diff --git a/.github/workflows/main.yml b/.github/workflows/main.yml index 7ef718bba..4fb654d4f 100644 --- a/.github/workflows/main.yml +++ b/.github/workflows/main.yml @@ -7,7 +7,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: [3.8, 3.9, '3.10', '3.11'] + python-version: [3.9, '3.10', '3.11', '3.12'] steps: - name: Checkout PyAutoConf uses: actions/checkout@v2 diff --git a/.gitignore b/.gitignore index 804d51ed2..c7e026207 100644 --- a/.gitignore +++ b/.gitignore @@ -125,13 +125,13 @@ venv.bak/ .idea workspace/output/ output -test/optimize/test_fit +test/mle/test_fit test/test_files/text/ test/ -test_autofit/optimize/test_fit/ +test_autofit/mle/test_fit/ test_autofit/test_files/text/psycopg2-binary==2.8.1 test_autofit/test_files/text/ -fit/test_autofit/optimize/test_fit +fit/test_autofit/mle/test_fit *.DS_Store test_autofit/config/priors/old @@ -157,7 +157,7 @@ test_autofit/samples.csv __MACOSX *.swp test/autofit/test_fit -# Byte-compiled / optimized / DLL files +# Byte-compiled / mled / DLL files __pycache__/ *.py[cod] *$py.class @@ -264,13 +264,13 @@ venv.bak/ .idea workspace/output/ output -test/optimize/test_fit +test/mle/test_fit test/test_files/text/ test/ -test_autofit/optimize/test_fit/ +test_autofit/mle/test_fit/ test_autofit/test_files/text/psycopg2-binary==2.8.1 test_autofit/test_files/text/ -fit/test_autofit/optimize/test_fit +fit/test_autofit/mle/test_fit *.DS_Store test_autofit/config/priors/old diff --git a/README.rst b/README.rst index c2bf20e40..79379c1d7 100644 --- a/README.rst +++ b/README.rst @@ -45,7 +45,7 @@ The following links are useful for new starters: - `The introduction Jupyter Notebook on Binder `_, where you can try **PyAutoFit** in a web browser (without installation). -- `The autofit_workspace GitHub repository `_, which includes example scripts and the `HowToFit Jupyter notebook lectures `_ which give new users a step-by-step introduction to **PyAutoFit**. +- `The autofit_workspace GitHub repository `_, which includes example scripts and the `HowToFit Jupyter notebook lectures `_ which give new users a step-by-step introduction to **PyAutoFit**. Support ------- @@ -72,7 +72,7 @@ API Overview To illustrate the **PyAutoFit** API, we use an illustrative toy model of fitting a one-dimensional Gaussian to noisy 1D data. Here's the ``data`` (black) and the model (red) we'll fit: -.. image:: https://raw.githubusercontent.com/rhayes777/PyAutoFit/master/files/toy_model_fit.png +.. image:: https://raw.githubusercontent.com/rhayes777/PyAutoFit/main/files/toy_model_fit.png :width: 400 We define our model, a 1D Gaussian by writing a Python class using the format below: diff --git a/autofit/__init__.py b/autofit/__init__.py index 7e055796c..c243f6edb 100644 --- a/autofit/__init__.py +++ b/autofit/__init__.py @@ -57,7 +57,7 @@ from .mapper.prior_model.annotation import AnnotationPriorModel from .mapper.prior_model.collection import Collection from .mapper.prior_model.prior_model import Model -from .mapper.prior_model.prior_model import Model +from .mapper.prior_model.array import Array from .non_linear.search.abstract_search import NonLinearSearch from .non_linear.analysis.visualize import Visualizer from .non_linear.analysis.analysis import Analysis @@ -66,7 +66,8 @@ from .non_linear.grid.sensitivity import Sensitivity from .non_linear.initializer import InitializerBall from .non_linear.initializer import InitializerPrior -from .non_linear.initializer import SpecificRangeInitializer +from .non_linear.initializer import InitializerParamBounds +from .non_linear.initializer import InitializerParamStartPoints from .non_linear.search.mcmc.auto_correlations import AutoCorrelationsSettings from .non_linear.search.mcmc.emcee.search import Emcee from .non_linear.search.mcmc.zeus.search import Zeus @@ -74,10 +75,11 @@ from .non_linear.search.nest.dynesty.search.dynamic import DynestyDynamic from .non_linear.search.nest.dynesty.search.static import DynestyStatic from .non_linear.search.nest.ultranest.search import UltraNest -from .non_linear.search.optimize.drawer.search import Drawer -from .non_linear.search.optimize.lbfgs.search import LBFGS -from .non_linear.search.optimize.pyswarms.search.globe import PySwarmsGlobal -from .non_linear.search.optimize.pyswarms.search.local import PySwarmsLocal +from .non_linear.search.mle.drawer.search import Drawer +from .non_linear.search.mle.bfgs.search import BFGS +from .non_linear.search.mle.bfgs.search import LBFGS +from .non_linear.search.mle.pyswarms.search.globe import PySwarmsGlobal +from .non_linear.search.mle.pyswarms.search.local import PySwarmsLocal from .non_linear.paths.abstract import AbstractPaths from .non_linear.paths import DirectoryPaths from .non_linear.paths import DatabasePaths @@ -132,4 +134,4 @@ def save_abc(pickler, obj): conf.instance.register(__file__) -__version__ = "2024.07.16.1" +__version__ = "2024.9.21.2" diff --git a/autofit/config/non_linear/README.rst b/autofit/config/non_linear/README.rst index 11774c406..9432dd852 100644 --- a/autofit/config/non_linear/README.rst +++ b/autofit/config/non_linear/README.rst @@ -6,4 +6,4 @@ Files - ``mcmc.yaml``: Settings default behaviour of MCMC non-linear searches (e.g. Emcee). - ``nest.yaml``: Settings default behaviour of nested sampler non-linear searches (e.g. Dynesty). -- ``optimizer.yaml``: Settings default behaviour of optimizer non-linear searches (e.g. PySwarms). \ No newline at end of file +- ``mle.yaml``: Settings default behaviour of maximum likelihood estimator (mle) searches (e.g. PySwarms). \ No newline at end of file diff --git a/autofit/config/non_linear/optimize.yaml b/autofit/config/non_linear/mle.yaml similarity index 66% rename from autofit/config/non_linear/optimize.yaml rename to autofit/config/non_linear/mle.yaml index 81a1640f5..3946b950d 100644 --- a/autofit/config/non_linear/optimize.yaml +++ b/autofit/config/non_linear/mle.yaml @@ -1,11 +1,11 @@ # Configuration files that customize the default behaviour of non-linear searches. -# **PyAutoFit** supports the following optimizer algorithms: +# **PyAutoFit** supports the following maximum likelihood estimator (MLE) algorithms: # - PySwarms: https://github.com/ljvmiranda921/pyswarms / https://pyswarms.readthedocs.io/en/latest/index.html # Settings in the [search], [run] and [options] entries are specific to each nested algorithm and should be -# determined by consulting that optimizers method's own readthedocs. +# determined by consulting that method's own readthedocs. PySwarmsGlobal: run: @@ -49,6 +49,44 @@ PySwarmsLocal: updates: iterations_per_update: 500 # The number of iterations of the non-linear search performed between every 'update', where an update performs tasks like outputting model.results. remove_state_files_at_end: true # Whether to remove the savestate of the seach (e.g. the Emcee hdf5 file) at the end to save hard-disk space (results are still stored as PyAutoFit pickles and loadable). +BFGS: + search: + tol: null + options: + disp: false + eps: 1.0e-08 + ftol: 2.220446049250313e-09 + gtol: 1.0e-05 + iprint: -1.0 + maxcor: 10 + maxfun: 15000 + maxiter: 15000 + maxls: 20 + initialize: # The method used to generate where walkers are initialized in parameter space {prior | ball}. + method: ball # priors: samples are initialized by randomly drawing from each parameter's prior. ball: samples are initialized by randomly drawing unit values from a narrow uniform distribution. + ball_lower_limit: 0.49 # The lower limit of the uniform distribution unit values are drawn from when initializing walkers using the ball method. + ball_upper_limit: 0.51 # The upper limit of the uniform distribution unit values are drawn from when initializing walkers using the ball method. + parallel: + number_of_cores: 1 # The number of cores the search is parallelized over by default, using Python multiprocessing. + printing: + silence: false # If True, the default print output of the non-linear search is silcened and not printed by the Python interpreter. + updates: + iterations_per_update: 500 # The number of iterations of the non-linear search performed between every 'update', where an update performs tasks like outputting model.results. + remove_state_files_at_end: true # Whether to remove the savestate of the seach (e.g. the Emcee hdf5 file) at the end to save hard-disk space (results are still stored as PyAutoFit pickles and loadable). +Drawer: + search: + total_draws: 50 + initialize: # The method used to generate where walkers are initialized in parameter space {prior | ball}. + method: ball # priors: samples are initialized by randomly drawing from each parameter's prior. ball: samples are initialized by randomly drawing unit values from a narrow uniform distribution. + ball_lower_limit: 0.49 # The lower limit of the uniform distribution unit values are drawn from when initializing walkers using the ball method. + ball_upper_limit: 0.51 # The upper limit of the uniform distribution unit values are drawn from when initializing walkers using the ball method. + parallel: + number_of_cores: 1 # The number of cores the search is parallelized over by default, using Python multiprocessing. + printing: + silence: false # If True, the default print output of the non-linear search is silcened and not printed by the Python interpreter. + updates: + iterations_per_update: 500 # The number of iterations of the non-linear search performed between every 'update', where an update performs tasks like outputting model.results. + remove_state_files_at_end: true # Whether to remove the savestate of the seach (e.g. the Emcee hdf5 file) at the end to save hard-disk space (results are still stored as PyAutoFit pickles and loadable). LBFGS: search: tol: null diff --git a/autofit/config/visualize/plots_search.yaml b/autofit/config/visualize/plots_search.yaml index b6e7331bb..7552865d3 100644 --- a/autofit/config/visualize/plots_search.yaml +++ b/autofit/config/visualize/plots_search.yaml @@ -2,5 +2,6 @@ nest: corner_anesthetic: true # Output corner figure (using anestetic) during a non-linear search fit? mcmc: corner_cornerpy: true # Output corner figure (using corner.py) during a non-linear search fit? -optimize: - corner_cornerpy: true # Output corner figure (using corner.py) during a non-linear search fit? \ No newline at end of file +mle: + subplot_parameters: true # Output a subplot of the best-fit parameters of the model? + log_likelihood_vs_iteration: true # Output a plot of the log likelihood versus iteration number? \ No newline at end of file diff --git a/autofit/database/aggregator/aggregator.py b/autofit/database/aggregator/aggregator.py index 4bfbebc9c..4c27f0c3d 100644 --- a/autofit/database/aggregator/aggregator.py +++ b/autofit/database/aggregator/aggregator.py @@ -1,5 +1,6 @@ import logging from abc import ABC, abstractmethod +from sqlalchemy import text from typing import Optional, List, Union, cast from ..sqlalchemy_ import sa @@ -370,7 +371,7 @@ def _fits_for_query(self, query: str) -> List[m.Fit]: query """ logger.debug(f"Executing query: {query}") - fit_ids = {row[0] for row in self.session.execute(query)} + fit_ids = {row[0] for row in self.session.execute(text(query))} logger.info(f"{len(fit_ids)} fit(s) found matching query") query = self.session.query(m.Fit).filter(m.Fit.id.in_(fit_ids)) diff --git a/autofit/database/migration/migration.py b/autofit/database/migration/migration.py index f540c43e1..c77c96515 100644 --- a/autofit/database/migration/migration.py +++ b/autofit/database/migration/migration.py @@ -1,14 +1,13 @@ import logging from abc import ABC, abstractmethod from hashlib import md5 +from sqlalchemy import text from typing import Union, Generator, Iterable, Optional from .session_wrapper import SessionWrapper from ..sqlalchemy_ import sa -logger = logging.getLogger( - __name__ -) +logger = logging.getLogger(__name__) class Identifiable(ABC): @@ -19,22 +18,13 @@ def id(self) -> str: A unique identifier generated by hashing a string """ - def __eq__( - self, - other: Union["Identifiable", str] - ) -> bool: + def __eq__(self, other: Union["Identifiable", str]) -> bool: """ Compares ids """ - if isinstance( - other, - Identifiable - ): + if isinstance(other, Identifiable): return self.id == other.id - if isinstance( - other, - str - ): + if isinstance(other, str): return self.id == other return False @@ -57,13 +47,7 @@ def id(self) -> str: """ Hash generated from underlying SQL statements """ - return md5( - ":".join( - self.strings - ).encode( - "utf-8" - ) - ).hexdigest() + return md5(":".join(self.strings).encode("utf-8")).hexdigest() def __str__(self): return "\n".join(self.strings) @@ -72,10 +56,7 @@ def __str__(self): class Revision(Identifiable): - def __init__( - self, - steps: Iterable[Step] - ): + def __init__(self, steps: Iterable[Step]): """ A specific revision of the database. This comprises a set of sequential steps and is uniquely identified @@ -95,12 +76,7 @@ def id(self) -> str: A unique identifier created by joining and hashing the identifiers of comprised steps. """ - return md5( - ":".join( - step.id for step - in self.steps - ).encode("utf-8") - ).hexdigest() + return md5(":".join(step.id for step in self.steps).encode("utf-8")).hexdigest() def __sub__(self, other: "Revision") -> "Revision": """ @@ -121,17 +97,11 @@ def __sub__(self, other: "Revision") -> "Revision": An object comprising steps required to move from the other revision to this revision. """ - return Revision(tuple( - step for step in self.steps - if step not in other.steps - )) + return Revision(tuple(step for step in self.steps if step not in other.steps)) class Migrator: - def __init__( - self, - *steps: Step - ): + def __init__(self, *steps: Step): """ Manages migration of an old database. @@ -153,14 +123,9 @@ def revisions(self) -> Generator[Revision, None, None]: starting on the first step and terminating on any step """ for i in range(1, len(self._steps) + 1): - yield Revision( - self._steps[:i] - ) + yield Revision(self._steps[:i]) - def get_steps( - self, - revision_id: Optional[str] = None - ) -> Iterable[Step]: + def get_steps(self, revision_id: Optional[str] = None) -> Iterable[Step]: """ Retrieve steps required to go from the specified revision to the latest revision. @@ -188,9 +153,7 @@ def latest_revision(self) -> Revision: The latest revision according to the steps passed to the Migrator """ - return Revision( - self._steps - ) + return Revision(self._steps) def migrate(self, session: sa.orm.Session): """ @@ -207,19 +170,11 @@ def migrate(self, session: sa.orm.Session): session A session pointing at some database. """ - wrapper = SessionWrapper( - session - ) + wrapper = SessionWrapper(session) revision_id = wrapper.revision_id - steps = list( - self.get_steps( - revision_id - ) - ) + steps = list(self.get_steps(revision_id)) if len(steps) == 0: - logger.info( - "Database already at latest revision" - ) + logger.info("Database already at latest revision") return latest_revision_id = self.latest_revision.id @@ -230,14 +185,10 @@ def migrate(self, session: sa.orm.Session): for step in steps: for string in step.strings: try: - session.execute( - string - ) + session.execute(text(string)) except sa.exc.OperationalError as e: logger.debug(e) wrapper.revision_id = self.latest_revision.id - logger.info( - f"revision_id updated to {wrapper.revision_id}" - ) + logger.info(f"revision_id updated to {wrapper.revision_id}") diff --git a/autofit/database/migration/session_wrapper.py b/autofit/database/migration/session_wrapper.py index 3f61071a8..3f2d52cda 100644 --- a/autofit/database/migration/session_wrapper.py +++ b/autofit/database/migration/session_wrapper.py @@ -1,16 +1,13 @@ import logging from functools import wraps +from sqlalchemy import text from typing import Optional from ..sqlalchemy_ import sa -logger = logging.getLogger( - __name__ -) +logger = logging.getLogger(__name__) -def needs_revision_table( - func -): +def needs_revision_table(func): """ Applies to functions that depend on the existence of the revision table. If the table does not exist @@ -59,11 +56,9 @@ def _init_revision_table(self): Creates the revision table with a single null entry """ self.session.execute( - "CREATE TABLE revision (revision_id VARCHAR PRIMARY KEY)" - ) - self.session.execute( - "INSERT INTO revision (revision_id) VALUES (null)" + text("CREATE TABLE revision (revision_id VARCHAR PRIMARY KEY)") ) + self.session.execute(text("INSERT INTO revision (revision_id) VALUES (null)")) @property def is_table(self) -> bool: @@ -71,9 +66,7 @@ def is_table(self) -> bool: Does the revision table exist? """ try: - self.session.execute( - "SELECT 1 FROM revision" - ) + self.session.execute(text("SELECT 1 FROM revision")) return True except sa.exc.OperationalError: return False @@ -85,9 +78,7 @@ def revision_id(self) -> Optional[str]: Describes the current revision of the database. None if no revisions have been made. """ - for row in self.session.execute( - "SELECT revision_id FROM revision" - ): + for row in self.session.execute(text("SELECT revision_id FROM revision")): return row[0] return None @@ -95,8 +86,6 @@ def revision_id(self) -> Optional[str]: @needs_revision_table def revision_id(self, revision_id: str): self.session.execute( - sa.text( - f"UPDATE revision SET revision_id = :revision_id" - ), - {"revision_id": revision_id} + text(f"UPDATE revision SET revision_id = :revision_id"), + {"revision_id": revision_id}, ) diff --git a/autofit/database/model/array.py b/autofit/database/model/array.py index 11bc29ddb..138901e92 100644 --- a/autofit/database/model/array.py +++ b/autofit/database/model/array.py @@ -29,7 +29,12 @@ def __init__(self, **kwargs): _shape = sa.Column(sa.String) fit_id = sa.Column(sa.String, sa.ForeignKey("fit.id")) - fit = sa.orm.relationship("Fit", uselist=False, foreign_keys=[fit_id]) + fit = sa.orm.relationship( + "Fit", + uselist=False, + foreign_keys=[fit_id], + back_populates="arrays", + ) @property def shape(self): @@ -86,6 +91,13 @@ class HDU(Array): "polymorphic_identity": "hdu", } + fit = sa.orm.relationship( + "Fit", + uselist=False, + foreign_keys=[Array.fit_id], + back_populates="hdus", + ) + @property def header(self): """ diff --git a/autofit/database/model/common.py b/autofit/database/model/common.py index 75b71b741..da7f2dc15 100644 --- a/autofit/database/model/common.py +++ b/autofit/database/model/common.py @@ -1,11 +1,7 @@ -from typing import Union +from typing import Union, ItemsView, Any, Iterable, Tuple from ..sqlalchemy_ import sa -from autofit.mapper.prior import abstract -from autofit.mapper.prior_model import prior_model -from autofit.mapper.prior_model import collection - from .model import Object @@ -26,7 +22,15 @@ class Dict(Object): __mapper_args__ = {"polymorphic_identity": "dict"} def __call__(self): - return {child.name: child() for child in self.children} + d = {} + for child in self.children: + instance = child() + if child.name != "": + d[child.name] = instance + else: + d[instance[0]] = instance[1] + + return d @classmethod def _from_object(cls, source: dict): @@ -34,3 +38,22 @@ def _from_object(cls, source: dict): instance._add_children(source.items()) instance.cls = dict return instance + + def _add_children( + self, items: Union[ItemsView[str, Any], Iterable[Tuple[str, Any]]] + ): + """ + Add database representations of child attributes + + Parameters + ---------- + items + Attributes such as floats or priors that are associated + with the real object + """ + self.children = [ + Object.from_object(value, name=key) + if isinstance(key, str) + else Object.from_object((key, value)) + for key, value in items + ] diff --git a/autofit/database/model/fit.py b/autofit/database/model/fit.py index e1cf0e78f..30b493a4c 100644 --- a/autofit/database/model/fit.py +++ b/autofit/database/model/fit.py @@ -1,6 +1,7 @@ import json import pickle from functools import wraps +from sqlalchemy.orm import Mapped from typing import List import numpy as np @@ -29,7 +30,11 @@ def __init__(self, **kwargs): name = sa.Column(sa.String) string = sa.Column(sa.String) fit_id = sa.Column(sa.String, sa.ForeignKey("fit.id")) - fit = sa.orm.relationship("Fit", uselist=False) + fit = sa.orm.relationship( + "Fit", + uselist=False, + back_populates="pickles", + ) @property def value(self): @@ -63,7 +68,11 @@ def __init__(self, **kwargs): name = sa.Column(sa.String) string = sa.Column(sa.String) fit_id = sa.Column(sa.String, sa.ForeignKey("fit.id")) - fit = sa.orm.relationship("Fit", uselist=False) + fit = sa.orm.relationship( + "Fit", + uselist=False, + back_populates="jsons", + ) @property def dict(self): @@ -110,7 +119,10 @@ class NamedInstance(Base): instance_id = sa.Column(sa.Integer, sa.ForeignKey("object.id")) __instance = sa.orm.relationship( - "Object", uselist=False, backref="named_instance", foreign_keys=[instance_id] + "Object", + uselist=False, + backref="named_instance", + foreign_keys=[instance_id], ) @property @@ -181,7 +193,10 @@ class Fit(Base): ) is_complete = sa.Column(sa.Boolean) - _named_instances: List[NamedInstance] = sa.orm.relationship("NamedInstance") + _named_instances: Mapped[List[NamedInstance]] = sa.orm.relationship( + "NamedInstance", + back_populates="fit", + ) @property @try_none @@ -203,7 +218,7 @@ def named_instances(self): def total_parameters(self): return self.model.prior_count if self.model else 0 - _info: List[Info] = sa.orm.relationship("Info") + _info: Mapped[List[Info]] = sa.orm.relationship("Info", back_populates="fit") def __init__(self, **kwargs): try: @@ -216,7 +231,7 @@ def __init__(self, **kwargs): parent_id = sa.Column(sa.String, sa.ForeignKey("fit.id")) - children: List["Fit"] = sa.orm.relationship( + children: Mapped[List["Fit"]] = sa.orm.relationship( "Fit", backref=sa.orm.backref("parent", remote_side=[id]) ) @@ -301,14 +316,22 @@ def model(self) -> AbstractPriorModel: def model(self, model: AbstractPriorModel): self.__model = Object.from_object(model) - pickles: List[Pickle] = sa.orm.relationship("Pickle", lazy="joined") - jsons: List[JSON] = sa.orm.relationship("JSON", lazy="joined") - arrays: List[Array] = sa.orm.relationship( + pickles: Mapped[List[Pickle]] = sa.orm.relationship( + "Pickle", + lazy="joined", + foreign_keys=[Pickle.fit_id], + ) + jsons: Mapped[List[JSON]] = sa.orm.relationship( + "JSON", + lazy="joined", + foreign_keys=[JSON.fit_id], + ) + arrays: Mapped[List[Array]] = sa.orm.relationship( "Array", lazy="joined", foreign_keys=[Array.fit_id], ) - hdus: List[HDU] = sa.orm.relationship( + hdus: Mapped[List[HDU]] = sa.orm.relationship( "HDU", lazy="joined", foreign_keys=[HDU.fit_id], @@ -332,7 +355,10 @@ def __getitem__(self, item: str): """ for p in self.jsons + self.arrays + self.hdus + self.pickles: if p.name == item: - return p.value + value = p.value + if item == "samples_summary": + value.model = self.model + return value return getattr(self, item) diff --git a/autofit/database/model/model.py b/autofit/database/model/model.py index 89e311335..a52b87228 100644 --- a/autofit/database/model/model.py +++ b/autofit/database/model/model.py @@ -1,13 +1,15 @@ import abc import inspect +from sqlalchemy.orm import Mapped from typing import List, Tuple, Any, Iterable, Union, ItemsView, Type import numpy as np from autoconf.class_path import get_class, get_class_path -from ..sqlalchemy_ import sa, declarative +from ..sqlalchemy_ import sa +from sqlalchemy.orm import declarative_base -Base = declarative.declarative_base() +Base = declarative_base() _schema_version = 1 @@ -32,7 +34,10 @@ class Object(Base): samples_for_id = sa.Column(sa.String, sa.ForeignKey("fit.id")) samples_for = sa.orm.relationship( - "Fit", uselist=False, foreign_keys=[samples_for_id] + "Fit", + uselist=False, + foreign_keys=[samples_for_id], + back_populates="_samples", ) latent_samples_for_id = sa.Column(sa.String, sa.ForeignKey("fit.id")) @@ -40,11 +45,13 @@ class Object(Base): "Fit", uselist=False, foreign_keys=[latent_samples_for_id], + back_populates="_latent_samples", ) - children: List["Object"] = sa.orm.relationship( + children: Mapped[List["Object"]] = sa.orm.relationship( "Object", uselist=True, + back_populates="parent", ) def __len__(self): diff --git a/autofit/example/analysis.py b/autofit/example/analysis.py index 836027d8c..581415228 100644 --- a/autofit/example/analysis.py +++ b/autofit/example/analysis.py @@ -116,7 +116,7 @@ def model_data_1d_from(self, instance: af.ModelInstance) -> np.ndarray: def save_attributes(self, paths: af.DirectoryPaths): """ Before the model-fit via the non-linear search begins, this routine saves attributes of the `Analysis` object - to the `pickles` folder such that they can be loaded after the analysis using PyAutoFit's database and + to the `files` folder such that they can be loaded after the analysis using PyAutoFit's database and aggregator tools. For this analysis the following are output: @@ -132,7 +132,7 @@ def save_attributes(self, paths: af.DirectoryPaths): Parameters ---------- paths - The PyAutoFit paths object which manages all paths, e.g. where the non-linear search outputs are stored, + The paths object which manages all paths, e.g. where the non-linear search outputs are stored, visualization, and the pickled objects used by the aggregator output by this function. """ paths.save_json(name="data", object_dict=self.data.tolist(), prefix="dataset") diff --git a/autofit/example/visualize.py b/autofit/example/visualize.py index a66f49824..7493021d6 100644 --- a/autofit/example/visualize.py +++ b/autofit/example/visualize.py @@ -33,7 +33,7 @@ def visualize_before_fit( analysis The analysis class used to perform the model-fit whose quantities are being visualized. paths - The PyAutoFit paths object which manages all paths, e.g. where the non-linear search outputs are stored, + The paths object which manages all paths, e.g. where the non-linear search outputs are stored, visualization, and the pickled objects used by the aggregator output by this function. model The model which is fitted to the data, which may be used to customize the visualization. @@ -84,7 +84,7 @@ def visualize( analysis The analysis class used to perform the model-fit whose quantities are being visualized. paths - The PyAutoFit paths object which manages all paths, e.g. where the non-linear search outputs are stored, + The paths object which manages all paths, e.g. where the non-linear search outputs are stored, visualization, and the pickled objects used by the aggregator output by this function. instance An instance of the model that is being fitted to the data by this analysis (whose parameters have been set @@ -163,7 +163,7 @@ def visualize_before_fit_combined( analyses A list of the analysis classes used to perform the model-fit whose quantities are being visualized. paths - The PyAutoFit paths object which manages all paths, e.g. where the non-linear search outputs are stored, + The paths object which manages all paths, e.g. where the non-linear search outputs are stored, visualization, and the pickled objects used by the aggregator output by this function. model The model which is fitted to the data, which may be used to customize the visualization. @@ -202,7 +202,7 @@ def visualize_combined( analyses A list of the analysis classes used to perform the model-fit whose quantities are being visualized. paths - The PyAutoFit paths object which manages all paths, e.g. where the non-linear search outputs are stored, + The paths object which manages all paths, e.g. where the non-linear search outputs are stored, visualization, and the pickled objects used by the aggregator output by this function. model The model which is fitted to the data, which may be used to customize the visualization. diff --git a/autofit/graphical/laplace/line_search.py b/autofit/graphical/laplace/line_search.py index 4f3e3e0f4..dc56d5952 100644 --- a/autofit/graphical/laplace/line_search.py +++ b/autofit/graphical/laplace/line_search.py @@ -11,7 +11,7 @@ from typing import Optional, Dict, Tuple import numpy as np -from scipy.optimize import linesearch +from scipy.optimize import _linesearch as linesearch from autoconf import cached_property from autofit.graphical.factor_graphs.abstract import ( @@ -188,7 +188,6 @@ def _next_state(self, stepsize): return next_state def step(self, stepsize): - if not stepsize: return self diff --git a/autofit/interpolator/covariance.py b/autofit/interpolator/covariance.py index fd090d287..4e7d0dcbd 100644 --- a/autofit/interpolator/covariance.py +++ b/autofit/interpolator/covariance.py @@ -195,8 +195,8 @@ def _relationships_for_value( """ analysis = self._analysis_for_value(value) model = self.model(path_relationship_map=path_relationship_map or {}) - optimizer = DynestyStatic() - result = optimizer.fit(model=model, analysis=analysis) + search = DynestyStatic() + result = search.fit(model=model, analysis=analysis) return result.instance def __getitem__(self, value: Equality) -> float: diff --git a/autofit/interpolator/linear.py b/autofit/interpolator/linear.py index 131381981..a923da9f6 100644 --- a/autofit/interpolator/linear.py +++ b/autofit/interpolator/linear.py @@ -1,4 +1,4 @@ -from scipy.stats import stats +from scipy.stats import linregress from .abstract import AbstractInterpolator @@ -9,5 +9,5 @@ class LinearInterpolator(AbstractInterpolator): @staticmethod def _interpolate(x, y, value): - slope, intercept, r, p, std_err = stats.linregress(x, y) + slope, intercept, r, p, std_err = linregress(x, y) return slope * value + intercept diff --git a/autofit/mapper/mock/mock_model.py b/autofit/mapper/mock/mock_model.py index fe8e299ba..f5811973b 100644 --- a/autofit/mapper/mock/mock_model.py +++ b/autofit/mapper/mock/mock_model.py @@ -41,7 +41,7 @@ def __init__(self, one=1, two=2, three=3): class MockClassx2Tuple: def __init__(self, one_tuple=(0.0, 0.0)): """Abstract MockParent, describing an object with y, x cartesian - coordinates """ + coordinates""" self.one_tuple = one_tuple def __eq__(self, other): @@ -94,7 +94,6 @@ def __init__(self, tup=(0.0, 0.0)): class MockOverload: - def __init__(self, one=1.0): self.one = one @@ -112,11 +111,11 @@ def two(self, two): class MockComponents: def __init__( - self, - components_0: list = None, - components_1: list = None, - parameter=None, - **kwargs + self, + components_0: list = None, + components_1: list = None, + parameter=None, + **kwargs ): self.parameter = parameter self.group_0 = components_0 @@ -134,7 +133,7 @@ def __eq__(self, other): class MockChildTuple(MockParent): def __init__(self, tup=(0.0, 0.0)): - """ Generic circular profiles class to contain functions shared by light and + """Generic circular profiles class to contain functions shared by light and mass profiles. Parameters @@ -147,7 +146,7 @@ def __init__(self, tup=(0.0, 0.0)): class MockChildTuplex2(MockChildTuple): def __init__(self, tup=(0.0, 0.0), one=1.0, two=0.0): - """ Generic elliptical profiles class to contain functions shared by light + """Generic elliptical profiles class to contain functions shared by light and mass profiles. Parameters @@ -166,7 +165,7 @@ def __init__(self, tup=(0.0, 0.0), one=1.0, two=0.0): class MockChildTuplex3(MockChildTuple): def __init__(self, tup=(0.0, 0.0), one=1.0, two=0.0, three=0.0): - """ Generic elliptical profiles class to contain functions shared by light + """Generic elliptical profiles class to contain functions shared by light and mass profiles. Parameters @@ -182,3 +181,13 @@ def __init__(self, tup=(0.0, 0.0), one=1.0, two=0.0, three=0.0): self.one = one self.two = two self.three = three + + +class Parameter: + def __init__(self, value: float = 0.5): + self.value = value + + +class WithString: + def __init__(self, arg: "Parameter"): + self.arg = arg diff --git a/autofit/mapper/model_object.py b/autofit/mapper/model_object.py index aea48e872..05be46ef7 100644 --- a/autofit/mapper/model_object.py +++ b/autofit/mapper/model_object.py @@ -150,6 +150,7 @@ def from_dict( from autofit.mapper.prior_model.collection import Collection from autofit.mapper.prior_model.prior_model import Model from autofit.mapper.prior.abstract import Prior + from autofit.mapper.prior.gaussian import GaussianPrior from autofit.mapper.prior.tuple_prior import TuplePrior from autofit.mapper.prior.arithmetic.compound import Compound from autofit.mapper.prior.arithmetic.compound import ModifiedPrior @@ -234,7 +235,10 @@ def get_class_path(): f"Could not find type for class path {class_path}. Defaulting to Instance placeholder." ) instance = ModelInstance() + elif type_ == "array": + from autofit.mapper.prior_model.array import Array + return Array.from_dict(d) else: try: return Prior.from_dict(d, loaded_ids=loaded_ids) @@ -276,6 +280,7 @@ def dict(self) -> dict: from autofit.mapper.prior_model.collection import Collection from autofit.mapper.prior_model.prior_model import Model from autofit.mapper.prior.tuple_prior import TuplePrior + from autofit.mapper.prior_model.array import Array if isinstance(self, Collection): type_ = "collection" @@ -285,6 +290,8 @@ def dict(self) -> dict: type_ = "model" elif isinstance(self, TuplePrior): type_ = "tuple_prior" + elif isinstance(self, Array): + type_ = "array" else: raise AssertionError( f"{self.__class__.__name__} cannot be serialised to dict" diff --git a/autofit/mapper/prior/abstract.py b/autofit/mapper/prior/abstract.py index 260ef5928..b5f2cd3f2 100644 --- a/autofit/mapper/prior/abstract.py +++ b/autofit/mapper/prior/abstract.py @@ -1,5 +1,4 @@ import itertools -import os import random from abc import ABC, abstractmethod from copy import copy @@ -96,16 +95,11 @@ def new(self): new.id = next(self._ids) return new + @abstractmethod def with_limits(self, lower_limit: float, upper_limit: float) -> "Prior": """ Create a new instance of the same prior class with the passed limits. """ - new = self.__class__( - lower_limit=max(lower_limit, self.lower_limit), - upper_limit=min(upper_limit, self.upper_limit), - ) - new.message = self.message - return new @property def factor(self): diff --git a/autofit/mapper/prior/arithmetic/compound.py b/autofit/mapper/prior/arithmetic/compound.py index 98952fc74..0355c36fb 100644 --- a/autofit/mapper/prior/arithmetic/compound.py +++ b/autofit/mapper/prior/arithmetic/compound.py @@ -17,7 +17,7 @@ def retrieve_name(var): first_name = None frame = inspect.currentframe() while frame is not None: - for name, value in frame.f_locals.items(): + for name, value in list(frame.f_locals.items()): if var is value: first_name = name frame = frame.f_back diff --git a/autofit/mapper/prior/log_gaussian.py b/autofit/mapper/prior/log_gaussian.py index d13f9d203..a02d77e7b 100644 --- a/autofit/mapper/prior/log_gaussian.py +++ b/autofit/mapper/prior/log_gaussian.py @@ -71,6 +71,37 @@ def __init__( id_=id_, ) + @classmethod + def with_limits(cls, lower_limit: float, upper_limit: float) -> "LogGaussianPrior": + """ + Create a new gaussian prior centred between two limits + with sigma distance between this limits. + + Note that these limits are not strict so exceptions will not + be raised for values outside of the limits. + + This function is typically used in prior passing, where the + result of a model-fit are used to create new Gaussian priors + centred on the previously estimated median PDF model. + + Parameters + ---------- + lower_limit + The lower limit of the new Gaussian prior. + upper_limit + The upper limit of the new Gaussian Prior. + + Returns + ------- + A new GaussianPrior + """ + return cls( + mean=(lower_limit + upper_limit) / 2, + sigma=upper_limit - lower_limit, + lower_limit=lower_limit, + upper_limit=upper_limit, + ) + def _new_for_base_message(self, message): """ Create a new instance of this wrapper but change the parameters used diff --git a/autofit/mapper/prior/uniform.py b/autofit/mapper/prior/uniform.py index 3aa799396..3cea04a90 100644 --- a/autofit/mapper/prior/uniform.py +++ b/autofit/mapper/prior/uniform.py @@ -64,6 +64,16 @@ def __init__( def tree_flatten(self): return (self.lower_limit, self.upper_limit), (self.id,) + def with_limits( + self, + lower_limit: float, + upper_limit: float, + ) -> "Prior": + return UniformPrior( + lower_limit=lower_limit, + upper_limit=upper_limit, + ) + def logpdf(self, x): # TODO: handle x as a numpy array if x == self.lower_limit: @@ -97,9 +107,9 @@ def value_for(self, unit: float, ignore_prior_limits: bool = False) -> float: physical_value = prior.value_for(unit=0.2) """ - return float(round( - super().value_for(unit, ignore_prior_limits=ignore_prior_limits), 14 - )) + return float( + round(super().value_for(unit, ignore_prior_limits=ignore_prior_limits), 14) + ) def log_prior_from_value(self, value): """ diff --git a/autofit/mapper/prior_model/abstract.py b/autofit/mapper/prior_model/abstract.py index 04695dc68..70acc8548 100644 --- a/autofit/mapper/prior_model/abstract.py +++ b/autofit/mapper/prior_model/abstract.py @@ -1334,7 +1334,9 @@ def _explode_path(path_): raise AssertionError(f"No path was found matching {name}") def instance_from_prior_name_arguments( - self, prior_name_arguments: Dict[str, float] + self, + prior_name_arguments: Dict[str, float], + ignore_assertions: bool = False, ): """ Instantiate the model from the names of priors and @@ -1346,6 +1348,8 @@ def instance_from_prior_name_arguments( The names of priors where names of models and the name of the prior have been joined by underscores, mapped to corresponding values. + ignore_assertions + If True, assertions will not be checked Returns ------- @@ -1355,10 +1359,15 @@ def instance_from_prior_name_arguments( { self.path_for_name(name): value for name, value in prior_name_arguments.items() - } + }, + ignore_assertions=ignore_assertions, ) - def instance_from_path_arguments(self, path_arguments: Dict[Tuple[str], float]): + def instance_from_path_arguments( + self, + path_arguments: Dict[Tuple[str], float], + ignore_assertions: bool = False, + ): """ Create an instance from a dictionary mapping paths to tuples to corresponding values. @@ -1370,6 +1379,8 @@ def instance_from_path_arguments(self, path_arguments: Dict[Tuple[str], float]): Note that, for linked priors, each path only needs to be specified once. If multiple paths for the same prior are specified then the value for the last path will be used. + ignore_assertions + If True, assertions will not be checked Returns ------- @@ -1378,7 +1389,10 @@ def instance_from_path_arguments(self, path_arguments: Dict[Tuple[str], float]): arguments = { self.object_for_path(path): value for path, value in path_arguments.items() } - return self._instance_for_arguments(arguments) + return self._instance_for_arguments( + arguments, + ignore_assertions=ignore_assertions, + ) @property def prior_count(self) -> int: @@ -1465,6 +1479,14 @@ def paths(self) -> List[Path]: """ return [path for path, _ in self.path_priors_tuples] + @property + def paths_formatted(self) -> List[Path]: + """ + A list of paths to all the priors in the model, ordered by their + ids + """ + return [path for path, _ in self.path_priors_tuples] + @property def composition(self): return [".".join(path) for path in self.paths] diff --git a/autofit/mapper/prior_model/array.py b/autofit/mapper/prior_model/array.py new file mode 100644 index 000000000..c37c786b2 --- /dev/null +++ b/autofit/mapper/prior_model/array.py @@ -0,0 +1,215 @@ +from typing import Tuple, Dict, Optional, Union + +from autoconf.dictable import from_dict +from .abstract import AbstractPriorModel +from autofit.mapper.prior.abstract import Prior +import numpy as np + +from autofit.jax_wrapper import register_pytree_node_class + + +@register_pytree_node_class +class Array(AbstractPriorModel): + def __init__( + self, + shape: Tuple[int, ...], + prior: Optional[Prior] = None, + ): + """ + An array of priors. + + Parameters + ---------- + shape : (int, int) + The shape of the array. + prior : Prior + The prior of every entry in the array. + """ + super().__init__() + self.shape = shape + self.indices = list(np.ndindex(*shape)) + + if prior is not None: + for index in self.indices: + self[index] = prior.new() + + @staticmethod + def _make_key(index: Tuple[int, ...]) -> str: + """ + Make a key for the prior. + + This is so an index (e.g. (1, 2)) can be used to access a + prior (e.g. prior_1_2). + + Parameters + ---------- + index + The index of an element in an array. + + Returns + ------- + The attribute name for the prior. + """ + if isinstance(index, int): + suffix = str(index) + else: + suffix = "_".join(map(str, index)) + return f"prior_{suffix}" + + def _instance_for_arguments( + self, + arguments: Dict[Prior, float], + ignore_assertions: bool = False, + ) -> np.ndarray: + """ + Create an array where the prior at each index is replaced with the + a concrete value. + + Parameters + ---------- + arguments + The arguments to replace the priors with. + ignore_assertions + Whether to ignore assertions in the priors. + + Returns + ------- + The array with the priors replaced. + """ + array = np.zeros(self.shape) + for index in self.indices: + value = self[index] + try: + value = value.instance_for_arguments( + arguments, + ignore_assertions, + ) + except AttributeError: + pass + + array[index] = value + return array + + def __setitem__( + self, + index: Union[int, Tuple[int, ...]], + value: Union[float, Prior], + ): + """ + Set the value at an index. + + Parameters + ---------- + index + The index of the prior. + value + The new value. + """ + setattr( + self, + self._make_key(index), + value, + ) + + def __getitem__( + self, + index: Union[int, Tuple[int, ...]], + ) -> Union[float, Prior]: + """ + Get the value at an index. + + Parameters + ---------- + index + The index of the value. + + Returns + ------- + The value at the index. + """ + return getattr( + self, + self._make_key(index), + ) + + @classmethod + def from_dict( + cls, + d, + reference: Optional[Dict[str, str]] = None, + loaded_ids: Optional[dict] = None, + ) -> "Array": + """ + Create an array from a dictionary. + + Parameters + ---------- + d + The dictionary. + reference + A dictionary of references. + loaded_ids + A dictionary of loaded ids. + + Returns + ------- + The array. + """ + arguments = d["arguments"] + shape = from_dict(arguments["shape"]) + array = cls(shape) + for key, value in arguments.items(): + if key.startswith("prior"): + setattr(array, key, from_dict(value)) + + return array + + def tree_flatten(self): + """ + Flatten this array model as a PyTree. + """ + members = [self[index] for index in self.indices] + return members, (self.shape,) + + @classmethod + def tree_unflatten(cls, aux_data, children): + """ + Unflatten a PyTree into an array model. + """ + (shape,) = aux_data + instance = cls(shape) + for index, child in zip(instance.indices, children): + instance[index] = child + + return instance + + @property + def prior_class_dict(self): + return { + **{ + prior: cls + for prior_model in self.direct_prior_model_tuples + for prior, cls in prior_model[1].prior_class_dict.items() + }, + **{prior: np.ndarray for _, prior in self.direct_prior_tuples}, + } + + def gaussian_prior_model_for_arguments(self, arguments: Dict[Prior, Prior]): + """ + Returns a new instance of model mapper with a set of Gaussian priors based on + tuples provided by a previous nonlinear search. + + Parameters + ---------- + arguments + Tuples providing the mean and sigma of gaussians + + Returns + ------- + A new model mapper populated with Gaussian priors + """ + new_array = Array(self.shape) + for index in self.indices: + new_array[index] = self[index].gaussian_prior_model_for_arguments(arguments) + + return new_array diff --git a/autofit/mapper/prior_model/prior_model.py b/autofit/mapper/prior_model/prior_model.py index 43f2dcc8c..ab2f99745 100644 --- a/autofit/mapper/prior_model/prior_model.py +++ b/autofit/mapper/prior_model/prior_model.py @@ -1,10 +1,9 @@ -import builtins import collections.abc import copy import inspect import logging -from typing import List import typing +from typing import * from autofit.jax_wrapper import register_pytree_node_class, register_pytree_node @@ -16,6 +15,7 @@ from autofit.mapper.prior.deferred import DeferredInstance from autofit.mapper.prior.tuple_prior import TuplePrior from autofit.mapper.prior_model.abstract import AbstractPriorModel +from autofit.mapper.prior_model.util import gather_namespaces from autofit.tools.namer import namer logger = logging.getLogger(__name__) @@ -112,19 +112,23 @@ def __init__( self.cls = cls - try: - annotations = inspect.getfullargspec(cls).annotations - for key, value in annotations.items(): - if isinstance(value, str): - annotations[key] = getattr(builtins, value) - except TypeError: - annotations = dict() + namespaces = gather_namespaces(cls) + + annotations = typing.get_type_hints( + cls.__init__, + namespaces, + namespaces, + ) try: arg_spec = inspect.getfullargspec(cls) defaults = dict( zip(arg_spec.args[-len(arg_spec.defaults) :], arg_spec.defaults) ) + defaults = { + key: value.default if hasattr(value, "default") else value + for key, value in defaults.items() + } except TypeError: defaults = {} @@ -170,12 +174,21 @@ def __init__( else: annotation = annotations[arg] - if ( - hasattr(annotation, "__origin__") - and issubclass( - annotation.__origin__, collections.abc.Collection + if isinstance(annotation, str): + continue + + if arg in defaults: + value = self._convert_value(defaults[arg]) + elif ( + ( + hasattr(annotation, "__origin__") + and issubclass( + annotation.__origin__, collections.abc.Collection + ) ) - ) or isinstance(annotation, collections.abc.Collection): + or isinstance(annotation, collections.abc.Collection) + or issubclass(annotation, collections.abc.Collection) + ): from autofit.mapper.prior_model.collection import Collection value = Collection() @@ -232,9 +245,17 @@ def instance_flatten(self, instance): """ Flatten an instance of this model as a PyTree. """ + attribute_names = [ + name + for name in self.direct_argument_names + if hasattr(instance, name) and name not in self.constructor_argument_names + ] return ( - [getattr(instance, name) for name in self.direct_argument_names], - None, + ( + [getattr(instance, name) for name in self.constructor_argument_names], + [getattr(instance, name) for name in attribute_names], + ), + (attribute_names,), ) def instance_unflatten(self, aux_data, children): @@ -250,7 +271,12 @@ def instance_unflatten(self, aux_data, children): ------- An instance of this model. """ - return self.cls(**dict(zip(self.direct_argument_names, children))) + constructor_arguments, other_arguments = children + attribute_names = aux_data[0] + instance = self.cls(*constructor_arguments) + for name, value in zip(attribute_names, other_arguments): + setattr(instance, name, value) + return instance def tree_flatten(self): """ @@ -279,7 +305,11 @@ def constructor_argument_names(self) -> List[str]: """ if self.cls not in class_args_dict: try: - class_args_dict[self.cls] = inspect.getfullargspec(self.cls).args[1:] + class_args_dict[self.cls] = [ + arg + for arg in inspect.getfullargspec(self.cls).args + if arg != "self" + ] except TypeError: class_args_dict[self.cls] = [] return class_args_dict[self.cls] @@ -320,6 +350,8 @@ def make_prior(self, attribute_name): If no configuration can be found """ cls = self.cls + if isinstance(cls, ConfigException): + return cls if not inspect.isclass(cls): # noinspection PyProtectedMember cls = inspect._findclass(cls) diff --git a/autofit/mapper/prior_model/util.py b/autofit/mapper/prior_model/util.py index b364ce9bf..822ce5035 100644 --- a/autofit/mapper/prior_model/util.py +++ b/autofit/mapper/prior_model/util.py @@ -1,3 +1,7 @@ +import inspect +from typing import Type, Dict +import typing + from autofit.mapper.prior_model.attribute_pair import AttributeNameValue @@ -5,3 +9,25 @@ class PriorModelNameValue(AttributeNameValue): @property def prior_model(self): return self.value + + +def gather_namespaces(cls: Type) -> Dict[str, Dict]: + """ + Recursively gather the globals and locals for a given class and its parent classes. + """ + namespaces = {} + + try: + for base in inspect.getmro(cls): + if base is object: + continue + + module = inspect.getmodule(base) + if module: + namespaces.update(vars(module)) + except AttributeError: + pass + + namespaces.update(vars(typing)) + + return namespaces diff --git a/autofit/mock.py b/autofit/mock.py index d29ae6aa2..98ddd262a 100644 --- a/autofit/mock.py +++ b/autofit/mock.py @@ -2,7 +2,7 @@ from autofit.non_linear.mock.mock_result import MockResult from autofit.non_linear.mock.mock_result import MockResultGrid from autofit.non_linear.mock.mock_search import MockSearch -from autofit.non_linear.mock.mock_search import MockOptimizer +from autofit.non_linear.mock.mock_search import MockMLE from autofit.non_linear.mock.mock_samples_summary import MockSamplesSummary from autofit.non_linear.mock.mock_samples import MockSamples from autofit.non_linear.mock.mock_samples import MockSamplesNest diff --git a/autofit/non_linear/analysis/analysis.py b/autofit/non_linear/analysis/analysis.py index 7d5b24a88..5e464bf56 100644 --- a/autofit/non_linear/analysis/analysis.py +++ b/autofit/non_linear/analysis/analysis.py @@ -1,3 +1,4 @@ +import inspect import logging from abc import ABC from typing import Optional, Dict @@ -9,8 +10,6 @@ from autofit.non_linear.result import Result from autofit.non_linear.samples.samples import Samples from autofit.non_linear.samples.sample import Sample -from autofit.mapper.prior_model.collection import Collection -from autofit.mapper.prior.gaussian import GaussianPrior from .visualize import Visualizer from ..samples.util import simple_model_for_kwargs @@ -42,6 +41,10 @@ def __getattr__(self, item: str): raise AttributeError(f"Analysis has no attribute {item}") def method(*args, **kwargs): + parameters = inspect.signature(_method).parameters + if "analyses" in parameters: + logger.debug(f"Skipping {item} as this is not a combined analysis") + return return _method(self, *args, **kwargs) return method @@ -69,7 +72,10 @@ def compute_latent_samples(self, samples: Samples) -> Optional[Samples]: log_prior=sample.log_prior, weight=sample.weight, kwargs=self.compute_latent_variables( - sample.instance_for_model(model) + sample.instance_for_model( + model, + ignore_assertions=True, + ) ), ) ) diff --git a/autofit/non_linear/analysis/combined.py b/autofit/non_linear/analysis/combined.py index af2f6fb2b..417b1f905 100644 --- a/autofit/non_linear/analysis/combined.py +++ b/autofit/non_linear/analysis/combined.py @@ -362,11 +362,15 @@ def make_result( paths=paths, samples=samples, search_internal=search_internal, - analysis=analysis + analysis=analysis, ) for analysis in self.analyses ] - return CombinedResult(child_results) + return CombinedResult( + child_results, + samples, + samples_summary, + ) def __len__(self): return len(self.analyses) diff --git a/autofit/non_linear/analysis/visualize.py b/autofit/non_linear/analysis/visualize.py index 7b0e74313..d0000c28a 100644 --- a/autofit/non_linear/analysis/visualize.py +++ b/autofit/non_linear/analysis/visualize.py @@ -32,7 +32,7 @@ def should_visualize( Parameters ---------- paths - The PyAutoFit paths object which manages all paths, e.g. where the non-linear search outputs are stored, + The paths object which manages all paths, e.g. where the non-linear search outputs are stored, visualization and the pickled objects used by the aggregator output by this function. Returns diff --git a/autofit/non_linear/fitness.py b/autofit/non_linear/fitness.py index e9857cf0a..0478de9d1 100644 --- a/autofit/non_linear/fitness.py +++ b/autofit/non_linear/fitness.py @@ -34,6 +34,7 @@ def __init__( fom_is_log_likelihood: bool = True, resample_figure_of_merit: float = -np.inf, convert_to_chi_squared: bool = False, + store_history: bool = False, ): """ Interfaces with any non-linear search to fit the model to the data and return a log likelihood via @@ -66,6 +67,11 @@ def __init__( instead. The appropriate value depends on the search, but is typically either `None`, `-np.inf` or `1.0e99`. All values indicate to the non-linear search that the model-fit should be resampled or ignored. + Many searches do not store the history of the parameters and log likelihood values, often to save + memory on large model-fits. However, this can be useful, for example to plot the results of a model-fit + versus iteration number. If the `store_history` bool is `True`, the parameters and log likelihoods are stored + in the `parameters_history_list` and `figure_of_merit_history_list` attribute of the fitness object. + Parameters ---------- analysis @@ -85,6 +91,8 @@ def __init__( convert_to_chi_squared If `True`, the figure of merit returned is the log likelihood multiplied by -2.0, such that it is a chi-squared value that is minimized. + store_history + If `True`, the parameters and log likelihood values of every model-fit are stored in lists. """ self.analysis = analysis @@ -93,11 +101,16 @@ def __init__( self.fom_is_log_likelihood = fom_is_log_likelihood self.resample_figure_of_merit = resample_figure_of_merit self.convert_to_chi_squared = convert_to_chi_squared + self.store_history = store_history + self._log_likelihood_function = None if self.paths is not None: self.check_log_likelihood(fitness=self) + self.parameters_history_list = [] + self.log_likelihood_history_list = [] + def __getstate__(self): state = self.__dict__.copy() del state["_log_likelihood_function"] @@ -153,6 +166,11 @@ def __call__(self, parameters, *kwargs): log_prior_list = self.model.log_prior_list_from_vector(vector=parameters) figure_of_merit = log_likelihood + sum(log_prior_list) + if self.store_history: + + self.parameters_history_list.append(parameters) + self.log_likelihood_history_list.append(log_likelihood) + if self.convert_to_chi_squared: figure_of_merit *= -2.0 diff --git a/autofit/non_linear/grid/grid_search/result.py b/autofit/non_linear/grid/grid_search/result.py index 3b1bef042..d85dcab27 100644 --- a/autofit/non_linear/grid/grid_search/result.py +++ b/autofit/non_linear/grid/grid_search/result.py @@ -1,4 +1,4 @@ -from typing import List, Optional, Union, Iterable +from typing import List, Optional, Union, Iterable, Tuple import numpy as np @@ -11,8 +11,66 @@ from autofit.non_linear.samples.interface import SamplesInterface +class AbstractGridSearchResult: + def __init__(self, samples: GridList): + self.samples = samples + + # noinspection PyTypeChecker + @as_grid_list + def physical_centres_lists_from( + self, + path: Union[str, Tuple[str, ...]], + ) -> GridList: + """ + Get the physical centres of the grid search from a path to an attribute of the instance in the samples. + + Parameters + ---------- + path + The path to the attribute to get from the instance + + Returns + ------- + A list of lists of physical values + """ + return self._physical_centres_lists_from(self.samples, path) + + # noinspection PyTypeChecker + @as_grid_list + def _physical_centres_lists_from( + self, + samples: GridList, + path: Union[str, Tuple[str, ...]], + ) -> GridList: + """ + Get the physical centres of the grid search from a path to an attribute of the instance in the samples. + + Parameters + ---------- + path + The path to the attribute to get from the instance + + Returns + ------- + A list of lists of physical values + """ + if isinstance(path, str): + path = path.split(".") + + def value_for_samples(samples): + return samples.model.object_for_path(path).mean + + else: + paths = [p.split(".") for p in path] + + def value_for_samples(samples): + return tuple(samples.model.object_for_path(p).mean for p in paths) + + return [value_for_samples(samples) for samples in samples] + + # noinspection PyTypeChecker -class GridSearchResult: +class GridSearchResult(AbstractGridSearchResult): def __init__( self, samples: List[SamplesInterface], @@ -34,13 +92,14 @@ def __init__( self.no_steps = len(lower_limits_lists) self.lower_limits_lists = GridList(lower_limits_lists, self.shape) - self.samples = GridList(samples, self.shape) if samples is not None else None self.side_length = int(self.no_steps ** (1 / self.no_dimensions)) self.step_size = 1 / self.side_length self.grid_priors = grid_priors self.parent = parent + super().__init__(GridList(samples, self.shape) if samples is not None else None) + @property @as_grid_list def physical_lower_limits_lists(self) -> GridList: @@ -209,7 +268,8 @@ def attribute_grid(self, attribute_path: Union[str, Iterable[str]]) -> GridList: @as_grid_list def log_likelihoods( - self, relative_to_value: float = 0.0, + self, + relative_to_value: float = 0.0, ) -> GridList: """ The maximum log likelihood of every grid search on a NumPy array whose shape is the native dimensions of the diff --git a/autofit/non_linear/grid/sensitivity/__init__.py b/autofit/non_linear/grid/sensitivity/__init__.py index 90c6107b1..ad83f2065 100644 --- a/autofit/non_linear/grid/sensitivity/__init__.py +++ b/autofit/non_linear/grid/sensitivity/__init__.py @@ -1,7 +1,7 @@ -import csv import logging import os from copy import copy +import numpy as np from pathlib import Path from typing import List, Generator, Callable, ClassVar, Optional, Union, Tuple @@ -9,11 +9,11 @@ from autofit.mapper.model import ModelInstance from autofit.mapper.prior_model.abstract import AbstractPriorModel from autofit.non_linear.grid.grid_search import make_lists, Sequential -from autofit.non_linear.grid.sensitivity.job import Job +from autofit.non_linear.grid.sensitivity.job import Job, MaskedJobResult from autofit.non_linear.grid.sensitivity.job import JobResult from autofit.non_linear.grid.sensitivity.result import SensitivityResult from autofit.non_linear.parallel import Process -from autofit.text.text_util import padding +from autofit.text.formatter import write_table class Sensitivity: @@ -29,6 +29,7 @@ def __init__( job_cls: ClassVar = Job, perturb_model_prior_func: Optional[Callable] = None, number_of_steps: Union[Tuple[int, ...], int] = 4, + mask: Optional[List[bool]] = None, number_of_cores: int = 2, limit_scale: int = 1, ): @@ -64,6 +65,10 @@ def __init__( The number of steps for each dimension of the sensitivity grid. If input as a float the dimensions are all that value. If input as a tuple of length the number of dimensions, each tuple value is the number of steps in that dimension. + mask + A mask to apply to the sensitivity grid, such that all `True` values are not included in the sensitivity + mapping. This is useful for removing regions of the sensitivity grid that are expected to have no + sensitivity, for example because they have no signal. number_of_cores How many cores does this computer have? limit_scale @@ -93,6 +98,22 @@ def __init__( self.job_cls = job_cls self.number_of_steps = number_of_steps + self.mask = None + + if mask is not None: + self.mask = np.asarray(mask) + if self.shape != self.mask.shape: + raise ValueError( + f""" + The mask of the Sensitivity object must have the same shape as the sensitivity grid. + + For your inputs, the shape of each are as follows: + + Sensitivity Grid: {self.shape} + Mask: {self.mask.shape} + """ + ) + self.number_of_cores = number_of_cores self.limit_scale = limit_scale @@ -116,9 +137,24 @@ def run(self) -> SensitivityResult: process_class = Process if self.number_of_cores > 1 else Sequential - results = list() + results = [] + jobs = [] + + for number in range(len(self._perturb_instances)): + if self._should_bypass(number=number): + model = self.model.copy() + model.perturb = self._perturb_models[number] + results.append( + MaskedJobResult( + number=number, + model=model, + ) + ) + else: + jobs.append(self._make_job(number)) + for result in process_class.run_jobs( - self._make_jobs(), number_of_cores=self.number_of_cores + jobs, number_of_cores=self.number_of_cores ): if isinstance(result, Exception): raise result @@ -127,32 +163,33 @@ def run(self) -> SensitivityResult: results = sorted(results) os.makedirs(self.paths.output_path, exist_ok=True) - with open(self.results_path, "w+") as f: - writer = csv.writer(f) - writer.writerow(headers) - for result_ in results: - values = physical_values[result_.number] - writer.writerow( - padding(item) - for item in [ - result_.number, - *values, - result_.log_evidence_increase, - result_.log_likelihood_increase, - ] - ) - result = SensitivityResult( + write_table( + headers=headers, + rows=[ + [ + result.number, + *physical_values[result.number], + result.log_evidence_increase, + result.log_likelihood_increase, + ] + for result in results + ], + filename=self.results_path, + ) + + sensitivity_result = SensitivityResult( samples=[result.result.samples_summary for result in results], perturb_samples=[ result.perturb_result.samples_summary for result in results ], shape=self.shape, + path_values=self.path_values, ) - self.paths.save_json("result", to_dict(result)) + self.paths.save_json("result", to_dict(sensitivity_result)) - return result + return sensitivity_result @property def shape(self) -> Tuple[int, ...]: @@ -176,6 +213,21 @@ def shape(self) -> Tuple[int, ...]: self.number_of_steps for _ in range(self.perturb_model.prior_count) ) + def shape_index_from_number(self, number: int) -> Tuple[int, ...]: + """ + Returns the index of the sensitivity grid from a number. + + Parameters + ---------- + number + The number of the sensitivity grid. + + Returns + ------- + The index of the sensitivity grid. + """ + return np.unravel_index(number, self.shape) + @property def step_size(self) -> Union[float, Tuple]: """ @@ -203,6 +255,17 @@ def _lists(self) -> List[List[float]]: """ return make_lists(self.perturb_model.prior_count, step_size=self.step_size) + @property + def path_values(self): + paths = [ + self.perturb_model.path_for_prior(prior) + for prior in self.perturb_model.priors_ordered_by_id + ] + + return { + path: list(values) for path, *values in zip(paths, *self._physical_values) + } + @property def _physical_values(self) -> List[List[float]]: """ @@ -227,31 +290,36 @@ def _headers(self) -> Generator[str, None, None]: yield path @property - def _labels(self) -> Generator[str, None, None]: + def _labels(self) -> List[str]: """ One label for each perturbation, used to distinguish fits for each perturbation by placing them in separate directories. """ + labels = [] for list_ in self._lists: strings = list() for value, prior_tuple in zip(list_, self.perturb_model.prior_tuples): path, prior = prior_tuple value = prior.value_for(value) strings.append(f"{path}_{value}") - yield "_".join(strings) + labels.append("_".join(strings)) + + return labels @property - def _perturb_instances(self) -> Generator[ModelInstance, None, None]: + def _perturb_instances(self) -> List[ModelInstance]: """ A list of instances each of which defines a perturbation to be applied to the image. """ - for list_ in self._lists: - yield self.perturb_model.instance_from_unit_vector(list_) + + return [ + self.perturb_model.instance_from_unit_vector(list_) for list_ in self._lists + ] @property - def _perturb_models(self) -> Generator[AbstractPriorModel, None, None]: + def _perturb_models(self) -> List[AbstractPriorModel]: """ A list of models representing a perturbation at each grid square. @@ -265,6 +333,8 @@ def _perturb_models(self) -> Generator[AbstractPriorModel, None, None]: step_sizes = (self.step_size,) * self.perturb_model.prior_count half_steps = [self.limit_scale * step_size / 2 for step_size in step_sizes] + + perturb_models = [] for list_ in self._lists: limits = [ ( @@ -277,38 +347,48 @@ def _perturb_models(self) -> Generator[AbstractPriorModel, None, None]: half_steps, ) ] - yield self.perturb_model.with_limits(limits) + perturb_models.append(self.perturb_model.with_limits(limits)) + return perturb_models + + def _should_bypass(self, number: int) -> bool: + shape_index = self.shape_index_from_number(number=number) + return self.mask is not None and np.asarray(self.mask)[shape_index] def _make_jobs(self) -> Generator[Job, None, None]: + for number, _ in enumerate(self._perturb_instances): + yield self._make_job(number) + + def _make_job(self, number) -> Generator[Job, None, None]: """ Create a list of jobs to be run on separate processes. Each job fits a perturb image with the original model and a model which includes a perturbation. """ - for number, (perturb_instance, perturb_model, label) in enumerate( - zip(self._perturb_instances, self._perturb_models, self._labels) - ): - if self.perturb_model_prior_func is not None: - perturb_model = self.perturb_model_prior_func( - perturb_instance=perturb_instance, perturb_model=perturb_model - ) + perturb_instance = self._perturb_instances[number] + perturb_model = self._perturb_models[number] + label = self._labels[number] - simulate_instance = copy(self.instance) - simulate_instance.perturb = perturb_instance - - paths = self.paths.for_sub_analysis( - label, + if self.perturb_model_prior_func is not None: + perturb_model = self.perturb_model_prior_func( + perturb_instance=perturb_instance, perturb_model=perturb_model ) - yield self.job_cls( - simulate_instance=simulate_instance, - model=self.model, - perturb_model=perturb_model, - base_instance=self.instance, - simulate_cls=self.simulate_cls, - base_fit_cls=self.base_fit_cls, - perturb_fit_cls=self.perturb_fit_cls, - paths=paths, - number=number, - ) + simulate_instance = copy(self.instance) + simulate_instance.perturb = perturb_instance + + paths = self.paths.for_sub_analysis( + label, + ) + + return self.job_cls( + simulate_instance=simulate_instance, + model=self.model, + perturb_model=perturb_model, + base_instance=self.instance, + simulate_cls=self.simulate_cls, + base_fit_cls=self.base_fit_cls, + perturb_fit_cls=self.perturb_fit_cls, + paths=paths, + number=number, + ) diff --git a/autofit/non_linear/grid/sensitivity/job.py b/autofit/non_linear/grid/sensitivity/job.py index c8f1b2ca3..ef14579c8 100644 --- a/autofit/non_linear/grid/sensitivity/job.py +++ b/autofit/non_linear/grid/sensitivity/job.py @@ -35,7 +35,10 @@ def log_evidence_increase(self) -> Optional[float]: if hasattr(self.result.samples, "log_evidence"): if self.result.samples.log_evidence is not None: - return float(self.perturb_result.samples.log_evidence - self.result.samples.log_evidence) + return float( + self.perturb_result.samples.log_evidence + - self.result.samples.log_evidence + ) @property def log_likelihood_increase(self) -> Optional[float]: @@ -48,6 +51,39 @@ def log_likelihood_increase(self) -> Optional[float]: return float(self.perturb_result.log_likelihood - self.result.log_likelihood) +class MaskedJobResult(AbstractJobResult): + """ + A placeholder result for a job that has been masked out. + """ + + def __init__(self, number, model): + super().__init__(number) + self.model = model + + @property + def result(self): + return self + + @property + def perturb_result(self): + return self + + def __getattr__(self, item): + return None + + @property + def samples_summary(self): + return self + + @property + def log_evidence(self): + return 0.0 + + @property + def log_likelihood(self): + return 0.0 + + class Job(AbstractJob): _number = count() @@ -92,6 +128,24 @@ def __init__( self.perturb_fit_cls = perturb_fit_cls self.paths = paths + @property + def base_paths(self): + return self.paths.for_sub_analysis("[base]") + + @property + def perturb_paths(self): + return self.paths.for_sub_analysis("[perturb]") + + @property + def is_complete(self) -> bool: + """ + Returns True if the job has been completed, False otherwise. + """ + return (self.base_paths.is_complete and self.perturb_paths.is_complete) or ( + (self.paths.output_path / "[base].zip").exists() + and (self.paths.output_path / "[perturb].zip").exists() + ) + def perform(self) -> JobResult: """ - Create one model with a perturbation and another without @@ -102,26 +156,33 @@ def perform(self) -> JobResult: An object comprising the results of the two fits """ - dataset = self.simulate_cls( - instance=self.simulate_instance, - simulate_path=self.paths.image_path.with_name("simulate"), - ) + if self.is_complete: + dataset = None + else: + dataset = self.simulate_cls( + instance=self.simulate_instance, + simulate_path=self.paths.image_path.with_name("simulate"), + ) result = self.base_fit_cls( model=self.model, dataset=dataset, - paths=self.paths.for_sub_analysis("[base]"), + paths=self.base_paths, + instance=self.simulate_instance, ) perturb_model = copy(self.model) - perturb_model.perturbation = self.perturb_model + perturb_model.perturb = self.perturb_model perturb_result = self.perturb_fit_cls( model=perturb_model, dataset=dataset, - paths=self.paths.for_sub_analysis("[perturb]"), + paths=self.perturb_paths, + instance=self.simulate_instance, ) return JobResult( - number=self.number, result=result, perturb_result=perturb_result + number=self.number, + result=result, + perturb_result=perturb_result, ) diff --git a/autofit/non_linear/grid/sensitivity/result.py b/autofit/non_linear/grid/sensitivity/result.py index 2b1a78093..8e5fd4839 100644 --- a/autofit/non_linear/grid/sensitivity/result.py +++ b/autofit/non_linear/grid/sensitivity/result.py @@ -1,29 +1,48 @@ -from typing import List, Tuple +from typing import List, Tuple, Union, Dict from autofit.non_linear.grid.grid_list import GridList, as_grid_list +from autofit.non_linear.grid.grid_search.result import AbstractGridSearchResult from autofit.non_linear.samples.interface import SamplesInterface + # noinspection PyTypeChecker -class SensitivityResult: +class SensitivityResult(AbstractGridSearchResult): def __init__( - self, - samples: List[SamplesInterface], - perturb_samples: List[SamplesInterface], - shape : Tuple[int, ...] + self, + samples: List[SamplesInterface], + perturb_samples: List[SamplesInterface], + shape: Tuple[int, ...], + path_values: Dict[Tuple[str, ...], List[float]], ): """ The result of a sensitivity mapping Parameters ---------- - results - The results of each sensitivity job. shape The shape of the sensitivity mapping grid. + path_values + A list of tuples of the path to the grid priors and the physical values themselves. """ - self.samples = GridList(samples, shape) + super().__init__(GridList(samples, shape)) self.perturb_samples = GridList(perturb_samples, shape) self.shape = shape + self.path_values = path_values + + def perturbed_physical_centres_list_from( + self, path: Union[str, Tuple[str, ...]] + ) -> GridList: + """ + Returns the physical centres of the perturbed model for each sensitivity fit + + Parameters + ---------- + path + The path to the physical centres in the samples + """ + if isinstance(path, str): + path = tuple(path.split(".")) + return self.path_values[path] def __getitem__(self, item): return self.samples[item] @@ -57,9 +76,10 @@ def log_evidence_differences(self) -> GridList: The log evidence differences between the base and perturbed models """ return [ - log_evidence_perturbed - log_evidence_base for - log_evidence_perturbed, log_evidence_base in - zip(self.log_evidences_perturbed, self.log_evidences_base) + log_evidence_perturbed - log_evidence_base + for log_evidence_perturbed, log_evidence_base in zip( + self.log_evidences_perturbed, self.log_evidences_base + ) ] @property @@ -85,13 +105,15 @@ def log_likelihood_differences(self) -> GridList: The log likelihood differences between the base and perturbed models """ return [ - log_likelihood_perturbed - log_likelihood_base for - log_likelihood_perturbed, log_likelihood_base in - zip(self.log_likelihoods_perturbed, self.log_likelihoods_base) + log_likelihood_perturbed - log_likelihood_base + for log_likelihood_perturbed, log_likelihood_base in zip( + self.log_likelihoods_perturbed, self.log_likelihoods_base + ) ] def figure_of_merits( - self, use_log_evidences: bool, + self, + use_log_evidences: bool, ) -> GridList: """ Convenience method to get either the log likelihoods difference or log evidence difference of the grid search. @@ -105,6 +127,3 @@ def figure_of_merits( if use_log_evidences: return self.log_evidence_differences return self.log_likelihood_differences - - - diff --git a/autofit/non_linear/initializer.py b/autofit/non_linear/initializer.py index 0de568cf6..a0966e531 100644 --- a/autofit/non_linear/initializer.py +++ b/autofit/non_linear/initializer.py @@ -25,6 +25,9 @@ class AbstractInitializer(ABC): def _generate_unit_parameter_list(self, model): pass + def info_from_model(self, model : AbstractPriorModel) -> str: + raise NotImplementedError + @staticmethod def figure_of_metric(args) -> Optional[float]: fitness, parameter_list = args @@ -175,7 +178,7 @@ def samples_in_test_mode(self, total_points: int, model: AbstractPriorModel): return unit_parameter_lists, parameter_lists, figure_of_merit_list -class SpecificRangeInitializer(AbstractInitializer): +class InitializerParamBounds(AbstractInitializer): def __init__( self, parameter_dict: Dict[Prior, Tuple[float, float]], @@ -183,14 +186,14 @@ def __init__( upper_limit=1.0, ): """ - Initializer that allows the range of possible starting points for each prior - to be specified explicitly. + Initializer which uses the bounds on input parameters as the starting point for the search (e.g. where + an MLE optimization starts or MCMC walkers are initialized). Parameters ---------- parameter_dict - A dictionary mapping priors to inclusive ranges of physical values that - the initial values for that dimension in the search may take + A dictionary mapping each parameter path to bounded ranges of physical values that + are where the search begins. lower_limit A default, unit lower limit used when a prior is not specified upper_limit @@ -226,7 +229,7 @@ def _generate_unit_parameter_list(self, model: AbstractPriorModel) -> List[float key = ".".join(model.path_for_prior(prior)) if key not in self._generated_warnings: logger.warning( - f"Range for {key} not set in the SpecificRangeInitializer. " + f"Range for {key} not set in the InitializerParamBounds. " f"Using defaults." ) self._generated_warnings.add(key) @@ -240,6 +243,84 @@ def _generate_unit_parameter_list(self, model: AbstractPriorModel) -> List[float return unit_parameter_list + def info_from_model(self, model : AbstractPriorModel) -> str: + """ + Returns a string showing the bounds of the parameters in the initializer. + """ + info = "Total Free Parameters = " + str(model.prior_count) + "\n" + info += "Total Starting Points = " + str(len(self.parameter_dict)) + "\n\n" + for prior in model.priors_ordered_by_id: + + key = ".".join(model.path_for_prior(prior)) + + try: + + value = self.info_value_from(self.parameter_dict[prior]) + + info += f"{key}: Start[{value}]\n" + + except KeyError: + + info += f"{key}: {prior})\n" + + return info + + def info_value_from(self, value : Tuple[float, float]) -> Tuple[float, float]: + """ + Returns the value that is used to display the bounds of the parameters in the initializer. + + This function simply returns the input value, but it can be overridden in subclasses for diffferent + initializers. + + Parameters + ---------- + value + The value to be displayed in the initializer info which is a tuple of the lower and upper bounds of the + parameter. + """ + return value + + +class InitializerParamStartPoints(InitializerParamBounds): + def __init__( + self, + parameter_dict: Dict[Prior, float], + ): + """ + Initializer which input values of the parameters as the starting point for the search (e.g. where + an MLE optimization starts or MCMC walkers are initialized). + + Parameters + ---------- + parameter_dict + A dictionary mapping each parameter path to the starting point physical values that + are where the search begins. + lower_limit + A default, unit lower limit used when a prior is not specified + upper_limit + A default, unit upper limit used when a prior is not specified + """ + parameter_dict_new = {} + + for key, value in parameter_dict.items(): + parameter_dict_new[key] = (value - 1.0e-8, value + 1.0e-8) + + super().__init__(parameter_dict=parameter_dict_new) + + def info_value_from(self, value : Tuple[float, float]) -> float: + """ + Returns the value that is used to display the starting point of the parameters in the initializer. + + This function returns the mean of the input value, as the starting point is a single value in the center of the + bounds. + + Parameters + ---------- + value + The value to be displayed in the initializer info which is a tuple of the lower and upper bounds of the + parameter. + """ + return (value[1] + value[0]) / 2.0 class Initializer(AbstractInitializer): def __init__(self, lower_limit: float, upper_limit: float): diff --git a/autofit/non_linear/mock/mock_search.py b/autofit/non_linear/mock/mock_search.py index 8f26e4c88..c68b5d9da 100644 --- a/autofit/non_linear/mock/mock_search.py +++ b/autofit/non_linear/mock/mock_search.py @@ -136,6 +136,7 @@ def _fit(self, model, analysis): ) self.paths.save_samples_summary(self.samples_summary) + self.paths.completed() return analysis.make_result( samples_summary=samples_summary, @@ -157,13 +158,13 @@ def perform_update(self, model, analysis, during_analysis, search_internal=None) ) -class MockOptimizer(MockSearch): +class MockMLE(MockSearch): def __init__(self, **kwargs): super().__init__(fit_fast=False, **kwargs) @property def samples_cls(self): - return MockOptimizer + return MockMLE def project( self, factor_approx: FactorApproximation, status: Status = Status() diff --git a/autofit/non_linear/paths/abstract.py b/autofit/non_linear/paths/abstract.py index 4ae9f4c37..d686a5d43 100644 --- a/autofit/non_linear/paths/abstract.py +++ b/autofit/non_linear/paths/abstract.py @@ -32,6 +32,7 @@ def __init__( parent: Optional["AbstractPaths"] = None, unique_tag: Optional[str] = None, identifier: str = None, + image_path_suffix: str = "", ): """ Manages the path structure for `NonLinearSearch` output, for analyses both not using and using the search @@ -63,6 +64,16 @@ def __init__( is_identifier_in_paths If True output path and symlink path terminate with an identifier generated from the search and model + parent + The parent paths object of this paths object. + unique_tag + A unique tag for the search, used to differentiate between searches with the same name. + identifier + A custom identifier for the search, if this is not None it will be used instead of the automatically + generated identifier + image_path_suffix + A suffix which is appended to the image path. This is used to differentiate between different + image outputs, for example the image of the starting point of an MLE. """ self.name = name or "" @@ -87,6 +98,8 @@ def __init__( except NoSectionError as e: logger.exception(e) + self.image_path_suffix = image_path_suffix + @property @abstractmethod def samples(self): @@ -211,10 +224,10 @@ def image_path(self) -> Path: The path to the image folder. """ - if not os.path.exists(self.output_path / "image"): - os.makedirs(self.output_path / "image") + if not os.path.exists(self.output_path / f"image{self.image_path_suffix}"): + os.makedirs(self.output_path / f"image{self.image_path_suffix}") - return self.output_path / "image" + return self.output_path / f"image{self.image_path_suffix}" @property def profile_path(self) -> Path: @@ -231,10 +244,11 @@ def output_path(self) -> Path: strings = list( filter( - len, + None, [ str(conf.instance.output_path), str(self.path_prefix), + self.unique_tag, str(self.name), ], ) diff --git a/autofit/non_linear/paths/directory.py b/autofit/non_linear/paths/directory.py index 69452603e..6f7d429b2 100644 --- a/autofit/non_linear/paths/directory.py +++ b/autofit/non_linear/paths/directory.py @@ -2,6 +2,7 @@ import dill import json +import numpy as np import os from os import path from pathlib import Path @@ -21,7 +22,6 @@ from ..samples import load_from_table from autofit.non_linear.samples.pdf import SamplesPDF from autofit.non_linear.samples.summary import SamplesSummary -import numpy as np from ...visualise import VisualiseGraph @@ -348,9 +348,17 @@ def save_all(self, search_config_dict=None, info=None): VisualiseGraph( model=self.model, ).save(str(self.output_path / "model_graph.html")) + if info: self.save_json("info", info) + self.save_json("search", to_dict(self.search)) + try: + info_start = self.search.initializer.info_from_model(model=self.model) + self._save_model_start_point(info=info_start) + except (NotImplementedError, AttributeError): + pass + self.save_json("model", to_dict(self.model)) self._save_metadata(search_name=type(self.search).__name__.lower()) @@ -458,6 +466,13 @@ def _save_model_info(self, model): with open_(self.output_path / "model.info", "w+") as f: f.write(model.info) + def _save_model_start_point(self, info): + """ + Save the model.start file, which summarizes the start point of every parameter. + """ + with open_(self.output_path / "model.start", "w+") as f: + f.write(info) + def _save_parameter_names_file(self, model): """ Create the param_names file listing every parameter's label and Latex tag, which is used for corner.py diff --git a/autofit/non_linear/plot/mle_plotters.py b/autofit/non_linear/plot/mle_plotters.py new file mode 100644 index 000000000..4e31e8a5d --- /dev/null +++ b/autofit/non_linear/plot/mle_plotters.py @@ -0,0 +1,134 @@ +import matplotlib.pyplot as plt +import logging +import numpy as np + +from autofit.non_linear.plot.samples_plotters import SamplesPlotter + +from autofit.non_linear.plot.samples_plotters import skip_plot_in_test_mode + +logger = logging.getLogger(__name__) + +class MLEPlotter(SamplesPlotter): + + def subplot_parameters(self, use_log_y : bool = False, use_last_50_percent : bool = False, **kwargs): + """ + Plots a subplot of every parameter against iteration number. + + The subplot extends over all free parameters in the model-fit, with the number of parameters per subplot + given by the total number of free parameters in the model-fit. + + This often produces a large dynamic range in the y-axis. Plotting the y-axis on a log-scale or only + plotting the last 50% of samples can make the plot easier to inspect. + + Parameters + ---------- + use_log_y + If True, the y-axis is plotted on a log-scale. + use_last_50_percent + If True, only the last 50% of samples are plotted. + kwargs + Additional key word arguments can be passed to the `plt.subplots` method. + + Returns + ------- + + """ + + parameter_lists = self.samples.parameters_extract + + plt.subplots(self.model.total_free_parameters, 1, figsize=(12, 3 * len(parameter_lists))) + + for i, parameters in enumerate(parameter_lists): + + iteration_list = range(len(parameter_lists[0])) + + plt.subplot(self.model.total_free_parameters, 1, i + 1) + + if use_last_50_percent: + + iteration_list = iteration_list[int(len(iteration_list) / 2) :] + parameters = parameters[int(len(parameters) / 2) :] + + if use_log_y: + plt.semilogy(iteration_list, parameters, c="k") + else: + plt.plot(iteration_list, parameters, c="k") + + plt.xlabel("Iteration", fontsize=16) + plt.ylabel(self.model.parameter_labels_with_superscripts_latex[i], fontsize=16) + plt.xticks(fontsize=16) + plt.yticks(fontsize=16) + + filename = "subplot_parameters" + + if use_log_y: + filename += "_log_y" + + if use_last_50_percent: + filename += "_last_50_percent" + + self.output.subplot_to_figure( + auto_filename=filename + ) + plt.close() + + @skip_plot_in_test_mode + def log_likelihood_vs_iteration(self, use_log_y : bool = False, use_last_50_percent : bool = False, **kwargs): + """ + Plot the log likelihood of a model fit as a function of iteration number. + + For a maximum likelihood estimate, the log likelihood should increase with iteration number. + + This often produces a large dynamic range in the y-axis. Plotting the y-axis on a log-scale or only + plotting the last 50% of samples can make the plot easier to inspect. + + Parameters + ---------- + use_log_y + If True, the y-axis is plotted on a log-scale. + """ + + log_likelihood_list = self.samples.log_likelihood_list + iteration_list = range(len(log_likelihood_list)) + + if use_last_50_percent: + + iteration_list = iteration_list[int(len(iteration_list) / 2) :] + log_likelihood_list = log_likelihood_list[int(len(log_likelihood_list) / 2) :] + + plt.figure(figsize=(12, 12)) + + if use_log_y: + plt.semilogy(iteration_list, log_likelihood_list, c="k") + else: + plt.plot(iteration_list, log_likelihood_list, c="k") + + plt.xlabel("Iteration", fontsize=16) + plt.ylabel("Log Likelihood", fontsize=16) + plt.xticks(fontsize=16) + plt.yticks(fontsize=16) + + title = "Log Likelihood vs Iteration" + + if use_log_y: + + title += " (Log Scale)" + + if use_last_50_percent: + + title += " (Last 50 Percent)" + + plt.title("Log Likelihood vs Iteration", fontsize=24) + + filename = "log_likelihood_vs_iteration" + + if use_log_y: + filename += "_log_y" + + if use_last_50_percent: + filename += "_last_50_percent" + + self.output.to_figure( + auto_filename=filename, + ) + plt.close() diff --git a/autofit/non_linear/plot/nest_plotters.py b/autofit/non_linear/plot/nest_plotters.py index 64d4cd43a..2d4c470cb 100644 --- a/autofit/non_linear/plot/nest_plotters.py +++ b/autofit/non_linear/plot/nest_plotters.py @@ -1,5 +1,3 @@ -from anesthetic.samples import NestedSamples -from anesthetic import make_2d_axes from functools import wraps import numpy as np import warnings @@ -59,6 +57,8 @@ def corner_anesthetic(self, **kwargs): config_dict = conf.instance["visualize"]["plots_settings"]["corner_anesthetic"] + from anesthetic.samples import NestedSamples + from anesthetic import make_2d_axes import matplotlib.pylab as pylab params = {'font.size' : int(config_dict["fontsize"])} diff --git a/autofit/non_linear/plot/optimize_plotters.py b/autofit/non_linear/plot/optimize_plotters.py deleted file mode 100644 index a4a6b8a31..000000000 --- a/autofit/non_linear/plot/optimize_plotters.py +++ /dev/null @@ -1,12 +0,0 @@ -import matplotlib.pyplot as plt -import logging -import numpy as np - -from autofit.non_linear.plot.samples_plotters import SamplesPlotter - -logger = logging.getLogger(__name__) - -class OptimizePlotter(SamplesPlotter): - - pass - diff --git a/autofit/non_linear/plot/samples_plotters.py b/autofit/non_linear/plot/samples_plotters.py index 37949ec14..885f01ef1 100644 --- a/autofit/non_linear/plot/samples_plotters.py +++ b/autofit/non_linear/plot/samples_plotters.py @@ -71,6 +71,7 @@ def log_posterior_list(self): def close(self): if plt.fignum_exists(num=1): + plt.clf() plt.close() def log_plot_exception(self, plot_name : str): diff --git a/autofit/non_linear/samples/sample.py b/autofit/non_linear/samples/sample.py index 45bf67092..60051a826 100644 --- a/autofit/non_linear/samples/sample.py +++ b/autofit/non_linear/samples/sample.py @@ -175,7 +175,11 @@ def from_lists( ) return samples - def instance_for_model(self, model: AbstractPriorModel): + def instance_for_model( + self, + model: AbstractPriorModel, + ignore_assertions: bool = False, + ): """ Create an instance from this sample for a model @@ -183,6 +187,8 @@ def instance_for_model(self, model: AbstractPriorModel): ---------- model The model the this sample was taken from + ignore_assertions + If True, do not check that the instance is valid Returns ------- @@ -190,13 +196,22 @@ def instance_for_model(self, model: AbstractPriorModel): """ try: if self.is_path_kwargs: - return model.instance_from_path_arguments(self.kwargs) + return model.instance_from_path_arguments( + self.kwargs, + ignore_assertions=ignore_assertions, + ) else: - return model.instance_from_prior_name_arguments(self.kwargs) + return model.instance_from_prior_name_arguments( + self.kwargs, + ignore_assertions=ignore_assertions, + ) except KeyError: # TODO: Does this get used? If so, why? - return model.instance_from_vector(self.parameter_lists_for_model(model)) + return model.instance_from_vector( + self.parameter_lists_for_model(model), + ignore_prior_limits=ignore_assertions, + ) @split_paths def with_paths(self, paths: List[Tuple[str, ...]]) -> "Sample": diff --git a/autofit/non_linear/samples/samples.py b/autofit/non_linear/samples/samples.py index 00165df5f..742270b38 100644 --- a/autofit/non_linear/samples/samples.py +++ b/autofit/non_linear/samples/samples.py @@ -407,7 +407,6 @@ def minimise(self) -> "Samples": A copy of this object with only important samples retained """ samples = copy(self) - samples.model = None samples.sample_list = list( {self.max_log_likelihood_sample, self.max_log_posterior_sample} ) diff --git a/autofit/non_linear/samples/summary.py b/autofit/non_linear/samples/summary.py index 3c9f59c74..920868549 100644 --- a/autofit/non_linear/samples/summary.py +++ b/autofit/non_linear/samples/summary.py @@ -12,8 +12,6 @@ class SamplesSummary(SamplesInterface): - __exclude_fields__ = ["model"] - def __init__( self, max_log_likelihood_sample: Sample, diff --git a/autofit/non_linear/search/abstract_search.py b/autofit/non_linear/search/abstract_search.py index 2ac13edec..3fdb70d25 100644 --- a/autofit/non_linear/search/abstract_search.py +++ b/autofit/non_linear/search/abstract_search.py @@ -29,6 +29,7 @@ from autofit.graphical.utils import Status from autofit.mapper.prior_model.abstract import AbstractPriorModel from autofit.mapper.prior_model.collection import Collection +from autofit.mapper.model import ModelInstance from autofit.non_linear.initializer import Initializer from autofit.non_linear.fitness import Fitness from autofit.non_linear.parallel import SneakyPool, SneakierPool @@ -181,9 +182,6 @@ def __init__( self._logger = None - if unique_tag is not None and path_prefix is not None: - path_prefix = path_prefix / unique_tag - self.unique_tag = unique_tag if paths: @@ -652,10 +650,11 @@ def pre_fit_output( if not self.disable_output: self.logger.info(f"The output path of this fit is {self.paths.output_path}") else: - self.logger.info("Output to hard-disk disabled, input a search name to enable.") + self.logger.info( + "Output to hard-disk disabled, input a search name to enable." + ) if not self.paths.is_complete or self.force_pickle_overwrite: - if not self.disable_output: self.logger.info( f"Outputting pre-fit files (e.g. model.info, visualization)." @@ -958,7 +957,6 @@ def perform_update( self.iterations += self.iterations_per_update if not self.disable_output: - if during_analysis: self.logger.info( f"""Fit Running: Updating results after {self.iterations} iterations (see output folder).""" @@ -1050,9 +1048,11 @@ def perform_update( def perform_visualization( self, model: AbstractPriorModel, - analysis: AbstractPriorModel, - samples_summary: SamplesSummary, + analysis: Analysis, during_analysis: bool, + samples_summary: Optional[SamplesSummary] = None, + instance: Optional[ModelInstance] = None, + paths_override: Optional[AbstractPaths] = None, search_internal=None, ): """ @@ -1075,32 +1075,54 @@ def perform_visualization( analysis Contains the data and the log likelihood function which fits an instance of the model to the data, returning the log likelihood the `NonLinearSearch` maximizes. + samples_summary + The summary of the samples of the non-linear search, which are used for visualization. during_analysis If the update is during a non-linear search, in which case tasks are only performed after a certain number of updates and only a subset of visualization may be performed. + instance + The instance of the model that is used for visualization. If not input, the maximum log likelihood + instance from the samples is used. """ self.logger.debug("Visualizing") - if analysis.should_visualize(paths=self.paths, during_analysis=during_analysis): + paths = paths_override or self.paths + + if instance is None and samples_summary is None: + raise AssertionError( + """The search's perform_visualization method has been called without an input instance or + samples_summary. + + This should not occur, please ensure one of these inputs is provided. + """ + ) + + if instance is None: + instance = samples_summary.instance + + if analysis.should_visualize(paths=paths, during_analysis=during_analysis): analysis.visualize( - paths=self.paths, - instance=samples_summary.instance, + paths=paths, + instance=instance, during_analysis=during_analysis, ) analysis.visualize_combined( - paths=self.paths, - instance=samples_summary.instance, + paths=paths, + instance=instance, during_analysis=during_analysis, ) - if analysis.should_visualize(paths=self.paths, during_analysis=during_analysis): - if not isinstance(self.paths, NullPaths): - samples = self.samples_from( - model=model, search_internal=search_internal - ) + if analysis.should_visualize(paths=paths, during_analysis=during_analysis): + if not isinstance(paths, NullPaths): + try: + samples = self.samples_from( + model=model, search_internal=search_internal + ) - self.plot_results(samples=samples) + self.plot_results(samples=samples) + except FileNotFoundError: + pass @property def samples_cls(self): diff --git a/autofit/non_linear/search/mcmc/emcee/search.py b/autofit/non_linear/search/mcmc/emcee/search.py index e8c2d81a1..59a511e81 100644 --- a/autofit/non_linear/search/mcmc/emcee/search.py +++ b/autofit/non_linear/search/mcmc/emcee/search.py @@ -4,6 +4,8 @@ import emcee import numpy as np +from autoconf import conf + from autofit.database.sqlalchemy_ import sa from autofit.mapper.model_mapper import ModelMapper from autofit.mapper.prior_model.abstract import AbstractPriorModel @@ -82,6 +84,11 @@ def __init__( self.logger.debug("Creating Emcee Search") + # TODO : Emcee visualization tools rely on the .hdf file and thus require that the search internal is + # TODO : On hard-disk, which this forces to occur. + + conf.instance["output"]["search_internal"] = True + def _fit(self, model: AbstractPriorModel, analysis): """ Fit a model using Emcee and the Analysis class which contains the data and returns the log likelihood from diff --git a/autofit/non_linear/search/mcmc/zeus/search.py b/autofit/non_linear/search/mcmc/zeus/search.py index ea1d3b11f..32e63242e 100644 --- a/autofit/non_linear/search/mcmc/zeus/search.py +++ b/autofit/non_linear/search/mcmc/zeus/search.py @@ -1,6 +1,8 @@ +import logging from typing import Dict, Optional import numpy as np +import os from autofit.database.sqlalchemy_ import sa from autofit.mapper.model_mapper import ModelMapper @@ -116,7 +118,7 @@ def _fit(self, model: AbstractPriorModel, analysis): "You are attempting to perform a model-fit using Zeus. \n\n" "However, the optional library Zeus (https://zeus-mcmc.readthedocs.io/en/latest/) is " "not installed.\n\n" - "Install it via the command `pip install zeus==3.5.5`.\n\n" + "Install it via the command `pip install zeus-mcmc==2.5.4`.\n\n" "----------------------" ) @@ -274,13 +276,39 @@ def samples_via_internal_from(self, model, search_internal=None): search_internal = search_internal or self.paths.load_search_internal() - auto_correlations = self.auto_correlations_from(search_internal=search_internal) + if os.environ.get("PYAUTOFIT_TEST_MODE") == "1": - discard = int(3.0 * np.max(auto_correlations.times)) - thin = int(np.max(auto_correlations.times) / 2.0) - samples_after_burn_in = search_internal.get_chain( - discard=discard, thin=thin, flat=True - ) + samples_after_burn_in = search_internal.get_chain( + discard=5, thin=5, flat=True + ) + + else: + auto_correlations = self.auto_correlations_from( + search_internal=search_internal + ) + + discard = int(3.0 * np.max(auto_correlations.times)) + thin = int(np.max(auto_correlations.times) / 2.0) + samples_after_burn_in = search_internal.get_chain( + discard=discard, thin=thin, flat=True + ) + + if len(samples_after_burn_in) == 0: + + logging.info( + """ + After thinnng the Zeus samples in order to remove burn-in, no samples were left. + + To create a samples object containing samples, so that the code can continue and results + can be inspected, the full list of samples before removing burn-in has been used. This may + indicate that the sampler has not converged and therefore your results may not be reliable. + + To fix this, run Zeus with more steps to ensure convergence is achieved or change the auto + correlation settings to be less aggressive in thinning samples. + """ + ) + + samples_after_burn_in = search_internal.get_chain(flat=True) parameter_lists = samples_after_burn_in.tolist() log_posterior_list = search_internal.get_log_prob(flat=True).tolist() @@ -306,7 +334,9 @@ def samples_via_internal_from(self, model, search_internal=None): sample_list=sample_list, samples_info=self.samples_info_from(search_internal=search_internal), auto_correlation_settings=self.auto_correlation_settings, - auto_correlations=auto_correlations, + auto_correlations=self.auto_correlations_from( + search_internal=search_internal + ), ) def auto_correlations_from(self, search_internal=None): diff --git a/autofit/non_linear/search/optimize/__init__.py b/autofit/non_linear/search/mle/__init__.py similarity index 100% rename from autofit/non_linear/search/optimize/__init__.py rename to autofit/non_linear/search/mle/__init__.py diff --git a/autofit/non_linear/search/mle/abstract_mle.py b/autofit/non_linear/search/mle/abstract_mle.py new file mode 100644 index 000000000..530fb51e8 --- /dev/null +++ b/autofit/non_linear/search/mle/abstract_mle.py @@ -0,0 +1,88 @@ +from abc import ABC +import copy +from typing import List + +from autoconf import conf +from autofit.mapper.prior_model.abstract import AbstractPriorModel +from autofit.non_linear.analysis import Analysis +from autofit.non_linear.search.abstract_search import NonLinearSearch +from autofit.non_linear.samples import Samples +from autofit.non_linear.plot.mle_plotters import MLEPlotter +from autofit.non_linear.plot.output import Output + + +class AbstractMLE(NonLinearSearch, ABC): + @property + def config_type(self): + return conf.instance["non_linear"]["mle"] + + @property + def samples_cls(self): + return Samples + + @property + def plotter_cls(self): + return MLEPlotter + + def plot_start_point( + self, + parameter_vector : List[float], + model: AbstractPriorModel, + analysis: Analysis, + ): + """ + Visualize the starting point of the non-linear search, using an instance of the model at the starting point + of the maximum likelihood estimator. + + Plots are output to a folder named `image_start` in the output path, so that the starting point model + can be compared to the final model inferred by the non-linear search. + + Parameters + ---------- + model + The model used by the non-linear search + analysis + The analysis which contains the visualization methods which plot the starting point model. + + Returns + ------- + + """ + + self.logger.info( + f"Visualizing Starting Point Model in image_start folder." + ) + + instance = model.instance_from_vector(vector=parameter_vector) + paths = copy.copy(self.paths) + paths.image_path_suffix = "_start" + + self.perform_visualization( + model=model, + analysis=analysis, + instance=instance, + during_analysis=False, + paths_override=paths, + ) + + def plot_results(self, samples): + + def should_plot(name): + return conf.instance["visualize"]["plots_search"]["mle"][name] + + plotter = self.plotter_cls( + samples=samples, + output=Output(path=self.paths.image_path / "search", format="png"), + ) + + if should_plot("subplot_parameters"): + + plotter.subplot_parameters() + plotter.subplot_parameters(use_log_y=True) + plotter.subplot_parameters(use_last_50_percent=True) + + if should_plot("log_likelihood_vs_iteration"): + + plotter.log_likelihood_vs_iteration() + plotter.log_likelihood_vs_iteration(use_log_y=True) + plotter.log_likelihood_vs_iteration(use_last_50_percent=True) \ No newline at end of file diff --git a/autofit/non_linear/search/optimize/drawer/__init__.py b/autofit/non_linear/search/mle/bfgs/__init__.py similarity index 100% rename from autofit/non_linear/search/optimize/drawer/__init__.py rename to autofit/non_linear/search/mle/bfgs/__init__.py diff --git a/autofit/non_linear/search/mle/bfgs/search.py b/autofit/non_linear/search/mle/bfgs/search.py new file mode 100644 index 000000000..e620fec36 --- /dev/null +++ b/autofit/non_linear/search/mle/bfgs/search.py @@ -0,0 +1,337 @@ +from typing import Optional + +from autoconf import cached_property +from autofit.database.sqlalchemy_ import sa + +from autofit.mapper.prior_model.abstract import AbstractPriorModel +from autofit.non_linear.search.mle.abstract_mle import AbstractMLE +from autofit.non_linear.analysis import Analysis +from autofit.non_linear.fitness import Fitness +from autofit.non_linear.initializer import AbstractInitializer +from autofit.non_linear.samples.sample import Sample +from autofit.non_linear.samples.samples import Samples + +import copy +from scipy import optimize +import numpy as np + + +class AbstractBFGS(AbstractMLE): + + method = None + + def __init__( + self, + name: Optional[str] = None, + path_prefix: Optional[str] = None, + unique_tag: Optional[str] = None, + initializer: Optional[AbstractInitializer] = None, + iterations_per_update: int = None, + session: Optional[sa.orm.Session] = None, + visualize : bool = False, + **kwargs + ): + """ + Abstract wrapper for the BFGS and L-BFGS scipy non-linear searches. + + See the docstrings of the `BFGS` and `LBFGS` classes for a description of the arguments of this class. + """ + + super().__init__( + name=name, + path_prefix=path_prefix, + unique_tag=unique_tag, + initializer=initializer, + iterations_per_update=iterations_per_update, + session=session, + **kwargs + ) + + self.visualize = visualize + + self.logger.debug(f"Creating {self.method} Search") + + @cached_property + def config_dict_options(self): + config_dict = copy.deepcopy(self._class_config["options"]) + + for key, value in config_dict.items(): + try: + config_dict[key] = self.kwargs[key] + except KeyError: + pass + + return config_dict + + def _fit( + self, + model: AbstractPriorModel, + analysis: Analysis, + ): + """ + Fit a model using the scipy L-BFGS method and the Analysis class which contains the data and returns the log + likelihood from instances of the model, which the `NonLinearSearch` seeks to maximize. + + Parameters + ---------- + model + The model which generates instances for different points in parameter space. + analysis + Contains the data and the log likelihood function which fits an instance of the model to the data, + returning the log likelihood the `NonLinearSearch` maximizes. + + Returns + ------- + A result object comprising the Samples object that inclues the maximum log likelihood instance and full + chains used by the fit. + """ + fitness = Fitness( + model=model, + analysis=analysis, + paths=self.paths, + fom_is_log_likelihood=False, + resample_figure_of_merit=-np.inf, + convert_to_chi_squared=True, + store_history=self.visualize, + ) + + try: + search_internal_dict = self.paths.load_search_internal() + + x0 = search_internal_dict["x0"] + total_iterations = search_internal_dict["total_iterations"] + + self.logger.info( + "Resuming LBFGS non-linear search (previous samples found)." + ) + + except (FileNotFoundError, TypeError): + + ( + unit_parameter_lists, + parameter_lists, + log_posterior_list, + ) = self.initializer.samples_from_model( + total_points=1, + model=model, + fitness=fitness, + paths=self.paths, + n_cores=self.number_of_cores, + ) + + x0 = np.asarray(parameter_lists[0]) + + total_iterations = 0 + + self.logger.info( + f"Starting new {self.method} non-linear search (no previous samples found)." + ) + + if self.visualize: + + self.plot_start_point( + parameter_vector=x0, + model=model, + analysis=analysis, + ) + + maxiter = self.config_dict_options.get("maxiter", 1e8) + + while total_iterations < maxiter: + iterations_remaining = maxiter - total_iterations + + iterations = min(self.iterations_per_update, iterations_remaining) + + if iterations > 0: + config_dict_options = self.config_dict_options + config_dict_options["maxiter"] = iterations + + search_internal = optimize.minimize( + fun=fitness.__call__, + x0=x0, + method=self.method, + options=config_dict_options, + **self.config_dict_search + ) + + total_iterations += search_internal.nit + + search_internal.log_posterior_list = -0.5 * fitness( + parameters=search_internal.x + ) + + if self.visualize: + + search_internal.parameters_history_list = fitness.parameters_history_list + search_internal.log_likelihood_history_list = fitness.log_likelihood_history_list + + self.paths.save_search_internal( + obj=search_internal, + ) + + x0 = search_internal.x + + if search_internal.nit < iterations: + return search_internal + + self.perform_update( + model=model, + analysis=analysis, + during_analysis=True, + search_internal=search_internal, + ) + + self.logger.info(f"{self.method} sampling complete.") + + return search_internal + + def samples_via_internal_from( + self, model: AbstractPriorModel, search_internal=None + ): + """ + Returns a `Samples` object from the LBFGS internal results. + + The samples contain all information on the parameter space sampling (e.g. the parameters, + log likelihoods, etc.). + + The internal search results are converted from the native format used by the search to lists of values + (e.g. `parameter_lists`, `log_likelihood_list`). + + Parameters + ---------- + model + Maps input vectors of unit parameter values to physical values and model instances via priors. + """ + + if search_internal is None: + search_internal = self.paths.load_search_internal() + + x0 = search_internal.x + total_iterations = search_internal.nit + + + if self.visualize: + + parameter_lists = search_internal.parameters_history_list + log_prior_list = model.log_prior_list_from(parameter_lists=parameter_lists) + log_likelihood_list = search_internal.log_likelihood_history_list + + else: + + parameter_lists = [list(x0)] + log_prior_list = model.log_prior_list_from(parameter_lists=parameter_lists) + log_posterior_list = np.array([search_internal.log_posterior_list]) + log_likelihood_list = [ + lp - prior for lp, prior in zip(log_posterior_list, log_prior_list) + ] + + weight_list = len(log_likelihood_list) * [1.0] + + sample_list = Sample.from_lists( + model=model, + parameter_lists=parameter_lists, + log_likelihood_list=log_likelihood_list, + log_prior_list=log_prior_list, + weight_list=weight_list, + ) + + samples_info = { + "total_iterations": total_iterations, + "time": self.timer.time if self.timer else None, + } + + return Samples( + model=model, + sample_list=sample_list, + samples_info=samples_info, + ) + + +class BFGS(AbstractBFGS): + """ + The BFGS non-linear search, which wraps the scipy Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm. + + See the docstrings of the `BFGS` and `LBFGS` classes for a description of the arguments of this class. + + For a full description of the scipy BFGS method, checkout its documentation: + + https://docs.scipy.org/doc/scipy/reference/optimize.minimize-bfgs.html#optimize-minimize-bfgs + + If you use `BFGS` as part of a published work, please cite the package via scipy following the instructions + under the *Attribution* section of the GitHub page. + + By default, the BFGS method scipy implementation does not store the history of parameter values and + log likelihood values during the non-linear search. This is because storing these values can require a large + amount of memory, in contradiction to the BFGS method's primary advantage of being memory efficient. + This means that it is difficult to visualize the BFGS method results (e.g. log likelihood vs iteration). + + **PyAutoFit** extends the class with the option of using visualize mode, which stores the history of parameter + values and log likelihood values during the non-linear search. This allows the results of the BFGS method to be + visualized after the search has completed, and it is enabled by setting the `visualize` flag to `True`. + + Parameters + ---------- + name + The name of the search, controlling the last folder results are output. + path_prefix + The path of folders prefixing the name folder where results are output. + unique_tag + The name of a unique tag for this model-fit, which will be given a unique entry in the sqlite database + and also acts as the folder after the path prefix and before the search name. + initializer + Generates the initialize samples of non-linear parameter space (see autofit.non_linear.initializer). + number_of_cores: int + The number of cores sampling is performed using a Python multiprocessing Pool instance. + session + An SQLalchemy session instance so the results of the model-fit are written to an SQLite database. + visualize + If True, visualization of the search is enabled, which requires storing the history of parameter values and + log likelihood values during the non-linear search. + """ + + method = "BFGS" + +class LBFGS(AbstractBFGS): + """ + The L-BFGS non-linear search, which wraps the scipy Limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) + algorithm. + + See the docstrings of the `BFGS` and `LBFGS` classes for a description of the arguments of this class. + + For a full description of the scipy L-BFGS method, checkout its documentation: + + https://docs.scipy.org/doc/scipy/reference/optimize.minimize-lbfgsb.html + + If you use `LBFGS` as part of a published work, please cite the package via scipy following the instructions + under the *Attribution* section of the GitHub page. + + By default, the L-BFGS method scipy implementation does not store the history of parameter values and + log likelihood values during the non-linear search. This is because storing these values can require a large + amount of memory, in contradiction to the L-BFGS method's primary advantage of being memory efficient. + This means that it is difficult to visualize the L-BFGS method results (e.g. log likelihood vs iteration). + + **PyAutoFit** extends the class with the option of using visualize mode, which stores the history of parameter + values and log likelihood values during the non-linear search. This allows the results of the L-BFGS method to be + visualized after the search has completed, and it is enabled by setting the `visualize` flag to `True`. + + Parameters + ---------- + name + The name of the search, controlling the last folder results are output. + path_prefix + The path of folders prefixing the name folder where results are output. + unique_tag + The name of a unique tag for this model-fit, which will be given a unique entry in the sqlite database + and also acts as the folder after the path prefix and before the search name. + initializer + Generates the initialize samples of non-linear parameter space (see autofit.non_linear.initializer). + number_of_cores: int + The number of cores sampling is performed using a Python multiprocessing Pool instance. + session + An SQLalchemy session instance so the results of the model-fit are written to an SQLite database. + visualize + If True, visualization of the search is enabled, which requires storing the history of parameter values and + log likelihood values during the non-linear search. + """ + + method = "L-BFGS-B" \ No newline at end of file diff --git a/autofit/non_linear/search/optimize/lbfgs/__init__.py b/autofit/non_linear/search/mle/drawer/__init__.py similarity index 100% rename from autofit/non_linear/search/optimize/lbfgs/__init__.py rename to autofit/non_linear/search/mle/drawer/__init__.py diff --git a/autofit/non_linear/search/optimize/drawer/search.py b/autofit/non_linear/search/mle/drawer/search.py similarity index 95% rename from autofit/non_linear/search/optimize/drawer/search.py rename to autofit/non_linear/search/mle/drawer/search.py index 43c481740..aa6567ee0 100644 --- a/autofit/non_linear/search/optimize/drawer/search.py +++ b/autofit/non_linear/search/mle/drawer/search.py @@ -5,12 +5,12 @@ from autofit.mapper.prior_model.abstract import AbstractPriorModel from autofit.non_linear.fitness import Fitness -from autofit.non_linear.search.optimize.abstract_optimize import AbstractOptimizer +from autofit.non_linear.search.mle.abstract_mle import AbstractMLE from autofit.non_linear.initializer import AbstractInitializer from autofit.non_linear.samples import Samples, Sample -class Drawer(AbstractOptimizer): +class Drawer(AbstractMLE): __identifier_fields__ = ("total_draws",) def __init__( diff --git a/autofit/non_linear/search/optimize/pyswarms/__init__.py b/autofit/non_linear/search/mle/pyswarms/__init__.py similarity index 100% rename from autofit/non_linear/search/optimize/pyswarms/__init__.py rename to autofit/non_linear/search/mle/pyswarms/__init__.py diff --git a/autofit/non_linear/search/optimize/pyswarms/search/__init__.py b/autofit/non_linear/search/mle/pyswarms/search/__init__.py similarity index 100% rename from autofit/non_linear/search/optimize/pyswarms/search/__init__.py rename to autofit/non_linear/search/mle/pyswarms/search/__init__.py diff --git a/autofit/non_linear/search/optimize/pyswarms/search/abstract.py b/autofit/non_linear/search/mle/pyswarms/search/abstract.py similarity index 92% rename from autofit/non_linear/search/optimize/pyswarms/search/abstract.py rename to autofit/non_linear/search/mle/pyswarms/search/abstract.py index 7c0acb7cd..0e8c18490 100644 --- a/autofit/non_linear/search/optimize/pyswarms/search/abstract.py +++ b/autofit/non_linear/search/mle/pyswarms/search/abstract.py @@ -7,7 +7,7 @@ from autofit.mapper.prior_model.abstract import AbstractPriorModel from autofit.non_linear.fitness import Fitness from autofit.non_linear.initializer import AbstractInitializer -from autofit.non_linear.search.optimize.abstract_optimize import AbstractOptimizer +from autofit.non_linear.search.mle.abstract_mle import AbstractMLE from autofit.non_linear.samples.sample import Sample from autofit.non_linear.samples.samples import Samples @@ -66,7 +66,7 @@ def __call__(self, parameters, *kwargs): return np.asarray(figure_of_merit_list) -class AbstractPySwarms(AbstractOptimizer): +class AbstractPySwarms(AbstractMLE): def __init__( self, name: Optional[str] = None, @@ -79,7 +79,7 @@ def __init__( **kwargs ): """ - A PySwarms Particle Swarm Optimizer global non-linear search. + A PySwarms Particle Swarm MLE global non-linear search. For a full description of PySwarms, checkout its Github and readthedocs webpages: @@ -215,9 +215,10 @@ def _fit(self, model: AbstractPriorModel, analysis): total_iterations += iterations - self.paths.save_search_internal( - obj=search_internal, - ) + # TODO : Running PySwarms in NoteBook raises + # TODO: TypeError: cannot pickle '_hashlib.HMAC' object + + self.output_search_internal(search_internal=search_internal) self.perform_update( model=model, @@ -230,6 +231,14 @@ def _fit(self, model: AbstractPriorModel, analysis): return search_internal + def output_search_internal(self, search_internal): + try: + self.paths.save_search_internal( + obj=search_internal, + ) + except TypeError: + pass + def samples_via_internal_from(self, model, search_internal=None): """ Returns a `Samples` object from the pyswarms internal results. diff --git a/autofit/non_linear/search/optimize/pyswarms/search/globe.py b/autofit/non_linear/search/mle/pyswarms/search/globe.py similarity index 92% rename from autofit/non_linear/search/optimize/pyswarms/search/globe.py rename to autofit/non_linear/search/mle/pyswarms/search/globe.py index ab16dadb5..ebf620cc1 100644 --- a/autofit/non_linear/search/optimize/pyswarms/search/globe.py +++ b/autofit/non_linear/search/mle/pyswarms/search/globe.py @@ -2,7 +2,7 @@ from autofit.database.sqlalchemy_ import sa from autofit.non_linear.initializer import AbstractInitializer -from autofit.non_linear.search.optimize.pyswarms.search.abstract import AbstractPySwarms +from autofit.non_linear.search.mle.pyswarms.search.abstract import AbstractPySwarms class PySwarmsGlobal(AbstractPySwarms): @@ -25,7 +25,7 @@ def __init__( **kwargs ): """ - A PySwarms Particle Swarm Optimizer global non-linear search. + A PySwarms Particle Swarm MLE global non-linear search. For a full description of PySwarms, checkout its Github and readthedocs webpages: diff --git a/autofit/non_linear/search/optimize/pyswarms/search/local.py b/autofit/non_linear/search/mle/pyswarms/search/local.py similarity index 92% rename from autofit/non_linear/search/optimize/pyswarms/search/local.py rename to autofit/non_linear/search/mle/pyswarms/search/local.py index ff1102a26..4432eae8a 100644 --- a/autofit/non_linear/search/optimize/pyswarms/search/local.py +++ b/autofit/non_linear/search/mle/pyswarms/search/local.py @@ -1,7 +1,7 @@ from typing import Optional from autofit.database.sqlalchemy_ import sa -from autofit.non_linear.search.optimize.pyswarms.search.abstract import AbstractPySwarms +from autofit.non_linear.search.mle.pyswarms.search.abstract import AbstractPySwarms class PySwarmsLocal(AbstractPySwarms): @@ -25,7 +25,7 @@ def __init__( **kwargs ): """ - A PySwarms Particle Swarm Optimizer global non-linear search. + A PySwarms Particle Swarm MLE global non-linear search. For a full description of PySwarms, checkout its Github and readthedocs webpages: diff --git a/autofit/non_linear/search/nest/abstract_nest.py b/autofit/non_linear/search/nest/abstract_nest.py index 307f06837..5e356891d 100644 --- a/autofit/non_linear/search/nest/abstract_nest.py +++ b/autofit/non_linear/search/nest/abstract_nest.py @@ -8,7 +8,7 @@ from autofit.non_linear.initializer import ( InitializerPrior, AbstractInitializer, - SpecificRangeInitializer, + InitializerParamBounds, ) from autofit.non_linear.samples import SamplesNest from autofit.non_linear.plot.nest_plotters import NestPlotter @@ -44,9 +44,9 @@ def __init__( session An SQLAlchemy session instance so the results of the model-fit are written to an SQLite database. """ - if isinstance(initializer, SpecificRangeInitializer): + if isinstance(initializer, InitializerParamBounds): raise ValueError( - "SpecificRangeInitializer cannot be used for nested sampling" + "InitializerParamBounds cannot be used for nested sampling" ) super().__init__( diff --git a/autofit/non_linear/search/nest/nautilus/search.py b/autofit/non_linear/search/nest/nautilus/search.py index 603236dc2..289f1a010 100644 --- a/autofit/non_linear/search/nest/nautilus/search.py +++ b/autofit/non_linear/search/nest/nautilus/search.py @@ -290,6 +290,23 @@ def call_search(self, search_internal, model, analysis): finished = False + minimum_iterations_per_updates = 3 * self.config_dict_search["n_live"] + if self.iterations_per_update < minimum_iterations_per_updates: + + self.iterations_per_update = minimum_iterations_per_updates + + logger.info( + f""" + The number of iterations_per_update is less than 3 times the number of live points, which can cause + issues where Nautilus loses sampling information due to stopping to output results. The number of + iterations per update has been increased to 3 times the number of live points, therefore a value + of {minimum_iterations_per_updates}. + + To remove this warning, increase the number of iterations_per_update to three or more times the + number of live points. + """ + ) + while not finished: iterations, total_iterations = self.iterations_from( diff --git a/autofit/non_linear/search/optimize/abstract_optimize.py b/autofit/non_linear/search/optimize/abstract_optimize.py deleted file mode 100644 index ad8a7db30..000000000 --- a/autofit/non_linear/search/optimize/abstract_optimize.py +++ /dev/null @@ -1,31 +0,0 @@ -from abc import ABC - -from autoconf import conf -from autofit.non_linear.search.abstract_search import NonLinearSearch -from autofit.non_linear.samples import Samples -from autofit.non_linear.plot.optimize_plotters import OptimizePlotter -from autofit.non_linear.plot.output import Output - - -class AbstractOptimizer(NonLinearSearch, ABC): - @property - def config_type(self): - return conf.instance["non_linear"]["optimize"] - - @property - def samples_cls(self): - return Samples - - @property - def plotter_cls(self): - return OptimizePlotter - - def plot_results(self, samples): - - def should_plot(name): - return conf.instance["visualize"]["plots_search"]["optimize"][name] - - plotter = OptimizePlotter( - samples=samples, - output=Output(path=self.paths.image_path / "search", format="png") - ) \ No newline at end of file diff --git a/autofit/non_linear/search/optimize/lbfgs/search.py b/autofit/non_linear/search/optimize/lbfgs/search.py deleted file mode 100644 index 963160d0e..000000000 --- a/autofit/non_linear/search/optimize/lbfgs/search.py +++ /dev/null @@ -1,240 +0,0 @@ -from typing import Optional - -from autoconf import cached_property -from autofit.database.sqlalchemy_ import sa - -from autofit.mapper.prior_model.abstract import AbstractPriorModel -from autofit.non_linear.search.optimize.abstract_optimize import AbstractOptimizer -from autofit.non_linear.analysis import Analysis -from autofit.non_linear.fitness import Fitness -from autofit.non_linear.initializer import AbstractInitializer -from autofit.non_linear.samples.sample import Sample -from autofit.non_linear.samples.samples import Samples - -import copy -from scipy import optimize -import numpy as np - - -class LBFGS(AbstractOptimizer): - __identifier_fields__ = () - - def __init__( - self, - name: Optional[str] = None, - path_prefix: Optional[str] = None, - unique_tag: Optional[str] = None, - initializer: Optional[AbstractInitializer] = None, - iterations_per_update: int = None, - session: Optional[sa.orm.Session] = None, - **kwargs - ): - """ - A L-BFGS scipy non-linear search. - - For a full description of the scipy L-BFGS method, checkout its documentation: - - https://docs.scipy.org/doc/scipy/reference/optimize.minimize-lbfgsb.html - - If you use `LBFGS` as part of a published work, please cite the package via scipy following the instructions - under the *Attribution* section of the GitHub page. - - Parameters - ---------- - name - The name of the search, controlling the last folder results are output. - path_prefix - The path of folders prefixing the name folder where results are output. - unique_tag - The name of a unique tag for this model-fit, which will be given a unique entry in the sqlite database - and also acts as the folder after the path prefix and before the search name. - initializer - Generates the initialize samples of non-linear parameter space (see autofit.non_linear.initializer). - number_of_cores: int - The number of cores sampling is performed using a Python multiprocessing Pool instance. - session - An SQLalchemy session instance so the results of the model-fit are written to an SQLite database. - """ - - super().__init__( - name=name, - path_prefix=path_prefix, - unique_tag=unique_tag, - initializer=initializer, - iterations_per_update=iterations_per_update, - session=session, - **kwargs - ) - - self.logger.debug("Creating LBFGS Search") - - @cached_property - def config_dict_options(self): - config_dict = copy.deepcopy(self._class_config["options"]) - - for key, value in config_dict.items(): - try: - config_dict[key] = self.kwargs[key] - except KeyError: - pass - - return config_dict - - def _fit( - self, - model: AbstractPriorModel, - analysis: Analysis, - ): - """ - Fit a model using the scipy L-BFGS method and the Analysis class which contains the data and returns the log - likelihood from instances of the model, which the `NonLinearSearch` seeks to maximize. - - Parameters - ---------- - model - The model which generates instances for different points in parameter space. - analysis - Contains the data and the log likelihood function which fits an instance of the model to the data, - returning the log likelihood the `NonLinearSearch` maximizes. - - Returns - ------- - A result object comprising the Samples object that inclues the maximum log likelihood instance and full - chains used by the fit. - """ - fitness = Fitness( - model=model, - analysis=analysis, - paths=self.paths, - fom_is_log_likelihood=False, - resample_figure_of_merit=-np.inf, - convert_to_chi_squared=True, - ) - - try: - search_internal_dict = self.paths.load_search_internal() - - x0 = search_internal_dict["x0"] - total_iterations = search_internal_dict["total_iterations"] - - self.logger.info( - "Resuming LBFGS non-linear search (previous samples found)." - ) - - except (FileNotFoundError, TypeError): - ( - unit_parameter_lists, - parameter_lists, - log_posterior_list, - ) = self.initializer.samples_from_model( - total_points=1, - model=model, - fitness=fitness, - paths=self.paths, - n_cores=self.number_of_cores, - ) - - x0 = np.asarray(parameter_lists[0]) - - total_iterations = 0 - - self.logger.info( - "Starting new LBFGS non-linear search (no previous samples found)." - ) - - maxiter = self.config_dict_options.get("maxiter", 1e8) - - while total_iterations < maxiter: - iterations_remaining = maxiter - total_iterations - - iterations = min(self.iterations_per_update, iterations_remaining) - - if iterations > 0: - config_dict_options = self.config_dict_options - config_dict_options["maxiter"] = iterations - - search_internal = optimize.minimize( - fun=fitness.__call__, - x0=x0, - method="L-BFGS-B", - options=config_dict_options, - **self.config_dict_search - ) - - total_iterations += search_internal.nit - - search_internal.log_posterior_list = -0.5 * fitness( - parameters=search_internal.x - ) - - self.paths.save_search_internal( - obj=search_internal, - ) - - x0 = search_internal.x - - if search_internal.nit < iterations: - return search_internal - - self.perform_update( - model=model, - analysis=analysis, - during_analysis=True, - search_internal=search_internal, - ) - - self.logger.info("L-BFGS sampling complete.") - - return search_internal - - def samples_via_internal_from( - self, model: AbstractPriorModel, search_internal=None - ): - """ - Returns a `Samples` object from the LBFGS internal results. - - The samples contain all information on the parameter space sampling (e.g. the parameters, - log likelihoods, etc.). - - The internal search results are converted from the native format used by the search to lists of values - (e.g. `parameter_lists`, `log_likelihood_list`). - - Parameters - ---------- - model - Maps input vectors of unit parameter values to physical values and model instances via priors. - """ - - if search_internal is None: - search_internal = self.paths.load_search_internal() - - x0 = search_internal.x - total_iterations = search_internal.nit - log_posterior_list = np.array([search_internal.log_posterior_list]) - - parameter_lists = [list(x0)] - log_prior_list = model.log_prior_list_from(parameter_lists=parameter_lists) - - log_likelihood_list = [ - lp - prior for lp, prior in zip(log_posterior_list, log_prior_list) - ] - weight_list = len(log_likelihood_list) * [1.0] - - sample_list = Sample.from_lists( - model=model, - parameter_lists=parameter_lists, - log_likelihood_list=log_likelihood_list, - log_prior_list=log_prior_list, - weight_list=weight_list, - ) - - samples_info = { - "total_iterations": total_iterations, - "time": self.timer.time if self.timer else None, - } - - return Samples( - model=model, - sample_list=sample_list, - samples_info=samples_info, - ) diff --git a/autofit/plot/__init__.py b/autofit/plot/__init__.py index ab35eb913..525a81891 100644 --- a/autofit/plot/__init__.py +++ b/autofit/plot/__init__.py @@ -1,5 +1,5 @@ from autofit.non_linear.plot.samples_plotters import SamplesPlotter from autofit.non_linear.plot.mcmc_plotters import MCMCPlotter -from autofit.non_linear.plot.optimize_plotters import OptimizePlotter +from autofit.non_linear.plot.mle_plotters import MLEPlotter from autofit.non_linear.plot.nest_plotters import NestPlotter from autofit.non_linear.plot.output import Output \ No newline at end of file diff --git a/autofit/text/formatter.py b/autofit/text/formatter.py index 7508f3974..751ebff23 100644 --- a/autofit/text/formatter.py +++ b/autofit/text/formatter.py @@ -1,6 +1,8 @@ import csv import logging -from typing import Tuple +from typing import Tuple, Union + +from pathlib import Path from autoconf import conf from autofit.tools.util import open_ @@ -221,7 +223,7 @@ def output_list_of_strings_to_file(file, list_of_strings): f.write("".join(list_of_strings)) -def write_table(headers, rows, filename: str): +def write_table(headers, rows, filename: Union[str, Path]): """ Write a table of parameters, posteriors, priors and likelihoods. diff --git a/docs/api/plot.rst b/docs/api/plot.rst index 84055e406..ac381ca1e 100644 --- a/docs/api/plot.rst +++ b/docs/api/plot.rst @@ -24,4 +24,4 @@ Plotters NestPlotter MCMCPlotter - OptimizePlotter \ No newline at end of file + MLEPlotter \ No newline at end of file diff --git a/docs/api/searches.rst b/docs/api/searches.rst index 6c08e98db..6d070bf86 100644 --- a/docs/api/searches.rst +++ b/docs/api/searches.rst @@ -4,8 +4,8 @@ Non-Linear Searches A non-linear search is an algorithm which fits a model to data. -**PyAutoFit** currently supports three types of non-linear search algorithms: nested samplers, -Markov Chain Monte Carlo (MCMC) and optimizers. +**PyAutoFit** currently supports three types of non-linear search algorithms: nested samplers (nest), +Markov Chain Monte Carlo (MCMC) and Maximum Likelihood Estimators (MLE). **Examples / Tutorials:** @@ -41,8 +41,8 @@ MCMC Emcee Zeus -Optimizers ----------- +Maximum Likelihood Estimators +----------------------------- .. currentmodule:: autofit @@ -51,6 +51,8 @@ Optimizers :template: custom-class-template.rst :recursive: + BFGS + LBFGS PySwarmsLocal PySwarmsGlobal diff --git a/docs/conf.py b/docs/conf.py index a77132aaf..fcb18313a 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -3,7 +3,7 @@ # # This file only contains a selection of the most common options. For a full # list see the documentation: -# https://www.sphinx-doc.org/en/master/usage/configuration.html +# https://www.sphinx-doc.org/en/main/usage/configuration.html # -- Path setup -------------------------------------------------------------- @@ -73,7 +73,7 @@ intersphinx_mapping = { "python": ("https://docs.python.org/3", None), - "sphinx": ("https://www.sphinx-doc.org/en/master", None), + "sphinx": ("https://www.sphinx-doc.org/en/main", None), } # -- Options for TODOs ------------------------------------------------------- diff --git a/docs/cookbooks/analysis.rst b/docs/cookbooks/analysis.rst index 8b766c3d6..de6ce6d0d 100644 --- a/docs/cookbooks/analysis.rst +++ b/docs/cookbooks/analysis.rst @@ -641,7 +641,7 @@ These files can then also be loaded via the database, as described in the databa Parameters ---------- paths - The PyAutoFit paths object which manages all paths, e.g. where + The paths object which manages all paths, e.g. where the non-linear search outputs are stored, visualization, and the pickled objects used by the aggregator output by this function. """ @@ -674,7 +674,7 @@ These files can then also be loaded via the database, as described in the databa Parameters ---------- paths - The PyAutoFit paths object which manages all paths, e.g. where the + The paths object which manages all paths, e.g. where the non-linear search outputs are stored, visualization and the pickled objects used by the aggregator output by this function. result diff --git a/docs/cookbooks/model.rst b/docs/cookbooks/model.rst index 1f9081d75..84590ccca 100644 --- a/docs/cookbooks/model.rst +++ b/docs/cookbooks/model.rst @@ -10,6 +10,8 @@ This cookbook provides an overview of basic model composition tools. **Contents:** +**Models:** + If first describes how to use the ``af.Model`` object to define models with a single model component from single Python classes, with the following sections: @@ -21,6 +23,8 @@ Python classes, with the following sections: - **Tuple Parameters (Model)**: Defining model components with parameters that are tuples. - **Json Output (Model)**: Output a model in human readable text via a .json file and loading it back again. +**Collections:** + It then describes how to use the ``af.Collection`` object to define models with many model components from multiple Python classes, with the following sections: @@ -31,6 +35,19 @@ Python classes, with the following sections: - **Json Output (Collection)**: Output a collection in human readable text via a .json file and loading it back again. - **Extensible Models (Collection)**: Using collections to extend models with new model components, including the use of Python dictionaries and lists. +**Arrays:** + +The cookbook next describes using NumPy arrays via tbe `af.Array` object to compose models, where each entry of the +array is a free parameters, therefore offering maximum flexibility with the number of free parameter. This has +the following sections: + + - **Model Composition (af.Array)**: Composing models using NumPy arrays and `af.Array`(). + - **Prior Customization (af.Array)**: How to customize the priors of a numpy array model. + - **Instances (af.Array)**: Create an instance of a numpy array model via input parameters. + - **Model Customization (af.Array):** Customize a numpy array model (e.g. fixing parameters or linking them to one another). + - **Json Output (af.Array)**: Output a numpy array model in human readable text via a .json file and loading it back again. + - **Extensible Models (af.Array)**: Using numpy arrays to compose models with a flexible number of parameters. + Python Class Template --------------------- @@ -718,7 +735,7 @@ Python dictionaries can easily be saved to hard disk as a ``.json`` file. This means we can save any **PyAutoFit** model to hard-disk. -Checkout the file ``autofit_workspace/*/model/jsons/model.json`` to see the model written as a .json. +Checkout the file ``autofit_workspace/*/model/jsons/collection.json`` to see the model written as a .json. .. code-block:: python @@ -910,6 +927,321 @@ This gives the following output: normalization (Gaussian) = 5.0 sigma (Gaussian) = 6.0 +Model Composition (af.Array) +---------------------------- + +Models can be composed using NumPy arrays, where each element of the array is a free parameter. + +This offers a lot more flexibility than using ``Model`` and ``Collection`` objects, as the number of parameters in the +model is chosen on initialization via the input of the ``shape`` attribute. + +For many use cases, this flexibility is key to ensuring model composition is as easy as possible, for example when +a part of the model being fitted is a matrix of parameters which may change shape depending on the dataset being +fitted. + +To compose models using NumPy arrays, we use the ``af.Array`` object. + +.. code-block:: python + + model = af.Array( + shape=(2, 2), + prior=af.GaussianPrior(mean=0.0, sigma=1.0), + ) + +Each element of the array is a free parameter, which for ``shape=(2,2)`` means the model has 4 free parameters. + +.. code-block:: python + + print(f"Model Total Free Parameters = {model.total_free_parameters}") + +The ``info`` attribute of the model gives information on all of the parameters and their priors. + +.. code-block:: python + + print(model.info) + +This gives the following output: + +.. code-block:: bash + + Total Free Parameters = 4 + + model Array (N=4) + indices list (N=0) + + shape (2, 2) + indices + 0 (0, 0) + 1 (0, 1) + 2 (1, 0) + 3 (1, 1) + prior_0_0 GaussianPrior [124], mean = 0.0, sigma = 1.0 + prior_0_1 GaussianPrior [125], mean = 0.0, sigma = 1.0 + prior_1_0 GaussianPrior [126], mean = 0.0, sigma = 1.0 + prior_1_1 GaussianPrior [127], mean = 0.0, sigma = 1.0 + +Prior Customization (af.Array) +------------------------------ + +The prior of every parameter in the array is set via the ``prior`` input above. + +NumPy array models do not currently support default priors via config files, so all priors must be manually specified. + +The prior of every parameter in the array can be customized by normal NumPy array indexing: + +.. code-block:: python + + model = af.Array(shape=(2, 2), prior=af.GaussianPrior(mean=0.0, sigma=1.0)) + + model.array[0, 0] = af.UniformPrior(lower_limit=0.0, upper_limit=1.0) + model.array[0, 1] = af.LogUniformPrior(lower_limit=1e-4, upper_limit=1e4) + model.array[1, 0] = af.GaussianPrior(mean=0.0, sigma=2.0) + +The ``info`` attribute shows the customized priors. + +.. code-block:: python + + print(model.info) + +The output is as follows: + +.. code-block:: bash + + Total Free Parameters = 4 + + model Array (N=4) + indices list (N=0) + + shape (2, 2) + indices + 0 (0, 0) + 1 (0, 1) + 2 (1, 0) + 3 (1, 1) + prior_0_0 UniformPrior [133], lower_limit = 0.0, upper_limit = 1.0 + prior_0_1 LogUniformPrior [134], lower_limit = 0.0001, upper_limit = 10000.0 + prior_1_0 GaussianPrior [135], mean = 0.0, sigma = 2.0 + prior_1_1 GaussianPrior [132], mean = 0.0, sigma = 1.0 + +Instances (af.Array) +-------------------- + +Instances of numpy array model components can be created, where an input ``vector`` of parameters is mapped to create +an instance of the Python class of the model. + +If the priors of the numpy array are not customized, ordering of parameters goes from element [0,0] to [0,1] to [1,0], +as shown by the ``paths`` attribute. + +.. code-block:: python + + model = af.Array( + shape=(2, 2), + prior=af.GaussianPrior(mean=0.0, sigma=1.0), + ) + + print(model.paths) + +The output is as follows: + +.. code-block:: bash + + ['prior_0_0', 'prior_0_1', 'prior_1_0', 'prior_1_1'] + +An instance can then be created by passing a vector of parameters to the model via the ``instance_from_vector`` method. + +The ``instance`` created is a NumPy array, where each element is the value passed in the vector. + +.. code-block:: python + + instance = model.instance_from_vector(vector=[0.0, 1.0, 2.0, 3.0]) + + print("\nModel Instance:") + print(instance) + +The output is as follows: + +.. code-block:: bash + + Model Instance: + [[0. 1.] + [2. 3.]] + +Prior customization changes the order of the parameters, therefore if you customize the priors of the numpy +array you must check the ordering of the parameters in the ``paths`` attribute before passing a vector to +the ``instance_from_vector`` + + +.. code-block:: python + + model[0, 0] = af.UniformPrior(lower_limit=0.0, upper_limit=1.0) + model[0, 1] = af.LogUniformPrior(lower_limit=1e-4, upper_limit=1e4) + model[1, 0] = af.GaussianPrior(mean=0.0, sigma=2.0) + + print(model.paths) + +The output is as follows: + +.. code-block:: bash + + [('prior_1_1',), ('prior_0_0',), ('prior_0_1',), ('prior_1_0',)] + +If we create a vector and print its values from this customized model: + +.. code-block:: python + + instance = model.instance_from_vector(vector=[0.0, 1.0, 2.0, 3.0]) + + print("\nModel Instance:") + print(instance) + +The output is as follows: + +.. code-block:: bash + + Model Instance: + [[1. 2.] + [3. 0.]] + +Model Customization (af.Array) +------------------------------ + +The model customization API for numpy array models is the same as for ``af.Model`` and ``af.Collection`` objects. + +.. code-block:: python + + model = af.Array( + shape=(2, 2), + prior=af.GaussianPrior(mean=0.0, sigma=1.0), + ) + + model[0,0] = 50.0 + model[0,1] = model[1,0] + model.add_assertion(model[1,1] > 0.0) + + print(model.info) + +The output is as follows: + +.. code-block:: bash + Total Free Parameters = 2 + + model Array (N=2) + indices list (N=0) + + shape (2, 2) + indices + 0 (0, 0) + 1 (0, 1) + 2 (1, 0) + 3 (1, 1) + prior_0_0 50.0 + prior_0_1 - prior_1_0 GaussianPrior [147], mean = 0.0, sigma = 1.0 + prior_1_1 GaussianPrior [148], mean = 0.0, sigma = 1.0 + + +JSon Outputs (af.Array) +------------------------ + +An ``Array`` has a ``dict`` attribute, which express all information about the model as a Python dictionary. + +By printing this dictionary we can therefore get a concise summary of the model. + +.. code-block:: python + + model = af.Array( + shape=(2, 2), + prior=af.GaussianPrior(mean=0.0, sigma=1.0), + ) + + print(model.dict()) + +Python dictionaries can easily be saved to hard disk as a ``.json`` file. + +This means we can save any **PyAutoFit** model to hard-disk. + +Checkout the file ``autofit_workspace/*/model/jsons/array.json`` to see the model written as a .json. + +.. code-block:: python + + model_path = path.join("scripts", "model", "jsons") + + os.makedirs(model_path, exist_ok=True) + + model_file = path.join(model_path, "array.json") + + with open(model_file, "w+") as f: + json.dump(model.dict(), f, indent=4) + +We can load the model from its ``.json`` file, meaning that one can easily save a model to hard disk and load it +elsewhere. + +.. code-block:: python + + model = af.Array.from_json(file=model_file) + + print(f"\n Model via Json Prior Count = {model.prior_count}") + +Extensible Models (af.Array) +---------------------------- + +For ``Model`` objects, the number of parameters is fixed to those listed in the input Python class when the model is +created. + +For ``Collection`` objects, the use of dictionaries and lists allows for the number of parameters to be extended, but it +was still tied to the input Python classes when the model was created. + +For ``Array`` objects, the number of parameters is fully customizable, you choose the shape of the array and therefore +the number of parameters in the model when you create it. + +This makes ``Array`` objects the most extensible and flexible way to compose models. + +You can also combine ``Array`` objects with ``Collection`` objects to create models with a mix of fixed and extensible +parameters. + +.. code-block:: python + + model = af.Collection( + gaussian=Gaussian, + array=af.Array(shape=(3, 2), prior=af.GaussianPrior(mean=0.0, sigma=1.0)) + ) + + model.gaussian.sigma = 2.0 + model.array[0, 0] = 1.0 + + print(model.info) + +The output is as follows: + +.. code-block:: python + + Total Free Parameters = 7 + + model Collection (N=7) + gaussian Gaussian (N=2) + array Array (N=5) + indices list (N=0) + + gaussian + centre UniformPrior [165], lower_limit = 0.0, upper_limit = 100.0 + normalization LogUniformPrior [166], lower_limit = 1e-06, upper_limit = 1000000.0 + sigma 2.0 + array + shape (3, 2) + indices + 0 (0, 0) + 1 (0, 1) + 2 (1, 0) + 3 (1, 1) + 4 (2, 0) + 5 (2, 1) + prior_0_0 1.0 + prior_0_1 GaussianPrior [160], mean = 0.0, sigma = 1.0 + prior_1_0 GaussianPrior [161], mean = 0.0, sigma = 1.0 + prior_1_1 GaussianPrior [162], mean = 0.0, sigma = 1.0 + prior_2_0 GaussianPrior [163], mean = 0.0, sigma = 1.0 + prior_2_1 GaussianPrior [164], mean = 0.0, sigma = 1.0 + + Wrap Up ------- diff --git a/docs/cookbooks/result.rst b/docs/cookbooks/result.rst index ccbb3f56e..7102e759f 100644 --- a/docs/cookbooks/result.rst +++ b/docs/cookbooks/result.rst @@ -553,7 +553,7 @@ as 1D numpy arrays, are converted to a suitable dictionary output format. This u Parameters ---------- paths - The PyAutoFit paths object which manages all paths, e.g. where the non-linear search outputs are stored, + The paths object which manages all paths, e.g. where the non-linear search outputs are stored, visualization, and the pickled objects used by the aggregator output by this function. """ from autoconf.dictable import to_dict @@ -575,7 +575,7 @@ as 1D numpy arrays, are converted to a suitable dictionary output format. This u Parameters ---------- paths - The PyAutoFit paths object which manages all paths, e.g. where the non-linear search outputs are stored, + The paths object which manages all paths, e.g. where the non-linear search outputs are stored, visualization and the pickled objects used by the aggregator output by this function. result The result of a model fit, including the non-linear search, samples and maximum likelihood model. diff --git a/docs/cookbooks/search.rst b/docs/cookbooks/search.rst index 8bb92caaf..32bada27f 100644 --- a/docs/cookbooks/search.rst +++ b/docs/cookbooks/search.rst @@ -241,7 +241,7 @@ We now define the start point of certain parameters in the model as follows. .. code-block:: python - initializer = af.SpecificRangeInitializer( + initializer = af.InitializerParamBounds( { model.centre: (49.0, 51.0), model.normalization: (4.0, 6.0), diff --git a/docs/index.rst b/docs/index.rst index a1ac9d908..58e29d7b9 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -22,7 +22,7 @@ The following links are useful for new starters: - `The introduction Jupyter Notebook on Binder `_, where you can try **PyAutoFit** in a web browser (without installation). -- `The autofit_workspace GitHub repository `_, which includes example scripts and the `HowToFit Jupyter notebook lectures `_ which give new users a step-by-step introduction to **PyAutoFit**. +- `The autofit_workspace GitHub repository `_, which includes example scripts and the `HowToFit Jupyter notebook lectures `_ which give new users a step-by-step introduction to **PyAutoFit**. Support ------- @@ -49,7 +49,7 @@ API Overview To illustrate the **PyAutoFit** API, we use an illustrative toy model of fitting a one-dimensional Gaussian to noisy 1D data. Here's the ``data`` (black) and the model (red) we'll fit: -.. image:: https://raw.githubusercontent.com/rhayes777/PyAutoFit/master/files/toy_model_fit.png +.. image:: https://raw.githubusercontent.com/rhayes777/PyAutoFit/main/files/toy_model_fit.png :width: 400 We define our model, a 1D Gaussian by writing a Python class using the format below: diff --git a/docs/installation/overview.rst b/docs/installation/overview.rst index c0c2c47d8..ad3853de2 100644 --- a/docs/installation/overview.rst +++ b/docs/installation/overview.rst @@ -3,7 +3,7 @@ Overview ======== -**PyAutoFit** requires Python 3.8 - 3.11 and support the Linux, MacOS and Windows operating systems. +**PyAutoFit** requires Python 3.9 - 3.12 and support the Linux, MacOS and Windows operating systems. **PyAutoFit** can be installed via the Python distribution `Anaconda `_ or using `Pypi `_ to ``pip install`` **PyAutoFit** into your Python distribution. diff --git a/docs/overview/scientific_workflow.rst b/docs/overview/scientific_workflow.rst index f6f1ed9be..ff670803a 100644 --- a/docs/overview/scientific_workflow.rst +++ b/docs/overview/scientific_workflow.rst @@ -525,7 +525,7 @@ For simpler scenarios, adjustments might include: In more intricate cases, models might involve numerous parameters and complex compositions of multiple model components. **PyAutoFit** offers a sophisticated model composition API designed to handle these complexities. It provides -tools for constructing elaborate models using lists of Python classes and hierarchical structures of Python classes. +tools for constructing elaborate models using lists of Python classes, NumPy arrays and hierarchical structures of Python classes. For a detailed exploration of these capabilities, you can refer to the `model cookbook `_, which provides comprehensive diff --git a/docs/overview/the_basics.rst b/docs/overview/the_basics.rst index a67e9302c..da3d213d5 100644 --- a/docs/overview/the_basics.rst +++ b/docs/overview/the_basics.rst @@ -298,13 +298,13 @@ Analysis We now tell **PyAutoFit** how to fit the model to the data. -We define an `Analysis` class, which includes: +We define an ``Analysis`` class, which includes: -- An `__init__` constructor that takes `data` and `noise_map` as inputs (this can be extended with additional elements necessary for fitting the model to the data). +- An ``__init__`` constructor that takes ``data`` and ``noise_map`` as inputs (this can be extended with additional elements necessary for fitting the model to the data). -- A `log_likelihood_function` that defines how to fit an `instance` of the model to the data and return a log likelihood value. +- A ``log_likelihood_function`` that defines how to fit an ``instance`` of the model to the data and return a log likelihood value. -Read the comments and docstrings of the `Analysis` class in detail for a full description of how the analysis works. +Read the comments and docstrings of the ``Analysis`` class in detail for a full description of how the analysis works. .. code-block:: python @@ -623,7 +623,7 @@ examples demonstrating more complex model-fitting tasks. This includes cookbooks, which provide a concise reference guide to the **PyAutoFit** API for advanced model-fitting: -- [Model Cookbook](https://pyautofit.readthedocs.io/en/latest/cookbooks/model.html): Learn how to compose complex models using multiple Python classes, lists, dictionaries, and customize their parameterization. +- [Model Cookbook](https://pyautofit.readthedocs.io/en/latest/cookbooks/model.html): Learn how to compose complex models using multiple Python classes, lists, dictionaries, NumPy arrays and customize their parameterization. - [Analysis Cookbook](https://pyautofit.readthedocs.io/en/latest/cookbooks/search.html): Customize the analysis with model-specific output and visualization to gain deeper insights into your model fits. diff --git a/docs/requirements.txt b/docs/requirements.txt index 701716a5c..10099ec48 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -18,7 +18,7 @@ astunparse==1.6.3 threadpoolctl>=3.1.0,<=3.2.0 autoconf timeout-decorator==0.5.0 -sphinx==5.2.3 +sphinx xxhash<=3.4.1 pyvis==0.3.2 furo diff --git a/docs/science_examples/astronomy.rst b/docs/science_examples/astronomy.rst index 9cf977585..faaa2642e 100644 --- a/docs/science_examples/astronomy.rst +++ b/docs/science_examples/astronomy.rst @@ -321,7 +321,7 @@ An example project on the **autofit_workspace** shows how to use **PyAutoFit** t lensing data, using **multi-level model composition**. If you'd like to perform the fit shown in this script, checkout the -`simple examples `_ on the +`simple examples `_ on the ``autofit_workspace``. We detail how **PyAutoFit** works in the first 3 tutorials of the `HowToFit lecture series `_. diff --git a/files/citations.bib b/files/citations.bib index cee63f8ce..84e01c822 100644 --- a/files/citations.bib +++ b/files/citations.bib @@ -76,6 +76,23 @@ @article{multinest volume = {398}, year = {2009} } +@ARTICLE{nautilus, + author = {{Lange}, Johannes U.}, + title = "{NAUTILUS: boosting Bayesian importance nested sampling with deep learning}", + journal = {\mnras}, + keywords = {methods: data analysis, methods: statistical, software: data analysis, Astrophysics - Instrumentation and Methods for Astrophysics, Astrophysics - Cosmology and Nongalactic Astrophysics, Astrophysics - Earth and Planetary Astrophysics, Astrophysics - Astrophysics of Galaxies, Computer Science - Machine Learning}, + year = 2023, + month = oct, + volume = {525}, + number = {2}, + pages = {3181-3194}, + doi = {10.1093/mnras/stad2441}, +archivePrefix = {arXiv}, + eprint = {2306.16923}, + primaryClass = {astro-ph.IM}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2023MNRAS.525.3181L}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} @ARTICLE{numpy, author={S. {van der Walt} and S. C. {Colbert} and G. {Varoquaux}}, doi={10.1109/MCSE.2011.37}, diff --git a/optional_requirements.txt b/optional_requirements.txt index 2ef2fbada..cb20d0857 100644 --- a/optional_requirements.txt +++ b/optional_requirements.txt @@ -2,6 +2,6 @@ astropy>=5.0 getdist==1.4 #jax>=0.4.13 #jaxlib>=0.4.13 -nautilus-sampler==1.0.2 -ultranest==3.6.2 +nautilus-sampler==1.0.4 +ultranest==4.3.2 zeus-mcmc==2.5.4 diff --git a/paper/paper.json b/paper/paper.json index ff2e0a5c8..c4bc2faa4 100644 --- a/paper/paper.json +++ b/paper/paper.json @@ -1,5 +1,5 @@ { - "@context": "https://raw.githubusercontent.com/codemeta/codemeta/master/codemeta.jsonld", + "@context": "https://raw.githubusercontent.com/codemeta/codemeta/main/codemeta.jsonld", "@type": "Code", "author": [ { diff --git a/readthedocs.yml b/readthedocs.yml index aa4ec0851..ca3049c2a 100644 --- a/readthedocs.yml +++ b/readthedocs.yml @@ -3,7 +3,7 @@ version: 2 build: os: ubuntu-20.04 tools: - python: "3.9" + python: "3.11" python: install: diff --git a/requirements.txt b/requirements.txt index f50f3246a..9e79b89f0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,18 +1,18 @@ -anesthetic==2.8.1 -corner==2.2.1 +anesthetic==2.8.14 +corner==2.2.2 decorator>=4.2.1 dill>=0.3.1.1 -dynesty==2.1.2 +dynesty==2.1.4 typing-inspect>=0.4.0 -emcee>=3.1.3 +emcee>=3.1.6 gprof2dot==2021.2.21 matplotlib numpydoc>=1.0.0 pyprojroot==0.2.0 pyswarms==1.3.0 -h5py>=2.10.0 -SQLAlchemy==1.3.20 -scipy<=1.11.3 +h5py>=3.11.0 +SQLAlchemy==2.0.32 +scipy<=1.14.0 astunparse==1.6.3 threadpoolctl>=3.1.0,<=3.2.0 timeout-decorator==0.5.0 diff --git a/test_autofit/analysis/test_free_parameter.py b/test_autofit/analysis/test_free_parameter.py index 3a8eb790c..5acb0fbb9 100644 --- a/test_autofit/analysis/test_free_parameter.py +++ b/test_autofit/analysis/test_free_parameter.py @@ -2,7 +2,7 @@ import autofit as af from autofit.non_linear.analysis import FreeParameterAnalysis -from autofit.non_linear.mock.mock_search import MockOptimizer +from autofit.non_linear.mock.mock_search import MockMLE def test_copy(): @@ -77,8 +77,8 @@ def make_result( combined_analysis, model, ): - optimizer = MockOptimizer() - return optimizer.fit(model, combined_analysis) + search = MockMLE() + return search.fit(model, combined_analysis) @pytest.fixture(autouse=True) diff --git a/test_autofit/config/non_linear/README.rst b/test_autofit/config/non_linear/README.rst index 11774c406..9432dd852 100644 --- a/test_autofit/config/non_linear/README.rst +++ b/test_autofit/config/non_linear/README.rst @@ -6,4 +6,4 @@ Files - ``mcmc.yaml``: Settings default behaviour of MCMC non-linear searches (e.g. Emcee). - ``nest.yaml``: Settings default behaviour of nested sampler non-linear searches (e.g. Dynesty). -- ``optimizer.yaml``: Settings default behaviour of optimizer non-linear searches (e.g. PySwarms). \ No newline at end of file +- ``mle.yaml``: Settings default behaviour of maximum likelihood estimator (mle) searches (e.g. PySwarms). \ No newline at end of file diff --git a/test_autofit/config/non_linear/optimize.yaml b/test_autofit/config/non_linear/mle.yaml similarity index 78% rename from test_autofit/config/non_linear/optimize.yaml rename to test_autofit/config/non_linear/mle.yaml index ccb50e959..b6587c4c4 100644 --- a/test_autofit/config/non_linear/optimize.yaml +++ b/test_autofit/config/non_linear/mle.yaml @@ -1,22 +1,3 @@ -DownhillSimplex: - initialize: - method: prior - printing: - silence: false - search: - disp: 1 - ftol: 0.0001 - full_output: 0 - maxfun: null - maxiter: null - retall: 0 - xtol: 0.0001 - updates: - iterations_per_update: 11 - log_every_update: 1 - model_results_every_update: 1 - remove_state_files_at_end: true - visualize_every_update: 1 Drawer: initialize: ball_lower_limit: 0.49 diff --git a/test_autofit/config/non_linear/mock.yaml b/test_autofit/config/non_linear/mock.yaml index 87ba40e65..7b78ee760 100644 --- a/test_autofit/config/non_linear/mock.yaml +++ b/test_autofit/config/non_linear/mock.yaml @@ -1,4 +1,4 @@ -MockOptimizer: +MockMLE: initialize: method: prior printing: diff --git a/test_autofit/config/priors/mock_model.yaml b/test_autofit/config/priors/mock_model.yaml index 74cbc268a..3c21fcfbf 100644 --- a/test_autofit/config/priors/mock_model.yaml +++ b/test_autofit/config/priors/mock_model.yaml @@ -1,3 +1,14 @@ +Parameter: + value: + gaussian_limits: + lower: 0.0 + upper: 1.0 + lower_limit: 0.0 + type: Uniform + upper_limit: 1.0 + width_modifier: + type: Absolute + value: 0.2 MockChildTuple: tup_0: gaussian_limits: diff --git a/test_autofit/config/priors/model.yaml b/test_autofit/config/priors/model.yaml index 0fda20896..16e704f99 100644 --- a/test_autofit/config/priors/model.yaml +++ b/test_autofit/config/priors/model.yaml @@ -66,8 +66,8 @@ PhysicalNFW: type: Absolute value: 0.2 gaussian_limits: - lower: '-1.0' - upper: '1.0' + lower: -1.0 + upper: 1.0 log10m: type: Uniform lower_limit: 6.0 diff --git a/test_autofit/database/mass_sie__source_sersic/phase_mass[sie]_source[bulge]/settings__imaging[grid_sub_2]__lens[pos_off]/dynesty_static[nlive_50__bound_multi_vol_dec_0.5_vol_check_2.0__enlarge_1.0__sample_rwalk_walks_5_facc_0.2]/pickles/settings.pickle b/test_autofit/database/mass_sie__source_sersic/phase_mass[sie]_source[bulge]/settings__imaging[grid_sub_2]__lens[pos_off]/dynesty_static[nlive_50__bound_multi_vol_dec_0.5_vol_check_2.0__enlarge_1.0__sample_rwalk_walks_5_facc_0.2]/pickles/settings.pickle index 370367a0d..b24dd9470 100644 Binary files a/test_autofit/database/mass_sie__source_sersic/phase_mass[sie]_source[bulge]/settings__imaging[grid_sub_2]__lens[pos_off]/dynesty_static[nlive_50__bound_multi_vol_dec_0.5_vol_check_2.0__enlarge_1.0__sample_rwalk_walks_5_facc_0.2]/pickles/settings.pickle and b/test_autofit/database/mass_sie__source_sersic/phase_mass[sie]_source[bulge]/settings__imaging[grid_sub_2]__lens[pos_off]/dynesty_static[nlive_50__bound_multi_vol_dec_0.5_vol_check_2.0__enlarge_1.0__sample_rwalk_walks_5_facc_0.2]/pickles/settings.pickle differ diff --git a/test_autofit/database/migration/test_integration.py b/test_autofit/database/migration/test_integration.py index b1c84d733..5336c6bc4 100644 --- a/test_autofit/database/migration/test_integration.py +++ b/test_autofit/database/migration/test_integration.py @@ -1,74 +1,35 @@ import pytest +from sqlalchemy import text from autofit.database.migration import SessionWrapper -@pytest.fixture( - autouse=True -) -def create_table( - session -): - session.execute( - "CREATE TABLE test (id INTEGER PRIMARY KEY)" - ) - - -def test_run_migration( - migrator, - session, - revision_2 -): - migrator.migrate( - session - ) - assert len(list( - session.execute( - "SELECT * FROM test" - ) - )) == 2 - - assert SessionWrapper( - session - ).revision_id == revision_2 - - -def test_apply_twice( - migrator, - session -): +@pytest.fixture(autouse=True) +def create_table(session): + session.execute(text("CREATE TABLE test (id INTEGER PRIMARY KEY)")) + + +def test_run_migration(migrator, session, revision_2): + migrator.migrate(session) + assert len(list(session.execute(text("SELECT * FROM test")))) == 2 + + assert SessionWrapper(session).revision_id == revision_2 + + +def test_apply_twice(migrator, session): for _ in range(2): - migrator.migrate( - session - ) - assert len(list( - session.execute( - "SELECT * FROM test" - ) - )) == 2 - - -def test_run_partial_migration( - migrator, - session, - revision_1, - revision_2 -): - wrapper = SessionWrapper( - session - ) + migrator.migrate(session) + assert len(list(session.execute(text("SELECT * FROM test")))) == 2 + + +def test_run_partial_migration(migrator, session, revision_1, revision_2): + wrapper = SessionWrapper(session) wrapper.revision_id = revision_1.id assert len(migrator.get_steps(revision_1.id)) == 1 - migrator.migrate( - session - ) - assert len(list( - session.execute( - "SELECT * FROM test" - ) - )) == 1 + migrator.migrate(session) + assert len(list(session.execute(text("SELECT * FROM test")))) == 1 assert wrapper.revision_id == revision_2 diff --git a/test_autofit/database/test_regression.py b/test_autofit/database/test_regression.py index bae58d94e..119731d8a 100644 --- a/test_autofit/database/test_regression.py +++ b/test_autofit/database/test_regression.py @@ -1,7 +1,8 @@ import pytest +import numpy as np import autofit as af -from autofit import database as db +from autofit import database as db, Fit @pytest.fixture(name="model") @@ -39,3 +40,32 @@ def test_model_with_parameterless_component(): def test_instance_in_collection(): collection = af.Collection(gaussian=af.Gaussian()) assert list(collection.items()) == [("gaussian", af.Gaussian())] + + +def test_samples_summary_model(): + fit = af.db.Fit() + model = af.Model(af.Gaussian) + fit["samples_summary"] = af.Samples(model=model, sample_list=[]) + fit.model = model + + assert fit["samples_summary"].model.cls == af.Gaussian + + +def test_dict_with_tuple_keys(): + d = {("a", "b"): 1} + assert db.Object.from_object(d)() == d + + +def test_persist_values(session): + fit = Fit(id=1) + + fit.set_pickle("pickle", "test") + fit.set_array("array", np.array([1, 2, 3])) + + session.add(fit) + session.commit() + + fit = session.query(Fit).first() + + assert fit["pickle"] == "test" + assert fit["array"].tolist() == [1, 2, 3] diff --git a/test_autofit/database/test_search.py b/test_autofit/database/test_search.py index 5f14e3dbb..d1873f9d0 100644 --- a/test_autofit/database/test_search.py +++ b/test_autofit/database/test_search.py @@ -4,17 +4,17 @@ def test_is_database_paths(session): - optimizer = af.m.MockOptimizer(session=session) - assert isinstance(optimizer.paths, af.DatabasePaths) + mle = af.m.MockMLE(session=session) + assert isinstance(mle.paths, af.DatabasePaths) # noinspection PyUnresolvedReferences - assert optimizer.paths.save_all_samples is False + assert mle.paths.save_all_samples is False @pytest.mark.parametrize("save_all_samples", [True, False]) def test_save_all_samples_boolean(session, save_all_samples): - optimizer = af.m.MockOptimizer(session=session, save_all_samples=save_all_samples) + mle = af.m.MockMLE(session=session, save_all_samples=save_all_samples) # noinspection PyUnresolvedReferences - assert optimizer.paths.save_all_samples is save_all_samples + assert mle.paths.save_all_samples is save_all_samples def test_unique_tag(session): @@ -23,11 +23,11 @@ def test_unique_tag(session): unique_tag = "unique" - optimizer = af.m.MockOptimizer(session=session, unique_tag=unique_tag) + mle = af.m.MockMLE(session=session, unique_tag=unique_tag) - assert optimizer.paths.unique_tag == unique_tag + assert mle.paths.unique_tag == unique_tag - optimizer.fit(model, analysis) + mle.fit(model, analysis) fit = session.query(af.db.Fit).one() diff --git a/test_autofit/graphical/functionality/test_projection.py b/test_autofit/graphical/functionality/test_projection.py index bdddea557..cd343923f 100644 --- a/test_autofit/graphical/functionality/test_projection.py +++ b/test_autofit/graphical/functionality/test_projection.py @@ -12,16 +12,16 @@ def make_q_cavity(): def test_integration(q_cavity, probit_factor): - x = np.linspace(-3, 3, 2 ** 10) + x = np.linspace(-3, 3, 2**10) probit = stats.norm(loc=0.0, scale=1.0).cdf(x) q = stats.norm(loc=q_cavity.mean, scale=q_cavity.sigma).pdf(x) tilted_distribution = probit * q - assert tilted_distribution.shape == (2 ** 10,) + assert tilted_distribution.shape == (2**10,) ni_0, ni_1, ni_2 = ( - integrate.trapz(x ** i * tilted_distribution, x) for i in range(3) + integrate.simpson(x**i * tilted_distribution, x=x) for i in range(3) ) q_numerical = autofit.messages.normal.NormalMessage.from_sufficient_statistics( diff --git a/test_autofit/graphical/gaussian/test_optimizer.py b/test_autofit/graphical/gaussian/test_optimizer.py index ea78e3d02..b7214c014 100644 --- a/test_autofit/graphical/gaussian/test_optimizer.py +++ b/test_autofit/graphical/gaussian/test_optimizer.py @@ -56,8 +56,8 @@ def test_optimisation(self, factor_model, laplace, dynesty): @pytest.mark.filterwarnings('ignore::RuntimeWarning') def test_null_paths(self, factor_model): - optimizer = af.DynestyStatic(maxcall=10) - result, status = optimizer.optimise( + search = af.DynestyStatic(maxcall=10) + result, status = search.optimise( factor_model.mean_field_approximation().factor_approximation(factor_model) ) diff --git a/test_autofit/graphical/hierarchical/test_optimise.py b/test_autofit/graphical/hierarchical/test_optimise.py index 06370f831..f500029ae 100644 --- a/test_autofit/graphical/hierarchical/test_optimise.py +++ b/test_autofit/graphical/hierarchical/test_optimise.py @@ -9,9 +9,9 @@ def make_factor(hierarchical_factor): def test_optimise(factor): - optimizer = af.DynestyStatic(maxcall=100, dynamic_delta=False, delta=0.1,) + search = af.DynestyStatic(maxcall=100, dynamic_delta=False, delta=0.1,) - _, status = optimizer.optimise( + _, status = search.optimise( factor.mean_field_approximation().factor_approximation(factor) ) assert status diff --git a/test_autofit/interpolator/test_covariance.py b/test_autofit/interpolator/test_covariance.py index a6d264842..0ab3e22ae 100644 --- a/test_autofit/interpolator/test_covariance.py +++ b/test_autofit/interpolator/test_covariance.py @@ -44,23 +44,23 @@ def test_covariance_matrix(interpolator): # ) -def n_effective(func): +def maxcall(func): return with_config( "non_linear", "nest", "DynestyStatic", "run", - "n_effective", - value=0, + "maxcall", + value=1, )(func) -@n_effective +@maxcall def test_interpolate(interpolator): assert isinstance(interpolator[interpolator.t == 0.5].gaussian.centre, float) -@n_effective +@maxcall def test_interpolate_other_field(interpolator): assert isinstance( interpolator[interpolator.gaussian.centre == 0.5].gaussian.centre, @@ -79,7 +79,7 @@ def test_model(interpolator): assert model.prior_count == 6 -@n_effective +@maxcall def test_single_variable(): samples_list = [ af.SamplesPDF( @@ -106,7 +106,7 @@ def test_single_variable(): assert interpolator[interpolator.t == 50.0].v == pytest.approx(50.0, abs=1.0) -@n_effective +@maxcall def test_variable_and_constant(): samples_list = [ af.SamplesPDF( diff --git a/test_autofit/mapper/functionality/__init__.py b/test_autofit/mapper/functionality/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/test_autofit/mapper/test_by_path.py b/test_autofit/mapper/functionality/test_by_path.py similarity index 100% rename from test_autofit/mapper/test_by_path.py rename to test_autofit/mapper/functionality/test_by_path.py diff --git a/test_autofit/mapper/test_explicit_width_modifier.py b/test_autofit/mapper/functionality/test_explicit_width_modifier.py similarity index 100% rename from test_autofit/mapper/test_explicit_width_modifier.py rename to test_autofit/mapper/functionality/test_explicit_width_modifier.py diff --git a/test_autofit/mapper/test_from_data_names.py b/test_autofit/mapper/functionality/test_from_data_names.py similarity index 100% rename from test_autofit/mapper/test_from_data_names.py rename to test_autofit/mapper/functionality/test_from_data_names.py diff --git a/test_autofit/mapper/test_has.py b/test_autofit/mapper/functionality/test_has.py similarity index 100% rename from test_autofit/mapper/test_has.py rename to test_autofit/mapper/functionality/test_has.py diff --git a/test_autofit/mapper/test_take_attributes.py b/test_autofit/mapper/functionality/test_take_attributes.py similarity index 100% rename from test_autofit/mapper/test_take_attributes.py rename to test_autofit/mapper/functionality/test_take_attributes.py diff --git a/test_autofit/mapper/model/test_regression.py b/test_autofit/mapper/model/test_regression.py index 1955cc30b..fa381fb6a 100644 --- a/test_autofit/mapper/model/test_regression.py +++ b/test_autofit/mapper/model/test_regression.py @@ -5,6 +5,7 @@ import autofit as af from autoconf.exc import ConfigException from autofit.example.model import PhysicalNFW +from autofit.mapper.mock.mock_model import WithString from autofit.mapper.model_object import Identifier @@ -158,7 +159,6 @@ def test_random_instance(model_with_assertion): class TestModel: - __test__ = False def __init__(self, items: List[float]): @@ -173,3 +173,8 @@ def test_typing_annotations(): def test_no_default_tuple_priors(): model = af.Model(PhysicalNFW) assert model.prior_count == 6 + + +def test_string_annotation(): + model = af.Model(WithString) + assert model.instance_from_prior_medians().arg.value == 0.5 diff --git a/test_autofit/mapper/test_array.py b/test_autofit/mapper/test_array.py new file mode 100644 index 000000000..836eb91f9 --- /dev/null +++ b/test_autofit/mapper/test_array.py @@ -0,0 +1,201 @@ +import pytest +import numpy as np + +import autofit as af +from autoconf.dictable import to_dict + + +@pytest.fixture +def array(): + return af.Array( + shape=(2, 2), + prior=af.GaussianPrior(mean=0.0, sigma=1.0), + ) + + +@pytest.fixture +def array_3d(): + return af.Array( + shape=(2, 2, 2), + prior=af.GaussianPrior(mean=0.0, sigma=1.0), + ) + + +def test_prior_count(array): + assert array.prior_count == 4 + + +def test_prior_count_3d(array_3d): + assert array_3d.prior_count == 8 + + +def test_instance(array): + instance = array.instance_from_prior_medians() + assert (instance == [[0.0, 0.0], [0.0, 0.0]]).all() + + +def test_instance_3d(array_3d): + instance = array_3d.instance_from_prior_medians() + assert ( + instance + == [ + [[0.0, 0.0], [0.0, 0.0]], + [[0.0, 0.0], [0.0, 0.0]], + ] + ).all() + + +def test_modify_prior(array): + array[0, 0] = 1.0 + assert array.prior_count == 3 + assert ( + array.instance_from_prior_medians() + == [ + [1.0, 0.0], + [0.0, 0.0], + ] + ).all() + + +def test_correlation(array): + array[0, 0] = array[1, 1] + array[0, 1] = array[1, 0] + + instance = array.random_instance() + + assert instance[0, 0] == instance[1, 1] + assert instance[0, 1] == instance[1, 0] + + +@pytest.fixture +def array_dict(): + return { + "arguments": { + "indices": { + "type": "list", + "values": [ + {"type": "tuple", "values": [0, 0]}, + {"type": "tuple", "values": [0, 1]}, + {"type": "tuple", "values": [1, 0]}, + {"type": "tuple", "values": [1, 1]}, + ], + }, + "prior_0_0": { + "lower_limit": float("-inf"), + "mean": 0.0, + "sigma": 1.0, + "type": "Gaussian", + "upper_limit": float("inf"), + }, + "prior_0_1": { + "lower_limit": float("-inf"), + "mean": 0.0, + "sigma": 1.0, + "type": "Gaussian", + "upper_limit": float("inf"), + }, + "prior_1_0": { + "lower_limit": float("-inf"), + "mean": 0.0, + "sigma": 1.0, + "type": "Gaussian", + "upper_limit": float("inf"), + }, + "prior_1_1": { + "lower_limit": float("-inf"), + "mean": 0.0, + "sigma": 1.0, + "type": "Gaussian", + "upper_limit": float("inf"), + }, + "shape": {"type": "tuple", "values": [2, 2]}, + }, + "type": "array", + } + + +def test_to_dict(array, array_dict, remove_ids): + assert remove_ids(to_dict(array)) == array_dict + + +def test_from_dict(array_dict): + array = af.AbstractPriorModel.from_dict(array_dict) + assert array.prior_count == 4 + assert ( + array.instance_from_prior_medians() + == [ + [0.0, 0.0], + [0.0, 0.0], + ] + ).all() + + +@pytest.fixture +def array_1d(): + return af.Array( + shape=(2,), + prior=af.GaussianPrior(mean=0.0, sigma=1.0), + ) + + +def test_1d_array(array_1d): + assert array_1d.prior_count == 2 + assert (array_1d.instance_from_prior_medians() == [0.0, 0.0]).all() + + +def test_1d_array_modify_prior(array_1d): + array_1d[0] = 1.0 + assert array_1d.prior_count == 1 + assert (array_1d.instance_from_prior_medians() == [1.0, 0.0]).all() + + +def test_tree_flatten(array): + children, aux = array.tree_flatten() + assert len(children) == 4 + assert aux == ((2, 2),) + + new_array = af.Array.tree_unflatten(aux, children) + assert new_array.prior_count == 4 + assert ( + new_array.instance_from_prior_medians() + == [ + [0.0, 0.0], + [0.0, 0.0], + ] + ).all() + + +class Analysis(af.Analysis): + def log_likelihood_function(self, instance): + return -float( + np.mean( + ( + np.array( + [ + [0.1, 0.2], + [0.3, 0.4], + ] + ) + - instance + ) + ** 2 + ) + ) + + +def test_optimisation(): + array = af.Array( + shape=(2, 2), + prior=af.UniformPrior( + lower_limit=0.0, + upper_limit=1.0, + ), + ) + result = af.DynestyStatic().fit(model=array, analysis=Analysis()) + + posterior = result.model + array[0, 0] = posterior[0, 0] + array[0, 1] = posterior[0, 1] + + result = af.DynestyStatic().fit(model=array, analysis=Analysis()) + assert isinstance(result.instance, np.ndarray) diff --git a/test_autofit/non_linear/grid/conftest.py b/test_autofit/non_linear/grid/conftest.py index a7f68379c..93e05fbb8 100644 --- a/test_autofit/non_linear/grid/conftest.py +++ b/test_autofit/non_linear/grid/conftest.py @@ -26,7 +26,7 @@ def make_sample_name_paths(): @pytest.fixture(name="grid_search_10_result") def make_grid_search_10_result(mapper, sample_name_paths): grid_search = af.SearchGridSearch( - search=af.m.MockOptimizer( + search=af.m.MockMLE( samples_summary=MockSamplesSummary( model=mapper, median_pdf_sample=Sample( diff --git a/test_autofit/non_linear/grid/test_optimizer_grid_search.py b/test_autofit/non_linear/grid/test_optimizer_grid_search.py index c407c6422..129177d2a 100644 --- a/test_autofit/non_linear/grid/test_optimizer_grid_search.py +++ b/test_autofit/non_linear/grid/test_optimizer_grid_search.py @@ -142,14 +142,14 @@ def test_raises_exception_for_bad_limits(self, grid_search, mapper): @pytest.fixture(name="grid_search_05") def make_grid_search_05(): - search = af.SearchGridSearch(search=af.m.MockOptimizer(), number_of_steps=2) + search = af.SearchGridSearch(search=af.m.MockMLE(), number_of_steps=2) search.search.paths = af.DirectoryPaths(name="sample_name") return search @pytest.fixture(autouse=True) def empty_args(): - af.m.MockOptimizer.init_args = list() + af.m.MockMLE.init_args = list() def test_csv_headers(grid_search_10_result, sample_name_paths): diff --git a/test_autofit/non_linear/grid/test_paths/test_database_run.py b/test_autofit/non_linear/grid/test_paths/test_database_run.py index b98434f4d..602b5d7e3 100644 --- a/test_autofit/non_linear/grid/test_paths/test_database_run.py +++ b/test_autofit/non_linear/grid/test_paths/test_database_run.py @@ -8,7 +8,7 @@ ) def make_search(session): return af.SearchGridSearch( - search=af.m.MockOptimizer( + search=af.m.MockMLE( session=session ), number_of_steps=2 diff --git a/test_autofit/non_linear/grid/test_paths/test_indicators.py b/test_autofit/non_linear/grid/test_paths/test_indicators.py index a8d5f112b..c9b9696c8 100644 --- a/test_autofit/non_linear/grid/test_paths/test_indicators.py +++ b/test_autofit/non_linear/grid/test_paths/test_indicators.py @@ -24,7 +24,7 @@ def make_database_parent_search(session): def _make_grid_search(mapper, session=None): search = af.SearchGridSearch( - search=af.m.MockOptimizer(session=session), number_of_steps=2 + search=af.m.MockMLE(session=session), number_of_steps=2 ) search.fit( model=mapper, diff --git a/test_autofit/non_linear/grid/test_quantities_for_path.py b/test_autofit/non_linear/grid/test_quantities_for_path.py new file mode 100644 index 000000000..539b35b10 --- /dev/null +++ b/test_autofit/non_linear/grid/test_quantities_for_path.py @@ -0,0 +1,76 @@ +import pytest + +import autofit as af + + +@pytest.fixture(name="result") +def make_result(): + model = af.Model( + af.Gaussian, + ) + grid_priors = [model.centre, model.normalization] + lower_limits_lists = [ + [0.0, 1.0], + [0.0, 2.0], + [2.0, 1.0], + [2.0, 2.0], + ] + + sample = af.Sample( + 1.0, + 1.0, + 1.0, + { + "centre": 1.0, + "normalization": 2.0, + "sigma": 3.0, + }, + ) + + def make_samples(centre, normalization): + return af.Samples( + model=af.Model( + af.Gaussian, + centre=af.UniformPrior( + lower_limit=centre, + upper_limit=centre + 2.0, + ), + normalization=af.UniformPrior( + lower_limit=normalization, + upper_limit=normalization + 1.0, + ), + ), + sample_list=[sample], + ) + + samples = [ + make_samples(centre, normalisation) + for centre, normalisation in lower_limits_lists + ] + + return af.GridSearchResult( + samples=samples, + lower_limits_lists=lower_limits_lists, + grid_priors=grid_priors, + ) + + +@pytest.mark.parametrize( + "name, expected", + [ + ("centre", [1.0, 1.0, 3.0, 3.0]), + ("normalization", [1.5, 2.5, 1.5, 2.5]), + ], +) +def test_physical_centres_from(result, name, expected): + assert result.physical_centres_lists_from(name) == expected + assert result.shape == (2, 2) + + +def test_two_physical_centres(result): + assert result.physical_centres_lists_from(("centre", "normalization")) == [ + (1.0, 1.5), + (1.0, 2.5), + (3.0, 1.5), + (3.0, 2.5), + ] diff --git a/test_autofit/non_linear/grid/test_sensitivity/conftest.py b/test_autofit/non_linear/grid/test_sensitivity/conftest.py index e807b2a16..e2a71d2eb 100644 --- a/test_autofit/non_linear/grid/test_sensitivity/conftest.py +++ b/test_autofit/non_linear/grid/test_sensitivity/conftest.py @@ -16,8 +16,8 @@ def __init__(self): def __call__(self, instance: af.ModelInstance, simulate_path: Optional[str]): image = instance.gaussian(x) - if hasattr(instance, "perturbation"): - image += instance.perturbation(x) + if hasattr(instance, "perturb"): + image += instance.perturb(x) return image @@ -38,10 +38,11 @@ class BaseFit: def __init__(self, analysis_cls): self.analysis_cls = analysis_cls - def __call__(self, dataset, model, paths): + def __call__(self, dataset, model, paths, instance): search = af.m.MockSearch( return_sensitivity_results=True, samples_summary=MockSamplesSummary(model=model), + paths=paths, ) analysis = self.analysis_cls(dataset=dataset) @@ -53,10 +54,11 @@ class PerturbFit: def __init__(self, analysis_cls): self.analysis_cls = analysis_cls - def __call__(self, dataset, model, paths): + def __call__(self, dataset, model, paths, instance): search = af.m.MockSearch( return_sensitivity_results=True, samples_summary=MockSamplesSummary(model=model), + paths=paths, ) analysis = self.analysis_cls(dataset=dataset) @@ -88,6 +90,37 @@ def make_sensitivity( ) +@pytest.fixture(name="masked_sensitivity") +def make_masked_sensitivity( + perturb_model, +): + # noinspection PyTypeChecker + instance = af.ModelInstance() + instance.gaussian = af.Gaussian() + return s.Sensitivity( + simulation_instance=instance, + base_model=af.Collection(gaussian=af.Model(af.Gaussian)), + perturb_model=perturb_model, + simulate_cls=Simulate(), + base_fit_cls=BaseFit(Analysis), + perturb_fit_cls=PerturbFit(Analysis), + paths=af.DirectoryPaths(), + number_of_steps=2, + mask=np.array( + [ + [ + [True, True], + [True, True], + ], + [ + [True, True], + [True, True], + ], + ] + ), + ) + + @pytest.fixture(name="job") def make_job( perturb_model, diff --git a/test_autofit/non_linear/grid/test_sensitivity/test_functionality.py b/test_autofit/non_linear/grid/test_sensitivity/test_functionality.py index 4d1330d1c..612f0d7c1 100644 --- a/test_autofit/non_linear/grid/test_sensitivity/test_functionality.py +++ b/test_autofit/non_linear/grid/test_sensitivity/test_functionality.py @@ -29,11 +29,27 @@ def test_labels(sensitivity): def test_perform_job(job): + assert not job.is_complete + + result = job.perform() + assert isinstance(result, s.JobResult) + assert isinstance(result.perturb_result, af.Result) + assert isinstance(result.result, af.Result) + + assert job.is_complete + + +def test_perform_twice(job): + job.perform() + assert job.is_complete + result = job.perform() assert isinstance(result, s.JobResult) assert isinstance(result.perturb_result, af.Result) assert isinstance(result.result, af.Result) + assert job.is_complete + class TestPerturbationModels: @pytest.mark.parametrize( @@ -122,11 +138,6 @@ def test_prior_with_limits(self): assert prior.lower_limit == 3 assert prior.upper_limit == 5 - def test_existing_limits(self): - prior = af.UniformPrior(2, 4).with_limits(3, 5) - assert prior.lower_limit == 3 - assert prior.upper_limit == 4 - @pytest.fixture(name="tuple_sensitivity") def make_tuple_sensitivity(sensitivity): diff --git a/test_autofit/non_linear/grid/test_sensitivity/test_masked_sensitivity.py b/test_autofit/non_linear/grid/test_sensitivity/test_masked_sensitivity.py new file mode 100644 index 000000000..bf9c33314 --- /dev/null +++ b/test_autofit/non_linear/grid/test_sensitivity/test_masked_sensitivity.py @@ -0,0 +1,71 @@ +from math import prod + +import pytest + +import autofit as af + + +@pytest.fixture(name="masked_result") +def make_masked_result(masked_sensitivity): + return masked_sensitivity.run() + + +def test_result_size(masked_sensitivity, masked_result): + number_elements = prod(masked_sensitivity.shape) + assert len(masked_result.samples) == number_elements + + +def test_sample(masked_result): + sample = masked_result.samples[0] + assert sample.model is not None + assert sample.model.perturb is not None + assert sample.log_evidence == 0.0 + assert sample.log_likelihood == 0.0 + + +@pytest.mark.parametrize( + "lower, upper, mean", + [ + (0.0, 1.0, 0.5), + (-1.0, 1.0, 0.0), + (-1.0, 0.0, -0.5), + (0.5, 1.0, 0.75), + ], +) +def test_mean_uniform_prior( + lower, + upper, + mean, +): + prior = af.UniformPrior( + lower_limit=0.0, + upper_limit=1.0, + ) + assert ( + prior.with_limits( + lower, + upper, + ).mean + == mean + ) + + +def test_path_value_dicts(masked_sensitivity): + assert masked_sensitivity.path_values == { + ("centre",): [0.25, 0.25, 0.25, 0.25, 0.75, 0.75, 0.75, 0.75], + ("normalization",): [0.25, 0.25, 0.75, 0.75, 0.25, 0.25, 0.75, 0.75], + ("sigma",): [0.25, 0.75, 0.25, 0.75, 0.25, 0.75, 0.25, 0.75], + } + + +def test_perturbed_physical_centres_list_from(masked_result): + assert masked_result.perturbed_physical_centres_list_from("centre") == [ + 0.25, + 0.25, + 0.25, + 0.25, + 0.75, + 0.75, + 0.75, + 0.75, + ] diff --git a/test_autofit/non_linear/grid/test_sensitivity/test_results.py b/test_autofit/non_linear/grid/test_sensitivity/test_results.py index 343adbfc8..655794eba 100644 --- a/test_autofit/non_linear/grid/test_sensitivity/test_results.py +++ b/test_autofit/non_linear/grid/test_sensitivity/test_results.py @@ -2,18 +2,25 @@ from autofit.non_linear.grid.sensitivity.job import JobResult from autofit.non_linear.grid.sensitivity.result import SensitivityResult +import autofit as af class Samples: - def __init__(self, log_likelihood): self.log_likelihood = log_likelihood + self.model = af.Model( + af.Gaussian, + centre=af.UniformPrior( + 0.0, + 1.0, + ), + ) def summary(self): return self -class Result: +class Result: def __init__(self, samples): self.samples = samples @@ -35,13 +42,24 @@ def test_job_result(job_result): assert job_result.log_likelihood_increase == 1.0 -def test_result(job_result): - result = SensitivityResult( +@pytest.fixture(name="sensitivity_result") +def make_sensitivity_result(job_result): + return SensitivityResult( samples=[job_result.result.samples.summary()], perturb_samples=[job_result.perturb_result.samples.summary()], - shape=(1,) + shape=(1,), + path_values={ + ("centre",): [0.5], + }, ) - assert result.log_likelihoods_base == [1.0] - assert result.log_likelihoods_perturbed == [2.0] - assert result.log_likelihood_differences == [1.0] + +def test_result(sensitivity_result): + assert sensitivity_result.log_likelihoods_base == [1.0] + assert sensitivity_result.log_likelihoods_perturbed == [2.0] + assert sensitivity_result.log_likelihood_differences == [1.0] + + +def test_physical_centres(sensitivity_result): + assert sensitivity_result.physical_centres_lists_from("centre") == [0.5] + assert sensitivity_result.perturbed_physical_centres_list_from("centre") == [0.5] diff --git a/test_autofit/non_linear/paths/test_paths.py b/test_autofit/non_linear/paths/test_paths.py index ffd784988..e123e897e 100644 --- a/test_autofit/non_linear/paths/test_paths.py +++ b/test_autofit/non_linear/paths/test_paths.py @@ -39,7 +39,9 @@ def test_paths_argument(self): self.assert_paths_as_expected(search.paths) def test_combination_argument(self): - search = af.m.MockSearch("other",) + search = af.m.MockSearch( + "other", + ) search.paths = af.DirectoryPaths(name="name") self.assert_paths_as_expected(search.paths) @@ -69,3 +71,8 @@ def test_serialize(model): pickled_paths = pickle.loads(pickle.dumps(paths)) assert pickled_paths.model is not None + + +def test_unique_tag(): + paths = af.DirectoryPaths(unique_tag="unique_tag") + assert "unique_tag" in paths.output_path.parts diff --git a/test_autofit/non_linear/test_initializer.py b/test_autofit/non_linear/test_initializer.py index c63397306..4d117ebe2 100644 --- a/test_autofit/non_linear/test_initializer.py +++ b/test_autofit/non_linear/test_initializer.py @@ -17,7 +17,8 @@ def __call__(self, parameters): return self.figure_of_merit -def test__priors__samples_from_model(): +@pytest.fixture +def model_and_samples(): model = af.Model(af.m.MockClassx4) model.one = af.UniformPrior(lower_limit=0.099, upper_limit=0.101) model.two = af.UniformPrior(lower_limit=0.199, upper_limit=0.201) @@ -26,34 +27,42 @@ def test__priors__samples_from_model(): initializer = af.InitializerPrior() - ( - unit_parameter_lists, - parameter_lists, - figure_of_merit_list, - ) = initializer.samples_from_model( + unit_parameter_lists, parameter_lists, _ = initializer.samples_from_model( total_points=2, model=model, fitness=MockFitness(), paths=af.DirectoryPaths(), ) - assert 0.0 < unit_parameter_lists[0][0] < 1.0 - assert 0.0 < unit_parameter_lists[1][0] < 1.0 - assert 0.0 < unit_parameter_lists[0][1] < 1.0 - assert 0.0 < unit_parameter_lists[1][1] < 1.0 - assert 0.0 < unit_parameter_lists[0][2] < 1.0 - assert 0.0 < unit_parameter_lists[1][2] < 1.0 - assert 0.0 < unit_parameter_lists[0][3] < 1.0 - assert 0.0 < unit_parameter_lists[1][3] < 1.0 - - assert 0.099 < parameter_lists[0][0] < 0.101 - assert 0.099 < parameter_lists[1][0] < 0.101 - assert 0.199 < parameter_lists[0][1] < 0.201 - assert 0.199 < parameter_lists[1][1] < 0.201 - assert 0.299 < parameter_lists[0][2] < 0.301 - assert 0.299 < parameter_lists[1][2] < 0.301 - assert 0.399 < parameter_lists[0][3] < 0.401 - assert 0.399 < parameter_lists[1][3] < 0.401 + return unit_parameter_lists, parameter_lists + +@pytest.mark.parametrize("index, param_index, lower, upper", [ + (0, 0, 0.0, 1.0), + (1, 0, 0.0, 1.0), + (0, 1, 0.0, 1.0), + (1, 1, 0.0, 1.0), + (0, 2, 0.0, 1.0), + (1, 2, 0.0, 1.0), + (0, 3, 0.0, 1.0), + (1, 3, 0.0, 1.0), +]) +def test_unit_parameter_lists(model_and_samples, index, param_index, lower, upper): + unit_parameter_lists, _ = model_and_samples + assert lower < unit_parameter_lists[index][param_index] < upper + +@pytest.mark.parametrize("index, param_index, lower, upper", [ + (0, 0, 0.099, 0.101), + (1, 0, 0.099, 0.101), + (0, 1, 0.199, 0.201), + (1, 1, 0.199, 0.201), + (0, 2, 0.299, 0.301), + (1, 2, 0.299, 0.301), + (0, 3, 0.399, 0.401), + (1, 3, 0.399, 0.401), +]) +def test_parameter_lists(model_and_samples, index, param_index, lower, upper): + _, parameter_lists = model_and_samples + assert lower < parameter_lists[index][param_index] < upper def test__priors__samples_from_model__raise_exception_if_all_likelihoods_identical(): @@ -234,8 +243,8 @@ def make_model(): ) -def test_starting_point_initializer(model): - initializer = af.SpecificRangeInitializer( +def test_initializer_bounds(model): + initializer = af.InitializerParamBounds( { model.centre: (1.0, 2.0), model.normalization: (2.0, 3.0), @@ -249,8 +258,53 @@ def test_starting_point_initializer(model): assert 0.0 <= parameter <= 1.0 +def test__initializer_bounds__info_from_model(model): + + initializer = af.InitializerParamBounds( + { + model.centre: (1.0, 2.0), + model.normalization: (2.0, 3.0), + model.sigma: (-2.0, -1.0), + } + ) + + info = initializer.info_from_model(model) + + assert "Total Free Parameters = 3" in info + assert "centre: Start[(1.0, 2.0)]" in info + + +def test__initializer_start_point(model): + initializer = af.InitializerParamStartPoints( + { + model.centre: 1.5, + model.normalization: 2.5, + model.sigma: -1.5, + } + ) + + parameter_list = initializer._generate_unit_parameter_list(model) + assert len(parameter_list) == 3 + for parameter in parameter_list: + assert 0.0 <= parameter <= 1.0 + + +def test__initializer_start_point__info_from_model(model): + initializer = af.InitializerParamStartPoints( + { + model.centre: 1.5, + model.normalization: 2.5, + model.sigma: -1.5, + } + ) + info = initializer.info_from_model(model) + + assert "Total Free Parameters = 3" in info + assert "centre: Start[1.5]" in info + + def test_offset(model): - initializer = af.SpecificRangeInitializer( + initializer = af.InitializerParamBounds( { model.centre: (1.5, 2.0), model.normalization: (2.5, 3.0), @@ -265,7 +319,7 @@ def test_offset(model): def test_missing_parameter(model): - initializer = af.SpecificRangeInitializer( + initializer = af.InitializerParamBounds( { model.centre: (1.5, 2.0), model.normalization: (2.5, 3.0), diff --git a/test_autofit/non_linear/test_regression.py b/test_autofit/non_linear/test_regression.py index 82591b64d..64876ce63 100644 --- a/test_autofit/non_linear/test_regression.py +++ b/test_autofit/non_linear/test_regression.py @@ -14,18 +14,18 @@ def test_no_priors(): search.fit(model, af.Analysis()) -@pytest.fixture(name="optimizer") -def make_optimizer(): +@pytest.fixture(name="search") +def make_search(): return af.DynestyStatic("name") -def test_serialize_optimiser(optimizer): - optimizer = pickle.loads(pickle.dumps(optimizer)) - assert optimizer.name == "name" +def test_serialize_optimiser(search): + search = pickle.loads(pickle.dumps(search)) + assert search.name == "name" -def test_serialize_grid_search(optimizer): - grid_search = af.SearchGridSearch(optimizer) +def test_serialize_grid_search(search): + grid_search = af.SearchGridSearch(search) assert grid_search.logger.name == "GridSearch (name)" assert "logger" not in grid_search.__getstate__() diff --git a/test_autofit/serialise/test_samples.py b/test_autofit/serialise/test_samples.py index d621b1cdb..438038fab 100644 --- a/test_autofit/serialise/test_samples.py +++ b/test_autofit/serialise/test_samples.py @@ -46,6 +46,7 @@ def test_summary(summary, model, sample): assert summary.max_log_likelihood_sample == sample + @pytest.fixture(name="summary_dict") def make_summary_dict(): return { @@ -69,7 +70,6 @@ def make_summary_dict(): }, }, }, - "log_evidence": None, "values_at_sigma_3": { "type": "list", "values": [ @@ -78,21 +78,26 @@ def make_summary_dict(): {"type": "tuple", "values": [2.0, 6.0]}, ], }, - "values_at_sigma_1": { - "type": "list", - "values": [ - {"type": "tuple", "values": [0.0, 2.0]}, - {"type": "tuple", "values": [1.0, 4.0]}, - {"type": "tuple", "values": [2.0, 6.0]}, - ], - }, - "errors_at_sigma_3": { - "type": "list", - "values": [ - {"type": "tuple", "values": [2.0, 0.0]}, - {"type": "tuple", "values": [3.0, 0.0]}, - {"type": "tuple", "values": [4.0, 0.0]}, - ], + "model": { + "class_path": "autofit.example.model.Gaussian", + "type": "model", + "arguments": { + "centre": { + "lower_limit": 0.0, + "upper_limit": 1.0, + "type": "Uniform", + }, + "normalization": { + "lower_limit": 0.0, + "upper_limit": 1.0, + "type": "Uniform", + }, + "sigma": { + "lower_limit": 0.0, + "upper_limit": 1.0, + "type": "Uniform", + }, + }, }, "max_log_likelihood_sample": { "type": "instance", @@ -111,6 +116,14 @@ def make_summary_dict(): }, }, }, + "values_at_sigma_1": { + "type": "list", + "values": [ + {"type": "tuple", "values": [0.0, 2.0]}, + {"type": "tuple", "values": [1.0, 4.0]}, + {"type": "tuple", "values": [2.0, 6.0]}, + ], + }, "errors_at_sigma_1": { "type": "list", "values": [ @@ -119,6 +132,15 @@ def make_summary_dict(): {"type": "tuple", "values": [4.0, 0.0]}, ], }, + "errors_at_sigma_3": { + "type": "list", + "values": [ + {"type": "tuple", "values": [2.0, 0.0]}, + {"type": "tuple", "values": [3.0, 0.0]}, + {"type": "tuple", "values": [4.0, 0.0]}, + ], + }, + "log_evidence": None, }, } diff --git a/test_autofit/test_correspondence.py b/test_autofit/test_correspondence.py index 94b44155f..b2dfb9805 100644 --- a/test_autofit/test_correspondence.py +++ b/test_autofit/test_correspondence.py @@ -109,4 +109,4 @@ def test_regression(): assert log_likelihood == numerical_log_likelihood assert float(gradient) == pytest.approx(numerical_gradient, rel=0.001) - assert float(approx_gradient) == pytest.approx(numerical_gradient, rel=0.001) + assert float(approx_gradient[0]) == pytest.approx(numerical_gradient, rel=0.001)