<commit_before>from __future__ import print_function
import pylab as plt
import numpy as np
from astrometry.util.ttime import Time

def _detmap(X):
    from scipy.ndimage.filters import gaussian_filter
    from legacypipe.survey import tim_get_resamp
    (tim, targetwcs, H, W, apodize) = X
    R = tim_get_resamp(tim, targetwcs)
    if R is None:
        return None,None,None,None,None
    assert(tim.psf_sigma > 0)
    psfnorm = 1./(2. * np.sqrt(np.pi) * tim.psf_sigma)
    ie = tim.getInvvar()
    detim = tim.getImage().copy()
    tim.getSky().addTo(detim, scale=-1.)
    detim = gaussian_filter(detim, tim.psf_sigma) / psfnorm**2
    detsig1 = tim.sig1 / psfnorm
    subh,subw = tim.shape
    detiv = np.zeros((subh,subw), np.float32) + (1. / detsig1**2)
    detiv[ie == 0] = 0.
    detiv = gaussian_filter(detiv, tim.psf_sigma)

    if apodize:
        apodize = int(apodize)
        ramp = np.arctan(np.linspace(-np.pi, np.pi, apodize+2))
        ramp = (ramp - ramp.min()) / (ramp.max()-ramp.min())
        # drop first and last (= 0 and 1)
        ramp = ramp[1:-1]
        detiv[:len(ramp),:] *= ramp[:,np.newaxis]
        detiv[:,:len(ramp)] *= ramp[np.newaxis,:]
        detiv[-len(ramp):,:] *= ramp[::-1][:,np.newaxis]
        detiv[:,-len(ramp):] *= ramp[::-1][np.newaxis,:]

    (Yo,Xo,Yi,Xi) = R
    if tim.dq is None:
        sat = None
    else:
        sat = ((tim.dq[Yi,Xi] & tim.dq_saturation_bits) > 0)
    return Yo, Xo, detim[Yi,Xi], detiv[Yi,Xi], sat

def detection_maps(tims, targetwcs, bands, mp, apodize=None):
    # Render the detection maps
    H,W = targetwcs.shape
    H,W = np.int(H), np.int(W)
    ibands = dict([(b,i) for i,b in enumerate(bands)])

    detmaps = [np.zeros((H,W), np.float32) for b in bands]
    detivs  = [np.zeros((H,W), np.float32) for b in bands]
    satmaps = [np.zeros((H,W), bool)       for b in bands]
    for tim, (Yo,Xo,incmap,inciv,sat) in zip(
        tims, mp.map(_detmap, [(tim, targetwcs, H, W, apodize) for tim in tims])):
        if Yo is None:
            continue
        ib = ibands[tim.band]
        detmaps[ib][Yo,Xo] += incmap * inciv
        detivs [ib][Yo,Xo] += inciv
        if sat is not None:
            satmaps[ib][Yo,Xo] |= sat
    for detmap,detiv in zip(detmaps, detivs):
        detmap /= np.maximum(1e-16, detiv)
    return detmaps, detivs, satmaps

def sed_matched_filters(bands):
    '''
    Determines which SED-matched filters to run based on the available
    bands.

    Returns
    -------
    SEDs : list of (name, sed) tuples
    '''
    # single-band filters
    SEDs = []
    for i,band in enumerate(bands):
        sed = np.zeros(len(bands))
        sed[i] = 1.
        SEDs.append((band, sed))
    # Reverse the order -- run z-band detection filter *first*.
    SEDs = list(reversed(SEDs))

    if len(bands) > 1:
        flat = dict(g=1., r=1., i=1., z=1.)
        SEDs.append(('Flat', [flat[b] for b in bands]))
        red = dict(g=2.5, r=1., i=0.4, z=0.4)
        SEDs.append(('Red', [red[b] for b in bands]))

    print('SED-matched filters:', SEDs)

    return SEDs

def run_sed_matched_filters(SEDs, bands, detmaps, detivs, omit_xy,
                            targetwcs, nsigma=5, saturated_pix=None,
                            exclusion_radius=4.,
                            veto_map=None,
                            plots=False, ps=None, mp=None):
    '''
    Runs a given set of SED-matched filters.

    Parameters
    ----------
    SEDs : list of (name, sed) tuples
        The SEDs to run.  The `sed` values are lists the same length
        as `bands`.
    bands : list of string
        The band names of `detmaps` and `detivs`.
    detmaps : numpy array, float
        Detection maps for each of the listed `bands`.
    detivs : numpy array, float
        Inverse-variances of the `detmaps`.
    omit_xy : None, or (xx,yy,rr) tuple
        Existing sources to avoid: x, y, radius.
    targetwcs : WCS object
        WCS object to use to convert pixel values into RA,Decs for the
        returned Tractor PointSource objects.
    nsigma : float, optional
        Detection threshold
    saturated_pix : None or list of numpy arrays, booleans
        Passed through to sed_matched_detection.
        A map of pixels that are always considered "hot" when
        determining whether a new source touches hot pixels of an
        existing source.
    exclusion_radius: int
        How many pixels around an existing source to veto
    plots : boolean, optional
        Create plots?
    ps : PlotSequence object
        Create plots?
    mp : multiproc object
        Multiprocessing
        
    Returns
    -------
    Tnew : fits_table
        Table of new sources detected
    newcat : list of PointSource objects
        Newly detected objects, with positions and fluxes, as Tractor
        PointSource objects.
    hot : numpy array of bool
        "Hot pixels" containing sources.
        
    See also
    --------
    sed_matched_detection : run a single SED-matched filter.
    
    '''
    from astrometry.util.fits import fits_table
    from tractor import PointSource, RaDecPos, NanoMaggies

    if omit_xy is not None:
        xx,yy,rr = omit_xy
        n0 = len(xx)
    else:
        xx,yy,rr = [],[],[]
        n0 = 0

    H,W = detmaps[0].shape
    hot = np.zeros((H,W), bool)

    peaksn = []
    apsn = []

    for sedname,sed in SEDs:
        #print('SED', sedname)
        if plots:
            pps = ps
        else:
            pps = None
        #t0 = Time()
        sedhot,px,py,peakval,apval = sed_matched_detection(
            sedname, sed, detmaps, detivs, bands, xx, yy, rr,
            nsigma=nsigma, saturated_pix=saturated_pix, veto_map=veto_map,
            ps=pps)
        #print('SED took', Time()-t0)
        if sedhot is None:
            continue
        print('SED', sedname, ':', len(px), 'new peaks')
        hot |= sedhot
        # With an empty xx, np.append turns it into a double!
        xx = np.append(xx, px).astype(int)
        yy = np.append(yy, py).astype(int)
        rr = np.append(rr, np.zeros_like(px) + exclusion_radius).astype(int)
        peaksn.extend(peakval)
        apsn.extend(apval)

    # New peaks:
    peakx = xx[n0:]
    peaky = yy[n0:]
    if len(peakx) == 0:
        return None,None,None

    # Add sources for the new peaks we found
    pr,pd = targetwcs.pixelxy2radec(peakx+1, peaky+1)
    print('Adding', len(pr), 'new sources')
    # Also create FITS table for new sources
    Tnew = fits_table()
    Tnew.ra  = pr
    Tnew.dec = pd
    Tnew.ibx = peakx
    Tnew.iby = peaky
    assert(len(peaksn) == len(Tnew))
    assert(len(apsn) == len(Tnew))
    Tnew.peaksn = np.array(peaksn)
    Tnew.apsn = np.array(apsn)
    newcat = []
    for r,d,x,y in zip(pr,pd,peakx,peaky):
        fluxes = dict([(band, detmap[y, x])
                       for band,detmap in zip(bands,detmaps)])
        newcat.append(PointSource(RaDecPos(r,d),
                                  NanoMaggies(order=bands, **fluxes)))
    return Tnew, newcat, hot

def plot_boundary_map(X, rgb=(0,255,0), extent=None, iterations=1):
    from scipy.ndimage.morphology import binary_dilation

    H,W = X.shape

    it = iterations
    #bounds = np.logical_xor(binary_dilation(X), X)
    padded = np.zeros((H+2*it, W+2*it), bool)
    padded[it:-it, it:-it] = X.astype(bool)
    bounds = np.logical_xor(binary_dilation(padded), padded)
    #rgba = np.zeros((H,W,4), np.uint8)
    rgba = np.zeros((H+2*it,W+2*it,4), np.uint8)
    rgba[:,:,0] = bounds*rgb[0]
    rgba[:,:,1] = bounds*rgb[1]
    rgba[:,:,2] = bounds*rgb[2]
    rgba[:,:,3] = bounds*255
    if extent is None:
        extent = [-it, W+it, -it, H+it]
    else:
        x0,x1,y0,y1 = extent
        extent = [x0-it, x1+it, y0-it, y1+it]

    plt.imshow(rgba, interpolation='nearest', origin='lower', extent=extent)<commit_msg>return 'hot' map even if no sources detected?<commit_after>from __future__ import print_function
import pylab as plt
import numpy as np
from astrometry.util.ttime import Time

def _detmap(X):
    from scipy.ndimage.filters import gaussian_filter
    from legacypipe.survey import tim_get_resamp
    (tim, targetwcs, H, W, apodize) = X
    R = tim_get_resamp(tim, targetwcs)
    if R is None:
        return None,None,None,None,None
    assert(tim.psf_sigma > 0)
    psfnorm = 1./(2. * np.sqrt(np.pi) * tim.psf_sigma)
    ie = tim.getInvvar()
    detim = tim.getImage().copy()
    tim.getSky().addTo(detim, scale=-1.)
    detim = gaussian_filter(detim, tim.psf_sigma) / psfnorm**2
    detsig1 = tim.sig1 / psfnorm
    subh,subw = tim.shape
    detiv = np.zeros((subh,subw), np.float32) + (1. / detsig1**2)
    detiv[ie == 0] = 0.
    detiv = gaussian_filter(detiv, tim.psf_sigma)

    if apodize:
        apodize = int(apodize)
        ramp = np.arctan(np.linspace(-np.pi, np.pi, apodize+2))
        ramp = (ramp - ramp.min()) / (ramp.max()-ramp.min())
        # drop first and last (= 0 and 1)
        ramp = ramp[1:-1]
        detiv[:len(ramp),:] *= ramp[:,np.newaxis]
        detiv[:,:len(ramp)] *= ramp[np.newaxis,:]
        detiv[-len(ramp):,:] *= ramp[::-1][:,np.newaxis]
        detiv[:,-len(ramp):] *= ramp[::-1][np.newaxis,:]

    (Yo,Xo,Yi,Xi) = R
    if tim.dq is None:
        sat = None
    else:
        sat = ((tim.dq[Yi,Xi] & tim.dq_saturation_bits) > 0)
    return Yo, Xo, detim[Yi,Xi], detiv[Yi,Xi], sat

def detection_maps(tims, targetwcs, bands, mp, apodize=None):
    # Render the detection maps
    H,W = targetwcs.shape
    H,W = np.int(H), np.int(W)
    ibands = dict([(b,i) for i,b in enumerate(bands)])

    detmaps = [np.zeros((H,W), np.float32) for b in bands]
    detivs  = [np.zeros((H,W), np.float32) for b in bands]
    satmaps = [np.zeros((H,W), bool)       for b in bands]
    for tim, (Yo,Xo,incmap,inciv,sat) in zip(
        tims, mp.map(_detmap, [(tim, targetwcs, H, W, apodize) for tim in tims])):
        if Yo is None:
            continue
        ib = ibands[tim.band]
        detmaps[ib][Yo,Xo] += incmap * inciv
        detivs [ib][Yo,Xo] += inciv
        if sat is not None:
            satmaps[ib][Yo,Xo] |= sat
    for detmap,detiv in zip(detmaps, detivs):
        detmap /= np.maximum(1e-16, detiv)
    return detmaps, detivs, satmaps

def sed_matched_filters(bands):
    '''
    Determines which SED-matched filters to run based on the available
    bands.

    Returns
    -------
    SEDs : list of (name, sed) tuples
    '''
    # single-band filters
    SEDs = []
    for i,band in enumerate(bands):
        sed = np.zeros(len(bands))
        sed[i] = 1.
        SEDs.append((band, sed))
    # Reverse the order -- run z-band detection filter *first*.
    SEDs = list(reversed(SEDs))

    if len(bands) > 1:
        flat = dict(g=1., r=1., i=1., z=1.)
        SEDs.append(('Flat', [flat[b] for b in bands]))
        red = dict(g=2.5, r=1., i=0.4, z=0.4)
        SEDs.append(('Red', [red[b] for b in bands]))

    print('SED-matched filters:', SEDs)

    return SEDs

def run_sed_matched_filters(SEDs, bands, detmaps, detivs, omit_xy,
                            targetwcs, nsigma=5, saturated_pix=None,
                            exclusion_radius=4.,
                            veto_map=None,
                            plots=False, ps=None, mp=None):
    '''
    Runs a given set of SED-matched filters.

    Parameters
    ----------
    SEDs : list of (name, sed) tuples
        The SEDs to run.  The `sed` values are lists the same length
        as `bands`.
    bands : list of string
        The band names of `detmaps` and `detivs`.
    detmaps : numpy array, float
        Detection maps for each of the listed `bands`.
    detivs : numpy array, float
        Inverse-variances of the `detmaps`.
    omit_xy : None, or (xx,yy,rr) tuple
        Existing sources to avoid: x, y, radius.
    targetwcs : WCS object
        WCS object to use to convert pixel values into RA,Decs for the
        returned Tractor PointSource objects.
    nsigma : float, optional
        Detection threshold
    saturated_pix : None or list of numpy arrays, booleans
        Passed through to sed_matched_detection.
        A map of pixels that are always considered "hot" when
        determining whether a new source touches hot pixels of an
        existing source.
    exclusion_radius: int
        How many pixels around an existing source to veto
    plots : boolean, optional
        Create plots?
    ps : PlotSequence object
        Create plots?
    mp : multiproc object
        Multiprocessing
        
    Returns
    -------
    Tnew : fits_table
        Table of new sources detected
    newcat : list of PointSource objects
        Newly detected objects, with positions and fluxes, as Tractor
        PointSource objects.
    hot : numpy array of bool
        "Hot pixels" containing sources.
        
    See also
    --------
    sed_matched_detection : run a single SED-matched filter.
    
    '''
    from astrometry.util.fits import fits_table
    from tractor import PointSource, RaDecPos, NanoMaggies

    if omit_xy is not None:
        xx,yy,rr = omit_xy
        n0 = len(xx)
    else:
        xx,yy,rr = [],[],[]
        n0 = 0

    H,W = detmaps[0].shape
    hot = np.zeros((H,W), bool)

    peaksn = []
    apsn = []

    for sedname,sed in SEDs:
        #print('SED', sedname)
        if plots:
            pps = ps
        else:
            pps = None
        #t0 = Time()
        sedhot,px,py,peakval,apval = sed_matched_detection(
            sedname, sed, detmaps, detivs, bands, xx, yy, rr,
            nsigma=nsigma, saturated_pix=saturated_pix, veto_map=veto_map,
            ps=pps)
        #print('SED took', Time()-t0)
        if sedhot is None:
            continue
        print('SED', sedname, ':', len(px), 'new peaks')
        hot |= sedhot
        # With an empty xx, np.append turns it into a double!
        xx = np.append(xx, px).astype(int)
        yy = np.append(yy, py).astype(int)
        rr = np.append(rr, np.zeros_like(px) + exclusion_radius).astype(int)
        peaksn.extend(peakval)
        apsn.extend(apval)

    # New peaks:
    peakx = xx[n0:]
    peaky = yy[n0:]

    Tnew = None
    newcat = []
    if len(peakx):
        # Add sources for the new peaks we found
        pr,pd = targetwcs.pixelxy2radec(peakx+1, peaky+1)
        print('Adding', len(pr), 'new sources')
        # Also create FITS table for new sources
        Tnew = fits_table()
        Tnew.ra  = pr
        Tnew.dec = pd
        Tnew.ibx = peakx
        Tnew.iby = peaky
        assert(len(peaksn) == len(Tnew))
        assert(len(apsn) == len(Tnew))
        Tnew.peaksn = np.array(peaksn)
        Tnew.apsn = np.array(apsn)
        for r,d,x,y in zip(pr,pd,peakx,peaky):
            fluxes = dict([(band, detmap[y, x])
                           for band,detmap in zip(bands,detmaps)])
            newcat.append(PointSource(RaDecPos(r,d),
                                      NanoMaggies(order=bands, **fluxes)))
    return Tnew, newcat, hot

def plot_boundary_map(X, rgb=(0,255,0), extent=None, iterations=1):
    from scipy.ndimage.morphology import binary_dilation

    H,W = X.shape

    it = iterations
    #bounds = np.logical_xor(binary_dilation(X), X)
    padded = np.zeros((H+2*it, W+2*it), bool)
    padded[it:-it, it:-it] = X.astype(bool)
    bounds = np.logical_xor(binary_dilation(padded), padded)
    #rgba = np.zeros((H,W,4), np.uint8)
    rgba = np.zeros((H+2*it,W+2*it,4), np.uint8)
    rgba[:,:,0] = bounds*rgb[0]
    rgba[:,:,1] = bounds*rgb[1]
    rgba[:,:,2] = bounds*rgb[2]
    rgba[:,:,3] = bounds*255
    if extent is None:
        extent = [-it, W+it, -it, H+it]
    else:
        x0,x1,y0,y1 = extent
        extent = [x0-it, x1+it, y0-it, y1+it]

    plt.imshow(rgba, interpolation='nearest', origin='lower', extent=extent)
{'fraction_non_alphanumeric': 0.13266535607838165, 'fraction_numerical': 0.011097723381317775, 'mean_word_length': 2.2375282282898787, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2330230', 'n_tokens_mistral': 5839, 'n_tokens_neox': 5623, 'n_words': 1447}