def test_value(self):
with self.test_session(use_gpu=True):
def _test_value(logits, n_experiments, given):
logits = np.array(logits, np.float32)
normalized_logits = logits - misc.logsumexp(
logits, axis=-1, keepdims=True)
given = np.array(given)
dist = Multinomial(logits, n_experiments)
log_p = dist.log_prob(given)
target_log_p = np.log(misc.factorial(n_experiments)) - \
np.sum(np.log(misc.factorial(given)), -1) + \
np.sum(given * normalized_logits, -1)
self.assertAllClose(log_p.eval(), target_log_p)
p = dist.prob(given)
target_p = np.exp(target_log_p)
self.assertAllClose(p.eval(), target_p)
_test_value([-50., -20., 0.], 4, [1, 0, 3])
_test_value([1., 10., 1000.], 1, [1, 0, 0])
_test_value([[2., 3., 1.], [5., 7., 4.]], 3,
np.ones([3, 1, 3], dtype=np.int32))
_test_value([-10., 10., 20., 50.], 100, [[0, 1, 99, 100],
[100, 99, 1, 0]])
python类factorial()的实例源码
def power_series(z, t, C, spatial_der_order=0):
if not all([isinstance(item, (Number, np.ndarray)) for item in [z, t]]):
raise TypeError
z = np.atleast_1d(z)
t = np.atleast_1d(t)
if not all([len(item.shape) == 1 for item in [z, t]]):
raise ValueError
x = np.nan*np.zeros((len(t), len(z)))
for i in range(len(z)):
sum_x = np.zeros(t.shape[0])
for j in range(len(C)-spatial_der_order):
sum_x += C[j+spatial_der_order][0, :]*z[i]**j/sm.factorial(j)
x[:, i] = sum_x
if any([dim == 1 for dim in x.shape]):
x = x.flatten()
return x
def __init__(self, nsingle, indexing='Lin', symmetry='n', nleads=0):
"""
Initialization of the StateIndexingDM class
Parameters
----------
nsingle : int
Number of single-particle states.
indexing : str
String determining type of the indexing. Possible values are 'Lin', 'charge', 'sz', 'ssq'.
symmetry : str
String determining that the states will be augmented by the symmetry.
"""
StateIndexing.__init__(self, nsingle, indexing, symmetry, nleads)
self.ndm0_tot = int(factorial(2*self.nsingle)/factorial(self.nsingle)**2)
self.ndm1_tot = int(self.nsingle/(self.nsingle+1)*self.ndm0_tot)
#
self.shiftlst0 = np.zeros(self.ncharge+1, dtype=longnp)
self.shiftlst1 = np.zeros(self.ncharge, dtype=longnp)
self.lenlst = np.zeros(self.ncharge, dtype=longnp)
self.dictdm = np.zeros(self.nmany, dtype=longnp)
#
self.statesdm, self.mapdm0, self.mapdm1 = None, None, None
self.set_statesdm(self.chargelst)
def __init__(self, nsingle, indexing='Lin', symmetry='n', nleads=0):
"""
Initialization of the StateIndexingDMc class
Parameters
----------
nsingle : int
Number of single-particle states.
indexing : str
String determining type of the indexing. Possible values are 'Lin', 'charge', 'sz', 'ssq'.
symmetry : str
String determining that the states will be augmented by the symmetry.
"""
StateIndexing.__init__(self, nsingle, indexing, symmetry, nleads)
self.ndm0_tot = int(factorial(2*self.nsingle)/factorial(self.nsingle)**2)
self.ndm1_tot = int(self.nsingle/(self.nsingle+1)*self.ndm0_tot)
#
self.shiftlst0 = np.zeros(self.ncharge+1, dtype=longnp)
self.shiftlst1 = np.zeros(self.ncharge, dtype=longnp)
self.lenlst = np.zeros(self.ncharge, dtype=longnp)
self.dictdm = np.zeros(self.nmany, dtype=longnp)
#
self.statesdm, self.mapdm0, self.mapdm1 = None, None, None
self.set_statesdm(self.chargelst)
def test_log_combination(self):
with self.test_session(use_gpu=True):
def _test_func(n, ks):
tf_n = tf.convert_to_tensor(n, tf.float32)
tf_ks = tf.convert_to_tensor(ks, tf.float32)
true_value = np.log(misc.factorial(n)) - \
np.sum(np.log(misc.factorial(ks)), axis=-1)
test_value = log_combination(tf_n, tf_ks).eval()
self.assertAllClose(true_value, test_value)
_test_func(10, [1, 2, 3, 4])
_test_func([1, 2], [[1], [2]])
_test_func([1, 4], [[1, 0], [2, 2]])
_test_func([[2], [3]], [[[0, 2], [1, 2]]])
def __call__(self, t, d=0):
for m in range(len(self.T)):
if t > self.T[m] and m != len(self.T) - 1:
t -= self.T[m]
else: break
P = 0
for n in range(d, self.order + 1):
P += self.p[m][n] * (factorial(n) / factorial(n - d)) * t**(n - d)
return P
phonemeDurationStat.py 文件源码
项目：jingjuSingingPhraseMatching
作者: ronggong
项目源码
文件源码
阅读 18
收藏 0
点赞 0
评论 0
def poisson(k, lamb):
return (lamb**k/factorial(k)) * np.exp(-lamb)
def function(x, lamb, offsX, offsY, scaleX, scaleY):
# poisson function, parameter lamb is the fit parameter
x = (x - offsX) / scaleX
y = (lamb**x / factorial(x)) * np.exp(-lamb)
return scaleY * y + offsY
def jk(D,Bs,i,j,frac,possible,ri):
Num=factorial(ri-1)
Den=factorial(N_ij(D,Bs,i,j,ri)+ri-1)
frac*=np.float128(Num)/np.float128(Den)
for k in range(0,ri):
frac*=factorial(N_ijk(D,Bs,i,j,k,ri))
if N_ijk(D,Bs,i,j,k,ri)!=0 :
possible=1
return frac,possible
def __init__(self, states, interval, method, differential_order=0):
"""
:param states: tuple of states in beginning and end of interval
:param interval: time interval (tuple)
:param method: method to use (``poly`` or ``tanh``)
:param differential_order: grade of differential flatness :math:`\\gamma`
"""
self.yd = states
self.t0 = interval[0]
self.t1 = interval[1]
self.dt = interval[1] - interval[0]
# setup symbolic expressions
if method == "tanh":
tau, sigma = sp.symbols('tau, sigma')
# use a gevrey-order of alpha = 1 + 1/sigma
sigma = 1.1
phi = .5*(1 + sp.tanh(2*(2*tau - 1)/((4*tau*(1-tau))**sigma)))
elif method == "poly":
gamma = differential_order # + 1 # TODO check this against notes
tau, k = sp.symbols('tau, k')
alpha = sp.factorial(2 * gamma + 1)
f = sp.binomial(gamma, k) * (-1) ** k * tau ** (gamma + k + 1) / (gamma + k + 1)
phi = alpha / sp.factorial(gamma) ** 2 * sp.summation(f, (k, 0, gamma))
else:
raise NotImplementedError("method {} not implemented!".format(method))
# differentiate phi(tau), index in list corresponds to order
dphi_sym = [phi] # init with phi(tau)
for order in range(differential_order):
dphi_sym.append(dphi_sym[-1].diff(tau))
# lambdify
self.dphi_num = []
for der in dphi_sym:
self.dphi_num.append(sp.lambdify(tau, der, 'numpy'))
def temporal_derived_power_series(z, C, up_to_order, series_termination_index, spatial_der_order=0):
"""
compute the temporal derivatives
q^{(n)}(z) = \sum_{k=0}^{series_termination_index} C[k][n,:] z^k / k!
from n=0 to n=up_to_order
:param z: scalar
:param C:
:param up_to_order:
:param series_termination_index:
:param spatial_der_order:
:return: Q = np.array( [q^{(0)}, ... , q^{(up_to_order)}] )
"""
if not isinstance(z, Number):
raise TypeError
if any([C[i].shape[0] - 1 < up_to_order for i in range(series_termination_index+1)]):
raise ValueError
len_t = C[0].shape[1]
Q = np.nan*np.zeros((up_to_order+1, len_t))
for i in range(up_to_order+1):
sum_Q = np.zeros(len_t)
for j in range(series_termination_index+1-spatial_der_order):
sum_Q += C[j+spatial_der_order][i, :]*z**(j)/sm.factorial(j)
Q[i, :] = sum_Q
return Q
def _pade_delay(p, q, c):
"""Numerically evaluated state-space using Pade approximants.
This may have numerical issues for large values of p or q.
"""
i = np.arange(1, p+q+1, dtype=np.float64)
taylor = np.append([1.0], (-c)**i / factorial(i))
num, den = pade(taylor, q)
return LinearSystem((num, den))
def time(self, t, s=1.0):
"""
Complex Paul wavelet, centred at zero.
Parameters
----------
t : float
Time. If s is not specified, i.e. set to 1, this can be
used as the non-dimensional time t/s.
s : float
Scaling factor. Default is 1.
Returns
-------
complex: value of the paul wavelet at the given time
The Paul wavelet is defined (in time) as::
(2 ** m * i ** m * m!) / (pi * (2 * m)!) \
* (1 - i * t / s) ** -(m + 1)
"""
m = self.m
x = t / s
const = (2 ** m * 1j ** m * factorial(m)) \
/ (np.pi * factorial(2 * m)) ** .5
functional_form = (1 - 1j * x) ** -(m + 1)
output = const * functional_form
return output
# Fourier wavelengths
def frequency(self, w, s=1.0):
"""Frequency representation of Paul.
Parameters
----------
w : float
Angular frequency. If s is not specified, i.e. set to 1,
this can be used as the non-dimensional angular
frequency w * s.
s : float
Scaling factor. Default is 1.
Returns
-------
complex: value of the paul wavelet at the given time
"""
m = self.m
x = w * s
# heaviside mock
Hw = 0.5 * (np.sign(x) + 1)
# prefactor
const = 2 ** m / (m * factorial(2 * m - 1)) ** .5
functional_form = Hw * (x) ** m * np.exp(-x)
output = const * functional_form
return output
def time(self, t, s=1.0):
"""
Complex Paul wavelet, centred at zero.
Parameters
----------
t : float
Time. If s is not specified, i.e. set to 1, this can be
used as the non-dimensional time t/s.
s : float
Scaling factor. Default is 1.
Returns
-------
complex: value of the paul wavelet at the given time
The Paul wavelet is defined (in time) as::
(2 ** m * i ** m * m!) / (pi * (2 * m)!) \
* (1 - i * t / s) ** -(m + 1)
"""
m = self.m
x = t / s
const = (2 ** m * 1j ** m * factorial(m)) \
/ (np.pi * factorial(2 * m)) ** .5
functional_form = (1 - 1j * x) ** -(m + 1)
output = const * functional_form
return output
# Fourier wavelengths
def frequency(self, w, s=1.0):
"""Frequency representation of Paul.
Parameters
----------
w : float
Angular frequency. If s is not specified, i.e. set to 1,
this can be used as the non-dimensional angular
frequency w * s.
s : float
Scaling factor. Default is 1.
Returns
-------
complex: value of the paul wavelet at the given time
"""
m = self.m
x = w * s
# heaviside mock
Hw = 0.5 * (np.sign(x) + 1)
# prefactor
const = 2 ** m / (m * factorial(2 * m - 1)) ** .5
functional_form = Hw * (x) ** m * np.exp(-x)
output = const * functional_form
return output
def rnm(n,m,rho):
"""
Return an array with the zernike Rnm polynomial calculated at rho points.
"""
Rnm=zeros(rho.shape)
S=(n-abs(m))/2
for s in range (0,S+1):
CR=pow(-1,s)*factorial(n-s)/ \
(factorial(s)*factorial(-s+(n+abs(m))/2)* \
factorial(-s+(n-abs(m))/2))
p=CR*pow(rho,n-2*s)
Rnm=Rnm+p
Rnm[rho>1.0] = 0
return Rnm
def time(self, t, s=1.0):
"""
Complex Paul wavelet, centred at zero.
Parameters
----------
t : float
Time. If `s` is not specified, i.e. set to 1, this can be
used as the non-dimensional time t/s.
s : float
Scaling factor. Default is 1.
Returns
-------
out : complex
Value of the Paul wavelet at the given time
The Paul wavelet is defined (in time) as::
(2 ** m * i ** m * m!) / (pi * (2 * m)!) \
* (1 - i * t / s) ** -(m + 1)
"""
m = self.m
x = t / s
const = (2 ** m * 1j ** m * factorial(m)) \
/ (np.pi * factorial(2 * m)) ** .5
functional_form = (1 - 1j * x) ** -(m + 1)
output = const * functional_form
return output
# Fourier wavelengths
def frequency(self, w, s=1.0):
"""Frequency representation of Paul.
Parameters
----------
w : float
Angular frequency. If `s` is not specified, i.e. set to 1,
this can be used as the non-dimensional angular
frequency w * s.
s : float
Scaling factor. Default is 1.
Returns
-------
out : complex
Value of the Paul wavelet at the given frequency
"""
m = self.m
x = w * s
# Heaviside mock
Hw = 0.5 * (np.sign(x) + 1)
# prefactor
const = 2 ** m / (m * factorial(2 * m - 1)) ** .5
functional_form = Hw * (x) ** m * np.exp(-x)
output = const * functional_form
return output
def gamma(n, tau, tsample, tol=0.01):
"""Returns the impulse response of `n` cascaded leaky integrators
This function calculates the impulse response of `n` cascaded
leaky integrators with constant of proportionality 1/`tau`:
y = (t/theta).^(n-1).*exp(-t/theta)/(theta*factorial(n-1))
Parameters
----------
n : int
Number of cascaded leaky integrators
tau : float
Decay constant of leaky integration (seconds).
Equivalent to the inverse of the constant of proportionality.
tsample : float
Sampling time step (seconds).
tol : float
Cut the kernel to size by ignoring function values smaller
than a fraction `tol` of the peak value.
"""
n = int(n)
tau = float(tau)
tsample = float(tsample)
if n <= 0 or tau <= 0 or tsample <= 0:
raise ValueError("`n`, `tau`, and `tsample` must be nonnegative.")
if tau <= tsample:
raise ValueError("`tau` cannot be smaller than `tsample`.")
# Allocate a time vector that is long enough for sure.
# Trim vector later on.
t = np.arange(0, 5 * n * tau, tsample)
# Calculate gamma
y = (t / tau) ** (n - 1) * np.exp(-t / tau)
y /= (tau * spm.factorial(n - 1))
# Normalize to unit area
y /= np.trapz(np.abs(y), dx=tsample)
# Cut off tail where values are smaller than `tol`.
# Make sure to start search on the right-hand side of the peak.
peak = y.argmax()
small_vals = np.where(y[peak:] < tol * y.max())[0]
if small_vals.size:
t = t[:small_vals[0] + peak]
y = y[:small_vals[0] + peak]
return t, y
def _power_series_flat_out(z, t, n, param, y, bound_cond_type):
""" Provide the power series approximation for x(z,t) and x'(z,t).
:param z: [0, ..., l] (numpy array)
:param t: [0, ... , t_end] (numpy array)
:param n: series termination index (integer)
:param param: [a2, a1, a0, alpha, beta] (list)
:param y: flat output with derivation np.array([[y],...,[y^(n/2)]])
:return: field variable x(z,t) and spatial derivation x'(z,t)
"""
# TODO: documentation
a2, a1, a0, alpha, beta = param
shape = (len(t), len(z))
x = np.nan*np.ones(shape)
d_x = np.nan*np.ones(shape)
# Actually power_series() is designed for robin boundary condition by z=0.
# With the following modification it can also used for dirichlet boundary condition by z=0.
if bound_cond_type is 'robin':
is_robin = 1.
elif bound_cond_type is 'dirichlet':
alpha = 1.
is_robin = 0.
else:
raise ValueError("Selected Boundary condition {0} not supported! Use 'robin' or 'dirichlet'".format(
bound_cond_type))
# TODO: flip iteration order: z <--> t, result: one or two instead len(t) call's
for i in range(len(t)):
sum_x = np.zeros(len(z))
for j in range(n):
sum_b = np.zeros(len(z))
for k in range(j+1):
sum_b += sm.comb(j, k)*(-a0)**(j-k)*y[k, i]
sum_x += (is_robin+alpha*z/(2.*j+1.))*z**(2*j)/sm.factorial(2*j)/a2**j*sum_b
x[i, :] = sum_x
for i in range(len(t)):
sum_x = np.zeros(len(z))
for j in range(n):
sum_b = np.zeros(len(z))
for k in range(j+2):
sum_b += sm.comb(j+1, k)*(-a0)**(j-k+1)*y[k, i]
if j == 0:
sum_x += alpha*y[0, i]
sum_x += (is_robin+alpha*z/(2.*(j+1)))*z**(2*j+1)/sm.factorial(2*j+1)/a2**(j+1)*sum_b
d_x[i, :] = sum_x
return x, d_x
