ratt-ru / pfb-imaging

Preconditioned forward/backward clean algorithm
MIT License
7 stars 5 forks source link

spi fitter crashes with "ZeroDivisionError: division by zero" with two freq planes #40

Closed o-smirnov closed 3 years ago

o-smirnov commented 3 years ago
Fitting 516042 components
Traceback (most recent call last):
  File "/home/oms/.venv/pfb/bin/simple_spi_fitter.py", line 7, in <module>
    exec(compile(f.read(), __file__, 'exec'))
  File "/home/oms/projects/pfb-clean/scripts/simple_spi_fitter.py", line 366, in <module>
    main(args)
  File "/home/oms/projects/pfb-clean/scripts/simple_spi_fitter.py", line 295, in main
    np.float64(ref_freq), beam=beam_comps).compute()
  File "/home/oms/.venv/pfb/lib/python3.6/site-packages/dask/base.py", line 279, in compute
    (result,) = compute(self, traverse=False, **kwargs)
  File "/home/oms/.venv/pfb/lib/python3.6/site-packages/dask/base.py", line 561, in compute
    results = schedule(dsk, keys, **kwargs)
  File "/home/oms/.venv/pfb/lib/python3.6/site-packages/dask/threaded.py", line 84, in get
    **kwargs
  File "/home/oms/.venv/pfb/lib/python3.6/site-packages/dask/local.py", line 487, in get_async
    raise_exception(exc, tb)
  File "/home/oms/.venv/pfb/lib/python3.6/site-packages/dask/local.py", line 317, in reraise
    raise exc
  File "/home/oms/.venv/pfb/lib/python3.6/site-packages/dask/local.py", line 222, in execute_task
    result = _execute_task(task, data)
  File "/home/oms/.venv/pfb/lib/python3.6/site-packages/dask/core.py", line 121, in _execute_task
    return func(*(_execute_task(a, cache) for a in args))
  File "/home/oms/.venv/pfb/lib/python3.6/site-packages/dask/optimization.py", line 963, in __call__
    return core.get(self.dsk, self.outkey, dict(zip(self.inkeys, args)))
  File "/home/oms/.venv/pfb/lib/python3.6/site-packages/dask/core.py", line 151, in get
    result = _execute_task(task, cache)
  File "/home/oms/.venv/pfb/lib/python3.6/site-packages/dask/core.py", line 121, in _execute_task
    return func(*(_execute_task(a, cache) for a in args))
  File "/home/oms/projects/codex-africanus/africanus/model/spi/dask.py", line 28, in _fit_spi_components_wrapper
    maxiter=maxiter)
  File "/home/oms/projects/codex-africanus/africanus/model/spi/component_spi.py", line 86, in fit_spi_components
    tol, maxiter, mindet)
ZeroDivisionError: division by zero

Files are in /net/young/home/oms/projects/Shapley/selfcal-shap1/spi-2band, running as

simple_spi_fitter.py -model a3562-model.fits -residual a3562-residual.fits -o a3562-spi2band-thr5 -bm a3562-beam.fits -th 5 -pf 4 
landmanbester commented 3 years ago

Ah, I see what's happening. The effective number of degrees of freedom is zero and I divide by it to set the chisq per deg of freedom to one (so that error bars have the right scale). I've fixed it here

https://github.com/landmanbester/codex-africanus/tree/nonzero_dof

Should be in master in about an hr or so

landmanbester commented 3 years ago

The fix is in africanus master branch now if you want to give it a try @o-smirnov

o-smirnov commented 3 years ago

Yep, works now, thanks!

o-smirnov commented 3 years ago

The error map is very funky though. Is it supposed to be reasonable?

landmanbester commented 3 years ago

Well its a hack for sure. The error map conflates goodness of fit with uncertainties coming from the data. Also, in this case, the effective number of degrees of freedom is zero (2 data points and 2 parameters) so I would have thought that it fits the data perfectly. Are the uncertainties higher or lower than you expect?

o-smirnov commented 3 years ago

Are the uncertainties higher or lower than you expect?

Both, lol... A few isolated 1e+16 pixels in a sea of ~1e-16.

landmanbester commented 3 years ago

That's a bit worrying. We could always fit a parametrised model to the dirty image using pfb iterations but that is a problem for future me. Instead. I'll issue the following disclaimer: take those error bars with a heap of salt!

o-smirnov commented 3 years ago

With a sea of sniffing salts, even.

My faith in the veracity of spectral index measurements of low surface brightness objects is currently at the same level as my faith in the incorruptibility of the ANC.... too bad this is our science case. :rofl:

landmanbester commented 3 years ago

Hahaha. They are slightly more reliable when fitting with nband > 2