Updating example on MWC480 (polychromatic)

demos/example_image_reconstruction_betaLyrae.jl

@@ -1,5 +1,5 @@
 #
-# Very Basic Image reconstruction code
+# Reconstruction of beta Lyrae from snapshot MIRC-6T data
 #
 using OITOOLS
 oifitsfile = "./data/betlyr6t.oifits"

demos/example_image_reconstruction_polychromatic_yso.jl → demos/example_image_reconstruction_polychromatic_MWC480.jl

@@ -8,15 +8,14 @@ oifitsfile = "./data/MWC480.oifits"
 # Step 1: reconstruct a "gray" image
 data = readoifits(oifitsfile, filter_bad_data = true, polychromatic = false)[1,1]
 nx = 64 #number of pixels (side)
-pixsize = 0.2 # mas/pixel
+pixsize = 0.1 # mas/pixel
 ft = setup_nfft(data, nx, pixsize);
-regularizers = [ [ ["centering", 1e4], ["tv", 1e5] ]] 
-
+regularizers = [["centering", 1e7]] 
 pointsource = zeros(nx,nx); 
 pointsource[div(nx+1,2), div(nx+1,2)] = 1.0;
 x = copy(pointsource);
-x = reconstruct(x, data, ft, regularizers = regularizers, maxiter = 200, verb=false);
-
+x = reconstruct(x, data, ft, regularizers = regularizers, maxiter = 500, verb=true);
+imdisp(x.^.2, pixsize=pixsize)
 
 # Step 2: polychromatic reconstruction
 data = vec(readoifits(oifitsfile, filter_bad_data = true, polychromatic = true)) # vec is to get rid of degenerate (temporal) dimension

src/modelfit.jl

@@ -6,8 +6,8 @@
 #       - black body spectral law
 #       - spectral lines (Gaussian, Voigt)
 #       - custom r/μ profiles
-using Statistics, LinearAlgebra, Parameters, PyCall, UltraNest, LsqFit, NLopt, Printf, FFTW, NFFT
-
+using Statistics, LinearAlgebra, Parameters, PyCall, UltraNest, NLopt, Printf, FFTW, NFFT
+import LsqFit
 @with_kw mutable struct OIparam
     name::String = "" # optional name of the compoment (e.g. "primary", "central source")
     val::Float64 = 0
@@ -533,8 +533,10 @@ function fit_model_ultranest(data::OIdata, model::OImodel; lbounds = Float64[],
     end
 
     if cornerplot == true
-        PyDict(pyimport("matplotlib")."rcParams")["font.size"]=[10];
-        pyimport("ultranest.plot").cornerplot(result);
+        PyDict(pyimport("matplotlib")."rcParams")["font.size"]=[8];
+        histogram_color = "black"
+        contour_colors = py"{'colors': ['#0072B2','#56B4E9','#009E73','#F0E442'], 'linestyles': ['-', '-', '-', '-']}"
+        pyimport("ultranest.plot").cornerplot(result, contour_kwargs=contour_colors, color=histogram_color);
     end
     cvis_model = [];
     if calculate_vis == true
@@ -600,15 +602,15 @@ function fit_model_levenberg(data::OIdata, model::OImodel; verbose = true, calcu
     lbounds, hbounds = get_model_bounds(model);
     pinit = get_model_params(model);
     m = (x,p)->lsqmodelobs(p, model, data, weights=weights);
-    fit = curve_fit(m, [], ydata, wt, pinit, lower=lbounds, upper=hbounds, show_trace=verbose);
+    fit = LsqFit.curve_fit(m, [], ydata, wt, pinit, lower=lbounds, upper=hbounds, show_trace=verbose);
     minx = fit.param
     minf = model_to_chi2(data, model, minx, weights=weights);
 
     if fit.converged == true
         println("Levenberg-Marquardt fit converged to chi2 = $(minf) for p=$(minx)\n")
     end
-    sigma = stderror(fit)
-    covar = estimate_covar(fit)
+    sigma = LsqFit.stderror(fit)
+    covar = LsqFit.estimate_covar(fit)
     if verbose==true
         println("Name       \t\tMinimum\t\tMaximum\t\tInit\t\tConverged ± Error");
         pnames = get_model_pnames(model);

src/oichi2.jl

@@ -14,6 +14,22 @@ function setup_dft(data::OIdata, nx, pixsize)
     return setup_dft(data.uv, nx, pixsize)
 end
 
+# function setup_dft(data::Matrix{OIdata}, nx, pixsize)
+#     if size(data) == (1,1)
+#         return setup_dft(data[1,1].uv, nx, pixsize);
+#     else 
+#         error("Multidimensional DFT not implemented yet");
+#     end
+# end
+
+# function setup_nfft(data::Matrix{OIdata}, nx, pixsize)
+#     if size(data) == (1,1)
+#         return setup_nfft(data[1,1], nx, pixsize);
+#     else 
+#         return setup_nfft_multi(data, nx, pixsize);
+#     end
+# end
+
 function setup_nfft(data::OIdata, nx, pixsize)
     scale_rad = pixsize * (pi / 180.0) / 3600000.0*[-1;1].*data.uv;
     fftplan_uv  = plan_nfft(scale_rad, (nx,nx), m=4, σ=2.0);
@@ -37,22 +53,18 @@ function setup_nfft(uv::Array{Float64,2}, indx_vis, indx_v2, indx_t3_1, indx_t3_
     return [fftplan_uv,fftplan_vis,fftplan_v2,fftplan_t3_1,fftplan_t3_2,fftplan_t3_3]
 end
 
-
-function setup_nfft_t4(data, nx, pixsize)
-    scale_rad = pixsize * (pi / 180.0) / 3600000.0*[-1;1].*data.uv
-    fftplan_uv  = plan_nfft(scale_rad, (nx,nx), m=4, σ=2.0);
-    fftplan_vis  = plan_nfft(scale_rad[:, data.indx_vis], (nx,nx), m=4, σ=2.0);
-    fftplan_v2   = plan_nfft(scale_rad[:, data.indx_v2], (nx,nx), m=4, σ=2.0);
-    fftplan_t3_1 = plan_nfft(scale_rad[:, data.indx_t3_1], (nx,nx), m=4, σ=2.0);
-    fftplan_t3_2 = plan_nfft(scale_rad[:, data.indx_t3_2], (nx,nx), m=4, σ=2.0);
-    fftplan_t3_3 = plan_nfft(scale_rad[:, data.indx_t3_3], (nx,nx), m=4, σ=2.0);
-    fftplan_t4_1 = plan_nfft(scale_rad[:, data.indx_t4_1], (nx,nx), m=4, σ=2.0);
-    fftplan_t4_2 = plan_nfft(scale_rad[:, data.indx_t4_2], (nx,nx), m=4, σ=2.0);
-    fftplan_t4_3 = plan_nfft(scale_rad[:, data.indx_t4_3], (nx,nx), m=4, σ=2.0);
-    fftplan_t4_4 = plan_nfft(scale_rad[:, data.indx_t4_4], (nx,nx), m=4, σ=2.0);
-    return [fftplan_uv,fftplan_vis,fftplan_v2,fftplan_t3_1,fftplan_t3_2,fftplan_t3_3, fftplan_t4_1,fftplan_t4_2,fftplan_t4_3, fftplan_t4_4]
-end
-
+# function setup_nfft_multi(data, nx, pixsize)
+#     nwavs = size(data,1);
+#     nepochs = size(data,2);
+#     scale_rad = pixsize * (pi / 180.0) / 3600000.0;
+#     fftplan_multi = Array{Array{NFFT.NFFTPlan{Float64, 2, 1}},2}(undef, nwavs, nepochs);
+#     for i=1:nepochs
+#         for j=1:nwavs
+#             fftplan_multi[j,i]=setup_nfft(data[j,i], nx, pixsize);
+#         end
+#     end
+#     return fftplan_multi
+# end
 
 function setup_nfft_multiepochs(data, nx, pixsize)
     nepochs = size(data,1);