diff --git a/.gitignore b/.gitignore index 487c453..c0ba1bd 100644 --- a/.gitignore +++ b/.gitignore @@ -1,8 +1,4 @@ *~ *.o *.prog -voz1b1 -jozov -vozinit -voztie .mydepend diff --git a/c_tools/analysis/voidOverlap.cpp b/c_tools/analysis/voidOverlap.cpp index ee46a4b..01c61fb 100644 --- a/c_tools/analysis/voidOverlap.cpp +++ b/c_tools/analysis/voidOverlap.cpp @@ -177,26 +177,72 @@ int main(int argc, char **argv) { newMatch.commonVolProg = 0; newMatch.dist = rdist; - if (rdist/catalog1.voids[iVoid1].radius > 1.5) continue; + //if (rdist > 1.5) continue; + if (rdist/catalog1.voids[iVoid1].radius > 2.0) continue; - 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; + // see if center is contained in void + match = false; + voidID1 = catalog1.voids[iVoid1].voidID; + for (iZ1 = 0; iZ1 < catalog1.void2Zones[voidID1].numZones; iZ1++) { + zoneID1 = catalog1.void2Zones[voidID1].zoneIDs[iZ1]; + for (p1 = 0; p1 < catalog1.zones2Parts[zoneID1].numPart; p1++) { + partID1 = catalog1.zones2Parts[zoneID1].partIDs[p1]; + + dist[0] = fabs(catalog1.part[partID1].x - + catalog2.voids[iVoid2].barycenter[0]); + dist[1] = fabs(catalog1.part[partID1].y - + catalog2.voids[iVoid2].barycenter[1]); + dist[2] = fabs(catalog1.part[partID1].z - + catalog2.voids[iVoid2].barycenter[2]); + + if (periodicX) dist[0] = fmin(dist[0], 1.0 - dist[0]); + if (periodicY) dist[1] = fmin(dist[1], 1.0 - dist[1]); + if (periodicZ) dist[2] = fmin(dist[2], 1.0 - dist[2]); + + rdist = sqrt(dist[0]*dist[0] + dist[1]*dist[1] + dist[2]*dist[2]); + + r1 = pow(3./4./M_PI*catalog1.part[partID1].volume / + catalog1.numPartTot, 1./3.); + if (rdist <= matchDist*4*r1) { + match = true; + break; + } + } } + + if (!match) continue; + //if (rdist/catalog1.voids[iVoid1].radius > 1.) continue; + + //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; + //} } } + // pick the closest matches to speed up computation + for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) { + int numPotential = catalog1.voids[iVoid1].matches.size(); + catalog1.voids[iVoid1].numMatches = numPotential; + + if (numPotential > MAX_MATCHES) { + sortMatches(catalog1.voids[iVoid1].matches, catalog2, 1); + } + + catalog1.voids[iVoid1].matches.resize(MAX_MATCHES); + } + printf(" Determining overlap...\n"); for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) { printf(" Working on void %d of %d...\n", iVoid1+1, catalog1.numVoids); @@ -267,9 +313,13 @@ int main(int argc, char **argv) { int matchID = catalog1.voids[iVoid1].matches[iMatch].matchID; catalog1.voids[iVoid1].matches[iMatch].merit = catalog1.voids[iVoid1].matches[iMatch].commonVolOrig * - catalog1.voids[iVoid1].matches[iMatch].commonVolProg / + catalog1.voids[iVoid1].matches[iMatch].commonVolOrig / catalog1.voids[iVoid1].vol / - catalog2.voids[matchID].vol; + catalog2.voids[iVoid1].vol; + //catalog1.voids[iVoid1].matches[iMatch].commonVolOrig * + //catalog1.voids[iVoid1].matches[iMatch].commonVolProg / + //catalog1.voids[iVoid1].vol / + //catalog2.voids[matchID].vol; //pow(catalog1.voids[iVoid1].matches[iMatch].commonVol,2) / //catalog1.voids[iVoid1].vol / //catalog2.voids[matchID].vol; @@ -296,7 +346,6 @@ int main(int argc, char **argv) { catalog1.voids[iVoid1].vol; if (commonVolRatio > 0.2) catalog1.voids[iVoid1].numBigMatches++; } - catalog1.voids[iVoid1].numMatches = catalog1.voids[iVoid1].matches.size(); } // output summary @@ -305,7 +354,7 @@ int main(int argc, char **argv) { filename = string(args.outfile_arg); filename = filename.append("summary.out"); fp = fopen(filename.c_str(), "w"); - fprintf(fp, "# void ID, radius, radius ratio, common volume ratio (to original), common volume ratio (to match), relative dist, num matches, num significant matches, match ID, merit, ellipticity ratio\n"); + fprintf(fp, "# void ID, radius, radius ratio, common volume ratio (to original), common volume ratio (to match), relative dist, num matches, num significant matches, match ID, merit, ellipticity ratio, density contrast\n"); for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) { int voidID = catalog1.voids[iVoid1].voidID; if (catalog1.voids[iVoid1].numMatches > 0) { @@ -322,7 +371,7 @@ int main(int argc, char **argv) { rdist = catalog1.voids[iVoid1].matches[0].dist; rdist /= catalog1.voids[iVoid1].radius; - fprintf(fp, "%d %.4f %.4f %e %e %.4f %d %d %d %e %e\n", voidID, + fprintf(fp, "%d %.4f %.4f %e %e %.4f %d %d %d %e %e %e\n", voidID, //fprintf(fp, "%d %.4f %.4f %.4f %.4f %.4f %d %d %d\n", voidID, catalog1.voids[iVoid1].radiusMpc, rRatio, @@ -333,10 +382,11 @@ int main(int argc, char **argv) { catalog1.voids[iVoid1].numBigMatches, catalog2.voids[iVoid2].voidID, catalog1.voids[iVoid1].matches[0].merit, - ellipRatio); + ellipRatio, + catalog1.voids[iVoid1].densCon); } else { - fprintf(fp, "%d %.4f 0.0 0.0 0.0 0.0 0.0 0 0 0 0.0\n", voidID, + fprintf(fp, "%d %.4f 0.0 0.0 0.0 0.0 0.0 0 0 0 0.0 0.0\n", voidID, catalog1.voids[iVoid1].radiusMpc); } } // end printing @@ -347,25 +397,26 @@ int main(int argc, char **argv) { filename = string(args.outfile_arg); filename = filename.append("detail.out"); fp = fopen(filename.c_str(), "w"); - int MAX_OUT = 10; - fprintf(fp, "# void ID, match common vol\n"); + fprintf(fp, "# match ID, merit, relative dist, relative radius\n"); for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) { int voidID = catalog1.voids[iVoid1].voidID; - fprintf(fp,"%d: ", voidID); - for (iMatch = 0; iMatch < MAX_OUT; iMatch++) { - if (iMatch < catalog1.voids[iVoid1].matches.size()) { - commonVolRatio = catalog1.voids[iVoid1].matches[iMatch].commonVolOrig / - //catalog1.voids[iVoid1].numPart; - catalog1.voids[iVoid1].vol; + fprintf(fp,"#%d (%.4f Mpc/h):\n", voidID, catalog1.voids[iVoid1].radiusMpc); + for (iMatch = 0; iMatch < catalog1.voids[iVoid1].matches.size(); iMatch++) { + commonVolRatio = catalog1.voids[iVoid1].matches[iMatch].commonVolOrig / + //catalog1.voids[iVoid1].numPart; + catalog1.voids[iVoid1].vol; - fprintf(fp, "%.3f %.2f ", catalog1.voids[iVoid1].matches[iMatch].dist, commonVolRatio); - //fprintf(fp, "%.2f ", commonVolRatio); - } else { - fprintf(fp, "0.00 "); - } + int matchID = catalog1.voids[iVoid1].matches[iMatch].matchID; + fprintf(fp, "%d %e %.4f %.4f\n", matchID, + catalog1.voids[iVoid1].matches[iMatch].merit, + catalog1.voids[iVoid1].matches[iMatch].dist/ + catalog1.voids[iVoid1].radius, + catalog2.voids[matchID].radius/ + catalog1.voids[iVoid1].radius); } - fprintf(fp, "\n"); + fprintf(fp, "\n "); + fprintf(fp, "\n "); } // end printing detail fclose(fp); diff --git a/c_tools/stacking/pruneVoids.cpp b/c_tools/stacking/pruneVoids.cpp index 57d5e26..49c2ec1 100644 --- a/c_tools/stacking/pruneVoids.cpp +++ b/c_tools/stacking/pruneVoids.cpp @@ -578,6 +578,8 @@ int main(int argc, char **argv) { // compute eigenvalues and vectors for orientation and shape double inertia[9]; + for (int i = 0; i < 9; i++) inertia[i] = 0.; + for (int p = 0; p < voids[iVoid].numPart; p++) { dist[0] = voidPart[p].x - voids[iVoid].barycenter[0]; dist[1] = voidPart[p].y - voids[iVoid].barycenter[1]; @@ -614,9 +616,11 @@ int main(int argc, char **argv) { float c = sqrt(2.5*(gsl_vector_get(voids[iVoid].eval,0) + gsl_vector_get(voids[iVoid].eval,1) - gsl_vector_get(voids[iVoid].eval,2))); - float ca = c/a; + float ca; float cb = c/b; - voids[iVoid].ellip = 1.0 - c/a; + if (a < c) ca = a/c; + if (a >= c) ca = c/a; + voids[iVoid].ellip = fabs(1.0 - ca); } // iVoid @@ -938,7 +942,7 @@ void outputVoids(string outputDir, string sampleName, string prefix, outVoid.radius, outVoid.voidID); - fprintf(fpShapes, "%d %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f %.2f\n", + fprintf(fpShapes, "%d %.6f %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e %.2e\n", outVoid.voidID, outVoid.ellip, gsl_vector_get(outVoid.eval, 0), diff --git a/c_tools/zobov2/findrtop.c b/c_tools/zobov2/jozov2/findrtop.c similarity index 100% rename from c_tools/zobov2/findrtop.c rename to c_tools/zobov2/jozov2/findrtop.c diff --git a/c_tools/zobov2/jozov2.cpp b/c_tools/zobov2/jozov2/jozov2.cpp similarity index 100% rename from c_tools/zobov2/jozov2.cpp rename to c_tools/zobov2/jozov2/jozov2.cpp diff --git a/c_tools/zobov2/jozov2.hpp b/c_tools/zobov2/jozov2/jozov2.hpp similarity index 100% rename from c_tools/zobov2/jozov2.hpp rename to c_tools/zobov2/jozov2/jozov2.hpp diff --git a/c_tools/zobov2/jozov2_io.cpp b/c_tools/zobov2/jozov2/jozov2_io.cpp similarity index 100% rename from c_tools/zobov2/jozov2_io.cpp rename to c_tools/zobov2/jozov2/jozov2_io.cpp diff --git a/c_tools/zobov2/jozov2_watershed.cpp b/c_tools/zobov2/jozov2/jozov2_watershed.cpp similarity index 100% rename from c_tools/zobov2/jozov2_watershed.cpp rename to c_tools/zobov2/jozov2/jozov2_watershed.cpp diff --git a/c_tools/zobov2/jozov2_zones.cpp b/c_tools/zobov2/jozov2/jozov2_zones.cpp similarity index 100% rename from c_tools/zobov2/jozov2_zones.cpp rename to c_tools/zobov2/jozov2/jozov2_zones.cpp diff --git a/c_tools/zobov2/zobov.hpp b/c_tools/zobov2/jozov2/zobov.hpp similarity index 100% rename from c_tools/zobov2/zobov.hpp rename to c_tools/zobov2/jozov2/zobov.hpp diff --git a/c_tools/zobov2/voz1b1/voz.h b/c_tools/zobov2/voz1b1/voz.h new file mode 100644 index 0000000..1d48e92 --- /dev/null +++ b/c_tools/zobov2/voz1b1/voz.h @@ -0,0 +1,29 @@ +#ifndef __VOZ_H +#define __VOZ_H + +#define MAXVERVER 100000 +#define NGUARD 84 /*Actually, the number of SPACES between guard points +##define NGUARD 42 /*Actually, the number of SPACES between guard points + in each dim */ + +typedef int pid_t; + +typedef struct Partadj { + pid_t nadj; + pid_t *adj; +} PARTADJ; + +#ifdef __cplusplus +extern "C" { +#endif + +int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs); +int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol); +int posread(char *posfile, float ***p, float fact); + +#ifdef __cplusplus +} +#endif + + +#endif diff --git a/c_tools/zobov2/voz1b1/voz1b1.cpp b/c_tools/zobov2/voz1b1/voz1b1.cpp new file mode 100644 index 0000000..0bad2ce --- /dev/null +++ b/c_tools/zobov2/voz1b1/voz1b1.cpp @@ -0,0 +1,395 @@ +#include +#include +#include +#include +#include +#include "libqhull/qhull_a.h" +#include "voz.h" +#include "voz_io.hpp" + +using CosmoTool::MINIARG_STRING; +using CosmoTool::MINIARG_DOUBLE; +using CosmoTool::MINIARG_INT; +using CosmoTool::MINIARG_NULL; +using boost::format; +using namespace std; + +#define DL for (d=0;d<3;d++) + +bool checkParameters(int *numdiv, int *b) +{ + for (int i = 0; i < 3; i++) + { + if (numdiv[i] < 0 || b[i] < 0 || b[i] >= numdiv[i]) + return false; + } + return true; +} + + +struct BoxData +{ + float width[3], totwidth[3]; + float width2[3], totwidth2[3]; + float bf, s[3], g[3], c[3]; + coordT *parts; + pid_t *orig; + pid_t nvpall, nvp, nvpbuf; + bool guard_added; + double xyz_min[3], xyz_max[3]; + + void prepareBox(const PositionData& pdata, int *b_id, int *numdivs); + + void checkParticle(double *xyz, bool& in_main, bool& in_buf); + void addGuardPoints(); + void findBoundary(); +}; + +void BoxData::checkParticle(double *xyz, bool& in_main, bool& in_buf) +{ + in_main = in_buf = true; + for (int d = 0; d < 3; d++) + { + xyz[d] -= (double)c[d]; + if (xyz[d] > width2[d]) + xyz[d] -= width[d]; + if (xyz[d] < -width2[d]) + xyz[d] += width[d]; + + in_buf = in_buf && (fabs(xyz[d]) < totwidth2[d]); + in_main = in_main && (fabs(xyz[d]) <= width2[d]); + } +} + +void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) +{ + guard_added = false; + + for (int i = 0; i < 3; i++) + { + width[i] = (pdata.xyz_max[i] - pdata.xyz_min[i])/numdivs[i]; + width2[i] = 0.5*width[i]; + + totwidth[i] = width[i]+2.*bf; + totwidth2[i] = width2[i] + bf; + + s[i] = width[i]/(float)NGUARD; + if ((bf*bf - 2.*s[i]*s[i]) < 0.) + { + printf("bf = %f, s = %f.\n",bf,s[i]); + printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.", + sqrt(2.)*width[i]/bf); + exit(0); + } + g[i] = (bf / 2.)*(1. + sqrt(1 - 2.*s[i]*s[i]/(bf*bf))); + cout << format("s[%d] = %f, bf = %f, g[%d] = %f.") % s[i] % i % bf % g[i] % i << endl; + + c[i] = b_id[i] * width[i]; + } + cout.flush(); + + cout << format("c: %f,%f,%f") % c[0] % c[1] % c[2] << endl; + + /* Assign temporary array*/ + nvpbuf = 0; /* Number of particles to tesselate, including buffer */ + nvp = 0; /* Without the buffer */ + + for (pid_t i = 0; i < pdata.np; i++) + { + double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] }; + bool is_it_in_buf, is_it_in_main; + + checkParticle(xyz, is_it_in_buf, is_it_in_main); + + if (is_it_in_buf) + nvpbuf++; + if (is_it_in_main) + nvp++; + } + + nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard points */ + + parts = new coordT[3*nvpbuf]; + orig = new pid_t[nvpbuf]; + + if (parts == 0) + { + cout << "Unable to allocate parts" << endl; + exit(0); + } + if (orig == 0) + { + cout << "Unable to allocate orig" << endl; + exit(0); + } + + nvp = 0; nvpall = 0; /* nvp = number of particles without buffer */ + + for (int j = 0; j < 3; j++) + { + xyz_min[j] = -std::numeric_limits::max(); + xyz_max[j] = std::numeric_limits::max(); + } + for (pid_t i = 0; i < pdata.np; i++) + { + bool is_it_in_main, is_it_in_buf; + double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] }; + + checkParticle(xyz, is_it_in_main, is_it_in_buf); + + if (is_it_in_main) { + for (int j = 0; j < 3; j++) + { + parts[3*nvp+j] = xyz[j]; + xyz_min[j] = min(xyz_min[j], xyz[j]); + xyz_max[j] = max(xyz_max[j], xyz[j]); + } + orig[nvp] = i; + nvp++; + } + } + cout << format("nvp = %d") %nvp << endl; + cout << format("x: %f,%f; y: %f,%f; z:%f,%f") % xyz_min[0] % xyz_max[0] % xyz_min[1] % xyz_max[1] % xyz_min[2] % xyz_max[2] << endl; + nvpbuf = nvp; + for (pid_t i = 0; i < pdata.np; i++) + { + bool is_it_in_main, is_it_in_buf; + double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] }; + + checkParticle(xyz, is_it_in_main, is_it_in_buf); + + if (is_it_in_buf && !is_it_in_main) + { + for (int j = 0; j < 3; j++) + { + parts[3*nvpbuf+j] = xyz[j]; + xyz_min[j] = min(xyz_min[j], xyz[j]); + xyz_max[j] = max(xyz_max[j], xyz[j]); + } + orig[nvpbuf] = i; + + nvpbuf++; + } + } + cout << format("nvpbuf = %d") % nvpbuf << endl; + cout << format("x: %f,%f; y: %f,%f; z:%f,%f\n") % xyz_min[0] % xyz_max[0] % xyz_min[1] % xyz_max[1] % xyz_min[2] % xyz_max[2] << endl; + nvpall = nvpbuf; + + double predict = 1; + for (int j = 0; j < 3; j++) + predict *= (width[j]+2*bf); + predict *= pdata.np; + cout << format("There should be ~ %g points; there are %d\n") % predict % nvpbuf << endl; +} + +void BoxData::addGuardPoints() +{ + int rot_g[3][3] = { {0, 0, 1}, {0,1,0}, {1,0,0} }; + int rot_si[3][3] = { {1, 0, 0}, {1,0,0}, {0,1,0} }; + int rot_sj[3][3] = { {0, 1, 0}, {0,0,1}, {0,0,1} }; + + for (int k = 0; k < 3; k++) + { + + /* Add guard points */ + for (int i=0; i <= NGUARD; i++) + { + for (int j=0; j <= NGUARD; j++) + { + /* Bottom */ + for (int a = 0; a < 3; a++) + parts[3*nvpall+a] = -width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] - rot_g[k][a] * g[a]; + + nvpall++; + + /* Top */ + for (int a = 0; a < 3; a++) + parts[3*nvpall+a] = (2*rot_g[k][a]-1)*width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] + rot_g[k][a] * g[a]; + + nvpall++; + } + } + } +} + +void BoxData::findBoundary() +{ + double BF = std::numeric_limits::max(); + + for (int j = 0; j < 3; j++) + { + xyz_min[j] = -BF; + xyz_max[j] = BF; + } + + for (pid_t i = nvpbuf; i < nvpall; i++) { + for (int j = 0; j < 3; j++) + { + xyz_min[j] = std::min(xyz_min[j], parts[3*i+j]); + xyz_max[j] = std::min(xyz_max[j], parts[3*i+j]); + } + } +} + +void saveTesselation(const string& outfile, PositionData& pdata, BoxData& boxdata, PARTADJ *adjs, float *vols) +{ + ofstream out(outfile.c_str()); + if (!out) + { + cout << format("Unable to open %s") % outfile << endl; + exit(0); + } + out.write((char *)&pdata.np, sizeof(pid_t)); + out.write((char *)&boxdata.nvp, sizeof(pid_t)); + cout << format("nvp = %d") % boxdata.nvp << endl; + + /* Tell us where the original particles were */ + out.write((char *)boxdata.orig, sizeof(pid_t)*boxdata.nvp); + /* Volumes*/ + out.write((char *)vols,sizeof(float)*boxdata.nvp); + /* Adjacencies */ + for (pid_t i = 0; i < boxdata.nvp; i++) + { + out.write((char*)&adjs[i].nadj, sizeof(pid_t)); + if (adjs[i].nadj > 0) + out.write((char *)adjs[i].adj, adjs[i].nadj*sizeof(pid_t)); + else + (cout << "0").flush(); + } + out.close(); +} + +int main(int argc, char *argv[]) { + PositionData pdata; + BoxData boxdata; + + coordT deladjs[3*MAXVERVER], points[3*MAXVERVER]; + pointT intpoints[3*MAXVERVER]; + string outfile; + ofstream out; + + char *suffix, *outDir, *posfile; + PARTADJ *adjs; + float *vols; + pid_t *orig; + double border, boxsize[3]; + int b[3], numdiv[3]; + double totalvol; + + CosmoTool::MiniArgDesc args[] = { + { "POSITION FILE", &posfile, MINIARG_STRING }, + { "BORDER SIZE", &border, MINIARG_DOUBLE }, + { "BOX_X", &boxsize[0], MINIARG_DOUBLE }, + { "BOX_Y", &boxsize[1], MINIARG_DOUBLE }, + { "BOX_Z", &boxsize[2], MINIARG_DOUBLE }, + { "SUFFIX", &suffix, MINIARG_STRING }, + { "NUM_DIVISION_X", &numdiv[0], MINIARG_INT }, + { "NUM_DIVISION_Y", &numdiv[1], MINIARG_INT }, + { "NUM_DIVISION_Z", &numdiv[2], MINIARG_INT }, + { "B0", &b[0], MINIARG_INT }, + { "B1", &b[1], MINIARG_INT }, + { "B2", &b[2], MINIARG_INT }, + { "OUTPUT DIRECTORY", &outDir, MINIARG_STRING }, + { 0, 0, MINIARG_NULL } + }; + + if (!CosmoTool::parseMiniArgs(argc, argv, args)) + return 1; + + if (!checkParameters(numdiv, b)) + return 2; + + /* Boxsize should be the range in r, yielding a range 0-1 */ + if (!pdata.readFrom(posfile)) + return 3; + + (cout << pdata.np << " particles" << endl).flush(); + + pdata.findExtrema(); + + (cout << boost::format("np: %d, x: %f,%f; y: %f,%f; z: %f,%f") + % pdata.np + % pdata.xyz_min[0] % pdata.xyz_max[0] + % pdata.xyz_min[1] % pdata.xyz_max[1] + % pdata.xyz_min[2] % pdata.xyz_max[2]).flush(); + + + if (border > 0.) + boxdata.bf = border; + else + boxdata.bf = 0.1; + + boxdata.prepareBox(pdata, b, numdiv); + pdata.destroy(); + boxdata.addGuardPoints(); + + adjs = new PARTADJ[boxdata.nvpall]; + if (adjs == 0) + { + cout << "Unable to allocate adjs" << endl; + return 0; + } + + boxdata.findBoundary(); + + cout << format("Added guard points to total %d points (should be %d)") + % boxdata.nvpall % (boxdata.nvpbuf + 6*(NGUARD+1)*(NGUARD+1)) << endl; + cout << format("x: %f,%f; y: %f,%f; z:%f,%f") % boxdata.xyz_min[0] % boxdata.xyz_max[0] % boxdata.xyz_min[1] % boxdata.xyz_max[1] % boxdata.xyz_min[2] % boxdata.xyz_max[2] << endl; + + /* Do tesselation*/ + printf("File read. Tessellating ...\n"); fflush(stdout); + int exitcode = delaunadj(boxdata.parts, boxdata.nvp, boxdata.nvpbuf, boxdata.nvpall, &adjs); + if (exitcode != 0) + { + printf("Error while tesselating. Stopping here."); fflush(stdout); + return 4; + } + + /* Now calculate volumes*/ + printf("Now finding volumes ...\n"); fflush(stdout); + vols = new float[boxdata.nvp]; + + for (pid_t i = 0; i < boxdata.nvp; i++) + { /* Just the original particles + * Assign adjacency coordinate array*/ + /* Volumes */ + for (int j = 0; j < adjs[i].nadj; j++) + { + for (int d = 0; d < 3; d++) + { + deladjs[3*j + d] = boxdata.parts[3*adjs[i].adj[j]+d] - boxdata.parts[3*i+d]; + + if (deladjs[3*j+d] < -boxdata.width2[d]) + deladjs[3*j+d] += boxdata.width[d]; + if (deladjs[3*j+d] > boxdata.width2[d]) + deladjs[3*j+d] -= boxdata.width[d]; + } + } + + exitcode = vorvol(deladjs, points, intpoints, adjs[i].nadj, &(vols[i])); + vols[i] *= pdata.np/pdata.V0; + if ((i % 1000) == 0) + cout << format("%d: %d, %f") % i % adjs[i].nadj % vols[i] << endl; + } + + /* Get the adjacencies back to their original values */ + for (pid_t i=0; i +#include +#include +#include "voz_io.hpp" +#include "CosmoTool/fortran.hpp" +#include + +using namespace CosmoTool; +using namespace std; + +bool PositionData::readFrom(const string& fname) +{ + try + { + UnformattedRead f(fname.c_str()); + + f.beginCheckpoint(); + np = f.readInt32(); + f.endCheckpoint(); + + for (int j = 0; j < 3; j++) + { + xyz[j] = new float[np]; + + f.beginCheckpoint(); + for (int p = 0; p < np; p++) + xyz[j][p] = f.readReal32(); + f.endCheckpoint(); + } + } + catch (const NoSuchFileException& e) + { + return false; + } + catch (const InvalidUnformattedAccess& e) + { + return false; + } + catch (const EndOfFileException& e) + { + return false; + } + + return true; +} + +void PositionData::findExtrema() +{ + const float BF = std::numeric_limits::max(); + + for (int i = 0; i < 3; i++) + { + xyz_min[i] = BF; + xyz_max[i] = -BF; + } + + V0=1; + for (int i = 0; i < 3; i++) + { + for (pid_t p = 0; p < np; p++) + { + xyz_min[i] = min(xyz_min[i], xyz[i][p]); + xyz_max[i] = max(xyz_max[i], xyz[i][p]); + } + V0 *= (xyz_max[i]-xyz_min[i]); + } +} diff --git a/c_tools/zobov2/voz1b1/voz_io.hpp b/c_tools/zobov2/voz1b1/voz_io.hpp new file mode 100644 index 0000000..77a1e5c --- /dev/null +++ b/c_tools/zobov2/voz1b1/voz_io.hpp @@ -0,0 +1,24 @@ +#ifndef __VOZ_IO_HPP +#define __VOZ_IO_HPP + +#include + +struct PositionData +{ + float *xyz[3]; + pid_t np; + float xyz_min[3], xyz_max[3]; + float V0; + + void destroy() + { + for (int j = 0; j < 3; j++) + delete[] xyz[j]; + } + + void findExtrema(); + + bool readFrom(const std::string& fname); +}; + +#endif diff --git a/c_tools/zobov2/voz1b1/vozutil.c b/c_tools/zobov2/voz1b1/vozutil.c new file mode 100644 index 0000000..c56da6f --- /dev/null +++ b/c_tools/zobov2/voz1b1/vozutil.c @@ -0,0 +1,245 @@ +#include "libqhull/qhull_a.h" +#include "voz.h" + +#define FOREACHvertex2_(vertices) FOREACHsetelement_(vertexT, vertices2,vertex2) + +/*char qh_version[] = "user_eg 3.1 2001/10/04"; */ + +int compar(const void * n1, const void * n2) { + int i1,i2; + + i1 = *(int *)n1; + i2 = *(int *)n2; + return 2*(i1 > i2) - 1 + (i1 == i2); +} + +/* Finds Delaunay adjacencies of a set of points */ +int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs) { + + int dim= 3; /* dimension of points */ + boolT ismalloc= False; /* True if qhull should free points in qh_freeqhull() or reallocation */ + char flags[250]; /* option flags for qhull, see qh_opt.htm */ + FILE *outfile= stdout; /* output from qh_produce_output() + use NULL to skip qh_produce_output() */ + FILE *errfile= stderr; /* error messages from qhull code */ + int exitcode; /* 0 if no error from qhull */ + int curlong, totlong; /* memory remaining after qh_memfreeshort */ + int i, ver, count; + int numfacets, numsimplicial, numridges, totneighbors, numneighbors, + numcoplanars, numtricoplanars; + setT *vertices, *vertices2, *vertex_points, *coplanar_points; + vertexT *vertex, **vertexp; + vertexT *vertex2, **vertex2p; + int vertex_i, vertex_n; + facetT *facet, **facetp, *neighbor, **neighborp; + pointT *point, **pointp; + int numdiv; + + PARTADJ adjst; + + int errorReported = 0; + + adjst.adj = (int *)malloc(MAXVERVER*sizeof(int)); + if (adjst.adj == NULL) { + printf("Unable to allocate adjst.adj\n"); + exit(0); + } + + /* Delaunay triangulation*/ + sprintf (flags, "qhull s d Qt"); + exitcode= qh_new_qhull (dim, nvpall, points, ismalloc, + flags, outfile, errfile); + + + if (!exitcode) { /* if no error */ + /* 'qh facet_list' contains the convex hull */ + + /* From qh_printvneighbors */ + qh_countfacets(qh facet_list, NULL, 0, &numfacets, &numsimplicial, + &totneighbors, &numridges, &numcoplanars, &numtricoplanars); + qh_vertexneighbors(); + vertices= qh_facetvertices (qh facet_list, NULL, 0); + vertex_points= qh_settemp (nvpall); + coplanar_points= qh_settemp (nvpall); + qh_setzero (vertex_points, 0, nvpall); + qh_setzero (coplanar_points, 0, nvpall); + FOREACHvertex_(vertices) + qh_point_add (vertex_points, vertex->point, vertex); + FORALLfacet_(qh facet_list) { + FOREACHpoint_(facet->coplanarset) + qh_point_add (coplanar_points, point, facet); + } + ver = 0; + FOREACHvertex_i_(vertex_points) { + (*adjs)[ver].nadj = 0; + if (vertex) { + /* Count the neighboring vertices, check that all are real + neighbors */ + adjst.nadj = 0; + FOREACHneighbor_(vertex) { + if ((*adjs)[ver].nadj > -1) { + if (neighbor->visitid) { + vertices2 = neighbor->vertices; + FOREACHvertex2_(vertices2) { + if (ver != qh_pointid(vertex2->point)) { + adjst.adj[adjst.nadj] = (int)qh_pointid(vertex2->point); + adjst.nadj ++; + if (adjst.nadj > MAXVERVER) { + printf("Increase MAXVERVER to at least %d!\n",adjst.nadj); + exit(0); + } + } + } + } else { + printf(" %d",ver); + (*adjs)[ver].nadj = -1; /* There are unreal vertices here */ + } + } + } + } else (*adjs)[ver].nadj = -2; + + /* Enumerate the unique adjacencies*/ + if (adjst.nadj >= 4) { + qsort((void *)adjst.adj, adjst.nadj, sizeof(int), &compar); + count = 1; + + for (i=1; i= nvpbuf && !errorReported) { + errorReported = 1; + printf("Guard point encountered. Increase border and/or nguard.\n"); + printf("P:(%f,%f,%f), G: (%f,%f,%f)\n",points[3*ver],points[3*ver+1],points[3*ver+2], + points[3*adjst.adj[i]],points[3*adjst.adj[i]+1],points[3*adjst.adj[i]+2]); + } + count++; + } + (*adjs)[ver].adj = (int *)malloc(count*sizeof(int)); + if ((*adjs)[ver].adj == NULL) { + printf("Unable to allocate (*adjs)[ver].adj\n"); + exit(0); + } + (*adjs)[ver].adj[0] = adjst.adj[0]; + count = 1; + for (i=1; i %d\n",adjst.nadj,ver,ver); + exit(0); + } + ver++; + if (ver == nvp) break; + } + qh_settempfree (&coplanar_points); + qh_settempfree (&vertex_points); + qh_settempfree (&vertices); + } + qh_freeqhull(!qh_ALL); /* free long memory */ + qh_memfreeshort (&curlong, &totlong); /* free short memory and memory allocator */ + if (curlong || totlong) + fprintf (errfile, "qhull internal warning (delaunadj): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong); + free(adjst.adj); + return exitcode; +} + +/* Calculates the Voronoi volume from a set of Delaunay adjacencies */ +int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol) { + int dim= 3; /* dimension of points */ + boolT ismalloc= False; /* True if qhull should free points in qh_freeqhull() or reallocation */ + char flags[250]; /* option flags for qhull, see qh_opt.htm */ + FILE *outfile= NULL; /* output from qh_produce_output() + use NULL to skip qh_produce_output() */ + FILE *errfile= stderr; /* error messages from qhull code */ + int exitcode; /* 0 if no error from qhull */ + facetT *facet; /* set by FORALLfacets */ + int curlong, totlong; /* memory remaining after qh_memfreeshort */ + + coordT *point, *normp, *coordp, *feasiblep, *deladj; + int i, j, k; + boolT zerodiv; + float runsum; + char region; + /*coordT *points; + pointT *intpoints;*/ + + /* make point array from adjacency coordinates (add offset)*/ + /*points = (coordT *)malloc(4*numpoints*sizeof(coordT)); + if (points == NULL) { + printf("Unable to allocate points\n"); + exit(0); + }*/ + for (i=0; inormal; + feasiblep= qh feasible_point; + if (facet->offset < -qh MINdenom) { + for (k= qh hull_dim; k--; ) + *(coordp++)= (*(normp++) / - facet->offset) + *(feasiblep++); + }else { + for (k= qh hull_dim; k--; ) { + *(coordp++)= qh_divzero (*(normp++), facet->offset, qh MINdenom_1, + &zerodiv) + *(feasiblep++); + if (zerodiv) { + qh_memfree (point, qh normal_size); + printf("LABELprintinfinite\n"); + exit(0); + } + } + } + } + } + qh_freeqhull (!qh_ALL); + qh_memfreeshort (&curlong, &totlong); + + /* Now we calculate the volume */ + sprintf (flags, "qhull FA"); + exitcode= qh_new_qhull (dim, numpoints, intpoints, ismalloc, + flags, outfile, errfile); + + qh_getarea(qh facet_list); + *vol = qh totvol; + + qh_freeqhull (!qh_ALL); + qh_memfreeshort (&curlong, &totlong); + if (curlong || totlong) + fprintf (errfile, "qhull internal warning (vorvol): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong); + /*free(points); free(intpoints);*/ + + return exitcode; +} diff --git a/python_tools/pipeline_source/defaults.py b/python_tools/pipeline_source/defaults.py index a41070a..360bfd6 100644 --- a/python_tools/pipeline_source/defaults.py +++ b/python_tools/pipeline_source/defaults.py @@ -23,6 +23,8 @@ import os # ----------------------------------------------------------------------------- # DEFAULT CONFIGURATION +datasetName = "" + startCatalogStage = 1 endCatalogStage = 3 diff --git a/python_tools/pipeline_source/prepareCatalogs.in.py b/python_tools/pipeline_source/prepareCatalogs.in.py index 5bcaa7e..bc59fc7 100644 --- a/python_tools/pipeline_source/prepareCatalogs.in.py +++ b/python_tools/pipeline_source/prepareCatalogs.in.py @@ -29,6 +29,7 @@ import sys import void_python_tools as vp import argparse import imp +import subprocess # ----------------------------------------------------------------------------- @@ -81,31 +82,39 @@ def getSampleName(setName, redshift, useVel, iSlice=-1, iVol=-1): return sampleName #------------------------------------------------------------------------------ -def getNickName(sampleName): +def getNickName(setName, sampleName): splitName = sampleName.split('_') - + + nickName = datasetName + if "ss" in splitName[1]: - nickName = "Subsample = " + splitName[1].replace("ss","") + nickName += " SS " + splitName[1].replace("ss","") + #nickName = "Subsample = " + splitName[1].replace("ss","") if "pv" in splitName[2]: - nickName += ", z = " + splitName[3].replace("z","") + ", Pev.Vel." + nickName += ", z = " + splitName[3].replace("z","") + " (w/ PV)" else: nickName += ", z = " + splitName[2].replace("z","") elif "hod" in splitName[1]: - nickName = "HOD = " + splitName[2] + nickName += " HOD " + splitName[2] if "pv" in splitName[3]: - nickName += ", z = " + splitName[4].replace("z","") + ", Pec.Vel." + nickName += ", z = " + splitName[4].replace("z","") + " (w/ PV)" else: nickName += ", z = " + splitName[3].replace("z","") elif "halos" in splitName[1]: - nickName = "Halos, min mass = " + splitName[2].replace("min","") + if "none" in splitName[2]: + nickName += " All Halos" + else: + nickName += " Halos > " + splitName[2].replace("min","") if "pv" in splitName[3]: - nickName += ", z = " + splitName[4].replace("z","") + ", Pec.Vel." + nickName += ", z = " + splitName[4].replace("z","") + " (w/ PV)" else: nickName += ", z = " + splitName[3].replace("z","") else: nickName = sampleName + if "ran" in setName: nickName = "Random" + nickName + return nickName @@ -274,7 +283,7 @@ newSample.addStack({zMin}, {zMax}, 2*{minRadius}+18, 2*{minRadius}+24, True, Fal sampleName = getSampleName(setName, sliceMin, useVel, iSlice=iSlice, iVol=mySubvolume) - nickName = getNickName(sampleName) + nickName = getNickName(setName, sampleName) scriptFile.write(sampleInfo.format(dataFile=dataFileName, @@ -428,7 +437,7 @@ for iSubSample in xrange(len(subSamples)): rescale_position = hubble/1000./scale shift = lbox/2. rescale_velocity = 3.08567802e16/3.1558149984e16 - command = "%s %s x y z vz vy vx mass | awk '{print $1*%g+%g, $2*%g+%g, $3*%g+%g, $4*%g, $5*%g, $6*%g, $7}' > %s" % (SDFcvt_PATH, dataFile, + command = "%s -a 200000 %s x y z vz vy vx mass | awk '{print $1*%g+%g, $2*%g+%g, $3*%g+%g, $4*%g, $5*%g, $6*%g, $7}' > %s" % (SDFcvt_PATH, dataFile, rescale_position, shift, rescale_position, @@ -521,38 +530,51 @@ if (args.script or args.all) and haloFileBase != "": sys.stdout.flush() # estimate number of halos to get density - #if haloFileDummy == '': - # dataFile = catalogDir+haloFileBase+fileNums[0] - #else: - # dataFile = catalogDir+haloFileBase.replace(haloFileDummy, fileNums[0]) - # - #inFile = open(dataFile, 'r') - #numPart = 0 - #for (iLine, line) in enumerate(inFile): - # if iLine < haloFileNumComLines: continue - # line = line.split(haloFileColSep) - # if minHaloMass == "none" or float(line[haloFileMCol]) > minHaloMass: - # numPart += 1 - #inFile.close() + if haloFileDummy == '': + dataFile = catalogDir+haloFileBase+fileNums[0] + else: + dataFile = catalogDir+haloFileBase.replace(haloFileDummy, + fileNums[0]) + numPart = 0 + if dataFormat == "sdf": + SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" + if minHaloMass == "none": + command = "%s -a 200000 %s mass | wc" % (SDFcvt_PATH, dataFile) + else: + command = "%s -a 200000 %s mass | awk '{if ($1>%g) print $1}' | wc" % (SDFcvt_PATH, dataFile, minHaloMass) + numPart = subprocess.check_output(command, shell=True) + numPart = int(numPart.split()[0]) + else: + inFile = open(dataFile, 'r') + for (iHalo,line) in enumerate(inFile): + if iHalo < haloFileNumComLines: continue + line = line.split(haloFileColSep) + if minHaloMass == "none" or float(line[haloFileMCol]) > minHaloMass: + numPart += 1 + inFile.close() - #minRadius = 2*int(np.ceil(lbox/numPart**(1./3.))) - minRadies = 10 + minRadius = int(np.ceil(lbox/numPart**(1./3.))) - setName = prefix+"halos_min"+str(minHaloMass) + if minHaloMass != "none": + strMinHaloMass = "%.2e" % minHaloMass + else: + strMinHaloMass = "none" + + setName = prefix+"halos_min"+strMinHaloMass fileList = [] for (iRedshift, redshift) in enumerate(redshifts): sampleName = getSampleName(setName, redshift, False) outFileName = sampleName+".dat" fileList.append(outFileName) - writeScript(setName, prefix+"halos_min"+str(minHaloMass)+"_z", "multidark", + writeScript(setName, prefix+"halos_min"+strMinHaloMass+"_z", "multidark", scriptDir, catalogDir, fileNums, redshifts, numSubvolumes, numSlices, False, lbox, minRadius, omegaM, dataFileNameList = fileList) if doPecVel: - writeScript(setName, prefix+"halos_min"+str(minHaloMass)+"_z", + writeScript(setName, prefix+"halos_min"+strMinHaloMass+"_z", "multidark", scriptDir, catalogDir, fileNums, redshifts, @@ -583,6 +605,9 @@ if (args.halos or args.all) and haloFileBase != "": if "nhalos" in line: numPart = int(line.split()[3].strip(';')) break + if "npart" in line: + numPart = int(line.split()[3].strip(';')) + break inFile.close() else: for (iLine, line) in enumerate(inFile): @@ -592,7 +617,12 @@ if (args.halos or args.all) and haloFileBase != "": numPart += 1 inFile.close() - sampleName = prefix+"halos_min"+str(minHaloMass)+"_z"+redshifts[iRedshift] + if minHaloMass != "none": + strMinHaloMass = "%.2e" % minHaloMass + else: + strMinHaloMass = "none" + + sampleName = prefix+"halos_min"+strMinHaloMass+"_z"+redshifts[iRedshift] outFileName = catalogDir+"/"+sampleName+".dat" outFile = open(outFileName, 'w') outFile.write("%f\n" %(lbox)) @@ -604,8 +634,10 @@ if (args.halos or args.all) and haloFileBase != "": if dataFormat == "sdf": SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" - if minHaloMass == "none": minHaloMass = 0.0 - command = "%s %s mass id x y z vz vy vx | awk '{if ($1>%g) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, minHaloMass, outFileName ) + if minHaloMass == "none": + command = "%s -a 200000 %s mass ident x y z vz vy vx | awk '{print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, outFileName ) + else: + command = "%s -a 200000 %s mass ident x y z vz vy vx | awk '{if ($1>%g) print $2, $3, $4, $5, $6, $7, $8, $1}'>>%s" % (SDFcvt_PATH, dataFile, minHaloMass, outFileName ) os.system(command) outFile = open(outFileName, 'a') outFile.write("-99 -99 -99 -99 -99 -99 -99 -99\n") @@ -688,15 +720,38 @@ if (args.script or args.all) and haloFileBase != "": fileList.append(outFileName) print " ", thisHod['name'] + + # estimate number of halos to get density + if haloFileDummy == '': + dataFile = catalogDir+haloFileBase+fileNums[0] + else: + dataFile = catalogDir+haloFileBase.replace(haloFileDummy, + fileNums[0]) + numPart = 0 + if dataFormat == "sdf": + SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" + command = "%s -a 200000 %s mass | awk '{if ($1>%g) print $1}' | wc" % (SDFcvt_PATH, dataFile, thisHod['Mcut']) + numPart = subprocess.check_output(command, shell=True) + numPart = int(numPart.split()[0]) + else: + inFile = open(dataFile, 'r') + for (iHalo,line) in enumerate(inFile): + if iHalo < haloFileNumComLines: continue + line = line.split(haloFileColSep) + if float(line[haloFileMCol]) > thisHod['Mcut']: numPart += 1 + inFile.close() + + minRadius = int(np.ceil(lbox/numPart**(1./3.))) + setName = prefix+"hod_"+thisHod['name'] writeScript(setName, prefix+"hod_"+thisHod['name']+"_z", "multidark", scriptDir, catalogDir, fileNums, redshifts, - numSubvolumes, numSlices, False, lbox, 15, omegaM, + numSubvolumes, numSlices, False, lbox, minRadius, omegaM, dataFileNameList = fileList) if doPecVel: writeScript(setName, prefix+"hod_"+thisHod['name']+"_z", "multidark", scriptDir, catalogDir, fileNums, redshifts, - numSubvolumes, numSlices, True, lbox, 15, omegaM, + numSubvolumes, numSlices, True, lbox, minRadius, omegaM, dataFileNameList = fileList) if (args.hod or args.all) and haloFileBase != "": @@ -716,7 +771,7 @@ if (args.hod or args.all) and haloFileBase != "": inFile = haloFile outFile = haloFile+"_temp" SDFcvt_PATH = "@CMAKE_BINARY_DIR@/external/libsdf/apps/SDFcvt/SDFcvt.x86_64" - command = "%s %s mass x y z vx vy vz>>%s" % (SDFcvt_PATH, inFile, outFile) + command = "%s -a 200000 %s mass x y z vx vy vz>>%s" % (SDFcvt_PATH, inFile, outFile) os.system(command) haloFile = outFile @@ -730,13 +785,11 @@ if (args.hod or args.all) and haloFileBase != "": hubble=hubble, redshift=redshift, Mmin=thisHod['Mmin'], - #Mmin=1.23e13, M1=thisHod['M1'], sigma_logM=thisHod['sigma_logM'], alpha=thisHod['alpha'], Mcut=thisHod['Mcut'], galden=thisHod['galDens'], - #galden=0.000225, haloFile=haloFile, haloFileFormat=dataFormat, numPartPerSide=numPart**(1/3.), @@ -744,7 +797,13 @@ if (args.hod or args.all) and haloFileBase != "": workDir=catalogDir)) parFile.close() - os.system(hodPath+" "+parFileName+">& /dev/null") + tempFile = "./hod.out" + os.system(hodPath+" "+parFileName+">& " + tempFile) + for line in open(tempFile): + if "MLO" in line: + print " (minimum halo mass = ", line.split()[1], ")" + break + os.unlink(tempFile) sampleName = getSampleName(prefix+"hod_"+thisHod['name'], redshift, False) outFileName = catalogDir+"/"+sampleName+".dat" diff --git a/python_tools/void_python_tools/plotting/plotDefs.py b/python_tools/void_python_tools/plotting/plotDefs.py index 9872fc8..1980cdc 100644 --- a/python_tools/void_python_tools/plotting/plotDefs.py +++ b/python_tools/void_python_tools/plotting/plotDefs.py @@ -20,9 +20,8 @@ LIGHT_SPEED = 299792.458 colorList = ['r', 'b', 'g', 'y', 'c', 'm', - 'brown', 'grey', - 'darkred', 'orange', 'pink', 'darkblue', - 'lightblue', 'chocolate', + 'darkred', 'grey', + 'orange', 'darkblue', 'indigo', 'lightseagreen', 'maroon', 'olive', 'royalblue', 'palevioletred', 'seagreen', 'tomato', 'aquamarine', 'darkslateblue', diff --git a/python_tools/void_python_tools/xcor/xcorlib.py b/python_tools/void_python_tools/xcor/xcorlib.py index d5b489a..d70ebad 100644 --- a/python_tools/void_python_tools/xcor/xcorlib.py +++ b/python_tools/void_python_tools/xcor/xcorlib.py @@ -70,4 +70,4 @@ def shellavg( f, dx, Nmesh, Nbin = 100, xmin = 0., xmax = 1., scale = 'lin' ): fm = np.histogram(x, bins = bins, weights = f)[0]/Nm fs = np.sqrt((np.histogram(x, bins = bins, weights = f**2)[0]/Nm - fm**2)/(Nm-1)) - return (Nm, xm, fm, fs) \ No newline at end of file + return (Nm, xm, fm, fs)