<repo_name>andrewasheridan/hera_sandbox<file_sep>/aas/Estimating Delays/nn/data_creation/Data_Creator_C.py
"""Summary
"""
# Data_Creator_C

from data_manipulation import *
from Data_Creator import Data_Creator

import numpy as np

class Data_Creator_C(Data_Creator):
    """Data_Creator_C C = Classification
    
    Creates one dimensional inputs and their labels for training
    Inputs are simulated angle data generated with true data as their base
        - see ????.ipynb for discussion of input data
    Labels are one-hot vectors (optionally blurred vectors)
        - see ????.ipynb for discussion of labels and blurring
    
    Each input is generated:
        - a pair of redundant baselines is randomly chosen
        - the ratio of visibilties is constructed (a complex flatness, see discussion)
        - a separate target cable delay is applied to each row of the flatness
            - targets in the range (-abs_min_max_delay, abs_min_max_delay)
        - the angle is computed of the flatness with the applied delay
    
    Each label is generated:
        - the smooth range of targets is converted to a series of steps of the desired precision
            - see Classifier_Class_Explanation.ipynb for discussion
        - each unique step value is assigned a one-hot vector (from min to max)
        - each target is assigned its label target
        - if blur != 0 blur is applied
            - Imagine each one-hot vector as a delta function with area = 1
            - blur converts the delta function into a gaussian with area = 1
            - see ????.ipynb for discussion
    
    Args:
        num_flatnesses (int): number of flatnesses used to generate data.
            Number of data samples = 60 * num_flatnesses
        bl_data (dict): data source. Output of get_seps_data()
        bl_dict (dict): Dictionary of seps with bls as keys. An output of get_or_gen_test_train_red_bls_dicts()
        gains (dict): Gains for this data. An output of load_relevant_data()
        abs_min_max_delay (float, optional): targets are generated in range(-abs_min_max_delay, abs_min_max_delay)
        precision (float, optional): size of the class steps must be from (0.005,0.001,0.0005,0.00025,0.0001)
        evaluation (bool, optional): experimental
        single_dataset (bool, optional): experimental
        singleset_path (None, optional): experimental
        blur (float, optional): blur value. Convert a one-hot vector delta function to a gaussian.
            Spreads the 1.0 prob of a class over a range of nearby classes.
            blur = 0.5 spreads over about 7 classes with a peak prob of 0.5 (varies)
            blue = 0.1 spreads over about 40 classes with peak prob ~ 0.06 (varies)
            blur = 0.0, no blur
    
    Attributes:
        blur (float): Description
        evaluation (bool): experimental
        precision (float): Description
        single_dataset (bool): experimental
        singleset_path (bool): experimental
    """
    __doc__ += Data_Creator.__doc__

    def __init__(self,
                 num_flatnesses,
                 bl_data,
                 bl_dict ,
                 gains,
                 abs_min_max_delay = 0.040,
                 precision = 0.00025,
                 evaluation = False,
                 single_dataset = False,
                 singleset_path = None,
                 blur = 0.): # try 1, 0.5, 0.10
        """Summary
        
        Args:
            num_flatnesses (int): number of flatnesses used to generate data.
                Number of data samples = 60 * num_flatnesses
            bl_data (dict): data source. Output of get_seps_data()
            bl_dict (dict): Dictionary of seps with bls as keys. An output of get_or_gen_test_train_red_bls_dicts()
            gains (dict): Gains for this data. An output of load_relevant_data()
            abs_min_max_delay (float, optional): targets are generated in range(-abs_min_max_delay, abs_min_max_delay)
            precision (float, optional): size of the class steps must be from (0.005,0.001,0.0005,0.00025,0.0001)
            evaluation (bool, optional): experimental
            single_dataset (bool, optional): experimental
            singleset_path (None, optional): experimental
            blur (float, optional): blur value. Convert a one-hot vector delta function to a gaussian.
                Spreads the 1.0 prob of a class over a range of nearby classes.
                blur = 0.5 spreads over about 7 classes with a peak prob of 0.5 (varies)
                blue = 0.1 spreads over about 40 classes with peak prob ~ 0.06 (varies)
                blur = 0.0, no blur
        """
        Data_Creator.__init__(self,
                              num_flatnesses = num_flatnesses,
                              bl_data = bl_data,
                              bl_dict = bl_dict,
                              gains = gains,
                              abs_min_max_delay = abs_min_max_delay)

        self.precision = precision
        self.evaluation = evaluation
        self.single_dataset = single_dataset
        self.singleset_path = singleset_path
        self.blur = blur
        
    def _load_singleset(self):
        """_load_singleset

        load a set of specific baseline-pairs
        """
        def loadnpz(filename):
            """loadnpz
            
            Args:
                filename (str): path
            
            Returns:
                numpy array: loaded npz
            """
            a = np.load(filename)
            d = dict(zip(("data1{}".format(k) for k in a), (a[k] for k in a)))
            return d['data1arr_0']

        return [[(s[0][0],s[0][1]), (s[1][0],s[1][1])] for s in loadnpz(self.singleset_path)]
    
    def _blur(self, labels):
        """_blur

        Convert a one-hot vector delta function to a gaussian.
        Spreads the 1.0 prob of a class over a range of nearby classes.
        blur = 0.5 spreads over about 7 classes with a peak prob of 0.5
        blur = 0.1 spreads over about 40 classes with peak prob ~ 0.06
        
        Args:
            labels (list of lists of ints): class labels
        
        Returns:
            list of lists of floats: class labels, blurred
        """
        for i, label in enumerate(labels):
            true_val = np.argmax(label)
            mean = true_val; std = self.blur; variance = std**2
            x = np.arange(len(label))
            f = np.exp(-np.square(x-mean)/2*variance)
            labels[i] = f/np.sum(f)
        return labels
             
    def _gen_data(self):            
        """_gen_data

        Generates artiicial flatnesses by combining two redundant visbilites and their gains.
        Applies different known delay to each row of the floatnesses.
        Scales data for network.

        Converts numerical delay valyue to classification label
        """
        # scaling tools
        # the NN likes data in the range (0,1)
        angle_tx  = lambda x: (np.asarray(x) + np.pi) / (2. * np.pi)
        angle_itx = lambda x: np.asarray(x) * 2. * np.pi - np.pi

        delay_tx  = lambda x: (np.array(x) + self._tau) / (2. * self._tau)
        delay_itx = lambda x: np.array(x) * 2. * self._tau - self._tau
        
        targets = np.random.uniform(low = -self._tau, high = self._tau, size = (self._num * 60, 1))
        applied_delay = np.exp(-2j * np.pi * (targets * self._nu + np.random.uniform()))



        assert type(self._bl_data) != None, "Provide visibility data"
        assert type(self._bl_dict) != None, "Provide dict of baselines"
        assert type(self._gains)   != None, "Provide antenna gains"

        if self._bl_data_c == None:
            self._bl_data_c = {key : self._bl_data[key].conjugate() for key in self._bl_data.keys()}

        if self._gains_c == None:
            self._gains_c = {key : self._gains[key].conjugate() for key in self._gains.keys()}


        def _flatness(seps):
            """_flatness

            Create a flatness from a given pair of seperations, their data & their gains.
            
            Args:
                seps (list of tuples of ints): two redundant separations 
            
            Returns:
                numpy array of complex floats: visibility ratio flatnesss
            """

            a, b = seps[0][0], seps[0][1]
            c, d = seps[1][0], seps[1][1]


            return self._bl_data[seps[0]]   * self._gains_c[(a,'x')] * self._gains[(b,'x')] * \
                   self._bl_data_c[seps[1]] * self._gains[(c,'x')]   * self._gains_c[(d,'x')]

        inputs = []
        seps = []

        singleset_seps = None if self.singleset_path == None else self._load_singleset()

        for i in range(self._num):

            unique_baseline = random.sample(self._bl_dict.keys(), 1)[0]
            if self.singleset_path == None:
                two_seps = [random.sample(self._bl_dict[unique_baseline], 2)][0]
            elif self.singleset_path != None:
                two_seps = singleset_seps[i]
            inputs.append(_flatness(two_seps))
            for t in range(60):
                seps.append(two_seps)

        inputs = np.angle(np.array(inputs).reshape(-1,1024) * applied_delay)
        
        permutation_index = np.random.permutation(np.arange(self._num * 60))
        
        
#         TODO: Quadruple check that rounded_targets is doing what you think it is doing

        # rounded_targets is supposed to take the true targets and round their value
        # to the value closest to the desired precision. and the absolute.
        # pretty sure i have it working for the values below
        # classes is supposed to have all the possible unique rounded values
        
        #0.0001 precision - 401 classes
        if self.precision == 0.0001:
            rounded_targets = np.asarray([np.round(abs(np.round(d * 100,2)/100), 5) for d in targets[permutation_index]]).reshape(-1)
            classes = np.arange(0,0.04 + self.precision, self.precision)
        
        #0.00025 precision - 161 classes
        if self.precision == 0.00025:
            rounded_targets = np.asarray([np.round(abs(np.round(d * 40,2)/40), 5) for d in targets[permutation_index]]).reshape(-1)
            classes = np.arange(0,0.04 + self.precision, self.precision)

        # 0.0005 precision 81 classes
        if self.precision == 0.0005:
            rounded_targets = np.asarray([np.round(abs(np.round(d * 20,2)/20), 5) for d in targets[permutation_index]]).reshape(-1)
            classes = np.arange(0,0.04 + self.precision, self.precision)

        # 0.001 precision - 41 classes
        if self.precision == 0.001:
            rounded_targets = np.asarray([np.round(abs(np.round(d * 10,2)/10), 5) for d in targets[permutation_index]]).reshape(-1)
            classes = np.arange(0,0.04 + self.precision, self.precision)

        # for 0.005 precision there should be 9 different classes
        if self.precision == 0.005:
            rounded_targets = np.asarray([np.round(abs(np.round(d * 2,2)/2), 5) for d in targets[permutation_index]]).reshape(-1)
            classes = np.arange(0,0.04 + self.precision, self.precision)

        eye = np.eye(len(classes), dtype = int)
        classes_labels = {}
        for i, key in enumerate(classes):
            classes_labels[np.round(key,5)] = eye[i].tolist()
            
            
        labels = [classes_labels[x] for x in rounded_targets]
        labels = labels if self.blur == 0 else self._blur(labels)
        
        if self.evaluation == True:
            self._epoch_batch.append((angle_tx(inputs[permutation_index]), labels, targets[permutation_index], np.asarray(seps)[permutation_index]))
        if self.single_dataset == True:
            self._epoch_batch.append((angle_tx(inputs[permutation_index]), labels, np.asarray(seps)[permutation_index]))
        else:
            self._epoch_batch.append((angle_tx(inputs[permutation_index]), labels))
        

<file_sep>/aas/Estimating Delays/nn/network_trainers/CNN_C_Trainer.py
"""Summary
"""
# CNN_C_Trainer

import sys, os
sys.path.insert(1, os.path.join(sys.path[0], '../modules'))
from NN_Trainer import NN_Trainer

import tensorflow as tf
from tensorflow.python.client import timeline

import numpy as np
from sklearn.metrics import confusion_matrix
import matplotlib.pyplot as plt
import io
import itertools

np.seterr(divide='ignore', invalid='ignore') # for cm div/zero (handled)

class CNN_C_Trainer(NN_Trainer):
    """CNN_C_Trainer
        - for training classifcation networks
        - see ????.ipynb for discussion
        
    Args:
        - sample_keep_prob - Probability of keeping a value from an input row
        - conv_keep_prob - Probability of keeping an output value from a convolution
        - num_classes - int - Number of classification classes
            - make sure this matches the value from the network
            - why isnt this number just loaded from the network...
        - single_dataset - bool - train on a single fixed dataset

    """
    __doc__ += NN_Trainer.__doc__

    def __init__(self,
                 network,
                 Data_Creator, # Class
                 num_classes,
                 num_epochs = 100,
                 batch_size = 32,
                 log_dir = '../logs/',
                 model_save_interval = 25,
                 pretrained_model_path = None,
                 metric_names = ['costs', 'accuracies'],
                 sample_keep_prob = 0.80,
                 conv_keep_prob = 0.9,
                 verbose = True,
                 single_dataset = False):
        """__init__"""
        NN_Trainer.__init__(self,
                            network = network,
                            Data_Creator = Data_Creator,
                            num_epochs = num_epochs,
                            batch_size = batch_size,
                            log_dir = log_dir,
                            model_save_interval = model_save_interval,
                            pretrained_model_path = pretrained_model_path,
                            metric_names = metric_names,
                            verbose = verbose)
        

        self.sample_keep_prob = sample_keep_prob
        self.conv_keep_prob = conv_keep_prob
        self.num_classes = num_classes
        self.single_dataset = single_dataset

    def add_data(self,train_info, test_info, gains, num_flatnesses = 10, abs_min_max_delay = 0.040, precision = 0.00025, blur = 0):
        """Adds data to the Trainer.

            Args

                train_info - (tuple : (visibility data, baselines dict)) - The training data
                test_info - (tuple : (visibility data, baselines dict)) - The testing data
                gains - (dict) - The gains for the visibilites
                num_flatnesses - (int) - Number of flatnesses for each epoch
                    - Number of samples is num_flatnesses * 60
                abs_min_max_delay - (float) - The absolute value of the min or max delay for this dataset.

        """
        self._abs_min_max_delay = abs_min_max_delay

        self._train_batcher = self._Data_Creator(num_flatnesses,
                                                 train_info[0],
                                                 train_info[1],
                                                 gains,
                                                 abs_min_max_delay,
                                                 precision,
                                                 single_dataset = self.single_dataset,
                                                 blur = blur)
        self._train_batcher.gen_data()


        self._test_batcher = self._Data_Creator(num_flatnesses,
                                                test_info[0],
                                                test_info[1],
                                                gains,
                                                abs_min_max_delay,
                                                precision,
                                                single_dataset = self.single_dataset,
                                                blur = blur)
        self._test_batcher.gen_data()

    def train(self):
        
        self.save_params()



        costs = []
        accuracies = []

        tf.reset_default_graph()
        
        self._network.create_graph()
        self._network.save_params()
        
        saver = tf.train.Saver()

        with tf.Session() as session:
            
            if self.pretrained_model_path == None:
                session.run(tf.global_variables_initializer())
                
            else:
                saver.restore(session, self.pretrained_model_path)

            archive_loc = self.log_dir + self._network.name
            training_writer = tf.summary.FileWriter(archive_loc + '/training', session.graph)
            testing_writer = tf.summary.FileWriter(archive_loc + '/testing', session.graph)
            self.model_save_location = archive_loc + '/trained_model.ckpt'
            
            
            self._msg = '\rtraining';self._vprint(self._msg)

            try:
                for epoch in range(self.num_epochs):
                    
                    if self.single_dataset == True and epoch == 0:
                        training_inputs, training_labels, training_bl_dict_singleset = self._train_batcher.get_data(); self._train_batcher.gen_data()
                        testing_inputs, testing_labels, testing_bl_dict_singleset = self._test_batcher.get_data(); self._test_batcher.gen_data()
                        np.savez(archive_loc + '/params/train_bl_dict_singleset', training_bl_dict_singleset)   
                        np.savez(archive_loc + '/params/test_bl_dict_singleset', testing_bl_dict_singleset)
                        
                    if self.single_dataset == False:
                        training_inputs, training_labels = self._train_batcher.get_data(); self._train_batcher.gen_data()
                        testing_inputs, testing_labels = self._test_batcher.get_data(); self._test_batcher.gen_data() 

                    training_labels = np.asarray(training_labels)
                    testing_labels = np.asarray(testing_labels)
                    
                    # if the division here has a remainde some values are just truncated
                    batch_size = self.batch_size
                    num_entries = training_inputs.shape[0]

                    for j in range(int(num_entries/batch_size)):
                        self._msg = '\repoch'
                        self._msg += '- {:5.0f}/{:5.0f}'.format(epoch + 1,self.num_epochs)
                        self._msg += ' - batch: {:4.0f}/{:4.0f}'.format(j + 1, int(num_entries/batch_size))


                        training_inputs_batch = training_inputs[j*batch_size:(j + 1)*batch_size].reshape(-1,1,1024,1)
                        training_labels_batch = training_labels[j*batch_size:(j + 1)*batch_size].reshape(-1,self.num_classes)

                        feed_dict = {self._network.X: training_inputs_batch,
                                     self._network.labels: training_labels_batch,
                                     self._network.sample_keep_prob : self.sample_keep_prob,
                                     self._network.conv_keep_prob : self.conv_keep_prob,
                                     self._network.is_training : True}

                        batch_cost, batch_acc, _ = session.run([self._network.cost, self._network.accuracy, self._network.optimizer], feed_dict = feed_dict)
                        if epoch != 0:
                            self._msg += " batch cost: {:0.4f}".format(batch_cost)
                            self._msg += " batch accs: {:2.2f}".format(batch_acc)
                        self._vprint(self._msg); 
                        
                            
                    train_feed_dict = {self._network.X: training_inputs.reshape(-1,1,1024,1),
                                       self._network.labels: training_labels.reshape(-1,self.num_classes),
                                       self._network.sample_keep_prob : 1.,
                                       self._network.conv_keep_prob : 1.,
                                       self._network.is_training : False}

                    train_predicts =  session.run([self._network.predictions], train_feed_dict)

                    train_feed_dict[self._network.image_buf] = self.plt_confusion_matrix(training_labels.reshape(-1,self.num_classes), train_predicts)

                    training_cost, training_acc, training_summary = session.run([self._network.cost,
                                                                                 self._network.accuracy,
                                                                                 self._network.summary],
                                                                                 feed_dict = train_feed_dict) 

                    training_writer.add_summary(training_summary, epoch)
                    training_writer.flush()  
                
                    
                    test_feed_dict = {self._network.X: testing_inputs.reshape(-1,1,1024,1),
                                      self._network.labels: testing_labels.reshape(-1,self.num_classes),
                                      self._network.sample_keep_prob : 1.,
                                      self._network.conv_keep_prob : 1.,
                                      self._network.is_training : False}

                    test_predicts =  session.run([self._network.predictions], test_feed_dict)

                    test_feed_dict[self._network.image_buf] = self.plt_confusion_matrix(testing_labels.reshape(-1,self.num_classes), test_predicts) 

                    testing_cost, testing_acc, testing_summary = session.run([self._network.cost,
                                                                              self._network.accuracy,
                                                                              self._network.summary],
                                                                              feed_dict = test_feed_dict)

                    testing_writer.add_summary(testing_summary, epoch)
                    testing_writer.flush()
                    
                    costs.append((training_cost, testing_cost))
                    accuracies.append((training_acc, testing_acc))
                    
                    if (epoch + 1) % self.model_save_interval == 0:
                        saver.save(session, self.model_save_location, epoch + 1)
                
                self.msg = ''
            except KeyboardInterrupt:
                self._msg = ' TRAINING INTERUPPTED' # this never prints I dont know why
                pass

            self._msg += '\rtraining ended'; self._vprint(self._msg)
            
            training_writer.close()
            testing_writer.close()


        session.close()
        self._msg += ' - session closed'; self._vprint(self._msg)
        self._msg = ''

        self._metrics = [costs, accuracies]
        self.save_metrics()
        
    def plt_confusion_matrix(self, labels, pred, normalize=False, title='Confusion matrix'):
        """
        Given one-hot encoded labels and preds, displays a confusion matrix.
        Arguments:
            `labels`:
                The ground truth one-hot encoded labels.
            `pred`:
                The one-hot encoded labels predicted by a model.
            `normalize`:
                If True, divides every column of the confusion matrix
                by its sum. This is helpful when, for instance, there are 1000
                'A' labels and 5 'B' labels. Normalizing this set would
                make the color coding more meaningful and informative.
        """
        labels = [label.argmax() for label in np.asarray(labels).reshape(-1,self.num_classes)] # bc
        pred = [label.argmax() for label in np.asarray(pred).reshape(-1,self.num_classes)] #bc

        #classes = ['pos','neg']#np.arange(len(set(labels)))
        if self.num_classes == 9:
            precision = 0.005
        if self.num_classes == 41:
            precision = 0.001
        if self.num_classes == 81:
            precision = 0.0005
        if self.num_classes == 161:
            precision = 0.00025
        if self.num_classes == 401:
            precision = 0.0001
        if self.num_classes == 2:
            classes = ['pos','neg']
            cm = confusion_matrix(labels, pred)
        else:
            classes = np.round(np.arange(0,0.04 + precision, precision), 5)
            cm = confusion_matrix(labels, pred, np.arange(len(classes)))

        #if normalize:
        cm = cm.astype('float')*100 / cm.sum(axis=1)[:, np.newaxis]
        cm = np.nan_to_num(cm, copy=True)
        cm = cm.astype('int')

        fig, ax = plt.subplots(figsize = (5,5), dpi = 320)
        #plt.figure(figsize=(15,10))
        im = ax.imshow(cm, interpolation='nearest', aspect='auto', cmap='Oranges', vmin = 0, vmax = 100)
        ax.set_title(title)
        cbar = fig.colorbar(im)
        if len(classes) <= 161:
            tick_marks = np.arange(len(classes))

            fs = 4 if self.num_classes <= 41 else 3
            ax.set_xticks(tick_marks)
            ax.set_xticklabels(classes, fontsize=fs, rotation=-90,  ha='center')
            ax.set_yticks(tick_marks)
            ax.set_yticklabels(classes, fontsize=fs)
            ax.xaxis.set_tick_params(width=0.5)
            ax.yaxis.set_tick_params(width=0.5)
        else:
            ax.set_xticks([])
            ax.set_yticks([])
        
        ax.set_ylabel('True label')
        ax.set_xlabel('Predicted label')
        if self.num_classes <= 41:
            for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
                #s = '{:2.0}'.format(cm[i, j]) if cm[i,j] >= 1 else '.'
                ax.text(j, i, format(cm[i, j], 'd') if cm[i,j]!=0 else '.',
                        horizontalalignment="center", fontsize=5, verticalalignment='center', color= "black")

        fig.set_tight_layout(True)
        # plt.show()
        buf = io.BytesIO()
        fig.savefig(buf, format='png', dpi = 320)
        plt.close(fig)
        buf.seek(0)

        return buf.getvalue()

<file_sep>/aas/Estimating Delays/nn/networks/CNN_DS_BN_C.py
"""Summary
"""
# CNN_DS_BN_C

import sys, os

from RestoreableComponent import RestoreableComponent

import tensorflow as tf
import numpy as np

class CNN_DS_BN_C(RestoreableComponent):
    """CNN_DS_BN_C() - Child of RestoreableComponent
    
    CNN: Convolutional Neural Network.
    DS: DownSampling. Each layer ends with a downsampling convolution
    BN: All non-linearalities have batch-normalization applied.
    C: Classification, this network classifies samples as having one of many labels.
    
    Network Structure:
        Incoming sample has dropout sample_keep_prob - set with Trainer
    
        Each layer has 4 convolutions, each with 4**(i+1) filters (i = zero based layer index).
            - Each convolution:
                - is a 1D convolution
                - feeds into the next
                - has a filter of size (1, fw)
                - has biases and a LeakyReLU activation
                - has dropout with conv_keep_prob - set with Trainer
                - has batch normalization
            - Filter widths (fw) of the four convolutions are [3,5,7,9]
            TODO: Optional filter widths? (different per layer?)
            - Fourth convolution will do 50% downsample (horizontal stride = 2)
    
        Final convolution feeds into fully connected layer as logits
        Cost function softmax_cross_entropy_with_logits_v2
        Optimizer is Adam with adam_initial_learning_rate
        Predictions are softmax probabilites
    
    Usage:
        - create object and set args (or load params)
    
        Training :
            - pass object into appropriate network Trainer object
            - trainer will run create_graph()
    
        Structure Check:
            - run create_graph()
            - call network._layers to see layer output dimensions
            TODO: Make structure check better
                - should directly access _layers
                - should add save graph for tensorboard ?
    
        Other:
            - run create_graph()
            - start tensorflow session
    
    Methods:
        - create_graph() - Contructs the network graph
        
    Attributes:
        accuracy (tensorflow obj): Running accuracy of predictions
        adam_initial_learning_rate (tensorflow obj): Adam optimizer initial learning rate
        conv_keep_prob (tensorflow obj): Keep prob rate for convolutions
        correct_prediction (tensorflow obj): compares predicted label to true, for predictions
        cost (tensorflow obj): cost function
        dtype (tensorflow obj): Type used for all computations
        image_buf (tensorflow obj): buffer for image summaries for tensorboard
        is_training (tensorflow obj): flag for batch normalization
        labels (tensorflow obj): tensor of sample labels
        num_classes (int): Number of different classes
        num_downsamples (int): Number of downsample convolutions
        optimizer (tensorflow obj): optimization function
        pred_cls (tensorflow obj): predicted class index
        predictions (tensorflow obj): probability predictions or each class
        sample_keep_prob (tensorflow obj): keep rate for incoming sample
        summary (tensorflow obj): summary operation for tensorboard
        true_cls (tensorflow obj): true class index
        X (tensorflow obj): incoming sample
    
    """
    __doc__ += RestoreableComponent.__doc__

    def __init__(self,
                 name,
                 num_downsamples,
                 num_classes,
                 log_dir = '../logs/',
                 dtype = tf.float32,
                 adam_initial_learning_rate = 0.0001,
                 verbose = True):
        """__init__
        
        Args:
            name (str): Name of this network
            num_downsamples (int): Number of downsample convolutions
            num_classes (int): Number of different classes
            log_dir (str, optional): Log directory
            dtype (tensorboard datatype, optional): datatype for all operations
            adam_initial_learning_rate (float, optional): Adam optimizer initial learning rate
            verbose (bool, optional): be verbose
        """
        RestoreableComponent.__init__(self, name=name, log_dir=log_dir, verbose=verbose)
                
        self.num_downsamples = num_downsamples
        self.dtype = dtype
        self.adam_initial_learning_rate = adam_initial_learning_rate
        self._num_freq_channels = 1024
        self._layers = []
        self.num_classes = num_classes #  classifier...

    def create_graph(self):
        """create_graph

        Create the network graph for use in a tensorflow session
        """

        
        self._msg = '\rcreating network graph '; self._vprint(self._msg)

        
        tf.reset_default_graph()
        self.is_training = tf.placeholder(dtype = tf.bool, shape = [], name = 'is_training')
        
        with tf.variable_scope('keep_probs'):
            self._msg += '.'; self._vprint(self._msg)

            self.sample_keep_prob = tf.placeholder(self.dtype, name = 'sample_keep_prob')
            self.conv_keep_prob = tf.placeholder(self.dtype, name = 'conv_keep_prob')
        
        with tf.variable_scope('sample'):
            self._msg += '.'; self._vprint(self._msg)
            
            self.X = tf.placeholder(self.dtype, shape = [None, 1, self._num_freq_channels, 1], name = 'X')
            self.X = tf.nn.dropout(self.X, self.sample_keep_prob)

        for i in range(self.num_downsamples):
            self._msg += '.'; self._vprint(self._msg)
            
            # foc - filter out channels - number of filters
            foc = 4**(i+1) # num filters grows with each downsample
            
            with tf.variable_scope('downsample_{}'.format(i)):

                layer = self.X if i == 0 else self._layers[-1]
                fitlter_widths = [3,5,7,9]

                for fw in fitlter_widths:
                    with tf.variable_scope('conv_1x{}'.format(fw)):
                        
                        layer = tf.layers.conv2d(inputs = layer,
                                                 filters = foc, 
                                                 kernel_size = (1,fw),
                                                 strides = (1,1) if fw != fitlter_widths[-1] else (1,2), # downscale last conv
                                                 padding = 'SAME',
                                                 activation = tf.nn.leaky_relu,
                                                 use_bias = True,
                                                 bias_initializer = tf.contrib.layers.xavier_initializer())
                        
                        layer = tf.nn.dropout(layer, self.conv_keep_prob)
                        layer = tf.layers.batch_normalization(layer, training = self.is_training)
                        
            self._layers.append(layer)              
                
        self._msg += ' '
        with tf.variable_scope('labels'):
            self._msg += '.'; self._vprint(self._msg)
            
            self.labels = tf.placeholder(dtype = self.dtype, shape = [None, self.num_classes], name = 'labels')
            self.true_cls = tf.argmax(self.labels, axis = 1)
            
        with tf.variable_scope('logits'):
            self._msg += '.'; self._vprint(self._msg)
            
            self._logits = tf.contrib.layers.fully_connected(tf.layers.flatten(self._layers[-1]), self.num_classes, activation_fn = None)

        with tf.variable_scope('predictions'):
            self._msg += '.'; self._vprint(self._msg)
            
            self.predictions = tf.nn.softmax(self._logits)
            self.pred_cls = tf.argmax(self.predictions, axis = 1)
            self.correct_prediction = tf.equal(self.pred_cls, self.true_cls)
            self.accuracy = tf.reduce_mean(tf.cast(self.correct_prediction, self.dtype))
        
        
        with tf.variable_scope('costs'):
            self._msg += '.'; self._vprint(self._msg)

            self._cross_entropy = tf.nn.softmax_cross_entropy_with_logits_v2(labels = self.labels, logits = self._logits)
            self.cost = tf.reduce_mean(self._cross_entropy)
        
        # batch normalization requires the update_ops variable and control dependency
        update_ops = tf.get_collection(tf.GraphKeys.UPDATE_OPS)
        with tf.control_dependencies(update_ops):
            with tf.variable_scope('train'):
                self._msg += '.'; self._vprint(self._msg)

                if self.dtype == tf.float16:
                    epsilon=1e-4 # optimizer outputs NaN otherwise :(
                else:
                    epsilon=1e-8

                self.optimizer = tf.train.AdamOptimizer(self.adam_initial_learning_rate, epsilon=epsilon).minimize(self.cost)
            
        with tf.variable_scope('logging'):
            self._msg += '.'; self._vprint(self._msg)

            with tf.variable_scope('image'):
                
                self.image_buf = tf.placeholder(tf.string, shape=[])
                epoch_image = tf.expand_dims(tf.image.decode_png(self.image_buf, channels=4), 0)
            
            tf.summary.scalar(name = 'cost', tensor = self.cost)
            tf.summary.scalar(name = 'accuracy', tensor = self.accuracy)
            tf.summary.image('confusion_matrix', epoch_image)

            self.summary = tf.summary.merge_all()
            
        num_trainable_params = np.sum([np.prod(v.get_shape().as_list()) for v in tf.trainable_variables()])
        self._msg = '\rnetwork Ready - {} trainable parameters'.format(num_trainable_params); self._vprint(self._msg)
<file_sep>/aas/Estimating Delays/estdel/estdel/__init__.py
from estdel import VratioDelaySign, VratioDelayMagnitude, VratioDelay, DelaySolver
<file_sep>/aas/Estimating Delays/nn/network_trainers/__init__.py
from NN_Trainer import NN_Trainer 
from CNN_BC_Trainer import CNN_BC_Trainer
from CNN_C_Trainer import CNN_C_Trainer
from CNN_DS_BN_R_Trainer import CNN_DS_BN_R_Trainer
from CNN_R_Trainer import CNN_R_Trainer
<file_sep>/aas/Estimating Delays/estdel/estdel/tests/test_estdel.py
import unittest
from estdel.estdel import _DelayPredict
import numpy as np

class test_DelayPredict(unittest.TestCase):

	def test_init(self):

		data = np.exp(-2j*np.pi*np.arange(1024)*0.01).reshape(-1,1024)

		delayPredict = _DelayPredict(data)

		self.assertEqual(delayPredict._n_freqs, 1024)
		np.testing.assert_array_equal(delayPredict._data, data)

		data_2 = np.exp(-2j*np.pi*np.arange(1024 + 1)*0.01).reshape(-1,1024 + 1)

		self.assertRaises(AssertionError, _DelayPredict, data_2)




if __name__ == '__main__':
	unittest.main()
<file_sep>/aas/Estimating Delays/nn/networks/CNN_DS_BN_R.py
"""CNN_DS_BN_R
"""

import sys, os

from RestoreableComponent import RestoreableComponent

import tensorflow as tf
import numpy as np

class CNN_DS_BN_R(RestoreableComponent):

    """CNN_DS_BN_R
    
    Attributes:
        accuracy_threshold (float): the value for a target to be considered 'good'
        adam_initial_learning_rate (float): Adam optimizer initial learning rate
        cost (str): name of cost function. 'MSE', 'MQE', 'MISG', 'PWT_weighted_MSE', 'PWT_weighted_MISG'
            - use 'MSE', others experimental
        downsample_keep_prob (tensorflow obj): keep prob for downsampling layers
        dtype (tensorflow obj): type used for all computations
        gaussian_shift_scalar (float): value to shift gaussian for 'MISG' cost
        image_buf (tensorflow obj): buffer for image summaries for tensorboard
        is_training (tensorflow obj): flag for batch normalization
        MISG (tensorflow obj): cost function, experimental. Mean Inverse Shifted Gaussian 
        MQE (tensorflow obj): cost function, experimental. Mean Quad Error
        MSE (tensorflow obj): cost function. Mean Squared Error
        num_downsamples (int): Number of downsampling layers
        optimizer (tensorflow obj): optimization function
        predictions (tensorflow obj): predicted outputs
        PWT (tensorflow obj): Percent of predictions within threshold
        sample_keep_prob (tensorflow obj): keep prob for input samples
        summary (tensorflow obj): summary operation for tensorboard
        targets (tensorflow obj): true values for optimizer
        X (tensorflow obj): input samples
    """
    __doc__ += RestoreableComponent.__doc__
    
    def __init__(self,
                 name,
                 num_downsamples,
                 cost = 'MSE',
                 log_dir = '../logs/',
                 dtype = tf.float32,
                 adam_initial_learning_rate = 0.0001,
                 accuracy_threshold = 0.00625,
                 gaussian_shift_scalar = 1e-5,
                 verbose = True):
        """Summary
        
        Args:
            name (str): name of network
            num_downsamples (int): number of downsampling layers
            cost (str, optional): name of cost function. 'MSE', 'MQE', 'MISG', 'PWT_weighted_MSE', 'PWT_weighted_MISG'
                    - use 'MSE', others experimental
            log_dir (str, optional): Directory to store network model and params
            dtype (tensorflow obj, optional): type used for all computations
            adam_initial_learning_rate (float, optional): Adam optimizer initial learning rate
            accuracy_threshold (float, optional): the value for a target to be considered 'good'
            gaussian_shift_scalar (float, optional): value to shift gaussian for 'MISG' cost
            verbose (bool, optional): Be verbose 
        """
        RestoreableComponent.__init__(self, name=name, log_dir=log_dir, verbose=verbose)
                
        self.num_downsamples = num_downsamples
        
        self.dtype = dtype
        self.adam_initial_learning_rate = adam_initial_learning_rate
        self.accuracy_threshold = accuracy_threshold
        self.gaussian_shift_scalar = gaussian_shift_scalar
        self.cost = cost
        

        self._num_freq_channels = 1024
        self._layers = []
        self._num_outputs = 1 # Regression


    def create_graph(self):
        """create_graph

        Create the network graph for use in a tensorflow session
        """
        self.save_params()
        self._msg = '\rcreating network graph '; self._vprint(self._msg)

        
        tf.reset_default_graph()
        self.is_training = tf.placeholder(dtype = tf.bool, shape = [], name = 'is_training')
        
        with tf.variable_scope('keep_probs'):
            
            self._msg += '.'; self._vprint(self._msg)

            self.sample_keep_prob = tf.placeholder(self.dtype, name = 'sample_keep_prob')
            self.downsample_keep_prob = tf.placeholder(self.dtype, name = 'downsample_keep_prob')
        
        
        with tf.variable_scope('sample'):
            
            self._msg += '.'; self._vprint(self._msg)
            
            self.X = tf.placeholder(self.dtype, shape = [None, 1, self._num_freq_channels, 1], name = 'X')
            self.X = tf.nn.dropout(self.X, self.sample_keep_prob)
            
        trainable = lambda shape, name : tf.get_variable(name = name,
                                                         dtype = self.dtype,
                                                         shape = shape,
                                                         initializer = tf.contrib.layers.xavier_initializer())
        for i in range(self.num_downsamples):
            self._msg += '.'; self._vprint(self._msg)
            
            with tf.variable_scope('conv_layer_{}'.format(i)):

                layer = self.X if i == 0 else self._layers[-1]
                
                # filter shape:
                fh = 1 # filter height = 1 for 1D convolution
                fw = 3 # filter width
                fic = 2**(i) # num in channels = number of incoming filters
                foc = 2**(i+1) # num out channels = number of outgoing filters
                filters = trainable([fh, fw, fic, foc], 'filters')
                
                # stride shape
                sN = 1 # batch = 1 (why anything but 1 ?)
                sH = 1 # height of stride = 1 for 1D conv
                sW = 1 # width of stride = downsampling factor = 1 for no downsampling or > 1 for downsampling
                sC = 1 # depth = number of channels the stride walks over = 1 (why anything but 1 ?)
                strides_no_downsample = [sN, sH, sW, sC]
                layer = tf.nn.conv2d(layer, filters, strides_no_downsample, 'SAME')
                
                # shape of biases = [num outgoing filters]
                biases = trainable([foc], 'biases')
                layer = tf.nn.bias_add(layer, biases)
                layer = tf.nn.leaky_relu(layer)
                layer = tf.contrib.layers.batch_norm(layer, is_training = self.is_training)
                
                # downsample
                with tf.variable_scope('downsample'):
                    fw = 5
                    filters = trainable([fh, fw, foc, foc], 'filters')

                    sW = 2
                    strides_no_downsample = [sN, sH, sW, sC]
                    layer = tf.nn.conv2d(layer, filters, strides_no_downsample, 'SAME')

                    # shape of biases = [num outgoing filters]
                    biases = trainable([foc], 'biases')
                    layer = tf.nn.bias_add(layer, biases)
                    layer = tf.nn.leaky_relu(layer)
                    layer = tf.nn.dropout(layer, self.downsample_keep_prob)
                    layer = tf.contrib.layers.batch_norm(layer, is_training = self.is_training)

                self._layers.append(layer)
                              
        
                
        self._msg += ' '
        with tf.variable_scope('targets'):
            self._msg += '.'; self._vprint(self._msg)
            
            self.targets = tf.placeholder(dtype = self.dtype, shape = [None, self._num_outputs], name = 'labels')

        with tf.variable_scope('predictions'):
            self._msg += '.'; self._vprint(self._msg)
            
            self.predictions = tf.contrib.layers.fully_connected(tf.layers.flatten(self._layers[-1]), self._num_outputs)
            
        with tf.variable_scope('costs'):
            self._msg += '.'; self._vprint(self._msg)

            error = tf.subtract(self.targets, self.predictions, name = 'error')
            squared_error = tf.square(error, name = 'squared_difference')
            quad_error = tf.square(squared_error, name = 'quad_error' )

            with tf.variable_scope('mean_inverse_shifted_gaussian'):
                self._msg += '.'; self._vprint(self._msg)

                normal_dist = tf.contrib.distributions.Normal(0.0, self.accuracy_threshold, name = 'normal_dist')
                gaussian_prob = normal_dist.prob(error, name = 'gaussian_prob')
                shifted_gaussian = tf.add(gaussian_prob, self.gaussian_shift_scalar, name = 'shifted_gaussian')   

                self.MISG = tf.reduce_mean(tf.divide(1.0, shifted_gaussian), name = 'mean_inverse_shifted_gaussian')
                
            with tf.variable_scope('mean_squared_error'):
                self._msg += '.'; self._vprint(self._msg)
                
                self.MSE = tf.reduce_mean(squared_error / self.gaussian_shift_scalar)
                
            with tf.variable_scope('mean_quad_error'):
                self._msg += '.'; self._vprint(self._msg)
                
                self.MQE = tf.reduce_mean(quad_error / self.gaussian_shift_scalar)

        with tf.variable_scope('logging'):  

            with tf.variable_scope('image'):
                self._msg += '.'; self._vprint(self._msg)
                
                self.image_buf = tf.placeholder(tf.string, shape=[])
                epoch_image = tf.expand_dims(tf.image.decode_png(self.image_buf, channels=4), 0)

            with tf.variable_scope('percent_within_threshold'):
                self._msg += '.'; self._vprint(self._msg)
                
                self.PWT = 100.*tf.reduce_mean(tf.cast(tf.less_equal(tf.abs(self.targets - self.predictions), self.accuracy_threshold), self.dtype) )


            tf.summary.histogram(name = 'targets', values = self.targets)
            tf.summary.histogram(name = 'predictions',values =  self.predictions)
            tf.summary.scalar(name = 'MSE', tensor = self.MSE)
            tf.summary.scalar(name = 'MISG', tensor = self.MISG)
            tf.summary.scalar(name = 'MQE', tensor = self.MQE)
            tf.summary.scalar(name = 'PWT', tensor = self.PWT)
            tf.summary.image('prediction_vs_actual', epoch_image)
            self.summary = tf.summary.merge_all()
            
        update_ops = tf.get_collection(tf.GraphKeys.UPDATE_OPS)
        with tf.control_dependencies(update_ops):
            with tf.variable_scope('train'):
                self._msg += '.'; self._vprint(self._msg)
            
                if self.cost == 'MSE':
                    cost = self.MSE
                if self.cost == 'MQE':
                    cost = tf.log(self.MQE)
                if self.cost == 'MISG':
                    cost = self.MISG
                if self.cost == 'PWT_weighted_MSE':
                    cost = self.MSE * (100. - self.PWT)
                if self.cost == 'PWT_weighted_MISG':
                    cost = self.MISG * (100. - self.PWT)


                self.optimizer = tf.train.AdamOptimizer(self.adam_initial_learning_rate, epsilon=1e-08).minimize(cost)
            
        
        num_trainable_params = np.sum([np.prod(v.get_shape().as_list()) for v in tf.trainable_variables()])
        self._msg = '\rNetwork Ready - {} trainable parameters'.format(num_trainable_params); self._vprint(self._msg)<file_sep>/nsk/scripts/sky_image.py
#!/usr/bin/env python2.7
"""
sky_image.py
-------------
sky-based calibration with CASA 5.1.1

run the script as:
casa -c sky_image.py <args>

<NAME>
<EMAIL>
Nov. 2017
"""
import sys
import os
import numpy as np
import argparse
import subprocess
import shutil
import glob

## Set Arguments
# Required Arguments
a = argparse.ArgumentParser(description="Run with casa as: casa -c sky_image.py <args>")
a.add_argument('--script', '-c', type=str, help='name of this script', required=True)
a.add_argument('--msin', default=None, type=str, help='path to a CASA measurement set. if fed a .uvfits, will convert to ms', required=True)
a.add_argument('--source', default=None, type=str, help='source name', required=True)
# IO Arguments
a.add_argument('--out_dir', default=None, type=str, help='output directory')
a.add_argument("--silence", default=False, action='store_true', help="turn off output to stdout")
a.add_argument('--source_ext', default=None, type=str, help="Extension to default source name in output image files")
a.add_argument("--im_stem", default=None, type=str, help="Imagename stem for output images. Default is basename of input MS.")
# Calibration Arguments
a.add_argument("--model_im", default=None, type=str, help="path to model image, if None will look for a {source}.cl file")
a.add_argument('--refant', default=None, type=str, help='reference antenna')
a.add_argument('--ex_ants', default=None, type=str, help='bad antennas to flag')
a.add_argument('--rflag', default=False, action='store_true', help='run flagdata(mode=rflag)')
a.add_argument('--unflag', default=False, action='store_true', help='start by unflagging data')
a.add_argument('--KGcal', default=False, action='store_true', help='perform K (dly) & G (phs) calibration')
a.add_argument("--KGsnr", default=2.0, type=float, help="KG calibration Signal-to-Noise cut")
a.add_argument('--Acal', default=False, action='store_true', help='perform G (amp) calibration')
a.add_argument("--Asnr", default=2.0, type=float, help="G-amplitude calibration Signal-to-Noise cut")
a.add_argument('--BPcal', default=False, action='store_true', help='perform BandPass calibration (phs & amp)')
a.add_argument("--BPsnr", default=2.0, type=float, help="bandpass calibration Signal-to-Noise cut")
a.add_argument('--uvrange', default="", type=str, help="uvrange in meters (baseline length) to use in calibration and imaging")
a.add_argument('--timerange', default=[""], type=str, nargs='*', help="calibration and clean timerange(s)")
a.add_argument('--bpoly', default=False, action='store_true', help="use BPOLY mode in bandpass")
a.add_argument('--degamp', default=4, type=int, help="amplitude polynomial degree for BPOLY")
a.add_argument('--degphase', default=1, type=int, help="phase polynomial degree for BPOLY")
a.add_argument('--calspw', default='0:100~924', type=str, help="Calibration spectral window selection")
a.add_argument('--smodel', default=[], type=float, nargs='*', help="Stokes source model as I Q U V")
# Imaging Arguments
a.add_argument('--image_mfs', default=False, action='store_true', help="make an MFS image across the band")
a.add_argument('--niter', default=50, type=int, help='number of clean iterations.')
a.add_argument('--pxsize', default=300, type=int, help='pixel (cell) scale in arcseconds')
a.add_argument('--imsize', default=500, type=int, help='number of pixels along a side of the output square image.')
a.add_argument('--cleanspw', default="0:100~924", type=str, help="spectral window selection for clean")
a.add_argument('--image_model', default=False, action='store_true', help="image model datacolumn instead of data datacolumn")
a.add_argument('--spec_cube', default=False, action='store_true', help="image spectral cube as well as MFS.")
a.add_argument('--spec_dchan', default=40, type=int, help="number of channel averaging for a single image in the spectral cube.")
a.add_argument('--spec_start', default=100, type=int, help='starting channel for spectral cube')
a.add_argument('--spec_end', default=924, type=int, help='ending channel for spectral cube')
a.add_argument("--stokes", default='I', type=str, help="Stokes parameters to image.")
# Plotting Arguments
a.add_argument("--plot_uvdist", default=False, action='store_true', help='make a uvdist plot')


def echo(message, type=0):
    if verbose:
        if type == 0:
            print(message)
        elif type == 1:
            print("\n" + message + "\n" + "-"*40)

if __name__ == "__main__":
    # parse args
    args = a.parse_args()

    # get vars
    if args.source_ext is None:
        args.source_ext = ''
    verbose = args.silence is False

    # check for .loc and .cl files
    if os.path.exists('{}.loc'.format(args.source)) is False:
        raise AttributeError("{}.loc file doesn't exist in working directory".format(args.source))

    # configure refant
    if args.refant is None and (args.KGcal is True or args.Acal is True or args.BPcal is True):
        raise AttributeError("if calibrating, refant needs to be specified")
    args.refant = "HH" + str(args.refant)

    # get phase center
    ra, dec = np.loadtxt('{}.loc'.format(args.source), dtype=str)
    ra, dec = ra.split(':'), dec.split(':')
    fixdir = 'J2000 {}h{}m{}s {}d{}m{}s'.format(*(ra+dec))
    msin = args.msin

    # get paths
    base_ms = os.path.basename(msin)
    if args.out_dir is None:
        out_dir = os.path.dirname(msin)
    else:
        out_dir = args.out_dir

    # check for uvfits
    if base_ms.split('.')[-1] == 'uvfits':
        echo("...converting uvfits to ms", type=1)
        uvfits = msin
        msin = os.path.join(out_dir, '.'.join(base_ms.split('.')[:-1] + ['ms']))
        base_ms = os.path.basename(msin)
        msfiles = glob.glob("{}*".format(msin))
        if len(msfiles) != 0:
            for i, msf in enumerate(msfiles):
                try:
                    os.remove(msf)
                except OSError:
                    shutil.rmtree(msf)
        importuvfits(uvfits, msin)
        echo("{}".format(msin))

    # rephase to source
    echo("...fix vis to {}".format(fixdir), type=1)
    fixvis(msin, msin, phasecenter=fixdir)

    # insert source model
    if (args.KGcal is True or args.Acal is True or args.BPcal is True) or args.image_model is True:
        if args.model_im is None:
            echo("...inserting {} as MODEL".format("{}.cl".format(args.source)), type=1)
            ft(msin, complist="{}.cl".format(args.source), usescratch=True)
        else:
            echo("...inserting {} as MODEL".format(args.model_im), type=1)
            ft(msin, model=args.model_im, usescratch=True)

    # unflag
    if args.unflag is True:
        echo("...unflagging", type=1)
        flagdata(msin, mode='unflag')

    # flag autocorrs
    echo("...flagging autocorrs", type=1)
    flagdata(msin, autocorr=True)

    # flag bad ants
    if args.ex_ants is not None:
        args.ex_ants = ','.join(map(lambda x: "HH"+x, args.ex_ants.split(',')))
        echo("...flagging bad ants: {}".format(args.ex_ants), type=1)
        flagdata(msin, mode='manual', antenna=args.ex_ants)

    # rflag
    if args.rflag is True:
        echo("...rfi flagging", type=1)
        flagdata(msin, mode='rflag')

    def KGCAL(msin, gaintables=[]):
        ## perform per-antenna delay and phase calibration ##
        # setup calibration tables     
        kc = os.path.join(out_dir, base_ms+'.{}.cal'.format('K'))
        gpc = os.path.join(out_dir, base_ms+'.{}.cal'.format('Gphs'))

        # perform initial K calibration (per-antenna delay)
        echo("...performing K gaincal", type=1)
        if os.path.exists(kc):
            shutil.rmtree(kc)
        if os.path.exists("{}.png".format(kc)):
            os.remove("{}.png".format(kc))
        gaincal(msin, caltable=kc, gaintype="K", solint='inf', refant=args.refant, minsnr=args.KGsnr, spw=args.calspw,
                gaintable=gaintables, timerange=cal_timerange, uvrange=args.uvrange)
        plotcal(kc, xaxis='antenna', yaxis='delay', figfile='{}.png'.format(kc), showgui=False)
        gaintables.append(kc)

        # write delays as npz file
        tb.open(kc)
        delays = tb.getcol('FPARAM')[0, 0]
        delay_ants = tb.getcol('ANTENNA1')
        delay_flags = tb.getcol('FLAG')[0, 0]
        tb.close()
        np.savez("{}.npz".format(kc), delay_ants=delay_ants, delays=delays, delay_flags=delay_flags)
        echo("...Saving delays to {}.npz".format(kc))
        echo("...Saving plotcal to {}.png".format(kc))

        # perform initial G calibration for phase (per-spw and per-pol gain)
        echo("...performing G gaincal for phase", type=1)
        if os.path.exists(gpc):
            shutil.rmtree(gpc)
        if os.path.exists("{}.png".format(gpc)):
            os.remove("{}.png".format(gpc))
        gaincal(msin, caltable=gpc, gaintype='G', solint='inf', refant=args.refant, minsnr=args.KGsnr, calmode='p',
                spw=args.calspw, gaintable=gaintables, timerange=cal_timerange, uvrange=args.uvrange)
        plotcal(gpc, xaxis='antenna', yaxis='phase', figfile='{}.png'.format(gpc), showgui=False)
        gaintables.append(gpc)

        # write phase to file
        tb.open(gpc)
        phases = np.angle(tb.getcol('CPARAM')[0, 0])
        phase_ants = tb.getcol('ANTENNA1')
        phase_flags = tb.getcol('FLAG')[0, 0]
        tb.close()
        np.savez("{}.npz".format(gpc), phase_ants=phase_ants, phases=phases, phase_flags=phase_flags)
        echo("...Saving phases to {}.npz".format(gpc))
        echo("...Saving plotcal to {}.png".format(gpc))

        return gaintables

    def ACAL(msin, gaintables=[]):
        # gaincal G amplitude
        echo("...performing G gaincal for amplitude", type=1)
        gac = msin+'.{}.cal'.format('Gamp')
        if os.path.exists(gac):
            shutil.rmtree(gac)
        if os.path.exists("{}.png".format(gac)):
            os.remove("{}.png".format(gac))
        gaincal(msin, caltable=gac, gaintype='G', solint='inf', refant=args.refant, minsnr=args.Asnr, calmode='a',
                spw=args.calspw, gaintable=gaintables, timerange=cal_timerange, uvrange=args.uvrange)
        plotcal(gac, xaxis='antenna', yaxis='amp', figfile='{}.png'.format(gac), showgui=False)
        gaintables.append(gac)

        # write amp to file
        tb.open(gac)
        amps = np.abs(tb.getcol('CPARAM')[0, 0])
        amp_ants = tb.getcol('ANTENNA1')
        amp_flags = tb.getcol('FLAG')[0, 0]
        tb.close()
        np.savez("{}.npz".format(gac), amp_ants=amp_ants, amps=amps, amp_flags=amp_flags)
        echo("...Saving amps to {}.npz".format(gac))
        echo('...Saving G amp plotcal to {}.png'.format(gac))

        return gaintables

    def BPCAL(msin, gaintables=[]):
        # calibrated bandpass
        echo("...performing B bandpass cal", type=1)
        bc = msin+'.{}.cal'.format('B')
        Btype = "B"
        if args.bpoly:
            Btype="BPOLY"
        if os.path.exists(bc):
            shutil.rmtree(bc)
        if os.path.exists("{}.amp.png".format(bc)):
            os.remove("{}.amp.png".format(bc))
        if os.path.exists("{}.phs.png".format(bc)):
            os.remove("{}.phs.png".format(bc))
        bandpass(vis=msin, spw="", minsnr=args.BPsnr, bandtype=Btype, degamp=args.degamp, degphase=args.degphase,
                caltable=bc, gaintable=gaintables, solint='inf', refant=args.refant, timerange=cal_timerange,
                uvrange=args.uvrange, smodel=args.smodel)
        plotcal(bc, xaxis='chan', yaxis='amp', figfile="{}.amp.png".format(bc), showgui=False)
        plotcal(bc, xaxis='chan', yaxis='phase', figfile="{}.phs.png".format(bc), showgui=False)
        gaintables.append(bc)

        # write bp to file
        if args.bpoly is False:
            # get flags and bandpass data
            tb.open(bc)
            bp = tb.getcol('CPARAM')[0]
            bp_ants = tb.getcol("ANTENNA1")
            bp_flags = tb.getcol('FLAG')[0]
            tb.close()
            # load spectral window data
            tb.open(bc+"/SPECTRAL_WINDOW")
            bp_freqs = tb.getcol("CHAN_FREQ")
            tb.close()
            # write to file
            np.savez("{}.npz".format(bc), bp=bp, bp_ants=bp_ants, bp_flags=bp_flags, bp_freqs=bp_freqs)
            echo("...Saving bandpass to {}.npz".format(bc))
            echo("...Saving amp plotcal to {}.amp.png".format(bc))
            echo("...Saving phs plotcal to {}.phs.png".format(bc))
        else:
            echo("couldn't access bandpass data. Note BPOLY solutions not currently compatible w/ caltable2calfits.py")

        return gaintables

    if (args.KGcal is True or args.Acal is True or args.BPcal is True):
        ## Begin Calibration ##
        # init cal_timerange
        cal_timerange = ','.join(args.timerange)
        # run through various calibration options
        gaintables = []
        if args.KGcal:
            gaintables = KGCAL(msin, gaintables)

        if args.Acal:
            gaintables = ACAL(msin, gaintables)

        if args.BPcal:
            gaintables = BPCAL(msin, gaintables)

        # apply calibration gaintables
        echo("...applying gaintables: \n {}".format('\n'.join(gaintables)), type=1)
        applycal(msin, gaintable=gaintables)
        ms_split = os.path.join(out_dir, "{}.split".format(base_ms))
        files = glob.glob("{}*".format(ms_split))
        for f in files:
            if os.path.exists(f):
                try:
                    shutil.rmtree(f)
                except OSError:
                    os.remove(f)

        split(msin, ms_split, datacolumn="corrected")
        fixvis(ms_split, ms_split, phasecenter=fixdir)

    else:
        echo("...no calibration performed", type=1)
        ms_split = msin

    if args.image_mfs is True or args.spec_cube is True:
        # remove paths
        if args.im_stem is None:
            im_stem = os.path.join(out_dir, base_ms + '.' + args.source + args.source_ext)
        else:
            im_stem = args.im_stem
        echo("...performing clean for output files:\n{}.*".format(im_stem), type=1)
        source_files = glob.glob(im_stem+'*')
        if len(source_files) > 0:
            for sf in source_files:
                if os.path.exists(sf):
                    try:
                        shutil.rmtree(sf)
                    except OSError:
                        os.remove(sf)

    if args.image_model:
        model_ms_split = ms_split + ".model"
        if args.im_stem is None:
            model_im_stem = os.path.join(out_dir, base_ms + '.model.' + args.source + args.source_ext)
        else:
            model_im_stem = args.im_stem + '.model'
        if args.model_im is None:
            echo("...inserting {} as MODEL".format("{}.cl".format(args.source)), type=1)
            ft(ms_split, complist="{}.cl".format(args.source), usescratch=True)
        else:
            echo("...inserting {} as MODEL".format(args.model_im), type=1)
            ft(ms_split, model=args.model_im, usescratch=True)
        split(ms_split, model_ms_split, datacolumn='model')

    # create mfs image
    if args.image_mfs:
        echo("...running MFS clean", type=1)
        def image_mfs(msin, im_stem, timerange, cleanspw):
            clean(vis=msin, imagename=im_stem, spw=cleanspw, niter=args.niter, weighting='briggs', robust=0, imsize=[args.imsize, args.imsize],
                  cell=['{}arcsec'.format(args.pxsize)], mode='mfs', timerange=timerange, uvrange=args.uvrange, stokes=args.stokes)
            exportfits(imagename='{}.image'.format(im_stem), fitsimage='{}.fits'.format(im_stem))
            print("...saving {}".format('{}.fits'.format(im_stem)))

    
        for i, tr in enumerate(args.timerange):
            if i == 0:
                image_mfs(ms_split, im_stem, tr, args.cleanspw)
                if args.image_model:
                    image_mfs(model_ms_split, model_im_stem, tr, args.cleanspw)
            else:
                image_mfs(ms_split, im_stem+'_tr{}'.format(i), tr, args.cleanspw)
                if args.image_model:
                    image_mfs(model_ms_split, model_im_stem+'_tr{}'.format(i), tr, args.cleanspw)

    # create spectrum
    if args.spec_cube:
        echo("...running spectral cube clean", type=1)
        def spec_cube(msin, im_stem, timerange):
            dchan = args.spec_dchan
            for i, chan in enumerate(np.arange(args.spec_start, args.spec_end, dchan)):
                clean(vis=msin, imagename=im_stem+'.spec{:04d}'.format(chan), niter=0, spw="0:{}~{}".format(chan, chan+dchan-1),
                        weighting='briggs', robust=0, imsize=[args.imsize, args.imsize], timerange=timerange, uvrange=args.uvrange,
                        cell=['{}arcsec'.format(args.pxsize)], mode='mfs', stokes=args.stokes)#, mask='circle[[{}h{}m{}s, {}d{}m{}s ], 7deg]'.format(*(ra+dec)))
                exportfits(imagename='{}.spec{:04d}.image'.format(im_stem, chan), fitsimage='{}.spec{:04d}.fits'.format(im_stem, chan))
                print("...saving {}".format('{}.spec{:04d}.fits'.format(im_stem, chan)))

        for i, tr in enumerate(args.timerange):
            if i == 0:
                spec_cube(ms_split, im_stem, tr)
                if args.image_model:
                    spec_cube(model_ms_split, model_im_stem, tr)
            else:
                spec_cube(ms_split, im_stem+'_tr{}'.format(i), tr)
                if args.image_model:
                    spec_cube(model_ms_split, model_im_stem+'_tr{}'.format(i), tr)

    # make uvdist plot
    if args.plot_uvdist:
        echo("...plotting uvdistance", type=1)
        # add model to ms_split
        if args.model_im is None:
            ft(ms_split, complist="{}.cl".format(args.source), usescratch=True)
        else:
            ft(ms_split, model=args.model_im, usescratch=True)
        # load visibility amplitudes
        ms.open(ms_split)
        data = ms.getdata(["amplitude", "antenna1", "antenna2", "uvdist", "axis_info", "flag", "model_amplitude"], ifraxis=True)
        amps = data['amplitude']
        mamps = data['model_amplitude']
        flags = data['flag']
        uvdist = data['uvdist']
        freqs = data['axis_info']['freq_axis']['chan_freq']
        ms.close()
        # get rid of autos
        select = []
        for i, a1 in enumerate(data['antenna1']):
            if a1 != data['antenna2'][i]:
                select.append(i)
        amps = amps[:, :, select, :].squeeze()
        mamps = mamps[:, :, select, :].squeeze()
        uvdist = uvdist[select, :]
        flags = flags[:, :, select, :].squeeze()
        # omit flagged data
        amps[flags] *= np.nan
        mamps[flags] *= np.nan
        # average across time
        amps = np.nanmean(amps, axis=2)
        mamps = np.nanmean(mamps, axis=2)
        uvdist = np.nanmean(uvdist, axis=1)
        # average into channel bins
        freq_bins = np.median(np.split(np.linspace(100., 200., 1024, endpoint=True)[100:924], 4), axis=1)
        amps = np.nanmedian(np.split(amps[100:924, :], 4), axis=1)
        mamps = np.nanmedian(np.split(mamps[100:924, :], 4), axis=1)
        # plot
        import matplotlib.pyplot as plt
        def plot_uvdist(amp, ext):
            fig, ax = plt.subplots(1, 1, figsize=(10,7))
            ax.grid(True)
            p = ax.plot(uvdist, amp.T, marker='o', markersize=6, alpha=0.8, ls='')
            ax.set_xlabel("baseline length (meters)", fontsize=16)
            ax.set_ylabel("amplitude (arbitrary units)", fontsize=16)
            ax.set_title("{} for {}".format(ext, base_ms))
            ax.legend(p, map(lambda x: '{:0.0f} MHz'.format(x), freq_bins))
            file_fname = os.path.join(out_dir, "{}.{}.png".format(base_ms, ext))
            echo("...saving {}".format(file_fname))
            fig.savefig(file_fname, bbox_inches='tight', pad=0.05)
            plt.close()

        plot_uvdist(amps, "data_uvdist")
        plot_uvdist(mamps, "model_uvdist")



<file_sep>/aas/Estimating Delays/nn/data_creation/Data_Creator.py
"""Data_Creator

base class for data creators
"""

import sys, os
from data_manipulation import *

import numpy as np

class Data_Creator(object):
    """Data_Creator - Parent for network style-specific creators

        Generates data, in an alternate thread, for training.
        Child MUST implement _gen_data()
    
        Args:
            num_flatnesses(int): number of flatnesses used to generate data.
                                -   Number of data samples = 60 * num_flatnesses
            bl_data : data source. Output of get_seps_data()
            bl_dict (dict) - Dictionary of seps with bls as keys. An output of get_or_gen_test_train_red_bls_dicts()
            gains(dict): - Gains for this data. An output of load_relevant_data()
    
    """

    def __init__(self,
                 num_flatnesses,
                 bl_data,
                 bl_dict,
                 gains,
                 abs_min_max_delay = 0.040):     
        """__init__
        
        Args:
            num_flatnesses (int): Number of flatnesses to generate
            bl_data (dict): keys are baselines, data is complex visibilities
            bl_dict (dict): keys are unique redundant baseline, values are baselines in that redudany group
            gains (dict): complex antenna gains
            abs_min_max_delay (float, optional): MinMax value of delay
        """
        self._num = num_flatnesses
                    
        self._bl_data = bl_data
        self._bl_data_c = None
        
        self._bl_dict = bl_dict
        
        self._gains = gains
        self._gains_c = None
        
        self._epoch_batch = []
        self._nu = np.arange(1024)
        self._tau = abs_min_max_delay
        
    def _gen_data(self):
        """_gen_data

        This method is called in a separate thread via gen_data().
        Must be overridden
        """
        print('Must implement _gen_data()')
 
    def gen_data(self):
        """gen_data

        Starts a new thread and generates data there.
        """
        
        self._thread = Thread(target = self._gen_data, args=())
        self._thread.start()

    def get_data(self, timeout = 10):
        """get_data

        Retrieves the data from the thread.
  
        Args:
            timeout (int, optional): data generation timeout in seconds
        
        Returns:
            (list of complex floats): shape = (num_flatnesses, 60, 1024)
             - needs to be reshaped for training

        """
        
        if len(self._epoch_batch) == 0:
            self._thread.join(timeout)
            
        return self._epoch_batch.pop(0)<file_sep>/aas/Estimating Delays/nn/network_trainers/CNN_R_Trainer.py
# CNN_R_Trainer

import sys, os
sys.path.insert(1, os.path.join(sys.path[0], '../modules'))
from NN_Trainer import NN_Trainer

import tensorflow as tf
from tensorflow.python.client import timeline

import numpy as np

class CNN_R_Trainer(NN_Trainer):
    
    def __init__(self,
                 network,
                 Data_Creator, # Class
                 num_epochs = 100,
                 batch_size = 32,
                 log_dir = '../logs/',
                 model_save_interval = 25,
                 pretrained_model_path = None,
                 metric_names = ['MISGs', 'MSEs', 'PWTs'],
                 sample_keep_prob = 0.80,
                 conv_keep_prob = 0.9,
                 pred_keep_prob = 0.50,
                 verbose = True):
    
        NN_Trainer.__init__(self,
                            network = network,
                            Data_Creator = Data_Creator,
                            num_epochs = num_epochs,
                            batch_size = batch_size,
                            log_dir = log_dir,
                            model_save_interval = model_save_interval,
                            pretrained_model_path = pretrained_model_path,
                            metric_names = metric_names,
                            verbose = verbose)
        

        self.sample_keep_prob = sample_keep_prob
        self.conv_keep_prob = conv_keep_prob
        self.pred_keep_prob = pred_keep_prob
        
    def train(self):
        
        self.save_params()
        
        

        MISGs = []
        MSEs = []
        PWTs = []
        
        tf.reset_default_graph()
        
        self._network.create_graph()
        saver = tf.train.Saver()
        
        itx =  lambda x: np.array(x) * 2. * self._abs_min_max_delay - self._abs_min_max_delay
        
        
        with tf.Session() as session:

            if self.pretrained_model_path == None:
                session.run(tf.global_variables_initializer())
                
            else:
                saver.restore(session, self.pretrained_model_path)
                
            options = tf.RunOptions(trace_level = tf.RunOptions.FULL_TRACE)
            run_metadata = tf.RunMetadata()

            archive_loc = self.log_dir + self._network.name 
            training_writer = tf.summary.FileWriter(archive_loc + '/training', session.graph)
            testing_writer = tf.summary.FileWriter(archive_loc + '/testing', session.graph)
            self._model_save_location = archive_loc + '/trained_model.ckpt'   
            
            try:
                for epoch in range(self.num_epochs):




                    training_inputs, training_targets = self._train_batcher.get_data(); self._train_batcher.gen_data()
                    testing_inputs, testing_targets = self._test_batcher.get_data(); self._test_batcher.gen_data()  

                    # if the division here has a remainde some values are just truncated
                    batch_size = self.batch_size
                    num_entries = training_inputs.shape[0]

                    for j in range(int(num_entries/batch_size)):
                        self._msg = '\repoch'
                        self._msg += '- {:5.0f}/{:5.0f}'.format(epoch + 1,self.num_epochs)
                        self._msg += ' - batch: {:4.0f}/{:4.0f}'.format(j + 1, int(num_entries/batch_size))
                        if epoch != 0:
                            self._msg += ' - (Training, Testing) - '.format(epoch)
                            self._msg += ' MISG: ({:0.4f}, {:0.4f})'.format(training_MISG, testing_MISG)
                            self._msg += ' MSE: ({:0.4f}, {:0.4f})'.format(training_MSE, testing_MSE)
                            self._msg += ' PWT: ({:2.2f}, {:2.2f})'.format(training_PWT, testing_PWT)



                        training_inputs_batch = training_inputs[j*batch_size:(j + 1)*batch_size].reshape(-1,1,1024,1)
                        training_targets_batch = training_targets[j*batch_size:(j + 1)*batch_size].reshape(-1,1)

                        feed_dict = {self._network.samples: training_inputs_batch,
                                     self._network.targets: training_targets_batch,
                                     self._network.sample_keep_prob : self.sample_keep_prob,
                                     self._network.conv_keep_prob : self.conv_keep_prob,
                                     self._network.pred_keep_prob : self.pred_keep_prob,
                                     self._network.is_training : True}



                        if j == 0 and (epoch + 1) % self.model_save_interval == 0:

                            session.run([self._network.optimizer], feed_dict = feed_dict,
                                           options = options, run_metadata = run_metadata) 

                            training_writer.add_run_metadata(run_metadata, 'epoch%d' %epoch)

                            fetched_timeline = timeline.Timeline(run_metadata.step_stats)
                            chrome_trace = fetched_timeline.generate_chrome_trace_format()
                            direc = archive_loc + '/timelines/'

                            if not os.path.exists(direc):
                                os.makedirs(direc)
                            with open(direc + 'timeline_{}.json'.format(epoch), 'w') as f:
                                f.write(chrome_trace)
                                self._msg += ' saving timeline & metadata'
                        else:
                            session.run([self._network.optimizer], feed_dict = feed_dict) 

                        self._vprint(self._msg)

                    # Prediction: Scaled Train(ing results)   
                    PST = session.run(self._network.predictions,
                                      feed_dict = {self._network.samples: training_inputs.reshape(-1,1,1024,1),
                                                   self._network.sample_keep_prob : 1.,
                                                   self._network.conv_keep_prob : 1.,
                                                   self._network.pred_keep_prob : 1.,
                                                   self._network.is_training : False}) 

                    train_feed_dict = {self._network.samples: training_inputs.reshape(-1,1,1024,1),
                                       self._network.targets: training_targets.reshape(-1,1),
                                       self._network.sample_keep_prob : 1.,
                                       self._network.conv_keep_prob : 1.,
                                       self._network.pred_keep_prob : 1.,
                                       self._network.image_buf: self.gen_plot(PST, training_targets, itx),
                                       self._network.is_training : False}


                    training_MISG, training_MSE, training_PWT, training_summary = session.run([self._network.MISG,
                                                                                               self._network.MSE,
                                                                                               self._network.PWT,
                                                                                               self._network.summary],
                                                                                              feed_dict = train_feed_dict) 

                    training_writer.add_summary(training_summary, epoch)
                    training_writer.flush()  

                    # Prediction: Scaled test(ing results)   
                    PST = session.run(self._network.predictions,
                                      feed_dict = {self._network.samples: testing_inputs.reshape(-1,1,1024,1),
                                                   self._network.sample_keep_prob : 1.,
                                                   self._network.conv_keep_prob : 1.,
                                                   self._network.pred_keep_prob : 1.,
                                                   self._network.is_training : False}) 

                    test_feed_dict = {self._network.samples: testing_inputs.reshape(-1,1,1024,1),
                                      self._network.targets: testing_targets.reshape(-1,1),
                                      self._network.sample_keep_prob : 1.,
                                      self._network.conv_keep_prob : 1.,
                                      self._network.pred_keep_prob : 1.,
                                      self._network.image_buf: self.gen_plot(PST,testing_targets, itx),
                                      self._network.is_training : False} 

                    testing_MISG, testing_MSE, testing_PWT, testing_summary = session.run([self._network.MISG,
                                                                                           self._network.MSE,
                                                                                           self._network.PWT,
                                                                                           self._network.summary],
                                                                                          feed_dict = test_feed_dict)




                    testing_writer.add_summary(testing_summary, epoch)
                    testing_writer.flush()  

                    MISGs.append((training_MISG, testing_MISG))
                    MSEs.append((training_MSE, testing_MSE))
                    PWTs.append((training_PWT, testing_PWT))

                    if (epoch + 1) % self.model_save_interval == 0:
                        saver.save(session, self._model_save_location, epoch)

            except KeyboardInterrupt:
                
                pass

            self._msg += '\rtraining ended'; self._vprint(self._msg)

            training_writer.close()
            testing_writer.close()

        session.close()
        self._msg += ' - session closed'; self._vprint(self._msg)
        
        self._metrics = [MISGs, MSEs, PWTs]
        self.save_metrics()
        <file_sep>/aas/Estimating Delays/nn/__init__.py
from networks import *
from network_trainers import *
from data_creation import *
<file_sep>/aas/Estimating Delays/nn/data_creation/Data_Creator_BC.py
# Data_Creator_BC

from data_manipulation import *
from Data_Creator import Data_Creator

import numpy as np

class Data_Creator_BC(Data_Creator):
    """Creates data in an alternate thread. BC for binary classifier.
       look at labels variable in _gen_data()
    
    ## usage:
    ## data_maker = Data_Creator_BC(num_flatnesses=250, mode = 'train')
    ## data_maker.gen_data() #before loop
    ## inputs, targets = data_maker.get_data() #start of loop
    ## data_maker.gen_data() #immediately after get_data()

    """

    def __init__(self,
                 num_flatnesses,
                 bl_data = None,
                 bl_dict = None,
                 gains = None,
                 abs_min_max_delay = 0.040):
        
        """
        Arguments
            num_flatnesses : int - number of flatnesses used to generate data.
                                   Number of data samples = 60 * num_flatnesses
            bl_data : data source. Output of get_seps_data()
            bl_dict : dict - Dictionary of seps with bls as keys. An output of get_or_gen_test_train_red_bls_dicts()
            gains : dict - Gains for this data. An output of load_relevant_data()
            
                                   
        """
        Data_Creator.__init__(self,
                              num_flatnesses = num_flatnesses,
                              bl_data = bl_data,
                              bl_dict = bl_dict,
                              gains = gains,
                              abs_min_max_delay = abs_min_max_delay)
        
             
    def _gen_data(self):
        
        # scaling tools
        # the NN likes data in the range (0,1)
        angle_tx  = lambda x: (np.asarray(x) + np.pi) / (2. * np.pi)
        angle_itx = lambda x: np.asarray(x) * 2. * np.pi - np.pi

        delay_tx  = lambda x: (np.array(x) + self._tau) / (2. * self._tau)
        delay_itx = lambda x: np.array(x) * 2. * self._tau - self._tau
        
        targets = np.random.uniform(low = -self._tau, high = self._tau, size = (self._num * 60, 1))
        applied_delay = np.exp(-2j * np.pi * (targets * self._nu + np.random.uniform()))



        assert type(self._bl_data) != None, "Provide visibility data"
        assert type(self._bl_dict) != None, "Provide dict of baselines"
        assert type(self._gains)   != None, "Provide antenna gains"

        if self._bl_data_c == None:
            self._bl_data_c = {key : self._bl_data[key].conjugate() for key in self._bl_data.keys()}

        if self._gains_c == None:
            self._gains_c = {key : self._gains[key].conjugate() for key in self._gains.keys()}


        def _flatness(seps):
            """Create a flatness from a given pair of seperations, their data & their gains."""

            a, b = seps[0][0], seps[0][1]
            c, d = seps[1][0], seps[1][1]


            return self._bl_data[seps[0]]   * self._gains_c[(a,'x')] * self._gains[(b,'x')] * \
                   self._bl_data_c[seps[1]] * self._gains[(c,'x')]   * self._gains_c[(d,'x')]

        inputs = []
        for _ in range(self._num):

            unique_baseline = random.sample(self._bl_dict.keys(), 1)[0]
            two_seps = [random.sample(self._bl_dict[unique_baseline], 2)][0]

            inputs.append(_flatness(two_seps))
            

        inputs = np.angle(np.array(inputs).reshape(-1,1024) * applied_delay)
        
        permutation_index = np.random.permutation(np.arange(self._num * 60))
        labels = [[1,0] if d >=0 else [0,1] for d in targets[permutation_index]]

        self._epoch_batch.append((angle_tx(inputs[permutation_index]), labels))

<file_sep>/aas/Estimating Delays/nn/data_creation/data_manipulation.py
"""data_manipulation

Collection of data helper files
"""
import numpy as np
from pyuvdata import UVData
import hera_cal as hc
import random
from threading import Thread
from glob import glob
import os

def load_relevant_data(miriad_path, calfits_path):
    """Loads redundant baselines, gains, and data.
    
    Arguments:
        miriad_path (str): path to a miriad file for some JD
        calfits_path (str): path to a calfits file for the same JD
    
    Returns:
        red_bls, gains, uvd
    
    """

    # read the data
    uvd = UVData()
    uvd.read_miriad(miriad_path)

    # get the redundancies for that data
    aa = hc.utils.get_aa_from_uv(uvd)
    info = hc.omni.aa_to_info(aa)
    red_bls = np.array(info.get_reds())

    # gains for same data 
    gains, _ = hc.io.load_cal(calfits_path)
    
    return red_bls, gains, uvd

def get_good_red_bls(red_bls, gain_keys, min_group_len = 4):
    """get_good_red_bls

    Select all the good antennas from red_bls
    
    Each baseline group in red_bls will have its bad separations removed.
    Groups with less than min_group_len are removed.
    
    Arguments:
        red_bls (list of lists of tuples of ints)
             - Each sublist is a group of redundant separations for a unique baseline.
        gain_keys (dict): gains.keys() from hc.io.load_cal()
        min_group_len (int, optional): Minimum number of separations in a 'good' sublist.
            - Default = 4, so that both training and testing can take two seps)
    
    Returns:
        list of lists: Each sublist is a len >=4 list of separations of good antennas
    
    """
    
    def ants_good(sep):
        """ants_good

        Returns True if both antennas are good.
        
        Because we are using data from firstcal
        (is this right? I have trouble rememebring the names
        and properties of the different data sources)
        we can check for known good or bad antennas by looking to see if the antenna
        is represented in gains.keys(). If the antenna is present, its good.
        
        Arguments:
            sep (tuple of ints): antenna indices
        
        Returns:
            bool: True if both antennas are in gain_keys else False
        """

        ants = [a[0] for a in gain_keys]

        if sep[0] in ants and sep[1] in ants:
            return True

        else:
            return False

    good_redundant_baselines = []
    
    for group in red_bls:
        
        new_group = []
        
        for sep in group:
            
            # only retain seps made from good antennas
            if ants_good(sep) == True:
                new_group.append(sep)
                
        new_group_len = len(new_group)
        
        # make sure groups are large enough that both the training set
        # and the testing set can take two seps 
        if new_group_len >= min_group_len:
            
            good_redundant_baselines.append(sorted(new_group))
            
    return good_redundant_baselines

def _train_test_split_red_bls(red_bls, training_percent = 0.80):
    """_train_test_split_red_bls

    Split a list of redundant baselines into a training set and a testing set. 
    
    Each of the two sets has at least one pair of baselines from every group from red_bls.
    However, separations from one set will not appear the other.
    
    Arguments:
        red_bls (list of lists): Each sublist is a group of redundant separations
                                  for a unique baseline.
                                  Each sublist must have at least 4 separations.
        training_percent (float, optional): - *Approximate* portion of the separations that will
                                   appear in the training set.
    
    Returns:
        dict, dict :train_red_bls_dict, test_red_bls_dict
    
    """
    
    
    training_redundant_baselines_dict = {}
    testing_redundant_baselines_dict = {}

    # make sure that each set has at least 2 seps from each group
    #
    thinned_groups_dict = {}
    for group in red_bls:

        # group key is the sep with the lowest antenna indicies
        key = sorted(group)[0]
        
        # pop off seps from group and append them to train or test groups
        random.shuffle(group)
        training_group = []
        training_group.append(group.pop())
        training_group.append(group.pop())

        testing_group = []
        testing_group.append(group.pop())
        testing_group.append(group.pop())
        
        # add the new train & test groups into the dicts
        training_redundant_baselines_dict[key] = training_group
        testing_redundant_baselines_dict[key] = testing_group
        
        # if there are still more seps in the group, save them into a dict for later assignment
        if len(group) != 0:
            thinned_groups_dict[key] = group

    # Shuffle and split the group keys into two sets using 
    #
    thinned_dict_keys = thinned_groups_dict.keys()
    random.shuffle(thinned_dict_keys)
    
    # Because we are ensuring that each set has some seps from every group,
    # the ratio of train / test gets reduced a few percent.
    # This (sort of) accounts for that with an arbitrary shift found by trial and error.
    #    
    #    Without this the a setting of training_percent = 0.80 results in a 65/35 split, not 80/20.
    #    
    #    I assume there is a better way...
    t_pct = np.min([0.95, training_percent + 0.15])

    
    # these 'extra' keys are the keys that each set will extract seps from thinned_groups_dict with
    training_red_bls_extra, testing_red_bls_extra = np.split(thinned_dict_keys,
                                                             [int(len(thinned_dict_keys)*t_pct)])

    # extract seps from thinned_groups_dict and apply to same key in training set
    for key in training_red_bls_extra:
        key = tuple(key)
        group = thinned_groups_dict[key]
        training_group = training_redundant_baselines_dict[key]
        training_group.extend(group)
        training_redundant_baselines_dict[key] = training_group

    # extract seps from thinned_groups_dict and apply to same key in testing set
    for key in testing_red_bls_extra:
        key = tuple(key)
        group = thinned_groups_dict[key]
        testing_group = testing_redundant_baselines_dict[key]
        testing_group.extend(group)
        testing_redundant_baselines_dict[key] = testing_group
        
    return training_redundant_baselines_dict, testing_redundant_baselines_dict

def _loadnpz(filename):
    """_loadnpz

    Loads up npzs. For dicts do loadnpz(fn)[()]
    
    Args:
        filename (str): path to npz to load
    
    Returns:
        numpy array: the data in the npz
    """
    
    a = np.load(filename)
    d = dict(zip(("data1{}".format(k) for k in a), (a[k] for k in a)))
    
    return d['data1arr_0']

def get_or_gen_test_train_red_bls_dicts(red_bls = None,
                                        gain_keys = None,
                                        training_percent = 0.80,
                                        training_load_path = None,
                                        testing_load_path = None):
    """get_or_gen_test_train_red_bls_dicts

    Create or generate baseline dicts for training and testing
    
    Args:
        red_bls (list of lists of tuples of ints): List of redundant baselines
        gain_keys (list of int): antenna gains
        training_percent (float, optional): Approximate portion of data to assign as training data
        training_load_path (None, optional): Path to prepared dict of baselines
        testing_load_path (None, optional): Path to prepared dict of baselines
    
    Returns:
        dict, dict: training_red_bls_dict, testing_red_bls_dict
    """
    training_fn = 'data/training_redundant_baselines_dict_{}'.format(int(100*training_percent))
    testing_fn = 'data/testing_redundant_baselines_dict_{}'.format(int(100*training_percent))

    if training_load_path != None and testing_load_path != None:
        
        training_red_bls_dict = _loadnpz(training_load_path)[()]
        testing_red_bls_dict = _loadnpz(testing_load_path)[()]

    else:

        # if the files exist dont remake them.
        if os.path.exists(training_fn + '.npz') and os.path.exists(testing_fn + '.npz'):
            training_red_bls_dict = _loadnpz(training_fn + '.npz')[()]
            testing_red_bls_dict = _loadnpz(testing_fn + '.npz')[()]

        else:
        
            assert type(red_bls) is not None, "Provide a list of redundant baselines"
            assert type(gain_keys) is not None, "Provide a list of gain keys"
            
            good_red_bls = get_good_red_bls(red_bls, gain_keys)
            training_red_bls_dict, testing_red_bls_dict = _train_test_split_red_bls(good_red_bls,
                                                                                    training_percent = training_percent)

            np.savez(training_fn, training_red_bls_dict)
            np.savez(testing_fn, testing_red_bls_dict)

    return training_red_bls_dict, testing_red_bls_dict

def get_seps_data(red_bls_dict, uvd):
    """get_seps_data

    Get the data for all the seps in a redundant baselines dictionary.
    
    Args:
        red_bls_dict (dict): leys are unique baselines, values are tuples of antenna indicies
        uvd (pyuvdata.UVData: 
    
    Returns:
        dict: keys are separations (int, int), values are complex visibilitoes
    """
    
    data = {}
    for key in red_bls_dict.keys():
        for sep in red_bls_dict[key]:
            data[sep] = uvd.get_data(sep)
    
    return data

<file_sep>/aas/Estimating Delays/nn/networks/__init__.py
from RestoreableComponent import RestoreableComponent
from CNN_DS_BN_R import CNN_DS_BN_R
from CNN_DS_BN_C import CNN_DS_BN_C
from CNN_DS_BN_BC import CNN_DS_BN_BC
from CNN_DSFC_BN_R import CNN_DSFC_BN_R
from CNN_QP_BN_R import CNN_QP_BN_R
from FNN_BN_R import FNN_BN_R<file_sep>/aas/Estimating Delays/nn/data_creation/__init__.py
from data_manipulation import *
from Data_Creator_C import Data_Creator_C
from Data_Creator_R import Data_Creator_R
from Data_Creator_BC import Data_Creator_BC
from Data_Creator_BC_Odd_Even import Data_Creator_BC_Odd_Even
<file_sep>/aas/Estimating Delays/estdel/setup.py
from setuptools import setup

setup(name='estdel',
      version='0.0.1',
      description='Estimate interferometer anteanna delays',
      url='https://github.com/andrewasheridan/hera_sandbox/tree/master/aas/Estimating%20Delays/estdel',
      author='<NAME>',
      author_email='<EMAIL>',
      license='MIT',
      package_dir = {'estdel' : 'estdel'},
      packages=['estdel'],
      install_requires=[
          'numpy>=1.2',
          'tensorflow>=1.8.0',
      ],
      zip_safe=False,
      include_package_data=True)

# AIPY/BLOB/V3 SETUP.PY from import to the function<file_sep>/aas/Estimating Delays/nn/networks/RestoreableComponent.py
"""RestoreableComponent
"""

from pprint import pprint
import numpy as np
import sys
import os


class RestoreableComponent(object):
    """RestoreableComponent
    
    Allows saving / loading of object parameters.
    Excludes saving / loading of parameters begining with '_' or parameters that refer to tensorflow objects.
    Used for networks and network-trainers to ease recordkeeping.
    
    Args:
        name (str): Name of the object
        log_dir (str, optional): location to store object parameters.
        verbose (bool, optional): Be verbose
    
    Methods:
        print_params(): Print the parameters for the object
        save_params(): Saves parameters in log_dir/name/params/object_class.npz
        load_params(path): Loads parameters from path into existing object
                            - Retains the name of the existing object       
    
    Attributes:
        log_dir (str): store params in log_dir/params/CLASS_NAME
        name (str): Name of object
    """
    
    def __init__(self, name, log_dir = 'logs/', verbose = True):
        """__init__
        
        Args:
            name (str): Name of object
            log_dir (str, optional): store params in log_dir/params/CLASS_NAME
            verbose (bool, optional): Be verbose
        """
        self.name = name
        self.log_dir = log_dir
        self._vprint = sys.stdout.write if verbose else lambda *a, **k: None # for verbose printing
        self._msg = ''
        
    def _gen_params_dict(self):
        """_gen_params_dict

        Generate dict of class attributes.
        Excludes parameters begininning with '_' or parameters that refer to tensorflow objects.
        
        Returns:
            dict: keys are attribute names, values are their values
        """
        d = self.__dict__
        return {key : d[key] for key in d.keys() if key[0] != '_' if 'tensorflow' not in str(type(d[key]))}
    
    def print_params(self):
        """Prints restoreable parameters"""
        pprint(self._gen_params_dict())

    def load_params(self, path):
        """load_params

        Load in parameters from npz file.
        Do not include the .npz suffix in path
        Keeps the current object name.
        
        Args:
            path (str): path to restoreabale parameters
        """
    
        self._msg = '\rloading parameters';self._vprint(self._msg)
        
        a = np.load(path + '.npz')
        d = dict(zip(("data1{}".format(k) for k in a), (a[k] for k in a)))
        
        name = self.name
        params = d['data1arr_0'][()]
        for key in params:
            setattr(self, key, params[key])
        self.name = name  
        
        self._msg = '\rparams loaded';self._vprint(self._msg)
        
    def save_params(self, direc = None):
        """save_params

        Saves the restoreable parameters.
        Setting direc will override the default location of log_dir/params/
        
        Args:
            direc (str, optional): storeage directory
        """
        
        self._msg += '\rsaving parameters';self._vprint(self._msg)
        
        direc = self.log_dir + self.name + '/params/' if direc is None else direc
        if not os.path.exists(direc):
            
            self._msg += '- creating new directory';self._vprint(self._msg)
            
            os.makedirs(direc)
        np.savez(direc + self.__class__.__name__, self._gen_params_dict())  
        
        self._msg += ' - params saved';self._vprint(self._msg)<file_sep>/nsk/scripts/source_extract.py
#!/usr/bin/env python2.7
"""
source_extract.py
========

get image statistics of FITS file(s)
on a specific source
"""
import astropy.io.fits as fits
from astropy import modeling as mod
import argparse
import os
import sys
import shutil
import glob
import numpy as np
import scipy.stats as stats

a = argparse.ArgumentParser(description="Get FITS image statistics around source at center of image")

a.add_argument("files", type=str, nargs='*', help="filename(s) or glob-parseable string of filename(s)")
a.add_argument("--source", type=str, help="source name, with a <source>.loc file in working directory")
a.add_argument("--radius", type=float, default=2, help="radius in degrees around estimated source position to get source peak")
a.add_argument("--rms_max_r", type=float, default=None, help="max radius in degrees around source to make rms calculation")
a.add_argument("--rms_min_r", type=float, default=None, help="min radius in degrees around source to make rms calculation")
a.add_argument("--outdir", type=str, default=None, help="output directory")
a.add_argument("--ext", type=str, default=".spectrum.tab", help='extension string for spectrum file')
a.add_argument("--overwrite", default=False, action='store_true', help='overwite output')
a.add_argument("--gaussfit_mult", default=1.0, type=float, help="gaussian fit mask area is gaussfit_mult * synthesized_beam")

def source_extract(imfile, source, radius=1, gaussfit_mult=1.0, rms_max_r=None, rms_min_r=None, **kwargs):

    # open fits file
    hdu = fits.open(imfile)

    # get header and data
    head = hdu[0].header
    data = hdu[0].data.squeeze()

    # determine if freq precedes stokes in header
    if head['CTYPE3'] == 'FREQ':
        freq_ax = 3
        stok_ax = 4
    else:
        freq_ax = 4
        stok_ax = 3

    # get axes info
    npix1 = head["NAXIS1"]
    npix2 = head["NAXIS2"]
    nstok = head["NAXIS{}".format(stok_ax)]
    nfreq = head["NAXIS{}".format(freq_ax)]

    # calculate beam area in degrees^2
    beam_area = (head["BMAJ"] * head["BMIN"] * np.pi / 4 / np.log(2))

    # calculate pixel area in degrees ^2
    pixel_area = np.abs(head["CDELT1"] * head["CDELT2"])
    Npix_beam = beam_area / pixel_area
    
    # get ra dec coordiantes
    ra_axis = np.linspace(head["CRVAL1"]-head["CDELT1"]*head["NAXIS1"]/2, head["CRVAL1"]+head["CDELT1"]*head["NAXIS1"]/2, head["NAXIS1"])
    dec_axis = np.linspace(head["CRVAL2"]-head["CDELT2"]*head["NAXIS2"]/2, head["CRVAL2"]+head["CDELT2"]*head["NAXIS2"]/2, head["NAXIS2"])
    RA, DEC = np.meshgrid(ra_axis, dec_axis)

    # get source coordinates
    ra, dec = np.loadtxt('{}.loc'.format(source), dtype=str)
    ra, dec = map(float, ra.split(':')), map(float,dec.split(':'))
    ra = (ra[0] + ra[1]/60. + ra[2]/3600.) * 15
    dec = (dec[0] + np.sign(dec[0])*dec[1]/60. + np.sign(dec[0])*dec[2]/3600.)

    # get radius coordinates
    R = np.sqrt((RA - ra)**2 + (DEC - dec)**2)

    # select pixels
    select = R < radius

    # get peak brightness within pixel radius
    peak = np.max(data[select])

    # get rms outside of pixel radius
    if rms_max_r is not None and rms_max_r is not None:
        rms_select = (R < rms_max_r) & (R > rms_min_r)
        rms = np.sqrt(np.mean(data[select]**2))
    else:
        rms = np.sqrt(np.mean(data[~select]**2))

    # get peak error
    peak_err = rms / np.sqrt(Npix_beam / 2.0)

    # get frequency of image
    freq = head["CRVAL3"]

    ## fit a 2D gaussian and get integrated and peak flux statistics ##
    # recenter R array by peak flux point and get thata T array
    peak_ind = np.argmax(data[select])
    peak_ra = RA[select][peak_ind]
    peak_dec = DEC[select][peak_ind]
    X = (RA - peak_ra)
    Y = (DEC - peak_dec)
    R = np.sqrt(X**2 + Y**2)
    X[np.where(np.isclose(X, 0.0))] = 1e-5
    T = np.arctan(Y / X)

    # use synthesized beam as data mask
    ecc = head["BMAJ"] / head["BMIN"]
    beam_theta = head["BPA"] * np.pi / 180 + np.pi/2
    EMAJ = R * np.sqrt(np.cos(T+beam_theta)**2 + ecc**2 * np.sin(T+beam_theta)**2)
    fit_mask = EMAJ < (head["BMAJ"] / 2 * gaussfit_mult)
    masked_data = data.copy()
    masked_data[~fit_mask] = 0.0

    # fit 2d gaussian
    gauss_init = mod.functional_models.Gaussian2D(peak, ra, dec, x_stddev=head["BMAJ"]/2, y_stddev=head["BMIN"]/2) 
    fitter = mod.fitting.LevMarLSQFitter()
    gauss_fit = fitter(gauss_init, RA[fit_mask], DEC[fit_mask], data[fit_mask])

    # get gaussian fit properties
    peak_gauss_flux = gauss_fit.amplitude.value
    P = np.array([X, Y]).T
    beam_theta -= np.pi/2
    Prot = P.dot(np.array([[np.cos(beam_theta), -np.sin(beam_theta)], [np.sin(beam_theta), np.cos(beam_theta)]]))
    gauss_cov = np.array([[gauss_fit.x_stddev.value**2, 0], [0, gauss_fit.y_stddev.value**2]])
    try:
        model_gauss = stats.multivariate_normal.pdf(Prot, mean=np.array([0,0]), cov=gauss_cov)
        model_gauss *= gauss_fit.amplitude.value / model_gauss.max()
        int_gauss_flux = np.nansum(model_gauss.ravel()) / Npix_beam
    except:
        int_gauss_flux = 0

    return peak, peak_err, rms, peak_gauss_flux, int_gauss_flux, freq

if __name__ == "__main__":

    # parse args
    args = a.parse_args()

    # sort files
    files = sorted(args.files)

    # get filename
    if args.outdir is None:
        args.outdir = os.path.dirname(os.path.commonprefix(files))
    output_fname = os.path.join(args.outdir, os.path.basename(os.path.splitext(os.path.commonprefix(files))[0] + args.ext))
    if os.path.exists(output_fname) and args.overwrite is False:
        raise IOError("file {} exists, not overwriting".format(output_fname))

    # iterate over files
    peak_flux = []
    peak_flux_err = []
    peak_gauss_flux = []
    int_gauss_flux = []
    freqs = []
    for i, fname in enumerate(files):
        output = source_extract(fname, **vars(args))
        peak_flux.append(output[0])
        peak_flux_err.append(output[2])
        peak_gauss_flux.append(output[3])
        int_gauss_flux.append(output[4])
        freqs.append(output[5])

    freqs = np.array(freqs)
    peak_flux = np.array(peak_flux)
    peak_flux_err = np.array(peak_flux_err)
    peak_gauss_flux = np.array(peak_gauss_flux)
    int_gauss_flux = np.array(int_gauss_flux)

    data = np.vstack([freqs/1e6, peak_flux, peak_flux_err, peak_gauss_flux, int_gauss_flux]).T

    # save spectrum
    print("...saving {}".format(output_fname))
    np.savetxt(output_fname, data, fmt="%8.5f", header="freqs (MHz)\t peak flux (Jy/beam)\t flux err\t peak gaussfit (Jy/beam)\t integ gaussfit (Jy)", delimiter='\t')





<file_sep>/aas/Estimating Delays/nn/networks/CNN_DS_BN_BC.py
# CNN_DS_BN_BC

import sys
import numpy as np
import tensorflow as tf
from CNN_DS_BN_C import CNN_DS_BN_C 

class CNN_DS_BN_BC(CNN_DS_BN_C):
    """ CNN_DS_BN_BC

    Downsampling binary classifier
    CNN: Convolutional Neural Network.
    DS: DownSampling. Each convolution is followed by a downsampling convolution
    BN: All non-linearalities have batch-normalization applied.
    BC: Binary Classification, this network classifies samples of having one of two labels.

    Each layer starts with a 1x3 convolution with 2**i filters (i = layer index) with stride of 1.
    This is fed into a 1x5 convolution with the same number of filters, but stride of 2 (50% downsample)

    Leaky_ReLU and batch normalization is applied after each convolution. Downsamples convolutions have dropout.
    Biases are added before activations.

    Output of last layer is fed to fully connected layer (with no activation)

    Cost function is softmax cross entropy

    """
    __doc__ += CNN_DS_BN_C.__doc__
    def __init__(self,
                 name,
                 num_downsamples,
                 log_dir = '../logs/',
                 dtype = tf.float32,
                 adam_initial_learning_rate = 0.0001,
                 verbose = True):
        
         CNN_DS_BN_C.__init__(self,
                              name=name,
                              num_downsamples=num_downsamples,
                              num_classes=2,
                              log_dir=log_dir,
                              dtype=dtype,
                              adam_initial_learning_rate=adam_initial_learning_rate,
                              verbose=verbose)
 
<file_sep>/aas/Estimating Delays/nn/network_trainers/CNN_BC_Trainer.py
# CNN_BC_Trainer

import sys, os
sys.path.insert(1, os.path.join(sys.path[0], '../modules'))
from NN_Trainer import NN_Trainer

import tensorflow as tf
from tensorflow.python.client import timeline

import numpy as np
from sklearn.metrics import confusion_matrix
import matplotlib.pyplot as plt
import io
import itertools

class CNN_BC_Trainer(NN_Trainer):

    def __init__(self,
                 network,
                 Data_Creator, # Class
                 num_epochs = 100,
                 batch_size = 32,
                 log_dir = '../logs/',
                 model_save_interval = 25,
                 pretrained_model_path = None,
                 metric_names = ['costs', 'accuracies'],
                 sample_keep_prob = 0.80,
                 conv_keep_prob = 0.9,
                 verbose = True,
                 class_names = ['pos', 'neg']):
    
        NN_Trainer.__init__(self,
                            network = network,
                            Data_Creator = Data_Creator,
                            num_epochs = num_epochs,
                            batch_size = batch_size,
                            log_dir = log_dir,
                            model_save_interval = model_save_interval,
                            pretrained_model_path = pretrained_model_path,
                            metric_names = metric_names,
                            verbose = verbose)
        

        self.sample_keep_prob = sample_keep_prob
        self.conv_keep_prob = conv_keep_prob
        self.class_names = class_names
        
    def train(self):
        
        self.save_params()

        costs = []
        accuracies = []

        tf.reset_default_graph()
        
        self._network.create_graph()
        saver = tf.train.Saver()

        with tf.Session() as session:
            
            if self.pretrained_model_path == None:
                session.run(tf.global_variables_initializer())
                
            else:
                saver.restore(session, self.pretrained_model_path)

            archive_loc = self.log_dir + self._network.name
            training_writer = tf.summary.FileWriter(archive_loc + '/training', session.graph)
            testing_writer = tf.summary.FileWriter(archive_loc + '/testing', session.graph)
            self.model_save_location = archive_loc + '/trained_model.ckpt'   
            
            
            self._msg = '\rtraining';self._vprint(self._msg)

            try:
                for epoch in range(self.num_epochs):
                    
                    

                    training_inputs, training_labels = self._train_batcher.get_data(); self._train_batcher.gen_data()
                    testing_inputs, testing_labels = self._test_batcher.get_data(); self._test_batcher.gen_data()  

                    training_labels = np.asarray(training_labels)
                    testing_labels = np.asarray(testing_labels)
                    
                    # if the division here has a remainde some values are just truncated
                    batch_size = self.batch_size
                    num_entries = training_inputs.shape[0]

                    for j in range(int(num_entries/batch_size)):
                        self._msg = '\repoch'
                        self._msg += '- {:5.0f}/{:5.0f}'.format(epoch + 1,self.num_epochs)
                        self._msg += ' - batch: {:4.0f}/{:4.0f}'.format(j + 1, int(num_entries/batch_size))
                        if epoch != 0:
                            self._msg += ' - (Training, Testing) - '.format(epoch)
                            self._msg += " costs: ({:0.4f}, {:0.4f})".format(training_cost, testing_cost)
                            self._msg += " accss: ({:2.2f}, {:2.2f})".format(training_acc, testing_acc)
                        self._vprint(self._msg); 

                        training_inputs_batch = training_inputs[j*batch_size:(j + 1)*batch_size].reshape(-1,1,1024,1)
                        training_labels_batch = training_labels[j*batch_size:(j + 1)*batch_size].reshape(-1,2)

                        feed_dict = {self._network.X: training_inputs_batch,
                                     self._network.labels: training_labels_batch,
                                     self._network.sample_keep_prob : self.sample_keep_prob,
                                     self._network.conv_keep_prob : self.conv_keep_prob,
                                     self._network.is_training : True}

                        session.run([self._network.optimizer], feed_dict = feed_dict) 
                            
                    train_feed_dict = {self._network.X: training_inputs.reshape(-1,1,1024,1),
                                       self._network.labels: training_labels.reshape(-1,2),
                                       self._network.sample_keep_prob : 1.,
                                       self._network.conv_keep_prob : 1.,
                                       self._network.is_training : False}

                    train_predicts =  session.run([self._network.predictions], train_feed_dict)

                    train_feed_dict[self._network.image_buf] = self.plt_confusion_matrix(training_labels.reshape(-1,2), train_predicts)

                    training_cost, training_acc, training_summary = session.run([self._network.cost,
                                                                                 self._network.accuracy,
                                                                                 self._network.summary],
                                                                                 feed_dict = train_feed_dict) 

                    training_writer.add_summary(training_summary, epoch)
                    training_writer.flush()  
                

                    test_feed_dict = {self._network.X: testing_inputs.reshape(-1,1,1024,1),
                                      self._network.labels: testing_labels.reshape(-1,2),
                                      self._network.sample_keep_prob : 1.,
                                      self._network.conv_keep_prob : 1.,
                                      self._network.is_training : False} 
                                      
                    test_predicts =  session.run([self._network.predictions], test_feed_dict)

                    test_feed_dict[self._network.image_buf] = self.plt_confusion_matrix(testing_labels.reshape(-1,2), test_predicts) 


                    testing_cost, testing_acc, testing_summary = session.run([self._network.cost,
                                                                              self._network.accuracy,
                                                                              self._network.summary],
                                                                              feed_dict = test_feed_dict)

                    testing_writer.add_summary(testing_summary, epoch)
                    testing_writer.flush()
                    
                    costs.append((training_cost, testing_cost))
                    accuracies.append((training_acc, testing_acc))
                    
                    if (epoch + 1) % self.model_save_interval == 0:
                        saver.save(session, self.model_save_location, epoch + 1)
                
                self.msg = ''
            except KeyboardInterrupt:
                self._msg = ' TRAINING INTERUPPTED' # this never prints I dont know why
                pass

            self._msg += '\rtraining ended'; self._vprint(self._msg)
            
            training_writer.close()
            testing_writer.close()


        session.close()
        self._msg += ' - session closed'; self._vprint(self._msg)
        self._msg = ''

        self._metrics = [costs, accuracies]
        self.save_metrics()



    def plt_confusion_matrix(self, labels, pred, normalize=False, title='Confusion matrix'):
        """
        Given one-hot encoded labels and preds, displays a confusion matrix.
        Arguments:
            `labels`:
                The ground truth one-hot encoded labels.
            `pred`:
                The one-hot encoded labels predicted by a model.
            `normalize`:
                If True, divides every column of the confusion matrix
                by its sum. This is helpful when, for instance, there are 1000
                'A' labels and 5 'B' labels. Normalizing this set would
                make the color coding more meaningful and informative.
        """
        labels = [label.argmax() for label in np.asarray(labels).reshape(-1,2)] # bc
        pred = [label.argmax() for label in np.asarray(pred).reshape(-1,2)] #bc

        classes = self.class_names

        cm = confusion_matrix(labels, pred)

        #if normalize:
        cm = cm.astype('float')*100 / cm.sum(axis=1)[:, np.newaxis]
        cm = np.nan_to_num(cm, copy=True)
        cm = cm.astype('int')

        fig, ax = plt.subplots(figsize = (5,5), dpi = 144)
        im = ax.imshow(cm, interpolation='nearest', aspect='auto', cmap=plt.cm.Oranges, vmin = 0, vmax = 100)
        ax.set_title(title)
        cbar = fig.colorbar(im)
        tick_marks = np.arange(len(classes))
        ax.set_xticks(tick_marks)
        ax.set_xticklabels(classes)
        ax.set_yticks(tick_marks)
        ax.set_yticklabels(classes)

        ax.set_ylabel('True label')
        ax.set_xlabel('Predicted label')
        for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):
            #s = '{:2.0}'.format(cm[i, j]) if cm[i,j] >= 1 else '.'
            ax.text(j, i, format(cm[i, j], 'd') if cm[i,j]!=0 else '.',
                    horizontalalignment="center", fontsize=15, verticalalignment='center', color= "black")

        # plt.show()
        buf = io.BytesIO()
        fig.savefig(buf, format='png', dpi = 144)
        plt.close(fig)
        buf.seek(0)

        return buf.getvalue()<file_sep>/aas/Estimating Delays/README.md
# Estimating Delays
## ML experiments for HERA RedCal

Can we use ML to aid redcal? We use tensorflow to find out.

See `wrap_unwrap_initial_exploration.ipynb` for a brief overview of the problem of computing cable delays from complex data.

## Prereqs

```
numpy
matplotlib
tensorflow
pyuvdata
hera_cal
```
- and all their prereqs

## Data

Data is built with NRAO IDR-2 miriad / firstcal files (is this description correct?).

 Specific files used cover one ten minute window:
`zen_data/zen.2458098.58037.xx.HH.uv` &
`zen_data/zen.2458098.58037.xx.HH.uv.abs.calfits`

TODO: Add description of how data is manipulated. For now see `estdel/nn/data_manipulation.py`

## Directory Contents

`estdel/` - package for estimating delays

`experiments/` - Experiments in training different networks, solving Ax = b, and doing other things.<file_sep>/aas/Estimating Delays/estdel/estdel/estdel.py
"""estdel - for estimating delays
 - <NAME> <EMAIL>

Estimate the overcall cable delay in visibility ratios.

For each list of 60 x 1024 complex visibility ratios produce 60 estimated cable delays.

Passes each row through one of (or both of) two trained neural networks. 

Sign network classifies the sign of the delay as being positive or negative.
Magnitude network classfies the magnitude of the delay as being in one of 401 classes.
    -   401 classes from 0.00 to 0.0400, each of width 0.0001
    -   default conversion to 401 classes from 0 to ~ 400 ns

VratioDelaySign estimates the sign of the delay (< 0 or >= 0)
VratioDelayMagnitude estimates the magnitude of the delay
VratioDelay provides Delay_Sign * Delay_Magnitude

tau = # slope in range -0.0400 to 0.0400
nu = np.arange(1024) # unitless frequency channels
phi = # phase in range 0 to 1
data = np.exp(-2j*np.pi*(nu*tau + phi))

estimator = estdel.VratioDelay(data)
prediction = estimator.predict()

# prediction should output tau
"""

import pkg_resources
import numpy as np
import tensorflow as tf

# suppress tensorflow INFO messages
tf.logging.set_verbosity(tf.logging.WARN)

_TRAINED_MODELS_DIR = 'trained_models'

# fn for best postive-negative classifier
_SIGN_PATH = 'sign_NN_frozen.pb'

# fn for best magnitude classifier
_MAG_PATH = 'mag_NN_frozen.pb'

# number of frequency channels
N_FREQS = 1024



class _DelayPredict(object):
    """_DelayPredict
    
    Base class for predictions.
    
    Data processing and predictions.
    
    """
    
    def __init__(self, data):
        """__init__
        
        Constants:
            _n_freqs = 1024
        
        Args:
            data (list of complex floats): Visibliity data
        """
        self._data = data
        self._n_freqs = N_FREQS

        assert(self._data.shape[-1] == self._n_freqs)


    def _angle_tx(self, x):
        """_angle_tx
        
        Scales (-pi, pi) to (0, 1)
        
        Args:
            x (numpy array of floats): Angle data in range -pi to pi
        
        Returns:
            numpy array of floats: scaled angle data
        """
        
        return (x + np.pi) / (2. * np.pi)
    

    def _preprocess_data(self):
        """_preprocess_data
        
            Converts data from complex to real (angle), scales, and reshapes
        
        Returns:
            numpy array of floats: Scaled and reshaped angle data
        """
        return self._angle_tx(np.angle(self._data)).reshape(-1,1,self._n_freqs,1)


    def _predict(self):
        """_predict
        
        Import frozen tensorflow network, activate graph, 
        feed data, and make prediction.
        
        """

        resource_package = __name__ 
        resource_path = '/'.join((_TRAINED_MODELS_DIR, self._model_path))
        path = pkg_resources.resource_filename(resource_package, resource_path)

        with tf.gfile.GFile(path, "rb") as f:
            restored_graph_def = tf.GraphDef()
            restored_graph_def.ParseFromString(f.read())

        with tf.Graph().as_default() as graph:
            tf.import_graph_def(restored_graph_def, input_map=None, return_elements=None, name="")

            sample_keep_prob = graph.get_tensor_by_name('keep_probs/sample_keep_prob:0')
            conv_keep_prob = graph.get_tensor_by_name('keep_probs/conv_keep_prob:0')
            is_training = graph.get_tensor_by_name('is_training:0')
            X = graph.get_tensor_by_name('sample/X:0')

            # add hook to output operation
            pred_cls = graph.get_tensor_by_name('predictions/ArgMax:0')

        with tf.Session(graph=graph) as sess:
            feed_dict = {sample_keep_prob : 1.,
                         conv_keep_prob : 1.,
                         is_training : False,
                         X: self.data}

            # collect prediction
            self._pred_cls = sess.run(pred_cls, feed_dict = feed_dict)

            sess.close()




class VratioDelaySign(_DelayPredict):
    """VratioDelaySign
    
    Estimates visibility ratio cable delay sign by using two pretrained neural networks.
    
    Methods:
        predict()
            - call to make sign prediction for data
    
    Attributes:
        data (numpy array of floats): Input data of redundant visibility ratios is processed for predictions
        predictions (numpy array of floats): The converted raw magnitude predictions (see predict())
        raw_predictions (list of floats): The raw magnitude predictions from the network
    """

    def __init__(self, data):
        """__init__
        
        Preprocesses data for prediction.
        
            - converts complex data to angle
            - scales angles to range preferred by networks
            - reshapes 2D data to 4D tensor
        
        
        Args:
            data (list of complex): shape = (N, 1024)
                - redundant visibility ratios
        
        """
        _DelayPredict.__init__(self, data = data)

        self.data = self._preprocess_data()
        self._model_path = _SIGN_PATH


    def _pred_cls_to_sign(self):
        """_pred_cls_to_sign
        
            Convert index of predicted class index to value
              1 if class is postive
             -1 if class is negative
        
        Returns:
            list of ints: 
        """
        return [1 if x == 0 else -1 for x in self._pred_cls]
         
    

    def predict(self):
        """predict
        
        Returns:
            numpy array of floats: sign predictions
        """
        self._predict()
        self.raw_predictions = self._pred_cls_to_sign()
        self.predictions = np.array(self.raw_predictions)

        return self.predictions


def _default_conversion_fn(x):
    """_default_conversion_fn
    
    Convert unitless predictions to nanoseconds
    Based on 1024 channels and 
    0.100 GHz - 0.200 GHz range
    
    Args:
        x (numpy array of floats): Predicted values in range 0.000 to 0.040
    
    Returns:
        numpy array of floats: Converted predicted value
    """
    
    freqs = np.linspace(0.100,0.200,N_FREQS) 
    channel_width_in_GHz = np.mean(np.diff(freqs))

    return np.array(x) / channel_width_in_GHz




class VratioDelayMagnitude(_DelayPredict):
    """VratioDelayMagnitude
    
    Estimates visibility ratio total cable delay by using two pretrained neural networks.
    
    Methods:
        predict()
            - call to make prediction
    
    Attributes:
        data (list of complex floats): Visibility ratios
        predictions (numpy array of floats): The converted raw magnitude predictions (see predict())
        raw_predictions (list of floats): The raw magnitude predictions from the network
    """

    def __init__(self, data):
        """Preprocesses data for prediction.
        
            - converts complex data to angle
            - scales angles to range preferred by networks
            - reshapes 2D data to 4D tensor
        
        
        Args:
            data (list of complex floats): shape = (N, 1024)
                - redundant visibility ratios
        
        """
        _DelayPredict.__init__(self, data = data)

        self.data = self._preprocess_data()
        self._model_path = _MAG_PATH


    def _pred_cls_to_magnitude(self):
        """_pred_cls_to_magnitude
        
        Convert predicted label index to magnitude
        
        Returns:
            list of floats: Magnitudes
        """
        magnitudes = np.arange(0,0.04 + 0.0001, 0.0001)
        return [magnitudes[x] for x in self._pred_cls]
    

    def predict(self, conversion_fn='default'):
        """predict
        
        Args:
            conversion_fn (None, str, or function): - None - Do no conversion, output predictions are the raw predictions
                - 'default' - convert raw predictions to ns by using frequencies 
                    with a 100MHz range over 1024 channels
                - OR provide your own function to do the conversion
                    - one required argument, the raw predictions
        
        Returns:
            numpy array of floats or list of floats: predictions
        """
        self._conversion_fn = conversion_fn
        self._predict()
        self.raw_predictions = self._pred_cls_to_magnitude()

        if self._conversion_fn is None:
            
            self.predictions = self.raw_predictions

        if self._conversion_fn == 'default':

            self.predictions = _default_conversion_fn(self.raw_predictions)

        else:

            self.predictions = self._conversion_fn(self.raw_predictions)

        return np.array(self.predictions)


class VratioDelay(object):
    """VratioDelay
    
    Estimates visibility ratio total cable delay by using two pretrained neural networks.
    
    Methods:
        predict()
            - call to make prediction
    
    Attributes:
        raw_predictions (list of floats): The raw predictions from the network
        predictions (numpy array of floats or list of floats) = The converted raw predictions

    """
    
    def __init__(self, data):
        """__init__
        
        Preprocesses data for prediction.
        
            - converts complex data to angle
            - scales angles to range preferred by networks
            - reshapes 2D data to 4D tensor
        
        
        Args:
            data (list of complex floats): shape = (N, 1024)
                - redundant visibility ratios
        
        """
        
        self._mag_evaluator = VratioDelayMagnitude(data)
        self._sign_evaluator = VratioDelaySign(data)
    

    def predict(self, conversion_fn='default'):
        """predict
        
            Make predictions
        
        Args:
            conversion_fn (str, optional): Description
            conversion_fn (None, 'default', or function
                - None - Do no conversion, output predictions are the raw predictions
                - 'default' - convert raw predictions to ns by using frequencies 
                    with a 100MHz range over 1024 channels
                - OR provide your own function to do the conversion
                    - takes in one argument, the raw predictions
        
        Returns:
            numpy array of floats or list of floats: Predicted values
        """
        signs = self._sign_evaluator.predict()
        mags = self._mag_evaluator.predict(conversion_fn=conversion_fn)
        self.raw_predictions = [self._mag_evaluator.raw_predictions[i]*signs[i] for i in range(len(signs))]
        self.predictions = signs*mags
        return self.predictions




class DelaySolver(object):
    """DelaySolver
    
    Args:
        list_o_sep_pairs (list of lists of tuples): shape = (N, 2, 2)
            ex: list_o_sep_pairs[0] = [(1, 2), (3, 4)]
            each length 2 sublist is made of two redundant separations, one in each tuple
    
        v_ratios (list of complex floats): shape = (N, 60, 1024)
    
            Complex visibility ratios made from the the corresponding redundant sep pairs in list_o_sep_pairs
    
        true_ant_delays (dict): dict of delays with antennas as keys,
            ex : true_ant_delays[143] = 1.2
            if conversion_fn == 'default', ant delays should be in ns
    
    Attributes:
        A (numpy array of ints): The matrix representing the redundant visibility ratios
        b (numpy array of floats): A times x
        unique_ants (numpy array of ints): All the unique antennas in list_o_sep_pairs
        v_ratio_row_predictions (numpy array of floats or list of floats): Predicted values
        v_ratio_row_predictions_raw (list of floats): Predicted values with no conversion
        x (list floats): True delays in order of antenna
    """
    
    def __init__(self,
                 list_o_sep_pairs,
                 v_ratios,
                 conversion_fn='default',
                 true_ant_delays={}, # dict {ant_idx : delay}
                ):
        """__init__

        Preprocess data, make predictions, covert data to ns, 
        construct A matrix.

        """
        self._list_o_sep_pairs = list_o_sep_pairs # shape = (N, 2, 2)
        self._v_ratios = v_ratios # complex, shape = (N, 60, 1024) # will be reshaped to (-1, 1, 1024, 1)
        self._true_ant_delays = true_ant_delays
        self._conversion_fn = conversion_fn

        self.unique_ants = np.unique(list_o_sep_pairs)
        self._max_ant_idx = np.max(self.unique_ants) + 1 
                                    
        self._make_A_from_list_o_sep_pairs()

        self._predictor = VratioDelay(v_ratios)
        self._predictor.predict(conversion_fn=conversion_fn)

        self.v_ratio_row_predictions = self._predictor.predictions
        self.v_ratio_row_predictions_raw = self._predictor.raw_predictions
        

    def _get_A_row(self, sep_pair):
        """_get_A_row

        constructs a single row of A from a sep pair

        Args:
            sep_pair (tuple of ints): Antennas for this visibility

        Returns;
            list of ints: a row of A
    
        """
        
        a = np.array(list(sum(sep_pair, () ))) # [(1,2) , (3, 4)] -> [1, 2, 3, 4]

        # construct the row
        # https://stackoverflow.com/a/29831596

        # row is 4 x _max_ant_idx, all zeros
        row = np.zeros((a.size, self._max_ant_idx), dtype = int)

        # for each element in sep_pair, got to the corresponding row
        # and assign the corresponding antenna the value 1
        row[np.arange(a.size),a] = 1 

        # flip the sign of the middle two rows
        row[1] *= -1
        row[2] *= -1

        # add the rows, row is now 1 x _max_ant_idx
        row = np.sum(row, axis = 0)

        return row


    def _make_A_from_list_o_sep_pairs(self):
        """_make_A_from_list_o_sep_pairs

        Construct A row by row
        """

        self.A = []
        for sep_pair in self._list_o_sep_pairs:

            # each visibility ratio of height 60 has one sep_pair
            # so make 60 idential rows in A for each visibility
            # so that A is the correct shape
            # (because the prediction will output a unique prediction 
            # for each row in the visibility ratio)
            self.A.append(np.tile(self._get_A_row(sep_pair), (60,1)))

        self.A =  np.asarray(self.A).reshape(-1, self._max_ant_idx)


    def true_b(self):
        """ true_b

        do A times x to find the true values for each visibility
        where x is a list of the true antenna delays in order of antenna

        Returns:
            numpy array of floats
        """

        
        self.x = [self._true_ant_delays[ant] for ant in self.unique_ants]
        self.b = np.matmul(self.A[:, self.unique_ants], self.x)

        return self.b


<file_sep>/aas/Estimating Delays/experiments/README.md
# experiments

The notebooks `CNN_*` and `FNN_*` are where I try training netowrks with different hyperparameters to see how they can learn the data.

`Classifier_Class_Explanation` is where I explain of how I convert this regression problem into a classification problem.

`Classification_Evaluator` is where I try see how well trained classifcation networks perform.

`Solver` is where I try to solve Ax=b<file_sep>/nsk/scripts/skynpz2calfits.py
#!/usr/bin/env python2.7
"""
skynpz2calfits.py
---------------

convert sky_image.py calibration
solution output in .npz format
(originally from CASA .cal/ tables)
into pyuvdata .calfits format

<NAME>
Dec. 2017
"""
import matplotlib
matplotlib.use('agg')
import matplotlib.pyplot as plt
from pyuvdata import UVCal, UVData
import numpy as np
import argparse
import os
import scipy.signal as signal
from sklearn import gaussian_process as gp
import copy
import hera_cal as hc

a = argparse.ArgumentParser(description="Turn CASA calibration solutions in {}.npz files from sky_image.py script into .calfits files")

# Required Parameters
a.add_argument("--fname", type=str, help="output calfits filename.", required=True)
a.add_argument("--uv_file", type=str, help="Path to original miriad uv file of data.", required=True)
# Delay Solution Parameters
a.add_argument("--dly_files", type=str, default=None, nargs='*', help="Path to .npz file(s) with antenna delay output from sky_image.py (CASA K cal)")
a.add_argument("--TTonly", default=False, action="store_true", help="only store Tip-Tilt slopes of delay solution.")
a.add_argument("--plot_dlys", default=False, action="store_true", help="plot delay solutions across array.")
# Phase Solution Parameters
a.add_argument("--phs_files", type=str, default=None, nargs='*', help="Path to .npz file(s) with phase output from sky_image.py (CASA G phscal)")
a.add_argument("--plot_phs", default=False, action="store_true", help="plot phase solutions across array.")
# Amplitude Solution Parameters
a.add_argument("--amp_files", type=str, default=None, nargs='*',  help="Path to .npz file(s) with amplitude output from sky_image.py (CASA G ampcal)")
a.add_argument("--plot_amp", default=False, action='store_true', help='plot amp solution across array.')
a.add_argument("--gain_amp_antavg", default=False, action='store_true', help="average gain amplitudes across antennas")
# Bandpass Solution Parameters
a.add_argument("--bp_files", type=str, default=None, nargs='*', help="Path to .npz file(s) with antenna complex bandpass output from sky_image.py (CASA bandpass)")
a.add_argument("--bp_flag_frac", type=float, default=0.5, help="at each freq bin, fraction of antennas flagged needed to broadcast flag to all ants.")
a.add_argument("--bp_broad_flags", default=False, action='store_true', help="broadcast flags at freq channels that satisfy bp_flag_frac")
a.add_argument("--noBPamp", default=False, action='store_true', help="set BP amplitude solutions to zero.")
a.add_argument("--noBPphase", default=False, action='store_true', help="set BP phase solutions to zero.")
a.add_argument('--bp_medfilt', default=False, action='store_true', help="median filter bandpass solutions")
a.add_argument('--medfilt_flag', default=False, action='store_true', help='use median filter to flag bad BP solutions in frequency')
a.add_argument('--medfilt_kernel', default=7, type=int, help="kernel size (channels) for BP median filter")
a.add_argument('--bp_amp_antavg', default=False, action='store_true', help="average bandpass amplitudes across antennas")
a.add_argument('--bp_TTonly', default=False, action='store_true', help="use only tip-tilt phase mode in bandpass solution")
a.add_argument('--plot_bp', default=False, action='store_true', help="plot final bandpass solutions")
a.add_argument('--bp_gp_smooth', default=False, action='store_true', help='smooth bandpass w/ gaussian process. Recommended to precede w/ bp_medfilt.')
a.add_argument('--bp_gp_max_dly', default=200.0, type=float, help="maximum delay in nanosec allowed for gaussian process fit")
a.add_argument('--bp_gp_nrestart', default=1, type=int, help='number of restarts for GP hyperparameter gradient descent.')
a.add_argument('--bp_gp_thin', default=2, type=int, help="thinning factor for freq bins in GP smooth fit")
# Misc
a.add_argument("--out_dir", default=None, type=str, help="output directory for calfits file. Default is working directory path")
a.add_argument("--overwrite", default=False, action="store_true", help="overwrite output calfits file if it exists")
a.add_argument("--multiply_gains", default=False, action="store_true", help="change gain_convention from divide to multiply.")
a.add_argument('--silence', default=False, action='store_true', help="silence output to stdout")

def echo(message, type=0, verbose=True):
    if verbose:
        if type == 0:
            print(message)
        elif type == 1:
            print('\n{}\n{}'.format(message, '-'*40))

def skynpz2calfits(fname, uv_file, dly_files=None, amp_files=None, bp_files=None, out_dir=None, phs_files=None, overwrite=False,
                   TTonly=True, plot_dlys=False, plot_phs=False, gain_convention='multiply', plot_bp=False, noBPphase=False,
                   noBPamp=False, bp_TTonly=False, bp_medfilt=False, medfilt_kernel=5, bp_amp_antavg=False, bp_flag_frac=0.3,
                   bp_broad_flags=False, medfilt_flag=False, bp_gp_smooth=False, bp_gp_max_dly=500.0, bp_gp_nrestart=1, bp_gp_thin=2,
                   plot_amp=False, gain_amp_antavg=False, verbose=True):
    """
    Convert *.npz output from sky_image.py into single-pol calfits file.
    uv_file must be single-pol.

    """
    # get out_dir
    if out_dir is None:
        out_dir = "./"

    # load UVData
    echo("...loading uv_file", verbose=verbose)
    uvd = UVData()
    uvd.read_miriad(uv_file)

    # get ants and antpos
    antpos, ants = uvd.get_ENU_antpos(center=True, pick_data_ants=True)
    ants = list(ants)

    # get freqs, times, jones
    freqs = np.unique(uvd.freq_array)
    times = np.unique(uvd.time_array)
    Nfreqs = len(freqs)
    Nants = len(ants)
    jones = uvd.polarization_array
    num2str = {-5: 'x', -6: 'y', -7: 'jxy', -8: 'jyx'}
    str2num = {'x': -5, 'y': -6, 'jxy': -7, 'jyx': -8}
    pols = np.array(map(lambda j: num2str[j], jones))

    # construct blank gains and flags
    gains = np.ones((Nants, Nfreqs, 1, 1), dtype=np.complex)
    flags = np.zeros((Nants, Nfreqs, 1, 1), dtype=np.bool)
    flagged_ants = np.zeros(Nants, dtype=np.bool)

    # process delays if available
    if dly_files is not None:
        echo("...processing delays", verbose=verbose, type=1)
        delays = np.zeros_like(ants, dtype=np.float)
        delay_flags = np.zeros_like(ants, dtype=np.bool)
        for dly_file in dly_files:
            echo("...processing {}".format(dly_file))
            # get CASA delays and antennas
            dly_data = np.load(dly_file)
            dly_ants = dly_data['delay_ants']
            dlys = dly_data['delays']
            dly_flags = dly_data['delay_flags']
            dly_ants = dly_ants.tolist()
            dlys *= 1e-9
            # reorder antennas to be consistant w/ ants array
            dlys = np.array(map(lambda a: dlys[dly_ants.index(a)] if a in dly_ants else 0.0, ants))
            dly_flags = np.array(map(lambda a: dly_flags[dly_ants.index(a)] if a in dly_ants else True, ants))
            # keep only limited information
            if TTonly:
                echo("...keeping only TT component of delays", verbose=verbose)
                A = np.vstack([antpos[:, 0], antpos[:, 1]]).T
                fit = np.linalg.pinv(A.T.dot(A)).dot(A.T).dot(dlys)
                dlys = A.dot(fit)
            # add to delay array
            delays += dlys
            delay_flags += dly_flags
            flagged_ants += dly_flags

        # turn into complex gains
        dly_gains = np.exp(2j*np.pi*(freqs-freqs.min())*delays.reshape(-1, 1))[:, :, np.newaxis, np.newaxis]
        dly_gain_flags = delay_flags[:, np.newaxis, np.newaxis, np.newaxis]

        # multiply into gains
        gains *= dly_gains

        # add into flags
        flags += dly_gain_flags

    # process overall phase if available
    if phs_files is not None:
        echo("...processing gain phase", verbose=verbose, type=1)
        phases = np.zeros_like(ants, dtype=np.float)
        phase_flags = np.zeros_like(ants, dtype=np.bool)
        for phs_file in phs_files:
            echo("...processing {}".format(phs_file))
            # get phase antennas and phases
            phs_data = np.load(phs_file)
            phs_ants = phs_data['phase_ants']
            phs = phs_data['phases']
            phs_flags = phs_data['phase_flags']
            phs_ants = phs_ants.tolist()
            # reorder to match ants
            phs = np.array([phs[phs_ants.index(a)] if a in phs_ants else 0.0 for a in ants])
            phs_flags = np.array(map(lambda a: phs_flags[phs_ants.index(a)] if a in phs_ants else True, ants))
            # add to phases
            phases += phs
            phase_flags += phs_flags
            flagged_ants += phs_flags

        # construct gains
        phase_gains = np.exp(1j * phases)
        phase_gains = phase_gains[:, np.newaxis, np.newaxis, np.newaxis]
        phase_gain_flags = phase_flags[:, np.newaxis, np.newaxis, np.newaxis]

        # mult into gains
        gains *= phase_gains

        # add into flags
        flags += phase_gain_flags

    # process overall amplitude if available
    if amp_files is not None:
        echo("...processing gain amp", verbose=verbose, type=1)
        amplitudes = np.ones_like(ants, dtype=np.float)
        amplitude_flags = np.zeros_like(ants, dtype=np.bool)
        for amp_file in amp_files:
            echo("...processing {}".format(amp_file))
            # get amp antenna and amps
            amp_data = np.load(amp_file)
            amp_ants = amp_data['amp_ants']
            amps = amp_data['amps']
            amp_flags = amp_data['amp_flags']
            amp_ants = amp_ants.tolist()
            # reorder to match ants
            amps = np.array([amps[amp_ants.index(a)] if a in amp_ants else 1.0 for a in ants])
            amp_flags = np.array(map(lambda a: amp_flags[amp_ants.index(a)] if a in amp_ants else True, ants))
            # add to amplitudes
            amplitudes *= amps
            amplitude_flags += amp_flags
            flagged_ants += amp_flags

        # average across ants if desired
        if gain_amp_antavg:
            echo("...averaging antenna gain amplitudes", verbose=verbose)
            avg_amp = np.median(amplitudes[~amplitude_flags])
            amplitudes *= avg_amp / amplitudes

        # construct gains
        amplitude_gains = amplitudes
        amplitude_gains = amplitude_gains[:, np.newaxis, np.newaxis, np.newaxis]
        amplitude_flags = amplitude_flags[:, np.newaxis, np.newaxis, np.newaxis]

        # mult into gains
        gains *= amplitude_gains

        # add into flags
        flags += amplitude_flags

    # process bandpass if available
    if bp_files is not None:
        echo("...processing bandpass", verbose=verbose, type=1)
        for ii, bp_file in enumerate(bp_files):
            echo("...processing {}".format(bp_file))
            # get bandpass and form complex gains
            bp_data = np.load(bp_file)
            bp_freqs = bp_data['bp_freqs']
            bp_Nfreqs = len(bp_freqs)
            bandp_gains = bp_data['bp']
            bandp_flags = bp_data['bp_flags']
            bp_ants = bp_data['bp_ants'].tolist()
            # reorder to match ants
            bandp_gains = np.array([bandp_gains[:, bp_ants.index(a)].squeeze() if a in bp_ants else np.ones((bp_Nfreqs), dtype=np.complex) for a in ants])
            bandp_flags = np.array([bandp_flags[:, bp_ants.index(a)].squeeze() if a in bp_ants else np.ones((bp_Nfreqs), dtype=np.bool) for a in ants])
            # broadcast flags to all antennas at freq channels that satisfy bp_flag_frac
            flag_broadcast = (np.sum(bandp_flags, axis=0) / float(bandp_flags.shape[0])) > bp_flag_frac
            if bp_broad_flags:
                bandp_flags += np.repeat(flag_broadcast.reshape(1, -1).astype(np.bool), Nants, 0)
            # configure gains and flags shapes
            bandp_gains = bandp_gains[:, :, np.newaxis, np.newaxis]
            bandp_flags = bandp_flags[:, :, np.newaxis, np.newaxis]
            if ii == 0:
                bp_gains = bandp_gains
                bp_flags = bandp_flags
            else:
                bp_gains *= bandp_gains
                bp_flags += bandp_flags

        # median filter if desired
        if bp_medfilt or medfilt_flag:
            echo("...median filtering", verbose=verbose)
            bp_gains_medfilt = signal.medfilt(bp_gains.real, kernel_size=(1, medfilt_kernel, 1, 1)) + 1j*signal.medfilt(bp_gains.imag, kernel_size=(1, medfilt_kernel, 1, 1))

            if medfilt_flag:
                echo("...solving for BP flags w/ medfilt data")
                # get residual and MAD from unfiltered and filtered data
                residual = (np.abs(bp_gains) - np.abs(bp_gains_medfilt))
                residual[bp_flags] *= np.nan
                residual = residual.squeeze()
                resid_std = np.nanmedian(np.abs(residual - np.nanmedian(residual, axis=1)[:, np.newaxis]), axis=1) * 1.5
                # identify outliers as greater than 10-sigma
                bad = np.array([np.abs(residual[i]) > resid_std[i]*10 for i in range(residual.shape[0])])
                bp_flags += bad[:, :, np.newaxis, np.newaxis]

            if bp_medfilt:
                echo("...bandpass is the median filtered bandpass")
                bp_gains = bp_gains_medfilt

        # average amplitude across antenna if desired
        if bp_amp_antavg:
            echo("...averaging bandpass amplitude across antennas", verbose=verbose)
            amp_avg = np.nanmedian(np.abs(bp_gains), axis=0).reshape(1, -1, 1, 1)
            bp_gains *= amp_avg / np.abs(bp_gains)

        # smooth bandpass w/ gaussian process if desired
        if bp_gp_smooth:
            echo("...smoothing with gaussian process", verbose=verbose)
            freq_lambda = 1. / (bp_gp_max_dly*1e-3) # MHz
            kernel = 1**2 * gp.kernels.RBF(freq_lambda + 10, (freq_lambda, 200.0)) + gp.kernels.WhiteKernel(1e-4, (1e-8, 1e0))

            # configure data
            X = bp_freqs / 1e6
            bp_gains_real = []
            bp_gains_imag = []
            # iterate over antennas
            for i, a in enumerate(ants):
                if i % 10 == 0: echo("{}/{} ants".format(i, len(ants)), verbose=verbose)
                # skip flagged ants
                if np.min(bp_flags[i, :]):
                    bp_gains_real.append(np.ones_like(X.squeeze()))
                    bp_gains_imag.append(np.zeros_like(X.squeeze()))
                    continue
                GP = gp.GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=bp_gp_nrestart)
                xdata = X[~bp_flags[i].squeeze(), :][::bp_gp_thin]
                yreal = bp_gains[i].squeeze()[~bp_flags[i].squeeze()][::bp_gp_thin].real
                yimag = bp_gains[i].squeeze()[~bp_flags[i].squeeze()][::bp_gp_thin].imag
                ydata = np.vstack([yreal, yimag]).T
                ydata_med = np.median(ydata, axis=0)
                ydata -= ydata_med
                GP.fit(xdata, ydata)
                ypred = GP.predict(X) + ydata_med
                bp_gains_real.append(ypred[:, 0])
                bp_gains_imag.append(ypred[:, 1])
            # reconstruct bp gains
            bp_gains = (np.array(bp_gains_real) + 1j*np.array(bp_gains_imag))[:, :, np.newaxis, np.newaxis]

        # take only tip-tilt if desired
        if bp_TTonly:
            echo("...distilling bandpass to only tip-tilt phase", verbose=verbose)
            # get bandpass phase array
            bp_phs = np.angle(bp_gains).reshape(Nants, bp_Nfreqs)
            # form least squares estimate of phase slopes along X and Y
            A = np.vstack([antpos[:, 0], antpos[:, 1]]).T
            fit = np.linalg.pinv(A.T.dot(A)).dot(A.T).dot(bp_phs)
            # median filter
            fit = signal.medfilt(fit, kernel_size=(1, 11))
            # smooth across frequency via FFT
            window = signal.windows.gaussian(fit.shape[1], fit.shape[1]/50.0)
            fit_smooth = signal.fftconvolve(fit, window.reshape(1, -1), mode='same') / np.sum(window)
            # make prediction to get tip-tilt estimates
            bp_TT_phs = A.dot(fit_smooth).reshape(Nants, bp_Nfreqs, 1, 1)
            # update bandpass gains
            bp_gains *= np.exp(1j*bp_TT_phs - 1j*np.angle(bp_gains))

        # suppress amp and/or phase if desired
        if noBPamp:
            echo("...eliminating BP amplitudes", verbose=verbose)
            bp_gains /= np.abs(bp_gains)

        if noBPphase:
            echo("...eliminating BP phases", verbose=verbose)
            bp_gains /= np.exp(1j*np.angle(bp_gains))

        # mult into gains
        bp_freq_select = np.array([np.argmin(np.abs(freqs-x)) for x in bp_freqs])
        gains[:, bp_freq_select, :, :] *= bp_gains
        flags[:, bp_freq_select, :, :] += bp_flags
        bp_flagged_ants = np.min(bp_flags, axis=1).squeeze()
        flagged_ants += bp_flagged_ants

    # check filename
    if fname.split('.')[-1] != "calfits":
        fname += ".calfits"

    # make into calfits
    gain_dict = {}
    flag_dict = {}
    for i, a in enumerate(ants):
        gain_dict[(a, pols[0])] = gains[i, :, :, 0].T
        flag_dict[(a, pols[0])] = flags[i, :, :, 0].T
        if flagged_ants[i]:
            flag_dict[(a, pols[0])] += True
    uvc = hc.io.write_cal(fname, gain_dict, freqs, times[:1], flags=flag_dict, outdir=out_dir,
                          overwrite=overwrite, gain_convention=gain_convention)

    # plot dlys
    if plot_dlys:
        fig, ax = plt.subplots(1, figsize=(8,6))
        ax.grid(True)
        dly_max = np.max(np.abs(delays*1e9))
        dly_min = -dly_max
        for i, a in enumerate(ants):
            if flagged_ants[i] == True:
                continue
            cax = ax.scatter(antpos[i, 0], antpos[i, 1], c=delays[i]*1e9, s=200, cmap='coolwarm', vmin=dly_min, vmax=dly_max)
            ax.text(antpos[i,0]+1, antpos[i,1]+2, str(a), fontsize=12)
        cbar = fig.colorbar(cax)
        cbar.set_label("delay [nanosec]", size=16)
        ax.set_xlabel("X [meters]", fontsize=14)
        ax.set_ylabel("Y [meters]", fontsize=14)
        ax.set_title("delay solutions for {}".format(os.path.basename(dly_files[0])), fontsize=10)
        fig.savefig(dly_files[0]+'.png', dpi=100, bbox_inches='tight', pad=0.05)
        plt.close()

    # plot phs
    if plot_phs:
        fig, ax = plt.subplots(1, figsize=(8,6))
        ax.grid(True)
        phs_max = np.pi
        phs_min = -np.pi
        for i, a in enumerate(ants):
            if flagged_ants[i] == True:
                continue
            cax = ax.scatter(antpos[i, 0], antpos[i, 1], c=phases[i], s=200, cmap='viridis', vmin=phs_min, vmax=phs_max)
            ax.text(antpos[i,0]+1, antpos[i,1]+2, str(a), fontsize=12) 
        cbar = fig.colorbar(cax)
        cbar.set_label("phase [radians]", size=16)
        ax.set_xlabel("X [meters]", fontsize=14)
        ax.set_ylabel("Y [meters]", fontsize=14)
        ax.set_title("phase solutions for {}".format(os.path.basename(phs_files[0])), fontsize=10)
        fig.savefig(phs_files[0]+'.png', dpi=100, bbox_inches='tight', pad=0.05)
        plt.close()

    # plot amp
    if plot_amp:
        fig, ax = plt.subplots(1, figsize=(8,6))
        ax.grid(True)
        amp_med = np.nanmedian(amplitudes[~amplitude_flags.squeeze()])
        amp_std = np.std(amplitudes[~amplitude_flags.squeeze()])
        amp_max = amp_med + amp_std * 2
        amp_min = amp_med - amp_std * 2
        for i, a in enumerate(ants):
            if flagged_ants[i] == True:
                continue
            cax = ax.scatter(antpos[i, 0], antpos[i, 1], c=amplitudes[i], s=200, cmap='rainbow', vmin=amp_min, vmax=amp_max)
            ax.text(antpos[i,0]+1, antpos[i,1]+2, str(a))
        cbar = fig.colorbar(cax)
        cbar.set_label("amplitude", size=16)
        ax.set_xlabel("X [meters]", fontsize=14)
        ax.set_ylabel("Y [meters]", fontsize=14)
        ax.set_title("amplitude solutions for {}".format(os.path.basename(amp_files[0])), fontsize=10)
        fig.savefig(amp_files[0]+'.png', dpi=100, bbox_inches='tight', pad=0.05)
        plt.close()

    # plot bandpass
    if plot_bp:
        fig, axes = plt.subplots(2, 1, figsize=(12,6))
        fig.subplots_adjust(hspace=0.3)
        # amplitude
        ax = axes[0]
        ax.grid(True)
        bp_gains[bp_flags] *= np.nan
        pls = []
        bp_ant_select = []
        ant_sort = np.argsort(ants)
        for i, a in enumerate(np.array(ants)[ant_sort]):
            if flagged_ants[ant_sort][i] == True:
                continue
            p, = ax.plot(bp_freqs / 1e6, np.abs(bp_gains).squeeze()[ants.index(a)], marker='.', ls='')
            pls.append(p)
        ax.set_xlabel("Frequency [MHz]", fontsize=12)
        ax.set_ylabel("Amplitude", fontsize=12)
        ax.set_title("bandpass for {}".format(bp_file), fontsize=10)
        # phase
        ax = axes[1]
        ax.grid(True)
        plot_ants = []
        for i, a in enumerate(np.array(ants)[ant_sort]):
            if flagged_ants[ant_sort][i] == True:
                continue
            plot_ants.append(a)
            ax.plot(bp_freqs / 1e6, np.angle(bp_gains).squeeze()[ants.index(a)], marker='.', ls='')
        ax.set_xlabel("Frequency [MHz]", fontsize=12)
        ax.set_ylabel("Phase [radians]", fontsize=12)
        lax = fig.add_axes([1.01, 0.1, 0.05, 0.8])
        lax.axis('off')
        lax.legend(pls, plot_ants, ncol=2)
        fig.savefig(bp_file+'.png', dpi=100, bbox_inches='tight', pad=0.05)
        plt.close()

if __name__ == "__main__":
    args = a.parse_args()

    gain_convention = 'divide'
    if args.multiply_gains:
        gain_convention = 'multiply'

    kwargs = copy.copy(args).__dict__
    kwargs.pop('fname')
    kwargs.pop('uv_file')
    kwargs.pop('multiply_gains')
    kwargs['verbose'] = args.silence == False
    kwargs.pop('silence')
    kwargs['gain_convention'] = gain_convention
    skynpz2calfits(args.fname, args.uv_file, **kwargs)

<file_sep>/aas/Estimating Delays/nn/network_trainers/NN_Trainer.py
import sys, os
from estdel.nn.networks import RestoreableComponent


# NN_Trainer
import matplotlib.pyplot as plt
import numpy as np
import os, io

class NN_Trainer(RestoreableComponent):
    """ Parent for restoreable neural network trainer (does not include training method).
        Add train method particular to network.

        Args:

            network - (one of the network classes from this repo) - The network to train
            Data_Creator - (one of the Data_Creator_* classes from this repo) - Creates data..
            num_epochs - (int) - Number of training epochs
            batch_size - (int) - Number of samples in each batch
                - Number of batches = int(number_of_samples/batch_size)
            log_dir - (string) - Directory to store logs & parameters
            model_save_interval - (int) - Save model every (this many) epochs
            pretrained_model_path - (string) - Path to previously trained model
            metric_names - (list of strings) - Names of the various metrics for this trainer
            verbose - (bool) - Be verbose.
    """
    __doc__ += RestoreableComponent.__doc__
    
    def __init__(self,
                 network,
                 Data_Creator, # Class
                 num_epochs,
                 batch_size,
                 log_dir,
                 model_save_interval,
                 pretrained_model_path,
                 metric_names,
                 verbose = True):
    
        RestoreableComponent.__init__(self, name=network.name, log_dir=log_dir, verbose=verbose)
        
        self._network = network   
        self._Data_Creator = Data_Creator
        self.num_epochs = num_epochs
        self.batch_size = batch_size
        self.model_save_interval = model_save_interval 
        self.pretrained_model_path = pretrained_model_path
        self._metrics = [] # list of lists
        self.metric_names = metric_names # list of strings, one per metric

    def add_data(self,train_info, test_info, gains, num_flatnesses = 10, abs_min_max_delay = 0.040):
        """Adds data to the Trainer.

            Args

                train_info - (tuple : (visibility data, baselines dict)) - The training data
                test_info - (tuple : (visibility data, baselines dict)) - The testing data
                gains - (dict) - The gains for the visibilites
                num_flatnesses - (int) - Number of flatnesses for each epoch
                    - Number of samples is num_flatnesses * 60
                abs_min_max_delay - (float) - The absolute value of the min or max delay for this dataset.

        """
        
        self._abs_min_max_delay = abs_min_max_delay

        self._train_batcher = self._Data_Creator(num_flatnesses,
                                                 train_info[0],
                                                 train_info[1],
                                                 gains,
                                                 abs_min_max_delay)
        self._train_batcher.gen_data()


        self._test_batcher = self._Data_Creator(num_flatnesses,
                                                test_info[0],
                                                test_info[1],
                                                gains,
                                                abs_min_max_delay)
        self._test_batcher.gen_data()
        
    def save_metrics(self):
        """Saves the recorded metrics to disk"""

        self._msg = '\rsaving metrics'; self._vprint(self._msg)
        direc = self.log_dir + self._network.name + '/'

        if not os.path.exists(direc): # should already exist by this point buuuuuut
            self._msg += ' - creating new directory ';self._vprint(self._msg)
            
            os.makedirs(direc)
        np.savez(direc + 'metrics', self.get_metrics()) 
        
        self._msg += ' - saved';self._vprint(self._msg)
        
        
    def get_metrics(self):
        """Returns the recorded metrics (list of lists)"""
        return {self.metric_names[i] : self._metrics[i] for i in range(len(self._metrics))}

    def plot_metrics(self,figsize = (8,6) ):
        """Plots the recorded metrics"""

        num_vals = np.min([len(metric) for metric in self._metrics])
        xvals = np.arange(num_vals)

        fig, axes = plt.subplots(len(self._metrics), 1, figsize = figsize, dpi = 144)

        for i, ax in enumerate(axes.reshape(-1)):

            if i == 0:
                ax.set_title('{}'.format(self._network.name))

            ax.plot(xvals[:num_vals], self._metrics[i][:num_vals], lw = 0.5)
            ax.set_xlabel('Epoch')
            ax.set_ylabel(self.metric_names[i])

        plt.tight_layout()
        plt.show()


        
    def gen_plot(self, predicted_values, actual_values, itx):
        """Create a prediction plot and save to byte string. For tensorboard images tab."""

        prediction_unscaled = itx(predicted_values)
        actual_unscaled = itx(actual_values)

        sorting_idx = np.argsort(actual_unscaled.T[0])

        fig, ax = plt.subplots(figsize = (5, 3), dpi = 144)

        ax.plot(prediction_unscaled.T[0][sorting_idx],
                linestyle = 'none', marker = '.', markersize = 1,
                color = 'darkblue')

        ax.plot(actual_unscaled.T[0][sorting_idx],
                linestyle = 'none', marker = '.', markersize = 1, alpha = 0.50,
                color = '#E50000')       

        ax.set_title('std: %.9f' %np.std(prediction_unscaled.T[0][sorting_idx] - actual_unscaled.T[0][sorting_idx]))

        buf = io.BytesIO()
        fig.savefig(buf, format='png', dpi = 144)
        plt.close(fig)
        buf.seek(0)

        return buf.getvalue()
    
    <file_sep>/aas/Estimating Delays/estdel/README.md
# estdel

## Prereqs

```
numpy
tensorflow
``` 
<file_sep>/aas/Estimating Delays/nn/data_creation/Data_Creator_BC_Odd_Even.py
# Data_Creator_BC_Odd_Even

from data_manipulation import *
from Data_Creator import Data_Creator

import numpy as np

class Data_Creator_BC_Odd_Even(Data_Creator):

    def __init__(self,
                 num_flatnesses,
                 bl_data = None,
                 bl_dict = None,
                 gains = None,
                 abs_min_max_delay = 0.040):

        Data_Creator.__init__(self,
                              num_flatnesses = num_flatnesses,
                              bl_data = bl_data,
                              bl_dict = bl_dict,
                              gains = gains,
                              abs_min_max_delay = abs_min_max_delay)
        
             
    def _gen_data(self):
        
        # scaling tools
        # the NN likes data in the range (0,1)
        angle_tx  = lambda x: (np.asarray(x) + np.pi) / (2. * np.pi)
        angle_itx = lambda x: np.asarray(x) * 2. * np.pi - np.pi

        delay_tx  = lambda x: (np.array(x) + self._tau) / (2. * self._tau)
        delay_itx = lambda x: np.array(x) * 2. * self._tau - self._tau
        
        targets = np.random.uniform(low = -self._tau, high = self._tau, size = (self._num * 60, 1))
        applied_delay = np.exp(-2j * np.pi * (targets * self._nu + np.random.uniform()))



        assert type(self._bl_data) != None, "Provide visibility data"
        assert type(self._bl_dict) != None, "Provide dict of baselines"
        assert type(self._gains)   != None, "Provide antenna gains"

        if self._bl_data_c == None:
            self._bl_data_c = {key : self._bl_data[key].conjugate() for key in self._bl_data.keys()}

        if self._gains_c == None:
            self._gains_c = {key : self._gains[key].conjugate() for key in self._gains.keys()}


        def _flatness(seps):
            """Create a flatness from a given pair of seperations, their data & their gains."""

            a, b = seps[0][0], seps[0][1]
            c, d = seps[1][0], seps[1][1]


            return self._bl_data[seps[0]]   * self._gains_c[(a,'x')] * self._gains[(b,'x')] * \
                   self._bl_data_c[seps[1]] * self._gains[(c,'x')]   * self._gains_c[(d,'x')]

        inputs = []
        for _ in range(self._num):

            unique_baseline = random.sample(self._bl_dict.keys(), 1)[0]
            two_seps = [random.sample(self._bl_dict[unique_baseline], 2)][0]

            inputs.append(_flatness(two_seps))
            

        inputs = np.angle(np.array(inputs).reshape(-1,1024) * applied_delay)
        
        permutation_index = np.random.permutation(np.arange(self._num * 60))
        
        # convert 0.0149 to 149 etc, for finding odd or even
        rounded_targets = 10000*np.asarray([np.round(abs(np.round(d * 100,2)/100), 5) for d in targets[permutation_index]]).reshape(-1)
        labels = [[1,0] if x%2 == 0 else [0,1] for x in rounded_targets]
        

        self._epoch_batch.append((angle_tx(inputs[permutation_index]), labels))

<file_sep>/aas/Estimating Delays/nn/networks/FNN_BN_R.py
"""FNN_BN_R
"""

import sys, os

from RestoreableComponent import RestoreableComponent

import tensorflow as tf
import numpy as np

class FNN_BN_R(RestoreableComponent):

    """FNN_BN_R

    A neural network of fully connected layers.
    
    Attributes:
        accuracy_threshold (float): the value for a target to be considered 'good'
        adam_initial_learning_rate (float): Adam optimizer initial learning rate
        cost (str): name of cost function. 'MSE', 'MQE', 'MISG', 'PWT_weighted_MSE', 'PWT_weighted_MISG'
            - use 'MSE', others experimental
        dtype (tensorflow obj): type used for all computations
        fcl_keep_prob (tensorflow obj): keep prob for fully connected layers
        gaussian_shift_scalar (float): value to shift gaussian for 'MISG' cost
        image_buf (tensorflow obj): buffer for image summaries for tensorboard
        is_training (tensorflow obj): flag for batch normalization
        layer_nodes (list of ints): numbers of nodes in each layer
        MISG (tensorflow obj): cost function, experimental. Mean Inverse Shifted Gaussian 
        MQE (tensorflow obj): cost function, experimental. Mean Quad Error
        MSE (tensorflow obj): cost function. Mean Squared Error
        optimizer (tensorflow obj): optimization function
        predictions (tensorflow obj): predicted outputs
        PWT (tensorflow obj): Percent of predictions within threshold
        sample_keep_prob (tensorflow obj): keep prob for input samples
        summary (tensorflow obj): summary operation for tensorboard
        targets (tensorflow obj): true values for optimizer
        X (tensorflow obj): input samples
    """
    __doc__ += RestoreableComponent.__doc__
    
    def __init__(self,
                 name,
                 layer_nodes,
                 cost = 'MSE',
                 log_dir = '../logs/',
                 dtype = tf.float32,
                 adam_initial_learning_rate = 0.0001,
                 accuracy_threshold = 0.00625,
                 gaussian_shift_scalar = 1e-5,
                 verbose = True):
        """__init__
        
        Args:
            name (str): name of network
            layer_nodes (list of ints): numbers of nodes in each layer
            cost (str, optional): name of cost function. 'MSE', 'MQE', 'MISG', 'PWT_weighted_MSE', 'PWT_weighted_MISG'
                    - use 'MSE', others experimental
            log_dir (str, optional): Directory to store network model and params
            dtype (tensorflow obj, optional): type used for all computations
            adam_initial_learning_rate (float, optional): Adam optimizer initial learning rate
            accuracy_threshold (float, optional): the value for a target to be considered 'good'
            gaussian_shift_scalar (float, optional): value to shift gaussian for 'MISG' cost
            verbose (bool, optional): Be verbose   
        """
        RestoreableComponent.__init__(self, name=name, log_dir=log_dir, verbose=verbose)
                
        self.layer_nodes = layer_nodes
        
        self.dtype = dtype
        self.adam_initial_learning_rate = adam_initial_learning_rate
        self.accuracy_threshold = accuracy_threshold
        self.gaussian_shift_scalar = gaussian_shift_scalar
        self.cost = cost
        

        self._num_freq_channels = 1024
        self._layers = []
        self._num_outputs = 1 # Regression


    def create_graph(self):
        """create_graph

        Create the network graph for use in a tensorflow session
        """
        self.save_params()
        self._msg = '\rcreating network graph '; self._vprint(self._msg)

        
        tf.reset_default_graph()
        self.is_training = tf.placeholder(dtype = tf.bool, shape = [], name = 'is_training')
        
        with tf.variable_scope('keep_probs'):
            
            self._msg += '.'; self._vprint(self._msg)

            self.sample_keep_prob = tf.placeholder(self.dtype, name = 'sample_keep_prob')
            self.fcl_keep_prob = tf.placeholder(self.dtype, name = 'fcl_keep_prob')
        
        
        with tf.variable_scope('sample'):
            
            self._msg += '.'; self._vprint(self._msg)
            
            self.X = tf.placeholder(self.dtype, shape = [None,self._num_freq_channels], name = 'X')
            self.X = tf.nn.dropout(self.X, self.sample_keep_prob)
            
        with tf.variable_scope('input_layer'):
            
            b = tf.get_variable(name = 'biases', shape = [self.layer_nodes[0]],
                                initializer = tf.contrib.layers.xavier_initializer())
            
            w = tf.get_variable(name = 'weights', shape  = [1024, self.layer_nodes[0]],
                                initializer = tf.contrib.layers.xavier_initializer())
            
            layer = tf.nn.leaky_relu(tf.matmul(self.X, w) + b)
            layer = tf.contrib.layers.batch_norm(layer, is_training = self.is_training)
            layer = tf.nn.dropout(layer, self.fcl_keep_prob)
            self._layers.append(layer)
            
        for i in range(len(self.layer_nodes)):
            if i > 0:
                with tf.variable_scope('layer_%d' %(i)):
                    layer = tf.contrib.layers.fully_connected(self._layers[i-1], self.layer_nodes[i])
                    layer = tf.contrib.layers.batch_norm(layer, is_training = self.is_training)
                    self._layers.append(layer)

                              
        
                
        self._msg += ' '
        with tf.variable_scope('targets'):
            self._msg += '.'; self._vprint(self._msg)
            
            self.targets = tf.placeholder(dtype = self.dtype, shape = [None, self._num_outputs], name = 'labels')

        with tf.variable_scope('predictions'):
            self._msg += '.'; self._vprint(self._msg)
            
            self.predictions = tf.contrib.layers.fully_connected(tf.layers.flatten(self._layers[-1]), self._num_outputs)
            
        with tf.variable_scope('costs'):
            self._msg += '.'; self._vprint(self._msg)

            error = tf.subtract(self.targets, self.predictions, name = 'error')
            squared_error = tf.square(error, name = 'squared_difference')
            quad_error = tf.square(squared_error, name = 'quad_error' )

            with tf.variable_scope('mean_inverse_shifted_gaussian'):
                self._msg += '.'; self._vprint(self._msg)

                normal_dist = tf.contrib.distributions.Normal(0.0, self.accuracy_threshold, name = 'normal_dist')
                gaussian_prob = normal_dist.prob(error, name = 'gaussian_prob')
                shifted_gaussian = tf.add(gaussian_prob, self.gaussian_shift_scalar, name = 'shifted_gaussian')   

                self.MISG = tf.reduce_mean(tf.divide(1.0, shifted_gaussian), name = 'mean_inverse_shifted_gaussian')
                
            with tf.variable_scope('mean_squared_error'):
                self._msg += '.'; self._vprint(self._msg)
                
                self.MSE = tf.reduce_mean(squared_error / self.gaussian_shift_scalar)
                
            with tf.variable_scope('mean_quad_error'):
                self._msg += '.'; self._vprint(self._msg)
                
                self.MQE = tf.reduce_mean(quad_error / self.gaussian_shift_scalar)

        with tf.variable_scope('logging'):  

            with tf.variable_scope('image'):
                self._msg += '.'; self._vprint(self._msg)
                
                self.image_buf = tf.placeholder(tf.string, shape=[])
                epoch_image = tf.expand_dims(tf.image.decode_png(self.image_buf, channels=4), 0)

            with tf.variable_scope('percent_within_threshold'):
                self._msg += '.'; self._vprint(self._msg)
                
                self.PWT = 100.*tf.reduce_mean(tf.cast(tf.less_equal(tf.abs(self.targets - self.predictions), self.accuracy_threshold), self.dtype) )


            tf.summary.histogram(name = 'targets', values = self.targets)
            tf.summary.histogram(name = 'predictions',values =  self.predictions)
            tf.summary.scalar(name = 'MSE', tensor = self.MSE)
            tf.summary.scalar(name = 'MISG', tensor = self.MISG)
            tf.summary.scalar(name = 'MQE', tensor = self.MQE)
            tf.summary.scalar(name = 'PWT', tensor = self.PWT)
            tf.summary.image('prediction_vs_actual', epoch_image)
            self.summary = tf.summary.merge_all()
            
        update_ops = tf.get_collection(tf.GraphKeys.UPDATE_OPS)
        with tf.control_dependencies(update_ops):
            with tf.variable_scope('train'):
                self._msg += '.'; self._vprint(self._msg)
            
                if self.cost == 'MSE':
                    cost = self.MSE
                if self.cost == 'MQE':
                    cost = tf.log(self.MQE)
                if self.cost == 'MISG':
                    cost = self.MISG
                if self.cost == 'PWT_weighted_MSE':
                    cost = self.MSE * (100. - self.PWT)
                if self.cost == 'PWT_weighted_MISG':
                    cost = self.MISG * (100. - self.PWT)


                self.optimizer = tf.train.AdamOptimizer(self.adam_initial_learning_rate, epsilon=1e-08).minimize(cost)
            
        
        num_trainable_params = np.sum([np.prod(v.get_shape().as_list()) for v in tf.trainable_variables()])
        self._msg = '\rNetwork Ready - {} trainable parameters'.format(num_trainable_params); self._vprint(self._msg)<file_sep>/aas/Estimating Delays/nn/networks/CNN_QP_BN_R.py
"""CNN_QP_BN_R
"""

import sys, os

from RestoreableComponent import RestoreableComponent

import tensorflow as tf
import numpy as np

class CNN_QP_BN_R(RestoreableComponent):
    """CNN_QP_BN_R
    CNN: Convolutional neural network.
   QP: Each computational layer is a quad-path layer.
   BN: All non-linearalities have batch-normalization applied.
   R: Regression, this network predicts a single value for each input.

    
    Attributes:
        accuracy (tensorflow obj): Running accuracy of predictions
        adam_initial_learning_rate (float): Adam optimizer initial learning rate
        conv_keep_prob (tensorflow obj): Keep prob rate for convolutions
        cost (str): name of cost function. 'MSE', 'MQE', 'MISG', 'PWT_weighted_MSE', 'PWT_weighted_MISG'
            - use 'MSE', others experimental
        dtype (tensorflow obj): Type used for all computations
        gaussian_shift_scalar (float): value to shift gaussian for 'MISG' cost
        image_buf (tensorflow obj): buffer for image summaries for tensorboard
        is_training (tensorflow obj): flag for batch normalization
        layer_downsampling_factors (list of ints): factors to downsample each layer
        MISG (tensorflow obj): cost function, experimental. Mean Inverse Shifted Gaussian 
        MSE (tensorflow obj): cost function. Mean Squared Error
        num_1x1_conv_filters_per_layer (list of ints): number of 1x1 convolutions for each layer
        num_freq_channels (int): Number of frequency channels
        optimizer (tensorflow obj): optimization function
        pred_keep_prob (TYPE): prediction keep prob
        predictions (tensorflow obj): predicted outputs
        PWT (tensorflow obj): Percent of predictions within threshold
        sample_keep_prob (tensorflow obj): keep prob for input samples
        samples (TYPE): input samples
        summary (tensorflow obj): summary operation for tensorboard
        targets (tensorflow obj): true values for optimizer
        wide_convolution_filter_widths (list of ints): widths of each laters wide convolutions
    """
    __doc__ += RestoreableComponent.__doc__
    def __init__(self,
                 name,
                 wide_convolution_filter_widths,
                 layer_downsampling_factors, 
                 num_1x1_conv_filters_per_layer,
                 log_dir = '../logs/',
                 dtype = tf.float32,
                 adam_initial_learning_rate = 0.0001,
                 cost = 'MSE',
                 accuracy_threshold = 0.00625,
                 gaussian_shift_scalar = 1e-5,
                 verbose = True):
        """__init__
        
        Args:
            name (str): name of network
            wide_convolution_filter_widths (list of ints): widths of each laters wide convolutions
            layer_downsampling_factors (list of ints): factors to downsample each layer
            num_1x1_conv_filters_per_layer (list of ints): number of 1x1 convolutions for each layer
            log_dir (str, optional): directory to store network model and params
            dtype (tensorflow obj): Type used for all computations
            adam_initial_learning_rate (tensorflow obj): Adam optimizer initial learning rate
            cost (str): name of cost function. 'MSE', 'MQE', 'MISG', 'PWT_weighted_MSE', 'PWT_weighted_MISG'
                - use 'MSE', others experimental
            accuracy_threshold (float, optional): Threshold to count a prediction as being on target
            gaussian_shift_scalar (float): value to shift gaussian for 'MISG' cost
            verbose (bool, optional): be verbose
        """
        RestoreableComponent.__init__(self, name=name, log_dir=log_dir, verbose=verbose)
                
        self.wide_convolution_filter_widths = wide_convolution_filter_widths
        self.layer_downsampling_factors = layer_downsampling_factors
        
        self.dtype = dtype
        self.adam_initial_learning_rate = adam_initial_learning_rate
        self.cost = cost
        self.accuracy_threshold = accuracy_threshold
        self.gaussian_shift_scalar = gaussian_shift_scalar
        self.num_1x1_conv_filters_per_layer = num_1x1_conv_filters_per_layer 
        
        self.num_freq_channels = 1024
        
    def create_graph(self):
        """create_graph

        Create the network graph for use in a tensorflow session
        """
        self.save_params()
        self._msg = '\rcreating network graph '; self._vprint(self._msg)
        
        tf.reset_default_graph()
        
        self.is_training = tf.placeholder(dtype = tf.bool, name = 'is_training', shape = [])
        
        with tf.variable_scope('keep_probs'):
            self._msg += '.'; self._vprint(self._msg)
            
            self.sample_keep_prob = tf.placeholder(self.dtype, name = 'sample_keep_prob')
            self.conv_keep_prob = tf.placeholder(self.dtype, name = 'conv_keep_prob') 
            self.pred_keep_prob = tf.placeholder(self.dtype, name = 'pred_keep_prob')    

        with tf.variable_scope('samples'):
            self._msg += '.'; self._vprint(self._msg)
            
            self.samples = tf.placeholder(self.dtype, shape = [None, 1, self.num_freq_channels, 1], name = 'samples')
            self.samples = tf.nn.dropout(self.samples, self.sample_keep_prob)

        self._layers = []
        num_layers = len(self.wide_convolution_filter_widths)
        layer_names = ['layer_{}'.format(i) for i in range(num_layers)]

        for i in range(num_layers):
            self._msg += '.'; self._vprint(self._msg)
            
            # previous layer is input for current layer
            layer = self.samples if i == 0 else self._layers[i - 1]
            layer = Quad_Path_Layer(layer,
                                    layer_names[i],
                                    self.wide_convolution_filter_widths[i],
                                    self.layer_downsampling_factors[i],
                                    self.num_1x1_conv_filters_per_layer[i],
                                    self.conv_keep_prob,
                                    self.is_training,
                                    self.dtype)
            layer = layer.process()
            
            self._layers.append(layer)
            
        with tf.variable_scope('prediction'):
            self._msg += '.'; self._vprint(self._msg)
            reshaped_final_layer = tf.contrib.layers.flatten(self._layers[-1])
            reshaped_final_layer = tf.nn.dropout(reshaped_final_layer, self.pred_keep_prob)
            prediction_weight = tf.get_variable(name = 'weight',
                                                shape = [reshaped_final_layer.get_shape()[-1], 1],
                                                dtype = self.dtype,
                                                initializer = tf.contrib.layers.xavier_initializer())
            
            self.predictions = tf.matmul(reshaped_final_layer, prediction_weight)

        with tf.variable_scope('targets'):
            self._msg += '.'; self._vprint(self._msg)
            
            self.targets = tf.placeholder(dtype = self.dtype, shape = [None, 1], name = 'targets')

        with tf.variable_scope('costs'):
            self._msg += '.'; self._vprint(self._msg)

            error = tf.subtract(self.targets, self.predictions, name = 'error')
            squared_error = tf.square(error, name = 'squared_difference')

            with tf.variable_scope('mean_inverse_shifted_gaussian'):
                self._msg += '.'; self._vprint(self._msg)

                normal_dist = tf.contrib.distributions.Normal(0.0, self.accuracy_threshold, name = 'normal_dist')
                gaussian_prob = normal_dist.prob(error, name = 'gaussian_prob')
                shifted_gaussian = tf.add(gaussian_prob, self.gaussian_shift_scalar, name = 'shifted_gaussian')   

                self.MISG = tf.reduce_mean(tf.divide(1.0, shifted_gaussian), name = 'mean_inverse_shifted_gaussian')
                
            with tf.variable_scope('mean_squared_error'):
                self._msg += '.'; self._vprint(self._msg)
                
                self.MSE = tf.reduce_mean(squared_error)


        with tf.variable_scope('logging'):  

            with tf.variable_scope('image'):
                self._msg += '.'; self._vprint(self._msg)
                
                self.image_buf = tf.placeholder(tf.string, shape=[])
                epoch_image = tf.expand_dims(tf.image.decode_png(self.image_buf, channels=4), 0)

            with tf.variable_scope('percent_within_threshold'):
                self._msg += '.'; self._vprint(self._msg)
                
                self.PWT = 100.*tf.reduce_mean(tf.cast(tf.less_equal(tf.abs(self.targets - self.predictions), self.accuracy_threshold), self.dtype) )


            tf.summary.histogram(name = 'targets', values = self.targets)
            tf.summary.histogram(name = 'predictions',values =  self.predictions)
            tf.summary.scalar(name = 'MSE', tensor = self.MSE)
            tf.summary.scalar(name = 'MISG', tensor = self.MISG)
            tf.summary.scalar(name = 'PWT', tensor = self.PWT)
            tf.summary.image('prediction_vs_actual', epoch_image)
            self.summary = tf.summary.merge_all()
            
        update_ops = tf.get_collection(tf.GraphKeys.UPDATE_OPS)
        with tf.control_dependencies(update_ops):
            with tf.variable_scope('train'):
                self._msg += '.'; self._vprint(self._msg)
            
                if self.cost == 'MSE':
                    cost = self.MSE
                if self.cost == 'MISG':
                    cost = self.MISG
                if self.cost == 'PWT_weighted_MSE':
                    cost = self.MSE * (100. - self.PWT)
                if self.cost == 'PWT_weighted_MISG':
                    cost = self.MISG * (100. - self.PWT)


                self.optimizer = tf.train.AdamOptimizer(self.adam_initial_learning_rate, epsilon=1e-08).minimize(cost)
            
        
        num_trainable_params = np.sum([np.prod(v.get_shape().as_list()) for v in tf.trainable_variables()])
        self._msg = '\rNetwork Ready - {} trainable parameters'.format(num_trainable_params); self._vprint(self._msg)



class Quad_Path_Layer(object):
    """A layer for a convolutional network. Layer is made of four paths: Average, Max, Wide, Narrow.
    
     Average - Average pool followed by a 1x1 convolution.
     Max - Max pool followed by a 1x1 convolution.
     Wide - 1x1 convolution followed by 1xWide convoltuion.
     Narrow - 1x1 convolution followed 1x(Wide/2) convolution.
    
    
    Attributes:
        conv_keep_prob (float): keep prob for convolutions
        dtype (tensorflow obj): type used for all computations
        is_training (tensorflow obj): flag for batch normalization
        layer_name (str): name of layer
        num_1x1_conv_filters (int): Number of 1x1 convolution filters
        strides (list of ints): stride for this layer
        t (tensor): tensor to process
        wide_convolution_filter_width (int): Width of the wide convolution filter
    
    
    """
    
    def __init__(self,
                 t,
                 layer_name,
                 wide_convolution_filter_width,
                 layer_downsampling_factor,
                 num_1x1_conv_filters = 4,
                 conv_keep_prob = 0.90,
                 is_training = True,
                 dtype = tf.float32):
        """Summary
        
        Args:
            t (tensor): tensor to process
            layer_name (str): name of layer
            wide_convolution_filter_width (int): Width of the wide convolution filter
            layer_downsampling_factor (int): Factor to downsample by
            num_1x1_conv_filters (int, optional): Number of 1x1 convolution filters
            conv_keep_prob (float, optional): keep prob for convolutions
            is_training (bool, optional): flag for batch normalization
            dtype (tensorflow obj, optional): type used for all computations
        """
        self.t = t
        self.layer_name = layer_name
        self.wide_convolution_filter_width = wide_convolution_filter_width
        self.strides = [1,1,layer_downsampling_factor, 1]
        self.num_1x1_conv_filters = num_1x1_conv_filters
        self.conv_keep_prob = conv_keep_prob
        self.is_training = is_training
        self.dtype = dtype

        
    def process(self):
        """process

        Creates the quad path layer
        
        Returns:
            tensor: concatenated filtes after processing
        """
        with tf.variable_scope(self.layer_name):

            narrow_convolution_filter_width = self.wide_convolution_filter_width / 2

            num_narrow_conv_filters = np.max([self.num_1x1_conv_filters / 2, 2])
            num_wide_conv_filters = np.max([num_narrow_conv_filters / 2, 1])


            path_A = self._avg_scope()
            path_B = self._max_scope()

            path_C_and_D = self._1x1_conv(self.t)

            path_C = self._conv_scope(path_C_and_D,
                                      [1,
                                       narrow_convolution_filter_width,
                                       path_C_and_D.get_shape().as_list()[3],
                                       num_narrow_conv_filters],
                                      self.strides,
                                      scope_name = 'narrow')

            path_D = self._conv_scope(path_C_and_D,
                                      [1,
                                       self.wide_convolution_filter_width,
                                       path_C_and_D.get_shape().as_list()[3],
                                       num_wide_conv_filters],
                                      self.strides,
                                      scope_name = 'wide')

            t = self._filter_cat_scope([path_A, path_B, path_C, path_D])

        return t

    def _trainable(self, name, shape):
        """_trainable

        Wrapper for tensorflow.get_variable(), xavier initialized
        
        Args:
            name (str): name of variable
            shape (int or list of ints): shape for vairable
        
        Returns:
            tensorflow obj: trainable variable
        """
        return tf.get_variable(name = name,
                               dtype = self.dtype,
                               shape = shape,
                               initializer = tf.contrib.layers.xavier_initializer())
    
    def _bias_add_scope(self, t, shape):
        """_bias_add_scope

        Creates a scope around a trainable bias and its addition to input
        
        Args:
            t (tensor): tensor to add biases to
            shape (int): number of biases to add
        
        Returns:
            tensor: tensor with biases added to it
        """
        with tf.variable_scope('add_bias'):

            biases = self._trainable('biases', shape)
            t = tf.nn.bias_add(t, biases)

        return t


    def _conv_scope(self,t, filter_shape, strides, scope_name = 'convolution'):
        """_conv_scope

        Creates a scope around a convolution layer.
        t = dropout( batch_norm( relu( conv2d(t) + bias )))
        Args:
            t (tensor): tensor to process
            filter_shape (list of ints): Shape of convolution filters
            strides (list of ints): Convolution strides
            scope_name (str, optional): name of scope (for graph organization)
        
        Returns:
            tensor: processed tensor
        """
        with tf.variable_scope(scope_name):

            t = tf.nn.conv2d(t, self._trainable('filters', filter_shape),strides,'SAME')
            t = self._bias_add_scope(t, [filter_shape[-1]])
            t = tf.nn.relu(t)
            t = tf.contrib.layers.batch_norm(t, is_training = self.is_training)
            t = tf.nn.dropout(t, self.conv_keep_prob)

        return t
    
    def _1x1_conv(self, t):
        """_1x1_conv

        1x1 convolution with strides = 1 and num_1x1_conv_filters filters
        
        Args:
            t (tensor): tensor to convolve
        
        Returns:
            tensor: convolved tensor
        """
        return self._conv_scope(t,
                                [1,1,t.get_shape().as_list()[3],self.num_1x1_conv_filters],
                                [1,1,1,1],
                                "1x1_conv")

    def _avg_scope(self):
        """_avg_scope

        Creates a scope around the average-pool path.
        
        Returns:
            tensor: tensor after average pool and 1x1 convolution
        """
        with tf.variable_scope('average'):
            t = tf.nn.avg_pool(self.t, self.strides, self.strides, padding = "SAME")
            t = self._1x1_conv(t)

        return t

    def _max_scope(self):
        """_max_scope

        Creates a scope around the max-pool path
        
        Returns:
            tensor: tensor after max pool and 1x1 convolution
        """
        with tf.variable_scope('max'):
            t = tf.nn.max_pool(self.t, self.strides, self.strides, padding = "SAME")
            t = self._1x1_conv(t)

        return t

    def _filter_cat_scope(self,t_list):
        """_filter_cat_scope

        Creates a scope around filter concatation (layer output)
        
        Args:
            t_list (list of tensors): filters to concatenate
        
        Returns:
            tensor: concatenated tensors
        """
        with tf.variable_scope('filter_cat'):
            t = tf.concat(t_list, 3)
        return t