2 files changed, 150 insertions(+), 0 deletions(-)

A clusterfunc.py
A test_clusterfunc.py

@@ 1,136 @@
+"""Clustering functions for geometric stretching with capped spacing."""
+import numpy as np
+from scipy.optimize import root_scalar
+
+class ClusteringException(BaseException):
+    """Allow catching of clustering errors specifically."""
+    pass
+
+def check_normalised(x, Dst, Dmax, ERmax, rtol=1e-9):
+    """Verify that a normalised clustering vector obeys the desired limits, raise if not."""
+
+
+    # Check x
+
+    if not x[0] == 0.:
+        raise ClusteringException(f'Normalised start value x[0]={x[0]}, expected 0.0.')
+
+    if not np.isclose(x[-1],1.):
+        raise ClusteringException(f'Normalised end value x[-1]={x[-1]}, expected 1.0.')
+
+    if not ((x>=0.) & (x<= 1. + rtol)).all():
+        raise ClusteringException(
+            f'Normalised values outside unit interval, min(x)={x.min()}, max(x)={x.max()}.'
+            )
+
+    # Check dx
+
+    dx = np.diff(x)
+
+    if not (dx > 0.).all():
+        raise ClusteringException(f'Normalised spacing min(dx)={dx.min()}, expected all > 0.0.')
+
+    if not np.isclose(dx[0], Dst):
+        raise ClusteringException(f'Normalised start spacing dx[0]={dx[0]}, expected {Dst}')
+
+    if not (dx <= Dmax).all():
+        raise ClusteringException(f'Normalised spacing max(dx)={dx.max()} exceeds target {Dmax}.')
+
+    # Check expansion ratio
+
+    ER = dx[1:]/dx[:-1]
+
+    if not (ER <= ERmax).all():
+        raise ClusteringException(f'Expansion ratio max(ER)={ER.max()} exceeds target {ERmax}.')
+
+def normalised_N_fixed(Dst, Dmax, ERmax, N, check=True):
+    """Single-sided clustering on the unit interval with specified number of point."""
+
+    # Number of points needed to expand Dst to Dmax at full stretch
+    # Expanding for more points will give D > Dmax, which we don't want
+    N1a = np.floor(np.log(Dmax / Dst) / np.log(ERmax)).astype(int)+1
+
+    # Length at which we will need to cap the spacing
+    L1a = Dst * (1.-ERmax**(N1a-1))/(1.-ERmax)
+
+    # Number of points for unit length at full stretch
+    # Any less than this, we cannot reach
+    N1b = np.ceil(np.log(1.-(1.-ERmax)/Dst)/np.log(ERmax)).astype(int) + 1
+
+    print(N1a, N1b, N, L1a)
+
+    if N < N1b:
+        raise ClusteringException(
+        f'Insufficient points to reach unit length with Dst={Dst}, ERmax={ERmax} (for any Dmax)'
+        )
+
+    if N <= N1a or L1a > 1.:
+        # No need to cap spacing
+        # Reduce expansion ratio to achieve unit length
+
+        # Closure to return L and dLdER as function of ER
+        def _iter(ERi):
+            # Length error
+            L = Dst * (1.0 - ERi ** (N - 1)) / (1.0 - ERi) - 1.0
+            # Derivative
+            dL = (
+                Dst
+                * ((N - 2.0) * ERi**N - (N - 1.0) * ERi ** (N - 1.0) + ERi)
+                / (1.0 - ERi) ** 2.0
+                / ERi
+            )
+            return L, dL
+
+        # Solve for the ER that gives unit length
+        ER = root_scalar(_iter, x0=ERmax, fprime=True).root
+        print(ER)
+
+        # Evaluate coordinates
+        dx = Dst * ER ** np.arange(0, N - 1)
+        x = cumsum0(dx)
+
+    else:
+        # If we used all N points, D would exceed Dmax, so uniform needed
+        raise Exception('Spacing needs to be capped')
+
+    if check:
+        if not len(x) == N:
+            raise ClusteringException(f'Incorrect number of points len(x)={len(x)}, expected {N}.')
+        check_normalised(x, Dst, Dmax, ERmax)
+
+    return x
+
+
+def cumsum0(x, axis=None):
+    return np.insert(np.cumsum(x, axis=axis), 0, 0.0, axis=axis)
+
+
+import matplotlib.pyplot as plt
+
+
+
+
+Dst = 0.001
+Dmax = 0.1
+ER = 1.2
+
+N = 31
+
+xx = normalised_N_fixed(Dst, Dmax, ER, N)
+
+fig, ax = plt.subplots()
+ax.plot(xx,'k-x')
+ax.set_title('x')
+plt.show()
+
+fig, ax = plt.subplots()
+ax.set_title('dx')
+ax.plot(np.diff(xx),'k-x')
+plt.show()
+
+# fig, ax = plt.subplots()
+# ax.plot(dx,'k-x')
+# plt.show()
+# fig, ax = plt.subplots()
+# ax.plot(x,'k-x')
+# plt.show()

@@ 1,14 @@
+import clusterfunc
+
+def test_normalised_N_fixed():
+
+    Dst = 0.001
+    Dmax = 0.5
+    ER = 1.2
+    for N in range(32,100):
+
+        clusterfunc.normalised_N_fixed(Dst, Dmax, ER, N)
+
+
+if __name__=='__main__':
+    test_normalised_N_fixed()