Python matplotlib.pyplot.semilogx() Examples

def test_mimo(self):
      # MIMO
      B = np.matrix('1,0;0,1')
      D = np.matrix('0,0')
      sysMIMO = ss(self.A,B,self.C,D)

      frqMIMO = sysMIMO.freqresp(self.omega)
      tfMIMO = tf(sysMIMO)

      #bode(sysMIMO) # - should throw not implemented exception
      #bode(tfMIMO) # - should throw not implemented exception

      #plt.figure(3)
      #plt.semilogx(self.omega,20*np.log10(np.squeeze(frq[0])))

      #plt.figure(4)
      #bode(sysMIMO,self.omega) 

def plot_series(series):
    plt.figure(1)
    # colors = [np.array([1, 0.1, 0.1]), np.array([0.1, 1, 0.1]), np.array([0.1, 0.1, 1])]
    colors = ['m', 'g', 'r', 'b', 'y']
    for i, s in enumerate(series):
        print(s['x'], s['y'], s['std'], s['label'])
        small_number = np.ones_like(s['x']) * (s['x'][1]*0.1)
        x_axis = np.where(s['x'] == 0, small_number, s['x'])
        plt.plot(x_axis, s['y'], color=colors[i], label=s['label'])
        plt.fill_between(x_axis, s['y'] - s['std'], s['y'] + s['std'], color=colors[i], alpha=0.2)
    plt.semilogx()
    plt.xlabel('MI reward bonus')
    plt.ylabel('Final intrinsic reward')
    plt.title('Final intrinsic reward in pointMDP with 10 good modes')
    plt.legend(loc='best')
    plt.show() 

def plot_losses(loss_vals, loss_names, filename, title, xlabel, ylabel, spacing=0):
    """
    Given a list of errors, plot the objectives of the training and show
    """
    plt.close('all')
    for li, lvals in enumerate(loss_vals):
        iterations = range(len(lvals))
        # lvals.insert(0, 0)
        if spacing == 0:
            plt.loglog(iterations, lvals, '-',label=loss_names[li])
            # plt.semilogx(iterations, lvals, 'x-')
        else:
            xvals = [ii*spacing for ii in iterations]
            plt.loglog( xvals, lvals, '-',label=loss_names[li])

    plt.grid()
    plt.legend(loc='upper left')
    plt.title(title)
    plt.xlabel(xlabel)
    plt.ylabel(ylabel)
    plt.savefig(filename)
    plt.close('all') 

def test_k_coeffs_binned(self):
        wavelengths = np.exp(np.arange(np.log(0.31e-6), np.log(29e-6), 1./20))
        wavelength_bins = np.array([wavelengths[0:-1], wavelengths[1:]]).T
        
        xsec_calc = TransitDepthCalculator(method="xsec")
        xsec_calc.change_wavelength_bins(wavelength_bins)
        ktab_calc = TransitDepthCalculator(method="ktables")
        ktab_calc.change_wavelength_bins(wavelength_bins)
                
        wavelengths, xsec_depths = xsec_calc.compute_depths(R_sun, M_jup, R_jup, 300, logZ=1, CO_ratio=1.5)
        wavelengths, ktab_depths = ktab_calc.compute_depths(R_sun, M_jup, R_jup, 300, logZ=1, CO_ratio=1.5)
        
        diffs = np.abs(ktab_depths - xsec_depths)

        '''plt.semilogx(wavelengths, xsec_depths)
        plt.semilogx(wavelengths, ktab_depths)
        plt.figure()
        plt.semilogx(wavelengths, 1e6 * diffs)
        plt.show()'''
        
        self.assertTrue(np.median(diffs) < 10e-6)
        self.assertTrue(np.percentile(diffs, 95) < 20e-6)
        self.assertTrue(np.max(diffs) < 30e-6) 

def sep_2D_plot(x, y, x_name, y_name, molecules, semilogx=False, savePic=False):
    """2D plot to compare all the molecules
    :param x: data for x-axis
    :param y: data for y-axis, y[i,j], where i is for a molecule, j is for a data point corresponding to a x value
    :return: if savePic is True, then .png image, otherwise a 'matplotlib.figure.Figure' object
    """
    fig = plt.figure()
    plt.xlabel(x_name)
    plt.ylabel(y_name)
    lineList = []
    labelList = []
    for i in range(len(molecules)):
        if semilogx:
            line, = plt.semilogx(x, y[i], 'o-')
        else:
            line, = plt.plot(x, y[i], 'o-')
        lineList.append(line)
        lineLabel = 'Comp {0}'.format(molecules[i])
        labelList.append(lineLabel)
        lgd = fig.legend(lineList, labelList, loc='upper center', bbox_to_anchor=(0.5, 1), ncol=len(molecules) / 2)
    if savePic:
        file_name = y_name.replace(' ', '_')
        fig.savefig('{}.png'.format(file_name), bbox_extra_artists=(lgd,), bbox_inches='tight')
    return fig 

def plot_summarize(prefix, trial, title=None):
    res = summarize(prefix, trial)
    err = res[:, 1:, 0]
    num = res[:, 1:, 3]
    marker = ('bo', 'r^', 'gs', 'mp', 'c*')
    for i in range(err.shape[1]):
        plt.semilogx(num[:, i], err[:, i], marker[i], ms=12)
    plt.legend(['Proposed', 'inTrees', 'NH', 'BATree', 'DTree2'], numpoints=1)
    plt.xlabel('# of Rules', fontsize=20)
    plt.ylabel('Test Error', fontsize=20)
    if title is None:
        title = prefix
    plt.title(title, fontsize=24)
    plt.show()
    if not os.path.exists('./result/fig/'):
        os.mkdir('./result/fig/')
    plt.savefig('./result/fig/compare_%s.pdf' % (prefix,), format="pdf", bbox_inches="tight")
    plt.close() 

def db_magnitude_distance_by_trt(db1, dist_type,
        figure_size=(7, 5), filename=None, filetype="png", dpi=300):
    """
    Plot magnitude-distance comparison by tectonic region
    """
    trts=[]
    for i in db1.records:
        trts.append(i.event.tectonic_region)
    trt_types=list(set(trts))
    selector = SMRecordSelector(db1)
    plt.figure(figsize=figure_size)
    for trt in trt_types:
        subdb = selector.select_trt_type(trt, as_db=True)
        mag, dists = get_magnitude_distances(subdb, dist_type)
        plt.semilogx(dists, mag, "o", mec='k', mew=0.5, label=trt)
    plt.xlabel(DISTANCE_LABEL[dist_type], fontsize=14)
    plt.ylabel("Magnitude", fontsize=14)
    plt.title("Magnitude vs Distance by Tectonic Region", fontsize=18)
    plt.legend(loc='lower right', numpoints=1)
    plt.grid()
    _save_image(filename, filetype, dpi)
    plt.show() 

def db_magnitude_distance_by_site(db1, dist_type, classification="NEHRP",
        figure_size=(7, 5), filename=None, filetype="png", dpi=300):
    """
    Plot magnitude-distance comparison by site NEHRP or Eurocode 8 Site class   
    """ 
    if classification == "NEHRP":
        site_bounds = NEHRP_BOUNDS
    elif classification == "EC8":
        site_bounds = EC8_BOUNDS
    else:
        raise ValueError("Unrecognised Site Classifier!")
    selector = SMRecordSelector(db1)
    plt.figure(figsize=figure_size)
    total_idx = []
    for site_class in site_bounds.keys():
        site_idx = _site_selection(db1, site_class, classification)
        if site_idx:
            site_db = selector.select_records(site_idx, as_db=True)
            mags, dists = get_magnitude_distances(site_db, dist_type)
            plt.plot(np.array(dists), np.array(mags), "o", mec='k',
                     mew=0.5, label="Site Class %s" % site_class)
            total_idx.extend(site_idx)
    unc_idx = set(range(db1.number_records())).difference(set(total_idx))
    unc_db = selector.select_records(unc_idx, as_db=True)
    mag, dists = get_magnitude_distances(site_db, dist_type)
    plt.semilogx(np.array(dists), np.array(mags), "o", mfc="None", mec='k',
                 mew=0.5, label="Unclassified", zorder=0)
    plt.xlabel(DISTANCE_LABEL[dist_type], fontsize=14)
    plt.ylabel("Magnitude", fontsize=14)
    plt.grid()
    plt.legend(ncol=2,loc="lower right", numpoints=1)
    plt.title("Magnitude vs Distance (by %s Site Class)" % classification,
              fontsize=18)
    _save_image(filename, filetype, dpi)
    plt.show() 

def plotBode2(zpk, n=200, f_range=None, f_logspace=True):
    """
    Bode plot of ZerosAndPoles object using matplotlib
    """
    (f, y) = zpk.frequencyResponse(n=n, f_range=f_range, f_logspace=f_logspace)

    y_A = numpy.abs(y)
    y_phi = Misc.to_deg(Misc.continuousAngle(y))

    plt.figure()
    plt.subplot(211)
    if f_logspace:
        plt.loglog(f, y_A)
    else:
        plt.plot(f, y_A)
    plt.grid(True, which="both")
    plt.ylabel("Amplitude")

    plt.subplot(212)
    if f_logspace:
        plt.semilogx(f, y_phi)
    else:
        plt.plot(f, y_A)
    plt.grid(True, which="both")
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Phase [deg]")

    plt.show() 

def drawValidationCurve(self):
        """
        To draw the validation curve
        :return:NA
        """
        X, y = self.X_train, self.y_train.ravel()
        indices = np.arange(y.shape[0])
        np.random.shuffle(indices)
        X, y = X[indices], y[indices]

        train_sizes = np.logspace(-5,5)
        train_scores, valid_scores = validation_curve(self.clf, X, y, "C",
                                              train_sizes, cv=5)

        train_scores_mean = np.mean(train_scores, axis=1)
        train_scores_std = np.std(train_scores, axis=1)
        valid_scores_mean = np.mean(valid_scores, axis=1)
        valid_scores_std = np.std(valid_scores, axis=1)

        plt.fill_between(train_sizes, train_scores_mean - train_scores_std,
                     train_scores_mean + train_scores_std, alpha=0.1,
                     color="r")
        plt.fill_between(train_sizes, valid_scores_mean - valid_scores_std,
                         valid_scores_mean + valid_scores_std, alpha=0.1, color="g")
        plt.semilogx(train_sizes, train_scores_mean, 'o-', color="r",
                 label="Training Precision")
        plt.semilogx(train_sizes, valid_scores_mean, '*-', color="g",
                 label="Cross-validation Precision")

        plt.legend(loc="best")

        plt.xlabel('Tradeoff C')
        plt.ylabel('Precision')
        plt.title('Validation Curve with Logistic Regression on the parameter of C')
        plt.grid(True)
        plt.show() 

def plot_t(EM, HS, title, i):
    plt.figure(title, figsize=(10, 8))
    plt.subplot(i)
    plt.semilogx(time, EM)
    plt.semilogx(time, HS, '--')


###############################################################################
# Impulse HS 

def plot_f(EM, HS, title, i):
    plt.figure(title, figsize=(10, 8))
    plt.subplot(i)
    plt.semilogx(1/time, EM.real)
    plt.semilogx(1/time, HS.real, '--')
    plt.semilogx(1/time, EM.imag)
    plt.semilogx(1/time, HS.imag, '--')


###############################################################################
# Halfspace 

def get_plot_options():
  """ Given a list of options, this reads them from the command line and returns
      a namespace with the values.
  """
  parser = argparse.ArgumentParser(description='Plotting.')
  parser.add_argument('--file', default='', help='File path of single plot file.')
  parser.add_argument('--filelist', default='', help='File name containing file paths.')
  parser.add_argument('--type', default='semilogy', help='Type of plot. Default is ' +
                      'semilogy, other options are plot, loglog, semilogx.')
  parser.add_argument('--title', help='Title of plot.')
  options = parser.parse_args()

  return options 

def plot_roc_curve(fpr, tpr, name='model', fig=None):
        if fig is None:
            fig = plt.figure()
            plt.semilogx(np.arange(0, 1, 0.01), np.arange(0, 1, 0.01), 'r', linestyle='--', label='Random guess')
        plt.semilogx(fpr, tpr, color=(random.uniform(0, 1), random.uniform(0, 1), random.uniform(0, 1)),
                     label='ROC curve with {}'.format(name))
        plt.title('Receiver Operating Characteristic')
        plt.xlabel('False Positive Rate')
        plt.ylabel('True Positive Rate')
        plt.legend(loc='best')
        return fig 

def plot_results(data):
    # plot the pc_attacks_detected vs threshold_vals
    plt.figure()
    plt.plot(data["threshold_vals"], data["pc_attacks_detected"], marker="o")
    plt.xlabel("Threshold Values")
    plt.ylabel("% of Total Attacks Detected")
    #    plt.title('% Total Attacks Detected vs. Threshold Value')

    # plot the total # false positives vs threshold values
    plt.figure()
    plt.semilogy(data["threshold_vals"], data["num_FP"], marker="o")
    plt.xlabel("Threshold Values")
    plt.ylabel("Total # False Positives")
    #    plt.title('Total # False Positives vs. Threshold Value')

    plt.figure()
    plt.semilogx(data["num_FP"], data["pc_attacks_detected"], marker="o", linewidth=5)
    plt.xlabel("Total # False Positives")
    plt.ylabel("% of Total Attacks Detected")
    plt.yticks([0, 10, 20, 30, 40, 50])
    plt.xlim(max(data["num_FP"]), 0)
    plt.subplots_adjust(bottom=0.2)

    font = {"family": "normal", "weight": "bold", "size": 35}
    matplotlib.rc("font", **font)

    plt.show() 

def plot(self, confidence_level: float=None, show: bool=True, file_name: str=None):
        """
        Plot the linear plot line.

        :confidence_level: the desired confidence level
        :show: True if the plot is to be shown
        :file_name: Save the plot as "file_name"
        """
        if confidence_level:
            self._set_confidence_level(confidence_level)

        plt.semilogx(self.cdf_x, _ftolnln(self.cdf[0]))
        axis = plt.gca()

        axis.grid(True, which='both')

        formatter = mpl.ticker.FuncFormatter(_weibull_ticks)
        axis.yaxis.set_major_formatter(formatter)
        yt_F = np.array([0.001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5,
                         0.6, 0.7, 0.8, 0.9, 0.95, 0.99])
        yt_lnF = _ftolnln(yt_F)
        plt.yticks(yt_lnF)

        plt.ylim(yt_lnF[1], yt_lnF[-1])
        plt.xlim(self.cdf_x.min(), self.cdf_x.max())

        self._plot_annotate()

        plt.ylabel('failure rate')
        plt.xlabel('cycles')

        if file_name:
            plt.savefig(file_name)
        if show:
            plt.show() 

def db_magnitude_distance(db1, dist_type, figure_size=(7, 5),
        figure_title=None,filename=None, filetype="png", dpi=300):
    """
    Creates a plot of magnitude verses distance for a strong motion database
    """
    plt.figure(figsize=figure_size)
    mags, dists = get_magnitude_distances(db1, dist_type)
    plt.semilogx(np.array(dists), np.array(mags), "o", mec='k', mew=0.5)
    plt.xlabel(DISTANCE_LABEL[dist_type], fontsize=14)
    plt.ylabel("Magnitude", fontsize=14)
    if figure_title:
        plt.title(figure_title, fontsize=18)
    _save_image(filename, filetype, dpi)
    plt.grid()
    plt.show() 

def test_isothermal(self):
        Ts = 5700
        Tp = 1500
        p = Profile()
        p.set_isothermal(Tp)
        calc = EclipseDepthCalculator()
        wavelengths, depths, info_dict = calc.compute_depths(p, R_sun, M_jup, R_jup, Ts, full_output=True)
                
        blackbody = np.pi * 2*h*c**2/wavelengths**5/(np.exp(h*c/wavelengths/k_B/Tp) - 1)

        rel_diffs = (info_dict["planet_spectrum"] - blackbody)/blackbody

        plt.loglog(1e6 * wavelengths, 1e-3 * blackbody, label="Blackbody")
        plt.loglog(1e6 * wavelengths, 1e-3 * info_dict["planet_spectrum"], label="PLATON")
        plt.xlabel("Wavelength (micron)", fontsize=12)
        plt.ylabel("Planet flux (erg/s/cm$^2$/micron)", fontsize=12)
        plt.legend()
        plt.tight_layout()
        plt.figure()
        plt.semilogx(1e6 * wavelengths, 100 * rel_diffs)
        plt.xlabel("Wavelength (micron)", fontsize=12)
        plt.ylabel("Relative difference (%)", fontsize=12)
        plt.tight_layout()
        plt.show()
        
        # Should be exact, but in practice isn't, due to our discretization
        self.assertLess(np.percentile(np.abs(rel_diffs), 50), 0.02)
        self.assertLess(np.percentile(np.abs(rel_diffs), 99), 0.05)
        self.assertLess(np.max(np.abs(rel_diffs)), 0.1)

        blackbody_star = np.pi * 2*h*c**2/wavelengths**5/(np.exp(h*c/wavelengths/k_B/Ts) - 1)        
        approximate_depths = blackbody / blackbody_star * (R_jup/R_sun)**2
        # Not expected to be very accurate because the star is not a blackbody
        self.assertLess(np.median(np.abs(approximate_depths - depths)/approximate_depths), 0.2) 

def test_ktables_binned(self):
        wavelengths = np.exp(np.arange(np.log(0.31e-6), np.log(29e-6), 1./20))
        wavelengths = np.append(wavelengths[0:20], wavelengths[50:90])
        wavelength_bins = np.array([wavelengths[0:-1], wavelengths[1:]]).T

        profile = Profile()
        profile.set_from_radiative_solution(
            5052, 0.75 * R_sun, 0.03142 * AU, 1.129 * M_jup, 1.115 * R_jup,
            0.983, -1.77, -0.44, -0.56, 0.23)
        xsec_calc = EclipseDepthCalculator(method="xsec")
        xsec_calc.change_wavelength_bins(wavelength_bins)
        ktab_calc = EclipseDepthCalculator(method="ktables")
        ktab_calc.change_wavelength_bins(wavelength_bins)
        xsec_wavelengths, xsec_depths = xsec_calc.compute_depths(
            profile, 0.75 * R_sun, 1.129 * M_jup, 1.115 * R_jup,
            5052)
        ktab_wavelengths, ktab_depths = ktab_calc.compute_depths(
            profile, 0.75 * R_sun, 1.129 * M_jup, 1.115 * R_jup,
            5052)
        rel_diffs = np.abs(ktab_depths - xsec_depths)/ ktab_depths
        self.assertTrue(np.median(rel_diffs) < 0.03)
        self.assertTrue(np.percentile(rel_diffs, 95) < 0.15)
        self.assertTrue(np.max(rel_diffs) < 0.2)
        
        '''print(np.median(rel_diffs), np.percentile(rel_diffs, 95), np.max(rel_diffs))
        plt.loglog(xsec_wavelengths, xsec_depths)
        plt.loglog(ktab_wavelengths, ktab_depths)
        plt.figure()
        plt.semilogx(ktab_wavelengths, rel_diffs)
        plt.show()''' 

def test_superimpose(self):
      # Test to make sure that multiple calls to plots superimpose their
      # data on the same axes unless told to do otherwise

      # Generate two plots in a row; should be on the same axes
      plt.figure(1); plt.clf()
      ctrl.bode_plot(ctrl.tf([1], [1,2,1]))
      ctrl.bode_plot(ctrl.tf([5], [1, 1]))

      # Check to make sure there are two axes and that each axes has two lines
      self.assertEqual(len(plt.gcf().axes), 2)
      for ax in plt.gcf().axes:
         # Make sure there are 2 lines in each subplot
         assert len(ax.get_lines()) == 2
      
      # Generate two plots as a list; should be on the same axes
      plt.figure(2); plt.clf();
      ctrl.bode_plot([ctrl.tf([1], [1,2,1]), ctrl.tf([5], [1, 1])])

      # Check to make sure there are two axes and that each axes has two lines
      self.assertEqual(len(plt.gcf().axes), 2)
      for ax in plt.gcf().axes:
         # Make sure there are 2 lines in each subplot
         assert len(ax.get_lines()) == 2

      # Generate two separate plots; only the second should appear
      plt.figure(3); plt.clf();
      ctrl.bode_plot(ctrl.tf([1], [1,2,1]))
      plt.clf()
      ctrl.bode_plot(ctrl.tf([5], [1, 1]))

      # Check to make sure there are two axes and that each axes has one line
      self.assertEqual(len(plt.gcf().axes), 2)
      for ax in plt.gcf().axes:
         # Make sure there is only 1 line in the subplot
         assert len(ax.get_lines()) == 1

      # Now add a line to the magnitude plot and make sure if is there
      for ax in plt.gcf().axes:
         if ax.get_label() == 'control-bode-magnitude':
            break
      ax.semilogx([1e-2, 1e1], 20 * np.log10([1, 1]), 'k-')
      self.assertEqual(len(ax.get_lines()), 2) 

def design():
    """Show results of designs"""
    # equal weighting on each output
    k1, gam1 = synth(0.25, 0.25)
    # increase "bandwidth" of output 2 by moving crossover weighting frequency 100 times higher
    k2, gam2 = synth(0.25, 25)
    # now weight output 1 more heavily
    # won't plot this one, just want gamma
    _, gam3 = synth(25, 0.25)

    print('design 1 gamma {:.3g} (Skogestad: 2.80)'.format(gam1))
    print('design 2 gamma {:.3g} (Skogestad: 2.92)'.format(gam2))
    print('design 3 gamma {:.3g} (Skogestad: 6.73)'.format(gam3))

    # do the designs
    g = plant()
    w = np.logspace(-2, 2, 101)
    I = ss([], [], [], np.eye(2))
    s1 = I.feedback(g*k1)
    s2 = I.feedback(g*k2)

    # frequency response
    sv1 = triv_sigma(s1, w)
    sv2 = triv_sigma(s2, w)

    plt.figure(2)

    plt.subplot(1, 2, 1)
    plt.semilogx(w, 20*np.log10(sv1[:, 0]), label=r'$\sigma_1(S_1)$')
    plt.semilogx(w, 20*np.log10(sv1[:, 1]), label=r'$\sigma_2(S_1)$')
    plt.semilogx(w, 20*np.log10(sv2[:, 0]), label=r'$\sigma_1(S_2)$')
    plt.semilogx(w, 20*np.log10(sv2[:, 1]), label=r'$\sigma_2(S_2)$')
    plt.ylim([-60, 10])
    plt.ylabel('magnitude [dB]')
    plt.xlim([1e-2, 1e2])
    plt.xlabel('freq [rad/s]')
    plt.legend()
    plt.title('Singular values of S')

    # time response

    # in design 1, both outputs have an inverse initial response; in
    # design 2, output 2 does not, and is very fast, while output 1
    # has a larger initial inverse response than in design 1
    time = np.linspace(0, 10, 301)
    t1 = (g*k1).feedback(I)
    t2 = (g*k2).feedback(I)

    y1 = step_opposite(t1, time)
    y2 = step_opposite(t2, time)

    plt.subplot(1, 2, 2)
    plt.plot(time, y1[0], label='des. 1 $y_1(t))$')
    plt.plot(time, y1[1], label='des. 1 $y_2(t))$')
    plt.plot(time, y2[0], label='des. 2 $y_1(t))$')
    plt.plot(time, y2[1], label='des. 2 $y_2(t))$')
    plt.xlabel('time [s]')
    plt.ylabel('response [1]')
    plt.legend()
    plt.title('c/l response to reference [1,-1]') 

def design():
    """Show results of designs"""
    # equal weighting on each output
    k1, gam1 = synth(0.25, 0.25)
    # increase "bandwidth" of output 2 by moving crossover weighting frequency 100 times higher
    k2, gam2 = synth(0.25, 25)
    # now weight output 1 more heavily
    # won't plot this one, just want gamma
    _, gam3 = synth(25, 0.25)

    print('design 1 gamma {:.3g} (Skogestad: 2.80)'.format(gam1))
    print('design 2 gamma {:.3g} (Skogestad: 2.92)'.format(gam2))
    print('design 3 gamma {:.3g} (Skogestad: 6.73)'.format(gam3))

    # do the designs
    g = plant()
    w = np.logspace(-2, 2, 101)
    I = ss([], [], [], np.eye(2))
    s1 = I.feedback(g*k1)
    s2 = I.feedback(g*k2)

    # frequency response
    sv1 = triv_sigma(s1, w)
    sv2 = triv_sigma(s2, w)

    plt.figure(2)

    plt.subplot(1, 2, 1)
    plt.semilogx(w, 20*np.log10(sv1[:, 0]), label=r'$\sigma_1(S_1)$')
    plt.semilogx(w, 20*np.log10(sv1[:, 1]), label=r'$\sigma_2(S_1)$')
    plt.semilogx(w, 20*np.log10(sv2[:, 0]), label=r'$\sigma_1(S_2)$')
    plt.semilogx(w, 20*np.log10(sv2[:, 1]), label=r'$\sigma_2(S_2)$')
    plt.ylim([-60, 10])
    plt.ylabel('magnitude [dB]')
    plt.xlim([1e-2, 1e2])
    plt.xlabel('freq [rad/s]')
    plt.legend()
    plt.title('Singular values of S')

    # time response

    # in design 1, both outputs have an inverse initial response; in
    # design 2, output 2 does not, and is very fast, while output 1
    # has a larger initial inverse response than in design 1
    time = np.linspace(0, 10, 301)
    t1 = (g*k1).feedback(I)
    t2 = (g*k2).feedback(I)

    y1 = step_opposite(t1, time)
    y2 = step_opposite(t2, time)

    plt.subplot(1, 2, 2)
    plt.plot(time, y1[0], label='des. 1 $y_1(t))$')
    plt.plot(time, y1[1], label='des. 1 $y_2(t))$')
    plt.plot(time, y2[0], label='des. 2 $y_1(t))$')
    plt.plot(time, y2[1], label='des. 2 $y_2(t))$')
    plt.xlabel('time [s]')
    plt.ylabel('response [1]')
    plt.legend()
    plt.title('c/l response to reference [1,-1]') 

def plot_results(nsweeps):
    """
    Plotting routine for iteration counts

    Args:
        nsweeps: number of fine sweeps used
    """

    fname = 'data/results_conv_diffusion_NS' + str(nsweeps) + '.pkl'
    file = open(fname, 'rb')
    results_diff = pickle.load(file)
    file.close()

    fname = 'data/results_conv_advection_NS' + str(nsweeps) + '.pkl'
    file = open(fname, 'rb')
    results_adv = pickle.load(file)
    file.close()

    xvalues_diff = sorted(list(results_diff.keys()))
    niter_diff = []
    for x in xvalues_diff:
        niter_diff.append(results_diff[x])

    xvalues_adv = sorted(list(results_adv.keys()))
    niter_adv = []
    for x in xvalues_adv:
        niter_adv.append(results_adv[x])

    # set up plotting parameters
    params = {'legend.fontsize': 20,
              'figure.figsize': (12, 8),
              'axes.labelsize': 20,
              'axes.titlesize': 20,
              'xtick.labelsize': 16,
              'ytick.labelsize': 16,
              'lines.linewidth': 3
              }
    plt.rcParams.update(params)

    # set up figure
    plt.figure()
    plt.xlabel(r'$\mu$')
    plt.ylabel('no. of iterations')
    plt.xlim(min(xvalues_diff + xvalues_adv) / 10.0, max(xvalues_diff + xvalues_adv) * 10.0)
    plt.ylim(min(niter_diff + niter_adv) - 1, max(niter_diff + niter_adv) + 1)
    plt.grid()

    # plot
    plt.semilogx(xvalues_diff, niter_diff, 'r-', marker='s', markersize=10, label='diffusion')
    plt.semilogx(xvalues_adv, niter_adv, 'b-', marker='o', markersize=10, label='advection')

    plt.legend(loc=1, ncol=1, numpoints=1)

    # plt.show()
    # save plot, beautify
    fname = 'data/conv_test_niter_NS' + str(nsweeps) + '.png'
    plt.savefig(fname, rasterized=True, bbox_inches='tight')

    assert os.path.isfile(fname), 'ERROR: plotting did not create file' 

def peaking(GdB, fc, Q=3.5, fs=44100.):
    """
    A second-order peaking filter having GdB gain at fc and approximately
    and 0 dB otherwise.
    
    The filter coefficients returns correspond to a biquadratic system function
    containing five parameters.
    
    Parameters
    ----------
    GdB : Lowpass gain in dB
    fc : Center frequency in Hz
    Q : Filter Q which is inversely proportional to bandwidth
    fs : Sampling frquency in Hz
    
    Returns
    -------
    b : ndarray containing the numerator filter coefficients
    a : ndarray containing the denominator filter coefficients
    
    Examples
    --------
    >>> import matplotlib.pyplot as plt
    >>> import numpy as np
    >>> from sk_dsp_comm.sigsys import peaking
    >>> from scipy import signal
    >>> b,a = peaking(2.0,500)
    >>> f = np.logspace(1,5,400)
    >>> w,H = signal.freqz(b,a,2*np.pi*f/44100)
    >>> plt.semilogx(f,20*np.log10(abs(H)))
    >>> plt.ylabel("Power Spectral Density (dB)")
    >>> plt.xlabel("Frequency (Hz)")
    >>> plt.show()

    >>> b,a = peaking(-5.0,500,4)
    >>> w,H = signal.freqz(b,a,2*np.pi*f/44100)
    >>> plt.semilogx(f,20*np.log10(abs(H)))
    >>> plt.ylabel("Power Spectral Density (dB)")
    >>> plt.xlabel("Frequency (Hz)")
    """
    mu = 10**(GdB/20.)
    kq = 4/(1 + mu)*np.tan(2*np.pi*fc/fs/(2*Q))
    Cpk = (1 + kq *mu)/(1 + kq)
    b1 = -2*np.cos(2*np.pi*fc/fs)/(1 + kq*mu)
    b2 = (1 - kq*mu)/(1 + kq*mu)
    a1 = -2*np.cos(2*np.pi*fc/fs)/(1 + kq)
    a2 = (1 - kq)/(1 + kq)
    b = Cpk*np.array([1, b1, b2])
    a = np.array([1, a1, a2])
    return b,a 

def ten_band_eq_resp(GdB,Q=3.5):
    """
    Create a frequency response magnitude plot in dB of a ten band equalizer
    using a semilogplot (semilogx()) type plot
    
    
    Parameters
    ----------
    GdB : Gain vector for 10 peaking filters [G0,...,G9]
    Q : Quality factor for each peaking filter (default 3.5)
    
    Returns
    -------
    Nothing : two plots are created
    
    Examples
    --------
    >>> import matplotlib.pyplot as plt
    >>> from sk_dsp_comm import sigsys as ss
    >>> ss.ten_band_eq_resp([0,10.0,0,0,-1,0,5,0,-4,0])
    >>> plt.show()
    """
    fs = 44100.0 # Hz
    NB = len(GdB)
    if not NB == 10:
        raise ValueError("GdB length not equal to ten")
    Fc = 31.25*2**np.arange(NB)
    B = np.zeros((NB,3));
    A = np.zeros((NB,3));
    
    # Create matrix of cascade coefficients
    for k in range(NB):
        b,a = peaking(GdB[k],Fc[k],Q,fs)
        B[k,:] = b
        A[k,:] = a
    # Create the cascade frequency response
    F = np.logspace(1,np.log10(20e3),1000)
    H = np.ones(len(F))*np.complex(1.0,0.0)
    for k in range(NB):
       w,Htemp = signal.freqz(B[k,:],A[k,:],2*np.pi*F/fs)
       H *= Htemp
    plt.figure(figsize=(6,4))
    plt.subplot(211)
    plt.semilogx(F,20*np.log10(abs(H)))
    plt.axis([10, fs/2, -12, 12])
    plt.grid()
    plt.title('Ten-Band Equalizer Frequency Response')
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Gain (dB)')
    plt.subplot(212)
    plt.stem(np.arange(NB),GdB,'b','bs')
    #plt.bar(np.arange(NB)-.1,GdB,0.2)
    plt.axis([0, NB-1, -12, 12])
    plt.xlabel('Equalizer Band Number')
    plt.ylabel('Gain Set (dB)')
    plt.grid() 

def plot(color=plt.get_cmap('tab10').colors[0]):
    n = FLAGS.num_test_pts
    hyperparams.announceEval()
    eval_pairs = hyperparams.getEvalDataGen()
    # Very special case lfnet:
    if FLAGS.baseline == 'lfnet':

        forward_pass_dict = baselines.parseLFNetOuts(eval_pairs, n)
        fps = []
        for pair_i in range(len(eval_pairs)):
            pair = eval_pairs[pair_i]
            folder, a, b = pair.name().split(' ')
            forward_passes = [forward_pass_dict['%s%s' % (folder, i)]
                              for i in [a, b]]
            matched_indices = system.match(forward_passes)
            fps.append([forward_passes[i][matched_indices[i]] for i in [0, 1]])
    else:
        pair_outs = cache_forward_pass.loadOrCalculateOuts()
        if FLAGS.num_scales > 1 and FLAGS.baseline == '':
            raise NotImplementedError
        else:
            fps = []
            full_fps = []
            for pair in pair_outs:
                reduced = [system.forwardPassFromHicklable(i).reducedTo(n)
                           for i in pair]
                full_fps.append(reduced)
                matched_indices = system.match(reduced)
                fps.append([reduced[i][matched_indices[i]] for i in [0, 1]])

    pairs_fps = zip(eval_pairs, fps)
    masks_errs = [evaluate.p3pMaskAndError(pair_fps[0], pair_fps[1])
                  for pair_fps in pairs_fps]

    if FLAGS.baseline == '':
        for mask_e, pair_fps, full in zip(masks_errs, pairs_fps, full_fps):
            mask, rerr, terr = mask_e
            pair, fps = pair_fps
            evaluate.renderMatching(pair, full, fps, mask)

    ninl = np.array([np.count_nonzero(i[0]) for i in masks_errs])
    rerrs = np.array([i[1] for i in masks_errs])
    rerrs[ninl < 10] = np.inf
    terrs = np.array([i[2] for i in masks_errs])
    terrs[ninl < 10] = np.inf

    if FLAGS.baseline != 'sift':
        rlabel, tlabel = hyperparams.label(), None
    else:
        rlabel, tlabel = '%s, rot.' % hyperparams.label(), \
                         '%s, transl.' % hyperparams.label()

    x, y = evaluate.lessThanCurve(rerrs)
    plt.semilogx(x, y, '--', color=color, label=rlabel)
    x, y = evaluate.lessThanCurve(terrs)
    plt.semilogx(x, y, ':', color=color, label=tlabel) 

def plot_cdf(self, n=(0, 1, 2, 3), d_min=None, d_max=None, pts=500,
                 method='logarithmic'):   # pragma: no cover
        r'''Plot the cumulative distribution function of the particle size 
        distribution. The plotted range can be specified using `d_min` and 
        `d_max`, or estimated automatically. One or more order can be plotted,
        by providing an iterable of ints as the value of `n` or just one int.
                
        Parameters
        ----------
        n : tuple(int) or int, optional
            None (for the `order` specified when the distribution was created),
            0 (number), 1 (length), 2 (area), 3 (volume/mass),
            or any integer; as many as desired may be specified, [-]
        d_min : float, optional
            Lower particle size diameter, [m]
        d_max : float, optional
            Upper particle size diameter, [m]
        pts : int, optional
            The number of points for values to be calculated, [-]
        method : str, optional
            Either 'linear', 'logarithmic', a Renard number like 'R10' or 'R5' 
            or'R2.5', or one of the sieve standards 'ISO 3310-1 R40/3', 
            'ISO 3310-1 R20', 'ISO 3310-1 R20/3', 'ISO 3310-1', 
            'ISO 3310-1 R10', 'ASTM E11', [-]
        '''
        try:
            import matplotlib.pyplot as plt
        except:  # pragma: no cover
            raise Exception(NO_MATPLOTLIB_MSG)

        ds = self.ds_discrete(d_min=d_min, d_max=d_max, pts=pts, method=method)
        try:
            for ni in n:
                cdfs = self.cdf_discrete(ds=ds, n=ni)
                plt.semilogx(ds, cdfs, label=_label_distribution_n(ni))
        except:
            cdfs = self.cdf_discrete(ds=ds, n=n)
            plt.semilogx(ds, cdfs, label=_label_distribution_n(n))
        if self.points:
            plt.plot(self.ds, self.fraction_cdf, '+', label='Volume/Mass points')
            
            if hasattr(self, 'area_fractions'):
                plt.plot(self.ds, cumsum(self.area_fractions), '+', label='Area points')
            if hasattr(self, 'length_fractions'):
                plt.plot(self.ds, cumsum(self.length_fractions), '+', label='Length points')
            if hasattr(self, 'number_fractions'):
                plt.plot(self.ds, cumsum(self.number_fractions), '+', label='Number points')
                
        plt.ylabel('Cumulative density function, [-]')
        plt.xlabel('Particle diameter, [m]')
        plt.title('Cumulative density function of %s distribution with '
                  'parameters %s' %(self.name, self.parameters))
        plt.legend()
        plt.show() 

def plot_pdf(self, n=(0, 1, 2, 3), d_min=None, d_max=None, pts=500,
                 normalized=False, method='linear'):  # pragma: no cover
        r'''Plot the probability density function of the particle size 
        distribution. The plotted range can be specified using `d_min` and 
        `d_max`, or estimated automatically. One or more order can be plotted,
        by providing an iterable of ints as the value of `n` or just one int.
                
        Parameters
        ----------
        n : tuple(int) or int, optional
            None (for the `order` specified when the distribution was created),
            0 (number), 1 (length), 2 (area), 3 (volume/mass),
            or any integer; as many as desired may be specified, [-]
        d_min : float, optional
            Lower particle size diameter, [m]
        d_max : float, optional
            Upper particle size diameter, [m]
        pts : int, optional
            The number of points for values to be calculated, [-]
        normalized : bool, optional
            Whether to display the actual probability density function, which
            may have a huge magnitude - or to divide each point by the sum
            of all the points. Doing this is a common practice, but the values
            at each point are dependent on the number of points being plotted,
            and the distribution of the points;
            [-]
        method : str, optional
            Either 'linear', 'logarithmic', a Renard number like 'R10' or 'R5' 
            or'R2.5', or one of the sieve standards 'ISO 3310-1 R40/3', 
            'ISO 3310-1 R20', 'ISO 3310-1 R20/3', 'ISO 3310-1', 
            'ISO 3310-1 R10', 'ASTM E11', [-]
        '''
        try:
            import matplotlib.pyplot as plt
        except:  # pragma: no cover
            raise Exception(NO_MATPLOTLIB_MSG)           
        ds = self.ds_discrete(d_min=d_min, d_max=d_max, pts=pts, method=method)
        try:
            for ni in n:
                fractions = [self.pdf(d, n=ni) for d in ds]
                if normalized:
                    fractions = normalize(fractions)
                plt.semilogx(ds, fractions, label=_label_distribution_n(ni))
        except Exception as e:
            fractions = [self.pdf(d, n=n) for d in ds]
            if normalized:
                fractions = normalize(fractions)
            plt.semilogx(ds, fractions, label=_label_distribution_n(n))
        plt.ylabel('Probability density function, [-]')
        plt.xlabel('Particle diameter, [m]')
        plt.title('Probability density function of %s distribution with '
                  'parameters %s' %(self.name, self.parameters))
        plt.legend()
        plt.show()
        return fractions 

def plot2():
    x = []
    y = []

    if FLAGS.baseline == '':
        ochan = FLAGS.chan
        for chan in [64, 128, 256]:
            FLAGS.chan = chan
            repsize, acc = repsizeAcc()
            x.append(repsize)
            y.append(acc)
        FLAGS.chan = ochan
    elif FLAGS.baseline == 'sift_pca':
        for pca_dim in [3, 5, 10]:
            FLAGS.pca_dim = pca_dim
            for num in [64, 128, 256, 500, 1000]:
                FLAGS.baseline_num_ips = num
                repsize, acc = repsizeAcc()
                x.append(repsize)
                y.append(acc)
            x.append(None)
            y.append(None)
    elif FLAGS.baseline == 'sift_pq':
        for m, k, i in [(4, 16, 10), (2, 256, 10), (2, 16, 10), (8, 4, 10),
                        (4, 4, 10), (1, 256, 10)]:
            FLAGS.pq_m = m
            FLAGS.pq_k = k
            FLAGS.pq_n_init = i
            for num in [64, 128, 256, 500, 1000]:
                FLAGS.baseline_num_ips = num
                repsize, acc = repsizeAcc()
                x.append(repsize)
                y.append(acc)
            x.append(None)
            y.append(None)
    else:
        for num in [64, 128, 256, 500, 1000]:
            FLAGS.baseline_num_ips = num
            repsize, acc = repsizeAcc()
            x.append(repsize)
            y.append(acc)

    plt.semilogx(x, y, '--o',
                 label='%s' % (label()),
                 color=hyperparams.methodColor())