# Looking for ways to debug mathematical algorithms

**URL:** <https://discourse.julialang.org/t/looking-for-ways-to-debug-mathematical-algorithms/95435>\
**Category:** General Usage\
**Tags:** debugging\
**Created:** [March 2, 2023, 11:29am UTC](https://discourse.julialang.org/t/looking-for-ways-to-debug-mathematical-algorithms/95435 "2023-03-02T11:29:43Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![stepanzh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanzh/32/34851_2.png) [@stepanzh](https://discourse.julialang.org/u/stepanzh)\
**Post date:** [March 2, 2023, 11:29am UTC](https://discourse.julialang.org/t/looking-for-ways-to-debug-mathematical-algorithms/95435/1 "2023-03-02T11:29:43Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [March 2, 2023, 12:35pm UTC](https://discourse.julialang.org/t/looking-for-ways-to-debug-mathematical-algorithms/95435/2 "2023-03-02T12:35:44Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![ffevotte](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ffevotte/32/6587_2.png) [@ffevotte](https://discourse.julialang.org/u/ffevotte)\
**Post date:** [March 2, 2023, 1:07pm UTC](https://discourse.julialang.org/t/looking-for-ways-to-debug-mathematical-algorithms/95435/3 "2023-03-02T13:07:12Z")

</div>

> [@stepanzh](#):
>
> 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 find [`Infiltrator`](https://juliadebug.github.io/Infiltrator.jl/stable/) very useful when I want to inspect the intermediate objects created by my code:

> **[GitHub - JuliaDebug/Infiltrator.jl: No-overhead breakpoints in Julia](https://github.com/JuliaDebug/Infiltrator.jl)**
>
> No-overhead breakpoints in Julia. Contribute to JuliaDebug/Infiltrator.jl development by creating an account on GitHub.

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.

---

<div class="post-metadata">

**Author:** ![stepanzh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanzh/32/34851_2.png) [@stepanzh](https://discourse.julialang.org/u/stepanzh)\
**Post date:** [March 2, 2023, 1:58pm UTC](https://discourse.julialang.org/t/looking-for-ways-to-debug-mathematical-algorithms/95435/4 "2023-03-02T13:58:25Z")

</div>

Seems promising for me, thanks!

---

<div class="post-metadata">

**Author:** ![stepanzh](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stepanzh/32/34851_2.png) [@stepanzh](https://discourse.julialang.org/u/stepanzh)\
**Post date:** [March 2, 2023, 2:32pm UTC](https://discourse.julialang.org/t/looking-for-ways-to-debug-mathematical-algorithms/95435/5 "2023-03-02T14:32:56Z")

</div>

Thanks. Could you clarify the advice on this example?

```julia
# 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)`?

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [March 2, 2023, 2:55pm UTC](https://discourse.julialang.org/t/looking-for-ways-to-debug-mathematical-algorithms/95435/6 "2023-03-02T14:55:28Z")

</div>

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