<reponame>jacobic/redpipes
import os
import numpy as np
import pandas as pd

from astropy import units as u
from astropy.coordinates import Angle, SkyCoord, concatenate
from astropy.table import Table, vstack

from src.utils import Str, wcs_pixel_info
import src.globals as glo


def match_combine_tables(table1, table2, sep_limit=Angle(0.25, u.arcsec),
                         primary_cut=250):
    """
    Combine tables after eliminating duplicate sources. A duplicate source is
    taken to be a source which appears in both `table1` and `table2` within a
    specified search radius.

    Read here for more information:
    https://outerspace.stsci.edu/display/PANSTARRS/PS1+Sky+tessellation+patterns

    The pixel size is 0.25 arcsec, which is by design very close to the native
    image scale of the detector (0.258'').

    Adjacent skycells within a projection cell overlap by 240 pixels (60
    arcsec) at each edge with the neighboring skycells.  Consequently,
    for objects and structures smaller than 2 arcmin it is generally not
    necessary to combine data from adjacent skycells.

    Parameters
    ----------
    table1 : astropy.table
    table2 :
    sep_limit : Angle, optional
        Search radius over which matching is carried out.

    Returns
    ----------
    table_combined: astropy.table
        A table containing only unique sources from the filtered and
        vertically stacked input tables.
    """

    #TODO: use TargetDataFrame class with sky and sex_keys, save 8 lines !!!!!

    ra1 = Angle(table1[glo.key_ra], u.deg)
    dec1 = Angle(table1[glo.key_dec], u.deg)
    ra2 = Angle(table2[glo.key_ra], u.deg)
    dec2 = Angle(table2[glo.key_dec], u.deg)
    sky1 = SkyCoord(ra=ra1.deg, dec=dec1.deg, unit='deg')
    sky2 = SkyCoord(ra=ra2.deg, dec=dec2.deg, unit='deg')
    idx1, idx2, _, _ = sky2.search_around_sky(sky1, sep_limit)

    s1x, s1y = table1['XWIN_IMAGE_I'], table1['YWIN_IMAGE_I']
    s2x, s2y = table2['XWIN_IMAGE_I'], table2['YWIN_IMAGE_I']

    dim_ra = np.absolute(sky2[idx2].ra.max() - sky2[idx2].ra.min())
    dim_dec = np.absolute(sky2[idx2].dec.max() - sky2[idx2].dec.min())

    try:
        if dim_ra > dim_dec:
            # on its side, therefore at combine it must be up to down
            # therefore trim in dec direction.
            assert (sky1.dec.mean() < sky2.dec.mean()), 'Not ordered correctly.'
            # sky1 is below sky2.
            s1, s2 = s1y, s2y
        else:
            # Stood up therefore combine must be side to side
            # therefore trim in ra direction.
            assert (sky1.ra.mean() < sky2.ra.mean()), 'Not ordered correctly.'
            # sky1 is to the left of sky2.
            # XWIN_IMAGE is mirrored w.r.t RA.
            s1, s2, = -s1x, -s2x

        cut1 = s1 < (np.max(s1) - primary_cut)
        cut2 = s2 > (np.min(s2) + primary_cut)
        table1 = table1[cut1]
        table2 = table2[cut2]

    except AssertionError as e:
        return match_combine_tables(table2, table1)

    #TODO: use TargetDataFrame class with sky and sex_keys, save 8 lines !!!!!
    ra1 = Angle(table1[glo.key_ra], u.deg)
    dec1 = Angle(table1[glo.key_dec], u.deg)
    ra2 = Angle(table2[glo.key_ra], u.deg)
    dec2 = Angle(table2[glo.key_dec], u.deg)

    sky1 = SkyCoord(ra=ra1.deg, dec=dec1.deg, unit='deg')
    sky2 = SkyCoord(ra=ra2.deg, dec=dec2.deg, unit='deg')
    idx1, idx2, _, _ = sky2.search_around_sky(sky1, sep_limit)

    # Eliminate matches with boolean mask. These matches represent duplicates.
    mask1 = np.ones(len(table1), dtype=bool)
    mask1[idx1] = False
    cleaned_combined = vstack([table1[mask1], table2])
    return cleaned_combined


def calibrate_table(table, transform='transform_PS1-DES'):
    """

    Parameters
    ----------
    table : astropy.table
        The uncalibrated table.

    Returns
    ----------
    table : astropy.table.Table
        The calibrated table.
    """
    mag_auto_i = table['MAG_AUTO_I']
    mag_auto_z = table['MAG_AUTO_Z']

    limit_mag_auto_i = (mag_auto_i > 14.0) & (mag_auto_i < 28.0)
    limit_mag_auto_z = (mag_auto_z > 3.0) & (mag_auto_z < 28.0)
    table = table[limit_mag_auto_i & limit_mag_auto_z]
    transform = os.path.join(os.path.dirname(__file__), transform)
    calibration = pd.read_table(transform, delim_whitespace=True, header=0)

    for column in calibration.columns:
        table[column] = table[column] - calibration[column].values[0]

    # Clean bad mag detmodel values with mag auto.
    for col in [Str('gr'), Str('ri'), Str('iz')]:
        mag_auto = table['MAG_AUTO_{0:u}'.format(col)]
        mag_detmodel = table['MAG_DETMODEL_{0:u}'.format(col)]
        limit_mag = mag_detmodel - mag_auto > 7
        mag_detmodel[limit_mag] = mag_auto[limit_mag]
        table['MAG_DETMODEL_{0:u}'.format(col)] = mag_detmodel
    return table


def combine_batches(batches, initial=False):
    """

    Designed to process batches, i.e. a nested list of Table
    objects. Batches on the first level of nesting should be ordered by RA
    and each batch within batches should be ordered by Dec (or the other way
    round).

    i.e. this function recursively merges 9 skycell tables which have been
    split into batches as illustrated below. The x marks the cluster centre.

    ────────────────
    │    │    │    │
    │    │    │    │    batches[2]  = [Table20, Table21, Table22]
    │────│────│────│
    │    │  x │    │
    │    │    │    │    batches[1]  = [Table10, Table11, Table12]
    │────│────│────│
    │    │    │    │
    │    │    │    │    batches[0]  = [Table00, Table01, Table02]
    │────│────│────│


    Parameters
    ----------
    batches

    Returns
    -------

    """
    if initial is True:
        combined = []
        for b in batches:
            combined.append(combine_batches(b))
        return combine_batches(combined)
    else:
        part = []
        for b in batches[:2], batches[2:]:
            if not b:
                continue
            elif len(b) > 2:
                part.append(combine_batches(b))
            elif len(b) == 2:
                part.append(match_combine_tables(b[0], b[1]))
            else:
                part.append(b[0])
        try:
            return match_combine_tables(part[0], part[1])
        except IndexError as e:
            return part[0]


def ordered_tables_to_combine(updated_catalogs):
    """

    Parameters
    ----------
    updated_catalogs

    Returns
    -------

    """
    ordered_cols = sorted(updated_catalogs, key=lambda x: x.sky_mid.dec.deg)
    binsize = 0.38
    data = np.array([cat.sky_mid.dec.deg for cat in updated_catalogs])
    bins = np.arange(np.min(data), np.max(data) + binsize, binsize)
    hist, edges = np.histogram(data, bins=bins)
    hist_reduced = [i for i in hist if i != 0]

    def row(l, indicies):
        for i, idx in enumerate(indicies):
            try:
                yield sorted(l[idx:indicies[i + 1]], key=lambda x: x.sky_mid.ra.deg)
            except IndexError as e:
                continue

    ord_cells = list(row(ordered_cols, indicies=np.cumsum([0] + hist_reduced)))
    ord_tables = list(map(lambda x: [y.table() for y in x], ord_cells))
    return ord_tables


def calibrate_combine(c, files=None, trim=Angle(0.05, u.deg), write=True):
    """
    Calibrate and combine catalogs.

    batch1
    334.751732842 30.9949808447 ----/
    334.759649924 31.3937820752 ----/ half_combined[0]

    batch2
    335.079177115 30.9933694087 ----/
    335.069144661 31.386782124 ----/ half_combined[1]

    Parameters
    ----------
    trim

    Returns
    ----------
    None
    """

    try:
        if not files:
            tables = ordered_tables_to_combine(updated_catalogs=c.files_updated)
        else:
            tables = ordered_tables_to_combine(updated_catalogs=files)
        c.log.info('Attempting to combine files for {0}'.format(c))
        # list of catalogs to calibrate and combine sorted by RA.
        assert tables is not None, 'Tables is None'

        final_table = combine_batches(tables, initial=True)
        # Calibrate final table to DES colours.
        final_table = calibrate_table(final_table)

        # Now trim the perimeter of the final catalog.
        ra = Angle(final_table[glo.key_ra], u.deg)
        dec = Angle(final_table[glo.key_dec], u.deg)
        ra_trim = (ra < np.max(ra) - trim) & (ra > np.min(ra) + trim)
        dec_trim = (dec < np.max(dec) - trim) & (dec > np.min(dec) + trim)
        idx_trim = np.where(ra_trim & dec_trim)[0]

        # Write out final combined, cut, calibrated and trimmed catalog.
        if write:
            name_out = '{0}{1}'.format(c, glo.final_name)
            data_out = os.path.join(c.entry, 'processed', name_out)
            final_table[idx_trim].write(data_out, format='fits', overwrite=True)
            c.log.info('Successful combine files for {0}'.format(c))

        return final_table[idx_trim]

    except TypeError:
        c.log.exception("{0} caused a type error.".format(c))

    except AttributeError as e:
        c.log.exception("{0} caused attribute error.".format(c))

    except AssertionError as e:
        c.log.exception("{0} failed an assertion test.".format(c))

    except Exception as e:
        c.log.exception("{0} Unknown reason for failing.".format(c))

import os
import numpy as np
import pandas as pd

from astropy import units as u
from astropy.coordinates import Angle, SkyCoord, concatenate
from astropy.table import Table, vstack

from src.utils import Str, wcs_pixel_info
import src.globals as glo


def match_combine_tables(table1, table2, sep_limit=Angle(0.25, u.arcsec),
                         primary_cut=250):
    """
    Combine tables after eliminating duplicate sources. A duplicate source is
    taken to be a source which appears in both `table1` and `table2` within a
    specified search radius.

    Read here for more information:
    https://outerspace.stsci.edu/display/PANSTARRS/PS1+Sky+tessellation+patterns

    The pixel size is 0.25 arcsec, which is by design very close to the native
    image scale of the detector (0.258'').

    Adjacent skycells within a projection cell overlap by 240 pixels (60
    arcsec) at each edge with the neighboring skycells.  Consequently,
    for objects and structures smaller than 2 arcmin it is generally not
    necessary to combine data from adjacent skycells.

    Parameters
    ----------
    table1 : astropy.table
    table2 :
    sep_limit : Angle, optional
        Search radius over which matching is carried out.

    Returns
    ----------
    table_combined: astropy.table
        A table containing only unique sources from the filtered and
        vertically stacked input tables.
    """

    #TODO: use TargetDataFrame class with sky and sex_keys, save 8 lines !!!!!

    ra1 = Angle(table1[glo.key_ra], u.deg)
    dec1 = Angle(table1[glo.key_dec], u.deg)
    ra2 = Angle(table2[glo.key_ra], u.deg)
    dec2 = Angle(table2[glo.key_dec], u.deg)
    sky1 = SkyCoord(ra=ra1.deg, dec=dec1.deg, unit='deg')
    sky2 = SkyCoord(ra=ra2.deg, dec=dec2.deg, unit='deg')
    idx1, idx2, _, _ = sky2.search_around_sky(sky1, sep_limit)

    s1x, s1y = table1['XWIN_IMAGE_I'], table1['YWIN_IMAGE_I']
    s2x, s2y = table2['XWIN_IMAGE_I'], table2['YWIN_IMAGE_I']

    dim_ra = np.absolute(sky2[idx2].ra.max() - sky2[idx2].ra.min())
    dim_dec = np.absolute(sky2[idx2].dec.max() - sky2[idx2].dec.min())

    try:
        if dim_ra > dim_dec:
            # on its side, therefore at combine it must be up to down
            # therefore trim in dec direction.
            assert (sky1.dec.mean() < sky2.dec.mean()), 'Not ordered correctly.'
            # sky1 is below sky2.
            s1, s2 = s1y, s2y
        else:
            # Stood up therefore combine must be side to side
            # therefore trim in ra direction.
            assert (sky1.ra.mean() < sky2.ra.mean()), 'Not ordered correctly.'
            # sky1 is to the left of sky2.
            # XWIN_IMAGE is mirrored w.r.t RA.
            s1, s2, = -s1x, -s2x

        cut1 = s1 < (np.max(s1) - primary_cut)
        cut2 = s2 > (np.min(s2) + primary_cut)
        table1 = table1[cut1]
        table2 = table2[cut2]

    except AssertionError as e:
        return match_combine_tables(table2, table1)

    #TODO: use TargetDataFrame class with sky and sex_keys, save 8 lines !!!!!
    ra1 = Angle(table1[glo.key_ra], u.deg)
    dec1 = Angle(table1[glo.key_dec], u.deg)
    ra2 = Angle(table2[glo.key_ra], u.deg)
    dec2 = Angle(table2[glo.key_dec], u.deg)

    sky1 = SkyCoord(ra=ra1.deg, dec=dec1.deg, unit='deg')
    sky2 = SkyCoord(ra=ra2.deg, dec=dec2.deg, unit='deg')
    idx1, idx2, _, _ = sky2.search_around_sky(sky1, sep_limit)

    # Eliminate matches with boolean mask. These matches represent duplicates.
    mask1 = np.ones(len(table1), dtype=bool)
    mask1[idx1] = False
    cleaned_combined = vstack([table1[mask1], table2])
    return cleaned_combined


def calibrate_table(table, transform='transform_PS1-DES'):
    """

    Parameters
    ----------
    table : astropy.table
        The uncalibrated table.

    Returns
    ----------
    table : astropy.table.Table
        The calibrated table.
    """
    mag_auto_i = table['MAG_AUTO_I']
    mag_auto_z = table['MAG_AUTO_Z']

    limit_mag_auto_i = (mag_auto_i > 14.0) & (mag_auto_i < 28.0)
    limit_mag_auto_z = (mag_auto_z > 3.0) & (mag_auto_z < 28.0)
    table = table[limit_mag_auto_i & limit_mag_auto_z]
    transform = os.path.join(os.path.dirname(__file__), transform)
    calibration = pd.read_table(transform, delim_whitespace=True, header=0)

    for column in calibration.columns:
        table[column] = table[column] - calibration[column].values[0]

    # Clean bad mag detmodel values with mag auto.
    for col in [Str('gr'), Str('ri'), Str('iz')]:
        mag_auto = table['MAG_AUTO_{0:u}'.format(col)]
        mag_detmodel = table['MAG_DETMODEL_{0:u}'.format(col)]
        limit_mag = mag_detmodel - mag_auto > 7
        mag_detmodel[limit_mag] = mag_auto[limit_mag]
        table['MAG_DETMODEL_{0:u}'.format(col)] = mag_detmodel
    return table


def combine_batches(batches, initial=False):
    """

    Designed to process batches, i.e. a nested list of Table
    objects. Batches on the first level of nesting should be ordered by RA
    and each batch within batches should be ordered by Dec (or the other way
    round).

    i.e. this function recursively merges 9 skycell tables which have been
    split into batches as illustrated below. The x marks the cluster centre.

    ────────────────
    │    │    │    │
    │    │    │    │    batches[2]  = [Table20, Table21, Table22]
    │────│────│────│
    │    │  x │    │
    │    │    │    │    batches[1]  = [Table10, Table11, Table12]
    │────│────│────│
    │    │    │    │
    │    │    │    │    batches[0]  = [Table00, Table01, Table02]
    │────│────│────│


    Parameters
    ----------
    batches

    Returns
    -------

    """
    if initial is True:
        combined = []
        for b in batches:
            combined.append(combine_batches(b))
        return combine_batches(combined)
    else:
        part = []
        for b in batches[:2], batches[2:]:
            if not b:
                continue
            elif len(b) > 2:
                part.append(combine_batches(b))
            elif len(b) == 2:
                part.append(match_combine_tables(b[0], b[1]))
            else:
                part.append(b[0])
        try:
            return match_combine_tables(part[0], part[1])
        except IndexError as e:
            return part[0]


def ordered_tables_to_combine(updated_catalogs):
    """

    Parameters
    ----------
    updated_catalogs

    Returns
    -------

    """
    ordered_cols = sorted(updated_catalogs, key=lambda x: x.sky_mid.dec.deg)
    binsize = 0.38
    data = np.array([cat.sky_mid.dec.deg for cat in updated_catalogs])
    bins = np.arange(np.min(data), np.max(data) + binsize, binsize)
    hist, edges = np.histogram(data, bins=bins)
    hist_reduced = [i for i in hist if i != 0]

    def row(l, indicies):
        for i, idx in enumerate(indicies):
            try:
                yield sorted(l[idx:indicies[i + 1]], key=lambda x: x.sky_mid.ra.deg)
            except IndexError as e:
                continue

    ord_cells = list(row(ordered_cols, indicies=np.cumsum([0] + hist_reduced)))
    ord_tables = list(map(lambda x: [y.table() for y in x], ord_cells))
    return ord_tables


def calibrate_combine(c, files=None, trim=Angle(0.05, u.deg), write=True):
    """
    Calibrate and combine catalogs.

    batch1
    334.751732842 30.9949808447 ----/
    334.759649924 31.3937820752 ----/ half_combined[0]

    batch2
    335.079177115 30.9933694087 ----/
    335.069144661 31.386782124 ----/ half_combined[1]

    Parameters
    ----------
    trim

    Returns
    ----------
    None
    """

    try:
        if not files:
            tables = ordered_tables_to_combine(updated_catalogs=c.files_updated)
        else:
            tables = ordered_tables_to_combine(updated_catalogs=files)
        c.log.info('Attempting to combine files for {0}'.format(c))
        # list of catalogs to calibrate and combine sorted by RA.
        assert tables is not None, 'Tables is None'

        final_table = combine_batches(tables, initial=True)
        # Calibrate final table to DES colours.
        final_table = calibrate_table(final_table)

        # Now trim the perimeter of the final catalog.
        ra = Angle(final_table[glo.key_ra], u.deg)
        dec = Angle(final_table[glo.key_dec], u.deg)
        ra_trim = (ra < np.max(ra) - trim) & (ra > np.min(ra) + trim)
        dec_trim = (dec < np.max(dec) - trim) & (dec > np.min(dec) + trim)
        idx_trim = np.where(ra_trim & dec_trim)[0]

        # Write out final combined, cut, calibrated and trimmed catalog.
        if write:
            name_out = '{0}{1}'.format(c, glo.final_name)
            data_out = os.path.join(c.entry, 'processed', name_out)
            final_table[idx_trim].write(data_out, format='fits', overwrite=True)
            c.log.info('Successful combine files for {0}'.format(c))

        return final_table[idx_trim]

    except TypeError:
        c.log.exception("{0} caused a type error.".format(c))

    except AttributeError as e:
        c.log.exception("{0} caused attribute error.".format(c))

    except AssertionError as e:
        c.log.exception("{0} failed an assertion test.".format(c))

    except Exception as e:
        c.log.exception("{0} Unknown reason for failing.".format(c))
Combine tables after eliminating duplicate sources. A duplicate source is     taken to be a source which appears in both `table1` and `table2` within a     specified search radius.      Read here for more information:     https://outerspace.stsci.edu/display/PANSTARRS/PS1+Sky+tessellation+patterns      The pixel size is 0.25 arcsec, which is by design very close to the native     image scale of the detector (0.258'').      Adjacent skycells within a projection cell overlap by 240 pixels (60     arcsec) at each edge with the neighboring skycells.  Consequently,     for objects and structures smaller than 2 arcmin it is generally not     necessary to combine data from adjacent skycells.      Parameters     ----------     table1 : astropy.table     table2 :     sep_limit : Angle, optional         Search radius over which matching is carried out.      Returns     ----------     table_combined: astropy.table         A table containing only unique sources from the filtered and         vertically stacked input tables. #TODO: use TargetDataFrame class with sky and sex_keys, save 8 lines !!!!! # on its side, therefore at combine it must be up to down # therefore trim in dec direction. # sky1 is below sky2. # Stood up therefore combine must be side to side # therefore trim in ra direction. # sky1 is to the left of sky2. # XWIN_IMAGE is mirrored w.r.t RA. #TODO: use TargetDataFrame class with sky and sex_keys, save 8 lines !!!!! # Eliminate matches with boolean mask. These matches represent duplicates. Parameters     ----------     table : astropy.table         The uncalibrated table.      Returns     ----------     table : astropy.table.Table         The calibrated table. # Clean bad mag detmodel values with mag auto. Designed to process batches, i.e. a nested list of Table     objects. Batches on the first level of nesting should be ordered by RA     and each batch within batches should be ordered by Dec (or the other way     round).      i.e. this function recursively merges 9 skycell tables which have been     split into batches as illustrated below. The x marks the cluster centre.      ────────────────     │    │    │    │     │    │    │    │    batches[2]  = [Table20, Table21, Table22]     │────│────│────│     │    │  x │    │     │    │    │    │    batches[1]  = [Table10, Table11, Table12]     │────│────│────│     │    │    │    │     │    │    │    │    batches[0]  = [Table00, Table01, Table02]     │────│────│────│       Parameters     ----------     batches      Returns     ------- Parameters     ----------     updated_catalogs      Returns     ------- Calibrate and combine catalogs.      batch1     334.751732842 30.9949808447 ----/     334.759649924 31.3937820752 ----/ half_combined[0]      batch2     335.079177115 30.9933694087 ----/     335.069144661 31.386782124 ----/ half_combined[1]      Parameters     ----------     trim      Returns     ----------     None # list of catalogs to calibrate and combine sorted by RA. # Calibrate final table to DES colours. # Now trim the perimeter of the final catalog. # Write out final combined, cut, calibrated and trimmed catalog.