Interpolate x, y using a quadratic algorithm
    https://en.wikipedia.org/wiki/Spline_(mathematics)
547848
def quadratic_interpolate(x, y, precision=250, **kwargs):
    
    n = len(x) - 1
    delta_x = [x2 - x1 for x1, x2 in zip(x, x[1:])]
    delta_y = [y2 - y1 for y1, y2 in zip(y, y[1:])]
    slope = [delta_y[i] / delta_x[i] if delta_x[i] else 1 for i in range(n)]

    # Quadratic spline: a + bx + cx²
    a = y
    b = [0] * (n + 1)
    c = [0] * (n + 1)

    for i in range(1, n):
        b[i] = 2 * slope[i - 1] - b[i - 1]

    c = [(slope[i] - b[i]) / delta_x[i] if delta_x[i] else 0 for i in range(n)]

    for i in range(n + 1):
        yield x[i], a[i]
        if i == n or delta_x[i] == 0:
            continue
        for s in range(1, precision):
            X = s * delta_x[i] / precision
            X2 = X * X
            yield x[i] + X, a[i] + b[i] * X + c[i] * X2
457881
function interpolate(t) {
    var A = Math.sin(d - (t *= d)) / sd,
        B = Math.sin(t) / sd,
        x = A * cy0 * cx0 + B * cy1 * cx1,
        y = A * cy0 * sx0 + B * cy1 * sx1,
        z = A * sy0       + B * sy1;
    return [
      Math.atan2(y, x) / d3_geo_radians,
      Math.atan2(z, Math.sqrt(x * x + y * y)) / d3_geo_radians
    ];
  }
507274
def interpolate_grid(self, in_lon, in_lat):
        
        out_data = np.zeros((self.data.shape[0], in_lon.shape[0], in_lon.shape[1]))
        for d in range(self.data.shape[0]):
            print("Loading ", d, self.variable, self.start_date)
            if self.data[d].max() > -999:
                step = self.data[d]
                step[step < 0] = 0
                if self.lat[-1] < self.lat[0]:
                    spline = RectBivariateSpline(self.lat[::-1], self.lon, step[::-1], kx=3, ky=3)
                else:
                    spline = RectBivariateSpline(self.lat, self.lon, step, kx=3, ky=3)
                print("Evaluating", d, self.variable, self.start_date)
                flat_data = spline.ev(in_lat.ravel(), in_lon.ravel())
                out_data[d] = flat_data.reshape(in_lon.shape)
                del spline
            else:
                print(d, " is missing")
                out_data[d] = -9999
        return out_data
142451
def smooth_differentiation(x, y, weigths=None, order=5, smoothness=3, derivation=1):
    '''
    '''
    if (len(x) != len(y)):
        raise ValueError("x, y must have the same length")
    f = splrep(x, y, w=weigths, k=order, s=smoothness)  # spline function
    return splev(x, f, der=derivation)
1078254
function xinterp (points) {
        // Scale the data.
        var f = [xscale(config.xlim[1]) - xscale(config.xlim[0]),
                 yscale(config.ylim[1]) - yscale(config.ylim[0])],
            z = [xscale(config.xlim[0]),
                 yscale(config.ylim[0])],
            scaled = points.map(function (p) {
                return [(p[0] - z[0]) / f[0], (p[1] - z[1]) / f[1]];
            });

        // Compute the distance along the path using a map-reduce.
        var dists = scaled.map(function (d, i) {
            if (i == 0) return 0.0;
            var dx = d[0] - scaled[i - 1][0],
                dy = d[1] - scaled[i - 1][1];
            return Math.sqrt(dx * dx + dy * dy);
        }),
            dist = dists.reduce(function (curr, d) { return d + curr; }, 0.0);

        // Choose the number of interpolation points based on this distance.
        var N = Math.round(200 * dist);

        // Re-sample the line.
        var resampled = [];
        dists.map(function (d, i) {
            if (i == 0) return;
            var n = Math.max(3, Math.round(d / dist * N)),
                spline = d3.interpolate(scaled[i - 1][1], scaled[i][1]),
                delta = (scaled[i][0] - scaled[i - 1][0]) / (n - 1);
            for (var j = 0, x = scaled[i - 1][0]; j < n; ++j, x += delta)
                resampled.push([x, spline(j / (n - 1))]);
        });

        // Compute the gradients.
        var gradients = resampled.map(function (a, i, d) {
            if (i == 0) return [d[1][0] - d[0][0], d[1][1] - d[0][1]];
            if (i == resampled.length - 1)
                return [d[i][0] - d[i - 1][0], d[i][1] - d[i - 1][1]];
            return [0.5 * (d[i + 1][0] - d[i - 1][0]),
                    0.5 * (d[i + 1][1] - d[i - 1][1])];
        });

        // Normalize the gradient vectors to be unit vectors.
        gradients = gradients.map(function (d) {
            var len = Math.sqrt(d[0] * d[0] + d[1] * d[1]);
            return [d[0] / len, d[1] / len];
        });

        // Generate some perturbations.
        var perturbations = smooth(resampled.map(d3.random.normal()), 3);

        // Add in the perturbations and re-scale the re-sampled curve.
        var result = resampled.map(function (d, i) {
            var p = perturbations[i],
                g = gradients[i];
            return [(d[0] + config.magnitude * g[1] * p) * f[0] + z[0],
                    (d[1] - config.magnitude * g[0] * p) * f[1] + z[1]];
        });

        return result.join("L");
    }
919234
def _interpolated_template(self, templateid):
        """"""
        phase, y = self._get_template_by_id(templateid)

        # double-check that phase ranges from 0 to 1
        assert phase.min() >= 0
        assert phase.max() <= 1

        # at the start and end points, we need to add ~5 points to make sure
        # the spline & derivatives wrap appropriately
        phase = np.concatenate([phase[-5:] - 1, phase, phase[:5] + 1])
        y = np.concatenate([y[-5:], y, y[:5]])

        # Univariate spline allows for derivatives; use this!
        return UnivariateSpline(phase, y, s=0, k=5)
1226304
def interpolate(p, axis_values, pixelgrid, order=1, mode='constant', cval=0.0):
    
    # convert requested parameter combination into a coordinate
    p_ = np.array([np.searchsorted(av_,val) for av_, val in zip(axis_values,p)])
    lowervals_stepsize = np.array([[av_[p__-1], av_[p__]-av_[p__-1]] \
                         for av_, p__ in zip(axis_values,p_)])
    p_coord = (p-lowervals_stepsize[:,0])/lowervals_stepsize[:,1] + p_-1
    
    # interpolate
    if order > 1:
        prefilter = False
    else:
        prefilter = False

    return [ndimage.map_coordinates(pixelgrid[...,i],p_coord, order=order,
                                    prefilter=prefilter, mode=mode, cval=cval) \
                for i in range(np.shape(pixelgrid)[-1])]
829410
function generateBezier (mX1, mY1, mX2, mY2) {
        var NEWTON_ITERATIONS = 4,
            NEWTON_MIN_SLOPE = 0.001,
            SUBDIVISION_PRECISION = 0.0000001,
            SUBDIVISION_MAX_ITERATIONS = 10,
            kSplineTableSize = 11,
            kSampleStepSize = 1.0 / (kSplineTableSize - 1.0),
            float32ArraySupported = "Float32Array" in window;

        /* Must contain four arguments. */
        if (arguments.length !== 4) {
            return false;
        }

        /* Arguments must be numbers. */
        for (var i = 0; i < 4; ++i) {
            if (typeof arguments[i] !== "number" || isNaN(arguments[i]) || !isFinite(arguments[i])) {
                return false;
            }
        }

        /* X values must be in the [0, 1] range. */
        mX1 = Math.min(mX1, 1);
        mX2 = Math.min(mX2, 1);
        mX1 = Math.max(mX1, 0);
        mX2 = Math.max(mX2, 0);

        var mSampleValues = float32ArraySupported ? new Float32Array(kSplineTableSize) : new Array(kSplineTableSize);

        function A (aA1, aA2) { return 1.0 - 3.0 * aA2 + 3.0 * aA1; }
        function B (aA1, aA2) { return 3.0 * aA2 - 6.0 * aA1; }
        function C (aA1)      { return 3.0 * aA1; }

        function calcBezier (aT, aA1, aA2) {
            return ((A(aA1, aA2)*aT + B(aA1, aA2))*aT + C(aA1))*aT;
        }

        function getSlope (aT, aA1, aA2) {
            return 3.0 * A(aA1, aA2)*aT*aT + 2.0 * B(aA1, aA2) * aT + C(aA1);
        }

        function newtonRaphsonIterate (aX, aGuessT) {
            for (var i = 0; i < NEWTON_ITERATIONS; ++i) {
                var currentSlope = getSlope(aGuessT, mX1, mX2);

                if (currentSlope === 0.0) return aGuessT;

                var currentX = calcBezier(aGuessT, mX1, mX2) - aX;
                aGuessT -= currentX / currentSlope;
            }

            return aGuessT;
        }

        function calcSampleValues () {
            for (var i = 0; i < kSplineTableSize; ++i) {
                mSampleValues[i] = calcBezier(i * kSampleStepSize, mX1, mX2);
            }
        }

        function binarySubdivide (aX, aA, aB) {
            var currentX, currentT, i = 0;

            do {
                currentT = aA + (aB - aA) / 2.0;
                currentX = calcBezier(currentT, mX1, mX2) - aX;
                if (currentX > 0.0) {
                  aB = currentT;
                } else {
                  aA = currentT;
                }
            } while (Math.abs(currentX) > SUBDIVISION_PRECISION && ++i < SUBDIVISION_MAX_ITERATIONS);

            return currentT;
        }

        function getTForX (aX) {
            var intervalStart = 0.0,
                currentSample = 1,
                lastSample = kSplineTableSize - 1;

            for (; currentSample != lastSample && mSampleValues[currentSample] <= aX; ++currentSample) {
                intervalStart += kSampleStepSize;
            }

            --currentSample;

            var dist = (aX - mSampleValues[currentSample]) / (mSampleValues[currentSample+1] - mSampleValues[currentSample]),
                guessForT = intervalStart + dist * kSampleStepSize,
                initialSlope = getSlope(guessForT, mX1, mX2);

            if (initialSlope >= NEWTON_MIN_SLOPE) {
                return newtonRaphsonIterate(aX, guessForT);
            } else if (initialSlope == 0.0) {
                return guessForT;
            } else {
                return binarySubdivide(aX, intervalStart, intervalStart + kSampleStepSize);
            }
        }

        var _precomputed = false;

        function precompute() {
            _precomputed = true;
            if (mX1 != mY1 || mX2 != mY2) calcSampleValues();
        }

        var f = function (aX) {
            if (!_precomputed) precompute();
            if (mX1 === mY1 && mX2 === mY2) return aX;
            if (aX === 0) return 0;
            if (aX === 1) return 1;

            return calcBezier(getTForX(aX), mY1, mY2);
        };

        f.getControlPoints = function() { return [{ x: mX1, y: mY1 }, { x: mX2, y: mY2 }]; };

        var str = "generateBezier(" + [mX1, mY1, mX2, mY2] + ")";
        f.toString = function () { return str; };

        return f;
    }
734026
def spline(x,knots,p,i=0.0):
    """"""
    assert(p+1<len(knots))
    return np.array([N(float(u),float(i),float(p),knots) for u in x])
21535
def interpolate_2d(values, method='pad', axis=0, limit=None, fill_value=None,
                   dtype=None):
    

    transf = (lambda x: x) if axis == 0 else (lambda x: x.T)

    # reshape a 1 dim if needed
    ndim = values.ndim
    if values.ndim == 1:
        if axis != 0:  # pragma: no cover
            raise AssertionError("cannot interpolate on a ndim == 1 with "
                                 "axis != 0")
        values = values.reshape(tuple((1,) + values.shape))

    if fill_value is None:
        mask = None
    else:  # todo create faster fill func without masking
        mask = mask_missing(transf(values), fill_value)

    method = clean_fill_method(method)
    if method == 'pad':
        values = transf(pad_2d(
            transf(values), limit=limit, mask=mask, dtype=dtype))
    else:
        values = transf(backfill_2d(
            transf(values), limit=limit, mask=mask, dtype=dtype))

    # reshape back
    if ndim == 1:
        values = values[0]

    return values
536144
def lagrange_interpolate(x, y, precision=250, **kwargs):
    
    n = len(x) - 1
    delta_x = [x2 - x1 for x1, x2 in zip(x, x[1:])]
    for i in range(n + 1):
        yield x[i], y[i]
        if i == n or delta_x[i] == 0:
            continue

        for s in range(1, precision):
            X = x[i] + s * delta_x[i] / precision
            s = 0
            for k in range(n + 1):
                p = 1
                for m in range(n + 1):
                    if m == k:
                        continue
                    if x[k] - x[m]:
                        p *= (X - x[m]) / (x[k] - x[m])
                s += y[k] * p
            yield X, s
187160
function quadratic2cubicBezier (a, b) {
  return [
    [
      (a / 3 + (a + b) / 3 - a) / (b - a),
      (a * a / 3 + a * b * 2 / 3 - a * a) / (b * b - a * a)
    ], [
      (b / 3 + (a + b) / 3 - a) / (b - a),
      (b * b / 3 + a * b * 2 / 3 - a * a) / (b * b - a * a)
    ]
  ]
}
401382
private double cubicInterpolate(final double Y0, final double Y1, final double Y2, final double Y3, final double MU) {
        final double A0;
        final double A1;
        final double A2;
        final double A3;
        final double MU2;

        MU2 = MU * MU;
        A0 = Y3 - Y2 - Y0 + Y1;
        A1 = Y0 - Y1 - A0;
        A2 = Y2 - Y0;
        A3 = Y1;

        return ((A0 * MU * MU2) + (A1 * MU2) + (A2 * MU) + A3);
    }
369439
private static double interpolateUsingXAndYTables(final double[] xArr, final double[] yArr,
      final int offset, final double x) {
    return cubicInterpolate(xArr[offset], yArr[offset], xArr[offset + 1], yArr[offset + 1],
        xArr[offset + 2], yArr[offset + 2], xArr[offset + 3], yArr[offset + 3], x);
  }
1341480
def _interpolate(self, kind='linear'):
        """"""
        # drop any existing non-dimension coordinates along the resampling
        # dimension
        dummy = self._obj.copy()
        for k, v in self._obj.coords.items():
            if k != self._dim and self._dim in v.dims:
                dummy = dummy.drop(k)
        return dummy.interp(assume_sorted=True, method=kind,
                            kwargs={'bounds_error': False},
                            **{self._dim: self._full_index})
101778
def cubic_interpolate(x, y, precision=250, **kwargs):
    
    n = len(x) - 1
    # Spline equation is a + bx + cx² + dx³
    # ie: Spline part i equation is a[i] + b[i]x + c[i]x² + d[i]x³
    a = y
    b = [0] * (n + 1)
    c = [0] * (n + 1)
    d = [0] * (n + 1)
    m = [0] * (n + 1)
    z = [0] * (n + 1)

    h = [x2 - x1 for x1, x2 in zip(x, x[1:])]
    k = [a2 - a1 for a1, a2 in zip(a, a[1:])]
    g = [k[i] / h[i] if h[i] else 1 for i in range(n)]

    for i in range(1, n):
        j = i - 1
        l = 1 / (2 * (x[i + 1] - x[j]) - h[j] * m[j]) if x[i + 1] - x[j] else 0
        m[i] = h[i] * l
        z[i] = (3 * (g[i] - g[j]) - h[j] * z[j]) * l

    for j in reversed(range(n)):
        if h[j] == 0:
            continue
        c[j] = z[j] - (m[j] * c[j + 1])
        b[j] = g[j] - (h[j] * (c[j + 1] + 2 * c[j])) / 3
        d[j] = (c[j + 1] - c[j]) / (3 * h[j])

    for i in range(n + 1):
        yield x[i], a[i]
        if i == n or h[i] == 0:
            continue
        for s in range(1, precision):
            X = s * h[i] / precision
            X2 = X * X
            X3 = X2 * X
            yield x[i] + X, a[i] + b[i] * X + c[i] * X2 + d[i] * X3
441477
def smooth(polylines):
    
    for c in polylines:
        if len(c) < 9:
            # smoothing wouldn't make sense here
            continue
        x = c[:, 0]
        y = c[:, 1]

        t = np.arange(x.shape[0], dtype=float)
        t /= t[-1]
        x = UnivariateSpline(t, x)(t)
        y = UnivariateSpline(t, y)(t)
        c[:, 0] = x
        c[:, 1] = y