Fermion annihilation operator in matrix representation for a indexed
       particle in a bounded N-particles fermion fock space
975057
def destruct(particles, index):
    """"""

    mat = np.zeros((2**particles, 2**particles))

    flipper = 2**index
    for i in range(2**particles):
        ispin = btest(i, index)
        if ispin == 1:
            mat[i ^ flipper, i] = phase(i, index)
    return csr_matrix(mat)
17374
def quasiparticle_weight(self):
        """"""
        return np.array([self.expected(op)**2 for op in self.oper['O']])
713083
private void computeA( DenseMatrix64F A , DenseMatrix64F H , double lambda )
    {
        final int numParam = param.getNumElements();

        A.set(H);
        for( int i = 0; i < numParam; i++ ) {
            A.set(i,i, A.get(i,i) + lambda);
        }
    }
548313
def objective_fun(theta, hamiltonian=None,
                  quantum_resource=QVMConnection(sync_endpoint='http://localhost:5000')):
    
    if hamiltonian is None:
        # Hamiltonian is Identity
        return 1.0

    if isinstance(hamiltonian, PauliSum):
        result = estimate_locally_commuting_operator(ucc_circuit(theta), hamiltonian,
                                                     1.0E-6, quantum_resource=quantum_resource)
        result = result[0][0].real  # first output is expected value, second is variance, third is shots
    elif isinstance(hamiltonian, np.ndarray) and isinstance(quantum_resource, QVMConnection):
        wf = quantum_resource.wavefunction(ucc_circuit(theta))
        wf = wf.amplitudes.reshape((-1, 1))
        result = np.conj(wf).T.dot(hamiltonian).dot(wf)[0, 0].real
        print(result)
    else:
        raise TypeError("type of hamiltonian or qvm is unrecognized")

    return result
1145532
def multidimensional_discretization(rho, sigma, N=3, method='rouwenhorst', m=2):
    

    # rho is assumed to be a scalar
    # sigma is a positive symmetric matrix
    # N number of points in each non-degenerate dimension
    # m : standard deviations to approximate

    import scipy.linalg
    from itertools import product

    d = sigma.shape[1]

    sigma = sigma.copy()

    zero_columns = np.where(sigma.sum(axis=0)==0)[0]
    for i in zero_columns:
        sigma[i,i] = 1

    L = scipy.linalg.cholesky(sigma)

    N = int(N)

    if method=='tauchen':
        [nodes_1d, probas_1d] = tauchen(N, 0, rho, 1, m=m)
    elif method=='rouwenhorst':
        [nodes_1d, probas_1d] = rouwenhorst(rho, 1, N)

    markov_nodes = np.array( list( product( *([nodes_1d]*d))) ).T
    markov_indices = np.array( list( product( *([range(N)]*d)  )  ) ).T

    markov_nodes = np.dot(L, markov_nodes)


    transition_matrix = 1
    for i in range(d):
        transition_matrix = np.kron(transition_matrix, probas_1d)

    markov_nodes = np.ascontiguousarray(markov_nodes.T)

    for i in zero_columns:
        markov_nodes[:,i] = 0

    return [markov_nodes, transition_matrix]
711476
def bethe_fermi(energy, quasipart, shift, hopping, beta):
    """"""
    return fermi_dist(quasipart * energy - shift, beta) \
            * bethe_lattice(energy, hopping)
714096
private void computePrincipleVectors(int numPrincipleVectors) {

        // Group all of the sspace's vectors into a Matrix so that we can
        // cluster them
        final int numTerms = sspace.getWords().size();
        final List<DoubleVector> termVectors = 
            new ArrayList<DoubleVector>(numTerms);
        String[] termAssignments = new String[numTerms];
        int k = 0;
        for (String term : sspace.getWords()) {
            termVectors.add(Vectors.asDouble(sspace.getVector(term)));
            termAssignments[k++] = term;
        }

        Random random = new Random();

        final DoubleVector[] principles = new DoubleVector[numPrincipleVectors];
        for (int i = 0; i < principles.length; ++i)
            principles[i] = new DenseVector(sspace.getVectorLength());

        // Randomly assign the data set to different intitial clusters 
        final MultiMap<Integer,Integer> clusterAssignment = 
            new HashMultiMap<Integer,Integer>();
        for (int i = 0; i < numTerms; ++i) 
            clusterAssignment.put(random.nextInt(numPrincipleVectors), i);
        
        int numIters = 1;
        for (int iter = 0; iter < numIters; ++iter) {
            verbose(LOGGER, "Computing principle vectors (round %d/%d)", iter+1, numIters);
            // Compute the centroids
            for (Map.Entry<Integer,Set<Integer>> e
                     : clusterAssignment.asMap().entrySet()) {
                int cluster = e.getKey();
                DoubleVector principle = new DenseVector(sspace.getVectorLength());
                principles[cluster] = principle;
                for (Integer row : e.getValue())
                    VectorMath.add(principle, termVectors.get(row));
            }

            // Reassign each element to the centroid to which it is closest
            clusterAssignment.clear();
            final int numThreads = workQueue.availableThreads();
            Object key = workQueue.registerTaskGroup(numThreads);

            for (int threadId_ = 0; threadId_ < numThreads; ++threadId_) {
                final int threadId = threadId_;
                workQueue.add(key, new Runnable() {
                        public void run() {
                            // Thread local cache of all the cluster assignments
                            MultiMap<Integer,Integer> clusterAssignment_ = 
                                new HashMultiMap<Integer,Integer>();
                            // For each of the vectors that this thread is
                            // responsible for, find the principle vector to
                            // which it is closest
                            for (int i = threadId; i < numTerms; i += numThreads) {
                                DoubleVector v = termVectors.get(i);
                                double highestSim = -Double.MAX_VALUE;
                                int pVec = -1;
                                for (int j = 0; j < principles.length; ++j) {
                                    DoubleVector principle = principles[j];
                                    double sim = Similarity.cosineSimilarity(
                                        v, principle);
                                    assert sim >= -1 && sim <= 1 : "similarity "
                                        + " to principle vector " + j + " is "
                                        + "outside the expected range: " + sim;
                                    if (sim > highestSim) {
                                        highestSim = sim;
                                        pVec = j;
                                    }
                                }
                                assert pVec != -1 : "Could not find match to "
                                    + "any of the " + principles.length + 
                                    " principle vectors";
                                clusterAssignment_.put(pVec, i);
                            }
                            // Once all the vectors for this thread have been
                            // mapped, update the thread-shared mapping.  Do
                            // this here, rather than for each comparison to
                            // minimize the locking
                            synchronized (clusterAssignment) {
                                clusterAssignment.putAll(clusterAssignment_);
                            }
                        }
                    });
            }
            workQueue.await(key);
        }

        double mean = numTerms / (double)numPrincipleVectors;
        double variance = 0d;
        for (Map.Entry<Integer,Set<Integer>> e
                 : clusterAssignment.asMap().entrySet()) {
            Set<Integer> rows = e.getValue();
            Set<String> terms = new HashSet<String>();
            for (Integer i : rows)
                terms.add(termAssignments[i]);
             verbose(LOGGER, "Principle vectod %d is closest to %d terms",
                     e.getKey(), terms.size());
            double diff = mean - terms.size();
            variance += diff * diff;
            principleVectorToNearestTerms.putMany(
                principles[e.getKey()], terms);
        }

        verbose(LOGGER, "Average number terms per principle vector: %f, "+
                "(%f stddev)", mean, Math.sqrt(variance / numPrincipleVectors));
    }
1360873
def H(self) -> 'QubitVector':
        """"""
        N = self.qubit_nb
        R = self.rank

        # (super) operator transpose
        tensor = self.tensor
        tensor = bk.reshape(tensor, [2**(N*R//2)] * 2)
        tensor = bk.transpose(tensor)
        tensor = bk.reshape(tensor, [2] * R * N)

        tensor = bk.conj(tensor)

        return QubitVector(tensor, self.qubits)
966146
def _apply_adj(self, a):
        # is projection the zero operator?
        if self.V.shape[1] == 0:
            return numpy.zeros(a.shape)
        ''''''
        c = inner(self.V, a, ip_B=self.ip_B)
        if self.Q is not None and self.R is not None:
            c = self.Q.dot(scipy.linalg.solve_triangular(self.R.T.conj(), c,
                                                         lower=True))
        return self.W.dot(c)
305167
def adjoint(self):
        
        return Laplacian(self.range, self.domain,
                         pad_mode=self.pad_mode, pad_const=0)
436715
def bosonic_constraints(a):
    
    substitutions = {}
    for i, ai in enumerate(a):
        substitutions[ai * Dagger(ai)] = 1.0 + Dagger(ai) * ai
        for aj in a[i+1:]:
            # substitutions[ai*Dagger(aj)] = -Dagger(ai)*aj
            substitutions[ai*Dagger(aj)] = Dagger(aj)*ai
            substitutions[Dagger(ai)*aj] = aj*Dagger(ai)
            substitutions[ai*aj] = aj*ai
            substitutions[Dagger(ai) * Dagger(aj)] = Dagger(aj) * Dagger(ai)

    return substitutions
545209
private static void eltran(DenseMatrix A, DenseMatrix V, int[] perm) {
        int n = A.nrows();
        for (int mp = n - 2; mp > 0; mp--) {
            for (int k = mp + 1; k < n; k++) {
                V.set(k, mp, A.get(k, mp - 1));
            }
            int i = perm[mp];
            if (i != mp) {
                for (int j = mp; j < n; j++) {
                    V.set(mp, j, V.get(i, j));
                    V.set(i, j, 0.0);
                }
                V.set(i, mp, 1.0);
            }
        }
    }
400752
def spectral_embedding(self, n):
        
        coords = spectral_embedding(self._affinity, n)
        return CoordinateMatrix(normalise_rows(coords))
1244574
def locality_preserving_projections(self, coordinates, num_dims=None):
    ''''''
    X = np.atleast_2d(coordinates)  # n x d
    L = self.laplacian(normed=True)  # n x n
    u,s,_ = np.linalg.svd(X.T.dot(X))
    Fplus = np.linalg.pinv(u * np.sqrt(s))  # d x d
    n, d = X.shape
    if n >= d:  # optimized order: F(X'LX)F'
      T = Fplus.dot(X.T.dot(L.dot(X))).dot(Fplus.T)
    else:  # optimized order: (FX')L(XF')
      T = Fplus.dot(X.T).dot(L.dot(X.dot(Fplus.T)))
    L = 0.5*(T+T.T)
    return _null_space(L, num_vecs=num_dims, overwrite=True)
857861
@Override
    public DMatrixRMaj getT(DMatrixRMaj T ) {
        T = UtilDecompositons_DDRM.checkZeros(T,N,N);

        T.data[0] = QT.data[0];

        for( int i = 1; i < N; i++ ) {
            T.set(i,i, QT.get(i,i));
            double a = QT.get(i-1,i);
            T.set(i-1,i,a);
            T.set(i,i-1,a);
        }

        if( N > 1 ) {
            T.data[(N-1)*N+N-1] = QT.data[(N-1)*N+N-1];
            T.data[(N-1)*N+N-2] = QT.data[(N-2)*N+N-1];
        }
            
        return T;
    }
461801
public static ZMatrixRMaj hermitianPosDef(int width, Random rand) {
        // This is not formally proven to work.  It just seems to work.
        ZMatrixRMaj a = RandomMatrices_ZDRM.rectangle(width,1,rand);
        ZMatrixRMaj b = new ZMatrixRMaj(1,width);
        ZMatrixRMaj c = new ZMatrixRMaj(width,width);

        CommonOps_ZDRM.transposeConjugate(a,b);
        CommonOps_ZDRM.mult(a, b, c);

        for( int i = 0; i < width; i++ ) {
            c.data[2*(i*width+i)] += 1;
        }

        return c;
    }
540051
def trace(self) -> bk.BKTensor:
        
        N = self.qubit_nb
        R = self.rank
        if R == 1:
            raise ValueError('Cannot take trace of vector')

        tensor = bk.reshape(self.tensor, [2**(N*R//2)] * 2)
        tensor = bk.trace(tensor)
        return tensor