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