M clusterfunc/check.py => clusterfunc/check.py +52 -32
@@ 2,84 2,104 @@ import numpy as np
from clusterfunc.exceptions import ClusteringException
import clusterfunc.util
+
def unit_single(x, Dmin, 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 x[0] == 0.0:
+ raise ClusteringException(f"Normalised start value x[0]={x[0]}, expected 0.0.")
- if not np.isclose(x[-1],1.,rtol=rtol):
- raise ClusteringException(f'Normalised end value x[-1]={x[-1]}, expected 1.0.')
+ if not np.isclose(x[-1], 1.0, rtol=rtol):
+ raise ClusteringException(f"Normalised end value x[-1]={x[-1]}, expected 1.0.")
- if not ((x>=0.) & (x<= 1. + rtol)).all():
+ if not ((x >= 0.0) & (x <= 1.0 + rtol)).all():
raise ClusteringException(
- f'Normalised values outside unit interval, min(x)={x.min()}, max(x)={x.max()}.'
- )
+ 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 (dx > 0.0).all():
+ raise ClusteringException(
+ f"Normalised spacing min(dx)={dx.min()}, expected all > 0.0."
+ )
if not np.isclose(dx[0], Dmin):
- raise ClusteringException(f'Normalised start spacing dx[0]={dx[0]}, expected {Dmin}')
+ raise ClusteringException(
+ f"Normalised start spacing dx[0]={dx[0]}, expected {Dmin}"
+ )
if not (dx <= Dmax).all():
- raise ClusteringException(f'Normalised spacing max(dx)={dx.max()} exceeds target {Dmax}.')
+ raise ClusteringException(
+ f"Normalised spacing max(dx)={dx.max()} exceeds target {Dmax}."
+ )
# Check expansion ratio
- ER = dx[1:]/dx[:-1]
+ ER = dx[1:] / dx[:-1]
if not (ER <= ERmax).all():
- raise ClusteringException(f'Expansion ratio max(ER)={ER.max()} exceeds target {ERmax}.')
+ raise ClusteringException(
+ f"Expansion ratio max(ER)={ER.max()} exceeds target {ERmax}."
+ )
- if not (ER >= 1. - rtol).all():
- raise ClusteringException(f'Expansion ratio min(ER)={ER.min()} less than unity.')
+ if not (ER >= 1.0 - rtol).all():
+ raise ClusteringException(
+ f"Expansion ratio min(ER)={ER.min()} less than unity."
+ )
def unit_symmetric(x, Dmin, Dmax, ERmax, rtol=1e-9):
-
# Check x
- if not x[0] == 0.:
- raise ClusteringException(f'Normalised start value x[0]={x[0]}, expected 0.0.')
+ if not x[0] == 0.0:
+ raise ClusteringException(f"Normalised start value x[0]={x[0]}, expected 0.0.")
- if not np.isclose(x[-1],1.,rtol=rtol):
- raise ClusteringException(f'Normalised end value x[-1]={x[-1]}, expected 1.0.')
+ if not np.isclose(x[-1], 1.0, rtol=rtol):
+ raise ClusteringException(f"Normalised end value x[-1]={x[-1]}, expected 1.0.")
- if not ((x>=0.) & (x<= 1. + rtol)).all():
+ if not ((x >= 0.0) & (x <= 1.0 + rtol)).all():
raise ClusteringException(
- f'Normalised values outside unit interval, min(x)={x.min()}, max(x)={x.max()}.'
- )
+ 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 (dx > 0.0).all():
+ raise ClusteringException(
+ f"Normalised spacing min(dx)={dx.min()}, expected all > 0.0."
+ )
if not np.isclose(dx[0], Dmin, rtol=rtol):
- raise ClusteringException(f'Normalised start spacing dx[0]={dx[0]}, expected {Dmin}')
+ raise ClusteringException(
+ f"Normalised start spacing dx[0]={dx[0]}, expected {Dmin}"
+ )
if not np.isclose(dx[-1], Dmin, rtol=rtol):
- raise ClusteringException(f'Normalised end spacing dx[-1]={dx[0]}, expected {Dmin}')
+ raise ClusteringException(
+ f"Normalised end spacing dx[-1]={dx[0]}, expected {Dmin}"
+ )
# Check expansion ratio
# Evaluate and correct for the shrinking part
- ER = dx[1:]/dx[:-1]
+ ER = dx[1:] / dx[:-1]
Na, Nb = clusterfunc.util.split_cells(len(x))
ER1 = ER.copy()
- ER1[ER1<1.] = 1./ER1[ER1<1.]
+ ER1[ER1 < 1.0] = 1.0 / ER1[ER1 < 1.0]
if not (ER1 <= ERmax).all():
- raise ClusteringException(f'Expansion ratio max(ER)={ER1.max()} exceeds target {ERmax}.')
+ raise ClusteringException(
+ f"Expansion ratio max(ER)={ER1.max()} exceeds target {ERmax}."
+ )
- if not (ER1 >= 1. - rtol).all():
- raise ClusteringException(f'Expansion ratio min(ER)={ER1.min()} less than unity.')
+ if not (ER1 >= 1.0 - rtol).all():
+ raise ClusteringException(
+ f"Expansion ratio min(ER)={ER1.min()} less than unity."
+ )
M clusterfunc/exceptions.py => clusterfunc/exceptions.py +1 -2
@@ 1,6 1,5 @@
"""Allow catching of clustering errors specifically."""
+
class ClusteringException(BaseException):
pass
-
-
M clusterfunc/plot.py => clusterfunc/plot.py +14 -11
@@ 1,26 1,29 @@
import matplotlib.pyplot as plt
import numpy as np
+
def plot(x):
- n = np.linspace(0, len(x)-1, len(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.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')
+ 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))
+ n = np.linspace(0, len(ER) - 1, len(ER))
fig, ax = plt.subplots()
- ax.plot(n, ER, 'k-x')
+ 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_xlim(n[0] - 1, n[-1])
# ax.set_ylim(x.min(), x.max())
- ax.set_xlabel('Index')
- ax.set_ylabel('Coordinate')
+ ax.set_xlabel("Index")
+ ax.set_ylabel("Coordinate")
+
import clusterfunc.single
import clusterfunc.symmetric
M clusterfunc/single.py => clusterfunc/single.py +134 -76
@@ 5,109 5,99 @@ from clusterfunc.exceptions import ClusteringException
from scipy.optimize import root_scalar
import clusterfunc.check
-def unit_fixed_N(Dmin, Dmax, ERmax, N, check=True):
- """Single-sided clustering on the unit interval with specified number of points.
- Parameters
- ----------
- Dmin: float
- Minimum spacing
-
- Returns
- -------
-
- """
-
- # There are four cases in decreasing order of grid points:
- # 1) Too many points for a uniform grid at Dmin.
- # 2) Can reach unit length without capping the spacing
- # 3) A section with capped spacing is required
- # 4) Not enough points to reach unit length at full stretch
- # """
+def _unit_fixed_N(Dmin, Dmax, ERmax, N, check=True):
+ """Single-sided clustering on the unit interval with specified number of points."""
if N < 2:
- raise ClusteringException(f'Need at least two points to cluster, N={N}')
+ raise ClusteringException(f"Need at least two points to cluster, N={N}")
if Dmax <= Dmin:
- raise ClusteringException(f'Dmax={Dmax} should be > Dmin={Dmin}')
+ raise ClusteringException(f"Dmax={Dmax} should be > Dmin={Dmin}")
if ERmax <= 1:
- raise ClusteringException(f'ERmax={ERmax} should be > 1')
+ raise ClusteringException(f"ERmax={ERmax} should be > 1")
# Rule out having too many points for a uniform spacing at Dmin
- Nunif = np.ceil(1./Dmin).astype(int)+1
+ Nunif = np.ceil(1.0 / Dmin).astype(int) + 1
if N > Nunif:
raise ClusteringException(
- f'Too many points N={N} for a uniform spacing of Dmin={Dmin}.'
+ f"Too many points N={N} for a uniform spacing of Dmin={Dmin}."
)
# Rule out too few points to reach unit length, regardless of capping
- Nfull = np.ceil(np.log(Dmax / Dmin) / np.log(ERmax)).astype(int)+1
- Lfull = Dmin * (1.-ERmax**(Nfull-1))/(1.-ERmax)
- if N < Nfull and Lfull < 1.:
+ Nfull = np.ceil(np.log(Dmax / Dmin) / np.log(ERmax)).astype(int) + 1
+ Lfull = Dmin * (1.0 - ERmax ** (Nfull - 1)) / (1.0 - ERmax)
+ if N < Nfull and Lfull < 1.0:
raise ClusteringException(
- f'Not enough points N={N} to reach unit length for ERmax={ERmax}, Dmin={Dmin}, regardless of Dmax., L={Lfull}, Nfull={Nfull}'
+ f"Not enough points N={N} to reach unit length for ERmax={ERmax}, Dmin={Dmin}, regardless of Dmax., L={Lfull}, Nfull={Nfull}"
)
# Assume we do not need to cap the spacing
# What expansion ratio is needed, and what is the max spacing?
- Nm1, Nm2, Nm3 = N-1, N-2, N-3
- Np1 = N+1
- Nm1Nm2 = Nm1*Nm2
+ Nm1, Nm2 = N - 1, N - 2
+ Np1 = N + 1
+
def _iter_ER(ERi):
- y = Dmin*(1.-ERi**Nm1)/(1.-ERi) - 1.
- dy = Dmin*(Nm2*ERi**Np1 - Nm1*ERi**N + ERi**2)/(1.-ERi)**2/ERi**2
+ y = Dmin * (1.0 - ERi**Nm1) / (1.0 - ERi) - 1.0
+ dy = (
+ Dmin
+ * (Nm2 * ERi**Np1 - Nm1 * ERi**N + ERi**2)
+ / (1.0 - ERi) ** 2
+ / ERi**2
+ )
return y, dy
# If we cannot reach irregardless of capping, give up
- if _iter_ER(ERmax)[0] < 0.:
+ if _iter_ER(ERmax)[0] < 0.0:
raise ClusteringException(
- f'Insufficient points N={N} to reach unit length with Dmin={Dmin}, ERmax={ERmax}, no capping to Dmax needed.'
+ f"Insufficient points N={N} to reach unit length with Dmin={Dmin}, ERmax={ERmax}, no capping to Dmax needed."
)
# Check for a uniform solution first
- err_unif = 1./(N-1) - Dmin
+ err_unif = 1.0 / (N - 1) - Dmin
rtol = 1e-9
if np.abs(err_unif) < rtol:
- ERnocap = 1.
- # Otherwise, find ER root
+ ERnocap = 1.0
+ # Otherwise, find ER root
else:
soln = root_scalar(_iter_ER, x0=ERmax, fprime=True, maxiter=1000)
assert soln.converged
ERnocap = soln.root
- Dnocap = Dmin*ERnocap**Nm2
+ Dnocap = Dmin * ERnocap**Nm2
# If capping constraint is satisfied, but ER is too big, we don't have enough points
- if (ERnocap > ERmax and Dnocap <= Dmax):
+ if ERnocap > ERmax and Dnocap <= Dmax:
raise ClusteringException(
- f'Insufficient points to reach unit length with Dmin={Dmin}, ERmax={ERmax}, no capping to Dmax needed.'
+ f"Insufficient points to reach unit length with Dmin={Dmin}, ERmax={ERmax}, no capping to Dmax needed."
)
# If ERnocap and Dnocap satisfy our constraints, we have the solution
elif ERnocap <= ERmax and Dnocap <= Dmax:
- dx = Dmin*ERnocap**np.arange(Nm1)
+ dx = Dmin * ERnocap ** np.arange(Nm1)
# We need to have a capped section
else:
-
# Length as a function of spacing cap, with ERmax constant
def _iter_Dmax(ERi, Dmaxi):
- Nstretch_i = np.ceil(np.log(Dmaxi / Dmin) / np.log(ERi)).astype(int)+1
- Lstretch_i = Dmin * (1.-ERi**(Nstretch_i-1))/(1.-ERi)
- dx1max_i = Dmin * ERi ** (Nstretch_i-2)
+ Nstretch_i = np.ceil(np.log(Dmaxi / Dmin) / np.log(ERi)).astype(int) + 1
+ Lstretch_i = Dmin * (1.0 - ERi ** (Nstretch_i - 1)) / (1.0 - ERi)
+ dx1max_i = Dmin * ERi ** (Nstretch_i - 2)
Ncap_i = N - Nstretch_i
Lcap_i = Ncap_i * 0.5 * (dx1max_i + Dmaxi)
- return Lstretch_i + Lcap_i - 1.
+ return Lstretch_i + Lcap_i - 1.0
# If we use all points at max ER and don't have enough length, will not work
LL = _iter_Dmax(ERmax, Dmax)
- if LL < 0.:
- raise ClusteringException(f'Not enough points to cluster with Dmin={Dmin}, ERmax={ERmax}, Dmax={Dmax} capping, total length only {LL+1}.')
+ if LL < 0.0:
+ raise ClusteringException(
+ f"Not enough points to cluster with Dmin={Dmin}, ERmax={ERmax}, Dmax={Dmax} capping, total length only {LL+1}."
+ )
# Binary search for cap distance that gives closes to unit length
Dmax_high = Dmax
Dmax_low = Dmin
- while Dmax_high-Dmax_low > 0.01*Dmax:
- Dmax_new = 0.5*(Dmax_high+Dmax_low)
+ while Dmax_high - Dmax_low > 0.01 * Dmax:
+ Dmax_new = 0.5 * (Dmax_high + Dmax_low)
err_new = _iter_Dmax(ERmax, Dmax_new)
if err_new > 0:
Dmax_high = Dmax_new
@@ 115,15 105,15 @@ def unit_fixed_N(Dmin, Dmax, ERmax, N, check=True):
Dmax_low = Dmax_new
# Fine-tune the expansion ratio for exactly unit length
- soln = root_scalar(_iter_Dmax, bracket=(1.+1e-9, ERmax), args=(Dmax_high,))
+ soln = root_scalar(_iter_Dmax, bracket=(1.0 + 1e-9, ERmax), args=(Dmax_high,))
assert soln.converged
ER = soln.root
# Evaluate the clustering
- Nstretch = np.ceil(np.log(Dmax_high / Dmin) / np.log(ER)).astype(int)+1
+ Nstretch = np.ceil(np.log(Dmax_high / Dmin) / np.log(ER)).astype(int) + 1
Ncap = N - Nstretch
- dx1 = Dmin * ER ** np.arange(0, Nstretch-1)
- dx2 = np.linspace(dx1[-1],Dmax_high, Ncap)
+ dx1 = Dmin * ER ** np.arange(0, Nstretch - 1)
+ dx2 = np.linspace(dx1[-1], Dmax_high, Ncap)
dx = np.concatenate((dx1, dx2))
x = clusterfunc.util.cumsum0(dx)
@@ 131,43 121,49 @@ def unit_fixed_N(Dmin, Dmax, ERmax, N, check=True):
x /= x[-1]
if not len(x) == N:
- raise ClusteringException(f'Incorrect number of points len(x)={len(x)}, expected {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
-def unit_free_N(Dmin, Dmax, ERmax, mult=8):
+def _unit_free_N(Dmin, Dmax, ERmax, mult=8):
"""Single-sided clustering on the unit interval with free number of points."""
# Find high and low guesses for n, N = mult*n + 1
- nlow = np.round(1./Dmax/mult).astype(int)
- nhigh = np.minimum(np.floor(1./Dmin/mult).astype(int),1024) # Limit to protect against overflow
+ nlow = np.round(1.0 / Dmax / mult).astype(int)
+ nhigh = np.minimum(
+ np.floor(1.0 / Dmin / mult).astype(int), 1024
+ ) # Limit to protect against overflow
if nlow > nhigh:
nlow = 1
if nhigh == nlow:
- raise ClusteringException(f'Unable to find distinct guesses for numbers of points with mult={mult}, try decreasing mult')
+ raise ClusteringException(
+ f"Unable to find distinct guesses for numbers of points with mult={mult}, try decreasing mult"
+ )
Nhigh = nhigh * mult + 1
try:
- x = unit_fixed_N(Dmin, Dmax, ERmax, Nhigh, check=True)
+ 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}')
+ raise ClusteringException(
+ f"Failed to cluster with high number of points guess Nhigh={Nhigh}"
+ )
# Loop until converged
- while nhigh-nlow > 1:
-
+ while nhigh - nlow > 1:
# Bisect high and low guesses
- nnew = (nhigh+nlow)//2
+ nnew = (nhigh + nlow) // 2
Nnew = nnew * mult + 1
# Attempt to cluster
try:
- x = unit_fixed_N(Dmin, Dmax, ERmax, Nnew, check=True)
+ x = _unit_fixed_N(Dmin, Dmax, ERmax, Nnew, check=True)
# New n is valid, use as high limiit
nhigh = nnew
except ClusteringException:
@@ 178,19 174,81 @@ def unit_free_N(Dmin, Dmax, ERmax, mult=8):
def fixed_N(x0, x1, Dmin, Dmax, ERmax, N):
- """Single-sided clustering between two values with with 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
+ """Single-sided clustering between two values with with fixed number of points.
+
+ Generate a grid vector x of length N from x0 to x1. Use geometric
+ stretching from a minimum spacing at x0, with spacings capped to a maximum
+ value.
+
+ Parameters
+ ----------
+ x0: float
+ Start value.
+ x1: float
+ End value.
+ Dmin: float
+ Minimum spacing at zero.
+ Dmax: float
+ Maximum spacing.
+ ERmax: float
+ Expansion ratio > 1.
+ N: int
+ Number of points in the grid vector.
+
+ Returns
+ -------
+ x: array
+ Grid vector of clustered 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)
+ x = x0 + dx * _unit_fixed_N(Dmin / dxa, Dmax / dxa, ERmax, N)
return x
+
def free_N(x0, x1, Dmin, Dmax, ERmax, mult=8):
- """Single-sided clustering between two values with with 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
+ """Single-sided clustering between two values with with free number of points.
+
+ Generate a grid vector x from x0 to x1. Use geometric stretching from a
+ minimum spacing at x0, with spacings capped to a maximum value. Use
+ the minimum number of points required to satsify the constraints.
+
+ Parameters
+ ----------
+ x0: float
+ Start value.
+ x1: float
+ End value.
+ Dmin: float
+ Minimum spacing at zero.
+ Dmax: float
+ Maximum spacing.
+ ERmax: float
+ Expansion ratio > 1.
+ mult: int
+ Choose a number of cells divisible by this factor.
+
+ Returns
+ -------
+ x: array
+ Grid vector of clustered points.
+
+ """
+
+ if np.isclose(x0, x1):
+ raise ClusteringException(
+ (
+ "Cannot distribute points without distinct start and end points, "
+ f"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)
+ x = x0 + dx * _unit_free_N(Dmin / dxa, Dmax / dxa, ERmax, mult)
+
return x
M clusterfunc/symmetric.py => clusterfunc/symmetric.py +48 -32
@@ 1,31 1,32 @@
"""Distribute points with symmetric clustering."""
import numpy as np
-from scipy.optimize import root_scalar
from clusterfunc.exceptions import ClusteringException
import clusterfunc.check
import clusterfunc.util
import clusterfunc.single
-def unit_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')
+ raise ClusteringException(
+ f"Not enough points N={N} to symmetrically cluster, need N>=3"
+ )
# Split number of cells as equally as possible
Na, Nb = clusterfunc.util.split_cells(N)
xj = 0.5
- def _guess_xj( Dmin, Dmax, ERmax, N, xj):
-
+ def _guess_xj(Dmin, Dmax, ERmax, N, xj):
# Initially assume exact halves
La = xj
- Lb = (1.-xj)
+ Lb = 1.0 - xj
# Evaluate the halves
- 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))
+ xa = clusterfunc.single.fixed_N(0.0, La, Dmin, Dmax, ERmax, Na)
+ xb1 = clusterfunc.single.fixed_N(0.0, Lb, Dmin, Dmax, ERmax, Nb)
+ xb = 1.0 - np.flip(xb1)
x = np.concatenate((xa, xb[1:]))
dxa = np.diff(xa)[-1]
@@ 33,59 34,74 @@ def unit_fixed_N(Dmin, Dmax, ERmax, N):
return x, dxa, dxb
- if np.mod(N,2):
- # Odd number exactly splits at 0.5
+ if np.mod(N, 2):
+ # Odd number exactly splits at 0.5
x, _, _ = _guess_xj(Dmin, Dmax, ERmax, N, xj)
else:
- # Move the join point if even number
+ # Move the join point if even number
maxiter = 20
for _ in range(maxiter):
-
x, dxa, dxb = _guess_xj(Dmin, Dmax, ERmax, N, xj)
if dxa > dxb:
- ERnow = dxa/dxb
+ ERnow = dxa / dxb
else:
- ERnow = dxb/dxa
+ ERnow = dxb / dxa
if ERnow < ERmax:
break
else:
- xj -= dxa-dxb
-
+ xj -= dxa - dxb
assert len(x) == N
return x
-def unit_free_N(Dmin, Dmax, ERmax, mult=8):
+
+def _unit_free_N(Dmin, Dmax, ERmax, mult=8):
""""""
- 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)
+ xa = clusterfunc.single.free_N(0.0, 0.5, Dmin, Dmax, ERmax, mult)
+ assert np.diff(xa).max() * (1.0 + 1e-9) <= Dmax
+ xb = 1.0 - np.flip(xa)
x = np.concatenate((xa, xb[1:]))
- assert not np.mod((len(x)-1)/mult,1.)
+ assert not np.mod((len(x) - 1) / mult, 1.0)
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
+ """Double-sided symmetric clustering with fixed number of points."""
+
+ if np.isclose(x0, x1):
+ raise ClusteringException(
+ (
+ "Cannot distribute points without distinct start and end, "
+ f"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)
+ 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
+ """Double-sided symmetric clustering with free number of points."""
+
+ if np.isclose(x0, x1):
+ raise ClusteringException(
+ (
+ "Cannot distribute points without distinct start and end, "
+ f"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
+ x = x0 + dx * _unit_free_N(Dmin / dxa, Dmax / dxa, ERmax, mult)
+ return x
D clusterfunc/unit.py => clusterfunc/unit.py +0 -73
@@ 1,73 0,0 @@
-"""Clustering on the unit interval."""
-import numpy as np
-from scipy.optimize import root_scalar
-from clusterfunc.exceptions import ClusteringException
-import clusterfunc.check
-import clusterfunc.util
-
-
-
-
-
-def symmetric_fixed_N(Dmin, Dmax, ERmax, N):
- """"""
- if N < 3:
- raise ClusteringException(f'Not enough points N={N} to symmetrically cluster, need N>=3')
-
- # Split number of cells as equally as possible
- Na, Nb = clusterfunc.util.split_cells(N)
-
- xj = 0.5
-
- def _guess_xj( Dmin, Dmax, ERmax, N, xj):
-
- # Initially assume exact halves
- La = xj
- 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
- x = np.concatenate((xa, xb[1:]))
-
- dxa = np.diff(xa)[-1]
- dxb = np.diff(xb)[0]
-
- return x, dxa, dxb
-
- if np.mod(N,2):
- # Odd number exactly splits at 0.5
-
- x, _, _ = _guess_xj(Dmin, Dmax, ERmax, N, xj)
-
- else:
- # Move the join point if even number
-
- maxiter = 20
- for _ in range(maxiter):
-
- x, dxa, dxb = _guess_xj(Dmin, Dmax, ERmax, N, xj)
-
- if dxa > dxb:
- ERnow = dxa/dxb
- else:
- ERnow = dxb/dxa
-
- if ERnow < ERmax:
- break
- else:
- xj -= dxa-dxb
-
-
- assert len(x) == N
-
- return x
-
-def symmetric_free_N(Dmin, Dmax, ERmax, mult=8):
- """"""
- xa = single_free_N(Dmin*2., Dmax*2., ERmax, mult)/2.
- xb = 1.-np.flip(xa)
- x = np.concatenate((xa, xb[1:]))
- assert not np.mod((len(x)-1)/mult,1.)
- return x
-
M clusterfunc/util.py => clusterfunc/util.py +5 -4
@@ 1,7 1,8 @@
import numpy as np
+
def split_cells(N):
- Na = (N-1)//2 + 1
+ Na = (N - 1) // 2 + 1
Nb = N - Na + 1
return Na, Nb
@@ 9,9 10,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.]
+ ER = dx[1:] / dx[:-1]
+ ER[ER < 1.0] = 1.0 / ER[ER < 1.0]
return ER
-
D tests/not_test_unit.py => tests/not_test_unit.py +0 -42
@@ 1,42 0,0 @@
-import clusterfunc.unit
-import clusterfunc.check
-import numpy as np
-
-def test_single():
-
- for Dmin in (2e-1, 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.unit.single_free_N(Dmin, Dmax, ER, mult=1))
-
- # Test without capping
- for N in range(Nmin,Nmax+1):
-
- # Evaluate
- x = clusterfunc.unit.single_fixed_N(Dmin, Dmax, ER, N)
-
- # Check
- clusterfunc.check.unit_single(x, Dmin, Dmax, ER, rtol=1e-9)
-
-def test_unif():
-
- N = np.arange(2, 101).astype(int)
- Dmax = 2.
- ER = 1.2
- for n in N:
- Dmin = 1./(n-1)
-
- # Evaluate
- x = clusterfunc.unit.single_fixed_N(Dmin, Dmax, ER, n)
-
- # Check
- clusterfunc.check.unit_single(x, Dmin, Dmax, ER, rtol=1e-9)
M tests/test_single.py => tests/test_single.py +32 -21
@@ 4,65 4,76 @@ import numpy as np
import pytest
from clusterfunc.exceptions import ClusteringException
+
def test_single():
+ """Verify clustering for single-sided unit interval."""
- for Dmin in (2e-1, 1e-1, 1e-2, 1e-3, 1e-4, 1e-5):
+ for Dmin in (2e-1, 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)
+ Nmax = np.minimum(np.floor(1.0 / Dmin).astype(int), 256)
- for Dmax in (0.05, 0.1, 0.2, 1.):
+ for Dmax in (0.05, 0.1, 0.2, 1.0):
if Dmax <= Dmin:
continue
# Minimum number of points with capping
- Nmin = len(clusterfunc.single.unit_free_N(Dmin, Dmax, ER, mult=1))
+ Nmin = len(
+ clusterfunc.single._unit_free_N(Dmin, Dmax, ER, mult=1)
+ )
# Test without capping
- for N in range(Nmin,Nmax+1):
+ for N in range(Nmin, Nmax + 1):
# Evaluate
- x = clusterfunc.single.unit_fixed_N(Dmin, Dmax, ER, N)
+ x = clusterfunc.single._unit_fixed_N(Dmin, Dmax, ER, N)
# Check
clusterfunc.check.unit_single(x, Dmin, Dmax, ER, rtol=1e-9)
+
def test_unif():
+ """Verify clustering for uniform cases on unit interval."""
N = np.arange(2, 101).astype(int)
- Dmax = 2.
+ Dmax = 2.0
ER = 1.2
for n in N:
- Dmin = 1./(n-1)
+ Dmin = 1.0 / (n - 1)
# Evaluate
- x = clusterfunc.single.unit_fixed_N(Dmin, Dmax, ER, n)
+ x = clusterfunc.single._unit_fixed_N(Dmin, Dmax, ER, n)
# Check
clusterfunc.check.unit_single(x, Dmin, Dmax, ER, rtol=1e-9)
+
def test_abs():
+ """Verify clustering for absolute coordinates."""
rtol = 1e-9
ER = 1.2
- for x0 in (-1., 0., 1., 2.):
- for x1 in (-1., 0., 1., 2.):
- for Dmin in (2e-1, 1e-2, 1e-4):
- for Dmax in (0.05, 0.1, 0.2, 1.):
-
+ for x0 in (-1.0, 0.0, 1.0, 2.0):
+ for x1 in (-1.0, 0.0, 1.0, 2.0):
+ for Dmin in (2e-1, 1e-2, 1e-4):
+ for Dmax in (0.05, 0.1, 0.2, 1.0):
if Dmax <= Dmin:
continue
- if np.isclose(x0,x1):
+ if np.isclose(x0, x1):
with pytest.raises(ClusteringException):
- clusterfunc.single.free_N(x0, x1, Dmin, Dmax, ER, mult=1)
+ clusterfunc.single.free_N(
+ x0, x1, Dmin, Dmax, ER, mult=1
+ )
else:
- x = clusterfunc.single.free_N(x0, x1, Dmin, Dmax, ER, mult=1)
+ x = clusterfunc.single.free_N(
+ x0, x1, Dmin, Dmax, ER, mult=1
+ )
- assert np.allclose(x[(0,-1),], (x0, x1))
+ 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.diff(np.sign(dx)).all() == 0.0
+ assert np.abs(dx).max() <= Dmax * (1.0 + rtol)
+ assert np.abs(dx).min() * (1.0 + rtol) >= Dmin
M tests/test_symmetric.py => tests/test_symmetric.py +30 -27
@@ 6,58 6,61 @@ from clusterfunc.exceptions import ClusteringException
def test_symmetric_unit():
-
- for Dmin in (1e-1, 1e-2, 1e-3, 1e-4, 1e-5):
+ 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.):
+ Nmax = np.minimum(np.floor(1.0 / Dmin).astype(int), 256)
+ for Dmax in (0.05, 0.1, 0.2, 1.0):
if Dmax <= Dmin:
continue
# Minimum number of points with capping
- Nmin = len(clusterfunc.symmetric.unit_free_N(Dmin, Dmax, ER, mult=1))
- print(Nmin)
+ Nmin = len(
+ clusterfunc.symmetric._unit_free_N(Dmin, Dmax, ER, mult=1)
+ )
# Test without capping
- for N in range(Nmin,Nmax+1):
+ for N in range(Nmin, Nmax + 1):
print(N)
# Evaluate
- x = clusterfunc.symmetric.unit_fixed_N(Dmin, Dmax, ER, N)
+ x = clusterfunc.symmetric._unit_fixed_N(Dmin, Dmax, ER, N)
# Check
- clusterfunc.check.unit_symmetric(x, Dmin, Dmax, ER, rtol=1e-9)
+ 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.
+ assert np.all(dx <= Dmax * (1.0 + rtol))
+ assert np.all(dx * (1.0 + rtol) >= Dmin)
+ assert np.diff(np.sign(dx)).all() == 0.0
-def test_symmetric_abs():
+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.):
-
+ for x0 in (-1.0, 0.0, 1.0, 2.0):
+ for x1 in (-1.0, 0.0, 1.0, 2.0):
+ for Dmin in (1e-1, 1e-2, 1e-4):
+ for Dmax in (0.05, 0.1, 0.2, 1.0):
if Dmax <= Dmin:
continue
- if np.isclose(x0,x1):
+ if np.isclose(x0, x1):
with pytest.raises(ClusteringException):
- clusterfunc.symmetric.free_N(x0, x1, Dmin, Dmax, ER, mult=1)
+ clusterfunc.symmetric.free_N(
+ x0, x1, Dmin, Dmax, ER, mult=1
+ )
else:
- x = clusterfunc.symmetric.free_N(x0, x1, Dmin, Dmax, ER, mult=1)
+ x = clusterfunc.symmetric.free_N(
+ x0, x1, Dmin, Dmax, ER, mult=1
+ )
- assert np.allclose(x[(0,-1),], (x0, x1))
+ 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.diff(np.sign(dx)).all() == 0.0
+ assert np.abs(dx).max() <= Dmax * (1.0 + rtol)
+ assert np.abs(dx).min() * (1.0 + rtol) >= Dmin