Spline

spline.rb

# Based on this video: https://www.youtube.com/watch?v=dxvmafuP9Wk
# Screenshot:          https://twitter.com/josh_cheek/status/1049913914585165826
# Animation:           https://twitter.com/josh_cheek/status/1050054519911063552
require 'graphics'
require 'matrix'

class Spline < Graphics::Simulation
  def initialize(points, original_fn)
    super 800, 600
    @points      = points
    @original_fn = original_fn
    (@xmin, @xmax), (@ymin, @ymax) = points.transpose.map &:minmax
    @∆x          = xmax - xmin
    @∆y          = ymax - ymin
    @xmid        = xmin + ∆x/2
    @ymid        = ymin + ∆y/2
    @margin      = 100 # px
    @num_samples = 300
    @bounds_to_coefficients = solve
    register_color :orange, 0xFF, 0x88, 0x00
  end

  def draw(*)
    draw_axes :white
    draw_fn "Original", 0, :yellow, &original_fn
    draw_fn "Spline",   1, :red do |x|
      a, b, c = coefficients_for x
      a*x*x + b*x + c
    end
    draw_points :orange, :white
  end

  private

  attr_reader :points, :num_samples, :original_fn, :margin
  attr_reader :∆x, :∆y
  attr_reader :xmin, :xmid, :xmax
  attr_reader :ymin, :ymid, :ymax

  def draw_axes(color)
    vline project_x(xmid), color, 0, h
    hline project_y(ymid), color, 0, w
  end

  def draw_fn(label, label_index, color, &fn)
    num_samples.succ.times.map do |i|
      x = xmin + ∆x*i/num_samples.to_f
      [x, fn[x]]
    end.each_cons(2) do |p1, p2|
      line *project(p1), *project(p2), color
    end
    text label, margin, h-margin-font.height*label_index, color
  end

  def draw_points(fill, stroke)
    points.each { |point| circle *project(point), 5, fill, true }
    points.each { |point| circle *project(point), 5, stroke, false }
  end

  def solve
    # I don't understand how linear algebra works, but I looked up how to solve
    # a system of equations and it told me to do this.
    #
    # Our equations are quadratic, ie:
    #
    #    ax^2 + bx + c = y
    #
    # And we have one of these equations between each pair of points.
    # So, we have 3 coefficients for each line. Our system of equations, then,
    # is a set of equations that can have any of these coefficients set.
    # So each row in our coefficients matrix contains all coefficients for
    # for all lines. We restrict which ones we're talking about by zeroing out
    # most coefficients, and multipying by zero is zero, so they are neutralized.
    # But, the nice thing is that we can represent any set of coefficients in the
    # same equation by just moving them all to one side and setting their
    # coefficients at a consistent offset.
    #
    # The constants matrix is the RHS of the equation above. In that equation,
    # we have a variable, `y`, but we plug a value into it, thus it can't vary.
    num_lines    = points.length.pred
    coefficients = []
    constants    = []

    # Translates into matrix positions
    emit_row = lambda do |const, start_index, _coefficients, &block|
      row = [0] * num_lines * 3
      _coefficients.each_with_index do |coefficient, coef_index|
        row[start_index*3+coef_index] = coefficient
      end
      block&.call(row)
      quadratic_coefficients << row.map(&:to_f)
      constants              << [const.to_f]
    end

    # Each line segment needs to end on the boundary points
    points.each_cons(2).with_index do |((x1, y1), (x2, y2)), i|
      emit_row[y1, i, [x1**2, x1, 1]]
      emit_row[y2, i, [x2**2, x2, 1]]
    end

    # Lines intersect at the interior points
    _first, *line_intersections, _last = points

    # At these locations, their derivatives must match
    # this is why two adjacent line segments are smooth ("continuous" in math speak)
    line_intersections.each_with_index do |(x, y), i|
      emit_row[0, i, [2*x, 1, 0, -2*x, -1]]
    end

    # We need 1 more equation in order to solve it. A simple enough one is to
    # choose that the two lines at the ends be zero at their second derivative.
    # However, I've somehow come up with setting index 0 and -1 to 1...
    # The 0 makes sense, that would set `a` of line 1 to zero, (which is the same
    # as saying the second derivative is zero), the -1 seems like that would
    # set `c` of the last line, rather than `a`... and then, setting the value to
    # 1 makes no sense, it's the second derivative, so:
    #
    #     ax^2 + bx + c = y
    #     2ax + b = y
    #     2a = y
    #
    # So, second derivative is `2a = y`, we set `y` to zero with the first arg to
    # emit_row, but then shouldn't I set the coefficients to `2` instead of `1`?
    #
    # Unfortunatley, when I do that, it breaks (ExceptionForMatrix::ErrNotRegular)
    # and I can't remember how I came up with the thing below, and it seems wrong...
    # except that the result of this actually looks pretty decent :/
    emit_row[0, 0, []] { |row| row[0] = row[-1] = 1 }

    solutions = Matrix.rows(quadratic_coefficients)
                      .inverse
                      .*(Matrix.rows(constants))
                      .to_a
                      .flatten

    # a line is a set of a,b,c to the quadratic equation ax^2 + bx + c = y
    lines  = solutions.each_slice(3).to_a
    bounds = points.transpose.first.each_cons(2)
    bounds.zip lines
  end

  def coefficients_for(x)
    # it's okay that they overlap on the knot b/c their values match up
    # if we don't do that, then we can't ask for the value at the last knot
    _, line = @bounds_to_coefficients.find { |(x1, x2), _| x1 <= x && x <= x2 }
    line || @lines.last
  end

  def project((x, y))
    [project_x(x), project_y(y)]
  end

  def project_x(x)
    (w-2*margin)*(x-xmin)/∆x+margin
  end

  def project_y(y)
    (h-2*margin)*(y-ymin)/∆y+margin
  end
end


require 'matrix'
include Math
π = PI
original_fn = -> x { sin x }
points = [-π, -π*2/3, -π/3, 0, π/3, π*2/3, π].map { |x| [x, original_fn[x]] }
# => [[-3.141592653589793, -1.2246467991473532e-16],
#     [-2.0943951023931953, -0.8660254037844388],
#     [-1.0471975511965976, -0.8660254037844386],
#     [0, 0.0],
#     [1.0471975511965976, 0.8660254037844386],
#     [2.0943951023931953, 0.8660254037844388],
#     [3.141592653589793, 1.2246467991473532e-16]]

Spline.new(points, original_fn).run