mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-06-11 15:15:15 +02:00
Remove symbolic import.
This commit is contained in:
Neil Lawrence
commit
4a93e75ce2
28 changed files with 291 additions and 576 deletions
|
|
@ -213,7 +213,7 @@ class Model(Parameterized):
|
|||
def optimize(self, optimizer=None, start=None, **kwargs):
|
||||
"""
|
||||
Optimize the model using self.log_likelihood and self.log_likelihood_gradient, as well as self.priors.
|
||||
|
||||
|
||||
kwargs are passed to the optimizer. They can be:
|
||||
|
||||
:param max_f_eval: maximum number of function evaluations
|
||||
|
|
@ -231,7 +231,7 @@ class Model(Parameterized):
|
|||
- 'lbfgsb': the l-bfgs-b method (see scipy.optimize.fmin_l_bfgs_b),
|
||||
- 'sgd': stochastic gradient decsent (see scipy.optimize.sgd). For experts only!
|
||||
|
||||
|
||||
|
||||
"""
|
||||
if self.is_fixed:
|
||||
raise RuntimeError, "Cannot optimize, when everything is fixed"
|
||||
|
|
@ -300,7 +300,7 @@ class Model(Parameterized):
|
|||
# just check the global ratio
|
||||
dx = np.zeros(x.shape)
|
||||
dx[transformed_index] = step * (np.sign(np.random.uniform(-1, 1, transformed_index.size)) if transformed_index.size != 2 else 1.)
|
||||
|
||||
|
||||
# evaulate around the point x
|
||||
f1 = self._objective(x + dx)
|
||||
f2 = self._objective(x - dx)
|
||||
|
|
|
|||
|
|
@ -272,7 +272,7 @@ class Param(Parameterizable, ObsAr):
|
|||
</tr>"""
|
||||
header = header_format.format(x=self.hierarchy_name(), c=__constraints_name__, i=__index_name__, t=__tie_name__, p=__priors_name__) # nice header for printing
|
||||
if not ties: ties = itertools.cycle([''])
|
||||
return "\n".join(['<table>'] + [header] + ["<tr><td>{i}</td><td>{x}</td><td>{c}</td><td>{p}</td><td>{t}</td></tr>".format(x=x, c=" ".join(map(str, c)), p=" ".join(map(str, p)), t=(t or ''), i=i) for i, x, c, t, p in itertools.izip(indices, vals, constr_matrix, ties, prirs)] + ["</table>"])
|
||||
return "\n".join(['<table>'] + [header] + ["<tr><td>{i}</td><td align=\"right\">{x}</td><td>{c}</td><td>{p}</td><td>{t}</td></tr>".format(x=x, c=" ".join(map(str, c)), p=" ".join(map(str, p)), t=(t or ''), i=i) for i, x, c, t, p in itertools.izip(indices, vals, constr_matrix, ties, prirs)] + ["</table>"])
|
||||
|
||||
def __str__(self, constr_matrix=None, indices=None, prirs=None, ties=None, lc=None, lx=None, li=None, lp=None, lt=None, only_name=False):
|
||||
filter_ = self._current_slice_
|
||||
|
|
|
|||
|
|
@ -370,7 +370,7 @@ class Parameterized(Parameterizable):
|
|||
cl = max([len(str(x)) if x else 0 for x in constrs + ["Constraint"]])
|
||||
tl = max([len(str(x)) if x else 0 for x in ts + ["Tied to"]])
|
||||
pl = max([len(str(x)) if x else 0 for x in prirs + ["Prior"]])
|
||||
format_spec = "<tr><td>{{name:<{0}s}}</td><td>{{desc:>{1}s}}</td><td>{{const:^{2}s}}</td><td>{{pri:^{3}s}}</td><td>{{t:^{4}s}}</td></tr>".format(nl, sl, cl, pl, tl)
|
||||
format_spec = "<tr><td>{{name:<{0}s}}</td><td align=\"right\">{{desc:>{1}s}}</td><td>{{const:^{2}s}}</td><td>{{pri:^{3}s}}</td><td>{{t:^{4}s}}</td></tr>".format(nl, sl, cl, pl, tl)
|
||||
to_print = []
|
||||
for n, d, c, t, p in itertools.izip(names, desc, constrs, ts, prirs):
|
||||
to_print.append(format_spec.format(name=n, desc=d, const=c, t=t, pri=p))
|
||||
|
|
|
|||
|
|
@ -10,6 +10,10 @@ from parameterization.variational import VariationalPosterior
|
|||
|
||||
import logging
|
||||
from GPy.inference.latent_function_inference.posterior import Posterior
|
||||
from GPy.inference.optimization.stochastics import SparseGPStochastics,\
|
||||
SparseGPMissing
|
||||
#no stochastics.py file added! from GPy.inference.optimization.stochastics import SparseGPStochastics,\
|
||||
#SparseGPMissing
|
||||
logger = logging.getLogger("sparse gp")
|
||||
|
||||
class SparseGP(GP):
|
||||
|
|
@ -37,12 +41,7 @@ class SparseGP(GP):
|
|||
|
||||
def __init__(self, X, Y, Z, kernel, likelihood, inference_method=None,
|
||||
name='sparse gp', Y_metadata=None, normalizer=False,
|
||||
missing_data=False):
|
||||
|
||||
self.missing_data = missing_data
|
||||
if self.missing_data:
|
||||
self.ninan = ~np.isnan(Y)
|
||||
|
||||
missing_data=False, stochastic=False, batchsize=1):
|
||||
#pick a sensible inference method
|
||||
if inference_method is None:
|
||||
if isinstance(likelihood, likelihoods.Gaussian):
|
||||
|
|
@ -56,6 +55,22 @@ class SparseGP(GP):
|
|||
self.num_inducing = Z.shape[0]
|
||||
|
||||
GP.__init__(self, X, Y, kernel, likelihood, inference_method=inference_method, name=name, Y_metadata=Y_metadata, normalizer=normalizer)
|
||||
self.missing_data = missing_data
|
||||
|
||||
if stochastic and missing_data:
|
||||
self.missing_data = True
|
||||
self.ninan = ~np.isnan(Y)
|
||||
self.stochastics = SparseGPStochastics(self, batchsize)
|
||||
elif stochastic and not missing_data:
|
||||
self.missing_data = False
|
||||
self.stochastics = SparseGPStochastics(self, batchsize)
|
||||
elif missing_data:
|
||||
self.missing_data = True
|
||||
self.ninan = ~np.isnan(Y)
|
||||
self.stochastics = SparseGPMissing(self)
|
||||
else:
|
||||
self.stochastics = False
|
||||
|
||||
logger.info("Adding Z as parameter")
|
||||
self.link_parameter(self.Z, index=0)
|
||||
if self.missing_data:
|
||||
|
|
@ -71,6 +86,7 @@ class SparseGP(GP):
|
|||
print message,
|
||||
print ''
|
||||
|
||||
self.posterior = None
|
||||
|
||||
def has_uncertain_inputs(self):
|
||||
return isinstance(self.X, VariationalPosterior)
|
||||
|
|
@ -156,7 +172,31 @@ class SparseGP(GP):
|
|||
|
||||
value_indices:
|
||||
dictionary holding indices for the update in full_values.
|
||||
if the key exists the update rule is:
|
||||
if the key exists the update rule is:def df(x):
|
||||
m.stochastics.do_stochastics()
|
||||
grads = m._grads(x)
|
||||
print '\r',
|
||||
message = "Lik: {: 6.4E} Grad: {: 6.4E} Dim: {} Lik: {} Len: {!s}".format(float(m.log_likelihood()), np.einsum('i,i->', grads, grads), m.stochastics.d, float(m.likelihood.variance), " ".join(["{:3.2E}".format(l) for l in m.kern.lengthscale.values]))
|
||||
print message,
|
||||
return grads
|
||||
|
||||
def grad_stop(threshold):
|
||||
def inner(args):
|
||||
g = args['gradient']
|
||||
return np.sqrt(np.einsum('i,i->',g,g)) < threshold
|
||||
return inner
|
||||
|
||||
def maxiter_stop(maxiter):
|
||||
def inner(args):
|
||||
return args['n_iter'] == maxiter
|
||||
return inner
|
||||
|
||||
def optimize(m, maxiter=1000):
|
||||
#opt = climin.RmsProp(m.optimizer_array.copy(), df, 1e-6, decay=0.9, momentum=0.9, step_adapt=1e-7)
|
||||
opt = climin.Adadelta(m.optimizer_array.copy(), df, 1e-2, decay=0.9)
|
||||
ret = opt.minimize_until((grad_stop(.1), maxiter_stop(maxiter)))
|
||||
print
|
||||
return ret
|
||||
full_values[key][value_indices[key]] += current_values[key]
|
||||
"""
|
||||
for key in current_values.keys():
|
||||
|
|
@ -206,22 +246,29 @@ class SparseGP(GP):
|
|||
def _outer_loop_for_missing_data(self):
|
||||
Lm = None
|
||||
dL_dKmm = None
|
||||
Kmm = None
|
||||
|
||||
self._log_marginal_likelihood = 0
|
||||
full_values = self._outer_init_full_values()
|
||||
|
||||
woodbury_inv = np.zeros((self.num_inducing, self.num_inducing, self.output_dim))
|
||||
woodbury_vector = np.zeros((self.num_inducing, self.output_dim))
|
||||
m_f = lambda i: "Inference with missing data: {: >7.2%}".format(float(i+1)/self.output_dim)
|
||||
message = m_f(-1)
|
||||
print message,
|
||||
for d in xrange(self.output_dim):
|
||||
if self.posterior is None:
|
||||
woodbury_inv = np.zeros((self.num_inducing, self.num_inducing, self.output_dim))
|
||||
woodbury_vector = np.zeros((self.num_inducing, self.output_dim))
|
||||
else:
|
||||
woodbury_inv = self.posterior._woodbury_inv
|
||||
woodbury_vector = self.posterior._woodbury_vector
|
||||
|
||||
if not self.stochastics:
|
||||
m_f = lambda i: "Inference with missing_data: {: >7.2%}".format(float(i+1)/self.output_dim)
|
||||
message = m_f(-1)
|
||||
print message,
|
||||
|
||||
for d in self.stochastics.d:
|
||||
ninan = self.ninan[:, d]
|
||||
|
||||
print ' '*(len(message)) + '\r',
|
||||
message = m_f(d)
|
||||
print message,
|
||||
if not self.stochastics:
|
||||
print ' '*(len(message)) + '\r',
|
||||
message = m_f(d)
|
||||
print message,
|
||||
|
||||
posterior, log_marginal_likelihood, \
|
||||
grad_dict, current_values, value_indices = self._inner_parameters_changed(
|
||||
|
|
@ -236,19 +283,50 @@ class SparseGP(GP):
|
|||
|
||||
Lm = posterior.K_chol
|
||||
dL_dKmm = grad_dict['dL_dKmm']
|
||||
Kmm = posterior._K
|
||||
woodbury_inv[:, :, d] = posterior.woodbury_inv
|
||||
woodbury_vector[:, d:d+1] = posterior.woodbury_vector
|
||||
self._log_marginal_likelihood += log_marginal_likelihood
|
||||
print ''
|
||||
if not self.stochastics:
|
||||
print ''
|
||||
|
||||
self.posterior = Posterior(woodbury_inv=woodbury_inv, woodbury_vector=woodbury_vector,
|
||||
K=Kmm, mean=None, cov=None, K_chol=Lm)
|
||||
if self.posterior is None:
|
||||
self.posterior = Posterior(woodbury_inv=woodbury_inv, woodbury_vector=woodbury_vector,
|
||||
K=posterior._K, mean=None, cov=None, K_chol=posterior.K_chol)
|
||||
self._outer_values_update(full_values)
|
||||
|
||||
def _outer_loop_without_missing_data(self):
|
||||
self._log_marginal_likelihood = 0
|
||||
|
||||
if self.posterior is None:
|
||||
woodbury_inv = np.zeros((self.num_inducing, self.num_inducing, self.output_dim))
|
||||
woodbury_vector = np.zeros((self.num_inducing, self.output_dim))
|
||||
else:
|
||||
woodbury_inv = self.posterior._woodbury_inv
|
||||
woodbury_vector = self.posterior._woodbury_vector
|
||||
|
||||
d = self.stochastics.d
|
||||
posterior, log_marginal_likelihood, \
|
||||
grad_dict, current_values, _ = self._inner_parameters_changed(
|
||||
self.kern, self.X,
|
||||
self.Z, self.likelihood,
|
||||
self.Y_normalized[:, d], self.Y_metadata)
|
||||
self.grad_dict = grad_dict
|
||||
|
||||
self._log_marginal_likelihood += log_marginal_likelihood
|
||||
|
||||
self._outer_values_update(current_values)
|
||||
|
||||
woodbury_inv[:, :, d] = posterior.woodbury_inv[:, :, None]
|
||||
woodbury_vector[:, d] = posterior.woodbury_vector
|
||||
if self.posterior is None:
|
||||
self.posterior = Posterior(woodbury_inv=woodbury_inv, woodbury_vector=woodbury_vector,
|
||||
K=posterior._K, mean=None, cov=None, K_chol=posterior.K_chol)
|
||||
|
||||
def parameters_changed(self):
|
||||
if self.missing_data:
|
||||
self._outer_loop_for_missing_data()
|
||||
elif self.stochastics:
|
||||
self._outer_loop_without_missing_data()
|
||||
else:
|
||||
self.posterior, self._log_marginal_likelihood, self.grad_dict, full_values, _ = self._inner_parameters_changed(self.kern, self.X, self.Z, self.likelihood, self.Y_normalized, self.Y_metadata)
|
||||
self._outer_values_update(full_values)
|
||||
|
|
|
|||
|
|
@ -34,7 +34,8 @@ class SparseGP_MPI(SparseGP):
|
|||
|
||||
"""
|
||||
|
||||
def __init__(self, X, Y, Z, kernel, likelihood, variational_prior=None, inference_method=None, name='sparse gp mpi', Y_metadata=None, mpi_comm=None, normalizer=False, missing_data=False):
|
||||
def __init__(self, X, Y, Z, kernel, likelihood, variational_prior=None, inference_method=None, name='sparse gp mpi', Y_metadata=None, mpi_comm=None, normalizer=False,
|
||||
missing_data=False, stochastic=False, batchsize=1):
|
||||
self._IN_OPTIMIZATION_ = False
|
||||
if mpi_comm != None:
|
||||
if inference_method is None:
|
||||
|
|
@ -42,7 +43,8 @@ class SparseGP_MPI(SparseGP):
|
|||
else:
|
||||
assert isinstance(inference_method, VarDTC_minibatch), 'inference_method has to support MPI!'
|
||||
|
||||
super(SparseGP_MPI, self).__init__(X, Y, Z, kernel, likelihood, inference_method=inference_method, name=name, Y_metadata=Y_metadata, normalizer=normalizer, missing_data=missing_data)
|
||||
super(SparseGP_MPI, self).__init__(X, Y, Z, kernel, likelihood, inference_method=inference_method, name=name, Y_metadata=Y_metadata, normalizer=normalizer,
|
||||
missing_data=missing_data, stochastic=stochastic, batchsize=batchsize)
|
||||
self.update_model(False)
|
||||
self.link_parameter(self.X, index=0)
|
||||
if variational_prior is not None:
|
||||
|
|
|
|||
|
|
@ -12,7 +12,7 @@ from sympy.utilities.iterables import numbered_symbols
|
|||
import scipy
|
||||
import GPy
|
||||
#from scipy.special import gammaln, gamma, erf, erfc, erfcx, polygamma
|
||||
from GPy.util.symbolic import normcdf, normcdfln, logistic, logisticln, erfcx, erfc, gammaln
|
||||
#@NDL you removed this file! #from GPy.util.symbolic import normcdf, normcdfln, logistic, logisticln, erfcx, erfc, gammaln
|
||||
def getFromDict(dataDict, mapList):
|
||||
return reduce(lambda d, k: d[k], mapList, dataDict)
|
||||
|
||||
|
|
|
|||
|
|
@ -208,12 +208,12 @@ def _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim=False):
|
|||
Q_signal = 4
|
||||
import GPy
|
||||
import numpy as np
|
||||
np.random.seed(0)
|
||||
np.random.seed(3000)
|
||||
|
||||
k = GPy.kern.Matern32(Q_signal, 1., lengthscale=np.random.uniform(1,6,Q_signal), ARD=1)
|
||||
k = GPy.kern.Matern32(Q_signal, 10., lengthscale=1+(np.random.uniform(1,6,Q_signal)), ARD=1)
|
||||
t = np.c_[[np.linspace(-1,5,N) for _ in range(Q_signal)]].T
|
||||
K = k.K(t)
|
||||
s1, s2, s3, sS = np.random.multivariate_normal(np.zeros(K.shape[0]), K, size=(4))[:,:,None]
|
||||
s2, s1, s3, sS = np.random.multivariate_normal(np.zeros(K.shape[0]), K, size=(4))[:,:,None]
|
||||
|
||||
Y1, Y2, Y3, S1, S2, S3 = _generate_high_dimensional_output(D1, D2, D3, s1, s2, s3, sS)
|
||||
|
||||
|
|
@ -360,7 +360,6 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1,
|
|||
):
|
||||
from GPy import kern
|
||||
from GPy.models import BayesianGPLVM
|
||||
from GPy.inference.latent_function_inference.var_dtc import VarDTCMissingData
|
||||
|
||||
D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 400, 3, 4
|
||||
_, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim)
|
||||
|
|
|
|||
|
|
@ -82,7 +82,7 @@ class Laplace(LatentFunctionInference):
|
|||
|
||||
#define the objective function (to be maximised)
|
||||
def obj(Ki_f, f):
|
||||
return -0.5*np.dot(Ki_f.flatten(), f.flatten()) + likelihood.logpdf(f, Y, Y_metadata=Y_metadata)
|
||||
return -0.5*np.dot(Ki_f.flatten(), f.flatten()) + np.sum(likelihood.logpdf(f, Y, Y_metadata=Y_metadata))
|
||||
|
||||
difference = np.inf
|
||||
iteration = 0
|
||||
|
|
@ -152,7 +152,7 @@ class Laplace(LatentFunctionInference):
|
|||
Ki_W_i = K - C.T.dot(C)
|
||||
|
||||
#compute the log marginal
|
||||
log_marginal = -0.5*np.dot(Ki_f.flatten(), f_hat.flatten()) + likelihood.logpdf(f_hat, Y, Y_metadata=Y_metadata) - np.sum(np.log(np.diag(L)))
|
||||
log_marginal = -0.5*np.dot(Ki_f.flatten(), f_hat.flatten()) + likelihood.logpdf(f_hat, Y, Y_metadata=Y_metadata).sum() - np.sum(np.log(np.diag(L)))
|
||||
|
||||
# Compute matrices for derivatives
|
||||
dW_df = -likelihood.d3logpdf_df3(f_hat, Y, Y_metadata=Y_metadata) # -d3lik_d3fhat
|
||||
|
|
|
|||
48
GPy/inference/optimization/stochastics.py
Normal file
48
GPy/inference/optimization/stochastics.py
Normal file
|
|
@ -0,0 +1,48 @@
|
|||
'''
|
||||
Created on 9 Oct 2014
|
||||
|
||||
@author: maxz
|
||||
'''
|
||||
|
||||
class StochasticStorage(object):
|
||||
'''
|
||||
This is a container for holding the stochastic parameters,
|
||||
such as subset indices or step length and so on.
|
||||
'''
|
||||
def __init__(self, model):
|
||||
"""
|
||||
Initialize this stochastic container using the given model
|
||||
"""
|
||||
|
||||
def do_stochastics(self):
|
||||
"""
|
||||
Update the internal state to the next batch of the stochastic
|
||||
descent algorithm.
|
||||
"""
|
||||
pass
|
||||
|
||||
class SparseGPMissing(StochasticStorage):
|
||||
def __init__(self, model, batchsize=1):
|
||||
self.d = xrange(model.Y_normalized.shape[1])
|
||||
|
||||
class SparseGPStochastics(StochasticStorage):
|
||||
"""
|
||||
For the sparse gp we need to store the dimension we are in,
|
||||
and the indices corresponding to those
|
||||
"""
|
||||
def __init__(self, model, batchsize=1):
|
||||
import itertools
|
||||
self.batchsize = batchsize
|
||||
if self.batchsize == 1:
|
||||
self.dimensions = itertools.cycle(range(model.Y_normalized.shape[1]))
|
||||
else:
|
||||
import numpy as np
|
||||
self.dimensions = lambda: np.random.choice(model.Y_normalized.shape[1], size=batchsize, replace=False)
|
||||
self.d = None
|
||||
self.do_stochastics()
|
||||
|
||||
def do_stochastics(self):
|
||||
if self.batchsize == 1:
|
||||
self.d = [self.dimensions.next()]
|
||||
else:
|
||||
self.d = self.dimensions()
|
||||
|
|
@ -1,30 +0,0 @@
|
|||
try:
|
||||
import sympy as sym
|
||||
sympy_available=True
|
||||
except ImportError:
|
||||
sympy_available=False
|
||||
|
||||
import numpy as np
|
||||
from symbolic import Symbolic
|
||||
|
||||
class Eq(Symbolic):
|
||||
"""
|
||||
The exponentiated quadratic covariance as a symbolic function.
|
||||
|
||||
"""
|
||||
def __init__(self, input_dim, output_dim=1, variance=1.0, lengthscale=1.0, name='Eq'):
|
||||
|
||||
parameters = {'variance' : variance, 'lengthscale' : lengthscale}
|
||||
x = sym.symbols('x_:' + str(input_dim))
|
||||
z = sym.symbols('z_:' + str(input_dim))
|
||||
variance = sym.var('variance',positive=True)
|
||||
lengthscale = sym.var('lengthscale', positive=True)
|
||||
dist_string = ' + '.join(['(x_%i - z_%i)**2' %(i, i) for i in range(input_dim)])
|
||||
from sympy.parsing.sympy_parser import parse_expr
|
||||
dist = parse_expr(dist_string)
|
||||
|
||||
# this is the covariance function
|
||||
f = variance*sym.exp(-dist/(2*lengthscale**2))
|
||||
# extra input dim is to signify the output dimension.
|
||||
super(Eq, self).__init__(input_dim=input_dim, k=f, output_dim=output_dim, parameters=parameters, name=name)
|
||||
|
||||
|
|
@ -1,30 +0,0 @@
|
|||
try:
|
||||
import sympy as sym
|
||||
sympy_available=True
|
||||
except ImportError:
|
||||
sympy_available=False
|
||||
|
||||
import numpy as np
|
||||
from symbolic import Symbolic
|
||||
|
||||
class Heat_eqinit(Symbolic):
|
||||
"""
|
||||
A symbolic covariance based on laying down an initial condition of the heat equation with an exponentiated quadratic covariance. The covariance then has multiple outputs which are interpreted as observations of the diffused process with different diffusion coefficients (or at different times).
|
||||
|
||||
"""
|
||||
def __init__(self, input_dim, output_dim=1, param=None, name='Heat_eqinit'):
|
||||
|
||||
x = sym.symbols('x_:' + str(input_dim))
|
||||
z = sym.symbols('z_:' + str(input_dim))
|
||||
scale = sym.var('scale_i scale_j',positive=True)
|
||||
lengthscale = sym.var('lengthscale_i lengthscale_j', positive=True)
|
||||
shared_lengthscale = sym.var('shared_lengthscale', positive=True)
|
||||
dist_string = ' + '.join(['(x_%i - z_%i)**2' %(i, i) for i in range(input_dim)])
|
||||
from sympy.parsing.sympy_parser import parse_expr
|
||||
dist = parse_expr(dist_string)
|
||||
|
||||
# this is the covariance function
|
||||
f = scale_i*scale_j*sym.exp(-dist/(2*(shared_lengthscale**2 + lengthscale_i*lengthscale_j)))
|
||||
# extra input dim is to signify the output dimension.
|
||||
super(Heat_eqinit, self).__init__(input_dim=input_dim+1, k=f, output_dim=output_dim, name=name)
|
||||
|
||||
|
|
@ -6,6 +6,7 @@ import numpy as np
|
|||
from ...core.parameterization.parameterized import Parameterized
|
||||
from kernel_slice_operations import KernCallsViaSlicerMeta
|
||||
from ...util.caching import Cache_this
|
||||
from GPy.core.parameterization.observable_array import ObsAr
|
||||
|
||||
|
||||
|
||||
|
|
@ -54,6 +55,20 @@ class Kern(Parameterized):
|
|||
|
||||
self._sliced_X = 0
|
||||
self.useGPU = self._support_GPU and useGPU
|
||||
self._return_psi2_n_flag = ObsAr(np.zeros(1)).astype(bool)
|
||||
|
||||
@property
|
||||
def return_psi2_n(self):
|
||||
"""
|
||||
Flag whether to pass back psi2 as NxMxM or MxM, by summing out N.
|
||||
"""
|
||||
return self._return_psi2_n_flag[0]
|
||||
@return_psi2_n.setter
|
||||
def return_psi2_n(self, val):
|
||||
def visit(self):
|
||||
if isinstance(self, Kern):
|
||||
self._return_psi2_n_flag[0]=val
|
||||
self.traverse(visit)
|
||||
|
||||
@Cache_this(limit=20)
|
||||
def _slice_X(self, X):
|
||||
|
|
@ -162,7 +177,7 @@ class Kern(Parameterized):
|
|||
def __mul__(self, other):
|
||||
""" Here we overload the '*' operator. See self.prod for more information"""
|
||||
return self.prod(other)
|
||||
|
||||
|
||||
def __imul__(self, other):
|
||||
""" Here we overload the '*' operator. See self.prod for more information"""
|
||||
return self.prod(other)
|
||||
|
|
|
|||
|
|
@ -10,7 +10,6 @@ import sslinear_psi_comp
|
|||
import linear_psi_comp
|
||||
|
||||
class PSICOMP_RBF(Pickleable):
|
||||
|
||||
@Cache_this(limit=2, ignore_args=(0,))
|
||||
def psicomputations(self, variance, lengthscale, Z, variational_posterior):
|
||||
if isinstance(variational_posterior, variational.NormalPosterior):
|
||||
|
|
@ -19,7 +18,7 @@ class PSICOMP_RBF(Pickleable):
|
|||
return ssrbf_psi_comp.psicomputations(variance, lengthscale, Z, variational_posterior)
|
||||
else:
|
||||
raise ValueError, "unknown distriubtion received for psi-statistics"
|
||||
|
||||
|
||||
@Cache_this(limit=2, ignore_args=(0,1,2,3))
|
||||
def psiDerivativecomputations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, lengthscale, Z, variational_posterior):
|
||||
if isinstance(variational_posterior, variational.NormalPosterior):
|
||||
|
|
@ -28,10 +27,10 @@ class PSICOMP_RBF(Pickleable):
|
|||
return ssrbf_psi_comp.psiDerivativecomputations(dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, lengthscale, Z, variational_posterior)
|
||||
else:
|
||||
raise ValueError, "unknown distriubtion received for psi-statistics"
|
||||
|
||||
|
||||
def _setup_observers(self):
|
||||
pass
|
||||
|
||||
|
||||
class PSICOMP_Linear(Pickleable):
|
||||
|
||||
@Cache_this(limit=2, ignore_args=(0,))
|
||||
|
|
@ -42,7 +41,7 @@ class PSICOMP_Linear(Pickleable):
|
|||
return sslinear_psi_comp.psicomputations(variance, Z, variational_posterior)
|
||||
else:
|
||||
raise ValueError, "unknown distriubtion received for psi-statistics"
|
||||
|
||||
|
||||
@Cache_this(limit=2, ignore_args=(0,1,2,3))
|
||||
def psiDerivativecomputations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, Z, variational_posterior):
|
||||
if isinstance(variational_posterior, variational.NormalPosterior):
|
||||
|
|
@ -51,6 +50,6 @@ class PSICOMP_Linear(Pickleable):
|
|||
return sslinear_psi_comp.psiDerivativecomputations(dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, Z, variational_posterior)
|
||||
else:
|
||||
raise ValueError, "unknown distriubtion received for psi-statistics"
|
||||
|
||||
|
||||
def _setup_observers(self):
|
||||
pass
|
||||
|
|
@ -139,7 +139,7 @@ def _psi2compDer(dL_dpsi2, variance, lengthscale, Z, mu, S):
|
|||
denom2 = np.square(denom)
|
||||
|
||||
_psi2 = _psi2computations(variance, lengthscale, Z, mu, S) # NxMxM
|
||||
Lpsi2 = dL_dpsi2[None,:,:]*_psi2
|
||||
Lpsi2 = dL_dpsi2*_psi2 # dL_dpsi2 is MxM, using broadcast to multiply N out
|
||||
Lpsi2sum = np.einsum('nmo->n',Lpsi2) #N
|
||||
Lpsi2Z = np.einsum('nmo,oq->nq',Lpsi2,Z) #NxQ
|
||||
Lpsi2Z2 = np.einsum('nmo,oq,oq->nq',Lpsi2,Z,Z) #NxQ
|
||||
|
|
|
|||
|
|
@ -8,9 +8,9 @@ from ...core.parameterization.transformations import Logexp
|
|||
from ...util.linalg import tdot
|
||||
from ... import util
|
||||
import numpy as np
|
||||
from scipy import integrate
|
||||
from ...util.caching import Cache_this
|
||||
from scipy import integrate, weave
|
||||
from ...util.config import config # for assesing whether to use weave
|
||||
from ...util.caching import Cache_this
|
||||
|
||||
class Stationary(Kern):
|
||||
"""
|
||||
|
|
@ -72,6 +72,13 @@ class Stationary(Kern):
|
|||
|
||||
@Cache_this(limit=5, ignore_args=())
|
||||
def K(self, X, X2=None):
|
||||
"""
|
||||
Kernel function applied on inputs X and X2.
|
||||
In the stationary case there is an inner function depending on the
|
||||
distances from X to X2, called r.
|
||||
|
||||
K(X, X2) = K_of_r((X-X2)**2)
|
||||
"""
|
||||
r = self._scaled_dist(X, X2)
|
||||
return self.K_of_r(r)
|
||||
|
||||
|
|
@ -125,10 +132,22 @@ class Stationary(Kern):
|
|||
return ret
|
||||
|
||||
def update_gradients_diag(self, dL_dKdiag, X):
|
||||
"""
|
||||
Given the derivative of the objective with respect to the diagonal of
|
||||
the covariance matrix, compute the derivative wrt the parameters of
|
||||
this kernel and stor in the <parameter>.gradient field.
|
||||
|
||||
See also update_gradients_full
|
||||
"""
|
||||
self.variance.gradient = np.sum(dL_dKdiag)
|
||||
self.lengthscale.gradient = 0.
|
||||
|
||||
def update_gradients_full(self, dL_dK, X, X2=None):
|
||||
"""
|
||||
Given the derivative of the objective wrt the covariance matrix
|
||||
(dL_dK), compute the gradient wrt the parameters of this kernel,
|
||||
and store in the parameters object as e.g. self.variance.gradient
|
||||
"""
|
||||
self.variance.gradient = np.einsum('ij,ij,i', self.K(X, X2), dL_dK, 1./self.variance)
|
||||
|
||||
#now the lengthscale gradient(s)
|
||||
|
|
@ -166,6 +185,7 @@ class Stationary(Kern):
|
|||
return 1./np.where(dist != 0., dist, np.inf)
|
||||
|
||||
def weave_lengthscale_grads(self, tmp, X, X2):
|
||||
"""Use scipy.weave to compute derivatives wrt the lengthscales"""
|
||||
N,M = tmp.shape
|
||||
Q = X.shape[1]
|
||||
if hasattr(X, 'values'):X = X.values
|
||||
|
|
@ -183,7 +203,6 @@ class Stationary(Kern):
|
|||
grads[q] = gradq;
|
||||
}
|
||||
"""
|
||||
from scipy import weave
|
||||
weave.inline(code, ['tmp', 'X', 'X2', 'grads', 'N', 'M', 'Q'], type_converters=weave.converters.blitz, support_code="#include <math.h>")
|
||||
return -grads/self.lengthscale**3
|
||||
|
||||
|
|
@ -383,7 +402,7 @@ class Matern52(Stationary):
|
|||
|
||||
class ExpQuad(Stationary):
|
||||
"""
|
||||
The Exponentiated quadratic covariance function.
|
||||
The Exponentiated quadratic covariance function.
|
||||
|
||||
.. math::
|
||||
|
||||
|
|
|
|||
|
|
@ -136,10 +136,8 @@ class Bernoulli(Likelihood):
|
|||
assert np.atleast_1d(inv_link_f).shape == np.atleast_1d(y).shape
|
||||
#objective = y*np.log(inv_link_f) + (1.-y)*np.log(inv_link_f)
|
||||
state = np.seterr(divide='ignore')
|
||||
# TODO check y \in {0, 1} or {-1, 1}
|
||||
objective = np.where(y==1, np.log(inv_link_f), np.log(1-inv_link_f))
|
||||
np.seterr(**state)
|
||||
return np.sum(objective)
|
||||
p = np.where(y==1, inv_link_f, 1.-inv_link_f)
|
||||
return np.log(p)
|
||||
|
||||
def dlogpdf_dlink(self, inv_link_f, y, Y_metadata=None):
|
||||
"""
|
||||
|
|
|
|||
|
|
@ -131,6 +131,40 @@ class Likelihood(Parameterized):
|
|||
|
||||
return z, mean, variance
|
||||
|
||||
def variational_expectations(self, Y, m, v, gh_points=None):
|
||||
"""
|
||||
Use Gauss-Hermite Quadrature to compute
|
||||
|
||||
E_p(f) [ log p(y|f) ]
|
||||
d/dm E_p(f) [ log p(y|f) ]
|
||||
d/dv E_p(f) [ log p(y|f) ]
|
||||
|
||||
where p(f) is a Gaussian with mean m and variance v. The shapes of Y, m and v should match.
|
||||
|
||||
if no gh_points are passed, we construct them using defualt options
|
||||
"""
|
||||
|
||||
if gh_points is None:
|
||||
gh_x, gh_w = np.polynomial.hermite.hermgauss(20)
|
||||
|
||||
shape = m.shape
|
||||
m,v,Y = m.flatten(), v.flatten(), Y.flatten()
|
||||
|
||||
#make a grid of points
|
||||
X = gh_x[None,:]*np.sqrt(2.*v[:,None]) + m[:,None]
|
||||
logp = self.logpdf(X,Y[:,None])
|
||||
|
||||
p = np.clip(p, 1e-9, 1.-1e-9) # for numerical stability
|
||||
N = stats.norm.pdf(X)
|
||||
F = np.log(p).dot(self.gh_w)
|
||||
NoverP = N/p
|
||||
dF_dm = (NoverP*self.Ysign[:,None]).dot(self.gh_w)
|
||||
dF_dv = -0.5*(NoverP**2 + NoverP*X).dot(self.gh_w)
|
||||
return F, dF_dm, dF_dv
|
||||
|
||||
|
||||
|
||||
|
||||
def predictive_mean(self, mu, variance, Y_metadata=None):
|
||||
"""
|
||||
Quadrature calculation of the predictive mean: E(Y_star|Y) = E( E(Y_star|f_star, Y) )
|
||||
|
|
|
|||
|
|
@ -3,7 +3,7 @@
|
|||
|
||||
|
||||
import sympy as sym
|
||||
from GPy.util.symbolic import normcdfln
|
||||
#from GPy.util.symbolic import normcdfln
|
||||
import numpy as np
|
||||
import link_functions
|
||||
from symbolic import Symbolic
|
||||
|
|
|
|||
|
|
@ -4,7 +4,7 @@
|
|||
|
||||
import sympy as sym
|
||||
from sympy.utilities.lambdify import lambdify
|
||||
from GPy.util.symbolic import gammaln
|
||||
# does not exist! JH from GPy.util.symbolic import gammaln
|
||||
|
||||
import numpy as np
|
||||
import link_functions
|
||||
|
|
|
|||
|
|
@ -5,10 +5,8 @@ import numpy as np
|
|||
from .. import kern
|
||||
from ..core.sparse_gp_mpi import SparseGP_MPI
|
||||
from ..likelihoods import Gaussian
|
||||
from ..inference.optimization import SCG
|
||||
from ..util import linalg
|
||||
from ..core.parameterization.variational import NormalPosterior, NormalPrior, VariationalPosterior
|
||||
from ..inference.latent_function_inference.var_dtc_parallel import update_gradients, VarDTC_minibatch
|
||||
from ..core.parameterization.variational import NormalPosterior, NormalPrior
|
||||
from ..inference.latent_function_inference.var_dtc_parallel import VarDTC_minibatch
|
||||
from ..inference.latent_function_inference.var_dtc_gpu import VarDTC_GPU
|
||||
import logging
|
||||
|
||||
|
|
@ -27,7 +25,7 @@ class BayesianGPLVM(SparseGP_MPI):
|
|||
def __init__(self, Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10,
|
||||
Z=None, kernel=None, inference_method=None, likelihood=None,
|
||||
name='bayesian gplvm', mpi_comm=None, normalizer=None,
|
||||
missing_data=False):
|
||||
missing_data=False, stochastic=False, batchsize=1):
|
||||
self.mpi_comm = mpi_comm
|
||||
self.__IN_OPTIMIZATION__ = False
|
||||
|
||||
|
|
@ -77,7 +75,8 @@ class BayesianGPLVM(SparseGP_MPI):
|
|||
name=name, inference_method=inference_method,
|
||||
normalizer=normalizer, mpi_comm=mpi_comm,
|
||||
variational_prior=self.variational_prior,
|
||||
missing_data=missing_data)
|
||||
missing_data=missing_data, stochastic=stochastic,
|
||||
batchsize=batchsize)
|
||||
|
||||
def set_X_gradients(self, X, X_grad):
|
||||
"""Set the gradients of the posterior distribution of X in its specific form."""
|
||||
|
|
@ -90,7 +89,12 @@ class BayesianGPLVM(SparseGP_MPI):
|
|||
def _inner_parameters_changed(self, kern, X, Z, likelihood, Y, Y_metadata, Lm=None, dL_dKmm=None, subset_indices=None):
|
||||
posterior, log_marginal_likelihood, grad_dict, current_values, value_indices = super(BayesianGPLVM, self)._inner_parameters_changed(kern, X, Z, likelihood, Y, Y_metadata, Lm=Lm, dL_dKmm=dL_dKmm, subset_indices=subset_indices)
|
||||
|
||||
log_marginal_likelihood -= self.variational_prior.KL_divergence(X)
|
||||
kl_fctr = 1.
|
||||
if self.missing_data:
|
||||
d = self.output_dim
|
||||
log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(X)/d
|
||||
else:
|
||||
log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(X)
|
||||
|
||||
current_values['meangrad'], current_values['vargrad'] = self.kern.gradients_qX_expectations(
|
||||
variational_posterior=X,
|
||||
|
|
@ -104,8 +108,12 @@ class BayesianGPLVM(SparseGP_MPI):
|
|||
X.variance.gradient[:] = 0
|
||||
|
||||
self.variational_prior.update_gradients_KL(X)
|
||||
current_values['meangrad'] += X.mean.gradient
|
||||
current_values['vargrad'] += X.variance.gradient
|
||||
if self.missing_data:
|
||||
current_values['meangrad'] += kl_fctr*X.mean.gradient/d
|
||||
current_values['vargrad'] += kl_fctr*X.variance.gradient/d
|
||||
else:
|
||||
current_values['meangrad'] += kl_fctr*X.mean.gradient
|
||||
current_values['vargrad'] += kl_fctr*X.variance.gradient
|
||||
|
||||
if subset_indices is not None:
|
||||
value_indices['meangrad'] = subset_indices['samples']
|
||||
|
|
|
|||
|
|
@ -353,7 +353,7 @@ class TestNoiseModels(object):
|
|||
#print model._get_params()
|
||||
np.testing.assert_almost_equal(
|
||||
model.pdf(f.copy(), Y.copy()),
|
||||
np.exp(model.logpdf(f.copy(), Y.copy()))
|
||||
np.exp(model.logpdf(f.copy(), Y.copy()).sum())
|
||||
)
|
||||
|
||||
@with_setup(setUp, tearDown)
|
||||
|
|
|
|||
|
|
@ -139,16 +139,16 @@ class MiscTests(unittest.TestCase):
|
|||
np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
m.kern.randomize()
|
||||
m2.kern[''] = m.kern[:]
|
||||
np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
np.testing.assert_almost_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
m.kern.randomize()
|
||||
m2.kern[:] = m.kern[:]
|
||||
np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
np.testing.assert_almost_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
m.kern.randomize()
|
||||
m2.kern[''] = m.kern['']
|
||||
np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
np.testing.assert_almost_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
m.kern.randomize()
|
||||
m2.kern[:] = m.kern[''].values()
|
||||
np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
np.testing.assert_almost_equal(m.log_likelihood(), m2.log_likelihood())
|
||||
|
||||
def test_big_model(self):
|
||||
m = GPy.examples.dimensionality_reduction.mrd_simulation(optimize=0, plot=0, plot_sim=0)
|
||||
|
|
|
|||
|
|
@ -25,5 +25,6 @@ try:
|
|||
except ImportError as e:
|
||||
_sympy_available = False
|
||||
|
||||
if _sympy_available:
|
||||
import symbolic
|
||||
#@NDL: you deleted the symbolic.py file!
|
||||
#if _sympy_available:
|
||||
#import symbolic
|
||||
|
|
|
|||
|
|
@ -188,4 +188,6 @@ class Cache_this(object):
|
|||
self.ignore_args = ignore_args
|
||||
self.force_args = force_kwargs
|
||||
def __call__(self, f):
|
||||
return Cacher_wrap(f, self.limit, self.ignore_args, self.force_args)
|
||||
newf = Cacher_wrap(f, self.limit, self.ignore_args, self.force_args)
|
||||
update_wrapper(newf, f)
|
||||
return newf
|
||||
|
|
|
|||
|
|
@ -964,7 +964,7 @@ def singlecell_rna_seq_deng(dataset='singlecell_deng'):
|
|||
if not data_available(dataset):
|
||||
download_data(dataset)
|
||||
|
||||
from pandas import read_csv
|
||||
from pandas import read_csv, isnull
|
||||
dir_path = os.path.join(data_path, dataset)
|
||||
|
||||
# read the info .soft
|
||||
|
|
@ -983,6 +983,21 @@ def singlecell_rna_seq_deng(dataset='singlecell_deng'):
|
|||
c[1:4] = ['strain', 'cross', 'developmental_stage']
|
||||
sample_info.columns = c
|
||||
|
||||
# get the labels right:
|
||||
rep = re.compile('\(.*\)')
|
||||
def filter_dev_stage(row):
|
||||
if isnull(row):
|
||||
row = "2-cell stage embryo"
|
||||
if row.startswith("developmental stage: "):
|
||||
row = row[len("developmental stage: "):]
|
||||
if row == 'adult':
|
||||
row += " liver"
|
||||
row = row.replace(' stage ', ' ')
|
||||
row = rep.sub(' ', row)
|
||||
row = row.strip(' ')
|
||||
return row
|
||||
labels = sample_info.developmental_stage.apply(filter_dev_stage)
|
||||
|
||||
# Extract the tar file
|
||||
filename = os.path.join(dir_path, 'GSE45719_Raw.tar')
|
||||
with tarfile.open(filename, 'r') as files:
|
||||
|
|
@ -1016,8 +1031,7 @@ def singlecell_rna_seq_deng(dataset='singlecell_deng'):
|
|||
sample_info.index = data.index
|
||||
|
||||
# get the labels from the description
|
||||
rep = re.compile('fibroblast|\d+-cell|embryo|liver|blastocyst|blastomere|zygote', re.IGNORECASE)
|
||||
labels = sample_info.developmental_stage.apply(lambda row: " ".join(rep.findall(row)))
|
||||
#rep = re.compile('fibroblast|\d+-cell|embryo|liver|early blastocyst|mid blastocyst|late blastocyst|blastomere|zygote', re.IGNORECASE)
|
||||
|
||||
sys.stdout.write(' '*len(message) + '\r')
|
||||
sys.stdout.flush()
|
||||
|
|
|
|||
File diff suppressed because one or more lines are too long
|
|
@ -13,7 +13,7 @@ def initialize_latent(init, input_dim, Y):
|
|||
p = pca(Y)
|
||||
PC = p.project(Y, min(input_dim, Y.shape[1]))
|
||||
Xr[:PC.shape[0], :PC.shape[1]] = PC
|
||||
var = p.fracs[:input_dim]
|
||||
var = .1*p.fracs[:input_dim]
|
||||
else:
|
||||
var = Xr.var(0)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,367 +0,0 @@
|
|||
import sys
|
||||
import numpy as np
|
||||
import sympy as sym
|
||||
from sympy import Function, S, oo, I, cos, sin, asin, log, erf, pi, exp, sqrt, sign, gamma, polygamma
|
||||
from sympy.matrices import Matrix
|
||||
########################################
|
||||
## Try to do some matrix functions: problem, you can't do derivatives
|
||||
## with respect to matrix functions :-(
|
||||
|
||||
class GPySymMatrix(Matrix):
|
||||
def __init__(self, indices):
|
||||
Matrix.__init__(self)
|
||||
def atoms(self):
|
||||
return [e2 for e in self for e2 in e.atoms()]
|
||||
|
||||
class selector(Function):
|
||||
"""A function that returns an element of a Matrix depending on input indices."""
|
||||
nargs = 3
|
||||
def fdiff(self, argindex=1):
|
||||
return selector(*self.args)
|
||||
@classmethod
|
||||
def eval(cls, X, i, j):
|
||||
if i.is_Number and j.is_Number:
|
||||
return X[i, j]
|
||||
|
||||
##################################################
|
||||
class logistic(Function):
|
||||
"""The logistic function as a symbolic function."""
|
||||
nargs = 1
|
||||
def fdiff(self, argindex=1):
|
||||
x = self.args[0]
|
||||
return logistic(x)*(1-logistic(x))
|
||||
|
||||
@classmethod
|
||||
def eval(cls, x):
|
||||
if x.is_Number:
|
||||
return 1/(1+exp(-x))
|
||||
|
||||
class logisticln(Function):
|
||||
"""The log logistic, which can often be computed with more precision than the simply taking log(logistic(x)) when x is small or large."""
|
||||
nargs = 1
|
||||
|
||||
def fdiff(self, argindex=1):
|
||||
x = self.args[0]
|
||||
return 1-logistic(x)
|
||||
|
||||
@classmethod
|
||||
def eval(cls, x):
|
||||
if x.is_Number:
|
||||
return -np.log(1+exp(-x))
|
||||
|
||||
class erfc(Function):
|
||||
"""The complementary error function, erfc(x) = 1-erf(x). Used as a helper function, particularly for erfcx, the scaled complementary error function. and the normal distributions cdf."""
|
||||
nargs = 1
|
||||
|
||||
@classmethod
|
||||
def eval(cls, arg):
|
||||
return 1-erf(arg)
|
||||
|
||||
class erfcx(Function):
|
||||
nargs = 1
|
||||
def fdiff(self, argindex=1):
|
||||
x = self.args[0]
|
||||
return x*erfcx(x)-2/sqrt(pi)
|
||||
|
||||
@classmethod
|
||||
def eval(cls, x):
|
||||
if x.is_Number:
|
||||
return exp(x**2)*erfc(x)
|
||||
|
||||
class gammaln(Function):
|
||||
"""The log of the gamma function, which is often needed instead of log(gamma(x)) for better accuracy for large x."""
|
||||
nargs = 1
|
||||
|
||||
def fdiff(self, argindex=1):
|
||||
x=self.args[0]
|
||||
return polygamma(0, x)
|
||||
|
||||
@classmethod
|
||||
def eval(cls, x):
|
||||
if x.is_Number:
|
||||
return log(gamma(x))
|
||||
|
||||
|
||||
class normcdfln(Function):
|
||||
"""The log of the normal cdf. Can often be computed with better accuracy than log(normcdf(x)), particulary when x is either small or large."""
|
||||
nargs = 1
|
||||
|
||||
def fdiff(self, argindex=1):
|
||||
x = self.args[0]
|
||||
#return -erfcx(-x/sqrt(2))/sqrt(2*pi)
|
||||
#return exp(-normcdfln(x) - 0.5*x*x)/sqrt(2*pi)
|
||||
return sqrt(2/pi)*1/erfcx(-x/sqrt(2))
|
||||
|
||||
@classmethod
|
||||
def eval(cls, x):
|
||||
if x.is_Number:
|
||||
return log(normcdf(x))
|
||||
|
||||
def _eval_is_real(self):
|
||||
return self.args[0].is_real
|
||||
|
||||
class normcdf(Function):
|
||||
"""The cumulative distribution function of the standard normal. Provided as a convenient helper function. It is computed throught -0.5*erfc(-x/sqrt(2))."""
|
||||
nargs = 1
|
||||
def fdiff(self, argindex=1):
|
||||
x = self.args[0]
|
||||
return gaussian(x)
|
||||
|
||||
@classmethod
|
||||
def eval(cls, x):
|
||||
if x.is_Number:
|
||||
return 0.5*(erfc(-x/sqrt(2)))
|
||||
|
||||
def _eval_is_real(self):
|
||||
return self.args[0].is_real
|
||||
|
||||
class normalln(Function):
|
||||
"""The log of the standard normal distribution."""
|
||||
nargs = 1
|
||||
def fdiff(self, argindex=1):
|
||||
x = self.args[0]
|
||||
return -x
|
||||
|
||||
@classmethod
|
||||
def eval(cls, x):
|
||||
if x.is_Number:
|
||||
return 0.5*sqrt(2*pi) - 0.5*x*x
|
||||
|
||||
|
||||
def _eval_is_real(self):
|
||||
return True
|
||||
|
||||
|
||||
class normal(Function):
|
||||
"""The standard normal distribution. Provided as a convenience function."""
|
||||
nargs = 1
|
||||
@classmethod
|
||||
def eval(cls, x):
|
||||
return 1/sqrt(2*pi)*exp(-0.5*x*x)
|
||||
|
||||
def _eval_is_real(self):
|
||||
return True
|
||||
|
||||
class differfln(Function):
|
||||
nargs = 2
|
||||
|
||||
def fdiff(self, argindex=2):
|
||||
if argindex == 2:
|
||||
x0, x1 = self.args
|
||||
return -2/(sqrt(pi)*(erfcx(x1)-exp(x1**2-x0**2)*erfcx(x0)))
|
||||
elif argindex == 1:
|
||||
x0, x1 = self.args
|
||||
return 2/(sqrt(pi)*(exp(x0**2-x1**2)*erfcx(x1)-erfcx(x0)))
|
||||
else:
|
||||
raise ArgumentIndexError(self, argindex)
|
||||
|
||||
@classmethod
|
||||
def eval(cls, x0, x1):
|
||||
if x0.is_Number and x1.is_Number:
|
||||
return log(erfc(x1)-erfc(x0))
|
||||
|
||||
class dh_dd_i(Function):
|
||||
nargs = 5
|
||||
@classmethod
|
||||
def eval(cls, t, tprime, d_i, d_j, l):
|
||||
if (t.is_Number
|
||||
and tprime.is_Number
|
||||
and d_i.is_Number
|
||||
and d_j.is_Number
|
||||
and l.is_Number):
|
||||
|
||||
diff_t = (t-tprime)
|
||||
l2 = l*l
|
||||
h = h(t, tprime, d_i, d_j, l)
|
||||
half_l_di = 0.5*l*d_i
|
||||
arg_1 = half_l_di + tprime/l
|
||||
arg_2 = half_l_di - (t-tprime)/l
|
||||
ln_part_1 = ln_diff_erf(arg_1, arg_2)
|
||||
arg_1 = half_l_di
|
||||
arg_2 = half_l_di - t/l
|
||||
sign_val = sign(t/l)
|
||||
ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l)
|
||||
|
||||
base = ((0.5*d_i*l2*(d_i+d_j)-1)*h
|
||||
+ (-diff_t*sign_val*exp(half_l_di*half_l_di
|
||||
-d_i*diff_t
|
||||
+ln_part_1)
|
||||
+t*sign_val*exp(half_l_di*half_l_di
|
||||
-d_i*t-d_j*tprime
|
||||
+ln_part_2))
|
||||
+ l/sqrt(pi)*(-exp(-diff_t*diff_t/l2)
|
||||
+exp(-tprime*tprime/l2-d_i*t)
|
||||
+exp(-t*t/l2-d_j*tprime)
|
||||
-exp(-(d_i*t + d_j*tprime))))
|
||||
return base/(d_i+d_j)
|
||||
|
||||
class dh_dd_j(Function):
|
||||
nargs = 5
|
||||
@classmethod
|
||||
def eval(cls, t, tprime, d_i, d_j, l):
|
||||
if (t.is_Number
|
||||
and tprime.is_Number
|
||||
and d_i.is_Number
|
||||
and d_j.is_Number
|
||||
and l.is_Number):
|
||||
diff_t = (t-tprime)
|
||||
l2 = l*l
|
||||
half_l_di = 0.5*l*d_i
|
||||
h = h(t, tprime, d_i, d_j, l)
|
||||
arg_1 = half_l_di + tprime/l
|
||||
arg_2 = half_l_di - (t-tprime)/l
|
||||
ln_part_1 = ln_diff_erf(arg_1, arg_2)
|
||||
arg_1 = half_l_di
|
||||
arg_2 = half_l_di - t/l
|
||||
sign_val = sign(t/l)
|
||||
ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l)
|
||||
sign_val = sign(t/l)
|
||||
base = tprime*sign_val*exp(half_l_di*half_l_di-(d_i*t+d_j*tprime)+ln_part_2)-h
|
||||
return base/(d_i+d_j)
|
||||
|
||||
class dh_dl(Function):
|
||||
nargs = 5
|
||||
@classmethod
|
||||
def eval(cls, t, tprime, d_i, d_j, l):
|
||||
if (t.is_Number
|
||||
and tprime.is_Number
|
||||
and d_i.is_Number
|
||||
and d_j.is_Number
|
||||
and l.is_Number):
|
||||
|
||||
diff_t = (t-tprime)
|
||||
l2 = l*l
|
||||
h = h(t, tprime, d_i, d_j, l)
|
||||
return 0.5*d_i*d_i*l*h + 2./(sqrt(pi)*(d_i+d_j))*((-diff_t/l2-d_i/2.)*exp(-diff_t*diff_t/l2)+(-tprime/l2+d_i/2.)*exp(-tprime*tprime/l2-d_i*t)-(-t/l2-d_i/2.)*exp(-t*t/l2-d_j*tprime)-d_i/2.*exp(-(d_i*t+d_j*tprime)))
|
||||
|
||||
class dh_dt(Function):
|
||||
nargs = 5
|
||||
@classmethod
|
||||
def eval(cls, t, tprime, d_i, d_j, l):
|
||||
if (t.is_Number
|
||||
and tprime.is_Number
|
||||
and d_i.is_Number
|
||||
and d_j.is_Number
|
||||
and l.is_Number):
|
||||
if (t is S.NaN
|
||||
or tprime is S.NaN
|
||||
or d_i is S.NaN
|
||||
or d_j is S.NaN
|
||||
or l is S.NaN):
|
||||
return S.NaN
|
||||
else:
|
||||
half_l_di = 0.5*l*d_i
|
||||
arg_1 = half_l_di + tprime/l
|
||||
arg_2 = half_l_di - (t-tprime)/l
|
||||
ln_part_1 = ln_diff_erf(arg_1, arg_2)
|
||||
arg_1 = half_l_di
|
||||
arg_2 = half_l_di - t/l
|
||||
sign_val = sign(t/l)
|
||||
ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l)
|
||||
|
||||
|
||||
return (sign_val*exp(half_l_di*half_l_di
|
||||
- d_i*(t-tprime)
|
||||
+ ln_part_1
|
||||
- log(d_i + d_j))
|
||||
- sign_val*exp(half_l_di*half_l_di
|
||||
- d_i*t - d_j*tprime
|
||||
+ ln_part_2
|
||||
- log(d_i + d_j))).diff(t)
|
||||
|
||||
class dh_dtprime(Function):
|
||||
nargs = 5
|
||||
@classmethod
|
||||
def eval(cls, t, tprime, d_i, d_j, l):
|
||||
if (t.is_Number
|
||||
and tprime.is_Number
|
||||
and d_i.is_Number
|
||||
and d_j.is_Number
|
||||
and l.is_Number):
|
||||
if (t is S.NaN
|
||||
or tprime is S.NaN
|
||||
or d_i is S.NaN
|
||||
or d_j is S.NaN
|
||||
or l is S.NaN):
|
||||
return S.NaN
|
||||
else:
|
||||
half_l_di = 0.5*l*d_i
|
||||
arg_1 = half_l_di + tprime/l
|
||||
arg_2 = half_l_di - (t-tprime)/l
|
||||
ln_part_1 = ln_diff_erf(arg_1, arg_2)
|
||||
arg_1 = half_l_di
|
||||
arg_2 = half_l_di - t/l
|
||||
sign_val = sign(t/l)
|
||||
ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l)
|
||||
|
||||
|
||||
return (sign_val*exp(half_l_di*half_l_di
|
||||
- d_i*(t-tprime)
|
||||
+ ln_part_1
|
||||
- log(d_i + d_j))
|
||||
- sign_val*exp(half_l_di*half_l_di
|
||||
- d_i*t - d_j*tprime
|
||||
+ ln_part_2
|
||||
- log(d_i + d_j))).diff(tprime)
|
||||
|
||||
|
||||
class h(Function):
|
||||
nargs = 5
|
||||
def fdiff(self, argindex=5):
|
||||
t, tprime, d_i, d_j, l = self.args
|
||||
if argindex == 1:
|
||||
return dh_dt(t, tprime, d_i, d_j, l)
|
||||
elif argindex == 2:
|
||||
return dh_dtprime(t, tprime, d_i, d_j, l)
|
||||
elif argindex == 3:
|
||||
return dh_dd_i(t, tprime, d_i, d_j, l)
|
||||
elif argindex == 4:
|
||||
return dh_dd_j(t, tprime, d_i, d_j, l)
|
||||
elif argindex == 5:
|
||||
return dh_dl(t, tprime, d_i, d_j, l)
|
||||
|
||||
|
||||
@classmethod
|
||||
def eval(cls, t, tprime, d_i, d_j, l):
|
||||
if (t.is_Number
|
||||
and tprime.is_Number
|
||||
and d_i.is_Number
|
||||
and d_j.is_Number
|
||||
and l.is_Number):
|
||||
if (t is S.NaN
|
||||
or tprime is S.NaN
|
||||
or d_i is S.NaN
|
||||
or d_j is S.NaN
|
||||
or l is S.NaN):
|
||||
return S.NaN
|
||||
else:
|
||||
half_l_di = 0.5*l*d_i
|
||||
arg_1 = half_l_di + tprime/l
|
||||
arg_2 = half_l_di - (t-tprime)/l
|
||||
ln_part_1 = ln_diff_erf(arg_1, arg_2)
|
||||
arg_1 = half_l_di
|
||||
arg_2 = half_l_di - t/l
|
||||
sign_val = sign(t/l)
|
||||
ln_part_2 = ln_diff_erf(half_l_di, half_l_di - t/l)
|
||||
|
||||
|
||||
return (sign_val*exp(half_l_di*half_l_di
|
||||
- d_i*(t-tprime)
|
||||
+ ln_part_1
|
||||
- log(d_i + d_j))
|
||||
- sign_val*exp(half_l_di*half_l_di
|
||||
- d_i*t - d_j*tprime
|
||||
+ ln_part_2
|
||||
- log(d_i + d_j)))
|
||||
|
||||
|
||||
# return (exp((d_j/2.*l)**2)/(d_i+d_j)
|
||||
# *(exp(-d_j*(tprime - t))
|
||||
# *(erf((tprime-t)/l - d_j/2.*l)
|
||||
# + erf(t/l + d_j/2.*l))
|
||||
# - exp(-(d_j*tprime + d_i))
|
||||
# *(erf(tprime/l - d_j/2.*l)
|
||||
# + erf(d_j/2.*l))))
|
||||
|
||||
|
||||
|
||||
Loading…
Add table
Add a link
Reference in a new issue