~jb753/clusterfunc

25e94837bb7291ef07b108ece20dba7c0d4d264b — James Brind 11 months ago 6f95b4a
Improvements
13 files changed, 471 insertions(+), 125 deletions(-)

A .gitignore
A clusterfunc/__init__.py
A clusterfunc/check.py
A clusterfunc/exceptions.py
R clusterfunc.py => clusterfunc/single.py
A clusterfunc/symmetric.py
A clusterfunc/unit.py
A clusterfunc/util.py
A pyproject.toml
A pytest.ini
D test_clusterfunc.py
A tests/not_test_unit.py
A tests/test_single.py
A .gitignore => .gitignore +12 -0
@@ 0,0 1,12 @@
tags
dist
doc/runs/
tags*
__pycache__
doc/_build
.idea
*.pyc
.python-version
*.pdf
*.npz
*.egg-info/

A clusterfunc/__init__.py => clusterfunc/__init__.py +0 -0
A clusterfunc/check.py => clusterfunc/check.py +85 -0
@@ 0,0 1,85 @@
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 np.isclose(x[-1],1.,rtol=rtol):
        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], 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}.')

    # 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}.')

    if not (ER >= 1. - 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 np.isclose(x[-1],1.,rtol=rtol):
        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], Dmin, rtol=rtol):
        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}')

    # Check expansion ratio

    # Evaluate and correct for the shrinking part
    ER = dx[1:]/dx[:-1]
    Na, Nb = clusterfunc.util.split_cells(len(x))

    ER1 = ER.copy()
    ER1[ER1<1.] = 1./ER1[ER1<1.]

    if not (ER1 <= ERmax).all():
        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.')

A clusterfunc/exceptions.py => clusterfunc/exceptions.py +6 -0
@@ 0,0 1,6 @@
"""Allow catching of clustering errors specifically."""

class ClusteringException(BaseException):
    pass



R clusterfunc.py => clusterfunc/single.py +55 -90
@@ 1,59 1,28 @@
"""Clustering functions for geometric stretching with capped spacing."""
"""Distribute points with single-sided clustering."""
import numpy as np
import clusterfunc.util
from clusterfunc.exceptions import ClusteringException
from scipy.optimize import root_scalar

class ClusteringException(BaseException):
    """Allow catching of clustering errors specifically."""
    pass
def unit_fixed_N(Dmin, Dmax, ERmax, N):
    """Single-sided clustering on the unit interval with specified number of points.

def check_normalised(x, Dmin, Dmax, ERmax, rtol=1e-9):
    """Verify that a normalised clustering vector obeys the desired limits, raise if not."""
    Parameters
    ----------
    Dmin: float
        Minimum spacing 

    # Check x
    Returns
    -------

    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.,rtol=rtol):
        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], 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}.')

    # 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}.')

    if not (ER >= 1. - rtol).all():
        raise ClusteringException(f'Expansion ratio min(ER)={ER.min()} less than unity.')

def unit_fixed_N(Dmin, Dmax, ERmax, N, check=True):
    """Single-sided clustering on the unit interval with specified number of point.

    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
        """
    # 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
    #     """

    if N < 2:
        raise ClusteringException(f'Need at least two points to cluster, N={N}')


@@ 85,25 54,27 @@ def unit_fixed_N(Dmin, Dmax, ERmax, N, check=True):
    Np1 = N+1
    Nm1Nm2 = Nm1*Nm2
    def _iter_ER(ERi):
        # y = ERi - 1. + Dmin*(1.-ERi**Nm1)
        # dy = 1. - Dmin*Nm1*ERi**Nm2
        # d2y = -Dmin*Nm1Nm2*ERi**Nm3
        y = Dmin*(1.-ERi**Nm1)/(1.-ERi) - 1.
        dy = Dmin*(Nm2*ERi**Np1 - Nm1*ERi**N + ERi**2)/(1.-ERi)**2/ERi**2
        return y, dy#, d2y
        return y, dy

    # If we cannot reach irregardless of capping, give up
    # print(N,_iter_ER(ERmax)[0])
    # quit()
    if _iter_ER(ERmax)[0] > 0.:
    if _iter_ER(ERmax)[0] < 0.:
        raise ClusteringException(
        f'Insufficient points N={N} to reach unit length with Dmin={Dmin}, ERmax={ERmax}, no capping to Dmax needed.'
        )

    # Find root 
    soln = root_scalar(_iter_ER, x0=ERmax, fprime=True, maxiter=200)
    assert soln.converged
    ERnocap = soln.root
    # Check for a uniform solution first
    err_unif = 1./(N-1) - Dmin
    rtol = 1e-9
    if np.abs(err_unif) < rtol:
        ERnocap = 1.
    # 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

    # If capping constraint is satisfied, but ER is too big, we don't have enough points


@@ 154,30 125,32 @@ def unit_fixed_N(Dmin, Dmax, ERmax, N, check=True):
        dx2 = np.linspace(dx1[-1],Dmax_high, Ncap)
        dx = np.concatenate((dx1, dx2))

    x = cumsum0(dx)
    x = clusterfunc.util.cumsum0(dx)

    x /= x[-1]

    if check:
        if not len(x) == N:
            raise ClusteringException(f'Incorrect number of points len(x)={len(x)}, expected {N}.')
        check_normalised(x, Dmin, Dmax, ERmax)
    if not len(x) == N:
        raise ClusteringException(f'Incorrect number of points len(x)={len(x)}, expected {N}.')

    return x

def unit_free_N(Dmin, Dmax, ERmax, mult=8, check=True):

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.round(1./Dmin/mult).astype(int),512) # Limit to protect against overflow
    nhigh = np.minimum(np.ceil(1./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')

    Nhigh = nhigh * mult + 1
    x = unit_fixed_N(Dmin, Dmax, ERmax, Nhigh, check=False)
    try:
        x = unit_fixed_N(Dmin, Dmax, ERmax, Nhigh)
    except ClusteringException:
        raise ClusteringException(f'Failed to cluster with high number of points guess Nhigh={Nhigh}')

    # Loop until converged
    while nhigh-nlow > 1:


@@ 188,38 161,30 @@ def unit_free_N(Dmin, Dmax, ERmax, mult=8, check=True):

        # Attempt to cluster
        try:
            x = unit_fixed_N(Dmin, Dmax, ERmax, Nnew, check=False)
            x = unit_fixed_N(Dmin, Dmax, ERmax, Nnew)
            # New n is valid, use as high limiit
            nhigh = nnew
        except ClusteringException:
            # New n is not valid, use as low limiit
            nlow = nnew

    if check:
        check_normalised(x, Dmin, Dmax, ERmax)

    return x


def fixed_N(x0, x1, Dmin, Dmax, ERmax, N, check=True):
def fixed_N(x0, x1, Dmin, Dmax, ERmax, N):
    """Single-sided clustering between two values with with fixed number of points."""
    return x0 + (x1-x0) * unit_fixed_N(Dmin, Dmax, ERmax, N, check)
    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, check=True):
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
    return x0 + (x1-x0) * unit_free_N(Dmin/dx, Dmax/dx, ERmax, mult, check)

def symmetric_unit_fixed_N(Dmin, Dmax, ERmax, N, check=True):
    Na = N//2
    Nb = N - Na + 1
    xa = fixed_N(0., 0.5, Dmin, Dmax, ERmax, Na, check)
    xb = np.flip(fixed_N(1., 0.5, Dmin, Dmax, ERmax, Na, check))
    x = np.concatenate((xa, xb[1:]))
    assert len(x) == N


def cumsum0(x, axis=None):
    return np.insert(np.cumsum(x, axis=axis), 0, 0.0, axis=axis)

# symmetric_unit_fixed_N(1e-3,1e-2, 1.2, 1000)
    dxa = np.abs(dx)
    x = x0 + dx * unit_free_N(Dmin/dxa, Dmax/dxa, ERmax, mult)
    return x

A clusterfunc/symmetric.py => clusterfunc/symmetric.py +70 -0
@@ 0,0 1,70 @@
"""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


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


A clusterfunc/unit.py => clusterfunc/unit.py +73 -0
@@ 0,0 1,73 @@
"""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


A clusterfunc/util.py => clusterfunc/util.py +11 -0
@@ 0,0 1,11 @@
import numpy as np

def split_cells(N):
    Na = (N-1)//2 + 1
    Nb = N - Na + 1
    return Na, Nb


def cumsum0(x, axis=None):
    return np.insert(np.cumsum(x, axis=axis), 0, 0.0, axis=axis)


A pyproject.toml => pyproject.toml +21 -0
@@ 0,0 1,21 @@
[project]
name = "clusterfunc"
version = "0.0.1"
requires-python = ">=3.9"
dependencies = [
    'numpy',
    'scipy',
]

[project.optional-dependencies]
tests = [
    'pytest~=6.2.5',
]

[build-system]
requires = ["setuptools"]
build-backend = "setuptools.build_meta"

# # Only package what is in the turbigen folder
# [tool.setuptools.packages.find]
# include = ["turbigen",]

A pytest.ini => pytest.ini +3 -0
@@ 0,0 1,3 @@
[pytest]
testpaths = tests
filterwarnings = error

D test_clusterfunc.py => test_clusterfunc.py +0 -35
@@ 1,35 0,0 @@
import clusterfunc
import numpy as np

def test_cluster():

    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), 512)

            for Dmax in (0.05, 0.1, 0.2, 1.):

                if Dmax <= Dmin:
                    continue

                # Minimum number of points with capping
                Nmin = len(clusterfunc.unit_free_N(Dmin, Dmax, ER, mult=1))

                # Test without capping
                for N in range(Nmin,Nmax+1):
                    clusterfunc.unit_fixed_N(Dmin, Dmax, ER, N)

def test_unif():

    N = np.arange(2, 101).astype(int)
    Dmax = 2.
    ER = 1.2
    for n in N:
        Dmin = 1./(n-1)
        clusterfunc.unit_fixed_N(Dmin, Dmax, ER, n)

if __name__=='__main__':
    test_cluster()
    # test_unif()

A tests/not_test_unit.py => tests/not_test_unit.py +67 -0
@@ 0,0 1,67 @@
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)

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)

A tests/test_single.py => tests/test_single.py +68 -0
@@ 0,0 1,68 @@
import clusterfunc.single
import clusterfunc.check
import numpy as np
import pytest
from clusterfunc.exceptions import ClusteringException

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.single.unit_free_N(Dmin, Dmax, ER, mult=1))

                # Test without capping
                for N in range(Nmin,Nmax+1):

                    # Evaluate
                    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():

    N = np.arange(2, 101).astype(int)
    Dmax = 2.
    ER = 1.2
    for n in N:
        Dmin = 1./(n-1)

        # Evaluate
        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():

    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.):

                    if Dmax <= Dmin:
                        continue

                    if np.isclose(x0,x1):
                        with pytest.raises(ClusteringException):
                            clusterfunc.single.free_N(x0, x1, Dmin, Dmax, ER, mult=1)
                    else:
                        x = clusterfunc.single.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