text:
# This module has been developed by Gabriel Fonseca and Steve Furber at The University of Manchester as part of the
# paper:
#
# "Using Stochastic Spiking Neural Networks on SpiNNaker to Solve Constraint Satisfaction Problems"
# Submitted to the journal Frontiers in Neuroscience| Neuromorphic Engineering
#-----------------------------------------------------------------------------------------------------------------------
"""Implement a framework to map a constraint satisfaction problem into a spiking neural network.

This module contains the CSP class, which stands for Constraint Satisfaction Problem. Its methods allow the creation of
a network of leaky integrate and fire spiking neurons whose connectivity represent the CSP problem, the connections
are either inhibitory or excitatory. The neurons are stochastically stimulated by spike sources implementing a Poisson
process causing the network dynamics to implement a stochastic search of the satisying configuration.
"""
import spynnaker7.pyNN as p  # simulator
from pyNN.random import RandomDistribution
import numpy as np
import os

msg = '%s \n'%('='*70) # a separator for readability of messages on standard output

class CSP:
    """Map a constraint satisfaction problem into a spiking neural network. """
    live = False
    run_time = 30000
    # Lists for counting populations to build report.
    var_pops = []
    stim_pops = [[]]
    diss_pops = []
    n_populations = 1
    d_populations = 1
    # Lists for counting synapses to build report.
    core_conns = []
    internal_conns = []
    stim_conns = []
    diss_conns = []
    stim_times = []
    diss_times = []
    constraint_conns = []
    # Whether set_clues populations should receive inhibition from other sources.
    clues_inhibition = False
    # Parameters for the leaky integrate and fire neurons.
    cell_params_lif = {'cm': 0.25,  # nF          membrane capacitance
                       'i_offset': 0.3,  # nA          bias current
                       'tau_m': 20.0,  # ms          membrane time constant
                       'tau_refrac': 2.0,  # ms          refractory period
                       'tau_syn_E': 5.0,  # ms          excitatory synapse time constant
                       'tau_syn_I': 5.0,  # ms          inhibitory synapse time constant
                       'v_reset': -70.0,  # mV          reset membrane potential
                       'v_rest': -65.0,  # mV          rest membrane potential
                       'v_thresh': -50.0  # mV          firing threshold voltage
                       }

    def __init__(self, variables_number=0, domain_size=0, constraints=[], exc_constraints=[],
                 core_size=25, directed=False, run_time=30000):
        """Initialize the constraint satisfaction problem spiking neural network.

        Args:
            variables_number: how many variables has the CSP (integer).
            domain_size: how many values can assume each variable (integer).
            constraints: a list of tuples of conflicting variables (list of tuples).
            exc_constraints: a list of tuples of variables taking the same value (list of tuples).
            core_size: number of neurons to represent each value that a variable can assume (integer).
            directed: if True applies the constraints  from source to target as defined by the list of tuples,
                if False constraints are applied also from target to source (undirected graph).
        """
        self.variables_number = variables_number
        self.domain_size = domain_size
        self.core_size = core_size
        self.size = domain_size * core_size
        self.directed = directed
        self.constraints = constraints
        self.exc_constraints = exc_constraints
        self.clues = [[]]
        self.run_time = run_time

    def set_clues(self, clues):
        """Take set_clues as an array of the form [[list of variables],[list of values]].

        Here clues are fixed and predetermined values for particular variables, These influence the constraints
        implementation.

        args:
            clues: an array of the form [[list of variable ids],[list of values taken by those variables]]
        """
        self.clues = clues

    def build_domains_pops(self):
        """Generate an array of pyNN population objects, one for each CSP variable.

        The population size will be self.size = domain_size*core_size.
        var_pops[i] is the population for variable i including all domain sub-populations, each of size core_size.
        """
        print msg, 'creating an array of %d neural populations'%(self.variables_number)
        var_pops = []
        for variable in range(self.variables_number):
            var_pops.append(p.Population(self.size,
                                         p.IF_curr_exp, self.cell_params_lif,
                                         label="var%d" % (variable + 1)))
        self.var_pops = var_pops

    def poisson_params(self, n_populations, full=False, stim_ratio=1.0, shrink=1.0, phase=0.0):
        """Define time intervals for activation of the pyNN Poisson noise sources.

        This method defines the temporal dependence of the stimulating noise, it is an internal method called by the
        build_stimulation_pops method. Here we use the word noise to refer to spike sources implementing a random
        Poisson process. In pyNN these objects connect with neurons using synapses as if they were neurons too.
        Currently each SpikeSourcePoisson object accepts only a start time and a duration time, thus to change the noise
        level through time one should create n different populations and activate them at different times.
        Here we uniformly distribute the start times of the n_populations from phase to run_time. Each population will
        be active for a period lapso = shrink * self.run_time / n_populations if full=False otherwise will stay active
        during all run_time. To avoid synchronization of all the noise sources and improve the stochasticity of the
        search a time interval delta is defined to randomly spread the activation
        and deactivation times of the SpikeSourcePoisson objects, see diagram.

                   |comienza-|-------Full Noise--------|--termina|
                   |--delta--|                         |--delta--|
                   0%-------100%                      100%------0%
                   |------Noise=(lapso - stim_ratio)---|--stim_ratio--------|
                   |------lapso=runtime*shrink/n_popultions--------------...|
        |--phase---|------noise interval = runtime*shrink----------------------------------------...|
        |-----------------------------------------------run_time----------------------------------------------------...|

        Other stochastic search strategies to solve the CSP may be implemented modifying this method and the
        build_stimulation_pops method below.

        args:
            n_populations: number of noise populations to stimulate each CSP variable.
            full: controls if the noise deactivations should all happen after run_time or at the lapso width.
            stim_ratio: defines the portion of the stimulation window in which the noise will be active. A value
                of 0.5 will mean that stimulation happens only during the first half of the interval.
            shrink: shrinks the time interval throughout wich the noise populations will be distributed.
                It defines fraction of the run time, so it should be between 0.0 and 1.0.
            phase: a waiting time before the first noise population activates.
        returns:
            list of starting times
            lists of random distributions for start and duration of noise stimulation populations.
        """
        lapso = shrink * self.run_time / n_populations
        delta = lapso / self.run_time
        comienza = [RandomDistribution("uniform", [lapso * i+phase, lapso * i + delta+phase]) for i in range(
            n_populations)]
        if full:
            termina = [RandomDistribution("uniform", [self.run_time - delta, self.run_time]) for i
                       in range(n_populations)]
        else:
            termina = [RandomDistribution("uniform", [lapso * stim_ratio, lapso * stim_ratio + delta]) for i
                       in range(n_populations)]
        stim_times = [lapso * i+phase for i in range(n_populations)]
        return stim_times, comienza, termina

    def build_stimulation_pops(self, n_populations=1, shrink=1.0, stim_ratio=1.0, rate=(20.0, 20.0), full=True,
                               phase=0.0, clue_size=None):
        """Generate noise sources for each neuron and creates additional stimulation sources for clues.

        The noise sources are pyNN population objects of the SpikeSourcePoisson type, which generate spikes at times
        described by a Poisson random process. In pyNN these objects connect with neurons using synapses as if they
        were neurons too, these will be excitatorily  connected to the variable populations as stimulating noise.

        Currently each SpikeSourcePoisson object accepts only a start time and a duration time, thus to change the noise
        level through time one should create n different populations and activate them at different times. This method
        passes the arguments to the poisson_params method and reads activation times and duration from its returns.
        Such times define the stochastic search.

        args:
            n_populations: number of noise populations to stimulate each CSP variable.
            shrink: shrinks the time interval throughout wich the noise populations will be distributed.
                It defines fraction of the run time, so it should be between 0.0 and 1.0.
            stim_ratio: defines the portion of the stimulation window in which the noise will be active. A value
                of 0.5 will mean that stimulation happens only during the first half of the interval.
            rate: a tuple of floating-point numbers defining the rate of the Poisson process for the noise
                  populations, the first value is used for all CSP variable populations and the second value for
                  the clues.
            full: controls if the noise deactivations should all happen after run_time or at the lapso width.
            phase: a waiting time before the first noise population activates.
            clue_size: optional, number of neurons to use to stimulate clues, default value is core_size.
        """
        print msg, 'creating %d populations of SpikeSourcePoisson noise sources for each variable'%(n_populations)
        stim_times, comienza, termina = self.poisson_params(n_populations, full=full, stim_ratio=stim_ratio,
                                                            shrink=shrink, phase=phase)
        if clue_size == None:
            clue_size = self.core_size
        stim_pops = [[] for k in range(n_populations)]
        clues_stim = []
        for stimulus in range(n_populations):
            for variable in range(self.variables_number):
                stim_pops[stimulus].append(p.Population(self.size, p.SpikeSourcePoisson,
                                                 {"rate": rate[0], "start": comienza[stimulus], "duration":
                                                     termina[stimulus]}, label="stim%d_var%d" % (stimulus + 1, variable
                                                                                + 1)))
                if variable in self.clues[0]:
                    clues_stim.append(p.Population(clue_size, p.SpikeSourcePoisson,
                                                   {"rate": rate[1], "start": 0, "duration": self.run_time},
                                                   label='clues_stim%d' % variable))
        self.stim_pops = stim_pops
        self.clues_stim = clues_stim
        self.n_populations = n_populations
        self.stims = stim_times
        self.clue_size = clue_size

    def build_dissipation_pops(self, d_populations=1, shrink=1.0, stim_ratio=1.0, rate=20.0, full=True, phase=0.0):
        """Generate noise sinks for each neuron: pyNN population objects of the type SpikeSourcePoisson.

        the Poisson neural populations will be inhibitorilly connected to the variable populations, creating a
        dissipative effect. This method passes the arguments to the poisson_params method and reads activation
        times from the returns. If the clues_inhibition argument is False the clues will not receive inhibition.

        args:
            d_populations: number of dissipative noise populations to depress each CSP variable.
            shrink: shrinks the time interval throughout wich the noise populations will be distributed.
                It defines fraction of the run time, so it should be between 0.0 and 1.0.
            stim_ratio: defines the portion of the depression window in which the noise will be active. A value
                of 0.5 will mean that depression happens only during the first half of the interval.
            rate: a floating-point number defining the rate of the Poisson process for the noise
                  populations.
            full: controls if the noise deactivations should all happen after run_time or at the lapso width.
            phase: a waiting time before the first dissipation population activates.
        """
        print msg, 'creating %d populations of dissipative noise sources for each variable' % (d_populations)
        diss_times, comienza, termina = self.poisson_params(d_populations, full=full, stim_ratio=stim_ratio,
                                                            shrink=shrink, phase=phase)
        diss_pops = [[] for k in range(d_populations)]
        for k in range(d_populations):
            for variable in range(self.variables_number):
                diss_pops[k].append(p.Population(self.size, p.SpikeSourcePoisson,
                                                 {"rate": rate, "start": comienza[k], "duration": termina[k]},
                                                 label="diss%d_var%d" % (k + 1, variable + 1)))
        #TODO: if self.clues_inibition = False do not create the populations for the clues.
        self.diss_pops = diss_pops
        self.d_populations = d_populations
        self.disss = diss_times

    def connect_cores(self, w_range=[0.6, 1.2], d_range=[1.0, 1.2]):
        """Create internal excitatory connections between the neurons of each domain subpopulation of each variable.

        In the network representing the CSP, each neural population representing a variable contains a subpopulation
        for each possible value on its domain. This method connects all-to-all the neurons of each domain population
        of each variable using escitatory synapses.

        args:
            w_range: range for the random distribution of synaptic weights in the form [w_min, w_max].
            d_range: range for the random distribution of synaptic delays in the form [d_min, d_max].
        """
        print msg, 'internally connnecting the neurons of each domain of each variable'
        delays = RandomDistribution('uniform', d_range)
        weights = RandomDistribution('uniform', w_range)
        connections = [(m, n, weights.next() if m // self.core_size == n // self.core_size and m != n else 0.0,
                        delays.next()) for n in range(self.domain_size * self.core_size) for m in
                       range(self.domain_size * self.core_size)]
        for variable in range(self.variables_number):
            synapses = p.Projection(self.var_pops[variable], self.var_pops[variable], p.FromListConnector(connections,
                                                                                safe=True),
                                    target="excitatory")
            self.core_conns.append(synapses)

    def internal_inhibition(self, w_range=[-0.2, 0.0], d_range=[2.0, 2.0]):
        """Connect the domains populations of the same variable using inhibitory synapses.

        the connectiviy establishes a lateral inhibition circuit over the domains of each variable, so that most of the
        time only the neurons from a single domain population are active.

        args:
            w_range: range for the random distribution of synaptic weights in the form [w_min, w_max].
            d_range: range for the random distribution of synaptic delays in the form [d_min, d_max].
        """
        print msg, 'Creating lateral inhibition between domains of each variable'
        delays = RandomDistribution('uniform', d_range)
        weights = RandomDistribution('uniform',w_range)
        connections = [(m, n, 0.0 if m // self.core_size == n // self.core_size else  weights.next(), delays.next()) for
                       n in range(self.size) for m in range(self.size)]
        for variable in range(self.variables_number):
            if self.clues_inhibition:
                synapses = p.Projection(self.var_pops[variable], self.var_pops[variable], p.FromListConnector(
                    connections, safe=True),
                                        target="inhibitory")
                self.internal_conns.append(synapses)
            elif variable not in self.clues:
                synapses = p.Projection(self.var_pops[variable], self.var_pops[variable], p.FromListConnector(
                    connections, safe=True),
                                        target="inhibitory")
                self.internal_conns.append(synapses)

    def stimulate_cores(self, w_range=[1.4, 1.4], d_range=[1.0, 1.0], w_clues=[1.4, 1.6]):  # w_clues=[0.0, 0.2]
        """Connect stimulating noise sources to variables populations.

        args:
            w_range: range for the random distribution of synaptic weights in the form [w_min, w_max].
            d_range: range for the random distribution of synaptic delays in the form [d_min, d_max].
            w_clues: clues specific range for the random distribution of synaptic weights in the form [w_min, w_max].
        """
        p.set_number_of_neurons_per_core(p.IF_curr_exp, 150)
        print msg, 'connecting Poisson noise sources to neural populations for stimulation'
        delays = RandomDistribution('uniform', d_range)
        weights = RandomDistribution('uniform', w_range)
        weight_clues = RandomDistribution("uniform", w_clues)
        for stimulus in range(self.n_populations):
            for variable in range(self.variables_number):
                counter = 0
                if variable in self.clues[0]:
                    shift = self.clues[1][self.clues[0].index(variable)] * self.core_size
                    connections = [(m, n + shift, weight_clues.next(), delays.next()) for m in range(self.core_size) for
                                   n in range(self.clue_size)]
                    synapses = p.Projection(self.clues_stim[counter], self.var_pops[variable],
                                            p.FromListConnector(connections, safe=True), target='excitatory')
                    counter += 1
                    self.stim_conns.append(synapses)
                else:
                    synapses = p.Projection(self.stim_pops[stimulus][variable], self.var_pops[variable],
                                            p.OneToOneConnector(weights=weights, delays=delays), target='excitatory')
                    self.stim_conns.append(synapses)
        self.stim_times += self.stims

    def depress_cores(self, w_range=[-2.0, -1.5], d_range=[2.0, 2.0]):
        """Connect depressing noise sources to variables populations.

        args:
            w_range: range for the random distribution of synaptic weights in the form [w_min, w_max].
            d_range: range for the random distribution of synaptic delays in the form [d_min, d_max].
        """
        print msg, 'connecting Poisson noise sources to neural populations for dissipation'
        delays = RandomDistribution('uniform', d_range)
        weights = RandomDistribution('uniform', w_range)
        for depressor in range(self.d_populations):
            for variable in range(self.variables_number):
                if variable not in self.clues[0]:
                    synapses = p.Projection(self.diss_pops[depressor][variable], self.var_pops[variable],
                                            p.OneToOneConnector(weights=weights, delays=delays), target='inhibitory')
                    self.diss_conns.append(synapses)
        self.diss_times += self.disss

    def apply_constraints(self, kind="inhibitory", w_range=[-0.2, -0.0], d_range=[2.0, 2.0], random_cons=False,
                          pAF=0.5):
        """Map constraints list to inhibitory or excitatory connections between neural populations.

        The clues_inhibition class variable determines whether clues should receive inhibitory connections or not.

        args:
            kind: whether constraints are inhibitory or excitatory.
            w_range: range for the random distribution of synaptic weights in the form [w_min, w_max].
            d_range: range for the random distribution of synaptic delays in the form [d_min, d_max].
            random_cons: whether constraints are randomly choosen to be inhibitory or excitatory with probability pAF.
            pAF: probability of inhibitory connections, as a probability it should be between 0.0 and 1.0. It only works
                when random_cons is True.
         """
        delays = RandomDistribution('uniform', d_range)
        weights = RandomDistribution('uniform', w_range)  # 1.8 2.0 spin_system
        if 'weight' in self.constraints[0]:
            print msg, '''creating constraints between CSP variables with specified weights and randomly distributed
                        delays'''
        else:
            print msg, '''creating constraints between CSP variables with random and  uniformelly distributed delays
                       and weights'''
        for constraint in self.constraints:
            source = constraint['source']
            target = constraint['target']
            if random_cons:
                kind = np.random.choice(['inhibitory', 'excitatory'], p=[pAF, 1-pAF])
            #TODO find a way of reducing the next two conditionals, they're equal except for conditioning on target...
            #TODO ... being a clue.
            if self.clues_inhibition:
                connections = []
                for n in range(self.size):
                    for m in range(self.size):
                        if 'weight' in constraint:
                            weight = constraint['weight']
                        else:
                            weight = weights.next()
                        connections.append(
                            (m, n, weight if m // self.core_size == n // self.core_size else 0.0, delays.next()))
                synapses = p.Projection(self.var_pops[source], self.var_pops[target],
                                        p.FromListConnector(connections, safe=True), target=kind)
                self.constraint_conns.append(synapses)
                if self.directed == False:
                    synapses = p.Projection(self.var_pops[target], self.var_pops[source],
                                            p.FromListConnector(connections, safe=True), target=kind)
                    self.constraint_conns.append(synapses)
            elif target not in self.clues[0]:
                connections = []
                for n in range(self.size):
                    for m in range(self.size):
                        if 'weight' in constraint:
                            weight = constraint['weight']
                        else:
                            weight = weights.next()
                        connections.append(
                            (m, n, weight if m // self.core_size == n // self.core_size else 0.0, delays.next()))
                synapses = p.Projection(self.var_pops[source], self.var_pops[target],
                                        p.FromListConnector(connections, safe=True), target=kind)
                self.constraint_conns.append(synapses)
                if self.directed == False:
                    synapses = p.Projection(self.var_pops[target], self.var_pops[source],
                                            p.FromListConnector(connections, safe=True), target=kind)
                    self.constraint_conns.append(synapses)

    def initialize(self, v_range=[-65.0, -55.0]):
        """Randomly initialize the membrane voltage of the neurons in the range v_range.

        args:
            v_range: range for the random distribution of membrane potentials in the form [v_min, v_max].
        """
        print msg, 'randomly setting the initial voltage for each variable population'
        for variable in self.var_pops:
            variable.initialize("v", RandomDistribution("uniform", v_range))

    def recording(self):
        """Record spikes from neural populations representing CSP variables.

        If live class variable is set to True this method also activate live output for the neural populations
        representing CSP variables.
        """
        print msg, 'activating recording for variable populations'
        for population in self.var_pops:
            population.record()
        # Activate live output to be used for example with the Sudoku visualiser.
        if self.live:
            print 'activating live output'
            p.external_devices.activate_live_output_for(self.var_pops)

    def record_stimulation(self):
        """Record spikes from stimulating noise sources."""
        print msg, 'activating recording for noise populations'
        for stimulus in self.stim_pops:
            for population in stimulus:
                population.record()

    def record_dissipation(self):
        """Record spikes from depressing noise sources."""
        print msg, 'activating recording for dissipation populations'
        for depressor in self.diss_pops:
            depressor.record()

    def save(self, filename, DAT=False):
        """Save spikes recorded from neural populations representing CSP variables.

        The recording() method should be called first in order to make spikes available for saving. All files will be
        saved into the results folder.

        args:
            filename: prefix of the file name where spikes will be saved in binary format. The full name will be:
                filename_spikes_binary.
            DAT: whether spikes should be saved also in .dat format on an additional file.
        """
        if not os.path.exists('results'):
            os.makedirs('results')
        if DAT:
            print msg, 'saving spikes from CSP variables to file results/%s_variable#.dat' % filename
            for var_index, population in enumerate(self.var_pops):
                population.printSpikes('results/%s_variable%d.dat' % (filename, var_index))
        with open('results/%s_spikes_binary' % filename, 'w+') as file:
            # TODO refactor params as dictionary with *kwargs
            params = [self.run_time,
                      self.variables_number,
                      self.size,
                      self.domain_size,
                      self.core_size,
                      self.constraints,
                      self.stim_times,
                      self.diss_times]
            np.save(file, params)
            for population in self.var_pops:
                spikes = population.getSpikes()
                np.save(file, spikes)
        self.spikes_file = filename

    def save_stimulation(self, filename, DAT=False):
        """Save spikes recorded from stimulating noise sources.

        The record_stimulation() method should be called first in order to make spikes available for saving.

        args:
            filename: prefix of the file name where spikes will be saved in .dat format. The full name will be:
                stim_#_filename_variables#.dat.
            DAT: whether spikes should be saved also in .dat format on an additional file.
        """
        if DAT:
            for pulse_index, pulse in enumerate(self.stim_pops):
                print msg, 'saving spikes from noise sources to file results/stim_%d_%s.dat' % (pulse_index, filename)
                for var_index, population in enumerate(pulse):
                    population.printSpikes('results/stim_%d_%s_variables%d.dat' % (pulse_index, filename, var_index))
        with open('results/%s_stim_spikes_binary' % filename, 'w+') as file:
            for pulse in self.stim_pops:
                for population in pulse:
                    spikes = population.getSpikes()
                    np.save(file, spikes)

    def save_dissipation(self, filename, DAT=False):
        """Save spikes recorded from depressing noise sources.

        The record_dissipation() method should be called first in order to make spikes available for saving.

        args:
            filename: prefix of the file name where spikes will be saved in .dat format. The full name will be:
                diss_#_filename_variables#.dat.
            DAT: whether spikes should be saved also in .dat format on an additional file.
        """
        if DAT:
            for pulse_index, pulse in enumerate(self.diss_pops):
                print msg, 'saving dissipation to file results/diss_%d_%s.dat' % (pulse_index, filename)
                for var_index, population in enumerate(pulse):
                    population.printSpikes('results/diss_%d_%s_variable%d.dat' % (pulse_index, filename, var_index))
        with open('results/%s_diss_spikes_binary' % filename, 'w+') as file:
            for pulse in self.diss_pops:
                for population in pulse:
                    spikes = population.getSpikes()
                    np.save(file, spikes)

    def report_network_params(self, filename=False):
        """ Report the network dimensions and parameters to standard output or file.

        args:
            filename: name of file in which report will be saved. if not specified the report will show up only in the
                standard output.
        """
        # Count populations.
        if not os.path.exists('results'):
            os.makedirs('results')
        var_pops_num = len(self.var_pops)
        diss_pops_num = len(self.diss_pops)
        stim_pops_num = len(self.stim_pops) * len(self.stim_pops[0])
        pops_number = var_pops_num + diss_pops_num + stim_pops_num
        # Count neurons.
        var_neurons = var_pops_num * self.size
        stim_neurons = stim_pops_num * self.size
        diss_neurons = diss_pops_num * self.size
        net_neurons = var_neurons + stim_neurons + diss_neurons

        # Create function to count synapses from projections list.
        def projections_counter(projections_pop):
            """ Count synapses from a list of projections.

            args:
                projections_pop: list of projections, collected when calling other methods.

            returns:
                synapse_number: the number of synapses generated by the list of projections.
            """
            # synapse_number = sum(map(lambda x: len(x.getWeights('list').flatten()), projections_pop))
            synapse_number = lambda x: x.getWeights('list').flatten()
            synapse_number = map(synapse_number, projections_pop)
            synapse_number = map(len, synapse_number)
            synapse_number = sum(synapse_number)
            return synapse_number

        # Count synapses created by methods that connected neural populations.
        core_conns = projections_counter(self.core_conns)
        internal_conns = projections_counter(self.internal_conns)
        stim_conns = projections_counter(self.stim_conns)
        diss_conns = projections_counter(self.diss_conns)
        constraint_conns = projections_counter(self.constraint_conns)
        net_conns = sum([core_conns, internal_conns, stim_conns, diss_conns, constraint_conns])

        # Report template.
        report = """
        |======== Network Parameters =========|
        |_____________________________________|
        |Total number of neurons:     %d
        |      variables neurons:     %d
        |    stimulation neurons:     %d
        |    dissipation neurons:     %d
        |_____________________________________
        |Total number of populations: %d
        |      variables populations: %d
        |    dissipation populations: %d
        |   stimulation poppulations: %d
        |                domain size: %d
        |                  core size: %d
        |_____________________________________
        |Total number of synapses:    %d
        |    stimulating synapses:    %d
        |    dissipating synapses:    %d
        |    constraints synapses:    %d
        |   var internal synapses:    %d
        |  core internal synapses:    %d
        |=====================================|
        """ % (net_neurons, var_neurons, stim_neurons, diss_neurons, pops_number, var_pops_num,
               diss_pops_num, stim_pops_num, self.domain_size, self.core_size, net_conns,
               stim_conns, diss_conns, constraint_conns, internal_conns, core_conns)
        print report
        # Print report to file.
        if filename:
            with open('results/%s.dat' % filename, 'w+') as file:
                file.write(report)

max_stars_repo_name: GAFonsecaGuerra/SpiNNakerCSPs

max_issues_repo_name: GAFonsecaGuerra/SpiNNakerCSPs

max_forks_repo_name: GAFonsecaGuerra/SpiNNakerCSPs


text:
"""
Utility functions and classes
"""

from functools import wraps
import traceback
import numpy as np

class once(dict):
    """identical to dict; used to yield something only once"""
    pass

class Bunch:
    """convert a dictionary into a class with data members equal to the dictionary keys"""
    def __init__(self, adict):
        self.__dict__.update(adict)

    def __getitem__(self, key):
        return self.__dict__[key]

    def __setitem__(self, key, value):
        self.__dict__[key] = value

def doublewrap(func):
    """
    a decorator decorator, can be used as @decorator(...) or @decorator
    """

    @wraps(func)
    def new_func(*args, **kwargs):
        if len(args) == 2 and len(kwargs) == 0 and callable(args[1]):
            return func(*args)
        else:
            return lambda f: func(args[0], f, *args[1:], **kwargs)

    return new_func

def doublewrap_1(func):
    """
    a decorator decorator, can be used as @decorator(...) or @decorator
    """

    @wraps(func)
    def new_func(*args, **kwargs):
        if len(args) == 1 and len(kwargs) == 0 and callable(args[1]):
            return func(*args)
        else:
            return lambda f: func(args[0], f, *args[1:], **kwargs)

    return new_func

def yield_traceback(func):
    """decorator to properly yield the traceback of a function in a parallel environment"""

    def new_func(*args, **kwargs):
        try:
            return f(*args, **kwargs)
        except:
            raise Exception("".join(traceback.format_exception(*sys.exc_info())))

    return new_func

class first_argument:
    # current_iteration = {}
    def __init__(self, name, num_iterations=None):
        self.name = name
        self.num_iterations = num_iterations
        if num_iterations is not None:
            current_iteration[name] = Value('i', 0)

    def iterations(self):
        def gen():
            for i in range(self.num_iterations):
                yield current_iteration[self.name].value
                current_iteration[self.name].value += 1
        return gen()

def flatten_along(arr, axis=None):
    """flatten an array along a give axis / axes"""
    if axis is None:
        return arr.flatten()
    else:
        if np.isscalar(axis):
            arr = np.moveaxis(arr, axis, 0)
        else:
            axis = sorted(axis)
            for i,ax in enumerate(axis):
                arr = np.moveaxis(arr, ax, i)
            shape = arr.shape[len(axis):]
            arr = arr.reshape((-1,) + shape)
        return arr

ext: py

max_stars_repo_name: johnaparker/numpipe

max_issues_repo_name: johnaparker/numpipe

max_forks_repo_name: johnaparker/numpipe

reason: import numpy


text:
import numpy as np
from cosymlib.file_io import read_generic_structure_file, get_file_xyz_txt
from cosymlib.molecule.geometry import Geometry

def read_normal_modes_gaussian(file_name):
    read = False
    frequencies = []
    normal_modes_matrices = []
    with open(file_name, mode='r') as lines:
        for line in lines:
            if read:
                if len(line.split()) == 3 or not line.strip():
                    read = False
                else:
                    for i in range(int(len(line.split()[2:])/3)):
                        atom_vibration = line.split()[2 + 3*i:3 + 2 + 3*i]
                        normal_modes_matrices[-3+i].append([float(x)for x in atom_vibration])

            if 'Atom  AN' in line:
                read = True
                for _ in range(line.split()[2:].count('X')):
                    normal_modes_matrices.append([])
            if 'Frequencies --' in line:
                frequencies.append(line.split()[2:])
    frequencies = [float(frequency) for sublist in frequencies for frequency in sublist]

    return frequencies, normal_modes_matrices

def get_nm_vibration_path(geometry, normal_mode, points=10, backward=False, k=1.5):

    nm_np_matrix = np.array(normal_mode)
    geom_np = np.array(geometry)
    if backward:
        x = np.linspace(-k, k, points)
    else:
        x = np.linspace(k/points, k, points)

    structures_path = []
    for dx in x:
        structures_path.append(geom_np + dx*nm_np_matrix)

    return structures_path

molecule = read_generic_structure_file('sf6.fchk')
freq, nm_martices = read_normal_modes_gaussian('../old_examples/sf6_freq.out')
k_points = 1.5
points = 10
n_freq = 1
path = get_nm_vibration_path(molecule.geometry.get_positions(), nm_martices[n_freq - 1], k=k_points)

x = np.linspace(k_points/points, k_points, points)
output = open('sf6_freq{}.xyz'.format(n_freq), 'w')
for ids, structure in enumerate(path):
    xi = '{:.2f}'.format(x[ids]).replace('.', '')
    txt = get_file_xyz_txt(Geometry(structure, symbols=molecule.geometry.get_symbols(), name='freq : ' + str(freq[0]) +
                                                                                ' ' + str(xi)))
    output.write(txt)


text:
import numpy as np
import torch
from torch.autograd import Variable
import sys
import os
sys.path.append(os.path.join(os.getcwd(), 'model/cffi'))
import matplotlib as mpl
mpl.use('Agg')
from utils.visualize import save_occupancy_fig, save_mesh_fig
from utils.config import parse_args
from model.table import get_accept_topology
from model.loss import Loss
from data.data_loader import load_data, get_batch

if torch.cuda.is_available():
    dtype = torch.cuda.FloatTensor
    dtype_long = torch.cuda.LongTensor
else:
    dtype = torch.FloatTensor
    dtype_long = torch.LongTensor

def run_val(model, loss_obj, data_val, args, phase='train'):
    """ Test with the trained model

    Input:
        model: the Deep Marching Cubes model
        loss_obj: the loss instance
        data_val: validation data blob, including points and voxel grid
        args: configuration arguments
        phase: 'train' or 'val'
    """

    max_prob = []
    # evaluation
    loss_eval = 0
    for itest in range(args.num_val):
        sys.stdout.write('.')
        sys.stdout.flush()

        itest_ = Variable(dtype_long([itest]))
        net_input, pts_rnd, _ = get_batch(data_val, itest_, args)

        offset, topology, occupancy = model(net_input)

        loss = loss_obj.loss_eval(offset[0], topology[0], pts_rnd[0])

        loss_eval += loss.data.cpu()[0]

        interval = args.verbose_interval
        if phase == 'val':
            interval = 1
        if np.mod(itest, interval) == 0:
            topology_fused = topology[-1].data.cpu().numpy()
            topology_fused = np.maximum(topology_fused[:, 0:128],
                                        topology_fused[:, 256:127:-1])
            topology_fused = topology_fused[:, get_accept_topology()]
            save_occupancy_fig(
                    pts_rnd[-1].data.cpu().numpy(),
                    occupancy[-1].data.cpu().numpy(),
                    loss_obj.x_grids,
                    loss_obj.y_grids,
                    loss_obj.z_grids,
                    itest, args, 'val')

            topology_vis = topology[:, :, loss_obj.visTopology]

            save_mesh_fig(
                    pts_rnd[-1].data.cpu().numpy(),
                    offset[-1],
                    topology_vis[-1],
                    loss_obj.x_grids,
                    loss_obj.y_grids,
                    loss_obj.z_grids,
                    itest, args, 'val')

    print('')
    return loss_eval

if __name__ == '__main__':

    # parse args
    args = parse_args()

    # load data
    args, data_val = load_data(args, dtype, 'val')

    # setup loss object
    loss_obj = Loss(args)

    # initialize the model
    assert(os.path.isfile(args.model))
    print "Validating with snapshotted model %s ..." % args.model
    deep_marching_cubes = torch.load(args.model)
    if torch.cuda.is_available():
        deep_marching_cubes.cuda()

    # validation
    loss = run_val(deep_marching_cubes, loss_obj, data_val, args, 'val')
    print('============== average loss:%f' % (loss/args.num_val))

    print 'Done!'


text:
import torch
import torch.nn as nn
from base_learner import BaseLearner
from model.pytorch_lstm import *
from torch.utils.data import DataLoader
from smoke_video_dataset import SmokeVideoDataset
import torch.optim as optim
from torch.optim import lr_scheduler
import torch.nn.functional as F
from torch.autograd import Variable
import uuid
from sklearn.metrics import classification_report
from sklearn.metrics import accuracy_score
import numpy as np
import random
from torchvision.transforms import Compose
from video_transforms import RandomResizedCrop, RandomHorizontalFlip, ColorJitter, RandomPerspective, RandomErasing, Resize, Normalize
from torch.utils.tensorboard import SummaryWriter
import json

# Conv LSTM learner
class LSTMLearner(BaseLearner):
    def __init__(self,
                 batch_size=8, # Batch 8 per GPU
                 lr=0.01,
                 max_steps=64e3,
                 momentum=0.9,
                 #milestones = [40, 500, 1000],
                 gamma = 0.1,
                 weight_decay = 0.000001,
                 num_workers=1,
                 num_of_action_classes=2,
                 num_steps_per_update=2,
                 num_steps_per_check=10,
                 #use_cuda=torch.cuda.is_available(),
                 use_cuda=True,
                 parallel=True,
                 augment=True,
                 save_model_path="../data/saved_lstm/",  # path for saving the models
                 train_writer_p="../data/ts_runs/lstm_train",
                 val_writer_p="../data/ts_runs/lstm_val"
                 ):
        super().__init__()
        self.create_logger(log_path="LSTMLearner_test.log")
        self.log("Use LSTM learner")

        self.batch_size = batch_size
        self.lr = lr
        self.max_steps = max_steps
        self.momentum = momentum
        #self.milestones = milestones
        self.gamma = gamma
        self.weight_decay = weight_decay
        self.num_workers = num_workers
        self.num_of_action_classes = num_of_action_classes
        self.num_steps_per_update = num_steps_per_update
        self.num_steps_per_check = num_steps_per_check
        self.use_cuda = use_cuda
        self.parallel = parallel
        self.augment = augment
        self.save_model_path = save_model_path
        self.train_writer_p = train_writer_p
        self.val_writer_p = val_writer_p

    def random_frames_from_batch(self, data):
        b, c, f, h, w = list(data.shape)
        new_batch = torch.zeros(b, c, h, w)
        for i in range(b):
            frame = random.randint(0, f - 1)
            new_batch[i, :, :, :] = data[i, :, frame, :, :]
        return new_batch

    def compress_videos_to_frames(self, n):
        b, c, t, h, w = n.shape
        ret = torch.zeros(b, c*t, h, w)
        frame = 0
        for channel in range(c*t):
            ret[:,channel,:,:] = n[:,channel%c, frame,:,:]
            if channel > 0 and channel % c == 0: frame += 1
        return ret

    def set_dataloader(self, metadata_path, p_vid, mode, tf):
        dataloader = {}
        for phase in metadata_path:
            self.log("Create dataloader for " + phase)
            dataset = SmokeVideoDataset(metadata_path=metadata_path[phase], root_dir=p_vid, transform=tf[phase]) #mode=mode, transform=tf[phase])
            dataloader[phase] = DataLoader(dataset, batch_size=self.batch_size, shuffle=True, num_workers=self.num_workers, pin_memory=True)
        return dataloader

    def labels_to_list(self, labels):
        # Convert labels obtained from the dataloader to a list of action classes
        return np.argmax(labels.numpy().max(axis=2), axis=1).tolist()

    def to_variable(self, v):
        if self.use_cuda:
            v = v.cuda()  # move to gpu
        return Variable(v)

    def make_pred(self, model, frames):
        # Upsample prediction to frame length (because we want prediction for each frame)
        return F.interpolate(model(frames), frames.size(2), mode="linear", align_corners=True)

    def fit(self,
            mode="rgb",
            p_metadata_train="../data/split/metadata_train_split_0_by_camera.json",
            p_metadata_validation="../data/split/metadata_validation_split_0_by_camera.json",
            p_vid="../data/"):

        self.log("="*60)
        self.log("="*60)
        self.log("Start training...")

        if mode == "rgb":
            model = CNN_LSTM()
            p_vid = p_vid + "rgb/"

        elif mode == "flow":
            model = CNN_LSTM()
            p_vid = p_vid + "flow/"

        #self.load(model, "../data/saved_lstm/b1ca514-lstm-rgb/12912.pt")

        if self.use_cuda:
            model.cuda()
            if self.parallel and torch.cuda.device_count() > 1:
                self.log("Let's use " + str(torch.cuda.device_count()) + " GPUs!")
                model = nn.DataParallel(model)

        writer_train = SummaryWriter(self.train_writer_p)
        writer_val = SummaryWriter(self.val_writer_p)

        # Load datasets
        metadata_path = {"train": p_metadata_train, "validation": p_metadata_validation}
        ts = self.get_transform()
        transform = {"train": ts, "validation": ts}
        if self.augment:
            transform["train"] = self.get_transform(phase="train")
        dataloader = self.set_dataloader(metadata_path, p_vid, mode, transform)

        # Set optimizer
        optimizer = optim.SGD(params=model.parameters(), lr=self.lr, momentum=self.momentum, weight_decay = self.weight_decay)
        # lr_sche = optim.lr_scheduler.MultiStepLR(optimizer, milestones=self.milestones, gamma=self.gamma)

        # Set Loss Fum nction
        criterion = nn.BCEWithLogitsLoss()

        # Set logging format
        log_fm = "%s step: %d lr: %r loss: %.4f"

        # Train and validate
        steps = 0
        epochs = 0
        nspu = self.num_steps_per_update
        nspc = self.num_steps_per_check
        nspu_nspc = nspu * nspc
        model_id = str(uuid.uuid4())[0:7] + "-lstm-" + mode
        accum = {} # counter for accumulating gradients
        tot_loss = {} # total loss
        tot_loc_loss = {} # total localization loss
        tot_cls_loss = {} # total classification loss
        pred_labels = {} # predicted labels
        true_labels = {} # true labels
        for phase in ["train", "validation"]:
            accum[phase] = 0
            tot_loss[phase] = 0.0
            tot_loc_loss[phase] = 0.0
            tot_cls_loss[phase] = 0.0
            pred_labels[phase] = []
            true_labels[phase] = []
        while steps < self.max_steps:
            self.log("-"*40)
            # Each epoch has a training and validation phase
            for phase in ["train", "validation"]:
                self.log("phase " + phase)
                if phase == "train":
                    epochs += 1
                    self.log("epochs: %d steps: %d/%d" % (epochs, steps, self.max_steps))
                    model.train(True) # set model to training mode
                else:
                    model.train(False) # set model to evaluate mode
                optimizer.zero_grad()
                # Iterate over data
                for d in dataloader[phase]:
                    accum[phase] += 1
                    # Get inputs
                    frames, labels = d["frames"], d["labels"]
                    #frames = self.compress_videos_to_frames(frames)
                    frames = self.to_variable(frames)
                    true_labels[phase] += self.labels_to_list(labels)
                    labels = self.to_variable(d["labels"])
                    pred = model(frames)
                    pred = pred.permute(0,2,1)
                    _, predicted = torch.max(pred, dim=1)
                    cpu_pred = predicted.cpu().detach().tolist()
                    for lst in cpu_pred:
                        pred_labels[phase].append(lst[-1])
                    # Compute classification loss (with max-pooling along time, batch x channel x time)
                    cls_loss = criterion(torch.max(pred, dim=2)[0], torch.max(labels, dim=2)[0])
                    tot_cls_loss[phase] += cls_loss.data
                    # Backprop
                    loss = cls_loss / nspu
                    tot_loss[phase] += loss.data
                    if phase == "train":
                        loss.backward()
                    # Accumulate gradients during training
                    if (accum[phase] == nspu) and phase == "train":
                        steps += 1
                        accum[phase] = 0
                        optimizer.step()
                        optimizer.zero_grad()
                        # lr_sche.step()
                        if steps % nspc == 0:
                            #lr = lr_sche.get_lr()[0]
                            lr = self.lr
                            tl = tot_loss[phase]/nspc
                            writer_train.add_scalar("Training Loss", tl, global_step=steps)
                            self.log(log_fm % (phase, steps, lr, tl))
                            tot_loss[phase] = 0.0
                if phase == "validation":
                    # lr = lr_sche.get_lr()[0]
                    lr = self.lr
                    tl = (tot_loss[phase]*nspu)/accum[phase]
                    writer_val.add_scalar("Validation Loss", tl, global_step=steps)
                    self.log(log_fm % (phase, steps, lr, tl))
                    tot_loss[phase] = 0.0
                    accum[phase] = 0
                    p_model = self.save_model_path + model_id + "/" + str(steps) + ".pt"
                    self.save(model, p_model)
                    for phase in ["train", "validation"]:
                        self.log("Performance for " + phase)
                        accuracy = accuracy_score(true_labels[phase], pred_labels[phase])
                        if phase == "train": writer_train.add_scalar("Train Accuracy", accuracy, global_step=steps)
                        elif phase == "validation": writer_val.add_scalar("Val Accuracy", accuracy, global_step=steps)
                        self.log(classification_report(true_labels[phase], pred_labels[phase]))
                        pred_labels[phase] = []
                        true_labels[phase] = []

        self.log("Done fit")

    def test(self,
                mode="rgb",
                p_metadata_test="../data/split/metadata_test_split_0_by_camera.json",
                p_vid="../data/",
                p_model=None):
        self.log("="*60)
        self.log("="*60)
        self.log("Start testing...")

        if mode == "rgb":
            model  = CNN_LSTM()
            p_vid = p_vid + "rgb/"

        elif mode == "flow":
            model  = CNN_LSTM()
            p_vid = p_vid + "flow/"

        self.load(model, p_model)

        if self.use_cuda:
            model.cuda()
            if self.parallel and torch.cuda.device_count() > 1:
                self.log("Let's use " + str(torch.cuda.device_count()) + " GPUs!")
                model = nn.DataParallel(model)

        # Load dataset
        metadata_path = {"test": p_metadata_test}
        transform = {"test": self.get_transform()}
        dataloader = self.set_dataloader(metadata_path, p_vid, mode, transform)

        #self.load(model, p_model)

        all_pred = {}
        # Test
        model.train(False)
        true_labels = []
        pred_labels = []
        counter = 0
        with torch.no_grad():
            for d in dataloader["test"]:
                if counter % 20 == 0:
                    self.log("Process batch " + str(counter))
                counter += 1
                # Get inputs
                frames, labels = d["frames"], d["labels"]
                names = d["file_name"]
                #frames = self.compress_videos_to_frames(frames)
                frames = self.to_variable(frames)
                l = self.labels_to_list(labels)
                true_labels += l
                labels = self.to_variable(d["labels"])
                pred = model(frames)
                pred = pred.permute(0,2,1)
                all_pred[counter] = {}
                p = pred.cpu().detach()
                _, predicted = torch.max(pred, dim=1)
                cpu_pred = predicted.cpu().detach().tolist()
                for i, lst in enumerate(cpu_pred):
                    all_pred[counter][i] = []
                    all_pred[counter][i].append("neg")
                    all_pred[counter][i].append(p[i,0,-1].item())
                    all_pred[counter][i].append("pos")
                    all_pred[counter][i].append(p[i,1,-1].item())
                    all_pred[counter][i].append("actual")
                    all_pred[counter][i].append(l[i])
                    all_pred[counter][i].append("name")
                    all_pred[counter][i].append(names[i])
                    pred_labels.append(lst[-1])
                #frames = self.to_variable(frames)
                #true_labels += self.labels_to_list(labels)
                #labels = self.to_variable(d["labels"])
                #pred = model(frames)
                #pred.permute(0,2,1)
                #_, predicted = torch.max(pred, dim=1)
                #cpu_pred = predicted.cpu().detach().tolist()
                #for lst in cpu_pred:
                #    pred_labels.append(lst[-1])
        self.log(classification_report(true_labels, pred_labels))

        with open ('lstm_pred.json', 'w') as f:
            json.dump(all_pred, f)

        self.log("Done test")

max_stars_repo_name: yenchiah/deep-smoke-machine

max_issues_repo_name: yenchiah/deep-smoke-machine

max_forks_repo_name: yenchiah/deep-smoke-machine


text:
"""Main functionality to calculate array-I-V-curves.

Every relevant paramter is contained inside the solar_array class to run the
actual calculations.

TODO:
    Expose Illumination to photocurrent function
"""

from __future__ import annotations
import io
import os
import requests
import numpy as np
import datetime as dt
import pandas as pd
import matplotlib.image as mpimg
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
from pathlib import Path
from typing import List, Dict, TYPE_CHECKING

from scipy.ndimage import zoom
from ..pvgis import Pvgis
from .solar import SunData
from .pv_panel import Pv_panel
# from .solar import *
if TYPE_CHECKING:
    from ..construction import Building

class Shading():
    """Class to contain manage rendered images.

    Attributes:
        images (Dict[str, np.ndarray]): Dictionary of rendered images.
        shape (tuple(int,int)): Shape of rendered images in px.
        load_images (type):
        dim_x (float): Size of image in X direction in m
        dim_y (float): Size of image in Y direction in m
        avg (np.ndarray): Non weighted average of all rendered images
            (For display purposes).
        building (Building): Root building for the Shading
        zero (np.ndarray): Zero direct irradiaton 'image'
        render_dir (str): Path to rendered images
        suns (SunData): Sun object for current location . (Includes horizon)
        scale (float): Image scale in px/m. Defaults to 20 px/m

    """

    def __init__(self, render_dir: str, suns: SunData, scale: float = 0):
        """Initialize Shading with the following arguments.

        Args:
            render_dir (str): Path to rendered images
            suns (SunData): Sun object for current location
            scale (float): Image scale in px/m. Defaults to 20 px/m

        Returns:
            Shading: Object instance

        """
        self.render_dir = Path(render_dir)
        self.images, self.shape = self.load_images(render_dir)
        self.suns = suns
        if scale is 0:
            self.scale = 20.0
        else:
            self.scale = scale
        self.dim_x = self.shape[1]/self.scale
        self.dim_y = self.shape[0]/self.scale
        self.avg = self.average_img()
        self.building = None
        self.zero = np.zeros(self.shape)

    def load_images(self, render_dir: str):
        """Load rendered images and append to dictionary.

        Args:
            render_dir (str): Directory of rendered images

        Returns:
            Dict[str, np.ndarray]: Dictionary of rendered images.
            shape (tuple(int,int)): Shape of rendered images in px.

        """
        render_path = Path(render_dir)
        image_paths = os.listdir(render_path)
        images: Dict['str', 'np.ndarray'] = {}
        for file_path in image_paths:
            if file_path[-3:] == 'png':
                images.update(self.fname_to_dict(
                    render_path.joinpath(Path(file_path))))
            else:
                print('Skipped {}. Not a png file'.format(file_path))
        return images, images[next(iter(images))].shape

    def fname_to_dict(self, file_path: Path, threshhold: float = .01):
        """Load an image from file_path and create a dict based on the
           filename.

        Args:
            file_path (Path): Path to image.
            threshhold (float): Binarisation threshhold

        Returns:
            Dict[str, np.ndarray]: .

        Todo:
            newer versions of matplotlib are compatbile with pathlib, making
            the string conversion useless

        """
        img = (mpimg.imread(str(file_path)) > threshhold).astype(np.bool)
        return({file_path.stem: img})

    def average_img(self):
        """Calculate the non weighted average of all rendered images.

        Returns:
            np.ndarray: Non weighted average

        """
        img_sum = np.zeros(self.shape)
        for ind, (_, image) in enumerate(self.images.items()):
            img_sum += image
        return img_sum/ind

    def __getitem__(self, index):
        """Return shading image corresponding to index.

        Args:
            index (List[int,int], dt.datetime): index to get the shading image.
                Can be a list containing an integer of the sun azimuth and
                elevation or a datetime. Using datetime requires a building.

        Returns:
            image np.ndarray: requested shading image as numpy array.

        """
        if isinstance(index, list):
            az, al = index
            return self.images['{}_{}'.format(az, al)]
        elif isinstance(index, dt.datetime):
            "Need to register a building to use datetime"
            assert not(
                self.building is None)
            sun = self.building.sundata[index]
            if sun['is_up']:
                az, al = int(sun['azimuth']), int(sun['elevation'])
                return self.images['{}_{}'.format(az, al)]
            else:
                return self.zero
        else:
            raise ValueError(
                'Index needs to be a list of azimuth and elevation' +
                ', or datetime')

class solar_array():
    """Main class for shaded solar array simulation.

    Attributes:
        shading (type): Shading object for direct irradiaton shadows.
        offset_x (float): Offset of the origin of the array along the x-Axis
            from the top left corner of the image. Defaults to 0.
        offset_y (float): Offset of the origin of the array along the y-Axis
            from the top left corner of the image. Defaults to 0.
        grid (type):
        timestep (int): Simulation timestep in seconds.
        building (Building): Building the array is attached to.
        pv_panel (Pv_panel): Pv_panel object, containing information about the
            used panels in the array
        azimuth (float): azimuth orientation of the array
        inclination (float): inclination of the array
        resolution_I (int): Current (A) resolution for the simulation:
            Defaults to 110.
        start_year (int): Initial year for the simulation. Defaults to 2007.
        end_year (int): Last year for the simulation. Defaults to 2008.
        output_dir (str): Path to save i-v-curves.

    """

    def __init__(self, shading: Shading, building: Building,
                 pv_panel: Pv_panel, render_dir: str, azimuth: float,
                 inclination: float, resolution_I: int = 110,
                 start_year: int = 2007, end_year: int = 2008):
        """Initialize the solar_array.

        Args:
            shading (Shading): Shading object for direct irradiaton shadows.
            building (Building): Building the array is attached to
            pv_panel (Pv_panel): Pv_panel object, containing information about
                the used panels in the array
            render_dir (str): path to rendered roof images
            azimuth (float): azimuth orientation of the array
            inclination (float): inclination of the array
            resolution_I (int): Current (A) resolution for the simulation:
                Defaults to 110.
            start_year (int): Initial year for the simulation.
                Defaults to 2007.
            end_year (int): Last year for the simulation. Defaults to 2008.

        """
        self.shading = shading
        self.building = building
        self.pv_panel = pv_panel
        self.shading.building = building
        self.azimuth = azimuth
        self.inclination = inclination

        self.resolution_I = 110
        self.offset_x = 0.0
        self.offset_y = 0.0
        self.grid = None
        self.start_year = start_year
        self.end_year = end_year
        self.timestep = building.timestep
        self.output_dir = None

    def set_offsets(self, offset_x: float = 0, offset_y: float = 0):
        """Set offsets of the solar array grid root.

        Sets the the offsets of the solar array with respect to the Shading
        images. With respect to the size of the used panel, the scale of the
        images and the offsets this function also creates an empty grid for
        possible panel locations

        Args:
            offset_x (float): Offset of the origin of the array along the
                x-Axis from the top left corner of the image. Defaults to 0.
            offset_y (float): Offset of the origin of the array along the
                y-Axis from the top left corner of the image. Defaults to 0.
        """
        self.offset_x = offset_x
        self.offset_y = offset_y
        self.grid = -1*np.ones(self.grid_shape).T

    def preprocess(self):
        """Preprocess the solar_array using provided information."""
        self.gen_panel_slices()
        self.get_resize_shading_fun()
        self.Pvgis = Pvgis(loc=self.building.location, azimuth=self.azimuth,
                           inclination=self.inclination,
                           startyear=self.start_year, endyear=self.end_year)

    def show_array(self,):
        """Preview Array configuration on roof using matplotlib."""
        shading = self.shading
        pv_panel = self.pv_panel
        offset_x, offset_y = self.offset_x, self.offset_y
        fig, ax = plt.subplots(figsize=(14, 10))
        ax.imshow(shading.avg, extent=(
            0, shading.dim_x, shading.dim_y, 0), alpha=.7)
        if self.grid is None:
            grid_shape = self.grid_shape
            grid = -1*np.ones(grid_shape).T
            coords = np.argwhere(grid > -2)
        else:
            grid = self.grid
        coords = np.argwhere(grid > -2)
        ax.plot(self.offset_x, self.offset_y, 'x', markersize=20, color='red')
        for x, y in coords:
            string = int(grid[x, y])
            alpha = .5+.5*int(string > -1)
            rect = Rectangle((pv_panel.dim_x*x+offset_x,
                              pv_panel.dim_y*y+offset_y),
                             pv_panel.dim_x, pv_panel.dim_y,
                             fill=None, alpha=alpha)
            bound = rect.get_extents()
            center = np.average(bound, axis=0)
            ax.add_patch(rect)
            ax.text(*center, '({},{})\nStr: {}'.format(x, y, string),
                    horizontalalignment='center',
                    verticalalignment='center',
                    fontsize=12, color='black', alpha=alpha)
        return fig, ax

    @property
    def grid_shape(self):
        """Return a maximum possible grid shape.

        Depends on size and scale of the image and panel.
        """
        return (int((self.shading.dim_y-self.offset_y)//self.pv_panel.dim_y),
                int((self.shading.dim_x-self.offset_x)//self.pv_panel.dim_x))

    @property
    def cell_grid(self):
        """Return the actual cell grid shape depending on defined strings."""
        pv_p = self.pv_panel
        coords = np.argwhere(self.grid > 0)
        x0, y0 = coords.min(axis=0)
        x_max, y_max = coords.max(axis=0)+1
        return(self.grid[x0:x_max, y0:y_max].repeat(pv_p.cells_x, axis=0
                                                    ).repeat(
                                                        pv_p.cells_y, axis=1))

    def get_resize_shading_fun(self):
        """Set the image resize function."""
        assert not(self.grid is None), 'Grid needs to be defined first'
        cell_grid_shape = self.cell_grid.shape
        ind_y_start = int(self.offset_y*self.shading.scale)
        ind_x_start = int(self.offset_x*self.shading.scale)
        ind_y_delta = int(
            self.cell_grid.shape[1]*self.pv_panel.dim_y *
            self.shading.scale/self.pv_panel.cells_y)
        ind_x_delta = int(
            self.cell_grid.shape[0]*self.pv_panel.dim_x *
            self.shading.scale/self.pv_panel.cells_x)

        def resize(img):
            return zoom(img[ind_y_start:ind_y_start+ind_y_delta,
                            ind_x_start:ind_x_start+ind_x_delta, ],
                        zoom=[cell_grid_shape[1]/ind_y_delta,
                              cell_grid_shape[0]/ind_x_delta],
                        order=0)
        self.resize = resize

    def gen_panel_slices(self):
        """Set list of slices for each cell on a substring."""
        module_indices = []
        panels = np.argwhere(self.grid > 0)
        pv_panel = self.pv_panel
        for pos_x, pos_y in panels:
            x_min = pos_x*pv_panel.cells_x
            x_max = x_min
            substrings = []
            for bypass_x in pv_panel.bypass_config[0]:
                x_max += bypass_x
                y_min = pos_y*pv_panel.cells_y
                y_max = y_min
                if isinstance(pv_panel.bypass_config[1], int):
                    for bypass_y in [pv_panel.bypass_config[1]]:
                        y_max += bypass_y
                        substrings.append(
                            [slice(x_min, x_max), slice(y_min, y_max)])
                        y_min = y_max
                else:
                    for bypass_y in pv_panel.bypass_config[1]:
                        y_max += bypass_y
                        substrings.append(
                            [slice(x_min, x_max), slice(y_min, y_max)])
                        y_min = y_max
                x_min = x_max
            module_indices.append(substrings)
        self.module_indices = module_indices

    def show_panel_slices(self):
        """Short summary.

        Returns:
            type: .

        """
        data = np.zeros(self.cell_grid.T.shape)
        for ind, mod in enumerate(self.module_indices[:]):
            for ind2, substring in enumerate(mod):
                data[substring[1], substring[0]] += (ind2+1+ind*.3)
        plt.imshow(data)
        return data

    def run(self, Pvgis, start_date=None, end_date=None, n_threads=8,
            batch_size=24*6, ):
        """Short summary.

        Args:
            Pvgis (Pvgis): Pvgis data for the simulation.
            start_date (dt.datetime): Simulation startdate. Defaults to None.
            end_date (dt.datetime): . Simulation enddate to None.
            n_threads (int): . Defaults to 8.
            batch_size (int): . Defaults to 24*6.

        Returns:
            type: .

        """
        if start_date is None:
            start_date = dt.datetime(self.start_year, 1, 2)
        if end_date is None:
            end_date = dt.datetime(self.end_year, 1, 1)

        timestep = self.timestep
        steps = int((end_date.timestamp()-start_date.timestamp()) /
                    timestep/batch_size)
        batches = []
        for time in range(steps):
            curr_time = start_date + \
                dt.timedelta(seconds=time*timestep*batch_size)
            end_time = start_date + \
                dt.timedelta(seconds=(time+1)*timestep*batch_size)
            batches.append((curr_time, end_time))
        self.results = self.process(batches[0])
        # with Pool(n_threads) as p:
        #    self.batch_data = p.map(self.process, self.batches)

    def process(self, args):
        """Process function for use in multiprocessing.Pool.

        Args:
            args (List[dt.datetime, dt.datetime]): startdate and enddate of the
                simulation.

        Returns:
            type: .

        """
        timestep = self.timestep
        (start_date, end_date) = args
        steps = int(end_date.timestamp()-start_date.timestamp())//timestep
        I_conf = np.linspace(self.pv_panel.I_min,
                             self.pv_panel.I_max, self.resolution_I)
        n_modules = len(self.module_indices)
        V_I_data = np.zeros((steps, n_modules, self.resolution_I))
        for time in range(int(end_date.timestamp() -
                              start_date.timestamp())//timestep):
            curr_time = start_date+dt.timedelta(seconds=time*timestep)
            environment = self.Pvgis[curr_time]
            L_mult = environment['Bi']
            L_add_di = environment['Di']+environment['Ri']+environment['As']
            T_am = environment['Tamb']
            # Define Iph/G func more accurately (Right now it's just
            # a factor of 1e-3)
            W_10 = environment['W10']
            if L_add_di == 0:
                continue
            roof = self.resize(self.shading[curr_time])
            for ind, module in enumerate(self.module_indices):
                L_ = [np.reshape(roof[substring[0],
                                      substring[1]], -1
                                 )*L_mult*1e-3 +
                      L_add_di*1e-3 for substring in module]
                V_I_data[time, ind] = self.pv_panel.v_from_i(
                    (I_conf, L_, T_am, W_10, -.5))
            np.save("/home/clemens/data/VI2/" +
                    start_date.strftime('%Y-%m-%d'), V_I_data)
        return [start_date, end_date,  V_I_data]

max_stars_repo_name: cbaretzky/BuiidingEnergySimulation

max_issues_repo_name: cbaretzky/BuiidingEnergySimulation

max_forks_repo_name: cbaretzky/BuiidingEnergySimulation

reason: import numpy,from scipy


text:
"""
Copyright (C) 2017 University of Massachusetts Amherst.
This file is part of "xcluster"
http://github.com/iesl/xcluster
Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
"""

"""
Utilities for creating delta separated data sets.
"""
import numpy as np

def gen_k_centers(k, dim):
    """Create a k cluster data set with required separation.

    For the purposes of validating a proof, generate each cluster center such
    that it is at least 4 * delta away from any other cluster for some value of
    delta > 0.

    Args:
        k - the number of clusters.
        dim - (optional) the dimension of the points.

    Returns:
        A list of 2 cluster centers and a value of delta such that the clusters
        centers are 4 * delta away form each other.
    """
    delta = abs(np.random.normal(0.0, 5.0))
    eps = 0.001
    centers = []
    for i in range(k):
        c = np.random.multivariate_normal(np.zeros(dim), np.identity(dim))
        if len(centers):
            c1 = centers[0]
            x = np.random.multivariate_normal(c1, np.identity(c1.size)) - c1
            direction = x / np.linalg.norm(x)
            centers.append(c1 + 2.0 * i * delta * direction + eps)
        else:
            centers.append(c)
    return centers, delta

def create_dataset(dims, size, num_clusters=20):
    """Create a delta separated data set.

    Generate a set of centers for the clusters and from each center draw size
    number of points that constitute the points in that cluster.  Then return
    a dataset of all points.

    Args:
        dims - (int) the dimention of all data points.
        size - (int) the number of points to generate for each cluster.
        num_clusters - (int) the number of clusters.
    """
    clusters, delta = gen_k_centers(num_clusters, dims)
    return _create_constrained_dataset(clusters, delta, size)

def _create_constrained_dataset(centers, delta, size):
    """Create a delta-separated dataset.

    For each of the centers draw size number of points. No two points may be
    farther than delta away form each other. Thus, to generate each point,
    choosea random direction and random distance from the center (of up to 0.5
    delta).

    Args:
      centers - a list of cluster centers.
      delta - the maximum distance between two points in the same cluster.
      size - the number of points to draw per cluster.

    Returns:
      A list of points that represents the dataset.
    """
    dataset = []
    count = 0
    for i, c in enumerate(centers):
        for j in range(size):
            x = np.random.multivariate_normal(c, np.identity(np.size(c))) - c
            direction = x / np.linalg.norm(x)
            magnitude = np.random.uniform(0.0, 0.5 * delta)
            # magnitude = np.random.uniform(0.0, delta) # NOT DEL-SEPARATED
            dataset.append((c + magnitude * direction, i, count))
            count += 1
    return dataset

def gen_4_normal():
    """Create 4 cluster centers.

    Create gaussians centered at (1,1), (1,-1), (-1,-1) and (-1,1).  Each has
    standard covariance.

    Args:
        None

    Returns:
        A list of the four cluster centers.
    """
    return [mn(mean=np.array([1.0, 1.0]),
               cov=np.array([[1.0, 0.0], [0.0, 1.0]])),
            mn(mean=np.array([1.0, -1.0]),
               cov=np.array([[1.0, 0.0], [0.0, 1.0]])),
            mn(mean=np.array([-1.0, -1.0]),
               cov=np.array([[1.0, 0.0], [0.0, 1.0]])),
            mn(mean=np.array([-1.0, 1.0]),
               cov=np.array([[1.0, 0.0], [0.0, 1.0]]))]

def _4_normal_spread():
  """Create 4 cluster centers.

  Create gaussians centered at (10,10), (10,-10), (-10,-10) and (-10,10).
  Each has standard covariance.

  Args:
    None

  Returns:
    A list of the four cluster centers.
  """
  return [mn(mean=np.array([10.0, 10.0]),
             cov=np.array([[1.0, 0.0], [0.0, 1.0]])),
          mn(mean=np.array([10.0, -10.0]),
             cov=np.array([[1.0, 0.0], [0.0, 1.0]])),
          mn(mean=np.array([-10.0, -10.0]),
             cov=np.array([[1.0, 0.0], [0.0, 1.0]])),
          mn(mean=np.array([-10.0, 10.0]),
             cov=np.array([[1.0, 0.0], [0.0, 1.0]]))]

def _5x5_grid_clusters():
  """Create a 5x5 grid of cluster centers.

  Create 25 cluster centers on the grid I^{[0, 4] x [0,4]}.  Each center is a
  gaussian with standard covariance

  Args:
    None

  Returns:
    A list of cluster centers.
  """
  return [mn(mean=np.array([i, j]), cov=np.array([[1.0, 0.0],
                                                  [0.0, 1.0]]))
          for i in range(5)
          for j in range(5)]

def _5x5_grid_clusters_spread():
  """Create a 5x5 grid of cluster centers.

  Create 25 cluster centers on the grid I^{[0, 4] x [0,4]}.  Each center is a
  gaussian with standard covariance

  Args:
    None

  Returns:
    A list of cluster centers.
  """
  return [mn(mean=np.array([i * 25, j * 25]), cov=np.array([[1.0, 0.0],
                                                            [0.0, 1.0]]))
          for i in range(5)
          for j in range(5)]

def _5x5_grid_clusters_close():
  """Create a 5x5 grid of cluster centers.

  Create 25 cluster centers on the grid I^{[0, 4] x [0,4]}.  Each center is a
  gaussian with standard covariance

  Args:
    None

  Returns:
    A list of cluster centers.
  """
  return [mn(mean=np.array([i * 5, j * 5]), cov=np.array([[1.0, 0.0],
                                                          [0.0, 1.0]]))
          for i in range(5)
          for j in range(5)]

def _2x3_grid_clusters_close():
  """Create a 3x3 grid of cluster centers.

  Create 25 cluster centers on the grid I^{[0, 4] x [0,4]}.  Each center is a
  gaussian with standard covariance

  Args:
    None

  Returns:
    A list of cluster centers.
  """
  return [mn(mean=np.array([i * 5, j * 5]), cov=np.array([[1.0, 0.0],
                                                          [0.0, 1.0]]))
          for i in range(2)
          for j in range(3)]

def _2x3_grid_clusters_spread():
  """Create a 3x3 grid of cluster centers.

  Create 25 cluster centers on the grid I^{[0, 4] x [0,4]}.  Each center is a
  gaussian with standard covariance

  Args:
    None

  Returns:
    A list of cluster centers.
  """
  return [mn(mean=np.array([i * 25, j * 25]), cov=np.array([[1.0, 0.0],
                                                            [0.0, 1.0]]))
          for i in range(2)
          for j in range(3)]

def _10x10_grid_clusters_close():
  """Create a 3x3 grid of cluster centers.

  Create 25 cluster centers on the grid I^{[0, 4] x [0,4]}.  Each center is a
  gaussian with standard covariance

  Args:
    None

  Returns:
    A list of cluster centers.
  """
  return [mn(mean=np.array([i * 5, j * 5]), cov=np.array([[1.0, 0.0],
                                                          [0.0, 1.0]]))
          for i in range(10)
          for j in range(10)]

def _10x10_grid_clusters_spread():
  """Create a 3x3 grid of cluster centers.

  Create 25 cluster centers on the grid I^{[0, 4] x [0,4]}.  Each center is a
  gaussian with standard covariance

  Args:
    None

  Returns:
    A list of cluster centers.
  """
  return [mn(mean=np.array([i * 25, j * 25]), cov=np.array([[1.0, 0.0],
                                                            [0.0, 1.0]]))
          for i in range(10)
          for j in range(10)]

def _random_standard_centers(n=100):
  """Create random cluster centers.

  Create n cluster centers randomly.  Each cluster center is a draw from a
  gaussian distribution centered at (0,0) with standard covariance.

  Args:
    n - optional; the number of centers to draw (default 100).

  Returns:
    A list of cluster centers.
  """
  generator = mn(mean=np.array([0, 0]),
                 cov=np.array([[1.0, 0.0], [0.0, 1.0]]))
  return [mn(mean=pt, cov=np.array([[1.0, 0.0], [0.0, 1.0]]))
          for pt in generator.rvs(size=n)]

def _from_file(filename):
  with open(filename, 'r') as f:
    clustering = []
    for line in f:
      splits = line.split('\t')
      l, vec = int(splits[0]), np.array([float(x) for x in splits[1:]])
      clustering.append((vec, l))
  return clustering

max_stars_repo_name: iesl/xcluster

max_issues_repo_name: iesl/xcluster

max_forks_repo_name: iesl/xcluster


text:
import os
import pickle
from typing import List

import numpy as np
from PIL import Image

def unpickle(file):
    with open(file, "rb") as fo:
        dict = pickle.load(fo, encoding="bytes")
    return dict

def parse_pickle(rawdata, rootdir) -> List[List[str]]:
    for i in range(10):
        directory = f"{rootdir}/{i}"
        os.makedirs(directory, exist_ok=True)
    class_to_filename_list = []
    for i in range(len(rawdata[b"filenames"])):
        filename = rawdata[b"filenames"][i].decode("utf-8")
        label = rawdata[b"labels"][i]
        data = rawdata[b"data"][i]
        data = data.reshape(3, 32, 32)
        data = np.swapaxes(data, 0, 2)
        data = np.swapaxes(data, 0, 1)
        with Image.fromarray(data) as img:
            image_path = f"{rootdir}/{label}/{filename}"
            img.save(image_path)
        class_to_filename_list.append([label, filename])
    return class_to_filename_list

max_stars_repo_name: sudabon/ml-system-in-actions

max_issues_repo_name: sudabon/ml-system-in-actions

max_forks_repo_name: sudabon/ml-system-in-actions