# Comparing equivalency of ForwardDiff.Jacobian

**URL:** <https://discourse.julialang.org/t/comparing-equivalency-of-forwarddiff-jacobian/132194>\
**Category:** General Usage\
**Created:** [September 8, 2025, 3:07pm UTC](https://discourse.julialang.org/t/comparing-equivalency-of-forwarddiff-jacobian/132194 "2025-09-08T15:07:36Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![JAY\_SHUKLA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jay_shukla/32/218577_2.png) [@JAY\_SHUKLA](https://discourse.julialang.org/u/JAY_SHUKLA)\
**Post date:** [September 8, 2025, 3:07pm UTC](https://discourse.julialang.org/t/comparing-equivalency-of-forwarddiff-jacobian/132194/1 "2025-09-08T15:07:36Z")

</div>

Hello,

I am trying to work on a neural network framework that identifies reaction pathways autonomously and currently trying to add an analysis tool called degree of rate control. I am comparing two approaches where I am defining the kinetic ODEs of the reactions explicitly (ODE.jl) and then using neural network definition via matrix multiplication (CRNN.jl). Ideally both the scripts should produce the same results after calculating the Jacobian matrix but for some reason CRNN.jl does not produce the same result as ODE.jl. Can anyone help me diagnose why?

Since I am a new user and can’t upload attachments here, please find the scripts here:

ODE.jl

using DifferentialEquations  
using Plots  
using DiffEqSensitivity, ForwardDiff, DelimitedFiles

ns = 2  
nr = 2  
k = [1e-5, 1]  
alg = RK4()  
b0 = 0  
lb = 1e-5  
ub = 1e1  
nt = 15

function trueODEfunc(dydt, y, k, t)  
aA = 1  
aB = 1  
kr = 0  
dydt[1] = -k[1] \* aA \* y[1] + kr \* y[2] + k[2] \* aB \* y[2];  
dydt[2] = -dydt[1]  
#dydt[2] = k[1] \* aA \* y[1] - k[2] \* y[2] - k3 \* aB \* y[2];  
end

#u0 = zeros(Float64, ns);  
free\_ic = 1 - 1 / (1 + 1.0e5);  
adsorbed\_ic = (1 / (1 + 1.0e5));  
u0 = [free\_ic, adsorbed\_ic]  
tspan = (0., 4.)  
tsteps = LinRange(0, 4, nt)

prob = ODEProblem(trueODEfunc, u0, tspan, k);  
ode\_sol = Array(solve(prob, alg, saveat=tsteps))  
data\_matrix = hcat(tsteps, ode\_sol[1, :], ode\_sol[2, :])

# Add column headers

headers = [“Time\_Steps” “Free\_Sites” “Adsorbed\_Sites”]

# Combine headers and data

output\_data = vcat(headers, data\_matrix)

# Save to CSV file

writedlm(“ODE\_ipynb.csv”, output\_data, ‘,’)

function target\_rate(y, k)  
aB = 1  
rate = y[2] \* k[2] \* aB  
return rate  
end

function rate\_wrapper(lnk)  
k = exp.(lnk)  
\_prob = remake(prob, p=k)  
sol = Array(solve(\_prob, alg, saveat=tsteps, sensealg=ForwardDiffSensitivity()))  
println(size(sol))  
k\_matrix = reshape(k, 1, size(k, 1))  
k\_repeat = repeat(k\_matrix, nt, 1)  
rate = Array{Real, 2}(undef, nt, 1)  
for i in 1:nt  
rate[i, 1] = target\_rate(sol[:, i], k\_repeat[i, :])  
end  
println(“Rate”)  
println(rate)  
return log.(rate)  
end

drc = ForwardDiff.jacobian(rate\_wrapper, log.(k))

plt = plot()  
plot!(plt, tsteps, drc[:, 1],  
linewidth=3, xlabel=“Time (s)”, ylabel=“Degree of Rate Control”,  
label=“DRC-1”)  
plot!(plt, tsteps, drc[:, 2], linewidth=3, label=“DRC-2”)  
png(plt, string(“DRC\_ODE”))

CRNN.jl

using DifferentialEquations  
using Plots  
using DiffEqSensitivity, ForwardDiff, DelimitedFiles

ns = 2  
nr = 2  
k = [1e-5, 1]  
alg = RK4()  
b0 = 0  
lb = 1e-5  
ub = 1e1  
nt = 15

#u0 = zeros(Float64, ns);  
free\_ic = 1 - 1 / (1 + 1.0e5);  
adsorbed\_ic = (1 / (1 + 1.0e5));  
u0 = [free\_ic, adsorbed\_ic]  
tspan = (0., 4.)  
tsteps = LinRange(0, 4, nt)

function p2vec(p)  
w\_b = p[1:nr] .+ b0;  
# More robust reshaping that works with dual numbers  
remaining\_params = p[nr + 1:end]  
w\_out = reshape(remaining\_params, ns, nr);  
# w\_out = clamp.(w\_out, -2.5, 2.5);  
w\_in = clamp.(-w\_out, 0, 2.5);  
return w\_in, w\_b, w\_out  
end

function display\_p(p)  
w\_in, w\_b, w\_out = p2vec(p);  
println(“species (column) reaction (row)”)  
println(“w\_in”)  
show(stdout, “text/plain”, round.(w\_in’, digits=3))

```
println("\nw_b")
show(stdout, "text/plain", round.(exp.(w_b'), digits=6))

println("\nw_out")
show(stdout, "text/plain", round.(w_out', digits=3))
println("\n\n")

```

end

function crnn(du, u, p, t)  
w\_in, w\_b, w\_out = p2vec(p);  
w\_in\_x = w\_in’ \* @. log(clamp(u, lb, ub));  
du .= w\_out \* @. exp(w\_in\_x + w\_b);  
end

p = [log(1e-5_1), log(1_1), -1, 1, 1, -1]  
display\_p(p)

prob = ODEProblem(crnn, u0, tspan, p)

function predict\_neuralode(prob, u0)  
sol = Array(solve(prob, alg, u0=u0, saveat=tsteps))  
return sol  
end

sol = predict\_neuralode(prob, u0)  
data\_matrix = hcat(tsteps, sol[1, :], sol[2, :])

# Add column headers

headers = [“Time\_Steps” “Free\_Sites” “Adsorbed\_Sites”]

# Combine headers and data

output\_data = vcat(headers, data\_matrix)

# Save to CSV file

writedlm(“CRNN\_ipynb.csv”, output\_data, ‘,’)

function target\_rate(u, p)  
w\_in, w\_b, w\_out = p2vec(p);  
w\_in\_x = w\_in’ \* @. log(clamp(u, lb, ub));  
rate\_all\_reaction = @. exp(w\_in\_x + w\_b);  
println(size(rate\_all\_reaction))  
target\_rate = rate\_all\_reaction[2]  
return target\_rate  
end

function rate\_wrapper(p\_new)  
\_prob = remake(prob, p=p\_new)  
sol = predict\_neuralode(\_prob, u0)  
rate = Vector{eltype(p\_new)}(undef, nt) # Use eltype to handle dual numbers  
for i in 1:nt  
rate[i] = target\_rate(sol[:, i], p\_new)  
end  
println(“Rate”)  
println(rate)  
return log.(rate)  
end

# Compute Jacobian with respect to all parameters

drc = ForwardDiff.jacobian(rate\_wrapper, p)

# Plot only the first two columns (corresponding to the rate constants)

plt = plot()  
plot!(plt, tsteps, drc[:, 1],  
linewidth=3, xlabel=“Time (s)”, ylabel=“Degree of Rate Control”,  
label=“DRC-1 (rate constant 1)”)  
plot!(plt, tsteps, drc[:, 2], linewidth=3,  
label=“DRC-2 (rate constant 2)”)  
png(plt, string(“DRC\_CRNN”))

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 8, 2025, 4:38pm UTC](https://discourse.julialang.org/t/comparing-equivalency-of-forwarddiff-jacobian/132194/2 "2025-09-08T16:38:01Z")

</div>

From what you showed… you just built the neural network but never trained the coefficients?

---

<div class="post-metadata">

**Author:** ![JAY\_SHUKLA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jay_shukla/32/218577_2.png) [@JAY\_SHUKLA](https://discourse.julialang.org/u/JAY_SHUKLA)\
**Post date:** [September 8, 2025, 5:42pm UTC](https://discourse.julialang.org/t/comparing-equivalency-of-forwarddiff-jacobian/132194/3 "2025-09-08T17:42:59Z")

</div>

Correct, I have not trained parameters because I am using pre-trained parameters via the following line:

p = [log(1e-5_1), log(1_1), -1, 1, 1, -1]

Basically I am trying to calculate a property called “Degree of Rate Control” or DRC. DRC is partial derivative of log of rate of the reaction with respect to the log of rate constant. This is the paper I am referring to: Evaluation of the degree of rate control via automatic differentiation

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [September 9, 2025, 6:02am UTC](https://discourse.julialang.org/t/comparing-equivalency-of-forwarddiff-jacobian/132194/4 "2025-09-09T06:02:27Z")

</div>

How much is the difference?

---

<div class="post-metadata">

**Author:** ![JAY\_SHUKLA](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jay_shukla/32/218577_2.png) [@JAY\_SHUKLA](https://discourse.julialang.org/u/JAY_SHUKLA)\
**Post date:** [September 9, 2025, 1:08pm UTC](https://discourse.julialang.org/t/comparing-equivalency-of-forwarddiff-jacobian/132194/5 "2025-09-09T13:08:21Z")

</div>

The pre-defined parameters when put in the crnn function will replicate the ODEs defined in trueODEfunc in ODE.jl. So there is not any difference as such.

However, I found the fix. The issue was clamping of the solution of the ODEs in CRNN.jl and incorrect definition in the function rate\_wrapper.

Thanks for your help!!