
/***********************************************************
 *Mediation Analysis                                       *
 *                                                         *
 *   Author: Zhiyong Zhang                                 *
 *                                                         *
 *   Purpose: Test mediation through bootstrap and EM      *
 *            algorithm for missing data                   *
 *                                                         *
 *   Usage:                                                *
 *      Run the program and the message appears            *
 ***********************************************************/

/***********************************************************
Change history:
2008/12/20 
   Version 1 based on the BootMed

 **********************************************************/
 

// The library newmat was used for the matrix operation    //
#define WANT_STREAM
#define WANT_MATH
#define WANT_FSTREAM

#include "newmatap.h"
#include "newmatio.h"
#include "newmatnl.h"
#include "newran.h"
#include "math.h"

#include
#include
#include
#include
#include "time.h"

#ifdef use_namespace
using namespace RBD_LIBRARIES;
#endif

/*------------------------------------------------------*\ 
 * Normal quantile calculation function                 *
\*------------------------------------------------------*/ 
double dfQNorm(double alpha){
	ColumnVector ck(100);
	ck=0.0;
	ck(1)=1;
	ck(2)=1;
	for (int k=3; k<=100; k++){
		for (int m=1; m> X(i,j);
		}
	fin.close();
	return X;
}

/*------------------------------------------------------*\ 
 * Function to organize the data using missing pattern  *
 * Return the matrix of data and vector of pattern size *
 * Input: matData: original data                        *
 *        iN: sample size                               *
 * Output: cvN: a vector with pattern size              *
 *         matPatData: sorted data with missing pattern *
\*------------------------------------------------------*/ 
void pattern(Matrix matData, int iN, ColumnVector &cvN, Matrix &matPatData){
	int n1=0,n2=0,n3=0,n4=0,n5=0,n6=0,n7=0,n8=0;
	double dmiss=99999.0;
	//Get the dimension of different patterns
	for (int i1=1; i1<=iN; i1++){
		if (matData(i1,1)!=dmiss){
			if (matData(i1,2)!=dmiss){
				if (matData(i1,3)!=dmiss){
					n1=n1+1;
				}else{
					n2=n2+1;
				}
			}else{
				if (matData(i1,3)!=dmiss){
					n3=n3+1;
				}else{
					n4=n4+1;
				}
			}
		}else{
			if (matData(i1,2)!=dmiss){
				if (matData(i1,3)!=dmiss){
					n5=n5+1;
				}else{
					n6=n6+1;
				}
			}else{
				if (matData(i1,3)!=dmiss){
					n7=n7+1;
				}else{
					n8=n8+1;
				}
			}
		}
	}//end of the sample size for each pattern

	//Define the matrix to store the data for different patterns
	int m1=n1, m2=n2, m3=n3, m4=n4, m5=n5, m6=n6, m7=n7, m8=n8;
		   
	if (n1==0){m1=1;}
	if (n2==0){m2=1;}
	if (n3==0){m3=1;}
	if (n4==0){m4=1;}
	if (n5==0){m5=1;}
	if (n6==0){m6=1;}
	if (n7==0){m7=1;}
	if (n8==0){m8=1;}

	Matrix matPat1(iN,3);
	Matrix matPat2(iN,3);
	Matrix matPat3(iN,3);
	Matrix matPat4(iN,3);
	Matrix matPat5(iN,3);
	Matrix matPat6(iN,3);
	Matrix matPat7(iN,3);
	Matrix matPat8(iN,3);

	matPat1=0.0;
	matPat2=0.0;
	matPat3=0.0;
	matPat4=0.0;
	matPat5=0.0;
	matPat6=0.0;
	matPat7=0.0;
	matPat8=0.0;

	//Get the dimension of different patterns
	int k1=1, k2=1, k3=1, k4=1, k5=1, k6=1, k7=1, k8=1;
	for (int i2=1; i2<=iN; i2++){
		if (matData(i2,1)!=dmiss){
			if (matData(i2,2)!=dmiss){
				if (matData(i2,3)!=dmiss){
					matPat1.row(k1)=matData.row(i2);
					k1++;
				}else{
					matPat2.row(k2)=matData.row(i2);
					k2++;
				}
			}else{
				if (matData(i2,3)!=dmiss){
					matPat3.row(k3)=matData.row(i2);
					k3++;
				}else{
					matPat4.row(k4)=matData.row(i2);
					k4++;
				}
			}
		}else{
			if (matData(i2,2)!=dmiss){
				if (matData(i2,3)!=dmiss){
					matPat5.row(k5)=matData.row(i2);
					k5++;
				}else{
					matPat6.row(k6)=matData.row(i2);
					k6++;
				}
			}else{
				if (matData(i2,3)!=dmiss){
					matPat7.row(k7)=matData.row(i2);
					k7++;
				}else{
					matPat8.row(k8)=matData.row(i2);
					k8++;
				}

			}
		}
	}//end of the data pattern
	//put the data into the pattern size vector
	cvN<0){
		for (i3=1; i3<=n1; i3++){
			matPatData.row(i3)=matPat1.row(i3);
	    }
	}
	if (n2>0){
		for (i3=1; i3<=n2; i3++){
			matPatData.row(n1+i3)=matPat2.row(i3);
	    }
	}
	if (n3>0){
		for (i3=1; i3<=n3; i3++){
			matPatData.row(l2+i3)=matPat3.row(i3);
	    }
	}
	if (n4>0){
		for (i3=1; i3<=n4; i3++){
			matPatData.row(l3+i3)=matPat4.row(i3);
	    }
	}
	if (n5>0){
		for (i3=1; i3<=n5; i3++){
			matPatData.row(l4+i3)=matPat5.row(i3);
	    }
	}
	if (n6>0){
		for (i3=1; i3<=n6; i3++){
			matPatData.row(l5+i3)=matPat6.row(i3);
	    }
	}
	if (n7>0){
		for (i3=1; i3<=n7; i3++){
			matPatData.row(l6+i3)=matPat7.row(i3);
	    }
	}
	if (n8>0){
		for (i3=1; i3<=n8; i3++){
			matPatData.row(l7+i3)=matPat8.row(i3);
	    }
	}

}

/*------------------------------------------------------*\ 
 * Function to calculate the mean and covariance matrix *
 * Output: mCov
           mCov[1:c,r] is the the means
		   mCov[1:c,1:r-1] is the covariance
\*------------------------------------------------------*/ 
Matrix mfCov(Matrix matX, int iRow, int iCol){
	RowVector vMean = matX.sum_columns()/iRow;
	Matrix mXbar(iRow,iCol);
	for (int i=1; i<=iRow; i++){
		mXbar.row(i)=vMean;
	}
    Matrix mDesign=matX-mXbar;
	Matrix mCov(iCol+1, iCol);
	mCov.submatrix(1,iCol,1,iCol)=mDesign.t()*mDesign/(iRow-1.0);
	mCov.row(iCol+1)=vMean;

	return mCov;
}

/*------------------------------------------------------*\ 
 * Function to estimate the coefficients                *
\*------------------------------------------------------*/

RowVector cvfParas(Matrix mCov)
{
	double d_a1, d_a2, d_b1, d_sig1, d_sig2, d_sigx, d_a2b1, dybar, dmbar, dxbar;
	double d_sy, d_sx, d_sm, d_sxy, d_sxm, d_sym;
	double d_mM, d_mY;
	RowVector cvPara(8);

	d_sy=mCov(3,3);
	d_sx=mCov(1,1);
	d_sm=mCov(2,2);
	d_sxy=mCov(1,3);
	d_sym=mCov(2,3);
	d_sxm=mCov(1,2);

	dybar=mCov(4,3);
	dxbar=mCov(4,1);
	dmbar=mCov(4,2);

	if (d_sxm*d_sxm-d_sm*d_sx==0){
		cout<<"Problem in calculation"<0.0000001){
		ssCov=mfEstep(matData, mTCov, iN);
		mT2Cov=mfMstep(ssCov,iN);
		mDiff=mT2Cov-mTCov;
		e=mDiff.sum_absolute_value();
		mTCov=mT2Cov;
	}
	cvPara=cvfParas(mTCov);

	return cvPara;
}

/*------------------------------------------------------*\ 
 * Delete a row from a matrix                           *
\*------------------------------------------------------*/ 
Matrix mDelARow(Matrix matData, int iR, int iN){
	Matrix mDel(iN-1, 3);
	if (iR==1){
		mDel=matData.submatrix(2,iN,1,3);
	}else if (iR==iN){
		mDel=matData.submatrix(1,iN-1,1,3);
	}else{
		mDel.submatrix(1,iR-1,1,3)=matData.submatrix(1,iR-1,1,3);
		mDel.submatrix(iR, iN-1,1,3)=matData.submatrix(iR+1,iN,1,3);
	}
	return mDel;
} 

/*------------------------------------------------------*\ 
 * Calculate the acceleration factor using Jackknife    *
\*------------------------------------------------------*/ 
RowVector rvJackA(Matrix matData, Matrix mCov,int iN){
	Matrix mJackPara(iN, 8), mDel(iN-1,3), mTemp1(iN,8), mTemp2(iN,8);
	RowVector dA(8), dMean(8), dSum1(8), dSum2(8);
	int i,j;
	for (i=1;i<=iN;i++){
		mDel = mDelARow(matData,i,iN);
		mJackPara.row(i)=cvfEM(mDel,mCov,iN-1);
	}

	dMean = mJackPara.sum_columns()/iN;
	
	for (i=1; i<=iN; i++){
		for (j=1; j<=8;j++){
			mTemp1(i,j) = (dMean(j)-mJackPara(i,j))*(dMean(j)-mJackPara(i,j))*(dMean(j)-mJackPara(i,j));
			mTemp2(i,j) = (dMean(j)-mJackPara(i,j))*(dMean(j)-mJackPara(i,j));
		}
	}

	dSum1=mTemp1.sum_columns();
	dSum2=mTemp2.sum_columns();
	for (i=1; i<=8; i++){
		dA(i)=dSum1(i)/(pow(dSum2(i),1.5)*6.0);
	}
	return dA;
}
/*------------------------------------------------------*\ 
 * Function to bootstrap                                *
\*------------------------------------------------------*/
Matrix mfBootOnce(Matrix matData, int iRow){
	Matrix mBoot(iRow, 3);
	
	Uniform U;
	int intvalue;

    for (int i=1; i<=iRow; i++){
		intvalue=int(U.Next()*iRow)+1;
		mBoot.row(i)=matData.row(intvalue);
	}
	return mBoot;
}

/*------------------------------------------------------*\ 
 * The stratified bootstrap                             *
\*------------------------------------------------------*/
Matrix mfBootstrap(Matrix matData, ColumnVector cvN, int iN){

	Matrix mStratBoot(iN,3);
	int n1=(int)cvN(1);
	int n2=(int)cvN(2);
	int n3=(int)cvN(3);
	int n4=(int)cvN(4);
	int n5=(int)cvN(5);
	int n6=(int)cvN(6);
	int n7=(int)cvN(7);
	int m2=n1+n2;
	int m3=m2+n3;
	int m4=m3+n4;
	int m5=m4+n5;
	int m6=m5+n6;
	int m7=m6+n7;

	if (n1>0){
		mStratBoot.submatrix(1,n1,1,3)=mfBootOnce(matData.submatrix(1,n1,1,3),n1);
	}
	if (n2>0){
		mStratBoot.submatrix(n1+1,m2,1,3)=mfBootOnce(matData.submatrix(n1+1,m2,1,3),n2);
	}
	if (n3>0){
		mStratBoot.submatrix(m2+1,m3,1,3)=mfBootOnce(matData.submatrix(m2+1,m3,1,3),n3);
	}
	if (n4>0){
		mStratBoot.submatrix(m3+1,m4,1,3)=mfBootOnce(matData.submatrix(m3+1,m4,1,3),n4);
	}
	if (n5>0){
		mStratBoot.submatrix(m4+1,m5,1,3)=mfBootOnce(matData.submatrix(m4+1,m5,1,3),n5);
	}
	if (n6>0){
		mStratBoot.submatrix(m5+1,m6,1,3)=mfBootOnce(matData.submatrix(m5+1,m6,1,3),n6);
	}
	if (n7>0){
		mStratBoot.submatrix(m6+1,m7,1,3)=mfBootOnce(matData.submatrix(m6+1,m7,1,3),n7);
	}

	return mStratBoot;
}

/*------------------------------------------------------*\ 
 * Function to calculate the CI for one parameter       *
\*------------------------------------------------------*/
RowVector rvfCI(ColumnVector rvBootRes, double dThetaHat, double dA, int iB, double alpha){
	RowVector rvCI(2);
	int i;
	double iPos=0.0;

	sort_ascending(rvBootRes);

	for (i=1; i<=iB; i++){
		if (rvBootRes(i)>output;
	char *fout;
	fout=(char*)output.c_str();
	std::ofstream out;
	
	out.open(fout);

// Print out the file information	
	out <<"----------------------------------------------------------\n";
	out <<"| Program name: MedCI.exe (V3.0)                         |\n";
	out <<"| By Zhiyong Zhang & Lijuan Wang                         |\n";
	out <<"| Email johnnyzhz@gmail.com for questions and comments   |\n";
	out <<"| See readme.txt for more information                    |\n";
	out <<"----------------------------------------------------------\n";
	out <<"\n";
	out <<"\n";
	out <<"******************************************************************\n";
	out <<"*                             NOTICE                             *\n";
	out <<"* Although this program has been tested to perform as expected   *\n";
	out <<"* in many cases, we cannot guarantee its performance under       *\n";
	out <<"* all possible circumstances.                                    *\n";
	out <<"* Users can use it for free at their own risks.                  *\n";
	out <<"******************************************************************\n";
	out <<"\n\n";

	//time
	time_t rawtime;
	time(&rawtime);
    out <<"The BootMed program is ran on "<> filename;
	out << filename <<"\n";
	char *fstr;
    fstr= (char*)filename.c_str();
	size(fstr,r1, c1);
    //Decide if number of variables > 3
	if (c1>3){
		cout<<"The number of variables >3 !";
		cout<<"Only three varialbes (Y, X, Mediation) are allowed in the data file!";
		out<<"The number of variables >3 !\n";
		out<<"Only three varialbes (Y, X, Mediation) are allowed in the data file!";
		exit(1);
	}


//Read the data into a matrix
	Matrix INPORG=readdata(fstr,r1,c1);

	

	//request the input of the alpha level

	std::string stralpha;

	cout <<"\nPlease specify the alpha level (a number between 0 and 1 only):"<>stralpha;
	
    double alpha= atof(stralpha.c_str());
	//cout <=1 || alpha<=0){
		cout << "Only a number bigger than 0 and less than 1 is accepted for the seed!"<>strseed;
	out <=1 || seed<=0){
		cout << "Only a number bigger than 0 and less than 1 is accepted for the seed!"<0){
		cout <<"What do you expect if no data observed for all three variables?"<>strBoot;

    int iBoot=atoi(strBoot.c_str());
	out << strBoot<<"\n";

	//Running time
	clock_t tStTime, tEndTime;
	tStTime=clock();

	//cout<