inv_solve_TV_pdipm

 INV_SOLVE_TV_PDIPM inverse solver for Andrea Borsic's
   Total Variation solver for use with difference EIT
 img= inv_solve_TV_pdipm( inv_model, data1, data2)
 img        => output image
 inv_model  => inverse model struct
 data1      => differential data at earlier time
 data2      => differential data at later time
 Parameters (see primaldual_tvrecon_lsearch for description)
   alpha1 (imdl.inv_solve_TV_pdipm.alpha1)
   alpha2 (set with imdl.hyperparameter.value)
   beta   (imdl.inv_solve_TV_pdipm.beta)
   want_dual_variable  (set to 1 if you want access to dual)
 Termination parameters
  max_iters =  inv_model.parameters.max_iteration (default 10)
      Max number of iterations before stopping
  min change = inv_model.parameters.min_change   (default 0)
      Min Change in objective fcn (norm(y-Jx)^2 + hp*TV(x)) before stopping

0001 function img= inv_solve_TV_pdipm( inv_model, data1, data2)
0002 % INV_SOLVE_TV_PDIPM inverse solver for Andrea Borsic's
0003 %   Total Variation solver for use with difference EIT
0004 % img= inv_solve_TV_pdipm( inv_model, data1, data2)
0005 % img        => output image
0006 % inv_model  => inverse model struct
0007 % data1      => differential data at earlier time
0008 % data2      => differential data at later time
0009 % Parameters (see primaldual_tvrecon_lsearch for description)
0010 %   alpha1 (imdl.inv_solve_TV_pdipm.alpha1)
0011 %   alpha2 (set with imdl.hyperparameter.value)
0012 %   beta   (imdl.inv_solve_TV_pdipm.beta)
0013 %   want_dual_variable  (set to 1 if you want access to dual)
0014 % Termination parameters
0015 %  max_iters =  inv_model.parameters.max_iteration (default 10)
0016 %      Max number of iterations before stopping
0017 %  min change = inv_model.parameters.min_change   (default 0)
0018 %      Min Change in objective fcn (norm(y-Jx)^2 + hp*TV(x)) before stopping
0019 
0020 % (C) 2002-2009 Andrea Borsic and Andy Adler. License: GPL version 2 or version 3
0021 % $Id: inv_solve_TV_pdipm.m 6247 2022-03-23 13:47:00Z aadler $
0022 
0023 
0024 p= get_params(inv_model);
0025 
0026 dva = calc_difference_data( data1, data2, inv_model.fwd_model);
0027 % TEST CODE -> Put elsewhere
0028 back_val = get_good_background(inv_model, data1);
0029 inv_model.jacobian_bkgnd.value= back_val;
0030 
0031 sol= [];
0032 for i=1:size(dva,2)
0033    [soln,dual_x]=primaldual_tvrecon_lsearch(inv_model, dva(:,i), ...
0034         p.maxiter,p.alpha1,p.alpha2, p.tol, p.beta, p.min_change);
0035 
0036    if ~p.keepiters
0037       soln=soln(:,end);
0038    end
0039 
0040    sol=[sol, soln];
0041 end
0042 
0043 img.name= 'solved by inv_solve_TV_pdipm';
0044 img.elem_data = sol;
0045 img.fwd_model= inv_model.fwd_model;
0046 try if inv_model.inv_solve_TV_pdipm.want_dual_variable
0047    img.dual_data = dual_x;
0048 end; end
0049 
0050 function p = get_params(inv_model);
0051    
0052    try   
0053        p.alpha1= inv_model.inv_solve_TV_pdipm.alpha1;
0054    catch
0055        p.alpha1= 1e-2;
0056    end
0057 
0058    try   
0059        p.beta= inv_model.inv_solve_TV_pdipm.beta;
0060    catch
0061        p.beta= 1e-4;
0062    end
0063 
0064    p.alpha2= calc_hyperparameter( inv_model);
0065 
0066    try   
0067        p.min_change = inv_model.parameters.min_change;
0068    catch
0069        p.min_change = 0;
0070    end
0071 
0072    try   
0073        p.maxiter = inv_model.parameters.max_iterations;
0074    catch
0075        p.maxiter= 10;
0076    end
0077 
0078    try   
0079        p.keepiters = inv_model.parameters.keep_iterations;
0080    catch
0081        p.keepiters= 0;
0082    end
0083 
0084    p.tol = 0; % TODO
0085 
0086 function back_val = get_good_background(inv_mdl, data1);
0087 
0088    % Create homogeneous model
0089    IM= eidors_obj('image','');
0090    IM.fwd_model= inv_mdl.fwd_model;
0091    s= ones(size(IM.fwd_model.elems,1),1);
0092    IM.elem_data= s;
0093 
0094    vsim= fwd_solve( IM);
0095    back_val=abs( data1\vsim.meas );
0096    back_val=1;
0097