Virtual RM#
mica2 provides a vmap mode to do reverberation mapping analysis with a virtual driving light curve. This mode applies
in cases where the dirving light curve cannot be chosen or the driving light curve has a poor qaulity that is not suitable to act
as the dirving one.
To this end, mica2 assume the virtual drving light curve follows a DRW process with a variation amplitude (\(\sigma\)) of 0.1
and has a time lag of zero with respect to the first light curve of the input data. The remaining analysis is trival and
similar to the normal modes.
The input data should have a format as:
# 1
# 0:171:269
56690.6100 3.4270 0.0640 % lc1, 171 lines
56691.5400 3.5450 0.0650
...
56864.8600 4.3310 0.0740
56865.9200 4.7080 0.0780
56698.3570 2.1900 0.0560 % lc2, 269 lines
56699.5590 2.2000 0.0580
...
56830.1490 2.3000 0.0650
56830.4200 2.2900 0.0660
Note the in the second line, the first number is set to 0, meaning that the light curve is virtual and does not have data.
For the exectuable binary version, the typical parameter file looks like:
#
# lines starting with "#" are regarded as comments and are neglected
# if want to turn on the line, remove the beginning "#"
#==============================================================
FileDir ./
DataFile data/sim_data.txt
TypeModel 2 # 0: general model
# 1: pmap, photometric RM
# 2: vmap, use a virtual driving light curve.
TypeTF 0 # 0: Gaussian
# 1: Top-hat
MaxNumberSaves 2000 # number of steps
FlagUniformVarParams 0 # if each data set has the same variability parameters
FlagUniformTranFuns 0 # if each data set has the same tf parameters
FlagLongtermTrend 0 # if long-term trending
FlagLagPositivity 0 # if enable tf at positive lags
NumCompLow 1
NumCompUpp 1
FlagConSysErr 0
FlagLineSysErr 0
TypeLagPrior 0 # type of lag prior for each Gaussians/tophats.
# 0, limit0 < lag0 < lag1 < lag2 <... < limit1
#
# 1, limit0 + 0*width < lag0 < limit0 + 1*width
# limit0 + 1*width < lag1 < limit0 + 2*width
# ...
# width = (limit1 - limit0)/num_comp
#
# 2, lags fixed at specific values, no limit on Guassian sigma/tophat width
# lag0 = limit0 + 0*dlag
# lag1 = limit0 + 1*dlag
# ...
# dlag = (limit1 - limit0)/(num_comp-1)
#
# 3, lags fixed at specific values
# Gaussian sigma ranges at (dlag/2, dlag), tophat wdith=dlag/2
# lag0 = limit0 + 0*dlag
# lag1 = limit0 + 1*dlag
# ...
# dlag = (limit1 - limit0)/(num_comp-1)
# better to set a large mumber of components
StrLagPrior [0:10:10:50] # valid if TypeLagPrior==4
# format: [lag1_1:lag1_2:lag2_1:lag2_2...]
# "LagLimitLow" and "LagLimitUpp" no longer applicable
For python version, mica provide a module vmap callable as follows.
from mpi4py import MPI
import numpy as np
import pymica
import matplotlib.pyplot as plt
# initiate MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
# load data
if rank == 0:
lc0 = np.empty(0) # virtual light curve, empty
lc1 = np.loadtxt("mcg08_g.txt")
lc2 = np.loadtxt("mcg08_r.txt")
# make a data dict
data_input = {"set1":[lc0, lc1, lc2]}
else:
data_input = None
data_input = comm.bcast(data_input, root=0)
#create a model
#there are two ways
#1) one way from the param file
#model = pymica.gmodel(param_file="param/param_input")
#2) the ohter way is through the setup function
model = pymica.vmap()
# use Gaussians
model.setup(data=data_input, type_tf='gaussian', lag_limit=[-2, 5], number_component=[1, 1], max_num_saves=2000)
# or use tophats
#model.setup(data=data_input, type_tf='tophat', lag_limit=[0, 100], number_component=[1, 1], max_num_saves=2000)
#the full arguments are
#model.setup(data_file=None, data=None,
# type_tf='gaussian', max_num_saves=2000,
# flag_uniform_var_params=False, flag_uniform_tranfuns=False,
# flag_trend=0, flag_lag_posivity=False,
# lag_limit=[0, 100], number_component=[1, 1],
# width_limit=[0.1, 100],
# flag_con_sys_err=False, flag_line_sys_err=False,
# type_lag_prior=0, lag_prior=[[0, 50]],
# num_particles=2, thread_steps_factor=2,
# new_level_interval_factor=2, save_interval_factor=2,
# lam=10, beta=100, ptol=0.1,
# max_num_levels=0)
#run mica
model.run()
#posterior run, only re-generate posterior samples, do not run MCMC
# model.post_run()
#do decomposition for the cases of multiple components
#model.decompose()
# plot results
if rank == 0:
model.plot_results() # plot results
model.post_process() # generate plots for the properties of MCMC sampling
Here is an example for vmap analysis. The data is extracted from Fausnaugh et al. 2018, ApJ, 854, 10.
An examplary result of MICA2 analysis with vmap mode. The right topmost panel shows the virtual driving light curve.#