Added optional keyword argument "order" and fixed some glitches in travel_time()

[Jason Furtney]

Daniel,

I have added the optional keyword argument "order" to distance(),
travel_time() and extension_velocities() for testing. The default is
2.

travel_time() (the Eikonal equation case) can use the second order
stencil in initializing the narrow band now.

I added tests based on the email you sent me on June 14th and
scikit-fmm is giving the same results as lsmlib for order 1 and 2 for
both distance() and travel_time(). Tests here:
https://github.com/scikit-fmm/scikit-fmm/blob/master/tests/test_dw.py

One of my original tests is still not passing, when I resolve that I
will release v0.0.3

Thanks again for all your help,
Jason

[Daniel Wheeler]

Thanks Jason!

--
You received this message because you are subscribed to the Google Groups "scikit-fmm" group.
To post to this group, send email to sciki...@googlegroups.com.
To unsubscribe from this group, send email to scikit-fmm+...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/scikit-fmm?hl=en.

--
Daniel Wheeler

[Daniel Wheeler]

Hi Jason,

I hope all is well. I finally got around to working with Scikit-fmm
again. I'm trying to use the "extension_velocities" functionality. I
noticed the following.

>>> import skfmm
>>> skfmm.extension_velocities([-1, -1, 1, 1], [0, 2, 1, 0])
(array([-1.5, -0.5, 0.5, 1.5]), array([ 1., 1., 2., 2.]))

This is different from lsmlib, which gives

>>> import pylsmlib
>>> pylsmlib.computeExtensionFields([-1, -1, 1, 1], [0, 2, 1, 0])
(array([-1.5, -0.5, 0.5, 1.5]), array([ 1.5, 1.5, 1.5, 1.5]))

It seems to me that the lsmlib behaviour is preferable. It uses the
value averaged at the zero level set to calculate the extension
velocities. What are your thoughts on this? If you agree that this
needs fixing, I can give it a shot if you like. Basically, the cell
adjacent to the interface is only using the neighbour cell to evaluate
it's extension value.

<https://github.com/scikit-fmm/scikit-fmm/blob/master/skfmm/extension_velocity_marcher.cpp#L49>

In LSMLIB the value is calculated as a combination of the neighbor and
the current cell's value.

<https://github.com/wd15/LSMLIB/blob/fipy/src/serial/lsm_FMM_field_extension.c#L546>

Anyway, let me know what you think. I will go ahead and modify it and
get a patch back to you. I also want to add in an extension mask that
only uses some of the interface values to evaluate the initial
guess. I'll submit a pull request when I'm done.

I've been working a bit on Pylsmlib. I have made the interface exactly
like Scikit-fmm and ported the fipy test suite over and I can do the
same to Scikit-fmm as well. They can then both have compatible
interfaces and test suites, which I think is nice.

Cheers

On Fri, Jul 13, 2012 at 5:50 PM, Jason Furtney <jkfu...@gmail.com> wrote:

[Jason Furtney]

Daniel,

Thanks for the reply and for your continued efforts on scikit-fmm.
Yes, please have a crack at improving the extension velocities.

I have done some experimenting on speeding the module up. Some simple
loop unrolling in the point update functions gave a 15% speed-up, much
to my surprise. I think some additional speed-up is possible if we
create a specialized version of the point update routines for 2D
arrays (or for 3D arrays).

Are you specifically interested in 2D analysis or do you run cases in 3D also?

Thanks,
Jason

> For more options, visit https://groups.google.com/groups/opt_out.
>
>



--
--
Jason K. Furtney
Itasca Consulting Group
111 3rd Ave. South, Suite 450
Minneapolis, MN 55401 USA
(612) 371-4711

[Daniel Wheeler]

On Tue, Sep 18, 2012 at 5:07 PM, Jason Furtney <jkfu...@gmail.com> wrote:
> Daniel,
>

> Thanks for the reply and for your continued efforts on scikit-fmm.
> Yes, please have a crack at improving the extension velocities.

Okay. I'll have a go.

> I have done some experimenting on speeding the module up. Some simple
> loop unrolling in the point update functions gave a 15% speed-up, much
> to my surprise. I think some additional speed-up is possible if we
> create a specialized version of the point update routines for 2D
> arrays (or for 3D arrays).
>
> Are you specifically interested in 2D analysis or do you run cases in 3D also?

Currently, FiPy only has 2D examples, but I need 3D for my own
research. However, the actual time spent doing level set related
calculations is minimal compared with solving advection diffusion type
equations at least in 2D, but any additional improvements can't hurt.

--
Daniel Wheeler

[Daniel Wheeler]

On Tue, Sep 18, 2012 at 5:07 PM, Jason Furtney <jkfu...@gmail.com> wrote:
> Daniel,
>

> Thanks for the reply and for your continued efforts on scikit-fmm.
> Yes, please have a crack at improving the extension velocities.

I figured out a couple of issues and things seem to be working better,
but I haven't pushed to github yet. Do you mind if I clean up the test
harness? It might be nice to run it with a single command "python
setup.py test" for example or "import skfmm; skfmm.test()". I prefer
the latter.

Also, I am planning on adding a number of additional tests. How would
you feel about shifting the tests over to the doctest format? I can
push some of the tests to appropriate functions were they are
descriptive and include a _test() for miscellaneous tests. I really
favour doctest over the alternatives.

Cheers

--
Daniel Wheeler

[Jason Furtney]

Yes, please go ahead and clean-up the test system. Could you change
the section on testing in the README.txt if you make any major
changes? Thanks for your contribution. After I merge your branch in I
will have a try at some speed-ups. Do you have any thoughts on the
difficulty of, or the value of converting this library to work with
Python 3?

Jason

[Daniel Wheeler]

On Wed, Sep 19, 2012 at 6:53 PM, Jason Furtney <jkfu...@gmail.com> wrote:
> Yes, please go ahead and clean-up the test system. Could you change
> the section on testing in the README.txt if you make any major
> changes?

Will do.

> Thanks for your contribution. After I merge your branch in I
> will have a try at some speed-ups. Do you have any thoughts on the
> difficulty of, or the value of converting this library to work with
> Python 3?

It'll have to be done at some point. I think the community will move
to 3 on mass when a linux or mac distro starts shipping with it as the
default. My colleague switched FiPy over to work with Python 3 with
relative ease, but it doesn't have C++ modules. I've CCed him to see
if he has any thoughts on this.

--
Daniel Wheeler

[Daniel Wheeler]

Jason,

I submitted a pull request for the changes I've made. See what you think.

Cheers

P.S. Do you try and avoid fast forward merges?

Daniel Wheeler

[Jason Furtney]

Daniel,

Thanks for your efforts on this. The extension velocity example
picture looks much better now! I dimly wondered why the lines of
constant velocity were not perpendicular to the interface.. All the
tests pass on windows and linux. I will test on OS X and push a new
version to pypi soon.

Jason

On Mon, Oct 15, 2012 at 1:57 PM, Daniel Wheeler

[Daniel Wheeler]

Nice one. Thanks.

[Daniel Wheeler]

One other thing. I noticed an issue that occurs in both LSMLIB and
Sckit-fmm that could adversely affect second order accuracy. I am not
100% certain I fully understand it, but I want to try and fix the
issue in both codes but wanted to discuss it first. See
<https://github.com/ktchu/LSMLIB/issues/2>. I've CCed the author of
LSMLIB as well as he might have some insights into this.

[Daniel Wheeler]

Here is the Skfmm bug report

<https://github.com/scikit-fmm/scikit-fmm/issues/5>

On Mon, Oct 15, 2012 at 5:06 PM, Daniel Wheeler

--
Daniel Wheeler

[Jason Furtney]

Daniel and Kevin,

Thanks for the note. I looked at this and am not 100% sure either.
I basically implemented what is described here:
http://www2.imm.dtu.dk/pubdb/views/edoc_download.php/841/pdf/imm841.pdf
which does not mention this subtlety.

I experimented with a little patch to update the far point in the
second order stencil. This is not 100% correct because it causes some
of the tests to fail. I will give this some more thought.

patch below.

From f7f5c35868ca4caa8edb5b8df23e46765eeb2ef6 Mon Sep 17 00:00:00 2001
From: Jason Furtney <jkfu...@gmail.com>
Date: Tue, 16 Oct 2012 09:35:58 -0500
Subject: [PATCH] testing second order fix

---
skfmm/base_marcher.cpp | 32 ++++++++++++++++++++++++++++----
1 file changed, 28 insertions(+), 4 deletions(-)

diff --git a/skfmm/base_marcher.cpp b/skfmm/base_marcher.cpp
index 9b7913c..45d5aea 100644
--- a/skfmm/base_marcher.cpp
+++ b/skfmm/base_marcher.cpp
@@ -158,11 +158,35 @@ void baseMarcher::solve()
}
}
}
- }
- }
- }
+ //==========================================================
+ // update the far point in the second order stencil
+ // "jump" over a Frozen point if needed
+ if (order_ == 2)
+ {
+ int naddr2 = _getN(addr,dim,j*2,Frozen);
+ if (naddr2!=-1 && flag_[naddr2]!=Frozen)
+ {
+ double d = updatePointOrderTwo(naddr2);
+ if (d && flag_[naddr2]==Narrow)
+ {
+ heap_->set(heapptr_[naddr2], fabs(d));
+ distance_[naddr2]=d;
+ }
+ else if (d && flag_[naddr2]==Far)
+ {
+ distance_[naddr2]=d;
+ flag_[naddr2]=Narrow;
+ heapptr_[naddr2] = heap_->push(naddr2,fabs(d));
+ }
+ }
+ }
+ //==========================================================
+ } // for each direction
+ } // for each dimension
+ } // main loop of Fast Marching Method
+
// add back mask here. The python wrapper will look for elements
- // equal to mexDouble and add the mask back
+ // equal to maxDouble and add the mask back
for (int i=0; i<size_; i++)
{
if (flag_[i] == Mask) distance_[i] = maxDouble;
--
1.7.11.msysgit.1

On Mon, Oct 15, 2012 at 4:06 PM, Daniel Wheeler

Jason K. Furtney, PhD

Itasca Consulting Group
111 3rd Ave. South, Suite 450
Minneapolis, MN 55401 USA
(612) 371-4711

www.itascacg.com

[Daniel Wheeler]

Hi Jason,

I've just pushed a branch to my git repository

<https://github.com/wd15/scikit-fmm/tree/second-order>

I applied your patch and modified it somewhat and also included a test
case that better highlights the issue. Using the modified patch, the
tests pass other than the new test, so obviously the problem hasn't
been fixed. Anyway, I am not sure whether it is possible to fix this.
The test I added is as follows:

>>> phi = np.array([[-1, 1, 1, 1, 1, -1],
... [-1, -1, -1, -1, -1, -1],
... [-1, -1, -1, -1, -1, -1]])

>>> phi = distance(phi)
>>> print phi[2, 2] == phi[2, 3]

One would expect the answer to be symmetric about the second axis, but
it evidently is not when using second order. I find this confusing.
Maybe you have some insight into this.

Cheers

[Jason Furtney]

Daniel,

Thanks again for your efforts, this is an interesting example.

I think I understand what you mean now. I merged your second-order
branch into the trunk and applied some changes, this extra loop should
be doing what it is supposed to do now.

Maybe there are two separate issues? On the symmetry problem: it seems
like the symmetry is broken when two points with the same trial
distance are processed in arbitrary order. I experimented with
freezing batches of points which are equidistant from the front. The
solution becomes symmetric and all the tests still pass. As far as I
can tell, this is OK because it still respects the causality of a
monotonically expanding front. Maybe there are some subtle reasons
this is not a good idea? There are possibly some floating point
precision issues.

This experiment and some logging code are in a branch on GitHub:
https://github.com/scikit-fmm/scikit-fmm/tree/logging

I wonder if this second-order problem is just inherent to the method?
My understanding is the first-order stencil does not have this
problem. If so, it might be worth looking into Multi-Stencil FMM. The
stencils are all first-order but the accuracy is better than the
vanilla second-order scheme:

http://www.cvip.uofl.edu/wwwcvip/research/publications/Pub_Pdf/2007_2/Sabry_Farag_TPAMI_2007.pdf

Let me know what you think. Below is some of the logging output for
your symmetry example. Thanks.

=============================================
Logging output from the first commit on the logging branch: ad81a
=========================================


initializing beginning narrow-band
updating: (1, 0)
old value 0.000000
new value -1.054700
updating: (1, 5)
old value 0.000000
new value -1.054700
updating: (2, 1)
old value 0.000000
new value -1.451184
updating: (2, 2)
old value 0.000000
new value -1.500000
updating: (2, 3)
old value 0.000000
new value -1.500000
updating: (2, 4)
old value 0.000000
new value -1.451184

starting solve loop

FFFFFF
NFFFFN
0NNNN0
freezing: (1, 0)
final value: -1.054700

updating: (2, 0)
old value 0.000000
new value -1.906267

FFFFFF
FFFFFN
NNNNN0
freezing: (1, 5)
final value: -1.054700

updating: (2, 5)
old value 0.000000
new value -1.906267

FFFFFF
FFFFFF
NNNNNN
freezing: (2, 1)
final value: -1.451184

updating: (2, 0)
old value -1.906267
new value -1.850740
updating: (2, 2)
old value -1.500000
new value -1.499230

FFFFFF
FFFFFF
NFNNNN
freezing: (2, 4)
final value: -1.451184

updating: (2, 3)
old value -1.500000
new value -1.499230
updating: (2, 5)
old value -1.906267
new value -1.850740

FFFFFF
FFFFFF
NFNNFN
freezing: (2, 2)
final value: -1.499230

updating: (2, 0)
old value -1.850740
new value -1.850740
updating: (2, 3)
old value -1.499230
new value -1.499822

FFFFFF
FFFFFF
NFFNFN
freezing: (2, 3)
final value: -1.499822

updating: (2, 5)
old value -1.850740
new value -1.850740

FFFFFF
FFFFFF
NFFFFN
freezing: (2, 5)
final value: -1.850740


FFFFFF
FFFFFF
NFFFFF
freezing: (2, 0)
final value: -1.850740

[[-0.5 0.35355339 0.5 0.5 0.35355339 -0.5 ]
[-1.0547002 -0.5 -0.5 -0.5 -0.5 -1.0547002 ]
[-1.85073968 -1.45118446 -1.49923009 -1.49982156 -1.45118446 -1.85073968]]
-1.49923009455 -1.49982155712

=============================================
Logging output from the second commit on the logging branch: 2b2607
=========================================


initializing beginning narrow-band
updating: (1, 0)
old value 0.000000
new value -1.054700
updating: (1, 5)
old value 0.000000
new value -1.054700
updating: (2, 1)
old value 0.000000
new value -1.451184
updating: (2, 2)
old value 0.000000
new value -1.500000
updating: (2, 3)
old value 0.000000
new value -1.500000
updating: (2, 4)
old value 0.000000
new value -1.451184

starting solve loop

FFFFFF
NFFFFN
0NNNN0
freezing: (1, 0)
final value: -1.054700


freezing: (1, 5)
final value: -1.054700

updating: (2, 0)
old value 0.000000
new value -1.906267
updating: (2, 5)
old value 0.000000
new value -1.906267

FFFFFF
FFFFFF
NNNNNN
freezing: (2, 1)
final value: -1.451184


freezing: (2, 4)
final value: -1.451184

updating: (2, 0)
old value -1.906267
new value -1.850740
updating: (2, 2)
old value -1.500000
new value -1.499230
updating: (2, 3)
old value -1.500000
new value -1.499230
updating: (2, 5)
old value -1.906267
new value -1.850740

FFFFFF
FFFFFF
NFNNFN
freezing: (2, 2)
final value: -1.499230


freezing: (2, 3)
final value: -1.499230

updating: (2, 0)
old value -1.850740
new value -1.850740
updating: (2, 5)
old value -1.850740
new value -1.850740

FFFFFF
FFFFFF
NFFFFN
freezing: (2, 5)
final value: -1.850740


freezing: (2, 0)
final value: -1.850740

[[-0.5 0.35355339 0.5 0.5 0.35355339 -0.5 ]
[-1.0547002 -0.5 -0.5 -0.5 -0.5 -1.0547002 ]
[-1.85073968 -1.45118446 -1.49923009 -1.49923009 -1.45118446 -1.85073968]]
-1.49923009455 -1.49923009455

[Daniel Wheeler]

On Sun, Oct 21, 2012 at 10:32 PM, Jason Furtney <jkfu...@gmail.com> wrote:
> Daniel,
>

> Thanks again for your efforts, this is an interesting example.
>
> I think I understand what you mean now. I merged your second-order
> branch into the trunk and applied some changes, this extra loop should
> be doing what it is supposed to do now.

I'm not sure it fixes anything, but I'll let you be the arbiter.

> Maybe there are two separate issues?

I think you're right.

> On the symmetry problem: it seems
> like the symmetry is broken when two points with the same trial
> distance are processed in arbitrary order.

Obviously, there is nothing to be done about that.

> I experimented with
> freezing batches of points which are equidistant from the front. The
> solution becomes symmetric and all the tests still pass. As far as I
> can tell, this is OK because it still respects the causality of a
> monotonically expanding front. Maybe there are some subtle reasons
> this is not a good idea? There are possibly some floating point
> precision issues.

What about not using information from a location that is equal to your
value. I actually tried messing with <= and < and didn't seem to help,
but I wasn't very systematic.

> This experiment and some logging code are in a branch on GitHub:
> https://github.com/scikit-fmm/scikit-fmm/tree/logging

I'll put this on my list of stuff to look at and get back to you.

> I wonder if this second-order problem is just inherent to the method?

It seems to me that if you only update using values that are less than
and update the whole stencil then there should be no asymmetries, but
I have to work through you're stuff below carefully to understand it.

I'll take a look at what you've done in the next few weeks. I have
some deadlines coming up and I'm going on vacation too.

Thanks for your effort!

--
Daniel Wheeler