Added SOBOL

Author: blankjul

.idea/pysample.iml

@@ -14,41 +14,44 @@ def correlation(X):
     return np.sum(np.tril(M, -1) ** 2)
 
 
-def centered_l2_discrepancy_slow(X, check_bounds=True):
-    if check_bounds and (np.any(np.min(X) < 0 or np.max(X) > 1)):
-        raise Exception("Make sure that all X are between 0 and 1.")
+def centered_l2_discrepancy(X):
 
     n_points, n_dim = X.shape
-
     acmh = np.abs(X - 0.5)
 
-    val = 0
-
-    for k in range(n_points):
-        temp = 1
-        for i in range(n_dim):
-            temp = temp * (1 + 0.5 * (acmh[:, i] + acmh[k, i] - np.abs(X[:, i] - X[k, i])))
-        val += np.sum(temp)
+    _X = np.abs(X.T[..., None] - X.T[:, None, :])
+    _acmh = np.abs(acmh.T[..., None] + acmh.T[:, None, :])
 
+    val = np.sum(np.prod((1 + 0.5 * (_acmh - _X)), axis=0))
     val = ((13 / 12) ** n_dim - 2 / n_points * np.sum(
         np.prod(1 + 0.5 * (acmh - acmh ** 2), axis=1)) + val / n_points ** 2) ** 0.5
 
     return val
 
 
-def centered_l2_discrepancy(X, check_bounds=True):
-    if check_bounds and (np.any(np.min(X) < 0 or np.max(X) > 1)):
-        raise Exception("Make sure that all X are between 0 and 1.")
 
-    n_points, n_dim = X.shape
+def wrapped_l2_discrepancy(X):
 
-    acmh = np.abs(X - 0.5)
+    """
 
-    _X = np.abs(X.T[..., None] - X.T[:, None, :])
-    _acmh = np.abs(acmh.T[..., None] + acmh.T[:, None, :])
+    [Npts,Ndim]=size(coord);
+    WDL2=0;
+    for k=1:Npts
+        yada=abs(coord(:,1)-coord(k,1));
+        temp=(1.5-yada.*(1-yada));
+        for i=2:Ndim
+            yada=abs(coord(:,i)-coord(k,i));
+            temp=temp.*(1.5-yada.*(1-yada));
+        end
+        WDL2=WDL2+sum(temp);
+    end
+    WDL2=(-(4/3)^Ndim+WDL2/Npts^2)^0.5; %^0.5;
 
-    val = np.sum(np.prod((1 + 0.5 * (_acmh - _X)), axis=0))
-    val = ((13 / 12) ** n_dim - 2 / n_points * np.sum(
-        np.prod(1 + 0.5 * (acmh - acmh ** 2), axis=1)) + val / n_points ** 2) ** 0.5
+    end
+    """
+    n_points, n_dim = X.shape
+    _X = np.abs(X.T[..., None] - X.T[:, None, :])
+    val = np.sum(np.prod((1.5 - _X * (1-_X)), axis=0))
+    val = (-(4 / 3) ** n_dim + val / n_points ** 2) ** 0.5
 
     return val

.idea/workspace.xml

@@ -0,0 +1,24 @@
+import numpy as np
+
+from pysample.methods.sampling import Sampling
+from pysample.util import calc_primes_until
+
+
+class HaltonSampling(Sampling):
+
+    def _sample(self, n_samples, n_dim):
+        bases = calc_primes_until(500)[:n_dim]
+        X = np.column_stack([self.halton_sequence(n_samples, b) for b in bases])
+        return X
+
+    def halton_sequence(self, n, b):
+        return np.array([self.halton_sequence_by_index(i, b) for i in range(n)])
+
+    def halton_sequence_by_index(self, i, b):
+        f = 1.0
+        x = 0.0
+        while i > 0:
+            f /= b
+            x += f * (i % b)
+            i = np.floor(i / b)
+        return x

pysample/indicators.py

@@ -1,7 +1,7 @@
 import numpy as np
 
 from pysample.methods.sampling import Sampling
-from pysample.indicators import minimum_distance, correlation
+from pysample.indicators import minimum_distance, correlation, centered_l2_discrepancy, wrapped_l2_discrepancy
 from pysample.optimization import random_search
 
 
@@ -17,12 +17,16 @@ def __init__(self, criterion="maxmin", optimizer="random_search", iterations=100
             self.fun_obj = minimum_distance
         elif criterion == "correlation":
             self.fun_obj = correlation
+        elif criterion == "centered_l2_discrepancy":
+            self.fun_obj = centered_l2_discrepancy
+        elif criterion == "wrapped_l2_discrepancy":
+            self.fun_obj = wrapped_l2_discrepancy
         else:
             raise Exception("Criterion not known.")
 
     def _sample(self, n_samples, n_dim):
 
-        if self.optimizer == "random_trials":
+        if self.optimizer == "random_search":
             fun_sample = lambda: _sample(n_samples, n_dim)
             X = random_search(fun_sample, self.fun_obj, self.iterations)
 
@@ -31,13 +35,10 @@ def _sample(self, n_samples, n_dim):
             pass
 
         elif self.optimizer == "ga":
-
             n_generations = 100
             pop_size = 20
 
 
-
-
         else:
             X = _sample(n_samples, n_dim)
 

pysample/methods/halton.py

@@ -0,0 +1,94 @@
+import numpy as np
+
+import re
+
+from pysample.methods.sampling import Sampling
+from pysample.util import path_to_resources
+
+
+class SobolSampling(Sampling):
+
+    def __init__(self, n_skip=-1, n_leap=0, setup="matlab") -> None:
+        super().__init__()
+        fname = "sobol_%s.dat" % setup
+        self.setup = parse_file(path_to_resources(fname))
+
+        self.n_skip = n_skip
+        self.n_leap = n_leap
+
+    def _sample(self, n_samples, n_dim):
+
+        # find out how long the sequence needs to be - skip or leap included
+        I = np.arange(0, n_samples, self.n_leap+1)
+        if self.n_skip == -1:
+            I += n_samples
+        else:
+            I += self.n_skip
+        n_sequence = np.max(I) + 1
+
+        # containts all the values of the sequences as integer
+        _X = np.zeros((n_sequence, n_dim), dtype=np.int)
+
+        # number of bits which will be necessary for this sequence
+        L = int(np.ceil(np.log2(n_sequence)))
+
+        # number of bits necessary for the equation for each of them
+        C = [highest_bit(i) for i in range(0, n_sequence)]
+
+        for j in range(n_dim):
+
+            if j == 0:
+                V = np.concatenate((np.array([0]), 2 ** np.arange(31, -1, -1)))
+
+            else:
+
+                s, a, m = self.setup[j]["s"], self.setup[j]["a"], self.setup[j]["m"]
+                V = np.zeros(L + 1, dtype=np.int)
+
+                if L <= s:
+                    for k in range(1, L + 1):
+                        V[k] = m[k - 1] << (32 - k)
+
+                else:
+
+                    for k in range(1, s + 1):
+                        V[k] = m[k - 1] << (32 - k)
+
+                    for i in range(s + 1, L + 1):
+                        V[i] = V[i - s] ^ int(V[i - s] >> s)
+                        for k in range(1, s):
+                            V[i] ^= (((a >> (s - 1 - k)) & 1) * V[i - k])
+
+            for i in range(1, n_sequence):
+                _X[i, j] = _X[i - 1, j] ^ V[C[i - 1]]
+
+        X = (_X / 2 ** 32)[I]
+        return X
+
+
+def highest_bit(i):
+    bit = 1
+    while i > 0 and i % 2 != 0:
+        i = np.floor(i / 2)
+        bit += 1
+    return bit
+
+
+def parse_file(path_to_file):
+    setup = [{}]
+
+    with open(path_to_file) as f:
+        content = f.read().splitlines()[1:]
+
+    for line in content:
+        entries = [int(e) for e in re.split("\s+|\t", line.strip())]
+
+        setup.append({
+            'd': entries[0],
+            's': entries[1],
+            'a': entries[2],
+            'm': np.array(entries[3:])
+
+        })
+
+    return setup

pysample/methods/lhs.py

@@ -12,6 +12,8 @@ def random_search(fun_sample, fun_obj, iterations):
             f = _f
             X = _X
 
+            print(f)
+
     return X
 
 

pysample/methods/sobol.py

@@ -0,0 +1,41 @@
+d	s	a	m_i     
+2   1  1  1
+3   1  3  1
+4   2  7  1 1
+5   3  11  1 3 7
+6   3  13  1 1 5
+7   4  19  1 3 1 1
+8   4  25  1 1 3 7
+9   5  37  1 3 3 9 9
+10   5  59  1 3 7 13 3
+11   5  47  1 1 5 11 27
+12   5  61  1 3 5 1 15
+13   5  55  1 1 7 3 29
+14   5  41  1 3 7 7 21
+15   6  67  1 1 1 9 23 37
+16   6  97  1 3 3 5 19 33
+17   6  91  1 1 3 13 11 7
+18   6  109  1 1 7 13 25 5
+19   6  103  1 3 5 11 7 11
+20   6  115  1 1 1 3 13 39
+21   7  131  1 3 1 15 17 63 13
+22   7  193  1 1 5 5 1 27 33
+23   7  137  1 3 3 3 25 17 115
+24   7  145  1 1 3 15 29 15 41
+25   7  143  1 3 1 7 3 23 79
+26   7  241  1 3 7 9 31 29 17
+27   7  157  1 1 5 13 11 3 29
+28   7  185  1 3 1 9 5 21 119
+29   7  167  1 1 3 1 23 13 75
+30   7  229  1 3 3 11 27 31 73
+31   7  171  1 1 7 7 19 25 105
+32   7  213  1 3 5 5 21 9 7
+33   7  191  1 1 1 15 5 49 59
+34   7  253  1 1 1 1 1 33 65
+35   7  203  1 3 5 15 17 19 21
+36   7  211  1 1 7 11 13 29 3
+37   7  239  1 3 7 5 7 11 113
+38   7  247  1 1 5 3 15 19 61
+39   8  285  1 3 1 1 9 27 89 7
+40   8  369  1 1 3 7 31 15 45 23
+41   8  299  1 3 3 9 9 25 107 39