Python matplotlib.pyplot.loglog() Examples

def plot_convergence_curve(data, info, dirname):
    # plot the convergence curve
    eps = 1e-6

    # compute the best pobj over all methods
    best_pobj = np.min([np.min(r['pobj']) for _, r in data])

    fig = plt.figure("convergence", figsize=(12, 12))
    plt.ticklabel_format(style='sci', axis='y', scilimits=(0, 0))

    color_cycle = itertools.cycle(COLORS)
    for (args, res), color in zip(data, color_cycle):
        times = list(np.cumsum(res['times']))
        plt.loglog(
            times, (res['pobj'] - best_pobj) / best_pobj + eps, '.-',
            label=get_label(info['grid_key'], args), color=color,
            linewidth=2)
    plt.xlabel('Time (s)', fontsize=24)
    plt.ylabel('Objective value', fontsize=24)
    ncol = int(np.ceil(len(data) / 10))
    plt.legend(ncol=ncol, fontsize=24)

    plt.gca().tick_params(axis='x', which='both', bottom=False, top=False)
    plt.gca().tick_params(axis='y', which='both', left=False, right=False)
    plt.tight_layout()
    plt.grid(True)
    figname = "{}/convergence.png".format(dirname)
    fig.savefig(figname, dpi=150) 

def plot_wcc_distribution(_g, _plot_img):
    """Plot weakly connected components size distributions
    :param _g: Transaction graph
    :param _plot_img: WCC size distribution image (log-log plot)
    :return:
    """
    all_wcc = nx.weakly_connected_components(_g)
    wcc_sizes = Counter([len(wcc) for wcc in all_wcc])
    size_seq = sorted(wcc_sizes.keys())
    size_hist = [wcc_sizes[x] for x in size_seq]

    plt.figure(figsize=(16, 12))
    plt.clf()
    plt.loglog(size_seq, size_hist, 'ro-')
    plt.title("WCC Size Distribution")
    plt.xlabel("Size")
    plt.ylabel("Number of WCCs")
    plt.savefig(_plot_img) 

def plot_powA(k,powNbody,powLPT,powRecon,LxN,RxN,label,c_i):
    c = plt.rcParams['axes.prop_cycle'].by_key()['color']

    ax1.loglog(k,powLPT,color = c[c_i],ls='--')
    ax1.loglog(k,powRecon,color=c[c_i],ls=':')
    l0=ax1.loglog(k,powNbody,color=c[c_i],ls='-',label=label)

    l1 = ax2.plot(k, powLPT/powNbody,color = c[c_i],ls='--')
    l2 = ax2.plot(k, powRecon/powNbody,label = label,color = c[c_i],ls=':')
    ax2.axhline(y=1, color='k', linestyle='--')


    ax3.loglog(k, 1-(LxN/np.sqrt(powLPT*powNbody))**2,color = c[c_i],ls='--')
    ax3.loglog(k, 1-(RxN/np.sqrt(powRecon*powNbody))**2,color = c[c_i],ls=':')
    return l0,l1,l2
#----------plot residual -----------------# 

def plot_distance_comparisons(self, distance1, distance2, logaxis=False,
        figure_size=(7, 5), filename=None, filetype="png", dpi=300):
        """
        Creates a plot comparing different distance metrics for the 
        specific rupture and target site combination
        """
        xdist = self._calculate_distance(distance1)       
        ydist = self._calculate_distance(distance2)       
        plt.figure(figsize=figure_size)

        if logaxis:
            plt.loglog(xdist, ydist, color='b', marker='o', linestyle='None')
        else:
            plt.plot(xdist, ydist, color='b', marker='o', linestyle='None')
        
        plt.xlabel("%s (km)" % distance1, size='medium')
        plt.ylabel("%s (km)" % distance2, size='medium')
        plt.title('Rupture: M=%6.1f, Dip=%3.0f, Ztor=%4.1f, Aspect=%5.2f'
                   % (self.magnitude, self.dip, self.ztor, self.aspect))
        _save_image(filename, filetype, dpi)
        plt.show() 

def fig_memory_usage():

    # FAST memory
    x = [1,3,7,14,30,90,180]
    y_fast = [0.653,1.44,2.94,4.97,9.05,19.9,35.2]
    # ConvNetQuake
    y_convnet = [6.8*1e-5]*7
    # Create figure
    plt.loglog(x,y_fast,"o-")
    plt.hold('on')
    plt.loglog(x,y_convnet,"o-")
    # plot markers
    plt.loglog(x,[1e-5,1e-5,1e-5,1e-5,1e-5,1e-5,1e-5],'o')
    plt.ylabel("Memory usage (GB)")
    plt.xlabel("Continous data duration (days)")
    plt.xlim(1,180)
    plt.grid("on")
    plt.savefig("./figures/memoryusage.eps")
    plt.close() 

def PlotOrthogonalResiduals(ModelX, ModelY, DataX, DataY):

    """
    """
    
    # setup the figure
    Fig = CreateFigure(AspectRatio=1.2)
    
    # Get residuals
    Residuals, OrthoX, OrthoY = OrthogonalResiduals(ModelX, ModelY, DataX, DataY)
    
    # plot model and data
    plt.loglog()
    plt.axis('equal')
    plt.plot(ModelX,ModelY,'k-', lw=1)
    plt.plot(DataX,DataY,'k.', ms=2)
    
    # plot orthogonals
    for i in range(0,len(DataX)):
        plt.plot([DataX[i],OrthoX[i]],[DataY[i],OrthoY[i]],'-',color=[0.5,0.5,0.5])
    
    plt.savefig(PlotDirectory+FilenamePrefix + "_ESRSOrthoResiduals.png", dpi=300) 

def example2():
    """
    Compute the GRADEV of a nonstationary white phase noise.
    """
    N=1000 # number of samples
    f = 1 # data samples per second
    s=1+5/N*np.arange(0,N)
    y=s*np.random.randn(1,N)[0,:]
    x = [xx for xx in np.linspace(1,len(y),len(y))]
    x_ax, y_ax, (err_l, err_h) , ns = allan.gradev(y,data_type='phase',rate=f,taus=x)
    plt.loglog(x_ax, y_ax,'b.',label="No gaps")
    y[int(0.4*N):int(0.6*N,)] = np.NaN  # Simulate missing data
    x_ax, y_ax, (err_l, err_h), ns = allan.gradev(y,data_type='phase',rate=f,taus=x)
    plt.loglog(x_ax, y_ax,'g.',label="With gaps")
    plt.grid()
    plt.legend()
    plt.xlabel('Tau / s')
    plt.ylabel('Overlapping Allan deviation')
    plt.show() 

def fig_run_time():
    # fast run time
    x_fast = [1,3,7,14,30,90,180]
    y_fast = [289,1.13*1e3,2.48*1e3,5.41*1e3,1.56*1e4,
              6.61*1e4,1.98*1e5]
    x_auto = [1,3]
    y_auto = [1.54*1e4, 8.06*1e5]
    x_convnet = [1,3,7,14,30]
    y_convnet = [9,27,61,144,291]
    # create figure
    plt.loglog(x_auto,y_auto,"o-")
    plt.hold('on')
    plt.loglog(x_fast[0:5],y_fast[0:5],"o-")
    plt.loglog(x_convnet,y_convnet,"o-")
    # plot x markers
    plt.loglog(x_convnet,[1e0]*len(x_convnet),'o')
    # plot y markers
    y_markers = [1,60,3600,3600*24]
    plt.plot([1]*4,y_markers,'ko')
    plt.ylabel("run time (s)")
    plt.xlabel("continous data duration (days)")
    plt.xlim(1,35)
    plt.grid("on")
    plt.savefig("./figures/runtimes.eps") 

def _scale_curve(self):
        """
        Puts the data-misfit and regularizing function values in the range
        [-10, 10].
        """
        if self.loglog:
            x, y = numpy.log(self.dnorm), numpy.log(self.mnorm)
        else:
            x, y = self.dnorm, self.mnorm

        def scale(a):
            vmin, vmax = a.min(), a.max()
            l, u = -10, 10
            return (((u - l) / (vmax - vmin)) *
                    (a - (u * vmin - l * vmax) / (u - l)))
        return scale(x), scale(y) 

def plot_cltv(time_locks):
    exponent_min = 0
    exponent_max = 3

    bin_factor = 10

    bins_log = 10**np.linspace(
        exponent_min, exponent_max, (exponent_max - exponent_min) * bin_factor + 1)

    fig, ax = plt.subplots(figsize=standard_figsize, dpi=300)

    ax.hist(time_locks, bins=bins_log)
    plt.loglog()
    ax.set_xlabel("CLTV bins [blocks]")
    ax.set_ylabel("Number of channels")
    plt.tight_layout()
    plt.show() 

def plot_base_fees(base_fees):
    exponent_min = 0
    exponent_max = 5

    bin_factor = 10

    bins_log = 10**np.linspace(
        exponent_min, exponent_max, (exponent_max - exponent_min) * bin_factor + 1)

    fig, ax = plt.subplots(figsize=standard_figsize, dpi=300)
    ax.hist(base_fees, bins=bins_log)
    ax.axvline(x=1E3, c='k', ls='--')
    plt.loglog()
    ax.set_xlabel("Base fee bins [msat]")
    ax.set_ylabel("Number of channels")
    plt.tight_layout()
    plt.show() 

def plot_fee_rates(fee_rates):
    exponent_min = -6
    exponent_max = 0

    bin_factor = 10

    bins_log = 10**np.linspace(
        exponent_min, exponent_max, (exponent_max - exponent_min) * bin_factor + 1)
    print(bins_log)
    fig, ax = plt.subplots(figsize=standard_figsize, dpi=300)
    ax.axvline(x=1E-6, c='k', ls='--')
    ax.hist(fee_rates, bins=bins_log)
    plt.loglog()
    ax.set_xlabel("Fee rate bins [sat per sat]")
    ax.set_ylabel("Number of channels")
    plt.tight_layout()
    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 __init__(self, datamisfit, regul, regul_params, loglog=True,
                 njobs=1):
        assert njobs >= 1, "njobs should be >= 1. {} given.".format(njobs)
        self.regul_params = regul_params
        self.datamisfit = datamisfit
        self.regul = regul
        self.objectives = None
        self.dnorm = None
        self.mnorm = None
        self.fit_method = None
        self.fit_args = None
        self.njobs = njobs
        self.loglog = loglog
        # Estimated parameters from the L curve
        self.corner_ = None 

def plot_loss(reg_ratio):
    n_trials = 1
    n_channels = 1
    n_times, n_atoms, n_times_atom = 100000, 10, 100

    rng = np.random.RandomState(0)
    X = rng.randn(n_trials, n_channels, n_times)
    D = rng.randn(n_atoms, n_channels, n_times_atom)

    reg = reg_ratio * get_lambda_max(X, D).max()

    results = []
    for func, max_iter in all_func:
        print(func.__name__)
        res = run_one(func, n_times, n_atoms, n_times_atom, reg, max_iter, X,
                      D)
        results.append(res)

    best = np.inf
    for res in results:
        func_name, n_times, n_atoms, n_times_atom, reg, times, pobj = res
        if pobj[-1] < best:
            best = pobj[-1]

    fig = plt.figure()
    for (func, max_iter), res in zip(all_func, results):
        style = '-' if 'cd' in func.__name__ else '--'
        func_name, n_times, n_atoms, n_times_atom, reg, times, pobj = res
        plt.loglog(times, np.array(pobj) - best, style, label=func.__name__)
    plt.legend()
    plt.xlim(1e-2, None)
    name = ('T=%s_K=%s_L=%s_reg=%.3f' % (n_times, n_atoms, n_times_atom,
                                         reg_ratio))
    plt.title(name)
    plt.xlabel('Time (s)')
    plt.ylabel('loss function')
    save_name = 'figures/bench_gcd/' + name + '.png'
    print('Saving %s' % (save_name, ))
    fig.savefig(save_name)
    # plt.show()
    plt.close(fig) 

def plot_runtimes(run_time_objects, problem_sizes, symbols):
    _plot(plt.loglog, problem_sizes, run_time_objects, symbols,
          ylabel='time $t$', loc=2, logx=True) 

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 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 plot_lcurve(self, ax=None, guides=True):
        """
        Make a plot of the data-misfit x regularization values.

        The estimated corner value is shown as a blue triangle.

        Parameters:

        * ax : matplotlib Axes
            If not ``None``, will plot the curve on this Axes instance.
        * guides : True or False
            Plot vertical and horizontal lines across the corner value.

        """
        if ax is None:
            ax = mpl.gca()
        else:
            mpl.sca(ax)
        x, y = self.dnorm, self.mnorm
        if self.loglog:
            mpl.loglog(x, y, '.-k')
        else:
            mpl.plot(x, y, '.-k')
        if guides:
            vmin, vmax = ax.get_ybound()
            mpl.vlines(x[self.corner_], vmin, vmax)
            vmin, vmax = ax.get_xbound()
            mpl.hlines(y[self.corner_], vmin, vmax)
        mpl.plot(x[self.corner_], y[self.corner_], '^b', markersize=10)
        mpl.xlabel('Data misfit(data norm)')
        mpl.ylabel('Regularization(model norm)') 

def plot_pow(k,powNbody,powLPT,powRecon,LxN,RxN,title):
    fig = plt.figure(figsize=(6,8))

    ax1 = plt.subplot2grid((4,1),(0,0),rowspan=2)
    plt.plot(k, powLPT,label = '2LPT')
    plt.plot(k, powRecon,label = 'U-Net')
    plt.plot(k, powNbody,label ='fastPM')
    plt.ylabel('P(k)')
    plt.yscale('log')
    plt.legend(loc='lower left')
    plt.title(title)
    plt.setp(ax1.get_xticklabels(),visible=False)

    ax2 = plt.subplot2grid((4,1),(2,0), rowspan = 1,sharex=ax1)
    plt.axhline(y=1, color='k', linestyle='--')
    plt.plot(k, powLPT/powNbody,label = 'LPT')
    plt.plot(k, powRecon/powNbody,label = 'Predict')
    plt.ylabel(r'$T(k)$')
    plt.setp(ax2.get_xticklabels(),visible=False)

    ax3 = plt.subplot2grid((4,1),(3,0),sharex=ax1)
    plt.loglog(k, 1-(LxN/np.sqrt(powLPT*powNbody))**2,label = 'LPTxNbody')
    plt.loglog(k, 1-(RxN/np.sqrt(powRecon*powNbody))**2,label = 'ReconxNbody')
    plt.xticks(np.round([0.06+0.01*i for i in range(0,4,2)]+ [0.1+0.1*i for i in range(0,7,2)],2),
               np.round([0.06+0.01*i for i in range(0,4,2)]+ [0.1+0.1*i for i in range(0,7,2)],2))
    plt.xticks(rotation=45)
    plt.ylabel(r'1-$r^2$')
    #plt.tight_layout() 

def plot_pancake(k,powNbody,powRecon,powInput,title):
    pos = np.intersect1d(np.where(np.nan_to_num(powNbody) > 1e-3)[0], np.where(np.nan_to_num(powRecon)> 1e-3)[0])
    plt.figure(figsize=(6,6))
    gs = GridSpec(2, 1, height_ratios=[2,1],width_ratios=[1])
    ax1 = plt.subplot(gs[0, 0])
    #plt.loglog(k[pos],powLPT[pos],'*',label='2LPT',c=c[0])
    plt.loglog(k[pos],powRecon[pos],'*',label='U-Net',c=c[1])
    plt.loglog(k[pos],powNbody[pos],'x',label='FastPM',c=c[2])
    plt.loglog(k,powInput,'^',label='1 mode input',c=c[3])
    plt.vlines(k[np.argmax(np.nan_to_num(powInput))],\
               ymin = 1e-5,\
               ymax = powInput[np.argmax(np.nan_to_num(powInput))],alpha=0.5,linestyles='dashed')
    plt.legend(loc='upper left')
    plt.title(r"displacement "+title)
    plt.ylabel(r'$P(k)}$')
    plt.ylim(ymin = 1e-5)

    plt.subplot(gs[1, 0],sharex=ax1)
    plt.loglog(k[pos],powRecon[pos]/powNbody[pos],'*',label='Recon',c=c[1])
    plt.xlabel('k '+r'[h/Mpc]')
    plt.ylim(ymin = 1e-5)
    plt.ylabel(r'$T(k)$')
    plt.axhline(y=1,color='black',ls='dashed')
    plt.setp(ax1.get_xticklabels(),visible=False)
    plt.xticks(np.round([0.06+0.01*i for i in range(0,4,2)]+ [0.1+0.1*i for i in range(0,7,2)],2),\
                   np.round([0.06+0.01*i for i in range(0,4,2)]+ [0.1+0.1*i for i in range(0,7,2)],2))
    plt.xticks(rotation=45)
    #plt.tight_layout()

#----------plot one slice demonstration--------------# 

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 test_radial_profile(self):
        r = np.logspace(start=-2, stop=2, num=100)
        c = 10
        logM = 13.
        #kappa = self.nfw.kappa(r, logM=logM, c=c)
        import matplotlib.pyplot as plt
        #plt.loglog(r, kappa)
        #plt.show()
        #assert 1 == 0 

def _fractal_dfa_plot(windows, fluctuations, dfa):
    fluctfit = 2 ** np.polyval(dfa, np.log2(windows))
    plt.loglog(windows, fluctuations, "bo")
    plt.loglog(windows, fluctfit, "r", label=r"$\alpha$ = %0.3f" % dfa[0])
    plt.title("DFA")
    plt.xlabel(r"$\log_{2}$(Window)")
    plt.ylabel(r"$\log_{2}$(Fluctuation)")
    plt.legend()
    plt.show()


# =============================================================================
#  Utils MDDFA
# ============================================================================= 

def plot_fourier_spectrum(time_series, time_step, figure_size=(7, 5),
        filename=None, filetype="png", dpi=300):
    """
    Plots the Fourier spectrum of a time series 
    """
    freq, amplitude = get_fourier_spectrum(time_series, time_step)
    plt.figure(figsize=figure_size)
    plt.loglog(freq, amplitude, 'b-')
    plt.xlabel("Frequency (Hz)", fontsize=14)
    plt.ylabel("Fourier Amplitude", fontsize=14)
    _save_image(filename, filetype, dpi)
    plt.show() 

def plot_gt_freqs(fp):
    """
    Draws a scatterplot of the empirical frequencies of the counted species
    versus their Simple Good Turing smoothed values, in rank order. Depends on
    pylab and matplotlib.
    """
    MLE = MLENGram(1, filter_punctuation=False, filter_stopwords=False)
    MLE.train(fp, encoding="utf-8-sig")
    counts = dict(MLE.counts[1])

    GT = GoodTuringNGram(1, filter_stopwords=False, filter_punctuation=False)
    GT.train(fp, encoding="utf-8-sig")

    ADD = AdditiveNGram(1, 1, filter_punctuation=False, filter_stopwords=False)
    ADD.train(fp, encoding="utf-8-sig")

    tot = float(sum(counts.values()))
    freqs = dict([(token, cnt / tot) for token, cnt in counts.items()])
    sgt_probs = dict([(tok, np.exp(GT.log_prob(tok, 1))) for tok in counts.keys()])
    as_probs = dict([(tok, np.exp(ADD.log_prob(tok, 1))) for tok in counts.keys()])

    X, Y = np.arange(len(freqs)), sorted(freqs.values(), reverse=True)
    plt.loglog(X, Y, "k+", alpha=0.25, label="MLE")

    X, Y = np.arange(len(sgt_probs)), sorted(sgt_probs.values(), reverse=True)
    plt.loglog(X, Y, "r+", alpha=0.25, label="simple Good-Turing")

    X, Y = np.arange(len(as_probs)), sorted(as_probs.values(), reverse=True)
    plt.loglog(X, Y, "b+", alpha=0.25, label="Laplace smoothing")

    plt.xlabel("Rank")
    plt.ylabel("Probability")
    plt.legend()
    plt.tight_layout()
    plt.savefig("img/rank_probs.png")
    plt.close("all") 

def test_plot(test):
    ax = plt.gca()
    ax.set_color_cycle(['b', 'g', 'r', 'c', 'm', 'y', 'k', '0.8'])
    kinds = order_kinds(sorted(args.kind or list(data[test])))
    for kind in kinds:
        kind_plot(test, kind)
    plt.ylim(ymin=9e-7)
    plt.loglog()
    plt.title(args.name + ' Performance: ' + test)
    plt.ylabel('Seconds')
    plt.xlabel('List Size')
    plt.legend(kinds, loc=2) 

def plotit(depth, a, n, res1, res2, res3, title):
    """Call `comp_appres` and plot result."""

    # Compute the three different models
    rho1, AB2 = comp_appres(depth, res1, a, n)
    rho2, _ = comp_appres(depth, res2, a, n)
    rho3, _ = comp_appres(depth, res3, a, n)

    # Create figure
    plt.figure()

    # Plot curves
    plt.loglog(AB2, rho1, label='Case 1')
    plt.plot(AB2, rho2, label='Case 2')
    plt.plot(AB2, rho3, label='Case 3')

    # Legend, labels
    plt.legend(loc='best')
    plt.title(title)
    plt.xlabel('AB/2 (m)')
    plt.ylabel(r'Apparent resistivity $\rho_a (\Omega\,$m)')

    plt.show()


###############################################################################
# Model 1: 2 layers
# ~~~~~~~~~~~~~~~~~
#
# +--------+---------------------+---------------------+
# |layer   | depth (m)           | resistivity (Ohm m) |
# +========+=====================+=====================+
# |air     | :math:`-\infty` - 0 | 2e14                |
# +--------+---------------------+---------------------+
# |layer 1 | 0 - 50              | 10                  |
# +--------+---------------------+---------------------+
# |layer 2 | 50 - :math:`\infty` | 100 / 10 / 1        |
# +--------+---------------------+---------------------+ 

def _fractal_correlation_plot(r_vals, corr, d2):
    fit = 2 ** np.polyval(d2, np.log2(r_vals))
    plt.loglog(r_vals, corr, "bo")
    plt.loglog(r_vals, fit, "r", label=r"$D2$ = %0.3f" % d2[0])
    plt.title("Correlation Dimension")
    plt.xlabel(r"$\log_{2}$(r)")
    plt.ylabel(r"$\log_{2}$(c)")
    plt.legend()
    plt.show() 

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()'''