update documentation

h4_experiment.py

@@ -1,7 +1,9 @@
+from typing import Dict, List, Tuple
 import os
 import pickle
 from tqdm import tqdm
 import numpy as np
+import scipy
 from multiprocessing import Pool
 import datetime
 
@@ -17,8 +19,15 @@
 from sampler import local_dists_optimal
 
 
-def _convert_spin_sector(ham_f):
-    "convert |↑↓↑↓...> representation to |↑↑...↓↓...> representation"
+def _convert_spin_sector(ham_f: FermionOperator) -> FermionOperator:
+    """convert |↑↓↑↓...> representation to |↑↑...↓↓...> representation.
+
+    Args:
+        ham_f (FermionOperator): Second-quantized molecular Hamiltonian in |↑↓↑↓...> representation.
+
+    Returns:
+        FermionOperator: Second-quantized molecular Hamiltonian in |↑↑...↓↓...> representation.
+    """
     n = count_qubits(ham_f)
     n_half = int(n / 2)
     index_map = {i: int(i / 2) if i % 2 == 0 else n_half + int((i - 1) / 2) for i in range(n)}
@@ -27,7 +36,15 @@ def _convert_spin_sector(ham_f):
     return ham_convert
 
 
-def _get_h4_hamiltonian(length):
+def _get_h4_hamiltonian(length: float) -> QubitOperator:
+    """Generate H4 Hamiltonian from given interatomic distance.
+
+    Args:
+        length (float): Interatomic distance (Å).
+
+    Returns:
+        QubitOperator: H4 Hamiltonian.
+    """
     geom = [
         ("H", (-length / 2 * 3, 0, 0)),
         ("H", (-length / 2, 0, 0)),
@@ -52,7 +69,13 @@ def _get_h4_hamiltonian(length):
     return ham
 
 
-def simulate_energy_vs_interatomic_distance():
+def simulate_energy_vs_interatomic_distance() -> Tuple[float, float, float]:
+    """Calculate energies with various interatomic distances.
+
+    Returns:
+        Tuple[float, float, float]: Rigorous energies, and energies estimated through the ground state and CISD state.
+    """
+
     energies_rigorous = []
     energies_2n1p_gs = []
     energies_2n1p_cisd = []
@@ -79,7 +102,17 @@ def simulate_energy_vs_interatomic_distance():
     return energies_rigorous, energies_2n1p_gs, energies_2n1p_cisd
 
 
-def _estimate_qse_matrix_elements_noise_simulator(beta_eff, params):
+def _estimate_qse_matrix_elements_noise_simulator(beta_eff: np.ndarray, params: Dict) -> Tuple[np.ndarray, np.ndarray]:
+    """Estimate operator matrices of QSE through LBCS using noise simulator.
+
+    Args:
+        beta_eff (np.ndarray): Weighted bias of measurement bases within LBCS.
+        params (Dict): Parameters for execution.
+
+    Returns:
+        Tuple[np.ndarray, np.ndarray]: Estimated matrix elements  \tilde{H}_ij = <\psi| O_i^\dag H O_j |\psi> and \tilde{S}_ij = <\psi| O_i^\dag O_j |\psi>
+    """
+
     def estimate_qse_matrix_element_noise_simulator(i, j, op_mat, result_mat):
         if j >= i:
             for term, coef in op_mat[i, j].terms.items():
@@ -116,7 +149,15 @@ def get_n_shots_per_pauli(term, beta, shots_total):
     return h_eff, s_mtrc
 
 
-def _iterative_run_lbcs(params):
+def _iterative_run_lbcs(params: Dict) -> Dict:
+    """Apply adaptive measurement optimization (LBCS)
+
+    Args:
+        params (Dict): Parameters for execution.
+
+    Returns:
+        Dict: QSE results in each iteration
+    """
     subspace_expansion = params["subspace_expansion"]
     hamiltonian = params["molecule"].hamiltonian
     n_qubit = subspace_expansion.n_qubit
@@ -171,7 +212,20 @@ def _iterative_run_lbcs(params):
     return dict_iter_lbcs
 
 
-def _get_pauli_mat_dict(file_name, n_qubits, h_ij, load=True):
+def _get_pauli_mat_dict(
+    file_name: str, n_qubits: int, h_ij: Dict[Tuple, QubitOperator], load=True
+) -> Dict[Tuple, scipy.sparse.csc.csc_matrix]:
+    """Generate dictionary of sparse matrix.
+
+    Args:
+        file_name (str): Filename of sparse matrix dictionary to loda/save.
+        n_qubits (int): Number of qubits.
+        h_ij (Dict[Tuple, QubitOperator]): Operator matrix
+        load (bool, optional): If True, precomputed sparse matrix dictionary is used. Defaults to True.
+
+    Returns:
+        Dict[Tuple, scipy.sparse]: _description_
+    """
     if load and file_name in os.listdir("h4_experiment_results"):
         with open(f"h4_experiment_results/{file_name}", "rb") as f:
             pauli_mat_dict = pickle.load(f)
@@ -185,7 +239,12 @@ def _get_pauli_mat_dict(file_name, n_qubits, h_ij, load=True):
     return pauli_mat_dict
 
 
-def simulate_qse_convergence():
+def simulate_qse_convergence() -> Tuple[float, float, List[float]]:
+    """Benchmark convergence of adaptive update of measurement strategy
+
+    Returns:
+        Tuple[float, float, List[float]]: rigorous and CISD energies, and energies for iterations.
+    """
     params_qse = {
         "molecule": "H4",
         "n_qubits": 8,

reproduce_experiment.ipynb

@@ -30,6 +30,22 @@
     "# Table 1, 3 and Fig. 2"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "e651db91",
+   "metadata": {},
+   "source": [
+    "Description of Parameters:\n",
+    "* n_trial [int]:  Number of trials to gather statistics for mean absolute error and standard deviation.\n",
+    "* n_lev[\"auto\" or int]: Regularization parameter for the general eigenvalue problem. The maximum value is the dimension of the subspace. If \"auto\" is selected, n_lev is chosen to minimize the absolute error of estimation.\n",
+    "* subspace [\"1n\" or \"2n1p\"]: Choosing \"1n\" yields a subspace $S = Span\\{a_i |\\psi\\rangle \\}$, while \"2n1p\" yields a subspace $S = Span\\{a_i a_j^\\dagger a_k |\\psi\\rangle \\}$.\n",
+    "* spin_supspace [\"up\", \"down\", or \"all\"]: Selecting \"up(down)\" yields a supspace $S = Span\\{O_{i,\\sigma=\\uparrow(\\downarrow) } |\\psi\\rangle \\} $, while \"all\" yields a subspace $S = Span\\{O_{i,\\sigma\\in\\{\\uparrow, \\downarrow \\}} |\\psi\\rangle \\} $.\n",
+    "* cpu_assigned [int]: Number of CPUs assigned for the calculation.\n",
+    "* verbose [0, 1, 2]: Verbosity level for the calculation progress. To minimize output, set it to 0.\n",
+    "* load [bool]: If True, precomputed values are utilized for some calculations.\n",
+    "* write_result_matrix [bool]: If True, $\\tilde{H}_{ij}$ and $\\tilde{S}_{ij}$ are saved after the calculation.\""
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 17,
@@ -40,8 +56,8 @@
    "outputs": [],
    "source": [
     "params_fixed = {\n",
-    "    \"n_trial\": 10,\n",
-    "    \"n_lev\": \"auto\",\n",
+    "    \"n_trial\": 10, \n",
+    "    \"n_lev\": \"auto\", # parameter \n",
     "    \"subspace\": \"1n\",\n",
     "    \"spin_supspace\": \"up\",\n",
     "    \"cpu_assigned\": 10,\n",

sampler.py

@@ -1,4 +1,4 @@
-from typing import Iterable, List, Optional
+from typing import Dict, Iterable, List, Optional
 
 import numpy as np
 import pandas as pd
@@ -35,19 +35,22 @@ def generate_random_measurement_axis(
         self,
         lbcs_beta: Optional[Iterable[Iterable]] = None,
         ogm_meas_set: Optional[Iterable[Iterable]] = None,
-    ):
-        """
-        Creates a list of random measurement axis.
+    ) -> np.ndarray:
+        """   Creates a list of random measurement axis.
         For 3-qubit system, this looks like, e.g.,
         [1, 3, 2],
         which tells you to measure in
         [X, Z, Y]
         basis.
 
-        return:
-            List[n_qubit]
+        Args:
+            lbcs_beta (Optional[Iterable[Iterable]], optional): Weighted bias of measurement bases for LBCS. Defaults to None.
+            ogm_meas_set (Optional[Iterable[Iterable]], optional): Probability distribution of basis for OGM. Defaults to None.
+
+        Returns:
+            np.ndarray: Optimized measurement bases.
         """
-        assert not ((lbcs_beta is not None) and (ogm_meas_set is not None)), "you must choose either of lbcs or ogm"
+        assert not ((lbcs_beta is not None) and (ogm_meas_set is not None)), "you must choose either of LBCS or OGM"
         if lbcs_beta is not None:
             meas_axes = np.vstack(
                 [np.random.multinomial(1, beta_i, size=self.Ntot).argmax(axis=1) + 1 for beta_i in lbcs_beta]
@@ -60,12 +63,15 @@ def generate_random_measurement_axis(
             meas_axes = np.random.randint(1, 4, size=(self.Ntot, self.n_qubit))
         return meas_axes
 
-    def _sample_digits(self, meas_axis, nshot_per_axis=1):
-        """
-        returns the measurement result at meas_axis.
-        attributes:
-            meas_axis: List of axes
-            nshot_per_axis: number of measurement per axis
+    def _sample_digits(self, meas_axis: np.ndarray, nshot_per_axis=1) -> List[List[int]]:
+        """ Returns the measurement result at meas_axis.
+
+        Args:
+            meas_axis (np.ndarray): Optimized measurement basis axis.
+            nshot_per_axis (int, optional): number of measurement per axis. Defaults to 1.
+
+        Returns:
+            List[List[int]]: List of measurement result in same format of meas_axis.
         """
         meas_state = QuantumState(self.n_qubit)
         meas_state.load(self.state)
@@ -89,9 +95,26 @@ def _sample_digits(self, meas_axis, nshot_per_axis=1):
         return digits
 
 
-def local_dists_optimal(ham, num_qubits, objective, method, β_initial=None, bitstring_HF=None):
+def local_dists_optimal(
+    ham: QubitOperator,
+    num_qubits: int,
+    objective: str,
+    method: str,
+    β_initial: Dict = None,
+    bitstring_HF: str = None,
+) -> np.ndarray:    
     """Find optimal probabilities beta_{i,P} and return as dictionary
-    attn: qiskit ordering"""
+
+    Args:
+        ham (QubitOperator): Target Hamiltonian
+        num_qubits (int): Number of qubits
+        objective (str): Objective fuction
+        method (str): Optimization method
+        bitstring_HF (str, optional): HF representation. Defaults to None.
+
+    Returns:
+        np.ndarray: _description_
+    """
     assert objective in ["diagonal", "mixed"]
     assert method in ["scipy", "lagrange"]
 
@@ -138,12 +161,13 @@ def estimate_exp(
     meas_axes: Iterable[Iterable] = None,
     samples: np.ndarray = None,
 ) -> float:
-    """
-    Estimate expectation value of Observable for Basic Classical Shadow
+    """Estimate expectation value of Observable for Basic Classical Shadow
 
     Args:
         operator (QubitOperator): Observable such as Hamiltonian
         sampler (LocalPauliShadowSampler_core): Sampler Class
+        meas_axes (Iterable[Iterable], optional): Precomputed measurement axes. Defaults to None.
+        samples (np.ndarray, optional): Precomputed sampling results for given measurement axes. Defaults to None.
 
     Returns:
         float: Expectation value
@@ -179,13 +203,15 @@ def estimate_exp_lbcs(
     meas_axes: Iterable[Iterable] = None,
     samples: np.ndarray = None,
 ) -> float:
-    """
-    Estimate expectation value of Observable for Locally Biased Classical Shadow
+    """Estimate expectation value of Observable for Locally Biased Classical Shadow
 
     Args:
         operator (QubitOperator): Observable such as Hamiltonian
         sampler (LocalPauliShadowSampler_core): Sampler Class
         beta (Iterable[Iterable]): weighted bias
+        meas_axes (Iterable[Iterable], optional): Precomputed measurement axes. Defaults to None.
+        samples (np.ndarray, optional): Precomputed sampling results for given measurement axes. Defaults to None.
+
 
     Returns:
         float: Expectation value
@@ -231,6 +257,9 @@ def estimate_exp_ogm(
         sampler (LocalPauliShadowSampler_core): Sampler Class
         meas_dist (Iterable[Iterable]): measurement set and its distribution
                                         input is format of QubitOperator.terms
+        meas_axes (Iterable[Iterable], optional): Precomputed measurement axes. Defaults to None.
+        samples (np.ndarray, optional): Precomputed sampling results for given measurement axes. Defaults to None.
+
 
     Returns:
         float: Expectation value
@@ -287,6 +316,8 @@ def estimate_exp_derand(
     Args:
         operator (QubitOperator): Observable such as Hamiltonian
         sampler (LocalPauliShadowSampler_core): Sampler Class
+        meas_axes (Iterable[Iterable], optional): Precomputed measurement axes. Defaults to None.
+        samples (np.ndarray, optional): Precomputed sampling results for given measurement axes. Defaults to None.
 
     Returns:
         float: Expectation value