From 3c8d7b2fb29d86a9421f074454a2714b7691606c Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 20 Dec 2012 13:38:14 -0600 Subject: [PATCH 01/16] significant speedups to void matching algorithm --- c_tools/analysis/voidOverlap.cpp | 122 ++++++++++++++----------------- 1 file changed, 55 insertions(+), 67 deletions(-) diff --git a/c_tools/analysis/voidOverlap.cpp b/c_tools/analysis/voidOverlap.cpp index 0d1df14..37f9a9f 100644 --- a/c_tools/analysis/voidOverlap.cpp +++ b/c_tools/analysis/voidOverlap.cpp @@ -88,6 +88,8 @@ int main(int argc, char **argv) { int closestMatchID; float closestMatchDist; float commonVolRatio; + MATCHPROPS newMatch; + int MAX_MATCHES = 20; CATALOG catalog1, catalog2; @@ -134,38 +136,44 @@ int main(int argc, char **argv) { printf("Will assume z-direction is periodic.\n"); } - printf(" Determining overlap...\n"); - -/* - // just use centers + // find closest voids + printf(" Finding nearest matches...\n"); for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) { - closestMatchDist = 1.e99; for (iVoid2 = 0; iVoid2 < catalog2.numVoids; iVoid2++) { rdist = getDist(catalog1, catalog2, iVoid1, iVoid2, periodicX, periodicY, periodicZ); - if (rdist < closestMatchDist) { - closestMatchID = iVoid2; - closestMatchDist = rdist; + + newMatch.matchID = iVoid2; + newMatch.commonVol = 0; + newMatch.dist = rdist; + + if (catalog1.voids[iVoid1].matches.size() < MAX_MATCHES) { + catalog1.voids[iVoid1].matches.push_back(newMatch); + } else { + // find the farthest match + float farthestMatchDist = 0; + int farthestMatchID = 0; + for (iMatch = 0; iMatch < MAX_MATCHES; iMatch++) { + if (catalog1.voids[iVoid1].matches[iMatch].dist > farthestMatchDist){ + farthestMatchDist = catalog1.voids[iVoid1].matches[iMatch].dist; + farthestMatchID = iMatch; + } + } + if (rdist < farthestMatchDist) + catalog1.voids[iVoid1].matches[farthestMatchID] = newMatch; } - - } - - MATCHPROPS newMatch; - newMatch.matchID = closestMatchID; - newMatch.commonVol = 1; - newMatch.dist = closestMatchDist; - catalog1.voids[iVoid1].matches.push_back(newMatch); + } } -*/ + printf(" Determining overlap...\n"); for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) { - printf(" Working on void %d of %d...\n", iVoid1, catalog1.numVoids); + printf(" Working on void %d of %d...\n", iVoid1+1, catalog1.numVoids); voidID1 = catalog1.voids[iVoid1].voidID; - for (iVoid2 = 0; iVoid2 < catalog2.numVoids; iVoid2++) { + for (iMatch = 0; iMatch < catalog1.voids[iVoid1].matches.size();iMatch++) { + iVoid2 = catalog1.voids[iVoid1].matches[iMatch].matchID; voidID2 = catalog2.voids[iVoid2].voidID; - match = false; for (iZ1 = 0; iZ1 < catalog1.void2Zones[voidID1].numZones; iZ1++) { zoneID1 = catalog1.void2Zones[voidID1].zoneIDs[iZ1]; @@ -180,6 +188,7 @@ int main(int argc, char **argv) { partID2 = catalog2.zones2Parts[zoneID2].partIDs[p2]; match = false; + if (args.useID_flag) { if (catalog1.part[partID1].uniqueID == catalog2.part[partID2].uniqueID) match = true; @@ -203,35 +212,12 @@ int main(int argc, char **argv) { r2 = pow(3./4./M_PI*catalog2.part[partID2].volume / catalog2.numPartTot, 1./3.); if (rdist <= 0.1*r1 || rdist <= 0.1*r2) match = true; - } + } - if (match) { - bool alreadyMatched = false; - for (iMatch = 0; - iMatch < catalog1.voids[iVoid1].matches.size(); iMatch++) { - if (catalog1.voids[iVoid1].matches[iMatch].matchID == - iVoid2) { - alreadyMatched = true; - break; - } - } - - if (alreadyMatched) { - //catalog1.voids[iVoid1].matches[iMatch].commonVol += 1; - catalog1.voids[iVoid1].matches[iMatch].commonVol += - catalog1.part[partID1].volume; - } else { - MATCHPROPS newMatch; - newMatch.matchID = iVoid2; - //newMatch.commonVol = 1; - newMatch.commonVol = catalog1.part[partID1].volume; - - newMatch.dist = getDist(catalog1, catalog2, iVoid1, iVoid2, - periodicX, periodicY, periodicZ); - catalog1.voids[iVoid1].matches.push_back(newMatch); - } // end if match - } + catalog1.voids[iVoid1].matches[iMatch].commonVol += + catalog2.part[partID2].volume; + } // end if match } // end p2 } // end iZ2 @@ -240,11 +226,13 @@ int main(int argc, char **argv) { } // end iVoid2 } // end match finding + printf(" Sorting matches...\n"); for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) { sortMatches(catalog1.voids[iVoid1].matches); } // count up significant matches + printf(" Categorizing matches...\n"); for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) { closestMatchDist = 0.; for (iMatch = 0; iMatch < catalog1.voids[iVoid1].matches.size(); iMatch++) { @@ -257,12 +245,12 @@ int main(int argc, char **argv) { } // output summary + printf(" Output...\n"); std::string filename; filename = string(args.outfile_arg); - filename = filename.append("_summary.out"); + filename = filename.append("summary.out"); fp = fopen(filename.c_str(), "w"); - //fp = fopen(args.outfile_arg, "w"); - fprintf(fp, "# void ID, radius, radius ratio, common volume ratio, relative dist, num matches, num significant matches\n"); + fprintf(fp, "# void ID, radius, radius ratio, common volume ratio, common volume ratio 2, relative dist, num matches, num significant matches\n"); for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) { int voidID = catalog1.voids[iVoid1].voidID; if (catalog1.voids[iVoid1].numMatches > 0) { @@ -272,27 +260,31 @@ int main(int argc, char **argv) { commonVolRatio = catalog1.voids[iVoid1].matches[0].commonVol / //catalog1.voids[iVoid1].numPart; catalog1.voids[iVoid1].vol; + float volRatio = catalog1.voids[iVoid1].matches[0].commonVol / + catalog2.voids[iVoid2].vol; rdist = catalog1.voids[iVoid1].matches[0].dist; rdist /= catalog1.voids[iVoid1].radius; - fprintf(fp, "%d %.2f %.2f %.2f %.2f %d %d\n", voidID, + fprintf(fp, "%d %.2f %.2f %.2f %.2f %.2f %d %d\n", voidID, catalog1.voids[iVoid1].radiusMpc, rRatio, commonVolRatio, + volRatio, rdist, catalog1.voids[iVoid1].numMatches, catalog1.voids[iVoid1].numBigMatches); } else { - fprintf(fp, "%d %f 0.0 0.0 0.0 0 0\n", voidID, - catalog1.voids[iVoid1].radius); + fprintf(fp, "%d %.2f 0.0 0.0 0.0 0.0 0 0\n", voidID, + catalog1.voids[iVoid1].radiusMpc); } } // end printing fclose(fp); // output detail + printf(" Output detail...\n"); filename = string(args.outfile_arg); - filename = filename.append("_detail.out"); + filename = filename.append("detail.out"); fp = fopen(filename.c_str(), "w"); int MAX_OUT = 10; fprintf(fp, "# void ID, match common vol\n"); @@ -300,14 +292,13 @@ int main(int argc, char **argv) { int voidID = catalog1.voids[iVoid1].voidID; fprintf(fp,"%d: ", voidID); for (iMatch = 0; iMatch < MAX_OUT; iMatch++) { - - if (iMatch < catalog1.voids[iVoid].matches.size()) { + if (iMatch < catalog1.voids[iVoid1].matches.size()) { commonVolRatio = catalog1.voids[iVoid1].matches[iMatch].commonVol / //catalog1.voids[iVoid1].numPart; catalog1.voids[iVoid1].vol; - fprintf(fp, "%.2f ", commonVolRatio); - + fprintf(fp, "%.3f ", catalog1.voids[iVoid1].matches[iMatch].dist); + //fprintf(fp, "%.2f ", commonVolRatio); } else { fprintf(fp, "0.00 "); } @@ -443,10 +434,6 @@ void loadCatalog(const char *partFile, const char *volFile, catalog.voids[i].voidProb = voidProb; catalog.voids[i].radius = pow(voidVol/catalog.numPartTot*3./4./M_PI, 1./3.); - //catalog.voids[i].biggestMatchID = -1; - //catalog.voids[i].biggestMatchVol = 0; - //catalog.voids[i].closestMatchID = -1; - //catalog.voids[i].closestMatchDist = 1.e99; catalog.voids[i].numMatches = 0; catalog.voids[i].numBigMatches = 0; @@ -548,14 +535,15 @@ void sortMatches(std::vector& matches) { MATCHPROPS tempMatch; bool swapped; - if (matches.size() == 1) return; + if (matches.size() <= 1) return; - swapped = false; - while (!swapped) { + swapped = true; + while (swapped) { swapped = false; for (int iMatch = 0; iMatch < matches.size() - 1; iMatch++) { - if (matches[iMatch].commonVol > matches[iMatch+1].commonVol) { - tempMatch = matches[iMatch]; + if (matches[iMatch].dist > matches[iMatch+1].dist) { + //if (matches[iMatch].commonVol < matches[iMatch+1].commonVol) { + tempMatch = matches[iMatch+1]; matches[iMatch+1] = matches[iMatch]; matches[iMatch] = tempMatch; swapped = true; From d4febc694247c10940187d35f9352e4b468857e1 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 20 Dec 2012 13:38:59 -0600 Subject: [PATCH 02/16] streamlined merger tree processing --- crossCompare/analysis/datasetsToAnalyze.py | 6 +++--- crossCompare/analysis/mergerTree.py | 18 +++++++++--------- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/crossCompare/analysis/datasetsToAnalyze.py b/crossCompare/analysis/datasetsToAnalyze.py index 82cf9ac..88452f6 100755 --- a/crossCompare/analysis/datasetsToAnalyze.py +++ b/crossCompare/analysis/datasetsToAnalyze.py @@ -2,14 +2,14 @@ workDir = "/home/psutter2/workspace/Voids/" -figDir = "./figs" -dataDir = "./data" +dataDir = "/home/psutter2/workspace/Voids/multidark/crossCompare/mergerTree/" CTOOLS_PATH = "../../c_tools/" +baseSampleDir = "multidark/md_ss0.000175/sample_md_ss0.000175_z0.56_d00/" + sampleDirList = [ - "multidark/md_ss0.000175/sample_md_ss0.000175_z0.56_d00/", "multidark/md_ss0.0004/sample_md_ss0.0004_z0.56_d00/", ] diff --git a/crossCompare/analysis/mergerTree.py b/crossCompare/analysis/mergerTree.py index 47a4179..89c5e13 100755 --- a/crossCompare/analysis/mergerTree.py +++ b/crossCompare/analysis/mergerTree.py @@ -30,30 +30,30 @@ if not os.access(filename, os.F_OK): parms = imp.load_source("name", filename) globals().update(vars(parms)) -if not os.access(figDir, os.F_OK): - os.makedirs(figDir) +if not os.access(dataDir, os.F_OK): + os.makedirs(dataDir) outFileName = dataDir + "/" + dataNameBase #+ ".dat" +with open(workDir+baseSampleDir+"/sample_info.dat", 'rb') as input: + baseSample = pickle.load(input) + for (iSample, sampleDir) in enumerate(sampleDirList): - if iSample == 0: continue with open(workDir+sampleDir+"/sample_info.dat", 'rb') as input: sample = pickle.load(input) - with open(workDir+sampleDirList[0]+"/sample_info.dat", 'rb') as input: - baseSample = pickle.load(input) - - print " Working with", sample.fullName, "..." + print " Working with", sample.fullName, "...", sampleName = sample.fullName binPath = CTOOLS_PATH+"/analysis/voidOverlap" logFile = os.getcwd()+"/mergerTree.log" - stepOutputFileName = outFileName + "_" + sampleName + "_" + stepOutputFileName = outFileName + "_" + baseSample.fullName + "_" + \ + sampleName + "_" #stepOutputFileName = os.getcwd()+"/data/overlap_" - launchVoidOverlap(baseSample, sample, workDir+sampleDirList[0], + launchVoidOverlap(baseSample, sample, workDir+baseSampleDir, workDir+sampleDir, binPath, thisDataPortion="central", logFile=logFile, continueRun=False, workDir=workDir, From a7cdb807047b7ca217bdaaa04057c0f2cf31eef1 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Sun, 23 Dec 2012 16:38:54 -0600 Subject: [PATCH 03/16] added observation flag in info file; stack voids now shifts particles based on redshift of subsample --- c_tools/libzobov/particleInfo.cpp | 10 ++++++++++ c_tools/mock/generateFromCatalog.cpp | 1 + c_tools/mock/generateMock.cpp | 1 + 3 files changed, 12 insertions(+) diff --git a/c_tools/libzobov/particleInfo.cpp b/c_tools/libzobov/particleInfo.cpp index 9711db6..fa6f467 100644 --- a/c_tools/libzobov/particleInfo.cpp +++ b/c_tools/libzobov/particleInfo.cpp @@ -10,6 +10,7 @@ bool loadParticleInfo(ParticleInfo& info, const std::string& extra_info) { int numpart; + int isObservation; NcFile f_info(extra_info.c_str()); @@ -23,6 +24,7 @@ bool loadParticleInfo(ParticleInfo& info, info.ranges[2][0] = f_info.get_att("range_z_min")->as_double(0); info.ranges[2][1] = f_info.get_att("range_z_max")->as_double(0); info.mask_index = f_info.get_att("mask_index")->as_int(0); //PMS + isObservation = f_info.get_att("is_observation")->as_int(0); //PMS for (int i = 0; i < 3; i++) info.length[i] = info.ranges[i][1] - info.ranges[i][0]; @@ -59,6 +61,14 @@ bool loadParticleInfo(ParticleInfo& info, for (int i = 0; i < numpart; i++) info.particles[i].z = mul*f.readReal32(); f.endCheckpoint(); + + if (!isObservation) { + for (int i = 0; i < numpart; i++) { + info.particles[i].x += info.ranges[0][0]; + info.particles[i].y += info.ranges[1][0]; + info.particles[i].z += info.ranges[2][0]; + } + } } catch (const NoSuchFileException& e) { diff --git a/c_tools/mock/generateFromCatalog.cpp b/c_tools/mock/generateFromCatalog.cpp index d5a0aca..54cbf87 100644 --- a/c_tools/mock/generateFromCatalog.cpp +++ b/c_tools/mock/generateFromCatalog.cpp @@ -466,6 +466,7 @@ void saveForZobov(ParticleData& pdata, const string& fname, const string& paramn fp.add_att("range_z_min", -Lmax/100.); fp.add_att("range_z_max", Lmax/100.); fp.add_att("mask_index", pdata.mask_index); // PMS + fp.add_att("is_observation", 1); // PMS int nOutputPart = pdata.mask_index; //int nOutputPart = pdata.pos.size(); diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index a7830dd..fe21928 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -338,6 +338,7 @@ void saveBox(SimuData *&boxed, const std::string& outbox) f.add_att("range_z_min", ranges[4]); f.add_att("range_z_max", ranges[5]); f.add_att("mask_index", -1); + f.add_att("is_observation", 0); NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart); NcDim *NumSnap_dim = f.add_dim("numsnap_dim", num_snapshots); From 481da6edd0dd0cc03276eb1821416257626794c3 Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Wed, 16 Jan 2013 14:06:36 -0600 Subject: [PATCH 04/16] fixed bug in barycenter calculation --- c_tools/stacking/pruneVoids.cpp | 75 +++++++++++++++++++++++++++------ 1 file changed, 61 insertions(+), 14 deletions(-) diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index cd88462..b0c284c 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -103,6 +103,9 @@ int main(int argc, char **argv) { args_info.periodic_arg); // check for periodic box + periodicX = 0; + periodicY = 0; + periodicZ = 0; if (!args_info.isObservation_flag) { if ( strchr(args_info.periodic_arg, 'x') != NULL) { periodicX = 1; @@ -143,9 +146,9 @@ int main(int argc, char **argv) { temp = (float *) malloc(numPartTot * sizeof(float)); volNorm = numPartTot/(boxLen[0]*boxLen[1]*boxLen[2]); - printf("VOL NORM = %f\n", volNorm); + printf(" VOL NORM = %f\n", volNorm); - printf("CENTRAL DEN = %f\n", args_info.maxCentralDen_arg); + printf(" CENTRAL DEN = %f\n", args_info.maxCentralDen_arg); fread(&dummy, 1, 4, fp); fread(temp, numPartTot, 4, fp); @@ -311,13 +314,16 @@ int main(int argc, char **argv) { voids[iVoid].barycenter[2] = 0.; for (p = 0; p < voids[iVoid].numPart; p++) { - dist[0] = fabs(voidPart[p].x - voids[iVoid].center[0]); - dist[1] = fabs(voidPart[p].y - voids[iVoid].center[1]); - dist[2] = fabs(voidPart[p].z - voids[iVoid].center[2]); + dist[0] = voidPart[p].x - voids[iVoid].center[0]; + dist[1] = voidPart[p].y - voids[iVoid].center[1]; + dist[2] = voidPart[p].z - voids[iVoid].center[2]; - if (periodicX) dist[0] = fmin(dist[0], boxLen[0]-dist[0]); - if (periodicY) dist[1] = fmin(dist[1], boxLen[1]-dist[1]); - if (periodicZ) dist[2] = fmin(dist[2], boxLen[2]-dist[2]); + if (periodicX && fabs(dist[0]) > boxLen[0]/2.) + dist[0] = dist[0] - copysign(boxLen[0], dist[0]); + if (periodicY && fabs(dist[1]) > boxLen[1]/2.) + dist[1] = dist[1] - copysign(boxLen[1], dist[1]); + if (periodicZ && fabs(dist[2]) > boxLen[2]/2.) + dist[2] = dist[2] - copysign(boxLen[2], dist[2]); voids[iVoid].barycenter[0] += voidPart[p].vol*(dist[0]); voids[iVoid].barycenter[1] += voidPart[p].vol*(dist[1]); @@ -331,12 +337,25 @@ int main(int argc, char **argv) { voids[iVoid].barycenter[1] += voids[iVoid].center[1]; voids[iVoid].barycenter[2] += voids[iVoid].center[2]; - if (periodicX) - voids[iVoid].barycenter[0] = fmod(voids[iVoid].barycenter[0], boxLen[0]); - if (periodicY) - voids[iVoid].barycenter[1] = fmod(voids[iVoid].barycenter[1], boxLen[1]); - if (periodicZ) - voids[iVoid].barycenter[2] = fmod(voids[iVoid].barycenter[2], boxLen[2]); + if (periodicX) { + if (voids[iVoid].barycenter[0] > boxLen[0]) + voids[iVoid].barycenter[0] = voids[iVoid].barycenter[0] - boxLen[0]; + if (voids[iVoid].barycenter[0] < 0) + voids[iVoid].barycenter[0] = boxLen[0] - voids[iVoid].barycenter[0]; + } + if (periodicY) { + if (voids[iVoid].barycenter[1] > boxLen[1]) + voids[iVoid].barycenter[1] = voids[iVoid].barycenter[1] - boxLen[1]; + if (voids[iVoid].barycenter[1] < 1) + voids[iVoid].barycenter[1] = boxLen[1] - voids[iVoid].barycenter[1]; + } + if (periodicZ) { + if (voids[iVoid].barycenter[2] > boxLen[2]) + voids[iVoid].barycenter[2] = voids[iVoid].barycenter[2] - boxLen[2]; + if (voids[iVoid].barycenter[2] < 2) + voids[iVoid].barycenter[2] = boxLen[2] - voids[iVoid].barycenter[2]; + } + // compute central density centralRad = voids[iVoid].radius/args_info.centralRadFrac_arg; @@ -450,6 +469,11 @@ int main(int argc, char **argv) { voids[iVoid].nearestEdge = nearestEdge; } // iVoid + int numWrong = 0; + int numHighDen = 0; + int numEdge = 0; + int numTooSmall = 0; + printf(" Picking winners and losers...\n"); for (iVoid = 0; iVoid < numVoids; iVoid++) { voids[iVoid].accepted = 1; @@ -462,36 +486,43 @@ int main(int argc, char **argv) { if (voids[iVoid].centralDen > args_info.maxCentralDen_arg) { voids[iVoid].accepted = -1; + numHighDen++; } // toss out voids that are obviously wrong if (voids[iVoid].densCon > 1.e3) { voids[iVoid].accepted = -4; + numWrong++; } if (strcmp(args_info.dataPortion_arg, "edge") == 0 && tolerance*voids[iVoid].maxRadius < voids[iVoid].nearestMock) { voids[iVoid].accepted = -3; + numEdge++; } if (strcmp(args_info.dataPortion_arg, "central") == 0 && tolerance*voids[iVoid].maxRadius > voids[iVoid].nearestMock) { voids[iVoid].accepted = -3; + numEdge++; } if (voids[iVoid].radius < args_info.rMin_arg) { voids[iVoid].accepted = -2; + numTooSmall++; } // *always* clean out near edges since there are no mocks there if (tolerance*voids[iVoid].maxRadius > voids[iVoid].nearestEdge) { voids[iVoid].accepted = -3; + numEdge++; } // assume the lower z-boundary is "soft" in observations if (args_info.isObservation_flag && voids[iVoid].redshift < args_info.zMin_arg) { voids[iVoid].accepted = -3; + numEdge++; } } @@ -501,6 +532,10 @@ int main(int argc, char **argv) { } printf(" Number kept: %d (out of %d)\n", numKept, numVoids); + printf(" Rejected %d near the edge\n", numEdge); + printf(" Rejected %d too small\n", numTooSmall); + printf(" Rejected %d obviously bad\n", numWrong); + printf(" Rejected %d too high central density\n", numHighDen); printf(" Output...\n"); fp = fopen(args_info.output_arg, "w"); @@ -531,6 +566,10 @@ int main(int argc, char **argv) { fprintf(fpBarycenter, "%d %e %e %e\n", voids[iVoid].voidID, +// TEST + //voids[iVoid].center[0], + //voids[iVoid].center[1], + //voids[iVoid].center[2]); voids[iVoid].barycenter[0], voids[iVoid].barycenter[1], voids[iVoid].barycenter[2]); @@ -540,11 +579,19 @@ int main(int argc, char **argv) { voids[iVoid].nearestMock); double outCenter[3]; +// TEST + //outCenter[0] = voids[iVoid].center[0]; + //outCenter[1] = voids[iVoid].center[1]; + //outCenter[2] = voids[iVoid].center[2]; outCenter[0] = voids[iVoid].barycenter[0]; outCenter[1] = voids[iVoid].barycenter[1]; outCenter[2] = voids[iVoid].barycenter[2]; if (args_info.isObservation_flag) { +// TEST + //outCenter[0] = (voids[iVoid].center[0]-boxLen[0]/2.)*100.; + //outCenter[1] = (voids[iVoid].center[1]-boxLen[1]/2.)*100.; + //outCenter[2] = (voids[iVoid].center[2]-boxLen[2]/2.)*100.; outCenter[0] = (voids[iVoid].barycenter[0]-boxLen[0]/2.)*100.; outCenter[1] = (voids[iVoid].barycenter[1]-boxLen[1]/2.)*100.; outCenter[2] = (voids[iVoid].barycenter[2]-boxLen[2]/2.)*100.; From 4338a08f447210de97ecdcbf178758744b8a7ea1 Mon Sep 17 00:00:00 2001 From: nhamaus Date: Thu, 17 Jan 2013 15:10:25 +0100 Subject: [PATCH 05/16] fixed particle selection problem in generateMock --- c_tools/mock/generateMock.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index ab51315..ad334e9 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -221,7 +221,12 @@ void generateOutput(SimuData *data, int axis, // cleared. That way new particles can be appended if this is a multi-file snapshot. void selectBox(SimuData *simu, std::vector& targets, generateMock_info& args_info) { - float subsample = args_info.subsample_given ? args_info.subsample_arg : 1.0; + float subsample; + if (args_info.subsample_given) { + subsample = args_info.subsample_arg; + } else { + subsample = 1.0; + } double ranges[3][2] = { { args_info.rangeX_min_arg, args_info.rangeX_max_arg }, { args_info.rangeY_min_arg, args_info.rangeY_max_arg }, From ea58ff2b837d0346bad0cae3e863fadec841f64c Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Thu, 17 Jan 2013 09:31:13 -0600 Subject: [PATCH 06/16] made default profile bin size 'auto' --- .../pipeline_source/prepareCatalogs.in.py | 24 ++++++++++++------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/python_tools/pipeline_source/prepareCatalogs.in.py b/python_tools/pipeline_source/prepareCatalogs.in.py index fc7c927..cb0cc63 100755 --- a/python_tools/pipeline_source/prepareCatalogs.in.py +++ b/python_tools/pipeline_source/prepareCatalogs.in.py @@ -79,7 +79,7 @@ startCatalogStage = 1 endCatalogStage = 4 startAPStage = 1 -endAPStage = 7 +endAPStage = 1 ZOBOV_PATH = "@CMAKE_BINARY_DIR@/zobov/" CTOOLS_PATH = "@CMAKE_BINARY_DIR@/c_tools/" @@ -124,7 +124,7 @@ newSample = Sample(dataFile = "{dataFile}", zBoundaryMpc = ({zMinMpc}, {zMaxMpc}), omegaM = {omegaM}, minVoidRadius = {minRadius}, - profileBinSize = 1.0, + profileBinSize = "auto", includeInHubble = True, partOfCombo = False, isCombo = False, @@ -135,12 +135,20 @@ newSample = Sample(dataFile = "{dataFile}", useLightCone = {useLightCone}, subsample = {subsample}) dataSampleList.append(newSample) -newSample.addStack({zMin}, {zMax}, 2*{minRadius} , 2*{minRadius}+2, True, False) -newSample.addStack({zMin}, {zMax}, 2*{minRadius} , 2*{minRadius}+4, True, False) -newSample.addStack({zMin}, {zMax}, 2*{minRadius}+2, 2*{minRadius}+6, True, False) -newSample.addStack({zMin}, {zMax}, 2*{minRadius}+6, 2*{minRadius}+10, True, False) -newSample.addStack({zMin}, {zMax}, 2*{minRadius}+10, 2*{minRadius}+18, True, False) -newSample.addStack({zMin}, {zMax}, 2*{minRadius}+18, 2*{minRadius}+24, True, False) +newSample.addStack(0.0, 5.0, 20, 25, False, False) +newSample.addStack(0.0, 5.0, 30, 35, False, False) +newSample.addStack(0.0, 5.0, 40, 45, False, False) +newSample.addStack(0.0, 5.0, 50, 55, False, False) +newSample.addStack(0.0, 5.0, 60, 65, False, False) +newSample.addStack(0.0, 5.0, 70, 75, False, False) +newSample.addStack(0.0, 5.0, 80, 85, False, False) +newSample.addStack(0.0, 5.0, 90, 95, False, False) +#newSample.addStack({zMin}, {zMax}, 2*{minRadius} , 2*{minRadius}+2, True, False) +#newSample.addStack({zMin}, {zMax}, 2*{minRadius} , 2*{minRadius}+4, True, False) +#newSample.addStack({zMin}, {zMax}, 2*{minRadius}+2, 2*{minRadius}+6, True, False) +#newSample.addStack({zMin}, {zMax}, 2*{minRadius}+6, 2*{minRadius}+10, True, False) +#newSample.addStack({zMin}, {zMax}, 2*{minRadius}+10, 2*{minRadius}+18, True, False) +#newSample.addStack({zMin}, {zMax}, 2*{minRadius}+18, 2*{minRadius}+24, True, False) """ for (iFile, redshift) in enumerate(redshifts): fileNum = fileNums[iFile] From ae9701e99c8e762c142323fd68aabcdafadf551a Mon Sep 17 00:00:00 2001 From: "P.M. Sutter" Date: Fri, 18 Jan 2013 09:56:31 -0600 Subject: [PATCH 07/16] replaced void shape fitter with rudimentary moment of inertia calculation --- c_tools/stacking/pruneVoids.cpp | 2 +- .../apTools/profiles/__init__.py | 1 + .../void_python_tools/backend/launchers.py | 48 +++++++++++-------- 3 files changed, 30 insertions(+), 21 deletions(-) diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index b0c284c..37f39a3 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -490,7 +490,7 @@ int main(int argc, char **argv) { } // toss out voids that are obviously wrong - if (voids[iVoid].densCon > 1.e3) { + if (voids[iVoid].densCon > 1.e4) { voids[iVoid].accepted = -4; numWrong++; } diff --git a/python_tools/void_python_tools/apTools/profiles/__init__.py b/python_tools/void_python_tools/apTools/profiles/__init__.py index b05e408..0f04a3e 100644 --- a/python_tools/void_python_tools/apTools/profiles/__init__.py +++ b/python_tools/void_python_tools/apTools/profiles/__init__.py @@ -1,6 +1,7 @@ from build import * from draw import * from fit import * +from inertia import * from mcmc import * from generateExpFigure import * from getSurveyProps import * diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 45c117a..28cd70c 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -314,7 +314,8 @@ def launchPrune(sample, binPath, thisDataPortion=None, def launchVoidOverlap(sample1, sample2, sample1Dir, sample2Dir, binPath, thisDataPortion=None, logFile=None, workDir=None, - continueRun=None, outputFile=None): + continueRun=None, outputFile=None, + matchMethod=None): sampleName1 = sample1.fullName sampleName2 = sample2.fullName @@ -358,7 +359,7 @@ def launchVoidOverlap(sample1, sample2, sample1Dir, sample2Dir, cmd += " --zonePartFile2=" + sample2Dir+"/voidPart_" + \ str(sampleName2)+".dat" - cmd += " --useID" + if matchMethod == "useID": cmd += " --useID" cmd += periodicLine cmd += " --outfile=" + outputFile cmd += " &> " + logFile @@ -600,13 +601,18 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None, os.system("mv %s %s" % ("tree.data", treeFile)) os.system("mv %s %s" % ("void_indexes.txt", voidDir+"/")) - os.system("mv %s %s" % ("posx.nc", voidDir+"/posx.nc")) - os.system("mv %s %s" % ("posy.nc", voidDir+"/posy.nc")) - os.system("mv %s %s" % ("posz.nc", voidDir+"/posz.nc")) - os.system("mv %s %s" % ("redshifts.nc", voidDir+"/redshifts.nc")) + os.system("mv %s %s" % ("posx.nc", voidDir+"/")) + os.system("mv %s %s" % ("posy.nc", voidDir+"/")) + os.system("mv %s %s" % ("posz.nc", voidDir+"/")) + os.system("mv %s %s" % ("z_void_indexes.txt", voidDir+"/")) + os.system("mv %s %s" % ("z_posx.nc", voidDir+"/")) + os.system("mv %s %s" % ("z_posy.nc", voidDir+"/")) + os.system("mv %s %s" % ("z_posz.nc", voidDir+"/")) + os.system("mv %s %s" % ("redshifts.nc", voidDir+"/")) os.system("mv %s %s" % ("indexes.nc", voidDir+"/")) os.system("mv %s %s" % ("kdtree_stackvoids.dat", voidDir+"/")) os.system("mv %s %s" % ("centers.txt", voidDir+"/")) + os.system("mv %s %s" % ("z_centers.txt", voidDir+"/")) os.system("mv %s %s" % ("sky_positions.txt", voidDir+"/")) os.system("mv %s %s" % ("check.txt", voidDir+"/")) os.system("mv %s %s" % ("tracer.txt", voidDir+"/")) @@ -899,22 +905,24 @@ def launchFit(sample, stack, logFile=None, voidDir=None, figDir=None, while badChain: Rexpect = (stack.rMin+stack.rMax)/2 Rtruncate = stack.rMin*3. + 1 # TEST - if sample.dataType == "observation": - ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect, - Niter=300000, - Nburn=100000, - Rextracut=Rtruncate) - else: - ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect, - Niter=300000, - Nburn=100000, - Rextracut=Rtruncate) - badChain = (args[0][0] > 0.5 or args[0][1] > stack.rMax or \ - args[0][2] > stack.rMax) and \ - (ntries < maxtries) + #if sample.dataType == "observation": + # ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect, + # Niter=300000, + # Nburn=100000, + # Rextracut=Rtruncate) + #else: + # ret,fits,args = vp.fit_ellipticity(voidDir,Rbase=Rexpect, + # Niter=300000, + # Nburn=100000, + # Rextracut=Rtruncate) + #badChain = (args[0][0] > 0.5 or args[0][1] > stack.rMax or \ + # args[0][2] > stack.rMax) and \ + # (ntries < maxtries) + ret,fits,args = vp.compute_inertia(voidDir, stack.rMax) + badChain = False ntries += 1 - np.save(voidDir+"/chain.npy", ret) + #np.save(voidDir+"/chain.npy", ret) np.savetxt(voidDir+"/fits.out", fits) plotTitle = "Sample: "+sample.nickName+\ From b410694f6ab40b45f933d24536e9545d0bac8176 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sat, 19 Jan 2013 12:14:08 -0500 Subject: [PATCH 08/16] Added an assertion check to trouble shoot potential SEGV --- c_tools/libzobov/loadZobov.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/c_tools/libzobov/loadZobov.cpp b/c_tools/libzobov/loadZobov.cpp index 7399874..2d5a2ff 100644 --- a/c_tools/libzobov/loadZobov.cpp +++ b/c_tools/libzobov/loadZobov.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -152,6 +153,7 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName, for (int j = 0; j < z.allVoids[volId].zId.size(); j++) { int zzid = z.allVoids[volId].zId[j]; + assert(zzid < z.allZones.size()); actualNumber += z.allZones[zzid].pId.size(); } From da9f7ff4addce27ffb4a0d6d53c45187a6bec43d Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sat, 19 Jan 2013 12:15:10 -0500 Subject: [PATCH 09/16] Reposition the particles according to ranges and not length/2 (as designed at the beginning) From 61bb0684ccf67bfca384af41b6f20b9dbb97de1f Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sat, 19 Jan 2013 12:16:04 -0500 Subject: [PATCH 10/16] Adjusted showSliceExact for new evolution of the file formats. Do proper translation/rescaling/rotation if requested From 9ce6b3f683f7047a2c01383d1fdd7def89be2251 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Sat, 19 Jan 2013 12:16:52 -0500 Subject: [PATCH 11/16] Recount the number of adjacencies, which is not fully preserved by the merging and can cause assertion failure in other software --- zobov/voztie.c | 33 +++++++++++++++++++++++++++++---- 1 file changed, 29 insertions(+), 4 deletions(-) diff --git a/zobov/voztie.c b/zobov/voztie.c index 77530cb..cc8c006 100644 --- a/zobov/voztie.c +++ b/zobov/voztie.c @@ -8,6 +8,7 @@ int main(int argc, char *argv[]) { FILE *part, *adj, *vol; char partfile[200], *suffix, adjfile[200], volfile[200], *outDir; float *vols, volstemp; + int *cnt_adj; PARTADJ *adjs; @@ -78,13 +79,17 @@ int main(int argc, char *argv[]) { if (mockIndex == -1) mockIndex = np; // END PMS + cnt_adj = (int *)malloc(np*sizeof(int)); + if (cnt_adj == NULL) + printf("Could not allocate cnt_adj.\n"); + adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ)); if (adjs == NULL) printf("Couldn't allocate adjs.\n"); vols = (float *)malloc(np*sizeof(float)); if (vols == NULL) printf("Couldn't allocate vols.\n"); orig = (int *)malloc(nvpmax*sizeof(int)); if (orig == NULL) printf("Couldn't allocate orig.\n"); - if ((vols == NULL) || (orig == NULL) || (adjs == NULL)) { + if ((cnt_adj==NULL) || (vols == NULL) || (orig == NULL) || (adjs == NULL)) { printf("Not enough memory to allocate. Exiting.\n"); exit(0); } @@ -238,11 +243,24 @@ printf("\n"); // END OMS /* Adjacencies: first the numbers of adjacencies, and the number we're actually going to write per particle */ + + // Recount the number of adjacencies after merge + for(i=0;i i) nout++; fwrite(&nout,1,sizeof(int),adj); - for (j=0;j i) + for (j=0;j i) fwrite(&(adjs[i].adj[j]),1,sizeof(int),adj); + } } fclose(adj); @@ -275,5 +295,10 @@ printf("\n"); fclose(vol); + free(vols); + free(cnt_adj); + free(adjs); + free(orig); + return(0); } From 64e226c78683074538d3bd2b1183886199302c0d Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 21 Jan 2013 21:38:38 -0500 Subject: [PATCH 12/16] Be compatible with older version of the parameter file (without is_observation) --- c_tools/libzobov/particleInfo.cpp | 46 +++++++++++++++++++++++++------ 1 file changed, 38 insertions(+), 8 deletions(-) diff --git a/c_tools/libzobov/particleInfo.cpp b/c_tools/libzobov/particleInfo.cpp index fa6f467..d3eb2e8 100644 --- a/c_tools/libzobov/particleInfo.cpp +++ b/c_tools/libzobov/particleInfo.cpp @@ -1,3 +1,4 @@ +#include #include #include #include "particleInfo.hpp" @@ -5,6 +6,34 @@ using namespace std; using namespace CosmoTool; +template +double retrieve_attr_safe_double(NcFile& f, const char *name, double defval) +{ + NcAtt *a = f.get_att(name); + if (a == 0) + { + if (failure) + abort(); + return defval; + } + return a->as_double(0); +} + +template +int retrieve_attr_safe_int(NcFile& f, const char *name, int defval) +{ + NcAtt *a = f.get_att(name); + if (a == 0) + { + if (failure) + abort(); + return defval; + } + return a->as_int(0); +} + + + bool loadParticleInfo(ParticleInfo& info, const std::string& particles, const std::string& extra_info) @@ -13,18 +42,19 @@ bool loadParticleInfo(ParticleInfo& info, int isObservation; NcFile f_info(extra_info.c_str()); + NcError nerr(NcError::verbose_nonfatal); if (!f_info.is_valid()) return false; - info.ranges[0][0] = f_info.get_att("range_x_min")->as_double(0); - info.ranges[0][1] = f_info.get_att("range_x_max")->as_double(0); - info.ranges[1][0] = f_info.get_att("range_y_min")->as_double(0); - info.ranges[1][1] = f_info.get_att("range_y_max")->as_double(0); - info.ranges[2][0] = f_info.get_att("range_z_min")->as_double(0); - info.ranges[2][1] = f_info.get_att("range_z_max")->as_double(0); - info.mask_index = f_info.get_att("mask_index")->as_int(0); //PMS - isObservation = f_info.get_att("is_observation")->as_int(0); //PMS + info.ranges[0][0] = retrieve_attr_safe_double(f_info, "range_x_min", 0); + info.ranges[0][1] = retrieve_attr_safe_double(f_info, "range_x_max", 0); + info.ranges[1][0] = retrieve_attr_safe_double(f_info, "range_y_min", 0); + info.ranges[1][1] = retrieve_attr_safe_double(f_info, "range_y_max", 0); + info.ranges[2][0] = retrieve_attr_safe_double(f_info, "range_z_min", 0); + info.ranges[2][1] = retrieve_attr_safe_double(f_info, "range_z_max", 0); + info.mask_index = retrieve_attr_safe_int(f_info, "mask_index", 0); + isObservation = retrieve_attr_safe_int(f_info, "is_observation", 0); for (int i = 0; i < 3; i++) info.length[i] = info.ranges[i][1] - info.ranges[i][0]; From 9bd0cbbb8820cd076cb75f8e0eab83788a4f76de Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 29 Jan 2013 13:58:12 -0600 Subject: [PATCH 13/16] Fixed regeneration step in generateMock. Fixed gadget loading (missing initialization for multiple files) --- c_tools/mock/generateMock.cpp | 8 +++++--- c_tools/mock/loaders/gadget_loader.cpp | 2 +- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index ad334e9..b3864c6 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -317,9 +317,9 @@ void buildBox(SimuData *simu, long num_targets, long loaded, for (int j = 0; j < 3; j++) { - boxed->Pos[j][loaded] = (simu->Pos[j][pid]-ranges[j*2])*mul[j]; - assert(boxed->Pos[j][loaded] > 0); - assert(boxed->Pos[j][loaded] < 1); + boxed->Pos[j][loaded] = max(min((simu->Pos[j][pid]-ranges[j*2])*mul[j], double(1)), double(0)); + assert(boxed->Pos[j][loaded] >= 0); + assert(boxed->Pos[j][loaded] <= 1); } uniqueID[loaded] = (simu_uniqueID != 0) ? simu_uniqueID[pid] : 0; expansion_fac[loaded] = efac[pid]; @@ -400,6 +400,7 @@ void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& a mul = new double[3]; ranges = new double[6]; snapshot_split = new long[*num_snapshots]; + expansion_fac = new double[boxed->NumPart]; boxed->new_attribute("uniqueID", uniqueID, delete_adaptor); @@ -408,6 +409,7 @@ void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& a boxed->new_attribute("particle_id", particle_id, delete_adaptor); boxed->new_attribute("num_snapshots", num_snapshots, delete_adaptor); boxed->new_attribute("snapshot_split", snapshot_split, delete_adaptor); + boxed->new_attribute("expansion_fac", expansion_fac, delete_adaptor); v_id->get(particle_id, boxed->NumPart); v_snap->get(snapshot_split, *num_snapshots); diff --git a/c_tools/mock/loaders/gadget_loader.cpp b/c_tools/mock/loaders/gadget_loader.cpp index 1c108ee..906f6c1 100644 --- a/c_tools/mock/loaders/gadget_loader.cpp +++ b/c_tools/mock/loaders/gadget_loader.cpp @@ -77,7 +77,7 @@ public: SimulationLoader *gadgetLoader(const std::string& snapshot, double Mpc_unitLength, int flags) { - bool singleFile; + bool singleFile = false; int num_files; SimuData *d; From 5fb3adfc89116c0bf6188e35ba04c39fe5c548cf Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 29 Jan 2013 13:59:25 -0600 Subject: [PATCH 14/16] Fixed potential memory corruption while loading particle information. Slightly reworked the source code --- c_tools/stacking/pruneVoids.cpp | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index 37f39a3..1e23830 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -87,7 +87,7 @@ int main(int argc, char **argv) { double nearestEdge, redshift; char line[500], junkStr[10]; int mask_index; - double ranges[2][3], boxLen[3], mul; + double ranges[3][2], boxLen[3], mul; double volNorm, radius; int clock1, clock2; int periodicX=0, periodicY=0, periodicZ=0; @@ -174,7 +174,7 @@ int main(int argc, char **argv) { part[p].y += ranges[1][0]; part[p].z += ranges[2][0]; } - } + } fclose(fp); printf(" Read %d particles...\n", numPartTot); @@ -341,19 +341,19 @@ int main(int argc, char **argv) { if (voids[iVoid].barycenter[0] > boxLen[0]) voids[iVoid].barycenter[0] = voids[iVoid].barycenter[0] - boxLen[0]; if (voids[iVoid].barycenter[0] < 0) - voids[iVoid].barycenter[0] = boxLen[0] - voids[iVoid].barycenter[0]; + voids[iVoid].barycenter[0] = boxLen[0] + voids[iVoid].barycenter[0]; } if (periodicY) { if (voids[iVoid].barycenter[1] > boxLen[1]) voids[iVoid].barycenter[1] = voids[iVoid].barycenter[1] - boxLen[1]; - if (voids[iVoid].barycenter[1] < 1) - voids[iVoid].barycenter[1] = boxLen[1] - voids[iVoid].barycenter[1]; + if (voids[iVoid].barycenter[1] < 0) + voids[iVoid].barycenter[1] = boxLen[1] + voids[iVoid].barycenter[1]; } if (periodicZ) { if (voids[iVoid].barycenter[2] > boxLen[2]) voids[iVoid].barycenter[2] = voids[iVoid].barycenter[2] - boxLen[2]; - if (voids[iVoid].barycenter[2] < 2) - voids[iVoid].barycenter[2] = boxLen[2] - voids[iVoid].barycenter[2]; + if (voids[iVoid].barycenter[2] < 0) + voids[iVoid].barycenter[2] = boxLen[2] + voids[iVoid].barycenter[2]; } @@ -397,8 +397,8 @@ int main(int argc, char **argv) { for (p = 0; p < voids[iVoid].numPart; p++) { dist[0] = fabs(voidPart[p].x - voids[iVoid].barycenter[0]); - dist[0] = fabs(voidPart[p].y - voids[iVoid].barycenter[1]); - dist[0] = fabs(voidPart[p].z - voids[iVoid].barycenter[2]); + dist[1] = fabs(voidPart[p].y - voids[iVoid].barycenter[1]); + dist[2] = fabs(voidPart[p].z - voids[iVoid].barycenter[2]); if (periodicX) dist[0] = fmin(dist[0], boxLen[0]-dist[0]); if (periodicY) dist[1] = fmin(dist[1], boxLen[1]-dist[1]); @@ -564,12 +564,12 @@ int main(int argc, char **argv) { voids[iVoid].densCon, voids[iVoid].voidProb); - fprintf(fpBarycenter, "%d %e %e %e\n", + fprintf(fpBarycenter, "%d %e %e %e\n", voids[iVoid].voidID, // TEST - //voids[iVoid].center[0], - //voids[iVoid].center[1], - //voids[iVoid].center[2]); +// voids[iVoid].center[0], +// voids[iVoid].center[1], +//` voids[iVoid].center[2], voids[iVoid].barycenter[0], voids[iVoid].barycenter[1], voids[iVoid].barycenter[2]); From f209b89565e0863ad672614be1d476fd4be7ea7a Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 29 Jan 2013 14:00:23 -0600 Subject: [PATCH 15/16] Reactivated the density threshold test From 116cd950e3a48e63518d31ead588147e45a20d9a Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 29 Jan 2013 14:02:00 -0600 Subject: [PATCH 16/16] Added support for regeneration (to finely test impact of peculiar velocities) --- pipeline/generateCatalog.py | 3 +- .../pipeline_source/prepareCatalogs.in.py | 1 + .../void_python_tools/backend/launchers.py | 44 ++++++++++++++----- 3 files changed, 35 insertions(+), 13 deletions(-) diff --git a/pipeline/generateCatalog.py b/pipeline/generateCatalog.py index 5941933..fdf3729 100755 --- a/pipeline/generateCatalog.py +++ b/pipeline/generateCatalog.py @@ -20,6 +20,7 @@ if (len(sys.argv) > 1): print " Cannot find parameter file %s!" % filename exit(-1) parms = imp.load_source("name", filename) + regenerateFlag = False globals().update(vars(parms)) else: print " Using default parameters" @@ -90,7 +91,7 @@ for sample in dataSampleList: launchGenerate(sample, GENERATE_PATH, workDir=workDir, inputDataDir=inputDataDir, zobovDir=zobovDir, figDir=figDir, logFile=logFile, useLCDM=useLCDM, - continueRun=continueRun) + continueRun=continueRun, regenerate=regenerateFlag) # -------------------------------------------------------------------------- if (startCatalogStage <= 2) and (endCatalogStage >= 2) and not sample.isCombo: diff --git a/python_tools/pipeline_source/prepareCatalogs.in.py b/python_tools/pipeline_source/prepareCatalogs.in.py index cb0cc63..0ea0bd1 100755 --- a/python_tools/pipeline_source/prepareCatalogs.in.py +++ b/python_tools/pipeline_source/prepareCatalogs.in.py @@ -81,6 +81,7 @@ endCatalogStage = 4 startAPStage = 1 endAPStage = 1 +regenerateFlag = False ZOBOV_PATH = "@CMAKE_BINARY_DIR@/zobov/" CTOOLS_PATH = "@CMAKE_BINARY_DIR@/c_tools/" freshStack = True diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index 28cd70c..d7e01bf 100755 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -24,11 +24,18 @@ ncFloat = 'f8' # Double precision # ----------------------------------------------------------------------------- def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, zobovDir=None, figDir=None, logFile=None, useLCDM=False, - continueRun=None): + continueRun=None,regenerate=False): if sample.dataType == "observation": sampleName = sample.fullName + if regenerate: + inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+sampleName+".par" + outputFile = zobovDir+"/regenerated_zobov_slice_" + sampleName + else: + inputParameterFlag = "" + outputFile = zobovDir+"/zobov_slice_" + sampleName + if sample.dataFile == "": datafile = inputDataDir+"/"+sampleName else: @@ -50,16 +57,17 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, zMax %g density_fake %g %s - """ % (datafile, maskFile, zobovDir+"/zobov_slice_"+sampleName, + %s + """ % (datafile, maskFile, outputFile, zobovDir+"/zobov_slice_"+sampleName+".par", sample.zBoundary[0], sample.zBoundary[1], sample.fakeDensity, - useLCDMFlag) + useLCDMFlag, inputParameterFlag) parmFile = os.getcwd()+"/generate_"+sample.fullName+".par" file(parmFile, mode="w").write(conf) - if not (continueRun and jobSuccessful(logFile, "Done!\n")): + if regenerate or not (continueRun and jobSuccessful(logFile, "Done!\n")): cmd = "%s --configFile=%s &> %s" % (binPath,parmFile,logFile) os.system(cmd) if jobSuccessful(logFile, "Done!\n"): @@ -97,6 +105,13 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, datafile = inputDataDir+"/"+sample.dataFile + if regenerate: + inputParameterFlag = "inputParameter " + zobovDir+"/zobov_slice_"+sampleName+".par" + outputFile = zobovDir+"/regenerated_zobov_slice_" + sampleName + else: + inputParameterFlag = "" + outputFile = zobovDir+"/zobov_slice_" + sampleName + if sample.usePecVel: includePecVelString = "peculiarVelocities" else: @@ -134,20 +149,21 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None, rangeZ_min %g rangeZ_max %g subsample %g - """ % (dataFileLine, zobovDir+"/zobov_slice_"+sampleName, + %s + """ % (dataFileLine, outputFile, zobovDir+"/zobov_slice_"+sampleName+".par", includePecVelString, useLightConeString, sample.dataUnit, xMin, xMax, yMin, yMax, sample.zBoundaryMpc[0], sample.zBoundaryMpc[1], - sample.subsample) + sample.subsample,inputParameterFlag) parmFile = os.getcwd()+"/generate_"+sample.fullName+".par" file(parmFile, mode="w").write(conf) - if not (continueRun and jobSuccessful(logFile, "Done!\n")): + if regenerate or not (continueRun and jobSuccessful(logFile, "Done!\n")): cmd = "%s --configFile=%s &> %s" % (binPath,parmFile,logFile) os.system(cmd) if jobSuccessful(logFile, "Done!\n"): @@ -258,6 +274,7 @@ def launchPrune(sample, binPath, thisDataPortion=None, sample.boxLen <= 1.e-1: periodicLine += "z" periodicLine += "' " + periodicLine = "" if not (continueRun and jobSuccessful(logFile, "NetCDF: Not a valid ID\n")): cmd = binPath @@ -297,6 +314,9 @@ def launchPrune(sample, binPath, thisDataPortion=None, str(thisDataPortion)+"_"+\ str(sampleName)+".out" cmd += " &> " + logFile + f=file("run_prune.sh",mode="w") + f.write(cmd) + f.close() os.system(cmd) if jobSuccessful(logFile, "NetCDF: Not a valid ID\n") or \ @@ -380,14 +400,14 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None, voidDir=None, freshStack=True, runSuffix=None, zobovDir=None, INCOHERENT=False, ranSeed=None, summaryFile=None, - continueRun=None, dataType=None): + continueRun=None, dataType=None, prefixRun=""): sampleName = sample.fullName runSuffix = getStackSuffix(stack.zMin, stack.zMax, stack.rMin, stack.rMax, thisDataPortion) - logFile = logDir+"/stack_"+sampleName+"_"+runSuffix+".out" + logFile = logDir+("/%sstack_"%prefixRun)+sampleName+"_"+runSuffix+".out" treeFile = voidDir+"/tree.data" @@ -414,7 +434,7 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None, maxDen = 0.2*float(maskIndex)/float(totalPart) else: maskIndex = 999999999 - maxDen = 0.2 + maxDen = -0.8 if INCOHERENT: nullTestFlag = "INCOHERENT" @@ -455,7 +475,7 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None, zobovDir+"/vol_"+sampleName+".dat", stack.rMin, stack.rMax, - zobovDir+"/zobov_slice_"+sampleName, + zobovDir+("/%szobov_slice_"%prefixRun)+sampleName, zobovDir+"/zobov_slice_"+sampleName+".par", maxDen, centralRadius, @@ -471,7 +491,7 @@ def launchStack(sample, stack, binPath, thisDataPortion=None, logDir=None, zobovDir+"/boundaryDistances_"+thisDataPortion+"_"+sampleName+".out", rescaleFlag) - parmFile = os.getcwd()+"/stack_"+sample.fullName+".par" + parmFile = os.getcwd()+("/%sstack_"%prefixRun)+sample.fullName+".par" fp = file(parmFile, mode="w") fp.write(conf)