% detect_binding_bleach.m
% 用法（在终端）:
%   octave detect_binding_bleach.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv

%--------------------------------------------
% 基本设置
%--------------------------------------------
arg_list = argv();
if numel(arg_list) < 1
    error("Usage: octave detect_binding_bleach.m 
<track_intensity_csv>");
end
input_file = arg_list{1};

% 加载 HMM 采样器
%addpath("sampler_SRC");

% 如果安装了 statistics 包，用于 kmeans
try
    pkg load statistics;
catch
    warning("Could not load 'statistics' package. Make sure it is installed if kmeans is missing.");
end

%--------------------------------------------
% 读入 track intensity 文件（分号分隔）
%--------------------------------------------
fprintf("Reading file: %s\n", input_file);

fid = fopen(input_file, "r");
if fid < 0
    error("Cannot open file: %s", input_file);
end

% 第一行是注释头 "# track index;time (seconds);track intensity (photon counts)"
header_line = fgetl(fid);  % 忽略内容，只是读掉这一行

% 后面每行: track_index;time_sec;intensity
data = textscan(fid, "%f%f%f", "Delimiter", ";");
fclose(fid);

track_idx = data{1};
time_sec  = data{2};
counts    = data{3};

% 排序（确保同一个 track 内按时间排序）
[~, order] = sortrows([track_idx, time_sec], [1, 2]);
track_idx = track_idx(order);
time_sec  = time_sec(order);
counts    = counts(order);

tracks = unique(track_idx);
n_tracks = numel(tracks);

fprintf("Found %d tracks.\n", n_tracks);

%--------------------------------------------
% 结果结构体
%--------------------------------------------
results = struct( ...
    "track_id", {}, ...
    "binding_indices", {}, ...
    "binding_times", {}, ...
    "bleach_indices", {}, ...
    "bleach_times", {} );

%--------------------------------------------
% 循环每个 track，跑 HMM + 找事件
%--------------------------------------------
for ti = 1:n_tracks
    tr = tracks(ti);
    fprintf("\n===== Track %d =====\n", tr);

    sel = (track_idx == tr);
    t = time_sec(sel);
    z = counts(sel);
    z = z(:);      % 列向量

    %-------------------------
    % 1) 设定 ICON HMM 的参数
    %-------------------------
    opts = struct();

    % 超参数（和原 MATLAB 代码一致）
    opts.a = 1;
    opts.g = 1;

    opts.Q(1) = mean(z);
    opts.Q(2) = 1 / (std(z)^2 + eps);   % 防止除零
    opts.Q(3) = 0.1;
    opts.Q(4) = 0.00001;

    opts.M = 10;

    % 采样参数
    opts.dr_sk  = 1;
    opts.K_init = 50;

    flag_stat = true;
    flag_anim = false;    % 建议在 Octave 里关掉动画更稳定

    R = 1000;             % 采样次数，可根据需要调小/调大

    %-------------------------
    % 2) 运行采样器
    %-------------------------
    fprintf("  Running HMM sampler...\n");
    chain = chainer_main(z, R, [], opts, flag_stat, flag_anim);

    %-------------------------
    % 3) 做后验分析，得到平滑后的光子水平轨迹
    %-------------------------
    fr = 0.25;        % burn-in 比例
    dr = 2;           % sample 步距

    m_min = 0;
    m_max = 1;
    m_num = 25;

    [m_mod, m_red] = chainer_analyze_means(chain, fr, dr, m_min, m_max, m_num, z);

    % m_mod 是一个和 z 同长度的向量（通常是 [0,1] 范围的归一化水平）
    m_mod = m_mod(:);

    %-------------------------
    % 4) 把 m_mod 聚类成 3 个光子强度状态
    %    状态 1 = 最低强度 (dark / unbound / bleached)
    %    状态 3 = 最高强度 (bound)
    %-------------------------
    K = 3;   % 你可以改成 2 或其他
    try
        [idx_raw, centers] = kmeans(m_mod, K);
    catch
        % 如果没 kmeans，就简单按数值分位数来分三类
        warning("kmeans not available, using simple quantile-based clustering.");
        q1 = quantile(m_mod, 1/3);
        q2 = quantile(m_mod, 2/3);
        idx_raw = ones(size(m_mod));
        idx_raw(m_mod > q1 & m_mod <= q2) = 2;
        idx_raw(m_mod > q2) = 3;
        centers = zeros(K,1);
        for kk = 1:K
            centers(kk) = mean(m_mod(idx_raw == kk));
        end
    end

    % 根据中心值从小到大重排状态编号
    [~, order_centers] = sort(centers);
    state_seq = zeros(size(idx_raw));
    for k = 1:K
        state_seq(idx_raw == order_centers(k)) = k;
    end

    low_state  = 1;
    high_state = K;

    %-------------------------
    % 5) 检测 binding / bleach 跳变
    %    binding:  low -> high
    %    bleach : high -> low
    %-------------------------
    s = state_seq(:);
    % 前一时刻和后一时刻
    s_prev = s(1:end-1);
    s_next = s(2:end);

    bind_idx   = find(s_prev == low_state  & s_next == high_state) + 1;
    bleach_idx = find(s_prev == high_state & s_next == low_state) + 1;

    bind_times   = t(bind_idx);
    bleach_times = t(bleach_idx);

    fprintf("  Found %d binding event(s) and %d bleaching event(s).\n", ...
            numel(bind_idx), numel(bleach_idx));

    % 存入结果
    results(ti).track_id        = tr;
    results(ti).binding_indices = bind_idx;
    results(ti).binding_times   = bind_times;
    results(ti).bleach_indices  = bleach_idx;
    results(ti).bleach_times    = bleach_times;
end

%--------------------------------------------
% 6) 把结果写成 CSV
%--------------------------------------------
out_csv = "binding_bleach_events.csv";
fid_out = fopen(out_csv, "w");
if fid_out < 0
    error("Cannot open output file: %s", out_csv);
end

fprintf(fid_out, "track_index,event_type,sample_index,time_seconds\n");
for ti = 1:numel(results)
    tr = results(ti).track_id;

    % binding
    for k = 1:numel(results(ti).binding_indices)
        fprintf(fid_out, "%d,binding,%d,%.6f\n", tr, ...
                results(ti).binding_indices(k), ...
                results(ti).binding_times(k));
    end

    % bleach
    for k = 1:numel(results(ti).bleach_indices)
        fprintf(fid_out, "%d,bleach,%d,%.6f\n", tr, ...
                results(ti).bleach_indices(k), ...
                results(ti).bleach_times(k));
    end
end
fclose(fid_out);

fprintf("\nAll done. Events written to: %s\n", out_csv);

• 运行：

octave detect_binding_bleach.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv

---

1. 保留原始 ICON 后验分析：
   • chainer_export 导出样本
   • chainer_analyze_means 得到 m_mod
   • chainer_analyze_transitions 得到 m_edges, p_mean, d_dist
   • chainer_analyze_drift 得到 y_mean, y_std
2. 使用 binding/bleach 时间点 计算 dwell time 分布并做单指数 rate 拟合，并直接画直方图 + 指数拟合曲线。
3. 输出一个 CSV 表，包含所有 binding/bleach 时间点（这是最重要的输出）。

说明：

• 代码仍然基于你已有的 detect_binding_bleach.m 结构，假定当前目录下有 sampler_SRC（含 chainer_main.m、chainer_analyze_means.m、chainer_analyze_transitions.m、chainer_analyze_drift.m 等）。¹²
• Dwell time 这里使用“连续 binding 事件之间的时间差”来估计 bound-state dwell time（简化版）；如果 bleaching 明确是轨迹末端终止，也可以用 bleach_time - last_binding_time 作为一个 right‑censored dwell（这里暂不做 censoring 修正，只输出原始分布和单指数拟合）。

---

完整脚本：detect_binding_bleach_dwell_simple.m

% detect_binding_bleach_dwell.m
%
% 用法（终端）:
%   octave detect_binding_bleach_dwell.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv
%
% 输入 CSV 格式（分号分隔）:
%   # track index;time (seconds);track intensity (photon counts)
%   1;0.00;123
%   1;0.05;118
%   2;0.00; 95
%   ...

%--------------------------------------------
% 0. 基本设置
%--------------------------------------------
arg_list = argv();

if numel(arg_list) < 1
    error("Usage: octave detect_binding_bleach_dwell.m 
<input_csv>");
end

input_file = arg_list{1};

% 加载 HMM 采样器代码（确保 sampler_SRC 在当前目录下）
addpath("sampler_SRC");

% 如果安装了 statistics 包，用于 kmeans
try
    pkg load statistics;
catch
    warning("Could not load 'statistics' package. Make sure it is installed if kmeans is missing.");
end

%--------------------------------------------
% 1. 读入 track intensity 文件（分号分隔）
%--------------------------------------------
fprintf("Reading file: %s\n", input_file);

fid = fopen(input_file, "r");
if fid < 0
    error("Cannot open file: %s", input_file);
end

% 第一行是注释头 "# track index;time (seconds);track intensity (photon counts)"
header_line = fgetl(fid); % 忽略内容，只是读掉这一行

% 后面每行: track_index;time_sec;intensity
data = textscan(fid, "%f%f%f", "Delimiter", ";");
fclose(fid);

track_idx = data{1};
time_sec  = data{2};
counts    = data{3};

% 排序（确保同一个 track 内按时间排序）
[~, order] = sortrows([track_idx, time_sec], [1, 2]);
track_idx = track_idx(order);
time_sec  = time_sec(order);
counts    = counts(order);

tracks   = unique(track_idx);
n_tracks = numel(tracks);
fprintf("Found %d tracks.\n", n_tracks);

%--------------------------------------------
% 2. 结果结构体（binding / bleach 时间点）
%--------------------------------------------
results = struct( ...
    "track_id",        {}, ...
    "binding_indices", {}, ...
    "binding_times",   {}, ...
    "bleach_indices",  {}, ...
    "bleach_times",    {} );

% 用于 dwell time 收集的数组
all_binding_times = [];   % 所有 track 的 binding time 合并
all_bleach_times  = [];   % 所有 track 的 bleach time 合并

%--------------------------------------------
% 3. 循环每个 track，跑 HMM + ICON 分析 + 事件检测
%--------------------------------------------
for ti = 1:n_tracks

    tr = tracks(ti);
    fprintf("\n===== Track %d =====\n", tr);

    sel = (track_idx == tr);
    t   = time_sec(sel);
    z   = counts(sel);
    z   = z(:); % 列向量

    %-------------------------
    % 3.1 设定 ICON HMM 的参数
    %-------------------------
    opts        = struct();
    % 超参数（与 level_analysis.m 一致）
    opts.a      = 1;
    opts.g      = 1;
    opts.Q(1)   = mean(z);
    opts.Q(2)   = 1 / (std(z)^2 + eps); % 防止 std(z)=0 除零
    opts.Q(3)   = 0.1;
    opts.Q(4)   = 0.00001;
    opts.M      = 10;

    % 采样参数
    opts.dr_sk  = 1;
    opts.K_init = 50;

    flag_stat   = true;
    flag_anim   = false;  % 建议在 Octave 里关掉动画
    R           = 1000;   % 采样次数，可按需求调整

    %-------------------------
    % 3.2 运行采样器（ICON HMM）
    %-------------------------
    fprintf(" Running HMM sampler...\n");
    chain = chainer_main(z, R, [], opts, flag_stat, flag_anim);

    %-------------------------
    % 3.3 ICON 后验分析（与 level_analysis.m 等价）
    %-------------------------
    fr    = 0.25; % burn-in 比例
    dr    = 2;    % sample 步距
    m_min = 0;
    m_max = 1;
    m_num = 25;

    % (1) 均值轨迹：m_mod
    [m_mod, m_red] = chainer_analyze_means(chain, fr, dr, m_min, m_max, m_num, z);
    m_mod = m_mod(:);

    % (2) 导出 samples（为每条 track 单独存一个文件）
    sample_file = sprintf("samples_track_%d", tr);
    chainer_export(chain, fr, dr, sample_file, "mat");

    % (3) 转移概率 / 跃迁统计
    [m_edges, p_mean, d_dist] = chainer_analyze_transitions( ...
        chain, fr, dr, m_min, m_max, m_num, true);

    % (4) 漂移分析
    [y_mean, y_std] = chainer_analyze_drift(chain, fr, dr, z);

    % 你可以根据需要保存这些 ICON 分析结果
    mat_out = sprintf("icon_analysis_track_%d.mat", tr);
    save(mat_out, "m_mod", "m_red", "m_edges", "p_mean", "d_dist", ...
                  "y_mean", "y_std", "t", "z");

    %-------------------------
    % 3.4 把 m_mod 聚类成 3 个光子强度状态（用于事件检测）
    %-------------------------
    K = 3; % 可按物理需求修改

    try
        [idx_raw, centers] = kmeans(m_mod, K);
    catch
        warning("kmeans not available, using simple quantile-based clustering.");
        q1 = quantile(m_mod, 1/3);
        q2 = quantile(m_mod, 2/3);

        idx_raw = ones(size(m_mod));
        idx_raw(m_mod > q1 & m_mod <= q2) = 2;
        idx_raw(m_mod > q2)              = 3;

        centers = zeros(K,1);
        for kk = 1:K
            centers(kk) = mean(m_mod(idx_raw == kk));
        end
    end

    % 根据中心值从小到大重排状态编号，使 state_seq = 1..K 对应 low->high
    [~, order_centers] = sort(centers);
    state_seq = zeros(size(idx_raw));
    for k = 1:K
        state_seq(idx_raw == order_centers(k)) = k;
    end

    low_state  = 1;
    high_state = K;

    %-------------------------
    % 3.5 检测 binding / bleaching
    %
    %   binding: low_state -> high_state
    %   bleach : high_state -> low_state
    %-------------------------
    s      = state_seq(:);
    s_prev = s(1:end-1);
    s_next = s(2:end);

    bind_idx   = find(s_prev == low_state  & s_next == high_state) + 1;
    bleach_idx = find(s_prev == high_state & s_next == low_state) + 1;

    bind_times   = t(bind_idx);
    bleach_times = t(bleach_idx);

    fprintf(" Found %d binding event(s) and %d bleaching event(s).\n", ...
            numel(bind_idx), numel(bleach_idx));

    % 存入结果结构体
    results(ti).track_id        = tr;
    results(ti).binding_indices = bind_idx;
    results(ti).binding_times   = bind_times;
    results(ti).bleach_indices  = bleach_idx;
    results(ti).bleach_times    = bleach_times;

    % 用于全局 dwell-time 统计
    all_binding_times = [all_binding_times; bind_times(:)];
    all_bleach_times  = [all_bleach_times;  bleach_times(:)];

end

%--------------------------------------------
% 4. 输出 binding/bleach 时间点的总表（最重要输出）
%--------------------------------------------
out_csv = "binding_bleach_events.csv";
fid_out = fopen(out_csv, "w");
if fid_out < 0
    error("Cannot open output file: %s", out_csv);
end

fprintf(fid_out, "track_index,event_type,sample_index,time_seconds\n");

for ti = 1:numel(results)
    tr = results(ti).track_id;

    % binding
    for k = 1:numel(results(ti).binding_indices)
        fprintf(fid_out, "%d,binding,%d,%.6f\n", tr, ...
                results(ti).binding_indices(k), ...
                results(ti).binding_times(k));
    end

    % bleach
    for k = 1:numel(results(ti).bleach_indices)
        fprintf(fid_out, "%d,bleach,%d,%.6f\n", tr, ...
                results(ti).bleach_indices(k), ...
                results(ti).bleach_times(k));
    end
end

fclose(fid_out);
fprintf("\n[INFO] Binding/bleach events written to: %s\n", out_csv);

%--------------------------------------------
% 5. Dwell-time 直方图 + rate 拟合
%
%   简化假设：
%   - 把每个 track 中连续 binding 事件之间的时间差视为
%     bound state dwell time（中间没有 bleaching，就重新 binding）。
%   - 如果有明确的最后一个 binding 到 bleaching 的时间，
%     你也可以扩展这里把它们当作 dwell time 加进去。
%--------------------------------------------

% 5.1 收集 bound-state dwell times（简单用 binding 时间差）
dwell_times = [];

for ti = 1:numel(results)
    bt = sort(results(ti).binding_times(:)); % 单个 track 的 binding times
    if numel(bt) >= 2
        dt = diff(bt);    % 相邻 binding 之间的时间差
        dwell_times = [dwell_times; dt];
    end

    % 可选：如果你想把最后一个 binding -> 第一个 bleach 也算入 dwell:
    % if ~isempty(results(ti).binding_times) && ~isempty(results(ti).bleach_times)
    %     last_binding  = max(results(ti).binding_times);
    %     first_bleach  = min(results(ti).bleach_times);
    %     if first_bleach > last_binding
    %         dwell_times = [dwell_times; first_bleach - last_binding];
    %     end
    % end
end

if isempty(dwell_times)
    fprintf("[WARN] No sufficient binding events to compute dwell-time distribution.\n");
else
    % 5.2 绘制 dwell-time 直方图
    figure;
    hold on;
    nbins = max(10, round(sqrt(numel(dwell_times)))); % 简单经验 bin 数
    [counts_hist, edges] = histcounts(dwell_times, nbins);
    centers = (edges(1:end-1) + edges(2:end)) / 2;

    bar(centers, counts_hist, "hist");
    xlabel("Dwell time (s)");
    ylabel("Counts");
    title("Bound-state dwell-time histogram");

    % 5.3 拟合单指数: p(t) ~ (1/tau) * exp(-t/tau)
    % 使用简单的最小二乘拟合 log(count) 对 t
    valid = counts_hist > 0;
    t_fit = centers(valid);
    y_fit = log(counts_hist(valid));

    % 线性拟合: log(N) = a + b * t, 其中 b ~= -1/tau
    p = polyfit(t_fit, y_fit, 1);
    a = p(1);
    b = p(2);

    tau_est = -1 / b;
    fprintf("[INFO] Fitted single-exponential dwell time tau = %.4f s\n", tau_est);

    % 画出拟合曲线（归一到 histogram 的总 counts 约束）
    t_plot = linspace(min(dwell_times), max(dwell_times), 200);
    % 预测密度 ~ exp(a + b t)，但我们要匹配 histogram 的离散计数
    y_model = exp(a + b * t_plot);

    plot(t_plot, y_model, "r-", "LineWidth", 2);
    legend("Histogram", "Single-exp fit");
    hold off;

    % 可选：保存图像
    print("-dpng", "dwell_time_histogram_fit.png");
    fprintf("[INFO] Dwell-time histogram and fit saved to dwell_time_histogram_fit.png\n");
end

fprintf("\nAll done.\n");

---

使用与扩展建议

• 运行：

octave detect_binding_bleach_dwell_simple.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv

会生成：

• binding_bleach_events.csv：所有 track 的 binding/bleach 时间点（你最关心的表格）。
• icon_analysis_track_X.mat：每个 track 的 ICON 分析结果（m_mod, m_edges, p_mean, d_dist, y_mean, y_std 等）。
• 如果有足够事件，生成 dwell_time_histogram_fit.png，包含 dwell-time 直方图和单指数拟合曲线。
  • dwell time 定义： 当前用的是“相邻 binding 时间点之差”作为 dwell；你如果有更明确的物理模型（例如每条轨迹只有 1 次 binding + 1 次 bleach），可以把注释里那段 last_binding -> first_bleach 的逻辑打开或进一步改写成更严格的两态 dwell 定义。
  • rate constant： 对于单指数 $p(t) \propto \exp(-t/\tau)$，tau_est 是平均 dwell time，解离 rate 大致是 $k{\text{off}} = 1/\tau{\text{est}}$。如果你需要可以在脚本里再打印一条 k_off_est。

如果你希望 dwell-time 统计严格区分 “bound state dwell vs unbound state dwell”，或者想直接从 p_mean / d_dist 中做 rate 矩阵估计，可以在这份脚本基础上再往前一步；可以一起设计一个更“教科书版”的 HMM‑rate 推断流程。

---

特性概览：

• 使用 ICON HMM（chainer_main）分析每个 track。¹
• 做原始 level_analysis.m 的四类后验分析：chainer_export、chainer_analyze_means、chainer_analyze_transitions、chainer_analyze_drift。¹
• 基于 m_mod 聚类成 3 个强度状态，检测 binding / bleaching 事件。²
• 输出 binding_bleach_events.csv，包含所有 binding/bleach 时间点。²
• 用 binding 时间点计算 dwell-time 分布，做单指数拟合；绘图部分只在内存中生成，用 print + close 保存 PNG（类似 R 的 dev.off()），并使用 hist 而不是 histcounts。²

将以下内容保存为 detect_binding_bleach_dwell.m 即可直接运行。

---

detect_binding_bleach_dwell.m 完整代码

% detect_binding_bleach_dwell.m
%
% 用法（终端）:
%   octave detect_binding_bleach_dwell.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv
%
% 输入 CSV 格式（分号分隔）:
%   # track index;time (seconds);track intensity (photon counts)
%   1;0.00;123
%   1;0.05;118
%   2;0.00; 95
%   ...

%--------------------------------------------
% 0. 基本设置
%--------------------------------------------
arg_list = argv();

if numel(arg_list) < 1
    error("Usage: octave detect_binding_bleach_dwell.m 
<input_csv>");
end

input_file = arg_list{1};

% 加载 HMM 采样器代码（确保 sampler_SRC 在当前目录下）
addpath("sampler_SRC");

% 如果安装了 statistics 包，用于 kmeans
try
    pkg load statistics;
catch
    warning("Could not load 'statistics' package. Make sure it is installed if kmeans is missing.");
end

%--------------------------------------------
% 1. 读入 track intensity 文件（分号分隔）
%--------------------------------------------
fprintf("Reading file: %s\n", input_file);

fid = fopen(input_file, "r");
if fid < 0
    error("Cannot open file: %s", input_file);
end

% 第一行是注释头 "# track index;time (seconds);track intensity (photon counts)"
header_line = fgetl(fid); % 忽略内容，只是读掉这一行

% 后面每行: track_index;time_sec;intensity
data = textscan(fid, "%f%f%f", "Delimiter", ";");
fclose(fid);

track_idx = data{1};
time_sec  = data{2};
counts    = data{3};

% 排序（确保同一个 track 内按时间排序）
[~, order] = sortrows([track_idx, time_sec], [1, 2]);
track_idx = track_idx(order);
time_sec  = time_sec(order);
counts    = counts(order);

tracks   = unique(track_idx);
n_tracks = numel(tracks);
fprintf("Found %d tracks.\n", n_tracks);

%--------------------------------------------
% 2. 结果结构体（binding / bleach 时间点）
%--------------------------------------------
results = struct( ...
    "track_id",        {}, ...
    "binding_indices", {}, ...
    "binding_times",   {}, ...
    "bleach_indices",  {}, ...
    "bleach_times",    {} );

% 用于 dwell time 收集的数组（跨所有 track）
all_binding_times = [];   % 所有 track 的 binding time 合并
all_bleach_times  = [];   % 所有 track 的 bleach time 合并

%--------------------------------------------
% 3. 循环每个 track，跑 HMM + ICON 分析 + 事件检测
%--------------------------------------------
for ti = 1:n_tracks

    tr = tracks(ti);
    fprintf("\n===== Track %d =====\n", tr);

    sel = (track_idx == tr);
    t   = time_sec(sel);
    z   = counts(sel);
    z   = z(:); % 列向量

    %-------------------------
    % 3.1 设定 ICON HMM 的参数
    %-------------------------
    opts        = struct();
    % 超参数（与 level_analysis.m 一致）
    opts.a      = 1;
    opts.g      = 1;
    opts.Q(1)   = mean(z);
    opts.Q(2)   = 1 / (std(z)^2 + eps); % 防止 std(z)=0 除零
    opts.Q(3)   = 0.1;
    opts.Q(4)   = 0.00001;
    opts.M      = 10;

    % 采样参数
    opts.dr_sk  = 1;
    opts.K_init = 50;

    flag_stat   = true;
    flag_anim   = false;  % 建议在 Octave 里关掉动画
    R           = 1000;   % 采样次数，可按需求调整

    %-------------------------
    % 3.2 运行采样器（ICON HMM）
    %-------------------------
    fprintf(" Running HMM sampler...\n");
    chain = chainer_main(z, R, [], opts, flag_stat, flag_anim);

    %-------------------------
    % 3.3 ICON 后验分析（与 level_analysis.m 等价）
    %-------------------------
    fr    = 0.25; % burn-in 比例
    dr    = 2;    % sample 步距
    m_min = 0;
    m_max = 1;
    m_num = 25;

    % (1) 均值轨迹：m_mod
    [m_mod, m_red] = chainer_analyze_means(chain, fr, dr, m_min, m_max, m_num, z);
    m_mod = m_mod(:);

    % (2) 导出 samples（为每条 track 单独存一个文件）
    sample_file = sprintf("samples_track_%d", tr);
    chainer_export(chain, fr, dr, sample_file, "mat");
    fprintf("%s --- Exported\n", [sample_file ".mat"]);

    % (3) 转移概率 / 跃迁统计
    [m_edges, p_mean, d_dist] = chainer_analyze_transitions( ...
        chain, fr, dr, m_min, m_max, m_num, true);

    % (4) 漂移分析
    [y_mean, y_std] = chainer_analyze_drift(chain, fr, dr, z);

    % 保存这些 ICON 分析结果（可选）
    mat_out = sprintf("icon_analysis_track_%d.mat", tr);
    save(mat_out, "m_mod", "m_red", "m_edges", "p_mean", "d_dist", ...
                  "y_mean", "y_std", "t", "z");

    %-------------------------
    % 3.4 把 m_mod 聚类成 3 个光子强度状态（用于事件检测）
    %-------------------------
    K = 3; % 可按物理需求修改

    try
        [idx_raw, centers] = kmeans(m_mod, K);
    catch
        warning("kmeans not available, using simple quantile-based clustering.");
        q1 = quantile(m_mod, 1/3);
        q2 = quantile(m_mod, 2/3);

        idx_raw = ones(size(m_mod));
        idx_raw(m_mod > q1 & m_mod <= q2) = 2;
        idx_raw(m_mod > q2)              = 3;

        centers = zeros(K,1);
        for kk = 1:K
            centers(kk) = mean(m_mod(idx_raw == kk));
        end
    end

    % 根据中心值从小到大重排状态编号，使 state_seq = 1..K 对应 low->high
    [~, order_centers] = sort(centers);
    state_seq = zeros(size(idx_raw));
    for k = 1:K
        state_seq(idx_raw == order_centers(k)) = k;
    end

    low_state  = 1;
    high_state = K;

    %-------------------------
    % 3.5 检测 binding / bleaching
    %
    %   binding: low_state -> high_state
    %   bleach : high_state -> low_state
    %-------------------------
    s      = state_seq(:);
    s_prev = s(1:end-1);
    s_next = s(2:end);

    bind_idx   = find(s_prev == low_state  & s_next == high_state) + 1;
    bleach_idx = find(s_prev == high_state & s_next == low_state) + 1;

    bind_times   = t(bind_idx);
    bleach_times = t(bleach_idx);

    fprintf("Found %d binding event(s) and %d bleaching event(s).\n", ...
            numel(bind_idx), numel(bleach_idx));

    % 存入结果结构体
    results(ti).track_id        = tr;
    results(ti).binding_indices = bind_idx;
    results(ti).binding_times   = bind_times;
    results(ti).bleach_indices  = bleach_idx;
    results(ti).bleach_times    = bleach_times;

    % 用于全局 dwell-time 统计
    all_binding_times = [all_binding_times; bind_times(:)];
    all_bleach_times  = [all_bleach_times;  bleach_times(:)];

end

%--------------------------------------------
% 4. 输出 binding/bleach 时间点的总表（最重要输出）
%--------------------------------------------
out_csv = "binding_bleach_events.csv";
fid_out = fopen(out_csv, "w");
if fid_out < 0
    error("Cannot open output file: %s", out_csv);
end

fprintf(fid_out, "track_index,event_type,sample_index,time_seconds\n");

for ti = 1:numel(results)
    tr = results(ti).track_id;

    % binding
    for k = 1:numel(results(ti).binding_indices)
        fprintf(fid_out, "%d,binding,%d,%.6f\n", tr, ...
                results(ti).binding_indices(k), ...
                results(ti).binding_times(k));
    end

    % bleach
    for k = 1:numel(results(ti).bleach_indices)
        fprintf(fid_out, "%d,bleach,%d,%.6f\n", tr, ...
                results(ti).bleach_indices(k), ...
                results(ti).bleach_times(k));
    end
end

fclose(fid_out);
fprintf("\n[INFO] Binding/bleach events written to: %s\n", out_csv);

%--------------------------------------------
% 5. Dwell-time 直方图 + rate 拟合
%
%   简化假设：
%   - 把每个 track 中连续 binding 事件之间的时间差视为
%     bound state dwell time（中间没有 bleaching，就重新 binding）。
%   - 如果有明确的最后一个 binding 到 bleaching 的时间，
%     可以扩展此处加入（当前先保留简单版本）。
%--------------------------------------------

% 5.1 收集 bound-state dwell times（简单用 binding 时间差）
dwell_times = [];

for ti = 1:numel(results)
    bt = sort(results(ti).binding_times(:)); % 单个 track 的 binding times
    if numel(bt) >= 2
        dt = diff(bt);    % 相邻 binding 之间的时间差
        dwell_times = [dwell_times; dt];
    end

    % 可选：如果你想把最后一个 binding -> 第一个 bleach 也算入 dwell:
    % if ~isempty(results(ti).binding_times) && ~isempty(results(ti).bleach_times)
    %     last_binding  = max(results(ti).binding_times);
    %     first_bleach  = min(results(ti).bleach_times);
    %     if first_bleach > last_binding
    %         dwell_times = [dwell_times; first_bleach - last_binding];
    %     end
    % end
end

if isempty(dwell_times)
    fprintf("[WARN] No sufficient binding events to compute dwell-time distribution.\n");
else
    % 5.2 绘制 dwell-time 直方图（只在内存中画，不弹窗）
    fig = figure("visible", "off");
    hold on;

    nbins = max(10, round(sqrt(numel(dwell_times))));

    % 使用 hist（Octave 自带）代替 histcounts
    [counts_hist, centers] = hist(dwell_times, nbins);

    bar(centers, counts_hist, "hist");
    xlabel("Dwell time (s)");
    ylabel("Counts");
    title("Bound-state dwell-time histogram");

    % 5.3 拟合单指数: p(t) ~ (1/tau) * exp(-t/tau)
    valid = counts_hist > 0;
    t_fit = centers(valid);
    y_fit = log(counts_hist(valid));

    p = polyfit(t_fit, y_fit, 1);
    a = p(1);
    b = p(2);

    tau_est = -1 / b;
    fprintf("[INFO] Fitted single-exponential dwell time tau = %.4f s\n", tau_est);

    t_plot = linspace(min(dwell_times), max(dwell_times), 200);
    y_model = exp(a + b * t_plot);

    plot(t_plot, y_model, "r-", "LineWidth", 2);
    legend("Histogram", "Single-exp fit");
    hold off;

    % 保存图像并关闭 figure（类似 R 的 dev.off()）
    out_png = "dwell_time_histogram_fit.png";
    print(fig, out_png, "-dpng");
    close(fig);

    fprintf("[INFO] Dwell-time histogram and fit saved to %s\n", out_png);
end

fprintf("\nAll done.\n");

运行方式：

octave detect_binding_bleach_dwell.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv

如有下一步你想基于 d_dist 做更“纯 HMM”的 dwell 分布 / rate 矩阵，也可以在这个脚本上再加一节专门分析。

---

level_analysis.m vs detect_binding_bleach.m

这两个脚本在核心 HMM 采样算法上是一致的，但整体流程和用途并不完全相同。level_analysis.m只是做“ICON HMM 拟合 + 若干通用后验分析”，而 detect_binding_bleach.m在此基础上，额外做了“多 track 输入 + 聚类成离散状态 + 自动识别 binding/bleach 跳变 + 输出 CSV”。因此不能说是完全相同算法的简单移植，而是“在相同 HMM 内核上扩展出的一个更完整事件检测流程”。¹²

下面分两部分说明：先列出关键差异，然后用中文系统讲解 detect_binding_bleach.m 的整个处理流程（包括和原始 level_analysis.m 共享的 HMM 部分）。²¹

---

一、算法差异概览

1. 输入数据形式不同
   • level_analysis.m：直接从 photoncount.txt 读入一个向量 z，视为单一时间序列（一个 DNA tether 的 photon counts）。¹
   • detect_binding_bleach.m：从命令行传入一个分号分隔的 CSV，结构为 track_index;time;intensity，先按 (track_index, time) 排序，然后对每个 track 单独取出 z 与时间 t 做分析，相当于一次运行中批量处理多条轨迹。²
2. HMM 部分是否一致
   • 两者都使用 chainer_main 这一 ICON HMM 采样器，参数设置非常接近：
     • opts.a = 1; opts.g = 1;（浓度参数）。¹²
     • opts.Q(1) = mean(z);（均值超分布的均值）。²¹
     • opts.Q(2) = 1/std(z)^2; vs 1 / (std(z)^2 + eps)（只是在 Octave 版本里加了 eps 防除零，对数值稳定性有轻微改进，但理论上等价）。¹²
     • opts.Q(3) = 0.1; opts.Q(4) = 0.00001; opts.M = 10; 完全一致。²¹
     • 采样设置：opts.dr_sk = 1; opts.K_init = 50; R = 1000; 也一致。¹²
   • 因此，给定同一 z，ICON HMM 的抽样和后验均值轨迹的生成逻辑是相同的，只是 Octave 版本做了防数值异常的小修正。整体可以认为 HMM 核心算法是等价的。²¹
3. 后验分析部分的差异
   • level_analysis.m：侧重“通用分析”，调用：
     • chainer_export 导出样本。¹
     • chainer_analyze_means 得到状态均值轨迹 m_mod。¹
     • chainer_analyze_transitions 得到转移概率矩阵、边界等。¹
     • chainer_analyze_drift 得到漂移轨迹 y_mean, y_std。¹
   • detect_binding_bleach.m：只调用 chainer_analyze_means 得到 m_mod，然后完全走另一条分析线：
     • 用 kmeans 或分位数把 m_mod 聚成 3 个强度状态（低、中、高）。²
     • 通过状态序列中的 low→high 和 high→low 跳变来定义 binding 和 bleaching 事件，并输出事件时间。²
   • 所以 HMM 后处理的目标和输出是完全不同的：一个是做统计分析、可视化；另一个是做自动事件检测。
4. 输出结果不同
   • level_analysis.m：
     • 导出 MCMC 样本 samples.mat，方便后续任意分析；同时有若干变量在 workspace 中（m_mod, p_mean, d_dist, y_mean 等）。¹
   • detect_binding_bleach.m：
     • 对每条 track 计算并存储 binding_indices/binding_times 和 bleach_indices/bleach_times，最后全部写入统一的 binding_bleach_events.csv，用于后续统计或可视化。²
5. 附加逻辑（聚类与事件判定）是新加的
   • kmeans（或分位数替代）、状态重标号、状态序列差分、事件 CSV 输出等逻辑，在 level_analysis.m 中完全不存在，是在 Octave 版本里新增的功能层。也就是说，detect_binding_bleach.m 不只是移植，而是在 HMM 内核基础上的一个“面向 binding/bleach 识别”的上层算法。²¹

结论：

• HMM 模型和采样部分与原 MATLAB 代码在理论上是相同的，只做了数值稳定性改动。
• 整体脚本逻辑不相同，Octave 版多了一整套“轨迹遍历 + 状态聚类 + binding/bleach 检测 + CSV 输出”的算法层。

---

二、用中文说明 detect_binding_bleach.m 的完整流程（含与原脚本共用的 HMM 思路）

下面把 detect_binding_bleach.m 当作主角，用中文分步骤解释其处理流程，并在相关步骤指出它与原 level_analysis.m 的对应关系。这样方便你检查是否满足你对“算法一致性”的期待。²¹

---

1. 命令行输入与文件读取

• 脚本通过 argv() 读取命令行参数，第一个参数是输入文件名，例如： octave detect_binding_bleach.m p853_250706_p502_10pN_ch5_0bar_b3_1_track_intensity_data_blue_5s.csv。²
• 输入文件是以分号分隔的文本，第一行为注释头，后面每行包含三列：
  • track_index：轨迹编号（一个 DNA tether 或一个粒子）。
  • time_seconds：该轨迹下某一帧的时间点。
  • track_intensity：该时间点的 photon counts。
• 读入后，把三列分别存成 track_idx, time_sec, counts。然后用 sortrows([track_idx, time_sec]) 按 “轨迹号优先，再按时间” 排好顺序，保证每条轨迹都是时间有序的。之后提取所有不同的 tracks，准备对每一条轨迹单独做分析。²

这一步在原 level_analysis.m 中要简单得多，那里只是 z = load('photoncount.txt');，默认只处理单一系列，没有 track 结构。¹

---

2. 对每条轨迹循环，构建输入序列 z

• 外层 for ti = 1:n_tracks 遍历每一个 track。
• 对于当前轨迹 tr，用 sel = (track_idx == tr) 挑出这一条轨迹对应的时间和强度，得到：
  • t = time_sec(sel)：该轨迹的时间轴。
  • z = counts(sel)：该轨迹的 photon count 序列。
• 把 z 转成列向量，作为后续 HMM 的观测数据。

到这里为止，z 的角色与原脚本中的 z = load('photoncount.txt') 完全一致，只是这里是“多轨迹版本”的 z。¹²

---

3. 设置 ICON HMM 的超参数和采样参数

这一段是两份脚本的核心重合部分，也是你最关心的“算法一致性”所在。¹²

• 先构建 opts 结构体，用于向 chainer_main 传参。
• 超参数设置为：
  • opts.a = 1;：控制 HMM 中状态转换的 Dirichlet 浓度（类似“转移稀疏度”的先验强度）。
  • opts.g = 1;：控制状态先验分布的基本参数（base measure 的浓度）。
  • opts.Q(1) = mean(z);：假设状态均值的先验均值等于观测平均值。
  • opts.Q(2) = 1 / (std(z)^2 + eps);：状态均值先验的精度，与观测标准差成反比；加上 eps 是为了避免在极端情形下 std(z) = 0 导致除零。原 MATLAB 版本是 1/std(z)^2，数学上是一致的设定。
  • opts.Q(3) = 0.1;：方差逆数的 Gamma 先验的 shape。
  • opts.Q(4) = 0.00001;：方差逆数的 Gamma 先验的 scale。
  • opts.M = 10;：插值节点数（ICON 模型内部数值近似需要）。
• 采样相关的参数：
  • opts.dr_sk = 1;：保存样本的间隔步长。
  • opts.K_init = 50;：初始假设的隐状态数（实际有效状态数由采样自动调整）。
  • R = 1000;：采样迭代次数。

这些设置在 level_analysis.m 中完全一致，只是缺少 eps。所以从贝叶斯 HMM 模型结构和先验设定来看，两者是同一个 ICON 模型。²¹

---

4. 运行 HMM 采样器，得到 MCMC 链

• 调用 chain = chainer_main(z, R, [], opts, flag_stat, flag_anim); 来运行 ICON HMM 的 Gibbs sampler 或类似的 MCMC 算法。
  • z 是观测序列。
  • R 是采样轮数。
  • 第三个参数为空 []，表示不提供初始链，由函数内部自行初始化。
  • opts 是前面设定好的超参数和控制参数。
  • flag_stat = true 打开进度输出，flag_anim = false 关闭动画（相比原 MATLAB 代码那里默认 flag_anim = true，这是在 Octave 里为了稳定性和性能做的小改动，对推断结果没有理论影响）。¹²

chainer_main 内部会对每个时间点的隐状态、每个状态的均值和方差、转移矩阵等进行采样，得到一条或多条后验链 chain，用来近似整个后验分布。这个步骤在两个脚本里完全相同，只是 Octave 版本对可视化做了简化。²¹

---

5. 后验均值轨迹分析：chainer_analyze_means

• 设置后处理参数：
  • fr = 0.25：前 25% 的样本视为 burn-in，丢弃不用。
  • dr = 2：之后每隔 2 个样本取一个，用来减弱自相关。
  • m_min = 0; m_max = 1; m_num = 25;：把状态均值的可能范围离散成 $[0,1]$ 区间上 25 个格点，用于后验平均的插值或投影。
• 调用： [m_mod, m_red] = chainer_analyze_means(chain, fr, dr, m_min, m_max, m_num, z);
• 输出 m_mod 可以理解为“在每一个时间点上的后验平均光子水平”，往往已被归一化到 $[0,1]$ 区间，是一个与 z 等长的一维轨迹。m_red 是某种降维或简化形式，这里没有用到。¹²

这一部分与 level_analysis.m 完全一致，区别在于：原脚本后面还继续做转移概率估计和漂移分析；Octave 版则转向“离散状态 + 事件检测”。²¹

---

6. 把连续水平轨迹聚类成几个离散状态

接下来是 Octave 脚本中新加的上层算法部分，用来从连续的 m_mod 中自动提取“低/中/高”强度状态：

• 设定状态数 K = 3（可手动改成 2 或其他）。
• 首选使用 kmeans(m_mod, K) 对 m_mod 做聚类，得到：
  • idx_raw：每个时间点所属的簇编号。
  • centers：每个簇的中心值（大致对应不同光子水平）。
• 如果 Octave 没有统计包或 kmeans 不可用，就退化为基于分位数的简单阈值法：
  • 用三分位数 q1 = quantile(m_mod, 1/3); q2 = quantile(m_mod, 2/3); 把轨迹分为低、中、高三段。
  • 然后为每类算一个均值当作 centers。
• 得到 centers 后，按照中心值从小到大排序，重排状态编号，使得：
  • 状态 1：最低强度（low_state）。
  • 状态 K（=3）：最高强度（high_state）。

这一段在 level_analysis.m 中没有任何对应逻辑，是 Octave 版本为 binding/bleach 检测专门设计的。其本质是：用连续的后验均值轨迹做无监督聚类，把连续空间映射为少数几个离散层级状态。¹²

---

7. 从状态序列中检测 binding 和 bleaching 事件

有了状态序列 state_seq 后，就可以做简单的“状态跳变检测”来定义事件：

• 令 s = state_seq(:)。
• 构造前后相邻状态：
  • s_prev = s(1:end-1);
  • s_next = s(2:end);
• binding 事件：
  • 当某一时间点出现 s_prev == low_state 且 s_next == high_state 时，即从低强度直接跳到高强度，视为一次“分子结合”（binding）。
  • 对应的索引为 find(...) + 1，因为跳变是在第二个点发生。
• bleach 事件：
  • 当出现 s_prev == high_state 且 s_next == low_state 时，即从高强度直接掉回低强度，视为一次“光漂白或解离”（bleaching/unbinding）。
  • 索引同样是 find(...) + 1。
• 最后用 t(bind_idx) 和 t(bleach_idx) 把这些索引转换成真实时间点（秒），分别作为 binding_times 和 bleach_times。

这一步把 HMM 的连续后验轨迹，最终转化成了一个离散事件时间列表。在原 level_analysis.m 中，没有这类应用导向的 event detection，只停留在“描述状态分布和转移概率”的层面。²¹

---

8. 汇总所有轨迹结果并写出 CSV

• 在 for ti = 1:n_tracks 循环内，把每条轨迹的事件信息存入 results(ti)：
  • track_id
  • binding_indices, binding_times
  • bleach_indices, bleach_times。²
• 循环结束后，统一写出一个 CSV 文件 binding_bleach_events.csv，格式为：
  • track_index,event_type,sample_index,time_seconds
  • 其中 event_type 是 "binding" 或 "bleach"。
• 最终在命令行打印一行提示，告诉你输出文件名。

相较之下，level_analysis.m 用的是 chainer_export 把整条 MCMC 链以 .mat 形式保存，供交互式分析和可视化使用，而不是面向“事件列表”的输出。¹²

---

总结（中文角度）

• 从 HMM 模型和抽样的角度看： 两个脚本在 ICON HMM 的参数设定和采样过程上是同一套算法，仅在 Octave 版加入了 eps 防止除零，以及关闭动画提高稳定性，这不会改变算法本质。
• 从整个脚本的“算法流程”看： level_analysis.m 是一个“通用 ICON HMM 演示与分析脚本”，只做单一序列的建模和一些统计性后验分析，不涉及具体事件检测逻辑。 detect_binding_bleach.m 在沿用同一 HMM 内核的基础上，增加了：多轨迹处理、对后验均值轨迹的聚类、将轨迹离散为若干强度状态、基于状态跳变的 binding/bleach 自动识别以及 CSV 输出。这是一层新的、面向应用的算法，不存在于原 MATLAB 脚本中。

如果你的目标是“在 Octave 中复用原 ICON HMM 的统计建模部分，并进一步实现自动 binding/bleaching 检测”，那么当前 detect_binding_bleach.m 的设计是符合这一思路的：核心推断算法相同，外层分析逻辑则是针对具体实验需求进行了扩展。