#include <math.h>
#include "mex.h"
#include "matrix.h"
#include <stdlib.h>
//#include<fstream>


unsigned int numnonzeros(const mxArray*);


void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray*prhs[]) { 

	double *phipara;
	double *gammaprior;

	double *betapara, *betapara_VE, *betapara_PR;
	int *ExpTransTktFile, *betaindxflag;
	double *VolumTran, *TRCtran;
	double *psigammarowsum;
	double *psigamma;
	double *purerandtrans;
	

	double *phiparanew, *gammaparanew, *term2, *term3, *term5;
	double *temp, *temp1, *tempbeta, *tempbeta0, *tempbeta1, *tempbeta2, *tempbeta3;
	mxArray *betawork;

	mwSize i, j, k, l, p, betaindx,m;
	mwSize TotalTransSize, phicolnsize, gammarowsize, gammacolnsize, betarowsize,NumExperts;
	mwSize *dims;

	phipara=(double *)mxGetPr(prhs[0]);
	gammaprior=(double *)mxGetPr(prhs[1]);
	ExpTransTktFile=(int *)mxGetPr(prhs[2]);
	VolumTran=(double *)mxGetPr(prhs[3]);
    betapara=(double *)mxGetPr(prhs[4]);
	psigamma=(double *)mxGetPr(prhs[5]);
	psigammarowsum=(double *)mxGetPr(prhs[6]);
	purerandtrans=(double *)mxGetPr(prhs[7]);
	betaindxflag=(int *)mxGetPr(prhs[8]);
	TRCtran=(double *)mxGetPr(prhs[9]);
	betapara_VE=(double *)mxGetPr(prhs[10]);
	betapara_PR=(double *)mxGetPr(prhs[11]);



	phicolnsize=mxGetN(prhs[0]);
	gammarowsize=mxGetM(prhs[1]);
	gammacolnsize=mxGetN(prhs[1]);
	TotalTransSize=mxGetM(prhs[2]);
	betarowsize=mxGetM(prhs[4]);
	NumExperts=mxGetM(prhs[3]);

	dims = (mwSize *) mxMalloc (2 * sizeof(mwSize));
    dims[0] = TotalTransSize;
    dims[1] = phicolnsize;
    plhs[0]= mxCreateNumericArray (2, dims, mxDOUBLE_CLASS, mxREAL);
	phiparanew=(double *) mxGetData (plhs[0]);

	dims[0] = gammarowsize;
    dims[1] = gammacolnsize;  // only for the oracle 
    plhs[1]= mxCreateNumericArray (2, dims, mxDOUBLE_CLASS, mxREAL);
	gammaparanew=(double *) mxGetData (plhs[1]);

	dims[0] = 1;
    dims[1] = 1;  // only for the oracle 
    plhs[2]= mxCreateNumericArray (2, dims, mxDOUBLE_CLASS, mxREAL);
	term2=(double *) mxGetData (plhs[2]);
	term2[0]=0;

	dims[0] = 1;
    dims[1] = 1;  // only for the oracle 
    plhs[3]= mxCreateNumericArray (2, dims, mxDOUBLE_CLASS, mxREAL);
	term3=(double *) mxGetData (plhs[3]);
	term3[0]=0;

	dims[0] = 1;
    dims[1] = 1;  // only for the oracle 
    plhs[4]= mxCreateNumericArray (2, dims, mxDOUBLE_CLASS, mxREAL);
	term5=(double *) mxGetData (plhs[4]);
	term5[0]=0;

	temp=mxCreateDoubleMatrix(1,1,mxREAL);
	temp1=mxCreateDoubleMatrix(1,1,mxREAL);
	tempbeta0=mxCreateDoubleMatrix(1,1,mxREAL);
	tempbeta1=mxCreateDoubleMatrix(1,1,mxREAL);
	tempbeta2=mxCreateDoubleMatrix(1,1,mxREAL);
	tempbeta3=mxCreateDoubleMatrix(1,1,mxREAL);

	//tempbeta2[0]=0.5;

	betawork=mxCreateDoubleMatrix(6,1,mxREAL);
	tempbeta=mxGetPr(betawork);
	

	temp[0]=0;
	//temp1[0]=0;

	for (i=0;i<gammarowsize;i++)
	{
		for (j=0;j<gammacolnsize;j++)
		{
			gammaparanew[j*gammarowsize+i]=0;
		}
	}

	for (i=0;i<TotalTransSize;i++)
	{
		phiparanew[i]=phipara[i];
		phiparanew[TotalTransSize+i]=phipara[TotalTransSize+i];
		temp1[0]=0;
		//tempbeta1[0]=VolumTran[NumExperts*(ExpTransTktFile[TotalTransSize+i]-1)+ExpTransTktFile[i]-1];
		tempbeta[3]=VolumTran[NumExperts*(ExpTransTktFile[TotalTransSize+i]-1)+ExpTransTktFile[i]-1];
		tempbeta[5]=purerandtrans[NumExperts*(ExpTransTktFile[TotalTransSize+i]-1)+ExpTransTktFile[i]-1];
		tempbeta[1]=TRCtran[NumExperts*(ExpTransTktFile[TotalTransSize+i]-1)+ExpTransTktFile[i]-1];

		/*for (j=0;j<betarowsize;j++)
		{
			if (((int)(betapara[betarowsize+j]))==ExpTransTktFile[TotalTransSize*2+i] && ((int)(betapara[2*betarowsize+j]))==ExpTransTktFile[i])
			{
				if (((int)(betapara[j]))==ExpTransTktFile[TotalTransSize+i])
				{
				//tempbeta0[0]=betapara[3*betarowsize+j];
				tempbeta[0]=betapara[3*betarowsize+j];
				//term2[0]=1;
				//tempbeta[0]=temp[0];
				break;
				}				
			    
			}			
		}*/

		tempbeta[0]=betapara[3*betarowsize+betaindxflag[TotalTransSize*2+i]-1];
		tempbeta[2]=betapara_VE[3*betarowsize+betaindxflag[TotalTransSize*2+i]-1];
		tempbeta[4]=betapara_PR[3*betarowsize+betaindxflag[TotalTransSize*2+i]-1];

		///*** speacial processing for this mixture model

		for (k=2;k<phicolnsize;k++)
		{
	    phiparanew[k*TotalTransSize+i]=tempbeta[k-2]*exp(psigamma[(k-2)*gammarowsize+ExpTransTktFile[i]-1]-psigammarowsum[ExpTransTktFile[i]-1]);

		temp1[0]+=phiparanew[k*TotalTransSize+i];
		}

		for (k=2;k<phicolnsize;k++)
		{
	    phiparanew[k*TotalTransSize+i]/=temp1[0];

		gammaparanew[(k-2)*gammarowsize+ExpTransTktFile[i]-1]+=phiparanew[k*TotalTransSize+i];

		term3[0]=term3[0]+log(tempbeta[k-2])*phipara[k*TotalTransSize+i];  // we use the old one to estimate
		term2[0]=term2[0]+phipara[k*TotalTransSize+i]*(psigamma[(k-2)*gammarowsize+ExpTransTktFile[i]-1]-psigammarowsum[ExpTransTktFile[i]-1]);
		term5[0]=term5[0]-phipara[k*TotalTransSize+i]*log(phipara[k*TotalTransSize+i]);
		}

		

	}
		//term2[0]=psigammarowsum[ExpTransTktFile[0]-1];
}