mod d;
mod p;
mod q;
mod r;

use strafe_type::{LogProbability64, Positive64, Probability64, Real64};

pub(crate) use self::{d::*, p::*, q::*, r::*};
use crate::traits::{Distribution, RNG};

/// # The Normal Distribution
///
/// ## Description
///
/// Density, distribution function, quantile function and random generation for the normal
/// distribution with mean equal to mean and standard deviation equal to sd.
///
/// ## Arguments
///
/// * mean: vector of means.
/// * sd: vector of standard deviations.
///
/// ## Details
///
/// If mean or sd are not specified they assume the default values of 0 and 1, respectively.
///
/// The normal distribution has density
///
/// $ f(x) = \frac{1}{\sigma \sqrt{2 \pi}} e^{-\frac{(x - \mu)^2}{2 \sigma^2}} $
///
/// where $ \mu $ is the mean of the distribution and $ \sigma $ the standard deviation.
///
/// ## Value
///
/// dnorm gives the density, pnorm gives the distribution function, qnorm gives the quantile
/// function, and rnorm generates random deviates.
///
/// The length of the result is determined by n for rnorm, and is the maximum of the lengths of
/// the numerical arguments for the other functions.
///
/// The numerical arguments other than n are recycled to the length of the result. Only the first
/// elements of the logical arguments are used.
///
/// For sd = 0 this gives the limit as sd decreases to 0, a point mass at mu. sd < 0 is an error
/// and returns NaN.
///
/// ## Density Plot
///
/// ```rust
/// # use r2rs_base::traits::StatisticalSlice;
/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
/// # use strafe_plot::prelude::{IntoDrawingArea, Line, Plot, PlotOptions, SVGBackend, BLACK};
/// # use strafe_type::FloatConstraint;
/// let norm = NormalBuilder::new().build();
/// let x = <[f64]>::sequence(-3.0, 3.0, 1000);
/// let y = x
///     .iter()
///     .map(|x| norm.density(x).unwrap())
///     .collect::<Vec<_>>();
///
/// let root = SVGBackend::new("density.svg", (1024, 768)).into_drawing_area();
/// Plot::new()
///     .with_options(PlotOptions {
///         x_axis_label: "x".to_string(),
///         y_axis_label: "density".to_string(),
///         ..Default::default()
///     })
///     .with_plottable(Line {
///         x,
///         y,
///         color: BLACK,
///         ..Default::default()
///     })
///     .plot(&root)
///     .unwrap();
/// # use std::fs::rename;
/// #     drop(root);
/// #     rename(
/// #             format!("density.svg"),
/// #             format!("src/distribution/norm/doctest_out/density.svg"),
/// #     )
/// #     .unwrap();
/// ```
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = embed_doc_image::embed_image!("density", "src/distribution/norm/doctest_out/density.svg")))]
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = "![Density][density]"))]
///
/// ## Source
///
/// For pnorm, based on
///
/// Cody, W. D. (1993) Algorithm 715: SPECFUN – A portable FORTRAN package of special function
/// routines and test drivers. ACM Transactions on Mathematical Software 19, 22–32.
///
/// For qnorm, the code is a C translation of
///
/// Wichura, M. J. (1988) Algorithm AS 241: The percentage points of the normal distribution.
/// Applied Statistics, 37, 477–484.
///
/// which provides precise results up to about 16 digits.
///
/// For rnorm, see RNG for how to select the algorithm and for references to the supplie
/// methods.
///
/// ## References
///
/// Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988) The New S Language. Wadsworth &
/// Brooks/Cole.
///
/// Johnson, N. L., Kotz, S. and Balakrishnan, N. (1995) Continuous Univariate Distributions,
/// volume 1, chapter 13. Wiley, New York.
///
/// ## See Also
///
/// Distributions for other standard distributions, including dlnorm for the Lognormal distribution.
///
/// ## Examples
///
/// ```rust
/// # use num_traits::FloatConst;
/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
/// let norm = NormalBuilder::new().build();
/// println!("{}", norm.density(0));
/// println!("{}", 1.0 / (2.0 * f64::PI()).sqrt());
/// # use std::{fs::File, io::Write};
/// # let mut f = File::create("src/distribution/norm/doctest_out/dens1.md").unwrap();
/// # writeln!(f, "```output").unwrap();
/// # writeln!(f, "{}", norm.density(0)).unwrap();
/// # writeln!(f, "{}", 1.0 / (2.0 * f64::PI()).sqrt()).unwrap();
/// # writeln!(f, "```").unwrap();
/// ```
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/dens1.md")))]
///
/// ```rust
/// # use num_traits::FloatConst;
/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
/// let norm = NormalBuilder::new().build();
/// println!("{}", norm.density(1));
/// println!("{}", (-0.5_f64).exp() / (2.0 * f64::PI()).sqrt());
/// println!("{}", 1.0 / (2.0 * f64::PI() * 1.0_f64.exp()).sqrt());
/// # use std::{fs::File, io::Write};
/// # let mut f = File::create("src/distribution/norm/doctest_out/dens2.md").unwrap();
/// # writeln!(f, "```output").unwrap();
/// # writeln!(f, "{}", norm.density(1)).unwrap();
/// # writeln!(f, "{}", (-0.5_f64).exp() / (2.0 * f64::PI()).sqrt()).unwrap();
/// # writeln!(f, "{}", 1.0 / (2.0 * f64::PI() * 1.0_f64.exp()).sqrt()).unwrap();
/// # writeln!(f, "```").unwrap();
/// ```
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = include_str!("doctest_out/dens2.md")))]
///
/// ```rust
/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
/// # use strafe_plot::{
/// #     plot::Plot,
/// #     prelude::{IntoDrawingArea, Line, PlotOptions, SVGBackend, BLACK, RED},
/// # };
/// # use strafe_type::FloatConstraint;
/// let norm = NormalBuilder::new().build();
///
/// let dx = (-60..50).map(|x| x as f64).collect::<Vec<_>>();
/// let dy1 = dx
///     .iter()
///     .map(|x| norm.log_density(x).unwrap())
///     .collect::<Vec<_>>();
/// let dy2 = dx
///     .iter()
///     .map(|x| norm.density(x).unwrap().ln())
///     .collect::<Vec<_>>();
///
/// let px = (-50..10).map(|x| x as f64).collect::<Vec<_>>();
/// let py1 = px
///     .iter()
///     .map(|x| norm.log_probability(x, true).unwrap())
///     .collect::<Vec<_>>();
/// let py2 = px
///     .iter()
///     .map(|x| norm.probability(x, true).unwrap().ln())
///     .collect::<Vec<_>>();
///
/// let root = SVGBackend::new("log_plots.svg", (1024, 768)).into_drawing_area();
///
/// Plot::new()
///     .with_options(PlotOptions {
///         x_axis_label: "x".to_string(),
///         y_axis_label: "density".to_string(),
///         title: "Log Normal Density".to_string(),
///         legend_x: 0.1,
///         legend_y: 0.1,
///         plot_bottom: 0.5,
///         ..Default::default()
///     })
///     .with_plottable(Line {
///         x: dx.clone(),
///         y: dy1,
///         legend: true,
///         label: "log_density()".to_string(),
///         color: BLACK,
///         ..Default::default()
///     })
///     .with_plottable(Line {
///         x: dx.clone(),
///         y: dy2,
///         legend: true,
///         label: "density().ln()".to_string(),
///         color: RED,
///         ..Default::default()
///     })
///     .plot(&root)
///     .unwrap();
///
/// Plot::new()
///     .with_options(PlotOptions {
///         x_axis_label: "x".to_string(),
///         y_axis_label: "probability".to_string(),
///         title: "Log Normal Cumulative".to_string(),
///         legend_x: 0.1,
///         legend_y: 0.1,
///         plot_top: 0.5,
///         ..Default::default()
///     })
///     .with_plottable(Line {
///         x: px.clone(),
///         y: py1,
///         legend: true,
///         label: "log_probability()".to_string(),
///         color: BLACK,
///         ..Default::default()
///     })
///     .with_plottable(Line {
///         x: px.clone(),
///         y: py2,
///         legend: true,
///         label: "probability().ln()".to_string(),
///         color: RED,
///         ..Default::default()
///     })
///     .plot(&root)
///     .unwrap();
/// # use std::fs::rename;
/// #     drop(root);
/// #     rename(
/// #             format!("log_plots.svg"),
/// #             format!("src/distribution/norm/doctest_out/log_plots.svg"),
/// #     )
/// #     .unwrap();
/// ```
#[cfg_attr(feature = "doc_outputs", cfg_attr(all(), doc = embed_doc_image::embed_image!("log_plots", "src/distribution/norm/doctest_out/log_plots.svg")))]
#[cfg_attr(
    feature = "doc_outputs",
    cfg_attr(all(), doc = "![Log Plots][log_plots]")
)]
///
/// If you want the so-called 'error function'
/// (see Abramowitz and Stegun 29.2.29)
/// ```rust
/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
/// # use strafe_type::FloatConstraint;
/// fn erf(x: f64) -> f64 {
///     let norm = NormalBuilder::new().build();
///     2.0 * norm.probability(x * 2.0_f64.sqrt(), true).unwrap() - 1.0
/// }
/// ```
///
/// and the so-called 'complementary error function'
/// ```rust
/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
/// # use strafe_type::FloatConstraint;
/// fn erfc(x: f64) -> f64 {
///     let norm = NormalBuilder::new().build();
///     2.0 * norm.probability(x * 2.0_f64.sqrt(), false).unwrap()
/// }
/// ```
///
/// and the inverses
/// ```rust
/// # use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
/// # use strafe_type::FloatConstraint;
/// fn erfinv(x: f64) -> f64 {
///     let norm = NormalBuilder::new().build();
///     norm.quantile((1.0 + x) / 2.0, true).unwrap() / 2.0_f64.sqrt()
/// }
///
/// fn erfcinv(x: f64) -> f64 {
///     let norm = NormalBuilder::new().build();
///     norm.quantile(x / 2.0, false).unwrap() / 2.0_f64.sqrt()
/// }
/// ```
pub struct Normal {
    mean: Real64,
    standard_deviation: Positive64,
}

impl Distribution for Normal {
    fn density<R: Into<Real64>>(&self, x: R) -> Real64 {
        dnorm(x, self.mean, self.standard_deviation, false)
    }

    fn log_density<R: Into<Real64>>(&self, x: R) -> Real64 {
        dnorm(x, self.mean, self.standard_deviation, true)
    }

    fn probability<R: Into<Real64>>(&self, q: R, lower_tail: bool) -> Probability64 {
        pnorm(q, self.mean, self.standard_deviation, lower_tail)
    }

    fn log_probability<R: Into<Real64>>(&self, q: R, lower_tail: bool) -> LogProbability64 {
        log_pnorm(q, self.mean, self.standard_deviation, lower_tail)
    }

    fn quantile<P: Into<Probability64>>(&self, p: P, lower_tail: bool) -> Real64 {
        qnorm(p, self.mean, self.standard_deviation, lower_tail)
    }

    fn log_quantile<LP: Into<LogProbability64>>(&self, p: LP, lower_tail: bool) -> Real64 {
        log_qnorm(p, self.mean, self.standard_deviation, lower_tail)
    }

    fn random_sample<R: RNG>(&self, rng: &mut R) -> Real64 {
        rnorm(self.mean, self.standard_deviation, rng)
    }
}

pub struct NormalBuilder {
    mean: Option<Real64>,
    standard_deviation: Option<Positive64>,
}

impl NormalBuilder {
    pub fn new() -> Self {
        Self {
            mean: None,
            standard_deviation: None,
        }
    }

    pub fn with_mean<R: Into<Real64>>(&mut self, mean: R) -> &mut Self {
        self.mean = Some(mean.into());
        self
    }

    pub fn with_standard_deviation<P: Into<Positive64>>(
        &mut self,
        standard_deviation: P,
    ) -> &mut Self {
        self.standard_deviation = Some(standard_deviation.into());
        self
    }

    pub fn build(&self) -> Normal {
        let mean = self.mean.unwrap_or(0.0.into());
        let standard_deviation = self.standard_deviation.unwrap_or(1.0.into());

        Normal {
            mean,
            standard_deviation,
        }
    }
}

#[cfg(test)]
mod tests;

#[cfg(all(test, feature = "enable_proptest"))]
mod proptests;

#[cfg(all(test, feature = "enable_covtest"))]
mod covtests;