text:
import numpy as np

# Dictionary of object flags
objFlagDict = {'TOO_FEW_OBS':2**0,
               'BAD_COLOR':2**1,
               'VARIABLE':2**2,
               'TEMPORARY_BAD_STAR':2**3,
               'RESERVED':2**4,
               'REFSTAR_OUTLIER': 2**5,
               'BAD_QUANTITY': 2**6}

# Dictionary of observation flags
obsFlagDict = {'NO_EXPOSURE': 2**0,
               'BAD_ERROR': 2**1,
               'SUPERSTAR_OUTLIER': 2**2,
               'NO_ZEROPOINT': 2**4,
               'BAD_MAG': 2**5,
               'BAD_AIRMASS': 2**6,
               'FOCALPLANE_OUTLIER': 2**7,
               'FOCALPLANE_OUTLIER_REF': 2**8}

# Dictionary of exposure flags
expFlagDict = {'TOO_FEW_STARS':2**0,
               'EXP_GRAY_TOO_NEGATIVE':2**1,
               'VAR_GRAY_TOO_LARGE':2**2,
               'TOO_FEW_EXP_ON_NIGHT':2**3,
               'NO_STARS':2**4,
               'BAND_NOT_IN_LUT':2**5,
               'TEMPORARY_BAD_EXPOSURE':2**6,
               'EXP_GRAY_TOO_POSITIVE':2**7,
               'BAD_ZPFLAG':2**8}

# Dictionary of zeropoint flags
zpFlagDict = {'PHOTOMETRIC_FIT_EXPOSURE':2**0,
              'PHOTOMETRIC_NOTFIT_EXPOSURE':2**1,
              'NONPHOTOMETRIC_FIT_NIGHT':2**2,
              'NOFIT_NIGHT':2**3,
              'CANNOT_COMPUTE_ZEROPOINT':2**4,
              'TOO_FEW_STARS_ON_CCD':2**5}

# Dictionary of retrieval flags
retrievalFlagDict = {'EXPOSURE_RETRIEVED':2**0,
                     'EXPOSURE_INTERPOLATED':2**1,
                     'EXPOSURE_STANDARD':2**2}

logDict = {'NONE': 0,
           'WARN': 1,
           'INFO': 2,
           'DEBUG': 3}

class MaxFitIterations(Exception):
    """
    Raise when we've hit the maximum number of iterations
    """

    pass

def getMemoryString(location):
    """
    Get a string for memory usage (current and peak) for logging.

    parameters
    ----------
    location: string
       A short string which denotes where in the code the memory was recorded.
    """

    status = None
    result = {'peak':0, 'rss':0}
    memoryString = ''
    try:
        with open('/proc/self/status') as status:
            for line in status:
                parts = line.split()
                key = parts[0][2:-1].lower()
                if key in result:
                    result[key] = int(parts[1])/1000

            memoryString = 'Memory usage at %s: %d MB current; %d MB peak.' % (
                location, result['rss'], result['peak'])
    except:
        memoryString = 'Could not get process status for memory usage at %s!' % (location)

    return memoryString

def dataBinner(x,y,binSize,xRange,nTrial=100,xNorm=-1.0,minPerBin=5):
    """
    Bin data and compute errors via bootstrap resampling.  All median statistics.

    parameters
    ----------
    x: float array
       x values
    y: float array
       y values
    binSize: float
       Bin size
    xRange: float array [2]
       x limits for binning
    nTrial: float, optional, default=100
       Number of bootstraps
    xNorm: float, optional, default=-1
       Set the y value == 0 when x is equan to xNorm.  if -1.0 then no norm.
    minPerBin: int, optional, default=5
       Minimum number of points per bin

    returns
    -------
    binStruct: recarray, length nbins
       'X_BIN': Left edge of each bin
       'X': median x in bin
       'X_ERR': width of x distribution in bin
       'X_ERR_MEAN': error on the mean of the x's in bin
       'Y': mean y in bin
       'Y_WIDTH': width of y distribution in bin
       'Y_ERR': error on the mean of the y's in the bin

    """

    import esutil

    hist,rev=esutil.stat.histogram(x,binsize=binSize,min=xRange[0],max=xRange[1]-0.0001,rev=True)
    binStruct=np.zeros(hist.size,dtype=[('X_BIN','f4'),
                                        ('X','f4'),
                                        ('X_ERR_MEAN','f4'),
                                        ('X_ERR','f4'),
                                        ('Y','f4'),
                                        ('Y_WIDTH','f4'),
                                        ('Y_ERR','f4'),
                                        ('N','i4')])
    binStruct['X_BIN'] = np.linspace(xRange[0],xRange[1],hist.size)

    for i in range(hist.size):
        if (hist[i] >= minPerBin):
            i1a=rev[rev[i]:rev[i+1]]

            binStruct['N'][i] = i1a.size

            medYs=np.zeros(nTrial,dtype='f8')
            medYWidths=np.zeros(nTrial,dtype='f8')
            medXs=np.zeros(nTrial,dtype='f8')
            medXWidths=np.zeros(nTrial,dtype='f8')

            for t in range(nTrial):
                r=(np.random.random(i1a.size)*i1a.size).astype('i4')

                medYs[t] = np.median(y[i1a[r]])
                medYWidths[t] = 1.4826*np.median(np.abs(y[i1a[r]] - medYs[t]))

                medXs[t] = np.median(x[i1a[r]])
                medXWidths[t] = 1.4826*np.median(np.abs(x[i1a[r]] - medXs[t]))

            binStruct['X'][i] = np.median(medXs)
            binStruct['X_ERR'][i] = np.median(medXWidths)
            binStruct['X_ERR_MEAN'][i] = 1.4826*np.median(np.abs(medXs - binStruct['X'][i]))
            binStruct['Y'][i] = np.median(medYs)
            binStruct['Y_WIDTH'][i] = np.median(medYWidths)
            binStruct['Y_ERR'][i] = 1.4826*np.median(np.abs(medYs - binStruct['Y'][i]))

    if (xNorm >= 0.0) :
        ind=np.clip(np.searchsorted(binStruct['X_BIN'],xnorm),0,binStruct.size-1)
        binStruct['Y'] = binStruct['Y'] - binStruct['Y'][ind]

    return binStruct

def gaussFunction(x, *p):
    """
    A gaussian function for fitting

    parameters
    ----------
    x: float array
       x values
    *p: Tuple (A, mu, sigma)
    """

    A, mu, sigma = p
    return A*np.exp(-(x-mu)**2./(2.*sigma**2))

def histoGauss(ax,array):
    """
    Plot a histogram and fit a Gaussian to it.  Modeled after IDL histogauss.pro.

    parameters
    ----------
    ax: Plot axis object
       If None, return coefficients but do not plot
    array: float array to plot

    returns
    -------
    coeff: tuple (A, mu, sigma, fail)
       Fail is 0 if the fit succeeded
    """

    import scipy.optimize
    import matplotlib.pyplot as plt
    import esutil

    if array.size < 5:
        return np.array([1.0, 0.0, 1.0])

    q13 = np.percentile(array,[25,75])
    binsize=2*(q13[1] - q13[0])*array.size**(-1./3.)

    hist=esutil.stat.histogram(array,binsize=binsize,more=True)

    p0=[array.size,
        np.median(array),
        np.std(array)]

    try:
        with np.warnings.catch_warnings():
            np.warnings.simplefilter("ignore")

            # This fit might throw a warning, which we don't need now.
            # Note that in the future if we use the output from this fit in an
            # automated way we'll need to tell the parent that we didn't converge
            coeff, varMatrix = scipy.optimize.curve_fit(gaussFunction, hist['center'],
                                                        hist['hist'], p0=p0)

        coeff = np.append(coeff, 0)
    except:
        # set to starting values...
        coeff = np.append(p0, 1)

    hcenter=hist['center']
    hhist=hist['hist']

    rangeLo = coeff[1] - 5*coeff[2]
    rangeHi = coeff[1] + 5*coeff[2]

    lo,=np.where(hcenter < rangeLo)
    ok,=np.where(hcenter > rangeLo)
    hhist[ok[0]] += np.sum(hhist[lo])

    hi,=np.where(hcenter > rangeHi)
    ok,=np.where(hcenter < rangeHi)
    hhist[ok[-1]] += np.sum(hhist[hi])

    if ax is not None:
        ax.plot(hcenter[ok],hhist[ok],'b-',linewidth=3)
        ax.set_xlim(rangeLo,rangeHi)

        if coeff[3] == 0:
            xvals=np.linspace(rangeLo,rangeHi,1000)
            yvals=gaussFunction(xvals,*coeff[: -1])

            ax.plot(xvals,yvals,'k--',linewidth=3)
        ax.locator_params(axis='x',nbins=6)  # hmmm

    return coeff

"""
def plotCCDMap(ax, ccdOffsets, values, cbLabel, loHi=None):

    import matplotlib.pyplot as plt
    import matplotlib.patches as patches
    import matplotlib.colors as colors
    import matplotlib.cm as cmx

    cm = plt.get_cmap('rainbow')

    plotRARange = [ccdOffsets['DELTA_RA'].min() - ccdOffsets['RA_SIZE'].max()/2.,
                   ccdOffsets['DELTA_RA'].max() + ccdOffsets['RA_SIZE'].max()/2.]
    plotDecRange = [ccdOffsets['DELTA_DEC'].min() - ccdOffsets['DEC_SIZE'].max()/2.,
                    ccdOffsets['DELTA_DEC'].max() + ccdOffsets['DEC_SIZE'].max()/2.]

    if (loHi is None):
        st=np.argsort(values)

        lo=values[st[int(0.02*st.size)]]
        hi=values[st[int(0.98*st.size)]]
    else:
        lo = loHi[0]
        hi = loHi[1]

    cNorm = colors.Normalize(vmin=lo, vmax=hi)
    scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)

    Z=[[0,0],[0,0]]
    levels=np.linspace(lo,hi,num=150)
    CS3=plt.contourf(Z,levels,cmap=cm)

    ax.clear()

    ax.set_xlim(plotRARange[0]-0.05,plotRARange[1]+0.05)
    ax.set_ylim(plotDecRange[0]-0.05,plotDecRange[1]+0.05)
    ax.set_xlabel(r'$\delta\,\mathrm{R.A.}$',fontsize=16)
    ax.set_ylabel(r'$\delta\,\mathrm{Dec.}$',fontsize=16)
    ax.tick_params(axis='both',which='major',labelsize=14)

    for k in range(values.size):
        off=[ccdOffsets['DELTA_RA'][k],
             ccdOffsets['DELTA_DEC'][k]]

        ax.add_patch(
            patches.Rectangle(
                (off[0]-ccdOffsets['RA_SIZE'][k]/2.,
                 off[1]-ccdOffsets['DEC_SIZE'][k]/2.),
                ccdOffsets['RA_SIZE'][k],
                ccdOffsets['DEC_SIZE'][k],
                edgecolor="none",
                facecolor=scalarMap.to_rgba(values[k]))
            )

    cb=None
    cb = plt.colorbar(CS3,ticks=np.linspace(lo,hi,5))

    cb.set_label('%s' % (cbLabel), fontsize=14)

    return None
"""

class Cheb2dField(object):
    """
    Chebyshev 2d Field class.
    """
    def __init__(self, xSize, ySize, pars):
        """
        Instantiate a Cheb2dField

        Parameters
        ----------
        xSize: `int`
           Size of bounding box in x direction
        ySize: `int`
           Size of bounding box in y direction
        pars: `np.array`
           Parameters may be 2d (order + 1, order + 1) or
           1d (order + 1) * (order + 1)
        """
        self.xSize = xSize
        self.ySize = ySize

        if len(pars.shape) == 1:
            # this is a 1-d flat pars
            self.order = int(np.sqrt(len(pars))) - 1
            self.pars = pars.reshape((self.order + 1, self.order + 1))
        else:
            # This is a 2-d pars
            self.order = pars.shape[0] - 1
            self.pars = pars

    @classmethod
    def fit(cls, xSize, ySize, order, x, y, value, valueErr=None, triangular=True):
        """
        Construct a Cheb2dField by fitting a field of x/y/value

        Parameters
        ----------
        xSize: `int`
           Size of bounding box in x direction
        ySize: `int`
           Size of bounding box in y direction
        order: `int`
           Chebyshev order of fit
        x: `np.array`
           Float array of x values
        y: `np.array`
           Float array of y values
        value: `np.array`
           Float array of dependent values to fit
        valueErr: `np.array`, optional
           Float array of dependent value errors to fit.
           Default is None (unweighted fit)
        triangular: `bool`, optional
           Fit should suppress high-order cross terms.  Default is True

        Returns
        -------
        cheb2dField: `fgcm.Cheb2dField`
           The Cheb2dField object
        """

        fit = np.zeros((order + 1) * (order + 1))

        if triangular:
            iind = np.repeat(np.arange(order + 1), order + 1)
            jind = np.tile(np.arange(order + 1), order + 1)
            lowInds, = np.where((iind + jind) <= order)
        else:
            lowInds = np.arange(fit.size)

        # We add a 0.5 here because of lsst stack compatibility
        xScaled = (x + 0.5 - xSize/2.) / (xSize / 2.)
        yScaled = (y + 0.5 - ySize/2.) / (ySize / 2.)

        V = np.polynomial.chebyshev.chebvander2d(yScaled, xScaled, [order, order])

        if triangular:
            V = V[:, lowInds]

        if valueErr is not None:
            w = 1./valueErr**2.
            Vprime = np.matmul(np.diag(w), V)
        else:
            w = np.ones(value.size)
            Vprime = V

        fit[lowInds] = np.matmul(np.matmul(np.linalg.inv(np.matmul(Vprime.T,
                                                                   Vprime)),
                                           Vprime.T), value * w)

        return cls(xSize, ySize, fit.reshape((order + 1, order + 1)))

    def evaluate(self, x, y, flatPars=None):
        """
        Evaluate the chebyshev field at a given position.  Optionally can
        substitute a set of fit parameters (used in curve fitter).

        Parameters
        ----------
        x: `np.array`
           Float array of x values
        y: `np.array`
           Float array of y values
        flatPars: 'list' or `np.array`, optional
           Replacement parameters, must have size (order + 1) * (order + 1)
           If None, use self.pars

        Returns
        -------
        values: `np.array`
           Float array of Chebyshev field evaluated at x, y
        """

        if flatPars is not None:
            c = np.array(flatPars).reshape(self.order + 1, self.order + 1)
        else:
            c = self.pars

        xScaled = (x + 0.5 - self.xSize/2.) / (self.xSize / 2.)
        yScaled = (y + 0.5 - self.ySize/2.) / (self.ySize / 2.)

        return np.polynomial.chebyshev.chebval2d(yScaled, xScaled, c)

    def evaluateCenter(self):
        """
        Evaluate the chebyshev field at the center.

        Returns
        -------
        value: `float`
           Float value of Chebyshev field evaluated at center.
        """

        return float(self.evaluate(self.xSize/2. - 0.5, self.ySize/2. - 0.5))

    def __call__(self, xy, *flatpars):
        """
        Evaluate for use in a fitter.

        Parameters
        ----------
        yx: `np.array`
           Numpy vstack (2, nvalues)
        *flatpars: `float`
           Chebyshev parameters.

        Returns
        -------
        values: `np.array`
           Value of function evaluated at x = xy[0, :], y = xy[1, :]
        """

        return self.evaluate(xy[0, :], xy[1, :], flatpars)

"""
def plotCCDMap2d(ax, ccdOffsets, parArray, cbLabel, loHi=None):
    import matplotlib.pyplot as plt
    import matplotlib.colors as colors
    import matplotlib.cm as cmx

    cm = plt.get_cmap('rainbow')
    plt.set_cmap('rainbow')

    plotRARange = [ccdOffsets['DELTA_RA'].min() - ccdOffsets['RA_SIZE'].max()/2.,
                   ccdOffsets['DELTA_RA'].max() + ccdOffsets['RA_SIZE'].max()/2.]
    plotDecRange = [ccdOffsets['DELTA_DEC'].min() - ccdOffsets['DEC_SIZE'].max()/2.,
                    ccdOffsets['DELTA_DEC'].max() + ccdOffsets['DEC_SIZE'].max()/2.]

    # compute central values...
    centralValues = np.zeros(ccdOffsets.size)

    for i in range(ccdOffsets.size):
        field = Cheb2dField(ccdOffsets['X_SIZE'][i], ccdOffsets['Y_SIZE'][i], parArray[i, :])
        centralValues[i] = -2.5 * np.log10(field.evaluateCenter()) * 1000.0

    if (loHi is None):
        st=np.argsort(centralValues)

        lo = centralValues[st[int(0.02*st.size)]]
        hi = centralValues[st[int(0.98*st.size)]]
    else:
        lo = loHi[0]
        hi = loHi[1]

    cNorm = colors.Normalize(vmin=lo, vmax=hi)
    scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)

    Z=[[0,0],[0,0]]
    levels=np.linspace(lo,hi,num=150)
    CS3=plt.contourf(Z,levels,cmap=cm)

    ax.clear()

    ax.set_xlim(plotRARange[0]-0.05,plotRARange[1]+0.05)
    ax.set_ylim(plotDecRange[0]-0.05,plotDecRange[1]+0.05)
    ax.set_xlabel(r'$\delta\,\mathrm{R.A.}$',fontsize=16)
    ax.set_ylabel(r'$\delta\,\mathrm{Dec.}$',fontsize=16)
    ax.tick_params(axis='both',which='major',labelsize=14)

    for k in range(ccdOffsets.size):
        xValues = np.linspace(0.0, ccdOffsets['X_SIZE'][i], 50)
        yValues = np.linspace(0.0, ccdOffsets['Y_SIZE'][i], 50)

        xGrid = np.repeat(xValues, yValues.size)
        yGrid = np.tile(yValues, xValues.size)

        field = Cheb2dField(ccdOffsets['X_SIZE'][i], ccdOffsets['Y_SIZE'][i], parArray[k, :])
        zGrid = -2.5 * np.log10(np.clip(field.evaluate(xGrid, yGrid), 0.1, None)) * 1000.0

        # This seems to be correct
        extent = [ccdOffsets['DELTA_RA'][k] -
                  ccdOffsets['RASIGN'][k]*ccdOffsets['RA_SIZE'][k]/2.,
                  ccdOffsets['DELTA_RA'][k] +
                  ccdOffsets['RASIGN'][k]*ccdOffsets['RA_SIZE'][k]/2.,
                  ccdOffsets['DELTA_DEC'][k] -
                  ccdOffsets['DECSIGN'][k]*ccdOffsets['DEC_SIZE'][k]/2.,
                  ccdOffsets['DELTA_DEC'][k] +
                  ccdOffsets['DECSIGN'][k]*ccdOffsets['DEC_SIZE'][k]/2.]

        zGridPlot = zGrid.reshape(xValues.size, yValues.size)
        if ccdOffsets['XRA'][k]:
            zGridPlot = zGridPlot.T

        plt.imshow(zGridPlot,
                   interpolation='bilinear',
                   origin='lower',
                   extent=extent,
                   norm=cNorm)

    cb=None
    cb = plt.colorbar(CS3,ticks=np.linspace(lo,hi,5))

    cb.set_label('%s' % (cbLabel), fontsize=14)

    return None
"""

def plotCCDMap(ax, deltaMapper, values, cbLabel, loHi=None):
    """
    Plot CCD map with single values for each CCD.

    Parameters
    ----------
    ax :
        plot axis object
    deltaMapper : `np.recarray`
        Delta x/y, ra/dec mapper.
    values : `np.ndarray`
        Values for each ccd.
    cbLabel :  `str`
        Color bar label
    loHi : tuple [2], optional
        (lo, hi) if set.  Otherwise, scaling is computed from data.
    """
    import matplotlib.pyplot as plt
    import matplotlib.patches as patches
    import matplotlib.colors as colors
    import matplotlib.cm as cmx

    cm = plt.get_cmap('rainbow')

    plotRaRange = [np.min(deltaMapper['delta_ra']) - 0.02,
                   np.max(deltaMapper['delta_ra']) + 0.02]
    plotDecRange = [np.min(deltaMapper['delta_dec']) - 0.02,
                    np.max(deltaMapper['delta_dec']) + 0.02]

    if loHi is None:
        st = np.argsort(values)

        lo = values[st[int(0.02*st.size)]]
        hi = values[st[int(0.98*st.size)]]
    else:
        lo = loHi[0]
        hi = loHi[1]

    cNorm = colors.Normalize(vmin=lo, vmax=hi)
    scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)

    Z = [[0, 0], [0, 0]]
    levels = np.linspace(lo, hi, num=150)
    CS3 = plt.contourf(Z, levels, cmap=cm)

    ax.clear()

    ax.set_xlim(plotRaRange[0], plotRaRange[1])
    ax.set_ylim(plotDecRange[0], plotDecRange[1])
    ax.set_xlabel(r'$\delta\,\mathrm{R.A.}$', fontsize=16)
    ax.set_ylabel(r'$\delta\,\mathrm{Dec.}$', fontsize=16)
    ax.tick_params(axis='both', which='major', labelsize=14)

    zGrid = np.zeros_like(deltaMapper['x'])
    for k in range(deltaMapper.size):
        zGrid[k, :] = values[k]

    plt.scatter(deltaMapper['delta_ra'].ravel(), deltaMapper['delta_dec'].ravel(),
                s=0.1,
                c=zGrid.ravel(),
                vmin=lo, vmax=hi)

    cb = None
    cb = plt.colorbar(CS3, ticks=np.linspace(lo, hi, 5))
    cb.set_label('%s' % (cbLabel), fontsize=14)

    return None

def plotCCDMap2d(ax, deltaMapper, parArray, cbLabel, loHi=None):
    """
    Plot CCD map with Chebyshev fits for each CCD.

    Parameters
    ----------
    ax :
        plot axis object
    deltaMapper : `np.recarray`
        Delta x/y, ra/dec mapper.
    parArray : `np.ndarray`
        Chebyshev parameter array (nCCD, nPars)
    cbLabel :  `str`
        Color bar label
    loHi : tuple [2], optional
        (lo, hi) if set.  Otherwise, scaling is computed from data.
    """
    import matplotlib.pyplot as plt
    import matplotlib.colors as colors
    import matplotlib.cm as cmx

    cm = plt.get_cmap('rainbow')
    plt.set_cmap('rainbow')

    plotRaRange = [np.min(deltaMapper['delta_ra']) - 0.02,
                   np.max(deltaMapper['delta_ra']) + 0.02]
    plotDecRange = [np.min(deltaMapper['delta_dec']) - 0.02,
                    np.max(deltaMapper['delta_dec']) + 0.02]

    centralValues = np.zeros(len(deltaMapper))
    for i in range(deltaMapper.size):
        field = Cheb2dField(deltaMapper['x_size'][i],
                            deltaMapper['y_size'][i],
                            parArray[i, :])
        centralValues[i] = -2.5*np.log10(field.evaluateCenter())*1000.0

    if loHi is None:
        st = np.argsort(centralValues)
        lo = centralValues[st[int(0.02*st.size)]]
        hi = centralValues[st[int(0.98*st.size)]]
    else:
        lo = loHi[0]
        hi = loHi[1]

    cNorm = colors.Normalize(vmin=lo, vmax=hi)
    scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=cm)

    Z = [[0, 0], [0, 0]]
    levels = np.linspace(lo, hi, num=150)
    CS3 = plt.contourf(Z, levels, cmap=cm)

    ax.clear()

    ax.set_xlim(plotRaRange[0], plotRaRange[1])
    ax.set_ylim(plotDecRange[0], plotDecRange[1])
    ax.set_xlabel(r'$\delta\,\mathrm{R.A.}$', fontsize=16)
    ax.set_ylabel(r'$\delta\,\mathrm{Dec.}$', fontsize=16)
    ax.tick_params(axis='both', which='major', labelsize=14)

    zGrid = np.zeros_like(deltaMapper['x'])
    for k in range(deltaMapper.size):
        field = Cheb2dField(deltaMapper['x_size'][k],
                            deltaMapper['y_size'][k],
                            parArray[k, :])
        zGrid[k, :] = -2.5*np.log10(np.clip(field.evaluate(deltaMapper['x'][k, :],
                                                           deltaMapper['y'][k, :]),
                                                           0.1,
                                                           None))*1000.0

    plt.scatter(deltaMapper['delta_ra'].ravel(), deltaMapper['delta_dec'].ravel(),
                s=0.1,
                c=zGrid.ravel(),
                vmin=lo, vmax=hi)

    cb = None
    cb = plt.colorbar(CS3, ticks=np.linspace(lo, hi, 5))
    cb.set_label('%s' % (cbLabel), fontsize=14)

    return None

def logFlaggedExposuresPerBand(log, fgcmPars, flagName, raiseAllBad=True):
    """
    Log how many exposures per band were flagged with a given flag.

    Parameters
    ----------
    log: fgcmLog
    fgcmPars: `fgcm.fgcmParameters`
    flagName: `str`
       Name of flag
    raiseAllBad: `bool`, optional
       Raise if a band has all bad.  Default is True.
    """

    for i, band in enumerate(fgcmPars.bands):
        if not fgcmPars.hasExposuresInBand[i]:
            continue
        inBand, = np.where(fgcmPars.expBandIndex == i)
        inBandBad, = np.where((fgcmPars.expFlag[inBand] &
                               expFlagDict[flagName]) > 0)
        log.info(' %s: %s band: %d of %d exposures' %
                 (flagName, band, inBandBad.size, inBand.size))
        if raiseAllBad and inBandBad.size == inBand.size:
            raise RuntimeError("FATAL: All observations in %s band have been cut with %s" % (band, flagName))

def checkFlaggedExposuresPerBand(log, fgcmPars):
    """
    Raise if any band has had all exposures removed.

    Parameters
    ----------
    log: fgcmLog
    fgcmPars: `fgcm.fgcmParameters`

    """

    for i, band in enumerate(fgcmPars.bands):
        if not fgcmPars.hasExposuresInBand[i]:
            continue
        inBand, = np.where(fgcmPars.expBandIndex == i)
        inBandBad, = np.where(fgcmPars.expFlag[inBand] > 0)
        if inBandBad.size >= (inBand.size - 1):
            raise RuntimeError("FATAL: All observations in %s band have been cut!")

def computeDeltaRA(a, b, dec=None, degrees=False):
    """
    Compute a - b for RAs.  Optionally scale by cos(dec).

    Parameters
    ----------
    a : `float` or `np.ndarray`
    b : `float` or `np.ndarray`
    dec : `float` or `np.ndarray`
    degrees : `bool`, optional
    """
    if degrees:
        rotVal = 360.0
    else:
        rotVal = 2.0*np.pi

    delta = (np.atleast_1d(a) - np.atleast_1d(b)) % rotVal
    test, = np.where(delta > rotVal/2.)
    delta[test] -= rotVal

    if dec is not None and degrees:
        delta *= np.cos(np.deg2rad(dec))
    elif dec is not None:
        delta *= np.cos(dec)

    return delta

class FocalPlaneProjectorFromOffsets(object):
    """Create a focalPlaneProjector that returns the formatted
    ccdOffsetsTable.

    This is for backwards compatibility and straight DECam processing.

    Parameters
    ----------
    ccdOffsets : `np.ndarray`
        CCD offset table
    """
    def __init__(self, ccdOffsets):
        # Convert from input to internal format.
        dtype = ccdOffsets.dtype.descr
        dtype.extend([('XRA', bool),
                      ('RASIGN', 'i2'),
                      ('DECSIGN', 'i2')])
        self.ccdOffsets = np.zeros(ccdOffsets.size, dtype=dtype)
        for name in ccdOffsets.dtype.names:
            self.ccdOffsets[name][:] = ccdOffsets[name][:]

        self._deltaMapper = None
        self._computedSigns = False

    def __call__(self, orientation, nstep=50):
        """
        Make a focal plane projection mapping.

        Parameters
        ----------
        orientation : `int`
            Camera "orientation".  Unused for FocalPlaneProjectorFromOffsets
            which does not support arbitrary orientations.
        nstep : `int`
            Number of steps in x/y per detector for the mapping.

        Returns
        -------
        projectionMapping : `np.darray`
            A projection mapping object with x, y, x_size, y_size,
            delta_ra_cent, delta_dec_cent, delta_ra, delta_dec for
            each detector id.
        """
        if self._deltaMapper is not None:
            # Return cached value
            return self._deltaMapper

        if not self._computedSigns:
            raise RuntimeError("Cannot make a projection mapping without first "
                               "running computeCCDOffsetSigns")

        deltaMapper = np.zeros(len(self.ccdOffsets), dtype=[('id', 'i4'),
                                                            ('x', 'f8', nstep**2),
                                                            ('y', 'f8', nstep**2),
                                                            ('x_size', 'i4'),
                                                            ('y_size', 'i4'),
                                                            ('delta_ra_cent', 'f8'),
                                                            ('delta_dec_cent', 'f8'),
                                                            ('delta_ra', 'f8', nstep**2),
                                                            ('delta_dec', 'f8', nstep**2)])
        deltaMapper['id'] = np.arange(len(self.ccdOffsets))

        for i in range(len(self.ccdOffsets)):
            xSize = int(self.ccdOffsets[i]['X_SIZE'])
            ySize = int(self.ccdOffsets[i]['Y_SIZE'])

            xValues = np.linspace(0.0, xSize, nstep)
            yValues = np.linspace(0.0, ySize, nstep)

            deltaMapper['x'][i, :] = np.repeat(xValues, yValues.size)
            deltaMapper['y'][i, :] = np.tile(yValues, xValues.size)
            deltaMapper['x_size'][i] = xSize
            deltaMapper['y_size'][i] = ySize

            deltaMapper['delta_ra_cent'][i] = self.ccdOffsets[i]['DELTA_RA']
            deltaMapper['delta_dec_cent'][i] = self.ccdOffsets[i]['DELTA_DEC']

            # And translate these into delta_ra, delta_dec
            raValues = np.linspace(self.ccdOffsets['DELTA_RA'][i] -
                                   self.ccdOffsets['RASIGN'][i]*self.ccdOffsets['RA_SIZE'][i]/2.,
                                   self.ccdOffsets['DELTA_RA'][i] +
                                   self.ccdOffsets['RASIGN'][i]*self.ccdOffsets['RA_SIZE'][i]/2.,
                                   nstep)
            decValues = np.linspace(self.ccdOffsets['DELTA_DEC'][i] -
                                    self.ccdOffsets['DECSIGN'][i]*self.ccdOffsets['DEC_SIZE'][i]/2.,
                                    self.ccdOffsets['DELTA_DEC'][i] +
                                    self.ccdOffsets['DECSIGN'][i]*self.ccdOffsets['DEC_SIZE'][i]/2.,
                                    nstep)

            if not self.ccdOffsets['XRA'][i]:
                # Swap axes
                deltaMapper['delta_ra'][i, :] = np.tile(raValues, decValues.size)
                deltaMapper['delta_dec'][i, :] = np.repeat(decValues, raValues.size)
            else:
                deltaMapper['delta_ra'][i, :] = np.repeat(raValues, decValues.size)
                deltaMapper['delta_dec'][i, :] = np.tile(decValues, raValues.size)
        self._deltaMapper = deltaMapper

        return deltaMapper

    def computeCCDOffsetSigns(self, fgcmStars):
        """Compute plotting signs for x/y to ra/dec conversions.

        This does not support rotations; you need a proper
        focalPlaneProjector for that (only supported via lsst dm).

        Parameters
        ----------
        fgcmStars : `fgcmStars`
        """
        import scipy.stats
        import esutil
        from .sharedNumpyMemManager import SharedNumpyMemManager as snmm

        obsObjIDIndex = snmm.getArray(fgcmStars.obsObjIDIndexHandle)
        obsCCDIndex = snmm.getArray(fgcmStars.obsCCDHandle) - fgcmStars.ccdStartIndex
        obsExpIndex = snmm.getArray(fgcmStars.obsExpIndexHandle)

        obsX = snmm.getArray(fgcmStars.obsXHandle)
        obsY = snmm.getArray(fgcmStars.obsYHandle)
        objRA = snmm.getArray(fgcmStars.objRAHandle)
        objDec = snmm.getArray(fgcmStars.objDecHandle)

        h, rev = esutil.stat.histogram(obsCCDIndex, rev=True)

        for i in range(h.size):
            if h[i] == 0: continue

            i1a = rev[rev[i]:rev[i+1]]

            cInd = obsCCDIndex[i1a[0]]

            if self.ccdOffsets['RASIGN'][cInd] == 0:
                # choose a good exposure to work with
                hTest, revTest = esutil.stat.histogram(obsExpIndex[i1a], rev=True)
                maxInd = np.argmax(hTest)
                testStars = revTest[revTest[maxInd]:revTest[maxInd + 1]]

                testRA = objRA[obsObjIDIndex[i1a[testStars]]]
                testDec = objDec[obsObjIDIndex[i1a[testStars]]]
                testX = obsX[i1a[testStars]]
                testY = obsY[i1a[testStars]]

                corrXRA,_ = scipy.stats.pearsonr(testX, testRA)
                corrYRA,_ = scipy.stats.pearsonr(testY, testRA)

                if (np.abs(corrXRA) > np.abs(corrYRA)):
                    self.ccdOffsets['XRA'][cInd] = True
                else:
                    self.ccdOffsets['XRA'][cInd] = False

                if self.ccdOffsets['XRA'][cInd]:
                    # x is correlated with RA
                    if corrXRA < 0:
                        self.ccdOffsets['RASIGN'][cInd] = -1
                    else:
                        self.ccdOffsets['RASIGN'][cInd] = 1

                    corrYDec,_ = scipy.stats.pearsonr(testY, testDec)
                    if corrYDec < 0:
                        self.ccdOffsets['DECSIGN'][cInd] = -1
                    else:
                        self.ccdOffsets['DECSIGN'][cInd] = 1
                else:
                    # y is correlated with RA
                    if corrYRA < 0:
                        self.ccdOffsets['RASIGN'][cInd] = -1
                    else:
                        self.ccdOffsets['RASIGN'][cInd] = 1

                    corrXDec,_ = scipy.stats.pearsonr(testX, testDec)
                    if corrXDec < 0:
                        self.ccdOffsets['DECSIGN'][cInd] = -1
                    else:
                        self.ccdOffsets['DECSIGN'][cInd] = 1

        self._computedSigns = True

max_stars_repo_name: erykoff/fgcm

max_issues_repo_name: erykoff/fgcm

max_forks_repo_name: lsst/fgcm

reason: import numpy,import scipy


text:
import numpy as np
import math

line = input().split()
x1 = int(line[0])
y1 = int(line[1])
x2 = int(line[2])
y2 = int(line[3])

a = x2-x1
b = y2-y1
r = math.sqrt(a*a + b*b)

x3 = x2 - b
y3 = y2 + a
x4 = x3 - a
y4 = y3 - b

print("{} {} {} {}".format(x3, y3, x4, y4))

ext: py

reason: import numpy


text:
#!~/anaconda3/bin/python3
# ******************************************************
# Author: Pengshuai Yang
# Last modified: 2021-01-28 14:33
# Email: yps18@mails.tsinghua.edu.cn
# Filename: linear_train.py
# Description:
#   train the linear classifier
# ******************************************************
import os
import numpy as np
import random
import argparse
from collections import defaultdict
import torch
from torch.utils.data import Dataset, DataLoader
from torch.optim.lr_scheduler import ReduceLROnPlateau
import torch.optim as optim
import torch.nn as nn

import sklearn
from sklearn.model_selection import KFold
import timm

from utils import EarlyStopping
from dataset import get_training_set, get_validating_set, get_testing_set

import sys
sys.path.append('..')
from model import Linear, Cs_co

#TRAIN_DATA = '../BioImage_2015_data/dataset/train/'
#TEST_DATA = '../BioImage_2015_data/dataset/init_test/'
TRAIN_DATA = '../../NCT_CRC_data/all_data/train/'
TEST_DATA = '../../NCT_CRC_data/all_data/test/'

def get_img_list(path, k=None):
    total_list = os.listdir(path)
    total_list.sort()
    random.seed(8)
    random.shuffle(total_list)
    if k is None:
        return total_list
    else:
        return total_list[:k]

def train_data_info(img_list):
    info = defaultdict(int)
    for img in img_list:
        label = img.split('-')[0]
        info[label] += 1
    print(info)

def train_test_embedding(method, backbone, total_train_list, cls, device):
    if method == 'cs':
        embedding_net = Cs_co(backbone, 1, return_embedding=True)
        embedding_state_dict = embedding_net.state_dict()
        pretrain_dict = torch.load('../../generative_or_discriminative/model/recon_pretrain/0_ckpt0.04339.pth')
        state_dict = {k:v for k,v in pretrain_dict.items() if
                      k in embedding_state_dict}
        embedding_net.load_state_dict(state_dict)
    elif method == 'cs-co':
        embedding_net = Cs_co(backbone, 1, return_embedding=True)
        embedding_state_dict = embedding_net.state_dict()
        pretrain_dict = torch.load('../checkpoint/csco/csco_co_Adam-step_10000_96_0.001_1e-06_1.0_100_0.06461.pth')
        #pretrain_dict = torch.load('../../generative_or_discriminative/model/fix_decoder_hp/recon_NODE_b+r_Adam-step_10000_0.001_1e-06_1_l2.pth')
        state_dict = {k:v for k,v in pretrain_dict.items() if
                      k in embedding_state_dict}
        embedding_net.load_state_dict(state_dict)
    elif method == 'unfix':
        embedding_net = Cs_co(backbone, 1, return_embedding=True)
        embedding_state_dict = embedding_net.state_dict()
        pretrain_dict = torch.load('./checkpoint/unfix_decoder/unfix_decoder_co_Adam-step_10000_96_0.001_1e-06_1.0_100_0.09165.pth')
        state_dict = {k:v for k,v in pretrain_dict.items() if
                      k in embedding_state_dict}
        embedding_net.load_state_dict(state_dict)
    elif method == 'no_svp':
        embedding_net = Cs_co(backbone, 1, return_embedding=True)
        embedding_state_dict = embedding_net.state_dict()
        pretrain_dict = torch.load('./checkpoint/no_svp/no_svp_co_Adam-step_10000_96_0.001_1e-06_1.0_100_0.04802.pth')
        state_dict = {k:v for k,v in pretrain_dict.items() if
                      k in embedding_state_dict}
        embedding_net.load_state_dict(state_dict)

    embedding_net = embedding_net.to(device)

    single_channel = True
    train_dataset = get_testing_set(TRAIN_DATA, cls, single_channel, total_train_list)
    test_dataset = get_testing_set(TEST_DATA, cls, single_channel)

    train_loader = DataLoader(train_dataset, batch_size=128,
                              shuffle=False, drop_last=False, num_workers=4)
    test_loader = DataLoader(test_dataset, batch_size=128, shuffle=False,
                             drop_last=False, num_workers=4)

    train_embedding, train_label = compute_embedding(train_loader,
                                                     embedding_net, device,
                                                     single_channel)
    test_embedding, test_label = compute_embedding(test_loader, embedding_net,
                                                   device, single_channel)

    print(train_embedding.shape, train_label.shape)
    print(test_embedding.shape, test_label.shape)

    del embedding_net
    torch.cuda.empty_cache()
    return train_embedding, train_label, test_embedding, test_label

def compute_embedding(data_loader, embedding_net, device, single_channel):
    with torch.no_grad():
        embedding_net.eval()
        if single_channel:
            for batch_idx, (h, e, y) in enumerate(data_loader):
                h = h.to(device)
                e = e.to(device)

                tmp_embeddings = embedding_net(h,e)

                tmp_y_true = y.numpy()
                tmp_embeddings = tmp_embeddings.detach().cpu().numpy()

                if batch_idx == 0:
                    y_true = tmp_y_true
                    all_embeddings = tmp_embeddings
                else:
                    y_true = np.concatenate([y_true, tmp_y_true])
                    all_embeddings = np.concatenate([all_embeddings, tmp_embeddings])
        else:
            for batch_idx, (patches, y) in enumerate(data_loader):
                patches = patches.to(device)

                tmp_embeddings = embedding_net(patches)

                tmp_y_true = y.numpy()
                tmp_embeddings = tmp_embeddings.detach().cpu().numpy()

                if batch_idx == 0:
                    y_true = tmp_y_true
                    all_embeddings = tmp_embeddings
                else:
                    y_true = np.concatenate([y_true, tmp_y_true])
                    all_embeddings = np.concatenate([all_embeddings, tmp_embeddings])
    return all_embeddings, y_true

def train_epoch(train_loader, model, optimizer, loss_fn, device, epoch):
    model.train()
    losses = []
    total_loss = 0

    for batch_idx, (embed, y) in enumerate(train_loader):

        embed = embed.to(device)
        y = y.to(device)

        optimizer.zero_grad()
        y_prob = model(embed)

        loss_outputs = loss_fn(y_prob, y)
        losses.append(loss_outputs.item())
        total_loss += loss_outputs.item()

        if loss_outputs.requires_grad is True:
            loss_outputs.backward()
            optimizer.step()

        tmp_y_true = y.detach().cpu().numpy()
        tmp_y_pred = torch.argmax(y_prob, dim=1).detach().cpu().numpy()

        if batch_idx==0:
            y_true = tmp_y_true
            y_pred = tmp_y_pred
        else:
            y_true = np.concatenate([y_true, tmp_y_true])
            y_pred = np.concatenate([y_pred, tmp_y_pred])

        if batch_idx % 50 == 0:
            message = 'Train: [{}/{} ({:.2f}%)]\tLoss: {:.6f}, ACC: {:.6f}'.format(int(batch_idx*y.size(0)), len(train_loader.dataset), 100.*batch_idx/len(train_loader), np.mean(losses), np.mean(y_true==y_pred))
            print(message)
            losses = []

    print('Epoch: {}, total_loss: {:.6f}, acc: {:.6f}'.format(epoch, total_loss/(batch_idx+1), np.mean(y_true==y_pred)))

def eval_epoch(eval_loader, model, loss_fn, device, epoch, eval_type,
               early_stopping=None):
    with torch.no_grad():
        model.eval()
        val_loss = 0
        for batch_idx, (embed, y) in enumerate(eval_loader):
            embed = embed.to(device)
            y = y.to(device)

            y_prob = model(embed)

            loss_outputs = loss_fn(y_prob, y)
            val_loss += loss_outputs.item()

            tmp_y_true = y.detach().cpu().numpy()
            tmp_y_pred = torch.argmax(y_prob, dim=1).detach().cpu().numpy()

            if batch_idx==0:
                y_true = tmp_y_true
                y_pred = tmp_y_pred
            else:
                y_true = np.concatenate([y_true, tmp_y_true])
                y_pred = np.concatenate([y_pred, tmp_y_pred])

    print('{} loss: {:.6f}, acc: {:.6f}'.format(eval_type, val_loss/(batch_idx+1), np.mean(y_true==y_pred)))
    if eval_type == 'valid':
        early_stopping(val_loss/(batch_idx+1), model)
        return val_loss/(batch_idx+1)
    elif eval_type == 'test':
        return np.mean(y_true==y_pred)

class Narray_Dataset(Dataset):
    def __init__(self, embeddings, labels):
        super(Narray_Dataset, self).__init__()
        self.embeddings = embeddings
        self.labels = labels

    def __getitem__(self, index):
        x = self.embeddings[index]
        y = self.labels[index]
        return torch.from_numpy(x), y

    def __len__(self):
        return len(self.labels)

if __name__ == '__main__':

    parser = argparse.ArgumentParser(description='script for linear evaluation.')
    parser.add_argument('-m', '--method', dest='method', type=str,
                        choices=['cs', 'cs-co', 'unfix', 'no_svp'], help='embedding method')
    parser.add_argument('-b', '--backbone', dest='backbone', type=str,
                        default='resnet18', choices=['resnet18', 'resnet50'],
                        help='backbone')
    parser.add_argument('-d', '--datasize', dest='train_data_size', type=int,
                        default=1000, help='train data size')
    # get embedding method and backbone type
    args = parser.parse_args()

    # hyper-parameters
    device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')
    batch_size = 32
    epochs = 500
    lr = 0.001
    train_data_size = args.train_data_size
    num_fold = 5

    # classes and index
    #cls = {'NOR': 0, 'BEN': 1, 'INS': 2, 'INV': 3}
    cls = {'ADI':0, 'DEB':1, 'LYM':2, 'MUC':3, 'MUS':4, 'NORM':5, 'STR':6, 'TUM':7}
    print(cls)

    total_train_list = np.array(get_img_list(TRAIN_DATA+'patches/',
                                             k=train_data_size))
    train_data_info(total_train_list)

    # compute data embeddings
    train_embedding, train_label, test_embedding, test_label = train_test_embedding(args.method, args.backbone, total_train_list, cls, device)
    embed_dim = train_embedding.shape[1]

    # linear classifier and cross-validation
    kf = KFold(n_splits=num_fold)
    results = []
    e_losses = []
    stop_epochs = []
    for i, (train_index, valid_index) in enumerate(kf.split(total_train_list)):
        model = Linear(embed_dim, len(cls))
        model = model.to(device)
        optimizer = optim.SGD(model.parameters(), lr=lr, momentum=0.9)
        #optimizer = optim.Adam(model.parameters(), lr=lr)

        train_dataset = Narray_Dataset(train_embedding[train_index],
                                       train_label[train_index])
        valid_dataset = Narray_Dataset(train_embedding[valid_index],
                                       train_label[valid_index])
        test_dataset = Narray_Dataset(test_embedding, test_label)

        train_loader = DataLoader(train_dataset, batch_size=batch_size,
                                  shuffle=True, drop_last=True, num_workers=4)
        valid_loader = DataLoader(valid_dataset, batch_size=batch_size,
                                  shuffle=False, drop_last=False, num_workers=4)
        test_loader = DataLoader(test_dataset, batch_size=batch_size,
                                 shuffle=False, drop_last=False, num_workers=4)

        loss_fn = nn.CrossEntropyLoss()
        scheduler = ReduceLROnPlateau(optimizer, 'min', verbose=True)
        ckpt_path = './model/{}_'.format(args.method)+str(i)+'.pth'
        early_stopping = EarlyStopping(patience=15, verbose=True, delta=0.001,
                                       path=ckpt_path)

        for epoch in range(1, epochs+1):
            train_epoch(train_loader, model, optimizer, loss_fn, device, epoch)
            val_loss = eval_epoch(valid_loader, model, loss_fn, device, epoch, 'valid',
                                  early_stopping)
            scheduler.step(val_loss)
            if early_stopping.early_stop:
                print('Early Stop!')
                break

        model.load_state_dict(torch.load(ckpt_path))
        test_acc = eval_epoch(test_loader, model, loss_fn, device, 0,
                              'test')
        results.append(test_acc)
        e_losses.append(-early_stopping.best_score)
        stop_epochs.append(epoch)

        #break
    results = np.array(results)
    print('acc: ',results)
    print('loss: ',e_losses)
    print('stop_epochs: ', stop_epochs)
    print('5-fold test_acc, mean: {:.6f}, std: {:.6f}\n'.format(results.mean(),
                                                                results.std()))

    f = open('ablation.txt', 'a')
    f.write('method: {}, train_data_size: {}\n'.format(args.method,
                                                       args.train_data_size))
    f.write('{}\n'.format(str(results)))
    f.write('test_acc, mean: {:.6f}, std: {:.6f}\n\n'.format(results.mean(),
                                                                results.std()))
    f.close()


text:
from torch.utils.data import DataLoader
from Datasets.utils import ToTensor, Compose, CropCenter, dataset_intrinsics, DownscaleFlow, plot_traj, visflow
from Datasets.tartanTrajFlowDataset import TrajFolderDataset
from Datasets.transformation import ses2poses_quat
from evaluator.tartanair_evaluator import TartanAirEvaluator
from TartanVO import TartanVO

import argparse
import numpy as np
import cv2
from os import mkdir
from os.path import isdir

def get_args():
    parser = argparse.ArgumentParser(description='HRL')

    parser.add_argument('--batch-size', type=int, default=1,
                        help='batch size (default: 1)')
    parser.add_argument('--worker-num', type=int, default=1,
                        help='data loader worker number (default: 1)')
    parser.add_argument('--image-width', type=int, default=640,
                        help='image width (default: 640)')
    parser.add_argument('--image-height', type=int, default=448,
                        help='image height (default: 448)')
    parser.add_argument('--model-name', default='',
                        help='name of pretrained model (default: "")')
    parser.add_argument('--euroc', action='store_true', default=False,
                        help='euroc test (default: False)')
    parser.add_argument('--kitti', action='store_true', default=False,
                        help='kitti test (default: False)')
    parser.add_argument('--kitti-intrinsics-file',  default='',
                        help='kitti intrinsics file calib.txt (default: )')
    parser.add_argument('--test-dir', default='',
                        help='test trajectory folder where the RGB images are (default: "")')
    parser.add_argument('--pose-file', default='',
                        help='test trajectory gt pose file, used for scale calculation, and visualization (default: "")')
    parser.add_argument('--save-flow', action='store_true', default=False,
                        help='save optical flow (default: False)')

    args = parser.parse_args()

    return args

if __name__ == '__main__':
    args = get_args()

    testvo = TartanVO(args.model_name)

    # load trajectory data from a folder
    datastr = 'tartanair'
    if args.kitti:
        datastr = 'kitti'
    elif args.euroc:
        datastr = 'euroc'
    else:
        datastr = 'tartanair'
    focalx, focaly, centerx, centery = dataset_intrinsics(datastr)
    if args.kitti_intrinsics_file.endswith('.txt') and datastr=='kitti':
        focalx, focaly, centerx, centery = load_kiiti_intrinsics(args.kitti_intrinsics_file)

    transform = Compose([CropCenter((args.image_height, args.image_width)), DownscaleFlow(), ToTensor()])

    testDataset = TrajFolderDataset(args.test_dir,  posefile = args.pose_file, transform=transform,
                                        focalx=focalx, focaly=focaly, centerx=centerx, centery=centery)
    testDataloader = DataLoader(testDataset, batch_size=args.batch_size,
                                        shuffle=False, num_workers=args.worker_num)
    testDataiter = iter(testDataloader)

    motionlist = []
    testname = datastr + '_' + args.model_name.split('.')[0]
    if args.save_flow:
        flowdir = 'results/'+testname+'_flow'
        if not isdir(flowdir):
            mkdir(flowdir)
        flowcount = 0
    while True:
        try:
            sample = testDataiter.next()
        except StopIteration:
            break

        motions, flow = testvo.test_batch(sample)
        motionlist.extend(motions)

        if args.save_flow:
            for k in range(flow.shape[0]):
                flowk = flow[k].transpose(1,2,0)
                np.save(flowdir+'/'+str(flowcount).zfill(6)+'.npy',flowk)
                flow_vis = visflow(flowk)
                cv2.imwrite(flowdir+'/'+str(flowcount).zfill(6)+'.png',flow_vis)
                flowcount += 1

    poselist = ses2poses_quat(np.array(motionlist))

    # calculate ATE, RPE, KITTI-RPE
    if args.pose_file.endswith('.txt'):
        evaluator = TartanAirEvaluator()
        results = evaluator.evaluate_one_trajectory(args.pose_file, poselist, scale=True, kittitype=(datastr=='kitti'))
        if datastr=='euroc':
            print("==> ATE: %.4f" %(results['ate_score']))
        else:
            print("==> ATE: %.4f,\t KITTI-R/t: %.4f, %.4f" %(results['ate_score'], results['kitti_score'][0], results['kitti_score'][1]))

        # save results and visualization
        plot_traj(results['gt_aligned'], results['est_aligned'], vis=False, savefigname='results/'+testname+'.png', title='ATE %.4f' %(results['ate_score']))
        np.savetxt('results/'+testname+'.txt',results['est_aligned'])
    else:
        np.savetxt('results/'+testname+'.txt',poselist)

max_stars_repo_name: DinoHub/tartanvo

max_issues_repo_name: Grandzxw/tartanvo

max_forks_repo_name: Grandzxw/tartanvo


text:
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: percent
#       format_version: '1.3'
#       jupytext_version: 1.3.2
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# %% [markdown]
# ### Bar plots global means with extra cases
#
#

# %%
from oas_erf.constants import get_plotpath
from oas_erf.util.naming_conventions.var_info import get_fancy_var_name
from oas_erf.util.practical_functions import make_folders
from oas_erf.util.slice_average import one_val_tab
# noinspection PyBroadException
from oas_erf.util.slice_average.one_val_tab import get_mean_std_by_type

import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
from oas_erf.data_info import simulation_types
import matplotlib as mpl
import numpy as np
from IPython import get_ipython

try:
    _ipython = get_ipython()
    _magic = _ipython.magic
    _magic('load_ext autoreload')
    _magic('autoreload 2')
except:
    pass

# %% [markdown]
# ## Filenames

# %%
from oas_erf.util.plot.colors import get_case_col

version = 'clean_elvoc'
plt_path = get_plotpath('one_value')

def create_filename(name):
    fn = f'{plt_path}_bar_{version}_{name}.'
    make_folders(fn)
    return fn

# %%
startyear = '0004-01'
endyear = '0005-12'
pmin = 850.  # minimum pressure level
avg_over_lev = True  # True#True#False#True
pressure_adjust = True  # Can only be false if avg_over_lev false. Plots particular hybrid sigma lev

# %% [markdown]
# ## Models and variables to include:

# %%
varl = ['N_AER', 'NCONC01', 'TGCLDCWP', 'CDNUMC', 'NCFT_Ghan', 'DIR_Ghan', 'LWDIR_Ghan', 'SWDIR_Ghan', 'SWCF_Ghan',
        'LWCF_Ghan', 'cb_SO4_NA', 'cb_SOA_NA', 'cb_NA', 'SOA_NA', 'SO4_NA',
        'ACTNL_incld', 'ACTREL_incld', 'SFisoprene', 'SFmonoterp']
varl_ex = ['FSNT', 'FSNT_DRF', 'FLNT', 'FLNT_DRF', 'FSNTCDRF']
varl = varl + varl_ex
cases_sec = [
    'NF1850_SECT_ctrl',  # _smax',
    'NF1850_aeroxid2014_SECT_ctrl',  # _smax',
    'NF1850_SECT_elvoc_smax',
    'NF1850_aeroxid2014_SECT_elvoc_smax',
    'NF1850_SECT_svoc_smax',
    'NF1850_aeroxid2014_SECT_svoc_smax',

]
cases_nsec = [
    'NF1850_noSECT_def',  # _smax',
    'NF1850_aeroxid2014_noSECT_def',  # _smax',
    'NF1850_aeroxid2014_noSECT_ox_ricc',
    'NF1850_noSECT_ox_ricc',

]

# %%

case_types = ['PI', 'PIaerPD']
model_types = ['OsloAeroSec', 'OsloAero$_{imp}$', 'OsloAero$_{def}$']
case_types = ['PI', 'PD']
model_types = ['OsloAeroSec', 'OsloAeroSec$_{elvoc}$', 'OsloAero$_{imp}$', 'OsloAero$_{def}$'][
              ::-1]
model_types2 = ['OsloAeroSec', 'OsloAeroSec$_{elvoc}$', 'OsloAeroSec$_{svoc}$', 'OsloAero$_{imp}$', 'OsloAero$_{def}$'][
               ::-1]

cdic = {key: get_case_col(key) for key in model_types2}  # , ['r','g','b'])}
# %% [markdown]
# ## Import data

# %%
# varl = ['N_AER']
df2, dic_vals = one_val_tab.get_tab_yearly_mean(varl,
                                                cases_sec + cases_nsec,
                                                startyear,
                                                endyear,
                                                pmin=pmin,
                                                pressure_adjust=pressure_adjust,
                                                average_over_lev=avg_over_lev,
                                                groupby='time.year',  # 'time',
                                                dims=None,
                                                area='Global',
                                                invert_dic=True
                                                )

# %%

# %%
print(model_types2)

# %% [markdown]
# ## Get difference

# %%
relative = False
dic_diff = simulation_types.get_diff_by_type(dic_vals, varl, case_types=['PI', 'PD'],
                                             relative=relative,
                                             mod_types=model_types2,
                                             ctrl='PI'
                                             )['PD-PI']

ls = []
for key in dic_diff.keys():
    print(key)
    _df = dic_diff[key]

    _df['case'] = key
    print(_df.keys())
    ls.append(_df.reset_index())

df_tot = pd.concat(ls)

# %%
trans_dic = {v: get_fancy_var_name(v) for v in varl}

rn_dic = {
    'NCRE$_{Ghan}$': 'ERF$_{aci}$',
    'SWCRE$_{Ghan}$': 'ERF$_{aci,SW}$',
    'LWCRE$_{Ghan}$': 'ERF$_{aci,LW}$',
    'DRE$_{Ghan}$': 'ERF$_{ari}$'

}

svarl = ['NCFT_Ghan', 'SWCF_Ghan', 'LWCF_Ghan', 'DIR_Ghan']

df1 = df_tot[[*svarl, 'case']]
df1 = df1.rename(trans_dic, axis=1)
df1 = df1.rename(rn_dic, axis=1)
df2 = pd.melt(df1, id_vars='case')

# %%
cols = [cdic[c] for c in df_tot['case'].unique()]

# %%
df2.groupby(['case', 'variable']).mean()

# %%

figsize = [5, 3]
figsize = [4.6, 4.6]
f, ax = plt.subplots(figsize=figsize, dpi=150)

pts = np.linspace(0, np.pi * 2, 24)
circ = np.c_[np.sin(pts) / 2, -np.cos(pts) / 2]
vert = np.r_[circ, circ[::-1] * .7]
open_circle = mpl.path.Path(vert)
sns.barplot(y='value', x='variable', hue='case', data=df2, errcolor='.55', ci=90, palette=cols, alpha=.6)
g = sns.stripplot(y='value', x='variable', hue='case', data=df2, dodge=True, marker=open_circle, palette=['none'] * 3,
                  jitter=.2, label='_nolegend_')  # , add_=False)#label='_nolabel_')

handles, labels = ax.get_legend_handles_labels()

# When creating the legend, only use the first two elements
# to effectively remove the last two.
l = plt.legend(handles[5:], labels[5:], frameon=False, loc='center right')

ax.axhline(0, linewidth=0.4, c='k')
sns.despine(f, bottom=True, trim=True, offset=10)
# ax.legend(frameon=False)
plt.tick_params(
    axis='x',  # changes apply to the x-axis
    which='both',  # both major and minor ticks are affected
    bottom=False,  # ticks along the bottom edge are off
    top=False,  # ticks along the top edge are off
    labeltop=True,
    labelbottom=False
)  # labels along the bottom edge are off
ax.set_ylabel('[Wm$^{-2}$]')
ax.set_xlabel('')

fn = create_filename('forcing')
plt.tight_layout()
f.savefig(fn + 'pdf', dpi=300)
print(fn)
plt.show()

print(fn)

# %%
figsize = [5, 3]
svarl = ['FSNT', 'FSNT_DRF', 'FLNT', 'FLNT_DRF', 'FSNTCDRF']  # '', 'SWCF_Ghan', 'LWCF_Ghan','DIR_Ghan']

mean, std, mean_nn, std_nn = get_mean_std_by_type(dic_vals, svarl,
                                                  case_types=case_types, model_types=model_types, ctrl=case_types[0])

fig, ax = plt.subplots(1, figsize=figsize)
mean_nn.plot.bar(alpha=0.5, ax=ax, color=[cdic[c] for c in pd.DataFrame(mean).columns],
                 yerr=std_nn)  # , colors={'OsloAeroSec':'b'})
ax.axhline(0, linewidth=0.4, c='k')
sns.despine(fig, bottom=True, trim=True, offset=10)
ax.legend(frameon=False)
plt.tick_params(
    axis='x',  # changes apply to the x-axis
    which='both',  # both major and minor ticks are affected
    bottom=False,  # ticks along the bottom edge are off
    top=False,  # ticks along the top edge are off
    labeltop=True,
    labelbottom=False
)  # labels along the bottom edge are off
ax.set_ylabel('PDaer-PI [Wm$^{-2}$]')
fn = create_filename('forcing')
plt.tight_layout()
# fig.savefig(fn+'pdf', dpi=300)
plt.show()
# %% [markdown]
# ### diff NCRE:

# %%

# %%
figsize = [5, 3]
svarl = ['N_AER', 'NCONC01', 'cb_SOA_NA', 'cb_SO4_NA', 'cb_NA']
rn_dic = {
    'c.b. SOA$_{NPF}$+SO4$_{NPF}$': 'c.b. SOA$_{NPF}$ \n+ SO4$_{NPF}$',
    'c.b. SOA$_{NPF}$': 'c.b. \n SOA$_{NPF}$',
    'c.b. SO4$_{NPF}$': 'c.b. \n SO4$_{NPF}$',
}
mean, std, mean_nn, std_nn = get_mean_std_by_type(dic_vals, svarl, relative=True, model_types=model_types,
                                                  case_types=case_types)
mean_nn = mean_nn.rename(rn_dic)
std_nn = std_nn.rename(rn_dic)
fig, ax = plt.subplots(1, figsize=figsize)
mean_nn.plot.bar(alpha=0.5, ax=ax, color=[cdic[c] for c in pd.DataFrame(mean).columns],
                 yerr=std_nn)  # , colors={'OsloAeroSec':'b'})
ax.axhline(0, linewidth=0.4, c='k')
sns.despine(fig, bottom=True, trim=True, offset=10)
ax.legend(frameon=False)
plt.tick_params(
    axis='x',  # changes apply to the x-axis
    which='both',  # both major and minor ticks are affected
    bottom=False,  # ticks along the bottom edge are off
    top=False,  # ticks along the top edge are off
    labeltop=True,
    labelbottom=False
)  # labels along the bottom edge are off

# ax.set_ylabel('(PDaer-PI)/PI [%]')
fn = create_filename('aerosol')
plt.tight_layout()
fig.savefig(fn + 'pdf', dpi=300)
plt.show()
print(fn)

# %%
figsize = [5, 3]
svarl = ['N_AER', 'NCONC01', 'cb_SOA_NA', 'cb_SO4_NA', 'cb_NA']
rn_dic = {
    'c.b. SOA$_{NPF}$+SO4$_{NPF}$': 'c.b. SOA$_{NPF}$ \n+ SO4$_{NPF}$',
    'c.b. SOA$_{NPF}$': 'c.b. \n SOA$_{NPF}$',
    'c.b. SO4$_{NPF}$': 'c.b. \n SO4$_{NPF}$',
}
mean, std, mean_nn, std_nn = get_mean_std_by_type(dic_vals, svarl, relative=True, model_types=model_types2,
                                                  case_types=case_types)
mean_nn = mean_nn.rename(rn_dic)
std_nn = std_nn.rename(rn_dic)
fig, ax = plt.subplots(1, figsize=figsize)
mean_nn.plot.bar(alpha=0.5, ax=ax, color=[cdic[c] for c in pd.DataFrame(mean).columns],
                 yerr=std_nn)  # , colors={'OsloAeroSec':'b'})
ax.axhline(0, linewidth=0.4, c='k')
sns.despine(fig, bottom=True, trim=True, offset=10)
ax.legend(frameon=False)
plt.tick_params(
    axis='x',  # changes apply to the x-axis
    which='both',  # both major and minor ticks are affected
    bottom=False,  # ticks along the bottom edge are off
    top=False,  # ticks along the top edge are off
    labeltop=True,
    labelbottom=False
)  # labels along the bottom edge are off

# ax.set_ylabel('(PDaer-PI)/PI [%]')
fn = create_filename('aerosol')
plt.tight_layout()
# fig.savefig(fn+'pdf', dpi=300)
plt.show()
print(fn)

# %%

svarl = ['CDNUMC', 'ACTREL_incld', 'TGCLDCWP']  # 'cb_SOA_NA', 'cb_SO4_NA']
relative = True
mean, std, mean_nn, std_nn = get_mean_std_by_type(dic_vals, svarl, case_types=case_types, model_types=model_types,
                                                  relative=relative)

fig, ax = plt.subplots(1, figsize=figsize)
mean_nn.plot.bar(alpha=0.5, ax=ax, color=[cdic[c] for c in pd.DataFrame(mean).columns],
                 yerr=std_nn)  # , colors={'OsloAeroSec':'b'})
ax.axhline(0, linewidth=0.4, c='k')
sns.despine(fig, bottom=True, trim=True, offset=10)
ax.legend(frameon=False)
plt.tick_params(
    axis='x',  # changes apply to the x-axis
    which='both',  # both major and minor ticks are affected
    bottom=False,  # ticks along the bottom edge are off
    top=False,  # ticks along the top edge are off
    labeltop=True,
    labelbottom=False
)  # labels along the bottom edge are off

ax.set_ylabel('(PDaer-PI)/PI [%]')
fn = create_filename('cld_props')
plt.tight_layout()
# fig.savefig(fn+'pdf', dpi=300)
plt.show()

# %%

svarl = ['SFmonoterp', 'SFisoprene']  # 'cb_SOA_NA', 'cb_SO4_NA']
relative = True
mean, std, mean_nn, std_nn = get_mean_std_by_type(dic_vals, svarl, case_types=case_types,
                                                  model_types=model_types,
                                                  relative=relative)

fig, ax = plt.subplots(1, figsize=figsize)
mean_nn.plot.bar(alpha=0.5, ax=ax, color=[cdic[c] for c in pd.DataFrame(mean).columns],
                 yerr=std_nn)  # , colors={'OsloAeroSec':'b'})
ax.axhline(0, linewidth=0.4, c='k')
sns.despine(fig, bottom=True, trim=True, offset=10)
ax.legend(frameon=False)
plt.tick_params(
    axis='x',  # changes apply to the x-axis
    which='both',  # both major and minor ticks are affected
    bottom=False,  # ticks along the bottom edge are off
    top=False,  # ticks along the top edge are off
    labeltop=True,
    labelbottom=False
)  # labels along the bottom edge are off

ax.set_ylabel('(PDaer-PI)/PI [%]')
fn = create_filename('bvoc')
plt.tight_layout()
# fig.savefig(fn+'pdf', dpi=300)
plt.show()

max_stars_repo_name: sarambl/OAS-ERF

max_issues_repo_name: sarambl/OAS-ERF

max_forks_repo_name: sarambl/OAS-ERF


text:
from functools import reduce
from numpy import random
print('------------------------------------------------------------------')
print(' ***** ¡CONTINUA LA AVENTURA PARA RECORRER LAS ISLAS! ****')
print('---------------/El viaje a Isla Torva/------------------\n')

 # Parte I-A: Crea un diccionario con los objetivos y los golpes que hay
 # que dar a Valrock
hit={ "brazaleteD": 8, "brazaleteI": 8, "cinturon": 14 }
l1,l2=[],[]
for key, valor in hit.items():
        l1.append(key)
        l2 .append(valor)
total_hits=reduce(lambda x, y: x+y, l2)
print('-------------Tripulación ataquemos al monstruo...!!!!-------------')
print('Total de golpes que habrá de darle al Valrock es: ', total_hits )
print('------------------------------------------------------------------')
 # Parte I-B: Crea un método llamado attack que por parámetro reciba el
# objetivo, target, en formato string y muestre por pantalla la frase:
# "Pyratilla ha atacado al target".
l1=sorted(l1)
def attack(l1):
  """
  Args:

  Salida:
  """

  target=l1
  return '¡Pyratilla ha atacado con el:'
hit= {"brazaleteD": 8, "brazaleteI": 8, "cinturon": 14}
msj="¡Oh, no! ¡Pyratilla ha fallado!"
x=0
for key in hit:
    n,m,dispF,dispV =0,0,0,0
    print(f'** El objetivo a destruir es el: {key} ')
    if key=="brazaleteD" or key=="brazaleteI" or key=="cinturon":
        while m < (hit[key]):
          ataq=random.randint(0,2)
          if (dispF+ dispV)==5:
             print("¡Un segundo! ¡Hay que recargar los 5 cañones!")
             dispF, dispV=0,0
          if ataq ==0:
            print(msj)
            dispF +=1
          else:
            m +=1
            print(attack(l1), key)
            dispV +=1
            if m==hit[key]:
              x +=1
              print(f'ATENCION:¡El {key} ha sido eliminado!\n')
              if x==3:
                print('"¡¡¡El Valrock ha sido eliminado!!!" \n')
print('-----------¡Juega con Pyratilla a palabras encadenadas------------')

# Parte II-a: Crear la función palabras_encadenadas() que te permita hacer
# una partida de este juego crea una función helper llamada valid_option()
# que compare la última letra de la palabra anterior con la primera de la
# nueva. Si # coinciden, devuelve True. De lo contrario, devuelve False.

from functools import reduce
def valid_option(w1, w2):
    if w1[-1] == w2[0]:
        return True
    else:
        return False
w1 = input("Introduce la 1ra. palabra: ")
w2 = input("Introduce la 2da. palabra: ")

def palabras_encadenadas(w1, w2):
    played_words = []
    if w1==w2:
        print('No puedes introducir la misma palabra \n')
    else:
        if valid_option(w1, w2):
            print('"Felcitaciones Ud ha ganado..!! \n')
            played_words.append(w1)
            played_words.append(w2)
        else:
            print('Lo sentimos, la palabra no es válida')
            played_words.append(w1)
            played_words.append(w2)
    return played_words
played_words=palabras_encadenadas(w1, w2)
print(f'Las palabras jugadas fueron:', palabras_encadenadas(w1, w2))

# Parte II-b: ¡Ayuda a Pyratilla a comprobar si alguna de las palabras
# en played_words supera las 7 letras utilizando funciones lambda y
# alguna de las funciones vistas en la sección 11.

print('La logitud de las palabras ingresadas son:')
aux=reduce(lambda x,y: (len(x), len(y)), played_words)
print(aux)
aux=list(filter(lambda x: x>7, aux))
print(f' El número de caracter mayor que 7 es:  {aux} ')
# Parte II-c:  Ayuda a Pyratilla a eliminar Isla Torvaa de unvisited_islands y
# añadirla al conjunto visited_islands de la entrada 08 de su diario.
print('-------------------Llegada a la Isla Torvaa---------------------')
print('Navegante Pyerce, muestre la lista de todas las Islas pendientes por visitar')
islands=set(['Golosa','Arrecife','Lejana','Torva','Verde'])
print(islands)
print('\n Borren la Isla Torva, que ya estamos desembarcando, muestre la lista de Islas visitadas')
visited_islands=set(['Alegre','Arrecife','Espesura','Torva' ])
unvisited_islands=islands.difference(visited_islands)
print(f'Lista de Islas visitadas: {visited_islands}\n')

reason: from numpy


text:
from fltk import Fl
from PyCommon.modules.GUI import hpSimpleViewer as hsv
from PyCommon.modules.Renderer import ysRenderer as yr
from PyCommon.modules.Resource import ysMotionLoader as yf
import numpy as np

from random import randrange

def main(filenames):
    motion = []
    for filename in filenames:
        motion.append(yf.readBvhFile(filename))

    # viewer settings
    rd_contact_positions = [None]
    rd_contact_forces = [None]
    viewer = hsv.hpSimpleViewer(title='Bvh Viewer', rect=(0, 0, 1200, 800), viewForceWnd=False)
    # viewer.doc.addRenderer('MotionModel', yr.DartRenderer(ppo.env.ref_world, (150,150,255), yr.POLYGON_FILL))
    for motion_idx in range(len(motion)):
        viewer.doc.addRenderer('motion'+str(motion_idx), yr.JointMotionRenderer(motion[motion_idx], (randrange(256), randrange(256), randrange(256))))

    def preCallback(frame):
        for _motion in motion:
            _motion.frame = min(frame, len(_motion)-1)

    viewer.setPostFrameCallback_Always(preCallback)
    viewer.setMaxFrame(max([len(_motion)-1 for _motion in motion]))
    viewer.startTimer(1./30.)
    viewer.show()

    Fl.run()

if __name__ == '__main__':
    # main(['wd2_1foot_contact_run2.bvh', 'wd2_1foot_contact_run2_edit.bvh'])
    main(['wd2_1foot_contact_run2_edit.bvh', 'wd2_fast_2foot_hop_edit.bvh', 'wd2_slow_2foot_hop_edit.bvh', 'wd2_long_broad_jump_edit.bvh', 'wd2_short_broad_jump_edit.bvh', 'wd2_n_kick_edit.bvh', 'wd2_boxing_round_round_girl1_edit.bvh'])
    # main(['segfoot_wd2_n_kick_edit.bvh'])

    # main('../../data/woody2/Motion/Balancing/wd2_1foot_contact_run2.bvh')
    # main('../../data/woody2/Motion/Balancing/wd2_fast_2foot_hop.bvh')
    # main('../../data/woody2/Motion/Balancing/wd2_slow_2foot_hop.bvh')
    # main('../../data/woody2/Motion/Balancing/wd2_long_broad_jump.bvh')
    # main('../../data/woody2/Motion/Balancing/wd2_short_broad_jump.bvh')
    # main('../../data/woody2/Motion/Balancing/wd2_n_kick.bvh')

    # main('../../data/woody2/Motion/samsung_boxing/round/boxing_round_round_girl1.bvh')[505:658]

max_stars_repo_name: hpgit/HumanFoot

max_issues_repo_name: hpgit/HumanFoot

max_forks_repo_name: hpgit/HumanFoot


text:
import tensorflow as tf
import numpy as np
import pdb

from model import REG
from preprocessing import Synth, GenMatrix
from os.path import join
from utils import search_wav
from sklearn.cross_validation import train_test_split

np.random.seed(1234567)

FLAGS = tf.flags.FLAGS
tf.flags.DEFINE_string('clean_dir', '/mnt/md1/user_jonlu/Github/DeepDenoisingAutoencoder/data/raw/clean',
                       'set clean data folder')
tf.flags.DEFINE_string('noise_dir', '/mnt/md1/user_jonlu/Github/DeepDenoisingAutoencoder/data/raw/noise',
                       'set noise data folder')
tf.flags.DEFINE_string('noisy_dir', '../data/noisy',
                       'set noisy(clean mixed with noise) data folder')
tf.flags.DEFINE_string('enhanced_dir', '../data/enhanced',
                       'set enhanced data folder')
tf.flags.DEFINE_string('training_files_dir', '../data/training_files',
                       'set training files folder')
tf.flags.DEFINE_string('tb_dir', '../model/tb_logs', 'save tensorboard logs')
tf.flags.DEFINE_string('saver_dir', '../model/saver', 'set folder for saver')
tf.flags.DEFINE_bool('TRAIN', True, 'train this model or not')
tf.flags.DEFINE_bool('TEST', True, 'test this model or not')
tf.flags.DEFINE_integer('n_cores', 20, 'set cpu cores')
tf.flags.DEFINE_integer('epochs', 100, 'epochs for training iterations')
tf.flags.DEFINE_integer('batch_size', 32, 'number of batch size')
tf.flags.DEFINE_float('learning_rate', 1e-3, 'learning rate')

def main():
    clean_dir = FLAGS.clean_dir
    noise_dir = FLAGS.noise_dir
    noisy_dir = FLAGS.noisy_dir
    enhanced_dir = FLAGS.enhanced_dir
    tb_dir = FLAGS.tb_dir
    saver_dir = FLAGS.saver_dir
    TRAIN = FLAGS.TRAIN
    TEST = FLAGS.TEST
    ncores = FLAGS.n_cores
    epochs = FLAGS.epochs
    batch_size = FLAGS.batch_size
    lr = FLAGS.learning_rate

    train_task = 'same_noise'  # set task name for noting your dataset

    # ===========================================================
    # ===========       Synthesize Noisy Data        ============
    # ===========================================================
    clean_file_list = search_wav(clean_dir)
    clean_train_list, clean_test_list = train_test_split(
        clean_file_list, test_size=0.2)
    noise_file_list = search_wav(noise_dir)[0:20]
    noise_train_list, noise_test_list = train_test_split(
        noise_file_list, test_size=0.2)
    noise_test_list = noise_train_list # test on the same noise

    print('--- Synthesize Training Noisy Data ---')
    train_noisy_dir = join(noisy_dir, 'train')
    sr_clean = 16000
    sr_noise = 44100
    snr_list = ['20dB', '10dB', '0dB']
    data_num = None  # set data_num to make training data numbers for different snr
    syn_train = Synth(clean_train_list, noise_train_list, sr_clean, sr_noise)
    syn_train.gen_noisy(snr_list, train_noisy_dir,
                        data_num=data_num, ADD_CLEAN=True, cpu_cores=ncores)
    print('--- Synthesize Testing Noisy Data ---')
    test_noisy_dir = join(noisy_dir, 'test')
    sr_clean = 16000
    sr_noise = 44100
    data_num = None # set data_num to make testing data numbers for different snr
    snr_list = ['15dB']
    syn_test = Synth(clean_test_list, noise_test_list, sr_clean, sr_noise)
    syn_test.gen_noisy(snr_list, test_noisy_dir,
                        data_num=data_num, ADD_CLEAN=True, cpu_cores=ncores)
    # ===========================================================
    # ===========       Create Training Matrix       ============
    # ===========================================================
    print('--- Generate Training Matrix ---')
    train_task = 'same_noise'  # set task name for noting your dataset
    training_files_dir = FLAGS.training_files_dir
    train_noisy_dir = join(noisy_dir, 'train')
    DEL_TRAIN_WAV = True
    gen_mat = GenMatrix(training_files_dir, train_task, train_noisy_dir)
    split_num = 50  # number of spliting files
    iter_num = 2  # set iter number to use multi-processing, cpu_cores = split_num/iter_num
    input_sequence = False  # set input data is sequence or not
    gen_mat.create_h5(split_num=split_num, iter_num=iter_num,
                      input_sequence=input_sequence,
                      DEL_TRAIN_WAV=DEL_TRAIN_WAV)

    # ===========================================================
    # ===========             Main Model             ============
    # ===========================================================
    print('--- Build Model ---')
    note = 'DDAE'
    date = '0720'
    split_num = 50
    training_files_dir = FLAGS.training_files_dir
    model = REG(tb_dir, saver_dir, train_task, date, gpu_num='3', note=note)
    model.build(init_learning_rate=1e-3, reuse=False)

    print('--- Train Model ---')
    model.train(training_files_dir, split_num, epochs, batch_size)

    print('--- Test Model ---')
    testing_data_dir = join(noisy_dir, 'test')
    result_dir = '../data/enhanced/{}_{}/'.format(note, date)
    num_test = 30 # Set this number to decide how many testing data you wanna use. (None => All)
    cpu_cores = 30
    test_saver = '{}_{}/{}/best_saver_{}'.format(
        saver_dir, note, date, train_task)
    model.test(testing_data_dir, result_dir,
               test_saver, cpu_cores, num_test)

if __name__ == '__main__':
    main()


text:
"""Viewer for aligned datasets"""

import glob
import os

from datetime import datetime

from matplotlib import pyplot as plt
from matplotlib.backends.backend_pdf import PdfPages
from matplotlib import markers

from astropy.io import fits
import numpy as np

from stwcs import wcsutil

from drizzlepac.hlautils import astrometric_utils as amutils

plt.ioff()

class Datasets:
    def __init__(self, parent_dir, num_levels=1):
        # determine how many datasets are available to review
        levels = os.path.sep.join('*' * num_levels)
        dirnames = os.path.join(parent_dir, levels)
        self.parent_dir = os.path.split(os.path.abspath(parent_dir))[1]
        self.dirnames = dirnames

        # Determine list of products to be reviewed
        drznames = glob.glob(os.path.join(dirnames, '*drz.fits'))
        drcnames = glob.glob(os.path.join(dirnames, '*drc.fits'))
        # Now, remove duplicates (keep drc, del drz)
        prodnames = [d if d.replace('drz.fits', 'drc.fits') not in drcnames else d.replace('drz.fits', 'drc.fits') for d in drznames]
        # Insure that pytest-temp directories are not included either
        self.prodnames = [d for d in prodnames if 'current' not in d]

        # Extract key information from products for statistical use
        #
        # Define values to be included in summary
        self.wcsnames = []
        self.instruments = []
        self.detectors = []
        self.targnames = []
        self.texptimes = []
        print("Processing {} products".format(len(self.prodnames)))
        for p in self.prodnames:
            with fits.open(p) as prod:
                phdu = prod[0].header
                sci = prod[1].header if len(prod) > 1 else phdu
                self.wcsnames.append(sci['wcsname'])

    def create_summary(self):
        font_size = 12
        first_page = plt.figure(figsize=(8.2, 11.69))
        first_page.clf()
        txt = 'Summary of Alignment Results from {}'.format(self.parent_dir)
        first_page.text(0.05, 0.95, txt, transform=first_page.transFigure, size=font_size, ha="center")
        txt = 'Total Datasets : {}'.format(len(self.prodnames))
        first_page.text(0.1, 0.85, txt, transform=first_page.transFigure, size=font_size, ha="left")

        # compute basic stats on final WCS solutions
        apost = len([w for w in self.wcsnames if 'FIT' in w])
        defwcs = len([w for w in self.wcsnames if '-' not in w])
        apri = len(self.wcsnames) - apost - defwcs
        relonly = len([w for w in self.wcsnames if 'NONE' in w])

        txt = 'Datasets with a posteriori WCS : {}'.format(apost)
        first_page.text(0.1, 0.8, txt, transform=first_page.transFigure, size=font_size, ha="left")
        txt = 'Datasets with a posteriori relative alignment : {}'.format(relonly)
        first_page.text(0.1, 0.75, txt, transform=first_page.transFigure, size=font_size, ha="left")
        txt = 'Datasets with a priori WCS : {}'.format(apri)
        first_page.text(0.1, 0.7, txt, transform=first_page.transFigure, size=font_size, ha="left")
        txt = 'Datasets with pipeline-default WCS : {}'.format(defwcs)
        first_page.text(0.1, 0.65, txt, transform=first_page.transFigure, size=font_size, ha="left")

        now = datetime.now()
        rtime = datetime.strftime(now, "b %d %Y  %H:%M:%S")
        txt = 'Report created at {}'.format(rtime)
        first_page.text(0.1, 0.1, txt, transform=first_page.transFigure, size=font_size, ha="center")

        return first_page

    def create(self, pdfname='multipage_pdf.pdf'):

        with PdfPages(pdfname) as pdf:
            first_page = self.create_summary()
            pdf.savefig(first_page)
            plt.close()

            # Now generate a separate page for each dataset
            for p, w in zip(self.prodnames, self.wcsnames):
                result = create_product_page(p, wcsname=w)
                pdf.savefig(result)
                plt.close()

def create_product_page(prodname, zoom_size=128, wcsname=""):

    # obtain image data to display
    with fits.open(prodname) as prod:
        data = prod[1].data
        phdr = prod[0].header
        targname = phdr['targname']
        inst = phdr['instrume']
        det = phdr['detector']
        texptime = phdr['texptime']
        inexp = phdr['d001data'].split('[')[0]
        wcstype = prod[1].header['wcstype']
        wcs = wcsutil.HSTWCS(prod, ext=1)
    center = (data.shape[0] // 2, data.shape[1] // 2)
    prod_path = os.path.split(prodname)[0]

    data = np.nan_to_num(data, 0.0)

    # Get GAIA catalog
    refcat = amutils.create_astrometric_catalog(prodname, existing_wcs=wcs,
                                                output=None)
    refx, refy = wcs.all_world2pix(refcat['RA'], refcat['DEC'], 0)
    # Remove points outside the full-size image area
    rx = []
    ry = []
    zx = []
    zy = []
    for x, y in zip(refx, refy):
        if 0 < x < data.shape[1] and 0 < y < data.shape[0]:
            rx.append(x)
            ry.append(y)
        if -zoom_size < x - center[1] < zoom_size and -zoom_size < y - center[0] < zoom_size:
            zx.append(x - center[1] + zoom_size)
            zy.append(y - center[0] + zoom_size)

    # Define subplot regions on page
    fig = plt.figure(constrained_layout=True, figsize=(7, 10))
    gs = fig.add_gridspec(ncols=4, nrows=5)

    # title plots
    img_title = "{} image of {} with WCSNAME={}".format(prodname, targname, wcsname)
    plt.title(img_title, loc='center', fontsize=8)

    # Define image display
    fig_img = fig.add_subplot(gs[:3, :])
    fig_zoom = fig.add_subplot(gs[3:, 2:])
    fig_summary = fig.add_subplot(gs[3:, :2])

    # Compute display range
    dmax = (data.max() // 10)
    dscaled = np.log10(np.clip(data, -0.9, dmax) + 1)
    # identify zoom region around center of data
    zoom = dscaled[center[0] - zoom_size:center[0] + zoom_size,
                   center[1] - zoom_size:center[1] + zoom_size]

    # display full image
    fig_img.imshow(dscaled, cmap='gray', origin='lower')
    # display zoomed section
    fig_zoom.imshow(zoom, cmap='gray', origin='lower')

    # define markerstyle
    mstyle = markers.MarkerStyle(marker='o')
    mstyle.set_fillstyle('none')
    # plot GAIA sources onto full size image
    fig_img.scatter(rx, ry, marker=mstyle, alpha=0.25, c='cyan', s=3)
    fig_zoom.scatter(zx, zy, marker=mstyle, alpha=0.25, c='cyan')

    # Print summary info
    fsize = 8
    pname = os.path.split(prodname)[1]
    fig_summary.text(0.01, 0.95, "Summary for {}".format(pname), fontsize=fsize)
    fig_summary.text(0.01, 0.9, "WCSNAME: {}".format(wcsname), fontsize=fsize)
    fig_summary.text(0.01, 0.85, "TARGET: {}".format(targname), fontsize=fsize)
    fig_summary.text(0.01, 0.8, "Instrument: {}/{}".format(inst, det), fontsize=fsize)
    fig_summary.text(0.01, 0.7, "Total Exptime: {}".format(texptime), fontsize=fsize)
    fig_summary.text(0.01, 0.65, "WCSTYPE: {}".format(wcstype), fontsize=fsize)
    fig_summary.text(0.01, 0.5, "# of GAIA sources: {}".format(len(rx)), fontsize=fsize)

    if 'FIT' in wcsname:
        # Look for FIT RMS and other stats from headerlet
        exp = fits.open(os.path.join(prod_path, inexp))
        for ext in exp:
            if 'extname' in ext.header and ext.header['extname'] == 'HDRLET' \
                and ext.header['wcsname'] == wcsname:
                hdrlet = ext.headerlet
                rms_ra = hdrlet[0].header['rms_ra']
                rms_dec = hdrlet[0].header['rms_dec']
                nmatch = hdrlet[0].header['nmatch']
                catalog = hdrlet[0].header['catalog']
                break
        exp.close()
        fig_summary.text(0.01, 0.4, "RMS: RA={:0.3}mas, DEC={:0.3}mas".format(rms_ra, rms_dec), fontsize=fsize)
        fig_summary.text(0.01, 0.35, "# matches: {}".format(nmatch), fontsize=fsize)
        fig_summary.text(0.01, 0.3, "Matched to {} catalog".format(catalog), fontsize=fsize)

    return fig

max_stars_repo_name: srodney/drizzlepac

max_issues_repo_name: srodney/drizzlepac

max_forks_repo_name: srodney/drizzlepac

reason: import numpy,from astropy