Commented out weave functions for Py3 support

This commit is contained in:
Mike Croucher 2015-03-01 10:17:21 +00:00
parent 560950466d
commit a0dc90596c

View file

@ -2,10 +2,9 @@
# Licensed under the GNU GPL version 3.0 # Licensed under the GNU GPL version 3.0
import numpy as np import numpy as np
from scipy import weave #from scipy import weave
from . import linalg from . import linalg
def safe_root(N): def safe_root(N):
i = np.sqrt(N) i = np.sqrt(N)
j = int(i) j = int(i)
@ -13,58 +12,58 @@ def safe_root(N):
raise ValueError("N is not square!") raise ValueError("N is not square!")
return j return j
def flat_to_triang(flat): #def flat_to_triang(flat):
"""take a matrix N x D and return a M X M x D array where # """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.
# """
# 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
N = M(M+1)/2 #def triang_to_flat(L):
# M, _, D = L.shape
the lower triangluar portion of the d'th slice of the result is filled by the d'th column of flat. #
""" # L = np.ascontiguousarray(L) # should do nothing if L was created by flat_to_triang
N, D = flat.shape #
M = (-1 + safe_root(8*N+1))/2 # N = M*(M+1)/2
ret = np.zeros((M, M, D)) # flat = np.empty((N, D))
flat = np.ascontiguousarray(flat) # code = """
# int count = 0;
code = """ # for(int m=0; m<M; m++)
int count = 0; # {
for(int m=0; m<M; m++) # for(int mm=0; mm<=m; mm++)
{ # {
for(int mm=0; mm<=m; mm++) # for(int d=0; d<D; d++)
{ # {
for(int d=0; d<D; d++) # flat[count] = L[d + m*D*M + mm*D];
{ # count++;
ret[d + m*D*M + mm*D] = flat[count]; # }
count++; # }
} # }
} # """
} # weave.inline(code, ['flat', 'L', 'D', 'M'])
""" # return flat
weave.inline(code, ['flat', 'ret', 'D', 'M'])
return ret
def triang_to_flat(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_cov(L): def triang_to_cov(L):
return np.dstack([np.dot(L[:,:,i], L[:,:,i].T) for i in xrange(L.shape[-1])]) return np.dstack([np.dot(L[:,:,i], L[:,:,i].T) for i in xrange(L.shape[-1])])
@ -93,9 +92,6 @@ def multiple_dpotri_old(Ls):
def multiple_dpotri(Ls): def multiple_dpotri(Ls):
return np.dstack([linalg.dpotri(np.asfortranarray(Ls[:,:,i]), lower=1)[0] for i in range(Ls.shape[-1])]) 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): 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