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

FIXED zobov bug when only using one subdomain; added zobov buffer as tunable parameter

Browse source
This commit is contained in:
Paul M. Sutter 2025-05-22 09:57:08 -04:00
parent 8730193e71
commit 23d665f7bd
6 changed files with 139 additions and 83 deletions
Download patch file Download diff file Expand all files Collapse all files

72
c_source/zobov/voz1b1.c
View file

@ -17,7 +17,7 @@ int main(int argc, char *argv[]) {
coordT rtemp[3], *parts;
coordT deladjs[3*MAXVERVER], points[3*MAXVERVER];
pointT intpoints[3*MAXVERVER];
FILE *pos, *out;
FILE *pos, *out, *testFile;
char *posfile, outfile[200], *suffix, *outDir;
PARTADJ *adjs;
float *vols;
@ -33,8 +33,11 @@ int main(int argc, char *argv[]) {
realT c[3];
int b[3];
double totalvol;
int isObs = 0;
if (argc != 10) {
printf("Routine: voz1b1\n");
if (argc != 11) {
printf("Wrong number of arguments.\n");
printf("arg1: position file\n");
printf("arg2: border size\n");
@ -43,6 +46,7 @@ int main(int argc, char *argv[]) {
printf("arg5: number of divisions\n");
printf("arg6-8: b[0-2]\n\n");
printf("arg9: output directory\n");
printf("arg10: observation flag\n");
exit(0);
}
posfile = argv[1];
@ -79,8 +83,15 @@ int main(int argc, char *argv[]) {
exit(0);
}
outDir = argv[9];
if (sscanf(argv[10],"%d",&isObs) != 1) {
printf("That's no observation mode; try again.\n");
exit(0);
}
printf("Working with subdomain b=[%d,%d,%d]\n", b[0], b[1], b[2]);
/* Boxsize should be the range in r, yielding a range 0-1 */
printf("Reading particle file...\n");
np = posread(posfile,&r,1./boxsize);
printf("%d particles\n",np);fflush(stdout);
@ -90,7 +101,8 @@ int main(int argc, char *argv[]) {
if (r[i][1]<ymin) ymin = r[i][1]; if (r[i][1]>ymax) ymax = r[i][1];
if (r[i][2]<zmin) zmin = r[i][2]; if (r[i][2]>zmax) zmax = r[i][2];
}
printf("np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np,xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout);
printf("Full box particle stats: np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np,
xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout);
width = 1./(float)numdiv;
width2 = 0.5*width;
@ -108,7 +120,7 @@ int main(int argc, char *argv[]) {
exit(0);
}
g = (bf / 2.)*(1. + sqrt(1 - 2.*s*s/(bf*bf)));
printf("s = %f, bf = %f, g = %f.\n",s,bf,g);
printf("Guard particle information: s = %f, bf = %f, g = %f.\n",s,bf,g);
fflush(stdout);
@ -118,11 +130,10 @@ int main(int argc, char *argv[]) {
exit(1);
}
DL c[d] = ((float)b[d])*width;
printf("c: %f,%f,%f\n",c[0],c[1],c[2]);
DL c[d] = ((float)b[d]+0.5)*width;
printf("Counting particles in subdomain...\n");
/* Assign temporary array*/
nvpbuf = 0; /* Number of particles to tesselate, including
buffer */
nvpbuf = 0; /* Number of particles to tesselate, includng ibuffer */
nvp = 0; /* Without the buffer */
for (i=0; i<np; i++) {
isitinbuf = 1;
@ -139,8 +150,7 @@ int main(int argc, char *argv[]) {
if (isitinmain) nvp++;
}
nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard
points */
nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard points */
parts = (coordT *)malloc(3*nvpbuf*sizeof(coordT));
orig = (pid_t *)malloc(nvpbuf*sizeof(pid_t));
@ -156,6 +166,7 @@ int main(int argc, char *argv[]) {
exit(1);
}
printf("Placing particles in subdomain...\n");
nvp = 0; nvpall = 0; /* nvp = number of particles without buffer */
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
for (i=0; i<np; i++) {
@ -180,8 +191,13 @@ int main(int argc, char *argv[]) {
if (rtemp[2] > zmax) zmax = rtemp[2];
}
}
printf("nvp = %d\n",nvp);
printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
printf("Subdomain: b=(%d,%d,%d), c=(%f,%f,%f)\n",
b[0],b[1],b[2],c[0],c[1],c[2]);
printf(" nvp=%d, nvpbuf=%d\n", nvp,nvpbuf);
printf(" Particle ranges: x: %f,%f; y: %f,%f; z:%f,%f\n",
xmin,xmax,ymin,ymax,zmin,zmax);
printf("Adding particles to buffer regions...\n");
nvpbuf = nvp;
for (i=0; i<np; i++) {
isitinbuf = 1;
@ -195,9 +211,6 @@ int main(int argc, char *argv[]) {
isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
}
if (isitinbuf && !isitinmain) {
/*printf("%3.3f ",sqrt(rtemp[0]*rtemp[0] + rtemp[1]*rtemp[1] +
rtemp[2]*rtemp[2]));
printf("|%2.2f,%2.2f,%2.2f,%f,%f",r[i][0],r[i][1],r[i][2],width2,totwidth2);*/
parts[3*nvpbuf] = rtemp[0];
parts[3*nvpbuf+1] = rtemp[1];
parts[3*nvpbuf+2] = rtemp[2];
@ -212,16 +225,19 @@ int main(int argc, char *argv[]) {
if (rtemp[2] > zmax) zmax = rtemp[2];
}
}
printf("nvpbuf = %d\n",nvpbuf);
printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
nvpall = nvpbuf;
predict = pow(width+2.*bf,3)*(float)np;
printf("There should be ~ %f points; there are %d\n",predict,nvpbuf);
printf("Total particles including buffer = %d (predicted ~%f)\n",
nvpbuf, predict);
printf(" Particle ranges including buffer: x: %f,%f; y: %f,%f; z:%f,%f\n",
xmin,xmax,ymin,ymax,zmin,zmax);
for (i=0;i<np;i++) free(r[i]);
free(r);
/* Add guard points */
printf("Adding guard points to subdomain faces...\n");
for (i=0; i<NGUARD+1; i++) {
for (j=0; j<NGUARD+1; j++) {
/* Bottom */
@ -262,6 +278,23 @@ int main(int argc, char *argv[]) {
nvpall++;
}
}
/*
// TEST TEMP
testFile = fopen("test.txt", "a");
for (i=0; i<nvpall; i++){
int itype = 0;
if (i > nvp) itype = 1; // it's a buffer particle
if (i > nvpbuf) itype = 2; // it's a guard particle
fprintf(testFile, "%f %f %f %d %f %f %f\n", c[0], c[1], c[2],
itype,
parts[3*i],
parts[3*i+1],
parts[3*i+2]);
}
fclose(testFile);
*/
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
for (i=nvpbuf;i<nvpall;i++) {
if (parts[3*i] < xmin) xmin = parts[3*i];
@ -274,7 +307,8 @@ int main(int argc, char *argv[]) {
printf("Added guard points to total %d points (should be %d)\n",nvpall,
nvpbuf + 6*(NGUARD+1)*(NGUARD+1));
printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax);
printf("New particle ranges including guards: x: %f,%f; y: %f,%f; z:%f,%f\n",
xmin,xmax,ymin,ymax,zmin,zmax);
/* Do tesselation*/
printf("File read. Tessellating ...\n"); fflush(stdout);

45
c_source/zobov/vozinit.c
View file

@ -28,8 +28,11 @@ int main(int argc, char *argv[]) {
int numGuards;
int b[3];
int numThreads;
int isObs = 0;
if (argc != 9) {
printf("Routine: vozinit\n");
if (argc != 10) {
printf("Wrong number of arguments.\n");
printf("arg1: position file\n");
printf("arg2: buffer size (default 0.1)\n");
@ -39,6 +42,7 @@ int main(int argc, char *argv[]) {
printf("arg6: number of parallel threads\n");
printf("arg7: location of voboz executables\n");
printf("arg8: output directory\n");
printf("arg9: observation flag (set to 1 for observation mode)\n");
exit(0);
}
posfile = argv[1];
@ -70,8 +74,13 @@ int main(int argc, char *argv[]) {
}
vobozPath = argv[7];
outDir = argv[8];
if (sscanf(argv[9],"%d",&isObs) != 1) {
printf("That's no observation mode; try again.\n");
exit(0);
}
/* Read the position file */
printf("Reading particle file and calculating box sizes...\n");
np = posread(posfile,&rfloat,1./boxsize);
/* Boxsize should be the range in r, yielding a range 0-1 */
@ -97,6 +106,8 @@ int main(int argc, char *argv[]) {
nvpmax = 0; nvpbufmax = 0; nvpmin = np; nvpbufmin = np;
// count particles in each sub-domain
printf("Counting particles in each subdomain...\n");
for (b[0] = 0; b[0] < numdiv; b[0]++) {
c[0] = ((float)b[0]+0.5)*width;
for (b[1] = 0; b[1] < numdiv; b[1]++) {
@ -105,21 +116,21 @@ int main(int argc, char *argv[]) {
c[2] = ((float)b[2]+0.5)*width;
nvp = 0; /* Number of particles excluding buffer */
nvpbuf = 0; /* Number of particles to tesselate, including
buffer */
nvpbuf = 0; /* Number of particles to tesselate, including buffer */
xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF;
for (i=0; i<np; i++) {
isitinbuf = 1; isitinmain = 1;
for (d=0; d<3; d++) {
rtemp[d] = rfloat[i][d] - c[d];
//if (isObs == 0) {
if (rtemp[d] > 0.5) rtemp[d] --;
if (rtemp[d] < -0.5) rtemp[d] ++;
//}
isitinbuf = isitinbuf && (fabs(rtemp[d]) < totwidth2);
isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
}
if (isitinbuf) {
nvpbuf++;
}
if (isitinbuf) { nvpbuf++; }
if (isitinmain) {
nvp++;
if (rtemp[0] < xmin) xmin = rtemp[0];
@ -130,13 +141,17 @@ int main(int argc, char *argv[]) {
if (rtemp[2] > zmax) zmax = rtemp[2];
}
}
if (nvp > nvpmax) nvpmax = nvp;
if (nvpbuf > nvpbufmax) nvpbufmax = nvpbuf;
if (nvp < nvpmin) nvpmin = nvp;
if (nvpbuf < nvpbufmin) nvpbufmin = nvpbuf;
printf("b=(%d,%d,%d), c=(%f,%f,%f), nvp=%d, nvpbuf=%d\n",
b[0],b[1],b[2],c[0],c[1],c[2],nvp,nvpbuf);
printf("Subdomain: b=(%d,%d,%d), c=(%f,%f,%f)\n",
b[0],b[1],b[2],c[0],c[1],c[2]);
printf(" nvp=%d, nvpbuf=%d\n", nvp,nvpbuf);
printf(" Particle ranges: x: %f,%f; y: %f,%f; z:%f,%f\n",
xmin,xmax,ymin,ymax,zmin,zmax);
}
}
}
@ -144,9 +159,9 @@ int main(int argc, char *argv[]) {
printf("Nvpbuf range: %d,%d\n",nvpbufmin,nvpbufmax);
numGuards = 6*(NGUARD+1)*(NGUARD+1);
printf("Num guards = %d\n", numGuards);
printf("Total max particles: %d\n" , nvpbufmax+numGuards);
if (nvpbufmax+numGuards >= QHULL_MAX_PARTICLES)
{
if (nvpbufmax+numGuards >= QHULL_MAX_PARTICLES) {
printf("Too many particles to tesselate per division (Qhull will overflow). Increase divisions or reduce buffer size.\n");
fflush(stdout);
exit(1);
@ -166,21 +181,17 @@ int main(int argc, char *argv[]) {
for (b[0]=0;b[0]<numdiv; b[0]++) {
for (b[1] = 0; b[1] < numdiv; b[1]++) {
for (b[2] = 0; b[2] < numdiv; b[2]++) {
//fprintf(scr,"%s/../c_tools/zobov2/voz1b1/voz1b1_2 %s %f %f %f %f %s %d %d %d %d %d %d %s&\n",
// vobozPath,
// posfile,border,boxsize,boxsize,boxsize,suffix,numdiv,numdiv, numdiv,b[0],b[1],b[2],
// outDir);
fprintf(scr,"%s/voz1b1 %s %f %f %s %d %d %d %d %s&\n",
fprintf(scr,"%s/voz1b1 %s %f %f %s %d %d %d %d %s %d\n",
vobozPath,
posfile,border,boxsize,suffix,numdiv,b[0],b[1],b[2],
outDir);
outDir, isObs);
p++;
if ((p == numThreads)) { fprintf(scr, "wait\n"); p = 0; }
}
}
}
fprintf(scr,"wait\n");
fprintf(scr,"%s/voztie %d %s %s %d\n", vobozPath, numdiv,suffix, outDir);
fprintf(scr,"%s/voztie %d %s %s %d\n", vobozPath, numdiv,suffix, outDir, isObs);
fclose(scr);
sprintf(systemstr,"chmod u+x %s",scrfile);

12
c_source/zobov/voztie.c
View file

@ -33,13 +33,14 @@ int main(int argc, char *argv[]) {
double avgnadj, avgvol;
int numRemoved = 0;
int isObs = 0;
if (argc != 5) {
printf("Wrong number of arguments.\n");
printf("Voztie: Wrong number of arguments.\n");
printf("arg1: number of divisions (default 2)\n");
printf("arg2: suffix describing this run\n");
printf("arg3: file directory\n");
printf("arg4: Beginning index of mock particles\n\n");
printf("arg5: observation flag\n");
exit(0);
}
if (sscanf(argv[1],"%d",&numdiv) != 1) {
@ -54,6 +55,11 @@ int main(int argc, char *argv[]) {
suffix = argv[2];
outDir = argv[3];
if (sscanf(argv[4],"%d",&isObs) != 1) {
printf("That's no observation flag; try again.\n");
exit(0);
}
np = -1; nvpmax = -1; nvpsum = 0;
for (i = 0; i < numdiv; i++) {
@ -204,7 +210,7 @@ int main(int argc, char *argv[]) {
npnotdone = 0; avgnadj = 0.; avgvol = 0.;
for (p=0;p<np;p++) {
// PMS
// PMS - TOTO KEEP THIS HERE???
if (vols[p] == 1.e-29) continue;
// END PMS
if (vols[p] == -1.) npnotdone++;

10
python_source/backend/launchers.py
View file

@ -353,6 +353,7 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
# -----------------------------------------------------------------------------
def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None,
numZobovDivisions=None, numZobovThreads=None,
zobovBuffer=0.1,
mergingThreshold=0.2):
startTime = time.time()
@ -382,9 +383,14 @@ def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None,
if os.access(outputDir+"/voidDesc_"+sampleName+".out", os.F_OK):
os.unlink(outputDir+"/voidDesc_"+sampleName+".out")
cmd = [binPath+"/vozinit", datafile, "0.1", "1.0", str(numZobovDivisions), \
isObs = 0
if sample.dataType == "observation":
isObs = 1
cmd = [binPath+"/vozinit", datafile, str(zobovBuffer), \
"1.0", str(numZobovDivisions), \
"_"+sampleName, str(numZobovThreads), \
binPath, outputDir]
binPath, outputDir, str(isObs)]
log = open(logFile, 'w')
subprocess.call(cmd, stdout=log, stderr=log)
log.close()

1
python_source/vide_pipeline/__main__.py
View file

@ -101,6 +101,7 @@ for sample in dataSampleList:
launchZobov(sample, ZOBOV_PATH, outputDir=outputDir, logDir=logDir,
continueRun=continueRun, numZobovDivisions=numZobovDivisions,
numZobovThreads=numZobovThreads,
zobovBuffer=zobovBuffer,
mergingThreshold=mergingThreshold)
# -------------------------------------------------------------------------

4
python_source/voidUtil/catalogUtil.py
View file

@ -258,9 +258,7 @@ def loadVoidCatalog(sampleDir,
# clearNearBoundaries: remove voids where the maximum extent is
# greater than the distance to nearest edge
# maxHighCentralDen: if != -1, filters based on based on central density
# NOTE: we are moving away from the cumbersome void catalog outputs of
# older versions, and will eventually just output a single catalog.
#
# To replicate the old "central" catalog, choose:
# clearEdges = True
# clearTree = True