Algorithm-Classifier-IsolationForest
view release on metacpan or search on metacpan
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
/* C99 VLAs -- n_trees is small (typ. 100) and fits on the stack. */
const double *node_ptrs[n_trees];
const int *idx_ptrs[n_trees];
const double *val_ptrs[n_trees];
/* forest_bytes totals every buffer the tree walks touch; it decides
* between the two loop shapes below. */
size_t forest_bytes = 0;
for (ti = 0; ti < n_trees; ti++) {
SV** np = av_fetch(nodes_av, ti, 0);
SV** ip = av_fetch(idx_av, ti, 0);
SV** vp = av_fetch(val_av, ti, 0);
if (!np || !*np || !ip || !*ip || !vp || !*vp) {
croak("score_all_xs: missing tree %d", ti);
}
node_ptrs[ti] = (const double*)SvPVbyte(*np, tl); forest_bytes += tl;
idx_ptrs[ti] = (const int*) SvPVbyte(*ip, tl); forest_bytes += tl;
val_ptrs[ti] = (const double*)SvPVbyte(*vp, tl); forest_bytes += tl;
}
xd = (const double*)SvPVbyte(x_sv, tl);
sm = (double*)SvPVbyte_force(sm_sv, tl);
/* Two loop shapes over the same per-point ascending-t additions --
* bit-identical results either way, so the size heuristic choosing
* between them can never change scores.
*
* Point-major (small forests): each point walks all trees with its
* path-length sum held in a register. Cheapest per walk, and the
* whole forest stays cache-resident across points anyway.
*
* Tree-blocked (large forests): once the forest outgrows L3, the
* point-major loop re-streams every tree's nodes and coefficients
* from memory for every point -- an extended-mode tree is ~56 KB
* at 16 features (24 KB nodes + 32 KB dense coefficients), and its
* per-tree scoring cost measured 2.2x worse at 400 trees than at
* 100. Walking a block of points through ONE tree at a time keeps
* that tree hot in L1/L2 while the block's rows stream through it
* (measured 3.1x faster at 400 extended trees, 20k points). The
* blocked shape pays an sm[i] load+store per walk instead of a
* register add, which measurably hurts cheap axis walks while the
* forest still fits in cache -- hence the byte threshold rather
* than always tiling. */
if (forest_bytes <= (size_t)4 * 1024 * 1024) {
#ifdef _OPENMP
#pragma omp parallel for schedule(static) if(use_openmp)
#endif
for (int i = 0; i < n_pts; i++) {
const double *xi = xd + (size_t)i * (size_t)n_feats;
double sum = 0.0;
for (int t = 0; t < n_trees; t++) {
sum += if_walk_tree(node_ptrs[t], idx_ptrs[t],
val_ptrs[t], xi, n_feats);
}
sm[i] = sum;
}
}
else {
/* 256 rows x 16 features x 8 bytes = 32 KB of input per block
* -- comfortable in L2 next to one tree. Each OpenMP thread
* owns whole blocks and therefore a unique slice of sm[], so
* there is still no synchronisation. For small batches the
* tile shrinks to keep ~4 blocks per thread available; losing
* per-block tree reuse there is fine, since a small batch
* never re-streams much anyway. */
int tile = 256;
#ifdef _OPENMP
if (use_openmp) {
int min_blocks = omp_get_max_threads() * 4;
if (min_blocks > 0 && (n_pts + tile - 1) / tile < min_blocks) {
tile = (n_pts + min_blocks - 1) / min_blocks;
if (tile < 1) tile = 1;
}
}
#endif
int n_blocks = (n_pts + tile - 1) / tile;
#ifdef _OPENMP
#pragma omp parallel for schedule(static) if(use_openmp)
#endif
for (int blk = 0; blk < n_blocks; blk++) {
const int i0 = blk * tile;
const int i1 = (i0 + tile < n_pts) ? i0 + tile : n_pts;
for (int i = i0; i < i1; i++) sm[i] = 0.0;
for (int t = 0; t < n_trees; t++) {
const double *nd = node_ptrs[t];
const int *ico = idx_ptrs[t];
const double *vco = val_ptrs[t];
for (int i = i0; i < i1; i++) {
sm[i] += if_walk_tree(nd, ico, vco,
xd + (size_t)i * (size_t)n_feats,
n_feats);
}
}
}
}
}
/* vote_all_xs(nodes_av, idx_av, val_av, x_sv, sm_sv,
* n_pts, n_feats, n_trees, depth_cut, min_votes, use_openmp)
*
* Majority-voting (MVIForest) tree walk: a tree votes a point anomalous
* when the point's path length in that tree is <= depth_cut -- the
* depth-domain image of the per-tree score cutoff (the Perl side
* precomputes depth_cut = -c(psi) * log2(threshold), so no per-tree
* exp()/log() runs in here).
*
* min_votes == 0: sm[i] = the point's full vote count over all n_trees
* trees (a small integer stored as a double, so the existing
* finalize_* helpers work on the buffer unchanged).
* min_votes > 0: sm[i] = 1.0/0.0 anomaly label, with per-point early
* exit: the walk stops as soon as the point has min_votes votes (the
* remaining trees can't change the outcome) or can no longer reach
* min_votes. This is MVIForest's "stop at majority" scoring loop.
*
* Always point-major, unlike score_all_xs's two loop shapes: the vote
* count / early exit is per-point state, so a tree-blocked loop would
* have to re-load it per walk and could never exit a point early.
* Votes are integer counts, so there is no summation-order concern
* either way. Thread-safety matches score_all_xs: the parallel region
* reads extracted pointers and writes a unique sm[i] per iteration. */
void vote_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,
double depth_cut, int min_votes, int use_openmp){
STRLEN tl;
AV *nodes_av, *idx_av, *val_av;
const double *xd;
double *sm;
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
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);
left = _build_node_c(aTHX_ x, nf, lidx, ln, depth + 1, limit,
mode_flag, ext_active);
right = _build_node_c(aTHX_ x, nf, ridx, rn, depth + 1, limit,
mode_flag, ext_active);
result = _mk_oblique(aTHX_ pool, coef, active, b, left, right);
free(pool); free(coef);
}
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;
}
lo = (double*)malloc(nf * sizeof(double));
hi = (double*)malloc(nf * sizeof(double));
for (f = 0; f < nf; f++) {
lo[f] = HUGE_VAL;
hi[f] = -HUGE_VAL;
}
for (int i = 0; i < size; i++) {
const double* row = x + (size_t)idxs[i] * (size_t)nf;
/* See the matching comment in _build_node_c: no isnan() guard
* needed, since NaN < x / NaN > x are always false already --
* that's what lets this vectorize as a plain min/max scan.
* omp simd here is thread-safe to call from inside the caller's
* omp parallel region: it's a per-thread vectorization hint,
* not a team construct, so it doesn't nest into anything. */
#ifdef _OPENMP
#pragma omp simd
#endif
for (int f2 = 0; f2 < nf; f2++) {
double v = row[f2];
if (v < lo[f2]) lo[f2] = v;
if (v > hi[f2]) hi[f2] = v;
}
}
varying = (int*)malloc(nf * sizeof(int));
nv = 0;
for (f = 0; f < nf; f++) {
if (lo[f] < hi[f]) varying[nv++] = f;
}
if (nv == 0) {
free(lo); free(hi); free(varying);
my_idx = tb_push_node(buf, 0.0, (double)size, 0.0, 0.0, 0.0, 0.0);
free(idxs);
return my_idx;
}
if (mode_flag == 0) {
int attr = varying[(int)(sm64_drand(rng) * nv)];
double split = lo[attr] + sm64_drand(rng) * (hi[attr] - lo[attr]);
int *lidx = (int*)malloc(size * sizeof(int));
int *ridx = (int*)malloc(size * sizeof(int));
int ln = 0, rn = 0, i, li, ri;
for (i = 0; i < size; i++) {
int row = idxs[i];
double v = x[(size_t)row * (size_t)nf + attr];
if (v < split) 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);
my_idx = tb_push_node(buf, 1.0, (double)attr, split,
(double)li, (double)ri, 0.0);
} else {
int active = ext_active + 1;
int *pool, *lidx, *ridx;
double *coef;
double b = 0.0;
int ln = 0, rn = 0, i, k, li, ri, coff;
if (active > nv) active = nv;
pool = (int*)malloc(nv * sizeof(int));
memcpy(pool, varying, nv * sizeof(int));
for (i = 0; i < active; i++) {
int j = i + (int)(sm64_drand(rng) * (nv - i));
int tmp = pool[i]; pool[i] = pool[j]; pool[j] = tmp;
}
coef = (double*)malloc(active * sizeof(double));
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
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
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
# 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,
LIBS => '-lm -lgomp',
@install,
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
my ($s, $l) = $iforest->score_predict_split(\@data, 0.6);
=head1 DESCRIPTION
Isolation Forest (Liu, Fei Tony & Ting, Kai & Zhou, Zhi-Hua, 2008) detects anomalies by random
partitioning rather than by modelling normal points. Each tree repeatedly
splits the data. Points that get isolated after only a few splits are likely
anomalies. The score is the average isolation depth across many trees,
normalised so values approach 1 for anomalies and stay below 0.5 for normal
points.
In extended mode the module implements the Extended Isolation Forest
variant. Each split is a random hyperplane instead of an axis-aligned cut,
which removes the rectangular, axis-aligned bias in the score field and
tends to help on elongated or multi-modal data.
With C<< voting => 'majority' >> the module implements the Majority Voting
Isolation Forest (MVIForest) aggregation: each tree votes a sample
anomalous or normal against the decision threshold and the label is the
majority of the votes, with prediction stopping early once the majority is
reached. Trees are built identically either way, so this composes with
both axis and extended mode, and an existing model can be flipped between
the two modes with L</set_voting> without refitting; see C<voting> under
L</new(%args)>.
psi referenced below is ψ or the pitchfork math symbol referenced in the paper,
Liu, Fei Tony & Ting, Kai & Zhou, Zhi-Hua. (2008). Isolation Forest. 413 - 422. 10.1109/ICDM.2008.17.
... or max samples.
L<https://www.researchgate.net/publication/224384174_Isolation_Forest>
=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
- seed :: optional integer to seed srand with for reproducible trees...
see perldoc -f srand for more info. This number is processed via abs(int()).
default :: undef
- mode :: if it should be IF or EIF
axis :: classic axis-parallel splits (IF)
extended :: oblique hyperplane splits (EIF)
default :: axis
- extension_level :: extended mode only... how many features take partin each
split. 0 behaves like a single-feature (axis) cut; the
maximum (n_features - 1) uses every varying feature. undef
=> maximum. Clamped to [0, n_features - 1] at fit time.
- contamination :: expected fraction of anomalies, in (0, 0.5]. When given,
fit() learns a score threshold that flags this fraction of
the training set, and predict() uses it by default. undef
=> no learned threshold (predict() falls back to 0.5).
default :: undef
- missing :: how fit() treats undef (missing) feature cells. Scoring always
tolerates undef regardless of this setting; it governs fit().
die :: croak from fit() if the training data contains any
undef cell. Scoring still maps undef to 0 (the
long-standing behaviour), so a model fitted on clean
data can still score rows with missing features.
zero :: treat a missing cell as the value 0, at fit and score.
impute :: replace a missing cell with the per-feature mean (or
median, see impute_with) learned from the present
values at fit time. The fill vector is stored on the
model and reused for scoring and persistence.
nan :: build feature ranges from present values only and route
a point missing the split feature to the right child,
consistently at fit and score time. Missingness is
preserved as signal rather than filled.
default :: die
- impute_with :: 'mean' or 'median'; the statistic used to compute the
per-feature fill under missing => 'impute'. Ignored otherwise.
default :: mean
- voting :: how the per-tree results are aggregated at scoring time.
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
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'"
unless $mode eq 'axis' || $mode eq 'extended';
# How fit() treats undef (missing) feature cells. Scoring always
# tolerates undef regardless of this setting -- it governs fit only.
# die :: croak if the training data contains any undef cell (default)
# zero :: treat a missing cell as the value 0
# impute :: replace a missing cell with the per-feature mean/median
# learned from the present values at fit time
# nan :: build ranges over present values only and route a point
# missing the split feature consistently to one branch, at
# both fit and score time
my $missing = $args{missing} // 'die';
croak "missing must be one of: die, zero, impute, nan"
unless $missing =~ /\A(?:die|zero|impute|nan)\z/;
my $impute_with = $args{impute_with} // 'mean';
croak "impute_with must be 'mean' or 'median'"
unless $impute_with =~ /\A(?:mean|median)\z/;
# How per-tree results are aggregated at scoring time. Trees are
# built identically either way -- this knob never touches fit()'s
# forest, only how score/predict combine the per-tree path lengths.
# mean :: classic IForest: average path length across trees,
# normalised into one score (the default)
# majority :: MVIForest (Chabchoub et al. 2022): each tree votes
# anomalous/normal against the decision threshold and
# the label is the majority of the tree votes
my $voting = $args{voting} // 'mean';
croak "voting must be 'mean' or 'majority'"
unless $voting =~ /\A(?:mean|majority)\z/;
if ( defined( $args{seed} ) ) {
$args{seed} = abs( int( $args{seed} ) );
}
# Clamp the accel knobs against what the build actually has. Passing
# use_c => 1 on a machine where Inline::C never compiled would otherwise
# leave score_samples() calling an undefined XS sub at first use.
# OpenMP is meaningless without the C tree walk, so force it off
# whenever the C backend is off -- matches the documented
# "Ignored when use_c is false" semantics.
my $use_c
= defined $args{use_c}
? ( $args{use_c} && $HAS_C ? 1 : 0 )
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
# Slot 2 carries c(size) precomputed, so the C scoring loop
# adds it straight to the depth instead of paying a log()
# per point per tree at every leaf hit. _c is the same
# function the pure-Perl scorer uses, so both backends keep
# producing bit-identical path lengths.
$node_data[$my_idx] = [ 0.0, $node->[1] + 0.0, _c( $node->[1] ), 0.0, 0.0, 0.0 ];
} elsif ( $node->[0] == _NODE_AXIS ) {
my $li = $assign->( $node->[3] );
my $ri = $assign->( $node->[4] );
$node_data[$my_idx] = [
1.0,
$node->[1] + 0.0, # attr
$node->[2] + 0.0, # split
$li + 0.0,
$ri + 0.0,
0.0,
];
} else { # _NODE_OBLIQUE
my ( $idx_arr, $coef_arr, $b ) = ( $node->[1], $node->[2], $node->[3] );
my $coef_off = scalar @coef_idx;
my $num = scalar @$idx_arr;
# Dense-pack opportunity: when this oblique split uses
# every feature (extension_level == n_features - 1 and
# all features vary), pack the coefficients in feature
# order so val[k] is the coefficient for feature k. The
# C scoring path then detects `nf == n_feats` and switches
# to a no-gather inner loop (dot += val[k] * xi[k]) that
# auto-vectorizes cleanly with FMA.
if ( defined $n_features && $num == $n_features ) {
my %coef_for;
@coef_for{@$idx_arr} = @$coef_arr;
for my $k ( 0 .. $n_features - 1 ) {
push @coef_idx, $k;
push @coef_val, $coef_for{$k} + 0.0;
}
} else {
for my $i ( 0 .. $num - 1 ) {
push @coef_idx, int( $idx_arr->[$i] );
push @coef_val, $coef_arr->[$i] + 0.0;
}
}
my $li = $assign->( $node->[4] );
my $ri = $assign->( $node->[5] );
$node_data[$my_idx] = [ 2.0, $coef_off + 0.0, $num + 0.0, $li + 0.0, $ri + 0.0, $b + 0.0, ];
} ## end else [ if ( $node->[0] == _NODE_LEAF ) ]
return $my_idx;
}; ## end $assign = sub
$assign->($root);
my $nodes_packed = pack( 'd*', map { @$_ } @node_data );
my $idx_packed = @coef_idx ? pack( 'l*', @coef_idx ) : pack('l*');
my $val_packed = @coef_val ? pack( 'd*', @coef_val ) : pack('d*');
return ( $nodes_packed, $idx_packed, $val_packed );
} ## end sub _pack_tree
# Build packed C-ready representations for all trees and store them in
# $self->{_c_nodes}, $self->{_c_coef_idx}, $self->{_c_coef_val}.
# Called after fit() and from_json() when _use_c is true. n_features is
# threaded through so _pack_tree can spot the dense-pack opportunity.
sub _rebuild_c_trees {
my ($self) = @_;
my ( @c_nodes, @c_coef_idx, @c_coef_val );
for my $tree ( @{ $self->{trees} } ) {
my ( $np, $ip, $vp ) = _pack_tree( $tree, $self->{n_features} );
push @c_nodes, $np;
push @c_coef_idx, $ip;
push @c_coef_val, $vp;
}
$self->{_c_nodes} = \@c_nodes;
$self->{_c_coef_idx} = \@c_coef_idx;
$self->{_c_coef_val} = \@c_coef_val;
} ## end sub _rebuild_c_trees
sub _check_fitted {
my ($self) = @_;
croak "model is not fitted yet; call fit() first"
unless ref $self->{trees} eq 'ARRAY' && @{ $self->{trees} };
}
# Memoised "does this perl have a real fork()?". False on Windows
# without Cygwin; true on every Unix-like platform.
{
my $cached;
sub _fork_supported {
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) ),
( run in 0.681 second using v1.01-cache-2.11-cpan-c966e8aa7e8 )