class ParameterSet(dict):
    """
    Class to manage hierarchical parameter sets.

    Usage example:

        > sim_params = ParameterSet({'dt': 0.1, 'tstop': 1000.0})
        > exc_cell_params = ParameterSet({'tau_m': 20.0, 'cm': 0.5})
        > inh_cell_params = ParameterSet({'tau_m': 15.0, 'cm': 0.5})
        > network_params = ParameterSet({'excitatory_cells': exc_cell_params, 'inhibitory_cells': inh_cell_params})
        > P = ParameterSet({'sim': sim_params, 'network': network_params})
        > P.sim.dt
        0.1
        > P.network.inhibitory_cells.tau_m
        15.0
    """

    @staticmethod
    def read_from_file(pth, filename):
        """
        Import parameter dictionary stored as a text file. The
        file must be in standard format, i.e. {'key1': value1, ...}
        :param pth: path
        :param filename:
        :return:
        """

        assert os.path.exists(pth), "Incorrect path..."

        with open('{}/{}'.format(pth, filename), 'r') as fi:
            contents = fi.read()
        d = eval(contents)

        return d

    def __init__(self, initializer, label='global'):
        """

        :param initializer: parameters dictionary, or string locating the full path of text file with
        parameters dictionary written in standard format
        :param label: label the parameters set
        :return: ParameterSet
        """

        def convert_dict(d, label):
            """
            Iterate through the dictionary d, replacing all items by ParameterSet objects
            :param d:
            :param label:
            :return:
            """
            for k, v in d.items():
                if isinstance(v, ParameterSet):
                    d[k] = v
                elif isinstance(v, dict):
                    d[k] = convert_dict(d, k)
                else:
                    d[k] = v
            return ParameterSet(d, label)

        # self._url = None
        if isinstance(initializer, basestring):  # url or str
            if os.path.exists(initializer):
                pth = os.path.dirname(initializer)
                if pth == '':
                    pth = os.path.abspath('.')
                filename = ntpath.basename(initializer)
                initializer = ParameterSet.read_from_file(pth, filename)

        # initializer is now a dictionary (if it was a path before), and we can iterate it and replace all items by
        # ParameterSets
        if isinstance(initializer, dict):
            for k, v in initializer.items():
                if isinstance(v, ParameterSet):
                    self[k] = v

                elif isinstance(v, dict):
                    self[k] = ParameterSet(v)

                else:
                    self[k] = v
        else:
            raise TypeError("initializer must be either a string specifying "
                            "the full path of the parameters file, "
                            "or a parameters dictionary")

        # set label
        if isinstance(initializer, dict):
            if 'label' in initializer:
                self.label = initializer['label']
            else:
                self.label = label
        elif hasattr(initializer, 'label'):
            self.label = label or initializer.label
        else:
            self.label = label

    def flat(self):
        __doc__ = nesteddict_walk.__doc__
        return nesteddict_walk(self)

    def flatten(self):
        __doc__ = nesteddict_flatten.__doc__
        return nesteddict_flatten(self)

    def __eq__(self, other):
        """
        Simple function for equality check of ParameterSet objects. Check is done implicitly by converting the objects
        to dictionaries.
        :param other:
        :return:
        """
        # compare instances
        if not isinstance(self, other.__class__):
            return False

        ad, od = self.as_dict(), other.as_dict()
        for key in ad:
            if ad[key] != od[key]:
                return False
        return True

    def __getattr__(self, name):
        """
        Allow accessing parameters using dot notation.
        """
        try:
            return self[name]
        except KeyError:
            return self.__getattribute__(name)

    def __setattr__(self, name, value):
        """
        Allow setting parameters using dot notation.
        """
        self[name] = value

    def __getitem__(self, name):
        """
        Modified get that detects dots '.' in the names and goes down the
        nested tree to find it
        """
        split = name.split('.', 1)
        if len(split) == 1:
            return dict.__getitem__(self, name)
        # nested get
        return dict.__getitem__(self, split[0])[split[1]]

    def flat_add(self, name, value):
        """
        Like `__setitem__`, but it will add `ParameterSet({})` objects
        into the namespace tree if needed.
        """

        split = name.split('.', 1)
        if len(split) == 1:
            dict.__setitem__(self, name, value)
        else:
            # nested set
            try:
                ps = dict.__getitem__(self, split[0])
            except KeyError:
                # setting nested name without parent existing
                # create parent
                ps = ParameterSet({})
                dict.__setitem__(self, split[0], ps)
            # and try again
            ps.flat_add(split[1], value)

    def __setitem__(self, name, value):
        """
        Modified set that detects dots '.' in the names and goes down the
        nested tree to set it
        """
        if isinstance(name, str):
            split = name.split('.', 1)
            if len(split) == 1:
                dict.__setitem__(self, name, value)
            else:
                # nested set
                dict.__getitem__(self, split[0])[split[1]] = value

    def update(self, E, **F):
        """
        Update ParameterSet with dictionary entries
        """
        if hasattr(E, "has_key"):
            for k in E:
                self[k] = E[k]
        else:
            for (k, v) in E:
                self[k] = v
        for k in F:
            self[k] = F[k]

    def __getstate__(self):
        """
        For pickling.
        """
        return self

    def save(self, url=None):
        """
        Write the ParameterSet to a text file
        The text format should be valid python code...
        :param url: locator (full path+ filename, str) - if None, data will be saved to self._url
        """

        if not url:
            ("Please provide url")
        # url = self._url
        assert url != ''
        # if not self._url:
        # 	self._url = url

        with open(url, 'w') as fp:
            fp.write(self.pretty())

    def pretty(self, indent='   '):
        """
        Return a unicode string representing the structure of the `ParameterSet`.
        evaluating the string should recreate the object.
        :param indent: indentation type
        :return: string
        """

        def walk(d, indent, ind_incr):
            s = []
            for k, v in d.items():
                if isinstance(v, list):
                    if k == 'randomize_neuron_pars' and np.mean([isinstance(x, dict) for x in v]) == 1.:
                        s.append('%s"%s": [' % (indent, k))
                        for x in v:
                            s.append('{')
                            s.append(walk(x, indent + ind_incr, ind_incr))
                            s.append('},\n')
                        s.append('],\n')
                        v_arr = np.array([])
                        continue
                    elif np.mean([isinstance(x, list) for x in v]):
                        v_arr = np.array(list(itertools.chain(*v)))
                    else:
                        v_arr = np.array(v)
                    mark = np.mean([isinstance(x, tuple) for x in v])
                    if (v_arr.dtype == object) and mark:
                        if len(v) == 1:
                            if v and isinstance(v[0], tuple) and isinstance(v[0][0], types.BuiltinFunctionType):
                                if isinstance(v_arr.any(), types.BuiltinFunctionType):
                                    if v_arr.any().__name__ in np.random.__all__:
                                        s.append('%s"%s": [(%s, %s)],' % (indent, k, 'np.random.{0}'.format(v_arr.any(
                                        ).__name__), str(v_arr.all())))
                                elif isinstance(v_arr.any(), types.MethodType) and v_arr.any().__str__().find('scipy'):
                                    s.append('%s"%s": [(%s, %s)],' % (indent, k, 'st.{0}.rvs'.format(str(v_arr.any().im_self.name)),
                                                                      str(v_arr.all())))
                        else:
                            if v and isinstance(v_arr.any(), tuple):
                                if isinstance(v_arr.any()[0], types.BuiltinFunctionType):
                                    list_idx = np.where(v_arr.any() in v)[0][0]
                                    string = '%s"%s": [' % (indent, k)
                                    for idx, nnn in enumerate(v):
                                        if idx == list_idx:
                                            tmp = np.array(nnn)
                                            string += '(%s, %s), ' % ('np.random.{0}'.format(tmp.any().__name__), str(tmp.all()))
                                        elif nnn is not None and isinstance(nnn[0], types.BuiltinFunctionType):
                                            tmp = np.array(nnn)
                                            string += '(%s, %s), ' % ('np.random.{0}'.format(tmp.any().__name__), str(tmp.all()))
                                        else:
                                            string += '%s, ' % (str(nnn))
                                    string += '], '
                                    s.append(string)
                            elif v and isinstance(v_arr.any(), types.BuiltinFunctionType) and isinstance(v_arr.all(),
                                                                                                         dict):
                                string = '%s"%s": [' % (indent, k)
                                for idx, nnn in enumerate(v):
                                    tmp = np.array(nnn)
                                    string += '(%s, %s), ' % (
                                        'np.random.{0}'.format(tmp.any().__name__), str(tmp.all()))
                                string += '], '
                                s.append(string)
                    else:
                        if np.mean([isinstance(x, np.ndarray) for x in v]):
                            string = '%s"%s": [' % (indent, k)
                            for idx, n in enumerate(v):
                                if isinstance(n, np.ndarray):
                                    string += 'np.' + repr(n) + ', '
                                else:
                                    string += str(n) + ', '
                            string += '], '
                            s.append(string)
                        elif hasattr(v, 'items'):
                            s.append('%s"%s": {' % (indent, k))
                            s.append(walk(v, indent + ind_incr, ind_incr))
                            s.append('%s},' % indent)
                        elif isinstance(v, basestring):
                            s.append('%s"%s": "%s",' % (indent, k, v))
                        else:
                            # what if we have a dict or ParameterSet inside a list? currently they are not expanded.
                            # Should they be?
                            s.append('%s"%s": %s,' % (indent, k, v))

                elif isinstance(v, types.BuiltinFunctionType):
                    if v.__name__ in np.random.__all__:
                        s.append('%s"%s"' % (indent, 'np.random.{0}'.format(v.__name__)))
                else:
                    if hasattr(v, 'items'):
                        s.append('%s"%s": {' % (indent, k))
                        s.append(walk(v, indent + ind_incr, ind_incr))
                        s.append('%s},' % indent)
                    elif isinstance(v, basestring):
                        s.append('%s"%s": "%s",' % (indent, k, v))
                    elif isinstance(v, np.ndarray):
                        s.append('%s"%s": %s,' % (indent, k, repr(v)[6:-1]))
                    elif isinstance(v, tuple) and (isinstance(v[0], types.MethodType)
                                                   or isinstance(v[0], types.BuiltinFunctionType)):
                        v_arr = np.array(v)
                        if isinstance(v[0], types.MethodType):
                            s.append('%s"%s": (%s, %s),' % (indent, k, 'st.{0}.rvs'.format(
                                str(v_arr.any().im_self.name)), str(v_arr.all())))
                        elif isinstance(v[0], types.BuiltinFunctionType):
                            s.append('%s"%s": (%s, %s),' % (indent, k, 'np.random.{0}'.format(v_arr.any().__name__),
                                                            str(v_arr.all())))
                    else:
                        # what if we have a dict or ParameterSet inside a list? currently they are not expanded.
                        #  Should they be?
                        s.append('%s"%s": %s,' % (indent, k, v))
            return '\n'.join(s)

        return '{\n' + walk(self, indent, indent) + '\n}'

    def tree_copy(self):
        """
        Return a copy of the `ParameterSet` tree structure.
        Nodes are not copied, but re-referenced.
        """

        tmp = ParameterSet({})
        for key in self:
            value = self[key]
            if isinstance(value, ParameterSet):
                tmp[key] = value.tree_copy()
            else:
                tmp[key] = value
        return tmp

    def as_dict(self):
        """
        Return a copy of the `ParameterSet` tree structure as a nested dictionary.
        """
        tmp = {}

        for key in self:
            if not isinstance(self[key], types.BuiltinFunctionType):
                value = self[key]
                if isinstance(value, ParameterSet):
                    tmp[key] = value.as_dict()
                else:
                    tmp[key] = value
        return tmp

    def replace_values(self, **args):
        """
        This expects its arguments to be in the form path=value, where path is a
        . (dot) delimited path to a parameter in the  parameter tree rooted in
        this ParameterSet instance.

        This function replaces the values of each parameter in the args with the
        corresponding values supplied in the arguments.
        """
        for k in args.keys():
            self[k] = args[k]

    def clean(self, termination='pars'):
        """
        Remove fields from ParameterSet that do not contain the termination
        This is mostly useful if, in the specification of the parameters
        additional auxiliary variables were set, and have no relevance for
        the experiment at hand...
        """
        accepted_keys = [x for x in self.iterkeys() if x[-len(termination):] == termination]
        new_dict = {k: v for k, v in self.iteritems() if k in accepted_keys}

        return ParameterSet(new_dict) class ParameterSpace:
    """
    A collection of `ParameterSets`, representing multiple points (combinations) in parameter space.
    Parses parameter scripts and runs experiments locally or creates job files to be run on a cluster.
    Can also harvest stored results from previous experiments and recreate the parameter space for post-analysis.
    """

    def __init__(self, initializer, keep_all=False):
        """
        Generate ParameterSpace containing a list of all ParameterSets

        :param initializer: file url
        :param keep_all: store all original parameters (??)
        :return: (tuple) param_sets, param_axes, global_label, size of parameter space
        """
        assert isinstance(initializer, str), "Filename must be provided"
        with open(initializer, 'r') as fp:
            self.parameter_file = fp.readlines()

        def validate_parameters_file(module):
            """
            Checks for any errors / incompatibilities in the structure of the parameter file. Function
            `build_parameters` is required, with or without arguments.
            :param module: imported parameter module object
            :return:
            """
            # TODO anything else?
            # build_parameters function must be defined!
            if not hasattr(module, "build_parameters"):
                raise ValueError("`build_parameters` function is missing!")

            # check parameter range and function arguments
            range_arguments = inspect.getargspec(module.build_parameters)[0]
            if not hasattr(module, "parameter_range"):
                if len(range_arguments) == 0:
                    return
                raise ValueError("`parameter_range` and arguments of `build_parameters` do not match!")

            for arg in range_arguments:
                if arg not in module.parameter_range:
                    raise ValueError('ParameterRange variable `%s` is not in `parameter_range` dictionary!' % arg)
                if not isiterable(module.parameter_range[arg]):
                    raise ValueError('ParameterRange variable `%s` is not iterable! Should be list!' % arg)

        # TODO this needs some rethinking, possibly unnecessary
        # def validate_parameter_sets(param_sets):
        #     """
        #
        #     :param param_sets:
        #     :return:
        #     """
        #     required_dicts = ["kernel_pars", "encoding_pars", "net_pars"]
        #     for p_set in param_sets:
        #         for d in required_dicts:
        #             if d not in p_set:
        #                 raise ValueError("Required parameter (dictionary) `%s` not found!" % d)
        #
        #         # `data_prefix` required
        #         if "data_prefix" not in p_set["kernel_pars"]:
        #             raise ValueError("`data_prefix` missing from `kernel_pars`!")

        def parse_parameters_file(url):
            """

            :param url:
            :return:
            """
            module_name = import_mod_file(url)
            module_obj 	= globals()[module_name]
            try:
                validate_parameters_file(module_obj)
            except ValueError as error:
                logger.error("Invalid parameter file! Error: %s" % error)
                exit(-1)

            range_args	= inspect.getargspec(module_obj.build_parameters)[0]  # arg names in build_parameters function
            n_ranges 	= len(range_args)
            n_runs 	 	= int(np.prod([len(module_obj.parameter_range[arg]) for arg in range_args]))  # nr combinations
            param_axes 	= dict()
            param_sets 	= []

            # build cross product of parameter ranges, sorting corresponds to argument ordering in build_parameters
            # empty if no range parameters are present: [()]
            range_combinations = list(itertools.product(*[module_obj.parameter_range[arg] for arg in range_args]))
            # call build_parameters for each range combination, and pack returned values into a list (set)
            # contains a single dictionary if no range parameters
            param_ranges = [module_obj.build_parameters( *elem ) for elem in range_combinations]
            global_label = param_ranges[0]['kernel_pars']['data_prefix']

            # currently only support parameter spaces of dimensions <=3
            if n_ranges <= 3:# and not emoo:
                # # verify parameter integrity / completeness
                # try:
                #     validate_parameter_sets(param_ranges)
                # except ValueError as error:
                #     print("Invalid parameter file! Error: %s" % error)

                # build parameter axes
                axe_prefixes = ['x', 'y', 'z']
                for range_index in range(n_runs):
                    params_label = global_label
                    for axe_index, axe in enumerate(axe_prefixes[:n_ranges]):  # limit axes to n_ranges
                        param_axes[axe + 'label'] 		= range_args[axe_index]
                        param_axes[axe + 'ticks'] 		= module_obj.parameter_range[param_axes[axe + 'label']]
                        param_axes[axe + 'ticklabels'] 	= [str(xx) for xx in param_axes[axe + 'ticks']]
                        params_label += '_' + range_args[axe_index]
                        params_label += '=' + str(range_combinations[range_index][axe_index])

                    p_set = ParameterSet(param_ranges[range_index], label=params_label)
                    if not keep_all:
                        p_set = p_set.clean(termination='pars')
                    p_set.update({'label': params_label})
                    param_sets.append(p_set)
            else:
                raise ValueError("Parameter spaces of >3 dimensions are currently not supported")

            return param_sets, param_axes, global_label, n_ranges

        def parse_parameters_dict(url):
            """
            Simple parser for dictionary (text) parameter scripts, not .py!
            :param url:
            :return:
            """
            # TODO is set label always global?
            param_set 	= ParameterSet(url, label='global')
            # try:
            #     validate_parameter_sets([param_set])
            # except ValueError as error:
            #     print("Invalid parameter file! Error: %s" % error)

            param_axes	= {}
            label 		= param_set.kernel_pars["data_prefix"]
            dim			= 1
            return [param_set], param_axes, label, dim

        if initializer.endswith(".py"):
            self.parameter_sets, self.parameter_axes, self.label, self.dimensions = parse_parameters_file(initializer)
        else:
            self.parameter_sets, self.parameter_axes, self.label, self.dimensions = parse_parameters_dict(initializer)

    def update_run_parameters(self, cluster=None):
        """
        Update run type and experiment specific paths in case a cluster template is specified.

        :param cluster: name of cluster template, e.g., Blaustein
        :return:
        """
        if cluster is not None:
            assert cluster in paths.keys(), "Default setting for cluster {0} not found!".format(cluster)
            for param_set in self.parameter_sets:
                param_set.kernel_pars['data_path'] = paths[cluster]['data_path']
                param_set.kernel_pars['mpl_path'] = paths[cluster]['matplotlib_rc']
                param_set.kernel_pars['print_time'] = False
                param_set.kernel_pars.system.local = False
                param_set.kernel_pars.system.system_label = cluster
                param_set.kernel_pars.system.jdf_template = paths[cluster]['jdf_template']
                param_set.kernel_pars.system.remote_directory = paths[cluster]['remote_directory']
                param_set.kernel_pars.system.queueing_system = paths[cluster]['queueing_system']

    def iter_sets(self):
        """
        An iterator which yields the ParameterSets in ParameterSpace
        """
        for val in self.parameter_sets:
            yield val

    def __eq__(self, other):
        """
        For testing purposes

        :param other:
        :return:
        """
        if not isinstance(other, self.__class__):
            return False

        for key in self.parameter_axes:
            if key not in other.parameter_axes:
                return False
            else:
                if isinstance(self.parameter_axes[key], np.ndarray):
                    if not np.array_equal(self.parameter_axes[key], other.parameter_axes[key]):
                        return False
                elif self.parameter_axes[key] != other.parameter_axes[key]:
                    return False

        if self.label != other.label or self.dimensions != other.dimensions:
            return False

        return self.parameter_sets == other.parameter_sets

    def __getitem__(self, item):
        return self.parameter_sets[item]

    def __len__(self):
        return len(self.parameter_sets)

    def save(self, target_full_path):
        """
        Save the full ParameterSpace by re-writing the parameter file
        :return:
        """
        with open(target_full_path, 'wb') as fp:
            fp.writelines(self.parameter_file)

    def compare_sets(self, parameter):
        """
        Determine whether a given parameter is common to all parameter sets
        :param parameter: parameter to compare
        """
        common = dict(pair for pair in self.parameter_sets[0].items() if all((pair in d.items() for d in
                                                                              self.parameter_sets[1:])))
        result = False
        if parameter in common.keys():
            result = True
        else:
            for k, v in common.items():
                if isinstance(v, dict):
                    if parameter in v.keys():
                        result = True
        return result

    def run(self, computation_function, project_dir=None, **parameters):
        """
        Run a computation on all the parameters.

        :param computation_function: function to execute
        :param project_dir: project directory
        :param parameters: kwarg arguments for the function
        """
        system = self.parameter_sets[0].kernel_pars.system

        if system['local']:
            logger.info("\nRunning {0} serially on {1} Parameter Sets".format(
                str(computation_function.__module__.split('.')[1]), str(len(self))))

            results = None
            for par_set in self.parameter_sets:
                logger.info("\n- Parameters: {0}".format(str(par_set.label)))
                results = computation_function(par_set, **parameters)
            return results
        else:
            logger.info("\nPreparing job description files...")
            export_folder 			= system['remote_directory']
            main_experiment_folder 	= export_folder + '{0}/'.format(self.label)

            try:
                os.makedirs(main_experiment_folder)
            except OSError as err:
                if err.errno == errno.EEXIST and os.path.isdir(main_experiment_folder):
                    logger.info("Path `{0}` already exists, will be overwritten!".format(main_experiment_folder))
                else:
                    raise OSError(err.errno, "Could not create exported experiment folder.", main_experiment_folder)

            remote_run_folder 	  = export_folder + self.label + '/'
            project_dir 		  = os.path.abspath(project_dir)
            network_dir 		  = os.path.abspath(project_dir + '/../../')
            py_file_common_header = ("import sys\nsys.path.append('{0}')\nsys.path.append('{1}')\nimport matplotlib"
                                     "\nmatplotlib.use('Agg')\nfrom modules.parameters import *\nfrom "
                                     "modules.analysis import *\nfrom computations import {2}\n\n") \
                .format(project_dir, network_dir, computation_function.__module__.split('.')[1])

            write_to_submit 	= []
            submit_jobs_file 	= main_experiment_folder + 'submit_jobs.py'
            job_list 			= main_experiment_folder + 'job_list.txt'
            cancel_jobs_file 	= main_experiment_folder + 'cancel_jobs.py'

            for par_set in self.parameter_sets:
                system2 			= par_set.kernel_pars.system
                template_file 		= system2['jdf_template']
                queueing_system     = system2['queueing_system']
                exec_file_name 		= remote_run_folder + par_set.label + '.sh'
                local_exec_file 	= main_experiment_folder + par_set.label + '.sh'
                computation_file 	= main_experiment_folder + par_set.label + '.py'
                remote_computation_file = remote_run_folder + par_set.label + '.py'

                par_set.save(main_experiment_folder + par_set.label+'.txt')
                with open(computation_file, 'w') as fp:
                    fp.write(py_file_common_header)
                    fp.write(computation_function.__module__.split('.')[1] + '.' + computation_function.__name__ +
                             "({0}, **{1})".format("'./" + par_set.label + ".txt'", str(parameters)))

                system2['jdf_fields'].update({'{{ computation_script }}': remote_computation_file,
                                              '{{ script_folder }}': remote_run_folder})

                io.process_template(template_file, system2['jdf_fields'].as_dict(), save_to=local_exec_file)
                write_to_submit.append("{0}\n".format(exec_file_name))

            with open(job_list, 'w') as fp:
                fp.writelines(write_to_submit)

            with open(submit_jobs_file, 'w') as fp:
                fp.write("import os\n")
                fp.write("import sys\n\n")
                fp.write("def submit_jobs(start_idx=0, stop_idx=None):\n")
                fp.write("\twith open('{0}') as fp:\n\t\tfor idx, line in enumerate(fp):\n".format('./job_list.txt'))
                fp.write("\t\t\tif stop_idx is not None:\n")
                fp.write("\t\t\t\tif (idx>=start_idx) and (idx<=stop_idx):\n")
                if queueing_system == 'slurm':
                    fp.write("\t\t\t\t\tos.system('sbatch {0}'.format(line))\n\t\t\telse:\n\t\t\t\tos.system('sbatch {0}'.format(line))")
                elif queueing_system == 'sge':
                    fp.write(
                        "\t\t\t\t\tos.system('qsub {0}'.format(line))\n\t\t\telse:\n\t\t\t\tos.system('qsub {"
                        "0}'.format(line))")
                # fp.write("\t\t\telse:\n")
                # fp.write("\t\t\t\tos.system('sbatch {0}'.format(line))")
                fp.write("\n\n")
                fp.write("if __name__=='__main__':\n\tif len(sys.argv)>1:\n\t\tsubmit_jobs(int(sys.argv[1]), "
                         "int(sys.argv[2]))")

            with open(cancel_jobs_file, 'w') as fp:
                if queueing_system == 'slurm':
                    fp.write(
                        "import os\nimport numpy as np\nimport sys\ndef cancel_range(init, end):\n\trang = np.arange("
                        "init, end)\n\tfor n in rang:\n\t\tos.system('scancel '+ str(n))\n\nif "
                        "__name__=='__main__':\n\tcancel_range(int(sys.argv[1]), int(sys.argv[2]))")
                elif queueing_system == 'sge':
                    fp.write(
                        "import os\nimport numpy as np\nimport sys\ndef cancel_range(init, end):\n\trang = np.arange("
                        "init, end)\n\tfor n in rang:\n\t\tos.system('qdel '+ str(n))\n\nif "
                        "__name__=='__main__':\n\tcancel_range(int(sys.argv[1]), int(sys.argv[2]))")

    def print_stored_keys(self, data_path):
        """
        Print all the nested keys in the results dictionaries for the current data set
        :param data_path: location of "Results" folder
        :return:
        """

        def pretty(d, indent=0):
            if isinstance(d, dict):
                for key, value in d.iteritems():
                    logger.info('  ' * indent + str(key))
                    if isinstance(value, dict):
                        pretty(value, indent + 1)

        pars_labels = [n.label for n in self]
        # open example data
        ctr = 0
        found_ = False
        while ctr < len(pars_labels) and not found_:
            try:
                with open(data_path + 'Results_' + pars_labels[ctr], 'r') as fp:
                    results = pickle.load(fp)
                # logger.info("Loading ParameterSet {0}".format(self.label))
                found_ = True
            except:
                logger.warning("Dataset {0} Not Found, skipping".format(pars_labels[ctr]))
                ctr += 1
                continue
        logger.info("\n\nResults dictionary structure:")
        pretty(results)

    def harvest(self, data_path, key_set=None, operation=None):
        """
        Gather stored results data and populate the space spanned by the Parameters with the corresponding results

        :param data_path: full path to global data folder
        :param key_set: specific result to extract from each results dictionary (nested keys should be specified as
        'key1/key2/...'
        :param operation: function to apply to result, if any
        :return: (parameter_labels, results_array)
        """
        if self.parameter_axes:
            l = ['x', 'y', 'z']
            domain_lens = [len(self.parameter_axes[l[idx]+'ticks']) for idx in range(self.dimensions)]
            domain_values = [self.parameter_axes[l[idx]+'ticks'] for idx in range(self.dimensions)]
            var_names = [self.parameter_axes[l[idx]+'label'] for idx in range(self.dimensions)]
            dom_len = np.prod(domain_lens)

            results_array = np.empty(tuple(domain_lens), dtype=object)
            parameters_array = np.empty(tuple(domain_lens), dtype=object)

            assert len(self) == dom_len, "Domain length inconsistent"
            pars_labels = [n.label for n in self]
            for n in range(int(dom_len)):
                params_label = self.label
                index = []
                if self.dimensions >= 1:
                    idx_x = n % domain_lens[0]
                    params_label += '_' + var_names[0] + '=' + str(domain_values[0][idx_x])
                    index.append(idx_x)
                if self.dimensions >= 2:
                    idx_y = (n // domain_lens[0]) % domain_lens[1]
                    params_label += '_' + var_names[1] + '=' + str(domain_values[1][idx_y])
                    index.append(idx_y)
                if self.dimensions == 3:
                    idx_z = ((n // domain_lens[0]) // domain_lens[1]) % domain_lens[2]
                    params_label += '_' + var_names[2] + '=' + str(domain_values[2][idx_z])
                    index.append(idx_z)

                parameters_array[tuple(index)] = pars_labels[pars_labels.index(params_label)]
                try:
                    with open(data_path+'Results_'+params_label, 'r') as fp:
                        results = pickle.load(fp)
                    logger.info("Loading ParameterSet {0}".format(params_label))
                except:
                    logger.warning("Dataset {0} not found, skipping".format(params_label))
                    continue
                if key_set is not None:
                    # print results
                    try:
                        nested_result = io.NestedDict(results)
                        if operation is not None:
                            results_array[tuple(index)] = operation(nested_result[key_set])
                        else:
                            # print key_set
                            results_array[tuple(index)] = nested_result[key_set]
                        # print nested_result[key_set]
                    except:
                        continue
                else:
                    results_array[tuple(index)] = results
        else:
            parameters_array = self.label
            results_array = []
            try:
                with open(data_path+'Results_'+self.label, 'r') as fp:
                    results = pickle.load(fp)
                logger.info("Loading Dataset {0}".format(self.label))
                if key_set is not None:
                    nested_result = io.NestedDict(results)
                    assert isinstance(results, dict), "Results must be dictionary"
                    if operation is not None:
                        results_array = operation(nested_result[key_set])
                    else:
                        results_array = nested_result[key_set]
                else:
                    results_array = results
            except IOError:
                logger.warning("Dataset {0} not found, skipping".format(self.label))

        return parameters_array, results_array def extract_nestvalid_dict(d, param_type='neuron'):
    """
    Verify whether the parameters dictionary are in the correct format, with adequate keys, in agreement with the nest
    parameters dictionaries so that they can later be passed as direct input to nest.

    :param d: parameter dictionary
    :param param_type: type of parameters - kernel, neuron, population, network, connections, topology
    :return: valid dictionary
    """
    if param_type == 'neuron' or param_type == 'synapse' or param_type == 'device':
        assert d['model'] in nest.Models(), "Model %s not currently implemented in %s" % (d['model'], nest.version())

        accepted_keys = nest.GetDefaults(d['model']).keys()
        accepted_keys.remove('model')
        nest_dict = {k: v for k, v in d.iteritems() if k in accepted_keys}
    elif param_type == 'kernel':
        accepted_keys = nest.GetKernelStatus().keys()
        nest_dict = {k: v for k, v in d.iteritems() if k in accepted_keys}
    else:
        # TODO
        logger.error("{!s} not implemented yet".format(param_type))
        exit(-1)

    return nest_dict class EncodingLayer:
    """
    Wrapper for all the encoders and generators involved in the input conversion process
    """
    # TODO create all instances using a prng; this means prng can't remain None... throw an error here?
    def __init__(self, initializer=None, signal=None, stim_seq=None, online=False, prng=None):
        """

        :param initializer: ParameterSet or dictionary with all parameters for Encoders and Generators
        :param signal: signal object to be converted
        :param stim_seq: [list of str] - stimulus sequence to be converted (only if input is a fixed spatiotemporal
        spike sequence)
        :param prng: random number generator for reproducibility (numpy.random)
        """
        self.prng = prng
        self.total_delay = 0.
        if initializer is not None:
            if isinstance(initializer, dict):
                initializer = parameters.ParameterSet(initializer)
            if not isinstance(initializer, parameters.ParameterSet):
                logger.info("Please provide the encoding parameters as dict or ParameterSet")
            if hasattr(initializer, 'encoder'):
                # if encoder exists, it is a population object
                self.n_encoders = initializer.encoder.N
                self.encoder_names = initializer.encoder.labels
            self.parameters = initializer
        else:
            self.n_encoders = 0
            self.encoder_names = []
            self.parameters = None

        self.n_generators = []
        if signal is not None and not online:
            if isinstance(signal, InputSignal):
                # generators convert the input signal
                self.n_generators.append(signal.dimensions)
            elif isinstance(signal, InputNoise):
                self.n_generators.append(signal.N)
            elif isinstance(signal, signals.AnalogSignalList) or isinstance(signal, signals.SpikeList):
                self.n_generators.append(len(signal))
            else:
                raise TypeError("signal must be InputSignal, AnalogSignalList or SpikeList")
        elif signal is not None and online:
            self.n_generators.append(signal.dimensions)
        elif stim_seq is not None:
            seq = list(signals.iterate_obj_list(stim_seq))
            assert isinstance(seq[0], str), "Provided stim_seq must be string sequence"
            self.n_generators.append(len(np.unique(seq)))
        else:
            self.n_generators.append(0)

        self.generator_names = []
        self.encoding_layer_gid = None
        self.generators = []
        self.encoders = []
        self.connections = []
        self.connection_types = []
        self.synaptic_weights = {}
        self.connection_delays = {}
        self.signal = signal

        def create_encoders(encoder_pars, signal=None):
            """
            Create the encoding layer
            :param encoder_pars: ParameterSet or dictionary with the encoder parameters
            :param signal: InputSignal object
            :return encoders, encoder_labels:
            """
            if isinstance(encoder_pars, dict):
                encoder_pars = parameters.ParameterSet(encoder_pars)
            assert isinstance(encoder_pars, parameters.ParameterSet), "Parameters must be provided as dictionary or ParameterSet"
            if signal is not None:
                assert (isinstance(signal, InputSignal) or isinstance(signal, InputNoise)), "Encoder Input must be InputSignal"

            encoder_labels = []
            encoders = []
            st_gen = StochasticGenerator(self.prng)
            for nn in range(encoder_pars.N):
                encoder_labels.append(encoder_pars.labels[nn])
                enc_pop_dict = {'pop_names': encoder_pars.labels[nn],
                                'n_neurons': encoder_pars.n_neurons[nn],
                                'topology': encoder_pars.topology[nn],
                                'topology_dict': encoder_pars.topology_dict[nn],
                                'gids': None,
                                'layer_gid': None,
                                'is_subpop': False}
                if hasattr(st_gen, encoder_pars.models[nn]):
                    enc_dict = parameters.copy_dict(enc_pop_dict, {})
                    encoders.append(Encoder(parameters.ParameterSet(enc_dict)))
                    logger.info("- {0} []".format(encoder_pars.models[nn]))

                elif encoder_pars.models[nn] == 'NEF':
                    neuron_dict = encoder_pars.neuron_pars[nn]
                    nest.CopyModel(neuron_dict['model'], encoder_pars.labels[nn])
                    nest.SetDefaults(encoder_pars.labels[nn], parameters.extract_nestvalid_dict(neuron_dict,
                                                                                                param_type='neuron'))

                    if encoder_pars.topology[nn]:
                        tp_dict = encoder_pars.topology_dict[nn]
                        tp_dict.update({'elements': encoder_pars.labels[nn]})
                        layer = tp.CreateLayer(tp_dict)
                        gids = nest.GetLeaves(layer)[0]
                        enc_dict = parameters.copy_dict(enc_pop_dict, {'gids': gids, 'layer_gid': layer})
                        encoders.append(Encoder(parameters.ParameterSet(enc_dict)))
                    else:
                        gids = nest.Create(encoder_pars.labels[nn], n=int(encoder_pars.n_neurons[nn]))
                        enc_dict = parameters.copy_dict(enc_pop_dict, {'gids': gids})
                        encoders.append(Encoder(parameters.ParameterSet(enc_dict)))

                    logger.info("- {0} Population of {1} neurons [{2}-{3}]".format(encoder_pars.models[nn],
                                                                             str(encoder_pars.n_neurons), str(min(gids)),
                                                                             str(max(gids))))
                    if encoder_pars.record_spikes[nn]:
                        if encoder_pars.spike_device_pars[nn].has_key('label'):
                            label = encoder_pars.spike_device_pars[nn]['label']
                        else:
                            label = 'encoder_spikes'
                        dev_gid = encoders[-1].record_spikes(parameters.extract_nestvalid_dict(
                            encoder_pars.spike_device_pars[nn], param_type='device'), label=label)
                        logger.info("- Connecting %s to %s, with label %s and id %s" % (
                            encoder_pars.spike_device_pars[nn]['model'],
                            encoders[-1].name, label, str(dev_gid)))

                    if encoder_pars.record_analogs[nn]:
                        if encoder_pars.analogue_device_pars[nn].has_key('label'):
                            label = encoder_pars.analogue_device_pars[nn]['label']
                        else:
                            label = 'encoder_analogs'
                        dev_gid = encoders[-1].record_analog(parameters.extract_nestvalid_dict(
                            encoder_pars.analogue_device_pars[nn], param_type='device'), label=label)
                        logger.info("- Connecting %s to %s, with label %s and id %s" % (
                            encoder_pars.analogue_device_pars[nn]['model'],
                            encoders[-1].name, label, str(dev_gid)))
                # TODO randomize initial variables...
                elif encoder_pars.models[nn] == 'parrot_neuron':
                    neuron_dict = encoder_pars.neuron_pars[nn]
                    nest.CopyModel(neuron_dict['model'], encoder_pars.labels[nn])
                    nest.SetDefaults(encoder_pars.labels[nn], parameters.extract_nestvalid_dict(neuron_dict,
                                                                                                param_type='neuron'))
                    if encoder_pars.topology[nn]:
                        tp_dict = encoder_pars.topology_dict[nn]
                        tp_dict.update({'elements': encoder_pars.labels[nn]})
                        layer = tp.CreateLayer(tp_dict)
                        gids = nest.GetLeaves(layer)[0]
                        enc_dict = parameters.copy_dict(enc_pop_dict, {'gids': gids, 'layer_gid': layer})
                        encoders.append(Encoder(parameters.ParameterSet(enc_dict)))
                    else:
                        gids = nest.Create(encoder_pars.labels[nn], n=int(encoder_pars.n_neurons[nn]))
                        enc_dict = parameters.copy_dict(enc_pop_dict, {'gids': gids})
                        encoders.append(Encoder(parameters.ParameterSet(enc_dict)))

                    logger.info("- {0} Population of {1} neurons [{2}-{3}]".format(encoder_pars.models[nn],
                                                                             str(encoder_pars.n_neurons[nn]),
                                                                             str(min(gids)),
                                                                             str(max(gids))))
                    if encoder_pars.record_spikes[nn]:
                        if encoder_pars.spike_device_pars[nn].has_key('label'):
                            label = encoder_pars.spike_device_pars[nn]['label']
                        else:
                            label = 'encoder_spikes'
                        dev_gid = encoders[-1].record_spikes(parameters.extract_nestvalid_dict(
                            encoder_pars.spike_device_pars[nn], param_type='device'), label=label)
                        logger.info("- Connecting %s to %s, with label %s and id %s" % (
                            encoder_pars.spike_device_pars[nn]['model'],
                            encoders[-1].name, label, str(dev_gid)))

                    if encoder_pars.record_analogs[nn]:
                        raise TypeError("Cannot record analogs from parrot neuron")
                else:
                    logger.warning("Not implemented yet!")

            return encoders, encoder_labels

        def create_generators(encoding_pars, signal_=None, input_dim=None):
            """
            Creates all necessary generator objects.
            :param encoding_pars: global encoding parameters
            :param signal_:
            :param input_dim:
            :return:
            """
            logger.info("\nCreating Generators: ")
            pars 			 = encoding_pars.generator
            generators 		 = []
            generator_labels = []

            for n in range(pars.N):
                # assess where the input to this generator comes from
                tmp = list(signals.iterate_obj_list(encoding_pars.connectivity.connections))
                targets = [x[0] for x in tmp]

                if pars.labels[n] in targets:
                    idx = targets.index(pars.labels[n])
                    source = tmp[idx][1]
                    if source in encoding_pars.encoder.labels:
                        src_idx = encoding_pars.encoder.labels.index(source)
                        tmp_input_dims = encoding_pars.encoder.n_neurons[src_idx]
                    else:
                        tmp_input_dims = None  # will be assessed from the signal
                # TODO list of possible generator types instead of label based
                elif pars.labels[n] == 'spike_pattern':
                    tmp_input_dims = encoding_pars.encoder.n_neurons[0]
                elif pars.labels[n] == 'X_noise' or 'noise' in pars.labels[n]:
                    tmp_input_dims = 1
                elif input_dim is not None:
                    tmp_input_dims = input_dim
                else:
                    tmp_input_dims = None

                # create a specific parameter set for each generator
                gen_pars_dict = {'label': pars.labels[n],
                                 'model': pars.models[n],
                                 'model_pars': pars.model_pars[n],
                                 'topology': pars.topology[n],
                                 'topology_pars': pars.topology_pars[n]}
                gen = Generator(gen_pars_dict, signal_, dims=tmp_input_dims)
                generators.append(gen)
                generator_labels.append(gen.name)
                logger.info("- {0} [{1}-{2}]".format(pars.labels[n], str(min(gen.gids)), str(max(gen.gids))))

            return generators, generator_labels

        ################################################################################
        if online:
            if hasattr(initializer, 'encoder') and initializer.encoder.N:
                self.encoders, self.encoder_names = create_encoders(initializer.encoder)
            if hasattr(initializer, 'generator') and signal is not None:
                self.generators, self.generator_names = create_generators(initializer,
                                                                          signal_=None, input_dim=signal.dimensions)
            elif hasattr(initializer, 'generator') and signal is None:
                self.generators, self.generator_names = create_generators(initializer)
        else:
            if hasattr(initializer, 'encoder'):
                self.encoders, self.encoder_names = create_encoders(initializer.encoder, signal)
            if hasattr(initializer, 'generator'):
                self.generators, self.generator_names = create_generators(initializer, signal)

    def _get_connection_type_and_parameters(self, encoding_pars, conn_pars, idx, populations, pop_objs):
        """
        Determine the type of the connecting populations and their parameters.
        :param encoding_pars:
        :param conn_pars:
        :param idx:
        :param populations:
        :param pop_objs:
        :return:
        """
        src_name = conn_pars.connections[idx][1]
        tget_name = conn_pars.connections[idx][0]

        if hasattr(encoding_pars, 'encoder') and (src_name in encoding_pars.encoder.labels):
            src_id = encoding_pars.encoder.labels.index(src_name)
            src_dims = self.encoders[src_id].size
            if conn_pars.topology_dependent[idx]:
                src_gids = self.encoders[src_id].layer_gid
            else:
                src_gids = self.encoders[src_id].gids
            src_tp = self.encoders[src_id].topology
        elif hasattr(encoding_pars, 'generator') and (src_name in encoding_pars.generator.labels):
            src_id = encoding_pars.generator.labels.index(src_name)
            src_dims = self.generators[src_id].input_dimension
            if conn_pars.topology_dependent[idx]:
                src_gids = self.generators[src_id].layer_gid
            else:
                src_gids = self.generators[src_id].gids
            src_tp = self.generators[src_id].topology
        else:
            raise TypeError("Source Population Label is not specified")

        if hasattr(encoding_pars, 'encoder') and (tget_name in encoding_pars.encoder.labels):
            tget_type = 'encoder'
            tget_id = encoding_pars.encoder.labels.index(tget_name)
            tget_dims = self.encoders[tget_id].size
            if conn_pars.topology_dependent[idx]:
                tget_gids = self.encoders[tget_id].layer_gid
            else:
                tget_gids = self.encoders[tget_id].gids
            tget_tp = self.encoders[tget_id].topology
        elif hasattr(encoding_pars, 'generator') and (tget_name in encoding_pars.generator.labels):
            tget_type = 'generator'
            tget_id = encoding_pars.generator.labels.index(tget_name)
            tget_dims = self.generators[tget_id].input_dimension
            if conn_pars.topology_dependent[idx]:
                tget_gids = self.generators[tget_id].layer_gid
            else:
                tget_gids = self.generators[tget_id].gids
            tget_tp = self.generators[tget_id].topology
        elif tget_name in populations:
            tget_type = 'population'
            tget_id = populations.index(tget_name)
            tget_dims = pop_objs[tget_id].size
            if conn_pars.topology_dependent[idx]:
                tget_gids = pop_objs[tget_id].layer_gid
            else:
                tget_gids = pop_objs[tget_id].gids
            tget_tp = pop_objs[tget_id].topology
        else:
            raise TypeError("Target Population Label is not specified or is not among the "
                            "populations: {0}".format(tget_name))

        return src_gids, src_id, src_dims, src_tp, tget_gids, tget_id, tget_dims, tget_tp

    def connect(self, encoding_pars, net_obj=None, progress=True):
        """
        Connect the generators to their target populations according to specifications
        :param encoding_pars: ParameterSet or dictionary
        :param net_obj: Network object containing the populations to connect to
        """
        logger.info("\nConnecting Encoding Layer: ")

        if isinstance(encoding_pars, dict):
            encoding_pars = parameters.ParameterSet(encoding_pars)
        assert isinstance(encoding_pars, parameters.ParameterSet), "Parameters must be dictionary or ParameterSet"
        if net_obj is not None:
            assert isinstance(net_obj, net_architect.Network), "Please provide the network object to connect to"
        else:
            logger.info("No network object provided, cannot connect input to network")

        conn_pars 	= encoding_pars.connectivity
        populations = list(signals.iterate_obj_list(net_obj.population_names))
        pop_objs 	= list(signals.iterate_obj_list(net_obj.populations))

        if not signals.empty(net_obj.merged_populations):
            merg_pop_names = [x.name for x in net_obj.merged_populations]
            merg_pops = [x for x in net_obj.merged_populations]
            populations.extend(merg_pop_names)
            pop_objs.extend(merg_pops)

        # iterate through all connections
        for idx, nn in enumerate(conn_pars.connections):
            src_name = nn[1]
            tget_name = nn[0]
            # determine the type of the connecting populations and their parameters
            try:
                src_gids, src_id, src_dims, src_tp, tget_gids, tget_id, tget_dims, tget_tp = \
                    self._get_connection_type_and_parameters(encoding_pars, conn_pars, idx, populations, pop_objs)
            except NameError:  # we just continue, nothing special or wrong here
                continue

            if "synapse_name" in conn_pars and conn_pars.synapse_name[idx] is not None:
                synapse_name = conn_pars.synapse_name[idx]
            else:
                synapse_name = src_name + '_' + tget_name

            self.connection_types.append(synapse_name)
            self.connections.append((tget_name, src_name))

            if conn_pars.models[idx] is not None and synapse_name not in nest.Models():
                nest.CopyModel(conn_pars.models[idx], synapse_name)
                nest.SetDefaults(synapse_name, conn_pars.model_pars[idx])

            if synapse_name.find('copy') > 0:  # re-connect the same neurons...
                start = time.time()
                logger.info("    - Connecting {0} (*)".format(synapse_name))

                device_models = ['spike_detector', 'spike_generator', 'multimeter']
                target_synapse_name = synapse_name[:synapse_name.find('copy')-1]
                conns = nest.GetConnections(synapse_model=target_synapse_name)

                iterate_steps = 100
                its = np.arange(0, len(conns)+1, iterate_steps).astype(int)

                for nnn, it in enumerate(its):
                    if nnn < len(its)-1:
                        con = conns[it:its[nnn+1]]
                        st = nest.GetStatus(con)
                        source_gids = [x['source'] for x in st if
                                       nest.GetStatus([x['target']])[0]['model'] not in device_models]
                        target_gids = [x['target'] for x in st if
                                       nest.GetStatus([x['target']])[0]['model'] not in device_models]
                        syn_dict = parameters.copy_dict(conn_pars.syn_specs[idx], {'model': synapse_name, 'weight':
                            conn_pars.weight_dist[idx], 'delay': conn_pars.delay_dist[idx]})

                        syn_dicts = [{'synapsemodel': list(np.repeat(synapse_name, len(source_gids)))[iddx],
                                      'source': source_gids[iddx],
                                      'target': target_gids[iddx],
                                      'weight': syn_dict['weight'],		              # TODO distributions??
                                      'delay': syn_dict['delay'],
                                      'receptor_type': syn_dict['receptor_type']} for iddx in range(len(target_gids))]
                        nest.DataConnect(syn_dicts)
                    if progress:
                        visualization.progress_bar(float(nnn) / float(len(its)))
                logger.info("\tElapsed time: {0} s".format(str(time.time()-start)))
            else:
                ## Set up connections
                # 1) pre-computed weight matrix, or function to specify w
                if conn_pars.preset_W[idx] is not None:
                    if isinstance(conn_pars.preset_W[idx], np.ndarray):
                        w = conn_pars.preset_W[idx]
                    elif isinstance(conn_pars.preset_W[idx], tuple):
                        args = parameters.copy_dict(conn_pars.preset_W[idx][1], {'size': (int(src_dims), int(tget_dims))})
                        w = conn_pars.preset_W[idx][0](**args)
                    elif isinstance(conn_pars.preset_W[idx], str):
                        w = np.load(conn_pars.preset_W[idx])
                    else:
                        raise ValueError("Incorrect W specifications")

                    logger.info("- Connecting {0} to population {1}".format(src_name, tget_name))

                    for pre_syn_matidx, pre_syn_gid in enumerate(src_gids):
                        post_syn_matidx = w[pre_syn_matidx, :].nonzero()[0]
                        post_syn_gid = list(post_syn_matidx + min(tget_gids))
                        weights = [w[pre_syn_matidx, x] for x in post_syn_matidx]

                        if isinstance(conn_pars.delay_dist[idx], dict):
                            delay = conn_pars.delay_dist[idx]
                        elif isinstance(conn_pars.delay_dist[idx], float) or \
                                isinstance(conn_pars.delay_dist[idx], int):
                            delay = np.repeat(conn_pars.delay_dist[idx], len(weights))
                        elif isinstance(conn_pars.delay_dist[idx], tuple):
                            args = parameters.copy_dict(conn_pars.delay_dist[idx][1], {'size': len(weights)})
                            delay = conn_pars.delay_dist[idx][0](**args)
                        else:
                            raise TypeError("Delays not provided in correct format")
                        visualization.progress_bar(float(pre_syn_matidx) / len(src_gids))
                        conn_dict = conn_pars.conn_specs[idx]
                        for idxx, tget in enumerate(post_syn_gid):
                            if isinstance(delay, dict):
                                d = delay
                            else:
                                d = delay[idxx]

                            syn_dict = parameters.copy_dict(conn_pars.syn_specs[idx], {'model': synapse_name,
                                                                                       'weight': weights[idxx],
                                                                                       'delay': d})

                            if conn_dict is not None:
                                nest.Connect(pre_syn_gid, [tget], conn_spec=conn_dict, syn_spec=syn_dict)
                            else:
                                nest.Connect(pre_syn_gid, [tget], 'all_to_all', syn_dict)
                    logger.info("- Connecting {0} to population {1}".format(src_name, tget_name))

                # topology connections only allowed between generator (with specified topology) and population
                elif conn_pars.topology_dependent[idx]:
                    assert (src_tp and tget_tp), "Topological connections are only possible if both pre- and post- " \
                                                 "population are topologically defined"

                    tp_dict = parameters.copy_dict(conn_pars.conn_specs[idx], {'synapse_model': synapse_name})
                    tp.ConnectLayers(src_gids, tget_gids, tp_dict)

                    logger.info("- Connecting Layer {0} to Layer {1}".format(src_name, tget_name))

                else:
                    # add extra numerical index to the src_name
                    if ((src_name + '0' in nest.Models()) and (tget_name in nest.Models())) \
                            or ((src_name in nest.Models()) and (tget_name in nest.Models())) \
                            or ((src_name + '0' in nest.Models()) and
                                (nest.GetStatus([tget_gids[0]])[0]['element_type'] == 'neuron')) \
                            or ((src_name in nest.Models()) and
                                (nest.GetStatus([tget_gids[0]])[0]['element_type'] == 'neuron')):

                        # need to use synapse with the same model for all connections from a device
                        if "synapse_name" in conn_pars and conn_pars.synapse_name[idx] is not None:
                            synapse_name = conn_pars.synapse_name[idx]
                        else:
                            synapse_name = src_name + '_' + tget_name
                        nest.SetDefaults(conn_pars.models[idx], conn_pars.model_pars[idx])
                        syn_dict = parameters.copy_dict(conn_pars.syn_specs[idx], {'model': synapse_name, 'weight':
                            conn_pars.weight_dist[idx], 'delay': conn_pars.delay_dist[idx]})
                        conn_dict = conn_pars.conn_specs[idx]
                        if isinstance(src_gids, list):
                            src_gids = tuple(x[0] for x in src_gids)
                        nest.Connect(src_gids, tget_gids, conn_spec=conn_dict, syn_spec=syn_dict)
                        logger.info("- Connecting {0} to population {1} [{2}]".format(src_name, tget_name, synapse_name))
                    else:
                        x = StochasticGenerator(self.prng)
                        if hasattr(x, encoding_pars.encoder.models[src_id]):
                            # src is a stochastic encoder... set spike times for spike generator
                            act_list = self.encoders[src_id].spiking_activity
                            if isinstance(tget_gids, list):
                                tget_gids = tuple(x[0] for x in tget_gids)

                            dt = nest.GetKernelStatus()['resolution']
                            d = decimal.Decimal(str(dt)).as_tuple()[1][0]
                            for idxx, nnn in enumerate(act_list.id_list):
                                spk_times = act_list.spiketrains[nnn].spike_times
                                spk_times[np.where(spk_times == 0.)] = dt
                                spk_t = [round(tt, d) for tt in spk_times]
                                nest.SetStatus([tget_gids[idxx]], {'spike_times': spk_t})

                            logger.info("- Connecting {0} to generator {1}".format(src_name, tget_name))

                        elif (src_name in nest.Models()) and (tget_name in nest.Models()):
                            if "synapse_name" in conn_pars and conn_pars.synapse_name[idx] is not None:
                                synapse_name = conn_pars.synapse_name[idx]
                            else:
                                synapse_name = src_name + '_' + tget_name
                            nest.SetDefaults(conn_pars.models[idx], conn_pars.model_pars[idx])
                            syn_dict = parameters.copy_dict(conn_pars.syn_specs[idx], {'model': synapse_name, 'weight':
                                conn_pars.weight_dist[idx], 'delay': conn_pars.delay_dist[idx]})
                            conn_dict = conn_pars.conn_specs[idx]
                            if isinstance(src_gids, list):
                                src_gids = tuple(x[0] for x in src_gids)
                            nest.Connect(src_gids, tget_gids, conn_spec=conn_dict, syn_spec=syn_dict)
                            logger.info("- Connecting {0} to population {1} [{2}]".format(src_name, tget_name, synapse_name))
        self.signal = None

    # def connect_clone(self, encoding_pars, network=None, clone=None): TODO safely remove!
    # 	"""
    # 	Replicate EncodingLayer so that both network and clone network receive the exact same input
    # 	:param encoding_pars:
    # 	:param network:
    # 	:param clone:
    # 	:return:
    # 	"""
    # 	native_population_names = list(signals.iterate_obj_list(network.population_names))
    # 	native_populations = list(signals.iterate_obj_list(network.populations))
    # 	target_populations = [n[0] for n in encoding_pars.connectivity.connections if n[0] in native_population_names]
    # 	target_populations.extend([n+'_clone' for n in target_populations])
    #
    # 	if hasattr(encoding_pars, "encoder") and encoding_pars.encoder.N:
    # 		encoding_pars.encoder.labels.extend(['{0}_parrots'.format(n) for n in
    # 		                                                    native_population_names if n in target_populations])
    # 		new_encoders = encoding_pars.encoder.labels
    # 		encoder_size = encoding_pars.encoder.n_neurons.extend([n.size for idx, n in enumerate(
    # 				native_populations) if native_population_names[idx] in target_populations])
    # 		encoder_pars = encoding_pars.encoder.neuron_pars.extend([{'model': 'parrot_neuron'} for n in native_population_names if n in target_populations])
    # 		enc_models = encoding_pars.encoder.models.extend(['parrots' for _ in range(len(new_encoders))])
    # 		enc_model_pars = encoding_pars.encoder.model_pars.extend([{} for _ in range(len(new_encoders))])
    # 		topologies = encoding_pars.encoder.topology.extend([False for _ in range(len(new_encoders))])
    # 		tp_dicts = encoding_pars.encoder.topology_dict.extend([None for _ in range(len(new_encoders))])
    # 		rc_spikes = encoding_pars.encoder.record_spikes.extend([False for _ in range(len(new_encoders))])
    # 		spk_dvc_pars = encoding_pars.encoder.spike_device_pars.extend([None for _ in range(len(new_encoders))])
    # 		rc_analogs = encoding_pars.encoder.record_analogs.extend([False for _ in range(len(new_encoders))])
    # 		analog_dev_pars = encoding_pars.encoder.analog_device_pars.extend([None for _ in range(len(new_encoders))])
    # 	else:
    # 		new_encoders = ['{0}_parrots'.format(n) for n in native_population_names if n in target_populations]
    # 		encoder_size = [n.size for n in native_populations]
    # 		encoder_pars = [{'model': 'parrot_neuron'} for n in native_population_names if n in target_populations]
    # 		enc_models = ['NEF' for _ in range(len(new_encoders))]
    # 		enc_model_pars = [{} for _ in range(len(new_encoders))]
    # 		topologies = [False for _ in range(len(new_encoders))]
    # 		tp_dicts = [None for _ in range(len(new_encoders))]
    # 		rc_spikes = [False for _ in range(len(new_encoders))]
    # 		spk_dvc_pars = [None for _ in range(len(new_encoders))]
    # 		rc_analogs = [False for _ in range(len(new_encoders))]
    # 		analog_dev_pars = [None for _ in range(len(new_encoders))]
    #
    # 	connections = []
    # 	connections_clone = []
    # 	conn_specs = []
    # 	conn_specs_clone = []
    # 	syn_names = []
    # 	syn_names_clone = []
    # 	syn_specs = []
    # 	syn_specs_clone = []
    # 	models = []
    # 	models_clone = []
    # 	model_pars = []
    # 	model_pars_clone = []
    # 	weights = []
    # 	weights_clone = []
    # 	delays = []
    # 	delays_clone = []
    # 	pre_w = []
    # 	pre_w_clone = []
    # 	for idx, n_connection in enumerate(encoding_pars.connectivity.connections):
    # 		if n_connection[0] in native_population_names:
    # 			connections.append(('{0}'.format(str(n_connection[0]+'_parrots')), '{0}'.format(str(n_connection[1]))))
    # 			#connections_clone.append(('{0}'.format(str(n_connection[0]+'_parrots')), '{0}'.format(str(
    # 			# n_connection[1]))))
    # 			conn_specs.append({'rule': 'all_to_all'})
    # 			syn_names.append('{0}Parrots'.format(str(n_connection[1])))
    # 			syn_specs.append({})
    # 			models.append('static_synapse')
    # 			model_pars.append({})
    # 			weights.append(1.)
    # 			delays.append(0.1)
    # 			pre_w.append(None)
    #
    # 			connections.append(('{0}'.format(str(n_connection[0])), '{0}'.format(str(n_connection[0]+'_parrots'))))
    # 			connections_clone.append(('{0}_clone'.format(str(n_connection[0])), '{0}'.format(str(n_connection[0]+'_parrots'))))
    # 			conn_specs.append({'rule': 'one_to_one'})
    # 			conn_specs_clone.append({'rule': 'one_to_one'})
    # 			syn_names.append(encoding_pars.connectivity.synapse_name[idx])
    # 			syn_names_clone.append(encoding_pars.connectivity.synapse_name[idx])
    # 			syn_specs.append(encoding_pars.connectivity.syn_specs[idx])
    # 			syn_specs_clone.append(encoding_pars.connectivity.syn_specs[idx])
    # 			models.append(encoding_pars.connectivity.models[idx])
    # 			models_clone.append(encoding_pars.connectivity.models[idx])
    # 			model_pars.append(encoding_pars.connectivity.model_pars[idx])
    # 			model_pars_clone.append(encoding_pars.connectivity.model_pars[idx])
    # 			weights.append(encoding_pars.connectivity.weight_dist[idx])
    # 			weights_clone.append(encoding_pars.connectivity.weight_dist[idx])
    # 			delays.append(encoding_pars.connectivity.delay_dist[idx])
    # 			delays_clone.append(encoding_pars.connectivity.delay_dist[idx])
    # 			pre_w.append(encoding_pars.connectivity.preset_W[idx])
    # 			pre_w_clone.append(encoding_pars.connectivity.preset_W[idx])
    # 		else:
    # 			connections.append(n_connection)
    # 			#connections_clone.append(n_connection)
    # 			conn_specs.append(encoding_pars.connectivity.conn_specs[idx])
    # 			syn_names.append(encoding_pars.connectivity.synapse_name[idx])
    # 			syn_specs.append(encoding_pars.connectivity.syn_specs[idx])
    # 			models.append(encoding_pars.connectivity.models[idx])
    # 			model_pars.append(encoding_pars.connectivity.model_pars[idx])
    # 			weights.append(encoding_pars.connectivity.weight_dist[idx])
    # 			delays.append(encoding_pars.connectivity.delay_dist[idx])
    # 			pre_w.append(encoding_pars.connectivity.preset_W[idx])
    #
    # 	tp = [False for _ in range(len(connections))]
    # 	extra_encoder = {
    # 		'N': len(new_encoders),
    # 		'labels': new_encoders,
    # 		'models': enc_models,
    # 		'model_pars': enc_model_pars,
    # 		'n_neurons': encoder_size,
    # 		'neuron_pars': encoder_pars,
    # 		'topology': topologies,
    # 		'topology_dict': tp_dicts,
    # 		'record_spikes': rc_spikes,
    # 		'spike_device_pars': spk_dvc_pars,
    # 		'record_analogs': rc_analogs,
    # 		'analog_device_pars': analog_dev_pars
    # 	}
    # 	extra_connectivity_native = {
    # 			'synapse_name': syn_names,
    # 			'connections': connections,
    # 			'topology_dependent': tp,
    # 			'conn_specs': conn_specs,
    # 			'syn_specs': syn_specs,
    # 			'models': models,
    # 			'model_pars': model_pars,
    # 			'weight_dist': weights,
    # 			'delay_dist': delays,
    # 			'preset_W': pre_w}
    # 	extra_connectivity_clone = {
    # 			'synapse_name': syn_names_clone,
    # 			'connections': connections_clone,
    # 			'topology_dependent': [False for _ in range(len(connections_clone))],
    # 			'conn_specs': conn_specs_clone,
    # 			'syn_specs': syn_specs_clone,
    # 			'models': models_clone,
    # 			'model_pars': model_pars_clone,
    # 			'weight_dist': weights_clone,
    # 			'delay_dist': delays_clone,
    # 			'preset_W': pre_w_clone}
    # 	native_encoding_pars = parameters.ParameterSet({'encoder': extra_encoder, 'generator': encoding_pars.generator.as_dict(),
    # 	                                     'connectivity': extra_connectivity_native})
    # 	clone_encoding_pars = parameters.ParameterSet({'encoder': extra_encoder, 'generator': encoding_pars.generator.as_dict(),
    # 	                                     'connectivity': extra_connectivity_clone})
    # 	self.__init__(native_encoding_pars)
    # 	self.connect(native_encoding_pars, network)
    # 	self.connect(clone_encoding_pars, clone)

    @staticmethod
    def replicate_connections(net, clone, progress=True): # TODO - correct or remove
        """
        Replicate the connectivity from the encoding layer to the clone network
        :param net:
        :param clone:
        :param progress:
        :return:
        """
        start = time.time()
        target_population_names = clone.population_names
        target_population_gids = [n.gids for n in clone.populations]
        source_population_names = [n.split('_')[0] for n in target_population_names]
        source_population_gids = [net.populations[net.population_names.index(n)].gids for n in source_population_names]
        device_models = ['spike_detector', 'multimeter'] # 'spike_generator',

        for idx, n_pop in enumerate(source_population_gids):
            logger.info("\n Replicating Encoding Layer connections to {0}".format(clone.population_names[idx]))
            conns 			= nest.GetConnections(target=n_pop)
            iterate_steps 	= 100
            its 			= np.arange(0, len(conns) + 1, iterate_steps).astype(int)
            if len(its) > 1:
                for nnn, it in enumerate(its):
                    if nnn < len(its) - 1:
                        con = conns[it:its[nnn + 1]]
                        st = nest.GetStatus(con)

                        source_gids = [x['source'] for x in st if nest.GetStatus([x['source']])[0]['model'] not in
                                       device_models]

                        target_gids = [target_population_gids[idx][0] for x in st
                                       if nest.GetStatus([x['source']])[0]['model'] not in device_models]

                        weights 	= [x['weight'] for x in st if nest.GetStatus([x['source']])[0]['model'] not in device_models]
                        delays 		= [x['delay'] for x in st if nest.GetStatus([x['source']])[0]['model'] not in device_models]
                        models 		= [str(x['synapse_model']) for x in st if nest.GetStatus([x['source']])[0]['model'] not in device_models]
                        receptors 	= [x['receptor'] for x in st if nest.GetStatus([x['source']])[0]['model'] not in device_models]
                        syn_dicts 	= [{'synapsemodel': models[iddx],
                                         'source': source_gids[iddx],
                                         'target': target_gids[iddx],
                                         'weight': weights[iddx],
                                         'delay': delays[iddx],
                                         'receptor_type': receptors[iddx]} for iddx in range(len(target_gids))]
                        nest.DataConnect(syn_dicts)
                    if progress:
                        visualization.progress_bar(float(nnn) / float(len(its)))

                logger.info("\tElapsed time: {0} s".format(str(time.time() - start)))
            else:
                st = nest.GetStatus(conns)
                logger.info([nest.GetStatus([x['source']])[0]['model'] for x in st])

                source_gids = [x['source'] for x in st if nest.GetStatus([x['source']])[0]['model'] not in
                               device_models]
                target_gids = [target_population_gids[idx][0] for x in st if nest.GetStatus([x['source']])[0][
                    'model'] not in device_models]
                weights 	= [x['weight'] for x in st if nest.GetStatus([x['source']])[0]['model'] not in device_models]
                delays 		= [x['delay'] for x in st if nest.GetStatus([x['source']])[0]['model'] not in device_models]
                models 		= [str(x['synapse_model']) for x in st if
                                 nest.GetStatus([x['source']])[0]['model'] not in device_models]
                receptors = [x['receptor'] for x in st if
                             nest.GetStatus([x['source']])[0]['model'] not in device_models]

                syn_dicts = [{'synapsemodel': models[iddx],
                              'source': source_gids[iddx],
                              'target': target_gids[iddx],
                              'weight': weights[iddx],
                              'delay': delays[iddx],
                              'receptor_type': receptors[iddx]} for iddx in range(len(target_gids))]
                nest.DataConnect(syn_dicts)
                logger.info("\tElapsed time: {0} s".format(str(time.time() - start)))

    def connect_decoders(self, decoding_pars):
        """
        Connect decoders to encoding population
        :param decoding_pars:
        :return:
        """
        if isinstance(decoding_pars, dict):
            decoding_pars = parameters.ParameterSet(decoding_pars)
        assert isinstance(decoding_pars, parameters.ParameterSet), "DecodingLayer must be initialized with ParameterSet or " \
                                                                   "dictionary"
        decoder_params = {}
        # initialize state extractors:
        if hasattr(decoding_pars, "state_extractor"):
            pars_st = decoding_pars.state_extractor
            logger.info("\nConnecting Decoders: ")
            for ext_idx, n_src in enumerate(pars_st.source_population):
                assert(n_src in self.encoder_names), "State extractor must connect to encoder population"
                if n_src in self.encoder_names:
                    pop_index = self.encoder_names.index(n_src)
                    src_obj = self.encoders[pop_index]
                else:
                    raise TypeError("No source populations in Encoding Layer")

                decoder_params.update({'state_variable': pars_st.state_variable,
                                       'state_specs': pars_st.state_specs,
                                       'reset_states': pars_st.reset_states,
                                       'average_states': pars_st.average_states,
                                       'standardize': pars_st.standardize})

                if hasattr(decoding_pars, "readout"):
                    pars_readout = decoding_pars.readout
                    assert(len(decoding_pars.readout) == decoding_pars.state_extractor.N), "Specify one readout dictionary " \
                                                                                           "per state extractor"
                    decoder_params.update({'readout': pars_readout})
                src_obj.connect_decoders(parameters.ParameterSet(decoder_params))
        else:
            raise IOError("DecodingLayer requires the specification of state extractors")

    def extract_synaptic_weights(self, src_gids=None, tget_gids=None, syn_name=None, progress=True):
        """
        Determine the connection weights between src_gids and tget_gids for synapses syn_name
        :param src_gids:
        :param tget_gids:
        :param syn_name:
        :param progress:
        :return:
        """
        # TODO remove the use of 'copy' keyword (this is now done with the multiport synapse in NEST)
        if src_gids is None and tget_gids is None:
            for con in list(np.unique(self.connection_types)):
                if con[-4:] != 'copy':
                    status_dict = nest.GetStatus(nest.GetConnections(synapse_model=con))
                    src_gids = [status_dict[n]['source'] for n in range(len(status_dict))]
                    tget_gids = [status_dict[n]['target'] for n in range(len(status_dict))]
                    self.synaptic_weights.update({con: net_architect.extract_weights_matrix(list(np.unique(src_gids)),
                                                                                            list(np.unique(tget_gids)),
                                                                                            progress=progress)})

        elif src_gids and tget_gids:
            if syn_name is None:
                syn_name = str(nest.GetStatus(nest.GetConnections([src_gids[0]], [tget_gids[0]]))[0]['synapse_model'])
            self.synaptic_weights.update({syn_name: net_architect.extract_weights_matrix(list(np.unique(src_gids)),
                                                                                         list(np.unique(tget_gids)),
                                                                                         progress=progress)})
        else:
            logger.warning("Provide gids!!")

    def extract_synaptic_delays(self, src_gids=None, tget_gids=None, syn_name=None, progress=True):
        """
        Determine the synaptic delays between src_gids and tget_gids for synapses syn_name
        :param src_gids: source NEST ids
        :param tget_gids: target NEST ids
        :param syn_name: name of synapse
        :param progress: show progress bar
        :return:
        """
        if src_gids is None:
            for con in list(np.unique(self.connection_types)):
                status_dict = nest.GetStatus(nest.GetConnections(synapse_model=con))
                src_gids = [status_dict[n]['source'] for n in range(len(status_dict))]
                tget_gids = [status_dict[n]['target'] for n in range(len(status_dict))]
                self.connection_delays.update({con: net_architect.extract_delays_matrix(list(np.unique(src_gids)),
                                                                                        list(np.unique(tget_gids)),
                                                                                        progress=progress)})
        elif src_gids and tget_gids:
            if syn_name is None:
                syn_name = str(nest.GetStatus(nest.GetConnections([src_gids[0]], [tget_gids[0]]))[0]['synapse_model'])
            self.connection_delays.update({syn_name: net_architect.extract_delays_matrix(list(np.unique(src_gids)),
                                                                                         list(np.unique(tget_gids)),
                                                                                         progress=progress)})
        else:
            logger.warning("Provide gids!!")

    def extract_encoder_activity(self, t_start=None, t_stop=None):
        """
        Read the spiking activity of the encoders present in the EncodingLayer
        :return:
        """
        if not signals.empty(self.encoders):
            logger.info("\nExtracting and storing recorded activity from encoders: ")
        for n_enc in self.encoders:
            if n_enc.attached_devices:
                logger.info("- Encoder {0}".format(n_enc.name))
                for n_dev in n_enc.attached_devices:
                    if nest.GetStatus(n_dev)[0]['to_memory']:
                        # status = nest.GetStatus(n_dev)[0]['events']
                        # t_start, t_stop = np.min(status['times']), np.max(status['times'])
                        n_enc.load_activity(n_dev, t_start=t_start, t_stop=t_stop)
                    elif nest.GetStatus(n_dev)[0]['to_file']:
                        n_enc.load_activity(list(nest.GetStatus(n_enc)[0]['filenames']), t_start=t_start, t_stop=t_stop)

    def flush_records(self, decoders=False):
        """
        Delete all data from all devices connected to the encoders
        :return:
        """
        if not signals.empty(self.encoders):
            logger.info("\nClearing device data: ")
        for n_enc in self.encoders:
            if n_enc.attached_devices:
                devices = n_enc.attached_device_names
                for idx, n_dev in enumerate(n_enc.attached_devices):
                    logger.info(" - {0} {1}".format(devices[idx], str(n_dev)))
                    nest.SetStatus(n_dev, {'n_events': 0})
                    if nest.GetStatus(n_dev)[0]['to_file']:
                        io.remove_files(nest.GetStatus(n_dev)[0]['filenames'])
        if decoders:
            for n_enc in self.encoders:
                if n_enc.decoding_layer is not None:
                    n_enc.decoding_layer.flush_records()

    def update_state(self, signal, prev_signal_ids=None, ids_to_update=None):
        """
        Update state of encoders and generators. If ids_to_update is given, then only the appropriate generators
        will be updated for efficiency reasons.

        :param signal: [AnalogSignalList] signal used to update state of encoders and generators
        :param prev_signal_ids: [list] with signal ids that were active/updated in the previous step
        :param ids_to_update: [list] with signal ids for which the generators should be updated
        :return:
        """
        for n_gen in self.generators:
            n_gen.update_state(signal, prev_signal_ids, ids_to_update)

    # TODO this could use a bit of rethinking
    def extract_connectivity(self, net, sub_set=False, progress=False):
        """
        Extract encoding layer connections.

        :param net: Network object
        :param sub_set:
        :param progress: display progress bar
        :return:
        """
        n_items = 10
        if sub_set:
            tgets = list(signals.iterate_obj_list([list(signals.iterate_obj_list(n.gids)) for n in
                                                   net.populations]))[:n_items]
            srces = list(itertools.chain(*[n.gids for n in self.generators][0]))[:n_items]
            self.extract_synaptic_weights(srces, tgets, syn_name='Gen_Net', progress=progress)
            self.extract_synaptic_delays(srces, tgets, syn_name='Gen_Net', progress=progress)

        # if there's only one connection type
        elif len(np.unique(self.connection_types)) == 1:
            tgets = list(signals.iterate_obj_list([list(signals.iterate_obj_list(n.gids)) for n in net.populations]))
            srces = list(itertools.chain(*[n.gids for n in self.generators][0]))
            self.extract_synaptic_weights(srces, tgets, syn_name='Gen_Net', progress=progress)
            self.extract_synaptic_delays(srces, tgets, syn_name='Gen_Net', progress=progress)

        elif len(self.connections) != len(np.unique(self.connection_types)):
            for con_idx, n_con in enumerate(self.connections):
                src_name = n_con[1]
                tget_name = n_con[0]
                syn_name = self.connection_types[con_idx] + str(con_idx)
                if src_name in net.population_names:
                    src_gids = net.populations[net.population_names.index(src_name)].gids
                elif src_name in list([n.name for n in net.merged_populations]):
                    merged_populations = list([n.name for n in net.merged_populations])
                    src_gids = net.merged_populations[merged_populations.index(src_name)].gids
                elif src_name in self.encoder_names:
                    src_gids = self.encoders[self.encoder_names.index(src_name)].gids
                else:
                    gen_names = [x.name[0].split('0')[0] for x in self.generators]
                    src_gids = list(itertools.chain(*self.generators[gen_names.index(src_name)].gids))
                if tget_name in net.population_names:
                    tget_gids = net.populations[net.population_names.index(tget_name)].gids
                elif tget_name in list([n.name for n in net.merged_populations]):
                    merged_populations = list([n.name for n in net.merged_populations])
                    tget_gids = net.merged_populations[merged_populations.index(tget_name)].gids
                elif tget_name in self.encoder_names:
                    tget_gids = self.encoders[self.encoder_names.index(tget_name)].gids
                else:
                    gen_names = [x.name[0].split('0')[0] for x in self.generators]
                    tget_gids = list(itertools.chain(*self.generators[gen_names.index(tget_name)].gids))
                self.extract_synaptic_weights(src_gids, tget_gids, syn_name=syn_name, progress=progress)
                self.extract_synaptic_delays(src_gids, tget_gids, syn_name=syn_name, progress=progress)
        else:
            self.extract_synaptic_weights(progress=progress)
            self.extract_synaptic_delays(progress=progress)

    def determine_total_delay(self):
        """
        Determine the connection delays involved in the encoding layer
        :return:
        """
        assert (not signals.empty(self.connection_delays)), "Please run extract_connectivity first..."

        for k, v in self.connection_delays.items():
            delay = np.unique(np.array(v[v.nonzero()].todense()))
            assert (len(delay) == 1), "Heterogeneous delays in encoding layer are not supported.."

        self.total_delay = float(delay)
        logger.info("\n- total delays in EncodingLayer: {0} ms".format(str(self.total_delay))) class InputNoise(StochasticGenerator):
    """
    Generate and store AnalogSignal object referring to the noise to add to the input signal u(t)
    """
    def __init__(self, initializer, rng=None, seed=None, start_time=None, stop_time=None):
        """

        :param initializer:
        :param rng:
        :param seed:
        :param stop_time:
        """
        StochasticGenerator.__init__(self, rng, seed)

        if isinstance(initializer, dict):
            initializer = parameters.ParameterSet(initializer)

        assert isinstance(initializer, parameters.ParameterSet), "Initializer must be a parameter dictionary or ParameterSet"
        self.N = initializer.N
        self.rectify = initializer.rectify
        self.source = initializer.noise_source
        not_allowed_keys = ['_url']
        tmp_pars = {k: v for k, v in initializer.noise_pars.as_dict().items() if k not in not_allowed_keys}
        self.parameters = tmp_pars
        self.dt = initializer.resolution
        if isinstance(initializer.start_time, list):
            self.onset_times = initializer.start_time
            self.global_start = initializer.start_time[0]
            self.offset_times = initializer.stop_time
            self.global_stop = initializer.stop_time[-1]
        else:
            self.global_start = initializer.start_time
            self.global_stop = initializer.stop_time

        if start_time is not None:
            self.global_start = start_time
        if stop_time is not None:
            self.global_stop = stop_time

        self.noise_signal = []
        self.time_data = np.arange(self.global_start, self.global_stop, self.dt)

    # TODO @comment
    def generate(self):
        """
        """
        for ii in range(self.N):
            if len(self.source) == 1:
                self.source = list(np.repeat(self.source[0], self.N))
            if not isinstance(self.source[ii], basestring):
                amplitude = self.parameters.pop('amplitude')
                self.parameters.pop('label')
                self.noise_signal.append(self.continuous_rv_generator(self.source[ii], amplitude=amplitude,
                                                                      t_start=self.global_start,
                                                                      t_stop=self.global_stop, dt=self.dt,
                                                                      rectify=self.rectify, array=False,
                                                                      **self.parameters))
            elif isinstance(self.source[ii], basestring) and self.source[ii] not in ['GWN', 'OU']:
                amplitude = self.parameters.pop('amplitude')
                self.parameters.pop('label')
                self.noise_signal.append(self.continuous_rv_generator(self.source[ii], amplitude=amplitude,
                                                                      t_start=self.global_start,
                                                                      t_stop=self.global_stop, dt=self.dt,
                                                                      rectify=self.rectify, array=False,
                                                                      **self.parameters))
            elif self.source[ii] == 'GWN':
                self.noise_signal.append(self.GWN_generator(amplitude=self.parameters['amplitude'],
                                                            mean=self.parameters['mean'],
                                                            std=self.parameters['std'],
                                                            t_start=self.global_start,
                                                            t_stop=self.global_stop,
                                                            dt=self.dt, rectify=self.rectify, array=False))
            elif self.source[ii] == 'OU':
                self.noise_signal.append(self.OU_generator(dt=self.parameters['dt'], tau=self.parameters['tau'],
                                                           sigma=self.parameters['sigma'],
                                                           y0=self.parameters['y0'],
                                                           t_start=self.global_start, t_stop=self.global_stop,
                                                           rectify=self.rectify, array=False))
            else:
                logger.warning("{0} Not currently implemented".format(self.source[ii]))

        channel_ids = []
        sigs = []
        for idx, nn in enumerate(self.noise_signal):
            channel_ids.append(idx * np.ones_like(nn.signal))
            sigs.append(nn.signal)
        sigs = list(itertools.chain(*sigs))
        channel_ids = list(itertools.chain(*channel_ids))

        # generate AnalogSignalList
        if self.N:
            tmp = [(channel_ids[n], sigs[n]) for n in range(len(channel_ids))]
            self.noise_signal = signals.AnalogSignalList(tmp, np.unique(channel_ids).tolist(), times=self.time_data)
        else:
            self.noise_signal = []

    @staticmethod
    def re_seed(global_seed):
        """
        Reset the rng seed (in case it is changed during the signal creation)
        """
        np.random.seed(global_seed) class Network(object):
    """
    Network class is one of the central classes, containing lists of elements specifying
    the properties of each subpopulation (including the population objects)...

    Input:
        - net_pars_set -> ParameterSet object containing all the network parameters

    Population objects are iteratively created when the Network object is initialized
    """
    @staticmethod
    def verify_consistent_dimensions(w, src_gids, tget_gids):
        """
        Verify whether the connection matrix complies with the (source and target) population sizes
        """

        src_dim, tget_dim = w.shape
        if len(src_gids) == src_dim and len(tget_gids) == tget_dim:
            return True
        else:
            return False

    def __init__(self, net_pars_set):
        """
        Initialize network object by creating all sub-populations
        :param net_pars_set: ParameterSet object containing all relevant parameters
        """

        def create_populations(net_pars_set_):
            populations = [[] for _ in range(net_pars_set_.n_populations)]
            # specify the keys not to be passed to the population objects
            not_allowed_keys = ['_url', 'label', 'n_populations', 'parameters', 'names', 'description', '',
                                'integration']

            # iterate through the population list
            logger.info("\nCreating populations:")
            for n in range(net_pars_set_.n_populations):
                if isinstance(net_pars_set_.pop_names[n], list):
                    # Create a composite population (with distinct sub-populations)
                    populations[n] = []
                    for nn in range(len(net_pars_set_.pop_names[n])):
                        # create sub-population dictionary
                        subpop_dict = {k: v[n][nn] for k, v in net_pars_set_.iteritems() if k not in not_allowed_keys
                                       and isinstance(v[n], list)}
                        # if input is not a list, it means that the properties apply to the parent population,
                        # so to all sub-populations:
                        subpop_dict.update({k: v[n] for k, v in net_pars_set_.iteritems() if k not in not_allowed_keys
                                            and (not isinstance(v[n], list))})
                        neuron_dict = subpop_dict['neuron_pars']
                        nest.CopyModel(subpop_dict['neuron_pars']['model'], net_pars_set_.pop_names[n][nn])
                        nest.SetDefaults(net_pars_set_.pop_names[n][nn],
                                         parameters.extract_nestvalid_dict(neuron_dict, param_type='neuron'))
                        if net_pars_set_.topology[n][nn]:
                            tp_dict = subpop_dict['topology_dict']
                            tp_dict.update({'elements': net_pars_set_.pop_names[n][nn]})
                            layer = tp.CreateLayer(tp_dict)
                            gids = nest.GetLeaves(layer)[0]
                            subpop_dict.update({'topology_dict': tp_dict, 'layer_gid': layer,
                                                'is_subpop': True, 'gids': gids})
                            populations[n].append(Population(parameters.ParameterSet(subpop_dict)))
                        else:
                            gids = nest.Create(net_pars_set_.pop_names[n][nn], n=int(net_pars_set_.n_neurons[n][nn]))
                            subpop_dict.update({'is_subpop': False, 'gids': gids})
                            populations[n].append(Population(parameters.ParameterSet(subpop_dict)))

                        logger.info("- Population {0!s}, with ids [{1!s}-{2!s}]".format(net_pars_set_.pop_names[n][nn],
                                                                                  min(gids), max(gids)))
                else:
                    # create a normal population
                    pop_dict = {k: v[n] for k, v in net_pars_set_.iteritems() if k not in not_allowed_keys}
                    neuron_dict = net_pars_set_.neuron_pars[n]

                    # create neuron model named after the population
                    nest.CopyModel(net_pars_set_.neuron_pars[n]['model'], net_pars_set_.pop_names[n])
                    # set default parameters
                    nest.SetDefaults(net_pars_set_.pop_names[n], parameters.extract_nestvalid_dict(neuron_dict,
                                                                                                   param_type='neuron'))

                    if net_pars_set_.topology[n]:
                        tp_dict = pop_dict['topology_dict']
                        tp_dict.update({'elements': net_pars_set_.pop_names[n]})
                        layer = tp.CreateLayer(tp_dict)
                        gids = nest.GetLeaves(layer)[0]
                        pop_dict.update({'topology_dict': tp_dict, 'layer_gid': layer,
                                         'is_subpop': False, 'gids': gids})
                        populations[n] = Population(parameters.ParameterSet(pop_dict))
                    else:
                        # create population
                        gids = nest.Create(net_pars_set_.pop_names[n], n=int(net_pars_set_.n_neurons[n]))
                        # set up population objects
                        pop_dict.update({'gids': gids, 'is_subpop': False})
                        populations[n] = Population(parameters.ParameterSet(pop_dict))
                    logger.info("- Population {0!s}, with ids [{1!s}-{2!s}]".format(net_pars_set_.pop_names[n],
                                                                              min(gids), max(gids)))

                    # store sequence id of population in network
                    populations[n].seq_id = n

            return populations

        # verify_pars_consistency(net_pars_set, ['_url', 'label', 'n_populations', 'parameters', 'names',
        # 									   'analog_device_pars', 'description'], net_pars_set.n_populations)

        self.populations = create_populations(net_pars_set)
        self.n_populations = net_pars_set.n_populations
        self.n_neurons = net_pars_set.n_neurons
        self.population_names = net_pars_set.pop_names
        self.record_spikes = net_pars_set.record_spikes
        self.record_analogs = net_pars_set.record_analogs
        self.spike_device_pars = net_pars_set.spike_device_pars
        self.analog_device_pars = net_pars_set.analog_device_pars
        self.device_names = [[] for _ in range(self.n_populations)]
        self.device_gids = [[] for _ in range(self.n_populations)]
        self.n_devices = [[] for _ in range(self.n_populations)]
        self.device_type = [[] for _ in range(self.n_populations)]
        self.devices_connected = False
        self.spiking_activity = [[] for _ in range(self.n_populations)]
        self.analog_activity = [[] for _ in range(self.n_populations)]
        self.connection_types = []
        self.connection_names = []
        self.synaptic_weights = {}
        self.synaptic_delays = {}
        #self.state_extractors = []
        #self.readouts = []
        self.merged_populations = []

    def get_population_by_name(self, name):
        """
        Returns a population given its name. Raises ValueError if not found.

        :param name: [str] name of Population (e.g. 'E1')
        :return: Population object
        """
        for pop in list(itertools.chain(*[self.merged_populations, self.populations])):
            if pop.name == name:
                return pop
        raise ValueError("Population named '{0}' not found!".format(name))

    def merge_subpopulations(self, sub_populations, name='', merge_activity=False, store=True):
        """
        Combine sub-populations into a main Population object.

        :param sub_populations: [list] - of Population objects to merge
        :param name: [str] - name of new population
        :param merge_activity: [bool] whether to merge previously extracted and currently stored spiking/analog activity
        :param store: [bool] if true, append the newly created population to the list of merged populations
        :return: new Population object
        """
        assert sub_populations, "No sub-populations to merge provided..."
        gids_list = [list(x.gids) for x in sub_populations]
        gids = list(itertools.chain.from_iterable(gids_list))

        subpop_names = [x.name for x in sub_populations]
        if signals.empty(name):
            name = ''.join(subpop_names)

        pop_dict = {'pop_names': name, 'n_neurons': len(gids), 'gids': gids, 'is_subpop': False}

        if all([x.topology for x in sub_populations]):
            positions = [list(x.topology_dict['positions']) for x in sub_populations]
            positions = list(itertools.chain.from_iterable(positions))

            elements = [x.topology_dict['elements'] for x in sub_populations]

            layer_ids = [x.layer_gid for x in sub_populations]

            tp_dict = {'elements': elements, 'positions': positions,
                       'edge_wrap': all(bool('edge_wrap' in x.topology_dict and x.topology_dict['edge_wrap'])
                                        for x in sub_populations)}

            pop_dict.update({'topology': True, 'layer_gid': layer_ids, 'topology_dict': tp_dict})
        else:
            pop_dict.update({'topology': False})

        new_population = Population(parameters.ParameterSet(pop_dict))

        if merge_activity:
            assert (nest.GetKernelStatus()['time']), "There is no activity to merge, run the simulation first and " \
                                                     "merge_population_activity at the network level"
            gids 				= []
            n_neurons 			= np.sum([x.size for x in sub_populations])
            spk_activity_list 	= [x.spiking_activity for x in sub_populations]
            analog_activity 	= [x.analog_activity for x in sub_populations]

            if not any([sl.empty() for sl in spk_activity_list]):
                t_start = round(np.min([x.t_start for x in spk_activity_list]))
                t_stop 	= round(np.max([x.t_stop for x in spk_activity_list]))
                if t_stop > t_start:
                    new_spike_list = signals.SpikeList([], [], t_start=t_start, t_stop=t_stop, dims=n_neurons)

                    for n in spk_activity_list:
                        if not isinstance(n, list):
                            gids.append(n.id_list)
                            for idd in n.id_list:
                                new_spike_list.append(idd, n.spiketrains[idd])
                        else:
                            logger.info("Merge specific spiking activity")   # TODO
                    new_population.spiking_activity = new_spike_list

            if not signals.empty(analog_activity):
                for n in analog_activity:
                    new_population.analog_activity.append(n)

        # store sequence ids of sub-populations which make up the new merged one
        new_population.children = [p.seq_id for p in sub_populations]
        new_population.is_merged = True

        if store:
            self.merged_populations.append(new_population)
        return new_population

    # TODO rewrite all output %s to .format()
    def connect_devices(self):
        """
        Connect recording devices to the populations, according to the parameter specifications

        NOTE: Should only be called once! Otherwise, connections are repeated, although a warning is printed.
        """
        if self.devices_connected:
             warnings.warn('Devices already connected to network, connections will now be repeated! Are you sure?')

        logger.info("\nConnecting Devices: ")
        for n in range(self.n_populations):

            if isinstance(self.record_spikes[n], list):
                # recorder connected to specific sub-populations
                self.n_devices[n] = 0  # number of devices in parent population
                for nn in range(len(self.populations[n])):
                    if self.record_spikes[n][nn]:
                        dev_dict = self.spike_device_pars[n][nn].copy()
                        self.device_type[n].append(dev_dict['model'])
                        dev_dict['label'] += self.population_names[n][nn] + '_' + dev_dict['model']
                        self.device_names[n].append(dev_dict['label'])
                        dev_gid = self.populations[n][nn].record_spikes(parameters.extract_nestvalid_dict(
                            dev_dict, param_type='device'))
                        self.device_gids[n].append(dev_gid)
                        self.n_devices[n] += 1
                        logger.info("- Connecting %s to %s, with label %s and id %s" % (
                            dev_dict['model'], self.population_names[n][nn],
                            dev_dict['label'], str(dev_gid)))
                    if self.record_analogs[n][nn]:
                        dev_dict = self.analog_device_pars[n][nn].copy()
                        self.device_type[n].append(dev_dict['model'])
                        dev_dict['label'] += self.population_names[n][nn] + '_' + dev_dict['model']
                        self.device_names[n].append(dev_dict['label'])
                        if dev_dict['record_n'] != self.populations[n][nn].size:
                            tmp = np.random.permutation(self.n_neurons[n][nn])[:dev_dict['record_n']]
                            ids = []
                            for i in tmp:
                                ids.append(self.populations[n][nn].gids[i])
                        else:
                            ids = None
                        dev_gid = self.populations[n][nn].record_analog(parameters.extract_nestvalid_dict(
                            dev_dict, param_type='device'), ids=ids, record=dev_dict['record_from'])
                        self.device_gids[n].append(dev_gid)
                        self.n_devices[n] += 1
                        logger.info("- Connecting %s to %s %s, with label %s and id %s" % (
                            dev_dict['model'], self.population_names[n][nn], str(ids),
                            dev_dict['label'], str(dev_gid)))

            elif isinstance(self.populations[n], list):
                # recorder connected to large population, including all its subpopulations
                new_pop = self.merge_subpopulations(self.populations[n])
                self.n_devices[n] = 0
                if self.record_spikes[n]:
                    dev_dict = self.spike_device_pars[n].copy()
                    self.device_type[n].append(dev_dict['model'])
                    dev_dict['label'] += new_pop.name + '_' + dev_dict['model']
                    self.device_names[n].append(dev_dict['label'])
                    dev_gid = new_pop.record_spikes(parameters.extract_nestvalid_dict(dev_dict, param_type='device'))
                    self.device_gids[n].append(dev_gid)
                    self.n_devices[n] += 1
                    logger.info("- Connecting %s to %s, with label %s and id %s" % (dev_dict['model'], new_pop.name,
                                                                              dev_dict['label'], str(dev_gid)))
                if self.record_analogs[n]:
                    dev_dict = self.analog_device_pars[n].copy()
                    self.device_type[n].append(dev_dict['model'])
                    dev_dict['label'] += new_pop.name + '_' + dev_dict['model']
                    self.device_names[n].append(dev_dict['label'])
                    if dev_dict['record_n'] != self.populations[n].size:
                        tmp = np.random.permutation(self.n_neurons[n])[:dev_dict['record_n']]
                        ids = []
                        for i in tmp:
                            ids.append(self.populations[n].gids[i])
                    else:
                        ids = None
                    dev_gid = new_pop.record_analog(parameters.extract_nestvalid_dict(dev_dict, param_type='device'),
                                                    ids=ids,
                                                    record=dev_dict['record_from'])
                    self.device_gids[n].append(dev_gid)

                    self.n_devices[n] += 1
                    logger.info("- Connecting %s to %s [%s], with label %s and id %s" %
                          (dev_dict['model'], new_pop.name, str(ids), dev_dict['label'], str(dev_gid)))
            else:
                # there are no sub-populations
                self.n_devices[n] = 0
                if self.record_spikes[n]:
                    dev_dict = self.spike_device_pars[n].copy()
                    self.device_type[n].append(dev_dict['model'])
                    dev_dict['label'] += self.population_names[n] + '_' + dev_dict['model']
                    self.device_names[n].append(dev_dict['label'])
                    dev_gid = self.populations[n].record_spikes(parameters.extract_nestvalid_dict(dev_dict,
                                                                                                  param_type='device'))
                    self.device_gids[n].append(dev_gid)
                    self.n_devices[n] += 1
                    logger.info("- Connecting %s to %s, with label %s and id %s" %
                          (dev_dict['model'], self.population_names[n], dev_dict['label'], str(dev_gid)))
                if self.record_analogs[n]:
                    if isinstance(self.record_analogs[n], bool):
                        dev_dict = self.analog_device_pars[n].copy()
                        self.device_type[n].append(dev_dict['model'])
                        dev_dict['label'] += self.population_names[n] + '_' + dev_dict['model']
                        self.device_names[n].append(dev_dict['label'])
                        if dev_dict['record_n'] != self.populations[n].size:
                            tmp = np.random.permutation(self.n_neurons[n])[:dev_dict['record_n']]
                            ids = []
                            for i in tmp:
                                ids.append(self.populations[n].gids[i])
                        else:
                            ids = None
                        dev_gid = self.populations[n].record_analog(parameters.extract_nestvalid_dict(dev_dict,
                                                                                                      param_type='device'),
                                                                    ids=ids, record=dev_dict['record_from'])
                        self.device_gids[n].append(dev_gid)

                        self.n_devices[n] += 1
                        if (ids is not None) and (len(ids) == 1):
                            logger.info("- Connecting %s to %s [%s], with label %s and id %s" % (dev_dict['model'],
                                                                                           self.population_names[n],
                                                                                           str(ids),
                                                                                           dev_dict['label'],
                                                                                           str(dev_gid)))
                        elif ids is not None:
                            logger.info("- Connecting %s to %s [%s-%s], with label %s and id %s" % (dev_dict['model'],
                                                                                              self.population_names[n],
                                                                                              str(min(ids)),
                                                                                              str(max(ids)),
                                                                                              dev_dict['label'],
                                                                                              str(dev_gid)))
                        else:
                            logger.info("- Connecting %s to %s [%s], with label %s and id %s" % (dev_dict['model'],
                                                                                           self.population_names[n],
                                                                                           str('all'),
                                                                                           dev_dict['label'],
                                                                                           str(dev_gid)))
                    else:
                        for nnn in range(int(self.record_analogs[n])):
                            dev_dict = self.analog_device_pars[n][nnn].copy()
                            self.device_type[n].append(dev_dict['model'])
                            dev_dict['label'] += self.population_names[n] + '_' + dev_dict['model']
                            self.device_names[n].append(dev_dict['label'])
                            if dev_dict['record_n'] != self.populations[n].size:
                                tmp = np.random.permutation(self.n_neurons[n])[:dev_dict['record_n']]
                                ids = []
                                for i in tmp:
                                    ids.append(self.populations[n].gids[i])
                            else:
                                ids = None

                            dev_pars = parameters.extract_nestvalid_dict(dev_dict, param_type='device')
                            dev_gid = self.populations[n].record_analog(dev_pars, ids=ids,
                                                                        record=dev_dict['record_from'])
                            self.device_gids[n].append(dev_gid)

                            self.n_devices[n] += 1
                            if len(ids) == 1:
                                logger.info("- Connecting %s to %s [%s], with label %s and id %s" % (dev_dict['model'],
                                                                                               self.population_names[n],
                                                                                               str(ids),
                                                                                               dev_dict['label'],
                                                                                               str(dev_gid)))
                            else:
                                logger.info("- Connecting %s to %s [%s-%s], with label %s and id %s" % (dev_dict['model'],
                                                                                                  self.population_names[n],
                                                                                                  str(min(ids)),
                                                                                                  str(max(ids)),
                                                                                                  dev_dict['label'],
                                                                                                  str(dev_gid)))

    # TODO maybe break up into several functions?
    def connect_populations(self, connect_pars_set, progress=False, batch_mode=False):
        """
        Connect the sub-populations according to the specifications.

        :param connect_pars_set: ParameterSet for establishing connections
        :param progress:
        :param batch_mode: if connection parameters allow it, connect one presynaptic to all its postsynaptic neurons
                           for efficiency purposes.
        :return:
        """
        logger.info("\nConnecting populations: ")
        # iterate over all synapse types
        for n in range(connect_pars_set.n_synapse_types):
            logger.info("    - %s [%s]" % (connect_pars_set.synapse_types[n], connect_pars_set.models[n]))

            # index of source and target populations in the population lists
            if connect_pars_set.synapse_types[n][1] in self.population_names:
                src_pop_idx = self.population_names.index(connect_pars_set.synapse_types[n][1])
            else:
                src_pop_idx = []
                for ii, nn in enumerate(self.population_names):
                    if isinstance(nn, list):
                        src_pop_idx.append(ii)
                        src_pop_idx.append(nn.index(connect_pars_set.synapse_types[n][1]))
            if connect_pars_set.synapse_types[n][0] in self.population_names:
                tget_pop_idx = self.population_names.index(connect_pars_set.synapse_types[n][0])
            else:
                tget_pop_idx = []
                for ii, nn in enumerate(self.population_names):
                    if isinstance(nn, list):
                        tget_pop_idx.append(ii)
                        tget_pop_idx.append(nn.index(connect_pars_set.synapse_types[n][0]))

            # corresponding global ids of each population to connect, or layer ids if topological connections are
            # present
            if not isinstance(src_pop_idx, int):
                if connect_pars_set.topology_dependent[n]:
                    assert self.populations[src_pop_idx[0]][src_pop_idx[1]].topology, "Source layer doesn't have " \
                                                                                      "topology"
                    src_gids = self.populations[src_pop_idx[0]][src_pop_idx[1]].layer_gid
                else:
                    src_gids = self.populations[src_pop_idx[0]][src_pop_idx[1]].gids
            else:
                if connect_pars_set.topology_dependent[n]:
                    assert self.populations[src_pop_idx].topology, "Source layer doesn't have topology"
                    src_gids = self.populations[src_pop_idx].layer_gid
                else:
                    src_gids = self.populations[src_pop_idx].gids
            if not isinstance(tget_pop_idx, int):
                if connect_pars_set.topology_dependent[n]:
                    assert self.populations[tget_pop_idx[0]][tget_pop_idx[1]].topology, "Target layer doesn't have " \
                                                                                        "topology"
                    tget_gids = self.populations[tget_pop_idx[0]][tget_pop_idx[1]].layer_gid
                else:
                    tget_gids = self.populations[tget_pop_idx[0]][tget_pop_idx[1]].gids
            else:
                if connect_pars_set.topology_dependent[n]:
                    assert self.populations[tget_pop_idx].topology, "Target layer doesn't have topology"
                    tget_gids = self.populations[tget_pop_idx].layer_gid
                else:
                    tget_gids = self.populations[tget_pop_idx].gids

            # copy and modify synapse model
            if hasattr(connect_pars_set, "synapse_names"):
                synapse_name = connect_pars_set.synapse_names[n]
            else:
                synapse_name = connect_pars_set.synapse_types[n][1] + '_' + connect_pars_set.synapse_types[n][0]

            nest.CopyModel(connect_pars_set.models[n], synapse_name)
            nest.SetDefaults(synapse_name, connect_pars_set.model_pars[n])
            self.connection_names.append(synapse_name)
            self.connection_types.append(connect_pars_set.synapse_types[n])

            #**** set up connections ****
            if 'copy' in synapse_name:  # re-connect the same neurons...
                start = time.time()
                logger.info("    - Connecting {0} (*)".format(synapse_name))

                device_models = ['spike_detector', 'spike_generator', 'multimeter']
                target_synapse_name = synapse_name[:synapse_name.find('copy')-1]
                conns = nest.GetConnections(synapse_model=target_synapse_name)

                iterate_steps = 100
                its = np.arange(0, len(conns)+1, iterate_steps).astype(int)

                for nnn, it in enumerate(its):
                    if nnn < len(its)-1:
                        con = conns[it:its[nnn+1]]
                        st = nest.GetStatus(con)
                        source_gids = [x['source'] for x in st if
                                       nest.GetStatus([x['target']])[0]['model'] not in device_models]
                        target_gids = [x['target'] for x in st if
                                       nest.GetStatus([x['target']])[0]['model'] not in device_models]

                        syn_dict = parameters.copy_dict(connect_pars_set.syn_specs[n], {'model': synapse_name, 'weight':
                            connect_pars_set.weight_dist[n], 'delay': connect_pars_set.delay_dist[n]})

                        syn_dicts = [{'synapsemodel': list(np.repeat(synapse_name, len(source_gids)))[iddx],
                                      'source': source_gids[iddx],
                                      'target': target_gids[iddx],
                                      'weight': syn_dict['weight'],		              # TODO distributions??
                                      'delay': syn_dict['delay'],
                                      'receptor_type': syn_dict['receptor_type']} for iddx in range(len(target_gids))]
                        nest.DataConnect(syn_dicts)
                    if progress:
                        visualization.progress_bar(float(nnn) / float(len(its)))
                logger.info("\tElapsed time: {0} s".format(str(time.time()-start)))
            else:
                # 1) if pre-computed weights matrices are given
                if (connect_pars_set.pre_computedW[n] is not None) and (not connect_pars_set.topology_dependent[n]):

                    if isinstance(connect_pars_set.pre_computedW[n], str):
                        w = np.load(connect_pars_set.pre_computedW[n])
                    else:
                        w = connect_pars_set.pre_computedW[n]

                    if self.verify_consistent_dimensions(w, src_gids, tget_gids):
                        for preSyn_matidx, preSyn_gid in enumerate(src_gids):
                            postSyn_matidx = w[preSyn_matidx, :].nonzero()[0]
                            postSyn_gids = list(postSyn_matidx + min(tget_gids))
                            weights = [w[preSyn_matidx, x] for x in postSyn_matidx]

                            if not isinstance(connect_pars_set.delay_dist[n], dict) and isinstance(
                                    connect_pars_set.delay_dist[n], list):
                                delays = connect_pars_set.delay_dist[n]
                            elif isinstance(connect_pars_set.delay_dist[n], dict):
                                delays = [parameters.copy_dict(connect_pars_set.delay_dist[n]) for _ in range(len(
                                    weights))]
                            else:
                                delays = np.repeat(connect_pars_set.delay_dist[n], len(weights))

                            if not batch_mode:
                                # we can do this more efficiently in batch mode if parameters allow it, but in this
                                # case iterate through each connection one-by-one
                                for idx, tget in enumerate(postSyn_gids):
                                    update_dict = {'model': synapse_name, 'weight': weights[idx], 'delay': delays[idx]}
                                    syn_dict = parameters.copy_dict(connect_pars_set.syn_specs[n], update_dict)
                                    nest.Connect([preSyn_gid], [tget], syn_spec=syn_dict)

                            # TODO is it not always possible to do batch mode? Weight and delay matrices are computed
                            # TODO above, separately, so this shouldn't affect
                            elif len(postSyn_gids):
                                # do this in batch mode, calling nest.Connect only once, if there's only one synapse
                                # model. Weight and delay matrices have been computed above. Reshaping is necessary
                                # for NEST..
                                update_dict = {'model': synapse_name,
                                               'weight': np.reshape(weights, (len(weights), 1)),
                                               'delay': np.reshape(delays, (len(delays), 1))}

                                syn_dict = parameters.copy_dict(connect_pars_set.syn_specs[n], update_dict)
                                nest.Connect([preSyn_gid], postSyn_gids, syn_spec=syn_dict)

                            if progress:
                                visualization.progress_bar(float(preSyn_matidx)/float(len(src_gids)))
                    else:
                        raise Exception("Dimensions of W are inconsistent with population sizes")

                # 2) if no pre-computed weights, and no topology in pre/post-synaptic populations, use dictionaries
                elif (connect_pars_set.pre_computedW[n] is None) and (not connect_pars_set.topology_dependent[n]):

                    syn_dict = parameters.copy_dict(connect_pars_set.syn_specs[n], {'model': synapse_name, 'weight':
                        connect_pars_set.weight_dist[n], 'delay': connect_pars_set.delay_dist[n]})
                    conn_dict = connect_pars_set.conn_specs[n]

                    nest.Connect(src_gids, tget_gids, conn_spec=conn_dict, syn_spec=syn_dict)

                # 3) if no precomputed weights, but topology in populations and topological connections
                elif connect_pars_set.topology_dependent[n]:
                    assert nest.is_sequence_of_gids(src_gids) and len(src_gids) == 1, \
                        "Source ids are not topology layer"
                    assert nest.is_sequence_of_gids(tget_gids) and len(tget_gids) == 1, \
                        "Target ids are not topology layer"
                    tp_dict = parameters.copy_dict(connect_pars_set.conn_specs[n], {'synapse_model': synapse_name})
                    tp.ConnectLayers(src_gids, tget_gids, tp_dict)

    def flush_records(self, decoders=False):
        """
        Delete all data from all devices connected to the network.

        :param decoders: [bool] if true, delete data stored in the decoding layer as well
        :return:
        """
        if not signals.empty(self.device_names):
            logger.info("\nClearing device data: ")
        devices = list(itertools.chain.from_iterable(self.device_names))

        for idx, n in enumerate(list(itertools.chain.from_iterable(self.device_gids))):
            logger.info(" - {0} {1}".format(devices[idx], str(n)))
            nest.SetStatus(n, {'n_events': 0})
            if nest.GetStatus(n)[0]['to_file']:
                io.remove_files(nest.GetStatus(n)[0]['filenames'])

        if decoders:
            for ctr, n_pop in enumerate(list(itertools.chain(*[self.merged_populations, self.populations]))):
                if n_pop.decoding_layer is not None:
                    n_pop.decoding_layer.flush_records()

    def copy(self):
        """
        Returns a copy of the entire network object. Doesn't create new NEST objects.
        :return:
        """
        return copy.deepcopy(self)

    # TODO comment
    def clone(self, original_parameter_set, devices=True, decoders=True, display=True):
        """
        Creates a new network object.

        :param original_parameter_set:
        :param devices:
        :param decoders:
        :param display:
        :return:
        """
        param_set = copy.deepcopy(original_parameter_set)

        def create_clone(net):
            """
            Create Network object
            :return:
            """
            neuron_pars = [0 for _ in range(len(list(signals.iterate_obj_list(net.populations))))]
            topology = [0 for _ in range(len(list(signals.iterate_obj_list(net.populations))))]
            topology_dict = [0 for _ in range(len(list(signals.iterate_obj_list(net.populations))))]
            pop_names = [0 for _ in range(len(list(signals.iterate_obj_list(net.populations))))]
            spike_device_pars = [0 for _ in range(len(list(signals.iterate_obj_list(net.populations))))]
            analog_dev_pars = [0 for _ in range(len(list(signals.iterate_obj_list(net.populations))))]
            status_elements = ['archiver_length', 'element_type', 'frozen', 'global_id',
                               'has_connections', 'local', 'recordables', 't_spike',
                               'thread', 'thread_local_id', 'vp', 'n_synapses',
                               'local_id', 'model', 'parent']
            for idx_pop, pop_obj in enumerate(list(signals.iterate_obj_list(net.populations))):
                # get generic neuron_pars (update individually later):
                neuron = parameters.extract_nestvalid_dict(nest.GetStatus([pop_obj.gids[0]])[0], param_type='neuron')
                d_tmp = {k: v for k, v in neuron.items() if k not in status_elements}
                neuron_pars[idx_pop] = parameters.copy_dict(d_tmp, {'model': nest.GetStatus([pop_obj.gids[0]])[0][
                    'model']})
                if isinstance(pop_obj.topology, dict):
                    topology[idx_pop] = True
                    topology_dict[idx_pop] = pop_obj.topology
                else:
                    topology[idx_pop] = False
                    topology_dict[idx_pop] = pop_obj.topology
                pop_names[idx_pop] = net.population_names[idx_pop] + '_clone'
                if net.record_spikes[idx_pop]:
                    spike_device_pars[idx_pop] = parameters.copy_dict(net.spike_device_pars[idx_pop],
                                                                      {'label': net.spike_device_pars[idx_pop]['label']+'_clone'})
                else:
                    spike_device_pars[idx_pop] = None
                if net.record_analogs[idx_pop]:
                    analog_dev_pars[idx_pop] = parameters.copy_dict(net.analog_device_pars[idx_pop],
                                                                    {'label': net.analog_device_pars[idx_pop]['label']+'_clone'})
                else:
                    analog_dev_pars[idx_pop] = None

            network_parameters = {'n_populations': net.n_populations,
                                  'pop_names': pop_names,
                                  'n_neurons': net.n_neurons,
                                  'neuron_pars': neuron_pars,
                                  'topology': topology,
                                  'topology_dict': topology_dict,
                                  'record_spikes': net.record_spikes,
                                  'spike_device_pars': spike_device_pars,
                                  'record_analogs': net.record_analogs,
                                  'analog_device_pars': analog_dev_pars}
            clone_net = Network(parameters.ParameterSet(network_parameters, label='clone'))
            #clone_net.connect_devices()

            for pop_idx, pop_obj in enumerate(clone_net.populations):
                for n_neuron in range(net.n_neurons[pop_idx]):
                    src_gid = net.populations[pop_idx].gids[n_neuron]
                    tget_gid = pop_obj.gids[n_neuron]
                    status_dict = nest.GetStatus([src_gid])[0]
                    st = {k: v for k, v in status_dict.items() if k not in status_elements}
                    nest.SetStatus([tget_gid], st)
            return clone_net

        def connect_clone(network, clone, progress=False):
            """
            Connect the populations in the clone network
            (requires iteration through all the connections in the mother network...)
            :param network:
            :param clone:
            :return:
            """
            devices = ['stimulator', 'structure']
            base_idx = min(list(itertools.chain(*[n.gids for n in clone.populations])))-1
            logger.info("\n Replicating connectivity in clone network (*)")
            for syn_idx, synapse in enumerate(network.connection_names):
                start = time.time()
                copy_synapse_name = synapse + '_clone'
                nest.CopyModel(synapse, copy_synapse_name)
                logger.info("\t- {0}".format(str(network.connection_types[syn_idx])))
                conns = nest.GetConnections(synapse_model=synapse)
                # ##
                iterate_steps = 100
                its = np.arange(0, len(conns)+iterate_steps, iterate_steps).astype(int)
                its[-1] = len(conns)
                # ##
                clone.connection_names.append(copy_synapse_name)
                conn_types = network.connection_types[syn_idx]
                connection_type = (conn_types[0]+'_clone', conn_types[1]+'_clone')
                clone.connection_types.append(connection_type)

                # src_idx = clone.population_names.index(connection_type[1])
                # tget_idx = clone.population_names.index(connection_type[0])
                # src_gids = clone.populations[src_idx].gids
                # tget_gids = clone.populations[tget_idx].gids
                # base_idx = min(list(itertools.chain(*[src_gids, tget_gids]))) - 1

                for nnn, it in enumerate(its):
                    if nnn < len(its) - 1:
                        con = conns[it:its[nnn+1]]
                        st = nest.GetStatus(con)
                        source_gids = [x['source']+base_idx for x in st if nest.GetDefaults(nest.GetStatus([x[
                                                                                                                'target']])[0]['model'])['element_type'] not in
                                       devices]
                        target_gids = [x['target']+base_idx for x in st if nest.GetDefaults(nest.GetStatus([x[
                                                                                                                'target']])[0]['model'])['element_type'] not in
                                       devices]
                        weights = [x['weight'] for x in st if nest.GetDefaults(nest.GetStatus([x['target']])[0][
                                                                                   'model'])['element_type'] not in devices]
                        delays = [x['delay'] for x in st if nest.GetDefaults(nest.GetStatus([x['target']])[0][
                                                                                 'model'])['element_type'] not in devices]

                        receptors = [x['receptor'] for x in st if nest.GetDefaults(nest.GetStatus([x['target']])[0][
                                                                                       'model']).has_key('element_type') and nest.GetDefaults(nest.GetStatus([x['target']])[0][
                                                                                                                                                  'model'])['element_type'] not in devices]
                        # modify target receptors in cases where they are used:
                        # for iddx, n_rec in enumerate(receptors):
                        # 	# check if neuron accepts different receptors
                        # 	if nest.GetStatus([target_gids[iddx]])[0].has_key('rec_type'):
                        # 		rec_types = list(nest.GetStatus([target_gids[iddx]])[0]['rec_type'])
                        # 		receptors[iddx] = int(nest.GetStatus([target_gids[iddx]])[0]['rec_type'][n_rec])
                        syn_dicts = [{'synapse_model': list(np.repeat(copy_synapse_name, len(source_gids)))[iddx],
                                      'source': source_gids[iddx],
                                      'target': target_gids[iddx],
                                      'weight': weights[iddx],
                                      'delay': delays[iddx],
                                      'receptor_type': receptors[iddx]} for iddx in range(len(target_gids))]
                        nest.DataConnect(syn_dicts)
                        if progress:
                            visualization.progress_bar(float(nnn)/float(len(its)))
                logger.info("\tElapsed time: {0} s".format(str(time.time() - start)))

        def connect_decoders(network, parameters):
            """
            :param parameters:
            :return:
            """
            target_populations = parameters.decoding_pars.state_extractor.source_population
            copy_targets = [n+'_clone' for n in target_populations]
            parameters.decoding_pars.state_extractor.source_population = copy_targets
            network.connect_decoders(parameters.decoding_pars)

        clone_network = create_clone(self)
        connect_clone(self, clone_network, display)
        if devices:
            clone_network.connect_devices()
        if decoders:
            connect_decoders(clone_network, param_set)

        return clone_network

    def mirror(self, copy_net=None, from_main=True, to_main=False):
        """
        Returns a network object equal to self and either connected
        to main network or receiving connections from main network
        :return:
        """
        if copy_net is None:
            copy_net = self.copy()
            cn = None
        else:
            assert isinstance(copy_net, Network), "copy_net must be Network object"
            cn = 1.

        if to_main:
            logger.info("Connecting CopyNetwork: ")
            device_models = ['spike_detector', 'spike_generator', 'multimeter']
            for pop_idx, pop_obj in enumerate(copy_net.populations):
                original_gid_range = [min(self.populations[pop_idx].gids),
                                      max(self.populations[pop_idx].gids)]
                copy_gid_range = [min(pop_obj.gids), max(pop_obj.gids)]

                start = time.time()
                logger.info("    - {0}, {1}".format(copy_net.population_names[pop_idx], self.population_names[pop_idx]))

                for n_neuron in range(self.n_neurons[pop_idx]):
                    src_gid = self.populations[pop_idx].gids[n_neuron]
                    tget_gid = pop_obj.gids[n_neuron]

                    conns = nest.GetConnections(source=[src_gid])

                    # for memory conservation, iterate:
                    iterate_steps = 100
                    its = np.arange(0, len(conns), iterate_steps).astype(int)

                    for nnn, it in enumerate(its):
                        if nnn < len(its)-1:
                            conn = conns[it:its[nnn+1]]
                            st = nest.GetStatus(conn)

                        target_gids = [x['target'] for x in st if
                                       nest.GetStatus([x['target']])[0]['model'] not in device_models]
                        weights = [x['weight'] for x in st if
                                   nest.GetStatus([x['target']])[0]['model'] not in device_models]
                        delays = [x['delay'] for x in st if nest.GetStatus([x['target']])[0]['model'] not in device_models]
                        models = [x['synapse_model'] for x in st if
                                  nest.GetStatus([x['target']])[0]['model'] not in device_models]
                        receptors = [x['receptor'] for x in st if
                                     nest.GetStatus([x['target']])[0]['model'] not in device_models]

                        tgets = [x + (copy_gid_range[0] - original_gid_range[0]) for x in target_gids]
                        syn_dicts = [{'synapsemodel': models[iddx], 'source': tget_gid,
                                      'target': tgets[iddx], 'weight': weights[iddx],
                                      'delay': delays[iddx], 'receptor_type': receptors[iddx]} for iddx in range(len(
                            target_gids))]
                        nest.DataConnect(syn_dicts)
                logger.info("Elapsed Time: {0}".format(str(time.time()-start)))

        elif from_main:
            logger.info("\nConnecting CopyNetwork: ")
            device_models = ['spike_detector', 'spike_generator', 'multimeter']
            for pop_idx, pop_obj in enumerate(copy_net.populations):
                original_gid_range = [min(self.populations[pop_idx].gids),
                                      max(self.populations[pop_idx].gids)]
                copy_gid_range = [min(pop_obj.gids), max(pop_obj.gids)]

                start = time.time()
                logger.info("\t    - {0}, {1}".format(copy_net.population_names[pop_idx], self.population_names[pop_idx]))

                for n_neuron in range(self.n_neurons[pop_idx]):
                    src_gid = self.populations[pop_idx].gids[n_neuron]
                    tget_gid = pop_obj.gids[n_neuron]

                    conns = nest.GetConnections(target=[src_gid])

                    # for memory conservation, iterate:
                    iterate_steps = 100
                    its = np.arange(0, len(conns), iterate_steps).astype(int)

                    for nnn, it in enumerate(its):
                        if nnn < len(its)-1:
                            conn = conns[it:its[nnn+1]]
                            st = nest.GetStatus(conn)
                            source_gids = [x['source'] for x in st if
                                           nest.GetStatus([x['source']])[0]['model'] not in device_models]
                            weights = [x['weight'] for x in st if
                                       nest.GetStatus([x['target']])[0]['model'] not in device_models]
                            delays = [x['delay'] for x in st if nest.GetStatus([x['target']])[0]['model'] not in device_models]
                            models = [x['synapse_model'] for x in st if
                                      nest.GetStatus([x['target']])[0]['model'] not in device_models]
                            receptors = [x['receptor'] for x in st if
                                         nest.GetStatus([x['target']])[0]['model'] not in device_models]

                            sources = [x + (copy_gid_range[0] - original_gid_range[0]) for x in source_gids]
                            syn_dicts = [{'synapsemodel': models[iddx], 'source': sources[iddx],
                                          'target': tget_gid, 'weight': weights[iddx],
                                          'delay': delays[iddx], 'receptor_type': receptors[iddx]} for iddx in range(len(
                                source_gids))]
                            nest.DataConnect(syn_dicts)
                logger.info("Elapsed Time: {0}".format(str(time.time() - start)))
        if cn is None:
            return copy_net

    def simulate(self, t=10, reset_devices=False, reset_network=False):
        """
        Simulate network
        :param t: total time
        :return:
        """
        nest.Simulate(t)
        if reset_devices:
            self.flush_records()
        if reset_network:
            nest.ResetNetwork()

    def extract_population_activity(self, t_start=None, t_stop=None):
        """
        Iterate through the populations in the network, verify which recording devices have been connected to them
        and extract this data. The activity is then converted in SpikeList or AnalogList objects and attached to
        the properties of the corresponding population.
        To merge the data from multiple populations see extract_network_activity()
        """
        if not signals.empty(self.device_names):
            logger.info("\nExtracting and storing recorded activity from devices:")

        for p_idx, p in enumerate(self.populations):
            if isinstance(p, list):
                for pop_idx, pop in enumerate(p):
                    if pop.attached_devices:
                        logger.info("- Population {0}".format(pop.name))
                        for nnn in pop.attached_devices:
                            if nest.GetStatus(nnn)[0]['to_memory']:
                                # initialize load_activity with device gid
                                pop.load_activity(nnn, t_start=t_start, t_stop=t_stop)
                            elif nest.GetStatus(nnn)[0]['to_file']:
                                # initialize load_activity with file paths
                                pop.load_activity(list(nest.GetStatus(nnn)[0]['filenames']),
                                                  t_start=t_start, t_stop=t_stop)
            else:
                if p.attached_devices:
                    logger.info("- Population {0}".format(p.name))
                    for pop in p.attached_devices:
                        if nest.GetStatus(pop)[0]['to_memory']:
                            # initialize load_activity with device gid
                            p.load_activity(pop, t_start=t_start, t_stop=t_stop)
                        elif nest.GetStatus(pop)[0]['to_file']:
                            p.load_activity(list(nest.GetStatus(pop)[0]['filenames']), t_start=t_start, t_stop=t_stop)

    def extract_network_activity(self):
        """
        Combine the activity lists attached to each population into a list of SpikeList or AnalogList objects
        corresponding to each sub-population
        NOTE: This function will copy the activity lists contained in the lower-level Population objects into the
        properties of the Network object (it may be redundant to maintain both activity sets, and it's only useful in
        certain situations)
        """
        for n in range(self.n_populations):
            if isinstance(self.populations[n], list) and self.n_devices[n]:
                for nn in range(len(self.populations[n])):
                    self.spiking_activity[n].append(self.populations[n][nn].spiking_activity)
                    self.analog_activity[n].append(self.populations[n][nn].analog_activity)
            elif self.n_devices[n]:
                self.spiking_activity[n] = self.populations[n].spiking_activity
                self.analog_activity[n] = self.populations[n].analog_activity

    def merge_population_activity(self, start=0., stop=1000., save=False):
        """
        Merge spike and analog data from the different populations, according to the already merged
        and registered populations. NO global merging here!

        :param start: start time of merging
        :param stop: end time of merging (included?)
        :param save: [deprecated] not used anymore, only for backwards compatibility purposes
        :return:
        """
        if not signals.empty(self.merged_populations):
            for m_pop in self.merged_populations:
                assert m_pop.is_merged

                m_spiking_activity = signals.SpikeList([], [], start, stop, m_pop.size)
                m_analog_activity = []

                # iterate over each parent population (ids) in the merged one and collect spiking and analog activities
                for child_seq_id in m_pop.children:
                    # only add spiking activity if it's not empty
                    if not self.populations[child_seq_id].spiking_activity.empty():
                        child_spiking_activity = self.populations[child_seq_id].spiking_activity
                        # truncate to correct start and end times if necessary
                        if not m_spiking_activity.time_parameters() == child_spiking_activity.time_parameters():
                            child_spiking_activity = child_spiking_activity.time_slice(start, stop)
                        m_spiking_activity.concatenate(child_spiking_activity)

                    # add analog activity
                    m_analog_activity.append(self.populations[child_seq_id].analog_activity)

                m_pop.spiking_activity = m_spiking_activity
                m_pop.analog_activity = m_analog_activity

    def extract_synaptic_weights(self, src_gids=None, tget_gids=None):
        """
        Read and store the weights for all the connected populations, or just for
        the provided sources and targets.

        :param src_gids: list of gids of source neurons
        :param tget_gids: list of gids of target neurons
        """

        if src_gids is not None and tget_gids is not None:
            syn_name = str(nest.GetStatus(nest.GetConnections([src_gids[0]], [tget_gids[0]]))[0]['synapse_model'])
            self.synaptic_weights.update({syn_name: extract_weights_matrix(src_gids, tget_gids)})
        else:
            for nn in self.connection_types:
                src_idx = self.population_names.index(nn[1])
                tget_idx = self.population_names.index(nn[0])
                src_gids = self.populations[src_idx].gids
                tget_gids = self.populations[tget_idx].gids
                self.synaptic_weights.update({nn: extract_weights_matrix(src_gids, tget_gids)})

    def extract_synaptic_delays(self, src_gids=None, tget_gids=None):
        """
        Read and store the synaptic delays for all the connected populations, or just for
        the provided sources and targets.

        :param src_gids: list of gids of source neurons
        :param tget_gids: list of gids of target neurons
        :return:
        """
        if src_gids and tget_gids:
            syn_name = str(nest.GetStatus(nest.GetConnections([src_gids[0]], [tget_gids[0]]))[0]['synapse_model'])
            self.synaptic_weights.update({syn_name: extract_weights_matrix(src_gids, tget_gids)})
        else:
            for nn in self.connection_types:
                src_idx = self.population_names.index(nn[1])
                tget_idx = self.population_names.index(nn[0])
                src_gids = self.populations[src_idx].gids
                tget_gids = self.populations[tget_idx].gids
                self.synaptic_delays.update({nn: extract_delays_matrix(src_gids, tget_gids)})

    def connect_decoders(self, decoding_pars):
        """
        Connect a Decoding Layer to the current network object.

        :param: decoding_pars: [ParameterSet or dict]
        :return:
        """
        if isinstance(decoding_pars, dict):
            decoding_pars = parameters.ParameterSet(decoding_pars)
        assert isinstance(decoding_pars, parameters.ParameterSet), "DecodingLayer must be initialized with " \
                                                                   "ParameterSet or dictionary"

        population_names 		= list(signals.iterate_obj_list(self.population_names))
        population_objs 		= list(signals.iterate_obj_list(self.populations))
        merged_population_names = [x.name for x in self.merged_populations]
        merged_population_objs 	= [x for x in self.merged_populations]
        decoder_params 			= {}

        # initialize state extractors:
        if hasattr(decoding_pars, "state_extractor"):
            pars_st = decoding_pars.state_extractor
            logger.info("\nConnecting Decoders: ")

            # group state_extractors by source population
            sources = []
            source_populations = []
            for ext_idx, n_src in enumerate(pars_st.source_population):
                if isinstance(n_src, list):
                    pop_label = ''.join(n_src)
                    if pop_label not in merged_population_names:
                        sub_population_idx = [population_names.index(x) for x in n_src]
                        sub_populations = [population_objs[x] for x in sub_population_idx]
                        source_population = self.merge_subpopulations(sub_populations=sub_populations, name=pop_label,
                                                                      merge_activity=True)
                    else:
                        source_population = merged_population_objs[merged_population_names.index(pop_label)]
                else:
                    pop_label = n_src
                    pop_index = population_names.index(n_src)
                    source_population = population_objs[pop_index]
                sources.append(pop_label)
                source_populations.append(source_population)
            unique_sources = set(sources)

            extractor_indices = {}
            for n_ext in list(unique_sources):
                extractor_indices.update({n_ext: np.where(np.array(sources) == n_ext)[0]})

            # create decoder parameters dictionary for each source population
            keys = ['state_specs', 'state_variable', 'reset_states', 'average_states', 'standardize']
            for k, v in extractor_indices.items():
                decoder_params.update({k: {}})
                for k1 in keys:
                    decoder_params[k].update({k1: [pars_st[k1][x] for x in v]})

            if hasattr(decoding_pars, "readout"):
                pars_rd = decoding_pars.readout
                for k, v in extractor_indices.items():
                    decoder_params[k].update({'readout': [pars_rd[x] for x in v]})
        else:
            raise IOError("DecodingLayer requires the specification of state extractors")

        pops = [source_populations[sources.index(x)] for x in unique_sources]
        for (population_name, population) in zip(unique_sources, pops):
            population.connect_decoders(parameters.ParameterSet(decoder_params[population_name])) def set_storage_locations(parameter_set, save=True):
    """
    Define paths to store data
    :param parameter_set: ParameterSet object
    :param save: [bool] is False, no paths are created
    :return save_paths: dictionary containing all relevant storage locations
    """
    if save:
        print("\nSetting storage paths...")
        save_label = parameter_set.label

        main_folder = parameter_set.kernel_pars.data_path + parameter_set.kernel_pars.data_prefix + '/'

        figures = main_folder + 'Figures/'
        inputs = main_folder + 'Inputs/'
        parameters = main_folder + 'Parameters/'
        results = main_folder + 'Results/'
        activity = main_folder + 'Activity/'
        others = main_folder + 'Other/'
        logs = main_folder + 'Logs/'

        dirs = {'main': main_folder, 'figures': figures, 'inputs': inputs, 'parameters': parameters, 'results':
            results, 'activity': activity, 'other': others, 'logs': logs}

        for d in dirs.values():
            try:
                os.makedirs(d)
            except OSError:
                pass

        dirs['label'] = save_label

        return dirs
    else:
        print("No data will be saved!")
        return {'label': False, 'figures': False, 'activity': False} def iterate_obj_list(obj_list):
    """
    Build an iterator to iterate through any nested list
    :obj_list: list of objects to iterate
    :return:
    """
    for idx, n in enumerate(obj_list):
        if isinstance(n, list):
            for idxx, nn in enumerate(obj_list[idx]):
                yield obj_list[idx][idxx]
        else:
            yield obj_list[idx] def set_global_rcParams(rc_config_file):
    """
    Replace the rcParams configuration of matplotlib, with predefined values in the dictionary rc_dict. To restore
    the default values, call mpl.rcdefaults()

    For more details, consult the documentation of rcParams

    :param rc_config_file: path to file containing detailed mpl configuration
    """
    assert os.path.isfile(rc_config_file), 'input must be path to config file'
    mpl.rc_file(rc_config_file) class InputPlots(object):
    """
    Class of plotting routines for input structures

    Input:
    ------
    connection_pars - ParameterSet containing all the connectivity parameters
    populations - list of Population objects in the network to plot
    colors      - color associated with each population
    """

    def __init__(self, stim_obj=None, input_obj=None, noise_obj=None):
        """
        """
        self.stim = stim_obj
        self.input = input_obj
        self.noise = noise_obj

    def plot_stimulus_matrix(self, set='train', ax=None, save=False, display=True):
        """

        """
        if (ax is not None) and (not isinstance(ax, mpl.axes.Axes)):
            raise ValueError('ax must be matplotlib.axes.Axes instance.')
        if ax is None:
            fig, ax = pl.subplots()
            fig.suptitle('{0} Stimulus Set'.format(set))

        data = getattr(self.stim, set + '_set')
        labels = getattr(self.stim, set + '_set_labels')
        plt1 = ax.imshow(1 - data.todense(), interpolation='nearest', aspect='auto', extent=None,
                         cmap='gray')
        ax.set_xlabel('Time Step')
        ax.set_ylabel('Stimulus')
        ax.set_xticks(np.arange(len(labels)))
        ax.set_yticks(np.arange(len(np.unique(np.array(labels)))))
        ax.set_yticklabels(np.unique(np.array(labels)))

        if display:
            pl.show(False)
        if save:
            assert isinstance(save, str), "Please provide filename"
            pl.savefig(save+'_StimulusMatrix_{0}.pdf'.format(str(set)))

    def plot_input_signal(self, ax=None, save=False, display=True, **kwargs):
        """
        """
        ax_props = {k: v for k, v in kwargs.iteritems() if k in ax.properties()}
        plot_props = {k: v for k, v in kwargs.iteritems() if k not in ax.properties()}

        if (ax is not None) and (not isinstance(ax, mpl.axes.Axes)):
            raise ValueError('ax must be matplotlib.axes.Axes instance.')
        if ax is None:
            fig = pl.figure()
            for ii in range(len(self.input.input_signal)):
                ax = pl.subplot2grid((len(self.input.input_signal), 1), (ii, 0), rowspan=1, colspan=1)
                ax.plot(self.input.time_data, self.input.input_signal[ii].signal, **plot_props)
                ax.set_ylabel(r'$\sigma_{u}$')
                if ii < len(self.input.input_signal) - 1:
                    ax.set_xticklabels('')
                    ax.set_xlim(left=min(self.input.time_data), right=max(self.input.time_data))
                else:
                    ax.set_xlabel('Time [ms]')
                    ax.set_xlim(left=min(self.input.time_data), right=max(self.input.time_data))
                ax.set_ylim([min(self.input.input_signal[ii].signal)-10., max(self.input.input_signal[ii].signal)+10.])
                ax.set(**ax_props)
        else:
            for ii in range(len(self.input.input_signal)):
                ax.plot(self.input.time_data, self.input.input_signal[ii].signal, **plot_props)
                ax.set_xlim([min(self.input.time_data), max(self.input.time_data)])
                ax.set_ylim([min(self.input.input_signal[ii].signal)-10., max(self.input.input_signal[ii].signal)+10.])
                ax.set_xlabel('Time [ms]')
                ax.set_ylabel(r'$\sigma_{u}$')
                ax.set(**ax_props)
        if display:
            pl.show(False)
        if save:
            assert isinstance(save, str), "Please provide filename"
            if not os.path.exists(save+'_InputTimeSeries.pdf'):
                pl.savefig(save+'_InputTimeSeries.pdf')
            else: # cases when the function is called twice
                pl.savefig(save+'_InputTimeSeries2.pdf')

    def plot_noise_component(self, ax=None, save=False, display=True):
        """
        """
        if (ax is not None) and (not isinstance(ax, mpl.axes.Axes)):
            raise ValueError('ax must be matplotlib.axes.Axes instance.')
        if ax is None:
            fig = pl.figure()
            for ii in range(len(self.noise.noise_signal)):
                ax = pl.subplot2grid((len(self.noise.noise_signal), 1), (ii, 0), rowspan=1,
                                     colspan=1)
                ax.plot(self.noise.noise_signal[ii].time_axis(), self.noise.noise_signal[
                    ii].signal)
                ax.set_ylabel(r'$\sigma_{u}$')
                if ii < len(self.noise.noise_signal) - 1:
                    ax.set_xticklabels('')
                    ax.set_xlim(min(self.noise.time_data), max(self.noise.time_data))
                else:
                    ax.set_xlabel('Time [ms]')
                    ax.set_xlim(min(self.noise.time_data), max(self.noise.time_data))
        else:
            for ii in range(len(self.noise.noise_signal)):
                ax.plot(self.noise.noise_signal[ii].time_axis(), self.noise.noise_signal[ii].signal)
            # if self.start and self.stop:
            # 	ax.set_xlim(left=self.start, right=self.stop)
            # else:
            ax.set_xlim(left=min(self.input.time_data), right=max(self.input.time_data))
        if display:
            pl.show(False)
        if save:
            assert isinstance(save, str), "Please provide filename"
            pl.savefig(save+'NoiseTimeSeries.pdf')

    def plot_signal_and_noise(self, channel_index=0, save=False, display=True):
        """
        """
        if self.noise is not None:
            combined = self.input.input_signal[channel_index].signal + self.noise.noise_signal[channel_index].signal
        else:
            return

        fig = pl.figure()
        ax1 = pl.subplot2grid((3, 1), loc=(0, 0), rowspan=1, colspan=1)
        ax1.set_title('Signal')
        ax1.set_xticklabels('')
        ax1.set_xlim(left=min(self.input.time_data), right=max(self.input.time_data))

        ax2 = pl.subplot2grid((3, 1), loc=(1, 0), rowspan=1, colspan=1)
        ax2.set_xticklabels('')
        snr = st.signaltonoise(combined)
        ax2.set_title(r'Signal [{0}] + Noise $\Rightarrow$ SNR = {1}'.format(str(channel_index), str(snr)))
        ax2.set_xlim(left=min(self.input.time_data), right=max(self.input.time_data))
        ax3 = pl.subplot2grid((3, 1), loc=(2, 0), rowspan=1, colspan=1)
        ax3.set_title('Noise')
        ax3.set_xlabel('Time [ms]')
        ax3.set_xlim(min(self.noise.time_data), max(self.noise.time_data))

        ax1.plot(self.input.time_data, self.input.input_signal[channel_index].signal)
        ax2.plot(self.input.time_data, self.input.input_signal[channel_index].signal + self.noise.noise_signal[
            channel_index].signal)
        ax3.plot(self.noise.time_data, self.noise.noise_signal[channel_index].signal)

        if display:
            pl.show(False)
        if save:
            assert isinstance(save, str), "Please provide filename"
            pl.savefig(save+'SNR.pdf') def characterize_population_activity(population_object, parameter_set, analysis_interval, epochs=None,
                                     plot=True, display=True, save=False, color_map="coolwarm", color_subpop=False,
                                     analysis_pars=None):
    """
    Compute all the relevant metrics of recorded activity (spiking and analog signals), providing
    a thorough characterization and quantification of population dynamics

    :return results: dict
    :param population_object: Population or Network object whose activity should be analyzed
    :param parameter_set: complete ParameterSet
    :param analysis_interval: list or tuple with [start_time, stop time] specifying the time interval to analyse
    :param prng: numpy.random object for precise experiment reproduction
    :param epochs:
    :param plot:
    :param display:
    :param save:
    :param color_map:
    :param color_subpop:
    :param analysis_pars:
    :return:
    """
    if analysis_pars is None:
        raise ValueError("Analysis parameters are required for characterizing population activity!")

    ap = analysis_pars
    pars_activity = ap.population_activity
    subpop_names = None

    if isinstance(population_object, na.Population):
        gids = None
        base_population_object = None
    elif isinstance(population_object, na.Network):
        # TODO following 3-4 lines should be checked again
        new_population = population_object.merge_subpopulations(sub_populations=population_object.populations,
                                                                merge_activity=True, store=True)
        population_object.merge_population_activity(start=analysis_interval[0], stop=analysis_interval[1])

        gids = [n.id_list for n in list(sg.iterate_obj_list(population_object.spiking_activity))]
        subpop_names = population_object.population_names

        if not gids:
            gids = [np.array(n.gids) for n in list(sg.iterate_obj_list(population_object.populations))]

        base_population_object 	= population_object
        population_object 		= new_population
    else:
        raise TypeError("Incorrect population object. Must be Population or Network object")

    results = {'spiking_activity': {}, 'analog_activity': {}, 'metadata': {'population_name': population_object.name}}

    ########################################################################################################
    # Spiking activity analysis
    if population_object.spiking_activity:
        spike_list = population_object.spiking_activity
        assert isinstance(spike_list, sg.SpikeList), "Spiking activity should be SpikeList object"
        spike_list = spike_list.time_slice(analysis_interval[0], analysis_interval[1])

        results['spiking_activity'].update(compute_spikelist_metrics(spike_list, population_object.name, ap))

        if plot and ap.depth % 2 == 0: # save all data
            vz.plot_isi_data(results['spiking_activity'][population_object.name],
                             data_label=population_object.name, color_map=color_map, location=0,
                             display=display, save=save)
            if has_pyspike:
                vz.plot_synchrony_measures(results['spiking_activity'][population_object.name],
                                           label=population_object.name, time_resolved=pars_activity.time_resolved,
                                           epochs=epochs, display=display, save=save)

        if pars_activity.time_resolved:
            # *** Averaged time-resolved metrics
            results['spiking_activity'][population_object.name].update(
                compute_time_resolved_statistics(spike_list, label=population_object.name,
                                                 time_bin=pars_activity.time_bin, epochs=epochs,
                                                 window_len=pars_activity.window_len, color_map=color_map,
                                                 display=display, plot=plot, save=save))
        if plot:
            results['metadata']['spike_list'] = spike_list

        if color_subpop and subpop_names:
            results['metadata'].update({'sub_population_names': subpop_names, 'sub_population_gids': gids,
                                        'spike_data_file': ''})

        if gids and ap.depth >= 3:
            if len(gids) == 2:
                locations = [-1, 1]
            else:
                locations = [0 for _ in range(len(gids))]

            for indice, name in enumerate(subpop_names):
                if base_population_object.spiking_activity[indice]:
                    sl = base_population_object.spiking_activity[indice]
                    results['spiking_activity'].update(compute_spikelist_metrics(sl, name, ap))

                if plot and ap.depth % 2 == 0: # save all data
                    vz.plot_isi_data(results['spiking_activity'][name], data_label=name, color_map=color_map,
                                     location=locations[indice], display=display, save=save)
                    if has_pyspike:
                        vz.plot_synchrony_measures(results['spiking_activity'][name], label=name,
                                                   time_resolved=pars_activity.time_resolved,
                                                   display=display, save=save)
                if pars_activity.time_resolved:
                    # *** Averaged time-resolved metrics
                    results['spiking_activity'][name].update(compute_time_resolved_statistics(spike_list,
                                                                                              label=population_object.name, time_bin=pars_activity.time_bin,
                                                                                              epochs=epochs, color_map=color_map, display=display,
                                                                                              plot=plot, save=save, window_len=pars_activity.window_len))
    else:
        logger.warning("Warning, the network is not spiking or no spike recording devices were attached.")

    # Analog activity analysis
    if population_object.analog_activity and base_population_object is not None:
        results['analog_activity'] = {}
        for pop_n, pop in enumerate(base_population_object.populations):
            if bool(pop.analog_activity):
                results['analog_activity'].update({pop.name: {}})
                pop_idx = parameter_set.net_pars.pop_names.index(pop.name)
                if parameter_set.net_pars.analog_device_pars[pop_idx] is None:
                    break
                variable_names = list(np.copy(parameter_set.net_pars.analog_device_pars[pop_idx]['record_from']))

                results['analog_activity'][pop.name].update(compute_analog_stats(pop, parameter_set, variable_names,
                                                                                 analysis_interval, plot))
    if plot:
        vz.plot_state_analysis(parameter_set, results, summary_only=bool(ap.depth % 2 != 0),
                               start=analysis_interval[0], stop=analysis_interval[1],
                               display=display, save=save)
    return results
from modules.parameters import ParameterSet, ParameterSpace, extract_nestvalid_dict
from modules.input_architect import EncodingLayer, InputNoise
from modules.net_architect import Network
from modules.io import set_storage_locations
from modules.signals import iterate_obj_list
from modules.visualization import set_global_rcParams, InputPlots
from modules.analysis import characterize_population_activity
import cPickle as pickle
import numpy as np
import nest
import time
__author__ = 'duarte'


def run(parameter_set, plot=False, display=False, save=True):
    """

    :param parameter_set:
    :param plot:
    :param display:
    :param save:
    :return:
    """
    if not isinstance(parameter_set, ParameterSet):
        if isinstance(parameter_set, basestring) or isinstance(parameter_set, dict):
            parameter_set = ParameterSet(parameter_set)
        else:
            raise TypeError("parameter_set must be ParameterSet, string with full path to parameter file or "
                            "dictionary")

    # ######################################################################################################################
    # Setup extra variables and parameters
    # ======================================================================================================================
    if plot:
        set_global_rcParams(parameter_set.kernel_pars['mpl_path'])
    paths = set_storage_locations(parameter_set, save)

    np.random.seed(parameter_set.kernel_pars['np_seed'])
    results = dict()

    # ######################################################################################################################
    # Set kernel and simulation parameters
    # ======================================================================================================================
    print('\nRuning ParameterSet {0}'.format(parameter_set.label))
    nest.ResetKernel()
    nest.set_verbosity('M_WARNING')
    nest.SetKernelStatus(extract_nestvalid_dict(parameter_set.kernel_pars.as_dict(), param_type='kernel'))

    # ######################################################################################################################
    # Build network
    # ======================================================================================================================
    net = Network(parameter_set.net_pars)

    # ######################################################################################################################
    # Randomize initial variable values
    # ======================================================================================================================
    for idx, n in enumerate(list(iterate_obj_list(net.populations))):
        if hasattr(parameter_set.net_pars, "randomize_neuron_pars"):
            randomize = parameter_set.net_pars.randomize_neuron_pars[idx]
            for k, v in randomize.items():
                n.randomize_initial_states(k, randomization_function=v[0], **v[1])

    # ######################################################################################################################
    # Build and connect input
    # ======================================================================================================================
    # Poisson input

    enc_layer = EncodingLayer(parameter_set.encoding_pars)