From bccb4f15e9e2cfc35225bc9ddaf47809f88a1b2d Mon Sep 17 00:00:00 2001
From: Alex Jones <alex.lemony@gmail.com>
Date: Tue, 1 Aug 2023 19:19:24 +0100
Subject: [PATCH 01/12] Add ArrayDifferentialOperators

---
 src/diff.jl | 120 ++++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 112 insertions(+), 8 deletions(-)

diff --git a/src/diff.jl b/src/diff.jl
index 11c523dab..c6eab9bea 100644
--- a/src/diff.jl
+++ b/src/diff.jl
@@ -1,4 +1,5 @@
-abstract type Operator <: Function end
+abstract type AbstractOperator <: Function end
+abstract type Operator <: AbstractOperator end
 
 """
 $(TYPEDEF)
@@ -33,17 +34,17 @@ struct Differential <: Operator
     x
     Differential(x) = new(value(x))
 end
-(D::Differential)(x) = Term{symtype(x)}(D, [x])
-(D::Differential)(x::Num) = Num(D(value(x)))
-(D::Differential)(x::Complex{Num}) = wrap(ComplexTerm{Real}(D(unwrap(real(x))), D(unwrap(imag(x)))))
+(D::Operator)(x) = Term{symtype(x)}(D, [x])
+(D::Operator)(x::Num) = Num(D(value(x)))
+(D::Operator)(x::Complex{Num}) = wrap(ComplexTerm{Real}(D(unwrap(real(x))), D(unwrap(imag(x)))))
 SymbolicUtils.promote_symtype(::Differential, x) = x
 
 is_derivative(x) = istree(x) ? operation(x) isa Differential : false
 
-Base.:*(D1, D2::Differential) = D1 ∘ D2
-Base.:*(D1::Differential, D2) = D1 ∘ D2
-Base.:*(D1::Differential, D2::Differential) = D1 ∘ D2
-Base.:^(D::Differential, n::Integer) = _repeat_apply(D, n)
+Base.:*(D1, D2::Operator) = D1 ∘ D2
+Base.:*(D1::Operator, D2) = D1 ∘ D2
+Base.:*(D1::Operator, D2::Operator) = D1 ∘ D2
+Base.:^(D::Operator, n::Integer) = _repeat_apply(D, n)
 
 Base.show(io::IO, D::Differential) = print(io, "Differential(", D.x, ")")
 
@@ -785,3 +786,106 @@ end
 function SymbolicUtils.substitute(op::Differential, dict; kwargs...)
     @set! op.x = substitute(op.x, dict; kwargs...)
 end
+
+#######################################################################################################################
+# Vector Calculus
+#######################################################################################################################
+abstract type ArrayOperator <: AbstractOperator end
+
+struct ArrayDifferentialOperator <: ArrayOperator
+    """The variables to differentiate with respect to."""
+    vars
+    differentials
+    name
+    ArrayDifferentialOperator(vars, differentials, name) = new(ArrayOp(vars), ArrayOp(differentials), name)
+    ArrayDifferentialOperator(vars::ArrayOp, diffs, name) = new(vars, ArrayOp(diffs), name)
+end
+Nabla(vars) = ArrayDifferentialOperator(ArrayOp(value.(vars)), map(Differential, value.(vars)), "∇")
+Div(vars) = (x) -> Nabla(vars) ⋅ x
+Curl(vars) = (x) -> Nabla(vars) × x
+Laplacian(vars) = Nabla(vars) ⋅ Nabla(vars)
+
+#? How to get transpose and Jac working?
+
+function (D::ArrayDifferentialOperator)(x::SymVec)
+    @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
+    _call(d, x) = d(x)
+    map(_call, zip(D.differentials, x))
+end
+(D::ArrayDifferentialOperator)(x::Arr) = Arr(D(value(x)))
+
+function (D1::ArrayDifferentialOperator)(D2::ArrayDifferentialOperator)
+    @assert all(x -> any(isequal.((x,), D2.vars)), D1.vars)
+
+    ArrayDifferentialOperator(D1.vars, scalarize(D1.differentials .∘ D2.differentials), "("*D1.name*"∘"*D2.name*")")
+end    
+
+function LinearAlgebra.dot(D::ArrayDifferentialOperator, x::SymVec)
+    @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
+    _call(d, x) = d(x)
+    sum(_call, zip(D.differentials, x))
+end 
+LinearAlgebra.dot(D::ArrayDifferentialOperator, x::Arr) = Arr(D ⋅ value(x))
+
+function LinearAlgebra.dot(x::SymVec, D::ArrayDifferentialOperator)
+    @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
+    (x) -> sum(D -> D(x), D.differentials)
+end
+LinearAlgebra.dot(x::Arr, D::ArrayDifferentialOperator) = value(x) ⋅ D
+
+function LinearAlgebra.dot(D1::ArrayDifferentialOperator, D2::ArrayDifferentialOperator)
+    @assert all(scalarize(isequal.(D1.vars, D2.vars))) "Operators have different variables and cannot be composed."
+    (x) -> sum(i -> (D1.differentials[i] ∘ D2.differentials[i])(x), eachindex(D1.vars))
+end
+
+function εijk_cond(i, j, k)
+    if (i, j, k) in [(1, 2, 3), (2, 3, 1), (3, 1, 2)]
+        1
+    elseif (i == j) || (j == k) || (k == i)
+        0
+    else 
+        -1
+    end
+end
+
+function crosscompose(a, b)
+    v1 = x -> (a[2] ∘ b[3])(x) - (a[3] ∘ b[2])(x)
+    v2 = x -> (a[3] ∘ b[1])(x) - (a[1] ∘ b[3])(x)
+    v3 = x -> (a[1] ∘ b[2])(x) - (a[2] ∘ b[1])(x) 
+    return [v1, v2, v3]
+end 
+
+function LinearAlgebra.cross(D::ArrayDifferentialOperator, x::SymVec)
+    @assert length(D.vars) == length(x) == 3 "Cross product is only defined in 3 dimensions."
+    ε = [εijk_cond(i, j, k) for i in 1:3, j in 1:3, k in 1:3]
+    curl(i) = sum(j -> sum(k -> expand_derivatives(ε[i, j, k]*D.differentials[j](x[k])), 1:3), 1:3)
+    
+    return map(curl, ArrayOp(1:3))
+end
+LinearAlgebra.cross(D::ArrayDifferentialOperator, x::Arr) = Arr(D × x)
+
+function LinearAlgebra.cross(D1::ArrayDifferentialOperator, D2::ArrayDifferentialOperator)
+    @assert length(D1.vars) == length(D2.vars) == 3 "Cross product is only defined in 3 dimensions."
+    @assert all(scalarize(isequal.(D1.vars, D2.vars))) "Operators have different variables and cannot be composed."
+
+   ArrayDifferentialOperator(D1.vars, crosscompose(D1.differentials, D2.differentials), "("*D1.name*"×"*D2.name*")")
+end
+
+SymbolicUtils.promote_symtype(::Nabla, x) = x
+
+is_derivative(x) = istree(x) ? operation(x) isa Differential : false
+
+Base.:*(D1, D2::Differential) = D1 ∘ D2
+Base.:*(D1::Differential, D2) = D1 ∘ D2
+Base.:*(D1::Differential, D2::Differential) = D1 ∘ D2
+Base.:^(D::Differential, n::Integer) = _repeat_apply(D, n)
+
+Base.show(io::IO, D::ArrayDifferentialOperator) = print(io, "(D.name)(", scalarize(D.vars), ")")
+
+function Base.:(==)(D1::ArrayDifferentialOperator, D2::ArrayDifferentialOperator) 
+    @variables x[1:length(D1.vars)]
+    all(scalarize(isequal.(D1.vars, D2.vars))) && all(scalarize(isequal.(D1(x), D2(x))))
+
+
+_isfalse(occ::Bool) = occ === false
+_isfalse(occ::Symbolic) = istree(occ) && _isfalse(operation(occ))

From 03ac562c3c8ca79abadbd59407639bc1c29ededb Mon Sep 17 00:00:00 2001
From: Alex Jones <alex.lemony@gmail.com>
Date: Tue, 1 Aug 2023 19:23:27 +0100
Subject: [PATCH 02/12] tidy up

---
 src/diff.jl | 12 ++----------
 1 file changed, 2 insertions(+), 10 deletions(-)

diff --git a/src/diff.jl b/src/diff.jl
index c6eab9bea..3996de897 100644
--- a/src/diff.jl
+++ b/src/diff.jl
@@ -873,19 +873,11 @@ end
 
 SymbolicUtils.promote_symtype(::Nabla, x) = x
 
-is_derivative(x) = istree(x) ? operation(x) isa Differential : false
-
-Base.:*(D1, D2::Differential) = D1 ∘ D2
-Base.:*(D1::Differential, D2) = D1 ∘ D2
-Base.:*(D1::Differential, D2::Differential) = D1 ∘ D2
-Base.:^(D::Differential, n::Integer) = _repeat_apply(D, n)
+is_derivative(x) = istree(x) ? operation(x) isa ArrayDifferentialOperator : false
 
 Base.show(io::IO, D::ArrayDifferentialOperator) = print(io, "(D.name)(", scalarize(D.vars), ")")
 
 function Base.:(==)(D1::ArrayDifferentialOperator, D2::ArrayDifferentialOperator) 
     @variables x[1:length(D1.vars)]
     all(scalarize(isequal.(D1.vars, D2.vars))) && all(scalarize(isequal.(D1(x), D2(x))))
-
-
-_isfalse(occ::Bool) = occ === false
-_isfalse(occ::Symbolic) = istree(occ) && _isfalse(operation(occ))
+end
\ No newline at end of file

From 69c126570e0bd531351567343a00b95b1166deb8 Mon Sep 17 00:00:00 2001
From: Alex Jones <alex.lemony@gmail.com>
Date: Tue, 1 Aug 2023 19:28:08 +0100
Subject: [PATCH 03/12] a fix

---
 src/diff.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/diff.jl b/src/diff.jl
index 3996de897..b03ac0b12 100644
--- a/src/diff.jl
+++ b/src/diff.jl
@@ -822,14 +822,13 @@ end
 
 function LinearAlgebra.dot(D::ArrayDifferentialOperator, x::SymVec)
     @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
-    _call(d, x) = d(x)
     sum(_call, zip(D.differentials, x))
 end 
 LinearAlgebra.dot(D::ArrayDifferentialOperator, x::Arr) = Arr(D ⋅ value(x))
 
 function LinearAlgebra.dot(x::SymVec, D::ArrayDifferentialOperator)
     @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
-    (x) -> sum(D -> D(x), D.differentials)
+    (x) -> sum((X, D) -> X*D(x), zip(x, D.differentials))
 end
 LinearAlgebra.dot(x::Arr, D::ArrayDifferentialOperator) = value(x) ⋅ D
 

From 72b78553fa77f7c1e747b6e8c70203430cdf8cdc Mon Sep 17 00:00:00 2001
From: Alex Jones <alex.lemony@gmail.com>
Date: Tue, 1 Aug 2023 19:29:17 +0100
Subject: [PATCH 04/12] another fix

---
 src/diff.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/diff.jl b/src/diff.jl
index b03ac0b12..6c2fbea95 100644
--- a/src/diff.jl
+++ b/src/diff.jl
@@ -822,13 +822,14 @@ end
 
 function LinearAlgebra.dot(D::ArrayDifferentialOperator, x::SymVec)
     @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
+    _call(d, x) = d(x)
     sum(_call, zip(D.differentials, x))
 end 
 LinearAlgebra.dot(D::ArrayDifferentialOperator, x::Arr) = Arr(D ⋅ value(x))
 
 function LinearAlgebra.dot(x::SymVec, D::ArrayDifferentialOperator)
     @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
-    (x) -> sum((X, D) -> X*D(x), zip(x, D.differentials))
+    (y) -> sum((X, D) -> X*D(y), zip(x, D.differentials))
 end
 LinearAlgebra.dot(x::Arr, D::ArrayDifferentialOperator) = value(x) ⋅ D
 

From dfc7d94f154a8f1199316b1120a3ad5c2d2245ca Mon Sep 17 00:00:00 2001
From: Alex Jones <alex.lemony@gmail.com>
Date: Fri, 4 Aug 2023 18:17:43 +0100
Subject: [PATCH 05/12] (almost) working

---
 src/arrays.jl |  1 +
 src/diff.jl   | 44 ++++++++++++++++++++++++--------------------
 2 files changed, 25 insertions(+), 20 deletions(-)

diff --git a/src/arrays.jl b/src/arrays.jl
index 218545b97..19aab684d 100644
--- a/src/arrays.jl
+++ b/src/arrays.jl
@@ -251,6 +251,7 @@ end
 # turn `f(x...)` into `term(f, x...)`
 #
 function call2term(expr, arrs=[])
+    (expr isa QuoteNode) && return expr
     !(expr isa Expr) && return :($unwrap($expr))
     if expr.head == :call
         if expr.args[1] == :(:)
diff --git a/src/diff.jl b/src/diff.jl
index 6c2fbea95..5f354b03b 100644
--- a/src/diff.jl
+++ b/src/diff.jl
@@ -787,20 +787,20 @@ function SymbolicUtils.substitute(op::Differential, dict; kwargs...)
     @set! op.x = substitute(op.x, dict; kwargs...)
 end
 
+
 #######################################################################################################################
 # Vector Calculus
 #######################################################################################################################
-abstract type ArrayOperator <: AbstractOperator end
+abstract type ArrayOperator end
 
 struct ArrayDifferentialOperator <: ArrayOperator
-    """The variables to differentiate with respect to."""
+    """The variables to differentiate with resp≈ect to."""
     vars
     differentials
     name
-    ArrayDifferentialOperator(vars, differentials, name) = new(ArrayOp(vars), ArrayOp(differentials), name)
-    ArrayDifferentialOperator(vars::ArrayOp, diffs, name) = new(vars, ArrayOp(diffs), name)
+    ArrayDifferentialOperator(vars, differentials, name) = new(vars, differentials, name)
 end
-Nabla(vars) = ArrayDifferentialOperator(ArrayOp(value.(vars)), map(Differential, value.(vars)), "∇")
+Nabla(vars) = ArrayDifferentialOperator(value.(vars), map(Differential, scalarize(value.(vars))), "∇")
 Div(vars) = (x) -> Nabla(vars) ⋅ x
 Curl(vars) = (x) -> Nabla(vars) × x
 Laplacian(vars) = Nabla(vars) ⋅ Nabla(vars)
@@ -809,8 +809,7 @@ Laplacian(vars) = Nabla(vars) ⋅ Nabla(vars)
 
 function (D::ArrayDifferentialOperator)(x::SymVec)
     @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
-    _call(d, x) = d(x)
-    map(_call, zip(D.differentials, x))
+    @arrayop (i,) unwrap((D.differentials)[i](x[i])) term=D(y) reduce=+
 end
 (D::ArrayDifferentialOperator)(x::Arr) = Arr(D(value(x)))
 
@@ -822,20 +821,20 @@ end
 
 function LinearAlgebra.dot(D::ArrayDifferentialOperator, x::SymVec)
     @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
-    _call(d, x) = d(x)
-    sum(_call, zip(D.differentials, x))
+    sum(D(x))
 end 
-LinearAlgebra.dot(D::ArrayDifferentialOperator, x::Arr) = Arr(D ⋅ value(x))
+LinearAlgebra.dot(D::ArrayDifferentialOperator, x::Arr) = D ⋅ value(x)
 
 function LinearAlgebra.dot(x::SymVec, D::ArrayDifferentialOperator)
     @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
-    (y) -> sum((X, D) -> X*D(y), zip(x, D.differentials))
+    (y) -> sum(@arrayop (i,) x[i]*D.differentials[i](y) term = (x⋅D)(y))
 end
 LinearAlgebra.dot(x::Arr, D::ArrayDifferentialOperator) = value(x) ⋅ D
 
 function LinearAlgebra.dot(D1::ArrayDifferentialOperator, D2::ArrayDifferentialOperator)
     @assert all(scalarize(isequal.(D1.vars, D2.vars))) "Operators have different variables and cannot be composed."
-    (x) -> sum(i -> (D1.differentials[i] ∘ D2.differentials[i])(x), eachindex(D1.vars))
+    lap = x -> sum((D1.differentials[i] ∘ D2.differentials[i])(x) for i in 1:length(D1.vars))
+    (x) -> @arrayop (i,) lap(x[i]) term=(D1⋅D2)(x) reduce=+
 end
 
 function εijk_cond(i, j, k)
@@ -855,14 +854,18 @@ function crosscompose(a, b)
     return [v1, v2, v3]
 end 
 
+function crosscall(a, b)
+    v1 = a[2](b[3]) - a[3](b[2])
+    v2 = a[3](b[1]) - a[1](b[3])
+    v3 = a[1](b[2]) - a[2](b[1])
+    return [v1, v2, v3]
+end
 function LinearAlgebra.cross(D::ArrayDifferentialOperator, x::SymVec)
     @assert length(D.vars) == length(x) == 3 "Cross product is only defined in 3 dimensions."
-    ε = [εijk_cond(i, j, k) for i in 1:3, j in 1:3, k in 1:3]
-    curl(i) = sum(j -> sum(k -> expand_derivatives(ε[i, j, k]*D.differentials[j](x[k])), 1:3), 1:3)
-    
-    return map(curl, ArrayOp(1:3))
+    curl = crosscall(D.differentials, x)
+    @arrayop (i,) curl[i] term=D×y 
 end
-LinearAlgebra.cross(D::ArrayDifferentialOperator, x::Arr) = Arr(D × x)
+LinearAlgebra.cross(D::ArrayDifferentialOperator, x::Arr) = Arr(D × value(x))
 
 function LinearAlgebra.cross(D1::ArrayDifferentialOperator, D2::ArrayDifferentialOperator)
     @assert length(D1.vars) == length(D2.vars) == 3 "Cross product is only defined in 3 dimensions."
@@ -871,13 +874,14 @@ function LinearAlgebra.cross(D1::ArrayDifferentialOperator, D2::ArrayDifferentia
    ArrayDifferentialOperator(D1.vars, crosscompose(D1.differentials, D2.differentials), "("*D1.name*"×"*D2.name*")")
 end
 
-SymbolicUtils.promote_symtype(::Nabla, x) = x
+SymbolicUtils.promote_symtype(::ArrayDifferentialOperator, x) = x
 
 is_derivative(x) = istree(x) ? operation(x) isa ArrayDifferentialOperator : false
 
-Base.show(io::IO, D::ArrayDifferentialOperator) = print(io, "(D.name)(", scalarize(D.vars), ")")
+Base.show(io::IO, D::ArrayDifferentialOperator) = print(io, D.name)
+Base.nameof(D::ArrayDifferentialOperator) = Symbol(D.name)
 
 function Base.:(==)(D1::ArrayDifferentialOperator, D2::ArrayDifferentialOperator) 
     @variables x[1:length(D1.vars)]
     all(scalarize(isequal.(D1.vars, D2.vars))) && all(scalarize(isequal.(D1(x), D2(x))))
-end
\ No newline at end of file
+end

From 73352b95595225daec9c73bcf5e89a5418e6545a Mon Sep 17 00:00:00 2001
From: Alex Jones <alex.lemony@gmail.com>
Date: Mon, 7 Aug 2023 17:10:21 +0100
Subject: [PATCH 06/12] final fixes and tests

---
 src/diff.jl  | 19 +++++--------------
 test/diff.jl | 38 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 43 insertions(+), 14 deletions(-)

diff --git a/src/diff.jl b/src/diff.jl
index 5f354b03b..131a1fe0d 100644
--- a/src/diff.jl
+++ b/src/diff.jl
@@ -809,7 +809,7 @@ Laplacian(vars) = Nabla(vars) ⋅ Nabla(vars)
 
 function (D::ArrayDifferentialOperator)(x::SymVec)
     @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
-    @arrayop (i,) unwrap((D.differentials)[i](x[i])) term=D(y) reduce=+
+    @arrayop (i,) (D.differentials)[i](x[i]) term=D(y) 
 end
 (D::ArrayDifferentialOperator)(x::Arr) = Arr(D(value(x)))
 
@@ -821,9 +821,10 @@ end
 
 function LinearAlgebra.dot(D::ArrayDifferentialOperator, x::SymVec)
     @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
-    sum(D(x))
+    @show D(x), scalarize(D(x))
+    sum(scalarize(D(x)))
 end 
-LinearAlgebra.dot(D::ArrayDifferentialOperator, x::Arr) = D ⋅ value(x)
+LinearAlgebra.dot(D::ArrayDifferentialOperator, x::Arr) = Num(D ⋅ value(x))
 
 function LinearAlgebra.dot(x::SymVec, D::ArrayDifferentialOperator)
     @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
@@ -837,16 +838,6 @@ function LinearAlgebra.dot(D1::ArrayDifferentialOperator, D2::ArrayDifferentialO
     (x) -> @arrayop (i,) lap(x[i]) term=(D1⋅D2)(x) reduce=+
 end
 
-function εijk_cond(i, j, k)
-    if (i, j, k) in [(1, 2, 3), (2, 3, 1), (3, 1, 2)]
-        1
-    elseif (i == j) || (j == k) || (k == i)
-        0
-    else 
-        -1
-    end
-end
-
 function crosscompose(a, b)
     v1 = x -> (a[2] ∘ b[3])(x) - (a[3] ∘ b[2])(x)
     v2 = x -> (a[3] ∘ b[1])(x) - (a[1] ∘ b[3])(x)
@@ -884,4 +875,4 @@ Base.nameof(D::ArrayDifferentialOperator) = Symbol(D.name)
 function Base.:(==)(D1::ArrayDifferentialOperator, D2::ArrayDifferentialOperator) 
     @variables x[1:length(D1.vars)]
     all(scalarize(isequal.(D1.vars, D2.vars))) && all(scalarize(isequal.(D1(x), D2(x))))
-end
+end
\ No newline at end of file
diff --git a/test/diff.jl b/test/diff.jl
index 7486a78e7..17dff5bc7 100644
--- a/test/diff.jl
+++ b/test/diff.jl
@@ -348,3 +348,41 @@ let
     @test isequal(expand_derivatives(Differential(t)(im*t)), im)
     @test isequal(expand_derivatives(Differential(t)(t^2 + im*t)), 2t + im)
 end
+
+## Vector Calculus
+
+
+@variables x[1:3] y[1:3]
+
+∇ = Nabla(x)
+
+out = ∇(y)
+@test all(isequal.(scalarize(out), [Differential(x[1])(y[1]), Differential(x[2])(y[2]), Differential(x[3])(y[3])]))
+
+
+out = ∇ ⋅ y
+@test isequal(out, Div(x)(y))
+@test isequal(out, Differential(x[1])(y[1]) + Differential(x[2])(y[2]) + Differential(x[3])(y[3]))
+
+out = ∇ × y
+@test isequal(out, Curl(x)(y))
+@test all(isequal.(scalarize(out), [Differential(x[2])(y[3]) - Differential(x[3])(y[2]),
+                    Differential(x[3])(y[1]) - Differential(x[1])(y[3]),
+                    Differential(x[1])(y[2]) - Differential(x[2])(y[1])]))
+
+
+
+out = (∇ ⋅ ∇)(y)
+@test isequal(out, Laplacian(x)(y))
+@test all(isequal.(scalarize(out), [Differential(x[1])(Differential(x[1])(y[1])) + Differential(x[2])(Differential(x[2])(y[1])) + Differential(x[3])(Differential(x[3])(y[1])), 
+                         Differential(x[1])(Differential(x[1])(y[2])) + Differential(x[2])(Differential(x[2])(y[2])) + Differential(x[3])(Differential(x[3])(y[2])), 
+                         Differential(x[1])(Differential(x[1])(y[3])) + Differential(x[2])(Differential(x[2])(y[3])) + Differential(x[3])(Differential(x[3])(y[3]))]))
+
+
+out = (∇ × ∇)(y)
+
+@test isequal(out, Curl(x)(Curl(x)(y)))
+@test all(isequal.(scalarize(out), [Differential(x[2])(Differential(x[3])(y[1])) - Differential(x[3])(Differential(x[2])(y[1])), 
+                         Differential(x[3])(Differential(x[1])(y[2])) - Differential(x[1])(Differential(x[3])(y[2])), 
+                         Differential(x[1])(Differential(x[2])(y[3])) - Differential(x[2])(Differential(x[1])(y[3]))]))
+

From 0559b065da576b7db5927e5dccd4a12c867fe724 Mon Sep 17 00:00:00 2001
From: Alex Jones <alex.lemony@gmail.com>
Date: Mon, 7 Aug 2023 17:11:11 +0100
Subject: [PATCH 07/12] import scalarize

---
 test/diff.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/diff.jl b/test/diff.jl
index 17dff5bc7..faabcf3b6 100644
--- a/test/diff.jl
+++ b/test/diff.jl
@@ -1,7 +1,7 @@
 using Symbolics
 using Test
 using IfElse
-using Symbolics: value
+using Symbolics: value, scalarize
 
 # Derivatives
 @variables t σ ρ β

From caa49cc791b84dacb8299c93865696e7d8339757 Mon Sep 17 00:00:00 2001
From: Alex Jones <alex.lemony@gmail.com>
Date: Mon, 7 Aug 2023 17:14:27 +0100
Subject: [PATCH 08/12] a comment

---
 src/diff.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/diff.jl b/src/diff.jl
index 131a1fe0d..f72bd27ac 100644
--- a/src/diff.jl
+++ b/src/diff.jl
@@ -875,4 +875,6 @@ Base.nameof(D::ArrayDifferentialOperator) = Symbol(D.name)
 function Base.:(==)(D1::ArrayDifferentialOperator, D2::ArrayDifferentialOperator) 
     @variables x[1:length(D1.vars)]
     all(scalarize(isequal.(D1.vars, D2.vars))) && all(scalarize(isequal.(D1(x), D2(x))))
-end
\ No newline at end of file
+end
+
+# TODO: Add simplification rules for dot and cross products to remove 0 terms and simplify
\ No newline at end of file

From 75fb18cb991e44c07e8466943acbc72e09c78850 Mon Sep 17 00:00:00 2001
From: Alex Jones <alex.lemony@gmail.com>
Date: Mon, 7 Aug 2023 17:19:51 +0100
Subject: [PATCH 09/12] remove duplicate method

---
 src/diff.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/diff.jl b/src/diff.jl
index f72bd27ac..664ece98c 100644
--- a/src/diff.jl
+++ b/src/diff.jl
@@ -867,8 +867,6 @@ end
 
 SymbolicUtils.promote_symtype(::ArrayDifferentialOperator, x) = x
 
-is_derivative(x) = istree(x) ? operation(x) isa ArrayDifferentialOperator : false
-
 Base.show(io::IO, D::ArrayDifferentialOperator) = print(io, D.name)
 Base.nameof(D::ArrayDifferentialOperator) = Symbol(D.name)
 

From 09943a48def7bbba04115ade7bb5864ff0e4a8a3 Mon Sep 17 00:00:00 2001
From: Alex Jones <alex.lemony@gmail.com>
Date: Wed, 9 Aug 2023 15:29:27 +0100
Subject: [PATCH 10/12] test fixes

---
 src/diff.jl  |   5 +-
 test/diff.jl | 295 +++++++++++++++++++++++++++------------------------
 2 files changed, 158 insertions(+), 142 deletions(-)

diff --git a/src/diff.jl b/src/diff.jl
index 664ece98c..56f9def12 100644
--- a/src/diff.jl
+++ b/src/diff.jl
@@ -801,6 +801,7 @@ struct ArrayDifferentialOperator <: ArrayOperator
     ArrayDifferentialOperator(vars, differentials, name) = new(vars, differentials, name)
 end
 Nabla(vars) = ArrayDifferentialOperator(value.(vars), map(Differential, scalarize(value.(vars))), "∇")
+const Grad = Nabla
 Div(vars) = (x) -> Nabla(vars) ⋅ x
 Curl(vars) = (x) -> Nabla(vars) × x
 Laplacian(vars) = Nabla(vars) ⋅ Nabla(vars)
@@ -809,7 +810,7 @@ Laplacian(vars) = Nabla(vars) ⋅ Nabla(vars)
 
 function (D::ArrayDifferentialOperator)(x::SymVec)
     @assert length(D.vars) == length(x) "Vector must be same length as vars in Operator $(D.name)."
-    @arrayop (i,) (D.differentials)[i](x[i]) term=D(y) 
+    @arrayop (i,) (D.differentials)[i](x[i]) term=D(x) 
 end
 (D::ArrayDifferentialOperator)(x::Arr) = Arr(D(value(x)))
 
@@ -854,7 +855,7 @@ end
 function LinearAlgebra.cross(D::ArrayDifferentialOperator, x::SymVec)
     @assert length(D.vars) == length(x) == 3 "Cross product is only defined in 3 dimensions."
     curl = crosscall(D.differentials, x)
-    @arrayop (i,) curl[i] term=D×y 
+    @arrayop (i,) curl[i] term=D×x 
 end
 LinearAlgebra.cross(D::ArrayDifferentialOperator, x::Arr) = Arr(D × value(x))
 
diff --git a/test/diff.jl b/test/diff.jl
index faabcf3b6..532f4cd5d 100644
--- a/test/diff.jl
+++ b/test/diff.jl
@@ -1,7 +1,7 @@
 using Symbolics
 using Test
 using IfElse
-using Symbolics: value, scalarize
+using Symbolics: value, scalarize, Nabla, Div, Curl, Laplacian
 
 # Derivatives
 @variables t σ ρ β
@@ -37,32 +37,32 @@ dcsch = D(csch(t))
 
 # Chain rule
 dsinsin = D(sin(sin(t)))
-test_equal(expand_derivatives(dsinsin), cos(sin(t))*cos(t))
+test_equal(expand_derivatives(dsinsin), cos(sin(t)) * cos(t))
 
-d1 = D(sin(t)*t)
-d2 = D(sin(t)*cos(t))
-@test isequal(expand_derivatives(d1), simplify(t*cos(t)+sin(t)))
-@test isequal(expand_derivatives(d2), cos(t)^2-sin(t)^2)
+d1 = D(sin(t) * t)
+d2 = D(sin(t) * cos(t))
+@test isequal(expand_derivatives(d1), simplify(t * cos(t) + sin(t)))
+@test isequal(expand_derivatives(d2), cos(t)^2 - sin(t)^2)
 
-eqs = [σ*(y-x),
-       x*(ρ-z)-y,
-       x*y - β*z]
+eqs = [σ * (y - x),
+	x * (ρ - z) - y,
+	x * y - β * z]
 jac = Symbolics.jacobian(eqs, [x, y, z])
-test_equal(jac[1,1], -1σ)
-test_equal(jac[1,2], σ)
-test_equal(jac[1,3], 0)
-test_equal(jac[2,1],  ρ - z)
-test_equal(jac[2,2], -1)
-test_equal(jac[2,3], -1x)
-test_equal(jac[3,1], y)
-test_equal(jac[3,2], x)
-test_equal(jac[3,3], -1β)
+test_equal(jac[1, 1], -1σ)
+test_equal(jac[1, 2], σ)
+test_equal(jac[1, 3], 0)
+test_equal(jac[2, 1], ρ - z)
+test_equal(jac[2, 2], -1)
+test_equal(jac[2, 3], -1x)
+test_equal(jac[3, 1], y)
+test_equal(jac[3, 2], x)
+test_equal(jac[3, 3], -1β)
 
 # issue #545
 z = t + t^2
 #test_equal(expand_derivatives(D(z)), 1 + t * 2)
 
-z = t-2t
+z = t - 2t
 #test_equal(expand_derivatives(D(z)), -1)
 
 # Variable dependence checking in differentiation
@@ -73,8 +73,8 @@ z = t-2t
 
 @variables x(t) y(t) z(t)
 
-@test isequal(expand_derivatives(D(x * y)), simplify(y*D(x) + x*D(y)))
-@test isequal(expand_derivatives(D(x * y)), simplify(D(x)*y + x*D(y)))
+@test isequal(expand_derivatives(D(x * y)), simplify(y * D(x) + x * D(y)))
+@test isequal(expand_derivatives(D(x * y)), simplify(D(x) * y + x * D(y)))
 
 @test isequal(expand_derivatives(D(2t)), 2)
 @test isequal(expand_derivatives(D(2x)), 2D(x))
@@ -88,16 +88,16 @@ z = t-2t
 @test all(iszero, Symbolics.gradient(42, [t, x, y, z]))
 @test all(iszero, Symbolics.hessian(42, [t, x, y, z]))
 @test isequal(Symbolics.jacobian([t, x, 42], [t, x]),
-              Num[1  0
-                  Differential(t)(x)           1
-                  0  0])
+	Num[1  0
+		Differential(t)(x)           1
+		0  0])
 
 # issue 252
 @variables beta, alpha, delta
 @variables x1, x2, x3
 
 # expression
-tmp = beta * (alpha * exp(x1) * x2 ^ (alpha - 1) + 1 - delta) / x3
+tmp = beta * (alpha * exp(x1) * x2^(alpha - 1) + 1 - delta) / x3
 # derivative w.r.t. x1 and x2
 t1 = Symbolics.gradient(tmp, [x1, x2])
 @test_nowarn Symbolics.gradient(t1[1], [beta])
@@ -111,7 +111,7 @@ using Symbolics
 @variables t x(t)
 ∂ₜ = Differential(t)
 ∂ₓ = Differential(x)
-L = .5 * ∂ₜ(x)^2 - .5 * x^2
+L = 0.5 * ∂ₜ(x)^2 - 0.5 * x^2
 @test isequal(expand_derivatives(∂ₓ(L)), -1 * x)
 @test isequal(expand_derivatives(Differential(x)(L) - ∂ₜ(Differential(∂ₜ(x))(L))), -1 * (∂ₜ(∂ₜ(x)) + x))
 @test isequal(expand_derivatives(Differential(x)(L) - ∂ₜ(Differential(∂ₜ(x))(L))), (-1 * x) - ∂ₜ(∂ₜ(x)))
@@ -123,9 +123,9 @@ L = .5 * ∂ₜ(x)^2 - .5 * x^2
 @variables u(..)
 Dy = Differential(y)
 Dx = Differential(x)
-dxyu = Dx(Dy(u(x,y)))
+dxyu = Dx(Dy(u(x, y)))
 @test isequal(expand_derivatives(dxyu), dxyu)
-dxxu = Dx(Dx(u(x,y)))
+dxxu = Dx(Dx(u(x, y)))
 @test isequal(expand_derivatives(dxxu), dxxu)
 
 using Symbolics, LinearAlgebra, SparseArrays
@@ -137,112 +137,114 @@ canonequal(a, b) = isequal(simplify(a), simplify(b))
 @variables t σ ρ β
 @variables x y z
 @test isequal(
-    (Differential(z) * Differential(y) * Differential(x))(t),
-    Differential(z)(Differential(y)(Differential(x)(t)))
+	(Differential(z) * Differential(y) * Differential(x))(t),
+	Differential(z)(Differential(y)(Differential(x)(t))),
 )
 
 @test canonequal(
-                 Symbolics.derivative(sin(cos(x)), x),
-                 -sin(x) * cos(cos(x))
-                )
+	Symbolics.derivative(sin(cos(x)), x),
+	-sin(x) * cos(cos(x)),
+)
 
 Symbolics.@register_symbolic no_der(x)
 @test canonequal(
-                 Symbolics.derivative([sin(cos(x)), hypot(x, no_der(x))], x),
-                 [
-                  -sin(x) * cos(cos(x)),
-                  x/hypot(x, no_der(x)) +
-                      no_der(x)*Differential(x)(no_der(x))/hypot(x, no_der(x))
-                 ]
-                )
+	Symbolics.derivative([sin(cos(x)), hypot(x, no_der(x))], x),
+	[
+		-sin(x) * cos(cos(x)),
+		x / hypot(x, no_der(x)) +
+		no_der(x) * Differential(x)(no_der(x)) / hypot(x, no_der(x)),
+	],
+)
 
 Symbolics.@register_symbolic intfun(x)::Int
 @test Symbolics.symtype(intfun(x).val) === Int
 
-eqs = [σ*(y-x),
-       x*(ρ-z)-y,
-       x*y - β*z]
+eqs = [σ * (y - x),
+	x * (ρ - z) - y,
+	x * y - β * z]
 
-∂ = Symbolics.jacobian(eqs,[x,y,z])
+∂ = Symbolics.jacobian(eqs, [x, y, z])
 for i in 1:3
-    ∇ = Symbolics.gradient(eqs[i],[x,y,z])
-    @test canonequal(∂[i,:],∇)
+	∇ = Symbolics.gradient(eqs[i], [x, y, z])
+	@test canonequal(∂[i, :], ∇)
 end
 
-@test all(canonequal.(Symbolics.gradient(eqs[1],[x,y,z]),[σ * -1,σ,0]))
-@test all(canonequal.(Symbolics.hessian(eqs[1],[x,y,z]),0))
+@test all(canonequal.(Symbolics.gradient(eqs[1], [x, y, z]), [σ * -1, σ, 0]))
+@test all(canonequal.(Symbolics.hessian(eqs[1], [x, y, z]), 0))
 
-du = [x^2, y^3, x^4, sin(y), x+y, x+z^2, z+x, x+y^2+sin(z)]
-reference_jac = sparse(Symbolics.jacobian(du, [x,y,z]))
+du = [x^2, y^3, x^4, sin(y), x + y, x + z^2, z + x, x + y^2 + sin(z)]
+reference_jac = sparse(Symbolics.jacobian(du, [x, y, z]))
 
-@test findnz(Symbolics.jacobian_sparsity(du, [x,y,z]))[[1,2]] == findnz(reference_jac)[[1,2]]
+@test findnz(Symbolics.jacobian_sparsity(du, [x, y, z]))[[1, 2]] == findnz(reference_jac)[[1, 2]]
 
 let
-    @variables t x(t) y(t) z(t)
-    @test Symbolics.exprs_occur_in([x,y,z], x^2*y) == [true, true, false]
+	@variables t x(t) y(t) z(t)
+	@test Symbolics.exprs_occur_in([x, y, z], x^2 * y) == [true, true, false]
 end
 
-@test isequal(Symbolics.sparsejacobian(du, [x,y,z]), reference_jac)
+@test isequal(Symbolics.sparsejacobian(du, [x, y, z]), reference_jac)
 
 
-function f!(res,u)
-    (x,y,z)=u
-    res.=[x^2, y^3, x^4, sin(y), x+y, x+z^2, z+x, x+y^2+sin(z)]
+function f!(res, u)
+	(x, y, z) = u
+	res .= [x^2, y^3, x^4, sin(y), x + y, x + z^2, z + x, x + y^2 + sin(z)]
 end
-function f1!(res,u,a,b;c)
-    (x,y,z)=u
-    res.=[a*x^2, y^3, b*x^4, sin(y), c*x+y, x+z^2, a*z+x, x+y^2+sin(z)]
+function f1!(res, u, a, b; c)
+	(x, y, z) = u
+	res .= [a * x^2, y^3, b * x^4, sin(y), c * x + y, x + z^2, a * z + x, x + y^2 + sin(z)]
 end
 
-input=rand(3)
-output=rand(8)
+input = rand(3)
+output = rand(8)
 
-findnz(Symbolics.jacobian_sparsity(f!, output, input))[[1,2]] == findnz(reference_jac)[[1,2]]
-findnz(Symbolics.jacobian_sparsity(f1!, output, input,1,2,c=3))[[1,2]] == findnz(reference_jac)[[1,2]]
+findnz(Symbolics.jacobian_sparsity(f!, output, input))[[1, 2]] == findnz(reference_jac)[[1, 2]]
+findnz(Symbolics.jacobian_sparsity(f1!, output, input, 1, 2, c = 3))[[1, 2]] == findnz(reference_jac)[[1, 2]]
 
-input = rand(2,2)
-function f2!(res,u,a,b,c)
-    (x,y,z)=u[1,1],u[2,1],u[3,1]
-    res.=[a*x^2, y^3, b*x^4, sin(y), c*x+y, x+z^2, a*z+x, x+y^2+sin(z)]
+input = rand(2, 2)
+function f2!(res, u, a, b, c)
+	(x, y, z) = u[1, 1], u[2, 1], u[3, 1]
+	res .= [a * x^2, y^3, b * x^4, sin(y), c * x + y, x + z^2, a * z + x, x + y^2 + sin(z)]
 end
 
-findnz(Symbolics.jacobian_sparsity(f!, output, input))[[1,2]] == findnz(reference_jac)[[1,2]]
+findnz(Symbolics.jacobian_sparsity(f!, output, input))[[1, 2]] == findnz(reference_jac)[[1, 2]]
 
 # Check for failures due to du[4] undefined
-function f_undef(du,u)
-    du[1] = u[1]
-    du[2] = u[2]
-    du[3] = u[3] + u[4]
+function f_undef(du, u)
+	du[1] = u[1]
+	du[2] = u[2]
+	du[3] = u[3] + u[4]
 end
 u0 = rand(4)
 du0 = similar(u0)
-sparsity_pattern = Symbolics.jacobian_sparsity(f_undef,du0,u0)
+sparsity_pattern = Symbolics.jacobian_sparsity(f_undef, du0, u0)
 udef_ref = sparse([1 0 0 0
-                    0 1 0 0
-                    0 0 1 1
-                    0 0 0 0])
-findnz(sparsity_pattern)[[1,2]] == findnz(udef_ref)[[1,2]]
+	 0 1 0 0
+	 0 0 1 1
+	 0 0 0 0])
+findnz(sparsity_pattern)[[1, 2]] == findnz(udef_ref)[[1, 2]]
 
 using Symbolics
 
-rosenbrock(X) = sum(1:length(X)-1) do i
-    100 * (X[i+1] - X[i]^2)^2 + (1 - X[i])^2
-end
-rosenbrock2(X,a;b) = sum(1:length(X)-1) do i
-    a * (X[i+1] - X[i]^2)^2 + (b - X[i])^2
-end
-
-@variables a,b
-X = [a,b]
+rosenbrock(X) =
+	sum(1:length(X)-1) do i
+		100 * (X[i+1] - X[i]^2)^2 + (1 - X[i])^2
+	end
+rosenbrock2(X, a; b) =
+	sum(1:length(X)-1) do i
+		a * (X[i+1] - X[i]^2)^2 + (b - X[i])^2
+	end
+
+@variables a, b
+X = [a, b]
 input = rand(2)
 
-spoly(x) = simplify(x, expand=true)
+spoly(x) = simplify(x, expand = true)
 rr = rosenbrock(X)
 
 reference_hes = Symbolics.hessian(rr, X)
 @test findnz(sparse(reference_hes))[1:2] == findnz(Symbolics.hessian_sparsity(rr, X))[1:2]
 @test findnz(sparse(reference_hes))[1:2] == findnz(Symbolics.hessian_sparsity(rosenbrock, input))[1:2]
-@test findnz(sparse(reference_hes))[1:2] == findnz(Symbolics.hessian_sparsity(rosenbrock2, input,100,b=1))[1:2]
+@test findnz(sparse(reference_hes))[1:2] == findnz(Symbolics.hessian_sparsity(rosenbrock2, input, 100, b = 1))[1:2]
 
 sp_hess = Symbolics.sparsehessian(rr, X)
 @test findnz(sparse(reference_hes))[1:2] == findnz(sp_hess)[1:2]
@@ -252,24 +254,24 @@ sp_hess = Symbolics.sparsehessian(rr, X)
 @variables t x(t)[1:4] ẋ(t)[1:4]
 expression = sin(x[1] + x[2] + x[3] + x[4]) |> Differential(t) |> expand_derivatives
 expression2 = substitute(expression, Dict(collect(Differential(t).(x) .=> ẋ)))
-@test isequal(expression2, (ẋ[1] + ẋ[2] + ẋ[3] + ẋ[4])*cos(x[1] + x[2] + x[3] + x[4]))
+@test isequal(expression2, (ẋ[1] + ẋ[2] + ẋ[3] + ẋ[4]) * cos(x[1] + x[2] + x[3] + x[4]))
 
 @test isequal(
-    Symbolics.derivative(IfElse.ifelse(signbit(b), b^2, sqrt(b)), b),
-    IfElse.ifelse(signbit(b), 2b,(SymbolicUtils.unstable_pow(2Symbolics.unwrap(sqrt(b)), -1)))
+	Symbolics.derivative(IfElse.ifelse(signbit(b), b^2, sqrt(b)), b),
+	IfElse.ifelse(signbit(b), 2b, (SymbolicUtils.unstable_pow(2Symbolics.unwrap(sqrt(b)), -1))),
 )
 
 # Chain rule
 #
 let
-    @variables t x(..) y(..) z(..)
-    Dt = Differential(t)
+	@variables t x(..) y(..) z(..)
+	Dt = Differential(t)
 
-    @test isequal(expand_derivatives(Dt(x(y(t)))),
-                  Dt(y(t))*Differential(y(t))(x(y(t))))
+	@test isequal(expand_derivatives(Dt(x(y(t)))),
+		Dt(y(t)) * Differential(y(t))(x(y(t))))
 
-    @test isequal(expand_derivatives(Dt(x(y(t), z(t)))),
-                  Dt(y(t))*Differential(y(t))(x(y(t), z(t))) + Dt(z(t))*Differential(z(t))(x(y(t), z(t))))
+	@test isequal(expand_derivatives(Dt(x(y(t), z(t)))),
+		Dt(y(t)) * Differential(y(t))(x(y(t), z(t))) + Dt(z(t)) * Differential(z(t))(x(y(t), z(t))))
 end
 
 @variables x y
@@ -282,71 +284,74 @@ D = Differential(t)
 
 eqs = [D(x) ~ x, D(y) ~ y + x]
 
-sub_eqs = substitute(eqs, Dict([D(x)=>D(x), x=>1]))
+sub_eqs = substitute(eqs, Dict([D(x) => D(x), x => 1]))
 @test sub_eqs == [D(x) ~ 1, D(y) ~ 1 + y]
 
 @variables x y
-@test substitute([x + y; x - y], Dict(x=>1, y=>2)) == [3, -1]
+@test substitute([x + y; x - y], Dict(x => 1, y => 2)) == [3, -1]
 
 
 # 530#discussion_r825125589
 let
-    using Symbolics
-    @variables u[1:2] y[1:1] t
-    u = collect(u)
-    y = collect(y)
-    @test isequal(Symbolics.jacobian([u;u[1]^2; y], u), Num[1 0
-                                                            0 1
-                                                            2u[1] 0
-                                                            0 0])
+	using Symbolics
+	@variables u[1:2] y[1:1] t
+	u = collect(u)
+	y = collect(y)
+	@test isequal(
+		Symbolics.jacobian([u; u[1]^2; y], u),
+		Num[1 0
+			0 1
+			2u[1] 0
+			0 0],
+	)
 end
 
 # make sure derivative(x[1](t), y) does not fail
 let
-    @variables t a(t)
-    vars = collect(@variables(x(t)[1:1])[1])
-    ps = collect(@variables(ps[1:1])[1])
-    @test Symbolics.derivative(ps[1], vars[1]) == 0
-    @test Symbolics.derivative(ps[1], a) == 0
-    @test Symbolics.derivative(x[1], a) == 0
+	@variables t a(t)
+	vars = collect(@variables(x(t)[1:1])[1])
+	ps = collect(@variables(ps[1:1])[1])
+	@test Symbolics.derivative(ps[1], vars[1]) == 0
+	@test Symbolics.derivative(ps[1], a) == 0
+	@test Symbolics.derivative(x[1], a) == 0
 end
 
 # 580
 let
-    @variables x[1:3]
-    y = [sin.(x); cos.(x)]
-    dj = Symbolics.jacobian(y, x)
-    @test !iszero(dj)
-    @test isequal(dj, Symbolics.jacobian(Symbolics.scalarize.(y), x))
-    sj = Symbolics.sparsejacobian(y, x)
-    @test !iszero(sj)
-    @test isequal(sj, Symbolics.jacobian(Symbolics.scalarize.(y), x))
+	@variables x[1:3]
+	y = [sin.(x); cos.(x)]
+	dj = Symbolics.jacobian(y, x)
+	@test !iszero(dj)
+	@test isequal(dj, Symbolics.jacobian(Symbolics.scalarize.(y), x))
+	sj = Symbolics.sparsejacobian(y, x)
+	@test !iszero(sj)
+	@test isequal(sj, Symbolics.jacobian(Symbolics.scalarize.(y), x))
 end
 
 # substituting iv of differentials
 @variables t t2 x(t)
 D = Differential(t)
 ex = D(x)
-ex2 = substitute(ex, [t=>t2])
+ex2 = substitute(ex, [t => t2])
 @test isequal(operation(Symbolics.unwrap(ex2)).x, t2)
-ex3 = substitute(D(x) * 2 + x / t, [t=>t2])
+ex3 = substitute(D(x) * 2 + x / t, [t => t2])
 xt2 = substitute(x, [t => t2])
 @test isequal(ex3, xt2 / t2 + 2Differential(t2)(xt2))
 
 # 581
 #
 let
-    @variables x(t)[1:3]
-    @test iszero(Symbolics.derivative(x[1], x[2]))
+	@variables x(t)[1:3]
+	@test iszero(Symbolics.derivative(x[1], x[2]))
 end
 
 #908
 #
 let
-    using Symbolics
-    @variables t
-    @test isequal(expand_derivatives(Differential(t)(im*t)), im)
-    @test isequal(expand_derivatives(Differential(t)(t^2 + im*t)), 2t + im)
+	using Symbolics
+	@variables t
+	@test isequal(expand_derivatives(Differential(t)(im * t)), im)
+	@test isequal(expand_derivatives(Differential(t)(t^2 + im * t)), 2t + im)
 end
 
 ## Vector Calculus
@@ -367,22 +372,32 @@ out = ∇ ⋅ y
 out = ∇ × y
 @test isequal(out, Curl(x)(y))
 @test all(isequal.(scalarize(out), [Differential(x[2])(y[3]) - Differential(x[3])(y[2]),
-                    Differential(x[3])(y[1]) - Differential(x[1])(y[3]),
-                    Differential(x[1])(y[2]) - Differential(x[2])(y[1])]))
+	Differential(x[3])(y[1]) - Differential(x[1])(y[3]),
+	Differential(x[1])(y[2]) - Differential(x[2])(y[1])]))
 
 
 
 out = (∇ ⋅ ∇)(y)
-@test isequal(out, Laplacian(x)(y))
-@test all(isequal.(scalarize(out), [Differential(x[1])(Differential(x[1])(y[1])) + Differential(x[2])(Differential(x[2])(y[1])) + Differential(x[3])(Differential(x[3])(y[1])), 
-                         Differential(x[1])(Differential(x[1])(y[2])) + Differential(x[2])(Differential(x[2])(y[2])) + Differential(x[3])(Differential(x[3])(y[2])), 
-                         Differential(x[1])(Differential(x[1])(y[3])) + Differential(x[2])(Differential(x[2])(y[3])) + Differential(x[3])(Differential(x[3])(y[3]))]))
+@test all(isequal.(scalarize(out), scalarize(Laplacian(x)(y))))
+@test all(
+	isequal.(
+		scalarize(out),
+		[Differential(x[1])(Differential(x[1])(y[1])) + Differential(x[2])(Differential(x[2])(y[1])) + Differential(x[3])(Differential(x[3])(y[1])),
+			Differential(x[1])(Differential(x[1])(y[2])) + Differential(x[2])(Differential(x[2])(y[2])) + Differential(x[3])(Differential(x[3])(y[2])),
+			Differential(x[1])(Differential(x[1])(y[3])) + Differential(x[2])(Differential(x[2])(y[3])) + Differential(x[3])(Differential(x[3])(y[3]))],
+	),
+)
 
 
 out = (∇ × ∇)(y)
 
-@test isequal(out, Curl(x)(Curl(x)(y)))
-@test all(isequal.(scalarize(out), [Differential(x[2])(Differential(x[3])(y[1])) - Differential(x[3])(Differential(x[2])(y[1])), 
-                         Differential(x[3])(Differential(x[1])(y[2])) - Differential(x[1])(Differential(x[3])(y[2])), 
-                         Differential(x[1])(Differential(x[2])(y[3])) - Differential(x[2])(Differential(x[1])(y[3]))]))
+@test_broken isequal(out, Curl(x)(Curl(x)(y)))
+@test all(
+	isequal.(
+		scalarize(out),
+		[Differential(x[2])(Differential(x[3])(y[1])) - Differential(x[3])(Differential(x[2])(y[1])),
+			Differential(x[3])(Differential(x[1])(y[2])) - Differential(x[1])(Differential(x[3])(y[2])),
+			Differential(x[1])(Differential(x[2])(y[3])) - Differential(x[2])(Differential(x[1])(y[3]))],
+	),
+)
 

From af60a3538b016cf6c32ff7be5ab63c5334a58df6 Mon Sep 17 00:00:00 2001
From: Alex Jones <alex.lemony@gmail.com>
Date: Fri, 25 Aug 2023 18:07:31 +0100
Subject: [PATCH 11/12] removed extra types

---
 src/diff.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/diff.jl b/src/diff.jl
index 56f9def12..363fc52d3 100644
--- a/src/diff.jl
+++ b/src/diff.jl
@@ -1,5 +1,4 @@
-abstract type AbstractOperator <: Function end
-abstract type Operator <: AbstractOperator end
+abstract type Operator <: Function end
 
 """
 $(TYPEDEF)
@@ -791,12 +790,13 @@ end
 #######################################################################################################################
 # Vector Calculus
 #######################################################################################################################
-abstract type ArrayOperator end
 
-struct ArrayDifferentialOperator <: ArrayOperator
-    """The variables to differentiate with resp≈ect to."""
+struct ArrayDifferentialOperator 
+    """The variables to differentiate with respect to."""
     vars
+    """The differentials, can be other fucntions if composite"""
     differentials
+    """name"""
     name
     ArrayDifferentialOperator(vars, differentials, name) = new(vars, differentials, name)
 end

From 068ad6455a04da9fab471c0a870fbfcad9b7a10d Mon Sep 17 00:00:00 2001
From: Alex Jones <alex.lemony@gmail.com>
Date: Mon, 4 Sep 2023 16:07:17 +0100
Subject: [PATCH 12/12] make ADO a subtype of Operator

---
 src/diff.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/diff.jl b/src/diff.jl
index 363fc52d3..3a5eb07bd 100644
--- a/src/diff.jl
+++ b/src/diff.jl
@@ -791,10 +791,10 @@ end
 # Vector Calculus
 #######################################################################################################################
 
-struct ArrayDifferentialOperator 
+struct ArrayDifferentialOperator <: Operator
     """The variables to differentiate with respect to."""
     vars
-    """The differentials, can be other fucntions if composite"""
+    """The differentials, can be other functions if composite"""
     differentials
     """name"""
     name