text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
    This module implements DeepQL - version 1

    the DeepQL uses an MLP Q-network
    it is retrained only after 'episodes' iterations
    the training uses replay memory

    bibliographic references:
    ========================

    VAN HASSELT, Hado; GUEZ, Arthur; SILVER, David. Deep reinforcement learning with double q-learning.
    In: Thirtieth AAAI conference on artificial intelligence. 2016.

    HASSELT, Hado V. Double Q-learning. In: Advances in Neural Information Processing Systems. 2010. p. 2613-2621.

    WANG, Ziyu et al. Dueling network architectures for deep reinforcement learning. arXiv preprint arXiv:1511.06581, 2015.

    CLEMENTE, Alfredo V.; CASTEJÓN, Humberto N.; CHANDRA, Arjun. Efficient parallel methods for deep reinforcement learning.
    arXiv preprint arXiv:1705.04862, 2017.

    MNIH, Volodymyr et al. Asynchronous methods for deep reinforcement learning.
    In: International conference on machine learning. 2016. p. 1928-1937.

"""
from __future__ import print_function
import time
from datetime import datetime
import numpy as np
from collections import deque
from copy import deepcopy
import pickle

import logging
logging.basicConfig(level=logging.INFO)

def softmax(z):
    """ returns the softmax function (probabilities) given an array z

        @param z: an 1D array of float
        @type z: np.array

        @return: softmax(x)
        @rtype: np.array
    """
    assert len(z.shape) == 1
    e_x = np.exp(z - np.max(z))
    return e_x / sum(e_x)

def softmax_2d(z):
    """
        @param z: an array (2, n)
        @type z: np.array
    """
    assert len(z.shape) == 2
    s = np.max(z, axis=1)
    s = s[:, np.newaxis]  # necessary step to do broadcasting
    e_x = np.exp(z - s)
    div = np.sum(e_x, axis=1)
    div = div[:, np.newaxis]  # to broadcast corretly
    return e_x / div

class DeepQL(object):

    """
        ref. https://keon.io/deep-q-learning/
             https://github.com/simoninithomas/deep_q_learning/blob/master/DeepQL%20Cartpole.ipynb
             https://medium.com/@gtnjuvin/my-journey-into-deep-q-learning-with-keras-and-gym-3e779cc12762
             https://medium.com/@awjuliani/simple-reinforcement-learning-with-tensorflow-part-4-deep-q-networks-and-beyond-8438a3e2b8df
    """

    TO_MUCH_ERROR_IN_A_ROW = 20  # if that much error occurs, stop the program. None disable this test.

    def __init__(self,
                 env,  # the environment object
                 model,  # Keras model
                 memory,
                 timesteps=1,
                 epsilon=0.01,  # exploration rate
                 epsilon_min=0.1,
                 epsilon_decay=0.995,
                 learning_rate=0.001,
                 gamma=0.95,
                 batch_size=32,  # size of the batch in the replay
                 episodes=30,  # number of iterations before replay occurs
                 epochs=1,
                 interaction_interval=30,  # wait 30 second, before another cycle
                 log_level=logging.DEBUG,
                 **kwargs):
        """
            @param env: the environment class
            @param model: the Keras model used to approximate the Q-function
            @param memory: the replay memory implementation
        """
        self.interaction_interval = interaction_interval

        self.log = logging.getLogger('DeepQL')
        self.log.setLevel(log_level)
        # self.log.addHandler(logging.StreamHandler(sys.stdout))
        self.log.info("Debug {} activated".format("is" if log_level == logging.DEBUG else "is **NOT**"))

        assert epsilon_decay is not None and epsilon_decay >= 0 and epsilon_decay <= 1, 'epsilon should be in [0,1]'
        self.epsilon = epsilon

        assert epsilon_min is None or (epsilon_min >= 0 and epsilon_min <= epsilon_decay), 'min epsilon should be in [0, epsilon]'
        self.epsilon_min = epsilon_min

        assert epsilon_decay is None or (epsilon_decay >= 0 and epsilon_decay <= 1), 'epsilon_decay should be in [0,1]'
        self.epsilon_decay = epsilon_decay

        self.log.info("Epsilon {} Eps min {} Decay {}".format(epsilon, epsilon_min, epsilon_decay))

        self.env = env

        self.action_size = env.action_size
        self.log.info("Action size {}".format(self.action_size))

        self.state_size = env.state_size
        self.state_dim = env.state_dim
        self.log.info("State dim {}".format(self.state_dim))

        self.memory = memory
        self.replay_counter = 0
        self.batch_size = batch_size if batch_size > 0 else 1
        self.log.info("Replay memory: batch size: {}".format(self.batch_size))

        self.episodes = episodes if episodes > 0 else 1
        self.log.info("Episodes: {}".format(self.episodes))

        self.learning_rate = learning_rate
        self.gamma = gamma    # discount rate
        self.model = model

        self.epochs = epochs
        self.timesteps = timesteps

        self.can_run = True
        self.runs = 0
        # create a container to save 'timesteps' states for each AP
        self.prev_states = [deque(maxlen=timesteps) for i in range(len(env.aps))]

        self.log.info("Agent created")

    def save_model(self, model_filename='model.json'):
        """ save the model to a json file and the weights to a h5 file

            @param model_filename: the filename with '.json' extension
        """
        model_json = self.model.to_json()
        with open(model_filename, "w") as json_file:
            json_file.write(model_json)
        self.model.save_weights(model_filename.replace('.json', '.h5'))
        opt = {'name': str(self.model.optimizer),
               'config': self.model.optimizer.get_config()
               }
        pickle.dump(opt, open(model_filename.replace('.json', '.p'), 'wb'), protocol=pickle.HIGHEST_PROTOCOL)

    def remember(self, states, actions, next_states, rewards):
        """ pushes s, a, s', r
            @param states: list of initial state, one for each AP
            @param actions: list of actions taken, one for each AP
            @param next_states: list of next state, one for each AP
            @param rewards: list of rewards, one for each AP
        """
        self.log.debug("Push to replay: s:{} a:{} s'{} r{}".format(states, actions, next_states, rewards))
        self.memory.push(states, actions, next_states, rewards)

    def format_state_to_predict(self, values, batch_size=1):
        """ formats to use in predict, because
            predict needs first dimension ==> entries
            followed by the other dimension in values

            @param values: list of values to convert to a numpy array
            @param batch_size: defines the size of the batch (first dimension size)
            @return: a numpu array with self.timesteps,
                     composed by the self.prev_states values (saved from previous runs) and
                     the value passed as parameter
        """
        if type(values) == int:
            values = [values, ]
        v = np.array(values)
        if self.timesteps == 1:
            if len(v.shape) == 1:
                v = v.reshape(batch_size, -1)
        else:
            v = v.reshape(batch_size, self.timesteps, -1)
        self.log.debug('format_state_to_predict batch_size {} - self.timesteps {} > {}'.format(batch_size, self.timesteps, v))
        return v

    def get_q_max(self, sprime):
        """ the Q_max is calculated using the model network
            notice that you don't need to call self.format_state_to_predict() for sprime
            sprime format depends on the number of time steps, thus
            dim(s') = (1, timesteps, num_features)
            ... num_features = self.state_dim

            @param next_state: the next state s'
            @return: the Q_max for the state s'

            Q_max = max Q(s', a')
                     a'
        """
        q_prediction = self.model.predict(sprime)
        q_max = np.amax(q_prediction, axis=1)
        self.log.debug("s': {} Q max: {}".format(sprime, q_max))
        return q_max

    def update_epsilon(self):
        """ perform epsilon decay.
            To prevent the decay to occur, just set epsilon_decay to None.
            If epsilon_min is None, then decays forever.
            Otherwise decays while epsilon > epsilon_min
        """
        if self.epsilon_decay is not None and (self.epsilon_min is None or self.epsilon > self.epsilon_min):
            self.epsilon *= self.epsilon_decay
            self.log.info("New epsilon: {}".format(self.epsilon))

    def replay(self):
        """
            decides if the replay will occur, if not just returns
            uses the memory to recover a mini-batch that will be used to train the model network

        @return: True if replay occured
        """
        self.log.info('Doing replay')
        minibatch = self.memory.sample(batch_size=self.batch_size)
        self.log.debug("Minibatch:\t{}".format(minibatch))
        X = []
        Y = []
        for sequence in minibatch:
            self.log.debug("Replaying the sequence {}".format(sequence))
            # retrain the sequence (notice the steps):
            # 1) get the s' (sequence)
            sprime = self.format_state_to_predict([s.next_state for s in sequence])
            self.log.debug("next_state {}".format(sprime))

            # 2) calculate the Q-value of this element
            reward = sequence[-1].reward  # reward in the last step
            new_value = reward + self.gamma * self.get_q_max(sprime)
            self.log.debug("old reward {} --> predicted reward {}".format(reward, new_value))

            # 3) make the agent to approximately map the current state to future discounted reward
            #    has to pass to the Q-function the sequence of states (seq_state)
            #   also get w
            s = self.format_state_to_predict([x.state for x in sequence])
            self.log.debug("state {}".format(s))

            action = sequence[-1].action
            self.log.debug("action {}".format(action))

            # 4) get what is predicted for the state s (sequence)
            #    and change the returned value for the action chosen to the calculated 'new_value'
            target_f = self.model.predict(s)

            target_f[0][action] = new_value

            self.log.debug('target_f {} - shape {}'.format(target_f, target_f.shape))
            X.append(s)
            Y.append(target_f[0])

        # 6) retrain the network with the new value for the action
        # s = (batch_size, timesteps, num_features)
        # target_f = (batch_size, 1)
        n = len(minibatch)
        if n > 0:
            X = np.array(X).reshape(n, self.timesteps, -1)
            Y = np.array(Y).reshape(n, -1)

            self.log.info("X{} = {}".format(X.shape, X))  # remove print
            self.log.info("Y{} = {}".format(Y.shape, Y))  # remove print

            self.model.fit(X, Y, epochs=self.epochs, verbose=0)

            # perform decay
            self.update_epsilon()

            self.replay_counter += 1
            return True  # replay occured
        else:
            return False

    def predict_action_output(self, curr_state):
        """ Predict the reward value based on the given state
            this method formats 'curr_state' using self.format_state_to_predict() in order to call the keras predict()

            @param curr_state: the current state (one for each device)
            @return: values for all the actions
            @rtype: list
        """

        # compose the current state with the previous states to predict
        seq_states = deepcopy(self.prev_states)  # prevent from changing the original data
        n = len(curr_state)
        for i in range(n):
            s = curr_state[i]
            seq_states[i].append(s)
        batch_size = len(seq_states)
        seq_states = self.format_state_to_predict(seq_states, batch_size=batch_size)

        # self.log.debug("batch_size {}".format(batch_size))
        # self.log.debug("self.prev_states {}".format(self.prev_states))
        # self.log.debug("seq_states {} {}".format(seq_states, seq_states.shape))

        # predict the action using the Q-network
        act_values = self.model.predict(seq_states)  # return an array of action_size elements
        self.log.info("Predicted action: {}".format(act_values))
        return act_values

    def get_action_eps_greedy(self, curr_state):
        """ select the action (one for each ap), epsilon greedy way
            @param curr_state: a list of current states, one for each AP

            @return: list[int]: each entry is a number that represents the action for that state
        """
        if np.random.rand() <= self.epsilon:
            # The agent acts randomly
            self.log.info("Random action")
            a = [np.random.choice(self.env.valid_actions(curr_state), 1)[0] for state in curr_state]
        else:
            # Pick the action based on the predicted reward
            act_values = self.predict_action_output(curr_state)
            self.log.info("act_values {}".format(act_values))  # remove print
            # todo: what happens if we have 2 or more actions with the same value in the same axis?
            #  you should select randomly between them
            a = np.argmax(act_values, axis=1)
        self.log.info('Selected action: {}'.format(a))
        return a

    def get_action_boltzmann(self, curr_state):
        """ select the action (one for each ap), using boltzmann
            @param curr_state: a list of current state, one for each AP

            @return: list[int]: each entry is a number that represents the action for that state
        """
        # Pick the action based on the predicted reward
        act_values = self.predict_action_output(curr_state)  # get the Q-value for each action in state

        # TODO: check the correctness of this code below
        #       probably error dealing with axis in softmax_2d
        p = softmax_2d(act_values)  # convert Q to probabilities)
        # select an action based on the probability
        a = [np.random.choice(self.n_actions, 1, p=p_action)[0] for p_action in p]
        self.log.info('Selected action: {}'.format(a))
        return a

    def get_action(self, states):
        """ overwrite this method to call self.get_action_eps_greedy() or self.get_action_boltzmann() to implement the search policy
            @param states: a list of states, one for each AP
        """
        return self.get_action_eps_greedy(states)

    def stop(self):
        """ change the stopping flag in the run(), so the program will stop at the end of the iteration
        """
        self.can_run = False

    @property
    def number_of_runs(self):
        """ number of times the program iteracted and acted upon the environment
            @rtype: int
        """
        return self.runs

    def stop_running(self):
        """
            change the flag that controls the while loop in run()
            so the agent stop at the end of that execution
        """
        self.can_run = False

    def run(self, run_id=1, wait_for_states=10, save_iterations=20):
        """ executes the control loop

            @param wait_for_states: how much time should sleep between get_states request
            @type wait_for_states: int

            @param save_iterations: every 'save_iterations' iterations, save the model and the weights

            @return: if the agents detected to much errors in a row
            @rtype: bool
        """
        assert wait_for_states > 0, "wait_for_states {} is not > 0".format(wait_for_states)

        # need to save (timesteps - 1) states to feed the network during prediction.
        states = None
        while states is None:
            # wait until a valid state is retrieved
            states = self.env.get_states()
            time.sleep(wait_for_states)
        # fill the previous state memory
        for i in range(len(self.env.aps)):
                s = states[i]
                self.prev_states[i] = deque([s for _ in range(self.timesteps)], maxlen=self.timesteps)
        self.log.info("len(self.prev_states) {}".format(len(self.prev_states)))
        self.log.info("len(self.env.aps) {}".format(len(self.env.aps)))
        self.log.info("")
        assert len(self.prev_states) == len(self.env.aps), "Something wrong with prev_states size"

        self.can_run = True
        self.runs_with_error = 0
        self.log.info("Starting running agent")

        none_states = 0
        while self.can_run and not self.env.done and (self.TO_MUCH_ERROR_IN_A_ROW is None or self.runs_with_error < self.TO_MUCH_ERROR_IN_A_ROW):
            self.log.info("Run #{}".format(self.runs))
            got_an_error = False
            t0 = datetime.now()
            self.log.info("Starting @{}".format(t0))

            # get states for each AP
            states = self.env.get_states()  # one for each AP
            self.log.info("States: {}".format(states))
            if states is None or None in states:
                got_an_error = True
                none_states += 1
            else:
                none_states = 0

                # the states are ok, so go on...
                # Decide action
                actions = self.get_action(states)
                self.log.info("Actions: {}".format(actions))

                # perform the selected actions, retrieve the next state and the reward
                next_states, rewards = self.env.make_step(actions)
                self.log.info("New state: {}".format(next_states))
                self.log.info("Rewards: {}".format(rewards))

                # check if an error occurred
                got_an_error = np.nan in rewards

            if not got_an_error:
                # save moves for posterior training
                self.remember(states, actions, next_states, rewards)

                # make next_state the new current state for the next frame.
                states = next_states

                self.runs += 1  # increment the number of iterations
                # train the model network
                # self.replay() decides if the replay will occur
                self.replay()

                # add actual state to prev_states
                for i in range(len(states)):
                    s = states[i]
                    self.prev_states[i].append(s)

                self.runs_with_error = 0  # reset the error counter
            else:
                # add another error in a row
                self.runs_with_error += 1

            t1 = datetime.now()
            _dur = (t1 - t0).seconds
            self.log.info("Ended @{} - duration {} s".format(t1, _dur))

            # wait interaction_interval seconds
            if self.interaction_interval is not None:
                _sleep = max(self.interaction_interval - _dur, 0)
                if _sleep > 0:
                    self.log.debug("Sleep {} s".format(_sleep))
                    time.sleep(_sleep)

            if self.runs % save_iterations == 0:
                # save the intermediate weights
                self.save_model()

        self.log.info("Total of {} steps".format(self.runs))
        if not self.can_run:
            self.log.info("Execution time expired.")

        if self.TO_MUCH_ERROR_IN_A_ROW is not None and self.runs_with_error > self.TO_MUCH_ERROR_IN_A_ROW:
            self.log.info("Out of tries. {} errors in a row during execution.".format(self.runs_with_error))

        if self.env.done:
            self.log.info("Game ended!!!")

        # return indicator if too much errors occured
        return self.TO_MUCH_ERROR_IN_A_ROW is not None and self.runs_with_error > self.TO_MUCH_ERROR_IN_A_ROW


text:
import argparse
import os
import time
import numpy as np
from KurtGroup.Kurt import output_processing as op
from functools import partial
from multiprocessing import Pool, cpu_count
import matplotlib.pyplot as plt
from matplotlib import rc
from types import FunctionType
import json
import copy

#*************************** INPUT PARSING ****************************

class TerminalInformation():

    def __init__(self, total_nr: int, max_filename_length: int = None):
        self.full_bar = "▇"
        self.empty_bar = " "
        self.total_nr = total_nr
        self.max_filename_length = max_filename_length
        self.terminal_width, _ = os.get_terminal_size()

        self.UP = lambda x : f"\x1B[{x}A"
        self.CLR = "\x1B[0K"
        self.CLRnl = "\x1B[0K\n"
        self.additionalLines = 0

        self.setProgressbarWidth(True)
        self.updateProgressbar(0, False, False)

    def setProgressbarWidth(self, timing: bool):
        self.progressbar_scale_factor = self.terminal_width
        self.progressbar_scale_factor -= 2 # Corresponds to removing the sidebars
        self.progressbar_scale_factor -= len(str(self.total_nr)) * 2 - 2 # Removes the counter
        if timing:
            self.progressbar_scale_factor -= 23 # Removes the timing parts except the average files pr. second time
            self.progressbar_scale_factor -= 5 # Estimate for average files pr. second time
        self.progressbar_scale_factor -= 10 # Additional just in case

    def updateProgressbar(self, nr: int, print_filename: bool, timing: bool, *, filename: str = None):
        self.progress = nr/self.total_nr*self.progressbar_scale_factor
        self.additionalLines = 0

        self.bar = int(self.progress)*self.full_bar + (self.progressbar_scale_factor-int(self.progress))*self.empty_bar

        self.progressbar = f"|{self.bar}| {nr}/{self.total_nr}"

        if timing:
            self.progressbar += self.timer(nr)

        if print_filename:
            if nr == 1:
                print("")
            self.additionalLines += 1
            self.progressbar = f"<<< {filename:{self.max_filename_length}} >>>{self.CLRnl}{self.progressbar}"

        lines = self.progressbar.split('\n')
        if isinstance(lines, list) and len(lines) > 1:
            self.progressbar = f"{self.UP(len(lines)-1)}{self.progressbar}"

        print(self.progressbar + self.CLR, end="\r")

    def start_timer(self):
        start_time = time.perf_counter()
        self.times = [start_time]
        self.time_differences = []

    def timer(self, nr):
        self.times.append(time.perf_counter())

        self.time_differences.append(self.times[-1] - self.times[-2])

        minutes = lambda x : x // 60
        seconds = lambda x : x - x // 60

        if len(self.times) > 1:
            time_average = np.mean(self.time_differences)
            time_spent = sum(self.time_differences)
            time_left = time_average*(self.total_nr - nr)
            return f"; {minutes(time_spent):02.0f}:{seconds(time_spent):02.0f}<{minutes(time_left):02.0f}:{seconds(time_left):02.0f}; {1/time_average:.2f} files/s"
        else:
            return "; 00:00<00:00; 0 files/s"

# Stolen from https://geekflare.com/flatten-list-python/ and modified slightly
def flatten_list(data, flat_list):
    # iterating over the data
    for element in data:
        # checking for list
        if type(element) == list:
            # calling the same function with current element as new argument
            flatten_list(element, flat_list)
        else:
            flat_list.append(element)

def Resize(array: list) -> None:
    """Takes an array of arrays with varying sizes and resizes them to the same size.
    This is done by appending 'NaN' to the smaller lists.
    If the first value says 'Not implemented' the remainder of the array will be filled with 'Not implemented

    Args:
        array (list): Array of arrays with varying sizes
    """
    max_size = 0
    for arr in array:
        if len(arr) > max_size:
            max_size = len(arr)

    for i, arr in enumerate(array):
        if arr == ['Not implemented']:
            array[i] *= max_size
        else:
            array[i] += ['NaN'] * (max_size - len(arr))

def Data_Extraction(infile, Needed_Values: dict, quiet: bool = False, Temperature: float = 298.15) -> dict:
    Extracted_values = dict()

    infile = op.OutputType(str(infile), Quiet=quiet, Temperature=Temperature)

    # Extracting data
    Extract_data(quiet, Needed_Values, infile.filename, infile.extract, infile.input)

    # Removing the file text from memory as this is no longer needed (and fills in memory)
    infile.extract.__delattr__('lines')

    # List of all dictionary keys for infile.extract
    dict_keys = [*infile.extract.__dict__.keys()]

    # Dictionary
    collection_dict = dict()

    # Collecting the data in dictionaries
    for i in dict_keys[1:]:
        collection_dict[i] =  infile.extract.__dict__[i]

    # Assigning to the Extracted_values dictionary with the filename as key so all data can be easily found in the future
    Extracted_values[infile.filename] = collection_dict

    return Extracted_values

def Extract_data(suppressed: bool, Wanted_Values: dict, infile: str, file_text: dict, input_type: str) -> None:
    # Loops over all requested values and runs the corresponding function
    # If the function has not been implemented it will print an error message
    for i in Wanted_Values:
        try:
            method = getattr(type(file_text),i)
            method(file_text)
        except AttributeError:
            if not(suppressed):
                with open("collect_data.log", "a") as logfile:
                    logfile.write(f'{infile}: {i} has not been implemented for {input_type}\n')

def Check_if_Implemented(input_file: dict, Set_of_values: dict, Extracted_values: dict) -> None:
    # Checks to see if the keys of a double dictionary exists
    # If they don't it is assumed that the function related to the data hasn't been implemented
    for infile in input_file:
        for values in Set_of_values.values():
            for val in values:
                try:
                    Extracted_values[infile][val]
                except KeyError:
                    Extracted_values[infile][val] = ['Not implemented']

def Collect_and_sort_data(input_file: str, Set_of_values: dict, Extracted_values: dict) -> dict:
    Final_arrays = dict()

    # All data from Extracted_values is sorted by the Set_of_values
    # This way excess values in Extracged_values are ignored
    for key in Set_of_values:
        for val in Set_of_values[key]:
            Final_arrays[val] = []
            for infile in input_file:
                Final_arrays[val].append(Extracted_values[infile][val])
    return Final_arrays

def Downsizing_variable_arrays(Outputs: dict, Variable_arrays: dict, count: int, Final_arrays: dict) -> None:
    # Downsizes arrays if requested
    # If not requested nothing happens
    for item in Variable_arrays.items():
        if item[1] > 0:
            for val in Outputs[item[0]]:
                for file in range(0,count):
                    Final_arrays[val][file] = Final_arrays[val][file][0:item[1]]

def Upsizing_variable_arrays(Outputs: dict, Variable_arrays: dict, count: int, Final_arrays: dict, Arguments: dict) -> None:
    # Upsizes arrays if requested
    # If not requested nothing happens
    for key, arg in Variable_arrays.items():
        for val in Outputs[key]:
            for file in range(0,count):
                Final_arrays[val][file] += ['NaN'] * (arg-len(Final_arrays[val][file]))

def Create_Header(Header_text: dict, Set_of_values: dict, Final_arrays: dict) -> list:
    header = ['File']
    # Adds to the header row all relevant headers for the data-points requested
    for key in Set_of_values.keys():
        for val in Set_of_values[key]:
            if len(Final_arrays[val][0]) > 1:
                for i in range(len(Final_arrays[val][0])):
                    header.append(f'{Header_text[val]} {i+1}')
            else:
                header.append(Header_text[val])
    return header

def Fill_output_array(Set_of_values: dict, array_input: dict, count: int, Final_arrays: dict, output_array: list) -> None:
    # Fills the output array with the header row
    if count == 1:
        output_array[1,0] = array_input[0][0]
    else:
        output_array[1:,0] = np.array(np.concatenate(array_input))

    # Fills the output array with the requested data
    # As it is not known the length of all lists in the Final_arrays the col is used to slowly move through the array
    col = 1
    for key in Set_of_values.keys():
        for val in Set_of_values[key]:
            output_array[1:,col:col+len(np.array(Final_arrays[val][0]))] = np.array(Final_arrays[val])
            col += len(np.array(Final_arrays[val][0]))

def Make_complex_propagator_spectrum(input_file: list, suppressed: bool, Format: str, Extracted_Values: dict, SAVE: bool = True) -> None:
    # A LOT OF PLOT SETUP
    rc('text', usetex=True)
    xlabel_font = ylabel_font = title_font = 16
    plt.rc('font', size=12) # x/y axis font size
    NA=6.02214199E23 #a vogadros number
    c=299792458 # speed of light
    eps0=8.8541878176E-12
    aufreq=4.134137334E16
    hartreetohz=6.579683920502E15
    Save_Dict = dict()
        # PLOT SETUP DONE
    for file in input_file:
        filename = file.replace('.out','')
        plotname = f'{filename}-complex-propagator.{Format}'
        title = filename.replace("_"," ")

        frequencies = np.array([x[0] for x in Extracted_Values[file]['complex_propagator'] if x != 'NaN' and x != 'Not implemented'])
        imaginary_part = np.array([x[2] for x in Extracted_Values[file]['complex_propagator'] if x != 'NaN' and x != 'Not implemented'])

        if len(frequencies) == 0:
            if not(suppressed):
                print(f'Polarizability dampening frequencies have either not been implemented for this output type, or none were found in the file {filename}')
            continue

        span = c/(frequencies*hartreetohz)*10**9
        graph = (frequencies*aufreq*10000)/(c*eps0) * (imaginary_part*1.64877727436*10**(-41)*NA)/(2.3*1000)

        plt.title(title)
        plt.plot(span, graph)
        plt.ylim((0,max(graph)*1.2))
        plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
        plt.xlabel('Wavelength $(nm)$', fontsize=xlabel_font)
        plt.ylabel('Extinction coefficient', fontsize=ylabel_font)
        plt.xticks(fontsize=16)
        plt.yticks(fontsize=16)
        plt.savefig(plotname, format=f'{Format}', dpi=600)
        plt.close()

        if SAVE:
            Save_Dict[file] = [span,graph]
    if SAVE:
        np.savez('complex_propagator.npz', **Save_Dict)
        print(f'Complex propagator spectrum data has been saved in complex_propagator.npz')

def UVVIS_Spectrum(t: list, l: list, f: list, k: float, sigmacm: float):
    lambda_tot = np.zeros(len(t))
    for x in range(1,len(t)):
        lambda_tot[x] = sum((k/sigmacm)*f*np.exp(-4*np.log(2)*((1/t[x]-1/l)/(1E-7*sigmacm))**2))
    return lambda_tot

def Make_uvvis_spectrum(input_file: list, suppressed: bool, UVVIS_Spectrum: FunctionType, Format: str, Extracted_Values: dict, SAVE: bool = True) -> None:
    # A LOT OF PLOT SETUP
    rc('text', usetex=True)
    xlabel_font = ylabel_font = title_font = 16
    plt.rc('font', size=12) # x/y axis font size
    N=1000 # number of calculated points in curve
    NA=6.02214199*10**23 #avogadros number
    c=299792458 #speed of light
    e=1.60217662*10**(-19) #electron charge
    me=9.10938*10**(-31) #electron mass
    pi=np.pi
    epsvac=8.8541878176*10**(-12)
    sigmacm=0.4*8065.544
    k=(NA*e**2)/(np.log(10)*2*me*c**2*epsvac)*np.sqrt(np.log(2)/pi)*10**(-1)
    inv_cm_to_au = 1/219474.63068
    Save_Dict = dict()
        # PLOT SETUP DONE
    for file in input_file:
        filename = file.replace('.out','')
        plotname = f'{filename}-uvvis.{Format}'
        title = filename.replace("_"," ")

        excitations = np.array([x for x in Extracted_Values[file]['exc_energies'] if x != 'NaN' and x != 'Not implemented'])
        oscillations = np.array([x for x in Extracted_Values[file]['osc_strengths'] if x != 'NaN' and x != 'Not implemented'])

        if len(excitations) == 0 and len(oscillations) == 0:
            if not(suppressed):
                print(f'Excitation energies and oscillator strengths have either not been implemented for this output type, or none were found in the file {filename}')
            continue
        elif len(excitations) == 0:
            if not(suppressed):
                print(f'Excitation energies have either not been implemented for this output type, or none were found in the file {filename}')
            continue
        elif len(oscillations) == 0:
            if not(suppressed):
                print(f'Oscillator strengths have either not been implemented for this output type, or none were found in the file {filename}')
            continue
        elif len(excitations) > len(oscillations):
            excitations = excitations[0:len(oscillations)]
        elif len(oscillations) > len(excitations):
            oscillations = oscillations[0:len(excitations)]

        excitations = 1E7/(excitations/ inv_cm_to_au)   # From a.u. to cm^-1 to nm
        span = np.linspace(min(excitations)-20, max(excitations)+20, N, endpoint=True) # exctinction coefficient (wavelength range)

        graph = UVVIS_Spectrum(span, excitations, oscillations, k, sigmacm)
        plt.title(title)
        plt.plot(span, graph)
        plt.ylim((0,max(graph)*1.2))
        plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
        plt.xlabel('Wavelength $(nm)$', fontsize=xlabel_font)
        plt.ylabel('Extinction coefficient', fontsize=ylabel_font)
        plt.xticks(fontsize=16)
        plt.yticks(fontsize=16)
        plt.savefig(plotname, format=f'{Format}', dpi=600)
        plt.close()

        if SAVE:
            Save_Dict[file] = [span,graph]
    if SAVE:
        np.savez('UVVIS.npz', **Save_Dict)
        print(f'UVVIS spectrum data has been saved in UVVIS.npz')

    del Save_Dict # Deletes dictionary from memory since it is no longer needed

def Spectra(args):
    """
    This function is used to do any methods related to the Spectra keyword
    """
    InputFiles = args.infile
    UVVis = args.uvvis
    ComplexPropagator = args.complex_propagator
    Format = args.format
    Save = args.save
    Quiet = args.quiet
    Multiprocessing = args.multiprocessing

    if UVVis:
        NeededArguments = {'_Excitation_energies': -1, '_Oscillator_strengths': -1}
        ArgumentsToValues = {'_Excitation_energies': ['exc_energies'], '_Oscillator_strengths': ['osc_strengths']}

        if Multiprocessing:
            with Pool(int(cpu_count()/2)) as pool:
                ExtractedValues = pool.map(partial(Data_Extraction, Needed_Values=NeededArguments, quiet=Quiet), InputFiles)
                ExtractedValues = {key: value for dictionary in ExtractedValues for key, value in dictionary.items()} # Reformatting Extracted_values
        else:
            ExtractedValues = dict()
            for infile in InputFiles:
                ExtractedValues[infile] = Data_Extraction(infile, NeededArguments, Quiet)[infile]

        Check_if_Implemented(InputFiles, ArgumentsToValues, ExtractedValues)   # Finding functions not implemented

        for OuterKey in ExtractedValues:   # Turn everything into lists
            for InnerKey, Value in ExtractedValues[OuterKey].items():
                if not isinstance(Value, list):
                    ExtractedValues[OuterKey][InnerKey] = [Value]

        Make_uvvis_spectrum(InputFiles, Quiet, UVVIS_Spectrum, Format, ExtractedValues, Save)
        return

    elif ComplexPropagator:
        NeededArguments = {'_Complex_propagator': True}
        ArgumentsToValues = {'_Complex_propagator': ['complex_propagator']}

        if Multiprocessing:
            with Pool(int(cpu_count()/2)) as pool:
                ExtractedValues = pool.map(partial(Data_Extraction, Needed_Values=NeededArguments, quiet=Quiet), InputFiles)
                ExtractedValues = {key: value for dictionary in ExtractedValues for key, value in dictionary.items()} # Reformatting Extracted_values
        else:
            ExtractedValues = dict()
            for infile in InputFiles:
                ExtractedValues[infile] = Data_Extraction(infile, NeededArguments, Quiet)[infile]

        Check_if_Implemented(InputFiles, ArgumentsToValues, ExtractedValues)   #Finding functions not implemented

        for OuterKey, Dict in ExtractedValues.items():   # Turn everything into lists
            for InnerKey, Value in Dict.items():
                if not isinstance(Value, list):
                    ExtractedValues[OuterKey][InnerKey] = [Value]

        Make_complex_propagator_spectrum(InputFiles, Quiet, Format, ExtractedValues, Save)
        return

    print('Doing nothing - please use --uvvis or --complex-propagator')

def Extract(args):
    """
    This function is used for any methods related to the extract keyword
    """

    # Setting the input_files variable to all files
    InputFiles = args.infile

    # These are all the possible arguments that extract data
    # They are here set to correspond the argument passed
    RequestedArguments = {
        '_Energy': args.energy,
        '_ZPV': args.zpv,
        '_Dipole_moments': args.dipole,
        '_Polarizabilities': args.polar,
        '_Excitation_energies': args.exc,
        '_Oscillator_strengths': args.osc,
        '_Frequencies': args.freq,
        '_Enthalpy': args.enthalpy,
        '_Entropy': args.entropy,
        '_Gibbs': args.gibbs,
        '_PartitionFunctions': args.partfunc,
        '_CPUS': args.cpu_time,
        '_Optimized_Geometry': args.optgeom
    }

    # These are the datapoints that will be extracted per argument
    # The value for each key is how that value is saved in the Classes
    Outputs = {
        '_Energy': ['tot_energy'],
        '_ZPV': ['zpv'],
        '_Enthalpy': ['enthalpy'],
        '_Entropy': ['entropy'],
        '_Gibbs': ['gibbs'],
        '_Dipole_moments': ['dipolex', 'dipoley', 'dipolez', 'total_dipole'],
        '_Polarizabilities': ['polx', 'poly', 'polz', 'iso_polar'],
        '_Excitation_energies': ['exc_energies'],
        '_Oscillator_strengths': ['osc_strengths'],
        '_Frequencies': ['freq'],
        '_PartitionFunctions': ['qTotal'],
        '_CPUS': ['total_cpu_time', 'wall_cpu_time'],
    }

    # These are what will be written in the header for each data-point
    # The keys are what the data-points are saved as in the classes
    # The values are what the data-points should be written under in the output header
    HeaderText = {
        'tot_energy': 'Energy',
        'zpv': 'ZPV EnergTotaly',
        'enthalpy': 'Enthalpy',
        'entropy': 'Entropy (kJ/(mol*K)',
        'gibbs': 'Gibbs Free Energy',
        'dipolex': 'Dipole x',
        'dipoley': 'Dipole y',
        'dipolez': 'Dipole z',
        'total_dipole': 'Total Dipole',
        'polx': 'Polarizability xx',
        'poly': 'Polarizability yy',
        'polz': 'Polarizability zz',
        'iso_polar': 'Isotropic Polarizability',
        'exc_energies': 'Exc. energy',
        'osc_strengths': 'Osc. strength',
        'freq': 'Frequency',
        'qTotal': 'Total molar partition function',
        'total_cpu_time': f'Total CPU time ({args.cpu_time})',
        'wall_cpu_time': f'Wall CPU time ({args.cpu_time})'
    }

    # Remainder of arguments
    T = args.temp
    Save = args.save
    Quiet = args.quiet
    Multiprocessing = args.multiprocessing
    ProgressBar = args.progressbar
    UnitTesting = args.unittest
    SaveName = args.savename

    # Making a copy of RequestedArguments
    # This is so arguments that are dependent on others can be called independently
    NeededArguments = RequestedArguments.copy()

    # Ensuring that excitation energies are calculated when oscillator strengths are
    if NeededArguments['_Oscillator_strengths'] and NeededArguments['_Excitation_energies'] == None:
        NeededArguments['_Oscillator_strengths'] = NeededArguments['_Excitation_energies'] = -1
    elif NeededArguments['_Oscillator_strengths']:
        NeededArguments['_Oscillator_strengths'] = NeededArguments['_Excitation_energies']

    # Ensuring that enthalpies and entropies are calculated as these are needed to calculate the Gibbs free energy
    if NeededArguments['_Gibbs']:
        NeededArguments['_Enthalpy'] = True
        NeededArguments['_Entropy'] = True

    # Ensuring that frequencies are calculated as these are needed to calculate the Partition Functions
    if NeededArguments['_PartitionFunctions']:
        NeededArguments['_Frequencies'] = -1

    # Ensuring that frequencies are calculated as these are needed to calculate the enthalpy
    if NeededArguments['_Enthalpy']:
        NeededArguments['_Frequencies'] = -1
        NeededArguments['_Energy'] = True

    # Ensuring that frequencies are calculated as these are needed to calculate the entropy
    if NeededArguments['_Entropy']:
        NeededArguments['_Frequencies'] = -1

    # Dictionary of data where the amount of values printed can be changed
    # Examples of this are the Excitation energies and the Frequencies
    VariableArrays = dict([item for item in RequestedArguments.items() if type(item[1]) == int])

    # List of arguments that have been requested
    WantedValues = [key for key, val in RequestedArguments.items() if not(val == None or val == False)]

    # List of arguments that are needed for those requested
    # This is done so you can request as an example the Gibbs free energies withou also having the enthalpies and entropies printed
    NeededValues = [key for key, val in NeededArguments.items() if not(val == None or val == False)]

    # Data-points that will be written to the terminal or save file
    # These are found from the Outputs dictionary by comparing with the Wanted_Values list
    ArgumentsToValues = {key: val for key, val in Outputs.items() if key in WantedValues}

    Values = []
    flatten_list([val for val in ArgumentsToValues.values()], Values)

    # How many files to run the script on
    Count = len(InputFiles)

    # If multiprocessing is enabled it will be run using half of the available CPUS
    # Else they will be run in a linear fashion
    if ProgressBar:
        max_filename_length = len(max(InputFiles, key=len))
        TerminalOutput = TerminalInformation(Count, max_filename_length)
        TerminalOutput.start_timer()
    if Multiprocessing:
        with Pool(int(cpu_count()/2)) as pool:
            ExtractedValues = []
            for i, result in enumerate(pool.imap(partial(Data_Extraction, Needed_Values=NeededValues, quiet=Quiet, Temperature=T), InputFiles), start=1):
                if ProgressBar:
                    TerminalOutput.updateProgressbar(i, False, True)
                ExtractedValues.append(result)
            ExtractedValues = {key: value for dictionary in ExtractedValues for key, value in dictionary.items()} # Reformatting Extracted_values
    else:
        ExtractedValues = dict()
        for i, file in enumerate(InputFiles, start=1):
            if ProgressBar:
                TerminalOutput.updateProgressbar(i, True, True, filename=file)
            ExtractedValues[file] = Data_Extraction(file, NeededValues, Quiet, T)[file]

    # Creating Input_Array where all values are put in lists
    InputArray = [[i] for i in ExtractedValues]

    # Checking if some functions have not been implemented for the relevant extraction types
    Check_if_Implemented(InputFiles, ArgumentsToValues, ExtractedValues)

    # If the CPU time has been requested in second or hours instead of minutes
    # then the calculation from minutes to the requested unit is done here
    if RequestedArguments['_CPUS'] == 's':
        for file in InputFiles:
            if ExtractedValues[file]['total_cpu_time'] != ['Not Implemented'] or ExtractedValues[file]['total_cpu_time'] != ['Nan']:
                continue
            ExtractedValues[file]['total_cpu_time'] *= 60
            ExtractedValues[file]['wall_cpu_time'] *= 60
    elif NeededArguments['_CPUS'] == 'h':
        for file in InputFiles:
            if ExtractedValues[file]['total_cpu_time'] != ['Not Implemented'] or ExtractedValues[file]['total_cpu_time'] != ['Nan']:
                continue
            ExtractedValues[file]['total_cpu_time'] /= 60
            ExtractedValues[file]['wall_cpu_time'] /= 60

    # This is done purely for the unittest script to work correctly
    if UnitTesting:
        SaveDict = copy.deepcopy(ExtractedValues)
        for key_outer, dictionary in ExtractedValues.items():
            for key_inner in dictionary:
                if key_inner not in Values:
                    SaveDict[key_outer].pop(key_inner)
        return SaveDict

    # If something is at this point not in a list somehow they will be after this
    for OuterKey, Dict in ExtractedValues.items():
        for InnerKey, Value in Dict.items():
            if not isinstance(Value, list):
                ExtractedValues[OuterKey][InnerKey] = [Value]

    # Collecting all values in arrays in a dictionary
    # Some values in the Extracted_Values dictionary may not have been requested, so these are removed here
    FinalArrays = Collect_and_sort_data(InputFiles, ArgumentsToValues, ExtractedValues)

    # Resizing arrays
    # An example is Excitation energies where there may be more of them in one output file than another
    # By doing this it fits properly in what is printed to the terminal
    for key in FinalArrays:
        if isinstance(FinalArrays[key], list):
            Resize(FinalArrays[key])

    # Fixing the size of variable size arrays so that they match what was requested
    Downsizing_variable_arrays(Outputs, VariableArrays, Count, FinalArrays)
    Upsizing_variable_arrays(Outputs, VariableArrays, Count, FinalArrays, RequestedArguments)

    # Creation of header row
    Header = Create_Header(HeaderText, ArgumentsToValues, FinalArrays)

    # Create output array from header row
    OutputArray = np.array([Header] * (Count + 1), dtype=object)

    # Filling the output array with the extracted data
    Fill_output_array(ArgumentsToValues, InputArray, Count, FinalArrays, OutputArray)

#   ------------ IF CHOSEN PRINTS THE OUTPUT IN A CSV FILE ------------
#   ---------- ELSE THE RESULTS ARE DUMPED INTO THE TERMINAL ----------

    if ProgressBar:
        print("")

    # If this statement is true, then only the filenames have been written to the Output_Array
    if len(OutputArray) == OutputArray.size:
       print("No data was extracted, therefore nothing more will be printed")
       return

    elif Save == 'return':
        SaveDict = copy.deepcopy(ExtractedValues)
        for key_outer, dictionary in ExtractedValues.items():
            for key_inner in dictionary:
                if key_inner not in Values:
                    SaveDict[key_outer].pop(key_inner)
        return SaveDict

    elif Save == 'csv':
        np.savetxt(f'{SaveName}.csv', OutputArray, delimiter=',', fmt='%s')
        print(f'Data has been saved in {SaveName}.csv')
        return

    elif Save == 'npz':
        SaveDict = {i[0]: i[1:] for i in OutputArray}
        np.savez(f'{SaveName}.npz', **SaveDict)
        print(f'Data has been saved in {SaveName}.npz')
        return

    elif Save == 'json':
        SaveDict = copy.deepcopy(ExtractedValues)
        for key_outer, dictionary in ExtractedValues.items():
            for key_inner in dictionary:
                if key_inner not in Values:
                    SaveDict[key_outer].pop(key_inner)
        json_object = json.dumps(SaveDict, indent=4)

        with open(f"{SaveName}.json", "w") as outfile:
            outfile.write(json_object)
        print(f"Data has been saved in {SaveName}.json")

        return

    print(OutputArray)

def main():
    #---------------------------
    # Creating main parser
    #---------------------------
    Parser = argparse.ArgumentParser(formatter_class=argparse.RawDescriptionHelpFormatter, description=f'''
    A script designed to make it easier to extract data from output files

    To extract data you have to use the keyword extract before any other keywords
    To make spectra you will have to use the keyword spectra before any other keywords

            Currently the output formats supported are
            ------------------------------------------
                -  ORCA
                -  DALTON
                -  GAUSSIAN
                -  LSDALTON
                -  VELOXCHEM
                -  AMSTERDAM MODELING SUITE
''', epilog=f'''
For help contact
    Theo Juncker von Buchwald
    fnc970@alumni.ku.dk

    Magnus Bukhave Johansen
    qhw298@alumni.ku.dk
''')

    #---------------------------
    # Creating spectra subparser
    #---------------------------
    subparser = Parser.add_subparsers(dest='pars')
    subparser.required = False

    SpectraSubparser = subparser.add_parser('spectra', formatter_class=argparse.RawDescriptionHelpFormatter, description=f'''
    This part of the script is used for making spectra from data collected from output files

            Currently the output formats supported are
            ------------------------------------------
                -  ORCA
                -  DALTON
                -  GAUSSIAN
                -  LSDALTON
                -  VELOXCHEM
                -  AMSTERDAM MODELING SUITE

    It is currently possible to make UVVIS spectra using excitation energies and complex propagator theory

    The following is not implemented for ORCA
    -  UVVIS based on complex propagator theory

    The following is not implemented for DALTON
    -  UVVIS based on complex propagator theory

    The following is not implemented for GAUSSIAN
    -  UVVIS based on complex propagator theory

    The following is not implemented for VELOXCHEM
    -  UVVIS based on excitation energies
    -  UVVIS based on complex propagator theory

    The following is not implemented for AMSTERDAM MODELING SUITE
    -  UVVIS based on excitation energies
    -  UVVIS based on complex propagator theory
''', help='Use to make spectra such as UVVIS from excitation energies or complex propagator theory')

    # Setting the Spectra function to be run if spectra is used
    SpectraSubparser.set_defaults(func=Spectra)

    # Adding arguments
    SpectraSubparser.add_argument('infile', type=str, nargs='+', help='The file(s) to make spectra from', metavar='File')

    SpectraGroup = SpectraSubparser.add_mutually_exclusive_group()
    SpectraGroup.add_argument('--uvvis', action='store_true', help='Include to make UVVIS spectra')
    SpectraGroup.add_argument('--complex-propagator', action='store_true', help='Include to use complex propagator theory to make the spectra')

    SpectraDataProcessingGroup = SpectraSubparser.add_argument_group('Data processing commands')
    SpectraDataProcessingGroup.add_argument('--format', default='png', const='png', type=str, help='Include this to change the picture format. Will use png as default. If \'--save\' is used together with this, the processed data will be saved in a .npz file', nargs='?', choices=['png', 'eps', 'pdf', 'svg', 'ps'])
    SpectraDataProcessingGroup.add_argument('-s', '--save', action='store_true', help='Saves extracted and processed data in a npz file')

    SpectraAdditionalCommandsGroup = SpectraSubparser.add_argument_group('Additional commands')
    SpectraAdditionalCommandsGroup.add_argument('-q', '--quiet', action='store_true', help='Include for the script to stay silent - This will not remove error messages or the printing of data')
    SpectraAdditionalCommandsGroup.add_argument('-mp','--multiprocessing', action='store_true', help='Include to use the multiprocessing library for data extraction')

    #---------------------------
    # Creating extract subparser
    #---------------------------
    ExtractionSubparser = subparser.add_parser('extract', formatter_class=argparse.RawDescriptionHelpFormatter, description=f'''
    This part of the script is for extracting data from output files and either printing it in the terminal or saving it to a file

            This script can currently extract the data
            ------------------------------------------
                -  Total energies
                -  Zero-Point Vibrational energies
                -  Enthalpies
                -  Entropies
                -  Gibbs Free energies
                -  Dipole moments
                -  Polarizability
                -  Excitation energies
                -  Oscillator strengths
                -  Frequencies
                -  Partition functions at a given temperature
                -  CPU time used
                -  Optimized geometries (or last geometry in file)

    Though not all data types have been implemented for all of the output formats

    All values are extracted as Atomic Units where applicable

    The following is not implemented for LSDALTON
    -  Zero-Point Vibrational energies
    -  Enthalpies
    -  Entropies
    -  Gibbs Free energies
    -  Frequencies
    -  Partition functions

    The following is not implemented for VELOXCHEM
    -  Zero-Point Vibrational energies
    -  Enthalpies
    -  Entropies
    -  Gibbs Free energies
    -  Dipole moments
    -  Polarizability
    -  Excitation energies
    -  Oscillator strengths
    -  Frequencies
    -  Partition functions at a given temperature
    -  CPU time used

    The following is not implemented for AMSTERDAM MODELING SUITE
    -  Zero-Point Vibrational energies
    -  Enthalpies
    -  Entropies
    -  Gibbs Free energies
    -  Dipole moments
    -  Polarizability
    -  Excitation energies
    -  Oscillator strengths
    -  Frequencies
    -  Partition functions at a given temperature
    -  CPU time used
'''
    , help='Use to extract data from output files')

    # Setting the Extract function to be run if extract is used
    ExtractionSubparser.set_defaults(func=Extract)

    # Adding arguments
    ExtractionSubparser.add_argument('infile', type=str, nargs='+', help='The file(s) to extract data from', metavar='File')

    ExtractionGroup = ExtractionSubparser.add_argument_group('Data extraction commands')
    ExtractionGroup.add_argument('-E', '--energy', action='store_true', help='Include to extract the Total Energy')
    ExtractionGroup.add_argument('-Z', '--zpv', action='store_true', help='Include to extract the Zero-Point Vibrational Energy')
    ExtractionGroup.add_argument('-H', '--enthalpy', action='store_true', help='Include to calculate molar enthalpies (kJ/mol)')
    ExtractionGroup.add_argument('-S', '--entropy', action='store_true', help='Include to calculate molar entropes (kJ/(mol*K)')
    ExtractionGroup.add_argument('-G', '--gibbs', action='store_true', help='Include to extract the Gibbs Free Energy')
    ExtractionGroup.add_argument('-D', '--dipole', action='store_true', help='Include to extract the Dipole Moment')
    ExtractionGroup.add_argument('-P', '--polar', action='store_true', help='Include to extract the Polarizability')
    ExtractionGroup.add_argument('-X', '--exc', const=-1, type=int, help='Include to extract the Excitation Energies. Add a number to extract that amount of Excitation Energies. It will extract all Excitation energies as default',nargs='?')
    ExtractionGroup.add_argument('-O', '--osc', action='store_true', help='Include to extract the Oscillator Strengths')
    ExtractionGroup.add_argument('-F', '--freq', const=-1, type=int, help='Include to extract the Frequencies. Add a number to extract that amount of Frequencies. It will extract all Frequencies as default', nargs='?')
    ExtractionGroup.add_argument('-Q', '--partfunc', action='store_true', help='Include to calculate molar partition functions.')
    ExtractionGroup.add_argument('-T', '--temp', const=298.15, default=298.15, type=float, help='Include to calculate at a different temperature. Default is 298.15 K', nargs='?')
    ExtractionGroup.add_argument('-C', '--cpu_time', const=['m'], help='Include to extract total cpu time and pr. cpu time. You can change the output from being in seconds, minutes and hours, where the default is minutes', nargs='?', choices=['s', 'm', 'h'])
    ExtractionGroup.add_argument('-geom', '--optgeom', action='store_true',help='Include to extract optimized geometries and save to \'filename_opt.xyz\'.')

    ExtractionDataProcessingGroup = ExtractionSubparser.add_argument_group('Data processing commands')
    ExtractionDataProcessingGroup.add_argument('-s', '--save', const='csv', type=str, help='Saves extracted and processed data. The extracted data is by default saved in a csv file', nargs='?', choices=['csv', 'npz', 'json', 'return'])
    ExtractionDataProcessingGroup.add_argument('--name', default='data', const='data', type=str, help='Define the name of the datafile where the extracted data is stored', nargs='?', dest='savename')

    ExtractionAdditionalCommandsGroup = ExtractionSubparser.add_argument_group('Additional commands')
    ExtractionAdditionalCommandsGroup.add_argument('-q', '--quiet', '--no-log', action='store_true', help="Include to not print error messages to the 'collect_data.log' file", dest='quiet')
    ExtractionAdditionalCommandsGroup.add_argument('-mp','--multiprocessing', action='store_true', help='Include to use the multiprocessing library for data extraction')
    ExtractionAdditionalCommandsGroup.add_argument('--no-progressbar', action='store_false', help='Include to deactivate progress bar', dest='progressbar')
    ExtractionAdditionalCommandsGroup.add_argument('--unittest', action='store_true', help=argparse.SUPPRESS)

    # Parses the arguments
    args = Parser.parse_args()

    # The arguments are sent to the correct function
    # The function may be one of Spectra, Extract, ...
    args.func(args)

if __name__ == "__main__":
    main()

ext: py

max_stars_repo_name: TheoBuchwald/UCPH-KVM

max_issues_repo_name: TheoBuchwald/UCPH-KVM

max_forks_repo_name: TheoBuchwald/UCPH-KVM

reason: import numpy


text:
import os
import sys

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd
import nevergrad as ng
from tqdm.notebook import tqdm

from src.graph import Node, Edge, Graph
from src.processing import process_embedding_from_words

class NoPrint:
    def __enter__(self):
        self._original_stdout = sys.stdout
        sys.stdout = open(os.devnull, 'w')

    def __exit__(self, exc_type, exc_val, exc_tb):
        sys.stdout.close()
        sys.stdout = self._original_stdout

def phi_edge(i, j, w, embedding, y=1):
    """
    i and j position of the pair of world in df_cross_embedding
    y = y_ij 1 or 0
    w = numpy array

    Interpretation physique: ?
    """
    if y:
        return np.exp(np.dot(w, embedding[i, j]))
    else:
        return 1

def build_graph(w, words, embedding):
    dic_nodes = {}
    edge_names = []

    # T is a node with special distribution
    # TODO : CHeck TUTTE 1984
    dic_nodes["T"] = Node(name="T", node_type="F", dist="MTT")

    for i, wordi in enumerate(words):
        for j, wordj in enumerate(words):
            # TODO: Rajouter les siblings
            if i == j:
                continue
            name_V = f"V_{wordi}_{wordj}"
            node_V = Node(name_V, "V", nb_states=2)
            dic_nodes[name_V] = node_V
            dist = np.array([phi_edge(i, j, w, embedding, y) for y in [0, 1]])
            name_F = f"F_{wordi}_{wordj}"
            node_F = Node(name_F, "F", dist=dist, dist_index=[name_V])
            dic_nodes[name_F] = node_F
            edge_names.append((name_V, name_F))
            edge_names.append((name_V, 'T'))

    graph = Graph(dic_nodes, edge_names)
    return graph

def compute_belief_from_w(w, words, embedding):
    graph = build_graph(w, words, embedding)
    graph.belief_propagation()
    belief = graph.compute_belief_nodes()
    return belief

def build_state(words, df_graph):
    state = {}
    for i, wordi in enumerate(words):
        for j, wordj in enumerate(words):
            # TODO: Rajouter les siblings
            if i == j:
                continue
            name_V = f"V_{wordi}_{wordj}"
            state[name_V] = df_graph.iloc[i, j]
    return state

def proba_state(w, state, words, embedding, bal_coef=1):
    graph = build_graph(w, words, embedding)
    graph.belief_propagation()
    graph.compute_belief_nodes()
    return graph.proba_state(state, bal_coef=bal_coef)

def gradient_from_w(w, state, words, embedding, eps=1e-3, verbose=True, bal_coef=1):
    """
    Compute the gradient of the log max_likelihood with respect to w.
    output:
        grad: np.array
        likelihood: float
    """
    grad = np.zeros_like(w, dtype=np.float32)
    likelihood = proba_state(w, state, words, embedding, bal_coef=bal_coef)

    for i, _ in enumerate(tqdm(grad, total=len(w), disable=not verbose, leave=False)):
        w_eps = w.copy()
        w_eps[i] += eps
        likelihood_eps = proba_state(w_eps, state, words, embedding, bal_coef=bal_coef)
        grad[i] = (likelihood_eps - likelihood)/eps
    return grad, likelihood

def plot_state_likelihood_and_gradient(state, words, embedding, W_X, W_Y, bal_coef=1):
    # plot offset embedding[0] and feature embedding[1]
    embedding = embedding[:, :, :2].copy()

    @np.vectorize
    def likelihood_from_coordinates(*w_i):
        w = np.array(w_i)
        likelihood = proba_state(w, state, words, embedding, bal_coef=bal_coef)
        return likelihood

    @np.vectorize
    def norm_of_gradient_from_coordinates(*w_i):
        w = np.array(w_i)
        grad, _ = gradient_from_w(w, state, words, embedding, verbose=False, bal_coef=bal_coef)
        return np.linalg.norm(grad)

    def plot_function(W_X, W_Y, fig, ax, func):
        if func == "likelihood":
            L = likelihood_from_coordinates(W_X, W_Y)
            title = 'Log-likelihood of the state as a function of w'
        elif func == "gradient":
            L = norm_of_gradient_from_coordinates(W_X, W_Y)
            title = 'Gradient of the log-likelihood of the graph'
        else:
            raise ValueError(
                "func: {} not in ['likelihood', 'gradient'].".format(func))
        L = L[:-1, :-1]
        left_x, right_x = W_X.min(), W_X.max()
        left_y, right_y = W_Y.min(), W_Y.max()
        left_l, right_l = -np.abs(L).max(), np.abs(L).max()
        left_l, right_l = L.min(), L.max()
        L = ax.pcolormesh(W_X, W_Y, L, cmap='hot', vmin=left_l, vmax=right_l)
        ax.set_title(title, fontsize=16)
        ax.set_xlabel("$W_0$: offset", fontsize=14)
        ax.set_ylabel("$W_1$: feature", fontsize=14)
        ax.axis([left_x, right_x, left_y, right_y])
        fig.colorbar(L, ax=ax)

    figure, axes = plt.subplots(1, 2, figsize=(14, 6))
    plot_function(W_X, W_Y, figure, axes[0], "likelihood")
    plot_function(W_X, W_Y, figure, axes[1], "gradient")
    plt.tight_layout()
    plt.show()

def plot_nx_graph(w, words, embedding, verbose=False):
    G = nx.Graph()
    G.add_nodes_from(words)
    belief = compute_belief_from_w(w, words, embedding)
    for i, wordi in enumerate(words):
        for j, wordj in enumerate(words):
            if i == j:
                continue
            name_V = f"V_{wordi}_{wordj}"
            proba_link = np.exp(belief[name_V][1])
            if verbose:
                print(name_V, proba_link)
            if proba_link > 0.5:
                G.add_edge(wordi, wordj)
    nx.draw(G, with_labels=True, pos=nx.spring_layout(G))
    return G

def check_process_embedding_from_words(embedd, words, verbose=False):
    if embedd is None:
        return process_embedding_from_words(words)
    else:
        if verbose:
            print("Embedding already charged")
        return embedd

class TaxonomyModule:
    def __init__(self, nb_features):
        self.w = np.random.randn(nb_features)   # already includes the offset
        self.train_embedd = None
        self.test_embedd = None
        self.nb_features = nb_features
        self.losses = []

    def train(self, train_df_graph, nb_epochs, lr, bal_coef=1, verbose=True):
        """
        if lr = None : uses nevergrad
        """
        train_words = list(train_df_graph.columns)

        with NoPrint():
            self.train_embedd = check_process_embedding_from_words(
                self.train_embedd, train_words)

        train_state = build_state(train_words, train_df_graph)

        if lr is not None:
            self.losses = []
            for epoch in tqdm(range(nb_epochs), disable=not verbose):
                grad, likelihood = gradient_from_w(
                    self.w, train_state, train_words, self.train_embedd,
                    bal_coef=bal_coef, verbose=verbose)
                self.w += lr*grad
                self.losses.append(likelihood)
                if verbose:
                    print("[Epoch {}]: {:.3f}".format(epoch + 1, likelihood))
        else:
            def f(w):
                return - proba_state(
                    w, train_state, train_words, self.train_embedd, bal_coef=bal_coef)
            with NoPrint():
                optimizer = ng.optimizers.NGOpt(
                    parametrization=self.nb_features, budget=nb_epochs, num_workers=4)
            try:
                verbosity = 2 if verbose else 0
                recommendation = optimizer.minimize(f, verbosity=verbosity)
                self.w = recommendation.value
            except KeyboardInterrupt:
                self.w = optimizer.provide_recommendation()

    def f1_score(self, test_df_graph, eps=1e-8):
        test_words = list(test_df_graph.columns)

        with NoPrint():
            self.test_embedd = check_process_embedding_from_words(
                self.test_embedd, test_words)

        belief = compute_belief_from_w(self.w, test_words, self.test_embedd)
        truth_graph = np.array(test_df_graph, dtype=bool)
        inferred_graph = np.zeros(
            (len(test_words), len(test_words)), dtype=bool)
        for i, wordi in enumerate(test_words):
            for j, wordj in enumerate(test_words):
                if i == j:
                    continue
                name_V = f"V_{wordi}_{wordj}"
                proba_link = np.exp(belief[name_V][1])
                if proba_link > 0.5:
                    inferred_graph[i, j] = 1
        tp = np.sum((inferred_graph == 1) & (truth_graph == 1))
        fp = np.sum((inferred_graph == 1) & (truth_graph == 0))
        fn = np.sum((inferred_graph == 0) & (truth_graph == 1))
        precision = tp/(tp+fp+eps)
        recall = tp/(tp+fn+eps)
        f1_score = 2 * (precision * recall) / (precision + recall + eps)
        return f1_score

    def infer_graph(self, test, verbose=False):
        """
        input: either a dataframe df_graph or a list of words
        verbose : output probabilities of each link
        """
        if type(test) is pd.DataFrame:
            test_words = list(test.columns)
        elif type(test) is list:
            test_words = test
        else:
            raise TypeError(
                "input must be either a dataframe or a list of words")

        with NoPrint():
            self.test_embedd = check_process_embedding_from_words(
                self.test_embedd, test_words)

        return plot_nx_graph(self.w, test_words, self.test_embedd, verbose=verbose)

max_stars_repo_name: clement-bonnet/belief-propagation

max_issues_repo_name: clement-bonnet/belief-propagation

max_forks_repo_name: clement-bonnet/belief-propagation

reason: import numpy,import networkx


text:
#%%

import cvxpy as cp
import numpy as np
import matplotlib.pyplot as plt

def loss_fn(X, Y, beta):
    return cp.norm2(cp.matmul(X, beta) - Y)**2

def regularizer(beta):
    return cp.norm1(beta)

def objective_fn(X, Y, beta, lambd):
    return loss_fn(X, Y, beta) + lambd * regularizer(beta)

def mse(X, Y, beta):
    return (1.0 / X.shape[0]) * loss_fn(X, Y, beta).value

def generate_data(m=100, n=20, sigma=5, density=0.2):
    "Generates data matrix X and observations Y."
    np.random.seed(1)
    beta_star = np.random.randn(n)
    idxs = np.random.choice(range(n), int((1-density)*n), replace=False)
    for idx in idxs:
        beta_star[idx] = 0
    X = np.random.randn(m,n)
    Y = X.dot(beta_star) + np.random.normal(0, sigma, size=m)
    return X, Y, beta_star

m = 100
n = 20
sigma = 5
density = 0.2

X, Y, _ = generate_data(m, n, sigma)
X_train = X[:50, :]
Y_train = Y[:50]
X_test = X[50:, :]
Y_test = Y[50:]

beta = cp.Variable(n)
lambd = cp.Parameter(nonneg=True)
problem = cp.Problem(cp.Minimize(objective_fn(X_train, Y_train, beta, lambd)))

lambd_values = np.logspace(-2, 3, 50)
train_errors = []
test_errors = []
beta_values = []
for v in lambd_values:
    lambd.value = v
    problem.solve()
    train_errors.append(mse(X_train, Y_train, beta))
    test_errors.append(mse(X_test, Y_test, beta))
    beta_values.append(beta.value)

# matplotlib inline
# config InlineBackend.figure_format = 'svg'

def plot_train_test_errors(train_errors, test_errors, lambd_values):
    plt.plot(lambd_values, train_errors, label="Train error")
    plt.plot(lambd_values, test_errors, label="Test error")
    plt.xscale("log")
    plt.legend(loc="upper left")
    plt.xlabel(r"$\lambda$", fontsize=16)
    plt.title("Mean Squared Error (MSE)")
    plt.show()

plot_train_test_errors(train_errors, test_errors, lambd_values)

print('done')

reason: import numpy,import cvxpy


text:
import numpy as np
from torch.nn import functional as F
from gym import Env, spaces
from dataset import get_dataset
from merge_strategy import strategy

class TimeSeriesDetector(Env):
    def __init__(self, config: dict):
        super(TimeSeriesDetector, self).__init__()
        self.observation_length = config['observation_length']
        self.num_channels = config['channels']
        self.action_stride = config.get('action_stride', 0)
        num_actions = config['num_event_classes']+1  # abstain
        self.action_space = spaces.Discrete(num_actions)
        self.observation_space = spaces.Box(
            low=config['low'],
            high=config['high'],
            shape=(self.num_channels, self.observation_length),
            dtype=np.float32,
        )
        # Load the data
        self.ds = get_dataset(**config['data'])

    def step(self, action: int):
        reward = 0.0
        y = self.y[:, self.current_step:self.current_step +
                   self.observation_length]
        #loss = F.nll_loss(prediction, target)

        self.current_step += 1
        if action != 0:
            self.current_step += self.action_stride
        obs = self.get_observation()
        done = self.current_step >= self.x.shape[1] - self.observation_length
        info = self.get_info_dict()
        return obs, reward, done, info

    def reset(self):
        self.current_step = 0
        i = np.random.randint(0, len(self.ds))
        self.x, self.y = self.ds[i]
        if self.x.shape[1] < self.observation_length:
            raise ValueError(f'Example {i} is shorter ({self.x.shape[1]}) '
                             f'than the observation length ({self.observation_length})')
        return self.get_observation()

    def get_observation(self):
        return self.x[:, self.current_step:self.current_step +
                      self.observation_length].numpy()

    def get_info_dict(self) -> dict:
        return dict(current_step=self.current_step)

envs = {
    'TimeSeriesDetector': TimeSeriesDetector,
}

def get_env(name: str):
    if name not in envs:
        raise ValueError(f'Environment "{name}" not found '
                         f'valid options are {envs}')
    return envs[name]

max_stars_repo_name: thavlik/machine-learning-portfolio

max_issues_repo_name: thavlik/machine-learning-portfolio

max_forks_repo_name: thavlik/machine-learning-portfolio


text:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Feb 13 19:05:23 2018
The DV flux time series data is stored in h5d format after resampling
at a coarser time sampling.  Also, some data
about the target from the header of the fits file is also kept

@author: Christopher J. Burke
"""

import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt
import h5py
import glob
import os
import math

def make_data_dirs(prefix, sector, epic):
    secDir = 'S{0:02d}'.format(sector)
    localDir = os.path.join(prefix,secDir)
    if not os.path.exists(localDir):
        os.mkdir(localDir)
    epcDir = '{0:04d}'.format(int(math.floor(epic/1000.0)))
    localDir = os.path.join(prefix,secDir,epcDir)
    if not os.path.exists(localDir):
        os.mkdir(localDir)
    return localDir

def idx_filter(idx, *array_list):
    new_array_list = []
    for array in array_list:
        new_array_list.append(array[idx])
    return new_array_list

def dvts_resamp(file, dirOut, RESAMP, SECTOR=None, overwrite=True):
    """ Resample TESS dv time series file and save as h5d format
        resamp - Resample factor just make it odd okay"""

    hdulist = fits.open(file)
    prihdr = hdulist[0].header
    nTces = prihdr['NUMTCES']

    dataSpanMax = 0.0
    # Get header information that we should keep
    if not SECTOR is None:
        keepprihdr = ['TICID','RA_OBJ', \
                  'DEC_OBJ','PMRA','PMDEC','PMTOTAL','TESSMAG','TEFF', \
                  'LOGG','RADIUS']
        formatprihdr = [np.uint32, \
                    np.float,np.float,np.float,np.float,np.float, \
                    np.float,np.float,np.float,np.float,np.float]
    else:
        keepprihdr = ['TICID','SECTOR','PXTABLE','RA_OBJ', \
                      'DEC_OBJ','PMRA','PMDEC','PMTOTAL','TESSMAG','TEFF', \
                      'LOGG','RADIUS']
        formatprihdr = [np.uint32, np.int, np.int, \
                        np.float,np.float,np.float,np.float,np.float, \
                        np.float,np.float,np.float,np.float,np.float]

    # make empty arrays
    kpCadenceNo = np.array([0], dtype=np.int)
    kpTimetbjd = np.array([0.0], dtype=np.float)
    kpQuality = np.array([0], dtype=np.int)
    kpPDC = np.array([0.0], dtype=np.float)
    for ii in range(nTces):
        # Check if already done
        epic = hdulist[0].header['TICID']
        if SECTOR is None:
            sec = hdulist[0].header['SECTOR']
        else:
            sec = SECTOR
        pn = ii+1
        fileoutput = os.path.join(make_data_dirs(dirOut,sec,epic), 'tess_dvts_{0:016d}_{1:02d}.h5d'.format(epic,pn))
        fileExists=os.path.isfile(fileoutput)
        if (not fileExists) or overwrite:

            extname = 'TCE_{0:d}'.format(ii+1)
            arr = hdulist[extname].data['TIME']
            nImage = len(arr)

            cadenceNo = hdulist[extname].data['CADENCENO']
            kpCadenceNo = cadenceNo
            timetbjd = hdulist[extname].data['TIME']
            kpTimetbjd = timetbjd
            lc_init = hdulist[extname].data['LC_INIT']
            lc_init_err = hdulist[extname].data['LC_INIT_ERR']
            lc_white = hdulist[extname].data['LC_WHITE']
            lc_med_detrend = hdulist[extname].data['LC_DETREND']
            lc_model = hdulist[extname].data['MODEL_INIT']
            lc_white_model = hdulist[extname].data['MODEL_WHITE']
            lc_phase = hdulist[extname].data['PHASE']
            pdc_flux = hdulist['statistics'].data['PDCSAP_FLUX']
            pdc_flux_err = hdulist['statistics'].data['PDCSAP_FLUX_ERR']
            deweights = hdulist['statistics'].data['DEWEIGHTS']
            kpQuality = hdulist['statistics'].data['QUALITY']
            kpPDC = hdulist['statistics'].data['PDCSAP_FLUX']

            # Fix the issue of having times close to zero
            # First get time stamps on valid data
            idx = np.where((np.isfinite(timetbjd)) & (np.isfinite(lc_init)) & (np.isfinite(pdc_flux)))[0]
            # minimum time on valid data
            minGdTime = np.min(timetbjd[idx])
            idx = np.where(timetbjd < minGdTime-1.0)[0]
            timetbjd[idx] = np.nan

            newNImage = int(np.floor(nImage / RESAMP))
            oldNImage = newNImage*RESAMP
            # trim off the excess images not integral into resamp
            idx = np.arange(0,oldNImage)
            cadenceNo, timetbjd, lc_init, lc_init_err, lc_white, lc_med_detrend, \
                lc_model, lc_white_model, lc_phase, pdc_flux, pdc_flux_err, \
                deweights = idx_filter(idx, cadenceNo, timetbjd, lc_init, lc_init_err, lc_white, lc_med_detrend, \
                                       lc_model, lc_white_model, lc_phase, pdc_flux, pdc_flux_err, \
                                       deweights)

            # Do downsampling of data stream
            cadenceNoBeg = np.min(np.reshape(cadenceNo, (newNImage, RESAMP)), axis=1)
            cadenceNoEnd = np.max(np.reshape(cadenceNo, (newNImage, RESAMP)), axis=1)
            cadenceNo = np.mean(np.reshape(cadenceNo, (newNImage, RESAMP)), axis=1, dtype=np.int)
            timetbjd = np.mean(np.reshape(timetbjd, (newNImage, RESAMP)), axis=1)
            lc_init = np.mean(np.reshape(lc_init, (newNImage, RESAMP)), axis=1)
            lc_init_err = np.mean(np.reshape(lc_init_err, (newNImage, RESAMP)), axis=1)
            lc_white = np.mean(np.reshape(lc_white, (newNImage, RESAMP)), axis=1)
            lc_med_detrend = np.mean(np.reshape(lc_med_detrend, (newNImage, RESAMP)), axis=1)
            lc_model = np.mean(np.reshape(lc_model, (newNImage, RESAMP)), axis=1)
            lc_white_model = np.mean(np.reshape(lc_white_model, (newNImage, RESAMP)), axis=1)
            lc_phase = np.mean(np.reshape(lc_phase, (newNImage, RESAMP)), axis=1)
            pdc_flux = np.sum(np.reshape(pdc_flux, (newNImage, RESAMP)), axis=1)
            pdc_flux_err = np.mean(np.reshape(pdc_flux_err, (newNImage, RESAMP)), axis=1)
            deweights = np.mean(np.reshape(deweights, (newNImage, RESAMP)), axis=1)

            # Identify data that is missing or NaN
            idx = np.where((np.isfinite(timetbjd)) & (np.isfinite(lc_init)) & (np.isfinite(pdc_flux)))[0]
            valid_data_flag = np.zeros((newNImage,), dtype=np.bool_)
            valid_data_flag[idx] = True

            # Trim all leading in-valid data
            if not valid_data_flag[0]:
                idx = np.where(valid_data_flag)[0]
                if not len(idx) == 0:
                    idx = idx[0]
                    cadenceNoBeg = cadenceNoBeg[idx:]
                    cadenceNoEnd = cadenceNoEnd[idx:]
                    cadenceNo = cadenceNo[idx:]
                    timetbjd = timetbjd[idx:]
                    lc_init = lc_init[idx:]
                    lc_init_err = lc_init_err[idx:]
                    lc_white = lc_white[idx:]
                    lc_med_detrend = lc_med_detrend[idx:]
                    lc_model = lc_model[idx:]
                    lc_white_model = lc_white_model[idx:]
                    lc_phase = lc_phase[idx:]
                    pdc_flux = pdc_flux[idx:]
                    pdc_flux_err = pdc_flux_err[idx:]
                    deweights = deweights[idx:]
                    valid_data_flag = valid_data_flag[idx:]
                else:
                    print('No Valid data? {0:d} {1:d}'.format(ii,epic))

            # Get the data span in days for use in federation steps
            idx = np.where(valid_data_flag == True)[0]
            if (len(idx)>0):
                dataSpan = np.max(timetbjd[idx]) - np.min(timetbjd[idx])
                if dataSpan > dataSpanMax:
                    dataSpanMax = dataSpan

            # Now save data as h5py
            f = h5py.File(fileoutput, 'w')
            tmp = f.create_dataset('cadenceNo', data=cadenceNo, compression='gzip')
            tmp = f.create_dataset('cadenceNoBeg', data=cadenceNoBeg, compression='gzip')
            tmp = f.create_dataset('cadenceNoEnd', data=cadenceNoEnd, compression='gzip')
            tmp = f.create_dataset('timetbjd', data=timetbjd, compression='gzip')
            tmp = f.create_dataset('lc_init', data=lc_init, compression='gzip')
            tmp = f.create_dataset('lc_init_err', data=lc_init_err, compression='gzip')
            tmp = f.create_dataset('lc_white', data=lc_white, compression='gzip')
            tmp = f.create_dataset('lc_med_detrend', data=lc_med_detrend, compression='gzip')
            tmp = f.create_dataset('lc_model', data=lc_model, compression='gzip')
            tmp = f.create_dataset('lc_white_model', data=lc_white_model, compression='gzip')
            tmp = f.create_dataset('lc_phase', data=lc_phase, compression='gzip')
            tmp = f.create_dataset('pdc_flux', data=pdc_flux, compression='gzip')
            tmp = f.create_dataset('pdc_flux_err', data=pdc_flux_err, compression='gzip')
            tmp = f.create_dataset('deweights', data=deweights, compression='gzip')
            tmp = f.create_dataset('valid_data_flag', data=valid_data_flag, compression='gzip')
            for i in range(len(keepprihdr)):
                curval = hdulist[0].header[keepprihdr[i]]
                if np.isscalar(curval):
                    tmp = f.create_dataset(keepprihdr[i], data=np.array([hdulist[0].header[keepprihdr[i]]], dtype=formatprihdr[i]))
                else:
                    tmp = f.create_dataset(keepprihdr[i], data=np.array([-1], dtype=formatprihdr[i]))

            f.close()
    return dataSpanMax, kpCadenceNo, kpTimetbjd, kpQuality, kpPDC

if __name__ == "__main__":

    dirInputs = '/pdo/spoc-data/sector-042/dv-time-series/'
    dirOutputs = '/pdo/users/cjburke/spocvet/sector42/'
    RESAMP = 5  ###  USE AN ODD NUMBER HELPS WITH CADENCE NO ###
    SECTOR_OVRRIDE = None # If NOT multisector set this to None ###
    overwrite = True # Set False to keep old results and only do files that dont exist

    fileList = glob.glob(os.path.join(dirInputs, '*dvt.fits*'))

    # For debugging a target filter fileList to wanttic string
    #wanttic='278139637'
    #newFileList = []
    #for fil in fileList:
    #    if fil.find(wanttic) != -1:
    #        newFileList.append(fil)
    #fileList = newFileList
    # END target debugging

    cnt = 0
    # Keep track of max data span for use in federation
    dataSpanMax = 0.0
    gdFracMax = 0.0
    for fil in fileList:
        cnt = cnt + 1
        if np.mod(cnt,10) == 0:
            print(cnt,' Data Span {0:f} GdFrac: {1:f}'.format(dataSpanMax, gdFracMax))
        dataSpan, cadno, timetjd, quality, pdc = dvts_resamp(fil, dirOutputs, RESAMP, SECTOR=SECTOR_OVRRIDE, overwrite=overwrite)
        nData = len(quality)
        gdFrac = float(len(np.where((quality==0) & (np.isfinite(pdc)==True))[0]))/float(nData)
        if dataSpan+2.0*gdFrac > dataSpanMax+2.0*gdFracMax:
            dataSpanMax = dataSpan
            gdFracMax = gdFrac
            fout = open('cadnoVtimemap.txt','w')
            for i, cad in enumerate(cadno):
                momdump = int(0)
                pdcfinite = int(0)
                if quality[i] & 32:
                    momdump = int(1)
                if np.isfinite(pdc[i]):
                    pdcfinite = int(1)
                fout.write('{0:d} {1:f} {2:d} {3:d} {4:d}\n'.format(int(cad), timetjd[i], int(quality[i]), momdump, pdcfinite))
            fout.close()

max_stars_repo_name: christopherburke/TESS-ExoClass

max_issues_repo_name: christopherburke/TESS-ExoClass

max_forks_repo_name: christopherburke/TESS-ExoClass

reason: import numpy,from astropy


text:
#!/usr/bin/env python
# encoding: utf-8
"""
observingPlan.py

Created by José Sánchez-Gallego on 11 Dec 2013.
Copyright (c) 2013. All rights reserved.
Licensed under a 3-clause BSD license.

"""

from __future__ import division
from __future__ import print_function
from astropy import table, time
from Totoro import exceptions
import warnings
from Totoro import config, log, readPath
import numpy as np
import os

def getScheduleFile():

    try:
        # Gets the master scheduler schedule file
        import autoscheduler
        autoschedulerDir = os.path.dirname(autoscheduler.__file__)
        schedule = os.path.join(autoschedulerDir, '../../schedules',
                                'Sch_base.6yrs.txt.frm.dat')
    except:
        # If the autoscheduler is not available, uses the local schedule file
        schedule = readPath(config['observingPlan']['fallBackSchedule'])
        warnings.warn('The master autoscheduler could not be found. '
                      'Using a local copy of the schedule. BE CAREFUL! '
                      'THIS FILE MIGHT BE OUTDATED!',
                      exceptions.TotoroUserWarning)

    return schedule

class ObservingPlan(object):
    """The survey-wide observing plan object.

    Parameters
    ----------
    plan : str or `astropy.tableTable`
        A file pth containing the observing plan or an astropy table
        with the data.
    format : str
        The input format of the plan. The options are 'autoscheduler',
        'autoscheduler.utc', 'autoscheduler.jd', 'sidereal', 'utc' or 'jd'.
    useOptimisedPlan : bool
        It True, tries to read the optimised file (in JD format) for the
        input plan. This notably improves performance.
    saveOptimisedPlan : bool
        If True, once the input observing plan has been converted to JD,
        saves it to the configuration directory. This optimised file will
        later be read if useOptimisedPlan is True.
    survey : str
        If None, the observing plan for all plans will be kept, but survey
        will need to be defined for some methods.

    """

    def __init__(self, schedule=config['observingPlan']['schedule'],
                 **kwargs):

        if isinstance(schedule, basestring) and schedule.lower() == 'none':
            schedule = getScheduleFile()
        else:
            schedule = readPath(schedule)

        self.scheduleFile = os.path.realpath(schedule)
        scheduleToPrint = self.scheduleFile[:15] + '...' + \
            self.scheduleFile[-45:] if len(self.scheduleFile) > 70 else \
            self.scheduleFile

        if not os.path.exists(self.scheduleFile):
            raise exceptions.TotoroError('schedule {0} not found'.format(
                                         scheduleToPrint))

        plan = table.Table.read(self.scheduleFile,
                                format='ascii.no_header')

        self.plan = plan.copy()
        self.plan.keep_columns(['col1', 'col11', 'col12'])
        self.plan.rename_column('col1', 'JD')
        self.plan.rename_column('col11', 'JD0')
        self.plan.rename_column('col12', 'JD1')

        self.plan.add_column(
            table.Column(np.zeros(len(self.plan)), name='Position', dtype=int))
        for ii, row in enumerate(plan):
            if row['col11'] == 0:
                continue
            if row['col11'] > row['col5'] and row['col11'] > row['col9']:
                self.plan['Position'][ii] = 2
            else:
                self.plan['Position'][ii] = 1

        log.debug('observing plan {0} loaded.'.format(scheduleToPrint))

        self.addRunDayCol()
        self.plan = self.plan[(self.plan['JD0'] > 0) & (self.plan['JD1'] > 0)]

    def addRunDayCol(self):
        """Adds a column with the night within the run."""

        nDay = 1
        nRun = 1
        ll = []
        run = []

        for row in self.plan:
            if row['JD0'] != 0.:
                ll.append(nDay)
                run.append(nRun)
                nDay += 1
            else:
                ll.append(-1)
                run.append(-1)
                nDay = 1
                nRun += 1

        run = np.array(run)
        for nn, value in enumerate(np.unique(run[run > 0])):
            run[run == value] = nn + 1

        self.plan.add_column(table.Column(data=run, name='RUN', dtype=int))
        self.plan.add_column(table.Column(data=ll, name='RUN_DAY', dtype=int))

    def getSurveyStart(self):
        """Gets the start of survey."""

        validDates = self.plan[self.plan['JD0'] > 0.0]
        startTime = validDates[0]['JD0']
        return startTime

    def getSurveyEnd(self):
        """Gets the end of survey."""

        validDates = self.plan[self.plan['JD1'] > 0.0]
        endTime = validDates[-1]['JD1']
        return endTime

    def getClosest(self, dd):

        tt = self.plan[(self.plan['JD0'] < dd) & (self.plan['JD1'] > dd)]

        if len(tt) == 1:
            return tt

        tt = self.plan[self.plan['JD'] <= dd]

        return tt[-1]

    def getJD(self, jd=None):

        if jd is None:
            jd = time.Time.now().jd

        jd = int(jd)

        if jd not in self.plan['JD']:
            warnings.warn('JD={0} not found in schedule'.format(jd),
                          exceptions.TotoroUserWarning)
            return (None, None)

        night = self.plan[self.plan['JD'] == jd]

        return (night['JD0'][0], night['JD1'][0])

    def getObservingBlocks(self, startDate, endDate):
        """Returns an astropy table with the observation dates
        for each night between startDate and endDate."""

        validDates = self.plan[(self.plan['JD1'] >= startDate) &
                               (self.plan['JD0'] <= endDate) &
                               (self.plan['JD0'] > 0.0) &
                               (self.plan['JD1'] > 0.0)]

        if (startDate is None and endDate is None) or (len(validDates) == 0):
            log.info('no observing blocks selected.')
            return validDates

        if startDate > validDates['JD0'][0]:
            validDates['JD0'][0] = startDate

        if endDate < validDates['JD1'][-1]:
            validDates['JD1'][-1] = endDate

        totalTime = np.sum([row['JD1'] -
                            row['JD0'] for row in validDates]) * 24

        log.info(('{0} blocks (days) selected, '
                  'making a total of {1:.2f} hours').format(len(validDates),
                                                            totalTime))

        return validDates

    def _createObservingBlock(self, startDate, endDate):
        """Creates an observing block from start and end dates. Mainly for the
        plugger."""

        tmpPlan = self.plan.copy()
        tmpPlan.add_row((int(startDate), startDate, endDate, -1, -1))

        totalTime = (endDate - startDate) * 24.
        log.info(('1 block (days) selected, '
                  'making a total of {0:.2f} hours').format(totalTime))

        return table.Table(tmpPlan[len(tmpPlan) - 1])

    def getRun(self, startDate=None):

        jd = int(startDate if startDate is not None else time.Time.now().jd)

        if jd not in self.plan['JD']:
            raise exceptions.TotoroError(
                'JD={0} not found in schedule'.format(jd))

        firstDay = self.plan[self.plan['JD'] == jd]
        runNumber = firstDay['RUN']
        lastDay = self.plan[self.plan['RUN'] == runNumber][-1]

        return (firstDay['JD0'], lastDay['JD1'])

    def getMJD(self, mjd):

        mjd = int(mjd)

        jd = 2400000. + int(mjd)
        row = self.plan[self.plan['JD'] == jd]

        if len(row) == 1:
            return (row['JD0'][0], row['JD1'][0])
        else:
            return None

    def getPosition(self, jd):
        """Returns the position of the block in which a certain JD falls."""

        try:
            pos = int(self.getClosest(jd)['Position'])
        except:
            warnings.warn('Cannot found block position for JD={0}'.format(jd),
                          exceptions.TotoroUserWarning)
            pos = None

        return pos

    def __repr__(self):
        return self.plan.__repr__()

    def __str__(self):
        return self.plan.__str__()

    def __getitem__(self, slice):
        return self.plan[slice]

max_stars_repo_name: ApachePointObservatory/Totoro

max_issues_repo_name: ApachePointObservatory/Totoro

max_forks_repo_name: albireox/Totoro


text:
import nose.tools as nt
import numpy as np
import theano
import theano.tensor as T

import treeano
import treeano.nodes as tn
from treeano.sandbox.nodes import inverse

fX = theano.config.floatX

def test_inverse_node_serialization():
    tn.check_serialization(inverse.InverseNode("a"))

def test_inverse_node():
    network = tn.SequentialNode(
        "s",
        [tn.InputNode("i", shape=(1, 1, 2, 2)),
         tn.MaxPool2DNode("m", pool_size=(2, 2)),
         tn.InputNode("i2", shape=(1, 1, 1, 1)),
         inverse.InverseNode("in", reference="m")]
    ).network()
    fn = network.function(["i", "i2"], ["in"])
    x = np.array([[[[1, 2],
                    [3, 4]]]],
                 dtype=fX)
    x2 = np.array(np.random.randn(), dtype=fX)
    ans = x2 * np.array([[[[0, 0],
                           [0, 1]]]],
                        dtype=fX)

    np.testing.assert_equal(ans, fn(x, x2.reshape(1, 1, 1, 1))[0])

reason: import numpy,import theano