openscad/BOSL2
1
0
Fork 0
mirror of https://github.com/revarbat/BOSL2.git synced 2025-08-26 23:44:32 +02:00
Code Issues Projects Releases Wiki Activity

Move linear algebra to linalg.scad

Browse Source 
columns->column because the multiindex case is handled by submatrix
and also it never occurs in the code.
This commit is contained in:
Adrian Mariano
2021-10-26 16:45:14 -04:00
parent b66b964188
commit a75a368dd5
15 changed files with 812 additions and 992 deletions
Download Patch File Download Diff File Expand all files Collapse all files

344
math.scad
View File

@@ -5,7 +5,6 @@
// include <BOSL2/std.scad>
//////////////////////////////////////////////////////////////////////
// Section: Math Constants
// Constant: PHI
@@ -669,7 +668,7 @@ function lcm(a,b=[]) =
function sum(v, dflt=0) =
v==[]? dflt :
assert(is_consistent(v), "Input to sum is non-numeric or inconsistent")
is_vector(v) || is_matrix(v) ? [for(i=v) 1]*v :
is_finite(v[0]) || is_vector(v[0]) ? [for(i=v) 1]*v :
_sum(v,v[0]*0);
function _sum(v,_total,_i=0) = _i>=len(v) ? _total : _sum(v,_total+v[_i], _i+1);
@@ -799,18 +798,6 @@ function _cumprod_vec(v,_i=0,_acc=[]) =
);
// Function: outer_product()
// Usage:
// x = outer_product(u,v);
// Description:
// Compute the outer product of two vectors, a matrix.
// Usage:
// M = outer_product(u,v);
function outer_product(u,v) =
assert(is_vector(u) && is_vector(v), "The inputs must be vectors.")
[for(ui=u) ui*v];
// Function: mean()
// Usage:
// x = mean(v);
@@ -877,335 +864,6 @@ function convolve(p,q) =
// Section: Matrix math
// Function: ident()
// Usage:
// mat = ident(n);
// Topics: Affine, Matrices
// Description:
// Create an `n` by `n` square identity matrix.
// Arguments:
// n = The size of the identity matrix square, `n` by `n`.
// Example:
// mat = ident(3);
// // Returns:
// // [
// // [1, 0, 0],
// // [0, 1, 0],
// // [0, 0, 1]
// // ]
// Example:
// mat = ident(4);
// // Returns:
// // [
// // [1, 0, 0, 0],
// // [0, 1, 0, 0],
// // [0, 0, 1, 0],
// // [0, 0, 0, 1]
// // ]
function ident(n) = [
for (i = [0:1:n-1]) [
for (j = [0:1:n-1]) (i==j)? 1 : 0
]
];
// Function: linear_solve()
// Usage:
// solv = linear_solve(A,b)
// Description:
// Solves the linear system Ax=b. If `A` is square and non-singular the unique solution is returned. If `A` is overdetermined
// the least squares solution is returned. If `A` is underdetermined, the minimal norm solution is returned.
// If `A` is rank deficient or singular then linear_solve returns `[]`. If `b` is a matrix that is compatible with `A`
// then the problem is solved for the matrix valued right hand side and a matrix is returned. Note that if you
// want to solve Ax=b1 and Ax=b2 that you need to form the matrix `transpose([b1,b2])` for the right hand side and then
// transpose the returned value.
function linear_solve(A,b,pivot=true) =
assert(is_matrix(A), "Input should be a matrix.")
let(
m = len(A),
n = len(A[0])
)
assert(is_vector(b,m) || is_matrix(b,m),"Invalid right hand side or incompatible with the matrix")
let (
qr = m<n? qr_factor(transpose(A),pivot) : qr_factor(A,pivot),
maxdim = max(n,m),
mindim = min(n,m),
Q = submatrix(qr[0],[0:maxdim-1], [0:mindim-1]),
R = submatrix(qr[1],[0:mindim-1], [0:mindim-1]),
P = qr[2],
zeros = [for(i=[0:mindim-1]) if (approx(R[i][i],0)) i]
)
zeros != [] ? [] :
m<n ? Q*back_substitute(R,transpose(P)*b,transpose=true) // Too messy to avoid input checks here
: P*_back_substitute(R, transpose(Q)*b); // Calling internal version skips input checks
// Function: matrix_inverse()
// Usage:
// mat = matrix_inverse(A)
// Description:
// Compute the matrix inverse of the square matrix `A`. If `A` is singular, returns `undef`.
// Note that if you just want to solve a linear system of equations you should NOT use this function.
// Instead use {{linear_solve()}}, or use {{qr_factor()}}. The computation
// will be faster and more accurate.
function matrix_inverse(A) =
assert(is_matrix(A) && len(A)==len(A[0]),"Input to matrix_inverse() must be a square matrix")
linear_solve(A,ident(len(A)));
// Function: rot_inverse()
// Usage:
// B = rot_inverse(A)
// Description:
// Inverts a 2d (3x3) or 3d (4x4) rotation matrix. The matrix can be a rotation around any center,
// so it may include a translation.
function rot_inverse(T) =
assert(is_matrix(T,square=true),"Matrix must be square")
let( n = len(T))
assert(n==3 || n==4, "Matrix must be 3x3 or 4x4")
let(
rotpart = [for(i=[0:n-2]) [for(j=[0:n-2]) T[j][i]]],
transpart = [for(row=[0:n-2]) T[row][n-1]]
)
assert(approx(determinant(T),1),"Matrix is not a rotation")
concat(hstack(rotpart, -rotpart*transpart),[[for(i=[2:n]) 0, 1]]);
// Function: null_space()
// Usage:
// x = null_space(A)
// Description:
// Returns an orthonormal basis for the null space of `A`, namely the vectors {x} such that Ax=0.
// If the null space is just the origin then returns an empty list.
function null_space(A,eps=1e-12) =
assert(is_matrix(A))
let(
Q_R = qr_factor(transpose(A),pivot=true),
R = Q_R[1],
zrow = [for(i=idx(R)) if (all_zero(R[i],eps)) i]
)
len(zrow)==0 ? [] :
transpose(columns(Q_R[0],zrow));
// Function: qr_factor()
// Usage:
// qr = qr_factor(A,[pivot]);
// Description:
// Calculates the QR factorization of the input matrix A and returns it as the list [Q,R,P]. This factorization can be
// used to solve linear systems of equations. The factorization is `A = Q*R*transpose(P)`. If pivot is false (the default)
// then P is the identity matrix and A = Q*R. If pivot is true then column pivoting results in an R matrix where the diagonal
// is non-decreasing. The use of pivoting is supposed to increase accuracy for poorly conditioned problems, and is necessary
// for rank estimation or computation of the null space, but it may be slower.
function qr_factor(A, pivot=false) =
assert(is_matrix(A), "Input must be a matrix." )
let(
m = len(A),
n = len(A[0])
)
let(
qr = _qr_factor(A, Q=ident(m),P=ident(n), pivot=pivot, column=0, m = m, n=n),
Rzero = let( R = qr[1]) [
for(i=[0:m-1]) [
let( ri = R[i] )
for(j=[0:n-1]) i>j ? 0 : ri[j]
]
]
) [qr[0], Rzero, qr[2]];
function _qr_factor(A,Q,P, pivot, column, m, n) =
column >= min(m-1,n) ? [Q,A,P] :
let(
swap = !pivot ? 1
: _swap_matrix(n,column,column+max_index([for(i=[column:n-1]) sqr([for(j=[column:m-1]) A[j][i]])])),
A = pivot ? A*swap : A,
x = [for(i=[column:1:m-1]) A[i][column]],
alpha = (x[0]<=0 ? 1 : -1) * norm(x),
u = x - concat([alpha],repeat(0,m-1)),
v = alpha==0 ? u : u / norm(u),
Qc = ident(len(x)) - 2*outer_product(v,v),
Qf = [for(i=[0:m-1]) [for(j=[0:m-1]) i<column || j<column ? (i==j ? 1 : 0) : Qc[i-column][j-column]]]
)
_qr_factor(Qf*A, Q*Qf, P*swap, pivot, column+1, m, n);
// Produces an n x n matrix that swaps column i and j (when multiplied on the right)
function _swap_matrix(n,i,j) =
assert(i<n && j<n && i>=0 && j>=0, "Swap indices out of bounds")
[for(y=[0:n-1]) [for (x=[0:n-1])
x==i ? (y==j ? 1 : 0)
: x==j ? (y==i ? 1 : 0)
: x==y ? 1 : 0]];
// Function: back_substitute()
// Usage:
// x = back_substitute(R, b, [transpose]);
// Description:
// Solves the problem Rx=b where R is an upper triangular square matrix. The lower triangular entries of R are
// ignored. If transpose==true then instead solve transpose(R)*x=b.
// You can supply a compatible matrix b and it will produce the solution for every column of b. Note that if you want to
// solve Rx=b1 and Rx=b2 you must set b to transpose([b1,b2]) and then take the transpose of the result. If the matrix
// is singular (e.g. has a zero on the diagonal) then it returns [].
function back_substitute(R, b, transpose = false) =
assert(is_matrix(R, square=true))
let(n=len(R))
assert(is_vector(b,n) || is_matrix(b,n),str("R and b are not compatible in back_substitute ",n, len(b)))
transpose
? reverse(_back_substitute(transpose(R, reverse=true), reverse(b)))
: _back_substitute(R,b);
function _back_substitute(R, b, x=[]) =
let(n=len(R))
len(x) == n ? x
: let(ind = n - len(x) - 1)
R[ind][ind] == 0 ? []
: let(
newvalue = len(x)==0
? b[ind]/R[ind][ind]
: (b[ind]-list_tail(R[ind],ind+1) * x)/R[ind][ind]
)
_back_substitute(R, b, concat([newvalue],x));
// Function: cholesky()
// Usage:
// L = cholesky(A);
// Description:
// Compute the cholesky factor, L, of the symmetric positive definite matrix A.
// The matrix L is lower triangular and `L * transpose(L) = A`. If the A is
// not symmetric then an error is displayed. If the matrix is symmetric but
// not positive definite then undef is returned.
function cholesky(A) =
assert(is_matrix(A,square=true),"A must be a square matrix")
assert(is_matrix_symmetric(A),"Cholesky factorization requires a symmetric matrix")
_cholesky(A,ident(len(A)), len(A));
function _cholesky(A,L,n) =
A[0][0]<0 ? undef : // Matrix not positive definite
len(A) == 1 ? submatrix_set(L,[[sqrt(A[0][0])]], n-1,n-1):
let(
i = n+1-len(A)
)
let(
sqrtAii = sqrt(A[0][0]),
Lnext = [for(j=[0:n-1])
[for(k=[0:n-1])
j<i-1 || k<i-1 ? (j==k ? 1 : 0)
: j==i-1 && k==i-1 ? sqrtAii
: j==i-1 ? 0
: k==i-1 ? A[j-(i-1)][0]/sqrtAii
: j==k ? 1 : 0]],
Anext = submatrix(A,[1:n-1], [1:n-1]) - outer_product(list_tail(A[0]), list_tail(A[0]))/A[0][0]
)
_cholesky(Anext,L*Lnext,n);
// Function: det2()
// Usage:
// d = det2(M);
// Description:
// Optimized function that returns the determinant for the given 2x2 square matrix.
// Arguments:
// M = The 2x2 square matrix to get the determinant of.
// Example:
// M = [ [6,-2], [1,8] ];
// det = det2(M); // Returns: 50
function det2(M) =
assert(is_matrix(M,2,2), "Matrix must be 2x2.")
M[0][0] * M[1][1] - M[0][1]*M[1][0];
// Function: det3()
// Usage:
// d = det3(M);
// Description:
// Optimized function that returns the determinant for the given 3x3 square matrix.
// Arguments:
// M = The 3x3 square matrix to get the determinant of.
// Example:
// M = [ [6,4,-2], [1,-2,8], [1,5,7] ];
// det = det3(M); // Returns: -334
function det3(M) =
assert(is_matrix(M,3,3), "Matrix must be 3x3.")
M[0][0] * (M[1][1]*M[2][2]-M[2][1]*M[1][2]) -
M[1][0] * (M[0][1]*M[2][2]-M[2][1]*M[0][2]) +
M[2][0] * (M[0][1]*M[1][2]-M[1][1]*M[0][2]);
// Function: determinant()
// Usage:
// d = determinant(M);
// Description:
// Returns the determinant for the given square matrix.
// Arguments:
// M = The NxN square matrix to get the determinant of.
// Example:
// M = [ [6,4,-2,9], [1,-2,8,3], [1,5,7,6], [4,2,5,1] ];
// det = determinant(M); // Returns: 2267
function determinant(M) =
assert(is_matrix(M, square=true), "Input should be a square matrix." )
len(M)==1? M[0][0] :
len(M)==2? det2(M) :
len(M)==3? det3(M) :
sum(
[for (col=[0:1:len(M)-1])
((col%2==0)? 1 : -1) *
M[col][0] *
determinant(
[for (r=[1:1:len(M)-1])
[for (c=[0:1:len(M)-1])
if (c!=col) M[c][r]
]
]
)
]
);
// Function: is_matrix()
// Usage:
// test = is_matrix(A, [m], [n], [square])
// Description:
// Returns true if A is a numeric matrix of height m and width n. If m or n
// are omitted or set to undef then true is returned for any positive dimension.
// Arguments:
// A = The matrix to test.
// m = Is given, requires the matrix to have the given height.
// n = Is given, requires the matrix to have the given width.
// square = If true, requires the matrix to have a width equal to its height. Default: false
function is_matrix(A,m,n,square=false) =
is_list(A)
&& (( is_undef(m) && len(A) ) || len(A)==m)
&& (!square || len(A) == len(A[0]))
&& is_vector(A[0],n)
&& is_consistent(A);
// Function: norm_fro()
// Usage:
// norm_fro(A)
// Description:
// Computes frobenius norm of input matrix. The frobenius norm is the square root of the sum of the
// squares of all of the entries of the matrix. On vectors it is the same as the usual 2-norm.
// This is an easily computed norm that is convenient for comparing two matrices.
function norm_fro(A) =
assert(is_matrix(A) || is_vector(A))
norm(flatten(A));
// Function: matrix_trace()
// Usage:
// matrix_trace(M)
// Description:
// Computes the trace of a square matrix, the sum of the entries on the diagonal.
function matrix_trace(M) =
assert(is_matrix(M,square=true), "Input to trace must be a square matrix")
[for(i=[0:1:len(M)-1])1] * [for(i=[0:1:len(M)-1]) M[i][i]];
// Section: Comparisons and Logic