Daniel@0: /* marg_sparse_table.c ../potential/tables*/
Daniel@0: 
Daniel@0: /******************************************/
Daniel@0: /* 5 input & 1 output                     */
Daniel@0: /* Big sparse table                       */
Daniel@0: /* Big domain                             */
Daniel@0: /* Big sizes                              */
Daniel@0: /* onto                                   */
Daniel@0: /* maximize, if missed, maximize=0        */
Daniel@0: /*                                        */
Daniel@0: /* small sparse table                     */
Daniel@0: /******************************************/
Daniel@0: 
Daniel@0: #include <math.h>
Daniel@0: #include <stdlib.h>
Daniel@0: #include "mex.h"
Daniel@0: 
Daniel@0: int compare(const void* src1, const void* src2){
Daniel@0: 	int i1 = *(int*)src1 ;
Daniel@0: 	int i2 = *(int*)src2 ;
Daniel@0: 	return i1-i2 ;
Daniel@0: }
Daniel@0: 
Daniel@0: void ind_subv(int index, const int *cumprod, 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: mxArray* convert_table_to_sparse(const double *Table, const int *sequence, const int nzCounts, const int N){
Daniel@0: 	mxArray *spTable;
Daniel@0: 	int     i, temp, *irs, *jcs, count=0;
Daniel@0: 	double  *sr;
Daniel@0: 
Daniel@0: 	spTable = mxCreateSparse(N, 1, nzCounts, mxREAL);
Daniel@0:     sr  = mxGetPr(spTable);
Daniel@0:     irs = mxGetIr(spTable);
Daniel@0:     jcs = mxGetJc(spTable);
Daniel@0: 
Daniel@0: 	jcs[0] = 0;
Daniel@0: 	jcs[1] = nzCounts;
Daniel@0: 
Daniel@0: 	for(i=0; i<nzCounts; i++){
Daniel@0: 		irs[i] = sequence[count];
Daniel@0: 		count++;
Daniel@0: 		temp = sequence[count];
Daniel@0: 		sr[i] = Table[temp];
Daniel@0: 		count++;
Daniel@0: 	}
Daniel@0: 	return spTable;
Daniel@0: }
Daniel@0: 
Daniel@0: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
Daniel@0: 	int        i, j, count, bdim, sdim, NS, NZB, position, bindex, sindex, maximize, nzCounts=0;
Daniel@0: 	int        *mask, *sequence, *result, *bir, *bjc, *ssize, *bCumprod, *sCumprod, *bsubv, *ssubv;
Daniel@0: 	double     *sTable, *pbDomain, *psDomain, *pbSize, *bpr, *spr;
Daniel@0: 	const char *field_names[] = {"domain", "T", "sizes"};
Daniel@0: 
Daniel@0: 	if(nrhs < 5) maximize = 0;
Daniel@0: 	else maximize = (int)mxGetScalar(prhs[4]);
Daniel@0: 
Daniel@0: 	bdim = mxGetNumberOfElements(prhs[1]);
Daniel@0: 	sdim = mxGetNumberOfElements(prhs[3]);
Daniel@0: 	pbSize = mxGetPr(prhs[2]);
Daniel@0: 	pbDomain = mxGetPr(prhs[1]);
Daniel@0: 	psDomain = mxGetPr(prhs[3]);
Daniel@0: 	bpr = mxGetPr(prhs[0]);
Daniel@0: 	bir = mxGetIr(prhs[0]);
Daniel@0: 	bjc = mxGetJc(prhs[0]);
Daniel@0: 	NZB = bjc[1];
Daniel@0: 
Daniel@0: 	if(sdim == 0){
Daniel@0: 		plhs[0] = mxCreateSparse(1, 1, 1, mxREAL);
Daniel@0: 		spr = mxGetPr(plhs[0]);
Daniel@0: 		bir = mxGetIr(plhs[0]);
Daniel@0: 		bjc = mxGetJc(plhs[0]);
Daniel@0: 		*spr = 0;
Daniel@0: 		*bir = 0;
Daniel@0: 		bjc[0] = 0;
Daniel@0: 		bjc[1] = 1;
Daniel@0: 		if(maximize){
Daniel@0: 			for(i=0; i<NZB; i++){
Daniel@0: 				*spr = (*spr < bpr[i])? bpr[i] : *spr;
Daniel@0: 			}
Daniel@0: 		}
Daniel@0: 		else{
Daniel@0: 			for(i=0; i<NZB; i++){
Daniel@0: 				*spr += bpr[i];
Daniel@0: 			}
Daniel@0: 		}	
Daniel@0: 		return;
Daniel@0: 	}
Daniel@0: 
Daniel@0: 	mask = malloc(sdim * sizeof(int));
Daniel@0: 	count = 0;
Daniel@0: 	for(i=0; i<sdim; i++){
Daniel@0: 		for(j=0; j<bdim; j++){
Daniel@0: 			if(psDomain[i] == pbDomain[j]){
Daniel@0: 				mask[count] = j;
Daniel@0: 				count++;
Daniel@0: 				break;
Daniel@0: 			}
Daniel@0: 		}
Daniel@0: 	}
Daniel@0: 		
Daniel@0: 	sTable = malloc(NZB * sizeof(double));
Daniel@0: 	sequence = malloc(NZB * 2 * sizeof(double));
Daniel@0: 	bCumprod = malloc(bdim * sizeof(int));
Daniel@0: 	sCumprod = malloc(sdim * sizeof(int));
Daniel@0: 	bsubv = malloc(bdim * sizeof(int));
Daniel@0: 	ssubv = malloc(sdim * sizeof(int));
Daniel@0: 	ssize = malloc(sdim * sizeof(int));
Daniel@0: 
Daniel@0: 	NS = 1;
Daniel@0: 	for(i=0; i<count; i++){
Daniel@0: 		ssize[i] = (int)pbSize[mask[i]];
Daniel@0: 		NS *= ssize[i];
Daniel@0: 	}
Daniel@0: 
Daniel@0: 	for(i=0; i<NZB; i++)sTable[i] = 0;
Daniel@0: 	
Daniel@0: 	bCumprod[0] = 1;
Daniel@0: 	for(i=0; i<bdim-1; i++){
Daniel@0: 		bCumprod[i+1] = bCumprod[i] * (int)pbSize[i];
Daniel@0: 	}
Daniel@0: 	sCumprod[0] = 1;
Daniel@0: 	for(i=0; i<sdim-1; i++){
Daniel@0: 		sCumprod[i+1] = sCumprod[i] * ssize[i];
Daniel@0: 	}
Daniel@0: 
Daniel@0: 	count = 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<sdim; j++){
Daniel@0: 			ssubv[j] = bsubv[mask[j]];
Daniel@0: 		}
Daniel@0: 		sindex = subv_ind(sdim, sCumprod, ssubv);
Daniel@0: 		result = (int *) bsearch(&sindex, sequence, nzCounts, sizeof(int)*2, compare);
Daniel@0: 		if(result){
Daniel@0: 			position = (result - sequence) / 2;
Daniel@0: 			if(maximize) 
Daniel@0: 				sTable[position] = (sTable[position] < bpr[i]) ? bpr[i] : sTable[position];
Daniel@0: 			else sTable[position] += bpr[i];
Daniel@0: 		}
Daniel@0: 		else {
Daniel@0: 			if(maximize) 
Daniel@0: 				sTable[nzCounts] = (sTable[nzCounts] < bpr[i]) ? bpr[i] : sTable[nzCounts];
Daniel@0: 			else sTable[nzCounts] += bpr[i];
Daniel@0: 			sequence[count] = sindex;
Daniel@0: 			count++;
Daniel@0: 			sequence[count] = nzCounts;
Daniel@0: 			nzCounts++;
Daniel@0: 			count++;
Daniel@0: 		}
Daniel@0: 	}
Daniel@0: 	
Daniel@0: 	qsort(sequence, nzCounts, sizeof(int) * 2, compare);
Daniel@0: 	plhs[0] = convert_table_to_sparse(sTable, sequence, nzCounts, NS);
Daniel@0: 
Daniel@0: 	free(sTable);
Daniel@0: 	free(sequence);
Daniel@0: 	free(mask);
Daniel@0: 	free(bCumprod);
Daniel@0: 	free(sCumprod);
Daniel@0: 	free(bsubv);
Daniel@0: 	free(ssubv);
Daniel@0: 	free(ssize);
Daniel@0: }