Python matplotlib.pyplot.loglog() Examples

The following are code examples for showing how to use matplotlib.pyplot.loglog(). They are from open source Python projects. You can vote up the examples you like or vote down the ones you don't like.

Example 1
Project: DiscEvolution   Author: rbooth200   File: plot_planet_M-R.py    GNU General Public License v3.0 6 votes vote down vote up
def make_plot_planets(planets):
    for p in planets:
        plt.subplot(211)
        plt.loglog(p.R, p.M / 317.8)
        plt.subplot(212)
        plt.loglog(p.R, p.M_core)

    plt.subplot(211)
    plt.title('t_0 = {:g}yr'.format(p.t_form[0]))
    plt.ylabel('$M\,[M_J]$')
    plt.plot([0.1, 300], [1,1], 'k--')
    plt.xlim(0.1, 300)
    plt.subplot(212)
    plt.ylabel('$M_c\,[M_\oplus]$')
    plt.xlabel('$R\,[\mathrm{au}]$')
    plt.xlim(0.1, 300) 
Example 2
Project: DiscEvolution   Author: rbooth200   File: plot_planet_R-t.py    GNU General Public License v3.0 6 votes vote down vote up
def make_plot_planets(planets):
    for p in planets:
        plt.subplot(211)
        plt.loglog(p.R, p.M / 317.8)
        plt.subplot(212)
        plt.loglog(p.R, p.M_core)

    plt.subplot(211)
    plt.title('t_0 = {:g}yr'.format(p.t_form[0]))
    plt.ylabel('$M\,[M_J]$')
    plt.plot([0.1, 300], [1,1], 'k--')
    plt.xlim(0.1, 300)
    plt.subplot(212)
    plt.ylabel('$M_c\,[M_\oplus]$')
    plt.xlabel('$R\,[\mathrm{au}]$')
    plt.xlim(0.1, 300) 
Example 3
Project: AMLSim   Author: IBM   File: plot_distributions.py    Apache License 2.0 6 votes vote down vote up
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) 
Example 4
Project: POET   Author: kevin218   File: plots.py    GNU General Public License v3.0 6 votes vote down vote up
def rmsplot(event, fit, fignum, savefile=None, istitle=True, stderr=None, normfactor=None):
    if stderr == None:
        stderr = fit.stderr
    if normfactor == None:
        normfactor = 1e-6
    plt.rcParams.update({'legend.fontsize':11})
    plt.figure(fignum, figsize=(8,6))
    plt.clf()
    if istitle:
        a = plt.suptitle(event.eventname + ' Correlated Noise', size=16)
    plt.loglog(fit.binsz, fit.rms/normfactor, color='black', lw=1.5, label='Fit RMS', zorder=3)    # our noise
    plt.loglog(fit.binsz, stderr/normfactor, color='red', ls='-', lw=2, label='Std. Err.', zorder=1) # expected noise
    plt.xlim(0, fit.binsz[-1]*2)
    plt.ylim(stderr[-1]/normfactor/2., stderr[0]/normfactor*2.)
    plt.xlabel("Bin Size", fontsize=14)
    plt.ylabel("RMS (ppm)", fontsize=14)
    plt.xticks(size=12)
    plt.yticks(size=12)
    plt.legend()
    if savefile != None:
        plt.savefig(savefile)
    return

# Plot RMS vs. bin size (with RMS uncertainties) looking for time-correlated noise 
Example 5
Project: POET   Author: kevin218   File: plots.py    GNU General Public License v3.0 6 votes vote down vote up
def rmsploterr(event, fit, fignum, savefile=None, istitle=True, stderr=None, normfactor=None):
    if stderr == None:
        stderr = fit.stderr
    if normfactor == None:
        normfactor = stderr[0]
    plt.rcParams.update({'legend.fontsize':11})
    plt.figure(fignum, figsize=(8,6))
    plt.clf()
    if istitle:
        a = plt.suptitle(event.eventname + ' Correlated Noise', size=16)
    plt.loglog(fit.binsz, fit.rms/normfactor, color='black', lw=1.5, label='Fit RMS')    # our noise
    plt.loglog(fit.binsz, stderr/normfactor, color='red', ls='-', lw=2, label='Std. Err.') # expected noise
    plt.xlim(0, fit.binsz[-1]*2)
    plt.ylim(stderr[-1]/normfactor/2., stderr[0]/normfactor*2.)
    plt.xlabel("Bin Size", fontsize=14)
    plt.ylabel("Normalized RMS", fontsize=14)
    plt.xticks(size=12)
    plt.yticks(size=12)
    plt.legend()
    if savefile != None:
        plt.savefig(savefile)
    return 
Example 6
Project: LSDMappingTools   Author: LSDtopotools   File: plot_hillslope_morphology.py    MIT License 6 votes vote down vote up
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) 
Example 7
Project: ves_ajoros   Author: ajoros   File: VESinverse_vectorized.py    GNU Lesser General Public License v3.0 6 votes vote down vote up
def error(): # simple rms error calc
    sumerror = 0.
    #pltanswer = [0]*64
    spline(m, one30, one30, asavl, rl, y2) # So this calculates the predicted fit?
    # and essentially operates on the list in place?
    for i in range(1, ndat): # So you always skip the value 0? due to -inf returns?
        ans = splint(m, adatl[i], asavl, rl, y2) # Then this calulates error?
        sumerror = sumerror + (rdatl[i] - ans) * (rdatl[i] - ans)
        #print(i,sum1,rdat[i],rdatl[i],ans)
        pltanswerl[i] = ans
        pltanswer[i] = np.power(10, ans)
    rms = np.sqrt(sumerror / (ndat - 1))

    # check the spline routine
#    for i in range(1,m+1,1):
#        anstest = splint(m, asavl[i],asavl,rl,y2)
#        print( asavl[i], rl[i], anstest)
    #print(' rms  =  ', rms)
# if you erally want to get a good idea of all perdictions from Montecarlo
# perform the following plot (caution - change iter to a smaller number)
    #plt.loglog(adat[1:ndat],pltanswer[1:ndat])
    return rms 
Example 8
Project: ms-analysis   Author: tlnagy   File: run_analysis.py    Mozilla Public License 2.0 6 votes vote down vote up
def plot_scatter_intensities(intensity_cols, experiment1, experiment2="Control_WCL", path=""):
    """
    Create a scatter plot of the per-protein raw intensities between
    two experiments. Experiment names should be the names of the columns
    minus the "Intensity " part.
    """
    if "Intensity %s"%experiment1 not in intensity_cols.columns:
        raise ValueError("Supplied column name not in dataframe")
    f = plt.figure(figsize=(3, 3))
    plt.loglog()
    tmp_prot_groups = intensity_cols.loc[:, ["Intensity %s"%experiment1, "Intensity %s"%experiment2]]
    sns.regplot(tmp_prot_groups.ix[:, 1], tmp_prot_groups.ix[:, 0], fit_reg=False)
    xs = np.linspace(tmp_prot_groups.min()[0], tmp_prot_groups.max()[0])
    plt.plot(xs, tmp_prot_groups.sum()[0]/tmp_prot_groups.sum()[1]*xs, "r", alpha=0.4)
    sns.axlabel("%s Raw Intensities"%experiment1.replace("_", " "),
                "%s Raw Intensities"%experiment2.replace("_", " "))
    sns.despine()
    plt.tight_layout()
    plt.savefig(os.path.join(path, "%sv%s_intensities.png"%(experiment1, experiment2)), dpi=300) 
Example 9
Project: poi   Author: jchluo   File: powerlaw.py    Apache License 2.0 6 votes vote down vote up
def plot(self, filename=None, marker='+', color='blue'):
        """Plot Power Law picture.
        filename: if assgin, picture will write to file,
                not show on screen, else show on screen.
        marker  : see http://matplotlib.org/api/markers_api.html
        color   : point color
        """
        plt.loglog(*self.points, linestyle='None', marker=marker, markeredgecolor=color)
        if self.line_ready:
            x = self.points[0]
            x = np.linspace(x[1], max(x), 10)
            y = [self.prob(i) for i in x]
            plt.loglog(x, y, linestyle='-')

        if filename is not None:
            plt.savefig(filename)
            log.debug("plot %s ok." % filename)
        else:
            plt.show() 
Example 10
Project: Emergence   Author: LennonLab   File: macroecotools.py    MIT License 6 votes vote down vote up
def plot_SARs(list_of_A_and_S):
    """Plot multiple SARs on a single plot. 
    
    Input: a list of lists, each sublist contains one vector for S and one vector for A.
    Output: a graph with SARs plotted on log-log scale, with colors spanning the spectrum.
    
    """
    N = len(list_of_A_and_S)
    HSV_tuples = [(x * 1.0 / N, 0.5, 0.5) for x in range(N)]
    RGB_tuples = map(lambda x: colorsys.hsv_to_rgb(*x), HSV_tuples)
    for i in range(len(list_of_A_and_S)):
        sublist = list_of_A_and_S[i]
        plt.loglog(sublist[0], sublist[1], color = RGB_tuples[i])
    plt.hold(False)
    plt.xlabel('Area')
    plt.ylabel('Richness') 
Example 11
Project: gamma_limits_sensitivity   Author: mahnen   File: __init__.py    MIT License 6 votes vote down vote up
def plot_effective_area(a_eff_interpol, style='k', label=''):
    '''
    fill a plot with the effective energy from the supplied
    interpolated data
    '''
    start = a_eff_interpol.x.min()
    stop = a_eff_interpol.x.max()
    samples = 1000

    energy_samples = np.linspace(start, stop, samples)
    area_samples = np.array([
        a_eff_interpol(energy)
        for energy
        in energy_samples
        ])
    plt.plot(np.power(10, energy_samples), area_samples/10000.,
             style,
             label=label)

    plt.loglog()
    plt.title('Effective Area')
    plt.xlabel('Energy / TeV')
    plt.ylabel('A$_{eff}$ / m$^2$')
    return 
Example 12
Project: gamma_limits_sensitivity   Author: mahnen   File: __init__.py    MIT License 6 votes vote down vote up
def plot_power_law(
        f_0,
        gamma,
        e_0,
        energy_range,
        fmt='k:',
        label='',
        alpha_plot=0.7
        ):
    '''
    This function generates a power law plot in
    the current figure
    '''
    e_x = 10**np.arange(
        np.log10(energy_range[0]),
        np.log10(energy_range[1])+0.05,
        0.05)
    e_y = power_law(e_x, f_0, gamma, e_0=e_0)

    plt.plot(e_x, e_y, fmt, label=label, alpha=alpha_plot)
    plt.loglog()

    plt.xlabel("E / TeV")
    plt.ylabel("dN/dE / [(cm$^2$ s TeV)$^{-1}$]") 
Example 13
Project: maxent   Author: TRIQS   File: plot_utils.py    GNU General Public License v3.0 6 votes vote down vote up
def _plotter(x, y, label=None, x_label=None, y_label=None,
             log_x=False, log_y=False, **kwargs):
    """ actually plotting a curve

    a small wrapper over matplotlib"""

    plot_command = plt.plot
    if log_x and log_y:
        plot_command = plt.loglog
    elif log_x:
        plot_command = plt.semilogx
    elif log_y:
        plot_command = plt.semilogy

    if np.any(np.iscomplex(y)):
        plot_command(x, y.real,
                     label='Re ' + label if label is not None else None)
        plot_command(x, y.imag,
                     label='Im ' + label if label is not None else None)
    else:
        plot_command(x, y, label=label)
    plt.xlabel(x_label)
    plt.ylabel(y_label) 
Example 14
Project: ANN-MD   Author: haakonvt   File: plot_tools.py    MIT License 6 votes vote down vote up
def plotTestVsTrainLoss(save_dir, list_of_rmse_train, list_of_rmse_test):
    if not list_of_rmse_test or not list_of_rmse_train:
        """
        No input gotten (or not enough), must read from file
        """
        list_of_rmse_test  = np.loadtxt(save_dir + "/testRMSE.txt")
        list_of_rmse_train = np.loadtxt(save_dir + "/trainRMSE.txt")

    plt.subplot(3,1,1)
    xTest_for_plot = np.linspace(0,1,len(list_of_rmse_test))
    xTrain_for_plot = np.linspace(0,1,len(list_of_rmse_train))
    plt.plot(xTrain_for_plot, list_of_rmse_train, label="train")
    plt.plot(xTest_for_plot, list_of_rmse_test, label="test") #, lw=2.0)
    plt.subplot(3,1,2)
    plt.semilogy(xTrain_for_plot, list_of_rmse_train, label="train")
    plt.semilogy(xTest_for_plot, list_of_rmse_test, label="test") #, lw=2.0)
    plt.subplot(3,1,3)
    plt.semilogy(xTrain_for_plot, list_of_rmse_train, label="train")
    plt.loglog(xTest_for_plot, list_of_rmse_test, label="test") #, lw=2.0)
    plt.savefig(save_dir+"/RMSE_evo.pdf")
    plt.show() 
Example 15
Project: flow   Author: nschloe   File: test_stokes.py    MIT License 6 votes vote down vote up
def show_errors(hmax, u_errors, p_errors):
    # plot order indicators
    for order in range(5):
        plt.loglog(
                [hmax[0], hmax[-1]],
                [u_errors[0], u_errors[0] * (hmax[-1] / hmax[0])**order],
                color='0.7'
                )

    plt.loglog(hmax, u_errors, linestyle='-', marker='.', label='||u - uh||')
    plt.loglog(hmax, p_errors, linestyle='-', marker='.', label='||p - ph||')

    plt.xlabel('hmax')
    plt.legend()
    plt.show()
    return 
Example 16
Project: ConvNetQuake   Author: tperol   File: fig_comparison.py    MIT License 6 votes vote down vote up
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() 
Example 17
Project: ConvNetQuake   Author: tperol   File: fig_comparison.py    MIT License 6 votes vote down vote up
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") 
Example 18
Project: rel_3d_pose   Author: matteorr   File: viz.py    MIT License 6 votes vote down vote up
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') 
Example 19
Project: extrapolation   Author: numerical-mathematics   File: basic_methods.py    MIT License 6 votes vote down vote up
def test_convergence(f, exact, t0, tf, y0, method, order, graph_title):
    hs = np.asarray([2**(-k) for k in range(1, 16)])
    err = np.zeros(len(hs))

    for i in range(len(hs)):
        ts, ys = method(f, y0, t0, tf, hs[i])
        err[i] = abs(ys[-1] - exact(tf))

    plt.hold('true')
    method_err,  = plt.loglog(hs, err, 's-')
    order_line,  = plt.loglog(hs, (hs**order)*(err[5]/hs[5]))
    plt.legend([method_err, order_line], ['Method Error', 'order ' + str(order)], loc=2)
    plt.title(graph_title)
    plt.ylabel('||error||')
    plt.xlabel('time step size')

    plt.show()

# Convergence Tests 
Example 20
Project: extrapolation   Author: numerical-mathematics   File: compare_test.py    MIT License 6 votes vote down vote up
def plot_results(methods,problem_name):
    import matplotlib
    matplotlib.use('agg')
    import matplotlib.pyplot as plt
    plt.hold('true')

    extrap, dopri5, dop853 = methods

    for method in [extrap, dopri5, dop853]:
        method['line'],   = plt.loglog(method['yerr'], method['runtime'], "s-")
    plt.legend([extrap['line'], dopri5['line'], dop853['line']], ["ParEx", "DOPRI5 (scipy)", "DOP853 (scipy)"], loc=1)
    plt.xlabel('Error')
    plt.ylabel('Wall clock time (seconds)')
    plt.title(problem_name)
    plt.show()
    plt.savefig('images/' + problem_name + '_err_vs_time.png')
    plt.close()


###############################################################
###################### TEST PROBLEMS ##########################
###############################################################

###### N-Body Problem ###### 
Example 21
Project: DiscEvolution   Author: rbooth200   File: run_model.py    GNU General Public License v3.0 5 votes vote down vote up
def _plot_grid(model):
    grid = model.disc.grid 

    try:
        eps = model.disc.dust_frac.sum(0)
        plt.subplot(222)
        plt.loglog(grid.Rc, eps)
        plt.xlabel('$R$')
        plt.ylabel('$\epsilon$')
        plt.ylim(ymin=1e-4)
        plt.subplot(223)
        plt.loglog(grid.Rc, model.disc.Stokes()[1])
        plt.xlabel('$R$')
        plt.ylabel('$St$')
        plt.subplot(224)
        plt.loglog(grid.Rc, model.disc.grain_size[1])
        plt.xlabel('$R$') 
        plt.ylabel('$a\,[\mathrm{cm}]$')

        plt.subplot(221)
        l, = plt.loglog(grid.Rc, model.disc.Sigma_D.sum(0), '--')
        c = l.get_color()
    except AttributeError:
        c = None

    plt.loglog(grid.Rc, model.disc.Sigma_G, c=c)
    plt.xlabel('$R$')
    plt.ylabel('$\Sigma_\mathrm{G, D}$')
    plt.ylim(ymin=1e-5) 
Example 22
Project: EIBalanceLSM   Author: zeqianli   File: lsm.py    MIT License 5 votes vote down vote up
def plot_gap_distribution(fire_collect, dt=0.1, title="gap distribution"):
    """Plot interval distribution between spikes. 
    
    SOC shows a power-law distribution.
    """
    have_spike = ~np.any(fire_collect, axis=0)
    ind_spike = np.where(have_spike)[0]
    gaps = ind_spike[1:] - ind_spike[:-1]
    gaps = gaps[gaps > 1]

    plt.figure()
    plt.title(title)

    if len(gaps) == 0:
        print("Error: len(gaps)==0")
    else:
        ps, bins = np.histogram(gaps, bins=50, density=True)
        xs = (bins[:-1] + bins[1:]) / 2
        _arg = ps > 0
        xs = xs[_arg]
        ys = ps[_arg]
        print(xs, ys)
        plt.scatter(xs, ys)
        plt.loglog()
        plt.ylim(ymax=1.2, ymin=min(ys) * 0.5)
        plt.xlim(xmin=1.5, xmax=max(xs * 1.5)) 
Example 23
Project: fastmat   Author: EMS-TU-Ilmenau   File: benchmark.py    Apache License 2.0 5 votes vote down vote up
def plotResult(self, nameResult, outPath='', addVersionTag=True,
                   strFormat='%.6e', chrDelimiter=',', chrNewline=os.linesep):
        '''
        Print data table to csv file. The result is organized as a list of
        dicts containing information in key=value pairs. The header will be
        extracted from the keys of the result data dictionaries. If the file
        does not exist it will be created.
        '''
        if nameResult not in self.results:
            raise ValueError("No benchmark result of such name found")

        result=self.results[nameResult]
        arrResult=result[BENCH.RESULT]

        filename=self.getFilename(nameResult, path=outPath,
                                  addVersionTag=addVersionTag)

        filename = filename.replace('csv', 'png')

        # force existance of output file
        dirname=os.path.dirname(filename)
        if dirname != '' and not os.path.exists(dirname):
            os.makedirs(dirname)

        import matplotlib.pyplot as plt

        plt.loglog(arrResult[:, 0], arrResult[:, 1:])

        plt.savefig(
            filename,
            dpi=300,
            transparent=True,
            bbox_inches='tight'
        )

        return filename 
Example 24
Project: ves_ajoros   Author: ajoros   File: VESinverse_vectorized.py    GNU Lesser General Public License v3.0 5 votes vote down vote up
def error(): # simple rms error calc
    sumerror = 0.
    #pltanswer = [0]*64
    spline(m, one30, one30, asavl, rl, y2) # So this calculates the predicted fit?
    # and essentially operates on the list in place?
    for i in range(1, ndat): # So you always skip the value 0? due to -inf returns?
        ans = splint(m, adatl[i], asavl, rl, y2) # Then this calulates error?
        sumerror = sumerror + (rdatl[i] - ans) * (rdatl[i] - ans)
        #print(i,sum1,rdat[i],rdatl[i],ans)
        pltanswerl[i] = ans
        pltanswer[i] = np.power(10, ans)
    rms = np.sqrt(sumerror / (ndat - 1))

    # check the spline routine
#    for i in range(1,m+1,1):
#        anstest = splint(m, asavl[i],asavl,rl,y2)
#        print( asavl[i], rl[i], anstest)
    #print(' rms  =  ', rms)
# if you erally want to get a good idea of all perdictions from Montecarlo
# perform the following plot (caution - change iter to a smaller number)
    #plt.loglog(adat[1:ndat],pltanswer[1:ndat])
    return rms
# my code to do a spline fit to predicted data at the nice spacing of Ghosh
# use splint to determine the spline interpolated prediction at the
# spacing where the measured resistivity was taken - to compare observation
# to prediction 
Example 25
Project: Nyles   Author: pvthinker   File: timing.py    MIT License 5 votes vote down vote up
def analyze_timing(path):

    mpl.rcParams['font.size'] = 14
    mpl.rcParams['lines.linewidth'] = 2

    filename = '%s/timing.pkl' % path
    pngtiming = '%s/timing.png' % path

    f = open(filename, 'br')
    timing = pickle.load(f)

    mean = []
    keys = []
    for k, vals in timing.items():
        mean += [np.mean(vals)]
        keys += [k]

    idx = np.argsort(mean)

    plt.figure(figsize=(10, 5))
    for k in idx[::-1]:
        vals = timing[keys[k]]
        plt.loglog(vals, label=keys[k])

    plt.xlabel('iteration')
    plt.ylabel('time [s]')
    plt.legend(loc='upper right', fontsize=10)
    plt.savefig(pngtiming) 
Example 26
Project: numpy-ml   Author: ddbourgin   File: plots.py    GNU General Public License v3.0 5 votes vote down vote up
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") 
Example 27
Project: pyFileFixity   Author: lrq3000   File: benchmark_plot.py    MIT License 5 votes vote down vote up
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) 
Example 28
Project: wlsice   Author: impaktor   File: fbm_msd_holtsberg.py    GNU General Public License v3.0 5 votes vote down vote up
def plot(t, n, msd, msd_disp):
    a = 2*H
    msd_theory = np.power(t[0:n],a)

    if isPlottingOn:
        plt.loglog(t[0:n],msd, 'b', label="fBm")
        #plt.loglog(t[0:n],mean_disp, 'r')
        plt.loglog(t[0:n],msd_theory, 'g--', label="theory")
        plt.xlabel('t')
        plt.ylabel('MSD')
        plt.savefig('fbm.png', bbox_inches='tight') 
Example 29
Project: empymod   Author: empymod   File: dc_rho-a_dip-dip.py    Apache License 2.0 5 votes vote down vote up
def plotit(depth, a, n, res1, res2, res3, title):
    """Call `calc_appres` and plot result."""

    # Calculate the three different models
    rho1, AB2 = calc_appres(depth, res1, a, n)
    rho2, _ = calc_appres(depth, res2, a, n)
    rho3, _ = calc_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        |
# +--------+---------------------+---------------------+ 
Example 30
Project: EvoDynamic   Author: SocratesNFR   File: cellular_automata_1d_criticality.py    Apache License 2.0 5 votes vote down vote up
def plot_distribution(distribution, title=""):
  import matplotlib.pyplot as plt
  plt.figure()
  x = np.linspace(1,len(distribution),len(distribution))
  plt.loglog(x, distribution, "o-", label="pdf")
  
  plt.loglog(x, [distribution[0]*xx**(-5.) for xx in x], label="slope=-5")
  plt.legend(loc='upper right')
  plt.title(title)
  plt.show() 
Example 31
Project: PyMimircache   Author: 1a1a11a   File: profilerUtils.py    GNU General Public License v3.0 5 votes vote down vote up
def draw2d(*args, **kwargs):


    figname = kwargs.get("figname", "2dPlot.png")

    if "plot_type" in kwargs:
        if kwargs['plot_type'] == "scatter":
            l = args[0]
            plt.scatter([i+1 for i in range(len(l))], l, label=kwargs.get("label", None))
    else:
        if 'logX' in kwargs and kwargs["logX"]:
            if 'logY' in kwargs and kwargs["logY"]:
                plt.loglog(*args, label=kwargs.get("label", None))
            else:
                plt.semilogx(*args, label=kwargs.get("label", None))
        else:
            if 'logY' in kwargs and kwargs["logY"]:
                plt.semilogy(*args, label=kwargs.get("label", None))
            else:
                plt.plot(*args, label=kwargs.get("label", None))

    set_fig(**kwargs)

    if not kwargs.get("no_save", False):
        # if folder does not exist, create the folder
        dname = os.path.dirname(figname)
        if dname and not os.path.exists(dname):
            os.makedirs(dname)
        plt.savefig(figname, dpi=600)
        if not kwargs.get("no_print_info", False):
            INFO("plot is saved as {}".format(figname))

    if not kwargs.get("no_show", False):
        try: plt.show()
        except: pass

    if not kwargs.get("no_clear", False):
        plt.clf() 
Example 32
Project: metaTOR   Author: koszullab   File: figures.py    GNU General Public License v3.0 5 votes vote down vote up
def draw_logplots(enrichment_files, output):

    proportional = "proportion" in enrichment_files or "proportion" in output

    enrichment_data = np.loadtxt(enrichment_files, usecols=(1, 2))
    sizes = enrichment_data[:, 0]
    hits = enrichment_data[:, 1]

    if SEABORN:
        sns.regplot(sizes, hits, fit_reg=False)
    else:
        plt.scatter(sizes, hits)

    plt.xlabel("Bin size")
    plt.ylabel("Number of hits")

    if proportional:
        plt.xscale("log")
        plt.title("Bin size vs relative hits")
    else:
        plt.loglog()
        plt.title("Bin size vs hits")
    plt.savefig(
        "{}_nhit.pdf".format(output), bbox_inches="tight", pad_inches=0.0
    )
    plt.close() 
Example 33
Project: lndmanage   Author: bitromortac   File: example_fee_market.py    MIT License 5 votes vote down vote up
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)

    plt.hist(fee_rates, bins=bins_log)
    plt.loglog()
    plt.xlabel("Fee rate bins")
    plt.ylabel("Number of channels")
    plt.show() 
Example 34
Project: lndmanage   Author: bitromortac   File: example_fee_market.py    MIT License 5 votes vote down vote up
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)
    print(bins_log)

    plt.hist(base_fees, bins=bins_log)
    plt.loglog()
    plt.xlabel("Base rate bins [msat]")
    plt.ylabel("Number of channels")
    plt.show() 
Example 35
Project: lndmanage   Author: bitromortac   File: example_fee_market.py    MIT License 5 votes vote down vote up
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)
    print(bins_log)

    plt.hist(time_locks, bins=bins_log)
    plt.loglog()
    plt.xlabel("CLTV bins [blocks]")
    plt.ylabel("Number of channels")
    plt.show() 
Example 36
Project: proptools   Author: mvernacc   File: turbopump.py    MIT License 5 votes vote down vote up
def turbine_efficiency_demo():
    plt.figure()
    uco = np.logspace(np.log10(0.04), np.log10(0.5))
    eta = [ssi_turbine_efficiency(u) for u in uco]
    plt.loglog(uco, eta)
    plt.xlabel('$u / c_o$')
    plt.ylabel('$\eta$')
    plt.grid(True) 
Example 37
Project: Emergence   Author: LennonLab   File: macroecotools.py    MIT License 5 votes vote down vote up
def plot_color_by_pt_dens(x, y, radius, loglog=0, plot_obj=None):
    """Plot bivariate relationships with large n using color for point density
    
    Inputs:
    x & y -- variables to be plotted
    radius -- the linear distance within which to count points as neighbors
    loglog -- a flag to indicate the use of a loglog plot (loglog = 1)
    
    The color of each point in the plot is determined by the logarithm (base 10)
    of the number of points that occur with a given radius of the focal point,
    with hotter colors indicating more points. The number of neighboring points
    is determined in linear space regardless of whether a loglog plot is
    presented.

    """
    plot_data = count_pts_within_radius(x, y, radius, loglog)
    sorted_plot_data = np.array(sorted(plot_data, key=lambda point: point[2]))
    
    if plot_obj == None:
        plot_obj = plt.axes()
        
    if loglog == 1:
        plot_obj.set_xscale('log')
        plot_obj.set_yscale('log')
        plot_obj.scatter(sorted_plot_data[:, 0], sorted_plot_data[:, 1],
                         c = np.sqrt(sorted_plot_data[:, 2]), edgecolors='none')
        plot_obj.set_xlim(min(x) * 0.5, max(x) * 2)
        plot_obj.set_ylim(min(y) * 0.5, max(y) * 2)
    else:
        plot_obj.scatter(sorted_plot_data[:, 0], sorted_plot_data[:, 1],
                    c = log10(sorted_plot_data[:, 2]), edgecolors='none')
    return plot_obj 
Example 38
Project: gamma_limits_sensitivity   Author: mahnen   File: __init__.py    MIT License 5 votes vote down vote up
def plot_ul_spectrum_figure(
        t_obs, lambda_lim, a_eff_interpol, e_0, n_points_to_plot, fmt='k', label=''
        ):
    '''
    fill a ul spectrum figure with the integral spectral exclusion zone plot
    '''
    gamma_range = [-5, -0.5]

    energy_range = [
        sensitive_energy(gamma_range[0], a_eff_interpol),
        sensitive_energy(gamma_range[1], a_eff_interpol)]
    energy_x = 10**np.linspace(
            np.log10(energy_range[0]),
            np.log10(energy_range[1]),
            n_points_to_plot
        )
    dn_de_y = [integral_spectral_exclusion_zone(
                energy,
                lambda_lim,
                a_eff_interpol,
                t_obs,
                e_0
                )
               for energy
               in energy_x
               ]
    dn_de_y = np.array(dn_de_y)

    plt.plot(energy_x, dn_de_y, fmt, label=label)
    plt.loglog()
    plt.title('Integral Spectral Exclusion Zone, t' +
              ('={0:1.1e} h'.format(t_obs/3600.)))
    plt.xlabel('E / TeV')
    plt.ylabel('dN/dE / [(cm$^2$ s TeV)$^{-1}$]')

    return energy_x, dn_de_y 
Example 39
Project: lenstronomy   Author: sibirrer   File: test_nfw_vir_trunc.py    MIT License 5 votes vote down vote up
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 
Example 40
Project: glosim2   Author: cosmo-epfl   File: LC.py    MIT License 5 votes vote down vote up
def plot_learning_curve(nTrains,scoreTest,err_scoreTest,Nfold,fout=None):
    with sns.plotting_context("notebook", font_scale=1.5):
        with sns.axes_style("ticks"):
            #plt.loglog(nTrains,scoreTest,'-ob',basex=10,basey=10);
            plt.errorbar(nTrains, scoreTest, yerr=err_scoreTest,fmt='o', capthick=1.5,capsize=3);
            plt.xscale('log')
            plt.yscale('log')
            plt.xlabel('Number of Training Samples')
            plt.ylabel('Test MAE [kJ/mol]')
            plt.title('Learning Curve testing on {:.0f}% of the dataset'.format(100*(1./Nfold)))
            if fout:
                plt.savefig(fout,dpi=300,bbox_inches='tight') 
Example 41
Project: pyodesys   Author: bjodah   File: _robertson.py    BSD 2-Clause "Simplified" License 5 votes vote down vote up
def plot(xout, yout, info):
    import matplotlib.pyplot as plt
    for idx in range(yout.shape[1]):
        plt.loglog(xout, yout[:, idx], label='ABC'[idx])
    plt.legend() 
Example 42
Project: alphacsc   Author: alphacsc   File: plot_output.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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) 
Example 43
Project: alphacsc   Author: alphacsc   File: greedy_coordinate_descent.py    BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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) 
Example 44
Project: pymcfost   Author: cpinte   File: SED.py    MIT License 5 votes vote down vote up
def verif(self):
        # Compares the SED from step 1 and step 2 to check if energy is properly conserved
        current_fig = plt.gcf().number
        plt.figure(20)
        plt.loglog(self._wl_th, np.sum(self._sed_th[0, :, :], axis=0), color="black")
        plt.loglog(self.wl, np.sum(self._sed_mc[0, 0, :, :], axis=0), color="red")
        plt.figure(current_fig) 
Example 45
Project: pymcfost   Author: cpinte   File: SED.py    MIT License 5 votes vote down vote up
def plot_Tr(self, h_r=0.05, log=True, symbol=None, **kwargs):

        grid = self.check_grid()

        T = self.T
        if grid.ndim > 2:
            Voronoi = False
            r_mcfost = grid[0, 0, :, :]
            z_mcfost = grid[1, 0, :, :]
        else:
            Voronoi = True
            r_mcfost = np.sqrt(grid[0, :] ** 2 + grid[1, :] ** 2)
            ou = r_mcfost > 1e-6  # Removing star
            T = T[ou]
            r_mcfost = r_mcfost[ou]
            z_mcfost = grid[2, ou]

        # Selecting data points
        ou = abs(z_mcfost) / r_mcfost < h_r

        # If we have cylindrical grid, we search for closest grid point
        # Todo

        r_mcfost = r_mcfost[ou]
        T = T[ou]

        if log:
            plt.loglog(r_mcfost, T, "o", linewidth=0.2, **kwargs)
        else:
            plt.plot(r_mcfost, T, "o", linewidth=0.2, **kwargs)
        plt.xlabel("z [au]")
        plt.ylabel("T [K]") 
Example 46
Project: extrapolation   Author: numerical-mathematics   File: tests.py    MIT License 5 votes vote down vote up
def tst_convergence(f, t0, tf, y0, order, exact, method, title="tst_convergence"):
    '''
        Runs a convergence test, integrating the system of initial value problems
        y'(t) = f(y, t) using a sequence of fixed step sizes with the provided
        extrapolation method.
        Creates a plot of the resulting errors versus step size with a reference
        line with the given order to compare with the method error 

        **Inputs**:
            - f         -- the right hand side function of the IVP.
                        Must output a non-scalar numpy.ndarray
            - [t0, tf]  -- the interval of integration
            - y0        -- the value of y(t0). Must be a non-scalar numpy.ndarray
            - order     -- the order of extrapolation
            - exact     -- the exact solution to the IVP.
                        Must output a non-scalar numpy.ndarray
            - method    -- the extrapolation method function
            - title     -- the title of the graph produced (optional)
    '''
    hs = np.asarray([2**(-k) for k in range(10)])
    err = np.zeros(len(hs))

    for i in range(len(hs)):
        y, _ = method(f, t0, tf, y0, p=order, step_size=(hs[i]), adaptive="fixed")
        err[i] = np.linalg.norm(y - exact(tf))

    plt.hold('true')
    method_err,  = plt.loglog(hs, err, 's-')
    order_line,  = plt.loglog(hs, (hs**order)*(err[5]/hs[5]))
    plt.legend([method_err, order_line], ['Method Error', 'order ' + str(order)], loc=2)
    plt.title(title)
    plt.ylabel('||error||')
    plt.xlabel('time step size')
    plt.show()
    return err 
Example 47
Project: extrapolation   Author: numerical-mathematics   File: tests.py    MIT License 5 votes vote down vote up
def tst_adaptive(f, t0, tf, y0, order, exact, method, title="tst_adaptive"):
    '''
        Runs a test, integrating a system of initial value problems y'(t) = f(y, t)
        with the given adaptive step size and adaptive order extrapolation method
        using a sequence of tolerances for local error.
        Creates a plot of the number of f evaluations versus the global error.

        **Inputs**:
            - f         -- the right hand side function of the IVP.
                        Must output a non-scalar numpy.ndarray
            - [t0, tf]  -- the interval of integration
            - y0        -- the value of y(t0). Must be a non-scalar numpy.ndarray
            - order     -- the order of extrapolation. 
            - exact     -- the exact solution to the IVP. 
                        Must output a non-scalar numpy.ndarray
            - method    -- the extrapolation method function
            - title     -- the title of the graph produced (optional)
    '''
    tol = np.asarray([2**(-k) for k in range(16)])
    err_step = np.zeros(len(tol))
    fe_step = np.zeros(len(tol))
    err_order = np.zeros(len(tol))
    fe_order = np.zeros(len(tol))

    for i in range(len(tol)):
        y, fe_step[i] = method(f, t0, tf, y0, p=order, atol=(tol[i]), rtol=(tol[i]), exact=exact, adaptive="step")
        err_step[i] = np.linalg.norm(y - exact(tf))
        y, fe_order[i], _, _ = method(f, t0, tf, y0, p=order, atol=(tol[i]), rtol=(tol[i]), exact=exact, adaptive="order")
        err_order[i] = np.linalg.norm(y - exact(tf))

    plt.hold('true')
    line_step, =plt.loglog(err_step, fe_step, 's-')
    line_order, =plt.loglog(err_order, fe_order, 's-')
    plt.legend([line_step, line_order], ["adaptive step", "adaptive order"], loc=2)

    plt.title(title + ' [order ' + str(order) + ']')
    plt.xlabel('||error||')
    plt.ylabel('fe')
    plt.show()
    return (err_step, err_order) 
Example 48
Project: extrapolation   Author: numerical-mathematics   File: tests.py    MIT License 5 votes vote down vote up
def test_dense():
    f = f_6
    exact = exact_6
    tol = [1.e-3,1.e-5,1.e-7,1.e-9,1.e-11,1.e-13]
    t0 = 0
    # t0 = random.random()*math.pi
    t_max = t0 + 2*math.pi
    y0 = exact(t0)
    t = np.linspace(t0,t_max,100)
    ys_exact = np.array([exact(t_) for t_ in t])
    err = np.zeros(len(tol))
    
    for i in range(len(tol)):
        print("tol = ", tol[i])
        ys = ex_p.ex_midpoint_parallel(f, y0, t, atol=tol[i], rtol=tol[i], adaptive="order")
        plt.hold(True)
        line_exact, = plt.plot(t, ys_exact, "s-")
        line_sol, = plt.plot(t, ys, "s-")
        plt.legend([line_exact, line_sol], ["exact", "sol w/ tol =" + str(tol[i])], loc=4)
        plt.show()
        err[i] = relative_error(ys, ys_exact)
        line_error, = plt.semilogy(t, abs(ys_exact - ys), "s-")
        plt.legend([line_error], ["error"], loc=4)
        plt.show()

    plt.loglog(tol, err, "s-")
    plt.show() 
Example 49
Project: cusumtools   Author: shadowk29   File: psdfit.py    GNU General Public License v3.0 5 votes vote down vote up
def main():
    root = tkinter.Tk()
    root.withdraw()
    psdfile = tkinter.filedialog.askopenfilename()
    psdfile = psdfile.replace('.bin','.psd')
    root.destroy()
    B = 100000
    N = 10000
    psd = pd.read_csv(psdfile,sep='\t',names=['f','S','integral','norm'])
    f = psd['f'].values[1:N]
    df = f[1]-f[0]
    S = psd['norm'].values[1:N]
    popt, pcov = curve_fit(fitfunc, f, np.log10(S), p0=[1.0,1,1000.0, 0.0001], sigma=np.sqrt(np.arange(1,N)+np.sqrt(3)/3), maxfev=100000)
    print(popt)
    print(np.sqrt(np.diag(pcov)))


    f0 = popt[0]
    alpha = popt[1]
    fstar = popt[2]
    offset = popt[3]
    
    pl.loglog(f,S)
    pl.loglog(f, 10**fitfunc(f, f0, alpha, fstar, offset))
    pl.show()

    print('Adjusted L: {0:.6f}'.format(corrected_L(f[:100], S[:100], f0, alpha, fstar, offset, df, B)))
    print('Old L: {0:.6f}'.format(old_L(f[:100], S[:100], df, B)))
    print('New alpha: {0:.6f}'.format(popt[1])) 
Example 50
Project: cusumtools   Author: shadowk29   File: noise-fit.py    GNU General Public License v3.0 5 votes vote down vote up
def fit_spectrum(self):
        self.p0 = [self.thermal, self.pink, self.brown]
        #pl.loglog(self.f,self.Pxx)
        #pl.show()
        popt, pcov = curve_fit(psd_fit, self.f, np.log10(self.Pxx), self.p0, sigma=np.sqrt(np.arange(1,len(self.f)+1)+np.sqrt(3)/3), maxfev=100000)
        self.thermal, self.pink, self.brown = popt
        print(self.pink)
        #pl.loglog(self.f,self.Pxx,self.f,10**psd_fit(self.f, self.thermal, self.pink, self.brown))
        #pl.show()