From 01271e38680a2a64f6f3a2c7c45300d86316c98c Mon Sep 17 00:00:00 2001
From: James Hensman <james.hensman@gmail.com>
Date: Fri, 21 Mar 2014 15:23:49 +0000
Subject: [PATCH] fixing the logexp (with MZ) and some stability issue in the
 stationary class

---
 GPy/core/parameterization/transformations.py | 12 ++++++------
 GPy/kern/_src/stationary.py                  |  7 +------
 2 files changed, 7 insertions(+), 12 deletions(-)

diff --git a/GPy/core/parameterization/transformations.py b/GPy/core/parameterization/transformations.py
index 5cda8d46..506d80cd 100644
--- a/GPy/core/parameterization/transformations.py
+++ b/GPy/core/parameterization/transformations.py
@@ -7,10 +7,10 @@ from domains import _POSITIVE,_NEGATIVE, _BOUNDED
 import weakref
 
 import sys
-#_lim_val = -np.log(sys.float_info.epsilon)
 
 _exp_lim_val = np.finfo(np.float64).max
-_lim_val = np.log(_exp_lim_val)
+_lim_val = 36.0
+epsilon = np.finfo(np.float64).resolution
 
 #===============================================================================
 # Fixing constants
@@ -54,19 +54,19 @@ class Transformation(object):
 class Logexp(Transformation):
     domain = _POSITIVE
     def f(self, x):
-        return np.where(x>_lim_val, x, np.log(1. + np.exp(np.clip(x, -_lim_val, _lim_val))))
+        return np.where(x>_lim_val, x, np.log(1. + np.exp(np.clip(x, -_lim_val, _lim_val)))) + epsilon
         #raises overflow warning: return np.where(x>_lim_val, x, np.log(1. + np.exp(x)))
     def finv(self, f):
         return np.where(f>_lim_val, f, np.log(np.exp(f+1e-20) - 1.))
     def gradfactor(self, f):
-        return np.where(f>_lim_val, 1., 1 - np.exp(-f))
+        return np.where(f>_lim_val, 1., 1. - np.exp(-f))
     def initialize(self, f):
         if np.any(f < 0.):
             print "Warning: changing parameters to satisfy constraints"
         return np.abs(f)
     def __str__(self):
         return '+ve'
-    
+
 
 class LogexpNeg(Transformation):
     domain = _POSITIVE
@@ -98,7 +98,7 @@ class NegativeLogexp(Transformation):
         return -self.logexp.initialize(f)  # np.abs(f)
     def __str__(self):
         return '-ve'
-    
+
 class LogexpClipped(Logexp):
     max_bound = 1e100
     min_bound = 1e-10
diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py
index a9e837a9..b6fea5ef 100644
--- a/GPy/kern/_src/stationary.py
+++ b/GPy/kern/_src/stationary.py
@@ -152,12 +152,7 @@ class Stationary(Kern):
         This term appears in derviatives.
         """
         dist = self._scaled_dist(X, X2).copy()
-        if X2 is None:
-            nondiag = util.diag.offdiag_view(dist)
-            nondiag[:] = 1./nondiag
-            return dist
-        else:
-            return 1./np.where(dist != 0., dist, np.inf)
+        return 1./np.where(dist != 0., dist, np.inf)
 
     def gradients_X(self, dL_dK, X, X2=None):
         """