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 /// 115 Result!Real integrate(Func, Real)(scope Func f, Real[] a, Real[] b, 116 Real epsRel = cast(Real) 1e-6, Real epsAbs = cast(Real) 0) 117 { 118 auto area = Area!Real( a, b ); 119 auto result = miser(f, area, epsRel, epsAbs, 200000*a.length); 120 return Result!Real( result.value, 121 result.error ); 122 } 123 124 /// The returned error is the expected variance in the result 125 Result!Real miser(Func, Real)(scope Func f, in Area!Real area, 126 Real epsRel = cast(Real) 1e-6, Real epsAbs = cast(Real) 0, 127 size_t npoints = 1000 ) 128 { 129 assert( volume(area) > 0, "Size of area is 0" ); 130 auto bounds = area.lower.zip(area.upper); 131 assert( bounds.all!((t) => t[1] > t[0] ) ); 132 133 auto minPoints = 15*area.dimension; 134 135 auto dim = max( 0.1*npoints, minPoints ).to!int; 136 auto leftOverPoints = npoints - dim; 137 auto points = 138 iota( 0, dim, 1 ) 139 .map!( (i) => bounds 140 .map!( (t) { 141 return uniform!"[]"( t[0], t[1] ).to!Real; 142 } 143 ).array 144 ); 145 auto values = points.map!((pnt) => f( pnt ) ).array; 146 147 auto result = values.meanAndVariance(area); 148 149 if ( npoints < minPoints 150 //|| result.error < epsAbs 151 || result.error/result.value < epsRel ) 152 return result; 153 154 // Try different subareas 155 Area!Real[] bestAreas; 156 auto bestEst = Real.max; 157 Result!Real[] bestResults; 158 foreach( j; 0..area.dimension ) 159 { 160 auto subAreas = area.splitArea( j ); 161 assert( volume(subAreas[0]) > 0, "Cannot divide the area further" ); 162 163 auto pntvs = zip( points, values ); 164 MeanSD[] msds = new MeanSD[2]; 165 foreach( pntv; zip( points, values ) ) 166 { 167 if (pntv[0].withinArea(subAreas[0])) 168 msds[0].put( pntv[1] ); 169 if (pntv[0].withinArea(subAreas[1])) 170 msds[1].put( pntv[1] ); 171 } 172 173 auto results = msds.zip(subAreas).map!((msd) => 174 meanAndVariance(msd[0], msd[1]) ); 175 Result!Real[] cacheResults; 176 // Optimize this by first only looking at first. Only if that 177 // is smaller than bestEstimate would we need to calculate second 178 Real runningError = 0; 179 while( !results.empty && runningError < bestEst ) 180 { 181 runningError += results.front.error; 182 cacheResults ~= results.front; 183 results.popFront; 184 } 185 186 if( results.empty && runningError < bestEst ) 187 { 188 bestEst = runningError; 189 bestResults = cacheResults; 190 bestAreas = subAreas; 191 } 192 } 193 assert( bestAreas.length == 2 ); 194 assert( bestResults.length == 2 ); 195 196 auto sdA = sqrt(bestResults[0].error); 197 auto sdB = sqrt(bestResults[1].error); 198 auto sumSd = sdA + sdB; 199 //assert( sumSd > 0, "Sum Errors to small" ); 200 if (sumSd == 0) 201 { 202 result = Result!Real( bestResults[0].value+bestResults[1].value, 203 sqrt( bestResults[0].error + bestResults[1].error ) ); 204 205 return result; 206 } 207 208 auto npntsl = round(leftOverPoints*sdA/sumSd).to!int; 209 210 auto rl = miser( f, bestAreas[0], 211 epsRel, epsAbs, npntsl ); 212 auto ru = miser( f, bestAreas[1], 213 epsRel, epsAbs, leftOverPoints-npntsl ); 214 215 result = Result!Real( rl.value+ru.value, 216 rl.error+ru.error ); 217 218 return result; 219 } 220 221 /// 222 unittest 223 { 224 import std.math : PI, pow; 225 import std.stdio : writeln; 226 auto func = function( double[] xs ) 227 { 228 if (pow(xs[0],2)+pow(xs[1],2)<= 1.0) 229 return 1.0; 230 return 0.0; 231 }; 232 233 auto result = integrate( func, [-1.0,-1], [1.0,1.0], 1e-5, 0 ); 234 result.writeln; 235 assert( result.value <= PI + 3*sqrt(result.error) ); 236 assert( result.value >= PI - 3*sqrt(result.error) ); 237 } 238 239 /// 240 unittest 241 { 242 import std.stdio : writeln; 243 auto func = function(double[] xs ) 244 { 245 return xs[0]*xs[1]; 246 }; 247 auto result = integrate( func, [0.0,0], [1.0,1] ); 248 result.writeln; 249 assert( result.value <= 0.25 + 3*sqrt(result.error) ); 250 assert( result.value >= 0.25 - 3*sqrt(result.error) ); 251 } 252 253 /// 254 unittest 255 { 256 import std.math : PI, cos; 257 import std.stdio : writeln; 258 auto func = function(real[] xs ) 259 { 260 return 1.0/(pow(PI,3)*(1-cos(xs[0])*cos(xs[1])*cos(xs[2]))); 261 }; 262 auto result = integrate( func, [0,0,0], [PI,PI,PI] ); 263 result.writeln; 264 assert( result.value <= 1.393204 + 3*sqrt(result.error) ); 265 assert( result.value >= 1.393204 - 3*sqrt(result.error) ); 266 }