eclipseplotter.py

from __future__ import absolute_import
from __future__ import print_function
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from matplotlib.ticker import MaxNLocator
from trm import roche
import sys
import configobj
import lfit
import emcee
import warnings
import GaussianProcess as GP
from mcmc_utils import *
import seaborn
from collections import MutableSequence
from model import Model
import time

sys.settrace

# parallellise with MPIPool
from emcee.utils import MPIPool
from six.moves import range

class LCModel(Model):
    """CV lightcurve model for multiple eclipses.
    
       Can be passed to routines for calculating model, chisq, prior, prob, etc.
       Can add eclipses at will with addEcl function. All eclipses share q, dphi, rwd.
       All other parameters vary. You cannot mix and match complex and simple bright 
       spot models for eclipses; all must use the same type of bright spot."""
       
    def __init__(self,parList,complex,nel_disc=1000,nel_donor=400):
        """Initialise model.
        
        Parameter list should be a 14 element (non-complex BS) or 18 element (complex BS)
        dictionary of Param objects. These are:
        wdFlux, dFlux, sFlux, rsFlux, q, dphi, rdisc, ulimb, rwd, scale, az, fis, dexp, phi0
        And additional params: exp1, exp2, tilt, yaw"""
        
        # Use of the super function allows abstract class in model.py to be referenced
        # Here the initialise function is referenced
        super(LCModel,self).__init__(parList)
        self.complex = complex
        
        # Need a way of checking number of parameters is correct
        if complex:
            assert len(parList)==18, "Wrong number of parameters"
        else:
            assert len(parList)==14, "Wrong number of parameters"
        
        # We need an LFIT CV object to do the calculations
        # First we create list of parameter names (first eclipse = 0)
        # Then we get values for each parameter through using getValue function from model.py
        # Finally, CV object calculated from these values
        parNames = ['wdFlux_0', 'dFlux_0', 'sFlux_0', 'rsFlux_0', 'q', 'dphi',\
            'rdisc_0', 'ulimb_0', 'rwd', 'scale_0', 'az_0', 'fis_0', 'dexp_0', 'phi0_0']
        if complex:
            parNames.extend(['exp1_0', 'exp2_0', 'tilt_0', 'yaw_0'])
        parVals = [self.getValue(name) for name in parNames]
        self.cv = lfit.CV(parVals)
        
        # How many eclipses?
        self.necl = 1

    def addEclipse(self,parList):
        """Allows additional eclipses to be added.
        
        Parameter list should include 11 or 15 Param objects (all params individual
        to each eclipse), depending on complexity of the bright spot model. These should be:
        wdFlux, dFlux, sFlux, rsFlux, rdisc, ulimb, scale, az, fis, dexp, phi0
        and additional params: exp1, exp2, tilt, yaw""" 
        
        # Need a way of checking number of parameters is correct  
        if self.complex:
            assert len(parList) == 15, "Wrong number of parameters"
        else:
            assert len(parList) == 11, "Wrong number of parameters"
        
        # How many eclipses?
        self.necl += 1
        # Add params from additional eclipses to existing parameter list
        self.plist.extend(parList)
        
    def calc(self,ecl,phi,width=None):
        """Extracts current parameter values for each eclipse and calculates CV flux."""

        # Template required for parameter names
        parNameTemplate = ['wdFlux_{0}', 'dFlux_{0}', 'sFlux_{0}', 'rsFlux_{0}', 'q', 'dphi',\
            'rdisc_{0}', 'ulimb_{0}', 'rwd', 'scale_{0}', 'az_{0}', 'fis_{0}', 'dexp_{0}', 'phi0_{0}']
        if complex:
            parNameTemplate.extend(['exp1_{0}', 'exp2_{0}', 'tilt_{0}', 'yaw_{0}'])
        # Template needs updating depending on eclipse number
        parNames = [template.format(ecl) for template in parNameTemplate]
        # List filled in with current parameter values 
        parVals = [self.getValue(name) for name in parNames]
        # CV Flux calculated from list of current parameter values
        try:
            return self.cv.calcFlux(parVals,phi,width)
        except:
            return -np.inf
        
    def chisq(self,phi,y,e,width=None):
        """Calculates chisq, which is required in ln_like"""
        retVal = 0.0
        for iecl in range(self.necl):
            if width:
                thisWidth=width[iecl]
            else:
                thisWidth=None
            # chisq calculation
            resids = (y[iecl] - self.calc(iecl,phi[iecl],thisWidth)) / e[iecl]
            # Check for bugs in model
            if np.any(np.isinf(resids)) or np.any(np.isnan(resids)):
                warnings.warn('model gave nan or inf answers')
                return -np.inf
            retVal += np.sum(resids**2)
        return retVal
        
    def ln_prior(self,verbose=False):
        """Returns the natural log of the prior probability of this model.
        
        Certain parameters (dphi, rdisc, scale, az) need to be treated as special cases,
        as the model contains more prior information than included in the parameter priors"""
        
        # Use of the super function allows abstract class in model.py to be referenced
        # Here the ln_prior function is referenced
        retVal = super(LCModel,self).ln_prior()
        
        # Remaining part of this function deals with special cases
        # dphi
        tol = 1.0e-6
        try:
            # Uses getParam function from model.py to get the objects of variable parameters
            q = self.getParam('q')
            dphi = self.getParam('dphi')
            # maxphi is dphi when i = 90
            maxphi = roche.findphi(q.currVal,90.0)
            # dphi cannot be greater than (or within a certain tolerance of) maxphi
            if dphi.currVal > maxphi-tol:
                if verbose:
                    print('Combination of q and dphi is invalid')
                retVal += -np.inf
            
        except:
            # We get here when roche.findphi raises error - usually invalid q
            retVal += -np.inf
            if verbose:
                print('Combination of q and dphi is invalid')
        
        # rdisc 
        try:
            xl1 = roche.xl1(q.currVal) # xl1/a
            maxrdisc = 0.46/xl1 # Maximum size disc can reach before precessing
            # rdisc is unique to each eclipse, so have to use slightly different method to 
            # obtain its object, compared to q and dphi which are shared parameters
            rdiscTemplate = 'rdisc_{0}'
            for iecl in range(self.necl):
                rdisc = self.getParam(rdiscTemplate.format(iecl))
                # rdisc cannot be greater than maxrdisc
                if rdisc.currVal > maxrdisc:
                    retVal += -np.inf
                
        except:
            # We get here when roche.findphi raises error - usually invalid q
            if verbose:
                print('Rdisc and q imply disc does not hit stream')
            retVal += -np.inf
        
        #BS scale
        rwd = self.getParam('rwd')
        minscale = rwd.currVal/3 # Minimum BS scale equal to 1/3 of rwd
        maxscale = rwd.currVal*3 # Maximum BS scale equal to 3x rwd
        scaleTemplate = 'scale_{0}'
        for iecl in range(self.necl):
            scale = self.getParam(scaleTemplate.format(iecl))
            # BS scale must be within allowed range 
            if scale.currVal < minscale or scale.currVal > maxscale:
                retVal += -np.inf
                if verbose:
                    print('BS Scale is not between 1/3 and 3 times WD size')
            
        #BS az
        slope = 80.0
        try:
            # Find position of bright spot where it hits disc
            azTemplate = 'az_{0}'
            for iecl in range(self.necl):
                rdisc = self.getParam(rdiscTemplate.format(iecl))
                rd_a = rdisc.currVal*xl1 # rdisc/a
                az = self.getParam(azTemplate.format(iecl))
                # Does stream miss disc? (rdisc/a < 0.2 or rdisc/a > 0.65 )
                # If yes, Tom's code will fail
                # Calculate position of BS
                x,y,vx,vy = roche.bspot(q.currVal,rd_a)
                # Find tangent to disc at this point
                alpha = np.degrees(np.arctan2(y,x))
                # Alpha is between -90 and 90. If negative, spot lags disc (i.e. alpha > 90)
                if alpha < 0: alpha = 90 - alpha
                tangent = alpha + 90 # Disc tangent
                # Calculate minimum and maximum az values using tangent and slope
                minaz = max(0,tangent-slope)
                maxaz = min(178,tangent+slope)
                # BS az must be within allowed range
                if az.currVal < minaz or az.currVal > maxaz:
                    retVal += -np.inf

        except:
            if verbose:
                print('Stream does not hit disc, or az is outside 80 degree tolerance')
            # We get here when roche.findphi raises error - usually invalid q
            retVal += -np.inf
           
        if complex:
            # BS exponential parameters
            # Total extent of bright spot is scale*(e1/e2)**(1/e2)
            # Limit this to half an inner lagrangian distance
            scaleTemplate = 'scale_{0}'
            exp1Template = 'exp1_{0}'
            exp2Template = 'exp2_{0}'
            for iecl in range(self.necl):
                sc = self.getParam(scaleTemplate.format(iecl))
                e1 = self.getParam(exp1Template.format(iecl))
                e2 = self.getParam(exp2Template.format(iecl))
                if sc.currVal*(e1.currVal/e2.currVal)**(1.0/e2.currVal) > 0.5:
                    retVal += -np.inf
    
        return retVal
         
    def ln_like(self,phi,y,e,width=None):
        """Calculates the natural log of the likelihood"""
        return -0.5*self.chisq(phi,y,e,width)
        
    def ln_prob(self,phi,y,e,width=None):
        """Calculates the natural log of the posterior probability (ln_prior + ln_like)"""
        # First calculate ln_prior
        lnp = self.ln_prior()
        # Then add ln_prior to ln_like
        if np.isfinite(lnp):
            try:
                return lnp + self.ln_like(phi,y,e,width)
            except:
                return -np.inf
        else:
            return lnp
            
class GPLCModel(LCModel):
    """CV lightcurve model for multiple eclipses, with added Gaussian process fitting"""
    
    def __init__(self,parList,complex,ampin_gp,ampout_gp,tau_gp,nel_disc=1000,nel_donor=400):
        """Initialise model.
        
        Parameter list should be a 17 element (non-complex BS) or 21 element (complex BS)
        dictionary of Param objects. These are:
        ampin_gp, ampout_gp, tau_gp, wdFlux, dFlux, sFlux, rsFlux, q, dphi, rdisc, ulimb, rwd, scale, az, fis, dexp, phi0
        And additional params: exp1, exp2, tilt, yaw"""
        
        super(GPLCModel,self).__init__(parList,complex,nel_disc,nel_donor)
        # Make sure GP parameters are variable when using this model
        self.plist.append(ampin_gp)
        self.plist.append(ampout_gp)
        self.plist.append(tau_gp)
        self._dist_cp = 10.0
        self._oldq = 10.0
        self._olddphi = 10.0
        self._oldrwd = 10.0
        
    def calcChangepoints(self,phi):
            
        # Also get object for dphi, q and rwd as this is required to determine changepoints
        dphi = self.getParam('dphi')
        q = self.getParam('q')
        rwd = self.getParam('rwd')
        phi0Template = 'phi0_{0}'
    
        dphi_change = np.fabs(self._olddphi - dphi.currVal)/dphi.currVal
        q_change    = np.fabs(self._oldq - q.currVal)/q.currVal
        rwd_change  = np.fabs(self._oldrwd - rwd.currVal)/rwd.currVal
    
        if (dphi_change > 1.2) or (q_change > 1.2) or (rwd_change > 1.2):
            # Calculate inclination
            inc = roche.findi(q.currVal,dphi.currVal)
            # Calculate wd contact phases 3 and 4
            phi3, phi4 = roche.wdphases(q.currVal, inc, rwd.currVal, ntheta=10)
            # Calculate length of wd egress
            dpwd = phi4 - phi3
            # Distance from changepoints to mideclipse
            dist_cp = (dphi.currVal+dpwd)/2.
        else:
            dist_cp = self._dist_cp
                        
        '''# Find location of all changepoints
        for iecl in range(self.necl):
            changepoints = []
            phi0 = self.getParam(phi0Template.format(iecl))
            # the following range construction gives a list
            # of all mid-eclipse phases within phi array
            for n in range (int( phi.min() ), int( phi.max() )+1, 1):
                changepoints.append(n+phi0.currVal-dist_cp)
                changepoints.append(n+phi0.currVal+dist_cp) '''
    
        # Find location of all changepoints
       
        changepoints = []
        # the following range construction gives a list
        # of all mid-eclipse phases within phi array
        for n in range (int( phi.min() ), int( phi.max() )+1, 1):
            changepoints.append(n-dist_cp)
            changepoints.append(n+dist_cp)       
        
        # save these values for speed
        if (dphi_change > 1.2) or (q_change > 1.2) or (rwd_change > 1.2):
            self._dist_cp = dist_cp
            self._oldq = q.currVal
            self._olddphi = dphi.currVal
            self._oldrwd = rwd.currVal
        
        return changepoints
     
    def createGP(self,phi):
        """Constructs a kernel, which is used to create Gaussian processes.
        
        Using values for the two hyperparameters (amp,tau), amp_ratio and dphi, this function:
        creates kernels for both inside and out of eclipse, works out the location of any 
        changepoints present, constructs a single (mixed) kernel and uses this kernel to create GPs"""
    
        # Get objects for ampin_gp, ampout_gp, tau_gp and find the exponential of their current values
        ln_ampin = self.getParam('ampin_gp')
        ln_ampout = self.getParam('ampout_gp')
        ln_tau = self.getParam('tau_gp')
        
        ampin = np.exp(ln_ampin.currVal)
        ampout = np.exp(ln_ampout.currVal)
        tau = np.exp(ln_tau.currVal)
       
        # Calculate kernels for both out of and in eclipse WD eclipse
        # Kernel inside of WD has smaller amplitude than that of outside eclipse,
        k_in  = ampin*GP.Matern32Kernel(tau)
        k_out = ampout*GP.Matern32Kernel(tau)
        #k_in  = ampin*GP.ExpKernel(tau)
        #k_out = ampout*GP.ExpKernel(tau)
        
        changepoints = self.calcChangepoints(phi)
        
        # Depending on number of changepoints, create kernel structure
        kernel_struc = [k_out]      
        for k in range (int( phi.min() ), int( phi.max() )+1, 1):
            kernel_struc.append(k_in)
            kernel_struc.append(k_out)
        
        # Create kernel with changepoints 
        kernel = GP.DrasticChangepointKernel(kernel_struc,changepoints)
        
        '''k1 = GP.Matern32Kernel(tau)
        
        gp_pars = np.array([ampout,ampin,ampout])
        changepoints = self.calcChangepoints(phi)
        
        k2 = GP.OutputScaleChangePointKernel(gp_pars,changepoints)
        
        kernel = k1*k2'''
        
        # Create GPs using this kernel
        gp = GP.GaussianProcess(kernel)
        return gp
            
    def ln_like(self,phi,y,e,width=None):
        """Calculates the natural log of the likelihood.
        
        This alternative ln_like function uses the createGP function to create Gaussian
        processes"""
        lnlike = 0.0
        # For each eclipse, create (and compute) Gaussian process and calculate the model
        for iecl in range(self.necl):
            gp = self.createGP(phi[iecl])
            gp.compute(phi[iecl],e[iecl])
            if width:
                thisWidth=width[iecl]
            else:
                thisWidth=None
            resids = y[iecl] - self.calc(iecl,phi[iecl],thisWidth)
                                
            # Check for bugs in model
            if np.any(np.isinf(resids)) or np.any(np.isnan(resids)):
                warnings.warn('model gave nan or inf answers')
                return -np.inf
                                
            # Calculate ln_like using lnlikelihood function from GaussianProcess.py             
            lnlike += gp.lnlikelihood(resids)         
        return lnlike
                
'''def parseInput(file):
    """Splits input file up making it easier to read"""
    # Reads in input file and splits it into lines
    blob = np.loadtxt(file,dtype='string',delimiter='\n')
    input_dict = {}
    for line in blob: 
        # Each line is then split at the equals sign
        k,v = line.split('=')
        input_dict[k.strip()] = v.strip()
    return input_dict'''
    
def parseInput(file):
    """Splits input file up making it easier to read"""
    input_dict = configobj.ConfigObj(file)
    return input_dict
            
if __name__ == "__main__":

    # Allows input file to be passed to code from argument line
    import argparse
    parser = argparse.ArgumentParser(description='Plot CV lightcurves')
    parser.add_argument('file',action='store',help='input file')
    args = parser.parse_args()
    
    # Use parseInput function to read data from input file
    input_dict = parseInput(args.file) 
    
    # Load in chain file
    file = input_dict['chain']
    
    # Read information about neclipses, plot ranges, complex bs, gps
    flat      = int( input_dict['flat'] )
    thin     = int( input_dict['nthin'] )
    neclipses = int(input_dict['neclipses'])
    start = float(input_dict['phi_start'])
    end = float(input_dict['phi_end'])
    complex   = bool(int(input_dict['complex']))
    useGP     = bool(int(input_dict['useGP']))
    corner = bool(int(input_dict['corner'])) 
    lc = bool(int(input_dict['lc'])) 
    
    if flat:
        chain = readflatchain(file)
    else:
        chain = readchain_dask(file)
        nwalkers, nsteps, npars = chain.shape
        chain = flatchain(chain,npars,thin=thin)
    
    # Read in file names containing eclipse data, as well as output plot names
    files = []
    output_plots = []
    for ecl in range(0,neclipses):
        files.append(input_dict['file_{0}'.format(ecl)])
        output_plots.append(input_dict['plot_{0}'.format(ecl)])
        
    # Create a model from the first eclipses (eclipse 0) parameters
    parNames = ['wdFlux_0', 'dFlux_0', 'sFlux_0', 'rsFlux_0', 'q', 'dphi',\
        'rdisc_0', 'ulimb_0', 'rwd', 'scale_0', 'az_0', 'fis_0', 'dexp_0', 'phi0_0']
    if complex:
        parNames.extend(['exp1_0', 'exp2_0', 'tilt_0', 'yaw_0'])
    # List of values obtained from input file using fromString function from mcmc_utils.py
    parList = [Param.fromString(name, input_dict[name]) for name in parNames]
    
    # Read in GP params using fromString function from mcmc_utils.py
    ampin_gp = Param.fromString('ampin_gp', input_dict['ampin_gp'])
    ampout_gp = Param.fromString('ampout_gp', input_dict['ampout_gp'])
    tau_gp = Param.fromString('tau_gp', input_dict['tau_gp'])
        
    # If fitting using GPs use GPLCModel, else use LCModel
    if useGP:
        model = GPLCModel(parList,complex,ampin_gp,ampout_gp,tau_gp)
    else:
        model = LCModel(parList,complex)
        
    # pickle is used for parallelisation
    # pickle cannot pickle methods of classes, so we wrap
    # the ln_prior, ln_like and ln_prob functions here to 
    # make something that can be pickled
    def ln_prior(parList):
        model.pars = parList
        return model.ln_prior()
    def ln_like(parList,phi,y,e,width=None):
        model.pars = parList
        #return model.ln_like(phi,y,e,width)
        lnlike = model.ln_like(phi,y,e,width)
        if np.isfinite(lnlike):
        	return lnlike
        else:
        	print(parList)
        	print(lnlike)
        	return lnlike
    def ln_prob(parList,phi,y,e,width=None):
        model.pars = parList
        return model.ln_prob(phi,y,e,width)
        
    # Add in additional eclipses as necessary
    parNameTemplate = ['wdFlux_{0}', 'dFlux_{0}', 'sFlux_{0}', 'rsFlux_{0}',\
        'rdisc_{0}', 'ulimb_{0}', 'scale_{0}', 'az_{0}', 'fis_{0}', 'dexp_{0}', 'phi0_{0}']
    if complex:
        parNameTemplate.extend(['exp1_{0}', 'exp2_{0}', 'tilt_{0}', 'yaw_{0}']) 
    for ecl in range(1,neclipses):
        # This line changes the eclipse number for each parameter name
        parNames = [template.format(ecl) for template in parNameTemplate]
        # List of values obtained from input file using fromString function from mcmc_utils.py
        parList = [Param.fromString(name, input_dict[name]) for name in parNames]
        # Use addEclipse function defined above to add eclipse parameters to parameter list
        model.addEclipse(parList)
        
    # Store your data in python lists, so that x[0] are the times for eclipse 0, etc.
    x = []
    y = []
    e = []
    w = []
    
    # Read in eclipse data
    for file in files:
        xt,yt,et = np.loadtxt(file,skiprows=16).T
        wt = np.mean(np.diff(xt))*np.ones_like(xt)/2.
        # Create mask
        mask = (xt > start) & (xt < end)
        x.append(xt[mask])
        y.append(yt[mask])
        e.append(et[mask])
        w.append(wt[mask])
       
    # How many parameters?  
    npars = model.npars
    # Current values of all parameters
    params = [par.currVal for par in model.pars] 
    
    # Print out and save individual parameters to file
    f = open("modparams.txt","w")
    f.close()
    params = []
    for i in range(npars):
        par = chain[:,i]
        lolim,best,uplim = np.percentile(par,[16,50,84])
        print("%s = %f +%f -%f" % (model.lookuptable[i],best,uplim-best,best-lolim))
        f = open("modparams.txt","a")
        f.write("%s = %f +%f -%f\n" % (model.lookuptable[i],best,uplim-best,best-lolim))
        f.close()
        params.append(best)
    # Plot cornerplot (include only shared and 1st eclipse parameters)
    if corner:
        # First extract shared and 1st eclipse parameters from chain
        if complex:
            pars_cp = 18
        else:
            pars_cp = 14
        if useGP:
            pars_cp += 3
        chain_cp = chain[:,0:pars_cp] 
        #chain_cp = chain[:,pars_cp:28]
        # Create corner plot
        fig = thumbPlot(chain_cp,model.lookuptable[0:pars_cp])
        #fig = thumbPlot(chain_cp,model.lookuptable[pars_cp:28])
        fig.savefig('cornerPlot.pdf')
        plt.close()
    # Update model with best parameters
    model.pars = params
    
    # Only plot lightcurves if instructed to do so
    if lc == 0:
    	sys.exit()
    
    # Print out chisq, ln prior, ln likelihood and ln probability for the model                
    print('\nFor this model:\n')
    # Size of data required in order to calculate degrees of freedom (D.O.F)
    dataSize = np.sum((xa.size for xa in x))
    print("Chisq          = %.2f (%d D.O.F)" % (model.chisq(x,y,e,w),dataSize - model.npars - 1))
    print("ln prior       = %.2f" % model.ln_prior())
    print("ln likelihood = %.2f" % model.ln_like(x,y,e,w))
    print("ln probability = %.2f" % model.ln_prob(x,y,e,w))
    print('\nFrom wrapper functions:\n')
    print("ln prior       = %.2f" % ln_prior(params))
    print("ln likelihood = %.2f" % ln_like(params,x,y,e,w))
    print("ln probability = %.2f" % ln_prob(params,x,y,e,w))
    
    # Save these to file
   
    f = open("modparams.txt","w")
    f.write("\nFor this model:\n\n")
    f.write("Chisq          = %.2f (%d D.O.F)\n" % (model.chisq(x,y,e,w),dataSize - model.npars - 1))
    f.write("ln prior       = %.2f\n" % model.ln_prior())
    f.write("ln likelihood = %.2f\n" % model.ln_like(x,y,e,w))
    f.write("ln probability = %.2f\n" % model.ln_prob(x,y,e,w))
    f.close()
    
    # Plot model & data
    # Use of gridspec to help with plotting
    gs = gridspec.GridSpec(2,1,height_ratios=[2,1])
    gs.update(hspace=0.0)
    seaborn.set(style='ticks')
    seaborn.set_style({"xtick.direction": "in","ytick.direction": "in"})

    for iecl in range(neclipses):
        xp = x[iecl]
        yp = y[iecl]
        ep = e[iecl]
        wp = w[iecl]
           
        xf = np.linspace(xp.min(),xp.max(),1000)
        wf = 0.5*np.mean(np.diff(xf))*np.ones_like(xf)
        yp_fit = model.calc(iecl,xp,wp)
        yf = model.calc(iecl,xf,wf)
        res = yp - yp_fit
        
        # Needed for plotting GP 
        if useGP:
            gp = model.createGP(xp)
            gp.compute(xp,ep)
            samples = gp.sample_conditional(res, xp, size = 300)
            mu = np.mean(samples,axis=0)
            std = np.std(samples,axis=0)
            fmu, _ = gp.predict(res, xf)
        
        ax1 = plt.subplot(gs[0,0])
        
        # CV model
        ax1.plot(xf,yf)
        ax1.plot(xf,model.cv.yrs)
        ax1.plot(xf,model.cv.ys)
        ax1.plot(xf,model.cv.ywd)
        ax1.plot(xf,model.cv.yd)
                      
        # Plot fill-between region
        #Read chain
        #First few lines used to determine number of parameters
        if complex:
            pars_fill = 15
        else:
            pars_fill = 11
        pars_fill_e1 = pars_fill + 3
        if useGP:
            pars_fill_e1 += 3
        if iecl == 0 and useGP == 1:
            # Read parameters of first eclipse from chain
            pars = chain[:,0:pars_fill_e1-3]
        if iecl == 0 and useGP == 0:
        	# Read parameters of first eclipse from chain
            pars = chain[:,0:pars_fill_e1]
        if iecl >> 0:
            # Parameters from additional eclipses are a little trickier
            # to read from chain
            pars_1 = chain[:,pars_fill_e1+pars_fill*(iecl-1):pars_fill_e1+pars_fill*(iecl-1)+4]
            pars_2 = chain[:,4:6]
            pars_3 = chain[:,pars_fill_e1+pars_fill*(iecl-1)+4:pars_fill_e1+pars_fill*(iecl-1)+6]
            pars_4 = chain[:,8]
            pars_5 = chain[:,pars_fill_e1+pars_fill*(iecl-1)+6:pars_fill_e1+pars_fill*iecl]
            pars = []
            for i in range(0,len(pars_1)):
                new_pars = np.append(pars_1[i],pars_2[i])
                new_pars = np.append(new_pars,pars_3[i])
                new_pars = np.append(new_pars,pars_4[i])
                new_pars = np.append(new_pars,pars_5[i])
                pars.append(new_pars)
            pars = np.array(pars)
                    
        # Create array of 1000 random numbers
        random_sample = np.random.randint(0,len(pars[0:]),1000)
                
        lcs = []
        for i in random_sample:
            CV = lfit.CV(pars[i])
                
            xf_2 = np.linspace(xp.min(),xp.max(),1000)
            wf_2 = 0.5*np.mean(np.diff(xf_2))*np.ones_like(xf_2)
            yf_2 = CV.calcFlux(pars[i],xf_2,wf_2)
            lcs.append(yf_2)
              
        # Plot filled area    
        lcs = np.array(lcs)
        mu_2 = lcs.mean(axis=0)
        std_2 = lcs.std(axis=0)
        ax1.fill_between(xf_2,mu_2+std_2,mu_2-std_2,color='b',alpha=0.2)
           
        # Data
        
        if useGP:
        	ax1.errorbar(xp,yp,yerr=ep,fmt='.',color='k',capsize=0,alpha=0.2,markersize=5,linewidth=1)
        	fmu, _ = gp.predict(res, xp)
        	ax1.errorbar(xp,yp-fmu,yerr=ep,fmt='.',color='k',capsize=0,alpha=0.6,markersize=5,linewidth=1)
        else:
        	ax1.errorbar(xp,yp,yerr=ep,fmt='.',color='k',capsize=0,alpha=0.6,markersize=5,linewidth=1)
        ax2 = plt.subplot(gs[1,0],sharex=ax1)
        ax2.errorbar(xp,yp-yp_fit,yerr=ep,color='k',fmt='.',capsize=0,alpha=0.6,markersize=5,linewidth=1)
        ax1.set_ylim(ymin=-0.0001,ymax=0.3601)
        #ax1.set_ylim(ymin=-0.0001)
        ax1.set_xlim(xmin=start,xmax=end)
        ax1.tick_params(top='on',right='on')
        ax2.tick_params(top='on',right='on')
        #ax2.set_xlim(ax1.get_xlim())
        #ax2.set_xlim(-0.1,0.12)
        if useGP:
            ax2.fill_between(xp,mu+1.0*std,mu-1.0*std,color='r',alpha=0.4)

        # Labels
        ax1.set_ylabel('Flux (mJy)', fontsize=15)
        ax2.set_ylabel('Residuals (mJy)', fontsize=15)
        ax1.set_xlabel('Orbital Phase', fontsize=15)
        ax2.set_xlabel('Orbital Phase', fontsize=15)
        ax2.yaxis.set_major_locator(MaxNLocator(4,prune='both'))
        ax1.tick_params(axis='both',labelbottom='off',labelsize=14,top='on',right='on')
        ax2.tick_params(axis='both',labelsize=14,top='on',right='on')
        
        for ax in plt.gcf().get_axes()[::2]:
            ax.yaxis.set_major_locator(MaxNLocator(prune='both'))
            
        plt.subplots_adjust(bottom=0.095, top=0.965, left=0.12, right=0.975)
        
        # Save plot images
        plt.savefig(output_plots[iecl])
        plt.close()