From 22d314e2fc5d1eba26ec0183fd67b8486e616d63 Mon Sep 17 00:00:00 2001 From: Zisis Eleftherios Date: Wed, 27 Mar 2024 11:40:36 +0100 Subject: [PATCH 1/6] Use voxcell ValueToIndexVoxels --- .../mtype_densities_from_map/create.py | 26 +++++++++---------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/atlas_densities/densities/mtype_densities_from_map/create.py b/atlas_densities/densities/mtype_densities_from_map/create.py index 936005d..50f1b80 100644 --- a/atlas_densities/densities/mtype_densities_from_map/create.py +++ b/atlas_densities/densities/mtype_densities_from_map/create.py @@ -18,6 +18,7 @@ from atlas_commons.typing import FloatArray from joblib import Parallel, delayed from tqdm import tqdm +from voxcell.voxel_data import ValueToIndexVoxels from atlas_densities.densities.mtype_densities_from_map.utils import ( _check_probability_map_consistency, @@ -95,34 +96,33 @@ def create_from_probability_map( # pylint: disable=too-many-arguments .set_index("acronym") .loc[probability_map.index.get_level_values("region")] .reset_index()[["region", "id"]] + .drop_duplicates(subset="region") + .reset_index(drop=True) ) - region_acronyms = set(region_info.region) - - region_masks = { - region_acronym: annotation.raw == region_id - for _, region_acronym, region_id in region_info.itertuples() - } + region_index = ValueToIndexVoxels(annotation.raw) # ensure output directory exists Path(output_dirpath).mkdir(exist_ok=True, parents=True) def _create_densities_for_metype(metype: str) -> Optional[Tuple[str, str]]: coefficients: Dict[str, Dict[str, Any]] = {} - for region_acronym in region_acronyms: + for region_acronym in region_info.region: coefficients[region_acronym] = { molecular_type: probability_map.at[(region_acronym, molecular_type), metype] for molecular_type in list(molecular_type_densities.keys()) if (region_acronym, molecular_type) in probability_map.index } - metype_density = np.zeros(annotation.shape, dtype=float) - for region_acronym in region_acronyms: - region_mask = region_masks[region_acronym] + metype_density = region_index.ravel(np.zeros(annotation.shape, dtype=float)) + for region_acronym, region_id in region_info.itertuples(index=False): + + region_indices = region_index.value_to_1d_indices(region_id) for molecular_type, coefficient in coefficients[region_acronym].items(): if coefficient <= 0.0: continue - density = molecular_type_densities[molecular_type] - metype_density[region_mask] += density[region_mask] * coefficient + density = region_index.ravel(molecular_type_densities[molecular_type]) + metype_density[region_indices] += density[region_indices] * coefficient + metype_density = np.reshape(metype_density, annotation.shape, order=region_index._order) if np.any(metype_density): # save density file @@ -137,7 +137,7 @@ def _create_densities_for_metype(metype: str) -> Optional[Tuple[str, str]]: returns = Parallel(n_jobs=n_jobs, return_as="generator")( delayed(_create_densities_for_metype)(metype) for metype in probability_map.columns ) - + #returns = [_create_densities_for_metype(metype) for metype in probability_map.columns[0:10]] # construct metadata output_legend: Dict[str, Dict[str, str]] = defaultdict(dict) for return_value in tqdm(returns, total=len(probability_map.columns)): From 535606cd426f317409188678db577dcae29c5747 Mon Sep 17 00:00:00 2001 From: Zisis Eleftherios Date: Wed, 27 Mar 2024 13:21:05 +0100 Subject: [PATCH 2/6] Remove comment --- atlas_densities/densities/mtype_densities_from_map/create.py | 1 - 1 file changed, 1 deletion(-) diff --git a/atlas_densities/densities/mtype_densities_from_map/create.py b/atlas_densities/densities/mtype_densities_from_map/create.py index 50f1b80..cd7fce0 100644 --- a/atlas_densities/densities/mtype_densities_from_map/create.py +++ b/atlas_densities/densities/mtype_densities_from_map/create.py @@ -137,7 +137,6 @@ def _create_densities_for_metype(metype: str) -> Optional[Tuple[str, str]]: returns = Parallel(n_jobs=n_jobs, return_as="generator")( delayed(_create_densities_for_metype)(metype) for metype in probability_map.columns ) - #returns = [_create_densities_for_metype(metype) for metype in probability_map.columns[0:10]] # construct metadata output_legend: Dict[str, Dict[str, str]] = defaultdict(dict) for return_value in tqdm(returns, total=len(probability_map.columns)): From 7783705a4fb3c5acbaeee9a4660c1ea987d6e752 Mon Sep 17 00:00:00 2001 From: Zisis Eleftherios Date: Wed, 3 Apr 2024 16:03:21 +0200 Subject: [PATCH 3/6] Use latest voxcell version with unravel method --- .../densities/mtype_densities_from_map/create.py | 14 +++++++++----- setup.py | 2 +- 2 files changed, 10 insertions(+), 6 deletions(-) diff --git a/atlas_densities/densities/mtype_densities_from_map/create.py b/atlas_densities/densities/mtype_densities_from_map/create.py index cd7fce0..a260398 100644 --- a/atlas_densities/densities/mtype_densities_from_map/create.py +++ b/atlas_densities/densities/mtype_densities_from_map/create.py @@ -99,7 +99,7 @@ def create_from_probability_map( # pylint: disable=too-many-arguments .drop_duplicates(subset="region") .reset_index(drop=True) ) - region_index = ValueToIndexVoxels(annotation.raw) + annotation_index = ValueToIndexVoxels(annotation.raw) # ensure output directory exists Path(output_dirpath).mkdir(exist_ok=True, parents=True) @@ -113,16 +113,20 @@ def _create_densities_for_metype(metype: str) -> Optional[Tuple[str, str]]: if (region_acronym, molecular_type) in probability_map.index } - metype_density = region_index.ravel(np.zeros(annotation.shape, dtype=float)) + # perform the manipulation in the 1d flat array + metype_density = np.zeros(np.prod(annotation.shape), dtype=float) + for region_acronym, region_id in region_info.itertuples(index=False): - region_indices = region_index.value_to_1d_indices(region_id) + region_indices = annotation_index.value_to_1d_indices(region_id) for molecular_type, coefficient in coefficients[region_acronym].items(): if coefficient <= 0.0: continue - density = region_index.ravel(molecular_type_densities[molecular_type]) + density = annotation_index.ravel(molecular_type_densities[molecular_type]) metype_density[region_indices] += density[region_indices] * coefficient - metype_density = np.reshape(metype_density, annotation.shape, order=region_index._order) + + # reshape the 1d metype_density array ack to the annotation's shape + metype_density = annotation_index.unravel(metype_density) if np.any(metype_density): # save density file diff --git a/setup.py b/setup.py index 308d6ee..7bf8bc6 100644 --- a/setup.py +++ b/setup.py @@ -28,7 +28,7 @@ # from the HiGHS library. We use the "highs" method in the densities module. "scipy>=1.6.0", "tqdm>=4.44.1", - "voxcell>=3.1.7", + "voxcell>=3.1.8", # ValueToIndexVoxels ravel/unravell ], extras_require={ "tests": [ From 72904b13d92831e0cf8a6de79192982897d5610a Mon Sep 17 00:00:00 2001 From: Zisis Eleftherios Date: Wed, 3 Apr 2024 16:05:41 +0200 Subject: [PATCH 4/6] Typo --- atlas_densities/densities/mtype_densities_from_map/create.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/atlas_densities/densities/mtype_densities_from_map/create.py b/atlas_densities/densities/mtype_densities_from_map/create.py index a260398..9995c0b 100644 --- a/atlas_densities/densities/mtype_densities_from_map/create.py +++ b/atlas_densities/densities/mtype_densities_from_map/create.py @@ -125,7 +125,7 @@ def _create_densities_for_metype(metype: str) -> Optional[Tuple[str, str]]: density = annotation_index.ravel(molecular_type_densities[molecular_type]) metype_density[region_indices] += density[region_indices] * coefficient - # reshape the 1d metype_density array ack to the annotation's shape + # reshape the 1d metype_density array back to the annotation's shape metype_density = annotation_index.unravel(metype_density) if np.any(metype_density): From 70135e530d423566c89d1ed1a2c7903bc625a877 Mon Sep 17 00:00:00 2001 From: Zisis Eleftherios Date: Wed, 3 Apr 2024 16:08:23 +0200 Subject: [PATCH 5/6] Fix typo --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 7bf8bc6..ee9c873 100644 --- a/setup.py +++ b/setup.py @@ -28,7 +28,7 @@ # from the HiGHS library. We use the "highs" method in the densities module. "scipy>=1.6.0", "tqdm>=4.44.1", - "voxcell>=3.1.8", # ValueToIndexVoxels ravel/unravell + "voxcell>=3.1.8", # ValueToIndexVoxels ravel/unravel ], extras_require={ "tests": [ From 576c7a3da985217fcb11a1fa6435efa8ac3e255f Mon Sep 17 00:00:00 2001 From: Zisis Eleftherios Date: Wed, 3 Apr 2024 16:24:45 +0200 Subject: [PATCH 6/6] Revert removed line --- atlas_densities/densities/mtype_densities_from_map/create.py | 1 + 1 file changed, 1 insertion(+) diff --git a/atlas_densities/densities/mtype_densities_from_map/create.py b/atlas_densities/densities/mtype_densities_from_map/create.py index 9995c0b..ff5eb02 100644 --- a/atlas_densities/densities/mtype_densities_from_map/create.py +++ b/atlas_densities/densities/mtype_densities_from_map/create.py @@ -141,6 +141,7 @@ def _create_densities_for_metype(metype: str) -> Optional[Tuple[str, str]]: returns = Parallel(n_jobs=n_jobs, return_as="generator")( delayed(_create_densities_for_metype)(metype) for metype in probability_map.columns ) + # construct metadata output_legend: Dict[str, Dict[str, str]] = defaultdict(dict) for return_value in tqdm(returns, total=len(probability_map.columns)):