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 }