Python Pgplot.plotxy方法代码示例

# 需要导入模块: import Pgplot [as 别名]
# 或者: from Pgplot import plotxy [as 别名]
 def plot_chi2_vs_sub(self, device='/xwin'):
     """
     plot_chi2_vs_sub(self, device='/xwin'):
         Plot (and return) an array showing the reduced-chi^2 versus
             the subband number.
     """
     # Sum the profiles in each subband
     profs = self.profs.sum(0)
     # Compute the averages and variances for the subbands
     avgs = profs.sum(1)/self.proflen
     vars = []
     for sub in range(self.nsub):
         var = 0.0
         if sub in self.killed_subbands:
             vars.append(var)
             continue
         for part in range(self.npart):
             if part in self.killed_intervals:
                 continue
             var += self.stats[part][sub][5] # foldstats prof_var
         vars.append(var)
     chis = Num.zeros(self.nsub, dtype='f')
     for ii in range(self.nsub):
         chis[ii] = self.calc_redchi2(prof=profs[ii], avg=avgs[ii], var=vars[ii])
     # Now plot it
     Pgplot.plotxy(chis, labx="Subband Number", laby="Reduced-\gx\u2\d",
                   rangey=[0.0, max(chis)*1.1], device=device)
     return chis

# 需要导入模块: import Pgplot [as 别名]
# 或者: from Pgplot import plotxy [as 别名]
def kuiper_uniform_test(data, output=0):
    """
    kuiper_uniform_test(data, output=0):
       Conduct a Kuiper test on the data.  The data must be values
       within [0,1) (e.g. phases from a periodicity search).  They
       will be compared to a uniform distribution.  The return value
       is the probability that the data is uniformly distributed.
    """
    sdata = num.asarray(data)
    N = sdata.size
    sdata.sort()
    f0 = num.arange(N, dtype=num.float64)/N
    fn = (num.arange(N, dtype=num.float64)+1.0)/N
    Dp = (fn - sdata).max()
    Dm = (sdata - f0).max()
    D = Dp + Dm
    P = kuiper_prob(D, N)
    if (output):
        xs = (num.arange(N+3, dtype=num.float64)/(N+2.0)).repeat(2)[1:-1]
        ys = num.concatenate((num.asarray([0.0]), sdata, num.asarray([1.0]))).repeat(2)
        Pgplot.plotxy(ys, xs, rangex=[-0.03, 1.03], rangey=[-0.03, 1.03], aspect=1.0, 
                      labx="Fraction of Data", laby="Cumulative Value", width=2)
        Pgplot.plotxy(num.asarray([0.0, 1.0]), num.asarray([0.0, 1.0]), width=1)
        Pgplot.closeplot()
        print("Max distance between the cumulative distributions (D) = %.5g" % D)
        print("Prob the data is from the specified distrbution   (P) = %.3g" % P)
    return (D, P)

# 需要导入模块: import Pgplot [as 别名]
# 或者: from Pgplot import plotxy [as 别名]
def estimate_rz(psr, T, show=0, device='/XWIN'):
    """
    estimate_rz(psr, T, show=0, device='/XWIN'):
        Return estimates of a pulsar's average Fourier freq ('r')
        relative to its nominal Fourier freq as well as its
        Fourier f-dot ('z') in bins, of a pulsar.
           'psr' is a psrparams structure describing the pulsar.
           'T' is the length of the observation in sec.
           'show' if true, displays plots of 'r' and 'z'.
           'device' if the device to plot to if 'show' is true.
    """
    startE = keplers_eqn(psr.orb.t, psr.orb.p, psr.orb.e, 1.0E-15)
    numorbpts = int(T / psr.orb.p + 1.0) * 1024 + 1
    dt = T / (numorbpts - 1)
    E = dorbint(startE, numorbpts, dt, psr.orb)
    z = z_from_e(E, psr, T)
    r = T/p_from_e(E, psr) - T/psr.p
    if show:
        times = np.arange(numorbpts) * dt
        Pgplot.plotxy(r, times, labx = 'Time', \
                      laby = 'Fourier Frequency (r)', device=device)
        if device=='/XWIN':
            print 'Press enter to continue:'
            i = raw_input()
        Pgplot.nextplotpage()
        Pgplot.plotxy(z, times, labx = 'Time',
                      laby = 'Fourier Frequency Derivative (z)', device=device)
        Pgplot.closeplot()
    return r.mean(), z.mean()

# 需要导入模块: import Pgplot [as 别名]
# 或者: from Pgplot import plotxy [as 别名]
 def plot_sumprof(self, device='/xwin'):
     """
     plot_sumprof(self, device='/xwin'):
         Plot the dedispersed and summed profile.
     """
     if not self.__dict__.has_key('subdelays'):
         print "Dedispersing first..."
         self.dedisperse()
     normprof = self.sumprof - min(self.sumprof)
     normprof /= max(normprof)
     Pgplot.plotxy(normprof, labx="Phase Bins", laby="Normalized Flux",
                   device=device)

# 需要导入模块: import Pgplot [as 别名]
# 或者: from Pgplot import plotxy [as 别名]
 def plot_chi2_vs_DM(self, loDM, hiDM, N=100, interp=0, device='/xwin'):
     """
     plot_chi2_vs_DM(self, loDM, hiDM, N=100, interp=0, device='/xwin'):
         Plot (and return) an array showing the reduced-chi^2 versus
             DM (N DMs spanning loDM-hiDM).  Use sinc_interpolation
             if 'interp' is non-zero.
     """
     # Sum the profiles in time
     sumprofs = self.profs.sum(0)
     if not interp:
         profs = sumprofs
     else:
         profs = Num.zeros(Num.shape(sumprofs), dtype='d')
     DMs = psr_utils.span(loDM, hiDM, N)
     chis = Num.zeros(N, dtype='f')
     subdelays_bins = self.subdelays_bins.copy()
     for ii, DM in enumerate(DMs):
         subdelays = psr_utils.delay_from_DM(DM, self.barysubfreqs)
         hifreqdelay = subdelays[-1]
         subdelays = subdelays - hifreqdelay
         delaybins = subdelays*self.binspersec - subdelays_bins
         if interp:
             interp_factor = 16
             for jj in range(self.nsub):
                 profs[jj] = psr_utils.interp_rotate(sumprofs[jj], delaybins[jj],
                                                     zoomfact=interp_factor)
             # Note: Since the interpolation process slightly changes the values of the
             # profs, we need to re-calculate the average profile value
             avgprof = (profs/self.proflen).sum()
         else:
             new_subdelays_bins = Num.floor(delaybins+0.5)
             for jj in range(self.nsub):
                 profs[jj] = psr_utils.rotate(profs[jj], int(new_subdelays_bins[jj]))
             subdelays_bins += new_subdelays_bins
             avgprof = self.avgprof
         sumprof = profs.sum(0)
         chis[ii] = self.calc_redchi2(prof=sumprof, avg=avgprof)
     # Now plot it
     Pgplot.plotxy(chis, DMs, labx="DM", laby="Reduced-\gx\u2\d", device=device)
     return (chis, DMs)

# 需要导入模块: import Pgplot [as 别名]
# 或者: from Pgplot import plotxy [as 别名]
     if tmpnumbins > numbins:  numbins = tmpnumbins
 # Powers averaged over orb.t as a function of orb.w
 pwrs_w = zeros((orbsperpt[ctype], numbins), Float32)
 for ct in range(orbsperpt[ctype]):
     wb = ct * 180.0 / orbsperpt[ctype]
     if debugout:  print('wb = '+repr(wb))
     psr = psrparams_from_list([pp, Pb, xb, ecc[ctype], wb, 0.0])
     for i in range(numffts):
         psr.orb.t = i * Tfft
         tmppwrs = spectralpower(gen_bin_response(0.0, numbetween,
                                                  psr.p, Tfft,
                                                  psr.orb, numbins))
         if debugout:  print('     tb = '+repr(psr.orb.t)+'  Max pow = '+\
            repr(max(tmppwrs)))
         if showplots:
             Pgplot.plotxy(tmppwrs)
             Pgplot.closeplot()
         pwrs_w[ct] = pwrs_w[ct] + tmppwrs
     if showsumplots:
         Pgplot.plotxy(pwrs_w[ct], title='power(w) averaged over orb.t')
         Pgplot.closeplot()
 pwrs_w = pwrs_w / numffts
 max_avg_pow = average(maximum.reduce(pwrs_w,1))
 if showsumplots:
     Pgplot.plotxy(add.reduce(pwrs_w), title='power(w) averaged over orb.t')
     Pgplot.closeplot()
 tim = clock() - stim
 if debugout:
     print('Time for this point was ',tim, ' s.')
 file.write('%8.6f  %10.5f  %10d  %13.9f\n' % \
            (pp, Tfft, int(Tfft/dt), max_avg_pow))

# 需要导入模块: import Pgplot [as 别名]
# 或者: from Pgplot import plotxy [as 别名]
     wb, tp = 0.0, ct * Pb / orbsperpt[ctype]
 else:
     (orbf, orbi)  = modf(ct / sqrt(orbsperpt[ctype]))
     orbi = orbi / sqrt(orbsperpt[ctype])
     wb, tp = orbf * 180.0, Pb * orbi
 if debugout:
     print 'T = '+`T`+'  ppsr = '+`ppsr[y]`+\
           ' Pb = '+`Pb`+' xb = '+`xb`+' eb = '+\
           `eb`+' wb = '+`wb`+' tp = '+`tp`
 psr = psrparams_from_list([ppsr[y], Pb, xb, eb, wb, tp])
 psr_numbins = 2 * bin_resp_halfwidth(psr.p, T, psr.orb)
 psr_resp = gen_bin_response(0.0, 1, psr.p, T, psr.orb,
                             psr_numbins)
 if showplots:
     print "The raw response:"
     Pgplot.plotxy(spectralpower(psr_resp))
     Pgplot.closeplot()
 # The following places the nominative psr freq
 # approx in bin len(data)/2
 datalen = next2_to_n(psr_numbins * 2)
 if datalen < 1024: datalen = 1024
 data = zeros(datalen, 'F')
 lo = (len(data) - len(psr_resp)) / 2
 hi = lo + len(psr_resp)
 data[lo:hi] = array(psr_resp, copy=1)
 (tryr, tryz) = estimate_rz(psr, T, show=showplots)
 tryr = tryr + len(data) / 2.0
 numr = 200
 numz = 200
 dr = 0.5
 dz = 1.0

# 需要导入模块: import Pgplot [as 别名]
# 或者: from Pgplot import plotxy [as 别名]
         template = sinc_interp.periodic_interp(template, numbins)[::oldlen]
 else:
     if gaussfitfile is not None:
         template = psr_utils.read_gaussfitfile(gaussfitfile, numbins)
     else:
         template = psr_utils.gaussian_profile(numbins, 0.0, gaussianwidth)
 # Normalize it
 template -= min(template)
 template /= max(template)
 # Rotate it so that it becomes a "true" template according to FFTFIT
 shift,eshift,snr,esnr,b,errb,ngood = measure_phase(template, template)
 template = psr_utils.fft_rotate(template, shift)
     
 # Determine the off-pulse bins
 if bkgd_vals is not None:
     Pgplot.plotxy(template, labx="Phase bins")
     Pgplot.plotxy(template[bkgd_vals], Num.arange(numbins)[bkgd_vals],
                   line=None, symbol=2, color='red')
     Pgplot.closeplot()
     offpulse_inds = bkgd_vals
     onpulse_inds = set(Num.arange(numbins)) - set(bkgd_vals)
 else:
     offpulse_inds = Num.compress(template<=bkgd_cutoff, Num.arange(numbins))
     onpulse_inds = Num.compress(template>bkgd_cutoff, Num.arange(numbins))
     Pgplot.plotxy(template)
     Pgplot.plotxy([bkgd_cutoff, bkgd_cutoff], [0.0, numbins], color='red')
     Pgplot.closeplot()
 # If the number of bins in the offpulse section is < 10% of the total
 # use the statistics in the .pfd file to set the RMS
 if (len(offpulse_inds) < 0.1*numbins):
     print "Number of off-pulse bins to use for RMS is too low.  Using .pfd stats."

# 需要导入模块: import Pgplot [as 别名]
# 或者: from Pgplot import plotxy [as 别名]
                    orbi = orbi / sqrt(orbsperpt[ctype])
                    wb, tp = orbf * 180.0, Pb * orbi

                # Generate the PSR response
                psr = psrparams_from_list([ppsr[y], Pb, xb, ecc[ctype], wb, tp])
                psr_numbins = 2 * bin_resp_halfwidth(psr.p, T, psr.orb)
                psr_resp = gen_bin_response(0.0, 1, psr.p, T, psr.orb,
                                            psr_numbins)
                if debugout:
                    print 'T = %9.3f  Pb = %9.3f  Ppsr = %9.7f' % \
                          (T, psr.orb.p, psr.p)

                newpows = slice_resp(psr, T, spectralpower(psr_resp))
                if showplots:
                    print "The raw response:"
                    Pgplot.plotxy(newpows)
                    Pgplot.closeplot()
                fftlen = len(newpows)
                noise = rng.sample(fftlen)
                tryamp[ct] = 500.0
                theo_sum_pow = powersum_at_sigma(detect_sigma,
                                                 int(T/psr.orb.p))
                if debugout:
                    print 'theo_sum_pow = ', theo_sum_pow
                newloop = 1
                tryamp[ct] = secant(mini_fft_sum_pows, tryamp[ct]/2,
                                    tryamp[ct], 0.01)
                # Pgplot.plotxy(spectralpower(fdata)[1:]/norm, \
                #              arange(len(fdata))*T/fftlen, \
                #              labx='Orbital Period (s))', \
                #              laby='Power')

# 需要导入模块: import Pgplot [as 别名]
# 或者: from Pgplot import plotxy [as 别名]
           presto.TWOPI*psr.orb.x/psr.p
     print ''
         
 # Create the data set
 cand = presto.orbitparams()
 m = 0
 comb = presto.gen_bin_response(0.0, 1, psr.p, T, psr.orb , 
                                presto.LOWACC, m)
 ind = len(comb)
 # The follwoing is performed automatically in gen_bin_resp() now
 # m = (ind / 2 + 10) * numbetween
 data = Numeric.zeros(3 * ind, 'F')
 data[ind:2*ind] = comb
 if showplots and not parallel:
     Pgplot.plotxy(presto.spectralpower(data), color='red',
                   title='Data', labx='Fourier Frequency',
                   laby='Relative Power')
     a = raw_input("Press enter to continue...")
     Pgplot.nextplotpage(1)
     
 # Perform the loops over the Keplerian parameters
 for job in range(numjobs):
     if parallel:
         myjob = work[myid]
     else:
         myjob = work[job]
     if myjob=='p':
         Dd = Dp
         psrref = psr.orb.p
     if myjob=='x':
         Dd = Dx

# 需要导入模块: import Pgplot [as 别名]
# 或者: from Pgplot import plotxy [as 别名]
     wb, tp = 0.0, ct * Pb / orbsperpt[ctype]
 else:
     (orbf, orbi)  = modf(ct / sqrt(orbsperpt[ctype]))
     orbi = orbi / sqrt(orbsperpt[ctype])
     wb, tp = orbf * 180.0, Pb * orbi
 if debugout:
     print('T = '+repr(T)+'  ppsr = '+repr(ppsr[y])+\
           ' Pb = '+repr(Pb)+' xb = '+repr(xb)+' eb = '+\
           repr(eb)+' wb = '+repr(wb)+' tp = '+repr(tp))
 psr = psrparams_from_list([ppsr[y], Pb, xb, eb, wb, tp])
 psr_numbins = 2 * bin_resp_halfwidth(psr.p, T, psr.orb)
 psr_resp = gen_bin_response(0.0, 1, psr.p, T, psr.orb,
                             psr_numbins)
 if showplots:
     print("The raw response:")
     Pgplot.plotxy(spectralpower(psr_resp))
     Pgplot.closeplot()
 if searchtype == 'ffdot':
     # The following places the nominative psr freq
     # approx in bin len(data)/2
     datalen = next2_to_n(psr_numbins * 2)
     if datalen < 1024: datalen = 1024
     data = zeros(datalen, 'F')
     lo = (len(data) - len(psr_resp)) / 2
     hi = lo + len(psr_resp)
     data[lo:hi] = array(psr_resp, copy=1)
     (tryr, tryz) = estimate_rz(psr, T, show=showplots)
     tryr = tryr + len(data) / 2.0
     numr = 200
     numz = 200
     dr = 0.5