Daniel@0: /* convert_to_sparse_table.c  convert a sparse discrete CPD with evidence into sparse table */
Daniel@0: /* convert_to_pot.m located in ../CPDs/discrete_CPD call it */
Daniel@0: /* 3 input */
Daniel@0: /* CPD      prhs[0] with 1D sparse CPT */
Daniel@0: /* domain   prhs[1]                    */
Daniel@0: /* evidence prhs[2]                    */
Daniel@0: /* 1 output */
Daniel@0: /* T        plhs[0] sparse table       */
Daniel@0: 
Daniel@0: #include <math.h>
Daniel@0: #include "mex.h"
Daniel@0: 
Daniel@0: void ind_subv(int index, const int *cumprod, const int n, int *bsubv){
Daniel@0: 	int i;
Daniel@0: 
Daniel@0: 	for (i = n-1; i >= 0; i--) {
Daniel@0: 		bsubv[i] = ((int)floor(index / cumprod[i]));
Daniel@0: 		index = index % cumprod[i];
Daniel@0: 	}
Daniel@0: }
Daniel@0: 
Daniel@0: int subv_ind(const int n, const int *cumprod, const int *subv){
Daniel@0: 	int i, index=0;
Daniel@0: 
Daniel@0: 	for(i=0; i<n; i++){
Daniel@0: 		index += subv[i] * cumprod[i];
Daniel@0: 	}
Daniel@0: 	return index;
Daniel@0: }
Daniel@0: 
Daniel@0: void reset_nzmax(mxArray *spArray, const int old_nzmax, const int new_nzmax){
Daniel@0: 	double *ptr;
Daniel@0: 	void   *newptr;
Daniel@0: 	int    *ir, *jc;
Daniel@0: 	int    nbytes;
Daniel@0: 
Daniel@0: 	if(new_nzmax == old_nzmax) return;
Daniel@0: 	nbytes = new_nzmax * sizeof(*ptr);
Daniel@0: 	ptr = mxGetPr(spArray);
Daniel@0: 	newptr = mxRealloc(ptr, nbytes);
Daniel@0: 	mxSetPr(spArray, newptr);
Daniel@0: 	nbytes = new_nzmax * sizeof(*ir);
Daniel@0: 	ir = mxGetIr(spArray);
Daniel@0: 	newptr = mxRealloc(ir, nbytes);
Daniel@0: 	mxSetIr(spArray, newptr);
Daniel@0: 	jc = mxGetJc(spArray);
Daniel@0: 	jc[0] = 0;
Daniel@0: 	jc[1] = new_nzmax;
Daniel@0: 	mxSetNzmax(spArray, new_nzmax);
Daniel@0: }
Daniel@0: 
Daniel@0: 
Daniel@0: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
Daniel@0: 	int     i, j, NS, NZB, count, bdim, match, domain, bindex, sindex, nzCounts=0;
Daniel@0: 	int     *observed, *bsubv, *ssubv, *bir, *sir, *bjc, *sjc, *mask, *ssize, *bcumprod, *scumprod;
Daniel@0: 	double  *pDomain, *pSize, *bpr, *spr;
Daniel@0: 	mxArray *pTemp;
Daniel@0: 
Daniel@0: 	pTemp = mxGetField(prhs[0], 0, "CPT");
Daniel@0: 	bpr = mxGetPr(pTemp);
Daniel@0: 	bir = mxGetIr(pTemp);
Daniel@0: 	bjc = mxGetJc(pTemp);
Daniel@0: 	NZB = bjc[1];
Daniel@0: 	pTemp = mxGetField(prhs[0], 0, "sizes");
Daniel@0: 	pSize = mxGetPr(pTemp);
Daniel@0: 
Daniel@0: 	pDomain = mxGetPr(prhs[1]);
Daniel@0: 	bdim = mxGetNumberOfElements(prhs[1]);
Daniel@0: 
Daniel@0: 	mask = malloc(bdim * sizeof(int));
Daniel@0: 	ssize = malloc(bdim * sizeof(int));
Daniel@0: 	observed = malloc(bdim * sizeof(int));
Daniel@0: 
Daniel@0: 	for(i=0; i<bdim; i++){
Daniel@0: 		ssize[i] = (int)pSize[i];
Daniel@0: 	}
Daniel@0: 
Daniel@0: 	count = 0;
Daniel@0: 	for(i=0; i<bdim; i++){
Daniel@0: 		domain = (int)pDomain[i] - 1;
Daniel@0: 		pTemp = mxGetCell(prhs[2], domain);
Daniel@0: 		if(pTemp){
Daniel@0: 			mask[count] = i;
Daniel@0: 			ssize[i] = 1;
Daniel@0: 			observed[count] = (int)mxGetScalar(pTemp) - 1;
Daniel@0: 			count++;
Daniel@0: 		}
Daniel@0: 	}
Daniel@0: 
Daniel@0: 	if(count == 0){
Daniel@0: 		pTemp = mxGetField(prhs[0], 0, "CPT");
Daniel@0: 		plhs[0] = mxDuplicateArray(pTemp);
Daniel@0: 		free(mask);
Daniel@0: 		free(ssize);
Daniel@0: 		free(observed);
Daniel@0: 		return;
Daniel@0: 	}
Daniel@0: 
Daniel@0: 	bsubv = malloc(bdim * sizeof(int));
Daniel@0: 	ssubv = malloc(count * sizeof(int));
Daniel@0: 	bcumprod = malloc(bdim * sizeof(int));
Daniel@0: 	scumprod = malloc(bdim * sizeof(int));
Daniel@0: 
Daniel@0: 	NS = 1;
Daniel@0: 	for(i=0; i<bdim; i++){
Daniel@0: 		NS *= ssize[i];
Daniel@0: 	}
Daniel@0: 
Daniel@0: 	plhs[0] = mxCreateSparse(NS, 1, NS, mxREAL);
Daniel@0: 	spr = mxGetPr(plhs[0]);
Daniel@0: 	sir = mxGetIr(plhs[0]);
Daniel@0: 	sjc = mxGetJc(plhs[0]);
Daniel@0: 	sjc[0] = 0;
Daniel@0: 	sjc[1] = NS;
Daniel@0: 
Daniel@0: 	bcumprod[0] = 1;
Daniel@0: 	scumprod[0] = 1;
Daniel@0: 	for(i=0; i<bdim-1; i++){
Daniel@0: 		bcumprod[i+1] = bcumprod[i] * (int)pSize[i];
Daniel@0: 		scumprod[i+1] = scumprod[i] * ssize[i];
Daniel@0: 	}
Daniel@0: 
Daniel@0: 	nzCounts = 0;
Daniel@0: 	for(i=0; i<NZB; i++){
Daniel@0: 		bindex = bir[i];
Daniel@0: 		ind_subv(bindex, bcumprod, bdim, bsubv);
Daniel@0: 		for(j=0; j<count; j++){
Daniel@0: 			ssubv[j] = bsubv[mask[j]];
Daniel@0: 		}
Daniel@0: 		match = 1;
Daniel@0: 		for(j=0; j<count; j++){
Daniel@0: 			if((ssubv[j]) != observed[j]){
Daniel@0: 				match = 0;
Daniel@0: 				break;
Daniel@0: 			}
Daniel@0: 		}
Daniel@0: 		if(match){
Daniel@0: 			spr[nzCounts] = bpr[i];
Daniel@0: 			sindex = subv_ind(bdim, scumprod, bsubv);
Daniel@0: 			sir[nzCounts] = sindex;
Daniel@0: 			nzCounts++;
Daniel@0: 		}
Daniel@0: 	}
Daniel@0: 
Daniel@0: 	reset_nzmax(plhs[0], NS, nzCounts);
Daniel@0: 	free(mask);
Daniel@0: 	free(ssize);
Daniel@0: 	free(observed);
Daniel@0: 	free(bsubv);
Daniel@0: 	free(ssubv);
Daniel@0: 	free(bcumprod);
Daniel@0: 	free(scumprod);
Daniel@0: }
Daniel@0: