MAJ terminée. Nous sommes passés en version 14.6.2 . Pour consulter les "releases notes" associées c'est ici :

https://about.gitlab.com/releases/2022/01/11/security-release-gitlab-14-6-2-released/
https://about.gitlab.com/releases/2022/01/04/gitlab-14-6-1-released/

FGenerateDistribution.hpp 9.48 KB
Newer Older
COULAUD Olivier's avatar
COULAUD Olivier committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
// ===================================================================================
// Copyright ScalFmm 2014 INRIA
//
// This software is a computer program whose purpose is to compute the FMM.
//
// This software is governed by the CeCILL-C and LGPL licenses and
// abiding by the rules of distribution of free software.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public and CeCILL-C Licenses for more details.
// "http://www.cecill.info".
// "http://www.gnu.org/licenses".
// ===================================================================================
16
17
18
#ifndef FGENERATEDISTRIBUTION_HPP
#define FGENERATEDISTRIBUTION_HPP

19
20
21
22
23
24
/**
 * \file
 * \brief Distribution generation implementations
 * \author O. Coulaud
 */

COULAUD Olivier's avatar
COULAUD Olivier committed
25

26
27
28
29
30
31
32
33
34
#include <cstdlib>
#include <ctime>
#include <iostream>
#include <fstream>
#include <algorithm>
//
#include "Utils/FMath.hpp"
#include "Utils/FParameters.hpp"

35
36
37
/**
 * \brief Seed the random number generator using current time
 */
38
void initRandom() {
39
40
    srand48(static_cast<long int>(time(nullptr)));
}
COULAUD Olivier's avatar
COULAUD Olivier committed
41

42
43
44
45
46
/**
 * \brief Generate a random number
 * \tparam FReal Floating point type
 * \return A random number in [0,1]
 */
47
template <class FReal>
48
49
50
FReal getRandom() {
    return static_cast<FReal>(drand48());
}
COULAUD Olivier's avatar
COULAUD Olivier committed
51

52
53
54
55
56
57
58
59
60
61
62
/**
 * \brief Generate points uniformly inside a cuboid
 *
 * \tparam FReal Floating point type
 *
 * \param N the number of points uniformly randomly sample on the unit cube
 * \param Lx the the X-length of the cuboid
 * \param Ly the the Y-length of the cuboid
 * \param Lz the the Z-length of the cuboid
 * \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0...
 */
63
template <class FReal>
64
65
66
67
68
69
70
71
72
73
void unifRandomPointsInCube(const FSize N, const FReal& Lx, const FReal& Ly,
                            const FReal& Lz, FReal* points)
{
    initRandom();
    for(FSize i = 0, j = 0 ; i< N ; ++i, j+=4)  {
        points[j]   = getRandom<FReal>() * Lx;
        points[j+1] = getRandom<FReal>() * Ly;
        points[j+2] = getRandom<FReal>() * Lz;
    }
}
74

75
76
77
78
79
80
81
82
83
84
85
/**
 * \brief Generate points uniformly inside a ball
 *
 * \tparam FReal Floating point type
 *
 * \param R the ball radius
 * \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0...
 */
template<class FReal>
void unifRandomPointsInBall(const FSize N, const FReal R, FReal* points) {
    initRandom();
86

87
88
89
    auto is_in_sphere = [&R](FReal* p) {
        return p[0]*p[0] + p[1]*p[1] + p[2]*p[2] < R*R;
    };
90

91
92
93
94
95
96
97
98
    for(FSize i = 0, j = 0 ; i< N ; ++i, j+=4)  {
        do {
            points[j]   = (getRandom<FReal>() - 0.5) * 2 * R;
            points[j+1] = (getRandom<FReal>() - 0.5) * 2 * R;
            points[j+2] = (getRandom<FReal>() - 0.5) * 2 * R;
        } while(! is_in_sphere(points + j));
    }
}
99

100
101
102
103
104
105
106
107
108
109
110
/**
 * \brief Generate N points non uniformly distributed on the ellipsoid of aspect ratio a:b:c
 *
 * \tparam FReal Floating point type
 *
 * \param N the number of points
 * \param a  the x semi-axe length
 * \param b  the y semi-axe length
 * \param c  the z semi-axe length
 * \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
 */
111
template <class FReal>
112
113
114
115
116
117
118
119
120
121
122
123
124
void nonunifRandomPointsOnElipsoid(const FSize N, const FReal& a, const FReal& b,
                                   const FReal& c, FReal* points)
{
    FReal u, v, cosu;
    for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4)  {
        u = FMath::FPi<FReal>() * getRandom<FReal>() - FMath::FPiDiv2<FReal>();
        v = FMath::FTwoPi<FReal>() * getRandom<FReal>() - FMath::FPi<FReal>();
        cosu = FMath::Cos(u);
        points[j]   = a * cosu * FMath::Cos(v);
        points[j+1] = b * cosu * FMath::Sin(v);
        points[j+2] = c * FMath::Sin(u);
    }
}
125
126


127
128
129
130
131
132
133
134
135
136
/**
 * \brief Generate N points uniformly distributed on the ellipsoid of aspect ratio a:a:c
 *
 * \tparam FReal Floating point type
 *
 * \param N the number of points
 * \param a  the x  semi-axe length
 * \param c  the z  semi-axe length
 * \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
*/
137
template <class FReal>
138
139
140
141
142
143
144
void unifRandomPointsOnProlate(const FSize N, const FReal& a, const FReal& c,
                               FReal* points)
{
    FReal u, w, v, ksi;
    FReal e = (a*a*a*a)/(c*c*c*c);
    bool isgood = false;
    FSize cpt = 0;
145

146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
    for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4)  {
        // Select a random point on the prolate
        do {
            ++cpt;
            u = 2.0 * getRandom<FReal>() - 1.0;
            v = FMath::FTwoPi<FReal>() * getRandom<FReal>();
            w = FMath::Sqrt(1 - u*u);
            points[j]	= a * w * FMath::Cos(v);
            points[j+1] = a * w * FMath::Sin(v);
            points[j+2] = c * u;
            // Accept the position ?
            ksi = a * getRandom<FReal>();
            isgood = (points[j]*points[j]
                      + points[j+1]*points[j+1]
                      + e*points[j+2]*points[j+2]) < ksi*ksi;
        } while(isgood);
    }
    std::cout.precision(4);
    std::cout << "Total tested points: " << cpt
              << " % of rejected points: "
              << 100 * static_cast<FReal>(cpt-N) / static_cast<FReal>(cpt) << " %"
              << std::endl;
}
169
170


171
172
173
174
175
176
177
178
179
180
181
/**
 * \brief  Generate N points uniformly distributed on the hyperbolic paraboloid of  aspect ratio a:b:c
 *
 * \tparam FReal Floating point type
 *
 * \param N the number of points
 * \param a  the x  semi-axe length
 * \param b  the y  semi-axe length
 * \param c  the z  semi-axe length
 * \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0...
 */
182
template <class FReal>
183
184
185
186
187
188
189
190
191
192
void unifRandomPointsOnHyperPara(const FSize N, const FReal &a, const FReal &b,
                                 const FReal &c, FReal * points)
{
    FReal u, v;
    for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4)  {
        u = 2.0 * getRandom<FReal>() - 1.0;
        v = 2.0 * getRandom<FReal>() - 1.0;
        points[j]   = a * u;
        points[j+1] = b * v;
        points[j+2] = c * (u*u - v*v);
193
194
195
196
    }
};


197
198
199
200
201
202
203
204
205
/**
 * \brief Generate N points uniformly distributed on the sphere of radius R
 *
 * \tparam FReal Floating point type
 *
 * \param N the number of points uniformly randomly sample on the sphere
 * \param R the radius of the sphere
 * \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0...
 */
206
template <class FReal>
207
208
209
210
211
212
213
214
215
216
217
218
219
220
void unifRandomPointsOnSphere(const FSize N, const FReal R, FReal* points) {
    initRandom();
    FReal u, v, theta, phi, sinPhi;
    for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4)  {
        u = getRandom<FReal>();
        v = getRandom<FReal>();
        theta  = FMath::FTwoPi<FReal>() * u;
        phi    = FMath::ACos(2*v - 1);
        sinPhi = FMath::Sin(phi);

        points[j]   = FMath::Cos(theta) * sinPhi * R;
        points[j+1] = FMath::Sin(theta) * sinPhi * R;
        points[j+2] = (2*v - 1) * R;
    }
221
222
223
};


224
225
226
227
228
229
230
231
232
/**
 * \brief Radial Plummer distribution
 *
 * \tparam FReal Floating point type
 *
 * \param cpt counter to know how many random selections we need to obtain a radius less than R
 * \param R   radius of the sphere that contains the particles
 * \return The radius according to the Plummer distribution
 */
233
template <class FReal>
234
235
236
237
238
239
240
241
242
243
FReal plummerDist(FSize& cpt, const FReal &R) {
    FReal radius, u;
    while(true) {
        u = FMath::pow(getRandom<FReal>(), 2.0/3.0);
        radius = FMath::Sqrt(u/(1.0-u));
        cpt++;
        if(radius <= R) {
            return static_cast<FReal>(radius);
        }
    }
244
245
}

246
247
248
249
/**
 * \brief Build N points following the Plummer distribution
 *
 * First we construct N points uniformly distributed on the unit sphere. Then
250
251
 * the radius in construct according to the Plummer distribution for
 * a  constant mass of 1/N
252
253
254
255
256
257
258
 *
 * \tparam FReal Floating point type
 *
 * \param N the number of points following the Plummer distribution
 * \param R the radius of the sphere that contains all the points
 * \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
 */
259
template <class FReal>
260
261
void unifRandomPlummer(const FSize N, const FReal R, FReal * points) {
    unifRandomPointsOnSphere<FReal>(N, 1, points);
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
    FReal mc = 1.0/static_cast<FReal>(N);
    for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4)  {
    	FReal m = getRandom<FReal>();
    	FReal r = FMath::Sqrt( 1.0/(FMath::pow(m, -2.0/3.0) - 1.0)) ;
        points[j]    *= r;
        points[j+1]  *= r;
        points[j+2]  *= r;
        points[j+3]   = mc;  // the mass
    }
}
/**
 * \brief Build N points following the Plummer like distribution
 *
 * First we construct N points uniformly distributed on the unit sphere. Then
 * the radius in construct according to the Plummer like distribution.
 *
 * \tparam FReal Floating point type
 *
 * \param N the number of points following the Plummer distribution
 * \param R the radius of the sphere that contains all the points
 * \param points array of size 4*N and stores data as follow x,y,z,0,x,y,z,0....
 */
template <class FReal>
void unifRandomPlummerLike(const FSize N, const FReal R, FReal * points) {
	FReal a = 1.0 ;
    unifRandomPointsOnSphere<FReal>(N, 1, points);
288
289
290
291
292
293
294
295
296
    FReal r, rm = 0.0;
    FSize cpt = 0;
    for (FSize i = 0, j = 0 ; i< N ; ++i, j+=4)  {
        r = plummerDist(cpt,R);
        rm = std::max(rm, r);
        points[j]    *= r;
        points[j+1]  *= r;
        points[j+2]  *= r;
    }
297

298
299
300
301
302
    std::cout << "Total tested points: " << cpt << " % of rejected points: "
              << 100 * static_cast<FReal>(cpt-N) / static_cast<FReal>(cpt)
              << " %"
              << std::endl;
}
COULAUD Olivier's avatar
COULAUD Olivier committed
303
//
304
#endif