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