inv_solve_abs_CG_logc

 inv_solve_CG_logc absolute solver using the conjugate gradient method with
 L-Curve criterion. Convert the element conductivity into their
 logarithm

0001 function img= inv_solve_abs_CG_logc( inv_model, data)
0002 % inv_solve_CG_logc absolute solver using the conjugate gradient method with
0003 % L-Curve criterion. Convert the element conductivity into their
0004 % logarithm
0005 
0006 % img= CG_lr_solve( inv_model, data)
0007 % img        => output image (or vector of images)
0008 % inv_model  => inverse model struct
0009 % data      => EIT data
0010 %
0011 % The conjugate gradient method is iterative and computes at each iteration
0012 % k the parameters:
0013 % p(k+1) = p(k) + alpha(k+1)* delta(k+1)
0014 % where alpha is the step length and delta the perturbation direction
0015 
0016 % Parameters:
0017 %     inv_model.parameters.max_iterations = N_max iter
0018 %     inv_model.parameters.perturb = vector with at least 3 variables to
0019 %     determine the step length
0020 %     inv_model.parameters.lambda = vector in logarithm space to estimate
0021 %     regularization defined by the L-curve criterion
0022 %     inv_model.parameters.normalisation = normalisation of the data sets,
0023 %     for instance the geometrical factor used in the apparent resistivity
0024 %     estimation
0025 
0026 % (C) 2012 Nolwenn Lesparre. License: GPL version 2 or version 3
0027 % $Id: inv_solve_abs_CG_logc.m 4039 2013-05-22 14:00:36Z aadler $
0028 
0029 
0030 % Possibility to use a different parameterisation between the inverse and
0031 % the forward problems. In that case, a coarse2fine matrix rely the
0032 % elements to reconstruct the medium conductivity
0033 
0034 if isstr(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0035 
0036 % Construct the image
0037 img = calc_jacobian_bkgnd( inv_model );
0038 
0039 if isfield(inv_model.fwd_model,'coarse2fine')
0040     nc = size(img.fwd_model.coarse2fine,2);
0041     img.elem_data = mean(img.elem_data)*ones(nc,1);
0042 end
0043 
0044 % Load the paramaters for the inversion
0045 if isfield(inv_model,'parameters')
0046     img.parameters= inv_model.parameters;
0047 else img.parameters.default= [];
0048 end
0049     
0050 if isfield(img.parameters,'max_iterations')
0051     iters = inv_model.parameters.max_iterations;
0052 else
0053     iters = 1;
0054 end
0055 
0056 if ~isfield(img.parameters,'perturb')
0057     img.parameters.perturb= [0 0.0001 0.001 0.01];
0058 end
0059 
0060 if ~isfield(img.parameters,'lambda')
0061     img.parameters.lambda= logspace(-5,0,500);
0062 end
0063 
0064 if ~isfield(img.parameters,'normalisation')
0065     img.parameters.normalisation= 1;
0066 end
0067 
0068 if isfield(img.parameters,'homogeneization')
0069     img = homogeneous_estimate( img, data );
0070 end
0071 
0072 
0073 residuals= zeros(size(data,1),iters+1);
0074 for k= 1:iters
0075     % Calculate Jacobian
0076     disp(['Begin Jacobian computation - Iteration ' num2str(k)]);
0077 
0078     J = calc_jacobian( img ); 
0079 
0080     % Convert Jacobian as the adjusted parameters are the logarithm of the
0081     % conductivity
0082     img.logCond= log(img.elem_data);
0083     dCond_dlogCond= img.elem_data;
0084     J = J.*repmat((dCond_dlogCond),1,size(data,1))';
0085     
0086     % Normalize the Jacobian
0087     J= img.parameters.normalisation*J;
0088        
0089     % Some of the parameters may not be adjusted, as for instance the
0090     % background for huge forward models
0091     if isfield(img.parameters,'fixed_background') && img.parameters.fixed_background==1
0092         J= J(:,1:nc-1);
0093     end
0094    
0095     if k==1 && size(J,1)<size(J,2)
0096 %         if isfield(img.parameters,'tol')
0097 %             tol= img.parameters.tol;
0098 %         elseif ~isfield(img.parameters,'tol')
0099             tol= svdAnalysisLcurvecrit(img.parameters.normalisation*data,img,J);
0100 %         end
0101     end
0102     
0103     % Estimate residuals between data and estimation
0104     vsim=  fwd_solve(img);
0105     res = data-vsim.meas;
0106     residuals(:,k)= img.parameters.normalisation*res;
0107     
0108     % Compute the step length for the Conjugate Gradient increment
0109     if size(J,1)<size(J,2)
0110         delta_params= pinv(J,tol)*(img.parameters.normalisation*(res));
0111     else
0112         delta_params= pinv(J)*(img.parameters.normalisation*(res));
0113     end
0114 %     figure; plot(delta_params)
0115 
0116     if isfield(inv_model.parameters,'fixed_background') && inv_model.parameters.fixed_background==1
0117         delta_params(nc)= 0;
0118     end
0119     % Compute the step length and adjust the parameters
0120     img = line_optimize(img, delta_params, data);
0121 end
0122 
0123 vsim=  fwd_solve(img);
0124 residuals(:,k+1) = img.parameters.normalisation*(vsim.meas-data);
0125 img.residuals= residuals;
0126 img.estimation= vsim.meas;
0127 end
0128 
0129 
0130 function  img = line_optimize(imgk, dx, data1)
0131 img = imgk;
0132 perturb= img.parameters.perturb;
0133 % Compute the forward model for eah perturbation step
0134 mlist= perturb*0;
0135 for i = 1:length(perturb); 
0136     img.logCond= imgk.logCond + perturb(i)*dx;
0137     img.elem_data= exp(img.logCond);
0138     vsim = fwd_solve( img );
0139     dv = vsim.meas-data1;
0140     dv= img.parameters.normalisation*dv;
0141     mlist(i) = norm(dv);
0142 end
0143 % Select the best fitting step
0144 pf= polyfit(log10(perturb(2:end)), mlist(2:end), 2);
0145 fmin= 10.^(-pf(2)/pf(1)/2); % poly minimum
0146 p= linspace(log10(perturb(2)),log10(perturb(end)),50);
0147 cp= pf(1)*p.^2+pf(2)*p+pf(3);
0148 
0149 
0150 if pf(1)*(log10(fmin))^2+pf(2)*log10(fmin)+pf(3) > min(mlist(2:end)) || fmin>max(perturb) || fmin<min(perturb(2:end)) 
0151     [mk,ik]= min(mlist(2:end));
0152     fmin= perturb(ik+1);
0153 end
0154 
0155 img.logCond= imgk.logCond + fmin*dx;
0156 img.elem_data= exp(img.logCond);
0157 vsim = fwd_solve( img );
0158 dv = vsim.meas-data1;
0159 dv= img.parameters.normalisation*dv;
0160 
0161 if norm(dv) > mlist(1)
0162     img.parameters.perturb= [0 logspace(-4,-2,5)];
0163     img.logCond= imgk.logCond;
0164 else
0165       img.parameters.perturb= [0 fmin/2 fmin fmin*2];
0166       img.logCond= imgk.logCond + fmin*dx;
0167       if size(imgk.fwd_model.elems,1) <= 200000
0168           img.parameters.perturb= [0 fmin/4 fmin/2 fmin fmin*2 fmin*4];
0169       else
0170           img.parameters.perturb= [0 fmin/2 fmin fmin*2];
0171       end
0172 end
0173 
0174 % figure; semilogx(perturb(2:end),mlist(2:end),'xk',fmin,pf(1)*log10(fmin)^2+pf(2)*log10(fmin)+pf(3),'or'); hold on
0175 % semilogx(10.^p,pf(1)*(p).^2+pf(2)*p+pf(3),'k','linewidth',2); axis tight
0176 % semilogx(10.^p,p*0+mlist(1),':k','linewidth',2); axis tight
0177 % xlabel('alpha','fontsize',20,'fontname','Times')
0178 % ylabel('Normalized residuals','fontsize',20,'fontname','Times')
0179 % title({['Best alpha = ' num2str(fmin,'%1.2e')] ; ...
0180 %     ['norm no move = ' num2str(mlist(1),4)]},'fontsize',30,'fontname','Times')
0181 % set(gca,'fontsize',20,'fontname','Times'); drawnow;
0182 %
0183 % Record the corresponding parameters
0184 img.elem_data= exp(img.logCond);
0185 img.res_data= exp(-img.logCond);
0186 end
0187 
0188 
0189 function tol= svdAnalysisLcurvecrit(data,img,J)
0190 % Compute the singular value decompisition
0191 [U,S,V]= svd(J); svj= diag(S);
0192 svj = svj(svj>eps);
0193 
0194 % figure; plot(svj)
0195 
0196 % Estimate the model parameters from a forward model linear approximation
0197 beta= U'*data;     %  adjust data
0198 beta= beta(1:length(svj));
0199 
0200 % Estimate the parameter of regularization
0201 lambda= img.parameters.lambda;
0202 % lambda= logspace(-4,-1,500);
0203 % lambda= logspace(-6,-2,500);
0204 
0205 % Determine the pinv tolerance following the L-curve criterion
0206 SVJ= repmat(svj,1,length(lambda));
0207 LAMBDA= repmat(lambda,length(svj),1);
0208 FI= SVJ.^2./(SVJ.^2+LAMBDA.^2);
0209 BETA= repmat(beta,1,length(lambda));
0210 XL2= sum((FI.*BETA./SVJ).^2,1);
0211 RES2= sum(((1-FI).*BETA).^2,1);
0212 resi= sqrt(RES2); xlambda = sqrt(XL2);
0213 %
0214 n= XL2; p= RES2;
0215 nnp= (1-FI).*(FI.^2).*(BETA.^2)./SVJ.^2;
0216 np= -(4./lambda).*sum(nnp,1);
0217 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);
0218 
0219 % [mk,ik]= min(kapa);
0220 % [m,ist]= min(abs(svj-lambda(ik)));
0221 % tol=svj(ist);
0222 
0223 n= length(lambda);
0224 
0225 pf1 = polyfit(log10(resi(1:round(n*0.2))), log10(xlambda(1:round(n*0.2))), 1);
0226 cp1= pf1(1)*log10(resi)+pf1(2);
0227 
0228 pf2 = polyfit(log10(resi(round(n*0.8):end)), log10(xlambda(round(n*0.8):end)), 1);
0229 cp2= pf2(1)*log10(resi)+pf2(2);
0230 %
0231 [mk,ik]= min(abs(cp1-cp2));
0232 [m,ist]= min(abs(svj-lambda(ik)));
0233 tol= svj(ist);
0234 disp(['Tolerance value = ' num2str(tol,3)])
0235 % tol= tol*5;
0236 
0237 if isfield(img.parameters,'tol')
0238     [mu,iu]= min(abs(img.parameters.tol-lambda));
0239     [m,istu]= min(abs(svj-lambda(iu)));
0240 else
0241     [mu,iu]= min(abs(tol-lambda));
0242     [m,istu]= min(abs(svj-lambda(iu)));
0243 end
0244 % [ms,ist1]= min(abs(svj-1));
0245 % % [ms,ist1]= min(abs(svj-1.5))
0246 % [ms,ist2]= min(abs(svj-2));
0247 % % [ms,ist1]= min(abs(svj-5))
0248 % [ist1 ist2];
0249 
0250 figure;
0251 loglog(resi,xlambda,'k','linewidth',2); axis tight; hold on
0252 % [x,y]= ginput;
0253 % [mk,ik]= min(abs(resi-x));
0254 % [m,ist]= min(abs(svj-lambda(ik)));
0255 % tol= svj(ist)
0256 % tol= lambda(ik)
0257 for ii= 1:100:length(lambda)
0258     plot(resi(ii),xlambda(ii),'x','MarkerSize',10,'linewidth',2)
0259     text(resi(ii),xlambda(ii),num2str(lambda(ii),3),'fontsize',14,'fontname','Times')
0260 end
0261 loglog(resi(ik),xlambda(ik),'or','linewidth',2,'MarkerSize',10)
0262 % if isfield(img.parameters,'tol')
0263 loglog(resi(iu),xlambda(iu),'pr','linewidth',2,'MarkerSize',10)
0264 % end
0265 loglog(resi,10.^cp1,'k:',resi,10.^cp2,'k:','linewidth',2)
0266 yl= get(gca,'ylim');
0267 ylabel('Roughness','fontsize',20,'fontname','Times')
0268 xlabel('Residuals','fontsize',20,'fontname','Times')
0269 ylim(yl);
0270 title({['Best solution lambda=' num2str(lambda(ik),'%1.2e')]; ...
0271     ['Number of Singular vector involved = ' num2str(ist) ';' num2str(istu) ]},'fontsize',20,'fontname','Times')
0272 set(gca,'fontsize',20,'fontname','Times'); drawnow;
0273 
0274 % figure;
0275 % semilogx(resi,kapa,'k','linewidth',2); axis tight; hold on
0276 % yl= get(gca,'ylim');
0277 % ylabel('Kapa','fontsize',30,'fontname','Times')
0278 % xlabel('Residuals','fontsize',30,'fontname','Times')
0279 % ylim(yl);
0280 % title(['Best solution lambda=' num2str(lambda(ik),'%1.2e')],'fontsize',30,'fontname','Times')
0281 % set(gca,'fontsize',30,'fontname','Times'); drawnow;
0282 
0283 
0284 if isfield(img.parameters,'tol')
0285     tol= img.parameters.tol;
0286 end
0287 
0288 end  
0289 
0290 
0291 function img = homogeneous_estimate( img, data )
0292 %    img = calc_jacobian_bkgnd( imdl );
0293 %    vs = fwd_solve(img);
0294 %
0295 %    if isstruct(data)
0296 %       data = data.meas;
0297 %    else
0298 %      meas_select = [];
0299 %      try
0300 %         meas_select = img.fwd_model.meas_select;
0301 %      end
0302 %      if length(data) == length(meas_select)
0303 %         data = data(meas_select);
0304 %      end
0305 %    end
0306 %
0307 %    pf = polyfit(data,vs.meas,1);
0308 %
0309 %    if isfield(img.fwd_model,'coarse2fine');
0310 % % TODO: the whole coarse2fine needs work here.
0311 % %   what happens if c2f doesn't cover the whole region
0312 %       nc = size(img.fwd_model.coarse2fine,2);
0313 %       img.elem_data = mean(img.elem_data)*ones(nc,1)*pf(1);
0314 %    else
0315 %       img.elem_data = mean(img.elem_data)*pf(1);
0316 %    end
0317 %     img.jacobian_bkgnd.value = mean(img.elem_data)*pf(1);
0318 %     disp(['Homogeneous resistivity = ' num2str(1/(mean(img.elem_data)*pf(1)),3) ' Ohm.m'])
0319 
0320     conductivity_background= (img.parameters.normalisation*data)\(ones(size(data)));
0321     img.jacobian_bkgnd.value = conductivity_background;
0322     img.elem_data= ones(size(img.elem_data))*conductivity_background;
0323     disp(['Homogeneous resistivity = ' num2str(1/(conductivity_background),3) ' Ohm.m'])
0324 end
0325 
0326 function do_unit_test
0327    unit_test_simdata
0328 end
0329 
0330 function unit_test_simdata
0331 shape_str = ['solid top    = plane(0,0,0;0,1,0);\n' ...
0332              'solid mainobj= top and orthobrick(-100,-200,-100;410,10,100) -maxh=20.0;\n'];
0333 e0 = linspace(0,310,64)';
0334 elec_pos = [e0,0*e0,0*e0,1+0*e0,0*e0,0*e0];
0335 elec_shape= [0.1,0.1,1];
0336 elec_obj = 'top';
0337 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0338 %fmdl.nodes = fmdl.nodes(:,[1,3,2]);
0339 
0340 % 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];
0341 % 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];
0342 % fmdl.stimulation= stim_pattern_geophys( 64, 'Schlumberger', {'spacings', spacing,'multiples',multiples});
0343 
0344 fmdl.stimulation= stim_pattern_geophys( 64, 'Wenner', {'spacings', 1:32} );
0345 
0346 cmdl= mk_grid_model([], 2.5+[-30,5,20,30:10:290,300,315,340], ...
0347                             -[0:5:10 17 30 50 75 100]);
0348 % cmdl = mk_grid_model([], 2.5+[-50,-20,0:10:310,330,360], ...
0349 %                              -[0:2.5:10, 15:5:25,30:10:80,100,120]);
0350 c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0351 S= sum(c2f,2); b= find(S<0.9999); a= find(S>=0.9999 & S<1);
0352 c2f(a,:)= c2f(a,:)./repmat(S(a),1,size(c2f,2));
0353 c2f(b,:)= 0; c2f(b,end+1)= 1;
0354 fmdl.coarse2fine= c2f;
0355 
0356 
0357 img = mk_image(fmdl,1);
0358 fm_pts = interp_mesh(fmdl);
0359 x_bary= fm_pts(:,1); z_bary= fm_pts(:,2);
0360 
0361 z_params= (min(fmdl.nodes(:,2)):max(fmdl.nodes(:,2)))';
0362 a = 0.36;
0363 b = 130;
0364 x_params= a*z_params+b;
0365 xlim=interp1(z_params,x_params,z_bary);
0366 img.elem_data(x_bary>xlim)= 0.01;
0367 
0368 % img2= mk_image(fmdl,img.elem_data);
0369 % figure; show_fem(img2);
0370 
0371 % img = mk_image(fmdl,0+ mk_c2f_circ_mapping(fmdl,[100;-30;0;50])*100);
0372 % img.elem_data(img.elem_data==0)= 0.1;
0373 dd  = fwd_solve(img);
0374 % show_fem(img);
0375 
0376 imdl= eidors_obj('inv_model','test');
0377 imdl.fwd_model= fmdl;
0378 imdl.rec_model= cmdl;
0379 imdl.solve = @inv_solve_abs_CG_logc;
0380 imdl.reconst_type = 'absolute';
0381 imdl.jacobian_bkgnd.value = 1;
0382 
0383 img1= mk_image(fmdl,1);
0384 vh1= fwd_solve(img1);
0385 normalisation= 1./vh1.meas;
0386 I= speye(length(normalisation));
0387 I(1:size(I,1)+1:size(I,1)*size(I,1))= normalisation;
0388 
0389 % imdl.parameters.lambda= logspace(-5.5,-2,1000);
0390 imdl.parameters.lambda= logspace(-3,2,1000);
0391 imdl.parameters.perturb= [0 logspace(-1,0,5)];
0392 
0393 imdl.parameters.max_iterations= 2;
0394 imdl.parameters.normalisation= I;
0395 imdl.parameters.homogeneization=1;
0396 imdl.parameters.fixed_background= 1;
0397 
0398 imgr= inv_solve(imdl, dd);
0399 
0400 imgCGd= imgr;
0401 imgCGd.fwd_model.coarse2fine= cmdl.coarse2fine;
0402 imgCGd.elem_data= log10(imgCGd.res_data(1:end-1));
0403 imgCGd.calc_colours.clim= 1.5;
0404 imgCGd.calc_colours.ref_level= 1.5;
0405 
0406 elec_posn= zeros(length(fmdl.electrode),3);
0407 for i=1:length(fmdl.electrode)
0408     elec_posn(i,:)= mean(fmdl.nodes(fmdl.electrode(1,i).nodes,:),1);
0409 end
0410    
0411 figure; show_fem(imgCGd,1);
0412 hold on; plot(elec_posn(:,1),elec_posn(:,3),'k*');
0413 axis tight; ylim([-100 0.5])
0414 xlabel('X (m)','fontsize',20,'fontname','Times')
0415 ylabel('Z (m)','fontsize',20,'fontname','Times')
0416 set(gca,'fontsize',20,'fontname','Times');
0417 
0418 img = mk_image( imdl );
0419 img.elem_data= imgr.elem_data;
0420 vCG= fwd_solve(img); vCG = vCG.meas;
0421 
0422 figure; plot(I*(dd.meas-vCG))
0423 figure; hist(I*(dd.meas-vCG),50)
0424 
0425 show_pseudosection( fmdl, I*dd.meas, '')
0426 show_pseudosection( fmdl, I*vCG, '')
0427 show_pseudosection( fmdl, (vCG-dd.meas)./dd.meas*100)
0428 % show_pseudosection( fmdl, I*vCG, '')
0429 end