/* used for calculating counterfactuals */

#include utils.gau;

proc counterfact(h1,h2,startt,maxt,res);
local prue,probs1,probs2,pp1,pp2,ff,fmax,fhat,tk,ptrue,i;

probs1=calprobs1(h1,startt,maxt);
pp1=probs1[.,1]|probs1[.,2];
pp1=pp1|(1-sumc(pp1));

probs2=calprobs1(h2,startt,maxt);
pp2=probs2[.,1]|probs2[.,2];
pp2=pp2|(1-sumc(pp2));

ptrue=pp1|pp2;

if maxc(res[.,1].eq minc(res[.,1]));
fhat=maxc(res[.,1]);
"fhat=";; fhat;
else;
"function values not identical";
endif;
tk=30;
fmax=-1;
for i (1,rows(res),1);
"ok";
print i~res[i,.];
print /flush " ";
ff=cf(res[i,2],res[i,3],res[i,2],1,maxt,ptrue,fhat,tk);
"ff_no_cancer_progress=";;ff;
if (ff.gt fmax); fmax=ff; endif;
ff=cf(res[i,2],res[i,3],1,1,maxt,ptrue,fhat,tk);
"ff70=";;ff;
ff=cf(res[i,2],res[i,3],res[i,2],res[i,3],maxt,ptrue,fhat,tk);
"ff00=";;ff;
ff=cf(res[i,2],res[i,3],1,res[i,3],maxt,ptrue,fhat,tk);
"ff_no_cvd_progress=";;ff;

"in data 70";;sumc(sumc(probs1[1:tk-1,1:2]));
"in data 00";;sumc(sumc(probs2[1:tk-1,1:2]));
endfor ;

retp(ff);
endp;

proc cf(a,b,acf,bcf,maxt,ptrue,fhat,tk);
/* 
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,arows,acols,xm,xn,tp,xmp,xnp,xm1,xm2,xm3,xica,eps;

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 */


   aa=bb~(-aa);
    
    arows=rows(aa); acols=cols(aa);
    aa=aa~-eye(arows);
    aa=aa|(fhat~zeros(1,acols-1)~(-ones(1,arows)));

    xm=rows(aa); xn=cols(aa)-1;
    
ppp=(points[.,1]*acf)~(points[.,2]*bcf);   

tp=minc(ppp');
cc=zeros(1,rows(ppp));
for i (1,rows(ppp),1);
if tp[i].lt tk;
    cc[i]=1;
endif;
endfor ;
    

    cc=(0~cc~zeros(1,arows));
  
    aa=cc|aa;  
    aa=aa|zeros(1,cols(aa));
    
    xmp=rows(aa); xnp=cols(aa); xm1=0; xm2=0; xm3=rows(bb)+1;
    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];
retp(fct);
endp;

load res=res_70_00wm1; 
rrr=selif(res,res[.,1].eq maxc(res[.,1]));
load h1=haz_70wm;
load h2=haz_80wm;
load h3=haz_90wm;
load h4=haz_00wm;

startt=45; maxt=35;
ff=counterfact(h1,h4,startt,maxt,rrr);