From 5edf464efa69f1457a0e992250da0e5c54a10062 Mon Sep 17 00:00:00 2001
From: James Hensman <james.hensman@gmail.com>
Date: Mon, 27 Apr 2015 20:15:36 +0100
Subject: [PATCH] help choleskies along a little

---
 GPy/core/parameterization/index_operations.py |   3 +
 GPy/util/choleskies.py                        | 113 +++---------------
 2 files changed, 18 insertions(+), 98 deletions(-)

diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py
index e4803f37..fd6438ce 100644
--- a/GPy/core/parameterization/index_operations.py
+++ b/GPy/core/parameterization/index_operations.py
@@ -5,6 +5,8 @@ import numpy
 from numpy.lib.function_base import vectorize
 from .lists_and_dicts import IntArrayDict
 from functools import reduce
+from transformations import Transformation
+from priors import Prior
 
 def extract_properties_to_index(index, props):
     prop_index = dict()
@@ -144,6 +146,7 @@ class ParameterIndexOperations(object):
         return prop_index
 
     def add(self, prop, indices):
+        #assert isinstance(prop, (Transformation, Prior)), "invalid key"
         self._properties[prop] = combine_indices(self._properties[prop], indices)
 
     def remove(self, prop, indices):
diff --git a/GPy/util/choleskies.py b/GPy/util/choleskies.py
index 37ac7211..2accfedc 100644
--- a/GPy/util/choleskies.py
+++ b/GPy/util/choleskies.py
@@ -5,10 +5,7 @@ import numpy as np
 from . import linalg
 from .config import config
 
-try:
-    from scipy import weave
-except ImportError:
-    config.set('weave', 'working', 'False')
+import choleskies_cython
 
 def safe_root(N):
     i = np.sqrt(N)
@@ -17,36 +14,6 @@ def safe_root(N):
         raise ValueError("N is not square!")
     return j
 
-def _flat_to_triang_weave(flat):
-    """take a matrix N x D and return a M X M x D array where
-
-    N = M(M+1)/2
-
-    the lower triangluar portion of the d'th slice of the result is filled by the d'th column of flat.
-    This is the weave implementation
-    """
-    N, D = flat.shape
-    M = (-1 + safe_root(8*N+1))/2
-    ret = np.zeros((M, M, D))
-    flat = np.ascontiguousarray(flat)
-
-    code = """
-    int count = 0;
-    for(int m=0; m<M; m++)
-    {
-      for(int mm=0; mm<=m; mm++)
-      {
-        for(int d=0; d<D; d++)
-        {
-          ret[d + m*D*M + mm*D] = flat[count];
-          count++;
-        }
-      }
-    }
-    """
-    weave.inline(code, ['flat', 'ret', 'D', 'M'])
-    return ret
-
 def _flat_to_triang_pure(flat_mat):
     N, D = flat_mat.shape
     M = (-1 + safe_root(8*N+1))//2
@@ -59,34 +26,11 @@ def _flat_to_triang_pure(flat_mat):
               count = count+1
     return ret
 
-if config.getboolean('weave', 'working'):
-	flat_to_triang =  _flat_to_triang_weave
-else:
-        flat_to_triang =  _flat_to_triang_pure
+def _flat_to_triang_cython(flat_mat):
+    N, D = flat_mat.shape
+    M = (-1 + safe_root(8*N+1))//2
+    return choleskies_cython.flat_to_triang(flat_mat, M)
 
-def _triang_to_flat_weave(L):
-    M, _, D = L.shape
-
-    L = np.ascontiguousarray(L) # should do nothing if L was created by flat_to_triang
-
-    N = M*(M+1)/2
-    flat = np.empty((N, D))
-    code = """
-    int count = 0;
-    for(int m=0; m<M; m++)
-    {
-      for(int mm=0; mm<=m; mm++)
-      {
-        for(int d=0; d<D; d++)
-        {
-          flat[count] = L[d + m*D*M + mm*D];
-          count++;
-        }
-      }
-    }
-    """
-    weave.inline(code, ['flat', 'L', 'D', 'M'])
-    return flat
 
 def _triang_to_flat_pure(L):
     M, _, D = L.shape
@@ -101,41 +45,19 @@ def _triang_to_flat_pure(L):
                 count = count +1
     return flat
 
-if config.getboolean('weave', 'working'):
-    triang_to_flat =  _triang_to_flat_weave
-else:
-    triang_to_flat =  _triang_to_flat_pure
+def _triang_to_flat_cython(L):
+    return choleskies_cython.triang_to_flat(L)
 
 def triang_to_cov(L):
     return np.dstack([np.dot(L[:,:,i], L[:,:,i].T) for i in range(L.shape[-1])])
 
-def multiple_dpotri_old(Ls):
-    M, _, D = Ls.shape
-    Kis = np.rollaxis(Ls, -1).copy()
-    [dpotri(Kis[i,:,:], overwrite_c=1, lower=1) for i in range(D)]
-    code = """
-    for(int d=0; d<D; d++)
-    {
-      for(int m=0; m<M; m++)
-      {
-        for(int mm=0; mm<m; mm++)
-        {
-          Kis[d*M*M + mm*M + m ] = Kis[d*M*M + m*M + mm];
-        }
-      }
-    }
-
-    """
-    weave.inline(code, ['Kis', 'D', 'M'])
-    Kis = np.rollaxis(Kis, 0, 3) #wtf rollaxis?
-    return Kis
-
 def multiple_dpotri(Ls):
     return np.dstack([linalg.dpotri(np.asfortranarray(Ls[:,:,i]), lower=1)[0] for i in range(Ls.shape[-1])])
 
 def indexes_to_fix_for_low_rank(rank, size):
     """
-    work out which indexes of the flatteneed array should be fixed if we want the cholesky to represent a low rank matrix
+    Work out which indexes of the flatteneed array should be fixed if we want
+    the cholesky to represent a low rank matrix
     """
     #first we'll work out what to keep, and the do the set difference.
 
@@ -153,15 +75,10 @@ def indexes_to_fix_for_low_rank(rank, size):
     return np.setdiff1d(np.arange((size**2+size)/2), keep)
 
 
+if config.getboolean('cython', 'working'):
+    triang_to_flat = _triang_to_flat_cython
+    flat_to_triang = _flat_to_triang_cython
+else:
+    triang_to_flat =  _triang_to_flat_pure
+    flat_to_triang = _flat_to_triang_pure
 
-#class cholchecker(GPy.core.Model):
-    #def __init__(self, L, name='cholchecker'):
-        #super(cholchecker, self).__init__(name)
-        #self.L = GPy.core.Param('L',L)
-        #self.link_parameter(self.L)
-    #def parameters_changed(self):
-        #LL = flat_to_triang(self.L)
-        #Ki = multiple_dpotri(LL)
-        #self.L.gradient = 2*np.einsum('ijk,jlk->ilk', Ki, LL)
-        #self._loglik = np.sum([np.sum(np.log(np.abs(np.diag()))) for i in range(self.L.shape[-1])])
-#