Work on plots

demo/post_mosaicing.py

@@ -80,28 +80,90 @@ def main():
     for ss in state:
         skysc.add(c.sky_domaintuple.force(ss))
 
-    sd_comparison(skysc.mean)
+    logskysc = ift.StatCalculator()
+    for ss in state:
+        logskysc.add(c.logsky.force(ss))
+
+    sc_sd_simdata = ift.StatCalculator()
+    op = c.responsesd.partial_insert(c.sky)
+    for ss in state:
+        sc_sd_simdata.add(op.force(ss))
+    sd_comparison(skysc.mean, logskysc.mean, sc_sd_simdata)
     dirty_images(skysc.mean)
 
 
-def sd_comparison(mean_sky):
+def sd_comparison(mean_sky, mean_logsky, sc_sd_simdata):
     fld, refval = rve.fits2field("../resources_multipointing/M83.spw17.I.sd.fits")
-    # FIXME Check with convert_alma_data.py script
     minfreq = 230172900000.0
     ind = int(np.round((minfreq - refval[0]) / fld.domain[0].distances[0]))
     extractor = ift.DomainTupleFieldInserter(fld.domain, 0, (ind,)).adjoint
     casasky = extractor(fld)
 
-    fig, axs = plt.subplots(1, 2)
-    axs = list(axs)
+    obs = c.sdobs["sd0"]
+    r_casa = (
+        rve.SingleDishResponse(
+            obs,
+            casasky.domain,
+            lambda x: rve.alma_beam_func(10.7, 0.75, np.mean(obs.freq), x),
+            c.global_phase_center,
+        )
+        @ ift.FieldAdapter(casasky.domain, "sky").adjoint
+    )
+
+    tmp = casasky.val_rw()
+    tmp[np.isnan(tmp)] = 0.0
+    tmp = ift.makeField(casasky.domain, tmp)
+
+    # FIXME Question: How can I compute data residuals for single dish data?
+    # FIXME Question: How do I arrive at the correct units?
+    d0 = obs.vis.val.ravel()
+    d1 = r_casa(tmp).val.ravel()
+    scaling_factor = np.vdot(d0, d1) / np.vdot(d1, d1)
+    print(
+        "Crazy scaling factor to make CASA reconstruction fit the data better",
+        scaling_factor,
+    )
+
+    tmp = tmp * scaling_factor
+    d1 = r_casa(tmp).val.ravel()
+
+    plt.hist(d0)
+    plt.savefig("data.png")
+    plt.close()
+
+    plt.hist(d1)
+    plt.savefig("data_sim.png")
+    plt.close()
+
+    fig, axs = plt.subplots(2, 2)
+    axs = list(axs.ravel())
+
+    axx = axs.pop(0)
+    imshow_symmetric(axx, r_casa.adjoint(obs.vis))
+    axx.set_title("Data")
 
     axx = axs.pop(0)
-    imshow_with_colorbar(axx, casasky)
+    imshow_symmetric(axx, r_casa.adjoint(obs.vis - r_casa(tmp)))
+    axx.set_title("CASA Residual")
+
+    plt.tight_layout()
+    plt.savefig("casa_sd_reconstruction.png")
+    plt.close()
+
+    fig, axs = plt.subplots(2, 2)
+    axs = list(axs.ravel())
+
+    axx = axs.pop(0)
+    imshow_with_colorbar(axx, casasky * scaling_factor, cmap="inferno")
     axx.set_title("CASA single dish")
 
     axx = axs.pop(0)
-    imshow_with_colorbar(axx, mean_sky)
-    axx.set_title("Resolve single dish posterior mean")
+    imshow_with_colorbar(axx, mean_sky, cmap="inferno")
+    axx.set_title("Resolve single dish")
+
+    axx = axs.pop(0)
+    imshow_with_colorbar(axx, mean_logsky.exp(), cmap="inferno")
+    axx.set_title("Resolve log. single dish")
 
     plt.tight_layout()
     fig.savefig("sd_reconstructions.png")