1function walnuts_sampler(N, dt, max_error, target)
2%WALNUTS_SAMPLER Interactive figure: WALNUTS sampling of a 2D target.
3% WALNUTS_SAMPLER(N, DT, MAX_ERROR, TARGET) draws N samples from TARGET
4% ('banana' | 'gaussian' | 'correlated' | 'donut') with the WALNUTS sampler
5% (helpers/walnuts.m) and opens a figure showing the density and the samples.
6%
7% This is the wiring: it selects the target, runs the sampler, evaluates the
8% target on a grid for the heatmap, loads the prebuilt figure app, and sends
9% both. It also handles requests from the figure (figure -> script): change
10% the target, resample, or record an orbit movie. Run walnuts_demo.m (which
11% addpath's helpers/).
13if nargin < 1 || isempty(N); N = 1000; end
14if nargin < 2 || isempty(dt); dt = 0.4; end
15if nargin < 3 || isempty(max_error); max_error = 0.8; end
16if nargin < 4 || isempty(target); target = 'banana'; end
18set_target(target);
19samples = run_chain(N, dt, max_error, target);
20dens = density_grid(target);
22html = fileread(fullfile('app', 'dist', 'index.html'));
23fig = figure;
24gl = uigridlayout(fig, [1 1], 'Padding', [0 0 0 0], ...
25 'RowHeight', {'1x'}, 'ColumnWidth', {'1x'});
26uihtml(gl, 'HTMLSource', html, 'Data', pack_data(samples, dens, N, dt, max_error, target), ...
27 'HTMLEventReceivedFcn', @(src, ev) on_event(src, ev));
28end
30function on_event(src, ev)
31% Figure -> script. The figure owns the current settings and passes them back:
32% 'resample' {n, dt, maxError, target} -> fresh chain; reply 'samples'.
33% 'setTarget' {target, ...} -> new target; reply full 'data'
34% (density + samples).
35% 'movie' {dt, maxError, target} -> record orbit trajectories; 'movie'.
36d = ev.HTMLEventData;
37N = 1000; dt = 0.4; max_error = 0.8; target = 'banana';
38if isstruct(d)
39 if isfield(d, 'n'); N = max(1, round(d.n)); end
40 if isfield(d, 'dt'); dt = d.dt; end
41 if isfield(d, 'maxError'); max_error = d.maxError; end
42 if isfield(d, 'target'); target = d.target; end
43end
44set_target(target);
45switch ev.HTMLEventName
46 case 'resample'
47 samples = run_chain(N, dt, max_error, target);
48 sendEventToHTMLSource(src, 'samples', ...
49 struct('x', samples(1, :), 'y', samples(2, :), ...
50 'n', N, 'dt', dt, 'maxError', max_error, 'target', target));
51 case 'setTarget'
52 samples = run_chain(N, dt, max_error, target);
53 dens = density_grid(target);
54 sendEventToHTMLSource(src, 'data', ...
55 pack_data(samples, dens, N, dt, max_error, target));
56 case 'movie'
57 n_steps = 12;
58 [~, traj] = walnuts(target_start(target), n_steps, dt, max_error, true);
59 sendEventToHTMLSource(src, 'movie', traj);
60end
61end
63function set_target(name)
64% Select the target the density functions evaluate (a process-global so
65% walnuts.m's leapfrog can stay target-agnostic).
66global WTARGET
67WTARGET = target_code(name);
68end
70function c = target_code(name)
71switch name
72 case 'gaussian'
73 c = 2;
74 case 'correlated'
75 c = 3;
76 case 'donut'
77 c = 4;
78 otherwise
79 c = 1; % banana
80end
81end
83function s = target_start(name)
84% A sensible interior starting point for each target (the ring's hole is a bad
85% start, so the donut starts on the ring).
86switch name
87 case 'donut'
88 s = [2.5; 0];
89 case 'banana'
90 s = [0; 1];
91 otherwise
92 s = [0; 0];
93end
94end
96function [xmin, xmax, ymin, ymax] = target_bounds(name)
97switch name
98 case 'banana'
99 xmin = -6; xmax = 6; ymin = -7; ymax = 3;
100 otherwise
101 xmin = -4; xmax = 4; ymin = -4; ymax = 4;
102end
103end
105function samples = run_chain(N, dt, max_error, target)
106burnin = 200;
107chain = walnuts(target_start(target), burnin + N, dt, max_error);
108samples = chain(:, burnin + 1:end);
109end
111function dens = density_grid(target)
112%DENSITY_GRID Evaluate the selected target on a grid for the heatmap.
113% Row-major flat values: index (iy-1)*nx + ix, iy from ymin (1) to ymax (ny).
114[xmin, xmax, ymin, ymax] = target_bounds(target);
115nx = 100; ny = 100;
116xs = linspace(xmin, xmax, nx);
117ys = linspace(ymin, ymax, ny);
118values = zeros(1, nx * ny);
119k = 1;
120for iy = 1:ny
121 for ix = 1:nx
122 values(k) = log_density([xs(ix); ys(iy)]);
123 k = k + 1;
124 end
125end
126dens = struct('values', values, 'nx', nx, 'ny', ny, ...
127 'xmin', xmin, 'xmax', xmax, 'ymin', ymin, 'ymax', ymax);
128end
130function data = pack_data(samples, dens, N, dt, max_error, target)
131data = struct();
132data.type = 'walnuts';
133data.samples = struct('x', samples(1, :), 'y', samples(2, :));
134data.density = dens;
135data.n = N;
136data.dt = dt;
137data.maxError = max_error;
138data.target = target;
139end