Algorithm-Classifier-IsolationForest
view release on metacpan or search on metacpan
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
AV* out;
int* all;
int t, i;
if (!SvROK(out_rv) || SvTYPE(SvRV(out_rv)) != SVt_PVAV) {
croak("build_forest_xs: out must be an arrayref");
}
x = (const double*)SvPVbyte(x_sv, tl);
out = (AV*)SvRV(out_rv);
av_clear(out);
if (n_trees > 0) av_extend(out, n_trees - 1);
all = (int*)malloc(n_pts * sizeof(int));
for (t = 0; t < n_trees; t++) {
int* sample;
for (i = 0; i < n_pts; i++) all[i] = i;
for (i = 0; i < psi; i++) {
int j = i + (int)(Drand01() * (n_pts - i));
int tmp = all[i]; all[i] = all[j]; all[j] = tmp;
}
sample = (int*)malloc(psi * sizeof(int));
memcpy(sample, all, psi * sizeof(int));
av_store(out, t,
_build_node_c(aTHX_ x, n_feats, sample, psi, 0, limit,
mode_flag, ext_level));
}
free(all);
}
/* ---------------------------------------------------------------------
* build_forest_openmp_xs -- OpenMP-parallel fit() tree builder.
*
* build_forest_xs (above) is bit-identical to the pure-Perl path
* because every random draw goes through Drand01(), the same
* generator Perl's rand()/srand() use -- but that generator is a
* single mutable struct shared by the whole interpreter, so calling
* it concurrently from multiple OpenMP threads would be a data race.
* The same is true of any Perl API call (newAV, newSViv, ...): Perl's
* SV allocator isn't safe to call from multiple OS threads sharing one
* interpreter without a lock that would just serialise everything
* anyway.
*
* So this builder trades the bit-identical guarantee for real thread
* parallelism: each tree gets its own splitmix64 PRNG stream, seeded
* from a tree index (not thread id or scheduling order), so results
* are still reproducible for a fixed seed and n_trees regardless of
* OMP_NUM_THREADS -- just different from what build_forest_xs or the
* pure-Perl path would produce for the same seed. The one Drand01()
* call in this function happens before the parallel region starts
* (single-threaded), so the result still varies with the model's
* `seed` the way every other code path does; it isn't used inside the
* parallel loop.
*
* Each tree is built entirely with plain C data (row-index int arrays,
* a growable TreeBuf of packed doubles/ints) -- no Perl API call
* happens anywhere inside the parallel region. Each node record in
* TreeBuf uses _pack_tree's 6-double SoA layout (see the file-top
* comment), but the node ORDER differs: records are appended
* post-order (a node is pushed after both its children, since child
* indices must be known first), so the root is the last record --
* _pack_tree's pre-order puts it at 0. _unpack_forest accounts for
* this. Oblique coefficients are also always stored sparse (in the
* random pool's order) -- the dense-pack fast path is skipped because
* its only purpose is speeding up score_all_xs, and _rebuild_c_trees
* reapplies it anyway once the caller unpacks these buffers back into
* the standard Perl tree shape and re-derives the scoring buffers.
*
* After the parallel region, each tree's TreeBuf is copied into a Perl
* string SV (one memcpy each, serially) and stored into nodes_rv /
* idx_rv / val_rv -- the caller unpacks these into ordinary nested
* Perl trees for $self->{trees} (so to_json/persistence/_rebuild_c_trees
* are unaffected). ------------------------------------------------ */
typedef struct {
double *nodes; size_t n_nodes, cap_nodes;
int *idx; size_t n_idx, cap_idx;
double *val; size_t n_val, cap_val;
} TreeBuf;
static void tb_init(TreeBuf *b) {
b->nodes = NULL; b->n_nodes = 0; b->cap_nodes = 0;
b->idx = NULL; b->n_idx = 0; b->cap_idx = 0;
b->val = NULL; b->n_val = 0; b->cap_val = 0;
}
static void tb_free(TreeBuf *b) {
free(b->nodes); free(b->idx); free(b->val);
}
static int tb_push_node(TreeBuf *b, double f0, double f1, double f2,
double f3, double f4, double f5) {
double *slot;
if (b->n_nodes == b->cap_nodes) {
size_t newcap = b->cap_nodes ? b->cap_nodes * 2 : 64;
b->nodes = (double*)realloc(b->nodes, newcap * 6 * sizeof(double));
b->cap_nodes = newcap;
}
slot = b->nodes + b->n_nodes * 6;
slot[0] = f0; slot[1] = f1; slot[2] = f2;
slot[3] = f3; slot[4] = f4; slot[5] = f5;
return (int)(b->n_nodes++);
}
/* Appends n (idx[k], val[k]) pairs and returns the offset they start
* at -- the `coff` an oblique node record stores. */
static int tb_push_coef(TreeBuf *b, const int *idx, const double *val,
int n) {
int off = (int)b->n_idx;
if (b->n_idx + (size_t)n > b->cap_idx) {
size_t newcap = b->cap_idx ? b->cap_idx * 2 : 64;
if (newcap < b->n_idx + (size_t)n) newcap = b->n_idx + (size_t)n;
b->idx = (int*)realloc(b->idx, newcap * sizeof(int));
b->cap_idx = newcap;
}
if (b->n_val + (size_t)n > b->cap_val) {
size_t newcap = b->cap_val ? b->cap_val * 2 : 64;
if (newcap < b->n_val + (size_t)n) newcap = b->n_val + (size_t)n;
b->val = (double*)realloc(b->val, newcap * sizeof(double));
b->cap_val = newcap;
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
* 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) {
lib/Algorithm/Classifier/IsolationForest.pm view on Meta::CPAN
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
#-------------------------------------------------------------------------------
# Packed input wrapper. pack_data() returns one of these so callers can
# score the same dataset many times without re-walking the AV/AV refs on
# every call -- a meaningful win at high feature counts where
# pack_input_xs is a non-trivial slice of total scoring time.
#
# It's a minimal blessed hashref: { packed, n_pts, n_feats }. The C
# scoring functions only need the packed bytes + dimensions.
#-------------------------------------------------------------------------------
sub pack_data {
my ( $self, $data ) = @_;
$self->_check_fitted;
croak "pack_data requires the Inline::C backend; install Inline::C"
unless $self->{_use_c};
croak "pack_data() expects an arrayref of samples"
unless ref $data eq 'ARRAY';
my $n_pts = scalar @$data;
my $nf = $self->{n_features};
my $x_packed = "\0" x ( $n_pts * $nf * 8 );
my ( $mode, $fill ) = $self->_pack_args;
pack_input_xs( $data, $x_packed, $n_pts, $nf, $mode, $fill );
return bless {
packed => $x_packed,
n_pts => $n_pts,
n_feats => $nf,
},
'Algorithm::Classifier::IsolationForest::PackedData';
} ## end sub pack_data
# Internal helper: given $data that may be a raw arrayref OR a PackedData
# instance, return the (n_pts, n_feats, x_packed) triple ready for
# score_all_xs. Called from every scoring fast path.
sub _resolve_input {
my ( $self, $data ) = @_;
if ( ref $data eq 'Algorithm::Classifier::IsolationForest::PackedData' ) {
croak "PackedData has $data->{n_feats} features but model expects " . $self->{n_features}
unless $data->{n_feats} == $self->{n_features};
return ( $data->{n_pts}, $data->{n_feats}, $data->{packed} );
}
my $n_pts = scalar @$data;
my $nf = $self->{n_features};
my $x_packed = "\0" x ( $n_pts * $nf * 8 );
my ( $mode, $fill ) = $self->_pack_args;
pack_input_xs( $data, $x_packed, $n_pts, $nf, $mode, $fill );
return ( $n_pts, $nf, $x_packed );
( run in 0.816 second using v1.01-cache-2.11-cpan-acf6aa7dc9e )