Python numpy 模块，pi() 实例源码
我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.pi()。
def deriveKernel(self, params, i):
self.checkParamsI(params, i)
ell = np.exp(params[0])
p = np.exp(params[1])
#compute d2
if (self.K_sq is None): d2 = sq_dist(self.X_scaled.T / ell) #precompute squared distances
else: d2 = self.K_sq / ell**2
#compute dp
dp = self.dp/p
K = np.exp(-d2 / 2.0)
if (i==0): return d2*K*np.cos(2*np.pi*dp)
elif (i==1): return 2*np.pi*dp*np.sin(2*np.pi*dp)*K
else: raise Exception('invalid parameter index:' + str(i))
def gelu(x):
return 0.5 * x * (1 + T.tanh(T.sqrt(2 / np.pi) * (x + 0.044715 * T.pow(x, 3))))
def get_local_wavenumbermesh(self, scaled=True, broadcast=False,
eliminate_highest_freq=False):
kx = fftfreq(self.N[0], 1./self.N[0])
ky = rfftfreq(self.N[1], 1./self.N[1])
if eliminate_highest_freq:
for i, k in enumerate((kx, ky)):
if self.N[i] % 2 == 0:
k[self.N[i]//2] = 0
Ks = np.meshgrid(kx, ky[self.rank*self.Np[1]//2:(self.rank*self.Np[1]//2+self.Npf)], indexing='ij', sparse=True)
if scaled is True:
Lp = 2*np.pi/self.L
Ks[0] *= Lp[0]
Ks[1] *= Lp[1]
K = Ks
if broadcast is True:
K = [np.broadcast_to(k, self.complex_shape()) for k in Ks]
return K
def _generate_data():
"""
?????
????u(k-1) ? y(k-1)?????y(k)
"""
# u = np.random.uniform(-1,1,200)
# y=[]
# former_y_value = 0
# for i in np.arange(0,200):
# y.append(former_y_value)
# next_y_value = (29.0 / 40) * np.sin(
# (16.0 * u[i] + 8 * former_y_value) / (3.0 + 4.0 * (u[i] ** 2) + 4 * (former_y_value ** 2))) \
# + (2.0 / 10) * u[i] + (2.0 / 10) * former_y_value
# former_y_value = next_y_value
# return u,y
u1 = np.random.uniform(-np.pi,np.pi,200)
u2 = np.random.uniform(-1,1,200)
y = np.zeros(200)
for i in range(200):
value = np.sin(u1[i]) + u2[i]
y[i] = value
return u1, u2, y
def ae(x):
if nonlinearity_name == 'relu':
f = tf.nn.relu
elif nonlinearity_name == 'elu':
f = tf.nn.elu
elif nonlinearity_name == 'gelu':
# def gelu(x):
# return tf.mul(x, tf.erfc(-x / tf.sqrt(2.)) / 2.)
# f = gelu
def gelu_fast(_x):
return 0.5 * _x * (1 + tf.tanh(tf.sqrt(2 / np.pi) * (_x + 0.044715 * tf.pow(_x, 3))))
f = gelu_fast
elif nonlinearity_name == 'silu':
def silu(_x):
return _x * tf.sigmoid(_x)
f = silu
# elif nonlinearity_name == 'soi':
# def soi_map(x):
# u = tf.random_uniform(tf.shape(x))
# mask = tf.to_float(tf.less(u, (1 + tf.erf(x / tf.sqrt(2.))) / 2.))
# return tf.cond(is_training, lambda: tf.mul(mask, x),
# lambda: tf.mul(x, tf.erfc(-x / tf.sqrt(2.)) / 2.))
# f = soi_map
else:
raise NameError("Need 'relu', 'elu', 'gelu', or 'silu' for nonlinearity_name")
h1 = f(tf.matmul(x, W['1']) + b['1'])
h2 = f(tf.matmul(h1, W['2']) + b['2'])
h3 = f(tf.matmul(h2, W['3']) + b['3'])
h4 = f(tf.matmul(h3, W['4']) + b['4'])
h5 = f(tf.matmul(h4, W['5']) + b['5'])
h6 = f(tf.matmul(h5, W['6']) + b['6'])
h7 = f(tf.matmul(h6, W['7']) + b['7'])
return tf.matmul(h7, W['8']) + b['8']
def score_samples(self, X):
"""Return the log-likelihood of each sample
See. "Pattern Recognition and Machine Learning"
by C. Bishop, 12.2.1 p. 574
or http://www.miketipping.com/papers/met-mppca.pdf
Parameters
----------
X: array, shape(n_samples, n_features)
The data.
Returns
-------
ll: array, shape (n_samples,)
Log-likelihood of each sample under the current model
"""
check_is_fitted(self, 'mean_')
X = check_array(X)
Xr = X - self.mean_
n_features = X.shape[1]
log_like = np.zeros(X.shape[0])
precision = self.get_precision()
log_like = -.5 * (Xr * (np.dot(Xr, precision))).sum(axis=1)
log_like -= .5 * (n_features * log(2. * np.pi)
- fast_logdet(precision))
return log_like
def getTrainTestKernel(self, params, Xtest):
self.checkParams(params)
ell = np.exp(params[0])
p = np.exp(params[1])
Xtest_scaled = Xtest/np.sqrt(Xtest.shape[1])
d2 = sq_dist(self.X_scaled.T/ell, Xtest_scaled.T/ell) #precompute squared distances
#compute dp
dp = np.zeros(d2.shape)
for d in xrange(self.X_scaled.shape[1]):
dp += (np.outer(self.X_scaled[:,d], np.ones((1, Xtest_scaled.shape[0]))) - np.outer(np.ones((self.X_scaled.shape[0], 1)), Xtest_scaled[:,d]))
dp /= p
K = np.exp(-d2 / 2.0)
return np.cos(2*np.pi*dp)*K
def reset(self,random_start_state=False, assign_state = False, n=None, k = None, \
perturb_params = False, p_LINK_LENGTH_1 = 0, p_LINK_LENGTH_2 = 0, \
p_LINK_MASS_1 = 0, p_LINK_MASS_2 = 0, **kw):
self.t = 0
self.state = np.random.uniform(low=-0.1,high=0.1,size=(4,))
self.LINK_LENGTH_1 = 1. # [m]
self.LINK_LENGTH_2 = 1. # [m]
self.LINK_MASS_1 = 1. #: [kg] mass of link 1
self.LINK_MASS_2 = 1.
if perturb_params:
self.LINK_LENGTH_1 += (self.LINK_LENGTH_1 * p_LINK_LENGTH_1) # [m]
self.LINK_LENGTH_2 += (self.LINK_LENGTH_2 * p_LINK_LENGTH_2) # [m]
self.LINK_MASS_1 += (self.LINK_MASS_1 * p_LINK_MASS_1) #: [kg] mass of link 1
self.LINK_MASS_2 += (self.LINK_MASS_2 * p_LINK_MASS_2) #: [kg] mass of link 2
# The idea here is that we can initialize our batch randomly so that we can get
# more variety in the state space that we attempt to fit a policy to.
if random_start_state:
self.state[:2] = np.random.uniform(-np.pi,np.pi,size=2)
if assign_state:
self.state[0] = wrap((2*k*np.pi)/(1.0*n),-np.pi,np.pi)
def calc_reward(self, action = None, state = None , **kw ):
'''Calculates the continuous reward based on the height of the foot (y position)
with a penalty applied if the hinge is moving (we want the acrobot to be upright
and stationary!), which is then normalized by the combined lengths of the links'''
t = self.target
if state is None:
s = self.state
else:
s = state
# Make sure that input state is clipped/wrapped to the given bounds (not guaranteed when coming from the BNN)
s[0] = wrap( s[0] , -np.pi , np.pi )
s[1] = wrap( s[1] , -np.pi , np.pi )
s[2] = bound( s[2] , -self.MAX_VEL_1 , self.MAX_VEL_1 )
s[3] = bound( s[3] , -self.MAX_VEL_1 , self.MAX_VEL_1 )
hinge, foot = self.get_cartesian_points(s)
reward = -0.05 * (foot[0] - self.LINK_LENGTH_1)**2
terminal = self.is_terminal(s)
return 10 if terminal else reward
def EStep(self):
P = np.zeros((self.M, self.N))
for i in range(0, self.M):
diff = self.X - np.tile(self.TY[i, :], (self.N, 1))
diff = np.multiply(diff, diff)
P[i, :] = P[i, :] + np.sum(diff, axis=1)
c = (2 * np.pi * self.sigma2) ** (self.D / 2)
c = c * self.w / (1 - self.w)
c = c * self.M / self.N
P = np.exp(-P / (2 * self.sigma2))
den = np.sum(P, axis=0)
den = np.tile(den, (self.M, 1))
den[den==0] = np.finfo(float).eps
self.P = np.divide(P, den)
self.Pt1 = np.sum(self.P, axis=0)
self.P1 = np.sum(self.P, axis=1)
self.Np = np.sum(self.P1)
def create_reference_image(size, x0=10., y0=-3., sigma_x=50., sigma_y=30., dtype='float64',
reverse_xaxis=False, correct_axes=True, sizey=None, **kwargs):
"""
Creates a reference image: a gaussian brightness with elliptical
"""
inc_cos = np.cos(0./180.*np.pi)
delta_x = 1.
x = (np.linspace(0., size - 1, size) - size / 2.) * delta_x
if sizey:
y = (np.linspace(0., sizey-1, sizey) - sizey/2.) * delta_x
else:
y = x.copy()
if reverse_xaxis:
xx, yy = np.meshgrid(-x, y/inc_cos)
elif correct_axes:
xx, yy = np.meshgrid(-x, -y/inc_cos)
else:
xx, yy = np.meshgrid(x, y/inc_cos)
image = np.exp(-(xx-x0)**2./sigma_x - (yy-y0)**2./sigma_y)
return image.astype(dtype)
def rotate_point_cloud(batch_data):
""" Randomly rotate the point clouds to augument the dataset
rotation is per shape based along up direction
Input:
BxNx3 array, original batch of point clouds
Return:
BxNx3 array, rotated batch of point clouds
"""
rotated_data = np.zeros(batch_data.shape, dtype=np.float32)
for k in range(batch_data.shape[0]):
rotation_angle = np.random.uniform() * 2 * np.pi
cosval = np.cos(rotation_angle)
sinval = np.sin(rotation_angle)
rotation_matrix = np.array([[cosval, 0, sinval],
[0, 1, 0],
[-sinval, 0, cosval]])
shape_pc = batch_data[k, ...]
rotated_data[k, ...] = np.dot(shape_pc.reshape((-1, 3)), rotation_matrix)
return rotated_data
def rotate_point_cloud_by_angle(batch_data, rotation_angle):
""" Rotate the point cloud along up direction with certain angle.
Input:
BxNx3 array, original batch of point clouds
Return:
BxNx3 array, rotated batch of point clouds
"""
rotated_data = np.zeros(batch_data.shape, dtype=np.float32)
for k in range(batch_data.shape[0]):
#rotation_angle = np.random.uniform() * 2 * np.pi
cosval = np.cos(rotation_angle)
sinval = np.sin(rotation_angle)
rotation_matrix = np.array([[cosval, 0, sinval],
[0, 1, 0],
[-sinval, 0, cosval]])
shape_pc = batch_data[k, ...]
rotated_data[k, ...] = np.dot(shape_pc.reshape((-1, 3)), rotation_matrix)
return rotated_data
def ac_solve(net):
"""
:param net:
:return:
"""
net.conductance_matrix()
net.dynamic_matrix()
net.rhs_matrix()
# frequency
f = float(net.analysis[-1])
# linear system definition
net.x = spsolve(net.G + 1j * 2 * np.pi * f* net.C, net.rhs)
def thinking(self):
"""Deliberate to avoid obstacles on the path."""
if self.motion.moveIsActive():
# Maneuver occurring. Let's finish it
# before taking any other measure.
pass
elif not self.sensors['proximity'][0].imminent_collision:
# Goes back to moving state.
self.behavior_ = self.BEHAVIORS.moving
elif all(s.imminent_collision for s in self.sensors['proximity']):
# There's nothing left to be done, only flag this is a dead-end.
self.behavior_ = self.BEHAVIORS.stuck
else:
peripheral_sensors = self.sensors['proximity'][1:]
for maneuver, sensor in zip(range(1, 4), peripheral_sensors):
if not sensor.imminent_collision:
# A sensor that indicates no obstacles were found.
# Move in that direction.
self.motion.post.moveTo(0, 0, np.pi / 2)
break
return self
def gaussian_nll(x, mus, sigmas):
"""
NLL for Multivariate Normal with diagonal covariance matrix
See:
wikipedia.org/wiki/Multivariate_normal_distribution#Likelihood_function
where \Sigma = diag(s_1^2,..., s_n^2).
x, mus, sigmas all should have the same shape.
sigmas (s_1,..., s_n) should be strictly positive.
Results in output shape of similar but without the last dimension.
"""
nll = lib.floatX(numpy.log(2. * numpy.pi))
nll += 2. * T.log(sigmas)
nll += ((x - mus) / sigmas) ** 2.
nll = nll.sum(axis=-1)
nll *= lib.floatX(0.5)
return nll
def tsukuba_load_poses(fn):
"""
Retrieve poses
X Y Z R P Y - > X -Y -Z R -P -Y
np.deg2rad(p[3]),-np.deg2rad(p[4]),-np.deg2rad(p[5]),
p[0]*.01,-p[1]*.01,-p[2]*.01, axes='sxyz') for p in P ]
"""
P = np.loadtxt(os.path.expanduser(fn), dtype=np.float64, delimiter=',')
return [ RigidTransform.from_rpyxyz(np.pi, 0, 0, 0, 0, 0) * \
RigidTransform.from_rpyxyz(
np.deg2rad(p[3]),np.deg2rad(p[4]),np.deg2rad(p[5]),
p[0]*.01,p[1]*.01,p[2]*.01, axes='sxyz') * \
RigidTransform.from_rpyxyz(np.pi, 0, 0, 0, 0, 0) for p in P ]
# return [ RigidTransform.from_rpyxyz(
# np.deg2rad(p[3]),-np.deg2rad(p[4]),-np.deg2rad(p[5]),
# p[0]*.01,-p[1]*.01,-p[2]*.01, axes='sxyz') for p in P ]
def __call__(self, z):
z1 = tf.reshape(tf.slice(z, [0, 0], [-1, 1]), [-1])
z2 = tf.reshape(tf.slice(z, [0, 1], [-1, 1]), [-1])
v1 = tf.sqrt((z1 - 5) * (z1 - 5) + z2 * z2) * 2
v2 = tf.sqrt((z1 + 5) * (z1 + 5) + z2 * z2) * 2
v3 = tf.sqrt((z1 - 2.5) * (z1 - 2.5) + (z2 - 2.5 * np.sqrt(3)) * (z2 - 2.5 * np.sqrt(3))) * 2
v4 = tf.sqrt((z1 + 2.5) * (z1 + 2.5) + (z2 + 2.5 * np.sqrt(3)) * (z2 + 2.5 * np.sqrt(3))) * 2
v5 = tf.sqrt((z1 - 2.5) * (z1 - 2.5) + (z2 + 2.5 * np.sqrt(3)) * (z2 + 2.5 * np.sqrt(3))) * 2
v6 = tf.sqrt((z1 + 2.5) * (z1 + 2.5) + (z2 - 2.5 * np.sqrt(3)) * (z2 - 2.5 * np.sqrt(3))) * 2
pdf1 = tf.exp(-0.5 * v1 * v1) / tf.sqrt(2 * np.pi * 0.25)
pdf2 = tf.exp(-0.5 * v2 * v2) / tf.sqrt(2 * np.pi * 0.25)
pdf3 = tf.exp(-0.5 * v3 * v3) / tf.sqrt(2 * np.pi * 0.25)
pdf4 = tf.exp(-0.5 * v4 * v4) / tf.sqrt(2 * np.pi * 0.25)
pdf5 = tf.exp(-0.5 * v5 * v5) / tf.sqrt(2 * np.pi * 0.25)
pdf6 = tf.exp(-0.5 * v6 * v6) / tf.sqrt(2 * np.pi * 0.25)
return -tf.log((pdf1 + pdf2 + pdf3 + pdf4 + pdf5 + pdf6) / 6)
def _evalfull(self, x):
fadd = self.fopt
curshape, dim = self.shape_(x)
# it is assumed x are row vectors
if self.lastshape != curshape:
self.initwithsize(curshape, dim)
# BOUNDARY HANDLING
# TRANSFORMATION IN SEARCH SPACE
x = x - self.arrxopt
x = monotoneTFosc(x)
idx = (x > 0)
x[idx] = x[idx] ** (1 + self.arrexpo[idx] * np.sqrt(x[idx]))
x = self.arrscales * x
# COMPUTATION core
ftrue = 10 * (self.dim - np.sum(np.cos(2 * np.pi * x), -1)) + np.sum(x ** 2, -1)
fval = self.noise(ftrue) # without noise
# FINALIZE
ftrue += fadd
fval += fadd
return fval, ftrue
def initwithsize(self, curshape, dim):
# DIM-dependent initialization
if self.dim != dim:
if self.zerox:
self.xopt = zeros(dim)
else:
self.xopt = compute_xopt(self.rseed, dim)
self.rotation = compute_rotation(self.rseed + 1e6, dim)
self.scales = (1. / self.condition ** .5) ** linspace(0, 1, dim) # CAVE?
self.linearTF = dot(compute_rotation(self.rseed, dim), diag(self.scales))
# decouple scaling from function definition
self.linearTF = dot(self.linearTF, self.rotation)
K = np.arange(0, 12)
self.aK = np.reshape(0.5 ** K, (1, 12))
self.bK = np.reshape(3. ** K, (1, 12))
self.f0 = np.sum(self.aK * np.cos(2 * np.pi * self.bK * 0.5)) # optimal value
# DIM- and POPSI-dependent initialisations of DIM*POPSI matrices
if self.lastshape != curshape:
self.dim = dim
self.lastshape = curshape
self.arrxopt = resize(self.xopt, curshape)
def pan(self, dx, dy, dz, relative=False):
"""
Moves the center (look-at) position while holding the camera in place.
If relative=True, then the coordinates are interpreted such that x
if in the global xy plane and points to the right side of the view, y is
in the global xy plane and orthogonal to x, and z points in the global z
direction. Distances are scaled roughly such that a value of 1.0 moves
by one pixel on screen.
"""
if not relative:
self.camera_center += QtGui.QVector3D(dx, dy, dz)
else:
cPos = self.cameraPosition()
cVec = self.camera_center - cPos
dist = cVec.length() ## distance from camera to center
xDist = dist * 2. * np.tan(0.5 * self.camera_fov * np.pi / 180.) ## approx. width of view at distance of center point
xScale = xDist / self.width()
zVec = QtGui.QVector3D(0,0,1)
xVec = QtGui.QVector3D.crossProduct(zVec, cVec).normalized()
yVec = QtGui.QVector3D.crossProduct(xVec, zVec).normalized()
self.camera_center = self.camera_center + xVec * xScale * dx + yVec * xScale * dy + zVec * xScale * dz
self.update()
def test_pitch_estimation(self):
"""
test pitch estimation algo with contrived small example
if pitch is within 5 Hz, then say its good (for this small example,
since the algorithm wasn't made for this type of synthesized signal)
"""
cfg = ExperimentConfig(pitch_strength_thresh=-np.inf)
# the next 3 variables are in Hz
tolerance = 5
fs = 48000
f = 150
# create a sine wave of f Hz freq sampled at fs Hz
x = np.sin(2*np.pi * f/fs * np.arange(2**10))
# estimate the pitch, it should be close to f
p, t, s = pest.pitch_estimation(x, fs, cfg)
self.assertTrue(np.all(np.abs(p - f) < tolerance))
def setFromQTransform(self, tr):
p1 = Point(tr.map(0., 0.))
p2 = Point(tr.map(1., 0.))
p3 = Point(tr.map(0., 1.))
dp2 = Point(p2-p1)
dp3 = Point(p3-p1)
## detect flipped axes
if dp2.angle(dp3) > 0:
#da = 180
da = 0
sy = -1.0
else:
da = 0
sy = 1.0
self._state = {
'pos': Point(p1),
'scale': Point(dp2.length(), dp3.length() * sy),
'angle': (np.arctan2(dp2[1], dp2[0]) * 180. / np.pi) + da
}
self.update()
def projectionMatrix(self, region=None):
# Xw = (Xnd + 1) * width/2 + X
if region is None:
region = (0, 0, self.width(), self.height())
x0, y0, w, h = self.getViewport()
dist = self.opts['distance']
fov = self.opts['fov']
nearClip = dist * 0.001
farClip = dist * 1000.
r = nearClip * np.tan(fov * 0.5 * np.pi / 180.)
t = r * h / w
# convert screen coordinates (region) to normalized device coordinates
# Xnd = (Xw - X0) * 2/width - 1
## Note that X0 and width in these equations must be the values used in viewport
left = r * ((region[0]-x0) * (2.0/w) - 1)
right = r * ((region[0]+region[2]-x0) * (2.0/w) - 1)
bottom = t * ((region[1]-y0) * (2.0/h) - 1)
top = t * ((region[1]+region[3]-y0) * (2.0/h) - 1)
tr = QtGui.QMatrix4x4()
tr.frustum(left, right, bottom, top, nearClip, farClip)
return tr
def pan(self, dx, dy, dz, relative=False):
"""
Moves the center (look-at) position while holding the camera in place.
If relative=True, then the coordinates are interpreted such that x
if in the global xy plane and points to the right side of the view, y is
in the global xy plane and orthogonal to x, and z points in the global z
direction. Distances are scaled roughly such that a value of 1.0 moves
by one pixel on screen.
"""
if not relative:
self.opts['center'] += QtGui.QVector3D(dx, dy, dz)
else:
cPos = self.cameraPosition()
cVec = self.opts['center'] - cPos
dist = cVec.length() ## distance from camera to center
xDist = dist * 2. * np.tan(0.5 * self.opts['fov'] * np.pi / 180.) ## approx. width of view at distance of center point
xScale = xDist / self.width()
zVec = QtGui.QVector3D(0,0,1)
xVec = QtGui.QVector3D.crossProduct(zVec, cVec).normalized()
yVec = QtGui.QVector3D.crossProduct(xVec, zVec).normalized()
self.opts['center'] = self.opts['center'] + xVec * xScale * dx + yVec * xScale * dy + zVec * xScale * dz
self.update()
def makeArrowPath(headLen=20, tipAngle=20, tailLen=20, tailWidth=3, baseAngle=0):
"""
Construct a path outlining an arrow with the given dimensions.
The arrow points in the -x direction with tip positioned at 0,0.
If *tipAngle* is supplied (in degrees), it overrides *headWidth*.
If *tailLen* is None, no tail will be drawn.
"""
headWidth = headLen * np.tan(tipAngle * 0.5 * np.pi/180.)
path = QtGui.QPainterPath()
path.moveTo(0,0)
path.lineTo(headLen, -headWidth)
if tailLen is None:
innerY = headLen - headWidth * np.tan(baseAngle*np.pi/180.)
path.lineTo(innerY, 0)
else:
tailWidth *= 0.5
innerY = headLen - (headWidth-tailWidth) * np.tan(baseAngle*np.pi/180.)
path.lineTo(innerY, -tailWidth)
path.lineTo(headLen + tailLen, -tailWidth)
path.lineTo(headLen + tailLen, tailWidth)
path.lineTo(innerY, tailWidth)
path.lineTo(headLen, headWidth)
path.lineTo(0,0)
return path
def setFromQTransform(self, tr):
p1 = Point(tr.map(0., 0.))
p2 = Point(tr.map(1., 0.))
p3 = Point(tr.map(0., 1.))
dp2 = Point(p2-p1)
dp3 = Point(p3-p1)
## detect flipped axes
if dp2.angle(dp3) > 0:
#da = 180
da = 0
sy = -1.0
else:
da = 0
sy = 1.0
self._state = {
'pos': Point(p1),
'scale': Point(dp2.length(), dp3.length() * sy),
'angle': (np.arctan2(dp2[1], dp2[0]) * 180. / np.pi) + da
}
self.update()
def projectionMatrix(self, region=None):
# Xw = (Xnd + 1) * width/2 + X
if region is None:
region = (0, 0, self.width(), self.height())
x0, y0, w, h = self.getViewport()
dist = self.opts['distance']
fov = self.opts['fov']
nearClip = dist * 0.001
farClip = dist * 1000.
r = nearClip * np.tan(fov * 0.5 * np.pi / 180.)
t = r * h / w
# convert screen coordinates (region) to normalized device coordinates
# Xnd = (Xw - X0) * 2/width - 1
## Note that X0 and width in these equations must be the values used in viewport
left = r * ((region[0]-x0) * (2.0/w) - 1)
right = r * ((region[0]+region[2]-x0) * (2.0/w) - 1)
bottom = t * ((region[1]-y0) * (2.0/h) - 1)
top = t * ((region[1]+region[3]-y0) * (2.0/h) - 1)
tr = QtGui.QMatrix4x4()
tr.frustum(left, right, bottom, top, nearClip, farClip)
return tr
def pan(self, dx, dy, dz, relative=False):
"""
Moves the center (look-at) position while holding the camera in place.
If relative=True, then the coordinates are interpreted such that x
if in the global xy plane and points to the right side of the view, y is
in the global xy plane and orthogonal to x, and z points in the global z
direction. Distances are scaled roughly such that a value of 1.0 moves
by one pixel on screen.
"""
if not relative:
self.opts['center'] += QtGui.QVector3D(dx, dy, dz)
else:
cPos = self.cameraPosition()
cVec = self.opts['center'] - cPos
dist = cVec.length() ## distance from camera to center
xDist = dist * 2. * np.tan(0.5 * self.opts['fov'] * np.pi / 180.) ## approx. width of view at distance of center point
xScale = xDist / self.width()
zVec = QtGui.QVector3D(0,0,1)
xVec = QtGui.QVector3D.crossProduct(zVec, cVec).normalized()
yVec = QtGui.QVector3D.crossProduct(xVec, zVec).normalized()
self.opts['center'] = self.opts['center'] + xVec * xScale * dx + yVec * xScale * dy + zVec * xScale * dz
self.update()
def make_wafer(self,wafer_r,frame,label,labelloc,labelwidth):
"""
Generate wafer with primary flat on the left. From https://coresix.com/products/wafers/ I estimated that the angle defining the wafer flat to arctan(flat/2 / radius)
"""
angled = 18
angle = angled*np.pi/180
circ = cad.shapes.Circle((0,0), wafer_r, width=self.boxwidth, initial_angle=180+angled, final_angle=360+180-angled, layer=self.layer_box)
flat = cad.core.Path([(-wafer_r*np.cos(angle),wafer_r*np.sin(angle)),(-wafer_r*np.cos(angle),-wafer_r*np.sin(angle))], width=self.boxwidth, layer=self.layer_box)
date = time.strftime("%d/%m/%Y")
if labelloc==(0,0):
labelloc=(-2e3,wafer_r-1e3)
# The label is added 100 um on top of the main cell
label_grid_chip = cad.shapes.LineLabel( self.name + " " +\
date,500,position=labelloc,
line_width=labelwidth,
layer=self.layer_label)
if frame==True:
self.add(circ)
self.add(flat)
if label==True:
self.add(label_grid_chip)
def evaluation(self, X_test, y_test):
# normalization
X_test = self.normalization(X_test)
# average over the output
pred_y_test = np.zeros([self.M, len(y_test)])
prob = np.zeros([self.M, len(y_test)])
'''
Since we have M particles, we use a Bayesian view to calculate rmse and log-likelihood
'''
for i in range(self.M):
w1, b1, w2, b2, loggamma, loglambda = self.unpack_weights(self.theta[i, :])
pred_y_test[i, :] = self.nn_predict(X_test, w1, b1, w2, b2) * self.std_y_train + self.mean_y_train
prob[i, :] = np.sqrt(np.exp(loggamma)) /np.sqrt(2*np.pi) * np.exp( -1 * (np.power(pred_y_test[i, :] - y_test, 2) / 2) * np.exp(loggamma) )
pred = np.mean(pred_y_test, axis=0)
# evaluation
svgd_rmse = np.sqrt(np.mean((pred - y_test)**2))
svgd_ll = np.mean(np.log(np.mean(prob, axis = 0)))
return (svgd_rmse, svgd_ll)
def nufft_scale1(N, K, alpha, beta, Nmid):
'''
calculate image space scaling factor
'''
# import types
# if alpha is types.ComplexType:
alpha = numpy.real(alpha)
# print('complex alpha may not work, but I just let it as')
L = len(alpha) - 1
if L > 0:
sn = numpy.zeros((N, 1))
n = numpy.arange(0, N).reshape((N, 1), order='F')
i_gam_n_n0 = 1j * (2 * numpy.pi / K) * (n - Nmid) * beta
for l1 in range(-L, L + 1):
alf = alpha[abs(l1)]
if l1 < 0:
alf = numpy.conj(alf)
sn = sn + alf * numpy.exp(i_gam_n_n0 * l1)
else:
sn = numpy.dot(alpha, numpy.ones((N, 1), dtype=numpy.float32))
return sn
def nufft_r(om, N, J, K, alpha, beta):
'''
equation (30) of Fessler's paper
'''
M = numpy.size(om) # 1D size
gam = 2.0 * numpy.pi / (K * 1.0)
nufft_offset0 = nufft_offset(om, J, K) # om/gam - nufft_offset , [M,1]
dk = 1.0 * om / gam - nufft_offset0 # om/gam - nufft_offset , [M,1]
arg = outer_sum(-numpy.arange(1, J + 1) * 1.0, dk)
L = numpy.size(alpha) - 1
# print('alpha',alpha)
rr = numpy.zeros((J, M), dtype=numpy.float32)
rr = iterate_l1(L, alpha, arg, beta, K, N, rr)
return (rr, arg)
def kaiser_bessel_ft(u, J, alpha, kb_m, d):
'''
Interpolation weight for given J/alpha/kb-m
'''
u = u * (1.0 + 0.0j)
import scipy.special
z = numpy.sqrt((2 * numpy.pi * (J / 2) * u) ** 2.0 - alpha ** 2.0)
nu = d / 2 + kb_m
y = ((2 * numpy.pi) ** (d / 2)) * ((J / 2) ** d) * (alpha ** kb_m) / \
scipy.special.iv(kb_m, alpha) * scipy.special.jv(nu, z) / (z ** nu)
y = numpy.real(y)
return y
def nufft_scale1(N, K, alpha, beta, Nmid):
'''
Calculate image space scaling factor
'''
# import types
# if alpha is types.ComplexType:
alpha = numpy.real(alpha)
# print('complex alpha may not work, but I just let it as')
L = len(alpha) - 1
if L > 0:
sn = numpy.zeros((N, 1))
n = numpy.arange(0, N).reshape((N, 1), order='F')
i_gam_n_n0 = 1j * (2 * numpy.pi / K) * (n - Nmid) * beta
for l1 in range(-L, L + 1):
alf = alpha[abs(l1)]
if l1 < 0:
alf = numpy.conj(alf)
sn = sn + alf * numpy.exp(i_gam_n_n0 * l1)
else:
sn = numpy.dot(alpha, numpy.ones((N, 1)))
return sn
def planetary_radius(mass, radius):
"""Calculate planetary radius if not given assuming a density dependent on
mass"""
if not isinstance(mass, (int, float)):
if isinstance(radius, (int, float)):
return radius
else:
return '...'
if mass < 0:
raise ValueError('Only positive planetary masses allowed.')
Mj = c.M_jup
Rj = c.R_jup
if radius == '...' and isinstance(mass, (int, float)):
if mass < 0.01: # Earth density
rho = 5.51
elif 0.01 <= mass <= 0.5:
rho = 1.64 # Neptune density
else:
rho = Mj/(4./3*np.pi*Rj**3) # Jupiter density
R = ((mass*Mj)/(4./3*np.pi*rho))**(1./3) # Neptune density
R /= Rj
else:
return radius
return R.value
def test_kbd():
M = 100
w = mdct.windows.kaiser_derived(M, beta=4.)
assert numpy.allclose(w[:M//2] ** 2 + w[-M//2:] ** 2, 1.)
with pytest.raises(ValueError):
mdct.windows.kaiser_derived(M + 1, beta=4.)
assert numpy.allclose(
mdct.windows.kaiser_derived(2, beta=numpy.pi/2)[:1],
[numpy.sqrt(2)/2])
assert numpy.allclose(
mdct.windows.kaiser_derived(4, beta=numpy.pi/2)[:2],
[0.518562710536, 0.855039598640])
assert numpy.allclose(
mdct.windows.kaiser_derived(6, beta=numpy.pi/2)[:3],
[0.436168993154, 0.707106781187, 0.899864772847])
def gaussian_kernel(kernel_shape, sigma=None):
"""
Get 2D Gaussian kernel
:param kernel_shape: kernel size
:param sigma: sigma of Gaussian distribution
:return: 2D Gaussian kernel
"""
kern = numpy.zeros((kernel_shape, kernel_shape), dtype='float32')
# get sigma from kernel size
if sigma is None:
sigma = 0.3*((kernel_shape-1.)*0.5 - 1.) + 0.8
def gauss(x, y, s):
Z = 2. * numpy.pi * s ** 2.
return 1. / Z * numpy.exp(-(x ** 2. + y ** 2.) / (2. * s ** 2.))
mid = numpy.floor(kernel_shape / 2.)
for i in xrange(0, kernel_shape):
for j in xrange(0, kernel_shape):
kern[i, j] = gauss(i - mid, j - mid, sigma)
return kern / kern.sum()
def is_grid(self, grid, image):
"""
Checks the "gridness" by analyzing the results of a hough transform.
:param grid: binary image
:return: wheter the object in the image might be a grid or not
"""
# - Distance resolution = 1 pixel
# - Angle resolution = 1° degree for high line density
# - Threshold = 144 hough intersections
# 8px digit + 3*2px white + 2*1px border = 16px per cell
# => 144x144 grid
# 144 - minimum number of points on the same line
# (but due to imperfections in the binarized image it's highly
# improbable to detect a 144x144 grid)
lines = cv2.HoughLines(grid, 1, np.pi / 180, 144)
if lines is not None and np.size(lines) >= 20:
lines = lines.reshape((lines.size / 2), 2)
# theta in [0, pi] (theta > pi => rho < 0)
# normalise theta in [-pi, pi] and negatives rho
lines[lines[:, 0] < 0, 1] -= np.pi
lines[lines[:, 0] < 0, 0] *= -1
criteria = (cv2.TERM_CRITERIA_EPS, 0, 0.01)
# split lines into 2 groups to check whether they're perpendicular
if cv2.__version__[0] == '2':
density, clmap, centers = cv2.kmeans(
lines[:, 1], 2, criteria, 5, cv2.KMEANS_RANDOM_CENTERS)
else:
density, clmap, centers = cv2.kmeans(
lines[:, 1], 2, None, criteria,
5, cv2.KMEANS_RANDOM_CENTERS)
if self.debug:
self.save_hough(lines, clmap)
# Overall variance from respective centers
var = density / np.size(clmap)
sin = abs(np.sin(centers[0] - centers[1]))
# It is probably a grid only if:
# - centroids difference is almost a 90° angle (+-15° limit)
# - variance is less than 5° (keeping in mind surface distortions)
return sin > 0.99 and var <= (5*np.pi / 180) ** 2
else:
return False
def save_hough(self, lines, clmap):
"""
:param lines: (rho, theta) pairs
:param clmap: clusters assigned to lines
:return: None
"""
height, width = self.image.shape
ratio = 600. * (self.step+1) / min(height, width)
temp = cv2.resize(self.image, None, fx=ratio, fy=ratio,
interpolation=cv2.INTER_CUBIC)
temp = cv2.cvtColor(temp, cv2.COLOR_GRAY2BGR)
colors = [(0, 127, 255), (255, 0, 127)]
for i in range(0, np.size(lines) / 2):
rho = lines[i, 0]
theta = lines[i, 1]
color = colors[clmap[i, 0]]
if theta < np.pi / 4 or theta > 3 * np.pi / 4:
pt1 = (rho / np.cos(theta), 0)
pt2 = (rho - height * np.sin(theta) / np.cos(theta), height)
else:
pt1 = (0, rho / np.sin(theta))
pt2 = (width, (rho - width * np.cos(theta)) / np.sin(theta))
pt1 = (int(pt1[0]), int(pt1[1]))
pt2 = (int(pt2[0]), int(pt2[1]))
cv2.line(temp, pt1, pt2, color, 5)
self.save2image(temp)
def is_grid(self, grid, image):
"""
Checks the "gridness" by analyzing the results of a hough transform.
:param grid: binary image
:return: wheter the object in the image might be a grid or not
"""
# - Distance resolution = 1 pixel
# - Angle resolution = 1° degree for high line density
# - Threshold = 144 hough intersections
# 8px digit + 3*2px white + 2*1px border = 16px per cell
# => 144x144 grid
# 144 - minimum number of points on the same line
# (but due to imperfections in the binarized image it's highly
# improbable to detect a 144x144 grid)
lines = cv2.HoughLines(grid, 1, np.pi / 180, 144)
if lines is not None and np.size(lines) >= 20:
lines = lines.reshape((lines.size/2), 2)
# theta in [0, pi] (theta > pi => rho < 0)
# normalise theta in [-pi, pi] and negatives rho
lines[lines[:, 0] < 0, 1] -= np.pi
lines[lines[:, 0] < 0, 0] *= -1
criteria = (cv2.TERM_CRITERIA_EPS, 0, 0.01)
# split lines into 2 groups to check whether they're perpendicular
if cv2.__version__[0] == '2':
density, clmap, centers = cv2.kmeans(
lines[:, 1], 2, criteria,
5, cv2.KMEANS_RANDOM_CENTERS)
else:
density, clmap, centers = cv2.kmeans(
lines[:, 1], 2, None, criteria,
5, cv2.KMEANS_RANDOM_CENTERS)
# Overall variance from respective centers
var = density / np.size(clmap)
sin = abs(np.sin(centers[0] - centers[1]))
# It is probably a grid only if:
# - centroids difference is almost a 90° angle (+-15° limit)
# - variance is less than 5° (keeping in mind surface distortions)
return sin > 0.99 and var <= (5*np.pi / 180) ** 2
else:
return False
def build_2D_cov_matrix(sigmax,sigmay,angle,verbose=True):
"""
Build a covariance matrix for a 2D multivariate Gaussian
--- INPUT ---
sigmax Standard deviation of the x-compoent of the multivariate Gaussian
sigmay Standard deviation of the y-compoent of the multivariate Gaussian
angle Angle to rotate matrix by in degrees (clockwise) to populate covariance cross terms
verbose Toggle verbosity
--- EXAMPLE OF USE ---
import tdose_utilities as tu
covmatrix = tu.build_2D_cov_matrix(3,1,35)
"""
if verbose: print ' - Build 2D covariance matrix with varinaces (x,y)=('+str(sigmax)+','+str(sigmay)+\
') and then rotated '+str(angle)+' degrees'
cov_orig = np.zeros([2,2])
cov_orig[0,0] = sigmay**2.0
cov_orig[1,1] = sigmax**2.0
angle_rad = (180.0-angle) * np.pi/180.0 # The (90-angle) makes sure the same convention as DS9 is used
c, s = np.cos(angle_rad), np.sin(angle_rad)
rotmatrix = np.matrix([[c, -s], [s, c]])
cov_rot = np.dot(np.dot(rotmatrix,cov_orig),np.transpose(rotmatrix)) # performing rot * cov * rot^T
return cov_rot
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def normalize_2D_cov_matrix(covmatrix,verbose=True):
"""
Calculate the normalization foctor for a multivariate gaussian from it's covariance matrix
However, not that gaussian returned by tu.gen_2Dgauss() is normalized for scale=1
--- INPUT ---
covmatrix covariance matrix to normaliz
verbose Toggle verbosity
"""
detcov = np.linalg.det(covmatrix)
normfac = 1.0 / (2.0 * np.pi * np.sqrt(detcov) )
return normfac
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def f(r,theta):
out = np.sin(theta)*np.cos(K*2*np.pi*(1./r))/r
out[-1] = 0
return out
def dfdr(r,theta):
out = (2*K*np.pi*np.sin(2*np.pi*K/r)
-r*np.cos(2*np.pi*K/r))*np.sin(theta)/(r**3)
out[-1] = 0
return out
def dfdrdtheta(r,theta):
out = (2*K*np.pi*np.sin(2*np.pi*K/r)
-r*np.cos(2*np.pi*K/r))*np.cos(theta)/(r**3)
out[-1] = 0
return out
def __init__(self,
order_X,r_h,
order_theta,
theta_min = 0,
theta_max = np.pi,
L=1):
"""Constructor.
Parameters
----------
order_X -- polynomial order in X direction
r_h -- physical minimum radius (uncompactified coordinates)
order_theta -- polynomial order in theta direction
theta_min -- minimum longitudinal value. Should be no less than 0.
theta_max -- maximum longitudinal value. Should be no greater than pi.
L -- Characteristic length scale of problem.
Needed for compactification
"""
self.order_X = order_X
self.order_theta = order_theta
self.r_h = r_h
self.theta_min = theta_min
self.theta_max = theta_max
self.L = L
super(PyballdDiscretization,self).__init__(order_X,
self.X_min,self.X_max,
order_theta,
theta_min,theta_max)
self.r = self.get_r_from_X(self.x)
self.R = self.get_r_from_X(self.X)
self.dRdX = self.get_drdX(self.X)
self.drdX = self.get_drdX(self.x)
self.dXdR = self.get_dXdr(self.X)
self.dXdr = self.get_dXdr(self.x)
self.d2XdR2 = self.get_d2Xdr2(self.X)
self.d2Xdr2 = self.get_d2Xdr2(self.x)
self.d2RdX2 = self.get_d2rdX2(self.X)
self.d2rdX2 = self.get_d2rdX2(self.x)
self.theta = self.y
self.THETA = self.Y
def get_integration_weights(order,nodes=None):
"""
Returns the integration weights for Gauss-Lobatto quadrature
as a function of the order of the polynomial we want to
represent.
See: https://en.wikipedia.org/wiki/Gaussian_quadrature
See: arXive:gr-qc/0609020v1
"""
if np.all(nodes == False):
nodes=get_quadrature_points(order)
if poly == polynomial.chebyshev.Chebyshev:
weights = np.empty((order+1))
weights[1:-1] = np.pi/order
weights[0] = np.pi/(2*order)
weights[-1] = weights[0]
return weights
elif poly == polynomial.legendre.Legendre:
interior_weights = 2/((order+1)*order*poly.basis(order)(nodes[1:-1])**2)
boundary_weights = np.array([1-0.5*np.sum(interior_weights)])
weights = np.concatenate((boundary_weights,
interior_weights,
boundary_weights))
return weights
else:
raise ValueError("Not a known polynomial type.")
return False
def gelu_fast(_x):
return 0.5 * _x * (1 + tf.tanh(tf.sqrt(2 / np.pi) * (_x + 0.044715 * tf.pow(_x, 3))))