diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json new file mode 100644 index 0000000..033cb5e --- /dev/null +++ b/.devcontainer/devcontainer.json @@ -0,0 +1,25 @@ +// For format details, see https://aka.ms/devcontainer.json. For config options, see the +// README at: https://github.com/devcontainers/templates/tree/main/src/ubuntu +{ + "name": "Ubuntu", + // Or use a Dockerfile or Docker Compose file. More info: https://containers.dev/guide/dockerfile + "image": "mcr.microsoft.com/devcontainers/base:jammy", + "features": { + "ghcr.io/devcontainers/features/python:1": {} + } + + // Features to add to the dev container. More info: https://containers.dev/features. + // "features": {}, + + // Use 'forwardPorts' to make a list of ports inside the container available locally. + // "forwardPorts": [], + + // Use 'postCreateCommand' to run commands after the container is created. + // "postCreateCommand": "uname -a", + + // Configure tool-specific properties. + // "customizations": {}, + + // Uncomment to connect as root instead. More info: https://aka.ms/dev-containers-non-root. + // "remoteUser": "root" +} diff --git a/.gitignore b/.gitignore index c590f49..56d281d 100644 --- a/.gitignore +++ b/.gitignore @@ -3,12 +3,6 @@ __pycache__/ *.py[cod] *$py.class - - - - - - # Distribution / packaging .Python build/ @@ -105,3 +99,6 @@ venv.bak/ # mypy .mypy_cache/ + +# macos +.DS_Store diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..6ba1afd --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,6 @@ +{ + "[python]": { + "editor.defaultFormatter": "ms-python.black-formatter" + }, + "python.formatting.provider": "none" +} \ No newline at end of file diff --git a/ouropy/genconnection.py b/ouropy/genconnection.py index 0c5a433..deea348 100644 --- a/ouropy/genconnection.py +++ b/ouropy/genconnection.py @@ -23,15 +23,15 @@ def __init__(self): def get_description(self): """Return a descriptive string for the connection""" - name = self.pre_pop.name + ' to ' + self.post_pop.name + '\n' + name = f'{self.pre_pop.name} to {self.post_pop.name}' + '\n' pre_cell_targets = '\n'.join([str(x) for x in self.pre_cell_targets]) return name + pre_cell_targets def get_name(self): if type(self.pre_pop) == str: - return self.pre_pop + ' to ' + str(self.post_pop) + return f'{self.pre_pop} to {str(self.post_pop)}' else: - return str(self.pre_pop) + ' to ' + str(self.post_pop) + return f'{str(self.pre_pop)} to {str(self.post_pop)}' def get_properties(self): """Get the and make them suitable for pickling""" @@ -128,12 +128,12 @@ def __init__(self, pre_pop, post_pop, for idx, curr_cell_pos in enumerate(pre_pop_pos): - curr_dist = [] - for post_cell_pos in post_pop_pos: - curr_dist.append(euclidian_dist(curr_cell_pos, post_cell_pos)) - + curr_dist = [ + euclidian_dist(curr_cell_pos, post_cell_pos) + for post_cell_pos in post_pop_pos + ] sort_idc = np.argsort(curr_dist) - closest_cells = sort_idc[0:target_pool] + closest_cells = sort_idc[:target_pool] picked_cells = np.random.choice(closest_cells, divergence, replace=False) @@ -248,20 +248,18 @@ def __init__(self, pre_pop, post_pop, for idx, curr_cell_pos in enumerate(pre_pop_pos): - curr_dist = [] - for post_cell_pos in post_pop_pos: - curr_dist.append(euclidian_dist(curr_cell_pos, post_cell_pos)) - + curr_dist = [ + euclidian_dist(curr_cell_pos, post_cell_pos) + for post_cell_pos in post_pop_pos + ] sort_idc = np.argsort(curr_dist) - closest_cells = sort_idc[0:target_pool] + closest_cells = sort_idc[:target_pool] picked_cells = np.random.choice(closest_cells, divergence, replace=False) pre_cell_target.append(picked_cells) for tar_c in picked_cells: - curr_syns = [] - curr_netcons = [] curr_conductances = [] curr_syn = h.pyr2pyr(post_pop[tar_c].soma(0.5)) @@ -293,14 +291,11 @@ def __init__(self, pre_pop, post_pop, curr_syn.Cdur_nmda = Cdur_nmda curr_syn.gbar_nmda = gbar_nmda - curr_syns.append(curr_syn) + curr_syns = [curr_syn] curr_netcon = h.NetCon(pre_pop[idx].soma(0.5)._ref_v, curr_syn, thr, Delay, weight, sec=pre_pop[idx].soma) - #curr_gvec = h.Vector() - #curr_gvec.record(curr_syn._ref_g) - #curr_conductances.append(curr_gvec) - curr_netcons.append(curr_netcon) + curr_netcons = [curr_netcon] netcons.append(curr_netcons) synapses.append(curr_syns) conductances.append(curr_conductances) @@ -1026,7 +1021,7 @@ def write_aps(self, directory='', fname=''): directory = os.getcwd() if not os.path.isdir(directory): os.mkdir(directory) - path = directory + '\\' + fname + '.npz' + path = os.path.join(directory, f'{fname}.npz') try: ap_list = [x[0].as_numpy() for x in self.ap_counters] except: diff --git a/ouropy/gennetwork.py b/ouropy/gennetwork.py index 1f7e0b7..cf0cb35 100644 --- a/ouropy/gennetwork.py +++ b/ouropy/gennetwork.py @@ -77,8 +77,10 @@ def __init__(self, celltypes=None, cellnums=None): if cellnums is None: return - for idx, cell_type in enumerate(celltypes): - self.populations.append(Population(cell_type, cellnums[idx], self)) + self.populations.extend( + Population(cell_type, cellnums[idx], self) + for idx, cell_type in enumerate(celltypes) + ) def mk_population(self, cell_type, n_cells): """Initialize instance empty or with cell populations. @@ -118,8 +120,6 @@ def set_numpy_seed(self, seed): def run_network(self, tstop=1000, dt=1): raise NotImplemented("run_network is not implemented yet") - h.tstop = tstop - h.run() def plot_aps(self, time=200): fig = plt.figure(figsize=(8.27, 11.69)) @@ -138,8 +138,7 @@ def plot_aps(self, time=200): plt.subplot(4, 1, idx+1) plt.eventplot(cells) - plt.ylabel(str(pop) + '\n' + str(pop.perc_active_cells())[0:4] + - '% active') + plt.ylabel((str(pop) + '\n' + str(pop.perc_active_cells())[:4] + '% active')) plt.xlim((0, time)) plt.xlabel("time (ms)") return fig @@ -149,24 +148,24 @@ def save_ap_fig(self, fig, directory=None, file_name=None): directory = os.getcwd() if not file_name: loc_time_str = '_'.join(time.asctime(time.localtime()).split(' ')) - file_name = str(self) + '_' + loc_time_str + file_name = f'{str(self)}_{loc_time_str}' if not os.path.isdir(directory): os.mkdir(directory) - full_file_path = directory + "\\" + file_name + full_file_path = os.path.join(directory, file_name) if os.path.isfile(full_file_path): raise ValueError("The file already exists.\n" + "shelve_network does not overwrite files.") - fig.savefig(full_file_path + ".pdf", dpi=300, format='pdf') - fig.savefig(full_file_path + ".eps", dpi=300, format='eps') + fig.savefig(f"{full_file_path}.pdf", dpi=300, format='pdf') + fig.savefig(f"{full_file_path}.eps", dpi=300, format='eps') plt.close() def get_properties(self): - properties = {'populations': [x.get_properties() - for x in self.populations], - 'init_params': self.init_params} - return properties + return { + 'populations': [x.get_properties() for x in self.populations], + 'init_params': self.init_params, + } def shelve_network(self, directory=None, file_name=None): """Saves the complete network information to a python shelve file. @@ -177,11 +176,11 @@ def shelve_network(self, directory=None, file_name=None): directory = os.getcwd() if not file_name: loc_time_str = '_'.join(time.asctime(time.localtime()).split(' ')) - file_name = str(self) + '_' + loc_time_str + file_name = f'{str(self)}_{loc_time_str}' if not os.path.isdir(directory): os.mkdir(directory) - full_file_path = directory + "\\" + file_name + '.pydd' + full_file_path = os.path.join(directory, f'{file_name}.pydd') if os.path.isfile(full_file_path): raise ValueError("The file already exists.\n" + "shelve_network does not overwrite files.") @@ -190,6 +189,26 @@ def shelve_network(self, directory=None, file_name=None): # BUG with paradigm_frequency_inhibition at 1Hz, possibly too long sim? curr_shelve[str(self)] = self.get_properties() curr_shelve.close() + + def shelve_aps(self, directory=None, file_name=None): + if not directory: + directory = os.getcwd() + if not file_name: + loc_time_str = '_'.join(time.asctime(time.localtime()).split(' ')) + file_name = f'{str(self)}_{loc_time_str}' + if not os.path.isdir(directory): + os.mkdir(directory) + + full_file_path = os.path.join(directory, f'{file_name}.pydd') + if os.path.isfile(full_file_path): + raise ValueError("The file already exists.\n" + + "shelve_aps does not overwrite files.") + + curr_shelve = shelve.open(full_file_path, flag='n') + # BUG with paradigm_frequency_inhibition at 1Hz, possibly too long sim? + curr_shelve['populations'] = [[(i, timestamps) for i, timestamps in enumerate(p.get_timestamps()) if len(timestamps) > 0] for p in self.populations] + curr_shelve.close() + def __str__(self): return str(self.__class__).split("'")[1] @@ -204,15 +223,15 @@ def __init__(self): def get_description(self): """Return a descriptive string for the connection""" - name = self.pre_pop.name + ' to ' + self.post_pop.name + '\n' + name = f'{self.pre_pop.name} to {self.post_pop.name}' + '\n' pre_cell_targets = '\n'.join([str(x) for x in self.pre_cell_targets]) return name + pre_cell_targets def get_name(self): if type(self.pre_pop) == str: - return self.pre_pop + ' to ' + str(self.post_pop) + return f'{self.pre_pop} to {str(self.post_pop)}' else: - return str(self.pre_pop) + ' to ' + str(self.post_pop) + return f'{str(self.pre_pop)} to {str(self.post_pop)}' def get_properties(self): """Get the and make them suitable for pickling""" @@ -310,12 +329,12 @@ def __init__(self, pre_pop, post_pop, for idx, curr_cell_pos in enumerate(pre_pop_pos): - curr_dist = [] - for post_cell_pos in post_pop_pos: - curr_dist.append(euclidian_dist(curr_cell_pos, post_cell_pos)) - + curr_dist = [ + euclidian_dist(curr_cell_pos, post_cell_pos) + for post_cell_pos in post_pop_pos + ] sort_idc = np.argsort(curr_dist) - closest_cells = sort_idc[0:target_pool] + closest_cells = sort_idc[:target_pool] picked_cells = np.random.choice(closest_cells, divergence, replace=False) @@ -439,10 +458,10 @@ def __init__(self, pre_pop, post_pop, for idx, curr_cell_pos in enumerate(pre_pop_pos): - curr_dist = [] - for post_cell_pos in post_pop_pos: - curr_dist.append(euclidian_dist(curr_cell_pos, post_cell_pos)) - + curr_dist = [ + euclidian_dist(curr_cell_pos, post_cell_pos) + for post_cell_pos in post_pop_pos + ] sort_idc = np.argsort(curr_dist) picked_cells = np.random.choice(sort_idc, divergence, replace=True, p=pdf) @@ -553,12 +572,12 @@ def __init__(self, pre_pop, post_pop, netcons = [] for idx, curr_cell_pos in enumerate(pre_pop_pos): - curr_dist = [] - for post_cell_pos in post_pop_pos: - curr_dist.append(euclidian_dist(curr_cell_pos, post_cell_pos)) - + curr_dist = [ + euclidian_dist(curr_cell_pos, post_cell_pos) + for post_cell_pos in post_pop_pos + ] sort_idc = np.argsort(curr_dist) - closest_cells = sort_idc[0:target_pool] + closest_cells = sort_idc[:target_pool] picked_cells = np.random.choice(closest_cells, divergence, replace=False) @@ -622,12 +641,12 @@ def __init__(self, pre_pop, post_pop, target_pool, target_segs, divergence, for idx, curr_cell_pos in enumerate(pre_pop_pos): - curr_dist = [] - for post_cell_pos in post_pop_pos: - curr_dist.append(euclidian_dist(curr_cell_pos, post_cell_pos)) - + curr_dist = [ + euclidian_dist(curr_cell_pos, post_cell_pos) + for post_cell_pos in post_pop_pos + ] sort_idc = np.argsort(curr_dist) - closest_cells = sort_idc[0:target_pool] + closest_cells = sort_idc[:target_pool] picked_cells = np.random.choice(closest_cells, divergence, replace=False) @@ -885,7 +904,7 @@ def make_cells(self, cell_type, n_cells): if not hasattr(self, 'cells'): self.cells = [] - for x in range(n_cells): + for _ in range(n_cells): self.cells.append(cell_type()) self.cells = np.array(self.cells, dtype=object) @@ -903,10 +922,7 @@ def get_cell_number(self): return len(self.cells) def record_aps(self): - counters = [] - for cell in self.cells: - counters.append(cell._AP_counter()) - + counters = [cell._AP_counter() for cell in self.cells] self.ap_counters = counters return counters @@ -933,13 +949,13 @@ def write_aps(self, directory='', fname=''): time_str = '_'.join(time_str.split(' ')) nw_name = self.parent_network.__class__.name pop_name = self.cell_type.name - fname = nw_name + '_' + pop_name + '_' + time_str + fname = f'{nw_name}_{pop_name}_{time_str}' fname = fname.replace(':', '-') if not directory: directory = os.getcwd() if not os.path.isdir(directory): os.mkdir(directory) - path = directory + '\\' + fname + '.npz' + path = os.path.join(directory, f'{fname}.npz') try: ap_list = [x[0].as_numpy() for x in self.ap_counters] except: @@ -972,8 +988,10 @@ def mk_current_clamp(self, cells, amp=0.3, dur=5, delays=3): self.cells[cell]._current_clamp_soma(amp=amp, dur=dur, delay=delay) def get_timestamps(self): - ap_list = [np.array(x[0]) for x in self.ap_counters] - return ap_list + try: + return [x[0].as_numpy() for x in self.ap_counters] + except: + return [np.array(x[0]) for x in self.ap_counters] def current_clamp_rnd(self, n_cells, amp=0.3, dur=5, delay=3): """DEPRECATE""" @@ -1003,10 +1021,7 @@ def add_connection(self, conn): def get_properties(self): """Get the properties of the network""" - try: - ap_time_stamps = [x[0].as_numpy() for x in self.ap_counters] - except: - ap_time_stamps = [np.array(x[0]) for x in self.ap_counters] + ap_time_stamps = self.get_timestamps() ap_numbers = [x[1].n for x in self.ap_counters] try: v_rec = [x.as_numpy() for x in self.VRecords] @@ -1014,21 +1029,19 @@ def get_properties(self): except: v_rec = [np.array(x) for x in self.VRecords] vclamp_i = [np.array(x) for x in self.VClamps_i] - properties = {'parent_network': str(self.parent_network), - 'cell_type': self.cell_type.name, - 'cell_number': self.get_cell_number(), - 'connections': [conn.get_properties() - for conn in self.connections], - 'ap_time_stamps': ap_time_stamps, - 'ap_number': ap_numbers, - 'v_records': v_rec, - 'VClamps_i': vclamp_i} - - properties - return properties + return { + 'parent_network': str(self.parent_network), + 'cell_type': self.cell_type.name, + 'cell_number': self.get_cell_number(), + 'connections': [conn.get_properties() for conn in self.connections], + 'ap_time_stamps': ap_time_stamps, + 'ap_number': ap_numbers, + 'v_records': v_rec, + 'VClamps_i': vclamp_i, + } def __str__(self): - return self.cell_type.name + 'Population' + return f'{self.cell_type.name}Population' def __iter__(self): return self diff --git a/ouropy/genpopulation.py b/ouropy/genpopulation.py index b280fcf..d08eeb0 100644 --- a/ouropy/genpopulation.py +++ b/ouropy/genpopulation.py @@ -112,7 +112,7 @@ def make_cells(self, cell_type, n_cells): if not hasattr(self, 'cells'): self.cells = [] - for x in range(n_cells): + for _ in range(n_cells): self.cells.append(cell_type()) self.cells = np.array(self.cells, dtype=object) @@ -122,10 +122,7 @@ def get_cell_number(self): return len(self.cells) def record_aps(self): - counters = [] - for cell in self.cells: - counters.append(cell._AP_counter()) - + counters = [cell._AP_counter() for cell in self.cells] self.ap_counters = counters return counters @@ -152,13 +149,13 @@ def write_aps(self, directory='', fname=''): time_str = '_'.join(time_str.split(' ')) nw_name = self.parent_network.__class__.name pop_name = self.cell_type.name - fname = nw_name + '_' + pop_name + '_' + time_str + fname = f'{nw_name}_{pop_name}_{time_str}' fname = fname.replace(':', '-') if not directory: directory = os.getcwd() if not os.path.isdir(directory): os.mkdir(directory) - path = directory + '\\' + fname + '.npz' + path = os.path.join(directory, f'{fname}.npz') try: ap_list = [x[0].as_numpy() for x in self.ap_counters] except: @@ -247,7 +244,7 @@ def get_properties(self): return properties def __str__(self): - return self.cell_type.name + 'Population' + return f'{self.cell_type.name}Population' def __iter__(self): return self diff --git a/paradigm_pattern_separation_baseline.py b/paradigm_pattern_separation_baseline.py index b77d9c6..48c3318 100644 --- a/paradigm_pattern_separation_baseline.py +++ b/paradigm_pattern_separation_baseline.py @@ -5,6 +5,7 @@ @author: DanielM """ + from neuron import h, gui # gui necessary for some parameters to h namespace import numpy as np from pydentate import net_tunedrev, neuron_tools @@ -65,7 +66,7 @@ dll_dir = os.path.join(dirname, 'x86_64', 'libnrnmech.so') print("DLL loaded from: " + dll_dir) # h.nrn_load_dll(dll_dir) -""" +""" neuron_tools.load_compiled_mechanisms(path='precompiled') # Start the runs of the model @@ -85,28 +86,28 @@ PP_to_GCs = [] for x in start_idc: - curr_idc = np.concatenate((GC_indices[x:2000], GC_indices[0:x])) + curr_idc = np.concatenate((GC_indices[x:2000], GC_indices[:x])) PP_to_GCs.append(np.random.choice(curr_idc, size=100, replace=False, p=pdf_gc)) PP_to_GCs = np.array(PP_to_GCs) - PP_to_GCs = PP_to_GCs[0:24] + PP_to_GCs = PP_to_GCs[:24] BC_indices = np.arange(24) start_idc = np.array(((start_idc/2000.0)*24), dtype=int) PP_to_BCs = [] for x in start_idc: - curr_idc = np.concatenate((BC_indices[x:24], BC_indices[0:x])) + curr_idc = np.concatenate((BC_indices[x:24], BC_indices[:x])) PP_to_BCs.append(np.random.choice(curr_idc, size=1, replace=False, p=pdf_bc)) PP_to_BCs = np.array(PP_to_BCs) - PP_to_BCs = PP_to_BCs[0:24] + PP_to_BCs = PP_to_BCs[:24] # Generate temporal patterns for the 100 PP inputs temporal_patterns = inhom_poiss(modulation_rate=input_frequency) - temporal_patterns[0:24] + temporal_patterns[:24] nw = net_tunedrev.TunedNetwork(nw_seed[0], temporal_patterns, PP_to_GCs, PP_to_BCs) @@ -121,21 +122,11 @@ print("Running model") neuron_tools.run_neuron_simulator() - tuned_save_file_name = (str(nw) + "-data-paradigm-local-pattern" + - "-separation_nw-seed_input-seed_input-frequency_scale_run_" + - str(nw_seed[0]) + '_' + - str(input_seed[0]) + '_' + - str(input_frequency[0]) + '_' + - str(input_scale).zfill(3) + '_' + - str(run).zfill(3) + '_') + tuned_save_file_name = f"{str(nw)}-data-paradigm-local-pattern-separation_nw-seed_input-seed_input-frequency_scale_run_{str(nw_seed[0])}_{str(input_seed[0])}_{str(input_frequency[0])}_{str(input_scale).zfill(3)}_{str(run).zfill(3)}_" - nw.shelve_network(savedir, tuned_save_file_name) + nw.shelve_aps(savedir, tuned_save_file_name) fig = nw.plot_aps(time=600) - tuned_fig_file_name = (str(nw) + "_spike-plot_paradigm_local-pattern" + - "-separation_run_scale_seed_input-seed_nw-seed_" + - str(run).zfill(3) + '_' + - str(input_scale).zfill(3) + '_' + str(10000) + - str(input_seed) + str(nw_seed)) + tuned_fig_file_name = f"{str(nw)}_spike-plot_paradigm_local-pattern-separation_run_scale_seed_input-seed_nw-seed_{str(run).zfill(3)}_{str(input_scale).zfill(3)}_10000{str(input_seed)}{str(nw_seed)}" nw.save_ap_fig(fig, savedir, tuned_fig_file_name) diff --git a/pydentate/inputs.py b/pydentate/inputs.py index 3da87de..ecb5db2 100644 --- a/pydentate/inputs.py +++ b/pydentate/inputs.py @@ -5,14 +5,15 @@ @author: Daniel & barisckuru """ +import pdb + import numpy as np +import quantities as pq from elephant import spike_train_generation as stg from neo.core import AnalogSignal -import quantities as pq +from scipy import stats from scipy.stats import skewnorm from skimage.measure import profile_line -from scipy import stats -import pdb def inhom_poiss(modulation_rate=10, max_rate=100, n_cells=400): @@ -26,26 +27,21 @@ def inhom_poiss(modulation_rate=10, max_rate=100, n_cells=400): t = np.arange(0, 0.5, sampling_interval.magnitude) - rate_profile = (np.sin(t*modulation_rate*np.pi*2-np.pi/2) + 1) * max_rate / 2 + rate_profile = (np.sin(t * modulation_rate * np.pi * 2 - np.pi / 2) + 1) * max_rate / 2 - rate_profile_as_asig = AnalogSignal(rate_profile, - units=1*pq.Hz, - t_start=0*pq.s, - t_stop=0.5*pq.s, - sampling_period=sampling_interval) + rate_profile_as_asig = AnalogSignal(rate_profile, units=1 * pq.Hz, t_start=0 * pq.s, t_stop=0.5 * pq.s, sampling_period=sampling_interval) spike_trains = [] - for x in range(n_cells): + for _ in range(n_cells): curr_train = stg.inhomogeneous_poisson_process(rate_profile_as_asig) # We have to make sure that there is sufficient space between spikes. # If there is not, we move the next spike by 0.1ms spike_trains.append(curr_train) - array_like = np.array([np.around(np.array(x.times)*1000, decimals=1) - for x in spike_trains], dtype=np.object) + array_like = np.array([np.around(np.array(x.times) * 1000, decimals=1) for x in spike_trains], dtype=object) for arr_idx in range(array_like.shape[0]): bad_idc = np.argwhere(np.diff(array_like[arr_idx]) == 0).flatten() - bad_idc = bad_idc+1 + bad_idc = bad_idc + 1 while bad_idc.any(): for bad_idx in bad_idc: array_like[arr_idx][bad_idx] = array_like[arr_idx][bad_idx] + 0.1 @@ -54,11 +50,20 @@ def inhom_poiss(modulation_rate=10, max_rate=100, n_cells=400): return array_like + def gaussian_connectivity_gc_bc(n_pre, n_gc, n_bc, n_syn_gc, n_syn_bc, scale_gc, scale_bc): """TODO""" pass -def gaussian_connectivity(n_pre, n_post, n_syn=[100,], scale=[1000, 12]): + +def gaussian_connectivity( + n_pre, + n_post, + n_syn=[ + 100, + ], + scale=[1000, 12], +): """TODO THIS IS A STUB FOR A GENERALIZED VERSION OF gaussian_connectivity_gc_bc Choose n_syn postsynaptic cells for each presynaptic cells. Clean up this function. It is not Pythonic. @@ -84,77 +89,79 @@ def gaussian_connectivity(n_pre, n_post, n_syn=[100,], scale=[1000, 12]): center = np.array(n_post) // 2 gauss = stats.norm(loc=center[0], scale=scale[0]) pdf = gauss.pdf(np.arange(n_post[0])) - pdf = pdf/pdf.sum() + pdf = pdf / pdf.sum() post_idc = np.arange(n_post[0]) - start_idc = np.random.randint(0, n_post[0]-1, size=n_pre) + start_idc = np.random.randint(0, n_post[0] - 1, size=n_pre) out_list = [] for idx, n_post_pop in enumerate(n_post): pre_to_post = [] # pdb.set_trace() for x in start_idc: - curr_idc = np.concatenate((post_idc[x:n_post_pop], post_idc[0:x])) + curr_idc = np.concatenate((post_idc[x:n_post_pop], post_idc[:x])) # pdb.set_trace() - pre_to_post.append(np.random.choice(curr_idc, size=n_syn[idx], replace=False, - p=pdf)) + pre_to_post.append(np.random.choice(curr_idc, size=n_syn[idx], replace=False, p=pdf)) out_list.append(pre_to_post) - if idx+1 < len(n_post): - gauss = stats.norm(loc=center[idx+1], scale=scale[idx+1]) - pdf = gauss.pdf(n_post[idx+1]) - pdf = pdf/pdf.sum() - post_idc = np.arange(n_post[idx+1]) - start_idc = np.array(((start_idc/n_post_pop)*n_post[idx+1]), dtype=int) - + if idx + 1 < len(n_post): + gauss = stats.norm(loc=center[idx + 1], scale=scale[idx + 1]) + pdf = gauss.pdf(n_post[idx + 1]) + pdf = pdf / pdf.sum() + post_idc = np.arange(n_post[idx + 1]) + start_idc = np.array(((start_idc / n_post_pop) * n_post[idx + 1]), dtype=int) return np.array(pre_to_post) -#Solstad 2006 Grid Model + +# Solstad 2006 Grid Model def _grid_maker(spacing, orientation, pos_peak, arr_size, sizexy, max_rate): - #define the params from input here, scale the resulting array for maxrate and sperate the xy for size and shift - arr_size = arr_size #200*200 dp was good enough in terms of resolution + # define the params from input here, scale the resulting array for maxrate and sperate the xy for size and shift + arr_size = arr_size # 200*200 dp was good enough in terms of resolution x, y = pos_peak - pos_peak = np.array([x,y]) + pos_peak = np.array([x, y]) max_rate = max_rate - lambda_spacing = spacing*(arr_size/100) #100 required for conversion - k = (4*np.pi)/(lambda_spacing*np.sqrt(3)) + lambda_spacing = spacing * (arr_size / 100) # 100 required for conversion + k = (4 * np.pi) / (lambda_spacing * np.sqrt(3)) degrees = orientation - theta = np.pi*(degrees/180) + theta = np.pi * (degrees / 180) meterx, metery = sizexy - arrx = meterx*arr_size # *arr_size for defining the 2d array size - arry = metery*arr_size - dims = np.array([arrx,arry]) + arrx = meterx * arr_size # *arr_size for defining the 2d array size + arry = metery * arr_size + dims = np.array([arrx, arry]) rate = np.ones(dims) - #implementation of grid function + # implementation of grid function # 3 k values for 3 cos gratings with different angles to generate grid fields - k1 = ((k/np.sqrt(2))*np.array((np.cos(theta+(np.pi)/12) + np.sin(theta+(np.pi)/12), - np.cos(theta+(np.pi)/12) - np.sin(theta+(np.pi)/12)))).reshape(2,) - k2 = ((k/np.sqrt(2))*np.array((np.cos(theta+(5*np.pi)/12) + np.sin(theta+(5*np.pi)/12), - np.cos(theta+(5*np.pi)/12) - np.sin(theta+(5*np.pi)/12)))).reshape(2,) - k3 = ((k/np.sqrt(2))*np.array((np.cos(theta+(9*np.pi)/12) + np.sin(theta+(9*np.pi)/12), - np.cos(theta+(9*np.pi)/12) - np.sin(theta+(9*np.pi)/12)))).reshape(2,) - - rate[i,j] = (np.cos(np.dot(k1, curr_dist))+ - np.cos(np.dot(k2, curr_dist))+ np.cos(np.dot(k3, curr_dist)))/3 - rate = max_rate*2/3*(rate+1/2) # arr is the resulting 2d grid out of 3 gratings + k1 = ((k / np.sqrt(2)) * np.array((np.cos(theta + (np.pi) / 12) + np.sin(theta + (np.pi) / 12), np.cos(theta + (np.pi) / 12) - np.sin(theta + (np.pi) / 12)))).reshape( + 2, + ) + k2 = ((k / np.sqrt(2)) * np.array((np.cos(theta + (5 * np.pi) / 12) + np.sin(theta + (5 * np.pi) / 12), np.cos(theta + (5 * np.pi) / 12) - np.sin(theta + (5 * np.pi) / 12)))).reshape( + 2, + ) + k3 = ((k / np.sqrt(2)) * np.array((np.cos(theta + (9 * np.pi) / 12) + np.sin(theta + (9 * np.pi) / 12), np.cos(theta + (9 * np.pi) / 12) - np.sin(theta + (9 * np.pi) / 12)))).reshape( + 2, + ) + + rate[i, j] = (np.cos(np.dot(k1, curr_dist)) + np.cos(np.dot(k2, curr_dist)) + np.cos(np.dot(k3, curr_dist))) / 3 + rate = max_rate * 2 / 3 * (rate + 1 / 2) # arr is the resulting 2d grid out of 3 gratings return rate - -def _grid_population(n_grid, max_rate, seed, arena_size=[1,1], arr_size=200): + + +def _grid_population(n_grid, max_rate, seed, arena_size=[1, 1], arr_size=200): # skewed normal distribution for grid spacings np.random.seed(seed) median_spc = 43 spc_max = 100 - skewness = 6 #Negative values are left skewed, positive values are right skewed. - grid_spc = skewnorm.rvs(a = skewness,loc=spc_max, size=n_grid) #Skewnorm function - grid_spc = grid_spc - min(grid_spc) #Shift the set so the minimum value is equal to zero. - grid_spc = grid_spc / max(grid_spc) #Standadize all the vlues between 0 and 1. - grid_spc = grid_spc * spc_max #Multiply the standardized values by the maximum value. + skewness = 6 # Negative values are left skewed, positive values are right skewed. + grid_spc = skewnorm.rvs(a=skewness, loc=spc_max, size=n_grid) # Skewnorm function + grid_spc = grid_spc - min(grid_spc) # Shift the set so the minimum value is equal to zero. + grid_spc = grid_spc / max(grid_spc) # Standadize all the vlues between 0 and 1. + grid_spc = grid_spc * spc_max # Multiply the standardized values by the maximum value. grid_spc = grid_spc + (median_spc - np.median(grid_spc)) - - grid_ori = np.random.randint(0, high=60, size=[n_grid,1]) #uniform dist for orientation btw 0-60 degrees - grid_phase = np.random.randint(0, high=(arr_size-1), size=[n_grid,2]) #uniform dist grid phase - + + grid_ori = np.random.randint(0, high=60, size=[n_grid, 1]) # uniform dist for orientation btw 0-60 degrees + grid_phase = np.random.randint(0, high=(arr_size - 1), size=[n_grid, 2]) # uniform dist grid phase + # create a 3d array with grids for n_grid - rate_grids = np.zeros((arr_size, arr_size, n_grid))#empty array + rate_grids = np.zeros((arr_size, arr_size, n_grid)) # empty array for i in range(n_grid): x = grid_phase[i][0] y = grid_phase[i][1] @@ -166,17 +173,11 @@ def _grid_population(n_grid, max_rate, seed, arena_size=[1,1], arr_size=200): def _inhom_poiss(arr, dur_s, poiss_seed=0, dt_s=0.025): np.random.seed(poiss_seed) n_cells = arr.shape[0] - spi_arr = np.zeros((n_cells), dtype = np.ndarray) + spi_arr = np.zeros((n_cells), dtype=np.ndarray) for grid_idc in range(n_cells): - np.random.seed(poiss_seed+grid_idc) - rate_profile = arr[grid_idc,:] - asig = AnalogSignal(rate_profile, - units=1*pq.Hz, - t_start=0*pq.s, - t_stop=dur_s*pq.s, - sampling_period=dt_s*pq.s, - sampling_interval=dt_s*pq.s) + np.random.seed(poiss_seed + grid_idc) + rate_profile = arr[grid_idc, :] + asig = AnalogSignal(rate_profile, units=1 * pq.Hz, t_start=0 * pq.s, t_stop=dur_s * pq.s, sampling_period=dt_s * pq.s, sampling_interval=dt_s * pq.s) curr_train = stg.inhomogeneous_poisson_process(asig) - spi_arr[grid_idc] = np.array(curr_train.times*1000) #time conv to ms + spi_arr[grid_idc] = np.array(curr_train.times * 1000) # time conv to ms return spi_arr -