我们从Python开源项目中，提取了以下21个代码示例，用于说明如何使用scipy.misc.factorial()。
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]])
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
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 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
