apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-15 06:52:39 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'devel' of github.com:SheffieldML/GPy into devel

Browse source
This commit is contained in:
Max Zwiessele 2014-10-16 12:43:19 +01:00
commit 3358d06e42
17 changed files with 184 additions and 522 deletions
Download patch file Download diff file Expand all files Collapse all files

10
GPy/core/model.py
View file

@ -381,6 +381,16 @@ class Model(Parameterized):
self.optimizer_array = x
return ret
def _repr_html_(self):
"""Representation of the model in html for notebook display."""
model_details = [['<b>Model</b>', self.name + '<br>'],
['<b>Log-likelihood</b>', '{}<br>'.format(float(self.log_likelihood()))],
["<b>Number of Parameters</b>", '{}<br>'.format(self.size)]]
from operator import itemgetter
to_print = [""] + ["{}: {}".format(name, detail) for name, detail in model_details] + ["<br><b>Parameters</b>:"]
to_print.append(super(Model, self)._repr_html_())
return "\n".join(to_print)
def __str__(self):
model_details = [['Name', self.name],
['Log-likelihood', '{}'.format(float(self.log_likelihood()))],

23
GPy/core/parameterization/param.py
View file

@ -251,6 +251,29 @@ class Param(Parameterizable, ObsAr):
if ind.size > 4: indstr = ','.join(map(str, ind[:2])) + "..." + ','.join(map(str, ind[-2:]))
else: indstr = ','.join(map(str, ind))
return name + '[' + indstr + ']'
def _repr_html_(self, constr_matrix=None, indices=None, prirs=None, ties=None):
"""Representation of the parameter in html for notebook display."""
filter_ = self._current_slice_
vals = self.flat
if indices is None: indices = self._indices(filter_)
ravi = self._raveled_index(filter_)
if constr_matrix is None: constr_matrix = self.constraints.properties_for(ravi)
if prirs is None: prirs = self.priors.properties_for(ravi)
if ties is None: ties = self._ties_for(ravi)
ties = [' '.join(map(lambda x: x, t)) for t in ties]
header_format = """
<tr>
<td><b>{i}</b></td>
<td><b>{x}</b></td>
<td><b>{c}</b></td>
<td><b>{p}</b></td>
<td><b>{t}</b></td>
</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 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_
vals = self.flat

29
GPy/core/parameterization/parameterized.py
View file

@ -357,6 +357,35 @@ class Parameterized(Parameterizable):
@property
def _ties_str(self):
return [','.join(x._ties_str) for x in self.flattened_parameters]
def _repr_html_(self, header=True):
"""Representation of the parameters in html for notebook display."""
name = adjust_name_for_printing(self.name) + "."
constrs = self._constraints_str;
ts = self._ties_str
prirs = self._priors_str
desc = self._description_str; names = self.parameter_names()
nl = max([len(str(x)) for x in names + [name]])
sl = max([len(str(x)) for x in desc + ["Value"]])
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 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))
sep = '-' * (nl + sl + cl + + pl + tl + 8 * 2 + 3)
if header:
header = """
<tr>
<td><b>{name}</b>
<td><b>Value</b></td>
<td><b>Constraint</b></td>
<td><b>Prior</b></td>
<td><b>Tied to</b></td>""".format(name=name)
to_print.insert(0, header)
return '<table>' + '\n'.format(sep).join(to_print) + '\n</table>'
def __str__(self, header=True):
name = adjust_name_for_printing(self.name) + "."
constrs = self._constraints_str;

12
GPy/core/parameterization/priors.py
View file

@ -52,7 +52,7 @@ class Gaussian(Prior):
self.constant = -0.5 * np.log(2 * np.pi * self.sigma2)
def __str__(self):
return "N(" + str(np.round(self.mu)) + ', ' + str(np.round(self.sigma2)) + ')'
return "N({:.2g}, {:.2g})".format(self.mu, self.sigma)
def lnpdf(self, x):
return self.constant - 0.5 * np.square(x - self.mu) / self.sigma2
@ -82,7 +82,7 @@ class Uniform(Prior):
self.upper = float(upper)
def __str__(self):
return "[" + str(np.round(self.lower)) + ', ' + str(np.round(self.upper)) + ']'
return "[{:.2g}, {:.2g}]".format(self.lower, self.upper)
def lnpdf(self, x):
region = (x>=self.lower) * (x<=self.upper)
@ -122,7 +122,7 @@ class LogGaussian(Prior):
self.constant = -0.5 * np.log(2 * np.pi * self.sigma2)
def __str__(self):
return "lnN(" + str(np.round(self.mu)) + ', ' + str(np.round(self.sigma2)) + ')'
return "lnN({:.2g}, {:.2g})".format(self.mu, self.sigma)
def lnpdf(self, x):
return self.constant - 0.5 * np.square(np.log(x) - self.mu) / self.sigma2 - np.log(x)
@ -220,7 +220,7 @@ class Gamma(Prior):
self.constant = -gammaln(self.a) + a * np.log(b)
def __str__(self):
return "Ga(" + str(np.round(self.a)) + ', ' + str(np.round(self.b)) + ')'
return "Ga({:.2g}, {:.2g})".format(self.a, self.b)
def summary(self):
ret = {"E[x]": self.a / self.b, \
@ -281,7 +281,7 @@ class InverseGamma(Prior):
self.constant = -gammaln(self.a) + a * np.log(b)
def __str__(self):
return "iGa(" + str(np.round(self.a)) + ', ' + str(np.round(self.b)) + ')'
return "iGa({:.2g}, {:.2g})".format(self.a, self.b)
def lnpdf(self, x):
return self.constant - (self.a + 1) * np.log(x) - self.b / x
@ -317,7 +317,7 @@ class HalfT(Prior):
self.constant = gammaln(.5*(self.nu+1.)) - gammaln(.5*self.nu) - .5*np.log(np.pi*self.A*self.nu)
def __str__(self):
return "hT(" + str(np.round(self.A)) + ', ' + str(np.round(self.nu)) + ')'
return "hT({:.2g}, {:.2g})".format(self.A, self.nu)
def lnpdf(self,theta):
return (theta>0) * ( self.constant -.5*(self.nu+1) * np.log( 1.+ (1./self.nu) * (theta/self.A)**2 ) )

2
GPy/core/symbolic.py
View file

@ -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)

4
GPy/inference/latent_function_inference/laplace.py
View file

@ -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

30
GPy/kern/_src/eq.py
View file

@ -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)

30
GPy/kern/_src/heat_eqinit.py
View file

@ -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)

74
GPy/kern/_src/stationary.py
View file

@ -10,6 +10,7 @@ from ... import util
import numpy as np
from scipy import integrate
from ...util.caching import Cache_this
from ...util.config import config # for assesing whether to use weave
class Stationary(Kern):
"""
@ -146,7 +147,17 @@ class Stationary(Kern):
#self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0)/self.lengthscale**3
tmp = dL_dr*self._inv_dist(X, X2)
if X2 is None: X2 = X
self.lengthscale.gradient = np.array([np.einsum('ij,ij,...', tmp, np.square(X[:,q:q+1] - X2[:,q:q+1].T), -1./self.lengthscale[q]**3) for q in xrange(self.input_dim)])
if config.getboolean('weave', 'working'):
try:
self.lengthscale.gradient = self.weave_lengthscale_grads(tmp, X, X2)
except:
print "\n Weave compilation failed. Falling back to (slower) numpy implementation\n"
config.set('weave', 'working', 'False')
self.lengthscale.gradient = np.array([np.einsum('ij,ij,...', tmp, np.square(X[:,q:q+1] - X2[:,q:q+1].T), -1./self.lengthscale[q]**3) for q in xrange(self.input_dim)])
else:
self.lengthscale.gradient = np.array([np.einsum('ij,ij,...', tmp, np.square(X[:,q:q+1] - X2[:,q:q+1].T), -1./self.lengthscale[q]**3) for q in xrange(self.input_dim)])
else:
r = self._scaled_dist(X, X2)
self.lengthscale.gradient = -np.sum(dL_dr*r)/self.lengthscale
@ -161,10 +172,43 @@ class Stationary(Kern):
dist = self._scaled_dist(X, X2).copy()
return 1./np.where(dist != 0., dist, np.inf)
def weave_lengthscale_grads(self, tmp, X, X2):
N,M = tmp.shape
Q = X.shape[1]
if hasattr(X, 'values'):X = X.values
if hasattr(X2, 'values'):X2 = X2.values
grads = np.zeros(self.input_dim)
code = """
double gradq;
for(int q=0; q<Q; q++){
gradq = 0;
for(int n=0; n<N; n++){
for(int m=0; m<M; m++){
gradq += tmp(n,m)*(X(n,q)-X2(m,q))*(X(n,q)-X2(m,q));
}
}
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
def gradients_X(self, dL_dK, X, X2=None):
"""
Given the derivative of the objective wrt K (dL_dK), compute the derivative wrt X
"""
if config.getboolean('weave', 'working'):
try:
return self.gradients_X_weave(dL_dK, X, X2)
except:
print "\n Weave compilation failed. Falling back to (slower) numpy implementation\n"
config.set('weave', 'working', 'False')
return self.gradients_X_(dL_dK, X, X2)
else:
return self.gradients_X_(dL_dK, X, X2)
def gradients_X_(self, dL_dK, X, X2=None):
invdist = self._inv_dist(X, X2)
dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
tmp = invdist*dL_dr
@ -184,6 +228,34 @@ class Stationary(Kern):
return ret
def gradients_X_weave(self, dL_dK, X, X2=None):
invdist = self._inv_dist(X, X2)
dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
tmp = invdist*dL_dr
if X2 is None:
tmp = tmp + tmp.T
X2 = X
ret = np.zeros(X.shape)
code = """
int n,q,d;
double retnd;
for(n=0;n<N;n++){
for(d=0;d<D;d++){
retnd = 0;
for(q=0;q<Q;q++){
retnd += tmp(n,q)*(X(n,d)-X2(q,d));
}
ret(n,d) = retnd;
}
}
"""
from scipy import weave
N,D = X.shape
Q = tmp.shape[1]
weave.inline(code, ['ret', 'N', 'D', 'Q', 'tmp', 'X', 'X2'], type_converters=weave.converters.blitz)
return ret/self.lengthscale**2
def gradients_X_diag(self, dL_dKdiag, X):
return np.zeros(X.shape)

6
GPy/likelihoods/bernoulli.py
View file

@ -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):
"""

34
GPy/likelihoods/likelihood.py
View file

@ -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) )

1
GPy/likelihoods/poisson.py
View file

@ -5,7 +5,6 @@ from __future__ import division
import numpy as np
from scipy import stats,special
import scipy as sp
from GPy.util.univariate_Gaussian import std_norm_pdf,std_norm_cdf
import link_functions
from likelihood import Likelihood

2
GPy/likelihoods/skew_exponential.py
View file

@ -5,8 +5,6 @@
import sympy as sym
from GPy.util.symbolic import normcdfln
import numpy as np
from ..util.univariate_Gaussian import std_norm_pdf, std_norm_cdf
#from GPy.util.functions import clip_exp
import link_functions
from symbolic import Symbolic
from scipy import stats

2
GPy/testing/likelihood_tests.py
View file

@ -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)

5
GPy/util/__init__.py
View file

@ -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

75
GPy/util/datasets/swiss_roll.pickle
View file

File diff suppressed because one or more lines are too long

367
GPy/util/symbolic.py
View file

@ -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))))