~jb753/clusterfunc

9b109bd9c53d7bdc753482adec6175d6ace650eb — James Brind 10 months ago 312638f
Cleanup and sample
M clusterfunc/check.py => clusterfunc/check.py +16 -7
@@ 4,19 4,24 @@ 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."""
    """Verify that unit clustering vector obeys the desired limits."""

    # Check x

    if not x[0] == 0.0:
        raise ClusteringException(f"Normalised start value x[0]={x[0]}, expected 0.0.")
        raise ClusteringException(
            f"Normalised start value x[0]={x[0]}, expected 0.0."
        )

    if not np.isclose(x[-1], 1.0, rtol=rtol):
        raise ClusteringException(f"Normalised end value x[-1]={x[-1]}, expected 1.0.")
        raise ClusteringException(
            f"Normalised end value x[-1]={x[-1]}, expected 1.0."
        )

    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, "
            f"min(x)={x.min()}, max(x)={x.max()}."
        )

    # Check dx


@@ 27,7 32,7 @@ def unit_single(x, Dmin, Dmax, ERmax, rtol=1e-9):
            f"Normalised spacing min(dx)={dx.min()}, expected all > 0.0."
        )

    if not np.isclose(dx[0], Dmin):
    if not np.isclose(dx[0], Dmin, rtol):
        raise ClusteringException(
            f"Normalised start spacing dx[0]={dx[0]}, expected {Dmin}"
        )


@@ 56,10 61,14 @@ def unit_symmetric(x, Dmin, Dmax, ERmax, rtol=1e-9):
    # Check x

    if not x[0] == 0.0:
        raise ClusteringException(f"Normalised start value x[0]={x[0]}, expected 0.0.")
        raise ClusteringException(
                f"Normalised start value x[0]={x[0]}, expected 0.0."
                )

    if not np.isclose(x[-1], 1.0, rtol=rtol):
        raise ClusteringException(f"Normalised end value x[-1]={x[-1]}, expected 1.0.")
        raise ClusteringException(
                f"Normalised end value x[-1]={x[-1]}, expected 1.0."
                )

    if not ((x >= 0.0) & (x <= 1.0 + rtol)).all():
        raise ClusteringException(

M clusterfunc/double.py => clusterfunc/double.py +140 -132
@@ 9,6 9,90 @@ import clusterfunc.plot
import warnings


def fixed(dx0, dx1, N, x0=0.0, x1=1.0):
    """Double-sided clustering with fixed number of points.

    Generate a grid vector x of length N, by default over the unit interval.
    Use Vinokur stretching from specified spacings at x0 and x1. Expansion
    ratio and maximum spacing are not controlled.

    Parameters
    ----------
    dx0 float
        Boundary spacing at x0.
    dx1 float
        Boundary spacing at x1.
    N: int
        Number of points in the grid vector.
    x0: float
        Start value.
    x1: float
        End value.

    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}"
        )

    # Scale the unit distribution between the given values
    Dx = x1 - x0
    Dxa = np.abs(Dx)
    x = x0 + Dx * _unit_fixed(dx0 / Dxa, dx1 / Dxa, N)

    return x


def free(dx0, dx1, dmax, ERmax, x0=0.0, x1=1.0, mult=8):
    """Double-sided clustering between two values with with fixed number of points.

    Generate a grid vector x from x0 to x1. Use Vinokur stretching from
    specified spacings at x0 and x1. Increase the number of points until
    maximum spacing and expansion ratio criteria are satisfied.

    Parameters
    ----------
    dx0 float
        Boundary spacing at x0.
    dx1 float
        Boundary spacing at x1.
    dmax: float
        Maximum spacing.
    ERmax: float
        Expansion ratio > 1.
    x0: float
        Start value.
    x1: float
        End value.
    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(
            f"Cannot distribute points without distinct start and end points, got x0={x0} and x1={x1}"
        )

    # Scale the unit distribution between the given values
    Dx = x1 - x0
    Dxa = np.abs(Dx)
    x = x0 + Dx * _unit_free(dx0 / Dxa, dx1 / Dxa, dmax / Dxa, ERmax, mult)

    return x


def _vinokur(ds, N):
    """Two sided analytic clustering function after Vinokur."""



@@ 50,58 134,63 @@ def _vinokur(ds, N):
def _invert_sinhx_x(y):
    """Return solution x for y = sinh(x)/x in Eqns. (62-67)."""

    y1 = y - 1.0
    x_low = np.sqrt(6.0 * y1) * (
        1.0
        - 0.15 * y1
        + 0.057321429 * y1**2.0
        - 0.024907295 * y1**3.0
        + 0.0077424461 * y1**4.0
        - 0.0010794123 * y1**5.0
    )

    v = np.log(y)
    w = 1.0 / y - 0.028527431
    x_high = (
        v
        + (1.0 + 1.0 / v) * np.log(2.0 * v)
        - 0.02041793
        + 0.24902722 * w
        + 1.9496443 * w**2.0
        - 2.6294547 * w**3.0
        + 8.56795911 * w**4.0
    )

    return np.where(y < 2.7829681, x_low, x_high)
    if y < 2.7829681:
        y1 = y - 1.0
        x2 = np.sqrt(6.0 * y1) * (
            1.0
            - 0.15 * y1
            + 0.057321429 * y1**2.0
            - 0.024907295 * y1**3.0
            + 0.0077424461 * y1**4.0
            - 0.0010794123 * y1**5.0
        )
    else:
        v = np.log(y)
        w = 1.0 / y - 0.028527431
        x2 = (
            v
            + (1.0 + 1.0 / v) * np.log(2.0 * v)
            - 0.02041793
            + 0.24902722 * w
            + 1.9496443 * w**2.0
            - 2.6294547 * w**3.0
            + 8.56795911 * w**4.0
        )

    # if not np.isclose(x, x2):
    #     print(x, x2, y)
    #     raise Exception('beans')

    return x2


def _invert_sinx_x(y):
    """Return solution x for y = sin(x)/x from Eqns. (68-71)."""
    if y < 0.26938972:
        x = np.pi * (
            1.0
            - y
            + y**2.0
            - (1.0 + np.pi**2.0 / 6.0) * y**3.0
            + 6.794732 * y**4.0
            - 13.205501 * y**5.0
            + 11.726095 * y**6.0
        )
    else:
        y1 = 1.0 - y
        x = np.sqrt(6.0 * y1) * (
            1.0
            + 0.15 * y1
            + 0.057321429 * y1**2.0
            + 0.048774238 * y1**3.0
            - 0.053337753 * y1**4.0
            + 0.075845134 * y1**5.0
        )

    return x

    x_low = np.pi * (
        1.0
        - y
        + y**2.0
        - (1.0 + np.pi**2.0 / 6.0) * y**3.0
        + 6.794732 * y**4.0
        - 13.205501 * y**5.0
        + 11.726095 * y**6.0
    )

    y1 = 1.0 - y
    x_high = np.sqrt(6.0 * y1) * (
        1.0
        + 0.15 * y1
        + 0.057321429 * y1**2.0
        + 0.048774238 * y1**3.0
        - 0.053337753 * y1**4.0
        + 0.075845134 * y1**5.0
    )

    return np.where(y < 0.26938972, x_low, x_high)


def _unit_fixed(dx0, dx1, N, rtol=1e-3):

def _unit_fixed(dx0, dx1, N, rtol=1e-2):
    """General double-side clustering with fixed number of points"""
    maxiter = 100
    dx = np.array([dx0, dx1])


@@ 117,7 206,7 @@ def _unit_fixed(dx0, dx1, N, rtol=1e-3):
    return x


def _unit_free(dx0, dx1, dmax, ERmax, mult=8, rtol=1e-3):
def _unit_free(dx0, dx1, dmax, ERmax, mult=8, rtol=1e-2):
    """General double-side clustering with free number of points"""
    n = 1
    maxiter = 1000


@@ 126,6 215,9 @@ def _unit_free(dx0, dx1, dmax, ERmax, mult=8, rtol=1e-3):
        warnings.filterwarnings("error")
        for _ in range(maxiter):
            N = mult * n + 1
            if N < 4:
                n += 1
                continue
            try:
                x = _unit_fixed(dx0, dx1, N, rtol)
            except RuntimeWarning:


@@ 143,87 235,3 @@ def _unit_free(dx0, dx1, dmax, ERmax, mult=8, rtol=1e-3):
        raise ClusteringException("Could not double cluster")

    return x


def fixed(x0, x1, dx0, dx1, N):
    """Double-sided clustering between two values with with fixed number of points.

    Generate a grid vector x of length N from x0 to x1. Use Vinokur stretching
    from specified spacings at x0 and x1. Expansion ratio and maximum spacing
    are not controlled.

    Parameters
    ----------
    x0: float
        Start value.
    x1: float
        End value.
    dx0 float
        Boundary spacing at x0.
    dx1 float
        Boundary spacing at x1.
    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}"
        )

    # Scale the unit distribution between the given values
    Dx = x1 - x0
    Dxa = np.abs(Dx)
    x = x0 + Dx * _unit_fixed(dx0 / Dxa, dx1 / Dxa, N)

    return x


def free(x0, x1, dx0, dx1, dmax, ERmax, mult=8):
    """Double-sided clustering between two values with with fixed number of points.

    Generate a grid vector x from x0 to x1. Use Vinokur stretching from
    specified spacings at x0 and x1. Increase the number of points until
    maximum spacing and expansion ratio criteria are satisfied.

    Parameters
    ----------
    x0: float
        Start value.
    x1: float
        End value.
    dx0 float
        Boundary spacing at x0.
    dx1 float
        Boundary spacing at x1.
    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(
            f"Cannot distribute points without distinct start and end points, got x0={x0} and x1={x1}"
        )

    # Scale the unit distribution between the given values
    Dx = x1 - x0
    Dxa = np.abs(Dx)
    x = x0 + Dx * _unit_free(dx0 / Dxa, dx1 / Dxa, dmax / Dxa, ERmax, mult)

    return x

M clusterfunc/single.py => clusterfunc/single.py +88 -87
@@ 6,6 6,94 @@ from scipy.optimize import root_scalar
import clusterfunc.check


def fixed(dmin, dmax, ERmax, N, x0=0.0, x1=1.0, check=True):
    """Single-sided clustering with fixed number of points.

    Generate a grid vector x of length N, by default over the unit interval.
    Use geometric stretching from a minimum start spacing, with spacings capped
    to a maximum value.

    Parameters
    ----------
    dmin: float
        Minimum spacing at start.
    dmax: float
        Maximum spacing.
    ERmax: float
        Expansion ratio > 1.
    x0: float
        Start value.
    x1: float
        End value.
    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}"
        )

    # Scale the unit distribution between the given values
    Dx = x1 - x0
    Dxa = np.abs(Dx)
    x = x0 + Dx * _unit_fixed(dmin / Dxa, dmax / Dxa, ERmax, N, check)

    return x


def free(dmin, dmax, ERmax, x0=0.0, x1=1.0, mult=8):
    """Single-sided clustering with free number of points.

    Generate a grid vector x, by default over the unit interval. Use geometric
    stretching from a minimum start spacing, with spacings capped to a maximum
    value. Use the minimum number of points required to satsify the
    constraints.

    Parameters
    ----------
    dmin: float
        Minimum spacing at zero.
    dmax: float
        Maximum spacing.
    ERmax: float
        Expansion ratio > 1.
    x0: float
        Start value.
    x1: float
        End value.
    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}"
            )
        )

    # Scale the unit distribution between the given values
    Dx = x1 - x0
    Dxa = np.abs(Dx)
    x = x0 + Dx * _unit_free(dmin / Dxa, dmax / Dxa, ERmax, mult)

    return x


def _unit_fixed(dmin, dmax, ERmax, N, check=True):
    """Single-sided clustering on the unit interval with specified number of points."""



@@ 171,90 259,3 @@ def _unit_free(dmin, dmax, ERmax, mult=8):
            nlow = nnew

    return x


def fixed(x0, x1, dmin, dmax, ERmax, N, check=True):
    """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 start.
    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}"
        )

    # Scale the unit distribution between the given values
    Dx = x1 - x0
    Dxa = np.abs(Dx)
    x = x0 + Dx * _unit_fixed(dmin / Dxa, dmax / Dxa, ERmax, N, check)

    return x


def free(x0, x1, dmin, dmax, ERmax, mult=8):
    """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}"
            )
        )

    # Scale the unit distribution between the given values
    Dx = x1 - x0
    Dxa = np.abs(Dx)
    x = x0 + Dx * _unit_free(dmin / Dxa, dmax / Dxa, ERmax, mult)

    return x

M clusterfunc/symmetric.py => clusterfunc/symmetric.py +59 -102
@@ 3,105 3,62 @@ import numpy as np
from clusterfunc.exceptions import ClusteringException
import clusterfunc.check
import clusterfunc.util
import clusterfunc.single


def _unit_fixed(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.0 - xj

        # Evaluate the halves
        xa = clusterfunc.single.fixed(0.0, La, Dmin, Dmax, ERmax, Na)
        xb1 = clusterfunc.single.fixed(0.0, Lb, Dmin, Dmax, ERmax, Nb)
        xb = 1.0 - np.flip(xb1)
        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 _unit_free(Dmin, Dmax, ERmax, mult=8):
    """"""
    xa = clusterfunc.single.free(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.0)
    return x


def fixed(x0, x1, Dmin, Dmax, ERmax, N):
    """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(Dmin / dxa, Dmax / dxa, ERmax, N)

    return x


def free(x0, x1, Dmin, Dmax, ERmax, mult=8):
    """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(Dmin / dxa, Dmax / dxa, ERmax, mult)

    return x
import clusterfunc.double


def fixed(dmin, N, x0=0.0, x1=1.0):
    """Double-sided clustering between two values with with fixed number of points.

    Generate a grid vector x of length N, defaulting to the unit interval. Use
    Vinokur stretching from specified minimum end spacing. Expansion ratio and
    maximum spacing are not controlled.

    Parameters
    ----------
    dmin float
        Boundary spacing.
    N: int
        Number of points in the grid vector.
    x0: float
        Start value.
    x1: float
        End value.

    Returns
    -------
    x: array
        Grid vector of clustered points.

    """
    return clusterfunc.double.fixed(dmin, dmin, N, x0, x1)


def free(dmin, dmax, ERmax, x0=0.0, x1=1.0, mult=8):
    """Symmetric clustering between two values with with fixed number of points.

    Generate a grid vector x, by default over the unit interval. Use Vinokur
    stretching from specified minimum spacing at both ends. Increase the number
    of points until maximum spacing and expansion ratio criteria are satisfied.

    Parameters
    ----------
    dmin float
        Boundary spacing at both ends.
    dmax: float
        Maximum spacing.
    ERmax: float
        Expansion ratio > 1.
    x0: float
        Start value.
    x1: float
        End value.
    mult: int
        Choose a number of cells divisible by this factor.

    Returns
    -------
    x: array
        Grid vector of clustered points.

    """
    return clusterfunc.double.free(dmin, dmin, dmax, ERmax, x0, x1, mult)

D compare.py => compare.py +0 -19
@@ 1,19 0,0 @@
import clusterfunc.symmetric
import clusterfunc.double
import clusterfunc.plot

dmin = 2e-3
dmax = 1e-1
ERmax = 1.2

xd = clusterfunc.double.free(0.0, 1.0, dmin, dmin, dmax, ERmax)
xs = clusterfunc.symmetric.free(0.0, 1.0, dmin, dmax, ERmax)
print(len(xd), len(xs))

clusterfunc.plot.plot(xd)
clusterfunc.plot.plot_ER(xd)
clusterfunc.plot.plot(xs)
clusterfunc.plot.plot_ER(xs)
import matplotlib.pyplot as plt

plt.show()

A sample.py => sample.py +18 -0
@@ 0,0 1,18 @@
import clusterfunc.symmetric
import clusterfunc.double
import clusterfunc.plot
import matplotlib.pyplot as plt

dmin = 2e-3
dmax = 2e-1
ERmax = 1.2

x1 = clusterfunc.single.free(dmin, dmax, ERmax)
x2 = clusterfunc.double.free(dmin, 5.*dmin, dmax, ERmax)
x3 = clusterfunc.symmetric.free(dmin, dmax, ERmax)

clusterfunc.plot.plot(x1)
clusterfunc.plot.plot(x2)
clusterfunc.plot.plot(x3)

plt.show()

M tests/test_double.py => tests/test_double.py +2 -2
@@ 56,11 56,11 @@ def test_double_abs():
                        if np.isclose(x0, x1):
                            with pytest.raises(ClusteringException):
                                clusterfunc.double.free(
                                    x0, x1, Dmin0, Dmin1, Dmax, ERmax, mult=1
                                    Dmin0, Dmin1, Dmax, ERmax, x0, x1, mult=1
                                )
                        else:
                            x = clusterfunc.double.free(
                                x0, x1, Dmin0, Dmin1, Dmax, ERmax, mult=1
                                Dmin0, Dmin1, Dmax, ERmax, x0, x1, mult=1
                            )

                            assert np.allclose(x[(0, -1),], (x0, x1))

M tests/test_single.py => tests/test_single.py +2 -2
@@ 59,9 59,9 @@ def test_single_abs():

                    if np.isclose(x0, x1):
                        with pytest.raises(ClusteringException):
                            clusterfunc.single.free(x0, x1, Dmin, Dmax, ER, mult=1)
                            clusterfunc.single.free(Dmin, Dmax, ER, x0, x1, mult=1)
                    else:
                        x = clusterfunc.single.free(x0, x1, Dmin, Dmax, ER, mult=1)
                        x = clusterfunc.single.free(Dmin, Dmax, ER, x0, x1, mult=1)

                        assert np.allclose(x[(0, -1),], (x0, x1))
                        dx = np.diff(x)

M tests/test_symmetric.py => tests/test_symmetric.py +6 -6
@@ 16,15 16,15 @@ def test_symmetric_unit():
                    continue

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

                # Test without capping
                for N in range(Nmin, Nmax + 1, 3):
                    # Evaluate
                    x = clusterfunc.symmetric._unit_fixed(Dmin, Dmax, ER, N)
                    x = clusterfunc.symmetric.fixed(Dmin, N)

                    # Check
                    clusterfunc.check.unit_symmetric(x, Dmin, Dmax, ER, rtol=1e-9)
                    clusterfunc.check.unit_symmetric(x, Dmin, Dmax, ER, rtol=1e-2)

                    dx = np.diff(x)
                    rtol = 1e-9


@@ 34,7 34,7 @@ def test_symmetric_unit():


def test_symmetric_abs():
    rtol = 1e-9
    rtol = 1e-3
    ER = 1.2
    for x0 in (-1.0, 0.0, 2.0):
        for x1 in (-1.0, 0.0, 2.0):


@@ 45,9 45,9 @@ def test_symmetric_abs():

                    if np.isclose(x0, x1):
                        with pytest.raises(ClusteringException):
                            clusterfunc.symmetric.free(x0, x1, Dmin, Dmax, ER, mult=1)
                            clusterfunc.symmetric.free(Dmin, Dmax, ER, x0, x1, mult=1)
                    else:
                        x = clusterfunc.symmetric.free(x0, x1, Dmin, Dmax, ER, mult=1)
                        x = clusterfunc.symmetric.free(Dmin, Dmax, ER, x0, x1, mult=1)

                        assert np.allclose(x[(0, -1),], (x0, x1))
                        dx = np.diff(x)