text:
# ---------------------------------------------------------------
# Copyright (c) 2020, NVIDIA CORPORATION. All rights reserved.
#
# This work is licensed under the NVIDIA Source Code License
# for NVAE. To view a copy of this license, see the LICENSE file.
# ---------------------------------------------------------------

import torch
import torch.nn as nn
import torch.nn.functional as F
import numpy as np

from .utils import one_hot

@torch.jit.script
def soft_clamp5(x: torch.Tensor):
    return x.div_(5.).tanh_().mul(5.)    #  5. * torch.tanh(x / 5.) <--> soft differentiable clamp between [-5, 5]

@torch.jit.script
def sample_normal_jit(mu, sigma):
    eps = mu.mul(0).normal_()
    z = eps.mul_(sigma).add_(mu)
    return z, eps

class Normal:
    def __init__(self, mu, log_sigma, temp=1.):
        self.mu = soft_clamp5(mu)
        log_sigma = soft_clamp5(log_sigma)
        self.sigma = torch.exp(log_sigma) + 1e-2      # we don't need this after soft clamp
        if temp != 1.:
            self.sigma *= temp

    def sample(self):
        return sample_normal_jit(self.mu, self.sigma)

    def sample_given_eps(self, eps):
        return eps * self.sigma + self.mu

    def log_p(self, samples):
        normalized_samples = (samples - self.mu) / self.sigma
        log_p = - 0.5 * normalized_samples * normalized_samples - 0.5 * np.log(2 * np.pi) - torch.log(self.sigma)
        return log_p

    def kl(self, normal_dist):
        term1 = (self.mu - normal_dist.mu) / normal_dist.sigma
        term2 = self.sigma / normal_dist.sigma

        return 0.5 * (term1 * term1 + term2 * term2) - 0.5 - torch.log(term2)

class DiscLogistic:
    def __init__(self, param):
        B, C, H, W = param.size()
        self.num_c = C // 2
        self.means = param[:, :self.num_c, :, :]                              # B, 3, H, W
        self.log_scales = torch.clamp(param[:, self.num_c:, :, :], min=-8.0)  # B, 3, H, W

    def log_prob(self, samples):
        assert torch.max(samples) <= 1.0 and torch.min(samples) >= 0.0
        # convert samples to be in [-1, 1]
        samples = 2 * samples - 1.0

        B, C, H, W = samples.size()
        assert C == self.num_c

        centered = samples - self.means                                         # B, 3, H, W
        inv_stdv = torch.exp(- self.log_scales)
        plus_in = inv_stdv * (centered + 1. / 255.)
        cdf_plus = torch.sigmoid(plus_in)
        min_in = inv_stdv * (centered - 1. / 255.)
        cdf_min = torch.sigmoid(min_in)
        log_cdf_plus = plus_in - F.softplus(plus_in)
        log_one_minus_cdf_min = - F.softplus(min_in)
        cdf_delta = cdf_plus - cdf_min
        mid_in = inv_stdv * centered
        log_pdf_mid = mid_in - self.log_scales - 2. * F.softplus(mid_in)

        log_prob_mid_safe = torch.where(cdf_delta > 1e-5,
                                        torch.log(torch.clamp(cdf_delta, min=1e-10)),
                                        log_pdf_mid - np.log(127.5))
        # woow the original implementation uses samples > 0.999, this ignores the largest possible pixel value (255)
        # which is mapped to 0.9922
        log_probs = torch.where(samples < -0.999, log_cdf_plus, torch.where(samples > 0.99, log_one_minus_cdf_min,
                                                                            log_prob_mid_safe))   # B, 3, H, W

        return log_probs

    def sample(self):
        u = torch.Tensor(self.means.size()).uniform_(1e-5, 1. - 1e-5).cuda()                        # B, 3, H, W
        x = self.means + torch.exp(self.log_scales) * (torch.log(u) - torch.log(1. - u))            # B, 3, H, W
        x = torch.clamp(x, -1, 1.)
        x = x / 2. + 0.5
        return x

class DiscMixLogistic:
    def __init__(self, param, num_mix=10, num_bits=8):
        B, C, H, W = param.size()
        self.num_mix = num_mix
        self.logit_probs = param[:, :num_mix, :, :]                                   # B, M, H, W
        l = param[:, num_mix:, :, :].view(B, 3, 3 * num_mix, H, W)                    # B, 3, 3 * M, H, W
        self.means = l[:, :, :num_mix, :, :]                                          # B, 3, M, H, W
        self.log_scales = torch.clamp(l[:, :, num_mix:2 * num_mix, :, :], min=-7.0)   # B, 3, M, H, W
        self.coeffs = torch.tanh(l[:, :, 2 * num_mix:3 * num_mix, :, :])              # B, 3, M, H, W
        self.max_val = 2. ** num_bits - 1

    def log_prob(self, samples):
        assert torch.max(samples) <= 1.0 and torch.min(samples) >= 0.0
        # convert samples to be in [-1, 1]
        samples = 2 * samples - 1.0

        B, C, H, W = samples.size()
        assert C == 3, 'only RGB images are considered.'

        samples = samples.unsqueeze(4)                                                  # B, 3, H , W
        samples = samples.expand(-1, -1, -1, -1, self.num_mix).permute(0, 1, 4, 2, 3)   # B, 3, M, H, W
        mean1 = self.means[:, 0, :, :, :]                                               # B, M, H, W
        mean2 = self.means[:, 1, :, :, :] + \
                self.coeffs[:, 0, :, :, :] * samples[:, 0, :, :, :]                     # B, M, H, W
        mean3 = self.means[:, 2, :, :, :] + \
                self.coeffs[:, 1, :, :, :] * samples[:, 0, :, :, :] + \
                self.coeffs[:, 2, :, :, :] * samples[:, 1, :, :, :]                     # B, M, H, W

        mean1 = mean1.unsqueeze(1)                          # B, 1, M, H, W
        mean2 = mean2.unsqueeze(1)                          # B, 1, M, H, W
        mean3 = mean3.unsqueeze(1)                          # B, 1, M, H, W
        means = torch.cat([mean1, mean2, mean3], dim=1)     # B, 3, M, H, W
        centered = samples - means                          # B, 3, M, H, W

        inv_stdv = torch.exp(- self.log_scales)
        plus_in = inv_stdv * (centered + 1. / self.max_val)
        cdf_plus = torch.sigmoid(plus_in)
        min_in = inv_stdv * (centered - 1. / self.max_val)
        cdf_min = torch.sigmoid(min_in)
        log_cdf_plus = plus_in - F.softplus(plus_in)
        log_one_minus_cdf_min = - F.softplus(min_in)
        cdf_delta = cdf_plus - cdf_min
        mid_in = inv_stdv * centered
        log_pdf_mid = mid_in - self.log_scales - 2. * F.softplus(mid_in)

        log_prob_mid_safe = torch.where(cdf_delta > 1e-5,
                                        torch.log(torch.clamp(cdf_delta, min=1e-10)),
                                        log_pdf_mid - np.log(self.max_val / 2))
        # the original implementation uses samples > 0.999, this ignores the largest possible pixel value (255)
        # which is mapped to 0.9922
        log_probs = torch.where(samples < -0.999, log_cdf_plus, torch.where(samples > 0.99, log_one_minus_cdf_min,
                                                                            log_prob_mid_safe))   # B, 3, M, H, W

        log_probs = torch.sum(log_probs, 1) + F.log_softmax(self.logit_probs, dim=1)  # B, M, H, W
        return torch.logsumexp(log_probs, dim=1)                                      # B, H, W

    def sample(self, t=1.):
        gumbel = -torch.log(- torch.log(torch.Tensor(self.logit_probs.size()).uniform_(1e-5, 1. - 1e-5).cuda()))  # B, M, H, W
        sel = one_hot(torch.argmax(self.logit_probs / t + gumbel, 1), self.num_mix, dim=1)          # B, M, H, W
        sel = sel.unsqueeze(1)                                                                 # B, 1, M, H, W

        # select logistic parameters
        means = torch.sum(self.means * sel, dim=2)                                             # B, 3, H, W
        log_scales = torch.sum(self.log_scales * sel, dim=2)                                   # B, 3, H, W
        coeffs = torch.sum(self.coeffs * sel, dim=2)                                           # B, 3, H, W

        # cells from logistic & clip to interval
        # we don't actually round to the nearest 8bit value when sampling
        u = torch.Tensor(means.size()).uniform_(1e-5, 1. - 1e-5).cuda()                        # B, 3, H, W
        x = means + torch.exp(log_scales) / t * (torch.log(u) - torch.log(1. - u))             # B, 3, H, W

        x0 = torch.clamp(x[:, 0, :, :], -1, 1.)                                                # B, H, W
        x1 = torch.clamp(x[:, 1, :, :] + coeffs[:, 0, :, :] * x0, -1, 1)                       # B, H, W
        x2 = torch.clamp(x[:, 2, :, :] + coeffs[:, 1, :, :] * x0 + coeffs[:, 2, :, :] * x1, -1, 1)  # B, H, W

        x0 = x0.unsqueeze(1)
        x1 = x1.unsqueeze(1)
        x2 = x2.unsqueeze(1)

        x = torch.cat([x0, x1, x2], 1)
        x = x / 2. + 0.5
        return x


text:
import sys
import re
import numpy as np
import matplotlib.pyplot as plt

# need one argument as input file
if len(sys.argv) != 2:
    print('Usage: %s input_log_file.txt' % sys.argv[0])
    sys.exit(1)

# opening the file
filename = sys.argv[1]
try:
    f = open(filename, 'r')
except:
    print('Please specify a valid file')
    exit(1)

# regular expression to match
pattern = '^([0-9]+)\s\+([0-9]+)\sms\s([0-9a-z]+)\s\[([^@]+)\@([^:]+):([^]]+)]\s(.*)$'
program = re.compile(pattern)
# 1 timestamp, 2 delta_ms, 3 thread, 4 function, 5 file, 6 line, 7 message

# transition array
execution_time = {}

# begin deltas
begin_times = {}

# ignoring first lines
ignore_lines = 0

# begin time for a task
block_begin = {}

# pattern of messages to analyze
analyzer_pattern = re.compile("^(\[ANALYZE]) ([^ ]+) (.*)$")

# all tids
all_threads = set()

for i, line in enumerate(f):
    # ignoring first lines
    if i < ignore_lines: continue

    # removing newline
    line = line.strip()

    # trying to match with re
    matching = program.split(line)

    # 9 is a magic number (see definition of pattern)
    if matching and len(matching) == 9:
        # parsing the array
        _, timestamp, delta_ms, thread, function, file_, line, message, _ = matching

        # splitting over analyze pattern
        split = analyzer_pattern.split(message)

        if split and len(split) == 5:
            _, _, status, block, _ = split
 #           print(timestamp, thread, block, status)
            if status == 'Begin':
                if block not in execution_time.keys(): execution_time[block] = []
                block_begin[block] = (timestamp, thread)
                if block not in begin_times: begin_times[block] = []
                begin_times[block].append(int(timestamp))
            elif status == 'End':
                if block not in block_begin.keys():
                    print('Stray block end: %s %s %s' % (block, timestamp, thread))
                    continue
                if thread != block_begin[block][1]:
                    print('Multistart for block ' + block + ' tid1 ' + block_begin[block][1] + ' tid2 ' + thread)
                    continue
                if block not in execution_time.keys(): execution_time[block] = []
                execution_time[block].append(int(timestamp) - int(block_begin[block][0]))
                all_threads.add(thread)
            else:
                assert(False)

# begin deltas for blocks
begin_deltas = {x: list(filter(lambda x: x < 1000, list(np.array(y[1:]) - np.array(y[:-1])))) for x, y in begin_times.items()}

# printing all successors
#print(execution_time)

# all tids
print('Used threads:')
print(all_threads)

# all blocks
print('Blocks:')
print(block_begin.keys())

# show only these blocks
#block_filter = {'Uchiya': 'Detector', 'QtCameraMain': 'Camera', 'HueThresholdManual': 'HSV threshold', 'FrameConvert': 'yuv2rgb'}
block_filter = {x: x for x in block_begin.keys()}

# printing the result
all_arrays = []
all_labels = []
all_blocks = []
for block, arr in execution_time.items():
    if block not in block_filter.keys(): continue
    block = block_filter[block]
    if len(arr) > 0:
        all_arrays.append(arr)
        all_labels.append('%s $\\rightarrow$\n$%.1f$ FPS, ms: $%.2f\pm %.2f$ med $%.2f$' % (block,  1000./ np.mean(arr), np.mean(arr), np.std(arr), np.median(arr)))
        all_blocks.append(block)
    arr1 = begin_deltas[block]
    all_arrays.append(arr1)
    all_labels.append('%s $\circlearrowright$\n$%.1f$ FPS, ms: $%.2f\pm %.2f$ med $%.2f$' % (block,  1000./ np.mean(arr1), np.mean(arr1), np.std(arr1), np.median(arr1)))
    all_blocks.append(block + " Delta")

for label in all_labels:
    print(' '.join(label.split()))

plt.figure(figsize=(15, 15))
plt.title('Processing time for thread loop')
plt.xlabel('ms')
#plt.yticks(rotation=25)
plt.boxplot(all_arrays, labels = all_labels, vert = False)
plt.savefig(filename + '.png', bbox_inches='tight')
print('Output: ' + filename + '.png')
plt.clf()

for arr, block in zip(all_arrays, all_blocks):
    plt.figure(figsize=(10, 10))
    plt.hist(arr, 20, facecolor='g', alpha=0.75)
    plt.xlabel('ms')
    plt.ylabel('Count')
    plt.title('Histogram of ' + block + ' time')
    plt.grid(True)
    fn_local = filename + '_' + block + '.png'
    plt.savefig(fn_local, bbox_inches='tight')
    print('Block ' + block + ': ' + fn_local)
    plt.clf()

ext: py

max_stars_repo_name: chili-epfl/qml-ar

max_issues_repo_name: chili-epfl/qml-ar

max_forks_repo_name: chili-epfl/qml-ar

reason: import numpy


text:
import pyflo.component
import suspect

from collections import deque
import numpy
import numbers

class BeginGroup:
    pass

class EndGroup:
    pass

class SVDChannelWeights(pyflo.component.Component):

    def __init__(self, metadata):
        inports = {
            "in": {},
        }
        outports = {
            "weights": {}
        }
        super().__init__(inports, outports)
        self.inports["in"].on('data', self.calculate_channel_weights)

    def calculate_channel_weights(self, data):
        weights = suspect.processing.channel_combination.svd_weighting(data)
        self.outports["weights"].send_data(weights)

class WeightedAverage(pyflo.component.Component):
    axis = None

    def __init__(self, metadata):
        inports = {
            "in": {},
            "weights": {},
            "axis": {},
        }
        outports = {
            "out": {}
        }
        super().__init__(inports, outports)
        self.inports["in"].on('data', self.data)
        self.queued_data = deque()
        self.inports["weights"].on('data', self.weights)
        self.queued_weights = deque()
        self.inports["axis"].on('data', self.set_axis)

    def data(self, data):
        self.queued_data.append(data)
        self.combine()

    def weights(self, data):
        self.queued_weights.append(data)
        self.combine()

    def combine(self):
        # options for combining data
        # 1: axis not attached, weights not attached, data in queue
        # 2: axis set, weights not attached, data in queue
        # 3: axis not attached, weights in queue, data in queue
        # 4: axis set, weights in queue, data in queue
        if self.inports["axis"].attached():
            if self.axis is None:
                return
            else:
                axis = int(self.axis)
        else:
            axis = 0
        if self.inports["weights"].attached():
            while len(self.queued_data) > 0 and len(self.queued_weights) > 0:
                data_to_combine = self.queued_data.popleft()
                weights = self.queued_weights.popleft()
                combined_data = data_to_combine.inherit(numpy.average(data_to_combine, axis=axis, weights=weights))
                self.outports["out"].send_data(combined_data)
        else:
            while len(self.queued_data) > 0:
                data_to_combine = self.queued_data.popleft()
                combined_data = data_to_combine.inherit(numpy.average(data_to_combine, axis=axis))
                self.outports["out"].send_data(combined_data)

    def set_axis(self, axis):
        self.axis = axis
        self.combine()

class WaterPeakAlignment(pyflo.component.Component):
    def __init__(self, metadata):
        inports = {
            "in": {}
        }
        outports = {
            "shift": {}
        }
        super().__init__(inports, outports)
        self.inports["in"].on("beginGroup", self.outports["shift"].begin_group)
        self.inports["in"].on("data", self.data)
        self.inports["in"].on("endGroup", self.outports["shift"].end_group)

    def data(self, data):
        frequency_shift = suspect.processing.frequency_correction.residual_water_alignment(data)
        self.outports["shift"].send_data(frequency_shift)

class FrequencyShift(pyflo.component.Component):
    def __init__(self, metadata):
        inports = {
            "in": {},
            "shift": {},
        }
        outports = {
            "out": {},
        }
        super().__init__(inports, outports)
        self.inports["in"].on("beginGroup", self.in_begin_group)
        self.inports["in"].on("data", self.data)
        self.inports["in"].on("endGroup", self.in_end_group)
        self.in_queue = deque()
        self.inports["shift"].on("beginGroup", self.shift_begin_group)
        self.inports["shift"].on("data", self.shift)
        self.inports["shift"].on("endGroup", self.shift_end_group)
        self.shift_queue = deque()

    def in_begin_group(self):
        self.in_queue.append(BeginGroup())

    def data(self, data):
        self.in_queue.append(data)
        self.process()

    def in_end_group(self):
        self.in_queue.append(EndGroup())
        self.process()

    def shift_begin_group(self):
        self.shift_queue.append(BeginGroup())

    def shift(self, shift):
        self.shift_queue.append(shift)
        self.process()

    def shift_end_group(self):
        self.shift_queue.append(EndGroup())
        self.process()

    def process(self):
        while len(self.in_queue) and len(self.shift_queue):
            in_val = self.in_queue.popleft()
            shift_val = self.shift_queue.popleft()
            if isinstance(in_val, BeginGroup) and isinstance(shift_val, BeginGroup):
                self.outports["out"].begin_group()
            elif isinstance(in_val, EndGroup) and isinstance(shift_val, EndGroup):
                self.outports["out"].end_group()
            elif type(in_val) == suspect.MRSData and isinstance(shift_val, numbers.Number):
                print("correcting frequency shift {}".format(shift_val))
                transformed_fid = suspect.processing.frequency_correction.transform_fid(in_val, shift_val, 0)
                self.outports["out"].send_data(transformed_fid)
            else:
                raise Exception("in {} and shift {} values didn't match".format(in_val, shift_val))

import matplotlib
matplotlib.use("Agg")
from matplotlib import pyplot

class TempPlot(pyflo.component.Component):
    def __init__(self, metadata):
        inports = {
            "in": {},
        }
        super().__init__(inports, {})
        self.inports["in"].on('data', self.data)
        self.inports["in"].on('beginGroup', self.begin)
        self.inports["in"].on('endGroup', self.end)

    def begin(self):
        pass

    def data(self, data):
        spectrum = numpy.fft.fftshift(numpy.fft.fft(data))
        pyplot.plot(numpy.abs(spectrum))
        pyplot.savefig("/home/ben/sample_plot.pdf")

    def end(self):
        print("saving figure")
        pyplot.savefig("/home/ben/sample_plot.pdf")

max_stars_repo_name: bennyrowland/pyflo-suspect

max_issues_repo_name: bennyrowland/pyflo-suspect

max_forks_repo_name: bennyrowland/pyflo-suspect


text:
"""
Copyright Matt DeMartino (Stravajiaxen)
Licensed under MIT License -- do whatever you want with this, just don't sue me!

This code attempts to solve Project Euler (projecteuler.net)
Problem #26 Reciprocal cycles

A unit fraction contains 1 in the numerator.
The decimal representation of the unit fractions with denominators 2 to 10 are given:

Where 0.1(6) means 0.166666..., and has a 1-digit recurring cycle.
It can be seen that 1/7 has a 6-digit recurring cycle.

Find the value of d < 1000 for which 1/d contains the longest recurring cycle in its
decimal fraction part.
"""
import time

# Based on some code from:
# https://stackoverflow.com/questions/25181642/how-set-numpy-floating-point-accuracy
import sympy as smp
import mpmath as mp

def main():
    mp.mp.dps = 3000  # Computation precision is 1000 digits

    num = mp.mpf("1.0")
    for i in range(1, 1001):

        den = mp.mpf(str(i))
        # print(f"1/{i} = ", num/den)
        val = str(num/den)
        for offset in range(1, 1000):
            for cycle_size in range(1, 1000):
                raise NotImplementedError("I was able to solve this one from my research!")

if __name__ == "__main__":
    start_time = time.time()
    main()
    elapsed_time = time.time() - start_time
    print("Elapsed Time: ", elapsed_time)

max_stars_repo_name: stravajiaxen/project-euler-solutions

max_issues_repo_name: stravajiaxen/project-euler-solutions

max_forks_repo_name: stravajiaxen/project-euler-solutions

reason: import sympy,import mpmath


text:
"""
This program will generate the first step in the DR16Q Superset file chain
from the spAll-v5_13_0.fits file. This input file is ~16Gb, so be careful.
Because the files in use are large, and some functions can take a long time
each function will output its own completed file, which the next
function will load in turn. Have some HDD space free. A lot.

Dependencies
----------
sdss_flagval : a pydl utility
github repository : https://github.com/bradlyke/utilities

Parameters
----------
None

Output
----------
Four catalog files output to the ../data directory
1) DR16Q_bitcut_{date}.fits: The initial spAll file subsampled on QSO bits
2) DR16Q_autoclass_{date}.fits: The catalog with AUTOCLASS_DR14Q populated
3) DR16Q_vdb_zw_{date}.fits: The classified catalog with previous
                             visual inspection data added
4) DR16Q_vdb_rem_{date}.fits: The vdb file with records removed that have
                              bad ZWARNING flags (this is the final file)

"""

from astropy.io import fits
import numpy as np
import sys
import os
import time
from pydl.pydlutils.sdss import sdss_flagval as flagval
import sdss_name_v2 as sdn
import cat_tools as ct
import progressBar as pb
import empty_cat_maker as ecm

#This function will select only the object records that were targeted as quasars
#within the spAll file. This will only include, by its nature, BOSS and eBOSS
#observations.
def cat_cut():
    #BOSS_TARGET1 bit flag names
    bt1_flag_names = np.array(['QSO_CORE','QSO_BONUS','QSO_KNOWN_MIDZ','QSO_KNOWN_LOHIZ',
                    'QSO_NN','QSO_UKIDSS','QSO_KDE_COADD','QSO_LIKE','QSO_FIRST_BOSS',
                    'QSO_KDE','QSO_CORE_MAIN','QSO_BONUS_MAIN','QSO_CORE_ED',
                    'QSO_CORE_LIKE','QSO_KNOWN_SUPPZ'])
    #EBOSS_TARGET0 bit flag names
    et0_flag_names = np.array(['QSO_EBOSS_CORE','QSO_PTF','QSO_REOBS','QSO_EBOSS_KDE',
                    'QSO_EBOSS_FIRST','QSO_BAD_BOSS','QSO_BOSS_TARGET','QSO_SDSS_TARGET',
                    'QSO_KNOWN','SPIDERS_RASS_AGN','SPIDERS_ERASS_AGN','TDSS_A',
                    'TDSS_FES_DE','TDSS_FES_NQHISN','TDSS_FES_MGII','TDSS_FES_VARBAL',
                    'SEQUELS_PTF_VARIABLE'])
    #EBOSS_TARGET1 bit flag names
    et1_flag_names = np.array(['QSO1_VAR_S82','QSO1_EBOSS_CORE','QSO1_PTF','QSO1_REOBS',
                    'QSO1_EBOSS_KDE','QSO1_EBOSS_FIRST','QSO1_BAD_BOSS','QSO_BOSS_TARGET',
                    'QSO_SDSS_TARGET','QSO_KNOWN','TDSS_TARGET','SPIDERS_TARGET','S82X_TILE1',
                    'S82X_TILE2','S82X_TILE3'])
    #EBOSS_TARGET2 bit flag names
    et2_flag_names = np.array(['SPIDERS_RASS_AGN','SPIDERS_ERASS_AGN','SPIDERS_XMMSL_AGN',
                    'TDSS_A','TDSS_FES_DE','TDSS_FES_NQHISN','TDSS_FES_MGII','TDSS_FES_VARBAL',
                    'TDSS_B','TDSS_FES_HYPQSO','TDSS_CP','S82X_BRIGHT_TARGET',
                    'S82X_XMM_TARGET','S82X_WISE_TARGET','S82X_SACLAY_VAR_TARGET',
                    'S82X_SACLAY_BDT_TARGET','S82X_SACLAY_HIZ_TARGET',
                    'S82X_RICHARDS15_PHOTOQSO_TARGET','S82X_PETERS15_COLORVAR_TARGET',
                    'S82X_LSSTZ4_TARGET','S82X_UNWISE_TARGET','S82X_GTRADMZ4_TARGET',
                    'S82X_CLAGN1_TARGET','S82X_CLAGN2_TARGET'])
    #ANCILLARY_TARGET1 bit flag names
    at1_flag_names = np.array(['BLAZGVAR','BLAZR','BLAZXR','BLAZXRSAM','BLAZXRVAR',
                    'XMMBRIGHT','XMMGRIZ','XMMHR','XMMRED','FBQSBAL','LBQSBAL','ODDBAL',
                    'OTBAL','PREVBAL','VARBAL','QSO_AAL','QSO_AALS','QSO_IAL','QSO_RADIO',
                    'QSO_RADIO_AAL','QSO_RADIO_IAL','QSO_NOAALS','QSO_GRI','QSO_HIZ',
                    'QSO_RIZ','BLAZGRFLAT','BLAZGRQSO','BLAZGX','BLAZGXQSO','BLAZGXR',
                    'CXOBRIGHT','CXORED'])
    #ANCILLARY_TARGET2 bit flag names
    at2_flag_names = np.array(['HIZQSO82','HIZQSOIR','KQSO_BOSS','QSO_VAR','QSO_VAR_FPG',
                    'RADIO_2LOBE_QSO','QSO_SUPPZ','QSO_VAR_SDSS','QSO_WISE_SUPP',
                    'QSO_WISE_FULL_SKY','DISKEMITTER_REPEAT','WISE_BOSS_QSO',
                    'QSO_XD_KDE_PAIR','TDSS_PILOT','SPIDERS_PILOT','TDSS_SPIDERS_PILOT',
                    'QSO_VAR_LF','QSO_EBOSS_W3_ADM','XMM_PRIME','XMM_SECOND',
                    'SEQUELS_TARGET','RM_TILE1','RM_TILE2','QSO_DEEP'])

    #Turn the bit flag name arrays into a bit mask value
    bt1 = flagval('BOSS_TARGET1',bt1_flag_names)
    et0 = flagval('EBOSS_TARGET0',et0_flag_names)
    et1 = flagval('EBOSS_TARGET1',et1_flag_names)
    et2 = flagval('EBOSS_TARGET2',et2_flag_names)
    at1 = flagval('ANCILLARY_TARGET1',at1_flag_names)
    at2 = flagval('ANCILLARY_TARGET2',at2_flag_names)

    loc = os.environ.get('BOSS_SPECTRO_REDUX') #Get the location for the spAll file.
    vnum = 'v5_13_0'
    full_file = '{0}/{1}/spAll-{2}.fits'.format(loc,vnum,vnum)
    spfull = fits.open(full_file)[1].data #Load the whole 16Gb file.

    #Cut down to only the records targeted as quasars using bitmasks.
    wq = np.where((spfull['BOSS_TARGET1']&bt1)|(spfull['EBOSS_TARGET0']&et0)|
            (spfull['EBOSS_TARGET1']&et1)|(spfull['EBOSS_TARGET2']&et2)|
            (spfull['ANCILLARY_TARGET1']&at1)|(spfull['ANCILLARY_TARGET2']&at2))[0]

    #If the user wants to track which version of the spAll file was used, this
    #will record that file version in the Primary (HDU0) header.
    fname = 'spAll-{}'.format(vnum)
    prim_hdr = fits.Header()
    prim_hdr['DATA_REL']=(fname,'spAll file used')
    prim_hdu = fits.PrimaryHDU(header=prim_hdr)

    #Output the cutdown file of only those records targeted as quasars.
    data_out = np.array(spfull[wq])
    data_hdu = fits.BinTableHDU.from_columns(data_out)
    data_of = fits.HDUList([prim_hdu,data_hdu])
    file_dtag = time.strftime('%Y%m%d')
    outfile_name = '../data/DR16Q_bitcut_{}'.format(file_dtag)
    #FET means (F)ITS writer (E)rror (T)rap. This will force the file to be written
    #just in case a file already exists with that file name.
    superset_name = ct.fet(data_of,outfile_name)

    return superset_name #To load this output in the next function, return output name.

#################################CLASSIFICATION#################################

################################################################################
#                                                                              #
#  SCRIPT: VI Catalogue Flagging Algorithm                                     #
#  AUTHOR: Brad Lyke, with base algorithm design by Isabelle Paris             #
#                                                                              #
#  PURPOSE: This script will analyze all objects from the latest spAll file    #
#           from SDSS and generate a DRQ flag for visual inspection followup.  #
#           The algorithm pseudocode is as follows:                            #
#           1)If the object 'CLASS' is 'STAR', then label as 'STAR'            #
#           2)If the object 'CLASS' is 'GALAXY' and 'Z' < 1, then 'GALAXY'     #
#           3)If the 'CLASS' is 'GALAXY' and 'Z' >= 1, check spZall next top 4 #
#             model fits and decide:                                           #
#             --If 1+ are 'GALAXY' then 'GALAXY'                               #
#           4)If the 'CLASS' is 'QSO', check spZall next top 4 fits and decide:#
#             --If 2+ are 'STAR' then 'STAR',                                  #
#             --Else If < 2 are 'STAR', ZWARNING=0 then 'QSO'                  #
#             --Else then 'QFLAG'                                              #
#           If, after these, the object has no label (or QFLAG which is an     #
#           internal variable for efficiency), then it will need VI followup.  #
#                                                                              #
#  VARIABLES: sp_cat - The name of the current spAll file to use, passed as a  #
#                      system argument.                                        #
#                                                                              #
#  INPUT: This file expects the input catalogue to only be items from the      #
#           spAll file that have OBJTYPE 'QSO'. Do not use full spAll file.    #
#         Runs from terminal. Command as follows                               #
#            python spAll_class_full.py [spAll filename]                       #
#            eg. python spAll_class_full.py spAll-v5_13_0.fits                 #
#  OUTPUT: -A new fits file with object classifications in AUTOCLASS_DR14Q.    #
#          -Prints to the terminal user-friendly feedback about progress on    #
#           number of objects matched. Runtimes not tracked as this will be a  #
#           run-once script for the newest spAll catalogue.                    #
#                                                                              #
# Note: REQUIRES PYTHON 3.5 or later.                                          #
################################################################################

################################################################################
#                                                                              #
#  FUNCTION: GALAXY CONFIRMATION                                               #
#  PURPOSE: This is the function that chooses a flag for a GALAXY that has a   #
#           Z >= 1. This is step 3 in the pseudocode above.                    #
#                                                                              #
#  ACCEPTS: fil  - File name of spZall file to check.                          #
#           p_in - The 'PLATE' value of the object currently being checked.    #
#           m_in - The 'MJD' value of the object currently being checked.      #
#           f_in - The 'FIBERID' value of the object currently being checked.  #
#  RETURNS: out_Gflag - The classification flag for the object, saved to the   #
#             mask_arr that holds all classifications for the file.            #
#  INPUTS:  None                                                               #
#  OUTPUTS: None                                                               #
#  CALLED BY: obj_class() in same file.                                        #
#                                                                              #
################################################################################
def gal3_confirm(fil,p_in,m_in,f_in):
    #OBJ_CLASS() for loop iterates through FIBERIDs for the objects in a file.
    #This function, then, only needs to check the current file for the matching
    #FIBERID. PLATE and MJD are not currently used, but can be in case of a
    #larger or combined spZall file that holds multiple plate-mjd combos.
    wg = np.where(fil['FIBERID'] == f_in)[0]

    #Check the 'CLASS' fields for the next top 4 fits and count them.
    t_gal,t_galC = np.unique(fil['CLASS'][wg[1:5]],return_counts=True)
    tGw = np.where(t_gal == 'GALAXY')[0]

    #The next 2 IF/ELSE statements are error-trapping in case one of the
    #counts in t_galC is 0. If it is, then it doesn't show up in t_gal, so
    #the np.where turns up empty (not 0). These manually force the value to
    #zero in such a case.
    if len(tGw) > 0:
        tG_check = t_galC[tGw][0]
    else:
        tG_check = 0

    #These now assign the classification flag based on the pseudocode in the
    #header above. If the number of galaxies is >= 1, mark it galaxy.
    #Otherwise leve it blank for future visual inspection.
    if tG_check >= 1:
        out_Gflag = 'GALAXY'
    else:
        out_Gflag = ''

    #Return the classification flag to obj_class()
    return out_Gflag

################################################################################
#                                                                              #
#  FUNCTION: QSO CONFIRMATION                                                  #
#  PURPOSE: This is the function that chooses a flag for QSOs from             #
#             step 4 in the pseudocode above.                                  #
#                                                                              #
#  ACCEPTS: fil   - File name of spZall file to check.                         #
#           p_in  - The 'PLATE' value of the object currently being checked.   #
#           m_in  - The 'MJD' value of the object currently being checked.     #
#           f_in  - The 'FIBERID' value of the object currently being checked. #
#           zw_in - The 'ZWARNING' value of the object currently being checked.#
#  RETURNS: out_Qflag - The classification flag for the object, saved to the   #
#             mask_arr that holds all classifications for the catalogue        #
#  INPUTS:  None                                                               #
#  OUTPUTS: None                                                               #
#  CALLED BY: obj_class() in same file.                                        #
#                                                                              #
################################################################################
def qso_confirm(fil,p_in,m_in,f_in,zw_in):
    #OBJ_CLASS() for loop iterates through FIBERIDs for the objects in a file.
    #This function, then, only needs to check the current file for the matching
    #FIBERID. PLATE and MJD are not currently used, but can be in case of a
    #larger or combined spZall file that holds multiple plate-mjd combos.
    wq = np.where(fil['FIBERID'] == f_in)[0]

    #Check the 'CLASS' field of the next top 4 fits. Find 'STAR', 'GALAXY',
    #and 'QSO'
    t_QSO,t_QSOC = np.unique(fil['CLASS'][wq[1:5]],return_counts=True)
    tQSw = np.where(t_QSO == 'STAR')[0]
    tQQw = np.where(t_QSO == 'QSO')[0]

    #The next 2 IF/ELSE statements are error-trapping in case one of the
    #counts in t_QSOC is 0. If it is, then it doesn't show up in the t_QSO array
    #so the np.where turns up empty (not 0). These manually force the value to
    #zero in such a case.
    if len(tQSw) > 0:
        tQS_check = t_QSOC[tQSw][0]
    else:
        tQS_check = 0
    if len(tQQw) > 0:
        tQQ_check = t_QSOC[tQQw][0]
    else:
        tQQ_check = 0

    #These now assign the classification flag based on the pseudocode in the
    #header above. The number of STARS is checked FIRST, then the number of stars
    #is compared along with ZWARNING and the number of other QSOs.
    #Failing these checks a QFLAG classification is assigned (meaning VI later).
    if tQS_check >= 2:
        out_Qflag = 'STAR'
    elif ((tQS_check < 2) & (zw_in == 0)):
        out_Qflag = 'QSO'
    else:
        out_Qflag = 'QFLAG'

    #Send the classification flag back to obj_class() loop that called it.
    return out_Qflag

################################################################################
#                                                                              #
#  FUNCTION: OBJECT CLASSIFICATION MAIN                                        #
#  PURPOSE: This is the function that actually iterates through all of the     #
#           objects in the cutdown catalog file input by the user. It will call#
#           gal3_confirm() and qso_confirm() in different for loops.           #
#                                                                              #
#  ACCEPTS: ifile - The name of the input cutdown catalog file, from cat_cut() #
#  RETURNS: None                                                               #
#  INPUTS:  None                                                               #
#  OUTPUTS: --A new FITS catalog. The spAll file is recreated with a new       #
#             column up front called 'AUTOCLASS_DR14Q' which holds the         #
#             classification flags from mask_arr.                              #
#           --User-friendly text output about current progress.                #
#  OUTPUTS TO: FITS file: Folder where the script was called from.             #
#              TEXT : Screen output.                                           #
#  CALLED BY: None                                                             #
#                                                                              #
################################################################################
def obj_class(ifile):
    #Open the cutdown catalog file and find the total number of objects.
    data = fits.open(ifile)[1].data
    dnum = len(data)

    #Counts the number of objects by type at the beginning for feedback and
    #benchmark purposes. Outputs this to the screen.
    nStar = len(np.where(data['CLASS'] == 'STAR')[0])
    nGal = len(np.where(data['CLASS'] == 'GALAXY')[0])
    nQSO = len(np.where(data['CLASS'] == 'QSO')[0])
    print('\n')
    print('----------------------------------------')
    print('Number of Starting Objects  : '+str(dnum))
    print('Number of Stars before      : '+str(nStar))
    print('Number of Galaxies before   : '+str(nGal))
    print('Number of QSOs before       : '+str(nQSO))
    print('----------------------------------------')
    print('\n')

    #This is a placeholder array to hold local flags based on the algorithm. It
    #will become a new column in the output FITS file called 'AUTOCLASS_DR14Q'.
    mask_arr = np.chararray(dnum,itemsize=6,unicode=True)

    #First Classification | STEP 1 in pseudocode.
    wStar = np.where(data['CLASS'] == 'STAR')[0]
    mask_arr[wStar] = 'STAR'

    #Second Classification | STEP 2 in pseudocode.
    wGal = np.where((data['CLASS'] == 'GALAXY')&(data['Z'] < 1))[0]
    mask_arr[wGal] = 'GALAXY'

    #spZall Classifications
    #Classifications 3 and 4 require the spZall files to be checked. To minimize
    #the number of times a file is loaded into memory, a list of necessary files
    #for the GALAXY classification (3) is made and downloaded in turn. While
    #these files are present, the QSO objects are also checked so that the files
    #won't need to be loaded again later.
    #GALAXY-CLASS objects
    wG2 = np.where((data['CLASS'] == 'GALAXY')&(data['Z'] >= 1))[0]
    count = len(wG2)

    #QSO-CLASS objects
    wQSO = np.where(data['CLASS'] == 'QSO')[0]
    countQ = len(wQSO)

    #Objects are identified by a unique combination of plate number, mjd, and
    #fiberid. This generates a list of single values for each object, so unique
    #combinations can be found. This allows the use of np.where() calls instead
    #of nested for loops later.
    pm = (data['PLATE']**2) + (data['MJD'] - 1)

    #This finds the unique combinations of plate and mjd which define the spZall
    #files. In each file each fiberid will have a number of pipeline fits.
    #These are the unique plate/mjd combos for GALAXIES.
    uniq,uniqindices = np.unique(pm[wG2], return_index=True)
    ugplates = data['PLATE'][wG2[uniqindices]]
    umjd = data['MJD'][wG2[uniqindices]]

    #This does the same as above, but for QSOs. Here this is used not to find
    #the files for loading, but to identify which objects will appear in the
    #files that are loaded for the GALAXY-3 classification.
    uQ,uQind = np.unique(pm[wQSO],return_index=True)
    uqplates = data['PLATE'][wQSO[uQind]]
    uqmjd = data['MJD'][wQSO[uQind]]

    #Makes a list of file names. This is passed to the for loop later.
    file_list = np.chararray(len(ugplates),itemsize=23,unicode=True)
    for j in range(len(ugplates)):
        file_list[j] = 'spZall-{}-{}.fits'.format(ugplates[j],umjd[j])

    #A counter used by the for loop.
    file_count = len(file_list)

    #Initializing the plate, mjd, fiber, run2d variables.
    plate = 0
    mjd = 0
    fiber = 0
    run2d = 0

    #Where are the spZall files stored
    uloc0 = os.environ.get('BOSS_SPECTRO_REDUX')

    #This loop will go through the files identified in file_list (from GAL-3
    #plate/mjd combos). It will iterate through all of the spZall files. While
    #doing so it outputs a single line about its progress that is overwritten.
    for i in range(file_count):
        #Define the PLATE and MJD we are working with on this iteration.
        plate = ugplates[i]
        mjd = umjd[i]
        run2d = data['RUN2D'][wG2[uniqindices[i]]]

        #Get the file name for the current spZall file.
        z_file = '{}/{}/{}/{}/{}'.format(uloc0,run2d,plate,run2d,file_list[i])

        #dall is the spZall file currently being checked.
        dall = fits.open(z_file)[1].data

        #Each object, identified by plate/mjd/fiber has 134 separate fits in the
        #spZall file. These will find only the fibers appearing in the spZall
        #file we are currently using.
        objs = np.where((data['PLATE'][wG2] == plate)&(data['MJD'][wG2] == mjd))[0]
        fibers,fiberind = np.unique(data['FIBERID'][wG2[objs]],return_index=True)

        #Each spZall file holds 1000 fibers. This iterates through the FIBERIDs
        #identified in the previous step, then passes the information to
        #the gal3_confirm() function above.
        for h in range(len(fibers)):
            fiber_t = fibers[h]
            mask_arr[wG2[objs[fiberind[h]]]] = gal3_confirm(dall,plate,mjd,fiber_t)

        #This does the same as above but for any QSOs that might be in the same
        #files.
        objsQ = np.where((data['PLATE'][wQSO] == plate)&(data['MJD'][wQSO] == mjd))[0]
        fibQ,fibQind = np.unique(data['FIBERID'][wQSO[objsQ]],return_index=True)

        #Again, this iterates through the FIBERIDs. ZWARNING is found first,
        #then everything is passed to the qso_confirm() function above.
        for k in range(len(fibQ)):
            fiber_q = fibQ[k]
            z_warn = data['ZWARNING'][wQSO[objsQ[fibQind[k]]]]
            mask_arr[wQSO[objsQ[fibQind[k]]]] = qso_confirm(dall,plate,mjd,fiber_q,z_warn)

        #This outputs the single line that is overwritten as files are completed
        st = '\rFile: {} | spZall Complete: {}/{}'.format(file_list[i],i+1,file_count)
        sys.stdout.write(st)
        sys.stdout.flush()

    #Leftover QSOs
    #This is for any classification 4 QSOs that weren't covered by the GAL-3
    #list of spZall files. It is the same as above, but only for QSOs.
    wlq = np.where((data['CLASS'] == 'QSO')&(mask_arr == ''))[0]
    ulq,ulqindex = np.unique(pm[wlq],return_index=True)
    ulplates = data['PLATE'][wlq[ulqindex]]
    ulmjd = data['MJD'][wlq[ulqindex]]

    #Creates a new list of spZall files, for the leftover QSOs.
    file_lq = np.chararray(len(ulplates),itemsize=23,unicode=True)
    for j in range(len(ulplates)):
        file_lq[j] = 'spZall-{}-{}.fits'.format(ulplates[j],ulmjd[j])

    file_lq_count = len(file_lq)

    #Reinitialize all of the object address information. Don't want old data
    #messing us up.
    plate = 0
    mjd = 0
    fiber = 0
    run2d = 0
    print('\n')
    print('----------------------------------------')
    print('LEFTOVER QSOs')
    print('----------------------------------------')

    for i in range(file_lq_count):
        plate = ulplates[i]
        mjd = ulmjd[i]
        run2d = data['RUN2D'][wlq[ulqindex[i]]]

        z_file = '{}/{}/{}/{}/{}'.format(uloc0,run2d,plate,run2d,file_lq[i])

        dall = fits.open(z_file)[1].data

        objsQ = np.where((data['PLATE'][wlq] == plate)&(data['MJD'][wlq] == mjd))[0]
        fibQ,fibQind = np.unique(data['FIBERID'][wlq[objsQ]],return_index=True)

        for k in range(len(fibQ)):
            fiber_q = fibQ[k]
            z_warn = data['ZWARNING'][wlq[objsQ[fibQind[k]]]]
            mask_arr[wlq[objsQ[fibQind[k]]]] = qso_confirm(dall,plate,mjd,fiber_q,z_warn)

        st1 = '\rFile: {} | spZall Complete: {}/{}'.format(file_lq[i],i+1,file_lq_count)
        sys.stdout.write(st1)
        sys.stdout.flush()

    #TESTING FEEDBACK
    #These lines are only for feedback before and after for quick reference on how
    #many objects need VI. We're hoping for no more than 7% neededing followup.
    wflagged = np.where((mask_arr == 'QFLAG') | (mask_arr == ''))[0]
    nflagged = len(wflagged)
    wt2 = np.where(mask_arr == 'STAR')[0]
    wt3 = np.where(mask_arr == 'GALAXY')[0]
    wt4 = np.where(mask_arr == 'QSO')[0]
    prct_vi = (float(nflagged) / float(dnum)) * 100

    print('\n')
    print('----------------------------------------')
    print('Number of Stars after      : '+str(len(wt2)))
    print('Number of Galaxies after   : '+str(len(wt3)))
    print('Number of QSOs after       : '+str(len(wt4)))
    print('----------------------------------------')
    print('Number to be Visually inspected : '+str(nflagged))
    print('Percent to be Visually inspected: {0:.2f}%'.format(prct_vi))
    print('\n')

    #After classification flags are completed, find any objects still left Unclassified
    #or that have a QFLAG (which is an unclassified QSO). Mark all objects for
    #visual inspection followup.
    mask_arr[wflagged] = 'VI'
    #We also want to visually inspect QSOs with pipeline Z > 3.5.
    wzu = np.where((mask_arr=='QSO')&(data['Z']>3.5))[0]
    mask_arr[wzu] = 'VI'

    min_mjd = np.amin(data['MJD'])
    max_mjd = np.amax(data['MJD'])
    mjd_range = '{0:05d}-{1:05d}'.format(min_mjd,max_mjd)
    spvname = 'spAll-v5_13_0'
    #Write out a fits file with this new classification column. Also modifies
    #the primary HDU header (extension 0) to keep information on what the Input
    #spAll file was, what DR this is for, and what program ran the damn thing.
    prim_hdrc = fits.Header()
    prim_hdrc['ALG_FILE']=('spAll_class_full.py','Algorithm file used')
    prim_hdrc['DATA_REL']=(spvname,'spAll file used')
    prim_hdrc['MJD_RNGE']=(mjd_range,'MJD Range of objects')
    prim_hduc = fits.PrimaryHDU(header=prim_hdrc)

    #Tack this new classification column onto the columns from the cutdown catalog.
    mask_col = fits.ColDefs([fits.Column(name='AUTOCLASS_DR14Q',format='6A',array=mask_arr)])
    data_cols = data.columns
    data_hduc = fits.BinTableHDU.from_columns(mask_col + data_cols)

    ofile_dtag = time.strftime('%Y%m%d')
    data_ofc = fits.HDUList([prim_hduc,data_hduc])
    out_file_namec = '../data/DR16Q_autoclass_{}'.format(ofile_dtag)
    classified_name = ct.fet(data_ofc,out_file_namec)

    return classified_name

#############################CATALOG COMBINE####################################

def cat_combine(inrec_name):

    #These are just the names of the columns that appear in the spAll files.
    spcolkeep = np.array(['RA','DEC','AUTOCLASS_DR14Q','THING_ID','PLATE','MJD','FIBERID',
                        'Z','ZWARNING','BOSS_TARGET1','EBOSS_TARGET0','EBOSS_TARGET1',
                        'EBOSS_TARGET2','ANCILLARY_TARGET1','ANCILLARY_TARGET2',
                        'OBJID','RUN','RERUN','CAMCOL','SKYVERSION','LAMBDA_EFF',
                        'FIELD','PSFFLUX','PSFFLUX_IVAR','PSFMAG','PSFMAGERR',
                        'EXTINCTION','SN_MEDIAN_ALL','ZOFFSET','XFOCAL','YFOCAL','CHUNK',
                        'TILE','PLATESN2','ID'])
    vdcopycols = np.array(['CLASS_PERSON','Z_CONF','Z_VI'])

    #Load the three files I need.
    virec_name = '../data/pre16_vi_database.fits'
    inrec = fits.open(inrec_name)[1].data #My classified FITS record
    virec = fits.open(virec_name)[1].data #The VI database dump from Patrick.

    vnum = 'v5_13_0'
    #Load the full current spAll file too.
    loc = os.environ.get('BOSS_SPECTRO_REDUX')
    full_file = '{0}/{1}/spAll-{2}.fits'.format(loc,vnum,vnum)
    sprec = fits.open(full_file)[1].data

    #Make the hash arrays.
    inhash = ct.mk_hash(inrec)
    vihash = ct.mk_hash(virec)
    sphash = ct.mk_hash(sprec)

    #Find the common objects from my FITS file and the vi database.
    inargs,viargs = ct.rec_match_srt(inhash,vihash)

    #Find the objects in the vi database that aren't in my FITS file.
    vi_mask = np.zeros(len(virec),dtype='i2')
    vi_mask[viargs] = 1
    wv0 = np.where(vi_mask==0)[0]
    vinhash = ct.mk_hash(virec[wv0])
    #Find those objects that aren't in my FITS file that ARE in the spAll file.
    vinargs,spargs = ct.rec_match_srt(vinhash,sphash)

    #The total number of records the output will hold. This is without duplicate
    #reduction.
    num_in = len(inrec)
    num_vin = len(vinhash)
    num_tot = len(inrec) + len(vinargs)

    #Create the structured array that will hold my record first.
    inarr = ecm.supercat_maker(num_in)
    vinarr = ecm.supercat_maker(num_vin)

    #Copy everything from my record to the structured array.
    for cname in spcolkeep:
        if cname == 'Z':
            cname1 = 'Z_PIPE'
        elif cname == 'RUN':
            cname1 = 'RUN_NUMBER'
        elif cname == 'RERUN':
            cname1 = 'RERUN_NUMBER'
        elif cname == 'CAMCOL':
            cname1 = 'CAMCOL_NUMBER'
        elif cname == 'FIELD':
            cname1 = 'FIELD_NUMBER'
        elif cname == 'ID':
            cname1 = 'ID_NUMBER'
        else:
            cname1 = cname
        inarr[cname1] = inrec[cname]

    #Copy the VI information from the VI database to the struct.
    for cname in vdcopycols:
        inarr[cname][inargs] = virec[cname][viargs]

    #Load everything from the spAll file for the VI objects that didn't match my record.
    for cname in spcolkeep:
        if cname == 'Z':
            cname1 = 'Z_PIPE'
        elif cname == 'RUN':
            cname1 = 'RUN_NUMBER'
        elif cname == 'RERUN':
            cname1 = 'RERUN_NUMBER'
        elif cname == 'CAMCOL':
            cname1 = 'CAMCOL_NUMBER'
        elif cname == 'FIELD':
            cname1 = 'FIELD_NUMBER'
        elif cname == 'ID':
            cname1 = 'ID_NUMBER'
        elif cname == 'AUTOCLASS_DR14Q':
            vinarr['AUTOCLASS_DR14Q'][vinargs] = 'UNK'
            continue
        else:
            cname1 = cname
        vinarr[cname1][vinargs] = sprec[cname][spargs]

    #Copy the VI information from the VI database.
    for cname in vdcopycols:
        vinarr[cname][vinargs] = virec[cname][wv0[vinargs]]

    #Copy both placeholder structs into one struct.
    dsarr = ecm.supercat_maker(num_tot)
    colnames = np.array(dsarr.dtype.names)
    for cname in colnames:
        dsarr[cname][0:num_in] = inarr[cname][:]
        dsarr[cname][num_in:num_tot] = vinarr[cname][vinargs]

    #Now we need to generate the SDSS names that go with each quasar
    name_tarr = sdn.s_name(dsarr['RA'],dsarr['DEC'])
    dsarr['SDSS_NAME'] = name_tarr

    #And sort the catalog by SDSS name, which is standard
    name_args = np.argsort(dsarr['SDSS_NAME'])
    dsarr = dsarr[name_args]

    #We might want all of the records, before any removals. This does that.
    prim_hdr = fits.open(inrec_name)[0].header
    prim_hdu = fits.PrimaryHDU(header=prim_hdr)

    #We need to make the output file name for the catalog that includes
    #the records with bad ZWARNING flags.
    outfile_dtag = time.strftime('%Y%m%d')
    outfile_name = '../data/DR16Q_vdb_zw_{}'.format(outfile_dtag)

    #Make the HDU1 catalog and write it out.
    data_hdu = fits.BinTableHDU.from_columns(dsarr)
    data_of = fits.HDUList([prim_hdu,data_hdu])
    vdb_out_name = ct.fet(data_of,outfile_name)

    return vdb_out_name

##################REMOVE BAD ZWARNING FLAGS#####################################

def zw_remove(ifile):
    #We also want to remove all records with bad ZWARNING flags.
    zwarn_arr = ct.file_load(ifile)

    zwarn_mask = np.zeros(len(zwarn_arr),dtype='i2')
    zwarn_mask[:] = 1
    unplug_val = flagval('ZWARNING',['UNPLUGGED','SKY','LITTLE_COVERAGE','NODATA',
                         'BAD_TARGET'])
    wunplugged = np.where(zwarn_arr['ZWARNING']&unplug_val)[0]
    zwarn_mask[wunplugged] = 0
    wplugged = np.where(zwarn_mask == 1)[0]

    #Again, we need to grab the HDU0 header information
    prim_hdr = fits.open(ifile)[0].header
    prim_hdu = fits.PrimaryHDU(header=prim_hdr)

    #We also need to make another output file name
    outfile_dtag = time.strftime('%Y%m%d')
    outfile_name = '../data/DR16Q_vdb_rem_{}'.format(outfile_dtag)

    #And make the output data HDU1 and HDU list
    data_out = np.array(zwarn_arr[wplugged])
    data_hdu = fits.BinTableHDU.from_columns(data_out)
    data_hdulist = fits.HDUList([prim_hdu,data_hdu])

    #Write out that file. This is the last write out from this program
    vdb_zw_name = ct.fet(data_hdulist,outfile_name)

    #Return the name of the file that contains VI database data.
    #This is not the superset yet.
    #Note to me: I have gotten this far in git-prep.
    return vdb_zw_name

#############COMMAND LINE FUNCTION CALLS########################################

#This requires the utilities scripts:
#    progressBar.py, sdss_name_v2.py, cat_tools.py, and empty_cat_maker.py
#This requires the following data file in ../data/
#    pre16_vi_database.fits

#Run from command line with:
#  python spAll_class_full.py
#  and wait.
if __name__=='__main__':
    supername = cat_cut() #This will run the superset cutdown.
    clsnamed = obj_class(supername) #This will run the Classification Algorithm.
    vdb_combo_name = cat_combine(clsnamed) #This combines the superset with the VI db.
    final_outname = zw_remove(vdb_combo_name) #This removes records with bad ZWARNINGs.

reason: import numpy,from astropy


text:
#!/usr/bin/env python2
# -*-coding:utf-8 -*-

import numpy as np
import os
import sys
import tensorflow as tf
import copy
import yaml
import cv2
import tarfile
import six.moves.urllib as urllib
from tensorflow.core.framework import graph_pb2
from collections import defaultdict
from io import StringIO
import matplotlib
from matplotlib import pyplot as plt

from object_detection.utils import label_map_util
from object_detection.utils import visualization_utils as vis_util
from stuff.helper import FPS2, WebcamVideoStream

# for ros
import rospy
from std_msgs.msg import String
from std_msgs.msg import Header
from sensor_msgs.msg import Image
from cv_bridge import CvBridge, CvBridgeError
from copy import deepcopy
from detection_2d.msg import ObjectInfoWithROI
from detection_2d.msg import ObjectInfoArray

CONFIG_PATH = os.path.dirname(sys.path[0]) + '/config'+'/config.yml'
print("config path : %r" % (CONFIG_PATH))

with open(CONFIG_PATH, 'r') as ymlfile:
    cfg = yaml.load(ymlfile)
video_input         = cfg['video_input']
visualize           = cfg['visualize']
vis_text            = cfg['vis_text']
max_frames          = cfg['max_frames']
width               = cfg['width']
height              = cfg['height']
fps_interval        = cfg['fps_interval']
allow_memory_growth = cfg['allow_memory_growth']
det_interval        = cfg['det_interval']
det_th              = cfg['det_th']
model_name          = cfg['model_name']
model_path          = cfg['model_path']
label_path          = cfg['label_path']
num_classes         = cfg['num_classes']
split_model         = cfg['split_model']
log_device          = cfg['log_device']

MODEL_PATH = os.path.dirname(sys.path[0]) + "/" + model_path
print("model path : %r" % (MODEL_PATH))

LABEL_PATH = os.path.dirname(sys.path[0]) + "/scripts" + "/" + label_path
print("label path : %r" % (LABEL_PATH))

# Download Model form TF's Model Zoo
def download_model():
    model_file = model_name + '.tar.gz'
    download_base = 'http://download.tensorflow.org/models/object_detection/'
    if not os.path.isfile(MODEL_PATH):
        print('Model not found. Downloading it now.')
        opener = urllib.request.URLopener()
        opener.retrieve(download_base + model_file, model_file)
        tar_file = tarfile.open(model_file)
        for file in tar_file.getmembers():
          file_name = os.path.basename(file.name)
          if 'frozen_inference_graph.pb' in file_name:
            tar_file.extract(file, os.getcwd() + '/models/')
        os.remove(os.getcwd() + '/' + model_file)
    else:
        print('Model found. Proceed.')

def _node_name(n):
  if n.startswith("^"):
    return n[1:]
  else:
    return n.split(":")[0]

def load_frozenmodel():
    print('Loading frozen model into memory')
    if not split_model:
        detection_graph = tf.Graph()
        with detection_graph.as_default():
          od_graph_def = tf.GraphDef()
          with tf.gfile.GFile(model_path, 'rb') as fid:
            serialized_graph = fid.read()
            od_graph_def.ParseFromString(serialized_graph)
            tf.import_graph_def(od_graph_def, name='')
        return detection_graph, None, None

    else:
        # load a frozen Model and split it into GPU and CPU graphs
        # need to add gpu_options
        input_graph = tf.Graph()
        config = tf.ConfigProto()
        config.gpu_options.allow_growth = allow_memory_growth
        with tf.Session(graph=input_graph, config=config):
            score = tf.placeholder(tf.float32, shape=(None, 1917, 90), name="Postprocessor/convert_scores")
            expand = tf.placeholder(tf.float32, shape=(None, 1917, 1, 4), name="Postprocessor/ExpandDims_1")
            for node in input_graph.as_graph_def().node:
                if node.name == "Postprocessor/convert_scores":
                    score_def = node
                if node.name == "Postprocessor/ExpandDims_1":
                    expand_def = node

        detection_graph = tf.Graph()
        with detection_graph.as_default():
          od_graph_def = tf.GraphDef()
          with tf.gfile.GFile(MODEL_PATH, 'rb') as fid:
            serialized_graph = fid.read()
            od_graph_def.ParseFromString(serialized_graph)
            dest_nodes = ['Postprocessor/convert_scores','Postprocessor/ExpandDims_1']

            edges = {}
            name_to_node_map = {}
            node_seq = {}
            seq = 0
            for node in od_graph_def.node:
              n = _node_name(node.name)
              name_to_node_map[n] = node
              edges[n] = [_node_name(x) for x in node.input]
              node_seq[n] = seq
              seq += 1

            for d in dest_nodes:
              assert d in name_to_node_map, "%s is not in graph" % d

            nodes_to_keep = set()
            next_to_visit = dest_nodes[:]
            while next_to_visit:
              n = next_to_visit[0]
              del next_to_visit[0]
              if n in nodes_to_keep:
                continue
              nodes_to_keep.add(n)
              next_to_visit += edges[n]

            nodes_to_keep_list = sorted(list(nodes_to_keep), key=lambda n: node_seq[n])

            nodes_to_remove = set()
            for n in node_seq:
              if n in nodes_to_keep_list: continue
              nodes_to_remove.add(n)
            nodes_to_remove_list = sorted(list(nodes_to_remove), key=lambda n: node_seq[n])

            keep = graph_pb2.GraphDef()
            for n in nodes_to_keep_list:
              keep.node.extend([copy.deepcopy(name_to_node_map[n])])

            remove = graph_pb2.GraphDef()
            remove.node.extend([score_def])
            remove.node.extend([expand_def])
            for n in nodes_to_remove_list:
              remove.node.extend([copy.deepcopy(name_to_node_map[n])])

            with tf.device('/gpu:0'):
              tf.import_graph_def(keep, name='')
            with tf.device('/cpu:0'):
              tf.import_graph_def(remove, name='')

        return detection_graph, score, expand

def load_labelmap():
    print('Loading label map')
    label_map = label_map_util.load_labelmap(LABEL_PATH)
    categories = label_map_util.convert_label_map_to_categories(label_map, max_num_classes=num_classes, use_display_name=True)
    category_index = label_map_util.create_category_index(categories)
    return category_index

class ObjectDetection(object):
    def __init__(self, detection_graph, category_index, score, expand):
        self.image_sub = rospy.Subscriber("/image", Image, self.imageCallback, queue_size=10)
        self.image_pub = rospy.Publisher("/object_image", Image, queue_size=10)
        self.bbox_pub = rospy.Publisher("/object_info", ObjectInfoArray, queue_size=10)
        self.detection_graph = detection_graph
        self.category_index = category_index
        self.score = score
        self.expand = expand
        self.flag = False

    def imageCallback(self, image_msg):
        try:
            self.cv_image = CvBridge().imgmsg_to_cv2(image_msg, "bgr8")
            self.frame_id = image_msg.header.frame_id
            self.flag = True
        except CvBridgeerror as e:
            print (e)

    def publish_boxes_and_labels(
                     self,
                     image,
                     boxes,
                     classes,
                     scores,
                     category_index,
                     max_boxes_to_draw=20,
                     min_score_thresh=.30):
        bbox_array = ObjectInfoArray()
        if not max_boxes_to_draw:
            max_boxes_to_draw = boxes.shape[0]
        for i in range(min(max_boxes_to_draw, boxes.shape[0])):
            if scores is None or scores[i] > min_score_thresh:
                if classes[i] in category_index.keys():
                    class_name = category_index[classes[i]]['name']
                else:
                    class_name = 'N/A'
                box = tuple(boxes[i].tolist())
                label = classes[i]
                score = scores[i]
                ymin, xmin, ymax, xmax = box
                im_height = image.shape[0]
                im_width  = image.shape[1]
                (left, right, top, dowm) = (xmin * im_width, xmax * im_width,
                        ymin * im_height, ymax * im_height)
                bbox = ObjectInfoWithROI()
                bbox.Class = class_name
                bbox.probability = score
                bbox.xmin = left
                bbox.xmax = right
                bbox.ymin = top
                bbox.ymax = dowm
                bbox_array.object_array.append(bbox)
                print('%s:%.2f left:%d right:%d top:%d dowm:%d' % (class_name, score, left, right, top, dowm))
        bbox_array.header.stamp = rospy.Time.now()
        bbox_array.header.frame_id = self.frame_id
        self.bbox_pub.publish(bbox_array)

    def object_detection(self, sess):
        image_np = self.cv_image
        image_np = cv2.resize(image_np, (height, width))
        image_np_expanded = np.expand_dims(image_np, axis=0)

        # actual Detection
        if not split_model:
            (boxes, scores, classes, num) = sess.run(
                [self.detection_boxes, self.detection_scores, self.detection_classes, self.num_detections],
                feed_dict={self.image_tensor: image_np_expanded})
        else:
            # Split Detection in two sessions.
            (score, expand) = sess.run([self.score_out, self.expand_out], feed_dict={self.image_tensor: image_np_expanded})
            (boxes, scores, classes, num) = sess.run(
                [self.detection_boxes, self.detection_scores, self.detection_classes, self.num_detections],
                feed_dict={self.score_in:score, self.expand_in: expand})

        # publish object info
        self.publish_boxes_and_labels(
            self.cv_image,
            np.squeeze(boxes),
            np.squeeze(classes).astype(np.int32),
            np.squeeze(scores),
            self.category_index)

        # Visualization of the results of a detection.
        if visualize:
            vis_util.visualize_boxes_and_labels_on_image_array(
                image_np,
                np.squeeze(boxes),
                np.squeeze(classes).astype(np.int32),
                np.squeeze(scores),
                self.category_index,
                use_normalized_coordinates=True,
                line_thickness=8)

        image_height = self.cv_image.shape[0]
        image_width  = self.cv_image.shape[1]
        resize_image = cv2.resize(image_np, (image_width, image_height))
        pub_image = CvBridge().cv2_to_imgmsg(resize_image, "bgr8")
        self.image_pub.publish(pub_image)

    def main(self):
        rospy.init_node('detection_node')
        rate = rospy.Rate(30)
        per = rospy.get_param('~per', 0.8)
        dev = rospy.get_param('~dev', "0")

        config = tf.ConfigProto(allow_soft_placement=True, log_device_placement=log_device)
        config.gpu_options.allow_growth=allow_memory_growth

        with self.detection_graph.as_default():
            with tf.Session(graph=self.detection_graph,config=config) as sess:
                # Define Input and Ouput tensors
                self.image_tensor        = self.detection_graph.get_tensor_by_name('image_tensor:0')
                self.detection_boxes     = self.detection_graph.get_tensor_by_name('detection_boxes:0')
                self.detection_scores    = self.detection_graph.get_tensor_by_name('detection_scores:0')
                self.detection_classes   = self.detection_graph.get_tensor_by_name('detection_classes:0')
                self.num_detections      = self.detection_graph.get_tensor_by_name('num_detections:0')
                if split_model:
                    self.score_out       = self.detection_graph.get_tensor_by_name('Postprocessor/convert_scores:0')
                    self.expand_out      = self.detection_graph.get_tensor_by_name('Postprocessor/ExpandDims_1:0')
                    self.score_in        = self.detection_graph.get_tensor_by_name('Postprocessor/convert_scores_1:0')
                    self.expand_in       = self.detection_graph.get_tensor_by_name('Postprocessor/ExpandDims_1_1:0')

                    while not rospy.is_shutdown():
                        if self.flag:
                            self.object_detection(sess)
                            self.flag = False
                        rate.sleep()

def main():
    download_model()
    graph, score, expand = load_frozenmodel()
    category = load_labelmap()
    detection = ObjectDetection(graph, category, score, expand)
    detection.main()

if __name__ == '__main__':
    main()


text:
# coding=utf-8
# Copyright 2019 The TensorFlow Datasets Authors.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""Cifar-10.1 dataset."""

from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import numpy as np
import tensorflow as tf

import tensorflow_datasets.public_api as tfds

_CITATION = """\
@article{recht2018cifar10.1,
  author = {Benjamin Recht and Rebecca Roelofs and Ludwig Schmidt and Vaishaal Shankar},
  title = {Do CIFAR-10 Classifiers Generalize to CIFAR-10?},
  year = {2018},
  note = {\\url{https://arxiv.org/abs/1806.00451}},
}

@article{torralba2008tinyimages,
  author = {Antonio Torralba and Rob Fergus and William T. Freeman},
  journal = {IEEE Transactions on Pattern Analysis and Machine Intelligence},
  title = {80 Million Tiny Images: A Large Data Set for Nonparametric Object and Scene Recognition},
  year = {2008},
  volume = {30},
  number = {11},
  pages = {1958-1970}
}
"""

_DESCRIPTION = """\
The CIFAR-10.1 dataset is a new test set for CIFAR-10. CIFAR-10.1 contains roughly 2,000 new test images
that were sampled after multiple years of research on the original CIFAR-10 dataset. The data collection
for CIFAR-10.1 was designed to minimize distribution shift relative to the original dataset. We describe
the creation of CIFAR-10.1 in the paper "Do CIFAR-10 Classifiers Generalize to CIFAR-10?".
The images in CIFAR-10.1 are a subset of the TinyImages dataset.
There are currently two versions of the CIFAR-10.1 dataset: v4 and v6.
"""

_DL_URL_IMAGES = "https://github.com/modestyachts/CIFAR-10.1/blob/master/datasets/cifar10.1_{}_data.npy?raw=true"
_DL_URL_LABELS = "https://github.com/modestyachts/CIFAR-10.1/blob/master/datasets/cifar10.1_{}_labels.npy?raw=true"

_DATA_OPTIONS = ["v4", "v6"]

class Cifar10_1Config(tfds.core.BuilderConfig):  # pylint: disable=invalid-name
  """BuilderConfig for Cifar-10.1."""

  @tfds.core.disallow_positional_args
  def __init__(self, data=None, **kwargs):
    """Constructs a Cifar10_1Config.

    Args:
      data: `str`, one of `_DATA_OPTIONS`.
      **kwargs: keyword arguments forwarded to super.
    """
    if data not in _DATA_OPTIONS:
      raise ValueError("data must be one of %s" % _DATA_OPTIONS)
    kwargs.setdefault("name", data)
    super(Cifar10_1Config, self).__init__(**kwargs)
    self.data = data

class Cifar10_1(tfds.core.GeneratorBasedBuilder):  # pylint: disable=invalid-name
  """Cifar-10.1 dataset."""

  BUILDER_CONFIGS = [
      Cifar10_1Config(
          description=(
              "It is the first version of our dataset on which we tested any classifier. As mentioned above, this "
              "makes the v4 dataset independent of the classifiers we evaluate. The numbers reported in the main "
              "sections of our paper use this version of the dataset. It was built from the top 25 TinyImages "
              "keywords for each class, which led to a slight class imbalance. The largest difference is that ships "
              "make up only 8% of the test set instead of 10%. v4 contains 2,021 images."
          ),
          version=tfds.core.Version("1.0.0"),
          data="v4",
      ),
      Cifar10_1Config(
          description=(
              "It is derived from a slightly improved keyword allocation that is exactly class balanced. This version "
              "of the dataset corresponds to the results in Appendix D of our paper. v6 contains 2,000 images."
          ),
          version=tfds.core.Version("1.0.0"),
          data="v6",
      )
  ]

  def _info(self):
    return tfds.core.DatasetInfo(
        builder=self,
        description=_DESCRIPTION,
        features=tfds.features.FeaturesDict({
            "image": tfds.features.Image(shape=(32, 32, 3)),
            "label": tfds.features.ClassLabel(num_classes=10),
        }),
        supervised_keys=("image", "label"),
        urls=["https://github.com/modestyachts/CIFAR-10.1"],
        citation=_CITATION,
    )

  def _split_generators(self, dl_manager):
    """Returns SplitGenerators."""

    image_url = _DL_URL_IMAGES.format(self.builder_config.data)
    label_url = _DL_URL_LABELS.format(self.builder_config.name)

    image_path, label_path = dl_manager.download([
        image_url,
        label_url,
    ])

    return [
        tfds.core.SplitGenerator(
            name=tfds.Split.TEST,
            gen_kwargs={
                "image_path": image_path,
                "label_path": label_path,
            }),
    ]

  def _generate_examples(self, image_path, label_path):
    with tf.io.gfile.GFile(image_path, "rb") as f:
      images = np.load(f)
    with tf.io.gfile.GFile(label_path, "rb") as f:
      labels = np.load(f)
    for i, (image, label) in enumerate(zip(images, labels)):
      record = {
          "image": image,
          "label": label,
      }
      yield i, record

max_stars_repo_name: SiegeLordEx/datasets

max_issues_repo_name: SiegeLordEx/datasets

max_forks_repo_name: SiegeLordEx/datasets


text:
# -----------------------------------------------------------------------------
#   @author:
#       Tingwu Wang
#   @brief:
#       Define the dynamic models for the system, which takes input two
#       adjacent states and output the predicted actions.
# -----------------------------------------------------------------------------
import tensorflow as tf
import numpy as np

from .base_dynamics import base_dynamics_network
from mbbl.config import init_path
# from mbbl.util import model_saver
from mbbl.util.common import logger
from mbbl.util.common import tf_networks
# from mbbl.util.common import tf_norm
# from mbbl.util.common import tf_utils

class dynamics_network(base_dynamics_network):
    '''
        @brief:
    '''

    def __init__(self, args, session, name_scope,
                 observation_size, action_size):
        '''
            @input:
                @ob_placeholder:
                    if this placeholder is not given, we will make one in this
                    class.

                @trainable:
                    If it is set to true, then the policy weights will be
                    trained. It is useful when the class is a subnet which
                    is not trainable
        '''
        super(dynamics_network, self).__init__(
            args, session, name_scope, observation_size, action_size
        )
        self._base_dir = init_path.get_abs_base_dir()
        self._debug_it = 0

    def build_network(self):
        # the placeholders
        self._build_ph()

        self._tensor = {}

        # construct the input to the forward network, we normalize the state
        # input, and concatenate with the action
        self._tensor['normalized_start_state'] = (
            self._input_ph['start_state'] -
            self._whitening_operator['state_mean']
        ) / self._whitening_operator['state_std']

        self._tensor['net_input'] = tf.concat(
            [self._tensor['normalized_start_state'],
             self._input_ph['action']], 1
        )

        # the setting of mlp network for every layer (by using a list), provide
        # the activation function, normalization function and initialzation
        network_shape = [self._observation_size + self._action_size] + \
            self.args.dynamics_network_shape + [self._observation_size]
        num_layer = len(network_shape) - 1
        act_type = \
            [self.args.dynamics_activation_type] * (num_layer - 1) + [None]
        norm_type = \
            [self.args.dynamics_normalizer_type] * (num_layer - 1) + [None]
        init_data = []
        for _ in range(num_layer):
            init_data.append(
                {'w_init_method': 'xavier', 'w_init_para': {'uniform': False},
                 'b_init_method': 'xavier', 'b_init_para': {'uniform': False}}
            )

        self._MLP = tf_networks.MLP(
            dims=network_shape, scope='dynamics_mlp', train=True,
            activation_type=act_type, normalizer_type=norm_type,
            init_data=init_data
        )

        self._tensor['net_output'] = self._MLP(self._tensor['net_input'])

        # get the predicted next state
        self._tensor['pred_output'] = self._tensor['net_output'] * \
            self._whitening_operator['diff_state_std'] + \
            self._whitening_operator['diff_state_mean'] + \
            self._input_ph['start_state']

        # fetch all the trainable variables
        self._set_var_list()

    def build_loss(self):
        self._update_operator = {}

        # get the groundtruth output
        self._tensor['normalized_state_diff'] = \
            (self._input_ph['end_state'] - self._input_ph['start_state'] -
             self._whitening_operator['diff_state_mean']) / \
            self._whitening_operator['diff_state_std']
        # get the loss and optimizer
        self._update_operator['pred_error'] = tf.square(
            self._tensor['net_output'] - self._tensor['normalized_state_diff']
        )
        # from util.common.fpdb import fpdb; fpdb().set_trace()
        self._update_operator['loss'] = \
            tf.reduce_mean(self._update_operator['pred_error'])

        self._update_operator['update_op'] = tf.train.AdamOptimizer(
            learning_rate=self.args.dynamics_lr,
            # beta1=0.5, beta2=0.99, epsilon=1e-4
        ).minimize(self._update_operator['loss'])

    def train(self, data_dict, replay_buffer, training_info={}):
        # update the whitening stats of the network
        self._set_whitening_var(data_dict['whitening_stats'])
        self._debug_it += 1

        # get the validation set
        new_data_id = list(range(len(data_dict['start_state'])))
        self._npr.shuffle(new_data_id)
        num_val = min(int(len(new_data_id) * self.args.dynamics_val_percentage),
                      self.args.dynamics_val_max_size)
        val_data = {
            key: data_dict[key][new_data_id][:num_val]
            for key in ['start_state', 'end_state', 'action']
        }

        # get the training set
        replay_train_data = replay_buffer.get_all_data()
        train_data = {
            key: np.concatenate(
                [data_dict[key][new_data_id][num_val:], replay_train_data[key]]
            ) for key in ['start_state', 'end_state', 'action']
        }

        for i_epochs in range(self.args.dynamics_epochs):
            # get the number of batches
            num_batches = len(train_data['action']) // \
                self.args.dynamics_batch_size
            # from util.common.fpdb import fpdb; fpdb().set_trace()
            assert num_batches > 0, logger.error('batch_size > data_set')
            avg_training_loss = []

            for i_batch in range(num_batches):
                # train for each sub batch
                feed_dict = {
                    self._input_ph[key]: train_data[key][
                        i_batch * self.args.dynamics_batch_size:
                        (i_batch + 1) * self.args.dynamics_batch_size
                    ] for key in ['start_state', 'end_state', 'action']
                }
                fetch_dict = {
                    'update_op': self._update_operator['update_op'],
                    'train_loss': self._update_operator['loss']
                }

                training_stat = self._session.run(fetch_dict, feed_dict)
                avg_training_loss.append(training_stat['train_loss'])

            val_loss = self.eval(val_data)

            logger.info('[dynamics]: Val Loss: {}, Train Loss: {}'.format(
                val_loss, np.mean(avg_training_loss))
            )

            '''
            if self._debug_it > 20:
                for i in range(2):
                    test_feed_dict = {
                        self._input_ph[key]: feed_dict[self._input_ph[key]][[i]]
                        for key in ['start_state', 'end_state', 'action']
                    }
                    t_end_state = self._session.run(self._tensor['pred_output'], test_feed_dict)
                    print('train_loss', training_stat)
                    # print('gt', test_feed_dict[self._input_ph['end_state']])
                    # print('pred', t_end_state)
                    diff = t_end_state - test_feed_dict[self._input_ph['end_state']]
                    print(i, 'pred_diff', diff, np.abs(diff).max())
                    print('end-start diff',
                          test_feed_dict[self._input_ph['end_state']] -
                          test_feed_dict[self._input_ph['start_state']])
                    print('pred_end-start diff',
                          t_end_state,
                          test_feed_dict[self._input_ph['start_state']])

                from util.common.fpdb import fpdb; fpdb().set_trace()
            '''
        training_stat['val_loss'] = val_loss
        training_stat['avg_train_loss'] = np.mean(avg_training_loss)
        return training_stat

    def get_weights(self):
        return self._get_network_weights()

    def set_weights(self, weight_dict):
        return self._set_network_weights(weight_dict)

    def eval(self, data_dict):
        feed_dict = {self._input_ph[key]: data_dict[key]
                     for key in ['start_state', 'end_state', 'action']}
        return self._session.run(self._update_operator['loss'], feed_dict)

    def pred(self, data_dict):
        feed_dict = {self._input_ph[key]: data_dict[key]
                     for key in ['start_state', 'action']}
        return self._session.run(self._tensor['pred_output'], feed_dict), -1, -1


text:
import numpy as np
import h5py

def get_data():
    A01T = h5py.File("./project_datasets/A01T_slice.mat", 'r')
    X = np.copy(A01T["image"])
    y = np.copy(A01T["type"])
    y = y[0, 0:X.shape[0]:1]
    y = np.asarray(y, dtype=np.int32)
    X = X[:, :22, :]

    for i in range(len(y)):
        if y[i] == 769:
            y[i] = 0
        elif y[i] == 770:
            y[i] = 1
        elif y[i] == 771:
            y[i] = 2
        elif y[i] == 772:
            y[i] = 3
        else:
            y[i] = 0
    X = X.reshape(-1, 22*1000)
    y = y[~np.isnan(X).any(axis=1)]
    X = X[~np.isnan(X).any(axis=1)]
    X = X.reshape(-1, 22, 1000)

    # Define training and testing data
    t_size = np.array(range(X.shape[0]))
    np.random.shuffle(t_size)
    test_mask = t_size[:50]
    train_mask = t_size[50:]
    X_test = X[test_mask]
    y_test = y[test_mask]
    X_train = X[train_mask]
    y_train = y[train_mask]

    return X_test, y_test, X_train, y_train

def get_all_data():

    A01T = h5py.File("./project_datasets/A01T_slice.mat", 'r')
    X = np.copy(A01T["image"])
    y = np.copy(A01T["type"])
    y = y[0, 0:X.shape[0]:1]
    y = np.asarray(y, dtype=np.int32)
    X = X[:, :22, :]

    for i in range(len(y)):
        if y[i] == 769:
            y[i] = 0
        elif y[i] == 770:
            y[i] = 1
        elif y[i] == 771:
            y[i] = 2
        elif y[i] == 772:
            y[i] = 3
        else:
            y[i] = 0
    X = X.reshape(-1, 22*1000)
    y = y[~np.isnan(X).any(axis=1)]
    X = X[~np.isnan(X).any(axis=1)]
    X = X.reshape(-1, 22, 1000)

    # Define training and testing data
    t_size = np.array(range(X.shape[0]))
    np.random.shuffle(t_size)
    test_mask = t_size[:50]
    train_mask = t_size[50:]
    X_test = X[test_mask]
    y_test = y[test_mask]
    X_train = X[train_mask]
    y_train = y[train_mask]

    for j in range(2, 10):
        string_name = "./project_datasets/A0" + str(j) + "T_slice.mat"
        A01T = h5py.File(string_name, 'r')
        X = np.copy(A01T["image"])
        y = np.copy(A01T["type"])
        y = y[0, 0:X.shape[0]:1]
        y = np.asarray(y, dtype=np.int32)
        X = X[:, :22, :]

        for i in range(len(y)):
            if y[i] == 769:
                y[i] = 0
            elif y[i] == 770:
                y[i] = 1
            elif y[i] == 771:
                y[i] = 2
            elif y[i] == 772:
                y[i] = 3
            else:
                y[i] = 0
        X = X.reshape(-1, 22 * 1000)
        y = y[~np.isnan(X).any(axis=1)]
        X = X[~np.isnan(X).any(axis=1)]
        X = X.reshape(-1, 22, 1000)

        # Define training and testing data
        t_size = np.array(range(X.shape[0]))
        np.random.shuffle(t_size)
        test_mask = t_size[:50]
        train_mask = t_size[50:]
        X_test = np.concatenate((X_test, X[test_mask]), axis=0)
        y_test = np.concatenate((y_test, y[test_mask]), axis=0)
        X_train = np.concatenate((X_train, X[train_mask]), axis=0)
        y_train = np.concatenate((y_train, y[train_mask]), axis=0)

    return X_test, y_test, X_train, y_train

def crop_data(X_test, y_test, X_train, y_train):
    X_test_out = []
    y_test_out = []
    X_train_out = []
    y_train_out = []
    X_test_out = X_test[:, :, :500]
    y_test_out = y_test
    X_train_out = X_train[:, :, :500]
    y_train_out = y_train

    for i in range(1, 501, 25):
        X_test_out = np.concatenate((X_test_out, X_test[:, :, i:i+500]), axis=0)
        # y_test_out = np.concatenate((y_test_out, y_test))
        X_train_out = np.concatenate((X_train_out, X_train[:, :, i:i+500]), axis=0)
        y_train_out = np.concatenate((y_train_out, y_train))
    print(X_test.shape)
    return X_test_out, y_test_out, X_train_out, y_train_out

def noise_data(X_test, y_test, X_train, y_train):
    X_test_out = X_test
    y_test_out = y_test
    X_train_out = X_train
    y_train_out = y_train

    for _ in range(4):
        X_train_out = np.concatenate((X_train_out, np.random.normal(size=X_train.shape) + X_train), axis=0)
        y_train_out = np.concatenate((y_train_out, y_train))

    return X_test_out, y_test_out, X_train_out, y_train_out

max_stars_repo_name: charliezjw/Neural-Signal-Decoder

max_issues_repo_name: charliezjw/Neural-Signal-Decoder

max_forks_repo_name: charliezjw/Neural-Signal-Decoder


text:
### TODO ###
### source activate alignment
import os
import sys
import numpy
import scipy
import math
import logging
import time
import rlcompleter
import readline

import datetime as dt
import epics as PyEpics
import matplotlib.pyplot as plt

from scipy.special import erf
from scipy.optimize import curve_fit
from pprint import pprint
from mpl_toolkits.mplot3d import axes3d

### PIN TRACKING PACKAGE
# sys.path.append('/home/beams/S1IDUSER/auto_sample_alignment/src')
from align import Alignment

def tracking(params, algorithm='pin', repeat=1):

    #### ENABLE TAB COMPLETION
    readline.parse_and_bind('tab: complete')

    #################################################
    pixel_size = 1.172; ## um / pixel

    slit_center_x=1000;
    slit_center_y=600;

    #################################################
    ### MOTOR SETTINGS
    #################################################
    mtr_samXE = PyEpics.Motor('1ide1:m34');
    mtr_samYE = PyEpics.Motor('1ide1:m35');
    mtr_samZE = PyEpics.Motor('1ide1:m36');
    mtr_samTh = PyEpics.Motor('1ide1:m86');
    mtr_samChi = PyEpics.Motor('1ide1:m87');
    mtr_samOme = PyEpics.Motor('1ide:m9');
    mtr_aeroXE = PyEpics.Motor('1ide1:m101');
    mtr_aeroZE = PyEpics.Motor('1ide1:m102');
    #################################################

    pname='/home/beams/S1IDUSER/mnt/s1c/mli_nov19/tomo';

    for i in range(repeat):
        #################################################
        ### ALIGN ROTATION AXIS TO THE HS
        #################################################
        ### MOVE STAGE TO 0
        mtr_samOme.move(0, wait=True)

        ### TAKE AN IMAGE AT OME 0
        PyEpics.caput('1idPG1:cam1:Acquire', 1)
        time.sleep(3)

        fname=PyEpics.caget('1idPG1:TIFF1:FileName_RBV', 'str') + "_%06d"%(PyEpics.caget('1idPG1:TIFF1:FileNumber_RBV')-1) + '.tif'
        pfname=os.path.join(pname, fname)

        print(pfname)

        align0 = Alignment(pfname)
        x0, y0 = align0.compute_center(algorithm, params)
        print(f"pin x,y position when omega is 0 : (x = {x0}, y = {y0}")

        if (x0 == -1) or (y0 == -1):
            print("Alignment failed at zero degress")
            break

        ### MOVE STAGE TO 180
        mtr_samOme.move(180, wait=True)

        ### TAKE AN IMAGE AT OME 180
        PyEpics.caput('1idPG1:cam1:Acquire', 1)
        time.sleep(3)

        fname=PyEpics.caget('1idPG1:TIFF1:FileName_RBV', 'str') + "_%06d"%(PyEpics.caget('1idPG1:TIFF1:FileNumber_RBV')-1) + '.tif'
        pfname=os.path.join(pname, fname)

        align180 = Alignment(pfname)
        x180, y180 = align180.compute_center(algorithm, params)
        if (x180 == -1) or (y180 == -1):
            print("Alignment failed at 180 degrees")
            break

        print(f"pin x,y position when omega is 180 : (X = {x180}, y = {180})")

        ### COMPUTE MOTIONS TO MAKE
        mid_x = (x180 + x0)/2;
        half_delta_x = (x180 - x0)/2;

        print(mid_x)
        print(half_delta_x)

        ### NEED TO CHECK / FIGURE OUT THE SIGNS ON THESE
        aeroXE_motion = ((mid_x - slit_center_x)*pixel_size)/1000;
        samXE_motion = -(half_delta_x*pixel_size)/1000;

        print('motions to execute')
        print(aeroXE_motion)
        print(samXE_motion)

        mtr_aeroXE.move(aeroXE_motion, relative=True, wait=True)
        mtr_samXE.move(samXE_motion, relative=True, wait=True)

reason: import numpy,import scipy,from scipy