1 // ******************************************************************************
2 //
3 // Title: Force Field X.
4 // Description: Force Field X - Software for Molecular Biophysics.
5 // Copyright: Copyright (c) Michael J. Schnieders 2001-2024.
6 //
7 // This file is part of Force Field X.
8 //
9 // Force Field X is free software; you can redistribute it and/or modify it
10 // under the terms of the GNU General Public License version 3 as published by
11 // the Free Software Foundation.
12 //
13 // Force Field X is distributed in the hope that it will be useful, but WITHOUT
14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
16 // details.
17 //
18 // You should have received a copy of the GNU General Public License along with
19 // Force Field X; if not, write to the Free Software Foundation, Inc., 59 Temple
20 // Place, Suite 330, Boston, MA 02111-1307 USA
21 //
22 // Linking this library statically or dynamically with other modules is making a
23 // combined work based on this library. Thus, the terms and conditions of the
24 // GNU General Public License cover the whole combination.
25 //
26 // As a special exception, the copyright holders of this library give you
27 // permission to link this library with independent modules to produce an
28 // executable, regardless of the license terms of these independent modules, and
29 // to copy and distribute the resulting executable under terms of your choice,
30 // provided that you also meet, for each linked independent module, the terms
31 // and conditions of the license of that module. An independent module is a
32 // module which is not derived from or based on this library. If you modify this
33 // library, you may extend this exception to your version of the library, but
34 // you are not obligated to do so. If you do not wish to do so, delete this
35 // exception statement from your version.
36 //
37 // ******************************************************************************
38 package ffx.numerics.spline;
39
40 import java.util.Arrays;
41
42 /**
43 * Static methods to generate and differentiate uniform b-Splines.
44 *
45 * @author Michael J. Schnieders
46 * @see <ul>
47 * <li><a href="http://www.springer.com/mathematics/analysis/book/978-0-387-95366-3"
48 * target="_blank"> C. de Boor, A Practical Guide to Splines. (Springer, New York, 2001)
49 * </a>
50 * <li><a href="http://www.wikipedia.org/wiki/B-spline" target="_blank">b-Splines at
51 * Wikipedia</a>
52 * <li><a href="http://mathworld.wolfram.com/B-Spline.html" target="_blank">b-Splines at
53 * MathWorld</a>
54 * </ul>
55 * @since 1.0
56 */
57 public class UniformBSpline {
58
59 /**
60 * Do not allow instantiation of UniformBSpline. All methods are static.
61 */
62 private UniformBSpline() {
63 }
64
65 /**
66 * Generate uniform b-Spline coefficients.
67 *
68 * @param x A double in the range [0.0, 1.0] where 0.5 is over a grid point (0.0 is half-way to
69 * the previous grid point, and 1.0 is half-way to the next grid point).
70 * @param order b-Spline order (degree + 1).
71 * @param coefficients b-Spline coefficients (n coefficients for order n).
72 * @since 1.0
73 */
74 public static void bSpline(final double x, final int order, final double[] coefficients) {
75
76 // Initialization to get to a linear b-Spline (degree 1).
77 coefficients[0] = 1.0 - x;
78 coefficients[1] = x;
79
80 // Apply b-Spline recursion to desired degree.
81 for (int k = 2; k < order; k++) {
82 bSplineRecursion(x, k, coefficients, coefficients);
83 }
84 }
85
86 /**
87 * Generate uniform b-Spline coefficients and their derivatives.
88 *
89 * @param x A double in the range [0.0, 1.0].
90 * @param order b-Spline order (degree + 1).
91 * @param deriveOrder Derivative order. <br> 0 = no derivative. <br> 1 = 1rst derivative. <br>
92 * It must not be greater than the b-Spline degree (order - 1). <br> The method is currently
93 * limited to deriveOrder .LE. 5. <br>
94 * @param coefficients The b-Spline coefficient array of size [order][deriveOrder + 1].
95 * @param work A work array of size [order][order].
96 * @since 1.0
97 */
98 public static void bSplineDerivatives(
99 final double x,
100 final int order,
101 final int deriveOrder,
102 final double[][] coefficients,
103 final double[][] work) {
104
105 assert (deriveOrder <= order - 1 && deriveOrder <= 5);
106
107 for (int k = 0; k < order; k++) {
108 Arrays.fill(work[k], 0.0);
109 }
110
111 // Initialization to get to 2nd order.
112 work[1][0] = 1.0 - x;
113 work[1][1] = x;
114
115 // Perform one pass to get to 3rd order.
116 work[2][0] = 0.5 * (1.0 - x) * work[1][0];
117 work[2][1] = 0.5 * ((x + 1.0) * work[1][0] + (2.0 - x) * work[1][1]);
118 work[2][2] = 0.5 * x * work[1][1];
119
120 // Compute standard B-spline recursion to desired order.
121 for (int k = 3; k < order; k++) {
122 bSplineRecursion(x, k, work[k - 1], work[k]);
123 }
124
125 int o1 = order - 1;
126 // do derivatives
127 try {
128 if (deriveOrder > 0) {
129 int o2 = order - 2;
130 bSplineDiff(work[o2], o1);
131 if (deriveOrder > 1) {
132 int o3 = order - 3;
133 bSplineDiff(work[o3], o2);
134 bSplineDiff(work[o3], o1);
135 if (deriveOrder > 2) {
136 int o4 = order - 4;
137 bSplineDiff(work[o4], o3);
138 bSplineDiff(work[o4], o2);
139 bSplineDiff(work[o4], o1);
140 if (deriveOrder > 3) {
141 int o5 = order - 5;
142 bSplineDiff(work[o5], o4);
143 bSplineDiff(work[o5], o3);
144 bSplineDiff(work[o5], o2);
145 bSplineDiff(work[o5], o1);
146 if (deriveOrder > 4) {
147 int o6 = order - 6;
148 bSplineDiff(work[o6], o5);
149 bSplineDiff(work[o6], o4);
150 bSplineDiff(work[o6], o3);
151 bSplineDiff(work[o6], o2);
152 bSplineDiff(work[o6], o1);
153 if (deriveOrder > 5) {
154 throw new Exception(" Unsupported option: dr_ord > 5");
155 }
156 }
157 }
158 }
159 }
160 }
161 } catch (Exception e) {
162 // Index out of bounds if deriveOrder is too high for order.
163 }
164
165 int deriveOrder1 = deriveOrder + 1;
166 for (int k = 0; k < order; k++) {
167 double[] tk = coefficients[k];
168 try {
169 for (int j = 0; j < deriveOrder1; j++) {
170 tk[j] = work[o1 - j][k];
171 }
172 } catch (Exception e) {
173 // Index out of bounds if deriveOrder is too high for order.
174 }
175 }
176 }
177
178 /**
179 * Differentiate a uniform b-Spline in place.
180 *
181 * @param coefficients B-Spline coefficients.
182 * @param order B-Spline order.
183 * @since 1.0
184 */
185 private static void bSplineDiff(final double[] coefficients, final int order) {
186 final int order1 = order - 1;
187 coefficients[order] = coefficients[order1];
188 for (int i = order1; i > 0; i--) {
189 coefficients[i] = coefficients[i - 1] - coefficients[i];
190 }
191 coefficients[0] = -coefficients[0];
192 }
193
194 /**
195 * Uniform b-Spline recursion.
196 *
197 * @param x A double in the range [0.0, 1.0].
198 * @param order Current b-Spline order.
199 * @param coefficients Current b-Spline coefficients.
200 * @param newCoefficients New b-Spline coefficients for order + 1.
201 * @since 1.0
202 */
203 private static void bSplineRecursion(
204 final double x,
205 final int order,
206 final double[] coefficients,
207 final double[] newCoefficients) {
208 final double div = 1.0 / (double) order;
209 final double k1mw = order + (1.0 - x);
210 final int km1 = order - 1;
211 newCoefficients[order] = div * x * coefficients[km1];
212 for (int i = 1; i < order; i++) {
213 final int kmi = order - i;
214 final int km1i = km1 - i;
215 final double x1 = x + (double) i;
216 final double k1mwi = k1mw - (double) i;
217 newCoefficients[kmi] = div * (x1 * coefficients[km1i] + k1mwi * coefficients[kmi]);
218 }
219 double oneX = 1.0 - x;
220 newCoefficients[0] = div * oneX * coefficients[0];
221 }
222 }