Added a variety of vector functions. Modified linear solver to use partial pivoting. Added a tridiagonal solver.

Author: adam-rumpf

GameMakerScripts.yyp

@@ -12,23 +12,22 @@ Common array functions (such as searching and counting).
 
 * [`_array_count`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_array_count/_array_count.gml): Counts the number of occurrences of a specified value in an array.
 * [`_array_index`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_array_index/_array_index.gml): Finds the first index at which a specified value occurs in an array (or -1 if not found). An optional argument allows the search to begin with a specified index.
-* [`_array_max`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_array_max/_array_max.gml): Returns the maximum value in an array.
-* [`_array_min`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_array_min/_array_min.gml): Returns the minimum value in an array.
-
-## Computational Mathematics Functions
-
-Algorithms from computational mathematics (such as regression and interpolation).
-
-* `...`
+* [`_array_max`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_array_max/_array_max.gml): Returns the maximum value (or the index of the maximum value) in an array.
+* [`_array_min`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_array_min/_array_min.gml): Returns the minimum value (or the index of the minimum value) in an array.
+* [`_invert_permutation`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_invert_permutation/_invert_permutation.gml): Generates a permutation array to invert a given permutation array.
+* [`_linspace`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_linspace/_linspace.gml): Generates an array with a specified number of equally-spaced values that cover a specified range.
+* [`_permute`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_permute/_permute.gml): Permutes the elements of an array according to a given permutation array.
+* [`_range`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_range/_range.gml): Generates an array of equally-spaced values over a specified range with a specified step size.
 
 ## Linear Algebra Functions
 
 Common linear algebra algorithms for dealing with matrices and vectors. In all functions, _vectors_ are considered to be 1D arrays while _matrices_ are 2D arrays.
 
-* [`_linear_solve`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_linear_solve/_linear_solve.gml): Solves a linear system given a coefficient matrix and righthand vector.
+* [`_linear_solve`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_linear_solve/_linear_solve.gml): Solves a linear system given a coefficient matrix and righthand vector (using Gaussian elimination with partial pivoting).
 * [`_matrix_multiply`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_matrix_multiply/_matrix_multiply.gml): Performs matrix-matrix, matrix-vector, and vector-vector multiplication.
 * [`_matrix_transpose`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_matrix_transpose/_matrix_transpose.gml): Transposes a 2D matrix.
-* [`_triangular_solve`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_triangular_solve/_triangular_solve.gml): Solves an upper or lower triangular system.
+* [`_tridiagonal_solve`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_tridiagonal_solve/_tridiagonal_solve.gml): Solves a tridiagonal system (using the Thomas algorithm).
+* [`_triangular_solve`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_triangular_solve/_triangular_solve.gml): Solves an upper or lower triangular system (using forward or back substitution).
 * [`_unit_vector`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_unit_vector/_unit_vector.gml): Defines a unit direction vector between two coordinates.
 * [`_vector_angle`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_vector_angle/_vector_angle.gml): Calculates the acute angle between two vectors.
 * [`_vector_distance`](https://github.com/adam-rumpf/game-maker-scripts/blob/master/scripts/_vector_distance/_vector_distance.gml): Calculates the Lp-distance between two vectors.

README.md

@@ -1,18 +1,29 @@
-/// @func _array_max(arr)
-/// @desc Finds the maximum value in an array.
+/// @func _array_max(arr[, val=true])
+/// @desc Finds the maximum value (or its index) in an array.
 /// @param {real[]} arr Array to search.
-/// @return {real} Maximum value in arr.
+/// @param {bool} [val=true] If true the maximum value is returned, if false the (first) index of the maximum value is returned.
+/// @return {real} Maximum value (or its index) in arr.
 
 function _array_max(arr)
 {
+	// Check for optional index argument
+	var val = (argument_count > 1 ? argument[1] : true);
+	
 	// Go through array to find largest value
-	var val = -infinity;
+	var mv = -infinity; // maximum value
+	var mi = -1; // index of maximum value
 	for (var i = 0; i < array_length(arr); i++)
 	{
-		if (arr[i] > val)
-			val = arr[i];
+		if (arr[i] > mv)
+		{
+			mv = arr[i];
+			mi = i;
+		}
 	}
 
-	// Return largest value
-	return val;
+	// Return largest value or its index
+	if (val == true)
+		return mv;
+	else
+		return mi;
 }

scripts/_array_max/_array_max.gml

@@ -1,18 +1,29 @@
-/// @func _array_min(arr)
-/// @desc Finds the minimum value in an array.
+/// @func _array_min(arr[, val=true])
+/// @desc Finds the minimum value (or its index) in an array.
 /// @param {real[]} arr Array to search.
-/// @return {real} Minimum value in arr.
+/// @param {bool} [val=true] If true the minimum value is returned, if false the (first) index of the minimum value is returned.
+/// @return {real} Minimum value (or its index) in arr.
 
 function _array_min(arr)
 {
+	// Check for optional index argument
+	var val = (argument_count > 1 ? argument[1] : true);
+	
 	// Go through array to find smallest value
-	var val = infinity;
+	var mv = infinity; // minimum value
+	var mi = -1; // index of minimum value
 	for (var i = 0; i < array_length(arr); i++)
 	{
-		if (arr[i] < val)
-			val = arr[i];
+		if (arr[i] < mv)
+		{
+			mv = arr[i];
+			mi = i;
+		}
 	}
 
-	// Return smallest value
-	return val;
+	// Return smallest value or its index
+	if (val == true)
+		return mv;
+	else
+		return mi;
 }

scripts/_array_min/_array_min.gml

@@ -0,0 +1,26 @@
+/// @func _invert_permutation(perm)
+/// @desc Creates the inverse of a given permutation array.
+/// @param {int[]} perm Permutation array, consisting of a sequence of indices 0, 1, 2, ...
+/// @return {int[]} Permutation array which inverts the permutation defined by perm (undefined in case of error).
+
+function _invert_permutation(perm)
+{
+	// Get array length
+	var n = array_length(perm);
+	
+	// Initialize inverse
+	var inv = array_create(n);
+	
+	// Determine values of inverse array
+	for (var i = 0; i < n; i++)
+	{
+		// Verify that index is within bounds
+		if ((perm[i] < 0) || (perm[i] >= n))
+			return undefined;
+		
+		// Reverse mapping
+		inv[perm[i]] = i;
+	}
+	
+	return inv;
+}

scripts/_invert_permutation/_invert_permutation.gml

@@ -1,8 +1,7 @@
 /// @func _linear_solve(a, b)
-/// @desc Solves a linear system of the form Ax = b (via Gaussian elimination) for a given coefficient matrix and constant vector.
-/// @requires _triangular_solve
-/// @param {real[]} a Square coefficient matrix (as a 2D array).
-/// @param {real[]} b Constraint vector (as a 1D array whose length matches a).
+/// @desc Solves a linear system of the form Ax = b (via Gaussian elimination with partial pivoting) for a given coefficient matrix and constant vector.
+/// @param {real[][]} a Square coefficient matrix.
+/// @param {real[]} b Constant vector (as an array whose length matches a).
 /// @return {real[]} Solution vector x of the linear system Ax = b (undefined in case of inconsistent dimensions or singular matrix).
 
 function _linear_solve(a, b)
@@ -12,31 +11,60 @@ function _linear_solve(a, b)
 		return undefined;
 
 	// Initialize system
-	var n, l, u;
+	var n, aa, bb, xx, perm;
 	n = array_length(b); // size of system
-	l = array_create(n, array_create(n, 0)); // lower triangular matrix
+	aa = a; // diagonalized matrix
+	bb = b; // constant vector
+	xx = array_create(n); // solution vector
+	perm = array_create(n); // row permutation array
 	for (var i = 0; i < n; i++)
-		l[i][i] = 1; // initialize as identity matrix
-	u = a; // upper triangular matrix
+		perm[i] = i;
 
-	// Perform Gaussian elimination to obtain LU decomposition
-	for (var i = 0; i < n-1; i++)
+	// Perform Gaussian elimination with partial pivoting
+	for (var i = 0; i < n; i++)
 	{
+		// Find maximum absolute value in current column's subdiagonal
+		var p = i; // pivot row index
 		for (var j = i+1; j < n; j++)
 		{
-			// Verify that we are not dividing by zero
-			if (u[i][i] == 0)
-				return undefined;
+			if (abs(aa[perm[j]][i]) > abs(aa[perm[p]][i]))
+				p = j;
+		}
+		
+		// Verify that the element is nonzero
+		if (aa[perm[p]][i] == 0)
+			return undefined;
+		
+		// Set permutation array to exchange the current row and the pivot row
+		if (i != p)
+		{
+			var temp = perm[i];
+			perm[i] = perm[p];
+			perm[p] = temp;
+		}
+		
+		// Perform elimination using pivot element
+		for (var j = 0; j < n; j++)
+		{
+			// Ignore pivot row
+			if (j == i)
+				continue;
+			
+			// Calculate multiplicative factor
+			var multi = aa[perm[j]][i]/aa[perm[i]][i];
 
-			l[j][i] = u[j][i]/u[i][i];
+			// Eliminate row
 			for (var k = i; k < n; k++)
-				u[j][k] = u[j][k] - l[j][i]*u[i][k];
+				aa[perm[j]][k] -= multi*aa[perm[i]][k];
+			
+			// Transform constant vector
+			bb[perm[j]] -= multi*bb[perm[i]];
 		}
 	}
 
-	// Solve Ly = b
-	var yy = _triangular_solve(l, b, false);
+	// Solve diagonalized system
+	for (var i = 0; i < n; i++)
+		xx[i] = bb[perm[i]]/aa[perm[i]][i];
 
-	// Solve Ux = y
-	return _triangular_solve(u, yy, true);
+	return xx;
 }

scripts/_linear_solve/_linear_solve.gml

@@ -0,0 +1,31 @@
+/// @func _linspace(start, stop, n)
+/// @desc Generates a specified number of equally-spaced values within a specified range.
+/// @param {real} start First value of array.
+/// @param {real} stop Last value of array.
+/// @param {int} n Number of values in array. Must be at least 2 (unless the endpoints are equal).
+/// @return {real[]} An array with first value start, final value stop, and a total of n elements (or empty array in case of error).
+
+function _linspace(start, stop, n)
+{
+	// Verify that number of elements is valid
+	if (n < 2)
+	{
+		// n=1 is allowed only if the endpoints are equal
+		if ((start == stop) && (n == 1))
+			return [start];
+		
+		// Otherwise we return an empty array
+		return [];
+	}
+	
+	// Determine step size
+	var step = (stop - start)/(n - 1);
+	
+	// Generate array
+	var arr = array_create(n);
+	for (var i = 0; i < n-1; i++)
+		arr[i] = start + i*step;
+	arr[n-1] = stop; // ensure that final value is exact
+	
+	return arr;
+}

scripts/_linspace/_linspace.gml

@@ -1,7 +1,7 @@
 /// @func _matrix_multiply(a, b)
 /// @desc Multiplies two matrices and/or vectors of compatible dimensions.
-/// @param {real[]} a 1D or 2D array to use as the left factor of the product.
-/// @param {real[]} b 1D or 2D array to use as the right factor of the product.
+/// @param {real[]*} a 1D or 2D array to use as the left factor of the product.
+/// @param {real[]*} b 1D or 2D array to use as the right factor of the product.
 /// @return {real[]} Product of matrix multiplication (undefiend if dimensions not compatible; scalar for vector-vector products).
 
 function _matrix_multiply(a, b)

scripts/_matrix_multiply/_matrix_multiply.gml

@@ -1,7 +1,7 @@
 /// @func _matrix_transpose(a)
 /// @desc Transposes a 2D array.
-/// @param {*[]} a Matrix to transpose.
-/// @return {*[]} The transpose of the input array a.
+/// @param {*[][]} a Matrix to transpose.
+/// @return {*[][]} The transpose of the input array a.
 
 function _matrix_transpose(a)
 {

scripts/_matrix_transpose/_matrix_transpose.gml

@@ -0,0 +1,28 @@
+/// @func _permute(arr, perm)
+/// @desc Permutes a given array according to a given permutation array.
+/// @param {*[]} arr Array to be permuted.
+/// @param {int[]} perm Permutation array (length must match arr and consist of consecutive integers 0, 1, 2, ...).
+/// @return {*[]} Array consisting of values in arr permutated according to perm (or original array if lengths mismatch).
+
+function _permute(arr, perm)
+{
+	// Verify that array lengths match
+	if (array_length(arr) != array_length(perm))
+		return arr;
+	
+	// Get array length
+	var n = array_length(arr);
+	
+	// Apply permutation
+	var arrp = array_create(n); // permuted array
+	for (var i = 0; i < n; i++)
+	{
+		// Verify that index is within bounds
+		if ((perm[i] < 0) || (perm[i] >= n))
+			return arr;
+		
+		arrp[i] = arr[perm[i]];
+	}
+	
+	return arrp;
+}