Algorithm-Classifier-IsolationForest

 view release on metacpan or  search on metacpan

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

 * 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;

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

            SV** lv = av_fetch(lo, f, 1);
            SV** hv = av_fetch(hi, f, 1);
            if (x[f] < SvNV(*lv)) sv_setnv(*lv, x[f]);
            if (x[f] > SvNV(*hv)) sv_setnv(*hv, x[f]);
        }
    }

    if (SvIV(slots[OL_TYPE]) == 0) {    /* leaf */
        if ((double)count >= _ol_split_threshold(eta, adaptive, depth) &&
            (double)depth < limit) {
            AV* lo      = (AV*)SvRV(AvARRAY(node)[OL_LO]);
            AV* hi      = (AV*)SvRV(AvARRAY(node)[OL_HI]);
            double* lod = (double*)malloc(nf * sizeof(double));
            double* hid = (double*)malloc(nf * sizeof(double));
            double* pts = (double*)malloc((size_t)count * (size_t)nf * sizeof(double));
            int* idx    = (int*)malloc((size_t)count * sizeof(int));
            SV* subtree;
            int i;

            for (f = 0; f < nf; f++) {
                SV** lv = av_fetch(lo, f, 0);
                SV** hv = av_fetch(hi, f, 0);
                lod[f] = (lv && *lv && SvOK(*lv)) ? SvNV(*lv) : 0.0;
                hid[f] = (hv && *hv && SvOK(*hv)) ? SvNV(*hv) : 0.0;
            }
            /* Synthetic points, point-major, one draw per feature with
             * width > 0 -- _sample_box's exact draw order. */
            for (i = 0; i < (int)count; i++) {
                for (f = 0; f < nf; f++) {
                    double w = hid[f] - lod[f];
                    pts[(size_t)i * (size_t)nf + f]
                        = (w > 0) ? lod[f] + Drand01() * w : lod[f];
                }
                idx[i] = i;
            }
            subtree = _ol_build(aTHX_ pts, idx, (int)count, nf, depth, limit,
                                eta, adaptive);
            free(lod);
            free(hid);
            free(pts);
            free(idx);
            return subtree;
        }
        return node_rv;
    }

    {
        int attr     = (int)SvIV(slots[OL_ATTR]);
        double split = SvNV(slots[OL_SPLIT]);
        int ci       = (x[attr] < split) ? OL_LEFT : OL_RIGHT;
        SV* child    = AvARRAY(node)[ci];
        SV* nc = _ol_node_learn(aTHX_ child, x, nf, depth + 1, limit, eta,
                                adaptive);
        if (nc != child) av_store(node, ci, nc);
        return node_rv;
    }
}

/* Union of two nodes' boxes into caller-provided arrays, matching
 * _box_union(): nodes without a box are skipped, the first boxed node
 * is copied and the second folded in.  Returns 0 when neither node has
 * a box. */
static int _ol_box_union(pTHX_ SV* a_rv, SV* b_rv, int nf,
                         double* lo, double* hi) {
    SV* boxed[2];
    int nb = 0, bi, f;

    if (SvOK(AvARRAY((AV*)SvRV(a_rv))[OL_LO])) boxed[nb++] = a_rv;
    if (SvOK(AvARRAY((AV*)SvRV(b_rv))[OL_LO])) boxed[nb++] = b_rv;
    if (nb == 0) return 0;

    for (bi = 0; bi < nb; bi++) {
        AV* node = (AV*)SvRV(boxed[bi]);
        AV* blo  = (AV*)SvRV(AvARRAY(node)[OL_LO]);
        AV* bhi  = (AV*)SvRV(AvARRAY(node)[OL_HI]);
        for (f = 0; f < nf; f++) {
            SV** lv  = av_fetch(blo, f, 0);
            SV** hv  = av_fetch(bhi, f, 0);
            double l = (lv && *lv && SvOK(*lv)) ? SvNV(*lv) : 0.0;
            double h = (hv && *hv && SvOK(*hv)) ? SvNV(*hv) : 0.0;
            if (bi == 0) {
                lo[f] = l;
                hi[f] = h;
            } else {
                if (l < lo[f]) lo[f] = l;
                if (h > hi[f]) hi[f] = h;
            }
        }
    }
    return 1;
}

/* Aggregate a subtree back into one leaf -- _collapse().  Children are
 * collapsed first so their boxes can be unioned; intermediate leaves
 * built while collapsing internal children are temporaries and dropped
 * here.  Returns a NEW leaf ref (refcount 1) for internal nodes, or the
 * node itself when it is already a leaf. */
static SV* _ol_collapse(pTHX_ SV* node_rv, int nf) {
    AV* node = (AV*)SvRV(node_rv);
    SV *l_rv, *r_rv, *cl, *cr, *leaf;
    SV *lo_rv = NULL, *hi_rv = NULL;
    double *lo, *hi;

    if (SvIV(AvARRAY(node)[OL_TYPE]) == 0) return node_rv;

    l_rv = AvARRAY(node)[OL_LEFT];
    r_rv = AvARRAY(node)[OL_RIGHT];
    cl   = _ol_collapse(aTHX_ l_rv, nf);
    cr   = _ol_collapse(aTHX_ r_rv, nf);

    lo = (double*)malloc(nf * sizeof(double));
    hi = (double*)malloc(nf * sizeof(double));
    if (_ol_box_union(aTHX_ cl, cr, nf, lo, hi)) {
        lo_rv = _ol_mk_box(aTHX_ lo, nf);
        hi_rv = _ol_mk_box(aTHX_ hi, nf);
    } else if (SvOK(AvARRAY(node)[OL_LO])) {
        /* Both children boxless: keep the node's own box (the Perl code
         * moves the same arrayrefs into the new leaf). */
        lo_rv = SvREFCNT_inc(AvARRAY(node)[OL_LO]);
        hi_rv = SvREFCNT_inc(AvARRAY(node)[OL_HI]);
    }

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

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|/load($path)>,
which dispatches on the stored format tag.

Throughout these docs, B<psi> is the paper's ψ: the number of points each
tree is built from, which C<sample_size> sets and which other
implementations often call I<max samples>.  See L</REFERENCES>.

=head1 NATIVE ACCELERATION (Inline::C and OpenMP)

Both the scoring hot path (C<score_samples>, C<predict>, C<path_lengths>,
C<score_predict_samples>, and C<score_predict_split>) and the C<fit()>
tree builder are automatically accelerated through
L<Inline::C> when it is installed and a working C compiler is reachable.
If the toolchain also accepts C<-fopenmp> and can link against
C<libgomp>, the per-point tree walk runs in parallel across all
available CPU cores using OpenMP, and the extended-mode oblique dot
product is vectorised via C<#pragma omp simd> -- which on modern x86
compilers translates to an unrolled FMA / AVX gather chain that's
substantially faster for high-feature-count extended models.

C<fit()>'s tree builder (subsampling plus the recursive axis/oblique
split search) runs in C the same way when C<use_c> is on, replacing the
per-node Perl arrayref copying with plain int-array partitioning --
typically an order of magnitude faster, and dramatically more so at
higher feature counts where the pure-Perl per-cell loop dominates. Its
random draws go through the same generator C<rand()>/C<srand()> use
internally, in the same call order the pure-Perl builder uses, so a
given C<seed> produces bit-identical trees whether C<use_c> is on or
off -- switching backends changes only how fast the model is built, not
the model itself. On perls whose NV is wider than a C double
(C<-Duselongdouble> / C<-Dusequadmath>) the pure-Perl builder rounds
each stored value to double precision to preserve this parity; axis
mode matches exactly, while extended mode can still differ on rare
libm rounding ties (double vs long-double transcendentals).

By default this C builder is single-threaded per call, because Perl's
RNG state isn't safe to share across OpenMP threads. Two ways to scale
fit() across cores are available (see below for why they don't compose):

=over 4

=item * C<parallel_fit> forks N worker processes, each building its
share of the trees with the (still single-threaded) C builder. Fixed
IPC/serialisation overhead per worker means this can cost more than it
saves once a fit already completes in milliseconds; it's most useful
once a single-process fit is large enough that the fork/Storable
overhead is small relative to the work being split.

=item * C<use_openmp_fit> builds trees across OpenMP threads within a
single process (one tree per thread), using a separate, thread-safe
PRNG seeded per tree index instead of Perl's C<rand()>. This means
trees built with C<use_openmp_fit> are I<not> bit-identical to the
default C<use_c> path for the same seed -- but a fixed seed and
C<n_trees> still reproduce the same trees regardless of
C<OMP_NUM_THREADS> or how OpenMP schedules the work. It's off by
default (unlike C<use_c>/C<use_openmp>, which only ever change speed,
this changes which trees get built) and only takes effect when C<use_c>
is also on and OpenMP is linked in.

=back

These two do NOT compose, despite both existing to parallelise fit().
A process that has run any OpenMP region -- including plain
C<score_samples()>/C<predict()> with the default C<use_openmp> -- and
then C<fork()>s (as C<parallel_fit> does) hands each child a copy of
libgomp's thread pool whose worker threads did not survive the fork. A
child that then starts its own C<#pragma omp parallel> region (as
C<use_openmp_fit> would) tries to reuse that now-invalid pool and
hangs. This is a general limitation of combining C<fork()> with OpenMP,
not something fixable from Perl, so C<parallel_fit>'s forked workers
always use the single-threaded C builder regardless of
C<use_openmp_fit> -- setting both just means C<parallel_fit> wins and
C<use_openmp_fit> has no effect for that call.

Detection happens once when the module is loaded.  When the
distribution was installed with C<Inline::C> available, the C backend
was already compiled during C<make> and the installed object is loaded
directly (see L</Compile at install time (the prebuilt object)> below);
otherwise the backend is compiled on first load and the artefact is
cached under C<_Inline/> and reused on subsequent runs.  Five package
variables report what the load picked up:

    $Algorithm::Classifier::IsolationForest::HAS_C       # 0/1
    $Algorithm::Classifier::IsolationForest::HAS_OPENMP  # 0/1
    $Algorithm::Classifier::IsolationForest::HAS_SIMD    # 0/1 (OpenMP 4.0+)
    $Algorithm::Classifier::IsolationForest::OPT_LEVEL   # e.g. "-O3 -march=native", '' if HAS_C is 0
    $Algorithm::Classifier::IsolationForest::C_SOURCE    # 'prebuilt' / 'runtime', '' if HAS_C is 0

Neither dependency is required.  Without C<Inline::C> the module falls
back to a pure-Perl implementation that produces identical results, just
slower; without OpenMP the C backend runs single-threaded.

The bundled C<iforest accel> subcommand performs a tiny fit + score and
prints which backend is active (including the build flags below), which
is the recommended way to verify the build picked up the optional
dependencies on a given machine.

=head2 Compile at install time (the prebuilt object)

When C<Inline::C> is usable while the distribution itself is being
built, C<perl Makefile.PL> arranges for the C backend to be compiled
once during C<make> and installed alongside the module like any XS
object.  At run time that object is loaded directly through
L<XSLoader>: no C compiler, no C<Inline> modules, and no C<_Inline/>
cache directory are needed on the machine the module ends up running

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

			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

sub schema_version { return $_[0]->{schema_version} }

=head2 schema_description

Returns the free-text description of the variable schema stored with the
model, or undef if none was recorded.

=cut

sub schema_description { return $_[0]->{schema_description} }

=head2 feature_descriptions

Returns the hashref of per-feature description strings stored with the
model, or undef if none were recorded.  Keys are feature names; coverage
may be partial.

=cut

sub feature_descriptions { return $_[0]->{feature_descriptions} }

=head2 fit

Trains the model on the specified data.

The data taken is an array of arrays. Each sub-array is one sample and must
contain one or more numeric features. All samples must have the same number
of features. There is no upper limit on dimensionality.

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

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.

    my $iforest = Algorithm::Classifier::IsolationForest->new(
        n_trees       => 100,
        sample_size   => 256,
        contamination => 0.01,
        seed          => 42,
    );
    $iforest->fit_from_csv('huge.csv', header => 1);

=cut

sub fit_from_csv {
	my ( $self, $path, %opt ) = @_;

	croak "fit_from_csv() requires a path to a CSV file"
		unless defined $path && length $path;
	croak "fit_from_csv(): '$path' is not a readable file"
		unless -f $path && -r _;
	croak "fit_from_csv() does not support munger models; " . "load the data through fit_tagged/fit instead"
		if _plan($self);

	# Decide once whether the first line is a header, so every pass skips the
	# same line and row indices stay aligned.  header => 1 forces it; otherwise
	# a first line carrying any non-numeric text (or matching stored
	# feature_names) is a header -- data rows hold only numbers and blanks.
	my $skip_first = $self->_detect_header( $path, $opt{header} );

	# ---- Pass 1: census.  Establish the true row count (which fixes psi) and
	# pin the feature width.  Neither census variant parses cells into numbers
	# or checks for missing values -- only the rows that actually train
	# (_numify_row, in gather) or get scored (the threshold pass) are validated,

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

	croak "fit_from_csv(): no data rows in '$path'" unless $n;

	# A stored feature_names schema fixes the expected column count.
	if ( ref $self->{feature_names} eq 'ARRAY' && @{ $self->{feature_names} } ) {
		my $want = scalar @{ $self->{feature_names} };
		croak "fit_from_csv(): $want feature_names but CSV has $n_features columns"
			unless $want == $n_features;
	}

	$self->{n_features} = $n_features;
	my ( $psi, $limit ) = $self->_resolve_geometry( $n, $n_features );

	# ---- Choose which row indices each tree trains on: an independent,
	# uniform psi-subset drawn without replacement.  Floyd's algorithm draws
	# each subset in O(psi) memory, never materialising the 0..n-1 index
	# vector _subsample() uses -- for an out-of-core n that vector would not
	# fit either.  %want maps each needed row index to the (tree, slot) pairs
	# awaiting it, so a row several trees picked is gathered once and shared.
	srand( $self->{seed} ) if defined $self->{seed};
	my @tree_idx;    # tree => [ chosen row indices ]
	my %want;        # row index => [ [tree, slot], ... ]
	for my $t ( 0 .. $self->{n_trees} - 1 ) {
		my @idx = _sample_indices_distinct( $n, $psi );
		$tree_idx[$t] = \@idx;
		push @{ $want{ $idx[$_] } }, [ $t, $_ ] for 0 .. $#idx;
	}

	# ---- Pass 2: gather the sampled rows (validated + coerced by _numify_row).
	# With an offset table this seeks straight to them; otherwise it re-scans.
	my $pool
		= $offsets
		? $self->_gather_indexed( $path, $offsets, \%want )
		: $self->_gather_stream( $path, $skip_first, \%want );

	# Baselines for ablation explanations, learned from the gathered
	# sub-sample -- like the impute fill below, the only rows bounded-
	# memory fitting ever holds.  Before _apply_missing_to_pool, which
	# densifies the pool in place (the medians must come from present
	# values only, matching fit()).  Median is an exact order statistic,
	# so the hash-order values() walk cannot change the result.
	$self->{feature_baselines} = $self->_compute_feature_baselines( [ values %$pool ] );

	# Apply the missing-value strategy to the retained rows (see fit()'s
	# _prepare_fit_data; the impute fill is learned from this training
	# sub-sample, the only value bounded-memory fitting can compute).
	$self->_apply_missing_to_pool($pool);

	# ---- Build. Each tree trains on its gathered sub-sample -- structurally
	# identical to fit()'s per-tree _build_tree(_subsample(...)).
	my @trees;
	for my $t ( 0 .. $self->{n_trees} - 1 ) {
		my $sample = [ @{$pool}{ @{ $tree_idx[$t] } } ];
		push @trees, $self->_build_tree( $sample, 0, $limit );
	}
	$self->{trees} = \@trees;

	# Repack the just-built trees for the C scorer up front.  fit() scores its
	# learned threshold pure-Perl because its training set is small; the
	# streaming contamination pass here may score millions of rows, so paying
	# one repack now lets that pass run through the C path instead (bit-
	# identical result, minutes -> seconds).  The delete first clears any stale
	# buffers from a prior fit so the repack reflects the new forest.
	delete @$self{qw(_c_nodes _c_coef_idx _c_coef_val)};
	$self->_rebuild_c_trees() if $self->{_use_c};

	# c_scan (default on): let the C scorer's packer coerce raw cells via SvNV
	# in the contamination pass instead of validating each with looks_like_number
	# in Perl -- a large saving over millions of rows.  It only takes effect on
	# the C mean-voting path; see _stream_scores.
	my $c_scan = exists $opt{c_scan} ? $opt{c_scan} : 1;
	$self->_learn_contamination_threshold_streaming( $path, $skip_first, $n, $c_scan )
		if defined $self->{contamination};

	return $self;
} ## end sub fit_from_csv

=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};

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

            $e->{score}, $top->{name} // $top->{index}, $top->{weight};
    }

Each hashref is

    {
        score    => 0.91,          # same value score_samples returns
        method   => 'ablation',    # which attribution method produced this
        features => [              # every feature, most responsible first
            { index => 1, name => 'bytes_log', weight => 0.62, value => 9.3,
              delta => 0.31, baseline => 4.1 },
            { index => 0, name => 'method',    weight => 0.21, value => 2.0,
              delta => 0.10, baseline => 1.0 },
            ...
        ],
    }

C<name> is taken from the model's stored C<feature_names> and is undef
when none were recorded.  C<weight> is each feature's normalised share
of the responsibility, in [0, 1] and summing to 1 (all weights are 0
in the degenerate case where no attribution information exists).
C<value> is the feature value the model actually scored.

Options:

  - method :: how per-feature responsibility is computed
        ablation :: (default) counterfactual substitution.  Each
            feature in turn is replaced by its baseline (the
            per-feature training-data median learned at fit time) and
            the sample is re-scored; the drop in anomaly score is that
            feature's C<delta>, and weights are the positive deltas
            normalised.  Features whose substitution does not
            de-anomalise the sample get 0.  Each feature entry
            additionally carries C<delta> (may be negative:
            substituting a correlated feature can make the sample MORE
            anomalous) and C<baseline> (the value swapped in).  Costs
            n_samples * (n_features + 1) rows through the ordinary
            scorer -- C-accelerated when available.
        path :: split attribution (local DIFFI; Carletti, Terzi &
            Susto -- see REFERENCES).  The sample walks every tree;
            each split node crossed credits its feature(s) with
            1/h_t - 1/h_max, where h_t is that tree's path length for
            the sample and h_max the longest across the forest -- so
            the trees that isolated the sample quickly speak loudest
            and trees that treated it as normal say nothing.  An
            oblique (extended-mode) split spreads its credit across
            its features in proportion to |coefficient * value|.
            Needs nothing stored beyond the trees (so it works on
            models saved before baseline support) and runs pure Perl
            over the tree walks.

Why ablation is the default: it directly answers "would this sample
still be an outlier with feature j at a normal value?", and that
question has an answer for ANY scored sample.  C<path> can only
apportion the splits the forest actually built, so it attributes well
for samples that were IN the training data (post-fit forensics of
flagged training rows) but degrades for unseen samples in the tails --
an out-of-range sample walks through boundary regions whose splits
mostly test OTHER features, and the attribution dilutes or even
misleads.  Prefer C<ablation> whenever the model has baselines;
C<path> is the fallback for old saved models and a second opinion on
training-set outliers.  Neither method untangles strongly correlated
features -- a sample anomalous only in the JOINT distribution of two
features may show weak attributions on both.

Under C<< voting => 'majority' >> the score (and ablation's deltas) are
in vote-fraction space, matching C<score_samples>' semantics in that
mode; the C<path> method is aggregation-independent and unchanged.

C<ablation> requires stored baselines: models fitted before explanation
support have none and croak (models with C<< missing => 'impute' >>
fall back to their stored fill vector); refitting stores them.

=cut

sub explain_samples {
	my ( $self, $data, %opts ) = @_;
	$self->_check_fitted;

	my $method = delete $opts{method} // 'ablation';
	croak "explain_samples: method must be 'path' or 'ablation'"
		unless $method =~ /\A(?:path|ablation)\z/;
	croak "explain_samples: unknown option(s): " . join( ', ', sort keys %opts )
		if %opts;

	my $rows = $self->_to_arrayref($data);
	croak "explain_samples() expects a non-empty arrayref of samples"
		unless @$rows;

	return $method eq 'ablation'
		? $self->_explain_ablation($rows)
		: $self->_explain_path($rows);
} ## end sub explain_samples

=head2 explain_sample_tagged(\%row, %opts)

Explains a single sample supplied as a hashref of named feature values
-- the tagged counterpart of
L<explain_samples|/explain_samples(\@data, %opts)>, taking the same
C<method> option and returning the single explanation hashref (with
C<name> filled from the stored feature names).

    my $e = $forest->explain_sample_tagged({ cpu => 0.9, mem => 0.4 });
    print "worst feature: $e->{features}[0]{name}\n";

Croaks under the same conditions as L<tagged_row_to_array|/tagged_row_to_array(\%row, $caller)>.

=cut

sub explain_sample_tagged {
	my ( $self, $row, %opts ) = @_;
	my $vec = $self->tagged_row_to_array( $row, 'explain_sample_tagged' );
	return $self->explain_samples( [$vec], %opts )->[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} ) {
			my ( $n_pts, $nf, $x_packed ) = $self->_resolve_input($data);
			my $votes_packed = "\0" x ( $n_pts * 8 );
			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} );
			my $result = [];
			vote_score_predict_xs( $votes_packed, $n_pts, $t + 0.0, $maj + 0.0, $result );
			return $result;
		}

		my $votes = $self->_vote_counts_perl( $self->_prepare_perl_input($data), $cut );
		return [ map { [ $_ / $t, ( $_ >= $maj ? 1 : 0 ) ] } @$votes ] if _NV_IS_DOUBLE;

		# Wide-NV perls: match vote_score_predict_xs's double division --
		# see _NV_IS_DOUBLE.
		return [ map { [ _to_double( $_ / $t ), ( $_ >= $maj ? 1 : 0 ) ] } @$votes ];
	} ## end if ( $self->{voting} eq 'majority' )

	# Fast path: build [score, label] pairs straight from the sum buffer
	# in one C call.  Avoids the intermediate scores arrayref + Perl
	# foreach that allocates ~3*n_pts SVs.  Gated identically to predict()
	# so the threshold conversion is valid.

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

			my $ri = $assign->( $node->[5] );
			$node_data[$my_idx] = [ 2.0, $coef_off + 0.0, $num + 0.0, $li + 0.0, $ri + 0.0, $b + 0.0, ];
		} ## end else [ if ( $node->[0] == _NODE_LEAF ) ]
		return $my_idx;
	}; ## end $assign = sub
	$assign->($root);

	my $nodes_packed = pack( 'd*', map { @$_ } @node_data );
	my $idx_packed   = @coef_idx ? pack( 'l*', @coef_idx ) : pack('l*');
	my $val_packed   = @coef_val ? pack( 'd*', @coef_val ) : pack('d*');
	return ( $nodes_packed, $idx_packed, $val_packed );
} ## end sub _pack_tree

# Build packed C-ready representations for all trees and store them in
# $self->{_c_nodes}, $self->{_c_coef_idx}, $self->{_c_coef_val}.
# Called after fit() and from_json() when _use_c is true.  n_features is
# threaded through so _pack_tree can spot the dense-pack opportunity.
#
# Args: none beyond the model itself, which must already hold its trees.
#
# Returns: nothing.  Leaves _c_nodes, _c_coef_idx and _c_coef_val holding
# one packed string per tree, in tree order -- the form score_all_xs and
# friends expect.
#
# Example:
#   $self->{trees} = $self->_build_forest_c( $train, $psi, $limit, 100 );
#   $self->_rebuild_c_trees;   # scoring can now take the C path
sub _rebuild_c_trees {
	my ($self) = @_;
	my ( @c_nodes, @c_coef_idx, @c_coef_val );
	for my $tree ( @{ $self->{trees} } ) {
		my ( $np, $ip, $vp ) = _pack_tree( $tree, $self->{n_features} );
		push @c_nodes,    $np;
		push @c_coef_idx, $ip;
		push @c_coef_val, $vp;
	}
	$self->{_c_nodes}    = \@c_nodes;
	$self->{_c_coef_idx} = \@c_coef_idx;
	$self->{_c_coef_val} = \@c_coef_val;
} ## end sub _rebuild_c_trees

# Guard for every method that needs a forest to work with.  Called at the
# top of the scoring and packing entry points so an unfitted model fails
# with one clear message instead of an obscure error deeper in.
#
# Args: none beyond the model itself.
#
# Returns: nothing on success.  Croaks when the model holds no trees --
# either never fitted, or constructed but not yet loaded from JSON.
#
# Example:
#   $self->_check_fitted;   # croaks with "model is not fitted yet"
sub _check_fitted {
	my ($self) = @_;
	croak "model is not fitted yet; call fit() first"
		unless ref $self->{trees} eq 'ARRAY' && @{ $self->{trees} };
}

# ---------------------------------------------------------------------------
# Optional Algorithm::ToNumberMunger integration -- see the MUNGERS POD
# section.  Both helpers are plain functions (not methods) so the Online
# class's delegating methods can hand them their own $self: the plan and
# spec live in the same hash slots on either class.
# ---------------------------------------------------------------------------

# Validate a feature_descriptions hash against the feature names: every
# described feature must exist (a description for one that does not is
# either a typo or a stale leftover from a schema change) and every
# description must be a plain string.  Partial coverage is fine.  A plain
# function, like the munger helpers, so both classes and the prototype
# validator can call it.
#
# Args:
#   $tags :: the model's feature_names, an arrayref of strings.  Must be
#            non-empty -- there is nothing to validate against otherwise.
#   $fd :: the feature_descriptions hashref, keyed by feature name with
#          plain strings as values.
#
# Returns: 1 when the descriptions are usable.  Croaks naming the offending
# key when one is not a known feature, when a value is a reference, or when
# either argument is the wrong shape.
#
# Example:
#   _validate_feature_descriptions( [ 'cpu', 'mem' ],
#       { cpu => 'load average over 1m' } );   # 1; partial coverage is fine
sub _validate_feature_descriptions {
	my ( $tags, $fd ) = @_;
	croak "feature_descriptions must be a hashref of 'feature name => description'"
		unless ref $fd eq 'HASH';
	croak "feature_descriptions requires feature_names to validate against"
		unless ref $tags eq 'ARRAY' && @$tags;
	my %known = map { $_ => 1 } @$tags;
	for my $k ( sort keys %$fd ) {
		croak "feature_descriptions describes '$k', which is not in feature_names"
			unless $known{$k};
		croak "feature_descriptions entry for '$k' must be a plain string"
			if ref $fd->{$k};
	}
	return 1;
} ## end sub _validate_feature_descriptions

# Compile a munger spec against the model's feature names.  Requires
# Algorithm::ToNumberMunger on demand -- it is an optional dependency --
# and lets its compile() croak on any spec problem.
#
# Args:
#   $tags :: the model's feature_names, an arrayref of strings in column
#            order.  Must be non-empty.
#   $mungers :: the munger spec hashref, shaped as
#               Algorithm::ToNumberMunger->compile expects it.
#
# Returns: the compiled plan object.  Croaks when feature_names are
# missing, when the module cannot be loaded, or when compile() rejects the
# spec.
#
# Example:
#   my $plan = _compile_mungers( [ 'method', 'path_len' ],
#       { method => { munger => 'http_method_enum', default => -1 } } );
sub _compile_mungers {
	my ( $tags, $mungers ) = @_;
	croak "this model has mungers but no feature_names to compile them against"

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

# ---------------------------------------------------------------------------

# Returns the training data to actually build trees on, after applying the
# missing-value strategy.
#
# Args:
#   $data :: the caller's raw training set, an arrayref of feature-value
#            arrayrefs.  Never modified -- zero/impute build a copy rather
#            than densifying in place.
#
# Returns: the arrayref of rows to build trees from.  That is $data itself
# under die (already dense) and nan (undef is meaningful there), and a
# dense copy under zero/impute -- except on the C tree-building path, which
# fills from missing_fill itself and so also gets $data through unchanged.
# Under die it croaks instead, naming the first missing cell's row and
# column.  Under impute it also sets $self->{missing_fill}.
#
# Example:
#   my $train = $self->_prepare_fit_data($data);
#   my $tree  = $self->_build_tree( _subsample( $train, $psi ), 0, $limit );
sub _prepare_fit_data {
	my ( $self, $data ) = @_;
	my $m  = $self->{missing};
	my $nf = $self->{n_features};

	if ( $m eq 'die' ) {

		# Locate the first missing cell, then report it.  The scan is over
		# 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;
	}



( run in 2.590 seconds using v1.01-cache-2.11-cpan-d01c6094234 )