notebooks/cookbooks/analysis.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Cookbook: Analysis\n",
        "==================\n",
        "\n",
        "The `Analysis` class is the interface between the data and model, whereby a `log_likelihood_function` is defined\n",
        "and called by the non-linear search to fit the model.\n",
        "\n",
        "This cookbook provides an overview of how to use and extend `Analysis` objects in **PyAutoFit**.\n",
        "\n",
        "__Contents__\n",
        "\n",
        " - Example: A simple example of an analysis class which can be adapted for you use-case.\n",
        " - Customization: Customizing an analysis class with different data inputs and editing the `log_likelihood_function`.\n",
        " - Visualization: Using a `visualize` method so that model-specific visuals are output to hard-disk.\n",
        " - Custom Result: Return a custom Result object with methods specific to your model fitting problem.\n",
        " - Latent Variables: Adding a `compute_latent_variables` method to the analysis to output latent variables to hard-disk.\n",
        " - Custom Output: Add methods which output model-specific results to hard-disk in the `files` folder (e.g. as .json\n",
        "   files) to aid in the interpretation of results."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "%matplotlib inline\n",
        "from pyprojroot import here\n",
        "workspace_path = str(here())\n",
        "%cd $workspace_path\n",
        "print(f\"Working Directory has been set to `{workspace_path}`\")\n",
        "\n",
        "import json\n",
        "import numpy as np\n",
        "from os import path\n",
        "from typing import Dict, Optional\n",
        "\n",
        "import autofit as af\n",
        "import autofit.plot as aplt"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Example__\n",
        "\n",
        "An example simple `Analysis` class, to remind ourselves of the basic structure and inputs.\n",
        "\n",
        "This can be adapted for your use case."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "\n",
        "\n",
        "class Analysis(af.Analysis):\n",
        "    def __init__(self, data: np.ndarray, noise_map: np.ndarray):\n",
        "        \"\"\"\n",
        "        The `Analysis` class acts as an interface between the data and model in **PyAutoFit**.\n",
        "\n",
        "        Its `log_likelihood_function` defines how the model is fitted to the data and it is called many times by\n",
        "        the non-linear search fitting algorithm.\n",
        "\n",
        "        In this example the `Analysis` `__init__` constructor only contains the `data` and `noise-map`, but it can be\n",
        "        easily extended to include other quantities.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        data\n",
        "            A 1D numpy array containing the data (e.g. a noisy 1D signal) fitted in the workspace examples.\n",
        "        noise_map\n",
        "            A 1D numpy array containing the noise values of the data, used for computing the goodness of fit\n",
        "            metric, the log likelihood.\n",
        "        \"\"\"\n",
        "        super().__init__()\n",
        "\n",
        "        self.data = data\n",
        "        self.noise_map = noise_map\n",
        "\n",
        "    def log_likelihood_function(self, instance) -> float:\n",
        "        \"\"\"\n",
        "        Returns the log likelihood of a fit of a 1D Gaussian to the dataset.\n",
        "\n",
        "        The data is fitted using an `instance` of the `Gaussian` class where its `model_data_from`\n",
        "        is called in order to create a model data representation of the Gaussian that is fitted to the data.\n",
        "        \"\"\"\n",
        "\n",
        "        xvalues = np.arange(self.data.shape[0])\n",
        "\n",
        "        model_data = instance.model_data_from(xvalues=xvalues)\n",
        "\n",
        "        residual_map = self.data - model_data\n",
        "        chi_squared_map = (residual_map / self.noise_map) ** 2.0\n",
        "        chi_squared = sum(chi_squared_map)\n",
        "        noise_normalization = np.sum(np.log(2 * np.pi * self.noise_map**2.0))\n",
        "        log_likelihood = -0.5 * (chi_squared + noise_normalization)\n",
        "\n",
        "        return log_likelihood\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "An instance of the analysis class is created as follows."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dataset_path = path.join(\"dataset\", \"example_1d\", \"gaussian_x1\")\n",
        "data = af.util.numpy_array_from_json(file_path=path.join(dataset_path, \"data.json\"))\n",
        "noise_map = af.util.numpy_array_from_json(\n",
        "    file_path=path.join(dataset_path, \"noise_map.json\")\n",
        ")\n",
        "\n",
        "analysis = Analysis(data=data, noise_map=noise_map)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Customization__\n",
        "\n",
        "The `Analysis` class can be fully customized to be suitable for your model-fit.\n",
        "\n",
        "For example, additional inputs can be included in the `__init__` constructor and used in the `log_likelihood_function`.\n",
        "if they are required for your `log_likelihood_function` to work.\n",
        "\n",
        "The example below includes three additional inputs:\n",
        "\n",
        " - Instead of inputting a `noise_map`, a `noise_covariance_matrix` is input, which means that corrrlated noise is \n",
        "   accounted for in the `log_likelihood_function`.\n",
        " \n",
        " - A `mask` is input which masks the data such that certain data points are omitted from the log likelihood\n",
        " \n",
        " - A `kernel` is input which can account for certain blurring operations during data acquisition."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "\n",
        "\n",
        "class Analysis(af.Analysis):\n",
        "    def __init__(\n",
        "        self,\n",
        "        data: np.ndarray,\n",
        "        noise_covariance_matrix: np.ndarray,\n",
        "        mask: np.ndarray,\n",
        "        kernel: np.ndarray,\n",
        "    ):\n",
        "        \"\"\"\n",
        "        The `Analysis` class which has had its inputs edited for a different model-fit.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        data\n",
        "            A 1D numpy array containing the data (e.g. a noisy 1D signal) fitted in the workspace examples.\n",
        "        noise_covariance_matrix\n",
        "            A 2D numpy array containing the noise values and their covariances for the data, used for computing the\n",
        "            goodness of fit whilst accounting for correlated noise.\n",
        "        mask\n",
        "            A 1D numpy array containing a mask, where `True` values mean a data point is masked and is omitted from\n",
        "            the log likelihood.\n",
        "        kernel\n",
        "            A 1D numpy array containing the blurring kernel of the data, used for creating the model data.\n",
        "        \"\"\"\n",
        "        super().__init__()\n",
        "\n",
        "        self.data = data\n",
        "        self.noise_covariance_matrix = noise_covariance_matrix\n",
        "        self.mask = mask\n",
        "        self.kernel = kernel\n",
        "\n",
        "    def log_likelihood_function(self, instance) -> float:\n",
        "        \"\"\"\n",
        "        The `log_likelihood_function` now has access to the  `noise_covariance_matrix`, `mask` and `kernel`\n",
        "        input above.\n",
        "        \"\"\"\n",
        "        print(self.noise_covariance_matrix)\n",
        "        print(self.mask)\n",
        "        print(self.kernel)\n",
        "\n",
        "        \"\"\"\n",
        "        We do not provide a specific example of how to use these inputs in the `log_likelihood_function` as they are\n",
        "        specific to your model fitting problem.\n",
        "        \n",
        "        The key point is that any inputs required to compute the log likelihood can be passed into the `__init__`\n",
        "        constructor of the `Analysis` class and used in the `log_likelihood_function`.\n",
        "        \"\"\"\n",
        "\n",
        "        log_likelihood = None\n",
        "\n",
        "        return log_likelihood\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "An instance of the analysis class is created as follows."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "dataset_path = path.join(\"dataset\", \"example_1d\", \"gaussian_x1\")\n",
        "data = af.util.numpy_array_from_json(file_path=path.join(dataset_path, \"data.json\"))\n",
        "\n",
        "noise_covariance_matrix = np.ones(shape=(data.shape[0], data.shape[0]))\n",
        "mask = np.full(fill_value=False, shape=data.shape)\n",
        "kernel = np.full(fill_value=1.0, shape=data.shape)\n",
        "\n",
        "analysis = Analysis(\n",
        "    data=data, noise_covariance_matrix=noise_covariance_matrix, mask=mask, kernel=kernel\n",
        ")"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Visualization__\n",
        "\n",
        "If a `name` is input into a non-linear search, all results are output to hard-disk in a folder.\n",
        "\n",
        "By overwriting the `Visualizer` object of an `Analysis` class with a custom `Visualizer` class, custom results of the\n",
        "model-fit can be visualized during the model-fit.\n",
        "\n",
        "The `Visualizer` below has the methods `visualize_before_fit` and `visualize`, which perform model specific \n",
        "visualization will also be output into an `image` folder, for example as `.png` files.\n",
        "\n",
        "This uses the maximum log likelihood model of the model-fit inferred so far.\n",
        "\n",
        "Visualization of the results of the search, such as the corner plot of what is called the \"Probability Density \n",
        "Function\", are also automatically output during the model-fit on the fly."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "\n",
        "\n",
        "class Visualizer(af.Visualizer):\n",
        "    @staticmethod\n",
        "    def visualize_before_fit(\n",
        "        analysis, paths: af.DirectoryPaths, model: af.AbstractPriorModel\n",
        "    ):\n",
        "        \"\"\"\n",
        "        Before a model-fit, the `visualize_before_fit` method is called to perform visualization.\n",
        "\n",
        "        The function receives as input an instance of the `Analysis` class which is being used to perform the fit,\n",
        "        which is used to perform the visualization (e.g. it contains the data and noise map which are plotted).\n",
        "\n",
        "        This can output visualization of quantities which do not change during the model-fit, for example the\n",
        "        data and noise-map.\n",
        "\n",
        "        The `paths` object contains the path to the folder where the visualization should be output, which is determined\n",
        "        by the non-linear search `name` and other inputs.\n",
        "        \"\"\"\n",
        "\n",
        "        import matplotlib.pyplot as plt\n",
        "\n",
        "        xvalues = np.arange(analysis.data.shape[0])\n",
        "\n",
        "        plt.errorbar(\n",
        "            x=xvalues,\n",
        "            y=analysis.data,\n",
        "            yerr=analysis.noise_map,\n",
        "            linestyle=\"\",\n",
        "            color=\"k\",\n",
        "            ecolor=\"k\",\n",
        "            elinewidth=1,\n",
        "            capsize=2,\n",
        "        )\n",
        "        plt.title(\"Maximum Likelihood Fit\")\n",
        "        plt.xlabel(\"x value of profile\")\n",
        "        plt.ylabel(\"Profile Normalization\")\n",
        "        plt.savefig(path.join(paths.image_path, f\"data.png\"))\n",
        "        plt.clf()\n",
        "\n",
        "    @staticmethod\n",
        "    def visualize(analysis, paths: af.DirectoryPaths, instance, during_analysis):\n",
        "        \"\"\"\n",
        "        During a model-fit, the `visualize` method is called throughout the non-linear search.\n",
        "\n",
        "        The function receives as input an instance of the `Analysis` class which is being used to perform the fit,\n",
        "        which is used to perform the visualization (e.g. it generates the model data which is plotted).\n",
        "\n",
        "        The `instance` passed into the visualize method is maximum log likelihood solution obtained by the model-fit\n",
        "        so far and it can be used to provide on-the-fly images showing how the model-fit is going.\n",
        "\n",
        "        The `paths` object contains the path to the folder where the visualization should be output, which is determined\n",
        "        by the non-linear search `name` and other inputs.\n",
        "        \"\"\"\n",
        "        xvalues = np.arange(analysis.data.shape[0])\n",
        "\n",
        "        model_data = instance.model_data_from(xvalues=xvalues)\n",
        "        residual_map = analysis.data - model_data\n",
        "\n",
        "        \"\"\"\n",
        "        The visualizer now outputs images of the best-fit results to hard-disk (checkout `visualizer.py`).\n",
        "        \"\"\"\n",
        "        import matplotlib.pyplot as plt\n",
        "\n",
        "        plt.errorbar(\n",
        "            x=xvalues,\n",
        "            y=analysis.data,\n",
        "            yerr=analysis.noise_map,\n",
        "            linestyle=\"\",\n",
        "            color=\"k\",\n",
        "            ecolor=\"k\",\n",
        "            elinewidth=1,\n",
        "            capsize=2,\n",
        "        )\n",
        "        plt.plot(xvalues, model_data, color=\"r\")\n",
        "        plt.title(\"Maximum Likelihood Fit\")\n",
        "        plt.xlabel(\"x value of profile\")\n",
        "        plt.ylabel(\"Profile Normalization\")\n",
        "        plt.savefig(path.join(paths.image_path, f\"model_fit.png\"))\n",
        "        plt.clf()\n",
        "\n",
        "        plt.errorbar(\n",
        "            x=xvalues,\n",
        "            y=residual_map,\n",
        "            yerr=analysis.noise_map,\n",
        "            linestyle=\"\",\n",
        "            color=\"k\",\n",
        "            ecolor=\"k\",\n",
        "            elinewidth=1,\n",
        "            capsize=2,\n",
        "        )\n",
        "        plt.title(\"Residuals of Maximum Likelihood Fit\")\n",
        "        plt.xlabel(\"x value of profile\")\n",
        "        plt.ylabel(\"Residual\")\n",
        "        plt.savefig(path.join(paths.image_path, f\"model_fit.png\"))\n",
        "        plt.clf()\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The `Analysis` class is defined following the same API as before, but now with its `Visualizer` class attribute\n",
        "overwritten with the `Visualizer` class above."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "\n",
        "\n",
        "class Analysis(af.Analysis):\n",
        "\n",
        "    \"\"\"\n",
        "    This over-write means the `Visualizer` class is used for visualization throughout the model-fit.\n",
        "\n",
        "    This `VisualizerExample` object is in the `autofit.example.visualize` module and is used to customize the\n",
        "    plots output during the model-fit.\n",
        "\n",
        "    It has been extended with visualize methods that output visuals specific to the fitting of `1D` data.\n",
        "    \"\"\"\n",
        "\n",
        "    Visualizer = Visualizer\n",
        "\n",
        "    def __init__(self, data, noise_map):\n",
        "        \"\"\"\n",
        "        An Analysis class which illustrates visualization.\n",
        "        \"\"\"\n",
        "        super().__init__()\n",
        "\n",
        "        self.data = data\n",
        "        self.noise_map = noise_map\n",
        "\n",
        "    def log_likelihood_function(self, instance):\n",
        "        \"\"\"\n",
        "        The `log_likelihood_function` is identical to the example above\n",
        "        \"\"\"\n",
        "        xvalues = np.arange(self.data.shape[0])\n",
        "\n",
        "        model_data = instance.model_data_from(xvalues=xvalues)\n",
        "        residual_map = self.data - model_data\n",
        "        chi_squared_map = (residual_map / self.noise_map) ** 2.0\n",
        "        chi_squared = sum(chi_squared_map)\n",
        "        noise_normalization = np.sum(np.log(2 * np.pi * noise_map**2.0))\n",
        "        log_likelihood = -0.5 * (chi_squared + noise_normalization)\n",
        "\n",
        "        return log_likelihood\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Custom Result__\n",
        "\n",
        "The `Result` object is returned by a non-linear search after running the following code:\n",
        "\n",
        "`result = search.fit(model=model, analysis=analysis)`\n",
        "\n",
        "The result can be can be customized to include additional information about the model-fit that is specific to your \n",
        "model-fitting problem.\n",
        "\n",
        "For example, for fitting 1D profiles, the `Result` could include the maximum log likelihood model 1D data: \n",
        "\n",
        "`print(result.max_log_likelihood_model_data_1d)`\n",
        "\n",
        "In other examples, this quantity has been manually computed after the model-fit has completed.\n",
        "\n",
        "The custom result API allows us to do this. First, we define a custom `Result` class, which includes the property\n",
        "`max_log_likelihood_model_data_1d`."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "\n",
        "\n",
        "class ResultExample(af.Result):\n",
        "    @property\n",
        "    def max_log_likelihood_model_data_1d(self) -> np.ndarray:\n",
        "        \"\"\"\n",
        "        Returns the maximum log likelihood model's 1D model data.\n",
        "\n",
        "        This is an example of how we can pass the `Analysis` class a custom `Result` object and extend this result\n",
        "        object with new properties that are specific to the model-fit we are performing.\n",
        "        \"\"\"\n",
        "        xvalues = np.arange(self.analysis.data.shape[0])\n",
        "\n",
        "        return self.instance.model_data_from(xvalues=xvalues)\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The custom result has access to the analysis class, meaning that we can use any of its methods or properties to \n",
        "compute custom result properties.\n",
        "\n",
        "To make it so that the `ResultExample` object above is returned by the search we overwrite the `Result` class attribute \n",
        "of the `Analysis` and define a `make_result` object describing what we want it to contain:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "\n",
        "\n",
        "class Analysis(af.Analysis):\n",
        "\n",
        "    \"\"\"\n",
        "    This overwrite means the `ResultExample` class is returned after the model-fit.\n",
        "    \"\"\"\n",
        "\n",
        "    Result = ResultExample\n",
        "\n",
        "    def __init__(self, data, noise_map):\n",
        "        \"\"\"\n",
        "        An Analysis class which illustrates custom results.\n",
        "        \"\"\"\n",
        "        super().__init__()\n",
        "\n",
        "        self.data = data\n",
        "        self.noise_map = noise_map\n",
        "\n",
        "    def log_likelihood_function(self, instance):\n",
        "        \"\"\"\n",
        "        The `log_likelihood_function` is identical to the example above\n",
        "        \"\"\"\n",
        "        xvalues = np.arange(self.data.shape[0])\n",
        "\n",
        "        model_data = instance.model_data_from(xvalues=xvalues)\n",
        "        residual_map = self.data - model_data\n",
        "        chi_squared_map = (residual_map / self.noise_map) ** 2.0\n",
        "        chi_squared = sum(chi_squared_map)\n",
        "        noise_normalization = np.sum(np.log(2 * np.pi * noise_map**2.0))\n",
        "        log_likelihood = -0.5 * (chi_squared + noise_normalization)\n",
        "\n",
        "        return log_likelihood\n",
        "\n",
        "    def make_result(\n",
        "        self,\n",
        "        samples_summary: af.SamplesSummary,\n",
        "        paths: af.AbstractPaths,\n",
        "        samples: Optional[af.SamplesPDF] = None,\n",
        "        search_internal: Optional[object] = None,\n",
        "        analysis: Optional[object] = None,\n",
        "    ) -> Result:\n",
        "        \"\"\"\n",
        "        Returns the `Result` of the non-linear search after it is completed.\n",
        "\n",
        "        The result type is defined as a class variable in the `Analysis` class (see top of code under the python code\n",
        "        `class Analysis(af.Analysis)`.\n",
        "\n",
        "        The result can be manually overwritten by a user to return a user-defined result object, which can be extended\n",
        "        with additional methods and attribute specific to the model-fit.\n",
        "\n",
        "        This example class does example this, whereby the analysis result has been overwritten with the `ResultExample`\n",
        "        class, which contains a property `max_log_likelihood_model_data_1d` that returns the model data of the\n",
        "        best-fit model. This API means you can customize your result object to include whatever attributes you want\n",
        "        and therefore make a result object specific to your model-fit and model-fitting problem.\n",
        "\n",
        "        The `Result` object you return can be customized to include:\n",
        "\n",
        "        - The samples summary, which contains the maximum log likelihood instance and median PDF model.\n",
        "\n",
        "        - The paths of the search, which are used for loading the samples and search internal below when a search\n",
        "        is resumed.\n",
        "\n",
        "        - The samples of the non-linear search (e.g. MCMC chains) also stored in `samples.csv`.\n",
        "\n",
        "        - The non-linear search used for the fit in its internal representation, which is used for resuming a search\n",
        "        and making bespoke visualization using the search's internal results.\n",
        "\n",
        "        - The analysis used to fit the model (default disabled to save memory, but option may be useful for certain\n",
        "        projects).\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        samples_summary\n",
        "            The summary of the samples of the non-linear search, which include the maximum log likelihood instance and\n",
        "            median PDF model.\n",
        "        paths\n",
        "            An object describing the paths for saving data (e.g. hard-disk directories or entries in sqlite database).\n",
        "        samples\n",
        "            The samples of the non-linear search, for example the chains of an MCMC run.\n",
        "        search_internal\n",
        "            The internal representation of the non-linear search used to perform the model-fit.\n",
        "        analysis\n",
        "            The analysis used to fit the model.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        Result\n",
        "            The result of the non-linear search, which is defined as a class variable in the `Analysis` class.\n",
        "        \"\"\"\n",
        "        return self.Result(\n",
        "            samples_summary=samples_summary,\n",
        "            paths=paths,\n",
        "            samples=samples,\n",
        "            search_internal=search_internal,\n",
        "            analysis=self,\n",
        "        )\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For the sake of brevity, we do not run the code below, but the following code would work:\n",
        "\n",
        "`result = search.fit(model=model, analysis=analysis)`\n",
        "`print(result.max_log_likelihood_model_data_1d)`\n",
        "\n",
        "__Latent Variables__\n",
        "\n",
        "A latent variable is not a model parameter but can be derived from the model. Its value and errors may be of interest \n",
        "and aid in the interpretation of a model-fit. \n",
        "  \n",
        "For example, for the simple 1D Gaussian example, it could be the full-width half maximum (FWHM) of the Gaussian. \n",
        "This is not included in the model but can be easily derived from the Gaussian's sigma value.\n",
        "\n",
        "By overwriting the Analysis class's `compute_latent_variables` method we can manually specify latent variables that \n",
        "are calculated. If the search has a `name`, these are output to a `latent.csv` file, which mirrors \n",
        "the `samples.csv` file.\n",
        "  \n",
        "There may also be a `latent.results` and `latent_summary.json` files output. The `output.yaml` config file contains\n",
        "settings customizing what files are output and how often."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "\n",
        "\n",
        "class Analysis(af.Analysis):\n",
        "    def __init__(self, data, noise_map):\n",
        "        \"\"\"\n",
        "        An Analysis class which illustrates latent variables.\n",
        "        \"\"\"\n",
        "        super().__init__()\n",
        "\n",
        "        self.data = data\n",
        "        self.noise_map = noise_map\n",
        "\n",
        "    def log_likelihood_function(self, instance):\n",
        "        \"\"\"\n",
        "        The `log_likelihood_function` is identical to the example above\n",
        "        \"\"\"\n",
        "        xvalues = np.arange(self.data.shape[0])\n",
        "\n",
        "        model_data = instance.model_data_from(xvalues=xvalues)\n",
        "        residual_map = self.data - model_data\n",
        "        chi_squared_map = (residual_map / self.noise_map) ** 2.0\n",
        "        chi_squared = sum(chi_squared_map)\n",
        "        noise_normalization = np.sum(np.log(2 * np.pi * noise_map**2.0))\n",
        "        log_likelihood = -0.5 * (chi_squared + noise_normalization)\n",
        "\n",
        "        return log_likelihood\n",
        "\n",
        "    def compute_latent_variables(self, instance) -> Dict[str, float]:\n",
        "        \"\"\"\n",
        "        A latent variable is not a model parameter but can be derived from the model. Its value and errors may be\n",
        "        of interest and aid in the interpretation of a model-fit.\n",
        "\n",
        "        For example, for the simple 1D Gaussian example, it could be the full-width half maximum (FWHM) of the\n",
        "        Gaussian. This is not included in the model but can be easily derived from the Gaussian's sigma value.\n",
        "\n",
        "        By overwriting this method we can manually specify latent variables that are calculated and output to\n",
        "        a `latent.csv` file, which mirrors the `samples.csv` file.\n",
        "\n",
        "        In the example below, the `latent.csv` file will contain one column with the FWHM of every Gausian model\n",
        "        sampled by the non-linear search.\n",
        "\n",
        "        This function is called for every non-linear search sample, where the `instance` passed in corresponds to\n",
        "        each sample.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        instance\n",
        "            The instances of the model which the latent variable is derived from.\n",
        "\n",
        "        Returns\n",
        "        -------\n",
        "        A dictionary mapping every latent variable name to its value.\n",
        "\n",
        "        \"\"\"\n",
        "        return {\"fwhm\": instance.fwhm}\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Outputting latent variables manually after a fit is complete is simple, just call \n",
        "the `analysis.compute_latent_variables()` function. \n",
        "\n",
        "For many use cases, the best set up disables autofit latent variable output during a fit via the `output.yaml`\n",
        "file and perform it manually after completing a successful model-fit. This will save computational run time by not \n",
        "computing latent variables during a any model-fit which is unsuccessful."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "analysis = Analysis(data=data, noise_map=noise_map)\n",
        "\n",
        "# Commented out because we do not run the search in this cookbook\n",
        "\n",
        "# latent_samples = analysis.compute_latent_variables(samples=result.samples)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Analysing and interpreting latent variables is described fully in the result cookbook.\n",
        "\n",
        "However, in brief, the `latent_samples` object is a `Samples` object and uses the same API as samples objects."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# print(latent_samples.median_pdf().fwhm)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "__Custom Output__\n",
        "\n",
        "When performing fits which output results to hard-disc, a `files` folder is created containing .json / .csv files of \n",
        "the model, samples, search, etc.\n",
        "\n",
        "These files are human readable and help one quickly inspect and interpret results. \n",
        "\n",
        "By extending an `Analysis` class with the methods `save_attributes` and `save_results`, \n",
        "custom files can be written to the `files` folder to further aid this inspection. \n",
        "\n",
        "These files can then also be loaded via the database, as described in the database cookbook."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "\n",
        "\n",
        "class Analysis(af.Analysis):\n",
        "    def __init__(self, data: np.ndarray, noise_map: np.ndarray):\n",
        "        \"\"\"\n",
        "        Standard Analysis class example used throughout PyAutoFit examples.\n",
        "        \"\"\"\n",
        "        super().__init__()\n",
        "\n",
        "        self.data = data\n",
        "        self.noise_map = noise_map\n",
        "\n",
        "    def log_likelihood_function(self, instance) -> float:\n",
        "        \"\"\"\n",
        "        Standard log likelihood function used throughout PyAutoFit examples.\n",
        "        \"\"\"\n",
        "\n",
        "        xvalues = np.arange(self.data.shape[0])\n",
        "\n",
        "        model_data = instance.model_data_from(xvalues=xvalues)\n",
        "\n",
        "        residual_map = self.data - model_data\n",
        "        chi_squared_map = (residual_map / self.noise_map) ** 2.0\n",
        "        chi_squared = sum(chi_squared_map)\n",
        "        noise_normalization = np.sum(np.log(2 * np.pi * self.noise_map**2.0))\n",
        "        log_likelihood = -0.5 * (chi_squared + noise_normalization)\n",
        "\n",
        "        return log_likelihood\n",
        "\n",
        "    def save_attributes(self, paths: af.DirectoryPaths):\n",
        "        \"\"\"\n",
        "        Before the non-linear search begins, this routine saves attributes of the `Analysis` object to the `files`\n",
        "        folder such that they can be loaded after the analysis using PyAutoFit's database and aggregator tools.\n",
        "\n",
        "        For this analysis, it uses the `AnalysisDataset` object's method to output the following:\n",
        "\n",
        "        - The dataset's data as a .json file.\n",
        "        - The dataset's noise-map as a .json file.\n",
        "\n",
        "        These are accessed using the aggregator via `agg.values(\"data\")` and `agg.values(\"noise_map\")`.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        paths\n",
        "            The paths object which manages all paths, e.g. where the non-linear search outputs are stored,\n",
        "            visualization, and the pickled objects used by the aggregator output by this function.\n",
        "        \"\"\"\n",
        "        # The path where data.json is saved, e.g. output/dataset_name/unique_id/files/data.json\n",
        "\n",
        "        file_path = (path.join(paths._json_path, \"data.json\"),)\n",
        "\n",
        "        with open(file_path, \"w+\") as f:\n",
        "            json.dump(self.data, f, indent=4)\n",
        "\n",
        "        # The path where noise_map.json is saved, e.g. output/noise_mapset_name/unique_id/files/noise_map.json\n",
        "\n",
        "        file_path = (path.join(paths._json_path, \"noise_map.json\"),)\n",
        "\n",
        "        with open(file_path, \"w+\") as f:\n",
        "            json.dump(self.noise_map, f, indent=4)\n",
        "\n",
        "    def save_results(self, paths: af.DirectoryPaths, result: af.Result):\n",
        "        \"\"\"\n",
        "        At the end of a model-fit,  this routine saves attributes of the `Analysis` object to the `files`\n",
        "        folder such that they can be loaded after the analysis using PyAutoFit's database and aggregator tools.\n",
        "\n",
        "        For this analysis it outputs the following:\n",
        "\n",
        "        - The maximum log likelihood model data as a .json file.\n",
        "\n",
        "        This is accessed using the aggregator via `agg.values(\"model_data\")`.\n",
        "\n",
        "        Parameters\n",
        "        ----------\n",
        "        paths\n",
        "            The paths object which manages all paths, e.g. where the non-linear search outputs are stored,\n",
        "            visualization and the pickled objects used by the aggregator output by this function.\n",
        "        result\n",
        "            The result of a model fit, including the non-linear search, samples and maximum likelihood model.\n",
        "        \"\"\"\n",
        "        xvalues = np.arange(self.data.shape[0])\n",
        "\n",
        "        instance = result.max_log_likelihood_instance\n",
        "\n",
        "        model_data = instance.model_data_from(xvalues=xvalues)\n",
        "\n",
        "        # The path where model_data.json is saved, e.g. output/dataset_name/unique_id/files/model_data.json\n",
        "\n",
        "        file_path = (path.join(paths._json_path, \"model_data.json\"),)\n",
        "\n",
        "        with open(file_path, \"w+\") as f:\n",
        "            json.dump(model_data, f, indent=4)\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finish."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [],
      "outputs": [],
      "execution_count": null
    }
  ],
  "metadata": {
    "anaconda-cloud": {},
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}