Merge branch 'enh/cz_sims' into feat/adaptive_sampling

examples/Adaptive sampling example using a mock device object.ipynb

@@ -39,8 +39,9 @@
     "\n",
     "mock_device = Mock_Device('mock_device')\n",
     "mock_device.mw_pow(-20)\n",
-    "mock_device.res_freq(7.4e9)\n",
-    "mock_device.cw_noise_level(.0005)\n"
+    "mock_device.res_freq(7.400023457e9)\n",
+    "mock_device.cw_noise_level(.0005)\n",
+    "mock_device.acq_delay(.05)\n"
    ]
   },
   {
@@ -76,6 +77,15 @@
     "This can also be done using an adaptive `Leaner1D` object, chosing 100 points optimally in the interval. "
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mock_device.acq_delay(.05)"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -113,13 +123,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "\n",
+    "d = det.Function_Detector(mock_device.S21, value_names=['Magn', 'Phase'], value_units=['V', 'deg'])\n",
     "MC.set_sweep_function(mock_device.mw_freq)\n",
     "MC.set_sweep_function_2D(mock_device.mw_pow)\n",
     "MC.set_detector_function(d)\n",
     "MC.set_adaptive_function_parameters({'adaptive_function': adaptive.Learner2D, \n",
-    "                                     'goal':lambda l: l.npoints>3000, \n",
-    "                                     'bounds':((7.39e9, 7.41e9), \n",
+    "                                     'goal':lambda l: l.npoints>20*20, \n",
+    "                                     'bounds':((7.398e9, 7.402e9), \n",
     "                                               (-20, -10))})\n",
     "dat = MC.run(mode='adaptive')\n",
     "\n"

pycqed/analysis/tools/plot_interpolation.py

@@ -0,0 +1,67 @@
+import numpy as np
+import logging
+from scipy import interpolate
+
+def areas(ip):
+    p = ip.tri.points[ip.tri.vertices]
+    q = p[:, :-1, :] - p[:, -1, None, :]
+    areas = abs(q[:, 0, 0] * q[:, 1, 1] - q[:, 0, 1] * q[:, 1, 0]) / 2
+    return areas
+
+
+def scale(points, xy_mean, xy_scale):
+    points = np.asarray(points, dtype=float)
+    return (points - xy_mean) / xy_scale
+
+
+def unscale(points, xy_mean, xy_scale):
+    points = np.asarray(points, dtype=float)
+    return points * xy_scale + xy_mean
+
+def interpolate_heatmap(x, y, z, n=None):
+    """
+    Args:
+
+    Returns:
+        x_grid : N*1 array of x-values of the interpolated grid
+        y_grid : N*1 array of x-values of the interpolated grid
+        z_grid : N*N array of z-values that form a grid.
+
+    The output of this method can directly be used for
+        plt.imshow(z_grid, extent=extent, aspect='auto')
+        where the extent is determined by the min and max of the x_grid and
+        y_grid
+    """
+
+    points = list(zip(x, y))
+    lbrt = np.min(points, axis=0), np.max(points, axis=0)
+    lbrt = lbrt[0][0], lbrt[0][1], lbrt[1][0], lbrt[1][1]
+
+    xy_mean = np.mean([lbrt[0], lbrt[2]]), np.mean([lbrt[1], lbrt[3]])
+    xy_scale = np.ptp([lbrt[0], lbrt[2]]), np.ptp([lbrt[1], lbrt[3]])
+
+    # interpolation needs to happen on a rescaled grid, this is somewhat akin to an
+    # assumption in the interpolation that the scale of the experiment is chosen sensibly.
+    ip = interpolate.LinearNDInterpolator(scale(points, xy_mean=xy_mean, xy_scale=xy_scale),
+                                          z)
+
+    if n is None:
+        # Calculate how many grid points are needed.
+        # factor from A=√3/4 * a² (equilateral triangle)
+        n = int(0.658 / np.sqrt(areas(ip).min()))
+        n = max(n, 10)
+        if n > 500:
+            logging.warning('n: {} larger than 500'.format(n))
+            n=500
+
+    x_lin = y_lin = np.linspace(-0.5, 0.5, n)
+    # Interpolation is evaulated linearly in the domain for interpolation
+    z_grid = ip(x_lin[:, None], y_lin[None, :]).squeeze()
+
+    # x and y grid points need to be rescaled from the linearly chosen points
+    points_grid = unscale(list(zip(x_lin, y_lin)), xy_mean=xy_mean, xy_scale=xy_scale)
+    x_grid = points_grid[:, 0]
+    y_grid = points_grid[:, 1]
+
+
+    return x_grid, y_grid, (z_grid).T

pycqed/analysis/tools/plotting.py

@@ -27,7 +27,7 @@ def set_xlabel(axis, label, unit=None, **kw):
         scale_factor, unit = SI_prefix_and_scale_factor(
             val=max(abs(xticks)), unit=unit)
         formatter = matplotlib.ticker.FuncFormatter(lambda x, pos:
-                                                    x*scale_factor)
+                                                    round(x*scale_factor, 3))
         axis.xaxis.set_major_formatter(formatter)
 
         axis.set_xlabel(label+' ({})'.format(unit), **kw)
@@ -52,7 +52,7 @@ def set_ylabel(axis, label, unit=None, **kw):
         scale_factor, unit = SI_prefix_and_scale_factor(
             val=max(abs(yticks)), unit=unit)
         formatter = matplotlib.ticker.FuncFormatter(lambda x, pos:
-                                                    x*scale_factor)
+                                                    round(x*scale_factor, 3))
         axis.yaxis.set_major_formatter(formatter)
 
         axis.set_ylabel(label+' ({})'.format(unit), **kw)

pycqed/instrument_drivers/virtual_instruments/mock_device.py

@@ -1,4 +1,5 @@
 import numpy as np
+import time
 from pycqed.analysis.fitting_models import hanger_func_complex_SI
 from qcodes.instrument.base import Instrument
 from qcodes.utils.validators import Numbers, Enum, Ints
@@ -54,13 +55,17 @@ def __init__(self, name, **kw):
                                           'magn_phase', 'magn'),
                            parameter_class=ManualParameter)
 
+        self.add_parameter('acq_delay', initial_value=0,
+                           unit='s',
+                           parameter_class=ManualParameter)
+
         self.add_parameter('cw_noise_level', initial_value=0,
                            parameter_class=ManualParameter)
 
         self.add_parameter('S21', unit='V', get_cmd=self.measure_transmission)
 
     def measure_transmission(self):
-
+        time.sleep(self.acq_delay())
         # TODO: add attenuation and gain
         transmission = dBm_to_Vpeak(self.mw_pow())
 

pycqed/measurement/measurement_control.py

@@ -14,6 +14,8 @@
 from pycqed.measurement.mc_parameter_wrapper import wrap_par_to_det
 from pycqed.analysis.tools.data_manipulation import get_generation_means
 
+from pycqed.analysis.tools.plot_interpolation import interpolate_heatmap
+
 from qcodes.instrument.base import Instrument
 from qcodes.instrument.parameter import ManualParameter
 from qcodes.utils import validators as vals
@@ -30,7 +32,7 @@
     print('Could not import msvcrt (used for detecting keystrokes)')
 
 try:
-    from qcodes.plots.pyqtgraph import QtPlot
+    from qcodes.plots.pyqtgraph import QtPlot, TransformState
 except Exception:
     print('pyqtgraph plotting not supported, '
           'try "from qcodes.plots.pyqtgraph import QtPlot" '
@@ -56,8 +58,7 @@ def __init__(self, name: str,
                            initial_value=datadir,
                            vals=vals.Strings(),
                            parameter_class=ManualParameter)
-        # Soft average is currently only available for "hard"
-        # measurements. It does not work with adaptive measurements.
+        # Soft average does not work with adaptive measurements.
         self.add_parameter('soft_avg',
                            label='Number of soft averages',
                            parameter_class=ManualParameter,
@@ -264,7 +265,6 @@ def measure_soft_adaptive(self, method=None):
         self.save_optimization_settings()
         self.adaptive_function = self.af_pars.pop('adaptive_function')
         if self.live_plot_enabled():
-            # self.initialize_plot_monitor()
             self.initialize_plot_monitor_adaptive()
         for sweep_function in self.sweep_functions:
             sweep_function.prepare()
@@ -449,6 +449,8 @@ def measurement_function(self, x):
         self.iteration += 1
         if self.mode != 'adaptive':
             self.print_progress(stop_idx)
+        else:
+            self.print_progress_adaptive()
         return vals
 
     def optimization_function(self, x):
@@ -703,6 +705,81 @@ def update_plotmon_2D(self, force_update=False):
             except Exception as e:
                 logging.warning(e)
 
+    def initialize_plot_monitor_2D_interp(self, ld=0):
+        """
+        Initialize a 2D plot monitor for interpolated (adaptive) plots
+        """
+        if self.live_plot_enabled() and len(self.sweep_function_names) ==2:
+            self.time_last_2Dplot_update = time.time()
+
+            # self.secondary_QtPlot.clear()
+            slabels = self.sweep_par_names
+            sunits = self.sweep_par_units
+            zlabels = self.detector_function.value_names
+            zunits = self.detector_function.value_units
+
+            self.im_plots = []
+            self.im_plot_scatters = []
+
+            for j in range(len(self.detector_function.value_names)):
+                self.secondary_QtPlot.add(x=[0, 1],
+                                          y=[0, 1],
+                                          z=np.zeros([2,2]),
+                                          xlabel=slabels[0], xunit=sunits[0],
+                                          ylabel=slabels[1], yunit=sunits[1],
+                                          zlabel=zlabels[j], zunit=zunits[j],
+                                          subplot=j+1,
+                                          cmap='viridis')
+
+                self.im_plots.append(self.secondary_QtPlot.traces[-1])
+                self.secondary_QtPlot.add(x=[0], y=[0],
+                                          pen=None,
+                                          color=1.0, width=0,
+                                          symbol='o', symbolSize=2,
+                                          subplot=j+1)
+                self.im_plot_scatters.append(self.secondary_QtPlot.traces[-1])
+
+    def update_plotmon_2D_interp(self, force_update=False):
+        '''
+        Updates the interpolated 2D heatmap
+        '''
+        if self.live_plot_enabled() and len(self.sweep_function_names) ==2:
+            try:
+                if (time.time() - self.time_last_2Dplot_update >
+                        self.plotting_interval() or force_update):
+                    # exists to force reset the x- and y-axis scale
+                    new_sc = TransformState(0, 1, True)
+
+                    x_vals = self.dset[:, 0]
+                    y_vals = self.dset[:, 1]
+                    for j in range(len(self.detector_function.value_names)):
+                        z_ind = len(self.sweep_functions) + j
+                        z_vals = self.dset[:, z_ind]
+
+                        # Interpolate points
+                        x_grid, y_grid, z_grid  = interpolate_heatmap(
+                            x_vals, y_vals, z_vals)
+                        # trace = self.secondary_QtPlot.traces[j]
+                        trace = self.im_plots[j]
+                        trace['config']['x'] = x_grid
+                        trace['config']['y'] = y_grid
+                        trace['config']['z'] = z_grid
+                        # force rescale the axes
+                        trace['plot_object']['scales']['x'] = new_sc
+                        trace['plot_object']['scales']['y'] = new_sc
+
+                        # Mark all measured points on which the interpolation
+                        # is based
+                        trace = self.im_plot_scatters[j]
+                        trace['config']['x'] = x_vals
+                        trace['config']['y'] = y_vals
+
+                    self.time_last_2Dplot_update = time.time()
+                    self.secondary_QtPlot.update_plot()
+            except Exception as e:
+                logging.warning(e)
+
+
     def initialize_plot_monitor_adaptive(self):
         '''
         Uses the Qcodes plotting windows for plotting adaptive plot updates
@@ -711,27 +788,43 @@ def initialize_plot_monitor_adaptive(self):
             return self.initialize_plot_monitor_adaptive_cma()
         else:
             self.initialize_plot_monitor()
-        self.time_last_ad_plot_update = time.time()
-        self.secondary_QtPlot.clear()
+            self.time_last_ad_plot_update = time.time()
+            self.secondary_QtPlot.clear()
+            self.initialize_plot_monitor_2D_interp()
+
+            zlabels = self.detector_function.value_names
+            zunits = self.detector_function.value_units
+
+            self.iter_traces = []
+
+            # Because of a bug in QCoDes pytqtgraph backend we don't
+            # want line plots and heatmaps in the same plotmon
+            # this if statement prevents that from happening
+            if len(self.sweep_functions) ==2:
+                iter_plotmon = self.main_QtPlot
+                iter_plotmon.win.nextRow()
+                iter_start_idx = len(self.sweep_functions)*len(zlabels)
+            else:
+                iter_plotmon = self.secondary_QtPlot
+                iter_start_idx=0
+
+            for j in range(len(zlabels)):
+                iter_plotmon.add(x=[0], y=[0],
+                                          xlabel='iteration',
+                                          ylabel=zlabels[j], yunit=zunits[j],
+                                          subplot=j+1+iter_start_idx,
+                                          symbol='o',
+                                          symbolSize=5)
+                self.iter_traces.append(iter_plotmon.traces[-1])
+
 
-        zlabels = self.detector_function.value_names
-        zunits = self.detector_function.value_units
 
-        for j in range(len(self.detector_function.value_names)):
-            self.secondary_QtPlot.add(x=[0],
-                                      y=[0],
-                                      xlabel='iteration',
-                                      ylabel=zlabels[j],
-                                      yunit=zunits[j],
-                                      subplot=j+1,
-                                      symbol='o', symbolSize=5)
 
     def update_plotmon_adaptive(self, force_update=False):
         if self.adaptive_function.__module__ == 'cma.evolution_strategy':
             return self.update_plotmon_adaptive_cma(force_update=force_update)
         else:
             self.update_plotmon(force_update=force_update)
-
         if self.live_plot_enabled():
             try:
                 if (time.time() - self.time_last_ad_plot_update >
@@ -740,12 +833,14 @@ def update_plotmon_adaptive(self, force_update=False):
                         y_ind = len(self.sweep_functions) + j
                         y = self.dset[:, y_ind]
                         x = range(len(y))
-                        self.secondary_QtPlot.traces[j]['config']['x'] = x
-                        self.secondary_QtPlot.traces[j]['config']['y'] = y
+                        self.iter_traces[j]['config']['x'] = x
+                        self.iter_traces[j]['config']['y'] = y
                         self.time_last_ad_plot_update = time.time()
                         self.secondary_QtPlot.update_plot()
             except Exception as e:
                 logging.warning(e)
+        self.update_plotmon_2D_interp(force_update=force_update)
+
 
     def initialize_plot_monitor_adaptive_cma(self):
         '''
@@ -1239,6 +1334,19 @@ def print_progress(self, stop_idx=None):
                 end_char = '\n'
             print('\r', progress_message, end=end_char)
 
+    def print_progress_adaptive(self):
+        if self.verbose():
+            acquired_points = self.dset.shape[0]
+
+            elapsed_time = time.time() - self.begintime
+            progress_message = \
+                "\rAcquired {acquired_points} points, \telapsed time: "\
+                "{t_elapsed}s".format(
+                    acquired_points = acquired_points,
+                    t_elapsed=round(elapsed_time, 1))
+            end_char = ''
+            print('\r', progress_message, end=end_char)
+
     def is_complete(self):
         """
         Returns True if enough data has been acquired.

pycqed/measurement/waveform_control_CC/waveform.py

@@ -384,7 +384,7 @@ def martinis_flux_pulse(length: float, lambda_2: float, lambda_3: float,
     scale = t[-1]/t_samples[-1]
     interp_wave = scipy.interpolate.interp1d(
         t/scale, theta_wave_clipped, bounds_error=False,
-        fill_value=0)(t_samples)
+        fill_value='extrapolate')(t_samples)
 
     # Return in the specified units
     if return_unit == 'theta':
@@ -414,6 +414,9 @@ def martinis_flux_pulse(length: float, lambda_2: float, lambda_3: float,
     # why sometimes the last sample is nan is not known,
     # but we will surely figure it out someday.
     # (Brian and Adriaan, 14.11.2017)
+    # This may be caused by the fill_value of the interp_wave (~30 lines up)
+    # that was set to 0 instead of extrapolate. This caused
+    # the np.tan(interp_wave) to divide by zero. (MAR 10-05-2018)
     voltage_wave = np.nan_to_num(voltage_wave)
 
     return voltage_wave