Source code for paramz.transformations

#===============================================================================
# Copyright (c) 2012 - 2014, GPy authors (see AUTHORS.txt).
# Copyright (c) 2015, Max Zwiessele
#
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
#
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * 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.
#
# * Neither the name of paramax 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.
#===============================================================================


import numpy as np
from .domains import _POSITIVE,_NEGATIVE, _BOUNDED
import weakref, logging

import sys

_log_lim_val = np.log(np.finfo(np.float64).max)
_exp_lim_val = np.finfo(np.float64).max
_lim_val = 36.0
epsilon = np.finfo(np.float64).resolution

#===============================================================================
# Fixing constants
__fixed__ = "fixed"
FIXED = False
UNFIXED = True
#===============================================================================
logger = logging.getLogger(__name__)

[docs]class Transformation(object):
    domain = None
    _instance = None
    def __new__(cls, *args, **kwargs):
        if not cls._instance or cls._instance.__class__ is not cls:
            cls._instance = super(Transformation, cls).__new__(cls, *args, **kwargs)
        return cls._instance
[docs]    def f(self, opt_param):
        raise NotImplementedError
[docs]    def finv(self, model_param):
        raise NotImplementedError
[docs]    def log_jacobian(self, model_param):
        """
        compute the log of the jacobian of f, evaluated at f(x)= model_param
        """
        logger.warning("log Jacobian for transformation {} not implemented, approximating by 0.".format(self.__class__))
        return 0.
[docs]    def log_jacobian_grad(self, model_param):
        """
        compute the drivative of the log of the jacobian of f, evaluated at f(x)= model_param
        """
        logger.warning("gradient of log Jacobian for transformation {} not implemented, approximating by 0.".format(self.__class__))
        return 0.
[docs]    def gradfactor(self, model_param, dL_dmodel_param):
        """ df(opt_param)_dopt_param evaluated at self.f(opt_param)=model_param, times the gradient dL_dmodel_param,

        i.e.:
        define

        .. math::

            \frac{\frac{\partial L}{\partial f}\left(\left.\partial f(x)}{\partial x}\right|_{x=f^{-1}(f)\right)}
        """
        raise NotImplementedError
[docs]    def gradfactor_non_natural(self, model_param, dL_dmodel_param):
        return self.gradfactor(model_param, dL_dmodel_param)
[docs]    def initialize(self, f):
        """ produce a sensible initial value for f(x)"""
        raise NotImplementedError
[docs]    def plot(self, xlabel=r'transformed $\theta$', ylabel=r'$\theta$', axes=None, *args,**kw):
        assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
        import matplotlib.pyplot as plt
        x = np.linspace(-8,8)
        plt.plot(x, self.f(x), *args, ax=axes, **kw)
        axes = plt.gca()
        axes.set_xlabel(xlabel)
        axes.set_ylabel(ylabel)
    def __str__(self):
        raise NotImplementedError
    def __repr__(self):
        return self.__class__.__name__

[docs]class Logexp(Transformation):
    domain = _POSITIVE
[docs]    def f(self, x):
        return np.where(x>_lim_val, x, np.log1p(np.exp(np.clip(x, -_log_lim_val, _lim_val)))) #+ epsilon
        #raises overflow warning: return np.where(x>_lim_val, x, np.log(1. + np.exp(x)))
[docs]    def finv(self, f):
        return np.where(f>_lim_val, f, np.log(np.expm1(f)))
[docs]    def gradfactor(self, f, df):
        return df*np.where(f>_lim_val, 1.,  - np.expm1(-f))
[docs]    def initialize(self, f):
        if np.any(f < 0.):
            logger.info("Warning: changing parameters to satisfy constraints")
        return np.abs(f)
[docs]    def log_jacobian(self, model_param):
        return np.where(model_param>_lim_val, model_param, np.log(np.expm1(model_param))) - model_param
[docs]    def log_jacobian_grad(self, model_param):
        return 1./(np.expm1(model_param))
    def __str__(self):
        return '+ve'

[docs]class Exponent(Transformation):
    domain = _POSITIVE
[docs]    def f(self, x):
        return np.where(x<_lim_val, np.where(x>-_lim_val, np.exp(x), np.exp(-_lim_val)), np.exp(_lim_val))
[docs]    def finv(self, x):
        return np.log(x)
[docs]    def gradfactor(self, f, df):
        return np.einsum('i,i->i', df, f)
[docs]    def initialize(self, f):
        if np.any(f < 0.):
            logger.info("Warning: changing parameters to satisfy constraints")
        return np.abs(f)
[docs]    def log_jacobian(self, model_param):
        return np.log(model_param)
[docs]    def log_jacobian_grad(self, model_param):
        return 1./model_param
    def __str__(self):
        return '+ve'

[docs]class NegativeLogexp(Transformation):
    domain = _NEGATIVE
    logexp = Logexp()
[docs]    def f(self, x):
        return -self.logexp.f(x)  # np.log(1. + np.exp(x))
[docs]    def finv(self, f):
        return self.logexp.finv(-f)  # np.log(np.exp(-f) - 1.)
[docs]    def gradfactor(self, f, df):
        return self.logexp.gradfactor(-f, -df)
[docs]    def initialize(self, f):
        return -self.logexp.initialize(-f)  # np.abs(f)
    def __str__(self):
        return '-ve'

[docs]class NegativeExponent(Exponent):
    domain = _NEGATIVE
[docs]    def f(self, x):
        return -Exponent.f(self, x)
[docs]    def finv(self, f):
        return -Exponent.finv(self, -f)
[docs]    def gradfactor(self, f, df):
        return -Exponent.gradfactor(self, f, df)
[docs]    def initialize(self, f):
        return -Exponent.initialize(self, f) #np.abs(f)
    def __str__(self):
        return '-ve'

[docs]class Square(Transformation):
    domain = _POSITIVE
[docs]    def f(self, x):
        return x ** 2
[docs]    def finv(self, x):
        return np.sqrt(x)
[docs]    def gradfactor(self, f, df):
        return np.einsum('i,i->i', df, 2 * np.sqrt(f))
[docs]    def initialize(self, f):
        return np.abs(f)
    def __str__(self):
        return '+sq'

[docs]class Logistic(Transformation):
    domain = _BOUNDED
    _instances = []
    def __new__(cls, lower=1e-6, upper=1e-6, *args, **kwargs):
        if cls._instances:
            cls._instances[:] = [instance for instance in cls._instances if instance()]
            for instance in cls._instances:
                if instance().lower == lower and instance().upper == upper:
                    return instance()
        newfunc = super(Transformation, cls).__new__
        if newfunc is object.__new__:
            o = newfunc(cls)
        else:
            o = newfunc(cls, lower, upper, *args, **kwargs)
        cls._instances.append(weakref.ref(o))
        return cls._instances[-1]()
    def __init__(self, lower, upper):
        assert lower < upper
        self.lower, self.upper = float(lower), float(upper)
        self.difference = self.upper - self.lower
[docs]    def f(self, x):
        if (x<-300.).any():
            x = x.copy()
            x[x<-300.] = -300.
        return self.lower + self.difference / (1. + np.exp(-x))
[docs]    def finv(self, f):
        return np.log(np.clip(f - self.lower, 1e-10, np.inf) / np.clip(self.upper - f, 1e-10, np.inf))
[docs]    def gradfactor(self, f, df):
        return np.einsum('i,i->i', df, (f - self.lower) * (self.upper - f) / self.difference)
[docs]    def initialize(self, f):
        if np.any(np.logical_or(f < self.lower, f > self.upper)):
            logger.info("Warning: changing parameters to satisfy constraints")
        #return np.where(np.logical_or(f < self.lower, f > self.upper), self.f(f * 0.), f)
        #FIXME: Max, zeros_like right?
        return np.where(np.logical_or(f < self.lower, f > self.upper), self.f(np.zeros_like(f)), f)
    def __str__(self):
        return '{},{}'.format(self.lower, self.upper)