From 4add86ee1f365a36a6c4fdc8b48f726e6186305c Mon Sep 17 00:00:00 2001 From: Thomas Limbacher Date: Fri, 10 Jul 2020 15:07:41 +0200 Subject: [PATCH] Add utils --- utils/eval_rewiring.py | 104 +++++++ utils/eval_rewiring_squential.py | 104 +++++++ utils/graph_rewiring_v1.py | 453 +++++++++++++++++++++++++++++ utils/graph_rewiring_v2.py | 420 ++++++++++++++++++++++++++ utils/graph_rewiring_v3.py | 440 ++++++++++++++++++++++++++++ utils/plot_sensitivity_analysis.py | 113 +++++++ utils/stats_rewiring.py | 438 ++++++++++++++++++++++++++++ 7 files changed, 2072 insertions(+) create mode 100644 utils/eval_rewiring.py create mode 100644 utils/eval_rewiring_squential.py create mode 100644 utils/graph_rewiring_v1.py create mode 100644 utils/graph_rewiring_v2.py create mode 100644 utils/graph_rewiring_v3.py create mode 100644 utils/plot_sensitivity_analysis.py create mode 100644 utils/stats_rewiring.py diff --git a/utils/eval_rewiring.py b/utils/eval_rewiring.py new file mode 100644 index 0000000..dc00c6c --- /dev/null +++ b/utils/eval_rewiring.py @@ -0,0 +1,104 @@ +# /usr/bin/env python + +import os + +import utils as utils +from spaghetti import spaghetti_global as spaghetti +from spaghetti.mc_lif_group import McLifGroup +from spaghetti.poisson_pattern_group import PoissonPatternGroup +from spaghetti.rewiring_connection import RewiringConnection +from spaghetti.spike_monitor import SpikeMonitor +from spaghetti.voltage_monitor import VoltageMonitor + + +def main(simulation_time, experiment, sim_date, assemblies, + assembly_neurons_idc, config, output_directory): + + input_params = config["input_parameters"] + connection_params = config["connection_parameters"] + neuron_params = config["neuron_parameters"] + simulation_time_per_pattern = (input_params["pattern_duration"] + + input_params["pattern_delay"]) + + # Initialize the simulation environment. + spaghetti.spaghetti_init(directory=output_directory) + + # Set the random seed. + spaghetti.kernel.set_master_seed(config["master_seed"]) + + # Create input neurons. + inp = PoissonPatternGroup(input_params["num_inputs"], input_params["rate"], + input_params["rate_bg"], params=input_params) + + # Create the neuron. + neuron = McLifGroup(1, neuron_params["num_branches"], neuron_params) + + # Connect input to neuron. + conn = RewiringConnection(inp, neuron, neuron.branch.syn_current, + connection_params) + conn.learn = False + + # Create some monitors which will record the simulation data. + SpikeMonitor(neuron, spaghetti.kernel.fn("test_output", "ras")) + SpikeMonitor(inp, spaghetti.kernel.fn("test_input", "ras")) + VoltageMonitor(neuron.soma, 0, spaghetti.kernel.fn("test_soma", "mem"), + paste_spikes=True, paste_spikes_from=neuron) + for i in range(neuron_params["num_branches"]): + VoltageMonitor(neuron.branch, (i, 0), + spaghetti.kernel.fn("test_branch", "mem", i)) + + # Now simulate the model. + conn.set_weights(weights_pre) + for assembly in assemblies: + inp.set_assemblies(assembly) + inp.set_assembly_neurons_idc(assembly_neurons_idc) + spaghetti.kernel.run_chunk( + simulation_time_per_pattern, 0, simulation_time) + + conn.set_weights(weights_post) + for assembly in assemblies: + inp.set_assemblies(assembly) + inp.set_assembly_neurons_idc(assembly_neurons_idc) + spaghetti.kernel.run_chunk( + simulation_time_per_pattern, 0, simulation_time) + + spaghetti.kernel.run_chunk(input_params["pattern_delay"], 0, + simulation_time) + + +if __name__ == '__main__': + import numpy as np + + experiment = "rewiring_ex1" + sim_date = "181228_163308-19" + master_seed = 13 + simulation_time = 3.2 + assemblies = [0, 1, 2] + assembly_neurons_idc = [i * 40 + np.arange(40) for i in assemblies[0:2]] + assembly_neurons_idc += [4 * 40 + np.arange(40)] + + input_directory = os.path.join("results", experiment, sim_date, "data") + + # Load the configuration file. + config = utils.load_configuration(os.path.join( + input_directory, "..", "config_" + experiment + ".yaml")) + num_branches = config["neuron_parameters"]["num_branches"] + config["master_seed"] = master_seed + if experiment == "rewiring_ex3": + assembly_neurons_idc = [] + assembly_neurons_idc = np.loadtxt(os.path.join(input_directory, + "assembly_neurons_idc"), + dtype=np.int) + else: + np.save(os.path.join(input_directory, "assembly_neurons_idc"), + assembly_neurons_idc) + + # Load the weights. + header_lenght = 3 + with open(os.path.join(input_directory, "weights.0.dat"), "rb") as f: + lines = f.readlines() + weights_pre = np.loadtxt(lines[header_lenght:num_branches + header_lenght]) + weights_post = np.loadtxt(lines[-num_branches:]) + + main(simulation_time, experiment, sim_date, assemblies, + assembly_neurons_idc, config, input_directory) diff --git a/utils/eval_rewiring_squential.py b/utils/eval_rewiring_squential.py new file mode 100644 index 0000000..8391d27 --- /dev/null +++ b/utils/eval_rewiring_squential.py @@ -0,0 +1,104 @@ +# /usr/bin/env python + +import os + +import utils as utils +from spaghetti import spaghetti_global as spaghetti +from spaghetti.mc_lif_group import McLifGroup +from spaghetti.poisson_pattern_group import PoissonPatternGroup +from spaghetti.rewiring_connection import RewiringConnection +from spaghetti.spike_monitor import SpikeMonitor +from spaghetti.voltage_monitor import VoltageMonitor + + +def main(experiment, sim_date, assemblies, weights, + assembly_neurons_idc, config, output_directory): + + input_params = config["input_parameters"] + connection_params = config["connection_parameters"] + neuron_params = config["neuron_parameters"] + simulation_time_per_pattern = input_params["pattern_duration"] + input_params["pattern_delay"] + + # Initialize the simulation environment. + spaghetti.spaghetti_init(directory=output_directory) + + # Set the random seed. + spaghetti.kernel.set_master_seed(config["master_seed"]) + + # Create input neurons. + inp = PoissonPatternGroup(input_params["num_inputs"], input_params["rate"], input_params["rate_bg"], + params=input_params) + + # Create the neuron. + neuron = McLifGroup(1, neuron_params["num_branches"], neuron_params) + + # Connect input to neuron. + conn = RewiringConnection(inp, neuron, neuron.branch.syn_current, connection_params) + conn.learn = False + + # Create some monitors which will record the simulation data. + SpikeMonitor(neuron, spaghetti.kernel.fn("test_output", "ras")) + SpikeMonitor(inp, spaghetti.kernel.fn("test_input", "ras")) + VoltageMonitor(neuron.soma, 0, spaghetti.kernel.fn("test_soma", "mem"), + paste_spikes=True, paste_spikes_from=neuron) + for i in range(neuron_params["num_branches"]): + VoltageMonitor(neuron.branch, (i, 0), spaghetti.kernel.fn("test_branch", "mem", i)) + + # Now simulate the model. + for w in weights: + conn.set_weights(w) + for assembly in assemblies: + inp.set_assemblies(assembly) + inp.set_assembly_neurons_idc(assembly_neurons_idc) + spaghetti.kernel.run_chunk( + simulation_time_per_pattern, 0, + len(assemblies) * len(weights) * simulation_time_per_pattern + + input_params["pattern_delay"]) + + spaghetti.kernel.run_chunk(input_params["pattern_delay"], 0, + len(assemblies) * len(weights) * + simulation_time_per_pattern + + input_params["pattern_delay"]) + + +if __name__ == '__main__': + import numpy as np + + experiment = "rewiring_ex3" + sim_date = "191204_132602/17" + master_seed = 4 + np.random.seed(master_seed) + assemblies = [0, 3, 6] + assembly_neurons_idc = [i * 40 + np.arange(40) for i in range(8)] + + # t = 1000s + # idc1 = 5424 + # idc2 = 21674 + # idc3 = 43345 + # t = 2000s + idc1 = 10845 + idc2 = 43345 + idc3 = 75845 + # t = 5000s + # idc1 = 27095 + # idc2 = 108345 + # idc3 = 189595 + + input_directory = os.path.join("results", experiment, sim_date, "data") + + # Load the configuration file. + config = utils.load_configuration(os.path.join( + input_directory, "..", "config_" + experiment + ".yaml")) + num_branches = config["neuron_parameters"]["num_branches"] + config["master_seed"] = master_seed + + # Load the weights. + with open(os.path.join(input_directory, "weights.0.dat"), "rb") as f: + lines = f.readlines() + + weights1 = np.loadtxt(lines[idc1:idc1 + num_branches]) + weights2 = np.loadtxt(lines[idc2:idc2 + num_branches]) + weights3 = np.loadtxt(lines[idc3:idc3 + num_branches]) + weights = [weights1, weights2, weights3] + + main(experiment, sim_date, assemblies, weights, assembly_neurons_idc, config, input_directory) diff --git a/utils/graph_rewiring_v1.py b/utils/graph_rewiring_v1.py new file mode 100644 index 0000000..bda032d --- /dev/null +++ b/utils/graph_rewiring_v1.py @@ -0,0 +1,453 @@ +#!/usr/bin/env python + +import os +import subprocess + +import ruamel_yaml as yaml + +import configure_seaborn as cs +import matplotlib.lines as mlines +import matplotlib.pyplot as plt +import networkx as nx +import numpy as np +import seaborn as sns +import utils as utils +from matplotlib.gridspec import GridSpec +from matplotlib.patches import Ellipse, FancyArrowPatch, Polygon + +sns.set(context='paper', style='ticks', rc=cs.rc_params) + + +def load_configuration(filepath): + with open(filepath, "r") as f: + config = yaml.safe_load(f) + + return config + + +def draw_connections(G, pos, node_colors, ax, edge_weights=None): + + alpha = 1.0 + for n in G: + c = Ellipse(pos[n], width=0.015, height=0.015, + alpha=alpha, color=node_colors[n], clip_on=False) + ax.add_patch(c) + G.nodes[n]["patch"] = c + x, y = pos[n] + seen = {} + alpha = 1.0 # 0.8 + for (u, v, d) in G.edges(data=True): + n1 = G.nodes[u]["patch"] + n2 = G.nodes[v]["patch"] + rad = 0.1 + if (u, v) in seen: + rad = seen.get((u, v)) + rad = (rad + np.sign(rad) * 0.1) * -1 + color = node_colors[u] + e = FancyArrowPatch(n1.center, + n2.center, + patchA=n1, + patchB=n2, + shrinkA=0, + shrinkB=0, + arrowstyle='-', + linewidth=0.5, + connectionstyle="arc3, rad=%s" % rad, + mutation_scale=10.0, + alpha=alpha, + color=color, + clip_on=False) + ax.add_patch(e) + seen[(u, v)] = rad + + +def draw_assemblies(G, assemblies, colors): + node_colors = [] + for i, assembly in enumerate(assemblies): + G.add_nodes_from(assembly) + node_colors += [colors[i]] * len(assembly) + + return node_colors + + +def draw_neuron(G, ax, branch_nodes, center_assemblies): + x_c = center_assemblies + x_branch_end = [x_c - 0.75, x_c - 0.2, x_c + 0.6] + + b1 = lambda x: -0.2666667 * x - 0.573 # noqa + b2 = lambda x: -3.18 * x - 1.35 # noqa + b3 = lambda x: 0.666667 * x - 0.423 # noqa + b_fun = [b1, b2, b3] + + # Branch nodes + pos_branch_nodes = [] + for xe, bf, branch_node in zip(x_branch_end, b_fun, branch_nodes): + dx = (x_c - xe) / 4 + pos_branch_nodes += [(x_c - dx, bf(x_c - dx)), + (x_c - 2 * dx, bf(x_c - 2 * dx)), + (x_c - 3 * dx, bf(x_c - 3 * dx))] + G.add_nodes_from(branch_node) + node_colors = ["w"] * len(pos_branch_nodes) + + # Neuron + xy = np.array([[x_c, -1.6], [x_c + 0.07, -1 + y_offset], + [x_c - 0.07, -1 + y_offset]]) + nrn = Polygon(xy, clip_on=False, fill=False, color="k", lw=1) + ax.add_patch(nrn) + + trunk = mlines.Line2D([x_c, x_c], [-1.6, -0.5], clip_on=False, color="k", linewidth=1) + ax.add_line(trunk) + + branch1 = mlines.Line2D([x_c - 0.01, x_c - 0.75], [-0.5, -0.3], clip_on=False, color="k", linewidth=1) + ax.add_line(branch1) + + branch2 = mlines.Line2D([x_c - 0.003, x_c - 0.2], [-0.486, 0.15], clip_on=False, color="k", linewidth=1) + ax.add_line(branch2) + + branch3 = mlines.Line2D([x_c + 0.01, x_c + 0.6], [-0.6, -0.2], clip_on=False, color="k", linewidth=1) + ax.add_line(branch3) + + return pos_branch_nodes, node_colors + + +def add_connections(experiment, weights, assemblies, assembly_idc, assembly_map, idc_other_assemblies, + num_neurons_per_assembly, idc_branch_nodes, min_plot_weight): + conn = [] + for i, w in enumerate(weights): + nrns = np.where(w > min_plot_weight)[0] + map_idx_a_idx_nb = {} + last_idx_nb = 0 + for nrn in nrns: + idx_a = np.random.choice( + map_neuron_id_to_assembly_id(nrn, assembly_idc)) + if experiment == "rewiring_ex3": + if idx_a in map_idx_a_idx_nb: + idx_nb = map_idx_a_idx_nb[idx_a] + else: + if last_idx_nb == 0: + idx_nb = 2 + last_idx_nb = 2 + else: + idx_nb = 0 + last_idx_nb = 0 + map_idx_a_idx_nb[idx_a] = idx_nb + else: + idx_nb = np.random.choice(idc_branch_nodes) + idx_nrn = np.random.choice(num_neurons_per_assembly) + if idx_a not in assembly_map.keys(): + idx_a = np.random.choice(idc_other_assemblies) + conn.append((assemblies[idx_a][idx_nrn], branch_nodes[i][idx_nb])) + else: + conn.append((assemblies[assembly_map[idx_a]][idx_nrn], branch_nodes[i][idx_nb])) + G.add_edges_from(conn) + + +def plot_input_spikes(input_spike_times, input_size, ax, xlim, pattern_labels, pattern_colors): + idc = np.where((input_spike_times[:, 0] >= xlim[0]) & (input_spike_times[:, 0] <= xlim[1])) + + ax.scatter(input_spike_times[idc, 0], input_spike_times[idc, 1], s=1.0, color="k", edgecolor="none") + + ax.set_xlim(xlim) + ax.set_ylim([None, input_size + 20]) + ax.set_xticks([]) + ax.set_yticks([1, input_size]) + ax.set_yticklabels("%d" % f for f in ax.get_yticks()) + ax.set_ylabel(r"Input", labelpad=3.2) + + for p, (pl, pc) in enumerate(zip(pattern_labels, pattern_colors)): + ax.plot([xlim[0] + p * 0.5 + 0.2, xlim[0] + (p + 1) * 0.5], + [355, 355], color=pc, alpha=1.0, linestyle='-', + linewidth=1, clip_on=False) + ax.text(xlim[0] + (p + 1) * 0.5 - 0.15, 400, pl, ha="center", color=pc, alpha=1.0) + + +def plot_soma_potential(mem_soma, ax, xlim): + idc = np.where((mem_soma[:, 0] >= xlim[0]) & (mem_soma[:, 0] <= xlim[1])) + ax.plot(mem_soma[idc][:, 0], mem_soma[idc][:, 1], color="k", linewidth=0.6) + + ax.set_ylabel(r"$V^{\mathrm{soma}}$ [mV]") + ax.set_xlabel(r"$t$ [s]") + ax.set_ylim([None, -25]) + ax.set_yticks([-70, -25]) + ax.set_yticklabels("%d" % f for f in [-70, -25]) + ax.set_xlim(xlim) + xticks = np.linspace(xlim[0], xlim[1], 3) + ax.set_xticks(xticks) + ax.set_xticklabels("%.1f" % (f - xlim[0]) for f in xticks) + + +def plot_branch_potential(mem_branch, branch_id, ax, xlim): + idc = np.where((mem_branch[:, 0] >= xlim[0]) & (mem_branch[:, 0] <= xlim[1])) + ax.plot(mem_branch[idc][:, 0] - 0 * xlim[0], mem_branch[idc][:, 1], color="k", linewidth=0.6) + ax.set_xlim(xlim) + ax.set_ylim([-72, -25]) + xticks = ax.get_xticks() + ax.set_xticks([]) + ax.set_yticks([]) + ax.set_ylabel(r"$V^{\mathrm{b}}_{" + str(branch_id) + "}$", rotation=0, va="center") + if branch_id >= 10: + ax.yaxis.set_label_coords(-0.044, 0.75) + else: + ax.yaxis.set_label_coords(-0.05, 0.75) + for spine in ["left", "bottom"]: + ax.spines[spine].set_visible(False) + + return xticks + + +def plot_scale(xlim): + line = mlines.Line2D([xlim[0] + 0.006, xlim[0] + 0.306], [-73.5, -73.5], clip_on=False, color=".15", + linewidth=0.7) + ax.add_line(line) + ax.text(np.mean([xlim[0] + 0.006, xlim[0] + 0.306]), -86.5, r"0.3 s", ha="center", fontsize=8) + + line = mlines.Line2D([xlim[0] - 0.03, xlim[0] - 0.03], [-68.9, -48.9], clip_on=False, color=".15", + linewidth=0.7) + ax.add_line(line) + ax.text(xlim[0] - 0.365, -64, r"20 mV", fontsize=8) + + +def map_neuron_id(gid, old_min=0, old_max=319, new_min=0, new_max=35): + old_range = (old_max - old_min) + new_range = (new_max - new_min) + return int(((((gid - old_min) * new_range) / old_range) + new_min)) + + +def map_neuron_id_to_assembly_id(gid, assembly_idc): + if not any(gid in x for x in assembly_idc): + return [-1] + + return np.where(assembly_idc == gid)[0] + + +# ------------------------------------------------------------------------------ +experiment = "rewiring_ex1" +sim_date = "191125_135334/1" +patterns = [0, 1, 4] +patterns_graph = [0, 1, 2] +branches = [1, 5, 11] +xlims = [[0.0, 1.7], [1.5, 3.2]] + +cs.set_figure_size(84 + 9, 87 + 8) + +# ------------------------------------------------------------------------------ +# Directory of simulation results and log files. +if experiment == "rewiring_ex4": + input_directory = os.path.join("results", experiment, "4", sim_date, "data") +else: + input_directory = os.path.join("results", experiment, sim_date, "data") + +# Directory for plots. +plots_directory = os.path.join(input_directory, "..", "plots") +if not os.path.exists(plots_directory): + os.makedirs(plots_directory) + +# Colors of patters. +c = sns.color_palette().as_hex() +c[4], c[5], c[6], c[7], c[8], c[9] = c[4], c[8], c[5], c[6], c[8], c[7] + +colors = [c[9]] +colors += [c[i] for i, p in enumerate(patterns)] +colors += [c[9]] +pattern_labels = [[r"$\mathrm{A}_%d$" % (p + 1) for p in patterns], + [r"$\mathrm{A}_%d$" % (p + 1) for p in patterns]] +pattern_colors = [[c[p] for p in patterns], [c[p] for p in patterns]] +assembly_map = {p: i + 1 for i, p in enumerate(patterns_graph)} +branches.sort() + +np.random.seed(0) +num_rows = 3 +node_colors = [] +num_branches = 3 +num_assemblies = 5 +min_plot_weight = 1.0 +idc_other_assemblies = [0, 4] +num_assemblies_real = 3 +num_branch_nodes = 3 +num_neurons_per_row = 7 +num_neurons_per_assembly = num_rows * num_neurons_per_row +pos_assemblies = [] + +x_offset = -0.54 +y_offset = -0.79 +for i in range(num_rows): + for x in np.linspace(-1 + x_offset, 1, 35): + pos_assemblies.append((x, 1.0 - i * 0.1)) +np.random.shuffle(pos_assemblies) + +assemblies = np.split(np.arange(num_neurons_per_assembly * num_assemblies), num_assemblies) + +branch1_nodes = np.max(assemblies) + 1 + np.arange(num_branch_nodes) +branch2_nodes = np.max(branch1_nodes) + 1 + np.arange(num_branch_nodes) +branch3_nodes = np.max(branch2_nodes) + 1 + np.arange(num_branch_nodes) +branch_nodes = [branch1_nodes, branch2_nodes, branch3_nodes] + +if experiment == "rewiring_ex5": + idc_branch_nodes = [0, 2] +else: + idc_branch_nodes = range(num_branch_nodes) + +# Load the configuration file. +config = utils.load_configuration(os.path.join( + input_directory, "..", "config_" + experiment + ".yaml")) +sim_simulation_time = config["simulation_time"] +sim_w_max = config["connection_parameters"]["w_max"] +sim_input_size = config["input_parameters"]["num_inputs"] +sim_num_assemblies = config["input_parameters"]["num_assemblies"] +sim_assembly_size = config["input_parameters"]["assembly_size"] +sim_num_branches = config["neuron_parameters"]["num_branches"] +sim_sampling_interval_weights = config["sampling_interval_weights"] +input_size = config["input_parameters"]["num_inputs"] + +if experiment == "rewiring_ex5": + assembly_idc = [] + assembly_idc = np.loadtxt( + os.path.join(input_directory, "assembly_neurons_idc"), dtype=np.int) + # for assembly_idx_low in np.sort(assembly_idc_low): + # assembly_idc.append(np.arange(assembly_idx_low, assembly_idx_low + + # sim_assembly_size)) +else: + assembly_idc = np.split(np.arange(sim_num_assemblies * sim_assembly_size), + sim_num_assemblies) + +# Load the simulation results. +mem_branch = [] +for b in branches: + mem_branch.append(np.loadtxt(os.path.join( + input_directory, "test_branch" + str(b) + ".0.mem"))) +mem_soma = np.loadtxt(os.path.join(input_directory, 'test_soma.0.mem')) +input_spike_times = np.loadtxt(os.path.join( + input_directory, 'test_input.0.ras')) +input_spike_times[:, 1] = np.random.randint(0, sim_input_size, + len(input_spike_times)) + +header_lenght = 3 +with open(os.path.join(input_directory, "weights.0.dat"), "rb") as f: + lines = f.readlines() +weights_pre = np.loadtxt(lines[header_lenght:sim_num_branches + header_lenght]) +weights_train_start = [weights_pre[b] for b in branches] +weights_post = np.loadtxt(lines[-sim_num_branches:]) +weights_train_end = [weights_post[b] for b in branches] + +# Before training. +# ------------------------------------------------------------------------------ +fig = plt.figure() +gs = GridSpec(6, 2) + +# Input spikes. +ax = plt.subplot(gs[0, :]) +plot_input_spikes(input_spike_times, input_size, ax, xlims[0], pattern_labels[0], pattern_colors[0]) + +# Create graph. +node_colors = [] +G = nx.MultiGraph() +ax = plt.subplot(gs[1:-1, 0]) + +# Draw the assemblies. +node_colors += draw_assemblies(G, assemblies, colors) + +# Draw the neuron. +xc = np.mean(pos_assemblies, axis=0)[0] +pos_branch_nodes, nc = draw_neuron(G, ax, branch_nodes, xc) +node_colors += nc + +# Add connections to graph. +add_connections(experiment, weights_train_start, assemblies, assembly_idc, + assembly_map, idc_other_assemblies, num_neurons_per_assembly, + idc_branch_nodes, min_plot_weight) + +# Draw connections. +draw_connections(G, (pos_assemblies + pos_branch_nodes), node_colors, ax) + +plt.axis('off') +ax.set_xlim([-1, 1]) +ax.set_ylim([-1 + y_offset, 1]) + +# Branch 3 potential. +ax = plt.subplot(gs[2, 1:]) +xticks = plot_branch_potential(mem_branch[2], branches[2] + 1, ax, xlims[0]) + +# Scale +plot_scale(xlims[0]) + +# Branch 2 potential. +ax = plt.subplot(gs[3, 1:]) +xticks = plot_branch_potential(mem_branch[1], branches[1] + 1, ax, xlims[0]) + +# Branch 1 potential. +ax = plt.subplot(gs[4, 1:]) +xticks = plot_branch_potential(mem_branch[0], branches[0] + 1, ax, xlims[0]) + +# Soma potential. +ax = plt.subplot(gs[-1, :]) +plot_soma_potential(mem_soma, ax, xlims[0]) + +plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=0.4) + +fname = os.path.join(plots_directory, "graph-before") +fig.savefig(fname + ".pdf", pad_inches=0.01) +subprocess.call(["pdftops", "-eps", fname + ".pdf", fname + ".eps"]) + +plt.close(fig) + + +# After training. +# ------------------------------------------------------------------------------ +fig = plt.figure() +gs = GridSpec(6, 2) + +# Input spikes. +ax = plt.subplot(gs[0, :]) +plot_input_spikes(input_spike_times, input_size, ax, xlims[1], pattern_labels[1], pattern_colors[1]) + +# Create graph. +node_colors = [] +G = nx.MultiGraph() +ax = plt.subplot(gs[1:-1, 0]) + +# Draw the assemblies. +node_colors += draw_assemblies(G, assemblies, colors) + +# Draw the neuron. +xc = np.mean(pos_assemblies, axis=0)[0] +pos_branch_nodes, nc = draw_neuron(G, ax, branch_nodes, xc) +node_colors += nc + +# Add connections to graph. +add_connections(experiment, weights_train_end, assemblies, assembly_idc, + assembly_map, idc_other_assemblies, num_neurons_per_assembly, + idc_branch_nodes, min_plot_weight) + +# Draw connections. +draw_connections(G, (pos_assemblies + pos_branch_nodes), node_colors, ax) + +plt.axis('off') +ax.set_xlim([-1, 1]) +ax.set_ylim([-1 + y_offset, 1]) + +# Branch 3 potential. +ax = plt.subplot(gs[2, 1:]) +xticks = plot_branch_potential(mem_branch[2], branches[2] + 1, ax, xlims[1]) + +# Scale +plot_scale(xlims[1]) + +# Branch 2 potential. +ax = plt.subplot(gs[3, 1:]) +xticks = plot_branch_potential(mem_branch[1], branches[1] + 1, ax, xlims[1]) + +# Branch 1 potential. +ax = plt.subplot(gs[4, 1:]) +xticks = plot_branch_potential(mem_branch[0], branches[0] + 1, ax, xlims[1]) + +# Soma potential. +ax = plt.subplot(gs[-1, :]) +plot_soma_potential(mem_soma, ax, xlims[1]) + +plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=None, hspace=0.4) + +fname = os.path.join(plots_directory, "graph-after") +fig.savefig(fname + ".pdf", pad_inches=0.01) +subprocess.call(["pdftops", "-eps", fname + ".pdf", fname + ".eps"]) + +plt.close(fig) diff --git a/utils/graph_rewiring_v2.py b/utils/graph_rewiring_v2.py new file mode 100644 index 0000000..01a72a9 --- /dev/null +++ b/utils/graph_rewiring_v2.py @@ -0,0 +1,420 @@ +#!/usr/bin/env python + +import os +import subprocess + +import ruamel_yaml as yaml + +import configure_seaborn as cs +import matplotlib.lines as mlines +import matplotlib.pyplot as plt +import networkx as nx +import numpy as np +import seaborn as sns +import utils as utils +from matplotlib.gridspec import GridSpec +from matplotlib.patches import Ellipse, FancyArrowPatch, Polygon + +sns.set(context='paper', style='ticks', rc=cs.rc_params) + + +def load_configuration(filepath): + with open(filepath, "r") as f: + config = yaml.safe_load(f) + + return config + + +def draw_connections(G, pos, node_colors, ax, edge_weights=None): + + alpha = 1.0 + for n in G: + c = Ellipse(pos[n], width=0.015, height=0.015, + alpha=alpha, color=node_colors[n], clip_on=False) + ax.add_patch(c) + G.nodes[n]["patch"] = c + x, y = pos[n] + seen = {} + alpha = 1.0 # 0.8 + for (u, v, d) in G.edges(data=True): + n1 = G.nodes[u]["patch"] + n2 = G.nodes[v]["patch"] + rad = 0.1 + if (u, v) in seen: + rad = seen.get((u, v)) + rad = (rad + np.sign(rad) * 0.1) * -1 + color = node_colors[u] + e = FancyArrowPatch(n1.center, + n2.center, + patchA=n1, + patchB=n2, + shrinkA=0, + shrinkB=0, + arrowstyle='-', + linewidth=0.5, + connectionstyle="arc3, rad=%s" % rad, + mutation_scale=10.0, + alpha=alpha, + color=color, + clip_on=False) + ax.add_patch(e) + seen[(u, v)] = rad + + +def draw_assemblies(G, assemblies, colors): + node_colors = [] + for i, assembly in enumerate(assemblies): + G.add_nodes_from(assembly) + node_colors += [colors[i]] * len(assembly) + + return node_colors + + +def draw_neuron(G, ax, branch_nodes, center_assemblies): + x_c = center_assemblies + x_branch_end = [x_c - 0.75, x_c - 0.2, x_c + 0.6, x_c + 1.0] + + def b1(x): return -0.2666667 * x - 0.573 # noqa + def b2(x): return -3.18 * x - 1.35 # noqa + def b3(x): return 0.666667 * x - 0.423 # noqa + def b4(x): return 0.201299 * x - 0.516 # noqa + b_fun = [b1, b2, b3, b4] + + # Branch nodes + pos_branch_nodes = [] + for i, (xe, bf, branch_node) in enumerate(zip(x_branch_end, b_fun, branch_nodes)): + dx = (x_c - xe) / 4 + pos_branch_nodes += [(x_c - dx, bf(x_c - dx)), + (x_c - 2 * dx, bf(x_c - 2 * dx)), + (x_c - 3 * dx, bf(x_c - 3 * dx))] + G.add_nodes_from(branch_node) + node_colors = ["w"] * len(pos_branch_nodes) + + # Neuron + xy = np.array([[x_c, -1.6], [x_c + 0.07, -1 + y_offset], + [x_c - 0.07, -1 + y_offset]]) + nrn = Polygon(xy, clip_on=False, fill=False, color="k", lw=1) + ax.add_patch(nrn) + + trunk = mlines.Line2D([x_c, x_c], [-1.6, -0.5], clip_on=False, color="k", linewidth=1) + ax.add_line(trunk) + + branch1 = mlines.Line2D([x_c - 0.01, x_c - 0.75], [-0.5, -0.3], clip_on=False, color="k", linewidth=1) + ax.add_line(branch1) + + branch2 = mlines.Line2D([x_c - 0.003, x_c - 0.2], [-0.486, 0.15], clip_on=False, color="k", linewidth=1) + ax.add_line(branch2) + + branch3 = mlines.Line2D([x_c + 0.01, x_c + 0.7], [-0.6, -0.11], clip_on=False, color="k", linewidth=1) + ax.add_line(branch3) + + branch4 = mlines.Line2D([x_c + 0.08, x_c + 0.85], [-0.555, -0.4], clip_on=False, color="k", linewidth=1) + ax.add_line(branch4) + + return pos_branch_nodes, node_colors + + +def add_connections(experiment, weights, assemblies, assembly_idc, assembly_map, idc_other_assemblies, + num_neurons_per_assembly, idc_branch_nodes, min_plot_weight): + conn = [] + for i, w in enumerate(weights): + nrns = np.where(w > min_plot_weight)[0] + print(len(nrns)) + map_idx_a_idx_nb = {} + last_idx_nb = 0 + for nrn in nrns: + idx_a = np.random.choice(map_neuron_id_to_assembly_id(nrn, assembly_idc)) + if experiment == "rewiring_ex5": + if idx_a in map_idx_a_idx_nb: + idx_nb = map_idx_a_idx_nb[idx_a] + else: + if last_idx_nb == 0: + idx_nb = 2 + last_idx_nb = 2 + else: + idx_nb = 0 + last_idx_nb = 0 + map_idx_a_idx_nb[idx_a] = idx_nb + else: + idx_nb = np.random.choice(idc_branch_nodes) + + idx_nrn = np.random.choice(num_neurons_per_assembly) + if idx_a not in assembly_map.keys(): + idx_a = np.random.choice(idc_other_assemblies) + conn.append((assemblies[idx_a][idx_nrn], branch_nodes[i][idx_nb])) + else: + conn.append((assemblies[assembly_map[idx_a]][idx_nrn], branch_nodes[i][idx_nb])) + + G.add_edges_from(conn) + + +def plot_input_spikes(input_spike_times, input_size, ax, xlim, pattern_labels, pattern_colors): + idc = np.where((input_spike_times[:, 0] >= xlim[0]) & (input_spike_times[:, 0] <= xlim[1])) + + ax.scatter(input_spike_times[idc, 0], input_spike_times[idc, 1], s=1.0, color="k", edgecolor="none") + + ax.set_xlim(xlim) + ax.set_ylim([None, input_size + 20]) + ax.set_xticks([]) + ax.set_yticks([1, input_size]) + ax.set_yticklabels("%d" % f for f in ax.get_yticks()) + ax.set_ylabel("Input", labelpad=3.2) + ax.yaxis.set_tick_params(width=0.8, length=2.5) + + for p, (pl, pc) in enumerate(zip(pattern_labels, pattern_colors)): + ax.plot([xlim[0] + p * 0.5 + 0.2, xlim[0] + (p + 1) * 0.5], + [355, 355], color=pc, alpha=1.0, linestyle='-', + linewidth=1, clip_on=False) + + ax.text(xlim[0] + (p + 1) * 0.5 - 0.15, 400, pl, ha="center", color=pc, alpha=1.0) + + +def plot_soma_potential(mem_soma, ax, xlim): + idc = np.where(mem_soma[:, 0] <= xlim[0] + 0.01) + mem_soma[:, 1][idc] = -70 + + idc = np.where((mem_soma[:, 0] >= xlim[0]) & (mem_soma[:, 0] <= xlim[1])) + ax.plot(mem_soma[idc][:, 0], mem_soma[idc][:, 1], color="k", linewidth=0.6) + + ax.set_ylabel(r"$V^{\mathrm{soma}}$", rotation=0, va="baseline") + ax.yaxis.set_label_coords(0.003, 0.86) + ax.set_yticks([]) + ax.set_ylim([-72, -25]) + ax.set_xlim(xlim) + ax.set_xticks([]) + for spine in ["left", "bottom"]: + ax.spines[spine].set_visible(False) + + +def plot_branch_potential(mem_branch, branch_id, ax, xlim): + idc = np.where((mem_branch[:, 0] >= xlim[0]) & (mem_branch[:, 0] <= xlim[1])) + ax.plot(mem_branch[idc][:, 0] - 0 * xlim[0], mem_branch[idc][:, 1], color="k", linewidth=0.6) + ax.set_xlim(xlim) + ax.set_ylim([-72, -25]) + xticks = ax.get_xticks() + ax.set_xticks([]) + ax.set_yticks([]) + ax.set_ylabel(r"$V^{\mathrm{b}}_{" + str(branch_id) + "}$", rotation=0, va="center") + if branch_id >= 10: + ax.yaxis.set_label_coords(-0.0101, 0.75) + else: + ax.yaxis.set_label_coords(-0.012, 0.75) + for spine in ["left", "bottom"]: + ax.spines[spine].set_visible(False) + + return xticks + + +def plot_scale(xlim): + line = mlines.Line2D([xlim[0] - 0.02, xlim[0] - 0.02], [-68.9, -48.9], clip_on=False, color=".15", + linewidth=0.7) + ax.add_line(line) + ax.text(xlim[0] - 0.19, -64, r"20 mV", fontsize=8) + + +def plot_scale_soma(xlim): + line = mlines.Line2D([xlim[0] + 0.006, xlim[0] + 0.306], [-76.0 - 2, -76.0 - 2], + clip_on=False, color=".15", linewidth=0.7) + ax.add_line(line) + ax.text(np.mean([xlim[0] + 0.006, xlim[0] + 0.306]), -93 - 2, r"0.3 s", ha="center", fontsize=8) + + line = mlines.Line2D([xlim[0] - 0.025, xlim[0] - 0.025], [-68.9 + 3, -44.5 + 3], clip_on=False, + color=".15", linewidth=0.7) + ax.add_line(line) + ax.text(xlim[0] - 0.19, -63 + 3, r"25 mV", fontsize=8) + + +def map_neuron_id(gid, old_min=0, old_max=319, new_min=0, new_max=35): + old_range = (old_max - old_min) + new_range = (new_max - new_min) + return int(((((gid - old_min) * new_range) / old_range) + new_min)) + + +def map_neuron_id_to_assembly_id(gid, assembly_idc): + if not any(gid in x for x in assembly_idc): + return [-1] + + return np.where(assembly_idc == gid)[0] + + +# ------------------------------------------------------------------------------ +experiment = "rewiring_ex2" +sim_date = "191125_114900/22" +patterns = [0, 1, 2, 3, 4, 5] +branches = [1, 2, 3, 4] +branch_names = [2, 3, 4, 5] +xlim = [3.0, 6.2] + +cs.set_figure_size(176 + 39, 68 + 23) + +# ------------------------------------------------------------------------------ +# Directory of simulation results and log files. +if experiment == "rewiring_ex4": + input_directory = os.path.join("results", experiment, "4", sim_date, "data") +else: + input_directory = os.path.join("results", experiment, sim_date, "data") + +# Directory for plots. +plots_directory = os.path.join(input_directory, "..", "plots") +if not os.path.exists(plots_directory): + os.makedirs(plots_directory) + +# Colors of patters. +# c = plt.rcParams['axes.prop_cycle'].by_key()['color'] +# c[6] = c[8] +# c[7] = c[9] +c = sns.color_palette().as_hex() +c2 = sns.hls_palette(8, l=.3, s=.8)[-3:] # noqa + +colors = [c[7]] +colors += [c[p] for i, p in enumerate(patterns[0:3])] +colors += [c[7]] +pattern_labels = [r"$\mathrm{A}_%d$" % (p + 1) for p in patterns[0:3]] +pattern_labels += [r"$\mathrm{R}_{%d}$" % (p + 1) for p in range(3)] +pattern_colors = [c[p] for p in patterns[0:3]] +pattern_colors += [c2[i] for i, p in enumerate(patterns[3:])] + +assembly_map = {p: i + 1 for i, p in enumerate(patterns[0:3])} +# branches.sort() + +np.random.seed(18) # 14, 16, 18 +num_rows = 3 +node_colors = [] +num_branches = 3 +num_assemblies = 5 +min_plot_weight = 0.9 +idc_other_assemblies = [0, 4] +num_assemblies_real = 3 +num_branch_nodes = 3 +num_neurons_per_row = 7 +num_neurons_per_assembly = num_rows * num_neurons_per_row +pos_assemblies = [] + +x_offset = -0.54 +y_offset = -0.79 +for i in range(num_rows): + for x in np.linspace(-1 + x_offset, 1, 35): + pos_assemblies.append((x, 1.0 - i * 0.1)) +np.random.shuffle(pos_assemblies) + +assemblies = np.split(np.arange(num_neurons_per_assembly * num_assemblies), num_assemblies) + +branch1_nodes = np.max(assemblies) + 1 + np.arange(num_branch_nodes) +branch2_nodes = np.max(branch1_nodes) + 1 + np.arange(num_branch_nodes) +branch3_nodes = np.max(branch2_nodes) + 1 + np.arange(num_branch_nodes) +branch4_nodes = np.max(branch3_nodes) + 1 + np.arange(num_branch_nodes) +# branch_nodes = [branch4_nodes, branch2_nodes, branch1_nodes, branch3_nodes] +branch_nodes = [branch1_nodes, branch2_nodes, branch3_nodes, branch4_nodes] + +if experiment == "rewiring_ex3": + idc_branch_nodes = [0, 2] +else: + idc_branch_nodes = range(num_branch_nodes) + +# Load the configuration file. +config = utils.load_configuration(os.path.join(input_directory, "..", "config_" + experiment + ".yaml")) +sim_simulation_time = config["simulation_time"] +sim_w_max = config["connection_parameters"]["w_max"] +sim_input_size = config["input_parameters"]["num_inputs"] +sim_num_assemblies = config["input_parameters"]["num_assemblies"] +sim_assembly_size = config["input_parameters"]["assembly_size"] +sim_num_branches = config["neuron_parameters"]["num_branches"] +sim_sampling_interval_weights = config["sampling_interval_weights"] +input_size = config["input_parameters"]["num_inputs"] + +if experiment == "rewiring_ex3": + assembly_idc = [] + assembly_idc_low = np.loadtxt(os.path.join(input_directory, "assembly_idc_low"), dtype=np.int) + for assembly_idx_low in np.sort(assembly_idc_low): + assembly_idc.append(np.arange(assembly_idx_low, assembly_idx_low + sim_assembly_size)) +else: + assembly_idc = np.split(np.arange(sim_num_assemblies * sim_assembly_size), sim_num_assemblies) + +# Load the simulation results. +mem_branch = [] +for b in branches: + mem_branch.append(np.loadtxt(os.path.join( + input_directory, "test_branch" + str(b) + ".0.mem"))) +mem_soma = np.loadtxt(os.path.join(input_directory, 'test_soma.0.mem')) +input_spike_times = np.loadtxt(os.path.join(input_directory, 'test_input.0.ras')) +input_spike_times[:, 1] = np.random.randint(0, sim_input_size, len(input_spike_times)) + +assembly_neurons_idc_random = np.load(os.path.join( + input_directory, 'assembly_neurons_idc.npy'))[3:] +assembly_neurons_idc = [i * 40 + np.arange(40) for i in range(8)] +overlaps = np.zeros((len(assembly_neurons_idc_random), len(assembly_neurons_idc))) +for i in range(len(assembly_neurons_idc_random)): + for j in range(len(assembly_neurons_idc)): + overlaps[i, j] = np.intersect1d(assembly_neurons_idc_random[i], assembly_neurons_idc[j]).size / 40 + +header_lenght = 3 +with open(os.path.join(input_directory, "weights.0.dat"), "rb") as f: + lines = f.readlines() +weights_post = np.loadtxt(lines[-sim_num_branches:]) +weights_train_end = [weights_post[b] for b in branches] + +# After training. +# ------------------------------------------------------------------------------ +fig = plt.figure() +gs = GridSpec(7, 4) + +# Input spikes. +ax = plt.subplot(gs[0, 1:]) +plot_input_spikes(input_spike_times, input_size, ax, xlim, pattern_labels, pattern_colors) + +# Create graph. +node_colors = [] +G = nx.MultiGraph() +ax = plt.subplot(gs[0:-2, 0]) + +# Draw the assemblies. +node_colors += draw_assemblies(G, assemblies, colors) + +# Draw the neuron. +xc = np.mean(pos_assemblies, axis=0)[0] +pos_branch_nodes, nc = draw_neuron(G, ax, branch_nodes, xc) +node_colors += nc + +# Add connections to graph. +add_connections(experiment, weights_train_end, assemblies, assembly_idc, assembly_map, idc_other_assemblies, + num_neurons_per_assembly, idc_branch_nodes, min_plot_weight) + +# Draw connections. +draw_connections(G, (pos_assemblies + pos_branch_nodes), node_colors, ax) + + +plt.axis('off') +ax.set_xlim([-1, 1]) +ax.set_ylim([-1 + y_offset + 0.00, 1.125]) + +# Branch 4 potential. +ax = plt.subplot(gs[1, 1:]) +xticks = plot_branch_potential(mem_branch[3], branch_names[3], ax, xlim) + +# Scale +# plot_scale(xlim) + +# Branch 3 potential. +ax = plt.subplot(gs[2, 1:]) +xticks = plot_branch_potential(mem_branch[2], branch_names[2], ax, xlim) + +# Branch 2 potential. +ax = plt.subplot(gs[3, 1:]) +xticks = plot_branch_potential(mem_branch[1], branch_names[1], ax, xlim) + +# Branch 1 potential. +ax = plt.subplot(gs[4, 1:]) +xticks = plot_branch_potential(mem_branch[0], branch_names[0], ax, xlim) + +# Soma potential. +ax = plt.subplot(gs[5, 1:]) +plot_soma_potential(mem_soma, ax, xlim) + +plot_scale_soma(xlim) + +plt.subplots_adjust(left=None, bottom=None, right=None, top=None, wspace=0.4, hspace=0.4) + +fname = os.path.join(plots_directory, "graph-v2") +fig.savefig(fname + ".pdf", pad_inches=0.01) +subprocess.call(["pdftops", "-eps", fname + ".pdf", fname + ".eps"]) +plt.close(fig) diff --git a/utils/graph_rewiring_v3.py b/utils/graph_rewiring_v3.py new file mode 100644 index 0000000..6e3f812 --- /dev/null +++ b/utils/graph_rewiring_v3.py @@ -0,0 +1,440 @@ +#!/usr/bin/env python + +import os +import subprocess + +import ruamel_yaml as yaml + +import configure_seaborn as cs +import matplotlib.lines as mlines +import matplotlib.pyplot as plt +import networkx as nx +import numpy as np +import seaborn as sns +import utils as utils +from matplotlib.gridspec import GridSpec +from matplotlib.patches import Ellipse, FancyArrowPatch, Polygon + +sns.set(context='paper', style='ticks', rc=cs.rc_params) + + +def load_configuration(filepath): + with open(filepath, "r") as f: + config = yaml.safe_load(f) + + return config + + +def draw_connections(G, pos, node_colors, ax, edge_weights=None): + + alpha = 1.0 + for n in G: + c = Ellipse(pos[n], width=0.013, height=0.013, + alpha=alpha, color=node_colors[n], clip_on=False) + ax.add_patch(c) + G.nodes[n]["patch"] = c + x, y = pos[n] + seen = {} + alpha = 1.0 # 0.8 + for (u, v, d) in G.edges(data=True): + n1 = G.nodes[u]["patch"] + n2 = G.nodes[v]["patch"] + rad = 0.1 + if (u, v) in seen: + rad = seen.get((u, v)) + rad = (rad + np.sign(rad) * 0.1) * -1 + color = node_colors[u] + e = FancyArrowPatch(n1.center, + n2.center, + patchA=n1, + patchB=n2, + shrinkA=0, + shrinkB=0, + arrowstyle='-', + connectionstyle="arc3, rad=%s" % rad, + mutation_scale=10.0, + alpha=alpha, + lw=0.5, + color=color, + clip_on=False) + ax.add_patch(e) + seen[(u, v)] = rad + + +def draw_assemblies(G, assemblies, colors): + node_colors = [] + for i, assembly in enumerate(assemblies): + G.add_nodes_from(assembly) + for i in range(len(assemblies)): + if i == 0: + node_colors += [colors[i]] * 14 + else: + node_colors += [colors[i]] * 13 + + return node_colors + + +def draw_neuron(G, ax, branch_nodes, center_assemblies): + x_c = center_assemblies + x_branch_end = [x_c - 0.75, x_c - 0.2, x_c + 0.6] + + b1 = lambda x: -0.2666667 * x - 0.573 # noqa + b2 = lambda x: -3.18 * x - 1.35 # noqa + b3 = lambda x: 0.666667 * x - 0.423 # noqa + b_fun = [b1, b2, b3] + + # Branch nodes + pos_branch_nodes = [] + for xe, bf, branch_node in zip(x_branch_end, b_fun, branch_nodes): + dx = (x_c - xe) / 4 + pos_branch_nodes += [(x_c - dx, bf(x_c - dx)), + (x_c - 2 * dx, bf(x_c - 2 * dx)), + (x_c - 3 * dx, bf(x_c - 3 * dx))] + G.add_nodes_from(branch_node) + node_colors = ["w"] * len(pos_branch_nodes) + + # Neuron + xy = np.array([[x_c, -1.6], [x_c + 0.07, -1 + y_offset], [x_c - 0.07, -1 + y_offset]]) + nrn = Polygon(xy, clip_on=False, fill=False, color="k", lw=1) + ax.add_patch(nrn) + + # Trunk + trunk = mlines.Line2D([x_c, x_c], [-1.6, -0.5], clip_on=False, color="k", linewidth=1) + ax.add_line(trunk) + + branch1 = mlines.Line2D([x_c - 0.01, x_c - 0.75], [-0.5, -0.3], clip_on=False, color="k", linewidth=1) + ax.add_line(branch1) + + branch2 = mlines.Line2D([x_c - 0.003, x_c - 0.2], [-0.486, 0.15], clip_on=False, color="k", linewidth=1) + ax.add_line(branch2) + + branch3 = mlines.Line2D([x_c + 0.01, x_c + 0.6], [-0.6, -0.2], clip_on=False, color="k", linewidth=1) + ax.add_line(branch3) + + return pos_branch_nodes, node_colors + + +def add_connections(experiment, weights, assemblies, assembly_idc, assembly_map, idc_other_assemblies, + num_neurons_per_assembly, idc_branch_nodes, min_plot_weight): + conn = [] + idc_a = [] + ii = [] + for i, w in enumerate(weights): + nrns = np.where(w > min_plot_weight)[0] + print(len(nrns)) + for nrn in nrns: + idc_a.append(np.random.choice(map_neuron_id_to_assembly_id(nrn, assembly_idc))) + ii.append(i) + for idx_a, i in zip(idc_a, ii): + idx_nb = np.random.choice(idc_branch_nodes) + idx_nrn = np.random.choice(num_neurons_per_assembly[idx_a]) + if idx_a not in assembly_map.keys(): + idx_a = np.random.choice(idc_other_assemblies) + conn.append((assemblies[idx_a][idx_nrn], branch_nodes[i][idx_nb])) + else: + conn.append((assemblies[assembly_map[idx_a]][idx_nrn], branch_nodes[i][idx_nb])) + G.add_edges_from(conn) + + +def plot_soma_potential(mem_soma, ax, xlim, pattern_labels, pattern_colors): + idc = np.where((mem_soma[:, 0] >= xlim[0]) & (mem_soma[:, 0] <= xlim[1])) + ax.plot(mem_soma[idc][:, 0], mem_soma[idc][:, 1], color="k", + clip_on=False, linewidth=0.6) + + for p, (pl, pc) in enumerate(zip(pattern_labels, pattern_colors)): + ax.plot([xlim[0] + p * 0.5 + 0.2, xlim[0] + (p + 1) * 0.5], [-12, -12], color=pc, alpha=1.0, + linestyle='-', linewidth=0.7, clip_on=False) + ax.text(xlim[0] + (p + 1) * 0.5 - 0.15, -5, pl, ha="center", color=pc, alpha=1.0, fontsize=8) + + ax.set_ylabel(r"$V^{\mathrm{soma}}$", rotation=0, va="center") + ax.yaxis.set_label_coords(-0.1, 1.4) + # ax.set_ylim([None, 20]) + ax.set_yticks([]) + ax.set_xticks([]) + ax.set_xlim(xlim) + + for spine in ["top", "right", "left", "bottom"]: + ax.spines[spine].set_visible(False) + + +def plot_scale(xlim): + line = mlines.Line2D([xlim[0] + 0.006, xlim[0] + 0.306], [-79.0, -79.0], clip_on=False, color="0.15", + linewidth=0.7) + ax.add_line(line) + ax.text(np.mean([xlim[0] + 0.006, xlim[0] + 0.306]), -106, r"0.3 s", ha="center", fontsize=8) + + line = mlines.Line2D([xlim[0] - 0.04, xlim[0] - 0.04], [-68.9, -43], clip_on=False, color="0.15", + linewidth=0.7) + ax.add_line(line) + ax.text(xlim[0] - 0.44, -66, r"25 mV", fontsize=8) + + +def map_neuron_id(gid, old_min=0, old_max=319, new_min=0, new_max=35): + old_range = (old_max - old_min) + new_range = (new_max - new_min) + return int(((((gid - old_min) * new_range) / old_range) + new_min)) + + +def map_neuron_id_to_assembly_id(gid, assembly_idc): + if not any(gid in x for x in assembly_idc): + return [-1] + + return np.where(assembly_idc == gid)[0] + + +# ------------------------------------------------------------------------------ +experiment = "rewiring_ex3" +sim_date = "191204_132602/17" +branches = [0, 11, 10] +patterns = [0, 3, 6] +master_seed = 5 +xlims = [[0.0, 1.7], [1.5, 3.2], [3.0, 4.7]] + +# t = 1000s +# idc1 = 5424 +# idc2 = 21674 +# idc3 = 43345 +# t = 2000s +idc1 = 10845 +idc2 = 43345 +idc3 = 75845 +# t = 5000s +# idc1 = 27095 +# idc2 = 108345 +# idc3 = 189595 + + +cs.set_figure_size(58 + 6, 45 + 8) + +# ------------------------------------------------------------------------------ +# Directory of simulation results and log files. +input_directory = os.path.join("results", experiment, sim_date, "data") + +# Directory for plots. +plots_directory = os.path.join(input_directory, "..", "plots") +if not os.path.exists(plots_directory): + os.makedirs(plots_directory) + +np.random.seed(master_seed) +num_rows = 3 +node_colors = [] +num_branches = 3 +num_assemblies = 8 +min_plot_weight = 1.3 +num_assemblies_real = 8 +num_branch_nodes = 3 +num_neurons_per_row = 7 +num_neurons_per_assembly = [14] +num_neurons_per_assembly += 7 * [13] +pos_assemblies = [] + +idc_other_assemblies1 = list(range(8)) +idc_other_assemblies2 = list(range(8)) +idc_other_assemblies3 = list(range(8)) +idc_other_assemblies1.remove(patterns[0]) +idc_other_assemblies2.remove(patterns[1]) +idc_other_assemblies3.remove(patterns[2]) + +# Colors of patters. +# c = plt.rcParams['axes.prop_cycle'].by_key()['color'] +# c[6] = c[8] +# c[7] = c[9] +c = sns.color_palette().as_hex() +c[4], c[5], c[6], c[7], c[8], c[9] = c[4], c[8], c[5], c[6], c[8], c[9] + +pattern_labels = [[r"$\mathrm{A}_%d$" % (p + 1) for p in patterns], + [r"$\mathrm{A}_%d$" % (p + 1) for p in patterns]] +pattern_colors = [[c[p] for p in patterns], [c[p] for p in patterns]] +assembly_map = {p: p for p in patterns} +branches.sort() + +x_offset = -0.54 +y_offset = -0.79 +for i in range(num_rows): + for x in np.linspace(-1 + x_offset, 1, 35): + pos_assemblies.append((x, 1.0 - i * 0.1)) +np.random.shuffle(pos_assemblies) + +assemblies = [] +for i in range(num_assemblies): + if i == 0: + assemblies.append(np.arange(14)) + else: + assemblies.append(np.max(assemblies[-1]) + 1 + np.arange(13)) + +branch1_nodes = np.max(assemblies[-1]) + 1 + np.arange(num_branch_nodes) +branch2_nodes = np.max(branch1_nodes) + 1 + np.arange(num_branch_nodes) +branch3_nodes = np.max(branch2_nodes) + 1 + np.arange(num_branch_nodes) +branch_nodes = [branch1_nodes, branch2_nodes, branch3_nodes] + +idc_branch_nodes = range(num_branch_nodes) + +# Load the configuration file. +config = utils.load_configuration(os.path.join( + input_directory, "..", "config_" + experiment + ".yaml")) +sim_simulation_time = config["simulation_time"] +sim_w_max = config["connection_parameters"]["w_max"] +sim_input_size = config["input_parameters"]["num_inputs"] +sim_num_assemblies = config["input_parameters"]["num_assemblies"] +sim_assembly_size = config["input_parameters"]["assembly_size"] +sim_num_branches = config["neuron_parameters"]["num_branches"] +sim_sampling_interval_weights = config["sampling_interval_weights"] +input_size = config["input_parameters"]["num_inputs"] + +assembly_idc = np.split(np.arange(sim_num_assemblies * sim_assembly_size), + sim_num_assemblies) + +# Load the simulation results. +mem_soma = np.loadtxt(os.path.join(input_directory, 'test_soma.0.mem')) + +with open(os.path.join(input_directory, "weights.0.dat"), "rb") as f: + lines = f.readlines() + +weights = np.loadtxt(lines[idc1:idc1 + sim_num_branches]) +weights1 = [weights[b] for b in branches] +weights = np.loadtxt(lines[idc2:idc2 + sim_num_branches]) +weights2 = [weights[b] for b in branches] +weights = np.loadtxt(lines[idc3:idc3 + sim_num_branches]) +weights3 = [weights[b] for b in branches] + +# After assembly 1. +# ------------------------------------------------------------------------------ +np.random.seed(master_seed) +fig = plt.figure() +gs = GridSpec(5, 2) + +# Create graph. +node_colors = [] +G = nx.MultiGraph() +ax = plt.subplot(gs[0:, 0]) + +# Draw the assemblies. +node_colors += draw_assemblies(G, assemblies, c) + +# Draw the neuron. +xc = np.mean(pos_assemblies, axis=0)[0] +pos_branch_nodes, nc = draw_neuron(G, ax, branch_nodes, xc) +node_colors += nc + +# Add connections to graph. +add_connections(experiment, weights1, assemblies, assembly_idc, + assembly_map, idc_other_assemblies1, num_neurons_per_assembly, + idc_branch_nodes, min_plot_weight) + +# Draw connections. +draw_connections(G, (pos_assemblies + pos_branch_nodes), node_colors, ax) + +plt.axis('off') +ax.set_xlim([-1, 1]) +ax.set_ylim([-1 + y_offset, 1]) + +# Soma potential. +ax = plt.subplot(gs[-1, 1:]) +plot_soma_potential(mem_soma, ax, xlims[0], pattern_labels[0], + pattern_colors[0]) + +# Scale +plot_scale(xlims[0]) + +plt.subplots_adjust(left=None, bottom=None, right=None, top=None, + wspace=-0.2, hspace=0.4) + +fname = os.path.join(plots_directory, "graph1") +fig.savefig(fname + ".pdf", pad_inches=0.01) +subprocess.call(["pdftops", "-eps", fname + ".pdf", fname + ".eps"]) +plt.close(fig) + + +# After assembly 2. +# ------------------------------------------------------------------------------ +np.random.seed(master_seed) +fig = plt.figure() +gs = GridSpec(5, 2) + +# Create graph. +node_colors = [] +G = nx.MultiGraph() +ax = plt.subplot(gs[0:, 0]) + +# Draw the assemblies. +node_colors += draw_assemblies(G, assemblies, c) + +# Draw the neuron. +xc = np.mean(pos_assemblies, axis=0)[0] +pos_branch_nodes, nc = draw_neuron(G, ax, branch_nodes, xc) +node_colors += nc + +# Add connections to graph. +add_connections(experiment, weights2, assemblies, assembly_idc, + assembly_map, idc_other_assemblies2, num_neurons_per_assembly, + idc_branch_nodes, min_plot_weight) + +# Draw connections. +draw_connections(G, (pos_assemblies + pos_branch_nodes), node_colors, ax) + +plt.axis('off') +ax.set_xlim([-1, 1]) +ax.set_ylim([-1 + y_offset, 1]) + +# Soma potential. +ax = plt.subplot(gs[-1, 1:]) +plot_soma_potential(mem_soma, ax, xlims[1], pattern_labels[0], + pattern_colors[0]) + +# Scale +plot_scale(xlims[1]) + +plt.subplots_adjust(left=None, bottom=None, right=None, top=None, + wspace=-0.2, hspace=0.4) + +fname = os.path.join(plots_directory, "graph2") +fig.savefig(fname + ".pdf", pad_inches=0.01) +subprocess.call(["pdftops", "-eps", fname + ".pdf", fname + ".eps"]) +plt.close(fig) + +# After assembly 3. +# ------------------------------------------------------------------------------ +np.random.seed(master_seed) +fig = plt.figure() +gs = GridSpec(5, 2) + +# Create graph. +node_colors = [] +G = nx.MultiGraph() +ax = plt.subplot(gs[0:, 0]) + +# Draw the assemblies. +node_colors += draw_assemblies(G, assemblies, c) + +# Draw the neuron. +xc = np.mean(pos_assemblies, axis=0)[0] +pos_branch_nodes, nc = draw_neuron(G, ax, branch_nodes, xc) +node_colors += nc + +# Add connections to graph. +add_connections(experiment, weights3, assemblies, assembly_idc, + assembly_map, idc_other_assemblies3, num_neurons_per_assembly, + idc_branch_nodes, min_plot_weight) + +# Draw connections. +draw_connections(G, (pos_assemblies + pos_branch_nodes), node_colors, ax) + +plt.axis('off') +ax.set_xlim([-1, 1]) +ax.set_ylim([-1 + y_offset, 1]) + +# Soma potential. +ax = plt.subplot(gs[-1, 1:]) +plot_soma_potential(mem_soma, ax, xlims[2], pattern_labels[0], + pattern_colors[0]) + +# Scale +plot_scale(xlims[2]) + +plt.subplots_adjust(left=None, bottom=None, right=None, top=None, + wspace=-0.2, hspace=0.4) + +fname = os.path.join(plots_directory, "graph3") +fig.savefig(fname + ".pdf", pad_inches=0.01) +subprocess.call(["pdftops", "-eps", fname + ".pdf", fname + ".eps"]) +plt.close(fig) diff --git a/utils/plot_sensitivity_analysis.py b/utils/plot_sensitivity_analysis.py new file mode 100644 index 0000000..04d642b --- /dev/null +++ b/utils/plot_sensitivity_analysis.py @@ -0,0 +1,113 @@ +#!/usr/bin/env python +"""TODO""" + +import os +import subprocess + +import plotting.configure_seaborn as cs +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import seaborn as sns + +sns.set(context='paper', style='ticks', rc=cs.rc_params) + +plots_directory = os.path.join(os.path.sep, "calc", os.getenv("USER"), + "data", "dendritic_rewiring", "rewiring_ex2/sensitivity_analysis") + +mean_best = 7.400 +std_best = 0.566 + +# These values have to be set to the values obtained via `stats_rewiring.py`. +mean_assemblies_neg = [7.320, + 7.080, + 7.240, + 7.280, + 7.080, + 7.440, + 7.200, + 6.360, + 7.240, + 7.280, + 7.240, + 7.400 + ] +std_assemblies_neg = [0.676, + 0.891, + 0.650, + 0.665, + 0.744, + 0.753, + 0.748, + 0.975, + 0.709, + 0.776, + 0.763, + 0.566 + ] + +mean_assemblies_pos = [7.400, + 7.280, + 7.200, + 7.040, + 7.080, + 7.040, + 7.160, + 6.400, + 7.600, + 7.280, + 7.360, + 7.400 + ] +std_assemblies_pos = [0.632, + 0.722, + 0.632, + 0.720, + 0.845, + 0.720, + 0.731, + 0.849, + 0.566, + 0.665, + 0.742, + 0.566 + ] + +mean_assemblies = np.array([mean_assemblies_pos, mean_assemblies_neg]) +std_assemblies = np.array([std_assemblies_pos, std_assemblies_neg]) + +diff_mean = np.abs(mean_best - mean_assemblies) +max_diff_mean = np.max(diff_mean, axis=0) +min_diff_mean = np.min(diff_mean, axis=0) +idc = np.argmax(diff_mean, axis=0) + +# sdi = max_diff_mean / np.diag(std_assemblies[idc]) +sdi = 100/10 * max_diff_mean / mean_best +# sdi = 100/10 * min_diff_mean / mean_best +# sdi = 100/10 * np.mean(diff_mean, axis=0) / mean_best + +print(diff_mean) +print(sdi) + +param = [r"$\eta$", r"$T$", r"$c_{\mathcal{L}}$", r"$\gamma$", r"$\lambda$", r"$c_{\mathrm{w}}$", + r"$c_{\mathcal{STDP}}$", r"$\mathcal{STDP}_{\mathrm{th}}$", r"$\tau_{\mathrm{x}}$", + r"$c_{\mathrm{ds}}$", r"$\Delta_{\mathrm{min}}^{\mathrm{ds}}$", + r"$\Delta_{\mathrm{max}}^{\mathrm{ds}}$"] + +df = pd.DataFrame(data={param[i]: sdi[i] for i in range(len(sdi))}, index=range(len(sdi))) + +colors = sns.color_palette().as_hex() +cs.set_figure_size(176 + 7, 60 + 6) + +fig = plt.figure() +ax = sns.barplot(data=df, ci=None, color=colors[0]) +ax.set_ylim([0, 1.4]) +ax.set_ylabel("$\mathrm{SI}\%$") +ax.set_xlabel("Parameter") + +plt.tight_layout() +fname = os.path.join(plots_directory, "sensitivity") +fig.savefig(fname + ".pdf", pad_inches=0.01) +subprocess.call(["pdftops", "-eps", fname + ".pdf", fname + ".eps"]) + +plt.close(fig) diff --git a/utils/stats_rewiring.py b/utils/stats_rewiring.py new file mode 100644 index 0000000..731f4fe --- /dev/null +++ b/utils/stats_rewiring.py @@ -0,0 +1,438 @@ +# /usr/bin/env python + +import glob +import os +import subprocess + +import plotting.configure_seaborn as cs +import matplotlib.colors as mcolors +import matplotlib.lines as mlines +import matplotlib.pyplot as plt +import numpy as np +import seaborn as sns +import utils as utils +from matplotlib.gridspec import GridSpec +from matplotlib.patches import FancyArrowPatch +from mpl_toolkits.axes_grid1 import ImageGrid, make_axes_locatable +from scipy.special import expit + +sns.set(context='paper', style='ticks', rc=cs.rc_params) + + +def _load_configuration(input_directory, experiment, keys): + config = utils.load_configuration(os.path.join( + input_directory, "..", "config_" + experiment + ".yaml")) + + if len(keys) > 1: + return config[keys[0]][keys[1]] + else: + return config[keys[0]] + + +def _load_simulation_results(input_directory, experiment, key, gid=None, + header_lenght=3): + if key == "input_spike_times": + data = np.loadtxt(os.path.join(input_directory, "input.0.ras")) + if key == "input_spike_times_test": + data = np.loadtxt(os.path.join(input_directory, "test_input.0.ras")) + elif key == "output_spike_times": + data = np.loadtxt(os.path.join(input_directory, "output.0.ras")) + elif key == "output_spike_times_test": + data = np.loadtxt(os.path.join(input_directory, "test_output.0.ras")) + elif key == "weights_end": + num_branches = _load_configuration( + input_directory, experiment, ["neuron_parameters", "num_branches"]) + with open(os.path.join(input_directory, "weights.0.dat"), "rb") as f: + lines = f.readlines() + data = np.loadtxt(lines[-num_branches:]) + elif key == "weights_start": + num_branches = _load_configuration( + input_directory, experiment, ["neuron_parameters", "num_branches"]) + with open(os.path.join(input_directory, "weights.0.dat"), "rb") as f: + lines = f.readlines() + data = np.loadtxt(lines[header_lenght:num_branches + header_lenght]) + elif key == "weights": + with open(os.path.join(input_directory, "weights.0.dat"), "rb") as f: + lines = f.readlines() + data = np.loadtxt(lines) + elif key == "branch_mem": + data = np.loadtxt(os.path.join(input_directory, "branch" + str(gid) + + ".0.mem")) + elif key == "branch_mem_test": + data = np.loadtxt(os.path.join(input_directory, "test_branch" + + str(gid) + ".0.mem")) + elif key == "soma_mem": + data = np.loadtxt(os.path.join(input_directory, "soma.0.mem")) + elif key == "soma_mem_test": + data = np.loadtxt(os.path.join(input_directory, "test_soma.0.mem")) + elif key == "plateau_duration": + data = np.loadtxt(os.path.join(input_directory, "branch0.0.pla")) + + return data + + +def _get_represented_assemblies(weights, assembly_idc, min_summed_weight, + min_num_synapses): + # w = utils.reject_outliers(weights) + w = weights + + num_synapses = [] + summed_weight = [] + for i, assembly_idx in enumerate(assembly_idc): + num_synapses.append(sum(np.heaviside(w[assembly_idx], 0))) + summed_weight.append( + sum(w[assembly_idx][w[assembly_idx] > 0], 0)) + + idc = np.where( + (np.asanyarray(summed_weight) >= min_summed_weight) & + (np.asanyarray(num_synapses) >= min_num_synapses))[0].tolist() + + return idc, num_synapses, summed_weight + + +def _get_start_and_stop_times_of_test_patterns(experiment, input_directory): + num_assemblies = _load_configuration(input_directory, experiment, ["input_parameters", "num_assemblies"]) + num_test_patterns_per_assembly = _load_configuration(input_directory, experiment, ["input_parameters", + "num_test_patterns_per_assembly"]) + num_patterns_per_assembly = _load_configuration(input_directory, experiment, ["input_parameters", + "num_patterns_per_assembly"]) + pattern_delay = _load_configuration(input_directory, experiment, ["input_parameters", "pattern_delay"]) + pattern_duration = _load_configuration(input_directory, experiment, ["input_parameters", + "pattern_duration"]) + + # Start and end times of test patterns. + duration_test_patterns = (pattern_delay + pattern_duration) * \ + num_test_patterns_per_assembly * num_assemblies + times_start = np.array( + [((pattern_delay + pattern_duration) * + num_patterns_per_assembly + + i * (pattern_delay + pattern_duration) * + (num_patterns_per_assembly + + (num_assemblies * num_test_patterns_per_assembly))) + for i in range(num_assemblies)]) + times_end = times_start + duration_test_patterns + + return times_start, times_end + + +def _compute_statistics(experiment, input_directory, logfile, min_summed_weight=50, min_num_synapses=10): + num_branches = _load_configuration(input_directory, experiment, ["neuron_parameters", "num_branches"]) + num_assemblies = _load_configuration(input_directory, experiment, ["input_parameters", "num_assemblies"]) + assembly_size = _load_configuration(input_directory, experiment, ["input_parameters", "assembly_size"]) + + weights_end = _load_simulation_results(input_directory, experiment, "weights_end") + + if experiment == "rewiring_ex5" or experiment == "rewiring_ex6": + logfile.write(" Overlap:\n".encode()) + + assembly_neurons_idc = np.loadtxt(os.path.join(input_directory, "assembly_neurons_idc"), + dtype=np.int) + overlap = np.zeros((len(assembly_neurons_idc), len(assembly_neurons_idc))) + for i in range(len(assembly_neurons_idc)): + for j in range(len(assembly_neurons_idc)): + overlap[i, j] = np.intersect1d( + assembly_neurons_idc[i], + assembly_neurons_idc[j]).size / assembly_size + + np.savetxt(logfile, overlap, fmt="%9.3f") + idc = np.triu_indices_from(overlap, k=1) + logfile.write("\n Mean overlap: {0:0.3f} var: {1:0.3f}\n\n".format( + np.mean(overlap[idc]), np.var(overlap[idc])).encode()) + + else: + assembly_neurons_idc = np.split(np.arange(num_assemblies * assembly_size), num_assemblies) + + assemblies = [] + num_synapses_assembly = [] + summed_weight_assembly = [] + num_synapses = [] + summed_weight = [] + num_assemblies_per_branch = np.zeros(num_branches) + for i in range(num_branches): + logfile.write(" Branch: {0}\n".format(i).encode()) + + idc, num_synapses, summed_weight = _get_represented_assemblies( + weights_end[i], assembly_neurons_idc, min_summed_weight, + min_num_synapses) + + if idc: + assemblies.append(idc) + num_assemblies_per_branch[i] = len(idc) + logfile.write(" Assembly: {0}, N_syn: {1}, sum(w): {2}\n" + .format(assemblies[-1], num_synapses, + [round(x) for x in summed_weight]) + .encode()) + num_synapses_assembly.append([num_synapses[a] + for a in assemblies[-1]]) + summed_weight_assembly.append([summed_weight[a] + for a in assemblies[-1]]) + else: + logfile.write(" Assembly: {0}, N_syn: {1}, sum(w): {2}\n" + .format(" ", num_synapses, + [round(x) for x in summed_weight]) + .encode()) + + num_represented_assemblies = len(np.unique([a for l in assemblies + for a in l])) + logfile.write("\n Represented assemblies {0}/{1}, {2}\n".format( + num_represented_assemblies, num_assemblies, assemblies).encode()) + + logfile.write(" Number of synapses per assembly" " {0:0.3f} SD:" + " {1:0.3f}\n".format( + np.mean([n for l in num_synapses_assembly for n in l]), + np.std([n for l in num_synapses_assembly for n in l])) + .encode()) + + logfile.write(" Summed weight per assembly" + " {0:0.3f} SD: {1:0.3f}\n\n".format( + np.mean([s for l in summed_weight_assembly for s in l]), + np.std([s for l in summed_weight_assembly for s in l])) + .encode()) + + if experiment == "rewiring_ex5" or experiment == "rewiring_ex6": + idc = np.triu_indices_from(overlap, k=1) + return (num_represented_assemblies, overlap, num_assemblies_per_branch) + else: + return num_represented_assemblies, num_assemblies_per_branch + + +def _plot_num_represented_assemblies_over_act_probability(a, m, s, plots_directory): + fig = plt.figure() + ax1 = fig.add_subplot(111) + e = ax1.errorbar(range(len(m)), m, s, fmt='o', ms=3, clip_on=False) + + for bar in e[1]: + bar.set_clip_on(False) + for bar in e[2]: + bar.set_clip_on(False) + + ax1.set_xlim([0, len(a) - 1]) + ax1.set_ylim([0, 8]) + ax1.set_xticklabels("%.1f" % f for f in np.linspace(np.max(a), np.min(a), len(a))) + ax1.set_yticklabels("%d" % f for f in range(0, 10, 2)) + ax1.set_ylabel(r"# of represented assemblies") + ax1.set_xlabel("Neuron activation probability") + + ax1.spines["bottom"].set_position(("data", -1.5)) + ax1.spines["left"].set_position(("data", -0.1)) + + plt.tight_layout() + fname = os.path.join(plots_directory, "..", "..", "num-represented-assemblies") + fig.savefig(fname + ".pdf", pad_inches=0.01) + subprocess.call(["pdftops", "-eps", fname + ".pdf", fname + ".eps"]) + plt.close(fig) + + +def _plot_weight_evolution_selected_branches(experiment, input_directory, plots_directory, branches=[0, 1, 2], + patterns=[0, 1, 2], subsample=20, w_thr=7.0, linewidth=0.5, + use_colors=True): + colors = sns.color_palette().as_hex() + if use_colors: + colors_p = ['#9db2d5', '#e7a783', '#a2d0ad', '#c4c5c4'] + else: + colors_p = ['#c4c5c4'] + colors = ['k'] + + num_branches = _load_configuration(input_directory, experiment, + ["neuron_parameters", "num_branches"]) + simulation_time = _load_configuration(input_directory, experiment, + ["simulation_time"]) + w_max = _load_configuration(input_directory, experiment, + ["connection_parameters", "w_max"]) + sampling_interval_weights = _load_configuration( + input_directory, experiment, ["sampling_interval_weights"]) + + weights = _load_simulation_results(input_directory, experiment, "weights") + + if (experiment == "rewiring_ex3" or experiment == "rewiring_ex6"): + t_start, t_end = _get_start_and_stop_times_of_test_patterns(experiment, input_directory) + + exclude_indc = [] + for ts, te in zip(t_start, t_end): + exclude_indc += list(range(int(ts / sampling_interval_weights), + int(te / sampling_interval_weights))) + effective_simulation_time = simulation_time - np.sum(np.array(t_end) - np.array(t_start)) + + else: + exclude_indc = [] + effective_simulation_time = simulation_time + + cs.set_figure_size(176 + 7, 71 + 7) + fig = plt.figure() + gs = GridSpec(len(branches), 1) + + for i, (b, p) in enumerate(zip(branches, patterns)): + ax = plt.subplot(gs[i]) + w = np.delete(weights[range(b, weights.shape[0], num_branches)], + exclude_indc, axis=0) + w[w < 0] = 0 + w = w[::subsample] + + ax.plot(w[:, w[-1, :] < w_thr], color=colors_p[-1], alpha=1.0, + linewidth=linewidth) + + if w[:, w[-1, :] >= w_thr].size: + ax.plot(w[:, w[-1, :] >= w_thr], color=colors_p[p % len(colors_p)], + alpha=1.0, linewidth=linewidth) + ax.plot(np.mean(w[:, w[-1, :] >= w_thr], axis=1), + color=colors[p % len(colors)], linewidth=1.0) + ax.set_xlim([0, effective_simulation_time / (subsample * + sampling_interval_weights)]) + ax.set_xticks(np.linspace(ax.get_xlim()[0], ax.get_xlim()[1], 5)) + ax.set_xticklabels("%.1f" + % (f * subsample * sampling_interval_weights) + for f in np.linspace(ax.get_xlim()[0], + ax.get_xlim()[1], 5)) + ax.set_ylim([0, w_max]) + ax.set_yticks(np.arange(0, w_max + 1, 4)) + ax.set_yticklabels("%d" % f for f in ax.get_yticks()) + ax.set_ylabel(r"$w_{%di}$ [nA]" % (b + 1)) + + if b == branches[-1]: + ax.set_xlabel(r"$t$ [s]") + ax.set_xticklabels("%.1f" + % (f * subsample * sampling_interval_weights) + for f in ax.get_xticks()) + else: + ax.set_xticks([]) + + plt.tight_layout() + name = "weights" + for b in branches: + name += "-b%d" % b + fname = os.path.join(plots_directory, name) + fig.savefig(fname + ".pdf", pad_inches=0.01) + subprocess.call(["pdftops", "-eps", fname + ".pdf", fname + ".eps"]) + plt.close(fig) + + +def main(experiment, sim_date): + + make_plots = False + # make_plots = True + comp_stats = not make_plots + + # Set this if `make_plots` is set to `True`. + branches = [0, 1, 2] # What branches to plot + patterns = [0, 1, 2] # This is just to get the color assigned to the assembly + + # Directories of simulation results and log files. + input_directories = sorted(glob.iglob(os.path.join("results", experiment, sim_date, "data"))) + + paths = {} + for input_directory in input_directories: + date = input_directory.split("/")[-3] + logfile_path = os.path.join("results", experiment, date, "stats.txt") + + if logfile_path in paths: + paths[logfile_path].append(input_directory) + else: + paths[logfile_path] = [input_directory] + + experiment = experiment.split("/")[0] + + neuron_activation_probability = [] + mean_num_represented_assemblies = [] + std_num_represented_assemblies = [] + for logfile_path, input_directories in paths.items(): + if comp_stats: + logfile = open(logfile_path, "wb") + + if experiment == "rewiring_ex4": + neuron_activation_probability.append(_load_configuration( + input_directories[0], experiment, + ["input_parameters", + "neuron_activation_probability"])) + + overlaps = [] + num_assemblies_per_branch = [] + num_represented_assemblies = [] + for trial, input_directory in enumerate(input_directories): + + if os.stat(os.path.join(input_directory, + "weights.0.dat")).st_size == 0: + continue + + plots_directory = os.path.join(input_directory, "..", "plots") + if not os.path.exists(plots_directory): + os.makedirs(plots_directory) + + if comp_stats and experiment != "test_plateau_duration": + # Compute some statistics. + logfile.write("Trial {0} (Data path: {1}):\n".format( + trial, input_directory).encode()) + + result = _compute_statistics(experiment, input_directory, + logfile) + if experiment == "rewiring_ex5" or experiment == "rewiring_ex6": + num_represented_assemblies.append(result[0]) + overlaps.append(result[1]) + num_assemblies_per_branch.append(result[2]) + else: + num_represented_assemblies.append(result[0]) + num_assemblies_per_branch.append(result[1]) + + if make_plots: + # Make some plots. + _plot_weight_evolution_selected_branches(experiment, + input_directory, + plots_directory, + branches=branches, + patterns=patterns, + subsample=1, use_colors=True) + + mean_num_represented_assemblies.append(np.mean(num_represented_assemblies)) + std_num_represented_assemblies.append(np.std(num_represented_assemblies)) + + if comp_stats: + num_assemblies = _load_configuration( + input_directories[0], experiment, ["input_parameters", "num_assemblies"]) + + logfile.write("Represented assemblies {0:0.3f} " + "SD: {1:0.3f}/{2} {3}\n" + .format(mean_num_represented_assemblies[-1], + std_num_represented_assemblies[-1], + num_assemblies, + num_represented_assemblies).encode()) + logfile.write("Mean num assemblies per branch {0:0.3f} " + "SD: {1:0.3f}\n" + .format(np.mean(num_assemblies_per_branch), + np.std(num_assemblies_per_branch)).encode()) + mm = [] + for i in range(num_assemblies + 1): + mm.append([list(num).count(i) for num in + num_assemblies_per_branch]) + m = np.mean(mm, axis=1) + s = np.std(mm, axis=1) + logfile.write("Num assemblies per branch\n" + "0:{0:0.3f} SD: {1:0.3f}\n1:{2:0.3f} SD: {3:0.3f}\n" + "2:{4:0.3f} SD: {5:0.3f}\n3:{6:0.3f} SD: {7:0.3f}\n" + "4:{8:0.3f} SD: {9:0.3f}\n" + "5:{10:0.3f} SD: {11:0.3f}\n" + "6:{12:0.3f} SD: {13:0.3f}\n" + "7:{14:0.3f} SD: {15:0.3f}\n" + "8:{16:0.3f} SD: {17:0.3f}\n" + .format(m[0], s[0], m[1], s[1], m[2], s[2], m[3], + s[3], m[4], s[4], m[5], s[5], m[6], s[6], + m[7], s[7], m[8], s[8]).encode()) + + if (experiment == "rewiring_ex5" or experiment == "rewiring_ex6") and comp_stats: + idc = np.triu_indices_from(overlaps[0], k=1) + triu = [overlaps[i][idc] for i in range(len(overlaps))] + logfile.write("Mean overlap {0:0.4f} SD: {1:0.4f}\n" + .format(np.mean(triu), + np.std(triu)).encode()) + + logfile.close() + + if experiment == "rewiring_ex4" and len(paths) > 1 and comp_stats: + a, m, s = zip(*sorted(zip(neuron_activation_probability, + mean_num_represented_assemblies, + std_num_represented_assemblies), + reverse=True)) + _plot_num_represented_assemblies_over_act_probability(a, m, s, plots_directory) + + +if __name__ == "__main__": + import sys + main(sys.argv[1], sys.argv[2])