From f81a30a268af0e302da49c5eb9c27290e3c07b3a Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 2 Sep 2022 12:55:25 -0600 Subject: [PATCH 01/32] Add in --cores and --gpu as options. --- LWA_EPIC/LWA_EPIC.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 3df3b5d..6df9bb7 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -2117,6 +2117,8 @@ def gen_args(return_parser=False): group2.add_argument("--tbffile", type=str, help="TBF Data Path") group3 = parser.add_argument_group("Processing Options") + group3.add_argument("--cores", type=str, default="3,4,5,6,7", help="Comma separated list of CPU cores to bind to") + group3.add_argument("--gpu", type=int, default=0, help="GPU to bind to") group3.add_argument("--imagesize", type=int, default=64, help="1-D Image Size") group3.add_argument( "--imageres", type=float, default=1.79057, help="Image pixel size in degrees" @@ -2231,8 +2233,8 @@ def main(args, parser): log.setLevel(logging.DEBUG) # Setup the cores and GPUs to use - cores = [3, 4, 5, 6, 7] - gpus = [0, 0, 0, 0, 0] + cores = [int(v) for v in args.cores.split(',')] + gpus = [args.gpu,]*len(cores) # Setup the signal handling ops = [] From 69aa8420902ae65458f713ee03c41aaf449b2343 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 2 Sep 2022 13:09:32 -0600 Subject: [PATCH 02/32] Add some caching of locations and phase corrections to MOFFCorrelatorOp. --- LWA_EPIC/LWA_EPIC.py | 98 +++++++++++++++++++++++++------------------- 1 file changed, 56 insertions(+), 42 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 6df9bb7..047afcb 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -889,15 +889,24 @@ def main(self): itshape = (self.ntime_gulp, nchan, nstand, npol) freq = (chan0 + np.arange(nchan)) * CHAN_BW - sampling_length, locs, sll = GenerateLocations( - self.locations, - freq, - self.ntime_gulp, - nchan, - npol, - grid_size=self.grid_size, - grid_resolution=self.grid_resolution, - ) + locname = "locations_%s_%i_%i_%i_%i_%i_%i_%.6f.npy" % (self.station.name, chan0, nchan, nstand, npol, self.ant_extent, self.grid_size, self.grid_resolution) + locname = os.path.join(os.path.dirname(__file__), locname) + try: + loc_data = np.loadz(locname) + sampling_length = loc_data['sampling_length'].item() + locs = loc_data['locs'][...] + sll = loc_data['sll'].item() + except OSError: + sampling_length, locs, sll = GenerateLocations( + self.locations, + freq, + self.ntime_gulp, + nchan, + npol, + grid_size=self.grid_size, + grid_resolution=self.grid_resolution, + ) + np.savez(locname, sampling_length=sampling_length, locs=locs, sll=sll) try: copy_array(self.locs, bifrost.ndarray(locs.astype(np.int32))) except AttributeError: @@ -942,39 +951,45 @@ def main(self): # Setup the kernels to include phasing terms for zenith # Phases are Ntime x Nchan x Nstand x Npol x extent x extent - freq.shape += (1, 1) - phases = np.zeros( - (self.ntime_gulp, nchan, nstand, npol, self.ant_extent, self.ant_extent), - dtype=np.complex64 - ) - for i in range(nstand): - # X - a = self.station.antennas[2 * i + 0] - delay = a.cable.delay(freq) - a.stand.z / speed_of_light.value - phases[:, :, i, 0, :, :] = np.exp(2j * np.pi * freq * delay) - phases[:, :, i, 0, :, :] /= np.sqrt(a.cable.gain(freq)) - if npol == 2: - # Y - a = self.station.antennas[2 * i + 1] + phasename = "phases_%s_%i_%i_%i_%i_%i.npy" % (self.station.name, chan0, nchan, nstand, npol, self.ant_extent) + phasename = os.path.join(os.path.dirname(__file__), phasename) + try: + phases = np.load(phasename) + except OSError: + freq.shape += (1, 1) + phases = np.zeros( + (self.ntime_gulp, nchan, nstand, npol, self.ant_extent, self.ant_extent), + dtype=np.complex64 + ) + for i in range(nstand): + # X + a = self.station.antennas[2 * i + 0] delay = a.cable.delay(freq) - a.stand.z / speed_of_light.value - phases[:, :, i, 1, :, :] = np.exp(2j * np.pi * freq * delay) - phases[:, :, i, 1, :, :] /= np.sqrt(a.cable.gain(freq)) - # Explicit bad and suspect antenna masking - this will - # mask an entire stand if either pol is bad - if ( - self.station.antennas[2 * i + 0].combined_status < 33 - or self.station.antennas[2 * i + 1].combined_status < 33 - ): - phases[:, :, i, :, :, :] = 0.0 - # Explicit outrigger masking - we probably want to do - # away with this at some point - if ( - (self.station == lwasv and a.stand.id == 256) - or (self.station == lwa1 and a.stand.id in (35, 257, 258, 259, 260)) - ): - phases[:, :, i, :, :, :] = 0.0 - phases = phases.conj() - phases = bifrost.ndarray(phases) + phases[:, :, i, 0, :, :] = np.exp(2j * np.pi * freq * delay) + phases[:, :, i, 0, :, :] /= np.sqrt(a.cable.gain(freq)) + if npol == 2: + # Y + a = self.station.antennas[2 * i + 1] + delay = a.cable.delay(freq) - a.stand.z / speed_of_light.value + phases[:, :, i, 1, :, :] = np.exp(2j * np.pi * freq * delay) + phases[:, :, i, 1, :, :] /= np.sqrt(a.cable.gain(freq)) + # Explicit bad and suspect antenna masking - this will + # mask an entire stand if either pol is bad + if ( + self.station.antennas[2 * i + 0].combined_status < 33 + or self.station.antennas[2 * i + 1].combined_status < 33 + ): + phases[:, :, i, :, :, :] = 0.0 + # Explicit outrigger masking - we probably want to do + # away with this at some point + if ( + (self.station == lwasv and a.stand.id == 256) + or (self.station == lwa1 and a.stand.id in (35, 257, 258, 259, 260)) + ): + phases[:, :, i, :, :, :] = 0.0 + phases = phases.conj() + phases = bifrost.ndarray(phases) + np.save(phasename, phases) try: copy_array(gphases, phases) except NameError: @@ -1434,7 +1449,6 @@ def main(self): # Setup the kernels to include phasing terms for zenith # Phases are Nchan x Nstand x Npol # freq.shape += (1,) - phases = np.zeros((nchan, nstand, npol), dtype=np.complex64) for i in range(nstand): # X From a806b6bd44050b7d59b2459fa0c3a8262d0c72e0 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 2 Sep 2022 13:15:45 -0600 Subject: [PATCH 03/32] No 'z' in np.load. --- LWA_EPIC/LWA_EPIC.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 047afcb..baa3db8 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -892,7 +892,7 @@ def main(self): locname = "locations_%s_%i_%i_%i_%i_%i_%i_%.6f.npy" % (self.station.name, chan0, nchan, nstand, npol, self.ant_extent, self.grid_size, self.grid_resolution) locname = os.path.join(os.path.dirname(__file__), locname) try: - loc_data = np.loadz(locname) + loc_data = np.load(locname) sampling_length = loc_data['sampling_length'].item() locs = loc_data['locs'][...] sll = loc_data['sll'].item() From a7e19e8d361fa23b77015e8f6d245aaa393fc7e6 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 2 Sep 2022 13:18:16 -0600 Subject: [PATCH 04/32] Convert to a bifrost.ndarray after we load. --- LWA_EPIC/LWA_EPIC.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index baa3db8..6a5029c 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -988,8 +988,9 @@ def main(self): ): phases[:, :, i, :, :, :] = 0.0 phases = phases.conj() - phases = bifrost.ndarray(phases) np.save(phasename, phases) + phases = bifrost.ndarray(phases) + try: copy_array(gphases, phases) except NameError: From 5d4192ffffbd1d8db4dc39eeefffa37e1c0c07c5 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 2 Sep 2022 13:21:06 -0600 Subject: [PATCH 05/32] np.savez leads to .npz. --- LWA_EPIC/LWA_EPIC.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 6a5029c..ebf2008 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -889,7 +889,7 @@ def main(self): itshape = (self.ntime_gulp, nchan, nstand, npol) freq = (chan0 + np.arange(nchan)) * CHAN_BW - locname = "locations_%s_%i_%i_%i_%i_%i_%i_%.6f.npy" % (self.station.name, chan0, nchan, nstand, npol, self.ant_extent, self.grid_size, self.grid_resolution) + locname = "locations_%s_%i_%i_%i_%i_%i_%i_%.6f.npz" % (self.station.name, chan0, nchan, nstand, npol, self.ant_extent, self.grid_size, self.grid_resolution) locname = os.path.join(os.path.dirname(__file__), locname) try: loc_data = np.load(locname) From de5b15c1b9edc6ae39aebc4fd51ccad434ab4f27 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 2 Sep 2022 13:33:07 -0600 Subject: [PATCH 06/32] Lower buffer factors for DecimationOp and MOFFCorrelatorOp. --- LWA_EPIC/LWA_EPIC.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index ebf2008..b93ad2a 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -746,8 +746,8 @@ def main(self): ishape = (self.ntime_gulp, nchan, nstand, npol) ogulp_size = self.ntime_gulp * self.nchan_out * nstand * self.npol_out * 1 # ci4 oshape = (self.ntime_gulp, self.nchan_out, nstand, self.npol_out) - self.iring.resize(igulp_size, buffer_factor= 8) - self.oring.resize(ogulp_size, buffer_factor= 128) # , obuf_size) + self.iring.resize(igulp_size, buffer_factor= 5) + self.oring.resize(ogulp_size, buffer_factor= 10) # , obuf_size) ohdr = ihdr.copy() ohdr["nchan"] = self.nchan_out @@ -998,8 +998,8 @@ def main(self): oshape = (1, nchan, npol ** 2, self.grid_size, self.grid_size) ogulp_size = nchan * npol ** 2 * self.grid_size * self.grid_size * 8 - self.iring.resize(igulp_size, buffer_factor=128) - self.oring.resize(ogulp_size, buffer_factor=256) + self.iring.resize(igulp_size, buffer_factor=10) + self.oring.resize(ogulp_size, buffer_factor=10) prev_time = time.time() with oring.begin_sequence(time_tag=iseq.time_tag, header=ohdr_str) as oseq: iseq_spans = iseq.read(igulp_size) From 6d8072f77f4e476527161d162a39a404cf52b092 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 2 Sep 2022 14:15:57 -0600 Subject: [PATCH 07/32] Also include the gulp size. --- LWA_EPIC/LWA_EPIC.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index b93ad2a..3d9ba10 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -889,7 +889,7 @@ def main(self): itshape = (self.ntime_gulp, nchan, nstand, npol) freq = (chan0 + np.arange(nchan)) * CHAN_BW - locname = "locations_%s_%i_%i_%i_%i_%i_%i_%.6f.npz" % (self.station.name, chan0, nchan, nstand, npol, self.ant_extent, self.grid_size, self.grid_resolution) + locname = "locations_%s_%i_%i_%i_%i_%i_%i_%i_%.6f.npz" % (self.station.name, chan0, self.ntime_gulp, nchan, nstand, npol, self.ant_extent, self.grid_size, self.grid_resolution) locname = os.path.join(os.path.dirname(__file__), locname) try: loc_data = np.load(locname) @@ -951,7 +951,7 @@ def main(self): # Setup the kernels to include phasing terms for zenith # Phases are Ntime x Nchan x Nstand x Npol x extent x extent - phasename = "phases_%s_%i_%i_%i_%i_%i.npy" % (self.station.name, chan0, nchan, nstand, npol, self.ant_extent) + phasename = "phases_%s_%i_%i_%i_%i_%i_%i.npy" % (self.station.name, chan0, self.ntime_gulp, nchan, nstand, npol, self.ant_extent) phasename = os.path.join(os.path.dirname(__file__), phasename) try: phases = np.load(phasename) From f14db098ab0eac3d0ea5a35fbbe918a8ad2514b5 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 10:26:09 -0600 Subject: [PATCH 08/32] Formatting. --- LWA_EPIC/LWA_EPIC.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 3d9ba10..7978cb2 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -1174,7 +1174,7 @@ def main(self): bf_auto.init(self.locs, polmajor=False) bf_auto.execute(udata, autocorrs) autocorrs = autocorrs.reshape( - self.ntime_gulp, nchan, nstand, npol ** 2 + self.ntime_gulp, nchan, nstand, npol ** 2 ) From 7c9dfd40cf0a077877f89a3fe14d6cdf3fc38e04 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 10:28:06 -0600 Subject: [PATCH 09/32] Put the caches locations/phases into a 'cache' directory. --- LWA_EPIC/LWA_EPIC.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 7978cb2..5ae2fcb 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -828,7 +828,11 @@ def __init__( elif self.station == lwa1: locations[[i for i, a in enumerate(self.station.antennas[::2]) if a.stand.id in (35, 257, 258, 259, 260)], :] = 0.0 self.locations = locations - + + self.cache_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), 'cache') + if not os.path.exists(self.cache_dir): + os.mkdir(self.cache_dir) + self.grid_size = grid_size self.grid_resolution = grid_resolution @@ -890,7 +894,7 @@ def main(self): freq = (chan0 + np.arange(nchan)) * CHAN_BW locname = "locations_%s_%i_%i_%i_%i_%i_%i_%i_%.6f.npz" % (self.station.name, chan0, self.ntime_gulp, nchan, nstand, npol, self.ant_extent, self.grid_size, self.grid_resolution) - locname = os.path.join(os.path.dirname(__file__), locname) + locname = os.path.join(self.cache_dir, locname) try: loc_data = np.load(locname) sampling_length = loc_data['sampling_length'].item() @@ -952,7 +956,7 @@ def main(self): # Setup the kernels to include phasing terms for zenith # Phases are Ntime x Nchan x Nstand x Npol x extent x extent phasename = "phases_%s_%i_%i_%i_%i_%i_%i.npy" % (self.station.name, chan0, self.ntime_gulp, nchan, nstand, npol, self.ant_extent) - phasename = os.path.join(os.path.dirname(__file__), phasename) + phasename = os.path.join(self.cache_dir, phasename) try: phases = np.load(phasename) except OSError: From 7cf290536779e7b59a8dae590ba119e41e86f854 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 10:30:36 -0600 Subject: [PATCH 10/32] Remove some comments plus some more cleanup. --- LWA_EPIC/LWA_EPIC.py | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 5ae2fcb..7293118 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -1065,15 +1065,9 @@ def main(self): bf_vgrid = VGrid() bf_vgrid.init(self.locs, gphases, self.grid_size, polmajor=False) bf_vgrid.execute(udata, gdata) - - #try: - # bf_romein.execute(udata, gdata) - #except NameError: - # bf_romein = Romein() - # bf_romein.init(self.locs, gphases, self.grid_size, polmajor=False) - # bf_romein.execute(udata, gdata) - gdata = gdata.reshape(self.ntime_gulp * nchan * npol, self.grid_size, self.grid_size) - # gdata = self.LinAlgObj.matmul(1.0, udata, bfantgridmap, 0.0, gdata) + gdata = gdata.reshape( + self.ntime_gulp * nchan * npol, self.grid_size, self.grid_size + ) if self.benchmark is True: timeg2 = time.time() print(" Grid time: %f" % (timeg2 - timeg1)) @@ -1097,7 +1091,6 @@ def main(self): timefft2 = time.time() print(" FFT time: %f" % (timefft2 - timefft1)) - # print("Accum: %d"%accum,end='\n') if self.newflag is True: try: crosspol = crosspol.reshape( From 7092910694ec926f6d3cab49a1b824e92946a1c4 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 10:35:17 -0600 Subject: [PATCH 11/32] Division fix. --- LWA_EPIC/LWA_EPIC.py | 20 +++++--------------- 1 file changed, 5 insertions(+), 15 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 7293118..aafa3ee 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -1060,8 +1060,10 @@ def main(self): if self.benchmark is True: timeg1 = time.time() try: + bf_vgrid.execute(udata, gdata) except NameError: + bf_vgrid = VGrid() bf_vgrid.init(self.locs, gphases, self.grid_size, polmajor=False) bf_vgrid.execute(udata, gdata) @@ -1146,7 +1148,7 @@ def main(self): np.ones( shape=(3, 1, nchan, nstand, npol ** 2), dtype=np.int32 - ) * self.grid_size / 2, + ) * self.grid_size // 2, space="cuda", ) autocorr_il = bifrost.ndarray( @@ -1157,13 +1159,6 @@ def main(self): space="cuda", ) - # Cross multiply to calculate autocorrs - #bifrost.map( - # "a(i,j,k,l) += (b(i,j,k,l/2) * b(i,j,k,l%2).conj())", - # {"a": autocorrs, "b": udata, "t": self.ntime_gulp}, - # axis_names=("i", "j", "k", "l"), - # shape=(self.ntime_gulp, nchan, nstand, npol ** 2), - #) try: bf_auto.execute(udata, autocorrs) except NameError: @@ -1174,16 +1169,11 @@ def main(self): self.ntime_gulp, nchan, nstand, npol ** 2 ) - - #bifrost.map( - # "a(i,j,p,k,l) += b(0,i,j,p/2,k,l)*b(0,i,j,p%2,k,l).conj()", - # {"a": crosspol, "b": gdata}, - # axis_names=("i", "j", "p", "k", "l"), - # shape=(self.ntime_gulp, nchan, npol ** 2, self.grid_size, self.grid_size), - #) try: + bf_gmul.execute(gdata, crosspol) except NameError: + bf_gmul = XGrid() bf_gmul.init(self.grid_size, polmajor=False) bf_gmul.execute(gdata, crosspol) From 0966222ff35b16a51ef05fb7ef3c86fde318c679 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 11:11:15 -0600 Subject: [PATCH 12/32] Tweak DecimationOp to do channel averaging when possible to preserve overall bandwidth. --- LWA_EPIC/LWA_EPIC.py | 33 +++++++++++++++++++++++++++------ 1 file changed, 27 insertions(+), 6 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index aafa3ee..6789bf4 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -550,7 +550,7 @@ def main(self): ohdr["chan0"] = chan0 ohdr["nchan"] = nchan ohdr["cfreq"] = (chan0 + 0.5 * (nchan - 1)) * srate - ohdr["bw"] = srate * nchan + ohdr["bw"] = nchan * srate ohdr["nstand"] = nstand ohdr["npol"] = npol ohdr["nbit"] = 4 @@ -749,11 +749,21 @@ def main(self): self.iring.resize(igulp_size, buffer_factor= 5) self.oring.resize(ogulp_size, buffer_factor= 10) # , obuf_size) + do_truncate = True + act_chan_bw = CHAN_BW + if nchan % self.nchan_out == 0: + do_truncate = False + act_chan_bw = CHAN_BW * (nchan // self.nchan_out) + self.log.info("Decimation: Running in averageing mode") + else: + self.log.info("Decimation: Running in truncation mode") + self.log.info("Decimation: Channel bandwidth is %.3f kHz", act_chan_bw/1e3) + ohdr = ihdr.copy() ohdr["nchan"] = self.nchan_out ohdr["npol"] = self.npol_out - ohdr["cfreq"] = (chan0 + 0.5 * (self.nchan_out - 1)) * CHAN_BW - ohdr["bw"] = self.nchan_out * CHAN_BW + ohdr["cfreq"] = chan0 * CHAN_BW + 0.5 * (self.nchan_out - 1) * act_chan_bw + ohdr["bw"] = self.nchan_out * act_chan_bw ohdr_str = json.dumps(ohdr) prev_time = time.time() @@ -773,7 +783,14 @@ def main(self): idata = ispan.data_view(np.uint8).reshape(ishape) odata = ospan.data_view(np.uint8).reshape(oshape) - sdata = idata[:, :self.nchan_out, :, :] + if do_truncate: + sdata = idata[:, :self.nchan_out, :, :] + else: + sdata = sdata.reshape( + self.ntime_gulp, -1, nchan // self.nchan_out, nstand, npol + ) + sdata = sdate.mean(axis=2) + if self.npol_out != npol: sdata = sdata[:, :, :, :self.npol_out] odata[...] = sdata @@ -884,6 +901,8 @@ def main(self): self.log.info("MOFFCorrelatorOp: Config - %s" % ihdr) chan0 = ihdr["chan0"] nchan = ihdr["nchan"] + bw = ihdr["bw"] + act_chan_bw = bw / nchan nstand = ihdr["nstand"] npol = ihdr["npol"] self.newflag = True @@ -892,7 +911,7 @@ def main(self): igulp_size = self.ntime_gulp * nchan * nstand * npol * 1 # ci4 itshape = (self.ntime_gulp, nchan, nstand, npol) - freq = (chan0 + np.arange(nchan)) * CHAN_BW + freq = chan0 * CHAN_BW + np.arange(nchan) * act_chan_bw locname = "locations_%s_%i_%i_%i_%i_%i_%i_%i_%.6f.npz" % (self.station.name, chan0, self.ntime_gulp, nchan, nstand, npol, self.ant_extent, self.grid_size, self.grid_resolution) locname = os.path.join(self.cache_dir, locname) try: @@ -1386,6 +1405,8 @@ def main(self): nchan = ihdr["nchan"] nstand = ihdr["nstand"] npol = ihdr["npol"] + bw = ihdr["bw"] + act_chan_bw = bw / nchan self.newflag = True accum = 0 @@ -1393,7 +1414,7 @@ def main(self): itshape = (self.ntime_gulp, nchan, nstand, npol) # Sample locations at right u/v/w values - freq = (chan0 + np.arange(nchan)) * CHAN_BW + freq = chan0 * CHAN_BW + np.arange(nchan) * act_chan_bw locs = Generate_DFT_Locations( self.locations, freq, self.ntime_gulp, nchan, npol ) From 3a9dd03b3f5294134572d45f1c6a2acaaec062e3 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 11:13:23 -0600 Subject: [PATCH 13/32] Make a copy. --- LWA_EPIC/LWA_EPIC.py | 1 + 1 file changed, 1 insertion(+) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 6789bf4..730b72c 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -786,6 +786,7 @@ def main(self): if do_truncate: sdata = idata[:, :self.nchan_out, :, :] else: + sdata = idata[...] sdata = sdata.reshape( self.ntime_gulp, -1, nchan // self.nchan_out, nstand, npol ) From 8abee0122ecaaad2039c89cc0a617d413041ef70 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 11:16:40 -0600 Subject: [PATCH 14/32] Make a copy. --- LWA_EPIC/LWA_EPIC.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 730b72c..77e79be 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -786,11 +786,10 @@ def main(self): if do_truncate: sdata = idata[:, :self.nchan_out, :, :] else: - sdata = idata[...] - sdata = sdata.reshape( + sdata = idata.reshape( self.ntime_gulp, -1, nchan // self.nchan_out, nstand, npol ) - sdata = sdate.mean(axis=2) + sdata = sdata.mean(axis=2) if self.npol_out != npol: sdata = sdata[:, :, :, :self.npol_out] From 8af2e6da06dcdfa1a1ced260d79904c79b116340 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 11:56:21 -0600 Subject: [PATCH 15/32] This should be a better way to average channels together. --- LWA_EPIC/LWA_EPIC.py | 37 +++++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 77e79be..b5d21a2 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -786,10 +786,36 @@ def main(self): if do_truncate: sdata = idata[:, :self.nchan_out, :, :] else: - sdata = idata.reshape( + # Fix the type + udata = bifrost.ndarray( + shape=itshape, + dtype="ci4", + native=False, + buffer=idata.ctypes.data, + ) + + try: + Unpack(udata, cdata) + except NameError: + cdata = bifrost.zeros( + shape=idata.shape + dtype=np.complex64 + ) + Unpack(udata, cdata) + + cdata = cdata.reshape( self.ntime_gulp, -1, nchan // self.nchan_out, nstand, npol ) - sdata = sdata.mean(axis=2) + cdata = cdata.mean(axis=2, dtype=np.complex64) + + try: + Quantize(cdata, sdata) + except NameError: + qdata = bifrost.zeros( + shape=cdata.shape + dtype='ci4' + ) + Quantize(cdata, sdata) if self.npol_out != npol: sdata = sdata[:, :, :, :self.npol_out] @@ -806,6 +832,13 @@ def main(self): } ) + if not do_truncate: + try: + del cdata + del sdata + except NameError: + pass + class CalibrationOp(object): From df5f363ab3dc7d04dfb5a353f8a31c8e25ed6870 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 11:56:57 -0600 Subject: [PATCH 16/32] Ugh, typo. --- LWA_EPIC/LWA_EPIC.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index b5d21a2..639d377 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -798,7 +798,7 @@ def main(self): Unpack(udata, cdata) except NameError: cdata = bifrost.zeros( - shape=idata.shape + shape=idata.shape, dtype=np.complex64 ) Unpack(udata, cdata) @@ -812,7 +812,7 @@ def main(self): Quantize(cdata, sdata) except NameError: qdata = bifrost.zeros( - shape=cdata.shape + shape=cdata.shape, dtype='ci4' ) Quantize(cdata, sdata) From 6e038254e68b92f71b0c9b4e6177b403498ff59a Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 12:54:47 -0600 Subject: [PATCH 17/32] Several fixes plus move the averaging over to the GPU. --- LWA_EPIC/LWA_EPIC.py | 35 ++++++++++++++++++++--------------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 639d377..af307e7 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -754,7 +754,7 @@ def main(self): if nchan % self.nchan_out == 0: do_truncate = False act_chan_bw = CHAN_BW * (nchan // self.nchan_out) - self.log.info("Decimation: Running in averageing mode") + self.log.info("Decimation: Running in averaging mode") else: self.log.info("Decimation: Running in truncation mode") self.log.info("Decimation: Channel bandwidth is %.3f kHz", act_chan_bw/1e3) @@ -784,42 +784,47 @@ def main(self): odata = ospan.data_view(np.uint8).reshape(oshape) if do_truncate: - sdata = idata[:, :self.nchan_out, :, :] + sdata2 = idata[:, :self.nchan_out, :, :] else: # Fix the type udata = bifrost.ndarray( - shape=itshape, + shape=ishape, dtype="ci4", native=False, buffer=idata.ctypes.data, ) + udata = udata.copy(space='cuda') try: - Unpack(udata, cdata) + cdata = cdata.reshape( + idata.shape + ) except NameError: cdata = bifrost.zeros( shape=idata.shape, - dtype=np.complex64 + dtype=np.complex64, + space='cuda' ) - Unpack(udata, cdata) - + Unpack(udata, cdata) cdata = cdata.reshape( self.ntime_gulp, -1, nchan // self.nchan_out, nstand, npol ) - cdata = cdata.mean(axis=2, dtype=np.complex64) + adata = cdata.mean(axis=2, dtype=np.complex64) try: - Quantize(cdata, sdata) + Quantize(adata, sdata) except NameError: - qdata = bifrost.zeros( - shape=cdata.shape, - dtype='ci4' + sdata = bifrost.zeros( + shape=adata.shape, + dtype='ci4', + space='cuda' ) - Quantize(cdata, sdata) + Quantize(adata, sdata) + sdata2 = sdata.copy(space='system') if self.npol_out != npol: - sdata = sdata[:, :, :, :self.npol_out] - odata[...] = sdata + sdata3 = sdata2[:, :, :, :self.npol_out] + odata[...] = sdata3.view(np.uint8) curr_time = time.time() process_time = curr_time - prev_time From 357611c62a875cf79c0f55463cdd1c97744d9a1f Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 14:49:39 -0600 Subject: [PATCH 18/32] A better way to do the averaging on the GPU. --- LWA_EPIC/LWA_EPIC.py | 78 +++++++++++++++++++++++++------------------- 1 file changed, 44 insertions(+), 34 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index af307e7..f021b64 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -701,6 +701,7 @@ def __init__( npol_out=2, guarantee=True, core=-1, + gpu=-1, ): self.log = log self.iring = iring @@ -710,6 +711,7 @@ def __init__( self.npol_out = npol_out self.guarantee = guarantee self.core = core + self.gpu = gpu self.bind_proclog = ProcLog(type(self).__name__ + "/bind") self.in_proclog = ProcLog(type(self).__name__ + "/in") @@ -725,8 +727,15 @@ def __init__( def main(self): if self.core != -1: bifrost.affinity.set_core(self.core) + if self.gpu != -1: + BFSetGPU(self.gpu) self.bind_proclog.update( - {"ncore": 1, "core0": bifrost.affinity.get_core()} + { + "ncore": 1, + "core0": bifrost.affinity.get_core(), + "ngpu": 1, + "gpu0": BFGetGPU(), + } ) with self.oring.begin_writing() as oring: @@ -785,46 +794,45 @@ def main(self): if do_truncate: sdata2 = idata[:, :self.nchan_out, :, :] + if self.npol_out != npol: + sdata = sdata[:, :, :, :self.npol_out] else: - # Fix the type - udata = bifrost.ndarray( - shape=ishape, - dtype="ci4", - native=False, - buffer=idata.ctypes.data, - ) - udata = udata.copy(space='cuda') - try: - cdata = cdata.reshape( - idata.shape - ) + copy_array(gdata, idata) except NameError: - cdata = bifrost.zeros( - shape=idata.shape, - dtype=np.complex64, + gdata = idata.copy(space='cuda') + adata = bifrost.zeros( + shape=(self.ntime_gulp, self.nchan_out, nstand, self.npol_out), + dtype='u8', space='cuda' ) - Unpack(udata, cdata) - cdata = cdata.reshape( - self.ntime_gulp, -1, nchan // self.nchan_out, nstand, npol - ) - adata = cdata.mean(axis=2, dtype=np.complex64) - + + bifrost.map(""" + signed char re, im; + #pragma unroll + for(int l=0; l<{npol_out}; l++) {{ + re = ((signed char) (b(i,j*{navg},k,l) & 0xF0)) / 16; + im = ((signed char) ((b(i,j*{navg},k,l) & 0x0F) << 4)) / 16; + #pragma unroll + for(int m=1; m<{navg}; m++) {{ + re += ((signed char) (b(i,j*{navg}+m,k,l) & 0xF0)) / 16; + im += ((signed char) ((b(i,j*{navg}+m,k,l) & 0x0F) << 4)) / 16; + }} + re /= {navg}; + im /= {navg}; + a(i,j,k,l) = ((re * 16) & 0xF0) | ((im * 16) >> 4); + }} + """.format(npol_out=self.npol_out, navg=nchan // self.nchan_out), + {'a': adata, 'b': gdata}, + axis_names=('i', 'j', 'k'), + shape=adata.shape[:3]) + try: - Quantize(adata, sdata) + copy_array(sdata, adata) except NameError: - sdata = bifrost.zeros( - shape=adata.shape, - dtype='ci4', - space='cuda' - ) - Quantize(adata, sdata) - sdata2 = sdata.copy(space='system') + sdata = adata.copy(space='system') - if self.npol_out != npol: - sdata3 = sdata2[:, :, :, :self.npol_out] - odata[...] = sdata3.view(np.uint8) + odata[...] = sdata curr_time = time.time() process_time = curr_time - prev_time @@ -2178,7 +2186,7 @@ def gen_args(return_parser=False): group2.add_argument("--tbffile", type=str, help="TBF Data Path") group3 = parser.add_argument_group("Processing Options") - group3.add_argument("--cores", type=str, default="3,4,5,6,7", help="Comma separated list of CPU cores to bind to") + group3.add_argument("--cores", type=str, default="2,3,4,5,6,7", help="Comma separated list of CPU cores to bind to") group3.add_argument("--gpu", type=int, default=0, help="GPU to bind to") group3.add_argument("--imagesize", type=int, default=64, help="1-D Image Size") group3.add_argument( @@ -2387,6 +2395,7 @@ def handle_signal_terminate(signum, frame): nchan_out=args.channels, npol_out=1 if args.singlepol else 2, core=cores.pop(0), + gpu=gpus.pop(0), ) ) else: @@ -2430,6 +2439,7 @@ def handle_signal_terminate(signum, frame): nchan_out=args.channels, npol_out=1 if args.singlepol else 2, core=cores.pop(0), + gpu=gpus.pop(0), ) ) From a69f84d36a25cda38f6857b26dc3ac0e295767f3 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 14:55:14 -0600 Subject: [PATCH 19/32] Yet another typo. --- LWA_EPIC/LWA_EPIC.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index f021b64..3ca1730 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -793,7 +793,7 @@ def main(self): odata = ospan.data_view(np.uint8).reshape(oshape) if do_truncate: - sdata2 = idata[:, :self.nchan_out, :, :] + sdata = idata[:, :self.nchan_out, :, :] if self.npol_out != npol: sdata = sdata[:, :, :, :self.npol_out] else: From 0196256414c0b7c577c642b5342fa673827bc677 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Tue, 6 Sep 2022 15:03:15 -0600 Subject: [PATCH 20/32] chan0 adjustment in DecimationOp. --- LWA_EPIC/LWA_EPIC.py | 1 + 1 file changed, 1 insertion(+) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 3ca1730..ca26015 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -763,6 +763,7 @@ def main(self): if nchan % self.nchan_out == 0: do_truncate = False act_chan_bw = CHAN_BW * (nchan // self.nchan_out) + chan0 = chan0 + 0.5 * (nchan // self.nchan_out - 1) self.log.info("Decimation: Running in averaging mode") else: self.log.info("Decimation: Running in truncation mode") From c6b2ec40b653bf63e10f973877bc588ada10f6ba Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Thu, 8 Sep 2022 01:56:37 -0600 Subject: [PATCH 21/32] Update chan0 inside DecimationOp. --- LWA_EPIC/LWA_EPIC.py | 1 + 1 file changed, 1 insertion(+) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index ca26015..9fb44cd 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -770,6 +770,7 @@ def main(self): self.log.info("Decimation: Channel bandwidth is %.3f kHz", act_chan_bw/1e3) ohdr = ihdr.copy() + ohdr["chan0"] = chan0 ohdr["nchan"] = self.nchan_out ohdr["npol"] = self.npol_out ohdr["cfreq"] = chan0 * CHAN_BW + 0.5 * (self.nchan_out - 1) * act_chan_bw From c9ce4945d689ae39c7a3ae350e2ac19dd410c9dd Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 28 Oct 2022 19:00:31 -0600 Subject: [PATCH 22/32] Have DecimationOp output ci8 values instead of ci4 to help preserve the phase. --- LWA_EPIC/LWA_EPIC.py | 130 +++++++++++++++++++------------------------ 1 file changed, 58 insertions(+), 72 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 9fb44cd..6943387 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -427,7 +427,7 @@ def main(self): igulp_size = self.ntime_gulp * 1 * nstand * npol * 8 # complex64 ishape = (self.ntime_gulp // nchan, nchan, nstand, npol) - ogulp_size = self.ntime_gulp * 1 * nstand * npol * 1 # ci4 + ogulp_size = self.ntime_gulp * 1 * nstand * npol * 2 # ci8 oshape = (self.ntime_gulp // nchan, nchan, nstand, npol) # self.iring.resize(igulp_size) self.oring.resize(ogulp_size, buffer_factor=5) @@ -462,7 +462,7 @@ def main(self): # Setup and load idata = ispan.data_view(np.complex64).reshape(ishape) - odata = ospan.data_view(np.int8).reshape(oshape) + odata = ospan.data_view('ci8').reshape(oshape) # FFT, shift, and phase fdata = fft(idata, axis=1) @@ -473,11 +473,11 @@ def main(self): try: Quantize(fdata, qdata, scale=1. / np.sqrt(nchan)) except NameError: - qdata = bifrost.ndarray(shape=fdata.shape, native=False, dtype="ci4") + qdata = bifrost.ndarray(shape=fdata.shape, dtype="ci8") Quantize(fdata, qdata, scale=1. / np.sqrt(nchan)) # Save - odata[...] = qdata.copy(space="cuda_host").view(np.int8).reshape(oshape) + odata[...] = qdata.copy(space="cuda_host").reshape(oshape) if self.profile: spani += 1 @@ -753,7 +753,7 @@ def main(self): igulp_size = self.ntime_gulp * nchan * nstand * npol * 1 # ci4 ishape = (self.ntime_gulp, nchan, nstand, npol) - ogulp_size = self.ntime_gulp * self.nchan_out * nstand * self.npol_out * 1 # ci4 + ogulp_size = self.ntime_gulp * self.nchan_out * nstand * self.npol_out * 2 # ci8 oshape = (self.ntime_gulp, self.nchan_out, nstand, self.npol_out) self.iring.resize(igulp_size, buffer_factor= 5) self.oring.resize(ogulp_size, buffer_factor= 10) # , obuf_size) @@ -762,10 +762,12 @@ def main(self): act_chan_bw = CHAN_BW if nchan % self.nchan_out == 0: do_truncate = False + navg = nchan // self.nchan_out act_chan_bw = CHAN_BW * (nchan // self.nchan_out) chan0 = chan0 + 0.5 * (nchan // self.nchan_out - 1) self.log.info("Decimation: Running in averaging mode") else: + navg = 1 self.log.info("Decimation: Running in truncation mode") self.log.info("Decimation: Channel bandwidth is %.3f kHz", act_chan_bw/1e3) @@ -791,51 +793,45 @@ def main(self): reserve_time = curr_time - prev_time prev_time = curr_time - idata = ispan.data_view(np.uint8).reshape(ishape) - odata = ospan.data_view(np.uint8).reshape(oshape) + idata = ispan.data_view('ci4').reshape(ishape) + odata = ospan.data_view('ci8').reshape(oshape) - if do_truncate: - sdata = idata[:, :self.nchan_out, :, :] - if self.npol_out != npol: - sdata = sdata[:, :, :, :self.npol_out] - else: - try: - copy_array(gdata, idata) - except NameError: - gdata = idata.copy(space='cuda') - adata = bifrost.zeros( - shape=(self.ntime_gulp, self.nchan_out, nstand, self.npol_out), - dtype='u8', - space='cuda' - ) - - bifrost.map(""" - signed char re, im; - #pragma unroll - for(int l=0; l<{npol_out}; l++) {{ - re = ((signed char) (b(i,j*{navg},k,l) & 0xF0)) / 16; - im = ((signed char) ((b(i,j*{navg},k,l) & 0x0F) << 4)) / 16; - #pragma unroll - for(int m=1; m<{navg}; m++) {{ - re += ((signed char) (b(i,j*{navg}+m,k,l) & 0xF0)) / 16; - im += ((signed char) ((b(i,j*{navg}+m,k,l) & 0x0F) << 4)) / 16; - }} - re /= {navg}; - im /= {navg}; - a(i,j,k,l) = ((re * 16) & 0xF0) | ((im * 16) >> 4); - }} - """.format(npol_out=self.npol_out, navg=nchan // self.nchan_out), - {'a': adata, 'b': gdata}, - axis_names=('i', 'j', 'k'), - shape=adata.shape[:3]) - - try: - copy_array(sdata, adata) - except NameError: - sdata = adata.copy(space='system') - - odata[...] = sdata + try: + copy_array(gdata, idata) + except NameError: + gdata = idata.copy(space='cuda') + adata = bifrost.zeros( + shape=(self.ntime_gulp, self.nchan_out, nstand, self.npol_out), + dtype='ci4', + space='cuda' + ) + bifrost.map(""" + int jF; + signed char sample, re, im; + + #pragma unroll + for(int l=0; l<{npol_out}; l++) {{ + re = im = 0; + + #pragma unroll + for(int m=0; m<{navg}; m++) {{ + jF = j*DECIM + m; + sample = b(i,jF,k,l).real_imag; + re += ((signed char) (sample & 0xF0)) / 16; + im += ((signed char) ((sample & 0x0F) << 4)) / 16; + }} + + // Save + a(i,j,k,l) = Complex(re, im); + }} + """.format(npol_out=self.npol_out, navg=navg), + {'a': adata, 'b': gdata}, + axis_names=('i', 'j', 'k'), + shape=adata.shape[:3]) + + copy_array(odata, adata) + curr_time = time.time() process_time = curr_time - prev_time prev_time = curr_time @@ -847,12 +843,11 @@ def main(self): } ) - if not do_truncate: - try: - del cdata - del sdata - except NameError: - pass + try: + del gdata + del adata + except NameError: + pass @@ -956,7 +951,7 @@ def main(self): self.newflag = True accum = 0 - igulp_size = self.ntime_gulp * nchan * nstand * npol * 1 # ci4 + igulp_size = self.ntime_gulp * nchan * nstand * npol * 2 # ci8 itshape = (self.ntime_gulp, nchan, nstand, npol) freq = chan0 * CHAN_BW + np.arange(nchan) * act_chan_bw @@ -1092,19 +1087,12 @@ def main(self): # Correlator # Setup and load - idata = ispan.data_view(np.uint8).reshape(itshape) - # Fix the type - udata = bifrost.ndarray( - shape=itshape, - dtype="ci4", - native=False, - buffer=idata.ctypes.data, - ) + idata = ispan.data_view('ci8').reshape(itshape) if self.benchmark is True: time1 = time.time() - udata = udata.copy(space="cuda") + udata = idata.copy(space="cuda") if self.benchmark is True: time1a = time.time() print(" Input copy time: %f" % (time1a - time1)) @@ -1458,7 +1446,7 @@ def main(self): self.newflag = True accum = 0 - igulp_size = self.ntime_gulp * nchan * nstand * npol * 1 # ci4 + igulp_size = self.ntime_gulp * nchan * nstand * npol * 2 # ci8 itshape = (self.ntime_gulp, nchan, nstand, npol) # Sample locations at right u/v/w values @@ -1579,16 +1567,14 @@ def main(self): # Correlator # Setup and load - idata = ispan.data_view(np.uint8).reshape(itshape) - # Fix the type - tdata = bifrost.ndarray(shape=itshape, dtype="ci4", native=False, buffer=idata.ctypes.data) - + idata = ispan.data_view('ci8').reshape(itshape) + if self.benchmark is True: time1 = time.time() # chan pol stand - tdata = tdata.transpose((1, 3, 2, 0)) - tdata = tdata.reshape(nchan * npol, nstand, self.ntime_gulp) - tdata = tdata.copy() + idata = idata.transpose((1, 3, 2, 0)) + idata = idata.reshape(nchan * npol, nstand, self.ntime_gulp) + tdata = idata.copy() tdata = tdata.copy(space="cuda") # Unpack From 237e29c9e0c1aa566a363583bc1739e102369e99 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 28 Oct 2022 19:06:34 -0600 Subject: [PATCH 23/32] Add in a new 'acphases' array for the autocorrelation removal so that cable gains are taken into account. --- LWA_EPIC/LWA_EPIC.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 6943387..6adcd26 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -763,7 +763,7 @@ def main(self): if nchan % self.nchan_out == 0: do_truncate = False navg = nchan // self.nchan_out - act_chan_bw = CHAN_BW * (nchan // self.nchan_out) + act_chan_bw = CHAN_BW * navg chan0 = chan0 + 0.5 * (nchan // self.nchan_out - 1) self.log.info("Decimation: Running in averaging mode") else: @@ -1055,12 +1055,16 @@ def main(self): phases[:, :, i, :, :, :] = 0.0 phases = phases.conj() np.save(phasename, phases) + acphases = np.abs(phases, dtype=np.complex64) phases = bifrost.ndarray(phases) + acphases = bifrost.ndarray(acphases) try: copy_array(gphases, phases) + copy_array(gacphases, acphases) except NameError: gphases = phases.copy(space="cuda") + gacphases = acphases.copy(space="cuda") oshape = (1, nchan, npol ** 2, self.grid_size, self.grid_size) ogulp_size = nchan * npol ** 2 * self.grid_size * self.grid_size * 8 @@ -1206,13 +1210,6 @@ def main(self): ) * self.grid_size // 2, space="cuda", ) - autocorr_il = bifrost.ndarray( - np.ones( - shape=(1, nchan, nstand, npol ** 2, self.ant_extent, self.ant_extent), - dtype=np.complex64 - ), - space="cuda", - ) try: bf_auto.execute(udata, autocorrs) @@ -1265,7 +1262,7 @@ def main(self): except NameError: bf_vgrid_autocorr = VGrid() bf_vgrid_autocorr.init( - autocorr_lo, autocorr_il, self.grid_size, polmajor=False + autocorr_lo, gacphases, self.grid_size, polmajor=False ) bf_vgrid_autocorr.execute(autocorrs_av, autocorr_g) autocorr_g = autocorr_g.reshape(1 * nchan * npol ** 2, self.grid_size, self.grid_size) From 84232ed24c1dbd001d2453152ff8dd47bab14ecf Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 28 Oct 2022 19:11:38 -0600 Subject: [PATCH 24/32] Revert to the Romein gridder for the auto-correlation removal. --- LWA_EPIC/LWA_EPIC.py | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 6adcd26..6c971d1 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -1249,22 +1249,22 @@ def main(self): autocorr_g = autocorr_g.reshape( 1, nchan, npol ** 2, self.grid_size, self.grid_size ) - #try: - # bf_romein_autocorr.execute(autocorrs_av, autocorr_g) - #except NameError: - # bf_romein_autocorr = Romein() - # bf_romein_autocorr.init( - # autocorr_lo, autocorr_il, self.grid_size, polmajor=False - # ) - # bf_romein_autocorr.execute(autocorrs_av, autocorr_g) try: - bf_vgrid_autocorr.execute(autocorrs_av, autocorr_g) + bf_romein_autocorr.execute(autocorrs_av, autocorr_g) except NameError: - bf_vgrid_autocorr = VGrid() - bf_vgrid_autocorr.init( + bf_romein_autocorr = Romein() + bf_romein_autocorr.init( autocorr_lo, gacphases, self.grid_size, polmajor=False ) - bf_vgrid_autocorr.execute(autocorrs_av, autocorr_g) + bf_romein_autocorr.execute(autocorrs_av, autocorr_g) + # try: + # bf_vgrid_autocorr.execute(autocorrs_av, autocorr_g) + # except NameError: + # bf_vgrid_autocorr = VGrid() + # bf_vgrid_autocorr.init( + # autocorr_lo, gacphases, self.grid_size, polmajor=False + # ) + # bf_vgrid_autocorr.execute(autocorrs_av, autocorr_g) autocorr_g = autocorr_g.reshape(1 * nchan * npol ** 2, self.grid_size, self.grid_size) # autocorr_g = romein_float(autocorrs_av,autocorr_g,autocorr_il,autocorr_lx,autocorr_ly,autocorr_lz,self.ant_extent,self.grid_size,nstand,nchan*npol**2) # Inverse FFT From aba6d6399dc77732f6c2f4b62570a53da56ff7f4 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 28 Oct 2022 19:14:40 -0600 Subject: [PATCH 25/32] Typo. --- LWA_EPIC/LWA_EPIC.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 6c971d1..8b4d8ed 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -816,7 +816,7 @@ def main(self): #pragma unroll for(int m=0; m<{navg}; m++) {{ - jF = j*DECIM + m; + jF = j*{avg} + m; sample = b(i,jF,k,l).real_imag; re += ((signed char) (sample & 0xF0)) / 16; im += ((signed char) ((sample & 0x0F) << 4)) / 16; From f12b16e8acbb1dabba883d047a89e2f644974459 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 28 Oct 2022 23:38:05 -0600 Subject: [PATCH 26/32] Typo. --- LWA_EPIC/LWA_EPIC.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 8b4d8ed..77fb1ee 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -816,7 +816,7 @@ def main(self): #pragma unroll for(int m=0; m<{navg}; m++) {{ - jF = j*{avg} + m; + jF = j*{navg} + m; sample = b(i,jF,k,l).real_imag; re += ((signed char) (sample & 0xF0)) / 16; im += ((signed char) ((sample & 0x0F) << 4)) / 16; From c2c2528fb3ed179262bb72cd93ece3e761ab12ca Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 28 Oct 2022 23:39:11 -0600 Subject: [PATCH 27/32] New type. --- LWA_EPIC/LWA_EPIC.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 77fb1ee..d617f1a 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -802,7 +802,7 @@ def main(self): gdata = idata.copy(space='cuda') adata = bifrost.zeros( shape=(self.ntime_gulp, self.nchan_out, nstand, self.npol_out), - dtype='ci4', + dtype='ci8', space='cuda' ) From af4ef385ecf320920c08771d3163756b56b9f5fc Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 28 Oct 2022 23:41:59 -0600 Subject: [PATCH 28/32] Fix the type on 'acphases'. --- LWA_EPIC/LWA_EPIC.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index d617f1a..ec40f40 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -1055,7 +1055,8 @@ def main(self): phases[:, :, i, :, :, :] = 0.0 phases = phases.conj() np.save(phasename, phases) - acphases = np.abs(phases, dtype=np.complex64) + acphases = np.abs(phases) + 0j + acphases = acphases.astype(np.complex64) phases = bifrost.ndarray(phases) acphases = bifrost.ndarray(acphases) From 3318446ea1d67cd780d2a6b34cf062eea98bd327 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Fri, 28 Oct 2022 23:44:14 -0600 Subject: [PATCH 29/32] Import Romein. --- LWA_EPIC/LWA_EPIC.py | 1 + 1 file changed, 1 insertion(+) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index ec40f40..fa99548 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -45,6 +45,7 @@ from bifrost.libbifrost import bf from bifrost.fft import Fft from bifrost.linalg import LinAlg +from bifrost.romein import Romein from bifrost.ndarray import memset_array, copy_array from bifrost.device import set_device as BFSetGPU, get_device as BFGetGPU, set_devices_no_spin_cpu as BFNoSpinZone From e2efe7a14dbf96740ef96256b0884080736b27f5 Mon Sep 17 00:00:00 2001 From: jaycedowell Date: Sat, 29 Oct 2022 00:00:46 -0600 Subject: [PATCH 30/32] 'acphases' has a different shape than I thought. --- LWA_EPIC/LWA_EPIC.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index fa99548..8eea85f 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -1056,7 +1056,12 @@ def main(self): phases[:, :, i, :, :, :] = 0.0 phases = phases.conj() np.save(phasename, phases) - acphases = np.abs(phases) + 0j + acphases = np.zeros( + (1, nchan, nstand, npol**2, self.ant_extent, self.ant_extent), + dtype=np.complex64 + ) + for i in range(npol**2): + acphases[:,:,:,i,:,:] = np.abs(phases[0,:,:,i//2,:,:]*phases[0,:,:,i%2,:,:].conj()) acphases = acphases.astype(np.complex64) phases = bifrost.ndarray(phases) acphases = bifrost.ndarray(acphases) From 23c521d8f179b4968eff64a589a04d64992e4992 Mon Sep 17 00:00:00 2001 From: Craig Taylor <73905895+ctaylor-physics@users.noreply.github.com> Date: Tue, 17 Oct 2023 13:37:36 -0600 Subject: [PATCH 31/32] Update LWA_EPIC.py Slew of updates from my work on the EPIC imaging comparisons. --- LWA_EPIC/LWA_EPIC.py | 70 ++++++++++++++++++++++++++------------------ 1 file changed, 42 insertions(+), 28 deletions(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index 8eea85f..e1369f1 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -21,6 +21,7 @@ from astropy.io import fits from astropy.time import Time, TimeDelta +from astropy import units as u from astropy.coordinates import SkyCoord, FK5 from astropy.constants import c as speed_of_light @@ -37,7 +38,7 @@ import bifrost.affinity from bifrost.address import Address as BF_Address from bifrost.udp_socket import UDPSocket as BF_UDPSocket -from bifrost.packet_capture import PacketCaptureCallback, UDPCapture +from bifrost.packet_capture import PacketCaptureCallback, UDPVerbsCapture as UDPCapture from bifrost.ring import Ring from bifrost.unpack import unpack as Unpack from bifrost.quantize import quantize as Quantize @@ -237,9 +238,9 @@ def GenerateLocations( lsl_locs = lsl_locs.T lsl_locsf = lsl_locs[:, np.newaxis, np.newaxis, :] / sample_grid[np.newaxis, np.newaxis, :, np.newaxis] - lsl_locsf -= np.min(lsl_locsf, axis=3, keepdims=True) - + # Centre locations slightly + lsl_locsf -= np.min(lsl_locsf, axis=3, keepdims=True) lsl_locsf += (grid_size - np.max(lsl_locsf, axis=3, keepdims=True)) / 2. # add ntime axis @@ -537,7 +538,7 @@ def main(self): # Setup the ring metadata and gulp sizes ntime = data.shape[2] - nstand, npol = data.shape[0] / 2, 2 + nstand, npol = data.shape[0] // 2, 2 oshape = (ntime, nchan, nstand, npol) ogulp_size = ntime * nchan * nstand * npol * 1 # ci4 self.oring.resize(ogulp_size) @@ -596,7 +597,7 @@ def main(self): # Quantization try: - Quantize(idata, qdata, scale=1. / np.sqrt(nchan)) + Quantize(idata, qdata, scale=1. / np.sqrt(nchan)) # TODO: check sqrt against main? I dont have this except NameError: qdata = bifrost.ndarray(shape=idata.shape, native=False, dtype="ci4") Quantize(idata, qdata, scale=1.0) @@ -963,6 +964,7 @@ def main(self): sampling_length = loc_data['sampling_length'].item() locs = loc_data['locs'][...] sll = loc_data['sll'].item() + locs = np.around(locs) except OSError: sampling_length, locs, sll = GenerateLocations( self.locations, @@ -974,6 +976,7 @@ def main(self): grid_resolution=self.grid_resolution, ) np.savez(locname, sampling_length=sampling_length, locs=locs, sll=sll) + locs=np.around(locs) try: copy_array(self.locs, bifrost.ndarray(locs.astype(np.int32))) except AttributeError: @@ -1014,12 +1017,24 @@ def main(self): raise ValueError( "Cannot write fits file without knowing polarization list" ) - ohdr_str = json.dumps(ohdr) + # Setup the kernels to include phasing terms for zenith # Phases are Ntime x Nchan x Nstand x Npol x extent x extent phasename = "phases_%s_%i_%i_%i_%i_%i_%i.npy" % (self.station.name, chan0, self.ntime_gulp, nchan, nstand, npol, self.ant_extent) phasename = os.path.join(self.cache_dir, phasename) + + # Zenith for LWA-SV + pce = np.deg2rad(88.0666283) + pca = np.deg2rad(90.9743763) + + # Create Phase Center Vector + phase_center_sv = np.array([np.cos(pce)*np.sin(pca), np.cos(pce)*np.cos(pca), np.sin(pce)]) + + ohdr["pc_elevation"] = pce + ohdr["pc_azimuth"] = pca + ohdr_str = json.dumps(ohdr) + try: phases = np.load(phasename) except OSError: @@ -1031,13 +1046,13 @@ def main(self): for i in range(nstand): # X a = self.station.antennas[2 * i + 0] - delay = a.cable.delay(freq) - a.stand.z / speed_of_light.value + delay = a.cable.delay(freq) - np.dot(phase_center_sv,[a.stand.x,a.stand.y, a.stand.z]) / speed_of_light.value phases[:, :, i, 0, :, :] = np.exp(2j * np.pi * freq * delay) phases[:, :, i, 0, :, :] /= np.sqrt(a.cable.gain(freq)) if npol == 2: # Y a = self.station.antennas[2 * i + 1] - delay = a.cable.delay(freq) - a.stand.z / speed_of_light.value + delay = a.cable.delay(freq) - np.dot(phase_center_sv,[a.stand.x,a.stand.y, a.stand.z]) / speed_of_light.value phases[:, :, i, 1, :, :] = np.exp(2j * np.pi * freq * delay) phases[:, :, i, 1, :, :] /= np.sqrt(a.cable.gain(freq)) # Explicit bad and suspect antenna masking - this will @@ -1104,6 +1119,7 @@ def main(self): time1 = time.time() udata = idata.copy(space="cuda") + if self.benchmark is True: time1a = time.time() print(" Input copy time: %f" % (time1a - time1)) @@ -1257,23 +1273,15 @@ def main(self): 1, nchan, npol ** 2, self.grid_size, self.grid_size ) try: - bf_romein_autocorr.execute(autocorrs_av, autocorr_g) + bf_vgrid_autocorr.execute(autocorrs_av, autocorr_g) except NameError: - bf_romein_autocorr = Romein() - bf_romein_autocorr.init( + bf_vgrid_autocorr = VGrid() + bf_vgrid_autocorr.init( autocorr_lo, gacphases, self.grid_size, polmajor=False ) - bf_romein_autocorr.execute(autocorrs_av, autocorr_g) - # try: - # bf_vgrid_autocorr.execute(autocorrs_av, autocorr_g) - # except NameError: - # bf_vgrid_autocorr = VGrid() - # bf_vgrid_autocorr.init( - # autocorr_lo, gacphases, self.grid_size, polmajor=False - # ) - # bf_vgrid_autocorr.execute(autocorrs_av, autocorr_g) + bf_vgrid_autocorr.execute(autocorrs_av, autocorr_g) autocorr_g = autocorr_g.reshape(1 * nchan * npol ** 2, self.grid_size, self.grid_size) - # autocorr_g = romein_float(autocorrs_av,autocorr_g,autocorr_il,autocorr_lx,autocorr_ly,autocorr_lz,self.ant_extent,self.grid_size,nstand,nchan*npol**2) + # Inverse FFT try: ac_fft.execute(autocorr_g, autocorr_g, inverse=True) @@ -1965,8 +1973,10 @@ def main(self): crit_pix_x = float(ihdr["grid_size_x"] / 2 + 1) # Need to correct for shift in center pixel when we flipped dec dimension # when writing npz, Only applies for even dimension size - crit_pix_y = float(ihdr["grid_size_y"] / 2 + 1) - (ihdr["grid_size_x"] + 1) % 2 - + crit_pix_y = float(ihdr["grid_size_y"] / 2 + 1) # - (ihdr["grid_size_x"] + 1) % 2 + + + delta_x = -dtheta_x * 180.0 / np.pi delta_y = dtheta_y * 180.0 / np.pi delta_f = ihdr["bw"] / ihdr["nchan"] @@ -1989,7 +1999,10 @@ def main(self): if nints >= self.ints_per_file: image = np.fft.fftshift(image, axes=(3, 4)) - image = image[:, :, :, ::-1, :] + # image = image[:, :, :, ::-1, :] + # TODO: contentious debate about whether the above line is correct OR the complete picture. + # maybe we need to flip x coords. + # maybe we need to offset ourselves by a pixel? # Restructure data in preparation to stuff into fits # Now (Ntimes, Npol, Nfreq, y, x) @@ -2016,12 +2029,13 @@ def main(self): ) time_array = t0 + np.arange(nints) * dt - lsts = time_array.sidereal_time("apparent") - coords = SkyCoord( - lsts.deg, ihdr["latitude"], obstime=time_array, unit="deg" - ).transform_to(FK5(equinox="J2000")) + # Set Observer + observinglocation = EarthLocation(lat = ihdr["latitude"]*u.deg , lon =ihdr["longitude"]*u.deg , height=1500*u.m) + aa = AltAz(location=observinglocation, obstime = time_array, az = ihdr["pc_azimuth"]*u.rad , alt = ihdr["pc_elevation"]*u.rad) + coords = SkyCoord(aa).transform_to(FK5(equinox="J2000")) + hdul = [] for im_num, d in enumerate(image): hdu = fits.ImageHDU(data=d) From b2a77ae8b8158b5f936436dd282aac68e2028961 Mon Sep 17 00:00:00 2001 From: Craig Taylor <73905895+ctaylor-physics@users.noreply.github.com> Date: Tue, 17 Oct 2023 13:48:03 -0600 Subject: [PATCH 32/32] Update LWA_EPIC.py small typo on import --- LWA_EPIC/LWA_EPIC.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/LWA_EPIC/LWA_EPIC.py b/LWA_EPIC/LWA_EPIC.py index e1369f1..9fcd6ae 100755 --- a/LWA_EPIC/LWA_EPIC.py +++ b/LWA_EPIC/LWA_EPIC.py @@ -22,7 +22,7 @@ from astropy.io import fits from astropy.time import Time, TimeDelta from astropy import units as u -from astropy.coordinates import SkyCoord, FK5 +from astropy.coordinates import SkyCoord, FK5, EarthLocation, AltAz from astropy.constants import c as speed_of_light import datetime