Add MEX support for inverting cell-to-block mappings (i.e.,
partition vectors) to create block-to-cell mappings.
This commit is contained in:
Bård Skaflestad
parent
a98073072f
commit
d12bc05d86
160
mex_partition_invert.c
Normal file
160
mex_partition_invert.c
Normal file
@ -0,0 +1,160 @@
|
|||||||
|
#include <stddef.h>
|
||||||
|
#include <string.h>
|
||||||
|
|
||||||
|
#include <mex.h>
|
||||||
|
|
||||||
|
#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);
|
||||||
|
}
|
||||||
|
}
|
||||||
51
mex_partition_invert.m
Normal file
51
mex_partition_invert.m
Normal file
@ -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
|
||||||
107
partition.c
107
partition.c
@ -129,3 +129,110 @@ partition_compress(int n, int *p)
|
|||||||
|
|
||||||
return ret;
|
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];
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|||||||
16
partition.h
16
partition.h
@ -11,4 +11,20 @@ partition_unif_idx(int ndims, int nc,
|
|||||||
int
|
int
|
||||||
partition_compress(int n, int *p);
|
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 */
|
#endif /* PARTITION_H_INLCUDED */
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user