Merge branch 'master' into matrix_power

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/circuit.py

@@ -1,5 +1,6 @@
 import collections
 import copy
+import sys
 from typing import Dict, List, Optional, Tuple, Union
 
 import numpy as np
@@ -1267,18 +1268,8 @@ def _update_draw_matrix(self, matrix, idx, gate, gate_symbol=None):
 
         return matrix, idx
 
-    def draw(self, line_wrap=70, legend=False) -> str:
-        """Draw text circuit using unicode symbols.
-
-        Args:
-            line_wrap (int): maximum number of characters per line. This option
-                split the circuit text diagram in chunks of line_wrap characters.
-            legend (bool): If ``True`` prints a legend below the circuit for
-                callbacks and channels. Default is ``False``.
-
-        Return:
-            String containing text circuit diagram.
-        """
+    def diagram(self, line_wrap: int = 70, legend: bool = False) -> str:
+        """Build the string representation of the circuit diagram."""
         # build string representation of gates
         matrix = [[] for _ in range(self.nqubits)]
         idx = [0] * self.nqubits
@@ -1369,3 +1360,43 @@ def chunkstring(string, length):
             output += table
 
         return output.rstrip("\n")
+
+    def __str__(self):
+        return self.diagram()
+
+    def draw(self, line_wrap: int = 70, legend: bool = False):
+        """Draw text circuit using unicode symbols.
+
+        Args:
+            line_wrap (int, optional): maximum number of characters per line. This option
+                split the circuit text diagram in chunks of line_wrap characters.
+                Defaults to :math:`70`.
+            legend (bool, optional): If ``True`` prints a legend below the circuit for
+                callbacks and channels. Defaults to ``False``.
+
+        Returns:
+            String containing text circuit diagram.
+        """
+        qibo.config.log.warning(
+            "Starting on qibo 0.2.13, ``Circuit.draw`` will work in-place. "
+            + "The in-place method is currently implemented as ``Circuit.display``, but "
+            + "will be renamed as ``Circuit.draw`` on release 0.2.13. "
+            + "In release 0.2.12, the in-place display of circuits is accessible as "
+            + "``Circuit.display``."
+        )
+        return self.diagram(line_wrap, legend)
+
+    def display(self, line_wrap: int = 70, legend: bool = False):
+        """Draw text circuit using unicode symbols.
+
+        Args:
+            line_wrap (int, optional): maximum number of characters per line. This option
+                split the circuit text diagram in chunks of line_wrap characters.
+                Defaults to :math:`70`.
+            legend (bool, optional): If ``True`` prints a legend below the circuit for
+                callbacks and channels. Defaults to ``False``.
+
+        Returns:
+            String containing text circuit diagram.
+        """
+        sys.stdout.write(self.diagram(line_wrap, legend))

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,47 @@ def sample_training_circuit_cdr(
     return sampled_circuit
 
 
+def _curve_fit(
+    backend, model, params, xdata, ydata, lr=1.0, max_iter=int(1e2), tolerance_grad=1e-5
+):
+    """
+    Fits a model with given parameters on the data points (x,y). This is generally based on the
+    `scipy.optimize.curve_fit` function, except for the `PyTorchBackend` which makes use of the
+    `torch.optim.LBFGS` optimizer.
+
+    Args:
+    backend (:class:`qibo.backends.Backend`): simulation engine, this is only useful for `pytorch`.
+    model (function): model to fit, it should be a callable ``model(x, *params)``.
+    params (ndarray): initial parameters of the model.
+    xdata (ndarray): x data, i.e. inputs to the model.
+    ydata (ndarray): y data, i.e. targets ``y = model(x, *params)``.
+    lr (float, optional): learning rate, defaults to ``1``. Used only in the `pytorch` case.
+    max_iter (int, optional): maximum number of iterations, defaults to ``100``. Used only in the `pytorch` case.
+    tolerance_grad (float, optional): gradient tolerance, optimization stops after reaching it, defaults to ``1e-5``. Used only in the `pytorch` case.
+
+    Returns:
+        ndarray: the optimal parameters.
+    """
+    if backend.name == "pytorch":
+        # pytorch has some problems with the `scipy.optim.curve_fit` function
+        # thus we use a `torch.optim` optimizer
+        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 +427,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 +463,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 +518,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 +533,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 +546,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 +577,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 +859,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 +1107,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)