Python numpy.eye() Examples

def betaseries_file(tmpdir_factory,
                    deriv_betaseries_fname=deriv_betaseries_fname):
    bfile = tmpdir_factory.mktemp("beta").ensure(deriv_betaseries_fname)
    np.random.seed(3)
    num_trials = 40
    tgt_corr = 0.1
    bs1 = np.random.rand(num_trials)
    # create another betaseries with a target correlation
    bs2 = minimize(lambda x: abs(tgt_corr - pearsonr(bs1, x)[0]),
                   np.random.rand(num_trials)).x

    # two identical beta series
    bs_data = np.array([[[bs1, bs2]]])

    # the nifti image
    bs_img = nib.Nifti1Image(bs_data, np.eye(4))
    bs_img.to_filename(str(bfile))

    return bfile 

def time_align_visualize(alignments, time, y, namespace='time_align'):
    plt.figure()
    heat = np.flip(alignments + alignments.T +
                   np.eye(alignments.shape[0]), axis=0)
    sns.heatmap(heat, cmap="YlGnBu", vmin=0, vmax=1)
    plt.savefig(namespace + '_heatmap.svg')

    G = nx.from_numpy_matrix(alignments)
    G = nx.maximum_spanning_tree(G)

    pos = {}
    for i in range(len(G.nodes)):
        pos[i] = np.array([time[i], y[i]])

    mst_edges = set(nx.maximum_spanning_tree(G).edges())
    
    weights = [ G[u][v]['weight'] if (not (u, v) in mst_edges) else 8
                for u, v in G.edges() ]
    
    plt.figure()
    nx.draw(G, pos, edges=G.edges(), width=10)
    plt.ylim([-1, 1])
    plt.savefig(namespace + '.svg') 

def test_givens_inverse():
    r"""
    The Givens rotation in OpenFermion is defined as

    .. math::

        \begin{pmatrix}
            \cos(\theta) & -e^{i \varphi} \sin(\theta) \\
            \sin(\theta) &     e^{i \varphi} \cos(\theta)
        \end{pmatrix}.

    confirm numerically its hermitian conjugate is it's inverse
    """
    a = numpy.random.random() + 1j * numpy.random.random()
    b = numpy.random.random() + 1j * numpy.random.random()
    ab_rotation = givens_matrix_elements(a, b, which='right')

    assert numpy.allclose(ab_rotation.dot(numpy.conj(ab_rotation).T),
                          numpy.eye(2))
    assert numpy.allclose(numpy.conj(ab_rotation).T.dot(ab_rotation),
                          numpy.eye(2)) 

def fit(self, Xs, Xt):
        '''
        Transform Xs and Xt
        :param Xs: ns * n_feature, source feature
        :param Xt: nt * n_feature, target feature
        :return: Xs_new and Xt_new after TCA
        '''
        X = np.hstack((Xs.T, Xt.T))
        X /= np.linalg.norm(X, axis=0)
        m, n = X.shape
        ns, nt = len(Xs), len(Xt)
        e = np.vstack((1 / ns * np.ones((ns, 1)), -1 / nt * np.ones((nt, 1))))
        M = e * e.T
        M = M / np.linalg.norm(M, 'fro')
        H = np.eye(n) - 1 / n * np.ones((n, n))
        K = kernel(self.kernel_type, X, None, gamma=self.gamma)
        n_eye = m if self.kernel_type == 'primal' else n
        a, b = np.linalg.multi_dot([K, M, K.T]) + self.lamb * np.eye(n_eye), np.linalg.multi_dot([K, H, K.T])
        w, V = scipy.linalg.eig(a, b)
        ind = np.argsort(w)
        A = V[:, ind[:self.dim]]
        Z = np.dot(A.T, K)
        Z /= np.linalg.norm(Z, axis=0)
        Xs_new, Xt_new = Z[:, :ns].T, Z[:, ns:].T
        return Xs_new, Xt_new 

def __init__(self):
        """Initialize variable used by Kalman Filter class
        Args:
            None
        Return:
            None
        """
        self.dt = 0.005  # delta time

        self.A = np.array([[1, 0], [0, 1]])  # matrix in observation equations
        self.u = np.zeros((2, 1))  # previous state vector

        # (x,y) tracking object center
        self.b = np.array([[0], [255]])  # vector of observations

        self.P = np.diag((3.0, 3.0))  # covariance matrix
        self.F = np.array([[1.0, self.dt], [0.0, 1.0]])  # state transition mat

        self.Q = np.eye(self.u.shape[0])  # process noise matrix
        self.R = np.eye(self.b.shape[0])  # observation noise matrix
        self.lastResult = np.array([[0], [255]]) 

def fit(self, Xs, Xt):
        '''
        Fit source and target using KMM (compute the coefficients)
        :param Xs: ns * dim
        :param Xt: nt * dim
        :return: Coefficients (Pt / Ps) value vector (Beta in the paper)
        '''
        ns = Xs.shape[0]
        nt = Xt.shape[0]
        if self.eps == None:
            self.eps = self.B / np.sqrt(ns)
        K = kernel(self.kernel_type, Xs, None, self.gamma)
        kappa = np.sum(kernel(self.kernel_type, Xs, Xt, self.gamma) * float(ns) / float(nt), axis=1)

        K = matrix(K)
        kappa = matrix(kappa)
        G = matrix(np.r_[np.ones((1, ns)), -np.ones((1, ns)), np.eye(ns), -np.eye(ns)])
        h = matrix(np.r_[ns * (1 + self.eps), ns * (self.eps - 1), self.B * np.ones((ns,)), np.zeros((ns,))])

        sol = solvers.qp(K, -kappa, G, h)
        beta = np.array(sol['x'])
        return beta 

def netflix(es, ps, e0, l=.0001):
    """Combine predictions with the optimal weights to minimize RMSE.

    Args:
        es (list of float): RMSEs of predictions
        ps (list of np.array): predictions
        e0 (float): RMSE of all zero prediction
        l (float): lambda as in the ridge regression

    Returns:
        (tuple):

            - (np.array): ensemble predictions
            - (np.array): weights for input predictions
    """
    m = len(es)
    n = len(ps[0])

    X = np.stack(ps).T
    pTy = .5 * (n * e0**2 + (X**2).sum(axis=0) - n * np.array(es)**2)

    w = np.linalg.pinv(X.T.dot(X) + l * n * np.eye(m)).dot(pTy)

    return X.dot(w), w 

def invertFast(A, d):
	"""
	given an array A of shape d x k and a d x 1 vector d, computes (A * A.T + diag(d)) ^{-1}
	Checked.
	"""
	assert(A.shape[0] == d.shape[0])
	assert(d.shape[1] == 1)

	k = A.shape[1]
	A = np.array(A)
	d_vec = np.array(d)
	d_inv = np.array(1 / d_vec[:, 0])

	inv_d_squared = np.dot(np.atleast_2d(d_inv).T, np.atleast_2d(d_inv))
	M = np.diag(d_inv) - inv_d_squared * np.dot(np.dot(A, np.linalg.inv(np.eye(k, k) + np.dot(A.T, mult_diag(d_inv, A)))), A.T)

	return M 

def _N(self,s,r):
        """
        Lagrange's interpolate function
        params:
            s,r:natural position of evalue point.2-array.
        returns:
            2x(2x4) shape function matrix.
        """
        la1=(1-s)/2
        la2=(1+s)/2
        lb1=(1-r)/2
        lb2=(1+r)/2
        N1=la1*lb1
        N2=la1*lb2
        N3=la2*lb1
        N4=la2*lb2

        N=np.hstack(N1*np.eye(2),N2*np.eye(2),N3*np.eye(2),N4*np.eye(2))
        return N 

def energy_from_opdm(opdm, constant, one_body_tensor, two_body_tensor):
    """
    Evaluate the energy of an opdm assuming the 2-RDM is opdm ^ opdm

    :param opdm: single spin-component of the full spin-orbital opdm.
    :param constant: constant shift to the Hamiltonian. Commonly this is the
                     nuclear repulsion energy.
    :param one_body_tensor: spatial one-body integrals
    :param two_body_tensor: spatial two-body integrals
    :return:
    """
    spin_opdm = np.kron(opdm, np.eye(2))
    spin_tpdm = 2 * wedge(spin_opdm, spin_opdm, (1, 1), (1, 1))
    molecular_hamiltonian = generate_hamiltonian(constant=constant,
                                                 one_body_integrals=one_body_tensor,
                                                 two_body_integrals=two_body_tensor)
    rdms = InteractionRDM(spin_opdm, spin_tpdm)
    return rdms.expectation(molecular_hamiltonian).real 

def test_output():
    shape = (2,2)
    ones = mx.nd.ones(shape)
    zeros = mx.nd.zeros(shape)
    out = mx.nd.zeros(shape)
    mx.nd.ones(shape, out=out)
    assert_almost_equal(out.asnumpy(), ones.asnumpy())
    mx.nd.zeros(shape, out=out)
    assert_almost_equal(out.asnumpy(), zeros.asnumpy())
    mx.nd.full(shape, 2, out=out)
    assert_almost_equal(out.asnumpy(), ones.asnumpy() * 2)
    arange_out = mx.nd.arange(0, 20, dtype='int64')
    assert_almost_equal(arange_out.asnumpy(), np.arange(0, 20))
    N_array = np.random.randint(1, high=8, size=10)
    M_array = np.random.randint(1, high=8, size=10)
    k_array = np.random.randint(-10, high=10, size=10)
    for i in range(10):
        N = N_array[i]
        M = M_array[i]
        k = k_array[i]
        assert_almost_equal(np.eye(N, M, k), mx.nd.eye(N, M, k).asnumpy())
        assert_almost_equal(np.eye(N, k=k), mx.nd.eye(N, k=k).asnumpy()) 

def mtx_updated_G(phi_recon, M, mtx_amp2visi_ri, mtx_fri2visi_ri):
    """
    Update the linear transformation matrix that links the FRI sequence to the
    visibilities by using the reconstructed Dirac locations.
    :param phi_recon: the reconstructed Dirac locations (azimuths)
    :param M: the Fourier series expansion is between -M to M
    :param p_mic_x: a vector that contains microphones' x-coordinates
    :param p_mic_y: a vector that contains microphones' y-coordinates
    :param mtx_freq2visi: the linear mapping from Fourier series to visibilities
    :return:
    """
    L = 2 * M + 1
    ms_half = np.reshape(np.arange(-M, 1, step=1), (-1, 1), order='F')
    phi_recon = np.reshape(phi_recon, (1, -1), order='F')
    mtx_amp2freq = np.exp(-1j * ms_half * phi_recon)  # size: (M + 1) x K
    mtx_amp2freq_ri = np.vstack((mtx_amp2freq.real, mtx_amp2freq.imag[:-1, :]))  # size: (2M + 1) x K
    mtx_fri2amp_ri = linalg.lstsq(mtx_amp2freq_ri, np.eye(L))[0]
    # projection mtx_freq2visi to the null space of mtx_fri2amp
    mtx_null_proj = np.eye(L) - np.dot(mtx_fri2amp_ri.T,
                                       linalg.lstsq(mtx_fri2amp_ri.T, np.eye(L))[0])
    G_updated = np.dot(mtx_amp2visi_ri, mtx_fri2amp_ri) + \
                np.dot(mtx_fri2visi_ri, mtx_null_proj)
    return G_updated 

def set_camera(self, fov_vertical, z_near, z_far, aspect):
    width = 2*np.tan(np.deg2rad(fov_vertical)/2.0)*z_near*aspect;
    height = 2*np.tan(np.deg2rad(fov_vertical)/2.0)*z_near;
    egl_program = self.egl_program
    c = np.eye(4, dtype=np.float32)
    c[3,3] = 0
    c[3,2] = -1
    c[2,2] = -(z_near+z_far)/(z_far-z_near)
    c[2,3] = -2.0*(z_near*z_far)/(z_far-z_near)
    c[0,0] = 2.0*z_near/width
    c[1,1] = 2.0*z_near/height
    c = c.T
    
    projection_matrix_o = glGetUniformLocation(egl_program, 'uProjectionMatrix')
    projection_matrix = np.eye(4, dtype=np.float32)
    projection_matrix[...] = c
    projection_matrix = np.reshape(projection_matrix, (-1))
    glUniformMatrix4fv(projection_matrix_o, 1, GL_FALSE, projection_matrix) 

def output_shrink(K, L):
    """
    shrink the convolution output to half the size.
    used when both the annihilating filter and the uniform samples of sinusoids satisfy
    Hermitian symmetric.
    :param K: the annihilating filter size: K + 1
    :param L: length of the (complex-valued) b vector
    :return:
    """
    out_len = L - K
    if out_len % 2 == 0:
        half_out_len = np.int(out_len / 2.)
        mtx_r = np.hstack((np.eye(half_out_len),
                           np.zeros((half_out_len, half_out_len))))
        mtx_i = mtx_r
    else:
        half_out_len = np.int((out_len + 1) / 2.)
        mtx_r = np.hstack((np.eye(half_out_len),
                           np.zeros((half_out_len, half_out_len - 1))))
        mtx_i = np.hstack((np.eye(half_out_len - 1),
                           np.zeros((half_out_len - 1, half_out_len))))
    return linalg.block_diag(mtx_r, mtx_i) 

def classical_mds(self, D):
        ''' 
        Classical multidimensional scaling

        Parameters
        ----------
        D : square 2D ndarray
            Euclidean Distance Matrix (matrix containing squared distances between points
        '''

        # Apply MDS algorithm for denoising
        n = D.shape[0]
        J = np.eye(n) - np.ones((n,n))/float(n)
        G = -0.5*np.dot(J, np.dot(D, J))

        s, U = np.linalg.eig(G)

        # we need to sort the eigenvalues in decreasing order
        s = np.real(s)
        o = np.argsort(s)
        s = s[o[::-1]]
        U = U[:,o[::-1]]

        S = np.diag(s)[0:self.dim,:]
        self.X = np.dot(np.sqrt(S),U.T) 

def rotate_camera_to_point_at(up_from, lookat_from, up_to, lookat_to):
  inputs = [up_from, lookat_from, up_to, lookat_to]
  for i in range(4):
    inputs[i] = normalize(np.array(inputs[i]).reshape((-1,)))
  up_from, lookat_from, up_to, lookat_to = inputs
  r1 = r_between(lookat_from, lookat_to)

  new_x = np.dot(r1, np.array([1, 0, 0]).reshape((-1, 1))).reshape((-1))
  to_x = normalize(np.cross(lookat_to, up_to))
  angle = np.arccos(np.dot(new_x, to_x))
  if angle > ANGLE_EPS:
    if angle < np.pi - ANGLE_EPS:
      ax = normalize(np.cross(new_x, to_x))
      flip = np.dot(lookat_to, ax)
      if flip > 0:
        r2 = get_r_matrix(lookat_to, angle)
      elif flip < 0:
        r2 = get_r_matrix(lookat_to, -1. * angle)
    else:
      # Angle of rotation is too close to 180 degrees, direction of rotation
      # does not matter.
      r2 = get_r_matrix(lookat_to, angle)
  else:
    r2 = np.eye(3)
  return np.dot(r2, r1) 

def get_training_data(num_samples):
    """Generates some training data."""

    # As (x, y) Cartesian coordinates.
    x = np.random.randint(0, 2, size=(num_samples, 2))

    y = x[:, 0] + 2 * x[:, 1]  # 2-digit binary to integer.
    y = np.cast['int32'](y)

    x = np.cast['float32'](x) * 1.6 - 0.8  # Scales to [-1, 1].
    x += np.random.uniform(-0.1, 0.1, size=x.shape)

    y_ohe = np.cast['float32'](np.eye(4)[y])
    y = np.cast['float32'](np.expand_dims(y, -1))

    return x, y, y_ohe 

def get_mnist_data(binarize=False):
    """Puts the MNIST data in the right format."""

    (X_train, y_train), (X_test, y_test) = mnist.load_data()

    if binarize:
        X_test = np.where(X_test >= 10, 1, -1)
        X_train = np.where(X_train >= 10, 1, -1)
    else:
        X_train = (X_train.astype(np.float32) - 127.5) / 127.5
        X_test = (X_test.astype(np.float32) - 127.5) / 127.5

    X_train = np.expand_dims(X_train, axis=-1)
    X_test = np.expand_dims(X_test, axis=-1)

    y_train = np.eye(10)[y_train]
    y_test = np.eye(10)[y_test]

    return (X_train, y_train), (X_test, y_test) 

def test_circuit_generation_and_accuracy():
    for dim in range(2, 10):
        qubits = cirq.LineQubit.range(dim)
        u_generator = numpy.random.random(
            (dim, dim)) + 1j * numpy.random.random((dim, dim))
        u_generator = u_generator - numpy.conj(u_generator).T
        assert numpy.allclose(-1 * u_generator, numpy.conj(u_generator).T)

        unitary = scipy.linalg.expm(u_generator)
        circuit = cirq.Circuit()
        circuit.append(optimal_givens_decomposition(qubits, unitary))

        fermion_generator = QubitOperator(()) * 0.0
        for i, j in product(range(dim), repeat=2):
            fermion_generator += jordan_wigner(
                FermionOperator(((i, 1), (j, 0)), u_generator[i, j]))

        true_unitary = scipy.linalg.expm(
            get_sparse_operator(fermion_generator).toarray())
        assert numpy.allclose(true_unitary.conj().T.dot(true_unitary),
                              numpy.eye(2 ** dim, dtype=complex))

        test_unitary = cirq.unitary(circuit)
        assert numpy.isclose(
            abs(numpy.trace(true_unitary.conj().T.dot(test_unitary))), 2 ** dim) 

def run(self):
        D = np.eye(self.variables_num,self.variables_num)
        x = self.init_variables.reshape(-1,1)
        iteration = 0
        for i in range(self.epochs):
            g1 = gradients(self.func,x.flatten()).reshape((self.variables_num,1))
            d = D.dot(g1) # D is approximately equals to Hessain Matrix
            best_step = self.find_best_step(x,d)
            s = best_step * d
            x -= s
            iteration += 1
            g2 = gradients(self.func,x.flatten()).reshape((self.variables_num,1))
            if np.sqrt((g2 ** 2).sum()) < self.eps:
                break
            y = g2 - g1
            # update D
            D += s.dot(s.T) / ((s.T).dot(y)) - D.dot(y).dot(y.T).dot(D) / ((y.T).dot(D).dot(y))
            self.generations_points.append(x.flatten())
            self.generations_targets.append(self.func(x.flatten()))


        if self.func_type == newton_config.func_type_min:
            self.global_best_target = np.min(np.array(self.generations_targets))
            self.global_best_index = np.argmin(np.array(self.generations_targets))
            self.global_best_point = self.generations_points[int(self.global_best_index)]
        else:
            self.global_best_target = np.max(np.array(self.generations_targets))
            self.global_best_index = np.argmax(np.array(self.generations_targets))
            self.global_best_point = self.generations_points[int(self.global_best_index)] 

def __init__(self):
        ndim, dt = 4, 1.

        # Create Kalman filter model matrices.
        self._motion_mat = np.eye(2 * ndim, 2 * ndim)
        for i in range(ndim):
            self._motion_mat[i, ndim + i] = dt
        self._update_mat = np.eye(ndim, 2 * ndim)

        # Motion and observation uncertainty are chosen relative to the current
        # state estimate. These weights control the amount of uncertainty in
        # the model. This is a bit hacky.
        self._std_weight_position = 1. / 20
        self._std_weight_velocity = 1. / 160 

def calcSTM(self,dt,j):
        #allocate
        u = 0
        deltaU = 0
        t = 0
        counter = 0
        
        #For elliptic orbits, calculate period effects
        if self.beta[j] >0:
            P = 2*np.pi*self.mu[j]*self.beta[j]**(-3./2.)
            n = np.floor((dt + P/2 - 2*self.nu0[j]/self.beta[j])/P)
            deltaU = 2*np.pi*n*self.beta[j]**(-5./2.)
        
        #loop until convergence of the time array to the time step
        while (np.max(np.abs(t-dt)) > self.epsmult*np.spacing(dt)) and (counter < 1000):
            q = self.beta[j]*u**2./(1+self.beta[j]*u**2.)
            U0w2 = 1. - 2.*q
            U1w2 = 2.*(1.-q)*u
            temp = self.contFrac(q)
            U = 16./15.*U1w2**5.*temp + deltaU
            U0 = 2.*U0w2**2.-1.
            U1 = 2.*U0w2*U1w2
            U2 = 2.*U1w2**2.
            U3 = self.beta[j]*U + U1*U2/3.
            r = self.r0norm[j]*U0 + self.nu0[j]*U1 + self.mu[j]*U2
            t = self.r0norm[j]*U1 + self.nu0[j]*U2 + self.mu[j]*U3
            u = u - (t-dt)/(4.*(1.-q)*r)
            counter += 1
        
        if (counter == 1000):
            raise ValueError('Failed to converge on t: %e/%e'%(np.max(np.abs(t-dt)), self.epsmult*np.spacing(dt)))
        
        #Kepler solution
        f = 1 - self.mu[j]/self.r0norm[j]*U2
        g = self.r0norm[j]*U1 + self.nu0[j]*U2
        F = -self.mu[j]*U1/r/self.r0norm[j]
        G = 1 - self.mu[j]/r*U2
        
        Phi = np.vstack((np.hstack((np.eye(3)*f, np.eye(3)*g)),np.hstack((np.eye(3)*F, np.eye(3)*G))))
       
        return Phi 

def compute(self, max_iter):
        ftol = 1.0e-7
        dim_k = self.p0.shape[1]
        g = np.eye(dim_k + 1, dtype=self.p0.dtype)
        p = np.copy(self.p1)

        self.g_series = np.zeros((max_iter + 1, dim_k + 1, dim_k + 1), dtype=g.dtype)
        self.g_series[0, :, :] = g

        itr = -1
        for itr in range(max_iter):
            neighbor_idx = self.nearest.query(p)[1]
            targets = self.p0[neighbor_idx]
            R, t = _icp_find_rigid_transform(p, targets)

            new_p = np.dot(R, p.T).T + t
            if np.sum(np.abs(p - new_p)) < ftol:
                break

            p = np.copy(new_p)
            dg = _icp_Rt_to_matrix(R, t)
            new_g = np.dot(dg, g)
            g = np.copy(new_g)
            self.g_series[itr + 1, :, :] = g

        self.g_series[(itr+1):, :, :] = g

        return g, p, (itr + 1) 

def do(self):

        X = []
        scalings = []

        for k, p in enumerate(self.partitions):

            if p > 1:
                val = np.linspace(0, 1, p + 1)[1:-1]
                _X = []
                for i in range(self.n_dim):
                    for j in range(i + 1, self.n_dim):
                        x = np.zeros((len(val), self.n_dim))
                        x[:, i] = val
                        x[:, j] = 1 - val
                        _X.append(x)

                X.append(np.row_stack(_X + [np.eye(self.n_dim)]))

            elif p == 1:
                X.append(np.eye(self.n_dim))
            else:
                X.append(np.full(self.n_dim, 1 / self.n_dim)[None, :])

            scalings.append(1 - k / len(self.partitions))

        scalings = np.array(scalings)
        X = self._solve(X, scalings)

        return X


# ---------------------------------------------------------------------------------------------------------
# Energy Functions
# --------------------------------------------------------------------------------------------------------- 

def run(self):
        B = np.eye(self.variables_num, self.variables_num)
        x = self.init_variables.reshape(-1, 1)
        iteration = 0
        for i in range(self.epochs):
            g1 = gradients(self.func, x.flatten()).reshape((self.variables_num, 1))
            d = linalg.inv(B).dot(g1)  # D is approximately equals to Hessain Matrix
            best_step = self.find_best_step(x, d)
            s = best_step * d
            x -= s
            iteration += 1
            g2 = gradients(self.func, x.flatten()).reshape((self.variables_num, 1))
            if np.sqrt((g2 ** 2).sum()) < self.eps:
                break
            y = g2 - g1
            # update D
            B += y.dot(y.T) / ((y.T).dot(s)) - (B.dot(s).dot(s.T).dot(B)) / (s.T.dot(B).dot(s))
            self.generations_points.append(x.flatten())
            self.generations_targets.append(self.func(x.flatten()))

        if self.func_type == newton_config.func_type_min:
            self.global_best_target = np.min(np.array(self.generations_targets))
            self.global_best_index = np.argmin(np.array(self.generations_targets))
            self.global_best_point = self.generations_points[int(self.global_best_index)]
        else:
            self.global_best_target = np.max(np.array(self.generations_targets))
            self.global_best_index = np.argmax(np.array(self.generations_targets))
            self.global_best_point = self.generations_points[int(self.global_best_index)] 

def convert_to_onehot(sca_label, class_num=31):
    return np.eye(class_num)[sca_label] 

def CORAL_map(Xs,Xt):
    Ds = Xs.copy()
    Dt = Xt.copy()
      
    cov_src = np.ma.cov(Ds.T) + np.eye(Ds.shape[1])
    cov_tar = np.ma.cov(Dt.T) + np.eye(Dt.shape[1])
    
    Cs = scipy.linalg.sqrtm(np.linalg.inv(np.array(cov_src)))
    Ct = scipy.linalg.sqrtm(np.array(cov_tar))
    A_coral = np.dot(Cs, Ct)
    
    Xs_new = np.dot(Ds, A_coral)
        
    return Xs_new 

def test_random_interaction_operator_term(order, real, seed):
    op = ofctr.random_interaction_operator_term(order, real, seed)

    assert openfermion.is_hermitian(op)

    assert op.constant == 0
    assert op.one_body_tensor.shape == (order,) * 2
    assert op.two_body_tensor.shape == (order,) * 4

    for tensor in (op.one_body_tensor, op.two_body_tensor):
        for indices in np.argwhere(tensor):
            assert len(set(indices)) == order

    op_2 = ofctr.random_interaction_operator_term(order, real, seed)
    assert op == op_2

    if order == 1:
        assert op.one_body_tensor != 0
        assert op.two_body_tensor != 0
    elif order == 2:
        assert np.all((op.one_body_tensor == 0) == np.eye(2))
    elif order == 3:
        assert np.all(op.one_body_tensor == 0)
    elif order == 4:
        assert np.all(op.one_body_tensor == 0)
    else:
        assert np.all(op.one_body_tensor == 0)
        assert np.all(op.two_body_tensor == 0) 

def test_zero_weights():
    for gate_type in [
            ofc.QuadraticFermionicSimulationGate,
            ofc.CubicFermionicSimulationGate, ofc.QuarticFermionicSimulationGate
    ]:
        weights = (0,) * gate_type.num_weights()
        gate = gate_type(weights)
        n_qubits = gate.num_qubits()

        assert np.allclose(cirq.unitary(gate), np.eye(2**n_qubits))
        cirq.testing.assert_decompose_is_consistent_with_unitary(gate)

        operator = openfermion.InteractionOperator.zero(n_qubits)
        assert gate_type.from_interaction_operator(operator=operator) is None 

def assert_fswap_consistent(gate):
    gate = gate.__copy__()
    n_qubits = gate.num_qubits()
    for i in range(n_qubits - 1):
        fswap = cirq.kron(np.eye(1 << i), cirq.unitary(ofc.FSWAP),
                          np.eye(1 << (n_qubits - i - 2)))
        assert fswap.shape == (1 << n_qubits,) * 2
        generator = gate.qubit_generator_matrix
        fswapped_generator = np.linalg.multi_dot([fswap, generator, fswap])
        gate.fswap(i)
        assert np.allclose(gate.qubit_generator_matrix, fswapped_generator)
    for i in (-1, n_qubits):
        with pytest.raises(ValueError):
            gate.fswap(i)