Faster LACosmic

I have implemented a fast version of the LACosmic algorithm in C with a python wrapper. The code is posted at https://github.com/cmccully/lacosmicx. I was wondering if it would be useful to incorporate into ccdproc/astropy. My implementation is multithreaded using OpenMP and also makes the approximation of using the separable 2d median filter. I ran some performance tests today (based solely on speed). My test image set included 16 4096x2048 images. I found that the current implementation of lacosmic in ccdproc took 261 seconds to run one iteration per image. My implementation with multithreading enabled took 52 seconds for 4 iterations per image, an improvement of a factor of ~20. Running it as a single thread took 180 seconds for 4 iterations per image, still a factor ~6 over the current implementation.

Hi Curtis-- Thanks so much for running the tests. A few quick questions:

How do the results compare? Are the same cosmic rays detected when using the same parameters?
Can you run yours without the multithreading and make a comparison?

It would also be interesting to compare the multithreading by using the built in python multithreading. If it can be all done in python, it will be much easier to maintain, but if the c-code is truly faster, than it is definitely something we will want to include.

@mwcraig

Just for reference here are astropy's information about c extenstions: https://astropy.readthedocs.org/en/stable/development/ccython.html

And @mwcraig has been handling much more of the package layout. There will be additional c-code that we want at some point, so it wold be good to set this up now for reference for future contributions

I can handle adding in C code, I think -- but will have very limited availability from Jun 27-Jul 5

How do the results compare? Are the same cosmic rays detected when using the same parameters?
I haven't done this yet, but I will. I have checked my results pretty extensively against cosmics.py: http://obswww.unige.ch/~tewes/cosmics_dot_py/ and got decent results. I will run some comparisons tomorrow.
Can you run yours without the multithreading and make a comparison?
I have done this. This is the second figure of 180 seconds for 4 iterations per image, still a factor of ~6 faster than the current implementation in ccdproc.

It may be possible to implement some of my optimizations (like multithreading) in pure python. I agree that is certainly easier to maintain. We would likely still take a 50% performance loss over pure C (I have run some basic tests on this), but that isn't really much.

Sorry I missed the line about the single thread. Then it definitely looks like the c-implemenation is something that we want to add. Do you want to take a shot at making a pull request and adding in your code? Matt and I can help you through it.

There are a couple of ways I could implement the C code into the python. I could treat it is a function library and write the driving logic in python. That might be easier to maintain. Currently the c code has all of the logic and python just calls the main routine (basically).

Another possibility is that I could implement some of my optimizations into the python version. Basically you do two median filter calls but only use a 1d kernel. Is it possible to multithread scipy functions without rewriting the base c/fortran code?

I don't know if multithreading is possible with scipy...I'd lean towards more python-y, and ease of maintenance matters a lot (my availability is going to go down quite a bit in a month or two...)

@cmccully

Oops, meant to add -- thanks very much @cmccully for taking this on!

Check out saltcrclean in pysalt. there is an example of python multi-threading in there.

Ok. I have looked at the multithreading in Python (which I have used for a few basic things before). I think adding the parallelization in Python will make the code more complicated and difficult to maintain. I think I have a good idea of how I want to implement things so I will take a stab at it. Thanks for all of the help!

Yeah, the little bit I've done with threading in python has not been fun. Seems too low level right now to be convenient.

assuming ndimage releases the GIL one could just run the ndimage.median_filter in parallel using python threads and slightly overlapping image slices. Should be quite simple to implement.
something along the line of this:

kh = kernel_size // 2
res[:kh,:] = handle_top(img[:kernel_size,:])[:kg,:]
for i in range(-kh, ny - kh, block_size):
    # filter overlapping blocks to get correct borders on the block_size
    threadpool.apply_async(handle(img[i - kh:i + block_size + kh,:]))
res[-kh:,:] = handle_bottom(img[:kernel_size,:])[-kh:,:]
for t in threads:
    res[a:b,:] = t.get()[kh:-kh,:] # slice of the borders where the blocking is invalid

Also interesting would be using a median-heap based filter instead of a quickselect based one.
A lacosmic implementation based on ESO's CPL library which has such a filter runs through four iterations on a 4kx2k image in 50 seconds on an old core2duo with only one thread.

Another option could be shrinking the image to lots of smaller regions around the detected cosmics at each iteration, as if I understand correctly the iterations following the first just eat themselves into already existing detections and don't add new ones. Though the implementation for that is possibly more complex.

Hi Julian. Sorry for the delay in responding to you. I am finally settled in my new postdoc position at LCOGT. I have looked at some of the heap sorting algorithms, but I think I read they are still slower than quickselect for small kernels (which lacosmic uses). Do you know if that is the case? I have C code to do most of this so it is just a matter of translating and cleaning it up a bit. I will try to get something uploaded soon so you can see what you think. I am happy if you want to compare the speed using a heap filter instead of quickselect.

Navigation Menu

Repository navigation

Faster LACosmic #125

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Faster LACosmic #125

Description

Activity

crawfordsm commented on Jun 24, 2014

crawfordsm commented on Jun 24, 2014

mwcraig commented on Jun 24, 2014

cmccully commented on Jun 24, 2014

crawfordsm commented on Jun 24, 2014

cmccully commented on Jun 25, 2014

mwcraig commented on Jun 26, 2014

mwcraig commented on Jun 26, 2014

crawfordsm commented on Jun 26, 2014

cmccully commented on Jun 26, 2014

mwcraig commented on Jul 2, 2014

juliantaylor commented on Aug 19, 2014

cmccully commented on Sep 3, 2014

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions