Tesllation seems stable with new boundary handling procedure. Rooted out strange adjacency record-keeping due to some galaxies connecting to guard particles. Zobov still generates warning about guard particle encounters, but this is fine for observational datasets because we don't trust the tesselation at the edge galaxies anyway.

This commit is contained in:
Paul M. Sutter 2025-01-16 18:40:54 +08:00
parent 091cf5d1dc
commit 5d93a8a737
8 changed files with 63 additions and 12 deletions

View file

@ -90,9 +90,10 @@ void readAdjacencyFile(const string& adjfile, PARTICLE*& p, pid_t& np)
}
else{
p[i].adj[p[i].ncnt] = j;
p[j].adj[p[j].ncnt] = i;
p[i].ncnt++;
p[j].ncnt++; }
p[j].adj[p[j].ncnt] = i;
p[j].ncnt++;
}
}
else
{

View file

@ -202,6 +202,44 @@ void placeGalaxiesInCube(NYU_VData& data, ParticleData& output_data,
% (1e-2*output_data.box[1][1]) % (1e-2*output_data.box[1][0])
% (1e-2*output_data.box[2][1]) % (1e-2*output_data.box[2][0]) << endl;
/*
// TEST DEBUG
int nNewPerSide = 100;
Position p;
int idx = data.size() + 1;
float xwidth = output_data.box[0][0] - output_data.box[0][1];
float ywidth = output_data.box[1][0] - output_data.box[1][1];
float zwidth = output_data.box[2][0] - output_data.box[2][1];
float dx = xwidth / nNewPerSide;
float dy = ywidth / nNewPerSide;
float dz = zwidth / nNewPerSide;
for (int i = 0; i < nNewPerSide; i++) {
for (int j = 0; j < nNewPerSide; j++) {
for (int k = 0; k < nNewPerSide; k++) {
p.xyz[0] = output_data.box[0][1] + i*dx;
p.xyz[1] = output_data.box[1][1] + j*dy;
p.xyz[2] = output_data.box[2][1] + k*dz;
cout << "INSERT " << p.xyz[0] << ", " << p.xyz[1] << ", " << p.xyz[2] << endl;
cout << "COMPARING " << output_data.pos[i].xyz[0] << endl;
output_data.pos.push_back(p);
output_data.id_gal.push_back(idx);
output_data.ra.push_back(-1);
output_data.dec.push_back(-1);
output_data.redshift.push_back(-1);
output_data.uniqueID.push_back(-1);
idx++;
}
}
}
*/
// END TEST DEBUG
}
void flagEdgeGalaxies(prepObservation_info& args ,
@ -230,8 +268,10 @@ void flagEdgeGalaxies(prepObservation_info& args ,
// write a small text file with galaxy position (for diagnostic purposes)
// TODO - remove this
FILE *fp;
fp = fopen("galaxies.txt", "w");
//for (int i = 0; i < output_data.id_gal.size(); i++) {
for (int i = 0; i < data.size(); i++) {
Position& p = output_data.pos[i];
fprintf(fp, "%e %e %e\n",
@ -326,7 +366,7 @@ void saveForZobov(ParticleData& pdata, const string& fname,
f.writeInt32(pdata.pos.size());
f.endCheckpoint();
cout << format("Saving %d galaxies...") % pdata.pos.size() << endl;
for (int j = 0; j < 3; j++)
{
cout << format("Writing %c components...") % axis[j] << endl;

View file

@ -162,7 +162,7 @@ int main(int argc,char **argv) {
/*if (p[i].ncnt != p[i].nadj) {*/
// END PMS
p[i].nadj = p[i].ncnt;
printf("We didn't get all of %d's adj's; %d != %d.\n",i,nin,p[i].nadj);
printf("jovoz: We didn't get all of %d's adj's; %d != %d.\n",i,nin,p[i].nadj);
/*exit(0);*/
}
}

View file

@ -305,8 +305,16 @@ int main(int argc, char *argv[]) {
exitcode = vorvol(deladjs, points, intpoints, adjs[i].nadj, &(vols[i]));
vols[i] *= (float)np;
if (i % 1000 == 0)
printf("%d: %d, %f\n",i,adjs[i].nadj,vols[i]);
//if (i % 1000 == 0)
// printf("%d: %d, %f\n",i,adjs[i].nadj,vols[i]);
}
// PMS - reset number of adjancies to not include links to border guards
for (i=0; i<nvp; i++) {
int nActual = 0;
for (j = 0; j < adjs[i].nadj; j++)
if (adjs[i].adj[j] < nvp) nActual++;
adjs[i].nadj = nActual;
}
/* Get the adjacencies back to their original values */

View file

@ -21,8 +21,8 @@
import os
from numpy import abs
import matplotlib as mpl
mpl.use('Agg')
#import matplotlib as mpl
#mpl.use('Agg')
LIGHT_SPEED = 299792.458

View file

@ -118,6 +118,7 @@ def launchPrep(sample, binPath, workDir=None, inputDataDir=None,
# flag edge galaxies with python routine for now
galFile = outputDir+"/zobov_slice_"+sampleName
os.system("mv %s %s" % ("galaxies.txt", outputDir))
#galFile = outputDir + "galaxies.txt"
edgeGalFile = outputDir + "/galaxy_edge_flags.txt"
contourFile = outputDir + "/contour_map.fits"
@ -481,7 +482,7 @@ def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None,
edgeFlags = np.loadtxt(edgeFile, dtype=np.int32)
# set edge galaxy volumes to nearly 0 (implying very high density)
vols[ edgeFlags>0 ] = 1.e-4
vols[ edgeFlags>0 ] = 1.e-10
volFile = outputDir+"/vol_weighted_"+sampleName+".dat"
with open(volFile, mode='wb') as File:
@ -493,7 +494,6 @@ def launchZobov(sample, binPath, outputDir=None, logDir=None, continueRun=None,
else:
volFileToUse = outputDir+"/vol_"+sampleName+".dat"
cmd = [binPath+"/jozov2", \
outputDir+"/adj_"+sampleName+".dat", \
volFileToUse, \

View file

@ -146,6 +146,7 @@ def findEdgeGalaxies(galFile, maskFile, edgeGalFile, contourFile,
# load in galaxies
#galPos = np.genfromtxt(galFile)
# TODO - WHY IS THIS FASTER THAN np.column_stack???
galPos = np.genfromtxt(outputDir+"/galaxies.txt")
with open(galFile, 'rb') as File:
chk = np.fromfile(File, dtype=np.int32,count=1)

View file

@ -484,6 +484,7 @@ def loadVoidCatalog(sampleDir, dataPortion="central", loadParticles=True,
catalog.part[p].adjs.append(pAdj)
catalog.part[pAdj].adjs.append(p)
print(" Sanity checking adjacenies...")
for p in range(numPart):
catalog.part[p].nadjs = len(catalog.part[p].adjs)