Add sequential CRHMC relative-volume updating

dingo/matlab/volume_updating/LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2026 GeomScale
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

dingo/matlab/volume_updating/default_params.m

@@ -0,0 +1,29 @@
+function params = default_params()
+    %DEFAULT_PARAMS Return the default parameter struct for relative_volume.
+    %
+    % Usage:
+    %   params = default_params();
+    %   params = Setfield(default_params(), user_overrides);
+
+    % ---- Annealing schedule parameters ----
+    params.r     = 0.1;   % target lower bound for each volume ratio
+    params.delta = 0.05;   % half-width of target ratio interval [r, r+delta]
+    params.alpha = 0.20;   % significance level for U/L t-tests
+    params.nu    = 10;     % number of sublists for U/L tests
+    params.N_utest = 320;  % points per sublist (total = nu * N_utest)
+    params.window_size = 5000;  % size of the sliding window (We now set it in the new_window()!!)
+
+    % ---- Convergence parameters ----
+    params.epsilon = 0.10;             % target relative error for final estimate
+    params.max_annealing_iter = 200;    % max Phase 1 construction iterations
+    params.max_bisect_iter = 100;       % max binary search iterations per body
+
+    % ---- Sampler parameters ----
+    params.simdLen = 8;                % number of parallel chains per step()
+
+    % ---- t_max estimation ----
+    params.t_max_safety = 1.01;        % inflate sampled t_max by this factor
+
+    % ---- Verbosity ----
+    params.verb = 1;   % 0 = quiet, 1 = normal, 2 = detailed
+end

dingo/matlab/volume_updating/init.m

@@ -0,0 +1,81 @@
+function simdLen = init()
+    %INIT Add PolytopeSamplerMatlab and src/ to the MATLAB path and compile solvers.
+    %
+    %   simdLen = init()
+    %
+    % Returns the maximum usable simdLen (number of parallel SIMD chains).
+    % Tries simdLen=8 first; falls back to simdLen=4 if the MEX solver for 8
+    % cannot be compiled (e.g. no C++ compiler is available).  simdLen=4 has
+    % pre-compiled binaries in PolytopeSamplerMatlab/bin/ for Linux, macOS,
+    % and Windows on x86-64.
+    %
+    % This function is idempotent — subsequent calls return the cached value
+    % without recompiling.  Called automatically by relative_volume() and
+    % run_relative_volume().
+
+    persistent cached_simdLen;
+    if ~isempty(cached_simdLen)
+        simdLen = cached_simdLen;
+        return;
+    end
+
+    src_dir = fileparts(mfilename('fullpath'));
+    root    = fileparts(src_dir);
+
+    % Add src/ so default_params.m, sliding_window.m, ul_test.m are findable
+    addpath(src_dir);
+
+    % The original repository keeps PolytopeSamplerMatlab beside src/.  The
+    % Dingo overlay can be installed in a different layout, so first honor an
+    % explicit environment variable and then check deterministic locations.
+    sampler_root = getenv('DINGO_POLYTOPE_SAMPLER_MATLAB');
+    candidates = {
+        sampler_root, ...
+        fullfile(root, 'PolytopeSamplerMatlab'), ...
+        fullfile(fileparts(fileparts(fileparts(src_dir))), 'PolytopeSamplerMatlab'), ...
+        fullfile(fileparts(fileparts(fileparts(fileparts(src_dir)))), 'PolytopeSamplerMatlab') ...
+    };
+    sampler_root = '';
+    for candidate_idx = 1:numel(candidates)
+        candidate = candidates{candidate_idx};
+        if ~isempty(candidate) && exist(fullfile(candidate, 'code'), 'dir')
+            sampler_root = candidate;
+            break;
+        end
+    end
+    if isempty(sampler_root)
+        error(['PolytopeSamplerMatlab was not found. Set ', ...
+               'DINGO_POLYTOPE_SAMPLER_MATLAB to its repository root.']);
+    end
+    sampler_code = fullfile(sampler_root, 'code');
+    sampler_bin  = fullfile(sampler_root, 'bin');
+
+    addpath(genpath(sampler_code));
+    addpath(sampler_bin);
+
+    % Scalar solver — pre-compiled binaries exist for all platforms
+    compile_solver(0);
+
+    % SIMD solver.  Try simdLen=8 first (best throughput); fall back to
+    % simdLen=4 if the MEX binary cannot be found or compiled.  simdLen=4
+    % has pre-compiled binaries shipped with PolytopeSamplerMatlab.
+    for trial = [8, 4]
+        try
+            compile_solver(trial);
+            cached_simdLen = trial;
+            if trial == 4
+                warning('init:fallback', ...
+                        ['Cannot compile MEX solver for simdLen=8 ' ...
+                         '(no C++ compiler?).  Falling back to simdLen=4 ' ...
+                         'which uses pre-compiled binaries.']);
+            end
+            simdLen = cached_simdLen;
+            return;
+        catch ME
+            if trial == 4
+                rethrow(ME);
+            end
+            % simdLen=8 failed — try simdLen=4
+        end
+    end
+end

dingo/matlab/volume_updating/new_window.m

@@ -0,0 +1,82 @@
+function window = new_window(params, dim, eps_i)
+    %NEW_WINDOW Create a sliding window convergence tracker.
+    %
+    % Input:
+    %   dim   - effective dimension of the polytope
+    %   eps_i - per-ratio error budget
+    %
+    % Output:
+    %   window - struct with fields:
+    %     .W          window size = ceil(4*dim^2 + 500)
+    %     .min_val    current minimum in the window
+    %     .max_val    current maximum in the window
+    %     .min_index  index of min in circular buffer
+    %     .max_index  index of max in circular buffer
+    %     .last_W     circular buffer of length W
+    %     .index      write pointer (1-based)
+    %     .converged  logical flag
+    %     .eps_i      per-ratio error budget
+    %
+    % Reference: Cousins & Vempala (2016), Section 3.3.
+    % Adapted from Volume-and-Sampling/Volume.m lines 246-266.
+
+    %window.W = ceil(4 * dim^2 + 500);
+    window.W = ceil(params.simdLen * dim^0.52)+400;
+    %window.W = params.window_size;
+    window.min_val = -realmax;
+    window.max_val = realmax;
+    window.min_index = window.W;
+    window.max_index = window.W;
+    window.last_W = zeros(window.W, 1);
+    window.index = 1;
+    window.converged = false;
+    window.eps_i = eps_i;
+    window.filled = false;
+end
+
+function window = update_window(window, val)
+    %UPDATE_WINDOW Push a new cumulative running average into the window.
+    %
+    % Input:
+    %   window - window struct from new_window
+    %   val    - current cumulative running average (scalar)
+    %
+    % Output:
+    %   window - updated struct with .converged set to true if converged
+    %
+    % Maintains O(1) amortized min/max tracking via lazy eviction:
+    % when the outgoing value was NOT the min or max, updates are O(1).
+    % When the outgoing value WAS the min or max, a full scan recomputes.
+    %
+    % Adapted from Volume-and-Sampling/Volume.m lines 268-301.
+
+    window.last_W(window.index) = val;
+
+    % Update minimum
+    if val <= window.min_val
+        window.min_val = val;
+        window.min_index = window.index;
+    elseif window.min_index == window.index
+        % The outgoing element was the min -- recompute
+        window.filled = true;
+        [window.min_val, window.min_index] = min(window.last_W);
+    end
+
+    % Update maximum
+    if val >= window.max_val
+        window.max_val = val;
+        window.max_index = window.index;
+    elseif window.max_index == window.index
+        % The outgoing element was the max -- recompute
+        window.filled = true;
+        [window.max_val, window.max_index] = max(window.last_W);
+    end
+
+    % Check convergence: relative spread within eps_i/2
+    if (window.max_val - window.min_val) / window.max_val <= window.eps_i / 2 && window.filled
+    %if (max(window.last_W) - min(window.last_W)) / max(window.last_W) <= window.eps_i / 2
+        window.converged = true;
+    end
+    % Advance circular buffer pointer
+    window.index = mod(window.index, window.W) + 1;
+end