Algorithm-Classifier-IsolationForest
view release on metacpan or search on metacpan
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
use Config ();
use List::Util qw(min);
use Scalar::Util qw(looks_like_number);
use POSIX qw(ceil);
use JSON::PP ();
use File::Slurp qw(read_file write_file);
our $VERSION = '0.7.0';
use constant EULER => 0.5772156649015329;
# Narrowed to C double precision so _randn() multiplies by the exact
# constant _c_randn() uses. A no-op on nvsize == 8 perls.
use constant TWO_PI => unpack( 'd', pack 'd', 6.283185307179586 );
# Node-type tags stored in index 0 of every tree node arrayref.
# 0 is falsy, so while ($node->[0]) acts as while (!leaf).
use constant _NODE_LEAF => 0;
use constant _NODE_AXIS => 1;
use constant _NODE_OBLIQUE => 2;
# The Inline::C tree builder computes everything in C doubles. On a perl
# whose NV is wider than a double (-Duselongdouble / -Dusequadmath) the
# pure-Perl builder keeps extra low bits at every step, so the two
# backends would stop producing bit-identical trees for the same seed
# (the parity t/03-fit-determinism.t checks). _NV_IS_DOUBLE guards
# narrowing statements wherever the pure-Perl builder computes a value
# that gets STORED in a tree (split points, hyperplane coefficients and
# offsets, impute fills), rounding at the same points the C builder
# rounds. It is compile-time true on nvsize == 8 perls, so there the
# guarded statements are optimised away and cost nothing.
#
# The row-partition loops (v < split, dot <= b) are deliberately NOT
# narrowed: with both operands already double-exact an axis comparison
# is identical anyway, and an oblique dot product accumulated in a wider
# NV flips a comparison only when |dot - b| falls inside the NV-vs-double
# rounding gap (~1e-19 relative) -- negligible, and those are the hot
# loops.
use constant _NV_IS_DOUBLE => $Config::Config{nvsize} == 8;
# Round an NV to C double precision. Only ever reached on wide-NV
# perls -- see _NV_IS_DOUBLE.
#
# Args:
# $nv :: any number. On a -Duselongdouble / -Dusequadmath perl it may
# carry more precision than a C double can hold.
#
# Returns: the same number rounded to the nearest IEEE 754 double. A no-op
# for a value that is already double-exact.
#
# Example:
# _to_double( $lo + rand() * ( $hi - $lo ) ); # what the C builder stores
sub _to_double { unpack 'd', pack 'd', $_[0] }
# ---------------------------------------------------------------------------
# Optional Inline::C accelerator for the scoring hot path.
#
# pack_input_xs(data_sv, out_sv, n_pts, n_feats, miss_mode, fill_sv)
# Walks the Perl arrayref-of-arrayrefs and writes a packed double buffer
# 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)
* fill_sv is only dereferenced when miss_mode == 1. */
void pack_input_xs(SV* data_sv, SV* out_sv, int n_pts, int n_feats,
int miss_mode, SV* fill_sv){
STRLEN tl;
double* out;
const double* fill = NULL;
double missval;
AV* outer;
int i, k;
if (!SvROK(data_sv) || SvTYPE(SvRV(data_sv)) != SVt_PVAV) {
croak("pack_input_xs: data must be an arrayref");
}
outer = (AV*)SvRV(data_sv);
out = (double*)SvPVbyte_force(out_sv, tl);
if (miss_mode == 1) {
STRLEN fl;
fill = (const double*)SvPVbyte(fill_sv, fl);
}
missval = (miss_mode == 2) ? NAN : 0.0;
for (i = 0; i < n_pts; i++) {
SV** row_pp = av_fetch(outer, i, 0);
double* dst = out + (size_t)i * (size_t)n_feats;
if (!row_pp || !*row_pp || !SvROK(*row_pp) ||
SvTYPE(SvRV(*row_pp)) != SVt_PVAV) {
for (k = 0; k < n_feats; k++)
dst[k] = (miss_mode == 1) ? fill[k] : missval;
continue;
}
{
AV* row = (AV*)SvRV(*row_pp);
for (k = 0; k < n_feats; k++) {
SV** v = av_fetch(row, k, 0);
if (v && *v && SvOK(*v)) {
dst[k] = SvNV(*v);
} else {
dst[k] = (miss_mode == 1) ? fill[k] : missval;
}
}
}
}
}
/* first_missing_xs(data_sv, n_pts, n_feats, out_rv)
*
* C replacement for _prepare_fit_data's missing => 'die' scan, which is a
* pure-Perl double loop over every cell of the training set and, being the
* default strategy, the single largest fixed cost of a fit on wide or long
* data. Walks the same arrayref-of-arrayrefs pack_input_xs walks, in the
* same row-major order, and stops at the first cell that is undef or
* absent, pushing (row, col) into out_rv. Leaves out_rv empty when every
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
}
sm = (const double*)SvPVbyte(sm_sv, tl);
out = (AV*)SvRV(out_rv);
av_clear(out);
if (n_pts > 0) av_extend(out, n_pts - 1);
for (i = 0; i < n_pts; i++) {
AV* row = newAV();
av_extend(row, 1);
AvARRAY(row)[0] = newSVnv(sm[i] / t);
AvARRAY(row)[1] = newSViv(sm[i] >= min_votes ? 1 : 0);
AvFILLp(row) = 1;
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;
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
dense_ok = (n_features > 0 && num == n_features);
for (k = 0; dense_ok && k < num; k++) {
if (s->ix[k] < 0 || s->ix[k] >= n_features) dense_ok = 0;
}
if (dense_ok) {
memset(s->dense, 0, (size_t)n_features * sizeof(double));
for (k = 0; k < num; k++) s->dense[s->ix[k]] = s->cv[k];
coff = tb_push_coef(b, s->order, s->dense, n_features);
} else {
coff = tb_push_coef(b, s->ix, s->cv, num);
}
li = _pack_node_xs(aTHX_ slots[4], n_features, b, s);
ri = _pack_node_xs(aTHX_ slots[5], n_features, b, s);
rec = b->nodes + (size_t)my_idx * IF_NZ;
rec[0] = 2.0; rec[1] = (double)coff; rec[2] = (double)num;
rec[3] = (double)li; rec[4] = (double)ri; rec[5] = bcoef;
}
return my_idx;
}
/* n_features may be -1 ("unknown"), which just disables the dense pack.
* The three output SVs are overwritten with the packed bytes. */
void pack_tree_xs(SV* root_rv, int n_features, SV* nodes_sv, SV* idx_sv,
SV* val_sv) {
dTHX;
TreeBuf b;
PackScratch s;
tb_init(&b);
ps_init(&s);
_pack_node_xs(aTHX_ root_rv, n_features, &b, &s);
ps_free(&s);
/* An axis-only tree never calls tb_push_coef, so idx/val stay NULL --
* pass "" so the Perl side always receives a defined string, matching
* what build_forest_openmp_xs hands back. */
sv_setpvn(nodes_sv, (char*)b.nodes, b.n_nodes * IF_NZ * sizeof(double));
sv_setpvn(idx_sv, b.n_idx ? (char*)b.idx : "", b.n_idx * sizeof(int));
sv_setpvn(val_sv, b.n_val ? (char*)b.val : "", b.n_val * sizeof(double));
tb_free(&b);
}
/* ---------------------------------------------------------------------
* 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;
double pivot;
int i, j;
if (a[mid] < a[lo]) _dswap(&a[lo], &a[mid]);
if (a[hi] < a[lo]) _dswap(&a[lo], &a[hi]);
if (a[hi] < a[mid]) _dswap(&a[mid], &a[hi]);
_dswap(&a[mid], &a[hi]);
pivot = a[hi];
i = lo;
for (j = lo; j < hi; j++) {
if (a[j] < pivot) { _dswap(&a[i], &a[j]); i++; }
}
_dswap(&a[i], &a[hi]);
return i;
}
/* Quickselect: returns the k-th smallest (0-indexed) of a[0..n-1],
* reordering a[] in the process (fine -- it's a private scratch copy).
* O(n) average case vs. a full O(n log n) sort. */
static double _kth_smallest(double *a, int n, int k) {
int lo = 0, hi = n - 1;
while (lo < hi) {
int p = _partition_lomuto(a, lo, hi);
if (p == k) return a[p];
if (p < k) lo = p + 1; else hi = p - 1;
}
return a[lo];
}
/* Median of a[0..n-1] (reorders a[]). Odd n: the single middle order
* statistic. Even n: quickselect finds the lower-median at k = n/2-1,
* which leaves every a[i > k] >= a[k] (the standard quickselect
* post-condition) -- so the upper-median is just the min of that
* remaining slice, one more linear scan instead of a second full
* selection pass. */
static double _median_select(double *a, int n) {
if (n % 2 == 1) {
return _kth_smallest(a, n, n / 2);
} else {
int k = n / 2 - 1;
double lower = _kth_smallest(a, n, k);
double upper = a[k + 1];
int i;
for (i = k + 2; i < n; i++) {
if (a[i] < upper) upper = a[i];
}
return (lower + upper) / 2.0;
}
}
void impute_fill_xs(SV* data_sv, int n_pts, int n_feats, int how,
SV* out_rv) {
dTHX;
AV *outer, *out;
DVec *cols;
int i, f;
if (!SvROK(data_sv) || SvTYPE(SvRV(data_sv)) != SVt_PVAV) {
croak("impute_fill_xs: data must be an arrayref");
}
if (!SvROK(out_rv) || SvTYPE(SvRV(out_rv)) != SVt_PVAV) {
croak("impute_fill_xs: out must be an arrayref");
}
outer = (AV*)SvRV(data_sv);
out = (AV*)SvRV(out_rv);
cols = (DVec*)calloc((size_t)n_feats, sizeof(DVec));
for (i = 0; i < n_pts; i++) {
SV** row_pp = av_fetch(outer, i, 0);
AV* row;
if (!row_pp || !*row_pp || !SvROK(*row_pp) ||
SvTYPE(SvRV(*row_pp)) != SVt_PVAV) {
continue;
}
row = (AV*)SvRV(*row_pp);
for (f = 0; f < n_feats; f++) {
SV** v = av_fetch(row, f, 0);
if (v && *v && SvOK(*v)) {
dvec_push(&cols[f], SvNV(*v));
}
}
}
/* Validate every column before freeing anything: croak() longjmps
* out of this function, so any cleanup loop reachable after a
* partial computation has already started (and already freed some
* cols[i].v) risks a double free on those same pointers. Checking
* all columns up front, before the computation loop below frees
* anything, avoids that entirely. Matches the Perl fallback's
* behaviour of reporting the first empty column in feature order. */
for (f = 0; f < n_feats; f++) {
if (cols[f].n == 0) {
int col = f;
for (i = 0; i < n_feats; i++) free(cols[i].v);
free(cols);
croak("impute: feature column %d has no present values", col);
}
}
av_clear(out);
if (n_feats > 0) av_extend(out, n_feats - 1);
for (f = 0; f < n_feats; f++) {
double result;
if (how == 0) {
double sum = 0.0;
for (i = 0; i < (int)cols[f].n; i++) sum += cols[f].v[i];
result = sum / (double)cols[f].n;
} else {
result = _median_select(cols[f].v, (int)cols[f].n);
}
av_store(out, f, newSVnv(result));
free(cols[f].v);
}
free(cols);
}
/* ---------------------------------------------------------------------
* Online Isolation Forest (Algorithm::Classifier::IsolationForest::
* Online) learn / unlearn / score-row accelerators.
*
* Unlike everything above, these operate directly on LIVE Perl
* arrayref trees: online trees mutate on every learned point, so there
* is no immutable packed form to walk during learning. Node layout
* (Online.pm's _N_* constants):
*
* leaf: [0, count, \@lo, \@hi]
* internal: [1, count, \@lo, \@hi, attr, split, left, right]
*
* and a tree record is a hashref { root, count, depth_limit }; root is
* undef until the tree learns its first point, and a leaf built from an
* empty synthetic partition has undef lo/hi until a real point reaches
* it.
*
* Random draws go through Drand01() in EXACTLY the order the pure-Perl
* learn path calls rand(): an optional per-tree subsample gate (drawn
* only when subsample < 1), and -- only when a leaf splits --
* count * nf box-sample draws that SKIP zero-width features (the Perl
* _sample_box never draws for those), then per synthetic internal node
* one draw for the split feature and one for the split value, recursing
* left before right. A learn() with a given seed therefore produces
* BIT-IDENTICAL trees whether use_c is on or off (on nvsize == 8 perls;
* wide-NV perls keep extra low bits in the pure-Perl path, as with
* fit()), which is what lets the online class reuse the existing use_c
* knob for learning instead of growing a new one.
* ------------------------------------------------------------------ */
#define OL_TYPE 0
#define OL_COUNT 1
#define OL_LO 2
#define OL_HI 3
#define OL_ATTR 4
#define OL_SPLIT 5
#define OL_LEFT 6
#define OL_RIGHT 7
/* ln(4) as the exact double Perl's compile-time log(4) produces --
* spelled as a literal (like TWO_PI on the Perl side) so a compiler
* that constant-folds log(4.0) differently from libm cannot introduce
* a one-ulp parity break in the depth budget. */
#define OL_LOG4 1.3862943611198906
/* Depth budget for n points -- the C image of Online.pm's _rpl(). */
static double _ol_rpl(double n, int eta) {
if (n < (double)eta) return 0.0;
return log(n / (double)eta) / OL_LOG4;
}
/* Points a node at `depth` needs before it may split (below which, on
* forgetting, it collapses back into a leaf) -- _split_threshold().
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
Makefile rule that performs the install-time compile. Not meant for
manual use.
=back
If the prebuilt object cannot be loaded for any reason (deleted, built
against a different perl, version mismatch after an upgrade), the
module quietly falls through the same chain as always: runtime
Inline::C build first, pure Perl last.
=head2 Tuning the C build
These environment variables are read once, the first time the module is
loaded, so they must be set before that -- e.g. in the shell before
running a script, not via C<%ENV> inside the script itself. They are
also read by C<perl Makefile.PL> to pick the flags baked into the
prebuilt object (see above); at run time they override the recorded
configure-time values, at the price of a runtime compile.
=over 4
=item * C<IF_NO_C=1> -- skip attempting to build the C backend entirely.
Equivalent to constructing every instance with C<use_c =E<gt> 0>, but
without needing to touch every call site; useful for a clean pure-Perl
timing baseline, or to avoid the compile attempt's overhead/noise on a
host known to lack a C compiler (the attempt already fails gracefully
without this, so it's a convenience, not a correctness fix).
=item * C<IF_OPT=-O2> (or C<-O0>/C<-O1>/C<-Os>/C<-Og>/C<-Oz>) -- override
the default C<-O3>, e.g. to shorten build time while iterating, or work
around a miscompile on an unusual toolchain. Invalid values are ignored
with a warning rather than passed through, since this string reaches a
compiler command line.
=item * C<IF_ARCH=E<lt>valueE<gt>> -- adds C<-march=E<lt>valueE<gt>> so the
compiler can target specific instruction-set extensions (AVX2 gather +
FMA, etc.) for the extended-mode oblique dot product and the fit-time
min/max scan's C<#pragma omp simd> loops. Accepts values like
C<x86-64-v3>, C<skylake>, or C<znver3> -- whatever your compiler's
C<-march=> accepts. Also validated (a restricted character set, not
passed through as-is) for the same reason as C<IF_OPT>. The special
value C<none> (or an empty string) opts out of any arch recorded at
configure time, yielding a plain build. Whenever a C<-march> is in
effect the build also adds C<-ffp-contract=off>: with FMA available
the compiler would otherwise contract C<a*b+c> into fused
multiply-adds whose different rounding breaks the guarantee that
C<use_c =E<gt> 1> and C<use_c =E<gt> 0> build bit-identical trees (the
C<-march> speedup comes from vectorization, not contraction, so this
costs essentially nothing).
=item * C<IF_NATIVE=1> -- shorthand for C<IF_ARCH=native>; ignored if
C<IF_ARCH> is also set. Prefer a specific C<IF_ARCH> value over this on
a machine you don't control exclusively (a shared build host, a
container base image): blanket C<-march=native> pulls in whatever
instruction sets the build host happens to have, including AVX-512 on
some Intel CPUs -- which is known to trigger clock throttling under
sustained heavy use and can make throughput I<worse> than a
conservative target like C<x86-64-v3> (AVX2, no AVX-512). If in doubt,
benchmark both before committing to one.
=item * C<IF_NO_OPENMP=1> -- build (or select) the serial C backend: the
OpenMP compile attempt is skipped entirely, so the resulting object has
no libgomp linkage and never starts an OpenMP runtime inside the
process. This differs from C<OMP_NUM_THREADS=1>, which merely runs the
parallel code on one thread but still loads libgomp. Set at
C<perl Makefile.PL> time it yields a serial prebuilt object; set at run
time against an OpenMP prebuilt install it triggers a runtime serial
build (needing a compiler). An explicit C<IF_NO_OPENMP=0> re-enables
OpenMP over a serial configure-time default.
=back
Whichever of these are used, the cached artefact under C<_Inline/> is
pinned to that build's instruction set -- delete C<_Inline/> (or use a
separate one per host) if the directory is shared across machines with
different CPUs, or a stale binary built for a narrower instruction set
than the current host will simply keep being reused.
=head2 Tuning the OpenMP runtime
These are standard OpenMP environment variables libgomp already reads
at run time (set before running your script, no module-specific
handling needed) -- listed here because they matter most for exactly
the workloads this module has: C<score_all_xs>'s per-point parallel
loop and C<use_openmp_fit>'s per-tree parallel loop.
=over 4
=item * C<OMP_NUM_THREADS=N> -- caps how many threads a parallel region
uses. Useful to leave headroom for other work sharing the machine, or
to pin down C<use_openmp_fit> reproducibility checks (see its docs
above: results don't depend on this, but it's a natural thing to vary
when confirming that).
=item * C<OMP_PROC_BIND=close> / C<OMP_PLACES=cores> -- on multi-socket
or otherwise NUMA machines, pins each thread to a core near where its
data already lives instead of letting the OS scheduler migrate threads
across sockets mid-run. Both C<score_all_xs> (each thread scans its own
slice of the packed query buffer) and C<use_openmp_fit> (each thread
builds one tree from packed training data) benefit from this when the
input is large enough to not fit comfortably in one socket's cache.
=back
These cost nothing to try -- unlike C<IF_ARCH>/C<IF_NATIVE>, they're
read fresh every run, not baked into a cached binary, so there's no
downside to experimenting per invocation.
=head1 GENERAL METHODS
=head2 new(%args)
Inits the object.
- n_trees :: number of isolation trees in the ensemble
default :: 100
- sample_size :: sub-sample size used to build each tree... max samples
default :: 256
- max_depth :: per-tree height limit... if not defined is set to ceil(log2(psi))
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
# trains on $train, never the raw $data.
my $train = $self->_prepare_fit_data($data);
# Per-feature medians of the raw training data, stored for ablation
# explanations (explain_samples). From $data, not $train: the pure-
# Perl path's densified $train would bake fill values into the
# medians while the C path's raw pass-through would not, and the
# baselines must not depend on use_c.
$self->{feature_baselines} = $self->_compute_feature_baselines($data);
my $n = scalar @$train;
# Resolve sub-sample size, extension level and height limit against the
# data's shape. Shared with fit_from_csv(), which learns n from a census
# pass before it ever holds the rows in RAM.
my ( $psi, $limit ) = $self->_resolve_geometry( $n, $n_features );
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|/tagged_row_to_array(\%row, $caller)> (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|/tagged_row_to_array(\%row, $caller)>, 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 fit_from_csv($path, %opts)
Trains the model directly from a CSV file B<without loading it into RAM>,
for datasets too large to slurp. An Isolation Forest never trains on more
than C<n_trees * sample_size> rows -- each tree uses one independent
sub-sample of C<sample_size> points -- so the working set is bounded by the
model's parameters, not the file size.
The file is read in two (or three) passes:
=over 4
=item 1.
B<Census> -- count the rows (this is the true C<n> that fixes C<psi =
min(sample_size, n)>) and pin the feature width. By default this is an
I<index> pass (see C<index> below): a fast block-scan that also records each
data row's byte offset.
=item 2.
B<Gather> -- retain only the rows the trees actually sampled (chosen in
bounded memory via Floyd's algorithm), then build each tree from its
sub-sample exactly as L</fit> would. With the offset table from pass 1 this
B<seeks straight to the sampled rows> instead of re-scanning the file.
=item 3.
B<Threshold> (only when C<contamination> is set) -- score every row against
the built forest and place the cut at the contamination boundary. The cut
is B<exact>, not estimated: a single min-heap of the C<k+1> largest scores
(C<k = round(contamination * n)>) captures everything
L</decision_threshold> selection needs, with a rare extra pass to resolve a
tie block straddling the boundary. Under C<< use_c => 1 >> (the default when
the C backend is present) this pass runs through the C scorer, so it stays
fast even over millions of rows.
=back
Neither census variant parses cells into numbers or checks for missing
values: only the rows that actually train (via gather) or are scored (the
threshold pass) are validated, so a malformed or missing cell in a
never-used row is not reported. In particular, under C<< missing => 'die' >>
a missing cell is rejected exactly when it appears in a sampled training row.
Because rows are addressed by their position across passes, C<$path> must be
a stable, re-readable file (not a pipe or C<STDIN>).
The first line is skipped as a header when it holds any non-numeric text (a
feature-name row) or exactly matches the model's stored C<feature_names> --
data rows contain only numbers and empty cells. Pass C<< header => 1 >> to
force the first line to be skipped even when it is all-numeric.
Options:
=over 4
=item * C<< header => 1 >> -- always skip the first line as a column header
(otherwise a header is auto-detected as described above).
=item * C<< index => 0 >> -- disable the offset-index pass and use the
streaming two-pass reader instead. The index is on by default: it makes
gather a random-access read of just the sampled rows rather than a second
full scan, but costs an C<8 * n>-byte offset table. When that table would
exceed C<index_max> it is dropped automatically and the fit falls back to
streaming. (In index mode a blank line is one that is empty or contains only
whitespace, judged cheaply; pass C<< index => 0 >> for files where that
matters.)
=item * C<< index_max => bytes >> -- ceiling on the offset table (default
256 MiB). Above it, the index is abandoned mid-scan and gather streams.
=item * C<< c_scan => 0 >> -- validate each scored cell with
C<looks_like_number> in the threshold pass instead of letting the C packer
coerce it. C<c_scan> is on by default and takes effect only on the C
mean-voting scoring path, where it is a large saving over millions of rows;
it produces the identical result on valid data, but coerces a non-numeric
scored cell to C<0.0> (nudging the threshold) rather than dying on it.
=back
Every column is a feature; an empty cell is treated as a missing value and
handled by the model's C<missing> strategy. The resulting model is
identical in shape to one built by L</fit> and is fully deterministic given
the same file, C<seed>, and parameters (though not bit-identical to
L</fit>, whose RNG stream interleaves sampling and tree-building
differently).
Current limitations: C<mungers> and tagged/named columns are not supported
(load such data through L<fit_tagged|/fit_tagged(\@rows)>), and the trees are always built
serially in pure Perl -- C<parallel_fit>/C<use_openmp_fit> are ignored for
the build, though scoring is still accelerated when C<use_c> is on.
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
#-------------------------------------------------------------------------------
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.
#
# Args:
# $theta :: the per-tree score cutoff, normally in (0, 1]. Degenerate
# values are handled rather than rejected.
# $c :: c(psi) for this forest, i.e. $self->{c_psi}. Only <= 0 for a
# psi <= 1 forest.
#
# Returns: the path length at or below which a tree votes anomalous, a
# float. +inf when every tree should vote, -1.0 when none should.
#
# Example:
# _depth_cut( 0.6, 10.24 ); # ~7.55: isolate in <= 7.55 edges and vote
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").
#
# Args:
# $n_trees :: the forest's tree count, an integer >= 1.
#
# Returns: the required vote count, an integer in [1, $n_trees].
#
# Example:
# _min_votes(100); # 51
# _min_votes(101); # 51
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.
#
# Args:
# $desc :: arrayref of the per-point statistic sorted DESCENDING -- the
# mean-mode anomaly scores, or the majority pivots under
# voting => 'majority'. Must be non-empty.
# $k :: how many points the contamination rate wants flagged, an integer
# in [1, scalar @$desc].
#
# Returns: the cutoff, a float positioned so that `statistic >= cutoff`
# selects the intended points. Never equal to a value in $desc.
#
# Example:
# _threshold_from_ranked( [ 0.71, 0.68, 0.52, 0.50 ], 2 ); # 0.60
# _threshold_from_ranked( [ 0.71, 0.68, 0.68, 0.50 ], 2 ); # 0.59, tie
#-------------------------------------------------------------------------------
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.
#
# Args:
# $data :: arrayref of rows already through _prepare_perl_input -- dense
# under impute, raw undef preserved under nan.
# $cut :: the depth cut from _depth_cut.
#
# Returns: arrayref of per-point vote counts, integers in [0, n_trees],
# positionally matching $data.
#
# Example:
# my $cut = _depth_cut( 0.6, $self->{c_psi} );
# my $votes = $self->_vote_counts_perl( $rows, $cut ); # [ 3, 97, 12, ... ]
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
# 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.
#
# Args:
# $data :: the training set, an arrayref of feature-value arrayrefs. Raw
# undef cells are fine.
#
# Returns: nothing. Sets $self->{threshold} as its whole purpose.
#
# Example:
# $self->{contamination} = 0.05;
# $self->_learn_contamination_threshold( \@training_rows );
# $self->decision_threshold; # the cutoff flagging ~5% of the training set
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
#-------------------------------------------------------------------------------
# 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.
#
# Args:
# $data :: arrayref of feature-value arrayrefs, raw or prepared -- it goes
# through _prepare_perl_input here either way.
#
# Returns: arrayref of per-point pivots, each a float in (0, 1],
# positionally matching $data.
#
# Example:
# my $pivots = $self->_majority_pivot_scores( \@training_rows );
# # row i is flagged at per-tree cutoff theta exactly when
# # $pivots->[$i] >= theta
#-------------------------------------------------------------------------------
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.
#
# Args: none. Draws two uniforms from Perl's rand(), so the caller controls
# reproducibility through srand().
#
# Returns: one float from N(0,1), typically within +/-4. Consumes exactly
# two rand() draws, which is what keeps the Perl and C builders in step.
#
# Example:
# srand(42);
# my $coef = _randn(); # a hyperplane coefficient for one feature
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
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).
#
# Args:
# $x :: one sample, an arrayref of feature values. undef cells are
# allowed and handled per $nan.
# $node :: the node to start walking from, normally a tree root.
# $depth :: the depth credited to $node, 0 for a root. Callers only pass
# anything else to resume a partial walk.
# $nan :: true to use the nan-strategy routing described above; false or
# omitted to coerce undef to 0.
#
# Returns: the path length, a float -- edges walked plus c(leaf size), so
# it is not an integer whenever the walk ends in a multi-point leaf.
#
# Example:
# _path_length( [ 0.9, 0.4 ], $tree, 0, 0 ); # e.g. 3.51
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
#-------------------------------------------------------------------------------
# Explanation (explain_samples) internals.
#-------------------------------------------------------------------------------
# Instrumented twin of _path_length: routes $x through one tree with
# EXACTLY _path_length's split logic while recording which feature(s)
# each crossed node tested. Returns (path_length, \@pairs) where
# path_length carries the usual c(leaf size) adjustment and each pair
# is [feature_index, share] -- an axis node contributes one pair with
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
# every cell of the training set, so with the C backend on it goes
# through first_missing_xs -- same row-major order, same cell
# reported, without the per-cell Perl loop overhead. Both paths
# read a row that is not an arrayref as missing at column 0, which
# is how pack_input_xs treats one downstream.
my $where = [];
if ( $self->{_use_c} ) {
first_missing_xs( $data, scalar @$data, $nf, $where );
} else {
my $i = 0;
for my $row (@$data) {
if ( ref $row ne 'ARRAY' ) { $where = [ $i, 0 ]; last }
my $f = 0;
$f++ while $f < $nf && defined $row->[$f];
if ( $f < $nf ) { $where = [ $i, $f ]; last }
$i++;
}
} ## end else [ if ( $self->{_use_c} ) ]
croak "fit(): undef feature value at sample $where->[0], column $where->[1]; "
. "construct with missing => 'zero', 'impute', or 'nan' "
. "to train on data with missing values"
if @$where;
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. The
# optional $how_override selects the statistic independently of the
# model's impute_with knob (used by _compute_feature_baselines, which
# always wants the median).
#
# Args:
# $data :: the rows to learn from, an arrayref of feature-value
# arrayrefs. undef cells are skipped rather than counted.
# $how_override :: 'mean' or 'median' to force the statistic, or undef to
# follow the model's impute_with.
#
# Returns: an arrayref of n_features fill values. Croaks naming the column
# when a feature has no present value anywhere in $data.
#
# Example:
# $self->_compute_impute_fill( \@rows ); # per impute_with
# $self->_compute_impute_fill( \@rows, 'median' ); # forced median
sub _compute_impute_fill {
my ( $self, $data, $how_override ) = @_;
my $nf = $self->{n_features};
my $how = $how_override // $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
# Per-feature median of the present values of the training data -- the
# "typical row" that ablation explanations (explain_samples with
# method => 'ablation') substitute against, one feature at a time.
# Always the median (not impute_with's statistic): a baseline should be
# a robustly central value, and outliers in the training data drag a
# mean around far more than a median.
#
# Computed from the RAW rows, never a densified copy, so the stored
# baselines are identical whether use_c is on or off (the pure-Perl fit
# path densifies $train before the trees are built; the C path does
# not). A column with no present value at all -- legal under
# missing => 'zero'/'nan' -- cannot yield a median, so the fast path's
# croak falls back to a tolerant pure-Perl pass that gives such columns
# the fill value scoring maps their undefs to anyway (0).
#
# Args:
# $data :: the raw training rows, an arrayref of feature-value
( run in 1.706 second using v1.01-cache-2.11-cpan-e623d60df62 )