matthiasm@8: /* C mex version of ind2subv.m in misc directory */
matthiasm@8: /* 2 input, 1 output         */
matthiasm@8: /* siz, ndx                  */
matthiasm@8: /* sub                       */
matthiasm@8: 
matthiasm@8: #include "mex.h"
matthiasm@8: #include <math.h>
matthiasm@8: 
matthiasm@8: void rbinary(int num, int n, double *rbits){
matthiasm@8: 	int i, mask;
matthiasm@8: 	num = num - 1;
matthiasm@8: 
matthiasm@8: 	mask = 1 << (n-1); /* mask = 00100...0 , where the 1 is in column n (rightmost = col 1) */
matthiasm@8: 	for (i = 0; i < n; i++) {
matthiasm@8: 		rbits[n-i-1] = ((num & mask) == 0) ? 1 : 2;
matthiasm@8: 		num <<= 1;
matthiasm@8: 	}
matthiasm@8: }
matthiasm@8: 
matthiasm@8: void ind_subv(int num, const double *sizes, int n, double *rbits){
matthiasm@8: 	int i;
matthiasm@8: 	int *cumprod;
matthiasm@8: 
matthiasm@8: 	cumprod  = malloc(n * sizeof(int));
matthiasm@8: 	num = num - 1;
matthiasm@8: 	cumprod[0] = 1;
matthiasm@8: 	for (i = 0; i < n-1; i++)
matthiasm@8: 	cumprod[i+1] = cumprod[i] * (int)sizes[i];
matthiasm@8: 	for (i = n-1; i >= 0; i--) {
matthiasm@8: 		rbits[i] = ((int)floor(num / cumprod[i])) + 1;
matthiasm@8: 		num = num % cumprod[i];
matthiasm@8: 	}
matthiasm@8: 	free(cumprod);
matthiasm@8: }
matthiasm@8: 
matthiasm@8: 
matthiasm@8: void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
matthiasm@8: 	int    i, j, k, nCol, nRow, nnRow, binary, count, temp, temp1, start;
matthiasm@8: 	double *pSize, *pNdx, *pr;
matthiasm@8: 	double ndx;
matthiasm@8: 	int    *subv, *cumprod, *templai;
matthiasm@8: 
matthiasm@8: 	pSize = mxGetPr(prhs[0]);
matthiasm@8: 	pNdx = mxGetPr(prhs[1]);
matthiasm@8: 	nCol = mxGetNumberOfElements(prhs[0]);
matthiasm@8: 	nnRow = mxGetNumberOfElements(prhs[1]);
matthiasm@8: 
matthiasm@8: 	nRow = 1;
matthiasm@8: 	for(i=0; i<nCol; i++){
matthiasm@8: 		nRow *= (int)pSize[i];
matthiasm@8: 	}
matthiasm@8: 
matthiasm@8: 	if(nCol == 0){
matthiasm@8: 		plhs[0] = mxDuplicateArray(prhs[1]);
matthiasm@8: 		return;
matthiasm@8: 	}
matthiasm@8: 
matthiasm@8: 	binary = 2;
matthiasm@8: 	for (i = 0; i < nCol; i++){
matthiasm@8: 		if (pSize[i] > 2.0){
matthiasm@8: 			binary = 0;
matthiasm@8: 			break;
matthiasm@8: 		}
matthiasm@8: 		else if((int)pSize[i] == 1){
matthiasm@8: 			binary = 1;
matthiasm@8: 		}
matthiasm@8: 	}
matthiasm@8: 
matthiasm@8: 	if(nnRow == 1){
matthiasm@8: 		ndx = mxGetScalar(prhs[1]);
matthiasm@8: 		plhs[0] = mxCreateDoubleMatrix(1, nCol, mxREAL);
matthiasm@8: 		pr = mxGetPr(plhs[0]);
matthiasm@8: 		if(binary == 2)rbinary((int)ndx, nCol, pr);
matthiasm@8: 		else ind_subv((int)ndx, pSize, nCol, pr);
matthiasm@8: 		return;
matthiasm@8: 	}
matthiasm@8: 
matthiasm@8: 	plhs[0] = mxCreateDoubleMatrix(nnRow, nCol, mxREAL);
matthiasm@8: 	pr = mxGetPr(plhs[0]);
matthiasm@8: 
matthiasm@8: 	subv = malloc(nRow * nCol * sizeof(int));
matthiasm@8: 
matthiasm@8: 	if (binary == 2) {
matthiasm@8: 		for(j=0; j<nCol; j++){
matthiasm@8: 			temp = (1 << j);
matthiasm@8: 			temp1 = j * nRow;
matthiasm@8: 			for(i=0; i<nRow; i++){
matthiasm@8: 				subv[temp1 + i] = ((i & temp) == 0) ? 1 : 2;
matthiasm@8: 			}
matthiasm@8: 		}
matthiasm@8: 	}
matthiasm@8: 	else if(binary == 1){
matthiasm@8: 		cumprod = (int *)malloc(nCol * sizeof(int));
matthiasm@8: 		templai = (int *)malloc(nCol * sizeof(int));
matthiasm@8: 		cumprod[0] = 1;
matthiasm@8: 		templai[0] = nRow;
matthiasm@8: 		for(i=1; i<nCol; i++){
matthiasm@8: 			k = (int)pSize[i-1] - 1;
matthiasm@8: 			cumprod[i] = cumprod[i-1] << k;
matthiasm@8: 			templai[i] = templai[i-1] >> k;
matthiasm@8: 		}
matthiasm@8: 		for(j=0; j<nCol; j++){
matthiasm@8: 			temp1 = j * nRow;
matthiasm@8: 			if(pSize[j] == 1.0){
matthiasm@8: 				for(i=0; i<nRow; i++){
matthiasm@8: 					subv[temp1 + i] = 1;
matthiasm@8: 				}
matthiasm@8: 			}
matthiasm@8: 			else{
matthiasm@8: 				temp = 1;
matthiasm@8: 				count = 0;
matthiasm@8: 				for(i=0; i<templai[j]; i++){
matthiasm@8: 					if(temp > 2) temp = 1;
matthiasm@8: 					for(k=0; k<cumprod[j]; k++){
matthiasm@8: 						subv[temp1 + count] = temp;
matthiasm@8:                         count++;
matthiasm@8: 					}
matthiasm@8: 					temp++;
matthiasm@8: 				}
matthiasm@8: 			}
matthiasm@8: 		}
matthiasm@8: 		free(templai);
matthiasm@8: 		free(cumprod);
matthiasm@8: 	}
matthiasm@8: 	else {
matthiasm@8: 		cumprod = (int *)malloc(nCol * sizeof(int));
matthiasm@8: 		templai = (int *)malloc(nCol * sizeof(int));
matthiasm@8: 		cumprod[0] = 1;
matthiasm@8: 		for(i=1; i<nCol; i++){
matthiasm@8: 			cumprod[i] = (int)(cumprod[i-1] * pSize[i-1]);
matthiasm@8: 		}
matthiasm@8: 		templai[0] = nRow;
matthiasm@8: 		for(i=1; i<nCol; i++){
matthiasm@8: 			templai[i] = nRow / cumprod[i];
matthiasm@8: 		}
matthiasm@8: 		for(j=0; j<nCol; j++){
matthiasm@8: 			temp1 = j * nRow;
matthiasm@8: 			if(pSize[j] == 1.0){
matthiasm@8: 				for(i=0; i<nRow; i++){
matthiasm@8: 					subv[temp1 + i] = 1;
matthiasm@8: 				}
matthiasm@8: 			}
matthiasm@8: 			else{
matthiasm@8: 				temp = 1;
matthiasm@8: 				count = 0;
matthiasm@8: 				for(i=0; i<templai[j]; i++){
matthiasm@8: 					if(temp > (int)pSize[j]) temp = 1;
matthiasm@8: 					for(k=0; k<cumprod[j]; k++){
matthiasm@8: 						subv[temp1 + count] = temp;
matthiasm@8: 						count++;
matthiasm@8: 					}
matthiasm@8: 					temp++;
matthiasm@8: 				}
matthiasm@8: 			}
matthiasm@8: 		}
matthiasm@8: 		free(cumprod);
matthiasm@8: 		free(templai);
matthiasm@8: 	}
matthiasm@8: 
matthiasm@8: 	count = 0;
matthiasm@8: 	for(j=0; j<nCol; j++){
matthiasm@8: 		temp1 = j * nRow;
matthiasm@8: 		for(i=0; i<nnRow; i++){
matthiasm@8: 			start = (int)pNdx[i] - 1;
matthiasm@8: 			pr[count] = subv[temp1 + start];
matthiasm@8: 			count++;
matthiasm@8: 		}
matthiasm@8: 	}
matthiasm@8: 	free(subv);
matthiasm@8: }
matthiasm@8: 
matthiasm@8: 
matthiasm@8: 
matthiasm@8: 
matthiasm@8: 
matthiasm@8: