From 5ea91f12f58e6abef735d332d1b06d64ea5a1be7 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 10 Oct 2013 13:00:05 -0500 Subject: [PATCH] temporary commit to check adjacencies in voids --- c_tools/stacking/pruneVoids.cpp | 102 +++++++++++++++++- c_tools/stacking/pruneVoids.ggo | 2 + .../void_python_tools/backend/launchers.py | 1 + 3 files changed, 102 insertions(+), 3 deletions(-) diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index 6f8f60f..d3f0f15 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -57,6 +57,8 @@ using namespace std; typedef struct partStruct { float x, y, z, vol; + int nadj, ncnt; + int *adj; } PART; typedef struct zoneStruct { @@ -71,8 +73,10 @@ typedef struct voidZoneStruct { typedef struct voidStruct { float vol, coreDens, zoneVol, densCon, voidProb, radius; + float rescaledCoreDens; int voidID, numPart, numZones, coreParticle, zoneNumPart; float maxRadius, nearestMock, centralDen, redshift, redshiftInMpc; + float nearestMockFromCore, nearestGalFromCore; float nearestEdge; float center[3], barycenter[3]; int accepted; @@ -393,7 +397,56 @@ int main(int argc, char **argv) { part[p].vol = temp[0]; } fclose(fp); - free(temp); + +// GUILHEM LOOK HERE + // and finally finally adjacencies + printf(" Loading particle adjacencies...\n"); + fp = fopen(args.partAdj_arg, "r"); + fread(&mask_index, 1, 4, fp); + if (mask_index != mockIndex) { + printf("NON-MATCHING MOCK INDICES!? %d %d\n", mask_index, mockIndex); + exit(-1); + } + int tempInt; + for (p = 0; p < mask_index; p++) { + fread(&tempInt, 1, 4, fp); + part[p].nadj = tempInt; + part[p].ncnt = 0; + if (part[p].nadj > 0) + part[p].adj = (int *) malloc(part[p].nadj * sizeof(int)); + } + for (p = 0; p < mask_index; p++) { + fread(&tempInt, 1, 4, fp); + int nin = tempInt; + if (nin > 0) { + for (int nAdj = 0; nAdj < nin; nAdj++) { + int tempAdj; + fread(&tempAdj, 1, 4, fp); + + // this bit has been readjusted just in case we are + // accidentally still linking to mock particles + + //if (tempAdj < mask_index) { + assert(p < tempAdj); + //if (part[p].ncnt == part[p].nadj) { + // printf("OVERFLOW %d\n", p); + //} else if (part[tempAdj].ncnt == part[tempAdj].nadj) { + // printf("OVERFLOW %d\n", tempAdj); + //} else { + part[p].adj[part[p].ncnt] = tempAdj; + part[p].ncnt++; + if (tempAdj < mask_index) { + part[tempAdj].adj[part[tempAdj].ncnt] = p; + part[tempAdj].ncnt++; + } + //} + //} + } +//printf("ADJ %d %d %d %d %d\n", p, nin, part[p].nadj, nAdj, tempInt); + } + } +// END GUILHEM LOOK HERE + fclose(fp); // load voids *again* using Guilhem's code so we can get tree if (!args.isObservation_flag) { @@ -466,6 +519,19 @@ int main(int argc, char **argv) { voidPart[i].z = part[partID].z; voidPart[i].vol = part[partID].vol; +// GUILHEM LOOK HERE + // testing for edge contamination + if (part[partID].vol < 1.e-27) { + printf("CONTAMINATED!! %d %d\n", iVoid, partID); + } else { + //printf("NORMAL!! %d %d %e\n", iVoid, partID, part[partID].vol); + } + for (int iAdj = 0; iAdj < part[partID].ncnt; iAdj++) { + if (part[partID].adj[iAdj] > mockIndex) { + printf("CONTAMINATED!! %d %d %d\n", iVoid, partID, iAdj); + } + } +// END GUILHEM LOOK HERE i++; } } @@ -538,6 +604,32 @@ int main(int argc, char **argv) { voids[iVoid].centralDen = numCentral / (volNorm*4./3. * M_PI * pow(centralRad, 3.)); + coreParticle = voids[iVoid].coreParticle; + voids[iVoid].rescaledCoreDens = 1./part[coreParticle].vol*volNorm; + // compute distance from core to nearest mock + minDist = 1.e99; + for (p = mockIndex; p < numPartTot; p++) { + dist[0] = part[coreParticle].x - part[p].x; + dist[1] = part[coreParticle].y - part[p].y; + dist[2] = part[coreParticle].z - part[p].z; + + dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); + if (dist2 < minDist) minDist = dist2; + } + voids[iVoid].nearestMockFromCore = sqrt(minDist); + + // compute distance from core to nearest mock + minDist = 1.e99; + for (p = 0; p < mockIndex; p++) { + dist[0] = part[coreParticle].x - part[p].x; + dist[1] = part[coreParticle].y - part[p].y; + dist[2] = part[coreParticle].z - part[p].z; + + dist2 = pow(dist[0],2) + pow(dist[1],2) + pow(dist[2],2); + if (dist2 < minDist && dist2 > 1.e-10) minDist = dist2; + } + voids[iVoid].nearestGalFromCore = sqrt(minDist); + // compute maximum extent /* if (args.isObservation_flag) { @@ -994,9 +1086,13 @@ void outputVoids(string outputDir, string sampleName, string prefix, outVoid.barycenter[1], outVoid.barycenter[2]); - fprintf(fpDistances, "%d %e\n", + fprintf(fpDistances, "%d %e %e %e %e %e\n", outVoid.voidID, - outVoid.nearestMock); + outVoid.nearestMock, + outVoid.radius, + outVoid.rescaledCoreDens, + outVoid.nearestMockFromCore, + outVoid.nearestGalFromCore); fprintf(fpCenters, "%.2f %.2f %.2f %.2f %.2f %.5f %.2f %d %f %d %d %d %d %.2f\n", outCenter[0], diff --git a/c_tools/stacking/pruneVoids.ggo b/c_tools/stacking/pruneVoids.ggo index a61f3a2..4d0c8ba 100644 --- a/c_tools/stacking/pruneVoids.ggo +++ b/c_tools/stacking/pruneVoids.ggo @@ -14,6 +14,8 @@ option "void2Zone" - "Zone file from ZOBOV" string required option "partVol" - "Particle volume file from ZOBOV" string required +option "partAdj" - "Adjacency file from ZOBOV" string required + option "zone2Part" - "Particle file from ZOBOV" string required option "mockIndex" - "Beginning index of mock particles" int required diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index b5e6a83..434cb91 100644 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -342,6 +342,7 @@ def launchPrune(sample, binPath, cmd += " --void2Zone="+zobovDir+"/voidZone_"+str(sampleName)+".dat" cmd += " --zone2Part=" + zobovDir+"/voidPart_"+str(sampleName)+".dat" cmd += " --partVol=" + zobovDir+"/vol_"+str(sampleName)+".dat" + cmd += " --partAdj=" + zobovDir+"/adj_"+str(sampleName)+".dat" cmd += " --extraInfo=" + zobovDir+"/zobov_slice_"+str(sampleName)+\ ".par" cmd += " --tolerance=1.0"