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

sva_samples/collapse/
truncate.rs
