%%%%%%%%%% adesc.m (Second Stage) %%%%%%%%%% 
% adesc.m  
% Ascent-descent algorithm for FIR Filter Design 
% Solve problem on sets of extremal points using a 
% descent method based on the work of Demjanov & Malozemov 
% and Wolfe. 
% The parameter acc_min can be set to the desired accuracy; the 
% default value is 0.012 
%
% @(#) adesc.m: version 5.0 created on 02/95                        
% (c) 1995 by Lina Karam and James McClellan
%
%main program
adesc_flg=crmz_input('Press <return> to continue and 0 to stop: ',1); 
if (adesc_flg)
  disp('           ***********************************************');
  disp('           *****    Invoking Second Ascent  Stage     ****');
  disp('           ***********************************************');

% Set the variable "bands".
 bands = zeros(size(vec_edges));
 for k = 1:length(vec_edges)
     if vec_edges(k)~=1
       bands(k) = find(vec_edges(k)==grd);
     else 
       bands(k) = length(grd);
     end
 end
 J = sqrt(-1);
 no_stp = 1;
 a = a(1:Lb);
 a = a(:);
 n2 = length(a);
 n = 2*length(a);
 mxl = 2*L+1;
 maj_it = 0;
 alpha = 1;
 nu = 0; 
 %%%%%%%%%%% Set initial accuracy parameters %%%%%%%%%%%%%
 acc_scale = max(abs(WT.*A)); 
 acc_min = 0.012;         %%% ACCURACY: USER SPECIFIED %%% 
 acc = acc_min;
 if (acc < acc_min) 
   acc = acc_min; 
 end;
 r_min = 0.5;
 r_o = acc_scale * r_min;
 if (r_o < r_min)
   r_o = r_min;
 end;
 r = r_o;
 epsi_min = 0.012;
 epsi_o = acc_scale * epsi_min;
 epsil = epsi_o; 
 epsilon = 1-0.005;
 %%%%%%%%%%% Construct initial subset %%%%%%%%%%%%%%%%%%%%%%
 HH_o = HH;
 e_max = max(abs(EE)); 
 [iext] = adesc_findset(EE,iext,bands,delta);
 sub_EE = EE(iext);
 sub_grd = grd(iext);
 sub_WT = WT(iext);
 sub_max = max(abs(sub_EE));
 fmax = iext;
 [jext] = adesc_findextr(EE,bands,epsilon); 
 %%%%%%%%%%% Check optimality of initial approximation %%%%%
 gext = adesc_grad(EE, jext, grd, WT, M_str, Lf, is_odd);
 [rext] = adesc_minpolytope(gext);
 no_stp = (norm(rext) > acc);
 %%%%%%%%%%% Begin major iteration %%%%%%%%%%%%%%%%%%%%%%%%%
 while (no_stp)
     jext = find((max(abs(sub_EE))-abs(sub_EE)) <= epsil);
     G = adesc_grad(sub_EE, jext, sub_grd, sub_WT, M_str, Lf, is_odd);
     [NrG] = adesc_minpolytope(G);
     norm_NrG = norm(NrG);
     if (norm_NrG <= r)  
      nu = nu + 1; 
      epsil = epsi_o./(2^nu);
      r = r_o./(2^nu);
      f_extr = find(abs(sub_EE)/sub_max >= epsilon);
      gext = adesc_grad(sub_EE,f_extr,sub_grd,sub_WT,M_str,Lf,is_odd);
      [rext] = adesc_minpolytope(gext);
      if (norm(rext) < acc) 
        maj_it = maj_it + 1;
        nu = 0;
        epsil = epsi_o; 
        r = r_o;  
        fprintf('iter: %3.0f   ',maj_it)
        fprintf('peak wgt error:  %5.2e  ',e_max)
        fprintf('subset peak wgt err: %5.2e\n',sub_max)
        [iext] = adesc_findset(EE,iext,bands,sub_max);
        sub_EE = EE(iext);
        sub_grd = grd(iext);
        sub_WT = WT(iext);
        sub_max = max(abs(sub_EE));
        if (length(iext) == length(fmax))
          no_stp = ~(all(iext==fmax));
        end;
        fmax = iext;
      end;
     else 
      d = - NrG./norm_NrG;
      d2 = d(1:n2)+J*d(n2+1:n);
      [HH_d] = adesc_reconst(d2,h_str,HH_str,M_str,tgrid,Ng,...
               L,Lf,Lc,Ls,Lb,is_odd,vec_edges,indx_edges);
      [HH_n,EE,e_max,alpha] = adesc_linsearch(alpha,iext,HH_o,HH_d,...
                              sub_max,A,WT,ifgrid);
      a = a + alpha*d2;
      HH_o = HH_n;
      sub_EE = EE(iext);
      sub_max = max(abs(sub_EE));
     end;             %%% if %%%            
 end;                 %%% while no_stp %%%    
 %%%%%%%%%%%   end major iteration %%%%%%%%%%%%%%%%%%%%%%%%%
 HH = HH_n;
 a = a(:); 
 h = eval(h_str);
 epsilon = 0.9*sub_max/e_max;
 iext = adesc_findextr(EE,bands,epsilon);
 jext = iext';
 delta = max(abs(EE));
 fprintf('Optimal Solution Obtained !\n')
end;       %%%%%%% end if %%%%%%% 
%%%%%%%%%%%   End Ascent-Descent Stage  %%%%%%%%%%%%%%%%%%%%%%%
