Compare commits
No commits in common. "bd1eee136662ef638628bca9907ae64f1b7538f9" and "a63e6246acf80346932449d1145d05211bb00aa4" have entirely different histories.
bd1eee1366
...
a63e6246ac
|
@ -28,7 +28,7 @@ combine the most useful features of each library.
|
|||
## Documentation
|
||||
|
||||
Most of the exported utility functions have docstrings. There is also
|
||||
a (currently unfinished) guide on how to use NumIdr [here](docs/Contents.md).
|
||||
a (currently unfinished) guide on how to use NumIdr [here](docs/Intro.md).
|
||||
|
||||
## Usage
|
||||
|
||||
|
|
|
@ -1,19 +0,0 @@
|
|||
# NumIdr Guide
|
||||
|
||||
Welcome to NumIdr's guide! This guide will serve as a basic introduction to the core ideas of NumIdr, as well as a reference for common functions and utilities that you might otherwise miss.
|
||||
|
||||
If you're familiar with the Python library [NumPy](https://numpy.org/) or anything based on it, then a lot of the early content will be familiar, as NumIdr is based on the same array type.
|
||||
|
||||
> [!NOTE]
|
||||
> This guide assumes a basic understanding of linear algebra, and it will not explain any linear algebra concepts necessary for understanding NumIdr's structure.
|
||||
|
||||
### Primary Features
|
||||
|
||||
1. [Fundamental Data Types](DataTypes.md)
|
||||
2. [Basic Operations on Arrays](Operations.md)
|
||||
3. [Working with Vectors and Matrices](VectorsMatrices.md)
|
||||
4. [Transforms](Transforms.md)
|
||||
|
||||
### Advanced Concepts
|
||||
|
||||
5. Array Representations (coming soon!)
|
|
@ -49,6 +49,9 @@ When determining the index of a value inside the array, the order of the indices
|
|||
> (as in "multi-dimensional array" in the previous section), or to the lengths of its
|
||||
> axes. Conventionally, NumIdr reserves "dimension" for the second meaning, and uses
|
||||
> "rank" for the first meaning.
|
||||
>
|
||||
> This guide has ignored this convention until now to be more understandable to newcomers,
|
||||
> but will follow it from this point onward.
|
||||
|
||||
## Types of Arrays
|
||||
|
||||
|
@ -76,6 +79,8 @@ Vector n a = Array [n] a
|
|||
|
||||
A vector's type signature and stored data is effectively identical to that of the standard library type `Vect`, whose elements are confusingly also called "vectors"; we often refer to those as "vects" to differentiate.
|
||||
|
||||
Indices are typically written as lists of integers, but for vectors it is occasionally acceptable to write the single index number without putting it inside a list. This is mostly the case for indexing, where each indexing function has an alternate definition specifically for vectors.
|
||||
|
||||
### Matrices
|
||||
|
||||
As mentioned before, a matrix is a rank-2 array:
|
||||
|
@ -116,18 +121,36 @@ These are useful for clarity when working with both homogeneous and non-homogene
|
|||
|
||||
### Transforms
|
||||
|
||||
A transform is a wrapper type for a square matrix with certain properties that can be used to transform points in space.
|
||||
A transform is a wrapper type for a matrix with certain properties that can be used to transform points in space.
|
||||
|
||||
```idris
|
||||
Transform : (ty : TransType) -> (n : Nat) -> (a : Type) -> Type
|
||||
```
|
||||
|
||||
The `TransType` parameter dictates what kind of transform it is. More information on transforms and their operations can be found in the [Transforms](Transforms.md) chapter.
|
||||
The `TransType` parameter dictates what kind of transform it is. These eight options are currently available:
|
||||
|
||||
**Linear Types:**
|
||||
- `Trivial` (always the identity transformation)
|
||||
- `Rotation`
|
||||
- `Orthonormal` (rotation + reflection)
|
||||
- `Linear`
|
||||
|
||||
**Affine Types:**
|
||||
- `Translation`
|
||||
- `Rigid` (rotation + translation)
|
||||
- `Isometry` (rotation + reflection + translation)
|
||||
- `Affine`
|
||||
|
||||
The `TransType` value is obtained by prepending a capital T to these names. For example, an isometry may have the type `Isometry 3 Double`, which is an alias for `Transform TIsometry 3 Double`.
|
||||
|
||||
#### The Point Type
|
||||
|
||||
Transforms behave differently from regular matrices when applied to a vector. When a non-linear transform is used, the transform is first linearized, so that vectors only have linear transformations applied to them. This is not a bug.
|
||||
|
||||
In order to properly apply these transforms, the `Point` type must be used, which is a wrapper around the `Vector` type that supports these transforms. This separation between points and vectors is intended to make working with affine transformations more convenient, as it mirrors the separation between points and vectors in affine algebra.
|
||||
|
||||
### Permutations
|
||||
|
||||
The type `Permutation n` represents a permutation of `n` elements. Permutations are mostly used internally for various algorithms, but they are also an input in some operations, such as permuting the rows or columns of a matrix.
|
||||
The type `Permutation n` represents a permutation of `n` elements. Permutations are mostly used internally for various algorithms, but they are also an input in various operations, such as those that permute the axes of an array.
|
||||
|
||||
Permutations can be composed using `(*.)`, and a permutation can be converted into a matrix using `permuteM`.
|
||||
|
||||
[Contents](Contents.md) | [Next](Operations.md)
|
||||
[Contents](Intro.md) | [Next](Operations.md)
|
||||
|
|
16
docs/Intro.md
Normal file
16
docs/Intro.md
Normal file
|
@ -0,0 +1,16 @@
|
|||
# NumIdr Guide
|
||||
|
||||
Welcome to NumIdr's guide!
|
||||
|
||||
If you're familiar with the Python library [NumPy](https://numpy.org/) or anything based on it, then a lot of the early content will be familiar, as NumIdr is based on the same array type.
|
||||
|
||||
#### Tutorial
|
||||
|
||||
1. [Fundamental Data Types](DataTypes.md)
|
||||
2. [Basic Operations on Arrays](Operations.md)
|
||||
3. Working with Vectors and Matrices
|
||||
4. Transforms
|
||||
|
||||
#### Advanced Concepts
|
||||
|
||||
5. Array Representations
|
|
@ -1,10 +1,5 @@
|
|||
# Basic Operations on Arrays
|
||||
|
||||
> [!CAUTION]
|
||||
> Arrays and their associated functions are not intended to be evaluated at compile-time. If you try to compute an array in the REPL, you will not get the output you expect!
|
||||
>
|
||||
> If you really need to use the REPL to test array code, use `:exec`.
|
||||
|
||||
## Constructing Arrays
|
||||
|
||||
The most important array constructor is `array`, which returns an array of the specified values:
|
||||
|
@ -13,7 +8,7 @@ The most important array constructor is `array`, which returns an array of the s
|
|||
array [[1, 2, 3], [4, 5, 6]]
|
||||
```
|
||||
|
||||
Scalars, vectors and matrices have their own constructors, used in exactly the same way (`scalar`, `vector`, and `matrix`). These should be used instead of `array` wherever possible, as they provide more information to the type-checker.
|
||||
Scalars, vectors and matrices have their own constructors, used in exactly the same way (`scalar`, `vector`, and `matrix`). These should be used instead of `array` whereever possible, as they provide more information to the type-checker.
|
||||
|
||||
There are also a few other, more basic constructors for convenience.
|
||||
|
||||
|
@ -32,13 +27,14 @@ ones [2, 2, 3]
|
|||
|
||||
There are a few simple functions for accessing basic properties of arrays: `shape` and `rank`, which are self-explanatory, and `size`, which returns the total number of elements in the array.
|
||||
|
||||
> [!TIP]
|
||||
> The `shape` accessor is sufficient for most uses, but it can cause problems with the type-checker, as for an array `arr : Array s a` the type checker does not know that `shape arr` and `s` are equal. To solve this problem, you can use the `ShapeView`:
|
||||
>
|
||||
> ```idris
|
||||
> example {s} arr with (viewShape arr)
|
||||
> _ | Shape s = ...
|
||||
> ```
|
||||
The `shape` accessor is sufficient for most uses, but it can cause problems with the type-checker, as for an array `arr : Array s a` the type checker does not know that `shape arr` and `s` are equal. To solve this problem, a view for accessing the shape is provided:
|
||||
|
||||
```idris
|
||||
example {s} arr with (viewShape arr)
|
||||
_ | Shape s = ...
|
||||
```
|
||||
|
||||
This will fully placate the type-checker, as the `s` returned by the view is proven to be equal to the shape of the array.
|
||||
|
||||
## Indexing Arrays
|
||||
|
||||
|
@ -66,7 +62,7 @@ Not all combinations of these categories are defined by the library. Here are th
|
|||
| **Update** | `indexUpdate` | `indexUpdateRange` | `indexUpdateNB` | | | |
|
||||
| **Set** | `indexSet` | `indexSetRange` | `indexSetNB` | | | |
|
||||
|
||||
The accessor functions have operator forms for convenience.
|
||||
The accessor functions have operator forms for convenience, also specified within the table.
|
||||
|
||||
### Specifying Coordinates
|
||||
|
||||
|
@ -81,23 +77,20 @@ arr !! [1, 0]
|
|||
|
||||
With ranged indexing, a sub-array of the original array is accessed or modified. This sub-array is given by a list of _range specifiers_, one for each axis, which can be one of the following:
|
||||
|
||||
- `Bounds x y` - Every index from `x` (inclusive) to `y` (exclusive)
|
||||
- `Bounds x y` - Every index from `x` to `y`
|
||||
- `StartBound x` - Every index from `x` to the end of the axis
|
||||
- `EndBound y` - Every index from the start of the axis to `y`
|
||||
- `One i`, `One' i` - The single index `i`
|
||||
- `All` - Every index in the axis
|
||||
- `Indices xs` - A list of indices in the axis
|
||||
- `Filter p` - A predicate specifying which indices to access or modify
|
||||
|
||||
The range specifier `One` is special, as it is the only range specifier that decreases the rank of the resulting array. For example, the following matrix access returns the first column as a vector due to `One` being used:
|
||||
There is also an extra specifier, `One i`, that selects exactly one index of the axis. It is only usable in safe indexing, and it is notable for being the only range specifier that decreases the rank of the sub-array. For example, the following matrix access returns the first column as a vector due to `One` being used:
|
||||
|
||||
```idris
|
||||
matrix [[2, 0], [1, 2]] !!.. [All, One 0]
|
||||
== vector [2, 1]
|
||||
```
|
||||
|
||||
`One'` does not decrease the rank of the result in this way.
|
||||
|
||||
## Standard Interface Methods
|
||||
|
||||
The array type implements a number of standard interfaces. Most of these implementations are unremarkable, but a few have caveats that are worth noting.
|
||||
|
@ -113,8 +106,8 @@ matrix [[1, 1], [2, 5]] + matrix [[2, 3], [-1, 3]]
|
|||
|
||||
This elementwise behavior holds for `(+)`, `(*)`, `(-)`, `(/)`, and `(<+>)`. **`(*)` is not matrix multiplication!** For the generalized multiplication operator, which includes matrix multiplication, see [the next chapter](VectorsMatrices.md).
|
||||
|
||||
> [!WARNING]
|
||||
> Due to unfortunate restrictions in Idris's standard `Num` interface, `(+)` and `(*)` can only be used when the array's shape is available at run-time. If this is not the case, you must use `zipwith (+)` or `zipWith (*)` instead.
|
||||
> [!NOTE]
|
||||
> Due to unfortunate restrictions in Idris's standard `Num` interface, the addition and multiplication operations can only be used when the array's shape is available at run-time. If this is not the case, you must use `zipwith (+)` or `zipWith (*)` instead.
|
||||
|
||||
### `Foldable` and `Traversable`
|
||||
|
||||
|
@ -124,7 +117,7 @@ When folding or traversing the elements of an array, these elements are ordered
|
|||
|
||||
### Concatenation and Stacking
|
||||
|
||||
Two arrays can be concatenated along an axis (`concat`), so long as all other axes have the same dimensions. Two matrices being concatenated along the row axis requires that they must have the same number of columns.
|
||||
Two arrays can be concatenated along an axis, so long as all other axes have the same dimensions. Two matrices being concatenated along the row axis requires that they must have the same number of columns.
|
||||
|
||||
```idris
|
||||
-- 0 is the first axis i.e. the row axis
|
||||
|
@ -135,7 +128,7 @@ concat 0 (matrix [[1, 2], [3, 4]]) (matrix [[5, 6], [7, 8]])
|
|||
[7, 8]]
|
||||
```
|
||||
|
||||
Stacking (`stack`) is similar to concatenation, but slightly different. Stacking combines arrays with the exact same shape into a single array that is one rank higher. For example, vectors can be stacked along the row axis to obtain a matrix whose rows are the original vectors.
|
||||
Stacking is similar to concatenation, but slightly different. Stacking combines arrays with the exact same shape into a single array that is one rank higher. For example, vectors can be stacked along the row axis to obtain a matrix whose rows are the original vectors.
|
||||
|
||||
```idris
|
||||
stack 0 [vector [1, 2], vector [3, 4]]
|
||||
|
@ -156,7 +149,7 @@ reshape [3, 2] (vector [1, 2, 3, 4, 5, 6])
|
|||
[5, 6]]
|
||||
```
|
||||
|
||||
Arrays can also be resized with `resize`, which changes their shape while keeping every element at the same index. A default element must be provided to fill any indices that did not exist in the original array.
|
||||
Arrays can also be resized, which changes their shape while keeping every element at the same index. A default element must be provided to fill any indices that did not exist in the original array.
|
||||
|
||||
```idris
|
||||
resize [2, 4] 10 (matrix [[1, 2],
|
||||
|
@ -170,10 +163,8 @@ Instead of the `resize` function, one can also use the `resizeLTE` function, whi
|
|||
|
||||
### Transpose
|
||||
|
||||
The `transpose` function reverses the axis order of an array: For `arr : Array [3,2,4] Int`, we have `transpose arr : Array [4,2,3] Int`. For matrices, this corresponds to the usual definition of switching rows and columns. There is also a postfix form `(.T)`.
|
||||
The `transpose` function reverses the axis order of an array. For matrices, this corresponds to the usual definition of switching rows and columns. There is also a postfix form `(.T)`.
|
||||
|
||||
For more fine-grained control when rearranging arrays, there are the `swapAxes` and `permuteAxes` functions, where the first swaps only two axes and the second takes an arbitrary [permutation](DataTypes.md#Permutations). There are also `swapInAxis` and `permuteInAxis`, which permute inside an axis, e.g. swapping rows or columns in a matrix.
|
||||
|
||||
Like with concatenation and stacking, the swap and permute functions have forms specific to vectors and matrices: `swapCoords` and `permuteCoords` for vectors, and `swapRows`, `permuteRows`, `swapColumns`, and `permuteColumns` for matrices.
|
||||
|
||||
[Previous](DataTypes.md) | [Contents](Contents.md) | [Next](VectorsMatrices.md)
|
||||
[Previous](DataTypes.md) | [Contents](Intro.md) | [Next](VectorsMatrices.md)
|
||||
|
|
|
@ -1,108 +1,4 @@
|
|||
# Transforms
|
||||
|
||||
In code that works with linear algebra, it is common to use matrices as _transformations_, functions that take in a vector and output a new vector. These matrices can often be divided into categories based on the operations they perform, such as a rotation matrix or an affine transformation matrix.
|
||||
|
||||
The `Transform` wrapper exists to encode these differing properties on the type level, as well as to provide extra utilities for working with matrices in this fashion.
|
||||
|
||||
## Types of Transforms
|
||||
|
||||
`Transform` has the following type signature:
|
||||
|
||||
```idris
|
||||
Transform : (ty : TransType) -> (n : Nat) -> (a : Type) -> Type
|
||||
```
|
||||
|
||||
A `Transform ty n a` is a wrapper over an `HMatrix' n a`. The `ty` parameter is the transform type, which dictates what properties the transform has. These eight options are currently available:
|
||||
|
||||
**Affine Types:**
|
||||
- `TAffine`
|
||||
- `TIsometry` (rotation + reflection + translation)
|
||||
- `TRigid` (rotation + translation)
|
||||
- `TTranslation`
|
||||
|
||||
**Linear Types:**
|
||||
- `TLinear`
|
||||
- `TOrthonormal` (rotation + reflection)
|
||||
- `TRotation`
|
||||
- `TTrivial` (always the `identity`)
|
||||
|
||||
The capital T at the beginning of each of these names identifies it as a `TransType` value. To make working with transforms smoother, NumIdr provides synonyms for transforms of each type: for example, `Isometry n a` is a synonym for `Transform TIsometry n a`.
|
||||
|
||||
### Linear and Affine
|
||||
|
||||
Transform types are divided into linear and affine types. Linear transforms must preserve the origin point, whereas affine transforms do not have this restriction.
|
||||
|
||||
Linear and affine transform types are in a one-to-one correspondence: a linear transform can be converted to and from an affine transform by adding or removing a translation component.
|
||||
|
||||
```
|
||||
Linear <-> Affine
|
||||
Orthonormal <-> Isometry
|
||||
Rotation <-> Rigid
|
||||
Trivial <-> Translation
|
||||
```
|
||||
|
||||
The `setTranslation` and `linearize` functions perform these conversions.
|
||||
|
||||
For simplicity, both categories of transform are wrappers over homogeneous matrices, even though linear transforms could be represented by non-homogeneous matrices.
|
||||
|
||||
### Transform Type Casting
|
||||
|
||||
Some transform types can be cast into other types. For example, a `Rotation` can be cast into an `Orthonormal`, as all rotation matrices are orthonormal.
|
||||
|
||||
```idris
|
||||
rt : Rotation 3 Double
|
||||
|
||||
cast rt : Orthonormal 3 Double
|
||||
```
|
||||
|
||||
In the diagram from the previous section, lower types can be cast into types higher up. Each linear type (on the left) can also be cast into the corresponding affine type (on the right).
|
||||
|
||||
## Constructing Transforms
|
||||
|
||||
There are multiple ways to construct transforms, either by wrapping a matrix or directly through constructor functions.
|
||||
|
||||
For each transform type, `fromHMatrix` can be used to test if a homogeneous matrix satisfies the right properties, and converts it into a transform if it does. The `Rotation`, `Orthonormal` and `Linear` types also have `fromMatrix` for non-homogeneous matrices.
|
||||
|
||||
> [!WARNING]
|
||||
> The `fromHMatrix` and `fromMatrix` constructors use exact equality comparisons when testing matrices, which can be a problem if your element type is `Double` or a similar inexact number type. To fix this issue, NumIdr provides a `WithEpsilon` named implementation that defines equality approximately and allows you to specify an epsilon value.
|
||||
>
|
||||
> ```idris
|
||||
> fromHMatrix @{WithEpsilon 1.0e-6} mat
|
||||
> ```
|
||||
|
||||
There are also direct constructors for transforms, which are often more convenient as they do not have the possibility of failing. There are too many of these constructors to exhaustively list here, so I encourage you to look through the functions in the `Data.NumIdr.Transform.*` modules to see what is available.
|
||||
|
||||
## Multiplication with Transforms
|
||||
|
||||
Like most objects in NumIdr, transforms multiply with the generalized multiplication operator `(*.)`, and `identity` and `inverse` can also be used with transforms.
|
||||
|
||||
> [!NOTE]
|
||||
> There is no `tryInverse` function for transforms. This is because all transforms are required to be invertible, so there isn't any need for it; you can safely use `inverse` in all circumstances.
|
||||
|
||||
Transforms of any types can be multiplied. When two transforms of different types are multiplied, the resulting transform type is determined by taking the most specific type that both original types can be cast to. For example, an `Orthonormal` transform multiplied by a `Translation` returns an `Isometry`.
|
||||
|
||||
### The Point Type
|
||||
|
||||
Transforms behave differently from regular matrices when applied to vectors. When an affine transform is applied in this way, it is first linearized, so that vectors only have linear transforms applied to them. **This is not a bug!**
|
||||
|
||||
In order to properly apply affine transforms, the `Point` type must be used, which is a wrapper around the `Vector` type that supports these transforms. A point can be constructed with the `point` function, which is used exactly the same as the `vector` constructor.
|
||||
|
||||
```idris
|
||||
point [4, 3, 6]
|
||||
```
|
||||
|
||||
Points support most basic operations that vectors do, including indexing operations and standard library methods. However, a point cannot be added to another point. Instead, a vector must be added to a point:
|
||||
|
||||
```idris
|
||||
(+.) : Vector n a -> Point n a -> Point n a
|
||||
|
||||
(.+) : Point n a -> Vector n a -> Point n a
|
||||
|
||||
(-.) : Point n a -> Point n a -> Vector n a
|
||||
```
|
||||
|
||||
To remember the distinction between the two addition operators, the dot is always on the side of the point, not the vector.
|
||||
|
||||
This separation between points and vectors is intended to make working with affine transformations more convenient, as it mirrors the separation between points and vectors in affine algebra. These may feel like arbitrary restrictions, but you might be surprised by how convenient they are to work with!
|
||||
|
||||
[Previous](VectorsMatrices.md) | [Contents](Contents.md)
|
||||
[Previous](VectorsMatrices.md) | [Contents](Intro.md)
|
||||
|
|
|
@ -1,114 +1,4 @@
|
|||
# Working with Vectors and Matrices
|
||||
|
||||
As linear algebra is one of the main concerns of NumIdr, most of its provided functions are dedicated to vectors (rank-1 arrays) and matrices (rank-2 arrays).
|
||||
|
||||
## The Generalized Multiplication Operator
|
||||
|
||||
A linear algebra library wouldn't be very useful without matrix multiplication! Since `(*)` is already used for element-wise multiplication, NumIdr defines a new interface `Mult` that can accept and return values of different types:
|
||||
|
||||
```idris
|
||||
interface Mult a b c where
|
||||
(*.) : a -> b -> c
|
||||
```
|
||||
|
||||
The generalized multiplication operator `(*.)` covers matrix multiplication, scalar-vector multiplication, and any other linear algebra operation that's vaguely multiplication-like.
|
||||
|
||||
## Vectors
|
||||
|
||||
### Algebraic Operations
|
||||
|
||||
Vectors can be added together with element-wise addition `(+)`. Scalar-vector multiplication is done with the generalized multiplication operator `(*.)`.
|
||||
|
||||
```idris
|
||||
2 *. (vector [1, 1] + vector [2, 3])
|
||||
== vector [6, 8]
|
||||
```
|
||||
|
||||
A few other basic linear algebra operations are available:
|
||||
|
||||
- `dot`, The dot product
|
||||
- `cross`, The cross product
|
||||
- `perp`, The perpendicular product (sometimes called the 2D cross product)
|
||||
- `triple`, The scalar triple product
|
||||
|
||||
### Indexing
|
||||
|
||||
NumIdr provides special versions of `index` and `indexNB` and their infix forms `(!!)` and `(!?)` for use with vectors. These take a single numeric index instead of a list.
|
||||
|
||||
```idris
|
||||
Vector.index 2 v == index [2] v
|
||||
|
||||
v !! 2 == v !! [2]
|
||||
```
|
||||
|
||||
For convenience, when working with two- or three-dimensional vectors, there are postfix accessors `(.x)`, `(.y)`, and `(.z)`:
|
||||
|
||||
```idris
|
||||
v = vector [5, 6, 2]
|
||||
|
||||
v.x == 5
|
||||
v.y == 6
|
||||
v.z == 2
|
||||
```
|
||||
|
||||
### Other Operations
|
||||
|
||||
- `toVect` - Convert a vector into a `Vect`
|
||||
- `dim` - Returns the vector's length
|
||||
- `(++)` - Concatenate two vectors
|
||||
|
||||
## Matrices
|
||||
|
||||
### Arithmetic Operations
|
||||
|
||||
Like vectors, matrices can be added together using `(+)`. Matrix multiplication, as well as matrix-vector and matrix-scalar multiplication, are performed using `(*.)`.
|
||||
|
||||
For the purposes of working with matrices and matrix-like objects, the sub-interfaces `MultMonoid` and `MultGroup` are defined:
|
||||
|
||||
```idris
|
||||
interface Mult' a => MultMonoid a where
|
||||
identity : a
|
||||
|
||||
interface MultMonoid a => MultGroup a where
|
||||
inverse : a -> a
|
||||
```
|
||||
|
||||
The `identity` function returns an identity matrix, and `inverse` calculates a matrix's inverse. Note that `inverse` cannot tell you if an inverse of your matrix does not exist; if you want to handle that possibility, use `tryInverse` instead.
|
||||
|
||||
```idris
|
||||
tryInverse : FieldCmp a => Matrix' n a -> Maybe (Matrix' n a)
|
||||
```
|
||||
|
||||
You can also use the `invertible` predicate to test if a matrix has an inverse.
|
||||
|
||||
#### LU and LUP Decomposition
|
||||
|
||||
The functions `decompLU` and `decompLUP` compute LU and LUP decomposition on a matrix.
|
||||
|
||||
```idris
|
||||
decompLU : Field a => (mat : Matrix m n a) -> Maybe (DecompLU mat)
|
||||
|
||||
decompLUP : FieldCmp a => (mat : Matrix m n a) -> DecompLUP mat
|
||||
```
|
||||
|
||||
`DecompLU` and `DecompLUP` are record types holding the results of the corresponding decomposition. The accessors `lower`, `upper` and `permute` can be applied to get each component of the decomposition; `lower` and `upper` return matrices, and `permute` returns a `Permutation` value.
|
||||
|
||||
#### Other Algebraic Operations
|
||||
|
||||
- `trace` - The sum of the matrix's diagonal
|
||||
- `outer` - The matrix-valued outer product (or tensor product) of two vectors
|
||||
- `det` - Determinant of the matrix
|
||||
- `solve` - Apply an inverse matrix to a vector, useful for solving linear equations
|
||||
|
||||
> [!TIP]
|
||||
> The `det` and `solve` operations require computing an LUP decomposition, which can be expensive. The variants `detWithLUP` and `solveWithLUP` allow an existing LUP decomposition to be passed in, which can make your code more efficient.
|
||||
|
||||
### Indexing
|
||||
|
||||
Aside from the usual array indexing functions, there are a few functions specialized to matrix indexing:
|
||||
|
||||
- `getRow` and `getColumn` - Returns a specific row or column of the matrix
|
||||
- `diagonal` - Returns the diagonal elements of the matrix as a vector
|
||||
- `minor` - Removes a single row and column from the matrix
|
||||
|
||||
[Previous](Operations.md) | [Contents](Contents.md) | [Next](Transforms.md)
|
||||
[Previous](Operations.md) | [Contents](Intro.md) | [Next](Transforms.md)
|
||||
|
|
|
@ -1,5 +1,5 @@
|
|||
package numidr
|
||||
version = 0.3.0
|
||||
version = 2.0.0
|
||||
brief = "Linear algebra and data science library"
|
||||
|
||||
authors = "Kiana Sheibani"
|
||||
|
|
|
@ -26,9 +26,6 @@ export
|
|||
rangeLenZ : (x : Nat) -> length (range 0 x) = x
|
||||
rangeLenZ x = rangeLen 0 x `trans` minusZeroRight x
|
||||
|
||||
export %unsafe
|
||||
assertFin : Nat -> Fin n
|
||||
assertFin n = natToFinLt n @{believe_me Oh}
|
||||
|
||||
--------------------------------------------------------------------------------
|
||||
-- Array coordinate types
|
||||
|
@ -66,7 +63,6 @@ namespace Strict
|
|||
public export
|
||||
data CRange : Nat -> Type where
|
||||
One : Fin n -> CRange n
|
||||
One' : Fin n -> CRange n
|
||||
All : CRange n
|
||||
StartBound : Fin (S n) -> CRange n
|
||||
EndBound : Fin (S n) -> CRange n
|
||||
|
@ -83,8 +79,6 @@ namespace Strict
|
|||
namespace NB
|
||||
public export
|
||||
data CRangeNB : Type where
|
||||
One : Nat -> CRangeNB
|
||||
One' : Nat -> CRangeNB
|
||||
All : CRangeNB
|
||||
StartBound : Nat -> CRangeNB
|
||||
EndBound : Nat -> CRangeNB
|
||||
|
@ -134,7 +128,6 @@ namespace Strict
|
|||
public export
|
||||
cRangeToList : {n : Nat} -> CRange n -> Either Nat (List Nat)
|
||||
cRangeToList (One x) = Left (cast x)
|
||||
cRangeToList (One' x) = Right [cast x]
|
||||
cRangeToList All = Right $ range 0 n
|
||||
cRangeToList (StartBound x) = Right $ range (cast x) n
|
||||
cRangeToList (EndBound x) = Right $ range 0 (cast x)
|
||||
|
@ -155,8 +148,8 @@ namespace Strict
|
|||
newShape : {s : _} -> (rs : CoordsRange s) -> Vect (newRank rs) Nat
|
||||
newShape [] = []
|
||||
newShape (r :: rs) with (cRangeToList r)
|
||||
_ | Left _ = newShape rs
|
||||
_ | Right xs = length xs :: newShape rs
|
||||
newShape (r :: rs) | Left _ = newShape rs
|
||||
newShape (r :: rs) | Right xs = length xs :: newShape rs
|
||||
|
||||
|
||||
getNewPos : {s : _} -> (rs : CoordsRange {rk} s) -> Vect rk Nat -> Vect (newRank rs) Nat
|
||||
|
@ -182,14 +175,6 @@ namespace NB
|
|||
validateCRange (d :: s) (r :: rs) = [| validate' d r :: validateCRange s rs |]
|
||||
where
|
||||
validate' : (n : Nat) -> CRangeNB -> Maybe (CRange n)
|
||||
validate' n (One i) =
|
||||
case isLT i n of
|
||||
Yes _ => Just (One (natToFinLT i))
|
||||
_ => Nothing
|
||||
validate' n (One' i) =
|
||||
case isLT i n of
|
||||
Yes _ => Just (One' (natToFinLT i))
|
||||
_ => Nothing
|
||||
validate' n All = Just All
|
||||
validate' n (StartBound x) =
|
||||
case isLTE x n of
|
||||
|
@ -212,44 +197,29 @@ namespace NB
|
|||
export %unsafe
|
||||
assertCRange : (s : Vect rk Nat) -> Vect rk CRangeNB -> CoordsRange s
|
||||
assertCRange [] [] = []
|
||||
assertCRange (d :: s) (r :: rs) = assert' r :: assertCRange s rs
|
||||
assertCRange (d :: s) (r :: rs) = assert' d r :: assertCRange s rs
|
||||
where
|
||||
assert' : forall n. CRangeNB -> CRange n
|
||||
assert' (One i) = One (assertFin i)
|
||||
assert' (One' i) = One' (assertFin i)
|
||||
assert' All = All
|
||||
assert' (StartBound x) = StartBound (assertFin x)
|
||||
assert' (EndBound x) = EndBound (assertFin x)
|
||||
assert' (Bounds x y) = Bounds (assertFin x) (assertFin y)
|
||||
assert' (Indices xs) = Indices (assertFin <$> xs)
|
||||
assert' (Filter f) = Filter (f . finToNat)
|
||||
assert' : (n : Nat) -> CRangeNB -> CRange n
|
||||
assert' n All = All
|
||||
assert' n (StartBound x) = StartBound (believe_me x)
|
||||
assert' n (EndBound x) = EndBound (believe_me x)
|
||||
assert' n (Bounds x y) = Bounds (believe_me x) (believe_me y)
|
||||
assert' n (Indices xs) = Indices (believe_me <$> xs)
|
||||
assert' n (Filter f) = Filter (f . finToNat)
|
||||
|
||||
public export
|
||||
cRangeNBToList : Nat -> CRangeNB -> Either Nat (List Nat)
|
||||
cRangeNBToList s (One i) = Left i
|
||||
cRangeNBToList s (One' i) = Right [i]
|
||||
cRangeNBToList s All = Right $ range 0 s
|
||||
cRangeNBToList s (StartBound x) = Right $ range x s
|
||||
cRangeNBToList s (EndBound x) = Right $ range 0 x
|
||||
cRangeNBToList s (Bounds x y) = Right $ range x y
|
||||
cRangeNBToList s (Indices xs) = Right $ nub xs
|
||||
cRangeNBToList s (Filter p) = Right $ filter p $ range 0 s
|
||||
|
||||
public export
|
||||
newRank : Vect rk Nat -> Vect rk CRangeNB -> Nat
|
||||
newRank _ [] = 0
|
||||
newRank (d :: s) (r :: rs) =
|
||||
case cRangeNBToList d r of
|
||||
Left _ => newRank s rs
|
||||
Right _ => S (newRank s rs)
|
||||
cRangeNBToList : Nat -> CRangeNB -> List Nat
|
||||
cRangeNBToList s All = range 0 s
|
||||
cRangeNBToList s (StartBound x) = range x s
|
||||
cRangeNBToList s (EndBound x) = range 0 x
|
||||
cRangeNBToList s (Bounds x y) = range x y
|
||||
cRangeNBToList s (Indices xs) = nub xs
|
||||
cRangeNBToList s (Filter p) = filter p $ range 0 s
|
||||
|
||||
||| Calculate the new shape given by a coordinate range.
|
||||
public export
|
||||
newShape : (s : Vect rk Nat) -> (is : Vect rk CRangeNB) -> Vect (newRank s is) Nat
|
||||
newShape [] [] = []
|
||||
newShape (d :: s) (r :: rs) with (cRangeNBToList d r)
|
||||
_ | Left _ = newShape s rs
|
||||
_ | Right xs = length xs :: newShape s rs
|
||||
newShape : Vect rk Nat -> Vect rk CRangeNB -> Vect rk Nat
|
||||
newShape = zipWith (length .: cRangeNBToList)
|
||||
|
||||
export
|
||||
getAllCoords' : Vect rk Nat -> List (Vect rk Nat)
|
||||
|
|
|
@ -31,16 +31,16 @@ export
|
|||
FieldCmp Double where
|
||||
abslt = (<) `on` abs
|
||||
|
||||
-- Alternative implementations of `Eq` and `FieldCmp` that compare approximately,
|
||||
-- useful when working with flating point numbers
|
||||
-- Alternative implementations of `Eq` and `FieldCmp` that compare floating
|
||||
-- point numbers approximately, useful when working with transforms
|
||||
namespace Eq
|
||||
export
|
||||
WithEpsilon : (Neg a, Abs a, Ord a) => a -> Eq a
|
||||
WithEpsilon : Double -> Eq Double
|
||||
WithEpsilon ep = MkEq (\x,y => abs (x - y) < ep) (\x,y => abs (x - y) >= ep)
|
||||
|
||||
namespace FieldCmp
|
||||
export
|
||||
WithEpsilon : (Neg a, Fractional a, Abs a, Ord a) => a -> FieldCmp a
|
||||
WithEpsilon : Double -> FieldCmp Double
|
||||
WithEpsilon ep = MkFieldCmp @{WithEpsilon ep} ((<) `on` abs)
|
||||
|
||||
--------------------------------------------------------------------------------
|
||||
|
|
|
@ -568,8 +568,9 @@ solveLowerTri' {n} mat b with (viewShape b)
|
|||
construct [] = []
|
||||
construct {i=S i} (b :: bs) =
|
||||
let xs = construct bs
|
||||
i' = assert_total $ case natToFin i n of Just i' => i'
|
||||
in (b - sum (zipWith (*) xs (reverse $ toVect $ replace {p = flip Array a} (believe_me $ Refl {x=()}) $
|
||||
mat !#.. [One i, EndBound i]))) / mat!#[i,i] :: xs
|
||||
mat !!.. [One i', EndBound (weaken i')]))) / mat!#[i,i] :: xs
|
||||
|
||||
|
||||
solveUpperTri' : Field a => Matrix' n a -> Vector n a -> Vector n a
|
||||
|
@ -580,8 +581,9 @@ solveUpperTri' {n} mat b with (viewShape b)
|
|||
construct _ [] = []
|
||||
construct i (b :: bs) =
|
||||
let xs = construct (S i) bs
|
||||
i' = assert_total $ case natToFin i n of Just i' => i'
|
||||
in (b - sum (zipWith (*) xs (toVect $ replace {p = flip Array a} (believe_me $ Refl {x=()}) $
|
||||
mat !#.. [One i, StartBound (S i)]))) / mat!#[i,i] :: xs
|
||||
mat !!.. [One i', StartBound (FS i')]))) / mat!#[i,i] :: xs
|
||||
|
||||
|
||||
||| Solve a linear equation, assuming the matrix is lower triangular.
|
||||
|
@ -626,7 +628,6 @@ solve mat = solveWithLUP mat (decompLUP mat)
|
|||
--------------------------------------------------------------------------------
|
||||
|
||||
|
||||
||| Determine whether a matrix has an inverse.
|
||||
export
|
||||
invertible : FieldCmp a => Matrix' n a -> Bool
|
||||
invertible {n} mat with (viewShape mat)
|
||||
|
|
|
@ -23,6 +23,10 @@ update : Coords s -> (a -> a) -> Vects s a -> Vects s a
|
|||
update [] f v = f v
|
||||
update (i :: is) f v = updateAt i (update is f) v
|
||||
|
||||
export %unsafe
|
||||
assertFin : Nat -> Fin n
|
||||
assertFin n = natToFinLt n @{believe_me Oh}
|
||||
|
||||
export
|
||||
indexRange : {s : _} -> (rs : CoordsRange s) -> Vects s a -> Vects (newShape rs) a
|
||||
indexRange [] v = v
|
||||
|
|
|
@ -176,7 +176,7 @@ Traversable (Point n) where
|
|||
|
||||
export
|
||||
Show a => Show (Point n a) where
|
||||
showPrec d (MkPoint v) = showCon d "point" $ showArg $ toVect v
|
||||
showPrec d (MkPoint v) = showCon d "point" $ showArg $ elements v
|
||||
|
||||
export
|
||||
Cast a b => Cast (Point n a) (Point n b) where
|
||||
|
|
|
@ -25,13 +25,12 @@ isTranslation : Eq a => Num a => HMatrix' n a -> Bool
|
|||
isTranslation {n} mat with (viewShape mat)
|
||||
_ | Shape [S n,S n] = isHMatrix mat && getMatrix mat == identity
|
||||
|
||||
||| Try to construct a translation from a homogeneous matrix.
|
||||
export
|
||||
fromHMatrix : Eq a => Num a => HMatrix' n a -> Maybe (Translation n a)
|
||||
fromHMatrix mat = if isTranslation mat then Just (unsafeMkTrans mat) else Nothing
|
||||
|
||||
|
||||
||| Construct a translation given a vector.
|
||||
export
|
||||
translate : Num a => Vector n a -> Translation n a
|
||||
translate v = unsafeMkTrans (translationH v)
|
||||
|
||||
||| Try to construct a translation from a homogeneous matrix.
|
||||
export
|
||||
fromHMatrix : Eq a => Num a => HMatrix' n a -> Maybe (Translation n a)
|
||||
fromHMatrix mat = if isTranslation mat then Just (unsafeMkTrans mat) else Nothing
|
||||
|
|
|
@ -186,8 +186,8 @@ perp a b = a.x * b.y - a.y * b.x
|
|||
||| Calculate the cross product of the two vectors.
|
||||
export
|
||||
cross : Neg a => Vector 3 a -> Vector 3 a -> Vector 3 a
|
||||
cross v1 v2 = let [a, b, c] = toVect v1
|
||||
[x, y, z] = toVect v2
|
||||
cross v1 v2 = let [a, b, c] = elements v1
|
||||
[x, y, z] = elements v2
|
||||
in vector [b*z - c*y, c*x - a*z, a*y - b*x]
|
||||
|
||||
||| Calculate the triple product of the three vectors.
|
||||
|
|
Loading…
Reference in a new issue