hpparvi / muscat2_transit_pipeline Goto Github PK
View Code? Open in Web Editor NEWMuSCAT2 transit analysis pipeline
License: GNU General Public License v3.0
MuSCAT2 transit analysis pipeline
License: GNU General Public License v3.0
For some datasets, ta.apply_normalized_limits cause an error only when we put apply=True in the argument. When we put apply=False, it shows plots without any problem.
The error message is as follows.
ValueError Traceback (most recent call last)
in
1 #this is run if obvious outliers are seen in raw flux plots above
----> 2 ta.apply_normalized_limits(2, lower=0.9, upper=1.1, apply=True, plot=True, npoly=3, iterations=5)
~/github/MuSCAT2_transit_pipeline/muscat2ta/transitanalysis.py in apply_normalized_limits(self, iapt, lower, upper, plot, apply, npoly, iterations, erosion)
148 def apply_normalized_limits(self, iapt: int, lower: float = -inf, upper: float = inf, plot: bool = True,
149 apply: bool = True, npoly: int = 0, iterations: int = 5, erosion: int = 0) -> None:
--> 150 self.lpf.apply_normalized_limits(iapt, lower, upper, plot, apply, npoly, iterations, erosion)
151
152 def downsample(self, exptime: float) -> None:
~/github/MuSCAT2_transit_pipeline/muscat2ta/m2lpf.py in apply_normalized_limits(self, iapt, lower, upper, plot, apply, npoly, iterations, erosion)
487 if apply:
488 fluxes.append(ft[mask])
--> 489 self.times[ipb] = self.times[ipb][mask]
490 self.ofluxes[ipb] = self.ofluxes[ipb][mask, :]
491 for icid in range(self.cids.shape[1]):
ValueError: could not broadcast input array from shape (209) into shape (212)
We have several options to choose from
Either a custom-built GP module or one of the two latter options is probably the best way to go, since it looks like the development of george has stalled.
TOI02431.01 observed on 210213
for ph,band in zip(ta.phs,ta.passbands):
fig = ph.plot_raw(8, # plot raw flux of star ids from 0 to 7
#ylim=(0.92, 1.05), # y-limit of the plots
figsize=(13,5), # no need to change
);
fig.suptitle(f'{band}-band')
returns
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-23-18c757bb5d25> in <module>
1 for ph,band in zip(ta.phs,ta.passbands):
----> 2 fig = ph.plot_raw(8, # plot raw flux of star ids from 0 to 7
3 #ylim=(0.92, 1.05), # y-limit of the plots
4 figsize=(13,5), # no need to change
5 );
~/github/MuSCAT2_transit_pipeline/muscat2ph/phdata.py in plot_raw(self, nstars, ylim, ncols, figsize)
369 fig, axs = subplots(nrows, ncols, figsize=figsize, sharex='all', sharey='all')
370 nflux = self.flux / self.flux.median('mjd')
--> 371 aids = abs(nflux.diff('mjd')).median('mjd').argmin('aperture')
372 for i,(ax,iapt) in enumerate(zip(axs.flat, aids[:nstars])):
373 apt = apertures[iapt]
~/miniconda3/envs/py3/lib/python3.8/site-packages/xarray/core/dataarray.py in argmin(self, dim, axis, keep_attrs, skipna)
3989 Dimensions without coordinates: y
3990 """
-> 3991 result = self.variable.argmin(dim, axis, keep_attrs, skipna)
3992 if isinstance(result, dict):
3993 return {k: self._replace_maybe_drop_dims(v) for k, v in result.items()}
~/miniconda3/envs/py3/lib/python3.8/site-packages/xarray/core/variable.py in argmin(self, dim, axis, keep_attrs, skipna)
2274 DataArray.argmin, DataArray.idxmin
2275 """
-> 2276 return self._unravel_argminmax("argmin", dim, axis, keep_attrs, skipna)
2277
2278 def argmax(
~/miniconda3/envs/py3/lib/python3.8/site-packages/xarray/core/variable.py in _unravel_argminmax(self, argminmax, dim, axis, keep_attrs, skipna)
2195 # Return int index if single dimension is passed, and is not part of a
2196 # sequence
-> 2197 return self.reduce(
2198 argminmax_func, dim=dim, axis=axis, keep_attrs=keep_attrs, skipna=skipna
2199 )
~/miniconda3/envs/py3/lib/python3.8/site-packages/xarray/core/variable.py in reduce(self, func, dim, axis, keep_attrs, keepdims, allow_lazy, **kwargs)
1638 )
1639 if axis is not None:
-> 1640 data = func(input_data, axis=axis, **kwargs)
1641 else:
1642 data = func(input_data, **kwargs)
~/miniconda3/envs/py3/lib/python3.8/site-packages/xarray/core/duck_array_ops.py in f(values, axis, skipna, **kwargs)
329
330 try:
--> 331 return func(values, axis=axis, **kwargs)
332 except AttributeError:
333 if not is_duck_dask_array(values):
~/miniconda3/envs/py3/lib/python3.8/site-packages/xarray/core/nanops.py in nanargmin(a, axis)
95
96 module = dask_array if isinstance(a, dask_array_type) else nputils
---> 97 return module.nanargmin(a, axis=axis)
98
99
~/miniconda3/envs/py3/lib/python3.8/site-packages/xarray/core/nputils.py in f(values, axis, **kwargs)
225 result = bn_func(values, axis=axis, **kwargs)
226 else:
--> 227 result = getattr(npmodule, name)(values, axis=axis, **kwargs)
228
229 return result
<__array_function__ internals> in nanargmin(*args, **kwargs)
~/miniconda3/envs/py3/lib/python3.8/site-packages/numpy/lib/nanfunctions.py in nanargmin(a, axis)
497 mask = np.all(mask, axis=axis)
498 if np.any(mask):
--> 499 raise ValueError("All-NaN slice encountered")
500 return res
501
ValueError: All-NaN slice encountered
Hi Hannu, here's the code snippet that fetches the most recent TOI releases table:
import os
import pandas as pd
def get_tois(clobber=True, outdir='~/data/', verbose=False,
remove_FP=True, remove_known_planets=False):
"""Download TOI list from TESS Alert/TOI Release.
Parameters
----------
clobber : bool
re-download table and save as csv file
outdir : str
download directory location
verbose : bool
print texts
Returns
-------
d : pandas.DataFrame
TOI table as dataframe
"""
dl_link = 'https://exofop.ipac.caltech.edu/tess/download_toi.php?sort=toi&output=csv'
fp = os.path.join(outdir,'TOIs.csv')
if not os.path.exists(outdir):
os.makedirs(outdir)
if not os.path.exists(fp) or clobber:
d = pd.read_csv(dl_link)
#remove False Positives
if remove_FP:
d = d[d['TFOPWG Disposition']!='FP']
print('TOIs with TFPWG disposition==FP are removed.\n')
if remove_known_planets:
planet_keys = ['WASP','SWASP','HAT','HATS','KELT','QATAR','K2','Kepler']
keys = []
for key in planet_keys:
idx = ~np.array(d['Comments'].str.contains(key).tolist(),dtype=bool)
d = d[idx]
if idx.sum()>0:
keys.append(key)
print(f'{keys} planets are removed.\n')
d.to_csv(fp, index=False)
print(f'Saved: {fp}\n')
else:
d = pd.read_csv(fp)
#remove False Positives
if remove_FP:
d = d[d['TFOPWG Disposition']!='FP']
if verbose:
print(f'Loaded: {fp}')
return d.sort_values('TOI')
def get_toi(toi=None, tic=None, clobber=True, outdir='~/data/', verbose=False):
"""Query TOI from TOI list
Parameters
----------
tic : int
TIC id
toi : float
TOI id
clobber : bool
re-download csv file
outdir : str
csv path
verbose : bool
print texts
Returns
-------
q : pandas.DataFrame
TOI match else None
"""
df = get_tois(clobber=clobber, verbose=verbose, outdir=outdir)
if toi:
if isinstance(toi, int):
toi = float(str(toi)+'.01')
else:
planet = str(toi).split('.')[1]
assert len(planet)==2, 'use pattern: TOI.01'
idx = df['TOI'].isin([toi])
elif tic:
idx = df['TIC ID'].isin([tic])
q=df.loc[idx]
#return if empty, else continue
if len(q)==0:
raise ValueError('TOI not found!')
q.index = q['TOI'].values
if verbose:
print('Data from TOI Release:\n')
columns = ['Period (days)','Epoch (BJD)','Duration (hours)','Depth (ppm)','Comments']
print(f'{q[columns].T}\n')
if q['TFOPWG Disposition'].isin(['FP']).any():
print('\nTFOPWG disposition is a False Positive!\n')
return q.sort_values(by='TOI', ascending=True)
MuSCAT2_transit_pipeline/muscat2ta/lpf.py
Line 71 in 9684de4
If the star hosts multiple planets, and the host star name is erroneously passed as "target" to positional argument in m2fit script, the exocat query will return an array. Line 71 will then throw an error.
Here's a simple example in Ipython. Hat-p-44 has 2 planets.
import exodata
exocat = exodata.load_db_from_url()
target='hat-p-44'
system=exocat.searchPlanet(target)
system
[Planet('HAT-P-44 b'), Planet('HAT-P-44 c')]
Probably a band-aid solution is to add:
assert isinstance(self.planet, exodata.astroclasses.Planet)
before L71, or maybe require the planet letter explicitly in the "target" positional argument in m2fit (and --target-catalog-name argument in m2photometry too?).
MuSCAT2_transit_pipeline/muscat2ph/phdata.py
Line 305 in e4cb993
I propose to have more informative text labels so apertures instead of aperture indexes are shown in the plot.
def plot_raw(self, nstars, ylim=(0.85, 1.15), ncols=4):
apertures = self._ds.aperture.data
date = self.mjd[0].iso.split(' ')[0]
nrows = int(ceil(nstars / ncols))
fig, axs = subplots(nrows, ncols, figsize=self.figsize, sharex=True, sharey=True)
nflux = self.flux / self.flux.median('mjd')
aids = abs(nflux.diff('mjd')).median('mjd').argmin('aperture')
for i,(ax,iapt) in enumerate(zip(axs.flat, aids[:nstars])):
apt = apertures[iapt]
ax.plot(self.mjd.value, nflux[:,i,iapt], 'k')
ax.text(0.05, 0.05, 'star id={}\nr={} pix'.format(i,int(apt)),
transform=ax.transAxes, color='b', size='large')
setp(axs, ylim=ylim)
sb.despine(fig)
pl.suptitle(date, fontsize='large')
fig.tight_layout()
return fig
Hi Hannu,
I cannot figure out where exactly and why the kernel crashes whenever I run plot_possible_blends in TFOPAnalysis.
I tried to debug line by line with pdb but each line in the code doesn't throw an error until after the last line when the kernel suddenly stops working.
Please note that we have not been able to run this part ever since even after updating with the most recent dev. Can you suggest a way how to better debug this issue?
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(320)plot_possible_blends()
-> m_excl = ones(self.distances.size, bool)
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(321)plot_possible_blends()
-> if self.excluded_stars:
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(324)plot_possible_blends()
-> m = (self.distances < max_separation) & m_excl
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(325)plot_possible_blends()
-> sids = where(m)[0]
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(326)plot_possible_blends()
-> stars = sids[argsort(self.distances[m])]
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(327)plot_possible_blends()
-> nstars = len(stars)
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(328)plot_possible_blends()
-> nrows = int(ceil(nstars / ncols))
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(329)plot_possible_blends()
-> naxs = nrows * ncols
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(330)plot_possible_blends()
-> nempty = naxs - nstars
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(332)plot_possible_blends()
-> phs = self.phs if pbs is None else [self.phs[i] for i in pbs]
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(333)plot_possible_blends()
-> passbands = self.passbands if pbs is None else [self.passbands[i] for i in pbs]
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(334)plot_possible_blends()
-> plotname = self._dres.joinpath(f"{self.ticname}_20{self.date}_MuSCAT2_raw_lcs.pdf")
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(335)plot_possible_blends()
-> pdf = PdfPages(plotname) if save else None
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(337)plot_possible_blends()
-> for pb, ph in zip(passbands, phs):
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(338)plot_possible_blends()
-> fig, axs = subplots(nrows, ncols, figsize=(figwidth, nrows * axheight), constrained_layout=True,
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(339)plot_possible_blends()
-> sharex='all')
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(340)plot_possible_blends()
-> for i, istar in enumerate(tqdm(stars, leave=False)):
(Pdb) n
Internal StopIteration
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(340)plot_possible_blends()
-> for i, istar in enumerate(tqdm(stars, leave=False)):
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(342)plot_possible_blends()
-> setp(axs[:,0], ylabel='Normalized flux')
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(343)plot_possible_blends()
-> [axs.flat[naxs - i - 1].remove() for i in range(nempty)]
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(344)plot_possible_blends()
-> fig.suptitle(f"{self.ticname} 20{self.date} {pb}-band MuSCAT2 relative normalised fluxes.",
(Pdb) n
> /home/muscat/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py(345)plot_possible_blends()
-> size='large')
(Pdb) n
Currently, ta.lpf.ps
prints:
0 |G| tc [-inf .. inf],
1 |G| pr [0.00 .. inf],
2 |G| rho [0.00 .. inf],
3 |G| b [0.00 .. 1.00],
It would be useful if it prints the actual respective prior and bounds (or have another method to do this). This is a common misconception in our group. I always have to remind them to check ta.lpf.ps[i].prior
for each i to just make sure what is the actual prior used.
Hello, Hannu.
I am trying the analysis of Gj436b on 180318.
Images of frameid=1 to 17 are observer's mistakes.
So, I renamed these files to *.fit_removed
.
and I executed m2astrometry --start-frame=18 --overwrite
, m2photometry --start-frame=18
, but it shows aperture of frameid=1 over FOV of frameid=18.
Do you know any solutions?
I analyzed data of TOI01634.01 observed on 210213.
#optimize the parameter set using differential evolution
ta.optimize(500)
fig, axs = ta.plot_fit('de');
returns
Optimizing the model: 0%
0/500 [00:00<?, ?it/s]
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-48-da8ec6d1d9c3> in <module>
1 #optimize the parameter set using differential evolution
----> 2 ta.optimize(500)
3
4 fig, axs = ta.plot_fit('de');
~/github/MuSCAT2_transit_pipeline/muscat2ta/transitanalysis.py in optimize(self, niter, pop, plot_convergence, plot_lc)
207
208 def optimize(self, niter: int = 1000, pop: ndarray = None, plot_convergence: bool = True, plot_lc: bool = False):
--> 209 self.lpf.optimize_global(niter, self.npop, pop, label='Optimizing the model')
210 self.pv = self.lpf.de.minimum_location
211
~/miniconda3/envs/py3/lib/python3.8/site-packages/pytransit/lpf/logposteriorfunction.py in optimize_global(self, niter, npop, population, label, leave, plot_convergence, use_tqdm, plot_parameters)
142 else:
143 self.de._population[:, :] = population
--> 144 for _ in tqdm(self.de(niter), total=niter, desc=label, leave=leave, disable=(not use_tqdm)):
145 pass
146
~/miniconda3/envs/py3/lib/python3.8/site-packages/tqdm/notebook.py in __iter__(self, *args, **kwargs)
251 def __iter__(self, *args, **kwargs):
252 try:
--> 253 for obj in super(tqdm_notebook, self).__iter__(*args, **kwargs):
254 # return super(tqdm...) will not catch exception
255 yield obj
~/miniconda3/envs/py3/lib/python3.8/site-packages/tqdm/std.py in __iter__(self)
1164
1165 try:
-> 1166 for obj in iterable:
1167 yield obj
1168 # Update and possibly print the progressbar.
~/miniconda3/envs/py3/lib/python3.8/site-packages/pytransit/utils/de.py in _eval_vfun(self, ngen)
221 popt, fitt = self._trial_pop, self._trial_fit
222
--> 223 fitc[:] = self.m * self.minfun(self._population)
224
225 for igen in range(ngen):
~/miniconda3/envs/py3/lib/python3.8/site-packages/pytransit/utils/de.py in __call__(self, x)
248
249 def __call__(self, x):
--> 250 return self.f(x, *self.args, **self.kwargs)
~/miniconda3/envs/py3/lib/python3.8/site-packages/pytransit/lpf/logposteriorfunction.py in lnposterior(self, pv)
118
119 def lnposterior(self, pv):
--> 120 lnp = self.lnprior(pv) + self.lnlikelihood(pv)
121 return where(isfinite(lnp), lnp, -inf)
122
~/miniconda3/envs/py3/lib/python3.8/site-packages/pytransit/lpf/lpf.py in lnlikelihood(self, pvp)
444
445 for lnlikelihood in self._lnlikelihood_models:
--> 446 lnl += lnlikelihood(pvp, fmodel)
447 return lnl
448
~/github/MuSCAT2_transit_pipeline/muscat2ta/m2lpf.py in __call__(self, pvp, model)
206 return lnlike_logistic_v1d(self.lpf.ofluxa, flux_m, err, self.lpf.lcids)
207 else:
--> 208 return lnlike_logistic_v(self.lpf.relative_flux(pvp), flux_m, err, self.lpf.lcids)
209
210 class M2LPF(BaseLPF):
~/github/MuSCAT2_transit_pipeline/muscat2ta/m2lpf.py in relative_flux(self, pv)
781 return self.ofluxa
782 else:
--> 783 return self.target_flux(pv) / self.reference_flux(pv)
784
785 def target_flux(self, pv):
~/github/MuSCAT2_transit_pipeline/muscat2ta/m2lpf.py in target_flux(self, pv)
790 off = zeros((p.shape[0], self.timea.size))
791 for i, sl in enumerate(self.lcslices):
--> 792 off[:, sl] = self.target_fluxes[i][:, p].T
793 return squeeze(off)
794
ValueError: could not broadcast input array from shape (200,0) into shape (200,239)
In the .nc file the distribution of the central time of the transit has values with a fraction of a day, the full BJD of the central time of the transit is not stored anywhere.
Hello, Hannu.
I would like to analyze MuSCAT1 data, however I couldn't execute m2organize
.
So I organized by my self. And then, I executed m2photometry
.
m2photometry
required FITS filenames starting from MCT2
, so I renamed.
After that, m2photometry
also required target-name_stars.fits
.
I couldn't make this FITS file by myself, so I gave up at here.
I would like to analyze in your pipeline. In addition to it, Norio wanted to analyze MuSCAT1 data for TFOP in your pipeline.
If you can, Could you make a new option for analysis of MuSCAT1?
I am sorry for my many issues.
Sincerely,
Taku
When I put with_contamination=True as an argument in ta=TFOPAnalysis(), every process goes well but at last part, ta.plot_final_fit() fails. The error message is as below.
ValueError Traceback (most recent call last)
in
----> 1 fig = ta.plot_final_fit();
~/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py in plot_final_fit(self, figwidth, save, close)
492 fig.add_axes((0.03, 0.32, 0.96, 0.001), facecolor='k', xticks=[], yticks=[])
493 lpf.plot_posteriors(fig=fig,
--> 494 gridspec=dict(top=0.30, bottom=0.05, left=0.1, right=0.95, wspace=0.03, hspace=0.3))
495 fig.add_axes((0.03, 0.01, 0.96, 0.001), facecolor='k', xticks=[], yticks=[])
496 plotname = self._dres.joinpath(f"{self.ticname}_20{self.date}_MuSCAT2_transit_fit.pdf")
~/github/MuSCAT2_transit_pipeline/muscat2ta/m2lpf.py in plot_posteriors(self, figsize, fig, gridspec, plot_k)
880
881 x = arange(self.nlc)
--> 882 plot_estimates(x, 1e6 * df[k2cols].quantile([0.5, 0.16, 0.84, 0.025, 0.975]).values, axs[0, 0])
883 if self.toi is not None:
884 axs[0, 0].axhline(self.toi.depth[0], c='k', ls='--')
~/github/MuSCAT2_transit_pipeline/muscat2ta/m2lpf.py in plot_estimates(x, p, ax, bwidth, color)
875 k2cols = [c for c in df.columns if 'k2' in c]
876 def plot_estimates(x, p, ax, bwidth=0.8, color='k'):
--> 877 ax.bar(x, p[4, :] - p[3, :], bwidth, p[3, :], alpha=0.25, fc=color)
878 ax.bar(x, p[2, :] - p[1, :], bwidth, p[1, :], alpha=0.25, fc=color)
879 [ax.plot((xx - 0.47 * bwidth, xx + 0.47 * bwidth), (pp[[0, 0]]), 'k') for xx, pp in zip(x, p.T)]
~/miniconda3/lib/python3.6/site-packages/matplotlib/init.py in inner(ax, data, *args, **kwargs)
1599 def inner(ax, *args, data=None, **kwargs):
1600 if data is None:
-> 1601 return func(ax, *map(sanitize_sequence, args), **kwargs)
1602
1603 bound = new_sig.bind(ax, *args, **kwargs)
~/miniconda3/lib/python3.6/site-packages/matplotlib/axes/_axes.py in bar(self, x, height, width, bottom, align, **kwargs)
2373 x, height, width, y, linewidth = np.broadcast_arrays(
2374 # Make args iterable too.
-> 2375 np.atleast_1d(x), height, width, y, linewidth)
2376
2377 # Now that units have been converted, set the tick locations.
~/.local/lib/python3.6/site-packages/numpy-1.16.3-py3.6-linux-x86_64.egg/numpy/lib/stride_tricks.py in broadcast_arrays(*args, **kwargs)
257 args = [np.array(_m, copy=False, subok=subok) for _m in args]
258
--> 259 shape = _broadcast_shape(*args)
260
261 if all(array.shape == shape for array in args):
~/.local/lib/python3.6/site-packages/numpy-1.16.3-py3.6-linux-x86_64.egg/numpy/lib/stride_tricks.py in _broadcast_shape(*args)
191 # use the old-iterator because np.nditer does not handle size 0 arrays
192 # consistently
--> 193 b = np.broadcast(*args[:32])
194 # unfortunately, it cannot handle 32 or more arguments directly
195 for pos in range(32, len(args), 31):
ValueError: shape mismatch: objects cannot be broadcast to a single shape
We have a lot of TFOP datasets that have very close neighboring stars. I wanted to check from which star the signal comes from by limiting the range of apertures in the pipeline. For example,
ta = TFOPAnalysis(target='toi06240.01',
date=230811,
tid=0,
cids=[4,5,6],
dataroot='./photometry',
aperture_lims=(0, 2), # small aper
# aperture_lims=(3, 7), # big aper
)
Here I superpose the selected apertures with the centroids to visualize if the aperture fully encapsulates the centroid spread.
import matplotlib.pyplot as pl
from matplotlib.patches import Circle
import astropy.units as u
import pandas as pd
def plot_apertures(starid, band_id=0, unit="arcsec", rlim=7, pix_scale=0.44, ax=None):
d = ta.distances[starid]*u.arcmin.to(u.arcsec)
apers_in_pix = ta.phs[band_id]._ds.aperture.data[:rlim]
centroids = ta.phs[band_id]._ds.centroid[:,starid,]
dx,dy = centroids.T*pix_scale
if unit=="arcmin":
dx*=u.arcsec.to(u.arcmin)
dy*=u.arcsec.to(u.arcmin)
d *=u.arcsec.to(u.arcmin)
if ax is None:
fig, ax = pl.subplots()
ax.set_title(f'star ID={starid} (d={d:.1f} {unit} = {round(d/pix_scale)} pix)');
ax.scatter(dx+d,dy+d)
#get centroid median and place apertures there
mx, my = np.median(dx+d), np.median(dy+d)
ax.plot(mx, my,'ko', label=f"dx,dy=({mx:.1f},{my:.1f}) {unit}", ms=10, lw=5)
ax.plot(d, d,'rx', ms=10, lw=5)
for index,p in enumerate(apers_in_pix):
ls = '--' if index%2 else '-'
r = p*pix_scale
if unit=="arcmin":
r *= u.arcsec.to(u.arcmin)
circle = Circle((mx, my), r, facecolor='none',
edgecolor='k', linewidth=3, alpha=0.5, ls=ls)
ax.add_patch(circle)
ax.set_xlabel(f'dx ({unit})')
ax.set_ylabel(f'dy ({unit})')
ax.legend()
return ax
#plot the target and the nearest neighbor
idx = pd.Series(ta.distances).sort_values().argsort().index.tolist()
fig, ax = pl.subplots()
plot_apertures(idx[0], ax=ax, rlim=3);
plot_apertures(idx[1], ax=ax, rlim=3);
In this case, is it correct that aperture_lims=(0, 2)
can be used to avoid flux contamination from the nearby star?
This is related to the email I sent to Felipe and Hannu on Aug 20. When I choose a narrow aperture_lims=(5,6)
for example, the parameter set specifies the target aperture as
13 |L| tap U(a = 0.0, b = 0.999) [ 0.00 .. 1.00],
However, after freezing the photometry, the frozen_apertures
property does not take into account the index offset in self.phs[0]._ds.aperture[self.taps]
so it will incorrectly report the aperture e.g. for index 0 instead of index 5.
Hi Hannu,
We have one TOI dataset that has only g-band data available. After optimization, ta.plot_fit('de') failed because the figure axes expect single index instead of a tuple, it seems. May be a work around is to flatten the axes when there is only a single band to plot. Same error occurs for other plotting methods such as ta.plot_final_fit()
--------------------------------------------------------------------------
IndexError Traceback (most recent call last)
<ipython-input-19-99f3968dcbc5> in <module>
----> 1 fig = ta.plot_fit('de');
~/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py in plot_fit(self, model, figsize, save, plot_priors)
295
296 def plot_fit(self, model: str = 'de', figsize: tuple = (13, 8), save=False, plot_priors=True):
--> 297 fig, axs = self.lpf.plot_light_curves(model=model, figsize=figsize)
298 npb = self.lpf.npb
299 if plot_priors:
~/github/MuSCAT2_transit_pipeline/muscat2ta/m2lpf.py in plot_light_curves(self, model, figsize, fig, gridspec, ylim_transit, ylim_residuals, bin_width)
653 for i, sl in enumerate(self.lcslices):
654 t = self.timea[sl]
--> 655 axs[1, i].plot(t, self.ofluxa[sl], '.', alpha=0.5)
656 axs[1, i].plot(t, fm[0][sl], 'k', lw=2)
657 axs[2, i].plot(t, self.ofluxa[sl] / bl[0][sl], '.', alpha=0.5)
IndexError: too many indices for array
MuSCAT2_transit_pipeline/muscat2ta/lpf.py
Line 23 in 9684de4
Hi Hannu, I'd like to request to get a copy of this module (currently in private repo?) in order to completely run the pipeline. Thank you.
Preliminary dumb version of aperture size optimisation.
Parallelisation scheme will depend on the chosen MCMC sampler. Can be either
Also, parallelisation (other than OpenMP) may not be necessary given the small amount of data per analysis.
i.e. TransitAnalysis.load() should parse the xarray.Dataset produced by TransitAnalysis.save() and populate the relevant TransitAnalysis attributes, instead of simply returning the xarray.Dataset
TFOPAnalysis
throws the error shown below when it is initialized without r
in the passband
argument, and also when no r
-band dataset is available in the photometry folder, as in the case of TOI 1605 taken on Sep. 2.
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-12-0efaa7fc2cab> in <module>
4 CIDS=[1,2,3,4] #comparison id(s); can be guesses for now
5
----> 6 ta = TFOPAnalysis(target=OBJ,
7 date=DATE,
8 tid=TID,
~/github/MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py in __init__(self, target, date, tid, cids, dataroot, nlegendre, npop, mjd_start, mjd_end, excluded_mjd_ranges, aperture_lims, passbands, use_opencl, with_transit, with_contamination, radius_ratio, excluded_stars, toi, klims, clear_field_only, check_saturation, contamination_model, contamination_reference_passband)
86 contamination_reference_passband: str = "r'"):
87
---> 88 super().__init__(target, date, tid, cids, dataroot=dataroot,
89 nlegendre=nlegendre, npop=npop, mjd_start=mjd_start, mjd_end=mjd_end,
90 excluded_mjd_ranges=excluded_mjd_ranges,
~/github/MuSCAT2_transit_pipeline/muscat2ta/transitanalysis.py in __init__(self, target, date, tid, cids, dataroot, exptime_min, nlegendre, npop, mjd_start, mjd_end, excluded_mjd_ranges, aperture_lims, passbands, use_opencl, with_transit, with_contamination, radius_ratio, klims, catalog_name, init_lpf, check_saturation, contamination_model, contamination_reference_passband, target_coordinates)
149
150 if self.init_lpf:
--> 151 self.lpf = M2LPF(target, self.phs, tid, cids, pbs, aperture_lims=aperture_lims, use_opencl=use_opencl,
152 with_transit=with_transit, with_contamination=with_contamination,
153 n_legendre=nlegendre, radius_ratio=radius_ratio, klims=klims,
~/github/MuSCAT2_transit_pipeline/muscat2ta/m2lpf.py in __init__(self, target, photometry, tid, cids, filters, aperture_lims, use_opencl, n_legendre, use_toi_info, with_transit, with_contamination, radius_ratio, noise_model, klims, contamination_model, contamination_reference_passband)
301 tm = QuadraticModel(interpolate=False, klims=klims, nk=1024, nz=1024)
302
--> 303 BaseLPF.__init__(self, target, filters, times, fluxes, wns, arange(len(photometry)), covariates,
304 arange(len(fluxes)), tm = tm, tref=floor(times[0].min()))
305 self._start_flares = len(self.ps)
~/miniconda3/envs/py3/lib/python3.8/site-packages/pytransit/lpf/lpf.py in __init__(self, name, passbands, times, fluxes, errors, pbids, covariates, wnids, tm, nsamples, exptimes, init_data, result_dir, tref, lnlikelihood)
208 # Inititalise the instrument
209 # --------------------------
--> 210 self._init_instrument()
211
212 self._init_lnlikelihood()
~/github/MuSCAT2_transit_pipeline/muscat2ta/m2lpf.py in _init_instrument(self)
425 filters = {'g': sdss_g, 'r': sdss_r, 'i':sdss_i, 'z_s':sdss_z}
426 self.instrument = Instrument('MuSCAT2', [filters[pb] for pb in self.passbands])
--> 427 self.cm = SMContamination(self.instrument, self.contamination_reference_passband)
428
429 def add_flare(self, loc, amp=(0, 0.2)):
~/miniconda3/envs/py3/lib/python3.8/site-packages/pytransit/contamination/contamination.py in __init__(self, instrument, ref_pb)
299 name of the reference passband
300 """
--> 301 super().__init__(instrument, ref_pb)
302
303 self._spectra: DataFrame = DataFrame(pd.read_hdf(resource_filename(__name__, join("data", "spectra.h5")), 'Z0'))
~/miniconda3/envs/py3/lib/python3.8/site-packages/pytransit/contamination/contamination.py in __init__(self, instrument, ref_pb)
85 """
86 self.instrument = instrument
---> 87 self._ri = self.instrument.pb_names.index(ref_pb)
88 self._rpb = ref_pb
89
ValueError: "r'" is not in list
Is it possible to use SG1 ephemeris instead of that from TOI catalog?
Recently, we have two cases (i.e. TOI 6077.01 and TOI 6254.01: http://research.iac.es/proyecto/muscat/stars/view/1563) were the ephemeris from TOI catalog is different/outdated compared to SG1. The wiki is able to fetch SG1 ephemeris automatically, so I wonder if this can be implemented in the pipeline as well.
I can clone but I can't install it
Basic multicolour transit modelling using a parametric systematics model.
Hi Hannu,
There were several cases when we cannot produce wcs files after running m2astrometry.
I asked Felipe before and he said in email, quote:
"when the astrometry fails we usually sets things by hand and I don't know any more details"
So could you instruct us how you do this manually so that at least we have an output for ta.create_example_frame()
for TFOP upload?
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.