def test_gmmmap_swap():
static_dim = 2
T = 10
windows = _get_windows_set()[-1]
np.random.seed(1234)
src_mc = np.random.rand(T, static_dim * len(windows))
tgt_mc = np.random.rand(T, static_dim * len(windows))
# pseudo parallel data
XY = np.concatenate((src_mc, tgt_mc), axis=-1)
gmm = GaussianMixture(n_components=4)
gmm.fit(XY)
paramgen = MLPG(gmm, windows=windows, swap=False)
swap_paramgen = MLPG(gmm, windows=windows, swap=True)
mc_converted1 = paramgen.transform(src_mc)
mc_converted2 = swap_paramgen.transform(tgt_mc)
src_mc = src_mc[:, :static_dim]
tgt_mc = tgt_mc[:, :static_dim]
assert norm(tgt_mc - mc_converted1) < norm(src_mc - mc_converted1)
assert norm(tgt_mc - mc_converted2) > norm(src_mc - mc_converted2)
python类norm()的实例源码
def kNN_entity(self, vec, topk=10, method=0, self_vec_id=None):
q = []
for i in range(len(self.vec_e)):
#skip self
if self_vec_id != None and i == self_vec_id:
continue
if method == 1:
dist = SP.distance.cosine(vec, self.vec_e[i])
else:
dist = LA.norm(vec - self.vec_e[i])
if len(q) < topk:
HP.heappush(q, self.index_dist(i, dist))
else:
#indeed it fetches the biggest
tmp = HP.nsmallest(1, q)[0]
if tmp.dist > dist:
HP.heapreplace(q, self.index_dist(i, dist) )
rst = []
while len(q) > 0:
item = HP.heappop(q)
rst.insert(0, (self.vocab_e[self.vec2e[item.index]], item.dist))
return rst
#given entity name, find kNN
def _ncc_c(x, y):
"""
>>> _ncc_c([1,2,3,4], [1,2,3,4])
array([ 0.13333333, 0.36666667, 0.66666667, 1. , 0.66666667,
0.36666667, 0.13333333])
>>> _ncc_c([1,1,1], [1,1,1])
array([ 0.33333333, 0.66666667, 1. , 0.66666667, 0.33333333])
>>> _ncc_c([1,2,3], [-1,-1,-1])
array([-0.15430335, -0.46291005, -0.9258201 , -0.77151675, -0.46291005])
"""
den = np.array(norm(x) * norm(y))
den[den == 0] = np.Inf
x_len = len(x)
fft_size = 1<<(2*x_len-1).bit_length()
cc = ifft(fft(x, fft_size) * np.conj(fft(y, fft_size)))
cc = np.concatenate((cc[-(x_len-1):], cc[:x_len]))
return np.real(cc) / den
def reshape(self, vertices):
"""Make snake smoother."""
arc_length = 0
# calculate the total arc-length
for i in range(1, self.length):
arc_length += norm(vertices[i] - vertices[i - 1])
# average length for each segment
seg_length = arc_length / (self.length - 1)
if self._bc == 'PBC':
arc_length += norm(vertices[0] - vertices[-1])
seg_length = arc_length / self.length
for i in range(1, self.length - 1 if self._bc == 'OBC' else self.length):
# normalized tangent direction at node i-1
tan_direction = vertices[i] - vertices[i - 1]
tan_direction /= norm(tan_direction)
# move node i
vertices[i] = vertices[i - 1] + tan_direction * seg_length
return vertices
def distribution_parameters(self, parameter_name):
if parameter_name=='trend':
E = dot(dot(self.derivative_matrix.T,inv(diag(self.parameters.list['omega'].current_value))),self.derivative_matrix)
mean = dot(inv(eye(self.size)+E),self.data)
cov = (self.parameters.list['sigma2'].current_value)*inv(eye(self.size)+E)
return {'mean' : mean, 'cov' : cov}
elif parameter_name=='sigma2':
E = dot(dot(self.derivative_matrix.T,inv(diag(self.parameters.list['omega'].current_value))),self.derivative_matrix)
pos = self.size
loc = 0
scale = 0.5*dot((self.data-dot(eye(self.size),self.parameters.list['trend'].current_value)).T,(self.data-dot(eye(self.size),self.parameters.list['trend'].current_value)))+0.5*dot(dot(self.parameters.list['trend'].current_value.T,E),self.parameters.list['trend'].current_value)
elif parameter_name=='lambda2':
pos = self.size-self.total_variation_order-1+self.alpha
loc = 0.5*(norm(dot(self.derivative_matrix,self.parameters.list['trend'].current_value),ord=1))/self.parameters.list['sigma2'].current_value+self.rho
scale = 1
elif parameter_name==str('omega'):
pos = [sqrt(((self.parameters.list['lambda2'].current_value**2)*self.parameters.list['sigma2'].current_value)/(dj**2)) for dj in dot(self.derivative_matrix,self.parameters.list['trend'].current_value)]
loc = 0
scale = self.parameters.list['lambda2'].current_value**2
return {'pos' : pos, 'loc' : loc, 'scale' : scale}
def initialize_vectors(self):
self.vec2e = np.zeros(len(self.vocab_e), dtype=np.int)
for i in range(len(self.vocab_e)):
w = self.vocab_e[i]
self.e2vec[w] = i
self.vec2e[i] = i
#randomize and normalize the initial vector
self.vec_e = RD.randn(len(self.vocab_e), self.dim)
for i in range(len(self.vec_e)):
self.vec_e[i] /= LA.norm(self.vec_e[i])
self.vec2r = np.zeros(len(self.vocab_r), dtype=np.int)
for i in range(len(self.vocab_r)):
w = self.vocab_r[i]
#randomize and normalize the initial vector
self.r2vec[w] = i
self.vec2r[i] = i
self.vec_r = RD.randn(len(self.vocab_r), self.dim)
for i in range(len(self.vec_r)):
self.vec_r[i] /= LA.norm(self.vec_r[i])
print "Initialized vectors for all entities and relations."
#GD of h + r - t
def kNN_relation(self, vec, topk=10, method=0, self_vec_id=None):
q = []
for i in range(len(self.vec_r)):
#skip self
if self_vec_id != None and i == self_vec_id:
continue
if method == 1:
dist = SP.distance.cosine(vec, self.vec_r[i])
else:
dist = LA.norm(vec - self.vec_r[i])
if len(q) < topk:
HP.heappush(q, self.index_dist(i, dist))
else:
#indeed it fetches the biggest
tmp = HP.nsmallest(1, q)[0]
if tmp.dist > dist:
HP.heapreplace(q, self.index_dist(i, dist) )
rst = []
while len(q) > 0:
item = HP.heappop(q)
rst.insert(0, (self.vocab_r[self.vec2r[item.index]], item.dist))
return rst
#given relation name, find kNN
def gradient_decent(self,tr,v_e1,v_e2,const_decay, L1=False):
# f = || v_e1 tr - v_e2 ||_2^2
assert len(v_e1.shape) == 1
assert v_e1.shape == v_e2.shape
assert tr.shape == (v_e1.shape[0], v_e2.shape[0])
f_res = np.dot(v_e1, tr) - v_e2
d_v_e1 = 2.0 * np.dot(tr, f_res)
d_v_e2 = - 2.0 * f_res
d_tr = 2.0 * np.dot(v_e1[:, np.newaxis], f_res[np.newaxis, :])
v_e1 -= d_v_e1 * self.rate
v_e2 -= d_v_e2 * self.rate
tr -= d_tr * self.rate
v_e1 /= LA.norm(v_e1)
v_e2 /= LA.norm(v_e2)
# don't touch tr
return LA.norm(np.dot(v_e1, tr) - v_e2)
def kNN_entity_name(self, name, src_lan='en', tgt_lan='fr', topk=10, method=0, replace=True):
model = self.models.get(src_lan)
if model == None:
print "Model for language", src_lan," does not exist."
return None
id = model.e2vec.get(name)
if id == None:
print name," is not in vocab"
return None
if src_lan != tgt_lan:#if you're not quering the kNN in the same language, then no need to get rid of the "self" point. However, transfer the vector.
pass_vec = np.dot(model.vec_e[id], self.transfer[src_lan + tgt_lan])
pass_vec /= LA.norm(pass_vec)
return self.kNN_entity(pass_vec, tgt_lan, topk, method, self_vec_id=None, replace_q=replace)
return self.kNN_entity(model.vec_e[id], tgt_lan, topk, method, self_vec_id=id, replace_q=replace)
#return k nearest relations to a given vector (np.array)
def kNN_relation_name(self, name, src_lan='en', tgt_lan='fr', topk=10, method=0):
model = self.models.get(src_lan)
if model == None:
print "Model for language", src_lan," does not exist."
return None
id = model.r2vec.get(name)
if id == None:
print name," is not in vocab"
return None
if src_lan != tgt_lan:#if you're not quering the kNN in the same language, then no need to get rid of the "self" point
pass_vec = np.dot(model.vec_r[id], self.transfer[src_lan + tgt_lan])
pass_vec /= LA.norm(pass_vec)
return self.kNN_relation(pass_vec, tgt_lan, topk, method, self_vec_id=None)
return self.kNN_relation(model.vec_r[id], tgt_lan, topk, method, self_vec_id=id)
#entity name to vector
def test_matrix_3x3(self):
# This test has been added because the 2x2 example
# happened to have equal nuclear norm and induced 1-norm.
# The 1/10 scaling factor accommodates the absolute tolerance
# used in assert_almost_equal.
A = (1 / 10) * \
np.array([[1, 2, 3], [6, 0, 5], [3, 2, 1]], dtype=self.dt)
assert_almost_equal(norm(A), (1 / 10) * 89 ** 0.5)
assert_almost_equal(norm(A, 'fro'), (1 / 10) * 89 ** 0.5)
assert_almost_equal(norm(A, 'nuc'), 1.3366836911774836)
assert_almost_equal(norm(A, inf), 1.1)
assert_almost_equal(norm(A, -inf), 0.6)
assert_almost_equal(norm(A, 1), 1.0)
assert_almost_equal(norm(A, -1), 0.4)
assert_almost_equal(norm(A, 2), 0.88722940323461277)
assert_almost_equal(norm(A, -2), 0.19456584790481812)
def test_bad_args(self):
# Check that bad arguments raise the appropriate exceptions.
A = array([[1, 2, 3], [4, 5, 6]], dtype=self.dt)
B = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4)
# Using `axis=<integer>` or passing in a 1-D array implies vector
# norms are being computed, so also using `ord='fro'`
# or `ord='nuc'` raises a ValueError.
assert_raises(ValueError, norm, A, 'fro', 0)
assert_raises(ValueError, norm, A, 'nuc', 0)
assert_raises(ValueError, norm, [3, 4], 'fro', None)
assert_raises(ValueError, norm, [3, 4], 'nuc', None)
# Similarly, norm should raise an exception when ord is any finite
# number other than 1, 2, -1 or -2 when computing matrix norms.
for order in [0, 3]:
assert_raises(ValueError, norm, A, order, None)
assert_raises(ValueError, norm, A, order, (0, 1))
assert_raises(ValueError, norm, B, order, (1, 2))
# Invalid axis
assert_raises(ValueError, norm, B, None, 3)
assert_raises(ValueError, norm, B, None, (2, 3))
assert_raises(ValueError, norm, B, None, (0, 1, 2))
def signal_by_db(x1, x2, snr, handle_method):
x1 = x1.astype(np.int32)
x2 = x2.astype(np.int32)
l1 = x1.shape[0]
l2 = x2.shape[0]
if l1 != l2:
if handle_method == 'cut':
ll = min(l1, l2)
x1 = x1[:ll]
x2 = x2[:ll]
elif handle_method == 'append':
ll = max(l1, l2)
if l1 < ll:
x1 = np.append(x1, x1[:ll-l1])
if l2 < ll:
x2 = np.append(x2, x2[:ll-l1])
from numpy.linalg import norm
x2 = x2 / norm(x2) * norm(x1) / (10.0 ** (0.05 * snr))
mix = x1 + x2
return mix
def calc_rotation_matrix(q1, q2, ref_q1, ref_q2):
ref_nv = np.cross(ref_q1, ref_q2)
q_nv = np.cross(q1, q2)
if min(norm(ref_nv), norm(q_nv)) == 0.: # avoid 0 degree including angle
return np.identity(3)
axis = np.cross(ref_nv, q_nv)
angle = rad2deg(acos(ref_nv.dot(q_nv) / (norm(ref_nv) * norm(q_nv))))
R1 = axis_angle_to_rotation_matrix(axis, angle)
rot_ref_q1, rot_ref_q2 = R1.dot(ref_q1), R1.dot(ref_q2) # rotate ref_q1,2 plane to q1,2 plane
cos1 = max(min(q1.dot(rot_ref_q1) / (norm(rot_ref_q1) * norm(q1)), 1.), -1.) # avoid math domain error
cos2 = max(min(q2.dot(rot_ref_q2) / (norm(rot_ref_q2) * norm(q2)), 1.), -1.)
angle1 = rad2deg(acos(cos1))
angle2 = rad2deg(acos(cos2))
angle = (angle1 + angle2) / 2.
axis = np.cross(rot_ref_q1, q1)
R2 = axis_angle_to_rotation_matrix(axis, angle)
R = R2.dot(R1)
return R
def get_lipschitz(data):
"""Get the Lipschitz constant for a specific loss function.
Only square loss implemented.
Parameters
----------
data : (n, d) float ndarray
data matrix
loss : string
the selected loss function in {'square', 'logit'}
Returns
----------
L : float
the Lipschitz constant
"""
n, p = data.shape
if p > n:
tmp = np.dot(data, data.T)
else:
tmp = np.dot(data.T, data)
return la.norm(tmp, 2)
def prox_l1(w, alpha):
r"""Proximity operator for l1 norm.
:math:`\\hat{\\alpha}_{l,m} = sign(u_{l,m})\\left||u_{l,m}| - \\alpha \\right|_+`
Parameters
----------
u : ndarray
The vector (of the n-dimensional space) on witch we want
to compute the proximal operator
alpha : float
regularisation parameter
Returns
-------
ndarray : the vector corresponding to the application of the
proximity operator to u
"""
return np.sign(w) * np.maximum(np.abs(w) - alpha, 0.)
def __convert_survey_to_sequence(self):
s = self.__beamline
if 'LENGTH' not in s:
s['LENGTH'] = np.nan
offset = s['ORBIT_LENGTH'][0] / 2.0
if pd.isnull(offset):
offset = 0
self.__beamline['AT_CENTER'] = pd.DataFrame(
npl.norm(
[
s['X'].diff().fillna(0.0),
s['Y'].diff().fillna(0.0)
],
axis=0
) - (
s['LENGTH'].fillna(0.0) / 2.0 - s['ORBIT_LENGTH'].fillna(0.0) / 2.0
) + (
s['LENGTH'].shift(1).fillna(0.0) / 2.0 - s['ORBIT_LENGTH'].shift(1).fillna(0.0) / 2.0
)).cumsum() / 1000.0 + offset
self.__converted_from_survey = True
def _initialize(self):
logger.debug("Initializing Policy.")
# check if policy is already initialized by the user
if self.policy.initialized:
logger.debug("Use pre-set policy parameters.")
return self.policy.parameters
# outerwise draw an element at random from the parameter space
parameter = self.parameter_space.sample()
for _ in range(1000):
self.policy.parameters = parameter
grad = self.estimator(self.policy)
if (norm(grad) >= 1000 * self.eps):
return parameter
parameter = self.parameter_space.sample()
logger.error('Unable to find non-zero gradient.')
def is_feasible(x, u):
# Reject going too fast
for i, v in enumerate(x[3:]):
if v > velmax_pos_plan[i] or v < velmax_neg_plan[i]:
return False
# Body to world
R = np.array([
[np.cos(x[2]), -np.sin(x[2])],
[np.sin(x[2]), np.cos(x[2])],
])
# Boat vertices in world frame
verts = x[:2] + R.dot(vps).T
# Check for collisions over all obstacles
for ob in obs:
if np.any(npl.norm(verts - ob[:2], axis=1) <= ob[2]):
return False
return True
################################################# HEURISTICS
def test_matrix_3x3(self):
# This test has been added because the 2x2 example
# happened to have equal nuclear norm and induced 1-norm.
# The 1/10 scaling factor accommodates the absolute tolerance
# used in assert_almost_equal.
A = (1 / 10) * \
np.array([[1, 2, 3], [6, 0, 5], [3, 2, 1]], dtype=self.dt)
assert_almost_equal(norm(A), (1 / 10) * 89 ** 0.5)
assert_almost_equal(norm(A, 'fro'), (1 / 10) * 89 ** 0.5)
assert_almost_equal(norm(A, 'nuc'), 1.3366836911774836)
assert_almost_equal(norm(A, inf), 1.1)
assert_almost_equal(norm(A, -inf), 0.6)
assert_almost_equal(norm(A, 1), 1.0)
assert_almost_equal(norm(A, -1), 0.4)
assert_almost_equal(norm(A, 2), 0.88722940323461277)
assert_almost_equal(norm(A, -2), 0.19456584790481812)
def test_bad_args(self):
# Check that bad arguments raise the appropriate exceptions.
A = array([[1, 2, 3], [4, 5, 6]], dtype=self.dt)
B = np.arange(1, 25, dtype=self.dt).reshape(2, 3, 4)
# Using `axis=<integer>` or passing in a 1-D array implies vector
# norms are being computed, so also using `ord='fro'`
# or `ord='nuc'` raises a ValueError.
assert_raises(ValueError, norm, A, 'fro', 0)
assert_raises(ValueError, norm, A, 'nuc', 0)
assert_raises(ValueError, norm, [3, 4], 'fro', None)
assert_raises(ValueError, norm, [3, 4], 'nuc', None)
# Similarly, norm should raise an exception when ord is any finite
# number other than 1, 2, -1 or -2 when computing matrix norms.
for order in [0, 3]:
assert_raises(ValueError, norm, A, order, None)
assert_raises(ValueError, norm, A, order, (0, 1))
assert_raises(ValueError, norm, B, order, (1, 2))
# Invalid axis
assert_raises(ValueError, norm, B, None, 3)
assert_raises(ValueError, norm, B, None, (2, 3))
assert_raises(ValueError, norm, B, None, (0, 1, 2))
def test_polygon(self):
pts = [[1,0], [1,1], [0,1], [0,2], [6,2]]
c = CurveFactory.polygon(pts)
expected_knots = [0,0,1,2,3,9,9]
actual_knots = c.knots(0,True)
self.assertEqual(len(c), 5)
self.assertEqual(c.order(0), 2)
self.assertEqual(c.dimension, 2)
self.assertAlmostEqual(np.linalg.norm(expected_knots - actual_knots), 0.0)
c = CurveFactory.polygon([0,0], [1,0], [0,1], [-1,0], relative=True)
self.assertEqual(len(c), 4)
self.assertAlmostEqual(c[2][0], 1)
self.assertAlmostEqual(c[2][1], 1)
self.assertAlmostEqual(c[3][0], 0)
self.assertAlmostEqual(c[3][1], 1)
def test_bezier(self):
crv = CurveFactory.bezier([[0,0], [0,1], [1,1], [1,0], [2,0], [2,1],[1,1]])
self.assertEqual(len(crv.knots(0)), 3)
self.assertTrue(np.allclose(crv(0), [0,0]))
t = crv.tangent(0)
self.assertTrue(np.allclose(t/norm(t), [0,1]))
t = crv.tangent(1, above=False)
self.assertTrue(np.allclose(t/norm(t), [0,-1]))
t = crv.tangent(1, above=True)
self.assertTrue(np.allclose(t/norm(t), [1,0]))
self.assertTrue(np.allclose(crv(1), [1,0]))
self.assertTrue(crv.order(0), 4)
# test the exact same curve, only with relative keyword
crv = CurveFactory.bezier([[0,0], [0,1], [1,0], [0,-1], [1,0], [0,1],[1,0]], relative=True)
self.assertEqual(len(crv.knots(0)), 3)
self.assertTrue(np.allclose(crv(0), [0,0]))
t = crv.tangent(0)
self.assertTrue(np.allclose(t/norm(t), [0,1]))
t = crv.tangent(1, above=False)
self.assertTrue(np.allclose(t/norm(t), [0,-1]))
t = crv.tangent(1, above=True)
self.assertTrue(np.allclose(t/norm(t), [1,0]))
self.assertTrue(np.allclose(crv(1), [1,0]))
self.assertTrue(crv.order(0), 4)
def test_volume_sphere(self):
# unit ball
ball = VolumeFactory.sphere(type='square')
# test 7x7 grid for radius = 1
for surf in ball.faces():
for u in np.linspace(surf.start(0), surf.end(0), 7):
for v in np.linspace(surf.start(1), surf.end(1), 7):
self.assertAlmostEqual(np.linalg.norm(surf(u, v), 2), 1.0)
self.assertAlmostEqual(surf.area(), 4*pi/6, places=3)
# unit ball
ball = VolumeFactory.sphere(type='radial')
# test 7x7 grid for radius = 1
for u in np.linspace(surf.start(0), surf.end(0), 7):
for v in np.linspace(surf.start(1), surf.end(1), 7):
self.assertAlmostEqual(np.linalg.norm(ball(u, v, 0), 2), 1.0) # w=0 is outer shell
self.assertAlmostEqual(np.linalg.norm(ball(u, v, 1), 2), 0.0) # w=1 is the degenerate core
self.assertAlmostEqual(ball.faces()[4].area(), 4*pi, places=3)
def test_cylinder_surface(self):
# unit cylinder
surf = SurfaceFactory.cylinder()
# test 7x7 grid for xy-radius = 1 and v=z
for u in np.linspace(surf.start()[0], surf.end()[0], 7):
for v in np.linspace(surf.start()[1], surf.end()[1], 7):
x = surf(u, v)
self.assertAlmostEqual(np.linalg.norm(x[:2], 2), 1.0) # (x,y) coordinates to z-axis
self.assertAlmostEqual(x[2], v) # z coordinate should be linear
self.assertAlmostEqual(surf.area(), 2*pi, places=3)
# cylinder with parameters
surf = SurfaceFactory.cylinder(r=2, h=5, center=(0,0,1), axis=(1,0,0))
for u in np.linspace(surf.start()[0], surf.end()[0], 7):
for v in np.linspace(surf.start()[1], surf.end()[1], 7):
x = surf(u, v)
self.assertAlmostEqual(x[1]**2+(x[2]-1)**2, 2.0**2) # distance to (z-1)=y=0
self.assertAlmostEqual(x[0], v*5) # x coordinate should be linear
self.assertAlmostEqual(surf.area(), 2*2*pi*5, places=3)
def test_cylinder_volume(self):
# unit cylinder
vol = VolumeFactory.cylinder()
# test 5x5x5 grid for xy-radius = w and v=z
for u in np.linspace(vol.start()[0], vol.end()[0], 5):
for v in np.linspace(vol.start()[1], vol.end()[1], 5):
for w in np.linspace(vol.start()[2], vol.end()[2], 5):
x = vol(u, v, w)
self.assertAlmostEqual(
np.linalg.norm(x[:2], 2), u) # (x,y) coordinates to z-axis
self.assertAlmostEqual(x[2], w) # z coordinate should be linear
self.assertAlmostEqual(vol.volume(), pi, places=3)
# cylinder with parameters
vol = VolumeFactory.cylinder(r=2, h=5, center=(0,0,1), axis=(1,0,0))
for u in np.linspace(vol.start()[0], vol.end()[0], 5):
for v in np.linspace(vol.start()[1], vol.end()[1], 5):
for w in np.linspace(vol.start()[2], vol.end()[2], 5):
x = vol(u, v, w)
self.assertAlmostEqual(x[1]**2+(x[2]-1)**2, u**2)
self.assertAlmostEqual(x[0], w*5)
self.assertAlmostEqual(vol.volume(), 2*2*pi*5, places=3)
def eval(self, coords, grad=False):
v1 = (coords[self.i]-coords[self.j])/bohr
v2 = (coords[self.k]-coords[self.j])/bohr
dot_product = np.dot(v1, v2)/(norm(v1)*norm(v2))
if dot_product < -1:
dot_product = -1
elif dot_product > 1:
dot_product = 1
phi = np.arccos(dot_product)
if not grad:
return phi
if abs(phi) > pi-1e-6:
grad = [
(pi-phi)/(2*norm(v1)**2)*v1,
(1/norm(v1)-1/norm(v2))*(pi-phi)/(2*norm(v1))*v1,
(pi-phi)/(2*norm(v2)**2)*v2
]
else:
grad = [
1/np.tan(phi)*v1/norm(v1)**2-v2/(norm(v1)*norm(v2)*np.sin(phi)),
(v1+v2)/(norm(v1)*norm(v2)*np.sin(phi)) -
1/np.tan(phi)*(v1/norm(v1)**2+v2/norm(v2)**2),
1/np.tan(phi)*v2/norm(v2)**2-v1/(norm(v1)*norm(v2)*np.sin(phi))
]
return phi, grad
def quadratic_step(g, H, w, trust, log=no_log):
ev = np.linalg.eigvalsh((H+H.T)/2)
rfo = np.vstack((np.hstack((H, g[:, None])),
np.hstack((g, 0))[None, :]))
D, V = np.linalg.eigh((rfo+rfo.T)/2)
dq = V[:-1, 0]/V[-1, 0]
l = D[0]
if norm(dq) <= trust:
log('Pure RFO step was performed:')
on_sphere = False
else:
def steplength(l):
return norm(np.linalg.solve(l*eye(H.shape[0])-H, g))-trust
l = Math.findroot(steplength, ev[0]) # minimization on sphere
dq = np.linalg.solve(l*eye(H.shape[0])-H, g)
on_sphere = False
log('Minimization on sphere was performed:')
dE = dot(g, dq)+0.5*dq.dot(H).dot(dq) # predicted energy change
log('* Trust radius: {:.2}'.format(trust))
log('* Number of negative eigenvalues: {}'.format((ev < 0).sum()))
log('* Lowest eigenvalue: {:.3}'.format(ev[0]))
log('* lambda: {:.3}'.format(l))
log('Quadratic step: RMS: {:.3}, max: {:.3}'.format(Math.rms(dq), max(abs(dq))))
log('* Predicted energy change: {:.3}'.format(dE))
return dq, dE, on_sphere
def XB(Y,index_PQ, index_P, n_PQ, n_P):
case_number, _ = np.shape(Y)
G, B = iterator.util.get_G_B(Y)
X = np.zeros((case_number, case_number))
B_p = np.zeros((n_P,n_P))
B_pp = np.zeros((n_PQ,n_PQ))
for i in xrange(case_number):
for j in xrange(case_number):
if LA.norm(Y[i][j]) > 1e-5 and i != j:
X[i][j] = np.reciprocal(np.imag(np.reciprocal(Y[i][j])))
for i in xrange(case_number):
for j in xrange(case_number):
if i != j:
X[i][i] -= X[i][j]
for i in xrange(0, n_P):
for j in xrange(0, n_P):
B_p[i][j] = X[index_P[i]][index_P[j]]
#---------------------------------------------------
for i in xrange(0, n_PQ):
for j in xrange(0, n_PQ):
B_pp[i][j] = B[index_PQ[i]][index_PQ[j]]
return B_p, B_pp
def predict(self, X):
X = self._pre_processing_x(X)
Y = np.zeros((X.shape[0], self.n_class))
A = self.A
# because X is (N x p), A is (K x p), we can to get the X_star (NxK)
X_star = X @ A.T
for k in range(self.n_class):
# mu_s_star shape is (p,)
mu_k_star = A @ self.Mu[k]
# Ref: http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.norm.html
# Ref: http://stackoverflow.com/questions/1401712/how-can-the-euclidean-distance-be-calculated-with-numpy
Y[:, k] = LA.norm(X_star - mu_k_star, axis=1) * 0.5 - log(self.Pi[k])
# Python index start from 0, transform to start with 1
y_hat = Y.argmin(axis=1).reshape((-1, 1)) + 1
return y_hat
