Lattice kronecker

docs/api/discrete_distributions.md

@@ -28,6 +28,10 @@ jupyter:
 
 ::: qmcpy.discrete_distribution.korobov.KorobovLattice
 
+## `lattice_vector_wssd_search`
+
+::: qmcpy.discrete_distribution.lattice.lattice_vector_wssd_search.lattice_vector_wssd_search
+
 ## `Halton`
 
 ::: qmcpy.discrete_distribution.digital_net_any_bases.halton.Halton
@@ -52,6 +56,10 @@ jupyter:
 
 ::: qmcpy.discrete_distribution.latin_hypercube.LatinHypercube
 
+## `kronecker_vector_search_mobius_transform`
+
+::: qmcpy.discrete_distribution.kronecker.kronecker_search_methods.kronecker_vector_search_mobius_transform
+
 ## `DummySampler`
 
 ::: qmcpy.discrete_distribution.dummy_sampler.DummySampler

pyproject.toml

@@ -179,6 +179,7 @@ includes = [
     "qmcpy/discrete_distribution/generating_params/korobov_p2_table.npz",
     "qmcpy/discrete_distribution/digital_net_b2/generating_matrices/*.npy",
     "qmcpy/discrete_distribution/lattice/generating_vectors/*.npy",
+    "qmcpy/discrete_distribution/kronecker/generating_vectors/*.txt",
     "qmcpy/util/qmcpy.mplstyle",
     ]
 excludes = []

qmcpy/discrete_distribution/__init__.py

@@ -1,10 +1,10 @@
 from .abstract_discrete_distribution import AbstractDiscreteDistribution
 from .iid_std_uniform import IIDStdUniform
-from .lattice import Lattice
+from .lattice import Lattice, lattice_vector_wssd_search
 from .digital_net_b2 import DigitalNetB2
 from .digital_net_any_bases import DigitalNetAnyBases,Halton,Faure,Hammersley
 from .mpmc import MPMC
-from .kronecker import Kronecker
+from .kronecker import Kronecker, kronecker_vector_search_mobius_transform
 from .korobov import KorobovLattice
 from .dummy_sampler import DummySampler
 from .latin_hypercube import LatinHypercube
@@ -15,4 +15,3 @@
 DigitalNet = DigitalNetB2
 Net = DigitalNetB2
 NetB2 = DigitalNetB2
-

qmcpy/discrete_distribution/kronecker/__init__.py

@@ -0,0 +1,2 @@
+from .kronecker import Kronecker
+from .kronecker_search_methods import kronecker_vector_search_mobius_transform

qmcpy/discrete_distribution/kronecker/generating_vectors/kron_vector_d-100_N-2exp20_2026_06_01.txt

@@ -0,0 +1,100 @@
+0.618033988749895
+0.3173225474723
+0.59332263014446
+0.20776441643926
+0.27373719258623
+0.649734278361753
+0.478954018631769
+0.86866022435182
+0.22845082022244
+0.581365429377986
+0.282365231829842
+0.0822850909119904
+0.223849641007295
+0.5770772201756
+0.51769659336634
+0.568025390904592
+0.156782234569368
+0.82246227056154
+0.805675312097409
+0.63877102813393
+0.358300563495856
+0.241741343018598
+0.705003192174204
+0.1931911954956
+0.261022001488623
+0.897938992038015
+0.46839743115877
+0.884022067965329
+0.752352896871505
+0.1601583600427
+0.10727599509739
+0.151478435512877
+0.163863657127101
+0.948303450359399
+0.80350943597439
+0.426371623468333
+0.435930910910882
+0.21329852459791
+0.661698149534002
+0.900679822160453
+0.122436710671457
+0.483663584095611
+0.928181067731583
+0.443143014606576
+0.74491332336194
+0.87948409225588
+0.0428242449803
+0.534576896789579
+0.24340042100879
+0.30424418245585
+0.574003104342617
+0.897289023268963
+0.541424476559586
+0.356895660350464
+0.507567280910795
+0.513983550428507
+0.0610821922457415
+0.183871471606587
+0.446015178033969
+0.455684287415085
+0.280817534817491
+0.115220095666085
+0.433740673279323
+0.515605957977756
+0.113076735656464
+0.733928297688305
+0.0597515651584137
+0.422268695684775
+0.0979181139173599
+0.213699261322352
+0.866811679881922
+0.0878569329036737
+0.678412735893121
+0.181093969536107
+0.128913741473518
+0.109341703717108
+0.289067270578427
+0.352218331663839
+0.303605902333137
+0.0613899204730832
+0.959535877660851
+0.475508309069064
+0.688698902674194
+0.657037932118495
+0.645555897563869
+0.720658665263604
+0.914423387894897
+0.425763295044487
+0.328825255006553
+0.892452975558004
+0.16973367306396
+0.912292406867098
+0.0923260018966512
+0.216301713289429
+0.147861410064151
+0.8600781655845
+0.752129792595509
+0.337431120990153
+0.542476014178907
+0.307279789725491

qmcpy/discrete_distribution/kronecker.py → qmcpy/discrete_distribution/kronecker/kronecker.py

@@ -1,5 +1,5 @@
-from .abstract_discrete_distribution import AbstractLDDiscreteDistribution
-from ..util import ParameterError
+from ..abstract_discrete_distribution import AbstractLDDiscreteDistribution
+from ...util import ParameterError
 import numpy as np
 import warnings
 
@@ -243,6 +243,7 @@ def __init__(self,
                 - `"CBC"`: uses the first $d$ components of a known good Component-by-Component (CBC) generating vector.
                 - `"RICHTMYER"`: uses $\boldsymbol{\alpha}_j = \sqrt{p_j} \bmod 1$, where $p_j$ are primes. This is the classical Richtmyer construction.
                 - `"SUZUKI"`: uses a deterministic construction $\boldsymbol{\alpha}_j = 2^{j/(d+1)}$.
+                - `"CBC_MT"`: uses the first $d$ components of a known good CBC generating vector obtained using the Mobius transformation method, which can be found in kronecker_search_methods.py.
                 - np.array: user-specified generating vector.
 
             shift (np.ndarray): Shift vector $\boldsymbol{\delta}$. If `randomize=True`, this is ignored and a random shift is generated. Otherwise, a fixed shift is used.
@@ -284,7 +285,118 @@ def __init__(self,
             gen_vec = _richtmyer_generating_vector(self.dvec.max()+1)
         elif isinstance(generating_vector, str) and generating_vector.lower() == "suzuki":
             self.gen_vec_source = "SUZUKI"
-            gen_vec = _suzuki_generating_vector(self.dvec.max()+1)
+            gen_vec = _suzuki_generating_vector(self.dvec.max()+1)       
+        elif isinstance(generating_vector, str) and generating_vector.lower() == "cbc_mt":
+            self.gen_vec_source = "CBC_MT"
+            CBC_MT = np.array([0.618033988749895,
+                0.3173225474723,
+                0.59332263014446,
+                0.20776441643926,
+                0.27373719258623,
+                0.649734278361753,
+                0.478954018631769,
+                0.86866022435182,
+                0.22845082022244,
+                0.581365429377986,
+                0.282365231829842,
+                0.0822850909119904,
+                0.223849641007295,
+                0.5770772201756,
+                0.51769659336634,
+                0.568025390904592,
+                0.156782234569368,
+                0.82246227056154,
+                0.805675312097409,
+                0.63877102813393,
+                0.358300563495856,
+                0.241741343018598,
+                0.705003192174204,
+                0.1931911954956,
+                0.261022001488623,
+                0.897938992038015,
+                0.46839743115877,
+                0.884022067965329,
+                0.752352896871505,
+                0.1601583600427,
+                0.10727599509739,
+                0.151478435512877,
+                0.163863657127101,
+                0.948303450359399,
+                0.80350943597439,
+                0.426371623468333,
+                0.435930910910882,
+                0.21329852459791,
+                0.661698149534002,
+                0.900679822160453,
+                0.122436710671457,
+                0.483663584095611,
+                0.928181067731583,
+                0.443143014606576,
+                0.74491332336194,
+                0.87948409225588,
+                0.0428242449803,
+                0.534576896789579,
+                0.24340042100879,
+                0.30424418245585,
+                0.574003104342617,
+                0.897289023268963,
+                0.541424476559586,
+                0.356895660350464,
+                0.507567280910795,
+                0.513983550428507,
+                0.0610821922457415,
+                0.183871471606587,
+                0.446015178033969,
+                0.455684287415085,
+                0.280817534817491,
+                0.115220095666085,
+                0.433740673279323,
+                0.515605957977756,
+                0.113076735656464,
+                0.733928297688305,
+                0.0597515651584137,
+                0.422268695684775,
+                0.0979181139173599,
+                0.213699261322352,
+                0.866811679881922,
+                0.0878569329036737,
+                0.678412735893121,
+                0.181093969536107,
+                0.128913741473518,
+                0.109341703717108,
+                0.289067270578427,
+                0.352218331663839,
+                0.303605902333137,
+                0.0613899204730832,
+                0.959535877660851,
+                0.475508309069064,
+                0.688698902674194,
+                0.657037932118495,
+                0.645555897563869,
+                0.720658665263604,
+                0.914423387894897,
+                0.425763295044487,
+                0.328825255006553,
+                0.892452975558004,
+                0.16973367306396,
+                0.912292406867098,
+                0.0923260018966512,
+                0.216301713289429,
+                0.147861410064151,
+                0.8600781655845,
+                0.752129792595509,
+                0.337431120990153,
+                0.542476014178907,
+                0.307279789725491], dtype=np.float64)
+            gen_vec = CBC_MT
+            if not (self.dvec.max() < len(gen_vec)):
+                if warn:
+                    warnings.warn(
+                        f"CBC_MT generating vector only supports dimension <= {len(CBC_MT)}; falling back to Richtmyer.",
+                        RuntimeWarning,
+                    )
+                self.gen_vec_source = "RICHTMYER"
+                gen_vec = _richtmyer_generating_vector(self.dvec.max()+1)
         else:
             self.gen_vec_source = "CUSTOM"
             gen_vec = np.asarray(generating_vector, dtype=float)
@@ -347,9 +459,9 @@ def periodic_discrepancy(self, n, k_tilde=None, gamma=None):
             k_tilde = (lambda x, gamma: np.prod(1 + (x * (x - 1) + 1/6) * gamma, axis=-1), 1)
 
         return np.sqrt(self._square_periodic_discrepancies(n, k_tilde, gamma))
+
 
-
-    def wssd_discrepancy(self, n, weights, k_tilde = None, gamma = None):
+    def wssd_discrepancy(self, n, sample_weights, k_tilde = None, gamma = None):
         # calculates the weighted sum of square discrepancy
         if gamma is None:
             gamma = np.ones(self.d)
@@ -358,12 +470,15 @@ def wssd_discrepancy(self, n, weights, k_tilde = None, gamma = None):
             k_tilde = (lambda x, gamma: np.prod(1 + (x * (x - 1) + 1/6) * gamma, axis=-1), 1)
 
         discrepancies = self._square_periodic_discrepancies(n, k_tilde, gamma)
-        return np.sum(weights * discrepancies, axis=-1)
-
+        return np.sum(sample_weights * discrepancies, axis=-1)
 
+
     def _square_periodic_discrepancies(self, n, k_tilde, gamma):
         n_array = np.arange(1, n + 1)
-        k_tilde_terms = k_tilde[0](self.gen_samples(n=n), gamma)
+        # we need the points without a random shift for the calculation, so we can't use self._gen_samples
+        i = np.arange(0, n)
+        points = (i[:,None] * self.gen_vec[:,None,:]) % 1
+        k_tilde_terms = k_tilde[0](points, gamma)
 
         left_sum = np.cumsum(k_tilde_terms[...,1:], axis=-1) * n_array[1:]
         right_sum = np.cumsum(n_array[:-1] * k_tilde_terms[...,1:], axis=-1)