multi: make faster

demo/multi_resolve/beam_to_solid_angle.py

+import numpy as np
+
+
+def pix_to_beam(beam_size):
+    beam_solid_angle = np.pi * (beam_size / 2) ** 2
+    cell_size = beam_size / 5
+    pixel_solid_angle = cell_size ** 2
+    pix_to_beam = beam_solid_angle / pixel_solid_angle
+    return pix_to_beam

demo/multi_resolve/mosaic_correct.py

@@ -60,10 +60,11 @@ sdom = sky.target[3]
 x_direction = coords(sdom.shape[0], sdom.distances[0])
 y_direction = coords(sdom.shape[1], sdom.distances[1])
 sky_coordinates = np.array(np.meshgrid(
-    x_direction, y_direction, indexing='xy'))
+    # -x_direction, y_direction, indexing='xy'
+    x_direction, y_direction, indexing='xy'
+))
 
-
-output_directory = f"output/combined_corrected_{npix}"
+output_directory = f"output/combined_corrected_{npix}_old"
 
 
 sky_center = SkyCoord(center_ra, center_dec, frame=center_frame)
@@ -84,12 +85,7 @@ for fldid, oo in enumerate(all_obs):
     beam = ift.makeField(sdom, beam)
     # beam_operator = ift.DiagonalOperator(beam)
 
-    beam_direction = f'fld{fldid}'
-    beam_directions[beam_direction] = dict(
-        dx=dx,
-        dy=dy,
-        beam=beam
-    )
+    beam_directions[f'fld{fldid}'] = dict(dx=dx, dy=dy, beam=beam)
 
 
 # Used for the dtypes
@@ -131,10 +127,7 @@ def build_response(field_key, dx, dy, obs, sky_dtype):
         domain_dtype=obs.vis.dtype
     )
 
-    # FIXME: This response operator makes unnecessary multiple beam pattern
-    # multiplications. As the SKY_BEAMER calculates the beam pattern for all fields.
-    # Hence, we throw away all other fields.
-    return UPCAST_TO_STOKES_I @ RADIO_RESPONSE @ FIELD_EXTRACTOR @ SKY_BEAMER @ REDUCER
+    return UPCAST_TO_STOKES_I @ RADIO_RESPONSE @ FIELD_EXTRACTOR
 
 
 def build_energy(response, obs):
@@ -149,7 +142,8 @@ responses = []
 energies = []
 for kk, vv, obs in zip(beam_directions.keys(), beam_directions.values(), all_obs):
     R = build_response(kk, vv['dx'], vv['dy'], obs, tmp_sky.dtype)
-    responses.append(R)
+
+    responses.append(R @ SKY_BEAMER @ REDUCER)
     energies.append(build_energy(R, obs))
 
 
@@ -159,12 +153,16 @@ for ii in range(1):
 
     tot_response_adjoint = np.zeros_like(f.val)
     for rr, oo in zip(responses, all_obs):
-        tot_response_adjoint += rr.adjoint(oo.vis).val
+        dirty = rr.adjoint(oo.vis).val
+        tot_response_adjoint += dirty
+        plt.imshow(dirty[0, 0, 0], origin='lower')
+        plt.show()
+
     plt.imshow(tot_response_adjoint[0, 0, 0], origin='lower')
     plt.show()
 
 lh = reduce(lambda x, y: x+y, energies)
-lh = lh @ sky
+lh = lh @ SKY_BEAMER @ REDUCER @ sky
 
 
 def callback(samples, i):

demo/multi_resolve/ms2rve.py

@@ -19,7 +19,11 @@ outdir = args.output_dir
 makedirs(outdir, exist_ok=True)
 
 
-name = input.split('/')[-1] + '_fld{field:02d}'
-obs = rve.ms2observations(input, "DATA", True, 0)
+spw = 1
+
+name = input.split('/')[-1] + '_fld{field:02d}_spw{spw}'
+obs = rve.ms2observations(input, "DATA", True, spw)
 for ii, o in enumerate(obs):
-    o.save(join(outdir, name.format(field=ii)), False)
+    if o is None:
+        continue
+    o.save(join(outdir, name.format(field=ii, spw=spw)), False)

demo/multi_resolve/multi_obs.cfg

 [sky]
 freq mode = single
 polarization=I
-space npix x = 256
-space npix y = 256
+space npix x = 128
+space npix y = 128
 space fov x = 25as
 space fov y = 25as
 # space fov x = 30.72as