Improve test coverage

src/fqe/bitstring.py

@@ -22,9 +22,10 @@
 from numpy import ndarray as Nparray
 from scipy import special
 
+import fqe
 from fqe.lib.bitstring import _lexicographic_bitstring_generator, _count_bits, \
                               _get_occupation
-from fqe.settings import use_accelerated_code, c_string_max_norb
+from fqe.settings import c_string_max_norb
 
 
 def gbit_index(str0: int) -> Generator[int, None, None]:
@@ -68,7 +69,7 @@ def integer_index(string: int) -> 'Nparray':
         Nparray: a list of integers indicating the index of each occupied \
             orbital
     """
-    if use_accelerated_code:
+    if fqe.settings.use_accelerated_code:
         return _get_occupation(string)
     else:
         return numpy.array(list(gbit_index(int(string)))).astype(numpy.int32)
@@ -107,7 +108,7 @@ def lexicographic_bitstring_generator(nele: int, norb: int) -> 'Nparray':
     if nele > norb:
         raise ValueError("nele cannot be larger than norb")
 
-    if use_accelerated_code and norb <= c_string_max_norb:
+    if fqe.settings.use_accelerated_code and norb <= c_string_max_norb:
         out = numpy.zeros((int(special.comb(norb, nele)),), dtype=numpy.uint64)
         _lexicographic_bitstring_generator(out, norb, nele)
         return out
@@ -127,7 +128,7 @@ def count_bits(string: int) -> int:
     Returns:
         int: the number of bits equal to 1
     """
-    if use_accelerated_code:
+    if fqe.settings.use_accelerated_code:
         return _count_bits(string)
     else:
         return bin(int(string)).count('1')

src/fqe/fci_graph.py

@@ -22,13 +22,14 @@
 import numpy
 from numpy import ndarray as Nparray
 
+import fqe
 from fqe.bitstring import integer_index, lexicographic_bitstring_generator, \
     get_bit, set_bit, unset_bit, count_bits_above, count_bits_between
 from fqe.lib.fci_graph import _build_mapping_strings, _map_deexc, \
                               _calculate_string_address, \
                               _c_map_to_deexc_alpha_icol, \
                               _make_mapping_each, _calculate_Z_matrix
-from fqe.settings import use_accelerated_code, c_string_max_norb
+from fqe.settings import c_string_max_norb
 
 Spinmap = Dict[Tuple[int, ...], Nparray]
 
@@ -54,7 +55,7 @@ def map_to_deexc(mappings: Spinmap, states: int, norbs: int,
     index = numpy.zeros((states,), dtype=numpy.uint32)
     for (i, j), values in mappings.items():
         idx = i * norbs + j
-        if use_accelerated_code and norbs <= c_string_max_norb:
+        if fqe.settings.use_accelerated_code and norbs <= c_string_max_norb:
             _map_deexc(dexc, values, index, idx)
         else:
             for state, target, parity in values:
@@ -82,7 +83,7 @@ def _get_Z_matrix(norb: int, nele: int) -> Nparray:
     if Z.size == 0:
         return Z
 
-    if use_accelerated_code and norb <= c_string_max_norb:
+    if fqe.settings.use_accelerated_code and norb <= c_string_max_norb:
         _calculate_Z_matrix(Z, norb, nele)
     else:
         for k in range(1, nele):
@@ -217,7 +218,7 @@ def _build_mapping(self, strings: Nparray, nele: int,
         """
         norb = self._norb
 
-        if use_accelerated_code and norb <= c_string_max_norb:
+        if fqe.settings.use_accelerated_code and norb <= c_string_max_norb:
             return _build_mapping_strings(strings, _get_Z_matrix(norb, nele),
                                           nele, norb)
         else:
@@ -315,7 +316,7 @@ def _build_strings(self, nele: int,
         norb = self._norb
         blist = lexicographic_bitstring_generator(nele, norb)
 
-        if use_accelerated_code and norb <= c_string_max_norb:
+        if fqe.settings.use_accelerated_code and norb <= c_string_max_norb:
             Z = _get_Z_matrix(norb, nele)
             string_list = _calculate_string_address(Z, nele, norb, blist)
         else:
@@ -470,7 +471,7 @@ def _map_to_deexc_alpha_icol(self):
             length2,
         ), dtype=numpy.int32)
         astrings = self.string_alpha_all()
-        if use_accelerated_code and norb <= c_string_max_norb:
+        if fqe.settings.use_accelerated_code and norb <= c_string_max_norb:
             _c_map_to_deexc_alpha_icol(exc, diag, index, astrings, norb,
                                        self._alpha_map)
         else:
@@ -538,7 +539,7 @@ def make_mapping_each(self, result: 'Nparray', alpha: bool, dag: List[int],
             strings = self.string_beta_all()
             length = self.lenb()
 
-        if use_accelerated_code:
+        if fqe.settings.use_accelerated_code:
             count = _make_mapping_each(result, strings, length,
                                        numpy.array(dag, dtype=numpy.int32),
                                        numpy.array(undag, dtype=numpy.int32))

src/fqe/fci_graph_set.py

@@ -23,7 +23,8 @@
 import scipy
 from scipy import special
 
-from fqe.settings import use_accelerated_code, c_string_max_norb
+import fqe
+from fqe.settings import c_string_max_norb
 from fqe.fci_graph import FciGraph
 from fqe.util import alpha_beta_electrons
 from fqe.bitstring import integer_index, count_bits_between
@@ -171,7 +172,7 @@ def _postprocess(spinmap, dnv, index0, index1):
 
         if dna != 0:
             (iasec, jasec) = (isec, jsec) if dna < 0 else (jsec, isec)
-            if use_accelerated_code and norb <= c_string_max_norb:
+            if fqe.settings.use_accelerated_code and norb <= c_string_max_norb:
                 ndowna, nupa = _make_mapping_each_set(iasec.string_alpha_all(),
                                                       abs(dna), norb,
                                                       iasec.nalpha())
@@ -192,7 +193,7 @@ def _postprocess(spinmap, dnv, index0, index1):
 
         if dnb != 0:
             (ibsec, jbsec) = (isec, jsec) if dnb < 0 else (jsec, isec)
-            if use_accelerated_code and norb <= c_string_max_norb:
+            if fqe.settings.use_accelerated_code and norb <= c_string_max_norb:
                 ndownb, nupb = _make_mapping_each_set(ibsec.string_beta_all(),
                                                       abs(dnb), norb,
                                                       ibsec.nbeta())

src/fqe/wavefunction.py

@@ -28,7 +28,7 @@
 import pickle
 import numpy
 from scipy import linalg
-from scipy.special import factorial, jv
+from scipy.special import jv
 
 from fqe.fqe_decorators import wrap_apply, wrap_apply_generated_unitary
 from fqe.fqe_decorators import wrap_time_evolve, wrap_rdm
@@ -555,12 +555,12 @@ def apply_generated_unitary(self,
 
         if algo == 'taylor':
             ham_arrays = hamil.iht(time)
-
             time_evol = copy.deepcopy(base)
             work = copy.deepcopy(base)
+            coeff = 1.0
             for order in range(1, max_expansion):
                 work = work.apply(ham_arrays)
-                coeff = 1.0 / factorial(order)
+                coeff /= order
                 time_evol.ax_plus_y(coeff, work)
                 if work.norm() * numpy.abs(coeff) < accuracy:
                     break
@@ -1231,7 +1231,7 @@ def _evolve_individual_nbody(self,
                     raise ValueError(
                         'Operators in _evolve_individual_nbody is not Hermitian'
                     )
-        else:
+        elif hamil.nterms() == 1:
             [
                 (coeff0, alpha0, beta0),
             ] = hamil.terms()
@@ -1247,6 +1247,9 @@ def _evolve_individual_nbody(self,
                         'Operators in _evolve_individual_nbody is not Hermitian'
                     )
             coeff0 *= 0.5
+        else:
+            assert hamil.nterms() == 0
+            return self if inplace else copy.deepcopy(self)
 
         daga = []
         dagb = []

tests/adapt_vqe_test.py

@@ -101,4 +101,4 @@ def test_adapt():
                   verbose=True)
     adapt.adapt_vqe(fqe_wf)
     assert np.isclose(adapt.energies[-1], -8.957417182801091)
-    assert np.isclose(adapt.energies[0], -8.95741717733075)
+    assert np.isclose(adapt.energies[0], -8.95741717733075)

tests/bitstring_test.py

@@ -21,9 +21,10 @@
 import fqe.settings
 
 
-def test_bit_integer_index_val():
+def test_bit_integer_index_val(c_or_python):
     """The index of bits should start at 0.
     """
+    fqe.settings.use_accelerated_code = c_or_python
     _gbiti = bitstring.gbit_index(1)
     assert next(_gbiti) == 0
     _gbiti = bitstring.gbit_index(11)

tests/fci_graph_test.py

@@ -15,10 +15,14 @@
 """
 
 import numpy
-import fqe
 import pytest
 from scipy import special
+
+import fqe
 from fqe import fci_graph
+from fqe.lib.fci_graph import _c_map_to_deexc_alpha_icol
+from fqe.settings import CodePath
+
 from tests.unittest_data.fci_graph_data import loader
 from tests.comparisons import compare_Spinmap
 
@@ -311,6 +315,14 @@ def test_map_to_deexc_alpha_icol(c_or_python, norb, nalpha, nbeta):
     assert numpy.array_equal(rexc, exc)
     assert numpy.array_equal(rdiag, diag)
 
+    if c_or_python == CodePath.PYTHON:
+        return
+
+    tmp = numpy.zeros((norb, norb))
+    test = {(i, i): numpy.zeros((1, 1)) for i in range(norb)}
+    with pytest.raises(ValueError):
+        _c_map_to_deexc_alpha_icol(tmp, tmp, tmp, tmp, norb, test)
+
 
 @pytest.mark.parametrize("nalpha,nbeta,norb", cases)
 def test_get_block_mappings(norb, nalpha, nbeta):

tests/linalg_test.py

@@ -16,10 +16,22 @@
 import pytest
 import numpy
 import fqe
+from fqe.lib import c_double_complex
 from fqe.lib.linalg import _zimatadd, _transpose
 from fqe.settings import CodePath
 
 
+def test_ctype_double_complex(c_or_python):
+    if c_or_python == CodePath.PYTHON:
+        # No test needed for python
+        return
+
+    rval = 1.0
+    ival = -1.0
+    test = c_double_complex(rval, ival)
+    assert test.value == rval + 1.j * ival
+
+
 def test_zimatadd(c_or_python):
     """Testing zimatadd"""
     if c_or_python == CodePath.PYTHON:

tests/unittest_data/wavefunction/generate_data.py

@@ -65,6 +65,7 @@ def generate_data(param):
         cr = rng.uniform(-0.5, 0.5, size=value.coeff.shape)
         ci = rng.uniform(-0.5, 0.5, size=value.coeff.shape)
         value.coeff = cr + 1.j * ci
+    wfn.normalize()
 
     data['wfn'] = wfn
     if not spin_conserving:
@@ -115,8 +116,8 @@ def generate_data(param):
 
     # array ham
     full = not (spin_conserving and number_conserving)
-    h1e = build_H1(norbs, full=full)
-    h2e = build_H2(norbs, full=full)
+    h1e = build_H1(norbs, full=full) / 10
+    h2e = build_H2(norbs, full=full) / 20
     e_0 = rng.uniform() + rng.uniform() * 1j
     if not full:
         hamil = get_restricted_hamiltonian((
@@ -139,6 +140,7 @@ def generate_data(param):
     except ValueError:
         # Don't generate an output
         evolved = None
+    #except RuntimeError
 
     data['apply_array'] = {
         'hamil': hamil,

tests/wavefunction_test.py

@@ -15,18 +15,16 @@
 """
 # pylint: disable=protected-access
 
-import os
-import sys
 import copy
-import numpy
-import pytest
-
+import os
 from io import StringIO
+import sys
 
-from scipy.special import binom
-
+import numpy
 from openfermion import FermionOperator
 from openfermion.utils import hermitian_conjugated
+import pytest
+from scipy.special import binom
 
 import fqe
 from fqe.wavefunction import Wavefunction
@@ -116,7 +114,7 @@ def test_empty_copy():
     test1 = test.empty_copy()
     assert test1._norb == test._norb
     assert len(test1._civec) == len(test._civec)
-    assert (2, 0) in test1._civec.keys()
+    assert (2, 0) in test1._civec
     assert not numpy.any(test1._civec[(2, 0)].coeff)
 
 
@@ -160,6 +158,10 @@ def test_apply_type_error():
     with pytest.raises(TypeError):
         wfn.time_evolve(0.1, hamil)
 
+    hamil2 = restricted_hamiltonian.RestrictedHamiltonian((data,))
+    with pytest.raises(TypeError, match="expansion must be an int"):
+        wfn.apply_generated_unitary(0.1, 'taylor', hamil2, expansion=0.5)
+
 
 def test_apply_value_error():
     data = numpy.zeros((2, 2), dtype=numpy.complex128)
@@ -299,6 +301,60 @@ def test_apply_empty_nbody():
     assert (wfn - out1).norm() < 1e-13
 
 
+def test_expansion_failure():
+    """Test exceptions when the max expansion limit is reached.
+    """
+    norb = 4
+    nalpha = 2
+    nbeta = 1
+    nele = nalpha + nbeta
+    time = 0.1
+
+    h1e = build_hamiltonian.build_H1(norb, full=True, asymmetric=True)
+
+    eig, _ = numpy.linalg.eigh(h1e)
+    hamil = fqe.get_general_hamiltonian((h1e,))
+
+    wfn = Wavefunction([[nele, nalpha - nbeta, norb]])
+    wfn.set_wfn(strategy='random')
+
+    err = "maximum chebyshev expansion limit reached"
+    with pytest.raises(RuntimeError, match=err):
+        wfn.apply_generated_unitary(time,
+                                    'chebyshev',
+                                    hamil,
+                                    expansion=1,
+                                    accuracy=1e-9,
+                                    spec_lim=(eig[0], eig[-1]))
+
+    err = "maximum taylor expansion limit reached"
+    with pytest.raises(RuntimeError, match=err):
+        wfn.apply_generated_unitary(time,
+                                    'taylor',
+                                    hamil,
+                                    expansion=1,
+                                    accuracy=1e-9)
+
+
+def test_empty_nbody_evolve():
+    """Check that '_evolve_individual_nbody' works with an empty
+    FermionOperator/SparseHamiltonian
+    """
+    norb = 4
+    nele = 4
+    time = 0.1
+    ops = FermionOperator()
+    sham = fqe.get_sparse_hamiltonian(ops, conserve_spin=False)
+
+    wfn = fqe.get_number_conserving_wavefunction(nele, norb)
+    wfn.set_wfn(strategy='random')
+    wfn.normalize()
+
+    out = wfn._evolve_individual_nbody(time, sham)
+
+    assert (out - wfn).norm() < 1.e-8
+
+
 def test_nbody_evolve():
     """Check that 'individual_nbody' is consistent with a general 1-electron op
     """
@@ -626,7 +682,6 @@ def test_apply(c_or_python, param, kind):
         assert Wavefunction_isclose(out, reference_data['wfn_out'])
 
 
-@pytest.mark.skip(reason="Test seems to be failing due to bad reference_data")
 @pytest.mark.parametrize("param,kind", [(c, k) for c in all_cases for k in [
     'apply_array', 'apply_sparse', 'apply_diagonal', 'apply_quadratic',
     'apply_dc'