wolffd@0: /* marg_table.c  ../potential/tables     */
wolffd@0: 
wolffd@0: 
wolffd@0: /******************************************/
wolffd@0: /* 5 input & 1 output                     */
wolffd@0: /* Big table                              */
wolffd@0: /* Big domain                             */
wolffd@0: /* Big sizes                              */
wolffd@0: /* onto                                   */
wolffd@0: /* maximize, if missed, maximize=0        */
wolffd@0: /*                                        */
wolffd@0: /* small table                            */
wolffd@0: /******************************************/
wolffd@0: 
wolffd@0: #include "mex.h"
wolffd@0: 
wolffd@0: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
wolffd@0: 	int     i, j, count, NB, NS, siz_b, siz_s, ndim, temp, maximize;
wolffd@0: 	int     *mask, *sx, *sy, *cpsy, *subs, *s, *cpsy2, *ssize;
wolffd@0: 	double  *pb, *ps, *bp, *sp, *pbd;
wolffd@0: 
wolffd@0: 
wolffd@0: 	siz_b = mxGetNumberOfElements(prhs[1]);
wolffd@0: 	siz_s = mxGetNumberOfElements(prhs[3]);
wolffd@0: 	pb = mxGetPr(prhs[1]);
wolffd@0: 	ps = mxGetPr(prhs[3]);
wolffd@0: 
wolffd@0: 	NB = mxGetNumberOfElements(prhs[0]);
wolffd@0: 	bp = mxGetPr(prhs[0]);
wolffd@0: 
wolffd@0: 	pbd = mxGetPr(prhs[2]);
wolffd@0: 
wolffd@0: 	if(nrhs < 5) maximize = 0;
wolffd@0: 	else maximize = (int)mxGetScalar(prhs[4]);
wolffd@0: 
wolffd@0: 	if(siz_s == 0){
wolffd@0: 		plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
wolffd@0: 		sp = mxGetPr(plhs[0]);
wolffd@0: 		if(maximize){
wolffd@0: 			for(i=0; i<NB; i++){
wolffd@0: 				*sp = (*sp < bp[i])? bp[i] : *sp;
wolffd@0: 			}
wolffd@0: 		}
wolffd@0: 		else{
wolffd@0: 			for(i=0; i<NB; i++){
wolffd@0: 				*sp += bp[i];
wolffd@0: 			}
wolffd@0: 		}
wolffd@0: 		return;
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	mask = malloc(siz_s * sizeof(int));
wolffd@0: 	ssize = 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(ps[i] == pb[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)pbd[i];
wolffd@0: 		sy[i] = 1;
wolffd@0: 	}
wolffd@0: 	for(i=0; i<siz_s; i++){
wolffd@0: 		temp = mask[i];
wolffd@0: 		sy[temp] = sx[temp];
wolffd@0: 		ssize[i] = sx[temp];
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	NS = 1;
wolffd@0: 	for(i=0; i<ndim; i++){
wolffd@0: 		NS *= sy[i];
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	plhs[0] = mxCreateNumericArray(siz_s, ssize, mxDOUBLE_CLASS, mxREAL);
wolffd@0: 	sp = mxGetPr(plhs[0]);
wolffd@0: 
wolffd@0: 	if(NS == 1){
wolffd@0: 		if(maximize){
wolffd@0: 			for(i=0; i<NB; i++){
wolffd@0: 				*sp = (*sp < bp[i])? bp[i] : *sp;
wolffd@0: 			}
wolffd@0: 		}
wolffd@0: 		else{
wolffd@0: 			for(i=0; i<NB; i++){
wolffd@0: 				*sp += bp[i];
wolffd@0: 			}
wolffd@0: 		}
wolffd@0: 		free(mask);
wolffd@0: 		free(sx);
wolffd@0: 		free(sy);
wolffd@0: 		free(ssize);
wolffd@0: 		return;
wolffd@0: 	}
wolffd@0: 
wolffd@0: 	if(NS == NB){
wolffd@0: 		for(i=0; i<NB; i++) *sp++ = *bp++;
wolffd@0: 		free(mask);
wolffd@0: 		free(sx);
wolffd@0: 		free(sy);
wolffd@0: 		free(ssize);
wolffd@0: 		return;
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: 	if(maximize){
wolffd@0: 		for(j=0; j<NB; j++){
wolffd@0: 			*sp = (*sp < *bp)? *bp : *sp;
wolffd@0: 			bp++;
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: 						sp -= cpsy2[i];
wolffd@0: 				}
wolffd@0: 				else{
wolffd@0: 					subs[i]++;
wolffd@0: 					if(sy[i])
wolffd@0: 						sp += cpsy[i];
wolffd@0: 					break;
wolffd@0: 				}
wolffd@0: 			}
wolffd@0: 		}
wolffd@0: 	}
wolffd@0: 	else{
wolffd@0: 		for(j=0; j<NB; j++){
wolffd@0: 			*sp += *bp++;
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: 						sp -= cpsy2[i];
wolffd@0: 				}
wolffd@0: 				else{
wolffd@0: 					subs[i]++;
wolffd@0: 					if(sy[i])
wolffd@0: 						sp += cpsy[i];
wolffd@0: 					break;
wolffd@0: 				}
wolffd@0: 			}
wolffd@0: 		}
wolffd@0: 	}
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(ssize);
wolffd@0: }
wolffd@0: