Skip to main content

sva_samples/collapse/
truncate.rs

1// Concern: truncates every series to the terms the profile leaves, once per collapse | Non-concern: evaluating what is left (point.rs) | IO: (&SpectralSum or &Body) -> the same, series-free
2
3use std::f64::consts::TAU;
4
5use sva_formula::affine::{Axis, axis, exact_constant};
6use sva_formula::closed_form::{Bound, Series, children, map_children};
7use sva_formula::series::{mentions_line, ratio, substitute};
8use sva_formula::spectral_sum::atom::{Exp, Factors, Singular, SpectralAtom};
9use sva_formula::table::series::{Shape, read};
10use sva_formula::{
11    Body, C64, Codomain, Env, Lane, NodeId, ParamId, Part, Run, SpectralSum, Ty, Unary, Var, lines,
12};
13
14use crate::error::CollapseError;
15use crate::profile::Profile;
16
17/// What a series is truncated against: the observation's own ceiling, the amplitude a
18/// dropped tail provably stays within, and the floor a tail no bound sums falls under.
19#[derive(Clone, Copy, Debug, PartialEq)]
20pub struct Audible {
21    ceiling: f64,
22    precision: f64,
23    floor_db: f64,
24}
25
26impl Audible {
27    pub fn of(profile: &Profile, rate: u32) -> Audible {
28        Audible::on(profile, crate::Grid::of(rate))
29    }
30
31    pub fn on(profile: &Profile, grid: crate::Grid) -> Audible {
32        let ceiling = profile.ceiling_hz.min(grid.sr() / 2.0);
33        Audible {
34            ceiling,
35            precision: profile.half_lsb(),
36            floor_db: profile.floor(ceiling),
37        }
38    }
39}
40
41/// A line series is placed analytically under `sva_formula`'s own bound; an expanded one
42/// becomes a formula the sample loop walks, which is what this caps.
43const MAX_EXPANDED_TERMS: usize = 1 << 13;
44
45/// FORMAT 6.2 truncates a series once, at collapse: every term the ceiling and the precision
46/// leave becomes an ordinary atom before the first sample is read.
47pub fn spectral_sum(n: &SpectralSum, band: Audible) -> Result<SpectralSum, CollapseError> {
48    if n.lanes.iter().all(|l| l.series.is_empty()) {
49        return Ok(n.clone());
50    }
51    let mut lanes = Vec::with_capacity(n.lanes.len());
52    for lane in &n.lanes {
53        let mut held = Lane {
54            series: Vec::new(),
55            ..lane.clone()
56        };
57        for s in &lane.series {
58            held.atoms.extend(atoms(s, band)?);
59        }
60        lanes.push(held);
61    }
62    Ok(SpectralSum::of(n.var, lanes))
63}
64
65/// The same truncation over a written closed form, whose series a spectral sum never reached.
66pub fn written(f: &Body, band: Audible) -> Result<Body, CollapseError> {
67    if let Body::Series(s) = f {
68        return expanded(s, band);
69    }
70    let mut found = None;
71    let out = map_children(f, |p| match written(&p.body, band) {
72        Ok(body) => Part::new(p.origin, body),
73        Err(e) => {
74            found = Some(e);
75            p.clone()
76        }
77    });
78    match found {
79        Some(e) => Err(e),
80        None => Ok(out),
81    }
82}
83
84fn atoms(s: &Series, band: Audible) -> Result<Vec<SpectralAtom>, CollapseError> {
85    let (body, window) = sva_formula::crop_peeled(&s.term.body);
86    let bare = Series {
87        term: Part::new(s.term.origin, body),
88        ..s.clone()
89    };
90    let taken = enumerated(&bare, band);
91    if !taken.is_empty() {
92        return Ok(taken
93            .into_iter()
94            .map(|l| {
95                SpectralAtom::new(
96                    l.amp,
97                    Factors {
98                        exp: Some(Exp::at(0.0, TAU * l.hz)),
99                        ind: window,
100                        ..Factors::NONE
101                    },
102                    Singular::Regular,
103                    s.term.origin,
104                )
105            })
106            .collect());
107    }
108    let body = expanded(s, band)?;
109    let held = sva_formula::normalize(&body, Var::T)
110        .map_err(|_| CollapseError::NotEvaluable("a series term"))?;
111    Ok(held.lanes.into_iter().flat_map(|l| l.atoms).collect())
112}
113
114/// A line series keeps its terms as runs; anything else is summed term by term.
115fn expanded(s: &Series, band: Audible) -> Result<Body, CollapseError> {
116    let taken = enumerated(s, band);
117    if !taken.is_empty() {
118        let runs = Run::of(&taken).into_iter();
119        return Ok(sum(runs.map(|r| Body::Run(Box::new(r))).collect()));
120    }
121    let count = terms(s, band)?;
122    let mut parts = Vec::with_capacity(count);
123    for i in 0..count {
124        let form = substitute(&s.term.body, s.index, (s.lo + i as i64) as f64);
125        parts.push(written(&form, band)?);
126    }
127    Ok(sum(parts))
128}
129
130fn enumerated(s: &Series, band: Audible) -> Vec<sva_formula::Line> {
131    match read(&s.term.body) {
132        Some(Shape::Lines(_)) => lines(s, band.ceiling, band.floor_db, band.precision).taken,
133        _ => Vec::new(),
134    }
135}
136
137fn sum(parts: Vec<Body>) -> Body {
138    match parts.len() {
139        0 => Body::Const(C64::ZERO),
140        1 => parts.into_iter().next().expect("one part"),
141        _ => Body::Add(parts.into_iter().map(Part::bare).collect()),
142    }
143}
144
145/// A nesting is priced whole before it expands: a written bound counts like the floor's.
146fn terms(s: &Series, band: Audible) -> Result<usize, CollapseError> {
147    let cost = cost(s, band).ok_or(CollapseError::NotEvaluable(
148        "a series whose term count no coefficient bounds",
149    ))?;
150    if cost.depth == 1 && cost.own > MAX_EXPANDED_TERMS {
151        return Err(CollapseError::NotEvaluable(
152            "a series of more terms than one instant expands",
153        ));
154    }
155    match cost.terms > MAX_EXPANDED_TERMS {
156        true => Err(CollapseError::NestedSeries {
157            depth: cost.depth,
158            terms: cost.terms,
159            bound: MAX_EXPANDED_TERMS,
160        }),
161        false => Ok(cost.own),
162    }
163}
164
165/// Series deep, terms the whole nesting takes, terms this level alone takes.
166struct Cost {
167    depth: usize,
168    terms: usize,
169    own: usize,
170}
171
172fn cost(s: &Series, band: Audible) -> Option<Cost> {
173    let own = counted(s, band)?;
174    let inside = substitute(&s.term.body, s.index, s.lo as f64);
175    let (depth, inner) = price(&inside, band)?;
176    Some(Cost {
177        depth: depth + 1,
178        terms: own.saturating_mul(inner),
179        own,
180    })
181}
182
183/// A sum of series costs their counts added; every other spelling multiplies.
184fn price(f: &Body, band: Audible) -> Option<(usize, usize)> {
185    if let Body::Series(s) = f {
186        let cost = cost(s, band)?;
187        return Some((cost.depth, cost.terms));
188    }
189    let summed = matches!(f, Body::Add(_));
190    let mut depth = 0;
191    let mut terms: usize = if summed { 0 } else { 1 };
192    for p in children(f) {
193        let (d, n) = price(&p.body, band)?;
194        depth = depth.max(d);
195        terms = match summed {
196            true => terms.saturating_add(n),
197            false => terms.saturating_mul(n),
198        };
199    }
200    Some((depth, terms.max(1)))
201}
202
203/// A geometric magnitude bound stops where its whole tail rounds away; any other stand-in at
204/// the profile's floor.
205fn counted(s: &Series, band: Audible) -> Option<usize> {
206    let hi = match s.hi {
207        Bound::Finite(n) => return usize::try_from((n - s.lo + 1).max(0)).ok(),
208        Bound::Infinite => MAX_EXPANDED_TERMS,
209    };
210    let bound = bound(&s.term.body, band)?;
211    let ratio = match bound.exact {
212        true => ratio(&bound.body, s.index).filter(|r| *r < 1.0),
213        false => None,
214    };
215    let floor = 10f64.powf(band.floor_db / 20.0);
216    let mut peak = 0.0f64;
217    for i in 0..hi {
218        let at = substitute(&bound.body, s.index, (s.lo + i as i64) as f64);
219        let held = exact_constant(&at)?.abs();
220        peak = peak.max(held);
221        let gone = match ratio {
222            Some(r) => held / (1.0 - r) <= band.precision,
223            None => peak > 0.0 && held < peak * floor,
224        };
225        if gone {
226            return Some(i.max(1));
227        }
228    }
229    None
230}
231
232/// A term's coefficient in the index. `exact` where it bounds the term's magnitude; elsewhere
233/// it only stands in, a turning factor or a name read as one.
234struct Bounded {
235    body: Body,
236    exact: bool,
237}
238
239fn bound(f: &Body, band: Audible) -> Option<Bounded> {
240    let held = |body: Body, exact: bool| Some(Bounded { body, exact });
241    let under = |p: &Part| bound(&p.body, band).map(|b| (Part::new(p.origin, b.body), b.exact));
242    let all = |parts: &[Part], wrap: &dyn Fn(Part) -> Part| {
243        let mut exact = true;
244        let mut kept = Vec::with_capacity(parts.len());
245        for p in parts {
246            let (part, known) = under(p)?;
247            exact &= known;
248            kept.push(wrap(part));
249        }
250        Some((kept, exact))
251    };
252    match f {
253        Body::Series(s) => match total(s, band) {
254            Some(sum) => held(Body::Const(C64::real(sum)), true),
255            None => held(Body::Const(C64::ONE), false),
256        },
257        Body::Line | Body::Node(_) => held(Body::Const(C64::ONE), false),
258        Body::Mul(parts) => {
259            let (kept, exact) = all(parts, &|p| p)?;
260            held(Body::Mul(kept), exact)
261        }
262        Body::Add(parts) => {
263            let (kept, exact) = all(parts, &|p| Part::bare(Body::Apply(Unary::Abs, p)))?;
264            held(Body::Add(kept), exact)
265        }
266        Body::Div(a, b) if !mentions_line(&b.body) => {
267            let (num, exact) = under(a)?;
268            held(Body::Div(num, b.clone()), exact)
269        }
270        Body::Div(a, b) => {
271            let ((num, _), (den, _)) = (under(a)?, under(b)?);
272            held(Body::Div(num, den), false)
273        }
274        Body::Apply(Unary::Sin | Unary::Cos | Unary::Tanh | Unary::Sat | Unary::Step, arg) => {
275            held(Body::Const(C64::ONE), real(&arg.body))
276        }
277        Body::Crop { of, .. } | Body::Shift { of, .. } | Body::Warp { of, .. } => {
278            bound(&of.body, band)
279        }
280        Body::Pow(base, n) => {
281            let (part, exact) = under(base)?;
282            held(Body::Pow(part, *n), exact && *n >= 0)
283        }
284        other if !mentions_line(other) => held(other.clone(), true),
285        _ => None,
286    }
287}
288
289/// A geometric series sums to its first bound over `1 - ratio`; any other to what it keeps.
290fn total(s: &Series, band: Audible) -> Option<f64> {
291    let bound = bound(&s.term.body, band).filter(|b| b.exact)?.body;
292    let at = |i: i64| Some(exact_constant(&substitute(&bound, s.index, i as f64))?.abs());
293    if let (Bound::Infinite, Some(r)) = (s.hi, ratio(&bound, s.index).filter(|r| *r < 1.0)) {
294        return Some(at(s.lo)? / (1.0 - r));
295    }
296    let kept = i64::try_from(counted(s, band).filter(|n| *n <= MAX_EXPANDED_TERMS)?).ok()?;
297    (s.lo..s.lo + kept).try_fold(0.0, |held, i| Some(held + at(i)?))
298}
299
300fn real(f: &Body) -> bool {
301    axis(f, &Unread) == Axis::Real
302}
303
304struct Unread;
305
306impl Env for Unread {
307    fn node(&self, _: NodeId) -> Ty {
308        Ty::form(Var::T, false, Codomain::Complex)
309    }
310
311    fn param(&self, _: ParamId) -> Ty {
312        Ty::form(Var::T, false, Codomain::Complex)
313    }
314}