inv_solve_abs_CG_params

 INV_SOLVE_ABS_CG_PARAMS absolute solver using the conjugate gradient method with
 Least Square criterion. This function operates with a mapping function
 linking the inverse parameters to the forward elements. The mapping
 function is defined in the inv_model.params_mapping structure.

0001 function img= inv_solve_abs_CG_params( inv_model, data)
0002 % INV_SOLVE_ABS_CG_PARAMS absolute solver using the conjugate gradient method with
0003 % Least Square criterion. This function operates with a mapping function
0004 % linking the inverse parameters to the forward elements. The mapping
0005 % function is defined in the inv_model.params_mapping structure.
0006 
0007 % img= CG_params_solve( inv_model, data)
0008 % img        => output image (or vector of images)
0009 % inv_model  => inverse model struct
0010 % data      => EIT data
0011 %
0012 %
0013 % The conjugate gradient method is iterative and computes at each iteration
0014 % k the parameters:
0015 % p(k+1) = p(k) + alpha(k+1)* delta(k+1)
0016 % where alpha is the step length and delta the perturbation direction
0017 
0018 % Parameters:
0019 %     inv_model.parameters.max_iterations = N_max iter
0020 %     inv_model.parameters.perturb = vector with at least 3 variables to
0021 %     determine the step length
0022 %     inv_model.parameters.lambda = vector in logarithm space to estimate
0023 %     regularization defined by the L-curve criterion
0024 %     inv_model.parameters.normalisation = normalisation of the data sets,
0025 %     for instance the geometrical factor used in the apparent resistivity
0026 %     estimation
0027 
0028 % inv_model.params_mapping structure may at least contain:
0029 % params_mapping.function
0030 % params_mapping.params
0031 % params_mapping.perturb
0032 
0033 % (C) 2012 Nolwenn Lesparre. License: GPL version 2 or version 3
0034 % $Id: inv_solve_abs_CG_params.m 4039 2013-05-22 14:00:36Z aadler $
0035 
0036 
0037 % Necessity to define a different parameterisation of the inverse problem
0038 % by respect to the forward problem one. In that case, a mapping function
0039 % rely the elements to reconstruct the medium conductivity
0040 
0041 if isstr(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0042 
0043 if isfield(inv_model,'params_mapping') && isfield(inv_model.params_mapping,'function')
0044     mapping_function= inv_model.params_mapping.function;
0045     img= feval(mapping_function,inv_model);
0046 else
0047     error('The inverse model must contain a field "params_mapping" where a mapping function links the forward and inverse parameters');
0048 end
0049 
0050 if isfield(inv_model.parameters,'max_iterations')
0051     iters = inv_model.parameters.max_iterations;
0052 else
0053     iters = 1;
0054 end
0055 
0056 if ~isfield(inv_model.parameters,'perturb')
0057     img.parameters.perturb= [0.0001 0.001 0.01];
0058 end
0059 
0060 if ~isfield(inv_model.parameters,'lambda')
0061     img.parameters.lambda= logspace(-5,0,500);
0062 end
0063 
0064 if ~isfield(inv_model.parameters,'normalisation')
0065     img.parameters.normalisation= 1;
0066 end
0067 
0068 np= size(img.params_mapping.params,1);
0069 residuals= zeros(size(data,1),iters+1);
0070 for k= 1:iters
0071     % Estimate residuals between data and estimation
0072     img= feval(mapping_function,img);
0073     vsim=  fwd_solve(img);
0074     res = data-vsim.meas;
0075     residuals(:,k)= img.parameters.normalisation*res;
0076 %     norm(residuals(:,k))
0077     % Calculate jacobian
0078     disp(['Begin Jacobian computation - Iteration ' num2str(k)]);
0079 %     eidors_cache clear jacobian
0080     J = calc_jacobian( img );
0081     % Normalize the Jacobian
0082     J= img.parameters.normalisation*J;
0083 %     figure; plot(J);
0084     if k==1;
0085     if np<=7
0086         [U,S,V]= svd(J); svj= diag(S);
0087         svj = svj(svj>eps);%figure; semilogy(svj)
0088 %         tol= svj(np);
0089 %         tol= svj(round(np/2));
0090         tol= svj(end)
0091     elseif np<=100
0092         [U,S,V]= svd(J); svj= diag(S);
0093         svj = svj(svj>eps);
0094         tol= svj(round(np/2));
0095     else
0096         tol= svdAnalysisLcurvecrit(img.parameters.normalisation*data,img,J);
0097 %         tol= svdAnalysisLcurvecrit(residuals(:,k),img,J);
0098     end
0099     end
0100 %     delta_params= pinv(J,tol)*residuals(:,k);
0101     delta_params= pinv(J,tol)*(img.parameters.normalisation*res);
0102     if k== 1; figure; plot(delta_params,'x-'); drawnow; end
0103     % Compute the step length and adjust the parameters
0104     img = line_optimize(img, delta_params, data);
0105     [(img.params_mapping.params(1)) ; img.params_mapping.params(2) ; exp(img.params_mapping.params(end-1:end))]
0106 end
0107 %     [(img.params_mapping.params(1)) ; img.params_mapping.params(2) ; ...
0108 %         exp(img.params_mapping.params(end-1:end))]
0109 % end
0110 img.residuals= residuals;
0111 end
0112 
0113 
0114 function  img = line_optimize(imgk, dx, data1)
0115 img = imgk;
0116 perturb= img.parameters.perturb;
0117 mapping_function= img.params_mapping.function;
0118 % Compute the forward model for eah perturbation step
0119 mlist= perturb*0;
0120 for i = 1:length(perturb);
0121     img.params_mapping.params= imgk.params_mapping.params + perturb(i)*dx;
0122     img= feval(mapping_function,img);
0123     vsim = fwd_solve( img );
0124     dv = calc_difference_data( vsim , data1, img.fwd_model);% vsim.meas-data1;
0125     dv= img.parameters.normalisation*dv;
0126     mlist(i) = norm(dv);
0127 end
0128 % Select the best fitting step
0129 pf= polyfit(log10(perturb(2:end)), mlist(2:end), 2);
0130 fmin= 10.^(-pf(2)/pf(1)/2); % poly minimum
0131 p= linspace(log10(perturb(2)),log10(perturb(end)),50);
0132 cp= pf(1)*p.^2+pf(2)*p+pf(3);
0133 
0134 if pf(1)*(log10(fmin))^2+pf(2)*log10(fmin)+pf(3) > min(mlist(2:end)) || log10(fmin)>max(perturb)
0135     [mk,ik]= min(cp);
0136     fmin= 10.^(p(ik));
0137 end
0138 
0139 img.params_mapping.params= imgk.params_mapping.params + perturb(i)*dx;
0140 img= feval(mapping_function,img);
0141 vsim = fwd_solve( img );
0142 dv = vsim.meas-data1;
0143 dv= img.parameters.normalisation*dv;
0144 norm(dv)
0145 
0146 % if pf(1)*(log10(fmin))^2+pf(2)*log10(fmin)+pf(3) > mlist(1)
0147 if norm(dv) > mlist(1)
0148     img.parameters.perturb= [0 logspace(-4,-2,5)];
0149     img.params_mapping.params= imgk.params_mapping.params;
0150 else
0151 %       img.parameters.perturb= [0 fmin/2 fmin fmin*2];
0152     img.params_mapping.params= imgk.params_mapping.params + perturb(i)*dx;
0153     img.parameters.perturb= [0 fmin/4 fmin/2 fmin fmin*2 fmin*4];
0154 end
0155 
0156 img= feval(mapping_function,img);
0157 
0158 
0159 figure;
0160 semilogx(perturb,mlist,'xk',fmin,pf(1)*log10(fmin)^2+pf(2)*log10(fmin)+pf(3),'or'); hold on
0161 semilogx(10.^p,pf(1)*(p).^2+pf(2)*p+pf(3),'k','linewidth',2); axis tight
0162 % semilogx(perturb,mlist,'xk',fmin,pf(1)*fmin^2+pf(2)*fmin+pf(3),'or'); hold on
0163 % semilogx(p,pf(1)*p.^2+pf(2)*p+pf(3),'k','linewidth',2); axis tight
0164 xlabel('alpha','fontsize',30,'fontname','Times')
0165 ylabel('Normalized residuals','fontsize',30,'fontname','Times')
0166 title({['Best alpha = ' num2str(fmin,'%1.2e')] ; ...
0167     ['norm no move = ' num2str(round(mlist(1)))]},'fontsize',30,'fontname','Times')
0168 set(gca,'fontsize',30,'fontname','Times'); drawnow;
0169 
0170 end
0171 
0172 
0173 function tol= svdAnalysisLcurvecrit(data,img,J)
0174 % Compute the singular value decompisition
0175 [U,S,V]= svd(J); svj= diag(S);
0176 svj = svj(svj>eps);
0177 
0178 % Estimate the model parameters from a forward model linear approximation
0179 beta= U'*data;     %  adjust data
0180 beta= beta(1:length(svj));
0181 
0182 % Estimate the parameter of regularization
0183 lambda= img.parameters.lambda;
0184 
0185 % Determine the pinv tolerance following the L-curve criterion
0186 SVJ= repmat(svj,1,length(lambda));
0187 LAMBDA= repmat(lambda,length(svj),1);
0188 FI= SVJ.^2./(SVJ.^2+LAMBDA.^2);
0189 BETA= repmat(beta,1,length(lambda));
0190 XL2= sum((FI.*BETA./SVJ).^2,1);
0191 RES2= sum(((1-FI).*BETA).^2,1);
0192 n= XL2; p= RES2;
0193 nnp= (1-FI).*(FI.^2).*(BETA.^2)./SVJ.^2;
0194 np= -(4./lambda).*sum(nnp,1);
0195 kapa= (2*n.*p./np).*(lambda.^2.*np.*p+2*lambda.*n.*p+lambda.^4.*n.*np)./(lambda.^2.*n.^2+p.^2).^(3/2);
0196 
0197 resi= sqrt(RES2); xlambda = sqrt(XL2);
0198 [mk,ik]= min(kapa);
0199 [m,ist]= min(abs(svj-lambda(ik)));
0200 tol=svj(ist);
0201 
0202 figure;
0203 loglog(resi,xlambda,'k',resi(ik),xlambda(ik),'or','linewidth',2); axis tight; hold on
0204 yl= get(gca,'ylim');
0205 ylabel('Roughness','fontsize',30,'fontname','Times')
0206 xlabel('Residuals','fontsize',30,'fontname','Times')
0207 ylim(yl);
0208 title(['Best solution lambda=' num2str(lambda(ik),'%1.2e')],'fontsize',30,'fontname','Times')
0209 set(gca,'fontsize',30,'fontname','Times'); drawnow;
0210 end   
0211 
0212 
0213 function do_unit_test
0214 shape_str = ['solid top    = plane(0,0,0;0,1,0);\n' ...
0215              'solid mainobj= top and orthobrick(-100,-200,-100;410,10,100) -maxh=20.0;\n'];
0216 e0 = linspace(0,310,64)';
0217 elec_pos = [e0,0*e0,0*e0,1+0*e0,0*e0,0*e0];
0218 elec_shape= [0.1,0.1,1];
0219 elec_obj = 'top';
0220 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0221 fmdl.jacobian= @jacobian_diff_map_func;
0222 %fmdl.nodes = fmdl.nodes(:,[1,3,2]);
0223 
0224 % spacing= [1 1 1 2 3 3 4 4 5 6 6 7 8 8 9 10 10 11 12 12 13 14 14 15 16 17];
0225 % multiples= [1 2 3 2 1 5/3 1 2  1 1 7/6 1 1 10/8 1 1 12/10 1 1 13/12 1 1 15/14 1 1 1];
0226 % fmdl.stimulation= stim_pattern_geophys( 64, 'Schlumberger', {'spacings', spacing,'multiples',multiples});
0227 
0228 fmdl.stimulation= stim_pattern_geophys( 64, 'Wenner', {'spacings', 1:32} );
0229 img= mk_image(fmdl,1/20);
0230 fm_pts= interp_mesh(fmdl);
0231 x_bary= fm_pts(:,1); z_bary= fm_pts(:,2);
0232 
0233 z_params= (min(fmdl.nodes(:,2)):max(fmdl.nodes(:,2)))';
0234 a = 0.36; b = 130;
0235 x_params= a*z_params+b;
0236 xlim=interp1(z_params,x_params,z_bary);
0237 % xlim= 130;
0238 img.elem_data(x_bary>xlim)= 1/120;
0239 
0240 dd  = fwd_solve(img);
0241 % figure;  show_fem(img);
0242 
0243 img1= mk_image(fmdl,1);
0244 vh1= fwd_solve(img1);
0245 normalisation= 1./vh1.meas;
0246 I= speye(length(normalisation));
0247 I(1:size(I,1)+1:size(I,1)*size(I,1))= normalisation;
0248 
0249 a = 0.3; b = 150;
0250 res_params= log([10 100]');
0251 
0252 z_params= (min(fmdl.nodes(:,2)):max(fmdl.nodes(:,2)))';
0253 x_params= a*z_params+b;
0254 
0255 imdl = eidors_obj('inv_model','test');
0256 imdl.fwd_model= fmdl;
0257 imdl.params_mapping.params= [a; b ; res_params];
0258 imdl.params_mapping.perturb= [0.005; 2; 0.05 ; 0.5];
0259 % imdl.params_mapping.params= [b; res_params];
0260 % imdl.params_mapping.perturb= [2; 0.05; 0.5];
0261 
0262 imdl.params_mapping.function = @border_mapping;
0263 imdl.params_mapping.data.x_bary = x_bary;
0264 imdl.params_mapping.data.z_bary = z_bary;
0265 imdl.params_mapping.data.res_params = res_params;
0266 imdl.params_mapping.data.x_params = x_params;
0267 imdl.params_mapping.data.z_params = z_params;
0268 imdl.params_mapping.data.a = a;
0269 imdl.params_mapping.data.b = b;
0270 imdl.reconst_type= 'absolute';
0271 imdl.solve = @inv_solve_abs_CG_params;
0272 imdl.normalize_measurements= 1;
0273 imdl.parameters.normalisation= I;
0274 imdl.parameters.lambda= logspace(-5,5,1000);
0275 imdl.parameters.max_iterations = 10;
0276 imdl.parameters.perturb= [0 logspace(-2,0,5)];
0277 imdl.jacobian_bkgnd.value = 1;
0278 imgr= inv_solve(imdl, dd);
0279 
0280 img = mk_image( imdl );
0281 img.elem_data= imgr.elem_data;
0282 vCG= fwd_solve(img); vCG = vCG.meas;
0283 figure; hist(I*(dd.meas-vCG),50)
0284 show_pseudosection( fmdl, I*dd.meas, 'HorizontalDownward')
0285 show_pseudosection( fmdl, I*vCG, 'HorizontalDownward')
0286 show_pseudosection( fmdl, (vCG-dd.meas)./dd.meas*100','HorizontalDownward')
0287 
0288 end
0289 
0290 
0291 function img = border_mapping(img)
0292 %% Function to be called to perform the mapping in the forward problem
0293 z= img.params_mapping.data.z_params;
0294 res= img.params_mapping.params(end-1:end);
0295 
0296 a= img.params_mapping.params(1);
0297 b= img.params_mapping.params(2);
0298 % b= img.params_mapping.params(1);
0299 %
0300 % a= img.params_mapping.data.a;
0301 % b= img.params_mapping.params(1);
0302 x= z*a+b; 
0303 
0304 
0305 % x2= img.params_mapping.data.x_params2;
0306 % z2= img.params_mapping.data.z_params2;
0307 
0308 xi= img.params_mapping.data.x_bary;
0309 zi= img.params_mapping.data.z_bary;
0310 xlim=interp1(z,x,zi);
0311 % zlim=interp1(x2,z2,xi);
0312 
0313 % figure; plot(xi,zlim,'*')
0314 % hold on, plot(xlim,zi,'x')
0315 
0316 
0317 vi= zeros(size(img.fwd_model.elems,1),1) + res(1);
0318 % vi(xi>b)= res(2);
0319 
0320 vi(xi>xlim)= res(2);
0321 % vi(zi>zlim&xi>xlim+15)= res(1);
0322 
0323 img.elem_data= exp(-vi);
0324 
0325 % img.params_mapping.data.x_params = x;
0326 % img.params_mapping.params= [x(2:end) ;  res];
0327 
0328 img.params_mapping.data.x_params = x;
0329 % img.params_mapping.params= [b ;  res];
0330 img.params_mapping.params= [a ; b ;  res];
0331 % img.params_mapping.params= [b ;  res];
0332 
0333 end