mintegrated source code

1 module mintegrated;
2 
3 import std.array : array;
4 import std.conv : to;
5 import std.range : iota, zip;
6 import std.algorithm : all, each, map, max, reduce, filter;
7 import std.random : uniform;
8 import std.math : pow, round, sqrt;
9 
10 import scid.types : Result;
11 import dstats.summary : meanStdev, MeanSD;
12 
13 private struct Area(Real)
14 {
15     Real[] lower;
16     Real[] upper;
17 }
18 
19 private Real volume(Real)( in Area!Real area )
20 {
21     /+
22         D way: fails to compile, so use for loop
23         return zip(area.lower, area.upper)
24         .map((t) => t[1] - t[0])
25         .reduce!((a,b) => a*b );
26     +/
27     Real s = 1;
28     foreach( t; zip(area.lower, area.upper) )
29         s *= t[1] - t[0];
30     return s;
31 }
32 
33 private size_t dimension(Real)( in Area!Real area )
34 {
35     return area.lower.length;
36 }
37 
38 unittest
39 {
40     assert( volume( Area!double([-1.0,-2.0],[2.0,0.0]) ) == 3.0*2.0 );
41 }
42 
43 private Area!Real[] splitArea(Real)( in Area!Real area, size_t dimension )
44 {
45     assert( dimension < area.lower.length );
46     auto div = (area.upper[dimension]-area.lower[dimension])
47         //*0.5
48         *uniform(0.4,0.6)
49         + area.lower[dimension];
50     auto newLower = area.lower.dup;
51     newLower[dimension] = div;
52     auto newUpper = area.upper.dup;
53     newUpper[dimension] = div;
54     return [ Area!Real( area.lower.dup, newUpper ),
55            Area!Real( newLower, area.upper.dup ) ];
56 }
57 
58 unittest
59 {
60     auto a = Area!double([-1.0,-2.0],[2.0,0.0]);
61     assert( splitArea(a,0)[0].volume < a.volume );
62     assert( splitArea(a,1)[0].volume < a.volume );
63     assert( splitArea(a,0)[1].volume < a.volume );
64     assert( splitArea(a,1)[1].volume < a.volume );
65 
66     auto splitted = splitArea(a,0);
67     assert( a.volume == splitted[0].volume + splitted[1].volume );
68 }
69 
70 private bool withinArea(Real)( in Real[] point, in Area!Real area )
71 {
72     return zip(point, area.lower, area.upper).all!(
73             (t) => t[0] >= t[1] && t[0] <= t[2] );
74 }
75 
76 unittest
77 {
78     auto a = Area!double([-1.0,-2.0],[2.0,0.0]);
79     assert( [0.0,0.0].withinArea( a ) );
80     assert( ![-2.0,0.0].withinArea( a ) );
81     assert( ![3.0,0.0].withinArea( a ) );
82     assert( ![0.0,-2.1].withinArea( a ) );
83     assert( ![0.0,0.1].withinArea( a ) );
84 }
85 
86 private Result!Real meanAndVariance(Real, Range : MeanSD)( in Range msd, in Area!Real area )
87 {
88     auto v = area.volume;
89     return Result!Real( v*msd.mean().to!Real,
90             pow(v,2)*msd.mse().to!Real );
91 }
92 
93 private Result!Real meanAndVariance(Real, Range)( in Range values, in Area!Real area )
94 {
95     auto msd = meanStdev( values );
96     return msd.meanAndVariance!Real( area );
97 }
98 
99 
100 unittest
101 {
102     auto a = Area!double([-1.0,-2.0],[0.0,-1.0]);
103     auto vs = [1.0,2.0,1.5];
104     auto res = vs.meanAndVariance( a );
105     assert( res.value == 1.5 );
106     assert( res.error == 0.5/3 );
107 
108     a = Area!double([-1.0,-2.0],[1.0,-1.0]);
109     res = vs.meanAndVariance( a );
110     assert( res.value == 2*1.5 );
111     assert( res.error == 4*0.5/3 );
112 }
113 
114 Result!Real integrate(Func, Real)(scope Func f, Real[] a, Real[] b,
115     Real epsRel = cast(Real) 1e-6, Real epsAbs = cast(Real) 0)
116 {
117     auto area = Area!Real( a, b );
118     auto result = miser(f, area, epsRel, epsAbs, 150*a.length);
119     return Result!Real( result.value*area.volume, 
120             result.error*area.volume ); 
121 }
122 
123 Result!Real miser(Func, Real)(scope Func f, in Area!Real area,
124     Real epsRel = cast(Real) 1e-6, Real epsAbs = cast(Real) 0, 
125     size_t npoints = 1000, MeanSD msd = MeanSD() )
126 {
127     assert( volume(area) > 0, "Size of area is 0" );
128     auto bounds = area.lower.zip(area.upper);
129     assert( bounds.all!((t) => t[1] > t[0] ) );
130     auto points =
131         iota( 0, npoints, 1 )
132         .map!( (i) => bounds
133                 .map!( (t) {
134                     return uniform!"[]"( t[0], t[1] ).to!Real; 
135                     } 
136                 ).array 
137              );
138     auto values = points.map!((pnt) => f( pnt ) ).array;
139 
140     foreach( v; values )
141         msd.put(v);
142     auto result = Result!Real( msd.mean, sqrt(msd.mse) );
143 
144     if (area.volume*result.error < epsAbs 
145             || result.error/(area.volume*result.value) < epsRel)
146         return result;
147 
148     // Try different subareas
149     Area!Real[] bestAreas;
150     auto bestEst = Real.max;
151     Result!Real[] bestResults;
152     MeanSD[] bestMSDs;
153     foreach( j; 0..area.lower.length ) 
154     {
155         auto subAreas = area.splitArea( j );
156         assert( volume(subAreas[0]) > 0, "Cannot divide the area further" );
157 
158         auto pntvs = zip( points, values );
159         MeanSD[] cacheMSDs;
160         auto results = subAreas.map!( (a) 
161                 {
162                     MeanSD msd;
163                     foreach( pntv; zip( points, values ) )
164                         {
165                         if (pntv[0].withinArea(a))
166                         msd.put( pntv[1] );
167                         }
168                     cacheMSDs ~= msd;
169                     assert( msd.N > 0 );
170                     return Result!Real(msd.mean, sqrt(msd.mse));
171                 } );
172         Result!Real[] cacheResults;
173         // Optimize this by first only looking at first. Only if that
174         // is smaller than bestEstimate would we need to calculate second
175         Real runningError = 0;
176         while( !results.empty && runningError < bestEst )
177         {
178             runningError += results.front.error;
179             cacheResults ~= results.front;
180             results.popFront;
181         }
182 
183         if( results.empty && runningError < bestEst )
184         {
185             bestEst = runningError;
186             bestResults = cacheResults;
187             bestAreas = subAreas;
188             bestMSDs = cacheMSDs;
189         }
190     }
191     assert( bestAreas.length == 2 );
192     assert( bestResults.length == 2 );
193 
194     auto sdA = sqrt(bestResults[0].error);
195     auto sdB = sqrt(bestResults[1].error);
196     auto sumSd = sdA + sdB;
197     if (sumSd == 0)
198     {
199         result = Result!Real( 0.5*(bestResults[0].value+bestResults[1].value), 
200             sqrt(0.25*pow(bestResults[0].error,2)+0.25*pow(bestResults[1].error,2) ) );
201 
202         return result;
203     }
204 
205     auto npntsl = round(150*area.dimension*sdA/sumSd).to!int; 
206     auto npntsu = 150*area.dimension-npntsl;
207 
208     auto rl = miser( f, bestAreas[0], 
209             epsRel, epsAbs,
210             max( 15*area.dimension, npntsl ),
211             bestMSDs[0]
212             );
213     auto ru = miser( f, bestAreas[1], 
214             epsRel, epsAbs,
215             max( 15*area.dimension, npntsu ),
216             bestMSDs[1] );
217 
218     result = Result!Real( 0.5*(rl.value+ru.value), 
219             sqrt(0.25*pow(rl.error,2)+0.25*pow(ru.error,2) ) );
220 
221     return result; 
222 }
223 
224 ///
225 unittest
226 {
227     import std.math : PI, pow;
228     import std.stdio : writeln;
229     auto func = function( double[] xs )
230     {
231         if (pow(xs[0],2)+pow(xs[1],2)<= 1.0)
232             return 1.0;
233         return 0.0;
234     };
235 
236     auto result = integrate( func, [-1.0,-1], [1.0,1.0], 1e-5, 1e-7 );
237     result.writeln;
238     assert( result.value <= PI + 1e-3 );
239     assert( result.value >= PI - 1e-3 );
240 }