compute_delta

 computes smallest positive delta such that
 delta1 = min_{Gamma^c} ( {(ak-pk)/(dk-bk)}_+, {(-ak-pk)/(dk+bk)}_+)
 delta2 = min_{Gamma} ( -delx/xk )
 delta = min(delta1,delta2)
 
 input: 
   required: pk, dk, ak, bk, 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(in)
0002 %
0003 % computes smallest positive delta such that
0004 % delta1 = min_{Gamma^c} ( {(ak-pk)/(dk-bk)}_+, {(-ak-pk)/(dk+bk)}_+)
0005 % delta2 = min_{Gamma} ( -delx/xk )
0006 % delta = min(delta1,delta2)
0007 %
0008 % input:
0009 %   required: pk, dk, ak, bk, 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 if length(in.bk) == 1
0036     bk = in.bk*ones(length(gamma_xc),1);
0037 else
0038     bk = in.bk(gamma_xc);
0039 end
0040 delx_vec = in.delx_vec; x_k = in.x;
0041 
0042 % select delta flag... determine the type of shrinkage
0043 if isfield(in,'shrinkage_flag'); shrinkage_flag = in.shrinkage_flag; else shrinkage_flag = 2; end
0044 % non-negativity constraint?
0045 if isfield(in,'nonneg'); nonneg = in.nonneg; else nonneg = 0; end
0046 
0047 out = [];
0048        
0049 % Constraint violation (new element to add)
0050 if shrinkage_flag ~= 0 
0051     % For different values of regularization parameters
0052     
0053     % violating constraints are positive
0054     delta1_constr = (ak-pk)./(dk-bk);
0055     % delta1_constr = delta1_constr(gamma_xc);
0056     delta1_pos_ind = find(delta1_constr>0);
0057     delta1_pos = delta1_constr(delta1_pos_ind);
0058     [delta1 i_delta1] = min(delta1_pos);
0059     if isempty(delta1)
0060         delta1 = inf;
0061     end
0062     
0063     % violating constraints are negative
0064     delta2_constr = (-ak-pk)./(dk+bk);
0065     % delta2_constr = delta2_constr(gamma_xc);
0066     delta2_pos_ind = find(delta2_constr>0);
0067     delta2_pos = delta2_constr(delta2_pos_ind);
0068     [delta2 i_delta2] = min(delta2_pos);
0069     if isempty(delta2)
0070         delta2 = inf;
0071     end
0072     
0073     if delta1>delta2 || nonneg 
0074         % infact if we want non-negative solution, then there is no need to
0075         % compute delta1 at all.
0076         delta = delta2;
0077         idelta = gamma_xc(delta2_pos_ind(i_delta2));
0078     else
0079         delta = delta1;
0080         idelta = gamma_xc(delta1_pos_ind(i_delta1));
0081     end
0082     
0083     out.delta_in = delta; out.idelta_in = idelta;
0084     if shrinkage_flag == 1
0085         out.delta = delta; out.idelta = idelta;
0086         out.flag = 1;
0087         return;
0088     end    
0089 end
0090 
0091 if shrinkage_flag ~= 1
0092     delta3_constr = (-x_k(gamma_x)./delx_vec(gamma_x));
0093     delta3_pos_index = find(delta3_constr>0);
0094     [delta3 i_delta3] = min(delta3_constr(delta3_pos_index));
0095     out_x = gamma_x(delta3_pos_index(i_delta3));
0096     if isempty(delta3)
0097         delta3 = inf;
0098     end
0099     
0100     out.delta_out = delta3; out.idelta_out = out_x;
0101     
0102     if shrinkage_flag == 0
0103         out.delta = delta3; out.idelta = out_x;
0104         out.flag = 0;
0105         return;
0106     end
0107 end
0108 
0109 flag = 1;
0110 if delta3 > 0 && delta3 <= delta
0111     flag = 0;
0112     delta = delta3;
0113     idelta = out_x;
0114 end
0115 
0116 out.delta = delta; out.idelta = idelta;
0117 out.flag = flag;
0118