1#   Licensed under the Apache License, Version 2.0 (the "License");
2#   you may not use this file except in compliance with the License.
3#   You may obtain a copy of the License at
4#
5#       http://www.apache.org/licenses/LICENSE-2.0
6#
7#   Unless required by applicable law or agreed to in writing, software
8#   distributed under the License is distributed on an "AS IS" BASIS,
9#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
10#   See the License for the specific language governing permissions and
11#   limitations under the License.
12"""Tests for sparse_tools.py."""
13import os
14
15import unittest
16
17import itertools
18
19import numpy
20from numpy.linalg import multi_dot
21
22from scipy.linalg import eigh, norm
23from scipy.sparse import csc_matrix
24import scipy
25from scipy.special import comb
26
27from openfermion.config import DATA_DIRECTORY
28from openfermion.chem import MolecularData
29from openfermion.hamiltonians import (fermi_hubbard, jellium_model,
30                                      wigner_seitz_length_scale,
31                                      number_operator)
32from openfermion.ops.operators import (FermionOperator, BosonOperator,
33                                       QuadOperator, QubitOperator,
34                                       MajoranaOperator)
35from openfermion.ops.representations import DiagonalCoulombHamiltonian
36from openfermion.transforms.opconversions import (jordan_wigner,
37                                                  get_fermion_operator,
38                                                  normal_ordered)
39from openfermion.transforms.repconversions import (fourier_transform,
40                                                   get_interaction_operator)
41
42from openfermion.testing.testing_utils import random_hermitian_matrix
43from openfermion.linalg.sparse_tools import (
44    get_sparse_operator, get_ground_state, eigenspectrum, expectation,
45    jw_number_restrict_state, inner_product, jw_sz_restrict_state,
46    jw_get_ground_state_at_particle_number, jw_sparse_givens_rotation,
47    pauli_matrix_map, sparse_eigenspectrum, jordan_wigner_sparse,
48    kronecker_operators, identity_csc, pauli_x_csc, qubit_operator_sparse,
49    get_linear_qubit_operator_diagonal, jw_number_indices,
50    get_number_preserving_sparse_operator, jw_configuration_state,
51    jw_hartree_fock_state, jw_sz_restrict_operator, jw_sz_indices,
52    jw_number_restrict_operator, variance,
53    expectation_computational_basis_state,
54    expectation_db_operator_with_pw_basis_state, get_gap, boson_ladder_sparse,
55    boson_operator_sparse, single_quad_op_sparse, _iterate_basis_)
56
57from openfermion.utils.operator_utils import (is_hermitian, count_qubits,
58                                              hermitian_conjugated)
59from openfermion.utils.indexing import up_index, down_index
60from openfermion.utils.grid import Grid
61from openfermion.hamiltonians.jellium_hf_state import (
62    lowest_single_particle_energy_states)
63from openfermion.linalg.linear_qubit_operator import LinearQubitOperator
64
65
66class EigenSpectrumTest(unittest.TestCase):
67
68    def setUp(self):
69        self.n_qubits = 5
70        self.majorana_operator = MajoranaOperator((1, 4, 9))
71        self.fermion_term = FermionOperator('1^ 2^ 3 4', -3.17)
72        self.fermion_operator = self.fermion_term + hermitian_conjugated(
73            self.fermion_term)
74        self.qubit_operator = jordan_wigner(self.fermion_operator)
75        self.interaction_operator = get_interaction_operator(
76            self.fermion_operator)
77
78    def test_eigenspectrum(self):
79        fermion_eigenspectrum = eigenspectrum(self.fermion_operator)
80        qubit_eigenspectrum = eigenspectrum(self.qubit_operator)
81        interaction_eigenspectrum = eigenspectrum(self.interaction_operator)
82        for i in range(2**self.n_qubits):
83            self.assertAlmostEqual(fermion_eigenspectrum[i],
84                                   qubit_eigenspectrum[i])
85            self.assertAlmostEqual(fermion_eigenspectrum[i],
86                                   interaction_eigenspectrum[i])
87
88        with self.assertRaises(TypeError):
89            _ = eigenspectrum(BosonOperator())
90
91        with self.assertRaises(TypeError):
92            _ = eigenspectrum(QuadOperator())
93
94
95class SparseOperatorTest(unittest.TestCase):
96
97    def test_kronecker_operators(self):
98
99        self.assertAlmostEqual(
100            0,
101            numpy.amax(
102                numpy.absolute(
103                    kronecker_operators(3 * [identity_csc]) -
104                    kronecker_operators(3 * [pauli_x_csc])**2)))
105
106    def test_qubit_jw_fermion_integration(self):
107
108        # Initialize a random fermionic operator.
109        fermion_operator = FermionOperator(((3, 1), (2, 1), (1, 0), (0, 0)),
110                                           -4.3)
111        fermion_operator += FermionOperator(((3, 1), (1, 0)), 8.17)
112        fermion_operator += 3.2 * FermionOperator()
113
114        # Map to qubits and compare matrix versions.
115        qubit_operator = jordan_wigner(fermion_operator)
116        qubit_sparse = get_sparse_operator(qubit_operator)
117        qubit_spectrum = sparse_eigenspectrum(qubit_sparse)
118        fermion_sparse = jordan_wigner_sparse(fermion_operator)
119        fermion_spectrum = sparse_eigenspectrum(fermion_sparse)
120        self.assertAlmostEqual(
121            0., numpy.amax(numpy.absolute(fermion_spectrum - qubit_spectrum)))
122
123
124class JordanWignerSparseTest(unittest.TestCase):
125
126    def test_jw_sparse_0create(self):
127        expected = csc_matrix(([1], ([1], [0])), shape=(2, 2))
128        self.assertTrue(
129            numpy.allclose(
130                jordan_wigner_sparse(FermionOperator('0^')).A, expected.A))
131
132    def test_jw_sparse_1annihilate(self):
133        expected = csc_matrix(([1, -1], ([0, 2], [1, 3])), shape=(4, 4))
134        self.assertTrue(
135            numpy.allclose(
136                jordan_wigner_sparse(FermionOperator('1')).A, expected.A))
137
138    def test_jw_sparse_0create_2annihilate(self):
139        expected = csc_matrix(([-1j, 1j], ([4, 6], [1, 3])),
140                              shape=(8, 8),
141                              dtype=numpy.complex128)
142        self.assertTrue(
143            numpy.allclose(
144                jordan_wigner_sparse(FermionOperator('0^ 2', -1j)).A,
145                expected.A))
146
147    def test_jw_sparse_0create_3annihilate(self):
148        expected = csc_matrix(
149            ([-1j, 1j, 1j, -1j], ([8, 10, 12, 14], [1, 3, 5, 7])),
150            shape=(16, 16),
151            dtype=numpy.complex128)
152        self.assertTrue(
153            numpy.allclose(
154                jordan_wigner_sparse(FermionOperator('0^ 3', -1j)).A,
155                expected.A))
156
157    def test_jw_sparse_twobody(self):
158        expected = csc_matrix(([1, 1], ([6, 14], [5, 13])), shape=(16, 16))
159        self.assertTrue(
160            numpy.allclose(
161                jordan_wigner_sparse(FermionOperator('2^ 1^ 1 3')).A,
162                expected.A))
163
164    def test_qubit_operator_sparse_n_qubits_too_small(self):
165        with self.assertRaises(ValueError):
166            qubit_operator_sparse(QubitOperator('X3'), 1)
167
168    def test_qubit_operator_sparse_n_qubits_not_specified(self):
169        expected = csc_matrix(([1, 1, 1, 1], ([1, 0, 3, 2], [0, 1, 2, 3])),
170                              shape=(4, 4))
171        self.assertTrue(
172            numpy.allclose(
173                qubit_operator_sparse(QubitOperator('X1')).A, expected.A))
174
175    def test_get_linear_qubit_operator_diagonal_wrong_n(self):
176        """Testing with wrong n_qubits."""
177        with self.assertRaises(ValueError):
178            get_linear_qubit_operator_diagonal(QubitOperator('X3'), 1)
179
180    def test_get_linear_qubit_operator_diagonal_0(self):
181        """Testing with zero term."""
182        qubit_operator = QubitOperator.zero()
183        vec_expected = numpy.zeros(8)
184
185        self.assertTrue(
186            numpy.allclose(
187                get_linear_qubit_operator_diagonal(qubit_operator, 3),
188                vec_expected))
189
190    def test_get_linear_qubit_operator_diagonal_zero(self):
191        """Get zero diagonals from get_linear_qubit_operator_diagonal."""
192        qubit_operator = QubitOperator('X0 Y1')
193        vec_expected = numpy.zeros(4)
194
195        self.assertTrue(
196            numpy.allclose(get_linear_qubit_operator_diagonal(qubit_operator),
197                           vec_expected))
198
199    def test_get_linear_qubit_operator_diagonal_non_zero(self):
200        """Get non zero diagonals from get_linear_qubit_operator_diagonal."""
201        qubit_operator = QubitOperator('Z0 Z2')
202        vec_expected = numpy.array([1, -1, 1, -1, -1, 1, -1, 1])
203
204        self.assertTrue(
205            numpy.allclose(get_linear_qubit_operator_diagonal(qubit_operator),
206                           vec_expected))
207
208    def test_get_linear_qubit_operator_diagonal_cmp_zero(self):
209        """Compare get_linear_qubit_operator_diagonal with
210            get_linear_qubit_operator."""
211        qubit_operator = QubitOperator('Z1 X2 Y5')
212        vec_expected = numpy.diag(
213            LinearQubitOperator(qubit_operator) * numpy.eye(2**6))
214
215        self.assertTrue(
216            numpy.allclose(get_linear_qubit_operator_diagonal(qubit_operator),
217                           vec_expected))
218
219    def test_get_linear_qubit_operator_diagonal_cmp_non_zero(self):
220        """Compare get_linear_qubit_operator_diagonal with
221            get_linear_qubit_operator."""
222        qubit_operator = QubitOperator('Z1 Z2 Z5')
223        vec_expected = numpy.diag(
224            LinearQubitOperator(qubit_operator) * numpy.eye(2**6))
225
226        self.assertTrue(
227            numpy.allclose(get_linear_qubit_operator_diagonal(qubit_operator),
228                           vec_expected))
229
230
231class ComputationalBasisStateTest(unittest.TestCase):
232
233    def test_computational_basis_state(self):
234        comp_basis_state = jw_configuration_state([0, 2, 5], 7)
235        self.assertAlmostEqual(comp_basis_state[82], 1.)
236        self.assertAlmostEqual(sum(comp_basis_state), 1.)
237
238
239class JWHartreeFockStateTest(unittest.TestCase):
240
241    def test_jw_hartree_fock_state(self):
242        hartree_fock_state = jw_hartree_fock_state(3, 7)
243        self.assertAlmostEqual(hartree_fock_state[112], 1.)
244        self.assertAlmostEqual(sum(hartree_fock_state), 1.)
245
246
247class JWNumberIndicesTest(unittest.TestCase):
248
249    def test_jw_sparse_index(self):
250        """Test the indexing scheme for selecting specific particle numbers"""
251        expected = [1, 2]
252        calculated_indices = jw_number_indices(1, 2)
253        self.assertEqual(expected, calculated_indices)
254
255        expected = [3]
256        calculated_indices = jw_number_indices(2, 2)
257        self.assertEqual(expected, calculated_indices)
258
259    def test_jw_number_indices(self):
260        n_qubits = numpy.random.randint(1, 12)
261        n_particles = numpy.random.randint(n_qubits + 1)
262
263        number_indices = jw_number_indices(n_particles, n_qubits)
264        subspace_dimension = len(number_indices)
265
266        self.assertEqual(subspace_dimension, comb(n_qubits, n_particles))
267
268        for index in number_indices:
269            binary_string = bin(index)[2:].zfill(n_qubits)
270            n_ones = binary_string.count('1')
271            self.assertEqual(n_ones, n_particles)
272
273
274class JWSzIndicesTest(unittest.TestCase):
275
276    def test_jw_sz_indices(self):
277        """Test the indexing scheme for selecting specific sz value"""
278
279        def sz_integer(bitstring):
280            """Computes the total number of occupied up sites
281            minus the total number of occupied down sites."""
282            n_sites = len(bitstring) // 2
283
284            n_up = len([
285                site for site in range(n_sites)
286                if bitstring[up_index(site)] == '1'
287            ])
288            n_down = len([
289                site for site in range(n_sites)
290                if bitstring[down_index(site)] == '1'
291            ])
292
293            return n_up - n_down
294
295        def jw_sz_indices_brute_force(sz_value, n_qubits):
296            """Computes the correct indices by brute force."""
297            indices = []
298            for bitstring in itertools.product(['0', '1'], repeat=n_qubits):
299                if (sz_integer(bitstring) == int(2 * sz_value)):
300                    indices.append(int(''.join(bitstring), 2))
301
302            return indices
303
304        # General test
305        n_sites = numpy.random.randint(1, 10)
306        n_qubits = 2 * n_sites
307        sz_int = ((-1)**numpy.random.randint(2) *
308                  numpy.random.randint(n_sites + 1))
309        sz_value = sz_int / 2.
310
311        correct_indices = jw_sz_indices_brute_force(sz_value, n_qubits)
312        subspace_dimension = len(correct_indices)
313
314        calculated_indices = jw_sz_indices(sz_value, n_qubits)
315
316        self.assertEqual(len(calculated_indices), subspace_dimension)
317
318        for index in calculated_indices:
319            binary_string = bin(index)[2:].zfill(n_qubits)
320            self.assertEqual(sz_integer(binary_string), sz_int)
321
322        # Test fixing particle number
323        n_particles = abs(sz_int)
324
325        correct_indices = [
326            index for index in correct_indices
327            if bin(index)[2:].count('1') == n_particles
328        ]
329        subspace_dimension = len(correct_indices)
330
331        calculated_indices = jw_sz_indices(sz_value,
332                                           n_qubits,
333                                           n_electrons=n_particles)
334
335        self.assertEqual(len(calculated_indices), subspace_dimension)
336
337        for index in calculated_indices:
338            binary_string = bin(index)[2:].zfill(n_qubits)
339            self.assertEqual(sz_integer(binary_string), sz_int)
340            self.assertEqual(binary_string.count('1'), n_particles)
341
342        # Test exceptions
343        with self.assertRaises(ValueError):
344            jw_sz_indices(3, 3)
345
346        with self.assertRaises(ValueError):
347            jw_sz_indices(3.1, 4)
348
349        with self.assertRaises(ValueError):
350            jw_sz_indices(1.5, 8, n_electrons=6)
351
352        with self.assertRaises(ValueError):
353            jw_sz_indices(1.5, 8, n_electrons=1)
354
355
356class JWNumberRestrictOperatorTest(unittest.TestCase):
357
358    def test_jw_restrict_operator(self):
359        """Test the scheme for restricting JW encoded operators to number"""
360        # Make a Hamiltonian that cares mostly about number of electrons
361        n_qubits = 4
362        target_electrons = 2
363        penalty_const = 10.
364        number_sparse = jordan_wigner_sparse(number_operator(n_qubits))
365        bias_sparse = jordan_wigner_sparse(
366            sum([
367                FermionOperator(((i, 1), (i, 0)), 1.0) for i in range(n_qubits)
368            ], FermionOperator()))
369        hamiltonian_sparse = penalty_const * (
370            number_sparse -
371            target_electrons * scipy.sparse.identity(2**n_qubits)
372        ).dot(number_sparse - target_electrons *
373              scipy.sparse.identity(2**n_qubits)) + bias_sparse
374
375        restricted_hamiltonian = jw_number_restrict_operator(
376            hamiltonian_sparse, target_electrons, n_qubits)
377        true_eigvals, _ = eigh(hamiltonian_sparse.A)
378        test_eigvals, _ = eigh(restricted_hamiltonian.A)
379
380        self.assertAlmostEqual(norm(true_eigvals[:6] - test_eigvals[:6]), 0.0)
381
382    def test_jw_restrict_operator_hopping_to_1_particle(self):
383        hop = FermionOperator('3^ 1') + FermionOperator('1^ 3')
384        hop_sparse = jordan_wigner_sparse(hop, n_qubits=4)
385        hop_restrict = jw_number_restrict_operator(hop_sparse, 1, n_qubits=4)
386        expected = csc_matrix(([1, 1], ([0, 2], [2, 0])), shape=(4, 4))
387
388        self.assertTrue(numpy.allclose(hop_restrict.A, expected.A))
389
390    def test_jw_restrict_operator_interaction_to_1_particle(self):
391        interaction = FermionOperator('3^ 2^ 4 1')
392        interaction_sparse = jordan_wigner_sparse(interaction, n_qubits=6)
393        interaction_restrict = jw_number_restrict_operator(interaction_sparse,
394                                                           1,
395                                                           n_qubits=6)
396        expected = csc_matrix(([], ([], [])), shape=(6, 6))
397
398        self.assertTrue(numpy.allclose(interaction_restrict.A, expected.A))
399
400    def test_jw_restrict_operator_interaction_to_2_particles(self):
401        interaction = (FermionOperator('3^ 2^ 4 1') +
402                       FermionOperator('4^ 1^ 3 2'))
403        interaction_sparse = jordan_wigner_sparse(interaction, n_qubits=6)
404        interaction_restrict = jw_number_restrict_operator(interaction_sparse,
405                                                           2,
406                                                           n_qubits=6)
407
408        dim = 6 * 5 // 2  # shape of new sparse array
409
410        # 3^ 2^ 4 1 maps 2**4 + 2 = 18 to 2**3 + 2**2 = 12 and vice versa;
411        # in the 2-particle subspace (1, 4) and (2, 3) are 7th and 9th.
412        expected = csc_matrix(([-1, -1], ([7, 9], [9, 7])), shape=(dim, dim))
413
414        self.assertTrue(numpy.allclose(interaction_restrict.A, expected.A))
415
416    def test_jw_restrict_operator_hopping_to_1_particle_default_nqubits(self):
417        interaction = (FermionOperator('3^ 2^ 4 1') +
418                       FermionOperator('4^ 1^ 3 2'))
419        interaction_sparse = jordan_wigner_sparse(interaction, n_qubits=6)
420        # n_qubits should default to 6
421        interaction_restrict = jw_number_restrict_operator(
422            interaction_sparse, 2)
423
424        dim = 6 * 5 // 2  # shape of new sparse array
425
426        # 3^ 2^ 4 1 maps 2**4 + 2 = 18 to 2**3 + 2**2 = 12 and vice versa;
427        # in the 2-particle subspace (1, 4) and (2, 3) are 7th and 9th.
428        expected = csc_matrix(([-1, -1], ([7, 9], [9, 7])), shape=(dim, dim))
429
430        self.assertTrue(numpy.allclose(interaction_restrict.A, expected.A))
431
432    def test_jw_restrict_jellium_ground_state_integration(self):
433        n_qubits = 4
434        grid = Grid(dimensions=1, length=n_qubits, scale=1.0)
435        jellium_hamiltonian = jordan_wigner_sparse(
436            jellium_model(grid, spinless=False))
437
438        #  2 * n_qubits because of spin
439        number_sparse = jordan_wigner_sparse(number_operator(2 * n_qubits))
440
441        restricted_number = jw_number_restrict_operator(number_sparse, 2)
442        restricted_jellium_hamiltonian = jw_number_restrict_operator(
443            jellium_hamiltonian, 2)
444
445        _, ground_state = get_ground_state(restricted_jellium_hamiltonian)
446
447        number_expectation = expectation(restricted_number, ground_state)
448        self.assertAlmostEqual(number_expectation, 2)
449
450
451class JWSzRestrictOperatorTest(unittest.TestCase):
452
453    def test_restrict_interaction_hamiltonian(self):
454        """Test restricting a coulomb repulsion Hamiltonian to a specified
455        Sz manifold."""
456        x_dim = 3
457        y_dim = 2
458
459        interaction_term = fermi_hubbard(x_dim, y_dim, 0., 1.)
460        interaction_sparse = get_sparse_operator(interaction_term)
461        sz_value = 2
462        interaction_restricted = jw_sz_restrict_operator(
463            interaction_sparse, sz_value)
464        restricted_interaction_values = set(
465            [int(value.real) for value in interaction_restricted.diagonal()])
466        # Originally the eigenvalues run from 0 to 6 but after restricting,
467        # they should run from 0 to 2
468        self.assertEqual(restricted_interaction_values, {0, 1, 2})
469
470
471class JWNumberRestrictStateTest(unittest.TestCase):
472
473    def test_jw_number_restrict_state(self):
474        n_qubits = numpy.random.randint(1, 12)
475        n_particles = numpy.random.randint(0, n_qubits)
476
477        number_indices = jw_number_indices(n_particles, n_qubits)
478        subspace_dimension = len(number_indices)
479
480        # Create a vector that has entry 1 for every coordinate with
481        # the specified particle number, and 0 everywhere else
482        vector = numpy.zeros(2**n_qubits, dtype=float)
483        vector[number_indices] = 1
484
485        # Restrict the vector
486        restricted_vector = jw_number_restrict_state(vector, n_particles)
487
488        # Check that it has the correct shape
489        self.assertEqual(restricted_vector.shape[0], subspace_dimension)
490
491        # Check that it has the same norm as the original vector
492        self.assertAlmostEqual(
493            inner_product(vector, vector),
494            inner_product(restricted_vector, restricted_vector))
495
496
497class JWSzRestrictStateTest(unittest.TestCase):
498
499    def test_jw_sz_restrict_state(self):
500        n_sites = numpy.random.randint(1, 10)
501        n_qubits = 2 * n_sites
502        sz_int = ((-1)**numpy.random.randint(2) *
503                  numpy.random.randint(n_sites + 1))
504        sz_value = sz_int / 2
505
506        sz_indices = jw_sz_indices(sz_value, n_qubits)
507        subspace_dimension = len(sz_indices)
508
509        # Create a vector that has entry 1 for every coordinate in
510        # the specified subspace, and 0 everywhere else
511        vector = numpy.zeros(2**n_qubits, dtype=float)
512        vector[sz_indices] = 1
513
514        # Restrict the vector
515        restricted_vector = jw_sz_restrict_state(vector, sz_value)
516
517        # Check that it has the correct shape
518        self.assertEqual(restricted_vector.shape[0], subspace_dimension)
519
520        # Check that it has the same norm as the original vector
521        self.assertAlmostEqual(
522            inner_product(vector, vector),
523            inner_product(restricted_vector, restricted_vector))
524
525
526class JWGetGroundStatesByParticleNumberTest(unittest.TestCase):
527
528    def test_jw_get_ground_state_at_particle_number_herm_conserving(self):
529        # Initialize a particle-number-conserving Hermitian operator
530        ferm_op = FermionOperator('0^ 1') + FermionOperator('1^ 0') + \
531            FermionOperator('1^ 2') + FermionOperator('2^ 1') + \
532            FermionOperator('1^ 3', -.4) + FermionOperator('3^ 1', -.4)
533        jw_hamiltonian = jordan_wigner(ferm_op)
534        sparse_operator = get_sparse_operator(jw_hamiltonian)
535        n_qubits = 4
536
537        num_op = get_sparse_operator(number_operator(n_qubits))
538
539        # Test each possible particle number
540        for particle_number in range(n_qubits):
541            # Get the ground energy and ground state at this particle number
542            energy, state = jw_get_ground_state_at_particle_number(
543                sparse_operator, particle_number)
544
545            # Check that it's an eigenvector with the correct eigenvalue
546            self.assertTrue(
547                numpy.allclose(sparse_operator.dot(state), energy * state))
548
549            # Check that it has the correct particle number
550            num = expectation(num_op, state)
551            self.assertAlmostEqual(num, particle_number)
552
553    def test_jw_get_ground_state_at_particle_number_hubbard(self):
554
555        model = fermi_hubbard(2, 2, 1.0, 4.0)
556        sparse_operator = get_sparse_operator(model)
557        n_qubits = count_qubits(model)
558        num_op = get_sparse_operator(number_operator(n_qubits))
559
560        # Test each possible particle number
561        for particle_number in range(n_qubits):
562            # Get the ground energy and ground state at this particle number
563            energy, state = jw_get_ground_state_at_particle_number(
564                sparse_operator, particle_number)
565
566            # Check that it's an eigenvector with the correct eigenvalue
567            self.assertTrue(
568                numpy.allclose(sparse_operator.dot(state), energy * state))
569
570            # Check that it has the correct particle number
571            num = expectation(num_op, state)
572            self.assertAlmostEqual(num, particle_number)
573
574    def test_jw_get_ground_state_at_particle_number_jellium(self):
575
576        grid = Grid(2, 2, 1.0)
577        model = jellium_model(grid, spinless=True, plane_wave=False)
578        sparse_operator = get_sparse_operator(model)
579        n_qubits = count_qubits(model)
580        num_op = get_sparse_operator(number_operator(n_qubits))
581
582        # Test each possible particle number
583        for particle_number in range(n_qubits):
584            # Get the ground energy and ground state at this particle number
585            energy, state = jw_get_ground_state_at_particle_number(
586                sparse_operator, particle_number)
587
588            # Check that it's an eigenvector with the correct eigenvalue
589            self.assertTrue(
590                numpy.allclose(sparse_operator.dot(state), energy * state))
591
592            # Check that it has the correct particle number
593            num = expectation(num_op, state)
594            self.assertAlmostEqual(num, particle_number)
595
596
597class JWSparseGivensRotationTest(unittest.TestCase):
598
599    def test_bad_input(self):
600        with self.assertRaises(ValueError):
601            jw_sparse_givens_rotation(0, 2, 1., 1., 5)
602        with self.assertRaises(ValueError):
603            jw_sparse_givens_rotation(4, 5, 1., 1., 5)
604
605
606class GroundStateTest(unittest.TestCase):
607
608    def test_get_ground_state_hermitian(self):
609        ground = get_ground_state(
610            get_sparse_operator(
611                QubitOperator('Y0 X1') + QubitOperator('Z0 Z1')))
612        expected_state = csc_matrix(([1j, 1], ([1, 2], [0, 0])),
613                                    shape=(4, 1),
614                                    dtype=numpy.complex128).A
615        expected_state /= numpy.sqrt(2.0)
616
617        self.assertAlmostEqual(ground[0], -2)
618        self.assertAlmostEqual(
619            numpy.absolute(expected_state.T.conj().dot(ground[1]))[0], 1.)
620
621
622class ExpectationTest(unittest.TestCase):
623
624    def test_expectation_correct_sparse_matrix(self):
625        operator = get_sparse_operator(QubitOperator('X0'), n_qubits=2)
626        vector = numpy.array([0., 1.j, 0., 1.j])
627        self.assertAlmostEqual(expectation(operator, vector), 2.0)
628
629        density_matrix = scipy.sparse.csc_matrix(
630            numpy.outer(vector, numpy.conjugate(vector)))
631        self.assertAlmostEqual(expectation(operator, density_matrix), 2.0)
632
633    def test_expectation_correct_linear_operator(self):
634        operator = LinearQubitOperator(QubitOperator('X0'), n_qubits=2)
635        vector = numpy.array([0., 1.j, 0., 1.j])
636        self.assertAlmostEqual(expectation(operator, vector), 2.0)
637
638    def test_expectation_handles_column_vector(self):
639        operator = get_sparse_operator(QubitOperator('X0'), n_qubits=2)
640        vector = numpy.array([[0.], [1.j], [0.], [1.j]])
641        self.assertAlmostEqual(expectation(operator, vector), 2.0)
642
643    def test_expectation_correct_zero(self):
644        operator = get_sparse_operator(QubitOperator('X0'), n_qubits=2)
645        vector = numpy.array([1j, -1j, -1j, -1j])
646        self.assertAlmostEqual(expectation(operator, vector), 0.0)
647
648    def test_execptions(self):
649        operator = LinearQubitOperator(QubitOperator('X0'), n_qubits=2)
650        vector = scipy.sparse.csc_matrix(numpy.array([1, 0, 0, 0]))
651        with self.assertRaises(ValueError):
652            _ = expectation(operator, vector)
653        with self.assertRaises(ValueError):
654            _ = expectation(vector, operator)
655
656
657class VarianceTest(unittest.TestCase):
658
659    def test_variance_row_vector(self):
660
661        X = pauli_matrix_map['X']
662        Z = pauli_matrix_map['Z']
663        zero = numpy.array([1., 0.])
664        plus = numpy.array([1., 1.]) / numpy.sqrt(2)
665        minus = numpy.array([1., -1.]) / numpy.sqrt(2)
666
667        self.assertAlmostEqual(variance(Z, zero), 0.)
668        self.assertAlmostEqual(variance(X, zero), 1.)
669
670        self.assertAlmostEqual(variance(Z, plus), 1.)
671        self.assertAlmostEqual(variance(X, plus), 0.)
672
673        self.assertAlmostEqual(variance(Z, minus), 1.)
674        self.assertAlmostEqual(variance(X, minus), 0.)
675
676    def test_variance_column_vector(self):
677
678        X = pauli_matrix_map['X']
679        Z = pauli_matrix_map['Z']
680        zero = numpy.array([[1.], [0.]])
681        plus = numpy.array([[1.], [1.]]) / numpy.sqrt(2)
682        minus = numpy.array([[1.], [-1.]]) / numpy.sqrt(2)
683
684        self.assertAlmostEqual(variance(Z, zero), 0.)
685        self.assertAlmostEqual(variance(X, zero), 1.)
686
687        self.assertAlmostEqual(variance(Z, plus), 1.)
688        self.assertAlmostEqual(variance(X, plus), 0.)
689
690        self.assertAlmostEqual(variance(Z, minus), 1.)
691        self.assertAlmostEqual(variance(X, minus), 0.)
692
693
694class ExpectationComputationalBasisStateTest(unittest.TestCase):
695
696    def test_expectation_fermion_operator_single_number_terms(self):
697        operator = FermionOperator('3^ 3', 1.9) + FermionOperator('2^ 1')
698        state = csc_matrix(([1], ([15], [0])), shape=(16, 1))
699
700        self.assertAlmostEqual(
701            expectation_computational_basis_state(operator, state), 1.9)
702
703    def test_expectation_fermion_operator_two_number_terms(self):
704        operator = (FermionOperator('2^ 2', 1.9) + FermionOperator('2^ 1') +
705                    FermionOperator('2^ 1^ 2 1', -1.7))
706        state = csc_matrix(([1], ([6], [0])), shape=(16, 1))
707
708        self.assertAlmostEqual(
709            expectation_computational_basis_state(operator, state), 3.6)
710
711    def test_expectation_identity_fermion_operator(self):
712        operator = FermionOperator.identity() * 1.1
713        state = csc_matrix(([1], ([6], [0])), shape=(16, 1))
714
715        self.assertAlmostEqual(
716            expectation_computational_basis_state(operator, state), 1.1)
717
718    def test_expectation_state_is_list_single_number_terms(self):
719        operator = FermionOperator('3^ 3', 1.9) + FermionOperator('2^ 1')
720        state = [1, 1, 1, 1]
721
722        self.assertAlmostEqual(
723            expectation_computational_basis_state(operator, state), 1.9)
724
725    def test_expectation_state_is_list_fermion_operator_two_number_terms(self):
726        operator = (FermionOperator('2^ 2', 1.9) + FermionOperator('2^ 1') +
727                    FermionOperator('2^ 1^ 2 1', -1.7))
728        state = [0, 1, 1]
729
730        self.assertAlmostEqual(
731            expectation_computational_basis_state(operator, state), 3.6)
732
733    def test_expectation_state_is_list_identity_fermion_operator(self):
734        operator = FermionOperator.identity() * 1.1
735        state = [0, 1, 1]
736
737        self.assertAlmostEqual(
738            expectation_computational_basis_state(operator, state), 1.1)
739
740    def test_expectation_bad_operator_type(self):
741        with self.assertRaises(TypeError):
742            expectation_computational_basis_state(
743                'never', csc_matrix(([1], ([6], [0])), shape=(16, 1)))
744
745    def test_expectation_qubit_operator_not_implemented(self):
746        with self.assertRaises(NotImplementedError):
747            expectation_computational_basis_state(
748                QubitOperator(), csc_matrix(([1], ([6], [0])), shape=(16, 1)))
749
750
751class ExpectationDualBasisOperatorWithPlaneWaveBasisState(unittest.TestCase):
752
753    def setUp(self):
754        grid_length = 4
755        dimension = 1
756        wigner_seitz_radius = 10.
757        self.spinless = True
758        self.n_spatial_orbitals = grid_length**dimension
759
760        n_qubits = self.n_spatial_orbitals
761        self.n_particles = 3
762
763        # Compute appropriate length scale and the corresponding grid.
764        length_scale = wigner_seitz_length_scale(wigner_seitz_radius,
765                                                 self.n_particles, dimension)
766
767        self.grid1 = Grid(dimension, grid_length, length_scale)
768        # Get the occupied orbitals of the plane-wave basis Hartree-Fock state.
769        hamiltonian = jellium_model(self.grid1, self.spinless, plane_wave=True)
770        hamiltonian = normal_ordered(hamiltonian)
771        hamiltonian.compress()
772
773        occupied_states = numpy.array(
774            lowest_single_particle_energy_states(hamiltonian, self.n_particles))
775        self.hf_state_index1 = numpy.sum(2**occupied_states)
776
777        self.hf_state1 = numpy.zeros(2**n_qubits)
778        self.hf_state1[self.hf_state_index1] = 1.0
779
780        self.orbital_occupations1 = [
781            digit == '1' for digit in bin(self.hf_state_index1)[2:]
782        ][::-1]
783        self.occupied_orbitals1 = [
784            index for index, occupied in enumerate(self.orbital_occupations1)
785            if occupied
786        ]
787
788        self.reversed_occupied_orbitals1 = list(self.occupied_orbitals1)
789        for i in range(len(self.reversed_occupied_orbitals1)):
790            self.reversed_occupied_orbitals1[i] = -1 + int(
791                numpy.log2(self.hf_state1.shape[0])
792            ) - self.reversed_occupied_orbitals1[i]
793
794        self.reversed_hf_state_index1 = sum(
795            2**index for index in self.reversed_occupied_orbitals1)
796
797    def test_1body_hopping_operator_1D(self):
798        operator = FermionOperator('2^ 0')
799        operator = normal_ordered(operator)
800        transformed_operator = normal_ordered(
801            fourier_transform(operator, self.grid1, self.spinless))
802
803        expected = expectation(get_sparse_operator(transformed_operator),
804                               self.hf_state1)
805        actual = expectation_db_operator_with_pw_basis_state(
806            operator, self.reversed_occupied_orbitals1, self.n_spatial_orbitals,
807            self.grid1, self.spinless)
808        self.assertAlmostEqual(expected, actual)
809
810    def test_1body_number_operator_1D(self):
811        operator = FermionOperator('2^ 2')
812        operator = normal_ordered(operator)
813        transformed_operator = normal_ordered(
814            fourier_transform(operator, self.grid1, self.spinless))
815
816        expected = expectation(get_sparse_operator(transformed_operator),
817                               self.hf_state1)
818        actual = expectation_db_operator_with_pw_basis_state(
819            operator, self.reversed_occupied_orbitals1, self.n_spatial_orbitals,
820            self.grid1, self.spinless)
821        self.assertAlmostEqual(expected, actual)
822
823    def test_2body_partial_number_operator_high_1D(self):
824        operator = FermionOperator('2^ 1^ 2 0')
825        operator = normal_ordered(operator)
826        transformed_operator = normal_ordered(
827            fourier_transform(operator, self.grid1, self.spinless))
828
829        expected = expectation(get_sparse_operator(transformed_operator),
830                               self.hf_state1)
831        actual = expectation_db_operator_with_pw_basis_state(
832            operator, self.reversed_occupied_orbitals1, self.n_spatial_orbitals,
833            self.grid1, self.spinless)
834        self.assertAlmostEqual(expected, actual)
835
836    def test_2body_partial_number_operator_mid_1D(self):
837        operator = FermionOperator('1^ 0^ 1 2')
838        operator = normal_ordered(operator)
839        transformed_operator = normal_ordered(
840            fourier_transform(operator, self.grid1, self.spinless))
841
842        expected = expectation(get_sparse_operator(transformed_operator),
843                               self.hf_state1)
844        actual = expectation_db_operator_with_pw_basis_state(
845            operator, self.reversed_occupied_orbitals1, self.n_spatial_orbitals,
846            self.grid1, self.spinless)
847        self.assertAlmostEqual(expected, actual)
848
849    def test_3body_double_number_operator_1D(self):
850        operator = FermionOperator('3^ 2^ 1^ 3 1 0')
851        operator = normal_ordered(operator)
852        transformed_operator = normal_ordered(
853            fourier_transform(operator, self.grid1, self.spinless))
854
855        expected = expectation(get_sparse_operator(transformed_operator),
856                               self.hf_state1)
857        actual = expectation_db_operator_with_pw_basis_state(
858            operator, self.reversed_occupied_orbitals1, self.n_spatial_orbitals,
859            self.grid1, self.spinless)
860        self.assertAlmostEqual(expected, actual)
861
862    def test_2body_adjacent_number_operator_1D(self):
863        operator = FermionOperator('3^ 2^ 2 1')
864        operator = normal_ordered(operator)
865        transformed_operator = normal_ordered(
866            fourier_transform(operator, self.grid1, self.spinless))
867
868        expected = expectation(get_sparse_operator(transformed_operator),
869                               self.hf_state1)
870        actual = expectation_db_operator_with_pw_basis_state(
871            operator, self.reversed_occupied_orbitals1, self.n_spatial_orbitals,
872            self.grid1, self.spinless)
873        self.assertAlmostEqual(expected, actual)
874
875    def test_1d5_with_spin_10particles(self):
876        dimension = 1
877        grid_length = 5
878        n_spatial_orbitals = grid_length**dimension
879        wigner_seitz_radius = 9.3
880
881        spinless = False
882        n_qubits = n_spatial_orbitals
883        if not spinless:
884            n_qubits *= 2
885        n_particles_big = 10
886
887        length_scale = wigner_seitz_length_scale(wigner_seitz_radius,
888                                                 n_particles_big, dimension)
889
890        self.grid3 = Grid(dimension, grid_length, length_scale)
891        # Get the occupied orbitals of the plane-wave basis Hartree-Fock state.
892        hamiltonian = jellium_model(self.grid3, spinless, plane_wave=True)
893        hamiltonian = normal_ordered(hamiltonian)
894        hamiltonian.compress()
895
896        occupied_states = numpy.array(
897            lowest_single_particle_energy_states(hamiltonian, n_particles_big))
898        self.hf_state_index3 = numpy.sum(2**occupied_states)
899
900        self.hf_state3 = csc_matrix(([1.0], ([self.hf_state_index3], [0])),
901                                    shape=(2**n_qubits, 1))
902
903        self.orbital_occupations3 = [
904            digit == '1' for digit in bin(self.hf_state_index3)[2:]
905        ][::-1]
906        self.occupied_orbitals3 = [
907            index for index, occupied in enumerate(self.orbital_occupations3)
908            if occupied
909        ]
910
911        self.reversed_occupied_orbitals3 = list(self.occupied_orbitals3)
912        for i in range(len(self.reversed_occupied_orbitals3)):
913            self.reversed_occupied_orbitals3[i] = -1 + int(
914                numpy.log2(self.hf_state3.shape[0])
915            ) - self.reversed_occupied_orbitals3[i]
916
917        self.reversed_hf_state_index3 = sum(
918            2**index for index in self.reversed_occupied_orbitals3)
919
920        operator = (FermionOperator('6^ 0^ 1^ 3 5 4', 2) +
921                    FermionOperator('7^ 6^ 5 4', -3.7j) +
922                    FermionOperator('3^ 3', 2.1) + FermionOperator('3^ 2', 1.7))
923        operator = normal_ordered(operator)
924        normal_ordered(fourier_transform(operator, self.grid3, spinless))
925
926        expected = 2.1
927        # Calculated from expectation(get_sparse_operator(
928        #    transformed_operator), self.hf_state3)
929        actual = expectation_db_operator_with_pw_basis_state(
930            operator, self.reversed_occupied_orbitals3, n_spatial_orbitals,
931            self.grid3, spinless)
932
933        self.assertAlmostEqual(expected, actual)
934
935    def test_1d5_with_spin_7particles(self):
936        dimension = 1
937        grid_length = 5
938        n_spatial_orbitals = grid_length**dimension
939        wigner_seitz_radius = 9.3
940
941        spinless = False
942        n_qubits = n_spatial_orbitals
943        if not spinless:
944            n_qubits *= 2
945        n_particles_big = 7
946
947        length_scale = wigner_seitz_length_scale(wigner_seitz_radius,
948                                                 n_particles_big, dimension)
949
950        self.grid3 = Grid(dimension, grid_length, length_scale)
951        # Get the occupied orbitals of the plane-wave basis Hartree-Fock state.
952        hamiltonian = jellium_model(self.grid3, spinless, plane_wave=True)
953        hamiltonian = normal_ordered(hamiltonian)
954        hamiltonian.compress()
955
956        occupied_states = numpy.array(
957            lowest_single_particle_energy_states(hamiltonian, n_particles_big))
958        self.hf_state_index3 = numpy.sum(2**occupied_states)
959
960        self.hf_state3 = csc_matrix(([1.0], ([self.hf_state_index3], [0])),
961                                    shape=(2**n_qubits, 1))
962
963        self.orbital_occupations3 = [
964            digit == '1' for digit in bin(self.hf_state_index3)[2:]
965        ][::-1]
966        self.occupied_orbitals3 = [
967            index for index, occupied in enumerate(self.orbital_occupations3)
968            if occupied
969        ]
970
971        self.reversed_occupied_orbitals3 = list(self.occupied_orbitals3)
972        for i in range(len(self.reversed_occupied_orbitals3)):
973            self.reversed_occupied_orbitals3[i] = -1 + int(
974                numpy.log2(self.hf_state3.shape[0])
975            ) - self.reversed_occupied_orbitals3[i]
976
977        self.reversed_hf_state_index3 = sum(
978            2**index for index in self.reversed_occupied_orbitals3)
979
980        operator = (FermionOperator('6^ 0^ 1^ 3 5 4', 2) +
981                    FermionOperator('7^ 2^ 4 1') +
982                    FermionOperator('3^ 3', 2.1) +
983                    FermionOperator('5^ 3^ 1 0', 7.3))
984        operator = normal_ordered(operator)
985        normal_ordered(fourier_transform(operator, self.grid3, spinless))
986
987        expected = 1.66 - 0.0615536707435j
988        # Calculated with expected = expectation(get_sparse_operator(
989        #    transformed_operator), self.hf_state3)
990        actual = expectation_db_operator_with_pw_basis_state(
991            operator, self.reversed_occupied_orbitals3, n_spatial_orbitals,
992            self.grid3, spinless)
993
994        self.assertAlmostEqual(expected, actual)
995
996    def test_3d2_spinless(self):
997        dimension = 3
998        grid_length = 2
999        n_spatial_orbitals = grid_length**dimension
1000        wigner_seitz_radius = 9.3
1001
1002        spinless = True
1003        n_qubits = n_spatial_orbitals
1004        n_particles_big = 5
1005
1006        length_scale = wigner_seitz_length_scale(wigner_seitz_radius,
1007                                                 n_particles_big, dimension)
1008
1009        self.grid3 = Grid(dimension, grid_length, length_scale)
1010        # Get the occupied orbitals of the plane-wave basis Hartree-Fock state.
1011        hamiltonian = jellium_model(self.grid3, spinless, plane_wave=True)
1012        hamiltonian = normal_ordered(hamiltonian)
1013        hamiltonian.compress()
1014
1015        occupied_states = numpy.array(
1016            lowest_single_particle_energy_states(hamiltonian, n_particles_big))
1017        self.hf_state_index3 = numpy.sum(2**occupied_states)
1018
1019        self.hf_state3 = csc_matrix(([1.0], ([self.hf_state_index3], [0])),
1020                                    shape=(2**n_qubits, 1))
1021
1022        self.orbital_occupations3 = [
1023            digit == '1' for digit in bin(self.hf_state_index3)[2:]
1024        ][::-1]
1025        self.occupied_orbitals3 = [
1026            index for index, occupied in enumerate(self.orbital_occupations3)
1027            if occupied
1028        ]
1029
1030        self.reversed_occupied_orbitals3 = list(self.occupied_orbitals3)
1031        for i in range(len(self.reversed_occupied_orbitals3)):
1032            self.reversed_occupied_orbitals3[i] = -1 + int(
1033                numpy.log2(self.hf_state3.shape[0])
1034            ) - self.reversed_occupied_orbitals3[i]
1035
1036        self.reversed_hf_state_index3 = sum(
1037            2**index for index in self.reversed_occupied_orbitals3)
1038
1039        operator = (FermionOperator('4^ 2^ 3^ 5 5 4', 2) +
1040                    FermionOperator('7^ 6^ 7 4', -3.7j) +
1041                    FermionOperator('3^ 7', 2.1))
1042        operator = normal_ordered(operator)
1043        normal_ordered(fourier_transform(operator, self.grid3, spinless))
1044        expected = -0.2625 - 0.4625j
1045        # Calculated with expectation(get_sparse_operator(
1046        #    transformed_operator), self.hf_state3)
1047        actual = expectation_db_operator_with_pw_basis_state(
1048            operator, self.reversed_occupied_orbitals3, n_spatial_orbitals,
1049            self.grid3, spinless)
1050
1051        self.assertAlmostEqual(expected, actual)
1052
1053    def test_3d2_with_spin(self):
1054        dimension = 3
1055        grid_length = 2
1056        n_spatial_orbitals = grid_length**dimension
1057        wigner_seitz_radius = 9.3
1058
1059        spinless = False
1060        n_qubits = n_spatial_orbitals
1061        if not spinless:
1062            n_qubits *= 2
1063        n_particles_big = 9
1064
1065        length_scale = wigner_seitz_length_scale(wigner_seitz_radius,
1066                                                 n_particles_big, dimension)
1067
1068        self.grid3 = Grid(dimension, grid_length, length_scale)
1069        # Get the occupied orbitals of the plane-wave basis Hartree-Fock state.
1070        hamiltonian = jellium_model(self.grid3, spinless, plane_wave=True)
1071        hamiltonian = normal_ordered(hamiltonian)
1072        hamiltonian.compress()
1073
1074        occupied_states = numpy.array(
1075            lowest_single_particle_energy_states(hamiltonian, n_particles_big))
1076        self.hf_state_index3 = numpy.sum(2**occupied_states)
1077
1078        self.hf_state3 = csc_matrix(([1.0], ([self.hf_state_index3], [0])),
1079                                    shape=(2**n_qubits, 1))
1080
1081        self.orbital_occupations3 = [
1082            digit == '1' for digit in bin(self.hf_state_index3)[2:]
1083        ][::-1]
1084        self.occupied_orbitals3 = [
1085            index for index, occupied in enumerate(self.orbital_occupations3)
1086            if occupied
1087        ]
1088
1089        self.reversed_occupied_orbitals3 = list(self.occupied_orbitals3)
1090        for i in range(len(self.reversed_occupied_orbitals3)):
1091            self.reversed_occupied_orbitals3[i] = -1 + int(
1092                numpy.log2(self.hf_state3.shape[0])
1093            ) - self.reversed_occupied_orbitals3[i]
1094
1095        self.reversed_hf_state_index3 = sum(
1096            2**index for index in self.reversed_occupied_orbitals3)
1097
1098        operator = (FermionOperator('4^ 2^ 3^ 5 5 4', 2) +
1099                    FermionOperator('7^ 6^ 7 4', -3.7j) +
1100                    FermionOperator('3^ 7', 2.1))
1101        operator = normal_ordered(operator)
1102        normal_ordered(fourier_transform(operator, self.grid3, spinless))
1103
1104        expected = -0.2625 - 0.578125j
1105        # Calculated from expected = expectation(get_sparse_operator(
1106        #    transformed_operator), self.hf_state3)
1107        actual = expectation_db_operator_with_pw_basis_state(
1108            operator, self.reversed_occupied_orbitals3, n_spatial_orbitals,
1109            self.grid3, spinless)
1110
1111        self.assertAlmostEqual(expected, actual)
1112
1113
1114class GetGapTest(unittest.TestCase):
1115
1116    def test_get_gap(self):
1117        operator = QubitOperator('Y0 X1') + QubitOperator('Z0 Z1')
1118        self.assertAlmostEqual(get_gap(get_sparse_operator(operator)), 2.0)
1119
1120    def test_get_gap_nonhermitian_error(self):
1121        operator = (QubitOperator('X0 Y1', 1 + 1j) +
1122                    QubitOperator('Z0 Z1', 1j) + QubitOperator((), 2 + 1j))
1123        with self.assertRaises(ValueError):
1124            get_gap(get_sparse_operator(operator))
1125
1126
1127class InnerProductTest(unittest.TestCase):
1128
1129    def test_inner_product(self):
1130        state_1 = numpy.array([1., 1.j])
1131        state_2 = numpy.array([1., -1.j])
1132
1133        self.assertAlmostEqual(inner_product(state_1, state_1), 2.)
1134        self.assertAlmostEqual(inner_product(state_1, state_2), 0.)
1135
1136
1137class BosonSparseTest(unittest.TestCase):
1138
1139    def setUp(self):
1140        self.hbar = 1.
1141        self.d = 5
1142        self.b = numpy.diag(numpy.sqrt(numpy.arange(1, self.d)), 1)
1143        self.bd = self.b.conj().T
1144        self.q = numpy.sqrt(self.hbar / 2) * (self.b + self.bd)
1145        self.p = -1j * numpy.sqrt(self.hbar / 2) * (self.b - self.bd)
1146        self.Id = numpy.identity(self.d)
1147
1148    def test_boson_ladder_noninteger_trunc(self):
1149        with self.assertRaises(ValueError):
1150            boson_ladder_sparse(1, 0, 0, 0.1)
1151
1152        with self.assertRaises(ValueError):
1153            boson_ladder_sparse(1, 0, 0, -1)
1154
1155        with self.assertRaises(ValueError):
1156            boson_ladder_sparse(1, 0, 0, 0)
1157
1158    def test_boson_ladder_destroy_one_mode(self):
1159        b = boson_ladder_sparse(1, 0, 0, self.d).toarray()
1160        self.assertTrue(numpy.allclose(b, self.b))
1161
1162    def test_boson_ladder_create_one_mode(self):
1163        bd = boson_ladder_sparse(1, 0, 1, self.d).toarray()
1164        self.assertTrue(numpy.allclose(bd, self.bd))
1165
1166    def test_boson_ladder_single_adjoint(self):
1167        b = boson_ladder_sparse(1, 0, 0, self.d).toarray()
1168        bd = boson_ladder_sparse(1, 0, 1, self.d).toarray()
1169        self.assertTrue(numpy.allclose(b.conj().T, bd))
1170
1171    def test_boson_ladder_two_mode(self):
1172        res = boson_ladder_sparse(2, 0, 0, self.d).toarray()
1173        expected = numpy.kron(self.b, self.Id)
1174        self.assertTrue(numpy.allclose(res, expected))
1175
1176        res = boson_ladder_sparse(2, 1, 0, self.d).toarray()
1177        expected = numpy.kron(self.Id, self.b)
1178        self.assertTrue(numpy.allclose(res, expected))
1179
1180    def test_single_quad_noninteger_trunc(self):
1181        with self.assertRaises(ValueError):
1182            single_quad_op_sparse(1, 0, 'q', self.hbar, 0.1)
1183
1184        with self.assertRaises(ValueError):
1185            single_quad_op_sparse(1, 0, 'q', self.hbar, -1)
1186
1187        with self.assertRaises(ValueError):
1188            single_quad_op_sparse(1, 0, 'q', self.hbar, 0)
1189
1190    def test_single_quad_q_one_mode(self):
1191        res = single_quad_op_sparse(1, 0, 'q', self.hbar, self.d).toarray()
1192        self.assertTrue(numpy.allclose(res, self.q))
1193        self.assertTrue(numpy.allclose(res, res.conj().T))
1194
1195    def test_single_quad_p_one_mode(self):
1196        res = single_quad_op_sparse(1, 0, 'p', self.hbar, self.d).toarray()
1197        self.assertTrue(numpy.allclose(res, self.p))
1198        self.assertTrue(numpy.allclose(res, res.conj().T))
1199
1200    def test_single_quad_two_mode(self):
1201        res = single_quad_op_sparse(2, 0, 'q', self.hbar, self.d).toarray()
1202        expected = numpy.kron(self.q, self.Id)
1203        self.assertTrue(numpy.allclose(res, expected))
1204
1205        res = single_quad_op_sparse(2, 1, 'p', self.hbar, self.d).toarray()
1206        expected = numpy.kron(self.Id, self.p)
1207        self.assertTrue(numpy.allclose(res, expected))
1208
1209    def test_boson_operator_sparse_trunc(self):
1210        op = BosonOperator('0')
1211        with self.assertRaises(ValueError):
1212            boson_operator_sparse(op, 0.1)
1213
1214        with self.assertRaises(ValueError):
1215            boson_operator_sparse(op, -1)
1216
1217        with self.assertRaises(ValueError):
1218            boson_operator_sparse(op, 0)
1219
1220    def test_boson_operator_invalid_op(self):
1221        op = FermionOperator('0')
1222        with self.assertRaises(ValueError):
1223            boson_operator_sparse(op, self.d)
1224
1225    def test_boson_operator_sparse_empty(self):
1226        for op in (BosonOperator(), QuadOperator()):
1227            res = boson_operator_sparse(op, self.d)
1228            self.assertEqual(res, numpy.array([[0]]))
1229
1230    def test_boson_operator_sparse_identity(self):
1231        for op in (BosonOperator(''), QuadOperator('')):
1232            res = boson_operator_sparse(op, self.d)
1233            self.assertEqual(res, numpy.array([[1]]))
1234
1235    def test_boson_operator_sparse_single(self):
1236        op = BosonOperator('0')
1237        res = boson_operator_sparse(op, self.d).toarray()
1238        self.assertTrue(numpy.allclose(res, self.b))
1239
1240        op = BosonOperator('0^')
1241        res = boson_operator_sparse(op, self.d).toarray()
1242        self.assertTrue(numpy.allclose(res, self.bd))
1243
1244        op = QuadOperator('q0')
1245        res = boson_operator_sparse(op, self.d, self.hbar).toarray()
1246        self.assertTrue(numpy.allclose(res, self.q))
1247
1248        op = QuadOperator('p0')
1249        res = boson_operator_sparse(op, self.d, self.hbar).toarray()
1250        self.assertTrue(numpy.allclose(res, self.p))
1251
1252    def test_boson_operator_sparse_number(self):
1253        op = BosonOperator('0^ 0')
1254        res = boson_operator_sparse(op, self.d).toarray()
1255        self.assertTrue(numpy.allclose(res, numpy.dot(self.bd, self.b)))
1256
1257    def test_boson_operator_sparse_multi_mode(self):
1258        op = BosonOperator('0^ 1 1^ 2')
1259        res = boson_operator_sparse(op, self.d).toarray()
1260
1261        b0 = boson_ladder_sparse(3, 0, 0, self.d).toarray()
1262        b1 = boson_ladder_sparse(3, 1, 0, self.d).toarray()
1263        b2 = boson_ladder_sparse(3, 2, 0, self.d).toarray()
1264
1265        expected = multi_dot([b0.T, b1, b1.T, b2])
1266        self.assertTrue(numpy.allclose(res, expected))
1267
1268        op = QuadOperator('q0 p0 p1')
1269        res = boson_operator_sparse(op, self.d, self.hbar).toarray()
1270
1271        expected = numpy.identity(self.d**2)
1272        for term in op.terms:
1273            for i, j in term:
1274                expected = expected.dot(
1275                    single_quad_op_sparse(2, i, j, self.hbar, self.d).toarray())
1276        self.assertTrue(numpy.allclose(res, expected))
1277
1278    def test_boson_operator_sparse_addition(self):
1279        op = BosonOperator('0^ 1')
1280        op += BosonOperator('0 0^')
1281        res = boson_operator_sparse(op, self.d).toarray()
1282
1283        b0 = boson_ladder_sparse(2, 0, 0, self.d).toarray()
1284        b1 = boson_ladder_sparse(2, 1, 0, self.d).toarray()
1285
1286        expected = numpy.dot(b0.T, b1) + numpy.dot(b0, b0.T)
1287        self.assertTrue(numpy.allclose(res, expected))
1288
1289
1290class GetNumberPreservingSparseOperatorIntegrationTestLiH(unittest.TestCase):
1291
1292    def setUp(self):
1293        # Set up molecule.
1294        geometry = [('Li', (0., 0., 0.)), ('H', (0., 0., 1.45))]
1295        basis = 'sto-3g'
1296        multiplicity = 1
1297        filename = os.path.join(DATA_DIRECTORY, 'H1-Li1_sto-3g_singlet_1.45')
1298        self.molecule = MolecularData(geometry,
1299                                      basis,
1300                                      multiplicity,
1301                                      filename=filename)
1302        self.molecule.load()
1303
1304        # Get molecular Hamiltonian.
1305        self.molecular_hamiltonian = self.molecule.get_molecular_hamiltonian()
1306
1307        self.hubbard_hamiltonian = fermi_hubbard(2,
1308                                                 2,
1309                                                 1.0,
1310                                                 4.0,
1311                                                 chemical_potential=.2,
1312                                                 magnetic_field=0.0,
1313                                                 spinless=False)
1314
1315    def test_exceptions(self):
1316        op = FermionOperator('1')
1317        with self.assertRaises(ValueError):
1318            _ = get_number_preserving_sparse_operator(op, 2, 1)
1319
1320    def test_number_on_reference(self):
1321        sum_n_op = FermionOperator()
1322        sum_sparse_n_op = get_number_preserving_sparse_operator(
1323            sum_n_op,
1324            self.molecule.n_qubits,
1325            self.molecule.n_electrons,
1326            spin_preserving=False)
1327
1328        space_size = sum_sparse_n_op.shape[0]
1329        reference = numpy.zeros((space_size))
1330        reference[0] = 1.0
1331
1332        for i in range(self.molecule.n_qubits):
1333            n_op = FermionOperator(((i, 1), (i, 0)))
1334            sum_n_op += n_op
1335
1336            sparse_n_op = get_number_preserving_sparse_operator(
1337                n_op,
1338                self.molecule.n_qubits,
1339                self.molecule.n_electrons,
1340                spin_preserving=False)
1341
1342            sum_sparse_n_op += sparse_n_op
1343
1344            expectation = reference.dot(sparse_n_op.dot(reference))
1345
1346            if i < self.molecule.n_electrons:
1347                assert expectation == 1.0
1348            else:
1349                assert expectation == 0.0
1350
1351        convert_after_adding = get_number_preserving_sparse_operator(
1352            sum_n_op,
1353            self.molecule.n_qubits,
1354            self.molecule.n_electrons,
1355            spin_preserving=False)
1356
1357        assert scipy.sparse.linalg.norm(convert_after_adding -
1358                                        sum_sparse_n_op) < 1E-9
1359
1360        assert reference.dot(sum_sparse_n_op.dot(reference)) - \
1361            self.molecule.n_electrons < 1E-9
1362
1363    def test_space_size_correct(self):
1364        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
1365
1366        sparse_ham = get_number_preserving_sparse_operator(
1367            hamiltonian_fop,
1368            self.molecule.n_qubits,
1369            self.molecule.n_electrons,
1370            spin_preserving=True)
1371
1372        space_size = sparse_ham.shape[0]
1373
1374        # Naive Hilbert space size is 2**12, or 4096.
1375        assert space_size == 225
1376
1377    def test_hf_energy(self):
1378        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
1379
1380        sparse_ham = get_number_preserving_sparse_operator(
1381            hamiltonian_fop,
1382            self.molecule.n_qubits,
1383            self.molecule.n_electrons,
1384            spin_preserving=True)
1385
1386        space_size = sparse_ham.shape[0]
1387        reference = numpy.zeros((space_size))
1388        reference[0] = 1.0
1389
1390        sparse_hf_energy = reference.dot(sparse_ham.dot(reference))
1391
1392        assert numpy.linalg.norm(sparse_hf_energy -
1393                                 self.molecule.hf_energy) < 1E-9
1394
1395    def test_one_body_hf_energy(self):
1396        one_body_part = self.molecular_hamiltonian
1397        one_body_part.two_body_tensor = numpy.zeros_like(
1398            one_body_part.two_body_tensor)
1399
1400        one_body_fop = get_fermion_operator(one_body_part)
1401        one_body_regular_sparse_op = get_sparse_operator(one_body_fop)
1402
1403        make_hf_fop = FermionOperator(((3, 1), (2, 1), (1, 1), (0, 1)))
1404        make_hf_sparse_op = get_sparse_operator(make_hf_fop, n_qubits=12)
1405
1406        hf_state = numpy.zeros((2**12))
1407        hf_state[0] = 1.0
1408        hf_state = make_hf_sparse_op.dot(hf_state)
1409
1410        regular_sparse_hf_energy = \
1411            (hf_state.dot(one_body_regular_sparse_op.dot(hf_state))).real
1412
1413        one_body_sparse_op = get_number_preserving_sparse_operator(
1414            one_body_fop,
1415            self.molecule.n_qubits,
1416            self.molecule.n_electrons,
1417            spin_preserving=True)
1418
1419        space_size = one_body_sparse_op.shape[0]
1420        reference = numpy.zeros((space_size))
1421        reference[0] = 1.0
1422
1423        sparse_hf_energy = reference.dot(one_body_sparse_op.dot(reference))
1424
1425        assert numpy.linalg.norm(sparse_hf_energy -
1426                                 regular_sparse_hf_energy) < 1E-9
1427
1428    def test_ground_state_energy(self):
1429        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
1430
1431        sparse_ham = get_number_preserving_sparse_operator(
1432            hamiltonian_fop,
1433            self.molecule.n_qubits,
1434            self.molecule.n_electrons,
1435            spin_preserving=True)
1436
1437        eig_val, _ = scipy.sparse.linalg.eigsh(sparse_ham, k=1, which='SA')
1438
1439        assert numpy.abs(eig_val[0] - self.molecule.fci_energy) < 1E-9
1440
1441    def test_doubles_are_subset(self):
1442        reference_determinants = [[
1443            True, True, True, True, False, False, False, False, False, False,
1444            False, False
1445        ],
1446                                  [
1447                                      True, True, False, False, False, False,
1448                                      True, True, False, False, False, False
1449                                  ]]
1450
1451        for reference_determinant in reference_determinants:
1452            reference_determinant = numpy.asarray(reference_determinant)
1453            doubles_state_array = numpy.asarray(
1454                list(
1455                    _iterate_basis_(reference_determinant,
1456                                    excitation_level=2,
1457                                    spin_preserving=True)))
1458            doubles_int_state_array = doubles_state_array.dot(
1459                1 << numpy.arange(doubles_state_array.shape[1])[::-1])
1460
1461            all_state_array = numpy.asarray(
1462                list(
1463                    _iterate_basis_(reference_determinant,
1464                                    excitation_level=4,
1465                                    spin_preserving=True)))
1466            all_int_state_array = all_state_array.dot(
1467                1 << numpy.arange(all_state_array.shape[1])[::-1])
1468
1469            for item in doubles_int_state_array:
1470                assert item in all_int_state_array
1471
1472        for reference_determinant in reference_determinants:
1473            reference_determinant = numpy.asarray(reference_determinant)
1474            doubles_state_array = numpy.asarray(
1475                list(
1476                    _iterate_basis_(reference_determinant,
1477                                    excitation_level=2,
1478                                    spin_preserving=True)))
1479            doubles_int_state_array = doubles_state_array.dot(
1480                1 << numpy.arange(doubles_state_array.shape[1])[::-1])
1481
1482            all_state_array = numpy.asarray(
1483                list(
1484                    _iterate_basis_(reference_determinant,
1485                                    excitation_level=4,
1486                                    spin_preserving=False)))
1487            all_int_state_array = all_state_array.dot(
1488                1 << numpy.arange(all_state_array.shape[1])[::-1])
1489
1490            for item in doubles_int_state_array:
1491                assert item in all_int_state_array
1492
1493        for reference_determinant in reference_determinants:
1494            reference_determinant = numpy.asarray(reference_determinant)
1495            doubles_state_array = numpy.asarray(
1496                list(
1497                    _iterate_basis_(reference_determinant,
1498                                    excitation_level=2,
1499                                    spin_preserving=False)))
1500            doubles_int_state_array = doubles_state_array.dot(
1501                1 << numpy.arange(doubles_state_array.shape[1])[::-1])
1502
1503            all_state_array = numpy.asarray(
1504                list(
1505                    _iterate_basis_(reference_determinant,
1506                                    excitation_level=4,
1507                                    spin_preserving=False)))
1508            all_int_state_array = all_state_array.dot(
1509                1 << numpy.arange(all_state_array.shape[1])[::-1])
1510
1511            for item in doubles_int_state_array:
1512                assert item in all_int_state_array
1513
1514    def test_full_ham_hermitian(self):
1515        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
1516
1517        sparse_ham = get_number_preserving_sparse_operator(
1518            hamiltonian_fop,
1519            self.molecule.n_qubits,
1520            self.molecule.n_electrons,
1521            spin_preserving=True)
1522
1523        assert scipy.sparse.linalg.norm(sparse_ham - sparse_ham.getH()) < 1E-9
1524
1525    def test_full_ham_hermitian_non_spin_preserving(self):
1526        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
1527
1528        sparse_ham = get_number_preserving_sparse_operator(
1529            hamiltonian_fop,
1530            self.molecule.n_qubits,
1531            self.molecule.n_electrons,
1532            spin_preserving=False)
1533
1534        assert scipy.sparse.linalg.norm(sparse_ham - sparse_ham.getH()) < 1E-9
1535
1536    def test_singles_simple_one_body_term_hermitian(self):
1537        fop = FermionOperator(((3, 1), (1, 0)))
1538        fop_conj = FermionOperator(((1, 1), (3, 0)))
1539
1540        sparse_op = get_number_preserving_sparse_operator(
1541            fop,
1542            self.molecule.n_qubits,
1543            self.molecule.n_electrons,
1544            spin_preserving=True,
1545            excitation_level=1)
1546
1547        sparse_op_conj = get_number_preserving_sparse_operator(
1548            fop_conj,
1549            self.molecule.n_qubits,
1550            self.molecule.n_electrons,
1551            spin_preserving=True,
1552            excitation_level=1)
1553
1554        assert scipy.sparse.linalg.norm(sparse_op -
1555                                        sparse_op_conj.getH()) < 1E-9
1556
1557    def test_singles_simple_two_body_term_hermitian(self):
1558        fop = FermionOperator(((3, 1), (8, 1), (1, 0), (4, 0)))
1559        fop_conj = FermionOperator(((4, 1), (1, 1), (8, 0), (3, 0)))
1560
1561        sparse_op = get_number_preserving_sparse_operator(
1562            fop,
1563            self.molecule.n_qubits,
1564            self.molecule.n_electrons,
1565            spin_preserving=True,
1566            excitation_level=1)
1567
1568        sparse_op_conj = get_number_preserving_sparse_operator(
1569            fop_conj,
1570            self.molecule.n_qubits,
1571            self.molecule.n_electrons,
1572            spin_preserving=True,
1573            excitation_level=1)
1574
1575        assert scipy.sparse.linalg.norm(sparse_op -
1576                                        sparse_op_conj.getH()) < 1E-9
1577
1578    def test_singles_repeating_two_body_term_hermitian(self):
1579        fop = FermionOperator(((3, 1), (1, 1), (5, 0), (1, 0)))
1580        fop_conj = FermionOperator(((5, 1), (1, 1), (3, 0), (1, 0)))
1581
1582        sparse_op = get_number_preserving_sparse_operator(
1583            fop,
1584            self.molecule.n_qubits,
1585            self.molecule.n_electrons,
1586            spin_preserving=True,
1587            excitation_level=1)
1588
1589        sparse_op_conj = get_number_preserving_sparse_operator(
1590            fop_conj,
1591            self.molecule.n_qubits,
1592            self.molecule.n_electrons,
1593            spin_preserving=True,
1594            excitation_level=1)
1595
1596        assert scipy.sparse.linalg.norm(sparse_op -
1597                                        sparse_op_conj.getH()) < 1E-9
1598
1599    def test_singles_ham_hermitian(self):
1600        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
1601
1602        sparse_ham = get_number_preserving_sparse_operator(
1603            hamiltonian_fop,
1604            self.molecule.n_qubits,
1605            self.molecule.n_electrons,
1606            spin_preserving=True,
1607            excitation_level=1)
1608
1609        assert scipy.sparse.linalg.norm(sparse_ham - sparse_ham.getH()) < 1E-9
1610
1611    def test_singles_ham_hermitian_non_spin_preserving(self):
1612        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
1613
1614        sparse_ham = get_number_preserving_sparse_operator(
1615            hamiltonian_fop,
1616            self.molecule.n_qubits,
1617            self.molecule.n_electrons,
1618            spin_preserving=False,
1619            excitation_level=1)
1620
1621        assert scipy.sparse.linalg.norm(sparse_ham - sparse_ham.getH()) < 1E-9
1622
1623    def test_cisd_energy(self):
1624        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
1625
1626        sparse_ham = get_number_preserving_sparse_operator(
1627            hamiltonian_fop,
1628            self.molecule.n_qubits,
1629            self.molecule.n_electrons,
1630            spin_preserving=True,
1631            excitation_level=2)
1632
1633        eig_val, _ = scipy.sparse.linalg.eigsh(sparse_ham, k=1, which='SA')
1634
1635        assert numpy.abs(eig_val[0] - self.molecule.cisd_energy) < 1E-9
1636
1637    def test_cisd_energy_non_spin_preserving(self):
1638        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
1639
1640        sparse_ham = get_number_preserving_sparse_operator(
1641            hamiltonian_fop,
1642            self.molecule.n_qubits,
1643            self.molecule.n_electrons,
1644            spin_preserving=False,
1645            excitation_level=2)
1646
1647        eig_val, _ = scipy.sparse.linalg.eigsh(sparse_ham, k=1, which='SA')
1648
1649        assert numpy.abs(eig_val[0] - self.molecule.cisd_energy) < 1E-9
1650
1651    def test_cisd_matches_fci_energy_two_electron_hubbard(self):
1652        hamiltonian_fop = self.hubbard_hamiltonian
1653
1654        sparse_ham_cisd = get_number_preserving_sparse_operator(
1655            hamiltonian_fop, 8, 2, spin_preserving=True, excitation_level=2)
1656
1657        sparse_ham_fci = get_sparse_operator(hamiltonian_fop, n_qubits=8)
1658
1659        eig_val_cisd, _ = scipy.sparse.linalg.eigsh(sparse_ham_cisd,
1660                                                    k=1,
1661                                                    which='SA')
1662        eig_val_fci, _ = scipy.sparse.linalg.eigsh(sparse_ham_fci,
1663                                                   k=1,
1664                                                   which='SA')
1665
1666        assert numpy.abs(eig_val_cisd[0] - eig_val_fci[0]) < 1E-9
1667
1668    def test_weird_determinant_matches_fci_energy_two_electron_hubbard(self):
1669        hamiltonian_fop = self.hubbard_hamiltonian
1670
1671        sparse_ham_cisd = get_number_preserving_sparse_operator(
1672            hamiltonian_fop,
1673            8,
1674            2,
1675            spin_preserving=True,
1676            excitation_level=2,
1677            reference_determinant=numpy.asarray(
1678                [False, False, True, True, False, False, False, False]))
1679
1680        sparse_ham_fci = get_sparse_operator(hamiltonian_fop, n_qubits=8)
1681
1682        eig_val_cisd, _ = scipy.sparse.linalg.eigsh(sparse_ham_cisd,
1683                                                    k=1,
1684                                                    which='SA')
1685        eig_val_fci, _ = scipy.sparse.linalg.eigsh(sparse_ham_fci,
1686                                                   k=1,
1687                                                   which='SA')
1688
1689        assert numpy.abs(eig_val_cisd[0] - eig_val_fci[0]) < 1E-9
1690
1691    def test_number_restricted_spectra_match_molecule(self):
1692        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
1693
1694        sparse_ham_number_preserving = get_number_preserving_sparse_operator(
1695            hamiltonian_fop,
1696            self.molecule.n_qubits,
1697            self.molecule.n_electrons,
1698            spin_preserving=False)
1699
1700        sparse_ham = get_sparse_operator(hamiltonian_fop,
1701                                         self.molecule.n_qubits)
1702
1703        sparse_ham_restricted_number_preserving = jw_number_restrict_operator(
1704            sparse_ham,
1705            n_electrons=self.molecule.n_electrons,
1706            n_qubits=self.molecule.n_qubits)
1707
1708        spectrum_from_new_sparse_method = sparse_eigenspectrum(
1709            sparse_ham_number_preserving)
1710
1711        spectrum_from_old_sparse_method = sparse_eigenspectrum(
1712            sparse_ham_restricted_number_preserving)
1713
1714        spectral_deviation = numpy.amax(
1715            numpy.absolute(spectrum_from_new_sparse_method -
1716                           spectrum_from_old_sparse_method))
1717        self.assertAlmostEqual(spectral_deviation, 0.)
1718
1719    def test_number_restricted_spectra_match_hubbard(self):
1720        hamiltonian_fop = self.hubbard_hamiltonian
1721
1722        sparse_ham_number_preserving = get_number_preserving_sparse_operator(
1723            hamiltonian_fop, 8, 4, spin_preserving=False)
1724
1725        sparse_ham = get_sparse_operator(hamiltonian_fop, 8)
1726
1727        sparse_ham_restricted_number_preserving = jw_number_restrict_operator(
1728            sparse_ham, n_electrons=4, n_qubits=8)
1729
1730        spectrum_from_new_sparse_method = sparse_eigenspectrum(
1731            sparse_ham_number_preserving)
1732
1733        spectrum_from_old_sparse_method = sparse_eigenspectrum(
1734            sparse_ham_restricted_number_preserving)
1735
1736        spectral_deviation = numpy.amax(
1737            numpy.absolute(spectrum_from_new_sparse_method -
1738                           spectrum_from_old_sparse_method))
1739        self.assertAlmostEqual(spectral_deviation, 0.)
1740
1741    def test_number_sz_restricted_spectra_match_molecule(self):
1742        hamiltonian_fop = get_fermion_operator(self.molecular_hamiltonian)
1743
1744        sparse_ham_number_sz_preserving = get_number_preserving_sparse_operator(
1745            hamiltonian_fop,
1746            self.molecule.n_qubits,
1747            self.molecule.n_electrons,
1748            spin_preserving=True)
1749
1750        sparse_ham = get_sparse_operator(hamiltonian_fop,
1751                                         self.molecule.n_qubits)
1752
1753        sparse_ham_restricted_number_sz_preserving = jw_sz_restrict_operator(
1754            sparse_ham,
1755            0,
1756            n_electrons=self.molecule.n_electrons,
1757            n_qubits=self.molecule.n_qubits)
1758
1759        spectrum_from_new_sparse_method = sparse_eigenspectrum(
1760            sparse_ham_number_sz_preserving)
1761
1762        spectrum_from_old_sparse_method = sparse_eigenspectrum(
1763            sparse_ham_restricted_number_sz_preserving)
1764
1765        spectral_deviation = numpy.amax(
1766            numpy.absolute(spectrum_from_new_sparse_method -
1767                           spectrum_from_old_sparse_method))
1768        self.assertAlmostEqual(spectral_deviation, 0.)
1769
1770    def test_number_sz_restricted_spectra_match_hubbard(self):
1771        hamiltonian_fop = self.hubbard_hamiltonian
1772
1773        sparse_ham_number_sz_preserving = get_number_preserving_sparse_operator(
1774            hamiltonian_fop, 8, 4, spin_preserving=True)
1775
1776        sparse_ham = get_sparse_operator(hamiltonian_fop, 8)
1777
1778        sparse_ham_restricted_number_sz_preserving = jw_sz_restrict_operator(
1779            sparse_ham, 0, n_electrons=4, n_qubits=8)
1780
1781        spectrum_from_new_sparse_method = sparse_eigenspectrum(
1782            sparse_ham_number_sz_preserving)
1783
1784        spectrum_from_old_sparse_method = sparse_eigenspectrum(
1785            sparse_ham_restricted_number_sz_preserving)
1786
1787        spectral_deviation = numpy.amax(
1788            numpy.absolute(spectrum_from_new_sparse_method -
1789                           spectrum_from_old_sparse_method))
1790        self.assertAlmostEqual(spectral_deviation, 0.)
1791
1792
1793class GetSparseOperatorQubitTest(unittest.TestCase):
1794
1795    def test_sparse_matrix_Y(self):
1796        term = QubitOperator(((0, 'Y'),))
1797        sparse_operator = get_sparse_operator(term)
1798        self.assertEqual(list(sparse_operator.data), [1j, -1j])
1799        self.assertEqual(list(sparse_operator.indices), [1, 0])
1800        self.assertTrue(is_hermitian(sparse_operator))
1801
1802    def test_sparse_matrix_ZX(self):
1803        coefficient = 2.
1804        operators = ((0, 'Z'), (1, 'X'))
1805        term = QubitOperator(operators, coefficient)
1806        sparse_operator = get_sparse_operator(term)
1807        self.assertEqual(list(sparse_operator.data), [2., 2., -2., -2.])
1808        self.assertEqual(list(sparse_operator.indices), [1, 0, 3, 2])
1809        self.assertTrue(is_hermitian(sparse_operator))
1810
1811    def test_sparse_matrix_ZIZ(self):
1812        operators = ((0, 'Z'), (2, 'Z'))
1813        term = QubitOperator(operators)
1814        sparse_operator = get_sparse_operator(term)
1815        self.assertEqual(list(sparse_operator.data),
1816                         [1, -1, 1, -1, -1, 1, -1, 1])
1817        self.assertEqual(list(sparse_operator.indices), list(range(8)))
1818        self.assertTrue(is_hermitian(sparse_operator))
1819
1820    def test_sparse_matrix_combo(self):
1821        qop = (QubitOperator(((0, 'Y'), (1, 'X')), -0.1j) + QubitOperator(
1822            ((0, 'X'), (1, 'Z')), 3. + 2.j))
1823        sparse_operator = get_sparse_operator(qop)
1824
1825        self.assertEqual(
1826            list(sparse_operator.data),
1827            [3 + 2j, 0.1, 0.1, -3 - 2j, 3 + 2j, -0.1, -0.1, -3 - 2j])
1828        self.assertEqual(list(sparse_operator.indices),
1829                         [2, 3, 2, 3, 0, 1, 0, 1])
1830
1831    def test_sparse_matrix_zero_1qubit(self):
1832        sparse_operator = get_sparse_operator(QubitOperator((), 0.0), 1)
1833        sparse_operator.eliminate_zeros()
1834        self.assertEqual(len(list(sparse_operator.data)), 0)
1835        self.assertEqual(sparse_operator.shape, (2, 2))
1836
1837    def test_sparse_matrix_zero_5qubit(self):
1838        sparse_operator = get_sparse_operator(QubitOperator((), 0.0), 5)
1839        sparse_operator.eliminate_zeros()
1840        self.assertEqual(len(list(sparse_operator.data)), 0)
1841        self.assertEqual(sparse_operator.shape, (32, 32))
1842
1843    def test_sparse_matrix_identity_1qubit(self):
1844        sparse_operator = get_sparse_operator(QubitOperator(()), 1)
1845        self.assertEqual(list(sparse_operator.data), [1] * 2)
1846        self.assertEqual(sparse_operator.shape, (2, 2))
1847
1848    def test_sparse_matrix_identity_5qubit(self):
1849        sparse_operator = get_sparse_operator(QubitOperator(()), 5)
1850        self.assertEqual(list(sparse_operator.data), [1] * 32)
1851        self.assertEqual(sparse_operator.shape, (32, 32))
1852
1853    def test_sparse_matrix_linearity(self):
1854        identity = QubitOperator(())
1855        zzzz = QubitOperator(tuple((i, 'Z') for i in range(4)), 1.0)
1856
1857        sparse1 = get_sparse_operator(identity + zzzz)
1858        sparse2 = get_sparse_operator(identity, 4) + get_sparse_operator(zzzz)
1859
1860        self.assertEqual(list(sparse1.data), [2] * 8)
1861        self.assertEqual(list(sparse1.indices), [0, 3, 5, 6, 9, 10, 12, 15])
1862        self.assertEqual(list(sparse2.data), [2] * 8)
1863        self.assertEqual(list(sparse2.indices), [0, 3, 5, 6, 9, 10, 12, 15])
1864
1865
1866class GetSparseOperatorFermionTest(unittest.TestCase):
1867
1868    def test_sparse_matrix_zero_n_qubit(self):
1869        sparse_operator = get_sparse_operator(FermionOperator.zero(), 4)
1870        sparse_operator.eliminate_zeros()
1871        self.assertEqual(len(list(sparse_operator.data)), 0)
1872        self.assertEqual(sparse_operator.shape, (16, 16))
1873
1874
1875class GetSparseOperatorBosonTest(unittest.TestCase):
1876
1877    def setUp(self):
1878        self.hbar = 1.
1879        self.d = 4
1880        self.b = numpy.diag(numpy.sqrt(numpy.arange(1, self.d)), 1)
1881        self.bd = self.b.conj().T
1882        self.q = numpy.sqrt(self.hbar / 2) * (self.b + self.bd)
1883
1884    def test_sparse_matrix_ladder(self):
1885        sparse_operator = get_sparse_operator(BosonOperator('0'), trunc=self.d)
1886        self.assertTrue(numpy.allclose(sparse_operator.toarray(), self.b))
1887        self.assertEqual(sparse_operator.shape, (self.d, self.d))
1888
1889    def test_sparse_matrix_quad(self):
1890        sparse_operator = get_sparse_operator(QuadOperator('q0'), trunc=self.d)
1891        self.assertTrue(numpy.allclose(sparse_operator.toarray(), self.q))
1892        self.assertEqual(sparse_operator.shape, (self.d, self.d))
1893
1894    def test_sparse_matrix_error(self):
1895        with self.assertRaises(TypeError):
1896            _ = get_sparse_operator(1)
1897
1898
1899class GetSparseOperatorDiagonalCoulombHamiltonianTest(unittest.TestCase):
1900
1901    def test_diagonal_coulomb_hamiltonian(self):
1902        n_qubits = 5
1903        one_body = random_hermitian_matrix(n_qubits, real=False)
1904        two_body = random_hermitian_matrix(n_qubits, real=True)
1905        constant = numpy.random.randn()
1906        op = DiagonalCoulombHamiltonian(one_body, two_body, constant)
1907
1908        op1 = get_sparse_operator(op)
1909        op2 = get_sparse_operator(jordan_wigner(get_fermion_operator(op)))
1910        diff = op1 - op2
1911        discrepancy = 0.
1912        if diff.nnz:
1913            discrepancy = max(abs(diff.data))
1914        self.assertAlmostEqual(discrepancy, 0.)
1915