inv_solve_abs_annealingSimplex_params

 INV_SOLVE_ABS_ANNEALINGSIMPLEX_PARAMS absolute solver using the simplex annealing method. 
 This function operates with a mapping function linking the inverse
 parameters 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_annealingSimplex_params(inv_model, data)
0002 % INV_SOLVE_ABS_ANNEALINGSIMPLEX_PARAMS absolute solver using the simplex annealing method.
0003 % This function operates with a mapping function linking the inverse
0004 % parameters linking the inverse parameters to the forward elements.
0005 % The mapping function is defined in the inv_model.params_mapping structure.
0006 
0007 % img= annealingSimplex( inv_model, data)
0008 % img        => output image (or vector of images)
0009 % inv_model  => inverse model struct
0010 % data      => EIT data
0011 %
0012 % Parameters:
0013 %     inv_model.parameters.tempInit = N_max iter
0014 %     inv_model.parameters.tempfinale = vector with at least 3 variables
0015 %     inv_model.parameters.cooldown =
0016 %     inv_model.parameters.normalisation
0017 
0018 % (C) 2012 Nolwenn Lesparre. License: GPL version 2 or version 3
0019 % $Id: inv_solve_abs_annealingSimplex_params.m 6037 2019-12-30 22:26:51Z aadler $
0020 
0021 
0022 % Necessity to define a different parameterisation of the inverse problem
0023 % by respect to the forward problem one. In that case, a mapping function
0024 % rely the elements to reconstruct the medium conductivity
0025 
0026 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0027 
0028 if isfield(inv_model,'params_mapping') &&  isfield(inv_model.params_mapping,'function')
0029     mapping_function= inv_model.params_mapping.function;
0030     img= feval(mapping_function,inv_model);
0031 else
0032     error('The inverse model must contain a field "params_mapping" where a mapping function links the forward and inverse parameters');
0033 end
0034 
0035 if isfield(inv_model.parameters,'temp')
0036     tempInit = inv_model.parameters.temp;
0037     tempfinale= inv_model.parameters.tempfinale;
0038     cooldown= inv_model.parameters.cooldown;
0039     nMetro= inv_model.parameters.nMetro;
0040 else
0041     tempInit = 1000;
0042     tempfinale= 0.01;
0043     cooldown= 0.95;
0044     nMetro= 3;
0045 end
0046 
0047 if ~isfield(inv_model.parameters,'normalisation')
0048     img.parameters.normalisation= 1;
0049 end
0050 temp= tempInit; k= 1;
0051 while temp>=tempfinale
0052     temp= temp*cooldown;
0053     k=k+1;
0054 end
0055 D= zeros(k,1); temperature= zeros(k,1);
0056 niter= k;
0057 temp= tempInit; k= 1;
0058 
0059 % Estimate the number of parameters to adjust
0060 np= size(img.params_mapping.params,1);
0061 
0062 % Generates np+1 models with a Gaussian law with a mean equals to
0063 % inv_model.params_mapping.params and a standard deviation of
0064 % inv_model.params_mapping.perturb
0065 if isfield(inv_model.params_mapping,'inital_model')
0066     modeles= inv_model.params_mapping.inital_model;
0067 else
0068     modeles= randn(np,np+1).*repmat(inv_model.params_mapping.perturb,1,np+1) + ...
0069         repmat(inv_model.params_mapping.params,1,np+1);
0070     modelesr= modeles(end,:);
0071     keep= ismember(modeles(end,:)<modeles(end-1,:),1);
0072     modeles(end,keep)= modeles(end-1,keep);
0073     modeles(end-1,keep)= modelesr(keep);
0074 end
0075 modelesT= [];
0076 costhhi= [];
0077 modelesLo= [];
0078 
0079 % Estimate the cost of each model
0080 cost= zeros(1,np+1);
0081 for j= 1:np+1
0082     img.params_mapping.params= modeles(:,j);
0083     cost(j)= objectiveFunction(img,data);
0084 end
0085 dist= 1;
0086 % Proceed to the downhill simplex regression while reducing the temperature
0087 while temp>=tempfinale && dist >= 1e-5
0088     for i= 1:nMetro
0089         [costlo,ilo]= min(cost);
0090         [costhi,ihi]= max(cost);
0091         [costnhi]= max(setdiff(cost,costhi));
0092         bary= sum(modeles,2);
0093         % Reflexion
0094         [modtryReflection,costtryReflection]= deformation(modeles,bary,np,ihi,-1,img,data);
0095         boolean1= annealingProbability(costtryReflection,costlo,temp);
0096         boolean2= annealingProbability(costtryReflection,costnhi,temp);
0097         if  boolean1 %costtryReflection <= costlo
0098              % Dilatation
0099             [modtryExpansion,costtryExpansion]= deformation(modeles,bary,np,ihi,-2,img,data);
0100             boolean= annealingProbability(costtryExpansion,costtryReflection,temp);
0101             if boolean %costtryExpansion <= costtryReflection
0102                 modeles(:,ihi)= modtryExpansion;
0103                 cost(ihi)= costtryExpansion;
0104             else
0105                 modeles(:,ihi)= modtryReflection;
0106                 cost(ihi)= costtryReflection; 
0107             end
0108         elseif boolean2 % costtryReflection<=costnhi
0109              % Outward  contraction
0110             [modtryOutContraction,costtryOutContraction]= deformation(modeles,bary,np,ihi,-0.7,img,data);
0111             boolean= annealingProbability(costtryOutContraction,costtryReflection,temp);
0112             if boolean % costtryOutContraction<=costtryReflection
0113                 modeles(:,ihi)= modtryOutContraction;
0114                 cost(ihi)= costtryOutContraction;
0115             else
0116                 modeles(:,ihi)= modtryReflection;
0117                 cost(ihi)= costtryReflection;     
0118             end     
0119         else 
0120             % Inward contraction
0121             [modtryContraction,costtryContraction]= deformation(modeles,bary,np,ihi,0.7,img,data);
0122             boolean= annealingProbability(costtryContraction,costhi,temp);
0123             if boolean % costtryContraction <= costhi
0124                 modeles(:,ihi)= modtryContraction;
0125                 cost(ihi)= costtryContraction;
0126              else   % Nothing better -> shrinkage
0127                 modeles= (modeles+repmat(modeles(:,ilo),1,np+1))/2;
0128                 for j= 1:np+1
0129                     img.params_mapping.params= modeles(:,j);
0130                     cost(j)= objectiveFunction(img,data);
0131                 end
0132             end
0133         end
0134         costhi= max(cost);
0135         modelesT= [modelesT  modeles];
0136         costhhi= [costhhi ; costhi]; 
0137         [costlo,ilo]= min(cost);
0138         modelesLo= [modelesLo modeles(:,ilo)];
0139     end 
0140 %      bary1= sum(modelesT(:,end-np:end),1)/(np+1);
0141 %      bary2= sum(modelesT(:,end-2*(np+1)+1:end-(np+1)),1)/(np+1);
0142 %
0143 %      dist= sqrt(sum((bary1-bary2).^2));
0144 %      D(k)= dist;
0145      temperature(k)= temp;
0146 
0147     temp= temp*cooldown;
0148     niter= niter-1;
0149     k=k+1;
0150 %     disp(['Remaining iterations= ' num2str(niter) ' - Temperature= ' num2str(temp)])
0151 
0152     
0153 %     disp(['Remaining iterations= ' num2str(niter) ' - Temperature= ' num2str(temp) ' - Distance from previous iteration= ' num2str(dist)])
0154    
0155     if k>2 && D(k)>D(k-1); break; end
0156 end
0157 % modelesT1= reshape(modelesT(:,1),3,[]);
0158 % modelesT2= reshape(modelesT(:,2),3,[]);
0159 
0160 
0161 % figure; loglog(temperature,D); axis tight
0162 % xlabel('Temperature','fontsize',20,'fontname','Times')
0163 % ylabel('Models contraction','fontsize',20,'fontname','Times')
0164 % set(gca,'fontsize',15,'fontname','Times','xdir','reverse')
0165 % drawnow
0166 
0167 [costlo,ilo]= min(cost);
0168 img.params_mapping.params= modeles(:,ilo);
0169 img= feval(mapping_function,img);
0170 
0171 
0172 modeles(1:end-2,ilo);
0173 exp(modeles(end-1:end,ilo));
0174 
0175 v= 1:size(modelesLo,2);
0176 
0177 figure; plot(v,exp(modelesLo),'linewidth',2);
0178 xlabel('Iteration number','fontsize',20,'fontname','Times')
0179 ylabel('Inversion parameter','fontsize',20,'fontname','Times');
0180 set(gca,'fontsize',15,'fontname','Times');
0181 axis tight;
0182 
0183 
0184 
0185 figure; semilogy(reshape(costhhi,[],1),'k','linewidth',2)
0186 xlabel('Iteration number','fontsize',20,'fontname','Times')
0187 ylabel('Cost function','fontsize',20,'fontname','Times')
0188 set(gca,'fontsize',15,'fontname','Times'); drawnow
0189 
0190 
0191 end
0192 
0193 function [boolean,Pr]= annealingProbability(costtry,costhi,temp)
0194 E= (costtry-costhi)/temp;
0195 Ptry= exp(-E);
0196 Pr= min([1 Ptry]);
0197 boolean= 0;
0198 if Pr > rand ; boolean= 1; end
0199 end
0200 
0201 
0202 function cost= objectiveFunction(img,data)
0203 vsim=  fwd_solve(img);
0204 residuals= img.parameters.normalisation*(vsim.meas-data);
0205 cost= sqrt(sum(residuals.^2));
0206 end
0207 
0208 function [modtry,costtry]= deformation(modeles,bary,np,ihi,fac,img,data)
0209 
0210 fac1= (1-fac)/np;  % division by np since it is the weight of the barycenter
0211 fac2= fac1-fac;      % no division as weight of a single model
0212 % modification of the worst model
0213 modtry= bary*fac1-modeles(:,ihi)*fac2;
0214 
0215 imgtry= img;
0216 imgtry.params_mapping.params= modtry;
0217 costtry= objectiveFunction(imgtry,data);    
0218 end
0219 
0220 
0221 
0222 function do_unit_test
0223 shape_str = ['solid top    = plane(0,0,0;0,1,0);\n' ...
0224              'solid mainobj= top and orthobrick(-100,-200,-100;410,10,100) -maxh=20.0;\n'];
0225 e0 = linspace(0,310,64)';
0226 elec_pos = [e0,0*e0,0*e0,1+0*e0,0*e0,0*e0];
0227 elec_shape= [0.1,0.1,1];
0228 elec_obj = 'top';
0229 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0230 fmdl.stimulation= stim_pattern_geophys( 64, 'Wenner', {'spacings', 1:32} );
0231 
0232 img= mk_image(fmdl,1/20);
0233 fm_pts= interp_mesh(fmdl);
0234 x_bary= fm_pts(:,1); z_bary= fm_pts(:,2);
0235 
0236 z_params= (min(fmdl.nodes(:,2)):max(fmdl.nodes(:,2)))';
0237 a = 0.36; b = 130;
0238 x_params= a*z_params+b;
0239 xlim=interp1(z_params,x_params,z_bary);
0240 img.elem_data(x_bary>xlim)= 1/120;
0241 
0242 dd  = fwd_solve(img);
0243 sig= sqrt(norm(dd.meas)); m= size(dd.meas,1);
0244 noise= .05;
0245 ddn= dd;
0246 ddn.meas = dd.meas + noise*sig*randn(m,1);
0247 
0248 img1= mk_image(fmdl,1);
0249 vh1= fwd_solve(img1);
0250 normalisation= 1./vh1.meas;
0251 I= speye(length(normalisation));
0252 I(1:size(I,1)+1:size(I,1)*size(I,1))= normalisation;
0253 
0254 a = 0.3; b = 150;
0255 res_params= log([10 100]');
0256 
0257 z_params= (min(fmdl.nodes(:,2)):max(fmdl.nodes(:,2)))';
0258 x_params= a*z_params+b;
0259 
0260 imdl = eidors_obj('inv_model','testNoisy');
0261 imdl.fwd_model= fmdl;
0262 imdl.params_mapping.params= [a; b ; res_params];
0263 imdl.params_mapping.perturb= [0.1; 50; 2 ; 2];
0264 
0265 imdl.params_mapping.function = @border_mapping;
0266 imdl.params_mapping.data.x_bary = x_bary;
0267 imdl.params_mapping.data.z_bary = z_bary;
0268 imdl.params_mapping.data.res_params = res_params;
0269 imdl.params_mapping.data.x_params = x_params;
0270 imdl.params_mapping.data.z_params = z_params;
0271 imdl.params_mapping.data.a = a;
0272 imdl.params_mapping.data.b = b;
0273 imdl.reconst_type= 'absolute';
0274 imdl.solve = @inv_solve_abs_annealingSimplex_params;
0275 imdl.normalize_measurements= 1;
0276 imdl.parameters.normalisation= I;
0277 
0278 imdl.parameters.temp= 1000;
0279 imdl.parameters.tempfinale= 1;
0280 imdl.parameters.cooldown= 0.97;% 0.95;
0281 imdl.parameters.nMetro= 1;
0282 imdl.jacobian_bkgnd.value = 1;
0283 % imgr= inv_solve(imdl, dd);
0284 imgr= inv_solve(imdl, ddn);
0285 img = mk_image( imdl );
0286 img.elem_data= imgr.elem_data;
0287 vAS= fwd_solve(img); vAS = vAS.meas;
0288 
0289 figure; hist(I*(dd.meas-vAS),50)
0290 show_pseudosection( fmdl, I*dd.meas, 'HorizontalDownward')
0291 show_pseudosection( fmdl, I*vAS, 'HorizontalDownward')
0292 show_pseudosection( fmdl, (vAS-dd.meas)./dd.meas*100,'HorizontalDownward')
0293 
0294 end
0295 
0296 
0297 function img = border_mapping(img)
0298 %% Function to be called to perform the mapping in the forward problem
0299 z= img.params_mapping.data.z_params;
0300 res= img.params_mapping.params(end-1:end);
0301 
0302 a= img.params_mapping.params(1);
0303 b= img.params_mapping.params(2);
0304 x= z*a+b; 
0305 
0306 xi= img.params_mapping.data.x_bary;
0307 zi= img.params_mapping.data.z_bary;
0308 xlim=interp1(z,x,zi);
0309 
0310 vi= zeros(size(img.fwd_model.elems,1),1) + res(1);
0311 vi(xi>xlim)= res(2);
0312 
0313 img.elem_data= exp(-vi);
0314 
0315 img.params_mapping.data.x_params = x;
0316 img.params_mapping.params= [a ; b ;  res];
0317 
0318 end
0319