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 / origin	Accumulated count + per-feature [N] stats (nullable before the first learn). sum/sum_square are f64 running sums shifted by origin (the first learned row): Σ(x−origin), Σ(x−origin)² — variance/TSS formulas are shift-invariant, and the shift keeps them cancellation-free for large-mean features (raw timestamps, cumulative registers). Raw stats recover as Σx = total·origin + sum.
learn(input)	Learn from a [batch x N] (or 1-D [N]) tensor; O(batch·N), no [N,N] matrix. Accumulates in f64 whatever the input dtype. Rejects NaN/±Inf with an error, leaving the profile untouched.
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
rescaled_min_max()	Derive the profile of the min-max-scaled space in O(F) (no data pass) — what learn over min_max_scaling(data) would produce. Pass it to Kmeans::learning after scaling; keep the raw receiver for inverse scaling.

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 on the ComputeEngine, with meta-learning (best-of-N restarts) and an elbow-method search for K.

// Low-level engine loop (the high-level entry points below do all of this for you)
var engine = ComputeEngine {};
var model = Kmeans::configure(clusters, features, TensorType::f64, true);
engine.configure(true); // forward-only BEFORE compile — skips gradient zones (~2x less arena at large K)
engine.compile(model, batchSize);
Kmeans::initialize(engine, 42);
for (var round = 0; round < rounds; round++) {
    Kmeans::init_round(engine);
    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 });
}
Kmeans::calculate_stats(engine);
var centroids = Kmeans::getClustersCentroids(engine); // [K, F]
var assignment = Kmeans::cluster(engine, newData);    // inference

API Reference — Kmeans (engine level)

Method	Description
Kmeans::configure(clusters, features, type, stats)	Build the compute model (materialized [batch, K] distance matrix)
Kmeans::configureFused(clusters, features, type, stats, fast_math)	FUSED model — euclidean+argmin in one op, no distance matrix; what learning/score run; bit-identical results for f64
Kmeans::initialize(engine, seed)	Configure forward-only, set seed, initialize
Kmeans::init_round(engine) / end_round(engine)	Reset round accumulators / apply the centroid update
Kmeans::learn(engine, batch)	Train on one mini-batch
Kmeans::calculate_stats(engine)	Compute cluster statistics (once, after all rounds)
Kmeans::cluster(engine, batch)	Assign a batch, returning the [batch] assignment
Kmeans::getClustersCentroids(engine)	Centroids [K, F]
Kmeans::getClustersCounts(engine)	Member count per cluster [K]
Kmeans::getSumOfDistances(engine)	Total loss
Kmeans::getDistances(engine)	Last batch’s distances to ALL centroids [batch, K] — MATERIALIZED engines only (configure/getInferenceEngine); throws on fused engines (no distance variable)
Kmeans::getAssignment(engine) / getBestDistances(engine)	Last batch’s assignment [batch] / distance to its closest centroid [batch]
Kmeans::getClustersSumOfDistances(engine)	Per-cluster Σ distance [K]
Kmeans::getClustersAvgOfDistances(engine)	Per-cluster avg distance [K] (after calculate_stats)
Kmeans::getClustersDistancesToEachOther(engine)	Inter-centroid distances [K, K] (after calculate_stats)
Kmeans::sortClustersWith(engine, scratch)	Deterministic cluster order (ascending centroid coordinate sum); scratch = caller-owned [K, F] tensor (model dtype), overwritten
Kmeans::getInferenceEngine(result, maxBatchSize, stats)	Materialized inference engine from a trained result (keeps getDistances) — costs maxBatchSize·K·elem arena bytes plus a full matrix write+read per batch
Kmeans::getFusedInferenceEngine(result, maxBatchSize, stats)	FUSED inference engine — the predict-only fast path: same cluster/getAssignment/getBestDistances results (bit-identical for f64), but the [batch, K] matrix is never allocated (per-row arena F·elem + 16 instead of (K + F)·elem + 16 bytes); getDistances throws on it

Every engine get* accessor returns a view into the engine arena: it is overwritten by the next pass and keeps the whole arena alive while referenced — clone() before storing beyond the engine’s lifetime. All tensors on a returned KmeansResult are already standalone clones.

High-level training

The entry points build/compile the engine, seed the centroids, run the rounds, and return a KmeansResult. tensors is any: a single Tensor (one batch), an Array<Tensor>, a graph-backed nodeList<Tensor>, or a GcbReader<Tensor> streaming mini-batches from a .gcb file. Per-run knobs live in one KmeansConfig record — build with KmeansConfig::of(k, rounds, seed) and mutate the 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). Parallel note: the join holds EVERY restart’s result — incl. its two [N] assignment/distances tensors (8·N B each, 16·N B per restart) — until the call returns; budget metaRounds at large N
Kmeans::meta_meta_learning(tensors, config, minClusters, maxClusters, step, stopRatio, metaRounds, parallel)	Elbow-method K search, minClusters..maxClusters inclusive in steps of step; ignores config.nb_clusters (per-K clone(config))
Kmeans::estimate(config, totalRows, nbFeatures, tensorType, batchRows, throughput)	Size a run BEFORE building anything (KmeansEstimate) — see Sizing a run
Kmeans::calibrate(msBudget)	Measure this machine’s GFLOP/s + GB/s on the actual kernels (KmeansThroughput); graph-saved — calibrate(0) reuses the save (measuring only if none exists), a positive budget forces a re-measure

gaussian (also on Kmeans::score) is a GaussianArray? — the lean per-feature (diagonal) profiler: per-feature min/max/avg/std + sum/sumSq in O(N·F) (GaussianND’s [F,F] covariance is O(N·F²) and never read by k-means). Pass null and learning/single_learning/score compute it locally via Kmeans::featureProfile(tensors) (one extra O(N·F) pass); pass a precomputed one when calling repeatedly over the same data (e.g. train then score). meta_learning/meta_meta_learning build it once internally across restarts/the K sweep (no parameter).

KmeansConfig fields (only nb_clusters/nb_rounds/seed required; the rest nullable → 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 (decorrelated 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(8·K, min(40·K, 40960), ceil(√N)) capped to [1, N]	Rows for seed_strategy == subsample — the 8·K tail keeps S ≪ N at large K
kmeans_parallel_rounds	clamp(ceil(log2 N), 1, 3)	Oversample rounds for kMeansParallel (each is a FULL data pass; 2–3 saturate coverage)
kmeans_parallel_oversample	2.0	Oversample factor l = oversample · k for kMeansParallel. Must be > 0: learning/estimate throw on a non-positive/NaN value (it would silently produce zero D² candidates)
metrics	minimal	Quality-metric cost tier (KmeansMetrics, one ordered dial) — see the tier table
max_engine_memory	null	Byte budget for the engine arena: learning estimates the arena BEFORE compiling and throws with the largest fitting batch size instead of OOM-ing the worker
tol	0.0	Lloyd early-stop on the ABSOLUTE max centroid shift (feature-space units, data-scale dependent); never fires on a round that replaced empty clusters; <= 0 runs all nb_rounds
rel_tol	off	Scale-free early-stop: stop once relative inertia improvement (SSE_prev − SSE_cur)/SSE_cur stays below it for 2 consecutive non-replaced rounds (one O(N) native Σd² scan per round while on; null/<= 0 = off)
fast_math	false	Opt-in FAST-MATH kernels for every LEARNING-time distance pass: the fused assignment kernel (register-blocked plain-FMA squared-distance accumulation, argmin on squares, sqrt deferred to the winner — measured 5.1× per pass vs the compensated kernel, AVX2, K=4096/F=1024) and the kMeansParallel seeding folds + candidate reduction. Trade-off: last-ulp differences (argmin/min ties can flip), platform/build-dependent. score, all reported metrics (silhouette included) and the kmeansPlusPlus init route always run the exact kernel
update_rule	lloyd	lloyd (full-batch) or miniBatch (Sculley) — see Mini-batch training
mini_batch_rows	1024 clamped [1, N]	Rows per mini-batch step (miniBatch); ignored for GcbReader inputs (stored batches ARE the steps)
mini_batch_steps	nb_rounds · ceil(N / step_rows)	Upper bound on update steps (nb_rounds keeps an “epochs” meaning); step_rows = mini_batch_rows, or the MEAN stored batch size (N / batchCount) for a GcbReader (whose steps consume whole stored batches); EWA early stop usually finishes sooner
distance_metric	euclidean	euclidean or cosine (spherical) — see Cosine clustering
report	false	Emit live-progress files — see Live progress & cancellation
report_prefix	""	Filename prefix for the report files + cancel flag (concurrent-run isolation). A filename token, NOT a path (files live under File::userDir(); avoid / and ..). REQUIRES report: true: the entry points throw on a non-empty prefix with report off (every write AND cancel-flag poll is gated on report, so the combo would run blind and uncancellable)
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 observer (KmeansProgress) for the sequential paths — see In-process progress
final_repair	true	Repair clusters left empty by the FINAL assignment: up to 2 extra replace+reassign passes (strictly loss-non-increasing; stops when a pass improves neither empty count nor loss; off-data replace modes swap to farthestPoint). Each pass appends to roundsLoss/centroidShift, bumps nIter/emptyReplacements (so nIter <= nb_rounds + 2). Skipped when nb_rounds == 0 and after a cancel. false disables
compact_empty_clusters	false	Drop still-empty clusters from the result (all per-cluster tensors), remap assignment, update nbClusters (see Kmeans::compactEmptyClusters). Only matters when empties survive Lloyd + final repair (K > distinct points)

A null init_mode/replace_mode means unpinned: single-run entry points resolve to the fixed default above; meta_learning/meta_meta_learning randomize per restart — the init slot over the three cheap modes (fromInput/randomUniform/randomNormal; pin init_mode = kmeansPlusPlus to pay its K seeding passes on every restart), the replace slot over all six (Kmeans::nextRandomReplaceMode), so the strong on-data replacements participate. farthestPoint/splitLargest are replace-only: learning throws when either is passed as init_mode.

Mini-batch training (update_rule = miniBatch)

Sculley (2010) mini-batch k-means — the scale switch for very large N: each step samples mini_batch_rows rows, assigns them on the fused engine, and moves each TOUCHED centroid toward the batch mean with a per-center 1/cum_count learning rate. O(b·K·F) per step instead of O(N·K·F) per Lloyd round, for typically a few % worse inertia.

• nb_rounds becomes an EPOCH bound (steps ≤ nb_rounds · ceil(N/b)); the loop stops early once the EWA-smoothed batch inertia fails to improve for 10 consecutive steps. tol/rel_tol do not apply; centroidShift stays empty (a shift is a learning-rate artefact — final-repair passes included).
• roundsLoss holds the per-step EWA inertia curve (nIter = steps run); KmeansProgress.onRound fires per step, so KmeansProgressTracker.current_loss is the live smoothed inertia.
• Dead centers are never repaired per batch (at large K most clusters legitimately see zero rows per step): every max(10, 2·ceil(K/b)) steps, centers untouched over the whole window are re-seeded by a D²-weighted draw of DISTINCT rows from the current batch (capped at max(1, K/20) per window), counted in emptyReplacements. config.replace_mode plays no role during the steps — under miniBatch it applies only to the final_repair passes.
• The returned assignment/distances/loss and ALL metrics come from ONE exact full pass after training (plus the usual final_repair), so results stay comparable with Lloyd runs.
• A GcbReader<Tensor> cannot random-access rows: its stored batches are consumed sequentially as the mini-batches, cycling at EOF (mini_batch_rows ignored — the effective step size b, which scales the default step budget, the EWA constant, and the recycle window, is the reader’s MEAN stored batch size N / batchCount).
• Kmeans::estimate models the mode: flops_per_round/bytes_streamed_per_round describe one STEP, rounds_bound the step bound, total_flops adds the exact final passes.

Sizing a run before you start — Kmeans::estimate / Kmeans::calibrate

estimate(config, total_rows, nb_features, tensor_type, batch_rows, throughput) returns a KmeansEstimate: memory (engine arena at the resolved batch, dataset, result tensors, metric extras, peak_bytes — a slightly CONSERVATIVE sum: the arena is released before the metrics tail, so it never truly coexists with the metric extras — and max_batch_rows, the largest batch fitting both the 2^31−1 elements-per-variable cap and config.max_engine_memory), work (FLOPs + streamed bytes for seeding / one pass / the metrics tier, resolved_seed_strategy after fall-throughs, rounds_bound = the true per-restart bound — Lloyd rounds, or mini-batch update STEPS under update_rule == miniBatch — plus 2 when final_repair can extend the run), and — given a KmeansThroughput — predicted seconds per pass and total (roofline: max of compute-bound and memory-bound; expect ±30% on the calibrated machine). Pure arithmetic: nothing is allocated or compiled; the formulas are the same private helpers the training path uses, mirroring seedClusters’/learning’s exact resolution — including the FAILURES: it throws exactly where the run would (replace-only init modes, non-positive kmeans_parallel_oversample, a miniBatch step batch tripping the memory budget / element cap), never silently capping to a feasible-looking configuration learning rejects — and the fall-throughs to direct (degenerate subsample S >= N, stat-only randomUniform/randomNormal init). batch_rows == null resolves to min(max_batch_rows, total_rows), echoed in KmeansEstimate.batch_rows. Every kernel-generation change updates the traffic model + its pinned tests in the same change — estimate describes the CURRENT kernels.

calibrate(ms_budget) runs two ~(ms_budget/2) micro-benchmarks: a cache-resident tile of the exact kernel (GFLOP/s) and a 256 MB native copy (GB/s, read+write counted) — ~512 MB transient allocation, may overshoot the budget by up to one copy pass. The result is saved in the graph: calibrate(0) returns the save (measuring once with a 250 ms default only when nothing is saved) — the everyday call; a positive budget is the explicit re-measure (hardware change / kernel-generation upgrade).

var cfg = KmeansConfig::of(20_000, 10, 42);          // 500k × 2000 × K=20000 f64 target
var tp  = Kmeans::calibrate(0);
var est = Kmeans::estimate(cfg, 500_000, 2_000, TensorType::f64, null, tp);
// est.max_batch_rows           = 1_073_741 (fused engines bind on [batch, 2000] alone)
// est.arena_bytes              ~ 8.1 GB at batch 500_000 (batch·(F·8+16) + K·F·16)
// est.bytes_streamed_per_round ~ 2e13 under the cache-tiled kernel
// est.est_seconds_per_round    = max(compute time, memory time)  — roofline
// est.rounds_bound             = 12 (10 rounds + up to 2 final-repair passes)

KmeansMetrics — one ordered cost tier. Each tier is a superset of the cheaper ones. loss, inertia, and per-cluster inertia/max/std are always computed (they only reuse tensors the engine already produced):

Tier	Adds	Extra cost
minimal (default)	—	none
variance	explainedVariance, calinskiHarabasz	O(F) reads of the sum/sumSq profile; nearly free in learning (it builds the GaussianArray anyway) — only 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	model-dtype [K,K] variable (K²·elemSize B — 4 for f32, 8 for f64) + O(K²·F) pass + O(K²) DB loop + two [K,F] f64 tensors. Keep OFF at large K

Multi-restart meta runs (meta_learning / meta_meta_learning) pay the tiers ≥ silhouette ONCE, on the winning restart after selection — non-best restarts (including the result handed to KmeansProgress.onRestart) carry variance-tier metrics only; read silhouette/DB/feature stats from KmeansMetaResult.bestResult.

Choosing the input carrier (numerics are identical; only where the batches live differs). For parallel == true meta-learning, prefer a graph-backed nodeList<Tensor> (a raw Array<Tensor> is copied into every worker):

Carrier	Batches live	Peak resident memory	Use when
Array<Tensor>	all in RAM	the whole dataset	it comfortably fits in memory
nodeList<Tensor>	graph (gcdata/)	one batch (paged)	batches stored/reused, or parallel meta-learning
GcbReader<Tensor>	flat .gcb on disk	one batch + [K,F] centroids	too big for RAM and not worth persisting in the graph

Streaming from a .gcb file (GcbReader<Tensor>)

Cluster data bigger than RAM without materializing an Array<Tensor> or persisting a nodeList: write the mini-batches once with a GcbWriter<Tensor>, then hand any entry point a reader.

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();

var reader = GcbReader<Tensor> { path: path };
var result = Kmeans::learning(reader, Kmeans::featureProfile(reader), KmeansConfig::of(k, rounds, seed), 0);

The training paths iterate the reader forward one tensor at a time and rewind it (pos = 0) at the start of every pass — seeding, each Lloyd round, the stats passes — so the file is re-read once per pass (× restarts × the K sweep). Peak memory is your batch size, not N. Results are bit-identical to the same batches passed as Array/nodeList (pinned by kmeans_gcbreader_learning_matches_array). Two constraints (the reader is forward-only, read-only):

1. The in-place helpers (scaleBatches/sphericalizeBatches/setBatch) reject a reader — pre-scale / pre-L2-normalize before writing the .gcb.
2. For parallel == true meta-learning give each job its own reader (a shared pos races) — or use a nodeList<Tensor>.

Sizing at large K. The engine compiles to the LARGEST mini-batch. learning/score run a fused engine that never materializes [batch, K]: the arena is batch·(F·elemSize + 16) + K·F·(elemSize + 8) bytes (the K·F·8 term is the always-f64 SumIf accumulator) and the 2^31−1 elements-per-variable cap binds on [batch, F] alone. Only getInferenceEngine (materialized, for getDistances()) still pays elemSize·K per row, and there batch · K must stay below 2^31−1 — mini-batch the input rather than passing one giant Tensor; predict-only callers should use getFusedInferenceEngine instead (same assignments, no matrix). Set max_engine_memory to fail fast with the largest fitting batch size instead of OOM-ing.

Rounds & early stop. nb_rounds is the maximum Lloyd iteration count: each round performs one real centroid update + re-assignment and appends one roundsLoss entry, so roundsLoss.size() == result.nIter — equal to nb_rounds unless tol/rel_tol stops early (nIter < nb_rounds) or final_repair extends it (nIter <= nb_rounds + 2). Exception: a cooperative cancel discards the aborted pass’s partial loss, so a cancelled run may return nIter == roundsLoss.size() + 1. The returned loss/centroids/assignment are always mutually consistent. tol acts on the max per-row centroid shift ‖Δcentroid‖: the non-squared loss is NOT monotone under Lloyd (the mean update minimizes Σd², not Σd), so the shift — which settles to 0 as the means converge — is the stop signal, not a loss delta.

Search diagnostics & quality metrics. KmeansResult additionally exposes (all nullable, schema-safe):

• nbClusters (= centroids.shape()[0]), seedStrategy (strategy ACTUALLY used; kMeansParallel ignores initMode), nIter, emptyReplacements (passes that had to replace an empty; repair passes included).
• inertia + clustersInertia — Σ squared distance (SSE, the standard objective), distinct from loss (Σ non-squared distances).
• explainedVariance, calinskiHarabasz (from variance tier); daviesBouldin (full tier only — needs [K,K]).
• silhouette + clustersSilhouette (from silhouette tier): the simplified centroid-based form with a = distance to the nearest centroid, so it lies in [0, 1] (never negative like the classical pairwise form) — read it as a normalized separation margin, compare relatively across K. SINGLETON clusters contribute 0 (Rousseeuw/sklearn convention — one-point clusters cannot inflate the mean); EMPTY clusters are excluded from the nearest-other search (mirroring Davies–Bouldin’s occupied-only rule).
• clustersMaxDistance (radius [K]) and clustersDistanceStd [K] — always set. The std is the POPULATION std (ddof=0), derived from the raw per-cluster Σd/Σd² moments — subject to catastrophic cancellation when member distances are near-identical (σ/μ ≲ 1e-8): treat tiny values as “tight”, not exact.
• clustersFeatureVariance / clustersFeatureMae [K, F] (full tier; f64-typed for both f32 and f64 models — the kernels accumulate in f64).
• centroidShift [nIter] — per-iteration max ‖Δcentroid‖, measured AFTER empty-cluster replacement (expect spikes on emptyReplacements passes; tol never stops on those).

meta_meta_learning returns the full elbow trace on KmeansMetaResult: clusterCounts, clusterLosses (Σd curve), clusterInertias (the SSE curve the elbow decides on), and bestNbClusters. Both meta entry points also populate restarts — one KmeansRestartConvergence per restart of the kept K (index, nIter, final loss/inertia, full roundsLoss curve, emptyClusters = clusters still empty at convergence; with default final_repair this is > 0 only when even repair could not fill them). All modes (sequential, parallel, auto-K) build these through the same Kmeans::restartRecord collapse point. Kmeans::emptyClusterCount(counts) is the underlying helper (zero entries in a [K] counts tensor; null → 0).

Per-K browsing (collect_per_k). By default meta_meta_learning keeps only the elbow K’s clustering. With config.collect_per_k = true, KmeansMetaResult.perK holds one KmeansPerK per swept K, in sweep order, aligned 1:1 with clusterCounts (including non-accepted Ks and the stopping K): k, accepted (passed the elbow test), chosen (k == bestNbClusters, exactly one), and result — that K’s full KmeansResult (loss/inertia live there). No extra compute, but each retained result holds O(N) tensors, so it is off by default; perK is null when off (and for plain meta_learning).

Cross-K seed contract. Every swept K runs from the same base config.seed; restart i uses seed + i · 1_000_003 (wide prime stride, decorrelated independently of nb_rounds). A sequential run is deterministic for a given (seed, data) — pinned by kmeans_meta_meta_cross_k_seed_determinism.

Selection criterion. Restart selection and the elbow rank by inertia (SSE — what Lloyd minimizes), not loss (whose minimizer is the geometric median). The elbow compares each K against the best SSE seen so far and never stops on a stochastic uptick, so a noisy curve cannot produce a false elbow. Single-run numerics are unchanged; only multi-restart/elbow selection uses inertia.

Modes (KmeansMode). Initialization: fromInput (sample real rows), randomUniform, randomNormal, kmeansPlusPlus (D²-weighted seeding). Replacement-only: farthestPoint, kmeansPlusPlus-as-replace (D²-resample against the running distances), splitLargest (split the highest-inertia cluster, placing the new centroid on its farthest member — useful when k > distinct samples). The first three work in both contexts.

Scalable seeding (KmeansSeedStrategy) — orthogonal to init_mode: it picks over how much data the recipe runs, not what a centroid is (distinct from train-on-a-sample-then-score, which is a separate workflow):

Strategy	What it does	When
direct (default)	Run init_mode over all N rows.	Small/medium N.
subsample	Seed on subsample_size rows, then run the full Lloyd loop over all data. Only meaningful for the row-consuming init modes (fromInput/kmeansPlusPlus): the stat-only randomUniform/randomNormal route straight to direct (they read the full-data profile, not sampled rows — no gather pass, better bounds). A resolved S >= N degenerates the same way: seeding falls through to direct (equivalent, cheaper), logged. In both cases KmeansResult.seedStrategy reports direct.	Large N / large k — the recommended default (S ≈ 8k at large k).
kMeansParallel	k-means‖: kmeans_parallel_rounds full-data passes drawing ~l = oversample·k D²-weighted candidates each (capped at max(k, min(N, 4k)) total), then a weighted k-means++ reduction to k.	Streamed inputs (GcbReader) where random row access is expensive; on in-memory data subsample is usually the better deal.

Cosine / shape clustering (distance_metric = cosine). To cluster by direction, L2-normalize the rows first with Kmeans::sphericalizeBatches(tensors) (on training data and anything you later score), then set the metric. This runs spherical k-means: learning re-projects centroids onto the unit sphere after each update so the Euclidean argmin equals the cosine argmax. Reported distances/loss are chord distances (d = sqrt(2·(1 − cosθ))); inertia equals Σ 2·(1 − cos similarity); result.centroids are unit-norm.

Kmeans::sphericalizeBatches(batches);
var g = Kmeans::featureProfile(batches);
var cfg = KmeansConfig::of(k, rounds, seed);
cfg.distance_metric = KmeansDistance::cosine;
var r = Kmeans::learning(batches, g, cfg, 0);  // learning() takes every carrier; single_learning() is the one-Tensor-only variant

Representative member (medoid). Kmeans::nearestMemberPerCluster(result) → [K] i64 tensor of the global row index (in batch-consumption order) of the member closest to each centroid; -1 = empty cluster.

Scaling to very large N·F·K (e.g. 1M × 1000 × 10000): default seeding (kmeansPlusPlus + direct) costs O(K·N·F) before Lloyd — use seed_strategy = subsample (or kMeansParallel) and a positive tol. Keep metrics = minimal/variance at large K (silhouette is O(N·K); full’s [K,K] matrix is O(K²) work/memory, ~800 MB at K=10000). Prefer f32 to halve data/centroid RAM (metrics stay f64 — see Precision). All hot loops are native single-threaded kernels; parallelism comes from the worker pool at the meta-learning (restart) level. Beyond-RAM data: stream a GcbReader<Tensor>.

Reusable helpers

The library owns the generic orchestration; the project keeps feature extraction, storage, and metric/persistence decisions. None of these persist anything.

Helper	What it does
Kmeans::score(trained, tensors, gaussian, metrics)	Assign every row to trained’s centroids and (re)fill assignment/distances/loss, per-cluster count/sum/avg distance, and the metrics-tier quality metrics; returns the loss. Train on a sample, score the population — or re-score a loaded result. gaussian must describe the SCORED population; null computes it locally when the tier needs it. Throws on an EMPTY input (0 rows), like learning — it never wipes a trained result to a zero-loss, all-empty shell. Not cooperatively cancellable (no progress/report context)
Kmeans::featureProfile(tensors)	Learn a GaussianArray (per-feature min/max/sum/sumSq, O(N·F)) over all batches; streams a GcbReader once
Kmeans::totalRows(tensors)	Total rows across all batches (streams a reader forward)
Kmeans::batchCount(tensors) / Kmeans::batchAt(tensors, i)	Random-access carrier helpers: number of mini-batches / mini-batch i (a nodeList<Tensor> resolves only the requested batch). Both throw on a GcbReader<Tensor> (forward-only)
Kmeans::setBatch(tensors, i, t)	Write-back twin of batchAt (replace batch i); throws on a bare Tensor and on read-only carriers (GcbReader)
Kmeans::recomputeClusterAggregates(result, k)	Recompute per-cluster count/sum/avg distance from raw assignment + distances (used by score)
Kmeans::normalizeBatches(tensors)	Learn a GaussianArray, min-max-scale every batch IN PLACE, return Tuple<GaussianArray, GaussianArray>: .x = raw-space scaler (for inverse_min_max_scaling of trained centroids), .y = scaled-space profile (pass THIS to learning/meta_learning; the raw profile would seed/repair off-data and poison the variance metrics)
Kmeans::scaleBatches(tensors, g)	Min-max-scale every batch IN PLACE with an existing GaussianArray
Kmeans::sphericalizeBatches(tensors)	L2-normalize every row IN PLACE (prep for cosine)
Kmeans::nearestMemberPerCluster(result)	[K] i64 medoid indices (-1 = empty cluster)
Kmeans::compactEmptyClusters(result)	Remove still-empty clusters IN PLACE (per-cluster tensors, assignment remap, nbClusters); returns the count removed. Auto-called by config.compact_empty_clusters; works on loaded results
Kmeans::emptyClusterCount(counts)	Zero entries in a [K] counts tensor (null → 0)

Train-on-a-sample, score-the-population (the canonical large-N flow):

var trained = Kmeans::meta_learning(sampleBatches, cfg, metaRounds, parallel).bestResult!!;
var fullG = Kmeans::featureProfile(fullBatches);
Kmeans::score(trained, fullBatches, fullG, KmeansMetrics::full);

To normalize first: var p = Kmeans::normalizeBatches(batches) before training (pass p.y as the gaussian), then p.x.inverse_min_max_scaling(trained.centroids!!) after.

Precision — f32 vs f64. The compute dtype comes from the input tensors (or configure’s tensor_type); both are supported end to end. f32 halves the memory of data, centroids AND the [batch, K] matrix (materialized engines), and speeds up the distance loops. The centroid-sum accumulation runs in f64 even for f32 models (an f32 running sum absorbs nothing past 2^24), so the trade-off is confined to stored coordinates: result.centroids is f32 (rounded from the exact f64 mean); a materialized distance matrix is stored f32 but computed with f64 accumulation (argmin ties can differ in the last ulp). Everything measured stays f64:

• assignment is always i64; distances, min_distance, loss, inertia/clustersInertia are always f64 (upcast exactly, so loss/inertia and k-means++ sampling are dtype-stable).
• All derived metrics (silhouettes, DB, CH, explainedVariance, per-cluster max/std) are f64.
• clustersFeatureVariance/clustersFeatureMae are f64-typed for both dtypes.

If you normalize with GaussianArray, learn the profile in the same dtype as the data: min_max_scaling/standard_scaling output the profile’s dtype and reject a mismatch.

Live progress & cancellation

With report == true, the entry points publish JSON files under File::userDir() for an out-of-process poller (e.g. an @exposed endpoint). report == false (default) writes nothing. Files are published atomically (write .tmp, rename). <prefix> = config.report_prefix — set a unique per-run id (e.g. "run_<uuid>_") so concurrent runs don’t clobber each other’s singleton paths (empty prefix = the fixed default names; it is a filename token, not a path).

File	Type	Written by	Contents
<prefix>kmeans_restart_<i>.json	KmeansRestartFile	each restart of learning/single_learning	k, restart (= reportIdx), round, rounds (TRUE bound: nb_rounds + up to 2 repair passes), loss, done, partial roundsLoss curve (full once done)
<prefix>kmeans_progress.json	KmeansProgressControl	meta_meta_learning	running/done/cancelled, k_min/k_max/k_step/k_cur, best_k (= bestNbClusters once finished), best_loss, live k_values/k_losses elbow trace (rewritten after each K — no separate elbow file needed), total
<prefix>kmeans_cancel.flag	(presence only)	the caller	requests a cooperative stop at the next mini-batch boundary / restart / K. A pass cancelled mid-way is DISCARDED: its partial loss never enters roundsLoss; the returned result is best-effort (assignment/distances may mix the last two rounds) and carries KmeansResult.cancelled == true — restart/elbow selection never lets a cancelled (truncated, inertia-understating) result displace a completed one

KmeansProgressControl.total is the planned step count (nK · restarts · rounds, incl. repair margin) — an UPPER BOUND; tol, the elbow stop, and cancellation run fewer steps, so drive a progress bar off done/cancelled, not steps/total. To cancel, create the flag file; the run deletes it on exit and sets cancelled = true. The top-level entry points (single_learning, meta_learning, meta_meta_learning) purge a stale flag (and stale restart files) at start — learning() itself only polls the flag; clear it yourself when driving learning() directly with report. The file flag is polled per mini-batch inside a pass, per restart, and per K (only while report == true). Parallel meta_learning: each restart writes its own <prefix>kmeans_restart_<reportIdx>.json (parallel-safe); read exactly restarts_total of them.

In-process progress (KmeansProgress)

For a same-runtime consumer (Task, service), set config.progress to a KmeansProgress subtype — no files, no polling. GCL has no method override, so all five methods are abstract: implement each (no-ops where unneeded; 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 */ }
    }
    fn onRound(restart: int, round: int, loss: float) { /* fine progress bar */ }
    fn onBatch(restart: int, round: int, batchIndex: int, batchCount: int) { /* one mini-batch */ }
    fn onK(k: int, best: KmeansResult, accepted: bool) { /* one swept K finished */ }
    fn cancelled(): bool { return false; } // true => stop cooperatively
}
var cfg = KmeansConfig::of(8, 20, 42);
cfg.progress = MyProgress {};
var meta = Kmeans::meta_learning(batches, cfg, 5, false);

Callback contract: onRound fires per Lloyd iteration (per step in miniBatch mode) — sequential paths only; onBatch fires per mini-batch of every assignment pass, including the initial pre-loop pass (round == 0) — same restriction; onRestart fires once per finished restart in restart order from both sequential and parallel meta_learning (parallel restarts run in Job workers that cannot call back, so they get no live onRound/onBatch; onRestart fires after the join — and at tiers ≥ silhouette its result carries variance-tier metrics only, the expensive tiers being filled once on the winner); onK fires once per swept K. cancelled() is polled at every mini-batch of the sequential assignment passes and between rounds, restarts, and Ks. The observer is never serialized to parallel workers or persisted, and is fully orthogonal to the report files — use either, both, or neither.

Ready-made observer: KmeansProgressTracker — a concrete KmeansProgress aggregating what a progress UI needs: live position (current_restart/current_round/current_batch/current_k), current_loss (last completed round), the current restart’s live rounds_loss convergence curve (chartable), best loss/inertia across finished restarts (meta paths), elapsed(), an EMA-paced remaining() upper-bound against the planned rounds, progressRatio(), and cancel_requested (set true to cancel at the next mini-batch). Works for learning/single_learning, sequential meta_learning, and meta_meta_learning:

var cfg = KmeansConfig::of(8, 20, 42);
var tracker = KmeansProgressTracker::of(cfg, 5, 1); // 5 restarts, no K sweep
cfg.progress = tracker;
var meta = Kmeans::meta_learning(batches, cfg, 5, false);
// tracker.current_loss / best_loss / best_inertia; elapsed() remaining() progressRatio();
// tracker.cancel_requested = true;  // stop at the next mini-batch

For a time prediction BEFORE the first pace sample, pair it with Kmeans::estimate (see Sizing a run).

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.configure(true);  // forward-only mode — set BEFORE compile so gradient zones are never allocated
engine.compile(model, 10);
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, EuclideanArgMin (fused distance+argmin — no [N,K] matrix, optional fast_math; bit-identical to Euclidean+ArgMin for f64)

Reduction semantics worth knowing (doc-comments on the types in compute.gcl have the full contracts): Sum and SumIf accumulate into their outputs across forward() calls — pair them with a Fill to zero the accumulators each pass; SumIf silently drops rows whose class id falls outside [0, classes) and always emits f64 sums (f32 inputs widen; counts stays i64); Avg overwrites its output, writing 0.0 rows where counts == 0. Sum/SumIf/Avg accept only f32/f64 inputs, and Avg’s counts / SumIf’s ifCondition must be i64 [rows] vectors — all rejected at compile time.

Utility: Fill (in-place constant fill of a var; its value: any is converted numerically to the var’s dtype at compile — lossy float→int values and i32 overflows are compile errors), Filter, Clip