inv_solve_abs_GN_prior

 INV_SOLVE_ABS_GN_PRIOR inverse solver (WITH DIFFERENT PRIOR AT ITERATION!!!!!)
 img= mc_inv_solve_abs_GN( inv_model, data1, data2)
 img        => output image (or vector of images)
 inv_model  => inverse model struct
 data1      => simulated   data 
 data2      => measurement data

 both data1 and data2 may be matrices (MxT) each of
  M measurements at T times
 if either data1 or data2 is a vector, then it is expanded
  to be the same size matrix
M Crabb - 29.06.2012
TODO - Figure a nice interface to 
       (i) Allow a best fitting homogeneous (user may not want this)
       (ii) Allow a reference conductivity at each iteration
       (iii) Allow a global reference conductivity

0001 function [img,img_iteration]= inv_solve_abs_GN_prior( inv_model, meas_data)
0002 % INV_SOLVE_ABS_GN_PRIOR inverse solver (WITH DIFFERENT PRIOR AT ITERATION!!!!!)
0003 % img= mc_inv_solve_abs_GN( inv_model, data1, data2)
0004 % img        => output image (or vector of images)
0005 % inv_model  => inverse model struct
0006 % data1      => simulated   data
0007 % data2      => measurement data
0008 %
0009 % both data1 and data2 may be matrices (MxT) each of
0010 %  M measurements at T times
0011 % if either data1 or data2 is a vector, then it is expanded
0012 %  to be the same size matrix
0013 %M Crabb - 29.06.2012
0014 %TODO - Figure a nice interface to
0015 %       (i) Allow a best fitting homogeneous (user may not want this)
0016 %       (ii) Allow a reference conductivity at each iteration
0017 %       (iii) Allow a global reference conductivity
0018 warning off backtrace
0019 warning('EIDORS:experimental','%s is experimental, handle with care!',...
0020                 upper('inv_solve_abs_GN_prior'));
0021 warning on backtrace
0022 
0023 %Get parameters (default : 1 maxiter, 1e-3 tol,2 show_iter, 2 backtrack)
0024 [maxiter, tol, show_iter,bls,best_homog,best_homog_ref] = get_parameters(inv_model);
0025 
0026 %Do not show iterations
0027 if(show_iter==1); img_iteration=0; end
0028 
0029 %Background image, save as current image and simulate some data
0030 img_bkgnd= calc_jacobian_bkgnd( inv_model );
0031 img_cur=img_bkgnd.elem_data;
0032 sim_data=fwd_solve(img_bkgnd);
0033 
0034 %"Best" fitting homogeneous?
0035 if(best_homog==2) %BEST FIT
0036     sigma_opt=1/((sim_data.meas'*meas_data)/(sim_data.meas'*sim_data.meas));
0037     img_cur=img_cur*sigma_opt;
0038     img_bkgnd.elem_data=img_cur;
0039     sim_data=fwd_solve(img_bkgnd);
0040 end
0041 %TODO - Generalise this
0042 img_ref=img_cur;
0043                          
0044 %Calculate the Jacobian, prior matrix and hyperparameter
0045 RtR = calc_RtR_prior( inv_model );
0046 hp= calc_hyperparameter( inv_model );
0047 
0048 %Calculate the voltage difference data (meas-sim)
0049 volt_diff_meas_sim = calc_difference_data( sim_data, meas_data, inv_model.fwd_model);
0050 
0051 %Start the Gauss Newton iteration
0052 for i=1:maxiter
0053     %Print to screen if we want error
0054     if(show_iter==2);% && mod(i,ceil(maxiter/10)) == 0 )
0055         img_iteration{i}.error = norm(volt_diff_meas_sim);
0056         img_iteration{i}.name= sprintf('solved by mc_GN_solve iteration_%i',i);
0057         img_iteration{i}.elem_data = img_cur;
0058         img_iteration{i}.fwd_model= inv_model.fwd_model;
0059         img_iteration{i}.type='image';    
0060         
0061         fprintf(1,'Error at iteration %i is %f\n',i,norm(volt_diff_meas_sim));
0062 %        figure;show_fem(img_iteration{i});
0063     end
0064     
0065     %Calculate the Jacobian
0066     J = calc_jacobian( img_bkgnd);
0067 
0068     %Gradient of objective function (regularization term not needed)
0069     %grad_obj = J'*W*(-volt_diff_meas_sim);
0070     grad_obj = J'*(-volt_diff_meas_sim);
0071     
0072     %%TODO Implement with generic conductivity (see above)
0073     if(best_homog_ref==2)
0074         grad_obj=grad_obj - hp^2*RtR*(img_ref-img_cur);
0075     end
0076     
0077     %Hessian of objective function
0078     hess_obj = J'*J + hp^2*RtR;
0079     
0080     %Compute search direction - negate gradient and do search
0081     grad_obj=-grad_obj; p_search = hess_obj \ grad_obj;
0082     
0083     %% Backtracking line search??
0084     if(bls==2) %No linesearch
0085         img_cur = img_cur + p_search; img_bkgnd.elem_data=img_cur;
0086     else
0087        %Line search parameters
0088        alpha=1.0; alpha_bls=0.1; beta_bls=0.5;
0089        
0090        %Create new candidate, forward solve and difference with measurements
0091        img_new = img_cur + alpha*p_search; img_bkgnd.elem_data=img_new; 
0092        sim_data_new=fwd_solve(img_bkgnd);
0093        volt_diff_meas_sim_new = calc_difference_data( sim_data_new, meas_data, inv_model.fwd_model);   
0094 
0095        %Calculate the functions for BLS
0096        beta_u_x_n_1 = beta_f(volt_diff_meas_sim_new);
0097        beta_u_x_n   = beta_f(volt_diff_meas_sim);
0098        grad_u_x_n   = p_search'*grad_obj;
0099        while(beta_u_x_n_1 > beta_u_x_n + alpha_bls*alpha*grad_u_x_n)
0100            %Shrink the linesearch parameter and create new candidate
0101            alpha=alpha*beta_bls;
0102            
0103            %Create new candidate, forward solve and difference with measurements
0104            img_new = img_cur + alpha*p_search;
0105        
0106            %Forward solve on new data and calc difference with measure
0107            img_bkgnd.elem_data=img_new; 
0108            sim_data_new=fwd_solve(img_bkgnd);
0109            volt_diff_meas_sim_new = calc_difference_data( sim_data_new, meas_data, inv_model.fwd_model);   
0110            
0111            %Calculate the functions for BLS
0112            beta_u_x_n_1 = beta_f(volt_diff_meas_sim_new);               
0113            if(alpha < 1e-16)
0114                 error('Line Search failed');
0115            end
0116        end
0117        
0118        %Update the solution from the descent, decrease barrier and assign
0119        img_cur = img_new; img_bkgnd.elem_data=img_cur;
0120     end
0121        
0122     %Resolve model, find difference data and test convergence
0123     sim_data_new=fwd_solve(img_bkgnd);
0124     volt_diff_meas_sim = calc_difference_data( sim_data_new, meas_data, inv_model.fwd_model);   
0125     
0126     if norm(volt_diff_meas_sim)<tol; break; end  % test tolerance
0127 end
0128 
0129 %Create a data structure to return
0130 img.name= 'solved by mc_GN_solve';
0131 img.elem_data = img_cur;
0132 img.fwd_model= inv_model.fwd_model;
0133 img.type='image';
0134 
0135 
0136 
0137 %Objective function - voltage 2-norm (for linesearch)
0138 function [beta]=beta_f(diff_volt)
0139     %Objective function
0140     beta = 0.5*norm(diff_volt,2)^2;
0141 end
0142 
0143 
0144 %Default parameters for the GN solver
0145 function [maxiter, tol,show_iter,bls,best_homog,best_homog_ref] = get_parameters(inv_model)
0146    try
0147      maxiter= inv_model.parameters.max_iterations;
0148    catch
0149      maxiter= 1;
0150    end
0151 
0152    try
0153      tol = inv_model.parameters.term_tolerance;
0154    catch
0155      tol= 1e-3;
0156    end
0157    
0158    try
0159      show_iter = inv_model.parameters.show_iterations;
0160    catch
0161      show_iter= 1;
0162    end
0163    
0164    try
0165       bls = inv_model.parameters.bls; 
0166    catch
0167       bls=1;
0168    end
0169    
0170    try
0171       best_homog = inv_model.parameters.best_homog; 
0172    catch
0173       best_homog=1;
0174    end
0175    
0176    try
0177       best_homog_ref = inv_model.parameters.best_homog_ref; 
0178    catch
0179       best_homog_ref=1;
0180    end
0181        
0182 end
0183 
0184 end