From 7a119f2335ed0d8009637dc4c0f3f6a25d00983b Mon Sep 17 00:00:00 2001 From: Dimitri RODARIE Date: Tue, 19 Mar 2024 11:09:25 +0100 Subject: [PATCH 1/3] Add check to the region filter in case the input dataset is all equal to 0 for the region of interest. (#74) --- atlas_densities/densities/utils.py | 9 ++++----- tests/densities/test_utils.py | 14 ++++++++++++++ 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/atlas_densities/densities/utils.py b/atlas_densities/densities/utils.py index bc6430b..cd917f9 100644 --- a/atlas_densities/densities/utils.py +++ b/atlas_densities/densities/utils.py @@ -68,11 +68,11 @@ def normalize_intensity( 3D array, the normalized marker intensity array. """ - outside_mean = np.mean( - marker_intensity[np.logical_and((annotation == region_id), (marker_intensity > 0.0))] - ) output_intensity = copy_array(marker_intensity, copy=copy) - output_intensity -= outside_mean * threshold_scale_factor + mask_low_filter = np.logical_and((annotation == region_id), (marker_intensity > 0.0)) + if np.any(mask_low_filter): + outside_mean = np.mean(marker_intensity[mask_low_filter]) + output_intensity -= outside_mean * threshold_scale_factor output_intensity[output_intensity < 0.0] = 0.0 max_intensity = np.max(output_intensity) if max_intensity > 0: @@ -405,7 +405,6 @@ def materialize(node): def get_group_ids(region_map: "RegionMap", config: dict, _skip_check=False) -> dict[str, set[int]]: - """ Get AIBS structure ids for several region groups of interest. diff --git a/tests/densities/test_utils.py b/tests/densities/test_utils.py index 51b9ba8..5680a7f 100644 --- a/tests/densities/test_utils.py +++ b/tests/densities/test_utils.py @@ -59,6 +59,20 @@ def test_normalize(): # In place modification actual = tested.normalize_intensity(marker, annotation_raw, copy=False) npt.assert_array_almost_equal(marker, expected) + # Test outside with no intensity + marker = np.array( + [ + [ + [0.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 0.0], + [0.0, 3.0, 5.0, 0.0], + [0.0, 0.0, 0.0, 0.0], + ] + ], + dtype=float, + ) + actual = tested.normalize_intensity(marker, annotation_raw, copy=True) + npt.assert_array_almost_equal(marker / np.max(marker), actual) def test_compensate_cell_overlap(): From c2ad1294985790bf223e79096b08a52d5f348402 Mon Sep 17 00:00:00 2001 From: MikeG Date: Thu, 21 Mar 2024 15:04:05 +0100 Subject: [PATCH 2/3] indent Test_mtype_densities_from_composition test right (#76) --- tests/app/test_mtype_densities.py | 493 +++++++++++++++--------------- 1 file changed, 247 insertions(+), 246 deletions(-) diff --git a/tests/app/test_mtype_densities.py b/tests/app/test_mtype_densities.py index 45610ed..4927a80 100644 --- a/tests/app/test_mtype_densities.py +++ b/tests/app/test_mtype_densities.py @@ -166,273 +166,274 @@ def test_output(self, tmp_path): assert "bAC" in metadata["density_files"]["BP"] assert "EXC" == metadata["synapse_class"] - class Test_mtype_densities_from_composition: - @pytest.fixture(scope="session") - def class_tmpdir(self, tmpdir_factory): - """Create a session scoped temp dir using the class name""" - return tmpdir_factory.mktemp(type(self).__name__) - - @pytest.fixture(scope="session") - def annotation_path(self, class_tmpdir): - path = class_tmpdir.join("annotation.nrrd") - VoxelData( - np.array([[[101, 102, 103, 104, 105, 106]]], dtype=np.int32), - voxel_dimensions=[25] * 3, - ).save_nrrd(str(path)) - return path - - @pytest.fixture(scope="session") - def hierarchy_path(self, class_tmpdir): - path = class_tmpdir.join("hierarchy.json") - - hierarchy = { - "id": 315, - "acronym": "Isocortex", - "name": "Isocortex", - "children": [ - { - "id": 500, - "acronym": "MO", - "name": "Somatomotor areas", - "children": [ - { - "id": 101, - "acronym": "MO1", - "name": "Somatomotor areas, Layer 1", - "children": [], - }, - { - "id": 102, - "acronym": "MO2", - "name": "Somatomotor areas, Layer 2", - "children": [], - }, - { - "id": 103, - "acronym": "MO3", - "name": "Somatomotor areas, Layer 3", - "children": [], - }, - { - "id": 104, - "acronym": "MO4", - "name": "Somatomotor areas, Layer 4", - "children": [], - }, - { - "id": 105, - "acronym": "MO5", - "name": "Somatomotor areas, layer 5", - "children": [], - }, - { - "id": 106, - "acronym": "MO6", - "name": "Somatomotor areas, layer 6", - "children": [], - }, - ], - }, - ], - } - with open(path, "w", encoding="utf-8") as jsonfile: - json.dump(hierarchy, jsonfile, indent=1, separators=(",", ": ")) - return path - - @pytest.fixture(scope="session") - def metadata_path(self, class_tmpdir): - path = class_tmpdir.join("metadata.json") - metadata = { - "region": { - "name": "Isocortex", - "query": "Isocortex", - "attribute": "acronym", - "with_descendants": True, - }, - "layers": { - "names": ["layer_1", "layer_2", "layer_3", "layer_4", "layer_5", "layer_6"], - "queries": [ - "@.*1[ab]?$", - "@.*2[ab]?$", - "@.*[^/]3[ab]?$", - "@.*4[ab]?$", - "@.*5[ab]?$", - "@.*6[ab]?$", + +class Test_mtype_densities_from_composition: + @pytest.fixture(scope="session") + def class_tmpdir(self, tmpdir_factory): + """Create a session scoped temp dir using the class name""" + return tmpdir_factory.mktemp(type(self).__name__) + + @pytest.fixture(scope="session") + def annotation_path(self, class_tmpdir): + path = class_tmpdir.join("annotation.nrrd") + VoxelData( + np.array([[[101, 102, 103, 104, 105, 106]]], dtype=np.int32), + voxel_dimensions=[25] * 3, + ).save_nrrd(str(path)) + return path + + @pytest.fixture(scope="session") + def hierarchy_path(self, class_tmpdir): + path = class_tmpdir.join("hierarchy.json") + + hierarchy = { + "id": 315, + "acronym": "Isocortex", + "name": "Isocortex", + "children": [ + { + "id": 500, + "acronym": "MO", + "name": "Somatomotor areas", + "children": [ + { + "id": 101, + "acronym": "MO1", + "name": "Somatomotor areas, Layer 1", + "children": [], + }, + { + "id": 102, + "acronym": "MO2", + "name": "Somatomotor areas, Layer 2", + "children": [], + }, + { + "id": 103, + "acronym": "MO3", + "name": "Somatomotor areas, Layer 3", + "children": [], + }, + { + "id": 104, + "acronym": "MO4", + "name": "Somatomotor areas, Layer 4", + "children": [], + }, + { + "id": 105, + "acronym": "MO5", + "name": "Somatomotor areas, layer 5", + "children": [], + }, + { + "id": 106, + "acronym": "MO6", + "name": "Somatomotor areas, layer 6", + "children": [], + }, ], - "attribute": "acronym", - "with_descendants": True, }, - } - with open(path, "w", encoding="utf-8") as jsonfile: - json.dump(metadata, jsonfile, indent=1) + ], + } + with open(path, "w", encoding="utf-8") as jsonfile: + json.dump(hierarchy, jsonfile, indent=1, separators=(",", ": ")) + return path + + @pytest.fixture(scope="session") + def metadata_path(self, class_tmpdir): + path = class_tmpdir.join("metadata.json") + metadata = { + "region": { + "name": "Isocortex", + "query": "Isocortex", + "attribute": "acronym", + "with_descendants": True, + }, + "layers": { + "names": ["layer_1", "layer_2", "layer_3", "layer_4", "layer_5", "layer_6"], + "queries": [ + "@.*1[ab]?$", + "@.*2[ab]?$", + "@.*[^/]3[ab]?$", + "@.*4[ab]?$", + "@.*5[ab]?$", + "@.*6[ab]?$", + ], + "attribute": "acronym", + "with_descendants": True, + }, + } + with open(path, "w", encoding="utf-8") as jsonfile: + json.dump(metadata, jsonfile, indent=1) + + return path + + @pytest.fixture(scope="session") + def density(self): + return VoxelData( + np.array([[[0.3, 0.3, 0.3, 0.3, 0.3, 0.3]]], dtype=np.float32), + voxel_dimensions=[25] * 3, + ) + + @pytest.fixture(scope="session") + def density_path(self, density, class_tmpdir): + path = class_tmpdir.join("density.nrrd") + density.save_nrrd(str(path)) + + return path + + @pytest.fixture(scope="session") + def taxonomy(self): + return pd.DataFrame( + { + "mtype": ["L3_TPC:A", "L3_TPC:B", "L23_MC", "L4_TPC", "L4_LBC", "L4_UPC"], + "mClass": ["PYR", "PYR", "INT", "PYR", "INT", "PYR"], + "sClass": ["EXC", "EXC", "INH", "EXC", "INH", "EXC"], + }, + columns=["mtype", "mClass", "sClass"], + ) + + @pytest.fixture(scope="session") + def taxonomy_path(self, taxonomy, class_tmpdir): + """Creates a taxonomy file and returns its path""" + path = class_tmpdir.join("neurons-mtype-taxonomy.tsv") + taxonomy.to_csv(path, sep="\t") + + return path + + @pytest.fixture(scope="session") + def composition(self): + return pd.DataFrame( + { + "density": [51750.099, 14785.743, 2779.081, 62321.137, 2103.119, 25921.181], + "layer": ["layer_3", "layer_3", "layer_3", "layer_4", "layer_4", "layer_4"], + "mtype": ["L3_TPC:A", "L3_TPC:B", "L23_MC", "L4_TPC", "L4_LBC", "L4_UPC"], + }, + columns=["density", "layer", "mtype"], + ) + + @pytest.fixture(scope="session") + def composition_path(self, composition, class_tmpdir): + path = class_tmpdir.join("composition.yaml") + out_dict = {"neurons": []} + for row in composition.itertuples(): + out_dict["neurons"].append( + { + "density": row.density, + "traits": { + "layer": int(row.layer.replace("layer_", "")), + "mtype": row.mtype, + }, + } + ) + with open(path, "w", encoding="utf-8") as yaml_file: + yaml.dump(out_dict, yaml_file) - return path + return path - @pytest.fixture(scope="session") - def density(self): - return VoxelData( - np.array([[[0.3, 0.3, 0.3, 0.3, 0.3, 0.3]]], dtype=np.float32), - voxel_dimensions=[25] * 3, - ) + def test_load_neuronal_mtype_taxonomy(self, taxonomy, taxonomy_path): + pdt.assert_frame_equal(tested._load_neuronal_mtype_taxonomy(taxonomy_path), taxonomy) - @pytest.fixture(scope="session") - def density_path(self, density, class_tmpdir): - path = class_tmpdir.join("density.nrrd") - density.save_nrrd(str(path)) + def test_validate_mtype_taxonomy(self, taxonomy): + tested._validate_mtype_taxonomy(taxonomy) - return path + wrong_taxonomy = taxonomy.rename(columns={"sClass": "John"}) + with pytest.raises(AtlasDensitiesError): + tested._validate_mtype_taxonomy(wrong_taxonomy) - @pytest.fixture(scope="session") - def taxonomy(self): - return pd.DataFrame( - { - "mtype": ["L3_TPC:A", "L3_TPC:B", "L23_MC", "L4_TPC", "L4_LBC", "L4_UPC"], - "mClass": ["PYR", "PYR", "INT", "PYR", "INT", "PYR"], - "sClass": ["EXC", "EXC", "INH", "EXC", "INH", "EXC"], - }, - columns=["mtype", "mClass", "sClass"], - ) + wrong_taxonomy = taxonomy.drop(columns=["mtype"]) + with pytest.raises(AtlasDensitiesError): + tested._validate_mtype_taxonomy(wrong_taxonomy) - @pytest.fixture(scope="session") - def taxonomy_path(self, taxonomy, class_tmpdir): - """Creates a taxonomy file and returns its path""" - path = class_tmpdir.join("neurons-mtype-taxonomy.tsv") - taxonomy.to_csv(path, sep="\t") + wrong_taxonomy = deepcopy(taxonomy) + wrong_taxonomy.loc[1, "sClass"] = "OTHER" + with pytest.raises(AtlasDensitiesError): + tested._validate_mtype_taxonomy(wrong_taxonomy) - return path + def test_load_neuronal_mtype_composition(self, composition, composition_path): + pdt.assert_frame_equal( + tested._load_neuronal_mtype_composition(composition_path), composition + ) - @pytest.fixture(scope="session") - def composition(self): - return pd.DataFrame( - { - "density": [51750.099, 14785.743, 2779.081, 62321.137, 2103.119, 25921.181], - "layer": ["layer_3", "layer_3", "layer_3", "layer_4", "layer_4", "layer_4"], - "mtype": ["L3_TPC:A", "L3_TPC:B", "L23_MC", "L4_TPC", "L4_LBC", "L4_UPC"], - }, - columns=["density", "layer", "mtype"], - ) + def test_validate_density(self, density): + tested._validate_density(density) - @pytest.fixture(scope="session") - def composition_path(self, composition, class_tmpdir): - path = class_tmpdir.join("composition.yaml") - out_dict = {"neurons": []} - for row in composition.itertuples(): - out_dict["neurons"].append( - { - "density": row.density, - "traits": { - "layer": int(row.layer.replace("layer_", "")), - "mtype": row.mtype, - }, - } - ) - with open(path, "w", encoding="utf-8") as yaml_file: - yaml.dump(out_dict, yaml_file) + wrong_density = density.with_data(np.zeros_like(density.raw)) + with pytest.raises(AtlasDensitiesError): + tested._validate_density(wrong_density) - return path + wrong_density = deepcopy(density) + wrong_density.raw[0, 0, 2] = -10.0 + with pytest.raises(AtlasDensitiesError): + tested._validate_density(wrong_density) - def test_load_neuronal_mtype_taxonomy(self, taxonomy, taxonomy_path): - pdt.assert_frame_equal(tested._load_neuronal_mtype_taxonomy(taxonomy_path), taxonomy) + def test_validate_neuronal_mtype_composition(self, composition): + tested._validate_neuronal_mtype_composition(composition) - def test_validate_mtype_taxonomy(self, taxonomy): - tested._validate_mtype_taxonomy(taxonomy) + wrong_composition = composition.copy(deep=True) + wrong_composition[["density", 2]] = -1.0 + with pytest.raises(AtlasDensitiesError): + tested._validate_neuronal_mtype_composition(wrong_composition) - wrong_taxonomy = taxonomy.rename(columns={"sClass": "John"}) - with pytest.raises(AtlasDensitiesError): - tested._validate_mtype_taxonomy(wrong_taxonomy) + def test_check_taxonomy_composition_congruency(self, taxonomy, composition): + tested._check_taxonomy_composition_congruency(taxonomy, composition) - wrong_taxonomy = taxonomy.drop(columns=["mtype"]) - with pytest.raises(AtlasDensitiesError): - tested._validate_mtype_taxonomy(wrong_taxonomy) + with pytest.raises(AtlasDensitiesError): + tested._check_taxonomy_composition_congruency(taxonomy.drop([1]), composition) - wrong_taxonomy = deepcopy(taxonomy) - wrong_taxonomy.loc[1, "sClass"] = "OTHER" - with pytest.raises(AtlasDensitiesError): - tested._validate_mtype_taxonomy(wrong_taxonomy) + with pytest.raises(AtlasDensitiesError): + tested._check_taxonomy_composition_congruency(taxonomy, composition.drop([2])) - def test_load_neuronal_mtype_composition(self, composition, composition_path): - pdt.assert_frame_equal( - tested._load_neuronal_mtype_composition(composition_path), composition - ) + def test_create_from_composition( + self, + annotation_path, + hierarchy_path, + metadata_path, + density_path, + taxonomy_path, + composition_path, + class_tmpdir, + ): + output_dir = class_tmpdir.mkdir("output") - def test_validate_density(self, density): - tested._validate_density(density) - - wrong_density = density.with_data(np.zeros_like(density.raw)) - with pytest.raises(AtlasDensitiesError): - tested._validate_density(wrong_density) - - wrong_density = deepcopy(density) - wrong_density.raw[0, 0, 2] = -10.0 - with pytest.raises(AtlasDensitiesError): - tested._validate_density(wrong_density) - - def test_validate_neuronal_mtype_composition(self, composition): - tested._validate_neuronal_mtype_composition(composition) - - wrong_composition = composition.copy(deep=True) - wrong_composition[["density", 2]] = -1.0 - with pytest.raises(AtlasDensitiesError): - tested._validate_neuronal_mtype_composition(wrong_composition) - - def test_check_taxonomy_composition_congruency(self, taxonomy, composition): - tested._check_taxonomy_composition_congruency(taxonomy, composition) - - with pytest.raises(AtlasDensitiesError): - tested._check_taxonomy_composition_congruency(taxonomy.drop([1]), composition) - - with pytest.raises(AtlasDensitiesError): - tested._check_taxonomy_composition_congruency(taxonomy, composition.drop([2])) - - def test_create_from_composition( - self, - annotation_path, - hierarchy_path, - metadata_path, - density_path, - taxonomy_path, - composition_path, - class_tmpdir, - ): - output_dir = class_tmpdir.mkdir("output") - - runner = CliRunner() - - with runner.isolated_filesystem(temp_dir=class_tmpdir): - result = runner.invoke( - tested.create_from_composition, - [ - "--annotation-path", - annotation_path, - "--hierarchy-path", - hierarchy_path, - "--metadata-path", - metadata_path, - "--excitatory-neuron-density-path", - density_path, - "--taxonomy-path", - taxonomy_path, - "--composition-path", - composition_path, - "--output-dir", - output_dir, - ], - ) + runner = CliRunner() - assert result.exit_code == 0, result.output + with runner.isolated_filesystem(temp_dir=class_tmpdir): + result = runner.invoke( + tested.create_from_composition, + [ + "--annotation-path", + annotation_path, + "--hierarchy-path", + hierarchy_path, + "--metadata-path", + metadata_path, + "--excitatory-neuron-density-path", + density_path, + "--taxonomy-path", + taxonomy_path, + "--composition-path", + composition_path, + "--output-dir", + output_dir, + ], + ) - expected_filenames = { - "L3_TPC:A_densities.nrrd", - "L3_TPC:B_densities.nrrd", - "L4_TPC_densities.nrrd", - "L4_UPC_densities.nrrd", - } + assert result.exit_code == 0, result.output + + expected_filenames = { + "L3_TPC:A_densities.nrrd", + "L3_TPC:B_densities.nrrd", + "L4_TPC_densities.nrrd", + "L4_UPC_densities.nrrd", + } - filenames = set(os.listdir(output_dir)) + filenames = set(os.listdir(output_dir)) - assert filenames == expected_filenames + assert filenames == expected_filenames - for filename in filenames: - data = VoxelData.load_nrrd(str(Path(output_dir, filename))) - assert data.shape == (1, 1, 6) - assert not np.allclose(data.raw, 0.0) + for filename in filenames: + data = VoxelData.load_nrrd(str(Path(output_dir, filename))) + assert data.shape == (1, 1, 6) + assert not np.allclose(data.raw, 0.0) From 2b89cfdbd95c99900b3bbafa0e696fa5ca871570 Mon Sep 17 00:00:00 2001 From: MikeG Date: Fri, 22 Mar 2024 09:40:02 +0100 Subject: [PATCH 3/3] Cleanup (#71) * move functions from `helper` that were only used in `inhibitory_neuron_densities_optimization` * removed unused `with_descendants` argument from `_compute_region_cell_counts` * made `compute_region_volumes` faster * rename descendant_id_set -> descendant_ids --- atlas_densities/densities/fitting.py | 10 +- .../inhibitory_neuron_densities_helper.py | 123 +-------------- ...nhibitory_neuron_densities_optimization.py | 144 ++++++++++++++---- .../densities/measurement_to_density.py | 4 +- .../refined_inhibitory_neuron_densities.py | 8 +- atlas_densities/densities/utils.py | 26 ++-- ...test_inhibitory_neuron_densities_helper.py | 82 ---------- ..._inhibitory_neuron_density_optimization.py | 93 ++++++++++- .../test_optimization_initialization.py | 8 +- tests/densities/test_utils.py | 2 +- 10 files changed, 231 insertions(+), 269 deletions(-) diff --git a/atlas_densities/densities/fitting.py b/atlas_densities/densities/fitting.py index 3595160..b170ace 100644 --- a/atlas_densities/densities/fitting.py +++ b/atlas_densities/densities/fitting.py @@ -147,11 +147,11 @@ def fill_in_homogenous_regions( inhibitory_mask = homogenous_regions["cell_type"] == "inhibitory" for region_name in homogenous_regions[inhibitory_mask]["brain_region"]: - desc_id_set = hierarchy_info.at[region_name, "descendant_id_set"] + desc_id_set = hierarchy_info.at[region_name, "descendant_ids"] id_mask = [id_ in desc_id_set for id_ in hierarchy_info["id"]] for child_region_name in hierarchy_info[id_mask].index: region_mask = np.isin( - annotation, list(hierarchy_info.at[child_region_name, "descendant_id_set"]) + annotation, list(hierarchy_info.at[child_region_name, "descendant_ids"]) ) density = np.mean(neuron_density[region_mask]) densities.at[child_region_name, "inhibitory_neuron"] = density @@ -165,7 +165,7 @@ def fill_in_homogenous_regions( excitatory_mask = homogenous_regions["cell_type"] == "excitatory" for region_name in homogenous_regions[excitatory_mask]["brain_region"]: - desc_id_set = hierarchy_info.at[region_name, "descendant_id_set"] + desc_id_set = hierarchy_info.at[region_name, "descendant_ids"] id_mask = [id_ in desc_id_set for id_ in hierarchy_info["id"]] for child_region_name in hierarchy_info[id_mask].index: densities.at[child_region_name, "inhibitory_neuron"] = 0.0 @@ -232,7 +232,7 @@ def compute_average_intensities( } where each intensity array is of shape (W, H, D) and where the items of each slice list range in [0, W - 1]. - hierarchy_info: data frame with colums "descendant_id_set" (set[int]) and "brain_region" + hierarchy_info: data frame with colums "descendant_ids" (set[int]) and "brain_region" (str) and whose index is a list of region ids. See :fun:`atlas_densities.densities.utils.get_hierarchy_info`. @@ -260,7 +260,7 @@ def compute_average_intensities( region_count, ) for region_name in tqdm(result.index): - region_mask = np.isin(annotation, list(hierarchy_info.at[region_name, "descendant_id_set"])) + region_mask = np.isin(annotation, list(hierarchy_info.at[region_name, "descendant_ids"])) for marker, intensity in gene_marker_volumes.items(): result.at[region_name, marker.lower()] = compute_average_intensity( intensity["intensity"], region_mask, intensity["slices"] diff --git a/atlas_densities/densities/inhibitory_neuron_densities_helper.py b/atlas_densities/densities/inhibitory_neuron_densities_helper.py index 40f8ae0..632c10d 100644 --- a/atlas_densities/densities/inhibitory_neuron_densities_helper.py +++ b/atlas_densities/densities/inhibitory_neuron_densities_helper.py @@ -2,16 +2,12 @@ """ import logging -from typing import List, Optional, Set, Tuple +from typing import List import numpy as np import pandas as pd -from atlas_commons.typing import FloatArray - -from atlas_densities.exceptions import AtlasDensitiesError L = logging.getLogger(__name__) -MinMaxPair = Tuple[float, Optional[float]] def average_densities_to_cell_counts( @@ -34,11 +30,6 @@ def average_densities_to_cell_counts( :func:`atlas_densities.densities.utils.compute_region_volumes`. The volumes to be used are stored in the "volume" column. - - Note: Volumes can refer to volumes of entire regions, including all descendants - subregions, or volumes of single region identifiers, depending on the value of the - option `with_descendants` used when creating `region_volumes`. - Returns: data frame containing the cell counts of brain regions and their associated standard deviations. The data frame index is `average_densities.index` (brain region names) and @@ -54,35 +45,6 @@ def average_densities_to_cell_counts( return result -def resize_average_densities( - average_densities: pd.DataFrame, hierarchy_info: pd.DataFrame -) -> pd.DataFrame: - """ - Resize the `average_densities` data frame in such a way that the resized index coincides with - `hierarchy_info["brain_region"]`. - - Missing entries are set to ``np.nan``. - - average_densities: a data frame whose columns are described in - :func:`atlas_densities.densities.fitting.linear_fitting` containing the average - cell densities of brain regions and their associated standard deviations. Columns are - labelled by T and T_standard_deviation for various cell types T. The index of - `average_densities` is a list of region names. - hierarchy_info: data frame returned by - :func:`atlas_densities.densities.utils.get_hierarchy_info`. - - Returns: a data frame containing all the entries of `average_densities` and whose index is - `hierarchy_info["brain_region"]`. New entries are set to ``np.nan``. - - """ - resized = pd.DataFrame( - np.nan, index=hierarchy_info["brain_region"], columns=average_densities.columns - ) - resized.loc[average_densities.index] = average_densities - - return resized - - def get_cell_types(data_frame: pd.DataFrame) -> List[str]: """ Extract cell types from column labels. @@ -96,86 +58,3 @@ def get_cell_types(data_frame: pd.DataFrame) -> List[str]: """ return sorted({column.replace("_standard_deviation", "") for column in data_frame.columns}) - - -def check_region_counts_consistency( - region_counts: pd.DataFrame, hierarchy_info: pd.DataFrame, tolerance: float = 0.0 -) -> None: - """ - Check that cell counts considered as certain are consistent across the brain regions hierarchy. - - Args: - region_counts: data frame returned by - :fun:`densities.inhibitory_densities_helper.average_densities_to_cell_counts`. - A region is understood as a region of the brain hierarchy and includes all descendant - subregions. - hierarchy_info: data frame returned by - :func:`atlas_densities.densities.utils.get_hierarchy_info`. - tolerance: non-negative float number used as tolerance when comparing counts. - Defaults to 0.0. - - Raises: - AtlasDensitiesError if the sum of the cell counts of some leaf regions, all given with - certainty, exceeds the cell count of an ancestor, also given with certainty. - - """ - - def _check_descendants_consistency( - region_counts, hierarchy_info, region_name: str, id_set: Set[int], cell_type: str - ): - if region_counts.at[region_name, f"{cell_type}_standard_deviation"] == 0.0: - count = region_counts.at[region_name, cell_type] - descendants_names = hierarchy_info.loc[list(id_set), "brain_region"] - zero_std_mask = ( - region_counts.loc[descendants_names, f"{cell_type}_standard_deviation"] == 0.0 - ) - mask = zero_std_mask & ( - region_counts.loc[descendants_names, cell_type] > count + tolerance - ) - descendants_counts = region_counts.loc[descendants_names, cell_type][mask] - if not descendants_counts.empty: - raise AtlasDensitiesError( - f"Counts of {cell_type} cells in regions {list(descendants_names)} all exceed " - f"the count of its ancestor region {region_name} and each count is given " - f"for certain. The counts {descendants_counts} are all larger than " - f"{count}." - ) - names_with_certainty = descendants_names[zero_std_mask.to_list()].to_list() - leaf_names = [ - hierarchy_info.at[id_, "brain_region"] - for id_ in id_set - if len(hierarchy_info.at[id_, "descendant_id_set"]) == 1 - and hierarchy_info.at[id_, "brain_region"] in names_with_certainty - ] - leaves_count_sum = np.sum(region_counts.loc[leaf_names, cell_type]) - if leaves_count_sum > count + tolerance: - raise AtlasDensitiesError( - f"The sum of the counts of {cell_type} cells in leaf regions", - f" which are given with certainty, exceeds the count of the ancestor" - f" region {region_name}, also given with certainty: " - f"{leaves_count_sum} > {count}.", - ) - - cell_types = get_cell_types(region_counts) - for region_name, id_set in zip( - hierarchy_info["brain_region"], hierarchy_info["descendant_id_set"] - ): - for cell_type in cell_types: - _check_descendants_consistency( - region_counts, hierarchy_info, region_name, id_set, cell_type - ) - - -def replace_inf_with_none(bounds: FloatArray) -> List[MinMaxPair]: - """ - Replace the upper bounds equal to ``np.inf`` by None so as to comply with - the `bounds` interface of scipy.optimize.linprog. - - Args: - bounds: float array of shape (N, 2). Values in `bounds[..., 1]` can be - ``np.inf`` to indicate the absence of a constraining upper bound. - - Return: - list of pairs (min, max) where min is a float and max either a float or None. - """ - return [(float(min_), None if np.isinf(max_) else float(max_)) for min_, max_ in bounds] diff --git a/atlas_densities/densities/inhibitory_neuron_densities_optimization.py b/atlas_densities/densities/inhibitory_neuron_densities_optimization.py index 64d5987..389f2a7 100644 --- a/atlas_densities/densities/inhibitory_neuron_densities_optimization.py +++ b/atlas_densities/densities/inhibitory_neuron_densities_optimization.py @@ -36,10 +36,11 @@ variable deltas used in this module corresponds to the decision variables ``delta_{r, m}`` of the pdf file. """ +from __future__ import annotations import logging import warnings -from typing import Dict, List, Tuple +from typing import Dict, List, Optional, Set, Tuple import numpy as np import pandas as pd @@ -51,10 +52,7 @@ from atlas_densities.densities import utils from atlas_densities.densities.inhibitory_neuron_densities_helper import ( average_densities_to_cell_counts, - check_region_counts_consistency, get_cell_types, - replace_inf_with_none, - resize_average_densities, ) from atlas_densities.exceptions import AtlasDensitiesError, AtlasDensitiesWarning @@ -62,6 +60,114 @@ SKIP = np.inf # Used to signify that a delta variable of the linear program should be removed KEEP = np.nan # Used to signify that a delta variable should be added to the linear program +MinMaxPair = Tuple[float, Optional[float]] + + +def _resize_average_densities( + average_densities: pd.DataFrame, hierarchy_info: pd.DataFrame +) -> pd.DataFrame: + """ + Resize the `average_densities` data frame in such a way that the resized index coincides with + `hierarchy_info["brain_region"]`. + + Missing entries are set to ``np.nan``. + + average_densities: a data frame whose columns are described in + :func:`atlas_densities.densities.fitting.linear_fitting` containing the average + cell densities of brain regions and their associated standard deviations. Columns are + labelled by T and T_standard_deviation for various cell types T. The index of + `average_densities` is a list of region names. + hierarchy_info: data frame returned by + :func:`atlas_densities.densities.utils.get_hierarchy_info`. + + Returns: a data frame containing all the entries of `average_densities` and whose index is + `hierarchy_info["brain_region"]`. New entries are set to ``np.nan``. + + """ + resized = pd.DataFrame( + np.nan, index=hierarchy_info["brain_region"], columns=average_densities.columns + ) + resized.loc[average_densities.index] = average_densities + + return resized + + +def _check_region_counts_consistency( + region_counts: pd.DataFrame, hierarchy_info: pd.DataFrame, tolerance: float = 0.0 +) -> None: + """ + Check that cell counts considered as certain are consistent across the brain regions hierarchy. + + Args: + region_counts: data frame returned by + :fun:`densities.inhibitory_densities_helper.average_densities_to_cell_counts`. + A region is understood as a region of the brain hierarchy and includes all descendant + subregions. + hierarchy_info: data frame returned by + :func:`atlas_densities.densities.utils.get_hierarchy_info`. + tolerance: non-negative float number used as tolerance when comparing counts. + Defaults to 0.0. + + Raises: + AtlasDensitiesError if the sum of the cell counts of some leaf regions, all given with + certainty, exceeds the cell count of an ancestor, also given with certainty. + + """ + + def _check_descendants_consistency(region_name: str, id_set: Set[int], cell_type: str): + standard_deviation = f"{cell_type}_standard_deviation" + + if region_counts.at[region_name, standard_deviation] != 0.0: + return + + # count is certain, since standard_deviation is 0 + count = region_counts.at[region_name, cell_type] + descendants_names = hierarchy_info.loc[list(id_set), "brain_region"] + zero_std_mask = region_counts.loc[descendants_names, standard_deviation] == 0.0 + mask = zero_std_mask & (region_counts.loc[descendants_names, cell_type] > count + tolerance) + descendants_counts = region_counts.loc[descendants_names, cell_type][mask] + if not descendants_counts.empty: + raise AtlasDensitiesError( + f"Counts of {cell_type} cells in regions {list(descendants_names)} all exceed " + f"the count of its ancestor region {region_name} and each count is given " + f"for certain. The counts {descendants_counts} are all larger than " + f"{count}." + ) + + names_with_certainty = descendants_names[zero_std_mask.to_list()].to_list() + leaf_names = [ + hierarchy_info.at[id_, "brain_region"] + for id_ in id_set + if len(hierarchy_info.at[id_, "descendant_ids"]) == 1 + and hierarchy_info.at[id_, "brain_region"] in names_with_certainty + ] + leaves_count_sum = np.sum(region_counts.loc[leaf_names, cell_type]) + if leaves_count_sum > count + tolerance: + raise AtlasDensitiesError( + f"The sum of the counts of {cell_type} cells in leaf regions", + " which are given with certainty, exceeds the count of the ancestor" + f" region {region_name}, also given with certainty: " + f"{leaves_count_sum} > {count}.", + ) + + for _, region_name, id_set in hierarchy_info[["brain_region", "descendant_ids"]].itertuples(): + for cell_type in get_cell_types(region_counts): + _check_descendants_consistency(region_name, id_set, cell_type) + + +def _replace_inf_with_none(bounds: FloatArray) -> List[MinMaxPair]: + """ + Replace the upper bounds equal to ``np.inf`` by None so as to comply with + the `bounds` interface of scipy.optimize.linprog. + + Args: + bounds: float array of shape (N, 2). Values in `bounds[..., 1]` can be + ``np.inf`` to indicate the absence of a constraining upper bound. + + Return: + list of pairs (min, max) where min is a float and max either a float or None. + """ + return [(float(min_), None if np.isinf(max_) else float(max_)) for min_, max_ in bounds] def _compute_region_cell_counts( @@ -69,7 +175,6 @@ def _compute_region_cell_counts( density: FloatArray, voxel_volume: float, hierarchy_info: "pd.DataFrame", - with_descendants: bool = True, ) -> "pd.DataFrame": """ Compute the cell count of every 3D region in `annotation` labeled by a unique @@ -85,10 +190,6 @@ def _compute_region_cell_counts( This is (25 * 1e-6) ** 3 for an AIBS atlas nrrd file with 25um resolution. hierarchy_info: data frame returned by :func:`atlas_densities.densities.utils.get_hierarchy_info`. - with_descendants: if True, a computed cell count refers to the entire 3D region - designated by a region name, including all descendants subregions. Otherwise, the - computed cell count is the cell count of the 3D region labeled by the unique region - identifier. Defaults to True. Returns: DataFrame of the following form (values are fake): @@ -108,13 +209,6 @@ def _compute_region_cell_counts( index=hierarchy_info.index, ) - if with_descendants: - counts = [] - for id_set in tqdm(hierarchy_info["descendant_id_set"]): - count = result.loc[list(id_set), "cell_count"].sum() - counts.append(count) - result["cell_count"] = counts - return result @@ -174,7 +268,7 @@ def _zero_descendants(id_, cell_type, x_result, deltas, hierarchy_info): Zeroes the count of `cell_type` cells in every descendant regions of `id_`, including `id_`. """ - desc_ids = list(hierarchy_info.at[id_, "descendant_id_set"]) + desc_ids = list(hierarchy_info.at[id_, "descendant_ids"]) x_result.loc[desc_ids, cell_type] = 0.0 desc_names = hierarchy_info.loc[desc_ids, "brain_region"] deltas.loc[desc_names, cell_type] = 0.0 @@ -369,9 +463,7 @@ def check_constraint_3e( b_ub = [] variable_count = len(x_map) + len(deltas_map) cell_types = get_cell_types(region_counts) - for id_, region_name, set_ in zip( - hierarchy_info.index, hierarchy_info["brain_region"], hierarchy_info["descendant_id_set"] - ): + for id_, region_name, set_ in hierarchy_info[["brain_region", "descendant_ids"]].itertuples(): for cell_type in cell_types: if (region_name, cell_type) in deltas_map: constraint = np.zeros((variable_count,), dtype=float) @@ -490,9 +582,7 @@ def _compute_initial_cell_counts( L.info("Computing cell count estimates in every 3D region ...") region_counts = average_densities_to_cell_counts(average_densities, volumes) - # Detect cell counts inconsistency between a region and its descendants when - # estimates are deemed as certain. - check_region_counts_consistency(region_counts, hierarchy_info) + _check_region_counts_consistency(region_counts, hierarchy_info) L.info("Computing cell count estimates in atomic 3D region ...") volumes.drop("volume", axis=1, inplace=True) @@ -535,7 +625,7 @@ def _check_variables_consistency( """ cell_count_tolerance = 1e-2 # absolute tolerance to rule out round-off errors for region_name, id_, id_set in zip( - deltas.index, hierarchy_info.index, hierarchy_info["descendant_id_set"] + deltas.index, hierarchy_info.index, hierarchy_info["descendant_ids"] ): for cell_type in cell_types: if np.isfinite(deltas.loc[region_name, cell_type]): @@ -662,7 +752,7 @@ def create_inhibitory_neuron_densities( # pylint: disable=too-many-locals """ hierarchy_info = utils.get_hierarchy_info(RegionMap.from_dict(hierarchy), root=region_name) - average_densities = resize_average_densities(average_densities, hierarchy_info) + average_densities = _resize_average_densities(average_densities, hierarchy_info) L.info("Initialization of the linear program: started") region_counts, id_counts = _compute_initial_cell_counts( @@ -671,7 +761,7 @@ def create_inhibitory_neuron_densities( # pylint: disable=too-many-locals L.info("Retrieving overall neuron counts in atomic 3D regions ...") neuron_counts = _compute_region_cell_counts( - annotation, neuron_density, voxel_volume, hierarchy_info, with_descendants=False + annotation, neuron_density, voxel_volume, hierarchy_info ) assert np.all(neuron_counts["cell_count"] >= 0.0) @@ -712,7 +802,7 @@ def create_inhibitory_neuron_densities( # pylint: disable=too-many-locals c=c_row, A_ub=a_ub, b_ub=b_ub, - bounds=replace_inf_with_none(bounds), + bounds=_replace_inf_with_none(bounds), method="highs", ) if not result.success: diff --git a/atlas_densities/densities/measurement_to_density.py b/atlas_densities/densities/measurement_to_density.py index 65170e4..d300304 100644 --- a/atlas_densities/densities/measurement_to_density.py +++ b/atlas_densities/densities/measurement_to_density.py @@ -72,7 +72,7 @@ def compute_region_densities( The index is the sorted list of all region identifiers. """ densities = [] - descendants = hierarchy_info["descendant_id_set"] + descendants = hierarchy_info["descendant_ids"] for iset_ in tqdm(range(len(descendants))): mask = np.isin(annotation, list(descendants.iloc[iset_])) densities.append(np.sum(cell_density[mask]) / np.count_nonzero(mask)) @@ -245,7 +245,7 @@ def cell_count_per_slice_to_density( cell_counts_per_slice = measurements[mask_50um] hierarchy_info = hierarchy_info.set_index("brain_region") for index, row in cell_counts_per_slice.iterrows(): - id_set = hierarchy_info.loc[row["brain_region"], "descendant_id_set"] + id_set = hierarchy_info.loc[row["brain_region"], "descendant_ids"] average_slice_volume = voxel_volume * get_average_voxel_count_per_slice( id_set, annotation, thickness ) diff --git a/atlas_densities/densities/refined_inhibitory_neuron_densities.py b/atlas_densities/densities/refined_inhibitory_neuron_densities.py index 284e9ed..f8552cc 100644 --- a/atlas_densities/densities/refined_inhibitory_neuron_densities.py +++ b/atlas_densities/densities/refined_inhibitory_neuron_densities.py @@ -12,6 +12,7 @@ comparisons with the former approach of "A Cell Atlas for the Mouse Brain" by C. Eroe et al., 2018. """ +from __future__ import annotations import copy import logging @@ -24,6 +25,7 @@ from atlas_densities.densities.inhibitory_neuron_densities_helper import ( average_densities_to_cell_counts, + get_cell_types, ) from atlas_densities.densities.utils import compute_region_volumes, get_hierarchy_info @@ -368,7 +370,7 @@ def create_inhibitory_neuron_densities( # pylint: disable=too-many-locals annotation: AnnotationT, voxel_volume: float, neuron_density: FloatArray, - average_densities: "pd.DataFrame", + average_densities: pd.DataFrame, region_name: str = "root", ) -> Dict[str, FloatArray]: """ @@ -428,7 +430,7 @@ def create_inhibitory_neuron_densities( # pylint: disable=too-many-locals # modify cell densities recursively starting from the leaves L.info("Computing region cell counts ...") region_counts = average_densities_to_cell_counts(average_densities, volumes) - cell_types = {column.replace("_standard_deviation", "") for column in region_counts.columns} + cell_types = get_cell_types(region_counts) density_helper = VolumetricDensityHelper( annotation, @@ -436,7 +438,7 @@ def create_inhibitory_neuron_densities( # pylint: disable=too-many-locals neuron_density, region_counts, "gad67+", - list(cell_types - {"gad67+"}), + list(set(cell_types) - {"gad67+"}), ) L.info("Making leaf region cell counts consistent ...") diff --git a/atlas_densities/densities/utils.py b/atlas_densities/densities/utils.py index cd917f9..8648351 100644 --- a/atlas_densities/densities/utils.py +++ b/atlas_densities/densities/utils.py @@ -12,7 +12,6 @@ import scipy.misc import scipy.ndimage from atlas_commons.typing import AnnotationT, BoolArray, FloatArray -from tqdm import tqdm from atlas_densities.exceptions import AtlasDensitiesError, AtlasDensitiesWarning from atlas_densities.utils import copy_array @@ -458,7 +457,7 @@ def get_hierarchy_info( root: str = "Basic cell groups and regions", ) -> "pd.DataFrame": """ - Returns the name and the descendant_id_set of each region that can be found by `region_map`. + Returns the name and the descendant_ids of each region that can be found by `region_map`. Note: We assume that the hierarchy file has unique brain region names. @@ -471,24 +470,22 @@ def get_hierarchy_info( Returns: returns a dataframe with index a list of unique region ids and two columns (values are fake) - descendant_id_set brain_region + descendant_ids brain_region 1 {1, 3} "Cerebellum" 2 {2, 4, 10} "Isocortex" ... ... The index consists in the sorted list of the identifiers of every region recorded in - `region_map` under root. The column `descendant_id_set` holds for each region the set of + `region_map` under root. The column `descendant_ids` holds for each region the set of identifiers of the descendants including the region itself. `brain region` is the list of every region name. """ region_ids = list(region_map.find(root, attr="name", with_descendants=True)) region_ids.sort() - descendant_id_sets = [ - region_map.find(id_, attr="id", with_descendants=True) for id_ in region_ids - ] + descendant_ids = [region_map.find(id_, attr="id", with_descendants=True) for id_ in region_ids] 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}, + {"brain_region": region_names, "descendant_ids": descendant_ids}, index=region_ids, ) @@ -524,21 +521,20 @@ def compute_region_volumes( The latter column is created only if `hierarchy_info` has no `id_set` column, in which case its index is made of unique integer identifiers. """ - id_volumes = [] - for id_ in tqdm(hierarchy_info.index): - id_volumes.append(np.count_nonzero(annotation == id_) * voxel_volume) - result = pd.DataFrame( { "brain_region": hierarchy_info["brain_region"], - "id_volume": id_volumes, + "id_volume": 0.0, }, index=hierarchy_info.index, ) + ids, counts = np.unique(annotation, return_counts=True) + result["id_volume"] = pd.Series(counts * voxel_volume, index=ids) + volumes = [] - for id_ in tqdm(hierarchy_info.index): - id_set = hierarchy_info.loc[id_, "descendant_id_set"] + for id_ in hierarchy_info.index: + id_set = hierarchy_info.loc[id_, "descendant_ids"] volume = result.loc[list(id_set), "id_volume"].sum() volumes.append(volume) result["volume"] = volumes diff --git a/tests/densities/test_inhibitory_neuron_densities_helper.py b/tests/densities/test_inhibitory_neuron_densities_helper.py index f0330d0..b7b352d 100644 --- a/tests/densities/test_inhibitory_neuron_densities_helper.py +++ b/tests/densities/test_inhibitory_neuron_densities_helper.py @@ -98,42 +98,6 @@ def test_average_densities_to_cell_counts_with_reverse_index_order(): ) -def test_resize_average_densities(): - average_densities = pd.DataFrame( - { - "gad67+": [9.0, 12.0], - "gad67+_standard_deviation": [0.9, 1.2], - "pv+": [1.0, 2.0], - "pv+_standard_deviation": [0.1, 0.2], - "sst+": [3.0, 4.0], - "sst+_standard_deviation": [0.3, 0.4], - "vip+": [5.0, 6.0], - "vip+_standard_deviation": [0.5, 0.6], - }, - index=["A", "B"], - ) - hierarchy_info = pd.DataFrame( - {"brain_region": ["A", "B", "C"], "descendant_id_set": [{1}, {2}, {1, 2, 3}]}, - index=[1, 2, 3], - ) - - expected = pd.DataFrame( - { - "gad67+": [9.0, 12.0, np.nan], - "gad67+_standard_deviation": [0.9, 1.2, np.nan], - "pv+": [1.0, 2.0, np.nan], - "pv+_standard_deviation": [0.1, 0.2, np.nan], - "sst+": [3.0, 4.0, np.nan], - "sst+_standard_deviation": [0.3, 0.4, np.nan], - "vip+": [5.0, 6.0, np.nan], - "vip+_standard_deviation": [0.5, 0.6, np.nan], - }, - index=hierarchy_info["brain_region"], - ) - actual = tested.resize_average_densities(average_densities, hierarchy_info) - pdt.assert_frame_equal(actual, expected) - - def test_get_cell_types(): average_densities = pd.DataFrame( { @@ -150,49 +114,3 @@ def test_get_cell_types(): ) assert tested.get_cell_types(average_densities) == ["gad67+", "pv+", "sst+", "vip+"] - - -def test_check_region_counts_consistency(): - region_counts = pd.DataFrame( - { - "gad67+": [9.0, 12.0, 2.0], - "gad67+_standard_deviation": [0.9, 1.2, 1.0], - "pv+": [1.0, 1.0, 1.0], - "pv+_standard_deviation": [0.1, 0.2, 0.3], - "sst+": [3.0, 4.0, 5.0], - "sst+_standard_deviation": [0.3, 0.4, 0.5], - "vip+": [5.0, 6.0, 7.0], - "vip+_standard_deviation": [0.5, 0.6, 0.7], - }, - index=["A", "B", "C"], - ) - hierarchy_info = pd.DataFrame( - {"brain_region": ["A", "B", "C"], "descendant_id_set": [{1, 2, 3}, {2}, {3}]}, - index=[1, 2, 3], - ) - - # Uncertain or consistent - tested.check_region_counts_consistency(region_counts, hierarchy_info) - region_counts.loc["A", "gad67+_standard_deviation"] = 0.0 - tested.check_region_counts_consistency(region_counts, hierarchy_info) - - # Inconsistent parent-child count inequality - region_counts.loc["B", "gad67+_standard_deviation"] = 0.0 - with pytest.raises(AtlasDensitiesError, match=r"gad67\+.*exceed"): - tested.check_region_counts_consistency(region_counts, hierarchy_info) - - # The sum of the cell counts of the children regions B and C exceeds the count - # of the parent region A. - region_counts.loc["B", "gad67+_standard_deviation"] = 1.2 # restores valid value - region_counts.loc["A", "pv+_standard_deviation"] = 0.0 - region_counts.loc["B", "pv+_standard_deviation"] = 0.0 - region_counts.loc["C", "pv+_standard_deviation"] = 0.0 - with pytest.raises(AtlasDensitiesError, match=r"sum of the counts.*pv\+.*exceeds.*region A"): - tested.check_region_counts_consistency(region_counts, hierarchy_info) - - -def test_replace_inf_with_none(): - bounds = np.array([[0.1, 1.0], [-1.0, np.inf]]) - result = [(0.1, 1.0), (-1.0, None)] - - assert repr(result) == repr(tested.replace_inf_with_none(bounds)) diff --git a/tests/densities/test_inhibitory_neuron_density_optimization.py b/tests/densities/test_inhibitory_neuron_density_optimization.py index c7424a7..ba5276b 100644 --- a/tests/densities/test_inhibitory_neuron_density_optimization.py +++ b/tests/densities/test_inhibitory_neuron_density_optimization.py @@ -13,6 +13,90 @@ from atlas_densities.exceptions import AtlasDensitiesError +def test_resize_average_densities(): + average_densities = pd.DataFrame( + { + "gad67+": [9.0, 12.0], + "gad67+_standard_deviation": [0.9, 1.2], + "pv+": [1.0, 2.0], + "pv+_standard_deviation": [0.1, 0.2], + "sst+": [3.0, 4.0], + "sst+_standard_deviation": [0.3, 0.4], + "vip+": [5.0, 6.0], + "vip+_standard_deviation": [0.5, 0.6], + }, + index=["A", "B"], + ) + hierarchy_info = pd.DataFrame( + { + "brain_region": ["A", "B", "C"], + }, + index=[1, 2, 3], + ) + + expected = pd.DataFrame( + { + "gad67+": [9.0, 12.0, np.nan], + "gad67+_standard_deviation": [0.9, 1.2, np.nan], + "pv+": [1.0, 2.0, np.nan], + "pv+_standard_deviation": [0.1, 0.2, np.nan], + "sst+": [3.0, 4.0, np.nan], + "sst+_standard_deviation": [0.3, 0.4, np.nan], + "vip+": [5.0, 6.0, np.nan], + "vip+_standard_deviation": [0.5, 0.6, np.nan], + }, + index=hierarchy_info["brain_region"], + ) + actual = tested._resize_average_densities(average_densities, hierarchy_info) + pdt.assert_frame_equal(actual, expected) + + +def test_check_region_counts_consistency(): + region_counts = pd.DataFrame( + { + "gad67+": [9.0, 12.0, 2.0], + "gad67+_standard_deviation": [0.9, 1.2, 1.0], + "pv+": [1.0, 1.0, 1.0], + "pv+_standard_deviation": [0.1, 0.2, 0.3], + "sst+": [3.0, 4.0, 5.0], + "sst+_standard_deviation": [0.3, 0.4, 0.5], + "vip+": [5.0, 6.0, 7.0], + "vip+_standard_deviation": [0.5, 0.6, 0.7], + }, + index=["A", "B", "C"], + ) + hierarchy_info = pd.DataFrame( + {"brain_region": ["A", "B", "C"], "descendant_ids": [{1, 2, 3}, {2}, {3}]}, + index=[1, 2, 3], + ) + + # Uncertain or consistent + tested._check_region_counts_consistency(region_counts, hierarchy_info) + region_counts.loc["A", "gad67+_standard_deviation"] = 0.0 + tested._check_region_counts_consistency(region_counts, hierarchy_info) + + # Inconsistent parent-child count inequality + region_counts.loc["B", "gad67+_standard_deviation"] = 0.0 + with pytest.raises(AtlasDensitiesError, match=r"gad67\+.*exceed"): + tested._check_region_counts_consistency(region_counts, hierarchy_info) + + # The sum of the cell counts of the children regions B and C exceeds the count + # of the parent region A. + region_counts.loc["B", "gad67+_standard_deviation"] = 1.2 # restores valid value + region_counts.loc["A", "pv+_standard_deviation"] = 0.0 + region_counts.loc["B", "pv+_standard_deviation"] = 0.0 + region_counts.loc["C", "pv+_standard_deviation"] = 0.0 + with pytest.raises(AtlasDensitiesError, match=r"sum of the counts.*pv\+.*exceeds.*region A"): + tested._check_region_counts_consistency(region_counts, hierarchy_info) + + +def test_replace_inf_with_none(): + bounds = np.array([[0.1, 1.0], [-1.0, np.inf]]) + result = [(0.1, 1.0), (-1.0, None)] + + assert repr(result) == repr(tested._replace_inf_with_none(bounds)) + + def test_create_volumetric_densities(): annotation = np.array([[[1, 1, 2]]], dtype=int) neuron_density = np.array([[[0.5, 0.5, 1.0]]]) @@ -342,13 +426,6 @@ def get_hierarchy_info(): "Lobule II, Purkinje layer", "Lobule II, molecular layer", ], - "descendant_id_set": [ - {920, 976, 10708, 10709, 10710}, - {976, 10708, 10709, 10710}, - {10708}, - {10709}, - {10710}, - ], }, index=[920, 976, 10708, 10709, 10710], ) @@ -363,7 +440,7 @@ def test_compute_region_cell_counts(): cell_counts = pd.DataFrame( { "brain_region": hierarchy_info["brain_region"], - "cell_count": 2 * np.array([5.0, 4.0, 1.0, 1.0, 1.0]), + "cell_count": 2 * np.array([1.0, 1.0, 1.0, 1.0, 1.0]), }, index=hierarchy_info.index, ) diff --git a/tests/densities/test_optimization_initialization.py b/tests/densities/test_optimization_initialization.py index abf2540..f15bdc3 100644 --- a/tests/densities/test_optimization_initialization.py +++ b/tests/densities/test_optimization_initialization.py @@ -16,7 +16,7 @@ @pytest.fixture def counts_data_1(): hierarchy_info = pd.DataFrame( - {"brain_region": ["A", "B"], "descendant_id_set": [{1}, {2}]}, index=[1, 2] + {"brain_region": ["A", "B"], "descendant_ids": [{1}, {2}]}, index=[1, 2] ) region_counts = pd.DataFrame( { @@ -162,7 +162,7 @@ def test_set_known_values_3(counts_data_1): def counts_data_2(): hierarchy_info = pd.DataFrame( # A is now a parent region of B - {"brain_region": ["A", "B"], "descendant_id_set": [{1, 2}, {2}]}, + {"brain_region": ["A", "B"], "descendant_ids": [{1, 2}, {2}]}, index=[1, 2], ) region_counts = pd.DataFrame( @@ -335,7 +335,7 @@ def neuron_counts(): def counts_data_3(): hierarchy_info = pd.DataFrame( - {"brain_region": ["A", "B", "C"], "descendant_id_set": [{1, 2, 3}, {2}, {3}]}, + {"brain_region": ["A", "B", "C"], "descendant_ids": [{1, 2, 3}, {2}, {3}]}, index=[1, 2, 3], ) region_counts = pd.DataFrame( @@ -571,7 +571,7 @@ def test_create_bounds_2(): @pytest.fixture def aub_and_bub_data(): hierarchy_info = pd.DataFrame( - {"brain_region": ["A", "B", "C"], "descendant_id_set": [{1, 2, 3}, {2}, {3}]}, + {"brain_region": ["A", "B", "C"], "descendant_ids": [{1, 2, 3}, {2}, {3}]}, index=[1, 2, 3], ) x_map, deltas_map = expected_maps() diff --git a/tests/densities/test_utils.py b/tests/densities/test_utils.py index 5680a7f..a370327 100644 --- a/tests/densities/test_utils.py +++ b/tests/densities/test_utils.py @@ -353,7 +353,7 @@ def get_hierarchy_info(): "Lobule II, Purkinje layer", "Lobule II, molecular layer", ], - "descendant_id_set": [ + "descendant_ids": [ {920, 976, 10708, 10709, 10710}, {976, 10708, 10709, 10710}, {10708},