import math
import sys

import numpy
import scipy.stats

import Likelihood
import cellDetect

try:

    import astropy.io.fits as pyfits

except:

    import pyfits

import logging

log = logging.getLogger("main")


class FixedBinSearch(object):

    def __init__(self, X, Y, t, timebin, eventfile, expomap, hwu, region=None):

        # Store the hardware unit
        self.hwu = hwu

        # Read the events
        self.X = X
        self.Y = Y
        self.time = t

        self.timebin = float(timebin)

        # Read the source catalog
        sources = []
        with pyfits.open(eventfile) as f:

            srcX = []
            srcY = []

            for row in f['REGION'].data:
                srcX.append(row.field("X")[0])
                srcY.append(row.field("Y")[0])

        sources = numpy.array(zip(srcX, srcY))

        # Transform RA,Dec in X,Y for the sources

        self.eventfile = eventfile

        # wcs = XMMWCS.XMMWCS(self.eventfile)

        # sources = numpy.array( wcs.sky2xy(zip(ras,decs)) )

        log.debug("Read %s sources from the source catalog" % (len(sources)))

        if region is None:

            self.xmin = self.X.min()
            self.xmax = self.X.max()

            self.ymin = self.Y.min()
            self.ymax = self.Y.max()

            # Now keep only the sources within this region
            idx = ((sources[:, 0] >= self.xmin) &
                   (sources[:, 0] <= self.xmax) &
                   (sources[:, 1] >= self.ymin) &
                   (sources[:, 1] <= self.ymax))

            self.sources = sources[idx]

        else:

            # Use the region to filter the sources

            filt = region.get_filter()

            self.sources = []

            for src in sources:

                if filt.inside1(src[0], src[1]):
                    self.sources.append(src)

            # Use the region to determine the extreme of X and Y

            shape = region[0]

            # NOTE: Assuming all units are in pixels

            if shape.name == 'circle':

                xc, yc, radius = shape.coord_list

                # sometimes for no reason the radius is negative. Fix that
                radius = abs(radius)

            elif shape.name == 'ellipse':

                xc, yc, r1, r2, ang = shape.coord_list

                radius = max(abs(r1), abs(r2))

            else:

                raise RuntimeError("Region of shape %s is not supported" % (shape.name))

            pass

            self.xmin = xc - radius
            self.xmax = xc + radius

            self.ymin = yc - radius
            self.ymax = yc + radius

        nsrc = len(self.sources)

        nmax = 100

        if nsrc > nmax:
            log.info("There were %s sources in the region, they are too many. Keeping only %s" % (nsrc, nmax))

            self.sources = self.sources[:nmax]

        log.debug("Kept %s source(s) within the region" % (len(self.sources)))

        self.expomap = expomap

    def scan(self, regfilt, threshold=5.0):

        bins = numpy.arange(self.time.min(), self.time.max() + 1e-4, self.timebin)

        xbins = numpy.arange(self.xmin, self.xmax, 25)
        ybins = numpy.arange(self.ymin, self.ymax, 25)