Use fused types in denoise, warp by jni · Pull Request #3486 · scikit-image/scikit-image

jni · 2018-10-19T08:06:47Z

Description

Supersedes #3469, #3253 (I think), #1977.

All credit to @AetherUnbound and @hmaarrfk here, I have just rebased their work on the latest master after @hmaarrfk fixed our negative indexing use in Cython.

Checklist

Clean style in the spirit of PEP8
Docstrings for all functions
~~Gallery example in ./doc/examples (new features only)~~
Benchmark in ./benchmarks, if your changes aren't covered by an
existing benchmark
Unit tests

For reviewers

Check that the PR title is short, concise, and will make sense 1 year
later.
Check that new functions are imported in corresponding __init__.py.
Check that new features, API changes, and deprecations are mentioned in
doc/release/release_dev.rst.
Consider backporting the PR with @meeseeksdev backport to v0.14.x

pep8speaks · 2018-10-19T08:06:56Z

Hello @jni! Thanks for updating this PR. We checked the lines you've touched for PEP 8 issues, and found:

In the file skimage/restoration/_denoise.py:

Line 18:68: E226 missing whitespace around arithmetic operator

Comment last updated at 2019-03-05 23:40:45 UTC

hmaarrfk · 2018-10-19T10:55:04Z

Why is OSX failing? https://travis-ci.org/scikit-image/scikit-image/jobs/443575582#L5300

AetherUnbound

Looks good to me - thanks for compiling this all together! And thanks @hmaarrfk for fixing the negative indexing stuff.

hmaarrfk · 2018-10-19T14:45:54Z

@jni, I'm going to say that though I appreciate that my original PR is taken seriously, I don't know how worthwhile it is to complicate the warping code if it doesn't really give a speedup. In my experience, templated code is so much harder to read and to adapt. By including this, we are forcing future developers to also develop warping algorithms that function for all numeric types. This is a huge potential burden.

I would look into why it was that any speedup was obtained for float64 types and maybe just apply that fix. It might have been passing a parameter by reference in a tight loop which was necessary to template the functions.

I think the assumption that integer arithmetic is faster than floatingpoint is seriously flawed due to massive hardware optimizations for the floating point units in modern processors, though I'm not an expert in this field.

float32/float64 seems to be a valid assumption if you can guarantee not upcasting. float16 only recently started to exist because GPUs are running out of memory for machine learning.

I can't speak specifically for @AetherUnbound's PR, I was mostly helping with syntax.

jni · 2018-10-19T23:49:22Z

@hmaarrfk but speed is not the only advantage, there is also RAM usage, right? I think that's very valuable. And personally I don't find the code much more complicated. And you know I'm quite nitpicky about complicated code! =D

It looks like the test failure is to do with random noise values. Same failure on that Mac build as on AppVeyor:
https://travis-ci.org/scikit-image/scikit-image/jobs/443575582#L5285-L5301

Edit: here's the failure for quick reference:

__________________________ test_denoise_bilateral_2d ___________________________
    def test_denoise_bilateral_2d():
        img = checkerboard_gray.copy()[:50, :50]
        # add some random noise
        img += 0.5 * img.std() * np.random.rand(*img.shape)
        img = np.clip(img, 0, 1)
    
        out1 = restoration.denoise_bilateral(img, sigma_color=0.1,
                                             sigma_spatial=10, multichannel=False)
        out2 = restoration.denoise_bilateral(img, sigma_color=0.2,
                                             sigma_spatial=20, multichannel=False)
    
        # make sure noise is reduced in the checkerboard cells
        assert_(img[30:45, 5:15].std() > out1[30:45, 5:15].std())
>       assert_(out1[30:45, 5:15].std() > out2[30:45, 5:15].std())
E       AssertionError
skimage/restoration/tests/test_denoise.py:183: AssertionError

Dunno why it would happen in this PR and not elsewhere, since that test appears to use float64 anyway... Any ideas?

hmaarrfk · 2018-10-19T23:51:03Z

Probably because we added some tests, the seed is different than before. Numpy (the repo) actually prints the values of the arrays that failed.

Would be useful if we had that....

jni · 2018-10-20T00:00:49Z

@hmaarrfk it's not that the seed is different, it's that it's not set. You can run

pytest skimage/restoration/tests/test_denoise.py::test_denoise_bilateral_2d

repeatedly to get different results each time.

hmaarrfk · 2018-10-20T00:04:47Z

@jni really? i thought the python seed was set to some consistent value everytime. I guess this isn't matlab again. It probably needs an environment variable.

jni · 2018-10-20T00:24:39Z

Anyway, perhaps it's not a good test, but it's a bit concerning that this PR brought that out.

hmaarrfk · 2018-10-20T04:21:52Z

The very first time I ran it, it failed, then all other times, it passed. I even installed pytest-repeat and ran it with

➤ pytest skimage/restoration/tests/test_denoise.py::test_denoise_bilateral_2d --count 100

still passed....

hmaarrfk · 2018-10-20T04:36:23Z

export PYTHONHASHSEED=1   
pytest skimage/restoration/tests/test_denoise.py::test_denoise_bilateral_2d

will fail, but I can't run it with --pdb

hmaarrfk · 2018-10-20T04:39:02Z

    def test_denoise_bilateral_2d():
        img = checkerboard_gray.copy()[:50, :50]
        # add some random noise
        noise = 0.5 * img.std() * np.random.rand(*img.shape)
        img += noise
        img = np.clip(img, 0, 1)
    
        out1 = restoration.denoise_bilateral(img, sigma_color=0.1,
                                             sigma_spatial=10, multichannel=False)
        out2 = restoration.denoise_bilateral(img, sigma_color=0.2,
                                             sigma_spatial=20, multichannel=False)
    
        # make sure noise is reduced in the checkerboard cells
>       assert img[30:45, 5:15].std() > out1[30:45, 5:15].std()
E       assert 0.07044972960234154 > nan
E        +  where 0.07044972960234154 = <built-in method std of numpy.ndarray object at 0x7fcbe9b3ce40>()
E        +    where <built-in method std of numpy.ndarray object at 0x7fcbe9b3ce40> = array([[ 0.12671524,  0.15149392,  0.21206663,  0.15423946,  0.00545805,\n         0.16724259,  0.21282153,  0.1807636 ...1418,  0.08659465,  0.017158  ,  0.18889751,\n         0.09346233,  0.08582789,  0.12415705,  0.2002917 ,  0.07256701]]).std
E        +  and   nan = <built-in method std of numpy.ndarray object at 0x7fcbe9b3cdf0>()
E        +    where <built-in method std of numpy.ndarray object at 0x7fcbe9b3cdf0> = array([[ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan],\n       [ nan,  nan,  nan,  nan,  nan,  nan,  nan,...  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan],\n       [ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan]]).std

skimage/restoration/_denoise.py

skimage/restoration/_denoise_cy.pyx

hmaarrfk · 2018-10-20T12:03:56Z

skimage/restoration/_denoise_cy.pyx

    cdef:
-        double[:, :, ::1] cimage = np.ascontiguousarray(image)
-        double[:, :, ::1] cu = u
+        np_floats[:, :, ::1] cu = u


Can cu be removed? u is already a horrible name

hmaarrfk · 2018-10-20T12:05:20Z

Do we have a benchmark? is it possible to see what happens when we change the looping order of the tight loop in _denoise_tv_bregman to C ordering?

skimage/restoration/_denoise_cy.pyx

skimage/restoration/_denoise.py

skimage/restoration/_denoise_cy.pyx

hmaarrfk · 2018-10-20T13:34:36Z

Finally, when all those comments are addressed, you can release the gil. You will need this line of code somewhere

                      color_lut_bin = min(
                          <Py_ssize_t>(dist * dist_scale), max_color_lut_bin)

hmaarrfk · 2018-10-20T13:40:59Z

skimage/restoration/_denoise.py

+    max_value /= bins
+
+    for b in range(bins):
+        color_lut[b] = _gaussian_weight(sigma, b * max_value)


Python code should use numpy and not cythonisms.

hmaarrfk · 2018-10-20T14:02:28Z

@jni, thanks for helping find the bug.

Can you please close this PR? it seems there is more work to do in the denoising PR.

jni · 2018-10-21T01:54:00Z

@hmaarrfk I don't see why the debugging needs to happens specifically in one or the other PR? I'm happy to take this on. Thanks for the detailed and specific review either way!

hmaarrfk · 2018-10-21T02:11:17Z

whatever works for you and @AetherUnbound \o/

skimage/_shared/interpolation.pxd

jni · 2018-11-29T14:01:15Z

@AetherUnbound I've given this another go. Sorry to miss your ping earlier this month. I think @hmaarrfk is on holiday in Japan atm. =D

Let's see how my latest commits fare in CI! =)

jni · 2018-11-29T14:01:55Z

One last thing I want to do is move all of the LUT code out of Cython. But for now it's bedtime.

jni · 2018-12-01T03:02:33Z

It may not look like it, but CI is actually happy! 😝 🎉 🎉 🎉

@scikit-image/core can we have some reviews please? =D This fixes a regression in 0.14.1 relative to 0.14.0 so it's relatively urgent (within scikit-image levels of urgency, anyway =P).

The failures are in --pre because, as far as I can tell, NumPy 1.16 changed the deprecation message for numpy.matrix, so all the work in #3242 is now obsolete. 🤦‍♂️ I'll raise a separate issue but I don't think it should hold up this PR. [edit: @sciunto already found it in #3571]

https://travis-ci.org/scikit-image/scikit-image/jobs/461638460#L3340

ValueError: Unexpected warning: the matrix subclass is not the recommended way to represent matrices or deal with linear algebra (see https://docs.scipy.org/doc/numpy/user/numpy-for-matlab-users.html). Please adjust your code to use regular ndarray.

AetherUnbound

Looks good to me! Beyond that test comment change

skimage/_shared/interpolation.pxd

skimage/restoration/_denoise_cy.pyx

skimage/restoration/tests/test_denoise.py

AetherUnbound · 2018-12-03T20:52:56Z

Ah, I don't have permission to push to your branch 😜

jni · 2018-12-04T07:18:19Z

@AetherUnbound you don't have permission, but for future reference, you can always make a PR to my fork/branch.

jni · 2018-12-04T07:29:04Z

Also, you can probably use the "suggest" button when making comments on code...?

jni · 2019-03-05T00:37:35Z

Also, thank you very much for running the benchmarks! I have a visitor this week so I have been focusing on napari. But I agree with you that it would be excellent to push both 0.14.3 and 0.15 out the door. I think it's feasible. If we get this merged I'll devote some time today for release notes!

hmaarrfk · 2019-03-05T00:59:28Z

lets just say that I used to "wait for scikit image to comile. Now I go get a coffee.

jni · 2019-03-05T01:31:05Z

@hmaarrfk this sounds like a win-win! 😂

stefanv · 2019-03-05T07:17:31Z

This changes the return from skimage.transform.warp substantially. I am concerned about returning ints from an interpolation operation: there is significant precision loss there. I can imagine a dtype argument to determine the output, but it seems dangerous to presume the user always wants the same dtype out as in (yes, I know ndimage does this).

Either way, this behavior change is big enough that it would require deprecation if it were to go through.

jni · 2019-03-05T08:02:59Z

@stefanv I think you've misread the PR. The output from warp is determined here:

scikit-image/skimage/transform/_warps_cy.pyx

Line 138 in f9046b4

cdef double[:, ::1] out = np.zeros((out_r, out_c), dtype=np.double)

and this line is not touched. As far as the warps module is concerned, and as far as I can tell, only internal functions are touched — which might explain the benchmark results 😂.

stefanv · 2019-03-05T08:11:14Z

OK, so maybe I also missed the intent / benefit of this PR then. Do you mind summarizing?

jni · 2019-03-05T13:43:13Z

@stefanv the most immediate benefit is that it fixes #3449, a bug that was introduced in 0.14.1 that is preventing CellProfiler from upgrading past 0.14.0. The warping changes lay the groundwork for allowing float32 warps. Really there wouldn't be any changes to warp but we needed fused types to fix the denoise error and @hmaarrfk had already implemented them in #3253, so we built off that.

stefanv · 2019-03-05T18:55:41Z

OK, I am +1; @jni @hmaarrfk do you mind confirming that all inline comments have been addressed? E.g., the pointer[0] syntax should probably be pointer*, but otherwise everything looks good to me.

skimage/_shared/interpolation.pxd

... Hold on to your butts! =P Co-Authored-By: jni <juan.nunez-iglesias@monash.edu>

hmaarrfk · 2019-03-05T23:42:56Z

I think everything is in order @stefanv

This reverts commit f0d48db.

* MNT: Add a fused numeric type to make fused_types more constent. * MNT: make the interpolation use fused types for any->any type interpolation of images. * MNT: Move even more numpy function calls from Cython to Python. * MNT: More explicit type specifying * BUG: Use fused floats in denoise Pull all denoise array allocation out into python Function used in python should be def not cdef Have the correct number of arguments * MNT: Add additional type check for denoise, more docs * PEP8 fixes * Remove redundant array u and view cu * Mutate output array in Cython and trim in Python * Replace height and width with existing rows and cols * Change iteration order to match array order * Move range_lut and color_lut properly from Cy to Py * Ravel range LUT which is expected to be 1D * Update comment in denoise bilateral tests * Fix relative import in interpolation.pyx * BENCH: Benchmark for warping with many types * Rename warp benchmark file * Update pointer syntax Co-Authored-By: Mark Harfouche <mark.harfouche@gmail.com>

stefanv · 2019-03-06T01:38:04Z

Don't ask :)

jni · 2019-03-06T06:50:12Z

I think there's very little chance that this will work but...

@meeseeksdev backport to v0.14.x

lumberbot-app · 2019-03-06T06:50:19Z

Owee, I'm MrMeeseeks, Look at me.

There seem to be a conflict, please backport manually. Here are approximate instructions:

Checkout backport branch and update it.

$ git checkout v0.14.x
$ git pull

Cherry pick the first parent branch of the this PR on top of the older branch:

$ git cherry-pick -m1 f0d48db4c246989182aa01c837d04903bc2330ae

You will likely have some merge/cherry-pick conflict here, fix them and commit:

$ git commit -am 'Backport PR #3486: Use fused types in denoise, warp'

Push to a named branch :

git push YOURFORK v0.14.x:auto-backport-of-pr-3486-on-v0.14.x

Create a PR against branch v0.14.x, I would have named this PR:

"Backport PR #3486 on branch v0.14.x"

And apply the correct labels and milestones.

Congratulation you did some good work ! Hopefully your backport PR will be tested by the continuous integration and merged soon!

If these instruction are inaccurate, feel free to suggest an improvement.

* MNT: Add a fused numeric type to make fused_types more constent. * MNT: make the interpolation use fused types for any->any type interpolation of images. * MNT: Move even more numpy function calls from Cython to Python. * MNT: More explicit type specifying * BUG: Use fused floats in denoise Pull all denoise array allocation out into python Function used in python should be def not cdef Have the correct number of arguments * MNT: Add additional type check for denoise, more docs * PEP8 fixes * Remove redundant array u and view cu * Mutate output array in Cython and trim in Python * Replace height and width with existing rows and cols * Change iteration order to match array order * Move range_lut and color_lut properly from Cy to Py * Ravel range LUT which is expected to be 1D * Update comment in denoise bilateral tests * Fix relative import in interpolation.pyx * BENCH: Benchmark for warping with many types * Rename warp benchmark file * Update pointer syntax Co-Authored-By: Mark Harfouche <mark.harfouche@gmail.com>

jni · 2019-03-06T07:15:46Z

I've submitted a manual backport.

hmaarrfk · 2019-05-01T02:56:44Z

@meeseeksdev backport to v0.14.x

lumberbot-app · 2019-05-01T02:56:48Z

Something went wrong ... Please have a look at my logs.

* MNT: Add a fused numeric type to make fused_types more constent. * MNT: make the interpolation use fused types for any->any type interpolation of images. * MNT: Move even more numpy function calls from Cython to Python. * MNT: More explicit type specifying * BUG: Use fused floats in denoise Pull all denoise array allocation out into python Function used in python should be def not cdef Have the correct number of arguments * MNT: Add additional type check for denoise, more docs * PEP8 fixes * Remove redundant array u and view cu * Mutate output array in Cython and trim in Python * Replace height and width with existing rows and cols * Change iteration order to match array order * Move range_lut and color_lut properly from Cy to Py * Ravel range LUT which is expected to be 1D * Update comment in denoise bilateral tests * Fix relative import in interpolation.pyx * BENCH: Benchmark for warping with many types * Rename warp benchmark file * Update pointer syntax Co-Authored-By: Mark Harfouche <mark.harfouche@gmail.com>

* Backport: use fused types in denoise, warp (#3486) * MNT: Add a fused numeric type to make fused_types more constent. * MNT: make the interpolation use fused types for any->any type interpolation of images. * MNT: Move even more numpy function calls from Cython to Python. * MNT: More explicit type specifying * BUG: Use fused floats in denoise Pull all denoise array allocation out into python Function used in python should be def not cdef Have the correct number of arguments * MNT: Add additional type check for denoise, more docs * PEP8 fixes * Remove redundant array u and view cu * Mutate output array in Cython and trim in Python * Replace height and width with existing rows and cols * Change iteration order to match array order * Move range_lut and color_lut properly from Cy to Py * Ravel range LUT which is expected to be 1D * Update comment in denoise bilateral tests * Fix relative import in interpolation.pyx * BENCH: Benchmark for warping with many types * Rename warp benchmark file * Update pointer syntax Co-Authored-By: Mark Harfouche <mark.harfouche@gmail.com> * forward port check_sdist

AetherUnbound approved these changes Oct 19, 2018

View reviewed changes

hmaarrfk requested changes Oct 20, 2018

View reviewed changes

skimage/restoration/_denoise_cy.pyx Show resolved Hide resolved

hmaarrfk requested changes Oct 20, 2018

View reviewed changes

skimage/restoration/_denoise.py Outdated Show resolved Hide resolved

skimage/restoration/_denoise_cy.pyx Outdated Show resolved Hide resolved

hmaarrfk reviewed Oct 20, 2018

View reviewed changes

hmaarrfk reviewed Oct 22, 2018

View reviewed changes

skimage/_shared/interpolation.pxd Show resolved Hide resolved

jni force-pushed the bugfix/float-cast branch from c73f259 to 5f78d05 Compare November 29, 2018 14:00

AetherUnbound approved these changes Dec 3, 2018

View reviewed changes

hmaarrfk reviewed Mar 5, 2019

View reviewed changes

skimage/_shared/interpolation.pxd Outdated Show resolved Hide resolved

Update pointer syntax

346f465

... Hold on to your butts! =P Co-Authored-By: jni <juan.nunez-iglesias@monash.edu>

stefanv merged commit f0d48db into scikit-image:master Mar 6, 2019

stefanv added a commit that referenced this pull request Mar 6, 2019

Revert "Use fused types in denoise, warp (#3486)"

f2bf40a

This reverts commit f0d48db.

lumberbot-app bot added the Still Needs Manual Backport MrMeeseeks-managed label label Mar 6, 2019

hmaarrfk mentioned this pull request Mar 9, 2019

Use float32 when warping when provided by the user. #3798

Closed

9 tasks

This was referenced Mar 14, 2019

Added fused type warping #1977

Closed

Use fused types for warping #1287

Closed

jni mentioned this pull request Apr 11, 2019

About geometry transform (dtype and performance) #3833

Open

hmaarrfk mentioned this pull request May 6, 2019

Backport: use fused types in denoise, warp (#3486) #3787

Closed

9 tasks

jni mentioned this pull request Mar 25, 2020

ValueError: Big-endian buffer not supported on little-endian compiler, when using skimage.transform.resize() #4525

Closed

Uh oh!

Conversation

jni commented Oct 19, 2018

Description

Checklist

For reviewers

Uh oh!

pep8speaks commented Oct 19, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Comment last updated at 2019-03-05 23:40:45 UTC

Uh oh!

hmaarrfk commented Oct 19, 2018

Uh oh!

AetherUnbound left a comment

Choose a reason for hiding this comment

Uh oh!

hmaarrfk commented Oct 19, 2018

Uh oh!

jni commented Oct 19, 2018 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

hmaarrfk commented Oct 19, 2018

Uh oh!

jni commented Oct 20, 2018

Uh oh!

hmaarrfk commented Oct 20, 2018

Uh oh!

jni commented Oct 20, 2018

Uh oh!

hmaarrfk commented Oct 20, 2018

Uh oh!

hmaarrfk commented Oct 20, 2018

Uh oh!

hmaarrfk commented Oct 20, 2018

Uh oh!

Uh oh!

Uh oh!

Uh oh!

hmaarrfk Oct 20, 2018

Choose a reason for hiding this comment

Uh oh!

hmaarrfk commented Oct 20, 2018

Uh oh!

Uh oh!

Uh oh!

Uh oh!

hmaarrfk commented Oct 20, 2018

Uh oh!

hmaarrfk Oct 20, 2018

Choose a reason for hiding this comment

Uh oh!

hmaarrfk commented Oct 20, 2018

Uh oh!

jni commented Oct 21, 2018

Uh oh!

hmaarrfk commented Oct 21, 2018

Uh oh!

Uh oh!

jni commented Nov 29, 2018

Uh oh!

jni commented Nov 29, 2018

Uh oh!

jni commented Dec 1, 2018

Uh oh!

AetherUnbound left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

AetherUnbound commented Dec 3, 2018

Uh oh!

jni commented Dec 4, 2018

Uh oh!

jni commented Dec 4, 2018

Uh oh!

jni commented Mar 5, 2019

Uh oh!

pep8speaks commented Oct 19, 2018 •

edited

Loading

jni commented Oct 19, 2018 •

edited

Loading