8 files changed, 144 insertions(+), 37 deletions(-)

A clusterfunc/plot.py
M clusterfunc/single.py
M clusterfunc/symmetric.py
M clusterfunc/util.py
A matplotlibrc
M tests/not_test_unit.py
M tests/test_single.py
A tests/test_symmetric.py

@@ 0,0 1,35 @@
+import matplotlib.pyplot as plt
+import numpy as np
+
+def plot(x):
+    n = np.linspace(0, len(x)-1, len(x))
+    fig, ax = plt.subplots()
+    ax.plot(n, x, 'k-x')
+    ax.plot(np.zeros_like(x), x, 'k-x')
+    ax.set_xlim(n[0]-1, n[-1])
+    ax.set_ylim(x.min(), x.max())
+    ax.set_xlabel('Index')
+    ax.set_ylabel('Coordinate')
+
+def plot_ER(x):
+    ER = clusterfunc.util.ER(x)
+    n = np.linspace(0, len(ER)-1, len(ER))
+    fig, ax = plt.subplots()
+    ax.plot(n, ER, 'k-x')
+    # ax.plot(np.zeros_like(x), x, 'k-x')
+    ax.set_xlim(n[0]-1, n[-1])
+    # ax.set_ylim(x.min(), x.max())
+    ax.set_xlabel('Index')
+    ax.set_ylabel('Coordinate')
+
+import clusterfunc.single
+import clusterfunc.symmetric
+
+Dmin = 1e-1
+Dmax = 2e-1
+ERmax = 1.2
+# x = clusterfunc.single.unit_free_N(Dmin, Dmax, ERmax, mult=1)
+x = clusterfunc.symmetric.unit_free_N(Dmin, Dmax, ERmax, 1)
+plot(x)
+# plot_ER(x)
+plt.show()

@@ 3,8 3,9 @@ import numpy as np
 import clusterfunc.util
 from clusterfunc.exceptions import ClusteringException
 from scipy.optimize import root_scalar
+import clusterfunc.check
 
-def unit_fixed_N(Dmin, Dmax, ERmax, N):
+def unit_fixed_N(Dmin, Dmax, ERmax, N, check=True):
     """Single-sided clustering on the unit interval with specified number of points.
 
     Parameters


@@ 132,6 133,10 @@ def unit_fixed_N(Dmin, Dmax, ERmax, N):
     if not len(x) == N:
         raise ClusteringException(f'Incorrect number of points len(x)={len(x)}, expected {N}.')
 
+    if check:
+        clusterfunc.check.unit_single(x, Dmin, Dmax, ERmax)
+
+
     return x
 
 


@@ 140,7 145,7 @@ def unit_free_N(Dmin, Dmax, ERmax, mult=8):
 
     # Find high and low guesses for n, N = mult*n + 1
     nlow = np.round(1./Dmax/mult).astype(int)
-    nhigh = np.minimum(np.ceil(1./Dmin/mult).astype(int),1024) # Limit to protect against overflow
+    nhigh = np.minimum(np.floor(1./Dmin/mult).astype(int),1024) # Limit to protect against overflow
     if nlow > nhigh:
         nlow = 1
     if nhigh == nlow:


@@ 148,8 153,9 @@ def unit_free_N(Dmin, Dmax, ERmax, mult=8):
 
     Nhigh = nhigh * mult + 1
     try:
-        x = unit_fixed_N(Dmin, Dmax, ERmax, Nhigh)
-    except ClusteringException:
+        x = unit_fixed_N(Dmin, Dmax, ERmax, Nhigh, check=True)
+    except ClusteringException as e:
+        print(e)
         raise ClusteringException(f'Failed to cluster with high number of points guess Nhigh={Nhigh}')
 
     # Loop until converged


@@ 161,7 167,7 @@ def unit_free_N(Dmin, Dmax, ERmax, mult=8):
 
         # Attempt to cluster
         try:
-            x = unit_fixed_N(Dmin, Dmax, ERmax, Nnew)
+            x = unit_fixed_N(Dmin, Dmax, ERmax, Nnew, check=True)
             # New n is valid, use as high limiit
             nhigh = nnew
         except ClusteringException:

@@ 4,9 4,10 @@ from scipy.optimize import root_scalar
 from clusterfunc.exceptions import ClusteringException
 import clusterfunc.check
 import clusterfunc.util
+import clusterfunc.single
 
 
-def symmetric_fixed_N(Dmin, Dmax, ERmax, N):
+def unit_fixed_N(Dmin, Dmax, ERmax, N):
     """"""
     if N < 3:
         raise ClusteringException(f'Not enough points N={N} to symmetrically cluster, need N>=3')


@@ 23,8 24,8 @@ def symmetric_fixed_N(Dmin, Dmax, ERmax, N):
         Lb = (1.-xj)
 
         # Evaluate the halves
-        xa = single_fixed_N(Dmin/La, Dmax/La, ERmax, Na)*La
-        xb = 1.-np.flip(single_fixed_N(Dmin/Lb, Dmax/Lb, ERmax, Nb))*Lb
+        xa = clusterfunc.single.fixed_N(0., La, Dmin, Dmax, ERmax, Na)
+        xb = 1.-np.flip(clusterfunc.single.fixed_N(0., Lb, Dmin, Dmax, ERmax, Nb))
         x = np.concatenate((xa, xb[1:]))
 
         dxa = np.diff(xa)[-1]


@@ 60,11 61,31 @@ def symmetric_fixed_N(Dmin, Dmax, ERmax, N):
 
     return x
 
-def symmetric_free_N(Dmin, Dmax, ERmax, mult=8):
+def unit_free_N(Dmin, Dmax, ERmax, mult=8):
     """"""
-    xa = single_free_N(Dmin*2., Dmax*2., ERmax, mult)/2.
+    xa = clusterfunc.single.free_N(0., 0.5, Dmin, Dmax, ERmax, mult)
+    assert np.diff(xa).max() * (1.+1e-9) <= Dmax
     xb = 1.-np.flip(xa)
     x = np.concatenate((xa, xb[1:]))
     assert not np.mod((len(x)-1)/mult,1.)
     return x
 
+def fixed_N(x0, x1, Dmin, Dmax, ERmax, N):
+    """Double-sided symmetric clustering between two values, fixed number of points."""
+    if np.isclose(x0,x1):
+        raise ClusteringException(f'Cannot distribute points without distinct start and end points, got x0={x0} and x1={x1}')
+    dx = x1-x0
+    dxa = np.abs(dx)
+    x = x0 + dx * unit_fixed_N(Dmin/dxa, Dmax/dxa, ERmax, N)
+    return x
+
+
+def free_N(x0, x1, Dmin, Dmax, ERmax, mult=8):
+    """Double-sided symmetrirc clustering between two values, free number of points."""
+    if np.isclose(x0,x1):
+        raise ClusteringException(f'Cannot distribute points without distinct start and end points, got x0={x0} and x1={x1}')
+    dx = x1-x0
+    dxa = np.abs(dx)
+    x = x0 + dx * unit_free_N(Dmin/dxa, Dmax/dxa, ERmax, mult)
+    return x
+

@@ 9,3 9,9 @@ def split_cells(N):
 def cumsum0(x, axis=None):
     return np.insert(np.cumsum(x, axis=axis), 0, 0.0, axis=axis)
 
+def ER(x):
+    dx = np.diff(x)
+    ER = dx[1:]/dx[:-1]
+    ER[ER<1.] = 1./ER[ER<1.]
+    return ER
+

@@ 0,0 1,1 @@
+text.usetex : False

@@ 40,28 40,3 @@ def test_unif():
 
         # Check
         clusterfunc.check.unit_single(x, Dmin, Dmax, ER, rtol=1e-9)
-
-def test_symmetric():
-
-    for Dmin in  (1e-2, 1e-3, 1e-4, 1e-5):
-        for ER in (1.02, 1.05, 1.1, 1.2):
-
-            # Number of points for uniform, limit to not be too huge
-            Nmax = np.minimum(np.floor(1./Dmin).astype(int), 256)
-
-            for Dmax in (0.05, 0.1, 0.2, 1.):
-
-                if Dmax <= Dmin:
-                    continue
-
-                # Minimum number of points with capping
-                Nmin = len(clusterfunc.unit.symmetric_free_N(Dmin, Dmax, ER, mult=1))
-
-                # Test without capping
-                for N in range(Nmin,Nmax+1):
-
-                    # Evaluate
-                    x = clusterfunc.unit.symmetric_fixed_N(Dmin, Dmax, ER, N)
-
-                    # Check
-                    clusterfunc.check.unit_symmetric(x, Dmin, Dmax, ER, rtol=1e-9)

@@ 64,5 64,5 @@ def test_abs():
                         assert np.allclose(x[(0,-1),], (x0, x1))
                         dx = np.diff(x)
                         assert np.diff(np.sign(dx)).all() == 0.
-                        assert np.abs(dx.max()) <= Dmax * (1.+rtol)
-                        assert np.abs(dx.min()) * (1.+rtol) >= Dmin
+                        assert np.abs(dx).max() <= Dmax * (1.+rtol)
+                        assert np.abs(dx).min() * (1.+rtol) >= Dmin

@@ 0,0 1,63 @@
+import clusterfunc.symmetric
+import clusterfunc.check
+import numpy as np
+import pytest
+from clusterfunc.exceptions import ClusteringException
+
+
+def test_symmetric_unit():
+
+    for Dmin in  (1e-1, 1e-2, 1e-3, 1e-4, 1e-5):
+        for ER in (1.02, 1.05, 1.1, 1.2):
+
+            # Number of points for uniform, limit to not be too huge
+            Nmax = np.minimum(np.floor(1./Dmin).astype(int), 256)
+
+            for Dmax in (0.05, 0.1, 0.2, 1.):
+
+                if Dmax <= Dmin:
+                    continue
+
+                # Minimum number of points with capping
+                Nmin = len(clusterfunc.symmetric.unit_free_N(Dmin, Dmax, ER, mult=1))
+                print(Nmin)
+
+                # Test without capping
+                for N in range(Nmin,Nmax+1):
+                    print(N)
+
+                    # Evaluate
+                    x = clusterfunc.symmetric.unit_fixed_N(Dmin, Dmax, ER, N)
+
+                    # Check
+                    clusterfunc.check.unit_symmetric(x, Dmin, Dmax, ER, rtol=1e-9)
+
+                    dx = np.diff(x)
+                    rtol = 1e-9
+                    assert np.all(dx<= Dmax * (1.+rtol))
+                    assert np.all(dx* (1.+rtol)>= Dmin )
+                    assert np.diff(np.sign(dx)).all() == 0.
+
+def test_symmetric_abs():
+
+    rtol = 1e-9
+    ER = 1.2
+    for x0 in (-1., 0., 1., 2.):
+        for x1 in (-1., 0., 1., 2.):
+            for Dmin in  (1e-1, 1e-2, 1e-4):
+                for Dmax in (0.05, 0.1, 0.2, 1.):
+
+                    if Dmax <= Dmin:
+                        continue
+
+                    if np.isclose(x0,x1):
+                        with pytest.raises(ClusteringException):
+                            clusterfunc.symmetric.free_N(x0, x1, Dmin, Dmax, ER, mult=1)
+                    else:
+                        x = clusterfunc.symmetric.free_N(x0, x1, Dmin, Dmax, ER, mult=1)
+
+                        assert np.allclose(x[(0,-1),], (x0, x1))
+                        dx = np.diff(x)
+                        assert np.diff(np.sign(dx)).all() == 0.
+                        assert np.abs(dx).max() <= Dmax * (1.+rtol)
+                        assert np.abs(dx).min() * (1.+rtol) >= Dmin