/*
**                         FILE BPM13.ARC
**
** This file contains the bivariate probit model procedures.
** Part 13.
*/



/* Declaration of globals. */

declare _sfile ?=0;
declare __miss !=1;
declare _eps ?=0.001;
declare _maxit ?=30;
declare _file ?=0;
declare _cases !=0;
declare __output !=0;

/* Setting of libraries needed. */

library a:\bpm,stat,pgraph,dutil;





/*
** UPROBSEL
**
** FORMAT    {indv,info}=uprobsel(dataset,indvar,depvar,sig);
**
** INPUT      dataset  : String, the name of the dataset that
**                       contains the variables.
**            indvar   : An Lx1 character vector containing
**                       the names of the independent variables.
**                       case a unit column is appended and L-1 is incremented).
**            depvar   : A 1x1 character vector containing the name
**                       of the dependent variables.
**            sig      : Bound for chi-square test.
**
** OUTPUT     indv     : Character vector of chosen variables.
**            info     : rows(indv)x3 vector. Line j contains the log-
**                       likelihood with variables indv[1:j] included.
**                       info[j,2 3] contain the chisquared and P-value
**                       comparing the loglikelihood for variables 1 through
**                       j with the loglikelihood for variables 1 through
**                       j-1.
**
** GLOBALS    _sfile  0: Output device is file (default).
**                    1: Output device is printer.
**            _eps     : Tolerance for the gradient (default 0.001).
**            _maxit   : Upper bound for the number of iterations
**                       (default 30).
**
*/


proc (2)=  uprobsel(dataset,indvar,depvar,sig);

local nin,inlbl,nout,outlbl,outfile,indv,tel,schat0,w0,test,handle,akkoord,fmt,
      tijd1,tijd2,pending,stoppen,step,chosen,aantal,hess,lvec,mask,cnt,grad,x,
      indexi,indexd,nvec,handle,hulp1,hulp2,cnt,test,t,hulp,sortv,volgorde,llt,
      y,i,j,k,g1,g2,cnt,s,minima,maxima,rval,r,teta0,tetah,xx,yy,l,hulp,param;


tijd1=hsec;

/* Creating/resetting the output file. */

if _sfile==0;
   output file=bprobsel.out reset;
else;
   output file="lpt1" reset;
endif;




/*
** Defining the working data set, and creating a dataset with the variables
** passed in indvar and depvar.
*/

if dataset $/="";
   outfile="bprobsel";
   {nin,inlbl,nout,outlbl}=dselvar(dataset,outfile,indvar|depvar);
   open handle=^outfile for read;
   aantal=rowsf(handle);
   handle=close(handle);
else;
   aantal=rows(indvar);
endif;



/* Printing the title and the names of te variables. */

cls;
print " ";
print "            UNIVARIATE  PROBIT  :  SELECTION";
print "            ================================";
print " ";
print " ";

format 14,4;
if dataset $/="";
   print "Variable :";
   print "----------";
   print $depvar';
   print " ";
   print " ";
   print "Candidate predictors :";
   print "----------------------";
   print $indvar';
   print " ";
   print " ";
endif;
format 10,4;
print "         Criterion :    Likelihood ratio";
print "Significance level :" sig;
print "   Number of cases :" aantal;
print " ";
print " ";
print "************************************************************"\
      "*****************";
print " ";


/* Putting the variables into local matrices. */

if dataset$=="";

   x=indvar;
   y=depvar;
   let depvar="y1" "y2";
   let indvar="intercep" "x1" "x2" "x3" "x4" "x5" "x6" "x7" "x8" "x9";
   if x[.,1]==ones(rows(x),1);
      indvar=indvar[1:cols(x)];
   else;
      indvar=indvar[1:cols(x)+1];
   endif;

else;

   {indvar,indexi}=indices(outfile,indvar);
   {depvar,indexd}=indices(outfile,depvar);
   nvec=setvars(outfile);
   open handle=^outfile for read;
   do until eof(handle);
      hulp1=readr(handle,100);
      i=1;
      do while i<=rows(indvar);
         hulp2=varget(indvar[i]);
         if rows(hulp2)>1;
            hulp2=hulp2|hulp1[.,indexi[i]];
         else;
            hulp2=hulp1[.,indexi[i]];
         endif;
         cnt=varput(hulp2,indvar[i]);
         i=i+1;
      endo;

      i=1;
      do while i<=rows(depvar);
         hulp2=varget(depvar[i]);
         if rows(hulp2)>1;
            hulp2=hulp2|hulp1[.,indexd[i]];
         else;
            hulp2=hulp1[.,indexd[i]];
         endif;
         cnt=varput(hulp2,depvar[i]);
         i=i+1;
      endo;

   endo;
   handle=close(handle);
  y=varget(depvar[1]);
   x=varget(indvar[1]);
   i=2;
   do while i<=rows(indvar);
      x=x~varget(indvar[i]);
      i=i+1;
   endo;

   if not(x[.,1]==ones(rows(x),1));
      indvar="intercep"|indvar;
   endif;

   nvec=setvars(outfile);

endif;


/* Checking for faulty data and deletion of rows with missing values. */

if (rows(x)/=rows(y));
   errorlog "x en y must have the same number of rows.";
   end;
elseif (cols(y)/=1);
   errorlog "y does not stand for a 0-1 variable.";
   end;
endif;

t=x~y;
_cases=rows(t);
g1=cols(x);
g2=cols(t);
t=packr(t);
_cases=_cases|rows(t);
y=t[.,g1+1:g2];
x=t[.,1:g1];


/* Determining the shape of the tables. */

minima=minc(y);
maxima=maxc(y);
i=minima;
rval=0;
do while i<=maxima;
   if not y/=i;
      rval=rval|i;
   endif;
   i=i+1;
endo;
rval=rval[2:rows(rval)];

r=rows(rval);


/* Transforming the dependent variables into contingency tables. */

hulp=zeros(rows(y),r);
i=1;
do while i<=r;
   hulp[.,i]=y.==rval[i];
   i=i+1;
endo;

y=hulp;
clear hulp;
g1=cols(x);
g2=g1+cols(y);


/* Multiple columns sort on the independent variable(s). */

t=x~y;
sortv=seqa(1,1,g1);
t=sortmc(t,sortv);


/* Summing over equal independent variable vectors. */

i=1;
k=1;
s=zeros(rows(t),g2);
do while i<=rows(t);
   cnt=t[i,1:g1];
   j=i;
   test=1;
   do while (j<=rows(t)) and test;
      j=j+1;
      if j<=rows(t);
         test=t[i,1:g1]==t[j,1:g1];
      else;
         test=0;
      endif;
   endo;
   j=j-1;
   if j<=rows(t);
      s[k,g1+1:g2]=sumc(t[i:j,g1+1:g2])';
      s[k,1:g1]=cnt;
      k=k+1;
   endif;
   i=j+1;
endo;

s=s[1:k-1,.];



/*
** Separating dependent and independent variables and adding a column
** of ones to the matrix of independent variables if necessary.
*/

y=s[.,g1+1:g2];
x=s[.,1:g1];
if not(x[.,1]==1);
   x=ones(rows(x),1)~x;
endif;
l=cols(x);


/* Vector of initial values. */

i=1;
hulp2=zeros(rows(y),cols(y)*cols(y));
do while i<=rows(y);
   hulp1=diagrv(eye(cols(y)),y[i,.]');
   hulp2[i,.]=reshape(hulp1,1,cols(y)*cols(y));
   i=i+1;
endo;
output off;
teta0=initrs(x,hulp2,r,r);
output on;
teta0=teta0[1:l*(r-1)];
teta0=reshape(teta0,r-1,l)';


/* Scanning of individual covariates. */

output off;
tetah=vec(teta0[1,.]);
{xx,yy}=shrink(x[.,1],y);
{schat0,llt,grad,hess}=uprobrs(tetah,xx,yy);
output on;


_file=2;
tel=2;
lvec=zeros(rows(indvar)-1,3);
do while tel<=rows(indvar);
   tetah=vec(teta0[1 tel,.]);
   output off;
   {xx,yy}=shrink(x[.,1 tel],y);
   {schat0,lvec[tel-1,1],grad,hess}=uprobrs(tetah,xx,yy);
   output on;
   tel=tel+1;
endo;

tel=1;
do while tel<=rows(lvec);
   lvec[tel,2]=lvec[tel,1]-llt;
   lvec[tel,3]=cdfchic(abs(2*lvec[tel,2]),1);
   tel=tel+1;
endo;


/* Ordering the independent variables. */

volgorde=seqa(2,1,rows(indvar)-1);
lvec=-lvec~volgorde;
lvec=sortc(lvec,1);
volgorde=lvec[.,4];
lvec=-lvec[.,1:3];



print "STEP  1";
print "=======";
print " ";
print " ";
mask=0~1~1~1;
let fmt[4,3]=
"*.*s "  14  4
"*.*lf " 14  4
"*.*lf " 14  4
"*.*lf " 14  4;
chosen="intercep"~llt~0~0;
print "Variables in model :";
print "--------------------";
format 10,4;
print $chosen[.,1]';
print " ";
print " ";
print "Candidates :";
print "------------";
print " ";
print "      variable        log-lik      chisquare        p-value";
print "      --------        -------      ---------        -------";
cnt=printfm(indvar[volgorde]~lvec[.,1]~(2*lvec[.,2])~lvec[.,3],mask,fmt);
print " ";
akkoord=(lvec[1,3]<=sig);
if akkoord;
   print $indvar[volgorde[1]] "will be kept in the model.";
else;
   print "No more variables will be taken in the model.";
endif;
print " ";
print "************************************************************"\
      "*****************";
print " ";
print " ";


/* Selection procedure. */

llt=0;
indvar=indvar[1]|indvar[volgorde];
x=x[.,1]~x[.,volgorde];
teta0=teta0';
teta0=(teta0[.,1]~teta0[.,volgorde])';
indv=indvar[1 2];
chosen=chosen|(indv[2]~lvec[1,.]);
if rows(indvar)>=3;
   pending=indvar[3:rows(indvar)];
   stoppen=not akkoord;
else;
   stoppen=1;
endif;
tel=1;
step=1;


do until stoppen;
   tel=1;
   step=step+1;
   w0=lvec[1,1];
   lvec=zeros(rows(pending),3);
   format 3,0;
   print "STEP" step;
   print "=======";
   print " ";
   print "Variables in model :";
   print "--------------------";
   format 10,4;
   print $chosen[.,1]';
   print " ";
   print "Candidates :";
   print "------------";
   param=zeros(r+step-1,rows(pending));
   do while tel<=rows(pending);
      output off;
      tetah=vec(teta0[1:rows(chosen) rows(chosen)+tel,.]);
      {xx,yy}=shrink(x[.,1:rows(chosen) rows(chosen)+tel],y);
      {param[.,tel],lvec[tel,1],grad,hess}=uprobrs(tetah,xx,yy);
      output on;
      lvec[tel,2]=lvec[tel,1]-w0;
      lvec[tel,3]=cdfchic(abs(2*lvec[tel,2]),1);
      tel=tel+1;
   endo;

   param=param';
   volgorde=seqa(rows(chosen)+1,1,rows(pending));
   lvec=-lvec~volgorde;
   lvec=sortc(lvec,1);
   volgorde=lvec[.,4];
   lvec=-lvec[.,1:3];
   x=x[.,1:rows(chosen)]~x[.,volgorde];
   teta0=teta0[1:rows(chosen),.]|teta0[volgorde,.];
   volgorde=volgorde-rows(chosen);
   pending=pending[volgorde];
   param=param[volgorde,.];
   param=transoe1(param[1,.],r);
   param=reshape(param,r-1,rows(chosen)+1)';
   teta0[1:rows(param),.]=param;


   print " ";
   print "      variable        log-lik      chisquare        p-value";
   print "      --------        -------      ---------        -------";
   cnt=printfm(pending~lvec[.,1]~(2*lvec[.,2])~lvec[.,3],mask,fmt);
   print " ";
   if rows(pending)==1;

      print " ";
   endif;
   akkoord=(lvec[1,3]<=sig);
   if akkoord;
      print $pending[1] "will be kept in the model.";
   else;
      print "No more variables will be taken in the model.";
   endif;
   print " ";
   print "************************************************************"\
         "*****************";
   print " ";


   if akkoord;
      chosen=chosen|(pending[1]~lvec[1,.]);
      indv=indv|pending[1];
      if rows(pending)==1;
         stoppen=1;
      else;
         pending=pending[2:rows(pending)];
      endif;
   else;
      stoppen=1;
   endif;

endo;


output on;

print "Summary :";
print "---------";
print " ";
print "      variable        log-lik      chisquare        p-value";
print "      --------        -------      ---------        -------";
if rows(chosen)==2 and chosen[2,4]>sig;
   cnt=printfm(chosen[1,1 2],mask[1 2],fmt[1 2,.]);
else;
   cnt=printfm(chosen[1,1 2],mask[1 2],fmt[1 2,.]);
   print " ";
   cnt=printfm(chosen[2:rows(chosen),1 2]~(2*chosen[2:rows(chosen),3])~
   chosen[2:rows(chosen),4],mask,fmt);
endif;
print " ";
print " ";
print "************************************************************"\
      "*****************";


print " ";
tijd2=hsec;
tijd2=round((tijd2-tijd1)/100);
format /rd 5,0;
print "Computation time :" round(tijd2/60) "minutes" tijd2%60 "seconds.";
print " ";
print "************************************************************"\
      "*****************";
print " ";
print " ";

output off;
retp(chosen[.,1],chosen[.,2:4]);
endp;





/*
** UPROBS2
**
** FORMAT    {indv,info}=uprobs2(dataset,indvar,depvar,sig);
**
** INPUT      dataset  : String, the name of the dataset that
**                       contains the variables.
**            indvar   : An Lx1 or (L-1)x1 character vector containing
**                       the names of the independent variables.
**                       case a unit column is appended and L-1 is incremented).
**            depvar   : A 1x1 character vector containing the name
**                       of the dependent variables.
**            sig      : Bound for chi-square test.
**
** OUTPUT     indv     : Character vector of chosen variables.
**            info     : rows(indvar)x3 vector. Line j contains the log-
**                       likelihood with variables indv[1:j] included.
**                       info[j,2 3] contain the chisquared and P-value
**                       comparing the loglikelihood for variables 1 through
**                       j with the loglikelihood for variables 1 through
**                       j-1.
**
** GLOBALS    _sfile  0: Output device is file (default).
**                    1: Output device is printer.
**            _eps     : Tolerance for the gradient (default 0.001).
**            _maxit   : Upper bound for the number of iterations
**                       (default 30).
**
** REMARK     Version of UPROBSEL, for 2x1 tables.
*/


proc (2)=  uprobs2(dataset,indvar,depvar,sig);

local nin,inlbl,nout,outlbl,outfile,indv,tel,schat0,w0,test,handle,akkoord,fmt,
      tijd1,tijd2,pending,stoppen,step,chosen,aantal,hess,lvec,mask,cnt,grad,x,
      indexi,indexd,nvec,handle,hulp1,hulp2,cnt,test,t,hulp,sortv,volgorde,llt,
      y,i,j,k,g1,g2,cnt,s;


tijd1=hsec;

/* Creating/resetting the output file. */

if _sfile==0;
   output file=bprobsel.out reset;
else;
   output file="lpt1" reset;
endif;




/*
** Defining the working data set, and creating a dataset with the variables
** passed in indvar and depvar.
*/

if dataset $/="";
   outfile="bprobsel";
   {nin,inlbl,nout,outlbl}=dselvar(dataset,outfile,indvar|depvar);
   open handle=^outfile for read;
   aantal=rowsf(handle);
   handle=close(handle);
else;
   aantal=rows(indvar);
endif;



/* Printing the title and the names of te variables. */

cls;
print " ";
print "            UNIVARIATE  PROBIT  :  SELECTION";
print "            ================================";
print " ";
print " ";

format 14,4;
if dataset $/="";
   print "Variable :";
   print "----------";
   print $depvar';
   print " ";
   print " ";
   print "Candidate predictors :";
   print "----------------------";
   print $indvar';
   print " ";
   print " ";
endif;
format 10,4;
print "         Criterion :    Likelihood ratio";
print "Significance level :" sig;
print "   Number of cases :" aantal;
print " ";
print " ";
print "************************************************************"\
      "*****************";
print " ";


/* Putting the variables into local matrices. */

if dataset$=="";

   x=indvar;
   y=depvar;
   let depvar="y1" "y2";
   let indvar="intercep" "x1" "x2" "x3" "x4" "x5" "x6" "x7" "x8" "x9";
   if x[.,1]==ones(rows(x),1);
      indvar=indvar[1:cols(x)];
   else;
      indvar=indvar[1:cols(x)+1];
   endif;

else;

   {indvar,indexi}=indices(outfile,indvar);
   {depvar,indexd}=indices(outfile,depvar);
   nvec=setvars(outfile);
   open handle=^outfile for read;
   do until eof(handle);
      hulp1=readr(handle,100);
      i=1;
      do while i<=rows(indvar);
         hulp2=varget(indvar[i]);
         if rows(hulp2)>1;
            hulp2=hulp2|hulp1[.,indexi[i]];
         else;
            hulp2=hulp1[.,indexi[i]];
         endif;
         cnt=varput(hulp2,indvar[i]);
         i=i+1;
      endo;

      i=1;
      do while i<=rows(depvar);
         hulp2=varget(depvar[i]);
         if rows(hulp2)>1;
            hulp2=hulp2|hulp1[.,indexd[i]];
         else;
            hulp2=hulp1[.,indexd[i]];
         endif;
         cnt=varput(hulp2,depvar[i]);
         i=i+1;
      endo;

   endo;
   handle=close(handle);
  y=varget(depvar[1]);
   x=varget(indvar[1]);
   i=2;
   do while i<=rows(indvar);
      x=x~varget(indvar[i]);
      i=i+1;
   endo;

   if not(x[.,1]==ones(rows(x),1));
      indvar="intercep"|indvar;
   endif;

   nvec=setvars(outfile);

endif;


/* Checking for faulty data and deletion of rows with missing values. */

if (rows(x)/=rows(y));
   errorlog "x en y must have the same number of rows.";
   end;
elseif (cols(y)/=1);
   errorlog "y does not stand for a 0-1 variable.";
   end;
endif;

t=x~y;
_cases=rows(t);
g1=cols(x);
g2=cols(t);
t=packr(t);
_cases=_cases|rows(t);
y=t[.,g1+1:g2];
x=t[.,1:g1];

if sumc(sumc((y[.,1] .ne 1) .and (y[.,1] .ne 0)))/=0;
   errorlog "The dependent variables are not of the desired form.";
   end;
endif;


/* Transforming the dependent variables into contingency tables. */

hulp=zeros(rows(y),2);
hulp[.,1]=(y[.,1] .== 0);
hulp[.,2]=(y[.,1] .== 1);

y=hulp;
clear hulp;
g2=g2+1;



/* Multiple columns sort on the independent variable(s). */

t=x~y;
sortv=seqa(1,1,g1);
t=sortmc(t,sortv);


/* Summing over equal independent variable vectors. */

i=1;
k=1;
s=zeros(rows(t),g2);
do while i<=rows(t);
   cnt=t[i,1:g1];
   j=i;
   test=1;
   do while (j<=rows(t)) and test;
      j=j+1;
      if j<=rows(t);
         test=t[i,1:g1]==t[j,1:g1];
      else;
         test=0;
      endif;
   endo;
   j=j-1;
   if j<=rows(t);
      s[k,g1+1:g2]=sumc(t[i:j,g1+1:g2])';
      s[k,1:g1]=cnt;
      k=k+1;
   endif;
   i=j+1;
endo;

s=s[1:k-1,.];



/*
** Separating dependent and independent variables and adding a column
** of ones to the matrix of independent variables if necessary.
*/

y=s[.,cols(s)-1:cols(s)];
x=s[.,1:cols(s)-2];
if not(x[.,1]==1);
   x=ones(rows(x),1)~x;
endif;


/* Scanning of individual independent variables. */

output off;
{schat0,llt,grad,hess}=uniprob(0,x[.,1],y);
output on;

_file=2;
tel=2;
lvec=zeros(rows(indvar)-1,3);
do while tel<=rows(indvar);
   output off;
   {schat0,lvec[tel-1,1],grad,hess}=uniprob(zeros(2,1),x[.,1 tel],y);
   output on;
   tel=tel+1;
endo;


tel=1;
do while tel<=rows(lvec);
   lvec[tel,2]=lvec[tel,1]-llt;
   lvec[tel,3]=cdfchic(abs(2*lvec[tel,2]),1);
   tel=tel+1;
endo;


/* Ordering the independent variables. */

volgorde=seqa(2,1,rows(indvar)-1);
lvec=-lvec~volgorde;
lvec=sortc(lvec,1);
volgorde=lvec[.,4];
lvec=-lvec[.,1:3];



print "STEP  1";
print "=======";
print " ";
print " ";
mask=0~1~1~1;
let fmt[4,3]=
"*.*s "  14  4
"*.*lf " 14  4
"*.*lf " 14  4
"*.*lf " 14  4;
chosen="intercep"~llt~0~0;
print "Variables in model :";
print "--------------------";
format 10,4;
print $chosen[.,1]';
print " ";
print " ";
print "Candidates :";
print "------------";
print " ";
print "      variable        log-lik      chisquare        p-value";
print "      --------        -------      ---------        -------";
cnt=printfm(indvar[volgorde]~lvec[.,1]~(2*lvec[.,2])~lvec[.,3],mask,fmt);
print " ";
akkoord=(lvec[1,3]<=sig);
if akkoord;
   print $indvar[volgorde[1]] "will be kept in the model.";
else;
   print "No more variables will be taken in the model.";
endif;
print " ";
print "************************************************************"\
      "*****************";
print " ";
print " ";



/* Selection procedure. */

llt=0;
indvar=indvar[1]|indvar[volgorde];
x=x[.,1]~x[.,volgorde];
indv=indvar[1 2];
chosen=chosen|(indv[2]~lvec[1,.]);
if rows(indvar)>=3;
   pending=indvar[3:rows(indvar)];
   stoppen=not akkoord;
else;
   stoppen=1;
endif;
tel=1;
step=1;


do until stoppen;
   tel=1;
   step=step+1;
   w0=lvec[1,1];
   lvec=zeros(rows(pending),3);
   format 3,0;
   print "STEP" step;
   print "=======";
   print " ";
   print "Variables in model :";
   print "--------------------";
   format 10,4;
   print $chosen[.,1]';
   print " ";
   print "Candidates :";
   print "------------";
   do while tel<=rows(pending);
      output off;
      {schat0,lvec[tel,1],grad,hess}=
      uniprob(zeros(rows(chosen)+1,1),x[.,1:rows(chosen) rows(chosen)+tel],y);
      output on;
      lvec[tel,2]=lvec[tel,1]-w0;
      lvec[tel,3]=cdfchic(abs(2*lvec[tel,2]),1);
format 10,4;
      tel=tel+1;
   endo;

   volgorde=seqa(rows(chosen)+1,1,rows(pending));
   lvec=-lvec~volgorde;
   lvec=sortc(lvec,1);
   volgorde=lvec[.,4];
   lvec=-lvec[.,1:3];
   x=x[.,1:rows(chosen)]~x[.,volgorde];
   volgorde=volgorde-rows(chosen);
   pending=pending[volgorde];


   print " ";
   print "      variable        log-lik      chisquare        p-value";
   print "      --------        -------      ---------        -------";
   cnt=printfm(pending~lvec[.,1]~(2*lvec[.,2])~lvec[.,3],mask,fmt);
   print " ";
   if rows(pending)==1;

      print " ";
   endif;
   akkoord=(lvec[1,3]<=sig);
   if akkoord;
      print $pending[1] "will be kept in the model.";
   else;
      print "No more variables will be taken in the model.";
   endif;
   print " ";
   print "************************************************************"\
         "*****************";
   print " ";


   if akkoord;
      chosen=chosen|(pending[1]~lvec[1,.]);
      indv=indv|pending[1];
      if rows(pending)==1;
         stoppen=1;
      else;
         pending=pending[2:rows(pending)];
      endif;
   else;
      stoppen=1;
   endif;

endo;


output on;

print "Summary :";
print "---------";
print " ";
print "      variable        log-lik      chisquare        p-value";
print "      --------        -------      ---------        -------";
if rows(chosen)==2 and chosen[2,4]>sig;
   cnt=printfm(chosen[1,1 2],mask[1 2],fmt[1 2,.]);
else;
   cnt=printfm(chosen[1,1 2],mask[1 2],fmt[1 2,.]);
   print " ";
   cnt=printfm(chosen[2:rows(chosen),1 2]~(2*chosen[2:rows(chosen),3])~
   chosen[2:rows(chosen),4],mask,fmt);
endif;
print " ";
print " ";
print "************************************************************"\
      "*****************";


print " ";
tijd2=hsec;
tijd2=round((tijd2-tijd1)/100);
format /rd 5,0;
print "Computation time :" round(tijd2/60) "minutes" tijd2%60 "seconds.";
print " ";
print "************************************************************"\
      "*****************";
print " ";
print " ";

output off;
retp(chosen[.,1],chosen[.,2:4]);
endp;


