latest updates

Change-Id: I6049b34ea984fcfa25f6e92b47c099ddad0ce40a

diameter_synthesis/build_diameters.py

@@ -32,7 +32,7 @@
 ##################################
 
 
-def build_diameters(models, models_params, morphologies_dict, neurite_types, new_morph_path, extra_params, morph_path, plot=True, n_cpu=1, n_samples=1):
+def build_diameters(models, models_params, morphologies_dict, neurite_types, new_morph_path, extra_params, morph_path, plot=True, n_cpu=1, n_samples=1, ext = '.png'):
     """ Building the diameters from the generated diameter models"""
 
     all_models = {}
@@ -50,30 +50,27 @@ def build_diameters(models, models_params, morphologies_dict, neurite_types, new
         neurons = []
         for mtype in morphologies_dict:
             for neuron in morphologies_dict[mtype]:
-                name, ext = os.path.splitext(neuron)
-                if ext in {'.h5', '.asc', '.swc'} and os.path.exists(morph_path + '/' + neuron):
+                name, neuron_ext = os.path.splitext(neuron)
+                if neuron_ext in {'.h5', '.asc', '.swc'} and os.path.exists(morph_path + '/' + neuron):
                     neurons.append([neuron, mtype])
 
         # set all parameters
         build_diam_poolf = partial(build_diam_pool, all_models, model, models_params,
-                                   neurite_types, extra_params, morph_path, new_morph_path, plot, n_samples)
+                                   neurite_types, extra_params, morph_path, new_morph_path, plot, n_samples, ext)
 
         # generate diameters in parallel
         with Pool(processes=n_cpu) as p_build:  # initialise the parallel computation
             list(tqdm(p_build.imap(build_diam_poolf, neurons), total=len(neurons)))
 
 
-def build_diam_pool(all_models, model, models_params, neurite_types, extra_params, morph_path, new_morph_path, plot, n_samples, neuron_input):
+def build_diam_pool(all_models, model, models_params, neurite_types, extra_params, morph_path, new_morph_path, plot, n_samples, ext, neuron_input):
     """ build a neuron diameters, save and plot it """
 
     fname = neuron_input[0]
     mtype = neuron_input[1]
-    name, ext = os.path.splitext(fname)
 
     filepath = os.path.join(morph_path, fname)
     neuron = io.load_morphology(filepath)
-    if neuron.name == 'mtC130201A_idB':
-        print(neuron.name)
 
     np.random.seed(extra_params[model]['seed'])
 
@@ -95,8 +92,8 @@ def build_diam_pool(all_models, model, models_params, neurite_types, extra_param
     io.save_neuron(neuron, model, new_morph_path)
     if plot:
         folder = 'shapes_' + os.path.basename(new_morph_path[:-1])
-        plotting.plot_diameter_diff(name, morph_path, new_morph_path,
-                                    model, neurite_types, folder=folder)
+        plotting.plot_diameter_diff(os.path.splitext(fname)[0], morph_path, new_morph_path,
+                                    model, neurite_types, folder=folder, ext = ext)
 
 
 def diametrize_model_generic(neuron, params, neurite_types, extra_params):

diameter_synthesis/build_models.py

@@ -22,7 +22,7 @@
 ##############################################
 
 
-def sampling_model_sibling_asymmetry(morphologies, neurite_types, extra_params, tqdm_disable=False):
+def sampling_model_sibling_asymmetry_trunk(morphologies, neurite_types, extra_params, tqdm_disable=False):
     """ test for sampling models """
 
     sibling_sequential = 'asymmetry_threshold'
@@ -120,6 +120,105 @@ def sampling_model_sibling_asymmetry(morphologies, neurite_types, extra_params,
     return all_models, all_data
 
 
+
+def sampling_model_sibling_asymmetry(morphologies, neurite_types, extra_params, tqdm_disable=False):
+    """ test for sampling models """
+
+    sibling_sequential = 'asymmetry_threshold'
+    rall_deviation_sequential = 'asymmetry_threshold'
+    terminal_diameters_sequential = None
+    trunk_diameters_sequential = None #'max_branch' 
+    tapers_sequential = None
+
+    # initialise dictionaries for collecting morphological quantities
+    sibling_ratios = {}
+    rall_deviations = {}
+    terminal_diameters = {}
+    trunk_diameters = {}
+    tapers = {}
+    for neurite_type in neurite_types:
+        sibling_ratios[neurite_type] = []
+        rall_deviations[neurite_type] = []
+        terminal_diameters[neurite_type] = []
+        trunk_diameters[neurite_type] = []
+        tapers[neurite_type] = []
+
+    # loop first over all morphologies (TODO: could be parallelized)
+    i = 0
+    for neuron in tqdm(morphologies, disable=tqdm_disable):
+        # for each neurite in the neuron
+        for neurite in neuron.neurites:
+            # for each type of neurite we consider
+            for neurite_type in neurite_types:
+
+                if neurite.type == STR_TO_TYPES[neurite_type]:
+
+                    # compute here all the morphological values from the neurite
+                    sibling_ratios[neurite_type] += morph_funcs.sibling_ratios(neurite, seq=sibling_sequential)
+                    rall_deviations[neurite_type] += morph_funcs.rall_deviations(neurite, seq=rall_deviation_sequential)
+                    terminal_diameters[neurite_type] += morph_funcs.terminal_diameters(neurite, threshold=extra_params['terminal_threshold'], seq=terminal_diameters_sequential)
+                    trunk_diameters[neurite_type] += morph_funcs.trunk_diameter(neurite, seq=trunk_diameters_sequential)
+                    tapers[neurite_type] += morph_funcs.taper(neurite, params=extra_params['taper'], seq=tapers_sequential)
+
+    # do the fits of each morphological values
+    sibling_ratio_models = {}
+    rall_deviation_models = {}
+    terminal_diameters_models = {}
+    trunk_diameters_models = {}
+    tapers_models = {}
+    for neurite_type in neurite_types:
+
+        # sibling ratio
+        sibling_ratio_models[neurite_type] = {}
+        sibling_ratio_models[neurite_type]['distribution'] = 'expon_rev'
+        sibling_ratio_models[neurite_type]['sequential'] = sibling_sequential
+        sibling_ratio_models[neurite_type]['params'] = fit_distribution(sibling_ratios[neurite_type], sibling_ratio_models[neurite_type]['distribution'], seq=sibling_sequential, extra_params=extra_params, name ='sibling', threshold=extra_params['threshold'][neurite_type] )
+
+        # Rall deviation
+        rall_deviation_models[neurite_type] = {}
+        rall_deviation_models[neurite_type]['distribution'] = 'exponnorm'
+        rall_deviation_models[neurite_type]['sequential'] = rall_deviation_sequential
+        rall_deviation_models[neurite_type]['params'] = fit_distribution(rall_deviations[neurite_type], rall_deviation_models[neurite_type]['distribution'], seq=rall_deviation_sequential, extra_params=extra_params, name ='Rall', threshold=extra_params['threshold'][neurite_type] )
+
+        # terminal diameters
+        terminal_diameters_models[neurite_type] = {}
+        terminal_diameters_models[neurite_type]['distribution'] = 'exponnorm'
+        terminal_diameters_models[neurite_type]['sequential'] = terminal_diameters_sequential
+        terminal_diameters_models[neurite_type]['params'] = fit_distribution(terminal_diameters[neurite_type], terminal_diameters_models[neurite_type]['distribution'], seq=terminal_diameters_sequential, extra_params=extra_params)
+
+        # trunk diameters
+        trunk_diameters_models[neurite_type] = {}
+        trunk_diameters_models[neurite_type]['distribution'] = 'exponnorm'
+        trunk_diameters_models[neurite_type]['sequential'] = trunk_diameters_sequential
+        trunk_diameters_models[neurite_type]['params'] = fit_distribution(trunk_diameters[neurite_type], trunk_diameters_models[neurite_type]['distribution'], seq=trunk_diameters_sequential, extra_params=extra_params)
+
+        # taper
+        tapers_models[neurite_type] = {}
+        tapers_models[neurite_type]['distribution'] = 'exponnorm'
+        tapers_models[neurite_type]['sequential'] = tapers_sequential
+        tapers_models[neurite_type]['params'] = fit_distribution(
+            tapers[neurite_type], tapers_models[neurite_type]['distribution'], seq=tapers_sequential, extra_params=extra_params)
+
+    # collect all models in one dictionary
+    all_models = {
+        'sibling_ratio': sibling_ratio_models,
+        'rall_deviation': rall_deviation_models,
+        'terminal_diameter': terminal_diameters_models,
+        'trunk_diameter': trunk_diameters_models,
+        'taper': tapers_models
+    }
+
+    all_data = {
+        'sibling_ratio': sibling_ratios,
+        'rall_deviation': rall_deviations,
+        'terminal_diameter': terminal_diameters,
+        'trunk_diameter': trunk_diameters,
+        'taper': tapers
+    }
+
+    return all_models, all_data
+
+
 def sampling_model_trunk_path(morphologies, neurite_types, extra_params, tqdm_disable=False):
     """ test for sampling models """
 
@@ -334,9 +433,9 @@ def build_models(models, morphologies, neurite_types, extra_params, fig_folder='
         if model == 'M0':
             all_models[model] = sampling_model_generic
         if model == 'M1':
-            all_models[model] = sampling_model_trunk_path
-        if model == 'M2':
             all_models[model] = sampling_model_sibling_asymmetry
+        if model == 'M2':
+            all_models[model] = sampling_model_sibling_asymmetry_trunk
 
     tqdm_1, tqdm_2 = utils.tqdm_disable(morphologies)  # to have a single progression bar
 

diameter_synthesis/main.py

@@ -83,5 +83,4 @@ def run_diameters(config_file):
         models_params = json.load(f)
 
     # generate diameters
-    models_params = build_diameters(config['models'], models_params, morphologies_dict, config['neurite_types'], config['new_morph_path'],
-                                    config['extra_params'], config['morph_path'], plot=config['plot'], n_cpu=config['n_cpu'], n_samples=config['n_samples'])
+    models_params = build_diameters(config['models'], models_params, morphologies_dict, config['neurite_types'], config['new_morph_path'], config['extra_params'], config['morph_path'], plot=config['plot'], n_cpu=config['n_cpu'], n_samples=config['n_samples'], ext = config['ext'])

diameter_synthesis/plotting.py

@@ -387,7 +387,7 @@ def draw_axis(obj, mode='2d', ax=None, **kwargs):
     return fig, ax
 
 
-def plot_diameter_diff(neuron_name, morph_path, new_morph_path, model, neurite_types, folder):
+def plot_diameter_diff(neuron_name, morph_path, new_morph_path, model, neurite_types, folder, ext = '.png'):
     """ plot original morphology, new one and differences """
 
     if not os.path.isdir(folder):
@@ -444,7 +444,7 @@ def plot_diameter_diff(neuron_name, morph_path, new_morph_path, model, neurite_t
     draw_axis(neuron_diff_neg, ax=axs[1, 1])
     axs[1, 1].set_title('Negative diameter differences')
 
-    fig.savefig(folder + '/' + neuron_name + '_' + folder + '_' + model + '.png', dpi=500)
+    fig.savefig(folder + '/' + neuron_name + '_' + folder + '_' + model + ext, dpi=500)
     plt.close('all')
 
 

examples/create_jsons.py

@@ -7,7 +7,7 @@
     #'morph_path': '/gpfs/bbp.cscs.ch/project/proj81/InputData/2017Release/OriginalData/05_RepairUnravel-asc/', 
     'mtypes_sort': 'mtypes', 
     'n_morphs_max': None, 
-    'n_mtypes_max': 10,
+    'n_mtypes_max': 60,
     'models': ['M2',], # 'M1', 'M2', 'M3', 'M4', 'M5', 'M6', 'M7', 'M8', 'M9'],
     'neurite_types': ['basal', 'apical'],
     'models_params_file': 'model_params_mtypes.json',

scripts/diameter-checks/mm-config.json

@@ -1,7 +1,7 @@
 {
     "version": "1.0",
     "circuitmvd3_path": "circuit.mvd3",
-    "morph_path": "new_morphologies_super_mtypes",
+    "morph_path": "../../examples/2017_emodels/2017_IN_cells/new_morphologies_all",
     "rep_morph_path": "/gpfs/bbp.cscs.ch/project/proj81/InputData/2017Release/OriginalData/05_RepairUnravel-asc/",
     "unrep_morph_path": "/gpfs/bbp.cscs.ch/project/proj68/home/kanari/DiamCheck/SYN_for_paper/Raw/04_ZeroDiameterFix-ASC/",
     "emodels_repo": "ssh://bbpcode.epfl.ch/project/proj38/singlecell-optimization",

scripts/diameter-checks/run.sh

@@ -10,7 +10,6 @@ module purge all
 
 module load neuron/7.6.8/python2/serial
 
-
 export OMP_NUM_THREADS=1
 
 diameter-check run mm-config.json final.json emodel_etype_map.json

scripts/split_PC_vs_IN/split_PC_VS_IN.py

@@ -11,12 +11,14 @@
 script to split cells in two types: pyramidal cells, and interneuron cells
 """
 
-input_folder = '../../examples/05_RepairUnravel-asc/'
-
-if not os.path.isdir('PC_neurons'):
-    os.mkdir('PC_neurons')
-if not os.path.isdir('Inter_neurons'):
-    os.mkdir('Inter_neurons')
+#input_folder = '../../examples/05_RepairUnravel-asc/'
+input_folder = '../extract_morphologies/selected_morphologies/'
+PC_folder = 'PC_neurons_selected'
+Inter_folder = 'Inter_neurons_selected'
+if not os.path.isdir(PC_folder):
+    os.mkdir(PC_folder)
+if not os.path.isdir(Inter_folder):
+    os.mkdir(Inter_folder)
 
 for fname in os.listdir(input_folder):
 
@@ -25,10 +27,10 @@
         IN = True
         for neurite in neuron.neurites:
             if neurite.type == tpes.STR_TO_TYPES['apical']:
-                shutil.copyfile(input_folder+fname, 'PC_neurons/'+fname)
+                shutil.copyfile(input_folder + fname, PC_folder +'/' + fname)
                 print(fname, '==> PC cell')
                 IN = False
                 break
         if IN:
             print(fname, '==> IN cell')
-            shutil.copyfile(input_folder+fname, 'Inter_neurons/'+fname)
+            shutil.copyfile(input_folder + fname, Inter_folder + '/' + fname)