#include <oxstd.h>
#include <oxprob.h>
#include <oxdraw.h>

like(theta,x,T){
  decl stage=theta[0],alpha1=theta[1],alpha2=theta[2],alpha3=theta[3],d1=theta[4],d2=theta[5];
  decl n=rows(x),i,sum1=0,sum2=0,sum3=0,n1=0,n2=0,n3=0;
  decl C,like;

 if (stage==3) {//&& 
			  sum1=0;sum2=0;sum3=0;n1=0;n2=0;n3=0;
			  for (i=0;i<n;i++){
				 if (x[i]<=d1)                {	sum1=sum1+log(1-(x[i]/T));n1=n1+1;}
				 if (x[i]>d1 && x[i]<=T-d2)   {	sum2=sum2+log(1-(x[i]/T));n2=n2+1;}
				 if (x[i]>T-d2 )              {	sum3=sum3+log(1-(x[i]/T));n3=n3+1;}	}
			  
			  C = (alpha1*alpha2*alpha3/T) /
			      (   alpha3*(alpha1-alpha2)*((1-d1/T)^alpha2)
				    + alpha2*alpha3*((1-d1/T)^(alpha2-alpha1))
				    + ((d2/T)^alpha2)*alpha1*(alpha2-alpha3)   );
			
			  like = n*log(C) + n1*(alpha2-alpha1)*log(1-d1/T)
			                  + n3*(alpha2-alpha3)*log(d2/T)
							  + (alpha1-1)*sum1
							  + (alpha2-1)*sum2
							  + (alpha3-1)*sum3;
			  }//$$
if (stage==2) {//&&
              sum1=0;sum2=0;sum3=0;n1=0;n2=0;n3=0;
			  for (i=0;i<n;i++){
				 if (x[i]<=T-d2)              {	sum1=sum1+log(1-x[i]/T);n1=n1+1;}
				 if (x[i]>T-d2)               {	sum2=sum2+log(1-x[i]/T);n2=n2+1;}  }
			  like = n*log(alpha2)
			        -n*log(T)
					+(alpha2-1)*sum1
					+(alpha3-1)*sum2
					+n2*(alpha2-alpha3)*log(d2/T)
					-n*log(1 - (1-alpha2/alpha3)*(d2/T)^alpha2);	
              }//$$
if (stage==1) {//&&
              sum1=0;sum2=0;sum3=0;
			  for (i=0;i<n;i++){sum1=sum1+log(1-x[i]/T);}
			  like = n*log(alpha3/T) +(alpha3-1)*sum1;	
              }//$$

return like; }

//************************************************************************************//
main(){
decl data = new array[12], vT;
decl T=7, x = loadmat("Palm7Days.txt"), n=rows(x), simdata=0;
decl stop,count,l,m,i,o;
decl theta=zeros(1,6),numparam=columns(theta),hlike;
decl parent,cut,cross1,cross2,mutepoint;
decl maxcount=500,maxcross=50,maxpop=maxcross*2,muterate=0.1;
decl offspring=zeros(maxpop,numparam),popfit=zeros(maxpop,numparam+1);
decl minalp1,maxalp1,minalp2,maxalp2,minalp3,maxalp3,mind1,maxd1,mind2,maxd2;
decl AIC,BIC;
decl outfit;
decl file=fopen("GA-Estimates-AllData.txt","w");

data[0] = loadmat("Palm3Days.txt");
data[1] = loadmat("Palm5Days.txt");
data[2] = loadmat("Palm7Days.txt");
data[3] = loadmat("Cartier3Days.txt");
data[4] = loadmat("Cartier5Days.txt");
data[5] = loadmat("Cartier7Days.txt");
data[6] = loadmat("HarryPotter3Days.txt");
data[7] = loadmat("HarryPotter5Days.txt");
data[8] = loadmat("HarryPotter7Days.txt");
data[9] = loadmat("XBox3Days.txt");
data[10] = loadmat("XBox5Days.txt");
data[11] = loadmat("XBox7Days.txt");

vT = <3;5;7;3;5;7;3;5;7;3;5;7>;

file = fopen("GA-Estimates-AllData.txt", "a");
fprint(file,"  Order of Data: Palm, Cartier, HarryPotter, XBox  ","\n");
fprint(file,"  Length / Like / AIC / BIC / Stages / Params  ","\n");
fclose(file);

//LOOP OVER DATASETS
for (o=0;o<rows(vT);o++){
x = data[o]; T = vT[o];	n=rows(x);

//Generic Values for all Data
//if (simdata==0){
//minalp1=4;maxalp1=6;
//minalp2=0.2;maxalp2=0.6;
//minalp3=0.8;maxalp3=1.2;
//mind1=1;maxd1=2.5;
//mind2=0;maxd2=0.001;	}
if (simdata==0){
minalp1=1;   maxalp1=15;
minalp2=0.1; maxalp2=1;
minalp3=0.5; maxalp3=15;
mind1=1;     maxd1=5;
mind2=0;     maxd2=0.01;	}

//--Randomly initialize regression param's--//
for (l=0;l<maxpop;l++){
 offspring[l][0]    =  ceil(3* ranu(1,1));
 offspring[l][1]    =  minalp1 + (maxalp1-minalp1) * ranu(1,1);
 offspring[l][2]    =  minalp2 + (maxalp2-minalp2) * ranu(1,1);
 offspring[l][3]    =  minalp3 + (maxalp3-minalp3) * ranu(1,1);
 offspring[l][4]    =  mind1   + (maxd1-mind1) * ranu(1,1);
 offspring[l][5]    =  mind2   + (maxd2-mind2) * ranu(1,1);
}
		
//---Inner Loop	---//
stop=1;count=0;
while ( count<maxcount ) {
//-- evaluate population fitness--//
for (l=1;l<=maxpop;l++){   
 theta  = offspring[l-1][];
 hlike  =  like(theta,x,T);
 //calculate AIC and BIC model selection criterion
 BIC    = (-hlike + (log(n)*numparam/2))/n;
 AIC    = (-hlike + numparam)/n;
 hlike  = - BIC; //Notice: want to MINIMIZE BIC or AIC; so we MAXIMIZE its negative value!!
 popfit[l-1][] = hlike~theta;}
//select two members of population with largest fitness
popfit = sortbyc(popfit,0);  
parent = popfit[maxpop-10:][1:];	
stop = norm( (popfit[maxpop-1][0]-popfit[0][0]) / popfit[maxpop-1][0] );
//perform crossover
for (i=0;i<=maxcross-1;i++){
 //choose cutting point for crossover
 cross1 = floor(10*ranu(1,1));
 cross2 = floor(10*ranu(1,1));
 cut    = floor((numparam-1)*ranu(1,1)); 
 offspring[2*i][0:cut]   = parent[cross1][0:cut];
 offspring[2*i+1][0:cut] = parent[cross2][0:cut];
 offspring[2*i][cut+1:]  = parent[cross1][cut+1:];
 offspring[2*i+1][cut+1:]= parent[cross2][cut+1:];		   
}  
//perform mutations
for (i=0;i<=maxpop-1;i++){
 if (ranu(1,1)<muterate){ 
   mutepoint = floor(numparam*ranu(1,1));
   if (mutepoint==0) {offspring[i][mutepoint]= ceil(3* ranu(1,1));}
   if (mutepoint==1) {offspring[i][mutepoint]= minalp1 + (maxalp1-minalp1) * ranu(1,1);}
   if (mutepoint==2) {offspring[i][mutepoint]= minalp2 + (maxalp2-minalp2) * ranu(1,1);}
   if (mutepoint==3) {offspring[i][mutepoint]= minalp3 + (maxalp3-minalp3) * ranu(1,1);}
   if (mutepoint==4) {offspring[i][mutepoint]= mind1   + (maxd1-mind1) * ranu(1,1);}
   if (mutepoint==5) {offspring[i][mutepoint]= mind2   + (maxd2-mind2) * ranu(1,1);}   
  }}
count  = count + 1;
 theta = popfit[maxpop-1][1:];
 hlike  =  like(theta,x,T);
 BIC    = (-hlike + (log(n)*numparam/2))/n;
 AIC    = (-hlike + numparam)/n;
outfit = T ~ hlike ~ AIC ~ BIC ~ theta;
 println("Iter ",count," Estimates ",outfit);
}//---End Inner Loop---//
file = fopen("GA-Estimates-AllData.txt", "a");
fprint(file,outfit,"\n");
fclose(file);
println("*********************************************************");
println("Estimates are ",outfit);
println("*********************************************************");
}//END OF LOOP OVER DATASETS
}//END OF PROGRAM