diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml new file mode 100644 index 000000000..2d580f542 --- /dev/null +++ b/.github/workflows/tests.yml @@ -0,0 +1,116 @@ +name: CI + +# for the moment only runs on "action* branches". Coverage/pypi not yet set up. + +on: + push: + branches: + - 'action*' + tags: + - '*' + pull_request: + branches: + - '*' + +env: + COBAYA_INSTALL_SKIP: polychord,planck_2015,CamSpec2021,2018_highl_CamSpec,unbinned,keck,classy + COBAYA_PACKAGES_PATH: ../packages + +jobs: + build: + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + include: + - name: "Anaconda: jammy LTS Python (fast; pip CAMB)" + os: ubuntu-latest + pydist: "ANACONDA" + - name: "Latest Python 3.12" + os: ubuntu-latest + python-version: 3.12 + mpi: openmpi + - name: "OS X Python 3.8" + os: macos-latest + python-version: 3.8 + mpi: openmpi + - name: "Windows Python 3.12" + os: windows-latest + python-version: 3.12 + mpi: intelmpi + steps: + - run: ln -s $(which gfortran-14) /usr/local/bin/gfortran + if: matrix.os == 'macos-latest' + + - run: gfortran --version + + - uses: actions/checkout@v4 + with: + fetch-depth: 0 + + - name: Set up Python ${{ matrix.python-version }} + if: matrix.pydist != 'ANACONDA' + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + + - name: Set up Miniconda + if: matrix.pydist == 'ANACONDA' + uses: conda-incubator/setup-miniconda@v3 + with: + auto-activate-base: false + activate-environment: tests-environment + environment-file: tests-environment.yml + + - name: Install mpi + if: matrix.pydist != 'ANACONDA' + uses: mpi4py/setup-mpi@v1 + with: + mpi: ${{ matrix.mpi }} + +# - name: Cache dependencies +# uses: actions/cache@v2 +# with: +# path: | +# ${{ github.workspace }}/packages/data/planck_supp_data_and_covmats +# ${{ github.workspace }}/packages/data/bao_data +# ${{ github.workspace }}/packages/data/sn_data +# ${{ github.workspace }}/packages/data/des_data +# ${{ github.workspace }}/packages/data/planck_2018_pliklite_native +# key: ${{ runner.os }}-build-${{ matrix.python-version }}} + + - name: Install dependencies (pip) + if: matrix.pydist != 'ANACONDA' + run: | + pip install mpi4py -i https://pypi.anaconda.org/mpi4py/simple + pip install -r requirements.txt pytest-xdist pytest-cov flaky matplotlib dill coverage flake8 iminuit numba camb + + - name: Run flake8 + shell: bash -el {0} + run: | + flake8 cobaya --select=E713,E704,E703,E714,E741,E10,E11,E20,E22,E23,E25,E27,E301,E302,E304,E9,F405,F406,F5,F6,F7,F8,W1,W2,W3,W6 --show-source --statistics + + - name: Run cobaya install and tests + shell: bash -el {0} + run: | + coverage run --parallel-mode -m cobaya.install polychord --debug + coverage run --parallel-mode -m pytest tests/ -n auto -k "not cosmo" --skip-not-installed --no-flaky-report + coverage run --parallel-mode -m cobaya.install cosmo-tests --no-progress-bars --debug --skip-global + pytest tests/ --cov -vv -s -k "cosmo" -n 2 --skip-not-installed --no-flaky-report + + - name: Run MPI tests + shell: bash -el {0} + run: | + mpiexec -np 2 coverage run --parallel-mode -m pytest -x -m mpi tests/ --no-flaky-report + + - name: Run external likelihood tests + shell: bash -el {0} + run: | + git clone --depth=1 https://github.com/CobayaSampler/example_external_likelihood + pip install ./example_external_likelihood --quiet + coverage run --parallel-mode -m unittest test_package.tests.test + +# - name: Upload coverage to Codecov +# uses: codecov/codecov-action@v1 + + diff --git a/.travis.yml b/.travis.yml index 17b97eeec..344422544 100644 --- a/.travis.yml +++ b/.travis.yml @@ -8,7 +8,7 @@ env: global: - COBAYA_INSTALL_SKIP=polychord,planck_2015,CamSpec2021,2018_highl_CamSpec,unbinned,keck -if: (type = pull_request) OR (branch = master) OR (branch =~ ^test.*) OR (tag IS present) +if: (type = pull_request) OR (branch = master) OR (branch =~ ^test.*) OR (branch =~ ^alltest.*) OR (tag IS present) #Large CamSpec folders tend to hang, so exclude non-base likelihoods from cache cache: @@ -21,8 +21,6 @@ cache: - /home/travis/build/CosmoPars/packages/data/sn_data - /home/travis/build/CosmoPars/packages/data/des_data - /home/travis/build/CosmoPars/packages/data/planck_2018_pliklite_native - - /home/travis/build/CosmoPars/packages/data/planck_2018_lowE_native - - /home/travis/build/CosmoPars/packages/data/planck_2018_lowT_native # (Pre)Installation jobs: @@ -46,22 +44,21 @@ jobs: apt: packages: - gfortran + install: + - pip install -r requirements.txt camb env: - GCC_VERSION="ubuntu" python: "3.10" - - if: branch !~ ^test.* - name: "Anaconda: jammy LTS Python 3.11" + - name: "Anaconda: jammy LTS Python (fast; pip CAMB)" dist: jammy - addons: - apt: - packages: - - gfortran env: - GCC_VERSION="ubuntu" - PYDIST="ANACONDA" - ANACONDA_CHANNEL="defaults" - python: "3.11" - - name: "Latest jammy Python 3.12" + - COBAYA_INSTALL_SKIP="$COBAYA_INSTALL_SKIP,classy" + language: minimal + - if: branch !~ ^test.* + name: "Latest jammy Python 3.12" dist: jammy addons: apt: @@ -69,6 +66,7 @@ jobs: - gfortran env: - GCC_VERSION="ubuntu" + - COBAYA_INSTALL_SKIP="$COBAYA_INSTALL_SKIP,classy" python: "3.12" @@ -81,27 +79,25 @@ before_install: ln -s /usr/bin/g++-$GCC_VERSION gcc-symlinks/g++; export PATH=$PWD/gcc-symlinks:$PATH; fi - - if [[ "$GCC_VERSION" == "11" ]]; then - export LD_PRELOAD=/usr/lib/x86_64-linux-gnu/libgfortran.so.5; - fi - - gfortran --version + - which gfortran >/dev/null 2>&1 && gfortran --version || echo "gfortran not installed" # Install rest of system requisites - - sudo apt install openmpi-bin openmpi-common libopenmpi-dev libopenblas-dev liblapack-dev + # - sudo apt install openmpi-bin openmpi-common libopenmpi-dev libopenblas-dev liblapack-dev # Python requisites - if [[ "$PYDIST" == "ANACONDA" ]]; then - wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh; + wget -q https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh; bash miniconda.sh -b -p $HOME/miniconda; export PATH="$HOME/miniconda/bin:$PATH"; hash -r; conda config --set always_yes yes --set changeps1 no; conda info -a; - conda create -q -n test-environment -c $ANACONDA_CHANNEL python=$TRAVIS_PYTHON_VERSION scipy matplotlib cython PyYAML dill coverage pytest; + conda create -q -n test-environment -c $ANACONDA_CHANNEL scipy matplotlib cython PyYAML dill coverage pytest pandas; source activate test-environment; conda install -c conda-forge mpi4py openmpi iminuit; - pip install flake8 flaky pytest-forked pytest-cov; + pip install -r requirements.txt flake8 flaky pytest-xdist pytest-cov camb; else python -m pip install --upgrade pip setuptools wheel; - pip install mpi4py "numpy<2" pytest-forked pytest-cov flaky matplotlib dill coverage flake8 iminuit; + pip install openmpi mpi4py -i https://pypi.anaconda.org/mpi4py/simple; + pip install pytest-xdist pytest-cov flaky matplotlib dill coverage flake8 iminuit numba; fi - python --version @@ -112,11 +108,9 @@ script: # General tests: - export COBAYA_PACKAGES_PATH="../packages" - coverage run --parallel-mode -m cobaya.install polychord --debug - - coverage run --parallel-mode -m pytest tests/ -k "not cosmo" --skip-not-installed --no-flaky-report - # numba can speed things, but makes for some messy logs - # - pip install numba + - coverage run --parallel-mode -m pytest tests/ -n auto -k "not cosmo" --skip-not-installed --no-flaky-report # Cosmology tests: - - coverage run --parallel-mode -m cobaya.install cosmo-tests --no-progress-bars --debug + - coverage run --parallel-mode -m cobaya.install cosmo-tests --no-progress-bars --debug --skip-global - if [ -n "${CAMB_BRANCH}" ]; then rm -rf $COBAYA_PACKAGES_PATH/code/CAMB ; git clone --recursive --depth 1 -b $CAMB_BRANCH https://github.com/cmbant/CAMB $COBAYA_PACKAGES_PATH/code/CAMB ; @@ -125,10 +119,10 @@ script: # mpi tests - mpiexec -np 2 --mca orte_base_help_aggregate 0 --mca btl ^openib --oversubscribe coverage run --parallel-mode -m pytest -x -m mpi tests/ --no-flaky-report ; - mkdir covers; mv .coverage.* covers; ls -ltra covers - - pytest tests/ --cov -vv -s -k "cosmo" --forked --skip-not-installed --no-flaky-report + - pytest tests/ --cov -vv -s -k "cosmo" -n 1 --skip-not-installed --no-flaky-report - mv .coverage .coverage.pytest; mv covers/.cov* . # Test external cosmological likelihoods - - pip install -e $COBAYA_PACKAGES_PATH/code/CAMB --quiet + #- pip install -e $COBAYA_PACKAGES_PATH/code/CAMB --quiet #- git clone --depth=1 https://github.com/CobayaSampler/planck_lensing_external #- pip install ./planck_lensing_external --quiet #- coverage run --parallel-mode -m unittest plancklensing.tests.test_likes diff --git a/CHANGELOG.md b/CHANGELOG.md index fc3bc3cf4..dc8488032 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,8 @@ +## 3.5.x + +- Detect and fix incomplete last lines when resuming or minimizing from existing runs (#306, #378) +- Added functions module and refactored some numerical functions into it + ## 3.5.4 - Allow classes to have both yaml and class attributes as long as no duplicate keys diff --git a/DEVEL.rst b/DEVEL.rst index e594f2af3..6b827e624 100644 --- a/DEVEL.rst +++ b/DEVEL.rst @@ -21,7 +21,7 @@ Development flow for contributors 1. Fork and clone the repo from github. 2. From its folder, install in editable mode: ``pip install -e .`` 3. Modify stuff. -4. Test with pytest (first "pip install pytest pytest-forked pytest-cov flaky dill") +4. Test with pytest (first "pip install pytest pytest-xdist pytest-cov flaky dill") 5. Make a pull requests and check (about about 15 minutes) if the tests have passed. 6. Iterate until tests pass, then wait for or request feedback/merge diff --git a/TODO.md b/TODO.md index 3f6403ade..7cd8a03ab 100644 --- a/TODO.md +++ b/TODO.md @@ -5,7 +5,7 @@ # cosmetic/consistency/speed ## min/max bounds enforced on derived parameters (more generally, "bounds" as well as priors) -## make portalocker/numba/dill a requirement? +## make numba a requirement? ## version attribute should be in all components not just theory (samplers can have versions) [done for samplers; missing: likelihoods] ## In the docs "Bases" (and UML diagram) not hyperlinked correctly (not sure how to fix) ## dump log info along with each chain file if saving to file (currently in stdout) diff --git a/cobaya/__main__.py b/cobaya/__main__.py index 23ce085c0..5d4a7f591 100644 --- a/cobaya/__main__.py +++ b/cobaya/__main__.py @@ -13,7 +13,7 @@ def run_command(): prefix = "cobaya-" console_scripts = ( metadata.entry_points().select(group="console_scripts") - if sys.version_info > (3, 9) + if sys.version_info >= (3, 10) else metadata.entry_points()["console_scripts"] ) for script in console_scripts: diff --git a/cobaya/collection.py b/cobaya/collection.py index 8da770d6c..e0c5d11ca 100644 --- a/cobaya/collection.py +++ b/cobaya/collection.py @@ -96,7 +96,7 @@ def compute_temperature(logpost, logprior, loglike, check=True, extra_tolerance= """ Returns the temperature of a sample. - If ``check=True`` and the log-probabilites passed are arrays, checks consistency + If ``check=True`` and the log-probabilities passed are arrays, checks consistency of the sample temperature, and raises ``AssertionError`` if inconsistent. """ temp = (logprior + loglike) / logpost @@ -286,8 +286,7 @@ def __init__(self, model, output=None, cache_size=_default_cache_size, name=None self.reset() # If loaded, check sample weights, consistent logp sums, # and temperature (ignores the given one) - samples_loaded = len(self) > 0 - if samples_loaded: + if len(self) > 0: try: try: self.temperature = self._check_logps(extra_tolerance=False) @@ -499,8 +498,9 @@ def _check_logps(self, temperature_only=False, extra_tolerance=False): check=True, extra_tolerance=extra_tolerance) except AssertionError as excpt: raise LoggedError( - self.log, "The sample seems to have an inconsistent temperature.") \ - from excpt + self.log, "The sample seems to have an inconsistent temperature. " + "This could be due to input file truncation on the last line " + "due to crash/being killed before complete.") from excpt if not temperature_only: tols = { "rtol": 1e-4 * (10 if extra_tolerance else 1), diff --git a/cobaya/component.py b/cobaya/component.py index acc81b517..98799cb20 100644 --- a/cobaya/component.py +++ b/cobaya/component.py @@ -444,9 +444,7 @@ def compare_versions(cls, version_a, version_b, equal=True): :return: bool """ va, vb = version.parse(version_a), version.parse(version_b) - if va >= vb if equal else va > vb: - return True - return False + return va >= vb if equal else va > vb def __exit__(self, exception_type, exception_value, traceback): if self.timer and self.timer.n: @@ -482,6 +480,7 @@ def dump_timing(self): def get_versions(self, add_version_field=False) -> InfoDict: """ Get version dictionary + :return: dictionary of versions for all components """ @@ -494,6 +493,7 @@ def format_version(x): def get_speeds(self, ignore_sub=False) -> InfoDict: """ Get speeds dictionary + :return: dictionary of versions for all components """ from cobaya.theory import HelperTheory @@ -695,8 +695,7 @@ def module_class_for_name(m, name): for cls in classes_in_module(m, subclass_of=CobayaComponent): if cls.__name__.lower() in valid_names: if result is not None: - raise ValueError('More than one class with same lowercase name %s', - name) + raise ValueError(f'More than one class with same lowercase name {name}') result = cls return result @@ -777,7 +776,7 @@ def _bare_load_external_module(name, path=None, min_version=None, reload=False, def load_external_module(module_name=None, path=None, install_path=None, min_version=None, get_import_path=None, reload=False, logger=None, - not_installed_level=None): + not_installed_level=None, default_global=False): """ Tries to load an external module at initialisation, dealing with explicit paths and Cobaya's installation path. @@ -807,12 +806,14 @@ def load_external_module(module_name=None, path=None, install_path=None, min_ver found. If this exception will be handled at a higher level, you may pass `not_installed_level='debug'` to prevent printing non-important messages at error-level logging. + + If default_global=True, always attempts to load from the global path if not + installed at path (e.g. pip install). """ if not logger: logger = get_logger(__name__) load_kwargs = {"name": module_name, "path": path, "get_import_path": get_import_path, "min_version": min_version, "reload": reload, "logger": logger} - default_global = False if isinstance(path, str): if path.lower() == "global": msg_tried = "global import (`path='global'` given)" diff --git a/cobaya/containers.py b/cobaya/containers.py index 8830dd723..9bbc161b9 100644 --- a/cobaya/containers.py +++ b/cobaya/containers.py @@ -47,7 +47,7 @@ # Python requisites -- LC_ALL=C: Necessary just for pip <= 8.1.2 (Xenial version) ENV LC_ALL C RUN python -m pip install --upgrade pip -RUN python -m pip install pytest-forked matplotlib cython astropy --upgrade +RUN python -m pip install pytest-xdist matplotlib cython astropy --upgrade # Prepare environment and tree for external packages ------------------------- ENV LD_LIBRARY_PATH $LD_LIBRARY_PATH:/usr/local/lib ENV CONTAINED TRUE diff --git a/cobaya/cosmo_input/autoselect_covmat.py b/cobaya/cosmo_input/autoselect_covmat.py index 625ee4c3e..90d45807c 100644 --- a/cobaya/cosmo_input/autoselect_covmat.py +++ b/cobaya/cosmo_input/autoselect_covmat.py @@ -205,7 +205,7 @@ def score_likes(_key: CovmatFileKey, covmat): best_p_l = get_best_score(best_p, score_likes) if is_debug(log): log.debug("Subset based on params + likes:\n - " + - "\n - ".join([b["name"] for b in best_p_l])) + "\n - ".join([b["name"] for b in best_p_l.values()])) if key_tuple: def score_left_params(_key, _covmat): diff --git a/cobaya/functions.py b/cobaya/functions.py new file mode 100644 index 000000000..791478d8e --- /dev/null +++ b/cobaya/functions.py @@ -0,0 +1,91 @@ +import numpy as np +import logging +import scipy + +try: + import numba +except (ImportError, SystemError): + # SystemError caused usually by incompatible numpy version + + numba = None + logging.debug("Numba not available, install it for better performance.") + + from scipy.stats import special_ortho_group + + random_SO_N = special_ortho_group.rvs + +else: + import warnings + + + def random_SO_N(dim, random_state): + """ + Draw random samples from SO(N). + Equivalent to scipy function but about 10x faster + Parameters + ---------- + dim : integer + Dimension of rotation space (N). + random_state: generator + Returns + ------- + rvs : Random size N-dimensional matrices, dimension (dim, dim) + + """ + dim = np.int64(dim) + H = np.eye(dim) + xx = random_state.standard_normal(size=(dim + 2) * (dim - 1) // 2) + _rvs(dim, xx, H) + return H + + + logging.getLogger('numba').setLevel(logging.ERROR) + with warnings.catch_warnings(): + warnings.filterwarnings("ignore") + + + @numba.njit("void(int64,float64[::1],float64[:,::1])") + def _rvs(dim, xx, H): + D = np.empty((dim,)) + ix = 0 + for n in range(dim - 1): + x = xx[ix:ix + dim - n] + ix += dim - n + norm2 = np.dot(x, x) + x0 = x[0].item() + D[n] = np.sign(x[0]) if x[0] != 0 else 1 + x[0] += D[n] * np.sqrt(norm2) + x /= np.sqrt((norm2 - x0 ** 2 + x[0] ** 2) / 2.) + # Householder transformation + tmp = np.dot(H[:, n:], x) + H[:, n:] -= np.outer(tmp, x) + D[-1] = (-1) ** (dim - 1) * D[:-1].prod() + H[:, :] = (D * H.T).T + + +def chi_squared(c_inv, delta): + """ + Compute chi squared, i.e. delta.T @ c_inv @ delta + + :param c_inv: symmetric positive definite inverse covariance matrix + :param delta: 1D array + :return: delta.T @ c_inv @ delta + """ + if len(delta) < 1500: + return c_inv.dot(delta).dot(delta) + else: + # use symmetry + return scipy.linalg.blas.dsymv(alpha=1.0, + a=c_inv if np.isfortran(c_inv) else c_inv.T, + x=delta, lower=0).dot(delta) + + +def inverse_cholesky(cov): + """ + Get inverse of Cholesky decomposition + + :param cov: symmetric positive definite matrix + :return: L^{-1} where cov = L L^T + """ + cholesky = np.linalg.cholesky(cov) + return scipy.linalg.lapack.dtrtri(cholesky, lower=True)[0] diff --git a/cobaya/install.py b/cobaya/install.py index a633a1b8c..548d2e9c7 100644 --- a/cobaya/install.py +++ b/cobaya/install.py @@ -425,7 +425,7 @@ def download_file(url, path, *, size=None, no_progress_bars=False, logger=None): if not no_progress_bars: bar = tqdm.tqdm(total=size, unit='iB', unit_scale=True, unit_divisor=1024) with open(filename_tmp_path, 'wb') as f: - for data in req.iter_content(chunk_size=1024): + for data in req.iter_content(chunk_size=32768): chunk_size = f.write(data) if not no_progress_bars: bar.update(chunk_size) @@ -455,7 +455,7 @@ def download_file(url, path, *, size=None, no_progress_bars=False, logger=None): if extension == "tgz": extension = "gz" kwargs = {} - if sys.version_info >= (12, 0): + if sys.version_info >= (3, 12): kwargs['filter'] = 'data' with tarfile.open(filename_tmp_path, "r:" + extension) as tar: tar.extractall(path, **kwargs) @@ -533,7 +533,7 @@ def pip_install(packages, upgrade=False, logger=None, options=(), **kwargs): """ if hasattr(packages, "split"): packages = [packages] - cmd = [sys.executable, '-m', 'pip', 'install'] + cmd = [sys.executable, '-m', 'pip', 'install', '-q'] if upgrade: cmd += ['--upgrade'] cmd += list(options) diff --git a/cobaya/likelihoods/base_classes/DataSetLikelihood.py b/cobaya/likelihoods/base_classes/DataSetLikelihood.py index cac4db05b..86ca331c8 100644 --- a/cobaya/likelihoods/base_classes/DataSetLikelihood.py +++ b/cobaya/likelihoods/base_classes/DataSetLikelihood.py @@ -14,8 +14,6 @@ from cobaya.typing import InfoDict from cobaya.conventions import packages_path_input from .InstallableLikelihood import InstallableLikelihood -# import _fast_chi_square for backwards compatibility -from .InstallableLikelihood import ChiSquaredFunctionLoader as _fast_chi_square # noqa class DataSetLikelihood(InstallableLikelihood): diff --git a/cobaya/likelihoods/base_classes/InstallableLikelihood.py b/cobaya/likelihoods/base_classes/InstallableLikelihood.py index 11e1ad9cf..42e3a7af4 100644 --- a/cobaya/likelihoods/base_classes/InstallableLikelihood.py +++ b/cobaya/likelihoods/base_classes/InstallableLikelihood.py @@ -17,27 +17,7 @@ from cobaya.install import _version_filename from cobaya.component import ComponentNotInstalledError from cobaya.tools import VersionCheckError, resolve_packages_path - - -class ChiSquaredFunctionLoader: - - def __init__(self, method_name='_fast_chi_squared'): - """ - On first use will replace method_name with direct reference to chi2 function - :param method_name: name of the property that _fast_chi_square() is assigned to - """ - self._method_name = method_name - - def __get__(self, instance, owner): - # delay testing active camb until run time - try: - from camb.mathutils import chi_squared as fast_chi_squared - except ImportError: - def fast_chi_squared(covinv, x): - return covinv.dot(x).dot(x) - - setattr(instance, self._method_name, fast_chi_squared) - return fast_chi_squared +from cobaya.functions import chi_squared class InstallableLikelihood(Likelihood): @@ -49,7 +29,7 @@ class InstallableLikelihood(Likelihood): # fast convenience function, to get chi-squared (exploiting symmetry); can call # self._fast_chi_squared(cov_inv, delta) - _fast_chi_squared = ChiSquaredFunctionLoader() + _fast_chi_squared = staticmethod(chi_squared) def __init__(self, *args, **kwargs): # Ensure check for install and version errors diff --git a/cobaya/likelihoods/base_classes/__init__.py b/cobaya/likelihoods/base_classes/__init__.py index 6955b94fd..dec992bf6 100644 --- a/cobaya/likelihoods/base_classes/__init__.py +++ b/cobaya/likelihoods/base_classes/__init__.py @@ -1,6 +1,6 @@ from .InstallableLikelihood import InstallableLikelihood from .bao import BAO -from .DataSetLikelihood import DataSetLikelihood, _fast_chi_square +from .DataSetLikelihood import DataSetLikelihood from .cmblikes import CMBlikes, make_forecast_cmb_dataset from .des import DES from .H0 import H0 diff --git a/cobaya/likelihoods/base_classes/des.py b/cobaya/likelihoods/base_classes/des.py index adb3a93c9..592e0d910 100644 --- a/cobaya/likelihoods/base_classes/des.py +++ b/cobaya/likelihoods/base_classes/des.py @@ -56,6 +56,7 @@ from cobaya.likelihoods.base_classes import DataSetLikelihood from cobaya.log import LoggedError from cobaya.conventions import Const +from cobaya.functions import numba # DES data types def_DES_types = ['xip', 'xim', 'gammat', 'wtheta'] @@ -129,26 +130,27 @@ def get_def_cuts(): # pragma: no cover return ranges -try: - import numba -except ImportError: - numba = None -else: +if numba: + import warnings - @numba.njit("void(float64[::1],float64[::1],float64[::1],float64[::1])") - def _get_lensing_dots(wq_b, chis, n_chi, dchis): - for i, chi in enumerate(chis): - wq_b[i] = np.dot(n_chi[i:], (1 - chi / chis[i:]) * dchis[i:]) + with warnings.catch_warnings(): + warnings.filterwarnings("ignore") - @numba.njit("void(float64[:,::1],float64[:,::1],float64[::1],float64)") - def _limber_PK_terms(powers, ks, dchifac, kmax): - for ix in range(powers.shape[0]): - for i in range(powers.shape[1]): - if 1e-4 <= ks[ix, i] < kmax: - powers[ix, i] *= dchifac[i] - else: - powers[ix, i] = 0 + @numba.njit("void(float64[::1],float64[::1],float64[::1],float64[::1])") + def _get_lensing_dots(wq_b, chis, n_chi, dchis): + for i, chi in enumerate(chis): + wq_b[i] = np.dot(n_chi[i:], (1 - chi / chis[i:]) * dchis[i:]) + + + @numba.njit("void(float64[:,::1],float64[:,::1],float64[::1],float64)") + def _limber_PK_terms(powers, ks, dchifac, kmax): + for ix in range(powers.shape[0]): + for i in range(powers.shape[1]): + if 1e-4 <= ks[ix, i] < kmax: + powers[ix, i] *= dchifac[i] + else: + powers[ix, i] = 0 class DES(DataSetLikelihood): diff --git a/cobaya/likelihoods/gaussian_mixture/gaussian_mixture.py b/cobaya/likelihoods/gaussian_mixture/gaussian_mixture.py index 385300bf1..a4ab4457f 100644 --- a/cobaya/likelihoods/gaussian_mixture/gaussian_mixture.py +++ b/cobaya/likelihoods/gaussian_mixture/gaussian_mixture.py @@ -16,6 +16,7 @@ from cobaya.log import LoggedError from cobaya.mpi import share_mpi, is_main_process from cobaya.typing import InputDict, Union, Sequence +from cobaya.functions import inverse_cholesky derived_suffix = "_derived" @@ -111,7 +112,7 @@ def initialize_with_params(self): self.weights = 1 / len(self.gaussians) # Prepare the transformation(s) for the derived parameters - self.inv_choleskyL = [np.linalg.inv(np.linalg.cholesky(cov)) for cov in self.covs] + self.inv_choleskyL = [inverse_cholesky(cov) for cov in self.covs] def logp(self, **params_values): """ diff --git a/cobaya/samplers/mcmc/mcmc.py b/cobaya/samplers/mcmc/mcmc.py index 00b8ecca8..7f7d51c79 100644 --- a/cobaya/samplers/mcmc/mcmc.py +++ b/cobaya/samplers/mcmc/mcmc.py @@ -24,6 +24,7 @@ from cobaya.samplers.mcmc.proposal import BlockedProposer from cobaya.log import LoggedError, always_stop_exceptions from cobaya.tools import get_external_function, NumberWithUnits, load_DataFrame +from cobaya.functions import inverse_cholesky from cobaya.yaml import yaml_dump_file from cobaya.model import LogPosterior from cobaya import mpi @@ -724,7 +725,7 @@ def check_convergence_and_learn_proposal(self): converged_means = False # Cholesky of (normalized) mean of covs and eigvals of Linv*cov_of_means*L try: - L = np.linalg.cholesky(norm_mean_of_covs) + Linv = inverse_cholesky(norm_mean_of_covs) except np.linalg.LinAlgError: self.log.warning( "Negative covariance eigenvectors. " @@ -732,7 +733,6 @@ def check_convergence_and_learn_proposal(self): "contain enough information at this point. " "Skipping learning a new covmat for now.") else: - Linv = np.linalg.inv(L) try: eigvals = np.linalg.eigvalsh(Linv.dot(corr_of_means).dot(Linv.T)) success_means = True diff --git a/cobaya/samplers/mcmc/proposal.py b/cobaya/samplers/mcmc/proposal.py index 8d95876cb..715ff2052 100644 --- a/cobaya/samplers/mcmc/proposal.py +++ b/cobaya/samplers/mcmc/proposal.py @@ -16,13 +16,12 @@ See https://arxiv.org/abs/1304.4473 """ -import logging from itertools import chain - import numpy as np from cobaya.log import LoggedError, HasLogger -from cobaya.tools import choleskyL +from cobaya.tools import choleskyL_corr +from cobaya.functions import random_SO_N class IndexCycler: @@ -55,63 +54,6 @@ def next(self): return self.indices[self.loop_index] -try: - import numba -except (ImportError, SystemError): - # SystemError caused usually by incompatible numpy version - from scipy.stats import special_ortho_group - - random_SO_N = special_ortho_group.rvs - numba = None -else: - import warnings - - - def random_SO_N(dim, random_state): - """ - Draw random samples from SO(N). - Equivalent to scipy function but about 10x faster - Parameters - ---------- - dim : integer - Dimension of rotation space (N). - random_state: generator - Returns - ------- - rvs : Random size N-dimensional matrices, dimension (dim, dim) - - """ - dim = np.int64(dim) - H = np.eye(dim) - xx = random_state.standard_normal(size=(dim + 2) * (dim - 1) // 2) - _rvs(dim, xx, H) - return H - - - logging.getLogger('numba').setLevel(logging.ERROR) - with warnings.catch_warnings(): - warnings.filterwarnings("ignore") - - - @numba.njit("void(int64,float64[::1],float64[:,::1])") - def _rvs(dim, xx, H): - D = np.empty((dim,)) - ix = 0 - for n in range(dim - 1): - x = xx[ix:ix + dim - n] - ix += dim - n - norm2 = np.dot(x, x) - x0 = x[0].item() - D[n] = np.sign(x[0]) if x[0] != 0 else 1 - x[0] += D[n] * np.sqrt(norm2) - x /= np.sqrt((norm2 - x0 ** 2 + x[0] ** 2) / 2.) - # Householder transformation - tmp = np.dot(H[:, n:], x) - H[:, n:] -= np.outer(tmp, x) - D[-1] = (-1) ** (dim - 1) * D[:-1].prod() - H[:, :] = (D * H.T).T - - class RandDirectionProposer(IndexCycler): def propose_vec(self, scale: float = 1): """ @@ -275,7 +217,7 @@ def set_covariance(self, propose_matrix): "symmetric square matrix.") self.propose_matrix = propose_matrix.copy() propose_matrix_j_sorted = self.propose_matrix[np.ix_(self.i_of_j, self.i_of_j)] - sigmas_diag, L = choleskyL(propose_matrix_j_sorted, return_scale_free=True) + sigmas_diag, L = choleskyL_corr(propose_matrix_j_sorted) # Store the basis as transformation matrices self.transform = [] for j_start, bp in zip(self.j_start, self.proposer): diff --git a/cobaya/theories/camb/camb.py b/cobaya/theories/camb/camb.py index fc3d09b1d..cfe0c4e18 100644 --- a/cobaya/theories/camb/camb.py +++ b/cobaya/theories/camb/camb.py @@ -263,7 +263,8 @@ class CAMB(BoltzmannBase): def initialize(self): """Importing CAMB from the correct path, if given.""" try: - install_path = (lambda p: self.get_path(p) if p else None)(self.packages_path) + install_path = self.get_path(self.packages_path) \ + if self.packages_path else None min_version = None if self.ignore_obsolete else self._min_camb_version self.camb = load_external_module( "camb", path=self.path, install_path=install_path, @@ -1017,8 +1018,7 @@ def get_allow_agnostic(self): def must_provide(self, **requirements): super().must_provide(**requirements) - opts = requirements.get('CAMB_transfers') - if opts: + if opts := requirements.get('CAMB_transfers'): self.non_linear_sources = opts['non_linear'] self.needs_perts = opts['needs_perts'] self.cobaya_camb.check_no_repeated_input_extra() diff --git a/cobaya/tools.py b/cobaya/tools.py index 9cb21ead7..7858a1091 100644 --- a/cobaya/tools.py +++ b/cobaya/tools.py @@ -457,6 +457,19 @@ def read_dnumber(n: Any, dim: int): return NumberWithUnits(n, "d", dtype=int, scale=dim).value +def truncate_to_end_line(file_name): + with open(file_name, "r+b") as inp: + # Find the last complete line + inp.seek(0, 2) # Go to the end of the file + pos = inp.tell() - 1 + while pos > 0 and inp.read(1) != b"\n": + pos -= 1 + inp.seek(pos, 0) + if pos > 0: + inp.seek(pos + 1, 0) + inp.truncate() + + def load_DataFrame(file_name, skip=0, root_file_name=None): """ Loads a `pandas.DataFrame` from a text file @@ -493,7 +506,17 @@ def load_DataFrame(file_name, skip=0, root_file_name=None): inp, sep=" ", header=None, names=cols, comment="#", skipinitialspace=True, skiprows=skip, index_col=False) - return data + if not data.empty: + # Check if the last row contains any NaNs + if data.iloc[-1].isna().any(): + log.warning("Last row of %s is incomplete or contains NaNs", file_name) + # If the second-to-last row exists and doesn't contain NaNs, + # delete the last row assuming this was due to crash on write + if len(data) > 1 and not data.iloc[-2].isna().any(): + data = data.iloc[:-1] + log.info(f"Saving {file_name} deleting last (in)complete line") + truncate_to_end_line(file_name) + return data def prepare_comment(comment): @@ -610,14 +633,14 @@ def fast_logpdf(x): def _KL_norm(m1, S1, m2, S2): - """Performs the Guassian KL computation, without input testing.""" + """Performs the Gaussian KL computation, without input testing.""" dim = S1.shape[0] S2inv = np.linalg.inv(S2) return 0.5 * (np.trace(S2inv.dot(S1)) + (m1 - m2).dot(S2inv).dot(m1 - m2) - - dim + np.log(np.linalg.det(S2) / np.linalg.det(S1))) + dim + np.linalg.slogdet(S2)[1] - np.linalg.slogdet(S1)[1]) -def KL_norm(m1=None, S1=np.array([]), m2=None, S2=np.array([]), symmetric=False): +def KL_norm(m1=None, S1=(), m2=None, S2=(), symmetric=False): """Kullback-Leibler divergence between 2 gaussians.""" S1, S2 = [np.atleast_2d(S) for S in [S1, S2]] assert S1.shape[0], "Must give at least S1" @@ -634,37 +657,34 @@ def KL_norm(m1=None, S1=np.array([]), m2=None, S2=np.array([]), symmetric=False) return _KL_norm(m1, S1, m2, S2) -def choleskyL(M, return_scale_free=False): +def choleskyL_corr(M): r""" Gets the Cholesky lower triangular matrix :math:`L` (defined as :math:`M=LL^T`) - of a given matrix ``M``. + for the matrix ``M``, in the form :math:`L = S L^\prime` where S is diagonal. Can be used to create an affine transformation that decorrelates a sample - :math:`x=\{x_i\}` as :math:`y=Lx`, where :math:`L` is extracted from the - covariance of the sample. + :math:`x=\{x_i\}` with covariance M, as :math:`x=Ly`, + where :math:`L` is extracted from M and y has identity covariance. - If ``return_scale_free=True`` (default: ``False``), returns a tuple of - a matrix :math:`S` containing the square roots of the diagonal of the input matrix - (the standard deviations, if a covariance is given), and the scale-free - :math:`L^\prime=S^{-1}L`. + Returns a tuple of a matrix :math:`S` containing the square roots of the diagonal + of the input matrix (the standard deviations, if a covariance is given), + and the scale-free :math:`L^\prime=S^{-1}L`. + (could just use Cholesky directly for proposal) """ std_diag, corr = cov_to_std_and_corr(M) - Lprime = np.linalg.cholesky(corr) - if return_scale_free: - return std_diag, Lprime - else: - return np.linalg.inv(std_diag).dot(Lprime) + return np.diag(std_diag), np.linalg.cholesky(corr) def cov_to_std_and_corr(cov): """ - Gets the standard deviations (as a diagonal matrix) + Gets the standard deviations (as a 1D array and the correlation matrix of a covariance matrix. """ - std_diag = np.diag(np.sqrt(np.diag(cov))) - invstd_diag = np.linalg.inv(std_diag) - corr = invstd_diag.dot(cov).dot(invstd_diag) - return std_diag, corr + std = np.sqrt(np.diag(cov)) + inv_std = 1 / std + corr = inv_std[:, np.newaxis] * cov * inv_std[np.newaxis, :] + np.fill_diagonal(corr, 1.0) + return std, corr def are_different_params_lists(list_A, list_B, name_A="A", name_B="B"): diff --git a/cobaya/yaml.py b/cobaya/yaml.py index 43564b209..e137ad981 100644 --- a/cobaya/yaml.py +++ b/cobaya/yaml.py @@ -121,6 +121,7 @@ def path_constructor(loader, node): return (env_val or '') + value[match.end():] + DefaultsLoader.add_implicit_resolver('!path', path_matcher, None) DefaultsLoader.add_constructor('!path', path_constructor) @@ -218,7 +219,7 @@ def _numpy_array_representer(dumper, data): CustomDumper.add_representer(np.ndarray, _numpy_array_representer) def _numpy_int_representer(dumper, data): - return dumper.represent_int(data) + return dumper.represent_int(int(data)) CustomDumper.add_representer(np.int64, _numpy_int_representer) diff --git a/docs/cluster_amazon.rst b/docs/cluster_amazon.rst index cb469b8c4..f134fc446 100644 --- a/docs/cluster_amazon.rst +++ b/docs/cluster_amazon.rst @@ -22,7 +22,7 @@ Now install the requisites with $ bash miniconda.sh -b -p $HOME/miniconda $ export PATH="$HOME/miniconda/bin:$PATH" $ conda config --set always_yes yes --set changeps1 no - $ conda create -q -n cobaya-env python=3.9 scipy matplotlib cython PyYAML pytest pytest-forked flaky + $ conda create -q -n cobaya-env python=3.9 scipy matplotlib cython PyYAML pytest pytest-xdist flaky $ source activate cobaya-env $ pip install mpi4py diff --git a/pyproject.toml b/pyproject.toml index dbf230dfa..891b216fc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -36,7 +36,8 @@ dependencies = [ ] [project.optional-dependencies] -test = ["pytest", "pytest-forked", "flaky", "mpi4py", "iminuit"] +test = ["pytest", "pytest-xdist", "flaky", "mpi4py", "iminuit"] +speed = ["numba", "mpi4py"] gui = ["pyside6", "matplotlib"] docs = [ "sphinx", "sphinx_rtd_theme>=1.0", "sphinxcontrib-jquery", diff --git a/tests-environment.yml b/tests-environment.yml new file mode 100644 index 000000000..dda785a81 --- /dev/null +++ b/tests-environment.yml @@ -0,0 +1,24 @@ +name: tests-environment +channels: + - defaults + - conda-forge +dependencies: + - scipy + - matplotlib + - cython + - pyyaml + - dill + - coverage + - pytest + - pandas + - iminuit + - mpi4py + - openmpi + - pip + - pip: + - flake8 + - flaky + - pytest-xdist + - pytest-cov + - camb + - -r requirements.txt diff --git a/tests/test_mcmc.py b/tests/test_mcmc.py index ffb8b73cd..3eebcdbab 100644 --- a/tests/test_mcmc.py +++ b/tests/test_mcmc.py @@ -236,7 +236,7 @@ def _test_overhead_timing(dim=15): from cProfile import Profile from io import StringIO # noinspection PyUnresolvedReferences - from cobaya.samplers.mcmc import proposal # one-time numba compile out of profiling + from cobaya.functions import random_SO_N # one-time numba compile out of profiling like_test = _make_gaussian_like(dim) info: InputDict = {'likelihood': {'like': like_test}, 'debug': False,