text:
import xquantum.blas as xqla
import xquantum.qdk as xqqdk
import xquantum.qiskit as xqqkit
from typing import NamedTuple, List, Mapping, Any
import numpy as np
from enum import Enum, IntEnum
import math, random
import pandas as pd

# Engine

class EngineType(Enum):
    BLAS = 0
    QDK = 1
    QISKIT = 2

    def __str__(self):
        return self.name

class Engine(NamedTuple):
    type: EngineType
    config: Mapping[str, Any]

    def __str__(self):
        return "{0} with config {1}".format(self.type, self.config)

# Network structure

class Pauli(IntEnum):
    X = 1
    Y = 2
    Z = 3

    def __repr__(self):
        return "{}".format(self.value)

class RotationBlock(NamedTuple):
    fixedQubitIndex: int
    rotationQubitIndex: int
    rotationParameterIndex: int
    rotationAxis: Pauli

    def __str__(self):
        return "RB({0},{1},{2},{3})".format(self.fixedQubitIndex, self.rotationQubitIndex, self.rotationParameterIndex,
                                            self.rotationAxis.name)

class CircuitModel(NamedTuple):
    network: List[RotationBlock]
    measureQubitIndex: int
    # measureParameterIndex : int
    biasParameterIndex: int

    def __str__(self):
        return "\tNetwork: {0}\n\tmeasureQubit: {1}\n\tbiasParameter: {2}".format(
            ['{0}'.format(r) for r in self.network], self.measureQubitIndex, self.biasParameterIndex)

# Hyper Params and Outcome

class HyperParams(NamedTuple):
    epochs: int
    batchSize: int
    shots: int
    learningRate: float
    decay: float
    momentum: float
    padding: bool
    pad: float

    def __str__(self):
        return "HyperParams:\n\tepochs: {0}\n\tbatchSize: {1}\n\tshots: {2}\n\tlearningRate: {3:.3f}\n\tdecay: {4}\n\tmomentum: {5:.3f}\n\tpad: {6:.3f}".format(
            self.epochs, self.batchSize, self.shots, self.learningRate, self.decay, self.momentum,
            self.pad if self.padding else 0)

class Outcome(NamedTuple):
    label: int
    probability: float
    pi: float

    def __str__(self):
        return "Label: {0}  Pr({0}): {1:.5f}  pi: {2:.5f}".format(self.label, self.probability, self.pi)

# Functions

def trainModel(engine, circuitModel, hyperParams, params, trainingFeatureSet, labelSet):
    learningHistory = pd.DataFrame([], columns=['Epoch', 'Error', 'Learning rate', 'Params'])
    if len(trainingFeatureSet) != len(labelSet):
        print("Training and Label set do not have the same length!")
        return (params, 1, learningHistory)
    normalizedFeatureSet = _normalizeFeatureSet(trainingFeatureSet, hyperParams,
                                                _circuitModelRegisterSize(circuitModel))
    bestEpoch = 0
    learningRate = hyperParams.learningRate
    bestParams = params
    bestError = _errorEstimation(hyperParams, labelSet,
                                 _batchClassify(engine, circuitModel, hyperParams, params, normalizedFeatureSet))
    learningHistory = learningHistory.append(
        {
            'Epoch': 0,
            'Error': bestError,
            'Learning rate': learningRate,
            'Params': params
        }, ignore_index=True)
    batches = math.ceil(len(normalizedFeatureSet) / hyperParams.batchSize)
    indexes = list(range(len(normalizedFeatureSet)))
    epochParams = bestParams
    epochError = bestError
    batchParams = epochParams
    batchCorrections = [0] * len(batchParams)
    for r in range(hyperParams.epochs):
        random.shuffle(indexes)
        trainingFeatureSet = list(map(lambda i: normalizedFeatureSet[i], indexes))
        trainingLabelSet = list(map(lambda i: labelSet[i], indexes))
        print(
            "Starting epoch {0}/{1} - last epoch error {3:.3f} - best error {4:.3f} in epoch {5} - learningRate {2:.3f}"
            .format(r + 1, hyperParams.epochs, learningRate, epochError, bestError, bestEpoch))
        for b in range(batches):
            print("\tTraining model - batch {0}/{1} - params {2}".format(b + 1, batches,
                                                                         [np.around(x, 3) for x in batchParams]))
            batchStart = b * hyperParams.batchSize
            batchEnd = min((b + 1) * hyperParams.batchSize, len(normalizedFeatureSet))
            (batchParams, batchCorrections) = _paramsOptimizaiton(engine, circuitModel, hyperParams, learningRate,
                                                                  batchParams, batchCorrections,
                                                                  trainingFeatureSet[batchStart:batchEnd],
                                                                  trainingLabelSet[batchStart:batchEnd])
            # batchError = _errorEstimation(hyperParams, labelSet, _batchClassify(engine, circuitModel, hyperParams, batchParams, normalizedFeatureSet))
        epochParams = batchParams
        epochError = _errorEstimation(
            hyperParams, labelSet, _batchClassify(engine, circuitModel, hyperParams, epochParams, normalizedFeatureSet))
        learningHistory = learningHistory.append(
            {
                'Epoch': (r + 1),
                'Error': epochError,
                'Learning rate': learningRate,
                'Params': batchParams
            },
            ignore_index=True)
        learningRate = _annealing(hyperParams, learningRate, r, (bestError - epochError))
        if epochError < bestError:
            bestError = epochError
            bestParams = epochParams
            bestEpoch = r + 1
            print("\t** New best epoch - error {0:.3f} - params {1}".format(bestError,
                                                                            [np.around(x, 3) for x in bestParams]))
    return (bestParams, bestError, learningHistory)

def _paramsOptimizaiton(engine, circuitModel, hyperParams, learningRate, params, corrections, normalizedFeatureSet,
                        labelSet):
    # calculate gradient of error function by params
    outcomes = _batchClassify(engine, circuitModel, hyperParams, params, normalizedFeatureSet)
    mse = [2 * (p - l) for p, l in zip(map(lambda o: o.pi, outcomes), labelSet)]  # MSE derivative coefficients
    gradient = [0] * len(params)
    for i in range(len(normalizedFeatureSet)):
        derr = _circuitDerivativeByParams(engine, circuitModel, hyperParams, params, normalizedFeatureSet[i])
        for p in range(len(params)):
            gradient[p] += mse[i] * derr[p] / len(normalizedFeatureSet)
    # update params
    newCorrections = [
        hyperParams.momentum * corrections[i] + (1 - hyperParams.momentum) * learningRate * gradient[i]
        for i in range(len(params))
    ]
    newParams = [params[i] - newCorrections[i] for i in range(len(params))]
    return (newParams, newCorrections)

def _annealing(hyperParams, learningRate, epoch, errorTrend):
    return hyperParams.learningRate / (1 + hyperParams.decay * epoch)

def _errorEstimation(hyperParams, labels, outcomes):
    return sum([(p - l)**2 for l, p in zip(labels, map(lambda o: o.label, outcomes))]) / len(labels)

def testModel(engine, circuitModel, hyperParams, params, testFeatureSet):
    normalizedFeatureSet = _normalizeFeatureSet(testFeatureSet, hyperParams, _circuitModelRegisterSize(circuitModel))
    outcomes = []
    batches = math.ceil(len(normalizedFeatureSet) / hyperParams.batchSize)
    for b in range(batches):
        batchStart = b * hyperParams.batchSize
        batchEnd = min((b + 1) * hyperParams.batchSize, len(normalizedFeatureSet))
        print("\tTesting model - batch {0}/{1}".format(b + 1, batches))
        outcomes = outcomes + _testModel(engine, circuitModel, hyperParams, params,
                                         normalizedFeatureSet[batchStart:batchEnd])
    return outcomes

def classify(engine, circuitModel, hyperParams, params, features):
    normalizedFeatures = _normalizeFeatures(features, hyperParams, _circuitModelRegisterSize(circuitModel))
    return _classify(engine, circuitModel, hyperParams, params, normalizedFeatures)

def _batchClassify(engine, circuitModel, hyperParams, params, normalizedFeatureSet):
    return _testModel(engine, circuitModel, hyperParams, params, normalizedFeatureSet)

# Common functions

def _circuitModelRegisterSize(circuit):
    size = circuit.measureQubitIndex
    for un in circuit.network:
        size = max(size, un.fixedQubitIndex, un.rotationQubitIndex)
    return size + 1

def _normalizeFeatures(features, hyperParams, circuitSize):
    f = len(features) + (1 if hyperParams.padding else 0)
    n = max((f - 1).bit_length(), circuitSize)
    values = [0] * (2**n)
    norm = 0
    for i in range(0, f):
        values[i] = features[i] if i < len(features) else hyperParams.pad if hyperParams.padding else 0
        norm += values[i]**2
    norm = math.sqrt(norm)
    return list(map(lambda x: x / norm, values))

def _normalizeFeatureSet(featureSet, hyperParams, circuitSize):
    return list(map(lambda f: _normalizeFeatures(f, hyperParams, circuitSize), featureSet))

def _defineOutcome(pr, pi):
    return Outcome(1 if pi > 0.5 else 0, pr, pi)

# Redirect to specific Engine methods

def _testModel(engine, circuitModel, hyperParams, params, normalizedFeatureSet):
    if engine.type == EngineType.BLAS:
        return [
            _defineOutcome(pr, pi)
            for (pr, pi) in xqla.testModel(engine, circuitModel, hyperParams, params, normalizedFeatureSet)
        ]
    if engine.type == EngineType.QDK:
        return [
            _defineOutcome(pr, pi)
            for (pr, pi) in xqqdk.testModel(engine, circuitModel, hyperParams, params, normalizedFeatureSet)
        ]
    if engine.type == EngineType.QISKIT:
        return [
            _defineOutcome(pr, pi)
            for (pr, pi) in xqqkit.testModel(engine, circuitModel, hyperParams, params, normalizedFeatureSet)
        ]

def _classify(engine, circuitModel, hyperParams, params, normalizedFeatures):
    (pr, pi) = (0, 0)
    if engine.type == EngineType.BLAS:
        (pr, pi) = xqla.classifiedLableProbability(engine, circuitModel, hyperParams, params, normalizedFeatures)
    if engine.type == EngineType.QDK:
        (pr, pi) = xqqdk.classifiedLableProbability(engine, circuitModel, hyperParams, params, normalizedFeatures)
    if engine.type == EngineType.QISKIT:
        (pr, pi) = xqqkit.classifiedLableProbability(engine, circuitModel, hyperParams, params, normalizedFeatures)
    return _defineOutcome(pr, pi)

def _circuitDerivativeByParams(engine, circuitModel, hyperParams, params, normalizedFeatures):
    if engine.type == EngineType.BLAS:
        return xqla.circuitDerivativeByParams(engine, circuitModel, hyperParams, params, normalizedFeatures)
    if engine.type == EngineType.QDK:
        return xqqdk.circuitDerivativeByParams(engine, circuitModel, hyperParams, params, normalizedFeatures)
    if engine.type == EngineType.QISKIT:
        return xqqkit.circuitDerivativeByParams(engine, circuitModel, hyperParams, params, normalizedFeatures)

max_stars_repo_name: caccio/CrossQuantumClassifier

max_issues_repo_name: caccio/CrossQuantumClassifier

max_forks_repo_name: caccio/CrossQuantumClassifier


text:
'''
This script is used to produce lists of CCDs or bricks, for production
purposes (building qdo queue, eg).

python legacypipe/queue-calibs.py  | qdo load cal -

eg, DR3:

* Staging to $SCRATCH to images required to run the EDR region:

module switch legacysurvey/dr3
python legacypipe/queue-calibs.py --region edr --write-ccds edr-ccds.fits --calibs --touching --nper 100
for x in $(tablist edr-ccds.fits"[col image_filename]" | awk '{print $2}' | sort | uniq); do
  rsync -LRarv $LEGACY_SURVEY_DIR/images/./$(echo $x | sed s/o.i/o??/g) /global/cscratch1/sd/desiproc/legacypipe-dir/images/;
done

* Running calibration preprocessing for EDR:

module switch legacysurvey/dr3-cori-scratch
python legacypipe/queue-calibs.py --region edr --calibs --touching --nper 100
qdo load cal jobs

# qdo launch cal 1 --cores_per_worker 1 --batchqueue shared --script "python legacypipe/run-calib.py --splinesky" --walltime 4:00:00 --keep_env --batchopts "-a 0-15"

qdo launch cal 1 --cores_per_worker 1 --batchqueue shared --script $(pwd)/cal.sh --walltime 4:00:00 --keep_env --batchopts "-a 0-15"
cal.sh:
  cd /global/cscratch1/sd/desiproc/code/legacypipe/py
  python legacypipe/run-calib.py --splinesky $*

* Running bricks for EDR:

python legacypipe/queue-calibs.py --region edr --brickq 0 > bricks
grep '^....[pm]...$' bricks | qdo load edr0 -

qdo launch edr0 16 --cores_per_worker 8 --batchqueue regular --walltime 12:00:00 --script ../bin/pipebrick.sh --keep_env

Brick stage 1:

python legacypipe/queue-calibs.py --region edr --brickq 1 > bricks1
grep '^....[pm]...$' bricks1 | awk '{print "/global/cscratch1/sd/desiproc/code/legacypipe/bin/pipebrick.sh",$1}' | qdo load jobs -

* Staging larger regions:

python legacypipe/queue-calibs.py --region northwest --write-ccds ccds-nw.fits --calibs --near --nper 100 --command --opt "--threads 8"

for x in $(tablist ccds-nw.fits"[col image_filename]" | awk '{print $2}' | sort | uniq); do   rsync -LRarv /project/projectdirs/desiproc/dr3/images/./$(echo $x | sed s/o.i/o??/g) /global/cscratch1/sd/desiproc/dr3/images/; done

qdo load --priority -10 jobs jobs

'''
from __future__ import print_function
import sys
import os
import numpy as np

from astrometry.util.fits import fits_table

from legacypipe.survey import LegacySurveyData, wcs_for_brick, ccds_touching_wcs

from astrometry.libkd.spherematch import match_radec

def log(*s):
    print(' '.join([str(ss) for ss in s]), file=sys.stderr)

def main(args):
    """Main program.
    """
    import argparse

    parser = argparse.ArgumentParser(description="This script is used to produce lists of CCDs or bricks, for production purposes (building qdo queue, eg).")
    parser.add_argument('--calibs', action='store_true',
                      help='Output CCDs that need to be calibrated.')

    parser.add_argument('--nper', type=int, default=None,
                      help='Batch N calibs per line')
    parser.add_argument('--byexp', action='store_true', default=False,
                        help='Run one whole exposure per job (not one CCD per job)')

    parser.add_argument('--forced', action='store_true',
                      help='Output forced-photometry commands')

    parser.add_argument('--lsb', action='store_true',
                      help='Output Low-Surface-Brightness commands')

    parser.add_argument('--stage', help='Stage image files to given directory')

    parser.add_argument('--touching', action='store_true',
                      help='Cut to only CCDs touching selected bricks')
    parser.add_argument('--near', action='store_true',
                      help='Quick cut to only CCDs near selected bricks')

    parser.add_argument('--check-coadd', action='store_true',
                      help='Check which coadds actually need to run.')
    parser.add_argument('--out', help='Output filename for calibs, default %(default)s',
                      default='jobs')
    parser.add_argument('--command', action='store_true',
                      help='Write out full command-line to run calib')
    parser.add_argument('--opt', help='With --command, extra options to add')

    parser.add_argument('--maxra', type=float, help='Maximum RA to run')
    parser.add_argument('--minra', type=float, help='Minimum RA to run')
    parser.add_argument('--maxdec', type=float, help='Maximum Dec to run')
    parser.add_argument('--mindec', type=float, help='Minimum Dec to run')

    parser.add_argument('--region', help='Region to select')

    parser.add_argument('--bricks', help='Set bricks.fits file to load')
    parser.add_argument('--ccds', help='Set ccds.fits file to load')
    parser.add_argument('--ignore_cuts', action='store_true',default=False,help='no photometric cuts')
    parser.add_argument('--save_to_fits', action='store_true',default=False,help='save cut brick,ccd to fits table')
    parser.add_argument('--name', action='store',default='dr3',help='save with this suffix, e.g. refers to ccds table')

    parser.add_argument('--delete-sky', action='store_true',
                      help='Delete any existing sky calibration files')

    parser.add_argument('--write-ccds', help='Write CCDs list as FITS table?')

    parser.add_argument('--nccds', action='store_true', default=False, help='Prints number of CCDs per brick')

    parser.add_argument('--bands', default='g,r,z', help='Set bands to keep: comma-separated list.')

    opt = parser.parse_args(args)

    want_ccds = (opt.calibs or opt.forced or opt.lsb)
    want_bricks = not want_ccds

    survey = LegacySurveyData()
    if opt.bricks is not None:
        B = fits_table(opt.bricks)
        log('Read', len(B), 'from', opt.bricks)
    else:
        B = survey.get_bricks()

    log('Bricks Dec range:', B.dec.min(), B.dec.max())

    if opt.ccds is not None:
        T = fits_table(opt.ccds)
        log('Read', len(T), 'from', opt.ccds)
    else:
        T = survey.get_ccds()
        log(len(T), 'CCDs')
    T.index = np.arange(len(T))

    if opt.ignore_cuts == False:
        log('Applying CCD cuts...')
        if 'ccd_cuts' in T.columns():
            T.cut(T.ccd_cuts == 0)
            log(len(T), 'CCDs survive cuts')

    bands = opt.bands.split(',')
    log('Filters:', np.unique(T.filter))
    T.cut(np.flatnonzero(np.array([f in bands for f in T.filter])))
    log('Cut to', len(T), 'CCDs in filters', bands)

    log('CCDs Dec range:', T.dec.min(), T.dec.max())

    # I,J,d,counts = match_radec(B.ra, B.dec, T.ra, T.dec, 0.2, nearest=True, count=True)
    # plt.clf()
    # plt.hist(counts, counts.max()+1)
    # plt.savefig('bricks.png')
    # B.cut(I[counts >= 9])
    # plt.clf()
    # plt.plot(B.ra, B.dec, 'b.')
    # #plt.scatter(B.ra[I], B.dec[I], c=counts)
    # plt.savefig('bricks2.png')

    # DES Stripe82
    #rlo,rhi = 350.,360.
    # rlo,rhi = 300., 10.
    # dlo,dhi = -6., 4.
    # TINY bit
    #rlo,rhi = 350.,351.1
    #dlo,dhi = 0., 1.1

    # EDR+
    # 860 bricks
    # ~10,000 CCDs
    #rlo,rhi = 239,246
    #dlo,dhi =   5, 13

    # DR1
    #rlo,rhi = 0, 360
    # part 1
    #dlo,dhi = 25, 40
    # part 2
    #dlo,dhi = 20,25
    # part 3
    #dlo,dhi = 15,20
    # part 4
    #dlo,dhi = 10,15
    # part 5
    #dlo,dhi = 5,10
    # the rest
    #dlo,dhi = -11, 5
    #dlo,dhi = 15,25.5

    dlo,dhi = -25, 40
    rlo,rhi = 0, 360

    # Arjun says 3x3 coverage area is roughly
    # RA=240-252 DEC=6-12 (but not completely rectangular)

    # COSMOS
    #rlo,rhi = 148.9, 151.2
    #dlo,dhi = 0.9, 3.5

    # A nice well-behaved region (EDR2/3)
    # rlo,rhi = 243.6, 244.6
    # dlo,dhi = 8.1, 8.6

    # 56 bricks, ~725 CCDs
    #B.cut((B.ra > 240) * (B.ra < 242) * (B.dec > 5) * (B.dec < 7))
    # 240 bricks, ~3000 CCDs
    #B.cut((B.ra > 240) * (B.ra < 244) * (B.dec > 5) * (B.dec < 9))
    # 535 bricks, ~7000 CCDs
    #B.cut((B.ra > 240) * (B.ra < 245) * (B.dec > 5) * (B.dec < 12))

    if opt.region in ['test1', 'test2', 'test3', 'test4']:
        nm = dict(test1='2446p115', # weird stuff around bright star
                  test2='1183p292', # faint sources around bright galaxy
                  test3='3503p005', # DES
                  test4='1163p277', # Pollux
                  )[opt.region]

        B.cut(np.flatnonzero(np.array([s == nm for s in B.brickname])))
        log('Cut to', len(B), 'bricks')
        log(B.ra, B.dec)
        dlo,dhi = -90,90
        rlo,rhi = 0, 360

    elif opt.region == 'edr':
        # EDR:
        # 535 bricks, ~7000 CCDs
        rlo,rhi = 240,245
        dlo,dhi =   5, 12

    elif opt.region == 'dr8-decam':
        rlo,rhi =   0, 360
        dlo,dhi = -70,  40
        log('DR8-DECam region')

    elif opt.region == 'edrplus':
        rlo,rhi = 235,248
        dlo,dhi =   5, 15

    elif opt.region == 'edr-south':
        rlo,rhi = 240,245
        dlo,dhi =   5, 10

    elif opt.region == 'cosmos1':
        # 16 bricks in the core of the COSMOS field.
        rlo,rhi = 149.75, 150.75
        dlo,dhi = 1.6, 2.6

    elif opt.region == 'pristine':
        # Stream?
        rlo,rhi = 240,250
        dlo,dhi = 10,15

    elif opt.region == 'des':
        dlo, dhi = -6., 4.
        rlo, rhi = 317., 7.

        T.cut(np.flatnonzero(np.array(['CPDES82' in fn for fn in T.cpimage])))
        log('Cut to', len(T), 'CCDs with "CPDES82" in filename')

    elif opt.region == 'subdes':
        rlo,rhi = 320., 360.
        dlo,dhi = -1.25, 1.25

    elif opt.region == 'northwest':
        rlo,rhi = 240,360
        dlo,dhi = 20,40
    elif opt.region == 'north':
        rlo,rhi = 120,240
        dlo,dhi = 20,40
    elif opt.region == 'northeast':
        rlo,rhi = 0,120
        dlo,dhi = 20,40
    elif opt.region == 'southwest':
        rlo,rhi = 240,360
        dlo,dhi = -20,0
    elif opt.region == 'south':
        rlo,rhi = 120,240
        dlo,dhi = -20,0
    elif opt.region == 'southeast':
        rlo,rhi = 0,120
        dlo,dhi = -20,0
    elif opt.region == 'southsoutheast':
        rlo,rhi = 0,120
        dlo,dhi = -20,-10
    elif opt.region == 'midwest':
        rlo,rhi = 240,360
        dlo,dhi = 0,20
    elif opt.region == 'middle':
        rlo,rhi = 120,240
        dlo,dhi = 0,20
    elif opt.region == 'mideast':
        rlo,rhi = 0,120
        dlo,dhi = 0,20

    elif opt.region == 'grz':
        # Bricks with grz coverage.
        # Be sure to use  --bricks survey-bricks-in-dr1.fits
        # which has_[grz] columns.
        B.cut((B.has_g == 1) * (B.has_r == 1) * (B.has_z == 1))
        log('Cut to', len(B), 'bricks with grz coverage')

    elif opt.region == 'nogrz':
        # Bricks without grz coverage.
        # Be sure to use  --bricks survey-bricks-in-dr1.fits
        # which has_[grz] columns.
        B.cut(np.logical_not((B.has_g == 1) * (B.has_r == 1) * (B.has_z == 1)))
        log('Cut to', len(B), 'bricks withOUT grz coverage')

    elif opt.region == 'deep2':
        rlo,rhi = 250,260
        dlo,dhi = 30,35

    elif opt.region == 'deep2f2':
        rlo,rhi = 251.4, 254.4
        dlo,dhi =  34.6,  35.3

    elif opt.region == 'deep2f3':
        rlo,rhi = 351.25, 353.75
        dlo,dhi = 0, 0.5

    elif opt.region == 'deep3':
        rlo,rhi = 214,216
        dlo,dhi = 52.25,53.25

    elif opt.region == 'virgo':
        rlo,rhi = 185,190
        dlo,dhi =  10, 15

    elif opt.region == 'virgo2':
        rlo,rhi = 182,192
        dlo,dhi =   8, 18

    elif opt.region == 'coma':
        # van Dokkum et al Coma cluster ultra-diffuse galaxies: 3x3 field centered on Coma cluster
        rc,dc = 195., 28.
        dd = 1.5
        cosdec = np.cos(np.deg2rad(dc))
        rlo,rhi = rc - dd/cosdec, rc + dd/cosdec
        dlo,dhi = dc - dd, dc + dd

    elif opt.region == 'lsb':
        rlo,rhi = 147.2, 147.8
        dlo,dhi = -0.4, 0.4

    elif opt.region == 'eboss-sgc':
        # generous boundaries to make sure get all relevant images
        # RA -45 to +45
        # Dec -5 to +7
        rlo,rhi = 310., 50.
        dlo,dhi = -6., 6.

    elif opt.region == 'eboss-ngc':
        # generous boundaries to make sure get all relevant images
        # NGC ELGs
        # RA 115 to 175
        # Dec 15 to  30
        # rlo,rhi = 122., 177.
        # dlo,dhi =  12.,  32.
        rlo,rhi = 126., 168.
        dlo,dhi =  18.,  33.

    elif opt.region == 'mzls':
        dlo,dhi = -10., 90. # -10: pull in Stripe 82 data too

    elif opt.region == 'dr4-bootes':
        # https://desi.lbl.gov/trac/wiki/DecamLegacy/DR4sched
        #dlo,dhi = 34., 35.
        #rlo,rhi = 209.5, 210.5
        dlo,dhi = 33., 36.
        rlo,rhi = 216.5, 219.5

    elif opt.region == 'des-sn-x3':
        #rlo,rhi = 36., 37.
        #dlo,dhi = -5., -4.
        rlo,rhi = 36., 36.5
        dlo,dhi = -4.5, -4.

    elif opt.region == 'ngc2632':
        # open cluster
        rlo,rhi = 129.0, 131.0
        dlo,dhi = 19.0, 20.5

    elif opt.region == 'dr8sky':
        rlo,rhi = 35.0, 37.0
        dlo,dhi = -3.0, -1.0

    # ADM DR8 test regions, see, e.g.:
    # https://desi.lbl.gov/trac/wiki/DecamLegacy/DR8#Testregions
    elif opt.region == 'dr8-test-s82':
        rlo, rhi = 0, 45
        dlo, dhi = -1.25, 1.25
    elif opt.region == 'dr8-test-hsc-sgc':
        rlo, rhi = 30, 40
        dlo, dhi = -6.5, -1.25
    elif opt.region == 'dr8-test-hsc-ngc':
        rlo, rhi = 177.5, 182.5
        dlo, dhi = -1, 1
    elif opt.region == 'dr8-test-edr':
        rlo, rhi = 240, 245
        dlo, dhi = 5, 12
    elif opt.region == 'dr8-test-hsc-north':
        rlo, rhi = 240, 250
        dlo, dhi = 42, 45
    elif opt.region == 'dr8-test-deep2-egs':
        rlo, rhi = 213, 216.5
        dlo, dhi = 52, 54
    elif opt.region == 'dr8-test-overlap':
        rlo, rhi = 132, 140.5
        dlo, dhi = 31.5, 35

    if opt.mindec is not None:
        dlo = opt.mindec
    if opt.maxdec is not None:
        dhi = opt.maxdec
    if opt.minra is not None:
        rlo = opt.minra
    if opt.maxra is not None:
        rhi = opt.maxra

    if rlo < rhi:
        B.cut((B.ra >= rlo) * (B.ra <= rhi) *
              (B.dec >= dlo) * (B.dec <= dhi))
    else: # RA wrap
        B.cut(np.logical_or(B.ra >= rlo, B.ra <= rhi) *
              (B.dec >= dlo) * (B.dec <= dhi))
    log(len(B), 'bricks in range; cut Dec range', B.dec.min(), B.dec.max())
    #for name in B.get('brickname'):
    #    print(name)
    #B.writeto('bricks-cut.fits')

    bricksize = 0.25
    # A bit more than 0.25-degree brick radius + Bok image radius ~ 0.57
    search_radius = 1.05 * np.sqrt(2.) * (bricksize +
                                          (0.455 * 4096 / 3600.))/2.

    log(len(T), 'CCDs')
    log(len(B), 'Bricks')
    I,J,d = match_radec(B.ra, B.dec, T.ra, T.dec, search_radius,
                        nearest=True)
    B.cut(I)
    log('Cut to', len(B), 'bricks near CCDs')
    log('Bricks Dec range:', B.dec.min(), B.dec.max())

    # plt.clf()
    # plt.plot(B.ra, B.dec, 'b.')
    # plt.title('DR3 bricks')
    # plt.axis([360, 0, np.min(B.dec)-1, np.max(B.dec)+1])
    # plt.savefig('bricks.png')

    if opt.touching:
        I,J,d = match_radec(T.ra, T.dec, B.ra, B.dec, search_radius,
                            nearest=True)
        # list the ones that will be cut
        # drop = np.ones(len(T))
        # drop[I] = False
        # for i in np.flatnonzero(drop):
        #     from astrometry.util.starutil_numpy import degrees_between
        #     dists = degrees_between(B.ra, B.dec, T.ra[i], T.dec[i])
        #     mindist = min(dists)
        #     print('Dropping:', T.ra[i], T.dec[i], 'min dist', mindist, 'search_radius', search_radius)

        T.cut(I)
        log('Cut to', len(T), 'CCDs near bricks')

    if opt.forced:
        log('Writing forced-photometry commands to', opt.out)
        f = open(opt.out,'w')
        log('Total of', len(T), 'CCDs')
        #T.cut(allI)
        camexp = set([(c,e) for c,e in zip(T.camera, T.expnum)])
        print(len(camexp), 'unique camera/exposure pairs')
        for cam,exp in camexp:
            expstr = '%08i' % exp
            outfn = os.path.join('forced', cam, expstr[:5], 'forced-%s-%s.fits' % (cam, exp))
            f.write('%s %s all %s\n' % (cam, exp, outfn))
        f.close()
        log('Wrote', opt.out)
        return 0

    # sort by RA increasing
    B.cut(np.argsort(B.ra))

    if opt.save_to_fits:
        assert(opt.touching)
        # Write cut tables to file
        for tab,typ in zip([B,T],['bricks','ccds']):
            fn='%s-%s-cut.fits' % (typ,opt.region)
            if os.path.exists(fn):
                os.remove(fn)
            tab.writeto(fn)
            log('Wrote %s' % fn)
        # Write text files listing ccd and filename names
        # nm1,nm2= 'ccds-%s.txt'% opt.region,'filenames-%s.txt' % opt.region
        # if os.path.exists(nm1):
        #     os.remove(nm1)
        # if os.path.exists(nm2):
        #     os.remove(nm2)
        # f1,f2=open(nm1,'w'),open(nm2,'w')
        # fns= list(set(T.get('image_filename')))
        # for fn in fns:
        #     f2.write('%s\n' % fn.strip())
        # for ti in T:
        #     f1.write('%s\n' % ti.get('image_filename').strip())
        # f1.close()
        # f2.close()
        # log('Wrote *-names.txt')

    if opt.touching:

        if want_bricks:
            # Shortcut the list of bricks that are definitely touching CCDs --
            # a brick-ccd pair within this radius must be touching.
            closest_radius = 0.95 * (bricksize + 0.262 * 2048 / 3600.) / 2.

            J1,nil,nil = match_radec(B.ra, B.dec, T.ra, T.dec, closest_radius, nearest=True)
            log(len(J1), 'of', len(B), 'bricks definitely touch CCDs')
            tocheck = np.ones(len(B), bool)
            tocheck[J1] = False
            J2 = []
            for j in np.flatnonzero(tocheck):
                b = B[j]
                wcs = wcs_for_brick(b)
                I = ccds_touching_wcs(wcs, T)
                log(len(I), 'CCDs for brick', b.brickname)
                if len(I) == 0:
                    continue
                J2.append(j)
            J = np.hstack((J1, J2))
            J = np.sort(J).astype(int)
            B.cut(J)
            log('Cut to', len(B), 'bricks touching CCDs')

        else:
            J = []
            allI = set()
            for j,b in enumerate(B):
                wcs = wcs_for_brick(b)
                I = ccds_touching_wcs(wcs, T)
                log(len(I), 'CCDs for brick', b.brickname)
                if len(I) == 0:
                    continue
                allI.update(I)
                J.append(j)
            allI = list(allI)
            allI.sort()
            B.cut(np.array(J))
            log('Cut to', len(B), 'bricks touching CCDs')

    elif opt.near:
        # Find CCDs near bricks
        allI,nil,nil = match_radec(T.ra, T.dec, B.ra, B.dec, search_radius, nearest=True)
        # Find bricks near CCDs
        J,nil,nil = match_radec(B.ra, B.dec, T.ra, T.dec, search_radius, nearest=True)
        B.cut(J)
        log('Cut to', len(B), 'bricks near CCDs')

    else:
        allI = np.arange(len(T))

    if opt.byexp:
        nil,eI = np.unique(T.expnum[allI], return_index=True)
        allI = allI[eI]
        print('Cut to', len(allI), 'expnums')

    if opt.nccds:
        from queue import Queue
        from threading import Thread

        log('Checking number of CCDs per brick')

        def worker():
            while True:
                i = q.get()
                if i is None:
                    break
                b = B[i]
                wcs = wcs_for_brick(b)
                I = ccds_touching_wcs(wcs, T)
                log(b.brickname, len(I))
                q.task_done()

        q = Queue()
        num_threads = 24
        threads = []

        for i in range(num_threads):
            t = Thread(target=worker)
            t.start()
            threads.append(t)

        for i in range(len(B)):
            q.put(i)

        q.join()
        for i in range(num_threads):
            q.put(None)
        for t in threads:
            t.join()

    if opt.write_ccds:
        T[allI].writeto(opt.write_ccds)
        log('Wrote', opt.write_ccds)

    if want_bricks:
        # Print the list of bricks and exit.
        for b in B:
            print(b.brickname)
        if opt.save_to_fits:
            B.writeto('bricks-%s-touching.fits' % opt.region)
        if not want_ccds:
            return 0

    ## Be careful here -- T has been cut; we want to write out T.index.
    ## 'allI' contains indices into T.

    if opt.stage is not None:
        cmd_pat = 'rsync -LRarv %s %s'
        fns = set()
        for iccd in allI:
            im = survey.get_image_object(T[iccd])
            fns.update([im.imgfn, im.wtfn, im.dqfn, im.psffn, im.merged_psffn,
                   im.merged_splineskyfn, im.splineskyfn])
        for i,fn in enumerate(fns):
            print('File', i+1, 'of', len(fns), ':', fn)
            if not os.path.exists(fn):
                print('No such file:', fn)
                continue
            base = survey.get_survey_dir()
            if base.endswith('/'):
                base = base[:-1]

            rel = os.path.relpath(fn, base)

            dest = os.path.join(opt.stage, rel)
            print('Dest:', dest)
            if os.path.exists(dest):
                print('Exists:', dest)
                continue

            cmd = cmd_pat % ('%s/./%s' % (base, rel), opt.stage)
            print(cmd)
            rtn = os.system(cmd)
            assert(rtn == 0)
        return 0

    if opt.lsb:
        log('Writing LSB commands to', opt.out)
        f = open(opt.out,'w')
        log('Total of', len(allI), 'CCDs')
        for j,i in enumerate(allI):
            exp = T.expnum[i]
            ext = T.ccdname[i].strip()
            outfn = 'lsb/lsb-%s-%s.fits' % (exp, ext)
            f.write('python legacyanalysis/lsb.py --expnum %i --extname %s --out %s -F -n > lsb/lsb-%s-%s.log 2>&1\n' % (exp, ext, outfn, exp, ext))
        f.close()
        log('Wrote', opt.out)
        return 0

    log('Writing calibs to', opt.out)
    f = open(opt.out,'w')
    log('Total of', len(allI), 'CCDs')

    batch = []

    def write_batch(f, batch, cmd):
        if cmd is None:
            cmd = ''
        f.write(cmd + ' '.join(batch) + '\n')

    cmd = None
    if opt.command:
        cmd = 'python legacypipe/run-calib.py '
        if opt.opt is not None:
            cmd += opt.opt + ' '

    for j,i in enumerate(allI):

        if opt.delete_sky:
            log(j+1, 'of', len(allI))
            im = survey.get_image_object(T[i])
            if opt.delete_sky and os.path.exists(im.skyfn):
                log('  deleting:', im.skyfn)
                os.unlink(im.skyfn)

        if opt.command:
            if opt.byexp:
                s = '--expnum %i' % (T.expnum[i])
            else:
                s = '%i-%s' % (T.expnum[i], T.ccdname[i])
            prefix = 'python legacypipe/run-calib.py '
            if opt.opt is not None:
                prefix = prefix + opt.opt
            #('python legacypipe/run-calib.py --expnum %i --ccdname %s' %
            #     (T.expnum[i], T.ccdname[i]))
        else:
            s = '%i' % T.index[i]
            prefix = ''

        if j < 10:
            print('Index', T.index[i], 'expnum', T.expnum[i], 'ccdname', T.ccdname[i],
                  'filename', T.image_filename[i])

        if not opt.nper:
            f.write(prefix + s + '\n')
        else:
            batch.append(s)
            if len(batch) >= opt.nper:
                write_batch(f, batch, cmd)
                batch = []

    if len(batch):
        write_batch(f, batch, cmd)

    f.close()
    log('Wrote', opt.out)
    return 0
#
#
#
if __name__ == '__main__':
    import sys
    main(sys.argv[1:])

ext: py

max_stars_repo_name: michaelJwilson/legacypipe

max_issues_repo_name: michaelJwilson/legacypipe

max_forks_repo_name: michaelJwilson/legacypipe

reason: import numpy


text:
import numpy as np
import os
from sklearn.preprocessing import MinMaxScaler
import pickle
import tensorflow as tf
import vae_
from sklearn.utils import shuffle
import collections

STEM_SAMPLE_POINT_NUM = 20

data_size = 40*3
batch_size = 50

# settings for VAE.
n_hidden = 100  # num of hidden units in mlp.
dim_z = 20  # dim of latent.

br_data_size = 4 * dim_z + 3

# data_path = "C:/Users/Administrator/Desktop/neuron_factor_tree_test/test/data/L4_BC/data"
#
# sess_save_dir = "C:/Users/Administrator/Desktop/neuron_factor_tree_test/test/data/L4_BC/data/save"
# save_dir = "C:/Users/Administrator/Desktop/neuron_factor_tree_test/test/data/L4_BC/data/save"

'''
Load both branch codes and branch id sequences.
'''

class SeqLoader():
    def __init__(self, test_path, seq_size, branch_size, batch_size, seq_type):
        self.data_path = os.path.join(test_path, "seq_data")
        self.save_path = os.path.join(test_path, "vae_save")
        self.seq_size = seq_size
        self.branch_size = branch_size
        self.batch_size = batch_size
        self.type = seq_type

        # load rnn data.
        training_data = os.path.join(self.data_path, "factorDB"+"_"+str(self.type)+".txt")
        if not (os.path.exists(training_data)):
            self.build_seq_data()
        else:
            self.load_processed(self.data_path)
        # self.build_seq_data()

        self.create_batches()
        self.reset_batch()

    # Collect all the branch codes(normalized).
    def branches_loader(self, branch_codes_file):
        '''
        A Branch : [angle, stm_scale_1, stm_scale_2, br_codes[sample_point_num, 3]].
        '''
        branches = []
        branch = []
        branch_codes = []
        br_angle = 0.0
        stm_scale = 1.0
        with open(branch_codes_file) as f:
            for line in f.readlines():
                if line.startswith("#BRANCH:"):
                    branch = []
                    branch_codes = []
                elif line.startswith("#ANGLE:"):
                    br_angle = float(line[:-1].split(" ")[1])
                    branch.append(br_angle)
                elif line.startswith("#STEM:"):
                    stm_scale = float(line[:-1].split(" ")[1])
                    branch.append(stm_scale)
                elif line.startswith("#DATA:"):
                    code_str = np.array(line[:-1].split(" ")[1:4])
                    code = list(map(float, code_str))
                    branch_codes.append(code)
                elif line.startswith("#BRANCHEND"):
                    branch.append(branch_codes)
                    # print(branch)
                    branches.append(branch)

        # print(len(branches))
        return branches

    # Load sequence as branch ids, only load a type each time.
    def sequences_loader(self, branch_database_path):
        seqs = []
        seq = []
        current_type = 0
        with open(branch_database_path) as f:
            for line in f.readlines():
                if line.startswith("#SEQ"):
                    if current_type==self.type and len(seq) != 0:
                        seqs.append(seq)
                    seq = []
                    current_type = int(np.array(line[:-1].split(" ")[1])) # update current seq type

                elif line.startswith("#BID"):
                    b_id_str = np.array(line[:-1].split(" ")[1])
                    b_id = float(b_id_str)
                    seq.append(b_id)

        if current_type==self.type and len(seq)!=0:
            seqs.append(seq)

        # print(len(seqs))
        return seqs

    def process_branches(self, branches):

        print("Processing branches.")

        # load scaler_file created in branch_generator traning step.
        scaler_file = os.path.join(self.data_path, "scaler.txt")
        with open(scaler_file, "rb") as f:
            min_max_scaler = pickle.load(f)
        print("Scaler file loaded from [%s]." % scaler_file)
        self.min_max_scaler = min_max_scaler  # record the scaler.

        br_count = np.array([[br[0], br[1], br[2]] for br in branches])
        # print(br_count.shape)

        angle_stm_scaler = MinMaxScaler()
        angle_stm_scaler.fit(br_count)

        angle_stm_scaler_file = os.path.join(self.data_path, "angle_stm_scaler.txt")
        with open(angle_stm_scaler_file, "wb") as f:
            pickle.dump(angle_stm_scaler, f)
        print("Angle stm scaler file saved at [%s]." % angle_stm_scaler_file)

        # Normalixation
        data = [[angle_stm_scaler.transform(np.array([br[0], br[1], br[2]]).reshape(-1, 3)),
                 min_max_scaler.transform(br[3]).reshape(STEM_SAMPLE_POINT_NUM*2, 3,  1)]
                for br in branches]
        # print(len(data))
        # print(data[0])

        return data

    def build_seq_data(self):

        print("Building sequence datas.")

        branch_seq_path = os.path.join(self.data_path, "branch_dataset.brd")
        branch_code_path = os.path.join(self.data_path, "branch_code.bcd")
        # Load branches.
        branches = self.branches_loader(branch_code_path)
        print("{} branches loaded.".format(len(branches)))
        self.branches = self.process_branches(branches)  # [None, [angle, scale_1, scale_2, codes]].

        # Load branch id sequences.
        self.sequences = self.sequences_loader(branch_seq_path)
        print("{} sequences loaded.".format(len(self.sequences)))
        if int(len(self.sequences) / self.batch_size) == 0:
            self.tensor = []
            return
        # print(self.sequences)

        # Sequence normalization.
        seq_ids = []
        for seq in self.sequences:
            seq.reverse()
            seq_n = []
            for id in [int(f) for f in seq]:
                if id!=-1: # empty branches.
                    seq_n.append(id)
            if(len(seq_n)==0):
                continue
            elif len(seq_n) < self.seq_size:
                for i in range(self.seq_size - len(seq_n)):
                    seq_n.append(-1)
            else:  # equals or shorter than seq_size.
                seq_n = seq_n[(len(seq_n) - self.seq_size) + 1:]
                seq_n.append(-1) # make sure last one is -1.

            seq_ids.append(seq_n)

        # print("shape")
        # print(np.array(seq_ids).shape)

        # x = tf.placeholder(tf.float32, shape=[None, data_size], name='input_img')
        x = tf.placeholder(tf.float32, shape=[None, 20, 3, 1], name='input')
        latent = vae_.gaussian_MLP_encoder(x, n_hidden, dim_z, 1.0)

        with tf.Session() as sess:

            tf.global_variables_initializer().run()
            saver = tf.train.Saver(tf.global_variables())
            ckpt = tf.train.get_checkpoint_state(self.save_path)
            if ckpt:
                print("VAE model restored for branch latent transform.")
                saver.restore(sess, ckpt.model_checkpoint_path)
            else:
                return

            latent_data = []
            for seq in seq_ids:
                seq_latents = []
                # seq_latents.append(np.array([1.0 for _ in range(br_data_size)])) # all 1.0 for soma.
                for br in seq:
                    if br != -1:
                        _st_mean, _st_stdv = \
                            sess.run(latent, feed_dict={x: [self.branches[br][1][:STEM_SAMPLE_POINT_NUM]]})# stem 1
                        st_latent = np.append(_st_mean[0], _st_stdv[0])
                        _st_mean, _st_stdv = \
                            sess.run(latent, feed_dict={x: [self.branches[br][1][STEM_SAMPLE_POINT_NUM:]]})
                        st_latent_ = np.append(_st_mean[0], _st_stdv[0])
                        br_latent = np.append(np.array(self.branches[br][0]), st_latent)
                        br_latent = np.append(br_latent, st_latent_)
                        # print(len(br_latent))
                        # print(br_latent)
                        seq_latents.append(br_latent)
                for i in range(self.seq_size - len(seq_latents)):
                    seq_latents.append(np.array([0.0 for _ in range(self.branch_size)])) # all 0 for leaf.

                latent_data.append(seq_latents)
            # print(len(latent_data))
            # print(latent_data[0][0])
            x_data_ = np.array(latent_data, dtype=np.float32)
            print(x_data_.shape)  # (data_size, seq_size, 4*dim_z+3)

        print("Training data normalized.")

        x_data = shuffle(x_data_, random_state=42)

        # make up a neuronFactor-database. Each nf in DB have a id.
        factors_ = x_data.reshape(-1, self.branch_size)  # reshape to count factors.
        factors_s = [str(f) for f in factors_.tolist()]
        counter = collections.Counter(factors_s)
        count_pairs = sorted(counter.items(), key=lambda x: -x[1])
        self.factors, _ = zip(*count_pairs)
        print(len(self.factors))
        # for f in self.factors:
        #     print(f)
        self.factorDB_size = len(self.factors)
        self.factorDB = dict(zip(self.factors, range(len(self.factors))))
        db_file = os.path.join(self.data_path, "factorDB" + "_" + str(self.type) + ".txt")
        print("Got a NeuronFactor Database of size %d." % self.factorDB_size)
        tensor = np.array(list(map(self.factorDB.get, factors_s)))
        tensor_file = os.path.join(self.data_path, "tensor" + "_" + str(self.type) + "")
        self.tensor = tensor.reshape(-1, self.seq_size)  # seq_num * seq_size

        # save data.
        with open(db_file, "wb") as f:
            pickle.dump(self.factors, f)
            print("NeuronFactorDB saved at [%s]." % db_file)
        np.save(tensor_file, self.tensor)
        print("Factor tensors saved at [%s]." % tensor_file)

        print("Data preprocessed.")

    def load_processed(self, data_base_dir):
        db_file = os.path.join(data_base_dir, "factorDB_"+str(self.type)+".txt")
        tensor_file = os.path.join(data_base_dir, "tensor_"+str(self.type)+".npy")

        # Loading data.
        with open(db_file, "rb") as f:
            self.factors = pickle.load(f)
        self.factorDB_size = len(self.factors)
        self.factorDB = dict(zip(self.factors, range(len(self.factors))))
        self.tensor = np.load(tensor_file)

        print("[%d] training data loaded from [%s]." % (len(self.tensor), tensor_file))

    def create_batches(self):
        self.num_batches = int(len(self.tensor) / self.batch_size)
        if self.num_batches== 0:
            print("Don't have enough data, finish.")
            return
        # use the first self.batch_size*self.num_batches data.
        tensor = self.tensor[:self.batch_size*self.num_batches]
        x_data = tensor
        y_data = np.copy(tensor)
        for y_ in y_data:
            y_[:-1] = y_[1:]
            y_[-1] = 0

        self.x_batches = np.split(x_data, self.num_batches)
        self.y_batches = np.split(y_data, self.num_batches)

    def reset_batch(self):
        self.index_id = 0

    def next_batch(self):  # in test use for range.
        x, y = self.x_batches[self.index_id], \
               self.y_batches[self.index_id]
        self.index_id += 1
        return x, y

    def get_seq_ter_prob(self):
        branch_seq_path = os.path.join(self.data_path, "branch_dataset.brd")
        seqs = self.sequences_loader(branch_seq_path)
        seq_len_list = [len(seq) for seq in seqs]
        seq_len_counter = collections.Counter(seq_len_list)

        t_prob = {}  # A dict for terminate prob.
        for l, _ in seq_len_counter.items():
            c = sum([c_ for l_, c_ in seq_len_counter.items() if l_ == l])
            # rest longer one
            r = 0  # rest longer one
            for l_, c_ in seq_len_counter.items():
                r_ = 0
                if l_ > l:  # for longer one
                    for i in range(l_ - l):
                        c_ = c_ / 2  # get rid of repetitive sub tree
                    r_ += c_
                r += r_

            t_prob[l] = float(c) / float(c + r)

        print("Got termination prob.")

        return t_prob

# loader = SeqLoader("C:/Users/Administrator/Desktop/TestFullNeuron", 20, 2*STEM_SAMPLE_POINT_NUM, 50, 3)
# x, y = loader.next_batch()
# print(x)
# print(y)


text:
import copy
import logging
import numpy as np
import os
import SimpleITK as sitk
import tensorflow as tf
import time
import functions.image.image_processing as ip
import functions.network as network
import functions.setting.setting_utils as su
from functions.setting.experiment_setting import load_network_multi_stage_from_predefined

def multi_stage_error(setting, landmark_dict, pair_info, overwrite=False):
    """
    :param setting:
    :param pair_info: information of the pair to be registered.
    :param overwrite:
    :return: The output moved images and dvf will be written to the disk.
             1: registration is performed correctly
             2: skip overwriting
             3: the dvf is available from the previous experiment [4, 2, 1]. Then just upsample it.
    """
    stage_list = setting['ImagePyramidSchedule']
    if setting['CNN_Mode']:
        lstm_mode = False
    else:
        lstm_mode = True
    dvf_error_address = su.address_generator(setting, 'dvf_error', current_experiment=setting['lstm_exp'], pair_info=pair_info, stage_list=stage_list, base_reg=setting['BaseReg'])
    if os.path.isfile(dvf_error_address):
        if not overwrite:
            logging.debug('overwrite=False, file ' + dvf_error_address + ' already exists, skipping .....')
            return 2
        else:
            logging.debug('overwrite=True, file ' + dvf_error_address + ' already exists, but overwriting .....')

    im_dict = dict()
    im_info_fixed = pair_info[0]
    im_info_moving = pair_info[1]
    base_reg = copy.copy(setting['BaseReg'])
    # im_info_moving['base_reg'] = copy.copy(base_reg)

    network_lstm_design = setting['network_lstm_dict']['NetworkDesign']
    setting['stage'] = 1
    setting = su.load_network_setting(setting, network_name=network_lstm_design)
    # setting['ImPad_S1'] = 50  # more safe selection
    if im_info_fixed['data'] == 'DIR-Lab_COPD':
        if setting['ImPad_S1'] < 42:
            setting['ImPad_S1'] = 42

    if im_info_fixed['data'] == 'DIR-Lab_4D':
        if setting['ImPad_S1'] < 40:
            setting['ImPad_S1'] = 40

    setting['ImPad_S2'] = setting['ImPad_S1'] + 7  # 5
    setting['ImPad_S4'] = setting['ImPad_S1'] + 10  # 8

    base_key_list = ['fixed_im_s', 'fixed_mask_s', 'moved_im_s', 'moved_mask_s']
    for i_stage, stage in enumerate(stage_list):
        mask_to_zero_stage = setting['network_dict']['stage' + str(stage)]['MaskToZero']
        fixed_im_stage_address = su.address_generator(setting, 'Im', stage=stage, **im_info_fixed)
        fixed_mask_stage_address = su.address_generator(setting, mask_to_zero_stage, stage=stage, **im_info_fixed)
        im_dict['fixed_im_s'+str(stage)] = sitk.GetArrayFromImage(sitk.ReadImage(fixed_im_stage_address))
        im_dict['fixed_mask_s' + str(stage)] = sitk.GetArrayFromImage(sitk.ReadImage(fixed_mask_stage_address))

        if i_stage == 0:
            check_downsampled_base_reg(setting, stage, base_reg=base_reg, pair_info=pair_info, mask_to_zero_stage=mask_to_zero_stage)
            moved_im_stage_address = su.address_generator(setting, 'MovedImBaseReg', pair_info=pair_info, stage=stage, base_reg=base_reg, **im_info_moving)
            moved_mask_stage_address = su.address_generator(setting, 'Moved'+mask_to_zero_stage+'BaseReg', pair_info=pair_info, stage=stage, base_reg=base_reg, **im_info_moving)
        else:
            if setting['UseRegisteredImages']:
                moved_im_stage_address = su.address_generator(setting, 'MovedIm', current_experiment=setting['exp_multi_reg'], pair_info=pair_info,
                                         stage=stage, stage_list=stage_list, base_reg=base_reg)
                moved_mask_stage_address = su.address_generator(setting, 'Moved'+mask_to_zero_stage, current_experiment=setting['exp_multi_reg'],
                                         pair_info=pair_info, stage=stage, stage_list=stage_list, base_reg=base_reg)

            else:
                check_downsampled_base_reg(setting, stage, base_reg=base_reg, pair_info=pair_info, mask_to_zero_stage=mask_to_zero_stage)
                moved_im_stage_address = su.address_generator(setting, 'MovedImBaseReg', pair_info=pair_info,
                                         stage=stage, base_reg=base_reg, **im_info_moving)
                moved_mask_stage_address = su.address_generator(setting, 'Moved'+mask_to_zero_stage+'BaseReg', pair_info=pair_info,
                                         stage=stage, base_reg=base_reg, **im_info_moving)

        im_dict['moved_im_s' + str(stage)] = sitk.GetArrayFromImage(sitk.ReadImage(moved_im_stage_address))
        im_dict['moved_mask_s' + str(stage)] = sitk.GetArrayFromImage(sitk.ReadImage(moved_mask_stage_address))
        logging.info(fixed_im_stage_address + ' loaded')
        logging.info(fixed_mask_stage_address + ' loaded')
        logging.info(moved_im_stage_address + ' loaded')
        logging.info(moved_mask_stage_address + ' loaded')

        default_pixel = setting['data'][pair_info[0]['data']]['DefaultPixelValue']
        im_dict['fixed_im_s'+str(stage)][im_dict['fixed_mask_s'+str(stage)] == 0] = default_pixel
        im_dict['moved_im_s'+str(stage)][im_dict['moved_mask_s'+str(stage)] == 0] = default_pixel

        if setting['ImPad_S'+str(stage)] > 0:
            for base_key in base_key_list:
                im_dict[base_key + str(stage)] = np.pad(im_dict[base_key + str(stage)],
                                                        setting['ImPad_S' + str(stage)],
                                                        'constant',
                                                        constant_values=(default_pixel,))

    indices_padded = copy.deepcopy(landmark_dict['FixedLandmarksIndex'])
    indices_padded = indices_padded + setting['ImPad_S1']

    batch_size = 15
    for i in range(3):
        landmark_dict['DVF_error_times'+str(i)+'_logits'] = np.zeros((np.shape(landmark_dict['FixedLandmarksIndex'])[0],) + (setting['NumberOfLabels'],))
        landmark_dict['DVF_error_times'+str(i)+'_label'] = np.zeros(np.shape(landmark_dict['FixedLandmarksIndex'])[0])
    multi_stage_network_design = load_network_multi_stage_from_predefined(setting['exp_multi_reg'])
    multi_stage_network_load = get_parameter_multi_stage_network(setting, multi_stage_network_design)
    multi_stage_network_address = dict()
    for stage in setting['ImagePyramidSchedule']:
        multi_stage_network_address['stage'+str(stage)] = su.address_generator(setting, 'saved_model_with_step',
                                                                               current_experiment=multi_stage_network_load['stage'+str(stage)]['NetworkLoad'],
                                                                               step=multi_stage_network_load['stage'+str(stage)]['GlobalStepLoad'])
    # Network
    tf.reset_default_graph()
    tf.set_random_seed(0)
    with tf.variable_scope('InputImages'):
        images_s1_tf = tf.placeholder(tf.float32, shape=[batch_size, 2 * setting['R'] + 1, 2 * setting['R'] + 1, 2 * setting['R'] + 1, 2], name="Images_S1")
        images_s2_tf = tf.placeholder(tf.float32, shape=[batch_size, 2 * setting['R'] + 1, 2 * setting['R'] + 1, 2 * setting['R'] + 1, 2], name="Images_S2")
        images_s4_tf = tf.placeholder(tf.float32, shape=[batch_size, 2 * setting['R'] + 1, 2 * setting['R'] + 1, 2 * setting['R'] + 1, 2], name="Images_S4")
    bn_training = tf.placeholder(tf.bool, name='bn_training')

    if lstm_mode:
        out0_tf, out1_tf, out2_tf, state0_tf, state1_tf = getattr(getattr(network, network_lstm_design), 'network')(images_s1_tf, images_s2_tf, images_s4_tf,
                                                                                bn_training,
                                                                                detailed_summary=False,
                                                                                use_keras=setting['use_keras'],
                                                                                num_of_classes=setting['NumberOfLabels'],
                                                                                multi_stage_network_address=multi_stage_network_address
                                                                                )
    else:
        out2_tf = getattr(getattr(network, network_lstm_design), 'network')(images_s1_tf, images_s2_tf, images_s4_tf,
                                                                            bn_training,
                                                                            detailed_summary=setting['DetailedNetworkSummary'],
                                                                            use_keras=setting['use_keras'],
                                                                            num_of_classes=setting['NumberOfLabels'],
                                                                            multi_stage_network_address=multi_stage_network_address
                                                                            )
        out0_tf, out1_tf, state0_tf, state1_tf = None, None, None, None

    sess = tf.Session()

    saver_loading = tf.train.Saver(tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES))
    saver_loading.restore(sess, su.address_generator(setting, 'saved_model_with_step',
                                                     current_experiment=setting['network_lstm_dict']['NetworkLoad'],
                                                     step=setting['network_lstm_dict']['GlobalStepLoad']))

    time_before_dvf = time.time()
    begin_landmark_i = 0
    end_landmark_i = 0
    fill_last_batch = 0
    while end_landmark_i < np.shape(landmark_dict['DVF_error_times2_logits'])[0]:
        end_landmark_i = begin_landmark_i + batch_size
        if end_landmark_i > np.shape(landmark_dict['DVF_error_times2_logits'])[0]:
            fill_last_batch = end_landmark_i - np.shape(landmark_dict['DVF_error_times2_logits'])[0]
            end_landmark_i = np.shape(landmark_dict['DVF_error_times2_logits'])[0]
        logging.info('begin_landmark_i={}, end_landmark_i={}'.format(begin_landmark_i, end_landmark_i))
        batch_im = extract_batch_from_patch_seq(setting, im_dict, indices_padded[begin_landmark_i: end_landmark_i, :])
        if fill_last_batch:
            batch_im = fill_junk_to_batch_size(setting, batch_size, batch_im)
        out_np = [None for _ in range(3)]
        if lstm_mode:
            [out_np[0], out_np[1], out_np[2]] = sess.run([out0_tf, out1_tf, out2_tf],
                                                         feed_dict={images_s1_tf: batch_im['stage1'], images_s2_tf: batch_im['stage2'],
                                                                    images_s4_tf: batch_im['stage4'], bn_training: 0})
        else:
            [out_np[2]] = sess.run([out2_tf],
                                                         feed_dict={images_s1_tf: batch_im['stage1'], images_s2_tf: batch_im['stage2'],
                                                                    images_s4_tf: batch_im['stage4'], bn_training: 0})

        for i in range(3):
            if lstm_mode or i == 2:
                out_np_center = out_np[i][:, setting['Ry']+1, setting['Ry']+1, setting['Ry']+1, :]
                out_np_center_label = np.argmax(out_np_center[:, setting['Labels_time'+str(i)]], axis=1)
                landmark_dict['DVF_error_times'+str(i)+'_label'][begin_landmark_i:end_landmark_i] = out_np_center_label[0: (batch_size - fill_last_batch)]
                landmark_dict['DVF_error_times'+str(i)+'_logits'][begin_landmark_i:end_landmark_i, :] = out_np_center[0: (batch_size - fill_last_batch), :]
        begin_landmark_i += 15

    time_after_dvf = time.time()
    logging.debug(pair_info[0]['data']+', CN{}, ImType{} is done in {:.2f}s '.format(
        pair_info[0]['cn'], pair_info[0]['type_im'], time_after_dvf - time_before_dvf))

    landmark_dict['DVF_nonrigidGroundTruth_magnitude'] = np.sqrt(landmark_dict['DVF_nonrigidGroundTruth'][:, 0]**2 +
                                                                 landmark_dict['DVF_nonrigidGroundTruth'][:, 1]**2 +
                                                                 landmark_dict['DVF_nonrigidGroundTruth'][:, 2]**2)

    result_landmarks_folder = su.address_generator(setting, 'result_landmarks_folder', stage_list=stage_list, current_experiment=setting['lstm_exp'])
    if not os.path.isdir(result_landmarks_folder):
        os.makedirs(result_landmarks_folder)

    pair_info_text = setting['BaseReg'] + '_' +  pair_info[0]['data'] + \
                     '_CN{}_TypeIm{},'.format(pair_info[0]['cn'], pair_info[0]['type_im']) + '  Moving:' + \
                     pair_info[1]['data'] + '_CN{}_TypeIm{}'.format(pair_info[1]['cn'], pair_info[1]['type_im'])

    # import matplotlib.pyplot as plt
    # plt.figure()
    # plt.plot(landmark_dict['DVF_error_times2_label'], landmark_dict['DVF_nonrigidGroundTruth_magnitude'], 'o')
    # plt.draw()
    # plt.savefig(result_landmarks_folder+pair_info_text+'.png')
    # plt.close()

    return landmark_dict

def get_parameter_multi_stage_network(setting, network_dict):
    for key in network_dict.keys():
        net = network_dict[key]
        # if net['R'] == 'Auto' and net['Ry'] == 'Auto':
        #     net['R'], net['Ry'] = getattr(getattr(network, net['NetworkDesign']), 'radius_test')()

        if net['Ry_erode'] == 'Auto':
            if net['NetworkDesign'] == 'crop5_connection':
                net['Ry_erode'] = 0
            else:
                net['Ry_erode'] = 2
        net['GlobalStepLoad'] = su.get_global_step(setting, net['GlobalStepLoad'], net['NetworkLoad'])

    return network_dict

def extract_batch_from_patch_seq(setting, im_dict, indices_padded):
    r = setting['R']
    stage_list = setting['ImagePyramidSchedule']

    ind = dict()
    ind['S1_no_padding'] = copy.deepcopy(indices_padded - setting['ImPad_S1'])

    for stage in stage_list:
        if stage != 1:
            ind['S'+str(stage)+'_no_padding'] = np.zeros(ind['S1_no_padding'].shape, dtype=np.int16)
            ind['S'+str(stage)+'_no_padding'] = np.zeros(ind['S1_no_padding'].shape, dtype=np.int16)
            for patch_i in range(ind['S1_no_padding'].shape[0]):
                ind['S'+str(stage)+'_no_padding'][patch_i, :] = np.round(ind['S1_no_padding'][patch_i, :]/stage)

        ind['S'+str(stage)] = ind['S'+str(stage)+'_no_padding'] + setting['ImPad_S'+str(stage)]

    batch_both = dict()
    batch_first_s1 = np.stack([im_dict['moved_im_s1'][
                            ind['S1'][i_batch][2] - r: ind['S1'][i_batch][2] + r + 1,
                            ind['S1'][i_batch][1] - r: ind['S1'][i_batch][1] + r + 1,
                            ind['S1'][i_batch][0] - r: ind['S1'][i_batch][0] + r + 1,
                            np.newaxis] for i_batch in range(np.shape(indices_padded)[0])])

    batch_second_s1 = np.stack([im_dict['fixed_im_s1'][
                            ind['S1'][i_batch][2] - r: ind['S1'][i_batch][2] + r + 1,
                            ind['S1'][i_batch][1] - r: ind['S1'][i_batch][1] + r + 1,
                            ind['S1'][i_batch][0] - r: ind['S1'][i_batch][0] + r + 1,
                            np.newaxis] for i_batch in range(np.shape(indices_padded)[0])])

    batch_both_s1 = np.concatenate((batch_first_s1, batch_second_s1), axis=4)
    batch_both['stage1'] = batch_both_s1

    batch_first_s2 = np.stack([im_dict['moved_im_s2'][
                            ind['S2'][i_batch][2] - r: ind['S2'][i_batch][2] + r + 1,
                            ind['S2'][i_batch][1] - r: ind['S2'][i_batch][1] + r + 1,
                            ind['S2'][i_batch][0] - r: ind['S2'][i_batch][0] + r + 1,
                            np.newaxis] for i_batch in range(np.shape(indices_padded)[0])])

    batch_second_s2 = np.stack([im_dict['fixed_im_s2'][
                            ind['S2'][i_batch][2] - r: ind['S2'][i_batch][2] + r + 1,
                            ind['S2'][i_batch][1] - r: ind['S2'][i_batch][1] + r + 1,
                            ind['S2'][i_batch][0] - r: ind['S2'][i_batch][0] + r + 1,
                            np.newaxis] for i_batch in range(np.shape(indices_padded)[0])])

    batch_both_s2 = np.concatenate((batch_first_s2, batch_second_s2), axis=4)
    batch_both['stage2'] = batch_both_s2

    batch_first_s4 = np.stack([im_dict['moved_im_s4'][
                            ind['S4'][i_batch][2] - r: ind['S4'][i_batch][2] + r + 1,
                            ind['S4'][i_batch][1] - r: ind['S4'][i_batch][1] + r + 1,
                            ind['S4'][i_batch][0] - r: ind['S4'][i_batch][0] + r + 1,
                            np.newaxis] for i_batch in range(np.shape(indices_padded)[0])])

    batch_second_s4 = np.stack([im_dict['fixed_im_s4'][
                            ind['S4'][i_batch][2] - r: ind['S4'][i_batch][2] + r + 1,
                            ind['S4'][i_batch][1] - r: ind['S4'][i_batch][1] + r + 1,
                            ind['S4'][i_batch][0] - r: ind['S4'][i_batch][0] + r + 1,
                            np.newaxis] for i_batch in range(np.shape(indices_padded)[0])])
    batch_both_s4 = np.concatenate((batch_first_s4, batch_second_s4), axis=4)
    batch_both['stage4'] = batch_both_s4

    return batch_both

def fill_junk_to_batch_size(setting, batch_size, batch_im):
    for stage in setting['ImagePyramidSchedule']:
        if np.shape(batch_im['stage'+str(stage)])[0] < batch_size:
            diff_size = batch_size - np.shape(batch_im['stage'+str(stage)])[0]
            diff_np = np.zeros((diff_size,) + np.shape(batch_im['stage'+str(stage)])[1:4] + (2,))
            batch_im['stage'+str(stage)] = np.concatenate((batch_im['stage'+str(stage)], diff_np), axis=0)
    return batch_im

def check_downsampled_base_reg(setting, stage, base_reg=None, pair_info=None, mask_to_zero_stage=None):
    if 'DownSamplingByGPU' not in setting.keys():
        setting['DownSamplingByGPU'] = False
    im_info_moving = pair_info[1]
    im_list_downsample = [{'Image': 'MovedImBaseReg',
                           'interpolator': 'BSpline',
                           'DefaultPixelValue': setting['data'][im_info_moving['data']]['DefaultPixelValue'],
                           'ImageByte': setting['data'][im_info_moving['data']]['ImageByte']},

                          {'Image': 'Moved' + mask_to_zero_stage + 'BaseReg',
                           'interpolator': 'NearestNeighbor',
                           'DefaultPixelValue': 0,
                           'ImageByte': sitk.sitkInt8},
                          ]
    #
    # im_stage_address = su.address_generator(setting, 'MovedImBaseReg', pair_info=pair_info,
    #                          stage=stage, base_reg=base_reg, **im_info_moving)
    # mask_stage_address = su.address_generator(setting, 'Moved' + mask_to_zero_stage + 'BaseReg', pair_info=pair_info,
    #                          stage=stage, base_reg=base_reg, **im_info_moving)

    for im_dict in im_list_downsample:
        im_stage_address = su.address_generator(setting, im_dict['Image'], stage=stage, pair_info=pair_info,
                              base_reg=base_reg, **im_info_moving)

        if not os.path.isfile(im_stage_address):
            im_s1_sitk = sitk.ReadImage(su.address_generator(setting, im_dict['Image'], stage=1, pair_info=pair_info,
                              base_reg=base_reg, **im_info_moving))
            if setting['DownSamplingByGPU'] and im_dict['Image'] == 'MovedImBaseReg':
                im_s1 = sitk.GetArrayFromImage(im_s1_sitk)
                im_stage = ip.downsampler_gpu(im_s1, stage, normalize_kernel=True,
                                              default_pixel_value=im_dict['DefaultPixelValue'])

                im_stage_sitk = ip.array_to_sitk(im_stage,
                                                 origin=im_s1_sitk.GetOrigin(),
                                                 spacing=tuple(i * stage for i in im_s1_sitk.GetSpacing()),
                                                 direction=im_s1_sitk.GetDirection())
            else:
                if im_dict['interpolator'] == 'NearestNeighbor':
                    interpolator = sitk.sitkNearestNeighbor
                elif im_dict['interpolator'] == 'BSpline':
                    interpolator = sitk.sitkBSpline
                else:
                    raise ValueError("interpolator should be in ['NearestNeighbor', 'BSpline']")

                if im_dict['Image'] in ['MovedTorsoBaseReg', 'MovedLungBaseReg']:
                    im_ref_sitk = sitk.ReadImage(su.address_generator(setting, 'MovedImBaseReg', pair_info=pair_info,
                             stage=stage, base_reg=base_reg, **im_info_moving))
                else:
                    im_ref_sitk = None
                im_stage_sitk = ip.resampler_sitk(im_s1_sitk,
                                                  scale=stage,
                                                  im_ref=im_ref_sitk,
                                                  default_pixel_value=im_dict['DefaultPixelValue'],
                                                  interpolator=interpolator)

            sitk.WriteImage(sitk.Cast(im_stage_sitk, im_dict['ImageByte']), im_stage_address)
    return 0

max_stars_repo_name: hsokooti/RegNet

max_issues_repo_name: hsokooti/RegNet

max_forks_repo_name: hsokooti/RegNet


text:
import importlib

import numpy as np
import torch
from skimage import measure
from skimage.metrics import adapted_rand_error, peak_signal_noise_ratio

from pytorch3dunet.unet3d.losses import compute_per_channel_dice
from pytorch3dunet.unet3d.seg_metrics import AveragePrecision, Accuracy
from pytorch3dunet.unet3d.utils import get_logger, expand_as_one_hot, convert_to_numpy

logger = get_logger('EvalMetric')

class DiceCoefficient:
    """Computes Dice Coefficient.
    Generalized to multiple channels by computing per-channel Dice Score
    (as described in https://arxiv.org/pdf/1707.03237.pdf) and theTn simply taking the average.
    Input is expected to be probabilities instead of logits.
    This metric is mostly useful when channels contain the same semantic class (e.g. affinities computed with different offsets).
    DO NOT USE this metric when training with DiceLoss, otherwise the results will be biased towards the loss.
    """

    def __init__(self, epsilon=1e-6, **kwargs):
        self.epsilon = epsilon

    def __call__(self, input, target):
        # Average across channels in order to get the final score
        return torch.mean(compute_per_channel_dice(input, target, epsilon=self.epsilon))

class MeanIoU:
    """
    Computes IoU for each class separately and then averages over all classes.
    """

    def __init__(self, skip_channels=(), ignore_index=None, **kwargs):
        """
        :param skip_channels: list/tuple of channels to be ignored from the IoU computation
        :param ignore_index: id of the label to be ignored from IoU computation
        """
        self.ignore_index = ignore_index
        self.skip_channels = skip_channels

    def __call__(self, input, target):
        """
        :param input: 5D probability maps torch float tensor (NxCxDxHxW)
        :param target: 4D or 5D ground truth torch tensor. 4D (NxDxHxW) tensor will be expanded to 5D as one-hot
        :return: intersection over union averaged over all channels
        """
        assert input.dim() == 5

        n_classes = input.size()[1]

        if target.dim() == 4:
            target = expand_as_one_hot(target, C=n_classes, ignore_index=self.ignore_index)

        assert input.size() == target.size()

        per_batch_iou = []
        for _input, _target in zip(input, target):
            binary_prediction = self._binarize_predictions(_input, n_classes)

            if self.ignore_index is not None:
                # zero out ignore_index
                mask = _target == self.ignore_index
                binary_prediction[mask] = 0
                _target[mask] = 0

            # convert to uint8 just in case
            binary_prediction = binary_prediction.byte()
            _target = _target.byte()

            per_channel_iou = []
            for c in range(n_classes):
                if c in self.skip_channels:
                    continue

                per_channel_iou.append(self._jaccard_index(binary_prediction[c], _target[c]))

            assert per_channel_iou, "All channels were ignored from the computation"
            mean_iou = torch.mean(torch.tensor(per_channel_iou))
            per_batch_iou.append(mean_iou)

        return torch.mean(torch.tensor(per_batch_iou))

    def _binarize_predictions(self, input, n_classes):
        """
        Puts 1 for the class/channel with the highest probability and 0 in other channels. Returns byte tensor of the
        same size as the input tensor.
        """
        if n_classes == 1:
            # for single channel input just threshold the probability map
            result = input > 0.5
            return result.long()

        _, max_index = torch.max(input, dim=0, keepdim=True)
        return torch.zeros_like(input, dtype=torch.uint8).scatter_(0, max_index, 1)

    def _jaccard_index(self, prediction, target):
        """
        Computes IoU for a given target and prediction tensors
        """
        return torch.sum(prediction & target).float() / torch.clamp(torch.sum(prediction | target).float(), min=1e-8)

class AdaptedRandError:
    """
    A functor which computes an Adapted Rand error as defined by the SNEMI3D contest
    (http://brainiac2.mit.edu/SNEMI3D/evaluation).

    This is a generic implementation which takes the input, converts it to the segmentation image (see `input_to_segm()`)
    and then computes the ARand between the segmentation and the ground truth target. Depending on one's use case
    it's enough to extend this class and implement the `input_to_segm` method.

    Args:
        use_last_target (bool): use only the last channel from the target to compute the ARand
    """

    def __init__(self, use_last_target=False, **kwargs):
        self.use_last_target = use_last_target

    def __call__(self, input, target):
        """
        Compute ARand Error for each input, target pair in the batch and return the mean value.

        Args:
            input (torch.tensor): 5D (NCDHW) output from the network
            target (torch.tensor): 4D (NDHW) ground truth segmentation

        Returns:
            average ARand Error across the batch
        """

        def _arand_err(gt, seg):
            n_seg = len(np.unique(seg))
            if n_seg == 1:
                return 0.
            return adapted_rand_error(gt, seg)[0]

        # converts input and target to numpy arrays
        input, target = convert_to_numpy(input, target)
        if self.use_last_target:
            target = target[:, -1, ...]  # 4D
        else:
            # use 1st target channel
            target = target[:, 0, ...]  # 4D

        # ensure target is of integer type
        target = target.astype(np.int)

        per_batch_arand = []
        for _input, _target in zip(input, target):
            n_clusters = len(np.unique(_target))
            # skip ARand eval if there is only one label in the patch due to the zero-division error in Arand impl
            # xxx/skimage/metrics/_adapted_rand_error.py:70: RuntimeWarning: invalid value encountered in double_scalars
            # precision = sum_p_ij2 / sum_a2
            logger.info(f'Number of ground truth clusters: {n_clusters}')
            if n_clusters == 1:
                logger.info('Skipping ARandError computation: only 1 label present in the ground truth')
                per_batch_arand.append(0.)
                continue

            # convert _input to segmentation CDHW
            segm = self.input_to_segm(_input)
            assert segm.ndim == 4

            # compute per channel arand and return the minimum value
            per_channel_arand = [_arand_err(_target, channel_segm) for channel_segm in segm]
            logger.info(f'Min ARand for channel: {np.argmin(per_channel_arand)}')
            per_batch_arand.append(np.min(per_channel_arand))

        # return mean arand error
        mean_arand = torch.mean(torch.tensor(per_batch_arand))
        logger.info(f'ARand: {mean_arand.item()}')
        return mean_arand

    def input_to_segm(self, input):
        """
        Converts input tensor (output from the network) to the segmentation image. E.g. if the input is the boundary
        pmaps then one option would be to threshold it and run connected components in order to return the segmentation.

        :param input: 4D tensor (CDHW)
        :return: segmentation volume either 4D (segmentation per channel)
        """
        # by deafult assume that input is a segmentation volume itself
        return input

class BoundaryAdaptedRandError(AdaptedRandError):
    """
    Compute ARand between the input boundary map and target segmentation.
    Boundary map is thresholded, and connected components is run to get the predicted segmentation
    """

    def __init__(self, thresholds=None, use_last_target=True, input_channel=None, invert_pmaps=True,
                 save_plots=False, plots_dir='.', **kwargs):
        super().__init__(use_last_target=use_last_target, save_plots=save_plots, plots_dir=plots_dir, **kwargs)
        if thresholds is None:
            thresholds = [0.3, 0.4, 0.5, 0.6]
        assert isinstance(thresholds, list)
        self.thresholds = thresholds
        self.input_channel = input_channel
        self.invert_pmaps = invert_pmaps

    def input_to_segm(self, input):
        if self.input_channel is not None:
            input = np.expand_dims(input[self.input_channel], axis=0)

        segs = []
        for predictions in input:
            for th in self.thresholds:
                # threshold probability maps
                predictions = predictions > th

                if self.invert_pmaps:
                    # for connected component analysis we need to treat boundary signal as background
                    # assign 0-label to boundary mask
                    predictions = np.logical_not(predictions)

                predictions = predictions.astype(np.uint8)
                # run connected components on the predicted mask; consider only 1-connectivity
                seg = measure.label(predictions, background=0, connectivity=1)
                segs.append(seg)

        return np.stack(segs)

class GenericAdaptedRandError(AdaptedRandError):
    def __init__(self, input_channels, thresholds=None, use_last_target=True, invert_channels=None, **kwargs):

        super().__init__(use_last_target=use_last_target, **kwargs)
        assert isinstance(input_channels, list) or isinstance(input_channels, tuple)
        self.input_channels = input_channels
        if thresholds is None:
            thresholds = [0.3, 0.4, 0.5, 0.6]
        assert isinstance(thresholds, list)
        self.thresholds = thresholds
        if invert_channels is None:
            invert_channels = []
        self.invert_channels = invert_channels

    def input_to_segm(self, input):
        # pick only the channels specified in the input_channels
        results = []
        for i in self.input_channels:
            c = input[i]
            # invert channel if necessary
            if i in self.invert_channels:
                c = 1 - c
            results.append(c)

        input = np.stack(results)

        segs = []
        for predictions in input:
            for th in self.thresholds:
                # run connected components on the predicted mask; consider only 1-connectivity
                seg = measure.label((predictions > th).astype(np.uint8), background=0, connectivity=1)
                segs.append(seg)

        return np.stack(segs)

class GenericAveragePrecision:
    def __init__(self, min_instance_size=None, use_last_target=False, metric='ap', **kwargs):
        self.min_instance_size = min_instance_size
        self.use_last_target = use_last_target
        assert metric in ['ap', 'acc']
        if metric == 'ap':
            # use AveragePrecision
            self.metric = AveragePrecision()
        else:
            # use Accuracy at 0.5 IoU
            self.metric = Accuracy(iou_threshold=0.5)

    def __call__(self, input, target):
        if target.dim() == 5:
            if self.use_last_target:
                target = target[:, -1, ...]  # 4D
            else:
                # use 1st target channel
                target = target[:, 0, ...]  # 4D

        input1 = input2 = input
        multi_head = isinstance(input, tuple)
        if multi_head:
            input1, input2 = input

        input1, input2, target = convert_to_numpy(input1, input2, target)

        batch_aps = []
        i_batch = 0
        # iterate over the batch
        for inp1, inp2, tar in zip(input1, input2, target):
            if multi_head:
                inp = (inp1, inp2)
            else:
                inp = inp1

            segs = self.input_to_seg(inp, tar)  # expects 4D
            assert segs.ndim == 4
            # convert target to seg
            tar = self.target_to_seg(tar)

            # filter small instances if necessary
            tar = self._filter_instances(tar)

            # compute average precision per channel
            segs_aps = [self.metric(self._filter_instances(seg), tar) for seg in segs]

            logger.info(f'Batch: {i_batch}. Max Average Precision for channel: {np.argmax(segs_aps)}')
            # save max AP
            batch_aps.append(np.max(segs_aps))
            i_batch += 1

        return torch.tensor(batch_aps).mean()

    def _filter_instances(self, input):
        """
        Filters instances smaller than 'min_instance_size' by overriding them with 0-index
        :param input: input instance segmentation
        """
        if self.min_instance_size is not None:
            labels, counts = np.unique(input, return_counts=True)
            for label, count in zip(labels, counts):
                if count < self.min_instance_size:
                    input[input == label] = 0
        return input

    def input_to_seg(self, input, target=None):
        raise NotImplementedError

    def target_to_seg(self, target):
        return target

class BlobsAveragePrecision(GenericAveragePrecision):
    """
    Computes Average Precision given foreground prediction and ground truth instance segmentation.
    """

    def __init__(self, thresholds=None, metric='ap', min_instance_size=None, input_channel=0, **kwargs):
        super().__init__(min_instance_size=min_instance_size, use_last_target=True, metric=metric)
        if thresholds is None:
            thresholds = [0.4, 0.5, 0.6, 0.7, 0.8]
        assert isinstance(thresholds, list)
        self.thresholds = thresholds
        self.input_channel = input_channel

    def input_to_seg(self, input, target=None):
        input = input[self.input_channel]
        segs = []
        for th in self.thresholds:
            # threshold and run connected components
            mask = (input > th).astype(np.uint8)
            seg = measure.label(mask, background=0, connectivity=1)
            segs.append(seg)
        return np.stack(segs)

class BlobsBoundaryAveragePrecision(GenericAveragePrecision):
    """
    Computes Average Precision given foreground prediction, boundary prediction and ground truth instance segmentation.
    Segmentation mask is computed as (P_mask - P_boundary) > th followed by a connected component
    """

    def __init__(self, thresholds=None, metric='ap', min_instance_size=None, **kwargs):
        super().__init__(min_instance_size=min_instance_size, use_last_target=True, metric=metric)
        if thresholds is None:
            thresholds = [0.3, 0.4, 0.5, 0.6, 0.7]
        assert isinstance(thresholds, list)
        self.thresholds = thresholds

    def input_to_seg(self, input, target=None):
        # input = P_mask - P_boundary
        input = input[0] - input[1]
        segs = []
        for th in self.thresholds:
            # threshold and run connected components
            mask = (input > th).astype(np.uint8)
            seg = measure.label(mask, background=0, connectivity=1)
            segs.append(seg)
        return np.stack(segs)

class BoundaryAveragePrecision(GenericAveragePrecision):
    """
    Computes Average Precision given boundary prediction and ground truth instance segmentation.
    """

    def __init__(self, thresholds=None, min_instance_size=None, input_channel=0, **kwargs):
        super().__init__(min_instance_size=min_instance_size, use_last_target=True)
        if thresholds is None:
            thresholds = [0.3, 0.4, 0.5, 0.6]
        assert isinstance(thresholds, list)
        self.thresholds = thresholds
        self.input_channel = input_channel

    def input_to_seg(self, input, target=None):
        input = input[self.input_channel]
        segs = []
        for th in self.thresholds:
            seg = measure.label(np.logical_not(input > th).astype(np.uint8), background=0, connectivity=1)
            segs.append(seg)
        return np.stack(segs)

class PSNR:
    """
    Computes Peak Signal to Noise Ratio. Use e.g. as an eval metric for denoising task
    """

    def __init__(self, **kwargs):
        pass

    def __call__(self, input, target):
        input, target = convert_to_numpy(input, target)
        return peak_signal_noise_ratio(target, input)

def get_evaluation_metric(config):
    """
    Returns the evaluation metric function based on provided configuration
    :param config: (dict) a top level configuration object containing the 'eval_metric' key
    :return: an instance of the evaluation metric
    """

    def _metric_class(class_name):
        m = importlib.import_module('pytorch3dunet.unet3d.metrics')
        clazz = getattr(m, class_name)
        return clazz

    assert 'eval_metric' in config, 'Could not find evaluation metric configuration'
    metric_config = config['eval_metric']
    metric_class = _metric_class(metric_config['name'])
    return metric_class(**metric_config)


text:
import json
import os
import time
from gensim.models import Word2Vec
import numpy as np
from sklearn.model_selection import train_test_split
from torch import optim
from torch.optim.lr_scheduler import ReduceLROnPlateau
from tqdm import trange
from model_rnn import *

use_cuda = torch.cuda.is_available()
FloatTensor = torch.cuda.FloatTensor if use_cuda else torch.FloatTensor
LongTensor = torch.cuda.LongTensor if use_cuda else torch.LongTensor
Tensor = FloatTensor

# poetry_idx

poetry_data = './data/sample.json'
with open(poetry_data, 'rb') as f:
    poetry_data = json.load(f)

poetry_data = [x for x in poetry_data if len(x[0])>1]  # remove empty data

train, dev = train_test_split(poetry_data, test_size=0.001, train_size=0.99, random_state=4)
print("Length of train {} and length of dev {}".format(len(train), len(dev)))

weights = np.load('./model_results/word2vec_skipgram_weights.pk.npy')
embed_weights = torch.from_numpy(weights)

max_topic_len = 20
max_length = 100
batch_size = 256
train_data = prepare_data_loader(train,max_topic_len, max_length, batch_size, shuffle=True)
dev_data = prepare_data_loader(dev, max_topic_len, max_length, batch_size=10, shuffle=False)

rnn_hidden_dim = 300
rnn_layer = 2
rnn_bidre = False  # setting it True seems to break the data structure
rnn_dropout = 0
dense_dim = 300
dense_h_dropout = 0.5

model = PoetryRNN(embed_weights, max_topic_len,
                  embed_weights,rnn_hidden_dim, rnn_layer,
                  rnn_bidre, rnn_dropout,
                  dense_dim, dense_h_dropout,
                  freeze_embed=True)

model_check_point = "./model_results/fix_embed_model_check_point.pk"
if os.path.isfile(model_check_point):
    model.load_state_dict(torch.load(model_check_point, map_location='cpu'))
    model.train(True)
    print("Load previous training parameters.")

if use_cuda:
    model = model.cuda()
    print("Dump to cuda")

model.train(True)

learning_rate = 5e-3
loss_fn = nn.NLLLoss()
model_params = list(filter(lambda p: p.requires_grad, model.parameters()))
optimizer = optim.Adam(model_params, lr=learning_rate, amsgrad=True,weight_decay=0)

if os.path.isfile("check_point_optim.pkl"):
    optimizer.load_state_dict(torch.load("check_point_optim.pkl"))
    print("Load previous optimizer.")

lr_scheme = ReduceLROnPlateau(optimizer, mode='min', factor=0.4, patience=10, min_lr=1e-7,
                              verbose=True)
counter = 0
all_loss = []
t_range = trange(1)
for iteration in t_range:
    print("starting...")
    losses = []
    for batch in train_data:
        batch_topic, topic_lens, batch_x, batch_y, sorted_lens = sort_batches(batch)
        model.zero_grad()
        target_score, hidden = model(batch_topic, topic_lens, batch_x, sorted_lens)        # batch*seq_len x vocab dim
        loss = loss_fn(target_score, batch_y)

        loss.backward()
        nn.utils.clip_grad_norm_(model.parameters(), max_norm=5)
        optimizer.step()
        # lr_scheme.step(loss.data[0])
        print("Current batch {}, loss {}".format(counter, loss.item()))
        losses.append(loss.item())
        counter += 1
        all_loss.append(loss.item())

    one_ite_loss = np.mean(losses)
    lr_scheme.step(one_ite_loss)
    print("One iteration loss {:.3f}".format(one_ite_loss))

torch.save(model.state_dict(), 'model_check_point.pk')
torch.save(model, 'rnn_model.pk')
torch.save(optimizer.state_dict(), 'optimizer_check_point.pk')
torch.save(all_loss, 'all_loss_check_point.pk')

import torch
a = torch.randn(10,5,2)
lengths = [5,4,3,2,1,1,2,3,4,5]
for i in range(len(lengths)):
    if lengths[i] < a[i,:].size(0):
        a[i, lengths[i]:, :] = 0

max_stars_repo_path: Negative Sampling and LSTM/train_rnn.py

max_issues_repo_path: Negative Sampling and LSTM/train_rnn.py

max_forks_repo_path: Negative Sampling and LSTM/train_rnn.py


text:
from keras import models
from keras import callbacks
from keras.engine import input_layer
from keras_preprocessing import image
import numpy as np
import matplotlib.pyplot as plt
import os
import math
import shutil
import glob
import tensorflow
import scipy

# Number of the Elements in Specific Cases (0-LGG 1-HGG)

ROOT_DIR = "145"
number_of_images = {}
for dir in os.listdir(ROOT_DIR):
    number_of_images[dir] = len(os.listdir(os.path.join(ROOT_DIR, dir)))
number_of_images.items()
#len(os.listdir("145"))

def dataFolder(p,split):

#########################
# Creating Train Folder #
#########################

    if not os.path.exists("./"+p):
        os.mkdir("./"+p)
        for dir in os.listdir(ROOT_DIR):
            os.makedirs("./"+p+"/"+dir)
            for img in np.random.choice(a=os.listdir(os.path.join(ROOT_DIR,dir)),size=(math.floor(split*number_of_images[dir])-2),replace=False):
                O=os.path.join(ROOT_DIR,dir,img)
                D=os.path.join("./"+p,dir)
                shutil.copy(O,D)
                os.remove(O)
    else:
        print(f"{p}Folder exists")
dataFolder("train",0.7)
dataFolder("val",0.15)
dataFolder("test",0.15)

#############
#MODEL BUILD#
#############

from keras.layers import Conv2D,MaxPool2D,Dropout,Flatten,Dense,BatchNormalization,GlobalAvgPool2D
from keras.models import Sequential
from keras.preprocessing.image import ImageDataGenerator
import keras

#Convolution Neural Networks Model

model=Sequential()

model.add(Conv2D(filters=16,kernel_size=(3,3),activation='relu',input_shape =(240,240,3)))
model.add(MaxPool2D(pool_size=(2,2))) #removing unnecessary features

model.add(Conv2D(filters=32,kernel_size=(3,3),activation='relu'))
model.add(MaxPool2D(pool_size=(2,2))) #removing unnecessary features

model.add(Conv2D(filters=64,kernel_size=(3,3),activation='relu'))
model.add(MaxPool2D(pool_size=(2,2))) #removing unnecessary features

model.add(Conv2D(filters=128,kernel_size=(3,3),activation='relu'))
model.add(MaxPool2D(pool_size=(2,2))) #removing unnecessary features

model.add(Dropout(rate=0.25)) # to avoid overfitting

model.add(Flatten()) #reduce to 1 dimension
model.add(Dense(units=64,activation='relu')) # multilayer perceptron with 64 outputs
model.add(Dropout(rate=0.25)) # to avoid overfitting
model.add(Dense(units=16,activation='relu')) # multilayer perceptron with 64 outputs
model.add(Dropout(rate=0.25)) # to avoid overfitting
model.add(Dense(units=1,activation='sigmoid')) # multilayer perceptron with 1 output
model.summary()
model.compile(optimizer='adam',loss=keras.losses.binary_crossentropy,metrics=['accuracy'])
print(model.get_weights())

####################
# Data Preparation #
####################

from keras.preprocessing.image import ImageDataGenerator
def preprocessingImages1(path):
    """
    input: Path
    output: Preprocessed Images
    """

    image_data=ImageDataGenerator(zoom_range=0.2,shear_range=0.2,rescale=1/255,horizontal_flip=True) # data augmentation
    image=image_data.flow_from_directory(directory=path,target_size=(240,240),batch_size=16,class_mode='binary')
    return image

path="train"
train_data=preprocessingImages1(path)
train_data.class_indices
def preprocessingImages2(path):
    """
    input: Path
    output: Preprocessed Images
    """

    image_data=ImageDataGenerator(rescale=1/255)
    image=image_data.flow_from_directory(directory=path,target_size=(240,240),batch_size=16,class_mode='binary')
    return image
path="test"
test_data=preprocessingImages2(path)
path="val"
val_data=preprocessingImages2(path)

###########################################
# Early Stop and Checking Model Situation #
###########################################

from keras.callbacks import ModelCheckpoint,EarlyStopping

# Early Stopping #

es=EarlyStopping(monitor="val_accuracy",min_delta=0.001,patience=10,verbose=1,mode='auto')

# Early Stopping Checkpoint #

mc=ModelCheckpoint(monitor="val_accuracy",filepath="./bestmodel.h5",verbose=1,save_best_only=True,mode='auto')
cd=[es,mc]

# Model Training #

hs= model.fit(x=train_data,steps_per_epoch=5,epochs=30,verbose=1,validation_data=val_data,validation_steps=16,callbacks=cd)

# Model Graphical Transformation #

h=hs.history
h.keys()

# Accuracy Graphing #

import matplotlib.pyplot as plt
plt.plot(h['accuracy'])
plt.plot(h['val_accuracy'], c="red")
plt.title("accuracy(blue) vs validation accuracy(red)")
plt.show()

# Loss Graphing #

plt.plot(h['loss'])
plt.plot(h['val_loss'], c="red")
plt.title("loss(blue) vs validation loss(red)")
plt.show()

# Model Accuracy #

from keras.models import load_model
model=load_model('bestmodel.h5')
acc=model.evaluate(test_data)[1]
print(f"The Accuracy of Our Model by function is {acc*100}%")

# Own Accuracy Function #

from keras.preprocessing.image import load_img, img_to_array
Root_direc="test\LGG"
LGG_err_counter=0
LGG_loop_counter=0
#train_data.class_indices
for dir in os.listdir(Root_direc):
    LGG_loop_counter=LGG_loop_counter+1
    path ="test\LGG"+"/"+dir
    img=load_img(path,target_size=(240,240))
    input_arr=img_to_array(img)/255
    #input_arr.shape
    input_arr=np.expand_dims(input_arr,axis=0)
    pred= ((model.predict(input_arr)>0.5).astype("int32"))[0][0]
    if pred == 0:
        #print("Patient is having  HGG")
        LGG_err_counter=LGG_err_counter+1
    else:
        #print("Patient is  having LGG")
        continue

Root_direc="test\HGG"
HGG_err_counter=0
HGG_loop_counter=0
for dir in os.listdir(Root_direc):
    HGG_loop_counter=HGG_loop_counter+1
    path ="test\HGG"+"/"+dir
    img=load_img(path,target_size=(240,240))
    input_arr=img_to_array(img)/255
    #input_arr.shape
    input_arr=np.expand_dims(input_arr,axis=0)
    pred= ((model.predict(input_arr)>0.5).astype("int32"))[0][0]
    if pred == 0:
        #print("Patient is having  HGG")
        continue
    else:
        #print("Patient is  having LGG")
        HGG_err_counter=HGG_err_counter+1

# Accuracy Calculation

LGG_accuracy=1-(LGG_err_counter/LGG_loop_counter)
HGG_accuracy=1-(HGG_err_counter/HGG_loop_counter)
print(f"LGG Accuracy:{LGG_accuracy}%")
print(f"HGG Accuracy:{HGG_accuracy}%")
print(f"Model Accuracy(own function):{((LGG_accuracy+HGG_accuracy)/2)*100}%")

# Manual Testing Part #
path ="145\HGG\BraTS20_Training_226_t1ce_z78.png"
#path ="145\HGG\Y14.jpg"
#path ="145\LGG\BraTS20_Training_296_t1ce_z84.png"
#path ="test\HGG\BraTS20_Training_001_t1ce_z76.png"
img=load_img(path,target_size=(240,240))
input_arr=img_to_array(img)/255
plt.imshow(input_arr)
input_arr.shape
input_arr=np.expand_dims(input_arr,axis=0)
pred= ((model.predict(input_arr)>0.5).astype("int32"))[0][0]
#train_data.class_indices
if pred == 0:
    print("Patient is having  HGG")
else:
    print("Patient is  having LGG")

# Plotting #

temp=0
temp1=0
value=[]
validation=[]
for epoch in range(len(hs.epoch)):
    temp=temp+h['accuracy'][epoch]
    temp1=temp1+h['val_accuracy'][epoch]
    value.append(temp/(epoch+1))
    validation.append(temp1/(epoch+1))

plt.plot(value)
plt.plot(validation, c="red")
plt.title("Accuracy(blue) vs Validation Accuracy(red)")
plt.show()
#print(model.get_weights())
#print((model.predict(input_arr)>0.5))
#model.summary()

reason: import numpy,import scipy