Periodic merge of develop into master

pycqed/analysis/analysis_toolbox.py

@@ -1950,21 +1950,135 @@ def plot_errorbars(x, y, ax=None, linewidth=2 ,markersize=2, marker='none'):
 
 
 def calculate_transmon_transitions(EC, EJ, asym=0, reduced_flux=0,
-                                   no_transitions=2, dim=None, ng=0):
+                                   no_transitions=2, dim=None, ng=0,
+                                   return_injs=False):
     '''
     Calculates transmon energy levels from the full transmon qubit Hamiltonian.
     '''
     if dim is None:
-        dim = no_transitions*20
+        dim = no_transitions*10
 
     EJphi = EJ*np.sqrt(asym**2 + (1-asym**2)*np.cos(np.pi*reduced_flux)**2)
     Ham = 4*EC*np.diag(np.arange(-dim-ng, dim-ng+1)**2) - EJphi/2 * \
         (np.eye(2*dim+1, k=+1) + np.eye(2*dim+1, k=-1))
-    HamEigs = np.linalg.eigvalsh(Ham)
-    HamEigs.sort()
-    transitions = HamEigs[1:]-HamEigs[:-1]
-    return transitions[:no_transitions]
 
+    if return_injs:
+        HamEigs, HamEigVs = np.linalg.eigh(Ham)
+        # HamEigs.sort()
+        transitions = HamEigs[1:]-HamEigs[:-1]
+        charge_number_operator = np.diag(np.arange(-dim-ng, dim-ng+1))
+        injs = np.zeros([dim, dim])
+        for i in range(dim):
+            for j in range(dim):
+                vect_i = np.matrix(HamEigVs[:,i])
+                vect_j = np.matrix(HamEigVs[:,j])
+                injs[i, j] = vect_i*(charge_number_operator*vect_j.getH())
+        return transitions[:no_transitions], injs
+
+    else:
+        HamEigs = np.linalg.eigvalsh(Ham)
+        HamEigs.sort()
+        transitions = HamEigs[1:]-HamEigs[:-1]
+        return transitions[:no_transitions]
+
+
+def calculate_transmon_and_resonator_transitions(EC, EJ, f_r, g_01,
+                                   dim=None, ng=0, f_01=None, f_12=None,
+                                            g_12_approximation=1):
+    '''
+    Calculates transmon energy levels and resonator from the full transmon qubit Hamiltonian.
+    '''
+
+    #calculate the bare transmon transitions, hardcoded to three levels only
+    [f_01, f_12], injs = calculate_transmon_transitions(EC, EJ, asym=0, reduced_flux=0,
+                                            no_transitions=2, dim=dim, ng=ng,
+                                            return_injs=True)
+
+    #problem can be cut up in th 0, 1 and 2-excitation manifold with E_ij, i excitations in the qubit and j of the resonator
+    try:
+        E00 = 0
+        g_01 = np.abs(g_01)
+        H1 = np.array([[f_01,   g_01],
+                       [g_01,   f_r]])
+        E10, E01 = np.linalg.eigvalsh(H1)
+        g_12_transmon = abs(g_01*injs[1,2]/injs[0,1])
+        g_12_resonator = abs(g_01*np.sqrt(2))
+        # print('g_12_correction',g_12_transmon/g_12_resonator)
+        H2 = np.array([[f_01+f_12,  g_12_transmon,       0],
+                       [g_12_transmon,       f_01+f_r,   g_12_resonator],
+                       [0,          g_12_resonator,       2*f_r]])
+        E20, E11, E02 = np.linalg.eigvalsh(H2)
+        #print(EC/1e6, EJ/1e9, f_r/1e9, g_01/1e6, f_01/1e9, f_12/1e9)
+    except np.linalg.LinAlgError:
+        print(EC/1e6, EJ/1e9, f_r/1e9, g_01/1e6, f_01/1e9, f_12/1e9)
+
+
+
+    f_01_d = E10 - E00
+    f_12_d = E20 - E10
+    f_r_d = E01 - E00
+    f_01_res_shifted = E11 - E01
+    f_r_qubit_shifted = E11 - E10
+
+    return f_01_d, f_12_d, f_r_d, f_01_res_shifted, f_r_qubit_shifted
+
+
+
+
+def fit_EC_EJ_g_f_res_ng(flux_01, f_01, flux_12, f_12, flux_r, f_r, ng=0, asym=0):
+    '''
+    Fits EC, EJ, g and f_res from a list of f01, f12 and f resonators 
+    as a function of thier respective flux settings by numerical optimization.
+    for initial guess it takes the maximum of the inputs
+    '''
+    from scipy import optimize
+    # initial guesses
+    g_01_ss_guess = 300e6
+    EC_guess = np.max(f_01)-np.max(f_12)
+    print('EC_guess',EC_guess/1e6)
+    EJmax_guess = (np.max(f_01)+EC_guess)**2/(8*EC_guess)
+
+    f_r_guess = np.min(f_r)
+    asym_guess = 0.1
+
+    def g01(g_01_ss, EC, EJ, EJmax):
+        return (EJ/EC)**(1/4)/(EJmax/EC)**(1/4)*g_01_ss
+
+    def penaltyfn (params):
+        EC, EJmax, f_r_bare, g_01_ss, asym = params
+        # calculate f01s
+        f01s = []
+        for fl in flux_01:
+            EJ = EJmax*np.sqrt(asym**2 + (1-asym**2)*np.cos(np.pi*fl)**2)
+            g_01 = g01(g_01_ss, EC, EJ, EJmax)
+            f01s.append(calculate_transmon_and_resonator_transitions(EC, EJ, f_r_bare, g_01, ng=ng)[0])
+        f01s = np.array(f01s)
+        # calculate f12s
+        f12s = []
+        for fl in flux_12:
+            EJ = EJmax*np.sqrt(asym**2 + (1-asym**2)*np.cos(np.pi*fl)**2)
+            g_01 = g01(g_01_ss, EC, EJ, EJmax)
+            f12s.append(calculate_transmon_and_resonator_transitions(EC, EJ, f_r_bare, g_01, ng=ng)[1])
+        f12s = np.array(f12s)
+        # calculate f_rs
+        f_rs = []
+        for fl in flux_r:
+            EJ = EJmax*np.sqrt(asym**2 + (1-asym**2)*np.cos(np.pi*fl)**2)
+            g_01 = g01(g_01_ss, EC, EJ, EJmax)
+            f_rs.append(calculate_transmon_and_resonator_transitions(EC, EJ, f_r_bare, g_01, ng=ng)[2])
+        f_rs = np.array(f_rs)
+
+        penalty_01 = f_01 - f01s
+        penalty_12 = f_12 - f12s
+        penalty_alpha = (f_01-f_12)-(f01s-f12s)
+        penalty_r = f_r - f_rs
+        return np.concatenate((penalty_01, penalty_12, penalty_alpha*20,  penalty_r*15))
+
+    (EC0, EJmax, fres, g_01_ss, asym), success = optimize.leastsq(penaltyfn, (EC_guess, EJmax_guess, f_r_guess, g_01_ss_guess, asym_guess))
+    print(success)
+    print('with', EC0, EJmax, fres, g_01_ss, asym)
+
+    return EC0, EJmax, fres, g_01_ss, asym
 
 def fit_EC_EJ(f01, f12):
     '''

pycqed/analysis/tools/plotting.py

@@ -9,6 +9,7 @@
 import matplotlib.colors as col
 import hsluv
 
+
 def set_xlabel(axis, label, unit=None, **kw):
     """
     Takes in an axis object and add a unit aware label to it.
@@ -61,12 +62,13 @@ def set_ylabel(axis, label, unit=None, **kw):
     return axis
 
 
+SI_PREFIXES = dict(zip(range(-24, 25, 3), 'yzafpnμm kMGTPEZY'))
+SI_PREFIXES[0] = ""
 
-SI_PREFIXES_2 = dict(zip(range(-24, 25, 3), 'yzafpnμm kMGTPEZY'))
-# SI_PREFIXES_2[0] = ""
-#SI_PREFIXES = 'yzafpnμm kMGTPEZY'
+# N.B. not all of these are SI units, however, all of these support SI prefixes
+SI_UNITS = 'm,s,g,W,J,V,A,F,T,Hz,Ohm,S,N,C,px,b,B,K,Bar,Vpeak,Vpp,Vp,Vrms'.split(
+    ',')
 
-SI_UNITS = 'm,s,g,W,J,V,A,F,T,Hz,Ohm,S,N,C,px,b,B,K,Bar,Vpeak,Vpp,Vp,Vrms'.split(',')
 
 def SI_prefix_and_scale_factor(val, unit=None):
     """
@@ -84,7 +86,12 @@ def SI_prefix_and_scale_factor(val, unit=None):
         try:
             with np.errstate(all="ignore"):
                 prefix_power = np.log10(abs(val))//3 * 3
-            return 10 ** -prefix_power, SI_PREFIXES_2[prefix_power] + unit
+                prefix = SI_PREFIXES[prefix_power]
+                # Greek symbols not supported in tex
+                if plt.rcParams['text.usetex'] and prefix == 'μ':
+                    prefix = r'$\mu$'
+
+            return 10 ** -prefix_power,  prefix + unit
         except (KeyError, TypeError):
             pass
 
@@ -276,7 +283,6 @@ def flex_colormesh_plot_vs_xy(xvals, yvals, zvals, ax=None,
     xvals = np.array(xvals)
     yvals = np.array(yvals)
 
-
     # First, we need to sort the data as otherwise we get odd plotting
     # artefacts. An example is e.g., plotting a fourier transform
     sorted_x_arguments = xvals.argsort()
@@ -346,7 +352,6 @@ def autolabel_barplot(ax, rects, rotation=90):
                 ha='center', va='bottom', rotation=rotation)
 
 
-
 def set_axeslabel_color(ax, color):
     '''
     Ad hoc function to set the labels, ticks, ticklabels and title to a color.
@@ -361,8 +366,7 @@ def set_axeslabel_color(ax, color):
     plt.setp(ax.title, color=color)
 
 
-
-##### generate custom colormaps
+# generate custom colormaps
 def make_segmented_cmap():
     white = '#ffffff'
     black = '#000000'
@@ -372,21 +376,24 @@ def make_segmented_cmap():
         'anglemap', [black, red, white, blue, black], N=256, gamma=1)
     return anglemap
 
-def make_anglemap( N = 256, use_hpl = True ):
-    h = np.ones(N) # hue
-    h[:N//2] = 11.6 # red
-    h[N//2:] = 258.6 # blue
-    s = 100 # saturation
-    l = np.linspace(0, 100, N//2) # luminosity
-    l = np.hstack( (l,l[::-1] ) )
 
-    colorlist = np.zeros((N,3))
+def make_anglemap(N=256, use_hpl=True):
+    h = np.ones(N)  # hue
+    h[:N//2] = 11.6  # red
+    h[N//2:] = 258.6  # blue
+    s = 100  # saturation
+    l = np.linspace(0, 100, N//2)  # luminosity
+    l = np.hstack((l, l[::-1]))
+
+    colorlist = np.zeros((N, 3))
     for ii in range(N):
         if use_hpl:
-            colorlist[ii,:] = hsluv.hpluv_to_rgb( (h[ii], s, l[ii]) )
+            colorlist[ii, :] = hsluv.hpluv_to_rgb((h[ii], s, l[ii]))
         else:
-            colorlist[ii,:] = hsluv.hsluv_to_rgb( (h[ii], s, l[ii]) )
-    colorlist[colorlist > 1] = 1 # correct numeric errors
+            colorlist[ii, :] = hsluv.hsluv_to_rgb((h[ii], s, l[ii]))
+    colorlist[colorlist > 1] = 1  # correct numeric errors
     colorlist[colorlist < 0] = 0
-    return col.ListedColormap( colorlist )
-hsluv_anglemap = make_anglemap( use_hpl = False )
+    return col.ListedColormap(colorlist)
+
+
+hsluv_anglemap = make_anglemap(use_hpl=False)

pycqed/analysis_v2/cryo_scope_analysis.py

@@ -196,31 +196,37 @@ def extract_data(self):
         self.raw_data_dict['amps'] = []
         self.raw_data_dict['data'] = []
 
-        a = ma_old.MeasurementAnalysis(
-            timestamp=self.timestamp, auto=False, close_file=False)
-        a.get_naming_and_values()
-
-        ch_amp = a.data_file[self.ch_amp_key].attrs['value']
-        if self.ch_range_key is None:
-            ch_range = 2  # corresponds to a scale factor of 1
-        else:
-            ch_range = a.data_file[self.ch_range_key].attrs['value']
-        amp = ch_amp*ch_range/2
-        # amp = ch_amp
-
-        # read conversion polynomial from the datafile if not provided as input
-        if isinstance(self.polycoeffs_freq_conv, str):
-            self.polycoeffs_freq_conv = np.array(
-                a.data_file[self.polycoeffs_freq_conv])
-            print(np.array(self.polycoeffs_freq_conv))
-
-        self.raw_data_dict['data'] = a.measured_values[self.ch_idx_cos] + 1j * \
-            a.measured_values[self.ch_idx_sin]
+        for i, timestamp in enumerate(self.timestamps):
+            a = ma_old.MeasurementAnalysis(
+                timestamp=timestamp, auto=False, close_file=False)
+            a.get_naming_and_values()
+            if i == 0:
+                ch_amp = a.data_file[self.ch_amp_key].attrs['value']
+                if self.ch_range_key is None:
+                    ch_range = 2  # corresponds to a scale factor of 1
+                else:
+                    ch_range = a.data_file[self.ch_range_key].attrs['value']
+                amp = ch_amp*ch_range/2
+                # read conversion polynomial from the datafile if not provided as input
+                if isinstance(self.polycoeffs_freq_conv, str):
+                    self.polycoeffs_freq_conv = np.array(
+                        a.data_file[self.polycoeffs_freq_conv])
+                    print(np.array(self.polycoeffs_freq_conv))
+
+                self.raw_data_dict['data'] =\
+                    a.measured_values[self.ch_idx_cos] + \
+                    1j * a.measured_values[self.ch_idx_sin]
+
+                # hacky but required for data saving
+                self.raw_data_dict['folder'] = a.folder
+                self.raw_data_dict['amps'].append(amp)
 
-        # hacky but required for data saving
-        self.raw_data_dict['folder'] = a.folder
-        self.raw_data_dict['amps'].append(amp)
-        a.finish()
+            else:
+                # If multiple datasets are used, shapes must match
+                self.raw_data_dict['data'] +=\
+                    a.measured_values[self.ch_idx_cos] + \
+                    1j * a.measured_values[self.ch_idx_sin]
+            a.finish()
 
         self.raw_data_dict['times'] = a.sweep_points
         self.raw_data_dict['timestamps'] = self.timestamps
@@ -243,39 +249,44 @@ def process_data(self):
 
     def prepare_plots(self):
         # pass
+        if len(self.timestamps)>1:
+            t_stamp_id = '{}-{}'.format(self.timestamps[0], self.timestamps[1])
+        else:
+            t_stamp_id = self.timestamps[0]
+
         self.plot_dicts['raw_data'] = {
             'plotfn': self.ca.plot_raw_data,
-            'title': self.timestamp+'\nRaw cryoscope data'}
+            'title': t_stamp_id+'\nRaw cryoscope data'}
 
         self.plot_dicts['demod_data'] = {
             'plotfn': self.ca.plot_demodulated_data,
-            'title': self.timestamp+'\nDemodulated data'}
+            'title': t_stamp_id+'\nDemodulated data'}
 
         self.plot_dicts['norm_data_circ'] = {
             'plotfn': self.ca.plot_normalized_data_circle,
-            'title': self.timestamp+'\nNormalized cryoscope data'}
+            'title': t_stamp_id+'\nNormalized cryoscope data'}
 
         self.plot_dicts['demod_phase'] = {
             'plotfn': self.ca.plot_phase,
-            'title': self.timestamp+'\nDemodulated phase'}
+            'title': t_stamp_id+'\nDemodulated phase'}
 
         self.plot_dicts['frequency_detuning'] = {
             'plotfn': self.ca.plot_frequency,
-            'title': self.timestamp+'\nDetuning frequency'}
+            'title': t_stamp_id+'\nDetuning frequency'}
 
         self.plot_dicts['cryoscope_amplitude'] = {
             'plotfn': self.ca.plot_amplitude,
-            'title': self.timestamp+'\nCryoscope amplitude'}
+            'title': t_stamp_id+'\nCryoscope amplitude'}
 
         self.plot_dicts['short_time_fft'] = {
             'plotfn': self.ca.plot_short_time_fft,
-            'title': self.timestamp+'\nShort time Fourier Transform'}
+            'title': t_stamp_id+'\nShort time Fourier Transform'}
 
         self.plot_dicts['zoomed_cryoscope_amplitude'] = {
             'plotfn': make_zoomed_cryoscope_fig,
             't': self.ca.time,
             'amp': self.ca.get_amplitudes(),
-            'title': self.timestamp+'\n Zoomed cryoscope amplitude'}
+            'title': t_stamp_id+'\n Zoomed cryoscope amplitude'}
 
 
 class SlidingPulses_Analysis(ba.BaseDataAnalysis):