test_inference hangs up with suggested parameters

[kelslund]

When running test case notebook test_Inference on a laptop, the code hangs up (1.5 hr+) when the suggested parameters (live_point=2000, max_calls=60000) are used in the step where the sampler is set up. Output includes:

home/lindsey/.local/lib/python3.8/site-packages/InferenceWorkflow/Likelihood.py:115: RuntimeWarning: divide by zero encountered in log
  likelihood = np.log(fx)
/home/lindsey/.local/lib/python3.8/site-packages/TOVsolver/solver_code.py:80: RuntimeWarning: overflow encountered in scalar multiply
  dpdr = dpdr / (r * (r - 2.0 * m))
/home/lindsey/.local/lib/python3.8/site-packages/TOVsolver/solver_code.py:79: RuntimeWarning: overflow encountered in scalar multiply
  dpdr = -(eps + pres) * (m + 4.0 * pi * r**3.0 * pres)
/home/lindsey/.local/lib/python3.8/site-packages/TOVsolver/solver_code.py:80: RuntimeWarning: invalid value encountered in scalar divide
  dpdr = dpdr / (r * (r - 2.0 * m))
/home/lindsey/.local/lib/python3.8/site-packages/scipy/integrate/_ode.py:1179: UserWarning: dopri5: step size becomes too small
  warnings.warn('{:s}: {:s}'.format(self.__class__.__name__,
/home/lindsey/.local/lib/python3.8/site-packages/TOVsolver/solver_code.py:284: RuntimeWarning: overflow encountered in scalar divide
  return stateTOV[1] * unit.c**2 / unit.G, sy.t

[ultranest] Sampling 2000 live points from prior ...

/home/lindsey/.local/lib/python3.8/site-packages/scipy/integrate/_ode.py:1179: UserWarning: dopri5: larger nsteps is needed
  warnings.warn('{:s}: {:s}'.format(self.__class__.__name__,
capi_return is NULL
Call-back cb_fcn_in___user__routines failed.
capi_return is NULL
Call-back cb_fcn_in___user__routines failed.

Changing max_calls to 10000 or 20000 removes these warnings, but still does not complete within 30 minutes,

[kelslund]

This problem has not been unique to the single notebok. @ChunHuangPhy, the text in the notebooks suggest these test cases were run on a laptop. I am also trying to run these on a laptop but the samplers keep hanging up and not giving results. This has been the case for sample_glitchBayesian and test_Bayesian_inference_SpeedofSound_EOS, so far.

[tuhinbits]

Key Improvements in the TOV Solver

1. Comprehensive Documentation

   • Extensively documented function with clear descriptions of parameters, return values, and methodology
   • Added detailed physics notes explaining the numerical approach and potential issues
   • Included explanations of the Taylor expansion used to avoid the central singularity

2. Adaptive Step Sizing Algorithm

   • Implemented the physically-motivated adaptive step sizing with citation:

     delta = 0.23  # optimal delta from https://articles.adsabs.harvard.edu/pdf/1971ApJ...170..299B
     dr = delta / ((1.0 / stateTOV[1]) * dmdr - (1.0 / stateTOV[0]) * dpdr)

3. Multiple Robust Stopping Conditions

   • Added layered stopping criteria to handle various edge cases:

     while (sy.successful() and 
            stateTOV[0] > Pmin and 
            sy.t < max_radius and
            rho_current > 0.01 * center_rho):

   • Density threshold condition to catch surface boundary more reliably
   • Maximum radius safety parameter to prevent runaway integration
   • Mass consistency check to detect numerical instabilities:

     if stateTOV[1] < prev_mass:
         break

4. Surface Region Handling

   • Special step size reduction near the stellar surface:

     if stateTOV[0] < Pmin * 10:
         dr = min(dr, 50.0 * unit.cm)

   • This directly addresses the warnings about step size becoming too small in steep gradient regions

5. Improved Initial Conditions

   • Carefully chosen initial radius (1e-18 * unit.cm) to balance numerical stability concerns
   • Taylor expansion for initial pressure to handle the coordinate singularity at r=0

6. Error Tolerance Configuration

   • Configured the ODE integrator with appropriate absolute and relative error tolerances:

     sy = ode(TOV, None).set_integrator("dopri5", atol=1e-8, rtol=1e-8)

   • Using a more suitable integrator (dopri5) for handling stiff differential equations

These improvements collectively enhance the numerical stability, accuracy, and robustness of the neutron star structure calculations across a wide range of equations of state.

Dynasty Sampler for Notebook Users

Added Dynasty sampler with multiprocessing pool support for efficient notebook execution
Key advantage: Utilizes all CPU cores without requiring MPI configuration
Performance: Achieves near-linear scaling with available cores on a single machine
Maintained UltraNest support for MPI-enabled cluster environments

[tuhinbits]

Scalability Analysis of Bayesian Inference for Neutron Star Equation of State

Executive Summary

This report presents a comprehensive scalability analysis of our Bayesian inference workflow for neutron star equation of state (EOS) investigations using the UltraNest sampler with MPI parallelization. The study was conducted on the Deucalion HPC platform in Portugal, utilizing up to 22 compute nodes and 2,816 CPU cores. Our results demonstrate strong but sub-linear scaling, with performance improvements diminishing as the node count increases.

Introduction

Modern astrophysical investigations of neutron star properties require sophisticated Bayesian inference frameworks to constrain the equation of state from observational data. These computationally intensive workflows benefit significantly from high-performance computing resources but understanding their scaling behavior is crucial for efficient resource allocation.

Our implementation uses a Relativistic Mean Field (RMF) approach with UltraNest sampling to infer parameters of the neutron star equation of state, constrained by mass-radius observations. This report analyzes the performance scaling of this workflow as a function of computational resources.

Methodology

Computational Setup

• HPC Platform: Deucalion HPC platform (Portugal)
• Nodes: 2 to 22 compute nodes
• CPU Cores: 128 cores per node (256 to 2,816 cores total)
• Inference Parameters: 8 physical parameters (g_sigma, g_omega, g_rho, kappa, lambda_0, zeta, Lambda_w, d1)
• Live Points: 2,000
• Maximum Function Calls: 500,000
• Test Code: test_Inference_ultranest_mpi.py

The workflow implements a Bayesian inference using UltraNest with MPI parallelization, calculating the neutron star mass-radius relationship for various parameter combinations and comparing against observational constraints.

Results

Table 1 presents the performance metrics across different node configurations.

Table 1: Scalability Performance Metrics

Nodes	CPUs	Runtime (s)	Runtime (min)	Speedup	Efficiency
2	256	10,938.00	182.30	1.00	1.00
6	768	4,139.00	68.98	2.64	0.88
10	1,280	2,642.00	44.03	4.14	0.83
14	1,792	1,993.00	33.22	5.49	0.78
18	2,304	1,948.00	32.47	5.61	0.62
22	2,816	1,611.00	26.85	6.79	0.62

Where:

• Speedup = Runtime(2 nodes) / Runtime(N nodes)
• Efficiency = Speedup / (N nodes / 2 nodes)

Discussion

Performance Scaling

Our analysis reveals several key insights:

1. Runtime Reduction: Increasing the node count from 2 to 22 reduced the runtime from approximately 3 hours to under 30 minutes, representing a 6.79× speedup.

2. Diminishing Returns: The performance scaling shows clear diminishing returns. While adding nodes consistently improves absolute runtime, the marginal improvement decreases significantly past 14 nodes.

3. Efficiency Drop: Parallel efficiency decreases from 100% (by definition) at the baseline to 62% at 22 nodes, indicating growing communication overhead and load imbalance as the system scales.

4. Performance Plateau: Between 14 and 18 nodes, we observe a performance plateau with minimal runtime improvement (only 45 seconds or 2.3% reduction), suggesting a potential scaling bottleneck.

5. Beyond 18 Nodes: Performance improvement resumes beyond 18 nodes, though still with reduced efficiency.

Computational Bottlenecks

The observed scaling behavior suggests several potential bottlenecks:

1. Communication Overhead: As the number of processes increases, the relative cost of inter-process communication grows, particularly for the MPI-based sampling coordination.

2. Load Imbalance: The dynamic nature of nested sampling can create imbalanced workloads across processors, particularly when some parameter combinations require more complex equation of state evaluations.

3. Memory Bandwidth: The calculation of neutron star structures is memory-intensive, and memory bandwidth limitations may constrain performance at higher node counts.

4. I/O Constraints: At scale, file system operations (reading the crust EOS data, writing posterior samples) may introduce serialization points that limit scaling.

Conclusions and Recommendations

This scalability study demonstrates that our Bayesian inference workflow achieves significant but sub-linear scaling on the Deucalion HPC platform. Based on our findings, we recommend:

1. Optimal Resource Allocation: For time-critical production runs, using 14-18 nodes appears to offer the best balance between absolute performance and resource efficiency. Beyond this point, the diminishing returns may not justify the additional computational resources.

2. Workload Optimization: Future development should focus on reducing inter-process communication and improving load balancing to maintain better efficiency at higher node counts.

3. Hybrid Parallelization: Implementing a hybrid MPI+OpenMP approach could potentially improve scaling by reducing the number of MPI ranks while utilizing all available cores.

4. Parameter Sensitivity Analysis: Investigating which model parameters contribute most to computational load could inform targeted optimizations.

Acknowledgments

This work was supported by the EURO HPC project (EHPC-DEV-2024D12-009). The computational component of this research was carried out on the Deucalion HPC platform in Portugal, whose resources and technical assistance are gratefully acknowledged.