fix: NaN stability, covariance regularization, and cube_samples support

.github/workflows/CI.yml

@@ -72,7 +72,7 @@ jobs:
     - name: Install dependencies
       run: |
         python -m pip install --upgrade pip
-        python -m pip install mpi4py pytest fortranformat anesthetic
+        python -m pip install mpi4py pytest anesthetic
 
     - name: Install pypolychord
       run: pip install -v '.[test]'

.gitignore

@@ -33,3 +33,6 @@ pypolychord/.ld_preload.sh
 # Chains folder (test results only)
 chains/*
 .ipynb_checkpoints/
+CLAUDE.md
+.mcp*
+uv.lock

Makefile

@@ -20,23 +20,33 @@ MPI=1
 # Whether to compile in debugging mode (default: false)
 DEBUG=0
 
-export MPI DEBUG
+# Whether to use LAPACK for eigendecomposition (default: true)
+# Set LAPACK=0 to use pure-Fortran Jacobi fallback
+LAPACK=1
+
+export MPI DEBUG LAPACK
 
 # We can autodetect the compiler type on unix systems via the shell.
 # if you want to override this then just run make with
 # make COMPILER_TYPE=<your type>
-# where <your type> is gnu or intel
+# where <your type> is gnu, intel, intel_llvm, or cray
 ifeq "$(shell which ftn >/dev/null 2>&1; echo $$?)" "0"
 COMPILER_TYPE=cray
-else ifeq "$(shell which ifort >/dev/null 2>&1; echo $$?)" "0" 
+else ifeq "$(shell which ifx >/dev/null 2>&1; echo $$?)" "0"
+# Detected Intel LLVM-based compilers (ifx, icx, icpx)
+COMPILER_TYPE=intel_llvm
+else ifeq "$(shell which ifort >/dev/null 2>&1; echo $$?)" "0"
+# Detected Intel classic compilers (ifort, icc, icpc)
 COMPILER_TYPE=intel
 else ifeq "$(shell which gfortran >/dev/null 2>&1; echo $$?)" "0"
 COMPILER_TYPE=gnu
 endif
 
-ifeq ($(COMPILER_TYPE),intel)
+ifeq ($(COMPILER_TYPE),intel_llvm)
+include Makefile_intel_llvm
+else ifeq ($(COMPILER_TYPE),intel)
 include Makefile_intel
-else ifeq ($(COMPILER_TYPE),gnu) 
+else ifeq ($(COMPILER_TYPE),gnu)
 include Makefile_gnu
 else ifeq ($(COMPILER_TYPE),cray)
 include Makefile_cray
@@ -105,6 +115,8 @@ print_CC:
 	@echo $(CC)
 print_CXX:
 	@echo $(CXX)
+print_FC:
+	@echo $(FC)
 
 CLEANDIRS = $(POLYCHORD_DIR) $(PYPOLYCHORD_DIR) $(LIKELIHOOD_DIR) $(BIN_DIR) $(LIB_DIR) $(DRIVERS_DIR) 
 .PHONY: clean veryclean print_CC print_CXX $(addsuffix clean,$(CLEANDIRS)) $(addsuffix veryclean,$(CLEANDIRS)) 
@@ -118,4 +130,3 @@ veryclean: clean $(addsuffix veryclean,$(CLEANDIRS))
 	$(RM) *~ build dist pypolychord.egg-info pypolychord/*.pyc pypolychord/__pycache__ __pycache__ pypolychord/lib/*.so
 $(addsuffix veryclean,$(CLEANDIRS))  : %veryclean: 
 	$(MAKE) -C $* veryclean
-

Makefile_cray

@@ -19,6 +19,14 @@ FFLAGS += -ffree-line-length-none -fallow-argument-mismatch
 CXXFLAGS += -ffree-line-length-none -fallow-argument-mismatch 
 endif
 
+# LAPACK support for eigendecomposition-based covariance regularization
+# Set LAPACK=0 to disable (will use pure-Fortran Jacobi fallback)
+# Cray systems typically have LAPACK via LibSci
+ifeq ($(LAPACK),1)
+FFLAGS += -DHAVE_LAPACK
+# LibSci is usually linked automatically on Cray systems
+endif
+
 ifdef IPO
 # Archive tool for compiling with ipo.
 AR = xiar rv

Makefile_gnu

@@ -43,6 +43,13 @@ else
 RPATH= 
 endif
 
+# LAPACK support for eigendecomposition-based covariance regularization
+# Set LAPACK=0 to disable (will use pure-Fortran Jacobi fallback)
+ifeq ($(LAPACK),1)
+FFLAGS += -DHAVE_LAPACK
+LDLIBS += -llapack -lblas
+endif
+
 ifeq ($(DEBUG),1)
 # Debugging mode
 # --------------
@@ -55,15 +62,15 @@ ifeq ($(DEBUG),1)
 # backtrace     : produce backtrace of error
 # fpe-trap      : search for floating point exceptions (dividing by zero etc)
 # fbounds-check : check array indices
-FFLAGS += -g -O0 -Wall -Wextra -pedantic -fcheck=all -fimplicit-none -fbacktrace -ffpe-trap=zero,overflow 
+FFLAGS += -g -O0 -Wall -Wextra -pedantic -fcheck=all -fimplicit-none -fbacktrace -ffpe-trap=zero,overflow -DDEBUG
 #
 CXXFLAGS += -g -O0 -Wall -Wextra -Wshadow -Weffc++
 CFLAGS += -g -O0 -Wall -Wextra -Wshadow -Weffc++
 else
 # Optimised mode
 # --------------
 # Ofast : maximum optimisation
-FFLAGS += -Ofast
+FFLAGS += -Ofast -fno-finite-math-only # no-finit-math-only is needed for detectinng NaNs
 CXXFLAGS += -Ofast
 CXFLAGS += -Ofast
 endif

Makefile_intel

@@ -28,6 +28,14 @@ CXXFLAGS += -std=c++11 -fpic
 LDFLAGS += -nofor-main
 LDLIBS += -cxxlib
 
+# LAPACK support for eigendecomposition-based covariance regularization
+# Set LAPACK=0 to disable (will use pure-Fortran Jacobi fallback)
+# Intel MKL provides LAPACK
+ifeq ($(LAPACK),1)
+FFLAGS += -DHAVE_LAPACK
+LDLIBS += -mkl
+endif
+
 ifeq ($(DEBUG),1)
 # Debugging mode
 # --------------
@@ -41,7 +49,7 @@ ifeq ($(DEBUG),1)
 # gen-interfaces : generate an interface block for each routine
 # warn-interfaces: warn on these interface blocks
 # fpe0           : halts program if dividing by zero, etc
-FFLAGS += -g -O0 -traceback -check all,noarg_temp_created -warn all -ftrapuv -debug all -gen-interfaces -warn-interfaces
+FFLAGS += -g -O0 -traceback -check all,noarg_temp_created -warn all -ftrapuv -debug all -gen-interfaces -warn-interfaces -DDEBUG
 CXXFLAGS += -g -O0 -traceback -Wcheck -Wremarks -Wall -Weffc++ -ftrapuv -debug all
 LDFLAGS += -g
 else

Makefile_intel_llvm

@@ -0,0 +1,71 @@
+# Makefile_intel_llvm (for Intel LLVM-based compilers: ifx, icx, icpx)
+
+ifeq ($(MPI),1)
+# Using specific MPI-enabled Intel LLVM compilers
+FC = mpiifx
+CC = mpiicx
+CXX = mpiicpx
+else
+# Using standard Intel LLVM compilers
+FC = ifx
+CXX = icpx
+CC = icx
+endif
+LDSHARED = $(LD) -shared
+
+# CPP uses the C compiler
+CPP = $(CC)
+# LD uses the C/C++ compiler as the primary linker driver for mixed-language projects
+LD = $(CC)
+
+# Archive tool
+AR = ar rv
+
+
+# default flags
+# --------------
+# fpp                    : perform preprocessing
+# fpic                   : shared object libraries
+# heap-arrays            : allocate local arrays on heap
+FFLAGS += -fpp -fpic -heap-arrays
+CXXFLAGS += -std=c++11 -fpic
+LDFLAGS +=
+
+# Link specific Intel runtime libraries
+LDLIBS += -lifcore -lifport -limf -lsvml -liomp5
+
+# LAPACK support for eigendecomposition-based covariance regularization
+# Set LAPACK=0 to disable (will use pure-Fortran Jacobi fallback)
+# Intel MKL provides LAPACK (use -qmkl for LLVM-based compilers)
+ifeq ($(LAPACK),1)
+FFLAGS += -DHAVE_LAPACK -qmkl
+LDLIBS += -qmkl
+endif
+
+ifeq ($(DEBUG),1)
+# Debugging mode
+# --------------
+# g              : enable gnu debugger compatibility
+# O0             : no optimisation
+# traceback      : create a backtrace on system failure
+# check all      : all checks (whilst compiling)
+# warn all       : all warnings (whilst running)
+# ftrapuv        : Traps uninitialized variables by setting them to very large values
+# debug all      : additional debugging information
+# gen-interfaces : generate an interface block for each routine
+# warn-interfaces: warn on these interface blocks
+# fpe0           : halts program if dividing by zero, etc
+FFLAGS += -g -O0 -traceback -check all,noarg_temp_created -warn all -ftrapuv -debug all -gen-interfaces -warn-interfaces -DDEBUG
+CXXFLAGS += -g -O0 -traceback -Wcheck -Wremarks -Wall -Weffc++ -ftrapuv -debug all
+LDFLAGS += -g
+else
+# Optimised mode
+# --------------
+# O2             : A high but safe level of optimisation.
+# fp-model precise: Disables value-unsafe floating-point optimisations. This is crucial for numerical stability and cross-compiler consistency.
+# xHost          : maximise architecture usage (from $(HOST) variable)
+# warn all       : Enable all warnings instead of suppressing them with -w.
+FFLAGS += -O2 -fp-model precise $(HOST) -warn all
+CXXFLAGS += -O2 -fp-model precise $(HOST) -Wall
+
+endif

pypolychord/polychord.py

@@ -521,6 +521,7 @@ def run(loglikelihood, nDims, **kwargs):
         default_comm = None
 
     paramnames = kwargs.pop('paramnames', None)
+    cube_samples = kwargs.pop('cube_samples', None)
 
     default_kwargs = {
         'nDerived': 0,
@@ -583,8 +584,8 @@ def run(loglikelihood, nDims, **kwargs):
                     (kwargs['file_root'] + ".paramnames"))
 
 
-    if 'cube_samples' in kwargs:
-        _make_resume_file(loglikelihood, kwargs['prior'], **kwargs)
+    if cube_samples is not None:
+        _make_resume_file(loglikelihood, cube_samples=cube_samples, **kwargs)
         read_resume = kwargs['read_resume']
         kwargs['read_resume'] = True
 
@@ -644,7 +645,7 @@ def wrap_prior(cube, theta):
                      kwargs['seed'],
                      kwargs['comm'])
 
-    if 'cube_samples' in kwargs:
+    if cube_samples is not None:
         kwargs['read_resume'] = read_resume
 
     try:
@@ -658,14 +659,52 @@ def wrap_prior(cube, theta):
 
 
 
+def _format_fortran_double(val, w=24, d=15, e=3):
+    """Format a float in Fortran E format (E24.15E3).
+
+    Produces output like: ' 0.123456789012345E+001'
+    """
+    import math
+    if math.isnan(val):
+        return f"{'NaN':>{w}}"
+    if math.isinf(val):
+        s = '-Infinity' if val < 0 else 'Infinity'
+        return f"{s:>{w}}"
+    if val == 0.0:
+        sign = '-' if math.copysign(1.0, val) < 0 else ' '
+        return f"{sign}0.{'0' * d}E+{'0' * e}".rjust(w)
+    sign = '-' if val < 0 else ' '
+    aval = abs(val)
+    # Get Python scientific notation with d digits of precision
+    # Python's :.{d}E gives d digits after the decimal with 1 digit before
+    py_str = f"{aval:.{d}E}"
+    # Parse mantissa and exponent from Python format: "d.dddE±ee"
+    mantissa_str, exp_str = py_str.split('E')
+    exp_val = int(exp_str)
+    # Python: d.ddd × 10^exp  →  Fortran: 0.dddd × 10^(exp+1)
+    # Shift decimal: remove '.', prepend '0.'
+    digits = mantissa_str.replace('.', '')  # e.g. "1234567890123456"
+    # digits has d+1 characters; we need d digits after '0.'
+    fortran_mantissa = '0.' + digits[:d]
+    fortran_exp = exp_val + 1
+    exp_sign = '+' if fortran_exp >= 0 else '-'
+    exp_digits = f"{abs(fortran_exp):0{e}d}"
+    result = f"{sign}{fortran_mantissa}E{exp_sign}{exp_digits}"
+    return result.rjust(w)
+
+
+def _format_fortran_int(val, w=12):
+    """Format an integer right-justified in width w."""
+    return f"{int(val):>{w}d}"
+
+
 def _make_resume_file(loglikelihood, **kwargs):
-    import fortranformat as ff
     resume_filename = Path(kwargs['base_dir']) / (kwargs['file_root']
                                                   + ".resume")
 
     try:
         from mpi4py import MPI
-        comm = MPI.COMM_WORLD
+        comm = kwargs.get('comm', None) or MPI.COMM_WORLD
         rank = comm.Get_rank()
         size = comm.Get_size()
     except ImportError:
@@ -675,43 +714,65 @@ def _make_resume_file(loglikelihood, **kwargs):
 
     lives = []
     logL_birth = kwargs['logzero']
-    for i in np.array_split(np.arange(len(kwargs['cube_samples'])), size)[rank]:
-        cube = kwargs['cube_samples'][i]
-        theta = kwargs['prior'](cube)
-        logL = loglikelihood(theta)
-        try:
-            logL, derived = logL
-        except TypeError:
-            derived = []
-        nDims = len(theta)
-        nDerived = len(derived)
-        lives.append(np.concatenate([cube,theta,derived,[logL_birth, logL]]))
-
-    if MPI:
-        sendbuf = np.array(lives).flatten()
+    n_samples = len(kwargs['cube_samples'])
+
+    # In MPI mode with multiple ranks, only workers (rank > 0) evaluate
+    # likelihoods, matching PolyChord's controller/worker split where rank 0
+    # is the administrator. In serial mode, rank 0 does everything.
+    if MPI and size > 1:
+        n_workers = size - 1
+        if rank > 0:
+            worker_index = rank - 1
+            for i in np.array_split(np.arange(n_samples), n_workers)[worker_index]:
+                cube = kwargs['cube_samples'][i]
+                theta = kwargs['prior'](cube)
+                logL = loglikelihood(theta)
+                try:
+                    logL, derived = logL
+                except TypeError:
+                    derived = []
+                nDims = len(theta)
+                nDerived = len(derived)
+                lives.append(np.concatenate([cube,theta,derived,[logL_birth, logL]]))
+
+        sendbuf = np.array(lives).flatten() if lives else np.array([], dtype=np.float64)
         sendcounts = np.array(comm.gather(len(sendbuf)))
         if rank == 0:
             recvbuf = np.empty(sum(sendcounts))
         else:
             recvbuf = None
         comm.Gatherv(sendbuf=sendbuf, recvbuf=(recvbuf, sendcounts), root=0)
 
+        if rank == 0:
+            # Infer row width: nDims (cube) + nDims (theta) + nDerived + 2 (logL_birth, logL)
+            nDims = len(kwargs['cube_samples'][0])
+            nDerived = kwargs.get('nDerived', 0)
+            row_width = 2 * nDims + nDerived + 2
+            lives = np.reshape(recvbuf, (n_samples, row_width))
+    else:
+        for i in range(n_samples):
+            cube = kwargs['cube_samples'][i]
+            theta = kwargs['prior'](cube)
+            logL = loglikelihood(theta)
+            try:
+                logL, derived = logL
+            except TypeError:
+                derived = []
+            nDims = len(theta)
+            nDerived = len(derived)
+            lives.append(np.concatenate([cube,theta,derived,[logL_birth, logL]]))
+        lives = np.array(lives)
+
     if rank == 0:
-        if MPI:
-            lives = np.reshape(recvbuf, (len(kwargs['cube_samples']), len(lives[0])))
-        else:
-            lives = np.array(lives)
         with open(resume_filename,"w") as f:
             def write(var):
                 var = np.atleast_1d(var)
                 if isinstance(var[0], np.integer):
-                    fmt = '(%iI12)' % var.size
-                elif isinstance(var[0], np.double):
-                    fmt = '(%iE24.15E3)' % var.size
+                    f.write(''.join(_format_fortran_int(v) for v in var) + '\n')
+                elif isinstance(var[0], (np.floating, float)):
+                    f.write(''.join(_format_fortran_double(v) for v in var) + '\n')
                 else:
-                    fmt = '(A)'
-                writer = ff.FortranRecordWriter(fmt)
-                f.write(writer.write(var) + '\n')
+                    f.write(str(var[0]) + '\n')
 
             write('=== Number of dimensions ===')
             write(nDims)