/*trial_estimation.cpp */
  
#include <math.h>
#include "trial_estimation.h"

/*define fixed parameters*/
void start_upcrmv(DVECTOR wdose1,DVECTOR wdose2,DVECTOR parm, DVECTOR priori_alpha,int *level1,int *level2,int n_r,int n_c)
{
    DVECTOR priori_p={-999,.05,.1,.2,.3,.5,.7};
	int i;

	for (i=1;i<=n_c;i++)
	{wdose2[i]=atanh(2*(1-priori_p[i])-1);
	 wdose1[i]=wdose2[i];
	}
		
	for (i=1;i<=3;i++)
		priori_alpha[i]=1;
}

/*define working model */
double hv(double dose1, double dose2, DVECTOR parm) /*start_upcrm*/
/*prior working model for probability of toxicity
  should be defined in starting up */
{ double p,p2,p3;
p2=parm[2];p3=parm[3];
p=1-pow((tanh (dose1) +1)/2,p2)*pow((tanh (dose2) +1)/2,p3);
if (p<0 || p>1)
{
	printf("\nerror from hv");
	printf("\ntanh dose1 %f, pow %f, tanh dose2 %f,pow %f",tanh (dose1),pow((tanh (dose1) +1)/2,p2),tanh (dose2),pow((tanh (dose2) +1)/2,p3));
	vprintd(parm,3);
	fprintf(stderr,"\nerror hv,dose1 %f,dose2 %f, p%f and parm",dose1,dose2,p);
	fprintf(stderr,"\np2 is %f and p2 is %f",p2,p3);
	exit(1);
}
return (p);
}

/***********************************************************/
void prep_crmv(DVECTOR crm_parm,DMATRIX toxicp,int *y,int *y_lev1,int *y_lev2,int tot_pid,int n_r,int n_c,double *r)
{
	static int info=0;
	static DVECTOR priori_alpha;
	static DVECTOR wdose1,wdose2;

	int lev1,lev2;
	if (info==0)
	{	start_upcrmv(wdose1,wdose2,crm_parm,priori_alpha,&lev1,&lev2,n_r,n_c);
		info=1;
	}
	crm_parm[1]=-999;//not used
//	printf("\nin prep crmv, alpha");
//	vprintd(priori_alpha,3);
//	vprintd(crm_parm,3);
	get_crmparm(y,y_lev1,y_lev2,wdose1,wdose2,tot_pid,priori_alpha,toxicp,n_r,n_c);
	*r=toxicp[y_lev1[tot_pid]][y_lev2[tot_pid]];	
	printf("\n\nAfter %d patients' responses are observed, the estimated toxicity matrix is:",tot_pid);
	mprintd(toxicp,n_r,n_c);
	fprintf(stderr,"\n\nAfter %d patients' responses are observed, the estimated toxicity matrix is:",tot_pid);
	mprintdout(toxicp,n_r,n_c,stderr);
}

/***********************************************************/
void get_crmparm(int *y, int *level1,int *level2,DVECTOR wdose1, DVECTOR wdose2,int pi,DVECTOR priori_alpha,DMATRIX toxicp,int n_r,int n_c)
/*	y: response 
	y_level : response dose level
	pi : patient id
	return next patient's dose level
*/
{ 
	compute_muparm(prod_phi,y,level1,level2,wdose1,wdose2,pi,priori_alpha,toxicp,n_r,n_c);
}

/***********************************************************/
void get_finaltoxicpv(DVECTOR parm,DVECTOR wdose1,DVECTOR wdose2,DMATRIX toxicp,int n_r, int n_c)
{ int i,j;
   for (i=1;i<=n_r;i++)
	for (j=1;j<=n_c;j++)
		toxicp[i][j]=hv(wdose1[i],wdose2[j],parm);
}

/***********************************************************/
void compute_muparm(long double f (IVECTOR,IVECTOR,IVECTOR,DVECTOR, DVECTOR,int,DVECTOR),int *y, \
int *level1,int *level2,DVECTOR wdose1,DVECTOR wdose2,int pi,DVECTOR priori_alpha,DMATRIX toxicp,int n_r,int n_c)
{
const int sample_n=5000; /*number of samples used in integration*/
//static double sample1[sample_n+1];
static double sample2 [sample_n +1];
static double sample3 [sample_n +1];

static int countsample=0,notready=1;
long double li,den;
long double num[4];
int i, j, k;
DVECTOR parm;
countsample+=1;
if (countsample%10==0)
	notready=1;
if (notready==1)
{	
	for (i=1;i<=sample_n;i++)
	{   
//		sample1[i]=RanExp(&seed,priori_alpha[1]);	
		sample2[i]=RanExp(&seed,priori_alpha[2]);	
		sample3[i]=RanExp(&seed,priori_alpha[3]);	
	}
	notready=0;
}

for (i=1;i<=4;i++)
   num[i]=0.0;

for (i=1;i<=n_r;i++)
   for (j=1;j<=n_c;j++)
	   toxicp[i][j]=0;
den=0;

for (i=1; i<=sample_n;i++)
{ 
//	parm[1]=sample1[i];
	parm[2]=sample2[i];
	parm[3]=sample3[i];
	
	li=f(y,level1,level2,wdose1,wdose2,pi,parm);
	den+=li;

	for (j=1;j<=n_r;j++)
		for (k=1;k<=n_c;k++)
			toxicp[j][k]+=li*hv(wdose1[j],wdose2[k],parm);
}
for (j=1;j<=n_r;j++)
		for (k=1;k<=n_c;k++)
          toxicp[j][k]=toxicp[j][k]/den;

if (den<=0)
{fprintf(stderr,"Crmv compute bottom is %f",den);
 exit(1);
}
}

/***********************************************************/
long double prod_phi(int *y,int *level1,int *level2,DVECTOR wdose1,DVECTOR wdose2,int pi, DVECTOR parm)
/* return product phi */
{
	int i;
	long double prod=30*log(10);
	for (i=1; i<=pi; i++)
		prod+=log(phi(y[i], wdose1[level1[i]],wdose2[level2[i]],parm));
	return ( expl(prod));
}
/***********************************************************/
double phi(int res, double dose1,double dose2, DVECTOR parm)
{
	if (res<=0 || res>=1)
	{
	if (res==0)
	return (1-hv(dose1,dose2,parm));
	else return (hv(dose1,dose2,parm));
	}
	else
	{fprintf(stderr,"\nerror in crmv, phi, res %d",res);
	exit(1);
	}

}

/***********************************************************/
double atanh(double x)
{
	if (x>1)
		x=1-1.0e-10;
	if (x<-1)
		x=-1+1.0e-10;
	return .5*log((1+x)/(1-x));
}