wolffd@0: /* C mex init_pot for in @jtree_sparse_inf_engine directory               */
wolffd@0: /* The file enter_evidence.m in directory @jtree_sparse_inf_engine call it*/
wolffd@0: 
wolffd@0: /**************************************/
wolffd@0: /* init_pot.c has 6 input & 2 output  */
wolffd@0: /* engine                             */
wolffd@0: /* clqs                               */
wolffd@0: /* pots                               */
wolffd@0: /* pot_type                           */
wolffd@0: /* onodes                             */
wolffd@0: /* ndx                                */
wolffd@0: /*                                    */
wolffd@0: /* clpot                              */
wolffd@0: /* seppot                             */
wolffd@0: /**************************************/
wolffd@0: #include <math.h>
wolffd@0: #include <search.h>
wolffd@0: #include "mex.h"
wolffd@0: 
wolffd@0: int compare(const void* src1, const void* src2){
wolffd@0: 	int i1 = *(int*)src1 ;
wolffd@0: 	int i2 = *(int*)src2 ;
wolffd@0: 	return i1-i2 ;
wolffd@0: }
wolffd@0: 
wolffd@0: void ind_subv(int index, const int *cumprod, int n, int *bsubv){
wolffd@0: 	int i;
wolffd@0: 
wolffd@0: 	for (i = n-1; i >= 0; i--) {
wolffd@0: 		bsubv[i] = ((int)floor(index / cumprod[i]));
wolffd@0: 		index = index % cumprod[i];
wolffd@0: 	}
wolffd@0: }
wolffd@0: 
wolffd@0: int subv_ind(const int n, const int *cumprod, const int *subv){
wolffd@0: 	int i, index=0;
wolffd@0: 
wolffd@0: 	for(i=0; i<n; i++){
wolffd@0: 		index += subv[i] * cumprod[i];
wolffd@0: 	}
wolffd@0: 	return index;
wolffd@0: }
wolffd@0: 
wolffd@0: void compute_fixed_weight(int *weight, const double *pbSize, const int *dmask, const int *bCumprod, const int ND, const int diffdim){
wolffd@0: 	int i, j;
wolffd@0: 	int *eff_cumprod, *subv, *diffsize, *diff_cumprod;
wolffd@0: 
wolffd@0: 	subv = malloc(diffdim * sizeof(int));
wolffd@0: 	eff_cumprod = malloc(diffdim * sizeof(int));
wolffd@0: 	diffsize = malloc(diffdim * sizeof(int));
wolffd@0: 	diff_cumprod = malloc(diffdim * sizeof(int));
wolffd@0: 	for(i=0; i<diffdim; i++){
wolffd@0: 		eff_cumprod[i] = bCumprod[dmask[i]];
wolffd@0: 		diffsize[i] = (int)pbSize[dmask[i]];
wolffd@0: 	}
wolffd@0: 	diff_cumprod[0] = 1;
wolffd@0: 	for(i=0; i<diffdim-1; i++){
wolffd@0: 		diff_cumprod[i+1] = diff_cumprod[i] * diffsize[i];
wolffd@0: 	}
wolffd@0: 	for(i=0; i<ND; i++){
wolffd@0: 		ind_subv(i, diff_cumprod, diffdim, subv);
wolffd@0: 		weight[i] = 0;
wolffd@0: 		for(j=0; j<diffdim; j++){
wolffd@0: 			weight[i] += eff_cumprod[j] * subv[j];
wolffd@0: 		}
wolffd@0: 	}
wolffd@0: 	free(eff_cumprod);
wolffd@0: 	free(subv);
wolffd@0: 	free(diffsize);
wolffd@0: 	free(diff_cumprod);
wolffd@0: }
wolffd@0: 
wolffd@0: mxArray* convert_to_sparse(const double *table, const int NB, const int counts){
wolffd@0: 	mxArray  *spTable;
wolffd@0:     int      i, k, *ir, *jc;
wolffd@0:     double   *sr;
wolffd@0:     
wolffd@0: 	spTable = mxCreateSparse(NB, 1, counts, mxREAL);
wolffd@0:     sr = mxGetPr(spTable);
wolffd@0:     ir = mxGetIr(spTable);
wolffd@0:     jc = mxGetJc(spTable);
wolffd@0: 
wolffd@0:     k = 0; 
wolffd@0: 	jc[0] = 0;
wolffd@0: 	jc[1] = counts;
wolffd@0: 	for(i=0; i<NB; i++){
wolffd@0: 		if(table[i] != 0.0){
wolffd@0: 			sr[k] = table[i];
wolffd@0: 			ir[k] = i;
wolffd@0: 			k++;
wolffd@0: 		}
wolffd@0:     }
wolffd@0: 
wolffd@0: 	return spTable;
wolffd@0: }
wolffd@0: 
wolffd@0: mxArray* convert_table_to_sparse(const double *bT, const int *index, const int nzCounts, const int NB){
wolffd@0: 	mxArray  *spTable;
wolffd@0:     int      i, *irs, *jcs;
wolffd@0:     double   *sr;
wolffd@0:     
wolffd@0: 	spTable = mxCreateSparse(NB, 1, nzCounts, mxREAL);
wolffd@0:     sr  = mxGetPr(spTable);
wolffd@0:     irs = mxGetIr(spTable);
wolffd@0:     jcs = mxGetJc(spTable);
wolffd@0: 
wolffd@0: 	jcs[0] = 0;
wolffd@0: 	jcs[1] = nzCounts;
wolffd@0: 
wolffd@0: 	for(i=0; i<nzCounts; i++){
wolffd@0: 			sr[i] = bT[i];
wolffd@0: 			irs[i] = index[i];
wolffd@0:     }
wolffd@0: 	return spTable;	
wolffd@0: }
wolffd@0: 
wolffd@0: mxArray* convert_ill_table_to_sparse(const double *bigTable, const int *sequence, const int nzCounts, const int NB){
wolffd@0: 	mxArray *spTable;
wolffd@0: 	int     i, temp, *irs, *jcs, count=0;
wolffd@0: 	double  *sr;
wolffd@0: 
wolffd@0: 	spTable = mxCreateSparse(NB, 1, nzCounts, mxREAL);
wolffd@0:     sr  = mxGetPr(spTable);
wolffd@0:     irs = mxGetIr(spTable);
wolffd@0:     jcs = mxGetJc(spTable);
wolffd@0: 
wolffd@0: 	jcs[0] = 0;
wolffd@0: 	jcs[1] = nzCounts;
wolffd@0: 
wolffd@0: 	for(i=0; i<nzCounts; i++){
wolffd@0: 		irs[i] = sequence[count];
wolffd@0: 		count++;
wolffd@0: 		temp = sequence[count];
wolffd@0: 		sr[i] = bigTable[temp];
wolffd@0: 		count++;
wolffd@0: 	}
wolffd@0: 	return spTable;
wolffd@0: }
wolffd@0: 
wolffd@0: void multiply_null_by_fuPot(mxArray *bigPot, const mxArray *smallPot){
wolffd@0: 	int     i, j, count, NB, NS, siz_b, siz_s, ndim, nzCounts=0;
wolffd@0: 	int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *jc;
wolffd@0: 	double  *pbDomain, *psDomain, *pbSize, *psSize, *bTable, *sTable, value;
wolffd@0: 	mxArray *pTemp, *pTemp1;
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "domain");
wolffd@0: 	pbDomain = mxGetPr(pTemp);
wolffd@0: 	siz_b = mxGetNumberOfElements(pTemp);
wolffd@0: 	pTemp = mxGetField(smallPot, 0, "domain");
wolffd@0: 	psDomain = mxGetPr(pTemp);
wolffd@0: 	siz_s = mxGetNumberOfElements(pTemp);
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "sizes");
wolffd@0: 	pbSize = mxGetPr(pTemp);
wolffd@0: 	pTemp = mxGetField(smallPot, 0, "sizes");
wolffd@0: 	psSize = mxGetPr(pTemp);
wolffd@0: 
wolffd@0: 	NB = 1;
wolffd@0: 	for(i=0; i<siz_b; i++){
wolffd@0: 		NB *= (int)pbSize[i];
wolffd@0: 	}
wolffd@0: 	NS = 1;
wolffd@0: 	for(i=0; i<siz_s; i++){
wolffd@0: 		NS *= (int)psSize[i];
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(smallPot, 0, "T");
wolffd@0: 	sTable = mxGetPr(pTemp);
wolffd@0: 	bTable = malloc(NB * sizeof(double));
wolffd@0: 	for(i=0; i<NB; i++){
wolffd@0: 		bTable[i] = 0;
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	if(NS == 1){
wolffd@0: 		value = *sTable;
wolffd@0: 		for(i=0; i<NB; i++){
wolffd@0: 			bTable[i] = value;
wolffd@0: 		}
wolffd@0: 		nzCounts = NB;
wolffd@0: 		pTemp = mxGetField(bigPot, 0, "T");
wolffd@0: 		if(pTemp)mxDestroyArray(pTemp);
wolffd@0: 		pTemp = convert_to_sparse(bTable, NB, NB);
wolffd@0: 		mxSetField(bigPot, 0, "T", pTemp);
wolffd@0: 		free(bTable);
wolffd@0: 		return;
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	if(NS == NB){
wolffd@0: 		for(i=0; i<NB; i++){
wolffd@0: 			bTable[i] = sTable[i];
wolffd@0: 			if(sTable[i] != 0) nzCounts++;
wolffd@0: 		}
wolffd@0: 		pTemp = mxGetField(bigPot, 0, "T");
wolffd@0: 		if(pTemp)mxDestroyArray(pTemp);
wolffd@0: 		pTemp = convert_to_sparse(bTable, NB, nzCounts);
wolffd@0: 		mxSetField(bigPot, 0, "T", pTemp);
wolffd@0: 		free(bTable);
wolffd@0: 		return;
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	mask = malloc(siz_s * sizeof(int));
wolffd@0: 	count = 0;
wolffd@0: 	for(i=0; i<siz_s; i++){
wolffd@0: 		for(j=0; j<siz_b; j++){
wolffd@0: 			if(psDomain[i] == pbDomain[j]){
wolffd@0: 				mask[count] = j;
wolffd@0: 				count++;
wolffd@0: 				break;
wolffd@0: 			}
wolffd@0: 		}
wolffd@0: 	}
wolffd@0: 	
wolffd@0: 	ndim = siz_b;
wolffd@0: 	sx = (int *)malloc(sizeof(int)*ndim);
wolffd@0: 	sy = (int *)malloc(sizeof(int)*ndim);
wolffd@0: 	for(i=0; i<ndim; i++){
wolffd@0: 		sx[i] = (int)pbSize[i];
wolffd@0: 		sy[i] = 1;
wolffd@0: 	}
wolffd@0: 	for(i=0; i<count; i++){
wolffd@0: 		sy[mask[i]] = sx[mask[i]];
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	s = (int *)malloc(sizeof(int)*ndim);
wolffd@0: 	*(cpsy = (int *)malloc(sizeof(int)*ndim)) = 1;
wolffd@0: 	subs =   (int *)malloc(sizeof(int)*ndim);
wolffd@0: 	cpsy2 =  (int *)malloc(sizeof(int)*ndim);
wolffd@0: 	for(i = 0; i < ndim; i++){
wolffd@0: 		subs[i] = 0;
wolffd@0: 		s[i] = sx[i] - 1;
wolffd@0: 	}
wolffd@0: 			
wolffd@0: 	for(i = 0; i < ndim-1; i++){
wolffd@0: 		cpsy[i+1] = cpsy[i]*sy[i]--;
wolffd@0: 		cpsy2[i] = cpsy[i]*sy[i];
wolffd@0: 	}
wolffd@0: 	cpsy2[ndim-1] = cpsy[ndim-1]*(--sy[ndim-1]);
wolffd@0: 
wolffd@0: 	for(j=0; j<NB; j++){
wolffd@0: 		bTable[j] = *sTable;
wolffd@0: 		if(*sTable != 0.0) nzCounts++;
wolffd@0: 		for(i = 0; i < ndim; i++){
wolffd@0: 			if(subs[i] == s[i]){
wolffd@0: 				subs[i] = 0;
wolffd@0: 				if(sy[i])
wolffd@0: 					sTable -= cpsy2[i];
wolffd@0: 			}
wolffd@0: 			else{
wolffd@0: 				subs[i]++;
wolffd@0: 				if(sy[i])
wolffd@0: 					sTable += cpsy[i];
wolffd@0: 				break;
wolffd@0: 			}
wolffd@0: 		}
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "T");
wolffd@0: 	if(pTemp)mxDestroyArray(pTemp);
wolffd@0: 	pTemp = convert_to_sparse(bTable, NB, nzCounts);
wolffd@0: 	mxSetField(bigPot, 0, "T", pTemp);
wolffd@0: 	pTemp1 = mxGetField(bigPot, 0, "T");
wolffd@0: 	jc = mxGetJc(pTemp1);
wolffd@0: 
wolffd@0: 	free(sx);
wolffd@0: 	free(sy);
wolffd@0: 	free(s);
wolffd@0: 	free(cpsy);
wolffd@0: 	free(subs);
wolffd@0: 	free(cpsy2);
wolffd@0:     free(mask);
wolffd@0: 	free(bTable);
wolffd@0: }
wolffd@0: 
wolffd@0: void multiply_null_by_spPot(mxArray *bigPot, const mxArray *smallPot){
wolffd@0: 	int     i, j, count, count1, match, temp, bdim, sdim, diffdim, NB, NS, ND, NZB, NZS, bindex, sindex, nzCounts=0;
wolffd@0: 	int     *samemask, *diffmask, *sir, *sjc, *bCumprod, *sCumprod, *ssubv, *sequence, *weight;
wolffd@0: 	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *spr;
wolffd@0: 	mxArray *pTemp, *pTemp1;
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "domain");
wolffd@0: 	pbDomain = mxGetPr(pTemp);
wolffd@0: 	bdim = mxGetNumberOfElements(pTemp);
wolffd@0: 	pTemp = mxGetField(smallPot, 0, "domain");
wolffd@0: 	psDomain = mxGetPr(pTemp);
wolffd@0: 	sdim = mxGetNumberOfElements(pTemp);
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "sizes");
wolffd@0: 	pbSize = mxGetPr(pTemp);
wolffd@0: 	pTemp = mxGetField(smallPot, 0, "sizes");
wolffd@0: 	psSize = mxGetPr(pTemp);
wolffd@0: 
wolffd@0: 	NB = 1;
wolffd@0: 	for(i=0; i<bdim; i++){
wolffd@0: 		NB *= (int)pbSize[i];
wolffd@0: 	}
wolffd@0: 	NS = 1;
wolffd@0: 	for(i=0; i<sdim; i++){
wolffd@0: 		NS *= (int)psSize[i];
wolffd@0: 	}
wolffd@0: 	ND = NB / NS;
wolffd@0: 
wolffd@0: 	if(ND == 1){
wolffd@0: 		pTemp = mxGetField(bigPot, 0, "T");
wolffd@0: 		if(pTemp)mxDestroyArray(pTemp);
wolffd@0: 		pTemp1 = mxGetField(smallPot, 0, "T");
wolffd@0: 		pTemp = mxDuplicateArray(pTemp1);
wolffd@0: 		mxSetField(bigPot, 0, "T", pTemp);
wolffd@0: 		return;
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(smallPot, 0, "T");
wolffd@0: 	spr = mxGetPr(pTemp);
wolffd@0: 	sir = mxGetIr(pTemp);
wolffd@0: 	sjc = mxGetJc(pTemp);
wolffd@0: 	NZS = sjc[1];
wolffd@0: 
wolffd@0: 	NZB = ND * NZS;
wolffd@0: 
wolffd@0: 	diffdim = bdim - sdim;
wolffd@0: 	sequence = malloc(NZB * 2 * sizeof(int));
wolffd@0: 	bigTable = malloc(NZB * sizeof(double));
wolffd@0: 	samemask = malloc(sdim * sizeof(int));
wolffd@0: 	diffmask = malloc(diffdim * sizeof(int));
wolffd@0: 	bCumprod = malloc(bdim * sizeof(int));
wolffd@0: 	sCumprod = malloc(sdim * sizeof(int));
wolffd@0: 	weight = malloc(ND * sizeof(int));
wolffd@0: 	ssubv = malloc(sdim * sizeof(int));
wolffd@0: 
wolffd@0: 	count = 0;
wolffd@0: 	count1 = 0;
wolffd@0: 	for(i=0; i<bdim; i++){
wolffd@0: 		match = 0;
wolffd@0: 		for(j=0; j<sdim; j++){
wolffd@0: 			if(pbDomain[i] == psDomain[j]){
wolffd@0: 				samemask[count] = i;
wolffd@0: 				match = 1;
wolffd@0: 				count++;
wolffd@0: 				break;
wolffd@0: 			}
wolffd@0: 		}
wolffd@0: 		if(match == 0){
wolffd@0: 			diffmask[count1] = i; 
wolffd@0: 			count1++;
wolffd@0: 		}
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	bCumprod[0] = 1;
wolffd@0: 	for(i=0; i<bdim-1; i++){
wolffd@0: 		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
wolffd@0: 	}
wolffd@0: 	sCumprod[0] = 1;
wolffd@0: 	for(i=0; i<sdim-1; i++){
wolffd@0: 		sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	count = 0;
wolffd@0: 	compute_fixed_weight(weight, pbSize, diffmask, bCumprod, ND, diffdim);
wolffd@0: 	for(i=0; i<NZS; i++){
wolffd@0: 		sindex = sir[i];
wolffd@0: 		ind_subv(sindex, sCumprod, sdim, ssubv);
wolffd@0: 		temp = 0;
wolffd@0: 		for(j=0; j<sdim; j++){
wolffd@0: 			temp += ssubv[j] * bCumprod[samemask[j]];
wolffd@0: 		}
wolffd@0: 		for(j=0; j<ND; j++){
wolffd@0: 			bindex = weight[j] + temp;
wolffd@0: 			bigTable[nzCounts] = spr[i];
wolffd@0: 			sequence[count] = bindex;
wolffd@0: 			count++;
wolffd@0: 			sequence[count] = nzCounts;
wolffd@0: 			nzCounts++;
wolffd@0: 			count++;
wolffd@0: 		}
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "T");
wolffd@0: 	if(pTemp)mxDestroyArray(pTemp);
wolffd@0: 	qsort(sequence, nzCounts, sizeof(int) * 2, compare);
wolffd@0: 	pTemp = convert_ill_table_to_sparse(bigTable, sequence, nzCounts, NB);
wolffd@0: 	mxSetField(bigPot, 0, "T", pTemp);
wolffd@0: 
wolffd@0: 	free(sequence); 
wolffd@0: 	free(bigTable);
wolffd@0: 	free(samemask);
wolffd@0: 	free(diffmask);
wolffd@0: 	free(bCumprod);
wolffd@0: 	free(sCumprod);
wolffd@0: 	free(weight);
wolffd@0: 	free(ssubv);
wolffd@0: }
wolffd@0: 
wolffd@0: void multiply_spPot_by_fuPot(mxArray *bigPot, const mxArray *smallPot){
wolffd@0: 	int     i, j, count, bdim, sdim, NB, NZB, bindex, sindex, nzCounts=0;
wolffd@0: 	int     *mask, *index, *bir, *bjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
wolffd@0: 	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value;
wolffd@0: 	mxArray *pTemp;
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "domain");
wolffd@0: 	pbDomain = mxGetPr(pTemp);
wolffd@0: 	bdim = mxGetNumberOfElements(pTemp);
wolffd@0: 	pTemp = mxGetField(smallPot, 0, "domain");
wolffd@0: 	psDomain = mxGetPr(pTemp);
wolffd@0: 	sdim = mxGetNumberOfElements(pTemp);
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "sizes");
wolffd@0: 	pbSize = mxGetPr(pTemp);
wolffd@0: 	pTemp = mxGetField(smallPot, 0, "sizes");
wolffd@0: 	psSize = mxGetPr(pTemp);
wolffd@0: 
wolffd@0: 	NB = 1;
wolffd@0: 	for(i=0; i<bdim; i++){
wolffd@0: 		NB *= (int)pbSize[i];
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "T");
wolffd@0: 	bpr = mxGetPr(pTemp);
wolffd@0: 	bir = mxGetIr(pTemp);
wolffd@0: 	bjc = mxGetJc(pTemp);
wolffd@0: 	NZB = bjc[1];
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(smallPot, 0, "T");
wolffd@0: 	spr = mxGetPr(pTemp);
wolffd@0: 
wolffd@0: 	bigTable = malloc(NZB * sizeof(double));
wolffd@0: 	index = malloc(NZB * sizeof(double));
wolffd@0: 	mask = malloc(sdim * sizeof(int));
wolffd@0: 	bCumprod = malloc(bdim * sizeof(int));
wolffd@0: 	sCumprod = malloc(sdim * sizeof(int));
wolffd@0: 	bsubv = malloc(bdim * sizeof(int));
wolffd@0: 	ssubv = malloc(sdim * sizeof(int));
wolffd@0: 
wolffd@0: 	for(i=0; i<NZB; i++){
wolffd@0: 		bigTable[i] = 0;
wolffd@0: 	}
wolffd@0: 	count = 0;
wolffd@0: 	for(i=0; i<sdim; i++){
wolffd@0: 		for(j=0; j<bdim; j++){
wolffd@0: 			if(psDomain[i] == pbDomain[j]){
wolffd@0: 				mask[count] = j;
wolffd@0: 				count++;
wolffd@0: 				break;
wolffd@0: 			}
wolffd@0: 		}
wolffd@0: 	}
wolffd@0: 	
wolffd@0: 	bCumprod[0] = 1;
wolffd@0: 	for(i=0; i<bdim-1; i++){
wolffd@0: 		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
wolffd@0: 	}
wolffd@0: 	sCumprod[0] = 1;
wolffd@0: 	for(i=0; i<sdim-1; i++){
wolffd@0: 		sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	for(i=0; i<NZB; i++){
wolffd@0: 		bindex = bir[i];
wolffd@0: 		ind_subv(bindex, bCumprod, bdim, bsubv);
wolffd@0: 		for(j=0; j<sdim; j++){
wolffd@0: 			ssubv[j] = bsubv[mask[j]];
wolffd@0: 		}
wolffd@0: 		sindex = subv_ind(sdim, sCumprod, ssubv);
wolffd@0: 		value = spr[sindex];
wolffd@0: 		if(value != 0){
wolffd@0: 			bigTable[nzCounts] = bpr[i] * value;
wolffd@0: 			index[nzCounts] = bindex;
wolffd@0: 			nzCounts++;
wolffd@0: 		}
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "T");
wolffd@0: 	if(pTemp)mxDestroyArray(pTemp);
wolffd@0: 	pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB);
wolffd@0: 	mxSetField(bigPot, 0, "T", pTemp);
wolffd@0: 
wolffd@0: 	free(bigTable);
wolffd@0: 	free(index);
wolffd@0: 	free(mask);
wolffd@0: 	free(bCumprod);
wolffd@0: 	free(sCumprod);
wolffd@0: 	free(bsubv);
wolffd@0: 	free(ssubv);
wolffd@0: }
wolffd@0: 
wolffd@0: void multiply_spPot_by_spPot(mxArray *bigPot, const mxArray *smallPot){
wolffd@0: 	int     i, j, count, bdim, sdim, NB, NZB, NZS, position, bindex, sindex, nzCounts=0;
wolffd@0: 	int     *mask, *index, *result, *bir, *sir, *bjc, *sjc, *bCumprod, *sCumprod, *bsubv, *ssubv;
wolffd@0: 	double  *bigTable, *pbDomain, *psDomain, *pbSize, *psSize, *bpr, *spr, value;
wolffd@0: 	mxArray *pTemp;
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "domain");
wolffd@0: 	pbDomain = mxGetPr(pTemp);
wolffd@0: 	bdim = mxGetNumberOfElements(pTemp);
wolffd@0: 	pTemp = mxGetField(smallPot, 0, "domain");
wolffd@0: 	psDomain = mxGetPr(pTemp);
wolffd@0: 	sdim = mxGetNumberOfElements(pTemp);
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "sizes");
wolffd@0: 	pbSize = mxGetPr(pTemp);
wolffd@0: 	pTemp = mxGetField(smallPot, 0, "sizes");
wolffd@0: 	psSize = mxGetPr(pTemp);
wolffd@0: 
wolffd@0: 	NB = 1;
wolffd@0: 	for(i=0; i<bdim; i++){
wolffd@0: 		NB *= (int)pbSize[i];
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "T");
wolffd@0: 	bpr = mxGetPr(pTemp);
wolffd@0: 	bir = mxGetIr(pTemp);
wolffd@0: 	bjc = mxGetJc(pTemp);
wolffd@0: 	NZB = bjc[1];
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(smallPot, 0, "T");
wolffd@0: 	spr = mxGetPr(pTemp);
wolffd@0: 	sir = mxGetIr(pTemp);
wolffd@0: 	sjc = mxGetJc(pTemp);
wolffd@0: 	NZS = sjc[1];
wolffd@0: 
wolffd@0: 	bigTable = malloc(NZB * sizeof(double));
wolffd@0: 	index = malloc(NZB * sizeof(double));
wolffd@0: 	mask = malloc(sdim * sizeof(int));
wolffd@0: 	bCumprod = malloc(bdim * sizeof(int));
wolffd@0: 	sCumprod = malloc(sdim * sizeof(int));
wolffd@0: 	bsubv = malloc(bdim * sizeof(int));
wolffd@0: 	ssubv = malloc(sdim * sizeof(int));
wolffd@0: 
wolffd@0: 	for(i=0; i<NZB; i++){
wolffd@0: 		bigTable[i] = 0;
wolffd@0: 	}
wolffd@0: 	count = 0;
wolffd@0: 	for(i=0; i<sdim; i++){
wolffd@0: 		for(j=0; j<bdim; j++){
wolffd@0: 			if(psDomain[i] == pbDomain[j]){
wolffd@0: 				mask[count] = j;
wolffd@0: 				count++;
wolffd@0: 				break;
wolffd@0: 			}
wolffd@0: 		}
wolffd@0: 	}
wolffd@0: 	
wolffd@0: 	bCumprod[0] = 1;
wolffd@0: 	for(i=0; i<bdim-1; i++){
wolffd@0: 		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
wolffd@0: 	}
wolffd@0: 	sCumprod[0] = 1;
wolffd@0: 	for(i=0; i<sdim-1; i++){
wolffd@0: 		sCumprod[i+1] = sCumprod[i] * (int)psSize[i];
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	for(i=0; i<NZB; i++){
wolffd@0: 		value = bpr[i];
wolffd@0: 		bindex = bir[i];
wolffd@0: 		ind_subv(bindex, bCumprod, bdim, bsubv);
wolffd@0: 		for(j=0; j<sdim; j++){
wolffd@0: 			ssubv[j] = bsubv[mask[j]];
wolffd@0: 		}
wolffd@0: 		sindex = subv_ind(sdim, sCumprod, ssubv);
wolffd@0: 		result = (int *) bsearch(&sindex, sir, NZS, sizeof(int), compare);
wolffd@0: 		if(result){
wolffd@0: 			position = result - sir;
wolffd@0: 			value *= spr[position];
wolffd@0: 			bigTable[nzCounts] = value;
wolffd@0: 			index[nzCounts] = bindex;
wolffd@0: 			nzCounts++;
wolffd@0: 		}
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	pTemp = mxGetField(bigPot, 0, "T");
wolffd@0: 	if(pTemp)mxDestroyArray(pTemp);
wolffd@0: 	pTemp = convert_table_to_sparse(bigTable, index, nzCounts, NB);
wolffd@0: 	mxSetField(bigPot, 0, "T", pTemp);
wolffd@0: 
wolffd@0: 	free(bigTable);
wolffd@0: 	free(index);
wolffd@0: 	free(mask);
wolffd@0: 	free(bCumprod);
wolffd@0: 	free(sCumprod);
wolffd@0: 	free(bsubv);
wolffd@0: 	free(ssubv);
wolffd@0: }
wolffd@0: 
wolffd@0: 
wolffd@0: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
wolffd@0: 	int     i, j, c, loop, nNodes, nCliques, ndomain, dims[2];
wolffd@0: 	double  *pClqs, *pr, *pt, *pSize;
wolffd@0: 	mxArray *pTemp, *pTemp1, *pStruct, *pCliques, *pBigpot, *pSmallpot;
wolffd@0: 	const char *field_names[] = {"domain", "T", "sizes"};
wolffd@0: 
wolffd@0: 	nNodes = mxGetNumberOfElements(prhs[1]);
wolffd@0: 	pCliques = mxGetField(prhs[0], 0, "cliques");
wolffd@0: 	nCliques = mxGetNumberOfElements(pCliques);
wolffd@0: 	pTemp = mxGetField(prhs[0], 0, "eff_node_sizes");
wolffd@0: 	pSize = mxGetPr(pTemp);
wolffd@0: 
wolffd@0: 	plhs[0] = mxCreateCellArray(1, &nCliques);
wolffd@0:     for(i=0; i<nCliques; i++){
wolffd@0:         pStruct = mxCreateStructMatrix(1, 1, 3, field_names);
wolffd@0: 		mxSetCell(plhs[0], i, pStruct);
wolffd@0: 		pTemp = mxGetCell(pCliques, i);
wolffd@0: 		ndomain = mxGetNumberOfElements(pTemp);
wolffd@0: 		pt = mxGetPr(pTemp);
wolffd@0: 		pTemp1 = mxDuplicateArray(pTemp);
wolffd@0: 		mxSetField(pStruct, 0, "domain", pTemp1);
wolffd@0: 		
wolffd@0: 		pTemp = mxCreateDoubleMatrix(1, ndomain, mxREAL);
wolffd@0: 		mxSetField(pStruct, 0, "sizes", pTemp);
wolffd@0: 		pr = mxGetPr(pTemp);
wolffd@0:         for(j=0; j<ndomain; j++){
wolffd@0:             pr[j] = pSize[(int)pt[j]-1];
wolffd@0:         }
wolffd@0:     }
wolffd@0: 
wolffd@0: 	pClqs = mxGetPr(prhs[1]);
wolffd@0: 	for(loop=0; loop<nNodes; loop++){
wolffd@0: 		c = (int)pClqs[loop] - 1;
wolffd@0: 		pSmallpot = mxGetCell(prhs[2], loop);
wolffd@0: 		pTemp = mxGetField(pSmallpot, 0, "T");
wolffd@0: 		pBigpot = mxGetCell(plhs[0], c);
wolffd@0: 		pTemp1 = mxGetField(pBigpot, 0, "T");
wolffd@0: 		if(pTemp1){
wolffd@0: 			if(mxIsSparse(pTemp))
wolffd@0: 				multiply_spPot_by_spPot(pBigpot, pSmallpot);
wolffd@0: 			else multiply_spPot_by_fuPot(pBigpot, pSmallpot);
wolffd@0: 		}
wolffd@0: 		else{
wolffd@0: 			if(mxIsSparse(pTemp))
wolffd@0: 				multiply_null_by_spPot(pBigpot, pSmallpot);
wolffd@0: 			else multiply_null_by_fuPot(pBigpot, pSmallpot);
wolffd@0: 		}		
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	dims[0] = nCliques;
wolffd@0: 	dims[1] = nCliques;
wolffd@0: 	plhs[1] = mxCreateCellArray(2, dims);
wolffd@0: }
wolffd@0: 
wolffd@0: