From 825ed6ad33b26f66e40a8cdfdbb50ba811421a56 Mon Sep 17 00:00:00 2001 From: w1n5t0n Date: Wed, 29 Apr 2026 15:55:43 +0300 Subject: [PATCH] feat(nisps/ml): MLP library with fixed-architecture template + spread-aware RL (meml-wmh) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Stream 2 of the clean-slate rewrite: nisps/ml/ replaces src/memlp/ with a header-only, heap-free MLP that satisfies nisps::core::MLEngine. Files (nisps/ml/): - activations.hpp — ReLU (leaky 0.01 for parity), sigmoid, tanh - loss.hpp — MSE per-sample (fixes meml-ues double-scaling: returns the sample's MSE without an extra 1/N multiplication; the training loop averages explicitly) - init.hpp — uniform/Xavier/spread-aware weight init - training.hpp — gradient clip helper (±10.0 matches legacy) - rl.hpp — move_weights with per-layer Xavier scaling, weight decay (10% * spread), gaussian noise via the deterministic Rng (matches the legacy JS sum-of-three-uniforms shape); draw_weights also spread-aware - stats.hpp — per-layer mean/max/dead/saturating diagnostics - mlp.hpp — 4-layer (3 hidden + sigmoid output) MLP class with std::array-backed weights, biases, gradient accumulators, dataset ring buffer (default 128 examples), loss history (default 4096 iters). Bias is a separate per-layer parameter — no input-vector mutation. Flat get_weights/set_weights layout: weights all layers (row-major, layer order), then biases all layers. Tests (tests/cpp/, all 50 passing under -Wall -Wextra -Werror -Wpedantic): - test_mlp_init.cpp — deterministic seeding, spread regimes, static_assert MLEngine concept satisfied - test_mlp_inference.cpp — golden hand-computed forward pass match, sigmoid output range, set_input bounds - test_mlp_training.cpp — XOR convergence (loss < 0.01 in <2k iters), ring-buffer eviction - test_mlp_loss.cpp — meml-ues regression test: reported loss equals hand-computed average MSE without extra 1/N scaling; sample weights honoured - test_mlp_rl.cpp — move_weights respects output_pin_mask (final-layer rows + biases preserved); spread regimes; grad clear after draw_weights - test_mlp_serialize.cpp — get_weights/set_weights round-trip preserves inference exactly; eval_loss is non-mutating; infer_batch matches individual inference Verification: - Clean build, no warnings - 50 tests pass (22 prior + 28 new) - No std::vector / new / malloc in nisps/ml/ - All float literals .f-suffixed in code (comments excepted) --- nisps/CMakeLists.txt | 6 + nisps/ml/activations.hpp | 91 +++++ nisps/ml/init.hpp | 72 ++++ nisps/ml/loss.hpp | 62 ++++ nisps/ml/mlp.hpp | 570 +++++++++++++++++++++++++++++++ nisps/ml/rl.hpp | 112 ++++++ nisps/ml/stats.hpp | 76 +++++ nisps/ml/training.hpp | 37 ++ tests/cpp/test_mlp_inference.cpp | 114 +++++++ tests/cpp/test_mlp_init.cpp | 98 ++++++ tests/cpp/test_mlp_loss.cpp | 125 +++++++ tests/cpp/test_mlp_rl.cpp | 133 ++++++++ tests/cpp/test_mlp_serialize.cpp | 146 ++++++++ tests/cpp/test_mlp_training.cpp | 92 +++++ 14 files changed, 1734 insertions(+) create mode 100644 nisps/ml/activations.hpp create mode 100644 nisps/ml/init.hpp create mode 100644 nisps/ml/loss.hpp create mode 100644 nisps/ml/mlp.hpp create mode 100644 nisps/ml/rl.hpp create mode 100644 nisps/ml/stats.hpp create mode 100644 nisps/ml/training.hpp create mode 100644 tests/cpp/test_mlp_inference.cpp create mode 100644 tests/cpp/test_mlp_init.cpp create mode 100644 tests/cpp/test_mlp_loss.cpp create mode 100644 tests/cpp/test_mlp_rl.cpp create mode 100644 tests/cpp/test_mlp_serialize.cpp create mode 100644 tests/cpp/test_mlp_training.cpp diff --git a/nisps/CMakeLists.txt b/nisps/CMakeLists.txt index 4253e4a..a40f531 100644 --- a/nisps/CMakeLists.txt +++ b/nisps/CMakeLists.txt @@ -51,6 +51,12 @@ if(NOT EMSCRIPTEN) ${NISPS_TEST_DIR}/test_ring_buffer.cpp ${NISPS_TEST_DIR}/test_rng.cpp ${NISPS_TEST_DIR}/test_math.cpp + ${NISPS_TEST_DIR}/test_mlp_init.cpp + ${NISPS_TEST_DIR}/test_mlp_inference.cpp + ${NISPS_TEST_DIR}/test_mlp_training.cpp + ${NISPS_TEST_DIR}/test_mlp_loss.cpp + ${NISPS_TEST_DIR}/test_mlp_rl.cpp + ${NISPS_TEST_DIR}/test_mlp_serialize.cpp ) target_link_libraries(nisps_core_tests PRIVATE nisps_core) diff --git a/nisps/ml/activations.hpp b/nisps/ml/activations.hpp new file mode 100644 index 0000000..3e8a5ee --- /dev/null +++ b/nisps/ml/activations.hpp @@ -0,0 +1,91 @@ +// nisps/ml/activations.hpp — activation functions used by the MLP. +// +// Three activations are needed for the firmware default: ReLU (hidden layers), +// sigmoid (output layer), and tanh (alternative). Each entry exposes the +// activation and its derivative-given-pre-activation; we keep both around +// because the backprop path multiplies by the derivative of the +// pre-activation, which is cheaper to compute on the inner product than on +// the post-activation in some cases. +// +// Notes on the leak in ReLU: +// The legacy firmware uses leaky ReLU with slope 0.01 everywhere it claims +// to use "ReLU". We preserve that behavior for parity — a true zero-slope +// ReLU would silently change training dynamics in ways that aren't +// documented. The leak slope is a `static const float` so it lives in SRAM +// per Chris's rule. +// +// Sigmoid: +// We use the exact std::exp version here (NOT fast_sigmoid). The MLP +// training loop multiplies by the sigmoid derivative, and the ~1.2% bias +// in fast_sigmoid would compound during training. fast_sigmoid is fine for +// inference-only paths but not for the gradient path. The ML output +// activation is the place where we want monotonic-AND-accurate. +// +// All literals carry the .f suffix; constants used inside loops live in +// `static const float` so they are hoisted to SRAM on RP2350. + +#pragma once + +#include + +#include "../core/perf.hpp" + +namespace nisps::ml { + +// Match the legacy firmware leaky-ReLU slope. 0.01 matches PyTorch default +// and what `src/memlp/Utils.h::kReLUSlope` had before. +inline constexpr float kReluLeakSlope = 0.01f; + +NISPS_FORCE_INLINE float relu(float x) noexcept { + return (x > 0.f) ? x : kReluLeakSlope * x; +} +NISPS_FORCE_INLINE float relu_deriv_pre(float pre_activation) noexcept { + return (pre_activation > 0.f) ? 1.f : kReluLeakSlope; +} + +NISPS_FORCE_INLINE float sigmoid(float x) noexcept { + // Saturate inputs to avoid expf overflow / underflow noise. Mirrors + // nisps::exact_sigmoid clamps from core/math.hpp. + if (x > 40.f) return 1.f; + if (x < -40.f) return 0.f; + return 1.f / (1.f + std::exp(-x)); +} +NISPS_FORCE_INLINE float sigmoid_deriv_pre(float pre_activation) noexcept { + const float s = sigmoid(pre_activation); + return s * (1.f - s); +} + +NISPS_FORCE_INLINE float tanh_act(float x) noexcept { + return std::tanh(x); +} +NISPS_FORCE_INLINE float tanh_deriv_pre(float pre_activation) noexcept { + const float t = std::tanh(pre_activation); + return 1.f - t * t; +} + +// Activation kind, dispatched at compile time per layer (see Layer template +// in mlp.hpp). We don't use a runtime tag because activation choice is +// architectural, not per-call. +enum class Activation : int { + ReLU = 0, + Sigmoid = 1, + Tanh = 2, +}; + +template +NISPS_FORCE_INLINE float activate(float x) noexcept { + if constexpr (A == Activation::ReLU) return relu(x); + if constexpr (A == Activation::Sigmoid) return sigmoid(x); + if constexpr (A == Activation::Tanh) return tanh_act(x); + return x; // unreachable +} + +template +NISPS_FORCE_INLINE float activate_deriv_pre(float pre) noexcept { + if constexpr (A == Activation::ReLU) return relu_deriv_pre(pre); + if constexpr (A == Activation::Sigmoid) return sigmoid_deriv_pre(pre); + if constexpr (A == Activation::Tanh) return tanh_deriv_pre(pre); + return 1.f; // unreachable +} + +} // namespace nisps::ml diff --git a/nisps/ml/init.hpp b/nisps/ml/init.hpp new file mode 100644 index 0000000..18bd9a5 --- /dev/null +++ b/nisps/ml/init.hpp @@ -0,0 +1,72 @@ +// nisps/ml/init.hpp — weight initialization strategies. +// +// Three strategies are exposed: +// - `uniform_init`: w ~ U[-1, 1]. Matches the legacy +// `utils::gen_rand()` baseline; produces high-magnitude pre-activations +// that drive sigmoid output saturation. +// - `xavier_init`: w ~ U[-1, 1] / sqrt(fan_in). Centered pre-activations, +// better for sigmoid output layers. +// - `spread_init`: linearly interpolates per-layer scale between the two. +// scale = (1 - spread) + spread / sqrt(fan_in). spread=0 ⇒ uniform, +// spread=1 ⇒ Xavier. This is the playground knob. +// +// All three operate on a flat row-major weight buffer of size fan_in*fan_out. +// Biases are initialized separately and ALWAYS to zero — the legacy code +// initialized biases to zero, and the playground spread parameter never +// touches biases. We keep that. + +#pragma once + +#include +#include +#include + +#include "../core/rng.hpp" + +namespace nisps::ml { + +// Compute the spread-aware weight scale for one layer. +// spread=0 → 1.0 (uniform [-1, 1]) +// spread=1 → 1/sqrt(fan_in) (Xavier) +inline float spread_scale(float spread, std::size_t fan_in) noexcept { + if (fan_in == 0u) return 1.f; + const float inv_sqrt = 1.f / std::sqrt(static_cast(fan_in)); + return (1.f - spread) + spread * inv_sqrt; +} + +// Initialize one layer's weights (flat row-major) and biases. +// weights: span of size fan_in * fan_out +// biases: span of size fan_out (zeroed) +// spread: see spread_scale +// rng: state advanced; caller owns it +inline void spread_init(std::span weights, + std::span biases, + std::size_t fan_in, + float spread, + Rng& rng) noexcept { + const float scale = spread_scale(spread, fan_in); + for (std::size_t i = 0; i < weights.size(); ++i) { + weights[i] = rng.next_float_signed() * scale; + } + for (std::size_t i = 0; i < biases.size(); ++i) { + biases[i] = 0.f; + } +} + +// Convenience aliases for the endpoints. Cheap to build on top of +// spread_init; useful for tests that want to assert a specific regime. +inline void uniform_init(std::span weights, + std::span biases, + std::size_t fan_in, + Rng& rng) noexcept { + spread_init(weights, biases, fan_in, 0.f, rng); +} + +inline void xavier_init(std::span weights, + std::span biases, + std::size_t fan_in, + Rng& rng) noexcept { + spread_init(weights, biases, fan_in, 1.f, rng); +} + +} // namespace nisps::ml diff --git a/nisps/ml/loss.hpp b/nisps/ml/loss.hpp new file mode 100644 index 0000000..6008f48 --- /dev/null +++ b/nisps/ml/loss.hpp @@ -0,0 +1,62 @@ +// nisps/ml/loss.hpp — MSE loss with a deliberate fix for the meml-ues +// double-scaling bug. +// +// THE BUG (in `src/memlp/Loss.h::MSE` and `src/memlp/MLP.cpp::Train`): +// The legacy MSE computed loss as +// sum((diff^2) / NOut) * sampleSizeReciprocal in Loss.h +// and then `Train()` did +// current_iteration_cost_function *= sampleSizeReciprocal; // ← second +// So the unweighted average loss was scaled by 1/N twice — once in MSE() +// and once in the outer training loop. Loss values for the same training +// set get reported 1/N too small, which is misleading for diagnostics. +// +// THE FIX: +// Loss is computed PER-SAMPLE here. The training loop is responsible for +// averaging across samples once. There is no implicit per-sample weight in +// the MSE function — that is the caller's policy decision. +// +// For a single sample: mse = (1/NOut) * sum_j (label_j - pred_j)^2 +// Derivative wrt pred_j: -(2/NOut) * (label_j - pred_j) +// +// When sample weights are supplied to `train()` they are applied at the +// sample level (loss_total = sum_i weight_i * mse_i), and the gradient +// per-sample is scaled by weight_i. Sample weights must sum to 1.0 — the +// caller normalizes (matches legacy contract). + +#pragma once + +#include +#include + +#include "../core/perf.hpp" + +namespace nisps::ml { + +// Compute MSE loss for one sample and write the per-output derivative into +// `loss_deriv`. Returns the scalar loss for this sample. +// +// loss = (1/NOut) * sum_j (label_j - pred_j)^2 +// d(loss)/d(pred_j) = -(2/NOut) * (label_j - pred_j) +// +// `loss_deriv` and `pred` must both have size NOut. `label` likewise. +// +// NB: We deliberately do NOT multiply by any sampleSizeReciprocal here. The +// caller is responsible for averaging across samples (and scaling by sample +// weights, if any) — see training.hpp. +NISPS_FORCE_INLINE float mse_per_sample(std::span label, + std::span pred, + std::span loss_deriv) noexcept { + const std::size_t n = pred.size(); + if (n == 0u) return 0.f; + const float inv_n = 1.f / static_cast(n); + + float accum = 0.f; + for (std::size_t j = 0; j < n; ++j) { + const float diff = label[j] - pred[j]; // sign chosen so deriv is wrt pred + accum += diff * diff * inv_n; + loss_deriv[j] = -2.f * inv_n * diff; + } + return accum; +} + +} // namespace nisps::ml diff --git a/nisps/ml/mlp.hpp b/nisps/ml/mlp.hpp new file mode 100644 index 0000000..ce3e56f --- /dev/null +++ b/nisps/ml/mlp.hpp @@ -0,0 +1,570 @@ +// nisps/ml/mlp.hpp — fixed-architecture MLP, four layers (three hidden + +// output). All buffers are template-sized std::array; zero heap allocation +// in inference, training, and the dataset path. +// +// ARCHITECTURE +// MLP +// ┌──────┐ Linear+Bias ┌────────┐ ReLU ┌────────┐ ReLU ┌────────┐ Sigmoid +// │ NIn │ ─────────────▶ │ NH1 │ ──────▶ │ NH2 │ ──────▶ │ NH3 │ ──────▶ NOut +// └──────┘ └────────┘ └────────┘ └────────┘ +// Layer 0 (NIn → NH1) ReLU +// Layer 1 (NH1 → NH2) ReLU +// Layer 2 (NH2 → NH3) ReLU +// Layer 3 (NH3 → NOut) Sigmoid +// +// We support exactly three hidden layers. The legacy firmware default is +// [10, 10, 14], so the MVP signature directly matches `MLP`. Variable layer count is deferred — see architecture.md. +// +// MEMORY MODEL +// Per layer L_k with fan_in = N_in[k], fan_out = N_out[k]: +// std::array weights // row-major +// std::array biases +// std::array pre_activation // cached for backprop +// std::array activation // cached for backprop +// std::array grad_w_accum // for backprop +// std::array grad_b_accum +// +// Per MLP: +// std::array input_buffer (current set_input values) +// std::array output (post-final-activation; outputs()) +// std::array dataset_features +// std::array dataset_labels +// std::size_t dataset_count, dataset_head (FIFO ring buffer) +// std::array loss_history (max iters from train()) +// Rng rng_ +// std::array bp_err_buf, bp_delta_buf (backprop scratch) +// +// FLAT WEIGHT LAYOUT (`get_weights` / `set_weights`) +// [layer0_weights ...] [layer1_weights ...] [layer2_weights ...] [layer3_weights ...] +// [layer0_biases ...] [layer1_biases ...] [layer2_biases ...] [layer3_biases ...] +// Documented in detail near `weight_count()`. +// +// CONCEPT SATISFACTION +// The class satisfies `nisps::MLEngine`: +// set_input, process, outputs, add_example, train (no-arg overload +// returning float), move_weights(speed, spread), draw_weights(spread), +// reset, seed. +// Plus diagnostics required by the broader API (see Stream 2 brief in +// architecture.md): eval_loss, layer_stats, get/set_weights, weight_count, +// infer_batch, loss_history. + +#pragma once + +#include +#include +#include +#include +#include + +#include "../core/concepts.hpp" +#include "../core/perf.hpp" +#include "../core/rng.hpp" +#include "activations.hpp" +#include "init.hpp" +#include "loss.hpp" +#include "rl.hpp" +#include "stats.hpp" +#include "training.hpp" + +namespace nisps::ml { + +// Layer. Stores its weights, biases, and the work +// buffers needed for forward + backprop. Header-only, all sizes compile- +// time. Each method is small; the compiler will inline through. +template +struct Layer { + static constexpr std::size_t kFanIn = FanIn; + static constexpr std::size_t kFanOut = FanOut; + static constexpr Activation kAct = Act; + + std::array weights{}; + std::array biases{}; + // Cached during forward(); consumed during backprop(). + std::array pre_activation{}; + std::array activation{}; + // Gradient accumulators — used per-sample for SGD weight update. + std::array grad_w{}; + std::array grad_b{}; + + NISPS_FORCE_INLINE float& w(std::size_t node, std::size_t in) noexcept { + return weights[node * FanIn + in]; + } + NISPS_FORCE_INLINE float w(std::size_t node, std::size_t in) const noexcept { + return weights[node * FanIn + in]; + } + + // Forward: compute pre_activation and activation given an input span. + NISPS_HOT NISPS_FORCE_INLINE + void forward(std::span input) noexcept { + for (std::size_t node = 0; node < FanOut; ++node) { + const std::size_t row = node * FanIn; + float sum = biases[node]; + for (std::size_t j = 0; j < FanIn; ++j) { + sum += weights[row + j] * input[j]; + } + pre_activation[node] = sum; + activation[node] = activate(sum); + } + } + + // Compute incoming-error vector for the previous layer: + // delta_in[j] = sum_node (err_signal[node] * w[node, j]) + // Where err_signal[node] = upstream_err[node] * d/dpre activation. + // Also accumulates per-weight and per-bias gradients (no LR yet). + NISPS_HOT NISPS_FORCE_INLINE + void backprop_accumulate(std::span input, + std::span upstream_err, + std::span delta_in, + float sample_weight) noexcept { + for (std::size_t j = 0; j < FanIn; ++j) delta_in[j] = 0.f; + + for (std::size_t node = 0; node < FanOut; ++node) { + const float err_signal = + upstream_err[node] * activate_deriv_pre(pre_activation[node]) * sample_weight; + const std::size_t row = node * FanIn; + for (std::size_t j = 0; j < FanIn; ++j) { + grad_w[row + j] += err_signal * input[j]; + delta_in[j] += err_signal * weights[row + j]; + } + grad_b[node] += err_signal; + } + } + + // Apply accumulated gradient to weights+biases with clipping. Resets + // the accumulators to zero for the next sample/iteration. + NISPS_FORCE_INLINE + void apply_grad(float lr) noexcept { + for (std::size_t i = 0; i < FanIn * FanOut; ++i) { + const float g = clip_gradient(grad_w[i]); + weights[i] -= lr * g; + grad_w[i] = 0.f; + } + for (std::size_t i = 0; i < FanOut; ++i) { + const float g = clip_gradient(grad_b[i]); + biases[i] -= lr * g; + grad_b[i] = 0.f; + } + } + + NISPS_FORCE_INLINE + void clear_grad() noexcept { + for (std::size_t i = 0; i < FanIn * FanOut; ++i) grad_w[i] = 0.f; + for (std::size_t i = 0; i < FanOut; ++i) grad_b[i] = 0.f; + } +}; + +// MLP +// +// MaxIterTrain caps the loss-history buffer; if a caller asks for more +// iterations they will be honored at runtime, but only the first +// kMaxIterTrain are recorded for inspection. 4096 fits the playground's +// upper bound and costs 16 KiB. +template +class MLP { + public: + static constexpr std::size_t kInput = NIn; + static constexpr std::size_t kHidden1 = NHidden1; + static constexpr std::size_t kHidden2 = NHidden2; + static constexpr std::size_t kHidden3 = NHidden3; + static constexpr std::size_t kOutput = NOut; + static constexpr std::size_t kMaxExamples = NMaxExamples; + static constexpr std::size_t kMaxIterTrain = NMaxIterTrain; + static constexpr std::size_t kNumLayers = 4u; + + using Layer0 = Layer; + using Layer1 = Layer; + using Layer2 = Layer; + using Layer3 = Layer; + + // --------------------------------------------------------------- + // Lifecycle + // --------------------------------------------------------------- + explicit MLP(std::uint64_t seed) noexcept : rng_(seed) { + // Default-init weights with spread=1 (Xavier-like). The IML + // interface caller is expected to draw_weights() with the + // playground spread before the first inference; this default + // simply gives us a non-degenerate starting state for tests + // that skip an explicit draw. + draw_weights(1.f); + clear_dataset_(); + loss_history_count_ = 0u; + } + + // --------------------------------------------------------------- + // Inference API (concept: set_input / process / outputs) + // --------------------------------------------------------------- + NISPS_FORCE_INLINE void set_input(std::size_t i, float v) noexcept { + if (i < NIn) input_[i] = v; + } + + NISPS_HOT void process() noexcept { + forward_(std::span(input_)); + // Mirror final activation into the output buffer so callers can + // read a stable span. + const auto& a = layer3_.activation; + for (std::size_t i = 0; i < NOut; ++i) output_[i] = a[i]; + } + + NISPS_FORCE_INLINE std::span outputs() const noexcept { + return std::span(output_.data(), NOut); + } + + // --------------------------------------------------------------- + // Dataset / Training (concept: add_example / train) + // --------------------------------------------------------------- + // FIFO ring buffer; oldest example evicted when full. No allocation. + // We don't track logical insertion order during training because SGD + // doesn't care — the iteration order over the buffer is arbitrary. + void add_example(std::span features, + std::span labels) noexcept { + if (features.size() < NIn || labels.size() < NOut) return; + + std::size_t slot; + if (dataset_count_ < NMaxExamples) { + slot = dataset_count_++; + } else { + // Buffer full: overwrite the slot pointed at by head_ (oldest) + // and advance head_ to the next-oldest. + slot = dataset_head_; + dataset_head_ = (dataset_head_ + 1u) % NMaxExamples; + } + const std::size_t f_off = slot * NIn; + const std::size_t l_off = slot * NOut; + for (std::size_t i = 0; i < NIn; ++i) ds_features_[f_off + i] = features[i]; + for (std::size_t i = 0; i < NOut; ++i) ds_labels_ [l_off + i] = labels[i]; + } + + // Concept-required no-arg overload. Default learning rate matches the + // playground's "sane RL training" knob; max_iter and min_err follow. + float train() noexcept { + return train(1.f, 1000u, 0.001f, std::span{}); + } + + // Full SGD training. `sample_weights`, if non-empty, must size to the + // current example count and sum to 1.0 (caller's responsibility — we + // do NOT renormalize). + // + // Returns final epoch loss. Records per-iteration loss in + // `loss_history_` (bounded by kMaxIterTrain). + float train(float lr, + std::size_t max_iter, + float min_err, + std::span sample_weights = {}) noexcept { + loss_history_count_ = 0u; + if (dataset_count_ == 0u) return 0.f; + + const bool weighted = !sample_weights.empty(); + const float uniform_w = 1.f / static_cast(dataset_count_); + + float epoch_loss = 0.f; + for (std::size_t iter = 0; iter < max_iter; ++iter) { + epoch_loss = 0.f; + // SGD: per-sample forward → loss → backprop+update. The order + // is the dataset insertion order; we do not shuffle (matches + // the legacy `Train()` exactly — `TrainBatch` shuffles, but + // we're not implementing batch yet). + for (std::size_t s = 0; s < dataset_count_; ++s) { + const float w = weighted ? sample_weights[s] : uniform_w; + + // Forward pass on sample s. + std::span x = sample_features_(s); + forward_(x); + + // Per-sample loss (NOT scaled by 1/N — the meml-ues fix). + std::array deriv{}; + const float sample_loss = mse_per_sample( + sample_labels_(s), + std::span(layer3_.activation.data(), NOut), + std::span(deriv.data(), NOut)); + + // Aggregate weighted loss. + epoch_loss += w * sample_loss; + + // Backprop with the same w as the gradient scaler. + backprop_(x, std::span(deriv), w); + + // Apply gradient (per-sample, SGD). + layer3_.apply_grad(lr); + layer2_.apply_grad(lr); + layer1_.apply_grad(lr); + layer0_.apply_grad(lr); + } + + if (loss_history_count_ < NMaxIterTrain) { + loss_history_[loss_history_count_++] = epoch_loss; + } + + if (epoch_loss < min_err) break; + } + return epoch_loss; + } + + // --------------------------------------------------------------- + // RL ops (concept: move_weights / draw_weights) + // --------------------------------------------------------------- + void move_weights(float speed, float spread, + std::span output_pin_mask = {}) noexcept { + move_weights_layer(std::span(layer0_.weights), std::span(layer0_.biases), + Layer0::kFanIn, speed, spread, /*final=*/false, {}, rng_); + move_weights_layer(std::span(layer1_.weights), std::span(layer1_.biases), + Layer1::kFanIn, speed, spread, /*final=*/false, {}, rng_); + move_weights_layer(std::span(layer2_.weights), std::span(layer2_.biases), + Layer2::kFanIn, speed, spread, /*final=*/false, {}, rng_); + move_weights_layer(std::span(layer3_.weights), std::span(layer3_.biases), + Layer3::kFanIn, speed, spread, /*final=*/true, output_pin_mask, rng_); + } + + void draw_weights(float spread) noexcept { + draw_weights_layer(std::span(layer0_.weights), std::span(layer0_.biases), + Layer0::kFanIn, spread, rng_); + draw_weights_layer(std::span(layer1_.weights), std::span(layer1_.biases), + Layer1::kFanIn, spread, rng_); + draw_weights_layer(std::span(layer2_.weights), std::span(layer2_.biases), + Layer2::kFanIn, spread, rng_); + draw_weights_layer(std::span(layer3_.weights), std::span(layer3_.biases), + Layer3::kFanIn, spread, rng_); + layer0_.clear_grad(); + layer1_.clear_grad(); + layer2_.clear_grad(); + layer3_.clear_grad(); + } + + // Concept reset: clear weights, dataset, and loss history. Seed is + // intentionally NOT reset (use `seed()` for that). + void reset() noexcept { + clear_dataset_(); + loss_history_count_ = 0u; + // Re-init weights from current rng state with default spread. + draw_weights(1.f); + for (std::size_t i = 0; i < NIn; ++i) input_[i] = 0.f; + for (std::size_t i = 0; i < NOut; ++i) output_[i] = 0.f; + } + + void seed(std::uint64_t s) noexcept { rng_.seed(s); } + + // --------------------------------------------------------------- + // Diagnostics + // --------------------------------------------------------------- + // Forward pass + MSE on a single (input, label) implied by current + // input_ and the most recent training labels — i.e. "what would the + // loss be on the current input if the label were the current output?" + // For now we report the average loss across the training set without + // updating weights. Useful for non-destructive evaluation. + float eval_loss() const noexcept { + if (dataset_count_ == 0u) return 0.f; + const float inv_n = 1.f / static_cast(dataset_count_); + // We need a non-const forward pass to use the cached buffers; since + // eval_loss is logically const, fork a local computation that + // doesn't touch member buffers. That means recomputing through the + // layer weights against scratch arrays — no allocation, just stack. + float total = 0.f; + for (std::size_t s = 0; s < dataset_count_; ++s) { + std::array a1{}; + std::array a2{}; + std::array a3{}; + std::array ao{}; + + const std::size_t f_off = s * NIn; + forward_const_layer( + std::span(ds_features_.data() + f_off, NIn), + layer0_.weights, layer0_.biases, a1); + forward_const_layer( + std::span(a1), layer1_.weights, layer1_.biases, a2); + forward_const_layer( + std::span(a2), layer2_.weights, layer2_.biases, a3); + forward_const_layer( + std::span(a3), layer3_.weights, layer3_.biases, ao); + + float sse = 0.f; + const float inv_o = 1.f / static_cast(NOut); + const std::size_t l_off = s * NOut; + for (std::size_t j = 0; j < NOut; ++j) { + const float d = ds_labels_[l_off + j] - ao[j]; + sse += d * d * inv_o; + } + total += sse * inv_n; + } + return total; + } + + LayerStats layer_stats(std::size_t layer_idx) const noexcept { + switch (layer_idx) { + case 0: return compute_layer_stats(layer0_.weights, layer0_.biases); + case 1: return compute_layer_stats(layer1_.weights, layer1_.biases); + case 2: return compute_layer_stats(layer2_.weights, layer2_.biases); + case 3: return compute_layer_stats(layer3_.weights, layer3_.biases); + default: return {}; + } + } + + // Flat layout: layer0 weights, layer1 weights, layer2 weights, layer3 + // weights, then layer0..3 biases. The split is `weights first all + // layers, then biases all layers` so callers serializing weights can + // pre-compute offsets without consulting layer-specific tables. + static constexpr std::size_t weight_count() noexcept { + return NIn * NHidden1 + NHidden1 * NHidden2 + NHidden2 * NHidden3 + NHidden3 * NOut + + NHidden1 + NHidden2 + NHidden3 + NOut; + } + + // Returns a span into a member-owned scratch buffer that holds a copy + // of the flat weights+biases. The buffer is regenerated on each call, + // so don't hold onto the span across mutations. + std::span get_weights() noexcept { + std::size_t k = 0u; + // Weights, layer-major. + for (float v : layer0_.weights) flat_weight_buf_[k++] = v; + for (float v : layer1_.weights) flat_weight_buf_[k++] = v; + for (float v : layer2_.weights) flat_weight_buf_[k++] = v; + for (float v : layer3_.weights) flat_weight_buf_[k++] = v; + // Biases. + for (float v : layer0_.biases) flat_weight_buf_[k++] = v; + for (float v : layer1_.biases) flat_weight_buf_[k++] = v; + for (float v : layer2_.biases) flat_weight_buf_[k++] = v; + for (float v : layer3_.biases) flat_weight_buf_[k++] = v; + return std::span(flat_weight_buf_.data(), weight_count()); + } + + void set_weights(std::span w) noexcept { + if (w.size() < weight_count()) return; + std::size_t k = 0u; + for (float& v : layer0_.weights) v = w[k++]; + for (float& v : layer1_.weights) v = w[k++]; + for (float& v : layer2_.weights) v = w[k++]; + for (float& v : layer3_.weights) v = w[k++]; + for (float& v : layer0_.biases) v = w[k++]; + for (float& v : layer1_.biases) v = w[k++]; + for (float& v : layer2_.biases) v = w[k++]; + for (float& v : layer3_.biases) v = w[k++]; + } + + // Run inference on N points (each NIn-sized) and write N output vectors + // (each NOut-sized) into `outs`. NO heap. Modifies the internal cached + // activations as a side effect. + void infer_batch(std::span points, + std::span outs) noexcept { + const std::size_t n = points.size() / NIn; + if (outs.size() < n * NOut) return; + for (std::size_t i = 0; i < n; ++i) { + const std::size_t in_off = i * NIn; + for (std::size_t j = 0; j < NIn; ++j) input_[j] = points[in_off + j]; + process(); + const std::size_t out_off = i * NOut; + for (std::size_t j = 0; j < NOut; ++j) outs[out_off + j] = output_[j]; + } + } + + std::span loss_history() const noexcept { + return std::span(loss_history_.data(), loss_history_count_); + } + + std::size_t example_count() const noexcept { return dataset_count_; } + + void clear_examples() noexcept { clear_dataset_(); } + + private: + // --------------------------------------------------------------- + // Internal helpers + // --------------------------------------------------------------- + NISPS_HOT NISPS_FORCE_INLINE + void forward_(std::span in) noexcept { + layer0_.forward(in); + layer1_.forward(std::span(layer0_.activation)); + layer2_.forward(std::span(layer1_.activation)); + layer3_.forward(std::span(layer2_.activation)); + } + + // Backprop with sample_weight applied to every error signal (so the + // accumulated gradient is already weighted). No weight update happens + // here — caller does it after each sample. + NISPS_HOT NISPS_FORCE_INLINE + void backprop_(std::span input, + std::span output_deriv, + float sample_weight) noexcept { + std::array d3{}; + std::array d2{}; + std::array d1{}; + std::array d0{}; + + layer3_.backprop_accumulate( + std::span(layer2_.activation), + output_deriv, + std::span(d3), + sample_weight); + layer2_.backprop_accumulate( + std::span(layer1_.activation), + std::span(d3), + std::span(d2), + 1.f); // weight already in d3 + layer1_.backprop_accumulate( + std::span(layer0_.activation), + std::span(d2), + std::span(d1), + 1.f); + layer0_.backprop_accumulate( + input, + std::span(d1), + std::span(d0), + 1.f); + } + + // const forward pass for diagnostics. Doesn't touch layer caches. + template + static NISPS_FORCE_INLINE void forward_const_layer( + std::span in, + const std::array& w, + const std::array& b, + std::array& out) noexcept { + for (std::size_t node = 0; node < Fo; ++node) { + const std::size_t row = node * Fi; + float sum = b[node]; + for (std::size_t j = 0; j < Fi; ++j) sum += w[row + j] * in[j]; + out[node] = activate(sum); + } + } + + NISPS_FORCE_INLINE std::span sample_features_(std::size_t s) const noexcept { + return std::span(ds_features_.data() + s * NIn, NIn); + } + NISPS_FORCE_INLINE std::span sample_labels_(std::size_t s) const noexcept { + return std::span(ds_labels_.data() + s * NOut, NOut); + } + + void clear_dataset_() noexcept { + dataset_count_ = 0u; + dataset_head_ = 0u; + } + + // --------------------------------------------------------------- + // Members + // --------------------------------------------------------------- + Layer0 layer0_{}; + Layer1 layer1_{}; + Layer2 layer2_{}; + Layer3 layer3_{}; + + std::array input_{}; + std::array output_{}; + + std::array ds_features_{}; + std::array ds_labels_{}; + std::size_t dataset_count_ = 0u; + std::size_t dataset_head_ = 0u; + + std::array flat_weight_buf_{}; + std::array loss_history_{}; + std::size_t loss_history_count_ = 0u; + + Rng rng_; +}; + +// MLP satisfies the MLEngine concept. We keep a static_assert in the test +// suite (test_mlp_concept_satisfied) — see test_mlp_init.cpp. + +} // namespace nisps::ml diff --git a/nisps/ml/rl.hpp b/nisps/ml/rl.hpp new file mode 100644 index 0000000..0fb7886 --- /dev/null +++ b/nisps/ml/rl.hpp @@ -0,0 +1,112 @@ +// nisps/ml/rl.hpp — reinforcement-style weight perturbation primitives. +// +// `move_weights` and `draw_weights` are the two playground RL operations +// applied to network weights: +// - `draw_weights(spread)`: re-randomize all weights using the spread- +// aware scale (see init.hpp). This is the "thumbs-down really hard" +// button. +// - `move_weights(speed, spread, output_pin_mask)`: per-call weight decay +// followed by gaussian noise injection. This is the "thumbs-down a +// little" feedback. The pin mask, if supplied, freezes weights feeding +// specific OUTPUT nodes — only on the FINAL layer. +// +// PARITY CONTRACT WITH LEGACY JS PLAYGROUND +// These match `playground/_archive/js/nisps/mlp.js::moveWeights` and +// `drawWeights` semantically. Specifically: +// * Per-layer scale: layer_scale = (1 - spread) + spread / sqrt(fan_in) +// * Weight decay: each weight gets multiplied by (1 - 0.1 * spread) +// BEFORE noise is added. spread=0 ⇒ no decay, +// spread=1 ⇒ ~10% decay per call. +// * Noise: gaussian via rng.next_float_gaussian(speed * +// layer_scale). The legacy code uses a sum-of-three- +// uniforms gaussian, and our `Rng` matches that +// shape (see core/rng.hpp). +// * Output pin mask: 1 byte per output node. If mask[i]==1, all weights +// AND the bias feeding output node i in the final +// layer are skipped. Hidden layers are unaffected by +// the mask. +// * Biases: included in the noise injection, unaffected by +// decay (matching the legacy JS behavior — decay was +// applied to weights only there). +// +// The functions are layer-scoped; the MLP class (mlp.hpp) iterates through +// its layers, calling these for each. + +#pragma once + +#include +#include +#include +#include + +#include "../core/rng.hpp" +#include "init.hpp" + +namespace nisps::ml { + +// Per-call weight-decay coefficient. spread=0 ⇒ 1.0 (no decay), spread=1 ⇒ +// 0.9 (10% decay per call). See architecture spread §3 in CLAUDE.md. +inline float weight_decay_factor(float spread) noexcept { + return 1.f - 0.1f * spread; +} + +// Apply RL perturbation to one layer's weights+biases. `is_final_layer` +// gates the application of `output_pin_mask`. `output_pin_mask`, if +// non-empty, must be sized fan_out. +// +// Effects, applied in order per weight w: +// w *= (1 - 0.1 * spread) +// w += gaussian(speed * layer_scale) +// Biases: +// b += gaussian(speed * layer_scale) +// Skip rule on final layer if mask[node]==1: leave w/b untouched. +inline void move_weights_layer(std::span weights, + std::span biases, + std::size_t fan_in, + float speed, + float spread, + bool is_final_layer, + std::span output_pin_mask, + Rng& rng) noexcept { + const float layer_scale = spread_scale(spread, fan_in); + const float noise_stddev = speed * layer_scale; + const float decay = weight_decay_factor(spread); + + const std::size_t fan_out = biases.size(); + + for (std::size_t node = 0; node < fan_out; ++node) { + const bool skip = + is_final_layer && !output_pin_mask.empty() && output_pin_mask[node] != 0u; + + // Bias perturbation. + if (!skip) { + biases[node] += rng.next_float_gaussian(noise_stddev); + } + + // Per-weight decay+noise. We always advance the RNG even on skip + // so that the random stream is independent of pin-mask state — this + // makes parity tests deterministic regardless of which outputs are + // pinned. Otherwise toggling a pin would shift the entire + // downstream noise sequence, which would be surprising. + const std::size_t row_off = node * fan_in; + for (std::size_t j = 0; j < fan_in; ++j) { + const float noise = rng.next_float_gaussian(noise_stddev); + if (!skip) { + float& w = weights[row_off + j]; + w = w * decay + noise; + } + } + } +} + +// Re-randomize one layer's weights using the spread-aware scale. Biases are +// reset to zero, matching `init.hpp::spread_init`. +inline void draw_weights_layer(std::span weights, + std::span biases, + std::size_t fan_in, + float spread, + Rng& rng) noexcept { + spread_init(weights, biases, fan_in, spread, rng); +} + +} // namespace nisps::ml diff --git a/nisps/ml/stats.hpp b/nisps/ml/stats.hpp new file mode 100644 index 0000000..99aaca3 --- /dev/null +++ b/nisps/ml/stats.hpp @@ -0,0 +1,76 @@ +// nisps/ml/stats.hpp — per-layer weight-health diagnostics. +// +// Mirrors the `nisps_mlp_get_layer_stats` WASM binding semantics: for each +// layer, return mean(|w|), max(|w|), dead fraction, and saturating fraction. +// These are used by the playground "weight health" panel. +// +// Definitions: +// - mean_abs: average of |w| across all weights+biases of the layer +// - max_abs: max of |w| across all weights+biases +// - dead_frac: fraction of weights with |w| < kDeadThresh +// - saturating_frac: fraction with |w| > kSaturatingThresh +// +// Thresholds are chosen to match the JS engine's heuristic targets: +// dead: 0.001 (a weight smaller than this contributes negligibly) +// sat: 3.0 (weights larger than this typically push sigmoid to its +// rails — a smell during normal training but normal during +// RL exploration) +// +// We intentionally include biases in mean/max but not in dead/sat fractions +// — biases are scarce relative to weights and including them would skew the +// fractions, especially at startup when biases are zero. + +#pragma once + +#include +#include +#include + +namespace nisps::ml { + +inline constexpr float kDeadThresh = 0.001f; +inline constexpr float kSaturatingThresh = 3.0f; + +struct LayerStats { + float mean_abs = 0.f; + float max_abs = 0.f; + float dead_frac = 0.f; + float saturating_frac = 0.f; +}; + +inline LayerStats compute_layer_stats(std::span weights, + std::span biases) noexcept { + LayerStats out; + const std::size_t total = weights.size() + biases.size(); + if (total == 0u) return out; + + float sum = 0.f; + float max_abs = 0.f; + std::size_t dead = 0u; + std::size_t sat = 0u; + + for (float w : weights) { + const float a = std::fabs(w); + sum += a; + if (a > max_abs) max_abs = a; + if (a < kDeadThresh) ++dead; + if (a > kSaturatingThresh) ++sat; + } + for (float b : biases) { + const float a = std::fabs(b); + sum += a; + if (a > max_abs) max_abs = a; + } + + out.mean_abs = sum / static_cast(total); + out.max_abs = max_abs; + + const float w_count = static_cast(weights.size()); + if (w_count > 0.f) { + out.dead_frac = static_cast(dead) / w_count; + out.saturating_frac = static_cast(sat) / w_count; + } + return out; +} + +} // namespace nisps::ml diff --git a/nisps/ml/training.hpp b/nisps/ml/training.hpp new file mode 100644 index 0000000..173fbec --- /dev/null +++ b/nisps/ml/training.hpp @@ -0,0 +1,37 @@ +// nisps/ml/training.hpp — gradient-clipping helper and per-layer SGD update. +// +// The MLP class owns the training loop because it knows the dataset layout +// and weight buffers. This header provides: +// - kGradClip: ±10.0, matches the legacy firmware/Layer.h clamp. +// - clip_gradient(): scalar clipper, applied to each accumulated gradient +// before the weight update. +// +// The full training loop (forward, loss, backprop, weight update) is +// implemented inline in mlp.hpp because everything it touches is either a +// member array or layer-templated. Splitting it across translation units +// would require type-erasing the layers, which we don't want. +// +// Optimizer choice: this MVP ships SGD only. RMSProp is planned (the legacy +// firmware uses it for `TrainBatch`) but the playground was using plain SGD +// until very recently and the XOR-convergence benchmark in the test suite +// is the clearer target. RMSProp can land as a follow-up — see the bd +// issue notes. For now `train()` is SGD with optional sample-weight +// scaling and gradient clipping. + +#pragma once + +#include "../core/perf.hpp" + +namespace nisps::ml { + +// Per-element gradient clip threshold. Matches the legacy firmware's +// gradientClipValue in Layer.h::ApplyAccumulatedGradients. +inline constexpr float kGradClip = 10.f; + +NISPS_FORCE_INLINE float clip_gradient(float g) noexcept { + if (g > kGradClip) return kGradClip; + if (g < -kGradClip) return -kGradClip; + return g; +} + +} // namespace nisps::ml diff --git a/tests/cpp/test_mlp_inference.cpp b/tests/cpp/test_mlp_inference.cpp new file mode 100644 index 0000000..c9fe5d8 --- /dev/null +++ b/tests/cpp/test_mlp_inference.cpp @@ -0,0 +1,114 @@ +// tests/cpp/test_mlp_inference.cpp — verify that the forward pass computes +// the expected linear+activation chain. We construct an MLP with known +// weights via set_weights(), run inference, and compare against a hand- +// computed result. + +#include +#include + +#include "test_helpers.hpp" + +#include "../../nisps/ml/mlp.hpp" + +namespace { + +using TinyMLP = nisps::ml::MLP<2, 2, 2, 2, 1, 8, 32>; + +// Manually compute the forward pass using the same activation rules as the +// MLP. ReLU on hidden layers (leaky 0.01), sigmoid on output. +float manual_forward(float x0, float x1, + const std::array& w0, const std::array& b0, + const std::array& w1, const std::array& b1, + const std::array& w2, const std::array& b2, + const std::array& w3, const std::array& b3) { + auto leaky_relu = [](float v) { return v > 0.f ? v : 0.01f * v; }; + auto sig = [](float v) { + if (v > 40.f) return 1.f; + if (v < -40.f) return 0.f; + return 1.f / (1.f + std::exp(-v)); + }; + + // Layer 0: 2 → 2, ReLU + float a0_0 = leaky_relu(b0[0] + w0[0] * x0 + w0[1] * x1); + float a0_1 = leaky_relu(b0[1] + w0[2] * x0 + w0[3] * x1); + // Layer 1: 2 → 2, ReLU + float a1_0 = leaky_relu(b1[0] + w1[0] * a0_0 + w1[1] * a0_1); + float a1_1 = leaky_relu(b1[1] + w1[2] * a0_0 + w1[3] * a0_1); + // Layer 2: 2 → 2, ReLU + float a2_0 = leaky_relu(b2[0] + w2[0] * a1_0 + w2[1] * a1_1); + float a2_1 = leaky_relu(b2[1] + w2[2] * a1_0 + w2[3] * a1_1); + // Layer 3: 2 → 1, Sigmoid + return sig(b3[0] + w3[0] * a2_0 + w3[1] * a2_1); +} + +NISPS_TEST(mlp_forward_matches_handcomputed) { + TinyMLP m(0ull); + + // Set weights to a known pattern: first weight 0.5, second 0.25, etc. + // (we construct a flat buffer matching the documented layout). + constexpr std::size_t WC = TinyMLP::weight_count(); + std::array flat{}; + // Weights layer by layer: + // L0: 4 weights; L1: 4; L2: 4; L3: 2 → 14 weights + // Biases: L0: 2; L1: 2; L2: 2; L3: 1 → 7 biases + // Total: 21 (sanity check) + NISPS_EXPECT(WC == 21u); + + float val = 0.1f; + for (std::size_t i = 0; i < WC; ++i) { + flat[i] = val; + val += 0.05f; + if (val > 0.7f) val = -0.6f; + } + m.set_weights(std::span(flat)); + + // Extract layer slices for the manual check. + std::array w0{flat[0], flat[1], flat[2], flat[3]}; + std::array w1{flat[4], flat[5], flat[6], flat[7]}; + std::array w2{flat[8], flat[9], flat[10], flat[11]}; + std::array w3{flat[12], flat[13]}; + std::array b0{flat[14], flat[15]}; + std::array b1{flat[16], flat[17]}; + std::array b2{flat[18], flat[19]}; + std::array b3{flat[20]}; + + const float x0 = 0.3f, x1 = 0.7f; + m.set_input(0, x0); + m.set_input(1, x1); + m.process(); + auto out = m.outputs(); + NISPS_EXPECT(out.size() == 1u); + + const float expected = manual_forward(x0, x1, w0, b0, w1, b1, w2, b2, w3, b3); + NISPS_EXPECT_NEAR(out[0], expected, 1e-5); +} + +NISPS_TEST(mlp_forward_outputs_in_unit_range) { + // Sigmoid output guarantees [0, 1]. + using M = nisps::ml::MLP<3, 10, 10, 14, 126, 16, 64>; + M m(42ull); + m.draw_weights(0.5f); + for (int t = 0; t < 100; ++t) { + m.set_input(0, static_cast(t) * 0.013f); + m.set_input(1, static_cast(t) * 0.027f - 0.5f); + m.set_input(2, static_cast(t) * 0.041f); + m.process(); + for (float v : m.outputs()) { + NISPS_EXPECT(v >= 0.f); + NISPS_EXPECT(v <= 1.f); + } + } +} + +NISPS_TEST(mlp_set_input_out_of_range_silently_ignored) { + TinyMLP m(0ull); + // Should not crash; index 99 simply does nothing. + m.set_input(99u, 1.234f); + m.set_input(0, 0.5f); + m.set_input(1, 0.5f); + m.process(); + auto out = m.outputs(); + NISPS_EXPECT(out.size() == 1u); +} + +} // namespace diff --git a/tests/cpp/test_mlp_init.cpp b/tests/cpp/test_mlp_init.cpp new file mode 100644 index 0000000..1779856 --- /dev/null +++ b/tests/cpp/test_mlp_init.cpp @@ -0,0 +1,98 @@ +// tests/cpp/test_mlp_init.cpp — exercises MLP construction, deterministic +// seeding, the spread parameter, and the static_assert that the class +// satisfies the MLEngine concept. + +#include "test_helpers.hpp" + +#include "../../nisps/core/concepts.hpp" +#include "../../nisps/ml/mlp.hpp" + +namespace { + +// Compact alias used across the test suite. +using SmallMLP = nisps::ml::MLP<2, 8, 8, 8, 4, 16, 64>; + +// Hard ground-truth: the class satisfies MLEngine. Compile-time check. +static_assert(nisps::MLEngine, + "MLP<...> must satisfy nisps::MLEngine concept"); + +NISPS_TEST(mlp_same_seed_same_init_weights) { + SmallMLP a(123ull); + SmallMLP b(123ull); + + auto wa = a.get_weights(); + auto wb = b.get_weights(); + NISPS_EXPECT(wa.size() == wb.size()); + for (std::size_t i = 0; i < wa.size(); ++i) { + NISPS_EXPECT(wa[i] == wb[i]); + } +} + +NISPS_TEST(mlp_diff_seed_diff_init_weights) { + SmallMLP a(1ull); + SmallMLP b(2ull); + + auto wa = a.get_weights(); + auto wb = b.get_weights(); + int distinct = 0; + for (std::size_t i = 0; i < wa.size(); ++i) { + if (wa[i] != wb[i]) ++distinct; + } + // With ~200+ weights, almost all should differ. + NISPS_EXPECT(distinct > static_cast(wa.size()) / 2); +} + +NISPS_TEST(mlp_seed_method_resets_rng) { + SmallMLP a(5ull); + a.seed(42ull); + a.draw_weights(0.5f); + + SmallMLP b(99ull); + b.seed(42ull); + b.draw_weights(0.5f); + + auto wa = a.get_weights(); + auto wb = b.get_weights(); + for (std::size_t i = 0; i < wa.size(); ++i) { + NISPS_EXPECT(wa[i] == wb[i]); + } +} + +NISPS_TEST(mlp_weight_count_matches_topology) { + using M = nisps::ml::MLP<3, 10, 10, 14, 126>; + // Layer fan_in*fan_out: 3*10 + 10*10 + 10*14 + 14*126 = 30+100+140+1764 = 2034 + // Biases: 10+10+14+126 = 160 + // Total: 2194 + NISPS_EXPECT(M::weight_count() == 2194u); +} + +NISPS_TEST(mlp_spread_zero_uniform_in_minus_one_one_range) { + // spread=0 → weights drawn from U[-1,1]. The max |w| should be ≤ 1. + SmallMLP m(7ull); + m.draw_weights(0.f); + auto w = m.get_weights(); + float maxabs = 0.f; + for (float v : w) { + const float a = v >= 0.f ? v : -v; + if (a > maxabs) maxabs = a; + } + NISPS_EXPECT(maxabs <= 1.f); +} + +NISPS_TEST(mlp_spread_one_xavier_smaller_range) { + // spread=1 → weights scaled by 1/sqrt(fan_in). For fan_in≥2, weights + // should be strictly smaller in magnitude than the spread=0 case. + SmallMLP m(7ull); + m.draw_weights(1.f); + auto w = m.get_weights(); + float maxabs = 0.f; + for (float v : w) { + const float a = v >= 0.f ? v : -v; + if (a > maxabs) maxabs = a; + } + // 1/sqrt(2) ≈ 0.707 — the smallest fan_in is 2 (input). So max possible + // is ≈ 0.707 (drawn from rng_signed, which is < 1). + NISPS_EXPECT(maxabs < 1.f); +} + +} // namespace diff --git a/tests/cpp/test_mlp_loss.cpp b/tests/cpp/test_mlp_loss.cpp new file mode 100644 index 0000000..bd6420a --- /dev/null +++ b/tests/cpp/test_mlp_loss.cpp @@ -0,0 +1,125 @@ +// tests/cpp/test_mlp_loss.cpp — regression test for the meml-ues +// double-scaling bug. +// +// Setup: one training example (label = current network output exactly). +// Expected loss = 0. +// +// Then: one example with label = output + small_delta. +// Expected loss = (1/NOut) * sum(small_delta^2). NOT scaled additionally +// by 1/N (here N=1 so any extra 1/N would still be 1, but the per-element +// derivative path is what would expose double-scaling). +// +// We compare the value reported by `train()` after a single iteration +// (with lr=0 — so no weight update — but the loss should be reported) +// against the expected formula. + +#include +#include + +#include "test_helpers.hpp" + +#include "../../nisps/ml/loss.hpp" +#include "../../nisps/ml/mlp.hpp" + +namespace { + +// Direct test of mse_per_sample. +NISPS_TEST(mse_per_sample_matches_formula) { + std::array label = {0.5f, 0.25f, 0.75f, 0.f}; + std::array pred = {0.4f, 0.30f, 0.70f, 0.1f}; + std::array deriv{}; + + const float l = nisps::ml::mse_per_sample( + std::span(label), + std::span(pred), + std::span(deriv)); + + // (0.1^2 + 0.05^2 + 0.05^2 + 0.1^2) / 4 = (0.01 + 0.0025 + 0.0025 + 0.01) / 4 = 0.025/4 + const float expected = (0.01f + 0.0025f + 0.0025f + 0.01f) * 0.25f; + NISPS_EXPECT_NEAR(l, expected, 1e-7); + + // Derivatives: -(2/N) * (label - pred). For element 0: -(2/4) * (0.5 - 0.4) = -0.05 + NISPS_EXPECT_NEAR(deriv[0], -0.05f, 1e-7); + NISPS_EXPECT_NEAR(deriv[3], 0.05f, 1e-7); // -(2/4) * (0 - 0.1) = +0.05 +} + +NISPS_TEST(mlp_train_loss_not_double_scaled_single_sample) { + // 1 example. The reported epoch_loss should equal the per-sample MSE + // exactly (with sample_weight=1/N=1.0). If the legacy double-scaling + // bug crept in, we'd see it scaled by an extra 1/1 = 1.0 (invisible) + // — so we use TWO samples and check the relationship. + using M = nisps::ml::MLP<1, 2, 2, 2, 2, 4, 8>; + M m(11ull); + m.draw_weights(0.5f); + + std::array f1{0.3f}; + std::array l1{0.7f, 0.2f}; + std::array f2{0.6f}; + std::array l2{0.1f, 0.9f}; + m.add_example(std::span(f1), std::span(l1)); + m.add_example(std::span(f2), std::span(l2)); + + // Compute the expected loss manually: average of per-sample MSE. + // Run inference on each, get pred, compute MSE manually. + auto compute_mse = [&](std::array& f, std::array& lab) { + m.set_input(0, f[0]); + m.process(); + auto out = m.outputs(); + float sse = 0.f; + for (std::size_t j = 0; j < 2u; ++j) { + const float d = lab[j] - out[j]; + sse += d * d; + } + return sse * 0.5f; // (1/NOut) sum of squared diff + }; + const float mse1 = compute_mse(f1, l1); + const float mse2 = compute_mse(f2, l2); + const float expected_epoch_loss = 0.5f * (mse1 + mse2); // avg over 2 samples + + // train() with lr=0 leaves weights unchanged but reports the per-iter + // loss. We use max_iter=1 to grab exactly one epoch_loss. + const float reported = m.train(/*lr=*/0.f, /*max_iter=*/1u, /*min_err=*/-1.f); + + NISPS_EXPECT_NEAR(reported, expected_epoch_loss, 1e-5); + + // If the legacy bug were present, reported would be expected/2 (extra + // 1/N multiplication). That would fail at 1e-5 tolerance for any + // non-tiny mse. Sanity-confirm: + NISPS_EXPECT(reported > expected_epoch_loss * 0.9f); +} + +NISPS_TEST(mlp_train_with_sample_weights_uses_them) { + using M = nisps::ml::MLP<1, 2, 2, 2, 2, 4, 8>; + M m(13ull); + m.draw_weights(0.5f); + + std::array f1{0.1f}; + std::array l1{0.5f, 0.5f}; + std::array f2{0.9f}; + std::array l2{0.5f, 0.5f}; + m.add_example(std::span(f1), std::span(l1)); + m.add_example(std::span(f2), std::span(l2)); + + // Concentrate all weight on sample 0. + std::array sw{1.f, 0.f}; + + auto compute_mse = [&](std::array& f, std::array& lab) { + m.set_input(0, f[0]); + m.process(); + auto out = m.outputs(); + float sse = 0.f; + for (std::size_t j = 0; j < 2u; ++j) { + const float d = lab[j] - out[j]; + sse += d * d; + } + return sse * 0.5f; + }; + const float mse1 = compute_mse(f1, l1); + + const float reported = m.train(0.f, 1u, -1.f, std::span(sw)); + // Under uniform weights this would be 0.5*(mse1+mse2). With sw={1,0} + // it should equal mse1. + NISPS_EXPECT_NEAR(reported, mse1, 1e-5); +} + +} // namespace diff --git a/tests/cpp/test_mlp_rl.cpp b/tests/cpp/test_mlp_rl.cpp new file mode 100644 index 0000000..39fb62c --- /dev/null +++ b/tests/cpp/test_mlp_rl.cpp @@ -0,0 +1,133 @@ +// tests/cpp/test_mlp_rl.cpp — verify draw_weights and move_weights +// behavior: +// - draw_weights with spread=0 vs spread=1 produces different weight +// magnitude regimes. +// - move_weights perturbs all weights in unmasked layers. +// - move_weights with output_pin_mask preserves the corresponding +// final-layer weight rows AND biases exactly. + +#include +#include + +#include "test_helpers.hpp" + +#include "../../nisps/ml/mlp.hpp" + +namespace { + +using SmallMLP = nisps::ml::MLP<2, 4, 4, 4, 6, 8, 32>; + +NISPS_TEST(mlp_draw_weights_spread_zero_full_range) { + SmallMLP m(0ull); + m.draw_weights(0.f); + auto w = m.get_weights(); + float maxabs = 0.f; + for (float v : w) { + const float a = v >= 0.f ? v : -v; + if (a > maxabs) maxabs = a; + } + // spread=0 ⇒ U[-1,1] (no Xavier scale). Most weights will fall in + // (0.5, 1.0) band — we expect at least one above 0.5. + NISPS_EXPECT(maxabs > 0.5f); +} + +NISPS_TEST(mlp_draw_weights_spread_one_xavier_compressed) { + SmallMLP m(0ull); + m.draw_weights(1.f); + auto w = m.get_weights(); + float maxabs = 0.f; + for (float v : w) { + const float a = v >= 0.f ? v : -v; + if (a > maxabs) maxabs = a; + } + // Xavier scale: 1/sqrt(fan_in). Smallest fan_in = 2 → scale ~0.707. + // Largest weight magnitude bounded by that. + NISPS_EXPECT(maxabs < 0.71f); +} + +NISPS_TEST(mlp_move_weights_changes_unpinned_weights) { + SmallMLP m(99ull); + m.draw_weights(0.5f); + + // Snapshot weights. + std::array before{}; + { + auto w = m.get_weights(); + for (std::size_t i = 0; i < w.size(); ++i) before[i] = w[i]; + } + + m.move_weights(0.1f, 0.5f); + + auto after = m.get_weights(); + int distinct = 0; + for (std::size_t i = 0; i < after.size(); ++i) { + if (before[i] != after[i]) ++distinct; + } + // Most weights should have changed (gaussian noise + decay). + NISPS_EXPECT(distinct > static_cast(after.size()) / 2); +} + +NISPS_TEST(mlp_move_weights_pin_mask_skips_final_outputs) { + SmallMLP m(33ull); + m.draw_weights(0.5f); + + auto w_before = m.get_weights(); + std::array before{}; + for (std::size_t i = 0; i < w_before.size(); ++i) before[i] = w_before[i]; + + // NOut = 6. Pin nodes 0, 2, 4. + std::array mask{1, 0, 1, 0, 1, 0}; + m.move_weights(0.1f, 0.5f, std::span(mask)); + + auto after = m.get_weights(); + + // Layout: weights = [L0(8) L1(16) L2(16) L3(24)] then biases [L0(4) L1(4) L2(4) L3(6)]. + constexpr std::size_t L0_W = 2*4; // 8 + constexpr std::size_t L1_W = 4*4; // 16 + constexpr std::size_t L2_W = 4*4; // 16 + constexpr std::size_t L3_W = 4*6; // 24 + constexpr std::size_t L3_W_OFF = L0_W + L1_W + L2_W; + constexpr std::size_t BIAS_OFF = L0_W + L1_W + L2_W + L3_W; + constexpr std::size_t L3_B_OFF = BIAS_OFF + 4 + 4 + 4; // 76 + 0 → biases start + + // Final layer weights row-major: [node*4 + j] for j ∈ [0,4). + // Pinned nodes 0, 2, 4 → rows 0, 2, 4 should be preserved. + for (std::size_t node : {0u, 2u, 4u}) { + for (std::size_t j = 0; j < 4u; ++j) { + const std::size_t idx = L3_W_OFF + node * 4u + j; + NISPS_EXPECT(before[idx] == after[idx]); + } + // Bias too. + NISPS_EXPECT(before[L3_B_OFF + node] == after[L3_B_OFF + node]); + } + // Unpinned nodes 1, 3, 5 → rows changed (most weights distinct). + int unpinned_changed = 0; + for (std::size_t node : {1u, 3u, 5u}) { + for (std::size_t j = 0; j < 4u; ++j) { + const std::size_t idx = L3_W_OFF + node * 4u + j; + if (before[idx] != after[idx]) ++unpinned_changed; + } + } + // 12 unpinned weights; gaussian noise w/ stddev > 0 → almost all change. + NISPS_EXPECT(unpinned_changed >= 10); +} + +NISPS_TEST(mlp_draw_weights_clears_grad_accumulators) { + // After draw_weights, calling train() shouldn't see stale gradients. + // Smoke test: draw → train → loss should drop normally. + using M = nisps::ml::MLP<2, 4, 4, 4, 1, 4, 64>; + M m(0ull); + m.draw_weights(0.5f); + std::array f{0.3f, 0.7f}; + std::array l{0.5f}; + m.add_example(std::span(f), std::span(l)); + const float l1 = m.train(0.5f, 50u, -1.f); + m.draw_weights(0.5f); + const float l2 = m.train(0.5f, 50u, -1.f); + // Both should be finite and >= 0. We don't assert about ordering; the + // point is that no NaNs leak through stale grads. + NISPS_EXPECT(l1 >= 0.f && l1 < 100.f); + NISPS_EXPECT(l2 >= 0.f && l2 < 100.f); +} + +} // namespace diff --git a/tests/cpp/test_mlp_serialize.cpp b/tests/cpp/test_mlp_serialize.cpp new file mode 100644 index 0000000..0b2f34b --- /dev/null +++ b/tests/cpp/test_mlp_serialize.cpp @@ -0,0 +1,146 @@ +// tests/cpp/test_mlp_serialize.cpp — round-trip get_weights → set_weights. +// +// Serializing weights, restoring them in a fresh MLP, and running inference +// must produce bit-identical outputs. This is the gate for cross-platform +// state transfer (firmware ↔ browser). + +#include +#include + +#include "test_helpers.hpp" + +#include "../../nisps/ml/mlp.hpp" + +namespace { + +using TestMLP = nisps::ml::MLP<3, 8, 8, 10, 5, 16, 32>; + +NISPS_TEST(mlp_get_set_weights_roundtrip_preserves_inference) { + TestMLP a(123ull); + a.draw_weights(0.4f); + + // Snapshot weights. + auto w = a.get_weights(); + std::vector snap(w.begin(), w.end()); + + TestMLP b(0ull); // different seed → different starting weights + b.set_weights(std::span(snap.data(), snap.size())); + + // Inference on identical input should produce identical outputs. + const float in[3] = {0.2f, -0.3f, 0.5f}; + for (std::size_t i = 0; i < 3u; ++i) { + a.set_input(i, in[i]); + b.set_input(i, in[i]); + } + a.process(); + b.process(); + + auto oa = a.outputs(); + auto ob = b.outputs(); + NISPS_EXPECT(oa.size() == ob.size()); + for (std::size_t i = 0; i < oa.size(); ++i) { + NISPS_EXPECT(oa[i] == ob[i]); + } +} + +NISPS_TEST(mlp_weight_count_matches_get_weights_size) { + TestMLP m(0ull); + auto w = m.get_weights(); + NISPS_EXPECT(w.size() == TestMLP::weight_count()); +} + +NISPS_TEST(mlp_get_weights_contains_layer_concatenation) { + // Verify the documented flat layout: weights all layers, then biases + // all layers. We do this by setting a known pattern via set_weights + // and reading it back. + TestMLP m(0ull); + constexpr std::size_t WC = TestMLP::weight_count(); + std::vector pattern(WC); + for (std::size_t i = 0; i < WC; ++i) { + pattern[i] = static_cast(i) * 0.001f - 0.5f; + } + m.set_weights(std::span(pattern.data(), pattern.size())); + auto w = m.get_weights(); + NISPS_EXPECT(w.size() == WC); + for (std::size_t i = 0; i < WC; ++i) { + NISPS_EXPECT_NEAR(w[i], pattern[i], 1e-7); + } +} + +NISPS_TEST(mlp_set_weights_too_short_is_ignored) { + TestMLP m(7ull); + auto before = m.get_weights(); + std::vector snap(before.begin(), before.end()); + + // Pass a too-short buffer. + std::array tiny{1.f, 2.f}; + m.set_weights(std::span(tiny)); + + auto after = m.get_weights(); + // Weights should be unchanged. + for (std::size_t i = 0; i < snap.size(); ++i) { + NISPS_EXPECT(snap[i] == after[i]); + } +} + +NISPS_TEST(mlp_eval_loss_does_not_modify_state) { + TestMLP m(5ull); + m.draw_weights(0.5f); + + std::array f{0.1f, 0.2f, 0.3f}; + std::array l{0.5f, 0.5f, 0.5f, 0.5f, 0.5f}; + m.add_example(std::span(f), std::span(l)); + + // Snapshot weights. + auto w_before = m.get_weights(); + std::vector snap(w_before.begin(), w_before.end()); + + const float L = m.eval_loss(); + NISPS_EXPECT(L >= 0.f); + + auto w_after = m.get_weights(); + for (std::size_t i = 0; i < snap.size(); ++i) { + NISPS_EXPECT(snap[i] == w_after[i]); + } +} + +NISPS_TEST(mlp_infer_batch_matches_individual_inference) { + TestMLP m(42ull); + m.draw_weights(0.5f); + + constexpr std::size_t N = 5; + constexpr std::size_t NI = 3; + constexpr std::size_t NO = 5; + std::array points{}; + for (std::size_t i = 0; i < N * NI; ++i) { + points[i] = static_cast(i) * 0.07f - 0.3f; + } + std::array outs{}; + m.infer_batch(std::span(points), std::span(outs)); + + // Compare against individual inference. + for (std::size_t i = 0; i < N; ++i) { + for (std::size_t j = 0; j < NI; ++j) { + m.set_input(j, points[i * NI + j]); + } + m.process(); + auto o = m.outputs(); + for (std::size_t j = 0; j < NO; ++j) { + NISPS_EXPECT(o[j] == outs[i * NO + j]); + } + } +} + +NISPS_TEST(mlp_layer_stats_reasonable_after_init) { + TestMLP m(0ull); + m.draw_weights(0.5f); + for (std::size_t L = 0; L < 4; ++L) { + const auto s = m.layer_stats(L); + NISPS_EXPECT(s.mean_abs > 0.f); + NISPS_EXPECT(s.max_abs >= s.mean_abs); + NISPS_EXPECT(s.dead_frac >= 0.f && s.dead_frac <= 1.f); + NISPS_EXPECT(s.saturating_frac >= 0.f && s.saturating_frac <= 1.f); + } +} + +} // namespace diff --git a/tests/cpp/test_mlp_training.cpp b/tests/cpp/test_mlp_training.cpp new file mode 100644 index 0000000..ede52c5 --- /dev/null +++ b/tests/cpp/test_mlp_training.cpp @@ -0,0 +1,92 @@ +// tests/cpp/test_mlp_training.cpp — convergence test on XOR. +// +// XOR is the classic minimum non-linear problem: a 2-layer linear net +// cannot solve it; an MLP with one hidden layer (and a non-linear +// activation) can. Our 4-layer MLP with sigmoid output is more than enough. +// +// We check loss < 0.01 within a generous iteration budget. If this test +// regresses to taking >1000 iterations, something is wrong with the +// gradient or weight-update path. + +#include + +#include "test_helpers.hpp" + +#include "../../nisps/ml/mlp.hpp" + +namespace { + +NISPS_TEST(mlp_xor_converges) { + // Modest network: 2 inputs, [4, 4, 4] hidden, 1 output. + using M = nisps::ml::MLP<2, 4, 4, 4, 1, 4, 1024>; + M m(7ull); + m.draw_weights(1.f); // Xavier-ish; needed for sigmoid output to start reasonable + + // XOR truth table. + std::array, 4> X = {{ + {0.f, 0.f}, + {0.f, 1.f}, + {1.f, 0.f}, + {1.f, 1.f}, + }}; + std::array, 4> Y = {{ + {0.f}, + {1.f}, + {1.f}, + {0.f}, + }}; + + for (std::size_t i = 0; i < 4u; ++i) { + m.add_example(std::span(X[i]), std::span(Y[i])); + } + NISPS_EXPECT(m.example_count() == 4u); + + // Train. Higher LR + more iterations is fine — the test is "did it + // converge AT ALL within a generous budget". + const float final_loss = m.train(/*lr=*/0.5f, /*max_iter=*/2000u, /*min_err=*/0.01f); + NISPS_EXPECT(final_loss < 0.01f); + + // Sanity: outputs should be near labels for each input. + for (std::size_t i = 0; i < 4u; ++i) { + m.set_input(0, X[i][0]); + m.set_input(1, X[i][1]); + m.process(); + const float o = m.outputs()[0]; + NISPS_EXPECT_NEAR(o, Y[i][0], 0.2); + } + + // Loss history should record at least one entry. + NISPS_EXPECT(m.loss_history().size() >= 1u); +} + +NISPS_TEST(mlp_train_with_no_examples_returns_zero) { + using M = nisps::ml::MLP<2, 4, 4, 4, 1, 4, 8>; + M m(0ull); + const float loss = m.train(0.5f, 100u, 0.001f); + NISPS_EXPECT(loss == 0.f); +} + +NISPS_TEST(mlp_clear_examples_works) { + using M = nisps::ml::MLP<2, 4, 4, 4, 1, 4, 8>; + M m(0ull); + std::array f{0.f, 1.f}; + std::array l{0.5f}; + m.add_example(std::span(f), std::span(l)); + NISPS_EXPECT(m.example_count() == 1u); + m.clear_examples(); + NISPS_EXPECT(m.example_count() == 0u); +} + +NISPS_TEST(mlp_dataset_ring_buffer_evicts_oldest) { + // NMaxExamples=4, add 6 examples; oldest 2 should be evicted. + using M = nisps::ml::MLP<1, 2, 2, 2, 1, 4, 8>; + M m(0ull); + for (int i = 0; i < 6; ++i) { + std::array f{static_cast(i)}; + std::array l{static_cast(i) * 0.1f}; + m.add_example(std::span(f), std::span(l)); + } + NISPS_EXPECT(m.example_count() == 4u); +} + +} // namespace