r2rs-stats 0.1.1

// "Whatever you do, work at it with all your heart, as working for the Lord,
// not for human masters, since you know that you will receive an inheritance
// from the Lord as a reward. It is the Lord Christ you are serving."
// (Col 3:23-24)

use r2rs_base::func::partial_sort;

pub fn lowess(x: &[f64], y: &[f64]) -> (Vec<f64>, Vec<f64>) {
    let mut x_ordered = vec![0.0; x.len()];
    let mut y_ordered = vec![0.0; x.len()];
    let mut indices = (0..x.len()).collect::<Vec<_>>();
    indices.sort_by(|&i1, &i2| x[i1].partial_cmp(&x[i2]).unwrap());
    for (old_index, new_index) in indices.iter().cloned().enumerate() {
        x_ordered[old_index] = x[new_index];
        y_ordered[old_index] = y[new_index];
    }
    let x = x_ordered;
    let y = y_ordered;

    let mut ys = vec![0.0; y.len()];
    let mut res = vec![0.0; y.len()];
    let mut rw = vec![0.0; y.len()];

    let n = x.len();
    let f = 2.0 / 3.0;
    let iter = 3;
    let delta = 0.01 * (x[n - 1] - x[0]);
    let nsteps = iter;

    let ns = 2.max(n.min((f * n as f64 + 1e-7) as usize));

    // robustness iterations

    let mut d1 = 0.0;
    let mut d2 = 0.0;

    let mut iter = 1;
    while iter < nsteps + 1 {
        let mut nleft = 1;
        let mut nright = ns;
        let mut last = 1; // index of prev estimated point
        let mut i = 1; // index of current point

        loop {
            if nright < n {
                // move nleft,  nright to right if radius decreases
                d1 = x[i - 1] - x[nleft - 1];
                d2 = x[nright + 1 - 1] - x[i - 1];

                // 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]

            let ok = lowest(
                &x,
                &y,
                n,
                &x,
                &mut ys,
                i,
                nleft,
                nright,
                &mut res,
                iter > 0,
                &rw,
            );
            if !ok {
                ys[i - 1] = y[i - 1];
            }

            // all weights zero
            // copy over value (all rw==0)

            if last < i {
                // Maybe last < i - 1
                let denom = x[i - 1] - x[last - 1];

                // skipped points -- interpolate
                // non-zero - proof?

                for j in last + 1..i {
                    let alpha = (x[j - 1] - x[last - 1]) / denom;
                    ys[j - 1] = alpha * ys[i - 1] + (1.0 - alpha) * ys[last - 1];
                }
            }

            // last point actually estimated
            last = i;

            // x coord of close points
            let cut = x[last - 1] + delta;
            for i2 in last + 1..=n {
                if x[i2 - 1] > cut {
                    i = i2;
                    break;
                }
                if x[i2 - 1] == x[last - 1] {
                    ys[i2 - 1] = ys[last - 1];
                    last = i2;
                }
            }
            i = (last + 1).max(i - 1);
            if last >= n {
                break;
            }
        }
        // residuals
        for i in 0..n {
            res[i] = y[i + 1 - 1] - ys[i + 1 - 1];
        }

        // overall scale estimate
        let mut sc = 0.0;
        for i in 0..n {
            sc += res[i].abs();
        }
        sc /= n as f64;

        // 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)
        for i in 0..n {
            rw[i] = res[i].abs();
        }

        // Compute   cmad := 6 * median(rw[], n)  ----
        // FIXME: We need C API in R for Median !
        let m1 = n / 2;
        // partial sort, for m1 & m2
        partial_sort(&mut rw, &[n - 1, m1 - 1]).expect("Error in partial sort");
        // 	rPsort(rw, n, m1);
        let cmad = if n % 2 == 0 {
            let m2 = n - m1 - 1;
            partial_sort(&mut rw, &[n - 1, m2 - 1]).expect("Error in partial sort");
            // rPsort(rw, n, m2);
            3.0 * (rw[m1] + rw[m2])
        } else {
            // n odd
            6.0 * rw[m1]
        };
        if cmad < 1e-7 * sc {
            // effectively zero
            break;
        }
        let c9 = 0.999 * cmad;
        let c1 = 0.001 * cmad;
        for i in 0..n {
            let r = res[i].abs();
            if r <= c1 {
                rw[i] = 1.0;
            } else if r <= c9 {
                rw[i] = fsquare(1.0 - fsquare(r / cmad));
            } else {
                rw[i] = 0.0;
            }
        }
        iter += 1;
    }

    (x, ys)
}

fn fcube(x: f64) -> f64 {
    x * x * x
}

fn fsquare(x: f64) -> f64 {
    x * x
}

fn lowest(
    x: &[f64],
    y: &[f64],
    n: usize,
    xs: &[f64],
    ys: &mut [f64],
    i: usize,
    nleft: usize,
    nright: usize,
    w: &mut [f64],
    userw: bool,
    rw: &[f64],
) -> bool {
    let range = x[n - 1] - x[0];
    let h = (xs[i - 1] - x[nleft - 1]).max(x[nright - 1] - xs[i - 1]);
    let h9 = 0.999 * h;
    let h1 = 0.001 * h;

    // sum of weights

    let mut a = 0.;
    let mut j = nleft;
    while j <= n {
        // compute weights (pick up all ties on right)

        w[j - 1] = 0.0;
        let r = (x[j - 1] - xs[i - 1]).abs();
        if r <= h9 {
            if r <= h1 {
                w[j - 1] = 1.0;
            } else {
                w[j - 1] = fcube(1.0 - fcube(r / h));
            }
            if userw {
                w[j - 1] *= rw[j - 1];
            }
            a += w[j - 1];
        } else if x[j - 1] > xs[i - 1] {
            break;
        }
        j = j + 1;
    }

    // rightmost pt (may be greater than nright because of ties)

    let nrt = j - 1;
    if a <= 0.0 {
        false
    } else {
        // weighted least squares
        // make sum of w[j] == 1

        for j in nleft..=nrt {
            w[j - 1] /= a;
        }
        if h > 0.0 {
            a = 0.0;

            // use linear fit weighted center of x values

            for j in nleft..=nrt {
                a += w[j - 1] * x[j - 1];
            }
            let mut b = xs[i - 1] - a;
            let mut c = 0.0;
            for j in nleft..=nrt {
                c += w[j - 1] * fsquare(x[j - 1] - a);
            }
            if c.sqrt() > 0.001 * range {
                b /= c;

                // points are spread out enough to compute slope

                for j in nleft..=nrt {
                    w[j - 1] *= b * (x[j - 1] - a) + 1.0;
                }
            }
        }
        ys[i - 1] = 0.0;
        for j in nleft..=nrt {
            ys[i - 1] += w[j - 1] * y[j - 1];
        }

        true
    }
}