astierCovFitParameters.py

Go to the documentation of this file.

# This file is part of cp_pipe.
#
# Developed for the LSST Data Management System.
# This product includes software developed by the LSST Project
# (https://www.lsst.org).
# See the COPYRIGHT file at the top-level directory of this distribution
# for details of code ownership.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
# Code developed at LPNHE ("saunerie", Paris)
 
import numpy as np
 
"""
The classes in this file were taken from https://github.com/PierreAstier/bfptc
by Pierre Astier (Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE),
Sorbonne Université, Paris, France).
 
File: bfptc/py/fitparameters.py
Commit hash: d46ba836fd5feb1c0065b61472c5f31b73b8480f
 
 
Notes from original code
-----------------------
 
This module contains utility classes to handle parameters in linear and
non-linear least square fits implemented in linearmodels and
nonlinearmodels. It provide 2 main features:
 
- `StructArray`: a derivative of numpy class ndarray to manage large
vectors organized into named subgroups.
 
- `FitParameters` : a class to manage large parameter vectors.  It
allows to easily fix/release specific parameters or entire subgroups,
and remap the remaining free parameters into a contiguous vector.
"""
 
 
class Structure(object):
    """ Collection of named slices
 
    Slices are specified by a name and a length. If omitted, the length
    default to one and the name to __0__ for the first unamed slice, __1__
    for the second and so on.
 
    Examples:
    ---------
    >>> s = Structure([('a', 7), 3])
    >>> print len(s)
    10
    >>> for name in s: print name
    a
    __0__
    >>> print len(Structure([('a', 3), 'b']))
    4
    """
    def __init__(self, groups):
        if isinstance(groups, Structure):
            self.slicesslices = groups.slices.copy()
            self._len_len = groups._len
        else:
            self.slicesslices = {}
            i = 0
            _n_unnamed = 0
            for group in groups:
                if isinstance(group, int):
                    name = "__%d__" % _n_unnamed
                    _len = group
                    _n_unnamed += 1
                elif isinstance(group, str):
                    name = group
                    _len = 1
                else:
                    try:
                        name, _len = group
                    except TypeError:
                        raise TypeError('Structure specification not understood: %s' % repr(group))
                self.slicesslices[name] = slice(i, i + _len)
                i += _len
            self._len_len = i
 
    def __getitem__(self, arg):
        if isinstance(arg, str):
            return self.slicesslices[arg]
        else:
            return arg
 
    def __iter__(self):
        return self.slicesslices.__iter__()
 
    def __len__(self):
        return self._len_len
 
 
class StructArray(np.ndarray):
    """Decorate numpy arrays with a collection of named slices.
 
    Array slices becomes accessible by their name. This is applicable to
    nd array, although the same `Structure` is shared between all
    dimensions.
 
    Examples:
    ---------
    >>> v = StructArray(np.zeros(10), [('a', 3), ('b', 7)])
    >>> print v['a']
    [ 0.  0.  0.]
 
    >>> C = StructArray(np.zeros((10,10)), [('a', 2), ('b', 8)])
    >>> print C['a', 'a']
    [[ 0.  0.]
     [ 0.  0.]]
    """
    def __new__(cls, array, struct=[]):
        obj = array.view(cls)
        obj.struct = Structure(struct)
        return obj
 
    def __array_finalize__(self, obj):
        if obj is None:
            return
        self.structstruct = getattr(obj, 'struct', [])
 
    def __getitem__(self, args):
        if not isinstance(args, tuple):
            args = args,
        newargs = tuple([self.structstruct[arg] for arg in args])
        return np.asarray(self)[newargs]
 
    def __setitem__(self, args, val):
        if not isinstance(args, tuple):
            args = args,
        newargs = tuple([self.structstruct[arg] for arg in args])
        np.asarray(self)[newargs] = val
 
    def __getslice__(self, start, stop):
        return self.__getitem____getitem__(slice(start, stop))
 
    def __reduce__(self):
        # Get the parent's __reduce__ tuple
        pickled_state = super(StructArray, self).__reduce__()
        # Create our own tuple to pass to __setstate__
        new_state = pickled_state[2] + (self.structstruct,)
        # Return a tuple that replaces the parent's __setstate__ tuple with our own
        return (pickled_state[0], pickled_state[1], new_state)
 
    def __setstate__(self, state):
        self.structstruct = state[-1]  # Set the info attribute
        # Call the parent's __setstate__ with the other tuple elements.
        super(StructArray, self).__setstate__(state[0:-1])
 
 
class FitParameters(object):
    """ Manages a vector of fit parameters with the possibility to mark a subset of
    them as fixed to a given value.
 
    The parameters can be organized in named slices (block of contiguous
    values) accessible through indexing by their name as in `StructArray`.
 
    >>> p_salt = FitParameters(['X0', 'X1', 'Color', 'Redshift'])
    >>> p_dice = FitParameters([('alpha', 2), ('S', 10), ('dSdT', 10), 'idark'])
 
    It is possible to modify the parameters in place. Using the indexing
    of slices by name simplifies somewhat the operations, as one does
    not need to care about the position of a slice within the entire
    parameter vector:
 
    >>> p_dice['idark'][0] = -1.0232E-12
    >>> p_dice['S'][4] = 25.242343E-9
    >>> p_dice['dSdT'][:] = 42.
    >>> print p_dice
    alpha: array([ 0.,  0.])
    S: array([  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
             0.00000000e+00,   2.52423430e-08,   0.00000000e+00,
             0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
             0.00000000e+00])
    dSdT: array([ 42.,  42.,  42.,  42.,  42.,  42.,  42.,  42.,  42.,  42.])
    idark: array([ -1.02320000e-12])
 
    It is also possible to mark parameters as fixed to a value.
    >>> p_dice.fix(0, 12.)
 
    Value is optional. The above is equivalent to:
    >>> p_dice[0] = 12.
    >>> p_dice.fix(0)
 
    Again named slices simplifies the operations:
    >>> p_dice['S'].fix([0, -1], 12.)
    >>> p_dice['dSdT'].fix([0, -1])
 
    One can fix entire slices at once:
    >>> p_dice['idark'].fix()
    >>> p_salt['Redshift'].fix(val=0.23432)
 
    The property ``full'' give access to the vector of parameters. The
    property "free" gives access to the free parameters:
    >>> print len(p_dice.free), len(p_dice.full)
    17 23
 
    Note that free relies on fancy indexing. Access thus trigger a
    copy. As a consequence, the following will not actually alter the
    data:
    >>> p_dice.free[0] = 12.
    >>> print p_dice.free[0]
    0.0
 
    It is still possible to set slices of free parameters as a
    contiguous vector. For example:
    >>> p_dice['S'].free = 12.
    >>> print p_dice['S'].free
    [ 12.  12.  12.  12.  12.  12.  12.  12.]
 
    >>> p_dice[:5].free = 4.
    >>> print p_dice[:5].free
    [ 4.  4.  4.]
 
    In particular, the typical use case which consists in updating the
    free parameters with the results of a fit works as expected:
    >>> p = np.arange(len(p_dice.free))
    >>> p_dice.free = p
 
    Last the class provide a convenience function that return the index
    of a subset of parameters in the global free parameters vector, and
    -1 for fixed parameters:
    >>> print p_dice['dSdT'].indexof()
    [-1  9 10 11 12 13 14 15 16 -1]
    >>> print p_dice['dSdT'].indexof([1,2])
    [ 9 10]
    """
    def __init__(self, groups):
        struct = Structure(groups)
        self._free_free = StructArray(np.ones(len(struct), dtype='bool'), struct)
        self._pars_pars = StructArray(np.zeros(len(struct), dtype='float'), struct)
        self._index_index = StructArray(np.zeros(len(struct), dtype='int'), struct)
        self._struct_struct = struct
        self._reindex_reindex()
        self._base_base = self
 
    def copy(self):
        cop = FitParameters(self._struct_struct)
        cop._free = StructArray(np.copy(self._free_free), cop._struct)
        cop._pars = StructArray(np.copy(self._pars_pars), cop._struct)
        cop._index = StructArray(np.copy(self._index_index), cop._struct)
        cop._reindex()
        cop._base = cop
        return cop
 
    def _reindex(self):
        self._index_index[self._free_free] = np.arange(self._free_free.sum())
        self._index_index[~self._free_free] = -1
 
    def fix(self, keys=slice(None), val=None):
        self._free_free[keys] = False
        if val is not None:
            self._pars_pars[keys] = val
        self._base_base._reindex()
 
    def release(self, keys=slice(None)):
        self._free_free[keys] = True
        self._base_base._reindex()
 
    def indexof(self, indices=slice(None)):
        return self._index_index[indices]
 
    def __getitem__(self, args):
        # Prevent conversion to scalars
        if isinstance(args, int):
            args = slice(args, args + 1)
        new = FitParameters.__new__(FitParameters)
        new._free = self._free_free[args]
        new._pars = self._pars_pars[args]
        new._index = self._index_index[args]
        new._base = self._base_base
        return new
 
    def __setitem__(self, args, val):
        self._pars_pars[args] = val
 
    @property
    def free(self):
        return self._pars_pars[self._free_free]
 
    @free.setter
    def free(self, val):
        self._pars_pars[self._free_free] = val
 
    @property
    def full(self):
        return self._pars_pars
 
    @full.setter
    def full(self, val):
        self._pars_pars = val
 
    def __repr__(self):
        if hasattr(self, '_struct'):
            s = "\n".join(['%s: %s' % (key, repr(self._pars_pars[key])) for key in self._struct_struct])
        else:
            s = repr(self._pars_pars)
        return s