Python numpy 模块，c_() 实例源码
我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用numpy.c_()。
def array2PIL(arr, size):
mode = 'RGBA'
arr = arr.reshape(arr.shape[0]*arr.shape[1], arr.shape[2])
if len(arr[0]) == 3:
arr = numpy.c_[arr, 255*numpy.ones((len(arr),1), numpy.uint8)]
return Image.frombuffer(mode, size, arr.tostring(), 'raw', mode, 0, 1)
def best_fit_plane(self):
"""Fits a plane to the point cloud using least squares.
Returns
-------
:obj:`tuple` of :obj:`numpy.ndarray` of float
A normal vector to and point in the fitted plane.
"""
X = np.c_[self.x_coords, self.y_coords, np.ones(self.num_points)]
y = self.z_coords
A = X.T.dot(X)
b = X.T.dot(y)
w = np.linalg.inv(A).dot(b)
n = np.array([w[0], w[1], -1])
n = n / np.linalg.norm(n)
n = Direction(n, self._frame)
x0 = self.mean()
return n, x0
def predict_bbox_regressor(model, feat, ex_boxes):
if ex_boxes.size == 0:
return np.array([]).reshape(-1, 4)
# predict regression targets
Y = np.dot(feat, model.Beta[:-1]) + model.Beta[-1]
# invert transformation
Y = dot(Y, model.T_inv)
# read out prediction
dst_size = Y[:, 2:]
dst_ctr = Y[:, 2:]
src_size = ex_boxes[:, 2:]
src_ctr = ex_boxes[:, :2] + 0.5 * src_size
pred_size = np.exp(dst_size) * src_size
pred_ctr = dst_ctr * src_ctr + src_ctr
pred = np.c_[pred_ctr - 0.5 * pred_size, pred_size]
return pred
def get_examples(bbox, gt):
# compute overlap ratio
O = overlap_ratio(bbox, gt)
# compute answer
src_size = bbox[:, 2:]
src_ctr = bbox[:, :2] + 0.5 * src_size
gt_size = gt[2:]
gt_ctr = gt[:2] + 0.5 * gt_size
dst_size = np.log(gt_size / src_size)
dst_ctr = (gt_ctr - src_ctr) * 1.0 / src_ctr
Y = np.c_[dst_ctr, dst_size]
return Y, O
def computeRotMatrix(self,Phi=False):
#######################################
# COMPUTE ROTATION MATRIX SUCH THAT m(t) = A*L(t)*A'*Hp
# Default set such that phi1,phi2 = 0 is UXO pointed towards North
if Phi is False:
phi1 = np.radians(self.phi[0])
phi2 = np.radians(self.phi[1])
phi3 = np.radians(self.phi[2])
else:
phi1 = np.radians(Phi[0]) # Roll (CCW)
phi2 = np.radians(Phi[1]) # Inclination (+ve is nose pointing down)
phi3 = np.radians(Phi[2]) # Declination (degrees CW from North)
# A1 = np.r_[np.c_[np.cos(phi1),-np.sin(phi1),0.],np.c_[np.sin(phi1),np.cos(phi1),0.],np.c_[0.,0.,1.]] # CCW Rotation about z-axis
# A2 = np.r_[np.c_[1.,0.,0.],np.c_[0.,np.cos(phi2),np.sin(phi2)],np.c_[0.,-np.sin(phi2),np.cos(phi2)]] # CW Rotation about x-axis (rotates towards North)
# A3 = np.r_[np.c_[np.cos(phi3),-np.sin(phi3),0.],np.c_[np.sin(phi3),np.cos(phi3),0.],np.c_[0.,0.,1.]] # CCW Rotation about z-axis (direction of head of object)
A1 = np.r_[np.c_[np.cos(phi1),np.sin(phi1),0.],np.c_[-np.sin(phi1),np.cos(phi1),0.],np.c_[0.,0.,1.]] # CW Rotation about z-axis
A2 = np.r_[np.c_[1.,0.,0.],np.c_[0.,np.cos(phi2),np.sin(phi2)],np.c_[0.,-np.sin(phi2),np.cos(phi2)]] # CW Rotation about x-axis (rotates towards North)
A3 = np.r_[np.c_[np.cos(phi3),np.sin(phi3),0.],np.c_[-np.sin(phi3),np.cos(phi3),0.],np.c_[0.,0.,1.]] # CW Rotation about z-axis (direction of head of object)
return np.dot(A3,np.dot(A2,A1))
def defineSensorLoc(self,XYZ):
#############################################
# DEFINE TRANSMITTER AND RECEIVER LOCATIONS
#
# XYZ: N X 3 array containing transmitter center locations
# **NOTE** for this sensor, we know where the receivers are relative to transmitters
self.TxLoc = XYZ
dx,dy = np.meshgrid([-0.8,-0.4,0.,0.4,0.8],[-0.8,-0.4,0.,0.4,0.8])
dx = mkvc(dx)
dy = mkvc(dy)
N = np.shape(XYZ)[0]
X = np.kron(XYZ[:,0],np.ones((25))) + np.kron(np.ones((N)),dx)
Y = np.kron(XYZ[:,1],np.ones((25))) + np.kron(np.ones((N)),dy)
Z = np.kron(XYZ[:,2],np.ones((25)))
self.RxLoc = np.c_[X,Y,Z]
self.TxID = np.kron(np.arange(1,np.shape(XYZ)[0]+1),np.ones((25)))
self.RxComp = np.kron(3*np.ones(np.shape(XYZ)[0]),np.ones((25)))
def updatePolarizations(self,r0,UB):
# Set operator and solution array
Hp = self.computeHp(r0=r0)
Brx = self.computeBrx(r0=r0)
P = self.computeP(Hp,Brx)
dunc = self.dunc
dobs = self.dobs
K = np.shape(dobs)[1]
q = np.zeros((6,K))
lb = np.zeros(6)
ub = UB*np.ones(6)
for kk in range(0,K):
LHS = P/np.c_[dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk]]
RHS = dobs[:,kk]/dunc[:,kk]
Sol = op.lsq_linear(LHS,RHS,bounds=(lb,ub),tol=1e-5)
q[:,kk] = Sol.x
self.q = q
def updatePolarizations(self,r0,UB):
# Set operator and solution array
Hp = self.computeHp(r0=r0)
Brx = self.computeBrx(r0=r0)
P = self.computeP(Hp,Brx)
dunc = self.dunc
dobs = self.dobs
K = np.shape(dobs)[1]
q = np.zeros((6,K))
lb = np.zeros(6)
ub = UB*np.ones(6)
for kk in range(0,K):
LHS = P/np.c_[dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk],dunc[:,kk]]
RHS = dobs[:,kk]/dunc[:,kk]
Sol = op.lsq_linear(LHS,RHS,bounds=(lb,ub),tol=1e-7)
q[:,kk] = Sol.x
self.q = q
def analytic_infinite_wire(obsloc,wireloc,orientation,I=1.):
"""
Compute the response of an infinite wire with orientation 'orientation'
and current I at the obsvervation locations obsloc
Output:
B: magnetic field [Bx,By,Bz]
"""
n,d = obsloc.shape
t,d = wireloc.shape
d = np.sqrt(np.dot(obsloc**2.,np.ones([d,t]))+np.dot(np.ones([n,d]),(wireloc.T)**2.)
- 2.*np.dot(obsloc,wireloc.T))
distr = np.amin(d, axis=1, keepdims = True)
idxmind = d.argmin(axis=1)
r = obsloc - wireloc[idxmind]
orient = np.c_[[orientation for i in range(obsloc.shape[0])]]
B = (mu_0*I)/(2*np.pi*(distr**2.))*np.cross(orientation,r)
return B
def predict_with_glm(X, y, model):
""" Predict number of mutation with GLM.
Args:
X (np.array): feature matrix.
y (pd.df): response.
model (dict): model meta-data.
Returns:
np.array: array of predictions.
"""
# Add const. to X
X = np.c_[X, np.ones(X.shape[0])]
if model['model_name'] == 'Binomial':
pred = np.array(model['model'].predict(X) * y.length * y.N)
elif model['model_name'] == 'NegativeBinomial':
pred = np.array(model['model'].predict(X, exposure=(y.length * y.N).values + 1))
else:
sys.stderr.write('Wrong model name in model info: {}. Need Binomial or NegativeBinomial.'.format(model['model_name']))
sys.exit(1)
return pred
def normals(self):
"""Face Normals
:rtype: numpy.array
:return: normals, (sum(nF), dim)
"""
if self.dim == 2:
nX = np.c_[
np.ones(self.nFx), np.zeros(self.nFx)
]
nY = np.c_[
np.zeros(self.nFy), np.ones(self.nFy)
]
return np.r_[nX, nY]
elif self.dim == 3:
nX = np.c_[
np.ones(self.nFx), np.zeros(self.nFx), np.zeros(self.nFx)
]
nY = np.c_[
np.zeros(self.nFy), np.ones(self.nFy), np.zeros(self.nFy)
]
nZ = np.c_[
np.zeros(self.nFz), np.zeros(self.nFz), np.ones(self.nFz)
]
return np.r_[nX, nY, nZ]
def tangents(self):
"""Edge Tangents
:rtype: numpy.array
:return: normals, (sum(nE), dim)
"""
if self.dim == 2:
tX = np.c_[
np.ones(self.nEx), np.zeros(self.nEx)
]
tY = np.c_[
np.zeros(self.nEy), np.ones(self.nEy)
]
return np.r_[tX, tY]
elif self.dim == 3:
tX = np.c_[
np.ones(self.nEx), np.zeros(self.nEx), np.zeros(self.nEx)
]
tY = np.c_[
np.zeros(self.nEy), np.ones(self.nEy), np.zeros(self.nEy)
]
tZ = np.c_[
np.zeros(self.nEz), np.zeros(self.nEz), np.ones(self.nEz)
]
return np.r_[tX, tY, tZ]
def test_invPropertyTensor2D(self):
M = discretize.TensorMesh([6, 6])
a1 = np.random.rand(M.nC)
a2 = np.random.rand(M.nC)
a3 = np.random.rand(M.nC)
prop1 = a1
prop2 = np.c_[a1, a2]
prop3 = np.c_[a1, a2, a3]
for prop in [4, prop1, prop2, prop3]:
b = invPropertyTensor(M, prop)
A = makePropertyTensor(M, prop)
B1 = makePropertyTensor(M, b)
B2 = invPropertyTensor(M, prop, returnMatrix=True)
Z = B1*A - sp.identity(M.nC*2)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
Z = B2*A - sp.identity(M.nC*2)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
def test_TensorType3D(self):
M = discretize.TensorMesh([6, 6, 7])
a1 = np.random.rand(M.nC)
a2 = np.random.rand(M.nC)
a3 = np.random.rand(M.nC)
a4 = np.random.rand(M.nC)
a5 = np.random.rand(M.nC)
a6 = np.random.rand(M.nC)
prop1 = a1
prop2 = np.c_[a1, a2, a3]
prop3 = np.c_[a1, a2, a3, a4, a5, a6]
for ii, prop in enumerate([4, prop1, prop2, prop3]):
self.assertTrue(TensorType(M, prop) == ii)
self.assertRaises(Exception, TensorType, M, np.c_[a1, a2, a3, a3])
self.assertTrue(TensorType(M, None) == -1)
def test_invPropertyTensor3D(self):
M = discretize.TensorMesh([6, 6, 6])
a1 = np.random.rand(M.nC)
a2 = np.random.rand(M.nC)
a3 = np.random.rand(M.nC)
a4 = np.random.rand(M.nC)
a5 = np.random.rand(M.nC)
a6 = np.random.rand(M.nC)
prop1 = a1
prop2 = np.c_[a1, a2, a3]
prop3 = np.c_[a1, a2, a3, a4, a5, a6]
for prop in [4, prop1, prop2, prop3]:
b = invPropertyTensor(M, prop)
A = makePropertyTensor(M, prop)
B1 = makePropertyTensor(M, b)
B2 = invPropertyTensor(M, prop, returnMatrix=True)
Z = B1*A - sp.identity(M.nC*3)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
Z = B2*A - sp.identity(M.nC*3)
self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < TOL)
def getError(self):
funR = lambda r, z: np.sin(2.*np.pi*r)
funZ = lambda r, z: np.sin(2.*np.pi*z)
sol = lambda r, t, z: (2*np.pi*r*np.cos(2*np.pi*r) + np.sin(2*np.pi*r))/r + 2*np.pi*np.cos(2*np.pi*z)
Fc = cylF2(self.M, funR, funZ)
Fc = np.c_[Fc[:, 0], np.zeros(self.M.nF), Fc[:, 1]]
F = self.M.projectFaceVector(Fc)
divF = self.M.faceDiv.dot(F)
divF_ana = call3(sol, self.M.gridCC)
err = np.linalg.norm((divF-divF_ana), np.inf)
return err
def getError(self):
funR = lambda r, z: np.sin(2.*np.pi*z) * np.cos(np.pi*r)
funZ = lambda r, z: np.sin(3.*np.pi*z) * np.cos(2.*np.pi*r)
Fc = cylF2(self.M, funR, funZ)
Fc = np.c_[Fc[:, 0], np.zeros(self.M.nF), Fc[:, 1]]
F = self.M.projectFaceVector(Fc)
aveF = self.M.aveF2CCV * F
aveF_anaR = funR(self.M.gridCC[:, 0], self.M.gridCC[:, 2])
aveF_anaZ = funZ(self.M.gridCC[:, 0], self.M.gridCC[:, 2])
aveF_ana = np.hstack([aveF_anaR, aveF_anaZ])
err = np.linalg.norm((aveF-aveF_ana), np.inf)
return err
def surface_points(self, grid_basis=True):
"""Returns the points on the surface.
Parameters
----------
grid_basis : bool
If False, the surface points are transformed to the world frame.
If True (default), the surface points are left in grid coordinates.
Returns
-------
:obj:`tuple` of :obj:`numpy.ndarray` of int, :obj:`numpy.ndarray` of float
The points on the surface and the signed distances at those points.
"""
surface_points = np.where(np.abs(self.data_) < self.surface_thresh_)
x = surface_points[0]
y = surface_points[1]
z = surface_points[2]
surface_points = np.c_[x, np.c_[y, z]]
surface_vals = self.data_[surface_points[:,0], surface_points[:,1], surface_points[:,2]]
if not grid_basis:
surface_points = self.transform_pt_grid_to_obj(surface_points.T)
surface_points = surface_points.T
return surface_points, surface_vals
def _check_freq(f):
"""Check the frequency definition."""
f = np.atleast_2d(np.asarray(f))
#
if len(f.reshape(-1)) == 1:
raise ValueError("The length of f should at least be 2.")
elif 2 in f.shape: # f of shape (N, 2) or (2, N)
if f.shape[1] is not 2:
f = f.T
elif np.squeeze(f).shape == (4,): # (fstart, fend, fwidth, fstep)
f = _pair_vectors(*tuple(np.squeeze(f)))
else: # Sequential
f = f.reshape(-1)
f.sort()
f = np.c_[f[0:-1], f[1::]]
return f
def convert_mask_to_locations(mask):
""" Return the converted Cartesian mask as sampling locations.
Parameters
----------
mask: np.ndarray, {0,1}
2D matrix, not necessarly a square matrix.
Returns
-------
samples_locations: np.ndarray
list of the samples between [-0.5, 0.5[.
"""
row, col = np.where(mask == 1)
row = row.astype("float") / mask.shape[0] - 0.5
col = col.astype("float") / mask.shape[1] - 0.5
return np.c_[row, col]
def _test_double_optimization():
"""Test double optimization on a simple example."""
# A simple sparse-sum function
X = [[1, 2], [3, 4], [5, 6]]
y = [sum(x) for x in X]
T = [[7, 8], [9, 10], [2, 1]]
# noisy variables
np.random.seed(0)
X = np.c_[X, np.random.random((3, 100))]
T = np.c_[T, np.random.random((3, 100))]
# Select the first 2 variables and calculate a linear model on them
dstep = DoubleStepEstimator(Lasso(tau=1.0), RidgeRegression(mu=0.0)).train(X, y)
# Coefficients
lasso = dstep.selector
ridge = dstep.estimator
assert_array_almost_equal([0.90635646, 0.90635646], lasso.beta[:2])
assert_array_almost_equal([1.0, 1.0], ridge.beta)
assert_array_almost_equal([1.0, 1.0], dstep.beta[:2])
# Prediction
y_ = dstep.predict(T)
assert_array_almost_equal([15., 19., 3.], y_)
def generate_hills(width, height, nhills):
'''
@param width float, terrain width
@param height float, terrain height
@param nhills int, #hills to gen. #hills actually generted is sqrt(nhills)^2
'''
# setup coordinate grid
xmin, xmax = -width/2.0, width/2.0
ymin, ymax = -height/2.0, height/2.0
x, y = np.mgrid[xmin:xmax:STEP, ymin:ymax:STEP]
pos = np.empty(x.shape + (2,))
pos[:, :, 0] = x; pos[:, :, 1] = y
# generate hilltops
xm, ym = np.mgrid[xmin:xmax:width/np.sqrt(nhills), ymin:ymax:height/np.sqrt(nhills)]
mu = np.c_[xm.flat, ym.flat]
sigma = float(width*height)/(nhills*8)
for i in range(mu.shape[0]):
mu[i] = multivariate_normal.rvs(mean=mu[i], cov=sigma)
# generate hills
sigma = sigma + sigma*np.random.rand(mu.shape[0])
rvs = [ multivariate_normal(mu[i,:], cov=sigma[i]) for i in range(mu.shape[0]) ]
hfield = np.max([ rv.pdf(pos) for rv in rvs ], axis=0)
return x, y, hfield
def generate_data(sample_size=200, pd=[[0.4, 0.4], [0.1, 0.1]]):
pd = np.array(pd)
pd /= pd.sum()
offset = 50
bins = np.r_[np.zeros((1,)), np.cumsum(pd)]
bin_counts = np.histogram(np.random.rand(sample_size), bins)[0]
data = np.empty((0, 2))
targets = []
for ((i, j), p), count in zip(np.ndenumerate(pd), bin_counts):
xs = np.random.uniform(low=0.0, high=50.0, size=count) + j * offset
ys = np.random.uniform(low=0.0, high=50.0, size=count) + -i * offset
data = np.vstack((data, np.c_[xs, ys]))
if i == j:
targets.extend([1] * count)
else:
targets.extend([-1] * count)
return np.c_[data, targets]
def show_classification_areas(X, Y, lr):
x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5
y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5
xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.02), np.arange(y_min, y_max, 0.02))
Z = lr.predict(np.c_[xx.ravel(), yy.ravel()])
Z = Z.reshape(xx.shape)
plt.figure(1, figsize=(30, 25))
plt.pcolormesh(xx, yy, Z, cmap=plt.cm.Pastel1)
# Plot also the training points
plt.scatter(X[:, 0], X[:, 1], c=np.abs(Y - 1), edgecolors='k', cmap=plt.cm.coolwarm)
plt.xlabel('X')
plt.ylabel('Y')
plt.xlim(xx.min(), xx.max())
plt.ylim(yy.min(), yy.max())
plt.xticks(())
plt.yticks(())
plt.show()
def plot_decision_boundary(X, Y, model):
# X - some data in 2dimensional np.array
x_min, x_max = X[:, 0].min() - 1, X[:, 0].max() + 1
y_min, y_max = X[:, 1].min() - 1, X[:, 1].max() + 1
xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.01),
np.arange(y_min, y_max, 0.01))
# here "model" is your model's prediction (classification) function
Z = model(np.c_[xx.ravel(), yy.ravel()])
# Put the result into a color plot
Z = Z.reshape(xx.shape)
plt.contourf(xx, yy, Z, cmap=plt.cm.Paired)
plt.axis('off')
for i in x:
print i
# Plot also the training points
plt.scatter(X[:, 0], X[:, 1], c=Y, cmap=plt.cm.Paired)
#???????
def get_double_integrator(m=1000, b=50, d=1):
N = 2
sys = LTISystem(
np.c_[[0, 1], [1, -b/m]], # A
np.r_[0, d/m], # B
# np.r_[0, 1], # C
)
def ref_func(*args):
if len(args) == 1:
x = np.zeros(N)
else:
x = args[1]
return np.r_[d/m, 0]-x
ref = SystemFromCallable(ref_func, N, N)
return sys, ref
def get_electromechanical(b=1, R=1, L=1, K=np.pi/5, M=1):
# TODO: determine good reference and/or initial_condition
# TODO: determine good default values for b, R, L, M
N = 3
sys = LTISystem(
np.c_[ # A
[0, 0, 0],
[1, -b/M, -K/L],
[0, K/M, -R/L]
],
np.r_[0, 0, 1/L], # B
# np.r_[1, 0, 0], # C
)
sys.initial_condition = np.ones(N)
def ref_func(*args):
if len(args) == 1:
x = np.zeros(N)
else:
x = args[1]
return np.r_[0, 0, 0]-x
ref = SystemFromCallable(ref_func, N, N)
return sys, ref
def get_cart_pendulum(m=1, M=3, L=0.5, g=9.81, pedant=False):
N = 4
sys = LTISystem(
np.c_[ # A
[0, 0, 0, 0],
[1, 0, 0, 0],
[0, m*g/M, 0, (-1)**(pedant)*(m+M)*g/(M*L)],
[0, 0, 1, 0]
],
np.r_[0, 1/M, 0, 1/(M*L)], # B
# np.r_[1, 0, 1, 0] # C
)
sys.initial_condition = np.r_[0, 0, np.pi/3, 0]
def ref_func(*args):
if len(args) == 1:
x = np.zeros(N)
else:
x = args[1]
return np.zeros(N)-x
ref = SystemFromCallable(ref_func, N, N)
return sys, ref
def update_equation_function(self, *args):
event_var = self.event_variable_equation_function(*args)
if self.condition_idx is None:
self.condition_idx = np.where(np.all(np.r_[
np.c_[[[True]], event_var >= self.event_bounds],
np.c_[event_var <= self.event_bounds, [[True]]]
], axis=0))[0][0]
else:
sq_dist = (event_var - self.event_bounds)**2
crossed_root_idx = np.where(sq_dist == np.min(sq_dist))[1][0]
if crossed_root_idx == self.condition_idx:
self.condition_idx += 1
elif crossed_root_idx == self.condition_idx-1:
self.condition_idx -= 1
else:
warnings.warn("SwitchedSystem did not cross a neighboring " +
"boundary. This may indicate an integration " +
"error. Continuing without updating " +
"condition_idx", UserWarning)
return self.state_update_equation_function(*args)
def plane2xyz(center, ij, plane):
"""
converts image pixel indices to xyz on the PLANE.
center : 2-tuple
ij : nx2 int array
plane : 4-tuple
return nx3 array.
"""
ij = np.atleast_2d(ij)
n = ij.shape[0]
ij = ij.astype('float')
xy_ray = (ij-center[None,:]) / DepthCamera.f
z = -plane[2]/(xy_ray.dot(plane[:2])+plane[3])
xyz = np.c_[xy_ray, np.ones(n)] * z[:,None]
return xyz
def viz_textbb(fignum,text_im, bb_list,alpha=1.0):
"""
text_im : image containing text
bb_list : list of 2x4xn_i boundinb-box matrices
"""
plt.close(fignum)
plt.figure(fignum)
plt.imshow(text_im)
plt.hold(True)
H,W = text_im.shape[:2]
for i in xrange(len(bb_list)):
bbs = bb_list[i]
ni = bbs.shape[-1]
for j in xrange(ni):
bb = bbs[:,:,j]
bb = np.c_[bb,bb[:,0]]
plt.plot(bb[0,:], bb[1,:], 'r', linewidth=2, alpha=alpha)
plt.gca().set_xlim([0,W-1])
plt.gca().set_ylim([H-1,0])
plt.show(block=False)
def bbox(self, image_coords=True):
"""
Return a 3 by 2 matrix, corresponding to the bounding box of the
annotation within the scan. If `scan_slice` is a numpy array
containing aslice of the scan, each slice of the annotation is
contained within the box:
bbox[0,0]:bbox[0,1]+1, bbox[1,0]:bbox[1,1]+1
If `image_coords` is `False` then each annotation slice is
instead contained within:
bbox[1,0]:bbox[1,1]+1, bbox[0,0]:bbox[0,1]+1
The last row of `bbox` give the inclusive lower and upper
bounds of the `image_z_position`.
"""
matrix = self.contours_to_matrix()
bbox = np.c_[matrix.min(axis=0), matrix.max(axis=0)]
return bbox if not image_coords else bbox[[1,0,2]]
def plot_and_save(x, func, xvv_inst, plot=False, fname=False):
mat = x.T
for m in range(xvv_inst.nsites):
for n in range(m+1):
if fname:
mat = np.c_[mat, func[:,m,n].T]
if plot:
plt.plot(x, func[:,m,n],
label='{}-{}'.format(xvv_inst.atom_names[m],
xvv_inst.atom_names[n]))
if fname:
np.savetxt(fname, mat)
if plot:
plt.legend()
plt.savefig('graph.png', dpi=300)
plt.show()
def read_csv(filename, skip_lines=0):
csvfile = file(filename, 'rb')
reader = csv.reader(csvfile)
data = np.empty(0, dtype=object)
last_count = np.NAN
for line in reader:
if skip_lines > 0:
skip_lines = skip_lines - 1
continue
if data.size > 0:
if len(line) != last_count:
raise Exception('unequal columes found')
data = np.c_[data, line]
last_count = len(line)
else:
data = np.array(line, dtype=object)
data = data.reshape(len(data), 1)
last_count = len(line)
csvfile.close()
return data.T
def print_evaluation_result(clf, bags_test, args):
pred_score = np.array([clf(B.data()) for B in bags_test])
pred_label = np.array([1 if score >= 0 else -1 for score in pred_score])
true_label = np.array([B.y for B in bags_test])
a = accuracy (pred_label, true_label) # accuracy
p = precision(pred_label, true_label) # precision
r = recall (pred_label, true_label) # recall
f = f_score (pred_label, true_label) # F-score
auc = metrics.roc_auc_score((true_label+1)/2, pred_score)
if not args.aucplot:
sys.stdout.write("""# accuracy,precision,recall,f-score,ROC-AUC
{:.3f},{:.3f},{:.3f},{:.3f},{:.3f}\n""".format(a, p, r, f, auc))
sys.stdout.flush()
else:
sys.stdout.write("""# accuracy,precision,recall,f-score,ROC-AUC
# {:.3f},{:.3f},{:.3f},{:.3f},{:.3f}\n""".format(a, p, r, f, auc))
sys.stdout.flush()
np.savetxt(sys.stdout.buffer, np.c_[pred_score, true_label])
def trotx(theta, unit="rad", xyz=[0, 0, 0]):
"""
TROTX Rotation about X axis
:param theta: rotation in radians or degrees
:param unit: "rad" or "deg" to indicate unit being used
:param xyz: the xyz translation, if blank defaults to [0,0,0]
:return: homogeneous transform matrix
trotx(THETA) is a homogeneous transformation (4x4) representing a rotation
of THETA radians about the x-axis.
trotx(THETA, 'deg') as above but THETA is in degrees
trotx(THETA, 'rad', [x,y,z]) as above with translation of [x,y,z]
"""
check_args.unit_check(unit)
tm = rotx(theta, unit)
tm = np.r_[tm, np.zeros((1, 3))]
mat = np.c_[tm, np.array([[xyz[0]], [xyz[1]], [xyz[2]], [1]])]
mat = np.asmatrix(mat.round(15))
return mat
# ---------------------------------------------------------------------------------------#
def troty(theta, unit="rad", xyz=[0, 0, 0]):
"""
TROTY Rotation about Y axis
:param theta: rotation in radians or degrees
:param unit: "rad" or "deg" to indicate unit being used
:param xyz: the xyz translation, if blank defaults to [0,0,0]
:return: homogeneous transform matrix
troty(THETA) is a homogeneous transformation (4x4) representing a rotation
of THETA radians about the y-axis.
troty(THETA, 'deg') as above but THETA is in degrees
troty(THETA, 'rad', [x,y,z]) as above with translation of [x,y,z]
"""
check_args.unit_check(unit)
tm = roty(theta, unit)
tm = np.r_[tm, np.zeros((1, 3))]
mat = np.c_[tm, np.array([[xyz[0]], [xyz[1]], [xyz[2]], [1]])]
mat = np.asmatrix(mat.round(15))
return mat
# ---------------------------------------------------------------------------------------#
def trotz(theta, unit="rad", xyz=[0, 0, 0]):
"""
TROTZ Rotation about Z axis
:param theta: rotation in radians or degrees
:param unit: "rad" or "deg" to indicate unit being used
:param xyz: the xyz translation, if blank defaults to [0,0,0]
:return: homogeneous transform matrix
trotz(THETA) is a homogeneous transformation (4x4) representing a rotation
of THETA radians about the z-axis.
trotz(THETA, 'deg') as above but THETA is in degrees
trotz(THETA, 'rad', [x,y,z]) as above with translation of [x,y,z]
"""
check_args.unit_check(unit)
tm = rotz(theta, unit)
tm = np.r_[tm, np.zeros((1, 3))]
mat = np.c_[tm, np.array([[xyz[0]], [xyz[1]], [xyz[2]], [1]])]
mat = np.asmatrix(mat.round(15))
return mat
# ---------------------------------------------------------------------------------------#
def trot2(theta, unit='rad'):
"""
TROT2 SE2 rotation matrix
:param theta: rotation in radians or degrees
:param unit: "rad" or "deg" to indicate unit being used
:return: homogeneous transform matrix (3x3)
TROT2(THETA) is a homogeneous transformation (3x3) representing a rotation of
THETA radians.
TROT2(THETA, 'deg') as above but THETA is in degrees.
Notes::
- Translational component is zero.
"""
tm = rot2(theta, unit)
tm = np.r_[tm, np.zeros((1, 2))]
mat = np.c_[tm, np.array([[0], [0], [1]])]
return mat
# ---------------------------------------------------------------------------------------#
def _scipy_bivariate_kde(x, y, bw, gridsize, cut, clip):
"""Compute a bivariate kde using scipy."""
data = np.c_[x, y]
kde = stats.gaussian_kde(data.T)
data_std = data.std(axis=0, ddof=1)
if isinstance(bw, string_types):
bw = "scotts" if bw == "scott" else bw
bw_x = getattr(kde, "%s_factor" % bw)() * data_std[0]
bw_y = getattr(kde, "%s_factor" % bw)() * data_std[1]
elif np.isscalar(bw):
bw_x, bw_y = bw, bw
else:
msg = ("Cannot specify a different bandwidth for each dimension "
"with the scipy backend. You should install statsmodels.")
raise ValueError(msg)
x_support = _kde_support(data[:, 0], bw_x, gridsize, cut, clip[0])
y_support = _kde_support(data[:, 1], bw_y, gridsize, cut, clip[1])
xx, yy = np.meshgrid(x_support, y_support)
z = kde([xx.ravel(), yy.ravel()]).reshape(xx.shape)
return xx, yy, z
def spiral(num_cls, dim, point_per_cls, rnd_state=1024):
np.random.seed(rnd_state)
points_per_cls = 100 # number of points per class
dim = 2 # dimensionality
num_cls = 3 # number of classes
X_data = np.zeros((points_per_cls * num_cls, dim))
y_data = np.zeros(points_per_cls * num_cls, dtype='uint8')
for j in range(num_cls):
ix = range(points_per_cls * j, points_per_cls * (j + 1))
r = np.linspace(0.0, 1, points_per_cls)
t = np.linspace(j * 4, (j + 1) * 4, points_per_cls) + np.random.randn(points_per_cls) * 0.2 # theta
X_data[ix] = np.c_[r * np.sin(t), r * np.cos(t)]
y_data[ix] = j
y_data_encoded = np.zeros((points_per_cls * num_cls, num_cls))
y_data_encoded[range(points_per_cls * num_cls), y_data] = 1
return X_data, y_data, y_data_encoded
def gradient(x0, X, y, alpha):
# gradient of the logistic loss
w, c = x0[1:137], x0[0]
#print("c is " + str(c))
z = X.dot(w) + c
z = phi(y * z)
z0 = (z - 1) * y
grad_w = np.matmul(z0,X) / X.shape[0] + alpha * w
grad_c = z0.sum() / X.shape[0]
grad_c = np.array(grad_c)
#print(grad_w[0,1:5])
return np.c_[([grad_c], grad_w)]
##### Stochastic Gradient Descent Optimiser ######
def average_ndcg(labels, query_ids, predicted_labels):
ndcg_list = np.zeros(len(set(query_ids)))
k = 0
for i in set(query_ids):
idx = [query_ids == i]
orders = np.c_[labels[idx],predicted_labels[idx]]
sorted_orders = orders[orders[:,1].argsort()[::-1]][:,0]
ndcg_list[k] = ndcg(sorted_orders)
k +=1
if k%2000 == 0:
print(str(k) + " queries calculated")
print("mean ndcg so far: " + str(np.mean(ndcg_list[0:k])))
return np.mean(ndcg_list)
# average ndcg is 0.26333
def average_ndcg(labels, query_ids, predicted_labels):
ndcg_list = np.zeros(len(set(query_ids)))
k = 0
for i in set(query_ids):
idx = [query_ids == i]
orders = np.c_[labels[idx],predicted_labels[idx]]
sorted_orders = orders[orders[:,1].argsort()[::-1]][:,0]
ndcg_list[k] = ndcg(sorted_orders)
k +=1
if k%2000 == 0:
print(str(k) + " queries calculated")
print("mean ndcg so far: " + str(np.mean(ndcg_list[0:k])))
return np.mean(ndcg_list)
# average ndcg is 0.26333
def computeGaussianWidthCandidates(self, referenceSamples=None, testSamples=None) :
"""
Compute a candidate list of Gaussian kernel widths. The best width will be
selected via cross-validation
"""
allSamples = numpy.c_[referenceSamples, testSamples]
medianDistance = self.getMedianDistanceBetweenSamples(allSamples.T)
return medianDistance * numpy.array([0.6, 0.8, 1, 1.2, 1.4])
def _fit(self,X,y=None):
if isinstance(X, pyisc.DataObject) and y is None:
assert y is None # Contained in the data object
self.class_column = X.class_column
if self.class_column >= 0:
self.classes_ = X.classes_
self._anomaly_detector._SetParams(
0,
-1 if X.class_column is None else X.class_column,
self.anomaly_threshold,
1 if self.is_clustering else 0
)
self._anomaly_detector._TrainData(X)
return self
if isinstance(X, ndarray):
class_column = -1
data_object = None
assert X.ndim <= 2
if X.ndim == 2:
max_class_column = X.shape[1]
else:
max_class_column = 1
if isinstance(y, list) or isinstance(y, ndarray):
assert len(X) == len(y)
class_column = max_class_column
data_object = pyisc.DataObject(numpy.c_[X, y], class_column=class_column)
elif y is None or int(y) == y and y > -1 and y <= max_class_column:
self.class_column = y
data_object = pyisc.DataObject(X,class_column=y)
if data_object is not None:
return self._fit(data_object)
raise ValueError("Unknown type of data to fit X, y:", type(X), type(y))
def _convert_to_data_object_in_scoring(self, X, y):
data_object = None
if isinstance(y, list) or isinstance(y, ndarray):
assert X.ndim == 2 and self.class_column == X.shape[1] or X.ndim == 1 and self.class_column == 1
data_object = pyisc.DataObject(numpy.c_[X, y], class_column=self.class_column,classes=self.classes_)
else:
assert self.class_column == y
data_object = pyisc.DataObject(X, class_column=self.class_column,classes=self.classes_ if y is not None else None)
return data_object
def loglikelihood(self,X,y=None):
assert isinstance(X, ndarray) and (self.class_column is None and y is None or len(y) == len(X))
if y is not None:
return self._anomaly_detector._LogProbabilityOfData(pyisc.DataObject(c_[X,y], class_column=len(X[0])), len(X)).sum()
else:
return self._anomaly_detector._LogProbabilityOfData(pyisc.DataObject(X), len(X)).sum()
def test_dataobject_set_column_values(self):
X = array([norm(1.0).rvs(10) for _ in range(1000)])
y = [None] * 1000
DO = DataObject(c_[X,y], class_column=len(X[0]))
assert_equal(len(X[0]), DO.class_column)
assert_equal(unique(y), DO.classes_)
classes=[None] + ['1', '2', '3', '4', '5']
DO = DataObject(c_[X,y], class_column=len(X[0]), classes=classes)
assert_equal(len(X[0]), DO.class_column)
assert_equal(classes, DO.classes_)
X2 = DO.as_2d_array()
assert_allclose(X2.T[:-1].T.astype(float), X)
assert_equal(X2.T[-1],y)
new_y = ["%i"%(divmod(i,5)[1]+1) for i in range(len(X))]
DO.set_column_values(len(X[0]), new_y)
assert_equal(len(X[0]), DO.class_column)
assert_equal([None]+list(unique(new_y)), DO.classes_)
X2 = DO.as_2d_array()
assert_allclose(X2.T[:-1].T.astype(float), X)
assert_equal(X2.T[-1], new_y)
def test_outlier_detection(self):
print "Start of test"
n_samples = 1000
norm_dist = stats.norm(0, 1)
truth = np.ones((n_samples,))
truth[-100:] = -1
X0 = norm_dist.rvs(n_samples)
X = np.c_[X0*5, X0+norm_dist.rvs(n_samples)*2]
uniform_dist = stats.uniform(-10,10)
X[-100:] = np.c_[uniform_dist.rvs(100),uniform_dist.rvs(100)]
outlier_detector = pyisc.SklearnOutlierDetector(
100.0/n_samples,
pyisc.P_Gaussian([0,1])
)
outlier_detector.fit(X, np.array([1]*len(X)))
self.assertLess(outlier_detector.threshold_, 0.35)
self.assertGreater(outlier_detector.threshold_, 0.25)
predictions = outlier_detector.predict(X, np.array([1]*len(X)))
accuracy = sum(truth == predictions)/float(n_samples)
print "accuracy", accuracy
self.assertGreater(accuracy, 0.85)