src/qinfer/abstract_model.py
#!/usr/bin/python
# -*- coding: utf-8 -*-
##
# abstract_model.py: Abstract interfaces for models with different levels of
# functionality.
##
# © 2017, Chris Ferrie (csferrie@gmail.com) and
# Christopher Granade (cgranade@cgranade.com).
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# 1. Redistributions of source code must retain the above copyright
# notice, this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# 3. Neither the name of the copyright holder nor the names of its
# contributors may be used to endorse or promote products derived from
# this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
# POSSIBILITY OF SUCH DAMAGE.
##
## FEATURES ##################################################################
from __future__ import absolute_import
from __future__ import division, unicode_literals
## EXPORTS ###################################################################
__all__ = [
'Simulatable',
'Model',
'FiniteOutcomeModel',
'DifferentiableModel'
]
## IMPORTS ###################################################################
from builtins import range, map
from future.utils import with_metaclass
import abc
# Python standard library package for specifying abstract classes.
import numpy as np
import warnings
from qinfer.utils import safe_shape
from qinfer.domains import IntegerDomain
## CLASSES ###################################################################
class Simulatable(with_metaclass(abc.ABCMeta, object)):
r"""
Represents a system which can be simulated according to
various model parameters and experimental control parameters
in order to produce representative data.
See :ref:`models_guide` for more details.
:param bool allow_identical_outcomes: Whether the method ``outcomes`` should
be allowed to return multiple identical outcomes for a given ``expparam``.
It will be more efficient to set this to ``True`` whenever it is likely
that multiple identical outcomes will occur.
"""
def __init__(self):
"""
:param bool always_resample_outcomes: Resample outcomes stochastically with
each outcome call.
:param :class:`~numpy.ndarray` initial_outcomes: Initial set of outcomes
that may be supplied. Otherwise initial outcomes default to
zeros.
"""
self._sim_count = 0
# Initialize a default scale matrix.
self._Q = np.ones((self.n_modelparams,))
## ABSTRACT PROPERTIES ##
@abc.abstractproperty
def n_modelparams(self):
"""
Returns the number of real model parameters admitted by this model.
This property is assumed by inference engines to be constant for
the lifetime of a :class:`Model` instance.
"""
pass
@abc.abstractproperty
def expparams_dtype(self):
"""
Returns the dtype of an experiment parameter array. For a
model with single-parameter control, this will likely be a scalar dtype,
such as ``"float64"``. More generally, this can be an example of a
record type, such as ``[('time', py.'float64'), ('axis', 'uint8')]``.
This property is assumed by inference engines to be constant for
the lifetime of a Model instance.
"""
pass
## CONCRETE PROPERTIES ##
@property
def is_n_outcomes_constant(self):
"""
Returns ``True`` if and only if both the domain and ``n_outcomes``
are independent of the expparam.
This property is assumed by inference engines to be constant for
the lifetime of a Model instance.
"""
return True
@property
def model_chain(self):
"""
Returns a tuple of models upon which this model is based,
such that properties and methods of underlying models for
models that decorate other models can be accessed. For a
standalone model, this is always the empty tuple.
"""
return ()
@property
def base_model(self):
"""
Returns the most basic model that this model depends on.
For standalone models, this property satisfies ``model.base_model is model``.
"""
return self
@property
def underlying_model(self):
"""
Returns the model that this model is based on (decorates)
if such a model exists, or ``None`` if this model is
independent.
"""
return self.model_chain[-1] if self.model_chain else None
@property
def sim_count(self):
"""
Returns the number of data samples that have been produced by
this simulator.
:rtype: int
"""
return self._sim_count
@property
def Q(self):
r"""
Returns the diagonal of the scale matrix :math:`\matr{Q}` that
relates the scales of each of the model parameters. In particular,
the quadratic loss for this Model is defined as:
.. math::
L_{\matr{Q}}(\vec{x}, \hat{\vec{x}}) =
(\vec{x} - \hat{\vec{x}})^\T \matr{Q} (\vec{x} - \hat{\vec{x}})
If a subclass does not explicitly define the scale matrix, it is taken
to be the identity matrix of appropriate dimension.
:return: The diagonal elements of :math:`\matr{Q}`.
:rtype: :class:`~numpy.ndarray` of shape ``(n_modelparams, )``.
"""
return self._Q
@property
def modelparam_names(self):
"""
Returns the names of the various model parameters admitted by this
model, formatted as LaTeX strings.
"""
return list(map("x_{{{}}}".format, range(self.n_modelparams)))
## CONCRETE METHODS ##
def _repr_html_(self, suppress_base=False):
s = r"""
<strong>{type.__name__}</strong> at 0x{id:0x}: {n_mp} model parameter{plural}
""".format(
id=id(self), type=type(self),
n_mp=self.n_modelparams,
plural="" if self.n_modelparams == 1 else "s"
)
if not suppress_base and self.model_chain:
s += r"""<br>
<p>Model chain:</p>
<ul>{}
</ul>
""".format(r"\n".join(
u"<li>{}</li>".format(model._repr_html_(suppress_base=True))
for model in reversed(self.model_chain)
))
return s
def are_expparam_dtypes_consistent(self, expparams):
"""
Returns ``True`` iff all of the given expparams
correspond to outcome domains with the same dtype.
For efficiency, concrete subclasses should override this method
if the result is always ``True``.
:param np.ndarray expparams: Array of expparamms
of type ``expparams_dtype``
:rtype: ``bool``
"""
if self.is_n_outcomes_constant:
# This implies that all domains are equal, so this must be true
return True
# otherwise we have to actually check all the dtypes
if expparams.size > 0:
domains = self.domain(expparams)
first_dtype = domains[0].dtype
return all(domain.dtype == first_dtype for domain in domains[1:])
else:
return True
## ABSTRACT METHODS ##
@abc.abstractmethod
def n_outcomes(self, expparams):
"""
Returns an array of dtype ``uint`` describing the number of outcomes
for each experiment specified by ``expparams``.
If the number of outcomes does not depend on expparams (i.e.
``is_n_outcomes_constant`` is ``True``), this method
should return a single number.
If there are an infinite (or intractibly large) number of outcomes,
this value specifies the number of outcomes to randomly sample.
:param numpy.ndarray expparams: Array of experimental parameters. This
array must be of dtype agreeing with the ``expparams_dtype``
property.
"""
pass
@abc.abstractmethod
def domain(self, exparams):
"""
Returns a list of :class:`Domain` objects, one for each input expparam.
:param numpy.ndarray expparams: Array of experimental parameters. This
array must be of dtype agreeing with the ``expparams_dtype``
property, or, in the case where ``n_outcomes_constant`` is ``True``,
``None`` should be a valid input.
:rtype: list of :class:`Domain`
"""
pass
@abc.abstractmethod
def are_models_valid(self, modelparams):
"""
Given a shape ``(n_models, n_modelparams)`` array of model parameters,
returns a boolean array of shape ``(n_models)`` specifying whether
each set of model parameters represents is valid under this model.
"""
pass
@abc.abstractmethod
def simulate_experiment(self, modelparams, expparams, repeat=1):
"""
Produces data according to the given model parameters and experimental
parameters, structured as a NumPy array.
:param np.ndarray modelparams: A shape ``(n_models, n_modelparams)``
array of model parameter vectors describing the hypotheses under
which data should be simulated.
:param np.ndarray expparams: A shape ``(n_experiments, )`` array of
experimental control settings, with ``dtype`` given by
:attr:`~qinfer.Model.expparams_dtype`, describing the
experiments whose outcomes should be simulated.
:param int repeat: How many times the specified experiment should
be repeated.
:rtype: np.ndarray
:return: A three-index tensor ``data[i, j, k]``, where ``i`` is the repetition,
``j`` indexes which vector of model parameters was used, and where
``k`` indexes which experimental parameters where used. If ``repeat == 1``,
``len(modelparams) == 1`` and ``len(expparams) == 1``, then a scalar
datum is returned instead.
"""
self._sim_count += modelparams.shape[0] * expparams.shape[0] * repeat
assert(self.are_expparam_dtypes_consistent(expparams))
## CONCRETE METHODS ##
def clear_cache(self):
"""
Tells the model to clear any internal caches used in computing
likelihoods and drawing samples. Calling this method should not cause
any different results, but should only affect performance.
"""
# By default, no cache to clear.
pass
def experiment_cost(self, expparams):
"""
Given an array of experimental parameters, returns the cost associated
with performing each experiment. By default, this cost is constant
(one) for every experiment.
:param expparams: An array of experimental parameters for which the cost
is to be evaluated.
:type expparams: :class:`~numpy.ndarray` of ``dtype`` given by
:attr:`~Model.expparams_dtype`
:return: An array of costs corresponding to the specified experiments.
:rtype: :class:`~numpy.ndarray` of ``dtype`` ``float`` and of the
same shape as ``expparams``.
"""
return np.ones(expparams.shape)
def distance(self, a, b):
r"""
Gives the distance between two model parameter vectors :math:`\vec{a}` and
:math:`\vec{b}`. By default, this is the vector 1-norm of the difference
:math:`\mathbf{Q} (\vec{a} - \vec{b})` rescaled by
:attr:`~Model.Q`.
:param np.ndarray a: Array of model parameter vectors having shape
``(n_models, n_modelparams)``.
:param np.ndarray b: Array of model parameters to compare to, having
the same shape as ``a``.
:return: An array ``d`` of distances ``d[i]`` between ``a[i, :]`` and
``b[i, :]``.
"""
return np.apply_along_axis(
lambda vec: np.linalg.norm(vec, 1),
1,
self.Q * (a - b)
)
def update_timestep(self, modelparams, expparams):
r"""
Returns a set of model parameter vectors that is the update of an
input set of model parameter vectors, such that the new models are
conditioned on a particular experiment having been performed.
By default, this is the trivial function
:math:`\vec{x}(t_{k+1}) = \vec{x}(t_k)`.
:param np.ndarray modelparams: Set of model parameter vectors to be
updated.
:param np.ndarray expparams: An experiment parameter array describing
the experiment that was just performed.
:return np.ndarray: Array of shape
``(n_models, n_modelparams, n_experiments)`` describing the update
of each model according to each experiment.
"""
return np.tile(modelparams, (expparams.shape[0],1,1)).transpose((1,2,0))
def canonicalize(self, modelparams):
r"""
Returns a canonical set of model parameters corresponding to a given
possibly non-canonical set. This is used for models in which there
exist model parameters :math:`\vec{x}_i` and :\math:`\vec{x}_j` such
that
.. math::
\Pr(d | \vec{x}_i; \vec{e}) = \Pr(d | \vec{x}_j; \vec{e})
for all outcomes :math:`d` and experiments :math:`\vec{e}`. For
models admitting such an ambiguity, this
method should then be overridden to return a consistent choice
out of such vectors, hence avoiding supurious model degeneracies.
Note that, by default, :class:`~qinfer.smc.SMCUpdater` will *not*
call this method.
"""
return modelparams
class Model(Simulatable):
"""
Represents a system which can be simulated according to
various model parameters and experimental control parameters
in order to produce the probability of a hypothetical data
record. As opposed to :class:`~qinfer.Simulatable`, instances
of :class:`~qinfer.Model` not only produce data consistent
with the description of a system, but also evaluate the probability
of that data arising from the system.
:param bool allow_identical_outcomes: Whether the method
``representative_outcomes`` should be allowed to return multiple
identical outcomes for a given ``expparam``.
:param float outcome_warning_threshold: Threshold value below which
``representative_outcomes``
will issue a warning about the representative outcomes
not adequately covering the domain with respect to
the relevant distribution.
See :ref:`models_guide` for more details.
"""
## INITIALIZERS ##
def __init__(self, allow_identical_outcomes=False,
outcome_warning_threshold=0.99):
super(Model, self).__init__()
self._call_count = 0
self._allow_identical_outcomes = allow_identical_outcomes
self._outcome_warning_threshold = outcome_warning_threshold
## CONCRETE PROPERTIES ##
@property
def call_count(self):
"""
Returns the number of points at which
the probability of this model has been evaluated,
where a point consists of a hypothesis about the model (a vector of
model parameters), an experimental control setting (expparams) and
a hypothetical or actual datum.
:rtype: int
"""
return self._call_count
## ABSTRACT METHODS ##
@abc.abstractmethod
def likelihood(self, outcomes, modelparams, expparams):
r"""
Calculates the probability of each given outcome, conditioned on each
given model parameter vector and each given experimental control setting.
:param np.ndarray modelparams: A shape ``(n_models, n_modelparams)``
array of model parameter vectors describing the hypotheses for
which the likelihood function is to be calculated.
:param np.ndarray expparams: A shape ``(n_experiments, )`` array of
experimental control settings, with ``dtype`` given by
:attr:`~qinfer.Simulatable.expparams_dtype`, describing the
experiments from which the given outcomes were drawn.
:rtype: np.ndarray
:return: A three-index tensor ``L[i, j, k]``, where ``i`` is the outcome
being considered, ``j`` indexes which vector of model parameters was used,
and where ``k`` indexes which experimental parameters where used.
Each element ``L[i, j, k]`` then corresponds to the likelihood
:math:`\Pr(d_i | \vec{x}_j; e_k)`.
"""
# Count the number of times the inner-most loop is called.
self._call_count += (
safe_shape(outcomes) * safe_shape(modelparams) * safe_shape(expparams)
)
@property
def allow_identical_outcomes(self):
"""
Whether the method ``representative_outcomes`` should be allowed to return multiple
identical outcomes for a given ``expparam``.
It will be more efficient to set this to ``True`` whenever it is likely
that multiple identical outcomes will occur.
:return: Flag state.
:rtype: ``bool``
"""
return self._allow_identical_outcomes
@allow_identical_outcomes.setter
def allow_identical_outcomes(self, value):
"""
Whether the method ``representative_outcomes`` should be allowed to return multiple
identical outcomes for a given ``expparam``.
It will be more efficient to set this to ``True`` whenever it is likely
that multiple identical outcomes will occur.
:param bool allow_identical_outcomes: Value of flag.
"""
self._allow_identical_outcomes = value
@property
def outcome_warning_threshold(self):
"""
Threshold value below which ``representative_outcomes``
will issue a warning about the representative outcomes
not adequately covering the domain with respect to
the relevant distribution.
:return: Threshold value.
:rtype: ``float``
"""
return self._outcome_warning_threshold
@outcome_warning_threshold.setter
def outcome_warning_threshold(self, value):
"""
Threshold value below which ``representative_outcomes``
will issue a warning about the representative outcomes
not adequately covering the domain with respect to
the relevant distribution.
:param float value: Threshold value.
"""
self._outcome_warning_threshold = value
## CONCRETE METHODS ##
# These methods depend on the abstract methods, and thus their behaviors
# change in each inheriting class.
def is_model_valid(self, modelparams):
"""
Returns True if and only if the model parameters given are valid for
this model.
"""
return self.are_models_valid(modelparams[np.newaxis, :])[0]
class LinearCostModelMixin(Model):
# FIXME: move this mixin to a new module.
# TODO: test this mixin.
"""
This mixin implements :meth:`Model.experiment_cost` by setting the
cost of an experiment equal to the value of a given field of each
``expparams`` element (by default, ``t``).
"""
_field = "t"
def experiment_cost(self, expparams):
return expparams[self._field]
class FiniteOutcomeModel(Model):
"""
Represents a system in the same way that :class:`~qinfer.Model`,
except that it is demanded that the number of outcomes for any
experiment be known and finite.
:param bool allow_identical_outcomes: Whether the method
``representative_outcomes`` should be allowed to return multiple
identical outcomes for a given ``expparam``.
:param float outcome_warning_threshold: Threshold value below which
``representative_outcomes``
will issue a warning about the representative outcomes
not adequately covering the domain with respect to
the relevant distribution.
:param int n_outcomes_cutoff: If ``n_outcomes`` exceeds this value,
``representative_outcomes`` will use this
value in its place. This is useful in the case
of a finite yet untractible number of outcomes. Use ``None``
for no cutoff.
See :class:`~qinfer.Model` and :ref:`models_guide` for more details.
"""
## INITIALIZERS ##
def __init__(self,
allow_identical_outcomes=False,
outcome_warning_threshold=0.99,
n_outcomes_cutoff=None
):
super(FiniteOutcomeModel, self).__init__(
outcome_warning_threshold=outcome_warning_threshold,
allow_identical_outcomes=allow_identical_outcomes)
self._n_outcomes_cutoff = n_outcomes_cutoff
if self.is_n_outcomes_constant:
# predefine if we can
self._domain = IntegerDomain(min=0,max=self.n_outcomes(None)-1)
## CONCRETE PROPERTIES ##
@property
def n_outcomes_cutoff(self):
"""
If ``n_outcomes`` exceeds this value for
some expparm, ``representative_outcomes`` will use this
value in its place. This is useful in the case
of a finite yet untractible number of outcomes.
:return: Cutoff value.
:rtype: ``int``
"""
return self._n_outcomes_cutoff
@n_outcomes_cutoff.setter
def n_outcomes_cutoff(self, value):
"""
If ``n_outcomes`` exceeds this value,
``representative_outcomes`` will use this
value in its place. This is useful in the case
of a finite yet untractible number of outcomes.
:param int value: Cutoff value.
"""
self.n_outcomes_cutoff = value
## ABSTRACT METHODS ##
## CONCRETE METHODS ##
# These methods depend on the abstract methods, and thus their behaviors
# change in each inheriting class.
def domain(self, expparams):
"""
Returns a list of :class:`Domain` objects, one for each input expparam.
:param numpy.ndarray expparams: Array of experimental parameters. This
array must be of dtype agreeing with the ``expparams_dtype``
property, or, in the case where ``n_outcomes_constant`` is ``True``,
``None`` should be a valid input.
:rtype: list of ``Domain``
"""
# As a convenience to most users, we define domain for them. If a
# fancier domain is desired, this method can easily be overridden.
if self.is_n_outcomes_constant:
return self._domain if expparams is None else [self._domain for ep in expparams]
else:
return [IntegerDomain(min=0,max=n_o-1) for n_o in self.n_outcomes(expparams)]
def simulate_experiment(self, modelparams, expparams, repeat=1):
# Call the superclass simulate_experiment, not recording the result.
# This is used to count simulation calls.
super(FiniteOutcomeModel, self).simulate_experiment(modelparams, expparams, repeat)
if self.is_n_outcomes_constant:
# In this case, all expparams have the same domain
all_outcomes = self.domain(None).values
probabilities = self.likelihood(all_outcomes, modelparams, expparams)
cdf = np.cumsum(probabilities, axis=0)
randnum = np.random.random((repeat, 1, modelparams.shape[0], expparams.shape[0]))
outcome_idxs = all_outcomes[np.argmax(cdf > randnum, axis=1)]
outcomes = all_outcomes[outcome_idxs]
else:
# Loop over each experiment, sadly.
# Assume all domains have the same dtype
assert(self.are_expparam_dtypes_consistent(expparams))
dtype = self.domain(expparams[0, np.newaxis])[0].dtype
outcomes = np.empty((repeat, modelparams.shape[0], expparams.shape[0]), dtype=dtype)
for idx_experiment, single_expparams in enumerate(expparams[:, np.newaxis]):
all_outcomes = self.domain(single_expparams).values
probabilities = self.likelihood(all_outcomes, modelparams, single_expparams)
cdf = np.cumsum(probabilities, axis=0)[..., 0]
randnum = np.random.random((repeat, 1, modelparams.shape[0]))
outcomes[:, :, idx_experiment] = all_outcomes[np.argmax(cdf > randnum, axis=1)]
return outcomes[0, 0, 0] if repeat == 1 and expparams.shape[0] == 1 and modelparams.shape[0] == 1 else outcomes
## STATIC METHODS ##
# These methods are provided as a convienence to make it easier to write
# simple models.
@staticmethod
def pr0_to_likelihood_array(outcomes, pr0):
"""
Assuming a two-outcome measurement with probabilities given by the
array ``pr0``, returns an array of the form expected to be returned by
``likelihood`` method.
:param numpy.ndarray outcomes: Array of integers indexing outcomes.
:param numpy.ndarray pr0: Array of shape ``(n_models, n_experiments)``
describing the probability of obtaining outcome ``0`` from each
set of model parameters and experiment parameters.
"""
pr0 = pr0[np.newaxis, ...]
pr1 = 1 - pr0
if len(np.shape(outcomes)) == 0:
outcomes = np.array(outcomes)[None]
return np.concatenate([
pr0 if outcomes[idx] == 0 else pr1
for idx in range(safe_shape(outcomes))
])
class DifferentiableModel(with_metaclass(abc.ABCMeta, Model)):
@abc.abstractmethod
def score(self, outcomes, modelparams, expparams, return_L=False):
r"""
Returns the score of this likelihood function, defined as:
.. math::
q(d, \vec{x}; \vec{e}) = \vec{\nabla}_{\vec{x}} \log \Pr(d | \vec{x}; \vec{e}).
Calls are represented as a four-index tensor
``score[idx_modelparam, idx_outcome, idx_model, idx_experiment]``.
The left-most index may be suppressed for single-parameter models.
If return_L is True, both `q` and the likelihood `L` are returned as `q, L`.
"""
pass
def fisher_information(self, modelparams, expparams):
"""
Returns the covariance of the score taken over possible outcomes,
known as the Fisher information.
The result is represented as the four-index tensor
``fisher[idx_modelparam_i, idx_modelparam_j, idx_model, idx_experiment]``,
which gives the Fisher information matrix for each model vector
and each experiment vector.
.. note::
The default implementation of this method calls
:meth:`~DifferentiableModel.score()` for each possible outcome,
which can be quite slow. If possible, overriding this method can
give significant speed advantages.
"""
if self.is_n_outcomes_constant:
outcomes = np.arange(self.n_outcomes(expparams))
scores, L = self.score(outcomes, modelparams, expparams, return_L=True)
assert len(scores.shape) in (3, 4)
if len(scores.shape) == 3:
scores = scores[np.newaxis, :, :, :]
# Note that E[score] = 0 by regularity assumptions, so we only
# need the expectation over the outer product.
return np.einsum("ome,iome,jome->ijme",
L, scores, scores
)
else:
# Indexing will be a major pain here, so we need to start
# by making an empty array, so that index errors will be raised
# when (not if!) we make mistakes.
fisher = np.empty((
self.n_modelparams, self.n_modelparams,
modelparams.shape[0], expparams.shape[0]
))
# Now we loop over experiments, since we cannot vectorize the
# expectation value over data.
for idx_experiment, experiment in enumerate(expparams):
experiment = experiment.reshape((1,))
n_o = self.n_outcomes(experiment)
outcomes = np.arange(n_o)
scores, L = self.score(outcomes, modelparams, experiment, return_L=True)
fisher[:, :, :, idx_experiment] = np.einsum("ome,iome,jome->ijme",
L, scores, scores
)
return fisher