text:
import numpy as np
import time  # You can use this to time how long different parts take to run--find inefficiencies
from test import X, y  # Data set and labels used for testing purposes.
#np.random.seed(0)
np.set_printoptions(threshold=100000)

def get_data_from_acath_csv():
    """
    Data found at: http://biostat.mc.vanderbilt.edu/wiki/Main/DataSets

    For description of data, see: http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/acath.html

    For label explanations, see: http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/Cacath.html
    I've edited the column names to better reflect what the values mean

    Outputs:
    [(training_data, training_sigdz_labels, training_tvdlm_labels),
    (validation_data, validation_sigdz_labels, validation_tvdlm_labels),
    (testing_data, testing_sigdz_labels, testing_tvdlm_labels)]

    Data is shuffled before being allocated to training/validation/testing.
    """
    lines = []
    for line in open('acath.csv'):
        line = (line.replace('\n', ''))
        line = line.split(',')
        lines.append(line)

    lines = lines[1:]
    np.random.shuffle(lines)  # shuffles the same way each time its run
    samples = []
    significant_coronary_disease_labels = []
    three_vessel_or_left_main_disease_labels = []
    for line in lines:
        if line[3] == '':
            line[3] = '0'  # If there's no cholesterol data, set it to 0. There are 1246 times we do this.
        if not line[5] == '':  # Excluding data points with no data for three vessel or left main disease
            # There are 3 empty values for this, so we are ignoring those 3 samples. See:
            # http://biostat.mc.vanderbilt.edu/wiki/pub/Main/DataSets/Cacath.html
            sample = [float(val) for val in line]  # Turning string values into actual floats
            samples.append(sample[0:4])  # Including sex, age, symptom duration, and cholesterol as features
            significant_coronary_disease_labels.append(int(sample[4]))
            three_vessel_or_left_main_disease_labels.append(int(sample[5]))

    # Since the data has already been shuffled, we'll just take out the pieces we want for the data sets

    training_data = np.array(samples[0:int(len(samples)*0.6)])
    training_sigdz_labels = np.array(significant_coronary_disease_labels[0:int(len(samples)*0.6)])
    training_tvdlm_labels = np.array(three_vessel_or_left_main_disease_labels[0:int(len(samples)*0.6)])

    validation_data = np.array(samples[round(len(samples)*0.6):round(len(samples)*0.8)])
    validation_sigdz_labels = np.array(significant_coronary_disease_labels
                                       [round(len(samples)*0.6):round(len(samples)*0.8)])
    validation_tvdlm_labels = np.array(three_vessel_or_left_main_disease_labels
                                       [round(len(samples)*0.6):round(len(samples)*0.8)])

    testing_data = np.array(samples[round(len(samples)*0.8):])
    testing_sigdz_labels = np.array(significant_coronary_disease_labels[round(len(samples)*0.8):])
    testing_tvdlm_labels = np.array(three_vessel_or_left_main_disease_labels[round(len(samples)*0.8):])

    # normalizing the data
    training_data = training_data/np.max(training_data, axis=0)
    validation_data = validation_data/np.max(validation_data, axis=0)
    testing_data = testing_data/np.max(testing_data, axis=0)

    train = (training_data, training_sigdz_labels, training_tvdlm_labels)
    validation = (validation_data, validation_sigdz_labels, validation_tvdlm_labels)
    test = (testing_data, testing_sigdz_labels, testing_tvdlm_labels)

    return [train, validation, test]

def get_data_from_falldetection_csv():
    """
    Data set found at https://www.kaggle.com/pitasr/falldata/version/1#
    """
    lines = []
    for line in open('falldetection.csv'):
        line = (line.replace('\n', ''))
        line = line.split(',')
        lines.append(line)
    column_names = lines.pop(0)
    np.random.shuffle(lines)
    labels = []
    samples = []
    for line in lines:
        labels.append(int(line[0]))
        samples.append([float(val) for val in line[1:]])

    all_data = [(label, sample) for label, sample in zip(labels, samples)]
    organized_data = {}
    for data_point in all_data:
        if data_point[0] in organized_data:
            organized_data[data_point[0]].append(data_point[1])
        else:
            organized_data[data_point[0]] = [data_point[1]]

    training_dict = {}
    rest_dict = {}
    for label in organized_data:
        if label != 3:
            training_dict[label] = organized_data[label][0:int(len(organized_data[3]) * 0.8/5)]
            rest_dict[label] = organized_data[label][int(len(organized_data[3]) * 0.8/5):]
    training_dict[3] = organized_data[3][0:int(len(organized_data[3]) * 0.8)]
    rest_dict[3] = organized_data[3][int(len(organized_data[3]) * 0.8):]

    training_data_and_labels = []
    rest_data_and_labels = []
    for key in training_dict:
        for sample in training_dict[key]:
            if key == 3:
                training_data_and_labels.append((1, sample))
            else:
                training_data_and_labels.append((0, sample))

        for sample2 in rest_dict[key]:
            if key == 3:
                rest_data_and_labels.append((1, sample2))
            else:
                rest_data_and_labels.append((0, sample2))

    np.random.shuffle(training_data_and_labels)
    np.random.shuffle(rest_data_and_labels)
    rest_data = [tup[1] for tup in rest_data_and_labels]
    rest_labels = [tup[0] for tup in rest_data_and_labels]

    training_data = np.array([tup[1] for tup in training_data_and_labels])
    training_labels = np.array([tup[0] for tup in training_data_and_labels])

    validation_data = np.array(rest_data[0:int(len(rest_data) * 0.3)])
    validation_labels = np.array(rest_labels[0:int(len(rest_data) * 0.3)])

    testing_data = np.array(rest_data[int(len(rest_data) * 0.3):])
    testing_labels = np.array(rest_labels[int(len(rest_data) * 0.3):])

    # normalizing the data
    #training_data = training_data / np.max(training_data, axis=0)
    #validation_data = validation_data / np.max(validation_data, axis=0)
    #testing_data = testing_data / np.max(testing_data, axis=0)

    train = (training_data, training_labels)
    validation = (validation_data, validation_labels)
    test = (testing_data, testing_labels)

    return [train, validation, test]

def get_data_from_heart_failure_clinical_records_dataset_csv():
    lines = []
    for line in open('heart_failure_clinical_records_dataset.csv'):
        lines.append(line.replace('\n', '').split(','))
    lines.pop(0)
    np.random.shuffle(lines)
    samples = []
    labels = []
    for line in lines:
        samples.append([float(val) for val in line[:-1]])
        labels.append(int(line[-1]))

    edited_samples = [sample[0:-1] for sample in samples]  # removes time
    assert len(labels) == len(samples)
    training_data = np.array(samples[0:int(len(samples)*0.7)])
    training_labels = np.array(labels[0:int(len(samples)*0.7)])

    validation_data = np.array(samples[int(len(samples)*0.7):int(len(samples)*0.8)])
    validation_labels = np.array(labels[int(len(samples)*0.7):int(len(samples)*0.8)])

    testing_data = np.array(samples[int(len(samples)*0.8):])
    testing_labels = np.array(labels[int(len(samples)*0.8):])

    training_data = training_data / np.max(training_data, axis=0)
    validation_data = validation_data / np.max(validation_data, axis=0)
    testing_data = testing_data / np.max(testing_data, axis=0)

    train = (training_data, training_labels)
    validate = (validation_data, validation_labels)
    test = (testing_data, testing_labels)

    return train, validate, test

class DenseLayer:
    def __init__(self, num_of_inputs, num_of_neurons, l2_weight_lambda=0.0,
                 l2_bias_lambda=0.0):
        self.weights = 0.01 * np.random.randn(num_of_inputs, num_of_neurons)
        # The parameters of randn are just the dimensions of the matrix it makes
        self.biases = np.zeros((1, num_of_neurons))
        # creates a matrix of height 1 and length num_of_neurons filled with 0's
        self.output = None

        # Factors used to penalize weights/biases of large magnitude. See the regularization_loss method in
        # the CE loss class
        self.l2_weight_lambda = l2_weight_lambda
        self.l2_bias_lambda = l2_bias_lambda

    def forward(self, inputs):
        self.inputs = inputs
        self.output = np.dot(self.inputs, self.weights) + self.biases

    def backward(self, d_error__d_denselayer):
        d_denselayer__d_weights = self.inputs.T
        d_denselayer__d_biases = np.ones((1, len(self.inputs)))
        # ^^ produces a matrix of 1x(number of samples) dimensions

        self.d_error__d_weights = np.dot(d_denselayer__d_weights,
                                         d_error__d_denselayer)
        if self.l2_weight_lambda > 0:
            self.d_error__d_weights += 2*self.l2_weight_lambda*self.weights
        self.d_error__d_biases = np.dot(d_denselayer__d_biases,
                                        d_error__d_denselayer)
        if self.l2_bias_lambda > 0:
            self.d_error__d_biases += 2*self.l2_bias_lambda*self.biases

        d_denselayer__d_input = self.weights.T
        self.d_error__d_inputs = np.dot(d_error__d_denselayer, d_denselayer__d_input)

class ActivationReLU:
    def __init__(self):
        self.output = None

    def forward(self, inputs):
        self.inputs = inputs
        self.output = np.maximum(0, inputs)

    def backward(self, d_error__d_relu):
        d_relu__d_inputs = np.ones((self.output.shape[0], self.output.shape[1]))

        d_relu__d_inputs[self.inputs <= 0] = 0
        # The above line is a fast implementation of this, below is a slower but more intuitive implementation

        '''
        for row_index in range(len(d_relu__d_inputs)):
            row = d_relu__d_inputs[row_index]
            for val_index in range(len(row)):
                if self.output[row_index][val_index] <= 0:
                    d_relu__d_inputs[row_index][val_index] = 0
        '''

        # element wise product
        self.d_error__d_inputs = np.multiply(d_error__d_relu, d_relu__d_inputs)

class ActivationSigmoid:
    def forward(self, inputs):
        self.inputs = inputs
        self.output = 1/(1+np.e**(-self.inputs))

    def backward(self, d_error__d_sigmoid):
        d_sigmoid__d_inputs = np.multiply(self.output, (1-self.output))
        # element wise product
        self.d_error__d_inputs = np.multiply(d_error__d_sigmoid, d_sigmoid__d_inputs)

class ActivationSoftmax:
    def __init__(self):
        self.output = None

    def forward(self, inputs):
        exponentiated_values = np.exp(inputs-np.max(inputs, axis=1,
                                                    keepdims=True))
        # the np.max part ensures all exponents are <=0, while actually
        # resulting in the same probability outputs!
        self.output = exponentiated_values / np.sum(exponentiated_values,
                                                    axis=1, keepdims=True)
        # we use axis=1 in both lines to make sure it does all these operations
        # based on rows.

    def backward(self, d_error__d_softmax):
        num_of_nodes = self.output.shape[1]

        gradient_matrices_for_d_softmax__d_inputs = []
        '''
        for sample in self.output:
            gradient_matrix = np.tile(sample, (len(sample), 1))  # creates a matrix of len(sample)xlen(sample)
                                                               # dim. where each row is = to sample
            for row_val_index in range(num_of_nodes):
                gradient_matrix[row_val_index] *= -sample[row_val_index]
                gradient_matrix[row_val_index][row_val_index] = sample[row_val_index]*(1-sample[row_val_index])

            gradient_matrices_for_d_softmax__d_inputs.append(gradient_matrix)

        '''  # An intuitive implementation, though less efficient

        for sample in self.output:
            gradient_matrix = np.zeros((num_of_nodes, num_of_nodes))
            for row_index in range(len(gradient_matrix)):
                for val_index in range(len(gradient_matrix[0])):
                    if row_index == val_index:
                        gradient_matrix[row_index][val_index] = sample[val_index] * (1 - sample[row_index])
                    else:
                        gradient_matrix[row_index][val_index] = -sample[row_index] * sample[val_index]

            gradient_matrices_for_d_softmax__d_inputs.append(gradient_matrix)

        d_softmax__d_inputs = np.array(gradient_matrices_for_d_softmax__d_inputs)

        gradient_matrices_for_d_error__d_inputs = []
        for row, sub_matrix in zip(d_error__d_softmax, d_softmax__d_inputs):
            gradient_matrices_for_d_error__d_inputs.append(
                np.dot(sub_matrix, row))

        self.d_error__d_inputs = np.array(gradient_matrices_for_d_error__d_inputs)

class LossCategoricalCrossEntropy:
    def __init__(self):
        # obligatory iS tHiS lOsS reference
        """
        │     │ |

        │ │   │ __
        """
        pass

    def forward(self, network_output, correct_output):
        self.network_output = np.clip(network_output, 1e-7, 1 - 1e-7)
        self.correct_output = correct_output
        # correct_output are the correct
        # classifications. If correct_output is [0, 1, 1], then that means the
        # first sample was category 0, the second sample was cat. 1, etc.

        neg_logs = []
        for output, actual in zip(self.network_output, self.correct_output):
            neg_logs.append(-np.log(output[actual]))

        self.error = np.mean(neg_logs)

    def backward(self):
        self.d_error__d_inputs = np.zeros((self.network_output.shape[0],
                                    self.network_output.shape[1]))
        for sample_index in range(len(self.network_output)):
            correct_output_class_for_this_sample = self.correct_output[sample_index]
            self.d_error__d_inputs[sample_index][correct_output_class_for_this_sample] = \
                -1/(self.network_output[sample_index][correct_output_class_for_this_sample] *
                    self.network_output.shape[0])

    def l2_regularization(self, layer):
        regularization_loss = 0
        if layer.l2_weight_lambda > 0:
            regularization_loss += layer.l2_weight_lambda * np.sum(layer.weights**2)
        if layer.l2_bias_lambda > 0:
            regularization_loss += layer.l2_bias_lambda * np.sum(layer.biases**2)

        return regularization_loss

class OptimizerVanillaSGD:
    def __init__(self, learning_rate=1.0, decay=0.1):
        self.starting_learning_rate = learning_rate
        self.current_learning_rate = learning_rate
        self.decay = decay
        self.iterations = 0

    def update_parameters(self, list_of_layers):
        self.current_learning_rate = self.starting_learning_rate / \
                                     (1+(self.decay*self.iterations))

        for layer_object in list_of_layers:
            update_weights = -(self.current_learning_rate * layer_object.d_error__d_weights)

            update_biases = -(self.current_learning_rate * layer_object.d_error__d_biases)

            layer_object.weights += update_weights
            layer_object.biases += update_biases

        self.iterations += 1

class OptimizerSGDWithMomentum:
    def __init__(self, learning_rate=1.0, decay=0.1, momentum=0.0):
        self.starting_learning_rate = learning_rate
        self.current_learning_rate = learning_rate
        self.decay = decay
        self.iterations = 0
        self.momentum = momentum

    def update_parameters(self, list_of_layers):
        self.current_learning_rate = self.starting_learning_rate / \
                                     (1+(self.decay*self.iterations))

        for layer_object in list_of_layers:
            if not hasattr(layer_object, 'previous_update_weights'):

                layer_object.previous_update_weights = np.zeros_like(layer_object.weights)
                layer_object.previous_update_biases = np.zeros_like(layer_object.biases)

            update_weights = (self.momentum * layer_object.previous_update_weights) +\
                               -(self.current_learning_rate * layer_object.d_error__d_weights)

            update_biases = (self.momentum * layer_object.previous_update_biases) +\
                            -(self.current_learning_rate * layer_object.d_error__d_biases)

            layer_object.previous_update_weights = update_weights
            layer_object.previous_update_biases = update_biases

            layer_object.weights += update_weights
            layer_object.biases += update_biases

        self.iterations += 1

class OptimizerAdaGrad:
    def __init__(self, learning_rate=1.0, decay=0.1, epsilon=1e-7):
        self.starting_learning_rate = learning_rate
        self.current_learning_rate = learning_rate
        self.decay = decay
        self.iterations = 0
        self.epsilon = epsilon

    def update_parameters(self, list_of_layers):
        self.current_learning_rate = self.starting_learning_rate / \
                                     (1+(self.decay*self.iterations))

        for layer_object in list_of_layers:
            if not (hasattr(layer_object, 'weight_cache') or hasattr(layer_object, 'bias_cache')):

                layer_object.weight_cache = np.zeros_like(layer_object.weights)
                layer_object.bias_cache = np.zeros_like(layer_object.biases)

            layer_object.weight_cache += layer_object.d_error__d_weights**2
            layer_object.bias_cache += layer_object.d_error__d_biases**2

            weight_learning_rates = self.current_learning_rate /\
                                   (np.sqrt(layer_object.weight_cache)+self.epsilon)
            bias_learning_rates = self.current_learning_rate /\
                                   (np.sqrt(layer_object.bias_cache)+self.epsilon)

            layer_object.weights += (-weight_learning_rates * layer_object.d_error__d_weights)

            layer_object.biases += (-bias_learning_rates * layer_object.d_error__d_biases)

        self.iterations += 1

class OptimizerRMSProp:
    def __init__(self, learning_rate=.001, decay=0.1, epsilon=1e-7, rho=0.9):
        # NOT SURE WHY WE SET SUCH A LOW STARTING LEARNING RATE
        self.starting_learning_rate = learning_rate
        self.current_learning_rate = learning_rate
        self.decay = decay
        self.iterations = 0
        self.epsilon = epsilon
        self.rho = rho

    def update_parameters(self, list_of_layers):
        self.current_learning_rate = self.starting_learning_rate / \
                                     (1+(self.decay*self.iterations))

        for layer_object in list_of_layers:
            if not (hasattr(layer_object, 'weight_cache') or hasattr(layer_object, 'bias_cache')):

                layer_object.weight_cache = np.zeros_like(layer_object.weights)
                layer_object.bias_cache = np.zeros_like(layer_object.biases)

            # Think of rho as the percent of the cache that we keep around. If we are changing the cache
            # less (i.e. higher rho value) then the cache value changes in a smoother way.
            # Smoother cache => smoother learning rates => smoother changes to weights/biases.
            layer_object.weight_cache = self.rho*layer_object.weight_cache + \
                                         (1-self.rho)*(layer_object.d_error__d_weights**2)

            layer_object.bias_cache = self.rho*layer_object.bias_cache + \
                                       (1-self.rho)*(layer_object.d_error__d_biases**2)

            weight_learning_rates = self.current_learning_rate /\
                                   (np.sqrt(layer_object.weight_cache)+self.epsilon)
            bias_learning_rates = self.current_learning_rate /\
                                   (np.sqrt(layer_object.bias_cache)+self.epsilon)

            layer_object.weights += (-weight_learning_rates * layer_object.d_error__d_weights)

            layer_object.biases += (-bias_learning_rates * layer_object.d_error__d_biases)

        self.iterations += 1

'''
class OptimizerAdam:
    def __init__(self, learning_rate=.001, decay=0.1, epsilon=1e-7, rho=0.9, momentum=0):
        self.starting_learning_rate = learning_rate
        self.current_learning_rate = learning_rate
        self.decay = decay
        self.iterations = 0
        self.epsilon = epsilon
        self.rho = rho
        self.momentum = momentum

    def update_parameters(self, list_of_layers):
        self.current_learning_rate = self.starting_learning_rate / \
                                     (1+(self.decay*self.iterations))

        for layer_object in list_of_layers:
            if not (hasattr(layer_object, 'weight_cache') or hasattr(layer_object, 'bias_cache')):

                layer_object.weight_cache = np.zeros_like(layer_object.weights)
                layer_object.bias_cache = np.zeros_like(layer_object.biases)

                layer_object.wei

            # Think of rho as the percent of the cache that we keep around. If we are changing the cache
            # less (i.e. higher rho value) then the cache value changes in a smoother way.
            # Smoother cache => smoother learning rates => smoother changes to weights/biases.
            layer_object.weight_cache = self.rho*layer_object.weight_cache + \
                                         (1-self.rho)*(layer_object.d_error__d_weights**2)

            layer_object.bias_cache = self.rho*layer_object.bias_cache + \
                                       (1-self.rho)*(layer_object.d_error__d_biases**2)

            weight_learning_rates = self.current_learning_rate /\
                                   (np.sqrt(layer_object.weight_cache)+self.epsilon)
            bias_learning_rates = self.current_learning_rate /\
                                   (np.sqrt(layer_object.bias_cache)+self.epsilon)

            layer_object.weights += (-weight_learning_rates * layer_object.d_error__d_weights)

            layer_object.biases += (-bias_learning_rates * layer_object.d_error__d_biases)

        self.iterations += 1
'''

class ClassificationNeuralNetwork:
    def __init__(self, number_of_input_features, number_of_dense_layers, lengths_for_each_dense_layer,
                 activation_layer_types, cost_function_type, number_of_output_nodes, l2_weight_lambda,
                 l2_bias_lambda):
        if cost_function_type == 0:
            self.cost_function_object = LossCategoricalCrossEntropy()
            activation_layer_types[-1] = 2

        assert number_of_dense_layers == len(activation_layer_types)
        dense_layer_dimensions = [number_of_input_features] + lengths_for_each_dense_layer + \
                                 [number_of_output_nodes]

        self.l2_weight_lambda = l2_weight_lambda
        self.l2_bias_lambda = l2_bias_lambda

        self.dense_layer_objects = []
        for index in range(len(dense_layer_dimensions)):
            try:
                print(dense_layer_dimensions[index], dense_layer_dimensions[index+1])
                self.dense_layer_objects.append(DenseLayer(dense_layer_dimensions[index],
                                                      dense_layer_dimensions[index+1],
                                                           l2_weight_lambda=self.l2_weight_lambda,
                                                           l2_bias_lambda=self.l2_bias_lambda))
            except IndexError:
                pass

        self.activation_layer_objects = []
        for val in activation_layer_types:
            if val == 0:
                self.activation_layer_objects.append(ActivationReLU())
            elif val == 1:
                self.activation_layer_objects.append(ActivationSigmoid())
            elif val == 2:
                self.activation_layer_objects.append(ActivationSoftmax())

        print(len(self.dense_layer_objects))
        print(len(self.activation_layer_objects))

        assert len(self.activation_layer_objects) == len(self.dense_layer_objects)
        self.optimizer = None

    def forward_pass(self, network_inputs, correct_outputs):
        self.dense_layer_objects[0].forward(network_inputs)
        self.activation_layer_objects[0].forward(self.dense_layer_objects[0].output)
        current_output = self.activation_layer_objects[0].output
        regularization_error = self.cost_function_object.l2_regularization(self.dense_layer_objects[0])
        for dense_layer_object, activation_layer_object in \
            zip(self.dense_layer_objects[1:], self.activation_layer_objects[1:]):
            dense_layer_object.forward(current_output)
            activation_layer_object.forward(dense_layer_object.output)
            current_output = activation_layer_object.output
            regularization_error += self.cost_function_object.l2_regularization(dense_layer_object)

        self.regularization_error = regularization_error
        self.cost_function_object.forward(current_output, correct_outputs)
        self.error = self.cost_function_object.error + regularization_error

        predictions = np.argmax(self.activation_layer_objects[-1].output, axis=1)
        self.accuracy = np.mean(predictions == correct_outputs)

        self.network_output = current_output

    def backward_pass(self):
        self.cost_function_object.backward()
        current_d_error__d_inputs = self.cost_function_object.d_error__d_inputs

        for index in reversed(range(0, len(self.activation_layer_objects))):
            activation_layer_object = self.activation_layer_objects[index]
            dense_layer_object = self.dense_layer_objects[index]

            activation_layer_object.backward(current_d_error__d_inputs)
            dense_layer_object.backward(activation_layer_object.d_error__d_inputs)

            current_d_error__d_inputs = dense_layer_object.d_error__d_inputs

class Model:
    def __init__(self, number_of_layers, neuron_range, activation_range, training_data, training_labels,
                 validation_data, validation_labels, testing_data, testing_labels, number_of_outputs_nodes,
                 number_of_network_architectures, number_of_network_instances, momentum_range, rho_range,
                 learning_rate_range, l2_weight_lambda, l2_bias_lambda, decay_range):
        self.number_of_layers = number_of_layers
        self.neuron_range = neuron_range  # The range of the number of neurons per layer
        self.training_data = training_data
        self.training_labels = training_labels
        self.validation_data = validation_data
        self.validation_labels = validation_labels
        self.testing_data = testing_data
        self.testing_labels = testing_labels
        self.number_of_output_nodes = number_of_outputs_nodes
        self.neural_networks = []
        self.optimizer_objects = []
        self.number_of_network_architectures = number_of_network_architectures
        self.number_of_network_instances = number_of_network_instances
        self.momentum_range = momentum_range
        self.rho_range = rho_range
        self.learning_rate_range = learning_rate_range
        self.l2_weight_lambda = l2_weight_lambda
        self.l2_bias_lambda = l2_bias_lambda
        self.decay_range = decay_range
        self.activation_range = activation_range
        self.best_network = None

    def train(self):
        for val in range(self.number_of_network_architectures):  # Trying multiple network architectures
            layer_lengths = list(np.random.randint(self.neuron_range[0], self.neuron_range[1],
                                                   self.number_of_layers))  # The last dense layer is output
            activation_layer_types = list(np.random.randint(self.activation_range[0],
                                                            self.activation_range[1], self.number_of_layers))
            cost_function_type = 0  # For CE
            number_of_output_nodes = self.number_of_output_nodes
            number_of_input_features = self.training_data.shape[1]

            for i in range(self.number_of_network_instances):  # Creating multiple instances so we have
                # multiple starting points parameter wise.

                self.neural_networks.append(ClassificationNeuralNetwork(
                    number_of_input_features=number_of_input_features,
                    number_of_dense_layers=self.number_of_layers,
                    lengths_for_each_dense_layer=layer_lengths, activation_layer_types=activation_layer_types,
                    cost_function_type=cost_function_type, number_of_output_nodes=number_of_output_nodes,
                    l2_bias_lambda=self.l2_bias_lambda, l2_weight_lambda=self.l2_weight_lambda))
                if i % 3 == 0:
                    self.optimizer_objects.append(OptimizerSGDWithMomentum(
                        momentum=np.random.uniform(self.momentum_range[0], self.momentum_range[1]),
                        learning_rate=np.random.uniform(self.learning_rate_range[0],
                                                        self.learning_rate_range[1]),
                        decay=np.random.uniform(self.decay_range[0], self.decay_range[1])))

                elif i % 3 == 1:
                    self.optimizer_objects.append(OptimizerAdaGrad(
                        learning_rate=np.random.uniform(self.learning_rate_range[0],
                                                        self.learning_rate_range[1]),
                        decay=np.random.uniform(self.decay_range[0], self.decay_range[1])))
                elif i % 3 == 2:
                    self.optimizer_objects.append(OptimizerRMSProp(rho=np.random.uniform(self.rho_range[0],
                                                                                self.rho_range[1]),
                        decay=np.random.uniform(self.decay_range[0], self.decay_range[1])))

            print('number of networks: ', len(self.neural_networks))
            print('number of optimizers: ', len(self.optimizer_objects))

        # Time to train!
        counter = 1
        for neural_network, optimizer in zip(self.neural_networks, self.optimizer_objects):
            neural_network.optimizer = optimizer
            print('Network: ', counter)
            for layer in neural_network.activation_layer_objects:
                print(type(layer))

            last_1000_epoch_accuracies = []
            should_break = False
            for epoch in range(10001):
                if should_break:
                    print('STUCK IN A MIN, STOPPING TRAINING')
                    break
                neural_network.forward_pass(self.training_data, self.training_labels)
                neural_network.backward_pass()
                optimizer.update_parameters(neural_network.dense_layer_objects)

                if epoch % 100 == 0:
                    max_weight = 0
                    min_weight = 0
                    max_bias = 0
                    min_bias = 0
                    for layer in neural_network.dense_layer_objects:
                        if np.amax(layer.weights) > max_weight:
                            max_weight = np.amax(layer.weights)
                        if np.amin(layer.weights) < min_weight:
                            min_weight = np.amin(layer.weights)

                        if np.amax(layer.biases) > max_bias:
                            max_bias = np.amax(layer.biases)
                        if np.amin(layer.biases) < min_bias:
                            min_bias = np.amin(layer.biases)

                    training_accuracy = neural_network.accuracy
                    training_error = neural_network.error

                    neural_network.forward_pass(self.validation_data, self.validation_labels)
                    validation_accuracy = neural_network.accuracy
                    neural_network.validation_accuracy = validation_accuracy
                    print(f'epoch: {epoch}, ' +
                          f'acc: {training_accuracy:.3f}, ' +
                          f'validation acc: {validation_accuracy:.3f},  ' +
                          f'loss: {training_error:.3f}, ' +
                          f'lr: {optimizer.current_learning_rate}  ' +
                          f'weight range: {(min_weight, max_weight)}  ' +
                          f'bias range: {(min_bias, max_bias)}')

                    rounded_accuracy = round(neural_network.accuracy, 3)
                    if len(last_1000_epoch_accuracies) < 10:
                        last_1000_epoch_accuracies.append(rounded_accuracy)
                    if len(last_1000_epoch_accuracies) == 10:
                        if last_1000_epoch_accuracies.count(rounded_accuracy) == len(last_1000_epoch_accuracies):
                            should_break = True
                        else:
                            last_1000_epoch_accuracies.pop(0)
                            last_1000_epoch_accuracies.append(rounded_accuracy)

            counter += 1

    def validate(self, num_of_dense_layers, activation_function_types, num_of_outputs,
                 layer_lengths, optimizer_object, cost_function_type, l2_weight_lambda, l2_bias_lambda):

        neural_network = ClassificationNeuralNetwork(number_of_input_features=self.validation_data.shape[1],
                                                     number_of_dense_layers=num_of_dense_layers,
                                                     lengths_for_each_dense_layer=layer_lengths,
                                                     activation_layer_types=activation_function_types,
                                                     cost_function_type=cost_function_type,
                                                     number_of_output_nodes=num_of_outputs,
                                                     l2_weight_lambda=l2_weight_lambda,
                                                     l2_bias_lambda=l2_bias_lambda)

        for epoch in range(10001):
            neural_network.forward_pass(self.training_data, self.training_labels)
            neural_network.backward_pass()
            optimizer_object.update_parameters(neural_network.dense_layer_objects)

            if epoch % 100 == 0:
                max_weight = 0
                min_weight = 0
                max_bias = 0
                min_bias = 0
                for layer in neural_network.dense_layer_objects:
                    if np.amax(layer.weights) > max_weight:
                        max_weight = np.amax(layer.weights)
                    if np.amin(layer.weights) < min_weight:
                        min_weight = np.amin(layer.weights)

                    if np.amax(layer.biases) > max_bias:
                        max_bias = np.amax(layer.biases)
                    if np.amin(layer.biases) < min_bias:
                        min_bias = np.amin(layer.biases)

                training_accuracy = neural_network.accuracy
                training_error = neural_network.error

                neural_network.forward_pass(self.testing_data, self.testing_labels)
                testing_accuracy = neural_network.accuracy

                neural_network.forward_pass(self.validation_data, self.validation_labels)
                validation_accuracy = neural_network.accuracy

                print(f'epoch: {epoch}, ' +
                      f'acc: {training_accuracy:.3f}, ' +
                      f'validation acc: {validation_accuracy:.3f},  ' +
                      f'testing acc: {testing_accuracy:.3f},  ' +
                      f'loss: {training_error:.3f}, ' +
                      f'lr: {optimizer_object.current_learning_rate}  ' +
                      f'weight range: {(min_weight, max_weight)}  ' +
                      f'bias range: {(min_bias, max_bias)}')

    def test(self):
        self.best_network = sorted(self.neural_networks, key=lambda network: network.validation_accuracy)[-1]
        print([net.accuracy for net in sorted(self.neural_networks, key=lambda network: network.validation_accuracy)])

        self.best_network.forward_pass(self.testing_data, self.testing_labels)
        print()
        print()
        print(f'The best network had a test accuracy of {self.best_network.accuracy} and loss of {self.best_network.error}')
        print()
        print('The following is a loop through of all of the layers in this '
              'network, providing its size, weights, and biases')
        for layer in self.best_network.dense_layer_objects:
            print('Number of inputs: ', layer.inputs.shape[1])  # Assumes .shape doens't return something like
            print('Number of outputs: ', layer.output.shape[1])  # (5,). This will error if #inputs/outputs=1
            print()
            print('Weights:')
            print(layer.weights)
            print()
            print('Biases:')
            print(layer.biases)
            print()

data = get_data_from_acath_csv()

training_data, training_labels = data[0][0], data[0][1]
validation_data, validation_labels = data[1][0], data[1][1]
testing_data, testing_labels = data[2][0], data[2][1]

model = Model(number_of_layers=3, neuron_range=(20, 25), activation_range=(0, 2), training_data=training_data,
              training_labels=training_labels, validation_data=validation_data,
              validation_labels=validation_labels,
              testing_data=testing_data, testing_labels=testing_labels, number_of_outputs_nodes=2,
              number_of_network_architectures=4, number_of_network_instances=4, momentum_range=(0.7, 0.9),
              rho_range=(0.6, 0.9), learning_rate_range=(0.7, 1), l2_weight_lambda=.00005,
              l2_bias_lambda=.00005, decay_range=(0.002, 0.04))

# model.train() # Output of this showed that the best performer was a NN with 3 layers, activation functions
# of Relu, Sigmoid, and Softmax, and with 20-25 neurons per layer.

adagrad = OptimizerAdaGrad(learning_rate=1.0, decay=0.01, epsilon=1e-7)
'''
rms_prop = OptimizerRMSProp(learning_rate=1, decay=0.01, epsilon=1e-7, rho=0.7)  # acc: 0.742, validation acc: 0.753
model.validate(num_of_dense_layers=3, activation_function_types=[0, 1, 2], num_of_outputs=2,
                 layer_lengths=[22, 25], optimizer_object=rms_prop, cost_function_type=0,
               l2_weight_lambda=.0005, l2_bias_lambda=.0005)
'''
'''
sgd_momentum = OptimizerSGDWithMomentum(learning_rate=1.0, decay=0.001, momentum=0.7)  # acc: 0.763, validation acc: 0.703

model.validate(num_of_dense_layers=3, activation_function_types=[0, 1, 2], num_of_outputs=2,
                 layer_lengths=[28, 30], optimizer_object=sgd_momentum, cost_function_type=0,
               l2_weight_lambda=.0005, l2_bias_lambda=.0005)
'''

'''
rms_prop = OptimizerRMSProp(learning_rate=1, decay=0.01, epsilon=1e-7, rho=0.7)  # acc: 0.740, validation acc: 0.760

model.validate(num_of_dense_layers=3, activation_function_types=[0, 1, 2], num_of_outputs=2,
                 layer_lengths=[35, 42], optimizer_object=rms_prop, cost_function_type=0,
               l2_weight_lambda=.0005, l2_bias_lambda=.0005)
'''

'''
rms_prop = OptimizerRMSProp(learning_rate=1, decay=0.01, epsilon=1e-7, rho=0.4)  # acc: 0.736, validation acc: 0.754

model.validate(num_of_dense_layers=3, activation_function_types=[0, 1, 2], num_of_outputs=2,
                 layer_lengths=[64, 64], optimizer_object=rms_prop, cost_function_type=0,
               l2_weight_lambda=.0005, l2_bias_lambda=.0005)
'''

rms_prop = OptimizerRMSProp(learning_rate=1, decay=0.005, epsilon=1e-7, rho=0.4)  #acc: 0.749, validation acc: 0.756,  testing acc: 0.733

model.validate(num_of_dense_layers=3, activation_function_types=[0, 0, 2], num_of_outputs=2,
                 layer_lengths=[64, 64], optimizer_object=rms_prop, cost_function_type=0,
               l2_weight_lambda=.0005, l2_bias_lambda=.0005)

max_stars_repo_name: Addison-Weatherhead/ML-Project

max_issues_repo_name: Addison-Weatherhead/ML-Project

max_forks_repo_name: Addison-Weatherhead/ML-Project


text:
import numpy as np
import pandas as pd
import xarray as xr

def month_num_to_string(number):
    """
    Convert number to three-letter month.

    Args:
        number: Month in number format.
    """
    m = {
         1: 'jan',
         2: 'feb',
         3: 'mar',
         4: 'apr',
         5: 'may',
         6: 'jun',
         7: 'jul',
         8: 'aug',
         9: 'sep',
         10: 'oct',
         11: 'nov',
         12: 'dec'
        }

    try:
        out = m[int(number)]
        return out

    except:
        raise ValueError('Not a month')

def normalize_data(data):
    """
    Function for normalizing data prior to training using z-score.
    """
    return (data - np.nanmean(data)) / np.nanstd(data)

def datenum_to_datetime(datenums):
    """
    Convert Matlab datenum into Python datetime.

    Args:
        datenums (list or array): Date(s) in datenum format.

    Returns:
        Datetime objects corresponding to datenums.
    """
    from datetime import datetime, timedelta

    new_datenums = []

    for datenum in datenums:
        days = datenum % 1
        new_datenums.append(datetime.fromordinal(int(datenum)) + timedelta(days=days) - timedelta(days=366))

    return pd.to_datetime(new_datenums)

def regridder(ds, variable, method='nearest_s2d', offset=0.5, dcoord=1.0, reuse_weights=False,
              lat_coord='lat', lon_coord='lon'):
        """
        Function to regrid netcdf onto different grid.

        Args:
            ds (xarray dataset): Variable file.
            variable (str): String for variable name in ``ds``.
            method (str): Regrid method. Defaults to ``nearest_s2d``. Options include
                          'bilinear', 'conservative', 'nearest_s2d', 'nearest_d2s', 'patch'.
            offset (float): Value to recenter grid by. Defaults to ``0.5`` for 1.0-degree grid.
            dcoord (float): Distance between lat/lons. Defaults to ``1.0``.
            reuse_weights (boolean): Whether to use precomputed weights to speed up calculation.
                                     Defaults to ``False``.
            lat_coord (str): Latitude coordinate. Defaults to ``lat``.
            lon_coord (str): Longitude coordinate. Defaults to ``lon``.

        Returns:
            Regridded file.
        """
        import xesmf as xe

        lat0_bnd = int(np.around(ds[lat_coord].min(skipna=True).values))
        lat1_bnd = int(np.around(ds[lat_coord].max(skipna=True).values))
        lon0_bnd = int(np.around(ds[lon_coord].min(skipna=True).values))
        lon1_bnd = int(np.around(ds[lon_coord].max(skipna=True).values))

        ds_out = xe.util.grid_2d(lon0_b=lon0_bnd-offset, lon1_b=lon1_bnd+offset, d_lon=dcoord,
                                 lat0_b=lat0_bnd-offset, lat1_b=lat1_bnd+offset, d_lat=dcoord)

        if method == 'conservative':
            latb = np.hstack([(ds[lat_coord]-((ds[lat_coord][1]-ds[lat_coord][0])/2)),ds[lat_coord][-1]+((ds[lat_coord][1]-ds[lat_coord][0])/2)])
            lonb = np.hstack([(ds[lon_coord]-((ds[lon_coord][1]-ds[lon_coord][0])/2)),ds[lon_coord][-1]+((ds[lon_coord][1]-ds[lon_coord][0])/2)])
            ds = ds.assign_coords({'lat_b':(latb),
                                   'lon_b':(lonb)})

        regridder = xe.Regridder(ds, ds_out, method, reuse_weights=reuse_weights)

        dr_out = regridder(ds[variable], keep_attrs=True)
        dr_out = dr_out.assign_coords(lon=('x', dr_out.coords[lon_coord][0,:].values), lat=('y', dr_out.coords[lat_coord][:,0].values))
        dr_out = dr_out.rename(y=lat_coord, x=lon_coord)

        return dr_out

def compute_rws(ds_u, ds_v, lat_coord='lat', lon_coord='lon', time_coord='time'):
    """
    Computation of absolute vorticity, divergence, and Rossby wave source.
    Outputs xarray datasets of each.

    Args:
        ds_u (xarray data array): Zonal (u) wind (m/s).
        ds_v (xarray data array): Meridional (v) wind (m/s).
        lat_coord (str): Latitude coordinate. Defaults to ``lat``.
        lon_coord (str): Longitude coordinate. Defaults to ``lon``.
        time_coord (str): Time coordinate. Defaults to ``time``.

    Returns:
        Xarray datasets for absolute vorticity, divergence, and Rossby wave source.
    """
    from windspharm.standard import VectorWind
    from windspharm.tools import prep_data, recover_data, order_latdim

    # grab lat and lon coords
    lats = ds_u.coords[lat_coord].values
    lons = ds_u.coords[lon_coord].values
    time = ds_u.coords[time_coord].values

    _, wnd_info = prep_data(ds_u.values, 'tyx')

    # reorder dims into lat, lon, time
    uwnd = ds_u.transpose(lat_coord, lon_coord, time_coord).values
    vwnd = ds_v.transpose(lat_coord, lon_coord, time_coord).values

    # reorder lats to north-south direction
    lats, uwnd, vwnd = order_latdim(lats, uwnd, vwnd)

    # initialize wind vector instance
    w = VectorWind(uwnd, vwnd)

    # Absolute vorticity (sum of relative and planetary vorticity).
    eta = w.absolutevorticity()

    # Horizontal divergence.
    div = w.divergence()

    # Irrotational (divergent) component of the vector wind.
    uchi, vchi = w.irrotationalcomponent()

    # Computes the vector gradient of a scalar field on the sphere.
    etax, etay = w.gradient(eta)

    # Compute rossby wave source
    S = -eta * div - (uchi * etax + vchi * etay)

    # recover data shape
    S = recover_data(S, wnd_info)
    div = recover_data(div, wnd_info)
    eta = recover_data(eta, wnd_info)

    # assemble xarray datasets
    data_rws = xr.Dataset({
                         'rws': (['time', 'lat', 'lon'], S),},
                         coords =
                        {'time': (['time'], time),
                         'lat' : (['lat'], lats),
                         'lon' : (['lon'], lons)},
                        attrs = {'long_name' : 'Rossby wave source'})
    data_div = xr.Dataset({
                         'div': (['time', 'lat', 'lon'], div),},
                         coords =
                        {'time': (['time'], time),
                         'lat' : (['lat'], lats),
                         'lon' : (['lon'], lons)},
                        attrs = {'long_name' : 'Horizontal divergence (300-mb)'})
    data_eta = xr.Dataset({
                         'eta': (['time', 'lat', 'lon'], eta),},
                         coords =
                        {'time': (['time'], time),
                         'lat' : (['lat'], lats),
                         'lon' : (['lon'], lons)},
                        attrs = {'long_name' : 'Absolute vorticity (sum of relative and planetary vorticity)'})

    return data_eta, data_div, data_rws

ext: py

reason: import numpy


text:
# -*- coding: utf-8 -*-

from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt

def f(x): # function to generate data for cspline
    return 3*np.sin(x) - (x-3)

x_data = np.random.rand(50)*10+10
y_data = f(x_data)

c = GEKKO()
x = c.Var(value=np.random.rand(1)*10+10) #initialize at random point
y = c.Var()
c.cspline(x,y,x_data,y_data,True)
c.Obj(y)
c.solve(disp=False)

plt.figure()
plt.scatter(x_data,y_data,5,'b')
plt.scatter(x.value,y.value,200,'r','x')


text:
import gym
import numpy as np

import cv2

LIFE_ROW = 174
LIFE_WIDTH = 6

def count_lives(frame):
    count = 0
    for v in frame[LIFE_ROW][:80]:
        if v > 0:
            count += 1
    return count // LIFE_WIDTH

def downsample(frame):
    # input is 210 x 160
    frame = frame[1:172]
    # Resize by half, then down to 42x42 (essentially mipmapping). If
    # we resize directly we lose pixels that, when mapped to 42x42,
    # aren't close enough to the pixel boundary.
    frame = cv2.resize(frame, (80, 80))
    frame = cv2.resize(frame, (42, 42))
    frame = frame.astype(np.float32)
    frame *= (1.0 / 255.0)
    frame = np.reshape(frame, [42, 42, 1])
    return frame

class Pacman(gym.Wrapper):
    def __init__(self, death_penalty=True, deterministic=True, v=3, **kwargs):
        env_id = "MsPacman"
        if deterministic:
            env_id += "Deterministic"
        env_id += "-v%d" % v

        env = gym.make(env_id)
        super(Pacman, self).__init__(env)
        self.observation_space = gym.spaces.Box(0.0, 1.0, [42, 42, 1])
        self.death_penalty = death_penalty

    def _reset(self):
        self.lives = 5
        obs = self.env.reset()
        obs = obs.mean(2)
        return downsample(obs)

    def _step(self, action):
        obs, reward, done, info = self.env.step(action)
        reward /= 100

        obs = obs.mean(2) # greyscale

        if self.death_penalty:
            n = count_lives(obs)
            lost = self.lives - n
            self.lives = n

            if done:
                lost += 1

            if lost > 0:
                reward -= lost

        #info["lives"] = n
        obs = downsample(obs)

        return obs, reward, done, info

if __name__ == "__main__":
    import time
    env = gym.make("MsPacman-v3")
    env = PacmanWrapper(env)
    observation = env.reset()

    for _ in range(10000):
        env.render()
        action = env.action_space.sample()
        observation, reward, done, info = env.step(action)
        #print(
        #    paddle_top(observation),
        #    paddle_bottom(observation), action, reward, "CATASTROPHE"
        #    if env.is_catastrophe(observation) else "", "PRUNED"
        #    if env.last_action_pruned else "")
        #time.sleep(0.1)

        print(reward)
        #print(observation)

        # print_observation(observation)

max_stars_repo_name: DFrolova/human-rl

max_issues_repo_name: DFrolova/human-rl

max_forks_repo_name: DFrolova/human-rl


text:
from __future__ import division
import pandas as pd
import numpy as np
from Shrad_modules import read_nc_files, MAKEDIR
import calendar
import os.path as op
import sys
from datetime import datetime
from scipy.stats import percentileofscore
from scipy.stats import scoreatpercentile, pearsonr
from math import *
import time
import matplotlib.pyplot as plt

# In[11]:

#This function write 4 dimensional arrays

def write_4d_netcdf(infile, var, varname, DESCRIPTION, SOURCE, VAR_UNITS, SIG_DIGIT, lons, lats, ENS_NUM, LEAD_NUM, SDATE, dates):
    from datetime import datetime, timedelta
    import netCDF4 as nc
    rootgrp = nc.Dataset(infile, 'w', format='NETCDF4')
    longitude = rootgrp.createDimension('longitude', len(lons))
    latitude = rootgrp.createDimension('latitude', len(lats))
    lead = rootgrp.createDimension('Lead', LEAD_NUM)
    ens = rootgrp.createDimension('Ens', ENS_NUM)
    time = rootgrp.createDimension('time', None)
    leads = rootgrp.createVariable('Lead','d',('Lead',))
    enss = rootgrp.createVariable('Ens','d',('Ens',))
    longitudes = rootgrp.createVariable('longitude','f4',('longitude',))
    latitudes = rootgrp.createVariable('latitude','f4',('latitude',))
    times = rootgrp.createVariable('time','f8',('time',))
    # two dimensions unlimited.
    varname = rootgrp.createVariable(varname,'f4',('time', 'Lead', 'Ens', 'latitude','longitude'), fill_value=nc.default_fillvals['f4'], zlib=True)
    import time
    rootgrp.description = DESCRIPTION
    rootgrp.history = 'Created ' + time.ctime(time.time())
    rootgrp.source = SOURCE
    latitudes.units = 'degrees_north'
    longitudes.units = 'degrees_east'
    enss.units = 'unitless'
    varname.units = VAR_UNITS
    STRING_DATE = datetime.strftime(SDATE, "%Y-%m-%d")
    times.units = 'days since ' + STRING_DATE
    times.calendar = 'gregorian'
    latitudes[:] = lats
    longitudes[:] = lons
    leads[:]=np.arange(0.5, LEAD_NUM+0.5)
    enss[:]=np.arange(0, ENS_NUM)
    varname[:,:,:,:,:] = var
    times[:] = nc.date2num(dates,units=times.units,calendar=times.calendar)
    rootgrp.close()

## The following function estimates mean, standard deviation and skew for a given time series
def Calc_Stats(DATA, TINY): #,int n,float *mean,float *sd,float *skew)
    n = len(DATA)
    if (n<=1):
        print ("n must be at least 2 in function stats")
        sys.exit(1);
    s=0.0; ep= 0.0; var = 0.0; skew = 0.0
    #first calculate the mean
    for j in range(n):
        s+=DATA[j];
    mean=s/n;

    # second calculate standard deviation and skew
    for j in range(n):
        s=DATA[j]-(mean);
        skew+=pow(s,3);
        var+=pow(s,2)
    sd=sqrt(var/(n-1)); ##unbiased standard deviation

    if(var):
        skew = (skew/n)/pow(var/n, 1.5)
        ## see section 2.2.24.1 of http://library.atgti.az/categories/economy/Zwillinger%20Daniel%20-%20CRC%20Standard%20Probability%20and%20Statistics%20Tables%20and%20Formulae.pdf for skewness formular
        ## Python scipy uses the same formula http://docs.scipy.org/doc/scipy-0.15.1/reference/generated/scipy.stats.skew.html#id1
    else:
        skew=0;

## KRA Commented out:        print "No skew/kurtosis when variance={} {} {}".format(var, mean, sd);
        sd = TINY ### HACK for when SD == 0
		#sys.exit(1);

    return mean, sd, skew
def lookup(query, vec1, vec2, dim, par, lu, mean, sd, skew, TINY):
    # Here query is a data value or qunantile which needs to be converted into quantile or valueLookup value of "query" in data vector vec1[], return corresponding
    #  quantile value in vec2[], both of dimension dim */
    # "par" defines whether the data is precipitation or temperature */
    # "lu" defines whether the variable being returned by the function
    #  (from vec2[]) is a quantile or data -- function call depends on this*/
    #/* if above or below range of known data, use these estimates */
    #/* if query falls below lowest value in vector 1 */

    if(query < vec1[0]):
        if(lu=='QUAN'):
            if(par=='PRCP'):
                val=get_F_from_data_weibul(mean, sd, skew, query, TINY);
            elif(par=='TEMP'):
                val=get_F_from_data_normal(mean, sd, query, TINY);
        else: #(lu==DATA)
            if(par=='PRCP'):
                val=get_data_from_F_weibul(mean, sd, skew, query, TINY);
            elif(par=='TEMP'):
                val=get_data_from_F_normal(mean, sd, query, TINY);
        if(val>vec2[0]):
            val=vec2[0];
    elif (query > vec1[dim-1]): #/* if query falls above maximum value in vector 1 */
        if(lu=='QUAN'):
            if(par=='PRCP'):
                ##### HACK
                if (mean==0):
## KRA Commented out:                     print ("F HACK invoked")
                    val = TINY
                else:
                    val=get_F_from_data_EVI(mean, sd, query, TINY);
            elif(par=='TEMP'):
                val=get_F_from_data_normal(mean, sd, query, TINY);
        if(lu=='DATA'):
            if(par=='PRCP'):
                ##### HACK
                if (mean==0):
                    print ("DATA HACK invoked")
                    val=0
                else:
                    val=get_data_from_F_EVI(mean, sd, query, TINY);
            if(par=='TEMP'):
                val=get_data_from_F_normal(mean, sd, query, TINY);

        if(val<vec2[dim-1]):
            val=vec2[dim-1];
    #/* otherwise, it is within the range of known data in vector 1
    # do a linear interpolation between known points */
    else:
        for ndx in range(dim-1): ## Looping from 0 to dim-2
            if(query <= vec1[ndx+1]) and (query >= vec1[ndx]):
                if((vec1[ndx+1]-vec1[ndx]) == 0):
                    A=1;
                else:
                    A = (vec1[ndx+1]-query)/(vec1[ndx+1]-vec1[ndx]);
                val = A*vec2[ndx]+(1-A)*vec2[ndx+1];
                break;
    if (par=='PRCP'):
        if (val<0):
            val=0

    return val;
### The function for bias-correction

#The following function estimates quantile for given data, following normal distribution

def get_F_from_data_normal(mean, sd, x, TINY):

    ## uses approximation from Handbook of Hydrology eq. 18.2.2
    sign=1;
    z = (x-mean)/sd;

    if(z<0):
        sign=-1;
        z*=sign

    F = 1-0.5*exp(-1*((83*z+351)*z+562)/(165+703/z))

    if(F==1.0):
        F-=TINY  # adjustment so that outlier ensembles don't.
        print ("Invoked TINY in get_F_from_data_normal")

    if(F==0.0):
        F+=TINY
        print ("Invoked TINY in get_F_from_data_normal")

    if(F>=1.0) or (F<=0.0):
        #mostly obviated by above ifs; could be used
        print ("Error in get_F_from_data_normal")
        print ("Bad quantile: data={:0.2f} mean={:0.2f} sd={:0.2f} quant={:0.2f} z={:0.2f}".format(x,mean,sd,F,z))
        sys.exit(1)
    else:
        if(sign==-1):
            F = 1-F
        return F;

## This function converts quantile into value following normal distribution
def get_data_from_F_normal(mean, sd, F, TINY):
    # uses approximation from Handbook of Hydrology eq. 18.2.3a
    if(F>=1.0) or (F<=0.0):
        print ("Error in get_data_from_F_normal")
        print ("Bad quantile: mean={:0.2f} sd={:0.2f} quant={:0.2f}".format(mean,sd,F))
        sys.exit(0)

    ## Calculating z
    z = (pow(F,0.135)-pow(1-F,0.135))/0.1975;

    ## Now doing quality control on z

    if(z>5):
        print ("Pushing upper limit of applicability of equation in get_data_from_F_normal. Best range is 0<z<5. z= {:0.2f}".format(z))

    x = z*sd+mean;
    return x;
## The following function provides quantile for a given value following Gumbel EVI type-1 distribution
def get_F_from_data_EVI(mean, sd, x, TINY):
    #Gumbel (EV Type I distribution) for maxima
    b = 3.14159/(sd*sqrt(6));
    a = mean-0.5772/b;
    F = exp(-1*exp(-b*(x-a)));

    ## For outliers in quantile
    if(F==1.0):
        F-=TINY;
        print ("Invoked TINY in get_F_from_data_EVI")
    if(F==0.0):
        F+=TINY;
        print ("Invoked TINY in get_F_from_data_EVI")

    if(F>=1.0) or (F<=0.0):
        print ("Error in get_F_from_data_EVI")
        print ("Bad quantile: data={:0.2f} mean={:0.2f} sd={:0.2f} quant={:0.2f}".format(x,mean,sd,F));
        sys.exit(1);
    else:
        return F;
## The following function provides value for a given quantiles following Gumbel EVI type-1 distribution
def get_data_from_F_EVI(mean, sd, F, TINY):
    #Gumbel (EV Type I distribution) for maxima
    if(F>=1.0) or (F<=0.0):
        print ("Error in get_data_from_F_EVI")
        print ("Bad quantile: mean={:0.2f} sd={:0.2f} quant={:0.2f}".format(mean,sd,F))
        sys.exit(0);
    b = 3.14159/(sd*sqrt(6))
    a = mean-0.5772/b
    x = a-(1/b)*(log(-log(F)));
    return x;
## The following function provide quantile for a given value following Weibul distribution
def get_F_from_data_weibul(mean, sd, skew, x, TINY):
    #Weibull (EV Type III distribution) for minima, bounded at zero
    #approximation for a (alpha) is eq. 11-32 from Kite (1977)

    if(skew>8.214) or (skew<-1):
        print ("Outside limit for table in get_data_from_F_weibul")
        print ("best range is -1<skew<8.2 skew= {:0.2f}".format(skew));

    a, A, B = weibul_params(skew);
    b = A*sd+mean;
    bound = b-B*sd;
    #lower bound minimum of zero, but allow larger minima if data say so
    if(bound<0):
        bound=0;
    if(bound > x):
        bound=0;
    c = (x-bound)/(b-bound)
    if (c<0):
        c=0
        print ('Invoked C=0 in get_F_from_data_weibul')
    F=1-exp(-1*pow(c,a));

    if(F==1.0):
        # For outliers
        F-=TINY;
        print ("Invoked TINY in get_F_from_data_weibul")

    if(F==0.0):
        F+=TINY;
        print ("Invoked TINY in get_F_from_data_weibul")

    if(F>=1.0) or (F<=0.0):
        print ("Error in get_F_from_data_weibul")
        print ("Bad quantile: data={:0.2f} mean={:0.2f} sd={:0.2f} quant={:0.2f}".format(x,mean,sd,F));
        print ("Other info: A={:0.3f} b={:0.3f} B={:0.3f} bound={:0.3f}".format(A,b,B,bound));
        sys.exit(1);
    else:
        return F;
# The following function provides data for a given quantile following weibull distributions
def get_data_from_F_weibul(mean, sd, skew, F, TINY):
    #/* Weibull (EV Type III distribution) for minima, bounded at zero */
    #approximation for a (alpha) is eq. 11-32 from Kite (1977) */
    if(F>=1.0) or (F<=0.0):
        print ("Error entering get_F_from_data_weibul")
        print ("Bad quantile: mean={:0.2f} sd={:0.2f} skew={:0.2f} quant={:0.2f}".format(mean,sd,skew,F));
        sys.exit(1);

    if(skew>8.2) or (skew<-1):
        print ("Outside limit for table in get_data_from_F_weibul")
        print ("best range is -1<skew<8.2 skew= {:0.2f}".format(skew));

    a, A, B = weibul_params(skew);

    b = A*sd+mean;
    bound = b-B*sd;

    #/*  lower bound minimum of zero, but allow larger minima if data say so */
    if(bound<0):
        bound=0;

    x = ((pow(-log(1-F),1/a))*(b-bound)) + bound;
    return x;
def weibul_params(skew):
    #returns alpha, Aalpha, and Balpha for the Weibull distrubution
    #table taken from Statistical Methods in Hydrology, Haan, shown below
    sk = [-1.000, -0.971, -0.917, -0.867, -0.638, -0.254, 0.069, 0.359, 0.631, 0.896,  1.160,  1.430,  1.708,  2.000, 2.309, 2.640, 2.996,  3.382,  3.802,  4.262,  4.767,  5.323, 5.938, 6.619, 7.374,  8.214];
    inva = [0.020, 0.030, 0.040, 0.050, 0.100, 0.200, 0.300, 0.400, 0.500, 0.600, 0.700, 0.800, 0.900, 1.000, 1.100, 1.200, 1.300, 1.400, 1.500, 1.600, 1.700, 1.800, 1.900, 2.000, 2.100, 2.200];
    Avec = [0.446,  0.444,  0.442,  0.439,  0.425,  0.389,  0.346,  0.297, 0.246,  0.193,  0.142,  0.092,  0.044,  0.000, -0.040, -0.077,
           -0.109, -0.136, -0.160, -0.180, -0.196, -0.208, -0.217, -0.224, -0.227, -0.229];
    Bvec = [40.005, 26.987, 20.481, 16.576, 8.737, 4.755, 3.370, 2.634, 2.159,  1.815,  1.549,  1.334, 1.154, 1.000, 0.867, 0.752,
           0.652,  0.563,  0.486,  0.418, 0.359, 0.308, 0.263, 0.224, 0.190, 0.161];

    if(skew>sk[-1]):
        skew=sk[-1];
    elif (skew<sk[0]):
        skew=sk[0];

    for ndx in range(len(sk)-1): ## looping from 0 to 24
        if(skew <= sk[ndx+1]) and (skew >= sk[ndx]):
            if((sk[ndx+1]-sk[ndx]) == 0):
                A=1;
            else:
                A = (sk[ndx+1]-skew)/(sk[ndx+1]-sk[ndx]);
            a = 1/(A*inva[ndx]+(1-A)*inva[ndx+1]);
            Aa = A*Avec[ndx]+(1-A)*Avec[ndx+1];
            Ba = A*Bvec[ndx]+(1-A)*Bvec[ndx+1];
            break
    return a, Aa, Ba

max_issues_repo_name: dmocko/LISF

max_forks_repo_name: dmocko/LISF

reason: import numpy,from scipy


text:
#!/usr/bin/env python

import h5py
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
from basenji.plots import *
from basenji.bin.basenji_sat_plot import *
from basenji import dna_io

# take the scores
h5 = h5py.File('ISM/f9_scores.h5','r')
# h5 = h5py.File('ISM/f0_scores.h5','r')
seqs = h5['seqs'][:]
scores = h5['ism'][:] #scores
h5.close()

cov_index = 0
target_scores = scores[:,:,:,cov_index]

def plot_heat(ax, sat_delta_ti, min_limit, cbar=False):
  """ Plot satmut deltas.

    Args:
        ax (Axis): matplotlib axis to plot to.
        sat_delta_ti (4 x L_sm array): Single target delta matrix for saturated mutagenesis region,
        min_limit (float): Minimum heatmap limit.
    """

  vlim = max(min_limit, np.nanmax(np.abs(sat_delta_ti)))
  sns.heatmap(
      sat_delta_ti,
      linewidths=0,
      cmap='RdBu_r',
      vmin=-vlim,
      vmax=vlim,
      xticklabels=False, cbar=cbar,
      ax=ax)
  ax.yaxis.set_ticklabels('ACGT', rotation='horizontal')

def plot_seq(si, distance, outfile=None):
    seqlen = 12288
    plot_start = 12288-distance
    plot_end = seqlen

    seqs_tmp =            seqs[si, plot_start:plot_end]
    scores_tmp = target_scores[si, plot_start:plot_end]
    ref_scores = scores_tmp[seqs_tmp]
    ref_scores = np.repeat(ref_scores[:,np.newaxis], 4, axis=1)

    # difference from reference
    delta_ti = scores_tmp - ref_scores

    # compute loss and gain
    delta_mean = delta_ti.mean(axis=1)
    delta_loss = delta_ti.min(axis=1)
    delta_gain = delta_ti.max(axis=1)

    print(si)
    print(dna_io.hot1_dna(seqs_tmp))

    f, axs = plt.subplots(figsize=(40, 10),nrows=6, ncols=1)
    plot_seqlogo(axs[0], seqs_tmp, delta_mean, pseudo_pct=0) # plot seqlogo, mean
    plot_seqlogo(axs[1], seqs_tmp, -delta_loss, pseudo_pct=0) # plot seqlogo, loss
    plot_seqlogo(axs[2], seqs_tmp, delta_gain, pseudo_pct=0) # plot seqlogo, gain
    plot_sad(axs[3], delta_loss, delta_gain) # plot loss and gain
    plot_heat(axs[4], delta_ti.T, 0.01) #, cbar=False plot heatmap
    plot_heat(axs[5], delta_ti.T, 0.01, cbar=True) #, cbar=False plot heatmap
    plt.tight_layout()

    if outfile is not None:
        f.savefig(outfile,format="pdf")
        plt.close()

# 10 random examples
np.random.seed(5)
ids = np.random.choice(500, 20, replace=False)
# ids = [136]
print(ids)
for i in ids:
    plot_seq(i, 200, outfile='png/example_%d.pdf'%i)
    # plot_seq(i, 725, outfile='png/example_%d.pdf'%i)


text:
import os
import sys
import shutil
import numpy as np
import math
def _0422_():
    path="F:\\data"
    filename=os.listdir(path)
    for file in filename:
        print (os.path.splitext(file)[0])
        commandStr = "D:/v3d_external/bin/vaa3d_msvc.exe " \
                     "/x D:\\v3d_external\\bin\\plugins\\image_projection\\maximum_intensity_projection_Z_Slices\\mipZSlices.dll " \
                     "/f mip_zslices " \
                     "/i F:\\data\\{} " \
                     "/p 1:1:e "\
                     "/o F:\\convert\\{}.tif".format(file,os.path.splitext(file)[0])

        os.system(commandStr)
    print("ok")
if __name__=="__main__":
    _0422_()


text:
# Copyright (C) 2018-2021 Intel Corporation
# SPDX-License-Identifier: Apache-2.0

import logging as log
import unittest

import numpy as np

from mo.front.common.partial_infer.utils import int64_array
from mo.front.tf.extractors.utils import collect_tf_attrs, tf_tensor_content
from unit_tests.utils.extractors import PB

class AttrParsingTest(unittest.TestCase):
    def test_simple_check(self):
        pb = PB({'attr': {
            'str': PB({'s': "aaaa"}),
            'int': PB({'i': 7}),
            'float': PB({'f': 2.0}),
            'bool': PB({'b': True}),
            'lisint': PB({'list': PB({'i': 5, 'i': 6})})}
        })

        res = collect_tf_attrs(pb.attr)

        # Reference results for given parameters
        ref = {
            'str': pb.attr['str'].s,
            'int': pb.attr['int'].i,
            'float': pb.attr['float'].f,
            'bool': pb.attr['bool'].b,
            'lisint': pb.attr['lisint'].list.i
        }
        for attr in ref:
            self.assertEqual(res[attr], ref[attr])

class TensorContentParsing(unittest.TestCase):
    def test_list_not_type_no_shape(self):
        pb_tensor = PB(dict(
            dtype=3,
            tensor_content=b'\x01\x00\x00\x00\x02\x00\x00\x00\x03\x00\x00\x00\x04\x00\x00\x00\x05\x00\x00\x00'
        ))
        tf_dtype = pb_tensor.dtype
        shape = np.array([5])
        ref = [1, 2, 3, 4, 5]
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(np.all(res == ref))

    def test_list_type_no_shape(self):
        pb_tensor = PB({
            'dtype': 3,
            'int_val': np.array([1, 2, 3, 4, 5], dtype=np.int32)
        })
        tf_dtype = pb_tensor.dtype
        shape = np.array([5])
        ref = [1, 2, 3, 4, 5]
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(np.all(res == ref))

    def test_list_not_type_shape(self):
        pb_tensor = PB({
            'dtype': 3,
            'tensor_content': b'\x01\x00\x00\x00\x02\x00\x00\x00\x03\x00\x00\x00\x04\x00\x00\x00\x05\x00\x00\x00'
        })
        tf_dtype = pb_tensor.dtype
        shape = np.array([10])
        ref = [1, 2, 3, 4, 5, 5, 5, 5, 5, 5]
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(np.all(res == ref))

    def test_list_type_shape(self):
        pb_tensor = PB({
            'dtype': 3,
            'int_val': np.array([1, 2, 3, 4, 5], dtype=np.int32)
        })
        tf_dtype = pb_tensor.dtype
        shape = np.array([10])
        ref = [1, 2, 3, 4, 5, 5, 5, 5, 5, 5]
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(np.all(res == ref))

    def test_0d_not_type_no_shape(self):
        pb_tensor = PB({
            'dtype': 3,
            'tensor_content': b'\x01\x00\x00\x00',
        })
        tf_dtype = pb_tensor.dtype
        shape = np.array([])
        ref = 1
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(res == ref)

    def test_0d_type_no_shape(self):
        pb_tensor = PB({
            'dtype': 3,
            'int_val': [5],
        })
        tf_dtype = pb_tensor.dtype
        shape = np.array([])
        ref = 5
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(res == ref)

    def test_0d_not_type_shape(self):
        pb_tensor = PB({
            'dtype': 3,
            'tensor_content': b'\x01\x00\x00\x00',
        })
        tf_dtype = pb_tensor.dtype
        shape = np.array([3])
        ref = [1, 1, 1]
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(np.all(res == ref))

    def test_0d_type_shape(self):
        pb_tensor = PB({
            'dtype': 3,
            'int_val': [5],
        })
        tf_dtype = pb_tensor.dtype
        shape = np.array([3])
        ref = [5, 5, 5]
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(np.all(res == ref))

    def test_0d_type_shape_1(self):
        pb_tensor = PB({
            'dtype': 3,
            'int_val': [5],
        })
        tf_dtype = pb_tensor.dtype
        shape = np.array([1])
        ref = [5]
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(np.all(res == ref))

    def test_nd_not_type_no_shape(self):
        pb_tensor = PB({
            'dtype': 3,
            'tensor_content':
                b'\x01\x00\x00\x00\x02\x00\x00\x00\x03\x00\x00\x00\x04\x00\x00\x00\x05\x00\x00\x00\x06\x00\x00\x00',
        })
        tf_dtype = pb_tensor.dtype
        shape = np.array([2, 3])
        ref = [[1, 2, 3], [4, 5, 6]]
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(np.all(res == ref))

    def test_nd_type_no_shape(self):
        pb_tensor = PB({
            'dtype': 3,
            'int_val': [[1, 2, 3], [4, 5, 6]],
        })
        tf_dtype = pb_tensor.dtype
        shape = np.array([2, 3])
        ref = [[1, 2, 3], [4, 5, 6]]
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(np.all(res == ref))

    def test_nd_not_type_shape(self):
        pb_tensor = PB({
            'dtype': 3,
            'tensor_content':
                b'\x01\x00\x00\x00\x02\x00\x00\x00\x03\x00\x00\x00\x04\x00\x00\x00\x05\x00\x00\x00\x06\x00\x00\x00',
        })
        tf_dtype = pb_tensor.dtype
        shape = np.array([2, 5, 2])
        ref = [[[1, 2], [3, 4], [5, 6], [6, 6], [6, 6]],
               [[6, 6], [6, 6], [6, 6], [6, 6], [6, 6]]]
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(np.all(res == ref))

    def test_nd_type_shape(self):
        pb_tensor = PB({
            'dtype': 3,
            'int_val': [[1, 2, 3], [4, 5, 6]],
        })
        tf_dtype = pb_tensor.dtype
        shape = np.array([2, 5, 2])
        ref = [[[1, 2], [3, 4], [5, 6], [6, 6], [6, 6]],
               [[6, 6], [6, 6], [6, 6], [6, 6], [6, 6]]]
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)
        self.assertTrue(np.all(res == ref))

    def test_str_decode(self):
        pb_tensor = PB({
            'dtype': 7,
            'string_val': [b"\037\000\036\000\002\000\303\237\035\000\002"]
        })
        tf_dtype = pb_tensor.dtype
        shape = int64_array([1])
        warning_message = 'ERROR:root:Failed to parse a tensor with Unicode characters. Note that Inference Engine ' \
                          'does not support string literals, so the string constant should be eliminated from the ' \
                          'graph.'
        ref_val = np.array([b'\x1f\x00\x1e\x00\x02\x00\xc3\x9f\x1d\x00\x02'])
        with self.assertLogs(log.getLogger(), level="ERROR") as cm:
            result = tf_tensor_content(tf_dtype, shape, pb_tensor)
            self.assertEqual([warning_message], cm.output)
            self.assertEqual(ref_val, result)

    def test_str_decode_list(self):
        pb_tensor = PB({
            'dtype': 7,
            'string_val': [b'\377\330\377\377\330\377'],
        })
        shape = int64_array([])
        warning_message = 'ERROR:root:Failed to parse a tensor with Unicode characters. Note that Inference Engine ' \
                          'does not support string literals, so the string constant should be eliminated from the ' \
                          'graph.'
        with self.assertLogs(log.getLogger(), level="ERROR") as cm:
            result = tf_tensor_content(pb_tensor.dtype, shape, pb_tensor)
            self.assertEqual([warning_message, warning_message], cm.output)

    def test_empty_value(self):
        pb_tensor = PB({
            'dtype': 1,
            'float_val': []
        })

        shape = int64_array([1, 1, 0])
        tf_dtype = pb_tensor.dtype
        ref = np.array([[[]]], dtype=np.float32)
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)

        self.assertEqual(res.shape, ref.shape)
        self.assertTrue(np.all(res == ref))

    def test_scalar_value(self):
        pb_tensor = PB({
            'dtype': 3,
            'int_val': 4
        })

        shape = int64_array([])
        tf_dtype = pb_tensor.dtype
        ref = np.array(4, dtype=np.int32)
        res = tf_tensor_content(tf_dtype, shape, pb_tensor)

        self.assertEqual(ref, res)

max_stars_repo_name: monroid/openvino

max_issues_repo_name: thomas-yanxin/openvino

max_forks_repo_name: thomas-yanxin/openvino


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

"""Run tabular behavioral cloning with representation learning."""
import os

from absl import app
from absl import flags
from dice_rl.data.dataset import Dataset
from dice_rl.data.tf_agents_onpolicy_dataset import TFAgentsOnpolicyDataset
import dice_rl.environments.gridworld.low_rank as low_rank
import dice_rl.environments.gridworld.navigation as navigation
import dice_rl.environments.gridworld.taxi as taxi
import dice_rl.environments.gridworld.tree as tree
from dice_rl.environments.infinite_frozenlake import InfiniteFrozenLake
from dice_rl.estimators import estimator as estimator_lib
import dice_rl.utils.common as common_utils
import numpy as np
import tensorflow.compat.v2 as tf
from tf_agents.environments import gym_wrapper
from tf_agents.environments import tf_py_environment

from rl_repr.contrastive_fourier.tabular_bc_energy import TabularBCEnergy
from rl_repr.contrastive_fourier.tabular_bc_sgd import TabularBCSGD
from rl_repr.contrastive_fourier.tabular_bc_svd import TabularBCSVD

FLAGS = flags.FLAGS

flags.DEFINE_string('env_name', 'lowrank_tree', 'Environment name.')
flags.DEFINE_integer('seed', 0, 'Initial random seed.')
flags.DEFINE_integer('num_trajectory', 500,
                     'Number of trajectories to collect.')
flags.DEFINE_integer('num_expert_trajectory', 5,
                     'Number of trajectories to collect.')
flags.DEFINE_integer('max_trajectory_length', 3,
                     'Cutoff trajectory at this step.')
flags.DEFINE_float('alpha', -0.125, 'How close to target policy.')
flags.DEFINE_float('alpha_expert', 1.125, 'How close to target policy.')
flags.DEFINE_bool('tabular_obs', True, 'Whether to use tabular observations.')
flags.DEFINE_string('load_dir', '../tests/testdata/',
                    'Directory to load dataset from.')
flags.DEFINE_string('save_dir', './output/', 'Directory to save result to.')
flags.DEFINE_float('gamma', 0.99, 'Discount factor.')

flags.DEFINE_string('embed_learner', None, 'One of [sgd, svd, energy] or None.')
flags.DEFINE_integer('embed_dim', 64, 'Embedder dimension.')
flags.DEFINE_integer('fourier_dim', None, 'Dimension of fourier features.')
flags.DEFINE_float('embed_learning_rate', 0.01, 'Learning rate.')
flags.DEFINE_float('learning_rate', 0.01, 'Learning rate.')
flags.DEFINE_integer('embed_pretraining_steps', 10_000,
                     'Number of pretraining training steps.')
flags.DEFINE_integer('num_steps', 100_000, 'Number of training steps.')
flags.DEFINE_integer('eval_interval', 1000,
                     'Number of pretraining training steps.')
flags.DEFINE_bool('finetune', False,
                  'Whether to finetune embedder during policy learning.')
flags.DEFINE_bool('latent_policy', True, 'Whether to learn a latent policy.')

flags.DEFINE_integer('batch_size', 1024, 'Batch size.')

def get_onpolicy_dataset(env_name, tabular_obs, policy_fn, policy_info_spec):
  """Gets target policy."""
  if env_name == 'taxi':
    env = taxi.Taxi(tabular_obs=tabular_obs)
  elif env_name == 'grid':
    env = navigation.GridWalk(tabular_obs=tabular_obs)
  elif env_name == 'lowrank_tree':
    env = tree.Tree(branching=2, depth=3, duplicate=10)
  elif env_name == 'frozenlake':
    env = InfiniteFrozenLake()
  elif env_name == 'low_rank':
    env = low_rank.LowRank()
  else:
    raise ValueError('Unknown environment: %s.' % env_name)

  tf_env = tf_py_environment.TFPyEnvironment(gym_wrapper.GymWrapper(env))
  tf_policy = common_utils.TFAgentsWrappedPolicy(
      tf_env.time_step_spec(),
      tf_env.action_spec(),
      policy_fn,
      policy_info_spec,
      emit_log_probability=True)

  return TFAgentsOnpolicyDataset(tf_env, tf_policy)

def main(argv):
  env_name = FLAGS.env_name
  seed = FLAGS.seed
  tabular_obs = FLAGS.tabular_obs
  num_trajectory = FLAGS.num_trajectory
  num_expert_trajectory = FLAGS.num_expert_trajectory
  max_trajectory_length = FLAGS.max_trajectory_length
  alpha = FLAGS.alpha
  alpha_expert = FLAGS.alpha_expert
  load_dir = FLAGS.load_dir
  save_dir = FLAGS.save_dir
  gamma = FLAGS.gamma
  assert 0 <= gamma < 1.
  embed_dim = FLAGS.embed_dim
  fourier_dim = FLAGS.fourier_dim
  embed_learning_rate = FLAGS.embed_learning_rate
  learning_rate = FLAGS.learning_rate
  finetune = FLAGS.finetune
  latent_policy = FLAGS.latent_policy
  embed_learner = FLAGS.embed_learner
  num_steps = FLAGS.num_steps
  embed_pretraining_steps = FLAGS.embed_pretraining_steps
  batch_size = FLAGS.batch_size

  hparam_str = ('{ENV_NAME}_tabular{TAB}_alpha{ALPHA}_seed{SEED}_'
                'numtraj{NUM_TRAJ}_maxtraj{MAX_TRAJ}').format(
                    ENV_NAME=env_name,
                    TAB=tabular_obs,
                    ALPHA=alpha,
                    SEED=seed,
                    NUM_TRAJ=num_trajectory,
                    MAX_TRAJ=max_trajectory_length)
  directory = os.path.join(load_dir, hparam_str)
  print('Loading dataset.')
  dataset = Dataset.load(directory)
  print('num loaded steps', dataset.num_steps)
  print('num loaded total steps', dataset.num_total_steps)
  print('num loaded episodes', dataset.num_episodes)
  print('num loaded total episodes', dataset.num_total_episodes)
  estimate = estimator_lib.get_fullbatch_average(dataset, gamma=gamma)
  print('data per step avg', estimate)

  hparam_str = ('{ENV_NAME}_tabular{TAB}_alpha{ALPHA}_seed{SEED}_'
                'numtraj{NUM_TRAJ}_maxtraj{MAX_TRAJ}').format(
                    ENV_NAME=env_name,
                    TAB=tabular_obs,
                    ALPHA=alpha_expert,
                    SEED=seed,
                    NUM_TRAJ=num_expert_trajectory,
                    MAX_TRAJ=max_trajectory_length)
  directory = os.path.join(load_dir, hparam_str)
  print('Loading expert dataset.')
  expert_dataset = Dataset.load(directory)
  print('num loaded expert steps', expert_dataset.num_steps)
  print('num loaded total expert steps', expert_dataset.num_total_steps)
  print('num loaded expert episodes', expert_dataset.num_episodes)
  print('num loaded total expert episodes', expert_dataset.num_total_episodes)
  expert_estimate = estimator_lib.get_fullbatch_average(
      expert_dataset, gamma=gamma)
  print('expert data per step avg', expert_estimate)

  hparam_dict = {
      'env_name': env_name,
      'alpha_expert': alpha_expert,
      'seed': seed,
      'num_trajectory': num_trajectory,
      'num_expert_trajectory': num_expert_trajectory,
      'max_trajectory_length': max_trajectory_length,
      'embed_learner': embed_learner,
      'embed_dim': embed_dim,
      'fourier_dim': fourier_dim,
      'embed_learning_rate': embed_learning_rate,
      'learning_rate': learning_rate,
      'latent_policy': latent_policy,
      'finetune': finetune,
  }
  hparam_str = ','.join(
      ['%s=%s' % (k, str(hparam_dict[k])) for k in sorted(hparam_dict.keys())])
  summary_writer = tf.summary.create_file_writer(
      os.path.join(save_dir, hparam_str, 'train'))

  if embed_learner == 'sgd' or not embed_learner:
    algo = TabularBCSGD(
        dataset.spec,
        gamma=gamma,
        embed_dim=embed_dim,
        embed_learning_rate=embed_learning_rate,
        learning_rate=learning_rate,
        finetune=finetune,
        latent_policy=latent_policy)
  elif embed_learner == 'svd':
    algo = TabularBCSVD(
        dataset.spec,
        gamma=gamma,
        embed_dim=embed_dim,
        learning_rate=learning_rate)
  elif embed_learner == 'energy':
    algo = TabularBCEnergy(
        dataset.spec,
        gamma=gamma,
        embed_dim=embed_dim,
        fourier_dim=fourier_dim,
        embed_learning_rate=embed_learning_rate,
        learning_rate=learning_rate)
  else:
    raise ValueError('embed learner %s not supported' % embed_learner)

  if embed_learner == 'svd':
    embed_dict = algo.solve(dataset)
    with summary_writer.as_default():
      for k, v in embed_dict.items():
        tf.summary.scalar(f'embed/{k}', v, step=0)
        print('embed', k, v)
  else:
    algo.prepare_datasets(dataset, expert_dataset)
    if embed_learner is not None:
      for step in range(embed_pretraining_steps):
        batch = dataset.get_step(batch_size, num_steps=2)
        embed_dict = algo.train_embed(batch)
        if step % FLAGS.eval_interval == 0:
          with summary_writer.as_default():
            for k, v in embed_dict.items():
              tf.summary.scalar(f'embed/{k}', v, step=step)
              print('embed', step, k, v)

  for step in range(num_steps):
    batch = expert_dataset.get_step(batch_size, num_steps=2)
    info_dict = algo.train_step(batch)
    if step % FLAGS.eval_interval == 0:
      with summary_writer.as_default():
        for k, v in info_dict.items():
          tf.summary.scalar(f'bc/{k}', v, step=step)
          print('bc', k, v)

      policy_fn, policy_info_spec = algo.get_policy()
      onpolicy_data = get_onpolicy_dataset(env_name, tabular_obs, policy_fn,
                                           policy_info_spec)
      onpolicy_episodes, _ = onpolicy_data.get_episode(
          100, truncate_episode_at=max_trajectory_length)
      with summary_writer.as_default():
        tf.print('eval/reward', np.mean(onpolicy_episodes.reward))
        tf.summary.scalar(
            'eval/reward', np.mean(onpolicy_episodes.reward), step=step)

if __name__ == '__main__':
  app.run(main)

max_stars_repo_name: DionysisChristopoulos/google-research

max_issues_repo_name: DionysisChristopoulos/google-research

max_forks_repo_name: admariner/google-research


text:
import os
import re
import sys
import random
import numpy as np
import tensorflow as tf

def setup_dirs(args):
  out_dir = args.OUTPUT_DIR
  os.makedirs(out_dir, exist_ok=True)

def seed_everything(seed=1234):
    random.seed(seed)
    os.environ['PYTHONHASHSEED'] = str(seed)
    np.random.seed(seed)
    tf.random.set_seed(seed)

# http://stackoverflow.com/questions/34950201/pycharm-print-end-r-statement-not-working
class Logger(object):
    def __init__(self):
        self.terminal = sys.stdout  #stdout
        self.file = None

    def open(self, file, mode=None):
        if mode is None: mode ='w'
        self.file = open(file, mode)

    def write(self, message, is_terminal=1, is_file=1 ):
        if '\r' in message: is_file=0

        if is_terminal == 1:
            self.terminal.write(message)
            self.terminal.flush()
            #time.sleep(1)

        if is_file == 1:
            self.file.write(message)
            self.file.flush()

    def flush(self):
        # this flush method is needed for python 3 compatibility.
        # this handles the flush command by doing nothing.
        # you might want to specify some extra behavior here.
        pass