r2rs-stats 0.1.1

Statistics programming for Rust based on R's stats package
Documentation
use crate::sexprec::{SEXP, SEXPREC, SEXPTYPE};
use ::libc;
extern "C" {
    #[no_mangle]
    fn R_rsort(_: *mut libc::c_double, _: libc::c_int);

    /*
     *  R : A Computer Language for Statistical Data Analysis
     *  Copyright (C) 1998-2005   The R Core Team
     *
     *  This header file is free software; you can redistribute it and/or modify
     *  it under the terms of the GNU Lesser General Public License as published by
     *  the Free Software Foundation; either version 2.1 of the License, or
     *  (at your option) any later version.
     *
     *  This file is part of R. R is distributed under the terms of the
     *  GNU General Public License, either Version 2, June 1991 or Version 3,
     *  June 2007. See doc/COPYRIGHTS for details of the copyright status of R.
     *
     *  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 Lesser General Public License for more details.
     *
     *  You should have received a copy of the GNU Lesser General Public License
     *  along with this program; if not, a copy is available at
     *  https://www.R-project.org/Licenses/
     */
    /* Included by R.h: API */
    #[no_mangle]
    fn error(_: *const libc::c_char, _: ...) -> !;
    #[no_mangle]
    fn LENGTH(x: SEXP) -> libc::c_int;
    #[no_mangle]
    fn REAL(x: SEXP) -> *mut libc::c_double;
    #[no_mangle]
    fn SET_STRING_ELT(x: SEXP, i: R_xlen_t, v: SEXP);
    #[no_mangle]
    fn SET_VECTOR_ELT(x: SEXP, i: R_xlen_t, v: SEXP) -> SEXP;
    #[no_mangle]
    static mut R_NamesSymbol: SEXP;
    #[no_mangle]
    fn asInteger(x: SEXP) -> libc::c_int;
    #[no_mangle]
    fn mkChar(_: *const libc::c_char) -> SEXP;
    #[no_mangle]
    fn setAttrib(_: SEXP, _: SEXP, _: SEXP) -> SEXP;
    #[no_mangle]
    fn allocVector(_: SEXPTYPE, _: R_xlen_t) -> SEXP;
    #[no_mangle]
    fn protect(_: SEXP) -> SEXP;
    #[no_mangle]
    fn unprotect(_: libc::c_int);
}
pub type ptrdiff_t = libc::c_long;
pub type R_xlen_t = ptrdiff_t;
/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 1997-2017   The R Core Team.
 *
 *  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/
 */
/* Speed up by `inlining' these (as macros) [since R version 1.2] : */
unsafe extern "C" fn line(
    mut x: *mut libc::c_double,
    mut y: *mut libc::c_double,
    mut z: *mut libc::c_double,
    mut w: *mut libc::c_double,
    mut n: libc::c_int,
    mut iter: libc::c_int,
    mut coef: *mut libc::c_double,
) {
    let mut i: libc::c_int = 0; /* z = ordered abscissae */
    let mut j: libc::c_int = 0;
    let mut k: libc::c_int = 0;
    let mut xb: libc::c_double = 0.;
    let mut x1: libc::c_double = 0.;
    let mut x2: libc::c_double = 0.;
    let mut xt: libc::c_double = 0.;
    let mut yt: libc::c_double = 0.;
    let mut yb: libc::c_double = 0.;
    let mut tmp1: libc::c_double = 0.;
    let mut tmp2: libc::c_double = 0.;
    i = 0 as libc::c_int;
    while i < n {
        *z.offset(i as isize) = *x.offset(i as isize);
        *w.offset(i as isize) = *y.offset(i as isize);
        i += 1
    }
    R_rsort(z, n);
    // All x-space computations can (and hence should) happen outside the iterations
    /* NB: Using quantiles is *not* quite correct (not following the reference).
           Rather need 3 groups with group sizes depending on n and (n % 3)
    */
    // x1 := quantile(x, 1/3) :
    tmp1 = *z.offset(
        ((n - 1 as libc::c_int) as libc::c_double * 1.0f64 / 3.0f64).floor() as libc::c_int
            as isize,
    );
    tmp2 = *z.offset(
        ((n - 1 as libc::c_int) as libc::c_double * 1.0f64 / 3.0f64).ceil() as libc::c_int as isize,
    );
    x1 = 0.5f64 * (tmp1 + tmp2);
    // x1 := quantile(x, 2/3) :
    tmp1 = *z.offset(
        ((n - 1 as libc::c_int) as libc::c_double * 2.0f64 / 3.0f64).floor() as libc::c_int
            as isize,
    );
    tmp2 = *z.offset(
        ((n - 1 as libc::c_int) as libc::c_double * 2.0f64 / 3.0f64).ceil() as libc::c_int as isize,
    );
    x2 = 0.5f64 * (tmp1 + tmp2);
    // xb := x_L := Median(x[i]; x[i] <= quantile(x, 1/3))
    i = 0 as libc::c_int;
    k = 0 as libc::c_int;
    while i < n {
        if *x.offset(i as isize) <= x1 {
            let fresh0 = k;
            k = k + 1;
            *z.offset(fresh0 as isize) = *x.offset(i as isize)
        }
        i += 1
    }
    R_rsort(z, k);
    xb = 0.5f64
        * (*z.offset(
            ((k - 1 as libc::c_int) as libc::c_double * 0.5f64).floor() as libc::c_int as isize,
        ) + *z.offset(
            ((k - 1 as libc::c_int) as libc::c_double * 0.5f64).ceil() as libc::c_int as isize,
        ));
    // xt := x_R := Median(x[i]; x[i] >= quantile(x, 2/3))
    i = 0 as libc::c_int;
    k = 0 as libc::c_int;
    while i < n {
        if *x.offset(i as isize) >= x2 {
            let fresh1 = k;
            k = k + 1;
            *z.offset(fresh1 as isize) = *x.offset(i as isize)
        }
        i += 1
    }
    R_rsort(z, k);
    xt = 0.5f64
        * (*z.offset(
            ((k - 1 as libc::c_int) as libc::c_double * 0.5f64).floor() as libc::c_int as isize,
        ) + *z.offset(
            ((k - 1 as libc::c_int) as libc::c_double * 0.5f64).ceil() as libc::c_int as isize,
        ));
    let mut slope: libc::c_double = 0.0f64;
    // "Polishing" iterations (incl first estimate) :
    j = 1 as libc::c_int; // iterations
    while j <= iter {
        // w[] = "current y"
        /* yb := Median(y[i]; x[i] <= quantile(x, 1/3) */
        i = 0 as libc::c_int;
        k = 0 as libc::c_int;
        while i < n {
            if *x.offset(i as isize) <= x1 {
                let fresh2 = k;
                k = k + 1;
                *z.offset(fresh2 as isize) = *w.offset(i as isize)
            }
            i += 1
        }
        R_rsort(z, k);
        yb = 0.5f64
            * (*z.offset(
                ((k - 1 as libc::c_int) as libc::c_double * 0.5f64).floor() as libc::c_int as isize,
            ) + *z.offset(
                ((k - 1 as libc::c_int) as libc::c_double * 0.5f64).ceil() as libc::c_int as isize
            ));
        /* yt := Median(y[i]; x[i] >= quantile(x, 2/3) */
        i = 0 as libc::c_int;
        k = 0 as libc::c_int;
        while i < n {
            if *x.offset(i as isize) >= x2 {
                let fresh3 = k;
                k = k + 1;
                *z.offset(fresh3 as isize) = *w.offset(i as isize)
            }
            i += 1
        }
        R_rsort(z, k);
        yt = 0.5f64
            * (*z.offset(
                ((k - 1 as libc::c_int) as libc::c_double * 0.5f64).floor() as libc::c_int as isize,
            ) + *z.offset(
                ((k - 1 as libc::c_int) as libc::c_double * 0.5f64).ceil() as libc::c_int as isize
            ));
        slope += (yt - yb) / (xt - xb);
        i = 0 as libc::c_int;
        while i < n {
            *w.offset(i as isize) = *y.offset(i as isize) - slope * *x.offset(i as isize);
            i += 1
        }
        j += 1
    }
    // intercept := median{ "residuals" } (= mean of the 3 median residuals)
    R_rsort(w, n);
    let mut yint: libc::c_double = 0.5f64
        * (*w.offset(
            ((n - 1 as libc::c_int) as libc::c_double * 0.5f64).floor() as libc::c_int as isize,
        ) + *w.offset(
            ((n - 1 as libc::c_int) as libc::c_double * 0.5f64).ceil() as libc::c_int as isize,
        ));
    i = 0 as libc::c_int;
    while i < n {
        *w.offset(i as isize) = yint + slope * *x.offset(i as isize);
        *z.offset(i as isize) = *y.offset(i as isize) - *w.offset(i as isize);
        i += 1
    }
    *coef.offset(0 as libc::c_int as isize) = yint;
    *coef.offset(1 as libc::c_int as isize) = slope;
}
// where is this used?
#[no_mangle]
pub unsafe extern "C" fn tukeyline0(
    mut x: *mut libc::c_double,
    mut y: *mut libc::c_double,
    mut z: *mut libc::c_double,
    mut w: *mut libc::c_double,
    mut n: *mut libc::c_int,
    mut coef: *mut libc::c_double,
) {
    line(x, y, z, w, *n, 1 as libc::c_int, coef);
}
#[no_mangle]
pub unsafe extern "C" fn tukeyline(
    mut x: SEXP,
    mut y: SEXP,
    mut iter: SEXP,
    mut call: SEXP,
) -> SEXP {
    let mut n: libc::c_int = LENGTH(x);
    if n < 2 as libc::c_int {
        error(b"insufficient observations\x00" as *const u8 as *const libc::c_char);
    }
    let mut ans: SEXP = 0 as *mut SEXPREC;
    ans = protect(allocVector(
        19 as libc::c_int as SEXPTYPE,
        4 as libc::c_int as R_xlen_t,
    ));
    let mut nm: SEXP = allocVector(16 as libc::c_int as SEXPTYPE, 4 as libc::c_int as R_xlen_t);
    setAttrib(ans, R_NamesSymbol, nm);
    SET_STRING_ELT(
        nm,
        0 as libc::c_int as R_xlen_t,
        mkChar(b"call\x00" as *const u8 as *const libc::c_char),
    );
    SET_STRING_ELT(
        nm,
        1 as libc::c_int as R_xlen_t,
        mkChar(b"coefficients\x00" as *const u8 as *const libc::c_char),
    );
    SET_STRING_ELT(
        nm,
        2 as libc::c_int as R_xlen_t,
        mkChar(b"residuals\x00" as *const u8 as *const libc::c_char),
    );
    SET_STRING_ELT(
        nm,
        3 as libc::c_int as R_xlen_t,
        mkChar(b"fitted.values\x00" as *const u8 as *const libc::c_char),
    );
    SET_VECTOR_ELT(ans, 0 as libc::c_int as R_xlen_t, call);
    let mut coef: SEXP = allocVector(14 as libc::c_int as SEXPTYPE, 2 as libc::c_int as R_xlen_t);
    SET_VECTOR_ELT(ans, 1 as libc::c_int as R_xlen_t, coef);
    let mut res: SEXP = allocVector(14 as libc::c_int as SEXPTYPE, n as R_xlen_t);
    SET_VECTOR_ELT(ans, 2 as libc::c_int as R_xlen_t, res);
    let mut fit: SEXP = allocVector(14 as libc::c_int as SEXPTYPE, n as R_xlen_t);
    SET_VECTOR_ELT(ans, 3 as libc::c_int as R_xlen_t, fit);
    line(
        REAL(x),
        REAL(y),
        REAL(res),
        REAL(fit),
        n,
        asInteger(iter),
        REAL(coef),
    );
    unprotect(1 as libc::c_int);
    return ans;
}