AM mini
In [1]:
Copied!
from IPython.display import Image
import astroemperor as emp
import numpy as np
import multiprocessing
np.random.seed(1234)
sim = emp.Simulation()
sim.load_data('HIP8923') # folder read from /datafiles/
sim.cores__ = multiprocessing.cpu_count()
from IPython.display import Image
import astroemperor as emp
import numpy as np
import multiprocessing
np.random.seed(1234)
sim = emp.Simulation()
sim.load_data('HIP8923') # folder read from /datafiles/
sim.cores__ = multiprocessing.cpu_count()
I couldnt grab the terminal size. Trying with pandas... Terminal size with pandas successful! ~~ Simulation Successfully Initialized ~~ Reading data from HIP8923_1_cor07.vels Reading data from HIP8923_1_cor14.vels Reading data from HIP8923_1_cor14T.vels
In [2]:
Copied!
ntemps, nwalkers, nsweeps, nsteps = 6, 128, 500, 1
sim.set_engine('reddemcee')
sim.engine_config['setup'] = [ntemps, nwalkers, nsweeps, nsteps]
# Custom adaptation options for reddemcee engine
sim.engine_config['tsw_history'] = True
sim.engine_config['smd_history'] = True
sim.engine_config['adapt_tau'] = 100
sim.engine_config['adapt_nu'] = 1.0
sim.engine_config['adapt_mode'] = 'SAR'
# Manual ladder
sim.engine_config['betas'] = list(np.geomspace(1, 0.001, ntemps))
# sim.engine_config['betas'][-1] = 0 # Optionally, for exact evidence estimation
sim.run_config['burnin'] = 0.5 # its in niter
if True:
sim.starmass = 1.0086627 # 0.9686578 1.0486716 # FLAME
sim.starmass_err = 0.0400069
sim.instrument_names_RV = ['Cor07', 'Cor14', 'Cor14T']
# Model options
sim.jitter_prargs = [0, 5] # Normal(0, 10)
sim.eccentricity_prargs = [0, 0.1] # Normal(0, 0.1)
sim.acceleration = 0 # no linear trend
# Offsets and Jitters manual boundaries
for i in range(1, 4):
sim.add_condition([f'Offset {i}', 'limits', [-50, 50]])
sim.add_condition([f'Jitter {i}', 'limits', [1e-4, 15]])
# we constrain a bit the planetary parameters
sim.add_condition(['Period 1', 'limits', [5100, 5200]])
sim.add_condition(['Amplitude 1', 'limits', [45, 50]])
sim.add_condition(['Eccentricity 1', 'limits', [0, 0.02]])
sim.add_condition(['Phase 1', 'limits', [4.9, 5.0]])
sim.add_condition(['Longitude 1', 'limits', [0.9, 0.94]])
sim.add_condition(['Inclination 1', 'limits', [0.40, 0.45]])
sim.add_condition(['Omega 1', 'limits', [1.70, 1.74]])
# We also constrain astrometric offsets
for offset in ['RA', 'DE', 'PLX', 'pm RA', 'pm DE']:
sim.add_condition([f'Offset {offset}', 'limits', [-1, 1]]) # usually -10, 10
# Finally, we set some sensible priors for Hipparcos and Gaia
sim.add_condition(['Jitter Gaia', 'limits', [1, 5]]) # fractional jitter
sim.add_condition(['Jitter Gaia', 'prior', 'Normal'])
sim.add_condition(['Jitter Gaia', 'prargs', [1, 0.1]])
sim.add_condition(['Jitter Hipp', 'limits', [0, 5]]) # additive jitter
if True:
sim.add_condition(['Period 1', 'init_pos', [5150, 5170]])
sim.add_condition(['Amplitude 1', 'init_pos', [48, 49.5]])
sim.add_condition(['Eccentricity 1', 'init_pos', [0.005, 0.014]])
sim.add_condition(['Phase 1', 'init_pos', [4.94, 4.96]])
sim.add_condition(['Longitude 1', 'init_pos', [0.91, 0.93]])
sim.add_condition(['Inclination 1', 'init_pos', [0.43, 0.44]])
sim.add_condition(['Omega 1', 'init_pos', [1.71, 1.73]])
sim.add_condition(['Offset 1', 'init_pos', [-30.3, -30.2]])
sim.add_condition(['Offset 2', 'init_pos', [-38.2, -38.1]])
sim.add_condition(['Offset 3', 'init_pos', [29.05, 29.10]])
sim.add_condition(['Jitter 1', 'init_pos', [1.710, 1.720]])
sim.add_condition(['Jitter 2', 'init_pos', [2.81, 2.83]])
sim.add_condition(['Jitter 3', 'init_pos', [0.93, 0.95]])
sim.add_condition(['Offset RA', 'init_pos', [0.600, 0.608]])
sim.add_condition(['Offset DE', 'init_pos', [0.400, 0.416]])
sim.add_condition(['Offset PLX', 'init_pos', [0.016, 0.018]])
sim.add_condition(['Offset pm RA', 'init_pos', [0.21, 0.23]])
sim.add_condition(['Offset pm DE', 'init_pos', [-0.240, -0.235]])
sim.add_condition(['Jitter Gaia', 'init_pos', [1.13, 1.14]]) # 0.01
sim.add_condition(['Jitter Hipparcos', 'init_pos', [1.01, 1.02]]) # 0.01
ntemps, nwalkers, nsweeps, nsteps = 6, 128, 500, 1
sim.set_engine('reddemcee')
sim.engine_config['setup'] = [ntemps, nwalkers, nsweeps, nsteps]
# Custom adaptation options for reddemcee engine
sim.engine_config['tsw_history'] = True
sim.engine_config['smd_history'] = True
sim.engine_config['adapt_tau'] = 100
sim.engine_config['adapt_nu'] = 1.0
sim.engine_config['adapt_mode'] = 'SAR'
# Manual ladder
sim.engine_config['betas'] = list(np.geomspace(1, 0.001, ntemps))
# sim.engine_config['betas'][-1] = 0 # Optionally, for exact evidence estimation
sim.run_config['burnin'] = 0.5 # its in niter
if True:
sim.starmass = 1.0086627 # 0.9686578 1.0486716 # FLAME
sim.starmass_err = 0.0400069
sim.instrument_names_RV = ['Cor07', 'Cor14', 'Cor14T']
# Model options
sim.jitter_prargs = [0, 5] # Normal(0, 10)
sim.eccentricity_prargs = [0, 0.1] # Normal(0, 0.1)
sim.acceleration = 0 # no linear trend
# Offsets and Jitters manual boundaries
for i in range(1, 4):
sim.add_condition([f'Offset {i}', 'limits', [-50, 50]])
sim.add_condition([f'Jitter {i}', 'limits', [1e-4, 15]])
# we constrain a bit the planetary parameters
sim.add_condition(['Period 1', 'limits', [5100, 5200]])
sim.add_condition(['Amplitude 1', 'limits', [45, 50]])
sim.add_condition(['Eccentricity 1', 'limits', [0, 0.02]])
sim.add_condition(['Phase 1', 'limits', [4.9, 5.0]])
sim.add_condition(['Longitude 1', 'limits', [0.9, 0.94]])
sim.add_condition(['Inclination 1', 'limits', [0.40, 0.45]])
sim.add_condition(['Omega 1', 'limits', [1.70, 1.74]])
# We also constrain astrometric offsets
for offset in ['RA', 'DE', 'PLX', 'pm RA', 'pm DE']:
sim.add_condition([f'Offset {offset}', 'limits', [-1, 1]]) # usually -10, 10
# Finally, we set some sensible priors for Hipparcos and Gaia
sim.add_condition(['Jitter Gaia', 'limits', [1, 5]]) # fractional jitter
sim.add_condition(['Jitter Gaia', 'prior', 'Normal'])
sim.add_condition(['Jitter Gaia', 'prargs', [1, 0.1]])
sim.add_condition(['Jitter Hipp', 'limits', [0, 5]]) # additive jitter
if True:
sim.add_condition(['Period 1', 'init_pos', [5150, 5170]])
sim.add_condition(['Amplitude 1', 'init_pos', [48, 49.5]])
sim.add_condition(['Eccentricity 1', 'init_pos', [0.005, 0.014]])
sim.add_condition(['Phase 1', 'init_pos', [4.94, 4.96]])
sim.add_condition(['Longitude 1', 'init_pos', [0.91, 0.93]])
sim.add_condition(['Inclination 1', 'init_pos', [0.43, 0.44]])
sim.add_condition(['Omega 1', 'init_pos', [1.71, 1.73]])
sim.add_condition(['Offset 1', 'init_pos', [-30.3, -30.2]])
sim.add_condition(['Offset 2', 'init_pos', [-38.2, -38.1]])
sim.add_condition(['Offset 3', 'init_pos', [29.05, 29.10]])
sim.add_condition(['Jitter 1', 'init_pos', [1.710, 1.720]])
sim.add_condition(['Jitter 2', 'init_pos', [2.81, 2.83]])
sim.add_condition(['Jitter 3', 'init_pos', [0.93, 0.95]])
sim.add_condition(['Offset RA', 'init_pos', [0.600, 0.608]])
sim.add_condition(['Offset DE', 'init_pos', [0.400, 0.416]])
sim.add_condition(['Offset PLX', 'init_pos', [0.016, 0.018]])
sim.add_condition(['Offset pm RA', 'init_pos', [0.21, 0.23]])
sim.add_condition(['Offset pm DE', 'init_pos', [-0.240, -0.235]])
sim.add_condition(['Jitter Gaia', 'init_pos', [1.13, 1.14]]) # 0.01
sim.add_condition(['Jitter Hipparcos', 'init_pos', [1.01, 1.02]]) # 0.01
In [3]:
Copied!
sim.plot_posteriors['temps'] = [0]
sim.plot_gaussian_mixtures['plot'] = False
sim.plot_trace['plot'] = False
sim.plot_histograms['plot'] = False
sim.plot_posteriors['temps'] = [0]
sim.plot_gaussian_mixtures['plot'] = False
sim.plot_trace['plot'] = False
sim.plot_histograms['plot'] = False
In [ ]:
Copied!
In [4]:
Copied!
sim.autorun(1, 1)
sim.autorun(1, 1)
Offset block added, OffsetBlock Jitter block added, JitterBlock Keplerian block added, AstrometryKeplerianBlock 1 AstrometryOffset block added, AstrometryOffsetBlock AstrometryJitter block added, AstrometryJitterBlock Condition applied: Parameter Period 1 attribute limits set to [5100, 5200] Condition applied: Parameter Period 1 attribute init_pos set to [5150, 5170] Condition applied: Parameter Amplitude 1 attribute limits set to [45, 50] Condition applied: Parameter Amplitude 1 attribute init_pos set to [48, 49.5] Condition applied: Parameter Phase 1 attribute limits set to [4.9, 5.0] Condition applied: Parameter Phase 1 attribute init_pos set to [4.94, 4.96] Condition applied: Parameter Eccentricity 1 attribute limits set to [0, 0.02] Condition applied: Parameter Eccentricity 1 attribute init_pos set to [0.005, 0.014] Condition applied: Parameter Longitude 1 attribute limits set to [0.9, 0.94] Condition applied: Parameter Longitude 1 attribute init_pos set to [0.91, 0.93] Condition applied: Parameter Inclination 1 attribute limits set to [0.4, 0.45] Condition applied: Parameter Inclination 1 attribute init_pos set to [0.43, 0.44] Condition applied: Parameter Omega 1 attribute limits set to [1.7, 1.74] Condition applied: Parameter Omega 1 attribute init_pos set to [1.71, 1.73] Condition applied: Parameter Offset 1 attribute limits set to [-50, 50] Condition applied: Parameter Offset 1 attribute init_pos set to [-30.3, -30.2] Condition applied: Parameter Offset 2 attribute limits set to [-50, 50] Condition applied: Parameter Offset 2 attribute init_pos set to [-38.2, -38.1] Condition applied: Parameter Offset 3 attribute limits set to [-50, 50] Condition applied: Parameter Offset 3 attribute init_pos set to [29.05, 29.1] Condition applied: Parameter Jitter 1 attribute limits set to [0.0001, 15] Condition applied: Parameter Jitter 1 attribute init_pos set to [1.71, 1.72] Condition applied: Parameter Jitter 2 attribute limits set to [0.0001, 15] Condition applied: Parameter Jitter 2 attribute init_pos set to [2.81, 2.83] Condition applied: Parameter Jitter 3 attribute limits set to [0.0001, 15] Condition applied: Parameter Jitter 3 attribute init_pos set to [0.93, 0.95] Condition applied: Parameter Offset DE attribute limits set to [-1, 1] Condition applied: Parameter Offset DE attribute init_pos set to [0.4, 0.416] Condition applied: Parameter Offset PLX attribute limits set to [-1, 1] Condition applied: Parameter Offset PLX attribute init_pos set to [0.016, 0.018] Condition applied: Parameter Offset pm DE attribute limits set to [-1, 1] Condition applied: Parameter Offset pm DE attribute init_pos set to [-0.24, -0.235] Condition applied: Parameter Jitter Hipp attribute limits set to [0, 5] Condition applied: Parameter Jitter Gaia attribute limits set to [1, 5] Condition applied: Parameter Jitter Gaia attribute prior set to Normal Condition applied: Parameter Jitter Gaia attribute prargs set to [1, 0.1] Condition applied: Parameter Jitter Gaia attribute init_pos set to [1.13, 1.14] ~~ Setup Info ~~ Current Engine is reddemcee 1.0 Number of cores is 24 Save location is datalogs/HIP8923/run_5/k1 Dynamical Criteria is None Posterior fit method is Gaussian Mixtures Limits constrain method is range Model Selection method is BIC ~~ Automatically Saving ~~ Logger : ✔ Samples : ✘ Posteriors : ✔ Likelihoods : ✔ Plots: Posteriors : ✔ Plots: Keplerian Model : ✔ Plots: Gaussian Mixture : ✘ Plots: Parameter Histograms : ✘ Plots: Corner : ✔ ~~ Pre-Run Info ~~ Parameter Prior Limits ---------------- ------------------ -------------- Period 1 ~𝓤 (5100, 5200) [5100, 5200] Amplitude 1 ~𝓤 (45, 50) [45, 50] Phase 1 ~𝓤 (4.9, 5.0) [4.9, 5] Eccentricity 1 ~𝓝 (0.0, 0.1) [0, 0.02] Longitude 1 ~𝓤 (0.9, 0.94) [0.9, 0.94] Inclination 1 ~I (0.4, 0.45) [0.4, 0.45] Omega 1 ~𝓤 (1.7, 1.74) [1.7, 1.74] ---------------- ------------------ -------------- Offset 1 ~𝓤 (-50, 50) [-50, 50] Offset 2 ~𝓤 (-50, 50) [-50, 50] Offset 3 ~𝓤 (-50, 50) [-50, 50] ---------------- ------------------ -------------- Jitter 1 ~𝓝 (0, 5) [0, 15] Jitter 2 ~𝓝 (0, 5) [0, 15] Jitter 3 ~𝓝 (0, 5) [0, 15] ---------------- ------------------ -------------- Offset RA* ~𝓤 (-10.0, 10.0) [-10, 10] Offset DE ~𝓤 (-1, 1) [-1, 1] Offset PLX ~𝓤 (-1, 1) [-1, 1] Offset pm RA* ~𝓤 (-10.0, 10.0) [-10, 10] Offset pm DE ~𝓤 (-1, 1) [-1, 1] ---------------- ------------------ -------------- Jitter Hipp ~𝓤 (0, 5) [0, 5] Jitter Gaia ~𝓝 (1.0, 0.1) [1, 5] Math for AstrometryKeplerianBlock 1: K⋅sin(I)⋅(cos(ν(t,P,𝜙,e)+𝜔 )+e⋅cos(𝜔 ))|₁ Math for OffsetBlock: γ₀|ᵢ Math for JitterBlock: 𝝈ᵢ Math for AstrometryOffsetBlock: γ₀|ᵢ Math for AstrometryJitterBlock: Generating Samples --------------------------------------------------------------------------- AttributeError Traceback (most recent call last) File /IronCrane/reddtea/pip_packages/astroemperor/docs/tutorials/temp_script_HIP8923_run_5.py:978 973 return p0 976 p1 = test_init() --> 978 sampler = reddemcee.PTSampler(nwalkers, 979 ndim, 980 my_likelihood, 981 my_prior, 982 ntemps=ntemps, 983 pool=mypool, 984 backend=my_backend, 985 betas=betas, 986 tsw_history=True, 987 smd_history=True, 988 adapt_tau=100, 989 adapt_nu=1.0, 990 adapt_mode="SAR" 991 ) 994 sampler.D_ = np.array([1.00000e+02, 5.00000e+00, 1.00000e-01, 2.00000e-02, 4.00000e-02, 995 5.00000e-02, 4.00000e-02, 1.00000e+02, 1.00000e+02, 1.00000e+02, 996 1.49999e+01, 1.49999e+01, 1.49999e+01, 2.00000e+01, 2.00000e+00, 997 2.00000e+00, 2.00000e+01, 2.00000e+00, 5.00000e+00, 4.00000e+00]) 999 def run_thing(): File ~/anaconda3/envs/emperor-test/lib/python3.12/site-packages/reddemcee/sampler.py:99, in PTSampler.__init__(self, nwalkers, ndim, log_like, log_prior, betas, ntemps, pool, moves, loglargs, loglkwargs, logpargs, logpkwargs, backend, vectorize, blobs_dtype, smd_history, tsw_history, adapt_tau, adapt_nu, adapt_mode, parameter_names) 97 self.config_adaptation_rate = adapt_nu 98 self.config_adaptation_mode = adapt_mode ---> 99 self.select_adjustment(self.config_adaptation_mode) 101 self.z_ngrid = 10001 102 self.z_nsim = 1000 File ~/anaconda3/envs/emperor-test/lib/python3.12/site-packages/reddemcee/sampler.py:38, in LadderAdaptation.select_adjustment(self, option) 37 def select_adjustment(self, option=0): ---> 38 self.adjustment_fn = getattr(self, f'calc_adjust{option}') AttributeError: 'PTSampler' object has no attribute 'calc_adjustSAR'
mv: cannot stat 'HIP8923_run_5.h5': No such file or directory mv: cannot stat 'HIP8923_run_5_0.h5': No such file or directory mv: cannot stat 'HIP8923_run_5_1.h5': No such file or directory mv: cannot stat 'HIP8923_run_5_2.h5': No such file or directory mv: cannot stat 'HIP8923_run_5_3.h5': No such file or directory mv: cannot stat 'HIP8923_run_5_4.h5': No such file or directory mv: cannot stat 'HIP8923_run_5_5.h5': No such file or directory
--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[4], line 1 ----> 1 sim.autorun(1, 1) File /IronCrane/reddtea/pip_packages/astroemperor/src/astroemperor/emp.py:2489, in Simulation.autorun(self, k_start, k_end) 2485 self._autorun_add_blocks() 2487 self._stat_holder_update() -> 2489 self.run() # also loads 2491 self.postprocess() 2493 self.run_plot_routines() File /IronCrane/reddtea/pip_packages/astroemperor/src/astroemperor/emp.py:2540, in Simulation.run(self) 2537 self._run_clean() 2539 # Load the sampler -> 2540 self._load_sampler() 2542 self.time_run = time.time() - time_run_init File /IronCrane/reddtea/pip_packages/astroemperor/src/astroemperor/emp.py:775, in emp_scribe._load_sampler(self) 773 def _load_sampler(self): 774 self.debug(f'run : _load_sampler() | {time.time()-self.time_init}') --> 775 getattr(self, f'_load_sampler_{self.engine__.__name__}')() File /IronCrane/reddtea/pip_packages/astroemperor/src/astroemperor/emp.py:790, in emp_scribe._load_sampler_reddemcee(self) 785 reader = PTHDFBackend(f'{self.saveplace}/restore/backends/{self.backend_name}.h5') 787 reader.backends = [HDFBackend_plus(f'{self.saveplace}/restore/backends/{self.backend_name}_{t}.h5') for t in range(ntemps)] --> 790 self.betas = list(reader.get_last_sample().betas) 791 self.engine_config["betas"] = self.betas 792 self.engine_config["ntemps"] = len(self.betas) File /IronCrane/reddtea/pip_packages/reddemcee/src/reddemcee/hdf.py:517, in PTHDFBackend.get_last_sample(self) 515 """Access the most recent sample in the chain""" 516 if (not self.initialized) or self.iteration <= 0: --> 517 raise AttributeError( 518 "you must run the sampler with " 519 "'store == True' before accessing the " 520 "results" 521 ) 523 states = PTState([b.get_last_sample() for b in self.backends]) 524 states.random_state = self.random_state AttributeError: you must run the sampler with 'store == True' before accessing the results
In [ ]:
Copied!
z_ti = np.round(sim.sampler.get_evidence_ti(discard=sim.reddemcee_discard)[:2], 3)
z_ss = np.round(sim.sampler.get_evidence_ss(discard=sim.reddemcee_discard), 3)
z_hy = np.round(sim.sampler.get_evidence_hybrid(discard=sim.reddemcee_discard), 3)
print(f'Thermodynamic Integration logZ = {z_ti[0]} ±{z_ti[1]}')
print(f'Stepping Stones logZ = {z_ss[0]} ±{z_ss[1]}')
print(f'Hybrid Method logZ = {z_hy[0]} ±{z_hy[1]}')
z_ti = np.round(sim.sampler.get_evidence_ti(discard=sim.reddemcee_discard)[:2], 3)
z_ss = np.round(sim.sampler.get_evidence_ss(discard=sim.reddemcee_discard), 3)
z_hy = np.round(sim.sampler.get_evidence_hybrid(discard=sim.reddemcee_discard), 3)
print(f'Thermodynamic Integration logZ = {z_ti[0]} ±{z_ti[1]}')
print(f'Stepping Stones logZ = {z_ss[0]} ±{z_ss[1]}')
print(f'Hybrid Method logZ = {z_hy[0]} ±{z_hy[1]}')
Thermodynamic Integration logZ = -289.327 ±20.011 Stepping Stones logZ = -286.358 ±0.317 Hybrid Method logZ = -299.35 ±20.098
Thermodynamic Integration logZ = -308.941 ±2.876
Stepping Stones logZ = -308.473 ±0.497
Hybrid Method logZ = -302.723 ±3.042
In [ ]:
Copied!
# inputs
input_theta = repr(np.asarray(sim.fit_max).tolist())
# inputs
input_theta = repr(np.asarray(sim.fit_max).tolist())
In [ ]:
Copied!
In [ ]:
Copied!
import matplotlib.pyplot as pl
from astroemperor.emp_model import ReddModel
from astroemperor.globals import _PLATFORM_SYSTEM, _CORES
from astroemperor.block_repo import *
from astroemperor.qol_utils import *
from astroemperor.math_utils import *
from reddcolors import Palette
import matplotlib.pyplot as pl
from astroemperor.emp_model import ReddModel
from astroemperor.globals import _PLATFORM_SYSTEM, _CORES
from astroemperor.block_repo import *
from astroemperor.qol_utils import *
from astroemperor.math_utils import *
from reddcolors import Palette
In [ ]:
Copied!
sim.model.get_attr_block('type_')
sim.model.get_attr_block('type_')
Out[ ]:
['Keplerian', 'Offset', 'Jitter', 'AstrometryOffset', 'AstrometryJitter']
In [ ]:
Copied!
Out[ ]:
True
In [ ]:
Copied!
options = {'saveloc':sim.saveplace,
'format':'pdf',
'figsize':(10, 8),
'label_fontsize':22,
'title_fontsize':12,
'tick_labelsize':18,
'line_width':2,
'major_tick_length':8,
'minor_tick_length':4,
'marker_size':8,
'scatter_size':20,
'legend_fs':14,
'fm_frame_lw':3,
}
def plot_astrometry(model_, input_theta, options={}):
if True:
saveplace_ = options['saveloc']
savefmt_ = options['format']
saveloc_plot_ = saveplace_ + f'/plots/models/astrometry_model.{savefmt_}'
figsize_ = options['figsize']
label_fontsize = options['label_fontsize']
title_fontsize = options['title_fontsize']
tick_labelsize = options['tick_labelsize']
line_width = options['line_width']
major_tick_length = options['major_tick_length']
minor_tick_length = options['minor_tick_length']
marker_size = options['marker_size']
scatter_size = options['scatter_size']
legend_fs = options['legend_fs']
fm_frame_lw = options['fm_frame_lw']
tail = ''
model_.model_script_no = 0
temp_script_name = f'temp_astrometry_plot{tail}.py'
temp_script_loc = f'{saveplace_}' # TODO: change
temp_plot_script = f'{temp_script_loc}/{temp_script_name}'
dependencies = model_.get_dependencies().tolist()
constants = model_.get_constants()
with open(temp_plot_script, 'w') as f:
if True:
# ASTROMETRY EXCLUSIVE
f.write('''
# dependencies
if True:
import matplotlib.gridspec as gridspec
import matplotlib.colors as plc
from matplotlib.ticker import MaxNLocator, FixedLocator
from matplotlib.cm import plasma
from matplotlib.colors import Normalize
from astropy.time import Time as AstroTime
from matplotlib.patches import Polygon, Patch
from matplotlib.lines import Line2D
from reddcolors import Palette
from astroemperor.canvas import hex2rgb, mk_cmap
''')
# --------
f.write(open(get_support('init_reddemcee.scr')).read())
for d in dependencies:
f.write(f'''
{d}''')
f.write('''
''')
for c in constants:
f.write(f'''{c} = {constants[c]}
''')
model_script_name = model_.write_model(loc=saveplace_,
tail=tail)
f.write(open(model_script_name).read())
# astrometry exclusive functions
f.write(open(get_support('astrometry/plot.00')).read())
f.write(f"""
theta = np.array({input_theta})
""")
kc = 1
for b in model_:
if b.type_=='Keplerian':
f.write(f"""
slice_kep{kc} = {repr(b.slice)}
""")
kc += 1
elif b.type_=='AstrometryOffset':
f.write(f"""
theta_am_off = theta[{repr(b.slice)}]
""")
elif b.type_=='AstrometryJitter':
f.write(f"""
theta_am_jit = theta[{repr(b.slice)}]
""")
# astrometry IAD
f.write("""
ca = 0
ca_err = 0
ca_hipp = 0
go = 0
nodes = 0
tt_cats = AM_catalogs_times[1:] - common_t
plx0 = AM_PLX_ref - theta_am_off[2]
df_hipp_b = data_iad_hipp.copy().groupby("BJD").agg({
'IORB': 'mean',
'EPOCH': 'mean',
'PARF': 'mean',
'CPSI': 'mean',
'SPSI': 'mean',
'RES': 'mean',
'SRES': 'mean'}).reset_index()
""")
kc = 1
for b in model_:
if b.type_=='Keplerian':
f.write(f"""
P{kc} = theta[slice_kep{kc}][0]
tt{kc} = np.linspace(0, P{kc}, 5000)
ca += calc_astro_new_t(theta[slice_kep{kc}], plx0, tt=tt{kc})
ca_hipp += calc_astro_hipp_binned(theta[slice_kep{kc}], plx0, df_hipp_b)
go += calc_astro_new_t(theta[slice_kep{kc}], plx0, tt=tt_cats)
nodes += calc_astro_nodes(theta[slice_kep{kc}], plx0)
""")
kc += 1
f.write(open(get_support('astrometry/plot.01')).read())
f.write("""
# astrometry inputs
def get_projections(iad_df):
ca_err = 0
t_hip = iad_df['BJD'].values
cpsi = iad_df['CPSI'].values
spsi = iad_df['SPSI'].values
r_al = iad_df['RES'].values
sig_al = iad_df['SRES'].values
hipp_iad_time = t_hip - common_t
""")
kc = 1
for b in model_:
if b.type_=='Keplerian':
f.write(f"""
ca_err += calc_astro_new_t(theta[slice_kep{kc}], plx0, tt=hipp_iad_time)
""")
kc += 1
f.write("""
x_mod = ca_err[0]
y_mod = ca_err[1]
u = np.column_stack([cpsi, spsi])
obs = np.column_stack([x_mod, y_mod]) + r_al[:, None]*u
sres = np.sqrt(sig_al**2 + theta_am_jit[0]**2)
seg_lo = obs - sres[:, None]*u
seg_hi = obs + sres[:, None]*u
segments_long = np.stack([seg_lo, seg_hi], axis=1)
seg_lo = obs - sig_al[:, None]*u
seg_hi = obs + sig_al[:, None]*u
segments_short = np.stack([seg_lo, seg_hi], axis=1)
return obs, segments_long, segments_short
proj = get_projections(data_iad_hipp0)
dra_hipp_iad = proj[0][:, 0]
dde_hipp_iad = proj[0][:, 1]
proj_b = get_projections(df_hipp_b)
dra_hipp_biad = proj_b[0][:, 0]
dde_hipp_biad = proj_b[0][:, 1]
binned_lines_long = proj_b[1]
binned_lines_shor = proj_b[2]
colors_bin = plasma(Normalize(vmin=df_hipp_b['BJD'].min(),
vmax=df_hipp_b['BJD'].max())(df_hipp_b['BJD']))
""")
f.write(f"""
figsize_ = {figsize_}
saveloc_ = '{saveloc_plot_}'
label_fontsize = {label_fontsize}
title_fontsize = {title_fontsize}
tick_labelsize = {tick_labelsize}
line_width = {line_width}
major_tick_length = {major_tick_length}
minor_tick_length = {minor_tick_length}
marker_size = {marker_size}
scatter_size = {scatter_size}
legend_fs = {legend_fs}
fm_frame_lw = {fm_frame_lw}
""")
f.write(open(get_support('astrometry/astrometry.plot')).read())
os.system(f'ipython {temp_script_loc}/{temp_script_name}')
pass
options = {'saveloc':sim.saveplace,
'format':'pdf',
'figsize':(10, 8),
'label_fontsize':22,
'title_fontsize':12,
'tick_labelsize':18,
'line_width':2,
'major_tick_length':8,
'minor_tick_length':4,
'marker_size':8,
'scatter_size':20,
'legend_fs':14,
'fm_frame_lw':3,
}
def plot_astrometry(model_, input_theta, options={}):
if True:
saveplace_ = options['saveloc']
savefmt_ = options['format']
saveloc_plot_ = saveplace_ + f'/plots/models/astrometry_model.{savefmt_}'
figsize_ = options['figsize']
label_fontsize = options['label_fontsize']
title_fontsize = options['title_fontsize']
tick_labelsize = options['tick_labelsize']
line_width = options['line_width']
major_tick_length = options['major_tick_length']
minor_tick_length = options['minor_tick_length']
marker_size = options['marker_size']
scatter_size = options['scatter_size']
legend_fs = options['legend_fs']
fm_frame_lw = options['fm_frame_lw']
tail = ''
model_.model_script_no = 0
temp_script_name = f'temp_astrometry_plot{tail}.py'
temp_script_loc = f'{saveplace_}' # TODO: change
temp_plot_script = f'{temp_script_loc}/{temp_script_name}'
dependencies = model_.get_dependencies().tolist()
constants = model_.get_constants()
with open(temp_plot_script, 'w') as f:
if True:
# ASTROMETRY EXCLUSIVE
f.write('''
# dependencies
if True:
import matplotlib.gridspec as gridspec
import matplotlib.colors as plc
from matplotlib.ticker import MaxNLocator, FixedLocator
from matplotlib.cm import plasma
from matplotlib.colors import Normalize
from astropy.time import Time as AstroTime
from matplotlib.patches import Polygon, Patch
from matplotlib.lines import Line2D
from reddcolors import Palette
from astroemperor.canvas import hex2rgb, mk_cmap
''')
# --------
f.write(open(get_support('init_reddemcee.scr')).read())
for d in dependencies:
f.write(f'''
{d}''')
f.write('''
''')
for c in constants:
f.write(f'''{c} = {constants[c]}
''')
model_script_name = model_.write_model(loc=saveplace_,
tail=tail)
f.write(open(model_script_name).read())
# astrometry exclusive functions
f.write(open(get_support('astrometry/plot.00')).read())
f.write(f"""
theta = np.array({input_theta})
""")
kc = 1
for b in model_:
if b.type_=='Keplerian':
f.write(f"""
slice_kep{kc} = {repr(b.slice)}
""")
kc += 1
elif b.type_=='AstrometryOffset':
f.write(f"""
theta_am_off = theta[{repr(b.slice)}]
""")
elif b.type_=='AstrometryJitter':
f.write(f"""
theta_am_jit = theta[{repr(b.slice)}]
""")
# astrometry IAD
f.write("""
ca = 0
ca_err = 0
ca_hipp = 0
go = 0
nodes = 0
tt_cats = AM_catalogs_times[1:] - common_t
plx0 = AM_PLX_ref - theta_am_off[2]
df_hipp_b = data_iad_hipp.copy().groupby("BJD").agg({
'IORB': 'mean',
'EPOCH': 'mean',
'PARF': 'mean',
'CPSI': 'mean',
'SPSI': 'mean',
'RES': 'mean',
'SRES': 'mean'}).reset_index()
""")
kc = 1
for b in model_:
if b.type_=='Keplerian':
f.write(f"""
P{kc} = theta[slice_kep{kc}][0]
tt{kc} = np.linspace(0, P{kc}, 5000)
ca += calc_astro_new_t(theta[slice_kep{kc}], plx0, tt=tt{kc})
ca_hipp += calc_astro_hipp_binned(theta[slice_kep{kc}], plx0, df_hipp_b)
go += calc_astro_new_t(theta[slice_kep{kc}], plx0, tt=tt_cats)
nodes += calc_astro_nodes(theta[slice_kep{kc}], plx0)
""")
kc += 1
f.write(open(get_support('astrometry/plot.01')).read())
f.write("""
# astrometry inputs
def get_projections(iad_df):
ca_err = 0
t_hip = iad_df['BJD'].values
cpsi = iad_df['CPSI'].values
spsi = iad_df['SPSI'].values
r_al = iad_df['RES'].values
sig_al = iad_df['SRES'].values
hipp_iad_time = t_hip - common_t
""")
kc = 1
for b in model_:
if b.type_=='Keplerian':
f.write(f"""
ca_err += calc_astro_new_t(theta[slice_kep{kc}], plx0, tt=hipp_iad_time)
""")
kc += 1
f.write("""
x_mod = ca_err[0]
y_mod = ca_err[1]
u = np.column_stack([cpsi, spsi])
obs = np.column_stack([x_mod, y_mod]) + r_al[:, None]*u
sres = np.sqrt(sig_al**2 + theta_am_jit[0]**2)
seg_lo = obs - sres[:, None]*u
seg_hi = obs + sres[:, None]*u
segments_long = np.stack([seg_lo, seg_hi], axis=1)
seg_lo = obs - sig_al[:, None]*u
seg_hi = obs + sig_al[:, None]*u
segments_short = np.stack([seg_lo, seg_hi], axis=1)
return obs, segments_long, segments_short
proj = get_projections(data_iad_hipp0)
dra_hipp_iad = proj[0][:, 0]
dde_hipp_iad = proj[0][:, 1]
proj_b = get_projections(df_hipp_b)
dra_hipp_biad = proj_b[0][:, 0]
dde_hipp_biad = proj_b[0][:, 1]
binned_lines_long = proj_b[1]
binned_lines_shor = proj_b[2]
colors_bin = plasma(Normalize(vmin=df_hipp_b['BJD'].min(),
vmax=df_hipp_b['BJD'].max())(df_hipp_b['BJD']))
""")
f.write(f"""
figsize_ = {figsize_}
saveloc_ = '{saveloc_plot_}'
label_fontsize = {label_fontsize}
title_fontsize = {title_fontsize}
tick_labelsize = {tick_labelsize}
line_width = {line_width}
major_tick_length = {major_tick_length}
minor_tick_length = {minor_tick_length}
marker_size = {marker_size}
scatter_size = {scatter_size}
legend_fs = {legend_fs}
fm_frame_lw = {fm_frame_lw}
""")
f.write(open(get_support('astrometry/astrometry.plot')).read())
os.system(f'ipython {temp_script_loc}/{temp_script_name}')
pass
In [ ]:
Copied!
plot_astrometry(sim.model, options)
plot_astrometry(sim.model, options)
I couldnt grab the terminal size. Trying with pandas... Terminal size with pandas successful!
In [ ]:
Copied!
--------------------------------------------------------------------------- NameError Traceback (most recent call last) Cell In[112], line 1 ----> 1 fdra NameError: name 'fdra' is not defined
In [ ]:
Copied!
Out[ ]:
'[5129.3759460102465, 48.52058448943847, 4.962099055382451, 0.009488241687722818, 0.9160063146697622, 0.42975136261506125, 1.7275335239960639, -29.725833513622447, -38.14485268785246, 29.135197557121206, 1.6771067552714802, 2.780029892978177, 1.0278469834717894, 0.627673375937684, 0.392951029972419, 0.017360874670579764, 0.20394204743578137, -0.25076696212102134, 1.3880360325317327, 1.0568185744846854]'
In [ ]:
Copied!
In [ ]:
Copied!
In [ ]:
Copied!
In [ ]:
Copied!