Looking for ways to debug mathematical algorithms

General Usage
1

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.
2

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.

3

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
4

Seems promising for me, thanks!

5

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)?
6

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[2], a[1]]
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