9 # include "discrepancy.h"
12 namespace pagmo{
namespace util {
namespace discrepancy {
28 double f = 1.0 / base;
31 retval += f * (i % base);
38 # define PRIME_MAX 1601
50 unsigned int npvec[PRIME_MAX] = { 1,
51 2, 3, 5, 7, 11, 13, 17, 19, 23, 29,
52 31, 37, 41, 43, 47, 53, 59, 61, 67, 71,
53 73, 79, 83, 89, 97, 101, 103, 107, 109, 113,
54 127, 131, 137, 139, 149, 151, 157, 163, 167, 173,
55 179, 181, 191, 193, 197, 199, 211, 223, 227, 229,
56 233, 239, 241, 251, 257, 263, 269, 271, 277, 281,
57 283, 293, 307, 311, 313, 317, 331, 337, 347, 349,
58 353, 359, 367, 373, 379, 383, 389, 397, 401, 409,
59 419, 421, 431, 433, 439, 443, 449, 457, 461, 463,
60 467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
61 547, 557, 563, 569, 571, 577, 587, 593, 599, 601,
62 607, 613, 617, 619, 631, 641, 643, 647, 653, 659,
63 661, 673, 677, 683, 691, 701, 709, 719, 727, 733,
64 739, 743, 751, 757, 761, 769, 773, 787, 797, 809,
65 811, 821, 823, 827, 829, 839, 853, 857, 859, 863,
66 877, 881, 883, 887, 907, 911, 919, 929, 937, 941,
67 947, 953, 967, 971, 977, 983, 991, 997, 1009, 1013,
68 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069,
69 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151,
70 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223,
71 1229, 1231, 1237, 1249, 1259, 1277, 1279, 1283, 1289, 1291,
72 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373,
73 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451,
74 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511,
75 1523, 1531, 1543, 1549, 1553, 1559, 1567, 1571, 1579, 1583,
76 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657,
77 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733,
78 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811,
79 1823, 1831, 1847, 1861, 1867, 1871, 1873, 1877, 1879, 1889,
80 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987,
81 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053,
82 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129,
83 2131, 2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213,
84 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287,
85 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357,
86 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423,
87 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531,
88 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593, 2609, 2617,
89 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687,
90 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741,
91 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819,
92 2833, 2837, 2843, 2851, 2857, 2861, 2879, 2887, 2897, 2903,
93 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999,
94 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
95 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181,
96 3187, 3191, 3203, 3209, 3217, 3221, 3229, 3251, 3253, 3257,
97 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331,
98 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407, 3413,
99 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511,
100 3517, 3527, 3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571,
101 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643,
102 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727,
103 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821,
104 3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907,
105 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967, 3989,
106 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057,
107 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139,
108 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231,
109 4241, 4243, 4253, 4259, 4261, 4271, 4273, 4283, 4289, 4297,
110 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409,
111 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483, 4493,
112 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583,
113 4591, 4597, 4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657,
114 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751,
115 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813, 4817, 4831,
116 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937,
117 4943, 4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003,
118 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087,
119 5099, 5101, 5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179,
120 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279,
121 5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387,
122 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443,
123 5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521,
124 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623, 5639,
125 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693,
126 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791,
127 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857,
128 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923, 5927, 5939,
129 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053,
130 6067, 6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133,
131 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221,
132 6229, 6247, 6257, 6263, 6269, 6271, 6277, 6287, 6299, 6301,
133 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367,
134 6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473,
135 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571,
136 6577, 6581, 6599, 6607, 6619, 6637, 6653, 6659, 6661, 6673,
137 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761,
138 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833,
139 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917,
140 6947, 6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997,
141 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103,
142 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207,
143 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297,
144 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411,
145 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487, 7489, 7499,
146 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561,
147 7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643,
148 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723,
149 7727, 7741, 7753, 7757, 7759, 7789, 7793, 7817, 7823, 7829,
150 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919,
151 7927, 7933, 7937, 7949, 7951, 7963, 7993, 8009, 8011, 8017,
152 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111,
153 8117, 8123, 8147, 8161, 8167, 8171, 8179, 8191, 8209, 8219,
154 8221, 8231, 8233, 8237, 8243, 8263, 8269, 8273, 8287, 8291,
155 8293, 8297, 8311, 8317, 8329, 8353, 8363, 8369, 8377, 8387,
156 8389, 8419, 8423, 8429, 8431, 8443, 8447, 8461, 8467, 8501,
157 8513, 8521, 8527, 8537, 8539, 8543, 8563, 8573, 8581, 8597,
158 8599, 8609, 8623, 8627, 8629, 8641, 8647, 8663, 8669, 8677,
159 8681, 8689, 8693, 8699, 8707, 8713, 8719, 8731, 8737, 8741,
160 8747, 8753, 8761, 8779, 8783, 8803, 8807, 8819, 8821, 8831,
161 8837, 8839, 8849, 8861, 8863, 8867, 8887, 8893, 8923, 8929,
162 8933, 8941, 8951, 8963, 8969, 8971, 8999, 9001, 9007, 9011,
163 9013, 9029, 9041, 9043, 9049, 9059, 9067, 9091, 9103, 9109,
164 9127, 9133, 9137, 9151, 9157, 9161, 9173, 9181, 9187, 9199,
165 9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283,
166 9293, 9311, 9319, 9323, 9337, 9341, 9343, 9349, 9371, 9377,
167 9391, 9397, 9403, 9413, 9419, 9421, 9431, 9433, 9437, 9439,
168 9461, 9463, 9467, 9473, 9479, 9491, 9497, 9511, 9521, 9533,
169 9539, 9547, 9551, 9587, 9601, 9613, 9619, 9623, 9629, 9631,
170 9643, 9649, 9661, 9677, 9679, 9689, 9697, 9719, 9721, 9733,
171 9739, 9743, 9749, 9767, 9769, 9781, 9787, 9791, 9803, 9811,
172 9817, 9829, 9833, 9839, 9851, 9857, 9859, 9871, 9883, 9887,
173 9901, 9907, 9923, 9929, 9931, 9941, 9949, 9967, 9973,10007,
174 10009,10037,10039,10061,10067,10069,10079,10091,10093,10099,
175 10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,
176 10181,10193,10211,10223,10243,10247,10253,10259,10267,10271,
177 10273,10289,10301,10303,10313,10321,10331,10333,10337,10343,
178 10357,10369,10391,10399,10427,10429,10433,10453,10457,10459,
179 10463,10477,10487,10499,10501,10513,10529,10531,10559,10567,
180 10589,10597,10601,10607,10613,10627,10631,10639,10651,10657,
181 10663,10667,10687,10691,10709,10711,10723,10729,10733,10739,
182 10753,10771,10781,10789,10799,10831,10837,10847,10853,10859,
183 10861,10867,10883,10889,10891,10903,10909,10937,10939,10949,
184 10957,10973,10979,10987,10993,11003,11027,11047,11057,11059,
185 11069,11071,11083,11087,11093,11113,11117,11119,11131,11149,
186 11159,11161,11171,11173,11177,11197,11213,11239,11243,11251,
187 11257,11261,11273,11279,11287,11299,11311,11317,11321,11329,
188 11351,11353,11369,11383,11393,11399,11411,11423,11437,11443,
189 11447,11467,11471,11483,11489,11491,11497,11503,11519,11527,
190 11549,11551,11579,11587,11593,11597,11617,11621,11633,11657,
191 11677,11681,11689,11699,11701,11717,11719,11731,11743,11777,
192 11779,11783,11789,11801,11807,11813,11821,11827,11831,11833,
193 11839,11863,11867,11887,11897,11903,11909,11923,11927,11933,
194 11939,11941,11953,11959,11969,11971,11981,11987,12007,12011,
195 12037,12041,12043,12049,12071,12073,12097,12101,12107,12109,
196 12113,12119,12143,12149,12157,12161,12163,12197,12203,12211,
197 12227,12239,12241,12251,12253,12263,12269,12277,12281,12289,
198 12301,12323,12329,12343,12347,12373,12377,12379,12391,12401,
199 12409,12413,12421,12433,12437,12451,12457,12473,12479,12487,
200 12491,12497,12503,12511,12517,12527,12539,12541,12547,12553,
201 12569,12577,12583,12589,12601,12611,12613,12619,12637,12641,
202 12647,12653,12659,12671,12689,12697,12703,12713,12721,12739,
203 12743,12757,12763,12781,12791,12799,12809,12821,12823,12829,
204 12841,12853,12889,12893,12899,12907,12911,12917,12919,12923,
205 12941,12953,12959,12967,12973,12979,12983,13001,13003,13007,
206 13009,13033,13037,13043,13049,13063,13093,13099,13103,13109,
207 13121,13127,13147,13151,13159,13163,13171,13177,13183,13187,
208 13217,13219,13229,13241,13249,13259,13267,13291,13297,13309,
209 13313,13327,13331,13337,13339,13367,13381,13397,13399,13411,
210 13417,13421,13441,13451,13457,13463,13469,13477,13487,13499 };
212 return npvec[PRIME_MAX-1];
213 }
else if ( n < PRIME_MAX ) {
218 pagmo_throw(value_error,
"n must be in [-1,1600]");
236 unsigned int ub = PRIME_MAX-1;
237 unsigned int min_p =
prime(lb);
238 unsigned int max_p =
prime(ub);
239 if (n<min_p || n>max_p) {
240 pagmo_throw(value_error,
"n is out of range");
242 if (n == min_p)
return min_p;
243 if (n == max_p)
return max_p;
245 unsigned int tmp = (ub+lb)/2;
246 unsigned int new_p =
prime(tmp);
261 std::vector<double> project_2_simplex::operator()(std::vector<double> retval)
const
263 if (retval.size() == m_dim-1) {
264 std::sort(retval.begin(),retval.end());
265 retval.insert(retval.begin(),0.0);
266 retval.push_back(1.0);
268 for (
unsigned int i = 0; i<retval.size()-1;++i) {
269 retval[i] = retval[i+1] - retval[i];
273 for (
unsigned int i = 0; i<retval.size();++i) {
279 pagmo_throw(value_error,
"To project on this simplex you need a point in dimension m_dim-1");
286 int *faure::binomial_table (
int qs,
int m,
int n )
325 coef =
new int[(m+1)*(n+1)];
327 for ( j = 0; j <= n; j++ )
329 for ( i = 0; i <= m; i++ )
338 for ( i = 1; i <= m; i++ )
343 for ( i = 1; i <= i4_min ( m, n ); i++ )
349 for( j = 1; j <= n; j++ )
351 for ( i = j + 1; i <= m; i++ )
353 coef[i+j*(m+1)] = ( coef[i-1+j*(m+1)] + coef[i-1+(j-1)*(m+1)] ) % qs;
361 void faure::faure_orig (
unsigned int dim_num,
unsigned int *seed,
double quasi[] )
432 if ( m_qs <= 0 || *seed <= 0 )
439 cout <<
"FAURE - Fatal error!\n";
440 cout <<
" PRIME_GE failed.\n";
452 hisum = i4_log_i4 ( *seed, m_qs );
457 if ( m_hisum_save != hisum )
459 if ( m_coef != NULL )
464 if ( m_ytemp != NULL )
469 m_hisum_save = hisum;
471 m_coef = binomial_table ( m_qs, hisum, hisum );
473 m_ytemp =
new int[hisum+1];
484 ktemp = i4_power ( m_qs, hisum + 1 );
487 for (
int i = hisum; 0 <= i; i-- )
489 ktemp = ktemp / m_qs;
490 mtemp = ltemp % ktemp;
491 m_ytemp[i] = ( ltemp - mtemp ) / ktemp;
501 r = ( ( double ) m_ytemp[hisum] );
502 for (
int i = hisum-1; 0 <= i; i-- )
504 r = ( ( double ) m_ytemp[i] ) + r / ( ( double ) m_qs );
507 quasi[0] = r / ( ( double ) m_qs );
511 for (
unsigned int k = 1; k < dim_num; k++ )
514 r = 1.0 / ( ( double ) m_qs );
516 for (
int j = 0; j <= hisum; j++ )
519 for (
int i = j; i <= hisum; i++ )
521 ztemp = ztemp + m_ytemp[i] * m_coef[i+j*(hisum+1)];
528 m_ytemp[j] = ztemp % m_qs;
529 quasi[k] = quasi[k] + ( ( double ) m_ytemp[j] ) * r;
530 r = r / ( ( double ) m_qs );
542 int faure::i4_log_i4 (
int i4,
int j4 )
621 while ( j4_abs <= i4_abs )
623 i4_abs = i4_abs / j4_abs;
632 int faure::i4_min (
int i1,
int i2 )
673 int faure::i4_power (
int i,
int j )
712 cout <<
"I4_POWER - Fatal error!\n";
713 cout <<
" I^J requested, with I = 0 and J negative.\n";
726 cout <<
"I4_POWER - Fatal error!\n";
727 cout <<
" I^J requested, with I = 0 and J = 0.\n";
742 for ( k = 1; k <= j; k++ )
752 int sobol::i8_bit_lo0 (
long long int n )
829 void sobol::i8_sobol (
unsigned int dim_num,
long long int *seed,
double quasi[ ] )
915 # define DIM_MAX2 1111
928 bool includ[LOG_MAX];
938 long long int seed_temp;
940 if ( !m_initialized || dim_num != m_dim_num_save )
942 m_initialized =
true;
943 long long int polyb[DIM_MAX2] =
945 1, 3, 7, 11, 13, 19, 25, 37, 59, 47,
946 61, 55, 41, 67, 97, 91, 109, 103, 115, 131,
947 193, 137, 145, 143, 241, 157, 185, 167, 229, 171,
948 213, 191, 253, 203, 211, 239, 247, 285, 369, 299,
949 301, 333, 351, 355, 357, 361, 391, 397, 425, 451,
950 463, 487, 501, 529, 539, 545, 557, 563, 601, 607,
951 617, 623, 631, 637, 647, 661, 675, 677, 687, 695,
952 701, 719, 721, 731, 757, 761, 787, 789, 799, 803,
953 817, 827, 847, 859, 865, 875, 877, 883, 895, 901,
954 911, 949, 953, 967, 971, 973, 981, 985, 995, 1001,
955 1019, 1033, 1051, 1063, 1069, 1125, 1135, 1153, 1163, 1221,
956 1239, 1255, 1267, 1279, 1293, 1305, 1315, 1329, 1341, 1347,
957 1367, 1387, 1413, 1423, 1431, 1441, 1479, 1509, 1527, 1531,
958 1555, 1557, 1573, 1591, 1603, 1615, 1627, 1657, 1663, 1673,
959 1717, 1729, 1747, 1759, 1789, 1815, 1821, 1825, 1849, 1863,
960 1869, 1877, 1881, 1891, 1917, 1933, 1939, 1969, 2011, 2035,
961 2041, 2053, 2071, 2091, 2093, 2119, 2147, 2149, 2161, 2171,
962 2189, 2197, 2207, 2217, 2225, 2255, 2257, 2273, 2279, 2283,
963 2293, 2317, 2323, 2341, 2345, 2363, 2365, 2373, 2377, 2385,
964 2395, 2419, 2421, 2431, 2435, 2447, 2475, 2477, 2489, 2503,
965 2521, 2533, 2551, 2561, 2567, 2579, 2581, 2601, 2633, 2657,
966 2669, 2681, 2687, 2693, 2705, 2717, 2727, 2731, 2739, 2741,
967 2773, 2783, 2793, 2799, 2801, 2811, 2819, 2825, 2833, 2867,
968 2879, 2881, 2891, 2905, 2911, 2917, 2927, 2941, 2951, 2955,
969 2963, 2965, 2991, 2999, 3005, 3017, 3035, 3037, 3047, 3053,
970 3083, 3085, 3097, 3103, 3159, 3169, 3179, 3187, 3205, 3209,
971 3223, 3227, 3229, 3251, 3263, 3271, 3277, 3283, 3285, 3299,
972 3305, 3319, 3331, 3343, 3357, 3367, 3373, 3393, 3399, 3413,
973 3417, 3427, 3439, 3441, 3475, 3487, 3497, 3515, 3517, 3529,
974 3543, 3547, 3553, 3559, 3573, 3589, 3613, 3617, 3623, 3627,
975 3635, 3641, 3655, 3659, 3669, 3679, 3697, 3707, 3709, 3713,
976 3731, 3743, 3747, 3771, 3791, 3805, 3827, 3833, 3851, 3865,
977 3889, 3895, 3933, 3947, 3949, 3957, 3971, 3985, 3991, 3995,
978 4007, 4013, 4021, 4045, 4051, 4069, 4073, 4179, 4201, 4219,
979 4221, 4249, 4305, 4331, 4359, 4383, 4387, 4411, 4431, 4439,
980 4449, 4459, 4485, 4531, 4569, 4575, 4621, 4663, 4669, 4711,
981 4723, 4735, 4793, 4801, 4811, 4879, 4893, 4897, 4921, 4927,
982 4941, 4977, 5017, 5027, 5033, 5127, 5169, 5175, 5199, 5213,
983 5223, 5237, 5287, 5293, 5331, 5391, 5405, 5453, 5523, 5573,
984 5591, 5597, 5611, 5641, 5703, 5717, 5721, 5797, 5821, 5909,
985 5913, 5955, 5957, 6005, 6025, 6061, 6067, 6079, 6081, 6231,
986 6237, 6289, 6295, 6329, 6383, 6427, 6453, 6465, 6501, 6523,
987 6539, 6577, 6589, 6601, 6607, 6631, 6683, 6699, 6707, 6761,
988 6795, 6865, 6881, 6901, 6923, 6931, 6943, 6999, 7057, 7079,
989 7103, 7105, 7123, 7173, 7185, 7191, 7207, 7245, 7303, 7327,
990 7333, 7355, 7365, 7369, 7375, 7411, 7431, 7459, 7491, 7505,
991 7515, 7541, 7557, 7561, 7701, 7705, 7727, 7749, 7761, 7783,
992 7795, 7823, 7907, 7953, 7963, 7975, 8049, 8089, 8123, 8125,
993 8137, 8219, 8231, 8245, 8275, 8293, 8303, 8331, 8333, 8351,
994 8357, 8367, 8379, 8381, 8387, 8393, 8417, 8435, 8461, 8469,
995 8489, 8495, 8507, 8515, 8551, 8555, 8569, 8585, 8599, 8605,
996 8639, 8641, 8647, 8653, 8671, 8675, 8689, 8699, 8729, 8741,
997 8759, 8765, 8771, 8795, 8797, 8825, 8831, 8841, 8855, 8859,
998 8883, 8895, 8909, 8943, 8951, 8955, 8965, 8999, 9003, 9031,
999 9045, 9049, 9071, 9073, 9085, 9095, 9101, 9109, 9123, 9129,
1000 9137, 9143, 9147, 9185, 9197, 9209, 9227, 9235, 9247, 9253,
1001 9257, 9277, 9297, 9303, 9313, 9325, 9343, 9347, 9371, 9373,
1002 9397, 9407, 9409, 9415, 9419, 9443, 9481, 9495, 9501, 9505,
1003 9517, 9529, 9555, 9557, 9571, 9585, 9591, 9607, 9611, 9621,
1004 9625, 9631, 9647, 9661, 9669, 9679, 9687, 9707, 9731, 9733,
1005 9745, 9773, 9791, 9803, 9811, 9817, 9833, 9847, 9851, 9863,
1006 9875, 9881, 9905, 9911, 9917, 9923, 9963, 9973,10003,10025,
1007 10043,10063,10071,10077,10091,10099,10105,10115,10129,10145,
1008 10169,10183,10187,10207,10223,10225,10247,10265,10271,10275,
1009 10289,10299,10301,10309,10343,10357,10373,10411,10413,10431,
1010 10445,10453,10463,10467,10473,10491,10505,10511,10513,10523,
1011 10539,10549,10559,10561,10571,10581,10615,10621,10625,10643,
1012 10655,10671,10679,10685,10691,10711,10739,10741,10755,10767,
1013 10781,10785,10803,10805,10829,10857,10863,10865,10875,10877,
1014 10917,10921,10929,10949,10967,10971,10987,10995,11009,11029,
1015 11043,11045,11055,11063,11075,11081,11117,11135,11141,11159,
1016 11163,11181,11187,11225,11237,11261,11279,11297,11307,11309,
1017 11327,11329,11341,11377,11403,11405,11413,11427,11439,11453,
1018 11461,11473,11479,11489,11495,11499,11533,11545,11561,11567,
1019 11575,11579,11589,11611,11623,11637,11657,11663,11687,11691,
1020 11701,11747,11761,11773,11783,11795,11797,11817,11849,11855,
1021 11867,11869,11873,11883,11919,11921,11927,11933,11947,11955,
1022 11961,11999,12027,12029,12037,12041,12049,12055,12095,12097,
1023 12107,12109,12121,12127,12133,12137,12181,12197,12207,12209,
1024 12239,12253,12263,12269,12277,12287,12295,12309,12313,12335,
1025 12361,12367,12391,12409,12415,12433,12449,12469,12479,12481,
1026 12499,12505,12517,12527,12549,12559,12597,12615,12621,12639,
1027 12643,12657,12667,12707,12713,12727,12741,12745,12763,12769,
1028 12779,12781,12787,12799,12809,12815,12829,12839,12857,12875,
1029 12883,12889,12901,12929,12947,12953,12959,12969,12983,12987,
1030 12995,13015,13019,13031,13063,13077,13103,13137,13149,13173,
1031 13207,13211,13227,13241,13249,13255,13269,13283,13285,13303,
1032 13307,13321,13339,13351,13377,13389,13407,13417,13431,13435,
1033 13447,13459,13465,13477,13501,13513,13531,13543,13561,13581,
1034 13599,13605,13617,13623,13637,13647,13661,13677,13683,13695,
1035 13725,13729,13753,13773,13781,13785,13795,13801,13807,13825,
1036 13835,13855,13861,13871,13883,13897,13905,13915,13939,13941,
1037 13969,13979,13981,13997,14027,14035,14037,14051,14063,14085,
1038 14095,14107,14113,14125,14137,14145,14151,14163,14193,14199,
1039 14219,14229,14233,14243,14277,14287,14289,14295,14301,14305,
1040 14323,14339,14341,14359,14365,14375,14387,14411,14425,14441,
1041 14449,14499,14513,14523,14537,14543,14561,14579,14585,14593,
1042 14599,14603,14611,14641,14671,14695,14701,14723,14725,14743,
1043 14753,14759,14765,14795,14797,14803,14831,14839,14845,14855,
1044 14889,14895,14909,14929,14941,14945,14951,14963,14965,14985,
1045 15033,15039,15053,15059,15061,15071,15077,15081,15099,15121,
1046 15147,15149,15157,15167,15187,15193,15203,15205,15215,15217,
1047 15223,15243,15257,15269,15273,15287,15291,15313,15335,15347,
1048 15359,15373,15379,15381,15391,15395,15397,15419,15439,15453,
1049 15469,15491,15503,15517,15527,15531,15545,15559,15593,15611,
1050 15613,15619,15639,15643,15649,15661,15667,15669,15681,15693,
1051 15717,15721,15741,15745,15765,15793,15799,15811,15825,15835,
1052 15847,15851,15865,15877,15881,15887,15899,15915,15935,15937,
1053 15955,15973,15977,16011,16035,16061,16069,16087,16093,16097,
1054 16121,16141,16153,16159,16165,16183,16189,16195,16197,16201,
1055 16209,16215,16225,16259,16265,16273,16299,16309,16355,16375,
1057 for ( i = 0; i < DIM_MAX2; i++ )
1060 for ( j = 0; j < LOG_MAX; j++ )
12787 v[1000][10] = 1051;
12788 v[1001][10] = 1261;
12790 v[1003][10] = 1555;
12791 v[1004][10] = 1757;
12792 v[1005][10] = 1777;
12794 v[1007][10] = 1583;
12795 v[1008][10] = 1957;
12798 v[1011][10] = 1163;
12801 v[1014][10] = 1963;
12803 v[1016][10] = 1905;
12805 v[1018][10] = 1677;
12809 v[1022][10] = 1723;
12811 v[1024][10] = 1885;
12812 v[1025][10] = 1249;
12814 v[1027][10] = 1803;
12819 v[1032][10] = 1767;
12822 v[1035][10] = 1035;
12824 v[1037][10] = 1263;
12825 v[1038][10] = 1881;
12826 v[1039][10] = 1779;
12827 v[1040][10] = 1565;
12832 v[1045][10] = 1419;
12834 v[1047][10] = 1889;
12836 v[1049][10] = 1871;
12837 v[1050][10] = 1869;
12840 v[1053][10] = 1547;
12841 v[1054][10] = 1799;
12843 v[1056][10] = 1441;
12845 v[1058][10] = 2021;
12846 v[1059][10] = 1303;
12847 v[1060][10] = 1505;
12848 v[1061][10] = 1735;
12849 v[1062][10] = 1619;
12850 v[1063][10] = 1065;
12851 v[1064][10] = 1161;
12852 v[1065][10] = 2047;
12856 v[1069][10] = 1447;
12859 v[1072][10] = 1065;
12863 v[1076][10] = 1257;
12864 v[1077][10] = 1833;
12866 v[1079][10] = 1841;
12867 v[1080][10] = 1733;
12868 v[1081][10] = 1179;
12869 v[1082][10] = 1191;
12870 v[1083][10] = 1025;
12871 v[1084][10] = 1639;
12872 v[1085][10] = 1955;
12873 v[1086][10] = 1423;
12874 v[1087][10] = 1685;
12875 v[1088][10] = 1711;
12879 v[1092][10] = 1653;
12884 v[1097][10] = 1505;
12886 v[1099][10] = 1449;
12887 v[1100][10] = 1573;
12888 v[1101][10] = 1297;
12889 v[1102][10] = 1821;
12890 v[1103][10] = 1691;
12893 v[1106][10] = 1187;
12895 v[1108][10] = 1535;
13562 v[1000][11] = 2977;
13563 v[1001][11] = 1979;
13564 v[1002][11] = 2271;
13565 v[1003][11] = 3247;
13566 v[1004][11] = 1267;
13567 v[1005][11] = 1747;
13571 v[1009][11] = 2001;
13572 v[1010][11] = 1195;
13573 v[1011][11] = 3065;
13575 v[1013][11] = 1499;
13576 v[1014][11] = 3529;
13577 v[1015][11] = 1081;
13578 v[1016][11] = 2877;
13579 v[1017][11] = 3077;
13581 v[1019][11] = 1793;
13582 v[1020][11] = 2409;
13583 v[1021][11] = 3995;
13584 v[1022][11] = 2559;
13585 v[1023][11] = 4081;
13586 v[1024][11] = 1195;
13587 v[1025][11] = 2955;
13588 v[1026][11] = 1117;
13589 v[1027][11] = 1409;
13592 v[1030][11] = 1521;
13593 v[1031][11] = 1607;
13595 v[1033][11] = 3055;
13596 v[1034][11] = 3123;
13597 v[1035][11] = 2533;
13598 v[1036][11] = 2329;
13599 v[1037][11] = 3477;
13601 v[1039][11] = 3683;
13602 v[1040][11] = 3715;
13604 v[1042][11] = 3139;
13605 v[1043][11] = 3311;
13607 v[1045][11] = 3511;
13608 v[1046][11] = 2299;
13610 v[1048][11] = 2941;
13611 v[1049][11] = 3067;
13612 v[1050][11] = 1331;
13613 v[1051][11] = 1081;
13614 v[1052][11] = 1097;
13615 v[1053][11] = 2853;
13616 v[1054][11] = 2299;
13618 v[1056][11] = 1745;
13620 v[1058][11] = 3819;
13622 v[1060][11] = 1059;
13623 v[1061][11] = 3559;
13625 v[1063][11] = 3743;
13628 v[1066][11] = 3501;
13630 v[1068][11] = 2599;
13631 v[1069][11] = 3983;
13632 v[1070][11] = 3961;
13634 v[1072][11] = 1899;
13636 v[1074][11] = 2493;
13637 v[1075][11] = 1795;
13641 v[1079][11] = 2361;
13642 v[1080][11] = 3093;
13643 v[1081][11] = 3119;
13644 v[1082][11] = 3679;
13645 v[1083][11] = 2367;
13646 v[1084][11] = 1701;
13647 v[1085][11] = 1445;
13648 v[1086][11] = 1321;
13649 v[1087][11] = 2397;
13650 v[1088][11] = 1241;
13651 v[1089][11] = 3305;
13652 v[1090][11] = 3985;
13653 v[1091][11] = 2349;
13654 v[1092][11] = 4067;
13655 v[1093][11] = 3805;
13656 v[1094][11] = 3073;
13657 v[1095][11] = 2837;
13658 v[1096][11] = 1567;
13659 v[1097][11] = 3783;
13661 v[1099][11] = 2441;
13662 v[1100][11] = 1181;
13665 v[1103][11] = 1201;
13666 v[1104][11] = 3735;
13667 v[1105][11] = 2517;
13669 v[1107][11] = 1535;
13670 v[1108][11] = 2175;
13671 v[1109][11] = 3613;
13672 v[1110][11] = 3019;
14193 v[1000][12] = 4939;
14194 v[1001][12] = 7671;
14195 v[1002][12] = 6059;
14196 v[1003][12] = 6275;
14197 v[1004][12] = 6517;
14198 v[1005][12] = 1931;
14199 v[1006][12] = 4583;
14200 v[1007][12] = 7301;
14201 v[1008][12] = 1267;
14202 v[1009][12] = 7509;
14203 v[1010][12] = 1435;
14204 v[1011][12] = 2169;
14205 v[1012][12] = 6939;
14206 v[1013][12] = 3515;
14207 v[1014][12] = 2985;
14208 v[1015][12] = 2787;
14209 v[1016][12] = 2123;
14210 v[1017][12] = 1969;
14211 v[1018][12] = 3307;
14213 v[1020][12] = 4359;
14214 v[1021][12] = 7059;
14215 v[1022][12] = 5273;
14216 v[1023][12] = 5873;
14217 v[1024][12] = 6657;
14218 v[1025][12] = 6765;
14219 v[1026][12] = 6229;
14220 v[1027][12] = 3179;
14221 v[1028][12] = 1583;
14222 v[1029][12] = 6237;
14223 v[1030][12] = 2155;
14226 v[1033][12] = 7491;
14227 v[1034][12] = 3309;
14228 v[1035][12] = 6805;
14229 v[1036][12] = 3015;
14230 v[1037][12] = 6831;
14231 v[1038][12] = 7819;
14233 v[1040][12] = 4747;
14234 v[1041][12] = 3935;
14235 v[1042][12] = 4109;
14236 v[1043][12] = 1311;
14238 v[1045][12] = 3089;
14239 v[1046][12] = 7059;
14240 v[1047][12] = 4247;
14241 v[1048][12] = 2989;
14242 v[1049][12] = 1509;
14243 v[1050][12] = 4919;
14244 v[1051][12] = 1841;
14245 v[1052][12] = 3045;
14246 v[1053][12] = 3821;
14247 v[1054][12] = 6929;
14248 v[1055][12] = 4655;
14249 v[1056][12] = 1333;
14250 v[1057][12] = 6429;
14251 v[1058][12] = 6649;
14252 v[1059][12] = 2131;
14253 v[1060][12] = 5265;
14254 v[1061][12] = 1051;
14256 v[1063][12] = 8057;
14257 v[1064][12] = 3379;
14258 v[1065][12] = 2179;
14259 v[1066][12] = 1993;
14260 v[1067][12] = 5655;
14261 v[1068][12] = 3063;
14262 v[1069][12] = 6381;
14263 v[1070][12] = 3587;
14264 v[1071][12] = 7417;
14265 v[1072][12] = 1579;
14266 v[1073][12] = 1541;
14267 v[1074][12] = 2107;
14268 v[1075][12] = 5085;
14269 v[1076][12] = 2873;
14270 v[1077][12] = 6141;
14272 v[1079][12] = 3537;
14273 v[1080][12] = 2157;
14275 v[1082][12] = 1999;
14276 v[1083][12] = 1465;
14277 v[1084][12] = 5171;
14278 v[1085][12] = 5651;
14279 v[1086][12] = 1535;
14280 v[1087][12] = 7235;
14281 v[1088][12] = 4349;
14282 v[1089][12] = 1263;
14283 v[1090][12] = 1453;
14284 v[1091][12] = 1005;
14285 v[1092][12] = 6893;
14286 v[1093][12] = 2919;
14287 v[1094][12] = 1947;
14288 v[1095][12] = 1635;
14289 v[1096][12] = 3963;
14292 v[1099][12] = 4569;
14294 v[1101][12] = 6737;
14295 v[1102][12] = 2995;
14296 v[1103][12] = 7235;
14297 v[1104][12] = 7713;
14299 v[1106][12] = 4821;
14300 v[1107][12] = 2377;
14301 v[1108][12] = 1673;
14303 v[1110][12] = 6541;
14307 if ( dim_num < 1 || DIM_MAX2 < dim_num )
14310 cout <<
"I8_SOBOL - Fatal error!\n";
14311 cout <<
" The spatial dimension DIM_NUM should satisfy:\n";
14312 cout <<
" 1 <= DIM_NUM <= " << DIM_MAX2 <<
"\n";
14313 cout <<
" But this input value is DIM_NUM = " << dim_num <<
"\n";
14317 m_dim_num_save = dim_num;
14333 for ( j = 0; j < m_maxcol; j++ )
14340 for ( i = 1; i < dim_num; i++ )
14365 for ( k = m-1; 0 <= k; k-- )
14368 includ[k] = ( j != ( 2 * j2 ) );
14377 for ( j = m; j < m_maxcol; j++ )
14382 for ( k = 0; k < m; k++ )
14388 newv = ( newv ^ ( l * v[i][j-k-1] ) );
14398 for ( j = m_maxcol - 2; 0 <= j; j-- )
14401 for ( i = 0; i < dim_num; i++ )
14403 v[i][j] = v[i][j] * l;
14409 recipd = 1.0E+00 / ( ( double ) ( 2 * l ) );
14416 for ( i = 0; i < dim_num; i++ )
14421 else if ( *seed == m_seed_save + 1 )
14423 l = i8_bit_lo0 ( *seed );
14425 else if ( *seed <= m_seed_save )
14429 for ( i = 0; i < dim_num; i++ )
14434 for ( seed_temp = m_seed_save; seed_temp <= (*seed)-1; seed_temp++ )
14437 l = i8_bit_lo0 ( seed_temp );
14439 for ( i = 0; i < dim_num; i++ )
14441 lastq[i] = ( lastq[i] ^ v[i][l-1] );
14444 l = i8_bit_lo0 ( *seed );
14446 else if ( m_seed_save+1 < *seed )
14448 for ( seed_temp = m_seed_save+1; seed_temp <= (*seed)-1; seed_temp++ )
14451 l = i8_bit_lo0 ( seed_temp );
14453 for ( i = 0; i < dim_num; i++ )
14455 lastq[i] = ( lastq[i] ^ v[i][l-1] );
14458 l = i8_bit_lo0 ( *seed );
14463 if ( m_maxcol < l )
14466 cout <<
"I8_SOBOL - Fatal error!\n";
14467 cout <<
" The value of SEED seems to be too large (too many sample points requested!)\n";
14468 cout <<
" SEED = " << *seed <<
"\n";
14469 cout <<
" MAXCOL = " << m_maxcol <<
"\n";
14470 cout <<
" L = " << l <<
"\n";
14477 for ( i = 0; i < dim_num; i++ )
14479 quasi[i] = ( ( double ) lastq[i] ) * recipd;
14481 lastq[i] = ( lastq[i] ^ v[i][l-1] );
14484 m_seed_save = *seed;
14495 unsigned int *lhs::perm_uniform (
unsigned int n)
14534 p =
new unsigned int[n];
14536 for ( i = 0; i < n; i++ )
14543 for ( i = 0; i < n; i++ )
14546 j = i+ (
unsigned int ) ( r * (
double) (n-i));
14556 std::vector<double> lhs::latin_random (
unsigned int dim_num,
unsigned int point_num)
14596 unsigned int *perm;
14598 std::vector<double> x(dim_num*point_num,0.0);
14607 for ( i = 0; i < dim_num; i++ )
14609 perm= perm_uniform ( point_num );
14611 for ( j = 0; j < point_num; j++ )
14614 x[k] = ( ( ( double ) perm[j] ) + r ) / ( (
double ) point_num );
14631 halton::halton(
unsigned int dim,
unsigned int count) :
base(dim,count), m_primes() {
14632 if (dim >10 || dim==0) {
14633 pagmo_throw(value_error,
"Halton sequences should not be used in dimension >10");
14636 pagmo_throw(value_error,
"The first element of the sequence has index 1");
14638 for (
size_t i=1; i<=dim; ++i) {
14639 m_primes.push_back(
prime(i));
14655 std::vector<double> retval;
14656 for (
size_t i=0; i<
m_dim; ++i) {
14671 pagmo_throw(value_error,
"Halton sequence first point id is 1");
14673 std::vector<double> retval;
14674 for (
size_t i=0; i<m_primes.size(); ++i) {
14690 faure::faure(
unsigned int dim,
unsigned int count) :
base(dim, count), m_coef(NULL), m_hisum_save(-1), m_qs(-1), m_ytemp(NULL) {
14691 if (dim >23 || dim <2) {
14692 pagmo_throw(value_error,
"Faure sequences can have dimension [2,23]");
14707 std::vector<double> retval(
m_dim,0.0);
14720 std::vector<double> retval(
m_dim,0.0);
14733 if (dim >10 || dim==0) {
14734 pagmo_throw(value_error,
"Halton sequences should not be used in dimension >10");
14737 pagmo_throw(value_error,
"The first element of the sequence has index 1");
14753 std::vector<double> tmp = m_generator();
14755 std::vector<double> retval = m_projector(tmp);
14766 std::vector<double> retval = m_projector(m_generator(n));
14780 sobol::sobol(
unsigned int dim,
unsigned int count) :
base(dim, count), m_dim_num_save(0), m_initialized(false), m_maxcol(62), m_seed_save(-1), recipd(0), lastq(), poly(), v(){
14781 if (dim >1111 || dim <1) {
14782 pagmo_throw(value_error,
"This Sobol sequence can have dimensions [1,1111]");
14797 std::vector<double> retval(
m_dim,0.0);
14798 signed long long int seed= (
signed long long int)
m_count;
14799 i8_sobol(
m_dim, &seed, &retval[0]);
14800 m_count= (
unsigned short int) seed;
14812 signed long long int seed= (
signed long long int)
m_count;
14813 std::vector<double> retval(
m_dim,0.0);
14814 i8_sobol(
m_dim, &seed, &retval[0]);
14815 m_count= (
unsigned short int) seed;
14825 lhs::lhs(
unsigned int dim,
unsigned int count) :
base(dim, count), m_initialised(false), m_set(), m_next(0) {
14839 std::vector<double> retval(
m_dim,0.0);
14840 if (!m_initialised){
14842 m_initialised=
true;
14845 for (
size_t i=0;i<
m_dim;i++){
14846 retval[i]=m_set[i+m_next*
m_dim];
14859 std::vector<double> retval(
m_dim,0.0);
14861 if (!m_initialised){
14863 m_initialised=
true;
14865 for (
size_t i=0;i<
m_dim;i++){
14866 retval[i]=m_set[i+m_next*
m_dim];
std::vector< double > operator()()
Operator ()
boost::shared_ptr< base > base_ptr
Smart pointer to the base discrepancy class.
base_ptr clone() const
Clone method.
std::vector< double > operator()()
Operator ()
unsigned int prime(int n)
Returns a prime number.
std::vector< double > operator()()
Operator ()
faure(unsigned int dim, unsigned int count=1)
Constructor.
double van_der_corput(unsigned int n, unsigned int base)
Van Der Corput sequence.
base_ptr clone() const
Clone method.
unsigned int prime_ge(unsigned int n)
Returns the smallest prime greater than or equal to n.
base_ptr clone() const
Clone method.
unsigned int m_count
Starting point of the sequence (can be used to skip initial values)
simplex(unsigned int dim, unsigned int count)
Constructor.
halton(unsigned int dim, unsigned int count=1)
Constructor.
std::vector< double > operator()()
Operator ()
base_ptr clone() const
Clone method.
sobol(unsigned int dim, unsigned int count)
Constructor.
base_ptr clone() const
Clone method.
Base low-discrepancy sequence class.
unsigned int m_dim
Hypercube dimension where sampling with low-discrepancy.
std::vector< double > operator()()
Operator ()
lhs(unsigned int dim, unsigned int count)
Constructor.