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-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.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     * Checks for equality to +/- 10 ulp.
67     *
68     * @param x1 a double.
69     * @param x2 a double.
70     * @return a boolean.
71     */
72    public static boolean approxEquals(double x1, double x2) {
73      return (approxEquals(x1, x2, 10.0));
74    }
75  
76    /**
77     * Compare two doubles to machine precision.
78     *
79     * @param x1 The first point.
80     * @param x2 The second point.
81     * @param ulpMultiple A multiple to apply to the ULP of the larger of the two values.
82     * @return True if the two values are within the specified ULPs of each other.
83     */
84    public static boolean approxEquals(double x1, double x2, double ulpMultiple) {
85      double diff = abs(x1 - x2);
86      // Ulp the larger of the two values.
87      double ulp = ulp(max(abs(x1), abs(x2)));
88      return (diff < (ulp * ulpMultiple));
89    }
90  
91    /**
92     * Evaluates the functions analytical integral over the entire range of points.
93     *
94     * @return Exact finite integral
95     */
96    public double analyticalIntegral() {
97      return analyticalIntegral(lowerBound(), upperBound());
98    }
99  
100   /**
101    * Evaluates the function's analytical integral over a range.
102    *
103    * @param lb Lower integration bound
104    * @param ub Upper integration bound
105    * @return Exact finite integral of range
106    */
107   public double analyticalIntegral(double lb, double ub) {
108     return integralAt(ub) - integralAt(lb);
109   }
110 
111   /** {@inheritDoc} */
112   @Override
113   public double binWidth() {
114     double divisor = halfWidthEnds() ? (double) (points.length - 2) : (double) (points.length - 1);
115     return (ub - lb) / divisor;
116   }
117 
118   /**
119    * Evaluates the function at x.
120    *
121    * @param x x
122    * @return f(x)
123    */
124   public abstract double fX(double x);
125 
126   /** {@inheritDoc} */
127   @Override
128   public double[] getAllFxPoints() {
129     int nPoints = points.length;
130     double[] retArray = new double[nPoints];
131     arraycopy(points, 0, retArray, 0, nPoints);
132     return retArray;
133   }
134 
135   /** {@inheritDoc} */
136   @Override
137   public double getFxPoint(int index) {
138     return points[index];
139   }
140 
141   /** {@inheritDoc} */
142   @Override
143   public double[] getX() {
144     double[] copyX = new double[x.length];
145     arraycopy(x, 0, copyX, 0, x.length);
146     return copyX;
147   }
148 
149   /** {@inheritDoc} */
150   @Override
151   public boolean halfWidthEnds() {
152     return halfWidthEnd;
153   }
154 
155   /**
156    * Analytical integral at a point.
157    *
158    * @param x Point
159    * @return Exact finite integral of 0 to this point
160    */
161   public abstract double integralAt(double x);
162 
163   /** {@inheritDoc} */
164   @Override
165   public double lowerBound() {
166     return lb;
167   }
168 
169   /** {@inheritDoc} */
170   @Override
171   public int numPoints() {
172     return points.length;
173   }
174 
175   /** {@inheritDoc} */
176   @Override
177   public String toString() {
178     StringBuilder sb =
179         new StringBuilder(
180             format(
181                 "Function f(x) curve with %d points from lower bound %9.3g and upper bound %9.3g",
182                 points.length, lb, ub));
183     if (halfWidthEnd) {
184       sb.append(" and half-width start/end bins");
185     }
186     sb.append(".");
187     return sb.toString();
188   }
189 
190   /** {@inheritDoc} */
191   @Override
192   public double upperBound() {
193     return ub;
194   }
195 
196   /**
197    * Used to check that x array is composed of equally-spaced points from lb to ub.
198    */
199   protected final void assertXIntegrity() {
200     assert ub > lb;
201     int nX = numPoints();
202     double sep = binWidth();
203     if (halfWidthEnd) {
204       assert x.length == nX;
205       assert lb == x[0];
206       assert ub == x[nX - 1];
207 
208       assert approxEquals(x[1], lb + 0.5 * sep);
209       assert approxEquals(x[nX - 2], (ub - 0.5 * sep));
210 
211       for (int i = 2; i < (nX - 2); i++) {
212         double target = lb + 0.5 * sep;
213         target += ((i - 1) * sep);
214         assert approxEquals(x[i], target);
215       }
216     } else {
217       for (int i = 0; i < x.length; i++) {
218         assert approxEquals(x[i], x[0] + i * sep);
219       }
220     }
221   }
222 }