Notes

demo/imaging_with_automatic_weighting.py

@@ -13,6 +13,8 @@ from os.path import isfile, splitext
 
 def main():
     parser = argparse.ArgumentParser()
+    # TODO Add modes: reconstruct, plot
+    # TODO OU process for automatic weighting
     parser.add_argument("-j", type=int, default=1)
     parser.add_argument("--use-cached", action="store_true")
     parser.add_argument("--use-wgridding", action="store_true")
@@ -40,16 +42,16 @@ def main():
     fov = np.array([rve.str2rad(args.xfov), rve.str2rad(args.yfov)])
     npix = np.array([args.xpix, args.ypix])
     rve.set_epsilon(1 / 10 / obs.max_snr())
-    ppos = []
-    for point in args.point:
-        ppos.append([rve.str2rad(point[0]), rve.str2rad(point[1])])
 
     dom = ift.RGSpace(npix, fov / npix)
     logsky = ift.SimpleCorrelatedField(
         dom, args.diffusefluxlevel, (1, 0.1), (5, 1), (1.2, 0.4), (0.2, 0.2), (-2, 0.5)
     )
     diffuse = logsky.exp()
-    if len(ppos) > 0:
+    if args.point is not None:
+        ppos = []
+        for point in args.point:
+            ppos.append([rve.str2rad(point[0]), rve.str2rad(point[1])])
         inserter = rve.PointInserter(dom, ppos)
         points = ift.InverseGammaOperator(
             inserter.domain, alpha=0.5, q=0.2 / dom.scalar_dvol
@@ -58,6 +60,7 @@ def main():
         sky = diffuse + points
     else:
         sky = diffuse
+    # TODO Add mode with independent noise learning
     npix = 2500
     effuv = np.linalg.norm(obs.effective_uv().T, axis=1)
     assert obs.nfreq == obs.npol == 1
@@ -77,6 +80,7 @@ def main():
     plotter.add("power spectrum bayesian weighting", logwgt.power_spectrum)
 
     if rve.mpi.master:
+        if args.point is not None:
             # MAP points with original weights
             lh = rve.ImagingLikelihood(obs, points)
             ham = ift.StandardHamiltonian(lh)
@@ -97,10 +101,13 @@ def main():
             "normalized residuals (original weights)", lh.normalized_residual
         )
         ham = ift.StandardHamiltonian(lh)
-        state = rve.MinimizationState(
-            ift.MultiField.union([0.1 * ift.from_random(diffuse.domain), state.mean]),
-            [],
+        if args.point is None:
+            fld = 0.1 * ift.from_random(diffuse.domain)
+        else:
+            fld = ift.MultiField.union(
+                [0.1 * ift.from_random(diffuse.domain), state.mean]
             )
+        state = rve.MinimizationState(fld, [])
         mini = ift.NewtonCG(
             ift.GradientNormController(name="newton", iteration_limit=20)
         )