Feature request: order of execution

[martijnende]

Section 4 of the paper of Treister & Haber details an efficient approach to computing the sensitivities of the Eikonal solution (corresponding Julia code: here). This approach requires the order in which the fast-marching algorithm was executed and the order of the finite difference operator (1st or 2nd order). Would it be possible to store these two quantities in eikonalfm so that they can be accessed for the sensitivity calculations?

[kevinganster]

When i wrote this package, i didn't look into the sensitivities part of the paper. What exactly would you need? The second array should contain the finite difference order (1 or 2) used, so for each grid-point we have a 1 if first order was used there and a 2 for second order. Correct? I don't really know what you mean with "order in which the fast-marching algorithm was executed". Is this the sequence which grid-point was evaluated when? So '0' for the source, '1' for a neighbour of the source and so on?

In general this would all be possible by simply creating two new arrays in the C-program and returning them at the end. But this also means that even if i don't want those two arrays, they would still need memory. I don't know how to avoid overhead when i don't want these arrays (i.e. in a "normal" run of the program). I would need to look into it, how this would be done as efficiently as possible in C++.

[martijnende]

Ideally, I'd need a function I can call to get a Jacobian-vector product for the "current" (= last run) forward simulation. For this one needs to construct a finite-difference matrix with the appropriate coefficients (forward 1st/2nd order or backward 1st/2nd order), and solve a system of equations in the order of the fast-marching operation (as you said: 0 at the source, 1 at the nearest neighbour, 2 the next nearest neighbour, etc.). Just having access to the order of execution and the appropriate derivative coefficients might be enough for me to write a Jacobian-vector routine in Python (this code doesn't look too challenging)

In the Julia code of Treister & Haber, they indeed allocate and fill two vectors during the main solver loop. In your case, you could consider using a flag. Since this flag is constant, evaluation of the conditionals should not incur much overhead if this output is not required (standard run).

It has been several years since I last touched a piece of C++ code, but if you're up for it we could look at this together?

[kevinganster]

I've succeeded so far to add an optional parameter to the C++-Interface called output_sensitivities, which when set to True, will store the execution sequence and orders (and return them as additional parameters). I'm still not sure if the "orders"-array i produce is the one you need. Since the finite differences are used in every dimension, i thought i should store the order of the FD for each dimension. The program

import numpy as np
import eikonalfm

c = np.ones((11, 11))
x_s = (5, 5)
dx = (0.1, 0.1)
order = 2

tau_fm, sequence, orders = eikonalfm.factored_fast_marching(c, x_s, dx, order, output_sensitivities=True)

print(sequence)
print(orders)

now produces the output

[[119 111 105  89  81  77  82  94 101 109 117]
 [115  97  71  66  49  45  50  64  72  98 112]
 [103  69  59  37  29  25  34  39  57  75 106]
 [ 91  61  41  24  17  11  19  21  44  65  93]
 [ 86  52  35  20   5   1   7  13  32  55  87]
 [ 80  48  28  12   3   0   2   9  27  47  78]
 [ 85  53  33  16   8   4   6  14  31  56  83]
 [ 95  62  38  22  15  10  18  23  43  63  96]
 [104  70  60  42  30  26  36  40  58  74 102]
 [114 100  73  67  54  46  51  68  76  99 116]
 [120 110 108  92  88  79  84  90 107 113 118]]
[[[ 2  2  2  2  2  2  2  2  2  2  2]
  [ 2  2  2  2  2  2  2  2  2  2  2]
  [ 2  2  2  2  2  2  2  2  2  2  2]
  [ 2  2  2  2  2  2  2  2  2  2  2]
  [ 1  1  1  1  1  1  1  1  1  1  1]
  [ 0  0  0  0  0  0  0  0  0  0  0]
  [-1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1]
  [-2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2]
  [-2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2]
  [-2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2]
  [-2 -2 -2 -2 -2 -2 -2 -2 -2 -2 -2]]

 [[ 2  2  2  2  1  0 -1 -2 -2 -2 -2]
  [ 2  2  2  2  1  0 -1 -2 -2 -2 -2]
  [ 2  2  2  2  1  0 -1 -2 -2 -2 -2]
  [ 2  2  2  2  1  0 -1 -2 -2 -2 -2]
  [ 2  2  2  2  1  0 -1 -2 -2 -2 -2]
  [ 2  2  2  2  1  0 -1 -2 -2 -2 -2]
  [ 2  2  2  2  1  0 -1 -2 -2 -2 -2]
  [ 2  2  2  2  1  0 -1 -2 -2 -2 -2]
  [ 2  2  2  2  1  0 -1 -2 -2 -2 -2]
  [ 2  2  2  2  1  0 -1 -2 -2 -2 -2]
  [ 2  2  2  2  1  0 -1 -2 -2 -2 -2]]]

Is this what you need?

[martijnende]

That's a quick implementation! I think this is indeed what is needed to compute the sensitivities w.r.t. to the velocity model. If you could push the code to a development branch, I could try it out and implement something in Python to see if it works.

[kevinganster]

Done. See the "dev" branch.

[martijnende]

Perfect, I'll start working on the Python implementation.