From a44ed76e3c92e6dce3d482466410a84cd6b0a391 Mon Sep 17 00:00:00 2001 From: Reinhold Willcox <33348824+reinhold-willcox@users.noreply.github.com> Date: Tue, 8 Oct 2024 23:52:53 +0200 Subject: [PATCH] More angles in output (#76) * updated version * added sum of omega_i and nu_i , and sum of omega_f and nu_f, as output parameters, which are better constrained for 0 eccentricity systems * fixed function calls to output \nu_f * added in extra nu_f parameters that break the code * removed problematic nu_f references * don't include revise in example dir * removed nu_f from runtests * fixed example plotting * fixed Project.toml for docs --- docs/Project.toml | 1 + examples/Project.toml | 1 + examples/plot_results_vfts243.jl | 69 ++++++++++++++++++-------------- src/CornerPlots.jl | 36 +++-------------- src/TuringModels.jl | 5 ++- 5 files changed, 48 insertions(+), 64 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index f85b128..af9c9d0 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,5 +1,6 @@ [deps] CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +CornerPlotting = "4a808bcd-0da0-4a1c-a768-2071018375da" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" diff --git a/examples/Project.toml b/examples/Project.toml index c37744a..f6de805 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -2,6 +2,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +CornerPlotting = "4a808bcd-0da0-4a1c-a768-2071018375da" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Revise = "295af30f-e4ad-537b-8983-00126c2a3abe" diff --git a/examples/plot_results_vfts243.jl b/examples/plot_results_vfts243.jl index 647f018..6ddc1fa 100644 --- a/examples/plot_results_vfts243.jl +++ b/examples/plot_results_vfts243.jl @@ -7,6 +7,7 @@ the previous example. We start by loading up using CairoMakie using SideKicks +using CornerPlotting using Distributions results, observations, priors, metadata = SideKicks.ExtractResults("vfts243_results.hdf5") @@ -17,23 +18,26 @@ As an initial check, we can plot the quantities that were used as observations. consistency check is to verify these are consistent with the MCMC samples. =# -plotting_props_obs_check = SideKicks.createPlottingProps([ - [:m1_f, m_sun, [15,40], L"M_1\;[M_{\odot}]"], - [:P_f, day, [10.35,10.45], L"P_f\;[\mathrm{days}]"], - [:e_f, 1, [0,0.1], L"e_f"], - [:K1, km_per_s, [77,90], L"K_1 \;[\mathrm{km s}^{-1}]"], -]) - -cp = create_corner_plot(results, plotting_props_obs_check, - supertitle="VFTS 243 - observables", - dists_to_plot = Dict( - :m1_f => Normal(25.0,2.3), - :e_f => Normal(0.017,0.012), - :K1 => Normal(81.4, 1.3), - :P_f => Normal(10.4031, 0.01) - ), - nbins = 20, nbins_contour = 10 - ) +names = [:m1_f, :P_f, :e_f, :K1] +scaling = Dict( :m1_f => m_sun, + :P_f => day, + :e_f => 1, + :K1 => km_per_s + ) +labels = Dict( :m1_f => L"M_{1f}\;[M_{\odot}]", + :P_f => L"P_f\;[\mathrm{days}]", + :e_f => L"e_f", + :K1 => L"K_1 \;[\mathrm{km s}^{-1}]", + ) + +set_theme!(CornerPlotting.default_theme()) +cp = CornerPlotting.CornerPlot(results, names, scaling=scaling, labels=labels) + +CornerPlotting.plot_extra_1D_distribution(cp, :m1_f, Normal(25.0,2.3)) +CornerPlotting.plot_extra_1D_distribution(cp, :e_f, Normal(0.017,0.012)) +CornerPlotting.plot_extra_1D_distribution(cp, :K1, Normal(81.4, 1.3)) +CornerPlotting.plot_extra_1D_distribution(cp, :P_f, Normal(10.4031, 0.01)) + save("vfts243_observables.png", cp.fig) cp.fig @@ -44,20 +48,23 @@ And, after verifying the samples do correspond to our observational constraints, consequences for explosion itself. =# -plotting_props = SideKicks.createPlottingProps([ - [:m2, m_sun, [0,25], L"M_2 \;[M_{\odot}]"], - [:dm2, m_sun, [0, 6], L"ΔM_2 \;[M_{\odot}]"], - [:P, day, [8,12], L"P \;[\mathrm{days}]"], - [:vkick, km_per_s, [0,50], L"v_{kick} \;[\mathrm{km s}^{-1}]"], - [:vsys, km_per_s, [0,50], L"v_{\mathrm{sys}} \;[\mathrm{km s}^{-1}]"], -]) - -cp = create_corner_plot(results, plotting_props, - show_CIs=true, - supertitle="VFTS 243 - derived quantities", - fraction_1D = 0.9, - nbins = 20, nbins_contour = 10 - ) +names = [:m2_i, :dm2, :P_i, :vkick, :vsys] +scaling = Dict( :m2_i => m_sun, + :dm2 => m_sun, + :P_i => day, + :vkick => km_per_s, + :vsys => km_per_s, + ) +labels = Dict( :m2_i => L"M_{2f}\;[M_{\odot}]", + :dm2 => L"ΔM_2 \;[M_{\odot}]", + :P_i => L"P_i \;[\mathrm{days}]", + :vkick => L"v_{kick} \;[\mathrm{km s}^{-1}]", + :vsys => L"v_{\mathrm{sys}} \;[\mathrm{km s}^{-1}]", + ) + + +set_theme!(CornerPlotting.default_theme()) +cp = CornerPlotting.CornerPlot(results, names, scaling=scaling, labels=labels) save("vfts243_derived.png", cp.fig) diff --git a/src/CornerPlots.jl b/src/CornerPlots.jl index 00de92c..a079c70 100644 --- a/src/CornerPlots.jl +++ b/src/CornerPlots.jl @@ -11,43 +11,17 @@ PlottingProps contains the props, units, ranges, and names (latex encouraged) for each property to plot. """ -@kwdef mutable struct PlottingProps - props::Vector{Symbol} - units::Vector{Float64} - ranges::Vector{Any} # RTW TODO - names::Vector{AbstractString} -end # TODO: need to make this return names, labels, ranges, scaling function createPlottingProps(props_matrix::Vector{Vector{Any}}) props_matrix = stack(props_matrix) # make into actual matrix - names = props_matrix[1,:] - #scaling = Dict(zip(names, props_matrix[2,:])) - #ranges = Dict(zip(names, props_matrix[3,:])) - #labels = Dict(zip(names, props_matrix[4,:])) - #return [names, scaling, ranges, labels] - return PlottingProps( - props = props_matrix[1,:], - units = props_matrix[2,:], - ranges = props_matrix[3,:], - names = props_matrix[4,:]) + names = convert(Vector{Symbol}, props_matrix[1,:]) + scaling = Dict(zip(names, props_matrix[2,:])) + ranges = Dict(zip(names, props_matrix[3,:])) + labels = Dict(zip(names, props_matrix[4,:])) + return [names, scaling, ranges, labels] end -function addPlottingProp(props_matrix::PlottingProps, new_prop::Vector{Any}) - props = props_matrix.props - units = props_matrix.units - ranges = props_matrix.ranges - names = props_matrix.names - push!(props, new_prop[1]) - push!(units , new_prop[2]) - push!(ranges, new_prop[3]) - push!(names, new_prop[4]) - return PlottingProps( - props = props, - units = units, - ranges = ranges, - names = names) -end """ create_corner_plot(results, plotting_props; diff --git a/src/TuringModels.jl b/src/TuringModels.jl index 5e9522e..26dd92b 100644 --- a/src/TuringModels.jl +++ b/src/TuringModels.jl @@ -360,14 +360,15 @@ function create_general_mcmc_model(; end dm2 = m2_i - m2_f + sum_ωi_νi = ω_i + ν_i - return ( m1_i, m2_i, P_i, e_i, a_i, Ω_i, ω_i, i_i, + return ( m1_i, m2_i, P_i, e_i, a_i, Ω_i, ω_i, i_i, sum_ωi_νi, vkick, θ, ϕ, dm2, frac, ν_i, m1_f, m2_f, P_f, e_f, a_f, Ω_f, ω_f, i_f, K1, K2, v_N, v_E, v_r, vsys, venv_N, venv_E, venv_r, vsys_N, vsys_E, vsys_r ) end - return_props = [:m1_i, :m2_i, :P_i, :e_i, :a_i, :Ω_i, :ω_i, :i_i, + return_props = [:m1_i, :m2_i, :P_i, :e_i, :a_i, :Ω_i, :ω_i, :i_i, :sum_ωi_νi, :vkick, :θ, :ϕ, :dm2, :frac, :ν_i, :m1_f, :m2_f, :P_f, :e_f, :a_f, :Ω_f, :ω_f, :i_f, :K1, :K2, :v_N, :v_E, :v_r, :vsys,