def rvs( self ):
if self.ss == 0:
x = random.random() * self.w
y = random.random() * self.h
self.set_point( x, y )
while self.qs:
x_idx = int( random.random() * self.qs )
s = self.queue[ x_idx ]
for y_idx in range( self.n ):
a = 2 * scipy.pi * random.random()
b = scipy.sqrt( self.A * random.random() + self.r2 )
x = s[0] + b*scipy.cos( a )
y = s[1] + b*scipy.sin( a )
if( x >= 0 )and( x < self.w ):
if( y >= 0 )and( y < self.h ):
if( self.distance( x, y ) ):
self.set_point( x, y )
del self.queue[x_idx]
self.qs -= 1
sample = list( filter( None, self.grid ) )
sample = scipy.asfarray( sample )
return sample
python类pi()的实例源码
poisson_disk_sampling.py 文件源码
项目:nodebox1-generative-tools
作者: x-raizor
项目源码
文件源码
阅读 41
收藏 0
点赞 0
评论 0
def fit(self, X, w):
if len(X) == 0:
raise NotEnoughParticles("Fitting not possible.")
self.X_arr = X.as_matrix()
ctree = cKDTree(X)
_, indices = ctree.query(X, k=min(self.k + 1, X.shape[0]))
covs, inv_covs, dets = list(zip(*[self._cov_and_inv(n, indices)
for n in range(X.shape[0])]))
self.covs = sp.array(covs)
self.inv_covs = sp.array(inv_covs)
self.determinants = sp.array(dets)
self.normalization = sp.sqrt(
(2 * sp.pi) ** self.X_arr.shape[1] * self.determinants)
if not sp.isreal(self.normalization).all():
raise Exception("Normalization not real")
self.normalization = sp.real(self.normalization)
def model(THETA, time, kplanets):
modelo = 0.0
if kplanets == 0:
return 0.0
for i in range(kplanets):
As, P, Ac, S, C = THETA[5*i:5*(i+1)]
A = As ** 2 + Ac ** 2
ecc = S ** 2 + C ** 2
w = sp.arccos(C / (ecc ** 0.5)) # longitude of periastron
phase = sp.arccos(Ac / (A ** 0.5))
### test
if S < 0:
w = 2 * sp.pi - sp.arccos(C / (ecc ** 0.5))
if As < 0:
phase = 2 * sp.pi - sp.arccos(Ac / (A ** 0.5))
###
per = sp.exp(P)
freq = 2. * sp.pi / per
M = freq * time + phase
E = sp.array([MarkleyKESolver().getE(m, ecc) for m in M])
f = (sp.arctan(((1. + ecc) ** 0.5 / (1. - ecc) ** 0.5) * sp.tan(E / 2.)) * 2.)
modelo += A * (sp.cos(f + w) + ecc * sp.cos(w))
return modelo
def kde(data, N=None, MIN=None, MAX=None):
# Parameters to set up the mesh on which to calculate
N = 2**14 if N is None else int(2**sci.ceil(sci.log2(N)))
if MIN is None or MAX is None:
minimum = min(data)
maximum = max(data)
Range = maximum - minimum
MIN = minimum - Range/10 if MIN is None else MIN
MAX = maximum + Range/10 if MAX is None else MAX
# Range of the data
R = MAX-MIN
# Histogram the data to get a crude first approximation of the density
M = len(data)
DataHist, bins = sci.histogram(data, bins=N, range=(MIN,MAX))
DataHist = DataHist/M
DCTData = scipy.fftpack.dct(DataHist, norm=None)
I = [iN*iN for iN in xrange(1, N)]
SqDCTData = (DCTData[1:]/2)**2
# The fixed point calculation finds the bandwidth = t_star
guess = 0.1
try:
t_star = scipy.optimize.brentq(fixed_point, 0, guess,
args=(M, I, SqDCTData))
except ValueError:
print 'Oops!'
return None
# Smooth the DCTransformed data using t_star
SmDCTData = DCTData*sci.exp(-sci.arange(N)**2*sci.pi**2*t_star/2)
# Inverse DCT to get density
density = scipy.fftpack.idct(SmDCTData, norm=None)*N/R
mesh = [(bins[i]+bins[i+1])/2 for i in xrange(N)]
bandwidth = sci.sqrt(t_star)*R
density = density/sci.trapz(density, mesh)
return bandwidth, mesh, density
def fixed_point(t, M, I, a2):
l=7
I = sci.float128(I)
M = sci.float128(M)
a2 = sci.float128(a2)
f = 2*sci.pi**(2*l)*sci.sum(I**l*a2*sci.exp(-I*sci.pi**2*t))
for s in range(l, 1, -1):
K0 = sci.prod(xrange(1, 2*s, 2))/sci.sqrt(2*sci.pi)
const = (1 + (1/2)**(s + 1/2))/3
time=(2*const*K0/M/f)**(2/(3+2*s))
f=2*sci.pi**(2*s)*sci.sum(I**s*a2*sci.exp(-I*sci.pi**2*time))
return t-(2*M*sci.sqrt(sci.pi)*f)**(-2/5)
def nLLeval(ldelta, Uy, S, REML=True):
"""
evaluate the negative log likelihood of a random effects model:
nLL = 1/2(n_s*log(2pi) + logdet(K) + 1/ss * y^T(K + deltaI)^{-1}y,
where K = USU^T.
Uy: transformed outcome: n_s x 1
S: eigenvectors of K: n_s
ldelta: log-transformed ratio sigma_gg/sigma_ee
"""
n_s = Uy.shape[0]
delta = scipy.exp(ldelta)
# evaluate log determinant
Sd = S + delta
ldet = scipy.sum(scipy.log(Sd))
# evaluate the variance
Sdi = 1.0 / Sd
Sdi=Sdi.reshape((Sdi.shape[0],1))
ss = 1. / n_s * (Uy*Uy*Sdi).sum()
# evalue the negative log likelihood
nLL = 0.5 * (n_s * scipy.log(2.0 * scipy.pi) + ldet + n_s + n_s * scipy.log(ss))
if REML:
pass
return nLL
def ellipse2bbox(a, b, angle, cx, cy):
a, b = max(a, b), min(a, b)
ca = sp.cos(angle)
sa = sp.sin(angle)
if sa == 0.0:
cta = 2.0 / sp.pi
else:
cta = ca / sa
if ca == 0.0:
ta = sp.pi / 2.0
else:
ta = sa / ca
x = lambda t: cx + a * sp.cos(t) * ca - b * sp.sin(t) * sa
y = lambda t: cy + b * sp.sin(t) * ca + a * sp.cos(t) * sa
# x = cx + a * cos(t) * cos(angle) - b * sin(t) * sin(angle)
# tan(t) = -b * tan(angle) / a
tx1 = sp.arctan(-b * ta / a)
tx2 = tx1 - sp.pi
x1, y1 = x(tx1), y(tx1)
x2, y2 = x(tx2), y(tx2)
# y = cy + b * sin(t) * cos(angle) + a * cos(t) * sin(angle)
# tan(t) = b * cot(angle) / a
ty1 = sp.arctan(b * cta / a)
ty2 = ty1 - sp.pi
x3, y3 = x(ty1), y(ty1)
x4, y4 = x(ty2), y(ty2)
minx, maxx = Util.minmax([x1, x2, x3, x4])
miny, maxy = Util.minmax([y1, y2, y3, y4])
return sp.floor(minx), sp.floor(miny), sp.ceil(maxx), sp.ceil(maxy)
multiple_independent_Utubes.py 文件源码
项目:pygfunction
作者: MassimoCimmino
项目源码
文件源码
阅读 30
收藏 0
点赞 0
评论 0
def _pipePositions(Ds, nPipes):
""" Positions pipes in an axisymetric configuration.
"""
dt = pi / float(nPipes)
pos = [(0., 0.) for i in range(2*nPipes)]
for i in range(nPipes):
pos[i] = (Ds*np.cos(2.0*i*dt+pi), Ds*np.sin(2.0*i*dt+pi))
pos[i+nPipes] = (Ds*np.cos(2.0*i*dt+pi+dt), Ds*np.sin(2.0*i*dt+pi+dt))
return pos
# Main function
def semimodel(self, params, time):
A, P, phase, w, ecc = params
freq = 2. * sp.pi / P
M = freq * time + phase
E = sp.array([MarkleyKESolver().getE(m, ecc) for m in M])
f = (sp.arctan(((1. + ecc) ** 0.5 / (1. - ecc) ** 0.5) * sp.tan(E / 2.)) * 2.)
modelo = A * (sp.cos(f + w) + ecc * sp.cos(w))
return modelo
def henshin(self, thetas, kplanets):
for i in range(kplanets):
Ask = thetas[:, i*5]
Pk = thetas[:, i*5 + 1]
Ack = thetas[:, i*5 + 2]
Sk = thetas[:, i*5 + 3]
Ck = thetas[:, i*5 + 4]
ecck = Sk ** 2 + Ck ** 2
Ak = Ask ** 2 + Ack ** 2
wk = sp.arccos(Ck / (ecck ** 0.5))
Phasek = sp.arccos(Ack / (Ak ** 0.5))
for j in range(len(Sk)):
if Sk[j] < 0:
wk[j] = 2 * sp.pi - sp.arccos(Ck[j] / (ecck[j] ** 0.5))
if Ask[j] < 0:
Phasek[j] = 2 * sp.pi - sp.arccos(Ack[j] / (Ak[j] ** 0.5))
thetas[:, i*5] = Ak
thetas[:, i*5 + 1] = sp.exp(Pk)
thetas[:, i*5 + 2] = Phasek
thetas[:, i*5 + 3] = wk
thetas[:, i*5 + 4] = ecck
return thetas
########################################
# los IC sirven para los mod_lims? NO!
# mod_lims sólo acotan el espacio donde buscar, excepto para periodo
def nLLeval(ldelta, Uy, S, REML=True):
"""
evaluate the negative log likelihood of a random effects model:
nLL = 1/2(n_s*log(2pi) + logdet(K) + 1/ss * y^T(K + deltaI)^{-1}y,
where K = USU^T.
Uy: transformed outcome: n_s x 1
S: eigenvectors of K: n_s
ldelta: log-transformed ratio sigma_gg/sigma_ee
"""
n_s = Uy.shape[0]
delta = scipy.exp(ldelta)
# evaluate log determinant
Sd = S + delta
ldet = scipy.sum(scipy.log(Sd))
# evaluate the variance
Sdi = 1.0 / Sd
Uy = Uy.flatten()
ss = 1. / n_s * (Uy * Uy * Sdi).sum()
# evalue the negative log likelihood
nLL = 0.5 * (n_s * scipy.log(2.0 * scipy.pi) + ldet + n_s + n_s * scipy.log(ss))
if REML:
pass
return nLL
def func_2vN(Ek, Ek_grid, l, eta, hfk):
"""
Linearly interpolate the value of hfk on Ek_grid at point Ek.
Parameters
----------
Ek : float
Energy value (not necessarily a grid point).
Ek_grid : array
Energy grid.
l : int
Lead label.
eta : int
Integer describing (+/-1) if infinitesimal eta is positive or negative.
hfk : array
Array containing Hilbert transform of Fermi function (or 1-Fermi).
Returns
-------
float
Interpolated value of hfk at Ek.
"""
if Ek<Ek_grid[0] or Ek>Ek_grid[-1]:
return 0
#
b_idx = int((Ek-Ek_grid[0])/(Ek_grid[1]-Ek_grid[0]))+1
if b_idx == len(Ek_grid): b_idx -= 1
a_idx = b_idx - 1
b, a = Ek_grid[b_idx], Ek_grid[a_idx]
#
fb = hfk[l, b_idx]
fa = hfk[l, a_idx]
rez = (fb-fa)/(b-a)*Ek + (b*fa-a*fb)/(b-a)
return pi*rez if eta+1 else pi*rez.conjugate()
def func_pauli(Ecb, mu, T, Dm, Dp, itype):
"""
Function used when generating Pauli master equation kernel.
Parameters
----------
Ecb : float
Energy.
mu : float
Chemical potential.
T : float
Temperature.
Dm,Dp : float
Bandwidth.
Returns
-------
array
| Array of two float numbers [cur0, cur1] containing
momentum-integrated current amplitudes.
| cur0 - particle current amplitude.
| cur1 - hole current amplitude.
"""
alpha = (Ecb-mu)/T
Rm, Rp = (Dm-mu)/T, (Dp-mu)/T
if itype == 1 or itype == 3 or (alpha < Rp and alpha > Rm):
cur0 = fermi_func(alpha)
cur1 = 1-cur0
rez = 2*pi*np.array([cur0, cur1])
else:
rez = np.zeros(2)
return rez
def logl(theta, time, rv, err, ins, staract, starflag, kplanets, nins, MOAV, totcornum):
i, lnl = 0, 0
ndat = len(time)
model_params = kplanets * 5
ins_params = nins * 2 * (MOAV + 1)
jitter, offset, macoef, timescale = sp.zeros(ndat), sp.zeros(ndat), sp.array([sp.zeros(ndat) for i in range(MOAV)]), sp.array([sp.zeros(ndat) for i in range(MOAV)])
ACC = theta[model_params] * (time - time[0])
residuals = sp.zeros(ndat)
for i in range(ndat):
jitpos = int(model_params + ins[i] * 2 * (MOAV+1) + 1)
jitter[i], offset[i] = theta[jitpos], theta[jitpos + 1] # jitt
for j in range(MOAV):
macoef[j][i], timescale[j][i] = theta[jitpos + 2*(j+1)], theta[jitpos + 2*(j+1) + 1]
a1 = (theta[:model_params])
if totcornum:
COR = sp.array([sp.array([sp.zeros(ndat) for k in range(len(starflag[i]))]) for i in range(len(starflag))])
SA = theta[model_params+ins_params+1:]
assert len(SA) == totcornum, 'error in correlations'
AR = 0.0 # just to remember to add this
counter = -1
for i in range(nins):
for j in range(len(starflag[i])):
counter += 1
passer = -1
for k in range(ndat):
if starflag[i][j] == ins[k]: #
passer += 1
COR[i][j][k] = SA[counter] * staract[i][j][passer]
FMC = 0
for i in range(len(COR)):
for j in range(len(COR[i])):
FMC += COR[i][j]
else:
FMC = 0
MODEL = model(a1, time, kplanets) + offset + ACC + FMC
for i in range(ndat):
residuals[i] = rv[i] - MODEL[i]
for c in range(MOAV):
if i > c:
MA = macoef[c][i] * sp.exp(-sp.fabs(time[i-1] - time[i]) / timescale[c][i]) * residuals[i-1]
residuals[i] -= MA
inv_sigma2 = 1.0 / (err**2 + jitter**2)
lnl = sp.sum(residuals ** 2 * inv_sigma2 - sp.log(inv_sigma2)) + sp.log(2*sp.pi) * ndat
return -0.5 * lnl
def get_at_k1(Ek, Ek_grid, l, cb, conj, phi1k, hphi1k):
"""
Linearly interpolate the values of phi1k, hphi1k on Ek_grid at point Ek.
Parameters
----------
Ek : float
Energy value (not necessarily a grid point).
Ek_grid : array
Energy grid.
l : int
Lead label.
cb : int
Index corresponding to Phi[1](k) matrix element.
conj : bool
If conj=True the term in the integral equation is conjugated.
phi1k : array
Numpy array with dimensions (len(Ek_grid), nleads, ndm1, ndm0),
cotaining energy resolved first order density matrix elements Phi[1](k).
hphi1k : array
Numpy array with dimensions (len(Ek_grid_ext), nleads, ndm1, ndm0),
containing Hilbert transform of phi1k.
Returns
-------
array, array
Interpolated values of phi1k and hphi1k at Ek.
"""
if Ek<Ek_grid[0] or Ek>Ek_grid[-1]:
return 0, 0
#
#b_idx = np.searchsorted(Ek_grid, Ek)
b_idx = int((Ek-Ek_grid[0])/(Ek_grid[1]-Ek_grid[0]))+1
if b_idx == len(Ek_grid): b_idx -= 1
a_idx = b_idx - 1
b, a = Ek_grid[b_idx], Ek_grid[a_idx]
#print(a_idx, b_idx, a, Ek, b)
#
fb = phi1k[b_idx, l, cb]
fa = phi1k[a_idx, l, cb]
u = (fb-fa)/(b-a)*Ek + (b*fa-a*fb)/(b-a)
u = u.conjugate() if conj else u
#
fb = hphi1k[b_idx, l, cb]
fa = hphi1k[a_idx, l, cb]
hu = (fb-fa)/(b-a)*Ek + (b*fa-a*fb)/(b-a)
hu = hu.conjugate() if conj else hu
return pi*hu, pi*u