/* C mex for collect_evidence.c in @jtree_sparse_inf_engine directory */
/* File enter_evidence.m in directory @jtree_sparse_inf_engine call it*/

/******************************************/
/* collect_evidence has 3 input & 2 output*/
/* engine                                 */
/* clpot                                  */
/* seppot                                 */
/*                                        */
/* clpot                                  */
/* seppot                                 */
/******************************************/

#include <math.h>
#include <stdlib.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);
}

void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){
	double *ptr;
	void   *newptr;
	int    *ir, *jc;
	int    nbytes;

	if(new_nzmax == old_nzmax) return;
	nbytes = new_nzmax * sizeof(*ptr);
	ptr = mxGetPr(spArray);
	newptr = mxRealloc(ptr, nbytes);
	mxSetPr(spArray, newptr);
	nbytes = new_nzmax * sizeof(*ir);
	ir = mxGetIr(spArray);
	newptr = mxRealloc(ir, nbytes);
	mxSetIr(spArray, newptr);
	jc = mxGetJc(spArray);
	jc[0] = 0;
	jc[1] = new_nzmax;
	mxSetNzmax(spArray, new_nzmax);
}

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

	spTable = mxCreateSparse(N, 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] = Table[temp];
		count++;
	}
	return spTable;
}

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, *bpr;
	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);

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

	NB = 1;
	for(i=0; i<bdim; i++){
		NB *= (int)pbSize[i];
	}

	if(sdim == 0){
		pTemp = mxCreateSparse(NB, 1, NB, mxREAL);
		mxSetField(bigPot, 0, "T", pTemp);
		bpr = mxGetPr(pTemp);
		sir = mxGetIr(pTemp);
		sjc = mxGetJc(pTemp);
		sjc[0] = 0;
		sjc[1] = NB;
		for(i=0; i<NB; i++){
			bpr[i] = *spr;
			sir[i] = i;
		}	
		return;
	}

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

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


	NZB = ND * NZS;

	diffdim = bdim - sdim;
	sequence = malloc(NZB * 2 * sizeof(int));
	bigTable = malloc(NZB * sizeof(double));
	samemask = malloc(sdim * sizeof(int));
	diffmask = malloc(diffdim * sizeof(int));
	bCumprod = malloc(bdim * sizeof(int));
	sCumprod = malloc(sdim * sizeof(int));
	weight = malloc(ND * sizeof(int));
	ssubv = malloc(sdim * sizeof(int));

	count = 0;
	count1 = 0;
	for(i=0; i<bdim; i++){
		match = 0;
		for(j=0; j<sdim; j++){
			if(pbDomain[i] == psDomain[j]){
				samemask[count] = i;
				match = 1;
				count++;
				break;
			}
		}
		if(match == 0){
			diffmask[count1] = i; 
			count1++;
		}
	}

	bCumprod[0] = 1;
	for(i=0; i<bdim-1; i++){
		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
	}
	sCumprod[0] = 1;
	for(i=0; i<sdim-1; i++){
		sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
	}

	count = 0;
	compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim);
	for(i=0; i<NZS; i++){
		sindex = sir[i];
		ind_subv(sindex, sCumprod, sdim, ssubv);
		temp = 0;
		for(j=0; j<sdim; j++){
			temp += ssubv[j] * bCumprod[samemask[j]];
		}
		for(j=0; j<ND; j++){
			bindex = weight[j] + temp;
			bigTable[nzCounts] = spr[i];
			sequence[count] = bindex;
			count++;
			sequence[count] = nzCounts;
			nzCounts++;
			count++;
		}
	}

	pTemp = mxGetField(bigPot, 0, "T");
	if(pTemp)mxDestroyArray(pTemp);
	qsort(sequence, nzCounts, sizeof(int) * 2, compare);
	pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB);
	mxSetField(bigPot, 0, "T", pTemp);

	free(sequence); 
	free(bigTable);
	free(samemask);
	free(diffmask);
	free(bCumprod);
	free(sCumprod);
	free(weight);
	free(ssubv);
}

void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){
	int     i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0;
	int     *mask, *result, *bir, *sir, *rir, *bjc, *sjc, *rjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
	double  *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, *rpr;
	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];
	}

	pTemp = mxGetField(bigPot, 0, "T");
	bpr = mxGetPr(pTemp);
	bir = mxGetIr(pTemp);
	bjc = mxGetJc(pTemp);
	NZB = bjc[1];

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

	if(sdim == 0){
		for(i=0; i<NZB; i++){
			bpr[i] *= *spr;
		}	
		return;
	}

	pTemp1 = mxCreateSparse(NB, 1, NZB, mxREAL);
	rpr = mxGetPr(pTemp1);
	rir = mxGetIr(pTemp1);
	rjc = mxGetJc(pTemp1);
	rjc[0] = 0;
	rjc[1] = NZB;

	mask = malloc(sdim * sizeof(int));
	bCumprod = malloc(bdim * sizeof(int));
	sCumprod = malloc(sdim * sizeof(int));
	bsubv = malloc(bdim * sizeof(int));
	ssubv = malloc(sdim * sizeof(int));

	count = 0;
	for(i=0; i<sdim; i++){
		for(j=0; j<bdim; j++){
			if(psDomain[i] == pbDomain[j]){
				mask[count] = j;
				count++;
				break;
			}
		}
	}
	
	bCumprod[0] = 1;
	for(i=0; i<bdim-1; i++){
		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
	}
	sCumprod[0] = 1;