inv_solve_abs_annealingMetropolis_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_annealingMetropolis_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 = Initial temperature
0014 %     inv_model.parameters.tempfinale = Final temperature
0015 %     inv_model.parameters.cooldown = Temperature decrease
0016 %     inv_model.parameters.normalisation
0017 
0018 % (C) 2012 Nolwenn Lesparre. License: GPL version 2 or version 3
0019 % $Id: inv_solve_abs_annealingMetropolis_params.m 5112 2015-06-14 13:00:41Z 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 % Estimate the number of parameters to adjust
0036 n_parameters= size(img.params_mapping.params,1);
0037 
0038 if isfield(inv_model.parameters,'temp')
0039     tempInit = inv_model.parameters.temp;
0040     tempfinale= inv_model.parameters.tempfinale;
0041     cooldown= inv_model.parameters.cooldown;
0042     nMetro= inv_model.parameters.nMetro;
0043 else
0044     tempInit = 1000;
0045     tempfinale= 0.01;
0046     cooldown= 0.95;
0047     nMetro= 300;
0048 end
0049 
0050 if ~isfield(inv_model.parameters,'normalisation')
0051     img.parameters.normalisation= 1;
0052 end
0053 
0054 if ~isfield(inv_model.params_mapping,'rangeParams')
0055     rangeParams= zeros(n_parameters,2);
0056     rangeParams(:,1)= -6; rangeParams(:,2)= 6;
0057 else
0058     rangeParams= inv_model.params_mapping.rangeParams;
0059 end
0060 
0061 if ~isfield(inv_model.params_mapping,'alpha')
0062     alpha= 1;
0063 else
0064     alpha= inv_model.params_mapping.alpha;
0065 end
0066 
0067 temp= tempInit; k= 1;
0068 while temp>=tempfinale
0069     temp= temp*cooldown;
0070     k=k+1;
0071 end
0072 
0073 disp(['Number of iteration ' num2str(k)])
0074 
0075 
0076 
0077 % Generates np+1 models with a Gaussian law with a mean equals to
0078 % inv_model.params_mapping.params and a standard deviation of
0079 % inv_model.params_mapping.perturb
0080 if isfield(inv_model.params_mapping,'inital_model')
0081     model= inv_model.params_mapping.inital_model;
0082 else
0083     np= numel(inv_model.params_mapping.params);
0084     model= 10.^((rangeParams(np,1).*(rangeParams(:,2)-rangeParams(:,1))+rangeParams(:,1)));    
0085 end
0086 
0087 % Estimate the cost of the initial model
0088 img.params_mapping.params= model;
0089 cost= objectiveFunction(img,data,alpha);
0090 costRef= cost;
0091 
0092 % Proceed to the Metropolis regression while reducing the temperature
0093 %% Run over the Simulated annealing loops
0094 
0095 temp= tempInit;
0096 n_temp= 1;
0097 while temp >= tempfinale % Loop over temperature
0098     for n= 1:nMetro % Metropolis loop
0099         model_try= model;
0100         idx_parametr= randi(n_parameters,1,1);
0101         model_try(idx_parametr)= (rand(1,1)*(rangeParams(idx_parametr,1))+rangeParams(idx_parametr,2));
0102         img.params_mapping.params= model_try;
0103         cost_try= objectiveFunction(img,data,alpha);
0104         if rand <= (exp(-cost_try/costRef)/exp(-cost/costRef))^(1/temp) % Metropolis test for acceptance of trial model
0105             likelihood= exp(-cost_try/costRef);
0106             cost= cost_try;
0107             model= model_try;
0108         end
0109     end % End of Metropolis loop
0110     resRec(n_temp)= cost;
0111     res_parRec(n_temp,:)= model';
0112     tempRec(n_temp)= temp;
0113     likelihoodRec(n_temp)= likelihood;
0114         
0115     n_temp= n_temp+1;
0116     temp = temp*cooldown;    % Decrease the temperature
0117 end % End of temperature loop
0118 res_parRec(:,3:4)= 10.^res_parRec(:,3:4);
0119 
0120 screenSize = get(0,'ScreenSize');
0121 h = figure;
0122 set(h,'Position',[0 0 screenSize(3)/2 screenSize(4)]);
0123 subplot(4,1,1)
0124 plot(tempRec,res_parRec(:,1),tempRec,res_parRec(:,2),'linewidth',2)
0125 set(gca,'fontsize',16,'fontname','Times','xDir','reverse','xScale','log')
0126 subplot(4,1,2)
0127 plot(tempRec,res_parRec(:,3),tempRec,res_parRec(:,4),'linewidth',2)
0128 set(gca,'fontsize',16,'fontname','Times','xDir','reverse','xScale','log','yScale','log')
0129 subplot(4,1,3)
0130 plot(tempRec,resRec,'linewidth',2)
0131 set(gca,'fontsize',16,'fontname','Times','xDir','reverse','xScale','log','yScale','log')
0132 subplot(4,1,4)
0133 plot(tempRec,likelihoodRec,'linewidth',2)
0134 set(gca,'fontsize',16,'fontname','Times','xDir','reverse','xScale','log')
0135 
0136 [res_parRec(end,1) res_parRec(end,2)]
0137 [res_parRec(end,3) res_parRec(end,4)]
0138 
0139 end
0140 
0141 function cost= objectiveFunction(img,data,alpha)
0142 vsim=  fwd_solve(img);
0143 residuals= img.parameters.normalisation*(vsim.meas-data);
0144 cost= sqrt(sum(residuals.^2))/alpha;
0145 end
0146 
0147 
0148 function do_unit_test
0149 shape_str = ['solid top    = plane(0,0,0;0,1,0);\n' ...
0150              'solid mainobj= top and orthobrick(-100,-200,-100;410,10,100) -maxh=20.0;\n'];
0151 e0 = linspace(0,310,64)';
0152 elec_pos = [e0,0*e0,0*e0,1+0*e0,0*e0,0*e0];
0153 elec_shape= [0.1,0.1,1];
0154 elec_obj = 'top';
0155 fmdl = ng_mk_gen_models(shape_str, elec_pos, elec_shape, elec_obj);
0156 fmdl.stimulation= stim_pattern_geophys( 64, 'Wenner', {'spacings', 1:32} );
0157 
0158 img= mk_image(fmdl,1/20);
0159 fm_pts= interp_mesh(fmdl);
0160 x_bary= fm_pts(:,1); z_bary= fm_pts(:,2);
0161 
0162 z_params= (min(fmdl.nodes(:,2)):max(fmdl.nodes(:,2)))';
0163 a = 0.36; b = 130;
0164 x_params= a*z_params+b;
0165 xlim=interp1(z_params,x_params,z_bary);
0166 img.elem_data(x_bary>xlim)= 1/120;
0167 
0168 dd  = fwd_solve(img);
0169 sig= sqrt(norm(dd.meas)); m= size(dd.meas,1);
0170 noise= .05;
0171 ddn= dd;
0172 ddn.meas = dd.meas + noise*sig*randn(m,1);
0173 
0174 img1= mk_image(fmdl,1);
0175 vh1= fwd_solve(img1);
0176 normalisation= 1./vh1.meas;
0177 I= speye(length(normalisation));
0178 I(1:size(I,1)+1:size(I,1)*size(I,1))= normalisation;
0179 
0180 a = 0.3; b = 150;
0181 res_params= log([10 100]');
0182 
0183 z_params= (min(fmdl.nodes(:,2)):max(fmdl.nodes(:,2)))';
0184 x_params= a*z_params+b;
0185 
0186 imdl = eidors_obj('inv_model','testNoisy');
0187 imdl.fwd_model= fmdl;
0188 imdl.params_mapping.params= [a; b ; res_params];
0189 imdl.params_mapping.perturb= [0.1; 50; 2 ; 2];
0190 
0191 imdl.params_mapping.function = @border_mapping;
0192 imdl.params_mapping.data.x_bary = x_bary;
0193 imdl.params_mapping.data.z_bary = z_bary;
0194 imdl.params_mapping.data.res_params = res_params;
0195 imdl.params_mapping.data.x_params = x_params;
0196 imdl.params_mapping.data.z_params = z_params;
0197 imdl.params_mapping.data.a = a;
0198 imdl.params_mapping.data.b = b;
0199 imdl.reconst_type= 'absolute';
0200 imdl.solve = @inv_solve_abs_annealingSimplex_params;
0201 imdl.normalize_measurements= 1;
0202 imdl.parameters.normalisation= I;
0203 
0204 imdl.parameters.temp= 1000;
0205 imdl.parameters.tempfinale= 1;
0206 imdl.parameters.cooldown= 0.97;% 0.95;
0207 imdl.parameters.nMetro= 1;
0208 imdl.jacobian_bkgnd.value = 1;
0209 % imgr= inv_solve(imdl, dd);
0210 imgr= inv_solve(imdl, ddn);
0211 img = mk_image( imdl );
0212 img.elem_data= imgr.elem_data;
0213 vAS= fwd_solve(img); vAS = vAS.meas;
0214 
0215 figure; hist(I*(dd.meas-vAS),50)
0216 show_pseudosection( fmdl, I*dd.meas, 'HorizontalDownward')
0217 show_pseudosection( fmdl, I*vAS, 'HorizontalDownward')
0218 show_pseudosection( fmdl, (vAS-dd.meas)./dd.meas*100,'HorizontalDownward')
0219 
0220 end
0221 
0222 
0223 function img = border_mapping(img)
0224 %% Function to be called to perform the mapping in the forward problem
0225 z= img.params_mapping.data.z_params;
0226 res= img.params_mapping.params(end-1:end);
0227 
0228 a= img.params_mapping.params(1);
0229 b= img.params_mapping.params(2);
0230 x= z*a+b; 
0231 
0232 xi= img.params_mapping.data.x_bary;
0233 zi= img.params_mapping.data.z_bary;
0234 xlim=interp1(z,x,zi);
0235 
0236 vi= zeros(size(img.fwd_model.elems,1),1) + res(1);
0237 vi(xi>xlim)= res(2);
0238 
0239 img.elem_data= exp(-vi);
0240 
0241 img.params_mapping.data.x_params = x;
0242 img.params_mapping.params= [a ; b ;  res];
0243 
0244 end
0245