1 /******************************************************************************* 2 3 Copyright: 4 Copyright (c) 2008. Fawzi Mohamed 5 Some parts copyright (c) 2009-2016 dunnhumby Germany GmbH. 6 All rights reserved. 7 8 License: 9 Tango Dual License: 3-Clause BSD License / Academic Free License v3.0. 10 See LICENSE_TANGO.txt for details. 11 12 Version: Initial release: July 2008 13 14 Authors: Fawzi Mohamed 15 16 *******************************************************************************/ 17 module ocean.math.random.Ziggurat; 18 import ocean.math.Bracket:findRoot; 19 import ocean.math.Math:abs; 20 import ocean.math.ErrorFunction:erfc; 21 import ocean.meta.traits.Basic /* : isRealType */; 22 import ocean.meta.types.Arrays /* : StripAllArrays */; 23 import ocean.core.Verify; 24 25 /// ziggurat method for decreasing distributions. 26 /// Marsaglia, Tsang, Journal of Statistical Software, 2000 27 /// If has negative is true the distribution is assumed to be symmetric with respect to 0, 28 /// otherwise it is assumed to be from 0 to infinity. 29 /// Struct based to avoid extra indirection when wrapped in a class (and it should be wrapped 30 /// in a class and not used directly). 31 /// Call style initialization avoided on purpose (this is a big structure, you don't want to return it) 32 struct Ziggurat(RandG,T,alias probDensityF,alias tailGenerator,bool hasNegative=true){ 33 static assert(isRealType!(T),T.stringof~" not acceptable, only floating point variables supported"); 34 enum int nBlocks=256; 35 T[nBlocks+1] posBlock; 36 T[nBlocks+1] fVal; 37 RandG r; 38 alias Ziggurat SourceT; 39 /// initializes the ziggurat 40 static Ziggurat create(alias invProbDensityF, alias cumProbDensityFCompl)(RandG rGenerator,real xLast=-1.0L,bool check_error=true){ 41 /// function to find xLast 42 real findXLast(real xLast){ 43 real v=xLast*probDensityF(xLast)+cumProbDensityFCompl(xLast); 44 real fMax=probDensityF(0.0L); 45 real pAtt=xLast; 46 real fAtt=probDensityF(xLast); 47 for (int i=nBlocks-2;i!=0;--i){ 48 fAtt+=v/pAtt; 49 if (fAtt>fMax) return fAtt+(i-1)*fMax; 50 pAtt=invProbDensityF(fAtt); 51 verify(pAtt>=0,"invProbDensityF is supposed to return positive values"); 52 } 53 return fAtt+v/pAtt-fMax; 54 } 55 void findBracket(ref real xMin,ref real xMax){ 56 real vMin=cumProbDensityFCompl(0.0L)/nBlocks; 57 real pAtt=0.0L; 58 for (int i=1;i<nBlocks;++i){ 59 pAtt+=vMin/probDensityF(pAtt); 60 } 61 real df=findXLast(pAtt); 62 if (df>0) { 63 // (most likely) 64 xMin=pAtt; 65 real vMax=cumProbDensityFCompl(0.0L); 66 xMax=pAtt+vMax/probDensityF(pAtt); 67 } else { 68 xMax=pAtt; 69 xMin=vMin/probDensityF(0.0L); 70 } 71 } 72 if (xLast<=0){ 73 real xMin,xMax; 74 findBracket(xMin,xMax); 75 xLast=findRoot(&findXLast,xMin,xMax); 76 // printf("xLast:%La => %La\n",xLast,findXLast(xLast)); 77 } 78 Ziggurat res; 79 with (res){ 80 r=rGenerator; 81 real v=probDensityF(xLast)*xLast+cumProbDensityFCompl(xLast); 82 real pAtt=xLast; 83 real fMax=probDensityF(0.0L); 84 posBlock[1]=cast(T)xLast; 85 real fAtt=probDensityF(xLast); 86 fVal[1]=cast(T)fAtt; 87 for (int i=2;i<nBlocks;++i){ 88 fAtt+=v/pAtt; 89 verify(fAtt<=fMax,"Ziggurat contruction shoot out"); 90 pAtt=invProbDensityF(fAtt); 91 verify(pAtt>=0,"invProbDensityF is supposed to return positive values"); 92 posBlock[i]=cast(T)pAtt; 93 fVal[i]=cast(T)fAtt; 94 } 95 posBlock[nBlocks]=0.0L; 96 fVal[nBlocks]=cast(T)probDensityF(0.0L); 97 real error=fAtt+v/pAtt-probDensityF(0.0L); 98 verify((!check_error) || error<real.epsilon*10000.0,"Ziggurat error larger than expected"); 99 posBlock[0]=cast(T)(xLast*(1.0L+cumProbDensityFCompl(xLast)/probDensityF(xLast))); 100 fVal[0]=0.0L; 101 for (int i=0;i<nBlocks;++i){ 102 verify(posBlock[i]>=posBlock[i+1],"decresing posBlock"); 103 verify(fVal[i]<=fVal[i+1],"increasing probabilty density function"); 104 } 105 } 106 return res; 107 } 108 /// returns a single value with the probability distribution of the current Ziggurat 109 /// and slightly worse randomness (in the normal case uses only 32 random bits). 110 /// Cannot be 0. 111 T getFastRandom() 112 { 113 static if (hasNegative){ 114 for (int iter=1000;iter!=0;--iter) 115 { 116 uint i0=r.uniform!(uint)(); 117 uint i=i0 & 0xFFU; 118 enum T scaleF=(cast(T)1)/(cast(T)uint.max+1); 119 T u= (cast(T)i0+cast(T)0.5)*scaleF; 120 T x = posBlock[i]*u; 121 if (x<posBlock[i+1]) return ((i0 & 0x100u)?x:-x); 122 if (i == 0) return tailGenerator(r,posBlock[1],x<0); 123 if ((cast(T)probDensityF(x))>fVal[i+1]+(fVal[i]-fVal[i+1])*((cast(T)r.uniform!(uint)+cast(T)0.5)*scaleF)) { 124 return ((i0 & 0x100u)?x:-x); 125 } 126 } 127 } else { 128 for (int iter=1000;iter!=0;--iter) 129 { 130 uint i0=r.uniform!(uint); 131 uint i= i0 & 0xFFU; 132 enum T scaleF=(cast(T)1)/(cast(T)uint.max+1); 133 T u= (cast(T)i0+cast(T)0.5)*scaleF; 134 T x = posBlock[i]*u; 135 if (x<posBlock[i+1]) return x; 136 if (i == 0) return tailGenerator(r,posBlock[1]); 137 if ((cast(T)probDensityF(x))>fVal[i+1]+(fVal[i]-fVal[i+1])*r.uniform!(T)) { 138 return x; 139 } 140 } 141 } 142 throw new Exception("max nr of iterations in Ziggurat, this should have probability<1.0e-1000"); 143 } 144 /// returns a single value with the probability distribution of the current Ziggurat 145 T getRandom() 146 { 147 static if (hasNegative){ 148 for (int iter=1000;iter!=0;--iter) 149 { 150 uint i0 = r.uniform!(uint); 151 uint i= i0 & 0xFF; 152 T u = r.uniform!(T)(); 153 T x = posBlock[i]*u; 154 if (x<posBlock[i+1]) return ((i0 & 0x100u)?x:-x); 155 if (i == 0) return tailGenerator(r,posBlock[1],x<0); 156 if ((cast(T)probDensityF(x))>fVal[i+1]+(fVal[i]-fVal[i+1])*r.uniform!(T)) { 157 return ((i0 & 0x100u)?x:-x); 158 } 159 } 160 } else { 161 for (int iter=1000;iter!=0;--iter) 162 { 163 uint i=r.uniform!(ubyte); 164 T u = r.uniform!(T)(); 165 T x = posBlock[i]*u; 166 if (x<posBlock[i+1]) return x; 167 if (i == 0) return tailGenerator(r,posBlock[1]); 168 if ((cast(T)probDensityF(x))>fVal[i+1]+(fVal[i]-fVal[i+1])*r.uniform!(T)) { 169 return x; 170 } 171 } 172 } 173 throw new Exception("max nr of iterations in Ziggurat, this should have probability<1.0e-1000"); 174 } 175 /// initializes the argument with the probability distribution given and returns it 176 /// for arrays this might potentially be faster than a naive loop 177 U randomize(U)(ref U a){ 178 static if(is(U S:S[])){ 179 size_t aL = a.length; 180 for (size_t i=0;i!=aL;++i){ 181 a[i]=cast(StripAllArrays!(U))getRandom(); 182 } 183 } else { 184 a=cast(U)getRandom(); 185 } 186 return a; 187 } 188 /// initializes the variable with the result of mapping op on the random numbers (of type T) 189 // unfortunately this (more efficent version) cannot use local delegates 190 template randomizeOp2(alias op){ 191 U randomizeOp2(U)(ref U a){ 192 static if(is(U S:S[])){ 193 alias StripAllArrays!(U) TT; 194 uint aL=a.length; 195 for (uint i=0;i!=aL;++i){ 196 static if(isComplexType!(TT)) { 197 a[i]=cast(TT)(op(getRandom())+1i*op(getRandom())); 198 } else static if (isImaginaryType!(TT)){ 199 a[i]=cast(TT)(1i*op(getRandom())); 200 } else { 201 a[i]=cast(TT)op(getRandom()); 202 } 203 } 204 } else { 205 static if(isComplexType!(U)) { 206 a=cast(U)(op(getRandom())+1i*op(getRandom())); 207 } else static if (isImaginaryType!(U)){ 208 el=cast(U)(1i*op(getRandom())); 209 } else { 210 a=cast(U)op(getRandom()); 211 } 212 } 213 return a; 214 } 215 } 216 /// initializes the variable with the result of mapping op on the random numbers (of type T) 217 U randomizeOp(U,S)(scope S delegate(T) op,ref U a){ 218 static if(is(U S:S[])){ 219 alias StripAllArrays!(U) TT; 220 size_t aL = a.length; 221 for (size_t i=0;i!=aL;++i){ 222 static if(isComplexType!(TT)) { 223 a[i]=cast(TT)(op(getRandom())+1i*op(getRandom())); 224 } else static if (isImaginaryType!(TT)){ 225 a[i]=cast(TT)(1i*op(getRandom())); 226 } else { 227 a[i]=cast(TT)op(getRandom()); 228 } 229 } 230 } else { 231 static if(isComplexType!(U)) { 232 a=cast(U)(op(getRandom())+1i*op(getRandom())); 233 } else static if (isImaginaryType!(U)){ 234 el=cast(U)(1i*op(getRandom())); 235 } else { 236 a=cast(U)op(getRandom()); 237 } 238 } 239 return a; 240 } 241 242 }