diff --git a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_LOAD_DATA-checkpoint.py b/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_LOAD_DATA-checkpoint.py
deleted file mode 100644
index 390dd475d8dd5f890ee4569bcf20da89b3caff8b..0000000000000000000000000000000000000000
--- a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_LOAD_DATA-checkpoint.py
+++ /dev/null
@@ -1,191 +0,0 @@
-import numpy as np
-import subprocess
-import matplotlib.pyplot as plt
-
-from multiarea_model import Analysis
-from config import data_path
-
-def load_and_create_data(M, A, raster_areas):
- """
- Calculate time-averaged population rates and store them in member pop_rates.
- If the rates had previously been stored with the same
- parameters, they are loaded from file.
-
- Parameters
- ----------
- t_min : float, optional
- Minimal time in ms of the simulation to take into account
- for the calculation. Defaults to 500 ms.
- t_max : float, optional
- Maximal time in ms of the simulation to take into account
- for the calculation. Defaults to the simulation time.
- compute_stat : bool, optional
- If set to true, the mean and variance of the population rate
- is calculated. Defaults to False.
- Caution: Setting to True slows down the computation.
- areas : list, optional
- Which areas to include in the calculcation.
- Defaults to all loaded areas.
- pops : list or {'complete'}, optional
- Which populations to include in the calculation.
- If set to 'complete', all populations the respective areas
- are included. Defaults to 'complete'.
- """
- A.create_pop_rates()
- # subprocess.run(['python3', './figures/Schmidt2018_dyn/compute_pop_rates.py'])
- # subprocess.run(['Rscript', '--vanilla', 'compute_bold_signal.R', fn, out_fn])
- # print("Computing population rates done")
-
-
- """
- Calculate poulation-averaged LvR (see Shinomoto et al. 2009) and
- store as member pop_LvR. Uses helper function LvR.
-
- Parameters
- ----------
- t_min : float, optional
- Minimal time in ms of the simulation to take into account
- for the calculation. Defaults to 500 ms.
- t_max : float, optional
- Maximal time in ms of the simulation to take into account
- for the calculation. Defaults to the simulation time.
- areas : list, optional
- Which areas to include in the calculcation.
- Defaults to all loaded areas.
- pops : list or {'complete'}, optional
- Which populations to include in the calculation.
- If set to 'complete', all populations the respective areas
- are included. Defaults to 'complete'.
- """
- A.create_pop_LvR()
- # print("Computing population LvR done")
-
-
- """
- Calculate time series of population- and area-averaged firing rates.
- Uses ah.pop_rate_time_series.
- If the rates have previously been stored with the
- same parameters, they are loaded from file.
-
-
- Parameters
- ----------
- t_min : float, optional
- Minimal time in ms of the simulation to take into account
- for the calculation. Defaults to 500 ms.
- t_max : float, optional
- Maximal time in ms of the simulation to take into account
- for the calculation. Defaults to the simulation time.
- areas : list, optional
- Which areas to include in the calculcation.
- Defaults to all loaded areas.
- pops : list or {'complete'}, optional
- Which populations to include in the calculation.
- If set to 'complete', all populations the respective areas
- are included. Defaults to 'complete'.
- kernel : {'gauss_time_window', 'alpha_time_window', 'rect_time_window'}, optional
- Specifies the kernel to be convolved with the spike histogram.
- Defaults to 'binned', which corresponds to no convolution.
- resolution: float, optional
- Width of the convolution kernel. Specifically it correponds to:
- - 'binned' : bin width of the histogram
- - 'gauss_time_window' : sigma
- - 'alpha_time_window' : time constant of the alpha function
- - 'rect_time_window' : width of the moving rectangular function
- """
- A.create_rate_time_series()
- # print("Computing rate time series done")
-
-
- # # Compute time series of firing rates convolved with a kernel
- # # data_path = sys.argv[1]
- # print(data_path)
- # label = M.simulation.label
- # method = 'auto_kernel'
- # for area in raster_areas:
- # subprocess.run(['python3', './figures/Schmidt2018_dyn/compute_rate_time_series.py', data_path, label, area, method])
- # print("Computing rate time series auto kernel done")
-
-
- """
- Calculate synchrony as the coefficient of variation of the population rate
- and store in member synchrony. Uses helper function synchrony.
- If the synchrony has previously been stored with the
- same parameters, they are loaded from file.
-
-
- Parameters
- ----------
- t_min : float, optional
- Minimal time in ms of the simulation to take into account
- for the calculation. Defaults to 500 ms.
- t_max : float, optional
- Maximal time in ms of the simulation to take into account
- for the calculation. Defaults to the simulation time.
- areas : list, optional
- Which areas to include in the calculcation.
- Defaults to all loaded areas.
- pops : list or {'complete'}, optional
- Which populations to include in the calculation.
- If set to 'complete', all populations the respective areas
- are included. Defaults to 'complete'.
- resolution : float, optional
- Resolution of the population rate. Defaults to 1 ms.
- """
- A.create_synchrony()
- # print("Computing synchrony done")
-
-
-# # Create corrrelation coefficient
-# # data_path = data_path
-# label = M.simulation.label
-
-# # Run the script with arguments
-# subprocess.run(['python', './figures/Schmidt2018_dyn/compute_corrcoeff.py', data_path, label])
-
-
-# """
-# Calculate synaptic input of populations and areas using the spike data.
-# Uses function ah.pop_synaptic_input.
-# If the synaptic inputs have previously been stored with the
-# same parameters, they are loaded from file.
-
-# Parameters
-# ----------
-# t_min : float, optional
-# Minimal time in ms of the simulation to take into account
-# for the calculation. Defaults to 500 ms.
-# t_max : float, optional
-# Maximal time in ms of the simulation to take into account
-# for the calculation. Defaults to the simulation time.
-# areas : list, optional
-# Which areas to include in the calculcation.
-# Defaults to all loaded areas.
-# pops : list or {'complete'}, optional
-# Which populations to include in the calculation.
-# If set to 'complete', all populations the respective areas
-# are included. Defaults to 'complete'.
-# kernel : {'gauss_time_window', 'alpha_time_window', 'rect_time_window'}, optional
-# Convolution kernel for the calculation of the underlying firing rates.
-# Defaults to 'binned' which corresponds to a simple histogram.
-# resolution: float, optional
-# Width of the convolution kernel. Specifically it correponds to:
-# - 'binned' : bin width of the histogram
-# - 'gauss_time_window' : sigma
-# - 'alpha_time_window' : time constant of the alpha function
-# - 'rect_time_window' : width of the moving rectangular function
-# """
- # A.create_synaptic_input()
- # # print("Computing synaptic input done")
-
- A.save()
-
- # """
- # Compute BOLD signal for a given area from the time series of
- # population-averaged spike rates of a given simulation using the
- # neuRosim package of R (see Schmidt et al. 2018 for more details).
- # """
- # try:
- # subprocess.run(['python3', './Schmidt2018_dyn/compute_bold_signal.py'])
- # except FileNotFoundError:
- # raise FileNotFoundError("Executing R failed. Did you install R?")
\ No newline at end of file
diff --git a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_instantaneous_and_mean_firing_rates-checkpoint.py b/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_instantaneous_and_mean_firing_rates-checkpoint.py
deleted file mode 100644
index c4e75fc5ff040e9b0e9e3a69b9d4f2587c3d8ae0..0000000000000000000000000000000000000000
--- a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_instantaneous_and_mean_firing_rates-checkpoint.py
+++ /dev/null
@@ -1,21 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as pl
-
-def plot_instan_mean_firing_rate(M):
- # load spike data and calculate instantaneous and mean firing rates
- data = np.loadtxt(M.simulation.data_dir + '/recordings/' + M.simulation.label + "-spikes-1-0.dat", skiprows=3)
- tsteps, spikecount = np.unique(data[:,1], return_counts=True)
- rate = spikecount / M.simulation.params['dt'] * 1e3 / np.sum(M.N_vec)
-
- # visualize calculate instantaneous and mean firing rates
- fig = pl.figure(figsize=(10, 5))
- ax = fig.add_subplot()
- ax.plot(tsteps, rate)
- ax.plot(tsteps, np.average(rate)*np.ones(len(tsteps)), label='mean')
- ax.set_title('Instantaneous and mean firing rate across all populations', fontsize=15)
- ax.set_xlabel('Time (ms)', fontsize=13)
- ax.set_ylabel('Firing rate (spikes / s)', fontsize=12)
- ax.set_xlim(0, M.simulation.params['t_sim'])
- ax.set_ylim(0, 50)
- ax.legend()
- pl.show()
\ No newline at end of file
diff --git a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_interareal_connectivity-checkpoint.py b/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_interareal_connectivity-checkpoint.py
deleted file mode 100644
index ebf1da13faeed741f3c614a6a267b802e2bd2b2c..0000000000000000000000000000000000000000
--- a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_interareal_connectivity-checkpoint.py
+++ /dev/null
@@ -1,353 +0,0 @@
-import json
-import numpy as np
-import matplotlib.pyplot as pl
-import os
-
-import sys
-sys.path.append('./figures/Schmidt2018')
-
-from helpers import area_list, datapath
-from matplotlib import gridspec
-from matplotlib.colors import LogNorm
-from matplotlib.ticker import FixedLocator
-from matplotlib import rc_file
-from multiarea_model import MultiAreaModel
-from plotcolors import myblue
-from scipy import stats
-
-# rc_file('plotstyle.rc')
-
-def visualize_interareal_connectivity(M):
- # full-scale model
- M_full_scale = MultiAreaModel({})
-
- """
- Figure layout
- """
- # nrows = 2
- nrows = 1
- ncols = 2
- # width = 6.8556
- width = 15
- panel_wh_ratio = 0.7 * (1. + np.sqrt(5)) / 2. # golden ratio
-
- height = width / panel_wh_ratio * float(nrows) / ncols
- # print(width, height)
- pl.rcParams['figure.figsize'] = (width, height)
-
- fig = pl.figure()
- fig.suptitle('Interareal connectivity for full-scale (left) and down-scale (right) multi-area model', fontsize=16, y=1.05)
- axes = {}
-
- # gs1 = gridspec.GridSpec(2, 2)
- gs1 = gridspec.GridSpec(1, 2)
- # gs1.update(left=0.06, right=0.95, top=0.95, bottom=0.1, wspace=0.1, hspace=0.3)
- gs1.update(left=0.06, right=0.95, top=0.95, bottom=0.1, wspace=0.3, hspace=0.3)
-
- # axes['A'] = pl.subplot(gs1[:1, :1])
- # axes['B'] = pl.subplot(gs1[:1, 1:2])
- axes['B'] = pl.subplot(gs1[:1, :1])
- # axes['D'] = pl.subplot(gs1[:1, 2:])
- axes['D'] = pl.subplot(gs1[:1, 1:2])
-
- # pos = axes['A'].get_position()
- pos2 = axes['D'].get_position()
- # axes['C'] = pl.axes([pos.x0 + 0.01, pos2.y0, pos.x1 - pos.x0 - 0.025, 0.23])
-
- # print(pos.x1 - pos.x0 - 0.025)
-
- # labels = ['A', 'B', 'C', 'D']
- labels = ['B', 'D']
- labels_display = ['Full-scale model', 'Down-scale model']
- # for label in labels:
- for i in range(len(labels)):
- label = labels[i]
- label_display = labels_display[i]
- if label in ['C']:
- label_pos = [-0.045, 1.18]
- else:
- label_pos = [-0.2, 1.04]
- # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label,
- # fontdict={'fontsize': 10, 'weight': 'bold',
- # 'horizontalalignment': 'left', 'verticalalignment':
- # 'bottom'}, transform=axes[label].transAxes)
- pl.text(label_pos[0], label_pos[1], label_display,
- fontdict={'fontsize': 12, 'weight': 'bold',
- 'horizontalalignment': 'left', 'verticalalignment':
- 'bottom'}, transform=axes[label].transAxes)
-
-# """
-# Load data
-# """
-# # M = MultiAreaModel({})
-# M_full_scale = MultiAreaModel({})
-
-# datapath = './multiarea_model/data_multiarea/'
-# with open(os.path.join(datapath, 'viscortex_processed_data.json'), 'r') as f:
-# proc = json.load(f)
-# with open(os.path.join(datapath, 'viscortex_raw_data.json'), 'r') as f:
-# raw = json.load(f)
-
-# FLN_Data_FV91 = proc['FLN_Data_FV91']
-
-# cocomac_data = raw['cocomac_data']
-# median_distance_data = raw['median_distance_data']
-
-# cocomac = np.zeros((32, 32))
-# conn_matrix = np.zeros((32, 32))
-# for i, area1 in enumerate(area_list[::-1]):
-# for j, area2 in enumerate(area_list):
-# # if M.K_areas[area1][area2] > 0. and area2 in cocomac_data[area1]:
-# if M_full_scale.K_areas[area1][area2] > 0. and area2 in cocomac_data[area1]:
-# cocomac[i][j] = 1.
-# if area2 in FLN_Data_FV91[area1]:
-# conn_matrix[i][j] = FLN_Data_FV91[area1][area2]
-
-# """
-# Panel A: CoCoMac Data
-# """
-# ax = axes['A']
-# # ax.yaxis.set_ticks_position("left")
-# # ax.xaxis.set_ticks_position("bottom")
-
-# ax.set_aspect(1. / ax.get_data_ratio())
-# ax.yaxis.set_ticks_position("none")
-# ax.xaxis.set_ticks_position("none")
-
-# masked_matrix = np.ma.masked_values(cocomac, 0.0)
-# cmap = pl.cm.binary
-# cmap.set_bad('w', 1.0)
-
-# x = np.arange(0, len(area_list) + 1)
-# y = np.arange(0, len(area_list[::-1]) + 1)
-# X, Y = np.meshgrid(x, y)
-
-# ax.set_xticks([i + 0.5 for i in np.arange(0, len(area_list) + 1, 1)])
-# ax.set_xticks([i + 0.5 for i in np.arange(0, len(area_list), 1)])
-# # ax.set_xticklabels(area_list, rotation=90, size=6.)
-# ax.set_xticklabels(area_list, rotation=90, size=10.)
-
-# # ax.set_yticks([i + 0.5 for i in np.arange(0, len(area_list) + 1, 1)])
-# ax.set_yticks([i + 0.5 for i in np.arange(0, len(area_list), 1)])
-# # ax.set_yticklabels(area_list[::-1], size=6.)
-# ax.set_yticklabels(area_list[::-1], size=6.)
-
-# ax.set_ylabel('Target area')
-# ax.set_xlabel('Source area')
-
-# im = ax.pcolormesh(masked_matrix, cmap=cmap,
-# edgecolors='None', vmin=0., vmax=1.)
-
-# t = FixedLocator([])
-# cbar = pl.colorbar(im, ticks=t, fraction=0.046, ax=ax)
-# cbar.set_alpha(0.)
- # cbar.remove()
-
-# """
-# Panel B: Data from Markov et al. (2014) "A weighted and directed
-# interareal connectivity matrix for macaque cerebral cortex."
-# Cerebral Cortex, 24(1), 17–36.
-# """
-# ax = axes['B']
-# ax.set_aspect(1. / ax.get_data_ratio())
-# ax.yaxis.set_ticks_position("none")
-# ax.xaxis.set_ticks_position("none")
-
-# masked_matrix = np.ma.masked_values(conn_matrix, 0.0)
-# cmap = pl.get_cmap('inferno')
-# cmap.set_bad('w', 1.0)
-
-# x = np.arange(0, len(area_list) + 1)
-# y = np.arange(0, len(area_list[::-1]) + 1)
-# X, Y = np.meshgrid(x, y)
-
-# ax.set_xticks([i + 0.5 for i in np.arange(0, len(area_list) + 1, 1)])
-# ax.set_xticklabels(area_list, rotation=90, size=6.)
-
-# ax.set_yticks([i + 0.5 for i in np.arange(0, len(area_list) + 1, 1)])
-# ax.set_yticklabels(area_list[::-1], size=6.)
-
-# im = ax.pcolormesh(masked_matrix, cmap=cmap,
-# edgecolors='None', norm=LogNorm(vmin=1e-6, vmax=1.))
-
-# t = FixedLocator([1e-6, 1e-4, 1e-2, 1])
-# cbar = pl.colorbar(im, ticks=t, fraction=0.046, ax=ax)
-# cbar.set_alpha(0.)
-
- """
- Panel B: Interareal connectivity of full-scaling multi-area model
- """
- conn_matrix_full_scale = np.zeros((32, 32))
- for i, area1 in enumerate(area_list[::-1]):
- for j, area2 in enumerate(area_list):
- conn_matrix_full_scale[i][j] = M_full_scale.K_areas[area1][
- area2] / np.sum(list(M_full_scale.K_areas[area1].values()))
-
- ax = axes['B']
- ax.yaxis.set_ticks_position("none")
- ax.xaxis.set_ticks_position("none")
-
- ax.set_aspect(1. / ax.get_data_ratio())
-
- masked_matrix_full_scale = np.ma.masked_values(conn_matrix_full_scale, 0.0)
- cmap = pl.get_cmap('inferno')
- cmap.set_bad('w', 1.0)
-
- x = np.arange(0, len(area_list) + 1)
- y = np.arange(0, len(area_list[::-1]) + 1)
- X, Y = np.meshgrid(x, y)
-
- # ax.set_xticks([i + 0.5 for i in np.arange(0, len(area_list) + 1, 1)])
- ax.set_xticks([i + 0.5 for i in np.arange(0, len(area_list), 1)])
- # ax.set_xticklabels(area_list, rotation=90, size=6.)
- ax.set_xticklabels(area_list, rotation=90, size=10.)
-
- # ax.set_yticks([i + 0.5 for i in np.arange(0, len(area_list) + 1, 1)])
- ax.set_yticks([i + 0.5 for i in np.arange(0, len(area_list), 1)])
- # ax.set_yticklabels(area_list[::-1], size=6.)
- ax.set_yticklabels(area_list[::-1], size=10.)
-
- # ax.set_ylabel('Target area')
- # ax.set_xlabel('Source area')
- ax.set_ylabel('Target area', fontsize=15)
- ax.set_xlabel('Source area', fontsize=15)
- im = ax.pcolormesh(masked_matrix_full_scale, cmap=cmap,
- edgecolors='None', norm=LogNorm(vmin=1e-6, vmax=1.))
-
- t = FixedLocator([1e-6, 1e-4, 1e-2, 1])
- cbar = pl.colorbar(im, ticks=t, fraction=0.046, ax=ax)
- cbar.set_alpha(0.)
-
- # """
- # Panel C: Exponential decay of FLN with distance
- # """
- # FLN_values_FV91 = np.array([])
- # distances_FV91 = np.array([])
-
- # for target_area in FLN_Data_FV91:
- # for source_area in FLN_Data_FV91[target_area]:
- # if target_area in median_distance_data and source_area in median_distance_data:
- # if FLN_Data_FV91[target_area][source_area]:
- # FLN_values_FV91 = np.append(FLN_values_FV91, FLN_Data_FV91[
- # target_area][source_area])
- # distances_FV91 = np.append(distances_FV91, median_distance_data[
- # target_area][source_area])
-
- # # Linear fit of distances vs. log FLN
- # print("\n \n Linear fit to logarithmic values")
- # gradient, intercept, r_value, p_value, std_err = stats.linregress(
- # distances_FV91, np.log(FLN_values_FV91))
- # print("Raw parameters: ", gradient, intercept)
- # print("Transformed parameters: ", -gradient, np.exp(intercept))
- # print('r_value**2', r_value ** 2)
- # print('p_value', p_value)
- # print('std_err', std_err)
-
- # ax = axes['C']
- # ax.yaxis.set_ticks_position("left")
- # ax.xaxis.set_ticks_position("bottom")
-
- # ax.yaxis.set_ticks_position("left")
- # ax.xaxis.set_ticks_position("bottom")
-
- # ax.spines['right'].set_color('none')
- # ax.spines['top'].set_color('none')
- # ax.yaxis.set_ticks_position("left")
- # ax.xaxis.set_ticks_position("bottom")
-
- # ax.plot(distances_FV91, np.log10(FLN_values_FV91), '.', color=myblue)
- # x = np.arange(np.min(distances_FV91), np.max(distances_FV91), 1)
- # ax.plot(x, (intercept + gradient * x) / np.log(10), linewidth=2.0,
- # color='Black', label='Linear regression fit')
-
- # ax.set_xlabel('Distance (mm)', labelpad=7)
- # ax.set_ylabel(r'$\log(FLN)$')
- # ax.set_yticks([-6, -4, -2, 0])
-
- # print("log fit")
- # print(np.corrcoef(gradient * distances_FV91 + intercept, np.log(FLN_values_FV91))[0][1])
-
- # """
- # Panel D: Resulting connectivity matrix
- # """
- # conn_matrix = np.zeros((32, 32))
- # for i, area1 in enumerate(area_list[::-1]):
- # for j, area2 in enumerate(area_list):
- # conn_matrix[i][j] = M.K_areas[area1][
- # area2] / np.sum(list(M.K_areas[area1].values()))
-
- # ax = axes['D']
- # ax.yaxis.set_ticks_position("none")
- # ax.xaxis.set_ticks_position("none")
-
- # ax.set_aspect(1. / ax.get_data_ratio())
-
- # masked_matrix = np.ma.masked_values(conn_matrix, 0.0)
- # cmap = pl.get_cmap('inferno')
- # cmap.set_bad('w', 1.0)
-
- # x = np.arange(0, len(area_list) + 1)
- # y = np.arange(0, len(area_list[::-1]) + 1)
- # X, Y = np.meshgrid(x, y)
-
- # ax.set_xticks([i + 0.5 for i in np.arange(0, len(area_list) + 1, 1)])
- # ax.set_xticklabels(area_list, rotation=90, size=6.)
-
- # ax.set_yticks([i + 0.5 for i in np.arange(0, len(area_list) + 1, 1)])
- # ax.set_yticklabels(area_list[::-1], size=6.)
-
- # ax.set_ylabel('Target area')
- # ax.set_xlabel('Source area')
- # im = ax.pcolormesh(masked_matrix, cmap=cmap,
- # edgecolors='None', norm=LogNorm(vmin=1e-6, vmax=1.))
-
- # t = FixedLocator([1e-6, 1e-4, 1e-2, 1])
- # cbar = pl.colorbar(im, ticks=t, fraction=0.046, ax=ax)
- # cbar.set_alpha(0.)
-
- """
- Panel D: Interareal connectivity of down-scaling multi-area model
- """
- conn_matrix_down_scale = np.zeros((32, 32))
- for i, area1 in enumerate(area_list[::-1]):
- for j, area2 in enumerate(area_list):
- conn_matrix_down_scale[i][j] = M.K_areas[area1][
- area2] / np.sum(list(M.K_areas[area1].values()))
-
- ax = axes['D']
- ax.yaxis.set_ticks_position("none")
- ax.xaxis.set_ticks_position("none")
-
- ax.set_aspect(1. / ax.get_data_ratio())
-
- masked_matrix_down_scale = np.ma.masked_values(conn_matrix_down_scale, 0.0)
- cmap = pl.get_cmap('inferno')
- cmap.set_bad('w', 1.0)
-
- x = np.arange(0, len(area_list) + 1)
- y = np.arange(0, len(area_list[::-1]) + 1)
- X, Y = np.meshgrid(x, y)
-
- # ax.set_xticks([i + 0.5 for i in np.arange(0, len(area_list) + 1, 1)])
- ax.set_xticks([i + 0.5 for i in np.arange(0, len(area_list), 1)])
- # ax.set_xticklabels(area_list, rotation=90, size=6.)
- ax.set_xticklabels(area_list, rotation=90, size=10.)
-
- # ax.set_yticks([i + 0.5 for i in np.arange(0, len(area_list) + 1, 1)])
- ax.set_yticks([i + 0.5 for i in np.arange(0, len(area_list), 1)])
- # ax.set_yticklabels(area_list[::-1], size=6.)
- ax.set_yticklabels(area_list[::-1], size=10.)
-
- ax.set_ylabel('Target area', fontsize=15)
- ax.set_xlabel('Source area', fontsize=15)
- im = ax.pcolormesh(masked_matrix_down_scale, cmap=cmap,
- edgecolors='None', norm=LogNorm(vmin=1e-6, vmax=1.))
-
- t = FixedLocator([1e-6, 1e-4, 1e-2, 1])
- cbar = pl.colorbar(im, ticks=t, fraction=0.046, ax=ax)
- cbar.set_alpha(0.)
-
- # """
- # Save figure
- # """
- # pl.savefig('Fig4_connectivity.eps')
\ No newline at end of file
diff --git a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_resting_state-checkpoint.py b/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_resting_state-checkpoint.py
deleted file mode 100644
index b19ae9c826974b05369fa769eb9598d24b3dcd86..0000000000000000000000000000000000000000
--- a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_resting_state-checkpoint.py
+++ /dev/null
@@ -1,474 +0,0 @@
-import json
-import numpy as np
-import os
-
-import sys
-sys.path.append('./figures/Schmidt2018_dyn')
-
-# from helpers import original_data_path, population_labels
-from helpers import population_labels
-from multiarea_model import MultiAreaModel
-from multiarea_model import Analysis
-from plotcolors import myred, myblue
-
-import matplotlib.pyplot as pl
-from matplotlib import gridspec
-# from matplotlib import rc_file
-# rc_file('plotstyle.rc')
-
-icolor = myred
-ecolor = myblue
-
-from M2E_LOAD_DATA import load_and_create_data
-
-def set_boxplot_props(d):
- for i in range(len(d['boxes'])):
- if i % 2 == 0:
- d['boxes'][i].set_facecolor(icolor)
- d['boxes'][i].set_color(icolor)
- else:
- d['boxes'][i].set_facecolor(ecolor)
- d['boxes'][i].set_color(ecolor)
- pl.setp(d['whiskers'], color='k')
- pl.setp(d['fliers'], color='k', markerfacecolor='k', marker='+')
- pl.setp(d['medians'], color='none')
- pl.setp(d['caps'], color='k')
- pl.setp(d['means'], marker='x', color='k',
- markerfacecolor='k', markeredgecolor='k', markersize=3.)
-
-def plot_resting_state(M, data_path, raster_areas=['V1', 'V2', 'FEF']):
- """
- Analysis class.
- An instance of the analysis class for the given network and simulation.
- Can be created as a member class of a multiarea_model instance or standalone.
-
- Parameters
- ----------
- network : MultiAreaModel
- An instance of the multiarea_model class that specifies
- the network to be analyzed.
- simulation : Simulation
- An instance of the simulation class that specifies
- the simulation to be analyzed.
- data_list : list of strings {'spikes', vm'}, optional
- Specifies which type of data is to load. Defaults to ['spikes'].
- load_areas : list of strings with area names, optional
- Specifies the areas for which data is to be loaded.
- Default value is None and leads to loading of data for all
- simulated areas.
- """
- # Instantiate an analysis class and load spike data
- A = Analysis(network=M,
- simulation=M.simulation,
- data_list=['spikes'],
- load_areas=None)
-
- # load data
- load_and_create_data(M, A, raster_areas)
-
- t_sim = M.simulation.params["t_sim"]
-
- """
- Figure layout
- """
-
- nrows = 4
- ncols = 4
- # width = 7.0866
- width = 12
- panel_wh_ratio = 0.7 * (1. + np.sqrt(5)) / 2. # golden ratio
-
- height = width / panel_wh_ratio * float(nrows) / ncols
- pl.rcParams['figure.figsize'] = (width, height)
-
-
- fig = pl.figure()
- axes = {}
-
- gs1 = gridspec.GridSpec(1, 3)
- gs1.update(left=0.06, right=0.72, top=0.95, wspace=0.4, bottom=0.35)
- # axes['A'] = pl.subplot(gs1[:-1, :1])
- # axes['B'] = pl.subplot(gs1[:-1, 1:2])
- # axes['C'] = pl.subplot(gs1[:-1, 2:])
- axes['A'] = pl.subplot(gs1[:1, :1])
- axes['B'] = pl.subplot(gs1[:1, 1:2])
- axes['C'] = pl.subplot(gs1[:1, 2:])
-
- gs2 = gridspec.GridSpec(3, 1)
- gs2.update(left=0.78, right=0.95, top=0.95, bottom=0.35)
- axes['D'] = pl.subplot(gs2[:1, :1])
- axes['E'] = pl.subplot(gs2[1:2, :1])
- axes['F'] = pl.subplot(gs2[2:3, :1])
-
-
- gs3 = gridspec.GridSpec(1, 1)
- gs3.update(left=0.1, right=0.95, top=0.3, bottom=0.075)
- axes['G'] = pl.subplot(gs3[:1, :1])
-
- # areas = ['V1', 'V2', 'FEF']
- area_list = ['V1', 'V2', 'VP', 'V3', 'V3A', 'MT', 'V4t', 'V4', 'VOT', 'MSTd',
- 'PIP', 'PO', 'DP', 'MIP', 'MDP', 'VIP', 'LIP', 'PITv', 'PITd',
- 'MSTl', 'CITv', 'CITd', 'FEF', 'TF', 'AITv', 'FST', '7a', 'STPp',
- 'STPa', '46', 'AITd', 'TH']
- if len(raster_areas) !=3:
- raise Exception("Error! Please give 3 areas to display as raster plots.")
- for area in raster_areas:
- if area not in area_list:
- raise Exception("Error! Given raster areas are either not from complete_area_list, please input correct areas to diaply the raster plots.")
- areas = raster_areas
-
- labels = ['A', 'B', 'C']
- for area, label in zip(areas, labels):
- label_pos = [-0.2, 1.01]
- # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label + ': ' + area,
- # fontdict={'fontsize': 10, 'weight': 'bold',
- # 'horizontalalignment': 'left', 'verticalalignment':
- # 'bottom'}, transform=axes[label].transAxes)
- pl.text(label_pos[0], label_pos[1], label + ': ' + area,
- fontdict={'fontsize': 10, 'weight': 'bold',
- 'horizontalalignment': 'left', 'verticalalignment':
- 'bottom'}, transform=axes[label].transAxes)
-
- label = 'G'
- label_pos = [-0.1, 0.92]
- # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label,
- # fontdict={'fontsize': 10, 'weight': 'bold',
- # 'horizontalalignment': 'left', 'verticalalignment':
- # 'bottom'}, transform=axes[label].transAxes)
- pl.text(label_pos[0], label_pos[1], label,
- fontdict={'fontsize': 10, 'weight': 'bold',
- 'horizontalalignment': 'left', 'verticalalignment':
- 'bottom'}, transform=axes[label].transAxes)
-
- labels = ['E', 'D', 'F']
- for label in labels:
- label_pos = [-0.2, 1.05]
- # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label,
- # fontdict={'fontsize': 10, 'weight': 'bold',
- # 'horizontalalignment': 'left', 'verticalalignment':
- # 'bottom'}, transform=axes[label].transAxes)
- pl.text(label_pos[0], label_pos[1], label,
- fontdict={'fontsize': 10, 'weight': 'bold',
- 'horizontalalignment': 'left', 'verticalalignment':
- 'bottom'}, transform=axes[label].transAxes)
-
-
- labels = ['A', 'B', 'C', 'D', 'E', 'F']
-
- for label in labels:
- axes[label].spines['right'].set_color('none')
- axes[label].spines['top'].set_color('none')
- axes[label].yaxis.set_ticks_position("left")
- axes[label].xaxis.set_ticks_position("bottom")
-
- for label in ['A', 'B', 'C']:
- axes[label].yaxis.set_ticks_position('none')
-
-
-# """
-# Load data
-# """
-# LOAD_ORIGINAL_DATA = True
-
-
-# if LOAD_ORIGINAL_DATA:
-# # use T=10500 simulation for spike raster plots
-# label_spikes = '3afaec94d650c637ef8419611c3f80b3cb3ff539'
-# # and T=100500 simulation for all other panels
-# label = '99c0024eacc275d13f719afd59357f7d12f02b77'
-# data_path = original_data_path
-# else:
-# from network_simulations import init_models
-# from config import data_path
-# models = init_models('Fig5')
-# label_spikes = models[0].simulation.label
-# label = models[1].simulation.label
-
- # """
- # Create MultiAreaModel instance to have access to data structures
- # """
- # M = MultiAreaModel({})
-
- # spike data
- # spike_data = {}
- # for area in areas:
- # spike_data[area] = {}
- # for pop in M.structure[area]:
- # spike_data[area][pop] = np.load(os.path.join(data_path,
- # label_spikes,
- # 'recordings',
- # '{}-spikes-{}-{}.npy'.format(label_spikes,
- # area, pop)))
- spike_data = A.spike_data
- label_spikes = M.simulation.label
- label = M.simulation.label
-
- # stationary firing rates
- fn = os.path.join(data_path, label, 'Analysis', 'pop_rates.json')
- with open(fn, 'r') as f:
- pop_rates = json.load(f)
-
- # time series of firing rates
- rate_time_series = {}
- for area in areas:
- # fn = os.path.join(data_path, label,
- # 'Analysis',
- # 'rate_time_series_full',
- # 'rate_time_series_full_{}.npy'.format(area))
- fn = os.path.join(data_path, label,
- 'Analysis',
- 'rate_time_series-{}.npy'.format(area))
- rate_time_series[area] = np.load(fn)
-
- # # time series of firing rates convolved with a kernel
- # rate_time_series_auto_kernel = {}
- # for area in areas:
- # fn = os.path.join(data_path, label,
- # 'Analysis',
- # 'rate_time_series_auto_kernel',
- # 'rate_time_series_auto_kernel_{}.npy'.format(area))
- # rate_time_series_auto_kernel[area] = np.load(fn)
-
- # local variance revised (LvR)
- fn = os.path.join(data_path, label, 'Analysis', 'pop_LvR.json')
- with open(fn, 'r') as f:
- pop_LvR = json.load(f)
-
- # correlation coefficients
- # fn = os.path.join(data_path, label, 'Analysis', 'corrcoeff.json')
- fn = os.path.join(data_path, label, 'Analysis', 'synchrony.json')
- with open(fn, 'r') as f:
- corrcoeff = json.load(f)
-
- """
- Plotting
- """
- # print("Raster plots")
-
- # t_min = 3000.
- # t_max = 3500.
- t_min = t_sim - 500
- t_max = t_sim
-
- icolor = myred
- ecolor = myblue
-
- # frac_neurons = 0.03
- frac_neurons = 1
-
- for i, area in enumerate(areas):
- ax = axes[labels[i]]
-
- if area in spike_data:
- n_pops = len(spike_data[area])
- # Determine number of neurons that will be plotted for this area (for
- # vertical offset)
- offset = 0
- n_to_plot = {}
- for pop in M.structure[area]:
- n_to_plot[pop] = int(M.N[area][pop] * frac_neurons)
- offset = offset + n_to_plot[pop]
- y_max = offset + 1
- prev_pop = ''
- yticks = []
- yticklocs = []
- for jj, pop in enumerate(M.structure[area]):
- if pop[0:-1] != prev_pop:
- prev_pop = pop[0:-1]
- yticks.append('L' + population_labels[jj][0:-1])
- yticklocs.append(offset - 0.5 * n_to_plot[pop])
- ind = np.where(np.logical_and(
- spike_data[area][pop][:, 1] <= t_max, spike_data[area][pop][:, 1] >= t_min))
- pop_data = spike_data[area][pop][ind]
- pop_neurons = np.unique(pop_data[:, 0])
- neurons_to_ = np.arange(np.min(spike_data[area][pop][:, 0]), np.min(
- spike_data[area][pop][:, 0]) + n_to_plot[pop], 1)
-
- if pop.find('E') > (-1):
- pcolor = ecolor
- else:
- pcolor = icolor
-
- for kk in range(n_to_plot[pop]):
- spike_times = pop_data[pop_data[:, 0] == neurons_to_[kk], 1]
-
- _ = ax.plot(spike_times, np.zeros(len(spike_times)) +
- offset - kk, '.', color=pcolor, markersize=1)
- offset = offset - n_to_plot[pop]
- y_min = offset
- ax.set_xlim([t_min, t_max])
- ax.set_ylim([y_min, y_max])
- ax.set_yticklabels(yticks)
- ax.set_yticks(yticklocs)
- ax.set_xlabel('Time (s)', labelpad=-0.1)
- ax.set_xticks([t_min, t_min + 250., t_max])
- # ax.set_xticklabels([r'$3.$', r'$3.25$', r'$3.5$'])
- l = t_min/1000
- m = (t_min + t_max)/2000
- r = t_max/1000
- ax.set_xticklabels([f'{l:.2f}', f'{m:.2f}', f'{r:.2f}'])
-
- # print("plotting Population rates")
-
- rates = np.zeros((len(M.area_list), 8))
- for i, area in enumerate(M.area_list):
- for j, pop in enumerate(M.structure[area][::-1]):
- rate = pop_rates[area][pop][0]
- if rate == 0.0:
- rate = 1e-5
- if area == 'TH' and j > 3: # To account for missing layer 4 in TH
- rates[i][j + 2] = rate
- else:
- rates[i][j] = rate
-
-
- rates = np.transpose(rates)
- masked_rates = np.ma.masked_where(rates < 1e-4, rates)
-
- ax = axes['D']
- d = ax.boxplot(np.transpose(rates), vert=False,
- patch_artist=True, whis=1.5, showmeans=True)
- set_boxplot_props(d)
-
- ax.plot(np.mean(rates, axis=1), np.arange(
- 1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
- ax.set_yticklabels(population_labels[::-1], size=8)
- ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
- ax.set_ylim((0., len(M.structure['V1']) + .5))
-
- x_max = 220.
- ax.set_xlim((-1., x_max))
- ax.set_xlabel(r'Rate (spikes/s)', labelpad=-0.1)
- ax.set_xticks([0., 50., 100.])
-
- # print("plotting Synchrony")
-
- syn = np.zeros((len(M.area_list), 8))
- for i, area in enumerate(M.area_list):
- for j, pop in enumerate(M.structure[area][::-1]):
- value = corrcoeff[area][pop]
- if value == 0.0:
- value = 1e-5
- if area == 'TH' and j > 3: # To account for missing layer 4 in TH
- syn[i][j + 2] = value
- else:
- syn[i][j] = value
-
-
- syn = np.transpose(syn)
- masked_syn = np.ma.masked_where(syn < 1e-4, syn)
-
- ax = axes['E']
- d = ax.boxplot(np.transpose(syn), vert=False,
- patch_artist=True, whis=1.5, showmeans=True)
- set_boxplot_props(d)
-
- ax.plot(np.mean(syn, axis=1), np.arange(
- 1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
-
- ax.set_yticklabels(population_labels[::-1], size=8)
- ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
- ax.set_ylim((0., len(M.structure['V1']) + .5))
- # ax.set_xticks(np.arange(0.0, 0.601, 0.2))
- ax.set_xticks(np.arange(0.0, 10., 2.0))
- # ax.set_xlabel('Correlation coefficient', labelpad=-0.1)
- ax.set_xlabel('Synchrony', labelpad=-0.1)
-
-
- # print("plotting Irregularity")
-
- LvR = np.zeros((len(M.area_list), 8))
- for i, area in enumerate(M.area_list):
- for j, pop in enumerate(M.structure[area][::-1]):
- value = pop_LvR[area][pop]
- if value == 0.0:
- value = 1e-5
- if area == 'TH' and j > 3: # To account for missing layer 4 in TH
- LvR[i][j + 2] = value
- else:
- LvR[i][j] = value
-
- LvR = np.transpose(LvR)
- masked_LvR = np.ma.masked_where(LvR < 1e-4, LvR)
-
- ax = axes['F']
- d = ax.boxplot(np.transpose(LvR), vert=False,
- patch_artist=True, whis=1.5, showmeans=True)
- set_boxplot_props(d)
-
- ax.plot(np.mean(LvR, axis=1), np.arange(
- 1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
- ax.set_yticklabels(population_labels[::-1], size=8)
- ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
- ax.set_ylim((0., len(M.structure['V1']) + .5))
-
-
- x_max = 2.9
- ax.set_xlim((0., x_max))
- ax.set_xlabel('Irregularity', labelpad=-0.1)
- ax.set_xticks([0., 1., 2.])
-
- axes['G'].spines['right'].set_color('none')
- axes['G'].spines['left'].set_color('none')
- axes['G'].spines['top'].set_color('none')
- axes['G'].spines['bottom'].set_color('none')
- axes['G'].yaxis.set_ticks_position("none")
- axes['G'].xaxis.set_ticks_position("none")
- axes['G'].set_xticks([])
- axes['G'].set_yticks([])
-
-
- # print("Plotting rate time series")
- pos = axes['G'].get_position()
- ax = []
- h = pos.y1 - pos.y0
- w = pos.x1 - pos.x0
- ax.append(pl.axes([pos.x0, pos.y0, w, 0.28 * h]))
- ax.append(pl.axes([pos.x0, pos.y0 + 0.33 * h, w, 0.28 * h]))
- ax.append(pl.axes([pos.x0, pos.y0 + 0.67 * h, w, 0.28 * h]))
-
- colors = ['0.5', '0.3', '0.0']
-
- # t_min = 500.
- # t_max = 10500.
- t_min = 500.
- t_max = t_sim
- # time = np.arange(500., t_max)
- time = np.arange(500, t_max)
- for i, area in enumerate(areas[::-1]):
- ax[i].spines['right'].set_color('none')
- ax[i].spines['top'].set_color('none')
- ax[i].yaxis.set_ticks_position("left")
- ax[i].xaxis.set_ticks_position("none")
-
- binned_spikes = rate_time_series[area][np.where(
- np.logical_and(time >= t_min, time < t_max))]
- ax[i].plot(time, binned_spikes, color=colors[0], label=area)
- # rate = rate_time_series_auto_kernel[area]
- # ax[i].plot(time, rate, color=colors[2], label=area)
- ax[i].set_xlim((t_min, t_max))
-
- ax[i].text(0.8, 0.7, area, transform=ax[i].transAxes)
-
- if i > 0:
- ax[i].spines['bottom'].set_color('none')
- ax[i].set_xticks([])
- ax[i].set_yticks([0., 30.])
- else:
- # ax[i].set_xticks([1000., 5000., 10000.])
- ax[i].set_xticks([t_min, (t_min+t_max)/2, t_max])
- l = t_min/1000
- m = (t_min + t_max)/2000
- r = t_max/1000
- # ax[i].set_xticklabels([r'$1.$', r'$5.$', r'$10.$'])
- ax[i].set_xticklabels([f'{l:.2f}', f'{m:.2f}', f'{r:.2f}'])
- ax[i].set_yticks([0., 5.])
- if i == 1:
- ax[i].set_ylabel(r'Rate (spikes/s)')
-
- ax[0].set_xlabel('Time (s)', labelpad=-0.05)
-
- fig.subplots_adjust(left=0.05, right=0.95, top=0.95,
- bottom=0.075, wspace=1., hspace=.5)
-
- # pl.savefig('Fig5_ground_state.eps')
\ No newline at end of file
diff --git a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_time_ave_pop_rates-checkpoint.py b/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_time_ave_pop_rates-checkpoint.py
deleted file mode 100644
index 89b991acfd40842b4bcb9a32839330f0a2a13a03..0000000000000000000000000000000000000000
--- a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_time_ave_pop_rates-checkpoint.py
+++ /dev/null
@@ -1,91 +0,0 @@
-import numpy as np
-
-from multiarea_model import Analysis
-import matplotlib.pyplot as plt
-from matplotlib.colors import LogNorm
-from matplotlib.ticker import FixedLocator
-
-def plot_time_averaged_population_rates(M, area_list=None, **keywords):
- """
- Plot overview over time-averaged population rates encoded in colors
- with areas along x-axis and populations along y-axis.
-
- Parameters
- ----------
- area_list : list, optional
- Specifies with areas are plotted in which order.
- Default to None, leading to plotting of all areas ordered by architectural type.
- output : {'pdf', 'png', 'eps'}, optional
- If given, the function stores the plot to a file of the given format.
- """
-
- A = Analysis(network=M,
- simulation=M.simulation,
- data_list=['spikes'],
- load_areas=None)
-
- A.create_pop_rates()
-
- if area_list is None:
- area_list = ['V1', 'V2', 'VP', 'V3', 'PIP', 'V3A', 'MT', 'V4t', 'V4',
- 'PO', 'VOT', 'DP', 'MIP', 'MDP', 'MSTd', 'VIP', 'LIP',
- 'PITv', 'PITd', 'AITv', 'MSTl', 'FST', 'CITv', 'CITd',
- '7a', 'STPp', 'STPa', 'FEF', '46', 'TF', 'TH', 'AITd']
-
- matrix = np.zeros((len(area_list), len(A.network.structure['V1'])))
-
- fig = plt.figure(figsize=(12, 4))
- fig.suptitle('Time-averaged population rates encoded in colors', fontsize=15, x=0.43)
- ax = fig.add_subplot(111)
-
- for i, area in enumerate(area_list):
- # print(i, area)
- # for j, pop in enumerate(A.network.structure_reversed['V1']):
- for j, pop in enumerate(A.network.structure['V1'][::-1]):
- if pop in A.network.structure[area]:
- rate = A.pop_rates[area][pop][0]
- if rate == 0.0:
- rate = 1e-5 # To distinguish zero-rate from non-existing populations
- else:
- rate = np.nan
- matrix[i][j] = rate
-
- cm = plt.cm.jet
- cm = cm.from_list('mycmap', [(0., 64./255., 192./255.), # custom dark blue
- (0., 128./255., 192./255.), # custom light blue
- 'white',
- (245./255., 157./255., 115./255.), # custom light red
- (192./255., 64./255., 0.)], N=256) # custom dark red
- cm.set_under('0.3')
- cm.set_bad('k')
-
- matrix = np.transpose(matrix)
- masked_matrix = np.ma.masked_where(np.isnan(matrix), matrix)
- ax.patch.set_hatch('x')
- im = ax.pcolormesh(masked_matrix, cmap=cm, edgecolors='None', norm=LogNorm(
- vmin=0.01, vmax=100.))
- ax.set_xlim(0, matrix[0].size)
-
- x_index = np.arange(4.5, 31.6, 5.0)
- x_ticks = [int(a + 0.5) for a in x_index]
- y_index = list(range(len(A.network.structure['V1'])))
- y_index = [a + 0.5 for a in y_index]
- # print(A.network.structure['V1'])
- # ax.set_xticks(x_index)
- ax.set_xticks([i + 0.5 for i in np.arange(0, len(area_list), 1)])
- # ax.set_xticklabels(x_ticks)
- ax.set_xticklabels(area_list, rotation=90, size=10.)
- ax.set_yticks(y_index)
- # ax.set_yticklabels(A.network.structure_reversed['V1'])
- ax.set_yticklabels(A.network.structure['V1'][::-1])
- ax.set_ylabel('Population', size=13)
- ax.set_xlabel('Area index', size=13)
- t = FixedLocator([0.01, 0.1, 1., 10., 100.])
-
- plt.colorbar(im, ticks=t)
-
- if 'output' in keywords:
- plt.savefig(os.path.join(A.output_dir, '{}_rates.{}'.format(A.simulation.label,
- keywords['output'])))
- else:
- fig.show()
\ No newline at end of file