`BoundsError` when using `expand_derivatives` on a particular expression in `y`, but not the same expression in `x`

[krsteffen]

MWE (Julia Version 1.10.2, Symbolics v5.28.0, SymbolicUtils v1.5.1):

using Symbolics

expr_gen = (var, fun) -> Differential(var)(Differential(var)(((-Differential(var)(Differential(var)(fun))) / g(var))))

@syms y f(y) g(y) h(y)

expr = expr_gen(y, f(y))
expand_derivatives(expr)
expr = expr_gen(y, g(y))
expand_derivatives(expr) # Error
expr = expr_gen(y, h(y))
expand_derivatives(expr) # Error

The error message:

> show(err)
1-element ExceptionStack:
BoundsError: attempt to access 1-element Vector{Any} at index [2]
Stacktrace:
 [1] getindex(A::Vector{Any}, i1::Int64)
   @ Base ./essentials.jl:13
 [2] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Number}, simplify::Bool; occurrences::SymbolicUtils.BasicSymbolic{Real}) (repeats 2 times)
   @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/diff.jl:245
 [3] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Number}, simplify::Bool; occurrences::Nothing)
   @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/diff.jl:245
 [4] expand_derivatives
   @ ~/.julia/packages/Symbolics/Eas9m/src/diff.jl:171 [inlined]
 [5] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Number})
   @ Symbolics ~/.julia/packages/Symbolics/Eas9m/src/diff.jl:171
 [6] top-level scope
   @ REPL[9]:1

For some reason all of the following run successfully. The only difference is replacing y with x.

using Symbolics

expr_gen = (var, fun) -> Differential(var)(Differential(var)(((-Differential(var)(Differential(var)(fun))) / g(var))))

@syms x f(x) g(x) h(x)

expr = expr_gen(x, f(x))
expand_derivatives(expr)
expr = expr_gen(x, g(x))
expand_derivatives(expr)
expr = expr_gen(x, h(x))
expand_derivatives(expr)

[bowenszhu]

This is still a problem of ordering.

• <ₑ not fully reliable SymbolicUtils.jl#483
• better ordering of terms SymbolicUtils.jl#549

As a result, inner_args[i] and arguments(occurrences)[i] sometimes does not actually match inside expand_derivatives.

Symbolics.jl/src/diff.jl

Lines 239 to 245 in 45de9b8

	inner_args = arguments(arg)
	l = length(inner_args)
	exprs = []
	c = 0
	for i in 1:l
	t2 = expand_derivatives(D(inner_args[i]),false, occurrences=arguments(occurrences)[i])

[orebas]

I am getting an error, and I am not sure if it is the same error as above. Is there any workaround? Of note, I only seem to get an error when there is quotient somewhere in the expression.

Here is the stacktrace:

ERROR: LoadError: BoundsError: attempt to access 3-element Vector{Any} at index [4]
Stacktrace:
  [1] throw_boundserror(A::Vector{Any}, I::Tuple{Int64})
    @ Base ./essentials.jl:14
  [2] getindex(A::Vector{Any}, i::Int64)
    @ Base ./essentials.jl:904
  [3] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Real}, simplify::Bool; occurrences::SymbolicUtils.BasicSymbolic{Real}) (repeats 4 times)
    @ Symbolics ~/.julia/packages/Symbolics/2UpZj/src/diff.jl:257
  [4] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Real}, simplify::Bool; occurrences::Nothing)
    @ Symbolics ~/.julia/packages/Symbolics/2UpZj/src/diff.jl:257
  [5] expand_derivatives
    @ ~/.julia/packages/Symbolics/2UpZj/src/diff.jl:183 [inlined]
  [6] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Real})
    @ Symbolics ~/.julia/packages/Symbolics/2UpZj/src/diff.jl:183

unfortunately, I am having some trouble making an MWE. the input to expand_derivatives is this expression, which when I paste it into an individual file, (and declare all the variables and make some replacements), does not seem to error

Differential(t)((-b*Differential(t)(Differential(t)(In(t)))*S(t)*Differential(t)(N(t)) - b*Differential(t)(Differential(t)(Differential(t)(N(t))))*S(t)*In(t) - 2b*Differential(t)(S(t))*Differential(t)(In(t))*Differential(t)(N(t)) - 2b*Differential(t)(S(t))*Differential(t)(Differential(t)(N(t)))*In(t) - 2b*Differential(t)(In(t))*Differential(t)(Differential(t)(N(t)))*S(t) - b*Differential(t)(Differential(t)(S(t)))*Differential(t)(N(t))*In(t)) / (N(t)^2) + (-b*d*Differential(t)(Differential(t)(Differential(t)(N(t))))*S(t)*Tr(t) - 2b*d*Differential(t)(S(t))*Differential(t)(Differential(t)(N(t)))*Tr(t) - 2b*d*Differential(t)(S(t))*Differential(t)(N(t))*Differential(t)(Tr(t)) - b*d*Differential(t)(Differential(t)(S(t)))*Tr(t)*Differential(t)(N(t)) - 2b*d*Differential(t)(Differential(t)(N(t)))*S(t)*Differential(t)(Tr(t)) - b*d*S(t)*Differential(t)(N(t))*Differential(t)(Differential(t)(Tr(t)))) / (N(t)^2) + (2b*d*Differential(t)(S(t))*Tr(t)*(Differential(t)(N(t))^2) + 4b*d*Differential(t)(Differential(t)(N(t)))*S(t)*Tr(t)*Differential(t)(N(t)) + 2b*d*S(t)*(Differential(t)(N(t))^2)*Differential(t)(Tr(t))) / (N(t)^3) + ((-b*Differential(t)(Differential(t)(In(t)))*S(t) - 2b*Differential(t)(S(t))*Differential(t)(In(t)) - b*Differential(t)(Differential(t)(S(t)))*In(t))*Differential(t)(N(t)) + (-b*Differential(t)(S(t))*In(t) - b*Differential(t)(In(t))*S(t))*Differential(t)(Differential(t)(N(t)))) / (N(t)^2) + (3b*Differential(t)(Differential(t)(In(t)))*Differential(t)(S(t)) + 3b*Differential(t)(In(t))*Differential(t)(Differential(t)(S(t))) + b*S(t)*Differential(t)(Differential(t)(Differential(t)(In(t)))) + b*Differential(t)(Differential(t)(Differential(t)(S(t))))*In(t)) / N(t) + (2b*Differential(t)(S(t))*(Differential(t)(N(t))^2)*In(t) + 2b*Differential(t)(In(t))*S(t)*(Differential(t)(N(t))^2) + 4b*Differential(t)(Differential(t)(N(t)))*S(t)*Differential(t)(N(t))*In(t)) / (N(t)^3) + ((-b*d*Differential(t)(S(t))*Tr(t) - b*d*S(t)*Differential(t)(Tr(t)))*Differential(t)(Differential(t)(N(t))) + (-2b*d*Differential(t)(S(t))*Differential(t)(Tr(t)) - b*d*Differential(t)(Differential(t)(S(t)))*Tr(t) - b*d*S(t)*Differential(t)(Differential(t)(Tr(t))))*Differential(t)(N(t))) / (N(t)^2) + (3b*d*Differential(t)(S(t))*Differential(t)(Differential(t)(Tr(t))) + 3b*d*Differential(t)(Differential(t)(S(t)))*Differential(t)(Tr(t)) + b*d*S(t)*Differential(t)(Differential(t)(Differential(t)(Tr(t)))) + b*d*Differential(t)(Differential(t)(Differential(t)(S(t))))*Tr(t)) / N(t) + (-a - g)*Differential(t)(Differential(t)(Differential(t)(In(t)))) - ((2b*d*Differential(t)(S(t))*Differential(t)(Tr(t)) + b*d*Differential(t)(Differential(t)(S(t)))*Tr(t) + b*d*S(t)*Differential(t)(Differential(t)(Tr(t)))) / (N(t)^2))*Differential(t)(N(t)) - Differential(t)(N(t))*((b*Differential(t)(Differential(t)(In(t)))*S(t) + 2b*Differential(t)(S(t))*Differential(t)(In(t)) + b*Differential(t)(Differential(t)(S(t)))*In(t)) / (N(t)^2)) - 2N(t)*((-b*d*Differential(t)(S(t))*Tr(t)*Differential(t)(N(t)) - b*d*Differential(t)(Differential(t)(N(t)))*S(t)*Tr(t) - b*d*S(t)*Differential(t)(N(t))*Differential(t)(Tr(t))) / (N(t)^4))*Differential(t)(N(t)) - 2N(t)*(((-b*d*Differential(t)(S(t))*Tr(t) - b*d*S(t)*Differential(t)(Tr(t)))*Differential(t)(N(t))) / (N(t)^4))*Differential(t)(N(t)) - 2N(t)*(((-b*Differential(t)(S(t))*In(t) - b*Differential(t)(In(t))*S(t))*Differential(t)(N(t))) / (N(t)^4))*Differential(t)(N(t)) - 2N(t)*((-b*Differential(t)(S(t))*Differential(t)(N(t))*In(t) - b*Differential(t)(In(t))*S(t)*Differential(t)(N(t)) - b*Differential(t)(Differential(t)(N(t)))*S(t)*In(t)) / (N(t)^4))*Differential(t)(N(t)) - 3((2b*d*S(t)*Tr(t)*(Differential(t)(N(t))^2)) / (N(t)^6))*(N(t)^2)*Differential(t)(N(t)) - 3(N(t)^2)*((2b*S(t)*(Differential(t)(N(t))^2)*In(t)) / (N(t)^6))*Differential(t)(N(t)))

[orebas]

Here's a shorter MWE for this. Weirdly, renaming x to t seems not to fail.
Related discourse question: https://discourse.julialang.org/t/julia-symbolics-error-when-taking-derivative-of-a-quotient/118007

using ModelingToolkit

function main()

	@variables t a(t) b(t)
	D = Differential(t)
	q = a / b

	for i in 1:5
        q2 = deepcopy(q)
        q3 = D(q2)
        q4 = expand_derivatives(q3)
        q = deepcopy(q4)
        display(q)
	end
end
main()

[orebas]

Here is a more explicit and readable MWE:

using Symbolics

function main()
	@variables t a(t) b(t)
	D = Differential(t)
	q = D(
		(-3D(a) * D(D(D(b))) - D(D(D(D(b)))) * a - D(b) * D(D(D(a))) - 3D(D(b)) * D(D(a))) / (b^2) + (-6D(a) * (D(b)^3) - 18(D(b)^2) * a * D(D(b))) / (b^4) + (-D(a) * D(D(D(b))) - D(b) * D(D(D(a))) - 2D(D(b)) * D(D(a))) / (b^2) + D(D(D(D(a)))) / b +
		(-2(-D(a) * D(b) - a * D(D(b))) * D(D(b)) - 2(-2D(a) * D(D(b)) - D(b) * D(D(a)) - a * D(D(D(b)))) * D(b)) / (b^3) + (-D(b) * D(D(D(a))) - D(D(b)) * D(D(a))) / (b^2) +
		(8D(a) * D(b) * D(D(b)) + 2(D(b)^2) * D(D(a)) + 4D(b) * a * D(D(D(b))) + 4a * (D(D(b))^2)) / (b^3) + (4D(a) * D(b) * D(D(b)) + 2(D(b)^2) * D(D(a))) / (b^3) - D(b) * (D(D(D(a))) / (b^2)) -
		2((-2D(a) * D(D(b)) - D(b) * D(D(a)) - a * D(D(D(b)))) / (b^4)) * D(b) * b - 2D(b) * ((-D(a) * D(D(b)) - D(b) * D(D(a))) / (b^4)) * b - 2D(b) * b * ((-D(b) * D(D(a))) / (b^4)) - 3((2D(a) * (D(b)^2)) / (b^6)) * D(b) * (b^2) -
		3D(b) * ((-2(-D(a) * D(b) - a * D(D(b))) * D(b)) / (b^6)) * (b^2) - 3D(b) * ((2D(a) * (D(b)^2) + 4D(b) * a * D(D(b))) / (b^6)) * (b^2) - 4D(b) * ((-6(D(b)^3) * a) / (b^8)) * (b^3),
	)
    display(q)
    expand_derivatives(q)
end
main()

with this error:

ERROR: LoadError: BoundsError: attempt to access 1-element Vector{Any} at index [2]
Stacktrace:
  [1] throw_boundserror(A::Vector{Any}, I::Tuple{Int64})
    @ Base ./essentials.jl:14
  [2] getindex(A::Vector{Any}, i::Int64)
    @ Base ./essentials.jl:915
  [3] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Real}, simplify::Bool; occurrences::SymbolicUtils.BasicSymbolic{Real}) (repeats 4 times)
    @ Symbolics ~/.julia/packages/Symbolics/zuLwn/src/diff.jl:257
  [4] expand_derivatives(O::SymbolicUtils.BasicSymbolic{Real}, simplify::Bool; occurrences::Nothing)
    @ Symbolics ~/.julia/packages/Symbolics/zuLwn/src/diff.jl:257
  [5] expand_derivatives(n::Num, simplify::Bool; occurrences::Nothing)
    @ Symbolics ~/.julia/packages/Symbolics/zuLwn/src/diff.jl:299
  [6] expand_derivatives
    @ ~/.julia/packages/Symbolics/zuLwn/src/diff.jl:298 [inlined]
  [7] expand_derivatives(n::Num)
    @ Symbolics ~/.julia/packages/Symbolics/zuLwn/src/diff.jl:298
  [8] main()
    @ Main ~/MWE-test/v5.jl:14
  [9] top-level scope
    @ ~/MWE-test/v5.jl:16
 [10] include(fname::String)
    @ Main ./sysimg.jl:38
 [11] top-level scope
    @ REPL[14]:1
in expression starting at /home/orebas/MWE-test/v5.jl:16

[karlwessel]

This bug keeps hitting me, especially for long but simple expressions, where CAS would actually start to be useful.

Is there some quick fix or workaround where we avoid using occurrences in expand_derivatives, since it seems to be used for optimization only?

I would prefer a slower but working method.

[karlwessel]

This is a Heisenbug, displaying Debug information changes the state of the system. The occurences tree/cache is built using the unsorted arguments list:

Symbolics.jl/src/diff.jl

Line 102 in c3f54a6

args = map(a->occursin_info(x, a, operation(expr) !== getindex), arguments(expr))

)

But as soon as the parent term of the arguments is displayed the arguments get sorted permanently:
https://github.com/JuliaSymbolics/SymbolicUtils.jl/blob/f9b0ade4a407b6c090ea8ce39fc9ca4ef238826e/src/types.jl#L139

This changes their order, while occurences still uses the old ordering. This definitely messes things up during debugging.

[karlwessel]

Is there some quick fix or workaround where we avoid using occurrences in expand_derivatives, since it seems to be used for optimization only?

I implemented a workaround in PR #1353.

[krsteffen]

The code runs successfully with Symbolics v6.27.0. I haven't tested earlier versions. Closing