1 /* R : A Computer Language for Statistical Data Analysis
2 * Copyright (C) 2010 The R Foundation
3 * Copyright (C) 2016 The R Core Team
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, a copy is available at
17 * https://www.R-project.org/Licenses/
18 */
19
20 #include "modreg.h"
21 #include <math.h>
22
23 /* To be "exported" (as part of R's C API): */
24 /**
25 * Modify the slopes m_k := s'(x_k) using Fritsch & Carlson (1980)'s algorithm
26 *
27 * @param m numeric vector of length n, the preliminary desired slopes s'(x_i), i = 1:n
28 * @param S the divided differences (y_{i+1} - y_i) / (x_{i+1} - x_i); i = 1:(n-1)
29 * @param n == length(m) == 1 + length(S)
30 * @return m*: the modified m[]'s: Note that m[] is modified in place
31 * @author Martin Maechler, Date: 19 Apr 2010
32 */
monoFC_mod(double * m,double S[],int n)33 void monoFC_mod(double *m, double S[], int n)
34 {
35 if(n < 2)
36 error(_("n must be at least two"));
37
38 for(int k = 0; k < n - 1; k++) {
39 /* modify both (m[k] & m[k+1]) if needed : */
40 double Sk = S[k];
41 int k1 = k + 1;
42 if(Sk == 0.) { /* or |S| < eps ?? FIXME ?? */
43 m[k] = m[k1] = 0.;
44 } else {
45 double
46 alpha = m[k ] / Sk,
47 beta = m[k1] / Sk, a2b3, ab23;
48 if((a2b3 = 2*alpha + beta - 3) > 0 &&
49 (ab23 = alpha + 2*beta - 3) > 0 &&
50 alpha * (a2b3 + ab23) < a2b3*a2b3) {
51 /* we are outside the monotonocity region ==> fix slopes */
52 double tauS = 3*Sk / sqrt(alpha*alpha + beta*beta);
53 m[k ] = tauS * alpha;
54 m[k1] = tauS * beta;
55 }
56 }
57 } /* end for */
58 }
59
monoFC_m(SEXP m,SEXP Sx)60 SEXP monoFC_m(SEXP m, SEXP Sx)
61 {
62 SEXP val;
63 int n = LENGTH(m);
64
65 if (isInteger(m))
66 val = PROTECT(coerceVector(m, REALSXP));
67 else {
68 if (!isReal(m))
69 error(_("Argument m must be numeric"));
70 val = PROTECT(duplicate(m));
71 }
72 if(n < 2) error(_("length(m) must be at least two"));
73 if(!isReal(Sx) || LENGTH(Sx) != n-1)
74 error(_("Argument Sx must be numeric vector one shorter than m[]"));
75
76 /* Fix up the slopes m[] := val[]: */
77 monoFC_mod(REAL(val), REAL(Sx), n);
78
79 UNPROTECT(1);
80 return(val);
81 }
82