view release on metacpan or  search on metacpan

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

    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;
    }
    memcpy(b->idx + b->n_idx, idx, (size_t)n * sizeof(int));
    memcpy(b->val + b->n_val, val, (size_t)n * sizeof(double));
    b->n_idx += n;
    b->n_val += n;
    return off;
}

/* splitmix64 -- fast, well-mixed, and per-stream state fits in one
 * uint64_t, which is all a thread-private PRNG needs here. Not
 * cryptographic; doesn't need to be. */
static uint64_t sm64_next(uint64_t *s) {
    uint64_t z = (*s += 0x9E3779B97F4A7C15ULL);
    z = (z ^ (z >> 30)) * 0xBF58476D1CE4E5B9ULL;
    z = (z ^ (z >> 27)) * 0x94D049BB133111EBULL;
    return z ^ (z >> 31);
}

static double sm64_drand(uint64_t *s) {
    return (double)(sm64_next(s) >> 11) * (1.0 / 9007199254740992.0);
}

static double _ts_randn(uint64_t *s) {
    double u1 = sm64_drand(s);
    double u2;
    if (u1 == 0.0) u1 = 1e-12;
    u2 = sm64_drand(s);
    return sqrt(-2.0 * log(u1)) * cos(6.283185307179586 * u2);
}

/* Thread-safe twin of _build_node_c: same split algorithm, but reads
 * randomness from a thread-private splitmix64 stream instead of
 * Drand01(), and writes into a TreeBuf instead of allocating Perl AVs
 * -- so it touches no interpreter-global state and is safe to call
 * concurrently from an OpenMP parallel region, one tree per thread. */
static int _build_node_packed(const double* x, int nf, int* idxs, int size,
                               int depth, int limit, int mode_flag,
                               int ext_active, TreeBuf *buf, uint64_t *rng) {
    double *lo, *hi;
    int *varying, nv, f, my_idx;

    if (depth >= limit || size <= 1) {
        my_idx = tb_push_node(buf, 0.0, (double)size, 0.0, 0.0, 0.0, 0.0);
        free(idxs);
        return my_idx;
    }

    lo = (double*)malloc(nf * sizeof(double));
    hi = (double*)malloc(nf * sizeof(double));
    for (f = 0; f < nf; f++) {
        lo[f] = HUGE_VAL;
        hi[f] = -HUGE_VAL;
    }
    for (int i = 0; i < size; i++) {
        const double* row = x + (size_t)idxs[i] * (size_t)nf;
        /* See the matching comment in _build_node_c: no isnan() guard
         * needed, since NaN < x / NaN > x are always false already --
         * that's what lets this vectorize as a plain min/max scan.
         * omp simd here is thread-safe to call from inside the caller's
         * omp parallel region: it's a per-thread vectorization hint,
         * not a team construct, so it doesn't nest into anything. */
        #ifdef _OPENMP
        #pragma omp simd
        #endif
        for (int f2 = 0; f2 < nf; f2++) {
            double v = row[f2];
            if (v < lo[f2]) lo[f2] = v;
            if (v > hi[f2]) hi[f2] = v;
        }
    }

    varying = (int*)malloc(nf * sizeof(int));
    nv      = 0;
    for (f = 0; f < nf; f++) {
        if (lo[f] < hi[f]) varying[nv++] = f;
    }

    if (nv == 0) {
        free(lo); free(hi); free(varying);
        my_idx = tb_push_node(buf, 0.0, (double)size, 0.0, 0.0, 0.0, 0.0);
        free(idxs);
        return my_idx;
    }

    if (mode_flag == 0) {
        int attr     = varying[(int)(sm64_drand(rng) * nv)];
        double split = lo[attr] + sm64_drand(rng) * (hi[attr] - lo[attr]);
        int *lidx = (int*)malloc(size * sizeof(int));
        int *ridx = (int*)malloc(size * sizeof(int));
        int ln = 0, rn = 0, i, li, ri;

        for (i = 0; i < size; i++) {
            int row  = idxs[i];
            double v = x[(size_t)row * (size_t)nf + attr];

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

# rather than the full file -- the only fill bounded-memory fitting can form.
#
# Args:
#   $pool :: hashref of row index => arrayref of cells, as returned by
#            _gather_stream or _gather_indexed.  Rewritten in place.
#
# Returns: nothing.  Under zero/impute every row in $pool comes back dense;
# under impute the learned fill is also stored in $self->{missing_fill}.
# die and nan return untouched (die data is already dense, nan wants its
# undefs kept).
#
# Example:
#   my $pool = $self->_gather_indexed( $path, $offsets, $want );
#   $self->_apply_missing_to_pool($pool);   # undef cells now carry the fill
sub _apply_missing_to_pool {
	my ( $self, $pool ) = @_;
	my $m = $self->{missing};
	return if $m eq 'die' || $m eq 'nan';    # die: already dense; nan: keep undef

	my $nf = $self->{n_features};
	my $fill;
	if ( $m eq 'impute' ) {
		$fill = $self->_compute_impute_fill( [ values %$pool ] );
		$self->{missing_fill} = $fill;
		delete $self->{_fill_packed};
	} else {                                 # zero
		$fill = [ (0) x $nf ];
	}
	for my $i ( keys %$pool ) {
		my $r = $pool->{$i};
		$pool->{$i} = [ map { defined $r->[$_] ? $r->[$_] : $fill->[$_] } 0 .. $nf - 1 ];
	}
	return;
} ## end sub _apply_missing_to_pool

# Streaming counterpart of _learn_contamination_threshold: learn the exact
# score cutoff for the contamination rate without holding every score.  k+1
# largest scores are enough for _threshold_from_ranked's boundary logic; a
# min-heap keeps them in one scoring pass, and (only under a boundary tie) one
# extra pass resolves the tie block's edges.
#
# Args:
#   $path :: path to the CSV file the model was just fitted from.
#   $skip_first :: true to drop the first non-blank line as a header.
#   $n :: the census row count, which fixes how many rows k stands for.
#   $c_scan :: true to let the C packer coerce cells, skipping the Perl
#              numeric validation on the scoring passes.  See _stream_scores.
#
# Returns: nothing.  Sets $self->{threshold} to the same value the in-RAM
# _learn_contamination_threshold would have produced for this data.
#
# Example:
#   $self->_learn_contamination_threshold_streaming( 'train.csv', 1, 1_000_000, 1 );
sub _learn_contamination_threshold_streaming {
	my ( $self, $path, $skip_first, $n, $c_scan ) = @_;

	my $k = int( $self->{contamination} * $n + 0.5 );
	$k = 1  if $k < 1;
	$k = $n if $k > $n;

	# Whole set flagged: sit the cut just below the global minimum score.
	if ( $k >= $n ) {
		my $min;
		$self->_stream_scores( $path, $skip_first, sub { $min = $_[0] if !defined $min || $_[0] < $min }, $c_scan );
		$self->{threshold} = $min - 1e-9;
		return;
	}

	# One scoring pass keeps the k+1 largest scores (a min-heap rooted at the
	# smallest kept).  contamination <= 0.5 bounds k at n/2, so this tail is
	# always the smaller side of the split.
	my @heap;
	my $cap = $k + 1;
	$self->_stream_scores(
		$path,
		$skip_first,
		sub {
			my $s = $_[0];
			if    ( @heap < $cap )  { _heap_push( \@heap, $s ) }
			elsif ( $s > $heap[0] ) { _heap_replace_root( \@heap, $s ) }
		},
		$c_scan
	);

	my @desc = sort { $b <=> $a } @heap;    # k+1 largest, descending
	my $v    = $desc[ $k - 1 ];             # k-th largest
	my $lo   = $desc[$k];                   # (k+1)-th largest

	# Clean gap at the boundary: cut midway, exactly as _threshold_from_ranked.
	if ( $lo < $v ) {
		$self->{threshold} = ( $v + $lo ) / 2.0;
		return;
	}

	# A tie block of value $v straddles rank k.  Reproduce _threshold_from_ranked's
	# tie branch: locate the block's edges (i = first index at $v, j = first
	# index below it) and the neighbouring scores with one more pass.
	my ( $cnt_gt, $cnt_eq, $above, $below ) = ( 0, 0, undef, undef );
	$self->_stream_scores(
		$path,
		$skip_first,
		sub {
			my $s = $_[0];
			if ( $s > $v ) {
				$cnt_gt++;
				$above = $s if !defined $above || $s < $above;    # smallest > $v
			} elsif ( $s == $v ) {
				$cnt_eq++;
			} else {
				$below = $s if !defined $below || $s > $below;    # largest < $v
			}
		},
		$c_scan
	);

	my $i = $cnt_gt;              # first rank holding $v
	my $j = $cnt_gt + $cnt_eq;    # first rank below $v
	if ( $i > 0 && ( $k - $i ) < ( $j - $k ) ) {
		$self->{threshold} = ( $above + $v ) / 2.0;    # exclude the block
	} elsif ( $j < $n ) {
		$self->{threshold} = ( $v + $below ) / 2.0;    # include the block