# Path: finitediff/util.py
# def interpolate_ahead(x, y, n, direction="fw"):
#     if not np.all(np.diff(x) > 0):
#         raise ValueError("x not strictly monotonic.")
# 
#     if direction == "both":
#         y1, slc1 = interpolate_ahead(x, y, n, "fw")
#         y2, slc2 = interpolate_ahead(x, y, n, "bw")
#         y = np.zeros_like(x)
#         y[slc1] = y1
#         y[slc2] += y2
#         y[slc1.start : slc2.stop] /= 2
#         return y, slice(None, None)
#     elif direction == "fw":
#         forward = True
#     elif direction == "bw":
#         forward = False
#     else:
#         raise ValueError("Unknown direction: %s" % direction)
# 
#     rev = slice(None, None, 1) if forward else slice(None, None, -1)
#     values = []
#     for idx, (xv, yv) in enumerate(zip(x[rev][n:], y[rev][n:])):
#         _x = np.ascontiguousarray(x[rev][idx : idx + n], dtype=np.float64)
#         _y = np.ascontiguousarray(y[rev][idx : idx + n], dtype=np.float64)
#         _v = np.array([xv])
#         values.append(
#             interpolate_by_finite_diff(_x, _y, _v, maxorder=0, ntail=n, nhead=0)
#         )
#     return np.array(values[rev]).squeeze(), slice(n, None) if forward else slice(
#         None, -n
#     )
# 
# def avg_stddev(arr, w):
#     """Calculates the average and standard deviation.
# 
#     Parameters
#     ----------
#     arr : array_like
#         Values.
#     w : array_like
#         Weights.
# 
#     Returns
#     -------
#     tuple of 2 floats (average & standard deviation)
# 
#     """
#     avg, wsum = np.average(arr, weights=w, returned=True)
#     res = arr - avg
#     stddev = np.sqrt(np.sum(np.dot(w, np.square(res)) / (res.size - 1) / wsum))
#     return avg, stddev


import numpy as np
from ..util import interpolate_ahead, avg_stddev
    dtg = np.diff(tg)
    err = np.zeros(y.size)
    for d in ("fw", "bw"):
        est, slc = interpolate_ahead(tg, y, ntrail, d)
        start = (ntrail - 1) if d == "fw" else 0
        stop = -(ntrail - 1) if d == "bw" else None
        err[slc] += np.abs(y[slc] - est) / dtg[start:stop] * dy[start:stop]
    imax = np.argsort(err)[-consider:][::-1]
    return [(tg[m], err[m]) for m in imax]


def pool_discontinuity_approx(loc_res, consistency_criterion=10):
    points = np.array(loc_res)
    w1 = np.abs(points[:, 1])
    return avg_stddev(points[:, 0], w1)


def grid_error(grid, y, ntrail=2):
    """Estimates error at each grid point from extrapolation.

    Extrapolates from left ('fw') or from right ('bw').

    Parameters
    ----------
    grid : array
    y : array
    ntrail : int
        Number of points to include in the look-ahead extrapolation.

    """

# Path: finitediff/util.py
# def interpolate_ahead(x, y, n, direction="fw"):
#     if not np.all(np.diff(x) > 0):
#         raise ValueError("x not strictly monotonic.")
# 
#     if direction == "both":
#         y1, slc1 = interpolate_ahead(x, y, n, "fw")
#         y2, slc2 = interpolate_ahead(x, y, n, "bw")
#         y = np.zeros_like(x)
#         y[slc1] = y1
#         y[slc2] += y2
#         y[slc1.start : slc2.stop] /= 2
#         return y, slice(None, None)
#     elif direction == "fw":
#         forward = True
#     elif direction == "bw":
#         forward = False
#     else:
#         raise ValueError("Unknown direction: %s" % direction)
# 
#     rev = slice(None, None, 1) if forward else slice(None, None, -1)
#     values = []
#     for idx, (xv, yv) in enumerate(zip(x[rev][n:], y[rev][n:])):
#         _x = np.ascontiguousarray(x[rev][idx : idx + n], dtype=np.float64)
#         _y = np.ascontiguousarray(y[rev][idx : idx + n], dtype=np.float64)
#         _v = np.array([xv])
#         values.append(
#             interpolate_by_finite_diff(_x, _y, _v, maxorder=0, ntail=n, nhead=0)
#         )
#     return np.array(values[rev]).squeeze(), slice(n, None) if forward else slice(
#         None, -n
#     )
# 
# def avg_stddev(arr, w):
#     """Calculates the average and standard deviation.
# 
#     Parameters
#     ----------
#     arr : array_like
#         Values.
#     w : array_like
#         Weights.
# 
#     Returns
#     -------
#     tuple of 2 floats (average & standard deviation)
# 
#     """
#     avg, wsum = np.average(arr, weights=w, returned=True)
#     res = arr - avg
#     stddev = np.sqrt(np.sum(np.dot(w, np.square(res)) / (res.size - 1) / wsum))
#     return avg, stddev


# Path: finitediff/grid/util.py
import numpy as np
from ..util import interpolate_ahead, avg_stddev

    dtg = np.diff(tg)
    err = np.zeros(y.size)
    for d in ("fw", "bw"):
        est, slc = interpolate_ahead(tg, y, ntrail, d)
        start = (ntrail - 1) if d == "fw" else 0
        stop = -(ntrail - 1) if d == "bw" else None
        err[slc] += np.abs(y[slc] - est) / dtg[start:stop] * dy[start:stop]
    imax = np.argsort(err)[-consider:][::-1]
    return [(tg[m], err[m]) for m in imax]


def pool_discontinuity_approx(loc_res, consistency_criterion=10):
    points = np.array(loc_res)
    w1 = np.abs(points[:, 1])
    return avg_stddev(points[:, 0], w1)


def grid_error(grid, y, ntrail=2):
    """Estimates error at each grid point from extrapolation.

    Extrapolates from left ('fw') or from right ('bw').

    Parameters
    ----------
    grid : array
    y : array
    ntrail : int
        Number of points to include in the look-ahead extrapolation.

    """

    est, _ = interpolate_ahead(grid, y, ntrail, direction="both")