tango.math.random.ExpSource source code

1 /*******************************************************************************
2         copyright:      Copyright (c) 2008. Fawzi Mohamed
3         license:        BSD style: $(LICENSE)
4         version:        Initial release: July 2008
5         author:         Fawzi Mohamed
6 *******************************************************************************/
7 module tango.math.random.ExpSource;
8 private import Integer = tango.text.convert.Integer;
9 import tango.math.Math:exp,log;
10 import tango.math.random.Ziggurat;
11 import tango.core.Traits:isRealType;
12 
13 /// class that returns exponential distributed numbers (f=exp(-x) for x>0, 0 otherwise)
14 final class ExpSource(RandG,T){
15     static assert(isRealType!(T),T.stringof~" not acceptable, only floating point variables supported");
16     /// probability distribution
17     static real probDensityF(real x){ return exp(-x); }
18     /// inverse probability distribution
19     static real invProbDensityF(real x){ return -log(x); }
20     /// complement of the cumulative density distribution (integral x..infinity probDensityF)
21     static real cumProbDensityFCompl(real x){ return exp(-x); }
22     /// tail for exponential distribution
23     static T tailGenerator(RandG r, T dMin) 
24     { 
25         return dMin-log(r.uniform!(T)());
26     }
27     alias Ziggurat!(RandG,T,probDensityF,tailGenerator,false) SourceT;
28     /// internal source of exp distribued numbers
29     SourceT source;
30     /// initializes the probability distribution
31     this(RandG r){
32         source=SourceT.create!(invProbDensityF,cumProbDensityFCompl)(r,0xf.64ec94bf5dc14bcp-1L);
33     }
34     /// chainable call style initialization of variables (thorugh a call to randomize)
35     ExpSource opCall(U,S...)(ref U a,S args){
36         randomize(a,args);
37         return this;
38     }
39     /// returns a exp distribued number
40     T getRandom(){
41         return source.getRandom();
42     }
43     /// returns a exp distribued number with the given beta (survival rate, average)
44     /// f=1/beta*exp(-x/beta)
45     T getRandom(T beta){
46         return beta*source.getRandom();
47     }
48     /// initializes the given variable with an exponentially distribued number
49     U randomize(U)(ref U x){
50         return source.randomize(x);
51     }
52     /// initializes the given variable with an exponentially distribued number with
53     /// scale parameter beta
54     U randomize(U,V)(ref U x,V beta){
55         return source.randomizeOp((T el){ return el*cast(T)beta; },x);
56     }
57     /// initializes the given variable with an exponentially distribued number and maps op on it
58     U randomizeOp(U,S)(scope S delegate(T)op,ref U a){
59         return source.randomizeOp(op,a);
60     }
61     /// exp distribution with different default scale parameter beta
62     /// f=1/beta*exp(-x/beta) for x>0, 0 otherwise
63     struct ExpDistribution{
64         T beta;
65         ExpSource source; // does not use Ziggurat directly to keep this struct small
66         /// constructor
67         static ExpDistribution create()(ExpSource source,T beta){
68             ExpDistribution res;
69             res.beta=beta;
70             res.source=source;
71             return res;
72         }
73         /// chainable call style initialization of variables (thorugh a call to randomize)
74         ExpDistribution opCall(U,S...)(ref U a,S args){
75             randomize(a,args);
76             return this;
77         }
78         /// returns a single number
79         T getRandom(){
80             return beta*source.getRandom();
81         }
82         /// initialize a
83         U randomize(U)(ref U a){
84             return source.randomizeOp((T x){return beta*x; },a);
85         }
86         /// initialize a
87         U randomize(U,V)(ref U a,V b){
88             return source.randomizeOp((T x){return (cast(T)b)*x; },a);
89         }
90     }
91     /// returns an exp distribution with a different beta
92     ExpDistribution expD(T beta){
93         return ExpDistribution.create(this,beta);
94     }
95 }