Skip to content

Sheshuk/fix model notebooks #402

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
JostMigenda merged 14 commits into from
Oct 7, 2025
Merged

Sheshuk/fix model notebooks #402

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
2bcc90a
Adding expand_dims, so Container would work with sdqueezed data arrays
Sheshuk Jul 4, 2025
effde21
Updated plot method
Sheshuk Jul 3, 2025
fc243c7
Skip integration over scalar dimension
Sheshuk Jul 4, 2025
3c9fdca
Updated the notebooks to use Flux.plot
Sheshuk Jul 4, 2025
9deea4d
Adding nbtest to run notebooks as tests separtely from doc building
Sheshuk Jul 7, 2025
b17bd24
Fix typo
Sheshuk Jul 7, 2025
78fa753
Allow errors in the notebooks, when building docs
Sheshuk Jul 7, 2025
238049a
Remove seaborn style
Sheshuk Jul 7, 2025
6b512eb
Fix Fornax: store additional parameters as attributes instead of pass…
Sheshuk Jul 9, 2025
c148c65
Fix Fornax_2021
Sheshuk Jul 9, 2025
3f03a48
Fix fornax_2022
Sheshuk Jul 9, 2025
2793e52
Update FlavorTransform notebook from Jpkneller
Sheshuk Oct 7, 2025
2f07e69
Remove nbsphinx_allow_errors
Sheshuk Oct 7, 2025
6fd113b
Update the plots
Sheshuk Oct 7, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 3 additions & 0 deletions .github/workflows/integration.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,9 @@ jobs:
run: |
python -m unittest python/snewpy/test/simplerate_integrationtest.py
pytest -m 'snowglobes'
- name: Run notebooks test
run: |
pytest --nbval-lax doc/source/nb
- uses: r-lib/actions/setup-pandoc@v2

- name: Build HTML docs
Expand Down
1 change: 1 addition & 0 deletions doc/source/conf.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,7 @@
# directories to ignore when looking for source files.
exclude_patterns = ['nb/dev']

nbsphinx_allow_errors = False
# The reST default role (used for this markup: `text`) to use for all
# documents.
#default_role = None
Expand Down
200 changes: 80 additions & 120 deletions doc/source/nb/FlavorTransformation.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,20 +19,19 @@
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"\n",
"from astropy import units as u \n",
"from snewpy.neutrino import MassHierarchy, MixingParameters, ThreeFlavorMixingParameters, FourFlavorMixingParameters\n",
"import snewpy.flavor_transformation as xforms\n",
"from snewpy.flavor import ThreeFlavor, FourFlavor\n",
"\n",
"from snewpy.neutrino import Flavor, MassHierarchy, MixingParameters\n",
"from snewpy.models.ccsn import Nakazato_2013\n",
"from snewpy.flavor_transformation import AdiabaticMSW, NonAdiabaticMSWH, \\\n",
" TwoFlavorDecoherence, ThreeFlavorDecoherence, \\\n",
" NeutrinoDecay, AdiabaticMSWes, NonAdiabaticMSWes\n",
"\n",
"mpl.rc('font', size=18)\n",
"%matplotlib inline\n",
"from astropy import units as u\n",
"from astropy import constants as c\n",
"from astropy.coordinates import SkyCoord, EarthLocation, AltAz\n",
"from astropy.time import Time \n",
"\n",
"# `np.trapz` was renamed in numpy 2.0 and now raises a DeprecationWarning\n",
"if np.__version__ < \"2.0\":\n",
" np.trapezoid = np.trapz"
"mpl.rc('font', size=18)\n",
"%matplotlib inline"
]
},
{
Expand Down Expand Up @@ -71,65 +70,47 @@
" A FlavorTransformation subclass; used to create an instance.\n",
" \"\"\"\n",
" \n",
" energies = np.linspace(0,60,121) * u.MeV \n",
" d = (10*u.kpc).to('cm').value # distance to SN\n",
" energies = np.linspace(1,60,119) * u.MeV\n",
" d = 10*u.kpc # distance to SN\n",
" \n",
" times = model.get_time()\n",
" burst_epoch = times <= 0.1*u.s\n",
" accretion_epoch = (times > 0.1*u.s) & (times <= 0.5*u.s)\n",
" cooling_epoch = (times > 0.5*u.s) & (times <= 10*u.s)\n",
" \n",
" ilum = {}\n",
" olum_nmo = {}\n",
" olum_imo = {}\n",
" ispec = model.get_flux(times, energies, d)\n",
" ospec_nmo = model.get_flux(times, energies, d, xform_nmo)\n",
" ospec_imo = model.get_flux(times, energies, d, xform_imo)\n",
" \n",
" ilum = ispec.integrate_or_sum('energy') \n",
" olum_nmo = ospec_nmo.integrate_or_sum('energy') \n",
" olum_imo = ospec_imo.integrate_or_sum('energy') \n",
" \n",
" for flavor in Flavor:\n",
" ilum[flavor] = np.zeros(len(times))\n",
" olum_nmo[flavor] = np.zeros(len(times))\n",
" olum_imo[flavor] = np.zeros(len(times))\n",
"\n",
" # Compute the transformed and untransformed flux at each time.\n",
" for i, t in enumerate(times):\n",
" ispec = model.get_initial_spectra(t, energies)\n",
" ospec_nmo = model.get_transformed_spectra(t, energies, xform_nmo)\n",
" ospec_imo = model.get_transformed_spectra(t, energies, xform_imo)\n",
"\n",
" for flavor in Flavor:\n",
" for j, E in enumerate(energies):\n",
" ispec[flavor][j] /= (4.*np.pi*d**2)\n",
" ospec_nmo[flavor][j] /= (4.*np.pi*d**2)\n",
" ospec_imo[flavor][j] /= (4.*np.pi*d**2)\n",
" \n",
" for flavor in Flavor:\n",
" ilum[flavor][i] = np.trapezoid(ispec[flavor].to('1/(erg*s)'), energies.to('erg')).value \n",
" olum_nmo[flavor][i] = np.trapezoid(ospec_nmo[flavor].to('1/(erg*s)'), energies.to('erg')).value\n",
" olum_imo[flavor][i] = np.trapezoid(ospec_imo[flavor].to('1/(erg*s)'), energies.to('erg')).value \n",
" \n",
" # make the figures \n",
" fig, axes = plt.subplots(3,3, figsize=(20,12), tight_layout=True)\n",
" \n",
" smax = [0.,0.,0.]\n",
" titles = ['Untransformed', 'Transformed (NMO)', 'Transformed (IMO)']\n",
" for i, spec in enumerate([ilum, olum_nmo, olum_imo]):\n",
" for i, spec in enumerate([ilum, olum_nmo, olum_imo]): \n",
" for j, phase in enumerate([burst_epoch, accretion_epoch, cooling_epoch]):\n",
" ax = axes[i,j]\n",
" timeunits = 'ms' if j==0 else 's'\n",
" \n",
" for flavor in Flavor:\n",
" x_values = times[phase].to(timeunits)\n",
" for flavor in ThreeFlavor:\n",
" y_values = spec.array.squeeze()[flavor][phase].to('1/(cm2 s)')\n",
" if i == 0:\n",
" smax[j] = np.maximum(smax[j], 1.1*np.max(spec[flavor][phase]))\n",
" smax[j] = np.maximum(smax[j], 1.1*np.max(y_values))\n",
" \n",
" ax.plot(times[phase].to(timeunits),\n",
" spec[flavor][phase], label=flavor.to_tex(), lw=3,\n",
" color='C0' if flavor.is_electron else 'C1',\n",
" ax.plot(x_values, y_values, label=flavor.to_tex(), lw=3,\n",
" color={'E':'C0','MU':'C1','TAU':'C2'}[flavor.lepton],\n",
" ls='-' if flavor.is_neutrino else ':')\n",
" \n",
" ax.set(xlim=(times[phase][0].to(timeunits).value, times[phase][-1].to(timeunits).value),\n",
" ylim=(0, smax[j]))\n",
" ax.set(xlim=(x_values[0].value, x_values[-1].value),\n",
" ylim=(0, smax[j].value))\n",
" \n",
" if j==0:\n",
" ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')\n",
" ax.legend(loc='upper right', ncol=1, fontsize=18)\n",
" ax.set(ylabel=f'flux [{y_values.unit._repr_latex_()}]')\n",
" ax.legend(loc='upper right', ncol=1, fontsize=18)\n",
" if j==1:\n",
" ax.set(title=titles[i])\n",
" if i < 2:\n",
Expand Down Expand Up @@ -171,72 +152,39 @@
" Time to compute the spectrum.\n",
" \"\"\"\n",
"\n",
" energies = np.linspace(0,60,121) * u.MeV \n",
" d = (10*u.kpc).to('cm').value # distance to SN\n",
" energies = np.linspace(1,60,119) * u.MeV\n",
" d = 10*u.kpc # distance to SN\n",
"\n",
" #get the spectra\n",
" ispec = model.get_initial_spectra(t, energies) \n",
" ospec_nmo = model.get_transformed_spectra(t, energies, xform_nmo)\n",
" ospec_imo = model.get_transformed_spectra(t, energies, xform_imo)\n",
"\n",
" for flavor in Flavor:\n",
" for j, E in enumerate(energies):\n",
" ispec[flavor][j] /= (4.*np.pi*d**2)\n",
" ospec_nmo[flavor][j] /= (4.*np.pi*d**2)\n",
" ospec_imo[flavor][j] /= (4.*np.pi*d**2)\n",
" ispec = model.get_flux(t, energies, d)\n",
" ospec_nmo = model.get_flux(t, energies, d, xform_nmo)\n",
" ospec_imo = model.get_flux(t, energies, d, xform_imo)\n",
" \n",
" fig, axes = plt.subplots(2,2, figsize=(18,15), sharex=True, sharey=True, tight_layout=True)\n",
"\n",
" for i, spec in enumerate([ispec, ospec_nmo]):\n",
" axes[0][0].plot(energies, spec[Flavor.NU_E]/1e16, \n",
" label='Untransformed '+Flavor.NU_E.to_tex() if i==0 else 'Transformed '+Flavor.NU_E.to_tex(),\n",
" color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n",
" axes[0][0].plot(energies, spec[Flavor.NU_X]/1e16, \n",
" label='Untransformed '+Flavor.NU_X.to_tex() if i==0 else 'Transformed '+Flavor.NU_X.to_tex(),\n",
" color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n",
"\n",
" axes[0][0].set(title='Neutrinos in the NMO: $t = ${:.1f}'.format(t))\n",
" axes[0][0].grid()\n",
" axes[0][0].legend(loc='upper right', ncol=2, fontsize=16)\n",
" \n",
" axes[0][1].plot(energies, spec[Flavor.NU_E_BAR]/1e16, \n",
" label='Untransformed '+Flavor.NU_E_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_E_BAR.to_tex(),\n",
" color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n",
" axes[0][1].plot(energies, spec[Flavor.NU_X_BAR]/1e16, \n",
" label='Untransformed '+Flavor.NU_X_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_X_BAR.to_tex(),\n",
" color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n",
"\n",
" axes[0][1].set(title='Antineutrinos in the NMO: $t = ${:.1f}'.format(t)) \n",
" axes[0][1].grid()\n",
" axes[0][1].legend(loc='upper right', ncol=2, fontsize=16) \n",
" \n",
" for i, spec in enumerate([ispec, ospec_imo]):\n",
" axes[1][0].plot(energies, spec[Flavor.NU_E]/1e16, \n",
" label='Untransformed '+Flavor.NU_E.to_tex() if i==0 else 'Transformed '+Flavor.NU_E.to_tex(),\n",
" color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n",
" axes[1][0].plot(energies, spec[Flavor.NU_X]/1e16, \n",
" label='Untransformed '+Flavor.NU_X.to_tex() if i==0 else 'Transformed '+Flavor.NU_X.to_tex(),\n",
" color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n",
"\n",
" axes[1][0].set(xlabel=r'$E$ [{}]'.format(energies.unit), title='Neutrinos in the IMO: $t = ${:.1f}'.format(t))\n",
" axes[1][0].grid()\n",
" axes[1][0].legend(loc='upper right', ncol=2, fontsize=16)\n",
" \n",
" axes[1][1].plot(energies, spec[Flavor.NU_E_BAR]/1e16, \n",
" label='Untransformed '+Flavor.NU_E_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_E_BAR.to_tex(),\n",
" color='C0', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n",
" axes[1][1].plot(energies, spec[Flavor.NU_X_BAR]/1e16, \n",
" label='Untransformed '+Flavor.NU_X_BAR.to_tex() if i==0 else 'Transformed '+Flavor.NU_X_BAR.to_tex(),\n",
" color='C1', ls='-' if i==0 else ':', lw=2, alpha=0.7)\n",
"\n",
" axes[1][1].set(xlabel=r'$E$ [{}]'.format(energies.unit), title='Antineutrinos in the IMO: $t = ${:.1f}'.format(t))\n",
" axes[1][1].grid()\n",
" axes[1][1].legend(loc='upper right', ncol=2, fontsize=18) \n",
"\n",
" ax = axes[0][0]\n",
" ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')\n",
" ax = axes[1][0]\n",
" ax.set(ylabel=r'flux [$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')\n",
" for idx_mixing, ospec in enumerate([ospec_nmo, ospec_imo]):\n",
" for idx_nu, flavors in enumerate([ThreeFlavor['NU'], ThreeFlavor['NU_BAR']]):\n",
" for idx_spec, spec in enumerate([ispec, ospec_nmo]):\n",
" spec_val = spec.array.squeeze()\n",
" spec_val = spec_val.to('1e16/(erg s cm2)')\n",
" for flavor in flavors:\n",
" ax = axes[idx_nu][idx_mixing]\n",
" ax.plot(energies, spec_val[flavor],\n",
" label = ('Untransformed','Transformed')[idx_spec]+' '+flavor.to_tex(),\n",
" color={'E':'C0','MU':'C1','TAU':'C2'}[flavor.lepton],\n",
" ls=('-',':')[idx_spec],\n",
" lw=2, alpha=0.7\n",
" )\n",
" \n",
" ax.set(title='{Nu}eutrinos in the {mixing}: $t = ${t:.1f}'.format(Nu=['N','Antin'][idx_nu],\n",
" mixing=['NMO','IMO'][idx_mixing],\n",
" t=t)\n",
" )\n",
" for ax in axes.flatten():\n",
" ax.grid()\n",
" ax.legend(loc='upper right', ncol=2, fontsize=16)\n",
" for ax in axes[:,0]:\n",
" ax.set(ylabel=f'flux {spec_val.unit._repr_latex_()}')#[$10^{16}$ erg$^{-1}$ cm$^{-2}$ s$^{-1}$]')\n",
" \n",
" return fig"
]
Expand All @@ -247,10 +195,16 @@
"metadata": {},
"outputs": [],
"source": [
"fig = plot_total_flux(model, AdiabaticMSW(MixingParameters()), AdiabaticMSW(MixingParameters('INVERTED')))\n",
"mixpars_nmo = MixingParameters('NORMAL')\n",
"mixpars_imo = MixingParameters('INVERTED')\n",
"\n",
"xform_nmo = xforms.AdiabaticMSW(mixpars_nmo)\n",
"xform_imo = xforms.AdiabaticMSW(mixpars_imo) \n",
"\n",
"fig = plot_total_flux(model, xform_nmo, xform_imo)\n",
"fig.show()\n",
"# fig.savefig('flux_adiabaticmsw.pdf')\n",
"fig = plot_spectra(model, AdiabaticMSW(MixingParameters()), AdiabaticMSW(MixingParameters('INVERTED')), 100*u.ms)\n",
"fig = plot_spectra(model, xform_nmo, xform_imo, 100*u.ms)\n",
"fig.show()\n",
"# fig.savefig('spectra_adiabaticmsw.pdf')"
]
Expand All @@ -261,13 +215,12 @@
"metadata": {},
"outputs": [],
"source": [
"mp_nmo = MixingParameters()\n",
"angles_nmo = mp_nmo.get_mixing_angles() + (10*u.deg,)\n",
"xf_nmo = AdiabaticMSWes(angles_nmo)\n",
"# set 1 degree as the mixing angle theta14\n",
"mixpars_nmo = FourFlavorMixingParameters(**MixingParameters('NORMAL'), theta14=1*u.deg)\n",
"xf_nmo = xforms.AdiabaticMSWes(mixpars_nmo)\n",
"\n",
"mp_imo = MixingParameters(mh=MassHierarchy.INVERTED)\n",
"angles_imo = mp_imo.get_mixing_angles() + (10*u.deg,)\n",
"xf_imo = AdiabaticMSWes(angles_imo)\n",
"mixpars_imo = FourFlavorMixingParameters(**MixingParameters('INVERTED'), theta14=1*u.deg)\n",
"xf_imo = xforms.AdiabaticMSWes(mixpars_imo)\n",
"\n",
"fig = plot_total_flux(model, xf_nmo, xf_imo)\n",
"fig.show()\n",
Expand All @@ -276,11 +229,18 @@
"fig.show()\n",
"# fig.savefig('spectra_adiabaticmswes.pdf')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "snewpy",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -294,9 +254,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.2"
"version": "3.13.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
"nbformat_minor": 4
}
21 changes: 8 additions & 13 deletions doc/source/nb/ccsn/Bollig_2016.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -141,28 +141,23 @@
"\n",
"for i, spec in enumerate([ispec, ospec_nmo]):\n",
" ax = axes[i]\n",
" for flavor in Flavor:\n",
" ax.plot(E, spec[flavor],\n",
" label=flavor.to_tex(),\n",
" color='C0' if flavor.is_electron else 'C1',\n",
" ls='-' if flavor.is_neutrino else ':', lw=2,\n",
" alpha=0.7)\n",
"\n",
" ax.set(xlabel=r'$E$ [{}]'.format(E.unit),\n",
" title='Initial Spectra: $t = ${:.1f}'.format(t) if i==0 else 'Oscillated Spectra: $t = ${:.1f}'.format(t))\n",
" plt.sca(ax)\n",
" spec.plot('energy')\n",
" \n",
" ax.set(title='Initial Spectra: $t = ${:.1f}'.format(t) if i==0 else 'Oscillated Spectra: $t = ${:.1f}'.format(t))\n",
" ax.grid()\n",
" ax.legend(loc='upper right', ncol=2, fontsize=16)\n",
"\n",
"ax = axes[0]\n",
"ax.set(ylabel=r'flux [erg$^{-1}$ s$^{-1}$]')\n",
"ax.set(ylabel=r'flux, MeV')\n",
"\n",
"fig.tight_layout();"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "snewpy",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -176,9 +171,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.2"
"version": "3.13.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
"nbformat_minor": 4
}
21 changes: 6 additions & 15 deletions doc/source/nb/ccsn/Fischer_2020.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,24 +121,15 @@
"\n",
"for i, spec in enumerate([ispec, ospec_nmo]):\n",
" ax = axes[i]\n",
" for flavor in Flavor:\n",
" if flavor.is_electron == True:\n",
" color='C0'\n",
" else:\n",
" color='C1'\n",
" ax.plot(E, spec[flavor],\n",
" label=flavor.to_tex(),\n",
" color=color,\n",
" ls='-' if flavor.is_neutrino else ':', lw=2,\n",
" alpha=0.7)\n",
"\n",
" ax.set(xlabel=r'$E$ [{}]'.format(E.unit),\n",
" title='Initial Spectra: $t = ${:.1f}'.format(t) if i==0 else 'Oscillated Spectra: $t = ${:.1f}'.format(t))\n",
" plt.sca(ax)\n",
" spec.plot('energy')\n",
" \n",
" ax.set(title='Initial Spectra: $t = ${:.1f}'.format(t) if i==0 else 'Oscillated Spectra: $t = ${:.1f}'.format(t))\n",
" ax.grid()\n",
" ax.legend(loc='upper right', ncol=2, fontsize=16)\n",
"\n",
"ax = axes[0]\n",
"ax.set(ylabel=r'flux [erg$^{-1}$ s$^{-1}$]')\n",
"ax.set(ylabel=r'flux, MeV')\n",
"\n",
"fig.tight_layout();"
]
Expand All @@ -160,7 +151,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.3"
"version": "3.13.5"
},
"vscode": {
"interpreter": {
Expand Down
Loading
Loading