do_rates_mle_wPSs.py

import numpy as np
import pandas as pd
from scipy import optimize, stats, interpolate
from astropy.io import fits
from astropy.wcs import WCS
import os
import argparse
import logging, traceback, time

from config import quad_dicts, EBINS0, EBINS1, drm_quad_dir, solid_angle_dpi_fname, fp_dir
from sqlite_funcs import get_conn, append_rate_tab
from dbread_funcs import get_info_tab, guess_dbfname, get_files_tab,\
                            get_twinds_tab, get_rate_fits_tab
from bkg_rate_estimation import get_quad_rate_objs_from_db, rate_obj_from_sqltab
from mle_rates_for_realtime import do_rate_mle, do_rate_mle_mp, get_abs_cor_rates,\
                            get_cnts_intp_obj, get_quad_cnts_tbins
from counting_and_quad_funcs import get_quad_cnts_tbins_fast, dmask2ndets_perquad
from drm_funcs import get_ebin_ind_edges, DRMs
from wcs_funcs import world2val
from event2dpi_funcs import det2dpis, mask_detxy
from models import Bkg_Model_wSA, Bkg_Model_wFlatA, CompoundModel,\
                    Point_Source_Model_Binned_Rates, Model
from ray_trace_funcs import RayTraces, FootPrints



def cli():
    parser = argparse.ArgumentParser()
    parser.add_argument('--evfname', type=str,\
            help="Event data file",
            default=None)
    parser.add_argument('--Njobs', type=int,\
            help="Total number of jobs",
            default=16)
    parser.add_argument('--job_id', type=int,\
            help="Which job this is",
            default=-1)
    parser.add_argument('--pix_fname', type=str,\
            help="Name of the file with good imx/y coordinates",\
            default='good_pix2scan.npy')
    parser.add_argument('--bkg_fname', type=str,\
            help="Name of the file with the bkg fits",\
            default='bkg_estimation.csv')
    parser.add_argument('--dbfname', type=str,\
            help="Name to save the database to",\
            default=None)
    parser.add_argument('--pcfname', type=str,\
            help="Name of the partial coding image",\
            default='pc_2.img')
    parser.add_argument('--min_pc', type=float,\
            help="Min partical coding fraction to use",\
            default=0.1)
    args = parser.parse_args()
    return args


def parse_bkg_csv(bkg_fname, solid_angle_dpi, ebins0, ebins1, bl_dmask, rt_dir):

    bkg_df = pd.read_csv(bkg_fname)
    col_names = bkg_df.columns
    nebins = len(ebins0)

    PSnames = []
    for name in col_names:
        if '_imx' in name:
            PSnames.append(name.split('_')[0])
    print PSnames
    Nsrcs = len(PSnames)
    if Nsrcs > 0:
        bkg_name = 'Background_'
    else:
        bkg_name = ''


    bkg_mod = Bkg_Model_wFlatA(bl_dmask, solid_angle_dpi, nebins, use_deriv=True)

    ps_mods = []

    if Nsrcs > 0:
        rt_obj = RayTraces(rt_dir)
        for i in range(Nsrcs):
            name = PSnames[i]
            imx = bkg_df[name+'_imx'][0]
            imy = bkg_df[name+'_imy'][0]
            mod = Point_Source_Model_Binned_Rates(imx, imy, 0.1,\
                                                  [ebins0,ebins1], rt_obj, bl_dmask,\
                                                  use_deriv=True, name=name)
            ps_mods.append(mod)

    return bkg_df, bkg_name, PSnames, bkg_mod, ps_mods


def im_dist(imx0, imy0, imx1, imy1):
    return np.hypot((imx1 - imx0), (imy1 - imy0))


def get_fp_arr(fp_dir):
    fnames = np.array(os.listdir(fp_dir))
    imxs = np.array([float(fn.split('_')[1]) for fn in fnames])
    imys = np.array([float(fn.split('_')[2][:-4]) for fn in fnames])
    dtp = [('imx', np.float), ('imy', np.float),
           ('fname', fnames.dtype)]
    fp_arr = np.empty(len(imxs), dtype=dtp)
    fp_arr['imx'] = imxs
    fp_arr['imy'] = imys
    fp_arr['fname'] = fnames
    return fp_arr

def gauss_sig_bkg_nllh(cnts, nsig, nbkg, bkg_err, sys_err=0.0):

    sigma2 = nbkg + nsig + bkg_err**2 + (sys_err*nsig)**2

    N_sub_bkg = cnts - nbkg

    nllh = -1*np.sum(stats.norm.logpdf(N_sub_bkg - nsig,\
                                   scale=np.sqrt(sigma2)))

    return nllh



def gauss_nllh2min(theta, data_counts, bkg_counts,\
               bkg_err, cnts_intp):

    Nsig = 10.**theta[0]
    gamma = theta[1]

    Nsigs = Nsig*cnts_intp(gamma)

    return gauss_sig_bkg_nllh(data_counts, Nsigs, bkg_counts,\
                       bkg_err)


def gauss_nllh2min_2regs(theta, data_counts, bkg_counts, bkg_err,\
                         data_counts2, bkg_counts2,\
                         bkg_err2, cnts_intp, Ndet_frac,\
                         Atrans):

    Nsig = 10.**theta[0]
    gamma = theta[1]

    Nsigs = Nsig*cnts_intp(gamma)

    sys_err = 0.05*np.ones_like(bkg_counts)

    nllh0 = gauss_sig_bkg_nllh(data_counts, Nsigs*Ndet_frac, bkg_counts,\
                       bkg_err, sys_err=sys_err)

    if Ndet_frac > 0.99:
        return nllh0

    sys_err2 = np.sqrt(sys_err**2 + (Atrans/4.**2))

    nllh1 = gauss_sig_bkg_nllh(data_counts2, Nsigs*(1.-Ndet_frac)*Atrans, bkg_counts2,\
                       bkg_err2, sys_err=sys_err2)

    return nllh0 + nllh1


def min_det_fp_nllh(data_counts, model, fp_bl,\
                    params, PSnames, bkg_name,\
                    dt, cnts_intp,\
                    Ndets, solid_ang_dmean, rt_sums,\
                    get_bkg_nllh=False):

    # bkg_rate, bkg_rate_err = bkg_obj.get_rate(t)

    nebins = model.nebins
    Nsrcs = len(PSnames)

    bkg_dpis = model.get_rate_dpis(params)
    bkg_cnts = np.array([np.sum(dpi[fp_bl])*dt for dpi in bkg_dpis])

    bkg_rate_errs = np.array([params['err_' + bkg_name + 'bkg_rate_' + str(j)] for j\
                              in range(nebins)])

    bkg_flatAs = np.array([params[bkg_name + 'flat_' + str(j)] for j in range(nebins)])
    bkg_diffAs = 1. - bkg_flatAs
    bkg_err = bkg_rate_errs*(bkg_flatAs*Ndets + bkg_diffAs*solid_ang_dmean)*dt

    tot_err2 = np.zeros(nebins)
    tot_err2 += bkg_err**2
    for i in range(Nsrcs):
        ps_rate_errs = np.array([params[PSnames[i] + '_rate_' + str(j)]\
                                 for j in range(nebins)])
        ps_err = ps_rate_errs*rt_sums[i]*dt
        tot_err2 += ps_err**2

    bkg_err = np.sqrt(tot_err2)

    args = (data_counts, bkg_cnts, bkg_err, cnts_intp)

    lowers = [-2., .25]
    uppers = [4., 2.25]
    bounds = optimize.Bounds(np.array(lowers), np.array(uppers))

#     x0s = [[1., 1.], [2., 1.],
#            [1., 2.], [2., 2.]]
    x0s = [[1., .725], [2., 1.105],
           [1., 1.605], [2., 1.995]]

#     x0 = [1., 1.]
    ress = []
    nlogls = np.zeros(len(x0s))

    for j, x0 in enumerate(x0s):
        res = optimize.minimize(gauss_nllh2min, x0, args=args,\
                            method='L-BFGS-B', bounds=bounds)

#         print res
        ress.append(res)
        nlogls[j] = res.fun

    if np.all(np.isnan(nlogls)):
        best_ind = 0
    else:
        best_ind = np.nanargmin(nlogls)

    bf_nsig = 10.**ress[best_ind].x[0]
    # bf_nsig = ress[best_ind].x[0]
    bf_ind = ress[best_ind].x[1]

    if get_bkg_nllh:

        bkg_nllh = gauss_sig_bkg_nllh(data_counts, 0., bkg_cnts, bkg_err)
        return bf_nsig, bf_ind, nlogls[best_ind], bkg_nllh

    return bf_nsig, bf_ind, nlogls[best_ind]


def min_det_fp_nllh2(data_counts0, data_counts1, model, fp_bl,\
                    params, PSnames, bkg_name,\
                    dt, cnts_intp,\
                    Ndets0, Ndets1, solid_ang_dmean0, solid_ang_dmean1,\
                    rt_sums0, rt_sums1, Atrans,\
                    get_bkg_nllh=False):

    # bkg_rate, bkg_rate_err = bkg_obj.get_rate(t)

    nebins = model.nebins
    Nsrcs = len(PSnames)

    bkg_dpis = model.get_rate_dpis(params)
    bkg_cnts0 = np.array([np.sum(dpi[fp_bl])*dt for dpi in bkg_dpis])
    bkg_cnts1 = np.array([np.sum(dpi[~fp_bl])*dt for dpi in bkg_dpis])


    bkg_rate_errs = np.array([params['err_' + bkg_name + 'bkg_rate_' + str(j)] for j\
                              in range(nebins)])

    bkg_flatAs = np.array([params[bkg_name + 'flat_' + str(j)] for j in range(nebins)])
    bkg_diffAs = 1. - bkg_flatAs
    bkg_err0 = bkg_rate_errs*(bkg_flatAs*Ndets0 + bkg_diffAs*solid_ang_dmean0)*dt
    bkg_err1 = bkg_rate_errs*(bkg_flatAs*Ndets1 + bkg_diffAs*solid_ang_dmean1)*dt

    tot_err02 = np.zeros(nebins)
    tot_err02 += bkg_err0**2
    tot_err12 = np.zeros(nebins)
    tot_err12 += bkg_err1**2
    for i in range(Nsrcs):
        ps_rate_errs = np.array([params['err_' + PSnames[i] + '_rate_' + str(j)]\
                                 for j in range(nebins)])
        ps_err0 = ps_rate_errs*rt_sums0[i]*dt
        tot_err02 += ps_err0**2
        ps_err1 = ps_rate_errs*rt_sums1[i]*dt
        tot_err12 += ps_err1**2

    bkg_err0 = np.sqrt(tot_err02)
    bkg_err1 = np.sqrt(tot_err12)

#     args = (data_counts, bkg_cnts, bkg_err, cnts_intp)
    Ndet_frac = Ndets0/float(Ndets0+Ndets1)
    args = (data_counts0, bkg_cnts0, bkg_err0, data_counts1,\
            bkg_cnts1, bkg_err1, cnts_intp, Ndet_frac, Atrans)

    lowers = [-2., .25]
    uppers = [10., 2.25]
    bounds = optimize.Bounds(np.array(lowers), np.array(uppers))

    # x0s = [[1.5, .725], [2., 1.105],
    #        [2.25, 1.605], [2.5, 1.995]]
    x0s = [[2.2, .725], [3.7, 1.605]]

#     x0 = [1., 1.]
    ress = []
    nlogls = np.zeros(len(x0s))

    for j, x0 in enumerate(x0s):
        res = optimize.minimize(gauss_nllh2min_2regs, x0, args=args,\
                            method='L-BFGS-B', bounds=bounds)

#         print res
        ress.append(res)
        nlogls[j] = res.fun

    if np.all(np.isnan(nlogls)):
        best_ind = 0
    else:
        best_ind = np.nanargmin(nlogls)

    bf_nsig = 10.**ress[best_ind].x[0]
    # bf_nsig = ress[best_ind].x[0]
    bf_ind = ress[best_ind].x[1]

    if get_bkg_nllh:

#         bkg_nllh0 = gauss_sig_bkg_nllh(data_counts0, 0., bkg_cnts0, bkg_err0)
#         bkg_nllh1 = gauss_sig_bkg_nllh(data_counts1, 0., bkg_cnts1, bkg_err1)
        bkg_nllh = gauss_nllh2min_2regs([-10.0, 1.5], *args)
        return bf_nsig, bf_ind, nlogls[best_ind], bkg_nllh

    return bf_nsig, bf_ind, nlogls[best_ind]


def get_cnts_from_tbins_ebins(evdata, tbins0, tbins1, ebins0, ebins1):

    nebins = len(ebins0)
    hist = np.zeros((len(tbins0),nebins), dtype=np.int)

    ebins = np.append(ebins0, [ebins1[-1]])
    ntbins = len(tbins0)
    tbin_size = tbins1[0] - tbins0[0]
    tstep = tbins0[1] - tbins0[0]
    dts = np.diff(tbins0)
    breaks = np.where(dts>tstep)
    tfreq = int(np.rint(tbin_size/tstep))
    i0 = 0

    for bi in range(len(breaks[0])+1):

        if bi >= len(breaks[0]):
            break_ = len(tbins0)
        else:
            break_ = breaks[0][bi]
        tbins0_ = tbins0[i0:break_+1]
        t_add = [tbins0_[-1] + (i+1)*tstep for i in range(tfreq)]
        tbins = np.append(tbins0_, t_add)
        h = np.histogramdd([evdata['TIME'],evdata['ENERGY']], bins=[tbins,ebins],\
                        range=((tbins[0],tbins[-1]),(ebins[0],ebins[-1])))[0]

        if tfreq <= 1:
            h2 = h
        else:
            h2 = np.zeros((h.shape[0]-(tfreq-1),h.shape[1]))
            for i in range(tfreq):
                i00 = i
                i11 = -tfreq + 1 + i
                if i11 < 0:
                    h2 += h[i00:i11]
                else:
                    h2 += h[i00:]

        hist[i0:break_+1] = h2
        i0 = break_+1

    return hist



class rates_fp_llh(object):

    def __init__(self, imxs, imys,\
                 ev_data, twind_tab,\
                 ebins0, ebins1, fp_dir,\
                 bl_dmask, drm_dir, bkg_fname,\
                 solid_ang_dpi, rt_dir):

        self.bkg_df, self.bkg_name, self.PSnames, self.bkg_mod, self.ps_mods =\
        parse_bkg_csv(bkg_fname, solid_ang_dpi, ebins0, ebins1, bl_dmask, rt_dir)
        self.Nsrcs = len(self.PSnames)
        if self.Nsrcs < 1:
            self.model = self.bkg_mod
            self.ray_traces = []
        else:
            self.model_list = [self.bkg_mod]
            self.model_list += self.ps_mods
            self.model = CompoundModel(self.model_list)
            self.ray_traces = [self.ps_mods[j].get_rt(self.ps_mods[j].imx,\
                                self.ps_mods[j].imy) for j in range(self.Nsrcs)]

        self.Nfps = len(imxs)
        self.fp_dir = fp_dir
        self.fp_obj = FootPrints(fp_dir)
        self.imxs = imxs
        self.imys = imys
        self.ebins0 = ebins0
        self.ebins1 = ebins1
        self.nebins = len(self.ebins0)
        self.twind_tab = twind_tab
        self.exp_groups = self.twind_tab.groupby('duration')
        self.Ndurs = len(self.exp_groups)

        # self.t_bins0 = t_bins0
        # self.t_bins1 = t_bins1
        self.ev_data = ev_data
        self.bl_dmask = bl_dmask
        self.ind_ax = np.linspace(-.5, 2.5, 20*3+1)
        self.drm_obj = DRMs(drm_dir)
        self.Ndets_tot = np.sum(bl_dmask)
        # self.bkg_mod = bkg_mod
        # self.bkg_df = bkg_df

        self.solid_ang_dpi = solid_ang_dpi
        self.solid_angle_tot = np.sum(solid_ang_dpi[bl_dmask])
        self.solid_angle_mean = np.mean(solid_ang_dpi[bl_dmask])
        self.solid_angs_dmean = self.solid_ang_dpi[bl_dmask]/self.solid_angle_mean
        self.solid_angs_dmean_sum = np.sum(self.solid_angs_dmean)

        # self.Atrans = np.array([.05, .06, .07, .08, .09, .1])
        self.Atrans = np.array([.05, .09, .125, .25, .4, .65])

    def get_fp_vals(self):

        self.fp_bls = []
        self.fpbls = []
        self.ndets = []
        self.solid_angs = []
        self.rt_sums = []
        self.rt_sums1 = []
        for i in range(self.Nfps):
#             fp_ind = np.argmin(im_dist(self.fp_arr['imx'], self.fp_arr['imy'],
#                                self.imxs[i], self.imys[i]))
#             fp = np.load(os.path.join(self.fp_dir,\
#                     self.fp_arr[fp_ind]['fname']))
            fp = self.fp_obj.get_fp(self.imxs[i], self.imys[i])
            self.fp_bls.append(mask_detxy(fp, self.ev_data))
            self.ndets.append(np.sum(self.bl_dmask&(fp==1)))
#             fpbl = (self.bl_dmask&(fp==1))
            fpbl = ((fp[self.bl_dmask]==1))
            self.fpbls.append(fpbl)
            self.solid_angs.append(np.sum(self.solid_angs_dmean[fpbl]))
            rtsums = []
            rtsums1 = []
            for j in range(self.Nsrcs):
                rtsums.append(np.sum(self.ray_traces[j][fpbl]))
                rtsums1.append(np.sum(self.ray_traces[j][~fpbl]))
            self.rt_sums.append(rtsums)
            self.rt_sums1.append(rtsums1)


#             self.ndets.append(np.sum(self.bl_dmask&(fp==0)))


    def get_cnts_tbins_ebins_fps(self, dur):

        # gaps in twinds might mess this up

        df_twind = self.exp_groups.get_group(dur)
        tbins0 = df_twind['time'].values
        tbins1 = df_twind['time_end'].values

        # tbins0 = self.t_bins0[dur_ind]
        # tbins1 = self.t_bins1[dur_ind]
        ntbins = len(tbins0)
#         tbin_size = tbins1[0] - tbins0[0]
#         tstep = tbins0[1] - tbins0[0]
#         tfreq = int(np.rint(tbin_size/tstep))
#         t_add = [tbins0[-1] + (i+1)*tstep for i in range(tfreq)]
#         tbins = np.append(tbins0, t_add)
#         ebins = np.append(self.ebins0, [self.ebins1[-1]])

        self.cnts_tot = get_cnts_from_tbins_ebins(self.ev_data, tbins0, tbins1,\
                                                  self.ebins0, self.ebins1)

        logging.debug("Done getting cnts_tot")

        self.cnts_fpte = np.zeros((self.Nfps,ntbins,self.nebins))

#         self.cnts_tot = np.zeros((ntbins,self.nebins))

#         h = np.histogramdd([self.ev_data['TIME'],\
#                     self.ev_data['ENERGY']],
#                bins=[tbins,ebins])[0]

#         if tfreq <= 1:
#             h2 = h
#         else:
#             h2 = np.zeros((h.shape[0]-(tfreq-1),h.shape[1]))
#             for i in range(tfreq):
#                 i0 = i
#                 i1 = -tfreq + 1 + i
#                 if i1 < 0:
#                     h2 += h[i0:i1]
#                 else:
#                     h2 += h[i0:]
#         self.cnts_tot = h2


        for ii in range(self.Nfps):

            fp_bl = self.fp_bls[ii]
            self.cnts_fpte[ii] = get_cnts_from_tbins_ebins(self.ev_data[fp_bl],\
                                                           tbins0, tbins1,\
                                                           self.ebins0, self.ebins1)
            logging.debug("Done with %d of %d cnts_fptes"%(ii+1,self.Nfps))

#             h = np.histogramdd([self.ev_data['TIME'][fp_bl],\
#                                 self.ev_data['ENERGY'][fp_bl]],
#                            bins=[tbins,ebins])[0]

#             if tfreq <= 1:
#                 h2 = h
#             else:
#                 h2 = np.zeros((h.shape[0]-(tfreq-1),h.shape[1]))
#                 for i in range(tfreq):
#                     i0 = i
#                     i1 = -tfreq + 1 + i
#                     if i1 < 0:
#                         h2 += h[i0:i1]
#                     else:
#                         h2 += h[i0:]
#             self.cnts_fpte[ii] = h2


    def get_drm_stuff(self):

        self.cnts_intps = []
        for i in range(self.Nfps):
            imx = self.imxs[i]
            imy = self.imys[i]
            drm = self.drm_obj.get_drm(imx, imy)
            ebin_ind_edges = get_ebin_ind_edges(drm, self.ebins0, self.ebins1)
            abs_cor = get_abs_cor_rates(imx, imy, drm)
            self.cnts_intps.append(get_cnts_intp_obj(self.ind_ax,\
                                drm, ebin_ind_edges, abs_cor))



    def run(self):
        t_0 = time.time()
        self.get_fp_vals()
        self.get_drm_stuff()
        logging.info("Done setting up footprints and drm stuff")
        logging.info("Took %.3f seconds" %(time.time()-t_0))

        res_dicts = []

        for ii, exp_group in enumerate(self.exp_groups):

            logging.info("Starting duration size %d of %d" %(ii+1, self.Ndurs))

            dur = exp_group[0]
            df_twind = exp_group[1]

            tbins0 = df_twind['time'].values
            tbins1 = df_twind['time_end'].values
            timeIDs = df_twind['timeID'].values

            ntbins = len(tbins0)
            tbin_size = tbins1[0] - tbins0[0]
            tstep = tbins0[1] - tbins0[0]

            logging.debug("ntbins: %d" %(ntbins))
            logging.debug("tbin_size: %.3f" %(tbin_size))
            logging.debug("tstep: %.3f" %(tstep))


            t_0 = time.time()
            self.get_cnts_tbins_ebins_fps(dur)
            logging.info("Done getting cnts_fpte")
            logging.info("Took %.3f seconds" %(time.time()-t_0))

            for jj in range(self.Nfps):

                cnts_intp = self.cnts_intps[jj]
                cnts_per_tbin = self.cnts_fpte[jj]
                cnts_per_tbin1 = self.cnts_tot - cnts_per_tbin
                Ndets = self.ndets[jj]
                solid_ang = self.solid_angs[jj]
                rt_sums = self.rt_sums[jj]
                rt_sums1 = self.rt_sums1[jj]
                fpbl = self.fpbls[jj]
                imx = self.imxs[jj]
                imy = self.imys[jj]

                bf_nsigs = np.zeros(ntbins)
                bf_inds = np.zeros(ntbins)
                nllhs = np.zeros(ntbins)
                bkg_nllhs = np.zeros(ntbins)


                t_0 = time.time()

                for kk in range(ntbins):

                    res_dict = {'dur':tbin_size,
                                'imx':imx, 'imy':imy,
                                'ndets':Ndets,
                                'solid_angle':solid_ang,
                                'timeID':timeIDs[kk]}
                    for Nps in range(self.Nsrcs):
                        res_dict[self.PSnames[Nps] + '_rt_sum'] = rt_sums[Nps]

                    res_dict['time'] = tbins0[kk]

                    bkg_ind = np.argmin(np.abs((tbins0[kk]+dur/2.) -\
                                            self.bkg_df['time']))
                    bkg_row = self.bkg_df.iloc[bkg_ind]
#                     bkg_diffuse = np.array([bkg_row['diffuse_'+str(i)] for i\
#                                     in range(self.nebins)])
#                     bkg_flat = np.array([bkg_row['flat_'+str(i)] for i\
#                                     in range(self.nebins)])

#                     res_dict['Nsig'], res_dict['Plaw_Ind'], res_dict['nllh'],\
#                     res_dict['bkg_nllh'] = min_det_fp_nllh(cnts_per_tbin[kk],\
#                                self.model, fpbl, bkg_row, self.PSnames, self.bkg_name,\
#                                 tbin_size,\
#                                 cnts_intp, Ndets, solid_ang, rt_sums,\
#                                 get_bkg_nllh=True)

                    res_dict['Nsig'], res_dict['Plaw_Ind'], res_dict['nllh'], res_dict['bkg_nllh'] =\
                                    min_det_fp_nllh2(cnts_per_tbin[kk], cnts_per_tbin1[kk],\
                                                            self.model, fpbl,\
                                    bkg_row, self.PSnames, self.bkg_name,\
                                    tbin_size, cnts_intp,\
                                    Ndets, self.Ndets_tot-Ndets, solid_ang,\
                                    self.solid_angs_dmean_sum-solid_ang,\
                                    rt_sums, rt_sums1, self.Atrans,\
                                    get_bkg_nllh=True)

                    TS = np.sqrt(2.*(res_dict['bkg_nllh'] - res_dict['nllh']))
                    if np.isnan(TS):
                        TS = 0.0
                    res_dict['TS'] = TS
                    res_dicts.append(res_dict)


#                 print "Done minimizing in loop of tbins"
#                 print "Took %.3f seconds" %(time.time()-t_0)

                logging.info("Done with %d of %d positions for duration %d of %d"%\
                        (jj+1,self.Nfps,ii+1,self.Ndurs))


        return res_dicts


def main(args):

    logging.basicConfig(filename='rates_llh_analysis_%d.log' %(args.job_id),\
                    level=logging.DEBUG,\
                    format='%(asctime)s-' '%(levelname)s- %(message)s')

    if args.dbfname is None:
        db_fname = guess_dbfname()
        if isinstance(db_fname, list):
            db_fname = db_fname[0]
    else:
        db_fname = args.dbfname

    logging.info('Connecting to DB')
    conn = get_conn(db_fname)

    info_tab = get_info_tab(conn)
    logging.info('Got info table')

    files_tab = get_files_tab(conn)
    logging.info('Got files table')

    trigtime = info_tab['trigtimeMET'][0]

    drm_dir = files_tab['drmDir'][0]
    rt_dir = files_tab['rtDir'][0]

    evfname = files_tab['evfname'][0]
    dmfname = files_tab['detmask'][0]
    ev_data = fits.open(evfname, memmap=False)[1].data
    logging.debug('Opened up event file')
    dmask = fits.open(dmfname)[0].data
    bl_dmask = (dmask==0)
    logging.debug('Opened up dmask file')

    ebins0 = np.array(EBINS0)
    ebins1 = np.array(EBINS1)
    nebins = len(ebins0)
    logging.debug("ebins0")
    logging.debug(ebins0)
    logging.debug("ebins1")
    logging.debug(ebins1)

    # probably get times from twind table

    twind_df = get_twinds_tab(conn)

    logging.info("Got TimeWindows table")

    logging.info("Getting bkg estimation from file")

    # rate_fits_df = get_rate_fits_tab(conn)
    while True:
        try:
            bkg_fits_df = pd.read_csv(args.bkg_fname)
            break
        except:
            time.sleep(10.0)
    # bkg_obj = rate_obj_from_sqltab(rate_fits_df, 0, 1)
    logging.info("Got bkg estimation")

    min_bin_size = np.min(twind_df['duration'])

    logging.info("Smallest duration to test is %.3fs" %(min_bin_size))

    exp_groups = twind_df.groupby('duration')

    nexps = len(exp_groups)

    # imxax = np.arange(-1.6, 1.61, 0.1)
    # imyax = np.arange(-.9, .91, 0.1)
    imxax = np.arange(-1.6, 1.61, 0.16)
    imyax = np.arange(-.9, .91, 0.05)

    imxg, imyg = np.meshgrid(imxax, imyax)
    imxs = imxg.ravel()
    imys = imyg.ravel()
    bl = np.isclose((imys%.1),0)
    imxs[bl] += .08

    PC = fits.open(args.pcfname)[0]
    pc = PC.data
    w_t = WCS(PC.header, key='T')

    pcs = world2val(w_t, pc, imxs, imys)
    logging.info("min, max pcs: %.4f, %.4f"%(np.min(pcs),np.max(pcs)))
    min_pc = max(args.min_pc - 0.25, 0.0099)
    logging.info("min_pc: %.4f"%(min_pc))
    logging.info("sum(pcs>min_pc: %d"%(np.sum(pcs>min_pc)))
    pc_bl = (pcs>min_pc)
    imxs = imxs[pc_bl]
    imys = imys[pc_bl]

    try:
        good_pix = np.load(args.pix_fname)
        im_dists = np.zeros_like(imxs)
        if len(good_pix) < 1:
            logging.info("pix2scan file is there are 0 pixels to scan")
            logging.info("Exiting")
            return
        for i in range(len(imxs)):
            im_dists[i] = np.min(im_dist(imxs[i], imys[i],\
                            good_pix['imx'], good_pix['imy']))
        bl = (im_dists<.1)
        imxs = imxs[bl]
        imys = imys[bl]
    except Exception as E:
        logging.error(E)
        logging.warning("Trouble reading the pix2scan file")
        logging.info("Using whole FoV")

    Npnts = len(imxs)
    logging.info("%d total grid points" %(Npnts))
    Nper_job = 1 + int(Npnts/float(args.Njobs))
    if args.job_id > -1:
        i0 = args.job_id*Nper_job
        i1 = i0 + Nper_job
        imxs = imxs[i0:i1]
        imys = imys[i0:i1]
    logging.info("%d grid points to do" %(len(imxs)))

    solid_angle_dpi = np.load(solid_angle_dpi_fname)

    # bkg_mod = Bkg_Model_wSA(bl_dmask, solid_angle_dpi, nebins)

    rate_llh_obj = rates_fp_llh(imxs, imys, ev_data, twind_df,\
                                ebins0, ebins1, fp_dir,\
                                bl_dmask, drm_dir, args.bkg_fname,\
                                solid_angle_dpi, rt_dir)

    res_dicts = rate_llh_obj.run()
    logging.info("Done with analysis")
    logging.info("%d results to write" %(len(res_dicts)))

    # append_rate_tab(conn, df_twind, quad_dict['id'], bkg_llh_tbins, llhs, bf_nsigs, bf_inds)
    #
    # logging.info("Appended rate results to DB")

    df = pd.DataFrame(res_dicts)
    logging.info("Done making results into DataFrame")
    save_fname = 'rates_llh_res_%d_.csv' %(args.job_id)
    df.to_csv(save_fname, index=False)



if __name__ == "__main__":

    args = cli()

    main(args)