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 }