From 7783705a4fb3c5acbaeee9a4660c1ea987d6e752 Mon Sep 17 00:00:00 2001
From: Zisis Eleftherios <eleftherios.zisis@epfl.ch>
Date: Wed, 3 Apr 2024 16:03:21 +0200
Subject: [PATCH] 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": [