Merge pull request #63 from Jammy2211/feature/save_delta_ellipticity

Feature/save delta ellipticity

autocti/__init__.py

@@ -62,6 +62,7 @@
 from .dataset_1d.model.result import ResultDataset1D
 from .charge_injection.model.analysis import AnalysisImagingCI
 from .charge_injection.model.result import ResultImagingCI
+from .model.analysis import AnalysisCTI
 from .model.model_util import CTI1D
 from .model.model_util import CTI2D
 from .model.settings import SettingsCTI1D

autocti/charge_injection/model/analysis.py

@@ -6,14 +6,14 @@
 
 import autoarray as aa
 import autofit as af
-from autoconf.dictable import output_to_json
 
 from autocti.charge_injection.imaging.imaging import ImagingCI
 from autocti.charge_injection.fit import FitImagingCI
 from autocti.charge_injection.model.visualizer import VisualizerImagingCI
 from autocti.charge_injection.model.result import ResultImagingCI
 from autocti.clocker.two_d import Clocker2D
 from autocti.charge_injection.hyper import HyperCINoiseCollection
+from autocti.model.analysis import AnalysisCTI
 from autocti.model.settings import SettingsCTI2D
 from autocti.preloads import Preloads
 
@@ -24,7 +24,7 @@
 logger.setLevel(level="INFO")
 
 
-class AnalysisImagingCI(af.Analysis):
+class AnalysisImagingCI(AnalysisCTI):
     def __init__(
         self,
         dataset: ImagingCI,
@@ -56,12 +56,12 @@ def __init__(
             The full dataset, which is visualized separate from the `dataset` that is fitted, which for example may
             not have the FPR masked and thus enable visualization of the FPR.
         """
-        super().__init__()
-
-        self.dataset = dataset
-        self.clocker = clocker
-        self.settings_cti = settings_cti
-        self.dataset_full = dataset_full
+        super().__init__(
+            dataset=dataset,
+            clocker=clocker,
+            settings_cti=settings_cti,
+            dataset_full=dataset_full,
+        )
 
         self.preloads = Preloads()
 
@@ -296,16 +296,6 @@ def output_dataset(dataset, prefix):
         if self.dataset_full is not None:
             output_dataset(dataset=self.dataset_full, prefix="dataset_full")
 
-    def in_ascending_fpr_order_from(self, quantity_list, fpr_value_list):
-        if not conf.instance["visualize"]["general"]["general"][
-            "subplot_ascending_fpr"
-        ]:
-            return quantity_list
-
-        indexes = sorted(range(len(fpr_value_list)), key=lambda k: fpr_value_list[k])
-
-        return [quantity_list[i] for i in indexes]
-
     def visualize_before_fit(self, paths: af.DirectoryPaths, model: af.Collection):
         if conf.instance["visualize"]["plots"]["combined_only"]:
             return

autocti/dataset_1d/model/analysis.py

@@ -1,20 +1,19 @@
 from typing import List, Optional
 
-from autoconf import conf
 from autoconf.dictable import to_dict
 
 import autofit as af
-from autoconf.dictable import output_to_json
 
 from autocti.dataset_1d.dataset_1d.dataset_1d import Dataset1D
 from autocti.dataset_1d.fit import FitDataset1D
 from autocti.dataset_1d.model.visualizer import VisualizerDataset1D
 from autocti.dataset_1d.model.result import ResultDataset1D
+from autocti.model.analysis import AnalysisCTI
 from autocti.model.settings import SettingsCTI1D
 from autocti.clocker.one_d import Clocker1D
 
 
-class AnalysisDataset1D(af.Analysis):
+class AnalysisDataset1D(AnalysisCTI):
     def __init__(
         self,
         dataset: Dataset1D,
@@ -46,12 +45,12 @@ def __init__(
             The full dataset, which is visualized separate from the `dataset` that is fitted, which for example may
             not have the FPR masked and thus enable visualization of the FPR.
         """
-        super().__init__()
-
-        self.dataset = dataset
-        self.clocker = clocker
-        self.settings_cti = settings_cti
-        self.dataset_full = dataset_full
+        super().__init__(
+            dataset=dataset,
+            clocker=clocker,
+            settings_cti=settings_cti,
+            dataset_full=dataset_full,
+        )
 
     def region_list_from(self) -> List:
         return ["fpr", "eper"]
@@ -182,16 +181,6 @@ def output_dataset(dataset, prefix):
             object_dict=to_dict(self.settings_cti),
         )
 
-    def in_ascending_fpr_order_from(self, quantity_list, fpr_value_list):
-        if not conf.instance["visualize"]["general"]["general"][
-            "subplot_ascending_fpr"
-        ]:
-            return quantity_list
-
-        indexes = sorted(range(len(fpr_value_list)), key=lambda k: fpr_value_list[k])
-
-        return [quantity_list[i] for i in indexes]
-
     def visualize_before_fit(self, paths: af.DirectoryPaths, model: af.Collection):
         region_list = self.region_list_from()
 

autocti/exc.py

@@ -1,4 +1,4 @@
-from autofit import exc
+from autofit.exc import *
 from autoarray.exc import ArrayException, MaskException, RegionException
 
 
@@ -30,7 +30,7 @@ class FittingException(Exception):
     pass
 
 
-class PriorException(exc.FitException):
+class PriorException(FitException):
     pass
 
 

autocti/model/analysis.py

@@ -0,0 +1,117 @@
+from typing import List, Optional, Union
+
+from autoconf import conf
+from autoconf.dictable import output_to_json
+
+import autoarray as aa
+import autofit as af
+
+from autocti.charge_injection.imaging.imaging import ImagingCI
+from autocti.charge_injection.model.result import ResultImagingCI
+from autocti.clocker.one_d import Clocker1D
+from autocti.clocker.two_d import Clocker2D
+from autocti.dataset_1d.dataset_1d.dataset_1d import Dataset1D
+from autocti.model.settings import SettingsCTI1D
+from autocti.model.settings import SettingsCTI2D
+
+from autocti import exc
+
+
+class AnalysisCTI(af.Analysis):
+    def __init__(
+        self,
+        dataset: Union[Dataset1D, ImagingCI],
+        clocker: Union[Clocker1D, Clocker2D],
+        settings_cti: Union[SettingsCTI1D, SettingsCTI2D],
+        dataset_full: Optional[aa.AbstractDataset] = None,
+    ):
+        """
+        Fits a CTI model to a CTI calibration dataset via a non-linear search.
+
+        The `Analysis` class defines the `log_likelihood_function` which fits the model to the dataset and returns the
+        log likelihood value defining how well the model fitted the data.
+
+        It handles many other tasks, such as visualization, outputting results to hard-disk and storing results in
+        a format that can be loaded after the model-fit is complete.
+
+        This class is used for model-fits which fit a CTI model via a `CTI2D` object to a charge injection
+        imaging dataset.
+
+        Parameters
+        ----------
+        dataset
+            The charge injection dataset that the model is fitted to.
+        clocker
+            The CTI arctic clocker used by the non-linear search and model-fit.
+        settings_cti
+            The settings controlling aspects of the CTI model in this model-fit.
+        dataset_full
+            The full dataset, which is visualized separate from the `dataset` that is fitted, which for example may
+            not have the FPR masked and thus enable visualization of the FPR.
+        """
+        super().__init__()
+
+        self.dataset = dataset
+        self.clocker = clocker
+        self.settings_cti = settings_cti
+        self.dataset_full = dataset_full
+
+    def region_list_from(self) -> List:
+        raise NotImplementedError
+
+    def save_results(self, paths: af.DirectoryPaths, result: ResultImagingCI):
+        """
+        At the end of a model-fit, this routine saves attributes of the `Analysis` object to the `files`
+        folder such that they can be loaded after the analysis using PyAutoFit's database and aggregator tools.
+
+        For this analysis it outputs the following:
+
+        - The Israel et al requirement on the spurious ellipticity based on the errors of the fit.
+
+        Parameters
+        ----------
+        paths
+            The PyAutoFit 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 and samples.
+        """
+
+        weight_list = []
+        delta_ellipticity_list = []
+
+        for sample in result.samples.sample_list:
+            try:
+                instance = sample.instance_for_model(model=result.samples.model)
+            except exc.FitException:
+                continue
+
+            weight_list.append(sample.weight)
+            delta_ellipticity_list.append(instance.cti.delta_ellipticity)
+
+        (
+            median_delta_ellipticity,
+            upper_delta_ellipticity,
+            lower_delta_ellipticity,
+        ) = af.marginalize(
+            parameter_list=delta_ellipticity_list,
+            sigma=2.0,
+            weight_list=weight_list,
+        )
+
+        delta_ellipticity = (upper_delta_ellipticity - lower_delta_ellipticity) / 2.0
+
+        output_to_json(
+            obj=delta_ellipticity,
+            file_path=paths._files_path / "delta_ellipticity.json",
+        )
+
+    def in_ascending_fpr_order_from(self, quantity_list, fpr_value_list):
+        if not conf.instance["visualize"]["general"]["general"][
+            "subplot_ascending_fpr"
+        ]:
+            return quantity_list
+
+        indexes = sorted(range(len(fpr_value_list)), key=lambda k: fpr_value_list[k])
+
+        return [quantity_list[i] for i in indexes]

test_autocti/charge_injection/model/test_analysis_ci.py

@@ -1,5 +1,7 @@
 import copy
+import os
 import pytest
+
 import autofit as af
 import autocti as ac
 

test_autocti/model/test_analysis.py

@@ -0,0 +1,52 @@
+import os
+import pytest
+
+import autofit as af
+import autocti as ac
+
+
+def test__save_results__delta_ellipyicity_output_to_json(
+    imaging_ci_7x7,
+    pre_cti_data_7x7,
+    traps_x1,
+    ccd,
+    parallel_clocker_2d,
+):
+    analysis = ac.AnalysisCTI(
+        dataset=imaging_ci_7x7,
+        clocker=parallel_clocker_2d,
+        settings_cti=ac.SettingsCTI2D(),
+    )
+
+    paths = af.DirectoryPaths()
+
+    model = af.Collection(
+        cti=af.Model(
+            ac.CTI2D, parallel_trap_list=[ac.TrapInstantCapture], parallel_ccd=ccd
+        ),
+    )
+
+    parameters = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0]]
+
+    sample_list = af.Sample.from_lists(
+        model=model,
+        parameter_lists=parameters,
+        log_likelihood_list=[1.0, 2.0, 3.0],
+        log_prior_list=[0.0, 0.0, 0.0],
+        weight_list=[0.2, 0.2, 0.2],
+    )
+
+    samples = ac.m.MockSamples(sample_list=sample_list, model=model)
+
+    analysis.save_results(
+        paths=paths,
+        result=ac.m.MockResult(samples=samples, model=model),
+    )
+
+    delta_ellipticity = ac.from_json(
+        file_path=paths._files_path / "delta_ellipticity.json"
+    )
+
+    assert delta_ellipticity == pytest.approx(-1.6145994035538491, 1.0e-4)
+
+    os.remove(paths._files_path / "delta_ellipticity.json")