Merge pull request #395 from Hjorthmedh/region_mesh_redux

Region mesh redux

snudda/core.py

@@ -216,6 +216,8 @@ def place_neurons(self,
         self.stop_parallel()
         self.close_log_file()
 
+        return sp
+
     ############################################################################
 
     def touch_detection_wrapper(self, args):
@@ -318,6 +320,8 @@ def touch_detection(self,
         self.stop_parallel()
         self.close_log_file()
 
+        return sd, sp
+
     ############################################################################
 
     def prune_synapses_wrapper(self, args):
@@ -390,6 +394,8 @@ def prune_synapses(self,
         self.stop_parallel()
         self.close_log_file()
 
+        return sp
+
     ############################################################################
 
     def setup_input_wrapper(self, args):
@@ -472,6 +478,8 @@ def setup_input(self,
         self.stop_parallel()
         self.close_log_file()
 
+        return si
+
     ############################################################################
 
     def export_to_SONATA_wrapper(self, args):
@@ -752,6 +760,7 @@ def simulate(self,
             print(f"Program run time: {stop - start:.1f}s")
 
         # sim.plot()
+        return sim
 
     ############################################################################
 

snudda/detect/detect.py

@@ -550,26 +550,33 @@ def get_hypervoxel_coords_and_section_id(self, neuron=None, neuron_info=None):
             hyper_voxel_id = self.hyper_voxel_id_lookup[tuple(hyper_voxel_coords[inside_idx, :].T)]
             section_type_id = subtree.section_data[inside_idx, :][:, [2, 0]]
 
-            hid_st_sid = np.hstack([hyper_voxel_id.reshape([hyper_voxel_id.shape[0], 1]), section_type_id])
-
-            # We also need to add parent points with the child branch's section id
-            # This is so we do not miss the first bit between the parent point and the first real point of the branch
-            parent_rows = []
-            if inside_idx.all():
-                # If inside_idx are all True then hyper_voxel_id is same length
-                # as subtree.geometry (should be valid for all but possibly Virtual Axons)
-
-                for tree_type in subtree.sections:
-                    for section in subtree.sections[tree_type].values():
-
-                        if section.section_type == 1:
-                            # Soma has no parent, skip
-                            continue
-
-                        parent_idx = section.point_idx[0]
-                        parent_rows.append([hyper_voxel_id[parent_idx], section.section_type, section.section_id])
-
-            tree_info[subtree_name] = np.unique(np.vstack([hid_st_sid, parent_rows]), axis=0)
+            tree_info[subtree_name] = np.hstack([hyper_voxel_id.reshape([hyper_voxel_id.shape[0], 1]), section_type_id])
+
+            # TODO: This should no longer be necessary! PARENT POINT should be included if
+            #       parent section type is the same
+            #
+            # # We also need to add parent points with the child branch's section id
+            # # This is to not miss the first bit between the parent point and the first real point of the branch
+            # parent_rows = []
+            # if inside_idx.all():
+            #     # If inside_idx are all True then hyper_voxel_id is same length
+            #     # as subtree.geometry (should be valid for all but possibly Virtual Axons)
+            #
+            #     for tree_type in subtree.sections:
+            #         for section in subtree.sections[tree_type].values():
+            #
+            #             if section.section_type == 1:
+            #                 # Soma has no parent, skip
+            #                 continue
+            #
+            #             parent_idx = section.point_idx[0]
+            #             parent_rows.append([hyper_voxel_id[parent_idx], section.section_type, section.section_id])
+            #
+            # if len(parent_rows) > 0:
+            #     tree_info[subtree_name] = np.unique(np.vstack([hid_st_sid, parent_rows]), axis=0)
+            # else:
+            #     # This is the case if we only have a soma
+            #     tree_info[subtree_name] = hid_st_sid
 
             # OBS, there is a rare case when a line segment starts in a hyper voxel, crosses a second hyper voxel
             # and ends up in a third hyper voxel. In this case the intermediate second hyper voxel will be missed if

snudda/neurons/morphology_data.py

@@ -227,6 +227,12 @@ def section_iterator_selective(self, section_type, section_id):
         for sid in section_id:
             yield self.sections[section_type][sid]
 
+    def has_axon(self):
+        return len(self.sections[2]) > 0
+
+    def has_dendrite(self):
+        return len(self.sections[3]) > 0
+
     def load_swc_file(self, swc_file=None, remapping_types={4: 3}, use_cache=True):
 
         """ Loads SWC morphology, not SNUDDA_DATA aware (file must exist).
@@ -252,6 +258,10 @@ def load_swc_file(self, swc_file=None, remapping_types={4: 3}, use_cache=True):
 
         data = np.loadtxt(swc_file)
 
+        if len(data.shape) == 1:
+            # This is to handle case when we only have a soma
+            data = data.reshape([1, 7])
+
         if any(np.diff(data[:, 0]) != 1):
             raise IndexError(f"SWC file has gaps in ID numbering ({swc_file})")
 

snudda/neurons/neuron_model_extended.py

@@ -277,8 +277,6 @@ def define_morphology(self, replace_axon=True, morph_file=None):
 
         assert (morph_file is not None)
 
-        # print("Using morphology: " + morph_file)
-
         return ephys.morphologies.NrnFileMorphology(morph_file, do_replace_axon=replace_axon)
 
     # OLD BUGFIX FOR segment pop
@@ -301,10 +299,6 @@ def build_section_lookup(self):
         # Soma is -1
         self.section_lookup[-1] = self.icell.soma[0]
 
-        # Dendrites are consecutive numbers starting from 1
-        # Ie neurons dend(0) is in pos 1, dend(99) is in pos 100
-        # This so we don't need to special treat soma (pos 0)
-
         for ic, c in enumerate(self.icell.dend):
             self.section_lookup[ic] = c
 

snudda/simulate/simulate.py

@@ -1234,6 +1234,9 @@ def set_resting_voltage(self, neuron_id, rest_volt=None):
 
             self.write_log(f"Neuron {self.neurons[neuron_id].name} ({neuron_id}) resting voltage = {rest_volt * 1e3}")
 
+            # import pdb
+            # pdb.set_trace()
+
             soma = [x for x in self.neurons[neuron_id].icell.soma]
             axon = [x for x in self.neurons[neuron_id].icell.axon]
             dend = [x for x in self.neurons[neuron_id].icell.dend]
@@ -1387,7 +1390,7 @@ def add_volt_recording_soma(self, cell_id=None):
         for cid in cell_id:
 
             if cid in self.neurons.keys() and not self.is_virtual_neuron[cid]:
-                self.add_volt_recording(cid, [0], [0.5])
+                self.add_volt_recording(cid, [-1], [0.5])  # Soma is sec_id=-1
 
     def add_volt_recording(self, cell_id: int, sec_id, sec_x):