From 1ee5959301755b05f3e49257546d531dcb634954 Mon Sep 17 00:00:00 2001 From: landmanbester Date: Tue, 14 May 2024 14:47:32 +0200 Subject: [PATCH 1/5] add S band to katbeam options --- spimple/apps/image_convolver.py | 9 +++++++-- spimple/apps/spi_fitter.py | 9 ++++++--- 2 files changed, 13 insertions(+), 5 deletions(-) diff --git a/spimple/apps/image_convolver.py b/spimple/apps/image_convolver.py index b229ae8..bb27ef9 100644 --- a/spimple/apps/image_convolver.py +++ b/spimple/apps/image_convolver.py @@ -36,7 +36,7 @@ def image_convolver(): "Use power_beam_maker to make power beam " "corresponding to image. ") parser.add_argument('-band', "--band", type=str, default='l', - help="Band to use with JimBeam. L or UHF") + help="Band to use with JimBeam. L, UHF or S") parser.add_argument('-pb-min', '--pb-min', type=float, default=0.05, help="Set image to zero where pb falls below this value") parser.add_argument('-pf', '--padding-frac', type=float, default=0.5, @@ -181,8 +181,13 @@ def image_convolver(): from katbeam import JimBeam if opts.band.lower() == 'l': beam = JimBeam('MKAT-AA-L-JIM-2020') - else: + elif opts.band.lower() == 'uhf': beam = JimBeam('MKAT-AA-UHF-JIM-2020') + elif opts.band.lower() == 's': + beam = JimBeam('MKAT-AA-S-JIM-2020') + else: + raise ValueError(f"Unknown beam model for katbeam in band {opts.band}") + beam_image = np.zeros(image.shape, dtype=opts.out_dtype) for v in range(freqs.size): beam_image[v] = beam.I(xx, yy, freqs[v]/1e6) # freqs in MHz diff --git a/spimple/apps/spi_fitter.py b/spimple/apps/spi_fitter.py index b9a054d..58492e8 100755 --- a/spimple/apps/spi_fitter.py +++ b/spimple/apps/spi_fitter.py @@ -94,7 +94,7 @@ def spi_fitter(): parser.add_argument('-ct', '--corr-type', type=str, default='linear', help="Correlation typ i.e. linear or circular. ") parser.add_argument('-band', "--band", type=str, default='l', - help="Band to use with JimBeam. L or UHF") + help="Band to use with JimBeam. L, UHF or S") parser.add_argument('-db', '--deselect-bands', default=None, nargs='+', type=int, help="Indices of subbands to exclude from the fitting \n" "By default, all the sub-bands are used for the residual image. \n" @@ -252,9 +252,12 @@ def spi_fitter(): from katbeam import JimBeam if opts.band.lower() == 'l': beam = JimBeam('MKAT-AA-L-JIM-2020') - else: + elif opts.band.lower() == 'uhf': beam = JimBeam('MKAT-AA-UHF-JIM-2020') - beam_image = np.zeros(model.shape, dtype=opts.out_dtype) + elif opts.band.lower() == 's': + beam = JimBeam('MKAT-AA-S-JIM-2020') + else: + raise ValueError(f"Unknown beam model for katbeam in band {opts.band}") for v in range(freqs.size): beam_image[v] = beam.I(xx, yy, freqs[v]/1e6) # freqs in MHz From f387af6bb20dde654f6198c01c87e73135b5d75f Mon Sep 17 00:00:00 2001 From: Athanaseus Javas Ramaila Date: Thu, 16 May 2024 07:49:37 +0200 Subject: [PATCH 2/5] debug --- spimple/apps/spi_fitter.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/spimple/apps/spi_fitter.py b/spimple/apps/spi_fitter.py index 58492e8..a96d2e1 100755 --- a/spimple/apps/spi_fitter.py +++ b/spimple/apps/spi_fitter.py @@ -250,6 +250,7 @@ def spi_fitter(): beam_image = load_fits(opts.beam_model, dtype=opts.out_dtype).squeeze() elif opts.beam_model == "JimBeam": from katbeam import JimBeam + beam_image = [] if opts.band.lower() == 'l': beam = JimBeam('MKAT-AA-L-JIM-2020') elif opts.band.lower() == 'uhf': @@ -259,7 +260,8 @@ def spi_fitter(): else: raise ValueError(f"Unknown beam model for katbeam in band {opts.band}") for v in range(freqs.size): - beam_image[v] = beam.I(xx, yy, freqs[v]/1e6) # freqs in MHz + beam_image.append(beam.I(xx, yy, freqs[v]/1e6)) # freqs in MHz + beam_image = np.stack(beam_image) else: beam_image = interpolate_beam(xx, yy, freqs, opts) From 99d185aaca6d537259bce4e75f568a516c8f670d Mon Sep 17 00:00:00 2001 From: landmanbester Date: Mon, 29 Jul 2024 15:25:24 +0200 Subject: [PATCH 3/5] add residual after fit as output --- spimple/apps/spi_fitter.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/spimple/apps/spi_fitter.py b/spimple/apps/spi_fitter.py index 58492e8..0a760f5 100755 --- a/spimple/apps/spi_fitter.py +++ b/spimple/apps/spi_fitter.py @@ -43,7 +43,7 @@ def spi_fitter(): "Default of zero means use all threads") parser.add_argument('-pb-min', '--pb-min', type=float, default=0.15, help="Set image to zero where pb falls below this value") - parser.add_argument('-products', '--products', default='aeikIcmrb', type=str, + parser.add_argument('-products', '--products', default='aeikIcmrbd', type=str, help="Outputs to write. Letter correspond to: \n" "a - alpha map \n" "e - alpha error map \n" @@ -54,6 +54,7 @@ def spi_fitter(): "m - convolved model \n" "r - convolved residual \n" "b - average power beam \n" + "d - difference between data and fitted model \n" "Default is to write all of them") parser.add_argument('-pf', "--padding-frac", default=0.5, type=float, help="Padding factor for FFT's.") @@ -452,15 +453,26 @@ def spi_fitter(): alpha_err_map[maskindices[:, 0], maskindices[:, 1]] = alpha_err i0map[maskindices[:, 0], maskindices[:, 1]] = Iref i0_err_map[maskindices[:, 0], maskindices[:, 1]] = i0_err + Irec_cube = i0map[None, :, :] * (freqs[:, None, None]/ref_freq)**alphamap[None, :, :] + fit_diff = np.zeros_like(model) + fit_diff[...] = np.nan + ix = maskindices[:, 0] + iy = maskindices[:, 1] + fit_diff[:, ix, iy] = model[:, ix, iy]/beam_image[:, ix, iy] + fit_diff[:, ix, iy] -= Irec_cube[:, ix, iy] if 'I' in opts.products: # get the reconstructed cube - Irec_cube = i0map[None, :, :] * \ - (freqs[:, None, None]/ref_freq)**alphamap[None, :, :] name = outfile + '.Irec_cube.fits' save_fits(name, np.expand_dims(Irec_cube, axis=4 - stokes_axis), mhdr, dtype=opts.out_dtype) print(f"Wrote reconstructed cube to {name}", file=log) + if 'd' in opts.products: + # get the reconstructed cube + name = outfile + '.fit_diff.fits' + save_fits(name, np.expand_dims(fit_diff, axis=4 - stokes_axis), mhdr, dtype=opts.out_dtype) + print(f"Wrote reconstructed cube to {name}", file=log) + # save alpha map if 'a' in opts.products: name = outfile + '.alpha.fits' From cc7433bbde46312733a27d2bad41c61db3b4be2e Mon Sep 17 00:00:00 2001 From: landmanbester Date: Wed, 31 Jul 2024 18:30:08 +0200 Subject: [PATCH 4/5] fix uninitialised beam image when using katbeam --- spimple/apps/spi_fitter.py | 1 + 1 file changed, 1 insertion(+) diff --git a/spimple/apps/spi_fitter.py b/spimple/apps/spi_fitter.py index 0a760f5..1e42e22 100755 --- a/spimple/apps/spi_fitter.py +++ b/spimple/apps/spi_fitter.py @@ -259,6 +259,7 @@ def spi_fitter(): beam = JimBeam('MKAT-AA-S-JIM-2020') else: raise ValueError(f"Unknown beam model for katbeam in band {opts.band}") + beam_image = np.zeros_like(model) for v in range(freqs.size): beam_image[v] = beam.I(xx, yy, freqs[v]/1e6) # freqs in MHz From ba3692c640be261ec9598b6322e12f4821e4ab7b Mon Sep 17 00:00:00 2001 From: landmanbester Date: Thu, 24 Oct 2024 09:02:01 +0200 Subject: [PATCH 5/5] np.float64 -> float in spimple-spifit parser --- spimple/apps/spi_fitter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/spimple/apps/spi_fitter.py b/spimple/apps/spi_fitter.py index b9a054d..5997426 100755 --- a/spimple/apps/spi_fitter.py +++ b/spimple/apps/spi_fitter.py @@ -67,7 +67,7 @@ def spi_fitter(): help="Per-channel freqs to use during fit to frequency axis\n"+ "These are automatically generated from the FITS header.\n"+ "Only specify if you have a non-standard FITS header.") - parser.add_argument('-rf', '--ref-freq', default=None, type=np.float64, + parser.add_argument('-rf', '--ref-freq', default=None, type=float, help='Reference frequency where the I0 map is sought. \n' "Will overwrite in fits headers of output.") parser.add_argument('-otype', '--out_dtype', default='f4', type=str,