/* mfem_cox.src


Written by Dan Powers (12/17/98),
Department of Sociology
University of Texas at Austin
E-mail: dpowers@mail.la.utexas.edu

This code is written and submitted for public, non-commercial use.  If you use
this code extensively in a research project, please cite the author.


mfEM_cox.src is a GAUSS program to estimate the multiplicative (gamma)
frailty model for clustered data of Guo and Rodriguez (JASA, 1992,
87:969-967) using the EM-Algorithm.  This routine differs from Guo and
Rodriguez insofar as the baseline hazard is unspecified (semi-parametric)
rather than piecewise constant. The code is an extensive rewrite and
modification of Naryan Satry's S-Plus procedure (repwe.s) for estimating
Guo and Rodiguez's model.

The data must contain a unique id number to identify clusters (i.e.,
families, sibships, communities, or individuals etc).  This program is
designed to estimate the variance in a gamma distributed cluster-level
random effect.  This program also works with individual-level data, in
which case all individuals are members of clusters of size one, i.e.,
cluster_id = individual_id.  If clustering is not an issue, then a simple
modification will reduce the computational burden of accumulating the
integrated hazard over individuals in clusters.  However, the program is
fast without these modifications.

An estimate of the variance in the cluster-level random effect is
given by 1/alpha (phi). Model results may be compared against the
conventional Cox model, which this program estimates after the
multiplicative frailty model.  The user can request that empirical
Bayes estimates of the cluster-level frailty and their variances be
written to an output file, for use in identifying high and low risk
clusters.

The program uses the EM accelerator used by Guo and Rodriguez, in
addition to the one proposed by Sastry (JASA, 1997, 92:426-435).  The
computational burden lies in the loops used to collect and append the
cluster-level information.

The underlying hazard model is left unspecified (A Cox proportional hazard
model) as suggested by Lancaster (1990).  It is possible to modify the
coxfit routine to accomodate nonproportional covariate effects (please
inform the author of any progress along these lines.)

Note that a calling program is used and this set of procedures is
"included".  I thank Paul Fackler for his set of gamma routines (used
below).  Note that in the calling program, maxit shoul be large (250) and
maxit1 should be small (less than 50).

REFERENCES:

Guo, G., and G. Rodriguez (1992) "Estimating a Multivariate Proportional
Hazards Model for Clustered Data using the EM-Algorithm, With an
Application to Child Survival in Guatemala,'' _Journal of the American
Statistical Association_ 87:969-976.

Lancaster, T. (1990) _The Econometric Analysis of Transition Data_
Cambridge:Cambridge Univ. Press.

Sastry, N. (1997) "A Nested Frailty Model for Survival Data, With an
Application to the Study of Child Survival in Northeast Brazil," _Journal
of the American Statistical Association_ 92:426-435.


*/


/*  Build Model */

ni=indcv(model,dlist);
 if ismiss(ni);
 "  MODEL VARIABLE NAME DOES NOT MATCH NAME IN DATA LIST";?;
stop;
 endif;
    k=1;
      do until k > rows(ni);
        if k==1;
          x=ndta[.,ni[k]];
           elseif k >1;
            x=x~ndta[.,ni[k]];
           endif;
         k=k+1;
      endo;

if rows(model) /= cols(x);

/* don't know how this could happen but.... */

"NUMBER OF NAMES DOES NOT CORRESPOND TO NUMBER OF VARIABLES";?;

stop;
endif;

/* initialization */
tot_ind_cl= zeros(n,1);
tot_ev_cl = zeros(n,1);
tot_event = zeros(m,1);
tot_indiv = zeros(m,1);
tot_e  = zeros(m,1);
u=ones(m,1);
tot_Ihaz_cl = zeros(n,1);
tot_Ihaz = zeros(m,1);
sum_event=0;
sum_Ihaz=0;
sum_indiv=0;
muhat1 = zeros(m,1);
muhat  = zeros(n,1);
xihat = zeros(n,1);
haz_b = zeros(cols(x),1);
old_b = zeros(cols(x),1);
xpwx = eye(cols(x));
pred_haz = zeros(m,1);
offs=ones(m,1);

"NUMBER OF CLUSTERS :                     ";; n;
"NUMBER OF (PSEUDO) OBSERVATIONS:         ";; m;?;


/* check for identified model */

        trap 1;
           vcinit = invpd(x'x);
        trap 0;
           if scalerr(vcinit);
?;

"MODEL IS NOT IDENTIFIED DUE TO LINEAR COMBINATIONS IN MODEL MATRIX";?;
"                      CHECK  MODEL SPECIFICATION                  ";?;
STOP;
endif;

format (12,3);
proc timer;
  local t,secs;
    t=time;
  secs=t[1]*3600+t[2]*60+t[3]+t[4]/100;
retp(secs);
endp;

tstart=timer;

/* loop thru sorted id numbers to get jj index for temp vars.  The index
   jj corresponds to Guo and Rodriguez // n_i //, the number of
   individuals in ith cluster.
*/


ii=0;
 j=0;
jj=zeros(n,1);
i=1;
  do until i > m;
        if i eq 1;    sum_event=d[i];sum_indiv=u[i];
           ii=1;
            j=1;
        else;
          if id[i] ne id[i-1];   sum_event=d[i]; sum_indiv=u[i];
           ii = ii+1;
            j=1;
        elseif id[i] eq id[i-1]; sum_event=sum_event+d[i];
                                 sum_indiv=sum_indiv+u[i];
            j=j+1;

           ii=ii;


        endif;

       endif;
       tot_ev_cl[ii]=sum_event;
       tot_ind_cl[ii]=sum_indiv;
       jj[ii]=j;
      i=i+1;
endo;

mc=maxc(jj);
"MAXIMUM NUMBER OF CLUSTER CONTRIBUTIONS: ";;mc;?;

/* initialize temp vars with -1 (for missing) */

tot_Ih_tmp=-ones(n,mc);
tot_Ev_tmp=-ones(n,mc);
tot_In_tmp=-ones(n,mc);
        mx=-ones(n,mc);

/* accumulate cluster level event counts */

i=1;
  do until i > n;
    /*   tot_ev_cl[i]=sumc( d.*(id.==c[i]) ); */
       j=1;
        do until j > jj[i];
          tot_ev_tmp[i,j]=tot_ev_cl[i];
          tot_in_tmp[i,j]=tot_ind_cl[i];
        j=j+1;
       endo;
  i=i+1;
endo;

/* vectorized -with missing deleted */

tot_event = packr(vecr(miss(tot_ev_tmp,mx)));
tot_indiv = packr(vecr(miss(tot_in_tmp,mx)));
clear tot_ev_tmp;
clear tot_in_tmp;

/* begin estimation loop  */

/* Initial Cox Fit */



{Ihaz,pred_haz,xpwx,haz_b} = coxfit(haz_b,t,d,x,offs,.0001,10);

alpha_new = 1;
iter = 0;
phidif = 1;

/* main iteration loop */

"Starting Alpha:";;alpha_st;

do while abs(phidif) > tol  and iter <= maxit1;
   if iter > 0;
         old_b = haz_b;
         muhat1 = (alpha_new + tot_event)./(alpha_new + tot_Ihaz);
   endif;

    alpha = (iter == 0) * alpha_st  + (iter > 0) * alpha_new;

 if iter >= 1;

   offs = muhat1;

 endif;

{Ihaz,pred_haz,xpwx,haz_b} = coxfit(haz_b,t,d,x,offs,tol,maxit1);

/* "phidif";;phidif; */

logl = n*(alpha.*ln(alpha)-ln(gamma(alpha)))
         + sumc(d.*pred_haz)
         + sumc(ln(gamma(alpha+tot_ev_cl)))
         - sumc((alpha+tot_ev_cl).*ln(alpha+tot_Ihaz_cl));

"LOG LIKELIHOOD   :";;logl;


/* accumulate cluster-level integrated hazard */

i=1;
  do until i > m;
        if i eq 1;    sum_Ihaz=Ihaz[i];
           ii=1;
            j=1;
        else;
          if id[i] ne id[i-1];   sum_Ihaz=Ihaz[i];
           ii = ii+1;
            j=1;
        elseif id[i] eq id[i-1]; sum_Ihaz=sum_Ihaz+Ihaz[i];
            j=j+1;

           ii=ii;
        endif;
       endif;

  tot_Ihaz_cl[ii]=sum_Ihaz;

       /*  ii;;j; */
       jj[ii]=j;
      i=i+1;
endo;



i=1;
  do until i > n;
      /* tot_Ihaz_cl[i] = sumc(Ihaz.*(id .== c[i])); */
            j=1;
           do until j > jj[i];
             tot_Ih_tmp[i,j] = tot_Ihaz_cl[i];
          j=j+1;
    endo;
  i=i+1;
endo;

/* vectorized with missing values removed */

tot_Ihaz = packr(vecr(miss(tot_Ih_tmp,mx)));


xihat = digamma(alpha + tot_ev_cl) - ln(alpha + tot_Ihaz_cl);
muhat = (alpha + tot_ev_cl)./(alpha + tot_Ihaz_cl);
ximudif = sumc(xihat - muhat);

{alpha_new} = Newt_Raph(alpha,ximudif,n,tol,maxit,0);

phidif3 = 1;
  iter3 = 0;

/* secondary loop */

do while abs(phidif3) > tol and iter3 <= maxit;
    xihat = digamma(alpha_new + tot_ev_cl) - ln(alpha_new + tot_Ihaz_cl);
    muhat = (alpha_new + tot_ev_cl)./(alpha_new + tot_Ihaz_cl);
    ximudif = sumc(xihat - muhat);
    alpha_old = alpha_new;

    {alpha_new} = Newt_Raph(alpha_old,ximudif,n,tol,maxit,0);

  phidif3 = 1/alpha_new - 1/alpha_old;
/*  iter3;;"phidif3";;phidif3; */
 iter3 = iter3 + 1;
endo; /* secondary loop: iter3 */

phidif = 1/alpha - 1/alpha_new;

iter = iter + 1;

if print_it == 1;

format /rd (12,0);
"EM-ITERATION     :";;iter;
"NR-ITERATIONS    :";;iter3;
format (12,6);
"ALPHA            :";;alpha;
"PHI              :";;1/alpha;
"TOLERANCE        :";;phidif;


endif;

endo; /* main loop: iter */

clear tot_ev_tmp;

tfinish = timer;

?;
"MODEL CONVERGED WITH TOLERANCE ";;phidif;;" IN ";; tfinish-tstart;;" SECONDS ";
?;

?;
"            C O M P U T I N G   S T A N D A R D   E R R O R S";
?;

/* variance estimation */

xihat = digamma(alpha + tot_ev_cl) - ln(alpha + tot_Ihaz_cl);
muhat = (alpha + tot_ev_cl)./(alpha + tot_Ihaz_cl);
d2Qdtheta2=xpwx;
p=rows(d2Qdtheta2);
EwD = (d2Qdtheta2|zeros(p,1)')~(zeros(p+1,1));
EwD[p+1,p+1] = -n * (1/alpha - trigamma(alpha));
var_w = muhat./(alpha + tot_Ihaz_cl);
var_logw = trigamma(alpha + tot_ev_cl);
cov_w_logw = muhat .* (digamma(alpha + tot_ev_cl + 1)
                          - digamma(alpha + tot_ev_cl));
var_logww = var_logw + var_w - 2.*cov_w_logw;
cov_w_logww = cov_w_logw - var_w;
A=zeros(p,1);
 dhaz_dbeta = x .* Ihaz;
 dhaz_dgam  = Ihaz;


/* loop over clusters and coefficient index */

i = 1;
 do until i > n;
 j = 1;
   do until j > p;
     A[j]=sumc(dhaz_dbeta[.,j].*(id.==c[i]));
   j=j+1;
endo;


/* build information matrix */


Ti = (-A~zeros(p,1))|(0~1);

VarWi = (var_w[i]~cov_w_logww[i])|(cov_w_logww[i]~var_logww[i]);
VarUi = Ti*VarWi*Ti';

if i == 1;
  VarU = VarUi;
else;
  VarU = VarU + VarUi;
endif;

i = i + 1;
endo;

/* either of these will work */

@
u=EwD-VarU;
var=sqrt(diag(invpd(u)));
@
{uu,dd,vv} = svdcusv(EwD-VarU);
var = sqrt(diag( (uu*  (vv' * 1/dd)  )' ));


logL=n*(alpha * ln(alpha)-ln(gamma(alpha))) + sumc(d.*pred_haz) +
sumc(ln(gamma(alpha+tot_ev_cl))) -
sumc((alpha + tot_ev_cl).*ln(alpha+tot_Ihaz_cl));

"Model Log_Likelihood";;logL;

phi    = 1/alpha;
se_phi = sqrt(alpha^(-4).*var[p+1]^2);
t_phi  = phi/se_phi;
tau    = phi/(2+phi);

"1/alpha (phi)    ";;1/alpha;
"se(1/alpha)      ";;sqrt(alpha^(-4).*var[p+1]^2);
"t(1/alpha)       ";;(1/alpha)./sqrt(alpha^(-4).*var[p+1]^2);
tv=haz_b./var[1:p];

?;
"                ESTIMATES FROM MULTIPLICATIVE FRAILTY MODEL
__________________________________________________________________________";
i=1;
"     VARIABLE    ESTIMATE     ST.ERROR      T-VAL      RISK RATIO
__________________________________________________________________________";?;
  do until i > p;
    $model[i];;haz_b[i];;var[i];;tv[i];;exp(haz_b[i]);
  i=i+1;
endo;

tb=tv;
hb=haz_b;
hb=zeros(p,1);
vb=var[1:p];

{Ihaz,pred_haz,xpwx,hb} = coxfit(hb,t,d,x,ones(rows(t),1),tol,maxit);

vb=sqrt(diag(invpd(xpwx)));
tb = hb./vb;

?;
"                       ESTIMATES FROM STANDARD MODEL
__________________________________________________________________________";
i=1;
"     VARIABLE    ESTIMATE     ST.ERROR      T-VAL      RISK RATIO
__________________________________________________________________________";?;
  do until i > p;
    $model[i];;hb[i];;vb[i];;tb[i];;exp(hb[i]);
  i=i+1;
endo;

/* output file=muhat.dat reset; muhat; */
if write_frail==1;
pause(5);;".....writing frailities to";$ffname;
output file=^ffname reset;
i=1;
  do until i > n;
c[i];;tot_ind_cl[i];;muhat[i];;var_w[i];
i=i+1;
endo;
endif;


/* procedures */

/*  Newton-Raphson Step for Alpha: Adapted from Naryan Sastray's S-Plus
    Program
*/

proc(1)=Newt_Raph(alpha_old,ximudif,n,tol,maxit,itprint);
   local dQdalpha, d2Qdalpha2, alpha, alpha_new, phidif,iter;

/*  Newton-Raphson Interations */

iter  = 0;
phidif = 1;
alpha = alpha_old;
do while abs(phidif) > tol and iter <= maxit;
  dQdalpha  = ximudif + n*(1 + ln(alpha) - digamma(alpha));
  d2Qdalpha2 = n * (1 / alpha - trigamma(alpha));
  alpha_new = alpha - dQdalpha/d2Qdalpha2;
  phidif = 1/alpha - 1/alpha_new;

if itprint == 1;

"     Alpha       PhiDif
____________________________";
alpha~phidif;

endif;

  iter = iter + 1;
alpha = alpha_new;
endo;
retp(alpha);
endp;


/* Proc COXFIT:  Gauss procedure to estimate the Cox proportional
                 hazards model using Breslow-Peto adjustment for
                 ties.  This procedure can be called from an EM
                 routine to estimate the variance of a gamma
                 distributed random effect.  The frailties are
                 estimated in the first part of the M-step and
                 treated as known in the second part of the
                 M-step. This procedure just does the coxfit
                 (second part of M-Step).

Inputs:    b  -  starting coeficient vector (kx1)
           t  -  waiting times              (nx1)
           d  -  event=1/censoring=0        (nx1)
           x  -  covariates                 (nxk)
        offs  -  offset                     (nx1)
        tol   -  convergence tolerance      (1x1)
      maxit   -  maximum iterations         (1x1)
       init   -  =1 if intitializing        (1x1)
                 =0 if subsequent call
Output
       Mod    -  column[1]=integrated hazards (Aalen's Estimate)
                 for use in computing gamma frailties.  (nx1)
              -  column[2]=original id number (after a few sort
                 procedures we need to check this).     (nx1)
         b    -  estimates                              (kx1)
         H    -  information                            (kxk)

*/

proc(4)=coxfit(b,t,d,x,offs,tol,maxit);


local dta, ndta, n_obs, n_fail, n_uni, id, i,
      ii, j, jj, k, tf, tf_u, df, xf, sx, dd, iter, db,
      new_b, U, H, lam, sum_d, sum_x, mc, inc, idx, R,
      sum_id,ntf,rsdta,tidx, tlam, haz, Ihaz, risk, LogL,
      modout;

/* first concatenate the input data */

dta = t~d~x~offs;
  k = cols(x);
 id = seqa(1,1,rows(dta));

/* reverse sort by failure-time for computational ease */

ndta = rev(sortc(dta,1));
   t = ndta[.,1];
   d = ndta[.,2];
   x = ndta[.,3:3+(k-1)];
offs = ndta[.,3+k];
n_obs  = rows(ndta);  /* number of observations    */
n_fail = sumc(d);     /* number of failures/events */


/*
Some calculations work only over failure times.
*/

tf     = selif(t,d.==1);   /* failure times                            */
df     = selif(d,d.==1);   /* ones                                     */
xf     = selif(x,d.==1);   /* covariate vector                         */
tf_u   = rev(unique(tf,1));/* unique failure times                     */
n_uni  = rows(tf_u);


/*    We need some quantities based on the
      unique failure times, such as:
dd  - number of ties at each unique failure time
sx  - the sum of covariate vector for tied obs.
*/


ii=0;
 j=0;
sx =zeros(n_uni,k);
dd =zeros(n_uni,1);
jj =zeros(n_uni,1);
sum_x  = zeros(1,k);
sum_d  = 0;

i=1;
  do until i > n_fail;
        if i eq 1;
  sum_x  = xf[1,.];
  sum_d  = df[1];
           ii=1;
            j=1;
        else;
          if tf[i] ne tf[i-1];
sum_x  = xf[i,.];
sum_d  = df[i];
           ii = ii+1;
            j=1;
        elseif tf[i] eq tf[i-1];
sum_x  = sum_x+xf[i,.];
sum_d  = sum_d+df[i];
            j=j+1;
           ii=ii;
        endif;
       endif;
       jj[ii]=j;
sx[ii,.] = sum_x;
dd[ii]   = sum_d;
i=i+1;
endo;

/* Get Risk Set Indicies */

idx=zeros(n_uni,1);
 i=1;
  do until i > n_uni;
     idx[i]=sumc(ones(n_obs,1).*(t .>= tf_u[i]));
  i=i+1;
endo;


mc=maxc(jj);
"NUMBER OF OBSERVATIONS                        : ";;n_obs;?;
"NUMBER OF FAILURES                            : ";;n_fail;?;
"PERCENT CENSORED                              ; ";;1-n_fail/n_obs;?;
"NUMBER OF UNIQUE FAILURE TIMES                : ";;n_uni;?;
"MAXIMUM NUMBER OF TIES PER UNIQUE FAILURE TIME: ";;mc;?;
"TOTAL NUMBER OF TIED FAILURE TIMES            : ";;
sumc(dd.*(dd.>1));?;
"PROPORTION OF TIES TO FAILURES                : ";;
sumc(dd.*(dd.>1))./n_fail;?;


/* if many ties call FIXUP before estimating this model */

/* NEWTON-STEP */

 db  = 1;
iter = 0;
do until abs(db) < tol or iter > maxit;

{lam,U,H,R} = CALC(b,x,sx,offs,k,n_uni,dd,idx);

    new_b = b  + solpd(U,H);
       db = new_b-b;
        b = new_b;
iter=iter+1;
endo;

"COX MODEL CONVERGED WITH TOLERANCE ";; format /rd (5,5); meanc(abs(db));;
"IN ";;format /rd (5,0);iter;;"ITERATIONS";
format /rd (12,8);

/*
   These are quantites to be output along with results.
   We need to get the times into ascending order to make
   this come out right.

   Sort by time and bring id and x along.  Next, sort by
  unique time, bring along estimated hazard.
*/

tidx  = sortc(t~d~id~x,1);
tlam  = sortc(tf_u~lam,1);

/*
   Get quantities needed to compute Aalen's estimate of the
   Integrated Hazard. We will also need this to compute the
   likelihood function for this model.
*/

risk  = exp(tidx[.,4:4+k-1]*b);

/* Loop over ALL observation to get Hazard and Cumulative Hazard */

Ihaz=zeros(n_obs,1);
 haz=zeros(n_obs,1);

i=1;
   do until i > n_obs;
          haz[i] = risk[i].*maxc(tlam[.,2].*(tlam[.,1] .== tidx[i,1]));
         Ihaz[i] = risk[i].*sumc(tlam[.,2].*(tlam[.,1] .<= tidx[i,1]));
   i=i+1;
endo;

/*
   Bookkeeping - rearrange by original id number so that Ihaz[i]
                 is matched up with the right observation in the
                 original data.
*/

LogL = sumc(sx*b - dd.*ln(R));

modout=sortc(Ihaz~haz~id,3);

"(Partial) Log-Likelihood ";;LogL;

retp(modout[.,1],modout[.,2],H,b);
endp;



/* PROC CALC:  risk-set loop to do main work */


proc(4)=CALC(b,x,sx,offs,npar,n_uni,dd,idx);

  local i,ii,j,xi, A, rf, rfx, rfxpx, xpx,
        rr,rx, rxx, U ,H ,lambda;

       xi = exp(x*b + ln(offs) );
    rf    = cumsumc(xi);
    rfx   = cumsumc(x.*xi);
    rfxpx = cumsumc(x*~x.*xi);


/* score and hessian */

rx=zeros(n_uni,npar);
rxx=zeros(n_uni,npar*npar);
rr=zeros(n_uni,1);

/* Score and Information  */

i = 1;
  do until i > n_uni;
ii=idx[i];
rx[i,.] = rfx[ii,.];
rxx[i,.]= rfxpx[ii,.];
rr[i]   = rf[ii];
i=i+1;
endo;

U = sumc(sx - dd.*rx./rr);

xpx=zeros(n_uni,npar*npar);

A=(rx./rr)*~(rx./rr);

j=1;
  do until j > npar*npar;
    xpx[.,j]=dd.*(rxx[.,j]./rr - A[.,j]);
   j=j+1;
  endo;

H=reshape(sumc(xpx),npar,npar);

/*
  Lambda is the estimated baseline hazard at each failure time -
  to be used to calculated the Integrated hazard and Log L.
*/

lambda = dd./rr;

retp(lambda,U,H,rr);
endp;

/*
   PROC fixup - replaces tied failure times with time + a small
   increment.  For use in situations with extreme ties.
   Another model would be preferred.
*/


Proc(1)=fixup(t);
  local n, nt, i, inc, ties;
n=rows(t);
nt=t;
ties=0;
   i=1;
 inc = .001;
     do until i > n;
           if t[i] ne t[i-1]; nt[i]=t[i];
       elseif t[i] eq t[i-1]; nt[i]=nt[i]+inc;
   ties=ties+1;
 endif;
   inc= inc + .001;
i=i+1;
endo;
"Ties " ;ties;
retp(nt);
endp;



/*
   Routines submitted by:
   Paul Fackler
   919-515-4535
   paul_fackler@ncsu.edu


digamma      (x) psi function (derivative of the log gamma function)

Calculates the value of the psi (digamma) function, the first
derivative of the loggamma function.  Accepts positive matrices.
>From Abromowitz and Stegun.  Uses formula 6.3.18 with recurrence formula
6.3.5. Requires workspace at least 5 times the size of X.

trigamma     (x) second derivative of the loggamma function

Calculates the value of the trigamma, the second, derivative of the
loggamma function. Accepts positive matrices. From Abromowitz and Stegun.
Uses formulas 6.4.11 and 6.4.12 with recurrence formula 6.4.6.
Each requires workspace at least 5 times the size of X.

*/

proc digamma(x);
  local p;
  x=x+6;
  p=1/(x.*x);
  p=(((0.004166666666667*p-0.003968253986254).*p+
      0.008333333333333).*p-0.083333333333333).*p;
  p=p+ln(x)-0.5/x-1/(x-1)-1/(x-2)-1/(x-3)-1/(x-4)-1/(x-5)-1/(x-6);
retp(p);
endp;


proc trigamma(x);
  local p,i;
  x=x+6;
  p=1/(x.*x);
  p=(((((0.075757575757576*p-0.033333333333333).*p+0.0238095238095238)
     .*p-0.033333333333333).*p+0.166666666666667).*p+1)./x+0.5*p;
  i=0; do until i>5;i=i+1;
    x=x-1;
    p=1/(x.*x)+p;
  endo;
  retp(p);
endp;







