PeriodicDetrend.py

import ipywidgets as widgets
from ipywidgets.widgets.widget_box import HBox
from ipywidgets.widgets.widget_description import DescriptionStyle
from ipywidgets.widgets.widget_string import Label
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import display
import os
from multiprocessing import Pool
from functools import partial
import matplotlib.patches as shp
from scipy import interpolate
from lmfit.models import VoigtModel, GaussianModel
import time

# Rodrigo Boufleur 2021 (c) 
# Função para calcular média móvel em dados não equiespaçados
# Retrona um com a media movel
# Não faz truncamento das bordas
def mmean(x, y, box):
    n = len(x)
    yout = np.zeros(n)
    for i in range(n):
        idx = (x >= (x[i] - box/2)) & (x < (x[i] + box/2))
        yout[i] = np.mean(y[idx])
    return yout




# Rodrigo Boufleur 2021 (c) 
# Busca período utilizando a técnica de plavchan

###############################################################################
# Description of the Plavchan algorithm taken from:
# https://exoplanetarchive.ipac.caltech.edu/docs/pgram/pgram_algo.html
#
# The Plavchan periodogram (Plavchan et al., 2008) is similar to a binless 
# variation of the "phase dispersion minimization" (PDM) algorithm 
# (Stellingwerf, 1978). In this method, the "basis" of periodic curves is 
# computed directly from the data. As in the BLS method, the time series is 
# folded to the candidate period. A dynamical prior is generated by box-car 
# smoothing the phased time series. The difference between the data and the 
# prior is squared and summed over a worst-fit subset of the data. When a 
# suitable period is found, the sum of the squared residuals from the smoothed 
# curve will be minimized. If no signal is present, the minimum sum of 
# squared errors will come from the model of no variability 
# (i.e., data values = constant). This is used as the normalization. 
# Periodogram power is defined as the normalization divided by the sum of 
# squared residuals to the smoothed curve. It will be greater than one 
# if the assumption of no variability is improved upon.
#
# * Number of outliers - The "number of outliers" parameter allows adjustment 
#   of the Plavchan power calculation. When comparing the time series to the 
#   dynamical prior, computation may be restricted to the N worst-fitting data 
#   points. The worst-fit data points may change for different candidate 
#   periods, as the prior also changes. This improves sensitivity in low 
#   signal-to-noise searches.
#
# * Phase-smoothing box size - The phase smoothing-box parameter specifies 
#   the width of the phase box over which to average the time-series data 
#   to compute the dynamical prior. A value of 0.05 is typical for ground-based 
#   transit surveys with a few thousand data points.


def _plavchan_statistics(test_period, x, y, linear_fit, n_outliers, box):
    '''Returns chi_ratio of period folded with linear model
    '''
    phase = (x/test_period) % 1
    idx = np.argsort(phase)
    average = mmean(phase[idx], y[idx], box)
    diff_model = ((average - y[idx])**2)
    idxpoor = np.argsort(diff_model)
    sum_model = np.sum( diff_model[idxpoor[-n_outliers:]] )
    return sum_model


def plavchan(x, y, per_ini, per_end, per_step=1/1400, window=0.06, n_outliers=500, n_threads=1, return_stats=False):
    '''Computes the plavchan statistics.
    '''
    
    n = len(x)   
    if n < n_outliers:
        n_outliers = n
        print('The number of outliers set exceeds the number of data points.')
        print('Using all data points as outliers.')
    
    period = np.arange(per_ini, per_end, per_step)
    # mean value linear fit
    linear_fit = np.repeat(np.mean(y), n)
    # run algorithm (_plavchan_statistics) in parallel
    pool = Pool(processes=n_threads)
    chipool=partial(_plavchan_statistics, x=x, y=y, linear_fit=linear_fit, n_outliers=n_outliers, box=window)
    statistics = np.asarray(pool.map(chipool, period))
    statistics = 1/statistics
    statistics = statistics/np.std(statistics, ddof=1)
    statistics = statistics-np.median(statistics)
    
    if return_stats:
        return [period, statistics]
    else:
        return period[np.argmax(statistics)]




class DetrendLightCurve():
    """Object containing the computational values
    """
    def __init__(self, time, flux, **kwargs):

        '''Constuctor method'''
        # allowed_kwargs = ['period']
        self.time = time
        self.flux = flux
        self.period = 0
        self.t = time
        self.f = flux
        self.solution = False
        self.detrend_log = []
        self.name = 'lightcurve'
        if 'period' in kwargs:
            self.period = kwargs['period']
        if 'name' in kwargs:
            self.name = kwargs['name']
        self.period_estimates = []
        self.period_err_estimates = []
        self.period_differences = []

    
    ### WIDGETS CREATION #######################################################
    def create_widgets(self):

        """generate widgets used in the app
        """
        ### EXECUTION SETTINGS WIDGETS
        # select the number or threads
        self.threads =  widgets.IntSlider(
            value=os.cpu_count()-1, 
            min=1, 
            max=os.cpu_count(), 
            description='THREADS: ',
            layout=widgets.Layout(width='30%')
        )

        # select fast or regular method (number of data points)
        self.method = widgets.RadioButtons(
            value='Regular',
            options=['Regular', 'Fast'], 
            description='METHOD: ',
            layout=widgets.Layout(width='30%')
        )

        # blank space
        self.blank_settings_tab_1 = widgets.Label(
            value='',
            layout=widgets.Layout(width='0%')
        )

        # create the 'set' button
        self.set_execution = widgets.Button(
            description="SET", 
            button_style='success',
            layout=widgets.Layout(width='10%')
        )
        # onclick behavion
        self.set_execution.on_click(self.on_set_execution_clicked)


        ### PERIOD SEARCH WIDGETS
        # define the minimum period to be searched
        self.minimum_period =  widgets.BoundedFloatText(
            value=0.5,
            min=0.0,
            description='MIN',
            step=0.5,
            layout=widgets.Layout(width='20%')

        )

        # define the maximum period to be searched
        self.maximum_period = widgets.BoundedFloatText(
            value=round((self.t.max()-self.t.min())/3.),
            min=0.0,
            description='MAX',
            step=0.5,
            layout=widgets.Layout(width='20%')

        )

        # define the period step of the search
        self.period_step = widgets.BoundedFloatText(
            value=0.05,
            min=0.0,
            step=0.01,
            description='STEP',
            layout=widgets.Layout(width='20%')

        )

        # define the window of the search
        self.window = widgets.BoundedFloatText(
            value=0.11,
            min=0.0,
            step=0.01,
            description='WINDOW',
            layout=widgets.Layout(width='20%')
        )

        # define the number of outliers
        self.outlier = widgets.BoundedIntText(
            value=int(len(self.t)*0.2),
            min=0,
            max=500,
            description='NUM OUTL',
            step=1,
            layout=widgets.Layout(width='20%')

        )

        self.blank_period_tab_1 = widgets.Label(
            value='',
            layout=widgets.Layout(width='50%')
        )

        # create the 'search' button
        self.period_lookup = widgets.Button(
            description="SEARCH", 
            button_style='success',
            layout = widgets.Layout(width='10%')
        )
        # onclick behavior
        self.period_lookup.on_click(self.on_period_lookup_clicked)
        



        ### CREATE DETREND WIDGETS

        # Mask 1
        self.mask1 = widgets.FloatRangeSlider(
            value=[0.0, 0.0],
            min=0,
            max=1,
            step=0.001,
            description='',
            disabled=False,
            continuous_update=True,
            orientation='horizontal',
            readout=False,
            readout_format='.2f',
            layout=widgets.Layout(width='auto', height='50%')
        )

        # Mask 2
        self.mask2 = widgets.FloatRangeSlider(
            value=[0.0, 0.0],
            min=0,
            max=1,
            step=0.001,
            description='',
            disabled=False,
            continuous_update=True,
            orientation='horizontal',
            readout=False,
            readout_format='.2f',
            layout=widgets.Layout(width='auto', height='50%')

        )

        # Mask 3
        self.mask3 = widgets.FloatRangeSlider(
            value=[0.0, 0.0],
            min=0,
            max=1,
            step=0.001,
            description='',
            disabled=False,
            continuous_update=True,
            orientation='horizontal',
            readout=False,
            readout_format='.2f',
            layout=widgets.Layout(width='auto', height='50%')

        )

        # range control widgets
        self.phasefolded_xaxis_slider = widgets.FloatRangeSlider(
            value=[0, 1],
            min=0,
            max=1,
            step=0.01,
            description='',
            disabled=False,
            continuous_update=True,
            orientation='horizontal',
            readout=False,
            readout_format='.2f',
            layout=widgets.Layout(width='auto')
        )

        self.phasefolded_yaxis_slider = widgets.FloatRangeSlider(
            value=[self.f.min()*0.95,self.f.max()*1.05],
            min=self.f.min()*0.95,
            max=self.f.max()*1.05,
            step=0.01,
            description='',
            disabled=False,
            continuous_update=True,
            orientation='vertical',
            readout=False,
            readout_format='.2f',
            layout=widgets.Layout(height='auto')
        )

        self.detrend_xaxis_slider = widgets.FloatRangeSlider(
            value=[self.t.min()-0.5, self.t.max()+0.5],
            min=self.t.min()-0.5,
            max=self.t.max()+0.5,
            step=0.1,
            description='',
            disabled=False,
            continuous_update=True,
            orientation='horizontal',
            readout=False,
            readout_format='.2f',
            layout=widgets.Layout(width='auto')
        )

        self.detrend_yaxis_slider = widgets.FloatRangeSlider(
            value=[self.f.min()*0.95,self.f.max()*1.05],
            min=self.f.min()*0.95,
            max=self.f.max()*1.05,
            step=0.01,
            description='',
            disabled=False,
            continuous_update=True,
            orientation='vertical',
            readout=False,
            readout_format='.2f',
            layout=widgets.Layout(height='auto')
        )
        
        # click buttons
        self.detrend_button = widgets.Button(
            description="DETREND", 
            button_style='success',
            layout=widgets.Layout(width='auto'),
            tooltip='Run the detrend algorithm.'
        )
        self.detrend_button.on_click(self.on_detrend_button_clicked)

        # dropdown lists
        self.interpolation_method_selection = widgets.Dropdown(
            options=[('Zero', 'zero'), ('Linear','slinear'), ('Quadratic','quadratic'), ('Cubic','cubic')],
            value='slinear',
            description='',
            layout=widgets.Layout(width='auto')
        )

        # define the window of the search
        self.moving_average_window = widgets.FloatText(
            value=0.33,
            min=0,
            step=0.01,
            layout=widgets.Layout(width='40%')        )
        
        self.reset_detrend = widgets.Button(
            description='RESET',
            button_style='danger',
            disabled=True,
            layout = widgets.Layout(width='auto'),
            tooltip='Reset and clear the log.'
        )
        self.reset_detrend.on_click(self.on_reset_detrend_clicked)        

 
        self.number_of_runs = widgets.IntText(
            value=1,
            description='NUM RUNS',
            disabled=False,
            layout=widgets.Layout(width='auto') 
        )

        self.optimize_period = widgets.Checkbox(
            value=True,
            description='Optimize Period'
        )
        
        self.save_results = widgets.Button(
            description='SAVE',
            button_style='info',
            disabled=True,
            layout = widgets.Layout(width='auto'),
            tooltip = 'Click to save the results.'
        )
        self.save_results.on_click(self.on_save_results_clicked)       

        # CREATE WIDGETS PLOT OUTPUTS
        self.plot_input_data = widgets.Output()
        self.plot_periodogram = widgets.Output()
        self.plot_folded_lightcurve = widgets.Output()
        self.plot_detrend = widgets.Output()
        self.detrend_log_out = widgets.Output()



    ### PLOTS GENERATION #######################################################
    # GENERATE TAB1 (SETTINGS) INPUT DATA PLOT
    def generate_plot_input_data(self):
        """Generates 'Execution Settings' tab plot
        """
        start, end = min(self.t)-0.5, max(self.t)+0.5
        xr = widgets.FloatRangeSlider(
            value=[start, end],
            min=start,
            max=end,
            step=1/1440,
            description=' ',
            disabled=False,
            continuous_update=True,
            orientation='horizontal',
            readout=False,
            readout_format='.2f',
            layout=widgets.Layout(width="90%")
        )

        plt.close()
        fig = plt.figure(figsize=(11, 4))
        ax = plt.axes()
        ax.grid(lw=.5)
        ax.set_xlabel('Time')
        ax.set_ylabel('Flux')
        fig.canvas.toolbar_visible  = False
        fig.canvas.header_visible=False
        fig.canvas.footer_visible=False

        @widgets.interact(xr=xr)
        def update(xr):
            xrl, xrr = xr
            [l.remove() for l in ax.lines]
            ax.plot(self.t, self.f, 'k.-')
            ax.set_xlim(xrl, xrr)

    # GENERATE TAB2 (PERIODOGRAM) PLOTS
    def generate_plot_plavchan_stats(self):
        """Plot periodogram from plavchan stats
        """
        start, end = min(self.plavchan_period), max(self.plavchan_period)
        xr = widgets.FloatRangeSlider(
            value=[start, end],
            min=start,
            max=end,
            step=1/1440,
            description=' ',
            disabled=False,
            continuous_update=True,
            orientation='horizontal',
            readout=False,
            readout_format='.2f',
            layout=widgets.Layout(width='90%')
        )

        logscale = widgets.Checkbox(
            value=False,
            description='Logaritmic Scale (X Axis)',
            continuous_update=True
        )

        plt.close()
        plavchan_plot = plt.figure(figsize=(11, 4))
        ax = plt.axes()
        ax.grid(lw=.5)
        ax.set_xlabel('Period')
        ax.set_ylabel('Plavchan Chi Ratio')
        plavchan_plot.canvas.toolbar_visible = False
        plavchan_plot.canvas.header_visible = False
        plavchan_plot.canvas.footer_visible = False

        @widgets.interact(xr=xr, logscale=logscale)
        def update(xr,logscale):
            xrl, xrr = xr
            [l.remove() for l in ax.lines]
            ax.plot(self.plavchan_period, self.plavchan_stats, 'r-')
            ax.set_xlim(xrl, xrr)
            if logscale:
                ax.set_xscale('log')
            else:
                ax.set_xscale('linear')

    # GENERATE TAB3 (FOLDED PHASE) PLOT
    def estimate_primary_transit_timing(self, phasevalue=False):
        '''Shift primary transit'''
        phase = (self.t/self.period) % 1
        idx = np.argsort(phase)
        mean_lc = mmean(phase[idx], self.f[idx], 0.01)
        reference = phase[idx[np.argmin(mean_lc)]]
        idx2 = (phase[idx] > (reference-0.1)) & (phase[idx] < (reference+0.1))

        mod = GaussianModel()
        x = phase[idx[idx2]]
        y = 1/self.f[idx[idx2]] - min(1/self.f[idx[idx2]])
        pars = mod.guess(y, x=x)
        out = mod.fit(y, pars, x=x)
        reference = out.params['center'].value
        if phasevalue:
            return reference
        if reference > phase[0]:
            return ( (reference - phase[0])*self.period + self.t[0] )
        if reference <= phase[0]:
            return ( (reference - phase[0] + 1)*self.period + self.t[0] )
        
    def shift_primary_transit_phase(self):
        '''Shift primary transit to 0.25'''
        defasage = 0.25 - self.estimate_primary_transit_timing(phasevalue=True)
        phase = (self.t/self.period) % 1
        phase2 = phase + defasage
        idx = phase2 > 1
        phase2[idx] = phase2[idx]-1
        idx = phase2 < 0
        phase2[idx] = phase2[idx]+1
        return phase2

    def generate_plot_folded_lightcurve(self):
        self.fig_phase_folded, self.ax_phase_folded = plt.subplots(constrained_layout=True, figsize=(5,3.71))
        self.ax_phase_folded = plt.axes()
        self.ax_phase_folded.grid(lw=.5)
        self.ax_phase_folded.set_xlabel('Orbital Phase')
        self.ax_phase_folded.set_ylabel('Flux')
        # self.fig_phase_folded.canvas.toolbar_position = 'bottom'
        self.fig_phase_folded.canvas.toolbar_visible  = False
        self.fig_phase_folded.canvas.header_visible = False
        self.fig_phase_folded.canvas.footer_visible = False

        if self.period > 0:
            self.phase = self.shift_primary_transit_phase()
            self.ax_phase_folded.plot(self.phase, self.f, '.', color='navy', alpha=0.75, label='Original')
        
        if self.solution:
            # plot detrended light curve
            self.phase = self.shift_primary_transit_phase()
            self.ax_phase_folded.plot(self.phase, self.f - self.trend , '.', color='orangered', label='Detrended' )
       
        # averaged folded light curve
        # self.ax_phase_folded.plot(self.phase[idx], mean_lc, 'k-', alpha=0.5)
       
        # Add making region 1
        self.m1 = plt.Rectangle(
            [self.mask1.value[0],self.f.min()*0.95],
            self.mask1.value[1]-self.mask1.value[0],
            self.f.max()*1.05-self.f.min()*0.95,
            edgecolor=None,
            facecolor='0.8')
        self.ax_phase_folded.add_artist(self.m1)
        # Add masking region 2
        self.m2 = plt.Rectangle(
            [self.mask2.value[0],self.f.min()*0.95],
            self.mask2.value[1]-self.mask2.value[0],
            self.f.max()*1.05-self.f.min()*0.95,
            edgecolor=None,
            facecolor='0.8')
        self.ax_phase_folded.add_artist(self.m2)
        # Add masking region 3
        self.m3 = plt.Rectangle(
            [self.mask3.value[0],self.f.min()*0.95],
            self.mask3.value[1]-self.mask3.value[0],
            self.f.max()*1.05-self.f.min()*0.95,
            edgecolor=None,
            facecolor='0.8')
        self.ax_phase_folded.add_artist(self.m3)
        # Plot legend
        if self.period > 0:
            plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.15), ncol=3, fancybox=False)
        # Monitor and observe updates
        self.phasefolded_xaxis_slider.observe(self.update_xaxisrange_folded_lightcurve, 'value')
        self.phasefolded_yaxis_slider.observe(self.update_yaxisrange_folded_lightcurve, 'value')
        self.mask1.observe(self.update_mask1, 'value')
        self.mask2.observe(self.update_mask2, 'value')
        self.mask3.observe(self.update_mask3, 'value') 
        
    # GENERATE TAB3 (DETREND MODEL) PLOT
    def generate_plot_detrend(self):
        self.fig_detrend, self.ax_detrend = plt.subplots(constrained_layout=True, figsize=(9.8,4))
        self.ax_detrend = plt.axes()
        self.ax_detrend.grid(lw=.5)
        self.ax_detrend.set_xlabel('Time')
        self.ax_detrend.set_ylabel('Flux')
        # self.fig_detrend.canvas.toolbar_position = 'bottom'
        self.fig_detrend.canvas.toolbar_visible  = False
        self.fig_detrend.canvas.header_visible = False
        self.fig_detrend.canvas.footer_visible = False
        # plot light curve
        self.ax_detrend.plot(self.t, self.f, '.-', color='grey', label='Original', alpha=0.5 )
        
        if self.solution:
                # plot detrended light curve
                self.ax_detrend.plot(self.t, self.f - self.trend , '-', color='navy', label='Detrended' )
                # plot trend
                self.ax_detrend.plot(self.t, self.trend + np.mean(self.f), '-', color='orangered', label='Trend Model' )
                
        # Monitor and observe updates
        self.detrend_xaxis_slider.observe(self.update_xaxisrange_detrend, 'value')
        self.detrend_yaxis_slider.observe(self.update_yaxisrange_detrend, 'value')
        plt.legend(loc='upper center', bbox_to_anchor=(0.5, 1.15), ncol=3, fancybox=False)
    
        
    ### UPDATE FUNCTIONS USED IN THE DETREND TAB ###############################
    def update_mask1(self, change):
        self.m1.set_xy([change.new[0],self.f.min()*0.95])
        self.m1.set_width(change.new[1]-change.new[0])

    def update_mask2(self, change):
        self.m2.set_xy([change.new[0],self.f.min()*0.95])
        self.m2.set_width(change.new[1]-change.new[0])

    def update_mask3(self, change):
        self.m3.set_xy([change.new[0],self.f.min()*0.95])
        self.m3.set_width(change.new[1]-change.new[0])

    def update_xaxisrange_folded_lightcurve(self, change):
        self.ax_phase_folded.set_xlim(change.new)

    def update_yaxisrange_folded_lightcurve(self, change):
        self.ax_phase_folded.set_ylim(change.new)

    def update_xaxisrange_detrend(self, change):
        self.ax_detrend.set_xlim(change.new)

    def update_yaxisrange_detrend(self, change):
        self.ax_detrend.set_ylim(change.new)
       
       


    ### THE DETRENDING FUNCTION ################################################
    def detrend(self):
        '''Detrends the light curve
        '''
        #subtract mean lc
        self.phase = (self.t/self.period) % 1
        index = np.argsort(self.phase)
        average_lightcurve = mmean(self.phase[index], self.f[index], 0.01)
        average_lightcurve = average_lightcurve[np.argsort(index)]
        self.trend = self.f - average_lightcurve

        # locate the index of the masked sections
        index_mask1 = (self.phase > self.mask1.value[0]) & (self.phase < self.mask1.value[1])
        index_mask2 = (self.phase > self.mask2.value[0]) & (self.phase < self.mask2.value[1])
        index_mask3 = (self.phase > self.mask3.value[0]) & (self.phase < self.mask3.value[1])

        self.mask = np.logical_or(index_mask1, index_mask2, index_mask3)

#         index = np.asarray(index)

        if np.sum([self.mask1.value,self.mask2.value,self.mask3.value]) > 0:
            self.optimize_period.value = False #d isables period optimization
                                                # need to be better implemented in the future 
            index = []
            i = 0
            while i < (len(self.t)-2):
                if self.mask[i]:
                    index_start = i
                    while self.mask[i]:
                        i = i + 1
                    index_end = i
                    if (index_start > 0) or (index_end < len(self.t)):
                        index.append([index_start,index_end])
                i = i + 1
            # transform list into numpy array
            
            for i in index:
            #print(i)
            ### find the center and the width for the ith segment
                center = self.t[i[0]] + (self.t[i[1]] - self.t[i[0]])/2 
                width = self.t[i[1]] - self.t[i[0]]

                ### set right and left limits for each segment
                ### the limits are 1.5 times the width
                left_lim = np.where(self.t > (center - 1.5*width))
                left_lim = left_lim[0][0]
                right_lim = np.where(self.t < (center + 1.5*width))
                right_lim = right_lim[0][-1]

                ### concatenate the correspondent time and flux extracted segments
                t_segment = np.concatenate([self.t[left_lim:i[0]],self.t[i[1]:right_lim]])
                f_segment = np.concatenate([self.trend[left_lim:i[0]],self.trend[i[1]:right_lim]])

                # interpolate using the selected algorithm
                # values are extrapolated in the edges when necessary
                f_int = interpolate.interp1d(
                    t_segment,
                    f_segment,
                    kind=self.interpolation_method_selection.value, 
                    fill_value="extrapolate"
                )

                new_segment = self.t[i[0]:i[1]]
                self.trend[i[0]:i[1]] = f_int(new_segment)


        # compute moving average of trend
        box_size = self.moving_average_window.value
        self.interpolated_trend = mmean(self.t, self.trend, box_size)
        self.trend = self.interpolated_trend

        if self.optimize_period.value: 
            # recompute period
            plavchan_period, plavchan_stats = plavchan(
                        self.t,
                        self.f-self.trend,
                        self.period-0.025,
                        self.period+0.025,
                        per_step=1/5760,
                        window=self.window.value,
                        n_outliers=self.outlier.value,
                        n_threads=self.threads.value,
                        return_stats=True
                    )


            mod = VoigtModel()
            x = plavchan_period
            y = plavchan_stats - min(plavchan_stats)
            pars = mod.guess(y, x=x)
            out = mod.fit(y, pars, x=x)
            oldperiod = self.period

#             self.period = out.params['center'].value
            self.period = plavchan_period[np.argmax(plavchan_stats)]
            self.period_err = out.params['fwhm'].value / 3.6013 # conversion of fwhm to sigma
            self.trend = self.interpolated_trend
            self.period_estimates.append(self.period)
            self.period_err_estimates.append(self.period_err)
            self.period_differences.append(self.period - out.params['center'].value)

#             ### determine new maskvalues
#             phase_correction = (self.t[0]/oldperiod)-(self.t[0]/self.period)

#             if np.sum([self.mask1.value,self.mask2.value,self.mask3.value]) > 0:
#                 # ordenate maskvalues and put zeros at the end if necessary
#                 maskvalues = np.concatenate([self.mask1.value,self.mask2.value,self.mask3.value])
#                 idx = np.argsort(maskvalues)
#                 maskvalues = maskvalues[idx]
#                 if np.sum(maskvalues == 0) > 1:
#                     maskvalues = np.roll(maskvalues,np.sum(maskvalues>0))
#                 print(maskvalues)
#                 # correct the new phase values
#                 NonZeroValues = (maskvalues > 0)
#                 maskvalues[NonZeroValues] = maskvalues[NonZeroValues] - phase_correction        
#                 # fix negative values
#                 check_negatives = (maskvalues<0)
#                 if np.sum(check_negatives) > 0:
#                     maskvalues[check_negatives] += 1
#                     if (np.sum(maskvalues == 0) == 2) and (np.sum(check_negatives) < 2):
#                         maskvalues[-1] = 1
#                 # fix values over 1
#                 check_over1 = (maskvalues>1)
#                 if np.sum(check_over1) > 0:
#                     maskvalues[check_over1] -= 1
#                     # case of 3
#                     if (np.sum(maskvalues == 0) == 2) and (np.sum(check_over1) < 2):
#                         maskvalues[-1] = 1

#                 # resort maskvalues
#                 idx = np.argsort(maskvalues)
#                 maskvalues = maskvalues[idx]
#                 if np.sum(maskvalues == 0) > 1:
#                     maskvalues = np.roll(maskvalues,np.sum(maskvalues>0))
#                 self.mask1.value = [maskvalues[0],maskvalues[1]]
#                 self.mask2.value = [maskvalues[2],maskvalues[3]]
#                 self.mask3.value = [maskvalues[4],maskvalues[4]]

        self.solution = True



    ### ON CLICK BUTTONS BEHAVIOR ##############################################
    def on_set_execution_clicked(self, change):
        self.set_execution.description='Setting...'
        self.set_execution.button_style='warning'

        if self.method.value=='Fast':
            idx = np.arange(len(self.time))
            new_idx = (idx % 5 == 0)
            self.t, self.f = self.time[new_idx], self.flux[new_idx]
        else:
            self.t, self.f = self.time, self.flux

        self.plot_input_data.clear_output()
        with self.plot_input_data:
            output = '{} method chosen. {} core(s) selected. \n{} data points in the light curve.'
            print(output.format(self.method.value, self.threads.value, len(self.t))) 
            self.generate_plot_input_data()
            self.set_execution.description='Set'
            self.set_execution.button_style='success'

        self.plot_folded_lightcurve.clear_output()
        with self.plot_folded_lightcurve:
            self.generate_plot_folded_lightcurve()

        # Replot detrended values
        self.plot_detrend.clear_output()        
        with self.plot_detrend:
            self.generate_plot_detrend()

    def on_period_lookup_clicked(self, change):
        # change button visual state
        self.period_lookup.button_style='warning'
        self.period_lookup.disabled=True
        self.period_lookup.description='Searching...'

        # lookup for the period
        self.plavchan_period, self.plavchan_stats = plavchan(
            self.t,
            self.f,
            self.minimum_period.value,
            self.maximum_period.value,
            per_step=self.period_step.value,
            window=self.window.value,
            n_outliers=self.outlier.value,
            n_threads=self.threads.value,
            return_stats=True
        )

        #clear previous outputs
        self.plot_periodogram.clear_output()

        with self.plot_periodogram:
            self.period = self.plavchan_period[np.argmax(self.plavchan_stats)]
            print('Best Period: {:.8f}'.format(self.period)) 
            self.generate_plot_plavchan_stats()

            #return 'lookup' button to original state
            self.period_lookup.button_style='success'
            self.period_lookup.disabled=False
            self.period_lookup.description='SEARCH'

        # update folded lightcurve
        self.plot_folded_lightcurve.clear_output()
        with self.plot_folded_lightcurve:
            self.generate_plot_folded_lightcurve()

    def on_reset_detrend_clicked(self, change):
        """Clears plavchan outputs
        """
        self.reset_detrend.button_style='warning'
        self.reset_detrend.disabled=True
        self.reset_detrend.description='Resetting...'
        # Replot folded light curve

        self.detrend_log_out.clear_output()

        self.plot_folded_lightcurve.clear_output()
        self.solution = False
        with self.plot_folded_lightcurve:
            self.generate_plot_folded_lightcurve()

        # Replot detrended values
        self.plot_detrend.clear_output()        
        with self.plot_detrend:
            self.generate_plot_detrend()

        self.number_of_runs.value=1

        self.reset_detrend.button_style='danger'
        self.reset_detrend.disabled=False
        self.reset_detrend.description='RESET'
        self.save_results.disabled=True
        self.mask1.value=[0,0]
        self.mask2.value=[0,0]
        self.mask3.value=[0,0]

    def on_save_results_clicked(self, change):
        try:
            os.remove(self.name+'_detrended.txt')
        except:
            pass
        
        with open(self.name+'_detrended.txt', "ab") as f:
            f.write(b"#Time Flux Flux-Trend Trend\n")
            np.savetxt(f, np.c_[self.time, self.flux, self.flux-self.trend, self.trend], fmt=('%.8f %.8f %.8f %.8f'))
            f.close()
        
        with self.detrend_log_out:
                self.detrend_log_out.append_stdout('Results saved in {}\n'.format(self.name+'_detrended.txt'))

    def on_detrend_button_clicked(self, change):
        # change button visual state
        self.detrend_button.button_style='warning'
        self.detrend_button.disabled=True

        # reset estimates
        self.period_estimates = []
        self.period_err_estimates = []
        self.period_differences = []
        
        if not self.optimize_period.value:
            self.number_of_runs.value = 1
            
        for nrun in range(self.number_of_runs.value):
            if self.number_of_runs.value > 1:
                self.detrend_button.description='Run {}...'.format(nrun+1)
            else:
                self.detrend_button.description='Detrending...'

            with self.detrend_log_out:
                self.detrend_log_out.append_stdout('Run {:03d} | '.format(nrun+1))            # lookup for the period
            self.detrend()
            # Update the detrend log
            if self.optimize_period.value:
                with self.detrend_log_out:
                    self.detrend_log_out.append_stdout('Period estimate: {:.8f} +-{:.8f}\n'.format(self.period,self.period_err))
                    if (self.number_of_runs.value - nrun < 2):
                        self.detrend_log_out.append_stdout('          Primary transit epoch estimate: {:.8f}\n'.format(self.estimate_primary_transit_timing()))
                        self.detrend_log_out.append_stdout('          Plot shifted by {:.5f} in phase.\n\n\n'.format(0.25-self.estimate_primary_transit_timing(phasevalue=True)))

            else:
                with self.detrend_log_out:
                    self.detrend_log_out.append_stdout('Detrend executed without period optimization\n\n\n')
            
        # Replot folded light curve
        self.plot_folded_lightcurve.clear_output()
        with self.plot_folded_lightcurve:
            self.generate_plot_folded_lightcurve()

        # Replot detrended values
        self.plot_detrend.clear_output()        
        with self.plot_detrend:
            self.generate_plot_detrend()

            
#         if self.number_of_runs.value > 1:
#             self.period = np.mean(self.period_estimates)
#             self.period_err = np.sqrt( np.mean(np.array(self.period_err_estimates)**2) + np.std(self.period_estimates, ddof=1)**2 )
            
#             if self.optimize_period.value:
#                 with self.detrend_log_out:
#                     self.detrend_log_out.append_stdout('Period estimate from {} runs: {:.8f} +-{:.8f}\n'.format(self.number_of_runs.value, self.period, self.period_err))
#                     self.detrend_log_out.append_stdout('Primary transit epoch estimate: {:.8f}\n\n\n'.format(self.estimate_primary_transit_timing()))
#             else:
#                 with self.detrend_log_out:
#                     self.detrend_log_out.append_stdout('Detrend executed without period optimization\n\n\n')
                    
        self.detrend_button.button_style='success'
        self.detrend_button.disabled=False
        self.detrend_button.description='DETREND'
        self.reset_detrend.disabled=False
        self.save_results.disabled=False


    def display(self):
        """Diplay the application
        """
        self.create_widgets()
        
        with self.plot_input_data:
            self.generate_plot_input_data()
            
        self.setting_display = widgets.VBox( [ widgets.HBox( [self.threads, 
                                                              self.method,
                                                              self.blank_settings_tab_1, 
                                                              self.set_execution ] ),
                                               widgets.HBox( [self.plot_input_data ] )
                                             ], width='1080px', height='auto')
        
        self.period_search_display =  widgets.VBox( [ widgets.HBox( [ self.minimum_period, 
                                                                      self.maximum_period, 
                                                                      self.period_step, 
                                                                      self.window ] ),
                                                      widgets.HBox( [ self.outlier,
                                                                      self.blank_period_tab_1, 
                                                                      self.period_lookup ] ),
                                                      widgets.HBox( [ self.plot_periodogram ] )
                                                    ] )
        
        
        # function to make labels 
        
        def make_label(value, width='auto'):
            ### Make label widget   
            label = widgets.Label(
                value=value, 
                layout=widgets.Layout(width=width)
            )
            return label

        with self.plot_folded_lightcurve:
            self.generate_plot_folded_lightcurve()

        with self.plot_detrend:
            self.generate_plot_detrend()
            
        self.detrend_display = widgets.GridspecLayout(25, 30, width='1080px', height='auto')
        # phase folded plot interaction
        self.detrend_display[0:9,0:1] = self.phasefolded_yaxis_slider
        self.detrend_display[0:9,1:20] = self.plot_folded_lightcurve
        self.detrend_display[9:10,1:20] = self.phasefolded_xaxis_slider
        # mask interaction widgets
        self.detrend_display[1:2,21:30] = make_label('MASK 1:')
        self.detrend_display[2:3,21:30] = self.mask1
        self.detrend_display[3:4,21:30] = make_label('MASK 2:')
        self.detrend_display[4:5,21:30] = self.mask2
        self.detrend_display[5:6,21:30] = make_label('MASK 3:')
        self.detrend_display[6:7,21:30] = self.mask3
        # interpolation degree selction
        self.detrend_display[7:8,21:30] = make_label('POLYNOMIAL DEGREE INTERPOLATION (only for masked regions):')
        self.detrend_display[8:9,25:30] = self.interpolation_method_selection

        # moving average windows size
        self.detrend_display[11:12,17:22] = make_label('MOV. AVG. WINDOW')
        self.detrend_display[11:12,22:23] = self.moving_average_window
        # event button
        self.detrend_display[11:12,23:27] = self.detrend_button
        self.detrend_display[11:12,27:30] = self.reset_detrend
        
        self.detrend_display[12:13,17:23] = self.optimize_period
        self.detrend_display[12:13,23:30] = self.number_of_runs  
        
        self.detrend_display[13:14,27:30] = self.save_results

        # detrend result plot
        self.detrend_display[15:24,0:1] = self.detrend_yaxis_slider
        self.detrend_display[15:24,1:30] = self.plot_detrend
        self.detrend_display[24:25,1:30] = self.detrend_xaxis_slider

        
        ### GENERATE THE OUTPUT TABS
        self.tab_contents = ['Execution Settings', 'Period Search', 'Detrend', 'Detrend Log']
        self.children = [ self.setting_display,
                     self.period_search_display,
                     self.detrend_display,
                     self.detrend_log_out]
        self.tab = widgets.Tab()
        self.tab.children = self.children
        self.tab.set_title(0, self.tab_contents[0])
        self.tab.set_title(1, self.tab_contents[1])
        self.tab.set_title(2, self.tab_contents[2])
        self.tab.set_title(3, self.tab_contents[3])
        display(self.tab)