r2rs-stats 0.1.1

use crate::sexprec::{SEXP, SEXPREC, SEXPTYPE};
use ::libc;
extern "C" {
    #[no_mangle]
    fn dcgettext(
        __domainname: *const libc::c_char,
        __msgid: *const libc::c_char,
        __category: libc::c_int,
    ) -> *mut libc::c_char;
    #[no_mangle]
    fn sqrt(_: libc::c_double) -> libc::c_double;

    /* General Support Functions */
    #[no_mangle]
    fn imax2(_: libc::c_int, _: libc::c_int) -> libc::c_int;
    #[no_mangle]
    fn imin2(_: libc::c_int, _: libc::c_int) -> libc::c_int;
    #[no_mangle]
    fn fmax2(_: libc::c_double, _: libc::c_double) -> libc::c_double;
    #[no_mangle]
    fn error(_: *const libc::c_char, _: ...) -> !;
    #[no_mangle]
    fn rPsort(_: *mut libc::c_double, _: libc::c_int, _: libc::c_int);
    #[no_mangle]
    fn TYPEOF(x: SEXP) -> libc::c_int;
    /* Vector Access Functions */
    #[no_mangle]
    fn LENGTH(x: SEXP) -> libc::c_int;
    #[no_mangle]
    fn REAL(x: SEXP) -> *mut libc::c_double;
    #[no_mangle]
    fn asInteger(x: SEXP) -> libc::c_int;
    #[no_mangle]
    fn asReal(x: SEXP) -> libc::c_double;
    /* NA_REAL: IEEE */
    #[no_mangle]
    static mut R_NaInt: libc::c_int;
    #[no_mangle]
    fn R_alloc(_: size_t, _: libc::c_int) -> *mut libc::c_char;
    #[no_mangle]
    fn allocVector(_: SEXPTYPE, _: R_xlen_t) -> SEXP;
    #[no_mangle]
    fn protect(_: SEXP) -> SEXP;
    #[no_mangle]
    fn unprotect(_: libc::c_int);
}
pub type Rboolean = libc::c_uint;
pub const TRUE: Rboolean = 1;
pub const FALSE: Rboolean = 0;
pub type size_t = libc::c_ulong;
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 */
/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 1996  Robert Gentleman and Ross Ihaka
 *  Copyright (C) 1999-2016 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/
 */
#[inline]
unsafe extern "C" fn fsquare(mut x: libc::c_double) -> libc::c_double {
    return x * x;
}
#[inline]
unsafe extern "C" fn fcube(mut x: libc::c_double) -> libc::c_double {
    return x * x * x;
}
unsafe extern "C" fn lowest(
    mut x: *mut libc::c_double,
    mut y: *mut libc::c_double,
    mut n: libc::c_int,
    mut xs: *mut libc::c_double,
    mut ys: *mut libc::c_double,
    mut nleft: libc::c_int,
    mut nright: libc::c_int,
    mut w: *mut libc::c_double,
    mut userw: Rboolean,
    mut rw: *mut libc::c_double,
    mut ok: *mut Rboolean,
) {
    let mut nrt: libc::c_int = 0;
    let mut j: libc::c_int = 0;
    let mut a: libc::c_double = 0.;
    let mut b: libc::c_double = 0.;
    let mut c: libc::c_double = 0.;
    let mut h: libc::c_double = 0.;
    let mut h1: libc::c_double = 0.;
    let mut h9: libc::c_double = 0.;
    let mut r: libc::c_double = 0.;
    let mut range: libc::c_double = 0.;
    x = x.offset(-1);
    y = y.offset(-1);
    w = w.offset(-1);
    rw = rw.offset(-1);
    range = *x.offset(n as isize) - *x.offset(1 as libc::c_int as isize);
    h = fmax2(
        *xs - *x.offset(nleft as isize),
        *x.offset(nright as isize) - *xs,
    );
    h9 = 0.999f64 * h;
    h1 = 0.001f64 * h;
    /* sum of weights */
    a = 0.0f64;
    j = nleft;
    while j <= n {
        /* compute weights */
        /* (pick up all ties on right) */
        *w.offset(j as isize) = 0.0f64;
        r = (*x.offset(j as isize) - *xs).abs();
        if r <= h9 {
            if r <= h1 {
                *w.offset(j as isize) = 1.0f64
            } else {
                *w.offset(j as isize) = fcube(1.0f64 - fcube(r / h))
            }
            if userw as u64 != 0 {
                *w.offset(j as isize) *= *rw.offset(j as isize)
            }
            a += *w.offset(j as isize)
        } else if *x.offset(j as isize) > *xs {
            break;
        }
        j = j + 1 as libc::c_int
    }
    /* rightmost pt (may be greater */
    /* than nright because of ties) */
    nrt = j - 1 as libc::c_int;
    if a <= 0.0f64 {
        *ok = FALSE
    } else {
        *ok = TRUE;
        /* weighted least squares */
        /* make sum of w[j] == 1 */
        j = nleft;
        while j <= nrt {
            *w.offset(j as isize) /= a;
            j += 1
        }
        if h > 0.0f64 {
            a = 0.0f64;
            /*  use linear fit */
            /* weighted center of x values */
            j = nleft;
            while j <= nrt {
                a += *w.offset(j as isize) * *x.offset(j as isize);
                j += 1
            }
            b = *xs - a;
            c = 0.0f64;
            j = nleft;
            while j <= nrt {
                c += *w.offset(j as isize) * fsquare(*x.offset(j as isize) - a);
                j += 1
            }
            if sqrt(c) > 0.001f64 * range {
                b /= c;
                /* points are spread out */
                /* enough to compute slope */
                j = nleft;
                while j <= nrt {
                    *w.offset(j as isize) *= b * (*x.offset(j as isize) - a) + 1.0f64;
                    j += 1
                }
            }
        }
        *ys = 0.0f64;
        j = nleft;
        while j <= nrt {
            *ys += *w.offset(j as isize) * *y.offset(j as isize);
            j += 1
        }
    };
}
unsafe extern "C" fn clowess(
    mut x: *mut libc::c_double,
    mut y: *mut libc::c_double,
    mut n: libc::c_int,
    mut f: libc::c_double,
    mut nsteps: libc::c_int,
    mut delta: libc::c_double,
    mut ys: *mut libc::c_double,
    mut rw: *mut libc::c_double,
    mut res: *mut libc::c_double,
) {
    let mut i: libc::c_int = 0;
    let mut iter: libc::c_int = 0;
    let mut j: libc::c_int = 0;
    let mut last: libc::c_int = 0;
    let mut m1: libc::c_int = 0;
    let mut m2: libc::c_int = 0;
    let mut nleft: libc::c_int = 0;
    let mut nright: libc::c_int = 0;
    let mut ns: libc::c_int = 0;
    let mut ok: Rboolean = FALSE;
    let mut alpha: libc::c_double = 0.;
    let mut c1: libc::c_double = 0.;
    let mut c9: libc::c_double = 0.;
    let mut cmad: libc::c_double = 0.;
    let mut cut: libc::c_double = 0.;
    let mut d1: libc::c_double = 0.;
    let mut d2: libc::c_double = 0.;
    let mut denom: libc::c_double = 0.;
    let mut r: libc::c_double = 0.;
    let mut sc: libc::c_double = 0.;
    if n < 2 as libc::c_int {
        *ys.offset(0 as libc::c_int as isize) = *y.offset(0 as libc::c_int as isize);
        return;
    }
    /* nleft, nright, last, etc. must all be shifted to get rid of these: */
    x = x.offset(-1);
    y = y.offset(-1);
    ys = ys.offset(-1);
    /* at least two, at most n points */
    ns = imax2(
        2 as libc::c_int,
        imin2(n, (f * n as libc::c_double + 1e-7f64) as libc::c_int),
    );
    /* robustness iterations */
    iter = 1 as libc::c_int; /* index of prev estimated point */
    while iter <= nsteps + 1 as libc::c_int {
        nleft = 1 as libc::c_int; /* index of current point */
        nright = ns;
        last = 0 as libc::c_int;
        i = 1 as libc::c_int;
        loop {
            if nright < n {
                /* move nleft,  nright to right */
                /* if radius decreases */
                d1 = *x.offset(i as isize) - *x.offset(nleft as isize);
                d2 = *x.offset((nright + 1 as libc::c_int) as isize) - *x.offset(i as isize);
                /* if d1 <= d2 with */
                /* x[nright+1] == x[nright], */
                /* lowest fixes */
                if d1 > d2 {
                    /* radius will not */
                    /* decrease by */
                    /* move right */
                    nleft += 1;
                    nright += 1;
                    continue;
                }
            }
            /* fitted value at x[i] */
            lowest(
                &mut *x.offset(1 as libc::c_int as isize),
                &mut *y.offset(1 as libc::c_int as isize),
                n,
                &mut *x.offset(i as isize),
                &mut *ys.offset(i as isize),
                nleft,
                nright,
                res,
                (iter > 1 as libc::c_int) as libc::c_int as Rboolean,
                rw,
                &mut ok,
            );
            if ok as u64 == 0 {
                *ys.offset(i as isize) = *y.offset(i as isize)
            }
            /* all weights zero */
            /* copy over value (all rw==0) */
            if last < i - 1 as libc::c_int {
                denom = *x.offset(i as isize) - *x.offset(last as isize);
                /* skipped points -- interpolate */
                /* non-zero - proof? */
                j = last + 1 as libc::c_int;
                while j < i {
                    alpha = (*x.offset(j as isize) - *x.offset(last as isize)) / denom;
                    *ys.offset(j as isize) = alpha * *ys.offset(i as isize)
                        + (1.0f64 - alpha) * *ys.offset(last as isize);
                    j += 1
                }
            }
            /* last point actually estimated */
            last = i;
            /* x coord of close points */
            cut = *x.offset(last as isize) + delta;
            i = last + 1 as libc::c_int;
            while i <= n {
                if *x.offset(i as isize) > cut {
                    break;
                }
                if *x.offset(i as isize) == *x.offset(last as isize) {
                    *ys.offset(i as isize) = *ys.offset(last as isize);
                    last = i
                }
                i += 1
            }
            i = imax2(last + 1 as libc::c_int, i - 1 as libc::c_int);
            if last >= n {
                break;
            }
        }
        /* residuals */
        i = 0 as libc::c_int;
        while i < n {
            *res.offset(i as isize) = *y.offset((i + 1 as libc::c_int) as isize)
                - *ys.offset((i + 1 as libc::c_int) as isize);
            i += 1
        }
        /* overall scale estimate */
        sc = 0.0f64;
        i = 0 as libc::c_int;
        while i < n {
            sc += (*res.offset(i as isize)).abs();
            i += 1
        }
        sc /= n as libc::c_double;
        /* compute robustness weights */
        /* except last time */
        if iter > nsteps {
            break;
        }
        /* Note: The following code, biweight_{6 MAD|Ri|}
        is also used in stl(), loess and several other places.
        --> should provide API here (MM) */
        i = 0 as libc::c_int;
        while i < n {
            *rw.offset(i as isize) = (*res.offset(i as isize)).abs();
            i += 1
        }
        /* Compute   cmad := 6 * median(rw[], n)  ---- */
        /* FIXME: We need C API in R for Median ! */
        m1 = n / 2 as libc::c_int;
        /* partial sort, for m1 & m2 */
        rPsort(rw, n, m1);
        if n % 2 as libc::c_int == 0 as libc::c_int {
            m2 = n - m1 - 1 as libc::c_int;
            rPsort(rw, n, m2);
            cmad = 3.0f64 * (*rw.offset(m1 as isize) + *rw.offset(m2 as isize))
        } else {
            /* n odd */
            cmad = 6.0f64 * *rw.offset(m1 as isize)
        }
        if cmad < 1e-7f64 * sc {
            break;
        }
        c9 = 0.999f64 * cmad;
        c1 = 0.001f64 * cmad;
        i = 0 as libc::c_int;
        while i < n {
            r = (*res.offset(i as isize)).abs();
            if r <= c1 {
                *rw.offset(i as isize) = 1.0f64
            } else if r <= c9 {
                *rw.offset(i as isize) = fsquare(1.0f64 - fsquare(r / cmad))
            } else {
                *rw.offset(i as isize) = 0.0f64
            }
            i += 1
        }
        iter += 1
    }
}
#[no_mangle]
pub unsafe extern "C" fn lowess(
    mut x: SEXP,
    mut y: SEXP,
    mut sf: SEXP,
    mut siter: SEXP,
    mut sdelta: SEXP,
) -> SEXP {
    if TYPEOF(x) != 14 as libc::c_int || TYPEOF(y) != 14 as libc::c_int {
        error(b"invalid input\x00" as *const u8 as *const libc::c_char);
    }
    let mut nx: libc::c_int = LENGTH(x);
    if nx == R_NaInt || nx == 0 as libc::c_int {
        error(b"invalid input\x00" as *const u8 as *const libc::c_char);
    }
    let mut f: libc::c_double = asReal(sf);
    if f.is_finite() as i32 == 0 || f <= 0 as libc::c_int as libc::c_double {
        error(dcgettext(
            b"stats\x00" as *const u8 as *const libc::c_char,
            b"\'f\' must be finite and > 0\x00" as *const u8 as *const libc::c_char,
            5 as libc::c_int,
        ));
    }
    let mut iter: libc::c_int = asInteger(siter);
    if iter == R_NaInt || iter < 0 as libc::c_int {
        error(dcgettext(
            b"stats\x00" as *const u8 as *const libc::c_char,
            b"\'iter\' must be finite and >= 0\x00" as *const u8 as *const libc::c_char,
            5 as libc::c_int,
        ));
    }
    let mut delta: libc::c_double = asReal(sdelta);
    let mut rw: *mut libc::c_double = 0 as *mut libc::c_double;
    let mut res: *mut libc::c_double = 0 as *mut libc::c_double;
    if delta.is_finite() as i32 == 0 || delta < 0 as libc::c_int as libc::c_double {
        error(dcgettext(
            b"stats\x00" as *const u8 as *const libc::c_char,
            b"\'delta\' must be finite and > 0\x00" as *const u8 as *const libc::c_char,
            5 as libc::c_int,
        ));
    }
    let mut ans: SEXP = 0 as *mut SEXPREC;
    ans = allocVector(14 as libc::c_int as SEXPTYPE, nx as R_xlen_t);
    protect(ans);
    rw = R_alloc(
        nx as size_t,
        ::std::mem::size_of::<libc::c_double>() as libc::c_ulong as libc::c_int,
    ) as *mut libc::c_double;
    res = R_alloc(
        nx as size_t,
        ::std::mem::size_of::<libc::c_double>() as libc::c_ulong as libc::c_int,
    ) as *mut libc::c_double;
    clowess(REAL(x), REAL(y), nx, f, iter, delta, REAL(ans), rw, res);
    unprotect(1 as libc::c_int);
    return ans;
}