/* C mex init_pot for in @jtree_sparse_inf_engine directory               */
/* The file enter_evidence.m in directory @jtree_sparse_inf_engine call it*/

/**************************************/
/* init_pot.c has 6 input & 2 output  */
/* engine                             */
/* clqs                               */
/* pots                               */
/* pot_type                           */
/* onodes                             */
/* ndx                                */
/*                                    */
/* clpot                              */
/* seppot                             */
/**************************************/
#include <math.h>
#include <search.h>
#include "mex.h"

int compare(const void* src1, const void* src2){
	int i1 = *(int*)src1 ;
	int i2 = *(int*)src2 ;
	return i1-i2 ;
}

void ind_subv(int index, const int *cumprod, int n, int *bsubv){
	int i;

	for (i = n-1; i >= 0; i--) {
		bsubv[i] = ((int)floor(index / cumprod[i]));
		index = index % cumprod[i];
	}
}

int subv_ind(const int n, const int *cumprod, const int *subv){
	int i, index=0;

	for(i=0; i<n; i++){
		index += subv[i] * cumprod[i];
	}
	return index;
}

void compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){
	int i, j;
	int *eff_cumprod, *subv, *diffsize, *diff_cumprod;

	subv = malloc(diffdim * sizeof(int));
	eff_cumprod = malloc(diffdim * sizeof(int));
	diffsize = malloc(diffdim * sizeof(int));
	diff_cumprod = malloc(diffdim * sizeof(int));
	for(i=0; i<diffdim; i++){
		eff_cumprod[i] = bCumprod[dmask[i]];
		diffsize[i] = (int)pbSize[dmask[i]];
	}
	diff_cumprod[0] = 1;
	for(i=0; i<diffdim-1; i++){
		diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i];
	}
	for(i=0; i<ND; i++){
		ind_subv(i, diff_cumprod, diffdim, subv);
		weight[i] = 0;
		for(j=0; j<diffdim; j++){
			weight[i] += eff_cumprod[j] * subv[j];
		}
	}
	free(eff_cumprod);
	free(subv);
	free(diffsize);
	free(diff_cumprod);
}

mxArray* convert_to_sparse(const double *table, const int NB, const int counts){
	mxArray  *spTable;
    int      i, k, *ir, *jc;
    double   *sr;
    
	spTable = mxCreateSparse(NB, 1, counts, mxREAL);
    sr = mxGetPr(spTable);
    ir = mxGetIr(spTable);
    jc = mxGetJc(spTable);

    k = 0; 
	jc[0] = 0;
	jc[1] = counts;
	for(i=0; i<NB; i++){
		if(table[i] != 0.0){
			sr[k] = table[i];
			ir[k] = i;
			k++;
		}
    }

	return spTable;
}

mxArray* convert_table_to_sparse(const double *bT, const int *index, const int nzCounts, const int NB){
	mxArray  *spTable;
    int      i, *irs, *jcs;
    double   *sr;
    
	spTable = mxCreateSparse(NB, 1, nzCounts, mxREAL);
    sr  = mxGetPr(spTable);
    irs = mxGetIr(spTable);
    jcs = mxGetJc(spTable);

	jcs[0] = 0;
	jcs[1] = nzCounts;

	for(i=0; i<nzCounts; i++){
			sr[i] = bT[i];
			irs[i] = index[i];
    }
	return spTable;	
}

mxArray* convert_ill_table_to_sparse(const double *bigTable, const int *sequence, const int nzCounts, const int NB){
	mxArray *spTable;
	int     i, temp, *irs, *jcs, count=0;
	double  *sr;

	spTable = mxCreateSparse(NB, 1, nzCounts, mxREAL);
    sr  = mxGetPr(spTable);
    irs = mxGetIr(spTable);
    jcs = mxGetJc(spTable);

	jcs[0] = 0;
	jcs[1] = nzCounts;

	for(i=0; i<nzCounts; i++){
		irs[i] = sequence[count];
		count++;
		temp = sequence[count];
		sr[i] = bigTable[temp];
		count++;
	}
	return spTable;
}

void multiply_null_by_fuPot(mxArray *bigPot, const mxArray *smallPot){
	int     i, j, count, NB, NS, siz_b, siz_s, ndim, nzCounts=0;
	int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *jc;
	double  *pbDomain, *psDomain, *pbSize, *psSize, *bTable, *sTable, value;
	mxArray *pTemp, *pTemp1;

	pTemp = mxGetField(bigPot, 0, "domain");
	pbDomain = mxGetPr(pTemp);
	siz_b = mxGetNumberOfElements(pTemp);
	pTemp = mxGetField(smallPot, 0, "domain");
	psDomain = mxGetPr(pTemp);
	siz_s = mxGetNumberOfElements(pTemp);

	pTemp = mxGetField(bigPot, 0, "sizes");
	pbSize = mxGetPr(pTemp);
	pTemp = mxGetField(smallPot, 0, "sizes");
	psSize = mxGetPr(pTemp);

	NB = 1;
	for(i=0; i<siz_b; i++){
		NB *= (int)pbSize[i];
	}
	NS = 1;
	for(i=0; i<siz_s; i++){
		NS *= (int)psSize[i];
	}

	pTemp = mxGetField(smallPot, 0, "T");
	sTable = mxGetPr(pTemp);
	bTable = malloc(NB * sizeof(double));
	for(i=0; i<NB; i++){
		bTable[i] = 0;
	}

	if(NS == 1){
		value = *sTable;
		for(i=0; i<NB; i++){
			bTable[i] = value;
		}
		nzCounts = NB;
		pTemp = mxGetField(bigPot, 0, "T");
		if(pTemp)mxDestroyArray(pTemp);
		pTemp = convert_to_sparse(bTable, NB, NB);
		mxSetField(bigPot, 0, "T", pTemp);
		free(bTable);
		return;
	}

	if(NS == NB){
		for(i=0; i<NB; i++){
			bTable[i] = sTable[i];
			if(sTable[i] != 0) nzCounts++;
		}
		pTemp = mxGetField(bigPot, 0, "T");
		if(pTemp)mxDestroyArray(pTemp);
		pTemp = convert_to_sparse(bTable, NB, nzCounts);
		mxSetField(bigPot, 0, "T", pTemp);
		free(bTable);
		return;
	}

	mask = malloc(siz_s * sizeof(int));
	count = 0;
	for(i=0; i<siz_s; i++){
		for(j=0; j<siz_b; j++){
			if(psDomain[i] == pbDomain[j]){
				mask[count] = j;
				count++;
				break;
			}
		}
	}
	
	ndim = siz_b;
	sx = (int *)malloc(sizeof(int)*ndim);
	sy = (int *)malloc(sizeof(int)*ndim);
	for(i=0; i<ndim; i++){
		sx[i] = (int)pbSize[i];
		sy[i] = 1;
	}
	for(i=0; i<count; i++){
		sy[mask[i]] = sx[mask[i]];
	}

	s = (int *)malloc(sizeof(int)*ndim);
	*(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
	subs =   (int *)malloc(sizeof(int)*ndim);
	cpsy2 =  (int *)malloc(sizeof(int)*ndim);
	for(i = 0; i < ndim; i++){
		subs[i] = 0;
		s[i] = sx[i] - 1;
	}
			
	for(i = 0; i < ndim-1; i++){
		cpsy[i+1] = cpsy[i]*sy[i]--;
		cpsy2[i] = cpsy[i]*sy[i];
	}
	cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);

	for(j=0; j<NB; j++){
		bTable[j] = *sTable;
		if(*sTable != 0.0) nzCounts++;
		for(i = 0; i < ndim; i++){
			if(subs[i] == s[i]){
				subs[i] = 0;
				if(sy[i])
					sTable -= cpsy2[i];
			}
			else{
				subs[i]++;
				if(sy[i])
					sTable += cpsy[i];
				break;
			}
		}
	}

	pTemp = mxGetField(bigPot, 0, "T");
	if(pTemp)mxDestroyArray(pTemp);
	pTemp = convert_to_sparse(bTable, NB, nzCounts);
	mxSetField(bigPot, 0, "T", pTemp);
	pTemp1 = mxGetField(bigPot, 0, "T");
	jc = mxGetJc(pTemp1);

	free(sx);
	free(sy);
	free(s);
	free(cpsy);
	free(subs);
	free(cpsy2);
    free(mask);
	free(bTable);
}

void multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){
	int     i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0;
	int     *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight;
	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr;
	mxArray *pTemp, *pTemp1;

	pTemp = mxGetField(bigPot, 0, "domain");
	pbDomain = mxGetPr(pTemp);
	bdim = mxGetNumberOfElements(pTemp);
	pTemp = mxGetField(smallPot, 0, "domain");
	psDomain = mxGetPr(pTemp);
	sdim = mxGetNumberOfElements(pTemp);

	pTemp = mxGetField(bigPot, 0, "sizes");
	pbSize = mxGetPr(pTemp);
	pTemp = mxGetField(smallPot, 0, "sizes");
	psSize = mxGetPr(pTemp);

	NB = 1;
	for(i=0; i<bdim; i++){
		NB *= (int)pbSize[i];
	}
	NS = 1;
	for(i=0; i<sdim; i++){
		NS *= (int)psSize[i];
	}
	ND = NB / NS;

	if(ND == 1){
		pTemp = mxGetField(bigPot, 0, "T");
		if(pTemp)mxDestroyArray(pTemp);
		pTemp1 = mxGetField(smallPot, 0, "T");
		pTemp = mxDuplicateArray(pTemp1);
		mxSetField(bigPot, 0, "T", pTemp);
		return;
	}

	pTemp = mxGetField(smallPot, 0, "T");
	spr = mxGetPr(pTemp);
	sir = mxGetIr(pTemp);
	sjc = mxGetJc(pTemp);
	NZS = sjc[1];

	NZB = ND * NZS;

	diffdim = bdim - sdim;