Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Hotfix for second order update in travel time #35

Merged
merged 2 commits into from
May 24, 2020
Merged

Hotfix for second order update in travel time #35

merged 2 commits into from
May 24, 2020

Conversation

f-fanni
Copy link
Contributor

@f-fanni f-fanni commented Feb 4, 2020

According to many sources, starting from wikipedia Eikonal equation to some research papers and PhD thesis, in case the discriminant is less than 0 one should try the updates with one dimension less and select the minimum. This is not happening on master branch, and in travel-time-fix only the 1D dimension is used as fallback.
The code with this pull request does the dimensionality decrease in the generic n-D case, falling back to (n-1)-D, and so on.
This commit is only partial and has two main issues.

  1. It only fix the second order update, but it should apply equally to the first order update (I just could not find the code for it in the files).
  2. The code is potentially inefficient as the recursion can potentially try to solve all the (n)-D, (n-1)-D, ..., (1-D) dimensional updates; From my understanding, the (n-k-1)-D should be checked only if ALL (n-k)-D updates have discriminant < 0. This is not a big problem in 3D, but it could be quite slow with high dimensionality inputs.

From my tests this fixes the issues #28 and #18 .

@jkfurtney
Copy link
Member

Thanks for this, I will have a look as soon as I can.

@jkfurtney
Copy link
Member

Thanks for this, I needed a reminder to get back to this, sorry for the dead air from me. I will find some time in this quarantine life to work on this pull request and get an update out.

Can you explain the rationale of the GIL release? Within the FMM algorithm, I think the points need to be calculated in specific order to maintain causality. So multi-threading would not speed up this algorithm as far as I understand.

Does you GIL modification here have another effect? Does it allow more than one simultaneous call to these functions? If so, did you test it and is it faster?

@f-fanni
Copy link
Contributor Author

f-fanni commented May 23, 2020

Hi, it's completely unrelated to the previous fix, I pushed on the wrong branch which updated the pull request. Sorry about that.

Anyway, the GIL release is not to parallelize a single call to the marcher, which is indeed not parallelizable, but is useful to allow n threads to call the library concurrently. Here a small example of what I mean, https://pastebin.com/ZG7Jybey. With this commit the threaded and the multiprocess versions take a similar amount of time (on my machine 10s v 30s for the serial calls), while your version should perform better only on multiprocess. (Can you please confirm this? I wrote this specific test after the commit and didn't test it on your version)

The real issue with GIL release is that we should not touch the python resources afterwards, as it could cause race condition or similar problems. Though, from what I saw in the code, this should not happen (the python variables passed as args to the functions are used as read only, and I release the GIL after the result allocation), and the extremely limited testing I've done seems to work just fine.

Again, I didn't intend to update the pull request, both because it is unreleated to the pull request itself, and because some more testing is probably necessary and I am not 100% sure of the inherent correctness of releasing the GIL.

@jkfurtney
Copy link
Member

OK, what you say makes sense, thanks for the detailed reply. I agree that it should be OK to call from multiple threads simultaneously. I will look at the original changes now and I can get to checking the threading situation next.

@jkfurtney jkfurtney merged commit 3e54d7f into scikit-fmm:travel-time-fix May 24, 2020
@jkfurtney
Copy link
Member

OK, thanks again for the time you put into fixing this and sorry that I sat on it for so long. I merged this into the branch then into master. I reverted the GIL stuff for now. It would be great if you could create another PR for that. I am glad to cross these issue of after many years. I will ask a few other people to test and I will get a new release out.

@f-fanni f-fanni mentioned this pull request May 25, 2020
Closed
@aleccarruthers
Copy link

I am still getting a negative discriminant error with the 'travel_time' function. Has there been any recent attempts to address this issue?

@aleccarruthers
Copy link

aleccarruthers commented Jun 4, 2020 •
edited

#phi = np.nan * np.ones((reciprocalLocalCostArray.shape))
phi = np.zeros((reciprocalLocalCostArray.shape))
#speed = np.ones((reciprocalLocalCostArray.shape))* np.nan
speed = np.zeros((reciprocalLocalCostArray.shape))
phi[maskedBasinFAC!=0] = 1
print(f'phi: {phi[~np.isnan(phi)]}')
speed[maskedBasinFAC!=0] = reciprocalLocalCostArray[maskedBasinFAC!=0]
print(np.max(speed[~np.isnan(speed)]))
phi[fastMarchingStartPointListFMM[0,i],fastMarchingStartPointListFMM[1,i]] = -1
print(len(phi[~np.isnan(phi)]))
try:
travelTimearray = skfmm.travel_time(phi, speed, dx=1)

Traceback (most recent call last): File "pygeonet_fast_marching.py", line 215, in <module> main() File "pygeonet_fast_marching.py", line 208, in main geodesicDistanceArray = Fast_Marching(fastMarchingStartPointListFMM, basinIndexArray, flowArray, reciprocalLocalCostArray) File "pygeonet_fast_marching.py", line 146, in Fast_Marching travelTimearray = skfmm.travel_time(phi, speed, dx=1) File "C:\Users\alecc\Anaconda3\envs\Spatial\lib\site-packages\skfmm\pfmm.py", line 173, in travel_time int(self_test), TRAVEL_TIME, order, narrow, periodic) RuntimeError: Negative discriminant in time marcher quadratic.

@jkfurtney
Copy link
Member

Hi Alec, thanks for the message, what version are you using? Have you tried the most recent master from github? We just put out some changes that may affect this behavior.

@aleccarruthers
Copy link

aleccarruthers commented Jun 4, 2020 via email

@jkfurtney
Copy link
Member

To test the latest you would need to build the module from source. I am planning to push a new release out soon and the anaconda packages should update automatically. In the example you posted I cannot run it to see what happens. Can you create a stand alone example with everything to reproduce the problem? I can test the latest version to see if it fixes the problem.

@aleccarruthers
Copy link

aleccarruthers commented Jun 4, 2020 via email

@aleccarruthers
Copy link

aleccarruthers commented Jun 4, 2020 via email

@jkfurtney
Copy link
Member

Hi, there are build instructions in the github readme. Good point that there is nothing in the docs about actually building the code. I will update that.

@aleccarruthers
Copy link

aleccarruthers commented Jun 5, 2020 via email

@jkfurtney
Copy link
Member

Hard to say without seeing any results. Can you post an example along with the MATLAB results?

@wmoebius
Copy link
Contributor

@aleccarruthers Might be able to help with the comparison to MATLAB. Found good agreement between scikit-fmm (this) and imsegfmm (MATLAB).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers
No reviews
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

Successfully merging this pull request may close these issues.

None yet
4 participants
@f-fanni @jkfurtney @aleccarruthers @wmoebius