From 26f6864f47d7c7838ab73ad428690d682129fa55 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Thu, 22 Nov 2012 10:51:50 -0500 Subject: [PATCH] Large rework of generateMock for handling partial load of snapshots and thus larger simulations. Work not yet complete. --- c_tools/mock/generateMock.cpp | 299 +++++++++++++++++++--------------- 1 file changed, 168 insertions(+), 131 deletions(-) diff --git a/c_tools/mock/generateMock.cpp b/c_tools/mock/generateMock.cpp index a6ac7dc..49d25e4 100644 --- a/c_tools/mock/generateMock.cpp +++ b/c_tools/mock/generateMock.cpp @@ -79,12 +79,12 @@ SimuData *doLoadRamses(const char *basename, int baseid, int velAxis, bool goRed SimuData *myLoadGadget(const char *fname, int id, int flags) { SimuData *sim = loadGadgetMulti(fname, id, flags); - sim->BoxSize *= gadgetUnit*1000; + sim->BoxSize *= gadgetUnit; for (int j = 0; j < 3; j++) { if (sim->Pos[j] != 0) { for (long i = 0; i < sim->NumPart; i++) - sim->Pos[j][i] *= gadgetUnit*1000; + sim->Pos[j][i] *= gadgetUnit; } } return sim; @@ -94,6 +94,7 @@ SimuData *doLoadSimulation(const char *gadgetname, int velAxis, bool goRedshift, { SimuData *d, *outd; bool singleFile = false; + double redshift_gravity; try { @@ -114,7 +115,7 @@ SimuData *doLoadSimulation(const char *gadgetname, int velAxis, bool goRedshift, outd = new SimuData; outd->NumPart = d->TotalNumPart; - outd->BoxSize = d->BoxSize/1000; + outd->BoxSize = d->BoxSize; outd->TotalNumPart = outd->NumPart; outd->Hubble = d->Hubble; outd->Omega_Lambda = d->Omega_Lambda; @@ -127,31 +128,34 @@ SimuData *doLoadSimulation(const char *gadgetname, int velAxis, bool goRedshift, outd->Id = new int[outd->NumPart]; delete d; + long *uniqueID = new long[outd->NumPart]; + + outd->new_attribute("uniqueID", uniqueID, delete_adaptor); + + redshift_gravity = goRedshift ? 1 : 0; + int curCpu = singleFile ? -1 : 0; cout << "loading file 0 " << endl; try { + long loaded = 0; while (1) { d = loadFunction(gadgetname, curCpu, NEED_POSITION|NEED_VELOCITY|NEED_GADGET_ID); - for (int k = 0; k < 3; k++) - for (int i = 0; i < d->NumPart; i++) - { - int pid = d->Id[i]-1; - assert(pid >= 0); - assert(pid < outd->TotalNumPart); - outd->Pos[k][pid] = d->Pos[k][i]/1000; - outd->Vel[2][pid] = d->Vel[velAxis][i]; - outd->Id[pid] = pid+1; - } + for (int i = 0; i < d->NumPart; i++) + { + for (int k = 0; k < 3; k++) + outd->Pos[k][loaded] = d->Pos[k][i]; + outd->Vel[2][loaded] = d->Vel[velAxis][i]; + outd->Id[loaded] = loaded; + uniqueID[loaded] = d->Id[i]; + assert(loaded < outd->TotalNumPart); + + outd->Pos[velAxis][loaded] += + redshift_gravity*d->Vel[velAxis][i]/100.; + loaded++; + } - if (goRedshift) - for (int i = 0; i < d->NumPart; i++) - { - int pid = d->Id[i]-1; - outd->Pos[velAxis][pid] += d->Vel[velAxis][i]/100.; - } - delete d; if (singleFile) break; @@ -425,30 +429,16 @@ void generateOutput(SimuData *data, int axis, } } -void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& args_info) +// This function prepares the list of targets for the specified snapshot. The target list is not +// cleared. That way new particles can be appended if this is a multi-file snapshot. +void selectBox(SimuData *simu, std::vector& targets, generateMock_info& args_info) { float subsample = args_info.subsample_given ? args_info.subsample_arg : 1.0; - uint32_t goodParticles = 0; double ranges[3][2] = { { args_info.rangeX_min_arg, args_info.rangeX_max_arg }, { args_info.rangeY_min_arg, args_info.rangeY_max_arg }, { args_info.rangeZ_min_arg, args_info.rangeZ_max_arg } }; - double mul[3]; - float minmax[2][3]; - int *particle_id; - bool *random_acceptance = 0; - - boxed = new SimuData; - boxed->Hubble = simu->Hubble; - boxed->Omega_M = simu->Omega_M; - boxed->Omega_Lambda = simu->Omega_Lambda; - boxed->time = simu->time; - boxed->BoxSize = simu->BoxSize; - - random_acceptance = new bool[simu->NumPart]; - - for (int j = 0; j < 3; j++) minmax[1][j] = minmax[0][j] = simu->Pos[j][0]; for (uint32_t i = 0; i < simu->NumPart; i++) { @@ -458,103 +448,140 @@ void makeBox(SimuData *simu, double *efac, SimuData *&boxed, generateMock_info& acceptance = acceptance && (simu->Pos[j][i] > ranges[j][0]) && - (simu->Pos[j][i] < ranges[j][1]); - minmax[0][j] = min(simu->Pos[j][i], minmax[0][j]); - minmax[1][j] = max(simu->Pos[j][i], minmax[1][j]); + (simu->Pos[j][i] < ranges[j][1]); } - random_acceptance[i] = acceptance && (drand48() <= subsample); - if (random_acceptance[i]) - goodParticles++; + + if (acceptance && (drand48() <= subsample)) + targets.push_back(i); } - cout << "Subsample fraction: " << subsample << endl; - cout << "Min range = " << ranges[0][0] << " " << ranges[1][0] << " " << ranges[2][0] << endl; - cout << "Max range = " << ranges[0][1] << " " << ranges[1][1] << " " << ranges[2][1] << endl; +} - cout << "Min position = " << minmax[0][0] << " " << minmax[0][1] << " " << minmax[0][2] << endl; - cout << "Max position = " << minmax[1][0] << " " << minmax[1][1] << " " << minmax[1][2] << endl; +void createBox(SimuData *simu, long num_targets, std::vector snapshot_split, SimuData *& boxed, generateMock_info& args_info) +{ + double *ranges = new double[6]; + double *mul = new double[3]; - cout << "Number of accepted particles: " << goodParticles << endl; + ranges[0] = args_info.rangeX_min_arg; + ranges[1] = args_info.rangeX_max_arg; + ranges[2] = args_info.rangeY_min_arg; + ranges[3] = args_info.rangeY_max_arg; + ranges[4] = args_info.rangeZ_min_arg; + ranges[5] = args_info.rangeZ_max_arg; + + boxed = new SimuData; + boxed->Hubble = simu->Hubble; + boxed->Omega_M = simu->Omega_M; + boxed->Omega_Lambda = simu->Omega_Lambda; + boxed->time = simu->time; + boxed->BoxSize = simu->BoxSize; for (int j = 0; j < 3; j++) { - boxed->Pos[j] = new float[goodParticles]; + boxed->Pos[j] = new float[num_targets]; boxed->Vel[j] = 0; - mul[j] = 1.0/(ranges[j][1] - ranges[j][0]); + mul[j] = 1.0/(ranges[2*j+1] - ranges[2*j+0]); } - long *uniqueID = new long[goodParticles]; - long *simu_uniqueID = simu->as("uniqueID"); - boxed->new_attribute("uniqueID", uniqueID, delete_adaptor); - + cout << "Min range = " << ranges[0] << " " << ranges[2] << " " << ranges[4] << endl; + cout << "Max range = " << ranges[1] << " " << ranges[3] << " " << ranges[5] << endl; + cout << "Number of accepted particles: " << num_targets << endl; cout << "Rescaling factors = " << mul[0] << " " << mul[1] << " " << mul[2] << endl; - boxed->NumPart = goodParticles; - particle_id = new int[goodParticles]; - double *expansion_fac = new double[goodParticles]; + long *uniqueID = new long[num_targets]; + long *simu_uniqueID = simu->as("uniqueID"); + long *particle_id = new long[num_targets]; + double *expansion_fac = new double[num_targets]; + long *snap_split = new long[snapshot_split.size()]; + int *numsnap_info = new int[1]; + copy(snapshot_split.begin(), snapshot_split.end(), snap_split); + *numsnap_info = snapshot_split.size(); + + boxed->NumPart = num_targets; + boxed->new_attribute("particle_id", particle_id, delete_adaptor); + boxed->new_attribute("expansion_fac", expansion_fac, delete_adaptor); + boxed->new_attribute("uniqueID", uniqueID, delete_adaptor); + boxed->new_attribute("mul", mul, delete_adaptor); + boxed->new_attribute("ranges", ranges, delete_adaptor); + boxed->new_attribute("snapshot_split", snap_split, delete_adaptor); + boxed->new_attribute("num_snapshots", numsnap_info, delete_adaptor); +} + +void buildBox(SimuData *simu, const std::vector& targets, long& loaded, + SimuData *boxed, double *efac) +{ uint32_t k = 0; - for (uint32_t i = 0; i < simu->NumPart; i++) + long *uniqueID = boxed->as("uniqueID"); + long *simu_uniqueID = simu->as("uniqueID"); + double *expansion_fac = boxed->as("expansion_fac"); + double *mul = boxed->as("mul"); + double *ranges = boxed->as("ranges"); + long *particle_id = boxed->as("particle_id"); + + for (uint32_t i = 0; i < targets.size(); i++) { - bool acceptance = random_acceptance[i]; + long pid = targets[i]; - if (acceptance) + for (int j = 0; j < 3; j++) { - for (int j = 0; j < 3; j++) - { - boxed->Pos[j][k] = (simu->Pos[j][i]-ranges[j][0])*mul[j]; - assert(boxed->Pos[j][k] > 0); - assert(boxed->Pos[j][k] < 1); - } - particle_id[k] = simu->Id[i]-1; - uniqueID[k] = (simu_uniqueID != 0) ? simu_uniqueID[i] : 0; - expansion_fac[k] = efac[i]; - k++; + boxed->Pos[j][loaded] = (simu->Pos[j][pid]-ranges[j*2])*mul[j]; + assert(boxed->Pos[j][loaded] > 0); + assert(boxed->Pos[j][loaded] < 1); } + particle_id[loaded] = pid; + uniqueID[loaded] = (simu_uniqueID != 0) ? simu_uniqueID[i] : 0; + expansion_fac[loaded] = efac[i]; + loaded++; + k++; } +} - delete[] random_acceptance; - +void saveBox(SimuData *&boxed, generateMock_info& args_info) +{ + double *ranges = boxed->as("ranges"); NcFile f(args_info.outputParameter_arg, NcFile::Replace); + long *particle_id = boxed->as("particle_id"); + double *expansion_fac = boxed->as("expansion_fac"); + long *snapshot_split = boxed->as("snapshot_split"); + int num_snapshots = *boxed->as("num_snapshots"); + long *uniqueID = boxed->as("uniqueID"); - f.add_att("range_x_min", ranges[0][0]); - f.add_att("range_x_max", ranges[0][1]); - f.add_att("range_y_min", ranges[1][0]); - f.add_att("range_y_max", ranges[1][1]); - f.add_att("range_z_min", ranges[2][0]); - f.add_att("range_z_max", ranges[2][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]); + f.add_att("range_y_max", ranges[3]); + f.add_att("range_z_min", ranges[4]); + f.add_att("range_z_max", ranges[5]); f.add_att("mask_index", -1); NcDim *NumPart_dim = f.add_dim("numpart_dim", boxed->NumPart); - NcVar *v = f.add_var("particle_ids", ncInt, NumPart_dim); + NcDim *NumSnap_dim = f.add_dim("numsnap_dim", num_snapshots); + NcVar *v = f.add_var("particle_ids", ncLong, NumPart_dim); NcVar *v2 = f.add_var("expansion", ncDouble, NumPart_dim); + NcVar *v3 = f.add_var("snapshot_split", ncLong, NumSnap_dim); v->put(particle_id, boxed->NumPart); - v2->put(expansion_fac, boxed->NumPart); + v2->put(expansion_fac, boxed->NumPart); + v3->put(snapshot_split, num_snapshots); + if (uniqueID != 0) + { + NcVar *v4 = f.add_var("unique_ids", ncLong, NumPart_dim); + v4->put(uniqueID, boxed->NumPart); + } delete[] particle_id; delete[] expansion_fac; - - -/* - FILE *fp = fopen("sample_info.txt", "w"); - fprintf(fp, "x_min = %f\n", ranges[0][0]); - fprintf(fp, "x_max = %f\n", ranges[0][1]); - fprintf(fp, "y_min = %f\n", ranges[1][0]); - fprintf(fp, "y_max = %f\n", ranges[1][1]); - fprintf(fp, "z_min = %f\n", ranges[2][0]); - fprintf(fp, "z_max = %f\n", ranges[2][1]); - fprintf(fp, "mask_index = -1\n"); - fprintf(fp, "total_particles = %d\n", boxed->NumPart); - fclose(fp); -*/ } -void makeBoxFromParameter(SimuData *simu, double *efac, SimuData* &boxed, generateMock_info& args_info) +void makeBoxFromParameter(SimuData *simu, SimuData* &boxed, generateMock_info& args_info) { NcFile f(args_info.inputParameter_arg); NcVar *v; - int *particle_id; + long *particle_id; double *expansion_fac; + long *uniqueID; + long *snapshot_split; + int *num_snapshots = new int[1]; boxed = new SimuData; boxed->Hubble = simu->Hubble; @@ -564,55 +591,56 @@ void makeBoxFromParameter(SimuData *simu, double *efac, SimuData* &boxed, genera boxed->BoxSize = simu->BoxSize; NcVar *v_id = f.get_var("particle_ids"); + NcVar *v_snap = f.get_var("snapshot_split"); long *edges1; - double ranges[3][2]; - double mul[3]; + long *dim_snap; + double *ranges; + double *mul; edges1 = v_id->edges(); + dim_snap = v_snap->edges(); assert(v_id->num_dims()==1); + assert(v_snap->num_dims()==1); boxed->NumPart = edges1[0]; + *num_snapshots = dim_snap[0]; + delete[] dim_snap; delete[] edges1; - particle_id = new int[boxed->NumPart]; + particle_id = new long[boxed->NumPart]; + uniqueID = new long[boxed->NumPart]; + mul = new double[3]; + ranges = new double[6]; + snapshot_split = new long[*num_snapshots]; + + + boxed->new_attribute("uniqueID", uniqueID, delete_adaptor); + boxed->new_attribute("mul", mul, delete_adaptor); + boxed->new_attribute("ranges", ranges, delete_adaptor); + boxed->new_attribute("particle_id", particle_id, delete_adaptor); + boxed->new_attribute("num_snapshots", num_snapshots, delete_adaptor); + boxed->new_attribute("snapshot_split", snapshot_split, delete_adaptor); v_id->get(particle_id, boxed->NumPart); + v_snap->get(snapshot_split, *num_snapshots); - ranges[0][0] = f.get_att("range_x_min")->as_double(0); - ranges[0][1] = f.get_att("range_x_max")->as_double(0); - ranges[1][0] = f.get_att("range_y_min")->as_double(0); - ranges[1][1] = f.get_att("range_y_max")->as_double(0); - ranges[2][0] = f.get_att("range_z_min")->as_double(0); - ranges[2][1] = f.get_att("range_z_max")->as_double(0); + ranges[0] = f.get_att("range_x_min")->as_double(0); + ranges[1] = f.get_att("range_x_max")->as_double(0); + ranges[2] = f.get_att("range_y_min")->as_double(0); + ranges[3] = f.get_att("range_y_max")->as_double(0); + ranges[4] = f.get_att("range_z_min")->as_double(0); + ranges[5] = f.get_att("range_z_max")->as_double(0); for (int j = 0; j < 3; j++) { boxed->Pos[j] = new float[boxed->NumPart]; boxed->Vel[j] = 0; - mul[j] = 1.0/(ranges[j][1] - ranges[j][0]); + mul[j] = 1.0/(ranges[2*j+1] - ranges[2*j+0]); } uint32_t k = 0; - uint32_t max_id = *std::max_element(simu->Id, simu->Id + simu->NumPart); - uint32_t min_id = *std::min_element(simu->Id, simu->Id + simu->NumPart); - uint32_t *id_to_part = new uint32_t[max_id-min_id+1]; - for (uint32_t i = 0; i < simu->NumPart; i++) - id_to_part[simu->Id[i]-min_id] = i; - - for (uint32_t i = 0; i < boxed->NumPart; i++) - { - int id = particle_id[i]; - assert(id >= min_id && id <= max_id); - int pid = id_to_part[id-min_id]; - - for (int j = 0; j < 3; j++) - { - boxed->Pos[j][i] = (simu->Pos[j][pid]-ranges[j][0])*mul[j]; - } - } - - delete[] id_to_part; - delete[] particle_id; + NcVar *v_uniq = f.get_var("unique_ids"); + v_uniq->get(uniqueID, boxed->NumPart); } int main(int argc, char **argv) @@ -707,13 +735,22 @@ int main(int argc, char **argv) args_info.peculiarVelocities_flag, expfact, args_info.cosmo_flag); - if (args_info.inputParameter_given) - makeBoxFromParameter(simu, expfact, simuOut, args_info); - else - makeBox(simu, expfact, simuOut, args_info); + vector targets, split; + long loaded = 0; + if (args_info.inputParameter_given) + makeBoxFromParameter(simu, simuOut, args_info); + else + { + selectBox(simu, targets, args_info); + split.push_back(targets.size()); + createBox(simu, targets.size(), split, simuOut, args_info); + } + + buildBox(simu, targets, loaded, simuOut, expfact); delete simu; + saveBox(simuOut, args_info); generateOutput(simuOut, args_info.axis_arg, args_info.output_arg); delete simuOut;