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  
42  import java.util.logging.Logger;
43  import java.util.stream.IntStream;
44  
45  /**
46   * This program integrates using four methods: rectangular integration, the trapezoidal method,
47   * Simpson's Three Point Integration, and Boole's Five Point Integration. Of these, Simpson's and
48   * Boole's integration methods generally perform the best. Parallelized methods exist and seem to
49   * function, but the performance of a parallelized method should not be necessary, and the
50   * parallelization introduces some complexity and inefficiency.
51   *
52   * <p>Primary methods: integrateData and integrateByBins, using some DataSet, which is often a
53   * DoublesDataSet.
54   *
55   * @author Claire O'Connell
56   * @author Jacob M. Litman
57   */
58  public class Integrate1DNumeric {
59  
60    private static final Logger logger = Logger.getLogger(Integrate1DNumeric.class.getName());
61    /** Constant for used for Simpson's rule. */
62    private static final double ONE_THIRD = (1.0 / 3.0);
63    /** Constant for used for Boole's rule. */
64    private static final double BOOLE_FACTOR = (2.0 / 45.0);
65  
66    /**
67     * Numerically integrates a data set using Boole's rule.
68     *
69     * @param data Data set to integrate
70     * @param side Side to integrate from
71     * @return Area of integral
72     */
73    public static double booles(DataSet data, IntegrationSide side) {
74      double area = 0;
75      int lb = 0;
76      int ub = data.numPoints() - 1;
77      if (data.halfWidthEnds()) {
78        area = trapezoidalEnds(data, side);
79        lb++;
80        ub--;
81      }
82      area += booles(data, side, lb, ub);
83      return area;
84    }
85  
86    /**
87     * Numerically integrates a data set, in bounds lb-ub inclusive, using Boole's rule.
88     *
89     * @param data Data set to integrate
90     * @param side Side to integrate from
91     * @param lb First index to integrate over
92     * @param ub Last index to integrate over
93     * @return Area of integral
94     */
95    public static double booles(DataSet data, IntegrationSide side, int lb, int ub) {
96      double area = 0;
97      double width = data.binWidth();
98      double[] points = data.getAllFxPoints();
99      int increment = 4;
100 
101     int nBins = (ub - lb) / increment;
102     int lowerNeglected = 0;
103     int upperNeglected = 0;
104 
105     switch (side) {
106       case RIGHT -> {
107         for (int i = ub; i > (lb + increment - 1); i -= increment) {
108           area += (7 * points[i] + 32 * points[i - 1] + 12 * points[i - 2]
109               + 32 * points[i - 3] + 7 * points[i - 4]);
110         }
111         lowerNeglected = lb;
112         upperNeglected = ub - (increment * nBins);
113       }
114       case LEFT -> {
115         for (int i = lb; i < (ub - increment + 1); i += increment) {
116           area += (7 * points[i] + 32 * points[i + 1] + 12 * points[i + 2]
117               + 32 * points[i + 3] + 7 * points[i + 4]);
118         }
119         lowerNeglected = lb + (increment * nBins);
120         upperNeglected = ub;
121       }
122     }
123     area *= BOOLE_FACTOR;
124     area *= width;
125 
126     area += finishIntegration(data, side, lowerNeglected, upperNeglected, IntegrationType.BOOLE);
127     return area;
128   }
129 
130   /**
131    * Numerically integrates a data set using Boole's rule. The sequential version is preferred unless
132    * necessary.
133    *
134    * @param data Data set to integrate
135    * @param side Side to integrate from
136    * @return Area of integral
137    */
138   public static double boolesParallel(DataSet data, IntegrationSide side) {
139     double area = 0;
140     int lb = 0;
141     int ub = data.numPoints() - 1;
142     if (data.halfWidthEnds()) {
143       area = trapezoidalEnds(data, side);
144       lb++;
145       ub--;
146     }
147     area += boolesParallel(data, side, lb, ub);
148     return area;
149   }
150 
151   /**
152    * Numerically integrates a data set, in bounds lb-ub inclusive, using Boole's rule. The sequential
153    * version is preferred unless necessary.
154    *
155    * @param data Data set to integrate
156    * @param side Side to integrate from
157    * @param lb First index to integrate over
158    * @param ub Last index to integrate over
159    * @return Area of integral
160    */
161   public static double boolesParallel(DataSet data, IntegrationSide side, int lb, int ub) {
162     double area = 0.0;
163     double width = data.binWidth();
164     double[] points = data.getAllFxPoints();
165     int increment = 4;
166 
167     int nBins = (ub - lb) / increment;
168     int lowerNeglected = 0;
169     int upperNeglected = 0;
170 
171     switch (side) {
172       case RIGHT -> {
173         lowerNeglected = lb;
174         upperNeglected = ub - (increment * nBins);
175         area =
176             IntStream.range(0, nBins).parallel().mapToDouble(
177                 (int i) -> {
178                   int fromPoint = ub - (increment * i);
179                   return (7 * points[fromPoint - 4]
180                       + 32 * points[fromPoint - 3]
181                       + 12 * points[fromPoint - 2]
182                       + 32 * points[fromPoint - 1]
183                       + 7 * points[fromPoint]);
184                 }).sum();
185       }
186       case LEFT -> {
187         lowerNeglected = lb + (increment * nBins);
188         upperNeglected = ub;
189         area = IntStream.range(0, nBins).parallel().mapToDouble(
190             (int i) -> {
191               int fromPoint = lb + (increment * i);
192               return (7 * points[fromPoint]
193                   + 32 * points[fromPoint + 1]
194                   + 12 * points[fromPoint + 2]
195                   + 32 * points[fromPoint + 3]
196                   + 7 * points[fromPoint + 4]);
197             }).sum();
198       }
199     }
200     area *= BOOLE_FACTOR;
201     area *= width;
202 
203     area += finishIntegration(data, side, lowerNeglected, upperNeglected, IntegrationType.BOOLE);
204     return area;
205   }
206 
207   /**
208    * Generates a set of points along x.
209    *
210    * @param lb Beginning value, inclusive
211    * @param ub Ending value, inclusive
212    * @param nPoints Total number of points
213    * @param halfWidthEnds If ends should have 1/2 regular separation
214    * @return an array of {@link double} objects.
215    */
216   public static double[] generateXPoints(double lb, double ub, int nPoints, boolean halfWidthEnds) {
217     if (lb >= ub) {
218       throw new IllegalArgumentException("ub must be greater than lb");
219     }
220     double[] points = new double[nPoints];
221     int sepDivisor = halfWidthEnds ? nPoints - 2 : nPoints - 1;
222     double sep = (ub - lb) / ((double) sepDivisor);
223 
224     if (halfWidthEnds) {
225       points[0] = lb;
226       points[nPoints - 1] = ub;
227       for (int i = 1; i < nPoints - 1; i++) {
228         points[i] = lb + ((double) i) * sep - 0.5 * sep;
229       }
230     } else {
231       for (int i = 0; i < nPoints; i++) {
232         points[i] = lb + ((double) i) * sep;
233       }
234     }
235     return points;
236   }
237 
238   /**
239    * Returns the contribution of each bin to the overall integral as an array; will be most accurate
240    * at break-points for the integration type. Overall integral is the sum of the array of doubles,
241    * plus or minus minor rounding errors.
242    *
243    * <p>If N is a breakpoint for Boole's rule: N+1 is a trapezoid from N to N+1. N+2 is Simpson's
244    * from N to N+2, minus the prior trapezoid (N+1). N+3 is 4-point integration from N to N+3, minus
245    * the N to N+2 parabola. N+4 is the full Boole's Rule from N to N+4, minus the N to N+3 4-point
246    * integral.
247    *
248    * @param data Data to integrate
249    * @param side Side to integrate from
250    * @param maxType Maximum rule to be used
251    * @return Per-bin contributions to integral.
252    */
253   public static double[] integrateByBins(
254       DataSet data, IntegrationSide side, IntegrationType maxType) {
255     boolean halfWide = data.halfWidthEnds();
256     double[] fX = data.getAllFxPoints();
257     int numPoints = data.numPoints();
258     double width = data.binWidth();
259     double[] values = new double[numPoints];
260 
261     int lb = halfWide ? 1 : 0;
262     int ub = halfWide ? numPoints - 2 : numPoints - 1;
263 
264     int increment = maxType.pointsNeeded() - 1;
265     increment = Math.max(1, increment); // Deal w/ rectangle integration.
266 
267     switch (side) {
268       case RIGHT -> {
269         values[ub] = width * fX[ub]; // Begin with rectangle.
270         /*
271           For each bin, its contribution to this sub-window's value is the integral from (start to
272           bin) minus the integral from (start to prior bin).
273          */
274         for (int i = ub - 1; i >= lb; i--) {
275           int fromUB = ub - i - 1;
276           fromUB /= increment;
277           fromUB *= increment;
278           int lastBegin = ub - fromUB;
279           double val = finishIntegration(data, side, i, lastBegin, maxType);
280           val -= finishIntegration(data, side, i + 1, lastBegin, maxType);
281           values[i] = val;
282         }
283         values[ub - 1] -= values[ub]; // Remove double-counting at start.
284         if (halfWide) {
285           if (maxType == IntegrationType.RECTANGULAR) {
286             values[1] += 0.5 * width * fX[1];
287             values[numPoints - 1] = (0.5 * width * fX[numPoints - 1]);
288           } else {
289             values[0] = 0.25 * width * fX[0];
290             values[1] += 0.25 * width * fX[1];
291             values[numPoints - 2] += (0.25 * width * fX[numPoints - 2]);
292             values[numPoints - 1] = (0.25 * width * fX[numPoints - 1]);
293           }
294         }
295       } // LEFT
296       case LEFT -> {
297         int shift = halfWide ? 1 : 0;
298         values[lb] = width * fX[lb]; // Begin with rectangle.
299 
300         /*
301           For each bin, its contribution to this sub-window's value is the integral from (start to
302           bin) minus the integral from (start to prior bin).
303          */
304         for (int i = lb + 1; i <= ub; i++) {
305           // Remove remainder via division-and-multiplication.
306           int lastBegin = ((i - 1 - shift) / increment);
307           lastBegin *= increment;
308           lastBegin += shift;
309           double val = finishIntegration(data, side, lastBegin, i, maxType);
310           val -= finishIntegration(data, side, lastBegin, i - 1, maxType);
311           values[i] = val;
312         }
313         values[lb + 1] -= values[lb]; // Remove double-counting at start.
314         if (halfWide) {
315           if (maxType == IntegrationType.RECTANGULAR) {
316             values[0] = 0.5 * width * fX[0];
317             values[numPoints - 2] += (0.5 * width * fX[numPoints - 2]);
318           } else {
319             values[0] = 0.25 * width * fX[0];
320             values[1] += 0.25 * width * fX[1];
321             values[numPoints - 2] += (0.25 * width * fX[numPoints - 2]);
322             values[numPoints - 1] = (0.25 * width * fX[numPoints - 1]);
323           }
324         }
325       }
326     }
327 
328     return values;
329   }
330 
331   /**
332    * Generic caller for 1D integration schemes given an IntegrationType.
333    *
334    * @param data To integrate
335    * @param side Integrate from side
336    * @param type Scheme to use
337    * @return Numeric integral
338    */
339   public static double integrateData(DataSet data, IntegrationSide side, IntegrationType type) {
340     switch (type) {
341       case RECTANGULAR -> {
342         return rectangular(data, side);
343       }
344       case TRAPEZOIDAL -> {
345         return trapezoidal(data, side);
346       }
347       case SIMPSONS -> {
348         return simpsons(data, side);
349       }
350       case BOOLE -> {
351         return booles(data, side);
352       }
353     }
354     logger.warning(
355         format(" Integration type %s not recognized! Defaulting to Simpson's integration", type));
356     return simpsons(data, side);
357   }
358 
359   /**
360    * Numerically integrates a data set using rectangular integration. Not recommended; preferably use
361    * at least trapezoidal integration.
362    *
363    * @param data Data set to integrate
364    * @param side Side to integrate from
365    * @return Area of integral
366    */
367   public static double rectangular(DataSet data, IntegrationSide side) {
368     double area = 0;
369     int lb = 0;
370     int ub = data.numPoints() - 1;
371     if (data.halfWidthEnds()) {
372       area = rectangularEnds(data, side);
373       ++lb;
374       --ub;
375     }
376     area += rectangular(data, side, lb, ub);
377     return area;
378   }
379 
380   /**
381    * Numerically integrates a data set, in bounds lb-ub inclusive, using rectangular integration. Not
382    * recommended; preferably use at least trapezoidal integration.
383    *
384    * @param data Data set to integrate
385    * @param side Side to integrate from
386    * @param lb First index to integrate over
387    * @param ub Last index to integrate over
388    * @return Area of integral
389    */
390   public static double rectangular(DataSet data, IntegrationSide side, int lb, int ub) {
391     double area = 0;
392     double width = data.binWidth();
393     double[] points = data.getAllFxPoints();
394     assert ub > lb;
395     assert ub < points.length;
396 
397     switch (side) {
398       case RIGHT -> {
399         for (int i = ub; i > lb; i--) {
400           area += (width * points[i]);
401         }
402       }
403       case LEFT -> {
404         for (int i = lb; i < ub; i++) {
405           area += width * points[i];
406         }
407       }
408     }
409 
410     return area;
411   }
412 
413   /**
414    * Treats half-width bins at the ends of a DataSet using rectangular integration. Not recommended,
415    * preferably use trapezoidal integration.
416    *
417    * @param data Data set to integrate
418    * @param side Side to integrate from
419    * @return Estimate of area in half-width start/end bins.
420    */
421   public static double rectangularEnds(DataSet data, IntegrationSide side) {
422     double width = 0.5 * data.binWidth();
423     double area = 0;
424     int nPoints = data.numPoints();
425     switch (side) {
426       case LEFT -> {
427         area = data.getFxPoint(0) * width;
428         area += (data.getFxPoint(nPoints - 2) * width);
429       }
430       case RIGHT -> {
431         area = data.getFxPoint(1) * width;
432         area += (data.getFxPoint(nPoints - 1) * width);
433       }
434     }
435     return area;
436   }
437 
438   /**
439    * Numerically integrates a data set using parallelized rectangular integration. Not recommended;
440    * preferably use at least trapezoidal integration, and avoid parallelized versions unless
441    * necessary.
442    *
443    * @param data Data set to integrate
444    * @param side Side to integrate from
445    * @return Area of integral
446    */
447   public static double rectangularParallel(DataSet data, IntegrationSide side) {
448     double area = 0;
449     int lb = 0;
450     int ub = data.numPoints() - 1;
451     if (data.halfWidthEnds()) {
452       area = rectangularEnds(data, side);
453       ++lb;
454       --ub;
455     }
456     area += rectangularParallel(data, side, lb, ub);
457     return area;
458   }
459 
460   /**
461    * Numerically integrates a data set, in bounds lb-ub inclusive, using rectangular integration. Not
462    * recommended; preferably use at least trapezoidal integration. Also, prefer parallelized versions
463    * unless necessary.
464    *
465    * @param data Data set to integrate
466    * @param side Side to integrate from
467    * @param lb First index to integrate over
468    * @param ub Last index to integrate over
469    * @return Area of integral
470    */
471   public static double rectangularParallel(DataSet data, IntegrationSide side, int lb, int ub) {
472     double width = data.binWidth();
473     double[] points = data.getAllFxPoints();
474     assert ub > lb;
475     assert ub < points.length;
476 
477     if (side == IntegrationSide.RIGHT) {
478       ++ub;
479       ++lb;
480     }
481     return IntStream.range(lb, ub).parallel().mapToDouble(
482         (int i) -> points[i] * width).sum();
483   }
484 
485   /**
486    * Numerically integrates a data set using Simpson's rule.
487    *
488    * @param data Data set to integrate
489    * @param side Side to integrate from
490    * @return Area of integral
491    */
492   public static double simpsons(DataSet data, IntegrationSide side) {
493     double area = 0;
494     int lb = 0;
495     int ub = data.numPoints() - 1;
496     if (data.halfWidthEnds()) {
497       area = trapezoidalEnds(data, side);
498       lb++;
499       ub--;
500     }
501     area += simpsons(data, side, lb, ub);
502     return area;
503   }
504 
505   /**
506    * Numerically integrates a data set, in bounds lb-ub inclusive, using Simpson's rule.
507    *
508    * @param data Data set to integrate
509    * @param side Side to integrate from
510    * @param lb First index to integrate over
511    * @param ub Last index to integrate over
512    * @return Area of integral
513    */
514   public static double simpsons(DataSet data, IntegrationSide side, int lb, int ub) {
515     double area = 0;
516     double width = data.binWidth();
517     double[] points = data.getAllFxPoints();
518     int increment = 2;
519 
520     int nBins = (ub - lb) / increment;
521     int lowerNeglected = 0;
522     int upperNeglected = 0;
523 
524     switch (side) {
525       case RIGHT -> {
526         for (int i = ub; i > (lb + increment - 1); i -= increment) {
527           area += points[i] + (4 * points[i - 1]) + points[i - 2];
528         }
529         lowerNeglected = lb;
530         upperNeglected = ub - (increment * nBins);
531       }
532       case LEFT -> {
533         for (int i = lb; i < (ub - increment + 1); i += increment) {
534           area += points[i] + (4 * points[i + 1]) + points[i + 2];
535         }
536         lowerNeglected = lb + (increment * nBins);
537         upperNeglected = ub;
538       }
539     }
540 
541     area *= ONE_THIRD;
542     area *= width;
543 
544     area += finishIntegration(data, side, lowerNeglected, upperNeglected, IntegrationType.SIMPSONS);
545     return area;
546   }
547 
548   /**
549    * Numerically integrates a data set using Boole's rule. The sequential version is preferred unless
550    * necessary.
551    *
552    * @param data Data set to integrate
553    * @param side Side to integrate from
554    * @return Area of integral
555    */
556   public static double simpsonsParallel(DataSet data, IntegrationSide side) {
557     double area = 0;
558     int lb = 0;
559     int ub = data.numPoints() - 1;
560     if (data.halfWidthEnds()) {
561       area = trapezoidalEnds(data, side);
562       lb++;
563       ub--;
564     }
565     area += simpsonsParallel(data, side, lb, ub);
566     return area;
567   }
568 
569   /**
570    * Numerically integrates a data set, in bounds lb-ub inclusive, using Simpson's rule. The
571    * sequential version is preferred unless necessary.
572    *
573    * @param data Data set to integrate
574    * @param side Side to integrate from
575    * @param lb First index to integrate over
576    * @param ub Last index to integrate over
577    * @return Area of integral
578    */
579   public static double simpsonsParallel(DataSet data, IntegrationSide side, int lb, int ub) {
580     double area = 0;
581     double width = data.binWidth();
582     double[] points = data.getAllFxPoints();
583     int increment = 2;
584     int nBins = (ub - lb) / increment;
585     int lowerNeglected = 0;
586     int upperNeglected = 0;
587 
588     switch (side) {
589       case RIGHT -> {
590         lowerNeglected = lb;
591         upperNeglected = ub - (increment * nBins);
592         area = IntStream.range(0, nBins).parallel().mapToDouble(
593             (int i) -> {
594               int fromPoint = ub - (increment * i);
595               return (points[fromPoint - 2] + 4 * points[fromPoint - 1] + points[fromPoint]);
596             }).sum();
597       }
598       case LEFT -> {
599         area = IntStream.range(0, nBins).parallel().mapToDouble(
600             (int i) -> {
601               int fromPoint = lb + (increment * i);
602               return (points[fromPoint] + 4 * points[fromPoint + 1] + points[fromPoint + 2]);
603             }).sum();
604         lowerNeglected = lb + (increment * nBins);
605         upperNeglected = ub;
606       }
607     }
608 
609     area *= ONE_THIRD;
610     area *= width;
611     area += finishIntegration(data, side, lowerNeglected, upperNeglected, IntegrationType.SIMPSONS);
612     return area;
613   }
614 
615   /**
616    * Numerically integrates a data set using trapezoidal integration. For most data sets, Simpson's
617    * rule or Boole's rule will out-perform trapezoidal integration.
618    *
619    * @param data Data set to integrate
620    * @param side Side to integrate from
621    * @return Area of integral
622    */
623   public static double trapezoidal(DataSet data, IntegrationSide side) {
624     double area = 0;
625     int lb = 0;
626     int ub = data.numPoints() - 1;
627     if (data.halfWidthEnds()) {
628       area = trapezoidalEnds(data, side);
629       lb++;
630       ub--;
631     }
632     area += trapezoidal(data, side, lb, ub);
633     return area;
634   }
635 
636   /**
637    * Numerically integrates a data set, in bounds lb-ub inclusive, using trapezoidal integration. For
638    * most data sets, Simpson's rule or Boole's rule will out-perform trapezoidal integration.
639    *
640    * @param data Data set to integrate
641    * @param side Side to integrate from
642    * @param lb First index to integrate over
643    * @param ub Last index to integrate over
644    * @return Area of integral
645    */
646   public static double trapezoidal(DataSet data, IntegrationSide side, int lb, int ub) {
647     double width = data.binWidth();
648     double[] points = data.getAllFxPoints();
649 
650     double area = 0.5 * points[lb];
651     area += 0.5 * points[ub];
652     for (int i = lb + 1; i < ub; i++) {
653       area += points[i];
654     }
655     area *= width;
656     return area;
657   }
658 
659   /**
660    * Treats half-width bins at the ends of a DataSet using trapezoidal integration.
661    *
662    * @param data Data set to integrate
663    * @param side Side to integrate from
664    * @return Estimate of area in half-width start/end bins.
665    */
666   public static double trapezoidalEnds(DataSet data, IntegrationSide side) {
667     double width = 0.5 * data.binWidth();
668     int nPts = data.numPoints();
669     double area = data.getFxPoint(0) + data.getFxPoint(1)
670         + data.getFxPoint(nPts - 2) + data.getFxPoint(nPts - 1);
671     area *= (0.5 * width);
672     return area;
673   }
674 
675   /**
676    * Numerically integrates a data set using trapezoidal integration. For most data sets, Simpson's
677    * rule or Boole's rule will out-perform trapezoidal integration. Prefer use of the sequential
678    * version unless necessary.
679    *
680    * @param data Data set to integrate
681    * @param side Side to integrate from
682    * @return Area of integral
683    */
684   public static double trapezoidalParallel(DataSet data, IntegrationSide side) {
685     double area = 0;
686     int lb = 0;
687     int ub = data.numPoints() - 1;
688     if (data.halfWidthEnds()) {
689       area = trapezoidalEnds(data, side);
690       lb++;
691       ub--;
692     }
693     area += trapezoidalParallel(data, side, lb, ub);
694     return area;
695   }
696 
697   /**
698    * Numerically integrates a data set, in bounds lb-ub inclusive, using trapezoidal integration. For
699    * most data sets, Simpson's rule or Boole's rule will out-perform trapezoidal integration. Prefer
700    * use of the sequential version unless necessary.
701    *
702    * @param data Data set to integrate
703    * @param side Side to integrate from
704    * @param lb First index to integrate over
705    * @param ub Last index to integrate over
706    * @return Area of integral
707    */
708   public static double trapezoidalParallel(DataSet data, IntegrationSide side, int lb, int ub) {
709     double width = data.binWidth();
710     double[] points = data.getAllFxPoints();
711 
712     double area = 0.5 * points[lb];
713     area += 0.5 * points[ub];
714     area += IntStream.range(lb + 1, ub).parallel().mapToDouble(
715         (int i) -> points[i]).sum();
716     area *= width;
717     return area;
718   }
719 
720   /**
721    * Integrates the remaining points after higher-order integration rules cannot evenly fit over the
722    * remaining data. For example, Boole's rule requires 5 points; if a data set has 7 points, it will
723    * integrate the first 5, but not handle the remaining 2; this method will, apply the highest-order
724    * integration rule that the remaining points (including the last one used previously) permit. In
725    * that case, it would be Simpson's rule.
726    *
727    * @param data Finish numerical integration on
728    * @param side Integration side
729    * @param lb Index of lower bound to complete integration on
730    * @param ub Index of upper bound to complete integration on
731    * @param type Highest-order rule permitted to finish integration
732    * @return Area of un-integrated region
733    */
734   private static double finishIntegration(
735       DataSet data, IntegrationSide side, int lb, int ub, IntegrationType type) {
736     int totPoints = (ub - lb);
737 
738     int perBin = type.pointsNeeded();
739     int increment = perBin - 1;
740     increment = Math.max(1, increment); // Needed for rectangular integration
741 
742     int nBins = totPoints / increment;
743     int remainder = totPoints % increment;
744 
745     IntegrateWindow intMode = switch (type) {
746       case BOOLE -> Integrate1DNumeric::booles;
747       case SIMPSONS -> Integrate1DNumeric::simpsons;
748       case RECTANGULAR -> Integrate1DNumeric::rectangular;
749       case TRAPEZOIDAL -> Integrate1DNumeric::trapezoidal;
750     };
751 
752     double area = 0.0;
753 
754     switch (side) {
755       case RIGHT -> {
756         for (int i = ub; i > (lb - 1 + increment); i -= increment) {
757           area += intMode.toArea(data, side, (i - increment), i);
758         }
759         ub -= (nBins * increment);
760       }
761       case LEFT -> {
762         for (int i = lb; i < (ub + 1 - increment); i += increment) {
763           area += intMode.toArea(data, side, i, (i + increment));
764         }
765         lb += (nBins * increment);
766       }
767     }
768 
769     assert remainder == (ub - lb);
770     switch (remainder) {
771       case 0:
772         break;
773       case 1:
774         area += trapezoidal(data, side, lb, ub);
775         break;
776       case 2:
777       case 3:
778         area += simpsons(data, side, lb, ub);
779         break;
780       case 4:
781       default:
782         throw new IllegalArgumentException("This should not be currently possible.");
783     }
784     return area;
785   }
786 
787   /**
788    * Left vs right-hand integration; left-hand integration will start from the first available point,
789    * run right as far as possible, and then clean up any remaining points using finishIntegration,
790    * while right-hand integration will start from the last available point, run left as far as
791    * possible, and then clean up any remaining points using finishIntegration.
792    *
793    * <p>Values: LEFT, RIGHT.
794    *
795    * <p>Not meaningful for trapezoidal integration; there will never be any unused points, and the
796    * method is symmetrical.
797    */
798   public enum IntegrationSide {
799     LEFT,
800     RIGHT
801   }
802 
803   /** Enumeration of implemented integration methods, and the number of points required by them. */
804   public enum IntegrationType {
805     RECTANGULAR(1),
806     TRAPEZOIDAL(2),
807     SIMPSONS(3),
808     BOOLE(5);
809 
810     private final int pointsNeeded;
811 
812     IntegrationType(int points) {
813       pointsNeeded = points;
814     }
815 
816     public final int pointsNeeded() {
817       return pointsNeeded;
818     }
819   }
820 
821   /**
822    * Functional interface used to select an integration method (primarily for finishIntegration).
823    */
824   @FunctionalInterface
825   private interface IntegrateWindow {
826 
827     /**
828      * Numerically integrates a range of x given f(x).
829      *
830      * @param data x and f(x) data to integrate
831      * @param side Side to integrate from
832      * @param lb Lower bound of x[] to integrate
833      * @param ub Upper bound of x[] to integrate
834      * @return Area of integral f(x) dx, from point lb to point ub
835      */
836     double toArea(DataSet data, IntegrationSide side, int lb, int ub);
837   }
838 }