From d12bc05d86faaee9217b4c94030ec4b01b534eab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 18 Aug 2010 15:04:33 +0000 Subject: [PATCH] Add MEX support for inverting cell-to-block mappings (i.e., partition vectors) to create block-to-cell mappings. --- mex_partition_invert.c | 160 +++++++++++++++++++++++++++++++++++++++++ mex_partition_invert.m | 51 +++++++++++++ partition.c | 107 +++++++++++++++++++++++++++ partition.h | 16 +++++ 4 files changed, 334 insertions(+) create mode 100644 mex_partition_invert.c create mode 100644 mex_partition_invert.m diff --git a/mex_partition_invert.c b/mex_partition_invert.c new file mode 100644 index 00000000..1296b3f8 --- /dev/null +++ b/mex_partition_invert.c @@ -0,0 +1,160 @@ +#include +#include + +#include + +#include "partition.h" + + +#define MAX(a,b) (((a) > (b)) ? (a) : (b)) + + +/* ---------------------------------------------------------------------- */ +static int +args_ok(int nlhs, int nrhs, const mxArray *prhs[]) +/* ---------------------------------------------------------------------- */ +{ + int ok; + + ok = (nlhs == 2) || (nlhs == 3); + ok = ok && (nrhs == 1); + ok = ok && (mxIsDouble(prhs[0]) || mxIsInt32(prhs[0])); + + return ok; +} + + +/* ---------------------------------------------------------------------- */ +static void +extract_rhs(const mxArray *M_p, int *p) +/* ---------------------------------------------------------------------- */ +{ + size_t e, ne; + + int *pi; + double *pd; + + ne = mxGetNumberOfElements(M_p); + + if (mxIsDouble(M_p)) { + pd = mxGetPr(M_p); + + for (e = 0; e < ne; e++) { p[e] = pd[e] - 1; } + } else { + pi = mxGetData(M_p); + + for (e = 0; e < ne; e++) { p[e] = pi[e] - 1; } + } +} + + +/* ---------------------------------------------------------------------- */ +static void +assign_int_vec(const int *v, mxArray *M_v) +/* ---------------------------------------------------------------------- */ +{ + size_t e, ne; + + int *pi; + double *pd; + + ne = mxGetNumberOfElements(M_v); + + if (mxIsDouble(M_v)) { + pd = mxGetPr(M_v); + + for (e = 0; e < ne; e++) { pd[e] = v[e] + 1; } + } else { + pi = mxGetData(M_v); + + for (e = 0; e < ne; e++) { pi[e] = v[e] + 1; } + } +} + + +/* ---------------------------------------------------------------------- */ +static int +max_block(int nc, const int *p) +/* ---------------------------------------------------------------------- */ +{ + int c, m; + + m = -1; + for (c = 0; c < nc; c++) + { + m = MAX(m, p[c]); + } + + return m; +} + + +/* ---------------------------------------------------------------------- */ +static void +adjust_one_based_idx(size_t n, int *v) +/* ---------------------------------------------------------------------- */ +{ + size_t i; + + for (i = 0; i < n; i++) { v[i] += 1; } +} + + +/* + * [pb2c, b2c] = mex_partition_invert(p) + * [pb2c, b2c, loc] = mex_partition_invert(p) + */ + +/* ---------------------------------------------------------------------- */ +void +mexFunction(int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) +/* ---------------------------------------------------------------------- */ +{ + int nc, max_blk, *p, *pb2c, *b2c, *loc; + char errmsg[1023 + 1]; + + if (args_ok(nlhs, nrhs, prhs)) { + nc = mxGetNumberOfElements(prhs[0]); + p = mxMalloc(nc * sizeof *p); + + extract_rhs(prhs[0], p); + max_blk = max_block(nc, p); + + if (partition_allocate_inverse(nc, max_blk, &pb2c, &b2c)) { + plhs[0] = mxCreateNumericMatrix(max_blk + 1 + 1, 1, mxINT32_CLASS, mxREAL); + plhs[1] = mxCreateNumericMatrix(nc, 1, mxINT32_CLASS, mxREAL); + if (nlhs == 3) { + plhs[2] = mxCreateNumericMatrix(nc, 1, mxINT32_CLASS, mxREAL); + } + + partition_invert(nc, p, pb2c, b2c); + + assign_int_vec(pb2c, plhs[0]); + assign_int_vec(b2c , plhs[1]); + + if (nlhs == 3) { + partition_localidx(max_blk + 1, pb2c, b2c, + mxGetData(plhs[2])); + adjust_one_based_idx(nc, mxGetData(plhs[2])); + } + } else { + plhs[0] = mxCreateDoubleScalar(mxGetNaN()); + plhs[1] = mxCreateDoubleScalar(mxGetNaN()); + if (nlhs == 3) { + plhs[2] = mxCreateDoubleScalar(mxGetNaN()); + } + } + + partition_deallocate_inverse(pb2c, b2c); + mxFree(p); /* p != NULL guaranteed here */ + } else { + sprintf(errmsg, + "Calling sequence is\n" + "\t[pb2c, b2c] = %s(p) %% or\n" + "\t[pb2c, b2c, loc] = %s(p)", + mexFunctionName(), mexFunctionName()); + + mexErrMsgTxt(errmsg); + } +} diff --git a/mex_partition_invert.m b/mex_partition_invert.m new file mode 100644 index 00000000..7d98a887 --- /dev/null +++ b/mex_partition_invert.m @@ -0,0 +1,51 @@ +function varargout = mex_partition_invert(varargin) +%Invert cell-to-block map (creating block-to-cell) using compiled C code. +% +% SYNOPSIS: +% [pb2c, b2c] = mex_partition_invert(p) +% [pb2c, b2c, loc] = mex_partition_invert(p) +% +% PARAMETERS: +% p - Partition vector. Should not contain any empty blocks. Use +% function 'mex_partition_compress' to remove empty blocks/bins. +% +% RETURNS: +% pb2c - Indirection map of size [MAX(p) + 1, 1] into the 'b2c' map +% array. Specifically, the cells of block 'b' are stored in +% +% b2c(pb2c(b) : pb2c(b + 1) - 1) +% +% b2c - Inverse cell map. The entries in pb2c(b):pb2c(b+1)-1 correspond +% to the result of FIND(p == b). +% +% loc - Local index within a block/bin. Specifically, +% +% loc(i) == FIND(b2c == i) - pb2c(p(i)) + 1 +% +% OPTIONAL. Only returned (and computed) if specifically +% requested. +% +% SEE ALSO: +% mex_partition_ui, mex_partition_compress. + +%{ +#COPYRIGHT# +%} + +% $Date$ +% $Revision$ + + buildmex CFLAGS="\$CFLAGS -Wall -Wextra -ansi -pedantic ... + -Wformat-nonliteral -Wcast-align -Wpointer-arith ... + -Wbad-function-cast -Wmissing-prototypes -Wstrict-prototypes ... + -Wmissing-declarations -Winline -Wundef -Wnested-externs ... + -Wcast-qual -Wshadow -Wconversion -Wwrite-strings ... + -Wno-conversion -Wchar-subscripts -Wredundant-decls" ... + ... + -O -largeArrayDims ... + ... + mex_partition_invert.c partition.c + + % Call MEX'ed edition. + [varargout{1:nargout}] = mex_partition_invert(varargin{:}); +end diff --git a/partition.c b/partition.c index 2381fcf3..954beecb 100644 --- a/partition.c +++ b/partition.c @@ -129,3 +129,110 @@ partition_compress(int n, int *p) return ret; } + + +/* ---------------------------------------------------------------------- */ +void +partition_deallocate_inverse(int *pi, int *inverse) +/* ---------------------------------------------------------------------- */ +{ + free(inverse); + free(pi); +} + + +/* ---------------------------------------------------------------------- */ +int +partition_allocate_inverse(int nc, int max_bin, + int **pi, int **inverse) +/* ---------------------------------------------------------------------- */ +{ + int nbin, ret, *ptr, *i; + + nbin = max_bin + 1; + + ptr = malloc((nbin + 1) * sizeof *ptr); + i = malloc(nc * sizeof *i ); + + if ((ptr == NULL) || (i == NULL)) { + partition_deallocate_inverse(ptr, i); + + *pi = NULL; + *inverse = NULL; + + ret = 0; + } else { + *pi = ptr; + *inverse = i; + + ret = nc; + } + + return ret; +} + + +/* ---------------------------------------------------------------------- */ +void +partition_invert(int nc, const *p, int *pi, int *inverse) +/* ---------------------------------------------------------------------- */ +{ + int nbin, b, i, j, tmp; + + nbin = 0; + for (i = 0; i < nc; i++) { + nbin = MAX(nbin, p[i]); + } + nbin += 1; /* Adjust for bin 0 */ + + /* Zero start pointers */ + for (b = 0; b < nbin; b++) { pi[b] = 0; } + + /* Count elements per bin */ + for (i = 0; i < nc ; i++) { pi[ p[i] ]++; } + + /* Derive start pointers for b=1:nbin */ + for (b = 1; b < nbin; b++) { pi[b] += pi[b - 1]; } + + /* Set end pointer in last bin */ + assert (pi[nbin - 1] == nc); + pi[nbin] = nc; + + /* Reverse insert bin elements whilst deriving start pointers */ + for (i = 0; i < nc; i++) { + inverse[-- pi[ p[i] ]] = i; + } + assert (pi[0] == 0); + + /* Reverse the reverse order, creating final inverse mapping */ + for (b = 0; b < nbin; b++) { + i = pi[b + 0] + 0; + j = pi[b + 1] - 1; + + while (i < j) { + /* Swap reverse (lower <-> upper) */ + tmp = inverse[i]; + inverse[i] = inverse[j]; + inverse[j] = tmp; + + i += 1; /* Increase lower bound */ + j -= 1; /* Decrease upper bound */ + } + } +} + + +/* ---------------------------------------------------------------------- */ +void +partition_localidx(int nbin, const int *pi, const int *inverse, + int *localidx) +/* ---------------------------------------------------------------------- */ +{ + int b, i; + + for (b = 0; b < nbin; b++) { + for (i = pi[b]; i < pi[b + 1]; i++) { + localidx[ inverse[i] ] = i - pi[b]; + } + } +} diff --git a/partition.h b/partition.h index bd5ceb69..783a1b35 100644 --- a/partition.h +++ b/partition.h @@ -11,4 +11,20 @@ partition_unif_idx(int ndims, int nc, int partition_compress(int n, int *p); + +int +partition_allocate_inverse(int nc, int max_blk, + int **pi, int **inverse); + +void +partition_deallocate_inverse(int *pi, int *inverse); + +void +partition_invert(int nc, const int *p, + int *pi, int *inverse); + +void +partition_localidx(int nblk, const int *pi, const int *inverse, + int *localidx); + #endif /* PARTITION_H_INLCUDED */