compute_delta_v2

 computes smallest positive delta such that
 delta1 = min_{Gamma^c} ( {(ak-pk)/(dk)}_+, {(-ak-pk)/(dk)}_+)
 delta2 = min_{Gamma} ( -delx/xk )
 delta = min(delta1,delta2)
 
 input: 
   required: pk, dk, ak, x, delx_vec
   optional: 
       shrinkage_flag (select stepsize by checking 
           0 -- only if an element in the active set shrinks to zero (outgoing)
           1 -- only if an inactive constraint becomes active (incoming)
           2 -- both shrinking elements in the active set
               and constraint violations of the inactive set
       nonneg (apply positivity constraint)
           0 -- no sign constraint on the solution
           1 -- consider only those constraints that become negative
           (i.e., causes signal amplitude to become positive)

0001 function out = compute_delta_v2(in)
0002 %
0003 % computes smallest positive delta such that
0004 % delta1 = min_{Gamma^c} ( {(ak-pk)/(dk)}_+, {(-ak-pk)/(dk)}_+)
0005 % delta2 = min_{Gamma} ( -delx/xk )
0006 % delta = min(delta1,delta2)
0007 %
0008 % input:
0009 %   required: pk, dk, ak, x, delx_vec
0010 %   optional:
0011 %       shrinkage_flag (select stepsize by checking
0012 %           0 -- only if an element in the active set shrinks to zero (outgoing)
0013 %           1 -- only if an inactive constraint becomes active (incoming)
0014 %           2 -- both shrinking elements in the active set
0015 %               and constraint violations of the inactive set
0016 %       nonneg (apply positivity constraint)
0017 %           0 -- no sign constraint on the solution
0018 %           1 -- consider only those constraints that become negative
0019 %           (i.e., causes signal amplitude to become positive)
0020 
0021 % output:
0022 %   delta   -- stepsize
0023 %   idelta  -- index of the element that enters or leave the support
0024 %   flag    -- 0    idelta leaves the support
0025 %              1    idelta enters the support
0026 %
0027 
0028 gamma_x = in.gamma; gamma_xc = in.gamma_c;
0029 pk = in.pk(gamma_xc); dk = in.dk(gamma_xc);
0030 if length(in.ak) == 1
0031     ak = in.ak*ones(length(gamma_xc),1);
0032 else
0033     ak = in.ak(gamma_xc);
0034 end
0035  
0036 delx_vec = in.delx_vec; x_k = in.x;
0037 
0038 % select delta flag... determine the type of shrinkage
0039 if isfield(in,'shrinkage_flag'); shrinkage_flag = in.shrinkage_flag; else shrinkage_flag = 2; end
0040 % non-negativity constraint?
0041 if isfield(in,'nonneg'); nonneg = in.nonneg; else nonneg = 0; end
0042 
0043 out = [];
0044        
0045 % Constraint violation (new element to add)
0046 if shrinkage_flag ~= 0 
0047     % For different values of regularization parameters
0048     
0049     % violating constraints are positive
0050     delta1_constr = (ak-pk)./(dk);
0051     % delta1_constr = delta1_constr(gamma_xc);
0052     delta1_pos_ind = find(delta1_constr>0);
0053     delta1_pos = delta1_constr(delta1_pos_ind);
0054     [delta1 i_delta1] = min(delta1_pos);
0055     if isempty(delta1)
0056         delta1 = inf;
0057     end
0058     
0059     % violating constraints are negative
0060     delta2_constr = (-ak-pk)./(dk);
0061     % delta2_constr = delta2_constr(gamma_xc);
0062     delta2_pos_ind = find(delta2_constr>0);
0063     delta2_pos = delta2_constr(delta2_pos_ind);
0064     [delta2 i_delta2] = min(delta2_pos);
0065     if isempty(delta2)
0066         delta2 = inf;
0067     end
0068     
0069     if delta1>delta2 || nonneg 
0070         % infact if we want non-negative solution, then there is no need to
0071         % compute delta1 at all.
0072         delta = delta2;
0073         idelta = gamma_xc(delta2_pos_ind(i_delta2));
0074     else
0075         delta = delta1;
0076         idelta = gamma_xc(delta1_pos_ind(i_delta1));
0077     end
0078     
0079     out.delta_in = delta; out.idelta_in = idelta;
0080     if shrinkage_flag == 1
0081         out.delta = delta; out.idelta = idelta;
0082         out.flag = 1;
0083         return;
0084     end    
0085 end
0086 
0087 if shrinkage_flag ~= 1
0088     delta3_constr = (-x_k(gamma_x)./delx_vec(gamma_x));
0089     delta3_pos_index = find(delta3_constr>0);
0090     [delta3 i_delta3] = min(delta3_constr(delta3_pos_index));
0091     out_x = gamma_x(delta3_pos_index(i_delta3));
0092     if isempty(delta3)
0093         delta3 = inf;
0094     end
0095     
0096     out.delta_out = delta3; out.idelta_out = out_x;
0097     
0098     if shrinkage_flag == 0
0099         out.delta = delta3; out.idelta = out_x;
0100         out.flag = 0;
0101         return;
0102     end
0103 end
0104 
0105 flag = 1;
0106 if delta3 > 0 && delta3 <= delta
0107     flag = 0;
0108     delta = delta3;
0109     idelta = out_x;
0110 end
0111 
0112 out.delta = delta; out.idelta = idelta;
0113 out.flag = flag;
0114