Add new smooth_templates.jl example (#35)

New `examples/templates/smooth_template.jl` added. This uses PyPlot rather than Plots.jl as I have been because I couldn't figure out how to get the color bars to scale how I wanted. This example required adding several new dependencies to the `docs/Project.toml`. 

This example is run on the docs build CI job and the plots from this example are included in `docs/src/fitting/fitting_intro.md`. This required adding steps to setup python and install matplotlib. I tried to use PyPlot.jl's default MiniConda installation but this failed sporadically so reverting to `setup-python` which hasn't failed yet. My plotting code uses tex rending so I also added a step to add those dependencies to the job. matplotlib v3.9 is currently not supported by PyPlot.jl so we are restricting matplotlib to 3.8.*.

.github/workflows/CI.yml

@@ -62,22 +62,32 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: '3.x'
+      - run: python -m pip install --upgrade pip
+      - run: pip install 'matplotlib==3.8.*'
       - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
       - uses: julia-actions/cache@v2
-      - run: |
+      - name: Instantiate docs project
+        run: |
           julia --project=docs -e '
             using Pkg
             Pkg.develop(PackageSpec(path=pwd()))
             Pkg.instantiate()'
-      - run: |
+      - name: Run doctests
+        run: |
           julia --project=docs -e '
             import Documenter: doctest, DocMeta
             import StarFormationHistories
             DocMeta.setdocmeta!(StarFormationHistories, :DocTestSetup, :(using StarFormationHistories); recursive=true)
             doctest(StarFormationHistories)'
-      - run: julia --project=docs docs/make.jl
+      - name: Install tex dependencies
+        run: sudo apt-get install dvipng texlive-latex-extra texlive-fonts-recommended cm-super
+      - name: Make and deploy docs
+        run: julia --project=docs docs/make.jl
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
           DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}

.gitignore

@@ -28,6 +28,8 @@ Manifest.toml
 
 # nbconvert generated static html in examples
 examples/*.html
+# rendered example figures
+examples/templates/*.svg
 
 # Some test files that don't need to be committed
 test/manual/madcash-2

docs/Project.toml

@@ -1,7 +1,19 @@
 [deps]
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+InitialMassFunctions = "e9251ff4-c148-4db3-bf46-89407750fae0"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 
 [compat]
+DelimitedFiles = "1"
+Distributions = "0.25"
 Documenter = "1"
+InitialMassFunctions = "0.1"
 Plots = "1"
+Printf = "<0.0.1, 1"
+PyPlot = "2"
+StatsBase = "0.34"

docs/make.jl

@@ -1,6 +1,18 @@
 using Documenter
 using StarFormationHistories
 
+
+# Run examples
+import PyPlot as plt
+plt.ioff()
+ENV["MPLBACKEND"] = "agg"
+# Run smooth_template.jl
+include("../examples/templates/smooth_template.jl")
+# Can't move yet as makedocs will clear the build folder.
+# Moving and showing in docs/src/fitting/fitting_intro.md.
+
+###########################################################
+
 # The `format` below makes it so that urls are set to "pretty" if you are pushing them to a hosting service, and basic if you are just using them locally to make browsing easier.
 
 # DocMeta.setdocmeta!(StarFormationHistories, :DocTestSetup, :(using StarFormationHistories; import Unitful; import UnitfulAstro); recursive=true)

docs/src/fitting/fitting_intro.md

@@ -40,6 +40,21 @@ Below we show a comparison of a smooth Hess diagram template constructed with [`
 
 ![Comparison of smooth Hess diagram template from `partial_cmd_smooth` and a Monte Carlo model made with `generate_stars_mass`.](figures/model_cmd.png)
 
+```@example
+mv("../../../examples/templates/template_compare.svg", "figures/template_compare.svg") # hide
+mv("../../../examples/templates/sigma_distribution.svg", "figures/sigma_distribution.svg") # hide
+nothing # hide
+```
+
+A worked example comparing a sampled stellar population with a smooth Hess diagram template is available in `examples/templates/smooth_template.jl`. The output figure is shown below. A distance modulus of 25 mag is used for this example, with photometric error and completeness functions roughly based on those we observe in the JWST/NIRCAM data of WLM (see [Weisz et al. 2024](https://ui.adsabs.harvard.edu/abs/2024ApJS..271...47W)). 
+
+![Comparison of CMD-sampled population with smooth Hess diagram template.](figures/template_compare.svg)
+
+At left is a population of stars sampled from an SSP with the methods described in the section of the documentation on [simulating CMDs](@ref simulate). The points from the isochrone are colored orange. The next figure shows the binned Hess diagram derived from these data. The next figure shows our smooth Hess diagram template calculated for this SSP. The final figure at right shows the residual between the data and model in units of standard deviations. These are sometimes called Pearson residuals. Below we show the distribution of these residuals, which should be Gaussian with mean 0 and standard deviation 1 if the model were perfect. We show a Gaussian PDF with standard deviation 1 and mean equal to the observed mean of the residuals for comparison.
+
+![Distribution of data - model residuals, in units of standard deviations.](figures/sigma_distribution.svg)
+
+The method used to create these smooth Hess diagram templates is [`partial_cmd_smooth`](@ref).
 
 ```@docs
 partial_cmd_smooth

examples/templates/smooth_template.jl

@@ -8,11 +8,20 @@ import PyPlot as plt
 import PyPlot: @L_str # For LatexStrings
 plt.rc("text", usetex=true)
 plt.rc("font", family="serif", serif=["Computer Modern"], size=14)
+# This gets close but not quite
+# plt.matplotlib.rcParams["axes.formatter.use_mathtext"] = true
+# plt.rc("font", family="serif", serif=["cmr10"], size=14)
 plt.rc("figure", figsize=(5,5))
-plt.rc("patch", linewidth=1, edgecolor="k", force_edgecolor=true) 
+plt.rc("patch", linewidth=1, edgecolor="k", force_edgecolor=true)
+# Disable interactive plotting when running on CI or building docs
+# if ("CI" in keys(ENV) && (ENV["CI"] == "true")) | (("DOCS_RUN" in keys(ENV)) && (ENV["DOCS_RUN"] == "true"))
+#     ENV["MPLBACKEND"] = "agg"
+#     plt.ioff()
+# end
 
 # Load example isochrone
-isochrone, mag_names = readdlm("../isochrone.txt", ' ', Float64, '\n'; header=true)
+# Path is relative to location of script, so use @__DIR__
+isochrone, mag_names = readdlm(joinpath(@__DIR__, "../isochrone.txt"), ' ', Float64, '\n'; header=true)
 # Unpack
 m_ini = isochrone[:,1]
 F090W = isochrone[:,2]
@@ -22,7 +31,7 @@ F150W = isochrone[:,3]
 distmod::Float64 = 25.0 # Distance modulus 
 
 # Set bins for Hess diagram
-edges = (range(-0.2, 1.2, length=75), range(distmod-6.0, distmod+3.0, length=75))
+edges = (range(-0.2, 1.2, length=75), range(distmod-6.0, distmod+5.0, length=75))
 
 # Set total stellar mass to normalize template to
 template_norm::Float64 = 1e7
@@ -77,55 +86,63 @@ signif[permutedims(obs_hess) .== 0] .= NaN
 # Plot
 fig,axs=plt.subplots(nrows=1,ncols=4,sharex=true,sharey=true,figsize=(20,5))
 fig.subplots_adjust(hspace=0.0,wspace=0.0)
-fig.suptitle(@sprintf("Stellar Mass: %.2e M\$_\\odot\$",template_norm))
+# fig.suptitle(@sprintf("Stellar Mass: %.2e M\$_\\odot\$",template_norm))
 
 axs[1].scatter(view(obs_mags,1,:) .- view(obs_mags,2,:), view(obs_mags,2,:), s=1, marker=".", c="k", alpha=0.05, rasterized=true, label="CMD-Sampled")
-axs[1].text(0.1,0.9,"Sampled CMD",transform=axs[1].transAxes)
+axs[1].text(0.05,0.95,@sprintf("Sampled CMD\nM\$_*\$ = %.2e M\$_\\odot\$", template_norm),transform=axs[1].transAxes,va="top",ha="left")
 
 im1 = axs[3].imshow(permutedims(template.weights), origin="lower", 
-      extent=(extrema(edges[1])..., extrema(edges[2])...), 
-      aspect="auto", cmap="Greys", norm=plt.matplotlib.colors.LogNorm(vmin=2.5 + log10(template_norm/1e7)), rasterized=true)
-axs[3].text(0.1,0.9,"Smooth Model",transform=axs[3].transAxes)
+                    extent=(extrema(edges[1])..., extrema(edges[2])...), 
+                    aspect="auto", cmap="Greys", norm=plt.matplotlib.colors.LogNorm(vmin=2.5 + log10(template_norm/1e7)), rasterized=true)
+axs[3].text(0.05,0.95,"Smooth Model",transform=axs[3].transAxes,va="top",ha="left")
 
 axs[2].imshow(permutedims(obs_hess), origin="lower", 
-    extent=(extrema(edges[1])..., extrema(edges[2])...), 
-    aspect="auto", cmap="Greys", norm=plt.matplotlib.colors.LogNorm(vmin=2.5 + log10(template_norm/1e7),vmax=im1.get_clim()[2]), rasterized=true, label="CMD-Sampled")
-axs[2].text(0.1,0.9,"Sampled Hess Diagram",transform=axs[2].transAxes)
+              extent=(extrema(edges[1])..., extrema(edges[2])...), 
+              aspect="auto", cmap="Greys", norm=plt.matplotlib.colors.LogNorm(vmin=2.5 + log10(template_norm/1e7),vmax=im1.get_clim()[2]), rasterized=true, label="CMD-Sampled")
+axs[2].text(0.05,0.95,"Sampled Hess Diagram",transform=axs[2].transAxes,va="top",ha="left")
 
-# im4 = axs[4].imshow( permutedims(obs_hess) .- permutedims(template.weights), origin="lower",
-#                      extent=(extrema(edges[1])..., extrema(edges[2])...), 
-#                      aspect="auto", cmap="Greys", clim=(-50,50), rasterized=true)
 im4 = axs[4].imshow( signif, 
                      origin="lower", extent=(extrema(edges[1])..., extrema(edges[2])...), 
                      aspect="auto", clim=(-2,2), rasterized=true)
-axs[4].text(0.1,0.9,L"(Obs - Model) / $\sigma$",transform=axs[4].transAxes)
-
-for ax in axs
-    ax.set_xlabel(L"F090W$-$F150W")
+axs[4].text(0.05,0.95,L"(Obs - Model) / $\sigma$",transform=axs[4].transAxes,va="top",ha="left")
+
+plot_isochrones::Bool = true
+for i in eachindex(axs)
+    axs[i].set_xlabel(L"F090W$-$F150W")
+    if plot_isochrones & (i != 4) # Don't plot on residual
+        axs[i].scatter(F090W .- F150W, F150W .+ distmod, marker=".", c="orange", s=1, alpha=0.3)
+    end
 end
 axs[1].set_ylabel("F150W")
 axs[1].set_ylim(reverse(extrema(edges[2]))) 
 axs[1].set_xlim(extrema(edges[1]))
 
 fig.colorbar(im1, ax=axs[1:3], pad=0.005, fraction=0.075) # fraction prevents too much padding on right
 fig.colorbar(im4, ax=axs[4], pad=0.015)
-# plt.savefig("test_cov.pdf", bbox_inches="tight")
-println( "sum(obs) - sum(template): ", sum(obs_hess) .- sum(template.weights))
-println( "Difference of sums / sum: ", (sum(obs_hess) .- sum(template.weights)) ./ sum(obs_hess))
-println( "Sum of residuals: ", sum(abs, permutedims(obs_hess) .- permutedims(template.weights)) )
+# println( "sum(obs) - sum(template): ", sum(obs_hess) .- sum(template.weights))
+# println( "Difference of sums / sum: ", (sum(obs_hess) .- sum(template.weights)) ./ sum(obs_hess))
+# println( "Sum of residuals: ", sum(abs, permutedims(obs_hess) .- permutedims(template.weights)) )
+plt.savefig(joinpath(@__DIR__,"template_compare.svg"), bbox_inches="tight")
 
 #################################
 # Distribution of σ discrepancies
 import Distributions: pdf, Normal, Poisson
-import StatsBase: mean
+import StatsBase: mean, std
 fig, ax1 = plt.subplots()
 hist1 = ax1.hist(filter(isfinite, signif), range=(-4,4), bins=25, density=true)
 ax1.set_xlim( extrema(hist1[2]) )
 ax1.set_xlabel("Residual / Standard Deviation")
 ax1.set_ylabel("PDF")
 let xplot = first(hist1[2]):0.01:last(hist1[2])
     # mean_σ = mean(filter(Base.Fix1(<,-5), filter(isfinite,signif)))
-    mean_σ = mean(filter(isfinite,signif))
-    ax1.plot( xplot, pdf.(Normal(mean_σ,1.0),xplot))
+    tmpgood = filter(x->first(hist1[2]) <= x <= last(hist1[2]), filter(isfinite,signif))
+    mean_σ = mean(tmpgood)
+    std_σ = std(tmpgood)
+    ax1.plot( xplot, pdf.(Normal(0.0,1.0),xplot), label="Ideal")
+    ax1.plot( xplot, pdf.(Normal(mean_σ,std_σ),xplot), label="Observed", c="magenta")
     ax1.axvline(mean_σ, c="k", ls="--")
+    ax1.text(0.05,0.95,(@sprintf("Mean: %.2f\n\$\\sigma\$: %.2f", mean_σ, std_σ)),transform=ax1.transAxes, va="top", ha="left")
 end
+ax1.legend(loc="upper right")
+plt.savefig(joinpath(@__DIR__,"sigma_distribution.svg"), bbox_inches="tight")
+