view release on metacpan or  search on metacpan

lib/Algorithm/Classifier/IsolationForest.pm  view on Meta::CPAN

use strict;
use warnings;
use Carp        qw(croak);
use Config      ();
use List::Util  qw(min);
use POSIX       qw(ceil);
use JSON::PP    ();
use File::Slurp qw(read_file write_file);

our $VERSION = '0.6.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.
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

lib/Algorithm/Classifier/IsolationForest.pm  view on Meta::CPAN

    const double *node_ptrs[n_trees];
    const int    *idx_ptrs[n_trees];
    const double *val_ptrs[n_trees];

    /* forest_bytes totals every buffer the tree walks touch; it decides
     * between the two loop shapes below. */
    size_t forest_bytes = 0;
    for (ti = 0; ti < n_trees; ti++) {
        SV** np = av_fetch(nodes_av, ti, 0);
        SV** ip = av_fetch(idx_av,   ti, 0);
        SV** vp = av_fetch(val_av,   ti, 0);
        if (!np || !*np || !ip || !*ip || !vp || !*vp) {
            croak("score_all_xs: missing tree %d", ti);
        }
        node_ptrs[ti] = (const double*)SvPVbyte(*np, tl); forest_bytes += tl;
        idx_ptrs[ti]  = (const int*)   SvPVbyte(*ip, tl); forest_bytes += tl;
        val_ptrs[ti]  = (const double*)SvPVbyte(*vp, tl); forest_bytes += tl;
    }

    xd = (const double*)SvPVbyte(x_sv, tl);
    sm = (double*)SvPVbyte_force(sm_sv, tl);

    /* Two loop shapes over the same per-point ascending-t additions --
     * bit-identical results either way, so the size heuristic choosing
     * between them can never change scores.
     *
     * Point-major (small forests): each point walks all trees with its
     * path-length sum held in a register.  Cheapest per walk, and the
     * whole forest stays cache-resident across points anyway.
     *
     * Tree-blocked (large forests): once the forest outgrows L3, the
     * point-major loop re-streams every tree's nodes and coefficients
     * from memory for every point -- an extended-mode tree is ~56 KB
     * at 16 features (24 KB nodes + 32 KB dense coefficients), and its
     * per-tree scoring cost measured 2.2x worse at 400 trees than at
     * 100.  Walking a block of points through ONE tree at a time keeps
     * that tree hot in L1/L2 while the block's rows stream through it
     * (measured 3.1x faster at 400 extended trees, 20k points).  The
     * blocked shape pays an sm[i] load+store per walk instead of a
     * register add, which measurably hurts cheap axis walks while the
     * forest still fits in cache -- hence the byte threshold rather
     * than always tiling. */
    if (forest_bytes <= (size_t)4 * 1024 * 1024) {
#ifdef _OPENMP
        #pragma omp parallel for schedule(static) if(use_openmp)
#endif
        for (int i = 0; i < n_pts; i++) {
            const double *xi = xd + (size_t)i * (size_t)n_feats;
            double sum = 0.0;
            for (int t = 0; t < n_trees; t++) {
                sum += if_walk_tree(node_ptrs[t], idx_ptrs[t],
                                    val_ptrs[t], xi, n_feats);
            }
            sm[i] = sum;
        }
    }
    else {
        /* 256 rows x 16 features x 8 bytes = 32 KB of input per block
         * -- comfortable in L2 next to one tree.  Each OpenMP thread
         * owns whole blocks and therefore a unique slice of sm[], so
         * there is still no synchronisation.  For small batches the
         * tile shrinks to keep ~4 blocks per thread available; losing
         * per-block tree reuse there is fine, since a small batch
         * never re-streams much anyway. */
        int tile = 256;
#ifdef _OPENMP
        if (use_openmp) {
            int min_blocks = omp_get_max_threads() * 4;
            if (min_blocks > 0 && (n_pts + tile - 1) / tile < min_blocks) {
                tile = (n_pts + min_blocks - 1) / min_blocks;
                if (tile < 1) tile = 1;
            }
        }
#endif
        int n_blocks = (n_pts + tile - 1) / tile;

#ifdef _OPENMP
        #pragma omp parallel for schedule(static) if(use_openmp)
#endif
        for (int blk = 0; blk < n_blocks; blk++) {
            const int i0 = blk * tile;
            const int i1 = (i0 + tile < n_pts) ? i0 + tile : n_pts;
            for (int i = i0; i < i1; i++) sm[i] = 0.0;
            for (int t = 0; t < n_trees; t++) {
                const double *nd  = node_ptrs[t];
                const int    *ico = idx_ptrs[t];
                const double *vco = val_ptrs[t];
                for (int i = i0; i < i1; i++) {
                    sm[i] += if_walk_tree(nd, ico, vco,
                                          xd + (size_t)i * (size_t)n_feats,
                                          n_feats);
                }
            }
        }
    }
}

/* vote_all_xs(nodes_av, idx_av, val_av, x_sv, sm_sv,
 *             n_pts, n_feats, n_trees, depth_cut, min_votes, use_openmp)
 *
 * Majority-voting (MVIForest) tree walk: a tree votes a point anomalous
 * when the point's path length in that tree is <= depth_cut -- the
 * depth-domain image of the per-tree score cutoff (the Perl side
 * precomputes depth_cut = -c(psi) * log2(threshold), so no per-tree
 * exp()/log() runs in here).
 *
 * min_votes == 0: sm[i] = the point's full vote count over all n_trees
 *   trees (a small integer stored as a double, so the existing
 *   finalize_* helpers work on the buffer unchanged).
 * min_votes > 0:  sm[i] = 1.0/0.0 anomaly label, with per-point early
 *   exit: the walk stops as soon as the point has min_votes votes (the
 *   remaining trees can't change the outcome) or can no longer reach
 *   min_votes.  This is MVIForest's "stop at majority" scoring loop.
 *
 * Always point-major, unlike score_all_xs's two loop shapes: the vote
 * count / early exit is per-point state, so a tree-blocked loop would
 * have to re-load it per walk and could never exit a point early.
 * Votes are integer counts, so there is no summation-order concern
 * either way.  Thread-safety matches score_all_xs: the parallel region
 * reads extracted pointers and writes a unique sm[i] per iteration. */
void vote_all_xs(SV* nodes_av_sv, SV* idx_av_sv, SV* val_av_sv,