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))
python类cos()的实例源码
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 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 _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 genplot(x, y, fit, xdata=None, ydata=None, maxpts=10000):
bin_range = (0, 360)
a = (np.arange(*bin_range))
f_a = nuth_func(a, fit[0], fit[1], fit[2])
nuth_func_str = r'$y=%0.2f*cos(%0.2f-x)+%0.2f$' % tuple(fit)
if xdata.size > maxpts:
import random
idx = random.sample(list(range(xdata.size)), 10000)
else:
idx = np.arange(xdata.size)
f, ax = plt.subplots()
ax.set_xlabel('Aspect (deg)')
ax.set_ylabel('dh/tan(slope) (m)')
ax.plot(xdata[idx], ydata[idx], 'k.', label='Orig pixels')
ax.plot(x, y, 'ro', label='Bin median')
ax.axhline(color='k')
ax.plot(a, f_a, 'b', label=nuth_func_str)
ax.set_xlim(*bin_range)
pad = 0.2 * np.max([np.abs(y.min()), np.abs(y.max())])
ax.set_ylim(y.min() - pad, y.max() + pad)
ax.legend(prop={'size':8})
return f
#Function copied from from openPIV pyprocess
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 ani_update(arg, ii=[0]):
i = ii[0] # don't ask...
if np.isclose(t_arr[i], np.around(t_arr[i], 1)):
fig2.suptitle('Evolution (Time: {})'.format(t_arr[i]), fontsize=24)
graphic_floor[0].set_data([-floor_lim*np.cos(incline_history[i]) + radius*np.sin(incline_history[i]), floor_lim*np.cos(incline_history[i]) + radius*np.sin(incline_history[i])], [-floor_lim*np.sin(incline_history[i])-radius*np.cos(incline_history[i]), floor_lim*np.sin(incline_history[i])-radius*np.cos(incline_history[i])])
graphic_wheel.center = (x_history[i], y_history[i])
graphic_ind[0].set_data([x_history[i], x_history[i] + radius*np.sin(w_history[i])],
[y_history[i], y_history[i] + radius*np.cos(w_history[i])])
graphic_pend[0].set_data([x_history[i], x_history[i] - cw_to_cm[1]*np.sin(q_history[i, 2])],
[y_history[i], y_history[i] + cw_to_cm[1]*np.cos(q_history[i, 2])])
graphic_dist[0].set_data([x_history[i] - cw_to_cm[1]*np.sin(q_history[i, 2]), x_history[i] - cw_to_cm[1]*np.sin(q_history[i, 2]) - pscale*p_history[i]*np.cos(q_history[i, 2])],
[y_history[i] + cw_to_cm[1]*np.cos(q_history[i, 2]), y_history[i] + cw_to_cm[1]*np.cos(q_history[i, 2]) - pscale*p_history[i]*np.sin(q_history[i, 2])])
ii[0] += int(1 / (timestep * framerate))
if ii[0] >= len(t_arr):
print("Resetting animation!")
ii[0] = 0
return [graphic_floor, graphic_wheel, graphic_ind, graphic_pend, graphic_dist]
# Run animation
def getZ2(el,recpos):
"""
sqrt( [a+h+s]^2 - [a*cos(el)]^2 ) - sqrt( [a+h]^2 - [a*cos(el)]^2 )
z(el) = ---------------------------------------------------------------------
s
a is height of observing station from earth center in km,
h = 300km is height of ionosphere slab
s = 200km is slab thickness
"""
a = np.linalg.norm(recpos)/1000
h=300
s = 200
Z = (np.sqrt((a+h+s)**2-(a*np.cos(np.radians(el[el>30])))**2)
-np.sqrt((a+h)**2-(a*np.cos(np.radians(el[el>30])))**2))/s
return Z
def great_circle_dist(p1, p2):
"""Return the distance (in km) between two points in
geographical coordinates.
"""
lon0, lat0 = p1
lon1, lat1 = p2
EARTH_R = 6372.8
lat0 = np.radians(float(lat0))
lon0 = np.radians(float(lon0))
lat1 = np.radians(float(lat1))
lon1 = np.radians(float(lon1))
dlon = lon0 - lon1
y = np.sqrt(
(np.cos(lat1) * np.sin(dlon)) ** 2
+ (np.cos(lat0) * np.sin(lat1)
- np.sin(lat0) * np.cos(lat1) * np.cos(dlon)) ** 2)
x = np.sin(lat0) * np.sin(lat1) + \
np.cos(lat0) * np.cos(lat1) * np.cos(dlon)
c = np.arctan2(y, x)
return EARTH_R * c
def x_axis_rotation(theta):
"""Generates a 3x3 rotation matrix for a rotation of angle
theta about the x axis.
Parameters
----------
theta : float
amount to rotate, in radians
Returns
-------
:obj:`numpy.ndarray` of float
A random 3x3 rotation matrix.
"""
R = np.array([[1, 0, 0,],
[0, np.cos(theta), -np.sin(theta)],
[0, np.sin(theta), np.cos(theta)]])
return R
def y_axis_rotation(theta):
"""Generates a 3x3 rotation matrix for a rotation of angle
theta about the y axis.
Parameters
----------
theta : float
amount to rotate, in radians
Returns
-------
:obj:`numpy.ndarray` of float
A random 3x3 rotation matrix.
"""
R = np.array([[np.cos(theta), 0, np.sin(theta)],
[0, 1, 0],
[-np.sin(theta), 0, np.cos(theta)]])
return R
def z_axis_rotation(theta):
"""Generates a 3x3 rotation matrix for a rotation of angle
theta about the z axis.
Parameters
----------
theta : float
amount to rotate, in radians
Returns
-------
:obj:`numpy.ndarray` of float
A random 3x3 rotation matrix.
"""
R = np.array([[np.cos(theta), -np.sin(theta), 0],
[np.sin(theta), np.cos(theta), 0],
[0, 0, 1]])
return R
def sph2cart(r, az, elev):
""" Convert spherical to cartesian coordinates.
Attributes
----------
r : float
radius
az : float
aziumth (angle about z axis)
elev : float
elevation from xy plane
Returns
-------
float
x-coordinate
float
y-coordinate
float
z-coordinate
"""
x = r * np.cos(az) * np.sin(elev)
y = r * np.sin(az) * np.sin(elev)
z = r * np.cos(elev)
return x, y, z
def mdct(x, odd=True):
""" Calculate modified discrete cosine transform of input signal in an
inefficient pure-Python method.
Use only for testing.
Parameters
----------
X : array_like
The input signal
odd : boolean, optional
Switch to oddly stacked transform. Defaults to :code:`True`.
Returns
-------
out : array_like
The output signal
"""
return trans(x, func=numpy.cos, odd=odd) * numpy.sqrt(2)
def imdct(X, odd=True):
""" Calculate inverse modified discrete cosine transform of input
signal in an inefficient pure-Python method.
Use only for testing.
Parameters
----------
X : array_like
The input signal
odd : boolean, optional
Switch to oddly stacked transform. Defaults to :code:`True`.
Returns
-------
out : array_like
The output signal
"""
return itrans(X, func=numpy.cos, odd=odd) * numpy.sqrt(2)
def cmdct(x, odd=True):
""" Calculate complex modified discrete cosine transform of input
inefficient pure-Python method.
Use only for testing.
Parameters
----------
X : array_like
The input signal
odd : boolean, optional
Switch to oddly stacked transform. Defaults to :code:`True`.
Returns
-------
out : array_like
The output signal
"""
return trans(x, func=lambda x: numpy.cos(x) - 1j * numpy.sin(x), odd=odd)
def icmdct(X, odd=True):
""" Calculate inverse complex modified discrete cosine transform of input
signal in an inefficient pure-Python method.
Use only for testing.
Parameters
----------
X : array_like
The input signal
odd : boolean, optional
Switch to oddly stacked transform. Defaults to :code:`True`.
Returns
-------
out : array_like
The output signal
"""
return itrans(X, func=lambda x: numpy.cos(x) + 1j * numpy.sin(x), odd=odd)
def test_with_fake_log_prob(self):
np.random.seed(42)
def grad_log_prob(x):
return -(x/2.0 + np.sin(x))*(1.0/2.0 + np.cos(x))
def fake_log_prob(x):
return -(x/5.0 + np.sin(x) )**2.0/2.0
generator = mh_generator(log_density=fake_log_prob,x_start=1.0)
tester = GaussianSteinTest(grad_log_prob,41)
selector = SampleSelector(generator, sample_size=1000,thinning=20,tester=tester, max_iterations=5)
data,converged = selector.points_from_stationary()
assert converged is False
def test_with_ugly(self):
np.random.seed(42)
def grad_log_prob(x):
return -(x/5.0 + np.sin(x))*(1.0/5.0 + np.cos(x))
def log_prob(x):
return -(x/5.0 + np.sin(x) )**2.0/2.0
generator = mh_generator(log_density=log_prob,x_start=1.0)
tester = GaussianSteinTest(grad_log_prob,41)
selector = SampleSelector(generator, sample_size=1000,thinning=20,tester=tester, max_iterations=5)
data,converged = selector.points_from_stationary()
tester = GaussianSteinTest(grad_log_prob,41)
assert tester.compute_pvalue(data)>0.05
assert converged is True
def random_points_in_circle(n,xx,yy,rr):
"""
get n random points in a circle.
"""
rnd = random(size=(n,3))
t = TWOPI*rnd[:,0]
u = rnd[:,1:].sum(axis=1)
r = zeros(n,'float')
mask = u>1.
xmask = logical_not(mask)
r[mask] = 2.-u[mask]
r[xmask] = u[xmask]
xyp = reshape(rr*r,(n,1))*column_stack( (cos(t),sin(t)) )
dartsxy = xyp + array([xx,yy])
return dartsxy
def get_data(filename,headers,ph_units):
# Importation des données .DAT
dat_file = np.loadtxt("%s"%(filename),skiprows=headers,delimiter=',')
labels = ["freq", "amp", "pha", "amp_err", "pha_err"]
data = {l:dat_file[:,i] for (i,l) in enumerate(labels)}
if ph_units == "mrad":
data["pha"] = data["pha"]/1000 # mrad to rad
data["pha_err"] = data["pha_err"]/1000 # mrad to rad
if ph_units == "deg":
data["pha"] = np.radians(data["pha"]) # deg to rad
data["pha_err"] = np.radians(data["pha_err"]) # deg to rad
data["phase_range"] = abs(max(data["pha"])-min(data["pha"])) # Range of phase measurements (used in NRMS error calculation)
data["Z"] = data["amp"]*(np.cos(data["pha"]) + 1j*np.sin(data["pha"]))
EI = np.sqrt(((data["amp"]*np.cos(data["pha"])*data["pha_err"])**2)+(np.sin(data["pha"])*data["amp_err"])**2)
ER = np.sqrt(((data["amp"]*np.sin(data["pha"])*data["pha_err"])**2)+(np.cos(data["pha"])*data["amp_err"])**2)
data["Z_err"] = ER + 1j*EI
# Normalization of amplitude
data["Z_max"] = max(abs(data["Z"])) # Maximum amplitude
zn, zn_e = data["Z"]/data["Z_max"], data["Z_err"]/data["Z_max"] # Normalization of impedance by max amplitude
data["zn"] = np.array([zn.real, zn.imag]) # 2D array with first column = real values, second column = imag values
data["zn_err"] = np.array([zn_e.real, zn_e.imag]) # 2D array with first column = real values, second column = imag values
return data
def get_data(filename,headers,ph_units):
# Importation des données .DAT
dat_file = np.loadtxt("%s"%(filename),skiprows=headers,delimiter=',')
labels = ["freq", "amp", "pha", "amp_err", "pha_err"]
data = {l:dat_file[:,i] for (i,l) in enumerate(labels)}
if ph_units == "mrad":
data["pha"] = data["pha"]/1000 # mrad to rad
data["pha_err"] = data["pha_err"]/1000 # mrad to rad
if ph_units == "deg":
data["pha"] = np.radians(data["pha"]) # deg to rad
data["pha_err"] = np.radians(data["pha_err"]) # deg to rad
data["phase_range"] = abs(max(data["pha"])-min(data["pha"])) # Range of phase measurements (used in NRMS error calculation)
data["Z"] = data["amp"]*(np.cos(data["pha"]) + 1j*np.sin(data["pha"]))
EI = np.sqrt(((data["amp"]*np.cos(data["pha"])*data["pha_err"])**2)+(np.sin(data["pha"])*data["amp_err"])**2)
ER = np.sqrt(((data["amp"]*np.sin(data["pha"])*data["pha_err"])**2)+(np.cos(data["pha"])*data["amp_err"])**2)
data["Z_err"] = ER + 1j*EI
# Normalization of amplitude
data["Z_max"] = max(abs(data["Z"])) # Maximum amplitude
zn, zn_e = data["Z"]/data["Z_max"], data["Z_err"]/data["Z_max"] # Normalization of impedance by max amplitude
data["zn"] = np.array([zn.real, zn.imag]) # 2D array with first column = real values, second column = imag values
data["zn_err"] = np.array([zn_e.real, zn_e.imag]) # 2D array with first column = real values, second column = imag values
return data
def genSphCoords():
""" Generates cartesian (x,y,z) and spherical (theta, phi) coordinates of a sphere
Returns
-------
coords : named tuple
holds cartesian (x,y,z) and spherical (theta, phi) coordinates
"""
coords = namedtuple('coords', ['x', 'y', 'z', 'az', 'el'])
az = _np.linspace(0, 2 * _np.pi, 360)
el = _np.linspace(0, _np.pi, 181)
coords.x = _np.outer(_np.cos(az), _np.sin(el))
coords.y = _np.outer(_np.sin(az), _np.sin(el))
coords.z = _np.outer(_np.ones(360), _np.cos(el))
coords.el, coords.az = _np.meshgrid(_np.linspace(0, _np.pi, 181),
_np.linspace(0, 2 * _np.pi, 360))
return coords
def sph2cartMTX(vizMTX):
""" Converts the spherical vizMTX data to named tuple contaibubg .xs/.ys/.zs
Parameters
----------
vizMTX : array_like
[180 x 360] matrix that hold amplitude information over phi and theta
Returns
-------
V : named_tuple
Contains .xs, .ys, .zs cartesian coordinates
"""
rs = _np.abs(vizMTX.reshape((181, -1)).T)
coords = genSphCoords()
V = namedtuple('V', ['xs', 'ys', 'zs'])
V.xs = rs * _np.sin(coords.el) * _np.cos(coords.az)
V.ys = rs * _np.sin(coords.el) * _np.sin(coords.az)
V.zs = rs * _np.cos(coords.el)
return V
def convert_cof_mag2mass(t0, te, u0, alpha, s, q):
"""
function to convert from center of magnification to center of mass
coordinates. Note that this function is for illustration only. It has
not been tested and may have sign errors.
"""
if s <= 1.0:
return t0, u0
else:
delta = q / (1. + q) / s
delta_u0 = delta * np.sin(alpha * np.pi / 180.)
delta_tau = delta * np.cos(alpha * np.pi / 180.)
t0_prime = t0 + delta_tau * te
u0_prime = u0 + delta_u0
return t0_prime, u0_prime
#Define model parameters in CoMAGN system
def _B_1_function(self, z, B_0=None):
"""
calculate B_1(z) function defined in:
Gould A. 1994 ApJ 421L, 71 "Proper motions of MACHOs
http://adsabs.harvard.edu/abs/1994ApJ...421L..71G
Yoo J. et al. 2004 ApJ 603, 139 "OGLE-2003-BLG-262: Finite-Source
Effects from a Point-Mass Lens"
http://adsabs.harvard.edu/abs/2004ApJ...603..139Y
"""
if B_0 is None:
B_0 = self._B_0_function(z)
function = (lambda r, theta: r * np.sqrt(1.-r**2) /
self.parameters.rho / np.sqrt(r**2+zz**2-2.*r*zz*np.cos(theta)))
lim_0 = lambda x: 0
lim_1 = lambda x: 1
W_1 = 0. * z
for (i, zz) in enumerate(z):
W_1[i] = integrate.dblquad(function, 0., 2.*np.pi, lim_0, lim_1)[0]
W_1 /= np.pi
return B_0 - 1.5 * z * self.parameters.rho * W_1
edge_canny_get_orientation_sector.py 文件源码
项目：Image-Processing-Algorithms
作者: machinelearninggod
项目源码
文件源码
阅读 31
收藏 0
点赞 0
评论 0
def get_orientation_sector(dx,dy):
# rotate (dx,dy) by pi/8
rotation = np.array([[np.cos(np.pi/8), -np.sin(np.pi/8)],
[np.sin(np.pi/8), np.cos(np.pi/8)]])
rotated = np.dot(rotation, np.array([[dx], [dy]]))
if rotated[1] < 0:
rotated[0] = -rotated[0]
rotated[1] = -rotated[1]
s_theta = -1
if rotated[0] >= 0 and rotated[0] >= rotated[1]:
s_theta = 0
elif rotated[0] >= 0 and rotated[0] < rotated[1]:
s_theta = 1
elif rotated[0] < 0 and -rotated[0] < rotated[1]:
s_theta = 2
elif rotated[0] < 0 and -rotated[0] >= rotated[1]:
s_theta = 3
return s_theta
def to_radec(coords, xc=0, yc=0):
"""
Convert the generated coordinates to (ra, dec) (unit: degree).
xc, yc: the center coordinate (ra, dec)
"""
results = []
for r, theta in coords:
# FIXME: spherical algebra should be used!!!
dx = r * np.cos(theta*np.pi/180)
dy = r * np.sin(theta*np.pi/180)
x = xc + dx
y = yc + dy
results.append((x, y))
if len(results) == 1:
return results[0]
else:
return results
