The following are code examples for showing how to use . They are extracted from open source Python projects. You can vote up the examples you like or vote down the exmaples you don’t like. You can also save this page to your account.
Example 1
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)) Example 2
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 Example 3
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 Example 4
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 Example 5
def monotoneTFosc(f):
"""Maps [-inf,inf] to [-inf,inf] with different constants
for positive and negative part.
"""
if np.isscalar(f):
if f > 0.:
f = np.log(f) / 0.1
f = np.exp(f + 0.49 * (np.sin(f) + np.sin(0.79 * f))) ** 0.1
elif f < 0.:
f = np.log(-f) / 0.1
f = -np.exp(f + 0.49 * (np.sin(0.55 * f) + np.sin(0.31 * f))) ** 0.1
return f
else:
f = np.asarray(f)
g = f.copy()
idx = (f > 0)
g[idx] = np.log(f[idx]) / 0.1
g[idx] = np.exp(g[idx] + 0.49 * (np.sin(g[idx]) + np.sin(0.79 * g[idx])))**0.1
idx = (f < 0)
g[idx] = np.log(-f[idx]) / 0.1
g[idx] = -np.exp(g[idx] + 0.49 * (np.sin(0.55 * g[idx]) + np.sin(0.31 * g[idx])))**0.1
return g Example 6
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)) Example 7
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) Example 8
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 Example 9
def solveIter(mu,e):
"""__solvIter returns an iterative solution for Ek
Mk = Ek - e sin(Ek)
"""
thisStart = np.asarray(mu-1.01*e)
thisEnd = np.asarray(mu + 1.01*e)
bestGuess = np.zeros(mu.shape)
for i in range(5):
minErr = 10000*np.ones(mu.shape)
for j in range(5):
thisGuess = thisStart + j*(thisEnd-thisStart)/10.0
thisErr = np.asarray(abs(mu - thisGuess + e*np.sin(thisGuess)))
mask = thisErr<minErr
minErr[mask] = thisErr[mask]
bestGuess[mask] = thisGuess[mask]
# reset for next loop
thisRange = thisEnd - thisStart
thisStart = bestGuess - thisRange/10.0
thisEnd = bestGuess + thisRange/10.0
return(bestGuess) Example 10
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 Example 11
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 Example 12
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 Example 13
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 Example 14
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 Example 15
def mdst(x, odd=True):
""" Calculate modified discrete sine 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.sin, odd=odd) * numpy.sqrt(2) Example 16
def imdst(X, odd=True):
""" Calculate inverse modified discrete sine 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.sin, odd=odd) * numpy.sqrt(2) Example 17
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) Example 18
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) Example 19
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 Example 20
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 Example 21
def bandpass(self, rin, sin, rout, sout):
''' To create a band pass two circle images are created, one inverted
and pasted into dthe other'''
# if radius zero dont create the inner circle
if rin != 0:
self.create_circle_mask(rin, sin)
else:
self.data = np.zeros(self.data.shape)
# create the outer circle
bigcircle = deepcopy(self)
bigcircle.create_circle_mask(rout, sout)
bigcircle.invert()
# sum the two pictures
m = (self + bigcircle)
# limit fro 0 to 1 and invert
m.limit(1)
m.invert()
self.data = m.data Example 22
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
Example 23
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 Example 24
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 Example 25
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 Example 26
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 Example 27
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 Example 28
def _B_0_function(self, z):
"""
calculate B_0(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
"""
out = 4. * z / np.pi
function = lambda x: (1.-value**2*np.sin(x)**2)**.5
for (i, value) in enumerate(z):
if value < 1.:
out[i] *= ellipe(value*value)
else:
out[i] *= integrate.quad(function, 0., np.arcsin(1./value))[0]
return out Example 29
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 Example 30
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 Example 31
def plotArc(start_angle, stop_angle, radius, width, **kwargs):
""" write a docstring for this function"""
numsegments = 100
theta = np.radians(np.linspace(start_angle+90, stop_angle+90, numsegments))
centerx = 0
centery = 0
x1 = -np.cos(theta) * (radius)
y1 = np.sin(theta) * (radius)
stack1 = np.column_stack([x1, y1])
x2 = -np.cos(theta) * (radius + width)
y2 = np.sin(theta) * (radius + width)
stack2 = np.column_stack([np.flip(x2, axis=0), np.flip(y2,axis=0)])
#add the first values from the first set to close the polygon
np.append(stack2, [[x1[0],y1[0]]], axis=0)
arcArray = np.concatenate((stack1,stack2), axis=0)
return patches.Polygon(arcArray, True, **kwargs), ((x1, y1), (x2, y2)) Example 32
def ct2lg(dX, dY, dZ, lat, lon):
n = dX.size
R = np.zeros((3, 3, n))
R[0, 0, :] = -np.multiply(np.sin(np.deg2rad(lat)), np.cos(np.deg2rad(lon)))
R[0, 1, :] = -np.multiply(np.sin(np.deg2rad(lat)), np.sin(np.deg2rad(lon)))
R[0, 2, :] = np.cos(np.deg2rad(lat))
R[1, 0, :] = -np.sin(np.deg2rad(lon))
R[1, 1, :] = np.cos(np.deg2rad(lon))
R[1, 2, :] = np.zeros((1, n))
R[2, 0, :] = np.multiply(np.cos(np.deg2rad(lat)), np.cos(np.deg2rad(lon)))
R[2, 1, :] = np.multiply(np.cos(np.deg2rad(lat)), np.sin(np.deg2rad(lon)))
R[2, 2, :] = np.sin(np.deg2rad(lat))
dxdydz = np.column_stack((np.column_stack((dX, dY)), dZ))
RR = np.reshape(R[0, :, :], (3, n))
dx = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0)
RR = np.reshape(R[1, :, :], (3, n))
dy = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0)
RR = np.reshape(R[2, :, :], (3, n))
dz = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0)
return dx, dy, dz Example 33
def distance(self, lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1 = lon1*pi/180
lat1 = lat1*pi/180
lon2 = lon2*pi/180
lat2 = lat2*pi/180
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
c = 2 * np.arcsin(np.sqrt(a))
km = 6371 * c
return km Example 34
def distance(self, lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1 = lon1*pi/180
lat1 = lat1*pi/180
lon2 = lon2*pi/180
lat2 = lat2*pi/180
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = numpy.sin(dlat/2)**2 + numpy.cos(lat1) * numpy.cos(lat2) * numpy.sin(dlon/2)**2
c = 2 * numpy.arcsin(numpy.sqrt(a))
km = 6371 * c
return km Example 35
def ct2lg(dX, dY, dZ, lat, lon):
n = dX.size
R = np.zeros((3, 3, n))
R[0, 0, :] = -np.multiply(np.sin(np.deg2rad(lat)), np.cos(np.deg2rad(lon)))
R[0, 1, :] = -np.multiply(np.sin(np.deg2rad(lat)), np.sin(np.deg2rad(lon)))
R[0, 2, :] = np.cos(np.deg2rad(lat))
R[1, 0, :] = -np.sin(np.deg2rad(lon))
R[1, 1, :] = np.cos(np.deg2rad(lon))
R[1, 2, :] = np.zeros((1, n))
R[2, 0, :] = np.multiply(np.cos(np.deg2rad(lat)), np.cos(np.deg2rad(lon)))
R[2, 1, :] = np.multiply(np.cos(np.deg2rad(lat)), np.sin(np.deg2rad(lon)))
R[2, 2, :] = np.sin(np.deg2rad(lat))
dxdydz = np.column_stack((np.column_stack((dX, dY)), dZ))
RR = np.reshape(R[0, :, :], (3, n))
dx = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0)
RR = np.reshape(R[1, :, :], (3, n))
dy = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0)
RR = np.reshape(R[2, :, :], (3, n))
dz = np.sum(np.multiply(RR, dxdydz.transpose()), axis=0)
return dx, dy, dz Example 36
def distance(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1 = lon1*pi/180
lat1 = lat1*pi/180
lon2 = lon2*pi/180
lat2 = lat2*pi/180
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2
c = 2 * np.arcsin(np.sqrt(a))
km = 6371 * c
return km Example 37
def todictionary(self, time_series=False):
# convert the ETM adjustment into a dirtionary
# optionally, output the whole time series as well
# start with the parameters
etm = dict()
etm['Linear'] = {'tref': self.Linear.tref, 'params': self.Linear.values.tolist()}
etm['Jumps'] = [{'type':jump.type, 'year': jump.year, 'a': jump.a.tolist(), 'b': jump.b.tolist(), 'T': jump.T} for jump in self.Jumps.table]
etm['Periodic'] = {'frequencies': self.Periodic.frequencies, 'sin': self.Periodic.sin.tolist(), 'cos': self.Periodic.cos.tolist()}
if time_series:
ts = dict()
ts['t'] = self.ppp_soln.t.tolist()
ts['x'] = self.ppp_soln.x.tolist()
ts['y'] = self.ppp_soln.y.tolist()
ts['z'] = self.ppp_soln.z.tolist()
etm['time_series'] = ts
return etm Example 38
def ct2lg(self, dX, dY, dZ, lat, lon):
n = dX.size
R = numpy.zeros((3, 3, n))
R[0, 0, :] = -numpy.multiply(numpy.sin(numpy.deg2rad(lat)), numpy.cos(numpy.deg2rad(lon)))
R[0, 1, :] = -numpy.multiply(numpy.sin(numpy.deg2rad(lat)), numpy.sin(numpy.deg2rad(lon)))
R[0, 2, :] = numpy.cos(numpy.deg2rad(lat))
R[1, 0, :] = -numpy.sin(numpy.deg2rad(lon))
R[1, 1, :] = numpy.cos(numpy.deg2rad(lon))
R[1, 2, :] = numpy.zeros((1, n))
R[2, 0, :] = numpy.multiply(numpy.cos(numpy.deg2rad(lat)), numpy.cos(numpy.deg2rad(lon)))
R[2, 1, :] = numpy.multiply(numpy.cos(numpy.deg2rad(lat)), numpy.sin(numpy.deg2rad(lon)))
R[2, 2, :] = numpy.sin(numpy.deg2rad(lat))
dxdydz = numpy.column_stack((numpy.column_stack((dX, dY)), dZ))
RR = numpy.reshape(R[0, :, :], (3, n))
dx = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0)
RR = numpy.reshape(R[1, :, :], (3, n))
dy = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0)
RR = numpy.reshape(R[2, :, :], (3, n))
dz = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0)
return dx, dy, dz Example 39
def distance(self, lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1 = lon1*pi/180
lat1 = lat1*pi/180
lon2 = lon2*pi/180
lat2 = lat2*pi/180
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = numpy.sin(dlat/2)**2 + numpy.cos(lat1) * numpy.cos(lat2) * numpy.sin(dlon/2)**2
c = 2 * numpy.arcsin(numpy.sqrt(a))
km = 6371 * c
return km Example 40
def distance(self, lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1 = lon1*pi/180
lat1 = lat1*pi/180
lon2 = lon2*pi/180
lat2 = lat2*pi/180
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = numpy.sin(dlat/2)**2 + numpy.cos(lat1) * numpy.cos(lat2) * numpy.sin(dlon/2)**2
c = 2 * numpy.arcsin(numpy.sqrt(a))
km = 6371 * c
return km Example 41
def ct2lg(self, dX, dY, dZ, lat, lon):
n = dX.size
R = numpy.zeros((3, 3, n))
R[0, 0, :] = -numpy.multiply(numpy.sin(numpy.deg2rad(lat)), numpy.cos(numpy.deg2rad(lon)))
R[0, 1, :] = -numpy.multiply(numpy.sin(numpy.deg2rad(lat)), numpy.sin(numpy.deg2rad(lon)))
R[0, 2, :] = numpy.cos(numpy.deg2rad(lat))
R[1, 0, :] = -numpy.sin(numpy.deg2rad(lon))
R[1, 1, :] = numpy.cos(numpy.deg2rad(lon))
R[1, 2, :] = numpy.zeros((1, n))
R[2, 0, :] = numpy.multiply(numpy.cos(numpy.deg2rad(lat)), numpy.cos(numpy.deg2rad(lon)))
R[2, 1, :] = numpy.multiply(numpy.cos(numpy.deg2rad(lat)), numpy.sin(numpy.deg2rad(lon)))
R[2, 2, :] = numpy.sin(numpy.deg2rad(lat))
dxdydz = numpy.column_stack((numpy.column_stack((dX, dY)), dZ))
RR = numpy.reshape(R[0, :, :], (3, n))
dx = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0)
RR = numpy.reshape(R[1, :, :], (3, n))
dy = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0)
RR = numpy.reshape(R[2, :, :], (3, n))
dz = numpy.sum(numpy.multiply(RR, dxdydz.transpose()), axis=0)
return dx, dy, dz Example 42
def orthogonalization_matrix(lengths, angles):
"""Return orthogonalization matrix for crystallographic cell coordinates.
Angles are expected in degrees.
The de-orthogonalization matrix is the inverse.
>>> O = orthogonalization_matrix([10, 10, 10], [90, 90, 90])
>>> numpy.allclose(O[:3, :3], numpy.identity(3, float) * 10)
True
>>> O = orthogonalization_matrix([9.8, 12.0, 15.5], [87.2, 80.7, 69.7])
>>> numpy.allclose(numpy.sum(O), 43.063229)
True
"""
a, b, c = lengths
angles = numpy.radians(angles)
sina, sinb, _ = numpy.sin(angles)
cosa, cosb, cosg = numpy.cos(angles)
co = (cosa * cosb - cosg) / (sina * sinb)
return numpy.array([
[ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0],
[-a*sinb*co, b*sina, 0.0, 0.0],
[ a*cosb, b*cosa, c, 0.0],
[ 0.0, 0.0, 0.0, 1.0]]) Example 43
def plotPlaceTxRxSphereXY(Ax,xtx,ytx,xrx,yrx,x0,y0,a):
Xlim = Ax.get_xlim()
Ylim = Ax.get_ylim()
FS = 20
Ax.scatter(xtx,ytx,s=100,color='k')
Ax.text(xtx-0.75,ytx+1.5,'$\mathbf{Tx}$',fontsize=FS+6)
Ax.scatter(xrx,yrx,s=100,color='k')
Ax.text(xrx-0.75,yrx-4,'$\mathbf{Rx}$',fontsize=FS+6)
xs = x0 + a*np.cos(np.linspace(0,2*np.pi,41))
ys = y0 + a*np.sin(np.linspace(0,2*np.pi,41))
Ax.plot(xs,ys,ls=':',color='k',linewidth=3)
Ax.set_xbound(Xlim)
Ax.set_ybound(Ylim)
return Ax Example 44
def calc_IndCurrent_cos_range(self,f,t):
"""Induced current over a range of times"""
Bpx = self.Bpx
Bpz = self.Bpz
a2 = self.a2
azm = np.pi*self.azm/180.
R = self.R
L = self.L
w = 2*np.pi*f
Ax = np.pi*a2**2*np.sin(azm)
Az = np.pi*a2**2*np.cos(azm)
Phi = (Ax*Bpx + Az*Bpz)
phi = np.arctan(R/(w*L))-np.pi # This is the phase and not phase lag
Is = -(w*Phi/(R*np.sin(phi) + w*L*np.cos(phi)))*np.cos(w*t + phi)
Ire = -(w*Phi/(R*np.sin(phi) + w*L*np.cos(phi)))*np.cos(w*t)*np.cos(phi)
Iim = (w*Phi/(R*np.sin(phi) + w*L*np.cos(phi)))*np.sin(w*t)*np.sin(phi)
return Ire,Iim,Is,phi Example 45
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 Example 46
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) Example 47
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 Example 48
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
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = Example 49
def residual(r,theta,u,d):
out = np.empty_like(u)
out[0] = (2*np.sin(theta)*r*d(u[0],1,0)
+ r*r*np.sin(theta)*d(u[0],2,0)
+ np.cos(theta)*d(u[0],0,1)
+ np.sin(theta)*d(u[1],0,2))
out[1] = (2*np.sin(theta)*r*d(u[1],1,0)
+ r*r*np.sin(theta)*d(u[1],2,0)
+ np.cos(theta)*d(u[1],0,1)
+ np.sin(theta)*d(u[1],0,2))
return out Example 50
def residual(r,theta,u,d):
u = u[0]
out = (2*np.sin(theta)*r*d(u,1,0)
+ r*r*np.sin(theta)*d(u,2,0)
+ np.cos(theta)*d(u,0,1)
+ np.sin(theta)*d(u,0,2))
out = out.reshape(tuple([1]) + out.shape)
return out 