def transfer_color(content, style):
import scipy.linalg as sl
# Mean and covariance of content
content_mean = np.mean(content, axis = (0, 1))
content_diff = content - content_mean
content_diff = np.reshape(content_diff, (-1, content_diff.shape[2]))
content_covariance = np.matmul(content_diff.T, content_diff) / (content_diff.shape[0])
# Mean and covariance of style
style_mean = np.mean(style, axis = (0, 1))
style_diff = style - style_mean
style_diff = np.reshape(style_diff, (-1, style_diff.shape[2]))
style_covariance = np.matmul(style_diff.T, style_diff) / (style_diff.shape[0])
# Calculate A and b
A = np.matmul(sl.sqrtm(content_covariance), sl.inv(sl.sqrtm(style_covariance)))
b = content_mean - np.matmul(A, style_mean)
# Construct new style
new_style = np.reshape(style, (-1, style.shape[2])).T
new_style = np.matmul(A, new_style).T
new_style = np.reshape(new_style, style.shape)
new_style = new_style + b
return new_style
python类sqrtm()的实例源码
def test_frechet_classifier_distance_value(self):
"""Test that `frechet_classifier_distance` gives the correct value."""
np.random.seed(0)
# Make num_examples > num_features to ensure scipy's sqrtm function
# doesn't return a complex matrix.
test_pool_real_a = np.float32(np.random.randn(512, 256))
test_pool_gen_a = np.float32(np.random.randn(512, 256))
fid_op = _run_with_mock(gan_metrics.frechet_classifier_distance,
test_pool_real_a, test_pool_gen_a,
classifier_fn=lambda x: x)
with self.test_session() as sess:
actual_fid = sess.run(fid_op)
expected_fid = _expected_fid(test_pool_real_a, test_pool_gen_a)
self.assertAllClose(expected_fid, actual_fid, 0.0001)
def test_trace_sqrt_product_value(self):
"""Test that `trace_sqrt_product` gives the correct value."""
np.random.seed(0)
# Make num_examples > num_features to ensure scipy's sqrtm function
# doesn't return a complex matrix.
test_pool_real_a = np.float32(np.random.randn(512, 256))
test_pool_gen_a = np.float32(np.random.randn(512, 256))
cov_real = np.cov(test_pool_real_a, rowvar=False)
cov_gen = np.cov(test_pool_gen_a, rowvar=False)
trace_sqrt_prod_op = _run_with_mock(gan_metrics.trace_sqrt_product,
cov_real, cov_gen)
with self.test_session() as sess:
# trace_sqrt_product: tsp
actual_tsp = sess.run(trace_sqrt_prod_op)
expected_tsp = _expected_trace_sqrt_product(cov_real, cov_gen)
self.assertAllClose(actual_tsp, expected_tsp, 0.01)
def renderStarGauss(image, cov, mu, first, scale = 5):
num_circles = 3
num_points = 64
cov = sqrtm(cov)
num = num_circles * num_points
pos = np.ones((num, 2))
for c in range(num_circles):
r = c + 1
for p in range(num_points):
angle = p / num_points * 2 * np.pi
index = c * num_points + p
x = r * np.cos(angle)
y = r * np.sin(angle)
pos[index, 0] = x * cov[0, 0] + y * cov[0, 1] + mu[0]
pos[index, 1] = x * cov[1, 0] + y * cov[1, 1] + mu[1]
#image = image.copy()
#image = cv2.cvtColor(image, cv2.COLOR_GRAY2BGR)
if first:
image = cv2.resize(image, (0, 0), None, scale, scale, cv2.INTER_NEAREST)
for c in range(num_circles):
pts = np.array(pos[c * num_points:(c + 1) * num_points, :] * scale + scale / 2, np.int32)
pts = pts.reshape((-1,1,2))
cv2.polylines(image, [pts], True, (255, 0, 0))
return image
def evd_r0(Vorg,
no_order=None,
no_samples=1000,
output_l=None):
"""
if output_l is not None, outputs will be saved to output_l
output_l = [generated_samples]
"""
n = Vorg - np.mean(Vorg, axis=0)
y0 = np.mean(Vorg, axis=0)
Rnn = np.dot(n.T, n) / n.shape[0]
if no_order is not None:
w, v = np.linalg.eig(Rnn)
w, v = w[:no_order], v[:,:no_order]
Rnn = np.dot(np.dot(v, np.diag(w)), v.T)
new_n = np.random.randn(no_samples, Vorg.shape[1])
y = y0.reshape(1,-1) + np.dot(new_n, linalg.sqrtm(Rnn))
if output_l is not None:
output_l.append(y)
def evd(Vorg,
no_order=None,
no_samples=1000):
n = Vorg - np.mean(Vorg, axis=0)
y0 = np.mean(Vorg, axis=0)
Rnn = np.dot(n.T, n) / n.shape[0]
if no_order is not None:
w, v = np.linalg.eig(Rnn)
w, v = w[:no_order], v[:,:no_order]
Rnn = np.dot(np.dot(v, np.diag(w)), v.T)
Rnn = np.abs(Rnn)
new_n = np.random.randn(no_samples, Vorg.shape[1])
y = y0.reshape(1,-1) + np.dot(new_n, linalg.sqrtm(Rnn))
y = np.real(y)
return y
def Y(W, cell=None):
"""Returns the normalized wave function `Y` (eq. 6 in Ps2).
Args:
W (numpy.ndarray): wave function sample points.
cell (pydft.geometry.Cell): describing the unit cell and sampling
points.
"""
from scipy.linalg import sqrtm
Usq = sqrtm(np.linalg.inv(U(W, cell)))
return np.dot(W, Usq)
def get_coral_params(ds, dt, lam=1e-3):
ms = ds.mean(axis=0)
ds = ds - ms
mt = dt.mean(axis=0)
dt = dt - mt
cs = np.cov(ds.T) + lam * np.eye(ds.shape[1])
ct = np.cov(dt.T) + lam * np.eye(dt.shape[1])
sqrt = splg.sqrtm
w = sqrt(ct).dot(np.linalg.inv(sqrt(cs)))
b = mt - w.dot(ms.reshape(-1, 1)).ravel()
return w, b
def _expected_fid(real_imgs, gen_imgs):
m = np.mean(real_imgs, axis=0)
m_v = np.mean(gen_imgs, axis=0)
sigma = np.cov(real_imgs, rowvar=False)
sigma_v = np.cov(gen_imgs, rowvar=False)
sqcc = scp_linalg.sqrtm(np.dot(sigma, sigma_v))
mean = np.square(m - m_v).sum()
trace = np.trace(sigma + sigma_v - 2 * sqcc)
fid = mean + trace
return fid
def _expected_trace_sqrt_product(sigma, sigma_v):
return np.trace(scp_linalg.sqrtm(np.dot(sigma, sigma_v)))
# A dummy GraphDef string with the minimum number of Ops.
def compute_induced_kernel_matrix_on_data(self,data_x,data_y):
'''Z follows the same distribution as X; W follows that of Y.
The current data generating methods we use
generate X and Y at the same time. '''
size_induced_set = max(self.num_inducex,self.num_inducey)
#print "size_induce_set", size_induced_set
if self.data_generator is None:
subsample_idx = np.random.randint(self.num_samples, size=size_induced_set)
self.data_z = data_x[subsample_idx,:]
self.data_w = data_y[subsample_idx,:]
else:
self.data_z, self.data_w = self.data_generator(size_induced_set)
self.data_z[[range(self.num_inducex)],:]
self.data_w[[range(self.num_inducey)],:]
print 'Induce Set'
if self.kernelX_use_median:
sigmax = self.kernelX.get_sigma_median_heuristic(data_x)
self.kernelX.set_width(float(sigmax))
if self.kernelY_use_median:
sigmay = self.kernelY.get_sigma_median_heuristic(data_y)
self.kernelY.set_width(float(sigmay))
Kxz = self.kernelX.kernel(data_x,self.data_z)
Kzz = self.kernelX.kernel(self.data_z)
#R = inv(sqrtm(Kzz))
R = inv(sqrtm(Kzz + np.eye(np.shape(Kzz)[0])*10**(-6)))
phix = Kxz.dot(R)
Kyw = self.kernelY.kernel(data_y,self.data_w)
Kww = self.kernelY.kernel(self.data_w)
#S = inv(sqrtm(Kww))
S = inv(sqrtm(Kww + np.eye(np.shape(Kww)[0])*10**(-6)))
phiy = Kyw.dot(S)
return phix, phiy
def randorth(shape):
from scipy.linalg import sqrtm, inv
assert len(shape) == 2
w = np.random.normal(0, size=shape)
w = w.dot(inv(sqrtm(w.T.dot(w))))
return G.sharedf(w)
# Softmax function
def fit(self, X):
"""
Compute the Dynamics Modes Decomposition to the input data.
:param X: the input snapshots.
:type X: numpy.ndarray or iterable
"""
self._snapshots, self._snapshots_shape = self._col_major_2darray(X)
n_samples = self._snapshots.shape[1]
X = self._snapshots[:, :-1]
Y = self._snapshots[:, 1:]
X, Y = self._compute_tlsq(X, Y, self.tlsq_rank)
Uy, sy, Vy = self._compute_svd(Y, self.svd_rank)
Ux, sx, Vx = self._compute_svd(X, self.svd_rank)
if len(sy) != len(sx):
raise ValueError(
'The optimal truncation produced different number of singular'
'values for the X and Y matrix, please specify different'
'svd_rank'
)
# Backward operator
bAtilde = self._build_lowrank_op(Uy, sy, Vy, X)
# Forward operator
fAtilde = self._build_lowrank_op(Ux, sx, Vx, Y)
self._Atilde = sqrtm(fAtilde.dot(np.linalg.inv(bAtilde)))
self._eigs, self._modes = self._eig_from_lowrank_op(
self._Atilde, Y, Ux, sx, Vx, self.exact
)
self._b = self._compute_amplitudes(
self._modes, self._snapshots, self._eigs, self.opt
)
self.original_time = {'t0': 0, 'tend': n_samples - 1, 'dt': 1}
self.dmd_time = {'t0': 0, 'tend': n_samples - 1, 'dt': 1}
return self
def _recursive_builder(self, operation, gate_name, control_qubits, target_qubit):
"""Helper function used to define the controlled gate recursively. It uses the algorithm in
the reference above. Namely it recursively constructs a controlled gate by applying a
controlled square root of the gate, followed by a toffoli gate with len(control_qubits) - 1
controls, applying the controlled adjoint of the square root, another toffoli with
len(control_qubits) - 1 controls, and finally another controlled copy of the gate.
:param numpy.ndarray operation: The matrix for the unitary to be controlled.
:param String gate_name: The name for the gate being controlled.
:param Sequence control_qubits: The qubits that are the controls.
:param Qubit or Int target_qubit: The qubit that the gate should be applied to.
:return: The intermediate Program being built.
:rtype: Program
"""
control_true = np.kron(ONE_PROJECTION, operation)
control_false = np.kron(ZERO_PROJECTION, np.eye(2, 2))
control_root_true = np.kron(ONE_PROJECTION, sqrtm(operation, disp=True))
controlled_gate = control_true + control_false
controlled_root_gate = control_root_true + control_false
sqrt_name = self.format_gate_name(SQRT_PREFIX, gate_name)
controlled_subprogram = pq.Program()
control_gate = pq.Program()
# For the base case, we check to see if there is just one control qubit.
if len(control_qubits) == 1:
control_name = self.format_gate_name(CONTROL_PREFIX, gate_name)
control_gate = self._defgate(control_gate, control_name, controlled_gate)
control_gate.inst((control_name, control_qubits[0], target_qubit))
return control_gate
else:
control_sqrt_name = self.format_gate_name(CONTROL_PREFIX, sqrt_name)
control_gate = self._defgate(control_gate, control_sqrt_name, controlled_root_gate)
control_gate.inst((control_sqrt_name, control_qubits[-1], target_qubit))
# Here we recurse to build a toffoli gate on n - 1 of the qubits.
n_minus_one_toffoli = self._recursive_builder(NOT_GATE,
NOT_GATE_LABEL,
control_qubits[:-1],
control_qubits[-1])
# We recurse to build a controlled sqrt of the target_gate, excluding the last control.
n_minus_one_controlled_sqrt = self._recursive_builder(sqrtm(operation, disp=True),
sqrt_name,
control_qubits[:-1],
target_qubit)
controlled_subprogram += control_gate
controlled_subprogram += n_minus_one_toffoli
controlled_subprogram += control_gate.dagger()
controlled_subprogram += n_minus_one_toffoli
controlled_subprogram += n_minus_one_controlled_sqrt
return controlled_subprogram
def compute_foes_svd_irls(matches_all_frames_, homographies_refined):
"""
Compute initial set of FoEs.
"""
matches_all_frames = filter_features(matches_all_frames_, homographies_refined)
def estimate_foe_svd_norm_irls(points1, points2, init=None):
"""
Estimate focus of expansion using least squares
"""
A,b = points2system(points1[:,0],points1[:,1],points2[:,0],points2[:,1],norm=False)
T = np.c_[A,-b]
weights = np.ones((T.shape[0],1))
mad = lambda x : 1.48 * np.median(np.abs(x - np.median(x)))
for it in range(100):
Tw = T * np.sqrt(weights)
# Compute correction matrix
M = np.array([[0.0,0.0,0.0],[0.0,0.0,0.0],[0.0,0.0,0.0]])
for i,p1 in enumerate(points1):
M += weights[i]**2 * np.array([[1.0,0.0,-p1[0]],[0.0,1.0,-p1[1]],[-p1[0],-p1[1],p1[0]**2 + p1[1]**2]])
Mcor = sqrtm(np.linalg.inv(M))
X = Mcor.dot(Tw.T.dot(Tw)).dot(Mcor)
U,S,V = np.linalg.svd(X)
err = T.dot(Mcor).dot(V[-1,:])
sigma = mad(err)
weights[:,0] = 1.0 / (1 + (err/sigma)**2)
weights /= weights.sum()
ep1 = V[2,:]
ep1 = Mcor.dot(ep1)
ep1 = ep1[:2] / ep1[2]
return ep1
epipoles = []
n_frames = matches_all_frames.shape[2]
for i in range(1,n_frames):
pt1 = matches_all_frames[:,:,0]
pt2 = matches_all_frames[:,:,i]
H = homographies_refined[i-1]
pt2_H = remove_homography(pt2,H)
epipoles.append(estimate_foe_svd_norm_irls(pt1,pt2_H))
return epipoles
def conv_vh_decomposition(model, args):
W = model.arg_params[args.layer+'_weight'].asnumpy()
N, C, y, x = W.shape
b = model.arg_params[args.layer+'_bias'].asnumpy()
W = W.transpose((1,2,0,3)).reshape((C*y, -1))
U, D, Q = np.linalg.svd(W, full_matrices=False)
sqrt_D = LA.sqrtm(np.diag(D))
K = args.K
V = U[:,:K].dot(sqrt_D[:K, :K])
H = Q.T[:,:K].dot(sqrt_D[:K, :K])
V = V.T.reshape(K, C, y, 1)
b_1 = np.zeros((K, ))
H = H.reshape(N, x, 1, K).transpose((0,3,2,1))
b_2 = b
W1, b1, W2, b2 = V, b_1, H, b_2
def sym_handle(data, node):
kernel = eval(node['param']['kernel'])
pad = eval(node['param']['pad'])
name = node['name']
name1 = name + '_v'
kernel1 = tuple((kernel[0], 1))
pad1 = tuple((pad[0], 0))
num_filter = W1.shape[0]
sym1 = mx.symbol.Convolution(data=data, kernel=kernel1, pad=pad1, num_filter=num_filter, name=name1)
name2 = name + '_h'
kernel2 = tuple((1, kernel[1]))
pad2 = tuple((0, pad[1]))
num_filter = W2.shape[0]
sym2 = mx.symbol.Convolution(data=sym1, kernel=kernel2, pad=pad2, num_filter=num_filter, name=name2)
return sym2
def arg_handle(arg_shape_dic, arg_params):
name1 = args.layer + '_v'
name2 = args.layer + '_h'
weight1 = mx.ndarray.array(W1)
bias1 = mx.ndarray.array(b1)
weight2 = mx.ndarray.array(W2)
bias2 = mx.ndarray.array(b2)
assert weight1.shape == arg_shape_dic[name1+'_weight'], 'weight1'
assert weight2.shape == arg_shape_dic[name2+'_weight'], 'weight2'
assert bias1.shape == arg_shape_dic[name1+'_bias'], 'bias1'
assert bias2.shape == arg_shape_dic[name2+'_bias'], 'bias2'
arg_params[name1 + '_weight'] = weight1
arg_params[name1 + '_bias'] = bias1
arg_params[name2 + '_weight'] = weight2
arg_params[name2 + '_bias'] = bias2
new_model = utils.replace_conv_layer(args.layer, model, sym_handle, arg_handle)
return new_model
def g2s_weighted(x, y, Zx, Zy, Lxx, Lxy, Lyx, Lyy, N=3):
"""
% Purpose : Computes the Global Weighted Least Squares reconstruction of a
% surface from its gradient field, whereby the weighting is defined by a
% weighted Frobenius norm
%
% Use (syntax):
% Z = g2sWeighted( Zx, Zy, x, y, N, Lxx, Lxy, Lyx, Lyy )
%
% Input Parameters :
% Zx, Zy := Components of the discrete gradient field
% x, y := support vectors of nodes of the domain of the gradient
% N := number of points for derivative formulas (default=3)
% Lxx, Lxy, Lyx, Lyy := Each matrix Lij is the covariance matrix of the
% gradient's i-component the in j-direction.
%
% Return Parameters :
% Z := The reconstructed surface
%
% Description and algorithms:
% The algorithm solves the normal equations of the Weighted Least Squares
% cost function, formulated by matrix algebra:
% e(Z) = || Lyy^(-1/2) * (D_y * Z - Zy) * Lyx^(-1/2) ||_F^2 +
% || Lxy^(-1/2) * ( Z * Dx' - Zx ) * Lxx^(-1/2) ||_F^2
% The normal equations are a rank deficient Sylvester equation which is
% solved by means of Householder reflections and the Bartels-Stewart
% algorithm.
"""
if Zx.shape != Zy.shape:
raise ValueError("Gradient components must be the same size")
if np.asmatrix(Zx).shape[1] != len(x) or np.asmatrix(Zx).shape[0] != len(y):
raise ValueError("Support vectors must have the same size as the gradient")
m, n = Zx.shape
Dx = diff_local(x, N, N)
Dy = diff_local(y, N, N)
Wxx = la.sqrtm(Lxx)
Wxy = la.sqrtm(Lxy)
Wyx = la.sqrtm(Lyx)
Wyy = la.sqrtm(Lyy)
# Solution for Zw (written here Z)
u = mldivide(Wxy, np.ones((m, 1)))
v = mldivide(Wyx, np.ones((n, 1)))
A = mldivide(Wyy, np.dot(Dy, Wxy))
B = mldivide(Wxx, np.dot(Dx, Wyx))
F = mldivide(Wyy, mrdivide(Zy, Wyx))
G = mldivide(Wxy, mrdivide(Zx, Wxx))
Z = _g2sSylvester(A, B, F, G, u, v)
# "Unweight" the solution
Z = Wxy * Z * Wyx
return Z
