Lacosmicx

[cmccully]

I have finished the optimized LA Cosmic implementation. It uses Cython for the main part of the algorithm and uses C for the individual routines. Everything is multithreaded using OpenMP (the Cython standard). For a 4K x 2K image, this implementation is ~35 times faster than cosmics.py written by Malte Tewes (which is already significantly faster than the original IRAF version).

Given the complexity of this code, should this be included in imageutils or should it be its own package?

Any other comments are welcome.

[cdeil]

@larrybradley and / or @crawfordsm Do you have time to review this?

[larrybradley]

I'll review this next week.

[astrofrog]

@cdeil @larrybradley - @cmccully raised a valid question as to whether this should be included in imageutils or should be its own package? Another option is to include it as part of ccdproc. I don't have a strong preference but just to be clear, the plan is to merge imageutils into astropy core within a couple of months, so this needs to be stable enough for merging API-wise.

[eteq]

I think this has to be a dependency for photutils and ccdproc, right? If so, then it can't be part of either separately. And I think a separate package just for this doesn't make sense because it will be a pain to maintain in lockstep with photutils and ccdproc. So that makes me think it should go into imageutils. But I agree with @astrofrog it has to be at least reasonably API-stable if we do that.

[eteq]

Also, @cmccully, am I reading right that it does not parallelize on clang? I ask because I know a lot of the user base are on macs...

[crawfordsm]

I'd recently reviewed an early version of this and a lot of my comments had already been integrated, but I suggested doing a pull request to get a few more eyes and thoughts about this as I'm not the most experienced person especially with astrop setup.py issues and also to take a look at what other build issues there might be.

There is an issue with clang as it does not include OpenMP and so there does need to be a switch to turn off all of the parrallelization code if it is compiling with compilers that it doesn't work with. There are some examples on google, but I wasn't sure if there was something astropy was already checking for or doing. So I think any thoughts about the best way to deal with that might be very helpful at this time.

[crawfordsm]

For the reasons @eteq states, I would support a version of this should be in imageutils, whether it is @cmccully or @larrybradley. Although the improved performance of @cmccully is a real benefit even if it does take a bit more work to get it stable.

[cmccully]

I am happy to have it in either location, but it sounds like it may be worth including in imageutils.

As Steve mentioned, clang does not support OpenMP as of yet. I have setup the code to build either way, (it doesn't fail with clang, but it only runs on one core). In principle, one could install http://clang-omp.github.io/ and then it would be parallel. Cython's parallel capabilities also rely on OpenMP so I am not sure that this should be an issue.

I see that the Travis Cl build failed. When I merged things, it seemed to work. Can someone give me a little direction on what that means? Thanks.

[cmccully]

I looked at the Travis Cl log. I ran into this issue on my own machine at one point. There is an undefined symbol: __pow_finite. This happened to me because the Anaconda python distribution ships with a bad libm. I don't know what the best way to protect against this issue.

[cmccully]

I removed the _astropy_init.py file as requested. I also updated the function that I believe is causing the Travis Cl build to fail. Most of the Travis Cl builds pass. The first five "egg_info" builds fail because they can't find Cython. from Cython.Build import cythonize ImportError: No module named Cython.Build Can someone let me know how to fix this?

[larrybradley]

@cmccully I can't get this to compile on my Mac (clang):

ld: library not found for -lgomp
clang: error: linker command failed with exit code 1 (use -v to see invocation)

This looks to be OpenMP related, but I thought the code would build without OpenMP.

[larrybradley]

Cython isn't installed for the egg_info builds. I'm not an expert on building C extensions, but the astropy documentation (http://docs.astropy.org/en/v0.4.2/development/ccython.html) says that get_extensions() should return a list of distutils.core.Extension objects. So I think (but I could be wrong!) that your setup_package.py should look something like:

import os
from distutils.core import Extension
from glob import glob
from astropy_helpers import setup_helpers

def get_extensions():
    cfg = setup_helpers.DistutilsExtensionArgs()
    cfg['sources'].extend(
        os.path.relpath(fname) for fname in
        glob(os.path.join(os.path.dirname(__file__), '*.pyx')))
    cfg['sources'].append('imageutils/lacosmicx/laxutils.c')
    cfg['include_dirs'].append('numpy')
    cfg['extra_compile_args'].extend(['-g', '-O3', '-fopenmp',
                                      '-funroll-loops', '-ffast-math'])
    cfg['extra_link_args'].extend(['-g', '-fopenmp'])
    return [Extension('imageutils.lacosmicx', **cfg)]

This will also pass for the egg_info test builds.

[cmccully]

Thanks @larrybradley for looking at this.

On the clang compiler: in the most recent clang update for Mac, the -fopenmp flag is no longer ignored, but actively causes the build to fail (even though OpenMP is not yet supported for clang). You can either remove the -fopenmp flag from the setup_package.py, use gcc, or install http://clang-omp.github.io/. I would very much appreciate any suggestions on how to check for this more robustly in the setup_package.py file.

As for compiling, I was not able to find an example Cython setup_package.py file so I wasn't sure how Cython was being called. The build file you posted did produce an error (which I am sure is related to my init.py file). Do you know how the init file should be updated to get rid of these issues: duplicate symbol _initlacosmicx in: build/temp.macosx-10.7-x86_64-2.7/imageutils/lacosmicx/lacosmicx.o build/temp.macosx-10.7-x86_64-2.7/imageutils/lacosmicx/laxwrappers.o duplicate symbol _pyx_module_is_main_imageutilslacosmicx in: build/temp.macosx-10.7-x86_64-2.7/imageutils/lacosmicx/lacosmicx.o build/temp.macosx-10.7-x86_64-2.7/imageutils/lacosmicx/laxwrappers.o

Thanks for all of your help.

[crawfordsm]

I did some checking around, but was hoping someone would have a better idea, but it looks like we should add something to detect clang to the setup.py and if it is, to then change so that openmp isn't used.

Here's one example: https://github.com/healpy/healpy/commit/ac5e2731e2a119bc291d370b19ea9a30496b854d

And here's another example that actually checks if openmp works by compiling a small script: https://github.com/gravitationalarry/bayestar-localization/blob/master/misc/setuptools_openmp.py

[larrybradley]

@cmccully I see those errors about duplicate symbols too. I suspect my suggested get_extensions() is in the right direction, but not completely correct. Perhaps @embray, @astrofrog, or @eteq can help with your setup_package.py and also suggest the best way to check for OpenMP so that those flags are turned off with clang?

[cmccully]

@larrybradley, @crawfordsm, @cdeil, @embray, @astrofrog, @eteq Any other suggestions for detecting clang to turn off OpenMP flags and how the setup.py package should be structured? Thanks.

[astrofrog]

@cmccully - I've put some fixes in https://github.com/cmccully/imageutils/pull/1 - this will allow travis to pass. We still need to figure out the openmp option but I'm working on it.

[astrofrog]

@cmccully - can you merge in https://github.com/cmccully/imageutils/pull/2 too? Should get Travis to pass.

[cdeil]

I think trying to compile a small test program with OpenMP (like here as suggested by @crawfordsm) with an option to turn it off (in case there's a problem) is much better than hard-coded checks which compiler is used (e.g. currently clang doesn't support OpenMP, but in a year or two it will).

We should probably file a new astropy-helpers issue for this build system question and keep the discussion on the code / method here.

I'm a bit worried about adding this to the Astropy core because of this build issue ... the fact that a add_openmp_flags_if_available function isn't available in a core package like setuptools or numpy and isn't widely used by compute-intensive package like scipy or scikit-learn could mean that it's not easy to maintain. Maybe worth asking on numpy-dev?

[eteq]

@cdeil is probably right that we should be cautious about OpenMP. It should be supported where possible given that it's the native Cython parallelization approach, but I share the concern that it might be fragile.

If we have to leave OpenMP off (or, better yet, as an "opt-in" option using the astropy configutation system), that's probably fine as a starting point, and subsquent PRs could try to integrate something into setup-helpers that that will actually check for OpenMP.

[astrofrog]

I wonder whether we could make this opt-in via an environment variable for now, which will be the easiest way for the user to change a setting in setup_package.py without having to change anything in setup_helpers for now? (so as a temporary solution). Something like ASTROPY_USE_OPENMP. If the user wants to compile with openmp they would do:

ASTROPY_USE_OPENMP=1 python setup.py install

or something like this?

[cdeil]

I also think an environment variable or command line option for setup.py so that OPENMP is off by default is the best solution for now. This way we can release 1.0 without worrying about build issues. In the future the behaviour could change so that this gets turned on automatically if / when someone implements checks in astropy-helpers to detect if OpenMP works with the user's compiler.

[crawfordsm]

So, @cmccully pointed me to here: https://github.com/pandegroup/IRMSD/blob/master/python/setup.py

which then led to a little exploring with distutils and there is a has_function built into distutils, so the following code snippet would actually detect if it is available or not:

from distutils.ccompiler import new_compiler
compiler = new_compiler()
compiler.add_library('gomp')
hasopenmp = compiler.has_function('omp_get_num_threads')

where hasopenmp would then be True if it is there. I'm sure this could all be wrapped in a try...except statement and it should be then easy to include.

[eteq]

I just checked now, and the snippet @crawfordsm gave correctly results in hasopenmp == False on OS X, but hasopenmp == True with gcc (at least on Red Hat, which I had easily accessible). So probably safe to just sub that in for USE_OPENMP = False? (with a try/except that results in False with a warning if it fails?)

[cdeil]

The check from @crawfordsm works for me as well with gcc and clang.

I see this warning printed with clang ... can this output be suppressed?

>>> hasopenmp = compiler.has_function('omp_get_num_threads')
/var/folders/sb/4qv5j4m90pz1rw7m70rj1b1r0000gn/T/omp_get_num_threads02r9yrie.c:1:1: warning: type specifier missing, defaults to 'int' [-Wimplicit-int]
main (int argc, char **argv) {
^~~~
/var/folders/sb/4qv5j4m90pz1rw7m70rj1b1r0000gn/T/omp_get_num_threads02r9yrie.c:2:5: warning: implicit declaration of function 'omp_get_num_threads' is invalid in C99
      [-Wimplicit-function-declaration]
    omp_get_num_threads();
    ^
2 warnings generated.
ld: library not found for -lgomp
clang: error: linker command failed with exit code 1 (use -v to see invocation)

It would be good to see if this check works with some more compilers. Maybe @yarikoptic has access to Debian testing machines and can help?

[crawfordsm]

@cdeil the output can be suppressed, but not easily from what I've seen so far. So it either means more code or leaving in the warnings. We could add a log message there to ignore these warnings.

[cmccully]

It is possible to suppress the warnings by doing the following:

ccompiler = new_compiler() try: ccompiler.add_library('gomp') ccompiler.set_excutable('compiler_so', ccompiler.executable('compiler_so')+' - w') USE_OPENMP=ccompiler.has_function('omp_get_num_threads') except: USE_OPENMP=False

While this suppresses the warnings, if OpenMP does not exist, the error message below will still print: ld: library not found for -lgomp clang: error: linker command failed with exit code 1 (use -v to see invocation)

[eteq]

@cmccully - I think you meant something like: ccompiler.set_executable('compiler_so', ccompiler.executables['compiler_so'][0] + ' -w'), is that right? The code you put above always gives USE_OPENMP==False, but that's because it's catching all the exceptions, including the errors from the above typos. With the change I'm showing to work, it works on Ubuntu, but with clang on osx it fails in a way that makes me think it's not even compiling (maybe -w isn't supposed to work on clang?).

[eteq]

Alright, I worked my way into the guts of distutils (I horrible thing to do, as always...), and I'm pretty sure there's no way to supress this output short of something drastic (and possibly bug-prone) like temporarily replacing stdout (and I'm not even sure that would work). So I think the best thing to do is just to do something like this:

ccompiler = new_compiler()
ccompiler.add_library('gomp')
USE_OPENMP=ccompiler.has_function('omp_get_num_threads')
if not USE_OPENMP:
    log('OpenMP was not found, continuing with non-parallelized version.  Any error messages you see immediately above this mentioning "omp_get_num_threads" can be ignored.')

Also, if it matters, the method above correctly (and silently) yields USE_OPENMP==True on Ubuntu (14.10), so that's another one that works.

[astrofrog]

@eteq - what about:

import sys
import tempfile
import subprocess

CODE = """
import sys
from distutils.ccompiler import new_compiler
ccompiler = new_compiler()
ccompiler.add_library('gomp')
sys.exit(int(ccompiler.has_function('omp_get_num_threads')))
"""

def has_openmp():
    with tempfile.NamedTemporaryFile(mode='w') as f:
        f.write(CODE)
        f.flush()
        s = subprocess.Popen([sys.executable, f.name],
                             stdout=subprocess.PIPE,
                             stderr=subprocess.PIPE)
        s.wait()
    return bool(s.returncode)

[eteq]

@astrofrog - oh, yeah, that works. Even safer might be to use stdin to pipe in the code instead of writing to a temp file. I'm not entirely sure how well that will behave when someone gives options to setup.py, but we can always try and see. The failure mode of "don't use OpenMP" should airways be safe after all.

[cmccully]

Ok. I pushed an update based on suggestions from @eteq and @astrofrog. What do you think?

[eteq]

@cmccully - see cmccully/imageutils#5 - that just adds some info that actually tells the user what's going on (and how they might debug it). Other than that, seems fine to me (assuming travis ends up happy).

[larrybradley]

@cmcully You'll need to rebase this PR on the current master before it can be merged.

Also, I did find some (small) differences compared to the version that was in photutils (which is now available at https://github.com/larrybradley/lacosmic for comparison testing), but I haven't had a chance to investigate the origin of the differences yet. I used sepmed=False, cleantype='medmask' for the comparison, which I think should emulate the "standard" lacosmic.

[cdeil]

python setup.py build_sphinx doesn't work for me locally: https://gist.github.com/cdeil/9ddcb835de9ea47bc8a1

@cmccully Can you git rebase master and then try python setup.py build_sphinx ? (if this doesn't solve the issue we might have to update astropy-helpers or investivate.)

[larrybradley]

The version of astropy-helpers in imageutils is up to date (0.4.3).

[cmccully]

Besides getting approval from Nicolas Devillard and Pieter van Dokkum, is there anything else that needs to be done on this pull request? Thanks again for everyone's help.

[ndevilla]

Hi guys:

Original author of the median code here. I hereby confirm that this code has been released in the public domain, allowing you to include it under whatever kind of license you might want to use, BSD, GNU, or whatever you please. As far as I know it has already been included in commercial software.

Fun fact: the research on median and the resulting code were conducted 15 years ago when I needed to speed up the routines I wrote for infrared image processing at the VLT. It is only fair that this code finds it way home into astropy.

Cheers Nicolas

[eteq]

@cmccully - I'll check with Pieter van Dokkum (he's just down the hall from my office...)

[astrofrog]

@ndevilla - thanks for confirming this!

@eteq - did you get a chance to speak to Pieter van Dokkum?

[astrofrog]

Just a general note - I don't want to hold this back, but at the moment this implements a fair bit of general purpose utilities (including some which overlap with other parts of imageutils and astropy). For example, there are functions for rebinning an array (which exists elsewhere in imageutils), convolution (which exists in astropy.convolution) and re-implements some of the models from astropy. I'm not suggesting that this can directly use what is in astropy because of performance considerations, but it does seem unecessary to duplicate things. So maybe in the long term we can try and improve astropy.convolution such that it exposes the functions needed here, and we can improve the rest of imageutils so that we don't need to duplicate things here.

I see two possible paths - we merge this asap and include it in astropy 1.0, but then over time try and reduce the amount of duplicated code. Or, we actually keep this open for a while while working on trying to reduce duplication.

[astrofrog]

I have specific suggestions which might help us reduce duplication in future:

• Move subsample and rebin to imageutils/sampling.pyx and give them more explicit names to indicate that they are specifically for factors of 2. It would be interesting to see fi the existing downsample and upsample can be used instead.
• Move all the math functions (specifically all the median calculation stuff) into a separate file in imageutils/lacosmicx, for example math_utils.pyx and math_utils.c.
• Similarly, move out all the convolution/filtering stuff to separate files.

Basically, if we split out things according to what kind of functionality they provide, it will be easier to try and refactor in future. At the moment, there are some very specific cosmic-ray utilities and very basic math functions in the same place.

[cdeil]

I think we should merge soon and improve refactor next year. It's unrealistic / painful to address all the points you mentioned in one Github PR. We can put a big fate note "astropy.image" is not stable at the top of the docs, like we did e.g. for modeling (and should have done for nddata / convolution).

[astrofrog]

@cdeil - ok, but I think it would still be worth splitting out the files as I suggested in my last comment, since this would be low effort for now and would make things a lot clearer to tackle in future.

[yarikoptic]

@cdeil I am sorry I have missed your original question about a Debian machine. If you don't need anything "exotic" then you could get Debian running via docker, virtualbox (e.g. we provide neurodebian appliance) or vagrant. If you need access to some specific platform (e.g. sparc), email me ;)

[cdeil]

@yarikoptic - The question is exactly if python setup.py install and python setup.test for this pull request works with exotic machines or compilers. It's probably easiest to ask for someone to do this testing after it's merge into the Astropy core repo for an Astropy 1.0 release candidate ... unless it's easy / little work for you to run this on a few test machines.

[yarikoptic]

;-) well -- seems not even an exotic one necessary (tested on my laptop and sparc box) to fail install:

In file included from /usr/lib/python2.7/dist-packages/numpy/core/include/numpy/ufuncobject.h:327:0,
                 from astropy/convolution/boundary_none.c:240:
/usr/lib/python2.7/dist-packages/numpy/core/include/numpy/__ufunc_api.h:241:1: warning: ‘_import_umath’ defined but not used [-Wunused-function]
 _import_umath(void)
 ^
error: Setup script exited with error: [Errno 2] No such file or directory: 'imageutils/cython_version.py'
PYTHONPATH=$PWD/INSTALL/lib/python2.7/site-packages/ python setup.py install   110.85s user 3.21s system 84% cpu 2:14.92 total

@cdeil email me your public ssh key and desired login name, paste here the fingerprint of your ssh key and I will give you access to an exotic box so you could try yourself

[astrofrog]

By the way, if people want to go ahead and merge this we can do that and I'm happy to do the splitting out of the files as I mentioned above.

[larrybradley]

+1 for splitting out the general purpose functions, but that can happen in a separate PR. I also wanted to investigate some small differences (when I get time), but that too can be in a separate PR.

[astrofrog]

Ok, so it seems like we're just waiting for @eteq to confirm back what Pieter van Dokkum says, and if all is well on that front, we can merge and do further work on this before starting to do a PR for astropy.