Comments (10)
Hi and thanks for getting in touch with this bug, I have not seen this one before.
It looks like this issue is similar to #15. In this case, the discriminant of the quadratic solved to update a point is negative and this condition is not checked for properly. I added a test for this in 5a1e1db and an exception is thrown. I will look in more detail at why this is happening this coming week.
from scikit-fmm.
It looks like this is related to the second order point update. Doing the calculation in first order seems to work:
print skfmm.travel_time(phi, speed, order=1)
array([[ 1.5, 0.5, 8.77192982, 9.03360658, 23.39632103, 12.82051282, 0.5 ],
[ 1.5, 0.5, 1.7766502, 1.2254902, 0.50150451, 0.51388574, 0.5 ],
[ 2.04532893, 1.20710678, 0.5, 0.5, 0.5, 0.5, 1.20710678],
[ 2.75243571, 2.04532893, 1.5, 1.5, 1.5, 1.5, 2.04532893],
[ 3.54804305, 2.94223041, 2.5, 2.5, 2.5, 2.5, 2.94223041]])
I will keep digging.
from scikit-fmm.
OK, a little more information. It looks like this bug goes back to the first release of this package. At least now we know it was not introduced recently.
from scikit-fmm.
Here is a simpler version which shows the problem
import numpy as np
import skfmm
phi = np.array([[1, -1, -1],
[0, 0, -1],
[0, 0, 1]])
speed = np.array([[ 1, 0.01, 0.1],
[ 0, 0, 0.1],
[ 0, 0, 1. ]])
phi = np.ma.MaskedArray(phi,phi==0)
print skfmm.travel_time(phi,speed)
from scikit-fmm.
Below, when value = 0.0660411041104
, we get an answer. When value = 0.0660408040804
, we get the negative discriminant error message. The arrival time in cell [0,2] goes asymptotically to zero between these values (which is non-physical).
phi = np.array([[1, -1, -1],
[0, 0, -1],
[0, 0, 1]])
speed = np.array([[ 1, value, 0.1],
[ 0, 0, 1. ],
[ 0, 0, 1. ]])
phi = np.ma.MaskedArray(phi,phi==0)
print skfmm.travel_time(phi,speed)
This only happens when two orthogonal parts of two points of the second-order stencil cross the zero level-set in the travel_time
function. I am thinking it is not physically correct to construct the gradient this way. Maybe there is a different way to construct a second-order gradient or maybe we need to revert to first-order in these cases?
from scikit-fmm.
The first-order method has the same problem, just at different values. Around value = 0.03405940594059409
in this case.
from scikit-fmm.
The zero level set of phi going through a region with a high contrast in speed seems to cause this problem. It looks like this problem is happening when points adjacent to the initially frozen points are updated for the first time, before the main marching algorithm starts. In this simple example,
import numpy as np
import skfmm
phi = np.array([[1, -1, -1],
[0, 0, -1],
[0, 0, 1]])
speed = np.array([[ 1, value, 1],
[ 0, 0, 1],
[ 0, 0, 1]])
phi = np.ma.MaskedArray(phi,phi==0)
tt = skfmm.travel_time(phi,speed,order=order)
as value
decreases from 1 the arrival time at point (0,2) increase up to a local maximum and then goes asymptotically to zero. For second-order this local maximum occurs at value=1/2.
and for first-order value=1/3
. For second-order, below value=1/2.
the up-wind finite difference approximation along the x direction, used to estimate the travel-time gradient, falls below 1/2 which may be physically incorrect. For first order, below value=1/3.
the finite-difference approximation in the x direction falls below zero.
I still do not understand exactly what is happening. The solution may be to detect when the gradient at a point meets these conditions and abandon one of the directions of the stencil? I will be away from the computer until August 12th but can have a look when I get back.
from scikit-fmm.
Hi, have you made any further progress with this issue, or do you have any further advice how to handle it? I have recently updated to 2019.1.30 (from pip) from 0.0.9 and am now unable to use skfmm. Does this suggest the travel times I calculated with the previous version are likely to contain erroneous values? I'm happy to provide details of my specific use case if that would be useful. Thanks.
from scikit-fmm.
Thanks for the note. If you could share some examples it would help to track down the problem. I have not looked at this recently but I will try to find some time soon. (I am a graduate of your department :) A student from the University of Exeter may be working on this module over the summer and they may be looking at this issue. @wmoebius
from scikit-fmm.
OK, thanks to @f-fanni this is also now fixed in master.
from scikit-fmm.
Related Issues (20)
- Unable to use wheel with older versions of numpy HOT 2
- PyPI build seems to be broken HOT 4
- No wheels for python3.9 HOT 3
- Support for Python 3.10 HOT 7
- GIL is held while marching HOT 3
- Pip install broken on CIs HOT 4
- about cfmm HOT 1
- ValueError: the array phi contains no zero contour (no zero level set) HOT 1
- cannot import skfmm due to arm64 incompatibility? HOT 1
- Support for python 3.11 HOT 2
- Installation went well, but unable to import HOT 2
- Marching issue in a skeleton HOT 2
- numpy.distutils is no longer available in Python 3.12 HOT 28
- Inaccurate travel paths HOT 2
- Add a less formal introduction of FMM and what it can be used for in the README HOT 1
- error: Multiple top-level packages discovered in a flat-layout: ['skfmm', 'profile']. HOT 5
- Negative curvature of signed distance HOT 9
- a relevant way to cite scikit-fmm
- [Enhancement] Release GIL HOT 3
- Exact representation of contour HOT 7
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from scikit-fmm.