function result = paper_code_bridge_study(stage, indices)
% Ports of actual code operations. No tuning of their numerical constants.
% Ranking: lowest worst active-case unsafe time. Break ties by mean error,
% then roughness. Show all metrics rather than hiding safety inside J.
file = fullfile(fileparts(mfilename('fullpath')),'paper_code_bridge_study.mat');
if isfile(file), s=load(file); result=s.result; else, result=struct(); end
cand = {
 'Power-level normalization',struct('c_norm',true)
 'Arithmetic confirmation gate',struct('c_gate',true)
 'Prior blend before inverse',struct('c_prior',true)
 'Inverse denominator floor',struct('c_den',true)
 'Median of identified margin',struct('c_idmed',true)
 'Median of power margin',struct('c_powmed',true)
 'Actual three reference snapshots',struct('c_ref',true)
 'Accumulate then recover',struct('both',true)
 'Activity-weighted accumulation',struct('c_activity',true)
 'Two-update drive kick',struct('c_kick',true)
 'Cap branch margins',struct('c_cap',true)
 'Output smoothing',struct('smooth',true)
 'Residual-weighted lambda guard',struct('c_guard',true)
 'Residual-weighted lambda deadband',struct('c_deadband',true)
 'Code upward lambda gain',struct('c_up',true)
 'Lambda gain times sqrt(lambda)',struct('c_rate',true)
 'Lambda step limit',struct('c_step',true)
 'Code frequency band to 100 Hz',struct('band',[2 100])
 'Exact spectral coefficient',struct('eta_s',1-exp(-0.125/0.25))
 'Level-dependent fit floor (requires gate)',struct('c_level',true)
 'Exact accumulation coefficient',struct('eta_g',1-exp(-0.125))
 'Exact recovery coefficient',struct('r_g',exp(-0.125/0.4))};
result.candidates=cand;
offsets=0; idle=false;
switch stage
 case 'single'
  methods=[{'Code as is',struct('code',true);'Paper as written',struct('guard',0)};cand];
 case {'two','three'}
  if strcmp(stage,'two'), source='single'; else, source='two'; end
  pool=result.(source);
  if strcmp(source,'single'), pool=pool(3:end); end
  rank=cell2mat(cellfun(@(e)[e.worstUnsafe,e.err,e.rough],pool,'UniformOutput',false)');
  [~,ix]=sortrows(rank,[1 2 3]); chosen=pool{ix(1)};
  result.([stage '_base'])=chosen;
  bases={chosen};
  if strcmp(stage,'three')
   [~,ixerr]=min(cellfun(@(e)e.err,pool));
   bases{2}=pool{ixerr}; result.three_accuracy_base=bases{2};
  end
  methods=cell(0,2);
  for ibase=1:numel(bases)
  chosen=bases{ibase};
  for j=1:size(cand,1)
   other=cand{j,2};
   if any(isfield(chosen.options,fieldnames(other))),continue,end
   opt=chosen.options;
   for f=fieldnames(other)',opt.(f{1})=other.(f{1});end
   if ~isempty(methods) && any(cellfun(@(o)isequaln(orderfields(o),orderfields(opt)),methods(:,2))),continue,end
   methods(end+1,:)={[chosen.name ' + ' cand{j,1}],opt}; %#ok<AGROW>
  end
  end
 case {'unseen','idle'}
  if ~isfield(result,'selected')
   pool=result.three;
   rank=cell2mat(cellfun(@(e)[e.worstUnsafe,e.err,e.rough],pool,'UniformOutput',false)');
   [~,ix]=sortrows(rank,[1 2 3]); result.selected=pool{ix(1)};
  end
  methods={'Code as is',struct('code',true);'Paper as written',struct('guard',0); ...
    result.two_base.name,result.two_base.options; ...
    result.three_base.name,result.three_base.options; ...
    result.selected.name,result.selected.options; ...
    'Previous rho-recovery recommendation',struct('recovery','drive','guard',0.02); ...
    'Rho recovery + slow power average',struct('recovery','drive','eta_q',0.12)};
  base=result.three_accuracy_base;
  methods(end+1,:)={base.name,base.options};
  for f={'c_guard','smooth','c_activity','c_up'}
   opt=base.options; opt.(f{1})=true;
   methods(end+1,:)={[base.name ' + ' f{1}],opt};
  end
  opt=base.options;opt.band=[2 100];
  methods(end+1,:)={[base.name ' + band 100 Hz'],opt};
  if strcmp(stage,'unseen'),offsets=[7000 8000 9000];else,idle=true;end
 case 'extended'
  % Explicitly labelled four/five-change diagnostics, never counted as three.
  base=struct('c_ref',true,'c_gate',true,'c_activity',true,'c_up',true);
  methods={'Reference + gate + activity + upward gain (4)',base};
  base.smooth=true;
  methods(end+1,:)={'Reference + gate + activity + upward gain + smoothing (5)',base};
  allopt=struct();
  for j=1:size(cand,1)
   other=cand{j,2};for f=fieldnames(other)',allopt.(f{1})=other.(f{1});end
  end
  methods(end+1,:)={'All code ports (diagnostic)',allopt};
  offsets=[0 7000 8000 9000];
end
result.([stage '_methods'])=methods;
if nargin<2,indices=1:size(methods,1);end
for j=indices
 opt=methods{j,2}; metrics=zeros(numel(offsets),12); cases=cell(size(offsets));
 for k=1:numel(offsets)
  evalc('R=run_code_bridge_benchmark(0.4,opt,offsets(k),idle);');
  t=R.table; errors=t.err;
  c3=R.cases(3); w=c3.tc>=45 & c3.tc<75;
  idle_bias=NaN; idle_unsafe=NaN;
  if idle
   c7=R.cases(7); wi=c7.tc>=60 & c7.tc<70;
   idle_bias=mean(c3.b_est(wi)-c3.b_true(wi));
   idle_unsafe=100*mean(c7.b_true(wi)<0.38);
  end
  metrics(k,:)=[mean(errors),mean(errors([1 2 5 6])),mean(errors([3 4 7 8])), ...
    t.unsafe(5:8)',mean(t.rough([1 2 5 6])), ...
    mean(c3.b_est(w)),mean(c3.b_true(w)),idle_bias,idle_unsafe];
  cases{k}=R.table;
 end
 m=mean(metrics,1);
 e=struct('name',methods{j,1},'options',opt,'err',m(1),'errW',m(2),'errN',m(3), ...
  'unsafe',m(4:7),'worstUnsafe',max(m(4:7)),'rough',m(8), ...
  'damageEstimate',m(9),'damageTruth',m(10),'idleBias',m(11),'idleUnsafe7',m(12), ...
  'J',m(1)+mean(m(4:7))/1000,'offsets',offsets,'metrics',metrics,'cases',{cases});
 result.(stage){j}=e;
 fprintf('%s %d/%d %-48s err %.4f unsafe [%5.1f %5.1f %5.1f %5.1f] rough %.2f J %.4f\n', ...
  stage,j,size(methods,1),e.name,e.err,e.unsafe,e.rough,e.J);
 save(file,'result');
end
end