Merge pull request #46 from didi-hou/Allow-simulating-some-aeas

Allow simulating only some of the areas

figures/MAM2EBRAINS/M2E_compute_corrcoeff.py

@@ -23,6 +23,7 @@ def compute_corrcoeff(M, data_path, label):
         sim_params = json.load(f)
     # T = sim_params['T']
     T = sim_params['sim_params']['t_sim']
+    areas_simulated = sim_params['sim_params']['areas_simulated']
 
     tmin = 500.
     subsample = 2000
@@ -35,7 +36,8 @@ def compute_corrcoeff(M, data_path, label):
 
     spike_data = {}
     cc_dict = {}
-    for area in M.area_list:
+    # for area in M.area_list:
+    for area in areas_simulated:
         cc_dict[area] = {}
         LvR_list = []
         N = []

figures/MAM2EBRAINS/M2E_visualize_resting_state.py

@@ -59,6 +59,7 @@ def plot_resting_state(M, data_path, raster_areas=['V1', 'V2', 'FEF']):
         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, 
@@ -123,11 +124,20 @@ def plot_resting_state(M, data_path, raster_areas=['V1', 'V2', 'FEF']):
                  '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.")
+
+    # if len(raster_areas) !=3:
+    #     raise Exception("Error! Please give 3 areas to display as raster plots.")
+
+    with open(os.path.join(data_path, M.simulation.label, 'custom_params_{}'.format(M.simulation.label)), 'r') as f:
+        sim_params = json.load(f)
+
+    areas_simulated = sim_params['sim_params']['areas_simulated']
+
     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.")
+            raise Exception("Error! Given raster areas are not from complete_area_list, please input correct areas to diaply the raster plots.")
+        if area not in areas_simulated:
+            raise Exception("Error! At least one of the given raster areas are not from the simulated areas, please input correct areas to diaply the raster plots.")
 
     areas = raster_areas
 
@@ -325,8 +335,19 @@ def plot_resting_state(M, data_path, raster_areas=['V1', 'V2', 'FEF']):
 
     # print("plotting Population rates")
 
-    rates = np.zeros((len(M.area_list), 8))
-    for i, area in enumerate(M.area_list):
+    # 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.zeros((len(areas), 8))
+    for i, area in enumerate(areas):
         for j, pop in enumerate(M.structure[area][::-1]):
             rate = pop_rates[area][pop][0]
             if rate == 0.0:
@@ -358,8 +379,19 @@ def plot_resting_state(M, data_path, raster_areas=['V1', 'V2', 'FEF']):
 
     # print("plotting Synchrony")
 
-    syn = np.zeros((len(M.area_list), 8))
-    for i, area in enumerate(M.area_list):
+    # 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.zeros((len(areas), 8))
+    for i, area in enumerate(areas):
         for j, pop in enumerate(M.structure[area][::-1]):
             value = corrcoeff[area][pop]
             if value == 0.0:
@@ -369,7 +401,6 @@ def plot_resting_state(M, data_path, raster_areas=['V1', 'V2', 'FEF']):
             else:
                 syn[i][j] = value
 
-
     syn = np.transpose(syn)
     masked_syn = np.ma.masked_where(syn < 1e-4, syn)
 
@@ -385,15 +416,26 @@ def plot_resting_state(M, data_path, raster_areas=['V1', 'V2', 'FEF']):
     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([0., 0.025, 0.05])
+    ax.set_xticks([0., 0.005, 0.01])
     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):
+    # 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.zeros((len(areas), 8))
+    for i, area in enumerate(areas):
         for j, pop in enumerate(M.structure[area][::-1]):
             value = pop_LvR[area][pop]
             if value == 0.0:

figures/MAM2EBRAINS/M2E_visualize_time_ave_pop_rates.py

@@ -1,11 +1,13 @@
 import numpy as np
+import os
+import json
 
 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):
+def plot_time_averaged_population_rates(M, data_path, area_list=None, **keywords):
     """
     Plot overview over time-averaged population rates encoded in colors
     with areas along x-axis and populations along y-axis.
@@ -31,14 +33,21 @@ def plot_time_averaged_population_rates(M, area_list=None, **keywords):
                      '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'])))
+
+    with open(os.path.join(data_path, M.simulation.label, 'custom_params_{}'.format(M.simulation.label)), 'r') as f:
+        sim_params = json.load(f)
+
+    areas_simulated = sim_params['sim_params']['areas_simulated']
+
+    # matrix = np.zeros((len(area_list), len(A.network.structure['V1'])))
+    matrix = np.zeros((len(areas_simulated), 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):
+    # for i, area in enumerate(area_list):
+    for i, area in enumerate(areas_simulated):
         # print(i, area)
         # for j, pop in enumerate(A.network.structure_reversed['V1']):
         for j, pop in enumerate(A.network.structure['V1'][::-1]):