Algorithm-Classifier-IsolationForest
view release on metacpan or search on metacpan
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
# into out_sv. Replaces the dominant per-call Perl map-pack loop.
# miss_mode selects how an undef cell is packed: 0 => 0.0, 1 => the
# per-feature fill from fill_sv (impute), 2 => NaN (nan strategy).
#
# score_all_xs(nodes_av, idx_av, val_av, x_sv, sm_sv,
# n_pts, n_feats, n_trees, use_openmp)
# Sums path lengths for all n_pts query points across all n_trees trees
# in one call. Outer loop over points is OpenMP-parallel when the
# module was built with OpenMP (each iteration writes to a unique sm[i],
# so no synchronisation is needed). Tree pointers are extracted from
# the AVs before the parallel region; the parallel region touches only
# raw int / double buffers.
#
# 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 (voting => 'majority') counterpart of score_all_xs:
# instead of summing path lengths it counts, per point, how many trees
# "vote anomalous" (path length <= depth_cut). min_votes == 0 writes
# the full vote count into sm[i]; min_votes > 0 writes a 0.0/1.0 label
# with per-point early exit once the majority outcome is decided --
# the MVIForest scoring loop. See the function's own comment.
#
# Node layout (6 doubles per node, "IF_NZ = 6"):
# leaf: [0, size, c(size), 0, 0, 0]
# axis: [1, attr, split, li, ri, 0]
# oblique: [2, coff, nf, li, ri, b]
#
# c(size) is the expected-path-length adjustment for a leaf holding
# `size` points, precomputed by _pack_tree (it involves a log(); doing
# it at pack time keeps transcendentals out of the per-point per-tree
# scoring loop). The fit-time TreeBuf writer leaves that slot 0 --
# its buffers are unpacked into Perl trees and re-packed by
# _pack_tree before score_all_xs ever sees them.
#
# Coefficient storage uses a Structure-of-Arrays layout: one int32 array
# per tree (feature indices, packed with 'l*') and one double array per
# tree (coefficients, packed with 'd*'). Both are indexed by `coff` --
# the same offset addresses paired entries in the two arrays. Splitting
# them this way halves index bandwidth, removes the per-element
# (int)<double> cast inside the SIMD loop, and lets the value loads be
# contiguous so the compiler emits a clean FMA chain over val[k] with
# the feature gather on xi[idx[k]] kept separate.
#
# Dense-pack fast path: when an oblique node uses every feature (the
# common case in extended mode with extension_level == n_features - 1),
# _pack_tree writes its coefficients in feature order so val[k] is the
# coefficient for feature k. score_all_xs detects this via `nf ==
# n_feats` and uses a no-gather dot product (dot += val[k] * xi[k])
# that vectorizes cleanly with FMA -- substantially faster than the
# sparse gather path on high-feature-count models.
# x: row-major doubles, n_pts rows of n_feats each.
# sums: out double array of length n_pts; score_all_xs writes once per i.
#
# OpenMP is enabled at module load when the toolchain accepts -fopenmp and
# libgomp is linkable; otherwise the same C code compiles to a serial loop
# (the #pragma is silently ignored without _OPENMP defined).
# ---------------------------------------------------------------------------
our $HAS_C = 0;
our $HAS_OPENMP = 0;
our $HAS_SIMD = 0;
our $OPT_LEVEL = ''; # the actual -O.../-march=... flags used to build, if any
our $C_SOURCE = ''; # 'prebuilt' (object installed at `make` time) or
# 'runtime' (compiled at first load into _Inline/);
# '' when $HAS_C is 0
{
my $C_CODE = <<'__INLINE_C__';
#include <math.h>
#include <string.h>
#include <stdint.h>
#ifdef _OPENMP
#include <omp.h>
#endif
#define IF_NZ 6
/* Data prefetch hint; a no-op on compilers without __builtin_prefetch.
* Purely a performance hint -- never affects results. */
#if defined(__GNUC__) || defined(__clang__)
#define IF_PREFETCH(p) __builtin_prefetch(p)
#else
#define IF_PREFETCH(p)
#endif
int has_openmp_xs(){
#ifdef _OPENMP
return 1;
#else
return 0;
#endif
}
/* SIMD on the extended-mode oblique dot product is enabled via
* `#pragma omp simd`, which OpenMP 4.0 (_OPENMP == 201307) introduced.
* Anything older silently ignores the pragma -- the loop still runs,
* just not auto-vectorised. So "simd available" really means the
* compiler is going to honour the pragma we put on that loop. */
int has_simd_xs(){
#if defined(_OPENMP) && _OPENMP >= 201307
return 1;
#else
return 0;
#endif
}
/* pack_input_xs(data_sv, out_sv, n_pts, n_feats, miss_mode, fill_sv)
*
* Walks a Perl arrayref-of-arrayrefs (n_pts rows of n_feats doubles each)
* directly in C and writes the packed double buffer into out_sv (which the
* caller pre-allocates with "\0" x (n_pts*n_feats*8)). Replaces
*
* pack('d*', map { my $r=$_; map { $r->[$_] // 0 } 0..$nf-1 } @$data)
*
* which was the dominant per-call overhead for high feature counts.
*
* miss_mode selects what an undef cell (or missing row) becomes:
* 0 => 0.0 (the 'die'/'zero' missing strategies)
* 1 => fill[k] (the 'impute' strategy; fill_sv is a packed
* double buffer of n_feats per-feature fill values)
* 2 => NaN (the 'nan' strategy; the C scorer's `<` / `<=`
* comparisons are both false for NaN, so a point
* missing the split feature falls to the right
* child -- matching how fit() routes it)
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
av_store(out, i, newRV_noinc((SV*)row));
}
}
/* vote_score_predict_split_xs(sm_sv, n_pts, t, min_votes,
* scores_rv, labels_rv)
*
* Parallel-arrays variant of vote_score_predict_xs, mirroring
* score_predict_split_xs's shape for the majority-voting path. */
void vote_score_predict_split_xs(SV* sm_sv, int n_pts, double t,
double min_votes,
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("vote_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(sm[i] / t));
av_store(labels, i, newSViv(sm[i] >= min_votes ? 1 : 0));
}
}
/* ---------------------------------------------------------------------
* build_forest_xs -- C-accelerated fit() tree builder.
*
* Replaces the pure-Perl _subsample + _build_tree + _axis_split /
* _oblique_split recursion with an equivalent C implementation that
* partitions plain `int` row-index arrays instead of copying arrayrefs
* of Perl SVs at every split. Random draws go through Drand01() --
* the exact generator Perl's own rand()/srand() use internally -- in
* the same call order the Perl code used, so a fit() with a given
* seed produces BIT-IDENTICAL trees whether use_c is on or off. This
* is what lets fit() reuse the existing `use_c` knob instead of a new
* one: switching backends never changes the model, only how fast it's
* built. (Verified by t/02-accel-selection.t's "identical seed =>
* identical trees" subtest, which exercises both backends.)
*
* Output trees are plain Perl arrayrefs in the same node shape
* _build_tree produces (leaf/axis/oblique -- see the file-top
* comment), so every downstream consumer (_pack_tree, to_json,
* from_json, the pure-Perl scorer) is unchanged.
*
* x_sv: packed row-major double buffer, n_pts rows of n_feats each
* (from pack_input_xs -- NaN marks a missing cell under the
* 'nan' missing-strategy).
* mode_flag: 0 => axis-parallel splits, 1 => oblique (extended).
* ext_level: extension_level_used (ignored when mode_flag == 0).
* out_rv: pre-existing arrayref; filled with n_trees tree roots.
* ------------------------------------------------------------------ */
/* Box-Muller normal draw, in the same rand() call order as _randn(). */
static double _c_randn(pTHX) {
double u1 = Drand01();
double u2;
if (u1 == 0.0) u1 = 1e-12;
u2 = Drand01();
return sqrt(-2.0 * log(u1)) * cos(6.283185307179586 * u2);
}
static SV* _mk_leaf(pTHX_ int size) {
AV* av = newAV();
av_extend(av, 1);
AvARRAY(av)[0] = newSVnv(0.0);
AvARRAY(av)[1] = newSViv(size);
AvFILLp(av) = 1;
return newRV_noinc((SV*)av);
}
static SV* _mk_axis(pTHX_ int attr, double split, SV* left, SV* right) {
AV* av = newAV();
av_extend(av, 4);
AvARRAY(av)[0] = newSVnv(1.0);
AvARRAY(av)[1] = newSViv(attr);
AvARRAY(av)[2] = newSVnv(split);
AvARRAY(av)[3] = left;
AvARRAY(av)[4] = right;
AvFILLp(av) = 4;
return newRV_noinc((SV*)av);
}
static SV* _mk_oblique(pTHX_ const int* idx, const double* coef, int n,
double b, SV* left, SV* right) {
AV *iav, *cav, *av;
int k;
iav = newAV();
cav = newAV();
if (n > 0) {
av_extend(iav, n - 1);
av_extend(cav, n - 1);
}
for (k = 0; k < n; k++) {
AvARRAY(iav)[k] = newSViv(idx[k]);
AvARRAY(cav)[k] = newSVnv(coef[k]);
}
AvFILLp(iav) = n - 1;
AvFILLp(cav) = n - 1;
av = newAV();
av_extend(av, 5);
AvARRAY(av)[0] = newSVnv(2.0);
AvARRAY(av)[1] = newRV_noinc((SV*)iav);
AvARRAY(av)[2] = newRV_noinc((SV*)cav);
AvARRAY(av)[3] = newSVnv(b);
AvARRAY(av)[4] = left;
AvARRAY(av)[5] = right;
AvFILLp(av) = 5;
return newRV_noinc((SV*)av);
}
/* Builds one node from the point set `idxs` (row indices into `x`,
* length `size`); recurses left-then-right, matching _build_tree's
* traversal order so nested splits draw random numbers in the same
* sequence the pure-Perl path would. Takes ownership of `idxs` --
* frees it before returning. */
static SV* _build_node_c(pTHX_ const double* x, int nf, int* idxs, int size,
int depth, int limit, int mode_flag,
int ext_active) {
double *lo, *hi;
int *varying, nv, f;
SV *result;
if (depth >= limit || size <= 1) {
SV* leaf = _mk_leaf(aTHX_ size);
free(idxs);
return leaf;
}
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;
/* No isnan() guard needed: NaN < x and NaN > x are always false
* under IEEE 754, so a NaN cell (the 'nan' missing strategy)
* already leaves lo/hi untouched without an explicit check --
* one less branch, and it's what lets this loop vectorize
* cleanly as a plain elementwise min/max scan. */
#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);
SV* leaf = _mk_leaf(aTHX_ size);
free(idxs);
return leaf;
}
if (mode_flag == 0) {
/* Axis-parallel split: one varying feature, one threshold. */
int attr = varying[(int)(Drand01() * nv)];
double split = lo[attr] + Drand01() * (hi[attr] - lo[attr]);
int *lidx = (int*)malloc(size * sizeof(int));
int *ridx = (int*)malloc(size * sizeof(int));
int ln = 0, rn = 0, i;
SV *left, *right;
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);
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_axis(aTHX_ attr, split, left, right);
} else {
/* Oblique split: a random hyperplane over `active` features. */
int active = ext_active + 1;
int *pool, *lidx, *ridx;
double *coef;
double b = 0.0;
int ln = 0, rn = 0, i, k;
SV *left, *right;
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)(Drand01() * (nv - i));
int tmp = pool[i]; pool[i] = pool[j]; pool[j] = tmp;
}
coef = (double*)malloc(active * sizeof(double));
for (k = 0; k < active; k++) {
int ff = pool[k];
double c = _c_randn(aTHX);
double p = lo[ff] + Drand01() * (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);
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;
}
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
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));
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) {
if (d->n == d->cap) {
size_t newcap = d->cap ? d->cap * 2 : 64;
d->v = (double*)realloc(d->v, newcap * sizeof(double));
d->cap = newcap;
}
d->v[d->n++] = x;
}
static void _dswap(double *a, double *b) { double t = *a; *a = *b; *b = t; }
/* Lomuto partition with a median-of-three pivot (avoids the O(n^2)
* worst case a fixed pivot hits on already-sorted or reverse-sorted
* input, which real feature columns -- timestamps, counters -- often
* are). Returns the pivot's final index. */
static int _partition_lomuto(double *a, int lo, int hi) {
int mid = lo + (hi - lo) / 2;
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
* adjustment) -- the C image of _score_row's loop, used by the
* prequential score_learn path where trees mutate between rows and the
* packed-snapshot scorer can never amortise. Draws nothing. */
double online_score_row_xs(SV* trees_av_sv, SV* row_sv, int nf, int eta) {
dTHX;
AV* trees;
double* x;
double sum = 0.0;
int t, n_trees;
if (!SvROK(trees_av_sv) || SvTYPE(SvRV(trees_av_sv)) != SVt_PVAV)
croak("online_score_row_xs: trees must be an arrayref");
trees = (AV*)SvRV(trees_av_sv);
n_trees = (int)av_len(trees) + 1;
x = (double*)malloc(nf * sizeof(double));
_ol_read_row(aTHX_ row_sv, x, nf);
for (t = 0; t < n_trees; t++) {
HV* tree = _ol_tree_hv(aTHX_ trees, t);
SV** rsv = hv_fetch(tree, "root", 4, 0);
AV* node;
SV** slots;
int depth = 0;
if (!rsv || !*rsv || !SvOK(*rsv)) continue;
node = (AV*)SvRV(*rsv);
slots = AvARRAY(node);
while (SvIV(slots[OL_TYPE]) != 0) {
int attr = (int)SvIV(slots[OL_ATTR]);
double split = SvNV(slots[OL_SPLIT]);
node = (AV*)SvRV(slots[(x[attr] < split) ? OL_LEFT : OL_RIGHT]);
slots = AvARRAY(node);
depth++;
}
sum += (double)depth + _ol_rpl((double)SvIV(slots[OL_COUNT]), eta);
}
free(x);
return sum;
}
__INLINE_C__
# IF_NO_C=1 skips even attempting to set up the C backend -- useful for
# forcing the pure-Perl path without touching every constructor call
# (use_c => 0), e.g. to get a clean timing baseline or to avoid the
# compile attempt's overhead/noise in a container known to lack a
# compiler. Everything below is skipped and $HAS_C stays 0.
unless ( $ENV{IF_NO_C} ) {
# Defaults recorded when `perl Makefile.PL` ran. Makefile.PL generates
# Algorithm::Classifier::IsolationForest::BuildFlags, capturing the
# IF_* values active at configure time plus whether a prebuilt object
# was scheduled for install (see "Compile at install time" in the POD
# below). From a plain source checkout the generated file is absent,
# the hard defaults here apply, and no prebuilt object is looked for.
my ( $def_opt, $def_arch, $def_no_omp, $prebuilt ) = ( '-O3', '', 0, 0 );
{
local $@;
my $rec = eval {
require Algorithm::Classifier::IsolationForest::BuildFlags;
Algorithm::Classifier::IsolationForest::BuildFlags::flags();
};
if ( ref $rec eq 'HASH' ) {
$def_opt = $rec->{opt} if defined $rec->{opt};
$def_arch = $rec->{arch} if defined $rec->{arch};
$def_no_omp = $rec->{no_openmp} ? 1 : 0;
$prebuilt = $rec->{prebuilt} ? 1 : 0;
}
}
# -O3 is the usual default: it's safe to enable unconditionally and
# matters here -- the extended-mode oblique dot product is wrapped in
# `#pragma omp simd`, but without aggressive optimization the compiler
# may still emit scalar code. Use OPTIMIZE (not CCFLAGS) -- CCFLAGS is
# prepended to the cc line and would be shadowed by Perl's own `-O2 -g`
# that ExtUtils::MakeMaker appends afterward (last `-O` wins in gcc).
# IF_OPT overrides the level itself (e.g. IF_OPT=-O2 to work around a
# miscompile, or to shorten build time while developing); it's
# 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,
LIBS => '-lm -lgomp',
@install,
);
$HAS_C = 1;
};
} ## end unless ($no_omp)
unless ($HAS_C) {
local $@;
eval {
require Inline;
Inline->import(
C => $serial_tag . $C_CODE,
OPTIMIZE => $opt_level,
LIBS => '-lm',
@install,
);
$HAS_C = 1;
};
} ## end unless ($HAS_C)
$C_SOURCE = 'prebuilt' if $HAS_C;
} else {
# Fast path: the object compiled at `make` time was installed
# under auto/ like any XS module, so plain XSLoader digs it out of
# @INC with no Inline involvement -- no compiler, no _Inline/
# directory, and a few ms instead of a first-run compile. Any
# failure (object deleted, different perl, version mismatch after
# an upgrade, libgomp since removed) just falls through to the
# runtime build.
if ($use_prebuilt) {
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
# checkout, IF_RUNTIME_BUILD=1, or IF_* values differing from the
# ones recorded at configure time. Try compiling with OpenMP
# first; on any failure (compiler doesn't accept -fopenmp, libgomp
# missing, etc.) fall back to a serial build.
unless ( $HAS_C or $no_omp ) {
local $@;
eval {
require Inline;
Inline->import(
C => $omp_tag . $C_CODE,
CCFLAGS => '-fopenmp',
OPTIMIZE => $opt_level,
LIBS => '-lm -lgomp',
);
$HAS_C = 1;
$C_SOURCE = 'runtime';
};
} ## end unless ( $HAS_C or $no_omp )
unless ($HAS_C) {
local $@;
eval {
require Inline;
Inline->import(
C => $serial_tag . $C_CODE,
OPTIMIZE => $opt_level,
LIBS => '-lm',
);
$HAS_C = 1;
$C_SOURCE = 'runtime';
};
} ## end unless ($HAS_C)
} ## end else [ if ( $ENV{IF_INSTALL_BUILD} ) ]
$OPT_LEVEL = $opt_level if $HAS_C;
} ## end unless ( $ENV{IF_NO_C} )
$HAS_OPENMP = ( $HAS_C && defined &has_openmp_xs && has_openmp_xs() ) ? 1 : 0;
$HAS_SIMD = ( $HAS_C && defined &has_simd_xs && has_simd_xs() ) ? 1 : 0;
}
=encoding UTF-8
=head1 NAME
Algorithm::Classifier::IsolationForest - unsupervised anomaly detection via Isolation Forest or Extended Isolation Forest
=head1 SYNOPSIS
use Algorithm::Classifier::IsolationForest;
my @data = ([0.1, -0.2], [0.0, 0.1], [5.0, 6.0], ...);
# Classic, axis-parallel Isolation Forest
my $iforest = Algorithm::Classifier::IsolationForest->new(
n_trees => 100,
sample_size => 256,
seed => 42,
);
$iforest->fit(\@data);
my $scores = $iforest->score_samples(\@data); # arrayref, each in (0,1]
my $flags = $iforest->predict(\@data, 0.6); # arrayref of 0/1
# Save and reload
$iforest->save('model.json');
my $reloaded = Algorithm::Classifier::IsolationForest->load('model.json');
# Extended Isolation Forest (oblique hyperplane splits)
my $eif = Algorithm::Classifier::IsolationForest->new(
mode => 'extended',
seed => 42,
);
$eif->fit(\@data);
# Parallel training (fork-based, Unix-like platforms): build the
# n_trees across several worker processes.
my $iforest = Algorithm::Classifier::IsolationForest->new(
n_trees => 200,
sample_size => 256,
seed => 42,
parallel_fit => 4, # 4 forked workers
);
$iforest->fit(\@data);
# Pre-pack a dataset to skip the per-call input-walk cost when the
# same data gets scored many times (interactive tuning, dashboards).
my $packed = $iforest->pack_data(\@data);
my $scores = $iforest->score_samples($packed);
my $flags = $iforest->predict($packed, 0.6);
# Get scores and labels as two flat arrayrefs in one call -- cheaper
# than score_predict_samples when you don't need the paired shape.
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)>.
For data that arrives as a stream and may drift over time, the companion
class L<Algorithm::Classifier::IsolationForest::Online> implements Online
Isolation Forest (Leveni et al. 2024): no C<fit()>, instead points are
learned as they arrive and forgotten once they age out of a sliding
window. Models saved by either class can be loaded through L</load>,
which dispatches on the stored format tag.
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
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
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
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
=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
- 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
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
for my $doc (qw(schema_version schema_description)) {
croak "$doc must be a plain string"
if defined $self->{$doc} && ref $self->{$doc};
}
_validate_feature_descriptions( $self->{feature_names}, $self->{feature_descriptions} )
if defined $self->{feature_descriptions};
# Optional Algorithm::ToNumberMunger integration: a declarative spec
# hash compiled into a plan that turns raw tagged values into numbers.
# Compiled eagerly so every spec error surfaces here rather than at
# first scoring; the module itself is only required when a spec is
# actually given, keeping it an optional dependency.
if ( defined $args{mungers} ) {
croak "mungers must be a hashref of 'tag => munger spec'"
unless ref $args{mungers} eq 'HASH';
croak "mungers requires feature_names (the munger plan compiles against them)"
unless ref $self->{feature_names} eq 'ARRAY' && @{ $self->{feature_names} };
$self->{mungers} = $args{mungers};
$self->{_munger_plan} = _compile_mungers( $self->{feature_names}, $self->{mungers} );
$self->{munger_module_version} = $Algorithm::ToNumberMunger::VERSION;
}
croak "n_trees must be >= 1" unless $self->{n_trees} >= 1;
croak "sample_size must be >= 1" unless $self->{sample_size} >= 1;
croak "extension_level must be >= 0"
if defined $self->{extension_level} && $self->{extension_level} < 0;
croak "contamination must be a number in (0, 0.5]"
if defined $self->{contamination}
&& !( $self->{contamination} > 0 && $self->{contamination} <= 0.5 );
croak "parallel_fit must be a positive integer"
if defined $self->{parallel_fit}
&& ( $self->{parallel_fit} !~ /^\d+$/ || $self->{parallel_fit} < 1 );
return bless $self, $class;
} ## end sub new
=head2 decision_threshold
The score cutoff C<predict> uses by default; undef unless C<contamination> was
set.
=cut
sub decision_threshold { return $_[0]->{threshold} }
=head2 set_voting
Switches the scoring-time aggregation between C<'mean'> and C<'majority'> on an
existing model and returns C<$self> (so it chains). The forest itself is
identical in both modes -- only the way per-tree results are combined changes
-- so this never rebuilds a single tree.
$iforest->set_voting('majority');
$iforest->set_voting('mean', \@training_data);
The one thing that does not carry over is a C<contamination>-learned
L</decision_threshold>. That cutoff is a quantile of whichever per-point
quantity the mode thresholds against -- the averaged anomaly score under
C<'mean'>, the per-tree majority pivot under C<'majority'> -- and those live in
different spaces, so a threshold learned in one mode flags the wrong fraction
in the other. When the model was fitted with C<contamination>, C<set_voting>
therefore relearns the threshold for the target mode, which requires the
original training data to be passed as the second argument (the model does not
retain it). Switching a model that had no C<contamination> needs no data:
C<predict> falls back to C<0.5>, which is meaningful in both modes.
Passing the current mode is a no-op (returns immediately, no data needed).
Calling this before L</fit> just records the mode for the eventual fit.
=cut
sub set_voting {
my ( $self, $voting, $data ) = @_;
croak "set_voting: voting must be 'mean' or 'majority'"
unless defined $voting && $voting =~ /\A(?:mean|majority)\z/;
return $self if $self->{voting} eq $voting;
# A learned threshold only exists once a contamination-fitted model has
# been fit(); that value is mode-specific and must be relearned against
# the same training set (see _learn_contamination_threshold). Everything
# else -- pre-fit models, and fitted models without contamination -- just
# flips the knob; predict()'s 0.5 fallback is valid in either mode.
my $fitted = ref $self->{trees} eq 'ARRAY' && @{ $self->{trees} };
my $recalibrate = $fitted && defined $self->{contamination};
if ($recalibrate) {
croak "set_voting: switching a contamination-fitted model requires "
. "the original training data as the second argument to "
. "recalibrate the decision threshold"
unless ref $data eq 'ARRAY' && @$data;
}
$self->{voting} = $voting;
$self->_learn_contamination_threshold($data) if $recalibrate;
return $self;
} ## end sub set_voting
=head2 feature_names
Returns the arrayref of feature name strings stored with the model, or undef
if none were provided at fit time.
my $names = $iforest->feature_names;
=cut
sub feature_names { return $_[0]->{feature_names} }
=head2 schema_version
Returns the user-owned schema version string stored with the model
(usually via a prototype -- see L</PROTOTYPES>), or undef if none was
recorded.
my $sv = $iforest->schema_version;
=cut
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
my $n = scalar @$train;
# The sub-sample cannot be larger than the data set itself.
my $psi = min( $self->{sample_size}, $n );
$self->{c_psi} = _c($psi);
$self->{psi_used} = $psi;
# Resolve the extension level against the data's dimensionality.
if ( $self->{mode} eq 'extended' ) {
my $max_ext = $n_features - 1;
my $ext
= defined $self->{extension_level}
? $self->{extension_level}
: $max_ext;
$ext = 0 if $ext < 0;
$ext = $max_ext if $ext > $max_ext;
$self->{extension_level_used} = $ext;
} else {
$self->{extension_level_used} = undef;
}
# Height limit: the average tree height ceil(log2(psi)). Past this depth the
# remaining points are scored using the c(size) adjustment instead.
my $limit
= defined $self->{max_depth}
? $self->{max_depth}
: ceil( log($psi) / log(2) );
$limit = 1 if $limit < 1;
$self->{max_depth_used} = $limit;
srand( $self->{seed} ) if defined $self->{seed};
my $workers = $self->{parallel_fit};
if ( defined $workers
&& $workers > 1
&& $self->{n_trees} > 1
&& _fork_supported() )
{
$self->{trees} = $self->_fit_trees_parallel( $train, $psi, $limit, $workers );
} elsif ( $self->{_use_c} && $self->{_use_openmp_fit} ) {
$self->{trees} = $self->_build_forest_openmp( $train, $psi, $limit, $self->{n_trees} );
} elsif ( $self->{_use_c} ) {
$self->{trees}
= $self->_build_forest_c( $train, $psi, $limit, $self->{n_trees} );
} else {
my @trees;
for ( 1 .. $self->{n_trees} ) {
my $sample = _subsample( $train, $psi );
push @trees, $self->_build_tree( $sample, 0, $limit );
}
$self->{trees} = \@trees;
}
# On a re-fit, packed scoring buffers from the previous fit are still
# sitting on the object; score_samples() below would pick them up and
# learn the contamination threshold against the OLD forest. Drop them
# so the training-set scoring runs pure-Perl against the trees just
# built; _rebuild_c_trees repacks from the new trees at the end.
delete @$self{qw(_c_nodes _c_coef_idx _c_coef_val)};
# If a contamination rate was requested, learn the score cutoff that flags
# that fraction of the training set. The threshold lands midway inside a
# real gap between flagged and unflagged training scores (ties at the
# k-boundary shift the cut to the nearest gap -- see
# _threshold_from_ranked), so it sits strictly between attainable values:
# unambiguous and robust to the tiny float rounding introduced by JSON
# serialisation.
#
# Under voting => 'majority' the value predict() thresholds against is
# the PER-TREE score, so the quantity to rank is each training point's
# majority pivot -- the per-tree cutoff at which that point loses its
# majority (see _majority_pivot_scores). A point is flagged iff its
# pivot >= threshold, exactly the relation the mean-mode score has, so
# the midpoint selection below serves both modes unchanged.
$self->_learn_contamination_threshold($train)
if defined $self->{contamination};
$self->_rebuild_c_trees() if $self->{_use_c};
return $self;
} ## end sub fit
=head2 fit_tagged(\@rows)
Trains the model on an arrayref of hashrefs of named feature values --
the tagged counterpart of L</fit>. Each row goes through
L</tagged_row_to_array> (and therefore through the munger plan when
C<mungers> is configured, which is the point: training data and scoring
data are munged by the identical plan), then the positional rows are
handed to C<fit>.
$iforest->fit_tagged([
{ cpu => 0.9, mem => 0.4, disk => 0.1 },
{ cpu => 0.2, mem => 0.3, disk => 0.2 },
...
]);
Requires stored C<feature_names>. Croaks under the same conditions as
L</tagged_row_to_array>, naming the offending row by index.
=cut
sub fit_tagged {
my ( $self, $data ) = @_;
croak "fit_tagged() expects a non-empty arrayref of hashref samples"
unless ref $data eq 'ARRAY' && @$data;
my @rows;
for my $i ( 0 .. $#$data ) {
push @rows, $self->tagged_row_to_array( $data->[$i], "fit_tagged (row $i)" );
}
return $self->fit( \@rows );
} ## end sub fit_tagged
=head2 pack_data(\@data)
Returns an opaque, blessed wrapper around the input dataset that the
scoring methods can use directly, skipping the per-call work of walking
the arrayref-of-arrayrefs and converting each cell into a double. At
high feature counts this is a meaningful win when the same dataset is
scored repeatedly (e.g. interactive threshold tuning, dashboards,
plotting that updates as parameters change).
Requires the Inline::C backend; croaks if C<use_c> is false.
my $packed = $forest->pack_data(\@data);
# Now any of these accept either an arrayref or the packed wrapper:
my $scores = $forest->score_samples($packed);
my $flags = $forest->predict($packed, 0.6);
my ($s, $l) = $forest->score_predict_split($packed);
The wrapper has C<n_pts> and C<n_feats> accessors for introspection.
The feature count is matched against the model on every call; passing a
packed dataset built for a different feature count is a fatal error.
=cut
=head2 path_lengths(\@data)
Returns an arrayref of the mean isolation depth per sample, for inspection.
my $lengths = $forest->path_lengths(\@data);
print "x, y, length\n";
my $int=0;
while (defined($data[$int])) {
print $data[$int][0].', '.$data[$int][1].', '.$lengths->[$int]."\n";
$int++;
}
=cut
sub path_lengths {
my ( $self, $data ) = @_;
$self->_check_fitted;
my $trees = $self->{trees};
my $t = scalar @$trees;
if ( $self->{_use_c} && $self->{_c_nodes} ) {
my ( $n_pts, $nf, $x_packed ) = $self->_resolve_input($data);
my $sums_packed = "\0" x ( $n_pts * 8 );
score_all_xs(
$self->{_c_nodes}, $self->{_c_coef_idx}, $self->{_c_coef_val},
$x_packed, $sums_packed, $n_pts,
$nf, $t, $self->{_use_openmp}
);
my $result = [];
finalize_path_lengths_xs( $sums_packed, $n_pts, $t + 0.0, $result );
return $result;
} ## end if ( $self->{_use_c} && $self->{_c_nodes} )
$data = $self->_prepare_perl_input($data);
my $nan = $self->{missing} eq 'nan' ? 1 : 0;
# Pure-Perl fallback (tree-outer, sample-inner for cache locality).
my @sums = (0) x @$data;
for my $tree (@$trees) {
for my $i ( 0 .. $#$data ) {
$sums[$i] += _path_length( $data->[$i], $tree, 0, $nan );
}
}
return [ map { $_ / $t } @sums ];
} ## end sub path_lengths
=head2 predict(\@data, $threshold)
Returns an arrayref of 0/1 labels for the specified data.
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
vote_all_xs( $self->{_c_nodes}, $self->{_c_coef_idx}, $self->{_c_coef_val},
$x_packed, $votes_packed, $n_pts, $nf, $t, $cut, 0, $self->{_use_openmp} );
# votes/t is exactly the "divide the sum buffer by t" shape
# finalize_path_lengths_xs implements, so reuse it.
my $result = [];
finalize_path_lengths_xs( $votes_packed, $n_pts, $t + 0.0, $result );
return $result;
} ## end if ( $self->{_use_c} && $self->{_c_nodes} )
my $votes = $self->_vote_counts_perl( $self->_prepare_perl_input($data), $cut );
return [ map { $_ / $t } @$votes ] if _NV_IS_DOUBLE;
# Wide-NV perls: the C finalizers divide in double, so narrow the
# vote fraction to match -- see _NV_IS_DOUBLE.
return [ map { _to_double( $_ / $t ) } @$votes ];
} ## end if ( $self->{voting} eq 'majority' )
if ( $self->{_use_c} && $self->{_c_nodes} ) {
my ( $n_pts, $nf, $x_packed ) = $self->_resolve_input($data);
my $sums_packed = "\0" x ( $n_pts * 8 );
score_all_xs(
$self->{_c_nodes}, $self->{_c_coef_idx}, $self->{_c_coef_val},
$x_packed, $sums_packed, $n_pts,
$nf, $t, $self->{_use_openmp}
);
if ( $c > 0 ) {
my $inv = log(2) / ( $c * $t );
my $result = [];
finalize_scores_xs( $sums_packed, $n_pts, $inv, $result );
return $result;
}
return [ (0.5) x $n_pts ];
} ## end if ( $self->{_use_c} && $self->{_c_nodes} )
$data = $self->_prepare_perl_input($data);
my $nan = $self->{missing} eq 'nan' ? 1 : 0;
# Pure-Perl fallback (tree-outer, sample-inner for cache locality).
my @sums = (0) x @$data;
for my $tree (@$trees) {
for my $i ( 0 .. $#$data ) {
$sums[$i] += _path_length( $data->[$i], $tree, 0, $nan );
}
}
# Precompute the single normalising factor; exp() is a direct FPU
# instruction and faster than Perl's general-purpose 2**x (pow).
# Derivation: 2**(-avg/c) = 2**(-(sum/t)/c) = exp(-sum * log(2)/(c*t))
if ( $c > 0 ) {
my $inv = log(2) / ( $c * $t );
return [ map { exp( -$_ * $inv ) } @sums ];
}
return [ (0.5) x @sums ];
} ## end sub score_samples
=head2 score_sample_tagged(\%row)
Scores a single sample supplied as a hashref of named feature values.
The model must have stored feature names (set via C<feature_names> in
C<new()> or the C<-t> CLI flag at fit time).
Returns a scalar anomaly score in (0, 1].
my $score = $forest->score_sample_tagged({ cpu => 0.9, mem => 0.4 });
Croaks if the model has no stored feature names, if the hashref contains a
key that is not a known feature name, or if a feature name is absent from the
hashref.
=cut
sub score_sample_tagged {
my ( $self, $row ) = @_;
my $vec = $self->tagged_row_to_array( $row, 'score_sample_tagged' );
my $result = $self->score_samples( [$vec] );
return $result->[0];
}
=head2 score_predict_samples
Returns an array ref of arrays. First value of each sub array is the score with the second being
0/1 for if it is a anomaly or not.
C<$threshold> defaults the same way as in C<predict>.
Under C<< voting => 'majority' >> the score is the anomaly vote fraction at
C<$threshold> (used as the per-tree cutoff) and the label is the majority
vote, matching C<score_samples>/C<predict> semantics in that mode.
my $results = $forest->score_predict_samples(\@data, $threshold);
print "x, y, score, result\n";
my $int=0;
while (defined($data[$int])) {
print $data[$int][0].', '.$data[$int][1].', '.$results->[$int][0].', '.$results->[$int][1]."\n";
$int++;
}
=cut
sub score_predict_samples {
my ( $self, $data, $threshold ) = @_;
$threshold
= defined $threshold ? $threshold
: defined $self->{threshold} ? $self->{threshold}
: 0.5;
$self->_check_fitted;
# Majority voting: [vote_fraction, majority_label] pairs. Needs the
# full vote counts (the fraction is part of the return shape), so no
# early exit here -- that saving is predict()-only.
if ( $self->{voting} eq 'majority' ) {
my $trees = $self->{trees};
my $t = scalar @$trees;
my $cut = _depth_cut( $threshold, $self->{c_psi} );
my $maj = _min_votes($t);
if ( $self->{_use_c} && $self->{_c_nodes} ) {
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
L<https://arxiv.org/abs/2505.09593>
L<https://github.com/ineveLoppiliF/Online-Isolation-Forest>
L<https://proceedings.mlr.press/v235/leveni24a.html>
=cut
###
###
### internal stuff below
###
###
#-------------------------------------------------------------------------------
# c(n): the expected path length of an unsuccessful search in a binary search
# tree of n nodes. Isolation Forest uses it (a) to adjust the path length when a
# leaf still holds more than one point (depth limit reached), and (b) to
# normalise the average path length into a 0..1 anomaly score.
#-------------------------------------------------------------------------------
sub _c {
my ($n) = @_;
return 0.0 if $n <= 1;
return 1.0 if $n == 2;
my $harmonic = log( $n - 1 ) + EULER; # H(n-1) ~= ln(n-1) + gamma
return 2.0 * $harmonic - ( 2.0 * ( $n - 1 ) / $n );
}
#-------------------------------------------------------------------------------
# Majority-voting (voting => 'majority') helpers. MVIForest -- Chabchoub,
# Togbe, Boly & Chiky 2022 (see REFERENCES) -- has each tree vote a point
# anomalous when the tree's own score 2**(-h/c(psi)) clears the decision
# threshold, and takes the majority of the votes as the label. Trees are
# untouched; only these scoring-time aggregation helpers differ from the
# classic mean-path-length pipeline.
#-------------------------------------------------------------------------------
# Depth-domain image of the per-tree score cutoff: a tree votes a point
# anomalous when 2**(-h/c) >= theta, i.e. h <= -c * log2(theta). Doing the
# log once here keeps exp/log out of the per-point per-tree loops (both C
# and Perl compare raw path lengths against this cut). Degenerate inputs
# pin the cut so `h <= cut` still behaves: theta <= 0 is cleared by every
# per-tree score (all in (0, 1]), so +inf lets every tree vote; c <= 0 only
# happens for psi <= 1 forests, whose score convention is a flat 0.5 (see
# score_samples), so all trees vote iff theta is at or below that pivot.
sub _depth_cut {
my ( $theta, $c ) = @_;
return ( $theta <= 0.5 ? 9**9**9 : -1.0 ) if $c <= 0;
return 9**9**9 if $theta <= 0;
return -$c * log($theta) / log(2);
}
# Smallest number of per-tree anomaly votes that constitutes a majority:
# int(t/2) + 1, i.e. strictly more than half the trees for both odd and
# even tree counts (the paper's "t/2 + 1").
sub _min_votes { return int( $_[0] / 2 ) + 1 }
#-------------------------------------------------------------------------------
# Contamination threshold selection: given the training scores ranked
# descending and the target flag count k, return a cutoff sitting midway
# inside the gap between the last flagged and the first unflagged score.
#
# Tied scores at the k-boundary make an exact count of k unattainable (the
# tie block can only go one way or the other) AND make the naive midpoint
# degenerate -- it equals the tied value, leaving predict()'s >= comparison
# balanced on exact float equality. Mean-mode scores are continuous enough
# that this practically never happens, but majority-mode pivots are
# structurally quantized (path lengths at the depth cap take few distinct
# values -- see _majority_pivot_scores), so ties there are the norm, and
# the score <-> depth-cut conversion adds an exp/log round trip that needs
# real slack around the cutoff rather than exact-equality behaviour. The
# whole tie block therefore goes to whichever side lands the flag count
# closest to k, preferring the flagging side on a dead heat.
#-------------------------------------------------------------------------------
sub _threshold_from_ranked {
my ( $desc, $k ) = @_;
my $n = scalar @$desc;
return $desc->[ $n - 1 ] - 1e-9 if $k >= $n; # flag everything
my $v = $desc->[ $k - 1 ];
return ( $v + $desc->[$k] ) / 2.0 if $desc->[$k] < $v; # clean gap at k
# Tie block straddling the k-boundary: locate its edges.
my $i = $k - 1;
$i-- while $i > 0 && $desc->[ $i - 1 ] == $v; # first index holding $v
my $j = $k;
$j++ while $j < $n && $desc->[$j] == $v; # first index below $v
if ( $i > 0 && ( $k - $i ) < ( $j - $k ) ) {
# Excluding the block lands closer to k: flag the $i points above it.
return ( $desc->[ $i - 1 ] + $v ) / 2.0;
}
return $j < $n
? ( $v + $desc->[$j] ) / 2.0 # include the block: flag $j
: $desc->[ $n - 1 ] - 1e-9; # block runs to the end
} ## end sub _threshold_from_ranked
# Pure-Perl vote counter: votes[i] = how many trees give point i a path
# length at or under the depth cut. Tree-outer / sample-inner for cache
# locality, mirroring the mean-mode fallback loops. $data must already
# be through _prepare_perl_input.
sub _vote_counts_perl {
my ( $self, $data, $cut ) = @_;
my $trees = $self->{trees};
my $nan = $self->{missing} eq 'nan' ? 1 : 0;
my @votes = (0) x @$data;
for my $tree (@$trees) {
for my $i ( 0 .. $#$data ) {
$votes[$i]++ if _path_length( $data->[$i], $tree, 0, $nan ) <= $cut;
}
}
return \@votes;
} ## end sub _vote_counts_perl
#-------------------------------------------------------------------------------
# Contamination support for majority voting: each training point's majority
# pivot -- the per-tree score threshold at which the point loses its
# majority. A point is flagged at cutoff theta iff at least min_votes of
# its per-tree path lengths h satisfy h <= -c*log2(theta), which holds iff
# its min_votes-th SMALLEST path length h_(maj) does, i.e. iff
# 2**(-h_(maj)/c) >= theta. So the pivot m = 2**(-h_(maj)/c) relates to
# the majority-mode threshold exactly as the mean-mode score relates to
# its threshold, and fit()'s midpoint selection works on either unchanged.
#
# Pure Perl by necessity: the per-tree path lengths never cross the C
# boundary individually (score_all_xs/vote_all_xs only return per-point
# aggregates), and fit() has already dropped any stale packed buffers when
# this runs -- the same situation as mean mode's training-set scoring pass.
#-------------------------------------------------------------------------------
# Learn the contamination cutoff for the CURRENT voting mode from a training
# set. Ranks the per-point quantity the active aggregation thresholds against
# -- the mean-mode anomaly score, or the majority pivot under
# voting => 'majority' -- and lands the cutoff midway inside a real gap between
# flagged and unflagged values (ties at the k-boundary shift it to the nearest
# gap; see _threshold_from_ranked), so it sits strictly between attainable
# values: unambiguous and robust to the float rounding JSON introduces. A
# point is flagged iff its statistic >= threshold in either mode, so the
# midpoint selection serves both unchanged. Shared by fit() (which passes the
# prepared training set after dropping any stale packed buffers) and
# set_voting() (which passes the caller-supplied training set against the
# live, fully packed forest); $data may hold raw undef cells either way, since
# the scorers below densify from missing_fill.
sub _learn_contamination_threshold {
my ( $self, $data ) = @_;
my $scores
= $self->{voting} eq 'majority'
? $self->_majority_pivot_scores($data)
: $self->score_samples($data);
my @desc = sort { $b <=> $a } @$scores;
my $n_pts = scalar @desc;
my $k = int( $self->{contamination} * $n_pts + 0.5 );
$k = 1 if $k < 1;
$k = $n_pts if $k > $n_pts;
$self->{threshold} = _threshold_from_ranked( \@desc, $k );
return;
} ## end sub _learn_contamination_threshold
sub _majority_pivot_scores {
my ( $self, $data ) = @_;
my $trees = $self->{trees};
my $t = scalar @$trees;
my $c = $self->{c_psi};
my $maj = _min_votes($t);
my $rows = $self->_prepare_perl_input($data);
my $nan = $self->{missing} eq 'nan' ? 1 : 0;
# psi <= 1 degenerate forest: every per-tree score is pinned at 0.5
# (matching score_samples' convention), so every pivot is too.
return [ (0.5) x @$rows ] unless $c > 0;
my $inv = log(2) / $c;
my @pivots;
for my $x (@$rows) {
my @paths = sort { $a <=> $b } map { _path_length( $x, $_, 0, $nan ) } @$trees;
push @pivots, exp( -$paths[ $maj - 1 ] * $inv );
}
return \@pivots;
} ## end sub _majority_pivot_scores
# One draw from the standard normal N(0,1) via Box-Muller. Used to pick the
# random hyperplane orientations in Extended Isolation Forest mode.
sub _randn {
my $u1 = rand() || 1e-12;
my $u2 = rand();
return sqrt( -2.0 * log($u1) ) * cos( TWO_PI * $u2 ) if _NV_IS_DOUBLE;
# Wide-NV perls: round after every operation _c_randn() performs in
# double, so both backends draw the same coefficient bit patterns
# (up to libm's own double-vs-long-double disagreements on rare
# rounding ties).
my $s = _to_double( sqrt( -2.0 * _to_double( log($u1) ) ) );
my $c = _to_double( cos( _to_double( TWO_PI * $u2 ) ) );
return _to_double( $s * $c );
} ## end sub _randn
#-------------------------------------------------------------------------------
# Draw $k samples without replacement via a partial Fisher-Yates shuffle of the
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
my ( @coef, $b );
$b = 0.0;
for my $f (@idx) {
my $c = _randn();
if (_NV_IS_DOUBLE) {
my $p = $lo->[$f] + rand() * ( $hi->[$f] - $lo->[$f] ); # point in the box
push @coef, $c;
$b += $c * $p;
} else {
# Round each op to double in the same order as the C builder's
# p = lo + rand() * (hi - lo); b += c * p;
# -- see _NV_IS_DOUBLE.
my $p = _to_double( rand() * _to_double( $hi->[$f] - $lo->[$f] ) );
$p = _to_double( $lo->[$f] + $p );
push @coef, $c;
$b = _to_double( $b + _to_double( $c * $p ) );
}
} ## end for my $f (@idx)
# A point missing any feature on the hyperplane (nan mode only) routes
# to the right child: in the C scorer the dot product becomes NaN and
# (NaN <= b) is false, so this keeps fit and score consistent. Under
# die/zero/impute every cell is defined, so the per-feature "defined"
# check and early-exit are dead weight there and skipped entirely.
my ( @left, @right );
if ($nan) {
for my $row (@$X) {
my $dot = 0.0;
my $missing = 0;
for ( 0 .. $#idx ) {
my $v = $row->[ $idx[$_] ];
if ( !defined $v ) { $missing = 1; last }
$dot += $coef[$_] * $v;
}
if ( !$missing && $dot <= $b ) { push @left, $row }
else { push @right, $row }
} ## end for my $row (@$X)
} else {
for my $row (@$X) {
my $dot = 0.0;
$dot += $coef[$_] * $row->[ $idx[$_] ] for 0 .. $#idx;
if ( $dot <= $b ) { push @left, $row }
else { push @right, $row }
}
}
return [ _NODE_OBLIQUE, \@idx, \@coef, $b, \@left, \@right ];
} ## end sub _oblique_split
#-------------------------------------------------------------------------------
# Path length of a single point in a single tree: edges traversed until a leaf,
# plus c(leaf size) when the leaf still holds several points.
#
# Node layout (arrayref, slot 0 = type):
# _NODE_LEAF [0, size]
# _NODE_AXIS [1, attr, split, left, right]
# _NODE_OBLIQUE [2, \@idx, \@coef, b, left, right]
#
# The type tag is also used as a loop sentinel: 0 (_NODE_LEAF) is falsy.
# No $self argument -- the node type encodes everything needed.
#-------------------------------------------------------------------------------
# The optional $nan flag selects the nan-strategy routing: a point missing
# the split feature goes to the right child (matching the C scorer, where
# the NaN comparison is false). Without it, undef is coerced to 0 -- the
# behaviour the die/zero/impute strategies rely on (their data is dense by
# the time it reaches here, so the "// 0" is normally a no-op).
sub _path_length {
my ( $x, $node, $depth, $nan ) = @_;
while ( $node->[0] ) { # false only for leaf (type 0)
if ( $node->[0] == _NODE_AXIS ) { # [1, attr, split, left, right]
if ($nan) {
my $v = $x->[ $node->[1] ];
$node = ( defined($v) && $v < $node->[2] ) ? $node->[3] : $node->[4];
} else {
$node = ( $x->[ $node->[1] ] // 0 ) < $node->[2] ? $node->[3] : $node->[4];
}
} else { # [2, \@idx, \@coef, b, left, right]
my ( $idx, $coef, $b ) = ( $node->[1], $node->[2], $node->[3] );
if ($nan) {
my $dot = 0.0;
my $missing = 0;
for ( 0 .. $#$idx ) {
my $v = $x->[ $idx->[$_] ];
if ( !defined $v ) { $missing = 1; last }
$dot += $coef->[$_] * $v;
}
$node = ( !$missing && $dot <= $b ) ? $node->[4] : $node->[5];
} else {
my $dot = 0.0;
$dot += $coef->[$_] * ( $x->[ $idx->[$_] ] // 0 ) for 0 .. $#$idx;
$node = $dot <= $b ? $node->[4] : $node->[5];
}
} ## end else [ if ( $node->[0] == _NODE_AXIS ) ]
$depth++;
} ## end while ( $node->[0] )
return $depth + _c( $node->[1] ); # leaf size at slot 1
} ## end sub _path_length
# Recursively convert a version-0 hash-based tree node to the version-1
# array format. Called by from_json when loading an old saved model.
sub _hash_node_to_array {
my ($node) = @_;
if ( $node->{leaf} ) {
return [ _NODE_LEAF, $node->{size} ];
} elsif ( exists $node->{attr} ) {
return [
_NODE_AXIS, $node->{attr},
$node->{split}, _hash_node_to_array( $node->{left} ),
_hash_node_to_array( $node->{right} ),
];
} else {
return [
_NODE_OBLIQUE, $node->{idx}, $node->{coef}, $node->{b},
_hash_node_to_array( $node->{left} ),
_hash_node_to_array( $node->{right} ),
];
}
} ## end sub _hash_node_to_array
# ---------------------------------------------------------------------------
# _pack_tree($root) -- flatten one tree into three packed buffers.
#
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
} 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 ] ],
[ @{$val}[ $coff .. $coff + $num - 1 ] ],
$b,
_unpack_node( $nodes, $idx, $val, int($li) ),
_unpack_node( $nodes, $idx, $val, int($ri) ),
];
} ## end else [ if ( $type == 0 ) ]
} ## end sub _unpack_node
# Unpacks every tree in the three per-tree packed-buffer arrayrefs
# build_forest_openmp_xs returns into the ordinary nested tree shape.
# The C builder pushes nodes post-order (a node is recorded after both
# of its children), so each tree's root is the LAST node record, not
# index 0 as in _pack_tree's pre-order layout.
sub _unpack_forest {
my ( $nodes_list, $idx_list, $val_list ) = @_;
my @trees;
for my $i ( 0 .. $#$nodes_list ) {
my @nodes = unpack( 'd*', $nodes_list->[$i] );
my @idx = unpack( 'l*', $idx_list->[$i] );
my @val = unpack( 'd*', $val_list->[$i] );
my $root = @nodes / 6 - 1;
push @trees, _unpack_node( \@nodes, \@idx, \@val, $root );
}
return \@trees;
} ## end sub _unpack_forest
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
if ( $m eq 'die' ) {
for my $i ( 0 .. $#$data ) {
my $row = $data->[$i];
for my $f ( 0 .. $nf - 1 ) {
next if defined $row->[$f];
croak "fit(): undef feature value at sample $i, column $f; "
. "construct with missing => 'zero', 'impute', or 'nan' "
. "to train on data with missing values";
}
}
return $data;
} ## end if ( $m eq 'die' )
# nan: leave undef in place -- _build_tree / the split routers handle it.
return $data if $m eq 'nan';
# zero / impute: undef has to become a real number somewhere before a
# split can look at it. The fill vector is computed either way (it's
# needed for persistence and for scoring later), but densifying $data
# into a second, fully separate Perl array here is only necessary for
# the pure-Perl tree builder (_build_tree assumes every cell is
# defined once missing != 'nan' -- see its lo/hi scan). The C
# tree-building path -- _build_forest_c/_build_forest_openmp, and
# every parallel_fit worker, all of which go through pack_input_xs --
# already fills undef cells itself from this same fill vector, so
# skip the redundant whole-dataset copy when that's the path fit()
# will actually take. Scoring the training set for a learned
# contamination threshold (below, in fit()) is unaffected: it always
# runs through the pure-Perl scorer regardless of use_c (fit() drops
# any previous fit's packed buffers before that scoring, and
# _rebuild_c_trees runs after), and that path already tolerates raw
# undef cells
# for both zero (_path_length's "// 0") and impute (_prepare_perl_input
# densifies on demand from missing_fill).
my $fill
= $m eq 'impute'
? $self->_compute_impute_fill($data)
: [ (0) x $nf ];
$self->{missing_fill} = $fill if $m eq 'impute';
delete $self->{_fill_packed};
return $data if $self->{_use_c};
return _densify( $data, $fill );
} ## end sub _prepare_fit_data
# Per-feature fill value (mean or median of the present values) for impute
# mode. Croaks if a feature has no present value to learn from.
sub _compute_impute_fill {
my ( $self, $data ) = @_;
my $nf = $self->{n_features};
my $how = $self->{impute_with};
# C fast path: walks the raw data directly and finds the median via
# quickselect (O(n) average) instead of the Perl fallback's full sort
# (O(n log n)). Produces the same fill values either way -- see
# impute_fill_xs's file-top comment -- so use_c only changes speed
# here, matching the rest of the module.
if ( $self->{_use_c} ) {
my $n = scalar @$data;
my $how_flag = $how eq 'median' ? 1 : 0;
my $fill = [];
impute_fill_xs( $data, $n, $nf, $how_flag, $fill );
return $fill;
}
my @fill;
for my $f ( 0 .. $nf - 1 ) {
my @vals = grep { defined } map { $_->[$f] } @$data;
croak "impute: feature column $f has no present values"
unless @vals;
if ( $how eq 'median' ) {
my @s = sort { $a <=> $b } @vals;
my $k = scalar @s;
$fill[$f]
= $k % 2
? $s[ int( $k / 2 ) ]
: ( $s[ $k / 2 - 1 ] + $s[ $k / 2 ] ) / 2.0;
} else { # mean
my $sum = 0;
if (_NV_IS_DOUBLE) {
$sum += $_ for @vals;
} else {
# impute_fill_xs accumulates the sum in double (over
# SvNV-truncated cells); match its rounding step for step.
$sum = _to_double( $sum + _to_double($_) ) for @vals;
}
$fill[$f] = $sum / scalar @vals;
} ## end else [ if ( $how eq 'median' ) ]
# The fill crosses into the C backend as a double (pack 'd' /
# SvNV), so on wide-NV perls store it already narrowed and both
# builders densify with the identical value.
$fill[$f] = _to_double( $fill[$f] ) unless _NV_IS_DOUBLE;
} ## end for my $f ( 0 .. $nf - 1 )
return \@fill;
} ## end sub _compute_impute_fill
# Return a dense copy of $data with every undef cell replaced by the
# matching per-feature fill value. Leaves present cells untouched.
sub _densify {
my ( $data, $fill ) = @_;
my $nf = scalar @$fill;
return [
map {
my $r = $_;
[ map { defined $r->[$_] ? $r->[$_] : $fill->[$_] } 0 .. $nf - 1 ]
} @$data
];
} ## end sub _densify
# (miss_mode, fill_packed) pair for pack_input_xs, per the active strategy.
# die/zero -> 0 (undef becomes 0.0); impute -> 1 (undef becomes fill[k]);
# nan -> 2 (undef becomes NaN, which the C scorer routes right).
sub _pack_args {
my ($self) = @_;
my $m = $self->{missing};
return ( 2, '' ) if $m eq 'nan';
if ( $m eq 'impute' ) {
my $fill = $self->{missing_fill};
croak "impute model is missing its fill vector"
unless ref $fill eq 'ARRAY' && @$fill == $self->{n_features};
$self->{_fill_packed} //= pack( 'd*', @$fill );
return ( 1, $self->{_fill_packed} );
}
return ( 0, '' ); # die, zero
} ## end sub _pack_args
# Pure-Perl fallback input prep: arrayref-ify, then fill for impute so the
# tree walk sees dense rows. zero/die rely on _path_length's "// 0"; nan
# keeps undef in place for _path_length to route. Returns the rows; the
# caller passes the nan flag to _path_length separately.
sub _prepare_perl_input {
my ( $self, $data ) = @_;
my $rows = $self->_to_arrayref($data);
if ( $self->{missing} eq 'impute' ) {
croak "impute model is missing its fill vector"
unless ref $self->{missing_fill} eq 'ARRAY';
$rows = _densify( $rows, $self->{missing_fill} );
}
return $rows;
} ## end sub _prepare_perl_input
# Minimal PackedData package: opaque token returned by pack_data().
# Exposes n_pts and n_feats accessors for users who want to introspect.
{
package Algorithm::Classifier::IsolationForest::PackedData;
sub n_pts { $_[0]->{n_pts} }
sub n_feats { $_[0]->{n_feats} }
}
1;
( run in 1.024 second using v1.01-cache-2.11-cpan-600a1bdf6e4 )