def compute_b_memory(X, omega, u):
    assert len(X.shape) == 2
    Phi2 = rff_feature_map_grad2(X, omega, u)
    return -np.mean(np.sum(Phi2, 0), 0) def compute_C_memory(X, omega, u):
    assert len(X.shape) == 2
    Phi2 = rff_feature_map_grad(X, omega, u)
    d = X.shape[1]
    N = X.shape[0]
    m = Phi2.shape[2]
    
#     # bottleneck! loop version is very slow
#     C = np.zeros((m, m))
#     t = time.time()
#     for i in range(N):
#         for ell in range(d):
#             phi2 = Phi2[ell, i]
#             C += np.outer(phi2, phi2)
#     print("loop", time.time()-t)
#      
#     # roughly 5x faster than the above loop
#     t = time.time()
#     Phi2_reshaped = Phi2.reshape(N*d, m)
#     C2=np.einsum('ij,ik->jk', Phi2_reshaped, Phi2_reshaped)
#     print("einsum", time.time()-t)
#
#     # cython implementation, is slowest
#     t = time.time()
#     Phi2_reshaped = Phi2.reshape(N*d, m)
#     C3 = outer_sum_cython(Phi2_reshaped)
#     print("cython", time.time()-t)
    
#     t = time.time()

    # fastest version using multicore: tensordot method
    Phi2_reshaped = Phi2.reshape(N * d, m)
    C4 = np.tensordot(Phi2_reshaped, Phi2_reshaped, [0, 0])
#     print("tensordot", time.time()-t)

    return C4 / N def _objective_sym_completely_manual(X, theta, omega, u):
    N = X.shape[0]
    D = X.shape[1]
    m = len(theta)
    
    J_manual = 0.
     
    for n in range(N):
        for d in range(D):
            b_term_manual = -np.sqrt(2. / m) * np.cos(np.dot(X[n], omega) + u) * (omega[d, :] ** 2)
            J_manual -= -np.dot(b_term_manual, theta)
             
            c_vec_manual = -np.sqrt(2. / m) * np.sin(np.dot(X[n], omega) + u) * omega[d, :]
            C_term_manual = np.outer(c_vec_manual, c_vec_manual)
            J_manual += 0.5 * np.dot(theta, np.dot(C_term_manual, theta))
    
    J_manual /= N
    return J_manual def _objective_sym_half_manual(X, theta, omega, u):
    N = X.shape[0]
    D = X.shape[1]
    
    J_manual = 0.
     
    for n in range(N):
        for d in range(D):
            b_term = -rff_feature_map_grad2_d(X[n], omega, u, d)
            J_manual -= np.dot(b_term, theta)

            c_vec = rff_feature_map_grad_d(X[n], omega, u, d)
            C_term_manual = np.outer(c_vec, c_vec)
            J_manual += 0.5 * np.dot(theta, np.dot(C_term_manual, theta))
    
    J_manual /= N
    return J_manual def update_b_single(x, b, n, omega, u):
    D = omega.shape[0]
    m = omega.shape[1]
    
    projections_sum = np.zeros(m)
    phi = rff_feature_map_single(x, omega, u)
    for d in range(D):
        # second derivative of feature map
        phi2 = -phi * (omega[d, :] ** 2)
        projections_sum -= phi2
    
    # Knuth's running average
    n += 1
    delta = projections_sum - b
    b += delta / n
    
    return b def update_L_C_single(x, L_C, n, omega, u):
    D = omega.shape[0]
    assert x.ndim == 1
    assert len(x) == D
    m = 1 if np.isscalar(u) else len(u)
    N = 1
    
    L_C *= np.sqrt(n)
    
    # since cholupdate works on transposed version
    L_C = L_C.T
    
    projection = np.dot(x[np.newaxis, :], omega) + u
    np.sin(projection, projection)
    projection *= -np.sqrt(2. / m)
    temp = np.zeros((N, m))
    for d in range(D):
        temp = -projection * omega[d, :]
        
        # here temp is 1xm, this costs O(m^2)
        cholupdate(L_C, temp[0])
        
    # since cholupdate works on transposed version
    L_C = L_C.T
    
    # since the new C has a 1/(n+1) term in it
    L_C /= np.sqrt(n + 1)
    
    return L_C def fit(X, omega, u, b=None, C=None):
    if b is None:
        b = compute_b(X, omega, u)
    
    if C is None:
        C = compute_C(X, omega, u)
    
    theta = np.linalg.solve(C, b)
    return theta def objective(X, theta, omega, u, b=None, C=None):
    if b is None:
        b = compute_b(X, omega, u)
        
    if C is None:
        C = compute_C(X, omega, u)
    
    return 0.5 * np.dot(theta, np.dot(C, theta)) - np.dot(theta, b) def compute_b(X, omega, u):
    assert len(X.shape) == 2
    m = 1 if np.isscalar(u) else len(u)
    D = X.shape[1]
    
    projections_sum = np.zeros(m)
    Phi2 = rff_feature_map(X, omega, u)
    for d in range(D):
        projections_sum += np.mean(-Phi2 * (omega[d, :] ** 2), axis=0)
        
    return -projections_sum def compute_C(X, omega, u):
    assert len(X.shape) == 2
    m = 1 if np.isscalar(u) else len(u)
    N = X.shape[0]
    D = X.shape[1]
    
    C = np.zeros((m, m))
    projection = np.dot(X, omega) + u
    np.sin(projection, projection)
    projection *= -np.sqrt(2. / m)
    temp = np.zeros((N, m))
    for d in range(D):
        temp = -projection * omega[d, :]
        C += np.tensordot(temp, temp, [0, 0])

    return C / N def update_C(x, C, n, omega, u):
    D = omega.shape[0]
    assert x.ndim == 1
    assert len(x) == D
    m = 1 if np.isscalar(u) else len(u)
    N = 1
    
    C_new = np.zeros((m, m))
    projection = np.dot(x[np.newaxis, :], omega) + u
    np.sin(projection, projection)
    projection *= -np.sqrt(2. / m)
    temp = np.zeros((N, m))
    for d in range(D):
        temp = -projection * omega[d, :]
        C_new += np.tensordot(temp, temp, [0, 0])
    
    # Knuth's running average
    n = n + 1
    delta = C_new - C
    C += delta / n
    
    return C class KernelExpFiniteGaussian(EstimatorBase):
    def __init__(self, sigma, lmbda, m, D):
        self.sigma = sigma
        self.lmbda = lmbda
        self.m = m
        self.D = D
        self.omega, self.u = rff_sample_basis(D, m, sigma)
        
        self._initialise_solution()
    
    def _initialise_solution(self):
        # components of linear system, stored for online updating
        # set in first iteration of update_fit, or in fit
        self.L_C = None
        self.b = None
        self.log_sum_weights = None
        
        # number of observed points via fit or update_fit
        self.n = 0
        
        # initial solution is just a flat function
        self.theta = np.zeros(self.m)
    
    def supports_update_fit(self):
        return True
    
    def supports_weights(self):
        return True
    
    def fit(self, X, log_weights=None):
        assert_array_shape(X, ndim=2, dims={1: self.D})
        N = len(X)

        # in any case, delete previous solution
        self._initialise_solution()
        
        if N <= 0:
            # dont do anything if no data observed
            return
        elif N == 1:
            # can get away with single update as not expensive
            self.update_fit(X, log_weights)
            return
        
        if log_weights is None:
            # b is the same as first x is used as b straight away in update_fit
            self.b = compute_b(X, self.omega, self.u)
            
            # remove first term from C computation as it is replaced with regulariser
            C = compute_C(X[1:], self.omega, self.u)
            
            # C so far consists of the average of N-1 terms
            # additional term is regulariser C (first in update_fit)
            # use Knuth online-update for new mean, which is average of N terms
            # effectively, this is equal to
            # C = (C * (N - 1) + np.eye(self.m) * self.lmbda) / N
            delta = np.eye(self.m) * self.lmbda - C
            C += delta / N
            self.L_C = np.linalg.cholesky(C)
            
            # as all weights are equal, this corresponds to repeated update_fit calls
            self.log_sum_weights = np.log(N)
            
            self.n = N
            
            self.theta = fit_L_C_precomputed(self.b, self.L_C)
        else:
            # weighted batch learning here corresponds to repeated online-learning
            assert_array_shape(log_weights, ndim=1, dims={0: N})
            self.update_fit(X, log_weights)
    
    def update_fit(self, X, log_weights=None):
        assert_array_shape(X, ndim=2, dims={1: self.D})
        N = len(X)
        
        # dont do anything if no data observed
        if N == 0:
            return
        
        if log_weights is None:
            log_weights = np.log(np.ones(N))
        assert_array_shape(log_weights, ndim=1, dims={0: N})
        
        # first update: use first of X and log_weights, and then discard
        if self.log_sum_weights is None:
            # assume have observed fake terms, which is needed for making the system well-posed
            # the L_C says that the fake terms had covariance self.lmbda, which is a regulariser
            self.L_C = np.eye(self.m) * np.sqrt(self.lmbda)
            self.log_sum_weights = log_weights[0]
            self.b = compute_b(X[0].reshape(1, self.D), self.omega, self.u)
            self.n = 1
            
            X = X[1:]
            log_weights = log_weights[1:]
            N -= 1
            
        # dont do anything if no data observed
        if N == 0:
            return
        
        old_L_C = np.array(self.L_C, copy=True)
        self.b, self.L_C = update_b_L_C_weighted(X, self.b, self.L_C,
                                                 self.log_sum_weights,
                                                 log_weights,
                                                 self.omega, self.u)
        
        if np.any(np.isnan(self.L_C)) or np.any(np.isinf(self.L_C)):
            logger.warning("Numerical error while updating Cholesky factor of C.\n"
                           "Before update:\n%s\n"
                           "After update:\n%s\n"
                           "Updating data:\n%s\n"
                           "Updating log weights:\n%s\n"
                           % (str(old_L_C), str(self.L_C), str(X), str(log_weights))
                           )
            raise RuntimeError("Numerical error while updating Cholesky factor of C.")
        
        # update terms and weights
        self.n += len(X)
        self.log_sum_weights = log_sum_exp(list(log_weights) + [self.log_sum_weights])
        
        # finally update solution
        self.theta = fit_L_C_precomputed(self.b, self.L_C)
    
    def log_pdf(self, x):
        if self.theta is None:
            raise RuntimeError("Model not fitted yet.")
        assert_array_shape(x, ndim=1, dims={0: self.D})
        
        phi = rff_feature_map_single(x, self.omega, self.u)
        return np.dot(phi, self.theta)
    
    def grad(self, x):
        if self.theta is None:
            raise RuntimeError("Model not fitted yet.")
        
        grad = rff_feature_map_grad_single(x, self.omega, self.u)
        return np.dot(grad, self.theta)
    
    if theano_available:
        def hessian(self, x):
            """
            Computes the Hessian of the learned log-density function.
            
            WARNING: This implementation slow, so don't call repeatedly.
            """
            assert_array_shape(x, ndim=1, dims={0: self.D})
            
            H = np.zeros((self.D, self.D))
            for i, theta_i in enumerate(self.theta):
                H += theta_i * rff_feature_map_comp_hessian_theano(x, self.omega[:, i], self.u[i])
        
            # RFF is a monte carlo average, so have to normalise by np.sqrt(m) here
            return H / np.sqrt(self.m)
        
        def third_order_derivative_tensor(self, x):
            """
            Computes the third order derivative tensor of the learned log-density function.
            
            WARNING: This implementation is slow, so don't call repeatedly.
            """
            assert_array_shape(x, ndim=1, dims={0: self.D})
            
            G3 = np.zeros((self.D, self.D, self.D))
            for i, theta_i in enumerate(self.theta):
                G3 += theta_i * rff_feature_map_comp_third_order_tensor_theano(x, self.omega[:, i], self.u[i])
        
            # RFF is a monte carlo average, so have to normalise by np.sqrt(m) here
            return G3 / np.sqrt(self.m)
    
    def log_pdf_multiple(self, X):
        assert_array_shape(X, ndim=2, dims={1: self.D})
        
        Phi = rff_feature_map(X, self.omega, self.u)
        return np.dot(Phi, self.theta)
    
    def objective(self, X):
        assert_array_shape(X, ndim=2, dims={1: self.D})
        
        # note we need to recompute b and C here
        return objective_L_C_precomputed(X, self.theta, self.omega, self.u, self.b, self.L_C)

    def get_parameter_names(self):
        return ['sigma', 'lmbda']
    
    def set_parameters_from_dict(self, param_dict):
        EstimatorBase.set_parameters_from_dict(self, param_dict)
        
        # update basis
        self.omega, self.u = rff_sample_basis(self.D, self.m, self.sigma) def update_b_L_C_weighted(X, b, L_C, log_sum_weights, log_weights, omega, u):
    m = 1 if np.isscalar(u) else len(u)
    N = X.shape[0]
    D = X.shape[1]
    
    # transform weights into (1-\lmbda)*old_mean+ \lmbda*new_term style updates
    lmbdas = log_weights_to_lmbdas(log_sum_weights, log_weights)
    
    # first and (negative) second derivatives of rff feature map
    projection = np.dot(X, omega) + u
    Phi = np.cos(projection) * np.sqrt(2. / m)
    Phi2 = np.sin(projection) * np.sqrt(2. / m)
    
    # not needed any longer
    del projection
    
    # work on upper triangular cholesky internally
    L_R = L_C.T
    
    b_new_term = np.zeros(m)
    for i in range(N):
        # downscale L_C once for every datum
        L_R *= np.sqrt(1 - lmbdas[i])
        
        b_new_term[:] = 0
        for d in range(D):
            b_new_term += Phi[i] * (omega[d, :] ** 2)
            
            # L_C is updated D times for every datum, each with fixed lmbda
            C_new_term = Phi2[i] * omega[d, :] * np.sqrt(lmbdas[i])
            cholupdate(L_R, C_new_term)
        
        # b is updated once per datum
        b = (1 - lmbdas[i]) * b + lmbdas[i] * b_new_term
    
    # transform back to lower triangular version
    L_C = L_R.T
    
    return b, L_C     def get_expr_gaussian_kernel(x, y, sigma):
    def get_expr_gaussian_kernel_grad(x, y, sigma):
    def get_expr_gaussian_kernel_hessian(x, y, sigma):
    def get_expr_gaussian_kernel_third_order_tensor(x, y, sigma):
    def get_expr_rff_feature_map_component(x, omega, u):
    def get_expr_rff_feature_map_component_grad(x, omega, u):
    def get_expr_rff_feature_map_component_hessian(x, omega, u):
    def get_expr_rff_feature_map_component_third_order_tensor(x, omega, u):
def gaussian_kernel(X, Y=None, sigma=1.):
def gaussian_kernel_grad(x, Y, sigma=1.):
def gaussian_kernel_dx_component(x, y, ell, sigma=1.):
def gaussian_kernel_hessian(x, y, sigma=1.):
def gaussian_kernel_hessian_entry(x, y, i, j, sigma=1.):
def gaussian_kernel_hessians(X, Y=None, sigma=1.0):
def gaussian_kernel_dx_dx(x, Y, sigma=1.):
def gaussian_kernel_dx_dx_component(x, y, ell, sigma=1.):
def gaussian_kernel_dx_dx_dy(x, y, sigma=1.):
def gaussian_kernel_dx_dx_dy_dy(x, y, sigma=1.):
def gaussian_kernel_dx_i_dx_j(x, y, sigma=1.):
def gaussian_kernel_dx_i_dx_j_component(x, y, ell, sigma=1.):
def gaussian_kernel_dx_i_dx_i_dx_j(x, y, sigma=1.):
def gaussian_kernel_dx_i_dx_i_dx_j_component(x, y, ell, sigma=1.):
def gaussian_kernel_dx_i_dx_i_dx_j_dx_j(x,y, sigma=1.):
def rff_sample_basis(D, m, sigma):
def rff_feature_map_single(x, omega, u):
def rff_feature_map(X, omega, u):
def rff_feature_map_grad_d(X, omega, u, d):
def rff_feature_map_grad2_d(X, omega, u, d):
def rff_feature_map_grad(X, omega, u):
def rff_feature_map_grad2(X, omega, u):
def rff_feature_map_grad_single(x, omega, u):
    K = np.exp(-(sq_dists) / sigma)
    G = (2.0 / sigma) * (k.T * differences)
    G = (2.0 / sigma) * (k * difference)
    H = k*(2*np.eye(d)/sigma - 4*np.outer(differences, differences)/sigma**2)
    H = k*(ridge - 4*(differences_i*differences_j)/sigma**2)
    N = X.shape[0]
    D = X.shape[1]
    N = X.shape[0]
    D = X.shape[1] def log_sum_exp(X):
    """
    Computes log sum_i exp(X_i).
    Useful if you want to solve log \int f(x)p(x) dx
    where you have samples from p(x) and can compute log f(x)
    """
    # extract minimum
    X0 = np.min(X)
    X_without_X0 = np.delete(X, np.argmin(X))
    
    return X0 + np.log(1 + np.sum(np.exp(X_without_X0 - X0)))
from nose import SkipTest
from nose.tools import assert_less_equal, assert_almost_equal, assert_equal
from numpy.ma.testutils import assert_close
from numpy.testing.utils import assert_allclose
from kernel_exp_family.estimators.finite.develop.gaussian import compute_b_memory, \
    compute_C_memory, _objective_sym_completely_manual, \
    _objective_sym_half_manual, update_b_single, update_L_C_single
from kernel_exp_family.estimators.finite.gaussian import fit, objective, \
    compute_b, compute_C, update_C, KernelExpFiniteGaussian, \
    update_b_L_C_weighted
from kernel_exp_family.kernels.kernels import rff_feature_map_grad2_d, \
    rff_feature_map_grad_d, theano_available, rff_sample_basis
from kernel_exp_family.tools.numerics import log_sum_exp
import numpy as np
    m = 10
    omega = np.random.randn(D, m)
    u = np.random.uniform(0, 2 * np.pi, m)
    X = np.random.randn(N, D)
    
    b = compute_b(X, omega, u)
    b_storage = compute_b_memory(X, omega, u)
    assert_allclose(b, b_storage)

def test_compute_C_equals_compute_C_memory():
    N = 100
    D = 3
    m = 10
    omega = np.random.randn(D, m)
    u = np.random.uniform(0, 2 * np.pi, m)
    X = np.random.randn(N, D)
    
    C = compute_C(X, omega, u)
    C_storage = compute_C_memory(X, omega, u)
    assert_allclose(C, C_storage, rtol=1e-4)

def test_update_b_single_equals_batch():
    N = 100
    D = 3
    m = 10
    omega = np.random.randn(D, m)
    u = np.random.uniform(0, 2 * np.pi, m)
    X = np.random.randn(N, D)
    x = np.random.randn(D)
    
    b = compute_b(X, omega, u)
    b = update_b_single(x, b, n=N, omega=omega, u=u)
    b_batch = compute_b(np.vstack((X, x)), omega, u)
    
    assert_allclose(b, b_batch)

def test_update_C_equals_batch():
    N = 100
    D = 3
    m = 10
    omega = np.random.randn(D, m)
    u = np.random.uniform(0, 2 * np.pi, m)
    X = np.random.randn(N, D)
    x = np.random.randn(D)
    
    C = compute_C(X, omega, u)
    C = update_C(x, C, n=N, omega=omega, u=u)
    C_batch = compute_C(np.vstack((X, x)), omega, u)
    
    assert_allclose(C, C_batch)

def test_update_b_L_C_weighted_equals_compute_b_and_compute_L_C_constant_weights():
    N = 2
    D = 2
    m = 2
    omega = np.random.randn(D, m)
    u = np.random.uniform(0, 2 * np.pi, m)
    X1 = np.random.randn(N, D)
    X2 = np.random.randn(N, D)
    log_weights1 = np.log(np.ones(N))
    log_weights2 = np.log(np.ones(N))
    
    stacked = np.vstack((X1, X2))
    b = compute_b(stacked, omega, u)
    L_C = np.linalg.cholesky(compute_C(stacked, omega, u))
    
    b_updated = compute_b(X1, omega, u)
    L_C_updated = np.linalg.cholesky(compute_C(X1, omega, u))
    log_sum_weights1 = log_sum_exp(log_weights1)
    
    b_updated, L_C_updated = update_b_L_C_weighted(X2, b_updated, L_C_updated, log_sum_weights1, log_weights2, omega, u)
    
    assert_allclose(b, b_updated)
    assert_allclose(L_C, L_C_updated)

def test_update_L_C_naive_equals_batch():
    N = 100
    D = 3
    m = 10
    omega = np.random.randn(D, m)
    u = np.random.uniform(0, 2 * np.pi, m)
    X = np.random.randn(N, D)
    x = np.random.randn(D)
    
    L_C = np.linalg.cholesky(compute_C(X, omega, u))
    L_C = update_L_C_single(x, L_C, n=N, omega=omega, u=u)
    L_C_batch = np.linalg.cholesky(compute_C(np.vstack((X, x)), omega, u))

    assert_allclose(L_C, L_C_batch)

def test_update_L_C_equals_batch():
    N = 100
    D = 3
    m = 10
    omega = np.random.randn(D, m)
    u = np.random.uniform(0, 2 * np.pi, m)
    X = np.random.randn(N, D)
    x = np.random.randn(D)
    
    L_C = np.linalg.cholesky(compute_C(X, omega, u))
    L_C = update_L_C_single(x, L_C, N, omega, u)
    L_C_batch = np.linalg.cholesky(compute_C(np.vstack((X, x)), omega, u))

    assert_allclose(L_C, L_C_batch)

def test_hessian_execute():
    if not theano_available:
        raise SkipTest("Theano not available.")
    sigma = 1.
    lmbda = 1.
    N = 100
    D = 2
    m = 10
    X = np.random.randn(N, D)
    
    est = KernelExpFiniteGaussian(sigma, lmbda, m, D)
    est.fit(X)
    est.hessian(X[0])

def test_third_order_derivative_tensor_execute():

    if not theano_available: