Skip to main content

secs_engine/
ba.rs

1//! Bass alignment, ported from `correction_engine.py:624–1095,1371–1428`.
2//! D1 is deliberately not applied: bus delay still penalizes coherent addition.
3
4use crate::chain::{BaRow, PreparedBaTerm, gather, sub_sum, sub_sum_prepared};
5use crate::decay::{BaLateGrid, ba_late_db};
6use crate::fit::{Fitter, xo_target_peq};
7use crate::grid::{BINS, SR_F, design_freqs};
8use crate::input::{BandType, BundleGain, DbaRow, LayoutKind, Peq, SubAllpass, SubMode};
9use crate::measurement::{Prepared, Spectra, onset_index, ramp_phase};
10use crate::notes::Note;
11use crate::setup::{Setup, SubFilter};
12use crate::{C, EngineConfig, EngineError, MeasurementSet, Tuning};
13use secs_compute::CancelToken;
14use secs_dsp::optimize::{DifferentialEvolutionOptions, differential_evolution};
15use secs_dsp::{
16    fft, pairwise_sum, response,
17    smooth::{self, Fraction},
18};
19use std::{f64::consts::PI, sync::Arc};
20
21pub struct BaInput<'a> {
22    pub measurements: &'a MeasurementSet,
23    pub prepared: &'a Prepared,
24    pub setup: &'a Setup,
25    pub config: &'a EngineConfig,
26    pub tuning: &'a Tuning,
27    pub xo_of: &'a [(String, f64)],
28    pub ref_level: f64,
29    pub cancel: &'a CancelToken,
30}
31
32#[derive(Debug, Clone, PartialEq)]
33pub struct ManualPlan {
34    pub fixed: Vec<String>,
35    pub bundles: Vec<(String, Vec<String>)>,
36    pub auto: Vec<String>,
37    pub gain_anchored: bool,
38}
39
40#[derive(Debug, Clone, PartialEq)]
41pub struct LayoutInfo {
42    pub kind: LayoutKind,
43    pub members: Vec<String>,
44    pub gains: Vec<(String, f64)>,
45    pub front: Vec<String>,
46    pub rear: Vec<String>,
47    pub rear_ms: f64,
48    pub auto_ms: Option<f64>,
49    pub trim_ms: f64,
50    pub min_ms: Option<f64>,
51    pub peq: Vec<Peq>,
52    #[cfg(test)]
53    candidates: Vec<(f64, f64)>,
54}
55
56pub struct SubAlignment {
57    pub rows: Vec<BaRow>,
58    pub sub_filters: Vec<(String, SubFilter)>,
59    pub h_subs: Vec<(u32, String, Arc<[C]>)>,
60    pub layout_rel: Vec<(String, Vec<(String, BaRow)>)>,
61    pub layout_info: Vec<(String, LayoutInfo)>,
62    pub sub_level_db: Option<f64>,
63    pub bus_gain: f64,
64    pub bus_peq: Vec<Peq>,
65    pub bus_rolloff_hz: Option<f64>,
66    pub sub_sum_by_pos: Vec<(u32, Arc<[C]>)>,
67    pub plan: Option<ManualPlan>,
68    pub notes: Vec<Note>,
69}
70
71fn check(input: &BaInput<'_>) -> Result<(), EngineError> {
72    if input.cancel.is_cancelled() {
73        Err(EngineError::Cancelled)
74    } else {
75        Ok(())
76    }
77}
78fn row(gain_db: f64, delay_ms: f64, polarity: i8, allpass: Vec<(f64, f64)>) -> BaRow {
79    BaRow {
80        gain_db,
81        delay_ms,
82        polarity,
83        allpass,
84    }
85}
86fn mean(x: &[f64]) -> f64 {
87    pairwise_sum(x) / x.len() as f64
88}
89fn median(mut x: Vec<f64>) -> f64 {
90    x.sort_by(f64::total_cmp);
91    let m = x.len() / 2;
92    if x.len().is_multiple_of(2) {
93        (x[m - 1] + x[m]) / 2.0
94    } else {
95        x[m]
96    }
97}
98fn mag(x: C) -> f64 {
99    x.re.hypot(x.im)
100}
101fn power(x: C) -> f64 {
102    let a = mag(x);
103    a * a
104}
105fn weight(pos: u32, main: u32) -> f64 {
106    if pos == main { 5.0 } else { 1.0 }
107}
108fn base(input: &BaInput<'_>, role: &str) -> f64 {
109    input
110        .prepared
111        .align_delay_ms
112        .iter()
113        .find(|(r, _)| r == role)
114        .unwrap()
115        .1
116}
117fn complete(h: &Spectra, roles: &[String]) -> Vec<u32> {
118    h.positions()
119        .iter()
120        .copied()
121        .filter(|&p| roles.iter().all(|r| h.get(p, r).is_some()))
122        .collect()
123}
124fn sums(
125    h: &Spectra,
126    roles: &[String],
127    rows: &[BaRow],
128    idx: &[usize],
129    positions: &[u32],
130) -> Vec<(u32, Arc<[C]>)> {
131    let f = gather(design_freqs(), idx);
132    positions
133        .iter()
134        .map(|&p| {
135            let cells: Vec<_> = roles
136                .iter()
137                .map(|r| gather(h.get(p, r).unwrap(), idx))
138                .collect();
139            let terms: Vec<_> = cells
140                .iter()
141                .zip(rows)
142                .map(|(h, r)| (h.as_slice(), r))
143                .collect();
144            (p, sub_sum(&terms, &f).into())
145        })
146        .collect()
147}
148fn grid(hi: f64) -> Vec<usize> {
149    let f = design_freqs();
150    let lo = f.partition_point(|f| *f < 20.0);
151    let hi = f.partition_point(|f| *f <= hi) - 1;
152    let start = (lo as f64).log10();
153    let stop = (hi as f64).log10();
154    let mut out: Vec<_> = (0..160)
155        .map(|i| {
156            let v = if i == 159 {
157                stop
158            } else {
159                start + i as f64 * ((stop - start) / 159.0)
160            };
161            10f64.powf(v).round_ties_even() as usize
162        })
163        .collect();
164    out.dedup();
165    out
166}
167
168/// NumPy polyfit: descending Vandermonde, strided column norm, scaled SVD,
169/// rcond = n * eps. A three-column one-sided Jacobi SVD avoids normal equations.
170/// LAPACK's divide-and-conquer SVD has different final bits; the oracle measures
171/// their effect on the objective, rather than claiming bit-identical coefficients.
172struct Quadratic {
173    x: Vec<f64>,
174    columns: [Vec<f64>; 3],
175    right: [[f64; 3]; 3],
176    scale: [f64; 3],
177    norm: [f64; 3],
178}
179impl Quadratic {
180    fn new(x: Vec<f64>) -> Self {
181        let mut columns = [
182            x.iter().map(|v| v * v).collect::<Vec<_>>(),
183            x.clone(),
184            vec![1.0; x.len()],
185        ];
186        // NumPy axis=0 on its C-order Vandermonde is a sequential reduction.
187        let scale = std::array::from_fn(|j| columns[j].iter().fold(0.0, |s, v| s + v * v).sqrt());
188        for j in 0..3 {
189            for v in &mut columns[j] {
190                *v /= scale[j];
191            }
192        }
193        let mut right = [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]];
194        for _ in 0..30 {
195            let mut changed = false;
196            for (p, q) in [(0, 1), (0, 2), (1, 2)] {
197                let a = columns[p].iter().map(|v| v * v).sum::<f64>();
198                let b = columns[q].iter().map(|v| v * v).sum::<f64>();
199                let c = columns[p]
200                    .iter()
201                    .zip(&columns[q])
202                    .map(|(v, w)| v * w)
203                    .sum::<f64>();
204                if c.abs() <= f64::EPSILON * (a * b).sqrt() {
205                    continue;
206                }
207                changed = true;
208                let z = (b - a) / (2.0 * c);
209                let t = z.signum() / (z.abs() + (1.0 + z * z).sqrt());
210                let co = 1.0 / (1.0 + t * t).sqrt();
211                let si = co * t;
212                let (left, right_column) = columns.split_at_mut(q);
213                for (a, b) in left[p].iter_mut().zip(&mut right_column[0]) {
214                    let old_a = *a;
215                    let old_b = *b;
216                    *a = co * old_a - si * old_b;
217                    *b = si * old_a + co * old_b;
218                }
219                for r in &mut right {
220                    let a = r[p];
221                    let b = r[q];
222                    r[p] = co * a - si * b;
223                    r[q] = si * a + co * b;
224                }
225            }
226            if !changed {
227                break;
228            }
229        }
230        let norm = std::array::from_fn(|j| {
231            pairwise_sum(&columns[j].iter().map(|v| v * v).collect::<Vec<_>>())
232        });
233        Self {
234            x,
235            columns,
236            right,
237            scale,
238            norm,
239        }
240    }
241    fn rough(&self, y: &[f64]) -> f64 {
242        let max = self.norm.iter().copied().fold(0.0, f64::max).sqrt();
243        let mut coef = [0.0; 3];
244        for j in 0..3 {
245            if self.norm[j].sqrt() <= self.x.len() as f64 * f64::EPSILON * max {
246                continue;
247            }
248            let dot = pairwise_sum(
249                &self.columns[j]
250                    .iter()
251                    .zip(y)
252                    .map(|(a, b)| a * b)
253                    .collect::<Vec<_>>(),
254            ) / self.norm[j];
255            for (i, c) in coef.iter_mut().enumerate() {
256                *c += self.right[i][j] * dot;
257            }
258        }
259        for (c, s) in coef.iter_mut().zip(self.scale) {
260            *c /= s;
261        }
262        mean(
263            &self
264                .x
265                .iter()
266                .zip(y)
267                .map(|(x, y)| {
268                    let d = y - ((coef[0] * x + coef[1]) * x + coef[2]);
269                    d * d
270                })
271                .collect::<Vec<_>>(),
272        )
273    }
274}
275fn variance(levels: &[Vec<f64>]) -> f64 {
276    if levels.len() <= 1 {
277        return 0.0;
278    }
279    let mut out = vec![0.0; levels[0].len()];
280    for i in 0..out.len() {
281        let m = levels.iter().fold(0.0, |s, l| s + l[i]) / levels.len() as f64;
282        out[i] = levels.iter().fold(0.0, |s, l| s + (l[i] - m).powi(2)) / levels.len() as f64;
283    }
284    mean(&out)
285}
286fn gd_variance(h: &[C], f: &[f64]) -> f64 {
287    let mut phase: Vec<_> = h.iter().map(|h| h.im.atan2(h.re)).collect();
288    let original = phase.clone();
289    let mut correction = 0.0;
290    for i in 1..phase.len() {
291        let d = original[i] - original[i - 1];
292        let mut dm = (d + PI).rem_euclid(2.0 * PI) - PI;
293        if dm == -PI && d > 0.0 {
294            dm = PI;
295        }
296        if d.abs() >= PI {
297            correction += dm - d;
298        }
299        phase[i] += correction;
300    }
301    let x: Vec<_> = f.iter().map(|f| 2.0 * PI * f).collect();
302    let n = x.len();
303    let mut gd = vec![0.0; n];
304    gd[0] = (phase[1] - phase[0]) / (x[1] - x[0]);
305    gd[n - 1] = (phase[n - 1] - phase[n - 2]) / (x[n - 1] - x[n - 2]);
306    for i in 1..n - 1 {
307        let a = x[i] - x[i - 1];
308        let b = x[i + 1] - x[i];
309        gd[i] = (-b / (a * (a + b))) * phase[i - 1]
310            + ((b - a) / (a * b)) * phase[i]
311            + (a / (b * (a + b))) * phase[i + 1];
312    }
313    for v in &mut gd {
314        *v = -*v * 1000.0;
315    }
316    let m = median(gd.clone());
317    mean(&gd.iter().map(|v| (v - m).powi(2)).collect::<Vec<_>>())
318}
319
320fn individual_db(power: &[Vec<f64>]) -> Vec<Vec<f64>> {
321    power
322        .iter()
323        .map(|row| row.iter().map(|p| 10.0 * (p + 1e-30).log10()).collect())
324        .collect()
325}
326
327struct BaSetup {
328    grid: Vec<usize>,
329    positions: Vec<u32>,
330    late: BaLateGrid,
331    raw_late: Vec<f64>,
332    h_grid: Vec<Vec<Vec<C>>>,
333    h_lin: Vec<Vec<Vec<C>>>,
334    f: Vec<f64>,
335    quadratic: Quadratic,
336}
337impl BaSetup {
338    fn new(input: &BaInput<'_>, h: &Spectra, hi: f64) -> Self {
339        let grid = grid(hi);
340        let positions = complete(h, &input.config.subs);
341        let f = gather(design_freqs(), &grid);
342        let late = BaLateGrid::new(hi);
343        let h_grid = positions
344            .iter()
345            .map(|&p| {
346                input
347                    .config
348                    .subs
349                    .iter()
350                    .map(|r| gather(h.get(p, r).unwrap(), &grid))
351                    .collect()
352            })
353            .collect();
354        let h_lin: Vec<Vec<Vec<C>>> = positions
355            .iter()
356            .map(|&p| {
357                input
358                    .config
359                    .subs
360                    .iter()
361                    .map(|r| h.get(p, r).unwrap()[..late.lin].to_vec())
362                    .collect()
363            })
364            .collect();
365        let rows: Vec<_> = input
366            .config
367            .subs
368            .iter()
369            .map(|r| row(0.0, base(input, r), 1, vec![]))
370            .collect();
371        let raw_late = h_lin
372            .iter()
373            .map(|hs| {
374                let terms: Vec<_> = hs
375                    .iter()
376                    .zip(&rows)
377                    .map(|(h, r)| (h.as_slice(), r))
378                    .collect();
379                ba_late_db(&late, &sub_sum(&terms, &design_freqs()[..late.lin]))
380            })
381            .collect();
382        let quadratic = Quadratic::new(f.iter().map(|f| f.ln()).collect());
383        Self {
384            grid,
385            positions,
386            late,
387            raw_late,
388            h_grid,
389            h_lin,
390            f,
391            quadratic,
392        }
393    }
394    fn indiv(&self, rows: Option<&[BaRow]>) -> Vec<Vec<f64>> {
395        self.h_grid
396            .iter()
397            .map(|hs| {
398                let mut total = vec![0.0; self.grid.len()];
399                for (j, h) in hs.iter().enumerate() {
400                    let g = rows.map_or(1.0, |r| 10f64.powf(r[j].gain_db / 10.0));
401                    for (s, &v) in total.iter_mut().zip(h) {
402                        *s += power(v) * g;
403                    }
404                }
405                total
406            })
407            .collect()
408    }
409    // Existing cost-oracle callers supply powers. Production precomputes the
410    // fixed reference in dB once before entering differential evolution.
411    #[cfg(test)]
412    fn cost(
413        &self,
414        input: &BaInput<'_>,
415        rows: &[BaRow],
416        fixed_indiv: Option<&[Vec<f64>]>,
417        bus: Option<&[C]>,
418    ) -> f64 {
419        let fixed_db = fixed_indiv.map(individual_db);
420        self.cost_with_db(input, rows, fixed_db.as_deref(), bus)
421    }
422
423    fn cost_with_db(
424        &self,
425        input: &BaInput<'_>,
426        rows: &[BaRow],
427        fixed_indiv: Option<&[Vec<f64>]>,
428        bus: Option<&[C]>,
429    ) -> f64 {
430        if input.cancel.is_cancelled() {
431            return 1e30;
432        }
433        let owned;
434        let indiv = if let Some(i) = fixed_indiv {
435            i
436        } else {
437            owned = individual_db(&self.indiv(Some(rows)));
438            &owned
439        };
440        let grid_terms: Vec<_> = rows
441            .iter()
442            .map(|r| PreparedBaTerm::new(&self.f, r))
443            .collect();
444        let lin_terms: Vec<_> = rows
445            .iter()
446            .map(|r| PreparedBaTerm::new(&design_freqs()[..self.late.lin], r))
447            .collect();
448        let mut levels = Vec::new();
449        let mut total = 0.0;
450        let mut weights = 0.0;
451        for (k, &pos) in self.positions.iter().enumerate() {
452            let sum = sub_sum_prepared(&self.h_grid[k], &grid_terms, self.f.len());
453            let db: Vec<_> = sum
454                .iter()
455                .map(|&h| 20.0 * (mag(h) + 1e-20).log10())
456                .collect();
457            let ind = &indiv[k];
458            let rough = self.quadratic.rough(&db);
459            let loss = mean(
460                &ind.iter()
461                    .zip(&db)
462                    .map(|(i, d)| (i - d - 3.0).max(0.0).powi(2))
463                    .collect::<Vec<_>>(),
464            );
465            let eff = mean(&db.iter().zip(ind).map(|(d, i)| d - i).collect::<Vec<_>>());
466            let with_bus: Vec<_> = sum
467                .iter()
468                .zip(&self.grid)
469                .map(|(&s, &i)| bus.map_or(s, |b| s * b[i]))
470                .collect();
471            let gd = gd_variance(&with_bus, &self.f);
472            let ls = sub_sum_prepared(&self.h_lin[k], &lin_terms, self.late.lin);
473            let worse = (ba_late_db(&self.late, &ls) - self.raw_late[k]).max(0.0);
474            let w = weight(pos, input.config.main_position);
475            weights += w;
476            total += w
477                * (rough + 0.5 * loss + 0.02 * gd + 0.5 * worse * worse
478                    - input.tuning.ba_eff_weight * eff);
479            levels.push(db);
480        }
481        total / weights + 0.5 * variance(&levels)
482    }
483}
484
485fn manual_value(setup: &Setup, role: &str) -> BaRow {
486    setup
487        .sub_manual
488        .iter()
489        .find(|(r, _)| r == role)
490        .map_or_else(
491            || row(0.0, 0.0, 1, vec![]),
492            |(_, m)| {
493                row(
494                    m.gain_db,
495                    m.delay_ms,
496                    if m.polarity < 0 { -1 } else { 1 },
497                    m.allpass.iter().take(2).copied().collect(),
498                )
499            },
500        )
501}
502fn manual_plan(input: &BaInput<'_>) -> Option<ManualPlan> {
503    let mut p = ManualPlan {
504        fixed: vec![],
505        bundles: vec![],
506        auto: vec![],
507        gain_anchored: false,
508    };
509    for r in &input.config.subs {
510        match input
511            .setup
512            .sub_manual
513            .iter()
514            .find(|(s, _)| s == r)
515            .map(|(_, m)| m)
516        {
517            Some(m) if m.mode == SubMode::Fixed => p.fixed.push(r.clone()),
518            Some(m) if m.mode == SubMode::Bundle => {
519                let name = if m.bundle.is_empty() { "A" } else { &m.bundle };
520                if let Some((_, rs)) = p.bundles.iter_mut().find(|(b, _)| b == name) {
521                    rs.push(r.clone());
522                } else {
523                    p.bundles.push((name.to_owned(), vec![r.clone()]));
524                }
525            }
526            _ => p.auto.push(r.clone()),
527        }
528    }
529    if p.fixed.is_empty() && p.bundles.is_empty() {
530        return None;
531    }
532    p.gain_anchored = !p.fixed.is_empty()
533        || (!p.bundles.is_empty() && input.tuning.bundle_gain == BundleGain::Fixed);
534    Some(p)
535}
536
537struct Unit {
538    roles: Vec<(String, BaRow)>,
539    center: f64,
540    gain_var: bool,
541    ap_var: bool,
542    auto: bool,
543}
544struct Search {
545    bounds: Vec<(f64, f64)>,
546    units: Vec<Unit>,
547    constants: Vec<(String, BaRow)>,
548    fixed_indiv: bool,
549    skip: bool,
550}
551impl Search {
552    fn new(input: &BaInput<'_>, alignment: &SubAlignment) -> Self {
553        let qmax = if input.tuning.sub_allpass == SubAllpass::LowQ {
554            1.0
555        } else {
556            3.0
557        };
558        let aps = [(20.0, 150.0), (0.3, qmax), (20.0, 150.0), (0.3, qmax)];
559        let mut out = Self {
560            bounds: vec![],
561            units: vec![],
562            constants: vec![],
563            fixed_indiv: true,
564            skip: false,
565        };
566        let Some(plan) = &alignment.plan else {
567            for r in &input.config.subs {
568                out.bounds.extend([
569                    (-8.0, 0.0),
570                    (-input.tuning.ba_delay_ms, input.tuning.ba_delay_ms),
571                    (0.0, 1.0),
572                ]);
573                out.bounds.extend(aps);
574                out.units.push(Unit {
575                    roles: vec![(r.clone(), row(0.0, 0.0, 1, vec![]))],
576                    center: base(input, r),
577                    gain_var: true,
578                    ap_var: true,
579                    auto: true,
580                });
581            }
582            out.skip = input.config.subs.len() == 1;
583            return out;
584        };
585        let extra = if plan.fixed.is_empty() {
586            0.0
587        } else {
588            input.config.sub_main_delay_ms
589        };
590        out.constants = plan
591            .fixed
592            .iter()
593            .map(|r| {
594                let mut v = manual_value(input.setup, r);
595                v.delay_ms -= extra;
596                (r.clone(), v)
597            })
598            .collect();
599        let off = if plan.fixed.is_empty() {
600            0.0
601        } else {
602            median(
603                out.constants
604                    .iter()
605                    .map(|(r, v)| v.delay_ms - base(input, r))
606                    .collect(),
607            )
608        };
609        let mut anchored = plan.fixed.clone();
610        if input.tuning.bundle_gain == BundleGain::Fixed {
611            anchored.extend(plan.bundles.iter().flat_map(|(_, rs)| rs).cloned());
612        }
613        let reference = (!anchored.is_empty()).then(|| {
614            mean(
615                &anchored
616                    .iter()
617                    .map(|r| manual_value(input.setup, r).gain_db)
618                    .collect::<Vec<_>>(),
619            )
620        });
621        out.fixed_indiv = reference.is_none();
622        let gain_bounds =
623            |rel: f64| reference.map_or((-rel - 8.0, -rel), |r| (r - rel - 8.0, r - rel + 8.0));
624        let delay_bounds = |center: f64| {
625            let lo = if plan.fixed.is_empty() {
626                -input.tuning.ba_delay_ms
627            } else {
628                (-input.tuning.ba_delay_ms).max(-extra - center)
629            };
630            (lo, input.tuning.ba_delay_ms.max(lo + 1.0))
631        };
632        for r in &plan.auto {
633            let center = base(input, r) + off;
634            out.bounds
635                .extend([gain_bounds(0.0), delay_bounds(center), (0.0, 1.0)]);
636            out.bounds.extend(aps);
637            out.units.push(Unit {
638                roles: vec![(r.clone(), row(0.0, 0.0, 1, vec![]))],
639                center,
640                gain_var: true,
641                ap_var: true,
642                auto: true,
643            });
644        }
645        for (name, rs) in &plan.bundles {
646            let layout = alignment.layout_rel.iter().find(|(b, _)| b == name);
647            let mut vals = layout.map_or_else(
648                || {
649                    rs.iter()
650                        .map(|r| (r.clone(), manual_value(input.setup, r)))
651                        .collect()
652                },
653                |(_, v)| v.clone(),
654            );
655            let dmin = vals
656                .iter()
657                .map(|(_, v)| v.delay_ms)
658                .fold(f64::INFINITY, f64::min);
659            let center = vals
660                .iter()
661                .map(|(r, v)| base(input, r) - (v.delay_ms - dmin))
662                .fold(f64::NEG_INFINITY, f64::max)
663                + off;
664            let gain_var = input.tuning.bundle_gain == BundleGain::Auto || layout.is_some();
665            let ap_var = layout.is_some() && input.tuning.sub_allpass != SubAllpass::None;
666            if gain_var {
667                out.bounds.push(gain_bounds(mean(
668                    &vals.iter().map(|(_, v)| v.gain_db).collect::<Vec<_>>(),
669                )));
670            }
671            out.bounds.extend([delay_bounds(center), (0.0, 1.0)]);
672            if ap_var {
673                out.bounds.extend(aps);
674            }
675            for (_, v) in &mut vals {
676                v.delay_ms -= dmin;
677            }
678            out.units.push(Unit {
679                roles: vals,
680                center,
681                gain_var,
682                ap_var,
683                auto: false,
684            });
685        }
686        out.skip = out.units.is_empty() || (out.units.len() == 1 && plan.fixed.is_empty());
687        out
688    }
689    fn unpack(&self, input: &BaInput<'_>, x: &[f64]) -> Vec<BaRow> {
690        let mut out = self.constants.clone();
691        let mut k = 0;
692        for u in &self.units {
693            let gain = if u.gain_var {
694                let v = x[k];
695                k += 1;
696                v
697            } else {
698                0.0
699            };
700            let d = x[k];
701            let pol = if x[k + 1] < 0.5 { 1 } else { -1 };
702            k += 2;
703            let mut aps = vec![];
704            if u.ap_var {
705                if input.tuning.sub_allpass != SubAllpass::None
706                    && (u.auto || (x[k] > 0.0 && x[k + 2] > 0.0))
707                {
708                    aps.extend([(x[k], x[k + 1]), (x[k + 2], x[k + 3])]);
709                }
710                k += 4;
711            }
712            for (r, v) in &u.roles {
713                let mut a = v.allpass.clone();
714                a.extend_from_slice(&aps);
715                out.push((
716                    r.clone(),
717                    row(
718                        v.gain_db + gain,
719                        u.center + v.delay_ms + d,
720                        v.polarity * pol,
721                        a,
722                    ),
723                ));
724            }
725        }
726        input
727            .config
728            .subs
729            .iter()
730            .map(|r| out.iter().find(|(s, _)| s == r).unwrap().1.clone())
731            .collect()
732    }
733}
734
735fn level(input: &BaInput<'_>, by_pos: &[(u32, Arc<[C]>)]) -> f64 {
736    let h = &by_pos
737        .iter()
738        .find(|(p, _)| *p == input.config.main_position)
739        .unwrap_or(&by_pos[0])
740        .1;
741    let indices: Vec<_> = design_freqs()
742        .iter()
743        .enumerate()
744        .filter(|(_, f)| **f >= 30.0 && **f <= 80.0)
745        .map(|(i, _)| i)
746        .collect();
747    let target = input
748        .config
749        .target
750        .sample(&gather(design_freqs(), &indices));
751    median(
752        indices
753            .iter()
754            .zip(target)
755            .map(|(&i, t)| 20.0 * (mag(h[i]) + 1e-20).log10() - t)
756            .collect(),
757    )
758}
759#[allow(clippy::too_many_arguments)]
760fn target_peq(
761    input: &BaInput<'_>,
762    by_pos: &[(u32, Arc<[C]>)],
763    lvl: f64,
764    hi: f64,
765    bands: usize,
766    boost: f64,
767    fmax: Option<f64>,
768    fitter: &mut Fitter<'_>,
769) -> (Vec<Peq>, Option<f64>, crate::low_protect::LimitOutcome) {
770    let rolloff = input
771        .config
772        .low_protect
773        .then(|| crate::low_protect::rolloff(by_pos, input.config.main_position, Some((0.0, 0.0))))
774        .flatten();
775    let by_pos: Vec<_> = by_pos
776        .iter()
777        .map(|(p, h)| (weight(*p, input.config.main_position), h.as_ref()))
778        .collect();
779    let mut peq = xo_target_peq(
780        &by_pos,
781        &input.config.target,
782        lvl,
783        rolloff.unwrap_or(20.0),
784        hi,
785        bands,
786        boost,
787        fmax,
788        fitter,
789    );
790    let outcome = rolloff.map_or_else(Default::default, |hz| {
791        crate::low_protect::limit(
792            &mut peq,
793            &vec![C::new(1.0, 0.0); BINS],
794            hz,
795            Default::default(),
796        )
797    });
798    (peq, rolloff, outcome)
799}
800fn bus_fit(
801    input: &BaInput<'_>,
802    by_pos: &[(u32, Arc<[C]>)],
803    fitter: &mut Fitter<'_>,
804) -> (Vec<Peq>, Option<f64>, crate::low_protect::LimitOutcome) {
805    if !input.tuning.xo_stage || input.xo_of.is_empty() {
806        let rolloff = input
807            .config
808            .low_protect
809            .then(|| {
810                crate::low_protect::rolloff(by_pos, input.config.main_position, Some((0.0, 0.0)))
811            })
812            .flatten();
813        return (vec![], rolloff, Default::default());
814    }
815    let xo: Vec<_> = input.xo_of.iter().map(|(_, v)| *v).collect();
816    let hi = if input.tuning.sub_bus_fix {
817        1.25 * xo.iter().copied().fold(f64::NEG_INFINITY, f64::max)
818    } else {
819        2.0 * mean(&xo)
820    };
821    target_peq(
822        input,
823        by_pos,
824        input.ref_level,
825        hi,
826        10,
827        6.0_f64.min(input.config.max_boost_db),
828        input.tuning.sub_bus_fix.then_some(hi),
829        fitter,
830    )
831}
832fn bus_response(
833    input: &BaInput<'_>,
834    h: &Spectra,
835    alignment: &SubAlignment,
836    rows: &[BaRow],
837    fitter: &mut Fitter<'_>,
838) -> Result<Option<Vec<C>>, EngineError> {
839    check(input)?;
840    if !input.tuning.xo_stage || input.xo_of.is_empty() {
841        return Ok(None);
842    }
843    let mut by_pos = sums(
844        h,
845        &input.config.subs,
846        rows,
847        &(0..BINS).collect::<Vec<_>>(),
848        &complete(h, &input.config.subs),
849    );
850    if !by_pos.iter().any(|(p, _)| *p == input.config.main_position) {
851        return Ok(None);
852    }
853    if !alignment.plan.as_ref().is_some_and(|p| p.gain_anchored) {
854        let gain = 10f64.powf((input.ref_level - level(input, &by_pos)) / 20.0);
855        for (_, h) in &mut by_pos {
856            *h = h.iter().map(|v| *v * C::new(gain, 0.0)).collect();
857        }
858    }
859    let (peq, _, _) = bus_fit(input, &by_pos, fitter);
860    check(input)?;
861    Ok((!peq.is_empty()).then(|| response::peq_response(design_freqs(), &peq)))
862}
863
864// The layout closure is kept separate from BaLateGrid: it is its own upstream
865// definition, even though its current window endpoints happen to be identical.
866fn layout_late(spec: &[C], band: &[f64], fs: f64) -> f64 {
867    let y = fft::irfft(
868        &spec
869            .iter()
870            .zip(band)
871            .map(|(s, b)| s * b)
872            .collect::<Vec<_>>(),
873        2 * (spec.len() - 1),
874    );
875    let top = y.iter().map(|v| v.abs()).fold(0.0, f64::max);
876    let k = y.iter().position(|v| v.abs() > 0.1 * top).unwrap_or(0);
877    let energy = |a: f64, b: f64| {
878        pairwise_sum(
879            &y[(k + (a * fs) as usize).min(y.len())..(k + (b * fs) as usize).min(y.len())]
880                .iter()
881                .map(|v| v * v)
882                .collect::<Vec<_>>(),
883        )
884    };
885    10.0 * (energy(0.1, 0.38) / (energy(0.0, 0.05) + 1e-30) + 1e-30).log10()
886}
887fn layout_stage(
888    input: &BaInput<'_>,
889    h: &mut Spectra,
890    out: &mut SubAlignment,
891    hi: f64,
892    fitter: &mut Fitter<'_>,
893) -> Result<(), EngineError> {
894    let f = design_freqs();
895    let grid = grid(hi);
896    let fg = gather(f, &grid);
897    let quadratic = Quadratic::new(fg.iter().map(|f| f.ln()).collect());
898    let lin = f.partition_point(|f| *f < 200.0);
899    let li: Vec<_> = (0..lin).collect();
900    let band: Vec<_> = f[..lin]
901        .iter()
902        .map(|f| ((f - 15.0) / 5.0).clamp(0.0, 1.0) * ((hi * 1.3 - f) / (hi * 0.3)).clamp(0.0, 1.0))
903        .collect();
904    let mask: Vec<_> = f
905        .iter()
906        .enumerate()
907        .filter(|(_, f)| **f >= 20.0 && **f <= hi)
908        .map(|(i, _)| i)
909        .collect();
910    for (name, lay) in &input.setup.layouts {
911        check(input)?;
912        let positions = complete(h, &lay.members);
913        if positions.is_empty() {
914            return Err(EngineError::Input(format!(
915                "bundle {name}: no complete position"
916            )));
917        }
918        let levels: Vec<_> = lay
919            .members
920            .iter()
921            .map(|r| {
922                let mut pw = vec![0.0; mask.len()];
923                let mut total = 0.0;
924                for &p in h.positions() {
925                    if let Some(h) = h.get(p, r) {
926                        let w = weight(p, input.config.main_position);
927                        total += w;
928                        for (v, &i) in pw.iter_mut().zip(&mask) {
929                            *v += w * power(h[i]);
930                        }
931                    }
932                }
933                median(
934                    pw.iter()
935                        .map(|v| 10.0 * (v / total + 1e-30).log10())
936                        .collect(),
937                )
938            })
939            .collect();
940        let low = levels.iter().copied().fold(f64::INFINITY, f64::min);
941        let gains: Vec<_> = lay
942            .members
943            .iter()
944            .cloned()
945            .zip(levels.iter().map(|v| low - v))
946            .collect();
947        let rear: Vec<_> = lay
948            .members
949            .iter()
950            .filter(|r| {
951                input
952                    .setup
953                    .sub_manual
954                    .iter()
955                    .any(|(s, m)| s == *r && m.dba_row == DbaRow::Rear)
956            })
957            .cloned()
958            .collect();
959        let front: Vec<_> = lay
960            .members
961            .iter()
962            .filter(|r| !rear.contains(r))
963            .cloned()
964            .collect();
965        let mut info = LayoutInfo {
966            kind: lay.kind,
967            members: lay.members.clone(),
968            gains: gains.clone(),
969            front: front.clone(),
970            rear: vec![],
971            rear_ms: 0.0,
972            auto_ms: None,
973            trim_ms: 0.0,
974            min_ms: None,
975            peq: vec![],
976            #[cfg(test)]
977            candidates: vec![],
978        };
979        let mut rows: Vec<_> = gains.iter().map(|(_, g)| row(*g, 0.0, 1, vec![])).collect();
980        if lay.kind == LayoutKind::Dba {
981            if front.is_empty() || rear.is_empty() {
982                info.kind = LayoutKind::Multi;
983                out.notes.push(Note::LayoutFallback(name.clone()));
984            } else {
985                let arrival = |p, r: &str| -> Result<f64, EngineError> {
986                    let m = input
987                        .measurements
988                        .get(p, r)
989                        .ok_or_else(|| EngineError::Input(format!("missing DBA IR {p}/{r}")))?;
990                    Ok(1000.0 * (m.t0 + onset_index(&m.samples) as f64 / SR_F))
991                };
992                let mut gap = 0.0_f64;
993                for &p in &positions {
994                    let r = rear
995                        .iter()
996                        .map(|r| arrival(p, r))
997                        .collect::<Result<Vec<_>, _>>()?;
998                    let fr = front
999                        .iter()
1000                        .map(|r| arrival(p, r))
1001                        .collect::<Result<Vec<_>, _>>()?;
1002                    gap = gap.max((mean(&r) - mean(&fr)).abs());
1003                }
1004                let lo = (gap.min(59.0) * 4.0).ceil() / 4.0;
1005                let mut best = (f64::INFINITY, 0.0);
1006                for k in 0..=((60.0 - lo) * 4.0) as usize {
1007                    check(input)?;
1008                    let delay = lo + k as f64 * 0.25;
1009                    let trial: Vec<_> = gains
1010                        .iter()
1011                        .map(|(r, g)| {
1012                            row(
1013                                *g,
1014                                if rear.contains(r) { delay } else { 0.0 },
1015                                if rear.contains(r) { -1 } else { 1 },
1016                                vec![],
1017                            )
1018                        })
1019                        .collect();
1020                    let gs = sums(h, &lay.members, &trial, &grid, &positions);
1021                    let ls = sums(h, &lay.members, &trial, &li, &positions);
1022                    let mut cost = 0.0;
1023                    let mut wsum = 0.0;
1024                    let mut dbs = vec![];
1025                    for ((p, g), (_, l)) in gs.iter().zip(&ls) {
1026                        let w = weight(*p, input.config.main_position);
1027                        wsum += w;
1028                        let db: Vec<_> =
1029                            g.iter().map(|v| 20.0 * (mag(*v) + 1e-20).log10()).collect();
1030                        cost += w * quadratic.rough(&db);
1031                        cost += w * 0.2 * layout_late(l, &band, 2.0 * f[lin - 1]);
1032                        dbs.push(db);
1033                    }
1034                    cost = cost / wsum + 0.5 * variance(&dbs);
1035                    #[cfg(test)]
1036                    info.candidates.push((delay, cost));
1037                    if cost < best.0 {
1038                        best = (cost, delay);
1039                    }
1040                }
1041                let delay = (best.1 + lay.rear_trim_ms).max(0.0);
1042                for ((r, _), v) in gains.iter().zip(&mut rows) {
1043                    if rear.contains(r) {
1044                        v.delay_ms = delay;
1045                        v.polarity = -1;
1046                    }
1047                }
1048                info.rear = rear;
1049                info.rear_ms = delay;
1050                info.auto_ms = Some(best.1);
1051                info.trim_ms = lay.rear_trim_ms;
1052                info.min_ms = Some(gap);
1053            }
1054        }
1055        let by_pos = sums(
1056            h,
1057            &lay.members,
1058            &rows,
1059            &(0..BINS).collect::<Vec<_>>(),
1060            &positions,
1061        );
1062        let (peq, rolloff, outcome) = target_peq(
1063            input,
1064            &by_pos,
1065            level(input, &by_pos),
1066            hi,
1067            6,
1068            6.0,
1069            None,
1070            fitter,
1071        );
1072        check(input)?;
1073        if let Some(hz) = rolloff.filter(|_| outcome.reduced.len() + outcome.cuts > 0) {
1074            out.notes.push(Note::LowProtect {
1075                role: lay.members.join(" + "),
1076                rolloff_hz: hz,
1077                hp_hz: None,
1078                limited: crate::low_protect::limited(&peq, &vec![C::new(1.0, 0.0); BINS], hz),
1079                reduced: outcome.reduced.len(),
1080                cuts: outcome.cuts,
1081            });
1082        }
1083        if !peq.is_empty() {
1084            let response = response::peq_response(f, &peq);
1085            for r in &lay.members {
1086                if !out.sub_filters.iter().any(|(s, _)| s == r) {
1087                    out.sub_filters.push((
1088                        r.clone(),
1089                        SubFilter {
1090                            hp_hz: 0.0,
1091                            hp_type: BandType::Lr24,
1092                            lp_hz: 0.0,
1093                            lp_type: BandType::Lr24,
1094                            peq: vec![],
1095                        },
1096                    ));
1097                }
1098                out.sub_filters
1099                    .iter_mut()
1100                    .find(|(s, _)| s == r)
1101                    .unwrap()
1102                    .1
1103                    .peq
1104                    .extend_from_slice(&peq);
1105                for p in h.positions().to_vec() {
1106                    if let Some(old) = h.get(p, r) {
1107                        let new: Arc<[C]> = old.iter().zip(&response).map(|(a, b)| a * b).collect();
1108                        h.set(p, r, new.clone());
1109                        out.h_subs.push((p, r.clone(), new));
1110                    }
1111                }
1112            }
1113        }
1114        info.peq = peq;
1115        out.layout_rel.push((
1116            name.clone(),
1117            lay.members.iter().cloned().zip(rows).collect(),
1118        ));
1119        out.layout_info.push((name.clone(), info));
1120        out.notes.push(Note::Layout(name.clone()));
1121    }
1122    Ok(())
1123}
1124
1125fn two_way(input: &BaInput<'_>, h: &Spectra, out: &mut SubAlignment) -> Result<usize, EngineError> {
1126    let filters = &out.sub_filters;
1127    let low: Vec<_> = input
1128        .config
1129        .subs
1130        .iter()
1131        .filter(|r| {
1132            filters
1133                .iter()
1134                .any(|(s, f)| s == *r && f.lp_hz > 0.0 && f.hp_hz == 0.0)
1135        })
1136        .cloned()
1137        .collect();
1138    let high: Vec<_> = input
1139        .config
1140        .subs
1141        .iter()
1142        .filter(|r| !low.contains(r))
1143        .collect();
1144    if low.is_empty() || high.is_empty() {
1145        return Ok(0);
1146    }
1147    let xo = filters.iter().find(|(r, _)| r == &low[0]).unwrap().1.lp_hz;
1148    if high.iter().any(|r| {
1149        (filters
1150            .iter()
1151            .find(|(s, _)| s == *r)
1152            .map_or(0.0, |(_, f)| f.hp_hz)
1153            - xo)
1154            .abs()
1155            > 0.5
1156    }) {
1157        return Ok(0);
1158    }
1159    let movable: Vec<_> = low
1160        .into_iter()
1161        .filter(|r| !out.plan.as_ref().is_some_and(|p| p.fixed.contains(r)))
1162        .collect();
1163    if movable.is_empty() {
1164        return Ok(0);
1165    }
1166    let f = design_freqs();
1167    let lo: Vec<_> = f
1168        .iter()
1169        .enumerate()
1170        .filter(|(_, f)| **f >= 12.0_f64.max(xo * 0.55) && **f <= xo * 0.85)
1171        .map(|(i, _)| i)
1172        .collect();
1173    let hi: Vec<_> = f
1174        .iter()
1175        .enumerate()
1176        .filter(|(_, f)| **f >= xo * 1.3 && **f <= 80.0_f64.min(xo * 3.2))
1177        .map(|(i, _)| i)
1178        .collect();
1179    if lo.is_empty() || hi.is_empty() {
1180        return Ok(0);
1181    }
1182    let target = input
1183        .config
1184        .target
1185        .sample(&f.iter().map(|f| f.max(1.0)).collect::<Vec<_>>());
1186    let mut total = 0.0;
1187    let mut iterations = 0;
1188    for _ in 0..4 {
1189        iterations += 1;
1190        check(input)?;
1191        let sum = sums(
1192            h,
1193            &input.config.subs,
1194            &out.rows,
1195            &(0..BINS).collect::<Vec<_>>(),
1196            &[input.config.main_position],
1197        );
1198        let p: Vec<_> = sum[0].1.iter().map(|h| power(*h)).collect();
1199        let sm = smooth::smooth_power(f, &p, Fraction::Scalar(3.0))?;
1200        let db: Vec<_> = sm
1201            .iter()
1202            .zip(&target)
1203            .map(|(p, t)| 10.0 * (p + 1e-30).log10() - t)
1204            .collect();
1205        let diff = median(gather(&db, &hi)) - median(gather(&db, &lo));
1206        let step = diff.clamp(-24.0 - total, input.config.infra_boost_max_db - total);
1207        if step.abs() < 0.1 {
1208            break;
1209        }
1210        for (r, v) in input.config.subs.iter().zip(&mut out.rows) {
1211            if movable.contains(r) {
1212                v.gain_db += step;
1213            }
1214        }
1215        total += step;
1216    }
1217    if total.abs() >= 0.1 {
1218        out.notes.push(Note::TwoWay {
1219            roles: movable,
1220            gain_db: total,
1221            xo_hz: xo,
1222        });
1223    }
1224    Ok(iterations)
1225}
1226
1227type DeSolver<'a> = dyn FnMut(
1228        &mut dyn FnMut(&[f64]) -> f64,
1229        &[(f64, f64)],
1230        Option<&[f64]>,
1231        DifferentialEvolutionOptions,
1232    ) -> Result<secs_dsp::optimize::DifferentialEvolutionResult, secs_dsp::DspError>
1233    + 'a;
1234
1235pub fn align(input: &BaInput<'_>, fitter: &mut Fitter<'_>) -> Result<SubAlignment, EngineError> {
1236    align_with(input, fitter, &mut |cost, bounds, x0, options| {
1237        differential_evolution(cost, bounds, x0, options)
1238    })
1239}
1240
1241fn align_with(
1242    input: &BaInput<'_>,
1243    fitter: &mut Fitter<'_>,
1244    solver: &mut DeSolver<'_>,
1245) -> Result<SubAlignment, EngineError> {
1246    check(input)?;
1247    let mut out = SubAlignment {
1248        rows: vec![],
1249        sub_filters: input.setup.sub_filters.clone(),
1250        h_subs: vec![],
1251        layout_rel: vec![],
1252        layout_info: vec![],
1253        sub_level_db: None,
1254        bus_gain: 1.0,
1255        bus_peq: vec![],
1256        bus_rolloff_hz: None,
1257        sub_sum_by_pos: vec![],
1258        plan: manual_plan(input),
1259        notes: vec![],
1260    };
1261    if input.config.subs.is_empty() {
1262        return Ok(out);
1263    }
1264    let hi = input
1265        .xo_of
1266        .iter()
1267        .map(|(_, x)| *x)
1268        .reduce(f64::max)
1269        .unwrap_or(100.0)
1270        * 1.5;
1271    let mut h = input.prepared.h.clone();
1272    layout_stage(input, &mut h, &mut out, hi, fitter)?;
1273    let search = Search::new(input, &out);
1274    if out.plan.is_some() {
1275        out.notes.push(Note::BaManual);
1276    }
1277    if out.plan.is_none() && input.config.subs.len() == 1 {
1278        out.rows = vec![row(0.0, base(input, &input.config.subs[0]), 1, vec![])];
1279    } else if search.skip {
1280        out.rows = search.unpack(input, &vec![0.0; search.bounds.len()]);
1281    } else {
1282        let s = BaSetup::new(input, &h, hi);
1283        if s.positions.is_empty() {
1284            return Err(EngineError::Input("BA has no complete position".into()));
1285        }
1286        let indiv = search.fixed_indiv.then(|| individual_db(&s.indiv(None)));
1287        let res = solver(
1288            &mut |x| s.cost_with_db(input, &search.unpack(input, x), indiv.as_deref(), None),
1289            &search.bounds,
1290            None,
1291            DifferentialEvolutionOptions::default(),
1292        )?;
1293        check(input)?;
1294        out.notes.push(Note::BaCost(res.fun));
1295        let mut x = res.x;
1296        if input.tuning.sub_bus_fix
1297            && let Some(p1) = bus_response(input, &h, &out, &search.unpack(input, &x), fitter)?
1298        {
1299            let x0: Vec<_> = x
1300                .iter()
1301                .zip(&search.bounds)
1302                .map(|(x, (a, b))| x.clamp(*a, *b))
1303                .collect();
1304            let res = solver(
1305                &mut |x| {
1306                    s.cost_with_db(input, &search.unpack(input, x), indiv.as_deref(), Some(&p1))
1307                },
1308                &search.bounds,
1309                Some(&x0),
1310                DifferentialEvolutionOptions {
1311                    maxiter: 40,
1312                    seed: 2,
1313                    ..Default::default()
1314                },
1315            )?;
1316            check(input)?;
1317            let b1 = bus_response(input, &h, &out, &search.unpack(input, &x), fitter)?;
1318            let c1 = s.cost_with_db(
1319                input,
1320                &search.unpack(input, &x),
1321                indiv.as_deref(),
1322                b1.as_deref(),
1323            );
1324            let b2 = bus_response(input, &h, &out, &search.unpack(input, &res.x), fitter)?;
1325            let c2 = s.cost_with_db(
1326                input,
1327                &search.unpack(input, &res.x),
1328                indiv.as_deref(),
1329                b2.as_deref(),
1330            );
1331            out.notes.push(Note::BaBusPass {
1332                first: c1,
1333                second: c2,
1334                kept_first: c2 >= c1,
1335            });
1336            if c2 < c1 {
1337                x = res.x;
1338            }
1339        }
1340        out.rows = search.unpack(input, &x);
1341    }
1342    two_way(input, &h, &mut out)?;
1343    check(input)?;
1344    out.sub_sum_by_pos = sums(
1345        &h,
1346        &input.config.subs,
1347        &out.rows,
1348        &(0..BINS).collect::<Vec<_>>(),
1349        &complete(&h, &input.config.subs),
1350    );
1351    if !out
1352        .sub_sum_by_pos
1353        .iter()
1354        .any(|(p, _)| *p == input.config.main_position)
1355    {
1356        return Err(EngineError::Input("missing BA reference position".into()));
1357    }
1358    let lvl = level(input, &out.sub_sum_by_pos);
1359    out.sub_level_db = Some(lvl - input.ref_level);
1360    if out.plan.as_ref().is_some_and(|p| p.gain_anchored) {
1361        out.notes.push(Note::BusGainAnchored(lvl - input.ref_level));
1362    } else {
1363        out.bus_gain = 10f64.powf((input.ref_level - lvl) / 20.0);
1364    }
1365    for (_, h) in &mut out.sub_sum_by_pos {
1366        *h = h.iter().map(|v| *v * C::new(out.bus_gain, 0.0)).collect();
1367    }
1368    let (peq, rolloff, outcome) = bus_fit(input, &out.sub_sum_by_pos, fitter);
1369    out.bus_peq = peq;
1370    out.bus_rolloff_hz = rolloff;
1371    if let Some(hz) = rolloff.filter(|_| outcome.reduced.len() + outcome.cuts > 0) {
1372        out.notes.push(Note::LowProtect {
1373            role: input.config.subs.join(" + "),
1374            rolloff_hz: hz,
1375            hp_hz: None,
1376            limited: crate::low_protect::limited(&out.bus_peq, &vec![C::new(1.0, 0.0); BINS], hz),
1377            reduced: outcome.reduced.len(),
1378            cuts: outcome.cuts,
1379        });
1380    }
1381    check(input)?;
1382    if input.tuning.xo_stage && !input.xo_of.is_empty() {
1383        let p = response::peq_response(design_freqs(), &out.bus_peq);
1384        for (_, h) in &mut out.sub_sum_by_pos {
1385            *h = h.iter().zip(&p).map(|(v, p)| v * p).collect();
1386        }
1387    }
1388    Ok(out)
1389}
1390
1391/// Main spectra include acoustic-XO PEQ, but not their electrical HP or delay.
1392/// Returns `Err(Cancelled)` once `cancelled()` reads true; it is polled once
1393/// per main position while the bands are gathered and once per candidate
1394/// delay, which leaves the search's arithmetic and order unchanged (TD-49).
1395#[allow(clippy::type_complexity)]
1396pub fn bus_delay(
1397    alignment: &SubAlignment,
1398    mains: &[(String, Vec<(u32, Arc<[C]>)>)],
1399    align_delay_ms: &[(String, f64)],
1400    xo_of: &[(String, f64)],
1401    lo_ms: f64,
1402    cancelled: &dyn Fn() -> bool,
1403) -> Result<(f64, i8), EngineError> {
1404    bus_delay_observe(
1405        alignment,
1406        mains,
1407        align_delay_ms,
1408        xo_of,
1409        lo_ms,
1410        cancelled,
1411        &mut |_, _, _| {},
1412    )
1413}
1414
1415#[allow(clippy::type_complexity)]
1416fn bus_delay_observe(
1417    alignment: &SubAlignment,
1418    mains: &[(String, Vec<(u32, Arc<[C]>)>)],
1419    align_delay_ms: &[(String, f64)],
1420    xo_of: &[(String, f64)],
1421    lo_ms: f64,
1422    cancelled: &dyn Fn() -> bool,
1423    observe: &mut dyn FnMut(i8, f64, f64),
1424) -> Result<(f64, i8), EngineError> {
1425    let check = || {
1426        if cancelled() {
1427            Err(EngineError::Cancelled)
1428        } else {
1429            Ok(())
1430        }
1431    };
1432    if alignment.plan.as_ref().is_some_and(|p| !p.fixed.is_empty()) {
1433        return Ok((0.0, 1));
1434    }
1435    let mut prepared = vec![];
1436    for (r, positions) in mains {
1437        let Some(&(_, xo)) = xo_of.iter().find(|(s, x)| s == r && *x != 0.0) else {
1438            continue;
1439        };
1440        let idx: Vec<_> = design_freqs()
1441            .iter()
1442            .enumerate()
1443            .filter(|(_, f)| **f >= xo / 2.0 && **f <= xo * 2.0)
1444            .map(|(i, _)| i)
1445            .collect();
1446        let f = gather(design_freqs(), &idx);
1447        let d = align_delay_ms.iter().find(|(s, _)| s == r).unwrap().1;
1448        let hp = response::lr4(&f, xo, "HP");
1449        let lp = response::lr4(&f, xo, "LP");
1450        for (pos, h) in positions {
1451            check()?;
1452            let sub = &alignment
1453                .sub_sum_by_pos
1454                .iter()
1455                .find(|(p, _)| p == pos)
1456                .expect("complete sub position")
1457                .1;
1458            let main: Vec<_> = idx
1459                .iter()
1460                .enumerate()
1461                .map(|(k, &i)| h[i] * C::from_polar(1.0, ramp_phase(f[k], d) * 0.001) * hp[k])
1462                .collect();
1463            let sub: Vec<_> = idx
1464                .iter()
1465                .enumerate()
1466                .map(|(k, &i)| sub[i] * lp[k])
1467                .collect();
1468            prepared.push((f.clone(), main, sub));
1469        }
1470    }
1471    if prepared.is_empty() {
1472        return Ok((0.0, 1));
1473    }
1474    let mut best = (1e18, 0.0, 1);
1475    for pol in [1, -1] {
1476        let mut k = 0;
1477        while lo_ms + k as f64 * 0.25 < 15.01 {
1478            check()?;
1479            let d = lo_ms + k as f64 * 0.25;
1480            k += 1;
1481            let mut cost = 0.0;
1482            for (f, main, sub) in &prepared {
1483                let errors: Vec<_> = f
1484                    .iter()
1485                    .zip(main)
1486                    .zip(sub)
1487                    .map(|((&f, &m), &s)| {
1488                        let s = s
1489                            * C::new(f64::from(pol), 0.0)
1490                            * C::from_polar(1.0, ramp_phase(f, d) * 0.001);
1491                        let total = power(m + s);
1492                        let reference = power(m) + power(s);
1493                        (10.0 * (reference + 1e-30).log10() - 10.0 * (total + 1e-30).log10())
1494                            .powi(2)
1495                    })
1496                    .collect();
1497                cost += mean(&errors);
1498            }
1499            observe(pol, d, cost);
1500            if cost < best.0 {
1501                best = (cost, d, pol);
1502            }
1503        }
1504    }
1505    Ok((best.1, best.2))
1506}
1507pub fn fold(rows: &[BaRow], bus_gain: f64, bus_delay_ms: f64, bus_polarity: i8) -> Vec<BaRow> {
1508    rows.iter()
1509        .map(|r| {
1510            row(
1511                r.gain_db + 20.0 * bus_gain.log10(),
1512                r.delay_ms + bus_delay_ms,
1513                r.polarity * bus_polarity,
1514                r.allpass.clone(),
1515            )
1516        })
1517        .collect()
1518}
1519pub fn sub_mode_label(role: &str, setup: &Setup, layout_info: &[(String, LayoutInfo)]) -> String {
1520    let Some((_, m)) = setup.sub_manual.iter().find(|(r, _)| r == role) else {
1521        return "auto".into();
1522    };
1523    match m.mode {
1524        SubMode::Auto => "auto".into(),
1525        SubMode::Fixed => "fixed".into(),
1526        SubMode::Bundle => {
1527            let name = if m.bundle.is_empty() { "A" } else { &m.bundle };
1528            if let Some((_, info)) = layout_info.iter().find(|(n, _)| n == name) {
1529                let kind = match info.kind {
1530                    LayoutKind::Multi => "multi",
1531                    LayoutKind::Dba => {
1532                        if info.rear.iter().any(|r| r == role) {
1533                            "dba-rear"
1534                        } else {
1535                            "dba-front"
1536                        }
1537                    }
1538                };
1539                format!("bundle:{name}:{kind}")
1540            } else {
1541                format!("bundle:{name}")
1542            }
1543        }
1544    }
1545}
1546
1547#[cfg(test)]
1548#[path = "tests/ba.rs"]
1549mod tests;