From e22a394a988b54ab503af221863f8310d74f2670 Mon Sep 17 00:00:00 2001
From: Maximilian Schmidt <max.schmidt@fz-juelich.de>
Date: Thu, 29 Mar 2018 15:20:17 +0900
Subject: [PATCH] Remove dependency from h5py_wrapper completely: Store data in
binary numpy data files, add test_analysis to test analysis class, fix many
bugs in analysis and analysis_helpers
---
README.md | 2 +-
environment.yaml | 4 -
multiarea_model/analysis.py | 198 ++++++++++++++--------------
multiarea_model/analysis_helpers.py | 121 ++++++++++++-----
multiarea_model/default_params.py | 14 +-
multiarea_model/simulation.py | 37 ++++--
multiarea_model/theory.py | 11 +-
tests/test_analysis.py | 33 +++++
tests/test_replace_cc.py | 12 +-
tests/test_simulate_V1.py | 10 +-
10 files changed, 276 insertions(+), 166 deletions(-)
create mode 100644 tests/test_analysis.py
diff --git a/README.md b/README.md
index 669f831..704be89 100644
--- a/README.md
+++ b/README.md
@@ -167,7 +167,7 @@ It can be conveniently run by executing `pytest` in the `tests/` folder:
## Requirements
-h5py\_wrapper, python\_dicthash ([https://github.com/INM-6/python-dicthash](https://github.com/INM-6/python-dicthash)),
+python\_dicthash ([https://github.com/INM-6/python-dicthash](https://github.com/INM-6/python-dicthash)),
correlation\_toolbox ([https://github.com/INM-6/correlation-toolbox](https://github.com/INM-6/correlation-toolbox)),
pandas, numpy, nested_dict, matplotlib (2.1.2), pyx, scipy, NEST 2.14.0
diff --git a/environment.yaml b/environment.yaml
index e4f6bdf..46ef9cf 100644
--- a/environment.yaml
+++ b/environment.yaml
@@ -1,16 +1,12 @@
name: multiarea_model
dependencies:
- - future
- numpy
- matplotlib
- pandas
- python==3.6
- scipy
- pip:
- - h5py
- - h5py_wrapper
- nested_dict
- - requests
- https://github.com/INM-6/python-dicthash/archive/master.zip
- pyx
diff --git a/multiarea_model/analysis.py b/multiarea_model/analysis.py
index 035a0b1..ea4b907 100644
--- a/multiarea_model/analysis.py
+++ b/multiarea_model/analysis.py
@@ -11,7 +11,7 @@ multi-area model of macaque visual cortex (Schmidt et al. 2017).
Classes
--------
-analysis : loads the data of the specified simulation and provides members
+Analysis : loads the data of the specified simulation and provides members
functions to post-process the data and plot it in various visualizations.
Authors
@@ -22,7 +22,6 @@ Sacha van Albada
"""
from . import analysis_helpers as ah
import glob
-import h5py_wrapper.wrapper as h5w
import inspect
import json
import matplotlib.pyplot as plt
@@ -30,6 +29,7 @@ import numpy as np
import os
import pandas as pd
+from copy import copy
from matplotlib.colors import LogNorm
from matplotlib.ticker import FixedLocator
from nested_dict import nested_dict
@@ -101,38 +101,35 @@ class Analysis:
rec_dir = os.path.join(self.simulation.data_dir, 'recordings')
for data_type in data_list:
if data_type == 'spikes':
- h5_file = os.path.join(rec_dir, '_'.join((self.simulation.label,
- 'spike_data.h5')))
- d = 'spike_dict'
columns = ['senders', 'times']
+ d = 'spike_dict'
elif data_type == 'vm':
assert(self.simulation.params['recording_dict']['record_vm']), "Trying to "
"load membrane potentials, but these data have not been recorded"
- h5_file = os.path.join(rec_dir, '_'.join((self.simulation.label,
- 'vm_data.h5')))
d = 'vm_dict'
columns = ['senders', 'times', 'V_m']
print('loading {}'.format(data_type))
data = {}
- # Check if the data has been merged into one common file.
- try:
- lazy_load = h5w.load(h5_file, lazy=True)
- except IOError:
- lazy_load = {}
- h5w.save(h5_file, {})
+ # Check if the data has already been stored in binary file
for area in self.areas_loaded:
- if area in lazy_load:
- data[area] = h5w.load(h5_file, path=area)
- else:
- area_dat = {}
+ data[area] = {}
+ fp = '-'.join((self.simulation.label,
+ self.simulation.params['recording_dict'][d]['label'],
+ area,
+ '*'))
+ files = glob.glob(os.path.join(rec_dir, fp))
+ if len(files) == 0:
+ continue
+ for pop in self.network.structure[area]:
fp = '-'.join((self.simulation.label,
self.simulation.params['recording_dict'][d]['label'],
area,
- '*'))
- files = glob.glob(os.path.join(rec_dir, fp))
- if len(files) == 0:
- continue
- for pop in self.network.structure[area]:
+ pop))
+ fn = os.path.join(rec_dir,
+ '.'.join((fp, 'npy')))
+ try:
+ data[area][pop] = np.load(fn)
+ except FileNotFoundError:
fp = '-'.join((self.simulation.label,
self.simulation.params['recording_dict'][d]['label'],
area,
@@ -145,10 +142,9 @@ class Analysis:
names=columns, sep='\t',
index_col=False),
ignore_index=True)
- area_dat[pop] = np.array(dat)
- data[area] = area_dat
+ np.save(fn, np.array(dat))
+ data[area][pop] = np.array(dat)
# Store spike data to single hdf5 file
- h5w.save(h5_file, data[area], path=area, write_mode='a')
if data_type == 'spikes':
self.spike_data = data
elif data_type == 'vm':
@@ -206,9 +202,13 @@ class Analysis:
default_dict = {'areas': self.areas_loaded,
'pops': 'complete', 'compute_stat': False}
params = ah._create_parameter_dict(default_dict, self.T, **keywords)
+ iterator = ah.model_iter(mode='single',
+ areas=params['areas'],
+ pops=params['pops'])
# Check if population rates have been stored with the same parameters
- self.pop_rates = ah._check_stored_data(os.path.join(self.output_dir, 'pop_rates.json'),
- params)
+ fp = os.path.join(self.output_dir, 'pop_rates.json')
+ self.pop_rates = ah._check_stored_data(fp,
+ copy(iterator), params)
if self.pop_rates is None:
print("Computing population rates")
@@ -231,9 +231,7 @@ class Analysis:
total_rates += rate[2]
d[area]['total'] = (np.mean(total_rates), np.std(total_rates))
else:
- for area, pop in ah.model_iter(mode='single',
- areas=params['areas'],
- pops=params['pops']):
+ for area, pop in iterator:
if pop in self.network.structure[area]:
spikes = self.spike_data[area][pop][:, 1]
indices = np.where(np.logical_and(spikes > params['t_min'],
@@ -256,7 +254,7 @@ class Analysis:
Calculate single neuron population rates and store them in member pop_rate_dists.
If the distributions had previously been stored with the
same parameters, they are loaded from file.
- Relies on helper function pop_rate_distribution.
+ Uses helper function pop_rate_distribution.
Parameters
----------
@@ -278,18 +276,20 @@ class Analysis:
default_dict = {'areas': self.areas_loaded, 'pops': 'complete'}
params = ah._create_parameter_dict(
default_dict, self.T, **keywords)
+ iterator = ah.model_iter(mode='single',
+ areas=params['areas'],
+ pops=params['pops'])
# Check if population rates have been stored with the same parameters
self.pop_rate_dists = ah._check_stored_data(os.path.join(self.output_dir,
- 'pop_rate_dists.h5'), params)
+ 'pop_rate_dists'),
+ copy(iterator), params)
if self.pop_rate_dists is None:
print("Computing population dists")
d = nested_dict()
d['Parameters'] = params
- for area, pop in ah.model_iter(mode='single',
- areas=params['areas'],
- pops=params['pops']):
+ for area, pop in iterator:
if pop in self.network.structure[area]:
res = list(ah.pop_rate_distribution(self.spike_data[area][pop],
params['t_min'],
@@ -303,7 +303,7 @@ class Analysis:
def create_synchrony(self, **keywords):
"""
Calculate synchrony as the coefficient of variation of the population rate
- and store in member synchrony. Relies on helper function synchrony.
+ 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.
@@ -331,19 +331,19 @@ class Analysis:
'pops': 'complete', 'resolution': 1.0}
params = ah._create_parameter_dict(
default_dict, self.T, **keywords)
-
+ iterator = ah.model_iter(mode='single',
+ areas=params['areas'],
+ pops=params['pops'])
# Check if synchrony values have been stored with the same parameters
self.synchrony = ah._check_stored_data(os.path.join(self.output_dir, 'synchrony.json'),
- params)
+ copy(iterator), params)
if self.synchrony is None:
print("Computing synchrony")
d = nested_dict()
d['Parameters'] = params
- for area, pop in ah.model_iter(mode='single',
- areas=params['areas'],
- pops=params['pops']):
- if pop in self.network.structure[area] > 0.:
+ for area, pop in iterator:
+ if pop in self.network.structure[area]:
d[area][pop] = ah.synchrony(self.spike_data[area][pop],
self.network.N[area][pop],
params['t_min'],
@@ -365,7 +365,7 @@ class Analysis:
def create_rate_time_series(self, **keywords):
"""
Calculate time series of population- and area-averaged firing rates.
- Relies on ah.pop_rate_time_series.
+ Uses ah.pop_rate_time_series.
If the rates have previously been stored with the
same parameters, they are loaded from file.
@@ -401,10 +401,16 @@ class Analysis:
default_dict, self.T, **keywords)
# Check if firing rates have been stored with the same parameters
- fn = os.path.join(self.output_dir, 'rate_time_series.h5')
- self.rate_time_series = ah._check_stored_data(fn, params)
- fn = os.path.join(self.output_dir, 'rate_time_series_pops.h5')
- self.rate_time_series_pops = ah._check_stored_data(fn, params)
+ fp = os.path.join(self.output_dir, 'rate_time_series')
+ iterator_areas = ah.model_iter(mode='single',
+ areas=params['areas'],
+ pops=None)
+ iterator_pops = ah.model_iter(mode='single',
+ areas=params['areas'],
+ pops=params['pops'])
+ self.rate_time_series = ah._check_stored_data(fp, copy(iterator_areas), params)
+ fp = os.path.join(self.output_dir, 'rate_time_series_pops')
+ self.rate_time_series_pops = ah._check_stored_data(fp, copy(iterator_pops), params)
if self.rate_time_series is None:
print('Computing rate time series')
@@ -415,10 +421,7 @@ class Analysis:
# population-averaged firing rates
d_pops = nested_dict()
d_pops['Parameters'] = params
-
- for area, pop in ah.model_iter(mode='single',
- areas=params['areas'],
- pops=params['pops']):
+ for area, pop in iterator_pops:
if pop in self.network.structure[area]:
time_series = ah.pop_rate_time_series(self.spike_data[area][pop],
self.network.N[area][pop],
@@ -444,7 +447,7 @@ class Analysis:
def create_synaptic_input(self, **keywords):
"""
Calculate synaptic input of populations and areas using the spike data.
- Relies on function ah.pop_synaptic_input.
+ Uses function ah.pop_synaptic_input.
If the synaptic inputs have previously been stored with the
same parameters, they are loaded from file.
@@ -479,10 +482,16 @@ class Analysis:
default_dict, self.T, **keywords)
# Check if synaptic inputs have been stored with the same parameters
- fn = os.path.join(self.output_dir, 'synaptic_input.h5')
- self.synaptic_input = ah._check_stored_data(fn, params)
- fn = os.path.join(self.output_dir, 'synaptic_input_pops.h5')
- self.synaptic_input_pops = ah._check_stored_data(fn, params)
+ iterator_areas = ah.model_iter(mode='single',
+ areas=params['areas'],
+ pops=None)
+ iterator_pops = ah.model_iter(mode='single',
+ areas=params['areas'],
+ pops=params['pops'])
+ fp = os.path.join(self.output_dir, 'synaptic_input')
+ self.synaptic_input = ah._check_stored_data(fp, copy(iterator_areas), params)
+ fp = os.path.join(self.output_dir, 'synaptic_input_pops')
+ self.synaptic_input_pops = ah._check_stored_data(fp, copy(iterator_pops), params)
if self.synaptic_input is None:
print('Computing rate time series')
@@ -491,32 +500,27 @@ class Analysis:
d_pops = nested_dict()
d_pops['Parameters'] = params
- for area, pop in ah.model_iter(mode='single',
- areas=params['areas'],
- pops=params['pops']):
+ for area, pop in copy(iterator_pops):
if pop in self.network.structure[area]:
if 'I' in pop:
- tau_syn = self.network.params['neuron_params']['single_neuron_dict']['tau_syn_in']
+ tau_syn = self.network.params['neuron_params'][
+ 'single_neuron_dict']['tau_syn_in']
else:
- tau_syn = self.network.params['neuron_params']['single_neuron_dict']['tau_syn_ex']
-# if pop in self.rate_time_series_pops[area]:
+ tau_syn = self.network.params['neuron_params'][
+ 'single_neuron_dict']['tau_syn_ex']
time_series = ah.synaptic_output(self.rate_time_series_pops[area][pop],
tau_syn, params['t_min'], params['t_max'],
resolution=params['resolution'])
- # else:
- # time_series = np.ones(params['t_max'] - params['t_min']) * np.nan
d_pops[area][pop] = time_series
self.synaptic_output_pops = d_pops.to_dict()
d_pops = nested_dict()
d_pops['Parameters'] = params
d_pops['Parameters'] = params
- for area, pop in ah.model_iter(mode='single',
- areas=params['areas'],
- pops=params['pops']):
+ for area, pop in iterator_pops:
if pop in self.network.structure[area]:
time_series = np.zeros(
- (params['t_max'] - params['t_min']) / params['resolution'])
+ int((params['t_max'] - params['t_min']) / params['resolution']))
for source_area, source_pop in ah.model_iter(mode='single',
areas=self.areas_loaded):
if source_pop in self.network.structure[source_area]:
@@ -530,7 +534,7 @@ class Analysis:
d['Parameters'] = params
for area in params['areas']:
d[area] = np.zeros(
- (params['t_max'] - params['t_min']) / params['resolution'])
+ int((params['t_max'] - params['t_min']) / params['resolution']))
for pop in self.network.structure[area]:
d[area] += d_pops[area][pop] * self.network.N[area][pop]
d[area] /= self.network.N[area]['total']
@@ -540,7 +544,7 @@ class Analysis:
def create_pop_cv_isi(self, **keywords):
"""
Calculate population-averaged CV ISI values and store as member pop_cv_isi.
- Relies on helper function cv_isi.
+ Uses helper function cv_isi.
If the CV ISI have previously been stored with the
same parameters, they are loaded from file.
@@ -565,17 +569,19 @@ class Analysis:
params = ah._create_parameter_dict(
default_dict, self.T, **keywords)
# Check if CV ISI have been stored with the same parameters
- self.pop_cv_isi = ah._check_stored_data(os.path.join(self.output_dir, 'pop_cv_isi.json'),
- params)
+ iterator = ah.model_iter(mode='single',
+ areas=params['areas'],
+ pops=params['pops'])
+ fp = os.path.join(self.output_dir, 'pop_cv_isi.json')
+ self.pop_cv_isi = ah._check_stored_data(fp,
+ copy(iterator), params)
if self.pop_cv_isi is None:
print("Computing population CV ISI")
d = nested_dict()
d['Parameters'] = params
- for area, pop in ah.model_iter(mode='single',
- areas=params['areas'],
- pops=params['pops']):
- if pop in self.network.structure[area] > 0.:
+ for area, pop in iterator:
+ if pop in self.network.structure[area]:
d[area][pop] = ah.pop_cv_isi(self.spike_data[area][pop],
params['t_min'],
params['t_max'])
@@ -584,7 +590,7 @@ class Analysis:
def create_pop_LvR(self, **keywords):
"""
Calculate poulation-averaged LvR (see Shinomoto et al. 2009) and
- store as member pop_LvR. Relies on helper function LvR.
+ store as member pop_LvR. Uses helper function LvR.
Parameters
----------
@@ -607,22 +613,24 @@ class Analysis:
default_dict, self.T, **keywords)
# Check if LvR have been stored with the same parameters
- self.pop_LvR = ah._check_stored_data(os.path.join(self.output_dir, 'pop_LvR.json'),
- params)
+ iterator = ah.model_iter(mode='single',
+ areas=params['areas'],
+ pops=params['pops'])
+ fp = os.path.join(self.output_dir, 'pop_LvR.json')
+ self.pop_LvR = ah._check_stored_data(fp,
+ copy(iterator), params)
if self.pop_LvR is None:
print("Computing population LvR")
d = nested_dict()
d['Parameters'] = params
- for area, pop in ah.model_iter(mode='single',
- areas=params['areas'],
- pops=params['pops']):
- if pop in self.network.structure[area] > 0.:
+ for area, pop in iterator:
+ if pop in self.network.structure[area]:
if self.network.N[area][pop] > 0.:
d[area][pop] = ah.pop_LvR(self.spike_data[area][pop],
2.0,
params['t_min'],
params['t_max'])[0]
- self.pop_LvR = d.to_dict()
+ self.pop_LvR = d.to_dict()
# ______________________________________________________________________________
# Function for plotting data
@@ -839,11 +847,11 @@ class Analysis:
if pop:
plt.savefig(os.path.join(self.output_dir,
'{}_power_spectrum_{}_{}.{}'.format(self.simulation.label,
- area, pop, keywords['output'])))
+ area, pop, keywords['output'])))
else:
plt.savefig(os.path.join(self.output_dir,
'{}_power_spectrum_{}.{}'.format(self.simulation.label,
- area, keywords['output'])))
+ area, keywords['output'])))
else:
fig.show()
@@ -928,19 +936,17 @@ class Analysis:
"""
members = inspect.getmembers(self)
save_list_json = ['structure', 'pop_rates', 'synchrony',
- 'func_conn', 'func_conn_pops', 'pop_cv_isi',
+ 'pop_cv_isi', 'pop_LvR',
'indegree_data', 'indegree_areas_data',
- 'outdegree_data', 'outdegree_areas_data',
- 'pop_LvR']
- save_list_hdf5 = ['activity_flux', 'activity_flux_pops',
- 'pop_rate_dists', 'rate_time_series',
- 'rate_time_series_pops', 'bold_signal',
- 'synaptic_input', 'synaptic_input_pops']
+ 'outdegree_data', 'outdegree_areas_data']
+ save_list_npy = ['pop_rate_dists', 'rate_time_series',
+ 'rate_time_series_pops', 'bold_signal',
+ 'synaptic_input', 'synaptic_input_pops']
for i in range(0, len(members)):
if members[i][0] in save_list_json:
f = open(self.output_dir + members[i][0] + '.json', 'w')
json.dump(members[i][1], f)
f.close()
- if members[i][0] in save_list_hdf5:
- f = self.output_dir + members[i][0] + '.h5'
- h5w.save(f, members[i][1], write_mode='w')
+ if members[i][0] in save_list_npy:
+ f = self.output_dir + members[i][0]
+ ah._save_dict_to_npy(f, members[i][1])
diff --git a/multiarea_model/analysis_helpers.py b/multiarea_model/analysis_helpers.py
index fc7b9b2..f8f3eaa 100644
--- a/multiarea_model/analysis_helpers.py
+++ b/multiarea_model/analysis_helpers.py
@@ -38,13 +38,12 @@ Sacha van Albada
"""
-import copy
+from copy import copy, deepcopy
+from nested_dict import nested_dict
import numpy as np
-import h5py_wrapper.wrapper as h5w
import json
from itertools import product
from scipy.signal import welch
-#################################
area_list = ['V1', 'V2', 'VP', 'V3', 'PIP', 'V3A', 'MT', 'V4t', 'V4',
@@ -73,8 +72,9 @@ def model_iter(mode='single', areas=None, pops='complete', areas2=None, pops2='c
Defaults to None, which correspons to taking all areas into account.
pops, pops2 : string or list, optional
If specified, loop only over these populations as target and source
- populations. Defaults to 'complete', which correspons to taking all
+ populations. Defaults to 'complete', which corresponds to taking all
areas into account.
+ If None, loop only over areas.
Returns
-------
@@ -83,6 +83,8 @@ def model_iter(mode='single', areas=None, pops='complete', areas2=None, pops2='c
"""
if mode == 'single':
assert((areas2 is None) and (pops2 is 'complete'))
+ if pops is None or pops2 is None:
+ assert((pops is None) and (pops2 is None) or mode == 'single')
if pops == 'complete':
pops = pop_list
if areas is None:
@@ -92,9 +94,15 @@ def model_iter(mode='single', areas=None, pops='complete', areas2=None, pops2='c
if areas2 is None:
areas2 = area_list
if mode == 'single':
- return product(areas, pops)
+ if pops is None:
+ return product(areas)
+ else:
+ return product(areas, pops)
elif mode == 'pairs':
- return product(areas, pops, areas2, pops2)
+ if pops is None:
+ return product(areas, areas2)
+ else:
+ return product(areas, pops, areas2, pops2)
def area_spike_train(spike_data):
@@ -117,7 +125,7 @@ def area_spike_train(spike_data):
data_array = np.array([])
for pop in spike_data:
data_array = np.append(data_array, spike_data[pop])
- data_array = np.reshape(data_array, (len(data_array) * 0.5, 2))
+ data_array = np.reshape(data_array, (-1, 2))
return data_array
@@ -131,7 +139,7 @@ def centralize(data, time=False, units=False):
"""
assert(time is not False or units is not False)
- res = copy.copy(data)
+ res = copy(data)
if time is True:
res = np.array([x - np.mean(x) for x in res])
if units is True:
@@ -179,8 +187,10 @@ def sort_gdf_by_id(data, idmin=None, idmax=None):
print('CT warning(sort_spiketrains_by_id): empty gdf data!')
return None, None
-# ______________________________________________________________________________
-# Helper functions for data loading
+
+"""
+Helper functions for data loading
+"""
def _create_parameter_dict(default_dict, T, **keywords):
@@ -210,7 +220,7 @@ def _create_parameter_dict(default_dict, T, **keywords):
return d
-def _check_stored_data(fn, param_dict):
+def _check_stored_data(fp, fn_iter, param_dict):
"""
Check if a data member of the data class has already
been computed with the same parameters.
@@ -223,21 +233,23 @@ def _check_stored_data(fn, param_dict):
Parameters of the calculation to compare with
the parameters of the stored data.
"""
- if 'json' in fn:
+ if 'json' in fp:
try:
- f = open(fn)
+ f = open(fp)
data = json.load(f)
f.close()
except IOError:
return None
- elif 'hdf5' in fn:
+ param_dict2 = data['Parameters']
+ else:
try:
- data = h5w.load(fn)
+ data = _load_npy_to_dict(fp, fn_iter)
except IOError:
return None
- param_dict2 = data['Parameters']
- param_dict_copy = copy.copy(param_dict)
- param_dict2_copy = copy.copy(param_dict2)
+ with open('-'.join((fp, 'parameters')), 'r') as f:
+ param_dict2 = json.load(f)
+ param_dict_copy = copy(param_dict)
+ param_dict2_copy = copy(param_dict2)
for k in param_dict:
if isinstance(param_dict_copy[k], list) or isinstance(param_dict_copy[k], np.ndarray):
param_dict_copy[k] = set(param_dict_copy[k])
@@ -251,8 +263,59 @@ def _check_stored_data(fn, param_dict):
"with different parameters")
return None
-# ______________________________________________________________________________
-# Analysis functions
+
+def _save_dict_to_npy(fp, data):
+ """
+ Save data dictionary to binary numpy files
+ by iteratively going through the dictionary.
+
+ Parameters
+ ----------
+ fp : str
+ File pattern to which the keys of the dictionary are attached.
+ data : dict
+ Dictionary containing the data
+ """
+ for key, val in data.items():
+ if key != 'Parameters':
+ fp_key = '-'.join((fp, key))
+ if isinstance(val, dict):
+ _save_dict_to_npy(fp_key, val)
+ else:
+ np.save(fp_key, val)
+ else:
+ fp_key = '-'.join((fp, 'parameters'))
+ with open(fp_key, 'w') as f:
+ json.dump(val, f)
+
+
+def _load_npy_to_dict(fp, fn_iter):
+ """
+ Load data stored in the files defined by fp
+ and fn_iter to a dictionary.
+
+ Parameters
+ ----------
+ fp : str
+ Base file pattern of the npy files
+ fn_iter : iterable
+ Iterable defining all the suffixes that are
+ appended to fp to form the file names.
+ """
+ data = nested_dict()
+ for it in fn_iter:
+ fp_it = (fp,) + it
+ fp_ = '{}.npy'.format('-'.join(fp_it))
+ if len(it) == 1:
+ data[it[0]] = np.load(fp_)
+ else:
+ data[it[0]][it[1]] = np.load(fp_)
+ return data
+
+
+"""
+Analysis functions
+"""
def online_hist(fname, tmin, tmax, resolution=1.):
@@ -383,8 +446,11 @@ def pop_rate_distribution(data_array, t_min, t_max, num_neur):
data_array[:, 1] < t_max))
neurons = data_array[:, 0][indices]
neurons = np.sort(neurons)
- n = neurons[0]
- rates = np.zeros(num_neur)
+ if len(neurons) > 0:
+ n = neurons[0]
+ else: # No spikes in [t_min, t_max]
+ n = None
+ rates = np.zeros(int(num_neur))
ii = 0
for i in range(neurons.size):
if neurons[i] == n:
@@ -474,8 +540,6 @@ def pop_rate_time_series(data_array, num_neur, t_min, t_max,
return time_series
-# Regularity measures
-
def pop_cv_isi(data_array, t_min, t_max):
"""
@@ -537,7 +601,7 @@ def ISI_SCC(data_array, t_min, t_max):
-------
bins : numpy.ndarray
ISI lags
- valyes : numpy.ndarray
+ values : numpy.ndarray
Serial correlation coefficient values
"""
indices = np.where(np.logical_and(data_array[:, 1] > t_min,
@@ -546,9 +610,8 @@ def ISI_SCC(data_array, t_min, t_max):
half = max(1000, 2 * (t_max - t_min)) / 2.0
if len(data_array) > 1 and len(indices) > 1:
for i in np.unique(data_array[:, 0]):
- print(i)
intervals = np.diff(data_array[indices][
- mlab.find(data_array[indices, 0] == i), 1])
+ np.where(data_array[indices, 0] == i), 1])
if intervals.size > 1:
mean = np.mean(intervals)
@@ -564,7 +627,6 @@ def ISI_SCC(data_array, t_min, t_max):
return 0.0
-
def pop_LvR(data_array, t_ref, t_min, t_max):
"""
Compute the LvR value of the given data_array.
@@ -594,7 +656,7 @@ def pop_LvR(data_array, t_ref, t_min, t_max):
LvR = np.array([])
for i in np.unique(data_array[:, 0]):
intervals = np.diff(data_array[i_min:i_max][
- mlab.find(data_array[i_min:i_max, 0] == i), 1])
+ np.where(data_array[i_min:i_max, 0] == i), 1])
if intervals.size > 1:
val = np.sum((1. - 4 * intervals[0:-1] * intervals[1:] / (intervals[0:-1] + intervals[
1:]) ** 2) * (1 + 4 * t_ref / (intervals[0:-1] + intervals[1:])))
@@ -604,7 +666,6 @@ def pop_LvR(data_array, t_ref, t_min, t_max):
return np.mean(LvR), LvR
-# Synchrony measures
def synchrony(data_array, num_neur, t_min, t_max, resolution=1.0):
"""
Compute the synchrony of an array of spikes as the coefficient
diff --git a/multiarea_model/default_params.py b/multiarea_model/default_params.py
index ccca642..d76b505 100644
--- a/multiarea_model/default_params.py
+++ b/multiarea_model/default_params.py
@@ -129,7 +129,7 @@ connection_params = {
# Poisson input with population-specific rates (het_poisson_stat)
# or by time-varying current input (het_current_nonstat)
# while still simulating all areas. In both cases, the data to replace
- # the cortico-cortical input is loaded from `replace_cc_input_file`.
+ # the cortico-cortical input is loaded from `replace_cc_input_source`.
'replace_cc': False,
# Whether to replace non-simulated areas by poisson sources
@@ -137,11 +137,17 @@ connection_params = {
# by specific rates ('het_poisson_stat')
# or by time-varying specific current ('het_current_nonstat')
# In the two latter cases, the data to replace the cortico-cortical
- # input is loaded from `replace_cc_input_file`
+ # input is loaded from `replace_cc_input_source`
'replace_non_simulated_areas': None,
- # File holding the input rates to replace cortico-cortical input
- 'replace_cc_input_file': None,
+ # Source of the input rates to replace cortico-cortical input
+ # Either a json file (has to end on .json) holding a scalar values
+ # for each population or
+ # a base name such that files with names
+ # $(replace_cc_input_source)-area-population.npy
+ # (e.g. '$(replace_cc_input_source)-V1-23E.npy')
+ # contain the time series for each population.
+ 'replace_cc_input_source': None,
# whether to redistribute CC synapse to meet literature value
# of E-specificity
diff --git a/multiarea_model/simulation.py b/multiarea_model/simulation.py
index 6417ee2..d0b3d47 100644
--- a/multiarea_model/simulation.py
+++ b/multiarea_model/simulation.py
@@ -8,9 +8,10 @@ Schmidt et al. (2018).
Classes
-------
-simulation : Loads a parameter file that specifies simulation parameters for a
+Simulation : Loads a parameter file that specifies simulation parameters for a
simulation of the instance of the model. A simulation is identified by a unique
hash label.
+
"""
import h5py_wrapper as h5w
@@ -22,6 +23,7 @@ import pprint
import shutil
import time
+from .analysis_helpers import _load_npy_to_dict, model_iter
from config import base_path, data_path
from copy import copy
from .default_params import nested_update, sim_params
@@ -123,8 +125,16 @@ class Simulation:
''.join(('custom_Data_Model_', self.network.label, '.json'))),
os.path.join('config_files',
'_'.join((self.network.label, 'config')))]
- if self.network.params['connection_params']['replace_cc_input_file'] is not None:
- files.append(self.network.params['connection_params']['replace_cc_input_file'])
+ if self.network.params['connection_params']['replace_cc_input_source'] is not None:
+ fs = self.network.params['connection_params']['replace_cc_input_source']
+ if '.json' in fs:
+ files.append(fs)
+ else: # Assume that the cc input is stored in one npy file per population
+ fn_iter = model_iter(mode='single', areas=self.network.area_list)
+ for it in fn_iter:
+ fp_it = (fs,) + it
+ fp_ = '{}.npy'.format('-'.join(fp_it))
+ files.append(fp_)
for f in files:
shutil.copy2(os.path.join(base_path, f),
self.data_dir)
@@ -170,18 +180,19 @@ class Simulation:
replace_cc = self.network.params['connection_params']['replace_cc']
replace_non_simulated_areas = self.network.params['connection_params'][
'replace_non_simulated_areas']
- if self.network.params['connection_params']['replace_cc_input_file'] is None:
- replace_cc_input_file = None
+ if self.network.params['connection_params']['replace_cc_input_source'] is None:
+ replace_cc_input_source = None
else:
- replace_cc_input_file = os.path.join(self.data_dir,
- self.network.params['connection_params'][
- 'replace_cc_input_file'])
+ replace_cc_input_source = os.path.join(self.data_dir,
+ self.network.params['connection_params'][
+ 'replace_cc_input_source'])
if not replace_cc and set(self.areas_simulated) != set(self.network.area_list):
if replace_non_simulated_areas == 'het_current_nonstat':
- non_simulated_cc_input = h5w.load(replace_cc_input_file)
+ fn_iter = model_iter(mode='single', areas=self.network.area_list)
+ non_simulated_cc_input = _load_npy_to_dict(replace_cc_input_source, fn_iter)
elif replace_non_simulated_areas == 'het_poisson_stat':
- fn = self.network.params['connection_params']['replace_cc_input_file']
+ fn = self.network.params['connection_params']['replace_cc_input_source']
with open(fn, 'r') as f:
non_simulated_cc_input = json.load(f)
elif replace_non_simulated_areas == 'hom_poisson_stat':
@@ -196,9 +207,11 @@ class Simulation:
" replace non-simulated areas.")
if replace_cc == 'het_current_nonstat':
- cc_input = h5w.load(replace_cc_input_file)
+ fn_iter = model_iter(mode='single', areas=self.network.area_list)
+ cc_input = _load_npy_to_dict(replace_cc_input_source, fn_iter)
elif replace_cc == 'het_poisson_stat':
- with open(self.network.params['connection_params']['replace_cc_input_file'], 'r') as f:
+ with open(self.network.params['connection_params'][
+ 'replace_cc_input_source'], 'r') as f:
cc_input = json.load(f)
elif replace_cc == 'hom_poisson_stat':
cc_input = {source_area_name:
diff --git a/multiarea_model/theory.py b/multiarea_model/theory.py
index b028421..e6ef7ac 100644
--- a/multiarea_model/theory.py
+++ b/multiarea_model/theory.py
@@ -10,16 +10,10 @@ Schmidt et al., 2017) to the network connectivity.
Classes
--------
-
-theory : provides functionality to predict the stable fixed point of
+Theory : provides functionality to predict the stable fixed point of
the model, perform further analysis and execute the stabilization
provedure.
-Authors
---------
-Maximilian Schmidt
-Jannis Schuecker
-
"""
import json
@@ -197,7 +191,8 @@ class Theory():
mu_CC = np.array([])
sigma2_CC = np.array([])
if self.network.params['connection_params']['replace_cc'] == 'het_poisson_stat':
- with open(self.network.params['connection_params']['replace_cc_input_file'], 'r') as f:
+ with open(self.network.params['connection_params'][
+ 'replace_cc_input_source'], 'r') as f:
rates = json.load(f)
self.cc_input_rates = dict_to_vector(rates,
self.network.area_list,
diff --git a/tests/test_analysis.py b/tests/test_analysis.py
new file mode 100644
index 0000000..fdb3f95
--- /dev/null
+++ b/tests/test_analysis.py
@@ -0,0 +1,33 @@
+import os
+from multiarea_model import MultiAreaModel
+
+"""
+Test analysis class:
+Run a small simulation, load data, compute all measures
+available in Analysis, save data and try load again.
+"""
+
+
+def test_analysis():
+ base_dir = os.getcwd()
+ fn = os.path.join(base_dir, 'fullscale_rates.json')
+ network_params = {'connection_params': {'replace_non_simulated_areas': 'het_poisson_stat',
+ 'replace_cc_input_source': fn},
+ 'N_scaling': 0.001,
+ 'K_scaling': 0.0001,
+ 'fullscale_rates': 'fullscale_rates.json'}
+ sim_params = {'t_sim': 500.,
+ 'areas_simulated': ['V1', 'V2']}
+ M = MultiAreaModel(network_params, simulation=True, sim_spec=sim_params, analysis=True)
+ M.simulation.simulate()
+ M = MultiAreaModel(network_params, simulation=True, sim_spec=sim_params, analysis=True)
+
+ M.analysis.create_pop_rates(t_min=100.)
+ M.analysis.create_pop_rate_dists(t_min=100.)
+ M.analysis.create_synchrony(t_min=100.)
+ M.analysis.create_rate_time_series(t_min=100.)
+ M.analysis.create_synaptic_input(t_min=100.)
+ M.analysis.create_pop_cv_isi(t_min=100.)
+ M.analysis.create_pop_LvR(t_min=100.)
+
+ M.analysis.save()
diff --git a/tests/test_replace_cc.py b/tests/test_replace_cc.py
index 43be22a..2b2791d 100644
--- a/tests/test_replace_cc.py
+++ b/tests/test_replace_cc.py
@@ -1,4 +1,3 @@
-import h5py_wrapper.wrapper as h5w
import json
import numpy as np
import os
@@ -6,6 +5,7 @@ from multiarea_model import MultiAreaModel
from multiarea_model.multiarea_helpers import vector_to_dict
from multiarea_model.multiarea_helpers import create_mask
from multiarea_model.default_params import complete_area_list, population_list
+from multiarea_model.analysis_helpers import _save_dict_to_npy
"""
Test replacing cortico-cortical connections.
@@ -23,7 +23,7 @@ def test_het_poisson_stat_mf():
json.dump(rates, f)
network_params = {'connection_params': {'replace_cc': 'het_poisson_stat',
- 'replace_cc_input_file': 'mf_rates.json'}}
+ 'replace_cc_input_source': 'mf_rates.json'}}
theory_params = {}
M = MultiAreaModel(network_params, theory=True, theory_spec=theory_params)
p, r = M.theory.integrate_siegert()
@@ -53,7 +53,7 @@ def test_het_poisson_stat_sim():
base_dir = os.getcwd()
fn = os.path.join(base_dir, 'fullscale_rates.json')
network_params = {'connection_params': {'replace_cc': 'het_poisson_stat',
- 'replace_cc_input_file': fn},
+ 'replace_cc_input_source': fn},
'N_scaling': 0.001,
'K_scaling': 0.0001,
'fullscale_rates': 'fullscale_rates.json'}
@@ -75,12 +75,12 @@ def test_hom_poisson_stat_sim():
def test_het_current_non_stat_sim():
curr = np.ones(10) * 10.
het_current = {area: {pop: curr for pop in population_list} for area in complete_area_list}
- h5w.save('het_current.h5', het_current, write_mode='w')
+ _save_dict_to_npy('het_current', het_current)
base_dir = os.getcwd()
- fn = os.path.join(base_dir, 'het_current.h5')
+ fs = os.path.join(base_dir, 'het_current')
network_params = {'connection_params': {'replace_cc': 'het_current_nonstat',
- 'replace_cc_input_file': fn},
+ 'replace_cc_input_source': fs},
'N_scaling': 0.001,
'K_scaling': 0.0001,
'fullscale_rates': 'fullscale_rates.json'}
diff --git a/tests/test_simulate_V1.py b/tests/test_simulate_V1.py
index 7754d6b..bf239cb 100644
--- a/tests/test_simulate_V1.py
+++ b/tests/test_simulate_V1.py
@@ -1,8 +1,8 @@
-import h5py_wrapper.wrapper as h5w
import numpy as np
import os
from multiarea_model import MultiAreaModel
from multiarea_model.default_params import complete_area_list, population_list
+from multiarea_model.analysis_helpers import _save_dict_to_npy
"""
Test simulating only V1
@@ -13,7 +13,7 @@ def test_het_poisson_stat_sim():
base_dir = os.getcwd()
fn = os.path.join(base_dir, 'fullscale_rates.json')
network_params = {'connection_params': {'replace_non_simulated_areas': 'het_poisson_stat',
- 'replace_cc_input_file': fn},
+ 'replace_cc_input_source': fn},
'N_scaling': 0.001,
'K_scaling': 0.0001,
'fullscale_rates': 'fullscale_rates.json'}
@@ -38,12 +38,12 @@ def test_hom_poisson_stat_sim():
def test_het_current_non_stat_sim():
curr = np.ones(10) * 10.
het_current = {area: {pop: curr for pop in population_list} for area in complete_area_list}
- h5w.save('het_current.h5', het_current, write_mode='w')
+ _save_dict_to_npy('het_current', het_current)
base_dir = os.getcwd()
- fn = os.path.join(base_dir, 'het_current.h5')
+ fs = os.path.join(base_dir, 'het_current')
network_params = {'connection_params': {'replace_non_simulated_areas': 'het_current_nonstat',
- 'replace_cc_input_file': fn},
+ 'replace_cc_input_source': fs},
'N_scaling': 0.001,
'K_scaling': 0.0001,
'fullscale_rates': 'fullscale_rates.json'}
--
GitLab