Is there any quicker way to get shortest path

[Peter Pan]

I used the code optimal_path_2d in wiki page, but it run very slow, the part of runge-kutta costs much, is there some other ways to finish it?

[Jason Furtney]

Dear Peter Pan,

Yes, I would imagine there would be a faster way to do this. Could you
send me the example that is running slowly please?

Thanks,
Jason

On Mon, Apr 4, 2016 at 10:06 AM, Peter Pan <panle...@gmail.com> wrote:
> I used the code optimal_path_2d in wiki page, but it run very slow, the part
> of runge-kutta costs much, is there some other ways to finish it?
>

> --
> You received this message because you are subscribed to the Google Groups
> "scikit-fmm" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to scikit-fmm+...@googlegroups.com.
> To post to this group, send email to sciki...@googlegroups.com.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/scikit-fmm/5d13f3fa-84a3-409d-bc61-87581fe5d098%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Lei Pan]

Dear Jason,

The attachment is my code, I use [`line_profiler`](https://github.com/rkern/line_profiler) to record time consuming. I added a break into the original optimal_path code, because it doesn't need too many iterations in Rung-Kutta, and now it needs 40-70 iterations for this example.

if add @profile to function get_kernel, the output:

```
6 @profile
7 def get_kernel(data, nx, ny, loc_source, loc_receiver, velocity, num_interval):
8 """Compute kernel function for ray theory
9 get kernel when given location of source and receiver based
10 on some velocity model
11
12
13 Parameter
14 =========
15
16 data: tuple
17 travel time data (num_of_receivers)
18 nx: int
19 number of grids along x axis
20 ny: int
21 number of grids along y axis
22 loc_source: tuple
23 (int, int), e.g. (2, 3)
24 loc_receiver: tuple
25 ((int, int), (int, int), ...), e.g. ((1, 2), (3, 4))
26 velocity: array
27 shape (ny, nx)
28 num_interval: int
29 number of sampling to get grids which shortest path goes across
30
31 Return
32 ======
33
34 kernel: array
35 shape (ny, nx)"""
36
37 # default set dx = 1
38 1 18 18.0 0.0 dx = 1
39
40 1 52 52.0 0.0 coordx = np.arange(nx)
41 1 6 6.0 0.0 coordy = np.arange(ny)
42 1 494 494.0 0.0 X, Y = np.meshgrid(coordx, coordy)
43
44 1 3 3.0 0.0 xs, ys = loc_source
45 1 148 148.0 0.0 phi = -1 * np.ones_like(X)
46 1 1168 1168.0 0.0 phi[np.logical_and(np.abs(X - xs) <= dx, np.abs(Y - ys) <= dx)] = 1
47 1 31392 31392.0 0.3 travel_time = skfmm.travel_time(phi, velocity, dx)
48
49 1 61 61.0 0.0 kernel = np.zeros_like(velocity)
50 1 4183 4183.0 0.0 t_rec = get_time(travel_time, coordx, coordy, loc_receiver)
51 1 1 1.0 0.0 i = 0
52 1 39 39.0 0.0 v2 = velocity**2
53 400 655 1.6 0.0 for rec in loc_receiver:
54 399 9107368 22825.5 98.0 yl, xl = optimal_path(travel_time, rec, loc_source, dx, coordx, coordy)
55 399 117388 294.2 1.3 grid_indx = get_curve(xl, yl, num_interval)
56 399 1170 2.9 0.0 dt = t_rec[i] - data[i]
57 399 4020 10.1 0.0 ix, iy = zip(*grid_indx)
58 399 27472 68.9 0.3 kernel[ix, iy] += -dt / v2[ix, iy]
59 399 727 1.8 0.0 i += 1
60 1 1 1.0 0.0 return kernel
```
The time that optimal_path costs is 98%, which is the main cost for the entire program.


The cost of time in optimal_path is below:
```
62 @profile
63 def optimal_path(travel_time,
64 starting_point,
65 ending_point,
66 dx,
67 coordx,
68 coordy,
69 N=10000):
70 """
71 Find the optimal path from starting_point to the zero contour
72 of travel_time. dx is the grid spacing
73 Solve the equation x_t = - grad t / | grad t |
74 """
75 399 278446 697.9 2.9 grad_t_y, grad_t_x = np.gradient(travel_time, dx)
76 399 2728 6.8 0.0 if isinstance(travel_time, np.ma.MaskedArray):
77 grad_t_y[grad_t_y.mask] = 0.0
78 grad_t_y = grad_t_y.data
79 grad_t_x[grad_t_x.mask] = 0.0
80 grad_t_x = grad_t_x.data
81
82 399 1562867 3917.0 16.3 gradx_interp = RectBivariateSpline(coordy, coordx, grad_t_x)
83 399 1511659 3788.6 15.8 grady_interp = RectBivariateSpline(coordy, coordx, grad_t_y)
84
85 399 1425 3.6 0.0 def get_velocity(position):
86 """ return normalized velocity at pos """
87 x, y = position
88 vel = np.array([gradx_interp(y, x)[0][0], grady_interp(y, x)[0][0]])
89 return vel / np.linalg.norm(vel)
90
91 399 885 2.2 0.0 def euler_point_update(pos, ds):
92 return pos - get_velocity(pos) * ds
93
94 399 814 2.0 0.0 def runge_kutta(pos, ds):
95 """ Fourth order Runge Kutta point update """
96 k1 = ds * get_velocity(pos)
97 k2 = ds * get_velocity(pos - k1 / 2.0)
98 k3 = ds * get_velocity(pos - k2 / 2.0)
99 k4 = ds * get_velocity(pos - k3)
100 return pos - (k1 + 2 * k2 + 2 * k3 + k4) / 6.0
101
102 399 211439 529.9 2.2 x = runge_kutta(starting_point, dx)
103 399 1115 2.8 0.0 xl, yl = [], []
104 399 870 2.2 0.0 xe, ye = ending_point
105 21687 46534 2.1 0.5 for i in range(N):
106 21687 59000 2.7 0.6 xl.append(x[1])
107 21687 51408 2.4 0.5 yl.append(x[0])
108 21687 5695582 262.6 59.4 x = runge_kutta(x, dx)
109 21687 112979 5.2 1.2 distance = ((x[1] - xl[-1])**2 + (x[0] - yl[-1])**2)**0.5
110 21687 55296 2.5 0.6 if distance < dx * 0.9999:
111 399 918 2.3 0.0 break
112 399 823 2.1 0.0 return yl, xl
```

Thanks

> You received this message because you are subscribed to a topic in the Google Groups "scikit-fmm" group.
> To unsubscribe from this topic, visit https://groups.google.com/d/topic/scikit-fmm/XvV_KhCUN5c/unsubscribe.
> To unsubscribe from this group and all its topics, send an email to scikit-fmm+...@googlegroups.com.

> To post to this group, send email to sciki...@googlegroups.com.

> To view this discussion on the web visit https://groups.google.com/d/msgid/scikit-fmm/CAOU2ZSW9Qx7oZhocJLqHHmFEbbGMC_k%3DpV1Ux%3DH4OfbPY1g7Nw%40mail.gmail.com.

[Lei Pan]

I'm sorry the last attachment missed a file. So I update a new one.

[Lei Pan]

Dear Jason,

I find the reason why it is slow, I call optimal_path for 100 times, and the record time is for the overall time. After applying only 1 time, it costs 60% in the entire program, and travel_time almost 30%.
I used FDTD to compute waveform and get arrival time, which is really slow for the purpose of testing inversion method. So I begin to try FMM, and I have thought it will be much faster. Yes, it is much faster to compute travel time, but it needs much more time if I need the shortest path, especially when the number of path is large.

Thanks,

> 在 2016年4月5日，下午2:19，Lei Pan <panle...@gmail.com> 写道：
>
> I'm sorry the last attachment missed a file. So I update a new one.<test_new.tar.gz>

>> 在 2016年4月4日，下午11:51，Jason Furtney <jkfu...@gmail.com> 写道：
>>

>> You received this message because you are subscribed to a topic in the Google Groups "scikit-fmm" group.
>> To unsubscribe from this topic, visit https://groups.google.com/d/topic/scikit-fmm/XvV_KhCUN5c/unsubscribe.
>> To unsubscribe from this group and all its topics, send an email to scikit-fmm+...@googlegroups.com.

>> To post to this group, send email to sciki...@googlegroups.com.

>> To view this discussion on the web visit https://groups.google.com/d/msgid/scikit-fmm/CAOU2ZSW9Qx7oZhocJLqHHmFEbbGMC_k%3DpV1Ux%3DH4OfbPY1g7Nw%40mail.gmail.com.

[Jason Furtney]

Dear Lei Pan,

Thank you for sending your example code. I ran this with the cProfile
module like this:

def gk():
return kernel.get_kernel(data, nx, ny, loc_src, loc_rec, speed0, 20)

import cProfile
cProfile.run('gk()', 'data')

import pstats
p=pstats.Stats('data')
p.strip_dirs().sort_stats('cumulative').print_stats(10)

and got:

ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 4.990 4.990 <string>:1(<module>)
1 0.000 0.000 4.990 4.990 test_kernel.py:41(gc)
1 0.023 0.023 4.990 4.990 kernel.py:6(get_kernel)
399 0.090 0.000 4.920 0.012 kernel.py:61(optimal_path)
799 2.306 0.003 2.343 0.003 fitpack2.py:1137(__init__)
22086 0.378 0.000 2.274 0.000 kernel.py:92(runge_kutta)
88344 0.515 0.000 1.896 0.000 kernel.py:83(get_velocity)
176689 0.465 0.000 0.847 0.000 fitpack2.py:774(__call__)
88344 0.229 0.000 0.433 0.000 linalg.py:1976(norm)
441722 0.163 0.000 0.424 0.000 numeric.py:406(asarray)

Just under half the time is spent in fitpack2. __init__.
It looks like you only need to calculate the interpolation weights
once for each travel time. I factored this out of the loop over the
rec points and got:

ncalls tottime percall cumtime percall filename:lineno(function)
1 0.000 0.000 2.487 2.487 <string>:1(<module>)
1 0.000 0.000 2.487 2.487 test_kernel.py:41(gc)
1 0.021 0.021 2.487 2.487 kernel.py:6(get_kernel)
399 0.086 0.000 2.414 0.006 kernel.py:76(optimal_path)
22086 0.377 0.000 2.301 0.000 kernel.py:98(runge_kutta)
88344 0.532 0.000 1.925 0.000 kernel.py:89(get_velocity)
176689 0.473 0.000 0.855 0.000 fitpack2.py:774(__call__)
88344 0.226 0.000 0.434 0.000 linalg.py:1976(norm)

My changes are attached, I did not test if the answer is the same
after this change.

These are the key lines:
22086 0.377 0.000 2.301 0.000 kernel.py:98(runge_kutta)
88344 0.532 0.000 1.925 0.000 kernel.py:89(get_velocity)
176689 0.473 0.000 0.855 0.000 fitpack2.py:774(__call__)

The call into fitpack is probably not going to be any faster. If you
want to go faster I would try to re-write get_velocity and runge_kutta
in Cython. NumPy arrays have (relatively) large overheads for the
length 2 arrays used in this example. In the Cython code try working
with the vector component doubles individually instead of using numpy.
Some other ideas to speed up the code could be to write your own
interpolation class which only calculates interpolation weights for
cells that are actually used. Another idea could be to use some kind
of numerical ODE solver to solve the shortest path.

Thanks,
Jason

> To view this discussion on the web visit https://groups.google.com/d/msgid/scikit-fmm/1A4B610F-EFB6-4863-B3B3-97230B28FF45%40gmail.com.

[Lei Pan]

Yes, you're right. The interpolation weights should be got out of the loop.  I run it again. The major part of consuming time is get_velocity now. 

Total time: 0.147592 s

File: kernel_new.py

Function: get_velocity at line 6

Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

     6                                               @profile

     7                                               def get_velocity(position):

     8                                                   """ return normalized velocity at pos """

     9      2000         5024      2.5      3.4          x, y = position

    10      2000        75148     37.6     50.9          vel = np.array([X_interp(y, x)[0][0], Y_interp(y, x)[0][0]])

    11      2000        67420     33.7     45.7          return vel / np.linalg.norm(vel)

I google np.linalg.norm and find a way to run it faster, the link is [](http://stackoverflow.com/questions/2850743/numpy-how-to-quickly-normalize-many-vectors), you can see the last one.

Total time: 0.081113 s

File: kernel_new.py

Function: get_velocity at line 7

Line #      Hits         Time  Per Hit   % Time  Line Contents

==============================================================

     7                                               @profile

     8                                               def get_velocity(position):

     9                                                   """ return normalized velocity at pos """

    10      2000         4221      2.1      5.2          x, y = position

    11      2000        64130     32.1     79.1          vel = np.array([X_interp(y, x)[0][0], Y_interp(y, x)[0][0]])

    12                                                   # return vel / np.linalg.norm(vel)

    13                                                   # return vel / np.sqrt(vel**2).sum(-1)

    14      2000        12762      6.4     15.7          return vel / np.sqrt(inner1d(vel, vel))

It reduces some time. However, I don't know how to proceed further. Because the number of loop may be very enormous, it will save a lot of time if reducing get_velocity a little time.

Could you please rewrite a Cython version for it?

Thanks,

在 2016年4月11日，下午10:17，Jason Furtney <jkfu...@gmail.com> 写道：

0.473

[Jason Furtney]

Dear Lei Pan,

Thanks for the reply and for your link, sorry for my late reply. Could
you put your latest version in a github repository so I can have a
look?

I will update the wiki page with your faster version.

You could reduce the number of calls to get_velocity by using lower
order integration. All the paths are independent so you could try the
multiprocessing module to calculate the paths in parallel if you have
many paths.

Have a look at the path finding algorithms in scikit-image:
http://scikit-image.org/docs/dev/api/skimage.graph.html
maybe these are better for your application?

Is the FMM based method faster than the scikit-image graph algorithms?
If you think there are applications for FMM based optimal path finding
maybe we should add it as a feature to scikit-fmm? It should be
possible to write a fast version of this algorithm in c or Cython.

Thanks,
Jason

> --
> You received this message because you are subscribed to the Google Groups
> "scikit-fmm" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to scikit-fmm+...@googlegroups.com.
> To post to this group, send email to sciki...@googlegroups.com.
> To view this discussion on the web visit

> https://groups.google.com/d/msgid/scikit-fmm/30255E38-B581-4F26-BBD8-97A9391512B5%40gmail.com.

[潘磊]

Thanks for your advises. I will look at the scikit-image. I put my codes on github,  this is the [address](https://github.com/panlei7/shortest_path)

You received this message because you are subscribed to a topic in the Google Groups "scikit-fmm" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/scikit-fmm/XvV_KhCUN5c/unsubscribe.
To unsubscribe from this group and all its topics, send an email to scikit-fmm+...@googlegroups.com.

To post to this group, send email to sciki...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/scikit-fmm/CAOU2ZSV9TGsg15Sch7smfJgh8%2BQ7pA-1q_LN47oc2wUFiKp8Kw%40mail.gmail.com.

[Jason Furtney]

Dear Lei Pan,

Thanks for posting your code on github. Are you still interested in
this many path finding problem? Have you made any improvements on your
example?

Thanks,
Jason

> https://groups.google.com/d/msgid/scikit-fmm/CANoUc%2BuYA2-%3D%3Dt3ZB52d%2BPNhFFBqUkr8RoJBZ8csUPx%3D3k1%3DJA%40mail.gmail.com.