Merge branch 'master' of bitbucket.org:cosmicvoids/void_identification

This commit is contained in:
P.M. Sutter 2013-03-20 07:19:34 -05:00
commit 7c172d5642
8 changed files with 77 additions and 56 deletions

View file

@ -274,7 +274,7 @@ void selectBox(SimuData *simu, std::vector<long>& targets, generateMock_info& ar
numAccepted++;
}
}
cout << "Num accepted here = " << numAccepted << " / input = " << simu->NumPart << endl;
cout << "SELECTBOX: Num accepted here = " << numAccepted << " / input = " << simu->NumPart << " (after resubsampling)" << endl;
}
class PreselectParticles: public SimulationPreprocessor
@ -385,7 +385,7 @@ void buildBox(SimuData *simu, long num_targets, long loaded,
}
}
void saveBox(SimuData *&boxed, const std::string& outbox)
void saveBox(SimuData *&boxed, const std::string& outbox, generateMock_info& args_info)
{
double *ranges = boxed->as<double>("ranges");
NcFile f(outbox.c_str(), NcFile::Replace, 0, 0, NcFile::Netcdf4);
@ -395,6 +395,11 @@ void saveBox(SimuData *&boxed, const std::string& outbox)
int num_snapshots = *boxed->as<int>("num_snapshots");
long *uniqueID = boxed->as<long>("uniqueID");
if (!f.is_valid())
{
cerr << "Could not create parameter file '" << outbox << "'. Aborting." << endl;
exit(1);
}
f.add_att("range_x_min", ranges[0]);
f.add_att("range_x_max", ranges[1]);
f.add_att("range_y_min", ranges[2]);
@ -403,6 +408,7 @@ void saveBox(SimuData *&boxed, const std::string& outbox)
f.add_att("range_z_max", ranges[5]);
f.add_att("mask_index", -1);
f.add_att("is_observation", 0);
f.add_att("data_subsampling", args_info.subsample_arg);
NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart);
NcDim *NumSnap_dim = f.add_dim("numsnap_dim", num_snapshots);
@ -446,6 +452,13 @@ void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& a
boxed->time = simu->time;
boxed->BoxSize = simu->BoxSize;
NcAtt *d_sub = f.get_att("data_subsampling");
if (d_sub == 0 || d_sub->as_double(0) != args_info.subsample_arg)
{
cerr << "Parameter file was not generated with the same simulation subsampling argument. Particles will be different. Stop here." << endl;
exit(1);
}
NcVar *v_id = f.get_var("particle_ids");
NcVar *v_snap = f.get_var("snapshot_split");
long *edges1;
@ -710,7 +723,7 @@ int main(int argc, char **argv)
delete[] efac;
}
saveBox(simuOut, args_info.outputParameter_arg);
saveBox(simuOut, args_info.outputParameter_arg, args_info);
generateOutput(simuOut, args_info.axis_arg,
args_info.output_arg);
delete preselector;

View file

@ -190,7 +190,7 @@ public:
if (load_flags & (NEED_POSITION | NEED_VELOCITY))
rescaleParticles(d, d->Hubble*1e-5, one_kpc/one_Gyr);
// enforceBoxSize(d);
enforceBoxSize(d);
applyTransformations(d);
basicPreprocessing(d, preproc);

View file

@ -234,6 +234,7 @@ ExternalProject_Add(cosmotool
-DGSLCBLAS_LIBRARY=${GSLCBLAS_LIBRARY}
-DNETCDF_LIBRARY=${NETCDF_LIBRARY}
-DNETCDFCPP_LIBRARY=${NETCDFCPP_LIBRARY}
-DENABLE_SHARP=OFF
)
SET(COSMOTOOL_LIBRARY ${CMAKE_BINARY_DIR}/ext_build/cosmotool/lib/libCosmoTool.a)
set(COSMOTOOL_INCLUDE_PATH ${CMAKE_BINARY_DIR}/ext_build/cosmotool/include)

View file

@ -163,6 +163,8 @@ def launchGenerate(sample, binPath, workDir=None, inputDataDir=None,
dataFileLine = "gadget " + datafile
elif sample.dataFormat == "sdf":
dataFileLine = "sdf " + datafile
else:
raise ValueError("unknown dataFormat '%s'" % sample.dataFormat)
iX = float(sample.mySubvolume[0])
iY = float(sample.mySubvolume[1])

View file

@ -85,7 +85,7 @@ int main(int argc,char **argv) {
printf("Bad density threshold.\n");
exit(0);
}
if (sscanf(argv[7],"%f",&mockIndex) == 0) {
if (sscanf(argv[7],"%d",&mockIndex) == 0) {
printf("Bad mock galaxy index.\n");
exit(0);
}
@ -108,29 +108,29 @@ int main(int argc,char **argv) {
p = (PARTICLE *)malloc(np*sizeof(PARTICLE));
/* Adjacencies*/
for (i=0;i<np;i++) {
fread(&p[i].nadj,1,sizeof(int),adj);
fread(&p[i].nadj,1,sizeof(pid_t),adj);
/* The number of adjacencies per particle */
if (p[i].nadj > 0)
p[i].adj = (int *)malloc(p[i].nadj*sizeof(int));
p[i].adj = (pid_t *)malloc(p[i].nadj*sizeof(pid_t));
else p[i].adj = 0;
p[i].ncnt = 0; /* Temporarily, it's an adj counter */
}
for (i=0;i<np;i++) {
fread(&nin,1,sizeof(int),adj);
fread(&nin,1,sizeof(pid_t),adj);
if (nin > 0)
for (k=0;k<nin;k++) {
fread(&j,1,sizeof(int),adj);
fread(&j,1,sizeof(pid_t),adj);
if (j < np) {
/* Set both halves of the pair */
assert(i < j);
if (p[i].ncnt == p[i].nadj)
{
p[i].adj = (int *)realloc(p[i].adj, (p[i].nadj+1)*sizeof(int));
p[i].nadj++;
p[i].adj = (pid_t *)realloc(p[i].adj, (p[i].nadj+1)*sizeof(pid_t));
p[i].nadj++;
}
if (p[j].ncnt == p[j].nadj)
{
p[j].adj = (int *)realloc(p[j].adj, (p[j].nadj+1)*sizeof(int));
p[j].adj = (pid_t *)realloc(p[j].adj, (p[j].nadj+1)*sizeof(pid_t));
p[j].nadj++;
}
p[i].adj[p[i].ncnt] = j;
@ -144,13 +144,13 @@ int main(int argc,char **argv) {
fclose(adj);
/* Check that we got all the pairs */
// adj = fopen(adjfile, "r");
// fread(&np,1, sizeof(int),adj);
/* adj = fopen(adjfile, "r");
fread(&np,1, sizeof(int),adj);*/
for (i=0;i<np;i++) {
// fread(&nin,1,sizeof(int),adj); /* actually nadj */
/* fread(&nin,1,sizeof(int),adj); /* actually nadj */
// PMS
if (p[i].ncnt != p[i].nadj && i < mockIndex) {
//if (p[i].ncnt != p[i].nadj) {
/*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);
@ -185,9 +185,9 @@ int main(int argc,char **argv) {
}
fclose(vol);
jumped = (int *)malloc(np*sizeof(int));
jumper = (int *)malloc(np*sizeof(int));
numinh = (int *)malloc(np*sizeof(int));
jumped = (pid_t *)malloc(np*sizeof(pid_t));
jumper = (pid_t *)malloc(np*sizeof(pid_t));
numinh = (pid_t *)malloc(np*sizeof(pid_t));
/* find jumper */
for (i = 0; i < np; i++) {
@ -252,7 +252,7 @@ int main(int argc,char **argv) {
}
for (h=0;h<nzones;h++) {
zt[h].adj = (int *)malloc(zt[h].nadj*sizeof(int));
zt[h].adj = (pid_t *)malloc(zt[h].nadj*sizeof(pid_t));
if (zt[h].adj == NULL) {
printf("Unable to allocate %d adj's of zone %d\n",zt[h].nadj,h);
exit(0);
@ -313,7 +313,7 @@ int main(int argc,char **argv) {
for (h=0;h<nzones;h++) {
/*printf("%d ",zt[h].nadj);*/
z[h].nadj = zt[h].nadj;
z[h].adj = (int *)malloc(zt[h].nadj*sizeof(int));
z[h].adj = (pid_t *)malloc(zt[h].nadj*sizeof(pid_t));
z[h].slv = (float *)malloc(zt[h].nadj*sizeof(float));
for (za = 0; za<zt[h].nadj; za++) {
z[h].adj[za] = zt[h].adj[za];
@ -326,12 +326,12 @@ int main(int argc,char **argv) {
free(zt);
free(numinh);
m = (int **)malloc(nzones*sizeof(int *));
m = (pid_t **)malloc(nzones*sizeof(pid_t *));
/* Not in the zone struct since it'll be freed up (contiguously, we hope)
soon */
nm = (int *)malloc(nzones*sizeof(int));
nm = (pid_t *)malloc(nzones*sizeof(pid_t));
for (h=0; h<nzones; h++) {
m[h] = (int *)malloc(z[h].np*sizeof(int));
m[h] = (pid_t *)malloc(z[h].np*sizeof(pid_t));
nm[h] = 0;
z[h].vol = 0.;
}
@ -354,14 +354,14 @@ int main(int argc,char **argv) {
printf("Problem opening zonefile %s.\n\n",zonfile);
exit(0);
}
fwrite(&np,1,4,zon);
fwrite(&nzones,1,4,zon);
fwrite(&np,1,sizeof(pid_t),zon);
fwrite(&nzones,1,sizeof(int),zon);
for (h=0; h<nzones; h++) {
// PMS
//printf("%d %d %d", &(z[h].np), m[h], z[h].np);
// END PMS
fwrite(&(z[h].np),1,4,zon);
fwrite(m[h],z[h].np,4,zon);
fwrite(&(z[h].np),1,sizeof(pid_t),zon);
fwrite(m[h],z[h].np,sizeof(pid_t),zon);
free(m[h]);
}
free(m);
@ -394,7 +394,7 @@ int main(int argc,char **argv) {
printf("Problem opening zonefile %s.\n\n",zonfile2);
exit(0);
}
fwrite(&nzones,1,4,zon2);
fwrite(&nzones,1,sizeof(int),zon2);
for (h = 0; h<nzones; h++) {
nhlcount = 0;
@ -532,8 +532,8 @@ int main(int argc,char **argv) {
z[h].nhl = nhl;
fwrite(&nhl,1,4,zon2);
fwrite(zonelist,nhl,4,zon2);
fwrite(&nhl,1,sizeof(int),zon2);
fwrite(zonelist,nhl,sizeof(int),zon2);
}
fclose(zon2);

View file

@ -3,7 +3,9 @@
##define NGUARD 42 /*Actually, the number of SPACES between guard points
in each dim */
typedef int pid_t;
typedef struct Partadj {
int nadj;
int *adj;
pid_t nadj;
pid_t *adj;
} PARTADJ;

View file

@ -12,7 +12,7 @@ int posread(char *posfile, float ***p, float fact);
int main(int argc, char *argv[]) {
int exitcode;
int i, j, np;
pid_t i, j, np;
float **r;
coordT rtemp[3], *parts;
coordT deladjs[3*MAXVERVER], points[3*MAXVERVER];
@ -22,11 +22,12 @@ int main(int argc, char *argv[]) {
PARTADJ *adjs;
float *vols;
float predict, xmin,xmax,ymin,ymax,zmin,zmax;
int *orig;
pid_t *orig;
int isitinbuf;
char isitinmain, d;
int numdiv, nvp, nvpall, nvpbuf;
int numdiv;
pid_t nvp, nvpall, nvpbuf;
float width, width2, totwidth, totwidth2, bf, s, g;
float border, boxsize;
float c[3];
@ -142,7 +143,7 @@ int main(int argc, char *argv[]) {
points */
parts = (coordT *)malloc(3*nvpbuf*sizeof(coordT));
orig = (int *)malloc(nvpbuf*sizeof(int));
orig = (pid_t *)malloc(nvpbuf*sizeof(pid_t));
if (parts == NULL) {
printf("Unable to allocate parts\n");
@ -182,15 +183,16 @@ int main(int argc, char *argv[]) {
nvpbuf = nvp;
for (i=0; i<np; i++) {
isitinbuf = 1;
isitinmain = 1;
DL {
rtemp[d] = r[i][d] - c[d];
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 > 0) &&
((fabs(rtemp[0])>width2)||(fabs(rtemp[1])>width2)||(fabs(rtemp[2])>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);*/
@ -323,14 +325,14 @@ int main(int argc, char *argv[]) {
printf("nvp = %d\n",nvp);
/* Tell us where the original particles were */
fwrite(orig,sizeof(int),nvp,out);
fwrite(orig,sizeof(pid_t),nvp,out);
/* Volumes*/
fwrite(vols,sizeof(float),nvp,out);
/* Adjacencies */
for (i=0;i<nvp;i++) {
fwrite(&(adjs[i].nadj),1,sizeof(int),out);
fwrite(&(adjs[i].nadj),1,sizeof(pid_t),out);
if (adjs[i].nadj > 0)
fwrite(adjs[i].adj,adjs[i].nadj,sizeof(int),out);
fwrite(adjs[i].adj,adjs[i].nadj,sizeof(pid_t),out);
else printf("0");
}
fclose(out);

View file

@ -14,10 +14,12 @@ int main(int argc, char *argv[]) {
int numdiv,np,np2,na;
int i,j,k,p,nout;
int nvp,npnotdone,nvpmax, nvpsum, *orig;
pid_t i,j,k,p,nout;
pid_t nvp,npnotdone,nvpmax, nvpsum, *orig;
double avgnadj, avgvol;
int numRemoved = 0;
// PMS
int mockIndex;
// END PMS
@ -87,7 +89,7 @@ int main(int argc, char *argv[]) {
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));
orig = (pid_t *)malloc(nvpmax*sizeof(pid_t));
if (orig == NULL) printf("Couldn't allocate orig.\n");
if ((cnt_adj==NULL) || (vols == NULL) || (orig == NULL) || (adjs == NULL)) {
printf("Not enough memory to allocate. Exiting.\n");
@ -96,6 +98,7 @@ int main(int argc, char *argv[]) {
for (p=0;p<np;p++)
vols[p] = -1.;
nvpsum = 0;
for (i = 0; i < numdiv; i++) {
for (j = 0; j < numdiv; j++) {
for (k = 0; k < numdiv; k++) {
@ -105,13 +108,13 @@ int main(int argc, char *argv[]) {
printf("Unable to open file %s.\n\n",partfile);
exit(0);
}
fread(&np2,1,sizeof(int),part);
fread(&nvp,1,sizeof(int),part);
fread(&np2,1,sizeof(pid_t),part);
fread(&nvp,1,sizeof(pid_t),part);
/*printf("nvp = %d\n",nvp);fflush(stdout);*/
nvpsum += nvp;
fread(orig,nvp,sizeof(int),part);
fread(orig,nvp,sizeof(pid_t),part);
for (p=0;p<nvp;p++) {
fread(&volstemp,1,sizeof(float),part);
if (vols[orig[p]] > -1.)
@ -131,12 +134,12 @@ int main(int argc, char *argv[]) {
fread(&na,1,sizeof(int),part);
if (na > 0) {
adjs[orig[p]].nadj = na;
adjs[orig[p]].adj = (int *)malloc(na*sizeof(int));
adjs[orig[p]].adj = (pid_t *)malloc(na*sizeof(pid_t));
if (adjs[orig[p]].adj == NULL) {
printf("Couldn't allocate adjs[orig[%d]].adj.\n",p);
exit(0);
}
fread(adjs[orig[p]].adj,na,sizeof(int),part);
fread(adjs[orig[p]].adj,na,sizeof(pid_t),part);
} else {
printf("0"); fflush(stdout);
}
@ -157,8 +160,6 @@ int main(int argc, char *argv[]) {
adjs[i].nadj = 0;
}
int numRemoved = 0;
// unlink particles adjacent to mock galaxies
for (i = 0; i < mockIndex; i++) {
for (j = 0; j < adjs[i].nadj; j++) {
@ -217,7 +218,7 @@ printf("\n");
avgvol += (double)(vols[p]);
}
if (npnotdone > 0)
printf("%d particles not done!\n");
printf("%d particles not done!\n", npnotdone);
printf("%d particles done more than once.\n",nvpsum-np);
avgnadj /= (double)np;
avgvol /= (double)np;
@ -272,9 +273,9 @@ printf("\n");
for (j=0;j<adjs[i].nadj; j++) if (adjs[i].adj[j] > i) nout++;
fwrite(&nout,1,sizeof(int),adj);
for (j=0;j<adjs[i].nadj; j++) {
int id = adjs[i].adj[j];
pid_t id = adjs[i].adj[j];
if (id > i)
fwrite(&(adjs[i].adj[j]),1,sizeof(int),adj);
fwrite(&id,1,sizeof(pid_t),adj);
}
}