Python scipy.linalg 模块，sqrtm() 实例源码

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

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