Density placement: Force scaling on voxels that are not full. Fix #27. (

#28)

Density placement: Force scaling on voxels that are not full. Fix #27.

* Minor fixes in docstrings.
* Fix warnings from divide by zero. Add test to exit density placement.
* update to pandas 2.0
* Fix warnings from nan values.

---------

Co-authored-by: Mike Gevaert <[email protected]>

atlas_densities/app/cell_densities.py

@@ -155,7 +155,7 @@ def cell_density(annotation_path, hierarchy_path, nissl_path, output_path):
 
     \b
     - the Nissl stain intensity, which is supposed to represent the overall cell density, up to
-    to region-dependent constant scaling factors.
+    region-dependent constant scaling factors.
     - cell counts from the scientific literature, which are used to determine a local
     linear dependency factor for each region where a cell count is available.
     - the optional soma radii, used to operate a correction.
@@ -417,7 +417,7 @@ def inhibitory_and_excitatory_neuron_densities(
     Note: optimization is not fully implemented and the current process only returns a
     feasible point.
 
-    The ouput densities are saved in the specified output directory under the following
+    The output densities are saved in the specified output directory under the following
     names:
 
     \b
@@ -490,7 +490,7 @@ def compile_measurements(
     Sexual Dimorphism' by Kim et al., 2017.
     https://ars.els-cdn.com/content/image/1-s2.0-S0092867417310693-mmc3.xlsx
     - `atlas_densities/app/data/measurements/gaba_papers.xlsx`, a compilation of measurements
-    from the scientifc literature made by Rodarie Dimitri (BBP).
+    from the scientific literature made by Rodarie Dimitri (BBP).
 
     This command extracts measurements from the above two files and gathers them into a unique
     CSV file with the following columns:
@@ -531,7 +531,7 @@ def compile_measurements(
     See `atlas_densities/densities/excel_reader.py` for more information.
 
     Note: This function should be deprecated once its output has been stored permanently as the
-    the unique source of density-related measurements for the AIBS mouse brain. New measurements
+    unique source of density-related measurements for the AIBS mouse brain. New measurements
     should be added to the stored file (Nexus).
     """
 
@@ -619,7 +619,7 @@ def measurements_to_average_densities(
     If several combinations of measurements yield several average cell densities for the same
     region, then the output data frame records the average of these measurements.
 
-    The ouput average cell densities are saved in under the CSV format as a dataframe with the same
+    The output average cell densities are saved in under the CSV format as a dataframe with the same
     columns as the input data frame specified via `--measurements-path`.
     See :mod:`atlas_densities.app.densities.compile_measurements`.
 
@@ -744,7 +744,7 @@ def fit_average_densities(
     gad67+) in every homogenous region of type "inhibitory".
 
     Our linear fitting of density values relies on the assumption that the average cell density
-    (number of cells per mm^3) of a cell type T in a brain region R depends linearily on the
+    (number of cells per mm^3) of a cell type T in a brain region R depends linearly on the
     average intensity of a gene marker of T. The conversion factor is a constant which depends only
     on T and on which of the following three groups R belongs to:
 
@@ -894,7 +894,7 @@ def inhibitory_neuron_densities(
     (inhibitory neurons) and those of the inhibitory neuron subtypes PV+, SST+ and VIP+.
 
     The function modifies the estimates in `average_densities` (the "first estimates" of the paper)
-    to make them consistent accross cell types: the average density of GAD67+ cells in a leaf
+    to make them consistent across cell types: the average density of GAD67+ cells in a leaf
     region L should be at most the sum of the average densities of its subtypes under scrutiny
     (e.g. PV+, SST+ and VIP+) and should not exceed the neuron density of L.
 

atlas_densities/densities/cell_density.py

@@ -107,7 +107,10 @@ def compute_cell_density(
     group_ids = get_group_ids(region_map)
     region_masks = get_region_masks(group_ids, annotation)
     fix_purkinje_layer_intensity(region_map, annotation, region_masks, nissl)
+    non_zero_nissl = nissl > 0
     for group, mask in region_masks.items():
-        nissl[mask] = nissl[mask] * (cell_counts()[group] / np.sum(nissl[mask]))
+        mask = np.logical_and(non_zero_nissl, mask)
+        if mask.any():
+            nissl[mask] = nissl[mask] * (cell_counts()[group] / np.sum(nissl[mask]))
 
     return nissl / voxel_volume

atlas_densities/densities/excel_reader.py

@@ -188,7 +188,7 @@ def _set_metadata_columns(
 
     In `dataframe`, the rows corresponding to the same article have no valid 'source title',
     'comment' nor 'specimen age' except the first one. This function makes sure every row is
-    set with approriate values in the aforementioned columns.
+    set with appropriate values in the aforementioned columns.
 
     Args:
         dataframe: DataFrame obtained when reading the worksheets 'GAD67 densities' or 'PV-SST-VIP'
@@ -352,16 +352,14 @@ def read_inhibitory_neuron_measurement_compilation(
     Read the neuron densities of the worksheet 'GAD67 densities' in gaba_papers.xlsx
 
     Args:
-        region_map: RegionMap object to navigate the brain regions hierarchy.
-            Assumed to be instantiated with AIBS 1.json.
         measurements_path: path to the file gaba_papers.xlsx
 
     Returns:
         a pd.DataFrame with the columns listed in the description at the top of
         this module.
     """
 
-    # Reading takes several seconds, possibly due to formulaes interpretation
+    # Reading takes several seconds, possibly due to formulas interpretation
     L.info("Loading excel worksheet ...")
     inhibitory_neurons_worksheet = pd.read_excel(
         str(measurements_path),
@@ -435,8 +433,6 @@ def read_pv_sst_vip_measurement_compilation(measurements_path: Union[str, "Path"
     Read the neuron densities of the worksheet 'PV-SST-VIP' in gaba_papers.xlsx
 
     Args:
-        region_map: RegionMap object to navigate the brain regions hierarchy.
-            Assumed to be instantiated with AIBS 1.json.
         measurements_path: path to the file gaba_papers.xlsx
 
     Returns:
@@ -484,7 +480,7 @@ def read_homogenous_neuron_type_regions(measurements_path: Union[str, "Path"]) -
 
     Returns:
         pd.DataFrame with two columns: 'brain_region' and 'cell_type'. A cell type value
-        is either 'ihibitory' or 'excitatory' and applies to every cell in the region.
+        is either 'inhibitory' or 'excitatory' and applies to every cell in the region.
     """
 
     homogenous_regions_worksheet = pd.read_excel(

atlas_densities/densities/excitatory_inhibitory_splitting.py

@@ -35,7 +35,7 @@ def scale_excitatory_densities(output, brain_regions, mapping, layer_ids, densit
         L.info("Performing layer: %s", layer)
         idx = np.nonzero(np.isin(brain_regions.raw, layer_ids[layer]))
 
-        for mtype, scale in df.iteritems():
+        for mtype, scale in df.items():
             if scale == 0:
                 continue
 

atlas_densities/densities/mtype_densities_from_composition.py

@@ -45,7 +45,7 @@ def create_from_composition(
 
     layer_masks = get_layer_masks(annotation.raw, region_map, metadata)
 
-    for layer, layer_data in excitatory_mtype_composition.groupby(["layer"]):
+    for layer, layer_data in excitatory_mtype_composition.groupby("layer"):
 
         layer_sum = layer_data["density"].sum()
 

atlas_densities/densities/utils.py

@@ -158,7 +158,7 @@ def optimize_distance_to_line(  # pylint: disable=too-many-arguments
     """
     Find inside a box the closest point to a line with prescribed coordinate sum.
 
-    This function aims at solving the following (convex quadratic) optmization problem:
+    This function aims at solving the following (convex quadratic) optimization problem:
 
     Given a sum S >= 0.0, a line D in the non-negative orthant of the Euclidean N-dimensional
     space and a box B in this orthant (an N-dimensional vector with non-negative coordinates),
@@ -176,26 +176,30 @@ def optimize_distance_to_line(  # pylint: disable=too-many-arguments
     Args:
         line_direction_vector: N-dimensional float vector with non-negative coordinates.
         upper_bounds: N-dimensional float vector with non-negative coordinates. Defines the
-            the box constraining the optimization process.
-        sum_constraints: non-negative float number. The coordinate sum constraints imposed on
+            box constraining the optimization process.
+        sum_constraint: non-negative float number. The coordinate sum constraints imposed on
             the point P we are looking for.
         threshold: non-negative float value. If the coordinate sum of the current point
             is below `threshold`, the function returns the current point.
         max_iter: maximum number of iterations.
-        copy: If True, the function makes a copy of the input `line_direction_vector`. Otherwise
+        copy: If True, the function makes a copy of the input `line_direction_vector`. Otherwise,
             `line_direction_vector` is modified in-place and holds the optimal value.
 
     Returns: N-dimensional float vector with non-negative coordinates. The solution point of the
-        optimization problem, if it exists, up to inacurracy due to threshold size or early
-        termination of the algorithm. Otherwise a point on the boundary of the box B defined by
+        optimization problem, if it exists, up to inaccuracy due to threshold size or early
+        termination of the algorithm. Otherwise, a point on the boundary of the box B defined by
         `upper_bounds`.
     """
     diff = float("inf")
     iter_ = 0
     point = line_direction_vector.copy() if copy else line_direction_vector
-    while diff > threshold and iter_ < max_iter:
-        point *= sum_constraint / np.sum(point)
+    scalable_voxels = point != 0
+    while diff > threshold and iter_ < max_iter and scalable_voxels.any():
+        point[scalable_voxels] *= (sum_constraint - np.sum(point[~scalable_voxels])) / np.sum(
+            point[scalable_voxels]
+        )
         point = np.min([point, upper_bounds], axis=0)
+        scalable_voxels = np.logical_and(point != 0, point < upper_bounds)
         diff = np.abs(np.sum(point) - sum_constraint)
         iter_ += 1
 
@@ -385,7 +389,7 @@ def get_group_names(region_map: "RegionMap", cleanup_rest: bool = False) -> Dict
 
     Args:
         region_map: object to navigate the mouse brain regions hierarchy
-            (instantied from AIBS 1.json).
+            (instantiated from AIBS 1.json).
         cleanup_rest: (Optional) If True, the name of any ascendant region of the Cerebellum and
             Isocortex groups are removed from the names of the Rest group. This makes sure that
             the set of names of the Rest group is closed under taking descendants.
@@ -485,7 +489,8 @@ def get_hierarchy_info(
     ]
     region_names = [region_map.get(id_, attr="name") for id_ in region_ids]
     data_frame = pd.DataFrame(
-        {"brain_region": region_names, "descendant_id_set": descendant_id_sets}, index=region_ids
+        {"brain_region": region_names, "descendant_id_set": descendant_id_sets},
+        index=region_ids,
     )
 
     return data_frame