github.com/gopherd/gonum@v0.0.4/stat/samplemv/halton.go (about)

     1  // Copyright ©2017 The Gonum Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package samplemv
     6  
     7  import (
     8  	"fmt"
     9  
    10  	"math/rand"
    11  
    12  	"github.com/gopherd/gonum/mat"
    13  	"github.com/gopherd/gonum/stat/distmv"
    14  )
    15  
    16  // Halton is a type for sampling using the Halton sequence from
    17  // the given distribution. The specific method for scrambling (or lack thereof)
    18  // is specified by the HaltonKind. If src is not nil, it will be used to generate
    19  // the randomness needed to scramble the sequence (if necessary), otherwise
    20  // the rand package will be used. Halton panics if the HaltonKind is unrecognized
    21  // or if q is nil.
    22  //
    23  // Halton sequence random number generation is a quasi-Monte Carlo procedure
    24  // where the samples are generated to be evenly spaced out across the distribution.
    25  // Note that this means the sample locations are correlated with one another.
    26  // The distmv.NewUnitUniform function can be used for easy sampling from the unit hypercube.
    27  type Halton struct {
    28  	Kind HaltonKind
    29  	Q    distmv.Quantiler
    30  	Src  rand.Source
    31  }
    32  
    33  // Sample generates rows(batch) samples using the Halton generation procedure.
    34  func (h Halton) Sample(batch *mat.Dense) {
    35  	halton(batch, h.Kind, h.Q, h.Src)
    36  }
    37  
    38  // HaltonKind specifies the type of algorithm used to generate Halton samples.
    39  type HaltonKind int
    40  
    41  const (
    42  	// Owen generates (scrambled) Halton samples using the Randomized van der Corput
    43  	// algorithm described in
    44  	//  A randomized Halton algorithm
    45  	//  Art Owen
    46  	//  https://arxiv.org/pdf/1706.02808.pdf
    47  	// Currently limited to 1000 dimensional inputs.
    48  	Owen = iota + 1
    49  )
    50  
    51  func halton(batch *mat.Dense, kind HaltonKind, q distmv.Quantiler, src rand.Source) {
    52  	// Code based from https://arxiv.org/pdf/1706.02808.pdf .
    53  	perm := rand.Perm
    54  	if src != nil {
    55  		perm = rand.New(src).Perm
    56  	}
    57  
    58  	n, d := batch.Dims()
    59  	// Each case should generate random numbers over the unit cube.
    60  	switch kind {
    61  	default:
    62  		panic("halton: unknown HaltonKind")
    63  	case Owen:
    64  		for j := 0; j < d; j++ {
    65  			b := nthPrime(j)
    66  			div := int64(1)
    67  			b2r := 1 / float64(b)
    68  			for 1-b2r < 1 {
    69  				p := perm(b)
    70  				for i := 0; i < n; i++ {
    71  					dig := (int64(i) / div) % int64(b)
    72  					pdig := float64(p[dig])
    73  					v := batch.At(i, j)
    74  					v += pdig * b2r
    75  					batch.Set(i, j, v)
    76  				}
    77  				div *= int64(b)
    78  				b2r /= float64(b)
    79  			}
    80  		}
    81  	}
    82  	p := make([]float64, d)
    83  	for i := 0; i < n; i++ {
    84  		copy(p, batch.RawRowView(i))
    85  		q.Quantile(batch.RawRowView(i), p)
    86  	}
    87  }
    88  
    89  // nthPrime returns the nth prime number (0 indexed).
    90  func nthPrime(n int) int {
    91  	if n > len(firstPrimes) {
    92  		panic(fmt.Sprintf("halton: dimension must be less than %d", len(firstPrimes)))
    93  	}
    94  	return firstPrimes[n]
    95  }
    96  
    97  var firstPrimes = []int{2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47,
    98  	53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131,
    99  	137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211,
   100  	223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293,
   101  	307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389,
   102  	397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479,
   103  	487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587,
   104  	593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673,
   105  	677, 683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773,
   106  	787, 797, 809, 811, 821, 823, 827, 829, 839, 853, 857, 859, 863, 877, 881,
   107  	883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991,
   108  	997, 1009, 1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069,
   109  	1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123, 1129, 1151, 1153, 1163, 1171,
   110  	1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277,
   111  	1279, 1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367,
   112  	1373, 1381, 1399, 1409, 1423, 1427, 1429, 1433, 1439, 1447, 1451, 1453, 1459,
   113  	1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553,
   114  	1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627,
   115  	1637, 1657, 1663, 1667, 1669, 1693, 1697, 1699, 1709, 1721, 1723, 1733, 1741,
   116  	1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847, 1861,
   117  	1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951,
   118  	1973, 1979, 1987, 1993, 1997, 1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053,
   119  	2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131, 2137, 2141,
   120  	2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267,
   121  	2269, 2273, 2281, 2287, 2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351,
   122  	2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417, 2423, 2437, 2441,
   123  	2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557,
   124  	2579, 2591, 2593, 2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677,
   125  	2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719, 2729, 2731, 2741, 2749,
   126  	2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851,
   127  	2857, 2861, 2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963,
   128  	2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037, 3041, 3049, 3061, 3067, 3079,
   129  	3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203,
   130  	3209, 3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313,
   131  	3319, 3323, 3329, 3331, 3343, 3347, 3359, 3361, 3371, 3373, 3389, 3391, 3407,
   132  	3413, 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527,
   133  	3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613,
   134  	3617, 3623, 3631, 3637, 3643, 3659, 3671, 3673, 3677, 3691, 3697, 3701, 3709,
   135  	3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821, 3823,
   136  	3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923,
   137  	3929, 3931, 3943, 3947, 3967, 3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027,
   138  	4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129, 4133, 4139,
   139  	4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243, 4253,
   140  	4259, 4261, 4271, 4273, 4283, 4289, 4297, 4327, 4337, 4339, 4349, 4357, 4363,
   141  	4373, 4391, 4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457, 4463, 4481, 4483,
   142  	4493, 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597,
   143  	4603, 4621, 4637, 4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703,
   144  	4721, 4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789, 4793, 4799, 4801, 4813,
   145  	4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943,
   146  	4951, 4957, 4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023,
   147  	5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101, 5107, 5113, 5119, 5147, 5153,
   148  	5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279,
   149  	5281, 5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, 5407,
   150  	5413, 5417, 5419, 5431, 5437, 5441, 5443, 5449, 5471, 5477, 5479, 5483, 5501,
   151  	5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623,
   152  	5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711,
   153  	5717, 5737, 5741, 5743, 5749, 5779, 5783, 5791, 5801, 5807, 5813, 5821, 5827,
   154  	5839, 5843, 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903, 5923,
   155  	5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067,
   156  	6073, 6079, 6089, 6091, 6101, 6113, 6121, 6131, 6133, 6143, 6151, 6163, 6173,
   157  	6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269, 6271, 6277,
   158  	6287, 6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367,
   159  	6373, 6379, 6389, 6397, 6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521,
   160  	6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599, 6607, 6619, 6637,
   161  	6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737,
   162  	6761, 6763, 6779, 6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, 6841, 6857,
   163  	6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 6947, 6949, 6959, 6961, 6967,
   164  	6971, 6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069,
   165  	7079, 7103, 7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, 7211,
   166  	7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283, 7297, 7307, 7309, 7321, 7331,
   167  	7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457, 7459, 7477, 7481,
   168  	7487, 7489, 7499, 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561,
   169  	7573, 7577, 7583, 7589, 7591, 7603, 7607, 7621, 7639, 7643, 7649, 7669, 7673,
   170  	7681, 7687, 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789,
   171  	7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907,
   172  	7919,
   173  }