Free cookie consent management tool by TermsFeed Policy Generator

source: stable/HeuristicLab.Random/3.3/NormalDistributedRandom.cs @ 17330

Last change on this file since 17330 was 17181, checked in by swagner, 5 years ago

#2875: Merged r17180 from trunk to stable

File size: 15.3 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 *
5 * This file is part of HeuristicLab.
6 *
7 * HeuristicLab is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * HeuristicLab is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
19 */
20#endregion
21
22using System;
23using HeuristicLab.Common;
24using HeuristicLab.Core;
25using HEAL.Attic;
26
27namespace HeuristicLab.Random {
28
29  /// <summary>
30  /// Normally distributed random variable.
31  /// Uses the Ziggurat transform to efficiently generate normal distributed random numbers.
32  /// See "The Ziggurat Method for Generating Random Variables" (G. Marsaglia and W.W. Tsang 2000).
33  /// </summary>
34  [Item("NormalDistributedRandom", "A pseudo random number generator which uses the Ziggurat method to create normally distributed random numbers.")]
35  [StorableType("240F7050-C05B-4E85-82C3-2871FE2138B1")]
36  public sealed class NormalDistributedRandom : Item, IRandom {
37    [Storable]
38    private double mu;
39    /// <summary>
40    /// Gets or sets the value for µ.
41    /// </summary>
42    public double Mu {
43      get { return mu; }
44      set { mu = value; }
45    }
46
47    [Storable]
48    private double sigma;
49    /// <summary>
50    /// Gets or sets the value for sigma.
51    /// </summary>
52    public double Sigma {
53      get { return sigma; }
54      set { sigma = value; }
55    }
56
57    [Storable]
58    private IRandom uniform;
59
60    private static double[] w = new double[] {
61      1.7290404664e-09,
62      1.2680928529e-10,
63      1.6897518107e-10,
64      1.9862687883e-10,
65      2.2232431174e-10,
66      2.4244936614e-10,
67      2.6016130916e-10,
68      2.7611987696e-10,
69      2.9073962682e-10,
70      3.0429969655e-10,
71      3.1699795566e-10,
72      3.2898020419e-10,
73      3.4035738117e-10,
74      3.5121602848e-10,
75      3.6162509098e-10,
76      3.7164057942e-10,
77      3.8130856805e-10,
78      3.9066758162e-10,
79      3.9975012189e-10,
80      4.0858399997e-10,
81      4.1719308563e-10,
82      4.2559822333e-10,
83      4.3381759296e-10,
84      4.4186720949e-10,
85      4.4976131153e-10,
86      4.5751258337e-10,
87      4.6513240481e-10,
88      4.7263104541e-10,
89      4.8001774777e-10,
90      4.8730097735e-10,
91      4.9448850570e-10,
92      5.0158732723e-10,
93      5.0860404777e-10,
94      5.1554460700e-10,
95      5.2241466708e-10,
96      5.2921933502e-10,
97      5.3596349581e-10,
98      5.4265170135e-10,
99      5.4928817050e-10,
100      5.5587695558e-10,
101      5.6242188684e-10,
102      5.6892646150e-10,
103      5.7539412124e-10,
104      5.8182819673e-10,
105      5.8823168558e-10,
106      5.9460769641e-10,
107      6.0095900478e-10,
108      6.0728838625e-10,
109      6.1359850534e-10,
110      6.1989202660e-10,
111      6.2617133700e-10,
112      6.3243904558e-10,
113      6.3869737277e-10,
114      6.4494881657e-10,
115      6.5119559745e-10,
116      6.5744004685e-10,
117      6.6368432972e-10,
118      6.6993072201e-10,
119      6.7618144417e-10,
120      6.8243871665e-10,
121      6.8870464887e-10,
122      6.9498151678e-10,
123      7.0127148533e-10,
124      7.0757677495e-10,
125      7.1389966161e-10,
126      7.2024242126e-10,
127      7.2660727435e-10,
128      7.3299660786e-10,
129      7.3941280876e-10,
130      7.4585826404e-10,
131      7.5233547170e-10,
132      7.5884698525e-10,
133      7.6539541372e-10,
134      7.7198347714e-10,
135      7.7861395109e-10,
136      7.8528972214e-10,
137      7.9201378789e-10,
138      7.9878920145e-10,
139      8.0561923799e-10,
140      8.1250728368e-10,
141      8.1945689123e-10,
142      8.2647166888e-10,
143      8.3355555791e-10,
144      8.4071272166e-10,
145      8.4794732347e-10,
146      8.5526402627e-10,
147      8.6266754851e-10,
148      8.7016316375e-10,
149      8.7775620106e-10,
150      8.8545243360e-10,
151      8.9325818964e-10,
152      9.0117996399e-10,
153      9.0922497309e-10,
154      9.1740082198e-10,
155      9.2571583732e-10,
156      9.3417884539e-10,
157      9.4279972718e-10,
158      9.5158891877e-10,
159      9.6055785548e-10,
160      9.6971930486e-10,
161      9.7908692265e-10,
162      9.8867602993e-10,
163      9.9850361313e-10,
164      1.0085882129e-09,
165      1.0189509236e-09,
166      1.0296150599e-09,
167      1.0406069340e-09,
168      1.0519566329e-09,
169      1.0636980186e-09,
170      1.0758701707e-09,
171      1.0885182755e-09,
172      1.1016947354e-09,
173      1.1154610569e-09,
174      1.1298901814e-09,
175      1.1450695947e-09,
176      1.1611052120e-09,
177      1.1781275955e-09,
178      1.1962995039e-09,
179      1.2158286600e-09,
180      1.2369856250e-09,
181      1.2601323318e-09,
182      1.2857697129e-09,
183      1.3146201905e-09,
184      1.3477839955e-09,
185      1.3870635751e-09,
186      1.4357403044e-09,
187      1.5008658760e-09,
188      1.6030947680e-09
189    };
190    private static long[] k = new long[] {
191      1991057938,
192      0,
193      1611602771,
194      1826899878,
195      1918584482,
196      1969227037,
197      2001281515,
198      2023368125,
199      2039498179,
200      2051788381,
201      2061460127,
202      2069267110,
203      2075699398,
204      2081089314,
205      2085670119,
206      2089610331,
207      2093034710,
208      2096037586,
209      2098691595,
210      2101053571,
211      2103168620,
212      2105072996,
213      2106796166,
214      2108362327,
215      2109791536,
216      2111100552,
217      2112303493,
218      2113412330,
219      2114437283,
220      2115387130,
221      2116269447,
222      2117090813,
223      2117856962,
224      2118572919,
225      2119243101,
226      2119871411,
227      2120461303,
228      2121015852,
229      2121537798,
230      2122029592,
231      2122493434,
232      2122931299,
233      2123344971,
234      2123736059,
235      2124106020,
236      2124456175,
237      2124787725,
238      2125101763,
239      2125399283,
240      2125681194,
241      2125948325,
242      2126201433,
243      2126441213,
244      2126668298,
245      2126883268,
246      2127086657,
247      2127278949,
248      2127460589,
249      2127631985,
250      2127793506,
251      2127945490,
252      2128088244,
253      2128222044,
254      2128347141,
255      2128463758,
256      2128572095,
257      2128672327,
258      2128764606,
259      2128849065,
260      2128925811,
261      2128994934,
262      2129056501,
263      2129110560,
264      2129157136,
265      2129196237,
266      2129227847,
267      2129251929,
268      2129268426,
269      2129277255,
270      2129278312,
271      2129271467,
272      2129256561,
273      2129233410,
274      2129201800,
275      2129161480,
276      2129112170,
277      2129053545,
278      2128985244,
279      2128906855,
280      2128817916,
281      2128717911,
282      2128606255,
283      2128482298,
284      2128345305,
285      2128194452,
286      2128028813,
287      2127847342,
288      2127648860,
289      2127432031,
290      2127195339,
291      2126937058,
292      2126655214,
293      2126347546,
294      2126011445,
295      2125643893,
296      2125241376,
297      2124799783,
298      2124314271,
299      2123779094,
300      2123187386,
301      2122530867,
302      2121799464,
303      2120980787,
304      2120059418,
305      2119015917,
306      2117825402,
307      2116455471,
308      2114863093,
309      2112989789,
310      2110753906,
311      2108037662,
312      2104664315,
313      2100355223,
314      2094642347,
315      2086670106,
316      2074676188,
317      2054300022,
318      2010539237
319   };
320    private static double[] f = new double[] {
321      1.0000000000e+00,
322      9.6359968185e-01,
323      9.3628269434e-01,
324      9.1304361820e-01,
325      8.9228165150e-01,
326      8.7324303389e-01,
327      8.5550057888e-01,
328      8.3878362179e-01,
329      8.2290720940e-01,
330      8.0773830414e-01,
331      7.9317700863e-01,
332      7.7914607525e-01,
333      7.6558417082e-01,
334      7.5244158506e-01,
335      7.3967725039e-01,
336      7.2725689411e-01,
337      7.1515148878e-01,
338      7.0333611965e-01,
339      6.9178915024e-01,
340      6.8049186468e-01,
341      6.6942769289e-01,
342      6.5858197212e-01,
343      6.4794182777e-01,
344      6.3749545813e-01,
345      6.2723249197e-01,
346      6.1714339256e-01,
347      6.0721951723e-01,
348      5.9745317698e-01,
349      5.8783704042e-01,
350      5.7836467028e-01,
351      5.6902998686e-01,
352      5.5982738733e-01,
353      5.5075180531e-01,
354      5.4179835320e-01,
355      5.3296267986e-01,
356      5.2424055338e-01,
357      5.1562821865e-01,
358      5.0712203979e-01,
359      4.9871864915e-01,
360      4.9041482806e-01,
361      4.8220765591e-01,
362      4.7409430146e-01,
363      4.6607214212e-01,
364      4.5813870430e-01,
365      4.5029163361e-01,
366      4.4252872467e-01,
367      4.3484783173e-01,
368      4.2724698782e-01,
369      4.1972434521e-01,
370      4.1227802634e-01,
371      4.0490642190e-01,
372      3.9760786295e-01,
373      3.9038079977e-01,
374      3.8322380185e-01,
375      3.7613546848e-01,
376      3.6911445856e-01,
377      3.6215949059e-01,
378      3.5526937246e-01,
379      3.4844297171e-01,
380      3.4167915583e-01,
381      3.3497685194e-01,
382      3.2833510637e-01,
383      3.2175290585e-01,
384      3.1522938609e-01,
385      3.0876362324e-01,
386      3.0235484242e-01,
387      2.9600214958e-01,
388      2.8970485926e-01,
389      2.8346219659e-01,
390      2.7727350593e-01,
391      2.7113807201e-01,
392      2.6505529881e-01,
393      2.5902456045e-01,
394      2.5304529071e-01,
395      2.4711695313e-01,
396      2.4123899639e-01,
397      2.3541094363e-01,
398      2.2963231802e-01,
399      2.2390270233e-01,
400      2.1822164953e-01,
401      2.1258877218e-01,
402      2.0700371265e-01,
403      2.0146611333e-01,
404      1.9597564638e-01,
405      1.9053204358e-01,
406      1.8513499200e-01,
407      1.7978426814e-01,
408      1.7447963357e-01,
409      1.6922089458e-01,
410      1.6400785744e-01,
411      1.5884037316e-01,
412      1.5371830761e-01,
413      1.4864157140e-01,
414      1.4361007512e-01,
415      1.3862377405e-01,
416      1.3368265331e-01,
417      1.2878671288e-01,
418      1.2393598258e-01,
419      1.1913054436e-01,
420      1.1437050998e-01,
421      1.0965602100e-01,
422      1.0498725623e-01,
423      1.0036443919e-01,
424      9.5787845552e-02,
425      9.1257803142e-02,
426      8.6774669588e-02,
427      8.2338899374e-02,
428      7.7950984240e-02,
429      7.3611505330e-02,
430      6.9321118295e-02,
431      6.5080583096e-02,
432      6.0890771449e-02,
433      5.6752663106e-02,
434      5.2667401731e-02,
435      4.8636294901e-02,
436      4.4660862535e-02,
437      4.0742866695e-02,
438      3.6884389818e-02,
439      3.3087886870e-02,
440      2.9356317595e-02,
441      2.5693291798e-02,
442      2.2103304043e-02,
443      1.8592102453e-02,
444      1.5167297795e-02,
445      1.1839478277e-02,
446      8.6244847625e-03,
447      5.5489949882e-03,
448      2.6696291752e-03
449    };
450
451    /// <summary>
452    /// Used by HeuristicLab.Persistence to initialize new instances during deserialization.
453    /// </summary>
454    /// <param name="deserializing">true, if the constructor is called during deserialization.</param>
455    [StorableConstructor]
456    private NormalDistributedRandom(StorableConstructorFlag _) : base(_) { }
457
458    /// <summary>
459    /// Initializes a new instance from an existing one (copy constructor).
460    /// </summary>
461    /// <param name="original">The original <see cref="NormalDistributedRandom"/> instance which is used to initialize the new instance.</param>
462    /// <param name="cloner">A <see cref="Cloner"/> which is used to track all already cloned objects in order to avoid cycles.</param>
463    private NormalDistributedRandom(NormalDistributedRandom original, Cloner cloner)
464      : base(original, cloner) {
465      uniform = cloner.Clone(original.uniform);
466      mu = original.mu;
467      sigma = original.sigma;
468    }
469
470    /// <summary>
471    /// Initializes a new instance of <see cref="NormalDistributedRandom"/> with µ = 0 and sigma = 1
472    /// and a new random number generator.
473    /// </summary>
474    public NormalDistributedRandom() {
475      this.mu = 0.0;
476      this.sigma = 1.0;
477      this.uniform = new MersenneTwister();
478    }
479
480    /// <summary>
481    /// Initializes a new instance of <see cref="NormalDistributedRandom"/> with the given parameters.
482    /// <note type="caution"> No CopyConstructor! The random number generator is not copied!</note>
483    /// </summary>   
484    /// <param name="uniformRandom">The random number generator.</param>
485    /// <param name="mu">The value for µ.</param>
486    /// <param name="sigma">The value for sigma.</param>
487    public NormalDistributedRandom(IRandom uniformRandom, double mu, double sigma) {
488      this.mu = mu;
489      this.sigma = sigma;
490      this.uniform = uniformRandom;
491    }
492
493    #region IRandom Members
494
495    /// <inheritdoc cref="IRandom.Reset()"/>
496    public void Reset() {
497      uniform.Reset();
498    }
499
500    /// <inheritdoc cref="IRandom.Reset(int)"/>
501    public void Reset(int seed) {
502      uniform.Reset(seed);
503    }
504
505    /// <summary>
506    /// This method is not implemented.
507    /// </summary>
508    public int Next() {
509      throw new NotImplementedException();
510    }
511
512    /// <summary>
513    /// This method is not implemented.
514    /// </summary>
515    public int Next(int maxVal) {
516      throw new NotImplementedException();
517    }
518
519    /// <summary>
520    /// This method is not implemented.
521    /// </summary>
522    public int Next(int minVal, int maxVal) {
523      throw new NotImplementedException();
524    }
525
526    /// <summary>
527    /// Generates a new double random number.
528    /// </summary>
529    /// <returns>A double random number.</returns>
530    public double NextDouble() {
531      return NormalDistributedRandom.NextDouble(uniform, mu, sigma);
532    }
533
534    #endregion
535
536    /// <summary>
537    /// Clones the current instance (deep clone).
538    /// </summary>
539    /// <remarks>Deep clone through <see cref="cloner.Clone"/> method of helper class
540    /// <see cref="Auxiliary"/>.</remarks>
541    /// <param name="clonedObjects">Dictionary of all already cloned objects. (Needed to avoid cycles.)</param>
542    /// <returns>The cloned object as <see cref="NormalDistributedRandom"/>.</returns>
543    public override IDeepCloneable Clone(Cloner cloner) {
544      return new NormalDistributedRandom(this, cloner);
545    }
546
547    public static double NextDouble(IRandom uniformRandom, double mu, double sigma) {
548      double signFactor = uniformRandom.Next() % 2 == 0 ? 1.0 : -1.0;
549      return sigma * signFactor * NextPositiveDouble(uniformRandom) + mu;
550    }
551
552    private static double NextPositiveDouble(IRandom uniformRandom) {
553      int j = uniformRandom.Next();
554      int i = (j & 127);
555      if (Math.Abs(j) < k[i]) {
556        return j * w[i];
557      } else {
558        double r = 3.442620;
559        double x, y;
560        x = j * w[i];
561        if (i == 0) {
562          do {
563            x = -Math.Log(ScaledUniform(uniformRandom)) * 0.2904764;
564            y = -Math.Log(ScaledUniform(uniformRandom));
565          } while (y + y < x * x);
566          return (j > 0) ? r + x : -r - x;
567        }
568        if (f[i] + ScaledUniform(uniformRandom) * (f[i - 1] - f[i]) < Math.Exp(-0.5 * x * x)) {
569          return x;
570        } else {
571          // recurse
572          return (NextPositiveDouble(uniformRandom));
573        }
574      }
575    }
576
577    private static double ScaledUniform(IRandom uniformRandom) {
578      return 0.5 + uniformRandom.Next() * 0.2328306e-9;
579    }
580  }
581}
Note: See TracBrowser for help on using the repository browser.