/* used for calculating counterfactuals when we eliminate a risk*/

#include utils.gau;

proc(2)=counterfact(h1,h2,startt,maxt,res);
local prue,probs1,probs2,pp1,pp2,ff,fmax,fhat,tk,ptrue,i,fmin, smin,smax,fct,ff2;

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;


smin=zeros(maxt,8); smax=zeros(maxt,8);

for tk (1,maxt,1); 

"TK = ";; tk;

fmax=-10+zeros(8,1);
fmin=10+zeros(8,1);
for i (1,rows(res),1);

print /flush  i;;
print res[i,.];;

fct=orig_lp(res[i,2],res[i,3],maxt,ptrue,tk,1,1);
fct;

fhat=fct;

/* "MAXIMUM"; */


ff=cf(res[i,2],res[i,3],-1,res[i,3],maxt,ptrue,fhat,tk,1);
/* "ff_no_cvd_progress=";;ff; */
if (ff.gt fmax[5]); fmax[5]=ff; endif;

ff=cf(res[i,2],res[i,3],res[i,2],-1,maxt,ptrue,fhat,tk,1);
/* "ff_no_cvd_progress=";;ff; */
if (ff.gt fmax[6]); fmax[6]=ff; endif;

ff=cf(res[i,2],res[i,3],-1,1,maxt,ptrue,fhat,tk,1);
/* "ff_no_cvd_progress=";;ff; */
if (ff.gt fmax[7]); fmax[7]=ff; endif;

ff=cf(res[i,2],res[i,3],1,-1,maxt,ptrue,fhat,tk,1);
/* "ff_no_cvd_progress=";;ff; */
if (ff.gt fmax[8]); fmax[8]=ff; endif;

ff=cf(res[i,2],res[i,3],1,1,maxt,ptrue,fhat,tk,1);
/* "ff_no_cancer_progress=";;ff; */
if (ff.gt fmax[1]); fmax[1]=ff; endif;


ff=cf(res[i,2],res[i,3],res[i,2],1,maxt,ptrue,fhat,tk,1);
/* "ff_no_cancer_progress=";;ff; */
if (ff.gt fmax[2]); fmax[2]=ff; endif;

ff=cf(res[i,2],res[i,3],1,res[i,3],maxt,ptrue,fhat,tk,1);
/* "ff_no_cvd_progress=";;ff; */
if (ff.gt fmax[3]); fmax[3]=ff; endif;

ff=cf(res[i,2],res[i,3],res[i,2],res[i,3],maxt,ptrue,fhat,tk,1);
/* "ff_no_cvd_progress=";;ff; */
if (ff.gt fmax[4]); fmax[4]=ff; endif;

/* min */
/* "MINIMUM"; */

ff=cf(res[i,2],res[i,3],1,1,maxt,ptrue,fhat,tk,2);
/* "ff_no_cancer_progress=";;ff; */
if (ff.lt fmin[1]); fmin[1]=ff; endif;

ff=cf(res[i,2],res[i,3],res[i,2],1,maxt,ptrue,fhat,tk,2);
/* "ff_no_cancer_progress=";;ff; */
if (ff.lt fmin[2]); fmin[2]=ff; endif;

ff=cf(res[i,2],res[i,3],1,res[i,3],maxt,ptrue,fhat,tk,2);
/* "ff_no_cvd_progress=";;ff; */
if (ff.lt fmin[3]); fmin[3]=ff; endif;

ff=cf(res[i,2],res[i,3],res[i,2],res[i,3],maxt,ptrue,fhat,tk,2);
/* "ff_no_cvd_progress=";;ff; */
if (ff.lt fmin[4]); fmin[4]=ff; endif;

ff=cf(res[i,2],res[i,3],-1,res[i,3],maxt,ptrue,fhat,tk,2);
/* "ff_no_cvd_progress=";;ff; */
if (ff.lt fmin[5]); fmin[5]=ff; endif;

ff=cf(res[i,2],res[i,3],res[i,2],-1,maxt,ptrue,fhat,tk,2);
/* "ff_no_cvd_progress=";;ff; */
if (ff.lt fmin[6]); fmin[6]=ff; endif;

ff=cf(res[i,2],res[i,3],-1,1,maxt,ptrue,fhat,tk,2);
/* "ff_no_cvd_progress=";;ff; */
if (ff.lt fmin[7]); fmin[7]=ff; endif;

ff=cf(res[i,2],res[i,3],1,-1,maxt,ptrue,fhat,tk,2);
/* "ff_no_cvd_progress=";;ff; */
if (ff.lt fmin[8]); fmin[8]=ff; endif;

endfor;

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

"fitted in 70";; fmin[1]~fmax[1];
"fitted with no_cancer_progress=";; fmin[2]~fmax[2];
"fitted with no_cvd_progress=";; fmin[3]~fmax[3];
"fitted in 00";; fmin[4]~fmax[4];
*/
format/rd 10,6;

"Min:";;fmin'; "MAX:";; fmax';
print /flush " ";

smin[tk,.]=1-fmax';
smax[tk,.]=1-fmin';

endfor;

retp(smin,smax);
endp;


proc cf(a,b,acf,bcf,maxt,ptrue,fhat,tk,max_or_min);
/* 
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. */

/* max_or_min = 1 or 2 depending on max or min, respectively */

local ppp,cc,fct,points,coef_x0,coef_x1,aa,bb,arows,acols,xm,xn,tp,xmp,xnp,xm1,xm2,xm3,xica,eps;

if .not (((a.eq acf ).or (acf.eq 1).or (acf.eq -1)).and ((bcf.eq 1).or (bcf.eq b).or (bcf.eq -1))); "something wrong with counerfactuals"; endif;

points= findpoints(a,b,maxt,tk,acf,bcf);

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]*abs(acf)+(tk+1)*(acf.le 0))~(points[.,2]*abs(bcf)+(tk+1)*(bcf.le 0));   

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

    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-11;
    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 -0.1);
        print "icase=";;
        print xica;;
        print " in lp_max_val";
        stop;
    endif; 
    fct=aa[1,1];
    if max_or_min.eq 2; fct=-fct; endif;
retp(fct);
endp;

proc orig_lp(a,b,maxt,ptrue,tk,acf,bcf);
/* This routine solves the original LP problem */



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,tk,acf,bcf);


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);

    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-11;
    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 -0.1);
        print "icase=";;
        print xica;;
        print " in lp_max_val";
        stop;
    endif; 
    fct=aa[1,1];
retp(fct);
endp;

proc(0)=wrilab(i);

if i.eq 1; "No Progress (fitted in 70)";
elseif i.eq 2; "Progress in CVD";
elseif i.eq 3; "Progress in Cancer";
elseif i.eq 4; "Progress in Both (fitted in end year)";
elseif i.eq 5; "Elim. CVD --- Cancer as in end year";
elseif i.eq 6; "Elim. Cancer --- CVD as in end year";
elseif i.eq 7; "Elim. CVD --- Cancer as in 70";
elseif i.eq 8; "Elim. Cancer --- CVD as in 70";
else; stop;
endif;
endp;

proc(0)=writeres(s1l,s2l,s3l,s1h,s2h,s3h,a1,a2,a3,labl);
local ii,i;

ii=30;
format /rd 9,3;

"Change in actual prob."; "&";;
"$";; a1[ii,1];; "-";; a1[ii,2];; "$&";;
"$";; a2[ii,1];; "-";; a2[ii,2];; "$&";;
"$";; a3[ii,1];; "-";; a3[ii,2];; "$";;
"\\\\";

for i (1,cols(s1l),1);
wrilab(i);; "&";;
"$(";; s1l[ii,i];; ",";; s1h[ii,i];; ")$&";;
"$(";; s2l[ii,i];; ",";; s2h[ii,i];; ")$&";;
"$(";; s3l[ii,i];; ",";; s3h[ii,i];; ")$";;
"\\\\";
endfor;
endp;



" COUNTERFACTUAL";



startt=45; maxt=35;

load h1=haz_70wm;
load h2=haz_80wm;
load h3=haz_90wm;
load h4=haz_00wm;

load res=res_70_00wm1; rrr=selif(res,res[.,1].eq maxc(res[.,1]));
{smin_wm_70_00,smax_wm_70_00}=counterfact(h1,h4,startt,maxt,rrr);
save smin_wm_70_00; save smax_wm_70_00;

load res=res_70_80wm1; rrr=selif(res,res[.,1].eq maxc(res[.,1]));
{smin_wm_70_80,smax_wm_70_80}=counterfact(h1,h2,startt,maxt,rrr);
save smin_wm_70_80; save smax_wm_70_80;

load res=res_70_90wm1; rrr=selif(res,res[.,1].eq maxc(res[.,1]));
{smin_wm_70_90,smax_wm_70_90}=counterfact(h1,h3,startt,maxt,rrr);
save smin_wm_70_90; save smax_wm_70_90;

load h1=haz_70wf;
load h2=haz_80wf;
load h3=haz_90wf;
load h4=haz_00wf;

load res=res_70_80wf1; rrr=selif(res,res[.,1].eq maxc(res[.,1]));
{smin_wf_70_80,smax_wf_70_80}=counterfact(h1,h2,startt,maxt,rrr);
save smin_wf_70_80; save smax_wf_70_80;

load res=res_70_90wf1; rrr=selif(res,res[.,1].eq maxc(res[.,1]));
{smin_wf_70_90,smax_wf_70_90}=counterfact(h1,h3,startt,maxt,rrr);
save smin_wf_70_90; save smax_wf_70_90;

load res=res_70_00wf1; rrr=selif(res,res[.,1].eq maxc(res[.,1]));
{smin_wf_70_00,smax_wf_70_00}=counterfact(h1,h4,startt,maxt,rrr);
save smin_wf_70_00; save smax_wf_70_00;

load h1=haz_70bm;
load h2=haz_80bm;
load h3=haz_90bm;
load h4=haz_00bm;

load res=res_70_80bm1; rrr=selif(res,res[.,1].eq maxc(res[.,1]));
{smin_bm_70_80,smax_bm_70_80}=counterfact(h1,h2,startt,maxt,rrr);
save smin_bm_70_80; save smax_bm_70_80;

load res=res_70_90bm1; rrr=selif(res,res[.,1].eq maxc(res[.,1]));
{smin_bm_70_90,smax_bm_70_90}=counterfact(h1,h3,startt,maxt,rrr);
save smin_bm_70_90; save smax_bm_70_90;

load res=res_70_00bm1; rrr=selif(res,res[.,1].eq maxc(res[.,1]));
{smin_bm_70_00,smax_bm_70_00}=counterfact(h1,h4,startt,maxt,rrr);
save smin_bm_70_00; save smax_bm_70_00;

load h1=haz_70bf;
load h2=haz_80bf;
load h3=haz_90bf;
load h4=haz_00bf;

load res=res_70_80bf1; rrr=selif(res,res[.,1].eq maxc(res[.,1]));
{smin_bf_70_80,smax_bf_70_80}=counterfact(h1,h2,startt,maxt,rrr);
save smin_bf_70_80; save smax_bf_70_80;

load res=res_70_90bf1; rrr=selif(res,res[.,1].eq maxc(res[.,1]));
{smin_bf_70_90,smax_bf_70_90}=counterfact(h1,h3,startt,maxt,rrr);
save smin_bf_70_90; save smax_bf_70_90;

load res=res_70_00bf1; rrr=selif(res,res[.,1].eq maxc(res[.,1]));
{smin_bf_70_00,smax_bf_70_00}=counterfact(h1,h4,startt,maxt,rrr);
save smin_bf_70_00; save smax_bf_70_00;



output file=Table_4_augmented.tex, reset;

cls;


labl="Adriana";

"% Created on ";; datestr(date);
"% by counter_tab.gau";
" "; " ";" "; " ";

"\\begin{center}";
"\\begin{tabular}{lccc}";
"\\multicolumn{4}{c}{\\textbf{TABLE: ";
"Counterfactual Probability of Surviving Age 75}} \\\\";
"\\\\";
"\\\\";
"\\multicolumn{4}{c}{\\textbf{Results for White Males}} \\\\";
"\\\\";
" & 1970--1980 & 1970--1990 & 1970--2000 \\\\";
"\\\\";

load smi00=smin_wm_70_00;
load sma00=smax_wm_70_00;
load smi90=smin_wm_70_90;
load sma90=smax_wm_70_90;
load smi80=smin_wm_70_80;
load sma80=smax_wm_70_80;

load act00=actual_wm_70_00;
load act80=actual_wm_70_80;
load act90=actual_wm_70_90;

writeres(smi80,smi90,smi00,sma80,sma90,sma00,act80,act90,act00,labl);

"\\\\";


"\\multicolumn{4}{c}{\\textbf{Results for White Females}} \\\\";
"\\\\";
" & 1970--1980 & 1970--1990 & 1970--2000 \\\\";
"\\\\";


load smi00=smin_wf_70_00;
load sma00=smax_wf_70_00;
load smi90=smin_wf_70_90;
load sma90=smax_wf_70_90;
load smi80=smin_wf_70_80;
load sma80=smax_wf_70_80;
load act00=actual_wf_70_00;
load act80=actual_wf_70_80;
load act90=actual_wf_70_90;

writeres(smi80,smi90,smi00,sma80,sma90,sma00,act80,act90,act00,labl);



"\\\\";
"\\\\";


"\\end{tabular}";
"\\end{center}";

"\\newpage";

"\\begin{center}";
"\\begin{tabular}{lccc}";
"\\multicolumn{4}{c}{\\textbf{TABLE: ";
"Counterfactual Probability of Surviving Past 75 (cont.)}} \\\\";
"\\\\";
"\\\\";


"\\multicolumn{4}{c}{\\textbf{Results for Black Males}} \\\\";
"\\\\";
" & 1970--1980 & 1970--1990 & 1970--2000 \\\\";
"\\\\";




load smi00=smin_bm_70_00;
load sma00=smax_bm_70_00;
load smi90=smin_bm_70_90;
load sma90=smax_bm_70_90;
load smi80=smin_bm_70_80;
load sma80=smax_bm_70_80;
load act00=actual_bm_70_00;
load act80=actual_bm_70_80;
load act90=actual_bm_70_90;

writeres(smi80,smi90,smi00,sma80,sma90,sma00,act80,act90,act00,labl);



"\\\\";


"\\multicolumn{4}{c}{\\textbf{Results for Black Females}} \\\\";
"\\\\";
" & 1970--1980 & 1970--1990 & 1970--2000 \\\\";
"\\\\";


load smi00=smin_bf_70_00;
load sma00=smax_bf_70_00;
load smi90=smin_bf_70_90;
load sma90=smax_bf_70_90;
load smi80=smin_bf_70_80;
load sma80=smax_bf_70_80;
load act00=actual_bf_70_00;
load act80=actual_bf_70_80;
load act90=actual_bf_70_90;

writeres(smi80,smi90,smi00,sma80,sma90,sma00,act80,act90,act00,labl);



"\\\\";
"\\\\";


"\\end{tabular}";
"\\end{center}";

" "; " ";
/* ******************************************** */
" "; " ";


output off;