inv_solve_abs_GN_logC

 INV_SOLVE_ABS_GNR absolute solver using Gauss Newton approximation
 img= gn_abs_solve( inv_model, data1, data2)
 img        => output image (or vector of images)
 inv_model  => inverse model struct
 data1      => EIT data

 Parameters:
   inv_model.parameters.max_iterations = N_max iter

0001 function img= inv_solve_abs_GN_logC( inv_model, data1);
0002 % INV_SOLVE_ABS_GNR absolute solver using Gauss Newton approximation
0003 % img= gn_abs_solve( inv_model, data1, data2)
0004 % img        => output image (or vector of images)
0005 % inv_model  => inverse model struct
0006 % data1      => EIT data
0007 %
0008 % Parameters:
0009 %   inv_model.parameters.max_iterations = N_max iter
0010 
0011 % (C) 2010 Andy Adler. License: GPL version 2 or version 3
0012 % $Id: inv_solve_abs_GN_logC.m 4039 2013-05-22 14:00:36Z aadler $
0013 
0014 if isstr(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0015 
0016 % Step 1: fit to background
0017 img = calc_jacobian_bkgnd( inv_model );
0018 % img = homogeneous_estimate( inv_model, data1 );
0019 if isfield(inv_model.fwd_model,'coarse2fine')
0020     nc = size(img.fwd_model.coarse2fine,2);
0021     img.elem_data = mean(img.elem_data)*ones(nc,1);
0022 end
0023 
0024 % Load the paramaters for the inversion
0025 if isfield(inv_model,'parameters')
0026     img.parameters= inv_model.parameters;
0027 else img.parameters.default= [];
0028 end
0029 
0030 if ~isfield(img.parameters,'normalisation')
0031     img.parameters.normalisation= 1;
0032 end
0033 
0034 if ~isfield(img.parameters,'perturb')
0035     img.parameters.perturb= [0 0.0001 0.001 0.01];
0036 end
0037 
0038 
0039 if isfield(img.parameters,'homogeneization')
0040     img = homogeneous_estimate( img, data1 );
0041 end
0042 
0043 
0044 hp  = calc_hyperparameter( inv_model );
0045 if isfield(img.parameters,'fixed_background') && img.parameters.fixed_background==1
0046     inv_model2= inv_model;
0047     inv_model2.fwd_model.coarse2fine= inv_model.fwd_model.coarse2fine(:,1:end-1);
0048     RtR = calc_RtR_prior( inv_model2 );
0049 else
0050     RtR = calc_RtR_prior( inv_model );
0051 end
0052 
0053 
0054 W   = calc_meas_icov( inv_model );
0055 hp2RtR= hp*RtR;
0056 
0057 iters = 1;
0058 if isfield(inv_model.parameters,'max_iterations')
0059    iters = inv_model.parameters.max_iterations; 
0060 end
0061 
0062 img0 = img;
0063 img0.logCond= log(img0.elem_data);
0064 
0065 residuals= zeros(size(data1,1),iters+1);
0066 
0067 for k = 1:iters  
0068     vsim=  fwd_solve(img);
0069     res = img.parameters.normalisation*(data1-vsim.meas);
0070     residuals(:,k)=res;
0071    
0072   % Calculate Jacobian
0073     disp(['Begin Jacobian computation - Iteration ' num2str(k)]);
0074     J = calc_jacobian( img ); 
0075 
0076     % Convert Jacobian as the adjusted parameters are the logarithm of the
0077     % conductivity
0078     img.logCond= log(img.elem_data);
0079     dCond_dlogCond= img.elem_data;
0080     J = J.*repmat((dCond_dlogCond),1,size(data1,1))';
0081     
0082     % Normalize the Jacobian
0083     J= img.parameters.normalisation*J;
0084     if isfield(img.parameters,'fixed_background') && img.parameters.fixed_background==1
0085         J= J(:,1:nc-1);
0086     end
0087 %     Jsum= sum(abs(J),1);
0088 %     Jmax= max(J); Jmax= (Jmax+abs(min(Jmax)));
0089 %     Jmax(Jmax== 0)= 1e-4;
0090 %     figure; hist(log10(Jmax))
0091 %     figure; hist(Jsum);
0092 %     img2= img; img2.fwd_model= inv_model.rec_model;
0093 %     img2.elem_data= Jsum;
0094 %     img2.calc_colours.clim= 25;
0095 %     img2.calc_colours.ref_level= 25;
0096 %     figure; show_fem(img2,3)
0097 %     img2.elem_data= log10(Jmax);
0098 %     img2.calc_colours.clim= 2;
0099 %     img2.calc_colours.ref_level= -1;
0100 %     figure; show_fem(img2,3)
0101 %
0102     if isfield(img.parameters,'fixed_background') && img.parameters.fixed_background==1
0103         RDx = hp2RtR*(img0.logCond(1:end-1) - img.logCond(1:end-1));
0104     else
0105         RDx = hp2RtR*(img0.logCond - img.logCond);
0106     end
0107 %     figure; plot(RDx);
0108     dx = (J'*W*J + hp2RtR)\(J'*res + RDx);
0109     if isfield(inv_model.parameters,'fixed_background') && inv_model.parameters.fixed_background==1
0110         dx(nc)= 0;
0111     end
0112     
0113     img= line_optimize(img, dx, data1);
0114 %     img0= img;
0115 end
0116 
0117 
0118 vsim=  fwd_solve(img);
0119 residuals(:,k+1) = img.parameters.normalisation*(vsim.meas-data1);
0120 img.residuals= residuals;
0121 img.estimation= vsim.meas;
0122 
0123 % Fit a parabola to the linefit and pick the best point
0124 % This is better than doing an exhaustive search
0125 % function  img = line_optimize(imgk, dx, data1);
0126 %   flist = [ 0.1,  0.5, 1.0];
0127 %   clim = mean(imgk.elem_data)/10; % prevent zero and negative conductivity
0128 %   img = imgk;
0129 %   for i = 1:length(flist);
0130 %      img.elem_data = imgk.elem_data + flist(i)*dx;
0131 %      img.elem_data(img.elem_data <= clim ) = clim;
0132 %      vsim = fwd_solve( img );
0133 %      dv = calc_difference_data( vsim , data1, img.fwd_model);
0134 %      mlist(i) = norm(dv);
0135 %   end
0136 %   pf = polyfit(flist, mlist, 2);
0137 %   fmin = -pf(2)/pf(1)/2; % poly minimum
0138 %   fmin(fmin>1) = 1; fmin(fmin<0) = 0;
0139 %
0140 %   img.elem_data = imgk.elem_data + flist(i)*dx;
0141 %   img.elem_data(img.elem_data <= clim ) = clim;
0142 function  img = line_optimize(imgk, dx, data1)
0143 img = imgk;
0144 perturb= img.parameters.perturb;
0145 % Compute the forward model for eah perturbation step
0146 mlist= perturb*0;
0147 for i = 1:length(perturb); 
0148     img.logCond= imgk.logCond + perturb(i)*dx;
0149     img.elem_data= exp(img.logCond);
0150     vsim = fwd_solve( img );
0151     dv = vsim.meas-data1;
0152     dv= img.parameters.normalisation*dv;
0153     mlist(i) = norm(dv);
0154 end
0155 % Select the best fitting step
0156 pf= polyfit(log10(perturb(2:end)), mlist(2:end), 2);
0157 fmin= 10.^(-pf(2)/pf(1)/2); % poly minimum
0158 p= linspace(log10(perturb(2)),log10(perturb(end)),50);
0159 cp= pf(1)*p.^2+pf(2)*p+pf(3);
0160 
0161 if pf(1)*(log10(fmin))^2+pf(2)*log10(fmin)+pf(3) > min(mlist(2:end)) || log10(fmin)>max(perturb) || log10(fmin)<min(perturb(2:end)) 
0162     [mk,ik]= min(mlist(2:end));
0163     fmin= perturb(ik+1);
0164 end
0165 
0166 img.logCond= imgk.logCond + fmin*dx;
0167 img.elem_data= exp(img.logCond);
0168 vsim = fwd_solve( img );
0169 dv = vsim.meas-data1;
0170 dv= img.parameters.normalisation*dv;
0171 
0172 if norm(dv) > mlist(1)
0173     img.parameters.perturb= [0 logspace(-4,-2,5)];
0174     img.logCond= imgk.logCond;
0175 else
0176       img.parameters.perturb= [0 fmin/2 fmin fmin*2];
0177       img.logCond= imgk.logCond + fmin*dx;
0178       if size(imgk.fwd_model.elems,1) <= 200000
0179           img.parameters.perturb= [0 fmin/4 fmin/2 fmin fmin*2 fmin*4];
0180       else
0181           img.parameters.perturb= [0 fmin/2 fmin fmin*2];
0182       end
0183 end
0184 
0185 figure; semilogx(perturb(2:end),mlist(2:end),'xk',fmin,pf(1)*log10(fmin)^2+pf(2)*log10(fmin)+pf(3),'or'); hold on
0186 semilogx(10.^p,pf(1)*(p).^2+pf(2)*p+pf(3),'k','linewidth',2); axis tight
0187 xlabel('alpha','fontsize',20,'fontname','Times')
0188 ylabel('Normalized residuals','fontsize',20,'fontname','Times')
0189 title({['Best alpha = ' num2str(fmin,'%1.2e')] ; ...
0190     ['norm no move = ' num2str(mlist(1),4)]},'fontsize',30,'fontname','Times')
0191 set(gca,'fontsize',20,'fontname','Times'); drawnow;
0192 
0193 % Record the corresponding parameters
0194 img.elem_data= exp(img.logCond);
0195 img.res_data= exp(-img.logCond);
0196 
0197 function img = homogeneous_estimate( img, data )
0198 %    img = calc_jacobian_bkgnd( imdl );
0199 %    vs = fwd_solve(img);
0200 %
0201 %    if isstruct(data)
0202 %       data = data.meas;
0203 %    else
0204 %      meas_select = [];
0205 %      try
0206 %         meas_select = imdl.fwd_model.meas_select;
0207 %      end
0208 %      if length(data) == length(meas_select)
0209 %         data = data(meas_select);
0210 %      end
0211 %    end
0212 %
0213 %    pf = polyfit(data,vs.meas,1);
0214 %
0215 %    if isfield(img.fwd_model,'coarse2fine');
0216 % % TODO: the whole coarse2fine needs work here.
0217 % %   what happens if c2f doesn't cover the whole region
0218 %       nc = size(img.fwd_model.coarse2fine,2);
0219 %       img.elem_data = mean(img.elem_data)*ones(nc,1)*pf(1);
0220 %    else
0221 %       img.elem_data = img.elem_data*pf(1);
0222 %    end
0223 
0224  conductivity_background= (img.parameters.normalisation*data)\(ones(size(data)));
0225  img.jacobian_bkgnd.value = conductivity_background;
0226  img.elem_data= ones(size(img.elem_data,1),1)*conductivity_background;
0227  disp(['Homogeneous resistivity = ' num2str(1/(conductivity_background),3) ' Ohm.m'])
0228  
0229     function do_unit_test
0230         unit_test_simdata
0231  
0232  
0233  function unit_test_simdata
0234 shape_str = ['solid top    = plane(0,0,0;0,1,0);\n' ...
0235              'solid mainobj= top and orthobrick(-100,-200,-100;410,10,100) -maxh=20.0;\n'];
0236 e0 = linspace(0,310,64)';
0237 elec_pos = [e0,0*e0,0*e0,1+0*e0,0*e0,0*e0];
0238 elec_shape= [0.1,0.1,1];
0239 elec_obj = 'top';
0240 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0241 %fmdl.nodes = fmdl.nodes(:,[1,3,2]);
0242 
0243 % 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];
0244 % 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];
0245 % fmdl.stimulation= stim_pattern_geophys( 64, 'Schlumberger', {'spacings', spacing,'multiples',multiples});
0246 
0247 fmdl.stimulation= stim_pattern_geophys( 64, 'Wenner', {'spacings', 1:32} );
0248 
0249 cmdl= mk_grid_model([], 2.5+[-30,5,20,30:10:290,300,315,340], ...
0250                             -[0:5:10 17 30 50 75 100]);
0251 % cmdl = mk_grid_model([], 2.5+[-50,-20,0:10:310,330,360], ...
0252 %                              -[0:2.5:10, 15:5:25,30:10:80,100,120]);
0253 c2f = mk_coarse_fine_mapping( fmdl, cmdl);
0254 S= sum(c2f,2); b= find(S<0.9999); a= find(S>=0.9999 & S<1);
0255 c2f(a,:)= c2f(a,:)./repmat(S(a),1,size(c2f,2));
0256 c2f(b,:)= 0; c2f(b,end+1)= 1;
0257 fmdl.coarse2fine= c2f;
0258 
0259 
0260 img = mk_image(fmdl,1);
0261 fm_pts = interp_mesh(fmdl);
0262 x_bary= fm_pts(:,1); z_bary= fm_pts(:,2);
0263 
0264 z_params= (min(fmdl.nodes(:,2)):max(fmdl.nodes(:,2)))';
0265 a = 0.36;
0266 b = 130;
0267 x_params= a*z_params+b;
0268 xlim=interp1(z_params,x_params,z_bary);
0269 img.elem_data(x_bary>xlim)= 0.01;
0270 
0271 % img2= mk_image(fmdl,img.elem_data);
0272 % figure; show_fem(img2);
0273 
0274 % img = mk_image(fmdl,0+ mk_c2f_circ_mapping(fmdl,[100;-30;0;50])*100);
0275 % img.elem_data(img.elem_data==0)= 0.1;
0276 dd  = fwd_solve(img);
0277 % show_fem(img);
0278 
0279 imdl= eidors_obj('inv_model','test');
0280 imdl.fwd_model= fmdl;
0281 imdl.rec_model= cmdl;
0282 imdl.fwd_model.normalize_measurements = 0;
0283 imdl.rec_model.normalize_measurements = 0;
0284 imdl.RtR_prior = @prior_laplace;
0285 imdl.solve = @inv_solve_abs_GN_logC;
0286 imdl.reconst_type = 'absolute';
0287 imdl.hyperparameter.value = 0.1;
0288 imdl.jacobian_bkgnd.value = 1;
0289 
0290 
0291 img1= mk_image(fmdl,1);
0292 vh1= fwd_solve(img1);
0293 normalisation= 1./vh1.meas;
0294 I= speye(length(normalisation));
0295 I(1:size(I,1)+1:size(I,1)*size(I,1))= normalisation;
0296 
0297 
0298 imdl.parameters.normalisation= I;
0299 imdl.parameters.homogeneization= 1;
0300 imdl.parameters.fixed_background= 1;
0301 imdl.parameters.perturb= [0 logspace(-5,-3,5)];
0302 imdl.parameters.max_iterations= 10;
0303 
0304 
0305 imgr= inv_solve(imdl, dd);
0306 
0307 imgGNd= imgr;
0308 imgGNd.fwd_model.coarse2fine= cmdl.coarse2fine;
0309 imgGNd.elem_data= log10(imgGNd.res_data(1:end-1));
0310 imgGNd.calc_colours.clim= 1.5;
0311 imgGNd.calc_colours.ref_level= 1.5;
0312 
0313 elec_posn= zeros(length(fmdl.electrode),3);
0314 for i=1:length(fmdl.electrode)
0315     elec_posn(i,:)= mean(fmdl.nodes(fmdl.electrode(1,i).nodes,:),1);
0316 end
0317    
0318 figure; show_fem(imgGNd,1);
0319 hold on; plot(elec_posn(:,1),elec_posn(:,3),'k*');
0320 axis tight; ylim([-100 0.5])
0321 xlabel('X (m)','fontsize',20,'fontname','Times')
0322 ylabel('Z (m)','fontsize',20,'fontname','Times')
0323 set(gca,'fontsize',20,'fontname','Times');
0324 
0325 img = mk_image( imdl );
0326 img.elem_data= imgr.elem_data;
0327 vCG= fwd_solve(img); vCG = vCG.meas;
0328 
0329 figure; plot(I*(dd.meas-vCG))
0330 figure; hist(I*(dd.meas-vCG),50)
0331 
0332 show_pseudosection( fmdl, I*dd.meas, '')
0333 show_pseudosection( fmdl, I*vCG, '')
0334 show_pseudosection( fmdl, (vCG-dd.meas)./dd.meas*100)
0335