use crate::alloc::borrow::ToOwned;
use core::fmt::Debug;
use core::{
fmt::Display,
ops::{Deref, DerefMut},
};
use alloc::string::{String, ToString};
use alloc::vec::Vec;
use errors::{DimensionsSource, ScheduleDecodingError};
use ndarray::{Array, Axis, Dimension, IntoDimension, Ix1};
use self::point_spread::PointSpread;
pub mod point_spread;
/// 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.
#[derive(Clone, PartialEq, Eq)]
pub struct Schedule<Dim: Dimension>(Array<bool, Dim>);
impl<Dim: Dimension> Debug for Schedule<Dim> {
fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
Display::fmt(self, f)
}
}
impl<Dim: Dimension> Display for Schedule<Dim> {
fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
if self.ndim() == 0 {
f.write_str("Zero dimensional schedule")?;
return Ok(());
}
let dim = self.shape();
let mut index = alloc::vec![0; self.ndim()];
for value in &**self {
f.write_str(if *value { "▩" } else { "_" })?;
index[0] += 1;
for i in 0..index.len() {
if index[i] == dim[i] {
// f.write_str("\n")?;
index[i] = 0;
if i + 1 < self.ndim() {
index[i + 1] += 1;
}
}
}
}
Ok(())
}
}
impl<Dim: Dimension> Deref for Schedule<Dim> {
type Target = Array<bool, Dim>;
fn deref(&self) -> &Self::Target {
&self.0
}
}
impl<Dim: Dimension> DerefMut for Schedule<Dim> {
fn deref_mut(&mut self) -> &mut Self::Target {
&mut self.0
}
}
/// When saving or loading a schedule, would the first sample be represented by zero or one?
#[derive(Clone, Copy, Debug)]
pub enum EncodingType {
/// The first sample is `0`
ZeroBased,
/// The first sample is `1`
OneBased,
}
impl<Dim: Dimension> Schedule<Dim> {
/// Create a new schedule from a list of booleans
pub const fn new(sched: Array<bool, Dim>) -> Schedule<Dim> {
Schedule(sched)
}
/// Unwrap the inner list of bools
pub fn into_inner(self) -> Array<bool, Dim> {
self.0
}
/// Count the number of samples in the schedule. Note that this is operation linear in the length of the schedule.
pub fn count(&self) -> usize {
self.0.iter().filter(|v| **v).count()
}
/// 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.
///
/// `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:
/// ```txt
/// 0
/// 100000000000000000000
/// ```
///
/// # Example
///
/// ```
/// # use nmr_schedule::*;
/// # use ndarray::*;
/// let sched = Schedule::new(Array::from_vec(vec![true, false, true, true, false, true, false, false]));
///
/// assert_eq!(
/// Schedule::<Ix1>::decode(
/// "0\n2\n3\n5",
/// EncodingType::ZeroBased,
/// |max| Ok(Ix1(max[0].next_power_of_two()))
/// ).unwrap(),
/// sched
/// );
///
/// assert!(
/// Schedule::<Ix1>::decode(
/// "1\n3\n4\n6",
/// EncodingType::OneBased,
/// |max| if max[0] > 4 { Err("Too big!".to_owned()) } else { Ok(max) },
/// ).is_err(),
/// );
///
/// let sched = Schedule::new(Array::from_shape_vec(IxDyn(&[8, 2]).strides(IxDyn(&[1, 8])), vec![
/// false, false, true, false, false, false, false, false,
/// true, false, false, true, false, true, false, false,
/// ]).unwrap());
///
/// assert_eq!(
/// Schedule::<IxDyn>::decode(
/// "1, 2\n3, 1\n4, 2\n6, 2",
/// EncodingType::OneBased,
/// |max| Ok(IxDyn(&max.as_array_view().iter().map(|v| v.next_power_of_two()).collect::<Vec<_>>())),
/// ).unwrap(),
/// sched,
/// );
///
/// let sched = Schedule::new(Array::from_vec(vec![true, false, true, true, false, true, false]));
/// assert_eq!(
/// Schedule::<Ix1>::decode(
/// "0\n2\n3\n5",
/// EncodingType::ZeroBased,
/// |_| Ok(Ix1(7)),
/// ).unwrap(),
/// sched,
/// );
/// ```
pub fn decode(
encoded: &str,
encoding_type: EncodingType,
dims_rule: impl FnOnce(Dim) -> Result<Dim, String>,
) -> Result<Schedule<Dim>, ScheduleDecodingError> {
let spots = encoded
.split('\n')
.filter(|s| !s.is_empty())
.map(|v| {
v.split(",")
.filter(|s| !s.is_empty())
.map(|v| {
v.trim().parse::<usize>().map_err(ScheduleDecodingError::NotANumber).and_then(|v| {
Ok(v - match encoding_type {
EncodingType::ZeroBased => 0,
EncodingType::OneBased => {
if v == 0 {
return Err(ScheduleDecodingError::ZeroSampleInOneEncodedData);
} else {
1
}
},
})
})
})
.collect::<Result<Vec<usize>, _>>()
})
.collect::<Result<Vec<Vec<usize>>, _>>()?;
let mut dim_amt = None;
for (i, spot) in spots.iter().enumerate() {
match dim_amt {
None => dim_amt = Some(spot.len()),
Some(existing_dims) => {
if existing_dims != spot.len() {
return Err(ScheduleDecodingError::DifferentDimensions {
expected: existing_dims,
line: i + 1,
found: spot.len(),
});
}
}
}
}
let dim_amt = dim_amt.unwrap_or(Dim::NDIM.unwrap_or(0));
if let Some(expected_dims) = Dim::NDIM {
if expected_dims != dim_amt {
return Err(ScheduleDecodingError::InconsistentDimensions {
found: dim_amt,
expected: expected_dims,
why: DimensionsSource::Ty(core::any::type_name::<Dim>()),
});
}
}
let mut min_dims = Dim::zeros(dim_amt);
for spot in spots.iter() {
for i in 0..dim_amt {
min_dims[i] = min_dims[i].max(spot[i] + 1);
}
}
let dims = match dims_rule(min_dims.to_owned()) {
Ok(v) => v,
Err(msg) => return Err(ScheduleDecodingError::CouldntFindDims(msg)),
};
if dims.ndim() != dim_amt && dim_amt != 0 {
return Err(ScheduleDecodingError::InconsistentDimensions {
found: dim_amt,
expected: dims.ndim(),
why: DimensionsSource::Rule,
});
}
for i in 0..dim_amt {
if dims[i] < min_dims[i] {
return Err(ScheduleDecodingError::DataTooBig {
given: dims.into_dyn(),
minimum: min_dims.into_dyn(),
});
}
}
let mut sched = Schedule(Array::from_elem(dims, false));
for spot in spots {
let mut dim = Dim::zeros(dim_amt);
for i in 0..dim_amt {
dim[i] = spot[i];
}
sched.0[dim] = true;
}
Ok(sched)
}
/// Writes the schedule to a string using the same format as is understood by [`Schedule::decode`].
///
/// # Example
///
/// ```
/// # use nmr_schedule::*;
/// # use ndarray::*;
/// let sched = Schedule::new(Array::from_vec(vec![true, false, false, true, true]));
/// assert_eq!(sched.encode(EncodingType::ZeroBased), "0\n3\n4");
/// assert_eq!(sched.encode(EncodingType::OneBased), "1\n4\n5");
/// ```
///
/// ```
/// # use nmr_schedule::*;
/// # use ndarray::*;
/// assert_eq!(
/// Schedule::decode(
/// "1, 2\n3, 4",
/// EncodingType::OneBased,
/// |v: Ix2| Ok(v),
/// )
/// .unwrap()
/// .encode(EncodingType::OneBased),
/// "1, 2\n3, 4",
/// );
/// ```
pub fn encode(&self, encoding_type: EncodingType) -> String {
self.indexed_iter()
.filter(|v| *v.1)
.map(|v| {
let mut spot = v.0.into_dimension();
spot.as_array_view_mut().iter_mut().for_each(|v| *v += match encoding_type {
EncodingType::ZeroBased => 0,
EncodingType::OneBased => 1,
});
spot.as_array_view().iter().map(|v| v.to_string()).collect::<Vec<_>>().join(", ")
})
.collect::<Vec<_>>()
.join("\n")
}
}
impl Schedule<Ix1> {
/// Calculates the Point Spread Function of the schedule.
///
/// The PSF is calculated as the discrete fourier transform of the schedule where `true` is replaced with `1.` and `false` is replaced with `0.`.
///
/// 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.
///
/// 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.
pub fn point_spread(&self) -> PointSpread {
PointSpread::calculate(self)
}
/// Calculate the Repeat Length Curve of the schedule for a given pattern.
///
/// 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.
///
/// 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`.
///
/// # Example
///
/// ```
/// # use nmr_schedule::*;
/// # use ndarray::*;
/// let sched = Schedule::new(Array::from_vec(vec![true, false, true, false, true, false, false, true, false, true]));
///
/// // The schedule has 3 instances of `T F T` and 1 instance of `T F T F T`
/// assert_eq!(sched.rlc([true, false]), vec![3, 1]);
/// ```
///
/// L. E. Cullen, A. Marchiori, D. Rovnyak, Magn Reson Chem 2023, 61(6), 337. <https://doi.org/10.1002/mrc.5340>
pub fn rlc(&self, seq: impl AsRef<[bool]>) -> Vec<u64> {
let mut counts = Vec::new();
let seq = seq.as_ref();
while counts.last() != Some(&0) {
counts.push(
self.axis_windows(Axis(0), seq.len() * (1 + counts.len()) + 1)
.into_iter()
.filter(|subsequence| {
subsequence
.iter()
.enumerate()
.all(|(i, v)| seq[i % seq.len()] == *v)
})
.count() as u64,
)
}
counts.pop();
counts
}
}
/// Errors that can be generated by the crate
pub mod errors {
use alloc::string::String;
use core::{fmt::Display, num::ParseIntError};
use ndarray::IxDyn;
/// The source of an amount of dimensions specified
#[derive(Clone, Copy, Debug)]
pub enum DimensionsSource {
/// The type defines a specific number of dimensions (two for `Ix2`)
Ty(&'static str),
/// The rule function passed in gave a specific number of dimensions
Rule,
}
/// An error that can occur while decoding a string as a schedule
#[derive(Clone, Debug)]
pub enum ScheduleDecodingError {
/// A value couldn't be decoded as a number
NotANumber(ParseIntError),
/// A sample was zero when the schedule was expected to be encoded with the first sample being one
ZeroSampleInOneEncodedData,
/// The dimensions expected of the data don't match the dimensions expected by the callee
InconsistentDimensions {
/// The dimensionality of the data
found: usize,
/// The dimensionality expected
expected: usize,
/// The reason the dimensionality was expected
why: DimensionsSource,
},
/// The length callback couldn't find a length
CouldntFindDims(String),
/// The data is too big to fit in the schedule
DataTooBig {
/// The dimensions specified by the rule function
given: IxDyn,
/// The minimum dimensionality of the schedule
minimum: IxDyn,
},
/// The schedule contains lines of different dimensions
DifferentDimensions {
/// The dimensions of previous lines
expected: usize,
/// The line number of the inconsistent line
line: usize,
/// The dimensionality of the line
found: usize,
},
}
impl Display for ScheduleDecodingError {
fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
match self {
ScheduleDecodingError::NotANumber(parse_error) => Display::fmt(parse_error, f),
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?"),
ScheduleDecodingError::InconsistentDimensions { found, expected, why } => {
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 "))?;
match why {
DimensionsSource::Ty(name) => f.write_fmt(format_args!("the choice of {name} as the dimensionality of the schedule.")),
DimensionsSource::Rule => f.write_fmt(format_args!("the definition of `dims_rule`.")),
}
},
ScheduleDecodingError::CouldntFindDims(message) => f.write_fmt(format_args!("Couldn't find a length for the schedule: {message}")),
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:?}.")),
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}")),
}
}
}
impl core::error::Error for ScheduleDecodingError {}
}
#[cfg(test)]
mod tests {
use ndarray::{Ix1, Ix2, IxDyn};
use crate::{ errors::{DimensionsSource, ScheduleDecodingError}, EncodingType, Schedule};
use alloc::borrow::ToOwned;
#[test]
fn decode_errors() {
assert!(matches!(
Schedule::<IxDyn>::decode("3\n2\noof", EncodingType::ZeroBased, Ok),
Err(ScheduleDecodingError::NotANumber(_))
));
assert!(matches!(
Schedule::<IxDyn>::decode("0\n1\n2", EncodingType::OneBased, Ok),
Err(ScheduleDecodingError::ZeroSampleInOneEncodedData)
));
assert!(matches!(
Schedule::<Ix1>::decode("1\n3\n500", EncodingType::OneBased, |max| if max[0] > 8 {
Err("Too big!".to_owned())
} else {
Ok(max)
}),
Err(ScheduleDecodingError::CouldntFindDims(_))
));
assert!(matches!(
Schedule::decode("2, 4\n5, 2\n245, 9", EncodingType::OneBased, |_| Ok(Ix2(
8, 8
))),
Err(ScheduleDecodingError::DataTooBig {
given: _,
minimum: _
})
));
assert!(matches!(
Schedule::decode("1, 2\n3, 4\n5, 6\n", EncodingType::ZeroBased, |_| Ok(Ix1(
8
))),
Err(ScheduleDecodingError::InconsistentDimensions {
found: 2,
expected: 1,
why: DimensionsSource::Ty(_)
}),
));
assert!(matches!(
Schedule::decode("1, 2\n3, 4\n5, 6\n", EncodingType::ZeroBased, |_| Ok(
IxDyn(&[8, 8, 8])
)),
Err(ScheduleDecodingError::InconsistentDimensions {
found: 2,
expected: 3,
why: DimensionsSource::Rule,
}),
));
}
}