In [ ]:
import numpy as np
import matplotlib.pyplot as plt
from finitediff.grid import adapted_grid, plot_convergence
from finitediff.grid.tests._common import g, g2
%matplotlib inline

In [ ]:
plot_convergence('grid_additions', [(32,), (16, 16), (8, 8, 8, 8)], g)

In [ ]:
plot_convergence('grid_additions', [(64,), (32, 32), (48, 12, 4)], g)

In [ ]:
blr = ((.3, .1), (.3, .1))
plot_convergence('grid_additions', [(64,), (32, 32), (48, 12, 4)], g, blurs=blr)

In [ ]:
plot_convergence('grid_additions', [(64,), (32, 32), (32, 24, 4, 4)], g, blurs=blr)

In [ ]:
plot_convergence('grid_additions', [(32,)*2, (16,)*4, (8,)*8], g, blurs=blr)

In [ ]:
plot_convergence('grid_additions', [(32,), (16,)*2, (8,)*4, (4,)*8], g, blurs=blr)

In [ ]:
plot_convergence('grid_additions', [(48,), (24,)*2, (12,)*4, (16,)*3], g)

In [ ]:


In [ ]:
plot_convergence('grid_additions', [(32,), (16,)*2, (8,)*4, (4,)*8], g2)

In [ ]:
plot_convergence('grid_additions', [(32,), (16,)*2, (8,)*4, (4,)*8], g2, ntrail=3)

In [ ]:
plot_convergence('grid_additions', [(32, 32), (16,)*4], g2)

In [ ]:
plot_convergence('grid_additions', [(32, 32), (16,)*4], g2, ntrail=3, blurs=blr)

In [ ]:
plot_convergence('grid_additions', [(64,)*2, (16,)*8], g2, ntrail=3, blurs=blr)

In [ ]:
def gs(x):
    return [(a, b) for a, b in zip(g(x), g2(x))]
plot_convergence('grid_additions', [(64,)*2, (16,)*8], gs, ntrail=3, blurs=blr,
                 metric=lambda x: x[0]+x[1])

In [ ]:
plot_convergence('ntrail', [2, 3, 4], g)

In [ ]:
plot_convergence('blurs', [((.5,), (.5,)), ((.5, .25), (.5, .25)), ((.125,), (.125,))], g)

In [ ]:
plot_convergence('blurs', [((.5,), (.5,)), ((.5, .25), (.5, .25)), ((.125,), (.125,))],
                 g, ntrail=4)

Grid refinement


In [ ]:
from scipy.interpolate import BPoly
from finitediff.grid import refine_grid

knot_arrays = [[0, 1, 2, 3, 7, 15], [0, 2.5, 5, 10, 11, 12, 13, 14, 15]]
knot_values = [[0., 0, 1, 1, 0, 0], [0., 0, 1, 1, 0, 0, 1, 1, 0]]
bpolys = []
for x, y in zip(knot_arrays, knot_values):
    bpolys.append(BPoly.from_derivatives(x, list(zip(y, *([np.zeros(len(y))]*2)))))
    
def plot_interpol(polys, xplt=None, ax=None):
    if ax is None:
        fig, ax = plt.subplots(1, 1, figsize=(16, 4))
    if xplt is None:
        xplt = np.linspace(0, 15, 500)
    for p in polys:
        ax.plot(xplt, p(xplt))
    for xx in knot_arrays:
        ax.vlines(xx, .1, .4, transform=ax.get_xaxis_transform(), color='k', linestyle='--', linewidth=0.5)
    ax.vlines(xplt, .6, .9, transform=ax.get_xaxis_transform(), lw=.5, alpha=.5, linestyle=':')
    return ax

In [ ]:
def demo_refine(cbs, ax=None, **kwargs):
    grid = np.linspace(0, 15, 160)
    for cb in cbs:
        grid, err = refine_grid(grid, cb, **kwargs)
    return grid, plot_interpol(bpolys, grid, ax=ax)

In [ ]:
fig, axes = plt.subplots(8, 1, figsize=(16, 32))
demo_refine(bpolys, grid_additions=[80]*2, ax=axes[0])
demo_refine(bpolys, grid_additions=[60]*3, ax=axes[1])
demo_refine(bpolys, grid_additions=[160], ax=axes[2])
demo_refine(bpolys[::-1], grid_additions=[80]*2, ax=axes[3])
demo_refine(bpolys, grid_additions=[80]*2, blurs=[(.5,)]*2, ax=axes[4])
demo_refine(bpolys, grid_additions=[80]*2, blurs=[(.5, .25, .125)]*2, ax=axes[5])
demo_refine(bpolys, grid_additions=[80]*2, ntrail=3, ax=axes[6])
grid, ax = demo_refine(bpolys + bpolys[::-1], grid_additions=[80], ntrail=3, ax=axes[7])

Grid pruning


In [ ]:
from finitediff.grid import grid_pruning_mask, grid_error
pruning_mask = grid_pruning_mask(grid, np.sum([grid_error(grid, cb(grid)) for cb in bpolys]))
new_grid = grid[pruning_mask]
plot_interpol(bpolys, new_grid)
grid.size, new_grid.size

Rebalancing grid


In [ ]:
%load_ext autoreload

In [ ]:
%autoreload 2

In [ ]:
from finitediff.grid import rebalanced_grid

In [ ]:
def mk_grid(grid, cbs, niter=1, pruning_factor=0.0, **kwargs):
    def _calc_err(grd):
        err = np.zeros_like(grid)
        for cb in cbs:
            err += np.abs(grid_error(grid, cb(grid)))
        return err
    grids, errs = [grid], []
    for _ in range(niter):
        errs.append(_calc_err(grid))
        grid = rebalanced_grid(grid, errs[-1], **kwargs)
        grids.append(grid)
    errs.append(_calc_err(grid))
    return grids, errs

In [ ]:
ngrids = 8
grids, errs = mk_grid(np.linspace(0, 15, 160), bpolys, niter=ngrids, base=0.2, smooth_fact=10.0)
fig, axes = plt.subplots(ngrids, 1, figsize=(16, 16))
for grid, ax in zip(grids, axes):
    plot_interpol(bpolys, grid, ax=ax)

In [ ]:
from finitediff.grid import combine_grids
combined_a = combine_grids(grids[:2], atol=1e-12)
combined_b = combine_grids(grids[2:4], atol=1e-12)
combined_c = combine_grids(grids[4:6], atol=1e-12)
combined_d = combine_grids(grids[6:8], atol=1e-12)

In [ ]:
fig, axes = plt.subplots(5, 1, figsize=(16, 16))
plot_interpol(bpolys, combined_a, ax=axes[0])
plot_interpol(bpolys, combined_b, ax=axes[1])
plot_interpol(bpolys, combined_c, ax=axes[2])
plot_interpol(bpolys, combined_d, ax=axes[3])
d_mask = grid_pruning_mask(combined_d, np.sum([grid_error(combined_d, cb(combined_d)) for cb in bpolys]))
plot_interpol(bpolys, combined_d[d_mask], ax=axes[4])

In [ ]:
def combine_and_prune(grids, errs, cbs):
    supergrid = combine_grids([grid[grid_pruning_mask(grid, err)] for grid, err in zip(grids, errs)], atol=1e-12)
    assert np.all(np.diff(supergrid) > 0)
    err = np.zeros_like(supergrid)
    for cb in cbs:
        err += np.abs(grid_error(supergrid, cb(supergrid)))
    ma = grid_pruning_mask(supergrid, err, supergrid.size - grids[0].size, grids[0].size//2)
    return supergrid[ma]

In [ ]:
plot_interpol(bpolys, combine_and_prune(grids[:4], errs[:4], bpolys))

In [ ]:
from itertools import chain
def cookbook_grid(cbs, n=160):
    ngrids = 4
    base_grids, base_errs = mk_grid(np.linspace(0, 15, n//4), cbs, niter=ngrids, base=0.2, smooth_fact=10.0)
    supergrid = combine_and_prune(base_grids, base_errs, cbs)
    for cb in chain(cbs, reversed(cbs)):
        supergrid, _errs = refine_grid(supergrid, cb, grid_additions=[n//4])
    prune_mask = grid_pruning_mask(supergrid, np.sum([grid_error(supergrid, cb(supergrid)) for cb in cbs]))
    return supergrid[prune_mask]

In [ ]:
plot_interpol(bpolys, cookbook_grid(bpolys))

In [ ]: