[WIP] GCI with shot noise

compile_with_cma.sh

@@ -1,10 +1,10 @@
 #!/bin/bash
-#SBATCH --job-name=boost
-#SBATCH --output=boost.log
+#SBATCH --job-name=gci_shots
+#SBATCH --output=gci_shots.log
 
 OPTIMIZATION_METHOD="cma"
-OPTIMIZATION_CONFIG="{ \"maxiter\": 20}"
+OPTIMIZATION_CONFIG="{ \"maxiter\": 5}"
 
-python3 compiling.py --backend numpy --path "./results/XXZ_5seeds/sgd_10q_6l_27/" \
-                     --epoch 500 --steps 2 --optimization_method $OPTIMIZATION_METHOD \
-                     --optimization_config "$OPTIMIZATION_CONFIG"
+python3 compiling.py --backend numpy --path "./results/moreonXXZ_shots/sgd_10q_8l_27/" \
+                     --epoch 200 --steps 2 --optimization_method $OPTIMIZATION_METHOD \
+                     --optimization_config "$OPTIMIZATION_CONFIG" --nshots 1000 \

compiling.py

@@ -97,8 +97,13 @@ def main(args):
         args.db_rotation,
     )
 
+
     # TODO: remove hardcoded magnetic field
     eo_d_type = getattr(double_bracket_evolution_oracles, args.eo_d)
+    params = [4 - np.sin(x / 3) for x in range(20)]
+
+    loss_value = gci.loss(0.01, eo_d_type.load(params), args.db_rotation, args.nshots)
+    print("LOSSONA: ", loss_value)
 
     print(
         f"The gci mode is {gci.double_bracket_rotation_type} rotation with {eo_d_type.__name__} as the oracle.\n"
@@ -151,6 +156,7 @@ def main(args):
                 eo_d_type=eo_d_type,
                 mode=args.db_rotation,
                 method=args.optimization_method,
+                nshots=args.nshots,
                 **opt_options,
             )
             best_s = optimized_params[0]
@@ -257,6 +263,9 @@ def print_report(report: dict):
     parser.add_argument(
         "--optimization_method", default="sgd", type=str, help="Optimization method"
     )
+    parser.add_argument(
+        "--nshots", type=int, help="Number of shots in case of shot-noisy GCI"
+    )
     parser.add_argument(
         "--optimization_config",
         type=str,

run.sh

@@ -3,6 +3,31 @@
 #SBATCH --output=10qVQE_%A_%a.log
 #SBATCH --array=0-249  # Creates 250 jobs in total
 
+# Define arrays of seeds and layers
+SEEDS=(1 2 3 4 42)
+LAYERS=(3 4 5 6 7 8 9 10)
+
+# Calculate the number of layers and seeds
+NUM_SEEDS=${#SEEDS[@]}
+NUM_LAYERS=${#LAYERS[@]}
+
+# Calculate the total number of combinations
+TOTAL_COMBINATIONS=$((NUM_SEEDS * NUM_LAYERS))
+
+# Check if the job array index is within bounds
+if [ ${SLURM_ARRAY_TASK_ID} -ge ${TOTAL_COMBINATIONS} ]; then
+  echo "Error: SLURM_ARRAY_TASK_ID=${SLURM_ARRAY_TASK_ID} exceeds number of combinations."
+  exit 1
+fi
+
+# Calculate the seed and layer for the current job
+SEED_INDEX=$((SLURM_ARRAY_TASK_ID / NUM_LAYERS))
+LAYER_INDEX=$((SLURM_ARRAY_TASK_ID % NUM_LAYERS))
+
+SEED=${SEEDS[$SEED_INDEX]}
+NLAYERS=${LAYERS[$LAYER_INDEX]}
+
+# Parameters
 NQUBITS=10
 
 DBI_STEPS=0

src/boostvqe/models/dbi/group_commutator_iteration_transpiler.py

@@ -1,5 +1,6 @@
 from dataclasses import dataclass
 from enum import Enum, auto
+import logging
 
 import hyperopt
 import numpy as np
@@ -205,7 +206,7 @@ def group_commutator(
         else:
             raise_error(ValueError, "Your EvolutionOracleType is not recognized")
 
-    def loss(self, step_duration: float, eo_d, mode_dbr):
+    def loss(self, step_duration: float, eo_d, mode_dbr, nshots:int=None):
         """
         Compute loss function distance between `look_ahead` steps.
 
@@ -218,17 +219,27 @@ def loss(self, step_duration: float, eo_d, mode_dbr):
         circ = self.get_composed_circuit()
         if step_duration is not None:
             circ = self.recursion_step_circuit(step_duration, eo_d, mode_dbr) + circ
-        return self.h.expectation(circ().state())
+
+        if nshots is None:
+            expval = self.h.expectation(circ().state())
+        else: 
+            circ.add(gates.M(*range(circ.nqubits)))
+            logging.info(circ.draw())
+            freq = circ(nshots=nshots).frequencies()
+            expval = self.h.expectation_from_samples(freq)
+
+        return expval
 
     def choose_step(
         self,
         d,
         step_grid=None,
         mode_dbr=None,
+        nshots=None,
     ):
         losses = []
         for s in step_grid:
-            losses.append(self.loss(s, d, mode_dbr))
+            losses.append(self.loss(s, d, mode_dbr, nshots))
         return step_grid[np.argmin(losses)], np.min(losses), losses
 
     def get_composed_circuit(self, step_duration=None, eo_d=None):

src/boostvqe/utils.py

@@ -354,9 +354,11 @@ def callback_D_optimization(params, gci, loss_history, params_history):
     loss_history.append(gci.loss(params[0]))
 
 
-def loss_function_D(gci_params, gci, eo_d_type, mode):
+def loss_function_D(gci_params, gci, eo_d_type, mode, nshots=None):
     """``params`` has shape [s0, b_list_0]."""
-    return gci.loss(gci_params[0], eo_d_type.load(gci_params[1:]), mode)
+    loss_value = gci.loss(gci_params[0], eo_d_type.load(gci_params[1:]), mode, nshots)
+    print(f"LOSS:{loss_value}")
+    return loss_value
 
 
 def optimize_D(
@@ -367,12 +369,14 @@ def optimize_D(
     method,
     s_bounds=(1e-4, 1e-1),
     b_bounds=(0.0, 9.0),
+    nshots=None,
     maxiter=100,
 ):
     """Optimize Ising GCI model using chosen optimization `method`."""
 
     # evolutionary strategy
     if method == "cma":
+        print("ciao")
         lower_bounds = s_bounds[0] + b_bounds[0] * (len(params) - 1)
         upper_bounds = s_bounds[1] + b_bounds[1] * (len(params) - 1)
         bounds = [lower_bounds, upper_bounds]
@@ -396,38 +400,38 @@ def optimize_D(
                 func=loss_function_D,
                 x0=params,
                 bounds=bounds,
-                args=(gci, eo_d_type, mode),
+                args=(gci, eo_d_type, mode, nshots),
                 maxiter=maxiter,
             )
         elif method == "differential_evolution":
             opt_results = optimize.differential_evolution(
                 func=loss_function_D,
                 x0=params,
                 bounds=bounds,
-                args=(gci, eo_d_type, mode),
+                args=(gci, eo_d_type, mode, nshots),
                 maxiter=maxiter,
             )
         elif method == "DIRECT":
             opt_results = optimize.direct(
                 func=loss_function_D,
                 bounds=bounds,
-                args=(gci, eo_d_type, mode),
+                args=(gci, eo_d_type, mode, nshots),
                 maxiter=maxiter,
             )
         elif method == "basinhopping":
             opt_results = optimize.basinhopping(
                 func=loss_function_D,
                 x0=params,
                 niter=maxiter,
-                minimizer_kwargs={"method": "Powell", "args": (gci, eo_d_type, mode)},
+                minimizer_kwargs={"method": "Powell", "args": (gci, eo_d_type, mode, nshots)},
             )
         # scipy local minimizers
         else:
             opt_results = optimize.minimize(
                 fun=loss_function_D,
                 x0=params,
                 bounds=bounds,
-                args=(gci, eo_d_type, mode),
+                args=(gci, eo_d_type, mode, nshots),
                 method=method,
                 options={"disp": 1, "maxiter": maxiter},
             )
@@ -440,10 +444,10 @@ def optimize_d_for_dbi(
     d_type,
     method,
     s_bounds=(-1e-1, 1e-1),
-    b_bounds=(0.0, 9.0),
+    b_bounds=(0.0, 9.0),            
     maxiter=100,
 ):
-    """Optimize Ising GCI model using chosen optimization `method`."""
+    """Optimize Ising D model using chosen optimization `method`."""
 
     # evolutionary strategy
     if method == "cma":