Merge pull request #383 from ptiede/ptiede-frequencysitearr

Multifrequency support for instrument modeling

Author: ptiede

Project.toml

@@ -89,7 +89,7 @@ ParameterHandling = "0.4, 0.5"
 Pigeons = "0.3, 0.4"
 PolarizedTypes = "0.1"
 PrettyTables = "1, 2"
-Pyehtim = "0.1"
+Pyehtim = "0.2"
 RecipesBase = "1"
 Reexport = "1"
 SpecialFunctions = "0.10, 1, 2"

docs/CondaPkg.toml

@@ -26,10 +26,14 @@ withenv("JULIA_DEBUG"=>"Literate") do
         jl_expr = "using Literate;"*
                   "preprocess(path, str) = replace(str, \"__DIR = @__DIR__\" => \"__DIR = \\\"\$(dirname(path))\\\"\");"*
                   "Literate.markdown(\"$(p_)\", \"$(joinpath(OUTPUT, d))\";"*
-                  "name=\"$name\", execute=false, flavor=Literate.DocumenterFlavor(),"*
+                  "name=\"$name\", execute=true, flavor=Literate.DocumenterFlavor(),"*
                   "preprocess=Base.Fix1(preprocess, \"$(p_)\"))"
         cm = `julia --project=$(@__DIR__) -e $(jl_expr)`
-        run(cm)
+        try
+            run(cm)
+        catch e
+            @warn "there was an issues with $cm\n $e"
+        end
 
     end
 end

docs/tutorials.jl

@@ -18,7 +18,7 @@ CairoMakie = "0.12"
 Distributions = "0.25"
 Optimization = "4"
 Plots = "1"
-Pyehtim = "0.1"
+Pyehtim = "0.2"
 StableRNGs = "1"
 StatsBase = "0.34"
 VLBIImagePriors = "0.9"

examples/advanced/HybridImaging/Project.toml

@@ -1,3 +1,14 @@
+import Pkg #hide
+__DIR = @__DIR__ #hide
+pkg_io = open(joinpath(__DIR, "pkg.log"), "w") #hide
+Pkg.activate(__DIR; io=pkg_io) #hide
+Pkg.develop(; path=joinpath(__DIR, "..", "..", ".."), io=pkg_io) #hide
+Pkg.instantiate(; io=pkg_io) #hide
+Pkg.precompile(; io=pkg_io) #hide
+close(pkg_io) #hide
+
+
+
 # # Hybrid Imaging of a Black Hole
 
 # In this tutorial, we will use **hybrid imaging** to analyze the 2017 EHT data.
@@ -16,18 +27,6 @@
 # This is the approach we will take in this tutorial to analyze the April 6 2017 EHT data
 # of M87.
 
-import Pkg #hide
-__DIR = @__DIR__ #hide
-pkg_io = open(joinpath(__DIR, "pkg.log"), "w") #hide
-Pkg.activate(__DIR; io=pkg_io) #hide
-Pkg.develop(; path=joinpath(__DIR, "..", "..", ".."), io=pkg_io) #hide
-Pkg.instantiate(; io=pkg_io) #hide
-Pkg.precompile(; io=pkg_io) #hide
-close(pkg_io) #hide
-
-ENV["GKSwstype"] = "nul" #hide
-
-
 # ## Loading the Data
 
 # To get started we will load Comrade
@@ -154,12 +153,10 @@ skym = SkyModel(sky, skyprior, g; metadata=skymetadata)
 using Enzyme
 post = VLBIPosterior(skym, intmodel, dvis; admode=set_runtime_activity(Enzyme.Reverse))
 
-# To sample from our prior we can do
-xrand = prior_sample(rng, post)
-
-# and then plot the results
+# We can sample from the prior to see what the model looks like
 using DisplayAs #hide
 import CairoMakie as CM
+xrand = prior_sample(rng, post)
 CM.activate!(type = "png", px_per_unit=1) #hide
 gpl = imagepixels(μas2rad(200.0), μas2rad(200.0), 128, 128)
 fig = imageviz(intensitymap(skymodel(post, xrand), gpl));
@@ -176,7 +173,7 @@ fig |> DisplayAs.PNG |> DisplayAs.Text #hide
 using Optimization
 using OptimizationOptimJL
 xopt, sol = comrade_opt(post, LBFGS(); 
-                        initial_params=xrand, maxiters=2000, g_tol=1e0)
+                        initial_params=xrand, maxiters=2000, g_tol=1e0);
 
 
 # First we will evaluate our fit by plotting the residuals
@@ -187,7 +184,8 @@ fig |> DisplayAs.PNG |> DisplayAs.Text #hide
 # These residuals suggest that we are substantially overfitting the data. This is a common
 # side effect of MAP imaging. As a result if we plot the image we see that there
 # is substantial high-frequency structure in the image that isn't supported by the data.
-imageviz(intensitymap(skymodel(post, xopt), gpl), figure=(;resolution=(500, 400),))
+fig = imageviz(intensitymap(skymodel(post, xopt), gpl), figure=(;resolution=(500, 400),));
+fig |> DisplayAs.PNG |> DisplayAs.Text #hide
 
 
 # To improve our results we will now move to Posterior sampling. This is the main method

examples/advanced/HybridImaging/main.jl

@@ -18,6 +18,6 @@ Distributions = "0.25"
 Optimization = "4"
 Pigeons = "0.4"
 Plots = "1"
-Pyehtim = "0.1"
+Pyehtim = "0.2"
 StableRNGs = "1"
 VLBIImagePriors = "0.9"

examples/beginner/GeometricModeling/Project.toml

@@ -1,3 +1,13 @@
+import Pkg; #hide
+__DIR = @__DIR__ #hide
+pkg_io = open(joinpath(__DIR, "pkg.log"), "w") #hide
+Pkg.activate(__DIR; io=pkg_io) #hide
+Pkg.develop(; path=joinpath(__DIR, "..", "..", ".."), io=pkg_io) #hide
+Pkg.instantiate(; io=pkg_io) #hide
+Pkg.precompile(; io=pkg_io) #hide
+close(pkg_io) #hide
+
+
 # # Geometric Modeling of EHT Data
 
 # `Comrade` has been designed to work with the EHT and ngEHT.
@@ -9,22 +19,15 @@
 # In this tutorial, we will construct a similar model and fit it to the data in under
 # 50 lines of code (sans comments). To start, we load Comrade and some other packages we need.
 
-import Pkg #hide
-__DIR = @__DIR__ #hide
-pkg_io = open(joinpath(__DIR, "pkg.log"), "w") #hide
-Pkg.activate(__DIR; io=pkg_io) #hide
-Pkg.develop(; path=joinpath(__DIR, "..", "..", ".."), io=pkg_io) #hide
-Pkg.instantiate(; io=pkg_io) #hide
-Pkg.precompile(; io=pkg_io) #hide
-close(pkg_io) #hide
-ENV["GKSwstype"] = "nul" #hide
-
+# To get started we load Comrade.
+#-
 
 using Comrade
 
 
+# Currently we use eht-imaging for data management, however this will soon be replaced
+# by a pure Julia solution.
 using Pyehtim
-
 # For reproducibility we use a stable random number genreator
 using StableRNGs
 rng = StableRNG(42)
@@ -81,7 +84,7 @@ prior = (
           ξτ= Uniform(0.0, π),
           f = Uniform(0.0, 1.0),
           σG = Uniform(μas2rad(1.0), μas2rad(100.0)),
-          τG = Uniform(0.0, 1.0),
+          τG = Exponential(1.0),
           ξG = Uniform(0.0, 1π),
           xG = Uniform(-μas2rad(80.0), μas2rad(80.0)),
           yG = Uniform(-μas2rad(80.0), μas2rad(80.0))
@@ -146,8 +149,8 @@ fpost = asflat(post)
 p = prior_sample(rng, post)
 
 # and then transform it to transformed space using T
-logdensityof(cpost, Comrade.TV.inverse(cpost, p))
-logdensityof(fpost, Comrade.TV.inverse(fpost, p))
+logdensityof(cpost, Comrade.inverse(cpost, p))
+logdensityof(fpost, Comrade.inverse(fpost, p))
 
 # note that the log densit is not the same since the transformation has causes a jacobian to ensure volume is preserved.
 
@@ -238,3 +241,4 @@ DisplayAs.Text(DisplayAs.PNG(p))
 # and the model, divided by the data's error:
 p = residual(post, chain[end]);
 DisplayAs.Text(DisplayAs.PNG(p))
+

examples/beginner/GeometricModeling/main.jl

@@ -6,4 +6,4 @@ Pyehtim = "3d61700d-6e5b-419a-8e22-9c066cf00468"
 
 [compat]
 Plots = "1"
-Pyehtim = "0.1"
+Pyehtim = "0.2"

examples/beginner/LoadingData/Project.toml

@@ -1,14 +1,3 @@
-# # Loading Data into Comrade
-
-# The VLBI field does not have a standardized data format, and the EHT uses a
-# particular uvfits format similar to the optical interferometry oifits format.
-# As a result, we reuse the excellent `eht-imaging` package to load data into `Comrade`.
-
-# Once the data is loaded, we then convert the data into the tabular format `Comrade`
-# expects. Note that this may change to a Julia package as the Julia radio
-# astronomy group grows.
-
-# To get started, we will load `Comrade` and `Plots` to enable visualizations of the data
 import Pkg #hide
 __DIR = @__DIR__ #hide
 pkg_io = open(joinpath(__DIR, "pkg.log"), "w") #hide
@@ -18,11 +7,18 @@ Pkg.instantiate(; io=pkg_io) #hide
 Pkg.precompile(; io=pkg_io) #hide
 close(pkg_io) #hide
 
-ENV["GKSwstype"] = "nul" #hide
+# # Loading Data into Comrade
 
+# The VLBI field does not have a standardized data format, and the EHT uses a
+# particular uvfits format similar to the optical interferometry oifits format.
+# As a result, we reuse the excellent `eht-imaging` package to load data into `Comrade`.
 
-using Comrade
+# Once the data is loaded, we then convert the data into the tabular format `Comrade`
+# expects. Note that this may change to a Julia package as the Julia radio
+# astronomy group grows.
 
+# To get started, we will load `Comrade` and `Plots` to enable visualizations of the data
+using Comrade
 using Plots
 
 # We also load Pyehtim since it loads eht-imaging into Julia using PythonCall and exports

examples/beginner/LoadingData/main.jl

@@ -24,6 +24,6 @@ Distributions = "0.25"
 Optimization = "4"
 Pkg = "1"
 Plots = "1"
-Pyehtim = "0.1"
+Pyehtim = "0.2"
 StableRNGs = "1"
 VLBIImagePriors = "0.9"