Python numpy.vdot() Examples

def solve(self, f, tol=0):
        dx = -self.alpha*f

        n = len(self.dx)
        if n == 0:
            return dx

        df_f = np.empty(n, dtype=f.dtype)
        for k in xrange(n):
            df_f[k] = vdot(self.df[k], f)

        try:
            gamma = solve(self.a, df_f)
        except LinAlgError:
            # singular; reset the Jacobian approximation
            del self.dx[:]
            del self.df[:]
            return dx

        for m in xrange(n):
            dx += gamma[m]*(self.dx[m] + self.alpha*self.df[m])
        return dx 

def trace(self, **kwargs):
        r"""Trace of the density operator corresponding to the state.

        For pure states the trace corresponds to the squared norm of the ket vector.

        For physical states this should always be 1, any deviations from this value are due
        to numerical errors and Hilbert space truncation artefacts.

        Returns:
          float: trace of the state
        """
        # pylint: disable=unused-argument
        if self.is_pure:
            return np.vdot(self.ket(), self.ket()).real  # <s|s>

        # need some extra steps to trace over multimode matrices
        eqn_indices = [
            [indices[idx]] * 2 for idx in range(self._modes)
        ]  # doubled indices [['i','i'],['j','j'], ... ]
        eqn = "".join(
            chain.from_iterable(eqn_indices)
        )  # flatten indices into a single string 'iijj...'
        return np.einsum(eqn, self.dm()).real 

def testVdotExecution(self):
        a_data = np.array([1 + 2j, 3 + 4j])
        b_data = np.array([5 + 6j, 7 + 8j])
        a = tensor(a_data, chunk_size=1)
        b = tensor(b_data, chunk_size=1)

        t = vdot(a, b)

        res = self.executor.execute_tensor(t)[0]
        expected = np.vdot(a_data, b_data)
        np.testing.assert_equal(res, expected)

        a_data = np.array([[1, 4], [5, 6]])
        b_data = np.array([[4, 1], [2, 2]])
        a = tensor(a_data, chunk_size=1)
        b = tensor(b_data, chunk_size=1)

        t = vdot(a, b)

        res = self.executor.execute_tensor(t)[0]
        expected = np.vdot(a_data, b_data)
        np.testing.assert_equal(res, expected) 

def solve(self, f, tol=0):
        dx = -self.alpha*f

        n = len(self.dx)
        if n == 0:
            return dx

        df_f = np.empty(n, dtype=f.dtype)
        for k in xrange(n):
            df_f[k] = vdot(self.df[k], f)

        try:
            gamma = solve(self.a, df_f)
        except LinAlgError:
            # singular; reset the Jacobian approximation
            del self.dx[:]
            del self.df[:]
            return dx

        for m in xrange(n):
            dx += gamma[m]*(self.dx[m] + self.alpha*self.df[m])
        return dx 

def matvec(self, f):
        dx = -f/self.alpha

        n = len(self.dx)
        if n == 0:
            return dx

        df_f = np.empty(n, dtype=f.dtype)
        for k in xrange(n):
            df_f[k] = vdot(self.df[k], f)

        b = np.empty((n, n), dtype=f.dtype)
        for i in xrange(n):
            for j in xrange(n):
                b[i,j] = vdot(self.df[i], self.dx[j])
                if i == j and self.w0 != 0:
                    b[i,j] -= vdot(self.df[i], self.df[i])*self.w0**2*self.alpha
        gamma = solve(b, df_f)

        for m in xrange(n):
            dx += gamma[m]*(self.df[m] + self.dx[m]/self.alpha)
        return dx 

def solve(self, f, tol=0):
        dx = -self.alpha*f

        n = len(self.dx)
        if n == 0:
            return dx

        df_f = np.empty(n, dtype=f.dtype)
        for k in xrange(n):
            df_f[k] = vdot(self.df[k], f)

        try:
            gamma = solve(self.a, df_f)
        except LinAlgError:
            # singular; reset the Jacobian approximation
            del self.dx[:]
            del self.df[:]
            return dx

        for m in xrange(n):
            dx += gamma[m]*(self.dx[m] + self.alpha*self.df[m])
        return dx 

def matvec(self, f):
        dx = -f/self.alpha

        n = len(self.dx)
        if n == 0:
            return dx

        df_f = np.empty(n, dtype=f.dtype)
        for k in xrange(n):
            df_f[k] = vdot(self.df[k], f)

        b = np.empty((n, n), dtype=f.dtype)
        for i in xrange(n):
            for j in xrange(n):
                b[i,j] = vdot(self.df[i], self.dx[j])
                if i == j and self.w0 != 0:
                    b[i,j] -= vdot(self.df[i], self.df[i])*self.w0**2*self.alpha
        gamma = solve(b, df_f)

        for m in xrange(n):
            dx += gamma[m]*(self.df[m] + self.dx[m]/self.alpha)
        return dx 

def test_sandwich(nr_sites, local_dim, rank, rgen, dtype):
    mps = factory.random_mpa(nr_sites, local_dim, rank,
                             randstate=rgen, dtype=dtype, normalized=True)
    mps2 = factory.random_mpa(nr_sites, local_dim, rank,
                              randstate=rgen, dtype=dtype, normalized=True)
    mpo = factory.random_mpa(nr_sites, [local_dim] * 2, rank,
                             randstate=rgen, dtype=dtype)
    mpo.canonicalize()
    mpo /= mp.trace(mpo)

    vec = mps.to_array().ravel()
    op = mpo.to_array_global().reshape([local_dim**nr_sites] * 2)
    res_arr = np.vdot(vec, np.dot(op, vec))
    res_mpo = mp.inner(mps, mp.dot(mpo, mps))
    res_sandwich = mp.sandwich(mpo, mps)
    assert_almost_equal(res_mpo, res_arr)
    assert_almost_equal(res_sandwich, res_arr)

    vec2 = mps2.to_array().ravel()
    res_arr = np.vdot(vec2, np.dot(op, vec))
    res_mpo = mp.inner(mps2, mp.dot(mpo, mps))
    res_sandwich = mp.sandwich(mpo, mps, mps2)
    assert_almost_equal(res_mpo, res_arr)
    assert_almost_equal(res_sandwich, res_arr) 

def vecangle(v1, v2):
    """angle between two vectors v1, v2 in radians"""
    v0 = pt(0, 0)
        
    if np.allclose(v1, v0) or np.allclose(v2, v0):
        return np.nan
    
    if np.allclose(v1, v2):
        return 0
    
    num = np.vdot(v1, v2)
    denom = (np.linalg.norm(v1) * np.linalg.norm(v2))
    
    if np.isclose(num, denom):
        return 0
    
    return math.acos(num / denom) 

def symm_inner_product_array(
    num_states, num_vecs, max_vecs_per_node, verbosity=1):
    """
    Computes symmetric inner product array from known vecs (as in POD).
    """
    vec_handles = [mr.VecHandlePickle(join(data_dir, row_vec_name%row_num))
        for row_num in mr.range(num_vecs)]

    generate_vecs(data_dir, num_states, vec_handles)

    my_VS = mr.VectorSpaceHandles(
        inner_product=np.vdot, max_vecs_per_node=max_vecs_per_node,
        verbosity=verbosity)

    prof = cProfile.Profile()
    start_time = time.time()
    prof.runcall(my_VS.compute_symm_inner_product_array, vec_handles)
    total_time = time.time() - start_time
    prof.dump_stats('IP_symm_array_r%d.prof'%mr.parallel.get_rank())

    return total_time 

def matvec(self, f):
        dx = -f/self.alpha

        n = len(self.dx)
        if n == 0:
            return dx

        df_f = np.empty(n, dtype=f.dtype)
        for k in xrange(n):
            df_f[k] = vdot(self.df[k], f)

        b = np.empty((n, n), dtype=f.dtype)
        for i in xrange(n):
            for j in xrange(n):
                b[i,j] = vdot(self.df[i], self.dx[j])
                if i == j and self.w0 != 0:
                    b[i,j] -= vdot(self.df[i], self.df[i])*self.w0**2*self.alpha
        gamma = solve(b, df_f)

        for m in xrange(n):
            dx += gamma[m]*(self.df[m] + self.dx[m]/self.alpha)
        return dx 

def setUp(self):
        if not os.access('.', os.W_OK):
            raise RuntimeError('Cannot write to current directory')
        self.test_dir = 'files_vectorspace_DELETE_ME'
        if not os.path.isdir(self.test_dir):
            parallel.call_from_rank_zero(os.mkdir, self.test_dir)

        self.max_vecs_per_proc = 10
        self.total_num_vecs_in_mem = (
            parallel.get_num_procs() * self.max_vecs_per_proc)

        self.vec_space = vspc.VectorSpaceHandles(
            inner_product=np.vdot, verbosity=0)
        self.vec_space.max_vecs_per_proc = self.max_vecs_per_proc

        # Default data members; set verbosity to 0 even though default is 1
        # so messages won't print during tests
        self.default_data_members = {
            'inner_product': np.vdot, 'max_vecs_per_node': 10000,
            'max_vecs_per_proc': (
                10000 * parallel.get_num_nodes() // parallel.get_num_procs()),
            'verbosity': 0, 'print_interval': 10, 'prev_print_time': 0.}
        parallel.barrier() 

def checkDeflate(self):
        if self.currentSize == self.maxM-1:
            self.deflated = 1
            for i in range(self.nEigen):
                self.sol[:self.currentSize] = self.evecs[:self.currentSize,i]
                self.constructSolV()            # Finds the "best" eigenvector for this eigenvalue
                self.Avlist[i,:] = self.cv.copy() # Puts this guess in self.Avlist rather than self.vlist for now...
                                                # since this would mess up self.constructSolV()'s solution
            for i in range(self.nEigen):
                self.cv = self.Avlist[i,:].copy() # This is actually the "best" eigenvector v, not A*v (see above)
                self.gramSchmidtCurrentVec(i)
                self.vlist[i,:] = self.cv.copy()

            orthog = 0.0
            for j in range(self.nEigen):
                for i in range(j):
                    orthog += np.vdot(self.vlist[j,:],self.vlist[i,:])**2.0

            for i in range(self.nEigen):
                self.cv = self.vlist[i].copy() # This is actually the "best" eigenvector v, not A*v (see above)
                self.hMult()                   # Use current vector cv to create cAv
                self.Avlist[i] = self.cAv.copy() 

def gramSchmidtCurrentVec(self,northo):
        for k in range(northo):
            fac = np.vdot( self.vlist[k,:], self.cv )
            self.cv -= fac * self.vlist[k,:] #/ np.vdot(self.vlist[i],self.vlist[i])
        cvnorm = np.linalg.norm(self.cv)
        if cvnorm < 1e-4:
            self.cv = np.random.rand(self.size)
            for k in range(northo):
                fac = np.vdot( self.vlist[k,:], self.cv )
                self.cv -= fac * self.vlist[k,:] #/ np.vdot(self.vlist[i],self.vlist[i])
            ########################################################################
            #  Sometimes some of the created vectors are linearly dependent.  i
            #  To get around this I'm not sure what to do other than to throw that
            #  solution out and start at that eigenvalue
            ########################################################################
            print(" ERROR!!!! ... restarting at eigenvalue #%" % \
                        (northo, cvnorm))
            self.ciEig = northo
        self.cv /= np.linalg.norm(self.cv) 

def evaluate_with_statevector(self, quantum_state):
        """

        Args:
            quantum_state (numpy.ndarray): quantum state

        Returns:
            float: the mean value
            float: the standard deviation
        Raises:
            AquaError: if Operator is empty
        """
        if self.is_empty():
            raise AquaError("Operator is empty, check the operator.")

        avg = np.vdot(quantum_state, self._matrix.dot(quantum_state))
        return avg, 0.0 

def guessInitial(self):
        nrm = np.linalg.norm(self.x0)
        self.x0 *= 1./nrm
        self.currentSize = self.nEigen
        for i in range(self.currentSize):
            self.vlist[i] *= 0.0
            self.vlist[i,:] += np.random.rand(self.size)
            self.vlist[i,i] /= np.linalg.norm(self.vlist[i,i])
            self.vlist[i,i] *= 12.
            for j in range(i):
                fac = np.vdot( self.vlist[j,:], self.vlist[i,:] )
                self.vlist[i,:] -= fac * self.vlist[j,:]
            self.vlist[i,:] /= np.linalg.norm(self.vlist[i,:])
        for i in range(self.currentSize):
            self.cv = self.vlist[i].copy()
            self.hMult()
            self.Avlist[i] = self.cAv.copy()
        self.constructSubspace() 

def evaluate_with_statevector(self, quantum_state):
        """
        Args:
            quantum_state (numpy.ndarray): a quantum state.

        Returns:
            float: the mean value
            float: the standard deviation
        Raises:
            AquaError: if Operator is empty
        """
        if self.is_empty():
            raise AquaError("Operator is empty, check the operator.")
        # convert to matrix first?
        # pylint: disable=import-outside-toplevel
        from .op_converter import to_matrix_operator
        mat_op = to_matrix_operator(self)
        avg = np.vdot(quantum_state, mat_op._matrix.dot(quantum_state))
        return avg, 0.0

    # pylint: disable=arguments-differ 

def test_adjoint_basis_vec_optional(self):
        """Test that adjoint modes default to direct modes"""
        no_adjoints_LTI_proj = lgp.LTIGalerkinProjectionHandles(
            np.vdot, self.basis_vec_handles, is_basis_orthonormal=True,
            verbosity=0)
        np.testing.assert_equal(
            no_adjoints_LTI_proj.adjoint_basis_vec_handles,
            self.basis_vec_handles) 

def get_AvgSim(alchemyAvgSoapA, alchemyAvgSoapB, chemicalSimGen):
    chemicalSim = chemicalSimGen(alchemyAvgSoapA, alchemyAvgSoapB)
    AvgSim = 0
    for spA in alchemyAvgSoapA:
        for spB in alchemyAvgSoapB:
            theta = chemicalSim[spA + spB]
            if theta == 0.:
                continue
            else:
                AvgSim += theta * np.vdot(alchemyAvgSoapA[spA], alchemyAvgSoapB[spB])

    return AvgSim 

def test_compute_decomp(self):
        """Test computation of the correlation array and SVD arrays."""
        rtol = 1e-10
        atol = 1e-12

        # Compute POD using modred
        POD = pod.PODHandles(inner_product=np.vdot, verbosity=0)
        eigvals, eigvecs = POD.compute_decomp(self.vec_handles)

        # Test correlation array values by simply recomputing them.  Here simply
        # take all inner products, rather than assuming a symmetric inner
        # product.
        np.testing.assert_allclose(
            POD.correlation_array,
            POD.vec_space.compute_inner_product_array(
                self.vec_handles, self.vec_handles),
            rtol=rtol, atol=atol)

        # Check POD eigenvectors and eigenvalues
        np.testing.assert_allclose(
            self.vecs_array.conj().T.dot(self.vecs_array.dot(eigvecs)),
            eigvecs.dot(np.diag(eigvals)), rtol=rtol, atol=atol)

        # Check that returned values match internal values
        np.testing.assert_equal(eigvals, POD.eigvals)
        np.testing.assert_equal(eigvecs, POD.eigvecs)


    #@unittest.skip('Testing something else.') 

def test_compute_modes(self):
        rtol = 1e-10
        atol = 1e-12

        # Compute POD using modred.  (The properties defining a POD mode require
        # manipulations involving the correct decomposition, so we cannot
        # isolate the mode computation from the decomposition step.)
        POD = pod.PODHandles(inner_product=np.vdot, verbosity=0)
        POD.compute_decomp(self.vec_handles)

        # Select a subset of modes to compute.  Compute at least half
        # the modes, and up to all of them.  Make sure to use unique
        # values.  (This may reduce the number of modes computed.)
        num_modes = parallel.call_and_bcast(
            np.random.randint,
            POD.eigvals.size // 2, POD.eigvals.size + 1)
        mode_idxs = np.unique(parallel.call_and_bcast(
            np.random.randint,
            0, POD.eigvals.size, num_modes))

        # Create handles for the modes
        mode_handles = [VecHandlePickle(self.mode_path % i) for i in mode_idxs]

        # Compute modes
        POD.compute_modes(mode_idxs, mode_handles, vec_handles=self.vec_handles)

        # Test modes
        np.testing.assert_allclose(
            POD.vec_space.compute_inner_product_array(
                mode_handles, self.vec_handles).dot(
                    POD.vec_space.compute_inner_product_array(
                        self.vec_handles, mode_handles)),
            np.diag(POD.eigvals[mode_idxs]),
            rtol=rtol, atol=atol)


    #@unittest.skip('Testing something else.') 

def get_AvgDeltaSim(AvgSoapA, AvgSoapB, chem_channels=False):
    if chem_channels:
        AvgDeltaSim = 0
        for spA in AvgSoapA:
            for spB in AvgSoapB:
                if np.all(spA != spB):
                    continue
                elif np.all(spA == spB):
                    AvgDeltaSim += np.vdot(AvgSoapA[spA], AvgSoapB[spB])
    else:
        AvgDeltaSim = np.vdot(AvgSoapA, AvgSoapB)

    return AvgDeltaSim 

def test_refcount_vdot(self):
        # Changeset #3443
        _assert_valid_refcount(np.vdot) 

def test_reduce_B(self):
        """Given modes, test reduced B array, orthogonal and non-orthogonal."""
        B_returned = self.LTI_proj.reduce_B(self.B_on_standard_basis_handles)
        np.testing.assert_allclose(B_returned, self.B_true)

        LTI_proj = lgp.LTIGalerkinProjectionHandles(
            np.vdot, self.basis_vec_handles,
            adjoint_basis_vec_handles=self.adjoint_basis_vec_handles,
            is_basis_orthonormal=False, verbosity=0)
        B_returned = LTI_proj.reduce_B(self.B_on_standard_basis_handles)
        np.testing.assert_allclose(B_returned, self.B_true_non_orth)


    #@unittest.skip('Testing something else') 

def setUp(self):
        if not os.access('.', os.W_OK):
            raise RuntimeError('Cannot write to current directory')

        self.test_dir ='file_LTIGalerkinProj_DELETE_ME'
        if parallel.is_rank_zero() and not os.path.exists(self.test_dir):
            os.mkdir(self.test_dir)
        parallel.barrier()

        self.basis_vec_path = join(self.test_dir, 'basis_vec_%02d.txt')
        self.adjoint_basis_vec_path = join(
            self.test_dir, 'adjoint_basis_vec_%02d.txt')
        self.A_on_basis_vec_path = join(self.test_dir, 'A_on_mode_%02d.txt')
        self.B_on_basis_path = join(self.test_dir, 'B_on_basis_%02d.txt')
        self.C_on_basis_vec_path = join(self.test_dir, 'C_on_mode_%02d.txt')

        self.num_basis_vecs = 10
        self.num_adjoint_basis_vecs = 10
        self.num_states = 11
        self.num_inputs = 3
        self.num_outputs = 2

        self.generate_data_set(
            self.num_basis_vecs, self.num_adjoint_basis_vecs,
            self.num_states, self.num_inputs, self.num_outputs)

        self.LTI_proj = lgp.LTIGalerkinProjectionHandles(
            np.vdot, self.basis_vec_handles,
            adjoint_basis_vec_handles=self.adjoint_basis_vec_handles,
            is_basis_orthonormal=True, verbosity=0) 

def lin_combine(
    num_states, num_bases, num_products, max_vecs_per_node, verbosity=1):
    """
    Computes linear combination of vecs from saved vecs and random coeffs

    num_bases is number of vecs to be linearly combined
    num_products is the resulting number of vecs
    """

    basis_handles = [mr.VecHandlePickle(join(data_dir, basis_name%basis_num))
        for basis_num in mr.range(num_bases)]
    product_handles = [mr.VecHandlePickle(join(data_dir,
        product_name%product_num))
        for product_num in mr.range(num_products)]

    generate_vecs(data_dir, num_states, basis_handles)
    my_VS = mr.VectorSpaceHandles(
        inner_product=np.vdot, max_vecs_per_node=max_vecs_per_node,
        verbosity=verbosity)
    coeff_array = np.random.random((num_bases, num_products))
    mr.parallel.barrier()

    prof = cProfile.Profile()
    start_time = time.time()
    prof.runcall(my_VS.lin_combine, *(product_handles, basis_handles,
        coeff_array))
    total_time = time.time() - start_time
    prof.dump_stats('lincomb_r%d.prof'%mr.parallel.get_rank())
    return total_time 

def inner_product_array(
    num_states, num_rows, num_cols, max_vecs_per_node, verbosity=1):
    """
    Computes inner products from known vecs.

    Remember that rows correspond to adjoint modes and cols to direct modes
    """
    col_vec_handles = [mr.VecHandlePickle(join(data_dir, col_vec_name%col_num))
        for col_num in mr.range(num_cols)]
    row_vec_handles = [mr.VecHandlePickle(join(data_dir, row_vec_name%row_num))
        for row_num in mr.range(num_rows)]

    generate_vecs(data_dir, num_states, row_vec_handles+col_vec_handles)

    my_VS = mr.VectorSpaceHandles(
        inner_product=np.vdot, max_vecs_per_node=max_vecs_per_node,
        verbosity=verbosity)

    prof = cProfile.Profile()
    start_time = time.time()
    prof.runcall(
        my_VS.compute_inner_product_array, *(col_vec_handles, row_vec_handles))
    total_time = time.time() - start_time
    prof.dump_stats('IP_array_r%d.prof'%mr.parallel.get_rank())

    return total_time 

def inner_product_array_uniform(vec1, vec2):
    """Takes inner product of numpy arrays without weighting. The first element
    is conjugated, i.e., IP = np.dot(vec1.conj().T, v2)
    """
    return np.vdot(vec1, vec2) 

def _calculate(self):
        min_, max_ = self._bounds
        n = self.reference.size
        f = n * (max_ - min_)**2

        diff = self.other - self.reference
        value = np.vdot(diff, diff) / f

        # calculate the gradient only when needed
        self._g_diff = diff
        self._g_f = f
        gradient = None
        return value, gradient 

def _analytic_prop_1q_gates(self, total_samples, omega_0, omega_a, omega_d0, phi):
        """Compute proportion for 0 and 1 states analytically for single qubit gates.
        Args:
            total_samples (int): length of pulses
            omega_0 (float): Q0 freq
            omega_a (float): Q0 drive amplitude
            omega_d0 (flaot): Q0 drive frequency
            phi (float): drive phase
        Returns:
            exp_prop (dict): expected value of 0 and 1 proportions from analytic computation
            """
        time = total_samples
        # write Hrot analytically
        h_rot = np.array([[
            (omega_d0 - omega_0) / 2,
            np.exp(1j * phi) * omega_a / 2
        ], [np.exp(-1j * phi) * omega_a / 2, -(omega_d0 - omega_0) / 2]])
        # exponentiate
        u_rot = expm(-1j * h_rot * time)
        state0 = np.array([1, 0])

        # compute analytic prob (proportion) of 0 state
        mat_elem0 = np.vdot(state0, np.dot(u_rot, state0))
        prop0 = np.abs(mat_elem0)**2

        # return expected proportion
        exp_prop = {'0': prop0, '1': 1 - prop0}
        return exp_prop