Skip to content

Instantly share code, notes, and snippets.

@jdmichaud

jdmichaud/index.html

Last active June 10, 2018 20:58
Show Gist options
  • Save jdmichaud/f6e30e46e6ea47a7380a6f2d5bb04a53 to your computer and use it in GitHub Desktop.
Save jdmichaud/f6e30e46e6ea47a7380a6f2d5bb04a53 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
index.html
<!DOCTYPE html>
<html>
<head>
<title>Least Squares</title>
<script src="https://cdnjs.cloudflare.com/ajax/libs/mathjs/4.4.2/math.js"></script>
<script src="main.js"></script>
<style>
.graph {
border: solid 1px black;
display: block;
}
</style>
</head>
<body>
<select class="orderSelector">
<option value=1>linear</option>
<option value=2>quadratic</option>
<option value=3>cubic</option>
<option value=4>quartic</option>
<option value=5>quintic</option>
</select>
<canvas class="graph" width=640 height=480></canvas>
</body>
</html>
Raw
least_square.py
from scipy import linalg
import numpy as np
# If A * x = b has no solution, x' can be approximated by multiplying by A.T:
# A.T * A * x' = A.T * b
# If the equation A * x = b represent the parameters of a set of equations that
# are supposed to match a list of data points.
# ex.:
# A list of point (0, 6) (1, 0) and (2, 0) should be matched by a affine
# function 'a * x + b = y'
# The set of equation would be:
# 0 * x + 1 * b = 6
# 1 * x + 1 * b = 0
# 2 * x + 1 * b = 0
# [ 0 1 ] [6]
# A =[ 1 1 ] x = [0]
# [ 2 1 ] [0]
# For complete proof, see:
# Introduction to Linear Algebra 4th ed., Ch. 4.3 (by Gilbert Strang)
# This function fits a function of order to the set of points
def least_square(points, order):
col1 = [p[0] for p in points]
A = np.mat([[1] * len(points)] + [[c**i for c in col1] for i in np.arange(1, order)]).T
b = np.mat([p[1] for p in points]).T
return linalg.inv(A.T * A) * A.T * b
Raw
linealg.md 

Ax = b

    |P..|  r = n = m
A = |.P.|  A is invertible
    |..P|  One solution for every b

    |P..|
    |.P.|  r = n < m
A = |...|  A is full column rank
    |..P|  One or no solution
    |...|

A = |.P....|  r = m < n
    |...P..|  A is full row rank
              One or infinity solutions

    |....|
A = |P...|  r < m & r < n
    |....|  Infinity of solutions
    |..P.|
Raw
main.js
let dots = [];
let line = [];
let order = 1;
// Default convertion function. Invert y and adds pan.
function convertToCanvas(canvas, pan, coord) {
return [coord[0] + pan[0], canvas.height - coord[1] - pan[1]];
}
function convertFromCanvas(canvas, pan, coord) {
const rect = canvas.getBoundingClientRect();
const x = coord[0] - rect.left;
const y = coord[1] - rect.top;
return [x - pan[0] - 3, canvas.height - y - pan[1] + 4];
}
// Draw the frame, dots and line on the canvas
// using the provided conversion function.
function draw(canvas, conv, dots, line) {
const ctx = canvas.getContext("2d");
function drawDot(d) {
ctx.fillRect(d[0], d[1], 4, 4);
}
function drawLine(from, to) {
ctx.beginPath();
ctx.moveTo(from[0], from[1]);
ctx.lineTo(to[0], to[1]);
ctx.stroke();
}
function drawFrame() {
drawLine(conv([0, 0]), conv([canvas.width - 20, 0]));
drawLine(conv([0, 0]), conv([0, canvas.height - 20]));
}
ctx.clearRect(0, 0, canvas.width, canvas.height);
drawFrame();
dots.forEach(d => drawDot(conv(d)));
// if (line.length == 2) drawLine(conv(line[0]), conv(line[1]));
line.forEach(p => ctx.fillRect(conv(p)[0], conv(p)[1], 1, 1));
}
// Compute the least square
// $$xh = (A^tA)^{-1}A^tb$$
function leastSquare(dots, order) {
const a = [];
for (let i = order; i >= 0; --i) {
a.push(dots.map(d => math.pow(d[0], i)));
}
const A = math.transpose(math.matrix(a));
const b = math.transpose(math.matrix(dots.map(d => d[1])));
const xh = math.multiply(
math.multiply(
math.inv(
math.multiply(
math.transpose(A),
A,
),
),
math.transpose(A),
),
b,
);
return xh.toArray();
}
// Apply the function as a polynomial of order params.length
function apply(canvas, pan, params) {
const left = -pan[0];
const right = canvas.width - pan[0];
const line = [];
for (let i = left; i < right; ++i) {
let res = 0;
for (let j = 0; j < params.length; ++j) {
res += params[j] * math.pow(i, params.length - j - 1);
}
line.push([i, res]);
}
console.log(line);
return line;
}
function mousedown(canvas, pan, toCanvas, fromCanvas, event) {
if (event.button == 0) {
// Right click, add dot
dots.push(fromCanvas([event.clientX, event.clientY]));
} else if (event.button == 2) {
// Left click, remove dot
const coord = fromCanvas([event.clientX, event.clientY]);
dots = dots.filter(d =>
// distance between a and b:
// (a - b)^t * (a - b)
math.multiply(
math.transpose(math.subtract(coord, d)),
math.subtract(coord, d),
) > 200
);
}
if (dots.length > 1) {
line = apply(canvas, pan, leastSquare(dots, order));
} else {
line = [];
}
draw(canvas, toCanvas, dots, line);
}
function orderChange(canvas, pan, toCanvas, event) {
order = parseInt(event.target.value, 10);
line = apply(canvas, pan, leastSquare(dots, order));
draw(canvas, toCanvas, dots, line);
}
// Retrieve canvas, define the pan, bind events and draw
function main() {
const pan = [10, 10];
const canvas = document.getElementsByTagName('canvas')[0];
const toCanvas = convertToCanvas.bind(this, canvas, pan);
const fromCanvas = convertFromCanvas.bind(this, canvas, pan);
const selector = document.getElementsByClassName('orderSelector')[0];
canvas.addEventListener('mousedown',
mousedown.bind(this, canvas, pan, toCanvas, fromCanvas));
canvas.oncontextmenu = (e) => e.preventDefault();
selector.addEventListener('change',
orderChange.bind(this, canvas, pan, toCanvas));
draw(canvas, toCanvas, dots, line);
}
window.onload = main;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment