dx=[dY, dX]?

[Geordie McBain]

Hello.  I was momentarily confused by how skfmm.distance interprets array-like dx.  In two dimensions, it appears to use the elements of dx as the increments of Y and X, if those arrays are returned by numpy.meshgrid (as in the examples).  This can be demonstrated using a circular contour on an anisotropic grid, with indexing in ['ij', 'xy'] and flip in [False, True]:

    x, y = [np.linspace(0, 1, n) for n in [nx, ny]]
    X, Y = np.meshgrid(x, y, indexing=indexing)
    dx = [x[1] - x[0], y[1] - y[0]]

    phi = np.ones_like(X)
    phi[X*X + Y*Y < 0.2**2] = -1

    d = distance(phi, dx=dx[::-1] if flip else dx)

    fig, ax = pl.subplots()
    C = pl.contour(X, Y, d)

Geometrically correct behaviour is obtained by passing meshgrid the nondefault indexing='ij' option or passing distance the reversed dx=[dY, dX]. (But not both.)

Is this behaviour intentional?  Just thought I'd check before proposing a pull-request.
 
What is the intended interpretation of dx for three dimensions?  I think numpy.meshgrid treats the third index differently to the first two.

[Jason Furtney]

Dear Geordie McBain,

Thank you for your enthusiasm for this module and for your recent pull
request. Sorry for my delay in getting back to you.

I added the feature of allowing the grid spacing to be different in
different directions as an after-thought when I first wrote this
module. I agree this is confusing and if you have a proposal to
improve it or to clarify the documentation go ahead. I need to spend
some more time looking to answer your question properly. The original
intention was when dx is an array-like it represent the grid spacing
in the x, y and z directions. Why the x and y are flipped in the
meshgrid call was confusing when I was first looking at it.

I have not used this "different grid-spacing" feature in anger; it is
on my list to add more tests of this capability. There is currently
only one test of different grid-spacing in the test suite.

There was some recent discussion on the FiPy mailing list about a related topic:

http://thread.gmane.org/gmane.comp.python.fipy/3536/focus=3548

Some people have asked that scikit-fmm be adapted to work on
non-uniform (but still structured) grids (where the grid spacing is
different in each cell). As part of looking into this I am going to
verify the current implementation for different grid-spacing.

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/32a588da-ef0b-4a07-ab3f-169bd33d86ca%40googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.



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

[Geordie McBain]

Le 21 déc. 2014 13:23, "Jason Furtney" <jkfu...@gmail.com> a écrit :
>
> Dear Geordie McBain,
>
> Thank you for your enthusiasm for this module and for your recent pull
> request. Sorry for my delay in getting back to you.
>
> I added the feature of allowing the grid spacing to be different in
> different directions as an after-thought when I first wrote this
> module. I agree this is confusing and if you have a proposal to
> improve it or to clarify the documentation go ahead. I need to spend
> some more time looking to answer your question properly. The original
> intention was when dx is an array-like it represent the grid spacing
> in the x, y and z directions. Why the x and y are flipped in the
> meshgrid call was confusing when I was first looking at it.

Hello.  I looked into this a bit more on Friday and did a few experiments on anisotropic two-dimensional grids.  I realized that the intention of sfkmm is obvious (elements of dx correspond to dimensions), that its behaviour is correct, and that it's the default behaviour of numpy.meshgrid (doubtless under the unfortunate influence of the MATLAB function of the same name) that's confusing.  At least numpy.meshgrid offers indexing='xy' as an option.  Provided I used that, skfmm worked perfectly.
   Thus I don't think that the source code needs any modification (except for the future generalization that you mention to nonuniform as opposed to merely anisotropic spacing).  I propose to modify the examples in the documentation to use indexing='ij' where necessary, and possibly also using a nonsquare array, e.g. ones((4, 3)) instead of ones ((3, 3)), to emphasize that the elements of phi should be indexed by geometrical dimension rather than as if it were a linear-algebraic matrix. 
   I think that these considerations also mean that the output of distance should be transposed (and maybe flipped?) for printing if it is to be correctly oriented.  Of course the first example should still be kept as simple as possible; skfmm's clean simple  interface is a real plus.