function result = paper_simple_search(stage, indices)
% Simple changes only. Select on original seeds, freeze before holdout.
% Run single 1:14, pair 1:52, holdout 1:15, develop 1:8, final 1:16.
here = fileparts(mfilename('fullpath'));
file = fullfile(here, 'paper_simple_search.mat');
if isfile(file), s = load(file); result = s.result; else, result = struct(); end
% Keep per-case metrics, not redundant full-rate simulation histories.
if isfield(result,'holdout_cases')
    for j = 1:numel(result.holdout_cases)
        for k = 1:numel(result.holdout_cases{j})
            if isstruct(result.holdout_cases{j}{k})
                result.holdout_cases{j}{k} = result.holdout_cases{j}{k}.table;
            end
        end
    end
end
candidates = {
    'Direct drive ratio', struct('drive', true)
    'Averaged drive ratio', struct('drive_avg', true)
    'Drive-ratio recovery', struct('recovery', 'drive')
    'Quiet-bin fit only', struct('quiet_fit', true)
    'Drive-ratio floor', struct('drive_floor', true)
    'Slow power average', struct('eta_q', 0.12)
    'Confirmation gate', struct('gate', true)
    'Lambda guard 0.02', struct('guard', 0.02)
    'Output smoothing', struct('smooth', true)
    'Lambda upward gain 0.32', struct('up_gain', 0.32)
    'Lambda deadband 0.008', struct('deadband', 0.008)
    'No power recovery', struct('recovery', 'none')};
result.candidates = candidates;
switch stage
    case 'single'
        methods = [{'Code as is', struct('code', true); 'Paper as written', struct('guard', 0)}; candidates];
    case 'pair'
        methods = cell(0,2);
        for a = 1:size(candidates,1)-1
            for b = a+1:size(candidates,1)
                % A replacement power branch makes these switches irrelevant.
                if a <= 2 && ismember(b,[2 3 4 5 6 7 12]), continue, end
                if a == 3 && b == 12, continue, end
                opt = candidates{a,2}; other = candidates{b,2};
                for f = fieldnames(other)', opt.(f{1}) = other.(f{1}); end
                methods(end+1,:) = {[candidates{a,1} ' + ' candidates{b,1}],opt}; %#ok<AGROW>
            end
        end
    case 'holdout'
        % Freeze selection using original seeds only, before any holdout runs.
        if ~isfield(result,'winner')
            [~, ix] = min(cellfun(@(x)x.J,result.pair));
            result.winner = result.pair{ix};
        end
        methods = [{'Code as is', struct('code', true); 'Paper as written', struct('guard', 0)}; ...
            candidates; {result.winner.name, result.winner.options}];
    case 'develop'
        % First holdout became development data after the original winner failed.
        ix = [10 12 13 14 15 16 32 42];
        methods = cell(numel(ix),2);
        for j = 1:numel(ix)
            e = result.pair{ix(j)}; methods(j,:) = {e.name,e.options};
        end
    case 'final'
        if ~isfield(result,'final_winner')
            pool = [result.holdout(3:end),result.develop];
            original = [result.single,result.pair];
            scores = zeros(size(pool));
            for j = 1:numel(pool)
                ix = find(cellfun(@(e)strcmp(e.name,pool{j}.name),original),1);
                scores(j) = (original{ix}.J + 3*pool{j}.J)/4;
                % Existing soft gate adds two constants (slope 2, floor 0.1).
                % Retain it as a comparator, exclude it under the one-new-parameter rule.
                if isfield(pool{j}.options,'gate'), scores(j) = Inf; end
            end
            [~,ix] = min(scores);
            result.final_winner = pool{ix};
            result.development_scores = scores;
        end
        methods = [result.holdout_methods; ...
            {result.final_winner.name,result.final_winner.options}];
end
result.([stage '_methods']) = methods;
if nargin < 2, indices = 1:size(methods,1); end
for j = indices
    opt = methods{j,2};
    if ismember(stage,{'holdout','develop'})
        offsets = [1000 2000 3000];
    elseif strcmp(stage,'final')
        offsets = [4000 5000 6000];
    else
        offsets = 0;
    end
    metrics = zeros(numel(offsets),6); cases = cell(size(offsets));
    for k = 1:numel(offsets)
        evalc('R = run_paper_benchmark_seeds(0.4,opt,offsets(k));');
        t = R.table;
        metrics(k,:) = [mean(t.err([1 2 5 6])),mean(t.err([3 4 7 8])), ...
            mean(t.unsafe(5:8)),mean(t.rough([1 2 5 6])),mean(t.travel(5:8)), ...
            mean(t.err)+mean(t.unsafe(5:8))/1000];
        cases{k} = R.table;
    end
    m = mean(metrics,1);
    entry = struct('name',methods{j,1},'options',opt,'errW',m(1),'errN',m(2), ...
        'unsafe',m(3),'rough',m(4),'travel',m(5),'J',m(6), ...
        'offsets',offsets,'metrics',metrics);
    result.(stage){j} = entry;
    result.([stage '_cases']){j} = cases;
    fprintf('%s %d/%d %-60s J %.5f W %.4f N %.4f unsafe %.2f rough %.2f\n', ...
        stage,j,size(methods,1),entry.name,entry.J,entry.errW,entry.errN,entry.unsafe,entry.rough);
    save(file,'result');
end
end