function simulate_single_channel()

% SIMULATE_SINGLE_CHANNEL (UIFigure; whole window scrollable; controls at bottom)
% Radios: 0, 1, 2, 3, 20  (white background, black text; at bottom)
% Guarantees:
%   • ≥ 2 episodes with exactly n simultaneous openings (n = radio button)
%   • All levels 1..n occur at least once
%   • Strictly decreasing occupancy probabilities: P(1) > P(2) > ... > P(n)
% Filtered trace auto-plays; histograms render AFTER playback.

% -------------------- CONFIG --------------------
cfg.SIM_DURATION_MS            = 250;     % total simulated duration (ms)
cfg.VIEW_WINDOW_MS             = 200;     % visible range (ms)
cfg.FS_HZ                      = 20e3;    % sampling rate (Hz)

cfg.PLAYBACK_SPEED             = 0.03;    % slower playback
cfg.TICK_MS                    = 10;      % animation step (ms revealed per tick)

cfg.SINGLE_CHANNEL_OPEN_PROB   = 0.08;    % Po for ONE channel
cfg.MEAN_OPEN_TIME_MS          = 5.0;     % mean open time (ms)

cfg.UNITARY_CURRENT_PA         = -2.0;    % pA (negative=inward)
cfg.NOISE_SD_PA                = 0.25;    % pA (pre-filter Gaussian noise)
cfg.LPF_HZ                     = 2e3;     % one-pole low-pass cutoff (Hz)
cfg.BASELINE_PA                = 0.0;     % baseline current (pA)

cfg.COOP_BETA                  = 1.4;     % negative cooperativity factor

cfg.AUTO_YLIM                  = true;
cfg.YLIM_FIXED                 = [-6 1];

% ---- Guarantees / shaping ----
cfg.EXACT_N_EPISODES_REQUIRED  = 2;       % ≥ two exact-n episodes
cfg.EXACT_N_MS                 = 0.7;     % ms per injected exact-n episode

% Default injection pacing (auto-tightened for large n inside runSim)
cfg.INJECT_DUR_MS_RANGE        = [0.8 2.0];   % ms per injected k<n episode
cfg.INJECT_MIN_SPACING_MS      = 3.0;         % ms spacing between injections
cfg.MAX_INJECT_TRIES_PER_LEVEL = 180;

% Enforce P(k) >= (1+ratio_step)*P(k+1) (stronger near n)
cfg.MONO_NEAR_BAND_WIDTH       = 4;                 % n-1..n-4
cfg.MONO_STEP_RATIOS           = [0.60 0.40 0.25 0.15];
cfg.MONO_BASE_RATIO            = 0.08;
% ------------------------------------------------

cfg.T_SEC = cfg.SIM_DURATION_MS/1000;
cfg.N     = round(cfg.FS_HZ * cfg.T_SEC);

% ===== UIFigure =====
fig = uifigure('Color','w','Position',[100 100 1200 800], 'Name','Ion Channel Simulator');
fig.AutoResizeChildren = 'off';   % we control layout explicitly

% ===== ONE scrollable panel that fills the whole figure =====
contentPanel = uipanel(fig, ...
    'Scrollable','on', ...
    'BorderType','none', ...
    'Units','normalized', ...
    'Position',[0 0 1 1]);

% Inner "canvas" panel; we size/position it in layout_bottom()
canvas = uipanel(contentPanel,'Tag','scrollCanvas','BorderType','none','Units','pixels');

% Axes (pixel-managed)
axTrace = uiaxes(canvas); axTrace.BackgroundColor = 'w';
axAmp   = uiaxes(canvas); axAmp.BackgroundColor   = 'w';
axDwell = uiaxes(canvas); axDwell.BackgroundColor = 'w';

% Trace visuals
cla(axTrace);
hLine   = plot(axTrace, NaN, NaN, 'k','LineWidth',1.0); hold(axTrace,'on');
hCursor = plot(axTrace, NaN, NaN, 'r:','LineWidth',1.0);
ylabel(axTrace,'Current (pA)'); xlabel(axTrace,'Time (ms)');
title(axTrace,'Neuron Patch Recording'); grid(axTrace,'on');
axTrace.XLimMode = 'manual';

% ===== Controls at the BOTTOM (white, black text) =====
bg = uibuttongroup(canvas, ...
    'Title','Number of Channels', ...
    'FontSize',18, 'FontWeight','bold', ...
    'BackgroundColor','w');

rb0   = uiradiobutton(bg, 'Text','0',  'FontSize',18,'FontWeight','bold');
rb1   = uiradiobutton(bg, 'Text','1',  'FontSize',18,'FontWeight','bold','Value',true);
rb2   = uiradiobutton(bg, 'Text','2',  'FontSize',18,'FontWeight','bold');
rb3   = uiradiobutton(bg, 'Text','3',  'FontSize',18,'FontWeight','bold');
rb20  = uiradiobutton(bg, 'Text','20', 'FontSize',18,'FontWeight','bold');
rbList = [rb0 rb1 rb2 rb3 rb20];

bg.SelectionChangedFcn = @selectionChanged_cb;

% Shared state (store with guidata)
app.cfg          = cfg;
app.contentPanel = contentPanel;
app.canvas       = canvas;
app.axTrace      = axTrace;
app.axAmp        = axAmp;
app.axDwell      = axDwell;
app.hLine        = hLine;
app.hCursor      = hCursor;
app.bg           = bg;
app.rbList       = rbList;
app.timer        = [];
app.anim         = [];
guidata(fig, app);

% Initial layout + resize handler
layout_bottom(fig);
fig.SizeChangedFcn = @(~,~) onResize(fig);

% Start with 1 channel
runSim(fig, 1);

% Closing behavior
fig.CloseRequestFcn = @close_cb;
end  % ===== end main =====


% ========================== SUBFUNCTIONS ==========================

function onResize(fig)
% Single-source resize: recompute canvas size and lay out everything
layout_bottom(fig);
end

function selectionChanged_cb(src, event)
fig = ancestor(src,'figure');
n_ch = str2double(event.NewValue.Text);  % uifigure radios use "Text"
runSim(fig, n_ch);
end

function layout_bottom(fig)
% Size the scrollable content to fill the figure when large,
% and only be taller (thus scroll) when the figure is small.

app = guidata(fig);

% Viewport is the scrollable panel
vp = getpixelposition(app.contentPanel);      % [x y w h]
W  = max(800, vp(3));                         % min useful width
H  = max(500, vp(4));                         % min useful height

% ---- Minimum comfortable sizes (one histogram ROW) ----
margin   = 24;
gap      = 16;
groupH   = 100;
traceMin = 100;
histMin  = 160;

% Only ONE histogram row vertically:
minContentH = margin + traceMin + gap + histMin + gap + groupH + margin;

% If viewport is big enough, don't force scroll: canvas == viewport height.
% If small, make canvas >= minContentH to allow scrolling.
canvasH = max(H, minContentH);

set(app.canvas, 'Units','pixels', 'Position',[0 0 W canvasH]);

% Lay out children now that canvas size is set
layout_content_bottom(app.canvas, app.axTrace, app.axAmp, app.axDwell, app.bg);
end

function layout_content_bottom(canvas, axTrace, axAmp, axDwell, bg)
% Place plots to fill width; put the radio group at the BOTTOM.
% If canvas is taller than the minimum, expand plots (no scroll needed).

cp = getpixelposition(canvas);  % [x y w h]
W  = cp(3);
H  = cp(4);

% Keep these constants in sync with layout_bottom
margin   = 40;
gap      = 30;
groupH   = 80;
traceMin = 100;
histMin  = 140;

% Minimum total content height (one hist row)
minContentH = margin + traceMin + gap + histMin + gap + groupH + margin;

% Space available for the two vertical blocks: [trace] and [hist row]
availForPlots = H - (2*margin + 2*gap + groupH);

% Make the trace roughly half as tall as before in spacious layouts
traceFrac = 0.30;              % was effectively 0.60 before
rowFrac   = 1 - traceFrac;     % 0.70

if H >= minContentH
    traceH = max(traceMin, round(traceFrac * availForPlots));
    histH  = max(histMin,  round(rowFrac   * availForPlots));
else
    traceH = traceMin;
    histH  = histMin;
end


% --- Bottom controls (span width) ---
bg.Position = [margin, margin, W - 2*margin, groupH];

% Equal-width radio buttons in numeric order (0,1,2,3,20)
kids = bg.Children;
vals = arrayfun(@(c) str2double(c.Text), kids);
[~, ix] = sort(vals, 'ascend');  kids = kids(ix);

n = numel(kids);
gapBtn = 14;
btnW = max(90, floor((bg.Position(3) - (n+1)*gapBtn)/n));
btnH = 34;  yBtn = (groupH - btnH)/2;
for i = 1:n
    x = gapBtn*i + btnW*(i-1);
    kids(i).Position = [x, yBtn, btnW, btnH];
end

% --- Plots (fill width above the controls; no big top blank) ---
% Histogram ROW sits directly above the group
yHist = margin + groupH + gap;
w2    = round((W - 3*margin)/2);
set(axAmp,  'Units','pixels', 'Position',[margin,         yHist, w2, histH]);
set(axDwell,'Units','pixels', 'Position',[2*margin + w2,  yHist, w2, histH]);

% Trace sits directly above the histogram row; top margin only
yTrace = yHist + histH + gap;
set(axTrace,'Units','pixels', 'Position',[margin, yTrace, W - 2*margin, traceH]);
end

function runSim(fig, n_ch)
app = guidata(fig);
cfg = app.cfg;
rng('shuffle');

% Kinetics from Po and mean open time
tau_open_s = cfg.MEAN_OPEN_TIME_MS/1000;
k_oc = 1/tau_open_s;                                   % O->C
Po   = max(0,min(0.999,cfg.SINGLE_CHANNEL_OPEN_PROB));
k_co = k_oc * Po/(1-Po);                               % C->O

% Per-sample transition probabilities
dt        = 1/cfg.FS_HZ;
t         = (0:cfg.N-1)*dt;
p_co_base = 1 - exp(-k_co*dt);
p_oc      = 1 - exp(-k_oc*dt);

% Dynamic pacing for large n (to fit everything into 250 ms)
spacing_eff   = cfg.INJECT_MIN_SPACING_MS;
dur_rng_eff   = cfg.INJECT_DUR_MS_RANGE;
if n_ch >= 100
    spacing_eff = 0.40;  dur_rng_eff = [0.15 0.40];
elseif n_ch >= 50
    spacing_eff = 0.60;  dur_rng_eff = [0.20 0.50];
elseif n_ch >= 20
    spacing_eff = 1.20;  dur_rng_eff = [0.40 1.00];
end

% ----- Simulate states with negative cooperativity -----
state = false(n_ch, cfg.N);
for k = 2:cfg.N
    if n_ch>0
        s_prev = state(:,k-1);
        m = sum(s_prev);
        p_co_eff = p_co_base * exp(-cfg.COOP_BETA * m);  % suppress high occupancy
        s_new = s_prev;
        if any(s_prev), s_new(s_prev) = rand(nnz(s_prev),1) >= p_oc; end % O->C
        if any(~s_prev), s_new(~s_prev) = rand(nnz(~s_prev),1) < p_co_eff; end % C->O
        state(:,k) = s_new;
    end
end

% ----- Guarantees: (a) ≥2 exact-n; (b) all k=1..n once -----
if n_ch > 0
    used = false(1,cfg.N);

    % (a) ≥2 exact-n episodes
    target_eq = sum(state,1) == n_ch;
    [s_eq, ~] = contiguous_runs(target_eq);
    need_eq = max(0, cfg.EXACT_N_EPISODES_REQUIRED - numel(s_eq));
    if need_eq > 0
        anchors = linspace(0.15, 0.40, need_eq+1); anchors(end) = [];
        for ii = 1:need_eq
            [state, used] = inject_k_sub(state, used, n_ch, cfg.FS_HZ, ...
                              n_ch, cfg.EXACT_N_MS, spacing_eff, anchors(ii));
        end
    end

    % (b) at least one for each k < n
    for k_occ = 1:max(0,n_ch-1)
        if ~any(sum(state,1) == k_occ)
            len_ms = mean(dur_rng_eff);
            [state, used] = inject_k_sub(state, used, n_ch, cfg.FS_HZ, ...
                              k_occ, len_ms, spacing_eff, []);
        end
    end
end

% ----- Enforce strictly decreasing occupancy probabilities -----
if n_ch > 1
    tries_per_k = cfg.MAX_INJECT_TRIES_PER_LEVEL;
    while true
        oc = sum(state,1);
        pk = zeros(1,n_ch);
        for k_occ = 1:n_ch
            pk(k_occ) = mean(oc == k_occ);
        end
        fixed = true;
        for k_occ = n_ch-1:-1:1
            band_idx = n_ch - k_occ; % 1 for n-1, 2 for n-2, ...
            if band_idx >=1 && band_idx <= cfg.MONO_NEAR_BAND_WIDTH
                ratio = cfg.MONO_STEP_RATIOS(min(band_idx, numel(cfg.MONO_STEP_RATIOS)));
            else
                ratio = cfg.MONO_BASE_RATIO;
            end
            if pk(k_occ) < (1+ratio) * pk(k_occ+1)
                fixed = false;
                ntries = tries_per_k;
                while ntries > 0 && pk(k_occ) < (1+ratio)*pk(k_occ+1)
                    len_ms = dur_rng_eff(1) + diff(dur_rng_eff)*rand;
                    [state, ~] = inject_k_sub(state, false(1,cfg.N), n_ch, cfg.FS_HZ, ...
                                        k_occ, len_ms, spacing_eff, []);
                    oc = sum(state,1);
                    pk(k_occ)   = mean(oc == k_occ);
                    pk(k_occ+1) = mean(oc == (k_occ+1));
                    ntries = ntries - 1;
                end
            end
        end
        if fixed || all(pk(1:end-1) > pk(2:end)), break; end
    end
end

% ----- Current & filtering -----
open_count = sum(state,1);
i_ideal    = cfg.BASELINE_PA + double(open_count)*cfg.UNITARY_CURRENT_PA;
i_raw      = i_ideal + cfg.NOISE_SD_PA*randn(1,cfg.N);

alpha  = exp(-2*pi*cfg.LPF_HZ*dt);
b0     = 1 - alpha; a1 = alpha;
i_filt = filter(b0, [1 -a1], i_raw);

% ----- Store animation state -----
anim.Fs       = cfg.FS_HZ;
anim.t_ms     = t*1e3;
anim.y        = i_filt;
anim.idx      = 1;
anim.step     = max(1, round(cfg.FS_HZ*(cfg.TICK_MS/1000)));
anim.speed    = cfg.PLAYBACK_SPEED;
anim.open_cnt = open_count;
anim.n_ch     = n_ch;
anim.view_ms  = min(cfg.VIEW_WINDOW_MS, cfg.SIM_DURATION_MS);
if n_ch==0
    anim.dwell_ch1 = [];
else
    anim.dwell_ch1 = dwell_times_from_state(state(1,:), 1, dt); % seconds
end
app.anim = anim;
guidata(fig, app);

% ----- Reset trace display & axes -----
set(app.hLine,'XData',NaN,'YData',NaN);
set(app.hCursor,'XData',NaN,'YData',NaN);
title(app.axTrace,'Neuron Patch Recording'); grid(app.axTrace,'on');
xlim(app.axTrace,[0 anim.view_ms]);

% Y limits
if cfg.AUTO_YLIM
    mask = anim.t_ms >= 0 & anim.t_ms <= anim.view_ms;
    if any(mask)
        yl = [min(anim.y(mask)), max(anim.y(mask))];
        if diff(yl)==0, yl = yl + [-1 1]; end
        pad = 0.1*range(yl);
        ylim(app.axTrace, [yl(1)-pad, yl(2)+pad]);
    else
        ylim(app.axTrace, cfg.YLIM_FIXED);
    end
else
    ylim(app.axTrace, cfg.YLIM_FIXED);
end

% Prepare hist axes ("Playing…")
cla(app.axAmp);  title(app.axAmp,'Open Channel Amplitude Distribution');
xlabel(app.axAmp,'Current (pA)'); ylabel(app.axAmp,'Counts'); grid(app.axAmp,'on');
text(app.axAmp,0.5,0.5,'Playing…','HorizontalAlignment','center','Units','normalized');

cla(app.axDwell); title(app.axDwell,'Open Channel Time Distribution');
xlabel(app.axDwell,'Open time (ms)'); ylabel(app.axDwell,'Counts'); grid(app.axDwell,'on');
text(app.axDwell,0.5,0.5,'Playing…','HorizontalAlignment','center','Units','normalized');

% Start animation
startTimer(fig);

% Reposition to monitor and re-layout (keeps “controls at bottom”)
% try
%     monPos = get(0, 'MonitorPositions');
%     if size(monPos,1) >= 2, mon = monPos(2,:); else, mon = monPos(1,:); end
%     figW = 1200; figH = 800;  x = mon(1) + 100;  y = mon(2) + mon(4) - figH - 100;
%     fig.Position = [x y figW figH];
% catch
% end
% Only do initial positioning once
app = guidata(fig);
if ~isfield(app,'didInitPos') || ~app.didInitPos
    try
        monPos = get(0, 'MonitorPositions');
        if size(monPos,1) >= 2, mon = monPos(2,:); else, mon = monPos(1,:); end
        figW = 1200; figH = 800;  x = mon(1) + 100;  y = mon(2) + mon(4) - figH - 100;
        fig.Position = [x y figW figH];
    catch
    end
    app.didInitPos = true;
    guidata(fig, app);
end

% Still re-layout content (safe and needed)
layout_bottom(fig);

layout_bottom(fig);
end

function startTimer(fig)
app = guidata(fig);
stopTimer(fig);
period = max( (app.cfg.TICK_MS/1000) / app.anim.speed, 0.01 );
t = timer('ExecutionMode','fixedSpacing', 'Period', period, 'TimerFcn', @onTick);
t.UserData = fig;
app.timer = t;
guidata(fig, app);
start(t);
end

function stopTimer(fig)
app = guidata(fig);
if ~isempty(app) && isfield(app,'timer') && ~isempty(app.timer) && isvalid(app.timer)
    try, stop(app.timer); catch, end
    try, delete(app.timer); catch, end
end
app.timer = [];
if ~isempty(app), guidata(fig, app); end
end

function onTick(t, ~)
fig = t.UserData;
if ~isvalid(fig), try, stop(t); delete(t); end, return; end
app = guidata(fig);
if isempty(app) || ~isfield(app,'anim') || isempty(app.anim)
    try, stop(t); delete(t); end
    return;
end
anim = app.anim;

anim.idx = min(anim.idx + anim.step, numel(anim.y));
set(app.hLine,'XData',anim.t_ms(1:anim.idx), 'YData',anim.y(1:anim.idx));
tcur = anim.t_ms(anim.idx);
yl = ylim(app.axTrace);
set(app.hCursor,'XData',[tcur tcur], 'YData',yl);

if anim.idx >= numel(anim.y)
    app.anim = anim; guidata(fig, app);
    stopTimer(fig);
    showHistograms(fig);
else
    app.anim = anim; guidata(fig, app);
end
drawnow limitrate;
end

function showHistograms(fig)
app = guidata(fig);
anim = app.anim;
n_ch = anim.n_ch;

cla(app.axAmp);  title(app.axAmp,'Open Channel Amplitude Distribution');
xlabel(app.axAmp,'Current (pA)'); ylabel(app.axAmp,'Counts'); grid(app.axAmp,'on');
if n_ch==0
    text(app.axAmp,0.5,0.5,'No open samples (0 channels)', ...
        'HorizontalAlignment','center','Units','normalized');
else
    open_mask = anim.open_cnt > 0;
    open_samples = anim.y(open_mask);
    if isempty(open_samples)
        text(app.axAmp,0.5,0.5,'No detected open samples', ...
            'HorizontalAlignment','center','Units','normalized');
    else
        histogram(app.axAmp, open_samples, 100);  % COUNTS
    end
end

cla(app.axDwell); title(app.axDwell,'Open Channel Time Distribution');
xlabel(app.axDwell,'Open time (ms)'); ylabel(app.axDwell,'Counts'); grid(app.axDwell,'on');
dwell_open = anim.dwell_ch1;  % seconds
if isempty(dwell_open)
    text(app.axDwell,0.5,0.5,'No openings detected', ...
        'HorizontalAlignment','center','Units','normalized');
else
    histogram(app.axDwell, 1e3*dwell_open, 40);
end
end

function close_cb(fig, ~)
stopTimer(fig);
if isvalid(fig), delete(fig); end
end

% ---------------------- helpers ----------------------
function [state, used] = inject_k_sub(state, used, n_ch, Fs, k_occ, len_ms, spacing_ms, anchor_frac)
% Inject a short window with exactly k_occ simultaneous openings
N = size(state,2);
L = max(1, round(Fs*(len_ms/1000)));
if nargin >= 8 && ~isempty(anchor_frac)
    pos = clamp(round(anchor_frac*N), 1, max(1,N-L+1));
else
    pos = randi([1, max(1,N-L+1)],1,1);
end
sep = max(1, round(Fs*(spacing_ms/1000)));
tries = 200;
while tries > 0
    lo = max(1, pos - sep);
    hi = min(N, pos + L - 1 + sep);
    if isempty(used), used = false(1,N); end
    if ~any(used(lo:hi)), break; end
    pos = randi([1, max(1,N-L+1)],1,1);
    tries = tries - 1;
end
lo = max(1, pos - sep);
hi = min(N, pos + L - 1 + sep);
which = randperm(n_ch, max(1,min(k_occ,n_ch)));
state(:, pos:pos+L-1) = false;
state(which, pos:pos+L-1) = true;
used(lo:hi) = true;
end

function [starts, ends_] = contiguous_runs(mask)
mask = mask(:)~=0;
d = diff([false; mask; false]);
starts = find(d==1);
ends_  = find(d==-1)-1;
end

function dts = dwell_times_from_state(state_vec, target, dt)
sv = logical(state_vec(:));
if target==0, sv = ~sv; end
d = diff([false; sv; false]);
run_starts = find(d==1);
run_ends   = find(d==-1)-1;
dts = (run_ends - run_starts + 1) * dt;
end

function v = clamp(v, a, b)
v = min(max(v, a), b);
end