/* PROJEC2.SRC - Projection Method (Minimum Weighted Residual Method)
** Version 1.1 May 1995
** (C) Copyright 1993-1995 by IHS Wien
** All Rights Reserved
**
**  Written by Michal Kejak
**
**=============================================================================
**
**  PROC PROJEC
**
**  FORMAT
**          { p,retcode } = PROJEC(&fres,mcons,xint,p0)
**
**  INPUT
**
**          &fres - pointer to a procedure that computes the residual function.
**                  The residual function embodies the Euler equations for the
**                  solved problem. This procedure must have two input
**                  arguments: a vector of values of variables, and parameters;
**                  and one output argument, the vector of values of the
**                  function evaluated at the input vectors of parameter
**                  and variables values.
**          mcons - row vector containig the basic constants of model:
**                  - number of state variables
**                  - vector of degrees of aproximation
**                  - number of control variables
**
**          xint  - vector of ranges of state variables for approximation
**
**          p0    - matrix of initial parameter values (see _prinit)
**
**  OUTPUT
**
**          p     - vector of parameters for approximation of optimal policy
**                  functions
**          retc  - return code:
**
**
**
**
**
**==============================================================================
**
**  GLOBAL VARIABLES
**
**  (default values in parentheses)
**  _prmeth - scalar, indicator of solution strategy:
**              = 1 NEWTON (Newton-Raphson)
**              = 2 start with steepest descent method and then continue
**                  with Newton method
**
**  _fres   - pointer to residual function
**
**  _prnst  - scalar, number of state variables (1)
**  _prncon - scalar, number of control variables (1)
**  _prdegs - vector (1,_prnst) of degrees of approximation (2)
**  _ch_coefs - matrix (n,n) of coefficients of Chebyshev polynomials,
**              n = max(_prdegs[i], i=1..._prnst) (0)
**  _prnpar - scalar, number of parameters (2)
**  _prvar  - scalar, variant of zeroes for 2nd degree approximation (2)
**  _prxint - vector (2,_prnst) of ranges of state variables (.5, 1.5)
**  _prtrd  - vector (1,_prnst), slope coefficient for linear transformation
**            of state variables
**  _prtre  - vector (1,_prnst), intercept coefficient for linear transformation
**            of state variables
**
**  _prinit - scalar, type of input of initial parameter values (0)
**              = 0     linear initial approx., matrix p0 (_prnst+1,_prcon)
**                      contains initial values for linear approximation,
**                      ie absolute term and linear coefficients for particular
**                      state variables
**
**              = 1     initial parameter values enter in the form of internal
**                      parameter representation, vector (_prnpar*_prncon,1).
**                      It is possible to read this information from file with
**                      name in the string variable _prldfn (extension .PAR).
**                      (It is useful to use solution results from
**                      the previous lower degree approximation which
**                      can be saved on disk (see _prsave and _prsvfn)).
**
**  _prldfn - string (max. 8 characters), the name of file where the initial
**                      values should be read from.
**
**  _prsave - scalar, indicator of saving results on disk (0)
**              = 0     computed parameters will not be saved on disk
**              = 1     computed parameters will be saved on disk (see _prsvfn)
**
**
**
**  _prsvfn - string, contains file name for saving computed parameters
**
**  _pritern- scalar, maximal number of iterations for Newton method (50)
**
**  _priters- scalar, maximal number of iterations for steepest descent
**            method (50)
**
**  _prtoln - scalar, the tolerance in convergence for Newton method (1e-10)
**
**  _prtols - scalar, the tolerance in convergence for steepest descent
**            method (.5) (after that switching to Newton method - _prmeth = 2)
**
**  _prretc - return code (0)
**                  = 0       normal convergence
**                  = 1       no convergence after _pritern (+ _priters)
**                            iterations
**
**  __output - scalar. If non-zero, output produced during iterations
**             is sent to the screen and/or output device.
**
**  __title  - string. A custom title to be printed at the top of
**             the iterations report. By default, only a generic
**             title will be printed.
**
**============================================================================*/

/* SOURCE CODE */

#include projec2.ext;

proc (2) = PROJEC(fres,mcons,xint,p0);
local a,p,retc,fres:proc,errmsg,prnst,prncon,prnpar,prdegs,save_vec,load_vec,
      t1,t2;
      
/* initialization of projection method */
t1      = time;
INIT_CON(mcons,xint);
_przero = ZERO_VEC;

p   = zeros(_prnpar*_prncon,1); a = p;
if (rows(p0)/=_prnst+1 OR cols(p0)/=_prncon) AND _prinit==0 OR
   (rows(p0)/=_prnpar*_prncon OR cols(p0)/=1) AND _prinit==1 AND _prldfn $=="";
    goto ERRORS("Wrong size of matrix for initial parameter values !!!");
endif;
if _prinit==0;
    p   = LINPAR(p0);
elseif _prinit==1;
    if _prldfn $/= "";
        load load_vec = ^_prldfn;
        prnst   = load_vec[1];
        prncon  = load_vec[2];
        prnpar  = load_vec[3];
        prdegs  = load_vec[4:3+prnst]';
        a       = load_vec[4+prnst:3+prnst+prnpar*prncon];
        if rows(p) /= _prnpar*_prncon;
            goto ERRORS("Wrong size of loaded matrix with initial parameters"\
                        " !!!");
        endif;
        p = PARAM(a,prnst,prncon,prnpar,prdegs);
/*        "a =" a;
        "p =" p;*/
    else;
        p = p0;
    endif;
else;
    goto ERRORS("Wrong value of variable _prinit !!!");
endif;
    if __output;
        call header("PROJEC","",_pr_ver);
/*        print chrs(45*ones(80,1));; */
        print "              Parameters";;
        print chrs(32*ones(32,1)) "      Residues";
        print chrs(45*ones(80,1));
        print;
    endif;

a       = COEF(p);
/* RESID(p); */
if _prmeth == 1;
    p = NEWTON(&RESID,p);
elseif _prmeth == 2;
    p = STEEP(&RESID,p);
    p = NEWTON(&RESID,p);
else;
endif;
a   = COEF(p);
OUTPUT_COEF(a);
/*----------------------------------
"residuum"; RESID(p);
----------------------------------*/
if _prsave;
    if _prsvfn $== "";
        goto ERRORS("Name of file for saving parameters is not given !!!");
    else;
        save_vec = _prnst|_prncon|_prnpar|_prdegs'|p;
        save ^_prsvfn = save_vec;
    endif;
endif;

t2    = time;
ELTIM(t1,t2);
retc  = _prretc;
retp(p,retc);

ERRORS:
    pop errmsg;
    errorlog errmsg;
    end;
endp;

proc (0) = ELTIM(t1,t2);
local t;
t=t2-t1;
if t[2]<0;
    t[1]=t[1]-1;
    t[2]=t[2]+60;
endif;
if t[3]<0;
    t[2]=t[2]-1;
    t[3]=t[3]+60;
endif;
format /rd 2,0;
if t[1]>0;
"Elapsed time: " t[1] " hours, " t[2] " minutes and " t[3] " seconds";
elseif t[2]>0;
"Elapsed time: " t[2] " minutes and " t[3] " seconds";
else;
"Elapsed time: " t[3] " seconds";
endif;
format /ro 16,8;
endp;

/* initialization  of model constants */
proc (0) = INIT_CON(mcons,xint);
local errmsg;
_prxint = xint;
_prnst  = mcons[1]; /* make check !!! */
_prncon = mcons[_prnst+2];
_prdegs = mcons[2:_prnst+1];
if not(cols(_prdegs)==_prnst);
    goto ERRORS("Wrong size of approximation degrees vector !!!");
endif;
_prnpar      = MULT(1,_prdegs);
_prtrd       = zeros(1,_prnst); _prtre = _prtrd;
_prtrd[1,.]  = 2/(_prxint[2,.]-_prxint[1,.]);
_prtre[1,.]  = -(_prxint[1,.]+_prxint[2,.])./(_prxint[2,.]-_prxint[1,.]);

/* --- computation of matrix of Chebyshev polynomial 
   --- coefficients later used throughout the program */

_ch_coefs = CH_COEF_INIT(maxc(_prdegs'));

retp;

ERRORS:
    pop errmsg;
    errorlog errmsg;
    end;
endp;

proc (0) = OUTPUT_COEF(coeff);
    local i,fmtstr,values,strings,str_num,all_degs,x_degs,y_degs,z_degs,order;
    print chrs(42*ones(80,1));
    "   Ordinary (external) representation of policy functions' coefficients:";
    print chrs(42*ones(80,1));

    fmtstr = {"-*.*s" 20 15,
              "-*.*s" 15 10};
    str_num = 0~0;
    strings = zeros(_prncon+1,1) $+ ("Element"$|"Contr #" $+ ftocv(seqa(1,1,_prncon),1,0));
    strings = strings';
    if _prncon>1;          @ expansion @
        fmtstr = fmtstr|(ones(_prncon-1,1)*fmtstr[2,.]);    
        str_num = str_num~(ones(_prncon-1,1)*str_num[2]); 
    endif;
    call printfm(strings,str_num,fmtstr);
    print;
    print chrs(42*ones(80,1));

    fmtstr = {"-*.*s" 20 15,
              "*.*lf" 15 4};
    str_num = 0~1;          
    if _prncon>1;          @ expansion @
        fmtstr = fmtstr|(ones(_prncon-1,1)*fmtstr[2,.]);    
        str_num = str_num~(ones(_prncon-1,1)*str_num[2]); 
    endif;

    @ output is not restricted in degrees of approximation @
    
    order = seqa(0,1,_prdegs[1]);
    x_degs = 0 $+ "x" $+ "^" $+ ftocv(order,1,0)$+" ";
    if _prnst>1;
        order = seqa(0,1,_prdegs[2]);
        y_degs = 0 $+ "y" $+ "^" $+ ftocv(order,1,0)$+" ";
    elseif _prnst>2;    
        order = seqa(0,1,_prdegs[3]);
        z_degs = 0 $+ "z" $+ "^" $+ ftocv(order,1,0)$+" ";
    endif;
   
    if _prnst==1;
      values = coeff';
      strings = x_degs;
    else;  
      strings = x_degs;
      {values,strings} = COEF_TRANSFORM(coeff',x_degs,y_degs,z_degs);
    endif;
    
    call printfm(strings~values,str_num,fmtstr);
    print chrs(42*ones(80,1));

endp;

/* rearrangement of the coefficients into the order for output */
/* only for _prnst = 2 or 3 */
proc (2) = COEF_TRANSFORM(a,xs,ys,zs);
@ n = _prdegs --> a = (_prnpar x _prncon) = (n[1]*n[2]*n[3] x _prncon) @
    local n,ordx,ordy,ordz,degx,degy,degz,degin,degord,i,last,lastdeg,counter,res,strin,resstr;
    n = ones(1,3);
    degin = zeros(_prnpar,1);
    degord = zeros(_prnpar,1);
    n[1] = _prdegs[1];
    n[2] = _prdegs[2];
    if _prnst==2;
        n[3] = 1;
    else;                       @ _prnst == 3 @
        n[3] = _prdegs[3];
    endif;
    ordx = seqa(0,1,n[1]);
    ordy = seqa(0,1,n[2]);
    ordz = seqa(0,1,n[3]);
    degx = ordx.*.ones(n[2]*n[3],1);
    xs   = xs.*.ones(n[2]*n[3],1);
    degy = ones(n[1],1).*.(ordy.*.ones(n[3],1));
    ys   = ones(n[1],1).*.(ys.*.ones(n[3],1));
    degz = ones(n[1]*n[2],1).*.ordz;
    if _prnst==3;
        zs   = ones(n[1]*n[2],1).*.zs;
    endif;
    degin = degx+degy+degz;  @ = (_prnpar x 1); total degree for all state vars @
    
    if _prnst==2;
        strin = xs $+ ys;
    else;                       @ _prnst == 3 @
        strin = xs $+ ys $+ zs;
    endif;
    
    last = rows(degin);     @ =n[1]n[2]n[3] @
    lastdeg = degin[last];  
    counter = 1;
    for i(0,lastdeg,1);     @ no. of degree @
     for j(last,1,-1);      @ no. of element @
        if degin[j]==i;
          degord[j] = counter;
          counter = counter + 1;
        endif;
     endfor; 
    endfor;
    
    for i(1,last,1);
      if cols(a)>1; 
       res[degord[i],.]=a[i,.]; @ rows @
      else;
       res[degord[i]]=a[i];     @ scalars @
      endif; 
      resstr[degord[i]]=strin[i];   @ strings @
    endfor;

    retp(res,resstr);
endp;

proc PARAM(a,prnst,prncon,prnpar,prdegs);
local npar,i,ic,ip,ip1,is,m,ip2,b;
npar = prnpar;
b    = zeros(_prnpar*_prncon,1);
ip   = 1;
do while ip<=prnpar*prncon;
    ip1 = ip-1;
    ic  = trunc(ip1/npar);
    ip1 = ip1%npar;
    ip2 = ic*MULT(1,_prdegs);
    i   = 1;
    do while i<= prnst-1;
        m   = MULT(i+1,prdegs);
        is  = trunc(ip1/m);
        ip1 = ip1%m;
        ip2 = ip2 + is*MULT(i+1,_prdegs);
        i   = i + 1;
    endo;
    ip2    = ip2 + ip1 + 1;
    b[ip2] = a[ip];
    ip     = ip + 1;
endo;
retp(b);
endp;

/* Parameter transformation */
proc LINPAR(k);             
@ k=(_prnst+1,_prcon) - matrix of init. values for lin. approx. @
@ its first row contains abs.term for a correspondent control var @
    local a,a_matrix,k_coef,aa,n,s,i_matrix,i_col,i;
    
    a    = zeros(_prnpar*_prncon,1);
    a_matrix = zeros(_prnpar,_prncon);
    k_coef = k[2:_prnst+1,.];               @ (_prnst,_prncon) @
    aa = k_coef./(_prtrd'*ones(1,_prncon)); @ (_prnst,_prncon) @
    a_matrix[1,.] = k[1,.]-_prtre*aa;       @ (1,_prncon); contains "residual" @
    n = _prnst;                             @ # of state vars @
    s = -1*seqa(n,-1,n);                    @ shift @
    i_matrix = shiftr(_prdegs'*ones(1,n),s,ones(n,1));
    i_col = prodc(i_matrix)+ones(n,1);      @ indices in 'a_matrix' @
    for i (1,n,1);
        a_matrix[i_col[i],.] = aa[i,.];
    endfor;
    a = vec(a_matrix);
    retp(a);
endp;

/* product of elements from i up to the last */
proc MULT(i,row);
    retp(prodc(row[i:cols(row)]'));
endp;

/* sum of elements from i up to the last */
proc SUM(i,row);
    retp(sumc(row[i:cols(row)]'));
endp;

/* calculation of the matrix with rows - all possible
   combinations of the Chebyshev zeros across state vars */
proc ZERO_VEC;
    local res,n,last,i;
    n = _prdegs;
    last = cols(n);         /* = number of state variables; max(last)=3 */
    res = ZERO_CH(n[last],last)';   /* n[last] x 1 */
    for i (1,2,1);
      if last>i;    
        res = vec(ones(rows(res),1)*ZERO_CH(n[last-i],last-i))~(ones(n[last-i],1).*.res);
        /* the 2nd component is vertical concatenation of matrix 'res' 'n[last-i]' times */
      endif;
    endfor;
    retp(res);          @ _prnpar x _prnst @
endp;

/* procedure for determination of Chebyshev zeros */
proc ZERO_CH(n,st);     @ n - # of degrees of appr, st - no. of state var @
local z,zz,order;
    z = zeros(1,n);
    zz = zeros(1,n);
    if n==2;
        if _prvar==1;
            z[1]=OUT_CH(1/sqrt(2),st);    /* a_1 */
            z[2]=OUT_CH(-1/sqrt(2),st);
        elseif _prvar==2;
            z[1]=OUT_CH(0,st);            /* a_2 */
            z[2]=OUT_CH(1/sqrt(2),st);
        else;
            z[1]=OUT_CH(0,st);            /* a_3 */
            z[2]=OUT_CH(-1/sqrt(2),st);
        endif;
    else;
        order = seqa(1,1,n)';                       @ = (1 2 ... n) @
        zz[1,.] = cos((2*order[1,.]-1)*pi/(2*n));   @ zeros @
        z = OUT_CH_ST(zz,st);                       @ transformed zeros @
    endif;
    retp(z);
endp;

/* transformation from Chebyshev approximation */
fn  OUT_CH(z,i)=(z+1)*(_prxint[2,i]-_prxint[1,i])/2+_prxint[1,i];

/* vector transformation from Chebyshev approximation across all state vars */
proc OUT_CHM(z);
    local y;
    y=zeros(1,_prnst);
    y[1,.]=(z[1,.]+1).*(_prxint[2,.]-_prxint[1,.])/2+_prxint[1,.];
    retp(y);
endp;

/* vector transformation from Chebyshev approximation 
   across all zeros for one (ith) state var; n = _prdegs[i] */
proc OUT_CH_ST(z,i);
    local y,n;
    n = cols(z);    @ z - row @
    y=zeros(1,n);
    y[1,.]=(z[1,.]+1).*(_prxint[2,i]-_prxint[1,i])/2+_prxint[1,i];
    retp(y);
endp;

/* transformation to Chebyshev approximation */
fn  IN_CH(z,i)=2*(z-_prxint[1,i])/(_prxint[2,i]-_prxint[1,i])-1;

/* vector transformation for Chebyshev approximation */
proc IN_CHM(z);
    local y;
    y=zeros(1,_prnst);
    y[1,.]=2*(z[1,.]-_prxint[1,.])./(_prxint[2,.]-_prxint[1,.])-1;
    retp(y);
endp;

/* computation of residual operator */
proc RESID(a);
    local i,b,bb;
    b = zeros(_prnpar*_prncon,1);
    bb = zeros(_prncon,_prnpar);
    for i (1,_prnpar,1);
        bb[.,i] = _FRES(_przero[i,.],a);   @ column @
    endfor;
    b = vec(bb);
    retp(b);
endp;

/* decomposition of parameter vector; result is (_prnpar,_prncon);
   this procedure is reversial for VEC(.) 
proc DECOMP(a);
    retp(reshape(a,_prncon,_prnpar)');  
endp;*/

/* procedure for modification of general power function */
proc __GPOW(x,pow);
    local y,a;
    a=zeros(3,1);
    if x>_prxeps;
        y=x^pow;
    else;
        a[1]=_prxeps^pow*(2-3*pow+pow*pow)/2;
        a[2]=pow*(2-pow)*_prxeps^(pow-1);
        a[3]=pow*(pow-1)*_prxeps^(pow-2)/2;
        y=a[1]+a[2]*x+a[3]*x*x;
    endif;
    retp(y);
endp;

/* computation of Chebyshev approximation */
/* horizontal direct product (*~) and tensor product (.*.) 
   give equal results for row vectors */
proc __APROX(x,a);              @ x - row, a - column(_prnpar,1) @
    local z,n,Ch_z,y;
    z = IN_CHM(x);            @ row (1x_prnst) @
    n = _prdegs;
    if (_prnst==1)or(_prnst==2)or(_prnst==3);
        Ch_z = CHP_ROW(z[1],n[1],0);
        if (_prnst==2)or(_prnst==3);
            Ch_z = Ch_z*~CHP_ROW(z[2],n[2],0);
        endif;
        if (_prnst==3);
            Ch_z = Ch_z*~CHP_ROW(z[3],n[3],0);
        endif;
    endif;
    y = Ch_z*a;                            
    retp(y);            @ 1x1 @
endp;

/* computation of the first derivation of Chebyshev approximation */
/* horizontal direct product (*~) and tensor product (.*.) 
   give equal results for row vectors */
proc __DAPROX(x,a,nder);              @ x - row, a - column, nder - scalar @
    local z,n,Ch_z,y,Ch_zz;
    z = IN_CHM(x);            @ row (1x_prnst) @
    n = _prdegs;
    if (_prnst==1)or(_prnst==2)or(_prnst==3);
        if nder==1;
            Ch_z = CHP_ROW(z[1],n[1],1);
        endif;    
        if (_prnst==2)or(_prnst==3);
            if nder==1;
                Ch_z = Ch_z*~CHP_ROW(z[2],n[2],0);
            else;
                Ch_zz = CHP_ROW(z[1],n[1],0);
                Ch_z = Ch_zz*~CHP_ROW(z[2],n[2],1);
            endif;
        endif;
        if (_prnst==3);
            if (nder==1)or(nder==2);
                Ch_z = Ch_z*~CHP_ROW(z[3],n[3],0);
            else;
                Ch_z = Ch_zz*~CHP_ROW(z[2],n[2],0)*~CHP_ROW(z[3],n[3],1);
            endif;    
        endif;
    endif;
    y = Ch_z*a;        
    y=y*_prtrd[nder];
    retp(y);            @ 1x1 @
endp;

/* computes matrix of Chebyshev polynomials coefficients - once and forever */
proc CH_COEF_INIT(n);            
    local i,coeff;          
    coeff = zeros(n,n);
    coeff[1,n] = 1;
    coeff[2,n-1] = 1;
    for i (3,n,1);
        coeff[i,.] = 2*shiftr(coeff[i-1,.],-1,0) - coeff[i-2,.];
    endfor;
    retp(coeff);   @ nxn @
endp;        

/* extracts matrix of Chebyshev polynomials coefficients - from the previously computed */
proc CH_COEF(m);  
    local n,coeff;
    n = cols(_ch_coefs);               @ = rows(_ch_coefs) @
    coeff = _ch_coefs[1:m,n-m+1:n];    @ the upper right square @
    retp(coeff);   @ nxn @
endp; 

/* computes matrix of coefficients for Chebyshev polynomials _derivatives_ */
proc DER_CH_COEF(n);            
    local der_coeff,coeff,rev_order;
    coeff = CH_COEF(n);
    rev_order = seqa(n,-1,n)';                  /* = (n n-1 ... 2 1); */
    der_coeff = shiftr(coeff,ones(n,1),zeros(n,1)).*(ones(n,1)*rev_order);
    retp(der_coeff);   @ nxn @
endp;

/* computes Chebyshev polynomials or its derivatives */
/* in x as a vector of values from T(0) to T(n-1) */
/* der = 1 - if it is derivatives, for simple polynomials der = 0 */
proc CHP_ROW(x,n,der);
    local coeff,pow_x,i,res;
    if der == 0;
        coeff = CH_COEF(n);
    else;
        coeff = DER_CH_COEF(n);    
    endif;
    pow_x = ones(n,1);
    for i (n-1,1,-1);           @ by recursion since so needs less memory @
        pow_x[i] = x*pow_x[i+1];
    endfor;
    res = coeff*pow_x;
    retp(res');      @ 1xn @
endp;

/* computation of coefficients of ordinary polynomial approximation */
/* TRANSFORMATION FROM THE INTERNAL REPRESENTATION INTO EXTERNAL */
proc COEF(a);                   @ a = (_prncon*_prnpar x 1); t = (_prnpar x _prnpar); @
    local i,tenprd,last,a_cols,res;
    last = _prnst;              /* = number of state variables; max(last)=3 */
    tenprd  = TRANS_CH(last);   
    for i (1,2,1);
        if last>i;
           tenprd = TRANS_CH(last-i).*.tenprd; 
        endif;
    endfor;
    a_cols = reshape(a,_prncon,_prnpar)'; /* every column ~ different control var */
    res = a_cols'tenprd;
    res = rev(res')';       /* reversing columns in a matrix */
    retp(res);              @ _prncon x _prnpar @
endp;

/* computation of one dimensional transformed Chebyshev polynomial */
proc TRANS_CH(ind);
    local n,br,order,up_j_i,up_n_j,up_n_i_fact,res;
    n = _prdegs[ind];
    br = zeros(n,n);
    order  = seqa(1,1,n);        @ =(1 2 ... n)' @
    up_j_i = shiftr(ones(n,1)*order',order,zeros(n,1));
    up_n_j = upmat(n*ones(n,n)-ones(n,1)*order');
    up_n_i_fact = upmat((n*ones(n,n)-order*ones(n,1)')!);  /* provides zeros in the lower part of 'br' */
    br = (up_n_i_fact)./((up_n_j!).*(up_j_i!)).*(_prtrd[ind]^up_n_j).*(_prtre[ind]^up_j_i);  /* D_x */
    res = CH_COEF(n)*br;
    retp(res);         @ nxn @              
endp;

/* the method of steepest descent */
proc STEEP(f,x0);
    local f:proc;
    local iter,grad,w,y,x,m;
    iter    = 0;
    grad    = 2*ones(rows(x0),1);
    x       = x0;
    "steep";
    do while (abs(maxc(grad)) > _prtols) AND (iter <= _priters);
        W       = Gradp(&f,x);
        y       = f(x);
        m       = y'W*W'y/((W*W'y)'W*W'y);
        grad    = W'y;
        x       = x - m*grad;
        iter    = iter + 1;
        "iter=" iter;
        "w=" w;
        "y=" y;
        "grad=" grad;
        "m=" m;
        "x=" x;
    endo;
    if iter > _priters; _prretc = 1; endif;
    retp(x);
endp;


/* Newton method */
proc NEWTON(f,x0);
    local f:proc;
    local iter,last_x,w,y,x,sequence,parlist,reslist,fmtstr;
    iter    = 0;
    last_x  = ones(rows(x0),1);
/*    "newton"; */
    x       = x0;
    do while (abs(sumc(x-last_x)) > _prtoln) AND (iter <= _pritern);
        y       = f(x);
        if __output;
            print "Iteration #" ftos(iter,"*.*lf",1,0);
            sequence    = seqa(1,1,_prnpar*_prncon);
            parlist     = zeros(_prnpar*_prncon,1) $+ "A[" $+
                         ftocv(sequence,floor(log(_prnpar*_prncon)+1),0) $+ "]";
            reslist     = zeros(_prnpar*_prncon,1) $+ "R[" $+
                         ftocv(sequence,floor(log(_prnpar*_prncon)+1),0) $+ "]";
            fmtstr      = {"-*.*lf" 15 8,
                            "*.*lf" 18 8,
                            "*.*lf" 22 8,
                            "*.*lf" 22 8};
            call printfm(parlist~x~reslist~y,0~1~0~1,fmtstr);
            print;
        endif;
        W       = Gradp(&f,x);
        y       = f(x);
        last_x  = x;
        x       = x - inv(w)*y;
        iter    = iter + 1;
/*        "w=" w;
        "y=" y;
        "x=" x;*/
    endo;
    if iter > _pritern; _prretc = 1; endif;
/*    f(x); */
    retp(x);
endp;

/* procedure for computing Gauss-Hermite quadrature */
proc GH_quad(f,x,a,n);
    local f:proc;
    local z,nmax,W,s,i,ff;
    nmax    = 10;
    z       = zeros(nmax,1);    /* abscissas */
    W       = z;                /* weights */
    if n==2;
        z[1]=.707106781186548;
        z[2]=-z[1];
        W[1]=.8862269254528;
        W[2]=W[1];
    elseif n==3;
        z[2]=1.224744871391589;
        z[3]=-z[2];
        W[1]=1.1816359006037;
        W[2]=.2954089751509;
        W[3]=W[2];
    elseif n==4;
        z[1]=.524647623275290;
        z[2]=-z[1];
        z[3]=1.650680123885785;
        z[4]=-z[3];
        W[1]=.8049140900055;
        W[2]=W[1];
        W[3]=.08131283544725;
        W[4]=W[3];
    elseif n==5;
        z[2]=.958572464613819;
        z[3]=-z[2];
        z[4]=2.020182870456086;
        z[5]=-z[4];
        W[1]=.9453087204829;
        W[2]=.3936193231522;
        W[3]=W[2];
        W[4]=.01995324205905;
        W[5]=W[4];
    elseif n==6;
        z[1]=.436077411927617;
        z[2]=-z[1];
        z[3]=1.335849074013697;
        z[4]=-z[3];
        z[5]=2.350604973674492;
        z[6]=-z[5];
        W[1]=.7246295952244;
        W[2]=W[1];
        W[3]=.1570673203229;
        W[4]=W[3];
        W[5]=.004530009905509;
        W[6]=W[5];
    else;
    endif;
    s = 0;
    i = 1;
    do while i<=n;
        ff = f(z[i],x,a);
        s  = s + W[i]*ff;
        i  = i + 1;
    endo;
    retp(s);
endp;

/*
**  PROC PROJSET
**
**  FORMAT
**          PROJSET;
**
**  INPUT
**          none
**
**  OUTPUT
**          none
**
**  REMARK:
**          PROJSET should be called before any calling of PROJEC. By setting
**          all global variables to their default values before calling PROJEC,
**          any unwanted side affects are avoided.
**
*/

proc (0) = PROJSET;

/*============================================================================*/
  _prmeth = 1;      /* indicator of solution strategy - Newton-Raphson method */
  _prnst  = 1;      /* number of state variables */
  _prncon = 1;      /* number of control variables */
  _prdegs = 2;      /* degrees of approximation */
  _ch_coefs = 0;     /* coefficients of Chebyshev polynomials */
  _prvar  = 2;      /* variant of zeroes for 2nd degree approximation */
  _prxint = {.5, 1.5};  /* range of state variables */
  _prinit = 0;      /* type of input of initial parameter values
                       (linear estimation) */
  _prldfn = "";     /* name of file with initial parameters values */
  _prsave = 0;      /* indicator of saving results on disk (no saving) */
  _prsvfn = "";     /* file name for saving computed parameters */
  _pr_ver = { 1.00, 1 };
  _pritern= 50;     /* max. number of iterations for Newton method */
  _priters= 50;     /* max. number of iterations for steepest descent method */
  _prtoln = 1e-10;  /* tolerance limit of convergence for Newton method */
  _prtols = .5;     /* tolerance of convergence for steepest descent method */
  _prretc = 0;      /* return code */
  _pr_ver = { 1.1, 1 };
  gausset;
/*============================================================================*/
endp;