hpparvi / muscat2_transit_pipeline Goto Github PK

View Code? Open in Web Editor NEW

4.0 4.0 0.0 10.62 MB

MuSCAT2 transit analysis pipeline

License: GNU General Public License v3.0

Jupyter Notebook 91.89% Python 8.11%

astronomy astrophysics exoplanet-transits exoplanets photometry

muscat2_transit_pipeline's People

Contributors

Stargazers

Watchers

muscat2_transit_pipeline's Issues

MuSCAT2 data IO

ta.apply_normalized_limits fails

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)

Need to decide the GP package

Choose the Gaussian process package

We have several options to choose from

Custom GP module to reduce dependencies
George (not sure if maintained any longer)
Scikit-learn
PyMC3

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.

# Plot raw light curves

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

ta.plot_basic_posteriors()

Hello, Hannu.
When I set radius_ratio='chromatic' in TransitAnalysis(), ta.plot_basic_posteriors() didin't work.
I attached the Error message.

fetching TOI table

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)

exocat query fix

MuSCAT2_transit_pipeline/muscat2ta/lpf.py

Line 71 in 9684de4

p = self.planet.P if self.planet else 5.

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?).

more informative raw_plot text labels

MuSCAT2_transit_pipeline/muscat2ph/phdata.py

Line 305 in e4cb993

def plot_raw(self, nstars, ylim=(0.85, 1.15), ncols=4, figsize=(11,11)):

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

jupyter kernel crashes after running ta.plot_possible_blends()

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

more informative prior description

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.

reference stars are not chose correctly

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?

ta.optimize says ValueError: could not broadcast input array from shape (200,0) into shape (200,239)

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)

Full central time of transit (in BJD) is not stored

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.

Additional option of analysis for MuSCAT1

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

with_contamination=True causes the error in ta.plot_final_fit()

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

visualizing appropriate photometric apertures

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?

Incorrect aperture size in `frozen_apertures`

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.

MuSCAT2_transit_pipeline/muscat2ta/tfopanalysis.py

Line 603 in 7fc2c73

tap, raps = self.lpf.frozen_apertures

ta.plot_fit() fails when single band is only available

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

missing modules

MuSCAT2_transit_pipeline/muscat2ta/lpf.py

Line 23 in 9684de4

 from pytransit.contamination import SMContamination, apparent_radius_ratio, true_radius_ratio 

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.

Simple aperture size optimisation

Preliminary dumb version of aperture size optimisation.

Parallelisation

Parallelisation scheme will depend on the chosen MCMC sampler. Can be either

OpenMP: using only PyTransit's OpenMP multi-threading inside a single computer
MPI: very easy if using emcee for MCMC, but useful only if we have a cluster to run the computations in
ipyparallel: possible, but needs to be studies. Also only useful if we have access to a cluster

Also, parallelisation (other than OpenMP) may not be necessary given the small amount of data per analysis.

implement TransitAnalysis.load() as inverse of TransitAnalysis.save()

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

error when r-band is missing

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

use SG1 catalog by default instead of TOI catalog

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.

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?