r2rs-stats 0.1.1

Statistics programming for Rust based on R's stats package
Documentation
use crate::sexprec::{SEXP, SEXPTYPE};
use ::libc;
extern "C" {
    #[no_mangle]
    fn sqrt(_: libc::c_double) -> libc::c_double;
    /* Vector Access Functions */
    #[no_mangle]
    fn LENGTH(x: SEXP) -> libc::c_int;
    #[no_mangle]
    fn INTEGER(x: SEXP) -> *mut libc::c_int;
    #[no_mangle]
    fn REAL(x: SEXP) -> *mut libc::c_double;
    #[no_mangle]
    fn SET_VECTOR_ELT(x: SEXP, i: R_xlen_t, v: SEXP) -> SEXP;
    #[no_mangle]
    fn asInteger(x: SEXP) -> libc::c_int;
    #[no_mangle]
    fn asReal(x: SEXP) -> libc::c_double;
    #[no_mangle]
    fn error(_: *const libc::c_char, _: ...) -> !;
    /*
     *  R : A Computer Language for Statistical Data Analysis
     *  Copyright (C) 1995, 1996  Robert Gentleman and Ross Ihaka
     *  Copyright (C) 1998--2016  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 */
    /*
       This used to define _BSD_SOURCE to make declarations of isfinite
       and isnan visible in glibc.  But that was deprecated in glibc 2.20,
       and --std=c99 suffices nowadays.
    */
    /* needed for isnan and isfinite, neither of which are used under C++ */
    /* implementation of these : ../../main/arithmetic.c */
    /* IEEE NaN */
    /* IEEE Inf */
    /* IEEE -Inf */
    /* NA_REAL: IEEE */
    /* NA_INTEGER:= INT_MIN currently */
    /* #define NA_FACTOR R_NaInt  unused */
    /* NA_STRING is a SEXP, so defined in Rinternals.h */
    /* True for R's NA only */
    /* True for special NaN, *not* for NA */
    #[no_mangle]
    fn R_finite(_: libc::c_double) -> libc::c_int;
    #[no_mangle]
    static mut R_NegInf: libc::c_double;
    #[no_mangle]
    static mut R_PosInf: libc::c_double;
    /* Defining NO_RINLINEDFUNS disables use to simulate platforms where
    this is not available */
    /* need remapped names here for use with R_NO_REMAP */
    /*
       These are the inlinable functions that are provided in Rinlinedfuns.h
       It is *essential* that these do not appear in any other header file,
       with or without the Rf_ prefix.
    */
    #[no_mangle]
    fn allocVector(_: SEXPTYPE, _: R_xlen_t) -> SEXP;
    #[no_mangle]
    fn ScalarReal(_: libc::c_double) -> SEXP;
    #[no_mangle]
    fn protect(_: SEXP) -> SEXP;
    #[no_mangle]
    fn unprotect(_: libc::c_int);
    #[no_mangle]
    fn dcgettext(
        __domainname: *const libc::c_char,
        __msgid: *const libc::c_char,
        __category: libc::c_int,
    ) -> *mut libc::c_char;
}
pub type ptrdiff_t = libc::c_long;
/* both config.h and Rconfig.h set SIZEOF_SIZE_T, but Rconfig.h is
skipped if config.h has already been included. */
pub type R_xlen_t = ptrdiff_t;
/* Fundamental Data Types:  These are largely Lisp
 * influenced structures, with the exception of LGLSXP,
 * INTSXP, REALSXP, CPLXSXP and STRSXP which are the
 * element types for S-like data objects.
 *
 *   --> TypeTable[] in ../main/util.c for  typeof()
 */
/* UUID identifying the internals version -- packages using compiled
code should be re-installed when this changes */
/*  These exact numeric values are seldom used, but they are, e.g., in
 *  ../main/subassign.c, and they are serialized.
*/
/* NOT YET using enum:
 *  1) The SEXPREC struct below has 'SEXPTYPE type : 5'
 * (making FUNSXP and CLOSXP equivalent in there),
 * giving (-Wall only ?) warnings all over the place
 * 2) Many switch(type) { case ... } statements need a final `default:'
 * added in order to avoid warnings like [e.g. l.170 of ../main/util.c]
 *   "enumeration value `FUNSXP' not handled in switch"
 */
/* real variables */
/* complex variables */
/* string vectors */
/* dot-dot-dot object */
/* make "any" args work.
Used in specifying types for symbol
registration to mean anything is okay  */
/* generic vectors */
/* expressions vectors */
/* byte code */
/* external pointer */
/* weak reference */
/* raw bytes */
/* S4, non-vector */
/* used for detecting PROTECT issues in memory.c */
/* fresh node created in new page */
/* node released by GC */
/* Closure or Builtin or Special */
/* NOT YET */
/* These are also used with the write barrier on, in attrib.c and util.c */
/* Avoid slow and possibly error-producing underflows by cutting off at
plus/minus sqrt(DELMAX) std deviations */
/* Formulae (6.67) and (6.69) of Scott (1992), the latter corrected. */
#[no_mangle]
pub unsafe extern "C" fn bw_ucv(mut sn: SEXP, mut sd: SEXP, mut cnt: SEXP, mut sh: SEXP) -> SEXP {
    let mut h: libc::c_double = asReal(sh);
    let mut d: libc::c_double = asReal(sd);
    let mut sum: libc::c_double = 0.0f64;
    let mut term: libc::c_double = 0.;
    let mut u: libc::c_double = 0.;
    let mut n: libc::c_int = asInteger(sn);
    let mut nbin: libc::c_int = LENGTH(cnt);
    let mut x: *mut libc::c_double = REAL(cnt);
    let mut i: libc::c_int = 0 as libc::c_int;
    while i < nbin {
        let mut delta: libc::c_double = i as libc::c_double * d / h;
        delta *= delta;
        if delta >= 1000 as libc::c_int as libc::c_double {
            break;
        }
        term = (-delta / 4.0f64).exp() - sqrt(8.0f64) * (-delta / 2.0f64).exp();
        sum += term * *x.offset(i as isize);
        i += 1
    }
    u = (0.5f64 + sum / n as libc::c_double)
        / (n as libc::c_double * h * 1.772453850905516027298167483341f64);
    // = 1 / (2 * n * h * sqrt(PI)) + sum / (n * n * h * sqrt(PI));
    return ScalarReal(u);
}
#[no_mangle]
pub unsafe extern "C" fn bw_bcv(mut sn: SEXP, mut sd: SEXP, mut cnt: SEXP, mut sh: SEXP) -> SEXP {
    let mut h: libc::c_double = asReal(sh);
    let mut d: libc::c_double = asReal(sd);
    let mut sum: libc::c_double = 0.0f64;
    let mut term: libc::c_double = 0.;
    let mut u: libc::c_double = 0.;
    let mut n: libc::c_int = asInteger(sn);
    let mut nbin: libc::c_int = LENGTH(cnt);
    let mut x: *mut libc::c_double = REAL(cnt);
    sum = 0.0f64;
    let mut i: libc::c_int = 0 as libc::c_int;
    while i < nbin {
        let mut delta: libc::c_double = i as libc::c_double * d / h;
        delta *= delta;
        if delta >= 1000 as libc::c_int as libc::c_double {
            break;
        }
        term = (-delta / 4 as libc::c_int as libc::c_double).exp()
            * (delta * delta - 12 as libc::c_int as libc::c_double * delta
                + 12 as libc::c_int as libc::c_double);
        sum += term * *x.offset(i as isize);
        i += 1
    }
    u = (1 as libc::c_int as libc::c_double + sum / (32.0f64 * n as libc::c_double))
        / (2.0f64 * n as libc::c_double * h * 1.772453850905516027298167483341f64);
    // = 1 / (2 * n * h * sqrt(PI)) + sum / (64 * n * n * h * sqrt(PI));
    return ScalarReal(u); /* add in diagonal */
}
#[no_mangle]
pub unsafe extern "C" fn bw_phi4(mut sn: SEXP, mut sd: SEXP, mut cnt: SEXP, mut sh: SEXP) -> SEXP {
    let mut h: libc::c_double = asReal(sh);
    let mut d: libc::c_double = asReal(sd);
    let mut sum: libc::c_double = 0.0f64;
    let mut term: libc::c_double = 0.;
    let mut u: libc::c_double = 0.;
    let mut n: libc::c_int = asInteger(sn);
    let mut nbin: libc::c_int = LENGTH(cnt);
    let mut x: *mut libc::c_double = REAL(cnt);
    let mut i: libc::c_int = 0 as libc::c_int;
    while i < nbin {
        let mut delta: libc::c_double = i as libc::c_double * d / h;
        delta *= delta;
        if delta >= 1000 as libc::c_int as libc::c_double {
            break;
        }
        term = (-delta / 2.0f64).exp() * (delta * delta - 6.0f64 * delta + 3.0f64);
        sum += term * *x.offset(i as isize);
        i += 1
    }
    sum = 2.0f64 * sum + n as libc::c_double * 3.0f64;
    u = sum / (n as libc::c_double * (n - 1 as libc::c_int) as libc::c_double * h.powf(5.0f64))
        * 0.398942280401432677939946059934f64;
    // = sum / (n * (n - 1) * h.powf(5.0) * sqrt(2 * PI));
    return ScalarReal(u); /* add in diagonal */
}
#[no_mangle]
pub unsafe extern "C" fn bw_phi6(mut sn: SEXP, mut sd: SEXP, mut cnt: SEXP, mut sh: SEXP) -> SEXP {
    let mut h: libc::c_double = asReal(sh);
    let mut d: libc::c_double = asReal(sd);
    let mut sum: libc::c_double = 0.0f64;
    let mut term: libc::c_double = 0.;
    let mut u: libc::c_double = 0.;
    let mut n: libc::c_int = asInteger(sn);
    let mut nbin: libc::c_int = LENGTH(cnt);
    let mut x: *mut libc::c_double = REAL(cnt);
    let mut i: libc::c_int = 0 as libc::c_int;
    while i < nbin {
        let mut delta: libc::c_double = i as libc::c_double * d / h;
        delta *= delta;
        if delta >= 1000 as libc::c_int as libc::c_double {
            break;
        }
        term = (-delta / 2 as libc::c_int as libc::c_double).exp()
            * (delta * delta * delta - 15 as libc::c_int as libc::c_double * delta * delta
                + 45 as libc::c_int as libc::c_double * delta
                - 15 as libc::c_int as libc::c_double);
        sum += term * *x.offset(i as isize);
        i += 1
    }
    sum = 2.0f64 * sum - 15.0f64 * n as libc::c_double;
    u = sum / (n as libc::c_double * (n - 1 as libc::c_int) as libc::c_double * h.powf(7.0f64))
        * 0.398942280401432677939946059934f64;
    // = sum / (n * (n - 1) * h.powf(7.0) * sqrt(2 * PI));
    return ScalarReal(u);
}
/*
   Use double cnt as from R 3.4.0, as counts can exceed INT_MAX for
   large n (65537 in the worse case but typically not at n = 1 million
   for a smooth distribution -- and this is by default no longer used
   for n > 500).
*/
#[no_mangle]
pub unsafe extern "C" fn bw_den(mut nbin: SEXP, mut sx: SEXP) -> SEXP {
    let mut nb: libc::c_int = asInteger(nbin);
    let mut n: libc::c_int = LENGTH(sx);
    let mut xmin: libc::c_double = 0.;
    let mut xmax: libc::c_double = 0.;
    let mut rang: libc::c_double = 0.;
    let mut dd: libc::c_double = 0.;
    let mut x: *mut libc::c_double = REAL(sx);
    xmin = R_PosInf;
    xmax = R_NegInf;
    let mut i: libc::c_int = 0 as libc::c_int;
    while i < n {
        if R_finite(*x.offset(i as isize)) == 0 {
            error(
                dcgettext(
                    b"stats\x00" as *const u8 as *const libc::c_char,
                    b"non-finite x[%d] in bandwidth calculation\x00" as *const u8
                        as *const libc::c_char,
                    5 as libc::c_int,
                ),
                i + 1 as libc::c_int,
            );
        }
        if *x.offset(i as isize) < xmin {
            xmin = *x.offset(i as isize)
        }
        if *x.offset(i as isize) > xmax {
            xmax = *x.offset(i as isize)
        }
        i += 1
    }
    rang = (xmax - xmin) * 1.01f64;
    dd = rang / nb as libc::c_double;
    let mut ans: SEXP = protect(allocVector(
        19 as libc::c_int as SEXPTYPE,
        2 as libc::c_int as R_xlen_t,
    ));
    let mut sc: SEXP = SET_VECTOR_ELT(
        ans,
        1 as libc::c_int as R_xlen_t,
        allocVector(14 as libc::c_int as SEXPTYPE, nb as R_xlen_t),
    );
    SET_VECTOR_ELT(ans, 0 as libc::c_int as R_xlen_t, ScalarReal(dd));
    let mut cnt: *mut libc::c_double = REAL(sc);
    let mut i_0: libc::c_int = 0 as libc::c_int;
    while i_0 < nb {
        *cnt.offset(i_0 as isize) = 0.0f64;
        i_0 += 1
    }
    let mut i_1: libc::c_int = 1 as libc::c_int;
    while i_1 < n {
        let mut ii: libc::c_int = (*x.offset(i_1 as isize) / dd) as libc::c_int;
        let mut j: libc::c_int = 0 as libc::c_int;
        while j < i_1 {
            let mut jj: libc::c_int = (*x.offset(j as isize) / dd) as libc::c_int;
            *cnt.offset((ii - jj).abs() as isize) += 1.0f64;
            j += 1
        }
        i_1 += 1
    }
    unprotect(1 as libc::c_int);
    return ans;
}
/* Input: counts for nb bins */
#[no_mangle]
pub unsafe extern "C" fn bw_den_binned(mut sx: SEXP) -> SEXP {
    let mut nb: libc::c_int = LENGTH(sx); // don't count distances to self
    let mut x: *mut libc::c_int = INTEGER(sx); // counts in the same bin got double-counted
    let mut ans: SEXP = protect(allocVector(14 as libc::c_int as SEXPTYPE, nb as R_xlen_t));
    let mut cnt: *mut libc::c_double = REAL(ans);
    let mut ib: libc::c_int = 0 as libc::c_int;
    while ib < nb {
        *cnt.offset(ib as isize) = 0.0f64;
        ib += 1
    }
    let mut ii: libc::c_int = 0 as libc::c_int;
    while ii < nb {
        let mut w: libc::c_int = *x.offset(ii as isize);
        *cnt.offset(0 as libc::c_int as isize) +=
            w as libc::c_double * (w as libc::c_double - 1.0f64);
        let mut jj: libc::c_int = 0 as libc::c_int;
        while jj < ii {
            *cnt.offset((ii - jj) as isize) += (w * *x.offset(jj as isize)) as libc::c_double;
            jj += 1
        }
        ii += 1
    }
    *cnt.offset(0 as libc::c_int as isize) *= 0.5f64;
    unprotect(1 as libc::c_int);
    return ans;
}