def treegauss_remove_row(
data_row,
tree_grid,
latent_row,
vert_ss,
edge_ss,
feat_ss, ):
# Update sufficient statistics.
for v in range(latent_row.shape[0]):
z = latent_row[v, :]
vert_ss[v, :, :] -= np.outer(z, z)
for e in range(tree_grid.shape[1]):
z1 = latent_row[tree_grid[1, e], :]
z2 = latent_row[tree_grid[2, e], :]
edge_ss[e, :, :] -= np.outer(z1, z2)
for v, x in enumerate(data_row):
if np.isnan(x):
continue
z = latent_row[v, :]
feat_ss[v] -= 1
feat_ss[v, 1] -= x
feat_ss[v, 2:] -= x * z # TODO Use central covariance.
python类outer()的实例源码
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
transformations.py 文件源码
项目：Neural-Networks-for-Inverse-Kinematics
作者: paramrajpura
项目源码
文件源码
阅读 38
收藏 0
点赞 0
评论 0
def reflection_matrix(point, normal):
"""Return matrix to mirror at plane defined by point and normal vector.
>>> v0 = numpy.random.random(4) - 0.5
>>> v0[3] = 1.
>>> v1 = numpy.random.random(3) - 0.5
>>> R = reflection_matrix(v0, v1)
>>> numpy.allclose(2, numpy.trace(R))
True
>>> numpy.allclose(v0, numpy.dot(R, v0))
True
>>> v2 = v0.copy()
>>> v2[:3] += v1
>>> v3 = v0.copy()
>>> v2[:3] -= v1
>>> numpy.allclose(v2, numpy.dot(R, v3))
True
"""
normal = unit_vector(normal[:3])
M = numpy.identity(4)
M[:3, :3] -= 2.0 * numpy.outer(normal, normal)
M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal
return M
def unscentedTransform(X, Wm, Wc, f):
Y = None
Ymean = None
fdim = None
N = np.shape(X)[1]
for j in range(0,N):
fImage = f(X[:,j])
if Y is None:
fdim = np.size(fImage)
Y = np.zeros((fdim, np.shape(X)[1]))
Ymean = np.zeros(fdim)
Y[:,j] = fImage
Ymean += Wm[j] * Y[:,j]
Ycov = np.zeros((fdim, fdim))
for j in range(0, N):
meanAdjustedYj = Y[:,j] - Ymean
Ycov += np.outer(Wc[j] * meanAdjustedYj, meanAdjustedYj)
return Y, Ymean, Ycov
def reflection_matrix(point, normal):
"""Return matrix to mirror at plane defined by point and normal vector.
>>> v0 = numpy.random.random(4) - 0.5
>>> v0[3] = 1.0
>>> v1 = numpy.random.random(3) - 0.5
>>> R = reflection_matrix(v0, v1)
>>> numpy.allclose(2., numpy.trace(R))
True
>>> numpy.allclose(v0, numpy.dot(R, v0))
True
>>> v2 = v0.copy()
>>> v2[:3] += v1
>>> v3 = v0.copy()
>>> v2[:3] -= v1
>>> numpy.allclose(v2, numpy.dot(R, v3))
True
"""
normal = unit_vector(normal[:3])
M = numpy.identity(4)
M[:3, :3] -= 2.0 * numpy.outer(normal, normal)
M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal
return M
def estimate_params(self, ess):
"""
Estimate Nomal param,
given expecteed sufficient stats
"""
n = len(ess)
mu = np.asarray([0,0])
for (m, c) in ess:
mu = mu + m
mu = mu * 1.0 / n
C = np.asarray([[0,0],[0,0]])
for (m, c) in ess:
C = C + c
C = C * 1.0 / n
C = C - np.outer(mu, mu)
if not self.full_cov:
C = reduce_cov(C)
return (mu, C)
def compute_pvalues_for_processes(self,U_matrix,chane_prob, num_bootstrapped_stats=100):
N = U_matrix.shape[0]
bootsraped_stats = np.zeros(num_bootstrapped_stats)
# orsetinW = simulate(N,num_bootstrapped_stats,corr)
for proc in range(num_bootstrapped_stats):
# W = np.sign(orsetinW[:,proc])
W = simulatepm(N,chane_prob)
WW = np.outer(W, W)
st = np.mean(U_matrix * WW)
bootsraped_stats[proc] = N * st
stat = N*np.mean(U_matrix)
return float(np.sum(bootsraped_stats > stat)) / num_bootstrapped_stats
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 _eig_local_op_mps(lv, ltens, rv):
"""Local operator contribution from an MPS"""
# MPS 1 / ltens: Interpreted as |psiXpsi| part of the operator
# MPS 2: The current eigvectector candidate
op = lv.T
# op axes: 0 mps2 bond, 1: mps1 bond
s = op.shape
op = op.reshape((s[0], 1, s[1]))
# op axes: 0 mps2 bond, 1: physical legs, 2: mps1 bond
for lt in ltens:
# op axes: 0: mps2 bond, 1: physical legs, 2: mps1 bond
op = np.tensordot(op, lt.conj(), axes=(2, 0))
# op axes: 0: mps2 bond, 1, 2: physical legs, 3: mps1 bond
s = op.shape
op = op.reshape((s[0], -1, s[3]))
# op axes: 0: mps2 bond, 1: physical legs, 2: mps1 bond
op = np.tensordot(op, rv, axes=(2, 0))
# op axes: 0: mps2 bond, 1: physical legs, 2: mps2 bond
op = np.outer(op.conj(), op)
# op axes:
# 0: (0a: left cc mps2 bond, 0b: physical row leg, 0c: right cc mps2 bond),
# 1: (1a: left mps2 bond, 1b: physical column leg, 1c: right mps2 bond)
return op
def test_mps_to_mpo(nr_sites, local_dim, rank, rgen):
mps = factory.random_mps(nr_sites, local_dim, rank, randstate=rgen)
# Instead of calling the two functions, we call mps_to_mpo(),
# which does exactly that:
# mps_as_puri = mp.mps_as_local_purification_mps(mps)
# mpo = mp.pmps_to_mpo(mps_as_puri)
mpo = mm.mps_to_mpo(mps)
# This is also a test of mp.mps_as_local_purification_mps() in the
# following sense: Local purifications are representations of
# mixed states. Therefore, compare mps and mps_as_puri by
# converting them to mixed states.
state = mps.to_array()
state = np.outer(state, state.conj())
state.shape = (local_dim,) * (2 * nr_sites)
state2 = mpo.to_array_global()
assert_array_almost_equal(state, state2)
def reflection_matrix(point, normal):
"""Return matrix to mirror at plane defined by point and normal vector.
>>> v0 = numpy.random.random(4) - 0.5
>>> v0[3] = 1.
>>> v1 = numpy.random.random(3) - 0.5
>>> R = reflection_matrix(v0, v1)
>>> numpy.allclose(2, numpy.trace(R))
True
>>> numpy.allclose(v0, numpy.dot(R, v0))
True
>>> v2 = v0.copy()
>>> v2[:3] += v1
>>> v3 = v0.copy()
>>> v2[:3] -= v1
>>> numpy.allclose(v2, numpy.dot(R, v3))
True
"""
normal = unit_vector(normal[:3])
M = numpy.identity(4)
M[:3, :3] -= 2.0 * numpy.outer(normal, normal)
M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal
return M
def vol(self):
"""Construct cell volumes of the 3D model as 1d array."""
if getattr(self, '_vol', None) is None:
vh = self.h
# Compute cell volumes
if self.dim == 1:
self._vol = utils.mkvc(vh[0])
elif self.dim == 2:
# Cell sizes in each direction
self._vol = utils.mkvc(np.outer(vh[0], vh[1]))
elif self.dim == 3:
# Cell sizes in each direction
self._vol = utils.mkvc(
np.outer(utils.mkvc(np.outer(vh[0], vh[1])), vh[2])
)
return self._vol
def test_minimummaximum_func(self):
a = np.ones((2, 2))
aminimum = minimum(a, a)
self.assertTrue(isinstance(aminimum, MaskedArray))
assert_equal(aminimum, np.minimum(a, a))
aminimum = minimum.outer(a, a)
self.assertTrue(isinstance(aminimum, MaskedArray))
assert_equal(aminimum, np.minimum.outer(a, a))
amaximum = maximum(a, a)
self.assertTrue(isinstance(amaximum, MaskedArray))
assert_equal(amaximum, np.maximum(a, a))
amaximum = maximum.outer(a, a)
self.assertTrue(isinstance(amaximum, MaskedArray))
assert_equal(amaximum, np.maximum.outer(a, a))
def test_TakeTransposeInnerOuter(self):
# Test of take, transpose, inner, outer products
x = arange(24)
y = np.arange(24)
x[5:6] = masked
x = x.reshape(2, 3, 4)
y = y.reshape(2, 3, 4)
assert_equal(np.transpose(y, (2, 0, 1)), transpose(x, (2, 0, 1)))
assert_equal(np.take(y, (2, 0, 1), 1), take(x, (2, 0, 1), 1))
assert_equal(np.inner(filled(x, 0), filled(y, 0)),
inner(x, y))
assert_equal(np.outer(filled(x, 0), filled(y, 0)),
outer(x, y))
y = array(['abc', 1, 'def', 2, 3], object)
y[2] = masked
t = take(y, [0, 3, 4])
assert_(t[0] == 'abc')
assert_(t[1] == 2)
assert_(t[2] == 3)
def test_testTakeTransposeInnerOuter(self):
# Test of take, transpose, inner, outer products
x = arange(24)
y = np.arange(24)
x[5:6] = masked
x = x.reshape(2, 3, 4)
y = y.reshape(2, 3, 4)
assert_(eq(np.transpose(y, (2, 0, 1)), transpose(x, (2, 0, 1))))
assert_(eq(np.take(y, (2, 0, 1), 1), take(x, (2, 0, 1), 1)))
assert_(eq(np.inner(filled(x, 0), filled(y, 0)),
inner(x, y)))
assert_(eq(np.outer(filled(x, 0), filled(y, 0)),
outer(x, y)))
y = array(['abc', 1, 'def', 2, 3], object)
y[2] = masked
t = take(y, [0, 3, 4])
assert_(t[0] == 'abc')
assert_(t[1] == 2)
assert_(t[2] == 3)
def test_4(self):
"""
Test of take, transpose, inner, outer products.
"""
x = self.arange(24)
y = np.arange(24)
x[5:6] = self.masked
x = x.reshape(2, 3, 4)
y = y.reshape(2, 3, 4)
assert self.allequal(np.transpose(y, (2, 0, 1)), self.transpose(x, (2, 0, 1)))
assert self.allequal(np.take(y, (2, 0, 1), 1), self.take(x, (2, 0, 1), 1))
assert self.allequal(np.inner(self.filled(x, 0), self.filled(y, 0)),
self.inner(x, y))
assert self.allequal(np.outer(self.filled(x, 0), self.filled(y, 0)),
self.outer(x, y))
y = self.array(['abc', 1, 'def', 2, 3], object)
y[2] = self.masked
t = self.take(y, [0, 3, 4])
assert t[0] == 'abc'
assert t[1] == 2
assert t[2] == 3
def outer(self, a, b):
"""
Return the function applied to the outer product of a and b.
"""
(da, db) = (getdata(a), getdata(b))
d = self.f.outer(da, db)
ma = getmask(a)
mb = getmask(b)
if ma is nomask and mb is nomask:
m = nomask
else:
ma = getmaskarray(a)
mb = getmaskarray(b)
m = umath.logical_or.outer(ma, mb)
if (not m.ndim) and m:
return masked
if m is not nomask:
np.copyto(d, da, where=m)
if not d.shape:
return d
masked_d = d.view(get_masked_subclass(a, b))
masked_d._mask = m
return masked_d
def reflection_matrix(point, normal):
"""Return matrix to mirror at plane defined by point and normal vector.
>>> v0 = numpy.random.random(4) - 0.5
>>> v0[3] = 1.
>>> v1 = numpy.random.random(3) - 0.5
>>> R = reflection_matrix(v0, v1)
>>> numpy.allclose(2, numpy.trace(R))
True
>>> numpy.allclose(v0, numpy.dot(R, v0))
True
>>> v2 = v0.copy()
>>> v2[:3] += v1
>>> v3 = v0.copy()
>>> v2[:3] -= v1
>>> numpy.allclose(v2, numpy.dot(R, v3))
True
"""
normal = unit_vector(normal[:3])
M = numpy.identity(4)
M[:3, :3] -= 2.0 * numpy.outer(normal, normal)
M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal
return M
def update_kl_loss(p, lambdas, T, Cs):
"""
Updates C according to the KL Loss kernel with the S Ts couplings calculated at each iteration
Parameters
----------
p : ndarray, shape (N,)
weights in the targeted barycenter
lambdas : list of the S spaces' weights
T : list of S np.ndarray(ns,N)
the S Ts couplings calculated at each iteration
Cs : list of S ndarray, shape(ns,ns)
Metric cost matrices
Returns
----------
C : ndarray, shape (ns,ns)
updated C matrix
"""
tmpsum = sum([lambdas[s] * np.dot(T[s].T, Cs[s]).dot(T[s])
for s in range(len(T))])
ppt = np.outer(p, p)
return np.exp(np.divide(tmpsum, ppt))
def __init__(self, n):
self.degree = 2*n - 2
a, A = numpy.polynomial.legendre.leggauss(n)
w = numpy.outer((1 + a) * A, A)
x = numpy.outer(1-a, numpy.ones(n)) / 2
y = numpy.outer(1+a, 1-a) / 4
self.weights = w.reshape(-1) / 4
self.points = numpy.stack([x.reshape(-1), y.reshape(-1)]).T
self.bary = numpy.array([
self.points[:, 0],
self.points[:, 1],
1 - numpy.sum(self.points, axis=1)
]).T
return
def rotation_matrix(u, theta):
'''Return matrix that implements the rotation around the vector :math:`u`
by the angle :math:`\\theta`, cf.
https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle.
:param u: rotation vector
:param theta: rotation angle
'''
# Cross-product matrix.
cpm = numpy.array([[0.0, -u[2], u[1]],
[u[2], 0.0, -u[0]],
[-u[1], u[0], 0.0]])
c = numpy.cos(theta)
s = numpy.sin(theta)
R = numpy.eye(3) * c \
+ s * cpm \
+ (1.0 - c) * numpy.outer(u, u)
return R
def reflection_matrix(point, normal):
"""Return matrix to mirror at plane defined by point and normal vector.
>>> v0 = numpy.random.random(4) - 0.5
>>> v0[3] = 1.
>>> v1 = numpy.random.random(3) - 0.5
>>> R = reflection_matrix(v0, v1)
>>> numpy.allclose(2, numpy.trace(R))
True
>>> numpy.allclose(v0, numpy.dot(R, v0))
True
>>> v2 = v0.copy()
>>> v2[:3] += v1
>>> v3 = v0.copy()
>>> v2[:3] -= v1
>>> numpy.allclose(v2, numpy.dot(R, v3))
True
"""
normal = unit_vector(normal[:3])
M = numpy.identity(4)
M[:3, :3] -= 2.0 * numpy.outer(normal, normal)
M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal
return M
def test_basic(self):
vecs = [[1, 2, 3],
[2],
np.array([1, 2, 3]).reshape(1, -1),
np.array([1, 2, 3]).reshape(-1, 1)]
num_ones_list = [4, 1]
for vec in vecs:
for num_ones in num_ones_list:
A = OnesOuterVec(num_ones=num_ones, vec=vec)
M = np.outer([1]*num_ones, vec)
V, v1, v2, U, u1, u2 = get_tst_mats(M.shape)
assert_allclose(A.dot(V), M.dot(V))
assert_allclose(A.dot(v1), M.dot(v1))
assert_allclose(A.dot(v2), M.dot(v2))
assert_allclose(A.T.dot(U), M.T.dot(U))
assert_allclose(A.T.dot(u1), M.T.dot(u1))
assert_allclose(A.T.dot(u2), M.T.dot(u2))
def test_basic(self):
shapes = [(50, 20), (1, 20), (50, 1)]
# sparse
for shape in shapes:
mats = self.get_Xs(shape)
m = mats[0].mean(axis=0).A1
ones = np.ones(shape[0])
M = mats[0].toarray() - np.outer(ones, m)
for X in mats:
A = col_mean_centered(X)
V, v1, v2, U, u1, u2 = get_tst_mats(M.shape)
assert_almost_equal(A.dot(V), M.dot(V))
assert_almost_equal(A.dot(v1), M.dot(v1))
assert_almost_equal(A.dot(v2), M.dot(v2))
assert_almost_equal(A.T.dot(U), M.T.dot(U))
assert_almost_equal(A.T.dot(u1), M.T.dot(u1))
assert_almost_equal(A.T.dot(u2), M.T.dot(u2))
def rotation_matrix(axis, angle):
"""
Calculate a three dimensional rotation matrix for a rotation around
the given angle and axis.
@type axis: (3,) numpy array
@param angle: angle in radians
@type angle: float
@rtype: (3,3) numpy.array
"""
axis = numpy.asfarray(axis) / norm(axis)
assert axis.shape == (3,)
c = math.cos(angle)
s = math.sin(angle)
r = (1.0 - c) * numpy.outer(axis, axis)
r.flat[[0,4,8]] += c
r.flat[[5,6,1]] += s * axis
r.flat[[7,2,3]] -= s * axis
return r
def gower_matrix(X):
"""
Gower, J.C. (1966). Some distance properties of latent root
and vector methods used in multivariate analysis.
Biometrika 53: 325-338
@param X: ensemble coordinates
@type X: (m,n,k) numpy.array
@return: symmetric dissimilarity matrix
@rtype: (n,n) numpy.array
"""
X = numpy.asarray(X)
B = sum(numpy.dot(x, x.T) for x in X) / float(len(X))
b = B.mean(1)
bb = b.mean()
return (B - numpy.add.outer(b, b)) + bb
def nextfastpower(n):
"""Return the next integral power of small factors greater than the given
number. Specifically, return m such that
m >= n
m == 2**x * 3**y * 5**z
where x, y, and z are integers.
This is useful for ensuring fast FFT sizes.
From https://gist.github.com/bhawkins/4479607 (Brian Hawkins)
"""
if n < 7:
return max (n, 1)
# x, y, and z are all bounded from above by the formula of nextpower.
# Compute all possible combinations for powers of 3 and 5.
# (Not too many for reasonable FFT sizes.)
def power_series (x, base):
nmax = ceil (log (x) / log (base))
return np.logspace (0.0, nmax, num=nmax+1, base=base)
n35 = np.outer (power_series (n, 3.0), power_series (n, 5.0))
n35 = n35[n35<=n]
# Lump the powers of 3 and 5 together and solve for the powers of 2.
n2 = nextpower (n / n35)
return int (min (n2 * n35))
def ests_ll_quad(self, params):
"""
Calculate the loglikelihood given model parameters `params`.
This method uses Gaussian quadrature, and thus returns an *approximate*
integral.
"""
mu0, gamma0, err0 = np.split(params, 3)
x = np.tile(self.z, (self.cfg.QCOUNT, 1, 1)) # (QCOUNTXnhospXnmeas)
loc = mu0 + np.outer(QC1, gamma0)
loc = np.tile(loc, (self.n, 1, 1))
loc = np.transpose(loc, (1, 0, 2))
scale = np.tile(err0, (self.cfg.QCOUNT, self.n, 1))
zs = lpdf_3d(x=x, loc=loc, scale=scale)
w2 = np.tile(self.w, (self.cfg.QCOUNT, 1, 1))
wted = np.nansum(w2 * zs, axis=2).T # (nhosp X QCOUNT)
qh = np.tile(QC1, (self.n, 1)) # (nhosp X QCOUNT)
combined = wted + norm.logpdf(qh) # (nhosp X QCOUNT)
return logsumexp(np.nan_to_num(combined), b=QC2, axis=1) # (nhosp)
def reflection_matrix(point, normal):
"""Return matrix to mirror at plane defined by point and normal vector.
>>> v0 = numpy.random.random(4) - 0.5
>>> v0[3] = 1.0
>>> v1 = numpy.random.random(3) - 0.5
>>> R = reflection_matrix(v0, v1)
>>> numpy.allclose(2., numpy.trace(R))
True
>>> numpy.allclose(v0, numpy.dot(R, v0))
True
>>> v2 = v0.copy()
>>> v2[:3] += v1
>>> v3 = v0.copy()
>>> v2[:3] -= v1
>>> numpy.allclose(v2, numpy.dot(R, v3))
True
"""
normal = unit_vector(normal[:3])
M = numpy.identity(4)
M[:3, :3] -= 2.0 * numpy.outer(normal, normal)
M[:3, 3] = (2.0 * numpy.dot(point[:3], normal)) * normal
return M
def forward_prop_random_thru_post_mm(self, model, mx, vx, mu, Su):
Kuu_noiseless = compute_kernel(
2 * model.ls, 2 * model.sf, model.zu, model.zu)
Kuu = Kuu_noiseless + np.diag(jitter * np.ones((self.M, )))
# TODO: remove inv
Kuuinv = np.linalg.inv(Kuu)
A = np.dot(Kuuinv, mu)
Smm = Su + np.outer(mu, mu)
B_sto = np.dot(Kuuinv, np.dot(Smm, Kuuinv)) - Kuuinv
psi0 = np.exp(2.0 * model.sf)
psi1, psi2 = compute_psi_weave(
2 * model.ls, 2 * model.sf, mx, vx, model.zu)
mout = np.einsum('nm,md->nd', psi1, A)
Bpsi2 = np.einsum('ab,nab->n', B_sto, psi2)[:, np.newaxis]
vout = psi0 + Bpsi2 - mout**2
return mout, vout
