1 /*
2 * Copyright (C) 2023 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical.interpolation;
17
18 /**
19 * Computes barycentric rational interpolation.
20 */
21 public class BarycentricRationalInterpolator extends BaseInterpolator {
22
23 /**
24 * Weights for barycentric rational interpolation.
25 */
26 private final double[] w;
27
28 /**
29 * Order of desired approximation.
30 */
31 private final int d;
32
33 /**
34 * Constructor.
35 *
36 * @param x x values to interpolate to. Values in x must be monotonic (either increasing or
37 * decreasing)
38 * @param y y values to interpolate to.
39 * @param d order of desired approximation.
40 */
41 public BarycentricRationalInterpolator(final double[] x, final double[] y, final int d) {
42 super(x, y, x.length);
43 w = new double[n];
44 this.d = d;
45
46 if (n <= d) {
47 throw new IllegalArgumentException("d too large for number of points");
48 }
49
50 for (int k = 0; k < n; k++) {
51 var imin = Math.max(k - d, 0);
52 var imax = k >= n - d ? n - d - 1 : k;
53 var temp = (imin & 1) != 0 ? -1.0 : 1.0;
54 var sum = 0.0;
55 for (int i = imin; i <= imax; i++) {
56 var jmax = Math.min(i + d, n - 1);
57 var term = 1.0;
58 for (var j = i; j <= jmax; j++) {
59 if (j == k) {
60 continue;
61 }
62 term *= (xx[k] - xx[j]);
63 }
64 term = temp / term;
65 temp = -temp;
66 sum += term;
67 }
68 w[k] = sum;
69 }
70 }
71
72 /**
73 * Gets order of desired approximation.
74 *
75 * @return order of desired approximation.
76 */
77 public int getD() {
78 return d;
79 }
80
81 /**
82 * Given a value x, returns an interpolated value, using data pointed to by {@link #xx} and
83 * {@link #yy}.
84 *
85 * @param x value to obtain interpolation for.
86 * @return interpolated value.
87 */
88 @Override
89 public double interpolate(final double x) {
90 return rawinterp(1, x);
91 }
92
93 /**
94 * Actual interpolation method.
95 *
96 * @param jlo index where value x to be interpolated in located in the array of xx.
97 * @param x value to obtain interpolation for.
98 * @return interpolated value.
99 */
100 @Override
101 public double rawinterp(int jlo, double x) {
102 var num = 0.0;
103 var den = 0.0;
104 for (var i = 0; i < n; i++) {
105 var h = x - xx[i];
106 if (h == 0.0) {
107 return yy[i];
108 } else {
109 var temp = w[i] / h;
110 num += temp * yy[i];
111 den += temp;
112 }
113 }
114 return num / den;
115 }
116 }