Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
#include <mex.h>
#include <math.h>
#define INDEX_TYPE_ID mxUINT8_CLASS
#define INDEX_TYPE uint8_T
void get_fourier_coefs(const double * X, const INDEX_TYPE * inds,
const double * mins, const double * ptps,
int num_pts, int dim, int num_bases, int M,
double * pcs, double * tmp, int stride) {
/* Gets projection coefficients onto the Fourier orthonormal basis.
*
* The basis on [0, 1] is:
* phi_1(x) = 1
* phi_k(x) = \sqrt{2} cos( (k-1) pi x )
*
* The basis on [0, 1]^d is their tensor product.
*
* Rescales the data points by subtracting mins and dividing by ptps, and
* then pushing in to [0, 1] if it's outside that range.
*
* X: a column-major (num_pts x dim) array
* inds: a column-major (dim x num_bases) array
* mins: a minimum value to subtract for each dimension
* ptps: a peak-to-peak value to divide by for each dimension
* num_pts: number of points in x
* dim: the dimension
* num_bases: the number of basis functions
* M: the max index of a basis function
* pcs: num_bases space for output
* tmp: num_pts * dim * M space for temporary values
* stride: how far to space things along in pcs (1 is normal)
**/
// i ranges over num_pts; j over dim; k over M; l over num_bases
unsigned int i, j, k, l;
double Xij_pi;
double sqrt_2 = sqrt(2);
double this_pc, this_prod;
#define tmp_for(i, j, k) ( (k) + M * ((j) + dim * (i)) )
// apply each phi func to each dataset point, in tmp
for (i = 0; i < num_pts; i++) {
for (j = 0; j < dim; j++) {
Xij_pi = M_PI *
fmax(0, fmin((X[i + num_pts * j] - mins[j]) / ptps[j], 1));
tmp[tmp_for(i, j, 0)] = 1;
for (k = 1; k < M; k++) {
tmp[tmp_for(i, j, k)] = sqrt_2 * cos(Xij_pi * k);
}
}
}
// compute projection coefficients:
// mean over data points of products of tmp values
for (l = 0; l < num_bases; l++) {
this_pc = 0;
for (i = 0; i < num_pts; i++) {
this_prod = 1;
for (j = 0; j < dim; j++) {
k = inds[j + dim * l] - 1;
this_prod *= tmp[tmp_for(i, j, k)];
}
this_pc += this_prod / num_pts;
}
pcs[l * stride] = this_pc;
}
#undef tmp_for
}
void mexFunction(int nlhs, mxArray ** plhs, int nrhs, const mxArray ** prhs) {
// function [pc] = fourier_coefs_mex(Xs, M, inds, rescale_mins, rescale_ptps)
mwSize i;
int do_checks = nrhs != 6;
if (do_checks && nrhs != 5)
mexErrMsgTxt("fourier_coefs_mex takes exactly five args");
if (do_checks && nlhs != 1)
mexErrMsgTxt("fourier_coefs_mex returns exactly one output");
// inputs
const mxArray * X_arrays = prhs[0];
mwSize num_xs, dim, num_pts, max_num_pts = 0;
int is_cell = mxIsCell(X_arrays);
if (is_cell) {
num_xs = mxGetNumberOfElements(X_arrays);
if (num_xs == 0) {
plhs[0] = mxCreateDoubleMatrix(mxGetN(prhs[2]), 0, mxREAL);
return;
}
for (i = 0; i < num_xs; i++) {
num_pts = mxGetM(mxGetCell(X_arrays, i));
if (num_pts > max_num_pts) {
max_num_pts = num_pts;
}
}
dim = mxGetN(mxGetCell(X_arrays, 0));
} else {
num_xs = 1;
dim = mxGetN(X_arrays);
}
const mxArray * M_array = prhs[1];
if (do_checks && mxGetNumberOfElements(M_array) != 1)
mexErrMsgTxt("M should be a scalar");
int M = mxGetScalar(M_array);
const mxArray * inds_array = prhs[2];
if (do_checks && mxGetNumberOfDimensions(inds_array) != 2)
mexErrMsgTxt("inds should be 2d");
if (do_checks && mxGetM(inds_array) != dim)
mexErrMsgTxt("size(inds, 1) should be size(X, 2)");
int num_bases = mxGetN(inds_array);
if (do_checks && mxGetClassID(inds_array) != INDEX_TYPE_ID)
mexErrMsgTxt("X should be of type uint8");
const INDEX_TYPE * inds = (INDEX_TYPE *) mxGetData(inds_array);
const mxArray * mins_array = prhs[3];
if (do_checks && mxGetNumberOfDimensions(mins_array) != 2)
mexErrMsgTxt("rescale_mins should be 2d");
if (do_checks && mxGetM(mins_array) != 1)
mexErrMsgTxt("size(rescale_mins, 1) should be 1");
if (do_checks && mxGetN(mins_array) != dim)
mexErrMsgTxt("size(rescale_mins, 2) should be size(X, 2)");
const double * mins = mxGetPr(mins_array);
const mxArray * ptps_array = prhs[4];
if (do_checks && mxGetNumberOfDimensions(ptps_array) != 2)
mexErrMsgTxt("rescale_ptps should be 2d");
if (do_checks && mxGetM(ptps_array) != 1)
mexErrMsgTxt("size(rescale_ptps, 1) should be 1");
if (do_checks && mxGetN(ptps_array) != dim)
mexErrMsgTxt("size(rescale_ptps, 2) should be size(X, 2)");
const double * ptps = mxGetPr(ptps_array);
// allocate for outputs
mxArray * pcs_array = mxCreateDoubleMatrix(num_xs, num_bases, mxREAL);
plhs[0] = pcs_array;
double * pcs = mxGetPr(pcs_array);
// allocate temporary space
double * tmp = mxMalloc(sizeof(double) * max_num_pts * dim * M);
// do it!
mxArray * X_array;
for (i = 0; i < num_xs; i++) {
X_array = is_cell ? mxGetCell(X_arrays, i) : X_array;
if (do_checks) {
if (mxGetNumberOfDimensions(X_array) != 2)
mexErrMsgTxt("X should be 2d");
if (!mxIsDouble(X_array))
mexErrMsgTxt("X should be of type double");
if (mxGetN(X_array) != dim)
mexErrMsgTxt("X should be of consistent dimension");
}
num_pts = mxGetM(X_array);
get_fourier_coefs(
mxGetPr(X_array), inds, mins, ptps, num_pts, dim, num_bases, M,
pcs + i, tmp, num_xs);
}
mxFree(tmp);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.