Unofficial Dub Package Documentation 
1 /** Authors:    Lars Tandle Kyllingstad
2     Copyright:  Copyright (c) 2010–2011, Lars T. Kyllingstad. All rights reserved.
3     License:    Boost License 1.0
4 */
5 module scid.internal.calculus.integrate_qng;
6 
7 
8 import std.conv;
9 import std.exception;
10 import std.math;
11 import std.numeric;
12 import std.traits;
13 
14 import scid.calculus: IntegrationException;
15 import scid.types;
16 import scid.util;
17 
18 version(unittest) import scid.core.testing;
19 
20 
21 
22 
23 Result!Real qng(Func, Real = ReturnType!Func)
24     (Func f, Real a, Real b, Real epsRel, Real epsAbs)
25 {
26     alias FPTemporary!Real T;
27 
28     enum epsRelMin = 50*T.epsilon;
29     enforce(
30         epsAbs > 0 || epsRel > epsRelMin,
31         text("Invalid accuracy request (must have epsAbs > 0 or epsRel > ",
32             epsRelMin, ")"));
33 
34 
35     // This function is used to calculate the absolute error.
36     static T errorCalc(T estimate, T absResult, T variance)
37     {
38         T error = estimate;
39         if (variance != 0 && error != 0)
40             error = variance * fmin(1, (200*error/variance)^^1.5);
41         if (absResult > T.min_normal/(50*T.epsilon))
42             error = fmax(error, 50*T.epsilon*absResult);
43         return error;
44     }
45 
46 
47     // Gauss-Kronrod-Patterson quadrature coefficients,
48     // calculated by L. W. Fullerton, Bell Labs, Nov. 1981.
49     static immutable T[5] x1 = staticArrayOf!T(
50         0.9739065285_1717172007_7964012084_452L,
51         0.8650633666_8898451073_2096688423_493L,
52         0.6794095682_9902440623_4327365114_874L,
53         0.4333953941_2924719079_9265943165_784L,
54         0.1488743389_8163121088_4826001129_720L);
55     static immutable T[5] w10 = staticArrayOf!T(
56         0.0666713443_0868813759_3568809893_332L,
57         0.1494513491_5058059314_5776339657_697L,
58         0.2190863625_1598204399_5534934228_163L,
59         0.2692667193_0999635509_1226921569_469L,
60         0.2955242247_1475287017_3892994651_338L);
61 
62 
63     static immutable T[5] x2 = staticArrayOf!T(
64         0.9956571630_2580808073_5527280689_003L,
65         0.9301574913_5570822600_1207180059_508L,
66         0.7808177265_8641689706_3717578345_042L,
67         0.5627571346_6860468333_9000099272_694L,
68         0.2943928627_0146019813_1126603103_866L);
69     static immutable T[5] w21a = staticArrayOf!T(
70         0.0325581623_0796472747_8818972459_390L,
71         0.0750396748_1091995276_7043140916_190L,
72         0.1093871588_0229764189_9210590325_805L,
73         0.1347092173_1147332592_8054001771_707L,
74         0.1477391049_0133849137_4841515972_068L);
75     static immutable T[6] w21b = staticArrayOf!T(
76         0.0116946388_6737187427_8064396062_192L,
77         0.0547558965_7435199603_1381300244_580L,
78         0.0931254545_8369760553_5065465083_366L,
79         0.1234919762_6206585107_7958109831_074L,
80         0.1427759385_7706008079_7094273138_717L,
81         0.1494455540_0291690566_4936468389_821L);
82 
83     static immutable T[11] x3 = staticArrayOf!T(
84         0.9993333609_0193208139_4099323919_911L,
85         0.9874334029_0808886979_5961478381_209L,
86         0.9548079348_1426629925_7919200290_473L,
87         0.9001486957_4832829362_5099494069_092L,
88         0.8251983149_8311415084_7066732588_520L,
89         0.7321483889_8930498261_2354848755_461L,
90         0.6228479705_3772523864_1159120344_323L,
91         0.4994795740_7105649995_2214885499_755L,
92         0.3649016613_4658076804_3989548502_644L,
93         0.2222549197_7660129649_8260928066_212L,
94         0.0746506174_6138332204_3914435796_506L);
95     static immutable T[10] w43a = staticArrayOf!T(
96         0.0162967342_8966656492_4281974617_663L,
97         0.0375228761_2086950146_1613795898_115L,
98         0.0546949020_5825544214_7212685465_005L,
99         0.0673554146_0947808607_5553166302_174L,
100         0.0738701996_3239395343_2140695251_367L,
101         0.0057685560_5976979618_4184327908_655L,
102         0.0273718905_9324884208_1276069289_151L,
103         0.0465608269_1042883074_3339154433_824L,
104         0.0617449952_0144256449_6240336030_883L,
105         0.0713872672_6869339776_8559114425_516L);
106     static immutable T[12] w43b = staticArrayOf!T(
107         0.0018444776_4021241410_0389106552_965L,
108         0.0107986895_8589165174_0465406741_293L,
109         0.0218953638_6779542810_2523123075_149L,
110         0.0325974639_7534568944_3882222526_137L,
111         0.0421631379_3519181184_7627924327_955L,
112         0.0507419396_0018457778_0189020092_084L,
113         0.0583793955_4261924837_5475369330_206L,
114         0.0647464049_5144588554_4689259517_511L,
115         0.0695661979_1235648452_8633315038_405L,
116         0.0728244414_7183320815_0939535192_842L,
117         0.0745077510_1417511827_3571813842_889L,
118         0.0747221475_1740300559_4425168280_423L);
119 
120     static immutable T[22] x4 = staticArrayOf!T(
121         0.9999029772_6272923449_0529830591_582L,
122         0.9979898959_8667874542_7496322365_960L,
123         0.9921754978_6068722280_8523352251_425L,
124         0.9813581635_7271277357_1916941623_894L,
125         0.9650576238_5838461912_8284110607_926L,
126         0.9431676131_3367059681_6416634507_426L,
127         0.9158064146_8550720959_1826430720_050L,
128         0.8832216577_7131650137_2117548744_163L,
129         0.8457107484_6241566660_5902011504_855L,
130         0.8035576580_3523098278_8739474980_964L,
131         0.7570057306_8549555832_8942793432_020L,
132         0.7062732097_8732181982_4094274740_840L,
133         0.6515894665_0117792253_4422205016_736L,
134         0.5932233740_5796108887_5273770349_144L,
135         0.5314936059_7083193228_5268948562_671L,
136         0.4667636230_4202284487_1966781659_270L,
137         0.3994248478_5921880473_2101665817_923L,
138         0.3298748771_0618828826_5053371824_597L,
139         0.2585035592_0216155180_2280975429_025L,
140         0.1856953965_6834665201_5917141167_606L,
141         0.1118422131_7990746817_2398359241_362L,
142         0.0373521233_9461987081_4998165437_704L);
143     static immutable T[21] w87a = staticArrayOf!T(
144         0.0081483773_8414917290_0002878448_190L,
145         0.0187614382_0156282224_3935059003_794L,
146         0.0273474510_5005228616_1582829741_283L,
147         0.0336777073_1163793004_6581056957_588L,
148         0.0369350998_2042790761_4589586742_499L,
149         0.0028848724_3021153050_1334156248_695L,
150         0.0136859460_2271270188_8950035273_128L,
151         0.0232804135_0288831112_3409291030_404L,
152         0.0308724976_1171335867_5466394126_442L,
153         0.0356936336_3941877071_9351355457_044L,
154         0.0009152833_4520224136_0843392549_948L,
155         0.0053992802_1930047136_7738743391_053L,
156         0.0109476796_0111893113_4327826856_808L,
157         0.0162987316_9678733526_2665703223_280L,
158         0.0210815688_8920383511_2433060188_190L,
159         0.0253709697_6925382724_3467999831_710L,
160         0.0291896977_5647575250_1446154084_920L,
161         0.0323732024_6720278968_5788194889_595L,
162         0.0347830989_5036514275_0781997949_596L,
163         0.0364122207_3135178756_2801163687_577L,
164         0.0372538755_0304770853_9592001191_226L);
165     static immutable T[23] w87b = staticArrayOf!T(
166         0.0002741455_6376207235_0016527092_881L,
167         0.0018071241_5505794294_8341311753_254L,
168         0.0040968692_8275916486_4458070683_480L,
169         0.0067582900_5184737869_9816577897_424L,
170         0.0095499576_7220164653_6053581325_377L,
171         0.0123294476_5224485369_4626639963_780L,
172         0.0150104473_4638895237_6697286041_943L,
173         0.0175489679_8624319109_9665352925_900L,
174         0.0199380377_8644088820_2278192730_714L,
175         0.0221949359_6101228679_6332102959_499L,
176         0.0243391471_2600080547_0360647041_454L,
177         0.0263745054_1483920724_1503786552_615L,
178         0.0282869107_8877120065_9968002987_960L,
179         0.0300525811_2809269532_2521110347_341L,
180         0.0316467513_7143992940_4586051078_883L,
181         0.0330504134_1997850329_0785944862_689L,
182         0.0342550997_0422606178_7082821046_821L,
183         0.0352624126_6015668103_3782717998_428L,
184         0.0360769896_2288870118_5500318003_895L,
185         0.0366986044_9845609449_8018047441_094L,
186         0.0371205492_6983257611_4119958413_599L,
187         0.0373342287_5193504032_1235449094_698L,
188         0.0373610737_6267902341_0321241766_599L);
189 
190 
191     // Some useful quantities.
192     immutable T halfLength = (b - a) / 2;
193     immutable T absHalfLength = abs(halfLength);
194 
195     immutable T center = (a + b) / 2;
196     immutable T fCenter = f(center);
197 
198     // In these arrays we save the result of function evaluations,
199     // so they can be reused for higher-order rules.
200     T[5] fSave1, fSave2, fSave3, fSave4;
201     T[21] fSave;
202 
203 
204     // Compute the integral using the 10- and 21-point formulae.
205     T result10 = 0;
206     T result21 = w21b[5] * fCenter;
207     T resultAbs = w21b[5] * abs(fCenter);
208 
209     foreach (k; 0 .. 5)
210     {
211         immutable x = halfLength * x1[k];
212 
213         immutable fValue1 = f(center + x);
214         immutable fValue2 = f(center - x);
215         immutable fValue = fValue1 + fValue2;
216 
217         result10 += w10[k] * fValue;
218         result21 += w21a[k] * fValue;
219         resultAbs += w21a[k] * (abs(fValue1) + abs(fValue2));
220 
221         fSave1[k] = fValue1;
222         fSave2[k] = fValue2;
223         fSave[k] = fValue;
224     }
225 
226     foreach (k; 0 .. 5)
227     {
228         immutable x = halfLength * x2[k];
229 
230         immutable fValue1 = f(center + x);
231         immutable fValue2 = f(center - x);
232         immutable fValue = fValue1 + fValue2;
233 
234         result21 += w21b[k] * fValue;
235         resultAbs += w21b[k] * (abs(fValue1) + abs(fValue2));
236 
237         fSave3[k] = fValue1;
238         fSave4[k] = fValue2;
239         fSave[k+5] = fValue;
240     }
241 
242     // resultAbs is the absolute value of the integral of the
243     // absolute value of the function.
244     resultAbs *= absHalfLength;
245 
246     // fMean is the mean value of f(x) in the interval [a, b],
247     // as evaluated by the 21-point rule.
248     immutable fMean = result21/2;
249 
250     // variance is the integral of f(x)-fMean on the interval [a, b],
251     // evaluated using the 21-point rule.
252     T variance = w21b[5] * abs(fCenter - fMean);
253     foreach (k; 0 .. 5)
254     {
255         variance +=
256             w21a[k] * (abs(fSave1[k] - fMean) + abs(fSave2[k] - fMean))
257             + w21b[k] * (abs(fSave3[k] - fMean) + abs(fSave4[k] - fMean));
258     }
259     variance *= absHalfLength;
260 
261     // Test for convergence.
262     {
263 
264         immutable result = result21 * halfLength;
265         immutable error = errorCalc(abs((result21 - result10)*halfLength),
266             resultAbs, variance);
267 
268         if (error <= fmax(epsAbs, epsRel * abs(result)))
269             return typeof(return)(result, error);
270     }
271 
272 
273     // Compute the integral using the 43-point formula.
274     T result43 = w43b[11] * fCenter;
275     foreach (k; 0 .. 10)
276         result43 += fSave[k] * w43a[k];
277 
278     foreach (k; 0 .. 11)
279     {
280         immutable x = halfLength * x3[k];
281         immutable fValue = f(center+x) + f(center-x);
282         result43 += fValue * w43b[k];
283         fSave[k+10] = fValue;
284     }
285 
286     // Test for convergence
287     {
288         immutable result = result43 * halfLength;
289         immutable error = errorCalc(abs((result43-result21)*halfLength),
290             resultAbs, variance);
291 
292         if (error <= fmax(epsAbs, epsRel * abs(result)))
293             return typeof(return)(result, error);
294     }
295 
296 
297     // Compute the integral using the 87-point formula.
298     T result87 = w87b[22] * fCenter;
299     foreach (k; 0 .. 21)
300         result87 += fSave[k] * w87a[k];
301 
302     foreach (k; 0 .. 22)
303     {
304         immutable x = halfLength * x4[k];
305         result87 += w87b[k] * (f(center+x) + f(center-x));
306     }
307 
308     // Test for convergence
309     {
310         immutable result = result87 * halfLength;
311         immutable error = errorCalc(abs((result87-result43)*halfLength),
312             resultAbs, variance);
313 
314         if (error > fmax(epsAbs, epsRel * abs(result)))
315             throw new IntegrationException(
316                 "The maximum number of steps has been executed",
317                 "The integral is probably too difficult to be calculated "
318                 ~"by the QNG algorithm.",
319                 result, error);
320 
321         return typeof(return)(result, error);
322     }
323 }
324 
325 
326 unittest
327 {
328     // Check that it compiles for different types.
329     alias qng!(float delegate(float)) fqng;
330     alias qng!(double delegate(double)) dqng;
331     alias qng!(double function(double)) dfqng;
332     alias qng!(real delegate(real)) rqng;
333 }
334 
335 
336 unittest
337 {
338     // QUADPACK book §4.4, integral 4
339     double alpha;
340     double f(double x) { return x^^alpha * log(1/x); }
341 
342     enum epsabs = 0.0;
343     enum epsrel = 1e-8;
344 
345     foreach (i; 1 .. 9)
346     {
347         alpha = 1 + i*0.2;
348         auto result = qng(&f, 0.0, 1.0, epsabs, epsrel);
349         auto expect = (alpha + 1)^^(-2);
350         assert (isAccurate(result.value, result.error, expect, epsrel, epsabs));
351     }
352 }
353 
354 
355 unittest
356 {
357     // QUADPACK book §4.4, integral 5
358     double alpha;
359     double f(double x) { return 4^^(-alpha) / ((x-PI_4)^^2 + 16^^(-alpha)); }
360 
361     enum epsabs = 0.0;
362     enum epsrel = 1e-8;
363 
364     foreach (i; 0 .. 2)
365     {
366         alpha = i;
367         auto result = qng(&f, 0.0, 1.0, epsabs, epsrel);
368         double expect = atan((4-PI) * 4^^(alpha-1)) + atan(PI * 4^^(alpha-1));
369         assert (isAccurate(result.value, result.error, expect, epsrel, epsabs));
370     }
371 }
372 
373 
374 unittest
375 {
376     // QUADPACK book §4.4, integral 6
377     double alpha;
378     double f(double x) { return cos(2^^alpha * sin(x)); }
379 
380     enum epsabs = 0.0;
381     enum epsrel = 1e-8;
382 
383     auto expect = [
384          2.4039394306344129983,
385          0.70337362695660089178,
386         -1.2476829250428461076,
387          0.53925691468609779719,
388         -0.54946164594662718058,
389          0.43378800263473354846,
390          0.29088010217372596783
391     ];
392 
393     foreach (i; 0 .. 6)
394     {
395         alpha = i;
396         auto result = qng(&f, 0.0, cast(double) PI, epsabs, epsrel);
397         //assert (isAccurate(result.value, result.error, expect[i], epsrel, epsabs));
398         assert (abs((result.value-expect[i])/result.value) <= epsrel);
399     }
400     // TODO:
401     // Some of the error estimates are a bit too low, which is why we
402     // can't use isAccurate.  This may not be a problem, but it's worth
403     // investigating closer.
404     //
405     // What's worse is that the routine should succeed for i=6 as well,
406     // but it doesn't.  Actually, it does converge to the right value,
407     // but it erroneously reports the error as ~0.8.
408 }

Page generated by adrdox