/* This version allows for kinks in the linear relationship. 
It is a modification of the old util 1, and we have added the new needed routines. */

#include optmum.prc;
clear bo_dat, bo_ptrue;

proc grids(h1,h2,startt,maxt,breakt,thmin1,thmax1,dth1);
local res, fct, a1, a2 ,b1,b2, adr,ngr1,ngr2,ptrue,probs,pp1,pp2, theta1,theta2,theta3, theta4,ia1,ia2,ib1,ib2,t1,t2;

format /rd 9,5;
t1=breakt;
t2=breakt;

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

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

ptrue=pp1|pp2;

ngr1=(thmax1-thmin1)./dth1;
ngr1=floor(ngr1);
ngr1=ngr1.*(ngr1.gt 0)+ (ngr1.le 0);
ngr1=ngr1+1;
"NGR     ";; ngr1';
"THMIN   ";; thmin1';
"THMAX   ";; thmax1';
"DTH     ";; dth1';
theta1=seqa(thmin1[1],dth1[1],ngr1[1]);
theta2=seqa(thmin1[2],dth1[2],ngr1[2]);
theta3=seqa(thmin1[3],dth1[3],ngr1[3]);
theta4=seqa(thmin1[4],dth1[4],ngr1[4]);

res=1000+zeros(1,6);

for ia1 (1,rows(theta1),1);
a1=theta1[ia1];
for ia2 (1,rows(theta2),1);
a2=theta2[ia2];
for ib1 (1,rows(theta3),1);
b1=theta3[ib1];
for ib2 (1,rows(theta4),1);
b2=theta4[ib2];

adr=hsec;
fct=func_2mult_wk(a1,a2,t1,b1,b2,t2,maxt,ptrue,2);
res=res|(fct~a1~a2~b1~b2~(hsec-adr));
print /flush maxc(res[2:rows(res),1])~fct~a1~a2~b1~b2~(hsec-adr); 


endfor;
endfor;
endfor;
endfor;

retp(res[2:rows(res),.]);
endp;


proc grid_seq(h1,h2,startt,maxt,breakt,dth1);
local res,fr,rrr,b,thmin1,thmax1,gridno;

format /rd 6,2;


b=getstart(h1,h2,startt,maxt,breakt);
thmin1=b-dth1*.9;
thmax1=b+dth1*.9;

output on;
" ";
"first grid";
" ";
thmin1~thmax1~dth1;
" ";
output off;
res=grids(h1,h2,startt,maxt,breakt,thmin1,thmax1,dth1);
fr=res[.,1];
fr=unique(fr,1);
fr=sortc(fr,1);
rrr=selif(res,res[.,1].ge fr[rows(fr)]);
output on;
"first round minimizers";
" ";
selif(res,res[.,1].ge maxc(res[.,1]));

gridno=1;

do until meanc(dth1).lt 0.005;
gridno=gridno+1;
dth1=dth1*0.66667;
thmin1=minc(rrr[.,2:5])-dth1; thmax1=maxc(rrr[.,2:5])+dth1; 

 " ";
"grid number ";;gridno;
" ";
thmin1~thmax1~dth1;
" ";
output off;

res=res|grids(h1,h2,startt,maxt,breakt,thmin1,thmax1,dth1);
fr=res[.,1];
fr=unique(fr,1);
fr=sortc(fr,1);
rrr=selif(res,res[.,1].ge fr[rows(fr)]);
output on;
"minimizers fr round ";;gridno;
" ";
selif(res,res[.,1].ge maxc(res[.,1]));
endo;

output off;
retp(res);
endp;


PROC kinked_line(s,a,b,t1);
RETP((s.lt t1).*s*a + (s.ge t1).*(t1*a+(s-t1).*b));
ENDP;

PROC kinked_line_inv(s,a,b,t1);
RETP((s.lt a*t1).*s/a + (s.ge a*t1).*(((s-a*t1)/b)+t1));
ENDP;

PROC fi_poi(s,a1,b1,t1,a2,b2,t2);

/* finds the basic points of support */

local points,np1,np2,i1,i2,xp1,xp2,ixp1,ixp2,ii1,ii2,c1,c2,c3,pa1,pa2,pa3,iter,p1,p2,p1k,p2k,tpa1,tpa2,dtpa,tau,tau1,tau2
    ,aa1,aa2,bb1,bb2,s1,s2,s2a,s2b,icross, itemp;

s=s|t1|t2;
p1=s|kinked_line_inv(s,a1,b1,t1);
p1=p1|(maxc(p1)+1);

p2=s|kinked_line_inv(s,a2,b2,t2);
p2=p2|(maxc(p2)+1);


p1=0|p1;
p2=0|p2;

p1=sortc(p1,1);
p2=sortc(p2,1);


p1=unique(p1,1);
p2=unique(p2,1);

p1=sortc(p1,1);
p2=sortc(p2,1);


p1k=kinked_line(p1,a1,b1,t1);
p2k=kinked_line(p2,a2,b2,t2);

np1=rows(p1)-1; np2=rows(p2)-1;

xp1=(p1[1:np1]+p1[2:np1+1])/2;
xp2=(p2[1:np2]+p2[2:np2+1])/2;

ixp1=int(xp1~kinked_line(xp1,a1,b1,t1));
/* cause 1 failure times in each period*/

ixp2=int(xp2~kinked_line(xp2,a2,b2,t2));
/* cause 2 failure times in each period*/

points=-100~-100;
for i1 (1,np1,1);
for i2 (1,np2,1);

ii1=(ixp1[i1,1].eq ixp2[i2,1]);
ii2=(ixp1[i1,2].eq ixp2[i2,2]);

if ii1*ii2; 

/* tie in both */
/* There are three situations dependning on whether 2,3 or 4 events can happen. We comhine the first two below*/


pa1=p1[i1]~p2[i2+1];
pa2=p1[i1+1]~p2[i2];

tpa1=p1k[i1]~p2k[i2+1];
tpa2=p1k[i1+1]~p2k[i2];


/*  
    (s1,s2) describes a point in the rectangle 0<s1<1, 0<s2<1
    the duration t1 when X=0 is given by pa[1]+s1*(pb[1]-pa[1])
        the duration t2 when X=0 is given by pa[2]+s2*(pb[2]-pa[2])
    t1=t2 when X=0 if  pa[1]+s1*(pb[1]-pa[1]) = pa[2]+s2*(pb[2]-pa[2]) 
                                 or s2 =  (pa[1]-pa[2])/ (pb[2]-pa[2]) + s1*(pb[1]-pa[1])/ (pb[2]-pa[2]) 
                       = aa1 +bb1*s1
    
    the duration t1 when X=1 is given by tpa[1]+s1*(tpb[1]-tpa[1])
        the duration t2 when X=1 is given by tpa[2]+s2*(tpb[2]-tpa[2])
    t1=t2 when X=1 if  tpa[1]+s1*(tpb[1]-tpa[1]) = tpa[2]+s2*(tpb[2]-tpa[2]) 
                                 or s2 =  (tpa[1]-tpa[2])/ (tpb[2]-tpa[2]) + s1*(tpb[1]-tpa[1])/ (tpb[2]-tpa[2])
                            = aa2 +bb2*s1 

    the two lines that give t1=t2 cross at 
    s1=  (aa1-aa2)/(bb2-bb1) (let s2 be the corresponding value. It should also be between 0 and 1 if s1 is)
    if s1 is between 0 and 1, then there are for cases. Otherwise, there are three.

    if s1 is between 0 and 1, the four values can be chosen as 
    (s1,s2/2), (s1,(s2+1)/2), (s1/2, ((aa1 +bb1*s1/2+aa2 +bb2*s1/2 ))/2),  ((1+s1)/2, ((aa1 +bb1*(1+s1)/2/2+aa2 +bb2*(1+s1)/2 ))/2)

    if s1 is not between 0 and 1 then the three points can be found by setting s1=1/2, 
    and then calculate the two correcponding values of s2. one then combines the values of s1 
    with the averages of the s2's, the average of 0 and the min s2, and teh average of 1 and the max s2.
    
*/


aa1=(pa1[1]-pa1[2])/ (pa2[2]-pa1[2]) ;
bb1=(pa2[1]-pa1[1])/ (pa2[2]-pa1[2]) ;
aa2=(tpa1[1]-tpa1[2])/ (tpa2[2]-tpa1[2]);
bb2=(tpa2[1]-tpa1[1])/ (tpa2[2]-tpa1[2]);

icross=1;
if bb2.eq bb1;
     icross =0; 
    s1=-9.99999e9;
else;
    s1=  (aa1-aa2)/(bb2-bb1) ;
    s2= aa1+bb1*s1;
endif;

if (s1.le 0).or (s1.ge 1); icross=0; endif;

if icross;


tau1=s1; tau2=s2/2;
points=points | ( (tau1*pa2[1]+(1-tau1)*pa1[1]) ~ (tau2*pa2[2]+(1-tau2)*pa1[2]) ); 

tau1=s1; tau2=(s2+1)/2;
points=points | ( (tau1*pa2[1]+(1-tau1)*pa1[1]) ~ (tau2*pa2[2]+(1-tau2)*pa1[2]) ); 

tau1=s1/2;
tau2= (aa1 +bb1*tau1+aa2 +bb2*tau1)/2;
points=points | ( (tau1*pa2[1]+(1-tau1)*pa1[1]) ~ (tau2*pa2[2]+(1-tau2)*pa1[2]) ); 

tau1=(1+s1)/2;
tau2= (aa1 +bb1*tau1+aa2 +bb2*tau1)/2;
points=points | ( (tau1*pa2[1]+(1-tau1)*pa1[1]) ~ (tau2*pa2[2]+(1-tau2)*pa1[2]) ); 

else;
/* add 2 or three points */
tau1=0.5;
s2a=aa1+bb1*tau1;
s2b=aa2+bb2*tau1;


if s2a.ne s2b;
tau2=(s2a+s2b)/2;

points=points | ( (tau1*pa2[1]+(1-tau1)*pa1[1]) ~ (tau2*pa2[2]+(1-tau2)*pa1[2]) ); 
endif;

tau2=(1+maxc(s2a|s2b))/2;

points=points | ( (tau1*pa2[1]+(1-tau1)*pa1[1]) ~ (tau2*pa2[2]+(1-tau2)*pa1[2]) ); 

tau2=(minc(s2a|s2b))/2;

points=points | ( (tau1*pa2[1]+(1-tau1)*pa1[1]) ~ (tau2*pa2[2]+(1-tau2)*pa1[2]) ); 

endif;
endif;

if ii1*(1-ii2); 
/* tie in the first period */

pa1=p1[i1]~p2[i2+1];
pa2=p1[i1+1]~p2[i2];


tpa1=pa1; tpa2=pa2;
dtpa=tpa1-tpa2;
tau=dtpa[1]/(dtpa[1]-dtpa[2]);

if tau.lt 0; "tau < 0"; stop; endif;
if tau.gt 1; "tau > 1"; stop; endif;
tau1=tau/2;
tau2=(tau+1)/2;
points=points|(tau1*pa1+(1-tau1)*pa2)|(tau2*pa1+(1-tau2)*pa2); 
endif;

if (1-ii1)*ii2; 
/* tie in the second period */
pa1=p1[i1]~p2[i2+1];
pa2=p1[i1+1]~p2[i2];


tpa1=p1k[i1]~p2k[i2+1];
tpa2=p1k[i1+1]~p2k[i2];

dtpa=tpa1-tpa2;
tau=dtpa[1]/(dtpa[1]-dtpa[2]);
if tau.lt 0; "tau < 0"; stop; endif;
if tau.gt 1; "tau > 1"; stop; endif;
tau1=tau/2;
tau2=(tau+1)/2;
points=points|(tau1*pa1+(1-tau1)*pa2)|(tau2*pa1+(1-tau2)*pa2); 

endif;

if (1-ii1)*(1-ii2); 
/* no tie */
points=points|(xp1[i1]~xp2[i2]); 
endif;

endfor;
endfor;
RETP(points[2:rows(points),.]);
ENDP;

PROC eli_poi(points,a1,b1,t1,a2,b2,t2,cen_t,pri);
/* eliminated unnecessary points */
local ppp, outc,tt1,tt2,tt,ii,outc1,outc2,cc;

/* x=0 */
tt1=points[.,1]; tt2=points[.,2];
tt=tt1~tt2;
if pri; "X=0";; tt; endif;
ii=(tt1.lt tt2)-(tt1.gt tt2);
tt=minc(tt');
cc=tt.ge cen_t;
tt=tt.*(1-cc);
ii=ii.*(1-cc);
tt=int(tt);
outc1=tt~ii~cc;


/* x=1 */
tt1=kinked_line(points[.,1],a1,b1,t1);
tt2=kinked_line(points[.,2],a2,b2,t2);
tt=tt1~tt2;
if pri; "X=0";; tt; endif;
ii=(tt1.lt tt2)-(tt1.gt tt2);
tt=minc(tt');
cc=tt.ge cen_t;
tt=tt.*(1-cc);
ii=ii.*(1-cc);
tt=int(tt);
outc2=tt~ii~cc;

outc=outc1~outc2;

if pri.eq 2; format /rd 7,3; points~outc; endif;
ppp=points[1,.];
outc2=outc[1,.];

for ii (2,rows(points),1);
outc1=outc2-outc[ii,.];
outc1=sumc((outc1').^2);
if minc(outc1).gt 0; 
ppp=ppp|points[ii,.];
outc2=outc2|outc[ii,.];
endif;
endfor;

RETP(ppp);
ENDP;

PROC find_poi_w_k(a1,a2,t1,b1,b2,t2,maxt);
/* 
a1 and a2 are the two slopes for waiting time 1
b1 and b2 are the two slopes for waiting time 2
*/
local pri,s,points;
pri=0;
s=seqa(1,1,maxt);
points=fi_poi(s,a1,a2,t1,b1,b2,t2);
points=eli_poi(points,a1,a2,t1,b1,b2,t2,maxt,pri);
RETP(points);
ENDP;

proc func_2mult_wk(a1,a2,t1,b1,b2,t2,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= find_poi_w_k(a1,a2,t1,b1,b2,t2,maxt);

ppp=points;
coef_x0=findcoefs(ppp,maxt);
ppp=kinked_line(points[.,1],a1,a2,t1)~kinked_line(points[.,2],b1,b2,t2);

coef_x1=findcoefs(ppp,maxt);

aa=coef_x0[.,1:cols(coef_x0)-1]~coef_x1[.,1:cols(coef_x1)-1]; 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 getstart(h1,h2,startt,maxt,breakt);
local probs,pp1,pp2, p,ftol,maxsec,maxit,b,f1,f2,f3,t1,t2,bopt,ib,ic;

format /rd 6,2;
t1=breakt;
t2=breakt;

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

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

bo_ptrue=pp1|pp2;
sumc(bo_ptrue);
bo_dat=zeros(3,1);
bo_dat[1]=t1; bo_dat[2]=t2; bo_dat[3]=maxt;

b=0|0|0|0;

bopt=b;
f1=func1(bopt);

format /rd 9,5;
for ic (1,4,1);
" IC";;
print /flush ic;;
print /flush f1;
for ib (1,60,1);
b[ic]=-0.2+ib/100;
f2=func1(b);

if f2.lt f1; bopt=b; f1=f2; print /flush f2~exp(b)'; endif;
endfor;
b=bopt;
endfor;



p=0.2*(eye(4)-0.25)|zeros(1,4);
p=p+bopt';
ftol=1.0e-3;
maxsec=240;
maxit=100;
{b,f1,f2,f3}= AMOEBA(p,ftol,maxsec,maxit,&func1);
b=b[1,.]';

bopt=b;
p=b;
f1=func1(bopt);

format /rd 9,5;
for ic (1,4,1);
" IC";;
print /flush ic;;
print /flush f1;
for ib (1,21,1);
b[ic]=p[ic]-0.15+3*(ib-1)/200;
f2=func1(b);

if f2.lt f1; bopt=b; f1=f2; print /flush f2~exp(b)'; endif;
endfor;
b=bopt;
endfor;



retp(exp(bopt));
endp;

proc func1(b);
b=exp(b);
retp(-func_2mult_wk(b[1],b[2],bo_dat[1],b[3],b[4],bo_dat[2],bo_dat[3],bo_ptrue,2));
endp;