1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
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
47
48
49
50
51
52
53
54
55
56
57
58 public class Integrate1DNumeric {
59
60 private static final Logger logger = Logger.getLogger(Integrate1DNumeric.class.getName());
61
62 private static final double ONE_THIRD = (1.0 / 3.0);
63
64 private static final double BOOLE_FACTOR = (2.0 / 45.0);
65
66
67
68
69
70
71
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
88
89
90
91
92
93
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
132
133
134
135
136
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
153
154
155
156
157
158
159
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
209
210
211
212
213
214
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
240
241
242
243
244
245
246
247
248
249
250
251
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);
266
267 switch (side) {
268 case RIGHT -> {
269 values[ub] = width * fX[ub];
270
271
272
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];
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 }
296 case LEFT -> {
297 int shift = halfWide ? 1 : 0;
298 values[lb] = width * fX[lb];
299
300
301
302
303
304 for (int i = lb + 1; i <= ub; i++) {
305
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];
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
333
334
335
336
337
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
361
362
363
364
365
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
382
383
384
385
386
387
388
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
415
416
417
418
419
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
440
441
442
443
444
445
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
462
463
464
465
466
467
468
469
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
487
488
489
490
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
507
508
509
510
511
512
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
550
551
552
553
554
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
571
572
573
574
575
576
577
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
617
618
619
620
621
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
638
639
640
641
642
643
644
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
661
662
663
664
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
677
678
679
680
681
682
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
699
700
701
702
703
704
705
706
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
722
723
724
725
726
727
728
729
730
731
732
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);
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
789
790
791
792
793
794
795
796
797
798 public enum IntegrationSide {
799 LEFT,
800 RIGHT
801 }
802
803
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
823
824 @FunctionalInterface
825 private interface IntegrateWindow {
826
827
828
829
830
831
832
833
834
835
836 double toArea(DataSet data, IntegrationSide side, int lb, int ub);
837 }
838 }