In [4]:

    

# %load NX01_master.py
#!/usr/bin/env python

"""
Created by stevertaylor
Copyright (c) 2014 Stephen R. Taylor

Code contributions by Rutger van Haasteren (piccard) and Justin Ellis (PAL/PAL2).

"""

import os, math, optparse, time, cProfile
from time import gmtime, strftime
from collections import OrderedDict
import h5py as h5

import numpy as np
from numpy import *

from scipy import integrate
from scipy import optimize
from scipy import constants
from numpy import random
from scipy import special as ss
from scipy import linalg as sl

import numexpr as ne
import ephem
from ephem import *

import libstempo as T2

import NX01_AnisCoefficients as anis
import NX01_utils as utils
import NX01_psr

import pyximport
pyximport.install(setup_args={"include_dirs":np.get_include()},
                  reload_support=True)

import NX01_jitter as jitter


    

In [36]:

    

%matplotlib inline
%config InlineBackend.figure_format = 'retina'
%load_ext autoreload
%autoreload 2
%load_ext line_profiler

from __future__ import division
import matplotlib
import matplotlib.pyplot as plt
matplotlib.rcParams['savefig.dpi'] = 1.5 * matplotlib.rcParams['savefig.dpi']


    

    




The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

In [64]:

    

plt.rcParams.update(plt.rcParamsDefault)
params = {'backend': 'pdf',
        'axes.labelsize': 10,
        'lines.markersize': 4,
        'font.size': 10,
        'xtick.major.size':6,
        'xtick.minor.size':3,  
        'ytick.major.size':6,
        'ytick.minor.size':3, 
        'xtick.major.width':0.5,
        'ytick.major.width':0.5,
        'xtick.minor.width':0.5,
        'ytick.minor.width':0.5,
        'lines.markeredgewidth':1,
        'axes.linewidth':1.2,
        'legend.fontsize': 7,
        'xtick.labelsize': 10,
        'ytick.labelsize': 10,
        'savefig.dpi':200,
        'path.simplify':True,
        'font.family': 'serif',
        'font.serif':'Times',
        'text.latex.preamble': [r'\usepackage{amsmath}', r'\usepackage{/Users/staylor/Research/NANOGrav/apjfonts}'],
        'text.usetex':True,
        'axes.color_cycle': ['b', 'lime', 'r', 'purple', 'g', 'c', 'm', 'orange', 'darkblue', \
                                'darkcyan', 'y','orangered','chartreuse','brown','deeppink','lightgreen', 'k'],
        #'font.serif':cm,
        #'figure.figsize': (3.39,2.1)}
        'figure.figsize': (3.39,2.5)}
plt.rcParams.update(params)


    

In [9]:

    

from_h5 = True
psrlist = './PsrListings.txt'
nmodes = 15
dmVar = False
ptmcmc = True
num_gwfreq_wins = 1
LMAX = 0
use_gpu = False
fix_slope = True
limit_or_detect_gwb = 'limit'
limit_or_detect_red = 'limit'
anis_modefile = None
fullN = True
num_psrs = 3
cadence = None


    

In [10]:

    

# Do you want to use GPU acceleration?
if use_gpu:
    import pycuda.autoinit
    from pycuda.compiler import SourceModule
    import pycuda.gpuarray as gpuarray
    import pycuda.driver as drv
    import pycuda.elementwise as el
    import pycuda.tools as tools
    import scikits.cuda.linalg as culinalg
    import scikits.cuda.misc as cumisc

    culinalg.init()


    

In [11]:

    

if nmodes:
    print "\n You've given me the number of frequencies to include in the low-rank time-frequency approximation, got it?\n"
else:
    print "\n You've given me the sampling cadence for the observations, which determines the upper frequency limit and the number of modes, got it?\n"

if ptmcmc:
    import PALInferencePTMCMC as PAL
else:
    import pymultinest


    

    




 You've given me the number of frequencies to include in the low-rank time-frequency approximation, got it?

In [12]:

    

################################################################################################################################
# PASSING THROUGH TEMPO2 VIA libstempo
################################################################################################################################

psr_pathinfo = np.genfromtxt(psrlist, dtype=str, skip_header=2) # name, hdf5-path, par-path, tim-path

if from_h5:

    tmp_psr = []
    for ii,tmp_name in enumerate(psr_pathinfo[:num_psrs,0]):
        tmp_psr.append(h5.File(psr_pathinfo[ii,1], 'r')[tmp_name])

    psr = [NX01_psr.PsrObjFromH5(p) for p in tmp_psr]
    
else:
    
    print 'Are you sure you do not want to use hdf5 files (recommended)?'
    
    t2psr=[]
    for ii in range(num_psrs):
        t2psr.append( T2.tempopulsar( parfile=psr_pathinfo[ii,2], timfile=psr_pathinfo[ii,3] ) )
        t2psr[ii].fit(iters=3)
        if np.any(np.isfinite(t2psr.residuals())==False)==True:
            t2psr = T2.tempopulsar( parfile=psr_pathinfo[ii,2], timfile=psr_pathinfo[ii,3] )

    psr = [NX01_psr.PsrObj(p) for p in t2psr]


# Grab all the pulsar quantities
[p.grab_all_vars() for p in psr]


    

    




--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 







    
Out[12]:





[None, None, None]

In [13]:

    

# Now, grab the positions and compute the ORF basis functions
psr_positions = [np.array([psr[ii].psr_locs[0],
                           np.pi/2. - psr[ii].psr_locs[1]])
                           for ii in range(len(psr))]
positions = np.array(psr_positions).copy()

CorrCoeff = np.array(anis.CorrBasis(positions,LMAX))       # Computing all the correlation basis-functions for the array.
harm_sky_vals = utils.SetupPriorSkyGrid(LMAX)              # Computing the values of the spherical-harmonics up to order
                                                                # LMAX on a pre-specified grid


    

In [14]:

    

if anis_modefile is None:
    gwfreqs_per_win = int(1.*nmodes/(1.*num_gwfreq_wins)) # getting the number of GW frequencies per window
    anis_modefreqs = np.arange(1,nmodes+1)
    anis_modefreqs = np.reshape(anis_modefreqs, (num_gwfreq_wins,gwfreqs_per_win))

    tmp_num_gwfreq_wins = num_gwfreq_wins
else:
    tmp_modefreqs = np.loadtxt(anis_modefile)
    tmp_num_gwfreq_wins = tmp_modefreqs.shape[0]
    anis_modefreqs = []
    for ii in range(tmp_num_gwfreq_wins):
        anis_modefreqs.append(np.arange(tmp_modefreqs[ii,0],tmp_modefreqs[ii,1]+1))

# Create a tag for evolving anisotropy searches
if (LMAX!=0) and (tmp_num_gwfreq_wins > 1):
    evol_anis_tag = 'EvAnis'
else:
    evol_anis_tag = ''


    

In [16]:

    

#############################################################################
# GETTING MAXIMUM TIME, COMPUTING FOURIER DESIGN MATRICES, AND GETTING MODES 
#############################################################################

Tmax = np.max([p.toas.max() - p.toas.min() for p in psr])

if nmodes:

    [p.makeTe(nmodes, Tmax, makeDM=dmVar) for p in psr]
    # get GW frequencies
    fqs = np.linspace(1/Tmax, nmodes/Tmax, nmodes)
    nmode = nmodes

else:

    nmode = int(round(0.5*Tmax/cadence))
    [p.makeTe(nmode, Tmax, makeDM=dmVar) for p in psr]
    # get GW frequencies
    fqs = np.linspace(1/Tmax, nmode/Tmax, nmode)


    

In [17]:

    

#######################################
# PRE-COMPUTING WHITE NOISE PROPERTIES 
#######################################

loglike1 = 0
TtNT = []
d = []
for ii,p in enumerate(psr):

    # compute ( T.T * N^-1 * T ) & log determinant of N
    new_err = (p.toaerrs).copy()
    if fullN==True:
        
        if len(p.ecorrs)>0:

            Jamp = np.ones(len(p.epflags))
            for jj,nano_sysname in enumerate(p.sysflagdict['nano-f'].keys()):
                Jamp[np.where(p.epflags==nano_sysname)] *= p.ecorrs[nano_sysname]**2.0

            Nx = jitter.cython_block_shermor_0D(p.res, new_err**2., Jamp, p.Uinds)
            d.append(np.dot(p.Te.T, Nx))
            logdet_N, TtNT_dummy = jitter.cython_block_shermor_2D(p.Te, new_err**2., Jamp, p.Uinds)
            TtNT.append(TtNT_dummy)
            det_dummy, dtNdt = jitter.cython_block_shermor_1D(p.res, new_err**2., Jamp, p.Uinds)

        else:
            
            d.append(np.dot(p.Te.T, p.res/( new_err**2.0 )))
        
            N = 1./( new_err**2.0 )
            right = (N*p.Te.T).T
            TtNT.append(np.dot(p.Te.T, right))
    
            logdet_N = np.sum(np.log( new_err**2.0 ))
        
            # triple product in likelihood function
            dtNdt = np.sum(p.res**2.0/( new_err**2.0 ))
        
    else:

        d.append(np.dot(p.Te.T, p.res/( new_err**2.0 )))
            
        N = 1./( new_err**2.0 )
        right = (N*p.Te.T).T
        TtNT.append(np.dot(p.Te.T, right))

        logdet_N = np.sum(np.log( new_err**2.0 ))
        
        # triple product in likelihood function
        dtNdt = np.sum(p.res**2.0/( new_err**2.0 ))

    loglike1 += -0.5 * (logdet_N + dtNdt)

d = np.concatenate(d)


    

In [18]:

    

##########################
# SETTING UP PRIOR RANGES
##########################

pmin = -20.0*np.ones(len(psr))
pmin = np.append(pmin,0.0*np.ones(len(psr)))
if dmVar==True:
    pmin = np.append(pmin,-20.0*np.ones(len(psr)))
    pmin = np.append(pmin,0.0*np.ones(len(psr)))
pmin = np.append(pmin,-18.0)
if fix_slope==False:
    pmin = np.append(pmin,0.0)
pmin = np.append(pmin,-10.0*np.ones( tmp_num_gwfreq_wins*(((LMAX+1)**2)-1) ))


pmax = -11.0*np.ones(len(psr))
pmax = np.append(pmax,7.0*np.ones(len(psr)))
if dmVar==True:
    pmax = np.append(pmax,-11.0*np.ones(len(psr)))
    pmax = np.append(pmax,7.0*np.ones(len(psr)))
pmax = np.append(pmax,-11.0)
if fix_slope==False:
    pmax = np.append(pmax,7.0)
pmax = np.append(pmax,10.0*np.ones( tmp_num_gwfreq_wins*(((LMAX+1)**2)-1) ))

##################################################################################


    

In [19]:

    

def my_prior(xx):
    logp = 0.
    
    if np.all(xx <= pmax) and np.all(xx >= pmin):
        logp = np.sum(np.log(1/(pmax-pmin)))
    else:
        logp = -np.inf
    
    return logp


    

In [59]:

    

def lnprob(xx):

    npsr = len(psr)
    #print npsr, nmodes

    if dmVar==True:
        Ared, gam_red, Adm, gam_dm, Agwb, gam_gwb, orf_coeffs = utils.masterSplitParams(xx, npsr, dmVar, fix_slope)
        mode_count = 4*nmode
    else:
        Ared, gam_red, Agwb, gam_gwb, orf_coeffs = utils.masterSplitParams(xx, npsr, dmVar, fix_slope)
        mode_count = 2*nmode


    # Reshaping freq-dependent anis coefficients,
    # and testing for power distribution physicality.

    orf_coeffs = orf_coeffs.reshape((tmp_num_gwfreq_wins,((LMAX+1)**2)-1))
    clm = np.array([[0.0]*((LMAX+1)**2) for ii in range(tmp_num_gwfreq_wins)])
    clm[:,0] = 2.0*np.sqrt(np.pi)

    physicality = 0.
    if LMAX!=0:

        for kk in range(tmp_num_gwfreq_wins):
            for ii in range(1,((LMAX+1)**2)):
                clm[kk,ii] = orf_coeffs[kk,ii-1]   

            # Testing for physicality of power distribution.
            if (utils.PhysPrior(clm[kk],harm_sky_vals) == 'Unphysical'):
                physicality += -10.0**7.0
            else:
                physicality += 0.

    
    # Computing frequency dependent overlap reduction functions.
    ORF=[]
    for ii in range(tmp_num_gwfreq_wins): # number of frequency windows
        for jj in range(len(anis_modefreqs[ii])): # number of frequencies in this window
            ORF.append( sum(clm[ii,kk]*CorrCoeff[kk] for kk in range(len(CorrCoeff))) )
    if dmVar==True:
        for ii in range(tmp_num_gwfreq_wins): # number of frequency windows
            for jj in range(len(anis_modefreqs[ii])): # number of frequencies in this window
                ORF.append( np.zeros((npsr,npsr)) )

    ORF = np.array(ORF)
    ORFtot = np.zeros((mode_count,npsr,npsr)) # shouldn't be applying ORF to dmfreqs,
                                                 # but the projection of GW spec onto dmfreqs
                                                 # is defined as zero below.
    ORFtot[0::2] = ORF
    ORFtot[1::2] = ORF
    
       
    # parameterize intrinsic red noise as power law
    Tspan = (1/fqs[0])*86400.0
    f1yr = 1/3.16e7
    rho = np.log10(Agwb**2/12/np.pi**2 * f1yr**(gam_gwb-3) * (fqs/86400.0)**(-gam_gwb)/Tspan)

    # parameterize intrinsic red-noise and DM-variations as power law
    kappa = [] 
    if dmVar==True:
        for ii in range(npsr):
            kappa.append(np.log10( np.append( Ared[ii]**2/12/np.pi**2 * f1yr**(gam_red[ii]-3) * (fqs/86400.0)**(-gam_red[ii])/Tspan,
                                        Adm[ii]**2/12/np.pi**2 * f1yr**(gam_dm[ii]-3) * (fqs/86400.0)**(-gam_dm[ii])/Tspan ) ))
    else:
        for ii in range(npsr):
            kappa.append(np.log10( Ared[ii]**2/12/np.pi**2 * f1yr**(gam_red[ii]-3) * (fqs/86400.0)**(-gam_red[ii])/Tspan ))
    

    # construct elements of sigma array
    sigdiag = []
    sigoffdiag = []
    
    if dmVar==True:
        gwbspec = np.append( 10**rho, np.zeros_like(rho) )
    else:
        gwbspec = 10**rho
        
    for ii in range(npsr):
        tot = np.zeros(mode_count)
        offdiag = np.zeros(mode_count)

        # off diagonal terms
        offdiag[0::2] = gwbspec
        offdiag[1::2] = gwbspec

        # diagonal terms
        tot[0::2] = ORF[:,ii,ii]*gwbspec + 10**kappa[ii]
        tot[1::2] = ORF[:,ii,ii]*gwbspec + 10**kappa[ii] 
                
        # fill in lists of arrays
        sigdiag.append(tot)
        sigoffdiag.append(offdiag)


    # compute Phi matrix
    smallMatrix = np.zeros((mode_count, npsr, npsr))
    for ii in range(npsr):
        for jj in range(ii,npsr):

            if ii == jj:
                smallMatrix[:,ii,jj] = sigdiag[jj] 
            else:
                smallMatrix[:,ii,jj] = ORFtot[:,ii,jj] * sigoffdiag[jj] 
                smallMatrix[:,jj,ii] = smallMatrix[:,ii,jj]

    
    # invert Phi matrix frequency-wise
    logdet_Phi = 0
    non_pos_def = 0
    for ii in range(mode_count):

        try:

            L = sl.cho_factor(smallMatrix[ii,:,:])
            smallMatrix[ii,:,:] = sl.cho_solve(L, np.eye(npsr))
            logdet_Phi += np.sum(2*np.log(np.diag(L[0])))

        except np.linalg.LinAlgError:

            print 'Cholesky Decomposition Failed!! Rejecting...'
            non_pos_def += 1

    # Break if we have non-positive-definiteness of Phi
    if non_pos_def > 0:

        return -np.inf

    else:

        bigTtNT = sl.block_diag(*TtNT)
        Phi = np.zeros_like( bigTtNT )

        # now fill in real covariance matrix
        ind = [0]
        ind = np.append(ind,np.cumsum([TtNT[ii].shape[0] for ii in range(len(psr))]))
        ind = [np.arange(ind[ii]+psr[ii].Gc.shape[1],ind[ii]+psr[ii].Gc.shape[1]+mode_count)
               for ii in range(len(ind)-1)]
        for ii in range(npsr):
            for jj in range(npsr):
                Phi[ind[ii],ind[jj]] = smallMatrix[:,ii,jj]
            
        # compute sigma
        Sigma = bigTtNT + Phi
            
        # cholesky decomp for second term in exponential
        if use_gpu:

            try:

                Sigma_gpu = gpuarray.to_gpu( Sigma.astype(np.float64).copy() )
                expval2_gpu = gpuarray.to_gpu( d.astype(np.float64).copy() )
                culinalg.cho_solve( Sigma_gpu, expval2_gpu ) # in-place linear-algebra:
                                                             # Sigma and expval2 overwritten
                logdet_Sigma = np.sum(2.0*np.log(np.diag(Sigma_gpu.get())))

            except cula.culaDataError:

                print 'Cholesky Decomposition Failed (GPU error!!)'
                return -np.inf

            logLike = -0.5 * (logdet_Phi + logdet_Sigma) + 0.5 * (np.dot(d, expval2_gpu.get() )) + loglike1
            
        else:

            try:

                cf = sl.cho_factor(Sigma)
                expval2 = sl.cho_solve(cf, d)
                logdet_Sigma = np.sum(2*np.log(np.diag(cf[0])))

            except np.linalg.LinAlgError:

                print 'Cholesky Decomposition Failed second time!! Using SVD instead'
                u,s,v = sl.svd(Sigma)
                expval2 = np.dot(v.T, 1/s*np.dot(u.T, d))
                logdet_Sigma = np.sum(np.log(s))


            logLike = -0.5 * (logdet_Phi + logdet_Sigma) + 0.5 * (np.dot(d, expval2)) + loglike1 


        # Multiplying likelihood to correct log-uniform
        # sampling thus making a uniform prior
        if limit_or_detect_gwb == 'limit':
            priorfac_gwb = np.log(Agwb * np.log(10.0))
        else:
            priorfac_gwb = 0.0

        if limit_or_detect_red == 'limit':
            priorfac_red = np.sum(np.log(Ared * np.log(10.0)))
        else:
            priorfac_red = 0.0


        return logLike + priorfac_gwb + priorfac_red + physicality


    

In [ ]:

    

nmodes_arr = [15,30,50] #[10,15,20,25,30,35,40,45,50]
num_psrs_arr = np.arange(1,19)

timings = np.zeros((len(nmodes_arr),len(num_psrs_arr)))


    

In [79]:

    

nmodes = 50

for mm in range(len(num_psrs_arr)):

    ################################################################################################################################
    # PASSING THROUGH TEMPO2 VIA libstempo
    ################################################################################################################################

    psr_pathinfo = np.genfromtxt(psrlist, dtype=str, skip_header=2) # name, hdf5-path, par-path, tim-path

    if from_h5:

        tmp_psr = []
        for ii,tmp_name in enumerate(psr_pathinfo[:num_psrs_arr[mm],0]):
            tmp_psr.append(h5.File(psr_pathinfo[ii,1], 'r')[tmp_name])

        psr = [NX01_psr.PsrObjFromH5(p) for p in tmp_psr]

    else:

        print 'Are you sure you do not want to use hdf5 files (recommended)?'

        t2psr=[]
        for ii in range(num_psrs_arr[mm]):
            t2psr.append( T2.tempopulsar( parfile=psr_pathinfo[ii,2], timfile=psr_pathinfo[ii,3] ) )
            t2psr[ii].fit(iters=3)
            if np.any(np.isfinite(t2psr.residuals())==False)==True:
                t2psr = T2.tempopulsar( parfile=psr_pathinfo[ii,2], timfile=psr_pathinfo[ii,3] )

        psr = [NX01_psr.PsrObj(p) for p in t2psr]


    # Grab all the pulsar quantities
    [p.grab_all_vars() for p in psr]


    # Now, grab the positions and compute the ORF basis functions
    psr_positions = [np.array([psr[ii].psr_locs[0],
                               np.pi/2. - psr[ii].psr_locs[1]])
                               for ii in range(len(psr))]
    positions = np.array(psr_positions).copy()

    CorrCoeff = np.array(anis.CorrBasis(positions,LMAX))       # Computing all the correlation basis-functions for the array.
    harm_sky_vals = utils.SetupPriorSkyGrid(LMAX)              # Computing the values of the spherical-harmonics up to order
                                                               # LMAX on a pre-specified grid


    if anis_modefile is None:
        gwfreqs_per_win = int(1.*nmodes/(1.*num_gwfreq_wins)) # getting the number of GW frequencies per window
        anis_modefreqs = np.arange(1,nmodes+1)
        anis_modefreqs = np.reshape(anis_modefreqs, (num_gwfreq_wins,gwfreqs_per_win))

        tmp_num_gwfreq_wins = num_gwfreq_wins
    else:
        tmp_modefreqs = np.loadtxt(anis_modefile)
        tmp_num_gwfreq_wins = tmp_modefreqs.shape[0]
        anis_modefreqs = []
        for ii in range(tmp_num_gwfreq_wins):
            anis_modefreqs.append(np.arange(tmp_modefreqs[ii,0],tmp_modefreqs[ii,1]+1))

    # Create a tag for evolving anisotropy searches
    if (LMAX!=0) and (tmp_num_gwfreq_wins > 1):
        evol_anis_tag = 'EvAnis'
    else:
        evol_anis_tag = ''



    #############################################################################
    # GETTING MAXIMUM TIME, COMPUTING FOURIER DESIGN MATRICES, AND GETTING MODES 
    #############################################################################

    Tmax = np.max([p.toas.max() - p.toas.min() for p in psr])

    if nmodes:

        [p.makeTe(nmodes, Tmax, makeDM=dmVar) for p in psr]
        # get GW frequencies
        fqs = np.linspace(1/Tmax, nmodes/Tmax, nmodes)
        nmode = nmodes

    else:

        nmode = int(round(0.5*Tmax/cadence))
        [p.makeTe(nmode, Tmax, makeDM=dmVar) for p in psr]
        # get GW frequencies
        fqs = np.linspace(1/Tmax, nmode/Tmax, nmode)


    #######################################
    # PRE-COMPUTING WHITE NOISE PROPERTIES 
    #######################################

    loglike1 = 0
    TtNT = []
    d = []
    for ii,p in enumerate(psr):

        # compute ( T.T * N^-1 * T ) & log determinant of N
        new_err = (p.toaerrs).copy()
        if fullN==True:

            if len(p.ecorrs)>0:

                Jamp = np.ones(len(p.epflags))
                for jj,nano_sysname in enumerate(p.sysflagdict['nano-f'].keys()):
                    Jamp[np.where(p.epflags==nano_sysname)] *= p.ecorrs[nano_sysname]**2.0

                Nx = jitter.cython_block_shermor_0D(p.res, new_err**2., Jamp, p.Uinds)
                d.append(np.dot(p.Te.T, Nx))
                logdet_N, TtNT_dummy = jitter.cython_block_shermor_2D(p.Te, new_err**2., Jamp, p.Uinds)
                TtNT.append(TtNT_dummy)
                det_dummy, dtNdt = jitter.cython_block_shermor_1D(p.res, new_err**2., Jamp, p.Uinds)

            else:

                d.append(np.dot(p.Te.T, p.res/( new_err**2.0 )))

                N = 1./( new_err**2.0 )
                right = (N*p.Te.T).T
                TtNT.append(np.dot(p.Te.T, right))

                logdet_N = np.sum(np.log( new_err**2.0 ))

                # triple product in likelihood function
                dtNdt = np.sum(p.res**2.0/( new_err**2.0 ))

        else:

            d.append(np.dot(p.Te.T, p.res/( new_err**2.0 )))

            N = 1./( new_err**2.0 )
            right = (N*p.Te.T).T
            TtNT.append(np.dot(p.Te.T, right))

            logdet_N = np.sum(np.log( new_err**2.0 ))

            # triple product in likelihood function
            dtNdt = np.sum(p.res**2.0/( new_err**2.0 ))

        loglike1 += -0.5 * (logdet_N + dtNdt)

    d = np.concatenate(d)


    ##########################
    # SETTING UP PRIOR RANGES
    ##########################

    pmin = -20.0*np.ones(len(psr))
    pmin = np.append(pmin,0.0*np.ones(len(psr)))
    if dmVar==True:
        pmin = np.append(pmin,-20.0*np.ones(len(psr)))
        pmin = np.append(pmin,0.0*np.ones(len(psr)))
    pmin = np.append(pmin,-18.0)
    if fix_slope==False:
        pmin = np.append(pmin,0.0)
    pmin = np.append(pmin,-10.0*np.ones( tmp_num_gwfreq_wins*(((LMAX+1)**2)-1) ))


    pmax = -11.0*np.ones(len(psr))
    pmax = np.append(pmax,7.0*np.ones(len(psr)))
    if dmVar==True:
        pmax = np.append(pmax,-11.0*np.ones(len(psr)))
        pmax = np.append(pmax,7.0*np.ones(len(psr)))
    pmax = np.append(pmax,-11.0)
    if fix_slope==False:
        pmax = np.append(pmax,7.0)
    pmax = np.append(pmax,10.0*np.ones( tmp_num_gwfreq_wins*(((LMAX+1)**2)-1) ))

    ##################################################################################


    #########################
    #########################

    # Set up the parameter list

    parameters=[]
    for ii in range(len(psr)):
        parameters.append('Ared_'+psr[ii].name)
    for ii in range(len(psr)):
        parameters.append('gam_red_'+psr[ii].name)
    if dmVar==True:
        for ii in range(len(psr)):
            parameters.append('Adm_'+psr[ii].name)
        for ii in range(len(psr)):
            parameters.append('gam_dm_'+psr[ii].name)
    parameters.append("Agwb")
    if fix_slope is False:
        parameters.append("gam_gwb")
        gamma_ext = 'GamVary'
    else:
        gamma_ext = 'Gam4p33'
    for ii in range( tmp_num_gwfreq_wins*(((LMAX+1)**2)-1) ):
        parameters.append('clm_{0}'.format(ii+1))


    print "\n You are searching for the following parameters: {0}\n".format(parameters)
    n_params = len(parameters)

    print "\n The total number of parameters is {0}\n".format(n_params)

    # Start the sampling off with some reasonable parameter choices
    x0 = np.log10(np.array([p.Redamp for p in psr]))
    x0 = np.append(x0,np.array([p.Redind for p in psr]))
    if dmVar==True:
        x0 = np.append(x0,np.log10(np.array([p.Redamp for p in psr])))
        x0 = np.append(x0,np.array([p.Redind for p in psr]))
    x0 = np.append(x0,-15.0)
    if fix_slope is False:
        x0 = np.append(x0,13./3.)
    x0 = np.append(x0,np.zeros( tmp_num_gwfreq_wins*(((LMAX+1)**2)-1) ))

    print "\n Your initial parameters are {0}\n".format(x0)

    # Make a reasonable covariance matrix to commence sampling
    cov_diag = 0.5*np.ones(len(psr))
    cov_diag = np.append(cov_diag,0.5*np.ones(len(psr)))
    if dmVar==True:
        cov_diag = np.append(cov_diag,0.5*np.ones(len(psr)))
        cov_diag = np.append(cov_diag,0.5*np.ones(len(psr)))
    cov_diag = np.append(cov_diag,0.5)
    if fix_slope is False:
        cov_diag = np.append(cov_diag,0.5)
    cov_diag = np.append(cov_diag,0.05*np.ones( tmp_num_gwfreq_wins*(((LMAX+1)**2)-1) ))



    result = %timeit -o lnprob(x0)
    timings[2,mm] = result.best


    

    




--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'gam_red_J1713+0747', 'Agwb']


 The total number of parameters is 3


 Your initial parameters are [-19.1114   1.7124 -15.    ]

100 loops, best of 3: 6.19 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'Agwb']


 The total number of parameters is 5


 Your initial parameters are [-19.1114  -15.1073    1.7124    2.88933 -15.     ]

100 loops, best of 3: 8.44 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'Agwb']


 The total number of parameters is 7


 Your initial parameters are [-19.1114  -15.1073  -17.0642    1.7124    2.88933   3.21169 -15.     ]

100 loops, best of 3: 11.8 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'Agwb']


 The total number of parameters is 9


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777    1.7124    2.88933   3.21169
   3.13292 -15.     ]

100 loops, best of 3: 16.1 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'Agwb']


 The total number of parameters is 11


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129     1.7124    2.88933
   3.21169   3.13292   3.25799 -15.     ]

10 loops, best of 3: 27.1 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'Agwb']


 The total number of parameters is 13


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325    1.7124
   2.88933   3.21169   3.13292   3.25799   3.21288 -15.     ]

10 loops, best of 3: 46.2 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1744-1134 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'Ared_J1744-1134', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'gam_red_J1744-1134', 'Agwb']


 The total number of parameters is 15


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325  -15.6471
   1.7124    2.88933   3.21169   3.13292   3.25799   3.21288   3.30492 -15.     ]

10 loops, best of 3: 54.3 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1744-1134 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0030+0451 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'Ared_J1744-1134', 'Ared_J0030+0451', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'gam_red_J1744-1134', 'gam_red_J0030+0451', 'Agwb']


 The total number of parameters is 17


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325  -15.6471
 -14.412     1.7124    2.88933   3.21169   3.13292   3.25799   3.21288
   3.30492   5.05806 -15.     ]

1 loops, best of 3: 90.8 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1744-1134 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0030+0451 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0613-0200 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'Ared_J1744-1134', 'Ared_J0030+0451', 'Ared_J0613-0200', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'gam_red_J1744-1134', 'gam_red_J0030+0451', 'gam_red_J0613-0200', 'Agwb']


 The total number of parameters is 19


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325  -15.6471
 -14.412   -13.6168    1.7124    2.88933   3.21169   3.13292   3.25799
   3.21288   3.30492   5.05806   3.16059 -15.     ]

1 loops, best of 3: 104 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1744-1134 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0030+0451 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0613-0200 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1614-2230 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'Ared_J1744-1134', 'Ared_J0030+0451', 'Ared_J0613-0200', 'Ared_J1614-2230', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'gam_red_J1744-1134', 'gam_red_J0030+0451', 'gam_red_J0613-0200', 'gam_red_J1614-2230', 'Agwb']


 The total number of parameters is 21


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325  -15.6471
 -14.412   -13.6168  -16.9988    1.7124    2.88933   3.21169   3.13292
   3.25799   3.21288   3.30492   5.05806   3.16059   3.31971 -15.     ]

10 loops, best of 3: 106 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1744-1134 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0030+0451 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0613-0200 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1614-2230 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting B1855+09 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'Ared_J1744-1134', 'Ared_J0030+0451', 'Ared_J0613-0200', 'Ared_J1614-2230', 'Ared_B1855+09', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'gam_red_J1744-1134', 'gam_red_J0030+0451', 'gam_red_J0613-0200', 'gam_red_J1614-2230', 'gam_red_B1855+09', 'Agwb']


 The total number of parameters is 23


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325  -15.6471
 -14.412   -13.6168  -16.9988  -13.8485    1.7124    2.88933   3.21169
   3.13292   3.25799   3.21288   3.30492   5.05806   3.16059   3.31971
   3.88472 -15.     ]

1 loops, best of 3: 127 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1744-1134 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0030+0451 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0613-0200 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1614-2230 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting B1855+09 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1853+1303 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'Ared_J1744-1134', 'Ared_J0030+0451', 'Ared_J0613-0200', 'Ared_J1614-2230', 'Ared_B1855+09', 'Ared_J1853+1303', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'gam_red_J1744-1134', 'gam_red_J0030+0451', 'gam_red_J0613-0200', 'gam_red_J1614-2230', 'gam_red_B1855+09', 'gam_red_J1853+1303', 'Agwb']


 The total number of parameters is 25


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325  -15.6471
 -14.412   -13.6168  -16.9988  -13.8485  -16.8149    1.7124    2.88933
   3.21169   3.13292   3.25799   3.21288   3.30492   5.05806   3.16059
   3.31971   3.88472   3.39717 -15.     ]

1 loops, best of 3: 141 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1744-1134 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0030+0451 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0613-0200 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1614-2230 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting B1855+09 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1853+1303 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2145-0750 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'Ared_J1744-1134', 'Ared_J0030+0451', 'Ared_J0613-0200', 'Ared_J1614-2230', 'Ared_B1855+09', 'Ared_J1853+1303', 'Ared_J2145-0750', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'gam_red_J1744-1134', 'gam_red_J0030+0451', 'gam_red_J0613-0200', 'gam_red_J1614-2230', 'gam_red_B1855+09', 'gam_red_J1853+1303', 'gam_red_J2145-0750', 'Agwb']


 The total number of parameters is 27


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325  -15.6471
 -14.412   -13.6168  -16.9988  -13.8485  -16.8149  -15.7129    1.7124
   2.88933   3.21169   3.13292   3.25799   3.21288   3.30492   5.05806
   3.16059   3.31971   3.88472   3.39717   2.67557 -15.     ]

1 loops, best of 3: 178 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1744-1134 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0030+0451 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0613-0200 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1614-2230 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting B1855+09 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1853+1303 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2145-0750 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1455-3330 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'Ared_J1744-1134', 'Ared_J0030+0451', 'Ared_J0613-0200', 'Ared_J1614-2230', 'Ared_B1855+09', 'Ared_J1853+1303', 'Ared_J2145-0750', 'Ared_J1455-3330', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'gam_red_J1744-1134', 'gam_red_J0030+0451', 'gam_red_J0613-0200', 'gam_red_J1614-2230', 'gam_red_B1855+09', 'gam_red_J1853+1303', 'gam_red_J2145-0750', 'gam_red_J1455-3330', 'Agwb']


 The total number of parameters is 29


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325  -15.6471
 -14.412   -13.6168  -16.9988  -13.8485  -16.8149  -15.7129  -16.752
   1.7124    2.88933   3.21169   3.13292   3.25799   3.21288   3.30492
   5.05806   3.16059   3.31971   3.88472   3.39717   2.67557   3.28721 -15.     ]

1 loops, best of 3: 205 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1744-1134 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0030+0451 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0613-0200 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1614-2230 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting B1855+09 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1853+1303 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2145-0750 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1455-3330 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1012+5307 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'Ared_J1744-1134', 'Ared_J0030+0451', 'Ared_J0613-0200', 'Ared_J1614-2230', 'Ared_B1855+09', 'Ared_J1853+1303', 'Ared_J2145-0750', 'Ared_J1455-3330', 'Ared_J1012+5307', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'gam_red_J1744-1134', 'gam_red_J0030+0451', 'gam_red_J0613-0200', 'gam_red_J1614-2230', 'gam_red_B1855+09', 'gam_red_J1853+1303', 'gam_red_J2145-0750', 'gam_red_J1455-3330', 'gam_red_J1012+5307', 'Agwb']


 The total number of parameters is 31


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325  -15.6471
 -14.412   -13.6168  -16.9988  -13.8485  -16.8149  -15.7129  -16.752
 -12.6189    1.7124    2.88933   3.21169   3.13292   3.25799   3.21288
   3.30492   5.05806   3.16059   3.31971   3.88472   3.39717   2.67557
   3.28721   1.02018 -15.     ]

1 loops, best of 3: 262 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1744-1134 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0030+0451 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0613-0200 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1614-2230 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting B1855+09 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1853+1303 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2145-0750 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1455-3330 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1012+5307 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1741+1351 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'Ared_J1744-1134', 'Ared_J0030+0451', 'Ared_J0613-0200', 'Ared_J1614-2230', 'Ared_B1855+09', 'Ared_J1853+1303', 'Ared_J2145-0750', 'Ared_J1455-3330', 'Ared_J1012+5307', 'Ared_J1741+1351', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'gam_red_J1744-1134', 'gam_red_J0030+0451', 'gam_red_J0613-0200', 'gam_red_J1614-2230', 'gam_red_B1855+09', 'gam_red_J1853+1303', 'gam_red_J2145-0750', 'gam_red_J1455-3330', 'gam_red_J1012+5307', 'gam_red_J1741+1351', 'Agwb']


 The total number of parameters is 33


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325  -15.6471
 -14.412   -13.6168  -16.9988  -13.8485  -16.8149  -15.7129  -16.752
 -12.6189  -16.9121    1.7124    2.88933   3.21169   3.13292   3.25799
   3.21288   3.30492   5.05806   3.16059   3.31971   3.88472   3.39717
   2.67557   3.28721   1.02018   3.38453 -15.     ]

1 loops, best of 3: 292 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1744-1134 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0030+0451 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0613-0200 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1614-2230 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting B1855+09 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1853+1303 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2145-0750 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1455-3330 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1012+5307 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1741+1351 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2010-1323 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'Ared_J1744-1134', 'Ared_J0030+0451', 'Ared_J0613-0200', 'Ared_J1614-2230', 'Ared_B1855+09', 'Ared_J1853+1303', 'Ared_J2145-0750', 'Ared_J1455-3330', 'Ared_J1012+5307', 'Ared_J1741+1351', 'Ared_J2010-1323', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'gam_red_J1744-1134', 'gam_red_J0030+0451', 'gam_red_J0613-0200', 'gam_red_J1614-2230', 'gam_red_B1855+09', 'gam_red_J1853+1303', 'gam_red_J2145-0750', 'gam_red_J1455-3330', 'gam_red_J1012+5307', 'gam_red_J1741+1351', 'gam_red_J2010-1323', 'Agwb']


 The total number of parameters is 35


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325  -15.6471
 -14.412   -13.6168  -16.9988  -13.8485  -16.8149  -15.7129  -16.752
 -12.6189  -16.9121  -16.4803    1.7124    2.88933   3.21169   3.13292
   3.25799   3.21288   3.30492   5.05806   3.16059   3.31971   3.88472
   3.39717   2.67557   3.28721   1.02018   3.38453   3.42768 -15.     ]

1 loops, best of 3: 346 ms per loop
--> Extracting J1713+0747 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1909-3744 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1640+2224 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1600-3053 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2317+1439 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1918-0642 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1744-1134 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0030+0451 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J0613-0200 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1614-2230 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting B1855+09 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1853+1303 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2145-0750 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1455-3330 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1012+5307 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1741+1351 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J2010-1323 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 

--> Extracting J1024-0719 from hdf5 file
--> Done extracting pulsar from hdf5 file :-) 


 You are searching for the following parameters: ['Ared_J1713+0747', 'Ared_J1909-3744', 'Ared_J1640+2224', 'Ared_J1600-3053', 'Ared_J2317+1439', 'Ared_J1918-0642', 'Ared_J1744-1134', 'Ared_J0030+0451', 'Ared_J0613-0200', 'Ared_J1614-2230', 'Ared_B1855+09', 'Ared_J1853+1303', 'Ared_J2145-0750', 'Ared_J1455-3330', 'Ared_J1012+5307', 'Ared_J1741+1351', 'Ared_J2010-1323', 'Ared_J1024-0719', 'gam_red_J1713+0747', 'gam_red_J1909-3744', 'gam_red_J1640+2224', 'gam_red_J1600-3053', 'gam_red_J2317+1439', 'gam_red_J1918-0642', 'gam_red_J1744-1134', 'gam_red_J0030+0451', 'gam_red_J0613-0200', 'gam_red_J1614-2230', 'gam_red_B1855+09', 'gam_red_J1853+1303', 'gam_red_J2145-0750', 'gam_red_J1455-3330', 'gam_red_J1012+5307', 'gam_red_J1741+1351', 'gam_red_J2010-1323', 'gam_red_J1024-0719', 'Agwb']


 The total number of parameters is 37


 Your initial parameters are [-19.1114  -15.1073  -17.0642  -16.9777  -17.129   -16.8325  -15.6471
 -14.412   -13.6168  -16.9988  -13.8485  -16.8149  -15.7129  -16.752
 -12.6189  -16.9121  -16.4803  -16.2659    1.7124    2.88933   3.21169
   3.13292   3.25799   3.21288   3.30492   5.05806   3.16059   3.31971
   3.88472   3.39717   2.67557   3.28721   1.02018   3.38453   3.42768
   3.338   -15.     ]

1 loops, best of 3: 514 ms per loop

In [86]:

    

plt.plot(np.arange(1,len(timings[0,:18])+1), timings[0,:18], lw=1.5, ls='solid', alpha=0.5, label=r'$15$ frequencies')
plt.plot(np.arange(1,len(timings[1,:18])+1), timings[1,:18], lw=1.5, ls='dashed', alpha=0.5, label=r'$30$ frequencies')
plt.plot(np.arange(1,len(timings[2,:18])+1), timings[2,:18], lw=1.5, ls='dotted', alpha=0.5, label=r'$50$ frequencies')

psr_pathinfo = np.genfromtxt(psrlist, dtype=str, skip_header=2) # name, hdf5-path, par-path, tim-path
plt.xticks(range(len(psr_pathinfo[:,0])), psr_pathinfo[:,0], rotation=45, size=8)
plt.yscale('log')
plt.ylim(1e-3,0.6)

plt.legend(loc='upper left')

plt.ylabel(r'Likelihood evaluation time [s]')

plt.savefig('nx01_timings.png',bbox_inches='tight')


    

In [37]:

    

%lprun -f lnprob lnprob(x0)