matlab_rng.h

Go to the documentation of this file.

/*
 * Copyright (c) 2012-2026, QORE Lab, Imperial College London
 * All rights reserved.
 */
#ifndef LINE_GEN_MATLAB_RNG_H
#define LINE_GEN_MATLAB_RNG_H
 
/**
 * @file
 * @ingroup line_gen
 * MATLAB-faithful random draws for the model generators: the C++ twin of the
 * JAR `jline.gen.MatlabRandom` and the native Python `line_solver.gen.matlab_rng`.
 *
 * A `MatlabRng` seeded with s yields the stream MATLAB yields after `rng(s)`
 * (MT19937 `init_genrand`, numpy's legacy `RandomState(s)`). `rand` is numpy's
 * `random_sample`, `randi(lo, hi)` is `lo-1+ceil((hi-lo+1)*rand)` (a degenerate
 * range still consumes a draw), `randperm(n)` is the stable sort order of
 * `rand(1,n)` and `randfixedsum` is Roger Stafford's original, draw for draw.
 * Checked against MATLAB R2026a on 2026-10-13.
 *
 * Unlike MATLAB's global generator the stream here belongs to its owner, so two
 * generators never share one; within one generator the draw sequence is MATLAB's.
 */
 
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <numeric>
#include <vector>
 
#include "line/opt/de/numpy_random_state.h"
 
namespace line {
namespace gen {
 
/** One MATLAB global random stream, as left by `rng(seed)`. */
class MatlabRng {
public:
    explicit MatlabRng(long long seed = 0) : state_(std::uint64_t(seed) & 0xffffffffULL) {}
 
    /** MATLAB `rng(seed)`. */
    void reseed(long long seed) { state_ = opt::de::NumpyRandomState(std::uint64_t(seed) & 0xffffffffULL); }
 
    /** One MATLAB `rand` draw. */
    double rand() { return state_.random_sample(); }
 
    /** MATLAB `randi([lo, hi])`. */
    int randi(int lo, int hi) { return lo - 1 + int(std::ceil(double(hi - lo + 1) * rand())); }
 
    /** MATLAB `randperm(n)`, returned 0-based. */
    std::vector<std::size_t> randperm(std::size_t n) {
        std::vector<double> u(n);
        for (std::size_t i = 0; i < n; ++i) u[i] = rand();
        std::vector<std::size_t> idx(n);
        std::iota(idx.begin(), idx.end(), std::size_t(0));
        std::stable_sort(idx.begin(), idx.end(),
                         [&u](std::size_t a, std::size_t b) { return u[a] < u[b]; });
        return idx;
    }
 
    /**
     * One column of `randfixedsum(n,1,s,a,b)` (matlab/lib/thirdparty/randfixedsum.m):
     * draws `rand(n-1,1)` twice, then `rand(n,1)` for the permutation.
     */
    std::vector<double> randfixedsum(std::size_t n, double s, double a, double b) {
        const double nn = double(n);
        s = (s - nn * a) / (b - a);
        const int k = std::max(std::min(int(std::floor(s)), int(n) - 1), 0);
        s = std::max(std::min(s, double(k + 1)), double(k));
        std::vector<double> s1(n), s2(n);
        for (std::size_t i = 0; i < n; ++i) {
            s1[i] = s - double(k - int(i));
            s2[i] = double(k + int(n) - int(i)) - s;
        }
        const double tiny = std::numeric_limits<double>::denorm_min();
        std::vector<std::vector<double>> w(n, std::vector<double>(n + 1, 0.0));
        w[0][1] = std::numeric_limits<double>::max();
        std::vector<std::vector<double>> t(n > 1 ? n - 1 : 0, std::vector<double>(n, 0.0));
        for (std::size_t i = 2; i <= n; ++i) {
            const double di = double(i);
            for (std::size_t c = 0; c < i; ++c) {
                const double tmp1 = w[i - 2][c + 1] * s1[c] / di;
                const double tmp2 = w[i - 2][c] * s2[n - i + c] / di;
                w[i - 1][c + 1] = tmp1 + tmp2;
                const double tmp3 = w[i - 1][c + 1] + tiny;
                const bool tmp4 = s2[n - i + c] > s1[c];
                t[i - 2][c] = tmp4 ? tmp2 / tmp3 : 1.0 - tmp1 / tmp3;
            }
        }
        std::vector<double> x(n, 0.0);
        std::vector<double> rt(n > 1 ? n - 1 : 0), rs(n > 1 ? n - 1 : 0);
        for (double& v : rt) v = rand();
        for (double& v : rs) v = rand();
        int j = k + 1;
        double sm = 0.0, pr = 1.0;
        for (std::size_t i = n - 1; i >= 1; --i) {
            const double e = rt[n - i - 1] <= t[i - 1][std::size_t(j - 1)] ? 1.0 : 0.0;
            const double sx = std::pow(rs[n - i - 1], 1.0 / double(i));
            sm += (1.0 - sx) * pr * s / double(i + 1);
            pr *= sx;
            x[n - i - 1] = sm + pr * e;
            s -= e;
            j -= int(e);
        }
        x[n - 1] = sm + pr * s;
        const std::vector<std::size_t> p = randperm(n);
        std::vector<double> out(n);
        for (std::size_t i = 0; i < n; ++i) out[i] = (b - a) * x[p[i]] + a;
        return out;
    }
 
private:
    opt::de::NumpyRandomState state_;
};
 
}  // namespace gen
}  // namespace line
 
#endif  // LINE_GEN_MATLAB_RNG_H