How to fix the period?

[AstroJLin]

Hi,

I want to fix the period, which should how to set up?

Thanks!

[adrn]

Hi @AstroJLin - here's a small example: https://gist.github.com/56707a370423d883a33bcc9049bdbf7a

[AstroJLin]

Hi,
Here, I fix the period with P = xu.with_unit(pm.Constant('P', P), u.day). However, it does not seen to run and the error is as follows. Should I how to deal with this error?

 File "/home/jlin/xinglong/RV/fixP.py", line 110, in
samples = joker.rejection_sample(data,prior_samples=100_000, max_posterior_samples=128)
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/thejoker/thejoker.py", line 248, in rejection_sample
samples = rejection_sample_helper(
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/thejoker/utils.py", line 294, in wrapper
raise e
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/thejoker/utils.py", line 292, in wrapper
func_return = func(*args, **kwargs)
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/thejoker/multiproc_helpers.py", line 231, in rejection_sample_helper
samples = make_full_samples(joker_helper, prior_samples_file, pool,
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/thejoker/multiproc_helpers.py", line 139, in make_full_samples
results = run_worker(make_full_samples_worker, pool, prior_samples_file,
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/thejoker/multiproc_helpers.py", line 48, in run_worker
for res in pool.map(worker, tasks):
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/thejoker/multiproc_helpers.py", line 126, in make_full_samples_worker
raw_samples, _ = joker_helper.batch_get_posterior_samples(batch,
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/thejoker/src/fast_likelihood.cpython-38-x86_64-linux-gnu.so", line 464, in thejoker.src.fast_likelihood.CJokerHelper.batch_get_posterior_samples
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/thejoker/src/fast_likelihood.cpython-38-x86_64-linux-gnu.so", line 537, in thejoker.src.fast_likelihood.CJokerHelper.batch_get_posterior_samples

builtins.ValueError: cannot reshape array of size 0 into shape (0,newaxis)

Thanks!

[adrn]

What versions of cython, numpy, and thejoker do you have installed?

[AstroJLin]

Hi,
Here, I use cython (0.29.21), numpy(1.21.1) and thejoker (1.1). In fact, If I run a small example that you give me, it can work. However, when I use my data to run with P = xu.with_unit(pm.Constant('P', P), u.day), this error will appear. So, I do not know what is happen？

Thanks!

[adrn]

What type is your variable P? You should probably name it something else, like period or _P since you redefine it as a model parameter. What if you try:

P = xu.with_unit(pm.Constant('P', float(_P)), u.day)

?

[AstroJLin]

Hi,

Although have named it something else, it still reports this error. Here I attach my data and code, could you help me to see what is happen?

Thanks!

from astropy.io import ascii
from astropy.time import Time
import astropy.units as u
import matplotlib.pyplot as plt
#%matplotlib inline
import numpy as np

import pymc3 as pm
import exoplanet.units as xu
import exoplanet as xo
from exoplanet.distributions import Angle
import corner
import arviz as az

import thejoker as tj
import emcee
seed = 42
rnd = np.random.default_rng(seed=seed)

tbl = ascii.read(
"""
2458183.36136016 -97.2 3.1
2458183.37733315 -37.7 4.9
2458183.39330614 24.5 3.8
2458187.34557366 -7.6 3
2458187.36154661 -89.7 2.5
2458187.37821403 -155.2 2
2458187.39418697 -193.7 2
2459295.33498398 141.1 2.2
2459295.35165123 153.2 2.1
2459295.36762401 147 1.9
2457091.33443744 -187.9 1.3
2458902.39203349 143 0.6
""",
names=['BJD', 'rv', 'rv_err'])
tbl['rv'].unit = u.km/u.s
tbl['rv_err'].unit = u.km/u.s

data = tj.RVData(
t=Time(tbl['BJD'], format='jd', scale='tcb'),
rv=u.Quantity(tbl['rv']),
rv_err=u.Quantity(tbl['rv_err']))

data.plot(markersize=3)
plt.show()
period = 0.2556686
with pm.Model() as model:
e = xu.with_unit(pm.Constant('e', 0),u.one)
P = xu.with_unit(pm.Constant('P', float(period)), u.day)
prior = tj.JokerPrior.default(
# Range of periods to consider
#P_min=0.1u.day, P_max=1.u.day,
sigma_K0=170u.km/u.s, # scale of the prior on semiamplitude, K
sigma_v=50u.km/u.s, # std dev of the prior on the systemic velocity, v0
#v0_offsets=[dv0_1],
pars={'e':e,'P':P}
)
joker = tj.TheJoker(prior)
samples = joker.rejection_sample(data,prior_samples=100_000, max_posterior_samples=128)
_ = tj.plot_rv_curves(samples, data=data)

[adrn]

Oh interesting - there is something weird about the way pm.Constant() is storing the value (it keeps it internally as an int, so the period becomes 0, which is invalid).

Try changing the definitions of e and P to:

    e = xu.with_unit(pm.Deterministic('e', tt.as_tensor_variable(0)), u.one)
    P = xu.with_unit(pm.Deterministic('P', tt.as_tensor_variable(period)), u.day)

Also added your example to the bottom here:
https://gist.github.com/56707a370423d883a33bcc9049bdbf7a

[AstroJLin]

Hi Adrian,
Thank you very much for your help and reply.

I have been able to fixed the period after using your suggestions. And then, when I try to run the pmx.sample() in the mcmc, it throw an error (builtins.TypeError: ('Wrong number of dimensions: expected 0, got 1 with shape (256,).', 'Container name "K"')).
I do not know what is causing this error. Could you give me some suggetions? Here I attach code and this error.

###
with pm.Model():
    e = xu.with_unit(pm.Deterministic('e', tt.constant(0)),
    P = xu.with_unit(pm.Deterministic('P', tt.constant(period)), u.day)
    prior_mcmc = tj.JokerPrior.default(
        sigma_K0=100*u.km/u.s,
        sigma_v=50*u.km/u.s,
        pars={'e': e, 'P':P}
    )
    joker_mcmc = tj.TheJoker(prior_mcmc, random_state=rnd)
    mcmc_init = joker_mcmc.setup_mcmc(data, samples)
    print (mcmc_init)


(an error??? )    trace = pmx.sample(
        tune=500, draws=1000,
        start=mcmc_init,
        random_seed=seed,
        cores=1, chains=2)  
    
print (az.summary(trace, var_names=prior_mcmc.par_names))
mcmc_samples = joker_mcmc.trace_to_samples(trace, data=data)
mcmc_samples.wrap_K()
print (mcmc_samples.mean().tbl['P', 'e', 'omega', 'M0', 's', 'K', 'v0'])
df = mcmc_samples.tbl.to_pandas()
df = df.drop(columns=['e', 's', 'omega','M0'])
_ = corner.corner(df,quantiles=(0.16,0.5,0.84),show_titles=True, title_kwargs={"fontsize": 12})
plt.show()
fig = plt.figure(figsize=(10, 8))
_ = tj.plot_phase_fold(mcmc_samples.median(), data)
plt.xlabel('Phase')
plt.ylabel('Radial Velocity (km/s)')
plt.show()
###

###The error???
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/theano/link/basic.py", line 82, in set
self.storage[0] = self.type.filter_inplace(
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/theano/graph/type.py", line 95, in filter_inplace
raise NotImplementedError()

builtins.NotImplementedError

During handling of the above exception, another exception occurred:

 File "/home/jlin/xinglong/RV/fixP.py", line 161, in
trace = pmx.sample(
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/pymc3_ext/sampling/sampling.py", line 97, in sample
return pm.sample(draws=draws, tune=tune, model=model, step=step, **kwargs)
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/pymc3/sampling.py", line 435, in sample
check_start_vals(start, model)
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/pymc3/util.py", line 234, in check_start_vals
initial_eval = model.check_test_point(test_point=elem)
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/pymc3/model.py", line 1384, in check_test_point
{RV.name: np.round(RV.logp(test_point), round_vals) for RV in self.basic_RVs},
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/pymc3/model.py", line 1384, in
{RV.name: np.round(RV.logp(test_point), round_vals) for RV in self.basic_RVs},
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/pymc3/model.py", line 1561, in call
return self.f(**point)
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/theano/compile/function/types.py", line 897, in call
self[k] = arg
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/theano/compile/function/types.py", line 559, in setitem
self.value[item] = value
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/theano/compile/function/types.py", line 515, in setitem
s.value = value
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/theano/link/basic.py", line 86, in set
self.storage[0] = self.type.filter(value, **kwargs)
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/theano/tensor/type.py", line 181, in filter
raise TypeError(

builtins.TypeError: ('Wrong number of dimensions: expected 0, got 1 with shape (256,).', 'Container name "K"')

[adrn]

For your data, thejoker returns many samples -- i.e., your samples variable has 256 posterior samples. When you run setup_mcmc(), you should only pass one sample. Here you can do that by just picking one at random, e.g.,

mcmc_init = joker_mcmc.setup_mcmc(data, samples[0])

or by taking the MAP or maximum likelihood sample, e.g.,

loglike = samples.ln_unmarginalized_likelihood(data)
mcmc_init = joker_mcmc.setup_mcmc(data, samples[loglike.argmax()])

[AstroJLin]

Hi Adrian,
Thank you very much for your reply again.

Here, I have two questions.

First, when I using the maximum likelihood sample, I do not use the median MCMC sample to fold the data and plot residuals relative to our inferred RV model. Here, an error is report as follow:

 File "/home/jlin/xinglong/RV/fixP.py", line 202, in
_ = tj.plot_phase_fold(mcmc_samples.median(), data)
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/thejoker/plot.py", line 272, in plot_phase_fold
orbit = sample.get_orbit()
File "/home/jlin/anaconda3/envs/thejokertest/lib/python3.8/site-packages/thejoker/samples.py", line 258, in get_orbit
raise ValueError("You must specify an index when the number "
builtins.ValueError: You must specify an index when the number of samples is >1 (here, it's 2000)

which means that I should set up the number of samples ( mcmc_samples.median()[?]). However, I do not know which samples corresponding to the median MCMC sample?

Second, If I do not try to fix the period, I do not need to set up the maximum likelihood sample as an example from the Thompson et al. 2019. So, I want to know the reason for this.

Thanks!