landmanbester / spimple

Spectral index fitting made simple
MIT License
3 stars 3 forks source link

Handle nan's more elegantly #11

Open francescaLoi opened 2 years ago

francescaLoi commented 2 years ago

I got this error:

# 2022-11-02 14:38:29: Running spimple-imconv --image /stimela_mount/output/continuum/polarization/images_true_all_mfs05-0006-Q-image.fits --output-filename /stimela_mount/output/continuum/polarization/images_true_all_mfs05-0006-Q-image_20-20-0 --psf-pars 20.0 20.0 0.0 --nthreads 50 --band l
# INFO      14:38:30 - IMCONV             | Input Options:
# INFO      14:38:30 - IMCONV             |                          image = /stimela_mount/output/continuum/polarization/images_true_all_mfs05-0006-Q-image.fits
# INFO      14:38:30 - IMCONV             |                output_filename = /stimela_mount/output/continuum/polarization/images_true_all_mfs05-0006-Q-image_20-20-0
# INFO      14:38:30 - IMCONV             |                       psf_pars = [20.0, 20.0, 0.0]
# INFO      14:38:30 - IMCONV             |                       nthreads = 50
# INFO      14:38:30 - IMCONV             |                       circ_psf = False
# INFO      14:38:30 - IMCONV             |                     beam_model = None
# INFO      14:38:30 - IMCONV             |                           band = l
# INFO      14:38:30 - IMCONV             |                         pb_min = 0.05
# INFO      14:38:30 - IMCONV             |                   padding_frac = 0.5
# INFO      14:38:30 - IMCONV             |                      out_dtype = f4
# INFO      14:38:30 - IMCONV             | Using emaj = 2.00e+01, emin = 2.00e+01, PA = 0.00e+00
# INFO      14:38:31 - IMCONV             | Convolving image
# /usr/local/lib/python3.8/dist-packages/spimple/utils.py:87: RuntimeWarning: divide by zero encountered in double_scalars
#   A = np.array([[1. / Smin ** 2, 0],
# /usr/local/lib/python3.8/dist-packages/spimple/utils.py:88: RuntimeWarning: divide by zero encountered in double_scalars
#   [0, 1. / Smaj ** 2]])
# Traceback (most recent call last):
#   File "/usr/local/bin/spimple-imconv", line 8, in <module>
#     sys.exit(image_convolver())
#   File "/usr/local/lib/python3.8/dist-packages/spimple/apps/image_convolver.py", line 130, in image_convolver
#     image, gausskern = convolve2gaussres(imagei, xx, yy, gaussparf,
#   File "/usr/local/lib/python3.8/dist-packages/spimple/utils.py", line 163, in convolve2gaussres
#     thiskern = Gaussian2D(xx, yy, gausspari[i], normalise=norm_kernel)
#   File "/usr/local/lib/python3.8/dist-packages/spimple/utils.py", line 100, in Gaussian2D
#     idx = ind[:, 0]
# IndexError: too many indices for array: array is 1-dimensional, but 2 were  indexed
# 2022-11-02 14:38:33: spimple-imconv exited with code 1

I think is because the input is a nan image... I'll include a check on CARACal for now but it would be better to fix it here

landmanbester commented 2 years ago

What I can do is raise an error if you pass an image containing only nan's. I don't think it is sensible to expect spimple to run through in this case, that should probably be handled by the pipeline (you probably want to abandon subsequent steps because they will just be a waste of time). Easy enough to add a check but please first try this with the current release. The error will be different, I don't see the above lines in the current code base

francescaLoi commented 2 years ago

Sorry @landmanbester , I was still using spimple through stimela. When I run the latest version this is the output:

/home/floi/.local/lib/python3.8/site-packages/numba/core/errors.py:175: UserWarning: Insufficiently recent colorama version found. Numba requires colorama >= 0.3.9
  warnings.warn(msg)
INFO      09:27:09 - IMCONV             | Input Options:
INFO      09:27:09 - IMCONV             |                          image = mfs05/output/continuum/polarization/images_true_all_mfs05-0006-Q-image.fits
INFO      09:27:09 - IMCONV             |                output_filename = /56TB_2/fornax/MFS-pol-imaging/prova
INFO      09:27:09 - IMCONV             |                       psf_pars = [0.00555, 0.00555, 0.0]
INFO      09:27:09 - IMCONV             |                       nthreads = 56
INFO      09:27:09 - IMCONV             |                       circ_psf = False
INFO      09:27:09 - IMCONV             |                         dilate = 1.05
INFO      09:27:09 - IMCONV             |                     beam_model = None
INFO      09:27:09 - IMCONV             |                           band = l
INFO      09:27:09 - IMCONV             |                         pb_min = 0.05
INFO      09:27:09 - IMCONV             |                   padding_frac = 0.5
INFO      09:27:09 - IMCONV             |                      out_dtype = f4
INFO      09:27:09 - IMCONV             | Using emaj = 5.55e-03, emin = 5.55e-03, PA = 0.00e+00
INFO      09:27:09 - IMCONV             | Convolving image
/home/floi/.local/lib/python3.8/site-packages/spimple/utils.py:92: RuntimeWarning: divide by zero encountered in double_scalars
  A = np.array([[1. / Smin ** 2, 0],
/home/floi/.local/lib/python3.8/site-packages/spimple/utils.py:93: RuntimeWarning: divide by zero encountered in double_scalars
  [0, 1. / Smaj ** 2]])
/home/floi/.local/lib/python3.8/site-packages/spimple/utils.py:168: RuntimeWarning: invalid value encountered in true_divide
  convkernhat = np.where(np.abs(thiskernhat)>1e-10, gausskernhat/thiskernhat, 0.0)
INFO      09:27:12 - IMCONV             | Wrote clean psf to /56TB_2/fornax/MFS-pol-imaging/prova.clean_psf.fits
INFO      09:27:14 - IMCONV             | Wrote convolved image to /56TB_2/fornax/MFS-pol-imaging/prova.convolved.fits
INFO      09:27:14 - IMCONV             | All done here

So all good, I just need an updated version of spimple on stimela. Can you please help me with this or shall we ask @Mulan-94 ? Thanks!

landmanbester commented 2 years ago

No problem @francescaLoi, I thought that may be the case. I have no experience with the old stimela framework so you would have to check with @Mulan-94. This would be much simpler if you could be persuaded to switch to stimela2. Maybe another busy week is in order

landmanbester commented 1 year ago

@Athanaseus has convinced me of the need to do this more intelligently. I hadn't quite appreciated that a typical mosaiced image inevitably contains nan's. The dropzweights branch will already look for and prune bands that are completely nan (should I also do this for bands that are all zeros?). I can also add a feature to deselect bands manually. Partial nan components (i.e. if some of the bands for a component are nan but the band still contains valid data elsewhere) will require a bit more thinking about but shouldn't be a showstopper

Mulan-94 commented 1 year ago

Does masking work in this case btw?

landmanbester commented 1 year ago

There is currently no option to pass in a mask but this should be trivial to add if there is a need for it. The mask for the fit is currently created by deselecting all components which either i) have a minimum flux along the frequency axis less than threshold * rms (or model/max()/maxDR if no residual is passed in) or ii) contain a nan in any frequency band. Is there a need for an external mask as well?