r2rs-stats 0.1.1

use std::{
    error::Error,
    fmt::{Display, Error as FmtError, Formatter},
};

use ::libc;
use r2rs_nmath::{distribution::NormalBuilder, traits::Distribution};
use strafe_trait::{
    Assumption, Concept, Conclusion, Normal, NotNormal, RejectionStatus, Statistic, StatisticalTest,
};
use strafe_type::{Alpha64, FloatConstraint};

#[derive(Copy, Clone, Debug)]
pub enum ShapiroWilkError {
    TooFew,
    TooMany,
    ZeroRange,
    Error7,
}

impl Display for ShapiroWilkError {
    fn fmt(&self, fmt: &mut Formatter) -> Result<(), FmtError> {
        match &self {
            ShapiroWilkError::TooFew => write!(fmt, "Too few observations"),
            ShapiroWilkError::TooMany => write!(fmt, "Too many observations"),
            ShapiroWilkError::ZeroRange => write!(fmt, "Observation range is zero"),
            ShapiroWilkError::Error7 => write!(fmt, "Error 7"),
        }
    }
}

impl Error for ShapiroWilkError {}

#[derive(Copy, Clone, Debug)]
pub struct ShapiroWilkTest {
    alpha: Alpha64,
}

impl StatisticalTest for ShapiroWilkTest {
    type Input = Vec<f64>;
    type Output = Result<ShapiroWilkStatistic, ShapiroWilkError>;

    fn assumptions() -> Vec<Box<dyn Assumption>> {
        Vec::new()
    }

    fn null_hypotheses() -> Vec<Box<dyn Conclusion>> {
        vec![Box::new(Normal::new())]
    }

    fn alternate_hypotheses() -> Vec<Box<dyn Conclusion>> {
        vec![Box::new(NotNormal::new())]
    }

    fn test(&mut self, x: &Self::Input) -> Self::Output {
        let mut x_sorted = x.to_owned();
        x_sorted.sort_by(|a1, a2| a1.partial_cmp(a2).unwrap());

        let mut w = 0.0;
        let mut pw = 0.0;
        let mut ifault = 0;
        swilk(&x_sorted, x.len() as i32, &mut w, &mut pw, &mut ifault)
            .expect("Error running swilk");

        match ifault {
            1 => Err(ShapiroWilkError::TooFew),
            2 => Err(ShapiroWilkError::TooMany),
            6 => Err(ShapiroWilkError::ZeroRange),
            7 => Err(ShapiroWilkError::Error7),
            _ => Ok(ShapiroWilkStatistic {
                alpha: self.alpha,
                w,
                p: pw,
            }),
        }
    }
}

#[derive(Copy, Clone, Debug)]
pub struct ShapiroWilkStatistic {
    pub(crate) alpha: Alpha64,
    pub(crate) w: f64,
    pub(crate) p: f64,
}

impl Statistic for ShapiroWilkStatistic {
    fn alpha(&self) -> Alpha64 {
        self.alpha
    }

    fn statistic(&self) -> f64 {
        self.w
    }

    fn probability_value(&self) -> f64 {
        self.p
    }

    fn conclusion(&self) -> RejectionStatus {
        if self.p < self.alpha.unwrap() {
            RejectionStatus::RejectInFavorOf(Box::new(NotNormal {}))
        } else {
            RejectionStatus::FailToReject(Box::new(Normal {}))
        }
    }
}

impl Default for ShapiroWilkTest {
    fn default() -> Self {
        Self { alpha: 0.05.into() }
    }
}

impl ShapiroWilkTest {
    pub fn new() -> Self {
        Self::default()
    }

    pub fn with_alpha<A: Into<Alpha64>>(self, alpha: A) -> Self {
        Self {
            alpha: alpha.into(),
        }
    }
}

pub fn swilk(
    mut x: &[libc::c_double],
    mut n: libc::c_int,
    mut w: &mut libc::c_double,
    mut pw: &mut libc::c_double,
    mut ifault: &mut libc::c_int,
) -> Result<(), Box<dyn Error>> {
    let mut nn2: libc::c_int = n / 2 as libc::c_int; /* 1-based */
    let vla = (nn2 + 1 as libc::c_int) as usize;
    let mut a: Vec<libc::c_double> = ::std::vec::from_elem(0., vla);
    /* ALGORITHM AS R94 APPL. STATIST. (1995) vol.44, no.4, 547-551.

        Calculates the Shapiro-Wilk W test and its significance level
    */
    let mut small: libc::c_double = 1e-19f64;
    /* polynomial coefficients */
    let mut g: [libc::c_double; 2] = [-2.273f64, 0.459f64];
    let mut c1: [libc::c_double; 6] = [
        0.0f64,
        0.221157f64,
        -0.147981f64,
        -2.07119f64,
        4.434685f64,
        -2.706056f64,
    ];
    let mut c2: [libc::c_double; 6] = [
        0.0f64,
        0.042981f64,
        -0.293762f64,
        -1.752461f64,
        5.682633f64,
        -3.582633f64,
    ];
    let mut c3: [libc::c_double; 4] = [0.544f64, -0.39978f64, 0.025054f64, -6.714e-4f64];
    let mut c4: [libc::c_double; 4] = [1.3822f64, -0.77857f64, 0.062767f64, -0.0020322f64];
    let mut c5: [libc::c_double; 4] = [-1.5861f64, -0.31082f64, -0.083751f64, 0.0038915f64];
    let mut c6: [libc::c_double; 3] = [-0.4803f64, -0.082676f64, 0.0030302f64];
    /* Local variables */
    let mut i: libc::c_int = 0;
    let mut j: libc::c_int = 0;
    let mut i1: libc::c_int = 0;
    let mut ssassx: libc::c_double = 0.;
    let mut summ2: libc::c_double = 0.;
    let mut ssumm2: libc::c_double = 0.;
    let mut gamma: libc::c_double = 0.;
    let mut range: libc::c_double = 0.;
    let mut a1: libc::c_double = 0.;
    let mut a2: libc::c_double = 0.;
    let mut an: libc::c_double = 0.;
    let mut m: libc::c_double = 0.;
    let mut s: libc::c_double = 0.;
    let mut sa: libc::c_double = 0.;
    let mut xi: libc::c_double = 0.;
    let mut sx: libc::c_double = 0.;
    let mut xx: libc::c_double = 0.;
    let mut y: libc::c_double = 0.;
    let mut w1: libc::c_double = 0.;
    let mut fac: libc::c_double = 0.;
    let mut asa: libc::c_double = 0.;
    let mut an25: libc::c_double = 0.;
    let mut ssa: libc::c_double = 0.;
    let mut sax: libc::c_double = 0.;
    let mut rsn: libc::c_double = 0.;
    let mut ssx: libc::c_double = 0.;
    let mut xsx: libc::c_double = 0.;
    *pw = 1.0f64;
    if n < 3 as libc::c_int {
        *ifault = 1 as libc::c_int;
        return Ok(());
    }
    an = n as libc::c_double;
    if n == 3 as libc::c_int {
        a[1] = std::f64::consts::FRAC_1_SQRT_2
    /* = sqrt(1/2) */
    } else {
        an25 = an + 0.25f64;
        summ2 = 0.0f64;
        i = 1 as libc::c_int;
        while i <= nn2 {
            let norm = NormalBuilder::new()
                .with_mean(0.0)
                .with_standard_deviation(1.0)
                .build();
            a[i as usize] = norm
                .quantile((i as libc::c_double - 0.375f64) / an25, true)
                .unwrap();
            let mut r__1: libc::c_double = a[i as usize];
            summ2 += r__1 * r__1;
            i += 1
        }
        summ2 *= 2.0f64;
        ssumm2 = summ2.sqrt();
        rsn = 1.0f64 / an.sqrt();
        a1 = poly(&c1, 6 as libc::c_int, rsn) - a[1] / ssumm2;
        /* Normalize a[] */
        if n > 5 as libc::c_int {
            i1 = 3 as libc::c_int;
            a2 = -a[2] / ssumm2 + poly(&c2, 6 as libc::c_int, rsn);
            fac = ((summ2 - 2.0f64 * (a[1] * a[1]) - 2.0f64 * (a[2] * a[2]))
                / (1.0f64 - 2.0f64 * (a1 * a1) - 2.0f64 * (a2 * a2)))
                .sqrt();
            a[2] = a2
        } else {
            i1 = 2 as libc::c_int;
            fac = ((summ2 - 2.0f64 * (a[1] * a[1])) / (1.0f64 - 2.0f64 * (a1 * a1))).sqrt()
        }
        a[1] = a1;
        i = i1;
        while i <= nn2 {
            a[i as usize] /= -fac;
            i += 1
        }
    }
    /* Check for zero range */
    range = x[(n - 1) as usize] - x[0];
    if range < small {
        *ifault = 6 as libc::c_int;
        return Ok(());
    }
    /* Check for correct sort order on range - scaled X */
    /* *ifault = 7; <-- a no-op, since it is changed below, in ANY CASE! */
    *ifault = 0 as libc::c_int;
    xx = x[0] / range;
    sx = xx;
    sa = -a[1];
    i = 1 as libc::c_int;
    j = n - 1 as libc::c_int;
    while i < n {
        xi = x[i as usize] / range;
        if xx - xi > small {
            /* Fortran had:  print *, "ANYTHING"
             * but do NOT; it *does* happen with sorted x (on Intel GNU/linux 32bit):
             *  shapiro.test(c(-1.7, -1,-1,-.73,-.61,-.5,-.24, .45,.62,.81,1))
             */
            *ifault = 7 as libc::c_int
        }
        sx += xi;
        i += 1;
        if i != j {
            sa += ((i - j) as libc::c_double).signum() * a[(if i > j { j } else { i }) as usize]
        }
        xx = xi;
        j -= 1
    }
    if n > 5000 as libc::c_int {
        *ifault = 2 as libc::c_int
    }
    /* Calculate W statistic as squared correlation
    between data and coefficients */
    sa /= n as libc::c_double;
    sx /= n as libc::c_double;
    sax = 0.0f64;
    ssx = sax;
    ssa = ssx;
    i = 0 as libc::c_int;
    j = n - 1 as libc::c_int;
    while i < n {
        if i != j {
            asa = ((i - j) as libc::c_double).signum()
                * a[(1 + (if i > j { j } else { i })) as usize]
                - sa
        } else {
            asa = -sa
        }
        xsx = x[i as usize] / range - sx;
        ssa += asa * asa;
        ssx += xsx * xsx;
        sax += asa * xsx;
        i += 1;
        j -= 1
    }
    /* W1 equals (1-W) calculated to avoid excessive rounding error
    for W very near 1 (a potential problem in very large samples) */
    ssassx = (ssa * ssx).sqrt();
    w1 = (ssassx - sax) * (ssassx + sax) / (ssa * ssx);
    *w = 1.0f64 - w1;
    /* Calculate significance level for W */
    if n == 3 as libc::c_int {
        /* exact P value : */
        let mut pi6: libc::c_double = 1.90985931710274f64; /* = asin(sqrt(3/4)) */
        let mut stqr: libc::c_double = std::f64::consts::FRAC_PI_3; /* n >= 12 */
        *pw = pi6 * (((*w).sqrt()).asin() - stqr); /* an "obvious" value, was 'small' which was 1e-19f */
        if *pw < 0.0f64 {
            *pw = 0.0f64
        }
        return Ok(());
    }
    y = w1.ln();
    xx = an.ln();
    if n <= 11 as libc::c_int {
        gamma = poly(&g, 2 as libc::c_int, an);
        if y >= gamma {
            *pw = 1e-99f64;
            return Ok(());
        }
        y = -(gamma - y).ln();
        m = poly(&c3, 4 as libc::c_int, an);
        s = (poly(&c4, 4 as libc::c_int, an)).exp()
    } else {
        m = poly(&c5, 4 as libc::c_int, xx);
        s = (poly(&c6, 3 as libc::c_int, xx)).exp()
    }
    /*DBG printf("c(w1=%g, w=%g, y=%g, m=%g, s=%g)\n",w1,*w,y,m,s); */
    let norm = NormalBuilder::new()
        .with_mean(m)
        .with_standard_deviation(s)
        .build();
    *pw = norm.probability(y, false).unwrap();

    Ok(())
}
/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 2000-2016   The R Core Team.
 *
 *  Based on Applied Statistics algorithms AS181, R94
 *    (C) Royal Statistical Society 1982, 1995
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, a copy is available at
 *  https://www.R-project.org/Licenses/
 */
/* swilk.f -- translated by f2c (version 19980913).
 * ------- and produced by f2c-clean,v 1.8 --- and hand polished: M.Maechler
 */
/* swilk */
fn poly(mut cc: &[libc::c_double], mut nord: libc::c_int, mut x: libc::c_double) -> libc::c_double {
    /* Algorithm AS 181.2 Appl. Statist. (1982) Vol. 31, No. 2

        Calculates the algebraic polynomial of order nord-1 with
        array of coefficients cc.  Zero order coefficient is cc(1) = cc[0]
    */
    let mut p: libc::c_double = 0.;
    let mut ret_val: libc::c_double = 0.;
    ret_val = cc[0];
    if nord > 1 as libc::c_int {
        p = x * cc[(nord - 1 as libc::c_int) as usize];
        let mut j: libc::c_int = nord - 2 as libc::c_int;
        while j > 0 as libc::c_int {
            p = (p + cc[j as usize]) * x;
            j -= 1
        }
        ret_val += p
    }
    ret_val
}
/* poly */