Skip to content

Commit

Permalink
formating using black
Browse files Browse the repository at this point in the history
  • Loading branch information
bjodah committed Aug 7, 2021
1 parent 244c3a5 commit bfb1940
Show file tree
Hide file tree
Showing 24 changed files with 520 additions and 299 deletions.
43 changes: 25 additions & 18 deletions examples/err.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,27 +14,34 @@ def demo_err():
"""
max_order = 7
n = 20
l = 0.25
fmt1 = '{0: <5s}\t{1: <21s}\t{2: >21s}\t{3: >21s}\t{4: >21s}'
fmt2 = '{0: <5d}\t{1:20.18f}\t{2: >21.18f}\t{3: >21.18f}\t{4: >21.18f}'
x = np.cumsum(np.random.rand(n)*l)
x = np.concatenate((x[::-1]*-1, x))
lp = 0.25
fmt1 = "{0: <5s}\t{1: <21s}\t{2: >21s}\t{3: >21s}\t{4: >21s}"
fmt2 = "{0: <5d}\t{1:20.18f}\t{2: >21.18f}\t{3: >21.18f}\t{4: >21.18f}"
x = np.cumsum(np.random.rand(n) * lp)
x = np.concatenate((x[::-1] * -1, x))
lst = []
derivs = np.zeros(n)
for order in range(max_order+1):
print('Order', order)
for m in range(1+order//2, n+1):
sub_x = x[n-m:n+m]
derivs[m-1] = derivatives_at_point_by_finite_diff(
sub_x, np.exp(sub_x), 0, order)[order]
print(fmt1.format('m', 'val', 'diff', 'analytical error',
'diff/analytical'))
for order in range(max_order + 1):
print("Order", order)
for m in range(1 + order // 2, n + 1):
sub_x = x[n - m : n + m]
derivs[m - 1] = derivatives_at_point_by_finite_diff(
sub_x, np.exp(sub_x), 0, order
)[order]
print(fmt1.format("m", "val", "diff", "analytical error", "diff/analytical"))
for m in range(1, n):
print(fmt2.format(
(m+1)*2, derivs[m], derivs[m]-derivs[m-1],
derivs[m]-1, (derivs[m]-derivs[m-1])/(derivs[m]-1)))
lst.append((derivs[-1], abs(derivs[-1]-derivs[-2])))
print(
fmt2.format(
(m + 1) * 2,
derivs[m],
derivs[m] - derivs[m - 1],
derivs[m] - 1,
(derivs[m] - derivs[m - 1]) / (derivs[m] - 1),
)
)
lst.append((derivs[-1], abs(derivs[-1] - derivs[-2])))
print(np.array(lst))

if __name__ == '__main__':

if __name__ == "__main__":
demo_err()
50 changes: 26 additions & 24 deletions examples/sine.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,9 +5,7 @@

import numpy as np

from finitediff import (
derivatives_at_point_by_finite_diff, interpolate_by_finite_diff
)
from finitediff import derivatives_at_point_by_finite_diff, interpolate_by_finite_diff


def demo_usage(n_data=50, n_fit=537, nhead=5, ntail=5, plot=False, alt=0):
Expand All @@ -28,64 +26,68 @@ def demo_usage(n_data=50, n_fit=537, nhead=5, ntail=5, plot=False, alt=0):

x0, xend = 0, 5
# shaky linspace -5% to +5% noise
x_data = (np.linspace(x0, xend, n_data) +
np.random.rand(n_data)*(xend-x0)/n_data/1.5)
y_data = np.sin(x_data) * (1.0+0.1*(np.random.rand(n_data)-0.5))
x_data = (
np.linspace(x0, xend, n_data)
+ np.random.rand(n_data) * (xend - x0) / n_data / 1.5
)
y_data = np.sin(x_data) * (1.0 + 0.1 * (np.random.rand(n_data) - 0.5))

x_fit = np.linspace(x0, xend, n_fit)

# Edges behave badly, work around:
x_fit[0] = x_fit[0] + (x_fit[1]-x_fit[0])/2
x_fit[-1] = x_fit[-2]+(x_fit[-1]-x_fit[-2])/2
x_fit[0] = x_fit[0] + (x_fit[1] - x_fit[0]) / 2
x_fit[-1] = x_fit[-2] + (x_fit[-1] - x_fit[-2]) / 2

if alt:
y_fit = np.empty(n_fit)
dydx_fit = np.empty(n_fit)
for i, xf in enumerate(x_fit):
# get index j of first data point beyond xf
j = np.where(x_data > xf)[0][0]
lower_bound = max(0, j-alt)
upper_bound = min(n_data-1, j+alt)
lower_bound = max(0, j - alt)
upper_bound = min(n_data - 1, j + alt)
y_fit[i] = derivatives_at_point_by_finite_diff(
x_data[lower_bound:upper_bound],
y_data[lower_bound:upper_bound], xf, 0)
x_data[lower_bound:upper_bound], y_data[lower_bound:upper_bound], xf, 0
)
dydx_fit[i] = derivatives_at_point_by_finite_diff(
x_data[lower_bound:upper_bound],
y_data[lower_bound:upper_bound], xf, 1)[1]
x_data[lower_bound:upper_bound], y_data[lower_bound:upper_bound], xf, 1
)[1]
else:
interp = interpolate_by_finite_diff(x_data, y_data, x_fit,
1, nhead, ntail)
interp = interpolate_by_finite_diff(x_data, y_data, x_fit, 1, nhead, ntail)
y_fit = interp[:, 0]
dydx_fit = interp[:, 1]

if plot:
import matplotlib.pyplot as plt

plt.subplot(221)
plt.plot(x_data, y_data, 'x', label='Data points (sin)')
plt.plot(x_fit, y_fit, '-', label='Fitted curve (order=0)')
plt.plot(x_data, np.sin(x_data), '-', label='Analytic sin(x)')
plt.plot(x_data, y_data, "x", label="Data points (sin)")
plt.plot(x_fit, y_fit, "-", label="Fitted curve (order=0)")
plt.plot(x_data, np.sin(x_data), "-", label="Analytic sin(x)")
plt.legend()

plt.subplot(222)
plt.plot(x_fit, y_fit-np.sin(x_fit), label='Error in order=0')
plt.plot(x_fit, y_fit - np.sin(x_fit), label="Error in order=0")
plt.legend()

plt.subplot(223)
plt.plot(x_fit, dydx_fit, '-', label='Fitted derivative (order=1)')
plt.plot(x_data, np.cos(x_data), '-', label='Analytic cos(x)')
plt.plot(x_fit, dydx_fit, "-", label="Fitted derivative (order=1)")
plt.plot(x_data, np.cos(x_data), "-", label="Analytic cos(x)")
plt.legend()

plt.subplot(224)
plt.plot(x_fit, dydx_fit-np.cos(x_fit), label='Error in order=1')
plt.plot(x_fit, dydx_fit - np.cos(x_fit), label="Error in order=1")
plt.legend()

plt.show()

if __name__ == '__main__':

if __name__ == "__main__":
try:
from argh import dispatch_command
except ImportError:

def dispatch_command(cb):
return cb()

dispatch_command(demo_usage)
6 changes: 3 additions & 3 deletions examples/tests/test_examples.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,11 +14,11 @@
import pytest


tests = glob.glob(os.path.join(os.path.dirname(__file__), '../*.py'))
tests = glob.glob(os.path.join(os.path.dirname(__file__), "../*.py"))


@pytest.mark.parametrize('pypath', tests)
@pytest.mark.parametrize("pypath", tests)
def test_examples(pypath):
py_exe = 'python3' if sys.version_info.major == 3 else 'python'
py_exe = "python3" if sys.version_info.major == 3 else "python"
p = subprocess.Popen([py_exe, pypath])
assert p.wait() == 0 # SUCCESS==0
17 changes: 11 additions & 6 deletions finitediff/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,19 +2,24 @@
"""
Finite difference weights for any derivative order on arbitrarily spaced grids.
"""
from __future__ import (absolute_import, division, print_function)
from __future__ import absolute_import, division, print_function

from ._release import __version__

from ._finitediff_c import (
derivatives_at_point_by_finite_diff, interpolate_by_finite_diff,
get_weights
derivatives_at_point_by_finite_diff,
interpolate_by_finite_diff,
get_weights,
)

__all__ = ['derivatives_at_point_by_finite_diff', 'interpolate_by_finite_diff', 'get_weights']
__all__ = [
"derivatives_at_point_by_finite_diff",
"interpolate_by_finite_diff",
"get_weights",
]


def get_include():
from pkg_resources import resource_filename, Requirement
return resource_filename(Requirement.parse(__name__),
'%s/include' % __name__)

return resource_filename(Requirement.parse(__name__), "%s/include" % __name__)
2 changes: 1 addition & 1 deletion finitediff/_config.py
Original file line number Diff line number Diff line change
@@ -1 +1 @@
env = {'USE_FORTRAN': '0'}
env = {"USE_FORTRAN": "0"}
2 changes: 1 addition & 1 deletion finitediff/_release.py
Original file line number Diff line number Diff line change
@@ -1 +1 @@
__version__ = '0.6.5.dev0+git'
__version__ = "0.6.5.dev0+git"
7 changes: 6 additions & 1 deletion finitediff/grid/__init__.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,9 @@
from .rebalance import rebalanced_grid, pre_pruning_mask, combine_grids, grid_pruning_mask
from .rebalance import (
rebalanced_grid,
pre_pruning_mask,
combine_grids,
grid_pruning_mask,
)
from .refine import refine_grid
from .make import adapted_grid
from .util import locate_discontinuity, pool_discontinuity_approx, grid_error
Expand Down
4 changes: 2 additions & 2 deletions finitediff/grid/make.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
# -*- coding: utf-8 -*-
from __future__ import (absolute_import, division, print_function)
from __future__ import absolute_import, division, print_function

import numpy as np
from .refine import refine_grid


def adapted_grid(xstart, xstop, cb, grid_additions=(50, 50), **kwargs):
"""" Creates an adapted (1D) grid by subsequent subgrid insertions.
""" " Creates an adapted (1D) grid by subsequent subgrid insertions.
Parameters
----------
Expand Down
26 changes: 17 additions & 9 deletions finitediff/grid/plotting.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# -*- coding: utf-8 -*-
from __future__ import (absolute_import, division, print_function)
from __future__ import absolute_import, division, print_function

import numpy as np
from .make import adapted_grid
Expand All @@ -10,22 +10,30 @@ def plot_convergence(key, values, cb, metric=None, xstart=0, xstop=2, **kwargs):

if key in kwargs:
raise ValueError("Cannot have key=%s when given in kwargs" % key)
fig, axes = plt.subplots(1, len(values), figsize=(16, 5),
sharey=True, gridspec_kw={'wspace': 0})
fig, axes = plt.subplots(
1, len(values), figsize=(16, 5), sharey=True, gridspec_kw={"wspace": 0}
)
scores, grid_sizes = [], []
if key is None and len(values) != 1:
raise ValueError("Not considering key")
for val, ax in zip(values, np.atleast_1d(axes)):
if key is not None:
kwargs[key] = val
grid, results = adapted_grid(xstart, xstop, cb, metric=metric, **kwargs)
y = results if metric is None else np.array(
[metric(r) for r in results])
ax.vlines(grid, 0, 1, transform=ax.get_xaxis_transform(),
linestyle='--', color='k', alpha=.3, linewidth=.5)
y = results if metric is None else np.array([metric(r) for r in results])
ax.vlines(
grid,
0,
1,
transform=ax.get_xaxis_transform(),
linestyle="--",
color="k",
alpha=0.3,
linewidth=0.5,
)
ax.plot(grid, y)
between_x = grid[:-1] + np.diff(grid)/2
between_y = y[:-1] + np.diff(y)/2
between_x = grid[:-1] + np.diff(grid) / 2
between_y = y[:-1] + np.diff(y) / 2
rbx = cb(between_x)
ybx = rbx if metric is None else np.array([metric(r) for r in rbx])
scores.append(np.sum(np.abs(between_y - ybx)))
Expand Down
45 changes: 25 additions & 20 deletions finitediff/grid/rebalance.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
# -*- coding: utf-8 -*-
from __future__ import (absolute_import, division, print_function)
from __future__ import absolute_import, division, print_function

import math
import numpy as np
Expand All @@ -10,31 +10,36 @@ def _avgdiff(x):
dx = np.diff(x)
dx2 = np.zeros_like(x)
dx2[0], dx2[-1] = dx[0], dx[-1]
dx2[1:-1] = 0.5*(dx[1:] + dx[:-1])
dx2[1:-1] = 0.5 * (dx[1:] + dx[:-1])
return dx2


def rebalanced_grid(grid, err, base=0.25, num=None, resolution_factor=10, smooth_fact=1.0):
def rebalanced_grid(
grid, err, base=0.25, num=None, resolution_factor=10, smooth_fact=1.0
):
if num is None:
num = grid.size

dx = np.diff(grid)
area_err = 0.5*np.dot(err[1:] + err[:-1], dx) # trapezoidal rule
area_err = 0.5 * np.dot(err[1:] + err[:-1], dx) # trapezoidal rule
dx2 = _avgdiff(grid)

def smooth_err(x):
tot = 0
for i, (gx, e) in enumerate(zip(grid, err)):
fwhm = dx2[i]*smooth_fact
tot += e*np.exp(-(x-gx)**2/(2*(fwhm/2.35482)**2))
fwhm = dx2[i] * smooth_fact
tot += e * np.exp(-((x - gx) ** 2) / (2 * (fwhm / 2.35482) ** 2))
return tot

finegrid = np.zeros((grid.size-1) * resolution_factor + 1)
for i in range(grid.size-1):
finegrid[i*resolution_factor:(i+1)*resolution_factor] = np.linspace(
grid[i], grid[i+1], resolution_factor+1)[:-1]
finegrid[-resolution_factor-1:] = np.linspace(grid[-2], grid[-1], resolution_factor + 1)
smoothed = smooth_err(finegrid) + base*area_err/(grid[-1] - grid[0])
finegrid = np.zeros((grid.size - 1) * resolution_factor + 1)
for i in range(grid.size - 1):
finegrid[i * resolution_factor : (i + 1) * resolution_factor] = np.linspace(
grid[i], grid[i + 1], resolution_factor + 1
)[:-1]
finegrid[-resolution_factor - 1 :] = np.linspace(
grid[-2], grid[-1], resolution_factor + 1
)
smoothed = smooth_err(finegrid) + base * area_err / (grid[-1] - grid[0])
assert np.all(smoothed > 0)
assert np.all(_avgdiff(finegrid) > 0)
interr = np.cumsum(smoothed * _avgdiff(finegrid))
Expand All @@ -43,7 +48,7 @@ def smooth_err(x):


def pre_pruning_mask(grid, rtol=1e-12, atol=0.0):
""" Returns a mask for grid pruning.
"""Returns a mask for grid pruning.
Any grid spacing smaller than ``rtol*gridvalue + atol`` will
be pruned. In general the value on the right is removed unless it is
Expand All @@ -62,10 +67,10 @@ def pre_pruning_mask(grid, rtol=1e-12, atol=0.0):
"""
if np.any(np.diff(grid) < 0):
raise ValueError("grid needs to be monotonic")
limit = grid[-1] - (atol + abs(rtol*grid[-1]))
limit = grid[-1] - (atol + abs(rtol * grid[-1]))
mask = np.empty(grid.size, dtype=np.bool_)
mask[grid.size - 1] = True # rightmost point included
for ridx in range(grid.size-2, -1, -1):
for ridx in range(grid.size - 2, -1, -1):
if grid[ridx] < limit:
mask[ridx] = True
break
Expand All @@ -74,18 +79,18 @@ def pre_pruning_mask(grid, rtol=1e-12, atol=0.0):
else:
raise ValueError("no grid-points left")
mask[0] = True # leftmost point included
limit = grid[0] + abs(rtol*grid[0]) + atol
limit = grid[0] + abs(rtol * grid[0]) + atol
for idx in range(1, ridx):
if grid[idx] < limit:
mask[idx] = False
else:
mask[idx] = True
limit = grid[idx] + abs(rtol*grid[idx]) + atol
limit = grid[idx] + abs(rtol * grid[idx]) + atol
return mask


def combine_grids(grids, **kwargs):
""" Combines multiple grids and prunes them using pre_pruning mask
"""Combines multiple grids and prunes them using pre_pruning mask
Parameters
----------
Expand All @@ -104,7 +109,7 @@ def combine_grids(grids, **kwargs):


def grid_pruning_mask(grid, err, ndrop=None, protect_sparse=None, pow_err=2, pow_dx=2):
""" Returns a mask for grid pruning.
"""Returns a mask for grid pruning.
Parameters
----------
Expand All @@ -126,7 +131,7 @@ def grid_pruning_mask(grid, err, ndrop=None, protect_sparse=None, pow_err=2, pow
protect_sparse = math.ceil(grid.size * 0.25)
dx = _avgdiff(grid)
protected = np.argsort(dx)[-protect_sparse:]
score = err**pow_err * dx**pow_dx
score = err ** pow_err * dx ** pow_dx
importance = np.argsort(score)
drop = []
for considered in importance:
Expand Down
Loading

0 comments on commit bfb1940

Please sign in to comment.