Unable to sample 2-planet GP model with MCMC

[hposborn]

I adapted this GP tutorial to work with 2 planets. However, the radvel.mcmc step fails when computing the GR statistic, and I'm not entirely sure why.

Here's some code to reproduce the error:

import radvel
import numpy as np

#generating fake data
t = np.sort(np.random.random(130)*100)
vel = np.random.normal(0.0,2.5,130)
vel+=1.25*np.sin(t*(np.pi*2)/2.49)
vel+=2.15*np.sin(t*(np.pi*2)/6.8)
vel+=4.25*np.sin(t*(np.pi*2)/21.3)
errvel = np.tile(1.5,130)

#initialising GP parameters:
hnames = [
'gp_amp', # eta_1; GP variability amplitude
'gp_explength', # eta_2; GP non-periodic characteristic length
'gp_per', # eta_3; GP variability period
'gp_perlength', # eta_4; GP periodic characteristic length
]

gp_explength_mean = 21*np.sqrt(2.) # sqrt(2)*tau in Dai+ 2017 [days]
gp_explength_unc = 1.0*np.sqrt(2.)
gp_perlength_mean = np.sqrt(1./(2.*3.32)) # sqrt(1/(2*gamma)) in Dai+ 2017
gp_perlength_unc = 0.019
gp_per_mean = 21 # T_bar in Dai+ 2017 [days]
gp_per_unc = 0.5

#planet parameters for use later:
Porb = np.array([2.5,6.98]) # orbital period [days]
Porb_unc = np.array([0.01,0.4])
Tc = np.array([2.4,17.3]) # [BJD]
Tc_unc = np.array([0.01,2.5])

#defining radvel.Parameters object
params2 = radvel.Parameters(2,basis='per tc secosw sesinw k')

params2['per1'] = radvel.Parameter(value=Porb[0])
params2['tc1'] = radvel.Parameter(value=Tc[0])
params2['sesinw1'] = radvel.Parameter(value=0.,vary=False) # fix eccentricity = 0
params2['secosw1'] = radvel.Parameter(value=0.,vary=False)
params2['k1'] = radvel.Parameter(value=4)
params2['per2'] = radvel.Parameter(value=Porb[1])
params2['tc2'] = radvel.Parameter(value=Tc[1])
params2['sesinw2'] = radvel.Parameter(value=0.,vary=False) # fix eccentricity = 0
params2['secosw2'] = radvel.Parameter(value=0.,vary=False)
params2['k2'] = radvel.Parameter(value=4)
params2['dvdt'] = radvel.Parameter(value=0.,vary=False)
params2['curv'] = radvel.Parameter(value=0.,vary=False)
params2['gp_amp'] = radvel.Parameter(value=25.0)
params2['gp_explength'] = radvel.Parameter(value=gp_explength_mean)
params2['gp_per'] = radvel.Parameter(value=gp_per_mean)
params2['gp_perlength'] = radvel.Parameter(value=gp_perlength_mean)

gpmodel2 = radvel.model.RVModel(params2)

like2 = radvel.likelihood.GPLikelihood(gpmodel2, t, vel,
                                        errvel, hnames, suffix='_harps',
                                        kernel_name="QuasiPer"
                                       )

params2['gamma_harps'] = radvel.Parameter(value=np.average(vel))
params2['jit_harps'] = radvel.Parameter(value=0.5)

gplike2 = radvel.likelihood.CompositeLikelihood([like2])
gppost2 = radvel.posterior.Posterior(gplike2)

#priors:
gppost2.priors += [radvel.prior.Gaussian('per1', Porb[0], Porb_unc[0])]
gppost2.priors += [radvel.prior.Gaussian('tc1', Tc[0], Tc_unc[0])]
gppost2.priors += [radvel.prior.Jeffreys('k1', 0.01, 10.)]
gppost2.priors += [radvel.prior.Gaussian('per2', Porb[1], Porb_unc[1])]
gppost2.priors += [radvel.prior.Gaussian('tc2', Tc[1], Tc_unc[1])]
gppost2.priors += [radvel.prior.Jeffreys('k2', 0.01, 10.)]

gppost2.priors += [radvel.prior.Jeffreys('gp_amp', 0.01, 100.)]
gppost2.priors += [radvel.prior.Jeffreys('jit_harps', 0.01,10.)]
gppost2.priors += [radvel.prior.Gaussian('gp_explength', gp_explength_mean, gp_explength_unc)]
gppost2.priors += [radvel.prior.Gaussian('gp_per', gp_per_mean, gp_per_unc)]
gppost2.priors += [radvel.prior.Gaussian('gp_perlength', gp_perlength_mean, gp_perlength_unc)]

#optimizing
from scipy import optimize

res2 = optimize.minimize(gppost2.neglogprob_array, gppost2.get_vary_params(), method='Powell')
print(gppost2)

#running the MCMC:
chains2 = radvel.mcmc(gppost2,nrun=100, savename='P2_chains.h5')

And the output I get:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-457-a6726dc8278f> in <module>
----> 1 chains2 = radvel.mcmc(gppost2,nrun=100, savename='P2_chains.h5')

~/transits/lib/python3.7/site-packages/radvel/mcmc.py in mcmc(post, nwalkers, nrun, ensembles, checkinterval, minAfactor, maxArchange, burnAfactor, burnGR, maxGR, minTz, minsteps, minpercent, thin, serial, save, savename, proceed, proceedname)
    380 
    381             convergence_check(minAfactor=minAfactor, maxArchange=maxArchange, maxGR=maxGR, minTz=minTz,
--> 382                               minsteps=minsteps, minpercent=minpercent)
    383 
    384             if save:

~/transits/lib/python3.7/site-packages/radvel/mcmc.py in convergence_check(minAfactor, maxArchange, maxGR, minTz, minsteps, minpercent)
    148             = convergence_calculate(statevars.chains,
    149                                     oldautocorrelation=statevars.oac, minAfactor=minAfactor, maxArchange=maxArchange,
--> 150                                     maxGR=maxGR, minTz=minTz)
    151         statevars.mintz = min(tz)
    152         statevars.maxgr = max(gr)

~/transits/lib/python3.7/site-packages/radvel/mcmc.py in convergence_calculate(chains, oldautocorrelation, minAfactor, maxArchange, minTz, maxGR)
    581     afactor = np.divide(chains.shape[0], autocorrelation)
    582 
--> 583     archange = np.divide(np.abs(np.subtract(autocorrelation, oldautocorrelation)), oldautocorrelation)
    584 
    585     # well-mixed criteria

ValueError: operands could not be broadcast together with shapes (12,) (9,) 

[hposborn]

Ah. This actually only happens when I have previously run some sort of radvel model with 1 planet (hence the 9 parameters) - the code I wrote to reproduce the error doesn't work if it's run "fresh". So clearly radvel.mcmc is remembering the 1-planet model, despite the fact it's initiated with a whole new model... Is there any way to "wipe" radvel's memory, or do I need to restart my Jupyter notebook kernel every time?

[bjfultn]

It's not easy to adjust the number of planets in a Posterior object after it is initialized with a given set of parameters. Try re-initializing the RVmodel, GPLikelihood, CompositeLikelihood, and Posterior objects.

[hposborn]

Thanks for the quick response! I'm not sure how to reinitialise those classes deeper than I have tried here. I guess my initial model isn't there, but params2 is a separate and disconnected dictionary from the initial params, and each of the created objects has a different name (gppost2, gplike2, etc). Yet, when fed into radvel.mcmc, it keeps some sort of memory! So is there some variable to feed to RVmodel, GPLikelihood, CompositeLikelihood, and Posterior (and maybe radvel.mcmc) which reinitialises completely?

Here's an example. If you run this snippet in a Jupyter notebook after running the above example, it fails:

params3 = radvel.Parameters(1,basis='per tc secosw sesinw k')

params3['per1'] = radvel.Parameter(value=6.9)
params3['tc1'] = radvel.Parameter(value=17.3)
params3['sesinw1'] = radvel.Parameter(value=0.,vary=False) # fix eccentricity = 0
params3['secosw1'] = radvel.Parameter(value=0.,vary=False)
params3['k1'] = radvel.Parameter(value=4)
params3['dvdt'] = radvel.Parameter(value=0.,vary=False)
params3['curv'] = radvel.Parameter(value=0.,vary=False)
params3['gp_amp'] = radvel.Parameter(value=25.0)
params3['gp_explength'] = radvel.Parameter(value=gp_explength_mean)
params3['gp_per'] = radvel.Parameter(value=gp_per_mean)
params3['gp_perlength'] = radvel.Parameter(value=gp_perlength_mean)

gpmodel3 = radvel.model.RVModel(params3)

like3 = radvel.likelihood.GPLikelihood(gpmodel3, t, vel,
                                        errvel, hnames, suffix='_harps',
                                        kernel_name="QuasiPer"
                                       )

params3['gamma_harps'] = radvel.Parameter(value=np.average(vel))
params3['jit_harps'] = radvel.Parameter(value=0.5)

gplike3 = radvel.likelihood.CompositeLikelihood([like3])
gppost3 = radvel.posterior.Posterior(gplike3)

gppost3.priors += [radvel.prior.Gaussian('per1', 6.9, 0.5)]
gppost3.priors += [radvel.prior.Gaussian('tc1', 17.3, 2.5)]
gppost3.priors += [radvel.prior.Jeffreys('k1', 0.01, 10.)]

gppost3.priors += [radvel.prior.Jeffreys('gp_amp', 0.01, 100.)]
gppost3.priors += [radvel.prior.Jeffreys('jit_harps', 0.01,10.)]
gppost3.priors += [radvel.prior.Gaussian('gp_explength', gp_explength_mean, gp_explength_unc)]
gppost3.priors += [radvel.prior.Gaussian('gp_per', gp_per_mean, gp_per_unc)]
gppost3.priors += [radvel.prior.Gaussian('gp_perlength', gp_perlength_mean, gp_perlength_unc)]

from scipy import optimize

res3 = optimize.minimize(gppost3.neglogprob_array, gppost3.get_vary_params(), method='Powell')
print(gppost3)

chains3 = radvel.mcmc(gppost3,nrun=100, savename='P2_chains.h5')

[bjfultn]

OK, I found the bug. The radvel.statevars namespace/class which is used to coordinate the MCMC run across multiple threads was not being reset between MCMC runs. For now you can add a radvel.statevars.reset() call between your mcmc runs.

I'll add this reset call into the radvel.mcmc function so that it is automatically called at the start of any new MCMC run. This change will be released with v1.3.8.

Please verify that the radvel.statevars.reset() trick works for you and close the issue if it does.

[hposborn]

Yep, that works. Thanks for the fix (and for a very useful bit of code in general)!