Expectation values as native backend arrays

src/qibo/backends/pytorch.py

@@ -120,7 +120,11 @@ def cast(
 
         if isinstance(x, self.np.Tensor):
             x = x.to(dtype)
-        elif isinstance(x, list) and all(isinstance(row, self.np.Tensor) for row in x):
+        elif (
+            isinstance(x, list)
+            and len(x) > 0
+            and all(isinstance(row, self.np.Tensor) for row in x)
+        ):
             x = self.np.stack(x)
         else:
             x = self.np.tensor(x, dtype=dtype, requires_grad=requires_grad)

src/qibo/hamiltonians/hamiltonians.py

@@ -564,7 +564,9 @@ def expectation_from_samples(self, freq, qubit_map=None):
             if len(term.factors) != len(set(term.factors)):
                 raise_error(NotImplementedError, "Z^k is not implemented since Z^2=I.")
         keys = list(freq.keys())
-        counts = np.array(list(freq.values())) / sum(freq.values())
+        counts = self.backend.cast(list(freq.values()), self.backend.precision) / sum(
+            freq.values()
+        )
         qubits = []
         for term in terms:
             qubits_term = []

src/qibo/models/error_mitigation.py

@@ -1,6 +1,7 @@
 """Error Mitigation Methods."""
 
 import math
+from inspect import signature
 from itertools import product
 
 import numpy as np
@@ -207,9 +208,12 @@ def ZNE(
         )
         expected_values.append(val)
 
-    gamma = get_gammas(noise_levels, analytical=solve_for_gammas)
+    gamma = backend.cast(
+        get_gammas(noise_levels, analytical=solve_for_gammas), backend.precision
+    )
+    expected_values = backend.cast(expected_values, backend.precision)
 
-    return np.sum(gamma * expected_values)
+    return backend.np.sum(gamma * expected_values)
 
 
 def sample_training_circuit_cdr(
@@ -302,6 +306,27 @@ def sample_training_circuit_cdr(
     return sampled_circuit
 
 
+def _curve_fit(
+    backend, model, params, xdata, ydata, lr=1, max_iter=int(1e2), tolerance_grad=1e-5
+):
+    if backend.name == "pytorch":
+        loss = lambda pred, target: backend.np.mean((pred - target) ** 2)
+        optimizer = backend.np.optim.LBFGS(
+            [params], lr=lr, max_iter=max_iter, tolerance_grad=tolerance_grad
+        )
+
+        def closure():
+            optimizer.zero_grad()
+            output = model(xdata, *params)
+            loss_val = loss(output, ydata)
+            loss_val.backward(retain_graph=True)
+            return loss_val
+
+        optimizer.step(closure)
+        return params
+    return curve_fit(model, xdata, ydata, p0=params)[0]
+
+
 def CDR(
     circuit,
     observable,
@@ -382,7 +407,17 @@ def CDR(
         )
         train_val["noisy"].append(val)
 
-    optimal_params = curve_fit(model, train_val["noisy"], train_val["noise-free"])[0]
+    nparams = (
+        len(signature(model).parameters) - 1
+    )  # first arg is the input and the *params afterwards
+    params = backend.cast(local_state.random(nparams), backend.precision)
+    optimal_params = _curve_fit(
+        backend,
+        model,
+        params,
+        backend.cast(train_val["noisy"], backend.precision),
+        backend.cast(train_val["noise-free"], backend.precision),
+    )
 
     val = get_expectation_val_with_readout_mitigation(
         circuit,
@@ -408,7 +443,7 @@ def vnCDR(
     noise_levels,
     noise_model,
     nshots: int = 10000,
-    model=lambda x, *params: (x * np.array(params).reshape(-1, 1)).sum(0),
+    model=None,
     n_training_samples: int = 100,
     insertion_gate: str = "CNOT",
     full_output: bool = False,
@@ -463,6 +498,9 @@ def vnCDR(
     """
     backend, local_state = _check_backend_and_local_state(seed, backend)
 
+    if model is None:
+        model = lambda x, *params: backend.np.sum(x * backend.np.vstack(params), axis=0)
+
     if readout is None:
         readout = {}
 
@@ -475,7 +513,7 @@ def vnCDR(
     for circ in training_circuits:
         result = backend.execute_circuit(circ, nshots=nshots)
         val = result.expectation_from_samples(observable)
-        train_val["noise-free"].append(val)
+        train_val["noise-free"].append(float(val))
         for level in noise_levels:
             noisy_c = get_noisy_circuit(circ, level, insertion_gate=insertion_gate)
             val = get_expectation_val_with_readout_mitigation(
@@ -488,12 +526,21 @@ def vnCDR(
                 seed=local_state,
                 backend=backend,
             )
-            train_val["noisy"].append(val)
-
-    noisy_array = np.array(train_val["noisy"]).reshape(-1, len(noise_levels))
+            train_val["noisy"].append(float(val))
 
-    params = local_state.random(len(noise_levels))
-    optimal_params = curve_fit(model, noisy_array.T, train_val["noise-free"], p0=params)
+    noisy_array = backend.cast(train_val["noisy"], backend.precision).reshape(
+        -1, len(noise_levels)
+    )
+    params = backend.cast(local_state.random(len(noise_levels)), backend.precision)
+    optimal_params = _curve_fit(
+        backend,
+        model,
+        params,
+        noisy_array.T,
+        backend.cast(train_val["noise-free"], backend.precision),
+        lr=1,
+        tolerance_grad=1e-7,
+    )
 
     val = []
     for level in noise_levels:
@@ -510,7 +557,10 @@ def vnCDR(
         )
         val.append(expval)
 
-    mit_val = model(np.array(val).reshape(-1, 1), *optimal_params[0])[0]
+    mit_val = model(
+        backend.cast(val, backend.precision).reshape(-1, 1),
+        *optimal_params,
+    )[0]
 
     if full_output:
         return mit_val, val, optimal_params, train_val
@@ -789,8 +839,7 @@ def get_expectation_val_with_readout_mitigation(
     exp_val = circuit_result.expectation_from_samples(observable)
 
     if "ncircuits" in readout:
-        exp_val /= circuit_result_cal.expectation_from_samples(observable)
-
+        return exp_val / circuit_result_cal.expectation_from_samples(observable)
     return exp_val
 
 
@@ -1038,8 +1087,9 @@ def ICS(
         data["noisy"].append(noisy_expectation)
         lambda_list.append(1 - noisy_expectation / expectation)
 
-    dep_param = np.mean(lambda_list)
-    dep_param_std = np.std(lambda_list)
+    lambda_list = backend.cast(lambda_list, backend.precision)
+    dep_param = backend.np.mean(lambda_list)
+    dep_param_std = backend.np.std(lambda_list)
 
     noisy_expectation = get_expectation_val_with_readout_mitigation(
         circuit,

tests/test_hamiltonians_symbolic.py

@@ -395,4 +395,4 @@ def test_symbolic_hamiltonian_with_constant(backend):
     h = hamiltonians.SymbolicHamiltonian(1e6 - Z(0), backend=backend)
 
     result = c.execute(nshots=10000)
-    assert result.expectation_from_samples(h) == approx(1e6, rel=1e-5, abs=0.0)
+    assert float(result.expectation_from_samples(h)) == approx(1e6, rel=1e-5, abs=0.0)

tests/test_models_error_mitigation.py

@@ -317,7 +317,7 @@ def test_readout_mitigation(backend, nqubits, nmeas, method, ibu_iters):
         c, obs, noise, nshots, readout, backend=backend
     )
 
-    assert np.abs(true_val - mit_val) <= np.abs(true_val - noisy_val)
+    assert backend.np.abs(true_val - mit_val) <= backend.np.abs(true_val - noisy_val)
 
 
 @pytest.mark.parametrize("nqubits", [3])