Sampling Galactic Binary parameters: a demo

import numpy as np
import matplotlib.pyplot as plt
import pickle
import m3c2.proposal as proposal
from m3c2.sampler import PTSampler
import m3c2.gb as gbtools
import corner
import scipy

import multiprocessing as mp
mp.set_start_method('fork') #needed for MacOS
import warnings
warnings.filterwarnings('ignore')

Galactic binary fast waveform is given by the LDC toolbox, as well as other common tools used here and in the mc3.gb submodule.

from ldc.waveform.fastGB import FastGB
from ldc.common.series import TDI
from ldc.lisa.noise import get_noise_model

Generate some data with a single source

src = {'Amplitude': 4.532e-23, 'EclipticLatitude': -0.61,'EclipticLongitude': 4.36,
       'Frequency': 0.0100680913,'FrequencyDerivative': 1.79e-15,
       'Inclination': 0.67,'InitialPhase': 5.48,'Polarization': 0.43}
fr_min = 0.010063
fr_max = 0.010073
dt = 15
Tobs = 365*24*60*60
GB = FastGB(delta_t=dt, T=Tobs)

X, Y, Z = GB.get_td_tdixyz(template=src)
AET = TDI(dict(zip(["X", "Y", "Z"], [X,Y,Z])))
AET.XYZ2AET()

plt.figure(figsize=(15,5))
plt.subplot(121);AET["A"].plot()
plt.subplot(122);AET["E"].plot()

[<matplotlib.lines.Line2D at 0x7fabfa656a90>]
png

Add some noise

def fftnoise(f):
    f = np.array(f, dtype='complex')
    Np = (len(f) - 1) // 2
    phases = np.random.rand(Np) * 2 * np.pi
    phases = np.cos(phases) + 1j * np.sin(phases)
    f[1:Np+1] *= phases
    f[-1:-1-Np:-1] = np.conj(f[1:Np+1])
    return np.fft.ifft(f).real

df = 1/Tobs
frange = np.linspace(0, 1/dt, int(1/(df*dt)))
noise = get_noise_model("MRDv1", frange)
for k in ["A", "E", "T"]:
    psd = noise.psd(option=k)
    psd[0] = 0
    AET[k] += fftnoise(np.sqrt(len(psd)*psd/dt/2))

plt.figure(figsize=(15,5))
plt.subplot(121);AET["A"].plot()
plt.subplot(122);AET["E"].plot()

[<matplotlib.lines.Line2D at 0x7fabfa542cd0>]
png 
f, psdX =  scipy.signal.welch(AET["A"].values, fs=1.0/dt, window='hamming', nperseg=256*256)
psd = noise.psd(option='A'); psd[0] = 0

plt.figure(figsize=(15,5))
plt.loglog(f, psdX, label="simulated data")
plt.loglog(frange, psd, label="noise model")
plt.ylabel("PSD A"); plt.xlabel("Freq [Hz]")
plt.legend()
plt.axis([1e-6, 1e-1, 1e-44, 1e-37])

(1e-06, 0.1, 1e-44, 1e-37)
png

Use the flat log-likelihood function

data = TDI(dict(zip(["A", "E", "T"],[AET["A"].ts.fft(),AET["E"].ts.fft(), AET["T"].ts.fft()])))
nsource = 1
fll = gbtools.FlatLogLik(data, GB, noise, [fr_min, fr_max], GB.T, nsource)

# Evaluate at true value
pars = fll.as_search_param(src)
l0 = fll.loglik(pars, plotIt=True)
print(f'Likelihood at true value is {l0}')

Likelihood at true value is -157.84451496017988
png

Define the sampler

# Parallel tempering parameters
Tmax = 10
Nchains = 5

# Set a starting point
priors = fll.make_prior()
#x0 = [gbtools.random_point(priors) for n in range(Nchains)]
x0 = [fll.as_search_param(src) for n in range(Nchains)] # start from truth for the demo

S = PTSampler(Nchains, fll.loglik, fll.logprior, fll.names, Tmax=Tmax, profiling=True)
S.set_starting_point(x0)

INFO:root:Tmax=10, range of temperature: [ 1.          1.77827941  3.16227766  5.62341325 10.        ]

Define the proposals

# Define proposals
SL = proposal.Slice(fll.names).slice
SC = proposal.SCAM(fll.names).SCAM
HS = gbtools.HarmonicShift(fll.names).shift
p_dict = [{SL:40, SC:70, HS:10}]*Nchains
S.set_proposals(p_dict)

Run the sampler

# Run mcmc
niter = 5000
c = S.run_mcmc(niter, pSwap=0.95, printN=1000, multiproc=True, n0_swap=2500)

INFO:root:iter 0 chain 0
INFO:root:current loglik: -157.8, best: -157.8, temp: 1.0, ratio: 1.0
INFO:root:iter 0 chain 1
INFO:root:current loglik: -157.8, best: -157.8, temp: 1.8, ratio: 1.0
INFO:root:iter 0 chain 2
INFO:root:current loglik: -157.8, best: -157.8, temp: 3.2, ratio: 1.0
INFO:root:iter 0 chain 3
INFO:root:current loglik: -157.8, best: -157.8, temp: 5.6, ratio: 1.0
INFO:root:iter 0 chain 4
INFO:root:current loglik: -157.8, best: -157.8, temp: 10.0, ratio: 1.0
INFO:root:iter 1000 chain 1
INFO:root:current loglik: -162.9, best: -157.4, temp: 1.8, ratio: 1.0
INFO:root:iter 1000 chain 3
INFO:root:current loglik: -184.6, best: -157.7, temp: 5.6, ratio: 1.0
INFO:root:iter 1000 chain 0
INFO:root:current loglik: -163.8, best: -157.3, temp: 1.0, ratio: 1.0
INFO:root:iter 1000 chain 2
INFO:root:current loglik: -161.8, best: -157.5, temp: 3.2, ratio: 1.0
INFO:root:iter 1000 chain 4
INFO:root:current loglik: -220.5, best: -157.8, temp: 10.0, ratio: 1.0
INFO:root:iter 2000 chain 1
INFO:root:current loglik: -157.8, best: -157.3, temp: 1.8, ratio: 1.0
INFO:root:iter 2000 chain 3
INFO:root:current loglik: -177.4, best: -157.7, temp: 5.6, ratio: 1.0
INFO:root:iter 2000 chain 0
INFO:root:current loglik: -159.5, best: -157.2, temp: 1.0, ratio: 1.0
INFO:root:iter 2000 chain 2
INFO:root:current loglik: -166.7, best: -157.5, temp: 3.2, ratio: 1.0
INFO:root:iter 2000 chain 4
INFO:root:current loglik: -199.6, best: -157.8, temp: 10.0, ratio: 1.0
INFO:root:iter 3000 chain 0
INFO:root:current loglik: -162.4, best: -157.2, temp: 1.0, ratio: 0.316008316008316
INFO:root:iter 3000 chain 1
INFO:root:current loglik: -168.5, best: -157.3, temp: 1.8, ratio: 0.3127659574468085
INFO:root:iter 3000 chain 2
INFO:root:current loglik: -192.4, best: -157.5, temp: 3.1, ratio: 0.3857442348008386
INFO:root:iter 3000 chain 3
INFO:root:current loglik: -195.8, best: -157.7, temp: 5.5, ratio: 0.6804123711340206
INFO:root:iter 3000 chain 4
INFO:root:current loglik: -201.1, best: -157.8, temp: 10.0, ratio: 1.0
INFO:root:iter 4000 chain 0
INFO:root:current loglik: -159.8, best: -157.0, temp: 1.0, ratio: 0.26586741889985893
INFO:root:iter 4000 chain 1
INFO:root:current loglik: -165.0, best: -157.3, temp: 1.8, ratio: 0.32723434201266716
INFO:root:iter 4000 chain 2
INFO:root:current loglik: -185.0, best: -157.5, temp: 3.1, ratio: 0.4188811188811189
INFO:root:iter 4000 chain 3
INFO:root:current loglik: -188.2, best: -157.7, temp: 5.3, ratio: 0.7016806722689075
INFO:root:iter 4000 chain 4
INFO:root:current loglik: -188.2, best: -157.8, temp: 10.0, ratio: 1.0
INFO:root:Temperature ladder adjustment processus is done
INFO:root:iter 4999 chain 0
INFO:root:current loglik: -158.4, best: -157.0, temp: 1.0, ratio: 0.30046550994498517
INFO:root:iter 4999 chain 1
INFO:root:current loglik: -164.0, best: -157.3, temp: 1.8, ratio: 0.36153198653198654
INFO:root:iter 4999 chain 2
INFO:root:current loglik: -179.7, best: -157.5, temp: 3.1, ratio: 0.46697517879680267
INFO:root:iter 4999 chain 3
INFO:root:current loglik: -209.1, best: -157.7, temp: 5.1, ratio: 0.7002114164904862
INFO:root:iter 4999 chain 4
INFO:root:current loglik: -209.1, best: -157.8, temp: 10.0, ratio: 1.0
INFO:root:Done !

Show likelihood evolution

L = np.load("chain_0.npy")['logL']
plt.figure()
plt.plot(L)
plt.axis([None, None, -1000, 0])

(-249.95000000000002, 5248.95, -1000.0, 0.0)
png

Show posterior

chn = np.load("chain_0.npy")
#chn_ = chn[:,:-1]
chn_ = np.array(chn[fll.names].tolist())
truth = fll.as_search_param(src)
labels = fll.names
fig = corner.corner(chn_,  bins=50, hist_kwargs={'density':True, 'lw':3},
                    plot_datapoints=False, fill_contours=False, quantiles=[0.05, 0.5, 0.95],
                    show_titles=True,
                    color='darkcyan', truths=truth, truth_color='k', use_math_test=True,
                    levels=[0.9], title_kwargs={"fontsize": 12}, labels=labels)
png

Log-likelihood at best value

chn = np.load("chain_0.npy")
loglik = chn['logL']
ibest = np.argmax(loglik)
chn = np.array(chn[fll.names].tolist())
pbest = chn[ibest,:]
l1 = fll.loglik(pbest, plotIt=True)
png