Algorithm-Classifier-IsolationForest
view release on metacpan or search on metacpan
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
SV* scores_rv, SV* labels_rv){
STRLEN tl;
const double* sm;
AV* scores;
AV* labels;
int i;
if (!SvROK(scores_rv) || SvTYPE(SvRV(scores_rv)) != SVt_PVAV ||
!SvROK(labels_rv) || SvTYPE(SvRV(labels_rv)) != SVt_PVAV) {
croak("score_predict_split_xs: scores/labels must be arrayrefs");
}
sm = (const double*)SvPVbyte(sm_sv, tl);
scores = (AV*)SvRV(scores_rv);
labels = (AV*)SvRV(labels_rv);
av_clear(scores);
av_clear(labels);
if (n_pts > 0) {
av_extend(scores, n_pts - 1);
av_extend(labels, n_pts - 1);
}
for (i = 0; i < n_pts; i++) {
av_store(scores, i, newSVnv(exp(-sm[i] * inv)));
av_store(labels, i, newSViv(sm[i] <= sum_threshold ? 1 : 0));
}
}
/* Walk one point through one tree; returns the path length (depth plus
* the precomputed c(leaf size) adjustment from the leaf record).
*
* Invariant: every feature index stored in a tree node is in
* [0, n_feats). fit() builds trees against n_features columns and
* pack_input_xs writes exactly that many doubles per row, and
* _resolve_input rejects PackedData with a mismatched feature count.
* So the loop can omit per-iteration bounds checks on attr / fi --
* this is what lets the oblique dot product vectorize cleanly under
* the omp-simd reductions below. */
#if defined(__GNUC__) || defined(__clang__)
__attribute__((always_inline))
#endif
static inline double if_walk_tree(const double *nd, const int *ico,
const double *vco, const double *xi,
int n_feats) {
int ni = 0, depth = 0;
for (;;) {
const double *node = nd + (size_t)ni * IF_NZ;
int type = (int)node[0];
if (type == 0) {
/* node[2] is c(leaf size), precomputed by _pack_tree; a
* log() here would otherwise run once per point per tree. */
return depth + node[2];
}
if (type == 1) {
double fv = xi[(int)node[1]];
ni = (fv < node[2]) ? (int)node[3] : (int)node[4];
} else {
int coff = (int)node[1], nf = (int)node[2];
double b = node[5], dot = 0.0;
const double *val_p = vco + (size_t)coff;
/* Both children are known before the dot product resolves
* which one gets taken, so start pulling their records in
* now and let the FMA loop below hide the latency. One of
* the two prefetches is always wasted -- affordable here
* on the oblique path, where there is real work to hide it
* under, but not on the axis path, whose single compare
* resolves immediately. */
const int li = (int)node[3], ri = (int)node[4];
IF_PREFETCH(nd + (size_t)li * IF_NZ);
IF_PREFETCH(nd + (size_t)ri * IF_NZ);
if (nf == n_feats) {
/* Dense oblique split: this node uses every feature,
* so _pack_tree laid the coefficients out in feature
* order. No gather -- the inner loop is a textbook
* FMA-vectorizable dot product over two contiguous
* double streams. Common case in extended mode at
* the default extension_level (== n_feats-1). */
#ifdef _OPENMP
#pragma omp simd reduction(+:dot)
#endif
for (int k = 0; k < n_feats; k++) {
dot += val_p[k] * xi[k];
}
} else {
/* Sparse oblique split: only nf < n_feats features
* participate, so we still need the gather on
* xi[idx_p[k]]. Storing idx as contiguous int32
* (rather than interleaved doubles) keeps the gather
* pattern clean and the val[] load contiguous. */
const int *idx_p = ico + (size_t)coff;
#ifdef _OPENMP
#pragma omp simd reduction(+:dot)
#endif
for (int k = 0; k < nf; k++) {
dot += val_p[k] * xi[idx_p[k]];
}
}
ni = (dot <= b) ? li : ri;
}
depth++;
}
}
/* score_all_xs(nodes_av, idx_av, val_av, x_sv, sm_sv,
* n_pts, n_feats, n_trees, use_openmp)
*
* Scores all points across all trees in one C call. See header comment
* above for the bigger picture. Writes sm[i] = sum_over_trees(path_len);
* the caller need not zero-init sm.
*
* idx_av holds per-tree packed int32 buffers of feature indices and
* val_av holds per-tree packed double buffers of coefficients (the SoA
* counterpart of the old interleaved layout). See the file-top
* comment for the rationale.
*
* Thread-safety: the parallel region only reads node/idx/val/x pointers
* (extracted before the region) and writes sm[i] for a unique i per
* iteration. No Perl API is called from inside the parallel region. */
void score_all_xs(SV* nodes_av_sv, SV* idx_av_sv, SV* val_av_sv,
SV* x_sv, SV* sm_sv,
int n_pts, int n_feats, int n_trees,
int use_openmp){
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
}
return result;
}
void build_forest_xs(SV* x_sv, int n_pts, int n_feats, int n_trees,
int psi, int limit, int mode_flag, int ext_level,
SV* out_rv) {
dTHX;
STRLEN tl;
const double* x;
AV* out;
int* all;
int t, i;
if (!SvROK(out_rv) || SvTYPE(SvRV(out_rv)) != SVt_PVAV) {
croak("build_forest_xs: out must be an arrayref");
}
x = (const double*)SvPVbyte(x_sv, tl);
out = (AV*)SvRV(out_rv);
av_clear(out);
if (n_trees > 0) av_extend(out, n_trees - 1);
all = (int*)malloc(n_pts * sizeof(int));
for (t = 0; t < n_trees; t++) {
int* sample;
for (i = 0; i < n_pts; i++) all[i] = i;
for (i = 0; i < psi; i++) {
int j = i + (int)(Drand01() * (n_pts - i));
int tmp = all[i]; all[i] = all[j]; all[j] = tmp;
}
sample = (int*)malloc(psi * sizeof(int));
memcpy(sample, all, psi * sizeof(int));
av_store(out, t,
_build_node_c(aTHX_ x, n_feats, sample, psi, 0, limit,
mode_flag, ext_level));
}
free(all);
}
/* ---------------------------------------------------------------------
* build_forest_openmp_xs -- OpenMP-parallel fit() tree builder.
*
* build_forest_xs (above) is bit-identical to the pure-Perl path
* because every random draw goes through Drand01(), the same
* generator Perl's rand()/srand() use -- but that generator is a
* single mutable struct shared by the whole interpreter, so calling
* it concurrently from multiple OpenMP threads would be a data race.
* The same is true of any Perl API call (newAV, newSViv, ...): Perl's
* SV allocator isn't safe to call from multiple OS threads sharing one
* interpreter without a lock that would just serialise everything
* anyway.
*
* So this builder trades the bit-identical guarantee for real thread
* parallelism: each tree gets its own splitmix64 PRNG stream, seeded
* from a tree index (not thread id or scheduling order), so results
* are still reproducible for a fixed seed and n_trees regardless of
* OMP_NUM_THREADS -- just different from what build_forest_xs or the
* pure-Perl path would produce for the same seed. The one Drand01()
* call in this function happens before the parallel region starts
* (single-threaded), so the result still varies with the model's
* `seed` the way every other code path does; it isn't used inside the
* parallel loop.
*
* Each tree is built entirely with plain C data (row-index int arrays,
* a growable TreeBuf of packed doubles/ints) -- no Perl API call
* happens anywhere inside the parallel region. Each node record in
* TreeBuf uses _pack_tree's 6-double SoA layout (see the file-top
* comment), but the node ORDER differs: records are appended
* post-order (a node is pushed after both its children, since child
* indices must be known first), so the root is the last record --
* _pack_tree's pre-order puts it at 0. _unpack_forest accounts for
* this. Oblique coefficients are also always stored sparse (in the
* random pool's order) -- the dense-pack fast path is skipped because
* its only purpose is speeding up score_all_xs, and _rebuild_c_trees
* reapplies it anyway once the caller unpacks these buffers back into
* the standard Perl tree shape and re-derives the scoring buffers.
*
* After the parallel region, each tree's TreeBuf is copied into a Perl
* string SV (one memcpy each, serially) and stored into nodes_rv /
* idx_rv / val_rv -- the caller unpacks these into ordinary nested
* Perl trees for $self->{trees} (so to_json/persistence/_rebuild_c_trees
* are unaffected). ------------------------------------------------ */
typedef struct {
double *nodes; size_t n_nodes, cap_nodes;
int *idx; size_t n_idx, cap_idx;
double *val; size_t n_val, cap_val;
} TreeBuf;
static void tb_init(TreeBuf *b) {
b->nodes = NULL; b->n_nodes = 0; b->cap_nodes = 0;
b->idx = NULL; b->n_idx = 0; b->cap_idx = 0;
b->val = NULL; b->n_val = 0; b->cap_val = 0;
}
static void tb_free(TreeBuf *b) {
free(b->nodes); free(b->idx); free(b->val);
}
static int tb_push_node(TreeBuf *b, double f0, double f1, double f2,
double f3, double f4, double f5) {
double *slot;
if (b->n_nodes == b->cap_nodes) {
size_t newcap = b->cap_nodes ? b->cap_nodes * 2 : 64;
b->nodes = (double*)realloc(b->nodes, newcap * 6 * sizeof(double));
b->cap_nodes = newcap;
}
slot = b->nodes + b->n_nodes * 6;
slot[0] = f0; slot[1] = f1; slot[2] = f2;
slot[3] = f3; slot[4] = f4; slot[5] = f5;
return (int)(b->n_nodes++);
}
/* Appends n (idx[k], val[k]) pairs and returns the offset they start
* at -- the `coff` an oblique node record stores. */
static int tb_push_coef(TreeBuf *b, const int *idx, const double *val,
int n) {
int off = (int)b->n_idx;
if (b->n_idx + (size_t)n > b->cap_idx) {
size_t newcap = b->cap_idx ? b->cap_idx * 2 : 64;
if (newcap < b->n_idx + (size_t)n) newcap = b->n_idx + (size_t)n;
b->idx = (int*)realloc(b->idx, newcap * sizeof(int));
b->cap_idx = newcap;
}
if (b->n_val + (size_t)n > b->cap_val) {
size_t newcap = b->cap_val ? b->cap_val * 2 : 64;
if (newcap < b->n_val + (size_t)n) newcap = b->n_val + (size_t)n;
b->val = (double*)realloc(b->val, newcap * sizeof(double));
b->cap_val = newcap;
}
memcpy(b->idx + b->n_idx, idx, (size_t)n * sizeof(int));
memcpy(b->val + b->n_val, val, (size_t)n * sizeof(double));
b->n_idx += n;
b->n_val += n;
return off;
}
/* splitmix64 -- fast, well-mixed, and per-stream state fits in one
* uint64_t, which is all a thread-private PRNG needs here. Not
* cryptographic; doesn't need to be. */
static uint64_t sm64_next(uint64_t *s) {
uint64_t z = (*s += 0x9E3779B97F4A7C15ULL);
z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL;
z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL;
return z ^ (z >> 31);
}
static double sm64_drand(uint64_t *s) {
return (double)(sm64_next(s) >> 11) * (1.0 / 9007199254740992.0);
}
static double _ts_randn(uint64_t *s) {
double u1 = sm64_drand(s);
double u2;
if (u1 == 0.0) u1 = 1e-12;
u2 = sm64_drand(s);
return sqrt(-2.0 * log(u1)) * cos(6.283185307179586 * u2);
}
/* Thread-safe twin of _build_node_c: same split algorithm, but reads
* randomness from a thread-private splitmix64 stream instead of
* Drand01(), and writes into a TreeBuf instead of allocating Perl AVs
* -- so it touches no interpreter-global state and is safe to call
* concurrently from an OpenMP parallel region, one tree per thread. */
static int _build_node_packed(const double* x, int nf, int* idxs, int size,
int depth, int limit, int mode_flag,
int ext_active, TreeBuf *buf, uint64_t *rng) {
double *lo, *hi;
int *varying, nv, f, my_idx;
if (depth >= limit || size <= 1) {
my_idx = tb_push_node(buf, 0.0, (double)size, 0.0, 0.0, 0.0, 0.0);
free(idxs);
return my_idx;
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
for (k = 0; k < active; k++) {
int ff = pool[k];
double c = _ts_randn(rng);
double p = lo[ff] + sm64_drand(rng) * (hi[ff] - lo[ff]);
coef[k] = c;
b += c * p;
}
lidx = (int*)malloc(size * sizeof(int));
ridx = (int*)malloc(size * sizeof(int));
for (i = 0; i < size; i++) {
int row = idxs[i];
double dot = 0.0;
for (k = 0; k < active; k++) {
dot += coef[k] * x[(size_t)row * (size_t)nf + pool[k]];
}
if (dot <= b) lidx[ln++] = row; else ridx[rn++] = row;
}
free(idxs); free(lo); free(hi); free(varying);
li = _build_node_packed(x, nf, lidx, ln, depth + 1, limit,
mode_flag, ext_active, buf, rng);
ri = _build_node_packed(x, nf, ridx, rn, depth + 1, limit,
mode_flag, ext_active, buf, rng);
coff = tb_push_coef(buf, pool, coef, active);
my_idx = tb_push_node(buf, 2.0, (double)coff, (double)active,
(double)li, (double)ri, b);
free(pool); free(coef);
}
return my_idx;
}
void build_forest_openmp_xs(SV* x_sv, int n_pts, int n_feats, int n_trees,
int psi, int limit, int mode_flag,
int ext_level, SV* nodes_rv, SV* idx_rv,
SV* val_rv, int use_openmp) {
dTHX;
STRLEN tl;
const double* x;
AV *nodes_av, *idx_av, *val_av;
TreeBuf *bufs;
uint64_t base_seed;
int t;
if (!SvROK(nodes_rv) || SvTYPE(SvRV(nodes_rv)) != SVt_PVAV ||
!SvROK(idx_rv) || SvTYPE(SvRV(idx_rv)) != SVt_PVAV ||
!SvROK(val_rv) || SvTYPE(SvRV(val_rv)) != SVt_PVAV) {
croak("build_forest_openmp_xs: nodes/idx/val must be arrayrefs");
}
x = (const double*)SvPVbyte(x_sv, tl);
nodes_av = (AV*)SvRV(nodes_rv);
idx_av = (AV*)SvRV(idx_rv);
val_av = (AV*)SvRV(val_rv);
av_clear(nodes_av); av_clear(idx_av); av_clear(val_av);
if (n_trees > 0) {
av_extend(nodes_av, n_trees - 1);
av_extend(idx_av, n_trees - 1);
av_extend(val_av, n_trees - 1);
}
/* Single Drand01() call, before the parallel region starts, so it's
* still a plain serial call into the interpreter's RNG state. */
base_seed = (uint64_t)(Drand01() * 18446744073709551615.0);
bufs = (TreeBuf*)malloc((size_t)n_trees * sizeof(TreeBuf));
for (t = 0; t < n_trees; t++) tb_init(&bufs[t]);
#ifdef _OPENMP
#pragma omp parallel for schedule(dynamic) if(use_openmp)
#endif
for (int t = 0; t < n_trees; t++) {
/* Seeded from the tree index, not thread id or iteration order,
* so the mapping from tree -> RNG stream is independent of
* OMP_NUM_THREADS / scheduling. sm64_next() mixes once more so
* adjacent tree indices (which differ by one golden-ratio step)
* don't start from too-similar states. */
uint64_t rng = base_seed + (uint64_t)t * 0x9E3779B97F4A7C15ULL;
rng = sm64_next(&rng);
int *all = (int*)malloc((size_t)n_pts * sizeof(int));
int *sample;
int i;
for (i = 0; i < n_pts; i++) all[i] = i;
for (i = 0; i < psi; i++) {
int j = i + (int)(sm64_drand(&rng) * (n_pts - i));
int tmp = all[i]; all[i] = all[j]; all[j] = tmp;
}
sample = (int*)malloc((size_t)psi * sizeof(int));
memcpy(sample, all, (size_t)psi * sizeof(int));
free(all);
_build_node_packed(x, n_feats, sample, psi, 0, limit, mode_flag,
ext_level, &bufs[t], &rng);
}
for (t = 0; t < n_trees; t++) {
/* newSVpvn(NULL, 0) makes an undef SV, not an empty-string one --
* axis-mode trees never call tb_push_coef, so idx/val stay NULL.
* Pass "" instead so the Perl side's unpack('...', $sv) always
* gets a defined (if empty) string, never undef. */
av_store(nodes_av, t, newSVpvn((char*)bufs[t].nodes,
bufs[t].n_nodes * 6 * sizeof(double)));
av_store(idx_av, t, bufs[t].n_idx
? newSVpvn((char*)bufs[t].idx, bufs[t].n_idx * sizeof(int))
: newSVpvn("", 0));
av_store(val_av, t, bufs[t].n_val
? newSVpvn((char*)bufs[t].val, bufs[t].n_val * sizeof(double))
: newSVpvn("", 0));
tb_free(&bufs[t]);
}
free(bufs);
}
/* ---------------------------------------------------------------------
* impute_fill_xs(data_sv, n_pts, n_feats, how, out_rv)
*
* C replacement for _compute_impute_fill's Perl loop: walks the raw
* arrayref-of-arrayrefs directly (like pack_input_xs), collecting each
* feature's present (defined) values, then reduces them to one fill
* value per feature -- mean (how == 0) or median (how == 1) -- and
* writes n_feats doubles into out_rv.
*
* Values are collected in row order (i = 0..n_pts-1), the same order
* the Perl version's `grep { defined } map { $_->[$f] } @data` walks
* them in, so the mean's left-to-right summation lands on the exact
* same float as the Perl path -- use_c toggles speed here, not the
* computed fill, matching the rest of the module.
*
* The median is an exact order statistic (not summation-dependent), so
* it matches the Perl path's sort-based median by definition regardless
* of which selection algorithm finds it. Croaks with the same message
* as the Perl fallback if a feature has no present values anywhere in
* the dataset. */
typedef struct { double *v; size_t n, cap; } DVec;
static void dvec_push(DVec *d, double x) {
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
return a[lo];
}
/* Median of a[0..n-1] (reorders a[]). Odd n: the single middle order
* statistic. Even n: quickselect finds the lower-median at k = n/2-1,
* which leaves every a[i > k] >= a[k] (the standard quickselect
* post-condition) -- so the upper-median is just the min of that
* remaining slice, one more linear scan instead of a second full
* selection pass. */
static double _median_select(double *a, int n) {
if (n % 2 == 1) {
return _kth_smallest(a, n, n / 2);
} else {
int k = n / 2 - 1;
double lower = _kth_smallest(a, n, k);
double upper = a[k + 1];
int i;
for (i = k + 2; i < n; i++) {
if (a[i] < upper) upper = a[i];
}
return (lower + upper) / 2.0;
}
}
void impute_fill_xs(SV* data_sv, int n_pts, int n_feats, int how,
SV* out_rv) {
dTHX;
AV *outer, *out;
DVec *cols;
int i, f;
if (!SvROK(data_sv) || SvTYPE(SvRV(data_sv)) != SVt_PVAV) {
croak("impute_fill_xs: data must be an arrayref");
}
if (!SvROK(out_rv) || SvTYPE(SvRV(out_rv)) != SVt_PVAV) {
croak("impute_fill_xs: out must be an arrayref");
}
outer = (AV*)SvRV(data_sv);
out = (AV*)SvRV(out_rv);
cols = (DVec*)calloc((size_t)n_feats, sizeof(DVec));
for (i = 0; i < n_pts; i++) {
SV** row_pp = av_fetch(outer, i, 0);
AV* row;
if (!row_pp || !*row_pp || !SvROK(*row_pp) ||
SvTYPE(SvRV(*row_pp)) != SVt_PVAV) {
continue;
}
row = (AV*)SvRV(*row_pp);
for (f = 0; f < n_feats; f++) {
SV** v = av_fetch(row, f, 0);
if (v && *v && SvOK(*v)) {
dvec_push(&cols[f], SvNV(*v));
}
}
}
/* Validate every column before freeing anything: croak() longjmps
* out of this function, so any cleanup loop reachable after a
* partial computation has already started (and already freed some
* cols[i].v) risks a double free on those same pointers. Checking
* all columns up front, before the computation loop below frees
* anything, avoids that entirely. Matches the Perl fallback's
* behaviour of reporting the first empty column in feature order. */
for (f = 0; f < n_feats; f++) {
if (cols[f].n == 0) {
int col = f;
for (i = 0; i < n_feats; i++) free(cols[i].v);
free(cols);
croak("impute: feature column %d has no present values", col);
}
}
av_clear(out);
if (n_feats > 0) av_extend(out, n_feats - 1);
for (f = 0; f < n_feats; f++) {
double result;
if (how == 0) {
double sum = 0.0;
for (i = 0; i < (int)cols[f].n; i++) sum += cols[f].v[i];
result = sum / (double)cols[f].n;
} else {
result = _median_select(cols[f].v, (int)cols[f].n);
}
av_store(out, f, newSVnv(result));
free(cols[f].v);
}
free(cols);
}
/* ---------------------------------------------------------------------
* Online Isolation Forest (Algorithm::Classifier::IsolationForest::
* Online) learn / unlearn / score-row accelerators.
*
* Unlike everything above, these operate directly on LIVE Perl
* arrayref trees: online trees mutate on every learned point, so there
* is no immutable packed form to walk during learning. Node layout
* (Online.pm's _N_* constants):
*
* leaf: [0, count, \@lo, \@hi]
* internal: [1, count, \@lo, \@hi, attr, split, left, right]
*
* and a tree record is a hashref { root, count, depth_limit }; root is
* undef until the tree learns its first point, and a leaf built from an
* empty synthetic partition has undef lo/hi until a real point reaches
* it.
*
* Random draws go through Drand01() in EXACTLY the order the pure-Perl
* learn path calls rand(): an optional per-tree subsample gate (drawn
* only when subsample < 1), and -- only when a leaf splits --
* count * nf box-sample draws that SKIP zero-width features (the Perl
* _sample_box never draws for those), then per synthetic internal node
* one draw for the split feature and one for the split value, recursing
* left before right. A learn() with a given seed therefore produces
* BIT-IDENTICAL trees whether use_c is on or off (on nvsize == 8 perls;
* wide-NV perls keep extra low bits in the pure-Perl path, as with
* fit()), which is what lets the online class reuse the existing use_c
* knob for learning instead of growing a new one.
* ------------------------------------------------------------------ */
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
# validated against a fixed set of GCC/Clang -O flags rather than
# interpolated as-is, since this string eventually reaches a shell
# command line via ExtUtils::MakeMaker.
my $opt = $def_opt;
if ( defined $ENV{IF_OPT} ) {
if ( $ENV{IF_OPT} =~ /\A-O[0123sgz]\z/ ) {
$opt = $ENV{IF_OPT};
} else {
warn "Algorithm::Classifier::IsolationForest: ignoring invalid "
. "IF_OPT value '$ENV{IF_OPT}' (expected one of -O0 -O1 -O2 "
. "-O3 -Os -Og -Oz); using $opt\n";
}
}
# -march=<value> lets the compiler target specific instruction-set
# extensions (AVX2 gather + FMA, etc.) for the oblique dot product
# and the fit-time min/max scan's `#pragma omp simd` loops.
#
# IF_ARCH=<value> sets it explicitly (e.g. "x86-64-v3", "skylake",
# "znver3") -- validated against a conservative identifier charset
# since, like IF_OPT, it flows into a compiler command line.
# IF_NATIVE=1 remains as shorthand for IF_ARCH=native and is used
# when IF_ARCH isn't set. Prefer a specific IF_ARCH value over
# IF_NATIVE on a machine you don't control exclusively: blanket
# -march=native pulls in whatever the build host has, including
# AVX-512 on some Intel CPUs, which is known to trigger clock
# throttling under sustained heavy use and can make throughput
# *worse* than a conservative target like x86-64-v3 (AVX2, no
# AVX-512). Either way, the cached artefact under _Inline/ is then
# pinned to that instruction set, so leave both unset if the
# directory is shared across machines with different CPUs.
my $arch = $def_arch;
if ( defined $ENV{IF_ARCH} ) {
if ( $ENV{IF_ARCH} eq '' or $ENV{IF_ARCH} eq 'none' ) {
# Explicit opt-out: overrides an arch recorded at configure
# time (there is no other way to request a plain build on
# an install configured with IF_ARCH).
$arch = '';
} elsif ( $ENV{IF_ARCH} =~ /\A[A-Za-z0-9_.+=-]+\z/ ) {
$arch = $ENV{IF_ARCH};
} else {
warn "Algorithm::Classifier::IsolationForest: ignoring invalid " . "IF_ARCH value '$ENV{IF_ARCH}'\n";
}
} elsif ( $ENV{IF_NATIVE} ) {
$arch = 'native';
}
# -ffp-contract=off rides along with any -march: once the target
# has FMA (x86-64-v3, most -march=native hosts), the compiler may
# otherwise contract a*b+c expressions into fused multiply-adds
# whose different rounding breaks the documented guarantee that
# use_c => 1 and use_c => 0 build bit-identical trees (one ulp in a
# split value cascades into a structurally different tree). The
# -march speedup comes from AVX2 vectorization, not contraction,
# so this costs little (verified against the fit-determinism and
# scoring-parity tests).
my $opt_level = $opt;
$opt_level .= " -march=$arch -ffp-contract=off" if length $arch;
# IF_NO_OPENMP=1 forces the serial C build: the OpenMP compile attempt
# is skipped, so the object has no libgomp linkage and never starts an
# OpenMP runtime in the process. Distinct from OMP_NUM_THREADS=1,
# which runs the parallel code on a single thread but still loads
# libgomp. An explicit IF_NO_OPENMP=0 re-enables OpenMP over a
# no-openmp configure-time default.
my $no_omp
= defined $ENV{IF_NO_OPENMP}
? ( $ENV{IF_NO_OPENMP} ? 1 : 0 )
: $def_no_omp;
# The prebuilt object is only trusted when the effective flags match
# what it was compiled with; any difference -- or an explicit
# IF_RUNTIME_BUILD=1 -- falls through to the classic runtime Inline::C
# build below, which honours the requested flags via the MD5-keyed
# _Inline/ cache exactly as before prebuilt support existed.
# IF_INSTALL_BUILD is the `make` rule driving the install-time compile
# (see Makefile.PL); it must never short-circuit into loading an
# older object.
my $use_prebuilt
= $prebuilt
&& !$ENV{IF_RUNTIME_BUILD}
&& !$ENV{IF_INSTALL_BUILD}
&& $opt eq $def_opt
&& $arch eq $def_arch
&& $no_omp == $def_no_omp;
# Inline::C hashes the C source to decide whether to rebuild but
# does NOT include CCFLAGS / OPTIMIZE in that hash. Without the
# tag below, toggling IF_NATIVE/IF_ARCH/IF_OPT (or editing the
# optimisation flags here) would silently reuse a cached binary
# built with stale flags. Embedding the active flags as a leading
# comment forces the hash to differ when they change. The OpenMP
# and serial builds get distinct tags so they cache to separate
# artefacts.
my $omp_tag = "/* if_build: openmp $opt_level */\n";
my $serial_tag = "/* if_build: serial $opt_level */\n";
if ( $ENV{IF_INSTALL_BUILD} ) {
# `make` is driving: the rule Makefile.PL appended runs this load
# with IF_INSTALL_BUILD=1 and @ARGV = (version, INST_ARCHLIB),
# which is where Inline's install mode reads them from. _INSTALL_
# makes Inline compile the backend and place the shared object
# under blib/arch so `make install` ships it; NAME/VERSION give
# the object a fixed identity XSLoader can find at run time
# (Inline's install mode also requires both and checks VERSION
# against $ARGV[0]). Same OpenMP-then-serial fallback as the
# runtime build below.
my @install = (
NAME => __PACKAGE__,
VERSION => $VERSION,
_INSTALL_ => 1,
);
unless ($no_omp) {
local $@;
eval {
require Inline;
Inline->import(
C => $omp_tag . $C_CODE,
CCFLAGS => '-fopenmp',
OPTIMIZE => $opt_level,
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
=head1 NATIVE ACCELERATION (Inline::C and OpenMP)
Both the scoring hot path (C<score_samples>, C<predict>, C<path_lengths>,
C<score_predict_samples>, and C<score_predict_split>) and the C<fit()>
tree builder are automatically accelerated through
L<Inline::C> when it is installed and a working C compiler is reachable.
If the toolchain also accepts C<-fopenmp> and can link against
C<libgomp>, the per-point tree walk runs in parallel across all
available CPU cores using OpenMP, and the extended-mode oblique dot
product is vectorised via C<#pragma omp simd> -- which on modern x86
compilers translates to an unrolled FMA / AVX gather chain that's
substantially faster for high-feature-count extended models.
C<fit()>'s tree builder (subsampling plus the recursive axis/oblique
split search) runs in C the same way when C<use_c> is on, replacing the
per-node Perl arrayref copying with plain int-array partitioning --
typically an order of magnitude faster, and dramatically more so at
higher feature counts where the pure-Perl per-cell loop dominates. Its
random draws go through the same generator C<rand()>/C<srand()> use
internally, in the same call order the pure-Perl builder uses, so a
given C<seed> produces bit-identical trees whether C<use_c> is on or
off -- switching backends changes only how fast the model is built, not
the model itself. On perls whose NV is wider than a C double
(C<-Duselongdouble> / C<-Dusequadmath>) the pure-Perl builder rounds
each stored value to double precision to preserve this parity; axis
mode matches exactly, while extended mode can still differ on rare
libm rounding ties (double vs long-double transcendentals).
By default this C builder is single-threaded per call, because Perl's
RNG state isn't safe to share across OpenMP threads. Two ways to scale
fit() across cores are available (see below for why they don't compose):
=over 4
=item * C<parallel_fit> forks N worker processes, each building its
share of the trees with the (still single-threaded) C builder. Fixed
IPC/serialisation overhead per worker means this can cost more than it
saves once a fit already completes in milliseconds; it's most useful
once a single-process fit is large enough that the fork/Storable
overhead is small relative to the work being split.
=item * C<use_openmp_fit> builds trees across OpenMP threads within a
single process (one tree per thread), using a separate, thread-safe
PRNG seeded per tree index instead of Perl's C<rand()>. This means
trees built with C<use_openmp_fit> are I<not> bit-identical to the
default C<use_c> path for the same seed -- but a fixed seed and
C<n_trees> still reproduce the same trees regardless of
C<OMP_NUM_THREADS> or how OpenMP schedules the work. It's off by
default (unlike C<use_c>/C<use_openmp>, which only ever change speed,
this changes which trees get built) and only takes effect when C<use_c>
is also on and OpenMP is linked in.
=back
These two do NOT compose, despite both existing to parallelise fit().
A process that has run any OpenMP region -- including plain
C<score_samples()>/C<predict()> with the default C<use_openmp> -- and
then C<fork()>s (as C<parallel_fit> does) hands each child a copy of
libgomp's thread pool whose worker threads did not survive the fork. A
child that then starts its own C<#pragma omp parallel> region (as
C<use_openmp_fit> would) tries to reuse that now-invalid pool and
hangs. This is a general limitation of combining C<fork()> with OpenMP,
not something fixable from Perl, so C<parallel_fit>'s forked workers
always use the single-threaded C builder regardless of
C<use_openmp_fit> -- setting both just means C<parallel_fit> wins and
C<use_openmp_fit> has no effect for that call.
Detection happens once when the module is loaded. When the
distribution was installed with C<Inline::C> available, the C backend
was already compiled during C<make> and the installed object is loaded
directly (see L</Compile at install time (the prebuilt object)> below);
otherwise the backend is compiled on first load and the artefact is
cached under C<_Inline/> and reused on subsequent runs. Five package
variables report what the load picked up:
$Algorithm::Classifier::IsolationForest::HAS_C # 0/1
$Algorithm::Classifier::IsolationForest::HAS_OPENMP # 0/1
$Algorithm::Classifier::IsolationForest::HAS_SIMD # 0/1 (OpenMP 4.0+)
$Algorithm::Classifier::IsolationForest::OPT_LEVEL # e.g. "-O3 -march=native", '' if HAS_C is 0
$Algorithm::Classifier::IsolationForest::C_SOURCE # 'prebuilt' / 'runtime', '' if HAS_C is 0
Neither dependency is required. Without C<Inline::C> the module falls
back to a pure-Perl implementation that produces identical results, just
slower; without OpenMP the C backend runs single-threaded.
The bundled C<iforest accel> subcommand performs a tiny fit + score and
prints which backend is active (including the build flags below), which
is the recommended way to verify the build picked up the optional
dependencies on a given machine.
=head2 Compile at install time (the prebuilt object)
When C<Inline::C> is usable while the distribution itself is being
built, C<perl Makefile.PL> arranges for the C backend to be compiled
once during C<make> and installed alongside the module like any XS
object. At run time that object is loaded directly through
L<XSLoader>: no C compiler, no C<Inline> modules, and no C<_Inline/>
cache directory are needed on the machine the module ends up running
on, and the first-load compile pause disappears entirely.
On x86-64 hardware from roughly the last decade,
C<IF_ARCH=x86-64-v3 perl Makefile.PL> is a reasonable configure line:
it bakes AVX2 + FMA (without AVX-512) into the prebuilt object, which
can speed up extended-mode scoring (how much is hardware-dependent --
benchmark with C<iforest bench> before assuming) while avoiding the
C<-march=native> caveats described under L</Tuning the C build>.
Bit-for-bit result parity with the pure-Perl backend is preserved
either way (see C<IF_ARCH> below).
The C<IF_*> build flags described below are captured when
C<perl Makefile.PL> runs -- set them in the environment of I<that>
command, not of C<make> -- and recorded in the generated
C<Algorithm::Classifier::IsolationForest::BuildFlags> module, which
thereby also fixes what the prebuilt object was compiled with. At run
time the recorded values serve as the defaults, so a process started
with no C<IF_*> variables set uses the prebuilt object as-is.
Setting C<IF_*> variables at run time keeps working exactly as in
releases without prebuilt support: if the requested flags differ from
the recorded ones, the prebuilt object (compiled with the wrong flags
for the request) is skipped and the module compiles at first load into
C<_Inline/> -- which does need C<Inline::C> and a compiler on that
machine. Two related knobs exist:
=over 4
=item * C<IF_RUNTIME_BUILD=1> -- ignore the prebuilt object
unconditionally and compile at first load even though the requested
flags match the recorded ones. Useful when the installed object is
suspect (built on a different CPU than it now runs on, linked against a
libgomp that has since changed) or to A/B a fresh local build against
the shipped one.
=item * C<IF_INSTALL_BUILD=1> -- internal; set by the generated
Makefile rule that performs the install-time compile. Not meant for
manual use.
=back
If the prebuilt object cannot be loaded for any reason (deleted, built
against a different perl, version mismatch after an upgrade), the
module quietly falls through the same chain as always: runtime
Inline::C build first, pure Perl last.
=head2 Tuning the C build
These environment variables are read once, the first time the module is
loaded, so they must be set before that -- e.g. in the shell before
running a script, not via C<%ENV> inside the script itself. They are
also read by C<perl Makefile.PL> to pick the flags baked into the
prebuilt object (see above); at run time they override the recorded
configure-time values, at the price of a runtime compile.
=over 4
=item * C<IF_NO_C=1> -- skip attempting to build the C backend entirely.
Equivalent to constructing every instance with C<use_c =E<gt> 0>, but
without needing to touch every call site; useful for a clean pure-Perl
timing baseline, or to avoid the compile attempt's overhead/noise on a
host known to lack a C compiler (the attempt already fails gracefully
without this, so it's a convenience, not a correctness fix).
=item * C<IF_OPT=-O2> (or C<-O0>/C<-O1>/C<-Os>/C<-Og>/C<-Oz>) -- override
the default C<-O3>, e.g. to shorten build time while iterating, or work
around a miscompile on an unusual toolchain. Invalid values are ignored
with a warning rather than passed through, since this string reaches a
compiler command line.
=item * C<IF_ARCH=E<lt>valueE<gt>> -- adds C<-march=E<lt>valueE<gt>> so the
compiler can target specific instruction-set extensions (AVX2 gather +
FMA, etc.) for the extended-mode oblique dot product and the fit-time
min/max scan's C<#pragma omp simd> loops. Accepts values like
C<x86-64-v3>, C<skylake>, or C<znver3> -- whatever your compiler's
C<-march=> accepts. Also validated (a restricted character set, not
passed through as-is) for the same reason as C<IF_OPT>. The special
value C<none> (or an empty string) opts out of any arch recorded at
configure time, yielding a plain build. Whenever a C<-march> is in
effect the build also adds C<-ffp-contract=off>: with FMA available
the compiler would otherwise contract C<a*b+c> into fused
multiply-adds whose different rounding breaks the guarantee that
C<use_c =E<gt> 1> and C<use_c =E<gt> 0> build bit-identical trees (the
C<-march> speedup comes from vectorization, not contraction, so this
costs essentially nothing).
=item * C<IF_NATIVE=1> -- shorthand for C<IF_ARCH=native>; ignored if
C<IF_ARCH> is also set. Prefer a specific C<IF_ARCH> value over this on
a machine you don't control exclusively (a shared build host, a
container base image): blanket C<-march=native> pulls in whatever
instruction sets the build host happens to have, including AVX-512 on
some Intel CPUs -- which is known to trigger clock throttling under
sustained heavy use and can make throughput I<worse> than a
conservative target like C<x86-64-v3> (AVX2, no AVX-512). If in doubt,
benchmark both before committing to one.
=item * C<IF_NO_OPENMP=1> -- build (or select) the serial C backend: the
OpenMP compile attempt is skipped entirely, so the resulting object has
no libgomp linkage and never starts an OpenMP runtime inside the
process. This differs from C<OMP_NUM_THREADS=1>, which merely runs the
parallel code on one thread but still loads libgomp. Set at
C<perl Makefile.PL> time it yields a serial prebuilt object; set at run
time against an OpenMP prebuilt install it triggers a runtime serial
build (needing a compiler). An explicit C<IF_NO_OPENMP=0> re-enables
OpenMP over a serial configure-time default.
=back
Whichever of these are used, the cached artefact under C<_Inline/> is
pinned to that build's instruction set -- delete C<_Inline/> (or use a
separate one per host) if the directory is shared across machines with
different CPUs, or a stale binary built for a narrower instruction set
than the current host will simply keep being reused.
=head2 Tuning the OpenMP runtime
These are standard OpenMP environment variables libgomp already reads
at run time (set before running your script, no module-specific
handling needed) -- listed here because they matter most for exactly
the workloads this module has: C<score_all_xs>'s per-point parallel
loop and C<use_openmp_fit>'s per-tree parallel loop.
=over 4
=item * C<OMP_NUM_THREADS=N> -- caps how many threads a parallel region
uses. Useful to leave headroom for other work sharing the machine, or
to pin down C<use_openmp_fit> reproducibility checks (see its docs
above: results don't depend on this, but it's a natural thing to vary
when confirming that).
=item * C<OMP_PROC_BIND=close> / C<OMP_PLACES=cores> -- on multi-socket
or otherwise NUMA machines, pins each thread to a core near where its
data already lives instead of letting the OS scheduler migrate threads
across sockets mid-run. Both C<score_all_xs> (each thread scans its own
slice of the packed query buffer) and C<use_openmp_fit> (each thread
builds one tree from packed training data) benefit from this when the
input is large enough to not fit comfortably in one socket's cache.
=back
These cost nothing to try -- unlike C<IF_ARCH>/C<IF_NATIVE>, they're
read fresh every run, not baked into a cached binary, so there's no
downside to experimenting per invocation.
=head1 GENERAL METHODS
=head2 new(%args)
Inits the object.
- n_trees :: number of isolation trees in the ensemble
default :: 100
- sample_size :: sub-sample size used to build each tree... max samples
default :: 256
- max_depth :: per-tree height limit... if not defined is set to ceil(log2(psi))
default :: undef
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
Trees are built identically in both settings -- only aggregation
changes -- so the knob composes with either mode (axis or
extended) and an existing model may switch it after the fact with
set_voting() (which relearns a contamination threshold for the
new mode).
mean :: classic Isolation Forest: a sample's path lengths
across all trees are averaged and normalised into
one anomaly score; predict() thresholds that score.
majority :: Majority Voting Isolation Forest (MVIForest;
Chabchoub, Togbe, Boly & Chiky 2022 -- see
REFERENCES). Each tree scores the sample on its own
(s_i = 2**(-h_i / c(psi))) and votes it anomalous
when s_i >= the decision threshold; predict() flags
the sample when more than half of the trees
(int(n_trees/2) + 1) vote anomalous, and stops
walking trees per sample as soon as the outcome is
decided. The threshold argument/default of the
predict methods is therefore the PER-TREE cutoff
here, not a forest-level score cutoff.
score_samples() returns the fraction of trees
voting anomalous -- still in [0, 1], but discrete
in steps of 1/n_trees. contamination composes: fit()
learns the per-tree cutoff that flags the requested
fraction of the training set.
default :: mean
- parallel_fit :: positive integer N => build the trees across N forked
worker processes during fit(). Each worker gets a derived seed
(parent seed + worker_id * 1009) so the parallel fit is
reproducible across runs at fixed worker count -- but the trees
produced are NOT bit-identical to a serial fit with the same
seed, because the RNG draws happen in a different order.
Inference is unaffected. Falls back silently to serial on
platforms without a real fork() (e.g. Windows without Cygwin).
default :: undef (serial)
- use_c :: boolean, override whether the Inline::C backend is used for
both scoring and fit()'s tree builder. When false the instance
falls back to pure Perl for both even if the C backend compiled
successfully. When true (or unset) the C backend is used if
available ($HAS_C). fit() with use_c on produces bit-identical
trees to use_c off for the same seed -- only build speed differs.
default :: $HAS_C
- use_openmp :: boolean, override whether OpenMP parallel scoring is
used inside score_all_xs(). When false the C tree walk runs
single-threaded even if OpenMP was linked in. Ignored when
use_c is false (pure Perl has no OpenMP path).
default :: $HAS_OPENMP
- use_openmp_fit :: boolean, build fit()'s trees across OpenMP threads
(one tree per thread) instead of the single-threaded C builder.
Opt-in and off by default: unlike use_c/use_openmp, this changes
which trees get built. Perl's RNG isn't safe to call from
multiple OS threads sharing one interpreter, so this path seeds
an independent PRNG per tree from the tree index rather than
Drand01() -- trees differ from the use_c (single-threaded)
and pure-Perl paths even with the same seed, though a fixed
seed and n_trees still reproduce the same trees regardless of
OMP_NUM_THREADS or scheduling. Does NOT compose with
parallel_fit: a forked child starting its own OpenMP region
after the parent process has used OpenMP for anything can
hang (a general fork()+libgomp limitation), so parallel_fit's
workers always use the single-threaded C builder regardless
of this setting -- setting both just means parallel_fit wins.
Ignored (clamped to 0) when use_c is false or OpenMP isn't
linked in.
default :: 0
- feature_names :: optional arrayref of per-feature labels enabling the
*_tagged methods (and required by mungers below).
default :: undef
- mungers :: optional hashref of declarative L<Algorithm::ToNumberMunger>
specs, keyed as that module's compile() expects (scalar mungers by
their output tag, expanding mungers by any label with an 'into'
list, combining mungers by their output tag with a 'from' list).
When set, every tagged row -- the *_tagged methods, fit_tagged,
and tagged_row_to_array -- is munged from raw values (strings,
timestamps, status codes, ...) into numbers through the compiled
plan, and munge_rows() applies the scalar mungers to positional
rows. Requires feature_names; the plan compiles against them, so
any spec error croaks here in new(). Algorithm::ToNumberMunger is
an optional dependency, required only when a spec is given (or a
loaded model carrying one is used with tagged data). The spec is
saved with the model, so a loaded model munges scoring input
exactly as it did training input. See L</MUNGERS> for details
and caveats.
default :: undef
- schema_version :: optional opaque string identifying the revision of
the variable schema this model was built against. Never parsed
or compared numerically; saved with the model and shown by
`iforest info`. Usually set from a prototype (see
L</PROTOTYPES>) rather than passed directly.
default :: undef
- schema_description :: optional opaque free-text description of what
the variable schema is. Same handling as schema_version.
default :: undef
- feature_descriptions :: optional hashref of 'feature name => free
text' describing individual features. Requires feature_names;
every key must name an entry there (a description for a feature
that does not exist croaks -- it is either a typo or a stale
leftover from a schema change). Partial coverage is fine.
Saved with the model and shown beside each tag by
`iforest info`.
default :: undef
Note: log2 under Perl is as below...
log($psi) / log(2)
=cut
sub new {
my ( $class, %args ) = @_;
my $mode = $args{mode} // 'axis';
croak "mode must be 'axis' or 'extended'"
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
return $cached if defined $cached;
require Config;
$cached
= ( ( $Config::Config{d_fork} || '' ) eq 'define' ) ? 1 : 0;
return $cached;
}
}
#-------------------------------------------------------------------------------
# Fork-based parallel tree builder. Used by fit() when parallel_fit > 1
# and the platform has a real fork(). Divides n_trees evenly among
# workers; each child seeds its own RNG ($seed + worker_id * 1009 so
# fixed-worker-count runs are reproducible), builds its share (via the
# C builder when _use_c is on, same as the non-parallel path), and
# returns the trees to the parent via Storable on a one-shot pipe.
#
# The trees that come back differ from a serial fit with the same seed
# because the RNG draws happen in a different order -- this is documented
# as part of the parallel_fit contract.
#-------------------------------------------------------------------------------
sub _fit_trees_parallel {
my ( $self, $data, $psi, $limit, $workers ) = @_;
require Storable;
require POSIX;
my $n_trees = $self->{n_trees};
$workers = $n_trees if $workers > $n_trees;
# Divide n_trees as evenly as possible across workers.
my @shares;
{
my $base = int( $n_trees / $workers );
my $extras = $n_trees - $base * $workers;
for my $w ( 0 .. $workers - 1 ) {
push @shares, $base + ( $w < $extras ? 1 : 0 );
}
}
my @procs; # { pid, rh, share }
for my $w ( 0 .. $workers - 1 ) {
my $share = $shares[$w];
next unless $share > 0;
pipe( my $rh, my $wh ) or croak "pipe failed: $!";
my $pid = fork();
croak "fork failed: $!" unless defined $pid;
if ( $pid == 0 ) {
# child
close $rh;
binmode $wh;
if ( defined $self->{seed} ) {
srand( $self->{seed} + $w * 1009 );
}
# Deliberately never _build_forest_openmp here, even when
# use_openmp_fit is on: if this process (or the parent that
# fork()ed us) already ran any OpenMP region before this
# fork -- including plain score_samples()/predict() with
# the default use_openmp -- libgomp's thread pool exists
# but its worker threads didn't survive the fork. A child
# starting its own #pragma omp parallel region then tries
# to reuse that now-invalid pool and hangs. This is a
# general fork()+libgomp limitation, not fixable from here,
# so forked workers always use the single-threaded C
# builder (or pure Perl) instead. See t/03-fit-determinism.t
# and the NATIVE ACCELERATION docs for the observed hang and
# why parallel_fit + use_openmp_fit isn't composed for real.
my $trees;
if ( $self->{_use_c} ) {
$trees = $self->_build_forest_c( $data, $psi, $limit, $share );
} else {
my @t;
for ( 1 .. $share ) {
my $sample = _subsample( $data, $psi );
push @t, $self->_build_tree( $sample, 0, $limit );
}
$trees = \@t;
}
print $wh Storable::freeze($trees);
close $wh;
# _exit so we don't run parent END/DESTROY in the child.
POSIX::_exit(0);
} ## end if ( $pid == 0 )
close $wh;
binmode $rh;
push @procs, { pid => $pid, rh => $rh, share => $share };
} ## end for my $w ( 0 .. $workers - 1 )
# Collect from each pipe in worker order so the canonical tree
# ordering is deterministic (worker 0's trees first, then 1's, ...).
my @all_trees;
for my $p (@procs) {
my $buf;
{
local $/;
$buf = readline( $p->{rh} );
}
close $p->{rh};
waitpid( $p->{pid}, 0 );
my $exit = $? >> 8;
croak "parallel_fit worker $p->{pid} exited with status $exit"
if $exit != 0;
my $trees = eval { Storable::thaw($buf) };
croak "parallel_fit worker $p->{pid} returned unparseable trees: $@"
if $@ || ref $trees ne 'ARRAY';
push @all_trees, @$trees;
} ## end for my $p (@procs)
return \@all_trees;
} ## end sub _fit_trees_parallel
#-------------------------------------------------------------------------------
# C-accelerated fit(): builds $n_trees trees against $data (a subset or
# the full training set) via build_forest_xs, which does its own
# per-tree subsampling internally. Random draws inside the C builder
# go through Drand01() -- the same generator Perl's rand() uses -- in
# the same call order _subsample/_build_tree used, so the returned
# trees are bit-identical to what the pure-Perl path would build from
# the same RNG state. That's what lets fit() switch backends on the
# existing `use_c` knob instead of a new one.
#-------------------------------------------------------------------------------
sub _build_forest_c {
my ( $self, $data, $psi, $limit, $n_trees ) = @_;
my $n = scalar @$data;
my $nf = $self->{n_features};
my $x_packed = "\0" x ( $n * $nf * 8 );
my ( $mode, $fill ) = $self->_pack_args;
pack_input_xs( $data, $x_packed, $n, $nf, $mode, $fill );
my $mode_flag = $self->{mode} eq 'extended' ? 1 : 0;
my $ext_level = $self->{extension_level_used} // 0;
my $trees = [];
build_forest_xs( $x_packed, $n, $nf, $n_trees, $psi, $limit, $mode_flag, $ext_level, $trees );
return $trees;
} ## end sub _build_forest_c
#-------------------------------------------------------------------------------
# OpenMP-parallel fit(): builds $n_trees trees across OpenMP threads (one
# tree per thread) via build_forest_openmp_xs. Unlike _build_forest_c,
# random draws come from a thread-private PRNG seeded per tree index
# rather than Drand01() -- Perl's RNG state can't be shared safely
# across OpenMP threads -- so the resulting trees are NOT bit-identical
# to the use_c (serial) or pure-Perl paths for the same seed, though a
# fixed seed + n_trees still reproduce the same trees regardless of
# OMP_NUM_THREADS. This is why it's gated by the separate, opt-in
# use_openmp_fit knob rather than reusing use_c/use_openmp.
#
# Only called from fit()'s non-forked branch. _fit_trees_parallel's
# workers never call this, even when use_openmp_fit is on: a forked
# child starting its own OpenMP region after the parent process has
# used OpenMP for anything (this includes plain score_samples()) can
# hang -- see the comment above that branch for the fork()+libgomp
# hazard this avoids.
#
# build_forest_openmp_xs hands back three arrayrefs of per-tree packed
# buffers (the same SoA layout _pack_tree produces) instead of Perl tree
# structures -- that's how it avoids any Perl API call inside its
# parallel region. _unpack_forest converts them back into the ordinary
# nested-arrayref tree shape so to_json/from_json/_rebuild_c_trees don't
# need to know this path exists.
#-------------------------------------------------------------------------------
sub _build_forest_openmp {
my ( $self, $data, $psi, $limit, $n_trees ) = @_;
my $n = scalar @$data;
my $nf = $self->{n_features};
my $x_packed = "\0" x ( $n * $nf * 8 );
my ( $mode, $fill ) = $self->_pack_args;
pack_input_xs( $data, $x_packed, $n, $nf, $mode, $fill );
my $mode_flag = $self->{mode} eq 'extended' ? 1 : 0;
my $ext_level = $self->{extension_level_used} // 0;
my ( @nodes, @idx, @val );
build_forest_openmp_xs( $x_packed, $n, $nf, $n_trees, $psi, $limit,
$mode_flag, $ext_level, \@nodes, \@idx, \@val, 1 );
return _unpack_forest( \@nodes, \@idx, \@val );
} ## end sub _build_forest_openmp
# Inverse of _pack_tree's SoA layout: given one tree's packed node
# buffer plus the shared idx/val coefficient buffers, reconstructs the
# ordinary nested-arrayref tree structure _build_tree/_build_node_c
# produce. li/ri fields hold the child's absolute node index, so this
# just follows them recursively from whatever index the caller says the
# root lives at. NOTE: _pack_tree numbers nodes DFS pre-order (root at
# 0), but build_forest_openmp_xs appends nodes post-order (children
# before parent), putting the root LAST -- the caller must pass the
# right root index for the buffer's origin.
sub _unpack_node {
my ( $nodes, $idx, $val, $node_i ) = @_;
my $off = $node_i * 6;
my $type = $nodes->[$off];
if ( $type == 0 ) {
return [ _NODE_LEAF, int( $nodes->[ $off + 1 ] ) ];
} elsif ( $type == 1 ) {
my ( $attr, $split, $li, $ri )
= @{$nodes}[ $off + 1 .. $off + 4 ];
return [
_NODE_AXIS, int($attr), $split,
_unpack_node( $nodes, $idx, $val, int($li) ),
_unpack_node( $nodes, $idx, $val, int($ri) ),
];
} else {
my ( $coff, $num, $li, $ri, $b ) = @{$nodes}[ $off + 1 .. $off + 5 ];
$coff = int($coff);
$num = int($num);
return [
_NODE_OBLIQUE,
[ @{$idx}[ $coff .. $coff + $num - 1 ] ],
( run in 0.879 second using v1.01-cache-2.11-cpan-9581c071862 )