inv_solve_abs_GN_constrain

Do Gauss Netwon Method with barrier to ensure positvity of the con
ductivity elements
INPUT v_h,v_i - simulated/measured voltages
      inv_model - inverse model structure
      tol - term tolerance (about 1e-5-1e-3 - default 1e-3)
      min_s - the minimal sigma per element (size(J,2))
      max_s - the maximal   ""        ""        ""   
      max_its - maximum iterations (default 1)

% PARAMETERS - MAKE AUTOMATIC

M Crabb and N Polydorides- 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_constrain(inv_model,meas_data)
0002 %Do Gauss Netwon Method with barrier to ensure positvity of the con
0003 %ductivity elements
0004 %INPUT v_h,v_i - simulated/measured voltages
0005 %      inv_model - inverse model structure
0006 %      tol - term tolerance (about 1e-5-1e-3 - default 1e-3)
0007 %      min_s - the minimal sigma per element (size(J,2))
0008 %      max_s - the maximal   ""        ""        ""
0009 %      max_its - maximum iterations (default 1)
0010 %
0011 %% PARAMETERS - MAKE AUTOMATIC
0012 %
0013 %M Crabb and N Polydorides- 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_constrain'));
0021 warning on backtrace
0022 
0023 %Convergence tolerance
0024 [maxiter, tol, min_s,max_s,rel_par,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 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 %Current image and the logistic equivalent
0045 n_cond=length(img_cur);
0046 log_img_cur=zeros(n_cond,1);
0047 for i=1:n_cond
0048     log_img_cur(i)=inv_logistic_f(img_cur(i),min_s,max_s,rel_par);
0049 end
0050                      
0051 %Calculate the Jacobian, prior matrix and hyperparameter
0052 RtR = calc_RtR_prior( inv_model );
0053 hp= calc_hyperparameter( inv_model );
0054 
0055 %Calculate the voltage difference data (meas-sim)
0056 volt_diff_meas_sim = calc_difference_data( sim_data, meas_data, inv_model.fwd_model);
0057 
0058 %Start the Gauss Newton iteration
0059 for i=1:maxiter
0060     %Print to screen if we want error
0061     if(show_iter==2);% && mod(i,ceil(maxiter/10))==0)
0062         fprintf(1,'Error at iteration %i is %f\n',i,norm(volt_diff_meas_sim));
0063     end
0064     
0065     %Calculate the Jacobian
0066     J = calc_jacobian( img_bkgnd);
0067                 
0068     %Compute the different vectors for method  (Polydorides 2012 Pg.10)
0069     d_s_d_m=zeros(n_cond,n_cond);
0070     for ii=1:length(log_img_cur)
0071         d_s_d_m(ii,ii)=d_logistic_f(log_img_cur(ii),min_s,max_s,rel_par);
0072     end
0073     
0074     %Multiply the partial derivative with Jacobian
0075     J=J*d_s_d_m;
0076     
0077     %Gradient of objective function (regularization term not needed)
0078     grad_obj = J'*(-volt_diff_meas_sim);
0079     
0080     %%TODO Implement with generic conductivity (see above)
0081     if(best_homog_ref==2)
0082         grad_obj=grad_obj - hp^2*RtR*(img_ref-img_cur);
0083     end
0084 
0085     %Hessian of objective function
0086     hess_obj = J'*J + hp^2*RtR;        
0087     
0088     %Compute search direction - negate gradient and do search
0089     grad_obj=-grad_obj; p_search = hess_obj \ grad_obj;
0090     
0091     %% Backtracking line search??
0092     if(bls==2) %No linesearch
0093         %Update the constrained conductivity
0094         log_img_cur = log_img_cur + p_search; 
0095         %Change variables to normal conductivity
0096         for iii=1:n_cond
0097             img_cur(iii)=logistic_f(log_img_cur(iii),min_s,max_s,rel_par);
0098         end        
0099         img_bkgnd.elem_data=img_cur;
0100     else
0101        %Line search parameters
0102        alpha=1.0; alpha_bls=0.1; beta_bls=0.5;
0103        
0104        %Create new candidate, forward solve and difference with measurements
0105        log_img_new = log_img_cur + alpha*p_search; 
0106        
0107        %Change variables to normal conductivity
0108        for iii=1:n_cond
0109            img_new(iii)=logistic_f(log_img_new(iii),min_s,max_s,rel_par);
0110        end               
0111        img_bkgnd.elem_data=img_new'; 
0112        sim_data_new=fwd_solve(img_bkgnd);
0113        volt_diff_meas_sim_new = calc_difference_data( sim_data_new, meas_data, inv_model.fwd_model);   
0114 
0115        %Calculate the functions for BLS
0116        beta_u_x_n_1 = beta_f(volt_diff_meas_sim_new);
0117        beta_u_x_n   = beta_f(volt_diff_meas_sim);
0118        grad_u_x_n   = p_search'*grad_obj;
0119        while(beta_u_x_n_1 > beta_u_x_n + alpha_bls*alpha*grad_u_x_n)
0120            %Shrink the linesearch parameter and create new candidate
0121            alpha=alpha*beta_bls;
0122            
0123            %Create new candidate, forward solve and difference with measurements
0124            log_img_new = log_img_cur + alpha*p_search;
0125        
0126            %Change variables to normal conductivity
0127            for iii=1:n_cond
0128                img_new(iii)=logistic_f(log_img_new(iii),min_s,max_s,rel_par);
0129            end               
0130            img_bkgnd.elem_data=img_new; 
0131            sim_data_new=fwd_solve(img_bkgnd);
0132            volt_diff_meas_sim_new = calc_difference_data( sim_data_new, meas_data, inv_model.fwd_model);   
0133            
0134            %Calculate the functions for BLS
0135            beta_u_x_n_1 = beta_f(volt_diff_meas_sim_new);               
0136            if(alpha < 1e-16)
0137                 error('Line Search failed');
0138            end
0139        end
0140        
0141        %Update the solution from the descent, decrease barrier and assign
0142        log_img_cur=log_img_new; img_cur = img_new'; img_bkgnd.elem_data=img_cur;
0143     end
0144        
0145     %Resolve model, find difference data and test convergence
0146     sim_data_new=fwd_solve(img_bkgnd);
0147     volt_diff_meas_sim = calc_difference_data( sim_data_new, meas_data, inv_model.fwd_model);   
0148     
0149     if norm(volt_diff_meas_sim)<tol; break; end  % test tolerance
0150 end
0151 
0152 %Create a data structure to return
0153 img.name= 'solved by mc_GN_box_solve';
0154 img.elem_data = img_cur;
0155 img.fwd_model= inv_model.fwd_model;
0156 img.type='image';
0157 
0158 
0159 
0160 %Objective function - voltage 2-norm (for linesearch)
0161 function [beta]=beta_f(diff_volt)
0162     %Objective function
0163     beta = 0.5*norm(diff_volt,2)^2;
0164 end
0165 
0166 %Logistic function, its inverse and partial derivatives
0167 function [logistic]=logistic_f(m_cur,min__s,max__s,relax_param)
0168     logistic = min__s + (max__s-min__s)/( 1 + exp(-m_cur/relax_param) );
0169 end
0170 
0171 function [d_logistic]=d_logistic_f(m_cur,min__s,max__s,relax_param)
0172     d_logistic = (max__s-min__s)/(   (1+exp(-m_cur/relax_param)) * ( (1+exp(m_cur/relax_param)) * relax_param));
0173 end
0174 
0175 function [inv_logistic]=inv_logistic_f(cond_cur,min__s,max__s,relax_param)
0176     inv_logistic = -relax_param*log( (cond_cur-max__s)/(min__s-cond_cur));
0177 end
0178 
0179 
0180 %Default parameters for the GN solver
0181 function [maxiter, tol,min_s,max_s,rel_par,show_iter,bls,best_homog,best_homog_ref] = get_parameters(inv_model)
0182    try
0183      maxiter= inv_model.parameters.max_iterations;
0184    catch
0185      maxiter= 1;
0186    end
0187 
0188    try
0189      tol = inv_model.parameters.term_tolerance;
0190    catch
0191      tol= 1e-3;
0192    end
0193    
0194    try
0195      min_s = inv_model.parameters.min_s;
0196    catch
0197      min_s= 1e-3;
0198    end
0199    
0200    try
0201      max_s = inv_model.parameters.max_s;
0202    catch
0203      tol= 1e3;
0204    end
0205    
0206    try
0207      rel_par = inv_model.parameters.rel_par;
0208    catch
0209      rel_par= 1.0;
0210    end
0211    
0212    try
0213      show_iter = inv_model.parameters.show_iterations;
0214    catch
0215      show_iter= 1;
0216    end
0217    
0218    try
0219       bls = inv_model.parameters.bls; 
0220    catch
0221       bls=1;
0222    end
0223    
0224    try
0225       best_homog = inv_model.parameters.best_homog; 
0226    catch
0227       best_homog=1;
0228    end
0229    
0230    try
0231       best_homog_ref = inv_model.parameters.best_homog_ref; 
0232    catch
0233       best_homog_ref=1;
0234    end
0235      
0236 end
0237 
0238 end