#ifos2win
    library lpmt;
    #include lpmt.sdf;
#else
#endif

proc(0)=showdim(a);
print /flush "the matrix has dimension ";;
print /flush rows(a)~cols(a);
endp;

proc findpoi1(gx,gy,s1,s2);
/* For two lines with slope s1 and s2, and sets of points on 
the x and y axes (gx and gy), consider all the resultiong polygoms 
in the first quadrant. Now vertigaly aggregate polygons above the 
steeper line. This program finds an interior point in all of the 
remaining polygons in such a way that the x--coordinated are the same
if two polygons are on top of each other. Each row in the output 
corresponds to a point */
local smi, sma, ppp, gxx, gyy, i;

smi=minc(s1|s2); 
sma=maxc(s1|s2); 

gxx=gx|(gy/sma)|(gy/smi);
gxx=0|gxx|(maxc(gxx)+1);
gxx=unique(gxx,1);
gxx=sortc(gxx,1);
gxx=(gxx[1:(rows(gxx)-1)]+gxx[2:rows(gxx)])/2;

for i (1,rows(gxx),1);
    gyy=0|gy|(gxx[i]*sma)|gxx[i]*smi;
    gyy=unique(gyy,1);
    gyy=selif(gyy,gyy.le (sma*gxx[i]));
    gyy=gyy|(maxc(gyy)+1);
    gyy=sortc(gyy,1);
    gyy=(gyy[1:(rows(gyy)-1)]+gyy[2:rows(gyy)])/2;
    if (i.eq 1);
        ppp=(gxx[i]*ones(rows(gyy),1))~gyy;
    else;
        ppp=ppp|((gxx[i]*ones(rows(gyy),1))~gyy);
    endif;
endfor;
retp(ppp);
endp;

proc findpoi(xp,yp,maxt,a1,a2);
local points,ii,jj,nxp,nyp,xx,yy,xx1,xx2,yy1,yy2,xl,xr,yl,yr,sl;

if(a1.lt a2); stop; "a1.lt a2"; endif;
nxp=rows(xp); nyp=rows(yp);
maxt=maxt+1;
points=maxt|maxt;

for ii (1,  nxp, 1);
if (ii.eq 1); xx=xp[1]/2; else; xx=(xp[ii]+xp[ii-1])/2; endif;
points=points~(xx|maxt);
endfor;

for ii (1,  nyp, 1);
if (ii.eq 1); yy=yp[1]/2; else; yy=(yp[ii]+yp[ii-1])/2; endif;
points=points~(maxt|yy);
endfor;


for ii (1,  nxp, 1);
for jj (1,  nyp, 1);

xr=xp[ii];
if (ii.eq 1); xl=0; else; xl=xp[ii-1]; endif;
yr=yp[jj];
if (jj.eq 1); yl=0; else; yl=yp[jj-1]; endif;

if a1*xr.le  yl; goto aaa1; endif;
if a2*xl.ge yr; goto aaa1; endif;



sl=(yr-yl)/(xr-xl);

xx1=(yr+sl*xl)/(a1+sl);
yy1=a1*xx1;
xx2=(yr+sl*xl)/(a2+sl);
yy2=a2*xx2;



if .not ( (xx1.ge xl) .and (xx1.le xr) .and (yy1.ge yl) .and (yy1.le yr)); 
    xx1=xl; yy1=yr; 
endif;
if .not ( (xx2.ge xl) .and (xx2.le xr) .and (yy2.ge yl) .and (yy2.le yr)); 
    xx2=xr; yy2=yl; 
endif;

xx1=(xx1+xx2)/2; yy1=(yy1+yy2)/2;
if (xx1*a2.ge yy1).or  (xx1*a1.le yy1); goto aaa1; endif;
points=points~(xx1|yy1);



aaa1:
endfor;
endfor;
retp(points);
endp;

proc findpoints(a,b,maxt);
/* finds the points of support */
local points,ss,xp,yp,maxs;

ss=seqa(1,1,maxt);
xp=ss|(ss/a);
yp=ss|(ss/b);
maxs=maxc(1|(1/a)|(1/b))*maxt;
xp=unique(xp,1);
yp=unique(yp,1);

if (a.le b);
points=findpoi(xp,yp,maxs,1,a/b);
else;
points=findpoi(xp,yp,maxs,a/b,1);
endif;

points=points';

retp(points);
endp;

proc findpoints1(a,maxt);
/* finds the points of support. Each row is a point */
local points,ss,xp,yp,maxs;

ss=seqa(1,1,maxt);
xp=ss|(ss/a);
yp=ss;

points=findpoi1(xp,yp,a,1);
retp(points);
endp;

/* ************************************************************** */



proc calprobs(h,startt,maxt);
/* turns a hazard rate into a distribution */
local probs, alive, i;

probs=zeros(maxt,2);
alive=1;
for i (1,maxt,1);
probs[i,1]=alive*h[i+startt,1];
probs[i,2]=alive*h[i+startt,2];
alive=alive*(1-h[i+startt,1]-h[i+startt,2]);
endfor;

retp(probs);
endp;

proc calprobs1(h,startt,maxt);
/* turns a hazard rate into a distribution */
/* same as calprobs */
local probs, alive, i,hh;

probs=zeros(maxt,2);
hh=h[.,1]+h[.,2];
hh=1-hh[startt+1:startt+maxt];
alive=cumprodc(hh);
alive=1|alive[1:maxt-1];

probs[.,1]=alive.*h[1+startt:maxt+startt,1];
probs[.,2]=alive.*h[1+startt:maxt+startt,2];

retp(probs);
endp;

proc findcoefs(ppp,maxt);
local tp, d1p, d2p, c, i, coef;

coef=zeros( rows(ppp), 2*maxt+1);
tp=minc(ppp');
d1p=ppp[.,1].lt ppp[.,2] + 0.5*(ppp[.,1].eq ppp[.,2]);
d2p=ppp[.,1].gt ppp[.,2] + 0.5*(ppp[.,1].eq ppp[.,2]);
c=tp.gt maxt;
d1p=d1p.*(1-c);
d2p=d2p.*(1-c);
tp=int(tp)+1;
for i (1,rows(ppp),1);
if (c[i].eq 0); 
   coef[i,tp[i]]=d1p[i]; 
   coef[i,maxt+tp[i]]=d2p[i]; 
endif;
  coef[i,2*maxt+1]=c[i];
endfor ;
retp(coef);
endp;

proc marg_eq(ppp);
/* given a set of points of support, ppp, this routine calculates the
matrix that corresponds to the restriction that the marginal distribution of the first 
elemement is the same for two distributions with those poinst of support */

local xp, i, x, constr, c;


xp=unique(ppp[.,1],1);


for i (1,rows(xp),1);
c=((ppp[.,1]'.eq xp[i])~(-(ppp[.,1]'.eq xp[i])));
if (i.eq 1);
    constr=c;
else;
    constr=constr|c;
endif;
endfor;
retp(constr);
endp;
proc lp_max_val(aa,bb,nmet);

/*  This routibe returns the function value for the LP problem that 
    determines whether aa*x=bb has a solution with x(i) >= 0 for all i */

local fct,arows,acols,xm,xn,xmp,xnp,xm1,xm2,xm3,xica,eps,cc,startv,ii;


if (nmet.eq 1);

    arows=rows(aa); acols=cols(aa);
    aa=eye(arows)~aa;
    /* here we add the v's */
    
    cc=-ones(arows,1)|zeros(acols,1);
    startv=bb|zeros(acols,1)|zeros(arows,1);
    #ifos2win

           
        struct LP lp0;
        lp0=createLP;
        lp0.a=aa;
        lp0.b=bb;
        lp0.c=cc;
        lp0.l=0;
        lp0.u=1e200;
        
        struct LPcontrol c0;
        c0=createLPcontrol;
        c0.output=0;
        c0.constraintType=3;
        c0.minimize=1;
        c0.start=startv;
        
        struct LPout out0;
        out0=lpmt(lp0,c0);
        fct=out0.optimumvalue;
        print "return code   ";; out0.returncode; 
        print "function value";; fct; 
   
        
    #else
       print /flush "OS not recognized"; 
       stop;
    #endif
elseif (nmet.eq 2);
    aa=bb~(-aa);
    
    arows=rows(aa); acols=cols(aa);
    aa=aa~-eye(arows);
    xm=rows(aa); xn=cols(aa)-1;
    
    aa=(0~zeros(1,acols-1)~(-ones(1,arows)))|aa;
    aa=aa|zeros(1,cols(aa));
    
    xmp=rows(aa); xnp=cols(aa); xm1=0; xm2=0; xm3=rows(bb);
    xica=0; 
    eps=1.0e-8;
    aa=aa';
    
    #IFUNIX
       dlibrary /home/honore/adriana/simplex1.so;
       dllcall simplex_(aa,xm,xn,xmp,xnp,xm1,xm2,xm3,xica);  
    #else
    #ifos2win
       dlibrary c:\research\forgau95\sim_new\release\sim_new.dll;
       dllcall isimpl(aa,xm,xn,xmp,xnp,xm1,xm2,xm3,xica,eps); 
    #else
       print /flush "OS not recognized"; 
       stop;
    #endif
    #endif
    if (xica.le -100);
        print "icase=";;
        print xica;;
        print " in lp_max_val";
        stop;
    endif; 
    fct=aa[1,1];
else;
    print /flush "nmet not valid in lp_max_val";
endif;
retp(fct);
endp;





proc func_2mult(a,b,maxt,ptrue,nmet);
/* This routine calculates the function value for the case 
where both T1 and T2 are modeled with multiplicative effects.
if nmet=1, then the routibe uses GAUSS's LP program. if nmet=2, 
then the routibe uses the LP program from Numerical Recipes. */
local ppp,cc,fct,points,coef_x0,coef_x1,aa,bb;

points= findpoints(a,b,maxt);

ppp=points;
coef_x0=findcoefs(ppp,maxt);
ppp=(points[.,1]*a)~(points[.,2]*b);
coef_x1=findcoefs(ppp,maxt);

aa=coef_x0~coef_x1; aa=aa~ones(rows(aa),1);

aa=aa';
bb=ptrue|1;
/* these are the real constraints */

fct=lp_max_val(aa,bb,nmet);
retp(fct);
endp;
 
proc func_1mult(a,maxt,ptrue,nmet);
/* This routine calculates the function value for the case 
where only T1 is modeled with multiplicative effects.
if nmet=1, then the routibe uses GAUSS's LP program. if nmet=2, 
then the routibe uses the LP program from Numerical Recipes. */
local ppp,cc,fct,points,coef_x0,coef_x1,aa,bb,mar;


points=findpoints1(a,maxt);
ppp=points;
coef_x0=findcoefs(ppp,maxt)';
ppp=(points[.,1]*a)~points[.,2];  
coef_x1=findcoefs(ppp,maxt)';
mar=marg_eq(ppp);
mar=mar[1:rows(mar)-1,.];
aa=(coef_x0~zeros(rows(coef_x0),cols(coef_x1)))|
   (zeros(rows(coef_x0),cols(coef_x1))~coef_x1); 
aa=aa|mar;
bb=ptrue|zeros(rows(mar),1);

/* these are the real constraints */

fct=lp_max_val(aa,bb,nmet);
retp(fct);
endp;