GreyCat Algebra Library

@library("algebra", "0.0.0");

Numerical computing and machine learning library for GreyCat. Provides statistical profiling, neural networks, clustering, signal processing, pattern detection, and polynomial regression — all running natively within the GreyCat runtime.

Features

• Statistical profiling — full multi-dimensional Gaussian analysis (GaussianND: min, max, avg, std, covariance, correlation) or the lean per-feature diagonal profiler (GaussianArray: same stats + scaling at O(N·F), no covariance matrix)
• Neural networks — regression, classification, and autoencoder architectures with Dense, Linear, LSTM layers
• Dimensionality reduction — PCA with automatic best-dimension detection
• Clustering — K-means with mini-batch support and meta-learning
• Signal processing — FFT, frequency analysis, low-pass filtering, extrapolation
• Pattern detection — Euclidean, DTW, FFT, and SAX-based time-series pattern matching
• Polynomial regression — curve fitting, prediction, and time-series compression
• Time-series decomposition — aggregate instant data into hourly, daily, weekly, monthly, yearly
• Kernel density estimation — 1D Gaussian KDE with scipy parity
• Climate — UTCI (Universal Thermal Climate Index) calculation

Statistical Profiling — GaussianND

Learn statistical properties from multi-dimensional data and apply normalization transforms.

// Create a profiler and learn from data
var profile = GaussianND {};
var data = Tensor {};
data.init(TensorType::f64, Array<int> { 0, 5 });  // [batch, 5 features]
data.append([0.67, -0.20, 0.19, -1.06, 0.46]);
data.append([-0.20, 3.82, -0.13, 1.06, -0.48]);
// ... add more observations
profile.learn(data);

// Access statistics
var avg = profile.avg();          // [5] averages
var std = profile.std();          // [5] standard deviations
var cov = profile.covariance();   // [5x5] covariance matrix
var corr = profile.correlation(); // [5x5] correlation matrix

// Normalize data
var normalized = profile.min_max_scaling(data);         // (x - min) / (max - min)
var original = profile.inverse_min_max_scaling(normalized);

var standardized = profile.standard_scaling(data);      // (x - avg) / std
var restored = profile.inverse_standard_scaling(standardized);

// Crop to subset of features
var sub_profile = profile.crop(0, 2);  // features 0 to 2

API Reference — GaussianND

Method	Description
learn(input)	Learn from a [batch x N] tensor
avg()	Returns [N] tensor of dimension averages
std()	Returns [N] tensor of standard deviations
covariance()	Returns [N x N] covariance matrix
correlation()	Returns [N x N] correlation matrix
dimensions()	Returns N (number of dimensions)
clear()	Reset all state
min_max_scaling(input)	Min-max normalization
inverse_min_max_scaling(input)	Inverse min-max normalization
standard_scaling(input)	Standard scaling (z-score)
inverse_standard_scaling(input)	Inverse standard scaling
crop(from, to)	Create sub-profile with feature subset

Per-feature Profiling — GaussianArray

GaussianArray is the lean, per-feature (diagonal) counterpart of GaussianND: an array of N independent 1-D Gaussians learned in one fused pass. It keeps exactly the per-feature statistics — min / max / sum / sum_square (each [N]) — and the same scaling transforms, but none of GaussianND’s [N, N] cross-feature machinery.

Use GaussianND when you actually need the covariance / correlation matrix or PCA. Use GaussianArray for everything else — per-feature avg / std, min-max / standard scaling, or feeding k-means. The difference is asymptotic, not constant: GaussianND.learn accumulates the full [N, N] outer product Σ xᵢxⱼ in O(batch·N²) time and N² memory, whereas GaussianArray.learn is O(batch·N) / O(N). At N = 2000 features that is the gap between ~2·10⁶ multiply-adds + 32 MB per learned batch and a couple thousand adds + 64 KB — the reason k-means profiles the feature space with GaussianArray (Kmeans::featureProfile).

var profile = GaussianArray {};
profile.learn(data);                              // [batch, N] — one fused O(batch·N) pass
var avg = profile.avg();                          // [N] averages
var std = profile.std();                          // [N] std (ddof=1)
var scaled = profile.standard_scaling(data);      // (x - avg) / std
var back = profile.inverse_standard_scaling(scaled);

API Reference — GaussianArray

Method / field	Description
total / min / max / sum / sum_square	Accumulated count + per-feature [N] stats (nullable before the first learn)
learn(input)	Learn from a [batch x N] (or 1-D [N]) tensor; O(batch·N), no [N,N] matrix. Accumulates.
avg() / std()	[N] per-feature average / standard deviation (ddof = 1)
dimensions()	Returns N
clear()	Reset all accumulators (reusable)
min_max_scaling / inverse_min_max_scaling	Min-max normalization and its inverse
standard_scaling / inverse_standard_scaling	Standard (z-score) scaling and its inverse

Dimensionality Reduction — PCA

Identify the most important dimensions using Principal Component Analysis.

// Learn PCA from a GaussianND profile
var profile = GaussianND {};
profile.learn(data);

var pca = PCA {};
pca.learn(profile.correlation()!!, profile.avg()!!, profile.std()!!, 0.95);
// pca.best_dimension now holds the number of dimensions retaining 95% variance

// Set target dimensionality and transform
pca.set_dimension(pca.best_dimension!!);
var reduced = pca.transform(data);         // [batch x dim] → [batch x best_dim]
var reconstructed = pca.inverse_transform(reduced);  // back to original space

API Reference — PCA

Method	Description
learn(correlation, avg, std, threshold?)	Learn eigenvectors from correlation matrix. Threshold (default 0.95) sets variance retention
set_dimension(dim)	Set number of output dimensions
transform(input)	Project from N to dim dimensions
inverse_transform(input)	Project back from dim to N dimensions
get_dimension(threshold)	Get number of dimensions for a given variance threshold

Neural Networks

High-level API for building, training, and evaluating neural networks.

Regression Network

var inputs = 7;
var outputs = 2;

// Create network
var nn = RegressionNetwork::new(
    inputs, outputs, TensorType::f64,
    false,  // inputs_gradients
    0,      // fixed_batch_size (0 = dynamic)
    42,     // seed
);

// Optional preprocessing
nn.setPreProcess(PreProcessType::standard_scaling, inputProfile);
nn.setPostProcess(PostProcessType::standard_scaling, outputProfile);

// Add layers
nn.addDenseLayer(5, true, ComputeActivationRelu {}, null);
nn.addDenseLayer(3, true, ComputeActivationSigmoid {}, null);
nn.addDenseLayer(outputs, true, ComputeActivationRelu {}, null);

// Configure loss and optimizer
nn.setLoss(ComputeRegressionLoss::square, ComputeReduction::auto);
nn.setOptimizer(ComputeOptimizerAdam {});

// Build and compile
var engine = ComputeEngine {};
var model = nn.build(true);
var batchSize = nn.initWithBatch(model, engine, null, batch);

// Training loop
for (var epoch = 0; epoch < 100; epoch++) {
    var inputTensor = nn.getInput(engine);
    var targetTensor = nn.getTarget(engine);
    // fill inputTensor and targetTensor with data...
    var loss = nn.train(engine);

    // Validation
    var valLoss = nn.validation(engine);
}

// Prediction
nn.getInput(engine)?.fill(newData);
var prediction = nn.predict(engine);

Classification Network

var inputs = 10;
var classes = 3;

var nn = ClassificationNetwork::new(
    inputs, classes, TensorType::f64,
    false,  // inputs_gradients
    0,      // fixed_batch_size
    42,     // seed
    true,   // calculate_probabilities
    true,   // from_logits
    false,  // has_class_weights
);

nn.addDenseLayer(5, true, ComputeActivationRelu {}, null);
nn.addDenseLayer(classes, true, ComputeActivationSigmoid {}, null);

nn.setLoss(ComputeClassificationLoss::sparse_categorical_cross_entropy, null);
nn.setOptimizer(ComputeOptimizerSgd {});

AutoEncoder Network

var nn = AutoEncoderNetwork::new(
    inputs, TensorType::f64,
    false, 0, 42,
);

// Encoder layers
nn.addDenseLayer(64, true, ComputeActivationRelu {}, null);
nn.addDenseLayer(16, true, ComputeActivationRelu {}, null);  // bottleneck
// Decoder layers
nn.addDenseLayer(64, true, ComputeActivationRelu {}, null);
nn.addDenseLayer(inputs, true, ComputeActivationSigmoid {}, null);

nn.setEncoderLayer(1);  // bottleneck layer index
nn.setLoss(ComputeRegressionLoss::square, null);

LSTM Layers

Add LSTM layers for sequence modeling:

var nn = RegressionNetwork::new(inputs, outputs, TensorType::f64, false, 0, 42);

nn.addDenseLayer(5, true, ComputeActivationRelu {}, null);
nn.addLSTMLayer(
    6,      // output size
    3,      // number of stacked LSTM layers
    10,     // sequence length
    true,   // use_bias
    true,   // return_sequences
    true,   // bidirectional
    null,   // initializer config
);
nn.addLSTMLayer(3, 3, 10, true, false, false, null);  // last LSTM: return_sequences=false
nn.addDenseLayer(outputs, true, ComputeActivationRelu {}, null);

Available Components

Activations: Relu, LeakyRelu, Sigmoid, Tanh, Softmax, Softplus, SoftSign, Selu, Elu, Celu, HardSigmoid, Exp

Optimizers:

Optimizer	Description
ComputeOptimizerAdam	Adam (default, lr=0.001)
ComputeOptimizerSgd	Stochastic Gradient Descent (lr=0.01)
ComputeOptimizerRmsProp	RMSprop
ComputeOptimizerAdaDelta	Adadelta
ComputeOptimizerAdaGrad	Adagrad
ComputeOptimizerAdaMax	Adamax
ComputeOptimizerNadam	Nadam
ComputeOptimizerFtrl	FTRL
ComputeOptimizerMomentum	SGD with momentum
ComputeOptimizerNesterov	SGD with Nesterov momentum

Layer types: Linear, Dense, LSTM, Activation, Filter

Loss functions:

• Regression: square, abs
• Classification: categorical_cross_entropy, sparse_categorical_cross_entropy

Preprocessing: min_max_scaling, standard_scaling, pca_scaling

Weight initializers: xavier, xavier_uniform, relu, relu_uniform, lecun_uniform, normal, uniform, pytorch, identity, constant, and more

K-Means Clustering

Mini-batch K-means clustering with meta-learning for optimal initialization.

var batchSize = 10;
var features = 6;
var clusters = 3;
var tensorType = TensorType::f64;
var rounds = 10;

// Configure and compile
var engine = ComputeEngine {};
var model = Kmeans::configure(clusters, features, tensorType, true);
engine.compile(model, batchSize);

// Initialize
Kmeans::initialize(engine, 42);

// Training loop
for (var round = 0; round < rounds; round++) {
    Kmeans::init_round(engine);

    // Feed mini-batches
    for (var mb = 0; mb < numBatches; mb++) {
        Kmeans::learn(engine, miniBatches[mb]);
    }

    Kmeans::end_round(engine);
    var loss = Kmeans::getSumOfDistances(engine).get(Array<int> { 0 });
}

// Compute statistics
Kmeans::calculate_stats(engine);
var centroids = Kmeans::getClustersCentroids(engine);
var counts = Kmeans::getClustersCounts(engine);
var avgDist = Kmeans::getClustersAvgOfDistances(engine);
var interDist = Kmeans::getClustersDistancesToEachOther(engine);

// Inference on new data
var assignment = Kmeans::cluster(engine, newData);

API Reference — Kmeans

Method	Description
Kmeans::configure(clusters, features, type, stats)	Build compute model
Kmeans::initialize(engine, seed)	Initialize engine
Kmeans::init_round(engine)	Reset counters for new round
Kmeans::learn(engine, batch)	Train on a mini-batch
Kmeans::end_round(engine)	Finalize training round
Kmeans::calculate_stats(engine)	Compute cluster statistics
Kmeans::cluster(engine, batch)	Assign clusters to data
Kmeans::getClustersCentroids(engine)	Get centroid positions [K x F]
Kmeans::getClustersCounts(engine)	Get sample count per cluster
Kmeans::getSumOfDistances(engine)	Get total loss
Kmeans::getDistances(engine)	Distance of the last mini-batch to all centroids [batch, nbCluster]
Kmeans::getAssignment(engine)	Cluster assignment of the last mini-batch [batch]
Kmeans::getBestDistances(engine)	Distance of each sample to its closest centroid [batch]
Kmeans::getClustersSumOfDistances(engine)	Sum of distances within each cluster [nbCluster]
Kmeans::getClustersAvgOfDistances(engine)	Get avg distance per cluster
Kmeans::getClustersDistancesToEachOther(engine)	Get inter-cluster distances [K x K]
Kmeans::sortClusters(engine)	Sort clusters deterministically
Kmeans::getInferenceEngine(result, batchSize, stats)	Create inference engine from trained result

High-level training

For most use cases, drive training through the high-level entry points rather than the manual loop above. They build and compile the engine, seed the centroids, run the rounds, and return a KmeansResult. The input data may be a single Tensor (treated as one batch), an Array<Tensor>, a graph-backed nodeList<Tensor> of mini-batches, or a GcbReader<Tensor> streaming mini-batches from a .gcb file one tensor at a time (see the note on tensors below).

All per-run knobs are bundled into a single KmeansConfig record (built with KmeansConfig { nb_clusters: …, nb_rounds: …, seed: … } or KmeansConfig::of(k, rounds, seed), then mutate optional fields):

Method	Description
Kmeans::learning(tensors, gaussian, config, reportIdx)	Train once over the input
Kmeans::single_learning(tensor, gaussian, config)	Train once over a single tensor
Kmeans::meta_learning(tensors, config, metaRounds, parallel)	Best-of-N independent runs (optionally parallel)
Kmeans::meta_meta_learning(tensors, config, minClusters, maxClusters, step, stopRatio, metaRounds, parallel)	Elbow-method search for the number of clusters

gaussian above (and on Kmeans::score below) is a GaussianArray? — the lean per-feature (diagonal) profiler, not GaussianND: k-means only needs per-feature min/max/avg/std and the sum/sumSq for the variance metrics, never the [F,F] covariance matrix, so GaussianArray computes the profile in O(N·F) instead of GaussianND’s O(N·F²). It is nullable: pass null and learning/single_learning/score compute it locally via Kmeans::featureProfile(tensors) (one extra O(N·F) pass). Pass a precomputed one instead when you call these repeatedly over the same data — e.g. train then score the same tensors — to avoid recomputing it (meta_learning / meta_meta_learning already do this internally across restarts / the K-sweep, which is why they don’t expose the parameter at all).

KmeansConfig fields (only nb_clusters / nb_rounds / seed are required; every other field is nullable and resolves to the noted default):

Field	Default	Meaning
nb_clusters	—	Number of clusters (ignored by meta_meta_learning, which sweeps K)
nb_rounds	—	Max Lloyd iterations per run
seed	—	Base random seed (incremented per restart by the meta entry points)
init_mode	kmeansPlusPlus (single) / randomized (meta)	Centroid initialization recipe (KmeansMode)
replace_mode	farthestPoint (single) / randomized (meta)	Empty-cluster replacement recipe (KmeansMode)
seed_strategy	direct	Initial-seeding scalability strategy (KmeansSeedStrategy)
subsample_size	max(k*40, ceil(sqrt(N)))	Subsample rows for seed_strategy == subsample
kmeans_parallel_rounds	clamp(ceil(log2 N), 1, 8)	Oversample rounds for seed_strategy == kMeansParallel
kmeans_parallel_oversample	2.0	Oversample factor l = oversample * k for kMeansParallel
metrics	minimal	Which quality metrics to compute, as ONE ordered cost tier (KmeansMetrics) — replaces the former calculate_* booleans. See the tier table below.
max_engine_memory	null	Optional byte budget for the compute-engine arena: when set, learning estimates the arena (dominated by the always-f64 [batch, K] distance matrix) before compiling and throws with the largest batch size that fits, instead of OOM-ing the worker
tol	0.0	Lloyd early-stop tolerance (<= 0 runs all nb_rounds)
distance_metric	euclidean	Distance metric (KmeansDistance): euclidean or cosine (spherical) — see Cosine / shape clustering below
report	false	Emit per-restart live-progress files
report_prefix	""	Filename prefix for the live report files + cancel flag, so concurrent runs don’t clobber each other’s singleton paths — see Live progress & cancellation
collect_per_k	false	(meta_meta_learning only) retain every swept K’s best result in KmeansMetaResult.perK — see Per-K browsing
progress	null	In-process progress observer for the sequential paths (KmeansProgress) — the callback alternative to the report files; see In-process progress (KmeansProgress)
final_repair	true	Repair clusters left empty by the FINAL assignment: up to two extra replace+reassign passes after the Lloyd loop (strictly loss-non-increasing; stops as soon as a pass improves neither the empty count nor the loss). Off-data replace modes (randomUniform/randomNormal) are swapped for farthestPoint during repair. Each pass appends to roundsLoss/centroidShift and bumps nIter/emptyReplacements, so nIter can exceed nb_rounds by up to 2. Skipped when nb_rounds == 0 and after a cooperative cancel. Set false for the historical no-repair behavior
compact_empty_clusters	false	Return the result WITHOUT still-empty clusters: drops every clusters_count == 0 row from all per-cluster tensors, remaps assignment ids and updates nbClusters (see Kmeans::compactEmptyClusters). Only matters when empties survive the Lloyd loop + final repair (e.g. K > distinct points)

A null init_mode / replace_mode means unpinned: the single-run entry points resolve it to the fixed default above, while meta_learning / meta_meta_learning randomize it per restart — the init slot over the three cheap seeding modes (fromInput / randomUniform / randomNormal; pin init_mode = kmeansPlusPlus to pay its K seeding passes on every restart), the replace slot over all six modes (Kmeans::nextRandomReplaceMode), so the strong on-data replacements (farthestPoint / D²-resample / splitLargest) participate in unpinned meta runs. KmeansSeedStrategy (direct / subsample / kMeansParallel) and the splitLargest replace mode are detailed under Scalable seeding below.

KmeansMetrics — one ordered cost tier for quality metrics. Each tier is a superset of the cheaper ones, so config.metrics is a single “how much am I willing to pay for diagnostics” dial — the right axis at large N/K/F, where every step up is an order of magnitude more work. loss, inertia, and the per-cluster inertia / max distance / distance std are always computed (they only reuse tensors the engine already produced); the tiers add the metrics that cost extra passes:

Tier	Adds on top of the previous	Extra cost
minimal (default)	— (loss, inertia, per-cluster inertia/max/std only)	none
variance	explainedVariance, calinskiHarabasz	the two scalars are O(F) reads of the per-feature sum/sumSq profile; learning builds that GaussianArray regardless of tier (seeding needs it), so this tier is nearly free there — only Kmeans::score (with a null gaussian) pays the O(N·F) profile pass from this tier up
silhouette	mean + per-cluster simplified silhouette	one O(N·K) assignment pass; no quadratic-in-K memory
full	[K,K] inter-centroid distances + daviesBouldin + per-cluster per-feature variance/MAE	always-f64 [K,K] variable (K²·8 B) + O(K²·F) pass + O(K²) DB loop + two [K,F] f64 tensors (2·K·F·8 B). Keep OFF at large K.

The tensors argument is any and accepts a Tensor, an Array<Tensor>, a nodeList<Tensor>, or a GcbReader<Tensor>. For parallel == true meta-learning, pass a graph-backed nodeList<Tensor> so the mini-batches can be shared across job workers (a raw Array<Tensor> is copied into every worker).

Choosing the input carrier. The three multi-batch carriers differ only in where the batches live — the numerics are identical:

Carrier	Where the batches live	Peak resident memory	Use when
Array<Tensor>	all in RAM at once	the whole dataset	it comfortably fits in memory
nodeList<Tensor>	persisted in the graph (gcdata/)	one batch at a time (paged)	you want the batches stored/reused, or need parallel meta-learning (batches shared across workers)
GcbReader<Tensor>	a flat .gcb file on disk	one batch + the [K,F] centroids	the data is too big for RAM and you don’t want to persist it in the graph

Streaming from a .gcb file (GcbReader<Tensor>) — big clustering without the memory bill

This is the key lever for large-scale clustering: you can cluster a dataset far bigger than RAM without ever materializing an Array<Tensor> in memory and without persisting it into a nodeList in the graph. Write the mini-batches once to a .gcb with a GcbWriter<Tensor>, then hand learning (or any entry point) a GcbReader<Tensor> over that file:

// 1. Write the mini-batches to a flat .gcb ONCE (e.g. straight out of your import/ETL,
//    one [batchRows, nbFeatures] tensor at a time — nothing accumulates in memory).
var path = "${File::userDir()}kmeans_input.gcb";
var w = GcbWriter<Tensor> { path: path, append: false };
for (var b = 0; b < numBatches; b++) {
    w.write(buildMiniBatch(b));          // pre-scaled / pre-L2-normalized if needed (see below)
}
w.flush();

// 2. Stream it through training. The reader is rewound (pos = 0) each pass, so the same
//    file feeds seeding, every Lloyd round, and the stats passes.
var reader = GcbReader<Tensor> { path: path };
var cfg = KmeansConfig::of(k, rounds, seed);
var result = Kmeans::learning(reader, Kmeans::featureProfile(reader), cfg, 0);

The training paths iterate the reader forward one tensor at a time (never holding more than a single mini-batch plus the [K,F] centroids), so peak memory is set by your batch size, not by N. Because k-means is multi-pass, the reader is rewound (pos = 0) at the start of every pass — seeding, each Lloyd round, and the stats passes — so the algorithm re-streams the file each pass. That trades memory for I/O: the file is re-read once per pass (× restarts × the K-sweep for the meta entry points). Results are bit-identical to the same batches passed as an Array/nodeList (locked by the kmeans_gcbreader_learning_matches_array test).

Two constraints follow from the reader being forward-only and read-only:

1. The in-place preprocessing helpers (Kmeans::scaleBatches / Kmeans::sphericalizeBatches / Kmeans::setBatch) do not accept a reader — there is no batch slot to write back. Pre-scale and (for cosine) pre-L2-normalize the tensors before writing them to the .gcb.
2. For parallel == true meta-learning give each job its own reader (a shared pos would race across workers) — or just prefer a nodeList<Tensor> there, which is built for sharing.

Sizing at large K. The engine compiles to the LARGEST mini-batch, and its dominant buffer is the [batch, K] distance matrix, which is always f64 regardless of the model dtype — so the arena is roughly batch·(8·K + F·elemSize + 16) + 2·K·F·elemSize bytes (f32 halves the data- and centroid-scaled terms only). Pick the batch size from your memory budget (e.g. at K = 20000 a 10000-row batch costs ≈ 1.6 GB of distance matrix alone), and set max_engine_memory to fail fast with the largest fitting batch size instead of OOM-ing. A single batched engine variable is hard capped at 2^31−1 elements, so batch · K must stay below that — mini-batch the input rather than passing one giant Tensor. All tensors stored on a returned KmeansResult are standalone copies (cloned at the result boundary), so holding results does not keep engine memory alive; the raw Kmeans::get* engine accessors, in contrast, return views into the engine’s arena — clone() those before storing them beyond the engine’s lifetime.

rounds is the maximum number of Lloyd iterations: each round that runs performs one real centroid update and a re-assignment and appends one entry to roundsLoss, so roundsLoss.size() == result.nIter — equal to rounds unless tol convergence or a cancellation stops the loop early (then nIter < rounds), or the default final_repair fills late empties (up to two extra passes, then nIter <= rounds + 2). The returned loss/centroids/assignment are always mutually consistent (rounds = 1 does one real update). tol is a convergence tolerance on the maximum per-row centroid shift ‖Δcentroid‖ (in feature-space units): training stops early once, on a round that did NOT relocate an empty cluster, the largest centroid moved less than tol. The non-squared-distance loss is not monotone under Lloyd (the mean update minimizes Σd², not Σd), so the centroid shift — which settles toward 0 as the means converge — is the stop signal, not a loss delta; pass tol <= 0.0 to always run all rounds. meta_meta_learning sweeps K from minClusters to maxClusters inclusive in increments of step (it ignores config.nb_clusters and runs a per-K clone(config)). The config.report flag (plus learning’s reportIdx argument) drives the opt-in live-progress files described in Live progress & cancellation below (report = false disables all of it; reportIdx is the per-restart file index used by meta_learning’s parallel restarts).

Search diagnostics & quality metrics. KmeansResult additionally exposes (all nullable): nbClusters (= centroids.shape()[0], a convenience copy of K so callers skip a shape read), seedStrategy (the seeding strategy actually used — kMeansParallel ignores initMode), nIter / emptyReplacements (process counters; final repair passes included), inertia + clustersInertia (Σ squared distance — the standard objective, distinct from loss which is the sum of non-squared distances), explainedVariance, calinskiHarabasz (both from the variance tier up), daviesBouldin (only at the full tier — it needs the [K,K] inter-centroid distances), silhouette (mean simplified silhouette over all samples; this centroid-based simplified form uses a = distance to the nearest centroid, so it is in [0, 1] and never goes negative like the classical pairwise silhouette — read it as a normalized separation margin and compare it relatively across K, not on the classical [-1, 1] scale; EMPTY clusters are excluded from its nearest-other search, mirroring daviesBouldin’s occupied-only rule; from the silhouette tier up), clustersSilhouette (per-cluster mean silhouette [K]; same tier), clustersMaxDistance (per-cluster max member distance / radius [K] — tightness and outlier signal, always set), clustersDistanceStd (per-cluster std of member distances [K], always set), clustersFeatureVariance (per-cluster per-feature variance [K, F] — mean squared deviation from the centroid; full tier, for both f32 and f64 models — see the precision note below), clustersFeatureMae (per-cluster per-feature mean absolute error [K, F]; same full tier), and centroidShift (per-iteration max Euclidean centroid movement ‖Δcentroid‖ [nIter] — always set; the natural convergence diagnostic behind a shift-based tol. It is measured AFTER empty-cluster replacement, so a pass that teleported an empty centroid records the large, REAL relocation — expect spikes on the passes counted by emptyReplacements; tol never stops on those passes). meta_meta_learning returns the full elbow trace on KmeansMetaResult: clusterCounts, clusterLosses (the Σ-distance curve, aligned), clusterInertias (the SSE curve), and the selected bestNbClusters. Both meta_learning and meta_meta_learning also populate restarts (nullable) — one KmeansRestartConvergence per restart of the kept K, each carrying that restart’s index, nIter, final loss, inertia, full per-round roundsLoss curve, and emptyClusters (the number of clusters still empty at convergence — with the default final_repair this is > 0 only when even the repair passes could not fill them, i.e. K outstrips the distinct/separable points; with final_repair = false a cluster starved late in the run can also stay empty; nullable for schema-safe evolution, always set on fresh results) — so callers can plot every restart’s convergence (and flag degenerate runs), not just the winning one. This count is populated for every mode (sequential meta, parallel meta, and auto-K), since all of them build their restart records through the same Kmeans::restartRecord collapse point. Kmeans::emptyClusterCount(counts: Tensor?) is the underlying helper (count of zero entries in a [K] clusters_count tensor; null → 0) if you want it on an arbitrary result.

Per-K browsing (collect_per_k). By default meta_meta_learning keeps only the elbow K’s clustering (bestResult) and discards every other swept K’s result. Set config.collect_per_k = true to retain them all: KmeansMetaResult.perK is then an Array<KmeansPerK> with one entry per swept K, in sweep order and aligned 1:1 with clusterCounts (including non-accepted K’s and the K the sweep stopped on). Each KmeansPerK carries k, loss, inertia, accepted (passed the elbow accept test), chosen (k == bestNbClusters — exactly one entry), and result — that K’s full KmeansResult (centroids, assignment, per-cluster stats, …), the same value meta_learning would return for that K. This lets a UI browse the best clustering at every K without re-running the sweep or re-implementing the elbow rule. It costs nothing extra to compute (the results are already produced), but each retained result holds O(N) tensors, so it is off by default and freed when the call returns; perK is null when the flag is off (and for plain meta_learning, which does not sweep K).

Cross-K seed contract. meta_meta_learning runs every swept K from the same base config.seed (the per-K clone(config) does not change the seed); inside each K, restart i uses seed + i * 1_000_003 (a wide prime stride that decorrelates restarts independently of nb_rounds). So restart i draws the same RNG span across all K, and a full sequential (parallel == false) run is deterministic for a given (seed, data) — pinned by the kmeans_meta_meta_cross_k_seed_determinism test.

Selection criterion. Restart selection (the kept bestResult in meta_learning) and the meta_meta_learning elbow rank candidates by inertia (Σ squared distance, the SSE objective Lloyd minimizes), not by loss (the sum of non-squared distances, whose minimizer is the geometric median). The elbow stop is robust to non-monotone loss: it compares each K against the best SSE seen so far (not just the previous K) and never stops on a stochastic uptick (a larger K landing a worse local optimum), so a noisy curve does not produce a false elbow. The single-run / seed_strategy = direct numerics are unchanged; only multi-restart/elbow selection uses inertia.

Initialization modes (KmeansMode): fromInput (sample real points), randomUniform, randomNormal, and kmeansPlusPlus (D²-weighted seeding — higher-quality starts that scale to large K). Replacement modes (for empty clusters during training): farthestPoint, kmeansPlusPlus (D²-resample against the running distances), and splitLargest (split the highest-inertia cluster) — these are not used for initialization. fromInput / randomUniform / randomNormal work in both contexts.

Scalable seeding (KmeansSeedStrategy). At large N, sequential kmeansPlusPlus over the full dataset costs K sequential passes. config.seed_strategy stages the chosen init_mode recipe over less data:

Strategy	What it does	When
direct (default)	Run init_mode over all N rows (historical behavior).	Small/medium N.
subsample	Seed init_mode on a row-subsample S ≪ N (subsample_size), then run the full Lloyd loop over all data. When the resolved S >= N (e.g. the k*40 default at large k), the strategy degenerates and seeding falls through to direct (no bootstrap copy, no second feature pass).	Large N, cheap strong baseline.
kMeansParallel	k-means‖: oversample O(log N) rounds (kmeans_parallel_rounds), drawing ~l = oversample·k D²-weighted candidates per round, then reduce the candidate set to k by a weighted k-means++.	Very large N, parallel-friendly seeding quality close to sequential k-means++.

seed_strategy is orthogonal to init_mode: it picks over how much data / via what staging the recipe runs, not what a centroid is. Note this is distinct from training on a sample and then scoring the full population (a separate workflow), whereas subsample only accelerates the seeding of a full-data train.

Empty-cluster fill at scale. replace_mode = kmeansPlusPlus re-seeds each empty cluster by D²-weighted sampling against the running per-sample distances; replace_mode = splitLargest splits the highest-inertia cluster (placing the new centroid on its farthest member) — useful when there are no spare points (k > distinct samples).

Cosine / shape clustering (distance_metric = cosine). To cluster by direction (shape) rather than magnitude, set config.distance_metric = KmeansDistance::cosine and L2-normalize the input rows first with Kmeans::sphericalizeBatches(tensors) (call it on the training data and on anything you later score). This runs spherical k-means: rows lie on the unit sphere and learning re-projects the centroids onto the sphere after each update, so the existing Euclidean kernel’s argmin equals the cosine argmax. The reported distances / loss are chord distances between unit vectors (d = sqrt(2·(1 − cosθ))); inertia (Σd²) is the quantity that equals Σ 2·(1 − cosine similarity). result.centroids are unit-norm direction vectors.

Kmeans::sphericalizeBatches(batches);                 // L2-normalize rows once
var g = Kmeans::featureProfile(batches);             // gaussian over the unit rows
var cfg = KmeansConfig::of(k, rounds, seed);
cfg.distance_metric = KmeansDistance::cosine;
var r = Kmeans::single_learning(batches, g, cfg);

Representative member (medoid). Kmeans::nearestMemberPerCluster(result) returns, per cluster, the global row index of the member closest to the centroid — a real representative sample (e.g. a representative entity) to show alongside the mean centroid. Indices are into the global row order the batches were consumed in; -1 marks an empty cluster.

Scaling to very large N · F · K (e.g. 1M rows × 1000 features × 10000 clusters): the default seeding (init_mode = kmeansPlusPlus, seed_strategy = direct) costs O(K·N·F) before the Lloyd loop — set config.seed_strategy = subsample (or kMeansParallel) so seeding runs over far less data, and a positive config.tol so converged runs stop early. Keep config.metrics = minimal (or at most variance) at large K: the silhouette tier’s pass is O(N·K) and the full tier’s [K, K] inter-centroid matrix behind Davies–Bouldin is O(K²) work and O(K²) memory (~800 MB at K = 10000). Prefer f32 to halve data/centroid RAM (all metrics stay f64-accurate — see the precision note). The per-element work in seeding, empty-cluster fill, all metrics, and the cosine projection runs in native single-threaded kernels; the engine’s worker pool provides parallelism at the meta-learning (restart) level. When the dataset itself exceeds RAM, stream it with a GcbReader<Tensor> (see Streaming from a .gcb file above) so peak memory is one mini-batch instead of the full Array<Tensor> — and without persisting it into a nodeList in the graph.

Reusable helpers

The library owns the generic orchestration around training so a consumer never re-implements it; the project keeps feature extraction, the storage container, and the decision of which metrics to compute / whether to persist. None of these helpers persist anything.

Helper	What it does
Kmeans::score(trained, tensors, gaussian, metrics)	Assign every row of tensors to trained’s centroids and (re)fill assignment / distances / loss, per-cluster count/sum/avg distance, and the metrics-tier quality metrics (KmeansMetrics) — in one call. Returns the loss. Use it to train on a sample and score the full population, or to re-score a loaded result on fresh data. gaussian (a GaussianArray?) must describe the scored population; pass null to have it computed locally over tensors when the tier needs it.
Kmeans::featureProfile(tensors)	Learn a GaussianArray (per-feature min/max/sum/sumSq, O(N·F)) over all mini-batches — the per-feature stats learning needs (no [F,F] covariance, unlike GaussianND). Accepts a GcbReader<Tensor> (streams it once, rewinding on open).
Kmeans::totalRows(tensors)	Total sample rows across all batches (streams a GcbReader<Tensor> forward, one batch at a time).
Kmeans::setBatch(tensors, i, t)	Write-back twin of batchAt (replace batch i); throws on a bare Tensor. Read-only carriers (a GcbReader<Tensor>) have no writable slot and are rejected — pre-transform before writing the .gcb.
Kmeans::recomputeClusterAggregates(result, k)	Recompute per-cluster count/sum/avg distance from a raw assignment + distances (used by score).
Kmeans::normalizeBatches(tensors)	Learn a GaussianArray and min-max-scale every batch IN PLACE to it; returns the GaussianArray so you can inverse_min_max_scaling the trained centroids back to raw space.
Kmeans::scaleBatches(tensors, g)	Min-max-scale every batch IN PLACE with an already-learned GaussianArray.
Kmeans::sphericalizeBatches(tensors)	L2-normalize every row IN PLACE (project onto the unit sphere) — the data-prep step for distance_metric == cosine.
Kmeans::nearestMemberPerCluster(result)	[K] i64 tensor of the global row index of the member nearest each centroid (approx. medoid / representative real sample; -1 for an empty cluster).
Kmeans::compactEmptyClusters(result)	Remove still-empty clusters from a finished result IN PLACE: drops their rows from every per-cluster tensor, remaps assignment ids, updates nbClusters; returns the number removed (0 = no-op). Called automatically when config.compact_empty_clusters is set; usable on loaded results too.
Kmeans::emptyClusterCount(counts)	Count of zero entries in a [K] clusters_count tensor (null → 0) — the “how many clusters ended empty” probe behind KmeansRestartConvergence.emptyClusters.

Train-on-a-sample, score-the-population (the canonical large-N flow): the project samples its own entities, builds the sample’s feature batches, trains, then scores the full population —

var sampleG = Kmeans::featureProfile(sampleBatches);
var trained = Kmeans::meta_learning(sampleBatches, cfg, metaRounds, parallel).bestResult!!;
// assign + score every row, refit the full metric tier over the full population:
var fullG = Kmeans::featureProfile(fullBatches);
Kmeans::score(trained, fullBatches, fullG, KmeansMetrics::full);

To normalize first, call Kmeans::normalizeBatches(batches) before training and g.inverse_min_max_scaling(trained.centroids!!) after — the project decides when and on which container (it already chose in-memory vs graph-backed).

Precision — f32 vs f64. The model’s compute dtype is taken from the input tensors (the tensor_type you pass to configure, or the dtype of the batches you hand to learning / meta_learning / meta_meta_learning). Both are fully supported end to end — k-means++ seeding, the Lloyd loop, empty-cluster replacement, and every quality metric run in either dtype. Choosing f32 halves the memory of the data and the centroids and speeds up the hot distance/accumulation loops; the trade-off is confined to the centroid coordinates themselves — result.centroids is f32, and the mean-update accumulates in f32, so centroid positions carry single-precision rounding (which compounds over very large N·F). Everything you measure stays double-precision regardless of the model dtype:

• result.assignment is always i64; result.distances, loss, inertia / clustersInertia, and the [batch, K] distance matrix are always f64 (the distance objective is computed and reduced in f64 even for an f32 model, so loss/inertia and k-means++ sampling are dtype-stable).
• All derived metrics — silhouette / clustersSilhouette, daviesBouldin, calinskiHarabasz, explainedVariance, clustersMaxDistance, clustersDistanceStd — are computed in f64.
• clustersFeatureVariance / clustersFeatureMae are produced for both dtypes: the native kernels read the batch/centroids at their element size but accumulate into f64 output tensors, so these are f64-typed and full-precision even when the model is f32.

In short: pick f32 for large N × F × K to halve data/centroid RAM and accelerate the distance pass — all reported metrics remain f64-accurate, and only the stored centroid coordinates are single-precision. If you normalize with GaussianArray (min_max_scaling / standard_scaling), learn the GaussianArray in the same dtype as the data: those methods output the profile’s dtype and reject a dtype mismatch, so a same-dtype profile keeps the whole pipeline in f32 (no silent upgrade to f64).

Live progress & cancellation

When report == true, the high-level entry points publish progress to JSON files under File::userDir() so a separate poller (e.g. an @exposed endpoint) can render a live convergence view of a long run. Passing report == false (the default for a one-shot call) disables all of it — no files are written. Each file is published atomically (write-to-.tmp then rename), so a reader never observes a half-written object.

The file names below are shown with the optional config.report_prefix token (<prefix>, empty by default). Set report_prefix to a unique per-run id (e.g. "run_<uuid>_") so several clusterings can stream concurrently without clobbering each other’s singleton files — the report layer threads the prefix through every path it touches, including the cancel flag. With an empty prefix the paths are the historical fixed names. report_prefix is a filename token, not a path: the files live directly under File::userDir(), so avoid / and ...

File	Type	Written by	Contents
<prefix>kmeans_restart_<i>.json	KmeansRestartFile	each restart of learning / single_learning	k, restart (= reportIdx), round, rounds, loss, done, and the partial roundsLoss curve so far (full once done)
<prefix>kmeans_progress.json	KmeansProgressControl	meta_meta_learning	overall sweep state: running / done / cancelled, k_min/k_max/k_step/k_cur, best_k (= the returned elbow bestNbClusters once finished), best_loss, the live k_values/k_losses elbow trace (rewritten after each K completes), and total
<prefix>kmeans_cancel.flag	(presence only)	the caller, to request cancellation	its existence asks the run to stop cooperatively at the next round / restart / K boundary

The live elbow trace needs no separate file: kmeans_progress.json already carries k_values / k_losses and is rewritten after each K completes, so a poller can animate the elbow during the sweep.

KmeansProgressControl.total is the planned step count (nK · restarts · rounds) — an UPPER BOUND. tol early-stop, the elbow stop, and cancellation all run fewer steps, so drive a progress bar off the done / cancelled flags rather than expecting steps / total to reach 100%. To cancel, create <prefix>kmeans_cancel.flag in File::userDir(); the run deletes it on exit and sets cancelled = true. The top-level entry points (single_learning, meta_learning, meta_meta_learning) also purge a stale flag at start, so a flag left over from a previous cancelled run can never silently truncate the next one (learning() itself only polls the flag — clear it yourself when driving learning() directly with report). Cancellation is checked at fine granularity — per Lloyd round inside a restart, per restart in meta_learning, and per K in meta_meta_learning (all only while report == true). For parallel meta_learning, each restart writes its own <prefix>kmeans_restart_<reportIdx>.json, so the per-restart files are parallel-safe; read exactly restarts_total of them.

In-process progress (KmeansProgress)

For an in-process consumer (same runtime — a Task, a service), the JSON files are unnecessary round-trips. Set config.progress to a KmeansProgress subtype and the library calls it directly. GCL has no method override, so every method is abstract — implement all five (make the ones you don’t need no-ops; cancelled returns false):

type MyProgress extends KmeansProgress {
    fn onRestart(restart: int, result: KmeansResult, conv: KmeansRestartConvergence, isBest: bool) {
        if (isBest) { /* publish a best-so-far snapshot the moment it improves */ }
    }
    fn onRound(restart: int, round: int, loss: float) { /* drive a fine progress bar */ }
    fn onBatch(restart: int, round: int, batchIndex: int, batchCount: int) { /* one mini-batch assigned */ }
    fn onK(k: int, best: KmeansResult, accepted: bool) { /* one swept K finished */ }
    fn cancelled(): bool { return false; } // return true to stop cooperatively
}

var cfg = KmeansConfig::of(8, 20, 42);
cfg.progress = MyProgress {};
var meta = Kmeans::meta_learning(batches, cfg, 5, false);

When cancelled, onRound fires per Lloyd iteration (sequential paths only), onBatch fires per mini-batch of every assignment pass — both the initial pre-loop assignment (round == 0) and every Lloyd round’s re-assignment (same restriction as onRound) — so a single pass over many mini-batches still reports smooth progress instead of jumping only between whole rounds, onRestart once per finished restart in restart order (both sequential and — after the join — parallel meta_learning; parallel restarts run in separate Job workers that cannot call back, so they get no live onRound/onBatch), and onK once per swept K from meta_meta_learning. cancelled() is polled between rounds, restarts, and Ks (not between individual mini-batches). The observer is never serialized to parallel workers or persisted. It is fully orthogonal to the report files — use either, both, or neither.

Signal Processing — FFT

Fast Fourier Transform for frequency analysis, filtering, and extrapolation.

var N = 1000;
var timeseries_complex = Tensor {};
timeseries_complex.init(TensorType::c128, Array<int> { N });

// Fill with a composite sine wave
var freq1 = 5.0;
var freq2 = 7.0;
var t = 0.0;
var dt = 1.0 / (freq1 * 200.0);
for (var i = 0; i < N; i++) {
    timeseries_complex.set(Array<int> { i }, sin(2 * MathConstants::pi * freq1 * t)
        + 0.3 * sin(2 * MathConstants::pi * freq2 * t));
    t = t + dt;
}

// Forward FFT: time → frequency domain
var frequency_complex = Tensor {};
var fft = FFT::new(N, false);
fft.transform(timeseries_complex, frequency_complex);

// Analyze frequency spectrum
var freq_table = FFT::get_frequency_table(frequency_complex, sampling_step);

// Apply low-pass filter
var filtered = Tensor {};
var cutoff = FFT::get_low_pass_filter_size(frequency_complex, 0.95);
FFT::apply_low_pass_filter(frequency_complex, filtered, cutoff);

// Inverse FFT: frequency → time domain
var fft_inv = FFT::new(N, true);
var reconstructed = Tensor {};
fft_inv.transform(reconstructed, frequency_complex);

// Extrapolation using frequency components
var value = FFT::extrapolate(frequency_complex, sampling_step, start_time, target_time, cutoff);

FFTModel — High-Level Time-Series Analysis

var model = FFTModel::train(myNodeTime, fromTime, toTime);

// Predict a single value
var predicted = model.extrapolate_value(futureTime, 0.95, null);

// Predict a range
var table = model.extrapolate(fromTime, toTime, 0.95, null, null);

API Reference — FFT

Method	Description
FFT::new(n, inverse)	Create FFT engine for n samples
fft.transform(time, freq)	Execute forward or inverse FFT
fft.transform_table(ts, time_c, freq_c)	Transform from Table, return frequency table
FFT::get_frequency_table(freq, step)	Get frequency analysis table
FFT::get_frequency_spectrum(freq, spec, db, filter)	Extract spectrum with optional dB conversion
FFT::apply_low_pass_filter(src, dst, cutoff)	Apply low-pass filter
FFT::get_low_pass_filter_size(freq, ratio)	Get cutoff for desired signal retention ratio
FFT::extrapolate(freq, step, start, t, filter)	Predict value at time t
FFT::extrapolate_table(time_c, step, start, from, to, skip)	Predict range of values
FFT::get_next_fast_size(n)	Get optimal FFT size >= n

Pattern Detection

Detect recurring patterns in time-series using multiple algorithms.

// Create time series
var ts = nodeTime<float> {};
for (var i = 0; i < 50; i++) {
    ts.setAt(time::new(i, DurationUnit::seconds), sin(MathConstants::pi * i / 10));
}

// Create detection engine (Euclidean, DTW, FFT, or SAX)
var engine = EuclideanPatternDetectionEngine::new(ts);
engine.state = PatternDetectionEngineState::new();

// Define reference patterns
engine.addPattern(
    time::new(10, DurationUnit::seconds),
    time::new(15, DurationUnit::seconds),
);

// Compute similarity scores
engine.initScoring();
engine.computeScores(null);

// Detect matches
engine.detect(PatternDetectionSensitivity {
    threshold: 0.0,   // minimum score threshold
    overlap: 1.0,     // allowed overlap ratio
}, null);

// Access results
for (timestamp, detection in engine.state.detections) {
    info("Match at ${timestamp}: score=${detection.score}, pattern=${detection.best_pattern}");
}

Available Detectors

Detector	Description
EuclideanPatternDetectionEngine	Euclidean distance-based matching
DTWPatternDetectionEngine	Dynamic Time Warping
FFTPatternDetectionEngine	FFT-based frequency matching
SaxPatternDetectionEngine	Symbolic Aggregate Approximation
RandomPatternDetectionEngine	Random baseline (for benchmarking)

Normalization Modes

Mode	Description
as_is	No normalization
shift	Vertical shift alignment
scaling	Vertical scaling alignment
shift_and_scaling	Both shift and scaling

Polynomial Regression

Fit polynomial curves for regression and time-series compression.

var N = 6;
var degree = 3;

var X = Tensor {};
X.init(TensorType::f64, Array<int> { N });
var Y = Tensor {};
Y.init(TensorType::f64, Array<int> { N });

for (var i = 0; i < N; i++) {
    var x = i * 10.0 + 1000;
    X.set(Array<int> { i }, x);
    Y.set(Array<int> { i }, 53.0 - 0.0002 * x + 0.00001 * x * x);
}

// Fit polynomial
var poly = Polynomial {};
var maxError = poly.learn(degree, X, Y);

// Predict
var predictions = poly.predict(X);
var singleValue = poly.predictValue(1050.0);

Time-Series Compression

Polynomial::compress(originalTS, polynomialTS, 5, 0.01, 1000);
Polynomial::decompress(originalTS, polynomialTS, 0.01, decompressedTS, errorTS);

API Reference — Polynomial

Method	Description
poly.learn(degree, X, Y)	Fit polynomial of given degree. Returns max error
poly.predict(X)	Predict Y values for tensor X
poly.predictValue(x)	Predict single Y value
Polynomial::compress(src, dst, maxDeg, maxErr, bufSize)	Compress time-series with adaptive polynomial fitting
Polynomial::decompress(src, poly, maxErr, dst, errTS)	Decompress and verify

Linear Solver

var weights = Solver::solve(X, Y);  // Solve X * w = Y for w

Time-Series Decomposition

Aggregate instant-level data into coarser time resolutions.

TimeSeriesDecomposition::calculateAll(
    instantTS,   // source
    hourlyTS,    // hourly aggregation (nullable)
    dailyTS,     // daily aggregation (nullable)
    weeklyTS,    // weekly aggregation (nullable)
    monthlyTS,   // monthly aggregation (nullable)
    yearlyTS,    // yearly aggregation (nullable)
    TimeZone::Europe_Luxembourg,
    null,        // lastUpdatedTime (null = full recalculation)
);

Supports incremental updates by passing lastUpdatedTime to recompute only from that point forward.

Calendar Profiling

The recurring-slot counterpart of TimeSeriesDecomposition (which down-samples to coarser instants): fold a time series onto fixed recurring calendar slots (“typical week”, “typical day”, “month-of-year”), each summarized by a Gaussian (mean / std / min / max / count). Each fold walks the series once; an empty/absent series yields an empty result.

// 7 days x 96 quarter-hours = 672 slots; slot = dayOfWeek*96 + (hour*96/24 + minute*96/1440)
var week  = CalendarProfile::weekShape(series, TimeZone::"UTC", 96);  // nodeList<Gaussian>
var day   = CalendarProfile::dayShape(series, tz, 96);                // 96 slots, folded across all days
var month = CalendarProfile::monthOfYear(series, tz);                 // 12 slots, slot = month - 1
var avgs  = GaussianSlots::means(week, 672, 0.0);  // dense [672] of slot means (absent slot -> 0.0)

Method	Description
CalendarProfile::weekShape(series, tz, slotsPerDay)	7*slotsPerDay slots: day-of-week × within-day
CalendarProfile::dayShape(series, tz, slotsPerDay)	slotsPerDay slots: a typical day folded across all days
CalendarProfile::monthOfYear(series, tz)	12 slots, slot = month - 1
GaussianSlots::means(slots, n, absent) / stds(...)	Dense [n] Array<float> of per-slot mean / sample-std from a sparse nodeList<Gaussian> (one ordered pass; absent or empty slot → absent). Works for any element type incl. time/duration (derived from the float accumulators).

CalendarPivot is the 2D form (e.g. a day-of-week × hour-of-day heatmap). GCL has no function-typed parameters, so the two axis extractors are abstract methods — subclass and implement rowOf / nRows / colOf / nCols, then call pivot(series, tz). Each cell is a Gaussian, so an empty cell is just count == null:

type WeekHourPivot extends CalendarPivot {
    fn rowOf(t: time, tz: TimeZone): int { return t.dayOfWeek(tz); }
    fn nRows(): int { return 7; }
    fn colOf(t: time, tz: TimeZone): int { return Date::from_time(t, tz).hour; }
    fn nCols(): int { return 24; }
}
var grid = WeekHourPivot {}.pivot(series, tz);  // Array<Array<Gaussian>> [7][24]

Sampling & vector utilities

Generic data-prep helpers used across the clustering / profiling pipelines.

// Seeded, deterministic sampling without replacement (partial Fisher–Yates):
var idx  = Sampling::sampleIndices(n, k, seed);     // k distinct indices in [0, n) (k clamped to [0, n])
var pick = Sampling::sample(items, k, seed);        // k distinct elements (input not mutated)

// In-place / pure ops on a plain Array<float>:
var sums   = Vectors::foldSum(v, group);            // sum each contiguous run of `group` (remainder dropped)
var phase  = Vectors::foldMeanByPhase(v, cycleLen); // mean across cycles at each phase (length == cycleLen)
Vectors::l2NormalizeVector(v);                      // scale to unit L2 length in place (no-op on zero)
var total  = Vectors::normalizeBySum(v);            // project onto the simplex in place; returns original sum

Method	Description
Sampling::sampleIndices(n, k, seed)	k distinct indices in [0, n), deterministic for (n, k, seed)
Sampling::sample(items, k, seed)	k distinct elements of items (input unchanged)
Vectors::foldSum(v, group)	Sum each contiguous run of group values; output length v.size()/group
Vectors::foldMeanByPhase(v, cycleLength)	Mean across cycles at each phase; output length cycleLength
Vectors::l2NormalizeVector(v)	Scale v to unit L2 length in place (no-op on a zero vector)
Vectors::normalizeBySum(v)	Project v onto the probability simplex in place; returns the original sum

Kernel Density Estimation

Estimate a smooth 1D probability density from samples with a Gaussian kernel. Matches scipy.stats.gaussian_kde (single-variable) to within 1e-9.

static native fn Kde::evaluate(
    data: Tensor,        // f64 sample points (1D)
    eval_points: Tensor, // f64 points to evaluate the density at (1D)
    bw_factor: float?,   // null => Scott's rule; else the scalar bandwidth factor
    bw_scale: float?,    // multiplies the Scott factor (default 1.0); ignored if bw_factor set
    max_samples: int?,   // null => use all data; else seeded subsample cap for the fit
    seed: int            // seed for the subsample selection (reproducible)
): Tensor;               // f64 density at each eval point, same length as eval_points

Bandwidth contract:

• bw_factor == null → Scott’s rule: factor = m^(-1/5) * bw_scale, where m is the fit size after any subsampling. bw_scale (default 1.0) scales the Scott factor.
• bw_factor != null → that value is used directly as the scalar bandwidth factor (scipy’s scalar bw_method); bw_scale is ignored.
• The kernel variance is factor² * var(fit), where var is the ddof=1 sample variance of the fit set.

Return value: the raw, un-normalized KDE density — exactly scipy’s model(x), i.e. Σ_i exp(-(x - fit_i)² / (2·kvar)) / (m·√(2π·kvar)). It is NOT rescaled to a max of 1.

Subsample determinism (deliberate scipy divergence): when max_samples is set and data is larger, the fit set is a seeded subsample of max_samples distinct points (a partial Fisher–Yates over the indices driven by a self-contained LCG keyed on seed). This is reproducible across runs/platforms but does not reproduce scipy’s unseeded np.random.choice draw — the goal is deterministic output, not byte-identical sample selection. Scott’s rule then uses the post-subsample size m = max_samples.

Kde::evaluate requires at least 2 data points (the ddof=1 variance is otherwise undefined) and raises a runtime error when fewer are supplied or when either tensor is null.

var data = Tensor {};
data.init(TensorType::f64, Array<int> { 5 });
var samples = Array<float> { -2.0, -1.0, 0.0, 1.0, 2.0 };
for (i: int, v: float in samples) { data.set(Array<int> { i }, v); }

var grid = Tensor {};
grid.init(TensorType::f64, Array<int> { 3 });
var gvals = Array<float> { -1.0, 0.0, 1.0 };
for (i: int, v: float in gvals) { grid.set(Array<int> { i }, v); }

// Scott's-rule bandwidth, no subsampling
var density = Kde::evaluate(data, grid, null, null, null, 0);

Climate Utilities

// Calculate Universal Thermal Climate Index
var utci_temp = utci(
    25.0,   // outdoor air temperature (°C)
    3.0,    // average wind speed (m/s)
    30.0,   // mean radiant temperature (°C)
    50.0,   // relative humidity (%)
);

ComputeEngine — Low-Level API

For advanced use cases, the ComputeEngine provides direct access to computational graphs.

// Define a compute model with custom operations
var model = ComputeModel {
    layers: Array<ComputeLayer> {
        ComputeLayerCustom {
            name: "ops",
            vars: Array<ComputeVariable> {
                ComputeVarInOut { name: "a", with_grad: false, shape: Array<int> { 3, 2 }, type: TensorType::f64 },
                ComputeVarInOut { name: "b", with_grad: false, shape: Array<int> { 3, 2 }, type: TensorType::f64 },
                ComputeVar { name: "c" },
            },
            ops: Array<ComputeOperation> {
                ComputeOperationAdd { input: "a", input2: "b", output: "c" },
            },
        },
    },
};

var engine = ComputeEngine {};
engine.compile(model, 10);
engine.configure(true);  // forward-only mode
engine.initialize();

// Set inputs and execute
engine.getVar("ops", "a")?.fill(1.5);
engine.getVar("ops", "b")?.fill(2.5);
engine.forward("ops");
var result = engine.getVar("ops", "c");  // 4.0

// State persistence
var state = ComputeState {};
engine.saveState(state);
// ... later ...
engine.loadState(state);

Available Operations

Arithmetic: Add, Sub, Mul, Div, Pow, MatMul, Scale, RaiseToPower, AddBias

Unary math: Abs, Neg, Sqrt, Exp, Log, Sign

Trigonometric: Sin, Cos, Tan, Asin, Acos, Atan, Sinh, Cosh, Tanh

Activation: Relu, LeakyRelu, Sigmoid, Softmax, Softplus, SoftSign, Selu, Elu, Celu, HardSigmoid, LogSoftmax, LeCunTanh

Reduction: Sum, Avg, ArgMin, ArgMax, SumIf, Euclidean

Utility: Fill, Filter, Clip