Compare commits

...

10 commits

15 changed files with 336 additions and 107 deletions

View file

@ -28,7 +28,7 @@ combine the most useful features of each library.
## Documentation ## Documentation
Most of the exported utility functions have docstrings. There is also Most of the exported utility functions have docstrings. There is also
a (currently unfinished) guide on how to use NumIdr [here](docs/Intro.md). a (currently unfinished) guide on how to use NumIdr [here](docs/Contents.md).
## Usage ## Usage

19
docs/Contents.md Normal file
View file

@ -0,0 +1,19 @@
# 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!)

View file

@ -49,9 +49,6 @@ 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 > (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 > axes. Conventionally, NumIdr reserves "dimension" for the second meaning, and uses
> "rank" for the first meaning. > "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 ## Types of Arrays
@ -79,8 +76,6 @@ 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. 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 ### Matrices
As mentioned before, a matrix is a rank-2 array: As mentioned before, a matrix is a rank-2 array:
@ -121,36 +116,18 @@ These are useful for clarity when working with both homogeneous and non-homogene
### Transforms ### Transforms
A transform is a wrapper type for a matrix with certain properties that can be used to transform points in space. A transform is a wrapper type for a square matrix with certain properties that can be used to transform points in space.
```idris ```idris
Transform : (ty : TransType) -> (n : Nat) -> (a : Type) -> Type Transform : (ty : TransType) -> (n : Nat) -> (a : Type) -> Type
``` ```
The `TransType` parameter dictates what kind of transform it is. These eight options are currently available: 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.
**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 ### 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 various operations, such as those that permute the axes of an array. 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.
[Contents](Intro.md) | [Next](Operations.md) Permutations can be composed using `(*.)`, and a permutation can be converted into a matrix using `permuteM`.
[Contents](Contents.md) | [Next](Operations.md)

View file

@ -1,16 +0,0 @@
# 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

View file

@ -1,5 +1,10 @@
# Basic Operations on Arrays # 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 ## Constructing Arrays
The most important array constructor is `array`, which returns an array of the specified values: The most important array constructor is `array`, which returns an array of the specified values:
@ -8,7 +13,7 @@ The most important array constructor is `array`, which returns an array of the s
array [[1, 2, 3], [4, 5, 6]] 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` whereever 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` wherever possible, as they provide more information to the type-checker.
There are also a few other, more basic constructors for convenience. There are also a few other, more basic constructors for convenience.
@ -27,14 +32,13 @@ 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. 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.
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: > [!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) > ```idris
_ | Shape s = ... > 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 ## Indexing Arrays
@ -62,7 +66,7 @@ Not all combinations of these categories are defined by the library. Here are th
| **Update** | `indexUpdate` | `indexUpdateRange` | `indexUpdateNB` | | | | | **Update** | `indexUpdate` | `indexUpdateRange` | `indexUpdateNB` | | | |
| **Set** | `indexSet` | `indexSetRange` | `indexSetNB` | | | | | **Set** | `indexSet` | `indexSetRange` | `indexSetNB` | | | |
The accessor functions have operator forms for convenience, also specified within the table. The accessor functions have operator forms for convenience.
### Specifying Coordinates ### Specifying Coordinates
@ -77,20 +81,23 @@ 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: 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` to `y` - `Bounds x y` - Every index from `x` (inclusive) to `y` (exclusive)
- `StartBound x` - Every index from `x` to the end of the axis - `StartBound x` - Every index from `x` to the end of the axis
- `EndBound y` - Every index from the start of the axis to `y` - `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 - `All` - Every index in the axis
- `Indices xs` - A list of indices in the axis - `Indices xs` - A list of indices in the axis
- `Filter p` - A predicate specifying which indices to access or modify - `Filter p` - A predicate specifying which indices to access or modify
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: 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:
```idris ```idris
matrix [[2, 0], [1, 2]] !!.. [All, One 0] matrix [[2, 0], [1, 2]] !!.. [All, One 0]
== vector [2, 1] == vector [2, 1]
``` ```
`One'` does not decrease the rank of the result in this way.
## Standard Interface Methods ## 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. The array type implements a number of standard interfaces. Most of these implementations are unremarkable, but a few have caveats that are worth noting.
@ -106,8 +113,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). 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).
> [!NOTE] > [!WARNING]
> 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. > 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.
### `Foldable` and `Traversable` ### `Foldable` and `Traversable`
@ -117,7 +124,7 @@ When folding or traversing the elements of an array, these elements are ordered
### Concatenation and Stacking ### Concatenation and Stacking
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. 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.
```idris ```idris
-- 0 is the first axis i.e. the row axis -- 0 is the first axis i.e. the row axis
@ -128,7 +135,7 @@ concat 0 (matrix [[1, 2], [3, 4]]) (matrix [[5, 6], [7, 8]])
[7, 8]] [7, 8]]
``` ```
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. 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.
```idris ```idris
stack 0 [vector [1, 2], vector [3, 4]] stack 0 [vector [1, 2], vector [3, 4]]
@ -149,7 +156,7 @@ reshape [3, 2] (vector [1, 2, 3, 4, 5, 6])
[5, 6]] [5, 6]]
``` ```
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. 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.
```idris ```idris
resize [2, 4] 10 (matrix [[1, 2], resize [2, 4] 10 (matrix [[1, 2],
@ -163,8 +170,10 @@ Instead of the `resize` function, one can also use the `resizeLTE` function, whi
### Transpose ### Transpose
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)`. 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)`.
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. 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.
[Previous](DataTypes.md) | [Contents](Intro.md) | [Next](VectorsMatrices.md) 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)

View file

@ -1,4 +1,108 @@
# Transforms # 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.
[Previous](VectorsMatrices.md) | [Contents](Intro.md) 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)

View file

@ -1,4 +1,114 @@
# Working with Vectors and Matrices # 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).
[Previous](Operations.md) | [Contents](Intro.md) | [Next](Transforms.md) ## 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)

View file

@ -1,5 +1,5 @@
package numidr package numidr
version = 2.0.0 version = 0.3.0
brief = "Linear algebra and data science library" brief = "Linear algebra and data science library"
authors = "Kiana Sheibani" authors = "Kiana Sheibani"

View file

@ -26,6 +26,9 @@ export
rangeLenZ : (x : Nat) -> length (range 0 x) = x rangeLenZ : (x : Nat) -> length (range 0 x) = x
rangeLenZ x = rangeLen 0 x `trans` minusZeroRight 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 -- Array coordinate types
@ -63,6 +66,7 @@ namespace Strict
public export public export
data CRange : Nat -> Type where data CRange : Nat -> Type where
One : Fin n -> CRange n One : Fin n -> CRange n
One' : Fin n -> CRange n
All : CRange n All : CRange n
StartBound : Fin (S n) -> CRange n StartBound : Fin (S n) -> CRange n
EndBound : Fin (S n) -> CRange n EndBound : Fin (S n) -> CRange n
@ -79,6 +83,8 @@ namespace Strict
namespace NB namespace NB
public export public export
data CRangeNB : Type where data CRangeNB : Type where
One : Nat -> CRangeNB
One' : Nat -> CRangeNB
All : CRangeNB All : CRangeNB
StartBound : Nat -> CRangeNB StartBound : Nat -> CRangeNB
EndBound : Nat -> CRangeNB EndBound : Nat -> CRangeNB
@ -128,6 +134,7 @@ namespace Strict
public export public export
cRangeToList : {n : Nat} -> CRange n -> Either Nat (List Nat) cRangeToList : {n : Nat} -> CRange n -> Either Nat (List Nat)
cRangeToList (One x) = Left (cast x) cRangeToList (One x) = Left (cast x)
cRangeToList (One' x) = Right [cast x]
cRangeToList All = Right $ range 0 n cRangeToList All = Right $ range 0 n
cRangeToList (StartBound x) = Right $ range (cast x) n cRangeToList (StartBound x) = Right $ range (cast x) n
cRangeToList (EndBound x) = Right $ range 0 (cast x) cRangeToList (EndBound x) = Right $ range 0 (cast x)
@ -148,8 +155,8 @@ namespace Strict
newShape : {s : _} -> (rs : CoordsRange s) -> Vect (newRank rs) Nat newShape : {s : _} -> (rs : CoordsRange s) -> Vect (newRank rs) Nat
newShape [] = [] newShape [] = []
newShape (r :: rs) with (cRangeToList r) newShape (r :: rs) with (cRangeToList r)
newShape (r :: rs) | Left _ = newShape rs _ | Left _ = newShape rs
newShape (r :: rs) | Right xs = length xs :: newShape rs _ | Right xs = length xs :: newShape rs
getNewPos : {s : _} -> (rs : CoordsRange {rk} s) -> Vect rk Nat -> Vect (newRank rs) Nat getNewPos : {s : _} -> (rs : CoordsRange {rk} s) -> Vect rk Nat -> Vect (newRank rs) Nat
@ -175,6 +182,14 @@ namespace NB
validateCRange (d :: s) (r :: rs) = [| validate' d r :: validateCRange s rs |] validateCRange (d :: s) (r :: rs) = [| validate' d r :: validateCRange s rs |]
where where
validate' : (n : Nat) -> CRangeNB -> Maybe (CRange n) 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 All = Just All
validate' n (StartBound x) = validate' n (StartBound x) =
case isLTE x n of case isLTE x n of
@ -197,29 +212,44 @@ namespace NB
export %unsafe export %unsafe
assertCRange : (s : Vect rk Nat) -> Vect rk CRangeNB -> CoordsRange s assertCRange : (s : Vect rk Nat) -> Vect rk CRangeNB -> CoordsRange s
assertCRange [] [] = [] assertCRange [] [] = []
assertCRange (d :: s) (r :: rs) = assert' d r :: assertCRange s rs assertCRange (d :: s) (r :: rs) = assert' r :: assertCRange s rs
where where
assert' : (n : Nat) -> CRangeNB -> CRange n assert' : forall n. CRangeNB -> CRange n
assert' n All = All assert' (One i) = One (assertFin i)
assert' n (StartBound x) = StartBound (believe_me x) assert' (One' i) = One' (assertFin i)
assert' n (EndBound x) = EndBound (believe_me x) assert' All = All
assert' n (Bounds x y) = Bounds (believe_me x) (believe_me y) assert' (StartBound x) = StartBound (assertFin x)
assert' n (Indices xs) = Indices (believe_me <$> xs) assert' (EndBound x) = EndBound (assertFin x)
assert' n (Filter f) = Filter (f . finToNat) assert' (Bounds x y) = Bounds (assertFin x) (assertFin y)
assert' (Indices xs) = Indices (assertFin <$> xs)
assert' (Filter f) = Filter (f . finToNat)
public export public export
cRangeNBToList : Nat -> CRangeNB -> List Nat cRangeNBToList : Nat -> CRangeNB -> Either Nat (List Nat)
cRangeNBToList s All = range 0 s cRangeNBToList s (One i) = Left i
cRangeNBToList s (StartBound x) = range x s cRangeNBToList s (One' i) = Right [i]
cRangeNBToList s (EndBound x) = range 0 x cRangeNBToList s All = Right $ range 0 s
cRangeNBToList s (Bounds x y) = range x y cRangeNBToList s (StartBound x) = Right $ range x s
cRangeNBToList s (Indices xs) = nub xs cRangeNBToList s (EndBound x) = Right $ range 0 x
cRangeNBToList s (Filter p) = filter p $ range 0 s 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)
||| Calculate the new shape given by a coordinate range. ||| Calculate the new shape given by a coordinate range.
public export public export
newShape : Vect rk Nat -> Vect rk CRangeNB -> Vect rk Nat newShape : (s : Vect rk Nat) -> (is : Vect rk CRangeNB) -> Vect (newRank s is) Nat
newShape = zipWith (length .: cRangeNBToList) newShape [] [] = []
newShape (d :: s) (r :: rs) with (cRangeNBToList d r)
_ | Left _ = newShape s rs
_ | Right xs = length xs :: newShape s rs
export export
getAllCoords' : Vect rk Nat -> List (Vect rk Nat) getAllCoords' : Vect rk Nat -> List (Vect rk Nat)

View file

@ -31,16 +31,16 @@ export
FieldCmp Double where FieldCmp Double where
abslt = (<) `on` abs abslt = (<) `on` abs
-- Alternative implementations of `Eq` and `FieldCmp` that compare floating -- Alternative implementations of `Eq` and `FieldCmp` that compare approximately,
-- point numbers approximately, useful when working with transforms -- useful when working with flating point numbers
namespace Eq namespace Eq
export export
WithEpsilon : Double -> Eq Double WithEpsilon : (Neg a, Abs a, Ord a) => a -> Eq a
WithEpsilon ep = MkEq (\x,y => abs (x - y) < ep) (\x,y => abs (x - y) >= ep) WithEpsilon ep = MkEq (\x,y => abs (x - y) < ep) (\x,y => abs (x - y) >= ep)
namespace FieldCmp namespace FieldCmp
export export
WithEpsilon : Double -> FieldCmp Double WithEpsilon : (Neg a, Fractional a, Abs a, Ord a) => a -> FieldCmp a
WithEpsilon ep = MkFieldCmp @{WithEpsilon ep} ((<) `on` abs) WithEpsilon ep = MkFieldCmp @{WithEpsilon ep} ((<) `on` abs)
-------------------------------------------------------------------------------- --------------------------------------------------------------------------------

View file

@ -568,9 +568,8 @@ solveLowerTri' {n} mat b with (viewShape b)
construct [] = [] construct [] = []
construct {i=S i} (b :: bs) = construct {i=S i} (b :: bs) =
let xs = construct 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=()}) $ in (b - sum (zipWith (*) xs (reverse $ toVect $ replace {p = flip Array a} (believe_me $ Refl {x=()}) $
mat !!.. [One i', EndBound (weaken i')]))) / mat!#[i,i] :: xs mat !#.. [One i, EndBound i]))) / mat!#[i,i] :: xs
solveUpperTri' : Field a => Matrix' n a -> Vector n a -> Vector n a solveUpperTri' : Field a => Matrix' n a -> Vector n a -> Vector n a
@ -581,9 +580,8 @@ solveUpperTri' {n} mat b with (viewShape b)
construct _ [] = [] construct _ [] = []
construct i (b :: bs) = construct i (b :: bs) =
let xs = construct (S i) 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=()}) $ in (b - sum (zipWith (*) xs (toVect $ replace {p = flip Array a} (believe_me $ Refl {x=()}) $
mat !!.. [One i', StartBound (FS i')]))) / mat!#[i,i] :: xs mat !#.. [One i, StartBound (S i)]))) / mat!#[i,i] :: xs
||| Solve a linear equation, assuming the matrix is lower triangular. ||| Solve a linear equation, assuming the matrix is lower triangular.
@ -628,6 +626,7 @@ solve mat = solveWithLUP mat (decompLUP mat)
-------------------------------------------------------------------------------- --------------------------------------------------------------------------------
||| Determine whether a matrix has an inverse.
export export
invertible : FieldCmp a => Matrix' n a -> Bool invertible : FieldCmp a => Matrix' n a -> Bool
invertible {n} mat with (viewShape mat) invertible {n} mat with (viewShape mat)

View file

@ -23,10 +23,6 @@ update : Coords s -> (a -> a) -> Vects s a -> Vects s a
update [] f v = f v update [] f v = f v
update (i :: is) f v = updateAt i (update is 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 export
indexRange : {s : _} -> (rs : CoordsRange s) -> Vects s a -> Vects (newShape rs) a indexRange : {s : _} -> (rs : CoordsRange s) -> Vects s a -> Vects (newShape rs) a
indexRange [] v = v indexRange [] v = v

View file

@ -176,7 +176,7 @@ Traversable (Point n) where
export export
Show a => Show (Point n a) where Show a => Show (Point n a) where
showPrec d (MkPoint v) = showCon d "point" $ showArg $ elements v showPrec d (MkPoint v) = showCon d "point" $ showArg $ toVect v
export export
Cast a b => Cast (Point n a) (Point n b) where Cast a b => Cast (Point n a) (Point n b) where

View file

@ -25,12 +25,13 @@ isTranslation : Eq a => Num a => HMatrix' n a -> Bool
isTranslation {n} mat with (viewShape mat) isTranslation {n} mat with (viewShape mat)
_ | Shape [S n,S n] = isHMatrix mat && getMatrix mat == identity _ | Shape [S n,S n] = isHMatrix mat && getMatrix mat == identity
||| 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. ||| Try to construct a translation from a homogeneous matrix.
export export
fromHMatrix : Eq a => Num a => HMatrix' n a -> Maybe (Translation n a) fromHMatrix : Eq a => Num a => HMatrix' n a -> Maybe (Translation n a)
fromHMatrix mat = if isTranslation mat then Just (unsafeMkTrans mat) else Nothing 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)

View file

@ -186,8 +186,8 @@ perp a b = a.x * b.y - a.y * b.x
||| Calculate the cross product of the two vectors. ||| Calculate the cross product of the two vectors.
export export
cross : Neg a => Vector 3 a -> Vector 3 a -> Vector 3 a cross : Neg a => Vector 3 a -> Vector 3 a -> Vector 3 a
cross v1 v2 = let [a, b, c] = elements v1 cross v1 v2 = let [a, b, c] = toVect v1
[x, y, z] = elements v2 [x, y, z] = toVect v2
in vector [b*z - c*y, c*x - a*z, a*y - b*x] in vector [b*z - c*y, c*x - a*z, a*y - b*x]
||| Calculate the triple product of the three vectors. ||| Calculate the triple product of the three vectors.