add docs setup

docs/Project.toml

@@ -0,0 +1,10 @@
+[deps]
+DualArrays = "429e4e16-f749-45f3-beec-30742fae38ce"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+
+[compat]
+Documenter = "1"
+DualArrays = "0.2.0"
+
+[sources]
+DualArrays = { path = ".." }

docs/make.jl

@@ -0,0 +1,22 @@
+using Pkg
+
+Pkg.activate(@__DIR__)
+Pkg.instantiate()
+
+using Documenter
+using DualArrays
+
+# Setup for doctests in docstrings
+DocMeta.setdocmeta!(DualArrays, :DocTestSetup, :(using DualArrays))
+
+makedocs(;
+    format = Documenter.HTML(
+        canonical = "https://github.com/dlfivefifty/DualArrays.jl",
+    ),
+    pages = [
+        "Home" => "index.md",
+        ],
+    sitename = "DualArrays.jl",
+)
+
+deploydocs(; repo = "github.com/dlfivefifty/DualArrays.jl")

docs/src/index.md

@@ -0,0 +1,90 @@
+# DualArrays.jl
+
+Documentation for `DualArrays.jl`.
+
+## Dual
+
+```@docs
+DualArrays.Dual
+```
+
+```@example
+using DualArrays
+
+a = Dual(2.0, [1.0, 2.0, 3.0])
+b = Dual(3.0, [4.0, 5.0, 6.0])
+
+a * b
+```
+
+## DualArray
+
+```@docs
+DualArrays.DualArray
+DualArrays.DualVector
+DualArrays.DualMatrix
+```
+
+```@example
+using DualArrays
+
+v = DualVector([1.0, 2.0, 3.0], [1 2 3; 4 5 6; 7 8 9])
+v[2]
+```
+
+## ArrayOperator
+
+```@docs
+DualArrays.ArrayOperator
+```
+
+### Broadcasting with ArrayOperators
+
+```@eval
+using DualArrays
+import Base.Docs
+import Markdown
+
+Markdown.parse(sprint(show, Docs.doc(DualArrays.ArrayOperatorBroadcastStyle)))
+```
+
+```@example
+using DualArrays
+
+t = ArrayOperator{1}([1 2 3; 4 5 6; 7 8 9])
+
+t .+ 1
+```
+
+### Multiplication with ArrayOperators
+
+```@eval
+using DualArrays
+import Base.Docs
+import Markdown
+
+Markdown.parse(sprint(show, Docs.doc(DualArrays._contract)))
+```
+
+```@example
+using DualArrays
+
+t = ArrayOperator{1}([1 2 3; 4 5 6; 7 8 9])
+
+t * t
+```
+
+## Utilities
+
+```@docs
+DualArrays.jacobian
+```
+
+```@example
+using DualArrays
+
+f(x) = sin.(x) .+ x .* x
+
+J = jacobian(f, [1.0, 2.0, 3.0])
+J.data
+```

src/DualArrays.jl

@@ -7,6 +7,8 @@ This package provides:
 - `Dual`: A dual number type for storing values and their derivatives
 - `DualVector`: A vector of dual numbers represented with a Jacobian matrix
 - `DualMatrix`: A matrix of dual numbers represented with a Jacobian tensor`
+- `ArrayOperator`: A data structure that stores an array equipped with an output and input dimension.
+    Useful for tensor computations involving higher order DualArrays
 Differentiation rules are mostly provided by ChainRules.jl.
 """
 module DualArrays

src/arithmetic.jl

@@ -3,35 +3,29 @@
 # Tensor transpose
 transpose(t::ArrayOperator{N, M}) where {N, M} = ArrayOperator{M}(transpose(t.data))
 
-"""
-Addition/Subtraction of DualVectors.
-"""
-
+# Addition/Subtraction of DualVectors.
 for op in (:+, :-)
     @eval $op(x::DualVector, y::DualVector) = DualVector($op(x.value, y.value), $op(x.jacobian, y.jacobian))
     @eval $op(x::DualVector, y::AbstractVector) = DualVector($op(x.value, y), x.jacobian)
     @eval $op(x::AbstractVector, y::DualVector) = DualVector($op(x, y.value), y.jacobian)
 end
 
-"""
-Matrix multiplication with a DualVector.
-"""
-*(x::AbstractMatrix, y::DualVector) = DualVector(x * y.value, x * y.jacobian)
-
-"""
-this section attempts to define broadcasting rules on DualVectors for functions
-that either:
 
-- take a single real argument (function applied to each element)
-- are binary operations (binary operation applied to scalar and each element)
-
-We use DiffRules.jl for symbolic derivatives and define our overloads accordingly.
-"""
+# Matrix multiplication with a DualVector.
+*(x::AbstractMatrix, y::DualVector) = DualVector(x * y.value, x * y.jacobian)
 
+###
+# this section attempts to define broadcasting rules on DualVectors for functions
+# that either:
+#
+# - take a single real argument (function applied to each element)
+# - are binary operations (binary operation applied to scalar and each element)
+#
+# We use DiffRules.jl for symbolic derivatives and define our overloads accordingly.
+###
+
+# Helper function: Given an n-ary function f, return its partial derivatives as a tuple of function.
 function diff_fn(f, n)
-    """
-    Helper function: Given an n-ary function f, return its partial derivatives as a tuple of function.
-    """
     syms = ntuple(_ -> gensym(), n)
     d = DiffRules.diffrule(:Base, f, syms...)
     d = d isa Tuple ? d : (d,)
@@ -121,12 +115,10 @@ LinearAlgebra.dot(x::AbstractVector, y::DualVector) = Dual(dot(x, y.value), tran
 # solve
 \(x::AbstractMatrix, y::DualVector) = DualVector(x \ y.value, ArrayOperator{1}(x \ y.jacobian.data))
 
-"""
------------------------
-Multiplication with ArrayOperators
------------------------
 
-* for ArrayOperators generalises from matrix/vector
+
+"""
+Multiplication (*) for ArrayOperators generalises from matrix/vector
 multiplication and uses tensor contraction provided by
 TensorOperations.jl. Let an (A, B) ArrayOperator denote
 an ArrayOperator with output dimension A and input dimension B
@@ -135,21 +127,22 @@ allow multiplication between an (A, B) ArrayOperator and a (B, C)
 ArrayOperator, with the result being an (A, C) ArrayOperator.
 This generalises from:
 
--An inner product: a row vector * a column vector -> a scalar:
-    (0, 1) * (1, 0) -> (0, 0)
+- An inner product: a row vector * a column vector -> a scalar:
+
+            (0, 1) * (1, 0) -> (0, 0)
 
--An outer product: a column vector * a row vector -> a matrix:
-    (1, 0) * (0, 1) -> (1, 1)
+- An outer product: a column vector * a row vector -> a matrix:
+
+            (1, 0) * (0, 1) -> (1, 1)
 
 - matrix/vector multiplication: a matrix * a column vector -> a column vector:
-    (1, 1) * (1, 0) -> (1, 0)
+
+            (1, 1) * (1, 0) -> (1, 0)
 
 We can treat multiplication between ArrayOperators and regular arrays as multiplication
 between two ArrayOperators: given an (A, B) ArrayOperator and an L-array, we can fix it as
 having input dimension B and infer C = L - B. We use this contraction rule to multiply the two.
 """
-
-# Helper function to perform tensor contraction between x and y
 function _contract(x, y, A, B, C)
     x_idx = (ntuple(i -> i, A), ntuple(i -> i + A, B))
     y_idx = (ntuple(i -> i, B), ntuple(i -> i + B, C))
@@ -158,6 +151,7 @@ function _contract(x, y, A, B, C)
     return TensorOperations.tensorcontract(x, x_idx, false, y, y_idx, false, ret_idx, 1)
 end
 
+
 function *(x::ArrayOperator{A, B}, y::ArrayOperator{B, C}) where {A, B, C}
     return ArrayOperator{A}(_contract(x.data, y.data, A, B, C))
 end

src/indexing.jl

@@ -4,8 +4,6 @@ using ArrayLayouts, FillArrays, LinearAlgebra, SparseArrays
 
 
 sparse_getindex(a...) = layout_getindex(a...)
-
-# TODO: should we move these?
 sparse_getindex(D::Diagonal, k::Integer, ::Colon) = OneElement(D.diag[k], k, size(D, 2))
 sparse_getindex(D::Diagonal, ::Colon, j::Integer) = OneElement(D.diag[j], j, size(D, 1))
 sparse_getindex(d::DualArray, args...) = d[args...]
@@ -49,16 +47,13 @@ getindex(t::ArrayOperator, i::Vararg{Int}) = getindex(t.data, i...)
 # Indexing with only slices/ranges returns a similar tensor
 getindex(t::ArrayOperator{N, M}, i::Vararg{Union{Colon, UnitRange}}) where {N, M} = ArrayOperator{N}(sparse_getindex(t.data, i...))
 
-"""
-Extract a single Dual number from a DualArray.
-"""
+
+# Extract a single Dual number from a DualArray.
 function getindex(x::DualArray, args::Vararg{Int})
     Dual(x.value[args...], x.jacobian[args, :])
 end
 
-"""
-Extract a sub-DualVector from a DualVector using a range.
-"""
+# Extract a sub-DualVector from a DualVector using a range.
 function getindex(x::DualVector, y::UnitRange)
     newval = x.value[y]
     newjac = x.jacobian[y, :]