math_test.ts

import * as matrix from './matrix';
import * as vec from './vec';

function expectMatrixSimilar(actual: matrix.Matrix, expected: matrix.Matrix) {
	expect(actual.length).toBe(expected.length);
	actual.forEach((row, i) => {
		row.forEach((v, k) => {
			expect(v).toBeCloseTo(expected[i][k]);
		});
	});
}

describe('matrix', () => {
	test('width', () => {
		expect(
			matrix.width([
				[1, 2],
				[2, 1],
			])
		).toEqual(2);
	});
	test('a checkerboard matrix', () => {
		expect(
			matrix.checkerboard([
				[1, 2],
				[3, 4],
			])
		).toEqual([
			[1, -2],
			[-3, 4],
		]);
	});
	test('.minors', () => {
		expect(
			matrix.minors([
				[3, 0, 2],
				[2, 0, -2],
				[0, 1, 1],
			] as const)
		).toEqual([
			[2, 2, 2],
			[-2, 3, 3],
			[-0, -10, 0], // dumb
		]);
	});
	describe('.identity', () => {
		test.each([
			[[1, 1], matrix.as<1, 1>([[1]] as const)],
			[
				[3, 3],
				matrix.as<3, 3>([
					[1, 0, 0],
					[0, 1, 0],
					[0, 0, 1],
				] as const),
			],

			[
				[2, 3],
				matrix.as<2, 3>([
					[1, 0],
					[0, 1],
					[0, 0],
				] as const),
			],
		] as const)('%#: (%p) => %p', ([i, j], b) => {
			expect(matrix.identity(i, j)).toEqual(b);
		});
	});

	describe('.transpose', () => {
		test.each([
			[
				matrix.as<3, 2>([
					[1, 2, 3],
					[4, 5, 6],
				] as const),

				matrix.as<2, 3>([
					[1, 4],
					[2, 5],
					[3, 6],
				] as const),
			],

			[
				matrix.as<3, 1>([[1, 2, 3]] as const),

				matrix.as<1, 3>([[1], [2], [3]] as const),
			],
		] as const)('%#: (%p) => %p', (a, b) => {
			expect(matrix.transpose(a)).toEqual(b);
		});
	});

	describe('.determinant', () => {
		test.each([
			[matrix.as<1, 1>([[1]] as const), 1],

			[
				matrix.as([
					[2, 0],
					[0, 5],
				] as const),
				10,
			],

			[
				matrix.as<2, 2>([
					[3, 8],
					[4, 6],
				] as const),

				-14,
			],

			[
				matrix.as<3, 3>([
					[6, 1, 1],
					[4, -2, 5],
					[2, 8, 7],
				] as const),

				-306,
			],

			[
				matrix.as([
					[6, -7],
					[-2, 4],
				] as const),
				10,
			],
		] as const)('%#: (%p) => %p', (a, b) => {
			expect(matrix.determinant(a)).toEqual(b);
		});
	});

	describe('.inverse', () => {
		test.each([
			//[matrix.zero, matrix.zero],

			[
				matrix.as([
					[2, 0],
					[0, 5],
				] as const),

				matrix.as([
					[1 / 2, 0],
					[0, 1 / 5],
				] as const),
			],

			/*[
				matrix.as<2, 2>([
					[3, 3.5],
					[3.2, 3.6],
				] as const),

				matrix.as<2, 2>([
					[-9, 8.75],
					[8, -7.5],
				] as const),

				// 1.12, -3.5
				// -3.2, 1.25
			],*/
			[
				matrix.as([
					[3, 0, 2],
					[2, 0, -2],
					[0, 1, 1],
				] as const),

				matrix.as([
					[0.2, 0.2, 0],
					[-0.2, 0.3, 1],
					[0.2, -0.3, 0],
				] as const),
			],

			[
				matrix.as<2, 2>([
					[4, 7],
					[2, 6],
				] as const),

				matrix.as<2, 2>([
					[0.6, -0.7],
					[-0.2, 0.4],
				] as const),
			],
		] as const)('%#: (%p) => %p', (a, b) => {
			expectMatrixSimilar(matrix.inverse(a), b);
		});

		test.each([
			/*[
				matrix.as<2, 2>([
					[3, 3.5],
					[3.2, 3.6],
				] as const),

				matrix.as<2, 2>([
					[-9, 8.75],
					[8, -7.5],
				] as const),

				// 1.12, -3.5
				// -3.2, 1.25
			],*/

			[
				matrix.as<2, 2>([
					[4, 7],
					[2, 6],
				] as const),
			],
			[
				matrix.as([
					[20, 100, 2],
					[0.2, -10, 1],
					[4, 3, 1],
				] as const),
			],
		] as const)('%#: %p * a => identity', a => {
			const [ij = 0] = matrix.size(a);
			expectMatrixSimilar(
				matrix.mul(a, matrix.inverse(a)),
				matrix.identity(ij, ij)
			);
		});
	});

	describe('.transpose', () => {
		test.each([
			[
				matrix.as<3, 2>([
					[1, 2, 3],
					[4, 5, 6],
				] as const),

				matrix.as<2, 3>([
					[1, 4],
					[2, 5],
					[3, 6],
				] as const),
			],

			[
				matrix.as<3, 1>([[1, 2, 3]] as const),

				matrix.as<1, 3>([[1], [2], [3]] as const),
			],
		] as const)('%#: (%p) => %p', (a, b) => {
			expect(matrix.transpose(a)).toEqual(b);
		});
	});

	describe('.col', () => {
		test.each([
			[
				matrix.as<2, 3>([
					[1, 2],
					[1, 4],
					[6, 5],
				] as const),
				0,
				[1, 1, 6],
			],
		])('%#: (%p, %p) => %p', (a, b, o) => {
			expect([...matrix.col(a, b)]).toEqual(o);
		});
	});

	describe('.row', () => {
		test.each([
			[matrix.as([[1]]), 0, [1]],
			[matrix.as([[1]]), 1, []],
			[matrix.as([[1]]), -1, []],
			[matrix.as([[1]]), Infinity, []],
		])('%#: (%p, %p) => %p', (a, b, o) => {
			expect([...matrix.row(a, b)]).toEqual(o);
		});

		test.each([
			[
				matrix.as<2, 3>([
					[1, 2],
					[1, 4],
					[6, 5],
				] as const),
				0,
				[1, 2],
			],
		])('(%p, %p) => %p', (a, b, o) => {
			expect([...matrix.row(a, b)]).toEqual(o);
		});
	});

	describe('.mul', () => {
		test.each([
			[
				matrix.as<2, 1>([[6, 9]] as const),

				matrix.as<2, 2>([
					[0, 1],
					[1, 0],
				] as const),

				matrix.as<2, 1>([[9, 6]] as const),
			],
			[
				matrix.as<2, 2>([
					[1, 0],
					[0, 1],
				] as const),

				matrix.as<2, 2>([
					[2, 4],
					[10, 9],
				] as const),

				matrix.as<2, 2>([
					[2, 4],
					[10, 9],
				] as const),
			],
		])('%#: (%p, %p) => %p', (a, b, o) => {
			expect(matrix.mul(a, b)).toEqual(o);
		});
	});

	describe('.add', () => {
		test.each([
			[
				matrix.as([
					[1, 2],
					[3, 4],
				] as const),

				matrix.as([
					[4, 3],
					[2, 1],
				] as const),

				matrix.as([
					[5, 5],
					[5, 5],
				]),
			],
		])('%#: (%p, %p) => %p', (a, b, o) => {
			expect(matrix.add(a, b)).toEqual(o);
		});
	});
});

describe('vec', () => {
	describe('.dot', () => {
		test.each([[[1, 2, 3], [3, 2, 1], 10]])(
			'.dot(%p, %p) => %p',
			(a, b, o) => {
				expect(vec.dot(a, b)).toEqual(o);
			}
		);
	});

	describe('.reverse', () => {
		test('[1,2,3] => [3,2,1]', () => {
			expect([...vec.reverse([1, 2, 3])]).toEqual([3, 2, 1]);
		});
	});

	describe('.mul', () => {
		test.each([[2, [3, 2, 1], [6, 4, 2]]])(
			'.mul(%p, %p) => %p',
			(a, b, o) => {
				expect(vec.mul(a, b)).toEqual(o);
			}
		);
	});
});

matrix.ts

import { Vector } from './vec';
import * as vec from './vec';

// J is effectively the number of ROWS and I is the number of COLUMNS
export type Matrix<
	I extends number = number,
	J extends number = number,
	T = number
> = Vector<J, Vector<I, T>>;

export type Square<IJ extends number, T = number> = Matrix<IJ, IJ, T>;

export const as: <
	I extends number = number,
	J extends number = number,
	T = number
>(
	v: readonly (readonly T[] & { length: I })[] & { length: J }
	/* eslint-disable-next-line @typescript-eslint/no-explicit-any */
) => Matrix<I, J, T> = v => v as any;

/* eslint-disable-next-line @typescript-eslint/no-explicit-any */
export const zero = as<0, 0>([] as any);

/**
 * Return a new matrix of given dimensions
 */
export function New<I extends number, J extends number>(
	i: I,
	j: J
): Matrix<I, J, undefined> {
	return vec.map(vec.New(j), () => vec.New(i));
}

export const add: <I extends number, J extends number>(
	m1: Matrix<I, J>,
	m2: Matrix<I, J>
) => Matrix<I, J> = <I extends number, J extends number>(
	m1: Matrix<I, J>,
	m2: Matrix<I, J>
) => vec.map(m1, (row, i) => vec.add(row, m2[i]));

/**
 * Returns a row of a given matrix as an Iterable.
 * Where the row does not exist, the Iterable is of length 0.
 */
export const row: <I extends number, J extends number, T>(
	v: Matrix<I, J, T>,
	r: number
) => Iterable<T> = function* (v, i) {
	const a = v[i];
	if (!a) return;
	for (let i = 0; i < a.length; i++) yield a[i];
};

export const rows: <I extends number, J extends number, T>(
	v: Matrix<I, J, T>,
	r: number
) => Iterable<Vector<I, T>> = v => v;

export const col: <I extends number, J extends number, T>(
	v: Matrix<I, J, T>,
	i: number
) => Iterable<T> = function* (v, i) {
	const [, jsize] = size(v);
	for (let j = 0; j < jsize; j++) yield v[j][i];
};

export const mul: <
	I1 extends number,
	J1 extends number,
	I2 extends number,
	J2 extends number
>(
	m1: Matrix<I1, J1>,
	m2: Matrix<I2, J2>
) => Multiply<Matrix<I1, J1>, Matrix<I2, J2>> = <
	I1 extends number,
	J1 extends number,
	I2 extends number,
	J2 extends number
>(
	m1: Matrix<I1, J1>,
	m2: Matrix<I2, J2>
) => {
	const [, /*i1*/ j1] = size(m1);
	const [i2 /*, j2*/] = size(m2);

	return vec.map(vec.New<J1>(j1), (_, i) =>
		vec.map(vec.New<I2>(i2), (_, j) => vec.dot(row(m1, i), col(m2, j)))
	);
};

export const map: <I extends number, J extends number, T, O>(
	m: Matrix<I, J, T>,
	f: (
		v: T,
		pos: readonly [i: number, j: number],
		matrix: Matrix<I, J, T>
	) => O
) => Matrix<I, J, O> = (m, f) =>
	vec.map(m, (row, j) => vec.map(row, (v, i) => f(v, [i, j], m)));

/**
 * Unsafely drops values in the matrix that do not make f return true.
 *
 * Unsafe in the sense that it may not return a valid matrix (since rows / cols might be incorrect)
 */
export const filter: <T>(
	m: Matrix<number, number, T>,
	f: (
		v: T,
		pos: readonly [i: number, j: number],
		matrix: Matrix<number, number, T>
	) => boolean
) => readonly (readonly T[])[] = (m, f) =>
	m
		.map((row, rowIndex) =>
			row.filter((v, colIndex) => f(v, [colIndex, rowIndex], m))
		)
		.filter(row => row.length !== 0);

export function is<T>(
	v: readonly (readonly T[])[]
): v is Matrix<number, number, T> {
	// each row must be of the same size

	return v.every((row, _, a) => row.length == a[0].length);
}

export function must<T>(
	v: readonly (readonly T[])[]
): asserts v is Matrix<number, number, T> {
	if (!is(v)) throw new Error(`${JSON.stringify(v)} is not a valid matrix`);

	return;
}

function isSquare<T>(m: readonly (readonly T[])[]): m is Square<number, T> {
	return m.every(row => row.length == m.length);
}

function mustIsSquare<T>(
	v: readonly (readonly T[])[]
): asserts v is Square<number, T> {
	if (!isSquare(v))
		throw new Error(`${JSON.stringify(v)} is not a valid square matrix`);
}

/**
 * `Multiply` gives the type of the matrix made by multiplying 2 given matricies.
 */
export type Multiply<
	A extends Matrix<number, number, unknown>,
	B extends Matrix<number, number, unknown>,
	O = number
> = [A, B] extends [
	/* eslint-disable-next-line @typescript-eslint/no-explicit-any */
	Matrix<any, infer J1, unknown>,
	/* eslint-disable-next-line @typescript-eslint/no-explicit-any */
	Matrix<infer I2, any, unknown>
]
	? Matrix<I2, J1, O>
	: never;

/**
 * `TransformTo` gives the type of the matrix that can be used to transform an input matrix
 * into an output matrix
 */
export type TransformTo<
	In extends Matrix<number, number, unknown>,
	Out extends Matrix<number, number, unknown>
	/*O = number*/
> = [In, Out] extends [
	/* eslint-disable-next-line @typescript-eslint/no-explicit-any */
	Matrix<any, infer J1, unknown>,
	Matrix<infer I2, infer J2, unknown>
]
	? J2 extends J1
		? /* eslint-disable-next-line @typescript-eslint/no-explicit-any */
		  Matrix<I2, any>
		: never // it isn't possible for these to be different
	: never;

export const size: <I extends number, J extends number>(
	/* eslint-disable-next-line @typescript-eslint/no-explicit-any */
	m: Matrix<I, J, any>
) => J extends 0 ? [undefined, J] : [I, J] = m =>
	/* eslint-disable-next-line @typescript-eslint/no-explicit-any */
	[m[0]?.length, m.length] as any;

export const transpose: <I extends number, J extends number>(
	m: Matrix<I, J>
) => Matrix<J, I> = <I extends number, J extends number>(m: Matrix<I, J>) => {
	const [i, j] = size(m);
	const rows = vec.New<I>(i);

	return vec.map(rows, (_, rj) =>
		vec.map(vec.New<J>(j), (__, vi) => m[vi][rj])
	);
};

/**
 * Returns an identity matrix of given dimensions
 */
export const identity: <I extends number, J extends number>(
	I: I,
	J: J
) => Matrix<I, J, 1 | 0> = <I extends number, J extends number>(i: I, j: J) => {
	return map(New<I, J>(i, j), (_, [i, j]) => +!(i - j) as 1 | 0);
};

// https://www.mathsisfun.com/algebra/matrix-determinant.html
export const determinant: <IJ extends number>(m: Square<IJ>) => number = m => {
	const [ij] = size(m);

	if ((ij ?? 0) == 0) return 0;
	if (ij == 1) return m[0][0];

	if (ij == 2) {
		const [[a, b], [c, d]] = m;
		return a * d - b * c;
	}

	const top = m[0];
	const rest = m.slice(1); // of length J - 1;
	return top.reduce((acc, cur, ind) => {
		const mul = ind % 2 == 0 ? 1 : -1;
		// remove every value in the same row
		const nm = rest.map(row => row.filter((_, ri) => ri !== ind));
		console.assert(nm[0].length == ij - 1);
		// this should give us a new matrix we can get the determinant for
		return acc + cur * determinant(nm) * mul;
	}, 0);
};

//https://www.mathsisfun.com/algebra/matrix-inverse-minors-cofactors-adjugate.html
export const minors: <IJ extends number>(s: Square<IJ>) => Square<IJ> = s =>
	map(s, (v, [column, row], m) => {
		// form a new matrix omitting the row and column of the selected value
		const smaller = filter(
			m,
			(_, [sColumn, sRow]) => column !== sColumn && row !== sRow
		);

		mustIsSquare(smaller);

		return determinant(smaller);
	});

export const width: (m: Matrix) => number = m => m?.[0]?.length ?? 0;

export const checkerboard: <I extends number, J extends number>(
	m: Matrix<I, J>
) => Matrix<I, J> = m =>
	map(m, (n, [col, row]) => ((col + row) % 2 == 0 ? n : -n));

/**
 * Returns the inverse of a matrix of given dimensions
 */
export const inverse: <IJ extends number>(m: Square<IJ>) => Square<IJ> = <
	IJ extends number
>(
	m: Square<IJ>
) => {
	const d = 1 / determinant(m);

	return map(transpose(checkerboard(minors(m))), n => d * n);
};

vec.ts

interface ReadonlyArrayOfLength<T, I extends number> extends ReadonlyArray<T> {
	length: I;
}

export interface Vector<I extends number = number, T = number>
	extends ReadonlyArrayOfLength<T, I> {
	length: I;
}

type mapFn<I extends number, T> = <U>(
	callbackFn: (value: T, index: number, array: Vector<I, T>) => U,
	thisArg?: unknown
) => Vector<I, U>;

/**
 * Map a Vector, returning a new Vector.
 */
export const map: <I extends number, T, U>(
	vec: Vector<I, T>,
	callbackFn: (value: T, index: number, array: Vector<I, T>) => U
	/* eslint-disable-next-line @typescript-eslint/no-explicit-any */
) => Vector<I, U> = (vec, c) => (vec.map as mapFn<any, any>)(c);

/**
 * Map an Iterable, returning a new Iterable.
 */
export const imap: <T, U>(
	v: Iterable<T>,
	f: (v: T, i: number) => U
) => Iterable<U> = function* (v, f) {
	let i = 0;
	for (const l of v) {
		yield f(l, i);
		i++;
	}
};

/**
 * Returns an iterable on a vector that iterates in reverse
 */
export const reverse: <I extends number, T>(v: Vector<I, T>) => Iterable<T> =
	function* (v) {
		for (let i = 0; i < v.length; i++) {
			yield v[v.length - i - 1];
		}
	};

export const as: <T, L extends number>(
	v: readonly T[] & { length: L }
	/* eslint-disable-next-line @typescript-eslint/no-explicit-any */
) => Vector<L, T> = v => v as any;

export const New: <N extends number>(n: number) => Vector<N, undefined> = n =>
	/* eslint-disable-next-line @typescript-eslint/no-explicit-any */
	[...Array(n)] as any as Vector<any, undefined>;

export const add: <I extends number>(
	v1: Vector<I>,
	v2: Vector<I>
) => Vector<I> = (v1, v2) => map(v1, (v, i) => v + v2[i]);

export const mul: <I extends number>(v1: number, v2: Vector<I>) => Vector<I> = (
	v1,
	v2
) => map(v2, (v /*i*/) => v * v1);

export const dot: (v1: Iterable<number>, v2: Iterable<number>) => number = (
	v1,
	v2
) =>
	sum(
		imap(
			zip(v1, v2, 0 as const),
			([a = 0 as const, b = 0 as const]) => a * b
		)
	);

export const sum: (v1: Iterable<number>) => number = v =>
	[...v].reduce((a, c) => a + c, 0);

const _zip: <T1, T2, T3>(
	v1: Iterable<T1>,
	v2: Iterable<T2>,
	fb: T3
) => Iterable<[T1 | T3, T2 | T3]> = function* (v1, v2, fb) {
	const [a, b] = [v1[Symbol.iterator](), v2[Symbol.iterator]()];
	for (let i = 0; ; i++) {
		const [ar, br] = [a.next(), b.next()];
		const left = ar.done ? fb : ar.value;
		const right = br.done ? fb : br.value;

		if (ar.done && br.done) break;

		yield [left, right];
	}
};
export const zip: {
	<T1, T2, L extends number>(
		v1: ReadonlyArrayOfLength<T1, L>,
		v2: ReadonlyArrayOfLength<T2, L>
	): Iterable<[T1, T2]>;
	<T1, T2, T3>(v1: Iterable<T1>, v2: Iterable<T2>, fb: T3): Iterable<
		[T1 | T3, T2 | T3]
	>;
	/* eslint-disable-next-line @typescript-eslint/no-explicit-any */
} = _zip as any;