fixed bug in adjacency reading

This commit is contained in:
Paul M. Sutter 2024-06-08 14:04:59 +01:00
parent ad752f21e8
commit 7f3afca2d7
2 changed files with 27 additions and 4 deletions

View file

@ -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++) {

View file

@ -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