TEST: Tests use 'lbfgsb' optimization function. Also some syntactic changes.

This commit is contained in:
Alexander Grigorievskiy 2016-03-15 17:01:05 +02:00
parent 7ecaf92ace
commit f0660dcde0
2 changed files with 1199 additions and 1159 deletions

View file

@ -5,16 +5,13 @@
Main functionality for state-space inference.
"""
import collections # for cheking whether a variable is iterable
import types # for cheking whether a variable is a function
import collections # for cheking whether a variable is iterable
import types # for cheking whether a variable is a function
import numpy as np
import scipy as sp
import scipy.linalg as linalg
#import pdb; pdb.set_trace()
try:
import state_space_setup
setup_available = True
@ -44,6 +41,7 @@ if print_verbose:
else:
print("state_space: cython is NOT used")
class Dynamic_Callables_Python(object):
def f_a(self, k, m, A):
@ -51,14 +49,16 @@ class Dynamic_Callables_Python(object):
p_a: function (k, x_{k-1}, A_{k}). Dynamic function.
k (iteration number), starts at 0
x_{k-1} State from the previous step
A_{k} Jacobian matrices of f_a. In the linear case it is exactly A_{k}.
A_{k} Jacobian matrices of f_a. In the linear case it is exactly
A_{k}.
"""
raise NotImplemented("f_a is not implemented!")
def Ak(self, k,m,P): # returns state iteration matrix
def Ak(self, k, m, P): # returns state iteration matrix
"""
function (k, m, P) return Jacobian of dynamic function, it is passed into p_a.
function (k, m, P) return Jacobian of dynamic function, it is passed
into p_a.
k (iteration number), starts at 0
m: point where Jacobian is evaluated
P: parameter for Jacobian, usually covariance matrix.
@ -75,8 +75,10 @@ class Dynamic_Callables_Python(object):
def Q_srk(self, k):
"""
function (k). Returns the square root of noise matrix of dynamic model on iteration k.
k (iteration number). starts at 0
function (k). Returns the square root of noise matrix of dynamic model
on iteration k.
k (iteration number). starts at 0
This function is implemented to use SVD prediction step.
"""
@ -97,9 +99,10 @@ class Dynamic_Callables_Python(object):
"""
raise NotImplemented("dQk is not implemented!")
def reset(self,compute_derivatives = False):
def reset(self, compute_derivatives=False):
"""
Return the state of this object to the beginning of iteration (to k eq. 0)
Return the state of this object to the beginning of iteration
(to k eq. 0).
"""
raise NotImplemented("reset is not implemented!")
@ -116,12 +119,13 @@ class Measurement_Callables_Python(object):
function (k, x_{k}, H_{k}). Measurement function.
k (iteration number), starts at 0
x_{k} state
H_{k} Jacobian matrices of f_h. In the linear case it is exactly H_{k}.
H_{k} Jacobian matrices of f_h. In the linear case it is exactly
H_{k}.
"""
raise NotImplemented("f_a is not implemented!")
def Hk(self, k,m_pred,P_pred): # returns state iteration matrix
def Hk(self, k, m_pred, P_pred): # returns state iteration matrix
"""
function (k, m, P) return Jacobian of measurement function, it is
passed into p_h.
@ -140,7 +144,6 @@ class Measurement_Callables_Python(object):
"""
raise NotImplemented("Rk is not implemented!")
def R_isrk(self, k):
"""
function (k). Returns the square root of the noise matrix of
@ -166,18 +169,21 @@ class Measurement_Callables_Python(object):
"""
raise NotImplemented("dQk is not implemented!")
def reset(self,compute_derivatives = False):
def reset(self, compute_derivatives=False):
"""
Return the state of this object to the beginning of iteration (to k eq. 0)
Return the state of this object to the beginning of iteration
(to k eq. 0)
"""
raise NotImplemented("reset is not implemented!")
if use_cython:
Measurement_Callables_Class = state_space_cython.Measurement_Callables_Cython
Measurement_Callables_Class = state_space_cython.\
Measurement_Callables_Cython
else:
Measurement_Callables_Class = Measurement_Callables_Python
class R_handling_Python(Measurement_Callables_Class):
"""
The calss handles noise matrix R.
@ -192,13 +198,15 @@ class R_handling_Python(Measurement_Callables_Class):
index - for each step of Kalman filter contains the corresponding index
in the array.
R_time_var_index - another index in the array R. Computed earlier and passed here.
R_time_var_index - another index in the array R. Computed earlier and
is passed here.
unique_R_number - number of unique noise matrices below which square roots
are cached and above which they are computed each time.
unique_R_number - number of unique noise matrices below which square
roots are cached and above which they are computed each time.
dR: 3D array[:, :, param_num]
derivative of R. Derivative is supported only when R do not change over time
derivative of R. Derivative is supported only when R do not change
over time
Output:
--------------
@ -208,7 +216,7 @@ class R_handling_Python(Measurement_Callables_Class):
"""
self.R = R
self.index = index
self.R_time_var_index = R_time_var_index
self.R_time_var_index = int(R_time_var_index)
self.dR = dR
if (len(np.unique(index)) > unique_R_number):
@ -218,38 +226,41 @@ class R_handling_Python(Measurement_Callables_Class):
self.R_square_root = {}
def Rk(self,k):
return self.R[:,:, self.index[self.R_time_var_index, k]]
def Rk(self, k):
return self.R[:, :, self.index[self.R_time_var_index, k]]
def dRk(self,k):
def dRk(self, k):
if self.dR is None:
raise ValueError("dR derivative is None")
return self.dR # the same dirivative on each iteration
return self.dR # the same dirivative on each iteration
def R_isrk(self, k):
"""
Function returns the inverse square root of R matrix on step k.
"""
ind = self.index[self.R_time_var_index, k]
R = self.R[:,:, ind ]
ind = int(self.index[self.R_time_var_index, k])
R = self.R[:, :, ind]
if (R.shape[0] == 1): # 1-D case handle simplier. No storage
# of the result, just compute it each time.
inv_square_root = np.sqrt( 1.0/R )
if (R.shape[0] == 1): # 1-D case handle simplier. No storage
# of the result, just compute it each time.
inv_square_root = np.sqrt(1.0/R)
else:
if self.svd_each_time:
(U,S,Vh) = sp.linalg.svd( R,full_matrices=False, compute_uv=True,
overwrite_a=False,check_finite=True)
(U, S, Vh) = sp.linalg.svd(R, full_matrices=False,
compute_uv=True, overwrite_a=False,
check_finite=True)
inv_square_root = U * 1.0/np.sqrt(S)
else:
if ind in self.R_square_root:
inv_square_root = self.R_square_root[ind]
else:
(U,S,Vh) = sp.linalg.svd( R,full_matrices=False, compute_uv=True,
overwrite_a=False,check_finite=True)
(U, S, Vh) = sp.linalg.svd(R, full_matrices=False,
compute_uv=True,
overwrite_a=False,
check_finite=True)
inv_square_root = U * 1.0/np.sqrt(S)
@ -262,25 +273,30 @@ if use_cython:
else:
R_handling_Class = R_handling_Python
class Std_Measurement_Callables_Python(R_handling_Class):
def __init__(self, H, H_time_var_index, R, index, R_time_var_index, unique_R_number, dH = None, dR=None):
super(Std_Measurement_Callables_Python,self).__init__(R, index, R_time_var_index, unique_R_number,dR)
def __init__(self, H, H_time_var_index, R, index, R_time_var_index,
unique_R_number, dH=None, dR=None):
super(Std_Measurement_Callables_Python,
self).__init__(R, index, R_time_var_index, unique_R_number, dR)
self.H = H
self.H_time_var_index = H_time_var_index
self.H_time_var_index = int(H_time_var_index)
self.dH = dH
def f_h(self, k, m, H):
"""
function (k, x_{k}, H_{k}). Measurement function.
k (iteration number), starts at 0
x_{k} state
H_{k} Jacobian matrices of f_h. In the linear case it is exactly H_{k}.
H_{k} Jacobian matrices of f_h. In the linear case it is exactly
H_{k}.
"""
return np.dot(H, m)
def Hk(self, k,m_pred,P_pred): # returns state iteration matrix
def Hk(self, k, m_pred, P_pred): # returns state iteration matrix
"""
function (k, m, P) return Jacobian of measurement function, it is
passed into p_h.
@ -289,22 +305,24 @@ class Std_Measurement_Callables_Python(R_handling_Class):
P: parameter for Jacobian, usually covariance matrix.
"""
return self.H[:,:, self.index[self.H_time_var_index, k]]
return self.H[:, :, self.index[self.H_time_var_index, k]]
def dHk(self,k):
def dHk(self, k):
if self.dH is None:
raise ValueError("dH derivative is None")
return self.dH # the same dirivative on each iteration
return self.dH # the same dirivative on each iteration
if use_cython:
Std_Measurement_Callables_Class = state_space_cython.Std_Measurement_Callables_Cython
Std_Measurement_Callables_Class = state_space_cython.\
Std_Measurement_Callables_Cython
else:
Std_Measurement_Callables_Class = Std_Measurement_Callables_Python
class Q_handling_Python(Dynamic_Callables_Class):
def __init__(self, Q, index, Q_time_var_index, unique_Q_number, dQ = None):
def __init__(self, Q, index, Q_time_var_index, unique_Q_number, dQ=None):
"""
Input:
---------------
@ -314,13 +332,15 @@ class Q_handling_Python(Dynamic_Callables_Class):
index - for each step of Kalman filter contains the corresponding index
in the array.
R_time_var_index - another index in the array R. Computed earlier and passed here.
R_time_var_index - another index in the array R. Computed earlier and
passed here.
unique_R_number - number of unique noise matrices below which square roots
are cached and above which they are computed each time.
unique_R_number - number of unique noise matrices below which square
roots are cached and above which they are computed each time.
dQ: 3D array[:, :, param_num]
derivative of Q. Derivative is supported only when Q do not change over time
derivative of Q. Derivative is supported only when Q do not
change over time
Output:
--------------
@ -341,39 +361,40 @@ class Q_handling_Python(Dynamic_Callables_Class):
self.Q_square_root = {}
def Qk(self,k):
def Qk(self, k):
"""
function (k). Returns noise matrix of dynamic model on iteration k.
k (iteration number). starts at 0
"""
return self.Q[:,:, self.index[self.Q_time_var_index, k]]
return self.Q[:, :, self.index[self.Q_time_var_index, k]]
def dQk(self,k):
def dQk(self, k):
if self.dQ is None:
raise ValueError("dQ derivative is None")
return self.dQ # the same dirivative on each iteration
return self.dQ # the same dirivative on each iteration
def Q_srk(self,k):
def Q_srk(self, k):
"""
function (k). Returns the square root of noise matrix of dynamic model on iteration k.
k (iteration number). starts at 0
function (k). Returns the square root of noise matrix of dynamic model
on iteration k.
k (iteration number). starts at 0
This function is implemented to use SVD prediction step.
"""
ind = self.index[self.Q_time_var_index, k]
Q = self.Q[:,:, ind]
Q = self.Q[:, :, ind]
if (Q.shape[0] == 1): # 1-D case handle simplier. No storage
# of the result, just compute it each time.
square_root = np.sqrt( Q )
if (Q.shape[0] == 1): # 1-D case handle simplier. No storage
# of the result, just compute it each time.
square_root = np.sqrt(Q)
else:
if self.svd_each_time:
(U,S,Vh) = sp.linalg.svd( Q,full_matrices=False, compute_uv=True,
overwrite_a=False,check_finite=True)
(U, S, Vh) = sp.linalg.svd(Q, full_matrices=False,
compute_uv=True,
overwrite_a=False,
check_finite=True)
square_root = U * np.sqrt(S)
else:
@ -381,8 +402,10 @@ class Q_handling_Python(Dynamic_Callables_Class):
if ind in self.Q_square_root:
square_root = self.Q_square_root[ind]
else:
(U,S,Vh) = sp.linalg.svd( Q,full_matrices=False, compute_uv=True,
overwrite_a=False,check_finite=True)
(U, S, Vh) = sp.linalg.svd(Q, full_matrices=False,
compute_uv=True,
overwrite_a=False,
check_finite=True)
square_root = U * np.sqrt(S)
@ -395,10 +418,13 @@ if use_cython:
else:
Q_handling_Class = Q_handling_Python
class Std_Dynamic_Callables_Python(Q_handling_Class):
def __init__(self, A, A_time_var_index, Q, index, Q_time_var_index, unique_Q_number, dA = None, dQ=None):
super(Std_Dynamic_Callables_Python,self).__init__(Q, index, Q_time_var_index, unique_Q_number,dQ)
def __init__(self, A, A_time_var_index, Q, index, Q_time_var_index,
unique_Q_number, dA=None, dQ=None):
super(Std_Dynamic_Callables_Python,
self).__init__(Q, index, Q_time_var_index, unique_Q_number, dQ)
self.A = A
self.A_time_var_index = A_time_var_index
@ -409,12 +435,12 @@ class Std_Dynamic_Callables_Python(Q_handling_Class):
f_a: function (k, x_{k-1}, A_{k}). Dynamic function.
k (iteration number), starts at 0
x_{k-1} State from the previous step
A_{k} Jacobian matrices of f_a. In the linear case it is exactly A_{k}.
A_{k} Jacobian matrices of f_a. In the linear case it is exactly
A_{k}.
"""
return np.dot(A, m)
return np.dot(A,m)
def Ak(self, k,m_pred,P_pred): # returns state iteration matrix
def Ak(self, k, m_pred, P_pred): # returns state iteration matrix
"""
function (k, m, P) return Jacobian of measurement function, it is
passed into p_h.
@ -423,23 +449,25 @@ class Std_Dynamic_Callables_Python(Q_handling_Class):
P: parameter for Jacobian, usually covariance matrix.
"""
return self.A[:,:, self.index[self.A_time_var_index, k]]
return self.A[:, :, self.index[self.A_time_var_index, k]]
def dAk(self,k):
def dAk(self, k):
if self.dA is None:
raise ValueError("dA derivative is None")
return self.dA # the same dirivative on each iteration
return self.dA # the same dirivative on each iteration
def reset(self,compute_derivatives = False):
def reset(self, compute_derivatives=False):
"""
Return the state of this object to the beginning of iteration (to k eq. 0)
Return the state of this object to the beginning of iteration
(to k eq. 0)
"""
return self
if use_cython:
Std_Dynamic_Callables_Class = state_space_cython.Std_Dynamic_Callables_Cython
Std_Dynamic_Callables_Class = state_space_cython.\
Std_Dynamic_Callables_Cython
else:
Std_Dynamic_Callables_Class = Std_Dynamic_Callables_Python
@ -460,19 +488,20 @@ class AddMethodToClass(object):
self.func = func
self.tp = tp
def __get__(self, obj, klass=None,*args, **kwargs):
def __get__(self, obj, klass=None, *args, **kwargs):
if self.tp == 'staticmethod':
return self.func
elif self.tp == 'normal':
def newfunc(obj, *args, **kwargs):
return self.func
return self.func
elif self.tp == 'classmethod':
def newfunc(klass, *args, **kwargs):
return self.func
return self.func
return newfunc
class DescreteStateSpaceMeta(type):
"""
Substitute necessary methods from cython.
@ -485,15 +514,23 @@ class DescreteStateSpaceMeta(type):
if use_cython:
if '_kalman_prediction_step_SVD' in attributes:
attributes['_kalman_prediction_step_SVD'] = AddMethodToClass(state_space_cython._kalman_prediction_step_SVD_Cython)
attributes['_kalman_prediction_step_SVD'] =\
AddMethodToClass(state_space_cython.
_kalman_prediction_step_SVD_Cython)
if '_kalman_update_step_SVD' in attributes:
attributes['_kalman_update_step_SVD'] = AddMethodToClass(state_space_cython._kalman_update_step_SVD_Cython)
attributes['_kalman_update_step_SVD'] =\
AddMethodToClass(state_space_cython.
_kalman_update_step_SVD_Cython)
if '_cont_discr_kalman_filter_raw' in attributes:
attributes['_cont_discr_kalman_filter_raw'] = AddMethodToClass(state_space_cython._cont_discr_kalman_filter_raw_Cython)
attributes['_cont_discr_kalman_filter_raw'] =\
AddMethodToClass(state_space_cython.
_cont_discr_kalman_filter_raw_Cython)
return super(DescreteStateSpaceMeta,
typeclass).__new__(typeclass, name, bases, attributes)
return super(DescreteStateSpaceMeta, typeclass).__new__(typeclass, name, bases, attributes)
class DescreteStateSpace(object):
"""
@ -539,32 +576,37 @@ class DescreteStateSpace(object):
along dimension 0, sample dimension - dimension 1, different
time series - dimension 2.
old_shape: tuple or None
If the shape has been modified, return old shape, otherwise None
If the shape has been modified, return old shape, otherwise
None.
"""
if (len(shape) > 3):
raise ValueError("Input array is not supposed to be more than 3 dimensional.")
raise ValueError("""Input array is not supposed to be more
than 3 dimensional.""")
if (len(shape) > desired_dim):
raise ValueError("Input array shape is more than desired shape.")
elif len(shape) == 1:
if (desired_dim==3):
return ((shape[0],1,1), shape) # last dimension is the time serime_series_no
elif (desired_dim==2):
return ((shape[0],1), shape)
if (desired_dim == 3):
return ((shape[0], 1, 1), shape) # last dimension is the
# time serime_series_no
elif (desired_dim == 2):
return ((shape[0], 1), shape)
elif len(shape) == 2:
if (desired_dim==3):
return ((shape[1],1,1), shape) if (shape[0] == 1) else ( (shape[0],shape[1],1), shape) # convert to column vector
elif (desired_dim==2):
return ((shape[1],1), shape) if (shape[0] == 1) else ( (shape[0],shape[1]), None) # convert to column vector
if (desired_dim == 3):
return ((shape[1], 1, 1), shape) if (shape[0] == 1) else\
((shape[0], shape[1], 1), shape) # convert to column
# vector
elif (desired_dim == 2):
return ((shape[1], 1), shape) if (shape[0] == 1) else\
((shape[0], shape[1]), None) # convert to column vector
else: # len(shape) == 3
return (shape,None) # do nothing
else: # len(shape) == 3
return (shape, None) # do nothing
@classmethod
def kalman_filter(cls,p_A, p_Q, p_H, p_R, Y, index = None, m_init=None,
def kalman_filter(cls, p_A, p_Q, p_H, p_R, Y, index=None, m_init=None,
P_init=None, p_kalman_filter_type='regular',
calc_log_likelihood=False,
calc_grad_log_likelihood=False, grad_params_no=None,

View file

@ -31,10 +31,8 @@ class StateSpaceKernelsTests(np.testing.TestCase):
if check_gradients:
self.assertTrue(m1.checkgrad())
#import pdb; pdb.set_trace()
if optimize:
m1.optimize(optimizer='bfgs',max_iters=optimize_max_iters)
m1.optimize(optimizer='lbfgsb',max_iters=optimize_max_iters)
if compare_with_GP and (predict_X is None):
predict_X = X
@ -45,7 +43,7 @@ class StateSpaceKernelsTests(np.testing.TestCase):
if compare_with_GP:
m2 = GPy.models.GPRegression(X,Y, gp_kernel)
m2.optimize(optimizer='bfgs', max_iters=optimize_max_iters)
m2.optimize(optimizer='lbfgsb', max_iters=optimize_max_iters)
#print(m2)
x_pred_reg_2 = m2.predict(predict_X)