mirror of
https://bitbucket.org/cosmicvoids/vide_public.git
synced 2025-07-04 15:21:11 +00:00
Large rework of generateMock for handling partial load of snapshots and thus larger simulations. Work not yet complete.
This commit is contained in:
Guilhem Lavaux
parent
60d77d3ae0
commit
26f6864f47
1 changed files with 168 additions and 131 deletions
|
|
@ -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<long>);
|
||||
|
||||
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<long>& 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<long> 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<long>("uniqueID");
|
||||
boxed->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
|
||||
|
||||
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<long>("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<long>);
|
||||
boxed->new_attribute("expansion_fac", expansion_fac, delete_adaptor<double>);
|
||||
boxed->new_attribute("uniqueID", uniqueID, delete_adaptor<long>);
|
||||
boxed->new_attribute("mul", mul, delete_adaptor<double>);
|
||||
boxed->new_attribute("ranges", ranges, delete_adaptor<double>);
|
||||
boxed->new_attribute("snapshot_split", snap_split, delete_adaptor<long>);
|
||||
boxed->new_attribute("num_snapshots", numsnap_info, delete_adaptor<int>);
|
||||
}
|
||||
|
||||
void buildBox(SimuData *simu, const std::vector<long>& targets, long& loaded,
|
||||
SimuData *boxed, double *efac)
|
||||
{
|
||||
uint32_t k = 0;
|
||||
for (uint32_t i = 0; i < simu->NumPart; i++)
|
||||
long *uniqueID = boxed->as<long>("uniqueID");
|
||||
long *simu_uniqueID = simu->as<long>("uniqueID");
|
||||
double *expansion_fac = boxed->as<double>("expansion_fac");
|
||||
double *mul = boxed->as<double>("mul");
|
||||
double *ranges = boxed->as<double>("ranges");
|
||||
long *particle_id = boxed->as<long>("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<double>("ranges");
|
||||
NcFile f(args_info.outputParameter_arg, NcFile::Replace);
|
||||
long *particle_id = boxed->as<long>("particle_id");
|
||||
double *expansion_fac = boxed->as<double>("expansion_fac");
|
||||
long *snapshot_split = boxed->as<long>("snapshot_split");
|
||||
int num_snapshots = *boxed->as<int>("num_snapshots");
|
||||
long *uniqueID = boxed->as<long>("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<long>);
|
||||
boxed->new_attribute("mul", mul, delete_adaptor<double>);
|
||||
boxed->new_attribute("ranges", ranges, delete_adaptor<double>);
|
||||
boxed->new_attribute("particle_id", particle_id, delete_adaptor<long>);
|
||||
boxed->new_attribute("num_snapshots", num_snapshots, delete_adaptor<int>);
|
||||
boxed->new_attribute("snapshot_split", snapshot_split, delete_adaptor<long>);
|
||||
|
||||
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<long> 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;
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue