1 /* --- Isotonic regression ---
2 * code simplified from VR_mds_fn() which is part of MASS.c,
3 * Copyright (C) 1995 Brian Ripley
4 * ---
5 * R : A Computer Language for Statistical Data Analysis
6 * Copyright (C) 2003 The R Foundation
7 *
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, a copy is available at
20 * https://www.R-project.org/Licenses/
21 */
22
23 #include "modreg.h"
24
isoreg(SEXP y)25 SEXP isoreg(SEXP y)
26 {
27 int n = LENGTH(y), i, ip, known, n_ip;
28 double tmp, slope;
29 SEXP yc, yf, iKnots, ans;
30 const char *anms[] = {"y", "yc", "yf", "iKnots", ""};
31
32 /* unneeded: y = coerceVector(y, REALSXP); */
33
34 PROTECT(ans = mkNamed(VECSXP, anms));
35
36 SET_VECTOR_ELT(ans, 0, y);
37 SET_VECTOR_ELT(ans, 1, yc = allocVector(REALSXP, n+1));
38 SET_VECTOR_ELT(ans, 2, yf = allocVector(REALSXP, n));
39 SET_VECTOR_ELT(ans, 3, iKnots= allocVector(INTSXP, n));
40
41 if (n == 0) {
42 UNPROTECT(1); /* ans */
43 return ans; /* avoid segfault below */
44 }
45
46 /* yc := cumsum(0,y) */
47 REAL(yc)[0] = 0.;
48 tmp = 0.;
49 for (i = 0; i < n; i++) {
50 tmp += REAL(y)[i];
51 REAL(yc)[i + 1] = tmp;
52 }
53 known = 0; ip = 0, n_ip = 0;
54 do {
55 slope = R_PosInf;/*1e+200*/
56 for (i = known + 1; i <= n; i++) {
57 tmp = (REAL(yc)[i] - REAL(yc)[known]) / (i - known);
58 if (tmp < slope) {
59 slope = tmp;
60 ip = i;
61 }
62 }/* tmp := max{i= kn+1,.., n} slope(p[kn] -> p[i]) and
63 * ip = argmax{...}... */
64 INTEGER(iKnots)[n_ip++] = ip;
65 for (i = known; i < ip; i++)
66 REAL(yf)[i] = (REAL(yc)[ip] - REAL(yc)[known]) / (ip - known);
67 } while ((known = ip) < n);
68
69 if (n_ip < n)
70 SET_VECTOR_ELT(ans, 3, lengthgets(iKnots, n_ip));
71 UNPROTECT(1);
72 return(ans);
73 }
74