rhayes777 / PyAutoFit

PyAutoFit: Classy Probabilistic Programming
https://pyautofit.readthedocs.io/
MIT License
59 stars 11 forks source link

Access `lower_limits_lists` and other quantities via model name #1026

Open Jammy2211 opened 1 month ago

Jammy2211 commented 1 month ago

A lower_limits_lists is passed into a GridSearchResult

class GridSearchResult:
    def __init__(
        self,
        samples: List[SamplesInterface],
        lower_limits_lists: Union[List, GridList],
        grid_priors: List[Prior],
        parent: Optional[NonLinearSearch] = None,
    ):

The lower_limits_lists is a list of list (of list of lists....) of the lower limit of the prior of every parameter that is searched over the grid search.

If the grid search has 3 dimensions, lower_limits_list will consist of 3 lists.

The quantities in lower_limits_list are used to compute other quantities, like the phyiscal value at the centre of each prior:

    @property
    @as_grid_list
    def physical_centres_lists(self) -> GridList:
        """
        The middle physical values for each grid square
        """
        return self._physical_values_for(self.centres_lists)

This is used in PyAutoLens to determine the centre of the cells where a subhalo search is performed, for visualization:

    def _array_2d_from(self, values) -> aa.Array2D:
        """
        Returns an `Array2D` where the input values are reshaped from list of lists to a 2D array, which is
        suitable for plotting.

        For example, this function may return the 2D array of the increases in log evidence for every lens model
        with a DM subhalo.

        The orientation of the 2D array and its values are chosen to ensure that when this array is plotted, DM
        subhalos with positive y and negative x `centre` coordinates appear in the top-left of the image.

        Parameters
        ----------
        values_native
            The list of list of values which are mapped to the 2D array (e.g. the `log_evidence` increases of every
            lens model with a DM subhalo).

        Returns
        -------
        The 2D array of values, where the values are mapped from the input list of lists.

        """
        values_reshaped = [value for values in values.native for value in values]

        return aa.Array2D.from_yx_and_values(
            y=[centre[0] for centre in self.physical_centres_lists],
            x=[centre[1] for centre in self.physical_centres_lists],
            values=values_reshaped,
            pixel_scales=self.physical_step_sizes,
            shape_native=self.shape,
        )

Up to now, this has always relied on a subhalo search being 2D (only being over the y and x coordinates). We may soon be doing higher dimensional grid searches and thus need to genelize the above method.

We basically just need a function which returns physical_centres_lists for the input model path, so something like:

centre_0_lists = self.physical_centres_lists_from("galaxies.subhalo.mass.centre.centre_0")

We also need support for multiple attribute inputs:

centre_lists = self.physical_centres_lists_from(["galaxies.subhalo.mass.centre.centre_0", "galaxies.subhalo.mass.centre.centre_1"])

It is crucial that the solution always maintains the structure of the physical_centres_lists used in the actual grid search.

However, the solution must work on all these different attribute,s as you can see above self.physical_step_sizes is also used in the function and would become ill defined if the grid search went to 3D.

Attribute Grid

A function exists in GridSearchResult which plays a similar role to the requested code above:

    @as_grid_list
    def attribute_grid(self, attribute_path: Union[str, Iterable[str]]) -> GridList:
        """
        Get a list of the attribute of the best instance from every search in a numpy array with the native dimensions
        of the grid search.

        Parameters
        ----------
        attribute_path
            The path to the attribute to get from the instance

        Returns
        -------
        A numpy array of the attribute of the best instance from every search in the grid search.
        """
        if isinstance(attribute_path, str):
            attribute_path = attribute_path.split(".")

        attribute_list = []
        for sample in self.samples:
            attribute = sample.instance
            for attribute_name in attribute_path:
                attribute = getattr(attribute, attribute_name)
            attribute_list.append(attribute)

        return attribute_list

Which is called via:

values=self.attribute_grid("galaxies.subhalo.mass.mass_at_200")

The difference for this function is that it returns the maximum log likelihood values of every input parameter based on model-path, as opposed to the values that define where the grid cell.

rhayes777 commented 1 month ago

We basically just need a function which returns physical_centres_lists for the input model path, so something like:

Are you talking about physical centres of just grid variables or any variable? Do we want this for the prior or posterior?

It is crucial that the solution always maintains the structure of the physical_centres_lists used in the actual grid search.

Presumably you mean that these values would be returned in an array the same shape as the grid search? From looking at the code it appears this is already supported for a grid search of arbitrary number of dimensions

Jammy2211 commented 1 month ago

Are you talking about physical centres of just grid variables or any variable? Do we want this for the prior or posterior?

I was only talking about grid variables, but having it work for any variable would be fine and actually preferable now that I think about it (and this appears to be what the current PR implementation does).

We want it to be for priors, e.g. the values used to set up the grid before the analysis runs. Doing the equivalent for the posterior requires the result.

rhayes777 commented 1 month ago

Cool ok so I think that's what this implementation does; as the model for each fit should depend on the grid limits the model alone should be sufficient. Asymmetric distributions are an exception as the mean will no longer be the centre of the grid square

Jammy2211 commented 1 month ago

Rather than using the mean of the samples, should the code not be using value_for to convert the values to physical values?

Something using this function I think?

    def _physical_values_for(self, unit_lists: GridList) -> List:
        """
        Compute physical values for lists of lists of unit hypercube
        values.

        Parameters
        ----------
        unit_lists
            A list of lists of hypercube values

        Returns
        -------
        A list of lists of physical values
        """
        return [
            [prior.value_for(limit) for prior, limit in zip(self.grid_priors, limits)]
            for limits in unit_lists
        ]

The physical values should be the centres of the priors, although maybe these both achieve the same thing.

rhayes777 commented 1 month ago

It could but that would be a much more complicated way of obtaining the same value for all except LogNormal priors

Jammy2211 commented 1 month ago

The subhalo mass will be a LogUniformPrior so autolens definitely neeeds it to work for that case.

rhayes777 commented 1 month ago

Cool will look at implementing it that way

rhayes777 commented 1 month ago

At the moment the grid priors are all converted to UniformPrior to peform the grid search.


    def make_arguments(self, values, grid_priors):
        arguments = {}
        for value, grid_prior in zip(values, grid_priors):
            if (
                float("-inf") == grid_prior.lower_limit
                or float("inf") == grid_prior.upper_limit
            ):
                raise exc.PriorException(
                    "Priors passed to the grid search must have definite limits"
                )
            lower_limit = grid_prior.lower_limit + value * grid_prior.width
            upper_limit = (
                grid_prior.lower_limit + (value + self.step_size) * grid_prior.width
            )
            prior = p.UniformPrior(lower_limit=lower_limit, upper_limit=upper_limit)
            arguments[grid_prior] = prior
        return arguments
rhayes777 commented 1 month ago

Also it looks like the implementation of mean is the physical value at half unit value:


>>> import autofit as af
>>> af.UniformPrior(1, 2).mean
1.5

>>> af.LogUniformPrior(1, 2).mean
1.4142135623730951

>>> af.LogUniformPrior(lower_limit=1.0, upper_limit=2.0).value_for(0.5)
1.4142135623730951
Jammy2211 commented 1 month ago

Ok, that will be why it works I guess.