guilhem_lavaux/vide_public
0
0
Fork 0
mirror of https://bitbucket.org/cosmicvoids/vide_public.git synced 2025-07-04 15:21:11 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

void overlap analysis appears to work

Browse source
This commit is contained in:
P.M. Sutter 2012-12-12 10:26:12 -06:00
parent 0b9673b24b
commit a5f88a8b3c
2 changed files with 166 additions and 68 deletions
Download patch file Download diff file Expand all files Collapse all files

226
c_tools/analysis/voidOverlap.cpp
View file

@ -16,6 +16,11 @@
#include "voidOverlap_conf.h" #include "voidOverlap_conf.h"
#include <vector> #include <vector>
#include <netcdfcpp.h> #include <netcdfcpp.h>
#include <iostream>
#include <sstream>
#include <fstream>
using namespace std;
typedef struct partStruct { typedef struct partStruct {
float x, y, z, volume; float x, y, z, volume;
@ -38,8 +43,10 @@ typedef struct voidStruct {
int voidID, numPart, numZones, coreParticle, zoneNumPart; int voidID, numPart, numZones, coreParticle, zoneNumPart;
float maxRadius, nearestMock, centralDen, redshift, redshiftInMpc; float maxRadius, nearestMock, centralDen, redshift, redshiftInMpc;
float nearestEdge; float nearestEdge;
float center[3], barycenter[3]; float barycenter[3];
std::vector<int> matches; std::vector<int> matches;
int biggestMatchID, numMatches;
float biggestMatchVol;
} VOID; } VOID;
typedef struct catalog { typedef struct catalog {
@ -53,11 +60,22 @@ typedef struct catalog {
void loadCatalog(const char *partFile, const char *volFile, void loadCatalog(const char *partFile, const char *volFile,
const char *voidFile, const char *zoneFile, const char *voidFile, const char *zoneFile,
const char *infoFile, const char *infoFile, const char *barycenterFile,
const char *zonePartFile, CATALOG& catalog); const char *zonePartFile, CATALOG& catalog);
float getDist(CATALOG& catalog1, CATALOG& catalog2, int& iVoid1, int& iVoid2);
int main(int argc, char **argv) { int main(int argc, char **argv) {
int p1, p2, iZ1, iZ2, iVoid1, iVoid2, iVoid, zoneID1, zoneID2;
int partID1, partID2;
int voidID1, voidID2;
bool periodicX=false, periodicY=false, periodicZ=false, match;
float dist[3], rdist, r1, r2;
FILE *fp;
CATALOG catalog1, catalog2;
// initialize arguments // initialize arguments
voidOverlap_info args; voidOverlap_info args;
voidOverlap_conf_params params; voidOverlap_conf_params params;
@ -79,25 +97,12 @@ int main(int argc, char **argv) {
return 1; return 1;
} }
char outVoid[200];
int p1, p2, iZ1, iZ2, iVoid1, iVoid2, iVoid, zoneID1, zoneID2;
int voidStartIndex, numVoids, numVoidsOut, targetVoidID, voidID1, voidID2;
FILE *fp;
float *temp, junk, voidVol;
long *temp2;
int junkInt, voidID, numPart, numZones, zoneID, partID;
char line[500], junkStr[10];
bool periodicX=false, periodicY=false, periodicZ=false, match;
float dist[3], rdist, r1, r2;
CATALOG catalog1, catalog2;
loadCatalog(args.partFile1_arg, args.volFile1_arg, args.voidFile1_arg, loadCatalog(args.partFile1_arg, args.volFile1_arg, args.voidFile1_arg,
args.zoneFile1_arg, args.infoFile1_arg, args.zoneFile1_arg, args.infoFile1_arg, args.barycenterFile1_arg,
args.zonePartFile1_arg, catalog1); args.zonePartFile1_arg, catalog1);
loadCatalog(args.partFile2_arg, args.volFile2_arg, args.voidFile2_arg, loadCatalog(args.partFile2_arg, args.volFile2_arg, args.voidFile2_arg,
args.zoneFile2_arg, args.infoFile2_arg, args.zoneFile2_arg, args.infoFile2_arg, args.barycenterFile2_arg,
args.zonePartFile2_arg, catalog2); args.zonePartFile2_arg, catalog2);
// check for periodic box // check for periodic box
@ -115,63 +120,111 @@ int main(int argc, char **argv) {
} }
printf(" Determining overlap...\n"); printf(" Determining overlap...\n");
//for (iVoid1 = 0; iVoid1 < 1; iVoid1++) {
for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) { 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, catalog1.numVoids);
voidID1 = catalog1.voids[iVoid1].voidID; 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++) {
for (iVoid2 = 0; iVoid2 < catalog2.numVoids; iVoid2++) { for (iVoid2 = 0; iVoid2 < catalog2.numVoids; iVoid2++) {
voidID2 = catalog2.voids[iVoid2].voidID; voidID2 = catalog2.voids[iVoid2].voidID;
match = false;
for (iZ1 = 0; iZ1 < catalog1.void2Zones[voidID1].numZones; iZ1++) {
if (match) break;
zoneID1 = catalog1.void2Zones[voidID1].zoneIDs[iZ1];
for (p1 = 0; p1 < catalog1.zones2Parts[zoneID1].numPart; p1++) {
if (match) break;
partID1 = catalog1.zones2Parts[zoneID1].partIDs[p1];
for (iZ2 = 0; iZ2 < catalog2.void2Zones[voidID2].numZones; iZ2++) { for (iZ2 = 0; iZ2 < catalog2.void2Zones[voidID2].numZones; iZ2++) {
if (match) break;
zoneID2 = catalog2.void2Zones[voidID2].zoneIDs[iZ2]; zoneID2 = catalog2.void2Zones[voidID2].zoneIDs[iZ2];
for (p2 = 0; p2 < catalog2.zones2Parts[zoneID2].numPart; p2++) {
match = false; for (p2 = 0; p2 < catalog2.zones2Parts[zoneID2].numPart; p2++) {
partID2 = catalog2.zones2Parts[zoneID2].partIDs[p2];
if (args.useID_flag) { if (args.useID_flag) {
if (catalog1.part[p1].uniqueID == if (catalog1.part[partID1].uniqueID ==
catalog2.part[p2].uniqueID) match = true; catalog2.part[partID2].uniqueID) match = true;
} else { } else {
dist[0] = catalog1.part[p1].x - catalog2.part[p2].x; dist[0] = fabs(catalog1.part[partID1].x -
dist[1] = catalog1.part[p1].y - catalog2.part[p2].y; catalog2.part[partID2].x);
dist[2] = catalog1.part[p1].z - catalog2.part[p2].z; dist[1] = fabs(catalog1.part[partID1].y -
catalog2.part[partID2].y);
dist[2] = fabs(catalog1.part[partID1].z -
catalog2.part[partID2].z);
if (periodicX) dist[0] = fmin(dist[0], if (periodicX) dist[0] = fmin(dist[0], 1.0 - dist[0]);
abs(catalog1.boxLen[0]-dist[0])); if (periodicY) dist[1] = fmin(dist[1], 1.0 - dist[1]);
if (periodicY) dist[1] = fmin(dist[1], if (periodicZ) dist[2] = fmin(dist[2], 1.0 - dist[2]);
abs(catalog1.boxLen[1]-dist[1]));
if (periodicZ) dist[2] = fmin(dist[2],
abs(catalog1.boxLen[2]-dist[2]));
rdist = sqrt(dist[0]*dist[0] + dist[1]*dist[1] + rdist = sqrt(dist[0]*dist[0] + dist[1]*dist[1] +
dist[2]*dist[2]); dist[2]*dist[2]);
r1 = pow(3./4./M_PI*catalog1.part[p1].volume, 1./3.); r1 = pow(3./4./M_PI*catalog1.part[partID1].volume /
r2 = pow(3./4./M_PI*catalog2.part[p2].volume, 1./3.); catalog1.numPartTot, 1./3.);
if (rdist <= r1 + r2) match = true; 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) { if (match) {
//printf("MATCHED TO %d\n", iVoid2);
catalog1.voids[iVoid1].matches.push_back(iVoid2); catalog1.voids[iVoid1].matches.push_back(iVoid2);
break; break;
} }
} // end p2 } // end p2
if (match) break; } // end iZ2
} // end iZ1 } // end p1
} // end iVoid2 } // end iZ1
} // end iVoid2
} // end p1
} // end iZ1
} // end iVoid1 } // end iVoid1
printf("HELLO: "); // output properties of matches
for (iVoid = 0; iVoid < catalog1.voids[0].matches.size(); iVoid++) { fp = fopen(args.outfile_arg, "w");
printf("%d ", catalog2.voids[iVoid].voidID); //for (iVoid1 = 0; iVoid1 < 1; iVoid1++) {
for (iVoid1 = 0; iVoid1 < catalog1.numVoids; iVoid1++) {
for (iVoid = 0; iVoid < catalog1.voids[iVoid1].matches.size(); iVoid++) {
iVoid2 = catalog1.voids[iVoid1].matches[iVoid];
rdist = getDist(catalog1, catalog2, iVoid1, iVoid2);
if (rdist/catalog1.voids[iVoid1].radius > 1 ||
rdist/catalog2.voids[iVoid2].radius > 1) {
catalog1.voids[iVoid].matches[iVoid] = -1;
continue;
}
if (catalog2.voids[iVoid2].vol >
catalog1.voids[iVoid1].biggestMatchVol) {
catalog1.voids[iVoid1].biggestMatchID = iVoid2;
catalog1.voids[iVoid1].biggestMatchVol = catalog2.voids[iVoid2].vol;
}
catalog1.voids[iVoid1].numMatches++;
//printf("CENTER %d %d %e %e %e\n", iVoid1, iVoid2, rdist, rdist/catalog1.voids[iVoid1].radius, rdist/catalog2.voids[iVoid2].radius);
}
int voidID = catalog1.voids[iVoid1].voidID;
if (catalog1.voids[iVoid1].numMatches > 0) {
iVoid2 = catalog1.voids[iVoid1].biggestMatchID;
float volRatio = catalog2.voids[iVoid2].vol /
catalog1.voids[iVoid1].vol;
rdist = getDist(catalog1, catalog2, iVoid1, iVoid2);
rdist /= catalog1.voids[iVoid1].radius;
fprintf(fp, "%d %e %e %d\n", voidID,
volRatio,
rdist,
catalog1.voids[iVoid1].numMatches);
} else {
fprintf(fp, "%d 0.0 0.0 0\n", voidID);
}
} }
fclose(fp);
printf("\n Well, bye.\n"); printf("\n Well, bye.\n");
return 0; return 0;
} // end main } // end main
@ -179,7 +232,7 @@ int main(int argc, char **argv) {
// ---------------------------------------------------------------------------- // ----------------------------------------------------------------------------
void loadCatalog(const char *partFile, const char *volFile, void loadCatalog(const char *partFile, const char *volFile,
const char *voidFile, const char *zoneFile, const char *voidFile, const char *zoneFile,
const char *infoFile, const char *infoFile, const char *barycenterFile,
const char *zonePartFile, CATALOG& catalog) { const char *zonePartFile, CATALOG& catalog) {
int i, p, numPartTot, numZonesTot, dummy, iVoid, iZ, numVolTot; int i, p, numPartTot, numZonesTot, dummy, iVoid, iZ, numVolTot;
@ -191,7 +244,7 @@ void loadCatalog(const char *partFile, const char *volFile,
char line[500], junkStr[10]; char line[500], junkStr[10];
float ranges[3][2]; float ranges[3][2];
printf("\n Loading info...\n"); printf("Loading info...\n");
NcFile f_info(infoFile); NcFile f_info(infoFile);
ranges[0][0] = f_info.get_att("range_x_min")->as_double(0); ranges[0][0] = f_info.get_att("range_x_min")->as_double(0);
ranges[0][1] = f_info.get_att("range_x_max")->as_double(0); ranges[0][1] = f_info.get_att("range_x_max")->as_double(0);
@ -206,13 +259,14 @@ void loadCatalog(const char *partFile, const char *volFile,
f_info.close(); f_info.close();
// read in all particle positions // read in all particle positions
printf("\n Loading particles...\n"); printf("Loading particles...\n");
fp = fopen(partFile, "r"); fp = fopen(partFile, "r");
fread(&dummy, 1, 4, fp); fread(&dummy, 1, 4, fp);
fread(&catalog.numPartTot, 1, 4, fp); fread(&numPartTot, 1, 4, fp);
fread(&dummy, 1, 4, fp); fread(&dummy, 1, 4, fp);
catalog.part.resize(numPartTot); catalog.part.resize(numPartTot);
catalog.numPartTot = numPartTot;
temp = (float *) malloc(numPartTot * sizeof(float)); temp = (float *) malloc(numPartTot * sizeof(float));
temp2 = (long *) malloc(numPartTot * sizeof(long)); temp2 = (long *) malloc(numPartTot * sizeof(long));
@ -275,24 +329,47 @@ void loadCatalog(const char *partFile, const char *volFile,
fgets(line, sizeof(line), fp); fgets(line, sizeof(line), fp);
catalog.voids.resize(catalog.numVoids); catalog.voids.resize(catalog.numVoids);
iVoid = 0; i = 0;
while (fgets(line, sizeof(line), fp) != NULL) { while (fgets(line, sizeof(line), fp) != NULL) {
scanf(line, "%d %d %d %f %f %d %d %f %d %f %f\n", &iVoid, &voidID, sscanf(line, "%d %d %d %f %f %d %d %f %d %f %f\n", &iVoid, &voidID,
&coreParticle, &coreDens, &zoneVol, &zoneNumPart, &numZones, &coreParticle, &coreDens, &zoneVol, &zoneNumPart, &numZones,
&voidVol, &numPart, &densCon, &voidProb); &voidVol, &numPart, &densCon, &voidProb);
catalog.voids[iVoid].coreParticle = coreParticle; catalog.voids[i].coreParticle = coreParticle;
catalog.voids[iVoid].zoneNumPart = zoneNumPart; catalog.voids[i].zoneNumPart = zoneNumPart;
catalog.voids[iVoid].coreDens = coreDens; catalog.voids[i].coreDens = coreDens;
catalog.voids[iVoid].zoneVol = zoneVol; catalog.voids[i].zoneVol = zoneVol;
catalog.voids[iVoid].voidID = voidID; catalog.voids[i].voidID = voidID;
catalog.voids[iVoid].vol = voidVol; catalog.voids[i].vol = voidVol;
catalog.voids[iVoid].numPart = numPart; catalog.voids[i].numPart = numPart;
catalog.voids[iVoid].numZones = numZones; catalog.voids[i].numZones = numZones;
catalog.voids[iVoid].densCon = densCon; catalog.voids[i].densCon = densCon;
catalog.voids[iVoid].voidProb = voidProb; catalog.voids[i].voidProb = voidProb;
catalog.voids[iVoid].radius = pow(voidVol/volNorm*3./4./M_PI, 1./3.); 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].numMatches = 0;
i++;
}
fclose(fp);
printf(" Read %d voids.\n", catalog.numVoids);
printf(" Loading barycenters\n");
fp = fopen(barycenterFile, "r");
float tempBary[3];
iVoid = 0;
while (fgets(line, sizeof(line), fp) != NULL) {
sscanf(line, "%d %f %f %f\n", &voidID, &tempBary[0], &tempBary[1],
&tempBary[2]);
tempBary[0] = (tempBary[0] - ranges[0][0])/catalog.boxLen[0];
tempBary[1] = (tempBary[1] - ranges[1][0])/catalog.boxLen[1];
tempBary[2] = (tempBary[2] - ranges[2][0])/catalog.boxLen[2];
catalog.voids[iVoid].barycenter[0] = tempBary[0];
catalog.voids[iVoid].barycenter[1] = tempBary[1];
catalog.voids[iVoid].barycenter[2] = tempBary[2];
iVoid++; iVoid++;
} }
fclose(fp); fclose(fp);
@ -304,7 +381,7 @@ void loadCatalog(const char *partFile, const char *volFile,
catalog.void2Zones.resize(catalog.numZonesTot); catalog.void2Zones.resize(catalog.numZonesTot);
for (iZ = 0; iZ < numZonesTot; iZ++) { for (iZ = 0; iZ < catalog.numZonesTot; iZ++) {
fread(&numZones, 1, 4, fp); fread(&numZones, 1, 4, fp);
catalog.void2Zones[iZ].numZones = numZones; catalog.void2Zones[iZ].numZones = numZones;
@ -323,6 +400,7 @@ void loadCatalog(const char *partFile, const char *volFile,
fread(&numZonesTot, 1, 4, fp); fread(&numZonesTot, 1, 4, fp);
catalog.zones2Parts.resize(numZonesTot); catalog.zones2Parts.resize(numZonesTot);
for (iZ = 0; iZ < numZonesTot; iZ++) { for (iZ = 0; iZ < numZonesTot; iZ++) {
fread(&numPart, 1, 4, fp); fread(&numPart, 1, 4, fp);
@ -337,3 +415,21 @@ void loadCatalog(const char *partFile, const char *volFile,
fclose(fp); fclose(fp);
} // end loadCatalog } // end loadCatalog
// ----------------------------------------------------------------------------
float getDist(CATALOG& catalog1, CATALOG& catalog2, int& iVoid1, int& iVoid2) {
float rdist, dist[3];
dist[0] = fabs(catalog1.voids[iVoid1].barycenter[0] -
catalog2.voids[iVoid2].barycenter[0]);
dist[1] = fabs(catalog1.voids[iVoid1].barycenter[1] -
catalog2.voids[iVoid2].barycenter[1]);
dist[2] = fabs(catalog1.voids[iVoid1].barycenter[2] -
catalog2.voids[iVoid2].barycenter[2]);
rdist = sqrt(dist[0]*dist[0] + dist[1]*dist[1] + dist[2]*dist[2]);
return rdist;
} // end getDist

4
c_tools/analysis/voidOverlap.ggo
View file

@ -10,6 +10,7 @@ option "voidFile1" - "Void info file for catalog 1" string yes
option "infoFile1" - "Extra info file for catalog 1" string yes option "infoFile1" - "Extra info file for catalog 1" string yes
option "zoneFile1" - "Zone file for catalog 1" string yes option "zoneFile1" - "Zone file for catalog 1" string yes
option "zonePartFile1" - "Zone-particle file for catalog 1" string yes option "zonePartFile1" - "Zone-particle file for catalog 1" string yes
option "barycenterFile1" - "Barycenter file for catalog 1" string yes
# void data for file 2 # void data for file 2
option "partFile2" - "Particle file for catalog 2" string yes option "partFile2" - "Particle file for catalog 2" string yes
@ -18,8 +19,9 @@ option "voidFile2" - "Void info file for catalog 2" string yes
option "infoFile2" - "Extra info file for catalog 2" string yes option "infoFile2" - "Extra info file for catalog 2" string yes
option "zoneFile2" - "Zone file for catalog 2" string yes option "zoneFile2" - "Zone file for catalog 2" string yes
option "zonePartFile2" - "Zone-particle file for catalog 2" string yes option "zonePartFile2" - "Zone-particle file for catalog 2" string yes
option "barycenterFile2" - "Barycenter file for catalog 2" string yes
# options # options
option "prefix" - "Prefix for output" string yes option "outfile" - "Output file" string yes
option "useID" - "Use unique catalog ID to match voids; otherwise use volumes" flag off option "useID" - "Use unique catalog ID to match voids; otherwise use volumes" flag off
option "periodic" - "Set of edges which are periodic" string optional default="xy" option "periodic" - "Set of edges which are periodic" string optional default="xy"