manual merging with AS

GPy/util/block_matrices.py

@@ -1,9 +1,37 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Copyright (c) 2014-2015, Alan Saul
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 import numpy as np
 
+def get_blocks_3d(A, blocksizes, pagesizes=None):
+    """
+    Given a 3d matrix, make a block matrix, where the first and second dimensions are blocked according
+    to blocksizes, and the pages are blocked using pagesizes
+    """
+    assert (A.shape[0]==A.shape[1]) and len(A.shape)==3, "can't blockify this non-square matrix, may need to use 2d version"
+    N = np.sum(blocksizes)
+    assert A.shape[0] == N, "bad blocksizes"
+    num_blocks = len(blocksizes)
+    if pagesizes == None:
+        #Assume each page of A should be its own dimension
+        pagesizes = range(A.shape[2])#[0]*A.shape[2]
+    num_pages = len(pagesizes)
+    B = np.empty(shape=(num_blocks, num_blocks, num_pages), dtype=np.object)
+    count_k = 0
+    #for Bk, k in enumerate(pagesizes):
+    for Bk in pagesizes:
+        count_i = 0
+        for Bi, i in enumerate(blocksizes):
+            count_j = 0
+            for Bj, j in enumerate(blocksizes):
+                #We want to have it count_k:count_k + k but its annoying as it makes a NxNx1 array is page sizes are set to 1
+                B[Bi, Bj, Bk] = A[count_i:count_i + i, count_j:count_j + j, Bk]
+                count_j += j
+            count_i += i
+        #count_k += k
+    return B
+
 def get_blocks(A, blocksizes):
-    assert (A.shape[0]==A.shape[1]) and len(A.shape)==2, "can;t blockify this non-square matrix"
+    assert (A.shape[0]==A.shape[1]) and len(A.shape)==2, "can't blockify this non-square matrix"
     N = np.sum(blocksizes)
     assert A.shape[0] == N, "bad blocksizes"
     num_blocks = len(blocksizes)
@@ -17,6 +45,11 @@ def get_blocks(A, blocksizes):
         count_i += i
     return B
 
+def get_block_shapes_3d(B):
+    assert B.dtype is np.dtype('object'), "Must be a block matrix"
+    #FIXME: This isn't general AT ALL...
+    return get_block_shapes(B[:,:,0]), B.shape[2]
+
 def get_block_shapes(B):
     assert B.dtype is np.dtype('object'), "Must be a block matrix"
     return [B[b,b].shape[0] for b in range(0, B.shape[0])]
@@ -35,7 +68,7 @@ def unblock(B):
         count_i += i
     return A
 
-def block_dot(A, B):
+def block_dot(A, B, diagonal=False):
     """
     Element wise dot product on block matricies
 
@@ -48,21 +81,30 @@ def block_dot(A, B):
     +-------------+   +------+------+    +-------+-------+
 
     ..Note
+        If any block of either (A or B) are stored as 1d vectors then we assume
+        that it denotes a diagonal matrix efficient dot product using numpy
+        broadcasting will be used, i.e. A11*B11
+
         If either (A or B) of the diagonal matrices are stored as vectors then a more
         efficient dot product using numpy broadcasting will be used, i.e. A11*B11
     """
     #Must have same number of blocks and be a block matrix
     assert A.dtype is np.dtype('object'), "Must be a block matrix"
     assert B.dtype is np.dtype('object'), "Must be a block matrix"
-    Ashape = A.shape
-    Bshape = B.shape
-    assert Ashape == Bshape
-    def f(A,B):
-        if Ashape[0] == Ashape[1] or Bshape[0] == Bshape[1]:
-            #FIXME: Careful if one is transpose of other, would make a matrix
-            return A*B
+    assert A.shape == B.shape
+    def f(C,D):
+        """
+        C is an element of A, D is the associated element of B
+        """
+        Cshape = C.shape
+        Dshape = D.shape
+        if diagonal and (len(Cshape) == 1 or len(Dshape) == 1\
+                or C.shape[0] != C.shape[1] or D.shape[0] != D.shape[1]):
+            print "Broadcasting, C: {} D:{}".format(C.shape, D.shape)
+            return C*D
         else:
-            return np.dot(A,B)
+            print "Dotting, C: {} C:{}".format(C.shape, D.shape)
+            return np.dot(C,D)
     dot = np.vectorize(f, otypes = [np.object])
     return dot(A,B)
 

GPy/util/linalg.py

@@ -20,7 +20,7 @@ try:
     from scipy import weave
 except ImportError:
     config.set('weave', 'working', 'False')
-    
+
 
 _scipyversion = np.float64((scipy.__version__).split('.')[:2])
 _fix_dpotri_scipy_bug = True
@@ -102,7 +102,6 @@ def jitchol(A, maxtries=5):
         num_tries = 1
         while num_tries <= maxtries and np.isfinite(jitter):
             try:
-                print(jitter)
                 L = linalg.cholesky(A + np.eye(A.shape[0]) * jitter, lower=True)
                 return L
             except:
@@ -115,7 +114,6 @@ def jitchol(A, maxtries=5):
     except:
         logging.warning('\n'.join(['Added jitter of {:.10e}'.format(jitter),
             '  in '+traceback.format_list(traceback.extract_stack(limit=2)[-2:-1])[0][2:]]))
-    import ipdb;ipdb.set_trace()
     return L
 
 # def dtrtri(L, lower=1):

GPy/util/misc.py

@@ -2,18 +2,37 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 import numpy as np
+from scipy.special import cbrt
 from .config import *
 
 _lim_val = np.finfo(np.float64).max
-
 _lim_val_exp = np.log(_lim_val)
 _lim_val_square = np.sqrt(_lim_val)
-_lim_val_cube = np.power(_lim_val, -3)
+#_lim_val_cube = cbrt(_lim_val)
+_lim_val_cube = np.nextafter(_lim_val**(1/3.0), -np.inf)
+_lim_val_quad = np.nextafter(_lim_val**(1/4.0), -np.inf)
+_lim_val_three_times = np.nextafter(_lim_val/3.0, -np.inf)
 
 def safe_exp(f):
     clip_f = np.clip(f, -np.inf, _lim_val_exp)
     return np.exp(clip_f)
 
+def safe_square(f):
+    f = np.clip(f, -np.inf, _lim_val_square)
+    return f**2
+
+def safe_cube(f):
+    f = np.clip(f, -np.inf, _lim_val_cube)
+    return f**3
+
+def safe_quad(f):
+    f = np.clip(f, -np.inf, _lim_val_quad)
+    return f**4
+
+def safe_three_times(f):
+    f = np.clip(f, -np.inf, _lim_val_three_times)
+    return 3*f
+
 def chain_1(df_dg, dg_dx):
     """
     Generic chaining function for first derivative
@@ -39,8 +58,8 @@ def chain_2(d2f_dg2, dg_dx, df_dg, d2g_dx2):
         return d2f_dg2
     if  len(d2f_dg2) > 1 and len(d2f_dg2.shape)>1 and d2f_dg2.shape[-1] > 1:
         raise NotImplementedError('Not implemented for matricies yet')
-    #dg_dx_2 = np.clip(dg_dx, 1e-12, _lim_val_square)**2
-    dg_dx_2 = dg_dx**2
+    dg_dx_2 = np.clip(dg_dx, -np.inf, _lim_val_square)**2
+    #dg_dx_2 = dg_dx**2
     return d2f_dg2*(dg_dx_2) + df_dg*d2g_dx2
 
 def chain_3(d3f_dg3, dg_dx, d2f_dg2, d2g_dx2, df_dg, d3g_dx3):
@@ -55,8 +74,8 @@ def chain_3(d3f_dg3, dg_dx, d2f_dg2, d2g_dx2, df_dg, d3g_dx3):
     if (  (len(d2f_dg2) > 1 and d2f_dg2.shape[-1] > 1)
            or (len(d3f_dg3) > 1 and d3f_dg3.shape[-1] > 1)):
         raise NotImplementedError('Not implemented for matricies yet')
-    #dg_dx_3 = np.clip(dg_dx, 1e-12, _lim_val_cube)**3
-    dg_dx_3 = dg_dx**3
+    dg_dx_3 = np.clip(dg_dx, -np.inf, _lim_val_cube)**3
+    #dg_dx_3 = dg_dx**3
     return d3f_dg3*(dg_dx_3) + 3*d2f_dg2*dg_dx*d2g_dx2 + df_dg*d3g_dx3
 
 def opt_wrapper(m, **kwargs):