Various guide edits

Is anyone actually using this?

README.md

@@ -27,14 +27,22 @@ combine the most useful features of each library.
 
 ## Documentation
 
-Most of the exported utility functions have docstrings. There are also plans for
-an in-depth tutorial on NumIdr's features, though that is not available yet.
+Most of the exported utility functions have docstrings. There is also
+a (currently unfinished) guide on how to use NumIdr [here](docs/Contents.md).
 
 ## Usage
 
-Basic package install:
+To install using idris2 directly:
 
-``` shell
+``` sh
+git clone https://github.com/kiana-S/numidr
+cd numidr
 idris2 --install numidr.ipkg
+
 ```
 
+Or you can install using pack:
+
+``` sh
+pack install numidr
+```

docs/Contents.md

@@ -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!)

docs/DataTypes.md

@@ -0,0 +1,133 @@
+# Fundamental Data Types
+
+NumIdr exports a number of different datatypes. The most important type and the cornerstone of the library is the _array_, but there are other useful types as well.
+
+## Arrays
+
+### What is an Array?
+
+In most programming languages, the word "array" is used to mean a one-dimensional list of values that is contiguous in memory. A typical array of integers may be written in list form like this:
+
+```
+[1, 4, 10, 2, -5, 18]
+```
+
+In this kind of array, elements are indexed by a single integer, starting at zero and increasing from left to right.
+
+NumIdr, however, uses the word a bit more generally: a NumIdr array is a multi-dimensional structure that can be indexed by any number of integers. NumIdr arrays are written as nested lists:
+
+```
+[[4, -9, -2],
+ [5, -6, 1]]
+```
+
+Unlike in other languages, however, this is not a nested structure. The above is a single array, and it is always manipulated as one object.
+
+### Properties of Arrays
+
+The `Array` datatype has the following parameters:
+
+```idris
+Array : (s : Vect rk Nat) -> (a : Type) -> Type
+```
+
+The first parameter is the _shape_, a list of numbers (the _dimensions_) where each dimension is the length of a particular _axis_ of the array. The second parameter is the _element type_, the type of the values inside the array.
+
+Let's return to the array example from earlier:
+
+```
+[[4, -9, -2],
+ [5, -6, 1]]
+```
+
+This is a rank-2 array, meaning that it has two axes. Rank-2 arrays are typically called matrices. To determine the dimensions of the array, we count the size of each nested list from the outside in, which in the case of matrices means the row axis comes before the column axis. This matrix has 2 rows and 3 columns, making its shape `[2, 3]`. Thus, a possible type for this array could be `Array [2, 3] Int`.
+
+When determining the index of a value inside the array, the order of the indices is the same as the order of the dimensions, and each index number counts from zero. For example, the index `[1, 0]` indicates the second row and first column, which contains `5`.
+
+> [!NOTE]
+> The word "dimensions" is often ambiguously used to either refer to the rank of an array
+> (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.
+
+## Types of Arrays
+
+Arrays are loosely divided into multiple subtypes mostly based on their rank. Each array subtype has an alias for convenience.
+
+### Scalars
+
+A scalar is a rank-0 array, meaning that it is indexed by 0 integers. Its alias is `Scalar`:
+
+```idris
+Scalar : (a : Type) -> Type
+Scalar a = Array [] a
+```
+
+A scalar has exactly one index, the empty list `[]`. This means that it is exactly the same as a single value and as such is largely pointless, but NumIdr still provides an alias for it just in case you need it.
+
+### Vectors
+
+A vector is a rank-1 array:
+
+```idris
+Vector : (n : Nat) -> (a : Type) -> Type
+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.
+
+### Matrices
+
+As mentioned before, a matrix is a rank-2 array:
+
+```idris
+Matrix : (m, n : Nat) -> (a : Type) -> Type
+Matrix m n a = Array [m, n] a
+```
+
+There is also an alias `Matrix'` for square matrices.
+
+```idris
+Matrix' : (n : Nat) -> (a : Type) -> Type
+Matrix' n a = Array [n, n] a
+```
+
+As a linear algebra library, the majority of the operations in NumIdr revolve around matrices.
+
+#### Homogeneous Matrices
+
+NumIdr also provides aliases for homogeneous matrices:
+
+```idris
+HMatrix : (m, n : Nat) -> (a : Type) -> Type
+HMatrix m n a = Array [S m, S n] a
+
+HMatrix' : (n : Nat) -> (a : Type) -> Type
+HMatrix' n a = Array [S n, S n] a
+
+-- To use with homogeneous matrices
+HVector : (n : Nat) -> (a : Type) -> Type
+HVector n a = Array [S n] a
+```
+
+These are useful for clarity when working with both homogeneous and non-homogeneous matrices.
+
+## Other Datatypes
+
+### Transforms
+
+A transform is a wrapper type for a square 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.
+
+### 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.
+
+Permutations can be composed using `(*.)`, and a permutation can be converted into a matrix using `permuteM`.
+
+[Contents](Contents.md) | [Next](Operations.md)

docs/Operations.md

@@ -0,0 +1,179 @@
+# 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:
+
+```idris
+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.
+
+There are also a few other, more basic constructors for convenience.
+
+```idris
+-- A 2x2x3 array filled with zeros
+repeat 0 [2, 2, 3]
+
+-- Same as previous
+zeros [2, 2, 3]
+
+-- A 2x2x3 array filled with ones
+ones [2, 2, 3]
+```
+
+## Accessing Array Properties
+
+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 = ...
+> ```
+
+## Indexing Arrays
+
+NumIdr provides multiple different indexing functions for different purposes. These functions can be grouped based on these categories:
+
+**Operation**
+- **Access** - Accesses and returns elements from the array.
+- **Update** - Returns a new array with the specified element or range updated by a function. These are indicated by an `Update` suffix.
+- **Set** - Returns a new array with the specified element or range set to new values. These are indicated by a `Set` suffix.
+
+**Range**
+- **Default** - Operates on a single array element.
+- **Range** - Operates on multiple elements at once. These are indicated by a `Range` suffix.
+
+**Safety**
+- **Safe** - Guarantees through its type that the index is within range by requiring each index to be a `Fin` value.
+- **Non-Bounded** - Does not guarantee through its type that the index is within range, and returns `Nothing` if the provided index is out of bounds. These are indicated by an `NB` suffix.
+- **Unsafe** - Does not perform any bounds checks at all. These are indicated by an `Unsafe` suffix. **Only use these if you really know what you are doing!**
+
+Not all combinations of these categories are defined by the library. Here are the currently provided indexing functions:
+
+|            | Safe            | Ranged Safe            | Non-Bounded       | Ranged Non-Bounded       | Unsafe                | Ranged Unsafe                |
+|------------|-----------------|------------------------|-------------------|--------------------------|-----------------------|------------------------------|
+| **Access** | `index`, `(!!)` | `indexRange`, `(!!..)` | `indexNB`, `(!?)` | `indexRangeNB`, `(!?..)` | `indexUnsafe`, `(!#)` | `indexRangeUnsafe`, `(!#..)` |
+| **Update** | `indexUpdate`   | `indexUpdateRange`     | `indexUpdateNB`   |                          |                       |                              |
+| **Set**    | `indexSet`      | `indexSetRange`        | `indexSetNB`      |                          |                       |                              |
+
+The accessor functions have operator forms for convenience.
+
+### Specifying Coordinates
+
+When indexing a single element, a coordinate is specified as a list of numbers, where each number is either a `Fin` value for safe indexing or a `Nat` value for non-bounded and unsafe indexing.
+
+```idris
+index [1, 0] arr
+
+-- Equivalent operator form
+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)
+- `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:
+
+```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.
+
+### Numeric Operations
+
+Arrays implement the numeric interfaces (`Num`, `Neg`, `Fractional`), as well as `Semigroup` and `Monoid`, if their element type supports those operations. These functions are computed elementwise.
+
+```idris
+matrix [[1, 1], [2, 5]] + matrix [[2, 3], [-1, 3]]
+ == matrix [[3, 4], [1, 8]]
+```
+
+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.
+
+### `Foldable` and `Traversable`
+
+When folding or traversing the elements of an array, these elements are ordered in row-major or "C-style" order, which corresponds to the order in which elements are written and displayed. This behavior should not be depended on, however, as it can change based on the internal [array representation](Representations.md); use the `elements` function if you specifically want row-major order.
+
+## Other Common Operations
+
+### 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.
+
+```idris
+-- 0 is the first axis i.e. the row axis
+concat 0 (matrix [[1, 2], [3, 4]]) (matrix [[5, 6], [7, 8]])
+  == matrix [[1, 2],
+             [3, 4],
+             [5, 6],
+             [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.
+
+```idris
+stack 0 [vector [1, 2], vector [3, 4]]
+  == matrix [[1, 2],
+             [3, 4]]
+```
+
+There are also specialized functions for operating on vectors and matrices: `(++)` for concatenating vectors; `vconcat` and `hconcat` for concatenating two matrices vertically (by row) and horizontally (by column) respectively, and `vstack` and `hstack` for stacking row vectors and column vectors respectively.
+
+### Reshaping and Resizing
+
+An array can be "reshaped" into any other shape, so long as the total number of elements is the same. This reshaping is done by arranging the elements into a linear order before inserting them into the new array. As with folding, the default order is row-major.
+
+```idris
+reshape [3, 2] (vector [1, 2, 3, 4, 5, 6])
+  == matrix [[1, 2],
+             [3, 4],
+             [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.
+
+```idris
+resize [2, 4] 10 (matrix [[1, 2],
+                          [3, 4],
+                          [5, 6]])
+  == matrix [[1, 2, 10, 10],
+             [3, 4, 10, 10]]
+```
+
+Instead of the `resize` function, one can also use the `resizeLTE` function, which does not require a default element, but only works if the new array would be provably smaller than the original one.
+
+### 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)`.
+
+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)

docs/Representations.md

@@ -0,0 +1,2 @@
+# Array Representations
+

docs/Transforms.md

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

docs/VectorsMatrices.md

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

numidr.ipkg

@@ -1,22 +1,17 @@
 package numidr
-version = 1.0.0
+version = 0.3.0
+brief = "Linear algebra and data science library"
 
 authors = "Kiana Sheibani"
 license = "MIT"
 
-langversion >= 0.5.1
+langversion >= 0.6.0
 
 sourcedir = "src"
 readme = "README.md"
 
-modules = Data.NP,
-          Data.Permutation,
+modules = Data.Permutation,
           Data.NumIdr,
-          Data.NumIdr.PrimArray,
-          Data.NumIdr.PrimArray.Bytes,
-          Data.NumIdr.PrimArray.Boxed,
-          Data.NumIdr.PrimArray.Linked,
-          Data.NumIdr.PrimArray.Delayed,
           Data.NumIdr.Array,
           Data.NumIdr.Array.Array,
           Data.NumIdr.Array.Coords,

sirdi.json

@@ -1,39 +0,0 @@
-[
-  {
-    "name": "numidr",
-    "version": "1.0.0",
-    "deps": [],
-    "modules": [
-      "Data.NP",
-      "Data.Permutation",
-      "Data.NumIdr",
-      "Data.NumIdr.Array",
-      "Data.NumIdr.Array.Array",
-      "Data.NumIdr.Array.Coords",
-      "Data.NumIdr.Array.Order",
-      "Data.NumIdr.Homogeneous",
-      "Data.NumIdr.Matrix",
-      "Data.NumIdr.Interfaces",
-      "Data.NumIdr.PrimArray",
-      "Data.NumIdr.Scalar",
-      "Data.NumIdr.Vector",
-      "Data.NumIdr.LArray",
-      "Data.NumIdr.Transform",
-      "Data.NumIdr.Transform.Affine",
-      "Data.NumIdr.Transform.Isometry",
-      "Data.NumIdr.Transform.Linear",
-      "Data.NumIdr.Transform.Orthonormal",
-      "Data.NumIdr.Transform.Point",
-      "Data.NumIdr.Transform.Rigid",
-      "Data.NumIdr.Transform.Rotation",
-      "Data.NumIdr.Transform.Transform",
-      "Data.NumIdr.Transform.Translation",
-      "Data.NumIdr.Transform.Trivial"
-    ],
-    "passthru": {
-      "authors": "Kiana Sheibani",
-      "license": "MIT",
-      "sourceloc": "https://github.com/kiana-S/numidr"
-    }
-  }
-]

src/Data/NP.idr

@@ -1,43 +0,0 @@
-module Data.NP
-
-import Data.Vect
-
-%default total
-
-
-||| A type constructor for heterogenous lists.
-public export
-data NP : (f : k -> Type) -> (ts : Vect n k) -> Type where
-  Nil  : NP {k} f []
-  (::) : f t -> NP {k} f ts -> NP {k} f (t :: ts)
-
-
-public export
-(all : NP (Eq . f) ts) => Eq (NP f ts) where
-  (==) {all=[]} [] [] = True
-  (==) {all=_::_} (x :: xs) (y :: ys) = (x == y) && (xs == ys)
-
-public export
-(all : NP (Show . f) ts) => Show (NP f ts) where
-  show xs = "[" ++ elems xs ++ "]"
-    where
-      elems : {0 f,ts : _} -> (all : NP (Show . f) ts) => NP f ts -> String
-      elems {all=[]} [] = ""
-      elems {all=[_]} [x] = show x
-      elems {all=_::_} (x :: xs) = show x ++ ", " ++ elems xs
-
--- use this when idris can't figure out the constraint above
-public export
-eqNP : forall f. (forall t. Eq (f t)) => (x, y : NP f ts) -> Bool
-eqNP [] [] = True
-eqNP (x :: xs) (y :: ys) = x == y && eqNP xs ys
-
-
-public export
-head : NP f (t :: ts) -> f t
-head (x :: _) = x
-
-public export
-index : (i : Fin n) -> NP {n} f ts -> f (index i ts)
-index FZ (x :: xs) = x
-index (FS i) (x :: xs) = index i xs

src/Data/NumIdr.idr

@@ -1,8 +1,8 @@
 module Data.NumIdr
 
-import public Data.NP
 import public Data.Permutation
 
+import Data.NumIdr.PrimArray
 import public Data.NumIdr.Interfaces
 import public Data.NumIdr.Array
 import public Data.NumIdr.Scalar

src/Data/NumIdr/Array.idr

@@ -1,6 +1,5 @@
 module Data.NumIdr.Array
 
-import public Data.NP
 import public Data.NumIdr.Array.Array
 import public Data.NumIdr.Array.Coords
 import public Data.NumIdr.Array.Rep

src/Data/NumIdr/Array/Array.idr

@@ -3,7 +3,6 @@ module Data.NumIdr.Array.Array
 import Data.List
 import Data.Vect
 import Data.Zippable
-import Data.NP
 import Data.Permutation
 import Data.NumIdr.Interfaces
 import Data.NumIdr.PrimArray
@@ -13,6 +12,8 @@ import Data.NumIdr.Array.Coords
 %default total
 
 
+||| The type of an array.
+|||
 ||| Arrays are the central data structure of NumIdr. They are an `n`-dimensional
 ||| grid of values, where `n` is a value known as the *rank* of the array. Arrays
 ||| of rank 0 are single values, arrays of rank 1 are vectors, and arrays of rank
@@ -23,22 +24,18 @@ import Data.NumIdr.Array.Coords
 ||| array's total size.
 |||
 ||| Arrays are indexed by row first, as in the standard mathematical notation for
-||| matrices. This is independent of the actual order in which they are stored; the
-||| default order is row-major, but this is configurable.
+||| matrices.
 |||
 ||| @ rk The rank of the array
 ||| @ s The shape of the array
 ||| @ a The type of the array's elements
 export
 data Array : (s : Vect rk Nat) -> (a : Type) -> Type where
-  ||| Internally, arrays are stored using Idris's primitive array type. This is
-  ||| stored along with the array's shape, and a vector of values called the
-  ||| *strides*, which determine how indexes into the internal array should be
-  ||| performed. This is how the order of the array is configurable.
+  ||| Internally, arrays are stored via one of a handful of representations.
   |||
-  ||| @ ord The order of the elements of the array
-  ||| @ sts The strides of the array
   ||| @ s   The shape of the array
+  ||| @ rep The internal representation of the array
+  ||| @ rc  A witness that the element type satisfies the representation constraint
   MkArray : (rep : Rep) -> (rc : RepConstraint rep a) => (s : Vect rk Nat) ->
                PrimArray rep s a @{rc} -> Array s a
 
@@ -56,24 +53,32 @@ unsafeMkArray = MkArray
 --------------------------------------------------------------------------------
 
 
-||| The shape of the array
+||| The shape of the array.
 export
 shape : Array {rk} s a -> Vect rk Nat
 shape (MkArray _ s _) = s
 
-||| The rank of the array
+||| The size of the array, i.e. the total number of elements.
+export
+size : Array s a -> Nat
+size = product . shape
+
+||| The rank of the array.
 export
 rank : Array s a -> Nat
 rank = length . shape
 
+||| The internal representation of the array.
 export
 getRep : Array s a -> Rep
 getRep (MkArray rep _ _) = rep
 
+||| The representation constraint of the array.
 export
 getRepC : (arr : Array s a) -> RepConstraint (getRep arr) a
 getRepC (MkArray _ @{rc} _ _) = rc
 
+||| Extract the primitive backend array.
 export
 getPrim : (arr : Array s a) -> PrimArray (getRep arr) s a @{getRepC arr}
 getPrim (MkArray _ _ pr) = pr
@@ -171,22 +176,12 @@ fromFunction : {default B rep : Rep} -> RepConstraint rep a =>
                   (s : Vect rk Nat) -> (Coords s -> a) -> Array s a
 fromFunction s f = MkArray rep s (PrimArray.fromFunction s f)
 
-||| Construct an array using a structure of nested vectors. The elements are arranged
-||| to the specified order before being written.
-|||
-||| @ s   The shape of the constructed array
-||| @ ord The order of the constructed array
-export
-array' : {default B rep : Rep} -> RepConstraint rep a =>
-          (s : Vect rk Nat) -> Vects s a -> Array s a
-array' s v = MkArray rep s (fromVects s v)
-
 ||| Construct an array using a structure of nested vectors.
 ||| To explicitly specify the shape and order of the array, use `array'`.
 export
 array : {default B rep : Rep} -> RepConstraint rep a =>
           {s : Vect rk Nat} -> Vects s a -> Array s a
-array {rep} = array' {rep} _
+array v = MkArray rep s (fromVects s v)
 
 
 --------------------------------------------------------------------------------
@@ -340,14 +335,26 @@ arr !#.. is = indexRangeUnsafe is arr
 -- Operations on arrays
 --------------------------------------------------------------------------------
 
+||| Map a function over an array.
+|||
+||| You should almost always use `map` instead; only use this function if you
+||| know what you are doing!
 export
 mapArray' : (a -> a) -> Array s a -> Array s a
 mapArray' f (MkArray rep _ arr) = MkArray rep _ (mapPrim f arr)
 
+||| Map a function over an array.
+|||
+||| You should almost always use `map` instead; only use this function if you
+||| know what you are doing!
 export
 mapArray : (a -> b) -> (arr : Array s a) -> RepConstraint (getRep arr) b => Array s b
 mapArray f (MkArray rep _ arr) @{rc} = MkArray rep @{rc} _ (mapPrim f arr)
 
+||| Combine two arrays of the same element type using a binary function.
+|||
+||| You should almost always use `zipWith` instead; only use this function if
+||| you know what you are doing!
 export
 zipWithArray' : (a -> a -> a) -> Array s a -> Array s a -> Array s a
 zipWithArray' {s} f a b with (viewShape a)
@@ -356,6 +363,10 @@ zipWithArray' {s} f a b with (viewShape a)
                           $ PrimArray.fromFunctionNB @{mergeRepConstraint (getRepC a) (getRepC b)} _
                               (\is => f (a !# is) (b !# is))
 
+||| Combine two arrays using a binary function.
+|||
+||| You should almost always use `zipWith` instead; only use this function if
+||| you know what you are doing!
 export
 zipWithArray : (a -> b -> c) -> (arr : Array s a) -> (arr' : Array s b) ->
                 RepConstraint (mergeRep (getRep arr) (getRep arr')) c => Array s c
@@ -377,6 +388,12 @@ export
 convertRep : (rep : Rep) -> RepConstraint rep a => Array s a -> Array s a
 convertRep rep (MkArray _ s arr) = MkArray rep s (convertRepPrim arr)
 
+||| Temporarily convert an array to a delayed representation to make modifying
+||| it more efficient, then convert it back to its original representation.
+export
+delayed : (Array s a -> Array s' a) -> Array s a -> Array s' a
+delayed f arr = convertRep (getRep arr) @{getRepC arr} $ f $ convertRep Delayed arr
+
 ||| Resize the array to a new shape, preserving the coordinates of the original
 ||| elements. New coordinates are filled with a default value.
 |||
@@ -394,7 +411,7 @@ resize s' def arr = fromFunction {rep=getRep arr} @{getRepC arr} s' (fromMaybe d
 export
 -- HACK: Come up with a solution that doesn't use `believe_me` or trip over some
 -- weird bug in the type-checker
-resizeLTE : (s' : Vect rk Nat) -> (0 ok : NP Prelude.id (zipWith LTE s' s)) =>
+resizeLTE : (s' : Vect rk Nat) -> (0 ok : All Prelude.id (zipWith LTE s' s)) =>
             Array {rk} s a -> Array s' a
 resizeLTE s' arr = resize s' (believe_me ()) arr
 
@@ -446,6 +463,7 @@ stack axis arrs = rewrite sym (lengthCorrect arrs) in
     getAxisInd FZ (i :: is) = (i, is)
     getAxisInd {s=_::_} (FS ax) (i :: is) = mapSnd (i::) (getAxisInd ax is)
 
+||| Join the axes of a nested array structure to form a single array.
 export
 joinAxes : {s' : _} -> Array s (Array s' a) -> Array (s ++ s') a
 joinAxes {s} arr with (viewShape arr)
@@ -460,6 +478,7 @@ joinAxes {s} arr with (viewShape arr)
     dropUpTo (x::xs) (y::ys) = dropUpTo xs ys
 
 
+||| Split an array into a nested array structure along the specified axes.
 export
 splitAxes : (rk : Nat) -> {0 rk' : Nat} -> {s : _} ->
             Array {rk=rk+rk'} s a -> Array (take {m=rk'} rk s) (Array (drop {m=rk'} rk s) a)
@@ -577,7 +596,8 @@ export
 
 
 -- Foldable and Traversable operate on the primitive array directly. This means
--- that their operation is dependent on the order of the array.
+-- that their operation is dependent on the internal representation of the
+-- array.
 
 export
 Foldable (Array s) where

src/Data/NumIdr/Array/Coords.idr

@@ -4,7 +4,8 @@ import Data.Either
 import Data.List
 import Data.List1
 import Data.Vect
-import Data.NP
+
+import public Data.Vect.Quantifiers
 
 %default total
 
@@ -25,6 +26,9 @@ 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
@@ -35,7 +39,16 @@ rangeLenZ x = rangeLen 0 x `trans` minusZeroRight x
 ||| values of `Fin dim`, where `dim` is the dimension of each axis.
 public export
 Coords : (s : Vect rk Nat) -> Type
-Coords = NP Fin
+Coords = All Fin
+
+
+-- Occasionally necessary for reasons of Idris not being great at
+-- resolving interface constraints
+public export
+[eqCoords] Eq (Coords s) where
+  [] == [] = True
+  (x :: xs) == (y :: ys) = x == y && xs == ys
+
 
 ||| Forget the shape of the array by converting each index to type `Nat`.
 export
@@ -53,6 +66,7 @@ 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
@@ -63,12 +77,14 @@ namespace Strict
 
   public export
   CoordsRange : (s : Vect rk Nat) -> Type
-  CoordsRange = NP CRange
+  CoordsRange = All CRange
 
 
 namespace NB
   public export
   data CRangeNB : Type where
+    One : Nat -> CRangeNB
+    One' : Nat -> CRangeNB
     All : CRangeNB
     StartBound : Nat -> CRangeNB
     EndBound : Nat -> CRangeNB
@@ -118,6 +134,7 @@ 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)
@@ -138,8 +155,8 @@ namespace Strict
   newShape : {s : _} -> (rs : CoordsRange s) -> Vect (newRank rs) Nat
   newShape [] = []
   newShape (r :: rs) with (cRangeToList r)
-    newShape (r :: rs) | Left _ = newShape rs
-    newShape (r :: rs) | Right xs = length xs :: newShape rs
+    _ | Left _ = newShape rs
+    _ | Right xs = length xs :: newShape rs
 
 
   getNewPos : {s : _} -> (rs : CoordsRange {rk} s) -> Vect rk Nat -> Vect (newRank rs) Nat
@@ -165,6 +182,14 @@ 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
@@ -187,29 +212,44 @@ namespace NB
   export %unsafe
   assertCRange : (s : Vect rk Nat) -> Vect rk CRangeNB -> CoordsRange s
   assertCRange [] [] = []
-  assertCRange (d :: s) (r :: rs) = assert' d r :: assertCRange s rs
+  assertCRange (d :: s) (r :: rs) = assert' r :: assertCRange s rs
     where
-      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)
+      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)
 
   public export
-  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
+  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)
 
   ||| Calculate the new shape given by a coordinate range.
   public export
-  newShape : Vect rk Nat -> Vect rk CRangeNB -> Vect rk Nat
-  newShape = zipWith (length .: cRangeNBToList)
+  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
 
 export
 getAllCoords' : Vect rk Nat -> List (Vect rk Nat)

src/Data/NumIdr/Array/Rep.idr

@@ -12,10 +12,10 @@ import Data.Buffer
 --------------------------------------------------------------------------------
 
 
-||| An order is an abstract representation of the way in which array
-||| elements are stored in memory. Orders are used to calculate strides,
-||| which provide a method of converting an array coordinate into a linear
-||| memory location.
+||| An order is an abstract representation of the way in which contiguous array
+||| elements are stored in memory. Orders are used to calculate strides, which
+||| provide a method of converting an array coordinate into a linear memory
+||| location.
 public export
 data Order : Type where
   ||| C-like order, or contiguous order. This order stores elements in a
@@ -49,24 +49,57 @@ calcStrides FOrder v@(_::_)  = scanl (*) 1 $ init v
 --------------------------------------------------------------------------------
 
 
+||| An internal representation of an array.
+|||
+||| Each array internally stores its values based on one of these
+||| representations.
 public export
 data Rep : Type where
-  Bytes : Order -> Rep
-  Boxed : Order -> Rep
+  ||| Byte-buffer array representation.
+  |||
+  ||| This representations stores elements by converting them into byte values
+  ||| storing them in a `Buffer`. Use of this representation is only valid if
+  ||| the element type implements `ByteRep`.
+  |||
+  ||| @ o The order to store array elements in
+  Bytes : (o : Order) -> Rep
+
+  ||| Boxed array representation.
+  |||
+  ||| This representation uses a boxed array type to store its elements.
+  |||
+  ||| @ o The order to store array elements in
+  Boxed : (o : Order) -> Rep
+
+  ||| Linked list array representation.
+  |||
+  ||| This representation uses Idris's standard linked list types to store its
+  ||| elements.
   Linked : Rep
+
+  ||| Delayed/Lazy array representation.
+  |||
+  ||| This representation delays the computation of the array's elements, which
+  ||| may be useful in writing efficient algorithms.
   Delayed : Rep
 
 %name Rep rep
 
 
+||| An alias for the representation `Boxed COrder`, the C-like boxed array
+||| representation.
+|||
+||| This representation is the default for all newly created arrays.
 public export
 B : Rep
 B = Boxed COrder
 
+||| An alias for the representation `Linked COrder`.
 public export
 L : Rep
 L = Linked
 
+||| An alias for the representation `Delayed`.
 public export
 D : Rep
 D = Delayed
@@ -80,6 +113,8 @@ Eq Rep where
   Delayed == Delayed = True
   _ == _ = False
 
+||| A predicate on representations for those that store their elements
+||| linearly.
 public export
 data LinearRep : Rep -> Type where
   BytesIsL : LinearRep (Bytes o)
@@ -140,11 +175,19 @@ forceRepNCCorrect Delayed = DelayedNC
 --------------------------------------------------------------------------------
 
 
+||| An interface for elements that can be converted into raw bytes.
+|||
+||| An implementation of this interface is required to use the `Bytes` array
+||| representation.
 public export
 interface ByteRep a where
+  ||| The number of bytes used to store each value.
   bytes : Nat
 
+  ||| Convert a value into a list of bytes.
   toBytes : a -> Vect bytes Bits8
+
+  ||| Convert a list of bytes into a value.
   fromBytes : Vect bytes Bits8 -> a
 
 export
@@ -294,6 +337,11 @@ export
     in  fromBytes bs1 :: fromBytes bs2
 
 
+||| The constraint that each array representation requires.
+|||
+||| Currently, only the `Bytes` representation has a constraint, requiring an
+||| implementation of `ByteRep`. All other representations can be used without
+||| constraint.
 public export
 RepConstraint : Rep -> Type -> Type
 RepConstraint (Bytes _) a = ByteRep a

src/Data/NumIdr/Interfaces.idr

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

src/Data/NumIdr/Matrix.idr

@@ -33,7 +33,7 @@ Matrix' n = Array [n,n]
 export
 matrix : {default B rep : Rep} -> RepConstraint rep a => {m, n : _} ->
           Vect m (Vect n a) -> Matrix m n a
-matrix x = array' {rep} [m,n] x
+matrix x = array {rep, s=[m,n]} x
 
 
 ||| Construct a matrix with a specific value along the diagonal.
@@ -168,7 +168,8 @@ diagonal {n} mat with (viewShape mat)
 export
 -- TODO: throw an actual proof in here to avoid the unsafety
 minor : Fin (S m) -> Fin (S n) -> Matrix (S m) (S n) a -> Matrix m n a
-minor i j mat = believe_me $ mat!!..[Filter (/=i), Filter (/=j)]
+minor i j mat = replace {p = flip Array a} (believe_me $ Refl {x = ()})
+  $ mat!!..[Filter (/=i), Filter (/=j)]
 
 
 filterInd : Num a => (Nat -> Nat -> Bool) -> Matrix m n a -> Matrix m n a
@@ -403,12 +404,12 @@ gaussStep i lu =
 export
 decompLU : Field a => (mat : Matrix m n a) -> Maybe (DecompLU mat)
 decompLU {m,n} mat with (viewShape mat)
-  _ | Shape [m,n] = map (MkLU . convertRep _ @{getRepC mat})
-    $ iterateN (minimum m n) (\i => (>>= gaussStepMaybe i)) (Just (convertRep Delayed mat))
+  _ | Shape [m,n] = map MkLU
+    $ iterateN (minimum m n) (\i => (>>= gaussStepMaybe i)) (Just mat)
   where
     gaussStepMaybe : Fin (minimum m n) -> Matrix m n a -> Maybe (Matrix m n a)
     gaussStepMaybe i mat = if mat!#[cast i,cast i] == 0 then Nothing
-                           else Just (gaussStep i mat)
+                           else Just $ gaussStep i mat
 
 
 ||| The LUP decomposition of a matrix.
@@ -519,7 +520,7 @@ decompLUP {m,n} mat with (viewShape mat)
     maxIndex x [] = x
     maxIndex _ [x] = x
     maxIndex x ((a,b)::(c,d)::xs) =
-      if abslt b d then assert_total $ maxIndex x ((c,d)::xs)
+      if abslt b d then maxIndex x ((c,d)::xs)
                    else assert_total $ maxIndex x ((a,b)::xs)
 
     gaussStepSwap : Fin (S $ minimum m n) -> DecompLUP mat -> DecompLUP mat
@@ -567,9 +568,8 @@ 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 $ believe_me $
-                  mat !!.. [One i', EndBound (weaken i')]))) / mat!#[i,i] :: xs
+      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
 
 
 solveUpperTri' : Field a => Matrix' n a -> Vector n a -> Vector n a
@@ -580,9 +580,8 @@ 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 $ believe_me $
-                  mat !!.. [One i', StartBound (FS i')]))) / mat!#[i,i] :: xs
+      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
 
 
 ||| Solve a linear equation, assuming the matrix is lower triangular.
@@ -627,6 +626,7 @@ 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)

src/Data/NumIdr/PrimArray.idr

@@ -3,7 +3,6 @@ module Data.NumIdr.PrimArray
 import Data.Buffer
 import Data.Vect
 import Data.Fin
-import Data.NP
 import Data.NumIdr.Array.Rep
 import Data.NumIdr.Array.Coords
 import public Data.NumIdr.PrimArray.Bytes
@@ -98,7 +97,7 @@ indexUpdate {rep = Linked} is f arr = update is f arr
 indexUpdate {rep = Delayed} is f arr =
   \is' =>
     let x = arr is'
-    in if eqNP is is' then f x else x
+    in if (==) @{eqCoords} is is' then f x else x
 
 export
 indexRange : {rep,s : _} -> RepConstraint rep a => (rs : CoordsRange s) -> PrimArray rep s a

src/Data/NumIdr/PrimArray/Boxed.idr

@@ -62,7 +62,7 @@ unsafeFromIns s ins = arrayAction s $ \arr =>
 export
 create : {o : _} -> (s : Vect rk Nat) -> (Nat -> a) -> PrimArrayBoxed o s a
 create s f = arrayAction s $ \arr =>
-    for_ [0..pred (product s)] $ \i => arrayDataSet i (f i) arr
+    for_ (range 0 (product s)) $ \i => arrayDataSet i (f i) arr
 
 ||| Index into a primitive array. This function is unsafe, as it performs no
 ||| boundary check on the index given.
@@ -88,7 +88,7 @@ export
 bytesToBoxed : {s : _} -> ByteRep a => PrimArrayBytes o s -> PrimArrayBoxed o s a
 bytesToBoxed (MkPABytes sts buf) = MkPABoxed sts $ unsafePerformIO $ do
   arr <- newArrayData (product s) (believe_me ())
-  for_ [0..pred (product s)] $ \i => arrayDataSet i !(readBuffer i buf) arr
+  for_ (range 0 (product s)) $ \i => arrayDataSet i !(readBuffer i buf) arr
   pure arr
 
 export
@@ -96,7 +96,7 @@ boxedToBytes : {s : _} -> ByteRep a => PrimArrayBoxed o s a -> PrimArrayBytes o
 boxedToBytes @{br} (MkPABoxed sts arr) = MkPABytes sts $ unsafePerformIO $ do
   Just buf <- newBuffer $ cast (product s * bytes @{br})
     | Nothing => die "Cannot create array"
-  for_ [0..pred (product s)] $ \i => writeBuffer i !(arrayDataGet i arr) buf
+  for_ (range 0 (product s)) $ \i => writeBuffer i !(arrayDataGet i arr) buf
   pure buf
 
 
@@ -115,7 +115,7 @@ foldl f z (MkPABoxed sts arr) =
   if product s == 0 then z
   else unsafePerformIO $ do
     ref <- newIORef z
-    for_ [0..pred $ product s] $ \n => do
+    for_ (range 0 (product s)) $ \n => do
         x <- readIORef ref
         y <- arrayDataGet n arr
         writeIORef ref (f x y)
@@ -127,7 +127,7 @@ foldr f z (MkPABoxed sts arr) =
   if product s == 0 then z
   else unsafePerformIO $ do
     ref <- newIORef z
-    for_ [pred $ product s..0] $ \n => do
+    for_ (range 0 (product s)) $ \n => do
          x <- arrayDataGet n arr
          y <- readIORef ref
          writeIORef ref (f x y)
@@ -137,99 +137,3 @@ export
 traverse : {o,s : _} -> Applicative f => (a -> f b) -> PrimArrayBoxed o s a -> f (PrimArrayBoxed o s b)
 traverse f = map (Boxed.fromList _) . traverse f . foldr (::) []
 
-{-
-export
-updateAt : Nat -> (a -> a) -> PrimArrayBoxed o s a -> PrimArrayBoxed o s a
-updateAt n f arr = if n >= length arr then arr else
-  unsafePerformIO $ do
-    let cpy = copy arr
-    x <- arrayDataGet n cpy.content
-    arrayDataSet n (f x) cpy.content
-    pure cpy
-
-export
-unsafeUpdateInPlace : Nat -> (a -> a) -> PrimArrayBoxed a -> PrimArrayBoxed a
-unsafeUpdateInPlace n f arr = unsafePerformIO $ do
-  x <- arrayDataGet n arr.content
-  arrayDataSet n (f x) arr.content
-  pure arr
-
-||| Convert a primitive array to a list.
-export
-toList : PrimArrayBoxed a -> List a
-toList arr = iter (length arr) []
-  where
-    iter : Nat -> List a -> List a
-    iter Z acc = acc
-    iter (S n) acc = let el = index n arr
-                     in  iter n (el :: acc)
-
-||| Construct a primitive array from a list.
-export
-fromList : List a -> PrimArrayBoxed a
-fromList xs = create (length xs)
-    (\n => assert_total $ fromJust $ getAt n xs)
-  where
-    partial
-    fromJust : Maybe a -> a
-    fromJust (Just x) = x
-
-
-
-export
-unsafeZipWith : (a -> b -> c) -> PrimArrayBoxed a -> PrimArrayBoxed b -> PrimArrayBoxed c
-unsafeZipWith f a b = create (length a) (\n => f (index n a) (index n b))
-
-export
-unsafeZipWith3 : (a -> b -> c -> d) ->
-                 PrimArrayBoxed a -> PrimArrayBoxed b -> PrimArrayBoxed c -> PrimArrayBoxed d
-unsafeZipWith3 f a b c = create (length a) (\n => f (index n a) (index n b) (index n c))
-
-export
-unzipWith : (a -> (b, c)) -> PrimArrayBoxed a -> (PrimArrayBoxed b, PrimArrayBoxed c)
-unzipWith f arr = (map (fst . f) arr, map (snd . f) arr)
-
-export
-unzipWith3 : (a -> (b, c, d)) -> PrimArrayBoxed a -> (PrimArrayBoxed b, PrimArrayBoxed c, PrimArrayBoxed d)
-unzipWith3 f arr = (map ((\(x,_,_) => x) . f) arr,
-                          map ((\(_,y,_) => y) . f) arr,
-                          map ((\(_,_,z) => z) . f) arr)
-
-
-export
-foldl : (b -> a -> b) -> b -> PrimArrayBoxed a -> b
-foldl f z (MkPABoxed size arr) =
-  if size == 0 then z
-  else unsafePerformIO $ do
-    ref <- newIORef z
-    for_ [0..pred size] $ \n => do
-        x <- readIORef ref
-        y <- arrayDataGet n arr
-        writeIORef ref (f x y)
-    readIORef ref
-
-export
-foldr : (a -> b -> b) -> b -> PrimArrayBoxed a -> b
-foldr f z (MkPABoxed size arr) =
-  if size == 0 then z
-  else unsafePerformIO $ do
-    ref <- newIORef z
-    for_ [pred size..0] $ \n => do
-         x <- arrayDataGet n arr
-         y <- readIORef ref
-         writeIORef ref (f x y)
-    readIORef ref
-
-export
-traverse : Applicative f => (a -> f b) -> PrimArrayBoxed a -> f (PrimArrayBoxed b)
-traverse f = map fromList . traverse f . toList
-
-
-||| Compares two primitive arrays for equal elements. This function assumes the
-||| arrays have the same length; it must not be used in any other case.
-export
-unsafeEq : Eq a => PrimArrayBoxed a -> PrimArrayBoxed a -> Bool
-unsafeEq a b = unsafePerformIO $
-                 map (concat @{All}) $ for [0..pred (arraySize a)] $
-                 \n => (==) <$> arrayDataGet n (content a) <*> arrayDataGet n (content b)
--}

src/Data/NumIdr/PrimArray/Bytes.idr

@@ -48,8 +48,7 @@ bufferAction @{br} s action = fromBuffer s $ unsafePerformIO $ do
 export
 constant : {o : _} -> ByteRep a => (s : Vect rk Nat) -> a -> PrimArrayBytes o s
 constant @{br} s x = bufferAction @{br} s $ \buf => do
-  let size = product s
-  for_ [0..pred size] $ \i => writeBuffer i x buf
+  for_ (range 0 (product s)) $ \i => writeBuffer i x buf
 
 export
 unsafeDoIns : ByteRep a => List (Nat, a) -> PrimArrayBytes o s -> IO ()
@@ -63,8 +62,7 @@ unsafeFromIns @{br} s ins = bufferAction @{br} s $ \buf =>
 export
 create : {o : _} -> ByteRep a => (s : Vect rk Nat) -> (Nat -> a) -> PrimArrayBytes o s
 create @{br} s f = bufferAction @{br} s $ \buf => do
-  let size = product s
-  for_ [0..pred size] $ \i => writeBuffer i (f i) buf
+  for_ (range 0 (product s)) $ \i => writeBuffer i (f i) buf
 
 export
 index : ByteRep a => Nat -> PrimArrayBytes o s -> a
@@ -102,7 +100,7 @@ foldl f z (MkPABytes sts buf) =
   if product s == 0 then z
   else unsafePerformIO $ do
     ref <- newIORef z
-    for_ [0..pred $ product s] $ \n => do
+    for_ (range 0 (product s)) $ \n => do
         x <- readIORef ref
         y <- readBuffer n buf
         writeIORef ref (f x y)

src/Data/NumIdr/PrimArray/Delayed.idr

@@ -2,9 +2,9 @@ module Data.NumIdr.PrimArray.Delayed
 
 import Data.List
 import Data.Vect
-import Data.NP
 import Data.NumIdr.Array.Rep
 import Data.NumIdr.Array.Coords
+import Data.NumIdr.PrimArray.Linked
 
 %default total
 
@@ -22,8 +22,8 @@ export
 indexRange : {s : _} -> (rs : CoordsRange s) -> PrimArrayDelayed s a -> PrimArrayDelayed (newShape rs) a
 indexRange [] v = v
 indexRange (r :: rs) v with (cRangeToList r)
-  _ | Left i = indexRange rs (\is' => v (believe_me i :: is'))
-  _ | Right is = \(i::is') => indexRange rs (\is'' => v (believe_me (Vect.index i (fromList is)) :: is'')) is'
+  _ | Left i = indexRange rs (\is' => v (assertFin i :: is'))
+  _ | Right is = \(i::is') => indexRange rs (\is'' => v (assertFin (Vect.index i (fromList is)) :: is'')) is'
 
 export
 indexSetRange : {s : _} -> (rs : CoordsRange s) -> PrimArrayDelayed (newShape rs) a

src/Data/NumIdr/PrimArray/Linked.idr

@@ -2,7 +2,6 @@ module Data.NumIdr.PrimArray.Linked
 
 import Data.Vect
 import Data.DPair
-import Data.NP
 import Data.NumIdr.Array.Rep
 import Data.NumIdr.Array.Coords
 
@@ -28,15 +27,17 @@ export
 indexRange : {s : _} -> (rs : CoordsRange s) -> Vects s a -> Vects (newShape rs) a
 indexRange [] v = v
 indexRange (r :: rs) v with (cRangeToList r)
-  _ | Left i = indexRange rs (Vect.index (believe_me i) v)
-  _ | Right is = believe_me $ map (\i => indexRange rs (Vect.index (believe_me i) v)) is
+  _ | Left i = indexRange rs (Vect.index (assertFin i) v)
+  _ | Right is = assert_total $
+        case toVect _ (map (\i => indexRange rs (Vect.index (assertFin i) v)) is) of
+          Just v => v
 
 export
 indexSetRange : {s : _} -> (rs : CoordsRange s) -> Vects (newShape rs) a -> Vects s a -> Vects s a
 indexSetRange {s=[]} [] rv _ = rv
 indexSetRange {s=_::_} (r :: rs) rv v with (cRangeToList r)
-  _ | Left i = updateAt (believe_me i) (indexSetRange rs rv) v
-  _ | Right is = foldl (\v,i => updateAt (believe_me i) (indexSetRange rs (Vect.index (believe_me i) rv)) v) v is
+  _ | Left i = updateAt (assertFin i) (indexSetRange rs rv) v
+  _ | Right is = foldl (\v,i => updateAt (assertFin i) (indexSetRange rs (Vect.index (assertFin i) rv)) v) v is
 
 
 export

src/Data/NumIdr/Transform/Point.idr

@@ -176,7 +176,7 @@ Traversable (Point n) where
 
 export
 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
 Cast a b => Cast (Point n a) (Point n b) where

src/Data/NumIdr/Transform/Translation.idr

@@ -25,12 +25,13 @@ 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
 
-||| 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
+
+
+||| Construct a translation given a vector.
+export
+translate : Num a => Vector n a -> Translation n a
+translate v = unsafeMkTrans (translationH v)

src/Data/NumIdr/Vector.idr

@@ -40,7 +40,8 @@ vector v = rewrite sym (lengthCorrect v)
 ||| Convert a vector into a `Vect`.
 export
 toVect : Vector n a -> Vect n a
-toVect v = believe_me $ Vect.fromList $ Prelude.toList v
+toVect v = replace {p = flip Vect a} (believe_me $ Refl {x=()}) $
+  Vect.fromList $ Prelude.toList v
 
 ||| Return the `i`-th basis vector.
 export
@@ -185,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] = elements v1
-                  [x, y, z] = elements v2
+cross v1 v2 = let [a, b, c] = toVect v1
+                  [x, y, z] = toVect v2
               in  vector [b*z - c*y, c*x - a*z, a*y - b*x]
 
 ||| Calculate the triple product of the three vectors.