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, 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, 100000*a.length); 120 return Result!Real( result.value, 121 result.error ); 122 } 123 124 /// 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 ); 136 auto points = 137 iota( 0, dim, 1 ) 138 .map!( (i) => bounds 139 .map!( (t) { 140 return uniform!"[]"( t[0], t[1] ).to!Real; 141 } 142 ).array 143 ); 144 auto values = points.map!((pnt) => f( pnt ) ).array; 145 146 auto result = values.meanAndVariance(area); 147 148 if ( npoints < minPoints 149 || result.error < epsAbs 150 || result.error/result.value < epsRel ) 151 return result; 152 153 // Try different subareas 154 Area!Real[] bestAreas; 155 auto bestEst = Real.max; 156 Result!Real[] bestResults; 157 foreach( j; 0..area.lower.length ) 158 { 159 auto subAreas = area.splitArea( j ); 160 assert( volume(subAreas[0]) > 0, "Cannot divide the area further" ); 161 162 auto pntvs = zip( points, values ); 163 MeanSD[] msds = new MeanSD[2]; 164 foreach( pntv; zip( points, values ) ) 165 { 166 if (pntv[0].withinArea(subAreas[0])) 167 msds[0].put( pntv[1] ); 168 if (pntv[0].withinArea(subAreas[1])) 169 msds[1].put( pntv[1] ); 170 } 171 172 auto results = msds.zip(subAreas).map!((msd) => 173 meanAndVariance(msd[0], msd[1]) ); 174 Result!Real[] cacheResults; 175 // Optimize this by first only looking at first. Only if that 176 // is smaller than bestEstimate would we need to calculate second 177 Real runningError = 0; 178 while( !results.empty && runningError < bestEst ) 179 { 180 runningError += results.front.error; 181 cacheResults ~= results.front; 182 results.popFront; 183 } 184 185 if( results.empty && runningError < bestEst ) 186 { 187 bestEst = runningError; 188 bestResults = cacheResults; 189 bestAreas = subAreas; 190 } 191 } 192 assert( bestAreas.length == 2 ); 193 assert( bestResults.length == 2 ); 194 195 auto sdA = bestResults[0].error; 196 auto sdB = bestResults[1].error; 197 auto sumSd = sdA + sdB; 198 //assert( sumSd > 0, "Sum Errors to small" ); 199 if (sumSd == 0) 200 { 201 result = Result!Real( bestResults[0].value+bestResults[1].value, 202 sqrt( bestResults[0].error + bestResults[1].error ) ); 203 204 return result; 205 } 206 207 auto leftOverPoints = round( 0.9*npoints ).to!int; 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 sqrt(pow(rl.error,2)+pow(ru.error,2) ) ); 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 + 1e-2 ); 236 assert( result.value >= PI - 1e-2 ); 237 }