Published
Edited
Feb 4, 2021
Importers
34 stars
// Takes in two arrays of knots and values, and returns a
// Float64Array of the slope at each knot.
natural_cubic_spline = function natural_cubic_spline(knots, values) {
const
t = knots, y = values, n = t.length - 1,
s = new Float64Array(n+1), // horizontal scales
d = new Float64Array(n+1), // diagonal matrix
m = new Float64Array(n+1); // result vector: slope at each knot
for (let i = 0; i < n; i++) {
s[i] = 1 / (t[i+1] - t[i]); // Note s[n] == 0
// Setup pass; at this point 'm' represents the right-hand side
m[i] += (m[i+1] = 3*s[i]*s[i]*(y[i+1] - y[i])); }
d[0] = 0.5 / s[0];
m[0] = d[0] * m[0];
for (let i = 1; i <= n; i++) {
d[i] = 1 / (2*(s[i] + s[i-1]) - s[i-1]*s[i-1]*d[i-1]);
// Forward pass; at this point 'm' represents Um.
m[i] = d[i]*(m[i] - s[i-1]*m[i-1]); }
for (let i = n-1; i >= 0; i--) {
m[i] -= d[i]*s[i]*m[i+1]; }
return m;
}
// Takes in an array of values at knots t = [0, 1, 2, ..., n]
// and returns a Float64Array of the slope at each knot.
natural_cubic_spline_equi = function natural_cubic_spline_equi(values) {
const
y = values, n = y.length - 1,
d = new Float64Array(n+1), // diagonal matrix
m = new Float64Array(n+1); // result vector: slope at each knot
for (let i = 1; i < n; i++)
m[i] = 3*(y[i+1] - y[i-1]);
d[0] = 0.5, m[0] = 1.5*(y[1] - y[0]), m[n] = 3*(y[n] - y[n-1]);
for (let i = 1; i < n; i++)
d[i] = 1 / (4 - d[i-1]), m[i] = (m[i] - m[i-1]) * d[i];
d[n] = 1 / (2 - d[n-1]), m[n] = (m[n] - m[n-1]) * d[n];
for (let i = n-1; i >= 0; i--)
m[i] -= d[i]*m[i+1];
return m;
}
// Takes in two arrays of knots and values, and returns a
// Float64Array of the slope at each knot.
// Assumes that the knots are in the periodic interval [0, 1].
cyclic_cubic_spline = function cyclic_cubic_spline(knots, values) {
const
t = knots, y = values, n = t.length - 1,
s = new Float64Array(n+1), // horizontal scales
d = new Float64Array(n+1), // diagonal matrix
z = new Float64Array(n+1), // correction vector
m = new Float64Array(n+1); // result vector: slope at each knot
for (let i = 0; i < n; i++) {
s[i] = 1 / (t[i+1] - t[i]);
m[i] += (m[i+1] += 3*s[i]*s[i]*(y[i+1] - y[i])); }
s[n] = 1 / (t[0] + 1 - t[n]);
const bn = 3*s[n]*s[n]*(y[0] - y[n]);
m[0] += bn; m[n] += bn;
d[0] = 1 / (2*s[0] + s[n]);
m[0] = d[0] * m[0];
z[0] = d[0] * s[n];
for (let i = 1; i < n; i++) {
d[i] = 1 / (2*(s[i] + s[i-1]) - s[i-1]*s[i-1]*d[i-1]);
m[i] = d[i]*(m[i] - s[i-1]*m[i-1]);
z[i] = d[i]*(0 - s[i-1]*z[i-1]); }
d[n] = 1 / (s[n] + 2*s[n-1] - s[n-1]*s[n-1]*d[n-1]);
m[n] = d[n]*(m[n] - s[n-1]*m[n-1]);
z[n] = d[n]*(s[n] - s[n-1]*z[n-1]);
for (let i = n-1; i >= 0; i--) {
m[i] -= d[i]*s[i]*m[i+1];
z[i] -= d[i]*s[i]*z[i+1]; }
const zscale = (m[0] + m[n]) / (1 + z[0] + z[n]);
for (let i = 0; i <= n; i++) {
m[i] -= zscale * z[i]; }
return m;
}
// Takes a list of knots, values, and derivatives, and constructs
// a cubic polynomial in Bernstein–Bézier basis on each subinterval,
// output as a Float64Array.
draw_hermite_spline = function draw_hermite_spline (knots, values, dvalues) {
const t = knots, y = values, m = dvalues, n = t.length - 1;
const bezpts = new Float64Array(8*n);
for (let i = 0; i < n; i++) {
const h = (t[i+1] - t[i]) * 0.3333333333333333;
bezpts.set([
t[i], y[i], t[i] + h, y[i] + m[i]*h,
t[i+1] - h, y[i+1] - m[i+1]*h, t[i+1], y[i+1]
], 8*i); }
return bezpts;
}
herm_to_cheb3 = function herm_to_cheb3(y0, y1, m0, m1) {
return [
0.5*(y0 + y1 + 0.125*(m0 - m1)),
0.03125*(18*(y1 - y0) - m0 - m1 ),
0.0625*(m1 - m0),
0.03125*(2*(y0 - y1) + m0 + m1)];
}
// Given a list of knots and a target point, finds the subinterval
// in which the target point falls. Points outside the interval are
// considered to be in the first or last interval.
binary_search = function binary_search(knots, target) {
let n = knots.length - 1, low = 0, half;
while ((half = (n/2)|0) > 0) {
low += half * (knots[low + half] <= target);
n -= half; }
return low;
}
// Given an array of knots and an array of values,
// return a function which evaluates a natural cubic spline
// extrapolated by linear functions on either end.
NaturalCubicSpline = function NaturalCubicSpline(knots, values) {
const
n = knots.length - 1,
t = new Float64Array(n+3), // knots
y = new Float64Array(n+3), // values
s = new Float64Array(n+2), // horizontal scales
m = new Float64Array(n+3), // slope at each knot
d = new Float64Array(n+2), // diagonal matrix
coeffs = new Float64Array(4*(n+2)); // Chebyshev coefficients
t.set(knots, 1), y.set(values, 1);
// Natural cubic spline algorithm
for (let i = 1; i < n+1; i++) {
s[i] = 1 / (t[i+1] - t[i]);
m[i] += (m[i+1] = 3*s[i]*s[i]*(y[i+1] - y[i])); }
d[1] = 0.5 / s[1];
m[1] = d[1] * m[1];
for (let i = 2; i <= n+1; i++) {
d[i] = 1 / (2*(s[i] + s[i-1]) - s[i-1]*s[i-1]*d[i-1]);
m[i] = d[i]*(m[i] - s[i-1]*m[i-1]); }
for (let i = n; i >= 1; i--) {
m[i] -= d[i]*s[i]*m[i+1]; }
// Linear extrapolation
t[0] = t[1] - 1; t[n+2] = t[n+1] + 1;
y[0] = y[1] - m[1]; y[n+2] = y[n+1] + m[n+1];
s[0] = s[n+1] = 1;
m[0] = m[1]; m[n+2] = m[n+1];
// Set up Chebyshev coefficients
for (let i = 0; i < n+2; i++) {
const h = t[i+1] - t[i];
const y0 = y[i], y1 = y[i+1], m0 = h*m[i], m1 = h*m[i+1], j = 4*i;
coeffs[j] = 0.5*(y0 + y1 + 0.125*(m0 - m1));
coeffs[j+1] = 0.03125*(18*(y1 - y0) - m0 - m1);
coeffs[j+2] = 0.0625*(m1 - m0);
coeffs[j+3] = 0.03125*(2*(y0 - y1) + m0 + m1); }
return Object.assign(
function eval_spline(u) {
// Binary search
let nn = t.length - 1, low = 0, half;
while ((half = (nn/2)|0) > 0) {
low += half * (t[low + half] <= u);
nn -= half; }
const i = low, j = 4*i;
u = (2*u - t[i] - t[i+1]) * s[i]; // scale to [–1, 1]
// Clenshaw's method.
let b1 = coeffs[j+3], b2 = coeffs[j+2] + 2 * u * b1;
b1 = coeffs[j+1] + 2 * u * b2 - b1;
return coeffs[j] + u * b1 - b2; },
{ knots: t, values: y, dvalues: m, scales: s, coeffs: coeffs });
}
// Takes in two arrays of knots and values, and returns a
// return a function which evaluates a cyclic cubic spline.
// Assumes that the knots are in the periodic interval [0, 1],
// or optionally some other interval.
CyclicCubicSpline = function CyclicCubicSpline(knots, values, domain=[0,1]) {
const
n = knots.length - 1,
scale = 1 / (domain[1] - domain[0]),
t = new Float64Array(n+3), // knots
y = new Float64Array(n+3), // values
s = new Float64Array(n+2), // horizontal scales
m = new Float64Array(n+3), // slope at each knot
d = new Float64Array(n+2), // diagonal matrix
z = new Float64Array(n+2), // correction vector
coeffs = new Float64Array(4*(n+2)); // Chebyshev coefficients
y.set(values, 1);
for (let i = n; i >= 0; i--) {
t[i+1] = (knots[i] - domain[0]) * scale; }
t[0] = t[n+1] - 1; t[n+2] = t[1] + 1;
y[0] = y[n+1]; y[n+2] = y[1];
for (let i = 1; i < n+1; i++) {
s[i] = 1 / (t[i+1] - t[i]);
m[i] += (m[i+1] += 3*s[i]*s[i]*(y[i+1] - y[i])); }
s[n+1] = 1 / (t[1] + 1 - t[n+1]);
const bn = 3*s[n+1]*s[n+1]*(y[1] - y[n+1]);
m[1] += bn; m[n+1] += bn;
d[1] = 1 / (2*s[1] + s[n+1]);
m[1] = d[1] * m[1];
z[1] = d[1] * s[n+1];
for (let i = 2; i < n+1; i++) {
d[i] = 1 / (2*(s[i] + s[i-1]) - s[i-1]*s[i-1]*d[i-1]);
m[i] = d[i]*(m[i] - s[i-1]*m[i-1]);
z[i] = d[i]*(0 - s[i-1]*z[i-1]); }
d[n+1] = 1 / (s[n+1] + 2*s[n] - s[n]*s[n]*d[n]);
m[n+1] = d[n+1]*(m[n+1] - s[n]*m[n]);
z[n+1] = d[n+1]*(s[n+1] - s[n]*z[n]);
for (let i = n; i >= 1; i--) {
m[i] -= d[i]*s[i]*m[i+1];
z[i] -= d[i]*s[i]*z[i+1]; }
const zscale = (m[1] + m[n+1]) / (1 + z[1] + z[n+1]);
for (let i = 1; i <= n+1; i++) {
m[i] -= zscale * z[i]; }
s[0] = s[n+1] = 1 / (t[1] - t[0]);
m[0] = m[n+1]; m[n+2] = m[1];
// Set up Chebyshev coefficients
for (let i = 0; i < n+2; i++) {
const h = t[i+1] - t[i];
const y0 = y[i], y1 = y[i+1], m0 = h*m[i], m1 = h*m[i+1], j = 4*i;
coeffs[j] = 0.5*(y0 + y1 + 0.125*(m0 - m1));
coeffs[j+1] = 0.03125*(18*(y1 - y0) - m0 - m1);
coeffs[j+2] = 0.0625*(m1 - m0);
coeffs[j+3] = 0.03125*(2*(y0 - y1) + m0 + m1); }
return Object.assign(
function eval_spline(u) {
u = ((u - domain[0]) * scale) % 1 // scale interval to [0, 1]
// Binary search
let nn = t.length - 1, low = 0, half;
while ((half = (nn/2)|0) > 0) {
low += half * (t[low + half] <= u);
nn -= half; }
const i = low, j = 4*i;
u = (2*u - t[i] - t[i+1]) * s[i]; // scale subinterval to [–1, 1]
// Clenshaw's method.
let b1 = coeffs[j+3], b2 = coeffs[j+2] + 2 * u * b1;
b1 = coeffs[j+1] + 2 * u * b2 - b1;
return coeffs[j] + u * b1 - b2; },
{ knots: t, values: y, dvalues: m, scales: s, coeffs: coeffs });
}
// Takes in two arrays of knots and values, and returns a
// Float64Array of cubic Bézier segments
parametric_NCS = function parametric_NCS(points, chordal=true) {
const
n = points.length - 1,
y0 = new Float64Array(n+1),
y1 = new Float64Array(n+1),
h = new Float64Array(n), // chord lengths
s = new Float64Array(n+1), // scale = 1/h
d = new Float64Array(n+1), // diagonal matrix
m0 = new Float64Array(n+1), // result vectors: parametric
m1 = new Float64Array(n+1), // slope at each knot
bezpts = new Float64Array(8*n);
for (let i = 0; i <= n; i++) {
[y0[i], y1[i]] = points[i]; }
if (chordal) {
for (let i = 0; i < n; i++) {
const dy0 = y0[i+1] - y0[i];
const dy1 = y1[i+1] - y1[i];
h[i] = Math.sqrt(dy0*dy0 + dy1*dy1)
s[i] = 1 / h[i]; } } // Note s[n] == 0
else {
h.fill(1, 0, n);
s.fill(1, 0, n); }
for (let i = 0; i < n; i++) {
// Setup pass; at this point 'm' represents the right-hand side
m0[i] += (m0[i+1] = 3*s[i]*s[i]*(y0[i+1] - y0[i]));
m1[i] += (m1[i+1] = 3*s[i]*s[i]*(y1[i+1] - y1[i])); }
d[0] = 0.5 / s[0];
m0[0] = d[0] * m0[0];
m1[0] = d[0] * m1[0];
for (let i = 1; i <= n; i++) {
d[i] = 1 / (2*(s[i] + s[i-1]) - s[i-1]*s[i-1]*d[i-1]);
// Forward pass; at this point 'm' represents Um.
m0[i] = d[i]*(m0[i] - s[i-1]*m0[i-1]);
m1[i] = d[i]*(m1[i] - s[i-1]*m1[i-1]); }
for (let i = n-1; i >= 0; i--) {
m0[i] -= d[i]*s[i]*m0[i+1];
m1[i] -= d[i]*s[i]*m1[i+1]; }
for (let i = 0; i < n; i++) {
const hi = h[i] * 0.3333333333333333;
bezpts.set([
y0[i], y1[i], y0[i] + hi*m0[i], y1[i] + hi*m1[i],
y0[i+1] - hi*m0[i+1], y1[i+1] - hi*m1[i+1], y0[i+1], y1[i+1],
], 8*i); }
return bezpts;
}
d3 = require("d3@5")
import {bezeval, bezpts_to_svgpath} from '@jrus/bezplot'
import {evaluate as chebeval} from '@jrus/cheb'
import {responsive_katex, math_css, $ as $tex} with {$css as $css} from "@jrus/misc"
$ = $tex.options({
macros: { "\\bt": "\\bar{t}" }
});
$$ = $tex.options({
displayMode: true,
macros: { "\\bt": "\\bar{t}" }
});
// version from D3
function controlPoints(x) {
var i, n = x.length - 1, m, a = new Array(n), b = new Array(n), r = new Array(n);
a[0] = 0, b[0] = 2, r[0] = x[0] + 2 * x[1];
for (i = 1; i < n - 1; i++) a[i] = 1, b[i] = 4, r[i] = 4 * x[i] + 2 * x[i + 1];
a[n - 1] = 2, b[n - 1] = 7, r[n - 1] = 8 * x[n - 1] + x[n];
for (i = 1; i < n; ++i) m = a[i] / b[i - 1], b[i] -= m, r[i] -= m * r[i - 1];
a[n - 1] = r[n - 1] / b[n - 1];
for (i = n - 2; i >= 0; --i) a[i] = (r[i] - a[i + 1]) / b[i];
b[n - 1] = (x[n] + a[n - 1]) / 2;
for (i = 0; i < n - 1; ++i) b[i] = 2 * x[i + 1] - a[i + 1];
return [a, b];
}
// Alternative to D3 version
control_points = function control_points(y) {
var i, n = y.length - 1, d = new Array(n+1), m = new Array(n+1);
for (i = 1; i < n; i++) m[i] = (y[i+1] - y[i-1]);
m[0] = (d[0] = .5)*(y[1] - y[0]), m[n] = (y[n] - y[n-1]);
for (i = 1; i < n; i++) m[i] = (m[i] - m[i-1]) * (d[i] = 1 / (4 - d[i-1]));
m[n] = (m[n] - m[n-1]) * (d[n] = 1 / (2 - d[n-1]));
for (i = n-1; i >= 0; i--) m[i] -= d[i]*m[i+1];
for (i = 0; i < n; i++) d[i] = y[i+1] - m[i+1], m[i] = y[i] + m[i];
d.pop(), m.pop();
return [m, d];
}
One platform to build and deploy the best data apps
Experiment and prototype by building visualizations in live JavaScript notebooks. Collaborate with your team and decide which concepts to build out.
Use Observable Framework to build data apps locally. Use data loaders to build in any language or library, including Python, SQL, and R.
Seamlessly deploy to Observable. Test before you ship, use automatic deploy-on-commit, and ensure your projects are always up-to-date.