nmr_schedule/schedule/mod.rs
1use crate::alloc::borrow::ToOwned;
2use core::fmt::Debug;
3use core::{
4 fmt::Display,
5 ops::{Deref, DerefMut},
6};
7
8use alloc::string::{String, ToString};
9use alloc::vec::Vec;
10use errors::{DimensionsSource, ScheduleDecodingError};
11use ndarray::{Array, Axis, Dimension, IntoDimension, Ix1};
12
13use self::point_spread::PointSpread;
14
15pub mod point_spread;
16
17/// Represents a sampling schedule. Schedules are represented as a list of booleans, where `true` means to take the sample at the position of the index and `false` means to not take the sample.
18#[derive(Clone, PartialEq, Eq)]
19pub struct Schedule<Dim: Dimension>(Array<bool, Dim>);
20
21impl<Dim: Dimension> Debug for Schedule<Dim> {
22 fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
23 Display::fmt(self, f)
24 }
25}
26
27impl<Dim: Dimension> Display for Schedule<Dim> {
28 fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
29 if self.ndim() == 0 {
30 f.write_str("Zero dimensional schedule")?;
31 return Ok(());
32 }
33
34 let dim = self.shape();
35 let mut index = alloc::vec![0; self.ndim()];
36
37 for value in &**self {
38 f.write_str(if *value { "▩" } else { "_" })?;
39
40 index[0] += 1;
41
42 for i in 0..index.len() {
43 if index[i] == dim[i] {
44 // f.write_str("\n")?;
45 index[i] = 0;
46
47 if i + 1 < self.ndim() {
48 index[i + 1] += 1;
49 }
50 }
51 }
52 }
53
54 Ok(())
55 }
56}
57
58impl<Dim: Dimension> Deref for Schedule<Dim> {
59 type Target = Array<bool, Dim>;
60
61 fn deref(&self) -> &Self::Target {
62 &self.0
63 }
64}
65
66impl<Dim: Dimension> DerefMut for Schedule<Dim> {
67 fn deref_mut(&mut self) -> &mut Self::Target {
68 &mut self.0
69 }
70}
71
72/// When saving or loading a schedule, would the first sample be represented by zero or one?
73#[derive(Clone, Copy, Debug)]
74pub enum EncodingType {
75 /// The first sample is `0`
76 ZeroBased,
77 /// The first sample is `1`
78 OneBased,
79}
80
81impl<Dim: Dimension> Schedule<Dim> {
82 /// Create a new schedule from a list of booleans
83 pub const fn new(sched: Array<bool, Dim>) -> Schedule<Dim> {
84 Schedule(sched)
85 }
86
87 /// Unwrap the inner list of bools
88 pub fn into_inner(self) -> Array<bool, Dim> {
89 self.0
90 }
91
92 /// Count the number of samples in the schedule. Note that this is operation linear in the length of the schedule.
93 pub fn count(&self) -> usize {
94 self.0.iter().filter(|v| **v).count()
95 }
96
97 /// Decodes a schedule from a string. The format is expected to be a list of line-break separated sample positions, where each sample position is a comma separated list of numbers for each index in each dimension.
98 ///
99 /// `len_rule` is a function that takes the minimum possible length to contain all of the values and returns either a length for the schedule or an error message that will be propagated. This is necessary because it is not known how many `false` values there there should be after the last `true` value, and it allows rejecting input that could cause resource exhaustion of the software. It would make sense to configure this function to return an error above a maximum value because it is possible to make a "schedule bomb" that would consume the entire RAM when decoded:
100 /// ```txt
101 /// 0
102 /// 100000000000000000000
103 /// ```
104 ///
105 /// # Example
106 ///
107 /// ```
108 /// # use nmr_schedule::*;
109 /// # use ndarray::*;
110 /// let sched = Schedule::new(Array::from_vec(vec![true, false, true, true, false, true, false, false]));
111 ///
112 /// assert_eq!(
113 /// Schedule::<Ix1>::decode(
114 /// "0\n2\n3\n5",
115 /// EncodingType::ZeroBased,
116 /// |max| Ok(Ix1(max[0].next_power_of_two()))
117 /// ).unwrap(),
118 /// sched
119 /// );
120 ///
121 /// assert!(
122 /// Schedule::<Ix1>::decode(
123 /// "1\n3\n4\n6",
124 /// EncodingType::OneBased,
125 /// |max| if max[0] > 4 { Err("Too big!".to_owned()) } else { Ok(max) },
126 /// ).is_err(),
127 /// );
128 ///
129 /// let sched = Schedule::new(Array::from_shape_vec(IxDyn(&[8, 2]).strides(IxDyn(&[1, 8])), vec![
130 /// false, false, true, false, false, false, false, false,
131 /// true, false, false, true, false, true, false, false,
132 /// ]).unwrap());
133 ///
134 /// assert_eq!(
135 /// Schedule::<IxDyn>::decode(
136 /// "1, 2\n3, 1\n4, 2\n6, 2",
137 /// EncodingType::OneBased,
138 /// |max| Ok(IxDyn(&max.as_array_view().iter().map(|v| v.next_power_of_two()).collect::<Vec<_>>())),
139 /// ).unwrap(),
140 /// sched,
141 /// );
142 ///
143 /// let sched = Schedule::new(Array::from_vec(vec![true, false, true, true, false, true, false]));
144 /// assert_eq!(
145 /// Schedule::<Ix1>::decode(
146 /// "0\n2\n3\n5",
147 /// EncodingType::ZeroBased,
148 /// |_| Ok(Ix1(7)),
149 /// ).unwrap(),
150 /// sched,
151 /// );
152 /// ```
153 pub fn decode(
154 encoded: &str,
155 encoding_type: EncodingType,
156 dims_rule: impl FnOnce(Dim) -> Result<Dim, String>,
157 ) -> Result<Schedule<Dim>, ScheduleDecodingError> {
158 let spots = encoded
159 .split('\n')
160 .filter(|s| !s.is_empty())
161 .map(|v| {
162 v.split(",")
163 .filter(|s| !s.is_empty())
164 .map(|v| {
165 v.trim().parse::<usize>().map_err(ScheduleDecodingError::NotANumber).and_then(|v| {
166 Ok(v - match encoding_type {
167 EncodingType::ZeroBased => 0,
168 EncodingType::OneBased => {
169 if v == 0 {
170 return Err(ScheduleDecodingError::ZeroSampleInOneEncodedData);
171 } else {
172 1
173 }
174 },
175 })
176 })
177 })
178 .collect::<Result<Vec<usize>, _>>()
179 })
180 .collect::<Result<Vec<Vec<usize>>, _>>()?;
181
182 let mut dim_amt = None;
183
184 for (i, spot) in spots.iter().enumerate() {
185 match dim_amt {
186 None => dim_amt = Some(spot.len()),
187 Some(existing_dims) => {
188 if existing_dims != spot.len() {
189 return Err(ScheduleDecodingError::DifferentDimensions {
190 expected: existing_dims,
191 line: i + 1,
192 found: spot.len(),
193 });
194 }
195 }
196 }
197 }
198
199 let dim_amt = dim_amt.unwrap_or(Dim::NDIM.unwrap_or(0));
200
201 if let Some(expected_dims) = Dim::NDIM {
202 if expected_dims != dim_amt {
203 return Err(ScheduleDecodingError::InconsistentDimensions {
204 found: dim_amt,
205 expected: expected_dims,
206 why: DimensionsSource::Ty(core::any::type_name::<Dim>()),
207 });
208 }
209 }
210
211 let mut min_dims = Dim::zeros(dim_amt);
212
213 for spot in spots.iter() {
214 for i in 0..dim_amt {
215 min_dims[i] = min_dims[i].max(spot[i] + 1);
216 }
217 }
218
219 let dims = match dims_rule(min_dims.to_owned()) {
220 Ok(v) => v,
221 Err(msg) => return Err(ScheduleDecodingError::CouldntFindDims(msg)),
222 };
223
224 if dims.ndim() != dim_amt && dim_amt != 0 {
225 return Err(ScheduleDecodingError::InconsistentDimensions {
226 found: dim_amt,
227 expected: dims.ndim(),
228 why: DimensionsSource::Rule,
229 });
230 }
231
232 for i in 0..dim_amt {
233 if dims[i] < min_dims[i] {
234 return Err(ScheduleDecodingError::DataTooBig {
235 given: dims.into_dyn(),
236 minimum: min_dims.into_dyn(),
237 });
238 }
239 }
240
241 let mut sched = Schedule(Array::from_elem(dims, false));
242
243 for spot in spots {
244 let mut dim = Dim::zeros(dim_amt);
245
246 for i in 0..dim_amt {
247 dim[i] = spot[i];
248 }
249
250 sched.0[dim] = true;
251 }
252
253 Ok(sched)
254 }
255
256 /// Writes the schedule to a string using the same format as is understood by [`Schedule::decode`].
257 ///
258 /// # Example
259 ///
260 /// ```
261 /// # use nmr_schedule::*;
262 /// # use ndarray::*;
263 /// let sched = Schedule::new(Array::from_vec(vec![true, false, false, true, true]));
264 /// assert_eq!(sched.encode(EncodingType::ZeroBased), "0\n3\n4");
265 /// assert_eq!(sched.encode(EncodingType::OneBased), "1\n4\n5");
266 /// ```
267 ///
268 /// ```
269 /// # use nmr_schedule::*;
270 /// # use ndarray::*;
271 /// assert_eq!(
272 /// Schedule::decode(
273 /// "1, 2\n3, 4",
274 /// EncodingType::OneBased,
275 /// |v: Ix2| Ok(v),
276 /// )
277 /// .unwrap()
278 /// .encode(EncodingType::OneBased),
279 /// "1, 2\n3, 4",
280 /// );
281 /// ```
282 pub fn encode(&self, encoding_type: EncodingType) -> String {
283 self.indexed_iter()
284 .filter(|v| *v.1)
285 .map(|v| {
286 let mut spot = v.0.into_dimension();
287 spot.as_array_view_mut().iter_mut().for_each(|v| *v += match encoding_type {
288 EncodingType::ZeroBased => 0,
289 EncodingType::OneBased => 1,
290 });
291
292 spot.as_array_view().iter().map(|v| v.to_string()).collect::<Vec<_>>().join(", ")
293 })
294 .collect::<Vec<_>>()
295 .join("\n")
296 }
297}
298
299impl Schedule<Ix1> {
300 /// Calculates the Point Spread Function of the schedule.
301 ///
302 /// The PSF is calculated as the discrete fourier transform of the schedule where `true` is replaced with `1.` and `false` is replaced with `0.`.
303 ///
304 /// If you take uniformly sampled data and set every sample not taken to zero, the [Convolution Theorem](https://en.wikipedia.org/wiki/Convolution_theorem) states that your true signal will be convolved with the schedule's PSF. Therefore, it is preferable to generate schedules where the PSF has minimal spikes because that minimizes sampling noise.
305 ///
306 /// The PSF can also be thought of as a plot of global bias towards sampling a particular frequency. For example, a PSF spike at wavenumber 1/3 indicates a bias towards sampling at positions that are at multiples of three, or some offset of a multiple of three, thereby undersampling frequencies that are out of phase with this bias. Therefore, the PSF can be thought of as a plot of how much a schedule undersamples different phases of a particular frequency.
307 pub fn point_spread(&self) -> PointSpread {
308 PointSpread::calculate(self)
309 }
310
311 /// Calculate the Repeat Length Curve of the schedule for a given pattern.
312 ///
313 /// This metric can be used to determine the number of repeating patterns in the sampling schedule. Schedules with fewer patterns are more condusive to signal reconstruction because patterns represent a local bias against sampling particular frequencies.
314 ///
315 /// The first argument is the pattern to be searched for. The function returns a list where the element at index `i` is number of places in the schedule where `seq` is repeated `i + 1` times followed by the first element of `seq`.
316 ///
317 /// # Example
318 ///
319 /// ```
320 /// # use nmr_schedule::*;
321 /// # use ndarray::*;
322 /// let sched = Schedule::new(Array::from_vec(vec![true, false, true, false, true, false, false, true, false, true]));
323 ///
324 /// // The schedule has 3 instances of `T F T` and 1 instance of `T F T F T`
325 /// assert_eq!(sched.rlc([true, false]), vec![3, 1]);
326 /// ```
327 ///
328 /// L. E. Cullen, A. Marchiori, D. Rovnyak, Magn Reson Chem 2023, 61(6), 337. <https://doi.org/10.1002/mrc.5340>
329 pub fn rlc(&self, seq: impl AsRef<[bool]>) -> Vec<u64> {
330 let mut counts = Vec::new();
331 let seq = seq.as_ref();
332
333 while counts.last() != Some(&0) {
334 counts.push(
335 self.axis_windows(Axis(0), seq.len() * (1 + counts.len()) + 1)
336 .into_iter()
337 .filter(|subsequence| {
338 subsequence
339 .iter()
340 .enumerate()
341 .all(|(i, v)| seq[i % seq.len()] == *v)
342 })
343 .count() as u64,
344 )
345 }
346
347 counts.pop();
348
349 counts
350 }
351}
352
353/// Errors that can be generated by the crate
354pub mod errors {
355 use alloc::string::String;
356 use core::{fmt::Display, num::ParseIntError};
357
358 use ndarray::IxDyn;
359
360 /// The source of an amount of dimensions specified
361 #[derive(Clone, Copy, Debug)]
362 pub enum DimensionsSource {
363 /// The type defines a specific number of dimensions (two for `Ix2`)
364 Ty(&'static str),
365 /// The rule function passed in gave a specific number of dimensions
366 Rule,
367 }
368
369 /// An error that can occur while decoding a string as a schedule
370 #[derive(Clone, Debug)]
371 pub enum ScheduleDecodingError {
372 /// A value couldn't be decoded as a number
373 NotANumber(ParseIntError),
374 /// A sample was zero when the schedule was expected to be encoded with the first sample being one
375 ZeroSampleInOneEncodedData,
376 /// The dimensions expected of the data don't match the dimensions expected by the callee
377 InconsistentDimensions {
378 /// The dimensionality of the data
379 found: usize,
380 /// The dimensionality expected
381 expected: usize,
382 /// The reason the dimensionality was expected
383 why: DimensionsSource,
384 },
385 /// The length callback couldn't find a length
386 CouldntFindDims(String),
387 /// The data is too big to fit in the schedule
388 DataTooBig {
389 /// The dimensions specified by the rule function
390 given: IxDyn,
391 /// The minimum dimensionality of the schedule
392 minimum: IxDyn,
393 },
394 /// The schedule contains lines of different dimensions
395 DifferentDimensions {
396 /// The dimensions of previous lines
397 expected: usize,
398 /// The line number of the inconsistent line
399 line: usize,
400 /// The dimensionality of the line
401 found: usize,
402 },
403 }
404
405 impl Display for ScheduleDecodingError {
406 fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
407 match self {
408 ScheduleDecodingError::NotANumber(parse_error) => Display::fmt(parse_error, f),
409 ScheduleDecodingError::ZeroSampleInOneEncodedData => f.write_str("A sample was zero when the schedule was expected to be encoded as one-based. Is this schedule really encoded as zero-based?"),
410 ScheduleDecodingError::InconsistentDimensions { found, expected, why } => {
411 f.write_fmt(format_args!("The dimensionality expected doesn't match the dimensionality of the data. The data is {found} dimensional whereas the schedule is expected to be {expected} dimensional. The dimensionality expected is from "))?;
412
413 match why {
414 DimensionsSource::Ty(name) => f.write_fmt(format_args!("the choice of {name} as the dimensionality of the schedule.")),
415 DimensionsSource::Rule => f.write_fmt(format_args!("the definition of `dims_rule`.")),
416
417 }
418 },
419 ScheduleDecodingError::CouldntFindDims(message) => f.write_fmt(format_args!("Couldn't find a length for the schedule: {message}")),
420 ScheduleDecodingError::DataTooBig { given: found, minimum } => f.write_fmt(format_args!("The schedule is too big to fit in the dimensions given. The data defines a minimum schedule of {minimum:?} but the dimensions given are {found:?}.")),
421 ScheduleDecodingError::DifferentDimensions { expected, line, found } => f.write_fmt(format_args!("The data contains sample positions of different dimensionalities. The dimensionality found so far is {expected} whereas the dimensionallity on line {line} is {found}")),
422 }
423 }
424 }
425
426 impl core::error::Error for ScheduleDecodingError {}
427}
428
429#[cfg(test)]
430mod tests {
431 use ndarray::{Ix1, Ix2, IxDyn};
432
433 use crate::{ errors::{DimensionsSource, ScheduleDecodingError}, EncodingType, Schedule};
434 use alloc::borrow::ToOwned;
435
436 #[test]
437 fn decode_errors() {
438
439 assert!(matches!(
440 Schedule::<IxDyn>::decode("3\n2\noof", EncodingType::ZeroBased, Ok),
441 Err(ScheduleDecodingError::NotANumber(_))
442 ));
443
444 assert!(matches!(
445 Schedule::<IxDyn>::decode("0\n1\n2", EncodingType::OneBased, Ok),
446 Err(ScheduleDecodingError::ZeroSampleInOneEncodedData)
447 ));
448
449 assert!(matches!(
450 Schedule::<Ix1>::decode("1\n3\n500", EncodingType::OneBased, |max| if max[0] > 8 {
451 Err("Too big!".to_owned())
452 } else {
453 Ok(max)
454 }),
455 Err(ScheduleDecodingError::CouldntFindDims(_))
456 ));
457
458 assert!(matches!(
459 Schedule::decode("2, 4\n5, 2\n245, 9", EncodingType::OneBased, |_| Ok(Ix2(
460 8, 8
461 ))),
462 Err(ScheduleDecodingError::DataTooBig {
463 given: _,
464 minimum: _
465 })
466 ));
467
468 assert!(matches!(
469 Schedule::decode("1, 2\n3, 4\n5, 6\n", EncodingType::ZeroBased, |_| Ok(Ix1(
470 8
471 ))),
472 Err(ScheduleDecodingError::InconsistentDimensions {
473 found: 2,
474 expected: 1,
475 why: DimensionsSource::Ty(_)
476 }),
477 ));
478
479 assert!(matches!(
480 Schedule::decode("1, 2\n3, 4\n5, 6\n", EncodingType::ZeroBased, |_| Ok(
481 IxDyn(&[8, 8, 8])
482 )),
483 Err(ScheduleDecodingError::InconsistentDimensions {
484 found: 2,
485 expected: 3,
486 why: DimensionsSource::Rule,
487 }),
488 ));
489 }
490}