.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/giese/giese_lisa_fig2.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_giese_giese_lisa_fig2.py: Giese LISA fig. 2 ================= Reproduction of :giese_2021:`\ `, fig. 2. These figures are used in Mika's M.Sc. thesis. .. GENERATED FROM PYTHON SOURCE LINES 8-264 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_examples/giese/images/sphx_glr_giese_lisa_fig2_001.png :alt: PTtools, $\alpha_{\bar{\theta}_n}$, PTtools, $\alpha_n$, Giese et al., $\alpha_{\bar{\theta}_n}$, Giese et al., $\alpha_n$ :srcset: /auto_examples/giese/images/sphx_glr_giese_lisa_fig2_001.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/giese/images/sphx_glr_giese_lisa_fig2_002.png :alt: $\alpha_{\bar{\theta}_n}$, $\alpha_n$ :srcset: /auto_examples/giese/images/sphx_glr_giese_lisa_fig2_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/giese/images/sphx_glr_giese_lisa_fig2_003.png :alt: PTtools, $\alpha_n$, Giese et al., $\alpha_n$ :srcset: /auto_examples/giese/images/sphx_glr_giese_lisa_fig2_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_examples/giese/images/sphx_glr_giese_lisa_fig2_004.png :alt: giese lisa fig2 :srcset: /auto_examples/giese/images/sphx_glr_giese_lisa_fig2_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/docs/checkouts/readthedocs.org/user_builds/pttools/checkouts/main/pttools/models/const_cs.py:681: RuntimeWarning: invalid value encountered in scalar multiply v_walls [0.2 0.21530612 0.23061224 0.24591837 0.26122449 0.27653061 0.29183673 0.30714286 0.32244898 0.3377551 0.35306122 0.36836735 0.38367347 0.39897959 0.41428571 0.42959184 0.44489796 0.46020408 0.4755102 0.49081633 0.50612245 0.52142857 0.53673469 0.55204082 0.56734694 0.58265306 0.59795918 0.61326531 0.62857143 0.64387755 0.65918367 0.6744898 0.68979592 0.70510204 0.72040816 0.73571429 0.75102041 0.76632653 0.78163265 0.79693878 0.8122449 0.82755102 0.84285714 0.85816327 0.87346939 0.88877551 0.90408163 0.91938776 0.93469388 0.95 ] $\alpha_{\bar{\theta}_n}$ css2=1/3, csb2=1/3 /home/docs/checkouts/readthedocs.org/user_builds/pttools/checkouts/main/examples/giese/giese_lisa_fig2.py:167: RuntimeWarning: All-NaN slice encountered [nan nan nan nan nan nan] css2=1/3, csb2=1/4 [nan nan nan nan nan nan] css2=1/4, csb2=1/3 [nan nan nan nan nan nan] css2=1/4, csb2=1/4 [nan nan nan nan nan nan] $\alpha_n$ css2=1/3, csb2=1/3 [nan nan nan nan nan nan] css2=1/3, csb2=1/4 [nan nan nan nan nan nan] css2=1/4, csb2=1/3 [nan nan nan nan nan nan] css2=1/4, csb2=1/4 [nan nan nan nan nan nan] $\alpha_n$ css2=1/3, csb2=1/3 [nan nan nan nan nan nan] css2=1/3, csb2=1/4 [nan nan nan nan nan nan] css2=1/4, csb2=1/3 [nan nan nan nan nan nan] css2=1/4, csb2=1/4 [nan nan nan nan nan nan] | .. code-block:: Python import logging import time import typing as tp import matplotlib.pyplot as plt import numpy as np from examples.utils import save from pttools.analysis.parallel import create_bubbles # from pttools.analysis.utils import A4_PAPER_SIZE from pttools.bubble.bubble_quantities import get_kappa_giese from pttools.bubble.giese import kappaNuMuModel from pttools.models import ConstCSModel from pttools.speedup import run_parallel, GITHUB_ACTIONS logger = logging.getLogger(__name__) def kappa_giese(params: np.ndarray, model: ConstCSModel) -> float: v_wall, alpha_tbn_giese = params try: kappa, v_arr, wow_arr, xi_arr, mode, vp, vm = kappaNuMuModel( # cs2s=model.cs2(model.w_crit, Phase.SYMMETRIC), # cs2b=model.cs2(model.w_crit, Phase.BROKEN), cs2s=model.css2, cs2b=model.csb2, al=alpha_tbn_giese, vw=v_wall ) except ValueError: return np.nan return kappa def kappas_giese( model: ConstCSModel, v_walls: np.ndarray, alpha_ns: np.ndarray, theta_bar: bool = False) -> np.ndarray: if theta_bar: alpha_tbns = alpha_ns else: alpha_tbns = np.empty((alpha_ns.size,)) for i, alpha_n in enumerate(alpha_ns): try: wn = model.wn(alpha_n, theta_bar=theta_bar) alpha_tbns[i] = model.alpha_theta_bar_n_from_alpha_n(alpha_n=alpha_n, wn=wn) except (ValueError, RuntimeError): alpha_tbns[i] = np.nan params = np.empty((alpha_ns.size, v_walls.size, 2)) for i_alpha_tbn, alpha_tbn in enumerate(alpha_tbns): for i_v_wall, v_wall in enumerate(v_walls): params[i_alpha_tbn, i_v_wall, 0] = v_wall params[i_alpha_tbn, i_v_wall, 1] = alpha_tbn kappas = run_parallel( func=kappa_giese, params=params, multiple_params=True, # output_dtypes=(float, ), # max_workers=max_workers, log_progress_percentage=20, kwargs={"model": model} ) return kappas def create_figure( axs: tp.Iterable[plt.Axes], models: tp.List[ConstCSModel], alpha_ns: np.ndarray, colors: tp.List[str], lss: tp.List[str], v_walls: np.ndarray, theta_bar: bool = False, giese: bool = False): kappas = np.empty((len(models), alpha_ns.size, v_walls.size)) for i_model, (model, ls) in enumerate(zip(models, lss)): if giese: if kappaNuMuModel is None: kappas[i_model, :, :] = np.nan else: kappas[i_model, :, :] = kappas_giese(model=model, v_walls=v_walls, alpha_ns=alpha_ns, theta_bar=theta_bar) else: bubbles, kappas[i_model, :, :] = create_bubbles( model=model, v_walls=v_walls, alpha_ns=alpha_ns, func=get_kappa_giese, bubble_kwargs={"theta_bar": theta_bar, "allow_invalid": False}, allow_bubble_failure=True ) for i_alpha_n, (alpha_n, color) in enumerate(zip(alpha_ns, colors)): try: i_max = np.nanargmax(kappas[i_model, i_alpha_n, :]) except ValueError: logger.error(f"Could not produce bubbles with alpha_n={alpha_n} for {model.label_unicode}") continue kwargs = {} # if ls == "-": # kwargs["label"] = rf"$\alpha={alpha_ns[i_alpha_n]}$" for ax in axs: ax.plot(v_walls, kappas[i_model, i_alpha_n, :], ls=ls, color=color, alpha=0.5, **kwargs) logger.info( f"alpha_n={alpha_n}, kappa_max={kappas[i_model, i_alpha_n, i_max]}, i_max={i_max}, " f"v_wall={v_walls[i_max]}, color={color}, ls={ls}, {model.label_unicode}" ) failed_inds = np.argwhere(np.isnan(kappas[i_model, i_alpha_n, :])) if failed_inds.size: logger.info( "Failed v_walls: %s", v_walls[failed_inds].flatten()) title = "" if giese: title += "Giese et al." else: title += "PTtools" title += ", " if theta_bar: title += r"$\alpha_{\bar{\theta}_n}$" else: title += r"$\alpha_n$" for ax in axs: ax.set_ylabel(r"$\kappa_{\bar{\theta}_n}$") ax.set_ylim(bottom=10 ** -2.5, top=1) ax.set_xlabel(r"$v_\text{wall}$") ax.set_yscale("log") ax.set_xlim(v_walls.min(), v_walls.max()) ax.set_title(title) ax.grid() return kappas def create_diff_figure( ax: plt.Axes, kappas_pttools: np.ndarray, kappas_giese: np.ndarray, models: tp.List[ConstCSModel], v_walls: np.ndarray, colors: tp.List[str], lss: tp.List[str], theta_bar: bool, title: bool = True): rel_diffs = np.abs(kappas_pttools - kappas_giese) / kappas_giese if theta_bar: title_str = r"$\alpha_{\bar{\theta}_n}$" else: title_str = r"$\alpha_n$" print(title_str) for i_model, (model, ls) in enumerate(zip(models, lss)): for i_alpha in range(kappas_pttools.shape[1]): ax.plot( v_walls, rel_diffs[i_model, i_alpha, :], color=colors[i_alpha], ls=ls, # label=f"Model {i_model}, alpha {i_alpha}", ) print(model.label_unicode) print(np.nanmax(rel_diffs[i_model, :, :], axis=1)) ax.set_xlabel(r"$v_\text{wall}$") ax.set_ylabel( r"$|\kappa_{\bar{\theta}_n,\text{PTtools}} - \kappa_{\bar{\theta}_n,\text{ref}}|" r" / \kappa_{\bar{\theta}_n,\text{ref}}$") ax.set_yscale("log") ax.set_ylim(10**(-6), 1) ax.set_xlim(v_walls.min(), v_walls.max()) ax.grid() if title: ax.set_title(title_str) def main(): alpha_ns = np.array([0.01, 0.03, 0.1, 0.3, 1, 3]) colors = ["b", "y", "r", "g", "purple", "grey"] n_v_walls = 20 if GITHUB_ACTIONS else 50 v_walls = np.linspace(0.2, 0.95, n_v_walls) lss = ["-", "--", ":", "-."] a_s = 5 a_b = 1 V_s = 1 models = [ ConstCSModel(css2=1/3, csb2=1/3, a_s=a_s, a_b=a_b, V_s=V_s, alpha_n_min=alpha_ns[0]), ConstCSModel(css2=1/3, csb2=1/4, a_s=a_s, a_b=a_b, V_s=V_s, alpha_n_min=alpha_ns[0]), ConstCSModel(css2=1/4, csb2=1/3, a_s=a_s, a_b=a_b, V_s=V_s, alpha_n_min=alpha_ns[0]), ConstCSModel(css2=1/4, csb2=1/4, a_s=a_s, a_b=a_b, V_s=V_s, alpha_n_min=alpha_ns[0]) ] logger.info(f"Minimum alpha_ns: %s", [model.alpha_n_min for model in models]) for model in models: logger.info("Model parameters: %s", model.params_str()) figsize_x = 8 fig1: plt.Figure = plt.figure(figsize=(figsize_x, 6)) fig2: plt.Figure = plt.figure(figsize=(figsize_x, 4)) fig3: plt.Figure = plt.figure(figsize=(figsize_x, 4)) fig4: plt.Figure = plt.figure(figsize=(figsize_x, 4)) axs1 = fig1.subplots(2, 2) axs2 = fig2.subplots(1, 2) axs3 = fig3.subplots(1, 2) ax4 = fig4.add_subplot() start_time = time.perf_counter() kappas_giese_atbn = create_figure( axs=(axs1[1, 0], ), models=models, alpha_ns=alpha_ns, colors=colors, lss=lss, v_walls=v_walls, theta_bar=True, giese=True ) kappas_giese_an = create_figure( axs=(axs1[1, 1], axs3[1]), models=models, alpha_ns=alpha_ns, colors=colors, lss=lss, v_walls=v_walls, theta_bar=False, giese=True ) giese_time = time.perf_counter() logger.info(f"Creating Giese kappa figures took {giese_time - start_time:.2f} s.") kappas_pttools_atbn = create_figure( axs=(axs1[0, 0], ), models=models, alpha_ns=alpha_ns, colors=colors, lss=lss, v_walls=v_walls, theta_bar=True ) kappas_pttools_an = create_figure( axs=(axs1[0, 1], axs3[0]), models=models, alpha_ns=alpha_ns, colors=colors, lss=lss, v_walls=v_walls, theta_bar=False ) print("v_walls") print(v_walls) create_diff_figure( ax=axs2[0], kappas_pttools=kappas_pttools_atbn, kappas_giese=kappas_giese_atbn, models=models, colors=colors, lss=lss, v_walls=v_walls, theta_bar=True ) create_diff_figure( ax=axs2[1], kappas_pttools=kappas_pttools_an, kappas_giese=kappas_giese_an, models=models, colors=colors, lss=lss, v_walls=v_walls, theta_bar=False ) create_diff_figure( ax=ax4, kappas_pttools=kappas_pttools_an, kappas_giese=kappas_giese_an, models=models, colors=colors, lss=lss, v_walls=v_walls, theta_bar=False, title=False ) logger.info(f"Creating PTtools kappa figures took {time.perf_counter() - giese_time:.2f} s.") fig1.tight_layout() fig2.tight_layout() fig3.tight_layout() fig4.tight_layout() return fig1, fig2, fig3, fig4 if __name__ == "__main__": fig, fig2, fig3, fig4 = main() save(fig, "giese_lisa_fig2") save(fig2, "giese_lisa_fig2_diff") save(fig3, "giese_lisa_fig2_alpha_n") save(fig4, "giese_lisa_fig2_alpha_n_diff") plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (2 minutes 2.074 seconds) **Estimated memory usage:** 267 MB .. _sphx_glr_download_auto_examples_giese_giese_lisa_fig2.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: giese_lisa_fig2.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: giese_lisa_fig2.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: giese_lisa_fig2.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_