primaldual_tvrecon_lsearch

PURPOSE ^

[rs,x]=primaldual_tvrecon_lsearch(inv_mdl, vmeas, ...

SYNOPSIS ^

function [rs,x]=primaldual_tvrecon_lsearch(inv_mdl, vmeas,maxiter,alpha1,alpha2, epsilon, beta, min_change)

DESCRIPTION ^

 [rs,x]=primaldual_tvrecon_lsearch(inv_mdl, vmeas, ...
             maxiter,alpha1,alpha2, epsilon, beta)
 11/02/01 By Andrea Borsic
 01/02/02 Modified: new steplength search on dual variable
 function rs=tvrecon(msh,c,vmeas,maxiter,alpha1,alpha2)

 ex. rs=pd_tvrecon(msh,c,vmeas,10,5e-3,1e-9); for simulated data
 ex. rs=pd_tvrecon(msh,c,vmeas,10,5e-3,5e-9); for tank data

 PARAMETERS
   alpha1 - hyperparameter for the initial Tikhonov step
   alpha2 - hyperparameter for the TV update
       A good alpha11 is 2e-5, a good alpha2 is 1e-9
   epsilon - termination tolerance (recommend 1e-6);
   beta    - initial value of approx: abs ~ sqrt(sqr+beta)
   min_change - minimum change in objective fcn
   rs - solution
   x  - dual solution

 PRIMAL DUAL IMPLEMENTATION

 Total variation reconstruction with norm-2 fitting of the measurements

 As reference see.. "A non linear primal-dual method for total variation-based image restoration" T.F. Chan, G.H. Golub, P. Mulet
 Note about notation: A is A' for us, and y is s, C is B, beta=mu^2

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [rs,x]=primaldual_tvrecon_lsearch(inv_mdl, vmeas, ...
0002               maxiter,alpha1,alpha2, epsilon, beta, min_change)
0003 
0004 % [rs,x]=primaldual_tvrecon_lsearch(inv_mdl, vmeas, ...
0005 %             maxiter,alpha1,alpha2, epsilon, beta)
0006 % 11/02/01 By Andrea Borsic
0007 % 01/02/02 Modified: new steplength search on dual variable
0008 % function rs=tvrecon(msh,c,vmeas,maxiter,alpha1,alpha2)
0009 %
0010 % ex. rs=pd_tvrecon(msh,c,vmeas,10,5e-3,1e-9); for simulated data
0011 % ex. rs=pd_tvrecon(msh,c,vmeas,10,5e-3,5e-9); for tank data
0012 %
0013 % PARAMETERS
0014 %   alpha1 - hyperparameter for the initial Tikhonov step
0015 %   alpha2 - hyperparameter for the TV update
0016 %       A good alpha11 is 2e-5, a good alpha2 is 1e-9
0017 %   epsilon - termination tolerance (recommend 1e-6);
0018 %   beta    - initial value of approx: abs ~ sqrt(sqr+beta)
0019 %   min_change - minimum change in objective fcn
0020 %   rs - solution
0021 %   x  - dual solution
0022 %
0023 % PRIMAL DUAL IMPLEMENTATION
0024 %
0025 % Total variation reconstruction with norm-2 fitting of the measurements
0026 %
0027 % As reference see.. "A non linear primal-dual method for total variation-based image restoration" T.F. Chan, G.H. Golub, P. Mulet
0028 % Note about notation: A is A' for us, and y is s, C is B, beta=mu^2
0029 %
0030 
0031 % (C) 2002-2006 Andrea Borsic. License: GPL version 2 or version 3
0032 % $Id: primaldual_tvrecon_lsearch.m 4838 2015-03-30 07:41:23Z aadler $
0033 
0034 % Initialisation
0035 fwd_model= inv_mdl.fwd_model;
0036 
0037 msh.TC = fwd_model.elems';
0038 msh.PC = fwd_model.nodes';
0039 
0040 decay_beta=0.7;
0041 
0042 % this is used for the line search procedure,
0043 % the last element, and the biggest must be one
0044 len=([0,1e-4,1e-3,1e-2,0.1,0.2,0.5,0.8,1]);
0045 
0046 % Inizialisation
0047 A=calc_R_prior( inv_mdl);
0048 
0049 n=size(A,1); % num_rows_L
0050 m=size(A,2); % num_elem
0051 
0052 % Create homogeneous model
0053 IM= eidors_obj('image','');
0054 IM.fwd_model= fwd_model;
0055 
0056 if 0 % static EIT - this code doesn't work yet
0057    scaling=vmeas\v_sim;
0058    s=s*scaling;
0059    %u=potentials(msh,s,c);
0060    %v_sim=measures(msh,u);
0061    v_sim= sim_measures( IM, s);
0062    de_v=vmeas-v_sim;
0063    %J=jacobian(msh,u);
0064    IM.s= s;
0065    J= calc_jacobian( IM );
0066    de_s=[J;alpha1*A]\[de_v;zeros(n,1)];
0067    s=s+de_s;
0068 else
0069    s= zeros(m,1); % solution
0070    x=zeros(n,1);  % dual variable
0071    J= calc_jacobian( calc_jacobian_bkgnd(inv_mdl) );
0072    IM.elem_data= s;
0073    IM.difference_rec = 1;
0074    IM.J = J;
0075    v_sim= sim_measures( IM, s);
0076    de_v=vmeas-v_sim;
0077    de_s=[J;alpha1*A]\[de_v;zeros(n,1)];
0078    s=s+de_s;
0079    scaling= 1;
0080 end
0081 
0082 %dispmsh(msh,s); colorbar; drawnow;
0083 
0084 % Iterative procedure
0085 
0086 terminate=0;
0087 iter=1;
0088 ind=1;
0089     rs(:,iter)=s;
0090 
0091 Obj_Fcn_old = inf;
0092 while (~terminate)&&(iter<maxiter)
0093     
0094 %   u=potentials(msh,s,c);
0095 %   v_sim=measures(msh,u);
0096 %   J=jacobian(msh,u);  - Jacobian is same in difference EIT
0097     v_sim= sim_measures( IM, s);
0098 %   J= calc_jacobian( IM );
0099 %   plot([v_sim, vmeas]);% pause
0100 
0101     z=A*s;    % This is an auxilliary variable
0102     
0103     eta=sqrt(z.^2+beta);
0104     
0105     grad=J'*(v_sim-vmeas);
0106     
0107     for i=1:n
0108         grad=grad+alpha2*A(i,:)'*A(i,:)*s/eta(i);
0109     end % for
0110     
0111     primal=sum(abs(z)); % we don't care here about 0.5*norm(de_v)
0112     dual=sum(x.*z);
0113 
0114     % TERMINATE ITERATIONS IF primal is no longer decreasing
0115     Obj_Fcn = norm(v_sim - vmeas)^2 + alpha2*primal;
0116     if abs(Obj_Fcn/Obj_Fcn_old - 1) < min_change
0117        eidors_msg('PDIPM: Breaking at iteration %d',iter,2);
0118        break
0119     else
0120        Obj_Fcn_old= Obj_Fcn;
0121     end
0122 
0123     
0124     eidors_msg('PDIPM: %2d & %1.3e & %1.3e & %1.3e & %1.3e & %1.3e & %1.3e & %1.3e & ',iter,primal,dual,primal-dual,norm(v_sim-vmeas),beta,len(ind),norm(grad),4);
0125         
0126     E=spdiags(eta,0,n,n);
0127     F=spdiags(ones(n,1)-(1./eta).*x.*z,0,n,n);
0128     
0129     B=(J'*J+alpha2*A'*inv(E)*F*A);
0130     
0131     %B=0.5*(B+B');
0132     
0133     de_s=-B\(alpha2*A'*inv(E)*z+J'*(v_sim-vmeas));
0134     
0135     ang=acos((dot(de_s,-grad)/(norm(de_s)*norm(-grad))))*(360/(2*pi));
0136     
0137     eidors_msg('PDIPM angle=%+3.1f deg',ang,4);
0138     
0139     % line search
0140     
0141     for k=1:length(len)
0142 %       meas_k = measures(msh,s+len(k)*de_s,c);
0143         meas_k= sim_measures( IM, s+len(k)*de_s);
0144         func_val(k)=0.5*norm( meas_k - vmeas )^2 + ...
0145                     alpha2*sum(abs(A*(s+len(k)*de_s)));
0146     end % for
0147     
0148     [temp,ind]=min(func_val);% disp(len(ind));
0149     
0150     % conductivity update
0151         
0152     s=s+len(ind)*de_s;
0153                
0154     de_x=-x+inv(E)*z+inv(E)*F*A*de_s;
0155     
0156     % dual step length rule
0157     
0158     lims=sign(de_x); % this are the limits (+1 or -1) torward x is pushed is de_x is added to it
0159     
0160     clearance=lims-x; % this is the signed distance to the limits
0161     
0162  % this protects against division, other values wil dominate, it doesn't affect the algorithm
0163     de_x(de_x==0)=1e-6;
0164         
0165  % stemps that will make on compunent of x exceed the limits of step*de_x is applied
0166     steps=clearance./de_x;
0167     
0168     idx= steps==0;
0169     steps(idx)=de_x(idx);
0170     
0171     % we need to pick up the smallest, and have some safety room
0172        
0173     x=x+min(1,0.99*min(steps))*de_x;
0174 
0175     if IM.difference_rec ==0 
0176     % Upper and lower limits enforcement
0177     s( s<0.01*scaling )=0.01*scaling;
0178     % and upper bounds, dynamic range=1e4
0179     s( s>100*scaling )=100*scaling;
0180     end
0181     
0182     beta=beta*decay_beta;       % beta is reduced
0183     decay_beta=decay_beta*0.8;  % the rate at wich beta is reduced is also adjusted
0184     if beta<1e-12
0185         beta=1e-12;
0186     end
0187     
0188 %    if norm(A*x)>gap(x,z) error('Rounding errors are spoiling the calculation, stopping.'); end % if
0189     
0190     % The primal-dual gap has been reduced and measures match
0191     if (sum(abs(z)-x.*z)<epsilon)&&(norm(v_sim-vmeas)<epsilon)
0192         terminate=1;
0193     end % if
0194     
0195     iter=iter+1;
0196 
0197     rs(:,iter)=s;
0198     
0199 end % while
0200 
0201 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0202 
0203 
0204 function v_sim= sim_measures( IM, s);
0205    if IM.difference_rec == 1
0206       v_sim = IM.J*s;
0207    else
0208       vh= fwd_solve( IM );
0209       IM.elem_data= s;
0210       vi= fwd_solve( IM );
0211 
0212       v_sim = calc_difference_data( vh, vi, IM.fwd_model);
0213    end

Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005