View Javadoc
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 }