prodef 0.1.0

//! A module that implements a multivariate normal probability density function.

use crate::{Density, Domain};
use itertools::{Itertools, zip_eq};
use nalgebra::{
    Const, DMatrix, DVector, DefaultAllocator, Dim, Dyn, MatrixView, OMatrix, OVector, RealField,
    Scalar, U1, VectorView, allocator::Allocator,
};
use rand::Rng;
use rand_distr::{Distribution, StandardNormal};
use serde::{Deserialize, Serialize};
use std::{
    cmp::Ordering,
    iter::Sum,
    ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign},
};

/// A multivariate normal probability density function.
#[derive(Clone, Debug, Deserialize, Serialize)]
#[serde(bound(
    serialize = "D: Serialize, B: Serialize, OVector<T, D>: Serialize, OMatrix<T, D, D>: Serialize"
))]
#[serde(bound(
    deserialize = "D: Serialize, B: Deserialize<'de>, OVector<T, D>: Deserialize<'de>, OMatrix<T, D, D>: Deserialize<'de>"
))]
pub struct MultiNormalDensity<T, D, B>
where
    T: Scalar,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
{
    matrix: OMatrix<T, D, D>,
    inverse: OMatrix<T, D, D>,
    ltm: OMatrix<T, D, D>,
    domain: B,
    mean: OVector<T, D>,
}

impl<T, D, B> MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    B: 'static + Domain<T, D>,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
{
    /// Returns the determinant of the covariance matrix.
    pub fn determinant(&self) -> T {
        self.ltm
            .diagonal()
            .iter()
            .fold(T::one(), |acc, next| {
                if *next != T::zero() {
                    acc * next.clone()
                } else {
                    acc
                }
            })
            .powi(2)
    }

    /// Returns the amount of Fisher information for derivatves of the mean `da` and `db`.
    pub fn fisher_information<RStride: Dim, CStride: Dim>(
        &self,
        da: &VectorView<T, D, RStride, CStride>,
        db: &VectorView<T, D, RStride, CStride>,
    ) -> T {
        (da.transpose() * &self.inverse * db)[(0, 0)].clone()
    }

    /// Create a [`MultiNormalDensity`] from a covariance matrix.
    pub fn from_matrix(matrix: OMatrix<T, D, D>, mean: OVector<T, D>, domain: B) -> Option<Self> {
        let n_dim = matrix.shape_generic().0;

        let inverse = {
            // Convert input matrix into a DMatrix to bypass annoying ToTypenum domain.
            let dmatrix =
                DMatrix::from_iterator(matrix.nrows(), matrix.ncols(), matrix.iter().cloned());

            let mut pinv = dmatrix
                .clone_owned()
                .pseudo_inverse(T::default_epsilon())
                .expect("failed to compute pseudo inverse");

            dmatrix
                .diagonal()
                .iter()
                .enumerate()
                .for_each(|(idx, value)| {
                    if matches!(
                        value
                            .partial_cmp(&T::zero())
                            .expect("covariance matrix contains NaN values"),
                        Ordering::Equal
                    ) {
                        pinv.set_column(idx, &DVector::<T>::zeros(matrix.ncols()));
                        pinv.set_row(idx, &DVector::<T>::zeros(matrix.ncols()).transpose());
                    }
                });

            let n_dim = matrix.shape_generic().0;

            OMatrix::<T, D, D>::from_iterator_generic(n_dim, n_dim, pinv.iter().cloned())
        };

        let mut d = OVector::<T, D>::zeros_generic(n_dim, Const::<1>);
        let mut l = OMatrix::<T, D, D>::zeros_generic(n_dim, n_dim);

        for cdx in 0..n_dim.value() {
            let mut d_j = matrix[(cdx, cdx)].clone();

            if cdx > 0 {
                for k in 0..cdx {
                    d_j -= d[k].clone() * l[(cdx, k)].clone().powi(2);
                }
            }

            d[cdx] = d_j;

            for rdx in cdx..n_dim.value() {
                let mut l_ij = matrix[(rdx, cdx)].clone();

                for k in 0..cdx {
                    l_ij -= d[k].clone() * l[(cdx, k)].clone() * l[(rdx, k)].clone();
                }

                if matches!(
                    d[cdx]
                        .partial_cmp(&T::zero())
                        .expect("matrix contains NaN values"),
                    Ordering::Equal
                ) {
                    l[(rdx, cdx)] = T::zero();
                } else {
                    l[(rdx, cdx)] = l_ij / d[cdx].clone();
                }
            }
        }

        let lsqrtd = l * OMatrix::from_diagonal(&OVector::from_iterator_generic(
            n_dim,
            Const::<1>,
            d.iter().map(|value| value.clone().sqrt()),
        ));

        Some(Self {
            matrix,
            inverse,
            ltm: lsqrtd,
            domain,
            mean,
        })
    }

    /// Create a [`MultiNormalDensity`] from a set of vectors, with optional weights.
    pub fn from_view<RStride: Dim, CStride: Dim>(
        vectors: &MatrixView<T, D, Dyn, RStride, CStride>,
        domain: B,
        opt_weights: Option<&[T]>,
    ) -> Option<Self>
    where
        T: Sum,
    {
        let n_dim = vectors.shape_generic().0;

        // Construct the covariance matrix.
        let mut matrix = OMatrix::<T, D, D>::from_iterator_generic(
            n_dim,
            n_dim,
            (0..(vectors.nrows().pow(2))).map(|idx| {
                let jdx = idx / n_dim.value();
                let kdx = idx % n_dim.value();

                if jdx <= kdx {
                    let x = vectors.row(jdx);
                    let y = vectors.row(kdx);

                    if !x.iter().all_equal() && !y.iter().all_equal() {
                        match opt_weights {
                            Some(w) => covariance_with_weights(x, y, w),
                            None => covariance(x, y),
                        }
                    } else {
                        T::zero()
                    }
                } else {
                    T::zero()
                }
            }),
        );

        // Fill up the other side of the covariance matrix.
        matrix += matrix.transpose() - OMatrix::from_diagonal(&matrix.diagonal());

        let mut mean = vectors.column_mean();

        // Set mean to first particle value if covariance is zero.
        // This fixes numerical issues where taking the mean over many
        // particles does not equal the constant value.
        matrix
            .diagonal()
            .iter()
            .zip(mean.iter_mut())
            .enumerate()
            .for_each(|(idx, (cov, value))| {
                if cov.eq(&T::zero()) {
                    *value = vectors[(idx, 0)].clone();
                }
            });

        Self::from_matrix(matrix, mean, domain)
    }

    /// Compute the Kullback-Leibler divergence between two [`MultiNormalDensity`]'s.
    pub fn kl_div(&self, other: &MultiNormalDensity<T, D, B>) -> Option<T>
    where
        T: Sum,
    {
        let mut l_0 = self.ltm.clone();
        let mu_0 = &self.mean;

        let mut l_1 = other.ltm.clone();
        let mu_1 = &other.mean;

        let mut n_dim = self.matrix.shape_generic().0.value();

        // Detect zero'd columns/rows that need to be modified.
        (0..l_0.nrows()).for_each(|idx| {
            if l_0[(idx, idx)].eq(&T::zero()) {
                l_0[(idx, idx)] = T::one() / T::zero();

                n_dim -= 1;

                // Set off diagonals to zero.
                for jdx in 0..l_0.ncols() {
                    if jdx != idx {
                        l_0[(idx, jdx)] = T::zero();
                        l_0[(jdx, idx)] = T::zero();
                    }
                }
            };
        });

        // Detect zero'd columns/rows that need to be modified.
        (0..l_1.nrows()).for_each(|idx| {
            if l_1[(idx, idx)].eq(&T::zero()) {
                l_1[(idx, idx)] = T::one() / T::zero();

                // Set off diagonals to zero.
                for jdx in 0..l_1.ncols() {
                    if jdx != idx {
                        l_1[(idx, jdx)] = T::zero();
                        l_1[(jdx, idx)] = T::zero();
                    }
                }
            };
        });

        let mut m = l_1.clone().solve_lower_triangular(&l_0).unwrap();

        // Detect NaN's and zero them out.
        m.iter_mut().for_each(|value| {
            if !value.is_finite() {
                *value = T::zero()
            }
        });

        let y = l_1.clone().solve_lower_triangular(&(mu_1 - mu_0)).unwrap();

        Some(
            (m.iter().cloned().sum::<T>() - T::from_usize(n_dim).unwrap()
                + y.norm()
                + T::from_usize(2).unwrap()
                    * l_1
                        .diagonal()
                        .iter()
                        .zip(l_0.diagonal().iter())
                        .map(|(a, b)| {
                            if a.is_finite() && b.is_finite() {
                                (a.clone() / b.clone()).ln()
                            } else {
                                T::zero()
                            }
                        })
                        .sum::<T>())
                / T::from_usize(2).unwrap(),
        )
    }

    /// Returns the log density value for a given sample.
    pub fn log_density<RStride: Dim, CStride: Dim>(
        &self,
        sample: &VectorView<T, D, RStride, CStride>,
    ) -> Option<T> {
        if !self.domain.contains(sample) {
            return None;
        }

        let diff = sample - self.mean.clone();
        let value = (diff.transpose() * &self.inverse * diff)[(0, 0)].clone();

        let p_nonzero: usize = self.matrix.diagonal().iter().fold(0, |acc, next| {
            if next.abs() > T::zero() { acc + 1 } else { acc }
        });

        Some(
            -(self.determinant().ln()
                + value
                + T::from_usize(p_nonzero).unwrap() * T::two_pi().ln())
                / T::from_usize(2).unwrap(),
        )
    }

    /// Returns the squared Mahalanobis distance.
    pub fn mahalanobis_distance_sq<RStride: Dim, CStride: Dim>(
        &self,
        x: &VectorView<T, D, RStride, CStride>,
    ) -> T {
        (&(x - &self.mean).transpose() * &self.inverse * (x - &self.mean))[(0, 0)].clone()
    }

    /// Returns the normalization factor for the multivariate distribution.
    pub fn normalization_factor(&self) -> T {
        T::one()
            / (T::two_pi()).powf(T::from_usize(self.rank()).unwrap() / T::from_usize(2).unwrap())
            / self.determinant().sqrt()
    }

    /// Returns the rank of the underlying covariance matrix.
    pub fn rank(&self) -> usize {
        self.ltm
            .diagonal()
            .fold(0, |acc, next| if next != T::zero() { acc + 1 } else { acc })
    }

    /// Sets the mean of the distribution.
    pub fn set_mean(&mut self, mean: OVector<T, D>) {
        self.mean = mean;
    }
}

impl<T, D, B> Density<T, D> for MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    B: 'static + Domain<T, D>,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
    StandardNormal: Distribution<T>,
{
    fn density<RStride: Dim, CStride: Dim>(
        &self,
        sample: &VectorView<T, D, RStride, CStride>,
    ) -> Option<T> {
        (&self).density(sample)
    }

    fn domain(&self) -> impl Domain<T, D> + 'static {
        self.domain.clone()
    }

    fn center(&self) -> OVector<T, D> {
        (&self).center()
    }

    fn is_constant(&self) -> OVector<bool, D> {
        (&self).is_constant()
    }

    fn sample(&self, rng: &mut impl Rng, max_attempts: usize) -> Option<OVector<T, D>> {
        (&self).sample(rng, max_attempts)
    }
}

impl<T, D, B> Density<T, D> for &MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    B: 'static + Domain<T, D>,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
    StandardNormal: Distribution<T>,
{
    fn density<RStride: Dim, CStride: Dim>(
        &self,
        sample: &VectorView<T, D, RStride, CStride>,
    ) -> Option<T> {
        if !self.domain.contains(sample) {
            return None;
        }

        let diff = sample - self.mean.clone();
        let value = (diff.transpose() * &self.inverse * diff)[(0, 0)].clone();

        let p_nonzero: usize = self.matrix.diagonal().iter().fold(0, |acc, next| {
            if next.abs() > T::zero() { acc + 1 } else { acc }
        });

        Some(
            (-value / T::from_usize(2).unwrap()).exp()
                / ((T::two_pi()).powi(p_nonzero as i32) * self.determinant()).sqrt(),
        )
    }

    fn domain(&self) -> impl Domain<T, D> + 'static {
        self.domain.clone()
    }

    fn center(&self) -> OVector<T, D> {
        self.mean.clone()
    }

    fn is_constant(&self) -> OVector<bool, D> {
        OVector::from_iterator_generic(
            self.matrix.shape_generic().0,
            U1,
            self.matrix
                .diagonal()
                .iter()
                .map(|value| value.eq(&T::zero())),
        )
    }

    fn sample(&self, rng: &mut impl Rng, max_attempts: usize) -> Option<OVector<T, D>> {
        let normal = StandardNormal;
        let n_dim = self.matrix.shape_generic().0;

        let mut proposal = self.mean.clone()
            + &self.ltm
                * OVector::<T, D>::from_iterator_generic(
                    n_dim,
                    U1,
                    (0..n_dim.value()).map(|_| rng.sample(normal)),
                );

        // Counter for rejected proposals.
        let mut attempts = 0;

        while !self.domain.contains(&proposal.as_view()) {
            proposal = self.mean.clone()
                + &self.ltm
                    * OVector::<T, D>::from_iterator_generic(
                        n_dim,
                        U1,
                        (0..n_dim.value()).map(|_| rng.sample(normal)),
                    );

            attempts += 1;

            if attempts > max_attempts {
                return None;
            }
        }

        Some(proposal)
    }
}

impl<T, D, B> Add<OVector<T, D>> for MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
{
    type Output = MultiNormalDensity<T, D, B>;

    fn add(self, rhs: OVector<T, D>) -> Self::Output {
        Self {
            matrix: self.matrix,
            inverse: self.inverse,
            ltm: self.ltm,
            domain: self.domain,
            mean: self.mean + rhs,
        }
    }
}

impl<T, D, B> Add<&OVector<T, D>> for MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
{
    type Output = MultiNormalDensity<T, D, B>;

    fn add(self, rhs: &OVector<T, D>) -> Self::Output {
        Self {
            matrix: self.matrix,
            inverse: self.inverse,
            ltm: self.ltm,
            domain: self.domain,
            mean: self.mean + rhs,
        }
    }
}

impl<T, D, B> AddAssign<OVector<T, D>> for MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
{
    fn add_assign(&mut self, rhs: OVector<T, D>) {
        self.mean += rhs
    }
}

impl<T, D, B> AddAssign<&OVector<T, D>> for MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
{
    fn add_assign(&mut self, rhs: &OVector<T, D>) {
        self.mean += rhs
    }
}

impl<T, D, B> Mul<T> for MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
{
    type Output = MultiNormalDensity<T, D, B>;

    fn mul(self, rhs: T) -> Self::Output {
        Self {
            matrix: self.matrix * rhs.clone(),
            inverse: self.inverse / rhs.clone(),
            ltm: self.ltm * rhs.sqrt(),
            domain: self.domain,
            mean: self.mean,
        }
    }
}

impl<T, D, B> MulAssign<T> for MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
{
    fn mul_assign(&mut self, rhs: T) {
        self.matrix *= rhs.clone();
        self.inverse /= rhs.clone();
        self.ltm *= rhs.sqrt();
    }
}

impl<T, D, B> Sub<OVector<T, D>> for MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
{
    type Output = MultiNormalDensity<T, D, B>;

    fn sub(self, rhs: OVector<T, D>) -> Self::Output {
        Self {
            matrix: self.matrix,
            inverse: self.inverse,
            ltm: self.ltm,
            domain: self.domain,
            mean: self.mean - rhs,
        }
    }
}

impl<T, D, B> Sub<&OVector<T, D>> for MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
{
    type Output = MultiNormalDensity<T, D, B>;

    fn sub(self, rhs: &OVector<T, D>) -> Self::Output {
        Self {
            matrix: self.matrix,
            inverse: self.inverse,
            ltm: self.ltm,
            domain: self.domain,
            mean: self.mean - rhs,
        }
    }
}

impl<T, D, B> SubAssign<OVector<T, D>> for MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
{
    fn sub_assign(&mut self, rhs: OVector<T, D>) {
        self.mean -= rhs
    }
}

impl<T, D, B> SubAssign<&OVector<T, D>> for MultiNormalDensity<T, D, B>
where
    T: RealField,
    D: Dim,
    DefaultAllocator: Allocator<D> + Allocator<U1, D> + Allocator<D, D>,
{
    fn sub_assign(&mut self, rhs: &OVector<T, D>) {
        self.mean -= rhs
    }
}

/// Computes the unbiased covariance over two slices.
///
/// The length of both iterators must be equal (panic).
pub fn covariance<'a, T, I>(x: I, y: I) -> T
where
    T: RealField + Sum,
    I: IntoIterator<Item = &'a T>,
    <I as IntoIterator>::IntoIter: Clone,
{
    let x_iter = x.into_iter();
    let y_iter = y.into_iter();

    let length = x_iter.clone().fold(0, |acc, _| acc + 1);

    let mu_x = x_iter.clone().cloned().sum::<T>() / T::from_usize(length).unwrap();
    let mu_y = x_iter.clone().cloned().sum::<T>() / T::from_usize(length).unwrap();

    zip_eq(x_iter, y_iter)
        .map(|(val_x, val_y)| (mu_x.clone() - val_x.clone()) * (mu_y.clone() - val_y.clone()))
        .sum::<T>()
        / T::from_usize(length - 1).unwrap()
}

/// Computes the unbiased covariance over two slices with weights.
///
/// The length of all three iterators must be equal (panic).
pub fn covariance_with_weights<'a, T, IV, IW>(x: IV, y: IV, w: IW) -> T
where
    T: RealField + Sum,
    IV: IntoIterator<Item = &'a T>,
    IW: IntoIterator<Item = &'a T>,
    <IV as IntoIterator>::IntoIter: Clone,
    <IW as IntoIterator>::IntoIter: Clone,
{
    let x_iter = x.into_iter();
    let y_iter = y.into_iter();
    let w_iter = w.into_iter();

    let wsum = w_iter.clone().cloned().sum::<T>();
    let wsumsq = w_iter.clone().map(|val_w| val_w.clone().powi(2)).sum::<T>();

    let wfac = wsum.clone() - wsumsq / wsum.clone();

    let mu_x = zip_eq(x_iter.clone(), w_iter.clone())
        .map(|(val_x, val_w)| val_x.clone() * val_w.clone())
        .sum::<T>()
        / wsum.clone();

    let mu_y = zip_eq(y_iter.clone(), w_iter.clone())
        .map(|(val_y, val_w)| val_y.clone() * val_w.clone())
        .sum::<T>()
        / wsum.clone();

    zip_eq(x_iter, zip_eq(y_iter, w_iter))
        .map(|(val_x, (val_y, val_w))| {
            (mu_x.clone() - val_x.clone()) * (mu_y.clone() - val_y.clone()) * val_w.clone()
        })
        .sum::<T>()
        / wfac
}

#[cfg(test)]
mod tests {
    use crate::domain::{MDomain, SDomain};

    use super::*;
    use approx::ulps_eq;
    use nalgebra::{Matrix, SVector, U3, VecStorage};
    use rand::{Rng, SeedableRng};
    use rand_xoshiro::Xoshiro256PlusPlus;

    #[test]
    fn test_multinormal_density() {
        let mut rng = Xoshiro256PlusPlus::seed_from_u64(1);
        let uniform = StandardNormal;

        let mut array = Matrix::<f64, U3, Dyn, VecStorage<f64, U3, Dyn>>::from_iterator(
            10000,
            (0..30000).map(|idx| {
                if idx % 3 == 1 {
                    0.0
                } else {
                    rng.sample::<f64, StandardNormal>(uniform)
                }
            }),
        );

        array.row_mut(0).iter_mut().for_each(|value| {
            *value += 0.1;
        });

        array.row_mut(2).iter_mut().for_each(|value| {
            *value += 0.25;
        });

        let mvnpdf = &MultiNormalDensity::from_view::<Dyn, U3>(
            &array.as_view(),
            MDomain::new(SVector::from([
                SDomain::new(Some(-0.75), Some(0.75)).unwrap(),
                SDomain::new(Some(-0.75), Some(0.75)).unwrap(),
                SDomain::new(Some(-0.75), Some(0.75)).unwrap(),
            ])),
            None,
        )
        .unwrap();

        assert!(ulps_eq!(
            mvnpdf
                .density::<U1, U3>(&SVector::from([0.2, 0.0, 0.35]).as_view())
                .unwrap(),
            0.16128485,
            epsilon = 1e-5,
            max_ulps = 5
        ));

        assert!(
            mvnpdf
                .density::<U1, U3>(&SVector::from([0.2f64, 1.0, 0.35]).as_view())
                .is_none()
        );

        assert!(ulps_eq!(
            mvnpdf.sample(&mut rng, 100).unwrap(),
            SVector::from([-0.418523, 0.0, 0.4995714]),
            epsilon = 1e-5,
            max_ulps = 5
        ));

        assert!(
            mvnpdf
                .domain()
                .contains::<U1, U3>(&mvnpdf.sample(&mut rng, 100).unwrap().as_view())
        );

        let mvpdf_ensbl = MultiNormalDensity::from_view::<Dyn, U3>(
            &array.as_view(),
            MDomain::new(SVector::from([
                SDomain::new(Some(-0.75), Some(0.75)).unwrap(),
                SDomain::new(Some(-0.75), Some(0.75)).unwrap(),
                SDomain::new(Some(-0.75), Some(0.75)).unwrap(),
            ])),
            None,
        )
        .unwrap();

        assert!(ulps_eq!(
            mvnpdf.ltm,
            mvpdf_ensbl.ltm,
            epsilon = 1e-5,
            max_ulps = 5
        ));
    }

    #[test]
    fn test_multinormal_density_kld() {
        let mut rng = Xoshiro256PlusPlus::seed_from_u64(1);
        let uniform = StandardNormal;

        let array_1 = Matrix::<f64, U3, Dyn, VecStorage<f64, U3, Dyn>>::from_iterator(
            10000,
            (0..30000).map(|idx| {
                if idx % 3 == 1 {
                    0.0
                } else {
                    rng.sample::<f64, StandardNormal>(uniform)
                }
            }),
        );

        let array_2 = Matrix::<f64, U3, Dyn, VecStorage<f64, U3, Dyn>>::from_iterator(
            10000,
            (0..30000).map(|idx| {
                if idx % 3 == 1 {
                    0.0
                } else {
                    0.25 + rng.sample::<f64, StandardNormal>(uniform)
                }
            }),
        );

        let mvpdf_1 = MultiNormalDensity::from_view::<Dyn, U3>(
            &array_1.as_view(),
            MDomain::new(SVector::from([
                SDomain::new(Some(-0.75), Some(0.75)).unwrap(),
                SDomain::new(Some(-0.75), Some(0.75)).unwrap(),
                SDomain::new(Some(-0.75), Some(0.75)).unwrap(),
            ])),
            None,
        )
        .unwrap();
        let mvpdf_2 = MultiNormalDensity::from_view::<Dyn, U3>(
            &array_2.as_view(),
            MDomain::new(SVector::from([
                SDomain::new(Some(-0.75), Some(0.75)).unwrap(),
                SDomain::new(Some(-0.75), Some(0.75)).unwrap(),
                SDomain::new(Some(-0.75), Some(0.75)).unwrap(),
            ])),
            None,
        )
        .unwrap();

        assert!(ulps_eq!(
            mvpdf_1.kl_div(&mvpdf_2).unwrap(),
            0.181914,
            epsilon = 1e-5,
            max_ulps = 5
        ));

        assert!(ulps_eq!(
            mvpdf_2.kl_div(&mvpdf_1).unwrap(),
            0.166584,
            epsilon = 1e-5,
            max_ulps = 5
        ));
    }
}