1 // ******************************************************************************
2 //
3 // Title: Force Field X.
4 // Description: Force Field X - Software for Molecular Biophysics.
5 // Copyright: Copyright (c) Michael J. Schnieders 2001-2025.
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.integrate;
39
40 import static java.lang.String.format;
41 import static java.lang.System.arraycopy;
42 import static org.apache.commons.math3.util.FastMath.abs;
43 import static org.apache.commons.math3.util.FastMath.max;
44 import static org.apache.commons.math3.util.FastMath.ulp;
45
46 /**
47 * A FunctionDataCurve represents a set of points along a 1-dimensional, analytically integrable
48 * function.
49 *
50 * @author Jacob M. Litman
51 */
52 public abstract class FunctionDataCurve implements DataSet {
53
54 /** Lower bound. */
55 protected double lb;
56 /** Upper bound. */
57 protected double ub;
58 /** Function values. */
59 protected double[] points;
60 /** Input points. */
61 protected double[] x;
62 /** If ends should have 1/2 regular separation. */
63 protected boolean halfWidthEnd;
64
65 /**
66 * Constructor.
67 */
68 public FunctionDataCurve() {
69 lb = 0.0;
70 ub = 0.0;
71 points = new double[0];
72 x = new double[0];
73 halfWidthEnd = false;
74 }
75
76 /**
77 * Checks for equality to +/- 10 ulp.
78 *
79 * @param x1 a double.
80 * @param x2 a double.
81 * @return a boolean.
82 */
83 public static boolean approxEquals(double x1, double x2) {
84 return (approxEquals(x1, x2, 10.0));
85 }
86
87 /**
88 * Compare two doubles to machine precision.
89 *
90 * @param x1 The first point.
91 * @param x2 The second point.
92 * @param ulpMultiple A multiple to apply to the ULP of the larger of the two values.
93 * @return True if the two values are within the specified ULPs of each other.
94 */
95 public static boolean approxEquals(double x1, double x2, double ulpMultiple) {
96 double diff = abs(x1 - x2);
97 // Ulp the larger of the two values.
98 double ulp = ulp(max(abs(x1), abs(x2)));
99 return (diff < (ulp * ulpMultiple));
100 }
101
102 /**
103 * Evaluates the functions analytical integral over the entire range of points.
104 *
105 * @return Exact finite integral
106 */
107 public double analyticalIntegral() {
108 return analyticalIntegral(lowerBound(), upperBound());
109 }
110
111 /**
112 * Evaluates the function's analytical integral over a range.
113 *
114 * @param lb Lower integration bound
115 * @param ub Upper integration bound
116 * @return Exact finite integral of range
117 */
118 public double analyticalIntegral(double lb, double ub) {
119 return integralAt(ub) - integralAt(lb);
120 }
121
122 /** {@inheritDoc} */
123 @Override
124 public double binWidth() {
125 double divisor = halfWidthEnds() ? (double) (points.length - 2) : (double) (points.length - 1);
126 return (ub - lb) / divisor;
127 }
128
129 /**
130 * Evaluates the function at x.
131 *
132 * @param x x
133 * @return f(x)
134 */
135 public abstract double fX(double x);
136
137 /** {@inheritDoc} */
138 @Override
139 public double[] getAllFxPoints() {
140 int nPoints = points.length;
141 double[] retArray = new double[nPoints];
142 arraycopy(points, 0, retArray, 0, nPoints);
143 return retArray;
144 }
145
146 /** {@inheritDoc} */
147 @Override
148 public double getFxPoint(int index) {
149 return points[index];
150 }
151
152 /** {@inheritDoc} */
153 @Override
154 public double[] getX() {
155 double[] copyX = new double[x.length];
156 arraycopy(x, 0, copyX, 0, x.length);
157 return copyX;
158 }
159
160 /** {@inheritDoc} */
161 @Override
162 public boolean halfWidthEnds() {
163 return halfWidthEnd;
164 }
165
166 /**
167 * Analytical integral at a point.
168 *
169 * @param x Point
170 * @return Exact finite integral of 0 to this point
171 */
172 public abstract double integralAt(double x);
173
174 /** {@inheritDoc} */
175 @Override
176 public double lowerBound() {
177 return lb;
178 }
179
180 /** {@inheritDoc} */
181 @Override
182 public int numPoints() {
183 return points.length;
184 }
185
186 /** {@inheritDoc} */
187 @Override
188 public String toString() {
189 StringBuilder sb =
190 new StringBuilder(
191 format(
192 "Function f(x) curve with %d points from lower bound %9.3g and upper bound %9.3g",
193 points.length, lb, ub));
194 if (halfWidthEnd) {
195 sb.append(" and half-width start/end bins");
196 }
197 sb.append(".");
198 return sb.toString();
199 }
200
201 /** {@inheritDoc} */
202 @Override
203 public double upperBound() {
204 return ub;
205 }
206
207 /**
208 * Used to check that x array is composed of equally-spaced points from lb to ub.
209 */
210 protected final void assertXIntegrity() {
211 assert ub > lb;
212 int nX = numPoints();
213 double sep = binWidth();
214 if (halfWidthEnd) {
215 assert x.length == nX;
216 assert lb == x[0];
217 assert ub == x[nX - 1];
218
219 assert approxEquals(x[1], lb + 0.5 * sep);
220 assert approxEquals(x[nX - 2], (ub - 0.5 * sep));
221
222 for (int i = 2; i < (nX - 2); i++) {
223 double target = lb + 0.5 * sep;
224 target += ((i - 1) * sep);
225 assert approxEquals(x[i], target);
226 }
227 } else {
228 for (int i = 0; i < x.length; i++) {
229 assert approxEquals(x[i], x[0] + i * sep);
230 }
231 }
232 }
233 }