imgal 0.3.0

use ndarray::{Array1, Array3, ArrayBase, AsArray, Ix1, ViewRepr, Zip};

use crate::filter::fft_convolve_1d;
use crate::prelude::*;
use crate::simulation::instrument;
use crate::statistics::sum;

/// Create a 1D Gaussian IRF convolved monoexponential or multiexponential decay
/// curve.
///
/// # Description
///
/// Creates a 1D Gaussian instrument response function (IRF) convolved with an
/// ideal monoexponential or multiexponential decay curve defined as the sum of
/// one or more exponential components, each characterized by a lifetime (τ) and
/// fractional intensity:
///
/// ```text
/// I(t) = [Σᵢ αᵢ × exp(-t/τᵢ)] ⊗ IRF(t)
/// ```
///
/// # Arguments
///
/// * `samples`: The number of discrete points that make up the decay curve.
/// * `period`: The period (*i.e.* time interval).
/// * `taus`: An array of lifetimes. For a monoexponential decay curve use a
///   single tau value and a fractional intensity of `1.0`. For a
///   multiexponential decay curve use two or more tau values, matched with
///   their respective fractional intensity. The `taus` and `fractions` arrays
///   must have the same length. Tau values set to `0.0` will be skipped.
/// * `fractions`: An array of fractional intensities for each tau in the `taus`
///   array. The `fractions` array must be the same length as the `taus` array
///   and sum to `1.0`. Fraction values set to `0.0` will be skipped.
/// * `total_counts`: The total intensity count (*e.g.* photon count) of the
///   decay curve.
/// * `irf_center`: The temporal position of the IRF peak within the time range.
/// * `irf_width`: The full width at half maximum (FWHM) of the IRF.
/// * `threads`: The requested number of threads to use for parallel execution.
///   If `None` or `Some(1)` sequential execution is used. If `Some(0)`, then
///   the maximum available parallelism is used. Thread counts are clamped to
///   the systems maximum.
///
/// # Returns
///
/// * `Ok(Vec<f64>)`: The 1D Gaussian IRF convolved monoexponential or
///   multiexponential decay curve.
/// * `Err(ImgalError)`: If `taus.len() != fractions.len()`. If
///   fractions array does not sum to `1.0`.
#[inline]
pub fn gaussian_exponential_decay_1d<'a, A>(
    samples: usize,
    period: f64,
    taus: A,
    fractions: A,
    total_counts: f64,
    irf_center: f64,
    irf_width: f64,
    threads: Option<usize>,
) -> Result<Array1<f64>, ImgalError>
where
    A: AsArray<'a, f64, Ix1>,
{
    let irf = instrument::gaussian_irf_1d(samples, period, irf_center, irf_width, threads);
    let i_arr =
        ideal_exponential_decay_1d(samples, period, taus, fractions, total_counts, threads)?;
    Ok(fft_convolve_1d(&i_arr, &irf, threads))
}

/// Create a 3D Gaussian IRF convolved monoexponential or multiexponential decay
/// curve.
///
/// # Description
///
/// Creates a 3D Gaussian instrument response function (IRF) convolved with an
/// ideal monoexponential or multiexponential decay curve defined as the sum of
/// one or more exponential components, each characterized by a lifetime (τ) and
/// fractional intensity:
///
/// ```text
/// I(t) = [Σᵢ αᵢ × exp(-t/τᵢ)] ⊗ IRF(t)
/// ```
///
/// # Arguments
///
/// * `samples`: The number of discrete points that make up the decay curve.
/// * `period`: The period (*i.e.* time interval).
/// * `taus`: An array of lifetimes. For a monoexponential decay curve use a
///   single tau value and a fractional intensity of `1.0`. For a
///   multiexponential decay curve use two or more tau values, matched with
///   their respective fractional intensity. The `taus` and `fractions` arrays
///   must have the same length. Tau values set to `0.0` will be skipped.
/// * `fractions`: An array of fractional intensities for each tau in the `taus`
///   array. The `fractions` array must be the same length as the `taus` array
///   and sum to `1.0`. Fraction values set to `0.0` will be skipped.
/// * `total_counts`: The total intensity count (*e.g.* photon count) of the
///   decay curve.
/// * `irf_center`: The temporal position of the IRF peak within the time range.
/// * `irf_width`: The full width at half maximum (FWHM) of the IRF.
/// * `shape`: The row and col shape to broadcast the decay curve into.
/// * `threads`: The requested number of threads to use for parallel execution.
///   If `None` or `Some(1)` sequential execution is used. If `Some(0)`, then
///   the maximum available parallelism is used. Thread counts are clamped to
///   the systems maximum.
///
/// # Returns
///
/// * `Ok(Array3<f64>)`: The 3D Gaussian IRF convolved monoexponential or
///   multiexponential decay curve with dimensions (row, col, t).
/// * `Err(ImgalError)`: If `taus.len() != fractions.len()`. If
///   fractions array does not sum to `1.0`.
#[inline]
pub fn gaussian_exponential_decay_3d<'a, A>(
    samples: usize,
    period: f64,
    taus: A,
    fractions: A,
    total_counts: f64,
    irf_center: f64,
    irf_width: f64,
    shape: (usize, usize),
    threads: Option<usize>,
) -> Result<Array3<f64>, ImgalError>
where
    A: AsArray<'a, f64, Ix1>,
{
    let i_arr = gaussian_exponential_decay_1d(
        samples,
        period,
        taus,
        fractions,
        total_counts,
        irf_center,
        irf_width,
        threads,
    )?;
    let dims = (shape.0, shape.1, samples);
    Ok(i_arr.broadcast(dims).unwrap().to_owned())
}

/// Create a 1D ideal monoexponential or multiexponential decay curve.
///
/// # Description
///
/// Creates a 1D ideal exponential decay curve by computing the sum of one or
/// more exponential components, each characterized by a lifetime (τ) and
/// fractional intensity as defined by:
///
/// ```text
/// I(t) = Σᵢ αᵢ × exp(-t/τᵢ)
/// ```
///
/// Where `αᵢ` are the pre-exponential factors derived from the fractional
/// intensities and lifetimes.
///
/// # Arguments
///
/// * `samples`: The number of discrete points that make up the decay curve.
/// * `period`: The period (*i.e.* time interval).
/// * `taus`: An array of lifetimes. For a monoexponential decay curve use a
///   single tau value and a fractional intensity of `1.0`. For a
///   multiexponential decay curve use two or more tau values, matched with
///   their respective fractional intensity. The `taus` and `fractions` arrays
///   must have the same length. Tau values set to `0.0` will be skipped.
/// * `fractions`: An array of fractional intensities for each tau in the `taus`
///   array. The `fractions` array must be the same length as the `taus` array
///   and sum to `1.0`. Fraction values set to `0.0` will be skipped.
/// * `total_counts`: The total intensity count (*e.g.* photon count) of the
///   decay curve.
/// * `threads`: The requested number of threads to use for parallel execution.
///   If `None` or `Some(1)` sequential execution is used. If `Some(0)`, then
///   the maximum available parallelism is used. Thread counts are clamped to
///   the systems maximum.
///
/// # Returns
///
/// * `Ok(Array1<f64>)`: The 1D monoexponential or multiexponential decay curve.
/// * `Err(ImgalError)`: If `taus.len() != fractions.len()`. If
///   fractions array does not sum to `1.0`.
///
/// # Reference
///
/// <https://doi.org/10.1111/j.1749-6632.1969.tb56231.x>
#[inline]
pub fn ideal_exponential_decay_1d<'a, A>(
    samples: usize,
    period: f64,
    taus: A,
    fractions: A,
    total_counts: f64,
    threads: Option<usize>,
) -> Result<Array1<f64>, ImgalError>
where
    A: AsArray<'a, f64, Ix1>,
{
    let taus: ArrayBase<ViewRepr<&'a f64>, Ix1> = taus.into();
    let fractions: ArrayBase<ViewRepr<&'a f64>, Ix1> = fractions.into();
    let tl = taus.len();
    let fl = fractions.len();
    if tl != fl {
        return Err(ImgalError::MismatchedArrayLengths {
            a_arr_name: "taus",
            a_arr_len: tl,
            b_arr_name: "fractions",
            b_arr_len: fl,
        });
    }
    let fs = sum(&fractions, threads);
    if fs != 1.0 {
        return Err(ImgalError::InvalidSum {
            expected: 1.0,
            got: fs,
        });
    }
    // compute the pre-exponential factors (alpha) and construct the decay curve
    // scaled to the total counts
    let alph_arr = &fractions / &taus;
    let mut i_arr = vec![0.0; samples];
    let time_arr = Array1::linspace(0.0, period, samples);
    alph_arr
        .iter()
        .zip(taus.iter())
        .filter(|&(&al, &ta)| al != 0.0 && ta != 0.0)
        .for_each(|(al, ta)| {
            Zip::from(&mut i_arr).and(&time_arr).for_each(|i, t| {
                *i += al * (-t / ta).exp();
            });
        });
    let scale = total_counts / sum(&i_arr, threads);
    i_arr.iter_mut().for_each(|v| *v *= scale);
    Ok(Array1::from_vec(i_arr))
}

/// Create a 3D ideal monoexponential or multiexponential decay curve.
///
/// # Description
///
/// Creates a 3D ideal exponential decay curve by computing the sum of one or
/// more exponential components, each characterized by a lifetime (τ) and
/// fractional intensity as defined by:
///
/// ```text
/// I(t) = Σᵢ αᵢ × exp(-t/τᵢ)
/// ```
///
/// Where `αᵢ` are the pre-exponential factors derived from the fractional
/// intensities and lifetimes.
///
/// # Arguments
///
/// * `samples`: The number of discrete points that make up the decay curve.
/// * `period`: The period (*i.e.* time interval).
/// * `taus`: An array of lifetimes. For a monoexponential decay curve use a
///   single tau value and a fractional intensity of `1.0`. For a
///   multiexponential decay curve use two or more tau values, matched with
///   their respective fractional intensity. The `taus` and `fractions` arrays
///   must have the same length. Tau values set to `0.0` will be skipped.
/// * `fractions`: An array of fractional intensities for each tau in the `taus`
///   array. The `fractions` array must be the same length as the `taus` array
///   and sum to `1.0`. Fraction values set to `0.0` will be skipped.
/// * `total_counts`: The total intensity count (*e.g.* photon count) of the
///   decay curve.
/// * `shape`: The row and col shape to broadcast the decay curve into.
/// * `threads`: The requested number of threads to use for parallel execution.
///   If `None` or `Some(1)` sequential execution is used. If `Some(0)`, then
///   the maximum available parallelism is used. Thread counts are clamped to
///   the systems maximum.
///
/// # Returns
///
/// * `Ok(Array3<f64>)`: The 3D monoexponential or multiexponential decay curve
///   with dimensions (row, col, t).
/// * `Err(ImgalError)`: If `taus.len() != fractions.len()`. If
///   fractions array does not sum to `1.0`.
///
/// # Reference
///
/// <https://doi.org/10.1111/j.1749-6632.1969.tb56231.x>
#[inline]
pub fn ideal_exponential_decay_3d<'a, A>(
    samples: usize,
    period: f64,
    taus: A,
    fractions: A,
    total_counts: f64,
    shape: (usize, usize),
    threads: Option<usize>,
) -> Result<Array3<f64>, ImgalError>
where
    A: AsArray<'a, f64, Ix1>,
{
    let i_arr =
        ideal_exponential_decay_1d(samples, period, taus, fractions, total_counts, threads)?;
    let dims = (shape.0, shape.1, samples);
    Ok(i_arr.broadcast(dims).unwrap().to_owned())
}

/// Create a 1D IRF convolved monoexponential or multiexponential decay curve.
///
/// # Description
///
/// Creates a 1D instrument response function (IRF) convolved with an ideal
/// monoexponential or multiexponential decay curve defined as the sum of one or
/// more exponential components, each characterized by a lifetime (τ) and
/// fractional intensity:
///
/// ```text
/// I(t) = [Σᵢ αᵢ × exp(-t/τᵢ)] ⊗ IRF(t)
/// ```
///
/// # Arguments
///
/// * `irf`: The IRF as a 1D array.
/// * `samples`: The number of discrete points that make up the decay curve.
/// * `period`: The period (*i.e.* time interval).
/// * `taus`: An array of lifetimes. For a monoexponential decay curve use a
///   single tau value and a fractional intensity of `1.0`. For a
///   multiexponential decay curve use two or more tau values, matched with
///   their respective fractional intensity. The `taus` and `fractions` arrays
///   must have the same length. Tau values set to `0.0` will be skipped.
/// * `fractions`: An array of fractional intensities for each tau in the `taus`
///   array. The `fractions` array must be the same length as the `taus` array
///   and sum to `1.0`. Fraction values set to `0.0` will be skipped.
/// * `total_counts`: The total intensity count (*e.g.* photon count) of the
///   decay curve.
/// * `threads`: The requested number of threads to use for parallel execution.
///   If `None` or `Some(1)` sequential execution is used. If `Some(0)`, then
///   the maximum available parallelism is used. Thread counts are clamped to
///   the systems maximum.
///
/// # Returns
///
/// * `Ok(Array1<f64>)`: The 1D IRF convolved monoexponential or
///   multiexponential decay curve.
/// * `Err(ImgalError)`: If `taus.len() != fractions.len()`. If fractions array
///   does not sum to `1.0`.
#[inline]
pub fn irf_exponential_decay_1d<'a, A, B>(
    irf: A,
    samples: usize,
    period: f64,
    taus: B,
    fractions: B,
    total_counts: f64,
    threads: Option<usize>,
) -> Result<Array1<f64>, ImgalError>
where
    A: AsArray<'a, f64, Ix1>,
    B: AsArray<'a, f64, Ix1>,
{
    let irf: ArrayBase<ViewRepr<&'a f64>, Ix1> = irf.into();
    let i_arr =
        ideal_exponential_decay_1d(samples, period, taus, fractions, total_counts, threads)?;
    Ok(fft_convolve_1d(i_arr.view(), irf, threads))
}

/// Create a 3D IRF convolved monoexponential or multiexponential decay curve.
///
/// # Description
///
/// Creates a 3D instrument response function (IRF) convolved with an ideal
/// monoexponential or multiexponential decay curve defined as the sum of one or
/// more exponential components, each characterized by a lifetime (τ) and
/// fractional intensity:
///
/// ```text
/// I(t) = [Σᵢ αᵢ × exp(-t/τᵢ)] ⊗ IRF(t)
/// ```
///
/// # Arguments
///
/// * `irf`: The IRF as a 1D array.
/// * `samples`: The number of discrete points that make up the decay curve.
/// * `period`: The period (*i.e.* time interval).
/// * `taus`: An array of lifetimes. For a monoexponential decay curve use a
///   single tau value and a fractional intensity of `1.0`. For a
///   multiexponential decay curve use two or more tau values, matched with
///   their respective fractional intensity. The `taus` and `fractions` arrays
///   must have the same length. Tau values set to `0.0` will be skipped.
/// * `fractions`: An array of fractional intensities for each tau in the `taus`
///   array. The `fractions` array must be the same length as the `taus` array
///   and sum to `1.0`. Fraction values set to `0.0` will be skipped.
/// * `total_counts`: The total intensity count (*e.g.* photon count) of the
///   decay curve.
/// * `shape`: The row and col shape to broadcast the decay curve into.
/// * `threads`: The requested number of threads to use for parallel execution.
///   If `None` or `Some(1)` sequential execution is used. If `Some(0)`, then
///   the maximum available parallelism is used. Thread counts are clamped to
///   the systems maximum.
///
/// # Returns
///
/// * `Ok(Array3<f64>)`: The 3D IRF convolved monoexponential or
///   multiexponential decay curve with dimensions (row, col, t).
/// * `Err(ImgalError)`: If `taus.len() != fractions.len()`. If fractions array
///   does not sum to `1.0`.
#[inline]
pub fn irf_exponential_decay_3d<'a, A, B>(
    irf: A,
    samples: usize,
    period: f64,
    taus: B,
    fractions: B,
    total_counts: f64,
    shape: (usize, usize),
    threads: Option<usize>,
) -> Result<Array3<f64>, ImgalError>
where
    A: AsArray<'a, f64, Ix1>,
    B: AsArray<'a, f64, Ix1>,
{
    let i_arr =
        irf_exponential_decay_1d(irf, samples, period, taus, fractions, total_counts, threads)?;
    let dims = (shape.0, shape.1, samples);
    Ok(i_arr.broadcast(dims).unwrap().to_owned())
}