Skip to content

Commit

Permalink
blackified everything
Browse files Browse the repository at this point in the history
  • Loading branch information
albop committed Aug 5, 2022
1 parent 14f6591 commit ce2420f
Show file tree
Hide file tree
Showing 14 changed files with 100 additions and 76 deletions.
4 changes: 2 additions & 2 deletions examples/example_mlinterp.py
Expand Up @@ -62,7 +62,7 @@ def vec_eval(u):

x1 = np.linspace(0, 1, 100) ** 2 # non-uniform points
x2 = np.linspace(0, 1, 100) ** 2 # non-uniform points
y = np.array([[np.sqrt(u1 ** 2 + u2 ** 2) for u2 in x2] for u1 in x1])
y = np.array([[np.sqrt(u1**2 + u2**2) for u2 in x2] for u1 in x1])
# (y[i,j] = sqrt(x1[i]**2+x2[j]**2)


Expand Down Expand Up @@ -112,7 +112,7 @@ def vec_eval(p):
x1 = np.linspace(0, 1, 100) ** 2 # non-uniform points for first dimensoin
x2 = (0.0, 1.0, 100) # uniform points for second dimension
grid = (x1, x2)
y = np.array([[np.sqrt(u1 ** 2 + u2 ** 2) for u2 in x2] for u1 in x1])
y = np.array([[np.sqrt(u1**2 + u2**2) for u2 in x2] for u1 in x1])


points = np.random.random((1000, 2))
Expand Down
2 changes: 1 addition & 1 deletion interpolation/smolyak/grid.py
Expand Up @@ -786,7 +786,7 @@ def __init__(self, d, mu, lb=None, ub=None):
def __repr__(self):
npoints = self.cube_grid.shape[0]
nz_pts = np.count_nonzero(self.B)
pct_nz = nz_pts / (npoints ** 2.0)
pct_nz = nz_pts / (npoints**2.0)

if isinstance(self.mu, int):
msg = "Smolyak Grid:\n\td: {0} \n\tmu: {1} \n\tnpoints: {2}"
Expand Down
4 changes: 2 additions & 2 deletions interpolation/smolyak/tests/test_interp.py
Expand Up @@ -5,13 +5,13 @@
import numpy as np

# func = lambda x, y: np.exp(x**2 - y**2)
func = lambda x, y: x ** 2 - y ** 2
func = lambda x, y: x**2 - y**2


func1 = lambda points: func(points[:, 0], points[:, 1])
func1_prime = lambda x: np.column_stack([2 * x[:, 0], -2 * x[:, 1]])

func2 = lambda x: np.sum(x ** 2, axis=1)
func2 = lambda x: np.sum(x**2, axis=1)
func2_prime = lambda x: 2 * x


Expand Down
18 changes: 9 additions & 9 deletions interpolation/splines/__init__.py
Expand Up @@ -12,24 +12,24 @@
def UCGrid(*args):
tt = numba.typeof((10.0, 1.0, 1))
for a in args:
assert(numba.typeof(a) == tt)
assert numba.typeof(a) == tt
min, max, n = a
assert(min<max)
assert(n>1)
assert min < max
assert n > 1

return tuple(args)


def CGrid(*args):
tt = numba.typeof((10.0, 1.0, 1))
for a in args:
if isinstance(a, np.ndarray):
assert(a.ndim==1)
assert(a.shape[0]>2)
elif (numba.typeof(a) == tt):
assert a.ndim == 1
assert a.shape[0] > 2
elif numba.typeof(a) == tt:
min, max, n = a
assert(min<max)
assert(n>1)
assert min < max
assert n > 1
else:
raise Exception(f"Unknown dimension specification: {a}")

Expand Down
84 changes: 52 additions & 32 deletions interpolation/splines/hermite.py
Expand Up @@ -3,23 +3,27 @@
import numpy.typing as npt
from typing import Tuple


class Vector:
pass


@jit(nopython=True)
def hermite_splines(lambda0: float)->Tuple[float, float, float, float]:
def hermite_splines(lambda0: float) -> Tuple[float, float, float, float]:
"""Computes the cubic Hermite splines in lambda0
Inputs: - float: lambda0
Output: - tuple: cubic Hermite splines evaluated in lambda0"""
h00 = 2*(lambda0**3) - 3*(lambda0**2) + 1
h10 = (lambda0**3) - 2*(lambda0**2) + lambda0
h01 = -2*(lambda0**3) + 3*(lambda0**2)
h00 = 2 * (lambda0**3) - 3 * (lambda0**2) + 1
h10 = (lambda0**3) - 2 * (lambda0**2) + lambda0
h01 = -2 * (lambda0**3) + 3 * (lambda0**2)
h11 = (lambda0**3) - (lambda0**2)
return (h00, h10, h01, h11)


@jit(nopython=True)
def hermite_interp(x0: float, xk: float, xkn: float, pk: float, pkn: float, mk: float, mkn: float)->float:
def hermite_interp(
x0: float, xk: float, xkn: float, pk: float, pkn: float, mk: float, mkn: float
) -> float:
"""Returns the interpolated value for x0.
Inputs: - float: x0, abscissa of the point to interpolate
- float: xk, abscissa of the nearest lowest point to x0 on the grid
Expand All @@ -30,9 +34,14 @@ def hermite_interp(x0: float, xk: float, xkn: float, pk: float, pkn: float, mk:
- float: mkn, tangent in xkn
Output: - float: interpolated value for x0
"""
t = (x0-xk)/(xkn-xk)
t = (x0 - xk) / (xkn - xk)
hsplines = hermite_splines(t)
return (pk*hsplines[0] + mk*(xkn-xk)*hsplines[1] + pkn*hsplines[2] + mkn*(xkn-xk)*hsplines[3])
return (
pk * hsplines[0]
+ mk * (xkn - xk) * hsplines[1]
+ pkn * hsplines[2]
+ mkn * (xkn - xk) * hsplines[3]
)


@jit(nopython=True)
Expand All @@ -48,12 +57,12 @@ def HermiteInterpolation(x0: float, x, y, yp):
return y[0]
elif x0 >= np.max(x):
return y[-1]

###### Interpolation case ######
indx = np.searchsorted(x, x0)
xk, xkn = x[indx-1], x[indx]
pk, pkn = y[indx-1], y[indx]
mk, mkn = yp[indx-1], yp[indx]
xk, xkn = x[indx - 1], x[indx]
pk, pkn = y[indx - 1], y[indx]
mk, mkn = yp[indx - 1], yp[indx]
return hermite_interp(x0, xk, xkn, pk, pkn, mk, mkn)


Expand All @@ -72,43 +81,54 @@ def HermiteInterpolationVect(xvect, x: Vector, y: Vector, yp: Vector):
out[i] = HermiteInterpolation(x0, x, y, yp)
return out


from numba import njit, types
from numba.extending import overload, register_jitable
from numba import generated_jit


def _hermite(x0,x,y,yp,out=None):
def _hermite(x0, x, y, yp, out=None):
pass


@overload(_hermite)
def _hermite(x0,x,y,yp,out=None):
def _hermite(x0,x,y,yp,out=None):
return HermiteInterpolation(x0,x,y,yp)
def _hermite(x0, x, y, yp, out=None):
def _hermite(x0, x, y, yp, out=None):
return HermiteInterpolation(x0, x, y, yp)

return _hermite


from numba.core.types.misc import NoneType as none


@generated_jit
def hermite(x0,x,y,yp,out=None):
def hermite(x0, x, y, yp, out=None):
try:
n = x0.ndim
if n==1:
input_type = 'vector'
elif n==2:
input_type = 'matrix'
if n == 1:
input_type = "vector"
elif n == 2:
input_type = "matrix"
else:
raise Exception("Invalid input type")
except:
# n must be a scalar
input_type = 'scalar'

if input_type == 'scalar':
def _hermite(x0,x,y,yp,out=None):
return HermiteInterpolation(x0,x,y,yp)
elif input_type == 'vector':
def _hermite(x0,x,y,yp,out=None):
return HermiteInterpolationVect(x0,x,y,yp)
elif input_type == 'matrix':
def _hermite(x0,x,y,yp,out=None):
return HermiteInterpolationVect(x0[:,0],x,y,yp)
return _hermite
input_type = "scalar"

if input_type == "scalar":

def _hermite(x0, x, y, yp, out=None):
return HermiteInterpolation(x0, x, y, yp)

elif input_type == "vector":

def _hermite(x0, x, y, yp, out=None):
return HermiteInterpolationVect(x0, x, y, yp)

elif input_type == "matrix":

def _hermite(x0, x, y, yp, out=None):
return HermiteInterpolationVect(x0[:, 0], x, y, yp)

return _hermite
16 changes: 6 additions & 10 deletions interpolation/splines/tests/test_derivatives.py
Expand Up @@ -4,7 +4,6 @@ def test_derivatives():

import numpy as np


grid = ((0.0, 1.0, 10),)

# grid = ((0.0, 1.0, 0.1),(0.0, 1.0, 0.1))
Expand All @@ -13,7 +12,6 @@ def test_derivatives():
Cx = np.concatenate([C[:, None], C[:, None] * 2])
points = np.random.random((10, 1))


eval_spline(
grid, C, (-0.1,), out=None, order=1, diff="None", extrap_mode="nearest"
) # no alloc
Expand All @@ -24,7 +22,6 @@ def test_derivatives():
grid, C, (-0.1,), out=None, order=1, diff="None", extrap_mode="linear"
) # no alloc


eval_spline(
grid, C, (1.1,), out=None, order=1, diff="None", extrap_mode="nearest"
) # no alloc
Expand All @@ -35,16 +32,13 @@ def test_derivatives():
grid, C, (1.1,), out=None, order=1, diff="None", extrap_mode="linear"
) # no alloc


eval_spline(
grid, Cx, points[0, :], out=None, order=1, diff="None", extrap_mode="linear"
)


eval_spline(grid, C, points, out=None, order=1, diff="None", extrap_mode="linear")
eval_spline(grid, Cx, points, out=None, order=1, diff="None", extrap_mode="linear")


orders = str(((0,), (1,)))

eval_spline(
Expand All @@ -56,18 +50,18 @@ def test_derivatives():
eval_spline(grid, C, points, out=None, order=1, diff=orders, extrap_mode="linear")
eval_spline(grid, Cx, points, out=None, order=1, diff=orders, extrap_mode="linear")


out = eval_spline(
grid, Cx, points, out=None, order=1, diff=orders, extrap_mode="linear"
)
out2 = np.zeros_like(out)
eval_spline(grid, Cx, points, out=out2, order=1, diff=orders, extrap_mode="linear")
print(abs(out - out2).max())


k = 3

eval_spline(grid, C, points[0, :], out=None, order=3, diff="None", extrap_mode="linear")
eval_spline(
grid, C, points[0, :], out=None, order=3, diff="None", extrap_mode="linear"
)

eval_spline(grid, C, points, out=None, order=k, diff="None", extrap_mode="linear")
eval_spline(
Expand All @@ -77,7 +71,9 @@ def test_derivatives():

orders = str(((0,), (1,)))

eval_spline(grid, C, points[0, :], out=None, order=k, diff=orders, extrap_mode="linear")
eval_spline(
grid, C, points[0, :], out=None, order=k, diff=orders, extrap_mode="linear"
)
eval_spline(grid, C, points, out=None, order=k, diff=orders, extrap_mode="linear")
eval_spline(
grid, Cx, points[0, :], out=None, order=k, diff=orders, extrap_mode="linear"
Expand Down
17 changes: 8 additions & 9 deletions interpolation/splines/tests/test_hermite_splines.py
@@ -1,22 +1,21 @@
def test_hermite_splines():

from interpolation.splines.hermite import HermiteInterpolationVect
import numpy as np
N = 10000 # Number of points in the initial dataset

N = 10000 # Number of points in the initial dataset
K = 1000 # Number of new points to interpolate

# Initial dataset
# grid = ((0.0, 1.0, K),) # Creation of an x-axis grid (xi)
grid = np.linspace(0.0, 1.0, N) # Creation of an x-axis grid (xi)
points = np.random.random((N)) # Random values for f(xi)
points = np.random.random((N)) # Random values for f(xi)
dpoints = np.random.random((N)) # Random derivatives for f'(xi)

# Generate new points
newgrid = np.random.random((K))

# Interpolation
out = HermiteInterpolationVect(newgrid, grid, points, dpoints)

print("OK")

2 changes: 1 addition & 1 deletion interpolation/splines/tests/test_multilinear_extrap.py
Expand Up @@ -14,7 +14,7 @@ def test_multilinear_extrap():

s = mlinspace(a, b, n)

f = lambda x: (x ** 2).sum(axis=1)
f = lambda x: (x**2).sum(axis=1)

x = f(s)
v = x.reshape(n)
Expand Down
2 changes: 1 addition & 1 deletion interpolation/tests/test_complete.py
Expand Up @@ -35,7 +35,7 @@ def f(x, y):
return x

def f2(x, y):
return x ** 3 - y
return x**3 - y

points = np.random.random((1000, 2))
vals = np.column_stack(
Expand Down
19 changes: 14 additions & 5 deletions interpolation/tests/test_derivs.py
Expand Up @@ -4,7 +4,7 @@


class Check1DDerivatives(unittest.TestCase):
"""
"""
Checks derivatives in a 1D interpolator
"""

Expand Down Expand Up @@ -35,7 +35,10 @@ def test_linear(self):
# 0-order must be the function
# 1-order must be the slope
result = np.vstack(
[y0 + slope * eval_points, np.ones_like(eval_points) * slope,]
[
y0 + slope * eval_points,
np.ones_like(eval_points) * slope,
]
).T

self.assertTrue(np.allclose(grad, result))
Expand All @@ -61,7 +64,10 @@ def test_nonlinear(self):
# 0-order must be the function
# 1-order must be + or - pi/2
result = np.vstack(
[np.array([0, -1, 0, 1, 0]), np.array([-1, -1, 1, 1, -1]) * 2 / np.pi,]
[
np.array([0, -1, 0, 1, 0]),
np.array([-1, -1, 1, 1, -1]) * 2 / np.pi,
]
).T

self.assertTrue(np.allclose(grad, result))
Expand All @@ -87,14 +93,17 @@ def test_nonlinear_approx(self):
# 0-order must be x^3
# 1-order must be close to 3x^2
result = np.vstack(
[np.power(eval_points, 3), np.power(eval_points, 2) * 3.0,]
[
np.power(eval_points, 3),
np.power(eval_points, 2) * 3.0,
]
).T

self.assertTrue(np.allclose(grad, result, atol=0.02))


class Check2DDerivatives(unittest.TestCase):
"""
"""
Checks derivatives in a 2D interpolator
"""

Expand Down

0 comments on commit ce2420f

Please sign in to comment.