nmr_schedule/pdf/
presets.rs

1use alloc::{borrow::ToOwned, vec::Vec};
2
3use crate::{InitialState, Monotonicity, find_zero};
4
5use super::Pdf;
6
7/// An unweighted PDF.
8///
9/// ```
10/// # use nmr_schedule::pdf::*;
11/// assert!(
12///     unweighted(256)
13///         .get_distribution()
14///         .iter()
15///         .all(|v| *v == 1. / 256.)
16/// );
17/// ```
18pub fn unweighted(len: usize) -> Pdf {
19    Pdf::from_integral(move |pos| pos / len as f64, len)
20}
21
22/// A linear PDF
23///
24/// ```
25/// # use nmr_schedule::pdf::*;
26/// assert!(
27///     linear(256)
28///         .get_distribution()
29///         .windows(2)
30///         .all(|v| v[0] - v[1] == 1. / 256. / 128.)
31/// );
32/// ```
33pub fn linear(len: usize) -> Pdf {
34    let inv_len = 1. / len as f64;
35    Pdf::from_integral(move |pos| inv_len * pos * (2. - inv_len * pos), len)
36}
37
38/// An exponential PDF. `scale` represents how quickly the PDF decays.
39///
40/// For QSched compatibility, `scale` can be set to `bias·time`
41///
42/// J.C.J Barna, E.D Laue, M.R Mayger, J Skilling, S.J.P Worrall, Exponential sampling, an alternative method for sampling in two-dimensional NMR experiments, Journal of Magnetic Resonance (1969), Volume 73, Issue 1, 1987, Pages 69-77, ISSN 0022-2364, https://doi.org/10.1016/0022-2364(87)90225-3.
43pub fn exponential(len: usize, scale: f64) -> Pdf {
44    let c = -scale / len as f64;
45    let normalize = 1. / ((-scale).exp() - 1.);
46
47    Pdf::from_integral(move |pos| normalize * (c * pos).exp(), len)
48}
49
50/// Options for the `bias` parameter used in QSched
51#[derive(Clone, Copy, Debug)]
52pub enum QSinBias {
53    /// Equivalent to setting 'bias' to '1'.
54    /// PDF is { 1 - sin(x) | 0 < x < time/2 }
55    Low,
56    /// Equivalent to setting 'bias' to '1.5'.
57    /// PDF is { (1 - sin(x))^1.8 | 0 < x < time/2 }
58    Med,
59    /// Equivalent to setting 'bias' to '2'.
60    /// PDF is { (0.85 - sin(x))^2 | 0 < x < time/2 }
61    High,
62}
63
64impl core::fmt::Display for QSinBias {
65    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
66        f.write_str(match self {
67            QSinBias::Low => "Low",
68            QSinBias::Med => "Med",
69            QSinBias::High => "High",
70        })
71    }
72}
73
74/// Equivalent to the `quant-sin` option in QSched
75///
76/// Palmer, M.R., Wenrich, B.R., Stahlfeld, P. et al. Performance tuning non-uniform sampling for sensitivity enhancement of signal-limited biological NMR. J Biomol NMR 58, 303–314 (2014). https://doi.org/10.1007/s10858-014-9823-5
77pub fn qsin(len: usize, bias: QSinBias, time: f64) -> Pdf {
78    let half_time = time / 2.;
79    let len_f = len as f64;
80
81    let s = half_time / len_f;
82
83    match bias {
84        QSinBias::Low => {
85            let f = |x: f64| x + x.cos();
86
87            let normalize = (f(half_time) - f(0.)).recip();
88
89            Pdf::from_integral(move |pos| f(pos * s) * normalize, len)
90        }
91        QSinBias::Med => {
92            let bars = 16;
93            let bars_f = bars as f64;
94
95            let mut discrete = alloc::vec![0.; len];
96
97            let mut sum = 0.;
98
99            for (i, v) in discrete.iter_mut().enumerate() {
100                let mut approximation = 0.;
101
102                for bar in 0..=bars {
103                    let x = i as f64 + bar as f64 / bars_f;
104                    approximation += (1. - (x * s).sin()).powf(1.8)
105                        * if bar == 0 || bar == bars { 0.5 } else { 1. };
106                }
107
108                *v = approximation;
109                sum += *v;
110            }
111
112            let scale = sum.recip();
113
114            discrete.iter_mut().for_each(|v| *v *= scale);
115
116            Pdf::from_discrete(discrete)
117        }
118        QSinBias::High => {
119            let f = |x: f64| 1.2225 * x + 1.7 * x.cos() - 0.25 * (2. * x).sin();
120
121            let normalize = (f(half_time) - f(0.)).recip();
122
123            Pdf::from_integral(move |pos| f(pos * s) * normalize, len)
124        }
125    }
126}
127
128/// Estimate the effective PDF of Sin Weighted Poisson Gap, if you would like to try it in other generators.
129///
130/// This has not been published and may be incorrect.
131#[allow(clippy::missing_panics_doc)] // Panicking is a bug
132pub fn sin_weighted_poisson_gap(len: usize, count: usize) -> Pdf {
133    // An expected gap length of 2 corresponds to a probability of (1/(2+1)) or 1/3:
134    // ▩__▩__▩__
135    // We want the sum of probabilities over the length to equal the count, so we can do a binary search for the correct value. Contributions welcome if you have a closed form for this.
136
137    let (_, mk_probabilities) = find_zero(
138        Monotonicity::Decreasing,
139        InitialState::new(len as f64 / count as f64, 0., f64::INFINITY),
140        crate::Precision::Image {
141            amount: 0.001,
142            debug: (&|_, _| -> () { panic!("Convergence failed") }) as &_,
143        },
144        |arbitrary| {
145            let probabilities = (0..len).map(move |i| {
146                1. / (arbitrary
147                    * ((i as f64 + 0.5) / (len as f64 + 1.) * core::f64::consts::FRAC_PI_2).sin()
148                    + 1.)
149            });
150
151            (
152                probabilities.to_owned().sum::<f64>() - count as f64,
153                move || {
154                    probabilities.collect() // Return a function to collect instead of allocating a vector every time (which would be slow)
155                },
156            )
157        },
158    );
159
160    let probabilities: Vec<f64> = mk_probabilities();
161    let probabilities_ref = &probabilities;
162
163    // Binary search the count parameter to turn the probabilities into a PDF
164    // Again, contributions welcome if anyone knows a closed form for this.
165    let (_, mk_pdf) = find_zero(
166        Monotonicity::Increasing,
167        InitialState::new(count as f64, 0., f64::INFINITY),
168        crate::Precision::Image {
169            amount: 0.00001,
170            debug: (&|_, _| -> () { panic!("Convergence failed") }) as &_,
171        },
172        |inverse_count| {
173            // This is actually the inverse operation of 1 - (1 - v)^(count)
174            // Not bothering to take 1/count, the binary search will do that
175            let pdf_gen = move |i| {
176                let complement: f64 = 1. - probabilities_ref[i];
177                1. - complement.powf(inverse_count)
178            };
179
180            ((0..len).map(pdf_gen).sum::<f64>() - 1., move || {
181                Pdf::from_discrete((0..len).map(pdf_gen).collect())
182            })
183        },
184    );
185
186    mk_pdf()
187}
nmr_schedule/pdf/presets.rs

nmr_schedule/pdf/
presets.rs
