# Looking for ways to debug mathematical algorithms

My research involves developing of thermodynamical solvers. They need equation solvers, minimization, automatic differentiation, etc. When I develop a new solver for a physical problem or face a bug in an existing one, I want to debug.

For example, before solving an equation `f(x) = 0`, I want to ensure that `f` has specific (for a problem or a solver) properties. In such situations, the ideal debug for me is to pause a program, plot `f` and inspect it.

I’m looking for handy ways for such types of debug.

I’ve tried

• logging of function’s arguments and values. The main problem here is that the argument’s precision is unknown a priori.
• saving function (closure, usually) as object and running a program in interactive mode (`-i` flag or run from REPL). When the program is finished, I used a plotting library. This approach is more flexible, but, in my implementation, solver must return debug-objects.

My (quite vague) piece of advice here would be: Encode as much program logic as possible in terms of data rather than procedure.

Data can be easily inspected and tested for its properties while procedures are much harder to analyse.

For example, turning a for loop with a nested if statement that reduces over some vector into a (sparse) matrix multiplication could make that algorithm much more testable.

I find `Infiltrator` very useful when I want to inspect the intermediate objects created by my code:

It is a very versatile and flexible tool, but one way of using it could be to exfiltrate `f` to the “safehouse” when your code is running, and inspect it afterwards.

2 Likes

Seems promising for me, thanks!

Thanks. Could you clarify the advice on this example?

``````# E.g. scalar function representing a physical model, say, an equation of state.
foo(x, y) = ...

# E.g. a solver of a physical problem.
function bar(xrange)
# ...
for x in xrange
target(y) = foo(x, y) - foo(2 * x, y)
root = find_root(solver, target, ymin, ymax)
# Do something with root...
end
# ...
end
``````

For me, the main point here is creating a throughaway closure `target(y)` which I want to debug (look at its plot in a “bad” case for `solver`). If I would apply your advice, I would

1. define in outter scope `target_closure_creator(x) -> target(y)`,
2. log necessary input,
3. go to REPL, recreate `target(y)` with `target_closure_creator` and log,
4. and, finally, inspect `target(y)`?

Let me illustrate my point with a simple example. Let’s say I want an algorithm to swap the first two elements of a vector. That operation should conserve the norm of said vector.

A procedural approach is

``````julia> f(a) = [a, a]
f (generic function with 1 method)

julia> f([1, 2])
2-element Vector{Int64}:
2
1

julia> using LinearAlgebra, Test

julia> a = rand(2)
2-element Vector{Float64}:
0.23595638406718
0.26032672255472455

julia> @test norm(a) == norm(f(a))
Test Passed
Expression: norm(a) == norm(f(a))
Evaluated: 0.35134800078859574 == 0.35134800078859574

``````

Instead, in the data-oriented approach I can test if the data, in this case a matrix, is unitary:

``````julia> using LinearAlgebra

julia> const swap_mat = [0 1; 1 0]
2×2 Matrix{Int64}:
0  1
1  0

julia> g(a) = swap_mat * a
g (generic function with 1 method)

julia> g([3, 4])
2-element Vector{Int64}:
4
3

julia> using Test

julia> isunitary(m) = iszero(norm(m' * m - I))
isunitary (generic function with 1 method)

julia> @test isunitary(swap_mat)
Test Passed
Expression: isunitary(swap_mat)
``````

As we scale up complexity, the data-oriented approach has a lot of benefits. Of course, it is absolutely not trivial to convert from procedure to data. But especially in scientific computing it can be quite helpful.

1 Like