/***********************************************************
 *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<iostream>
#include<fstream>
#include<ostream>
#include<string>
#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<k; m++){
			ck(k)=ck(k)+ck(m)*ck(k-m)/((m)*(2.0*m-1.0));
		}
	}
 
	double x;
	x=2.0*alpha-1.0;
 
	double pi=3.14159265358979323846;
 
	double erf=0.0;
	for (int i=1; i<=100; i++){
	   erf=erf+ck(i)/(2.0*i-1.0)*pow((sqrt(pi)/2.0*x),(2.0*i-1.0));
	}
 
	erf=sqrt(2.0)*erf;
 
	return erf;
} 
 
// factorial fucntion
double ifFactor(int N){
	double iFactor=1.0;
	for (int i=1; i<=N; i++){
		iFactor = iFactor * i;
	}
	return iFactor;
}
/*------------------------------------------------------*\ 
 * Normal cumulative distribution function              *
\*------------------------------------------------------*/ 
double dfPNorm(double dX){
 
	double dZ = dX/sqrt(2.0);
 
	double pi=3.14159265358979323846;
 
    double erf=dZ;
	for (int i=1; i<=100; i++){
	   erf=erf+pow(-1.0,i)*pow(dZ, 2.0*i+1.0)/(ifFactor(i)*(2.0*i+1.0));
	}
	erf=2.0*erf/sqrt(pi);
	erf=(1.0+erf)/2.0;
 
	return erf;
} 
 
/*------------------------------------------------------*\ 
 * Function to read in a data file and decide the size. *
 * Read in the file name and output the row and column  *
 * Numbers.                                             *
\*------------------------------------------------------*/ 
 
void size(char *filename, int &rn, int &cn)
{
	Tracer tr("size"); // for tracking exceptions
    //Open file and decide the size of the file
	std::ifstream in;
	in.open(filename);
	if (!in) { std::cout << "cannot open the data file\n"; exit(1); }
 
	std::string str;
    std::string buf;
 
	 // Decide the number of rows
	int row=0;
 
    while(std::getline(in,str))
    {
        if(row==0)
            buf=str;
        row++;
	//	std::cout<<str<<"\n";
    }
 
	// Number of columns
	//c_str() convert a string to char*
 
    int column=0;
    const char *c=buf.c_str();
    bool tag=true;
    while(*c!='\0')
    {
        if(*c!=' '&&tag==true)
    {
        c++;
        tag=false;
        column++;
    }
    else if(*c++==' ')
    {
            tag=true;
    }
    }
 
	rn=row;
    cn=column;
 
	in.close();
 
}
 
/*------------------------------------------------------*\ 
 * Read the data from a file into a matrix              *
\*------------------------------------------------------*/ 
 
Matrix readdata(char *filename, int row, int col){
	std::ifstream fin;
	fin.open(filename);
 
	if (!fin) { std::cout << "cannot open the data file\n"; exit(1); }
 
	Matrix X(row,col);
 
	for (int i=1; i<=row; i++)
		for (int j=1; j<=col; j++)
		{
			fin >> 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<<n1<<n2<<n3<<n4<<n5<<n6<<n7<<n8;
 
	//put the data into the pattern matrix
	int l2=n1+n2;
	int l3=l2+n3;
	int l4=l3+n4;
	int l5=l4+n5;
	int l6=l5+n6;
	int l7=l6+n7;
 
	int i3;
	if (n1>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"<<endl;
		exit(1);
	}else{
 
	d_a1=(d_sxm*d_sym-d_sm*d_sxy)/(d_sxm*d_sxm-d_sm*d_sx);
	d_a2=(d_sx*d_sym-d_sxm*d_sxy)/(d_sm*d_sx-d_sxm*d_sxm);
	d_b1=d_sxm/d_sx;
	d_a2b1=d_a2*d_b1;
 
	d_sigx=d_sx;
	d_sig1=d_sm-d_sxm*d_sxm/d_sx;
	d_sig2=(d_sm*d_sx*d_sy-d_sxm*d_sxm*d_sy-d_sx*d_sym*d_sym-d_sm*d_sxy*d_sxy+2*d_sxm*d_sym*d_sxy)/(d_sm*d_sx-d_sxm*d_sxm);
 
	d_mM = dmbar - dxbar*d_b1;
	d_mY = dybar - d_a1*dxbar - d_a2*dmbar;
 
	cvPara(1)=d_mM;
	cvPara(2)=d_mY;
	cvPara(3)=d_b1;
	cvPara(4)=d_a2;
	cvPara(5)=d_a1; //variance for residuals of M~X
	cvPara(6)=d_sig1; //variance for residuals of Y~X+M
	cvPara(7)=d_sig2;
	cvPara(8)=d_a2b1;
	}
 
	return cvPara;
}
 
/*------------------------------------------------------*\ 
 * EM algorithm - The M-step
 * Output: mCov -- covariance matrix and means
\*------------------------------------------------------*/ 
Matrix mfMstep(Matrix mSS, int iN){
	Matrix mSSV(3,3);
	RowVector mSSM(3);
	mSSV = mSS.submatrix(1,3,1,3);
	mSSM = mSS.row(4);
 
	RowVector mMean=mSSM/iN;
	Matrix mCov(4,3);
	mCov.submatrix(1,3,1,3) = mSSV/(iN-1.0) - iN/(iN-1.0)*mMean.t()*mMean;
	mCov.row(4)=mMean;
 
	return mCov;
}
 
/*------------------------------------------------------*\ 
 * EM algorithm - The E-step
 * Output: mCov -- covariance matrix and means
\*------------------------------------------------------*/ 
 
Matrix mfEstep(Matrix matData, Matrix mCov, int iN){
	ColumnVector sX2(iN), sM2(iN),sY2(iN),sXM(iN),sXY(iN),sMY(iN),mX(iN),mM(iN),mY(iN);
	int i;
	double dmiss=99999.0;
	//Get the dimension of different patterns
	for (i=1; i<=iN; i++){
		if (matData(i,1)!=dmiss){
			if (matData(i,2)!=dmiss){
				if (matData(i,3)!=dmiss){
					mX(i)=matData(i,1);
					mM(i)=matData(i,2);
					mY(i)=matData(i,3);
					sX2(i)=mX(i)*mX(i);
					sY2(i)=mY(i)*mY(i);
					sM2(i)=mM(i)*mM(i);
					sXM(i)=mX(i)*mM(i);
					sXY(i)=mX(i)*mY(i);
					sMY(i)=mM(i)*mY(i);
				}else{
					mX(i)=matData(i,1);
					mM(i)=matData(i,2);
					mY(i)=mCov(4,3) + ( (mCov(3,1)*mCov(2,2)-mCov(3,2)*mCov(1,2))*(mX(i)-mCov(4,1))+
						                (mCov(3,2)*mCov(1,1)-mCov(3,1)*mCov(1,2))*(mM(i)-mCov(4,2)))/(
										 mCov(1,1)*mCov(2,2)-mCov(1,2)*mCov(1,2));
					sX2(i)=mX(i)*mX(i);
					sY2(i)=mY(i)*mY(i)+ mCov(3,3)-( (mCov(3,1)*mCov(2,2)-mCov(3,2)*mCov(1,2))*mCov(3,1)+
						                (mCov(3,2)*mCov(1,1)-mCov(3,1)*mCov(1,2))*mCov(3,2))/(
										 mCov(1,1)*mCov(2,2)-mCov(1,2)*mCov(1,2));
					sM2(i)=mM(i)*mM(i);
					sXM(i)=mX(i)*mM(i);
					sXY(i)=mX(i)*mY(i);
					sMY(i)=mM(i)*mY(i);
				}
			}else{
				if (matData(i,3)!=dmiss){
					mX(i)=matData(i,1);
					mY(i)=matData(i,3);
					mM(i)=mCov(4,2)+( (mCov(2,1)*mCov(3,3)-mCov(3,2)*mCov(1,3))*(mX(i)-mCov(4,1))+
						                (mCov(3,2)*mCov(1,1)-mCov(3,1)*mCov(1,2))*(mY(i)-mCov(4,3)))/(
										 mCov(1,1)*mCov(3,3)-mCov(1,3)*mCov(1,3));
					sX2(i)=mX(i)*mX(i);
					sY2(i)=mY(i)*mY(i);
					sM2(i)=mM(i)*mM(i)+mCov(2,2)-( (mCov(2,1)*mCov(3,3)-mCov(3,2)*mCov(1,3))*mCov(2,1)+
						                (mCov(3,2)*mCov(1,1)-mCov(3,1)*mCov(1,2))*mCov(2,3))/(
										 mCov(1,1)*mCov(3,3)-mCov(1,3)*mCov(1,3));
					sXM(i)=mX(i)*mM(i);
					sXY(i)=mX(i)*mY(i);
					sMY(i)=mM(i)*mY(i);
				}else{
					mX(i)=matData(i,1);
					mM(i)=mCov(4,2)+mCov(1,2)/mCov(1,1)*(mX(i)-mCov(4,1));
					mY(i)=mCov(4,3)+mCov(1,3)/mCov(1,1)*(mX(i)-mCov(4,1));
					sX2(i)=mX(i)*mX(i);
					sY2(i)=mY(i)*mY(i)+mCov(3,3)-mCov(1,3)*mCov(1,3)/mCov(1,1);
					sM2(i)=mM(i)*mM(i)+mCov(2,2)-mCov(1,2)*mCov(1,2)/mCov(1,1);
					sXM(i)=mX(i)*mM(i);
					sXY(i)=mX(i)*mY(i);
					sMY(i)=mM(i)*mY(i)+mCov(2,3)-mCov(1,2)*mCov(1,3)/mCov(1,1);
				}
			}
		}else{
			if (matData(i,2)!=dmiss){
				if (matData(i,3)!=dmiss){
					mM(i)=matData(i,2);
					mY(i)=matData(i,3);
					mX(i)=mCov(4,1) + ( (mCov(2,1)*mCov(3,3)-mCov(3,2)*mCov(1,3))*(mM(i)-mCov(4,2))+
						                (mCov(1,3)*mCov(2,2)-mCov(3,2)*mCov(1,2))*(mY(i)-mCov(4,3)))/(
										 mCov(2,2)*mCov(3,3)-mCov(2,3)*mCov(2,3));
					sX2(i)=mX(i)*mX(i) + mCov(1,1)-( (mCov(2,1)*mCov(3,3)-mCov(3,2)*mCov(1,3))*(mCov(1,2))+
						                (mCov(1,3)*mCov(2,2)-mCov(3,2)*mCov(1,2))*(mCov(1,3)))/(
										 mCov(2,2)*mCov(3,3)-mCov(2,3)*mCov(2,3));
					sY2(i)=mY(i)*mY(i);
					sM2(i)=mM(i)*mM(i);
					sXM(i)=mX(i)*mM(i);
					sXY(i)=mX(i)*mY(i);
					sMY(i)=mM(i)*mY(i);
				}else{
					mM(i)=matData(i,2);
					mX(i)=mCov(4,1)+mCov(1,2)/mCov(2,2)*(mM(i)-mCov(4,2));
					mY(i)=mCov(4,3)+mCov(2,3)/mCov(2,2)*(mM(i)-mCov(4,2));
					sX2(i)=mX(i)*mX(i)+mCov(1,1)-mCov(1,2)*mCov(1,2)/mCov(2,2);
					sY2(i)=mY(i)*mY(i)+mCov(3,3)-mCov(3,2)*mCov(3,2)/mCov(2,2);
					sM2(i)=mM(i)*mM(i);
					sXM(i)=mX(i)*mM(i);
					sXY(i)=mX(i)*mY(i)+mCov(1,3)-mCov(1,2)*mCov(2,3)/mCov(2,2);
					sMY(i)=mM(i)*mY(i);
				}
			}else{
				if (matData(i,3)!=dmiss){
					mY(i)=matData(i,3);
					mX(i)=mCov(4,1)+mCov(1,3)/mCov(3,3)*(mY(i)-mCov(4,3));
					mM(i)=mCov(4,2)+mCov(2,3)/mCov(3,3)*(mY(i)-mCov(4,3));
					sX2(i)=mX(i)*mX(i)+mCov(1,1)-mCov(1,3)*mCov(1,3)/mCov(3,3);
					sY2(i)=mY(i)*mY(i);
					sM2(i)=mM(i)*mM(i)+mCov(2,2)-mCov(2,3)*mCov(2,3)/mCov(3,3);
					sXM(i)=mX(i)*mM(i)+mCov(1,2)-mCov(1,3)*mCov(2,3)/mCov(3,3);
					sXY(i)=mX(i)*mY(i);
					sMY(i)=mM(i)*mY(i);
				}
			}
		}
	}//end of the sample size for each pattern
	//Get the sum of each vector and put them into a matrix
	Matrix mSS(4,3);
	mSS(1,1)=sX2.sum();
	mSS(1,2)=sXM.sum();
	mSS(1,3)=sXY.sum();
	mSS(2,1)=mSS(1,2);
	mSS(2,2)=sM2.sum();
	mSS(2,3)=sMY.sum();
	mSS(3,1)=mSS(1,3);
	mSS(3,2)=mSS(2,3);
	mSS(3,3)=sY2.sum();
 
	mSS(4,1)=mX.sum();
	mSS(4,2)=mM.sum();
	mSS(4,3)=mY.sum();
 
	return mSS;
}
 
/*------------------------------------------------------*\ 
 * EM algorithm                                         *
\*------------------------------------------------------*/ 
RowVector cvfEM(Matrix matData, Matrix mCov, int iN){
	RowVector cvPara;
	double e=1.0;
	Matrix ssCov(4,3);
	Matrix mTCov(4,3);
	mTCov=mCov;
	Matrix mT2Cov(4,3),mDiff(4,3);
 
	while (e>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)<dThetaHat){
			iPos=iPos+1.0;
		}
	}
	double dZ01=dfQNorm(iPos/iB);
	double zalpha1=dfQNorm((1.0-alpha)/2.0);
	double tt1 = dfPNorm(dZ01+(dZ01+zalpha1)/(1-dA*(dZ01+zalpha1)));
 
	int iLow = int(tt1*iB+.5);
	rvCI(1)=rvBootRes(iLow);
 
	double zalpha2=dfQNorm((1.0+alpha)/2.0);
	double tt2 = dfPNorm(dZ01+(dZ01+zalpha2)/(1-dA*(dZ01+zalpha2)));
 
	int iUp = int(tt2*iB+.5);
	rvCI(2)=rvBootRes(iUp);
 
	return rvCI;
}
/*------------------------------------------------------*\ 
 * Function to calculate the CIs                        *
\*------------------------------------------------------*/ 
Matrix mBootCI(Matrix matData, Matrix mCov, ColumnVector cvN, int iN, int iB,double alpha){
	RowVector rvTrue=cvfEM(matData,mCov,iN);
	RowVector rvTemp(2), tempres(8);
	Matrix mBootData(iN,3),mBootRes(iB,8), mCIres(8,3);
    int k=int(iB/10);
	int i,j=1;
	cout<<"Running: ";
	for (i=1; i<=iB; i++){
		mBootData=mfBootstrap(matData,cvN,iN);
		tempres=cvfEM(mBootData,mCov,iN);
		mBootRes.row(i)=tempres;
		if (i==j*k){
			cout<<j*10<<"% ";
			j++;
		}
	}
	cout<<"DONE"<<"\n";
	//output the mean of the results
	//cout<<mBootRes.sum_columns()/iB;
 
	RowVector rvA=rvJackA(matData,mCov,iN);
 
	for (i=1; i<=8; i++){
		rvTemp = rvfCI(mBootRes.column(i),rvTrue(i),rvA(i),iB,alpha);
		mCIres(i,1)=rvTrue(i);
		mCIres(i,2)=rvTemp(1);
		mCIres(i,3)=rvTemp(2);	
	}
 
	return mCIres;
}
 
 
/*------------------------------------------------------*\ 
 * The main program                                     *
\*------------------------------------------------------*/ 
int main()
{
// Print the program information
    cout <<"\n";
	cout <<"\n";
	cout <<"---------------------------------------------------------"<<endl;
	cout <<"| Program name: MedCI.exe (V3.0)                        |"<<endl;
	cout <<"| By Zhiyong Zhang & Lijuan Wang                        |"<<endl;
	cout <<"| Email johnnyzhz@gmail.com for questions and comments  |"<<endl;
	cout <<"| See readme.txt for more information                   |"<<endl;
	cout <<"---------------------------------------------------------"<<endl;
	cout <<"\n";
	cout <<"\n";
 
	cout <<"******************************************************************"<<endl;
	cout <<"*                             NOTICE                             *"<<endl;
	cout <<"* Although this program has been tested to perform as expected   *"<<endl;
	cout <<"* in many cases, we cannot guarantee its performance under       *"<<endl;
	cout <<"* all possible circumstances.                                    *"<<endl;
	cout <<"* Users can use it for free at their own risks.                  *"<<endl;
	cout <<"******************************************************************"<<endl;
	cout <<endl<<endl;
 
// Ask the input of the output file
	cout <<"Please input the output file name: "<<endl;
	std::string output;
	std::cin >>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 "<<ctime(&rawtime)<<"\n";
 
	out << "\nPlease input the output file name: "<<"\n";
	out << output << "\n";
 
//Input the data file name and check the data	
	int c1;
	int r1;
	std::string filename;
    cout <<"\nPlease input the file name of the data:"<<endl;
	out <<"\nPlease input the file name of the data:"<<"\n";
    std::cin >> 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):"<<endl;
	out <<"\nPlease specify the alpha level (a number between 0 and 1 only):\n";
 
	std::cin >>stralpha;
 
    double alpha= atof(stralpha.c_str());
	//cout <<alpha <<endl;
	out <<alpha<<"\n";
	if (alpha>=1 || alpha<=0){
		cout << "Only a number bigger than 0 and less than 1 is accepted for the seed!"<<endl;
		out <<"Only a number bigger than 0 and less than 1 is accepted for the seed!\n";
		exit(1);
	}
 
 
	//request the input of seed for random number generation
 
	std::string strseed;
 
	cout <<"\nPlease specify the random number seed (a number between 0 and 1 only):"<<endl;
	out <<"\nPlease specify the random number seed (a number between 0 and 1 only):\n";
	std::cin >>strseed;
	out <<strseed<<"\n";
 
	double seed= atof(strseed.c_str());
	//cout <<seed;
 
	if (seed>=1 || seed<=0){
		cout << "Only a number bigger than 0 and less than 1 is accepted for the seed!"<<endl;
		out <<"Only a number bigger than 0 and less than 1 is accepted for the seed!\n";
		exit(1);
	}
	MotherOfAll urng(seed);                    // declare uniform random number generator
    Random::Set(urng); 
    ColumnVector cvN(8);
    Matrix matPatData(r1,c1);
    pattern(INPORG,r1,cvN,matPatData);
 
	if (cvN(1)==0){
		cout <<"Sorry, you asked too much for missing data analysis! There is no single case with complete data!"<<endl;
		out <<"Sorry, you asked too much for missing data analysis!\n There is no single case with complete data!\n"<<endl;
		exit(1);
	}
 
	if (cvN(8)>0){
		cout <<"What do you expect if no data observed for all three variables?"<<endl;
		out <<"What do you expect if no data observed for all three variables?\n"<<endl;
		exit(1);
	}
 
	std::string strBoot;
 
	cout <<"Please input the bootstrap sample size:"<<endl;
	out <<"Please input the bootstrap sample size:"<< "\n";
	std::cin>>strBoot;
 
    int iBoot=atoi(strBoot.c_str());
	out << strBoot<<"\n";
 
	//Running time
	clock_t tStTime, tEndTime;
	tStTime=clock();
 
	//cout<<matPatData;
	Matrix mCov=mfCov(matPatData.submatrix(1,cvN(1),1,3),cvN(1),3);
 
	Matrix test=mBootCI(matPatData,mCov,cvN,r1,iBoot,alpha);
 
	//end time
	tEndTime=clock();
	out <<"\n The total running time is "<<(double)(tEndTime-tStTime)/CLOCKS_PER_SEC<<" seconds. \n";
	cout<<test;
 
   out.close();
 
   return 0;
 
}