From b5d0fbd314314fb5f35405c38b6f302f82ed3037 Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Thu, 25 Sep 2025 20:01:19 +0200
Subject: [PATCH 01/19] Add norm symbol (as a multivariate operator). To be
 moved

It should be moved to univariate operators, but this has to be changed in MOI first
---
 src/reverse_mode.jl | 32 ++++++++++++++++++++++++++++++++
 src/sizes.jl        |  2 ++
 test/ArrayDiff.jl   | 41 ++++++++++++++++++++++++++++++++++++++++-
 3 files changed, 74 insertions(+), 1 deletion(-)

diff --git a/src/reverse_mode.jl b/src/reverse_mode.jl
index baab908..0106b0c 100644
--- a/src/reverse_mode.jl
+++ b/src/reverse_mode.jl
@@ -247,6 +247,22 @@ function _forward_eval(
                     tmp_dot += v1 * v2
                 end
                 @s f.forward_storage[k] = tmp_dot
+            elseif node.index == 14 # norm 
+                ix = children_arr[children_indices[1]]
+                tmp_norm_squared = zero(T)
+                for j in _eachindex(f.sizes, ix)
+                    v = @j f.forward_storage[ix]
+                    tmp_norm_squared += v * v
+                end
+                @s f.forward_storage[k] = sqrt(tmp_norm_squared)
+                for j in _eachindex(f.sizes, ix)
+                    v = @j f.forward_storage[ix]
+                    if tmp_norm_squared == 0
+                        @j f.partials_storage[ix] = zero(T)
+                    else
+                        @j f.partials_storage[ix] = v / @s f.forward_storage[k]
+                    end
+                end
             else # atan, min, max
                 f_input = _UnsafeVectorView(d.jac_storage, N)
                 ∇f = _UnsafeVectorView(d.user_output_buffer, N)
@@ -379,6 +395,22 @@ function _reverse_eval(f::_SubexpressionStorage)
                         end
                     end
                     continue
+                elseif op == :norm
+                    # Node `k` is scalar, the jacobian w.r.t. the vectorized input
+                    # child is a row vector whose entries are stored in `f.partials_storage`
+                    rev_parent = @s f.reverse_storage[k]
+                    for j in
+                        _eachindex(f.sizes, children_arr[children_indices[1]])
+                        ix = children_arr[children_indices[1]]
+                        partial = @j f.partials_storage[ix]
+                        val = ifelse(
+                            rev_parent == 0.0 && !isfinite(partial),
+                            rev_parent,
+                            rev_parent * partial,
+                        )
+                        @j f.reverse_storage[ix] = val  
+                    end
+                    continue
                 end
             end
         elseif node.type != MOI.Nonlinear.NODE_CALL_UNIVARIATE
diff --git a/src/sizes.jl b/src/sizes.jl
index f5ca9dd..9fcb562 100644
--- a/src/sizes.jl
+++ b/src/sizes.jl
@@ -181,6 +181,8 @@ function _infer_sizes(
                 _add_size!(sizes, k, (N,))
             elseif op == :dot
                 # TODO assert all arguments have same size
+            elseif op == :norm
+                # TODO actually norm should be moved to univariate
             elseif op == :+ || op == :-
                 # TODO assert all arguments have same size
                 _copy_size!(sizes, k, children_arr[first(children_indices)])
diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index 6dd1f6a..a909425 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -64,6 +64,45 @@ function test_objective_dot_bivariate()
     return
 end
 
+function test_objective_norm_univariate()
+    model = Nonlinear.Model()
+    x = MOI.VariableIndex(1)
+    Nonlinear.set_objective(model, :(norm([$x])))
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 1, 0]
+    @test sizes.size_offset == [0, 0, 0]
+    @test sizes.size == [1]
+    @test sizes.storage_offset == [0, 1, 2, 3]
+    x = [1.2]
+    @test MOI.eval_objective(evaluator, x) == abs(x[1])
+    g = ones(1)
+    MOI.eval_objective_gradient(evaluator, g, x)
+    @test g[1] == sign(x[1])
+    return
+end
+
+function test_objective_norm_bivariate()
+    model = Nonlinear.Model()
+    x = MOI.VariableIndex(1)
+    y = MOI.VariableIndex(2)
+    Nonlinear.set_objective(model, :(norm([$x, $y])))
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 1, 0, 0]
+    @test sizes.size_offset == [0, 0, 0, 0]
+    @test sizes.size == [2]
+    @test sizes.storage_offset == [0, 1, 3, 4, 5]
+    x = [3.0, 4.0]
+    @test MOI.eval_objective(evaluator, x) == 5.0
+    g = ones(2)
+    MOI.eval_objective_gradient(evaluator, g, x)
+    @test g == x / 5.0
+    return
+end
+
 end  # module
 
-TestArrayDiff.runtests()
+TestArrayDiff.runtests()
\ No newline at end of file

From 65fa6eff71b8bc07f7d05d5d7c2c598d54e9903d Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Thu, 25 Sep 2025 21:35:39 +0200
Subject: [PATCH 02/19] Add a test falling in the norm=0 - branch

Plus, does norm really need to be moved? I've changed my mind, it's multivar
---
 test/ArrayDiff.jl | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index a909425..d7b8efe 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -100,6 +100,11 @@ function test_objective_norm_bivariate()
     g = ones(2)
     MOI.eval_objective_gradient(evaluator, g, x)
     @test g == x / 5.0
+    y = [0.0, 0.0]
+    @test MOI.eval_objective(evaluator, y) == 0.0
+    g = ones(2)
+    MOI.eval_objective_gradient(evaluator, g, y)
+    @test g == [0.0, 0.0]
     return
 end
 

From fe207a5ca5a3fda9de8ce2faa65913214c0a2be2 Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Sat, 27 Sep 2025 18:35:08 +0200
Subject: [PATCH 03/19] Implement hcat evaluation WIP

Notes: * Gradient is not correct yet (test fails)
* hcat cannot be evaluated as it is, since the final output must be scalar and we do not have any operation to reduce the matrix to scalar at the moment
---
 src/reverse_mode.jl | 53 +++++++++++++++++++++++++++++++++++++++++++++
 src/sizes.jl        | 10 +++++++++
 test/ArrayDiff.jl   | 32 ++++++++++++++++++++++++++-
 3 files changed, 94 insertions(+), 1 deletion(-)

diff --git a/src/reverse_mode.jl b/src/reverse_mode.jl
index baab908..268ecd1 100644
--- a/src/reverse_mode.jl
+++ b/src/reverse_mode.jl
@@ -247,6 +247,24 @@ function _forward_eval(
                     tmp_dot += v1 * v2
                 end
                 @s f.forward_storage[k] = tmp_dot
+            elseif node.index == 12 # hcat
+                idx1, idx2 = children_indices
+                ix1 = children_arr[idx1]
+                ix2 = children_arr[idx2]
+                nb_cols1 = f.sizes.ndims[ix1] <= 1 ? 1 : _size(f.sizes, ix1, 2)
+                col_size = _size(f.sizes, k, 1)
+                for j in _eachindex(f.sizes, k)
+                    col = (j - 1) ÷ col_size + 1
+                    if col <= nb_cols1
+                        @j f.partials_storage[ix1] = one(T)
+                        val = @j f.forward_storage[ix1]
+                        @j f.forward_storage[k] = val
+                    else
+                        @j f.partials_storage[ix2] = one(T)
+                        val = @j f.forward_storage[ix2]
+                        @j f.forward_storage[k] = val
+                    end
+                end
             else # atan, min, max
                 f_input = _UnsafeVectorView(d.jac_storage, N)
                 ∇f = _UnsafeVectorView(d.user_output_buffer, N)
@@ -323,6 +341,18 @@ function _forward_eval(
             f.partials_storage[rhs] = zero(T)
         end
     end
+    # This function is written assuming that the final output is scalar.
+    # Therefore cannot return the matrix, so I guess I return it's first entry only, 
+    # as long as sum or matx-vect products are not implemented.
+
+    #println("Last node ", f.nodes[1].index)
+    #if f.nodes[1].index == 12
+    #    mtx = reshape(
+    #        f.forward_storage[_storage_range(f.sizes, 1)],
+    #        f.sizes.size[1:f.sizes.ndims[1]]...,
+    #    )
+    #    return mtx
+    #end
     return f.forward_storage[1]
 end
 
@@ -379,6 +409,29 @@ function _reverse_eval(f::_SubexpressionStorage)
                         end
                     end
                     continue
+                elseif op == :hcat
+                    total_cols = 0
+                    for c_idx in children_indices
+                        total_cols += f.sizes.ndims[children_arr[c_idx]] <= 1 ?
+                            1 : _size(f.sizes, children_arr[c_idx], 2)
+                    end
+                    row_size = _size(f.sizes, k, 1)
+                    col_start = 1
+                    for c_idx in children_indices
+                        child = children_arr[c_idx]
+                        child_col_size = f.sizes.ndims[child] <= 1 ?
+                            1 : _size(f.sizes, child, 2)
+                        for i in 1:row_size
+                            for j in 1:child_col_size
+                                rev_parent_j =
+                                    @j f.reverse_storage[k]
+                                # partial is 1 so we can ignore it
+                                @j f.reverse_storage[child] = rev_parent_j
+                            end
+                        end
+                        col_start += child_col_size
+                    end
+                    continue
                 end
             end
         elseif node.type != MOI.Nonlinear.NODE_CALL_UNIVARIATE
diff --git a/src/sizes.jl b/src/sizes.jl
index f5ca9dd..f5bd67a 100644
--- a/src/sizes.jl
+++ b/src/sizes.jl
@@ -184,6 +184,16 @@ function _infer_sizes(
             elseif op == :+ || op == :-
                 # TODO assert all arguments have same size
                 _copy_size!(sizes, k, children_arr[first(children_indices)])
+            elseif op == :hcat
+                total_cols = 0
+                for c_idx in children_indices
+                    total_cols += sizes.ndims[children_arr[c_idx]] <= 1 ?
+                        1 : _size(sizes, children_arr[c_idx], 2)
+                end
+                child_shape = _size(sizes, children_arr[first(children_indices)])
+                shape = sizes.ndims[children_arr[first(children_indices)]] <= 2 ?
+                    (child_shape[1], total_cols) : (child_shape[1], total_cols, child_shape[3:end]...)
+                _add_size!(sizes, k, tuple(shape...))
             elseif op == :*
                 # TODO assert compatible sizes and all ndims should be 0 or 2
                 first_matrix = findfirst(children_indices) do i
diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index 6dd1f6a..7fedbb6 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -64,6 +64,36 @@ function test_objective_dot_bivariate()
     return
 end
 
+function test_objective_matrix()
+    model = Nonlinear.Model()
+    x1 = MOI.VariableIndex(1)
+    x2 = MOI.VariableIndex(2)
+    x3 = MOI.VariableIndex(3)
+    x4 = MOI.VariableIndex(4)
+    Nonlinear.set_objective(
+        model,
+        :(hcat([$x1, $x3], [$x2, $x4]))
+    )
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [2, 1, 0, 0, 1, 0, 0]
+    @test sizes.size_offset == [2, 1, 0, 0, 0, 0, 0]
+    @test sizes.size == [2, 2, 2, 2]
+    @test sizes.storage_offset == [0, 4, 6, 7, 8, 10, 11, 12]
+    x1 = 1.0
+    x2 = 2.0
+    x3 = 3.0
+    x4 = 4.0
+    println(MOI.eval_objective(evaluator, [x1, x2, x3, x4]))
+    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) ==
+          1.0
+    g = ones(4)
+    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
+    @test g == [1.0, 0.0, 0.0, 0.0]
+    return
+end
+
 end  # module
 
-TestArrayDiff.runtests()
+TestArrayDiff.runtests()
\ No newline at end of file

From b59f8e6a80dfb9dc00d23ff208b0329e3e1ad196 Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Sun, 28 Sep 2025 16:23:05 +0200
Subject: [PATCH 04/19] Fix hcat eval and gradient

0dim still to be fixed (hcat scalars into a row)
---
 src/reverse_mode.jl | 62 ++++++++++++++++++++++-----------------------
 test/ArrayDiff.jl   | 16 ++++++------
 2 files changed, 39 insertions(+), 39 deletions(-)

diff --git a/src/reverse_mode.jl b/src/reverse_mode.jl
index 268ecd1..632a9db 100644
--- a/src/reverse_mode.jl
+++ b/src/reverse_mode.jl
@@ -252,18 +252,16 @@ function _forward_eval(
                 ix1 = children_arr[idx1]
                 ix2 = children_arr[idx2]
                 nb_cols1 = f.sizes.ndims[ix1] <= 1 ? 1 : _size(f.sizes, ix1, 2)
-                col_size = _size(f.sizes, k, 1)
-                for j in _eachindex(f.sizes, k)
-                    col = (j - 1) ÷ col_size + 1
-                    if col <= nb_cols1
-                        @j f.partials_storage[ix1] = one(T)
-                        val = @j f.forward_storage[ix1]
-                        @j f.forward_storage[k] = val
-                    else
-                        @j f.partials_storage[ix2] = one(T)
-                        val = @j f.forward_storage[ix2]
-                        @j f.forward_storage[k] = val
-                    end
+                col_size = f.sizes.ndims[ix1] == 0 ? 1 : _size(f.sizes, k, 1)
+                for j in _eachindex(f.sizes, ix1)
+                    @j f.partials_storage[ix1] = one(T)
+                    val = @j f.forward_storage[ix1]
+                    @j f.forward_storage[k] = val
+                end
+                for j in _eachindex(f.sizes, ix2)
+                    @j f.partials_storage[ix2] = one(T)
+                    val = @j f.forward_storage[ix2]
+                    _setindex!(f.forward_storage, val, f.sizes, k, j + nb_cols1 * col_size)
                 end
             else # atan, min, max
                 f_input = _UnsafeVectorView(d.jac_storage, N)
@@ -410,26 +408,28 @@ function _reverse_eval(f::_SubexpressionStorage)
                     end
                     continue
                 elseif op == :hcat
-                    total_cols = 0
-                    for c_idx in children_indices
-                        total_cols += f.sizes.ndims[children_arr[c_idx]] <= 1 ?
-                            1 : _size(f.sizes, children_arr[c_idx], 2)
+                    idx1, idx2 = children_indices
+                    ix1 = children_arr[idx1]
+                    ix2 = children_arr[idx2]
+                    nb_cols1 = f.sizes.ndims[ix1] <= 1 ? 1 : _size(f.sizes, ix1, 2)
+                    col_size = f.sizes.ndims[ix1] == 0 ? 1 : _size(f.sizes, k, 1)
+                    for j in _eachindex(f.sizes, ix1)
+                        partial = @j f.partials_storage[ix1]
+                        val = ifelse(
+                            _getindex(f.reverse_storage, f.sizes, k, j) == 0.0 && !isfinite(partial),
+                            _getindex(f.reverse_storage, f.sizes, k, j),
+                            _getindex(f.reverse_storage, f.sizes, k, j) * partial,
+                        )
+                        @j f.reverse_storage[ix1] = val
                     end
-                    row_size = _size(f.sizes, k, 1)
-                    col_start = 1
-                    for c_idx in children_indices
-                        child = children_arr[c_idx]
-                        child_col_size = f.sizes.ndims[child] <= 1 ?
-                            1 : _size(f.sizes, child, 2)
-                        for i in 1:row_size
-                            for j in 1:child_col_size
-                                rev_parent_j =
-                                    @j f.reverse_storage[k]
-                                # partial is 1 so we can ignore it
-                                @j f.reverse_storage[child] = rev_parent_j
-                            end
-                        end
-                        col_start += child_col_size
+                    for j in _eachindex(f.sizes, ix2)
+                        partial = @j f.partials_storage[ix2]
+                        val = ifelse(
+                            _getindex(f.reverse_storage, f.sizes, k, j + nb_cols1 * col_size) == 0.0 && !isfinite(partial),
+                            _getindex(f.reverse_storage, f.sizes, k, j + nb_cols1 * col_size),
+                            _getindex(f.reverse_storage, f.sizes, k, j + nb_cols1 * col_size) * partial,
+                        )
+                        @j f.reverse_storage[ix2] = val
                     end
                     continue
                 end
diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index 7fedbb6..6b16ea9 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -64,7 +64,7 @@ function test_objective_dot_bivariate()
     return
 end
 
-function test_objective_matrix()
+function test_objective_hcat_1dim()
     model = Nonlinear.Model()
     x1 = MOI.VariableIndex(1)
     x2 = MOI.VariableIndex(2)
@@ -72,25 +72,25 @@ function test_objective_matrix()
     x4 = MOI.VariableIndex(4)
     Nonlinear.set_objective(
         model,
-        :(hcat([$x1, $x3], [$x2, $x4]))
+        :(dot(hcat([$x1], [$x3]), hcat([$x2], [$x4])))
     )
     evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
     MOI.initialize(evaluator, [:Grad])
     sizes = evaluator.backend.objective.expr.sizes
-    @test sizes.ndims == [2, 1, 0, 0, 1, 0, 0]
-    @test sizes.size_offset == [2, 1, 0, 0, 0, 0, 0]
-    @test sizes.size == [2, 2, 2, 2]
-    @test sizes.storage_offset == [0, 4, 6, 7, 8, 10, 11, 12]
+    @test sizes.ndims == [0, 2, 1, 0, 1, 0, 2, 1, 0, 1, 0]
+    @test sizes.size_offset == [0, 6, 5, 0, 4, 0, 2, 1, 0, 0, 0]
+    @test sizes.size == [1, 1, 1, 2, 1, 1, 1, 2]
+    @test sizes.storage_offset == [0, 1, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13]
     x1 = 1.0
     x2 = 2.0
     x3 = 3.0
     x4 = 4.0
     println(MOI.eval_objective(evaluator, [x1, x2, x3, x4]))
     @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) ==
-          1.0
+          14.0
     g = ones(4)
     MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
-    @test g == [1.0, 0.0, 0.0, 0.0]
+    @test g == [2.0, 1.0, 4.0, 3.0]
     return
 end
 

From 0e90752fe1e67e66af6104fa4f2a7792ab940624 Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Sun, 28 Sep 2025 18:19:21 +0200
Subject: [PATCH 05/19] Fix 0dim

Meaning a row vector can now be written with [1 2] or hcat(1, 2)
---
 src/sizes.jl      | 11 ++++++++---
 test/ArrayDiff.jl | 30 ++++++++++++++++++++++++++++++
 2 files changed, 38 insertions(+), 3 deletions(-)

diff --git a/src/sizes.jl b/src/sizes.jl
index f5bd67a..4a44332 100644
--- a/src/sizes.jl
+++ b/src/sizes.jl
@@ -190,9 +190,14 @@ function _infer_sizes(
                     total_cols += sizes.ndims[children_arr[c_idx]] <= 1 ?
                         1 : _size(sizes, children_arr[c_idx], 2)
                 end
-                child_shape = _size(sizes, children_arr[first(children_indices)])
-                shape = sizes.ndims[children_arr[first(children_indices)]] <= 2 ?
-                    (child_shape[1], total_cols) : (child_shape[1], total_cols, child_shape[3:end]...)
+                if sizes.ndims[children_arr[first(children_indices)]] == 0
+                    shape = (1, total_cols)
+                elseif sizes.ndims[children_arr[first(children_indices)]] <= 2
+                    shape = ( _size(sizes, children_arr[first(children_indices)], 1), total_cols)
+                else 
+                    child_shape = _size(sizes, children_arr[first(children_indices)])
+                    shape = (child_shape[1], total_cols, child_shape[3:end]...)
+                end
                 _add_size!(sizes, k, tuple(shape...))
             elseif op == :*
                 # TODO assert compatible sizes and all ndims should be 0 or 2
diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index 6b16ea9..c388443 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -64,6 +64,36 @@ function test_objective_dot_bivariate()
     return
 end
 
+function test_objective_hcat_0dim()
+    model = Nonlinear.Model()
+    x1 = MOI.VariableIndex(1)
+    x2 = MOI.VariableIndex(2)
+    x3 = MOI.VariableIndex(3)
+    x4 = MOI.VariableIndex(4)
+    Nonlinear.set_objective(
+        model,
+        :(dot([$x1 $x3], [$x2 $x4]))
+    )
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 2, 0, 0, 2, 0, 0]
+    @test sizes.size_offset == [0, 2, 0, 0, 0, 0, 0]
+    @test sizes.size == [1, 2, 1, 2]
+    @test sizes.storage_offset == [0, 1, 3, 4, 5, 7, 8, 9]
+    x1 = 1.0
+    x2 = 2.0
+    x3 = 3.0
+    x4 = 4.0
+    println(MOI.eval_objective(evaluator, [x1, x2, x3, x4]))
+    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) ==
+          14.0
+    g = ones(4)
+    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
+    @test g == [2.0, 1.0, 4.0, 3.0]
+    return
+end
+
 function test_objective_hcat_1dim()
     model = Nonlinear.Model()
     x1 = MOI.VariableIndex(1)

From 7842545a2e933448edf31d8ca1845d8f78569206 Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Mon, 29 Sep 2025 09:20:32 +0200
Subject: [PATCH 06/19] Add Revise to toml (temporary)

---
 test/Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/Project.toml b/test/Project.toml
index 7ed96af..0e679fe 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,5 +2,6 @@
 ArrayDiff = "c45fa1ca-6901-44ac-ae5b-5513a4852d50"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

From 7576b515beaaf3db573b7cc504778bf3b817b0f8 Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Mon, 29 Sep 2025 09:42:10 +0200
Subject: [PATCH 07/19] Fix quick pb I introduced when merging

---
 src/reverse_mode.jl | 2 ++
 test/ArrayDiff.jl   | 4 +++-
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/reverse_mode.jl b/src/reverse_mode.jl
index 2238614..c12e4be 100644
--- a/src/reverse_mode.jl
+++ b/src/reverse_mode.jl
@@ -262,6 +262,7 @@ function _forward_eval(
                     @j f.partials_storage[ix2] = one(T)
                     val = @j f.forward_storage[ix2]
                     _setindex!(f.forward_storage, val, f.sizes, k, j + nb_cols1 * col_size)
+                end
             elseif node.index == 14 # norm 
                 ix = children_arr[children_indices[1]]
                 tmp_norm_squared = zero(T)
@@ -445,6 +446,7 @@ function _reverse_eval(f::_SubexpressionStorage)
                             _getindex(f.reverse_storage, f.sizes, k, j + nb_cols1 * col_size) * partial,
                         )
                         @j f.reverse_storage[ix2] = val
+                    end
                 elseif op == :norm
                     # Node `k` is scalar, the jacobian w.r.t. the vectorized input
                     # child is a row vector whose entries are stored in `f.partials_storage`
diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index 3b4776d..093f0c9 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -121,6 +121,8 @@ function test_objective_hcat_1dim()
     g = ones(4)
     MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
     @test g == [2.0, 1.0, 4.0, 3.0]
+    return
+end
   
 function test_objective_norm_univariate()
     model = Nonlinear.Model()
@@ -168,4 +170,4 @@ end
 
 end  # module
 
-TestArrayDiff.runtests()
\ No newline at end of file
+TestArrayDiff.runtests() 
\ No newline at end of file

From 5fe52642dbf1a0b9d1ab91da1c66f52650d2871b Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Mon, 29 Sep 2025 10:00:37 +0200
Subject: [PATCH 08/19] Add missing continue

To prevent checking if hcat has same size as children
---
 src/reverse_mode.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/reverse_mode.jl b/src/reverse_mode.jl
index c12e4be..8c53731 100644
--- a/src/reverse_mode.jl
+++ b/src/reverse_mode.jl
@@ -447,6 +447,7 @@ function _reverse_eval(f::_SubexpressionStorage)
                         )
                         @j f.reverse_storage[ix2] = val
                     end
+                    continue
                 elseif op == :norm
                     # Node `k` is scalar, the jacobian w.r.t. the vectorized input
                     # child is a row vector whose entries are stored in `f.partials_storage`

From bc3c92483d57f458631698862332c93f15e9f812 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Beno=C3=AEt=20Legat?= <benoit.legat@gmail.com>
Date: Wed, 1 Oct 2025 11:21:19 +0200
Subject: [PATCH 09/19] Update test/ArrayDiff.jl

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

diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index 093f0c9..0b8529a 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -170,4 +170,4 @@ end
 
 end  # module
 
-TestArrayDiff.runtests() 
\ No newline at end of file
+TestArrayDiff.runtests()
\ No newline at end of file

From b4770da6faf0a296574d45ffe25392b69efecc6d Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Thu, 2 Oct 2025 07:07:10 +0200
Subject: [PATCH 10/19] Fix format

---
 src/reverse_mode.jl | 43 ++++++++++++++++++++++++++++++++++---------
 src/sizes.jl        | 15 ++++++++++-----
 test/ArrayDiff.jl   | 17 ++++++-----------
 3 files changed, 50 insertions(+), 25 deletions(-)

diff --git a/src/reverse_mode.jl b/src/reverse_mode.jl
index 8c53731..6cdc70f 100644
--- a/src/reverse_mode.jl
+++ b/src/reverse_mode.jl
@@ -261,7 +261,13 @@ function _forward_eval(
                 for j in _eachindex(f.sizes, ix2)
                     @j f.partials_storage[ix2] = one(T)
                     val = @j f.forward_storage[ix2]
-                    _setindex!(f.forward_storage, val, f.sizes, k, j + nb_cols1 * col_size)
+                    _setindex!(
+                        f.forward_storage,
+                        val,
+                        f.sizes,
+                        k,
+                        j + nb_cols1 * col_size,
+                    )
                 end
             elseif node.index == 14 # norm 
                 ix = children_arr[children_indices[1]]
@@ -427,23 +433,42 @@ function _reverse_eval(f::_SubexpressionStorage)
                     idx1, idx2 = children_indices
                     ix1 = children_arr[idx1]
                     ix2 = children_arr[idx2]
-                    nb_cols1 = f.sizes.ndims[ix1] <= 1 ? 1 : _size(f.sizes, ix1, 2)
-                    col_size = f.sizes.ndims[ix1] == 0 ? 1 : _size(f.sizes, k, 1)
+                    nb_cols1 =
+                        f.sizes.ndims[ix1] <= 1 ? 1 : _size(f.sizes, ix1, 2)
+                    col_size =
+                        f.sizes.ndims[ix1] == 0 ? 1 : _size(f.sizes, k, 1)
                     for j in _eachindex(f.sizes, ix1)
                         partial = @j f.partials_storage[ix1]
                         val = ifelse(
-                            _getindex(f.reverse_storage, f.sizes, k, j) == 0.0 && !isfinite(partial),
+                            _getindex(f.reverse_storage, f.sizes, k, j) ==
+                            0.0 && !isfinite(partial),
                             _getindex(f.reverse_storage, f.sizes, k, j),
-                            _getindex(f.reverse_storage, f.sizes, k, j) * partial,
+                            _getindex(f.reverse_storage, f.sizes, k, j) *
+                            partial,
                         )
                         @j f.reverse_storage[ix1] = val
                     end
                     for j in _eachindex(f.sizes, ix2)
                         partial = @j f.partials_storage[ix2]
                         val = ifelse(
-                            _getindex(f.reverse_storage, f.sizes, k, j + nb_cols1 * col_size) == 0.0 && !isfinite(partial),
-                            _getindex(f.reverse_storage, f.sizes, k, j + nb_cols1 * col_size),
-                            _getindex(f.reverse_storage, f.sizes, k, j + nb_cols1 * col_size) * partial,
+                            _getindex(
+                                f.reverse_storage,
+                                f.sizes,
+                                k,
+                                j + nb_cols1 * col_size,
+                            ) == 0.0 && !isfinite(partial),
+                            _getindex(
+                                f.reverse_storage,
+                                f.sizes,
+                                k,
+                                j + nb_cols1 * col_size,
+                            ),
+                            _getindex(
+                                f.reverse_storage,
+                                f.sizes,
+                                k,
+                                j + nb_cols1 * col_size,
+                            ) * partial,
                         )
                         @j f.reverse_storage[ix2] = val
                     end
@@ -461,7 +486,7 @@ function _reverse_eval(f::_SubexpressionStorage)
                             rev_parent,
                             rev_parent * partial,
                         )
-                        @j f.reverse_storage[ix] = val  
+                        @j f.reverse_storage[ix] = val
                     end
                     continue
                 end
diff --git a/src/sizes.jl b/src/sizes.jl
index 518b040..002d290 100644
--- a/src/sizes.jl
+++ b/src/sizes.jl
@@ -189,15 +189,20 @@ function _infer_sizes(
             elseif op == :hcat
                 total_cols = 0
                 for c_idx in children_indices
-                    total_cols += sizes.ndims[children_arr[c_idx]] <= 1 ?
-                        1 : _size(sizes, children_arr[c_idx], 2)
+                    total_cols +=
+                        sizes.ndims[children_arr[c_idx]] <= 1 ? 1 :
+                        _size(sizes, children_arr[c_idx], 2)
                 end
                 if sizes.ndims[children_arr[first(children_indices)]] == 0
                     shape = (1, total_cols)
                 elseif sizes.ndims[children_arr[first(children_indices)]] <= 2
-                    shape = ( _size(sizes, children_arr[first(children_indices)], 1), total_cols)
-                else 
-                    child_shape = _size(sizes, children_arr[first(children_indices)])
+                    shape = (
+                        _size(sizes, children_arr[first(children_indices)], 1),
+                        total_cols,
+                    )
+                else
+                    child_shape =
+                        _size(sizes, children_arr[first(children_indices)])
                     shape = (child_shape[1], total_cols, child_shape[3:end]...)
                 end
                 _add_size!(sizes, k, tuple(shape...))
diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index 0b8529a..4bd270b 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -70,10 +70,7 @@ function test_objective_hcat_0dim()
     x2 = MOI.VariableIndex(2)
     x3 = MOI.VariableIndex(3)
     x4 = MOI.VariableIndex(4)
-    Nonlinear.set_objective(
-        model,
-        :(dot([$x1 $x3], [$x2 $x4]))
-    )
+    Nonlinear.set_objective(model, :(dot([$x1 $x3], [$x2 $x4])))
     evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
     MOI.initialize(evaluator, [:Grad])
     sizes = evaluator.backend.objective.expr.sizes
@@ -86,8 +83,7 @@ function test_objective_hcat_0dim()
     x3 = 3.0
     x4 = 4.0
     println(MOI.eval_objective(evaluator, [x1, x2, x3, x4]))
-    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) ==
-          14.0
+    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == 14.0
     g = ones(4)
     MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
     @test g == [2.0, 1.0, 4.0, 3.0]
@@ -102,7 +98,7 @@ function test_objective_hcat_1dim()
     x4 = MOI.VariableIndex(4)
     Nonlinear.set_objective(
         model,
-        :(dot(hcat([$x1], [$x3]), hcat([$x2], [$x4])))
+        :(dot(hcat([$x1], [$x3]), hcat([$x2], [$x4]))),
     )
     evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
     MOI.initialize(evaluator, [:Grad])
@@ -116,14 +112,13 @@ function test_objective_hcat_1dim()
     x3 = 3.0
     x4 = 4.0
     println(MOI.eval_objective(evaluator, [x1, x2, x3, x4]))
-    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) ==
-          14.0
+    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == 14.0
     g = ones(4)
     MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
     @test g == [2.0, 1.0, 4.0, 3.0]
     return
 end
-  
+
 function test_objective_norm_univariate()
     model = Nonlinear.Model()
     x = MOI.VariableIndex(1)
@@ -170,4 +165,4 @@ end
 
 end  # module
 
-TestArrayDiff.runtests()
\ No newline at end of file
+TestArrayDiff.runtests()

From 8e32bf01eaa5a9d1b0fbc3ba1246b6c2f51ecf12 Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Thu, 2 Oct 2025 09:04:00 +0200
Subject: [PATCH 11/19] Remove Revise and simplify hcat

* Remove Revise from toml as mentionned in review
* Remove hcat case 'more than 2 dims', as we do not pretend to handle it right now
---
 src/sizes.jl      | 9 ++++-----
 test/Project.toml | 1 -
 2 files changed, 4 insertions(+), 6 deletions(-)

diff --git a/src/sizes.jl b/src/sizes.jl
index 002d290..47aefbd 100644
--- a/src/sizes.jl
+++ b/src/sizes.jl
@@ -195,15 +195,14 @@ function _infer_sizes(
                 end
                 if sizes.ndims[children_arr[first(children_indices)]] == 0
                     shape = (1, total_cols)
-                elseif sizes.ndims[children_arr[first(children_indices)]] <= 2
+                else
+                    @assert sizes.ndims[children_arr[first(
+                        children_indices,
+                    )]] <= 2 "Hcat with ndims > 2 is not supported yet"
                     shape = (
                         _size(sizes, children_arr[first(children_indices)], 1),
                         total_cols,
                     )
-                else
-                    child_shape =
-                        _size(sizes, children_arr[first(children_indices)])
-                    shape = (child_shape[1], total_cols, child_shape[3:end]...)
                 end
                 _add_size!(sizes, k, tuple(shape...))
             elseif op == :*
diff --git a/test/Project.toml b/test/Project.toml
index 0e679fe..7ed96af 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,6 +2,5 @@
 ArrayDiff = "c45fa1ca-6901-44ac-ae5b-5513a4852d50"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
-Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

From 8983f4df7051509911408d5d97fa8f28aa987308 Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Thu, 2 Oct 2025 10:04:33 +0200
Subject: [PATCH 12/19] Delete personal comments

---
 src/reverse_mode.jl | 12 ------------
 1 file changed, 12 deletions(-)

diff --git a/src/reverse_mode.jl b/src/reverse_mode.jl
index 6cdc70f..16ef670 100644
--- a/src/reverse_mode.jl
+++ b/src/reverse_mode.jl
@@ -361,18 +361,6 @@ function _forward_eval(
             f.partials_storage[rhs] = zero(T)
         end
     end
-    # This function is written assuming that the final output is scalar.
-    # Therefore cannot return the matrix, so I guess I return it's first entry only, 
-    # as long as sum or matx-vect products are not implemented.
-
-    #println("Last node ", f.nodes[1].index)
-    #if f.nodes[1].index == 12
-    #    mtx = reshape(
-    #        f.forward_storage[_storage_range(f.sizes, 1)],
-    #        f.sizes.size[1:f.sizes.ndims[1]]...,
-    #    )
-    #    return mtx
-    #end
     return f.forward_storage[1]
 end
 

From 6c01d85ab5cbdc29698cf0768329183a48d385ce Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Sun, 5 Oct 2025 19:22:07 +0200
Subject: [PATCH 13/19] Add vcat and row - evaluation and derivs

---
 src/reverse_mode.jl |  91 +++++++++++++++++++++++++++++++++++++
 src/sizes.jl        |  27 +++++++++++
 test/ArrayDiff.jl   | 106 ++++++++++++++++++++++++++++++++++++++++++--
 3 files changed, 221 insertions(+), 3 deletions(-)

diff --git a/src/reverse_mode.jl b/src/reverse_mode.jl
index 16ef670..b4b8fdb 100644
--- a/src/reverse_mode.jl
+++ b/src/reverse_mode.jl
@@ -269,6 +269,35 @@ function _forward_eval(
                         j + nb_cols1 * col_size,
                     )
                 end
+            elseif node.index == 13 # vcat
+                idx1, idx2 = children_indices
+                ix1 = children_arr[idx1]
+                ix2 = children_arr[idx2]
+                nb_rows1 = f.sizes.ndims[ix1] <= 1 ? 1 : _size(f.sizes, ix1, 1)
+                nb_rows2 = f.sizes.ndims[ix2] <= 1 ? 1 : _size(f.sizes, ix2, 1)
+                nb_rows = nb_rows1 + nb_rows2
+                for j in _eachindex(f.sizes, ix1)
+                    @j f.partials_storage[ix1] = one(T)
+                    val = @j f.forward_storage[ix1]
+                    _setindex!(
+                        f.forward_storage,
+                        val,
+                        f.sizes,
+                        k,
+                        div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1
+                    )
+                end
+                for j in _eachindex(f.sizes, ix2)
+                    @j f.partials_storage[ix2] = one(T)
+                    val = @j f.forward_storage[ix2]
+                    _setindex!(
+                        f.forward_storage,
+                        val,
+                        f.sizes,
+                        k,
+                        div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1 + nb_rows1
+                    )
+                end
             elseif node.index == 14 # norm 
                 ix = children_arr[children_indices[1]]
                 tmp_norm_squared = zero(T)
@@ -285,6 +314,13 @@ function _forward_eval(
                         @j f.partials_storage[ix] = v / @s f.forward_storage[k]
                     end
                 end
+            elseif node.index == 16 # row
+                for j in _eachindex(f.sizes, k)
+                    ix = children_arr[children_indices[j]]
+                    @s f.partials_storage[ix] = one(T)
+                    val = @s f.forward_storage[ix]
+                    @j f.forward_storage[k] = val
+                end
             else # atan, min, max
                 f_input = _UnsafeVectorView(d.jac_storage, N)
                 ∇f = _UnsafeVectorView(d.user_output_buffer, N)
@@ -461,6 +497,53 @@ function _reverse_eval(f::_SubexpressionStorage)
                         @j f.reverse_storage[ix2] = val
                     end
                     continue
+                elseif op == :vcat
+                    idx1, idx2 = children_indices
+                    ix1 = children_arr[idx1]
+                    ix2 = children_arr[idx2]
+                    nb_rows1 =
+                        f.sizes.ndims[ix1] <= 1 ? 1 : _size(f.sizes, ix1, 1)
+                    nb_rows2 = 
+                        f.sizes.ndims[ix2] <= 1 ? 1 : _size(f.sizes, ix2, 1)
+                    nb_rows = nb_rows1 + nb_rows2
+                    row_size =
+                        f.sizes.ndims[ix1] == 0 ? 1 : _size(f.sizes, k, 2)
+                    for j in _eachindex(f.sizes, ix1)
+                        partial = @j f.partials_storage[ix1]
+                        val = ifelse(
+                            _getindex(f.reverse_storage, f.sizes, k, div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1) ==
+                            0.0 && !isfinite(partial),
+                            _getindex(f.reverse_storage, f.sizes, k, div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1),
+                            _getindex(f.reverse_storage, f.sizes, k, div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1) *
+                            partial,
+                        )
+                        @j f.reverse_storage[ix1] = val
+                    end
+                    for j in _eachindex(f.sizes, ix2)
+                        partial = @j f.partials_storage[ix2]
+                        val = ifelse(
+                            _getindex(
+                                f.reverse_storage,
+                                f.sizes,
+                                k,
+                                div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1 + nb_rows1,
+                            ) == 0.0 && !isfinite(partial),
+                            _getindex(
+                                f.reverse_storage,
+                                f.sizes,
+                                k,
+                                div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1 + nb_rows1,
+                            ),
+                            _getindex(
+                                f.reverse_storage,
+                                f.sizes,
+                                k,
+                                div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1 + nb_rows1,
+                            ) * partial,
+                        )
+                        @j f.reverse_storage[ix2] = val
+                    end
+                    continue
                 elseif op == :norm
                     # Node `k` is scalar, the jacobian w.r.t. the vectorized input
                     # child is a row vector whose entries are stored in `f.partials_storage`
@@ -477,6 +560,14 @@ function _reverse_eval(f::_SubexpressionStorage)
                         @j f.reverse_storage[ix] = val
                     end
                     continue
+                elseif op == :row
+                    for j in _eachindex(f.sizes, k)
+                        ix = children_arr[children_indices[j]]
+                        rev_parent_j = @j f.reverse_storage[k]
+                        # partial is 1 so we can ignore it
+                        @s f.reverse_storage[ix] = rev_parent_j
+                    end
+                    continue
                 end
             end
         elseif node.type != MOI.Nonlinear.NODE_CALL_UNIVARIATE
diff --git a/src/sizes.jl b/src/sizes.jl
index 47aefbd..60fa7fc 100644
--- a/src/sizes.jl
+++ b/src/sizes.jl
@@ -179,6 +179,14 @@ function _infer_sizes(
                     op,
                 )
                 _add_size!(sizes, k, (N,))
+            elseif op == :row
+                _assert_scalar_children(
+                    sizes,
+                    children_arr,
+                    children_indices,
+                    op,
+                )
+                _add_size!(sizes, k, (1,N))
             elseif op == :dot
                 # TODO assert all arguments have same size
             elseif op == :norm
@@ -205,6 +213,25 @@ function _infer_sizes(
                     )
                 end
                 _add_size!(sizes, k, tuple(shape...))
+            elseif op == :vcat
+                total_rows = 0
+                for c_idx in children_indices
+                    total_rows +=
+                        sizes.ndims[children_arr[c_idx]] <= 1 ? 1 :
+                        _size(sizes, children_arr[c_idx], 1)
+                end
+                if sizes.ndims[children_arr[first(children_indices)]] == 0
+                    shape = (total_rows, 1)
+                else
+                    @assert sizes.ndims[children_arr[first(
+                        children_indices,
+                    )]] <= 2 "Hcat with ndims > 2 is not supported yet"
+                    shape = (
+                        total_rows,
+                        _size(sizes, children_arr[first(children_indices)], 2),
+                    )
+                end
+                _add_size!(sizes, k, tuple(shape...))
             elseif op == :*
                 # TODO assert compatible sizes and all ndims should be 0 or 2
                 first_matrix = findfirst(children_indices) do i
diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index 4bd270b..04ce2aa 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -82,7 +82,6 @@ function test_objective_hcat_0dim()
     x2 = 2.0
     x3 = 3.0
     x4 = 4.0
-    println(MOI.eval_objective(evaluator, [x1, x2, x3, x4]))
     @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == 14.0
     g = ones(4)
     MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
@@ -111,7 +110,6 @@ function test_objective_hcat_1dim()
     x2 = 2.0
     x3 = 3.0
     x4 = 4.0
-    println(MOI.eval_objective(evaluator, [x1, x2, x3, x4]))
     @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == 14.0
     g = ones(4)
     MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
@@ -163,6 +161,108 @@ function test_objective_norm_bivariate()
     return
 end
 
+function test_objective_norm_of_row()
+    model = Nonlinear.Model()
+    x1 = MOI.VariableIndex(1)
+    x2 = MOI.VariableIndex(2)
+    Nonlinear.set_objective(model, :(norm(row($x1, $x2))))
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 2, 0, 0]
+    @test sizes.size_offset == [0, 0, 0, 0]
+    @test sizes.size == [1, 2]
+    @test sizes.storage_offset == [0, 1, 3, 4, 5]
+    x1 = 1.0
+    x2 = 2.0
+    @test MOI.eval_objective(evaluator, [x1, x2]) == sqrt(5.0)
+    g = ones(2)
+    MOI.eval_objective_gradient(evaluator, g, [x1, x2])
+    @test g == [1.0 / sqrt(5.0), 2.0 / sqrt(5.0)]
+    return
+end
+
+function test_objective_vcat_0dim()
+    model = Nonlinear.Model()
+    x1 = MOI.VariableIndex(1)
+    x2 = MOI.VariableIndex(2)
+    x3 = MOI.VariableIndex(3)
+    x4 = MOI.VariableIndex(4)
+    Nonlinear.set_objective(model, :(norm(vcat($x1, $x3))))
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 2, 0, 0]
+    @test sizes.size_offset == [0, 0, 0, 0]
+    @test sizes.size == [2, 1]
+    @test sizes.storage_offset == [0, 1, 3, 4, 5]
+    x1 = 1.0
+    x2 = 2.0
+    x3 = 3.0
+    x4 = 4.0
+    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(10.0)
+    g = ones(4)
+    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
+    @test g == [1.0 / sqrt(10.0), 0.0, 3.0 / sqrt(10.0), 0.0]
+    return
+end  
+
+function test_objective_vcat_1dim()
+    model = Nonlinear.Model()
+    x1 = MOI.VariableIndex(1)
+    x2 = MOI.VariableIndex(2)
+    x3 = MOI.VariableIndex(3)
+    x4 = MOI.VariableIndex(4)
+    Nonlinear.set_objective(
+        model,
+        :(norm(vcat([$x1 $x3], [$x2 $x4])))
+    )
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 2, 2, 0, 0, 2, 0, 0]
+    @test sizes.size_offset == [0, 4, 2, 0, 0, 0, 0, 0]
+    @test sizes.size == [1, 2, 1, 2, 2, 2]
+    @test sizes.storage_offset == [0, 1, 5, 7, 8, 9, 11, 12, 13]
+    x1 = 1.0
+    x2 = 2.0
+    x3 = 3.0
+    x4 = 4.0
+    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(30.0)
+    g = ones(4)
+    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
+    @test g == [1.0 / sqrt(30.0), 2.0 / sqrt(30.0), 3.0 / sqrt(30.0), 4.0 / sqrt(30.0)]
+    return
+end
+
+function test_objective_norm_of_matrix()
+    model = Nonlinear.Model()
+    x1 = MOI.VariableIndex(1)
+    x2 = MOI.VariableIndex(2)
+    x3 = MOI.VariableIndex(3)
+    x4 = MOI.VariableIndex(4)
+    Nonlinear.set_objective(
+        model,
+        :(norm([$x1 $x2; $x3 $x4])),
+    )
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 2, 2, 0, 0, 2, 0, 0]
+    @test sizes.size_offset == [0, 4, 2, 0, 0, 0, 0, 0]
+    @test sizes.size == [1, 2, 1, 2, 2, 2]
+    @test sizes.storage_offset == [0, 1, 5, 7, 8, 9, 11, 12, 13]
+    x1 = 1.0
+    x2 = 2.0
+    x3 = 3.0
+    x4 = 4.0
+    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(30.0)
+    g = ones(4)
+    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
+    @test g == [1.0 / sqrt(30.0), 2.0 / sqrt(30.0), 3.0 / sqrt(30.0), 4.0 / sqrt(30.0)]
+    return
+end
+
 end  # module
 
-TestArrayDiff.runtests()
+TestArrayDiff.runtests()
\ No newline at end of file

From 98ae6553fb83a7d3f15de7efeb346d141929ed38 Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Thu, 13 Nov 2025 15:59:13 +0100
Subject: [PATCH 14/19] Add test | scalar mult to be fixed

* Add test Frobenius norm of matrix
* !! Multiplication by scalar nb still to be fixed, SL
---
 test/ArrayDiff.jl | 56 +++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 56 insertions(+)

diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index 04ce2aa..7dc27d6 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -261,8 +261,64 @@ function test_objective_norm_of_matrix()
     MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
     @test g == [1.0 / sqrt(30.0), 2.0 / sqrt(30.0), 3.0 / sqrt(30.0), 4.0 / sqrt(30.0)]
     return
+end 
+
+function test_objective_norm_of_matrix_with_sum()
+    model = Nonlinear.Model()
+    x1 = MOI.VariableIndex(1)
+    x2 = MOI.VariableIndex(2)
+    x3 = MOI.VariableIndex(3)
+    x4 = MOI.VariableIndex(4)
+    Nonlinear.set_objective(
+        model,
+        :(norm([$x1 $x2; $x3 $x4] - [1 1; 1 1])),
+    )
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 2, 2, 2, 0, 0, 2, 0, 0, 2, 2, 0, 0, 2, 0, 0]
+    @test sizes.size_offset == [0, 12, 10, 8, 0, 0, 6, 0, 0, 4, 2, 0, 0, 0, 0, 0]
+    @test sizes.size == [1, 2, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 2]
+    @test sizes.storage_offset == [0, 1, 5, 9, 11, 12, 13, 15, 16, 17, 21, 23, 24, 25, 27, 28, 29]
+    x1 = 1.0
+    x2 = 2.0
+    x3 = 3.0
+    x4 = 4.0
+    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(14.0)
+    g = ones(4)
+    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
+    @test g == [0.0, 1.0 / sqrt(14.0), 2.0 / sqrt(14.0), 3.0 / sqrt(14.0)]
+    return
 end
 
+#function test_objective_norm_of_matrix_with_scalar_mult()
+#    model = Nonlinear.Model()
+#    x1 = MOI.VariableIndex(1)
+#    x2 = MOI.VariableIndex(2)
+#    x3 = MOI.VariableIndex(3)
+#    x4 = MOI.VariableIndex(4)
+#    Nonlinear.set_objective(
+#        model,
+#        :(norm(2 * [$x1 $x2; $x3 $x4] - [1 2; 3 4])),
+#    )
+#    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
+#    MOI.initialize(evaluator, [:Grad])
+#    sizes = evaluator.backend.objective.expr.sizes
+#    @test sizes.ndims == [0, 2, 2, 2, 0, 0, 2, 0, 0, 2, 2, 0, 0, 2, 0, 0]
+#    @test sizes.size_offset == [0, 12, 10, 8, 0, 0, 6, 0, 0, 4, 2, 0, 0, 0, 0, 0]
+#    @test sizes.size == [1, 2, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 2]
+#    @test sizes.storage_offset == [0, 1, 5, 9, 11, 12, 13, 15, 16, 17, 21, 23, 24, 25, 27, 28, 29]
+#    x1 = 2.0
+#    x2 = 4.0
+#    x3 = 6.0
+#    x4 = 8.0
+#    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(30.0)
+#    g = ones(4)
+#    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
+#    @test g == [1.0 / sqrt(30.0), 2.0 / sqrt(30.0), 3.0 / sqrt(30.0), 4.0 / sqrt(30.0)]
+#    return
+#end
+
 end  # module
 
 TestArrayDiff.runtests()
\ No newline at end of file

From 750c49955269ac3b23fa0c8f8d46f3ca84da4143 Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Tue, 9 Dec 2025 11:58:06 +0100
Subject: [PATCH 15/19] Format and rename some tests

Summary: tests cover dot and norm, between vectors or matrices constructed with vect, hcat, vcat or row
---
 src/reverse_mode.jl |  53 ++++++++++++++++++-----
 src/sizes.jl        |  28 ++++++++----
 test/ArrayDiff.jl   | 101 +++++++++++++++++++++-----------------------
 3 files changed, 110 insertions(+), 72 deletions(-)

diff --git a/src/reverse_mode.jl b/src/reverse_mode.jl
index b4b8fdb..3456ee7 100644
--- a/src/reverse_mode.jl
+++ b/src/reverse_mode.jl
@@ -284,7 +284,7 @@ function _forward_eval(
                         val,
                         f.sizes,
                         k,
-                        div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1
+                        div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1,
                     )
                 end
                 for j in _eachindex(f.sizes, ix2)
@@ -295,7 +295,10 @@ function _forward_eval(
                         val,
                         f.sizes,
                         k,
-                        div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1 + nb_rows1
+                        div(j-1, nb_rows1) * nb_rows +
+                        1 +
+                        (j-1) % nb_rows1 +
+                        nb_rows1,
                     )
                 end
             elseif node.index == 14 # norm 
@@ -503,7 +506,7 @@ function _reverse_eval(f::_SubexpressionStorage)
                     ix2 = children_arr[idx2]
                     nb_rows1 =
                         f.sizes.ndims[ix1] <= 1 ? 1 : _size(f.sizes, ix1, 1)
-                    nb_rows2 = 
+                    nb_rows2 =
                         f.sizes.ndims[ix2] <= 1 ? 1 : _size(f.sizes, ix2, 1)
                     nb_rows = nb_rows1 + nb_rows2
                     row_size =
@@ -511,11 +514,30 @@ function _reverse_eval(f::_SubexpressionStorage)
                     for j in _eachindex(f.sizes, ix1)
                         partial = @j f.partials_storage[ix1]
                         val = ifelse(
-                            _getindex(f.reverse_storage, f.sizes, k, div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1) ==
-                            0.0 && !isfinite(partial),
-                            _getindex(f.reverse_storage, f.sizes, k, div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1),
-                            _getindex(f.reverse_storage, f.sizes, k, div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1) *
-                            partial,
+                            _getindex(
+                                f.reverse_storage,
+                                f.sizes,
+                                k,
+                                div(j-1, nb_rows1) * nb_rows +
+                                1 +
+                                (j-1) % nb_rows1,
+                            ) == 0.0 && !isfinite(partial),
+                            _getindex(
+                                f.reverse_storage,
+                                f.sizes,
+                                k,
+                                div(j-1, nb_rows1) * nb_rows +
+                                1 +
+                                (j-1) % nb_rows1,
+                            ),
+                            _getindex(
+                                f.reverse_storage,
+                                f.sizes,
+                                k,
+                                div(j-1, nb_rows1) * nb_rows +
+                                1 +
+                                (j-1) % nb_rows1,
+                            ) * partial,
                         )
                         @j f.reverse_storage[ix1] = val
                     end
@@ -526,19 +548,28 @@ function _reverse_eval(f::_SubexpressionStorage)
                                 f.reverse_storage,
                                 f.sizes,
                                 k,
-                                div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1 + nb_rows1,
+                                div(j-1, nb_rows1) * nb_rows +
+                                1 +
+                                (j-1) % nb_rows1 +
+                                nb_rows1,
                             ) == 0.0 && !isfinite(partial),
                             _getindex(
                                 f.reverse_storage,
                                 f.sizes,
                                 k,
-                                div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1 + nb_rows1,
+                                div(j-1, nb_rows1) * nb_rows +
+                                1 +
+                                (j-1) % nb_rows1 +
+                                nb_rows1,
                             ),
                             _getindex(
                                 f.reverse_storage,
                                 f.sizes,
                                 k,
-                                div(j-1, nb_rows1) * nb_rows + 1 + (j-1) % nb_rows1 + nb_rows1,
+                                div(j-1, nb_rows1) * nb_rows +
+                                1 +
+                                (j-1) % nb_rows1 +
+                                nb_rows1,
                             ) * partial,
                         )
                         @j f.reverse_storage[ix2] = val
diff --git a/src/sizes.jl b/src/sizes.jl
index 60fa7fc..a747656 100644
--- a/src/sizes.jl
+++ b/src/sizes.jl
@@ -186,7 +186,7 @@ function _infer_sizes(
                     children_indices,
                     op,
                 )
-                _add_size!(sizes, k, (1,N))
+                _add_size!(sizes, k, (1, N))
             elseif op == :dot
                 # TODO assert all arguments have same size
             elseif op == :norm
@@ -241,14 +241,24 @@ function _infer_sizes(
                     last_matrix = findfirst(children_indices) do i
                         return !iszero(sizes.ndims[children_arr[i]])
                     end
-                    _add_size!(
-                        sizes,
-                        k,
-                        (
-                            _size(sizes, first_matrix, 1),
-                            _size(sizes, last_matrix, sizes.ndims[last_matrix]),
-                        ),
-                    )
+                    if sizes.ndims[last_matrix] == 0 ||
+                       sizes.ndims[first_matrix] == 0
+                        _add_size!(sizes, k, (1, 1))
+                        continue
+                    else
+                        _add_size!(
+                            sizes,
+                            k,
+                            (
+                                _size(sizes, first_matrix, 1),
+                                _size(
+                                    sizes,
+                                    last_matrix,
+                                    sizes.ndims[last_matrix],
+                                ),
+                            ),
+                        )
+                    end
                 end
             elseif op == :^ || op == :/
                 @assert N == 2
diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index 7dc27d6..732ed0a 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -64,7 +64,7 @@ function test_objective_dot_bivariate()
     return
 end
 
-function test_objective_hcat_0dim()
+function test_objective_hcat_scalars()
     model = Nonlinear.Model()
     x1 = MOI.VariableIndex(1)
     x2 = MOI.VariableIndex(2)
@@ -89,7 +89,7 @@ function test_objective_hcat_0dim()
     return
 end
 
-function test_objective_hcat_1dim()
+function test_objective_hcat_vectors()
     model = Nonlinear.Model()
     x1 = MOI.VariableIndex(1)
     x2 = MOI.VariableIndex(2)
@@ -117,6 +117,28 @@ function test_objective_hcat_1dim()
     return
 end
 
+function test_objective_dot_bivariate_on_rows()
+    model = Nonlinear.Model()
+    x = MOI.VariableIndex(1)
+    y = MOI.VariableIndex(2)
+    Nonlinear.set_objective(model, :(dot([$x $y] - [1 2], -[1 2] + [$x $y])))
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x, y])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 2, 2, 0, 0, 2, 0, 0, 2, 2, 2, 0, 0, 2, 0, 0]
+    @test sizes.size_offset ==
+          [0, 12, 10, 0, 0, 8, 0, 0, 6, 4, 2, 0, 0, 0, 0, 0]
+    @test sizes.size == [1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2]
+    @test sizes.storage_offset ==
+          [0, 1, 3, 5, 6, 7, 9, 10, 11, 13, 15, 17, 18, 19, 21, 22, 23]
+    x = [5, -1]
+    @test MOI.eval_objective(evaluator, x) ≈ 25
+    g = ones(2)
+    MOI.eval_objective_gradient(evaluator, g, x)
+    @test g == 2(x - [1, 2])
+    return
+end
+
 function test_objective_norm_univariate()
     model = Nonlinear.Model()
     x = MOI.VariableIndex(1)
@@ -165,7 +187,7 @@ function test_objective_norm_of_row()
     model = Nonlinear.Model()
     x1 = MOI.VariableIndex(1)
     x2 = MOI.VariableIndex(2)
-    Nonlinear.set_objective(model, :(norm(row($x1, $x2))))
+    Nonlinear.set_objective(model, :(norm([$x1 $x2])))
     evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2])
     MOI.initialize(evaluator, [:Grad])
     sizes = evaluator.backend.objective.expr.sizes
@@ -182,7 +204,7 @@ function test_objective_norm_of_row()
     return
 end
 
-function test_objective_vcat_0dim()
+function test_objective_norm_of_vcat_vector()
     model = Nonlinear.Model()
     x1 = MOI.VariableIndex(1)
     x2 = MOI.VariableIndex(2)
@@ -205,18 +227,15 @@ function test_objective_vcat_0dim()
     MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
     @test g == [1.0 / sqrt(10.0), 0.0, 3.0 / sqrt(10.0), 0.0]
     return
-end  
+end
 
-function test_objective_vcat_1dim()
+function test_objective_norm_of_vcat_matrix()
     model = Nonlinear.Model()
     x1 = MOI.VariableIndex(1)
     x2 = MOI.VariableIndex(2)
     x3 = MOI.VariableIndex(3)
     x4 = MOI.VariableIndex(4)
-    Nonlinear.set_objective(
-        model,
-        :(norm(vcat([$x1 $x3], [$x2 $x4])))
-    )
+    Nonlinear.set_objective(model, :(norm(vcat([$x1 $x3], [$x2 $x4]))))
     evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
     MOI.initialize(evaluator, [:Grad])
     sizes = evaluator.backend.objective.expr.sizes
@@ -231,7 +250,12 @@ function test_objective_vcat_1dim()
     @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(30.0)
     g = ones(4)
     MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
-    @test g == [1.0 / sqrt(30.0), 2.0 / sqrt(30.0), 3.0 / sqrt(30.0), 4.0 / sqrt(30.0)]
+    @test g == [
+        1.0 / sqrt(30.0),
+        2.0 / sqrt(30.0),
+        3.0 / sqrt(30.0),
+        4.0 / sqrt(30.0),
+    ]
     return
 end
 
@@ -241,10 +265,7 @@ function test_objective_norm_of_matrix()
     x2 = MOI.VariableIndex(2)
     x3 = MOI.VariableIndex(3)
     x4 = MOI.VariableIndex(4)
-    Nonlinear.set_objective(
-        model,
-        :(norm([$x1 $x2; $x3 $x4])),
-    )
+    Nonlinear.set_objective(model, :(norm([$x1 $x2; $x3 $x4])))
     evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
     MOI.initialize(evaluator, [:Grad])
     sizes = evaluator.backend.objective.expr.sizes
@@ -259,9 +280,14 @@ function test_objective_norm_of_matrix()
     @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(30.0)
     g = ones(4)
     MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
-    @test g == [1.0 / sqrt(30.0), 2.0 / sqrt(30.0), 3.0 / sqrt(30.0), 4.0 / sqrt(30.0)]
+    @test g == [
+        1.0 / sqrt(30.0),
+        2.0 / sqrt(30.0),
+        3.0 / sqrt(30.0),
+        4.0 / sqrt(30.0),
+    ]
     return
-end 
+end
 
 function test_objective_norm_of_matrix_with_sum()
     model = Nonlinear.Model()
@@ -269,17 +295,16 @@ function test_objective_norm_of_matrix_with_sum()
     x2 = MOI.VariableIndex(2)
     x3 = MOI.VariableIndex(3)
     x4 = MOI.VariableIndex(4)
-    Nonlinear.set_objective(
-        model,
-        :(norm([$x1 $x2; $x3 $x4] - [1 1; 1 1])),
-    )
+    Nonlinear.set_objective(model, :(norm([$x1 $x2; $x3 $x4] - [1 1; 1 1])))
     evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
     MOI.initialize(evaluator, [:Grad])
     sizes = evaluator.backend.objective.expr.sizes
     @test sizes.ndims == [0, 2, 2, 2, 0, 0, 2, 0, 0, 2, 2, 0, 0, 2, 0, 0]
-    @test sizes.size_offset == [0, 12, 10, 8, 0, 0, 6, 0, 0, 4, 2, 0, 0, 0, 0, 0]
+    @test sizes.size_offset ==
+          [0, 12, 10, 8, 0, 0, 6, 0, 0, 4, 2, 0, 0, 0, 0, 0]
     @test sizes.size == [1, 2, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 2]
-    @test sizes.storage_offset == [0, 1, 5, 9, 11, 12, 13, 15, 16, 17, 21, 23, 24, 25, 27, 28, 29]
+    @test sizes.storage_offset ==
+          [0, 1, 5, 9, 11, 12, 13, 15, 16, 17, 21, 23, 24, 25, 27, 28, 29]
     x1 = 1.0
     x2 = 2.0
     x3 = 3.0
@@ -291,34 +316,6 @@ function test_objective_norm_of_matrix_with_sum()
     return
 end
 
-#function test_objective_norm_of_matrix_with_scalar_mult()
-#    model = Nonlinear.Model()
-#    x1 = MOI.VariableIndex(1)
-#    x2 = MOI.VariableIndex(2)
-#    x3 = MOI.VariableIndex(3)
-#    x4 = MOI.VariableIndex(4)
-#    Nonlinear.set_objective(
-#        model,
-#        :(norm(2 * [$x1 $x2; $x3 $x4] - [1 2; 3 4])),
-#    )
-#    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
-#    MOI.initialize(evaluator, [:Grad])
-#    sizes = evaluator.backend.objective.expr.sizes
-#    @test sizes.ndims == [0, 2, 2, 2, 0, 0, 2, 0, 0, 2, 2, 0, 0, 2, 0, 0]
-#    @test sizes.size_offset == [0, 12, 10, 8, 0, 0, 6, 0, 0, 4, 2, 0, 0, 0, 0, 0]
-#    @test sizes.size == [1, 2, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 2]
-#    @test sizes.storage_offset == [0, 1, 5, 9, 11, 12, 13, 15, 16, 17, 21, 23, 24, 25, 27, 28, 29]
-#    x1 = 2.0
-#    x2 = 4.0
-#    x3 = 6.0
-#    x4 = 8.0
-#    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(30.0)
-#    g = ones(4)
-#    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
-#    @test g == [1.0 / sqrt(30.0), 2.0 / sqrt(30.0), 3.0 / sqrt(30.0), 4.0 / sqrt(30.0)]
-#    return
-#end
-
 end  # module
 
-TestArrayDiff.runtests()
\ No newline at end of file
+TestArrayDiff.runtests()

From 1b6d9c299d85eb2c0d6c452512017e90387ec63f Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Tue, 9 Dec 2025 15:03:01 +0100
Subject: [PATCH 16/19] Try to pass github tests before going further

---
 test/ArrayDiff.jl | 112 +++++++++++++++++++++++-----------------------
 1 file changed, 56 insertions(+), 56 deletions(-)

diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index 732ed0a..a775ea3 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -259,62 +259,62 @@ function test_objective_norm_of_vcat_matrix()
     return
 end
 
-function test_objective_norm_of_matrix()
-    model = Nonlinear.Model()
-    x1 = MOI.VariableIndex(1)
-    x2 = MOI.VariableIndex(2)
-    x3 = MOI.VariableIndex(3)
-    x4 = MOI.VariableIndex(4)
-    Nonlinear.set_objective(model, :(norm([$x1 $x2; $x3 $x4])))
-    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
-    MOI.initialize(evaluator, [:Grad])
-    sizes = evaluator.backend.objective.expr.sizes
-    @test sizes.ndims == [0, 2, 2, 0, 0, 2, 0, 0]
-    @test sizes.size_offset == [0, 4, 2, 0, 0, 0, 0, 0]
-    @test sizes.size == [1, 2, 1, 2, 2, 2]
-    @test sizes.storage_offset == [0, 1, 5, 7, 8, 9, 11, 12, 13]
-    x1 = 1.0
-    x2 = 2.0
-    x3 = 3.0
-    x4 = 4.0
-    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(30.0)
-    g = ones(4)
-    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
-    @test g == [
-        1.0 / sqrt(30.0),
-        2.0 / sqrt(30.0),
-        3.0 / sqrt(30.0),
-        4.0 / sqrt(30.0),
-    ]
-    return
-end
-
-function test_objective_norm_of_matrix_with_sum()
-    model = Nonlinear.Model()
-    x1 = MOI.VariableIndex(1)
-    x2 = MOI.VariableIndex(2)
-    x3 = MOI.VariableIndex(3)
-    x4 = MOI.VariableIndex(4)
-    Nonlinear.set_objective(model, :(norm([$x1 $x2; $x3 $x4] - [1 1; 1 1])))
-    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
-    MOI.initialize(evaluator, [:Grad])
-    sizes = evaluator.backend.objective.expr.sizes
-    @test sizes.ndims == [0, 2, 2, 2, 0, 0, 2, 0, 0, 2, 2, 0, 0, 2, 0, 0]
-    @test sizes.size_offset ==
-          [0, 12, 10, 8, 0, 0, 6, 0, 0, 4, 2, 0, 0, 0, 0, 0]
-    @test sizes.size == [1, 2, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 2]
-    @test sizes.storage_offset ==
-          [0, 1, 5, 9, 11, 12, 13, 15, 16, 17, 21, 23, 24, 25, 27, 28, 29]
-    x1 = 1.0
-    x2 = 2.0
-    x3 = 3.0
-    x4 = 4.0
-    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(14.0)
-    g = ones(4)
-    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
-    @test g == [0.0, 1.0 / sqrt(14.0), 2.0 / sqrt(14.0), 3.0 / sqrt(14.0)]
-    return
-end
+#function test_objective_norm_of_matrix()
+#    model = Nonlinear.Model()
+#    x1 = MOI.VariableIndex(1)
+#    x2 = MOI.VariableIndex(2)
+#    x3 = MOI.VariableIndex(3)
+#    x4 = MOI.VariableIndex(4)
+#    Nonlinear.set_objective(model, :(norm([$x1 $x2; $x3 $x4])))
+#    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
+#    MOI.initialize(evaluator, [:Grad])
+#    sizes = evaluator.backend.objective.expr.sizes
+#    @test sizes.ndims == [0, 2, 2, 0, 0, 2, 0, 0]
+#    @test sizes.size_offset == [0, 4, 2, 0, 0, 0, 0, 0]
+#    @test sizes.size == [1, 2, 1, 2, 2, 2]
+#    @test sizes.storage_offset == [0, 1, 5, 7, 8, 9, 11, 12, 13]
+#    x1 = 1.0
+#    x2 = 2.0
+#    x3 = 3.0
+#    x4 = 4.0
+#    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(30.0)
+#    g = ones(4)
+#    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
+#    @test g == [
+#        1.0 / sqrt(30.0),
+#        2.0 / sqrt(30.0),
+#        3.0 / sqrt(30.0),
+#        4.0 / sqrt(30.0),
+#    ]
+#    return
+#end
+#
+#function test_objective_norm_of_matrix_with_sum()
+#    model = Nonlinear.Model()
+#    x1 = MOI.VariableIndex(1)
+#    x2 = MOI.VariableIndex(2)
+#    x3 = MOI.VariableIndex(3)
+#    x4 = MOI.VariableIndex(4)
+#    Nonlinear.set_objective(model, :(norm([$x1 $x2; $x3 $x4] - [1 1; 1 1])))
+#    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
+#    MOI.initialize(evaluator, [:Grad])
+#    sizes = evaluator.backend.objective.expr.sizes
+#    @test sizes.ndims == [0, 2, 2, 2, 0, 0, 2, 0, 0, 2, 2, 0, 0, 2, 0, 0]
+#    @test sizes.size_offset ==
+#          [0, 12, 10, 8, 0, 0, 6, 0, 0, 4, 2, 0, 0, 0, 0, 0]
+#    @test sizes.size == [1, 2, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 2]
+#    @test sizes.storage_offset ==
+#          [0, 1, 5, 9, 11, 12, 13, 15, 16, 17, 21, 23, 24, 25, 27, 28, 29]
+#    x1 = 1.0
+#    x2 = 2.0
+#    x3 = 3.0
+#    x4 = 4.0
+#    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(14.0)
+#    g = ones(4)
+#    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
+#    @test g == [0.0, 1.0 / sqrt(14.0), 2.0 / sqrt(14.0), 3.0 / sqrt(14.0)]
+#    return
+#end
 
 end  # module
 

From db533ae072635a715e286ccb7f353b4b43b025d0 Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Tue, 9 Dec 2025 15:20:26 +0100
Subject: [PATCH 17/19] Small fix

---
 src/Coloring/Coloring.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Coloring/Coloring.jl b/src/Coloring/Coloring.jl
index 818b6e4..1ef8dc5 100644
--- a/src/Coloring/Coloring.jl
+++ b/src/Coloring/Coloring.jl
@@ -206,7 +206,7 @@ function acyclic_coloring(g::UndirectedGraph)
     firstVisitToTree = fill(_Edge(0, 0, 0), _num_edges(g))
     color = fill(0, _num_vertices(g))
     # disjoint set forest of edges in the graph
-    S = DataStructures.IntDisjointSets(_num_edges(g))
+    S = DataStructures.IntDisjointSet{Int}(_num_edges(g))
     @inbounds for v in 1:_num_vertices(g)
         n_neighbor = _num_neighbors(v, g)
         start_neighbor = _start_neighbors(v, g)

From 4644804fc87cd8d60975eb2aa95672649dbd6635 Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Tue, 9 Dec 2025 15:26:06 +0100
Subject: [PATCH 18/19] Add test to improve code coverage

---
 test/ArrayDiff.jl | 23 ++++++++++++++++++++++-
 1 file changed, 22 insertions(+), 1 deletion(-)

diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index a775ea3..f16f781 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -183,7 +183,7 @@ function test_objective_norm_bivariate()
     return
 end
 
-function test_objective_norm_of_row()
+function test_objective_norm_of_row_vector()
     model = Nonlinear.Model()
     x1 = MOI.VariableIndex(1)
     x2 = MOI.VariableIndex(2)
@@ -259,6 +259,27 @@ function test_objective_norm_of_vcat_matrix()
     return
 end
 
+function test_objective_norm_of_row()
+    model = Nonlinear.Model()
+    x1 = MOI.VariableIndex(1)
+    x2 = MOI.VariableIndex(2)
+    Nonlinear.set_objective(model, :(norm(row($x1, $x2))))
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 2, 0, 0]
+    @test sizes.size_offset == [0, 0, 0, 0]
+    @test sizes.size == [1, 2]
+    @test sizes.storage_offset == [0, 1, 3, 4, 5]
+    x1 = 1.0
+    x2 = 2.0
+    @test MOI.eval_objective(evaluator, [x1, x2]) == sqrt(5.0)
+    g = ones(2)
+    MOI.eval_objective_gradient(evaluator, g, [x1, x2])
+    @test g == [1.0 / sqrt(5.0), 2.0 / sqrt(5.0)]
+    return
+end
+
 #function test_objective_norm_of_matrix()
 #    model = Nonlinear.Model()
 #    x1 = MOI.VariableIndex(1)

From f0fad53f93f5cf149f11cf5d5fb3611ffea1c1ba Mon Sep 17 00:00:00 2001
From: Sophie L <sophie.lequeu@student.uclouvain.be>
Date: Tue, 9 Dec 2025 17:24:06 +0100
Subject: [PATCH 19/19] Add tests

---
 test/ArrayDiff.jl | 131 ++++++++++++++++++++++++++--------------------
 1 file changed, 75 insertions(+), 56 deletions(-)

diff --git a/test/ArrayDiff.jl b/test/ArrayDiff.jl
index f16f781..5eab3cc 100644
--- a/test/ArrayDiff.jl
+++ b/test/ArrayDiff.jl
@@ -40,6 +40,25 @@ function test_objective_dot_univariate()
     return
 end
 
+function test_objective_dot_univariate_and_scalar_mult()
+    model = Nonlinear.Model()
+    x = MOI.VariableIndex(1)
+    Nonlinear.set_objective(model, :(2*(dot([$x], [$x]))))
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 0, 0, 1, 0, 1, 0]
+    @test sizes.size_offset == [0, 0, 0, 1, 0, 0, 0]
+    @test sizes.size == [1, 1]
+    @test sizes.storage_offset == [0, 1, 2, 3, 4, 5, 6, 7]
+    x = [1.2]
+    @test MOI.eval_objective(evaluator, x) == 2*x[1]^2
+    g = ones(1)
+    MOI.eval_objective_gradient(evaluator, g, x)
+    @test g[1] == 4x[1]
+    return
+end
+
 function test_objective_dot_bivariate()
     model = Nonlinear.Model()
     x = MOI.VariableIndex(1)
@@ -280,62 +299,62 @@ function test_objective_norm_of_row()
     return
 end
 
-#function test_objective_norm_of_matrix()
-#    model = Nonlinear.Model()
-#    x1 = MOI.VariableIndex(1)
-#    x2 = MOI.VariableIndex(2)
-#    x3 = MOI.VariableIndex(3)
-#    x4 = MOI.VariableIndex(4)
-#    Nonlinear.set_objective(model, :(norm([$x1 $x2; $x3 $x4])))
-#    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
-#    MOI.initialize(evaluator, [:Grad])
-#    sizes = evaluator.backend.objective.expr.sizes
-#    @test sizes.ndims == [0, 2, 2, 0, 0, 2, 0, 0]
-#    @test sizes.size_offset == [0, 4, 2, 0, 0, 0, 0, 0]
-#    @test sizes.size == [1, 2, 1, 2, 2, 2]
-#    @test sizes.storage_offset == [0, 1, 5, 7, 8, 9, 11, 12, 13]
-#    x1 = 1.0
-#    x2 = 2.0
-#    x3 = 3.0
-#    x4 = 4.0
-#    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(30.0)
-#    g = ones(4)
-#    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
-#    @test g == [
-#        1.0 / sqrt(30.0),
-#        2.0 / sqrt(30.0),
-#        3.0 / sqrt(30.0),
-#        4.0 / sqrt(30.0),
-#    ]
-#    return
-#end
-#
-#function test_objective_norm_of_matrix_with_sum()
-#    model = Nonlinear.Model()
-#    x1 = MOI.VariableIndex(1)
-#    x2 = MOI.VariableIndex(2)
-#    x3 = MOI.VariableIndex(3)
-#    x4 = MOI.VariableIndex(4)
-#    Nonlinear.set_objective(model, :(norm([$x1 $x2; $x3 $x4] - [1 1; 1 1])))
-#    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
-#    MOI.initialize(evaluator, [:Grad])
-#    sizes = evaluator.backend.objective.expr.sizes
-#    @test sizes.ndims == [0, 2, 2, 2, 0, 0, 2, 0, 0, 2, 2, 0, 0, 2, 0, 0]
-#    @test sizes.size_offset ==
-#          [0, 12, 10, 8, 0, 0, 6, 0, 0, 4, 2, 0, 0, 0, 0, 0]
-#    @test sizes.size == [1, 2, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 2]
-#    @test sizes.storage_offset ==
-#          [0, 1, 5, 9, 11, 12, 13, 15, 16, 17, 21, 23, 24, 25, 27, 28, 29]
-#    x1 = 1.0
-#    x2 = 2.0
-#    x3 = 3.0
-#    x4 = 4.0
-#    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(14.0)
-#    g = ones(4)
-#    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
-#    @test g == [0.0, 1.0 / sqrt(14.0), 2.0 / sqrt(14.0), 3.0 / sqrt(14.0)]
-#    return
-#end
+function test_objective_norm_of_matrix()
+    model = Nonlinear.Model()
+    x1 = MOI.VariableIndex(1)
+    x2 = MOI.VariableIndex(2)
+    x3 = MOI.VariableIndex(3)
+    x4 = MOI.VariableIndex(4)
+    Nonlinear.set_objective(model, :(norm([$x1 $x2; $x3 $x4])))
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 2, 2, 0, 0, 2, 0, 0]
+    @test sizes.size_offset == [0, 4, 2, 0, 0, 0, 0, 0]
+    @test sizes.size == [1, 2, 1, 2, 2, 2]
+    @test sizes.storage_offset == [0, 1, 5, 7, 8, 9, 11, 12, 13]
+    x1 = 1.0
+    x2 = 2.0
+    x3 = 3.0
+    x4 = 4.0
+    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(30.0)
+    g = ones(4)
+    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
+    @test g == [
+        1.0 / sqrt(30.0),
+        2.0 / sqrt(30.0),
+        3.0 / sqrt(30.0),
+        4.0 / sqrt(30.0),
+    ]
+    return
+end
+
+function test_objective_norm_of_matrix_with_sum()
+    model = Nonlinear.Model()
+    x1 = MOI.VariableIndex(1)
+    x2 = MOI.VariableIndex(2)
+    x3 = MOI.VariableIndex(3)
+    x4 = MOI.VariableIndex(4)
+    Nonlinear.set_objective(model, :(norm([$x1 $x2; $x3 $x4] - [1 1; 1 1])))
+    evaluator = Nonlinear.Evaluator(model, ArrayDiff.Mode(), [x1, x2, x3, x4])
+    MOI.initialize(evaluator, [:Grad])
+    sizes = evaluator.backend.objective.expr.sizes
+    @test sizes.ndims == [0, 2, 2, 2, 0, 0, 2, 0, 0, 2, 2, 0, 0, 2, 0, 0]
+    @test sizes.size_offset ==
+          [0, 12, 10, 8, 0, 0, 6, 0, 0, 4, 2, 0, 0, 0, 0, 0]
+    @test sizes.size == [1, 2, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2, 2, 2]
+    @test sizes.storage_offset ==
+          [0, 1, 5, 9, 11, 12, 13, 15, 16, 17, 21, 23, 24, 25, 27, 28, 29]
+    x1 = 1.0
+    x2 = 2.0
+    x3 = 3.0
+    x4 = 4.0
+    @test MOI.eval_objective(evaluator, [x1, x2, x3, x4]) == sqrt(14.0)
+    g = ones(4)
+    MOI.eval_objective_gradient(evaluator, g, [x1, x2, x3, x4])
+    @test g == [0.0, 1.0 / sqrt(14.0), 2.0 / sqrt(14.0), 3.0 / sqrt(14.0)]
+    return
+end
 
 end  # module