def test_mode_all_but_economic(self):
a = array([[1, 2], [3, 4]])
b = array([[1, 2], [3, 4], [5, 6]])
for dt in "fd":
m1 = a.astype(dt)
m2 = b.astype(dt)
self.check_qr(m1)
self.check_qr(m2)
self.check_qr(m2.T)
self.check_qr(matrix(m1))
for dt in "fd":
m1 = 1 + 1j * a.astype(dt)
m2 = 1 + 1j * b.astype(dt)
self.check_qr(m1)
self.check_qr(m2)
self.check_qr(m2.T)
self.check_qr(matrix(m1))
python类matrix()的实例源码
def test_basic(self):
A = np.array([[1, 2], [3, 4]])
mA = matrix(A)
assert_(np.all(mA.A == A))
B = bmat("A,A;A,A")
C = bmat([[A, A], [A, A]])
D = np.array([[1, 2, 1, 2],
[3, 4, 3, 4],
[1, 2, 1, 2],
[3, 4, 3, 4]])
assert_(np.all(B.A == D))
assert_(np.all(C.A == D))
E = np.array([[5, 6], [7, 8]])
AEresult = matrix([[1, 2, 5, 6], [3, 4, 7, 8]])
assert_(np.all(bmat([A, E]) == AEresult))
vec = np.arange(5)
mvec = matrix(vec)
assert_(mvec.shape == (1, 5))
def test_sum(self):
"""Test whether matrix.sum(axis=1) preserves orientation.
Fails in NumPy <= 0.9.6.2127.
"""
M = matrix([[1, 2, 0, 0],
[3, 4, 0, 0],
[1, 2, 1, 2],
[3, 4, 3, 4]])
sum0 = matrix([8, 12, 4, 6])
sum1 = matrix([3, 7, 6, 14]).T
sumall = 30
assert_array_equal(sum0, M.sum(axis=0))
assert_array_equal(sum1, M.sum(axis=1))
assert_equal(sumall, M.sum())
assert_array_equal(sum0, np.sum(M, axis=0))
assert_array_equal(sum1, np.sum(M, axis=1))
assert_equal(sumall, np.sum(M))
def test_basic(self):
import numpy.linalg as linalg
A = np.array([[1., 2.],
[3., 4.]])
mA = matrix(A)
assert_(np.allclose(linalg.inv(A), mA.I))
assert_(np.all(np.array(np.transpose(A) == mA.T)))
assert_(np.all(np.array(np.transpose(A) == mA.H)))
assert_(np.all(A == mA.A))
B = A + 2j*A
mB = matrix(B)
assert_(np.allclose(linalg.inv(B), mB.I))
assert_(np.all(np.array(np.transpose(B) == mB.T)))
assert_(np.all(np.array(np.transpose(B).conj() == mB.H)))
def skew(v):
return np.matrix([[0,-v[2,0],v[1,0]], [v[2,0],0,-v[0,0]], [-v[1,0],v[0,0],0]])
def vector3(x, y, z):
return np.matrix([[x],[y],[z]])
def vector6(a, b, c, x, y, z):
return np.matrix([[a],[b],[c],[x],[y],[z]])
def col(v):
col = [[x] for x in v]
return np.matrix(col)
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 analytic_convolution_gaussian(mu1,covar1,mu2,covar2):
"""
The analytic vconvolution of two Gaussians is simply the sum of the two mean vectors
and the two convariance matrixes
--- INPUT ---
mu1 The mean of the first gaussian
covar1 The covariance matrix of of the first gaussian
mu2 The mean of the second gaussian
covar2 The covariance matrix of of the second gaussian
"""
muconv = mu1+mu2
covarconv = covar1+covar2
return muconv, covarconv
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
def get_facial_landmarks_from_mask(img, pts):
rect = cv2.boundingRect(pts)
rect = dlib.rectangle(rect[0], rect[1], rect[0] + rect[2], rect[1] + rect[3])
return np.matrix([list(pt) for pt in pts]), rect
def get_facial_landmarks(img):
# No need to upsample
rects = face_detector(img, 0)
if len(rects) == 0:
print "No faces"
return None
rect = rects[0]
shape = shape_predictor(img, rect)
return np.matrix([[pt.x, pt.y] for pt in shape.parts()]), rect
def get_tm_opp(pts1, pts2):
# Transformation matrix - ( Translation + Scaling + Rotation )
# using Procuster analysis
pts1 = np.float64(pts1)
pts2 = np.float64(pts2)
m1 = np.mean(pts1, axis = 0)
m2 = np.mean(pts2, axis = 0)
# Removing translation
pts1 -= m1
pts2 -= m2
std1 = np.std(pts1)
std2 = np.std(pts2)
std_r = std2/std1
# Removing scaling
pts1 /= std1
pts2 /= std2
U, S, V = np.linalg.svd(np.transpose(pts1) * pts2)
# Finding the rotation matrix
R = np.transpose(U * V)
return np.vstack([np.hstack((std_r * R,
np.transpose(m2) - std_r * R * np.transpose(m1))), np.matrix([0.0, 0.0, 1.0])])
def as_float_array(X, copy=True, force_all_finite=True):
"""Converts an array-like to an array of floats
The new dtype will be np.float32 or np.float64, depending on the original
type. The function can create a copy or modify the argument depending
on the argument copy.
Parameters
----------
X : {array-like, sparse matrix}
copy : bool, optional
If True, a copy of X will be created. If False, a copy may still be
returned if X's dtype is not a floating point type.
force_all_finite : boolean (default=True)
Whether to raise an error on np.inf and np.nan in X.
Returns
-------
XT : {array, sparse matrix}
An array of type np.float
"""
if isinstance(X, np.matrix) or (not isinstance(X, np.ndarray)
and not sp.issparse(X)):
return check_array(X, ['csr', 'csc', 'coo'], dtype=np.float64,
copy=copy, force_all_finite=force_all_finite,
ensure_2d=False)
elif sp.issparse(X) and X.dtype in [np.float32, np.float64]:
return X.copy() if copy else X
elif X.dtype in [np.float32, np.float64]: # is numpy array
return X.copy('F' if X.flags['F_CONTIGUOUS'] else 'C') if copy else X
else:
return X.astype(np.float32 if X.dtype == np.int32 else np.float64)
def get_precision(self):
"""Compute data precision matrix with the generative model.
Equals the inverse of the covariance but computed with
the matrix inversion lemma for efficiency.
Returns
-------
precision : array, shape=(n_features, n_features)
Estimated precision of data.
"""
n_features = self.components_.shape[1]
# handle corner cases first
if self.n_components_ == 0:
return np.eye(n_features) / self.noise_variance_
if self.n_components_ == n_features:
return linalg.inv(self.get_covariance())
# Get precision using matrix inversion lemma
components_ = self.components_
exp_var = self.explained_variance_
exp_var_diff = np.maximum(exp_var - self.noise_variance_, 0.)
precision = np.dot(components_, components_.T) / self.noise_variance_
precision.flat[::len(precision) + 1] += 1. / exp_var_diff
precision = np.dot(components_.T,
np.dot(linalg.inv(precision), components_))
precision /= -(self.noise_variance_ ** 2)
precision.flat[::len(precision) + 1] += 1. / self.noise_variance_
return precision
def diffD1D2(Flags):
# To check if the difference between D1 and D2 amplifies downstream
# First decide which model to use
delay = 1
if Flags == "Allsym":
(d1,d2,fsi,ti,ta,stn,gpi,ipctx,A,B,params) = calcRates(Flags,delay)
# D = Direct pathway, ID = Indirect pathway, HD = Hyperdirect pathway
# Reducing a full recurrent matrix leads to postive and negative contributions in ID and HD instead of pure just positive contributions
D = params['gpid1']
print "Direct",D
de1 = 1. + (params['d1d1']*params['fsifsi']) - (params['stnti']*params['stnstn'])
Ex1 = (params['stnti']*params['d1d1']*params['fsifsi']*params['gpistn'])/de1
Ex2 = (params['gpid1']*params['stnstn']*params['d1ta']*params['fsifsi'])/de1
IDpos = params['stnti']*params['tid2']*(params['d1ta']*params['gpid1']*params['tistn']+ Ex1 + Ex2)
print "IDpos",IDpos
Ex3 = (params['d1ta']+params['d1ta']*params['tata']+((params['stnta']*params['fsiti']*params['d1fsi'])/de1))
IDneg = params['gpid1']+params['gpiti']*params['tid2']+params['gpid1']*params['stnstn']*params['tid2']*Ex3
print "IDneg",IDneg
HDpos = (params['jstnctx']*params['gpid1']*params['stnstn']*params['fsiti']*params['d1fsi'])/de1
print "HDpos",HDpos
Ex4 = params['jstnctx']*params['d1d1']*params['fsifsi']*params['stnti']*params['gpistn']
Ex5 = params['jstnctx']*params['gpid1']*params['stnstn']*params['d1ta']*params['fsifsi']
HDneg = (Ex4 + Ex5)/de1
print "HDneg",HDneg
d1fix = np.mean(d1[:-10])
d2fix = np.mean(d2[:-10])
DelMSN = d1fix - d2fix
DelGpi = (D*d1fix) + ((IDpos - IDneg)*d2fix)
return (DelMSN,DelGpi)
def diffD1D2(Flags):
# To check if the difference between D1 and D2 amplifies downstream
# First decide which model to use
delay = 1
if Flags == "Allsym":
(d1,d2,fsi,ti,ta,stn,gpi,ipctx,A,B,params) = calcRates(Flags,delay)
# D = Direct pathway, ID = Indirect pathway, HD = Hyperdirect pathway
# Reducing a full recurrent matrix leads to postive and negative contributions in ID and HD instead of pure just positive contributions
D = params['gpid1']
print "Direct",D
de1 = 1. + (params['d1d1']*params['fsifsi']) - (params['stnti']*params['stnstn'])
Ex1 = (params['stnti']*params['d1d1']*params['fsifsi']*params['gpistn'])/de1
Ex2 = (params['gpid1']*params['stnstn']*params['d1ta']*params['fsifsi'])/de1
IDpos = params['stnti']*params['tid2']*(params['d1ta']*params['gpid1']*params['tistn']+ Ex1 + Ex2)
print "IDpos",IDpos
Ex3 = (params['d1ta']+params['d1ta']*params['tata']+((params['stnta']*params['fsiti']*params['d1fsi'])/de1))
IDneg = params['gpid1']+params['gpiti']*params['tid2']+params['gpid1']*params['stnstn']*params['tid2']*Ex3
print "IDneg",IDneg
HDpos = (params['jstnctx']*params['gpid1']*params['stnstn']*params['fsiti']*params['d1fsi'])/de1
print "HDpos",HDpos
Ex4 = params['jstnctx']*params['d1d1']*params['fsifsi']*params['stnti']*params['gpistn']
Ex5 = params['jstnctx']*params['gpid1']*params['stnstn']*params['d1ta']*params['fsifsi']
HDneg = (Ex4 + Ex5)/de1
print "HDneg",HDneg
d1fix = np.mean(d1[:-10])
d2fix = np.mean(d2[:-10])
DelMSN = d1fix - d2fix
DelGpi = (D*d1fix) + ((IDpos - IDneg)*d2fix)
return (DelMSN,DelGpi)
def diffD1D2(Flags):
# To check if the difference between D1 and D2 amplifies downstream
# First decide which model to use
delay = 1
if Flags == "Allsym":
(d1,d2,fsi,ti,ta,stn,gpi,ipctx,A,B,params) = calcRates(Flags,delay)
# D = Direct pathway, ID = Indirect pathway, HD = Hyperdirect pathway
# Reducing a full recurrent matrix leads to postive and negative contributions in ID and HD instead of pure just positive contributions
D = params['gpid1']
print "Direct",D
de1 = 1. + (params['d1d1']*params['fsifsi']) - (params['stnti']*params['stnstn'])
Ex1 = (params['stnti']*params['d1d1']*params['fsifsi']*params['gpistn'])/de1
Ex2 = (params['gpid1']*params['stnstn']*params['d1ta']*params['fsifsi'])/de1
IDpos = params['stnti']*params['tid2']*(params['d1ta']*params['gpid1']*params['tistn']+ Ex1 + Ex2)
print "IDpos",IDpos
Ex3 = (params['d1ta']+params['d1ta']*params['tata']+((params['stnta']*params['fsiti']*params['d1fsi'])/de1))
IDneg = params['gpid1']+params['gpiti']*params['tid2']+params['gpid1']*params['stnstn']*params['tid2']*Ex3
print "IDneg",IDneg
HDpos = (params['jstnctx']*params['gpid1']*params['stnstn']*params['fsiti']*params['d1fsi'])/de1
print "HDpos",HDpos
Ex4 = params['jstnctx']*params['d1d1']*params['fsifsi']*params['stnti']*params['gpistn']
Ex5 = params['jstnctx']*params['gpid1']*params['stnstn']*params['d1ta']*params['fsifsi']
HDneg = (Ex4 + Ex5)/de1
print "HDneg",HDneg
d1fix = np.mean(d1[:-10])
d2fix = np.mean(d2[:-10])
DelMSN = d1fix - d2fix
DelGpi = (D*d1fix) + ((IDpos - IDneg)*d2fix)
return (DelMSN,DelGpi)
def diffD1D2(Flags):
# To check if the difference between D1 and D2 amplifies downstream
# First decide which model to use
delay = 1
if Flags == "Allsym":
(d1,d2,fsi,ti,ta,stn,gpi,ipctx,A,B,params) = calcRates(Flags,delay)
# D = Direct pathway, ID = Indirect pathway, HD = Hyperdirect pathway
# Reducing a full recurrent matrix leads to postive and negative contributions in ID and HD instead of pure just positive contributions
D = params['gpid1']
print "Direct",D
de1 = 1. + (params['d1d1']*params['fsifsi']) - (params['stnti']*params['stnstn'])
Ex1 = (params['stnti']*params['d1d1']*params['fsifsi']*params['gpistn'])/de1
Ex2 = (params['gpid1']*params['stnstn']*params['d1ta']*params['fsifsi'])/de1
IDpos = params['stnti']*params['tid2']*(params['d1ta']*params['gpid1']*params['tistn']+ Ex1 + Ex2)
print "IDpos",IDpos
Ex3 = (params['d1ta']+params['d1ta']*params['tata']+((params['stnta']*params['fsiti']*params['d1fsi'])/de1))
IDneg = params['gpid1']+params['gpiti']*params['tid2']+params['gpid1']*params['stnstn']*params['tid2']*Ex3
print "IDneg",IDneg
HDpos = (params['jstnctx']*params['gpid1']*params['stnstn']*params['fsiti']*params['d1fsi'])/de1
print "HDpos",HDpos
Ex4 = params['jstnctx']*params['d1d1']*params['fsifsi']*params['stnti']*params['gpistn']
Ex5 = params['jstnctx']*params['gpid1']*params['stnstn']*params['d1ta']*params['fsifsi']
HDneg = (Ex4 + Ex5)/de1
print "HDneg",HDneg
d1fix = np.mean(d1[:-10])
d2fix = np.mean(d2[:-10])
DelMSN = d1fix - d2fix
DelGpi = (D*d1fix) + ((IDpos - IDneg)*d2fix)
return (DelMSN,DelGpi)
