From 18dd384f49519cce0e35b4e611ca6c54c330de43 Mon Sep 17 00:00:00 2001
From: Max Vassiliev <vassiliev.max@gmail.com>
Date: Tue, 4 Nov 2025 20:45:37 +0000
Subject: [PATCH 1/7] initial commit: create readme, organise code, fix bugs

---
 Project.toml               |  23 +++++++-
 README.md                  |  39 +++++++++++++-
 examples/newtonpendulum.jl |   4 +-
 src/DualArrays.jl          | 108 +++++--------------------------------
 src/arithmetic.jl          |  45 ++++++++++++++++
 src/indexing.jl            |  41 ++++++++++++++
 src/types.jl               |  61 +++++++++++++++++++++
 src/utilities.jl           |  64 ++++++++++++++++++++++
 8 files changed, 286 insertions(+), 99 deletions(-)
 create mode 100644 src/arithmetic.jl
 create mode 100644 src/indexing.jl
 create mode 100644 src/types.jl
 create mode 100644 src/utilities.jl

diff --git a/Project.toml b/Project.toml
index 0375150..815faea 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,15 +1,34 @@
 name = "DualArrays"
 uuid = "429e4e16-f749-45f3-beec-30742fae38ce"
-authors = ["Sheehan Olver <solver@mac.com>"]
 version = "0.1.0"
+authors = ["Sheehan Olver <solver@mac.com>"]
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
+DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+
+[compat]
+ChainRules = "1.72.6"
+DiffRules = "1.15.1"
+DifferentialEquations = "7.17.0"
+ForwardDiff = "1.2.2"
+Plots = "1.41.1"
+
+[extras]
+BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
+DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+
+[targets]
+dev = ["Revise"]
+examples = ["Plots", "DifferentialEquations", "ForwardDiff", "BandedMatrices"]
+test = ["Test", "ForwardDiff", "BandedMatrices"]
diff --git a/README.md b/README.md
index 8cc9eeb..a3976fc 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,39 @@
 # DualArrays.jl
- A package for working with arrays of dual numbers with sparsity
+DualArrays.jl is a package that provides the `DualVector` structure for use in forward mode automatic differentiation (autodiff). Existing forward mode autodiff implementations such as `ForwardDiff.jl` make use of a `Dual` structure representing Dual numbers as follows:
+
+```julia
+struct Dual{T,V<:Real,N} <: Real
+    value::V
+    partials::Partials{N,V}
+end
+```
+
+where partials is defined as follows:
+
+```julia
+struct Partials{N,V} <: AbstractVector{V}
+    values::NTuple{N,V}
+end
+```
+
+There are some limitations of this when differentiating vector valued functions, as having a vector of `Dual`s cannot properly exploit sparse Jacobian structures such as banded structures that appear when solving systems of ODEs. `DualArrays.jl` introduces a new `DualVector` type as follows:
+
+```julia
+struct DualVector{T, M <: AbstractMatrix{T}} <: AbstractVector{Dual{T}}
+    value::Vector{T}
+    jacobian::M
+end
+```
+
+with a similarly defined `Dual` type. This allows for jacobians to be any sparse matrix structure defined in the Julia linear algebra ecosystem. An efficient implementation of differentiating a function might be as follows:
+
+```julia
+using FillArrays
+
+function gradient(f::Function, x::Vector)
+    dx = DualVector(x, Eye(length(x)))
+    return f(dx).jacobian
+end
+```
+
+See the examples folder for more use cases
diff --git a/examples/newtonpendulum.jl b/examples/newtonpendulum.jl
index 81d727b..7802ef8 100644
--- a/examples/newtonpendulum.jl
+++ b/examples/newtonpendulum.jl
@@ -48,8 +48,8 @@ end
 x0 = zeros(Float64, N)
 
 #Solve and plot both solution and LHS ('deviation' from system)
-@time sol = newton_method_forwarddiff(f, x0, 100)
-@time sol = newton_method_dualvector(f, x0, 100)
+@time sol = newton_method_forwarddiff(f, x0, 10)
+@time sol = newton_method_dualvector(f, x0, 10)
 plot(0:ts:Tmax, [a; sol; b])
 plot!(0:ts:Tmax, [0; f(sol); 0])
 
diff --git a/src/DualArrays.jl b/src/DualArrays.jl
index e4074d5..c2c4ed7 100644
--- a/src/DualArrays.jl
+++ b/src/DualArrays.jl
@@ -1,103 +1,23 @@
-module DualArrays
-export DualVector
-
-import Base: +, ==, getindex, size, broadcast, axes, broadcasted, show, sum,
-             vcat, convert, *
-using LinearAlgebra, ArrayLayouts, BandedMatrices, FillArrays
-
-struct Dual{T, Partials <: AbstractVector{T}} <: Real
-    value::T
-    partials::Partials
-end
-
-==(a::Dual, b::Dual) = a.value == b.value && a.partials == b.partials
-
-sparse_getindex(a...) = layout_getindex(a...)
-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))
-
 """
-reprents a vector of duals given by
-    
-    values + jacobian * [ε_1,…,ε_n].
-
-For now the entries just return the values.
-"""
-
-struct DualVector{T, M <: AbstractMatrix{T}} <: AbstractVector{Dual{T}}
-    value::Vector{T}
-    jacobian::M
-    function DualVector(value::Vector{T},jacobian::M) where {T, M <: AbstractMatrix{T}}
-        if(size(jacobian)[1] != length(value))
-            x,y = length(value),size(jacobian)[1]
-            throw(ArgumentError("vector length must match number of rows in jacobian.\n
-            vector length: $x \n
-            no. of jacobian rows: $y"))
-        end
-        new{T,M}(value,jacobian)
-    end
-end
-
-function DualVector(value::AbstractVector, jacobian::AbstractMatrix)
-    T = promote_type(eltype(value), eltype(jacobian))
-    DualVector(convert(Vector{T}, value), convert(AbstractMatrix{T}, jacobian))
-end
-
+    DualArrays
 
+A Julia package for efficient automatic differentiation using dual numbers and dual arrays.
 
-function getindex(x::DualVector, y::Int)
-    Dual(x.value[y], sparse_getindex(x.jacobian,y,:))
-end
-
-function getindex(x::DualVector, y::UnitRange)
-    newval = x.value[y]
-    newjac = sparse_getindex(x.jacobian,y,:)
-    DualVector(newval, newjac)
-end
-size(x::DualVector) = length(x.value)
-axes(x::DualVector) = axes(x.value)
-+(x::DualVector,y::DualVector) = DualVector(x.value + y.value, x.jacobian + y.jacobian)
-*(x::AbstractMatrix, y::DualVector) = DualVector(x * y.value, x * y.jacobian)
-
-broadcasted(::typeof(sin),x::DualVector) = DualVector(sin.(x.value),Diagonal(cos.(x.value))*x.jacobian)
-
-function broadcasted(::typeof(*),x::DualVector,y::DualVector)
-    newval = x.value .* y.value
-    newjac = x.value .* y.jacobian + y.value .* x.jacobian
-    DualVector(newval,newjac)
-end
-
-function sum(x::DualVector)
-    n = length(x.value)
-    Dual(sum(x.value), vec(sum(x.jacobian; dims=1)))
-end
-
-_jacobian(d::Dual) = permutedims(d.partials)
-_jacobian(d::DualVector, ::Int) = d.jacobian
-_jacobian(x::Number, N::Int) = zeros(typeof(x), 1, N)
+This package provides:
+- `Dual`: A dual number type for storing values and their derivatives
+- `DualVector`: A vector of dual numbers represented efficiently with a Jacobian matrix
+"""
+module DualArrays
 
-_value(d::DualVector) = d.value
-_value(x::Number) = x
+export DualVector, Dual
 
-function vcat(x::Union{Dual, DualVector}...)
-    if length(x) == 1
-        return x[1]
-    end
-    value = vcat((d.value for d in x)...)
-    jacobian = vcat((_jacobian(d) for d in x)...)
-    DualVector(value,jacobian)
-end
+import Base: +, ==, getindex, size, axes, broadcasted, show, sum, vcat, convert, *
 
-_size(x::Real) = 1
-_size(x::DualVector) = size(x)
+using LinearAlgebra, ArrayLayouts, FillArrays
 
-function vcat(x::Union{Real, DualVector}...)
-    cols = max((_size(i) for i in x)...)
-    val = vcat((_value(i) for i in x)...)
-    jac = vcat((_jacobian(i, cols) for i in x)...)
-    DualVector(val, jac)
+include("types.jl")
+include("indexing.jl")
+include("arithmetic.jl")
+include("utilities.jl")
 end
 
-show(io::IO,::MIME"text/plain", x::DualVector) = (print(io,x.value); print(io," + "); print(io,x.jacobian);print("𝛜"))
-end
-# module DualArrays
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
new file mode 100644
index 0000000..7bfc690
--- /dev/null
+++ b/src/arithmetic.jl
@@ -0,0 +1,45 @@
+# Arithmetic operations for DualArrays.jl
+
+using LinearAlgebra, ChainRules
+
+"""
+Addition of DualVectors.
+"""
+Base.:+(x::DualVector, y::DualVector) = DualVector(x.value + y.value, x.jacobian + y.jacobian)
+Base.:+(x::DualVector, y::AbstractVector) = DualVector(x.value + y, x.jacobian)
+Base.:+(x::AbstractVector, y::DualVector) = DualVector(x + y.value, y.jacobian)
+
+"""
+    *(x::AbstractMatrix, y::DualVector)
+
+Matrix multiplication with a DualVector.
+"""
+Base.:*(x::AbstractMatrix, y::DualVector) = DualVector(x * y.value, x * y.jacobian)
+
+"""
+    broadcasted(::typeof(sin), x::DualVector)
+
+Broadcasted sine function for DualVectors using the chain rule.
+"""
+
+Base.broadcasted(::typeof(sin), x::DualVector) = 
+    DualVector(sin.(x.value), Diagonal(cos.(x.value)) * x.jacobian)
+
+"""
+    broadcasted(::typeof(*), x::DualVector, y::DualVector)
+
+Broadcasted multiplication of two DualVectors using the product rule.
+"""
+function Base.broadcasted(::typeof(*), x::DualVector, y::DualVector)
+    newval = x.value .* y.value
+    newjac = Diagonal(x.value) * y.jacobian + Diagonal(y.value) * x.jacobian
+    DualVector(newval, newjac)
+end
+
+for op in (:-, :/)
+    fdef = quote
+        function (op)(args::Vararg{Union{DualVector, AbstractVector}, 2})
+            _args = (NoTangent(), )
+        end
+    end
+end
\ No newline at end of file
diff --git a/src/indexing.jl b/src/indexing.jl
new file mode 100644
index 0000000..e2556ca
--- /dev/null
+++ b/src/indexing.jl
@@ -0,0 +1,41 @@
+# Indexing operations for DualArrays.jl
+
+using ArrayLayouts, FillArrays
+
+sparse_getindex(a...) = layout_getindex(a...)
+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))
+
+"""
+    getindex(x::DualVector, y::Int)
+
+Extract a single Dual number from a DualVector at position y.
+"""
+function Base.getindex(x::DualVector, y::Int)
+    Dual(x.value[y], sparse_getindex(x.jacobian, y, :))
+end
+
+"""
+    getindex(x::DualVector, y::UnitRange)
+
+Extract a sub-DualVector from a DualVector using a range.
+"""
+function Base.getindex(x::DualVector, y::UnitRange)
+    newval = x.value[y]
+    newjac = sparse_getindex(x.jacobian, y, :)
+    DualVector(newval, newjac)
+end
+
+"""
+    size(x::DualVector)
+
+Return the size of the DualVector (length of the value vector).
+"""
+Base.size(x::DualVector) = (length(x.value),)
+
+"""
+    axes(x::DualVector)
+
+Return the axes of the DualVector.
+"""
+Base.axes(x::DualVector) = axes(x.value)
\ No newline at end of file
diff --git a/src/types.jl b/src/types.jl
new file mode 100644
index 0000000..09af0c8
--- /dev/null
+++ b/src/types.jl
@@ -0,0 +1,61 @@
+# Core type definitions for DualArrays.jl
+
+"""
+    Dual{T, Partials <: AbstractVector{T}} <: Real
+
+A dual number type that stores a value and its partials (derivatives).
+
+# Fields
+- `value::T`: The primal value
+- `partials::Partials`: The partial derivatives as a vector
+"""
+struct Dual{T, Partials <: AbstractVector{T}} <: Real
+    value::T
+    partials::Partials
+end
+
+"""
+    DualVector{T, M <: AbstractMatrix{T}} <: AbstractVector{Dual{T}}
+
+Represents a vector of dual numbers given by:
+    
+    values + jacobian * [ε₁, …, εₙ]
+
+For now the entries just return the values when indexed.
+
+# Fields
+- `value::Vector{T}`: The primal values
+- `jacobian::M`: The Jacobian matrix containing partial derivatives
+
+# Constructor
+    DualVector(value::Vector{T}, jacobian::M) where {T, M <: AbstractMatrix{T}}
+
+Constructs a DualVector, ensuring that the vector length matches the number of rows in the Jacobian.
+"""
+struct DualVector{T, M <: AbstractMatrix{T}} <: AbstractVector{Dual{T}}
+    value::Vector{T}
+    jacobian::M
+    
+    function DualVector(value::Vector{T}, jacobian::M) where {T, M <: AbstractMatrix{T}}
+        if size(jacobian, 1) != length(value)
+            x, y = length(value), size(jacobian, 1)
+            throw(ArgumentError("vector length must match number of rows in jacobian.\n" *
+                               "vector length: $x\n" *
+                               "no. of jacobian rows: $y"))
+        end
+        new{T, M}(value, jacobian)
+    end
+end
+
+"""
+    DualVector(value::AbstractVector, jacobian::AbstractMatrix)
+
+Convenience constructor that promotes types to be compatible.
+"""
+function DualVector(value::AbstractVector, jacobian::AbstractMatrix)
+    T = promote_type(eltype(value), eltype(jacobian))
+    DualVector(convert(Vector{T}, value), convert(AbstractMatrix{T}, jacobian))
+end
+
+# Basic equality for Dual numbers
+==(a::Dual, b::Dual) = a.value == b.value && a.partials == b.partials
\ No newline at end of file
diff --git a/src/utilities.jl b/src/utilities.jl
new file mode 100644
index 0000000..5217b16
--- /dev/null
+++ b/src/utilities.jl
@@ -0,0 +1,64 @@
+# Utility functions for DualArrays.jl
+
+"""
+    sum(x::DualVector)
+
+Sum all elements of a DualVector, returning a single Dual number.
+"""
+function Base.sum(x::DualVector)
+    n = length(x.value)
+    Dual(sum(x.value), vec(sum(x.jacobian; dims=1)))
+end
+
+# Helper functions for vcat operations
+_jacobian(d::Dual) = permutedims(d.partials)
+_jacobian(d::DualVector) = d.jacobian
+_jacobian(d::DualVector, ::Int) = d.jacobian
+_jacobian(x::Number, N::Int) = zeros(typeof(x), 1, N)
+
+_value(v::AbstractVector) = v
+_value(d::DualVector) = d.value
+_value(d::Dual) = d.value
+_value(x::Number) = x
+
+_size(x::Real) = 1
+_size(x::DualVector) = length(x.value)
+
+"""
+    vcat(x::Union{Dual, DualVector}...)
+
+Vertically concatenate Dual numbers and DualVectors.
+"""
+function Base.vcat(x::Union{Dual, DualVector}...)
+    if length(x) == 1
+        return x[1]
+    end
+    value = vcat((d.value for d in x)...)
+    jacobian = vcat((_jacobian(d) for d in x)...)
+    DualVector(value, jacobian)
+end
+
+"""
+    vcat(x::Union{Real, DualVector}...)
+
+Vertically concatenate Real numbers and DualVectors.
+"""
+
+function Base.vcat(x::Union{Real, DualVector}...)
+    # Avoid stack overflow
+    if all(i -> i isa Real, x)
+        return Base.cat(x..., dims=1)
+    end
+    cols = maximum((_size(i) for i in x))
+    val = vcat((_value(i) for i in x)...)
+    jac = vcat((_jacobian(i, cols) for i in x)...)
+    DualVector(val, jac)
+end
+
+"""
+    show(io::IO, ::MIME"text/plain", x::DualVector)
+
+Custom display method for DualVectors.
+"""
+Base.show(io::IO, ::MIME"text/plain", x::DualVector) = 
+    (print(io, x.value); print(io, " + "); print(io, x.jacobian); print(io, "𝛜"))
\ No newline at end of file

From 5bde0a81f0de4d7c6a4e9ab97bd13ab6def22412 Mon Sep 17 00:00:00 2001
From: Max Vassiliev <vassiliev.max@gmail.com>
Date: Mon, 17 Nov 2025 12:00:34 +0000
Subject: [PATCH 2/7] define broadcasting, clean up code

---
 Project.toml      |   6 +--
 src/DualArrays.jl |   4 +-
 src/arithmetic.jl | 134 ++++++++++++++++++++++++++++++++++++++++------
 src/indexing.jl   |   8 ---
 src/types.jl      |  10 ++--
 src/utilities.jl  |   8 ---
 6 files changed, 129 insertions(+), 41 deletions(-)

diff --git a/Project.toml b/Project.toml
index 815faea..81c651a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,13 +1,13 @@
 name = "DualArrays"
 uuid = "429e4e16-f749-45f3-beec-30742fae38ce"
-version = "0.1.0"
+version = "0.1.1"
 authors = ["Sheehan Olver <solver@mac.com>"]
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
-DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
+ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -15,7 +15,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 
 [compat]
 ChainRules = "1.72.6"
-DiffRules = "1.15.1"
+ChainRulesCore = "1.26.0"
 DifferentialEquations = "7.17.0"
 ForwardDiff = "1.2.2"
 Plots = "1.41.1"
diff --git a/src/DualArrays.jl b/src/DualArrays.jl
index c2c4ed7..7feb3d5 100644
--- a/src/DualArrays.jl
+++ b/src/DualArrays.jl
@@ -5,7 +5,9 @@ A Julia package for efficient automatic differentiation using dual numbers and d
 
 This package provides:
 - `Dual`: A dual number type for storing values and their derivatives
-- `DualVector`: A vector of dual numbers represented efficiently with a Jacobian matrix
+- `DualVector`: A vector of dual numbers represented with a Jacobian matrix
+
+Differentiation rules are mostly provided by ChainRules.jl.
 """
 module DualArrays
 
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index 7bfc690..f84838c 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -10,36 +10,138 @@ Base.:+(x::DualVector, y::AbstractVector) = DualVector(x.value + y, x.jacobian)
 Base.:+(x::AbstractVector, y::DualVector) = DualVector(x + y.value, y.jacobian)
 
 """
-    *(x::AbstractMatrix, y::DualVector)
-
 Matrix multiplication with a DualVector.
 """
 Base.:*(x::AbstractMatrix, y::DualVector) = DualVector(x * y.value, x * y.jacobian)
 
 """
-    broadcasted(::typeof(sin), x::DualVector)
+Broadcasted multiplication of two DualVectors using the product rule.
+"""
+function Base.broadcasted(::typeof(*), x::DualVector, y::DualVector)
+    newval = x.value .* y.value
+    newjac = Diagonal(x.value) * y.jacobian + Diagonal(y.value) * x.jacobian
+    DualVector(newval, newjac)
+end
 
-Broadcasted sine function for DualVectors using the chain rule.
 """
 
-Base.broadcasted(::typeof(sin), x::DualVector) = 
-    DualVector(sin.(x.value), Diagonal(cos.(x.value)) * x.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)
+
+This implementation is loosely based on
+https://juliadiff.org/ChainRulesOverloadGeneration.jl/dev/examples/forward_mode.html
 
 """
-    broadcasted(::typeof(*), x::DualVector, y::DualVector)
 
-Broadcasted multiplication of two DualVectors using the product rule.
 """
-function Base.broadcasted(::typeof(*), x::DualVector, y::DualVector)
-    newval = x.value .* y.value
-    newjac = Diagonal(x.value) * y.jacobian + Diagonal(y.value) * x.jacobian
-    DualVector(newval, newjac)
+Defines how a frule in ChainRules.jl for a scalar function f should be broadcasted over a DualVector.
+"""
+function broadcast_rule(f, d::DualVector)
+    val = similar(d.value)
+    jac = similar(d.jacobian)
+
+    @inbounds for (i, x) in pairs(d.value)
+        y, dy = ChainRules.frule((ChainRules.NoTangent(), d.jacobian[i, :]), f, x[i])
+        val[i] = y
+        jac[i, :] = dy
+    end
+
+    return DualVector(val, jac)
 end
 
-for op in (:-, :/)
-    fdef = quote
-        function (op)(args::Vararg{Union{DualVector, AbstractVector}, 2})
-            _args = (NoTangent(), )
+"""
+Defines how a frule in ChainRules.jl for a binary
+operation a (+) b on reals a, b should be broadcasted over:
+- a (+) d (d a DualVector)
+- d (+) a (d a DualVector)
+"""
+
+function broadcast_rule(f, d::DualVector, x::Real)
+    val = similar(d.value)
+    jac = similar(d.jacobian)
+
+    @inbounds for (i, y) in pairs(d.value)
+        yval, dy = ChainRules.frule(
+            (ChainRules.NoTangent(), d.jacobian[i, :], zero(d.jacobian[i, :])),
+             f, y, x)
+        val[i] = yval
+        jac[i, :] = dy
+    end
+
+    return DualVector(val, jac)
+end
+
+function broadcast_rule(f, x::Real, d::DualVector)
+    val = similar(d.value)
+    jac = similar(d.jacobian)
+
+    @inbounds for (i, y) in pairs(d.value)
+        yval, dy = ChainRules.frule(
+            (ChainRules.NoTangent(), zero(d.jacobian[i, :]), d.jacobian[i, :]),
+             f, x, y)
+        val[i] = yval
+        jac[i, :] = dy
+    end
+
+    return DualVector(val, jac)
+end
+
+# Define a lookup table typeof(function) => function
+function_lookup = Dict{DataType, Function}()
+
+for n in names(Base)
+    f = getfield(Base, n)
+    if f isa Function
+        function_lookup[typeof(f)] = f
+    end
+end
+
+# Get all applicable frules defined in ChainRules
+# and define broadcasted versions for DualVector using vector_rule
+
+frules = methods(ChainRules.frule)
+for f in frules
+    # get signatures for each function with a frule
+    sig = Base.unwrap_unionall(Base.tuple_type_tail(f.sig))
+    
+    # split into operation and args, filter for
+    # single argument functions and binary operations that can act on real numbers
+    op, args = sig.parameters[2], sig.parameters[3:end]
+
+    isconcretetype(op) || continue
+
+    # if it is a single argument function...
+    if length(args) == 1
+        args[1] isa Type || continue
+        Real <: args[1] || continue
+
+        try
+            fn = function_lookup[op]
+
+            @eval Base.broadcasted(x::$op, d::DualVector) = broadcast_rule($fn, d)
+        catch e
+            # skip over erroneous definitions. Uncomment the below for debugging.
+            # println("Failed to define broadcasted for $op")
+        end
+    end
+
+    # if it is a binary operation...
+    if length(args) == 2
+        args[1] isa Type || continue
+        args[2] isa Type || continue
+        Real <: args[1] && Real <: args[2] || continue
+
+        try
+            fn = function_lookup[op]
+            
+            @eval Base.broadcasted(x::$op, d::DualVector, y::Real) = broadcast_rule($fn, d, y)
+            @eval Base.broadcasted(x::$op, y::Real, d::DualVector) = broadcast_rule($fn, y, d)
+        catch
+            # skip over erroneous definitions. Uncomment the below for debugging.
+            # println("Failed to define broadcasted for $op")
         end
     end
 end
\ No newline at end of file
diff --git a/src/indexing.jl b/src/indexing.jl
index e2556ca..7036091 100644
--- a/src/indexing.jl
+++ b/src/indexing.jl
@@ -7,8 +7,6 @@ sparse_getindex(D::Diagonal, k::Integer, ::Colon) = OneElement(D.diag[k], k, siz
 sparse_getindex(D::Diagonal, ::Colon, j::Integer) = OneElement(D.diag[j], j, size(D, 1))
 
 """
-    getindex(x::DualVector, y::Int)
-
 Extract a single Dual number from a DualVector at position y.
 """
 function Base.getindex(x::DualVector, y::Int)
@@ -16,8 +14,6 @@ function Base.getindex(x::DualVector, y::Int)
 end
 
 """
-    getindex(x::DualVector, y::UnitRange)
-
 Extract a sub-DualVector from a DualVector using a range.
 """
 function Base.getindex(x::DualVector, y::UnitRange)
@@ -27,15 +23,11 @@ function Base.getindex(x::DualVector, y::UnitRange)
 end
 
 """
-    size(x::DualVector)
-
 Return the size of the DualVector (length of the value vector).
 """
 Base.size(x::DualVector) = (length(x.value),)
 
 """
-    axes(x::DualVector)
-
 Return the axes of the DualVector.
 """
 Base.axes(x::DualVector) = axes(x.value)
\ No newline at end of file
diff --git a/src/types.jl b/src/types.jl
index 09af0c8..a0a0c43 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -32,18 +32,18 @@ For now the entries just return the values when indexed.
 
 Constructs a DualVector, ensuring that the vector length matches the number of rows in the Jacobian.
 """
-struct DualVector{T, M <: AbstractMatrix{T}} <: AbstractVector{Dual{T}}
-    value::Vector{T}
+struct DualVector{T, V <: AbstractVector{T},M <: AbstractMatrix{T}} <: AbstractVector{Dual{T}}
+    value::V
     jacobian::M
-    
-    function DualVector(value::Vector{T}, jacobian::M) where {T, M <: AbstractMatrix{T}}
+
+    function DualVector(value::V, jacobian::M) where {T, V <: AbstractVector{T}, M <: AbstractMatrix{T}}
         if size(jacobian, 1) != length(value)
             x, y = length(value), size(jacobian, 1)
             throw(ArgumentError("vector length must match number of rows in jacobian.\n" *
                                "vector length: $x\n" *
                                "no. of jacobian rows: $y"))
         end
-        new{T, M}(value, jacobian)
+        new{T,V, M}(value, jacobian)
     end
 end
 
diff --git a/src/utilities.jl b/src/utilities.jl
index 5217b16..544a756 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -1,8 +1,6 @@
 # Utility functions for DualArrays.jl
 
 """
-    sum(x::DualVector)
-
 Sum all elements of a DualVector, returning a single Dual number.
 """
 function Base.sum(x::DualVector)
@@ -25,8 +23,6 @@ _size(x::Real) = 1
 _size(x::DualVector) = length(x.value)
 
 """
-    vcat(x::Union{Dual, DualVector}...)
-
 Vertically concatenate Dual numbers and DualVectors.
 """
 function Base.vcat(x::Union{Dual, DualVector}...)
@@ -39,8 +35,6 @@ function Base.vcat(x::Union{Dual, DualVector}...)
 end
 
 """
-    vcat(x::Union{Real, DualVector}...)
-
 Vertically concatenate Real numbers and DualVectors.
 """
 
@@ -56,8 +50,6 @@ function Base.vcat(x::Union{Real, DualVector}...)
 end
 
 """
-    show(io::IO, ::MIME"text/plain", x::DualVector)
-
 Custom display method for DualVectors.
 """
 Base.show(io::IO, ::MIME"text/plain", x::DualVector) = 

From 16f6a2f06902d22813bcbab99f267fc950418fde Mon Sep 17 00:00:00 2001
From: Max Vassiliev <vassiliev.max@gmail.com>
Date: Mon, 17 Nov 2025 14:28:24 +0000
Subject: [PATCH 3/7] fix bugs and add unit tests

---
 Project.toml           |  5 ++-
 src/arithmetic.jl      | 53 ++++++++++++++++++++++++++++--
 src/types.jl           |  4 +--
 src/utilities.jl       |  2 +-
 test/broadcast_test.jl | 74 ++++++++++++++++++++++++++++++++++++++++++
 test/runtests.jl       | 72 +++++++++++++++++++++++-----------------
 6 files changed, 172 insertions(+), 38 deletions(-)
 create mode 100644 test/broadcast_test.jl

diff --git a/Project.toml b/Project.toml
index 81c651a..e9efb10 100644
--- a/Project.toml
+++ b/Project.toml
@@ -12,6 +12,7 @@ FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 ChainRules = "1.72.6"
@@ -19,6 +20,7 @@ ChainRulesCore = "1.26.0"
 DifferentialEquations = "7.17.0"
 ForwardDiff = "1.2.2"
 Plots = "1.41.1"
+SparseArrays = "1.12.0"
 
 [extras]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
@@ -27,8 +29,9 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [targets]
 dev = ["Revise"]
 examples = ["Plots", "DifferentialEquations", "ForwardDiff", "BandedMatrices"]
-test = ["Test", "ForwardDiff", "BandedMatrices"]
+test = ["Test", "ForwardDiff", "BandedMatrices", "SparseArrays"]
diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index f84838c..f615e7c 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -44,7 +44,7 @@ function broadcast_rule(f, d::DualVector)
     jac = similar(d.jacobian)
 
     @inbounds for (i, x) in pairs(d.value)
-        y, dy = ChainRules.frule((ChainRules.NoTangent(), d.jacobian[i, :]), f, x[i])
+        y, dy = ChainRules.frule((ChainRules.NoTangent(), d.jacobian[i, :]), f, x)
         val[i] = y
         jac[i, :] = dy
     end
@@ -62,10 +62,11 @@ operation a (+) b on reals a, b should be broadcasted over:
 function broadcast_rule(f, d::DualVector, x::Real)
     val = similar(d.value)
     jac = similar(d.jacobian)
+    z = zero(jac[1, :])
 
     @inbounds for (i, y) in pairs(d.value)
         yval, dy = ChainRules.frule(
-            (ChainRules.NoTangent(), d.jacobian[i, :], zero(d.jacobian[i, :])),
+            (ChainRules.NoTangent(), d.jacobian[i, :], z),
              f, y, x)
         val[i] = yval
         jac[i, :] = dy
@@ -77,10 +78,11 @@ end
 function broadcast_rule(f, x::Real, d::DualVector)
     val = similar(d.value)
     jac = similar(d.jacobian)
+    z = zero(jac[1, :])
 
     @inbounds for (i, y) in pairs(d.value)
         yval, dy = ChainRules.frule(
-            (ChainRules.NoTangent(), zero(d.jacobian[i, :]), d.jacobian[i, :]),
+            (ChainRules.NoTangent(), z, d.jacobian[i, :]),
              f, x, y)
         val[i] = yval
         jac[i, :] = dy
@@ -89,9 +91,46 @@ function broadcast_rule(f, x::Real, d::DualVector)
     return DualVector(val, jac)
 end
 
+"""
+Extend for broadcasting Dual and DualVector
+"""
+
+function broadcast_rule(f, d::DualVector, x::Dual)
+    val = similar(d.value)
+    jac = similar(d.jacobian)
+    z = x.partials
+
+    @inbounds for (i, y) in pairs(d.value)
+        yval, dy = ChainRules.frule(
+            (ChainRules.NoTangent(), d.jacobian[i, :], z),
+             f, y, x.value)
+        val[i] = yval
+        jac[i, :] = dy
+    end
+
+    return DualVector(val, jac)
+end
+
+function broadcast_rule(f, x::Dual, d::DualVector)
+    val = similar(d.value)
+    jac = similar(d.jacobian)
+    z = x.partials
+
+    @inbounds for (i, y) in pairs(d.value)
+        yval, dy = ChainRules.frule(
+            (ChainRules.NoTangent(), z, d.jacobian[i, :]),
+             f, x.value, y)
+        val[i] = yval
+        jac[i, :] = dy
+    end
+
+    return DualVector(val, jac)
+end
+
 # Define a lookup table typeof(function) => function
 function_lookup = Dict{DataType, Function}()
 
+# todo: extend to LinearAlgebra and FastMath functions
 for n in names(Base)
     f = getfield(Base, n)
     if f isa Function
@@ -99,6 +138,9 @@ for n in names(Base)
     end
 end
 
+# a set of defined broadcasts to avoid duplicate definitions
+defined = Set{DataType}()
+
 # Get all applicable frules defined in ChainRules
 # and define broadcasted versions for DualVector using vector_rule
 
@@ -112,6 +154,7 @@ for f in frules
     op, args = sig.parameters[2], sig.parameters[3:end]
 
     isconcretetype(op) || continue
+    op in defined && continue
 
     # if it is a single argument function...
     if length(args) == 1
@@ -122,6 +165,7 @@ for f in frules
             fn = function_lookup[op]
 
             @eval Base.broadcasted(x::$op, d::DualVector) = broadcast_rule($fn, d)
+            push!(defined, op)
         catch e
             # skip over erroneous definitions. Uncomment the below for debugging.
             # println("Failed to define broadcasted for $op")
@@ -139,6 +183,9 @@ for f in frules
             
             @eval Base.broadcasted(x::$op, d::DualVector, y::Real) = broadcast_rule($fn, d, y)
             @eval Base.broadcasted(x::$op, y::Real, d::DualVector) = broadcast_rule($fn, y, d)
+            @eval Base.broadcasted(x::$op, d::DualVector, du::Dual) = broadcast_rule($fn, d, du)
+            @eval Base.broadcasted(x::$op, du::Dual, d::DualVector) = broadcast_rule($fn, du, d)
+            push!(defined, op)
         catch
             # skip over erroneous definitions. Uncomment the below for debugging.
             # println("Failed to define broadcasted for $op")
diff --git a/src/types.jl b/src/types.jl
index a0a0c43..144d713 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -48,9 +48,7 @@ struct DualVector{T, V <: AbstractVector{T},M <: AbstractMatrix{T}} <: AbstractV
 end
 
 """
-    DualVector(value::AbstractVector, jacobian::AbstractMatrix)
-
-Convenience constructor that promotes types to be compatible.
+Constructor that forces type compatibility
 """
 function DualVector(value::AbstractVector, jacobian::AbstractMatrix)
     T = promote_type(eltype(value), eltype(jacobian))
diff --git a/src/utilities.jl b/src/utilities.jl
index 544a756..44818c3 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -1,4 +1,4 @@
-# Utility functions for DualArrays.jl
+# Miscellaneous functions for DualArrays.jl
 
 """
 Sum all elements of a DualVector, returning a single Dual number.
diff --git a/test/broadcast_test.jl b/test/broadcast_test.jl
new file mode 100644
index 0000000..df0062b
--- /dev/null
+++ b/test/broadcast_test.jl
@@ -0,0 +1,74 @@
+using DualArrays, Test, SparseArrays, LinearAlgebra
+
+@testset "Broadcasting" begin
+    d = DualVector([1.0, 2.0, 3.0], [1 0 0; 0 1 0; 0 0 1])
+
+    # Test broadcasting with single argument functions
+    s = sin.(d)
+    c = cos.(d)
+    e = exp.(d)
+    l = log.(d)
+    a = tanh.(d)
+
+    @test s isa DualVector
+    @test c isa DualVector
+    @test e isa DualVector
+    @test l isa DualVector
+    @test a isa DualVector
+
+    @test s.value ≈ sin.([1.0, 2.0, 3.0])
+    @test c.value ≈ cos.([1.0, 2.0, 3.0])
+    @test e.value ≈ exp.([1.0, 2.0, 3.0])
+    @test l.value ≈ log.([1.0, 2.0, 3.0])
+    @test a.value ≈ tanh.([1.0, 2.0, 3.0])
+
+    @test s.jacobian ≈ Diagonal(cos.([1.0, 2.0, 3.0]))
+    @test c.jacobian ≈ Diagonal(-sin.([1.0, 2.0, 3.0]))
+    @test e.jacobian ≈ Diagonal(exp.([1.0, 2.0, 3.0]))
+    @test l.jacobian ≈ Diagonal(1.0 ./ [1.0, 2.0, 3.0])
+    @test a.jacobian ≈ Diagonal(1.0 .- tanh.([1.0, 2.0, 3.0]).^2)
+
+    # Test broadcasting with binary operations
+    a = d .+ 2.0
+    s = 3.0 .- d
+    m = d .* 4.0
+    div = 8.0 ./ d
+
+    @test a isa DualVector
+    @test s isa DualVector 
+    @test m isa DualVector
+    @test div isa DualVector
+
+    @test a.value ≈ [3.0, 4.0, 5.0]
+    @test s.value ≈ [2.0, 1.0, 0.0]
+    @test m.value ≈ [4.0, 8.0, 12.0]
+    @test div.value ≈ [8.0, 4.0, 8/3.0]
+
+    @test a.jacobian ≈ d.jacobian
+    @test s.jacobian ≈ -d.jacobian
+    @test m.jacobian ≈ 4.0 * d.jacobian
+    @test div.jacobian ≈ -8.0 ./ (d.value .^ 2) .* d.jacobian
+
+    # Test broadcasting with binary operations on a Dual and DualVector
+    x = Dual(2.0, [1.0, 2.0, 3.0])
+
+    a = x .+ d
+    s = d .- x
+    m = x .* d
+    div = d ./ x
+
+    @test a isa DualVector
+    @test s isa DualVector
+    @test m isa DualVector
+    @test div isa DualVector
+
+    @test a.value ≈ [3.0, 4.0, 5.0]
+    @test s.value ≈ [-1.0, 0.0, 1.0]
+    @test m.value ≈ [2.0, 4.0, 6.0]
+    @test div.value ≈ [0.5, 1.0, 1.5]
+
+    @test a.jacobian ≈ [2.0 2.0 3.0; 1.0 3.0 3.0; 1.0 2.0 4.0]
+    @test s.jacobian ≈ [0.0 -2.0 -3.0; -1.0 -1.0 -3.0; -1.0 -2.0 -2.0]
+    @test m.jacobian ≈ [3.0 2.0 3.0; 2.0 6.0 6.0; 3.0 6.0 11.0]
+    @test div.jacobian ≈ [0.25 -0.5 -0.75; -0.5 -0.5 -1.5; -0.75 -1.5 -1.75]
+end
\ No newline at end of file
diff --git a/test/runtests.jl b/test/runtests.jl
index d1cba3c..12a6fc9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2,41 +2,53 @@ using DualArrays, Test, LinearAlgebra, ForwardDiff, BandedMatrices
 using DualArrays: Dual
 
 @testset "DualArrays" begin
-    @test_throws ArgumentError DualVector([1,2],I(3))
     
-    v = DualVector([1,2, 3], [1 2 3; 4 5 6;7 8 9])
-    @test v[1] isa Dual
-    @test v[1] == Dual(1,[1,2,3])
-    @test v[2] == Dual(2,[4,5,6])
-    @test v[3] == Dual(3,[7,8,9])
-    @test_throws BoundsError v[4]
-    @test v == DualVector([1,2, 3], [1 2 3; 4 5 6;7 8 9])
+    @testset "Type Definition" begin
+        @test_throws ArgumentError DualVector([1,2],I(3))
+    end
+    
+    @testset "Indexing" begin
+        v = DualVector([1, 2, 3], [1 2 3; 4 5 6;7 8 9])
+        @test v[1] isa Dual
+        @test v[1] == Dual(1,[1,2,3])
+        @test v[2] == Dual(2,[4,5,6])
+        @test v[3] == Dual(3,[7,8,9])
+        @test_throws BoundsError v[4]
+        @test v == DualVector([1,2, 3], [1 2 3; 4 5 6;7 8 9])
 
-    w = v + v
-    @test w == DualVector([2,4,6],[2 4 6;8 10 12;14 16 18])
-    @test w.jacobian == 2v.jacobian
+        x,y = v[1:2],v[2:3]
+        @test x == DualVector([1,2],[1 2 3;4 5 6])
+        @test y == DualVector([2,3],[4 5 6;7 8 9])
 
-    @test sin.(v) isa DualVector
-    @test sin.(v).value == [sin(1), sin(2), sin(3)]
-    @test sin.(v).jacobian == Diagonal(cos.(v.value)) * v.jacobian
+        n = 10
+        v = DualVector(1:n, I(n))
+        @test v[2:end].jacobian isa BandedMatrix
 
-    x,y = v[1:2],v[2:3]
-    @test x == DualVector([1,2],[1 2 3;4 5 6])
-    @test y == DualVector([2,3],[4 5 6;7 8 9])
-    @test x .* y == DualVector([2,6],[6 9 12;26 31 36])
+        @test sum(v[1:end-1] .* v[2:end]).partials == ForwardDiff.gradient(v -> sum(v[1:end-1] .* v[2:end]), 1:n)
+    end
     
-    @test sum(x .* y) isa Dual
-    @test sum(x .* y) == Dual(8,[32,40,48])
-
-    n = 10
-    v = DualVector(1:n, I(n))
-    @test v[2:end].jacobian isa BandedMatrix
+    @testset "Arithmetic" begin
+        v = DualVector([1, 2, 3], [1 2 3; 4 5 6;7 8 9])
+        w = v + v
+        @test w == DualVector([2,4,6],[2 4 6;8 10 12;14 16 18])
+        @test w.jacobian == 2v.jacobian
 
-    @test sum(v[1:end-1] .* v[2:end]).partials == ForwardDiff.gradient(v -> sum(v[1:end-1] .* v[2:end]), 1:n)
+        x = Dual(1, [1, 2, 3])
+        y = DualVector([2, 3], [4 5 6;7 8 9])
 
-    x = Dual(1, [1, 2, 3])
-    y = DualVector([2, 3], [4 5 6;7 8 9])
-    @test vcat(x) == x
-    @test vcat(x, x) == DualVector([1, 1], [1 2 3;1 2 3])
-    @test vcat(x, y) == DualVector([1, 2, 3], [1 2 3;4 5 6;7 8 9])
+        @test x .* y == DualVector([2,3],[6 9 12;10 14 18])
+        
+        @test sum(x .* y) isa Dual
+        @test sum(x .* y) == Dual(5,[16,23,30])
+    end
+    
+    @testset "vcat" begin
+        x = Dual(1, [1, 2, 3])
+        y = DualVector([2, 3], [4 5 6;7 8 9])
+        @test vcat(x) == x
+        @test vcat(x, x) == DualVector([1, 1], [1 2 3;1 2 3])
+        @test vcat(x, y) == DualVector([1, 2, 3], [1 2 3;4 5 6;7 8 9])
+    end
+    
+    include("broadcast_test.jl")
 end
\ No newline at end of file

From e5a3b10467b8ba0fbe8b33042af46f81c7d07ac0 Mon Sep 17 00:00:00 2001
From: Max Vassiliev <vassiliev.max@gmail.com>
Date: Mon, 17 Nov 2025 14:36:02 +0000
Subject: [PATCH 4/7] fix compatibility issue

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index e9efb10..a2b5947 100644
--- a/Project.toml
+++ b/Project.toml
@@ -20,7 +20,7 @@ ChainRulesCore = "1.26.0"
 DifferentialEquations = "7.17.0"
 ForwardDiff = "1.2.2"
 Plots = "1.41.1"
-SparseArrays = "1.12.0"
+SparseArrays = "1.10"
 
 [extras]
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"

From 7ca2cc5226beb8341cffd4edb8377e72d51c8c74 Mon Sep 17 00:00:00 2001
From: Maxim Vassiliev <76599693+max-vassili3v@users.noreply.github.com>
Date: Wed, 19 Nov 2025 10:18:04 +0000
Subject: [PATCH 5/7] correct version

Co-authored-by: Sheehan Olver <solver@mac.com>
---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index a2b5947..4a4e901 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "DualArrays"
 uuid = "429e4e16-f749-45f3-beec-30742fae38ce"
-version = "0.1.1"
+version = "0.2.0"
 authors = ["Sheehan Olver <solver@mac.com>"]
 
 [deps]

From 16d3482ac321cb3acb3bdbc37484b899f51aae21 Mon Sep 17 00:00:00 2001
From: Maxim Vassiliev <vassiliev.max@gmail.com>
Date: Wed, 19 Nov 2025 10:50:54 +0000
Subject: [PATCH 6/7] update readme

---
 README.md | 32 ++++++--------------------------
 1 file changed, 6 insertions(+), 26 deletions(-)

diff --git a/README.md b/README.md
index a3976fc..078a9f9 100644
--- a/README.md
+++ b/README.md
@@ -1,31 +1,9 @@
 # DualArrays.jl
-DualArrays.jl is a package that provides the `DualVector` structure for use in forward mode automatic differentiation (autodiff). Existing forward mode autodiff implementations such as `ForwardDiff.jl` make use of a `Dual` structure representing Dual numbers as follows:
+DualArrays.jl is a package that provides the `DualVector` structure for use in forward mode automatic differentiation (autodiff). Existing forward mode autodiff implementations such as `ForwardDiff.jl` make use of a `Dual` structure with a vector of dual parts.
 
-```julia
-struct Dual{T,V<:Real,N} <: Real
-    value::V
-    partials::Partials{N,V}
-end
-```
-
-where partials is defined as follows:
-
-```julia
-struct Partials{N,V} <: AbstractVector{V}
-    values::NTuple{N,V}
-end
-```
+There are some limitations of this when differentiating vector valued functions, as the dual components of each element will be each treated as a separate (dense) vector rather than a Jacobian. This misses the opportunity to exploit sparse Jacobian structures such as banded structures that appear when solving systems of ODEs. 
 
-There are some limitations of this when differentiating vector valued functions, as having a vector of `Dual`s cannot properly exploit sparse Jacobian structures such as banded structures that appear when solving systems of ODEs. `DualArrays.jl` introduces a new `DualVector` type as follows:
-
-```julia
-struct DualVector{T, M <: AbstractMatrix{T}} <: AbstractVector{Dual{T}}
-    value::Vector{T}
-    jacobian::M
-end
-```
-
-with a similarly defined `Dual` type. This allows for jacobians to be any sparse matrix structure defined in the Julia linear algebra ecosystem. An efficient implementation of differentiating a function might be as follows:
+`DualArrays.jl` provides a new structure, `DualVector`, consisting of a vector of real parts and a jacobian that can be any matrix structure in the Julia ecosystem. This carries over many of the optimisations provided by sparse matrix structures to the forward-mode autodiff process. The package also comes with its own `Dual` type for elementwise indexing or vector -> scalar functions. An efficient implementation of differentiating a function might be as follows:
 
 ```julia
 using FillArrays
@@ -36,4 +14,6 @@ function gradient(f::Function, x::Vector)
 end
 ```
 
-See the examples folder for more use cases
+See the examples folder for more use cases.
+
+Differentiation rules are mostly provided by the `ChainRules.jl` autodiff backend.

From 1cc20def40ef73915267c6aa566ee1d8ed0cc424 Mon Sep 17 00:00:00 2001
From: Maxim Vassiliev <vassiliev.max@gmail.com>
Date: Wed, 19 Nov 2025 13:14:39 +0000
Subject: [PATCH 7/7] remove unnecessary lookup table

---
 src/arithmetic.jl | 40 ++++++++--------------------------------
 1 file changed, 8 insertions(+), 32 deletions(-)

diff --git a/src/arithmetic.jl b/src/arithmetic.jl
index f615e7c..1fb94f7 100644
--- a/src/arithmetic.jl
+++ b/src/arithmetic.jl
@@ -127,17 +127,6 @@ function broadcast_rule(f, x::Dual, d::DualVector)
     return DualVector(val, jac)
 end
 
-# Define a lookup table typeof(function) => function
-function_lookup = Dict{DataType, Function}()
-
-# todo: extend to LinearAlgebra and FastMath functions
-for n in names(Base)
-    f = getfield(Base, n)
-    if f isa Function
-        function_lookup[typeof(f)] = f
-    end
-end
-
 # a set of defined broadcasts to avoid duplicate definitions
 defined = Set{DataType}()
 
@@ -161,15 +150,8 @@ for f in frules
         args[1] isa Type || continue
         Real <: args[1] || continue
 
-        try
-            fn = function_lookup[op]
-
-            @eval Base.broadcasted(x::$op, d::DualVector) = broadcast_rule($fn, d)
-            push!(defined, op)
-        catch e
-            # skip over erroneous definitions. Uncomment the below for debugging.
-            # println("Failed to define broadcasted for $op")
-        end
+        @eval Base.broadcasted(fn::$op, d::DualVector) = broadcast_rule(fn, d)
+        push!(defined, op)
     end
 
     # if it is a binary operation...
@@ -178,17 +160,11 @@ for f in frules
         args[2] isa Type || continue
         Real <: args[1] && Real <: args[2] || continue
 
-        try
-            fn = function_lookup[op]
-            
-            @eval Base.broadcasted(x::$op, d::DualVector, y::Real) = broadcast_rule($fn, d, y)
-            @eval Base.broadcasted(x::$op, y::Real, d::DualVector) = broadcast_rule($fn, y, d)
-            @eval Base.broadcasted(x::$op, d::DualVector, du::Dual) = broadcast_rule($fn, d, du)
-            @eval Base.broadcasted(x::$op, du::Dual, d::DualVector) = broadcast_rule($fn, du, d)
-            push!(defined, op)
-        catch
-            # skip over erroneous definitions. Uncomment the below for debugging.
-            # println("Failed to define broadcasted for $op")
-        end
+        
+        @eval Base.broadcasted(fn::$op, d::DualVector, y::Real) = broadcast_rule(fn, d, y)
+        @eval Base.broadcasted(fn::$op, y::Real, d::DualVector) = broadcast_rule(fn, y, d)
+        @eval Base.broadcasted(fn::$op, d::DualVector, du::Dual) = broadcast_rule(fn, d, du)
+        @eval Base.broadcasted(fn::$op, du::Dual, d::DualVector) = broadcast_rule(fn, du, d)
+        push!(defined, op)
     end
 end
\ No newline at end of file