1 // =============================================================================
2 // PROJECT CHRONO - http://projectchrono.org
3 //
4 // Copyright (c) 2014 projectchrono.org
5 // All rights reserved.
6 //
7 // Use of this source code is governed by a BSD-style license that can be found
8 // in the LICENSE file at the top level of the distribution and at
9 // http://projectchrono.org/license-chrono.txt.
10 //
11 // =============================================================================
12 // Authors: Alessandro Tasora, Radu Serban
13 // =============================================================================
14
15 #include "chrono/core/ChDistribution.h"
16 #include "chrono/core/ChMathematics.h"
17
18 namespace chrono {
19
ChMinMaxDistribution(double mmin,double mmax)20 ChMinMaxDistribution::ChMinMaxDistribution(double mmin, double mmax) : min(mmin), max(mmax) {}
21
GetRandom()22 double ChMinMaxDistribution::GetRandom() {
23 return min + (::chrono::ChRandom() * (max - min));
24 }
25
26 // -----------------------------------------------------------------------------
27
ChNormalDistribution(double mvariance,double mmean)28 ChNormalDistribution::ChNormalDistribution(double mvariance, double mmean)
29 : variance(mvariance), mean(mmean), hasSpare(false) {}
30
GetRandom()31 double ChNormalDistribution::GetRandom() {
32 if (hasSpare) {
33 hasSpare = false;
34 return (mean + sqrt(variance * rand1) * sin(rand2));
35 }
36
37 hasSpare = true;
38
39 rand1 = ::chrono::ChRandom();
40 if (rand1 < 1e-100)
41 rand1 = 1e-100;
42 rand1 = -2.0 * log(rand1);
43 rand2 = ::chrono::ChRandom() * CH_C_2PI;
44
45 return (mean + sqrt(variance * rand1) * cos(rand2));
46 }
47
48 // -----------------------------------------------------------------------------
49
ChWeibullDistribution(double mlambda,double mk)50 ChWeibullDistribution::ChWeibullDistribution(double mlambda, double mk) : lambda(mlambda), k(mk) {}
51
GetRandom()52 double ChWeibullDistribution::GetRandom() {
53 double rand = ::chrono::ChRandom();
54 if (rand < 1e-100)
55 rand = 1e-100;
56 return lambda * pow(-log(rand), (1.0 / k));
57 }
58
59 // -----------------------------------------------------------------------------
60
ChZhangDistribution(double average_size,double minimum_size)61 ChZhangDistribution::ChZhangDistribution(double average_size, double minimum_size) : minsize(minimum_size) {
62 lambdar = 1.0 / (average_size - minimum_size);
63 }
64
GetRandom()65 double ChZhangDistribution::GetRandom() {
66 double rand = ::chrono::ChRandom();
67 if (rand < 1e-100)
68 rand = 1e-100;
69 return (minsize + (1.0 / lambdar) * (-log(rand)));
70 }
71
72 // -----------------------------------------------------------------------------
73
ChContinuumDistribution(ChVectorDynamic<> & mx,ChVectorDynamic<> & my)74 ChContinuumDistribution::ChContinuumDistribution(ChVectorDynamic<>& mx, ChVectorDynamic<>& my) : x(mx), y(my) {
75 if (mx.size() != my.size())
76 throw ChException("Probability curve must have same number of elements in abscysse and ordinates");
77
78 cdf_x = mx;
79 cdf_y = my;
80
81 // compute CDF
82 double integral = 0;
83 for (int i = 0; i < x.size() - 1; i++) {
84 integral += 0.5 * (y(i) + y(i + 1)) * (x(i + 1) - x(i));
85 cdf_y(i) = integral;
86 cdf_x(i) = 0.5 * (x(i + 1) + x(i));
87 }
88
89 // normalize if P(x) had not unit integral
90 double totintegral = cdf_y(x.size() - 2);
91 if (totintegral != 1.0) {
92 for (int i = 0; i < x.size() - 1; i++) {
93 cdf_y(i) *= 1. / totintegral;
94 }
95 }
96
97 cdf_x(x.size() - 1) = x(x.size() - 1);
98 cdf_y(x.size() - 1) = 1.0;
99 }
100
GetRandom()101 double ChContinuumDistribution::GetRandom() {
102 double mx1 = x(0);
103 double mx2 = cdf_x(0);
104 double my1 = 0;
105 double my2 = cdf_y(0);
106
107 double rand = ChRandom();
108 for (int i = 0; i < x.size() - 1; i++) {
109 if ((rand <= cdf_y(i + 1)) && (rand > cdf_y(i))) {
110 mx1 = cdf_x(i);
111 mx2 = cdf_x(i + 1);
112 my1 = cdf_y(i);
113 my2 = cdf_y(i + 1);
114 break;
115 }
116 }
117 // linear interp
118 double val = mx1 + ((rand - my1) / (my2 - my1)) * (mx2 - mx1);
119 return val;
120 }
121
122 // -----------------------------------------------------------------------------
123
ChDiscreteDistribution(ChVectorDynamic<> & mx,ChVectorDynamic<> & my)124 ChDiscreteDistribution::ChDiscreteDistribution(ChVectorDynamic<>& mx, ChVectorDynamic<>& my) : x(mx), y(my) {
125 if (mx.size() != my.size())
126 throw ChException("Probability values and percentages must have the same size");
127
128 cdf_y = my;
129
130 // compute CDF
131 double integral = 0;
132 for (int i = 0; i < x.size(); i++) {
133 integral += y(i);
134 cdf_y(i) = integral;
135 }
136 // normalize if P(x) had not unit integral
137 double totintegral = cdf_y(x.size() - 1);
138 if (totintegral != 1.0) {
139 for (int i = 0; i < x.size(); i++) {
140 cdf_y(i) *= 1. / totintegral;
141 }
142 }
143 cdf_y(x.size() - 1) = 1.0;
144 }
145
GetRandom()146 double ChDiscreteDistribution::GetRandom() {
147 double rand = ChRandom();
148 double lastval = 0;
149 for (int i = 0; i < x.size(); i++) {
150 if ((rand <= cdf_y(i)) && (rand > lastval)) {
151 return x(i);
152 }
153 lastval = cdf_y(i);
154 }
155 return 0;
156 }
157
158 } // end namespace chrono
159