From 7f3afca2d7fc698841590d31806898fb835d5af4 Mon Sep 17 00:00:00 2001
From: "Paul M. Sutter" <pmsutter@gmail.com>
Date: Sat, 8 Jun 2024 14:04:59 +0100
Subject: [PATCH] fixed bug in adjacency reading

---
 c_source/zobov/voztie.c               |  2 ++
 python_source/voidUtil/catalogUtil.py | 29 +++++++++++++++++++++++----
 2 files changed, 27 insertions(+), 4 deletions(-)

diff --git a/c_source/zobov/voztie.c b/c_source/zobov/voztie.c
index 0f67b98..ae21042 100644
--- a/c_source/zobov/voztie.c
+++ b/c_source/zobov/voztie.c
@@ -179,6 +179,7 @@ int main(int argc, char *argv[]) {
   }
   printf("\n");
 
+/*
   // PMS : remove mock galaxies and anything adjacent to a mock galaxy
   printf("\nRemoving mock galaxies...\n"); 
  
@@ -221,6 +222,7 @@ int main(int argc, char *argv[]) {
   printf("There are %d galaxies remaining.\n", mockIndex-numRemoved);
 
   // END PMS
+*/
 
   npnotdone = 0; avgnadj = 0.; avgvol = 0.;
   for (p=0;p<np;p++) {
diff --git a/python_source/voidUtil/catalogUtil.py b/python_source/voidUtil/catalogUtil.py
index dcfb9a8..85a6e3e 100644
--- a/python_source/voidUtil/catalogUtil.py
+++ b/python_source/voidUtil/catalogUtil.py
@@ -435,7 +435,8 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True,
                                 y = partData[i][1],
                                 z = partData[i][2],
                                 volume = 0,
-                                adjs = 0,
+                                nadjs = 0,
+                                adjs = [],
                                 ra = extraData[i][0],
                                 dec = extraData[i][1],
                                 redshift = extraData[i][2],
@@ -467,10 +468,29 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True,
     print(" Loading adjacencies...")
     adjFile = sampleDir+"adj_"+sample.fullName+".dat"
     with open(adjFile, mode="rb") as File:
-      numPart = np.fromfile(adjFile, dtype=np.int32,count=1)[0]
+      numPart = np.fromfile(File, dtype=np.int32,count=1)[0]
+
+      # this the total number of adjancies per particle
+      nadjPerPart = np.fromfile(File, dtype=np.int32,count=numPart)
+      
+      # but the file only stores one half of each pair, so we need to match
       for p in range(numPart):
-        numAdjs = np.fromfile(File, dtype=np.int32,count=1)[0]
-        catalog.part[p].adjs = np.fromfile(File, dtype=np.int32, count=numAdjs)
+        nin = np.fromfile(File, dtype=np.int32, count=1)[0]
+        
+        for n in range(nin):
+          pAdj = np.fromfile(File, dtype=np.int32, count=1)[0]
+
+          if (p < pAdj):
+            catalog.part[p].adjs.append(pAdj)
+            catalog.part[pAdj].adjs.append(p)
+
+    print(" Sanity checking adjacenies...")
+    for p in range(numPart):
+      catalog.part[p].nadjs = len(catalog.part[p].adjs)
+      nHave = len(catalog.part[p].adjs)
+      nExpected = nadjPerPart[p]
+      if (nHave != nExpected):
+        print("   Error for particle %d: Have %d adj, expected %d" % (p, nHave, nExpected)) 
 
     print(" Loading zone-void membership info...")
     zoneFile = sampleDir+"/voidZone_"+sample.fullName+".dat"
@@ -502,6 +522,7 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True,
           partID = np.fromfile(File, dtype=np.int32,count=1)[0]
           catalog.zones2Parts[iZ].partIDs.append(partID)
 
+    print("Done loading catalog.")
   return catalog