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


unsigned int numnonzeros(const mxArray*);


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


	int  *ExpTransTktFile, *betaindxflag, *NumExperts, *edgeknowindx;
	double *betapara, *phipara, *delta, *sigmasquare, *mupara;

	double *d_expertnew, *wcost, *testindx, *TRCtran, *Cplussumfij;

	double  *tempprobsum, *temp,*tempprobsumsquare, *temp1, *temp2, *temp3, *temp_trc, *TRC_ij, *tempsigmasum;

	double expfactor;

	mwSize i, j, k, l, p;
	mwSize TotalTransSize, betarowsize, NumEdge, exp_row, exp_coln, NumTkts, NumExps, tktlength,explength, edgeknowcoln;
	mwSize *dims; 
	int tktindx, nextexpert, edgeknow_flag, maxtest, tempindx11, edgeknow_flag_ir;

	
			
	ExpTransTktFile=(int *)mxGetPr(prhs[0]); 	
	betapara=(double *)mxGetPr(prhs[1]); 
	phipara=(double *)mxGetPr(prhs[2]);	
	betaindxflag=(int *)mxGetPr(prhs[3]); 
	NumExperts=(int *)mxGetPr(prhs[4]); 
	TRCtran=(double *)mxGetPr(prhs[5]); 
	delta=(double *)mxGetPr(prhs[6]);  // after log transform
	sigmasquare=(double *)mxGetPr(prhs[7]);
	mupara=(double *)mxGetPr(prhs[8]);
	TRC_ij=(double *)mxGetPr(prhs[9]); 
	
	TotalTransSize=mxGetM(prhs[0]);
	betarowsize=mxGetM(prhs[1]); //the row size of betapara	
	

	tempprobsum=mxCreateDoubleMatrix(1,1,mxREAL);
	temp=mxCreateDoubleMatrix(1,1,mxREAL);
	temp1=mxCreateDoubleMatrix(1,1,mxREAL);
	temp2=mxCreateDoubleMatrix(1,1,mxREAL);
	temp3=mxCreateDoubleMatrix(1,1,mxREAL);
	temp_trc=mxCreateDoubleMatrix(1,1,mxREAL);
	Cplussumfij=mxCreateDoubleMatrix(1,1,mxREAL);
	tempsigmasum=mxCreateDoubleMatrix(1,1,mxREAL);

	tempprobsumsquare=mxCreateDoubleMatrix(1,1,mxREAL);



	dims = (mwSize *) mxMalloc (2 * sizeof(mwSize));  // This is very important !!!!	

	dims[0] = 2;
    dims[1] = 1;
    plhs[0]= mxCreateNumericArray (2, dims, mxDOUBLE_CLASS, mxREAL);
	d_expertnew=(double *) mxGetData (plhs[0]);

	dims[0] = 1;
    dims[1] = 1;
    plhs[1]= mxCreateNumericArray (2, dims, mxDOUBLE_CLASS, mxREAL);
	wcost=(double *) mxGetData (plhs[1]);

	

	//for (i=0;i<TotalTransSize;i++)
		for (i=0;i<1;i++)
	{
		tempprobsum[0]=0;
        temp_trc[0]=0;

		tempprobsumsquare[0]=0;
		tempsigmasum[0]=0;

		tktindx=ExpTransTktFile[2*TotalTransSize+i];

		if (betaindxflag[4*TotalTransSize+i]==0 && betaindxflag[3*TotalTransSize+i]!=0)
		{
			temp3[0]=betapara[3*betarowsize+betaindxflag[2*TotalTransSize+i]-1];
			wcost[0]=wcost[0]-phipara[2*TotalTransSize+i]*log(temp3[0]);
			//d_expertnew : no need to update mu and sigma
		}
		else if (betaindxflag[4*TotalTransSize+i]!=0)  
		{
		
		for (j=0;j<betaindxflag[TotalTransSize+i];j++)
		{
			if (betapara[4*betarowsize+betaindxflag[i]+j-1]==1)
			{
			edgeknow_flag=((int)betapara[2*betarowsize+betaindxflag[i]+j-1]-1)*NumExperts[0]+(int)betapara[betaindxflag[i]+j-1]-1;
			edgeknow_flag_ir=((int)betapara[2*betarowsize+betaindxflag[2*TotalTransSize+i]-1]-1)*NumExperts[0]+(int)betapara[betaindxflag[2*TotalTransSize+i]-1]-1;

			temp[0]=TRC_ij[edgeknow_flag];   //diffmatrix is of  the same size as betaparanew
			
			temp_trc[0]+=temp[0]*(delta[edgeknow_flag_ir]-delta[edgeknow_flag])/sigmasquare[0];

			tempprobsum[0]+=temp[0];

			tempsigmasum[0]+=temp[0]*(1/(2*sigmasquare[0]*sigmasquare[0]))*(delta[edgeknow_flag_ir]*delta[edgeknow_flag_ir]-delta[edgeknow_flag]*delta[edgeknow_flag]+2*mupara[0]*delta[edgeknow_flag]-2*mupara[0]*delta[edgeknow_flag_ir]);


			if(j==(betaindxflag[2*TotalTransSize+i]-betaindxflag[i]))
			{
				//betaindxflag[i]+j-1 the jth neighbor stored in betapara
				temp2[0]=temp[0];
				temp3[0]=betapara[3*betarowsize+betaindxflag[i]+j-1];  // 
			}
			}
			
		}

		tempprobsumsquare[0]=tempprobsum[0]*tempprobsum[0];
		temp_trc[0]=TRC_ij[edgeknow_flag_ir]*temp_trc[0]/tempprobsumsquare[0];
		tempsigmasum[0]=TRC_ij[edgeknow_flag_ir]*tempsigmasum[0]/tempprobsumsquare[0];

		d_expertnew[0]-=(phipara[2*TotalTransSize+i]/temp3[0])*temp_trc[0];
		d_expertnew[1]-=(phipara[2*TotalTransSize+i]/temp3[0])*tempsigmasum[0]*sigmasquare[0];
		//d_expertnew[1]=(delta[edgeknow_flag_ir]*delta[edgeknow_flag_ir]-delta[edgeknow_flag]*delta[edgeknow_flag]+2*mupara[0]*delta[edgeknow_flag]-2*mupara[0]*delta[edgeknow_flag_ir]);


		wcost[0]=wcost[0]-phipara[2*TotalTransSize+i]*log(temp3[0]);
		
		}
		else if (betaindxflag[3*TotalTransSize+i]==0)
		{		
		for (j=0;j<betaindxflag[TotalTransSize+i];j++)
		{
			
			edgeknow_flag=((int)betapara[2*betarowsize+betaindxflag[i]+j-1]-1)*NumExperts[0]+(int)betapara[betaindxflag[i]+j-1]-1;
			edgeknow_flag_ir=((int)betapara[2*betarowsize+betaindxflag[2*TotalTransSize+i]-1]-1)*NumExperts[0]+(int)betapara[betaindxflag[2*TotalTransSize+i]-1]-1;

			temp[0]=TRC_ij[edgeknow_flag];   //diffmatrix is of  the same size as betaparanew
			
			temp_trc[0]+=temp[0]*(delta[edgeknow_flag_ir]-delta[edgeknow_flag])/sigmasquare[0];

			tempprobsum[0]+=temp[0];
			tempsigmasum[0]+=temp[0]*(1/(2*sigmasquare[0]*sigmasquare[0]))*(delta[edgeknow_flag_ir]*delta[edgeknow_flag_ir]-delta[edgeknow_flag]*delta[edgeknow_flag]+2*mupara[0]*delta[edgeknow_flag]-2*mupara[0]*delta[edgeknow_flag_ir]);

			if(j==(betaindxflag[2*TotalTransSize+i]-betaindxflag[i]))
			{
				//betaindxflag[i]+j-1 the jth neighbor stored in betapara
				temp2[0]=temp[0];
				temp3[0]=betapara[3*betarowsize+betaindxflag[i]+j-1];  // 
			}		
			
		}

		tempprobsumsquare[0]=tempprobsum[0]*tempprobsum[0];
		temp_trc[0]=TRC_ij[edgeknow_flag_ir]*temp_trc[0]/tempprobsumsquare[0];
		
		tempsigmasum[0]=TRC_ij[edgeknow_flag_ir]*tempsigmasum[0]/tempprobsumsquare[0];

		d_expertnew[0]-=(phipara[2*TotalTransSize+i]/temp3[0])*temp_trc[0];
		d_expertnew[1]-=(phipara[2*TotalTransSize+i]/temp3[0])*tempsigmasum[0]*sigmasquare[0];

		wcost[0]=wcost[0]-phipara[2*TotalTransSize+i]*log(temp3[0]);
		
		}

		// the following is updating based on TRC only
		
		tempprobsum[0]=0;
		temp_trc[0]=0;
		tempsigmasum[0]=0;
		for (j=0;j<betaindxflag[TotalTransSize+i];j++)
		{
			
			edgeknow_flag=((int)betapara[2*betarowsize+betaindxflag[i]+j-1]-1)*NumExperts[0]+(int)betapara[betaindxflag[i]+j-1]-1;
			edgeknow_flag_ir=((int)betapara[2*betarowsize+betaindxflag[2*TotalTransSize+i]-1]-1)*NumExperts[0]+(int)betapara[betaindxflag[2*TotalTransSize+i]-1]-1;


			temp[0]=TRC_ij[edgeknow_flag];   //diffmatrix is of  the same size as betaparanew
			
			temp_trc[0]+=temp[0]*(delta[edgeknow_flag_ir]-delta[edgeknow_flag])/sigmasquare[0];

			tempprobsum[0]+=temp[0];

			tempsigmasum[0]+=temp[0]*(1/(2*sigmasquare[0]*sigmasquare[0]))*(delta[edgeknow_flag_ir]*delta[edgeknow_flag_ir]-delta[edgeknow_flag]*delta[edgeknow_flag]+2*mupara[0]*delta[edgeknow_flag]-2*mupara[0]*delta[edgeknow_flag_ir]);

			

			if(j==(betaindxflag[2*TotalTransSize+i]-betaindxflag[i]))
			{
				
				temp2[0]=temp[0];
				temp3[0]= TRCtran[edgeknow_flag_ir]; // 
			}		
			
		}

		tempprobsumsquare[0]=tempprobsum[0]*tempprobsum[0];
		temp_trc[0]=TRC_ij[edgeknow_flag_ir]*temp_trc[0]/tempprobsumsquare[0];
		tempsigmasum[0]=TRC_ij[edgeknow_flag_ir]*tempsigmasum[0]/tempprobsumsquare[0];

		d_expertnew[0]-=(phipara[3*TotalTransSize+i]/temp3[0])*temp_trc[0];
		d_expertnew[1]-=(phipara[3*TotalTransSize+i]/temp3[0])*tempsigmasum[0]*sigmasquare[0];
		
	    

		wcost[0]=wcost[0]-phipara[3*TotalTransSize+i]*log(temp3[0]);


		
	}

}