ocean.math.random.ExpSource source code

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.ExpSource;
18 import Integer = ocean.text.convert.Integer_tango;
19 import ocean.math.Math:exp,log;
20 import ocean.math.random.Ziggurat;
21 import ocean.meta.traits.Basic /* : isRealType */;
22 
23 /// class that returns exponential distributed numbers (f=exp(-x) for x>0, 0 otherwise)
24 final class ExpSource(RandG,T){
25     static assert(isRealType!(T),T.stringof~" not acceptable, only floating point variables supported");
26     /// probability distribution
27     static real probDensityF(real x){ return exp(-x); }
28     /// inverse probability distribution
29     static real invProbDensityF(real x){ return -log(x); }
30     /// complement of the cumulative density distribution (integral x..infinity probDensityF)
31     static real cumProbDensityFCompl(real x){ return exp(-x); }
32     /// tail for exponential distribution
33     static T tailGenerator(RandG r, T dMin)
34     {
35         return dMin-log(r.uniform!(T));
36     }
37     alias Ziggurat!(RandG,T,probDensityF,tailGenerator,false) SourceT;
38     /// internal source of exp distribued numbers
39     SourceT source;
40     /// initializes the probability distribution
41     this(RandG r){
42         source=SourceT.create!(invProbDensityF,cumProbDensityFCompl)(r,0xf.64ec94bf5dc14bcp-1L);
43     }
44     /// chainable call style initialization of variables (thorugh a call to randomize)
45     ExpSource opCall(U,S...)(ref U a,S args){
46         randomize(a,args);
47         return this;
48     }
49     /// returns a exp distribued number
50     T getRandom(){
51         return source.getRandom();
52     }
53     /// returns a exp distribued number with the given beta (survival rate, average)
54     /// f=1/beta*exp(-x/beta)
55     T getRandom(T beta){
56         return beta*source.getRandom();
57     }
58     /// initializes the given variable with an exponentially distribued number
59     U randomize(U)(ref U x){
60         return source.randomize(x);
61     }
62     /// initializes the given variable with an exponentially distribued number with
63     /// scale parameter beta
64     U randomize(U,V)(ref U x,V beta){
65         return source.randomizeOp((T el){ return el*cast(T)beta; },x);
66     }
67     /// initializes the given variable with an exponentially distribued number and maps op on it
68     U randomizeOp(U,S)(scope S delegate(T)op,ref U a){
69         return source.randomizeOp(op,a);
70     }
71     /// exp distribution with different default scale parameter beta
72     /// f=1/beta*exp(-x/beta) for x>0, 0 otherwise
73     struct ExpDistribution{
74         T beta;
75         ExpSource source; // does not use Ziggurat directly to keep this struct small
76         /// constructor
77         static ExpDistribution create()(ExpSource source,T beta){
78             ExpDistribution res;
79             res.beta=beta;
80             res.source=source;
81             return res;
82         }
83         /// chainable call style initialization of variables (thorugh a call to randomize)
84         ExpDistribution opCall(U,S...)(ref U a,S args){
85             randomize(a,args);
86             return this;
87         }
88         /// returns a single number
89         T getRandom(){
90             return beta*source.getRandom();
91         }
92         /// initialize a
93         U randomize(U)(ref U a){
94             return source.randomizeOp((T x){return beta*x; },a);
95         }
96         /// initialize a
97         U randomize(U,V)(ref U a,V b){
98             return source.randomizeOp((T x){return (cast(T)b)*x; },a);
99         }
100     }
101     /// returns an exp distribution with a different beta
102     ExpDistribution expD(T beta){
103         return ExpDistribution.create(this,beta);
104     }
105 }