From 65eb5c4af27b355d1eccfc8d8248920cdeb75718 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Thu, 18 Apr 2013 16:09:08 +0200 Subject: [PATCH 1/9] Revert "Indicate at which point of the text stream is the failure" This reverts commit 4f9fd340f0f76aedba431491118addd1372ecb7c. --- c_tools/libzobov/loadZobov.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/c_tools/libzobov/loadZobov.cpp b/c_tools/libzobov/loadZobov.cpp index e34ef45..a1a6cbb 100644 --- a/c_tools/libzobov/loadZobov.cpp +++ b/c_tools/libzobov/loadZobov.cpp @@ -134,15 +134,11 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName, } cout << "Loading description" << endl; - int lineid = 1; string line; getline(descFile, line); - lineid++; getline(descFile, line); - lineid++; getline(descFile, line); - lineid++; while (!descFile.eof()) { istringstream lineStream(line.c_str()); @@ -164,7 +160,7 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName, >> probability; if (!lineStream) { - cerr << "Error in text stream at line " << lineid << endl; + cerr << "Error in text stream" << endl; abort(); } @@ -191,7 +187,6 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName, << actualNumber << ")" << endl; } getline(descFile, line); - lineid++; } cout << "Done loading" << endl; From d05d5f241b86f897295746fe60deff71630a4152 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Thu, 18 Apr 2013 16:09:29 +0200 Subject: [PATCH 2/9] Revert "Revert "Indicate at which point of the text stream is the failure"" This reverts commit 65eb5c4af27b355d1eccfc8d8248920cdeb75718. --- c_tools/libzobov/loadZobov.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/c_tools/libzobov/loadZobov.cpp b/c_tools/libzobov/loadZobov.cpp index a1a6cbb..e34ef45 100644 --- a/c_tools/libzobov/loadZobov.cpp +++ b/c_tools/libzobov/loadZobov.cpp @@ -134,11 +134,15 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName, } cout << "Loading description" << endl; + int lineid = 1; string line; getline(descFile, line); + lineid++; getline(descFile, line); + lineid++; getline(descFile, line); + lineid++; while (!descFile.eof()) { istringstream lineStream(line.c_str()); @@ -160,7 +164,7 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName, >> probability; if (!lineStream) { - cerr << "Error in text stream" << endl; + cerr << "Error in text stream at line " << lineid << endl; abort(); } @@ -187,6 +191,7 @@ bool loadZobov(const char *descName, const char *adjName, const char *voidsName, << actualNumber << ")" << endl; } getline(descFile, line); + lineid++; } cout << "Done loading" << endl; From 15d6825f5d8a8f7fd941f61cba0334351aaeb11a Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 10 Jun 2013 03:48:56 -0400 Subject: [PATCH 3/9] Preliminary reimplementation of voz1b1 in C++ with more modularity --- .gitignore | 4 - c_tools/zobov2/voz1b1/voz.h | 29 +++ c_tools/zobov2/voz1b1/voz1b1.cpp | 337 +++++++++++++++++++++++++++++++ c_tools/zobov2/voz1b1/voz_io.cpp | 64 ++++++ c_tools/zobov2/voz1b1/voz_io.hpp | 22 ++ c_tools/zobov2/voz1b1/vozutil.c | 245 ++++++++++++++++++++++ 6 files changed, 697 insertions(+), 4 deletions(-) create mode 100644 c_tools/zobov2/voz1b1/voz.h create mode 100644 c_tools/zobov2/voz1b1/voz1b1.cpp create mode 100644 c_tools/zobov2/voz1b1/voz_io.cpp create mode 100644 c_tools/zobov2/voz1b1/voz_io.hpp create mode 100644 c_tools/zobov2/voz1b1/vozutil.c diff --git a/.gitignore b/.gitignore index 487c453..c0ba1bd 100644 --- a/.gitignore +++ b/.gitignore @@ -1,8 +1,4 @@ *~ *.o *.prog -voz1b1 -jozov -vozinit -voztie .mydepend diff --git a/c_tools/zobov2/voz1b1/voz.h b/c_tools/zobov2/voz1b1/voz.h new file mode 100644 index 0000000..1d48e92 --- /dev/null +++ b/c_tools/zobov2/voz1b1/voz.h @@ -0,0 +1,29 @@ +#ifndef __VOZ_H +#define __VOZ_H + +#define MAXVERVER 100000 +#define NGUARD 84 /*Actually, the number of SPACES between guard points +##define NGUARD 42 /*Actually, the number of SPACES between guard points + in each dim */ + +typedef int pid_t; + +typedef struct Partadj { + pid_t nadj; + pid_t *adj; +} PARTADJ; + +#ifdef __cplusplus +extern "C" { +#endif + +int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs); +int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol); +int posread(char *posfile, float ***p, float fact); + +#ifdef __cplusplus +} +#endif + + +#endif diff --git a/c_tools/zobov2/voz1b1/voz1b1.cpp b/c_tools/zobov2/voz1b1/voz1b1.cpp new file mode 100644 index 0000000..7912082 --- /dev/null +++ b/c_tools/zobov2/voz1b1/voz1b1.cpp @@ -0,0 +1,337 @@ +#include +#include +#include +#include "libqhull/qhull_a.h" +#include "voz.h" +#include "voz_io.hpp" + +using namespace std; + +#define DL for (d=0;d<3;d++) + +bool checkParameters(int *numdiv, int *b) +{ + for (int i = 0; i < 3; i++) + { + if (numdiv[i] < 0 || b[i] < 0 || b[i] >= numdiv[i]) + return false; + } + return true; +} + + +int main(int argc, char *argv[]) { + const float BF = std::numeric_limits::max(); + + int exitcode; + pid_t i, j, np; + PositionData pdata; + coordT rtemp[3], *parts; + coordT deladjs[3*MAXVERVER], points[3*MAXVERVER]; + pointT intpoints[3*MAXVERVER]; + FILE *pos, *out; + char *posfile, outfile[200], *suffix, *outDir; + PARTADJ *adjs; + float *vols; + realT predict, xmin,xmax,ymin,ymax,zmin,zmax; + pid_t *orig; + + int isitinbuf; + char isitinmain, d; + pid_t nvp, nvpall, nvpbuf; + realT width, width2, totwidth, totwidth2, bf, s, g; + double border, boxsize[3]; + realT c[3]; + int b[3]; + int numdiv[3]; + double totalvol; + + CosmoTool::MiniArgDesc args[] = { + { "POSITION FILE", MINIARG_STRING, &posfile }, + { "BORDER SIZE", MINIARG_DOUBLE, &border }, + { "BOX_X", MINIARG_DOUBLE, &boxsize[0] }, + { "BOX_Y", MINIARG_DOUBLE, &boxsize[1] }, + { "BOX_Z", MINIARG_DOUBLE, &boxsize[2] }, + { "SUFFIX", MINIARG_STRING, &suffix }, + { "NUM_DIVISION_X", MINIARG_INT, &numdiv[0] }, + { "NUM_DIVISION_Y", MINIARG_INT, &numdiv[1] }, + { "NUM_DIVISION_Z", MINIARG_INT, &numdiv[2] }, + { "B0", MINIARG_INT, &b[0] }, + { "B1", MINIARG_INT, &b[1] }, + { "B2", MINIARG_INT, &b[2] }, + { "OUTPUT DIRECTORY", MINIARG_STRING, &outDir }, + { 0, MINIARG_NULL, 0 } + }; + + if (!CosmoTool::parseMiniArgs(argc, argv, args)) + return 1; + + if (!checkParameters(numdiv, b)) + return 2; + + /* Boxsize should be the range in r, yielding a range 0-1 */ + if (!pdata.readFrom(posfile)) + return 3; + + (cout << pdata.np << " particles" << endl).flush(); + + pdata.findExtrema(); + + (cout << boost::format("np: %d, x: %f,%f; y: %f,%f; z: %f,%f") + % pdata.np + % xyz_min[0] % xyz_max[0] + % xyz_min[1] % xyz_max[1] + % xyz_min[2] % xyz_max[2]).flush(); + + + width = 1./(float)numdiv; + width2 = 0.5*width; + if (border > 0.) bf = border; + else bf = 0.1; + /* In units of 0-1, the thickness of each subregion's buffer*/ + totwidth = width+2.*bf; + totwidth2 = width2 + bf; + + s = width/(float)NGUARD; + if ((bf*bf - 2.*s*s) < 0.) { + printf("bf = %f, s = %f.\n",bf,s); + printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.", + sqrt(2.)*width/bf); + exit(0); + } + g = (bf / 2.)*(1. + sqrt(1 - 2.*s*s/(bf*bf))); + printf("s = %f, bf = %f, g = %f.\n",s,bf,g); + + fflush(stdout); + + adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ)); + if (adjs == NULL) { + printf("Unable to allocate adjs\n"); + exit(0); + } + + DL c[d] = ((float)b[d])*width; + printf("c: %f,%f,%f\n",c[0],c[1],c[2]); + /* Assign temporary array*/ + nvpbuf = 0; /* Number of particles to tesselate, including + buffer */ + nvp = 0; /* Without the buffer */ + for (i=0; i 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 (isitinmain) nvp++; + } + + 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)); + + if (parts == NULL) { + printf("Unable to allocate parts\n"); + fflush(stdout); + } + if (orig == NULL) { + printf("Unable to allocate orig\n"); + fflush(stdout); + } + + 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 0.5) rtemp[d] --; + if (rtemp[d] < -0.5) rtemp[d] ++; + isitinmain = isitinmain && (fabs(rtemp[d]) <= width2); + } + if (isitinmain) { + parts[3*nvp] = rtemp[0]; + parts[3*nvp+1] = rtemp[1]; + parts[3*nvp+2] = rtemp[2]; + orig[nvp] = i; + nvp++; + if (rtemp[0] < xmin) xmin = rtemp[0]; + if (rtemp[0] > xmax) xmax = rtemp[0]; + if (rtemp[1] < ymin) ymin = rtemp[1]; + if (rtemp[1] > ymax) ymax = rtemp[1]; + if (rtemp[2] < zmin) zmin = rtemp[2]; + 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); + nvpbuf = nvp; + for (i=0; i 0.5) rtemp[d] --; + if (rtemp[d] < -0.5) rtemp[d] ++; + isitinbuf = isitinbuf && (fabs(rtemp[d]) xmax) xmax = rtemp[0]; + if (rtemp[1] < ymin) ymin = rtemp[1]; + if (rtemp[1] > ymax) ymax = rtemp[1]; + if (rtemp[2] < zmin) zmin = rtemp[2]; + 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); + + pdata.destroy(); + + /* Add guard points */ + for (i=0; i xmax) xmax = parts[3*i]; + if (parts[3*i+1] < ymin) ymin = parts[3*i+1]; + if (parts[3*i+1] > ymax) ymax = parts[3*i+1]; + if (parts[3*i+2] < zmin) zmin = parts[3*i+2]; + if (parts[3*i+2] > zmax) zmax = parts[3*i+2]; + } + + 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); + + /* Do tesselation*/ + printf("File read. Tessellating ...\n"); fflush(stdout); + exitcode = delaunadj(parts, nvp, nvpbuf, nvpall, &adjs); + if (exitcode != 0) + { + printf("Error while tesselating. Stopping here."); fflush(stdout); + exit(1); + } + + /* Now calculate volumes*/ + printf("Now finding volumes ...\n"); fflush(stdout); + vols = (float *)malloc(nvp*sizeof(float)); + + for (i=0; i 0.5) deladjs[3*j+d]--; + } + + 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]); + } + + /* Get the adjacencies back to their original values */ + + for (i=0; i 0) + fwrite(adjs[i].adj,adjs[i].nadj,sizeof(pid_t),out); + else printf("0"); + } + fclose(out); + + return(0); +} diff --git a/c_tools/zobov2/voz1b1/voz_io.cpp b/c_tools/zobov2/voz1b1/voz_io.cpp new file mode 100644 index 0000000..6bf033f --- /dev/null +++ b/c_tools/zobov2/voz1b1/voz_io.cpp @@ -0,0 +1,64 @@ +#include +#include +#include "voz_io.hpp" +#include "CosmoTool/fortran.hpp" +#include + +using namespace CosmoTool; +using namespace std; + +bool PositionData::readFrom(const string& fname) +{ + try + { + UnformattedRead f(fname.c_str()); + + f.beginCheckpoint(); + np = f.readInt32(); + f.endCheckPoint(); + + for (int j = 0; j < 3; j++) + { + xyz[j] = new float[np]; + + f.beginCheckpoint(); + for (int p = 0; p < np; p++) + xyz[j][p] = f.readReal32(); + f.endCheckpoint(); + } + } + catch (const NoSuchFileException& e) + { + return false; + } + catch (const InvalidUnformattedAccess& e) + { + return false; + } + catch (const EndOfFileException& e) + { + return false; + } + + return true; +} + +void PositionData::findExtrema() +{ + const float BF = std::numeric_limits::max(); + + for (int i = 0; i < 3; i++) + { + xyz_min[i] = BF; + xyz_max[i] = -BF; + } + + for (int i = 0; i < 3; i++) + { + for (pid_t p = 0; p < np; p++) + { + xyz_min[p] = min(xyz_min[p], xyz[p][i]); + xyz_max[p] = max(xyz_max[p], xyz[p][i]); + } + } +} diff --git a/c_tools/zobov2/voz1b1/voz_io.hpp b/c_tools/zobov2/voz1b1/voz_io.hpp new file mode 100644 index 0000000..dbeb632 --- /dev/null +++ b/c_tools/zobov2/voz1b1/voz_io.hpp @@ -0,0 +1,22 @@ +#ifndef __VOZ_IO_HPP +#define __VOZ_IO_HPP + +#include + +struct PositionData +{ + float *xyz[3]; + pid_t np; + + void destroy() + { + for (int j = 0; j < 3; j++) + delete[] xyz[j]; + } + + void findExtrema(); + + bool readFrom(const std::string& fname); +}; + +#endif diff --git a/c_tools/zobov2/voz1b1/vozutil.c b/c_tools/zobov2/voz1b1/vozutil.c new file mode 100644 index 0000000..c56da6f --- /dev/null +++ b/c_tools/zobov2/voz1b1/vozutil.c @@ -0,0 +1,245 @@ +#include "libqhull/qhull_a.h" +#include "voz.h" + +#define FOREACHvertex2_(vertices) FOREACHsetelement_(vertexT, vertices2,vertex2) + +/*char qh_version[] = "user_eg 3.1 2001/10/04"; */ + +int compar(const void * n1, const void * n2) { + int i1,i2; + + i1 = *(int *)n1; + i2 = *(int *)n2; + return 2*(i1 > i2) - 1 + (i1 == i2); +} + +/* Finds Delaunay adjacencies of a set of points */ +int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs) { + + int dim= 3; /* dimension of points */ + boolT ismalloc= False; /* True if qhull should free points in qh_freeqhull() or reallocation */ + char flags[250]; /* option flags for qhull, see qh_opt.htm */ + FILE *outfile= stdout; /* output from qh_produce_output() + use NULL to skip qh_produce_output() */ + FILE *errfile= stderr; /* error messages from qhull code */ + int exitcode; /* 0 if no error from qhull */ + int curlong, totlong; /* memory remaining after qh_memfreeshort */ + int i, ver, count; + int numfacets, numsimplicial, numridges, totneighbors, numneighbors, + numcoplanars, numtricoplanars; + setT *vertices, *vertices2, *vertex_points, *coplanar_points; + vertexT *vertex, **vertexp; + vertexT *vertex2, **vertex2p; + int vertex_i, vertex_n; + facetT *facet, **facetp, *neighbor, **neighborp; + pointT *point, **pointp; + int numdiv; + + PARTADJ adjst; + + int errorReported = 0; + + adjst.adj = (int *)malloc(MAXVERVER*sizeof(int)); + if (adjst.adj == NULL) { + printf("Unable to allocate adjst.adj\n"); + exit(0); + } + + /* Delaunay triangulation*/ + sprintf (flags, "qhull s d Qt"); + exitcode= qh_new_qhull (dim, nvpall, points, ismalloc, + flags, outfile, errfile); + + + if (!exitcode) { /* if no error */ + /* 'qh facet_list' contains the convex hull */ + + /* From qh_printvneighbors */ + qh_countfacets(qh facet_list, NULL, 0, &numfacets, &numsimplicial, + &totneighbors, &numridges, &numcoplanars, &numtricoplanars); + qh_vertexneighbors(); + vertices= qh_facetvertices (qh facet_list, NULL, 0); + vertex_points= qh_settemp (nvpall); + coplanar_points= qh_settemp (nvpall); + qh_setzero (vertex_points, 0, nvpall); + qh_setzero (coplanar_points, 0, nvpall); + FOREACHvertex_(vertices) + qh_point_add (vertex_points, vertex->point, vertex); + FORALLfacet_(qh facet_list) { + FOREACHpoint_(facet->coplanarset) + qh_point_add (coplanar_points, point, facet); + } + ver = 0; + FOREACHvertex_i_(vertex_points) { + (*adjs)[ver].nadj = 0; + if (vertex) { + /* Count the neighboring vertices, check that all are real + neighbors */ + adjst.nadj = 0; + FOREACHneighbor_(vertex) { + if ((*adjs)[ver].nadj > -1) { + if (neighbor->visitid) { + vertices2 = neighbor->vertices; + FOREACHvertex2_(vertices2) { + if (ver != qh_pointid(vertex2->point)) { + adjst.adj[adjst.nadj] = (int)qh_pointid(vertex2->point); + adjst.nadj ++; + if (adjst.nadj > MAXVERVER) { + printf("Increase MAXVERVER to at least %d!\n",adjst.nadj); + exit(0); + } + } + } + } else { + printf(" %d",ver); + (*adjs)[ver].nadj = -1; /* There are unreal vertices here */ + } + } + } + } else (*adjs)[ver].nadj = -2; + + /* Enumerate the unique adjacencies*/ + if (adjst.nadj >= 4) { + qsort((void *)adjst.adj, adjst.nadj, sizeof(int), &compar); + count = 1; + + for (i=1; i= nvpbuf && !errorReported) { + errorReported = 1; + printf("Guard point encountered. Increase border and/or nguard.\n"); + printf("P:(%f,%f,%f), G: (%f,%f,%f)\n",points[3*ver],points[3*ver+1],points[3*ver+2], + points[3*adjst.adj[i]],points[3*adjst.adj[i]+1],points[3*adjst.adj[i]+2]); + } + count++; + } + (*adjs)[ver].adj = (int *)malloc(count*sizeof(int)); + if ((*adjs)[ver].adj == NULL) { + printf("Unable to allocate (*adjs)[ver].adj\n"); + exit(0); + } + (*adjs)[ver].adj[0] = adjst.adj[0]; + count = 1; + for (i=1; i %d\n",adjst.nadj,ver,ver); + exit(0); + } + ver++; + if (ver == nvp) break; + } + qh_settempfree (&coplanar_points); + qh_settempfree (&vertex_points); + qh_settempfree (&vertices); + } + qh_freeqhull(!qh_ALL); /* free long memory */ + qh_memfreeshort (&curlong, &totlong); /* free short memory and memory allocator */ + if (curlong || totlong) + fprintf (errfile, "qhull internal warning (delaunadj): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong); + free(adjst.adj); + return exitcode; +} + +/* Calculates the Voronoi volume from a set of Delaunay adjacencies */ +int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol) { + int dim= 3; /* dimension of points */ + boolT ismalloc= False; /* True if qhull should free points in qh_freeqhull() or reallocation */ + char flags[250]; /* option flags for qhull, see qh_opt.htm */ + FILE *outfile= NULL; /* output from qh_produce_output() + use NULL to skip qh_produce_output() */ + FILE *errfile= stderr; /* error messages from qhull code */ + int exitcode; /* 0 if no error from qhull */ + facetT *facet; /* set by FORALLfacets */ + int curlong, totlong; /* memory remaining after qh_memfreeshort */ + + coordT *point, *normp, *coordp, *feasiblep, *deladj; + int i, j, k; + boolT zerodiv; + float runsum; + char region; + /*coordT *points; + pointT *intpoints;*/ + + /* make point array from adjacency coordinates (add offset)*/ + /*points = (coordT *)malloc(4*numpoints*sizeof(coordT)); + if (points == NULL) { + printf("Unable to allocate points\n"); + exit(0); + }*/ + for (i=0; inormal; + feasiblep= qh feasible_point; + if (facet->offset < -qh MINdenom) { + for (k= qh hull_dim; k--; ) + *(coordp++)= (*(normp++) / - facet->offset) + *(feasiblep++); + }else { + for (k= qh hull_dim; k--; ) { + *(coordp++)= qh_divzero (*(normp++), facet->offset, qh MINdenom_1, + &zerodiv) + *(feasiblep++); + if (zerodiv) { + qh_memfree (point, qh normal_size); + printf("LABELprintinfinite\n"); + exit(0); + } + } + } + } + } + qh_freeqhull (!qh_ALL); + qh_memfreeshort (&curlong, &totlong); + + /* Now we calculate the volume */ + sprintf (flags, "qhull FA"); + exitcode= qh_new_qhull (dim, numpoints, intpoints, ismalloc, + flags, outfile, errfile); + + qh_getarea(qh facet_list); + *vol = qh totvol; + + qh_freeqhull (!qh_ALL); + qh_memfreeshort (&curlong, &totlong); + if (curlong || totlong) + fprintf (errfile, "qhull internal warning (vorvol): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong); + /*free(points); free(intpoints);*/ + + return exitcode; +} From 2ad09d3a269f429c27c4af3be8a09b1d03734a28 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 10 Jun 2013 16:36:16 +0200 Subject: [PATCH 4/9] Moved jozov2 into its own directory. More code reorganization in voz1b1. --- c_tools/zobov2/{ => jozov2}/findrtop.c | 0 c_tools/zobov2/{ => jozov2}/jozov2.cpp | 0 c_tools/zobov2/{ => jozov2}/jozov2.hpp | 0 c_tools/zobov2/{ => jozov2}/jozov2_io.cpp | 0 .../zobov2/{ => jozov2}/jozov2_watershed.cpp | 0 c_tools/zobov2/{ => jozov2}/jozov2_zones.cpp | 0 c_tools/zobov2/{ => jozov2}/zobov.hpp | 0 c_tools/zobov2/voz1b1/voz1b1.cpp | 511 ++++++++++-------- c_tools/zobov2/voz1b1/voz_io.hpp | 1 + 9 files changed, 275 insertions(+), 237 deletions(-) rename c_tools/zobov2/{ => jozov2}/findrtop.c (100%) rename c_tools/zobov2/{ => jozov2}/jozov2.cpp (100%) rename c_tools/zobov2/{ => jozov2}/jozov2.hpp (100%) rename c_tools/zobov2/{ => jozov2}/jozov2_io.cpp (100%) rename c_tools/zobov2/{ => jozov2}/jozov2_watershed.cpp (100%) rename c_tools/zobov2/{ => jozov2}/jozov2_zones.cpp (100%) rename c_tools/zobov2/{ => jozov2}/zobov.hpp (100%) diff --git a/c_tools/zobov2/findrtop.c b/c_tools/zobov2/jozov2/findrtop.c similarity index 100% rename from c_tools/zobov2/findrtop.c rename to c_tools/zobov2/jozov2/findrtop.c diff --git a/c_tools/zobov2/jozov2.cpp b/c_tools/zobov2/jozov2/jozov2.cpp similarity index 100% rename from c_tools/zobov2/jozov2.cpp rename to c_tools/zobov2/jozov2/jozov2.cpp diff --git a/c_tools/zobov2/jozov2.hpp b/c_tools/zobov2/jozov2/jozov2.hpp similarity index 100% rename from c_tools/zobov2/jozov2.hpp rename to c_tools/zobov2/jozov2/jozov2.hpp diff --git a/c_tools/zobov2/jozov2_io.cpp b/c_tools/zobov2/jozov2/jozov2_io.cpp similarity index 100% rename from c_tools/zobov2/jozov2_io.cpp rename to c_tools/zobov2/jozov2/jozov2_io.cpp diff --git a/c_tools/zobov2/jozov2_watershed.cpp b/c_tools/zobov2/jozov2/jozov2_watershed.cpp similarity index 100% rename from c_tools/zobov2/jozov2_watershed.cpp rename to c_tools/zobov2/jozov2/jozov2_watershed.cpp diff --git a/c_tools/zobov2/jozov2_zones.cpp b/c_tools/zobov2/jozov2/jozov2_zones.cpp similarity index 100% rename from c_tools/zobov2/jozov2_zones.cpp rename to c_tools/zobov2/jozov2/jozov2_zones.cpp diff --git a/c_tools/zobov2/zobov.hpp b/c_tools/zobov2/jozov2/zobov.hpp similarity index 100% rename from c_tools/zobov2/zobov.hpp rename to c_tools/zobov2/jozov2/zobov.hpp diff --git a/c_tools/zobov2/voz1b1/voz1b1.cpp b/c_tools/zobov2/voz1b1/voz1b1.cpp index 7912082..2b996f9 100644 --- a/c_tools/zobov2/voz1b1/voz1b1.cpp +++ b/c_tools/zobov2/voz1b1/voz1b1.cpp @@ -1,10 +1,13 @@ +#include #include +#include #include #include #include "libqhull/qhull_a.h" #include "voz.h" #include "voz_io.hpp" +using boost::format; using namespace std; #define DL for (d=0;d<3;d++) @@ -19,33 +22,208 @@ bool checkParameters(int *numdiv, int *b) return true; } +struct BoxData +{ + float width[3], totwidth[3]; + float width2[3], totwidth2[3]; + float bf, s[3], g[3], c[3]; + coordT *parts; + pid_t *orig; + pid_t nvpall, nvp, nvpbuf; + bool guard_added; + double xyz_min[3], xyz_max[3]; + + void prepareBox(const PositionData& pdata); + + void checkParticle(float *xyz, bool& in_main, bool& in_buf); +}; + +void BoxData::checkParticle(double *xyz, bool& in_main, bool& in_buf) +{ + in_main = in_buf = true; + for (int d = 0; d < 3; d++) + { + xyz[d] -= (double)c[d]; + if (xyz[d] > width2[d]) + xyz[d] -= width[d]; + if (xyz[d] < -width2[d]) + xyz[d] += width[d]; + + in_buf = in_buf && (fabs(xyz[d]) < totwidth2[d]); + in_main = in_main && (fabs(xyz[d]) <= width2[d]); + } +} + +void BoxData::prepareBox(const PositionData& pdata, int *b_id) +{ + guard_added = false; + + for (int i = 0; i < 3; i++) + { + float s; + width[i] = (pdata.xyz_max[i] - pdata.xyz_min[i])/(float)numdiv; + width2[i] = 0.5*width; + + totwidth[i] = width+2.*bf; + totwidth2[i] = width2 + bf; + + s[i] = width[i]/(float)NGUARD; + if ((bf*bf - 2.*s[i]*s[i]) < 0.) + { + printf("bf = %f, s = %f.\n",bf,s[i]); + printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.", + sqrt(2.)*width/bf); + exit(0); + } + g[i] = (bf / 2.)*(1. + sqrt(1 - 2.*s[i]*s[i]/(bf*bf))); + cout << format("s[%d] = %f, bf = %f, g[%d] = %f.") % s[i] % i % bf % g[i] % i << endl; + + c[i] = b_id[i] * width[i]; + } + cout.flush(); + + cout << format("c: %f,%f,%f") % c[0] % c[1] % c[2] << endl; + + /* Assign temporary array*/ + nvpbuf = 0; /* Number of particles to tesselate, including buffer */ + nvp = 0; /* Without the buffer */ + + for (pid_t i=0; i < np; i++) + { + double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] }; + bool is_it_in_buf, is_it_in_main; + + checkParticle(xyz, is_it_in_buf, is_it_in_main); + + if (isitinbuf) + nvpbuf++; + if (isitinmain) + nvp++; + } + + nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard points */ + + parts = new coordT[3*nvpbuf]; + orig = new pid_t[nvpuf]; + + if (parts == 0) + { + cout << "Unable to allocate parts" << endl; + exit(0); + } + if (orig == 0) + { + cout << "Unable to allocate orig" << endl; + exit(0); + } + + nvp = 0; nvpall = 0; /* nvp = number of particles without buffer */ + + for (int j = 0; j < 3; j++) + { + xyz_min[j] = -std::numeric_limits::max(); + xyz_max[j] = std::numeric_limits::max(); + } + for (pid_t i = 0; i < np; i++) + { + bool is_it_in_main, is_it_in_buf; + double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] }; + + checkParticle(xyz, is_it_in_main, is_it_in_buf); + + if (is_it_in_main) { + for (int j = 0; j < 3; j++) + { + parts[3*nvp+j] = xyz[j]; + xyz_min[j] = min(xyz_min[j], xyz[j]); + xyz_max[j] = max(xyz_max[j], xyz[j]); + } + orig[nvp] = i; + nvp++; + } + } + printf("nvp = %d\n",nvp); + printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax); + nvpbuf = nvp; + for (pid_t i = 0; i < np; i++) + { + bool is_it_in_main, is_it_in_buf; + double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] }; + + checkParticle(xyz, is_it_in_main, is_it_in_buf); + + if (isitinbuf && !isitinmain) + { + for (int j = 0; j < 3; j++) + { + parts[3*nvpbuf+j] = xyz[j]; + xyz_min[j] = min(xyz_min[j], xyz[j]); + xyz_max[j] = max(xyz_max[j], xyz[j]); + } + orig[nvpbuf] = i; + + nvpbuf++; + } + } + cout << format("nvpbuf = %d") % nvpbuf << endl; + cout << format("x: %f,%f; y: %f,%f; z:%f,%f\n") % xyz_min[0] % xyz_max[0] % xyz_min[1] % xyz_max[1] % xyz_min[2] % xyz_max[2] << endl; + nvpall = nvpbuf; + + double predict = 1; + for (int j = 0; j < 3; j++) + predict *= (width[j]+2*bf); + predict *= np; + cout << format("There should be ~ %g points; there are %d\n") % predict % nvpbuf << endl; +} + +void BoxData::addGuardPoints() +{ + int rot_g[3][3] = { {0, 0, 1}, {0,1,0}, {1,0,0} }; + int rot_si[3][3] = { {1, 0, 0}, {1,0,0}, {0,1,0} }; + int rot_sj[3][3] = { {0, 1, 0}, {0,0,1}, {0,0,1} }; + + for (int k = 0; k < 3; k++) + { + + /* Add guard points */ + for (i=0; i<=NGUARD; i++) + { + for (j=0; j<=NGUARD; j++) + { + /* Bottom */ + for (int a = 0; a < 3; a++) + parts[3*nvpall+a] = -width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] - rot_g[k][a] * g[a]; + + nvpall++; + + /* Top */ + for (int a = 0; a < 3; a++) + parts[3*nvpall+a] = (2*rot_g[k][a]-1)*width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] + rot_g[k][a] * g[a]; + + nvpall++; + } + } + } +} + int main(int argc, char *argv[]) { - const float BF = std::numeric_limits::max(); - - int exitcode; - pid_t i, j, np; PositionData pdata; - coordT rtemp[3], *parts; + BoxData boxdata; + coordT deladjs[3*MAXVERVER], points[3*MAXVERVER]; pointT intpoints[3*MAXVERVER]; - FILE *pos, *out; - char *posfile, outfile[200], *suffix, *outDir; + string outfile; + ofstream out; + + char *suffix, *outDir; PARTADJ *adjs; float *vols; - realT predict, xmin,xmax,ymin,ymax,zmin,zmax; pid_t *orig; - - int isitinbuf; - char isitinmain, d; - pid_t nvp, nvpall, nvpbuf; - realT width, width2, totwidth, totwidth2, bf, s, g; double border, boxsize[3]; - realT c[3]; - int b[3]; - int numdiv[3]; + int b[3], numdiv[3]; double totalvol; - + CosmoTool::MiniArgDesc args[] = { { "POSITION FILE", MINIARG_STRING, &posfile }, { "BORDER SIZE", MINIARG_DOUBLE, &border }, @@ -68,7 +246,7 @@ int main(int argc, char *argv[]) { if (!checkParameters(numdiv, b)) return 2; - + /* Boxsize should be the range in r, yielding a range 0-1 */ if (!pdata.readFrom(posfile)) return 3; @@ -77,180 +255,29 @@ int main(int argc, char *argv[]) { pdata.findExtrema(); - (cout << boost::format("np: %d, x: %f,%f; y: %f,%f; z: %f,%f") + (cout << boost::format("np: %d, x: %f,%f; y: %f,%f; z: %f,%f") % pdata.np % xyz_min[0] % xyz_max[0] % xyz_min[1] % xyz_max[1] % xyz_min[2] % xyz_max[2]).flush(); - width = 1./(float)numdiv; - width2 = 0.5*width; - if (border > 0.) bf = border; - else bf = 0.1; - /* In units of 0-1, the thickness of each subregion's buffer*/ - totwidth = width+2.*bf; - totwidth2 = width2 + bf; - - s = width/(float)NGUARD; - if ((bf*bf - 2.*s*s) < 0.) { - printf("bf = %f, s = %f.\n",bf,s); - printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.", - sqrt(2.)*width/bf); - exit(0); - } - g = (bf / 2.)*(1. + sqrt(1 - 2.*s*s/(bf*bf))); - printf("s = %f, bf = %f, g = %f.\n",s,bf,g); - - fflush(stdout); - - adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ)); - if (adjs == NULL) { - printf("Unable to allocate adjs\n"); - exit(0); - } - - DL c[d] = ((float)b[d])*width; - printf("c: %f,%f,%f\n",c[0],c[1],c[2]); - /* Assign temporary array*/ - nvpbuf = 0; /* Number of particles to tesselate, including - buffer */ - nvp = 0; /* Without the buffer */ - for (i=0; i 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 (isitinmain) nvp++; - } - - 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)); - - if (parts == NULL) { - printf("Unable to allocate parts\n"); - fflush(stdout); - } - if (orig == NULL) { - printf("Unable to allocate orig\n"); - fflush(stdout); - } - - 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 0.5) rtemp[d] --; - if (rtemp[d] < -0.5) rtemp[d] ++; - isitinmain = isitinmain && (fabs(rtemp[d]) <= width2); - } - if (isitinmain) { - parts[3*nvp] = rtemp[0]; - parts[3*nvp+1] = rtemp[1]; - parts[3*nvp+2] = rtemp[2]; - orig[nvp] = i; - nvp++; - if (rtemp[0] < xmin) xmin = rtemp[0]; - if (rtemp[0] > xmax) xmax = rtemp[0]; - if (rtemp[1] < ymin) ymin = rtemp[1]; - if (rtemp[1] > ymax) ymax = rtemp[1]; - if (rtemp[2] < zmin) zmin = rtemp[2]; - 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); - nvpbuf = nvp; - for (i=0; i 0.5) rtemp[d] --; - if (rtemp[d] < -0.5) rtemp[d] ++; - isitinbuf = isitinbuf && (fabs(rtemp[d]) xmax) xmax = rtemp[0]; - if (rtemp[1] < ymin) ymin = rtemp[1]; - if (rtemp[1] > ymax) ymax = rtemp[1]; - if (rtemp[2] < zmin) zmin = rtemp[2]; - 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); + if (border > 0.) + boxdata.bf = border; + else + boxdata.bf = 0.1; + boxdata.prepareBox(pdata); pdata.destroy(); - - /* Add guard points */ - for (i=0; i zmax) zmax = parts[3*i+2]; } - - 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); - + + cout << format("Added guard points to total %d points (should be %d)") + % boxdata.nvpall % (boxdata.nvpbuf + 6*(NGUARD+1)*(NGUARD+1)) << endl; + cout << format("x: %f,%f; y: %f,%f; z:%f,%f") % boxdata.xyz_min[0] % boxdata.xyz_max[0] % boxdata.xyz_min[1] % boxdata.xyz_max[1] % boxdata.xyz_min[2] % boxdata.xyz_max[2] << endl; + /* Do tesselation*/ printf("File read. Tessellating ...\n"); fflush(stdout); - exitcode = delaunadj(parts, nvp, nvpbuf, nvpall, &adjs); + exitcode = delaunadj(boxdata.parts, boxdata.nvp, boxdata.nvpbuf, boxdata.nvpall, &adjs); if (exitcode != 0) { printf("Error while tesselating. Stopping here."); fflush(stdout); - exit(1); + return 4; } - + /* Now calculate volumes*/ printf("Now finding volumes ...\n"); fflush(stdout); - vols = (float *)malloc(nvp*sizeof(float)); - - for (i=0; i 0.5) deladjs[3*j+d]--; - } - - 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]); - } + vols = new float[nvp]; + + for (pid_t i = 0; i < nvp; i++) + { /* Just the original particles + * Assign adjacency coordinate array*/ + /* Volumes */ + for (int j = 0; j < adjs[i].nadj; j++) + { + for (int d = 0; d < 3; d++) + { + deladjs[3*j + d] = parts[3*adjs[i].adj[j]+d] - parts[3*i+d]; + + if (deladjs[3*j+d] < -boxdata.width2[d]) + deladjs[3*j+d] += boxdata.width[d]; + if (deladjs[3*j+d] > boxdata.width2[d]) + deladjs[3*j+d] -= boxdata.width[d]; + } + } + + exitcode = vorvol(deladjs, points, intpoints, adjs[i].nadj, &(vols[i])); + vols[i] *= np/V0; + if ((i % 1000) == 0) + cout << format("%d: %d, %f") % i % adjs[i].nadj % vols[i] << endl; + } /* Get the adjacencies back to their original values */ - - for (i=0; i 0) - fwrite(adjs[i].adj,adjs[i].nadj,sizeof(pid_t),out); - else printf("0"); - } - fclose(out); - + for (pid_t i = 0; i < nvp; i++) + { + out.write((char*)&adjs[i].nadj, sizeof(pid_t)); + if (adjs[i].nadj > 0) + out.write((char *)adjs[i].adj, adjs[i].nadj*sizeof(pid_t)); + else + (cout << "0").flush(); + } + out.close(); + delete[] adjs; + return(0); } diff --git a/c_tools/zobov2/voz1b1/voz_io.hpp b/c_tools/zobov2/voz1b1/voz_io.hpp index dbeb632..01205bb 100644 --- a/c_tools/zobov2/voz1b1/voz_io.hpp +++ b/c_tools/zobov2/voz1b1/voz_io.hpp @@ -7,6 +7,7 @@ struct PositionData { float *xyz[3]; pid_t np; + float xyz_min[3], xyz_max[3]; void destroy() { From 51efae2e0912736f7f22d725ef415987d29d4007 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 10 Jun 2013 16:48:29 +0200 Subject: [PATCH 5/9] Tons of fixes --- c_tools/zobov2/voz1b1/voz1b1.cpp | 77 +++++++++++++++++--------------- 1 file changed, 40 insertions(+), 37 deletions(-) diff --git a/c_tools/zobov2/voz1b1/voz1b1.cpp b/c_tools/zobov2/voz1b1/voz1b1.cpp index 2b996f9..63348dd 100644 --- a/c_tools/zobov2/voz1b1/voz1b1.cpp +++ b/c_tools/zobov2/voz1b1/voz1b1.cpp @@ -7,6 +7,10 @@ #include "voz.h" #include "voz_io.hpp" +using CosmoTool::MINIARG_STRING; +using CosmoTool::MINIARG_DOUBLE; +using CosmoTool::MINIARG_INT; +using CosmoTool::MINIARG_NULL; using boost::format; using namespace std; @@ -33,9 +37,9 @@ struct BoxData bool guard_added; double xyz_min[3], xyz_max[3]; - void prepareBox(const PositionData& pdata); + void prepareBox(const PositionData& pdata, int *b_id, int *numdivs); - void checkParticle(float *xyz, bool& in_main, bool& in_buf); + void checkParticle(double *xyz, bool& in_main, bool& in_buf); }; void BoxData::checkParticle(double *xyz, bool& in_main, bool& in_buf) @@ -54,25 +58,24 @@ void BoxData::checkParticle(double *xyz, bool& in_main, bool& in_buf) } } -void BoxData::prepareBox(const PositionData& pdata, int *b_id) +void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) { guard_added = false; for (int i = 0; i < 3; i++) { - float s; - width[i] = (pdata.xyz_max[i] - pdata.xyz_min[i])/(float)numdiv; - width2[i] = 0.5*width; + width[i] = (pdata.xyz_max[i] - pdata.xyz_min[i])/numdivs[i]; + width2[i] = 0.5*width[i]; - totwidth[i] = width+2.*bf; - totwidth2[i] = width2 + bf; + totwidth[i] = width[i]+2.*bf; + totwidth2[i] = width2[i] + bf; s[i] = width[i]/(float)NGUARD; if ((bf*bf - 2.*s[i]*s[i]) < 0.) { printf("bf = %f, s = %f.\n",bf,s[i]); printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.", - sqrt(2.)*width/bf); + sqrt(2.)*width[i]/bf); exit(0); } g[i] = (bf / 2.)*(1. + sqrt(1 - 2.*s[i]*s[i]/(bf*bf))); @@ -88,23 +91,23 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id) nvpbuf = 0; /* Number of particles to tesselate, including buffer */ nvp = 0; /* Without the buffer */ - for (pid_t i=0; i < np; i++) + for (pid_t i = 0; i < pdata.np; i++) { double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] }; bool is_it_in_buf, is_it_in_main; checkParticle(xyz, is_it_in_buf, is_it_in_main); - if (isitinbuf) + if (is_it_in_buf) nvpbuf++; - if (isitinmain) + if (is_it_in_main) nvp++; } nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard points */ parts = new coordT[3*nvpbuf]; - orig = new pid_t[nvpuf]; + orig = new pid_t[nvpbuf]; if (parts == 0) { @@ -124,7 +127,7 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id) xyz_min[j] = -std::numeric_limits::max(); xyz_max[j] = std::numeric_limits::max(); } - for (pid_t i = 0; i < np; i++) + for (pid_t i = 0; i < pdata.np; i++) { bool is_it_in_main, is_it_in_buf; double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] }; @@ -142,17 +145,17 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id) nvp++; } } - printf("nvp = %d\n",nvp); - printf("x: %f,%f; y: %f,%f; z:%f,%f\n",xmin,xmax,ymin,ymax,zmin,zmax); + cout << format("nvp = %d") %nvp << endl; + cout << format("x: %f,%f; y: %f,%f; z:%f,%f") % xyz_min[0] % xyz_max[0] % xyz_min[1] % xyz_max[1] % xyz_min[2] % xyz_max[2] << endl; nvpbuf = nvp; - for (pid_t i = 0; i < np; i++) + for (pid_t i = 0; i < pdata.np; i++) { bool is_it_in_main, is_it_in_buf; double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] }; checkParticle(xyz, is_it_in_main, is_it_in_buf); - if (isitinbuf && !isitinmain) + if (is_it_in_buf && !is_it_in_main) { for (int j = 0; j < 3; j++) { @@ -172,7 +175,7 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id) double predict = 1; for (int j = 0; j < 3; j++) predict *= (width[j]+2*bf); - predict *= np; + predict *= pdata.np; cout << format("There should be ~ %g points; there are %d\n") % predict % nvpbuf << endl; } @@ -216,7 +219,7 @@ int main(int argc, char *argv[]) { string outfile; ofstream out; - char *suffix, *outDir; + char *suffix, *outDir, *posfile; PARTADJ *adjs; float *vols; pid_t *orig; @@ -225,20 +228,20 @@ int main(int argc, char *argv[]) { double totalvol; CosmoTool::MiniArgDesc args[] = { - { "POSITION FILE", MINIARG_STRING, &posfile }, - { "BORDER SIZE", MINIARG_DOUBLE, &border }, - { "BOX_X", MINIARG_DOUBLE, &boxsize[0] }, - { "BOX_Y", MINIARG_DOUBLE, &boxsize[1] }, - { "BOX_Z", MINIARG_DOUBLE, &boxsize[2] }, - { "SUFFIX", MINIARG_STRING, &suffix }, - { "NUM_DIVISION_X", MINIARG_INT, &numdiv[0] }, - { "NUM_DIVISION_Y", MINIARG_INT, &numdiv[1] }, - { "NUM_DIVISION_Z", MINIARG_INT, &numdiv[2] }, - { "B0", MINIARG_INT, &b[0] }, - { "B1", MINIARG_INT, &b[1] }, - { "B2", MINIARG_INT, &b[2] }, - { "OUTPUT DIRECTORY", MINIARG_STRING, &outDir }, - { 0, MINIARG_NULL, 0 } + { "POSITION FILE", &posfile, MINIARG_STRING }, + { "BORDER SIZE", &border, MINIARG_DOUBLE }, + { "BOX_X", &boxsize[0], MINIARG_DOUBLE }, + { "BOX_Y", &boxsize[1], MINIARG_DOUBLE }, + { "BOX_Z", &boxsize[2], MINIARG_DOUBLE }, + { "SUFFIX", &suffix, MINIARG_STRING }, + { "NUM_DIVISION_X", &numdiv[0], MINIARG_INT }, + { "NUM_DIVISION_Y", &numdiv[1], MINIARG_INT }, + { "NUM_DIVISION_Z", &numdiv[2], MINIARG_INT }, + { "B0", &b[0], MINIARG_INT }, + { "B1", &b[1], MINIARG_INT }, + { "B2", &b[2], MINIARG_INT }, + { "OUTPUT DIRECTORY", &outDir, MINIARG_STRING }, + { 0, 0, MINIARG_NULL } }; if (!CosmoTool::parseMiniArgs(argc, argv, args)) @@ -257,9 +260,9 @@ int main(int argc, char *argv[]) { (cout << boost::format("np: %d, x: %f,%f; y: %f,%f; z: %f,%f") % pdata.np - % xyz_min[0] % xyz_max[0] - % xyz_min[1] % xyz_max[1] - % xyz_min[2] % xyz_max[2]).flush(); + % pdata.xyz_min[0] % pdata.xyz_max[0] + % pdata.xyz_min[1] % pdata.xyz_max[1] + % pdata.xyz_min[2] % pdata.xyz_max[2]).flush(); if (border > 0.) From 26132586c056367c8289c2e2e398864af1e1199e Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 11 Jun 2013 10:26:50 +0200 Subject: [PATCH 6/9] voz1b1_2 compiles and links now. --- c_tools/zobov2/voz1b1/voz1b1.cpp | 108 ++++++++++++++++++------------- c_tools/zobov2/voz1b1/voz_io.cpp | 13 ++-- c_tools/zobov2/voz1b1/voz_io.hpp | 3 +- 3 files changed, 72 insertions(+), 52 deletions(-) diff --git a/c_tools/zobov2/voz1b1/voz1b1.cpp b/c_tools/zobov2/voz1b1/voz1b1.cpp index 63348dd..0bad2ce 100644 --- a/c_tools/zobov2/voz1b1/voz1b1.cpp +++ b/c_tools/zobov2/voz1b1/voz1b1.cpp @@ -26,6 +26,7 @@ bool checkParameters(int *numdiv, int *b) return true; } + struct BoxData { float width[3], totwidth[3]; @@ -40,6 +41,8 @@ struct BoxData void prepareBox(const PositionData& pdata, int *b_id, int *numdivs); void checkParticle(double *xyz, bool& in_main, bool& in_buf); + void addGuardPoints(); + void findBoundary(); }; void BoxData::checkParticle(double *xyz, bool& in_main, bool& in_buf) @@ -189,9 +192,9 @@ void BoxData::addGuardPoints() { /* Add guard points */ - for (i=0; i<=NGUARD; i++) + for (int i=0; i <= NGUARD; i++) { - for (j=0; j<=NGUARD; j++) + for (int j=0; j <= NGUARD; j++) { /* Bottom */ for (int a = 0; a < 3; a++) @@ -209,6 +212,52 @@ void BoxData::addGuardPoints() } } +void BoxData::findBoundary() +{ + double BF = std::numeric_limits::max(); + + for (int j = 0; j < 3; j++) + { + xyz_min[j] = -BF; + xyz_max[j] = BF; + } + + for (pid_t i = nvpbuf; i < nvpall; i++) { + for (int j = 0; j < 3; j++) + { + xyz_min[j] = std::min(xyz_min[j], parts[3*i+j]); + xyz_max[j] = std::min(xyz_max[j], parts[3*i+j]); + } + } +} + +void saveTesselation(const string& outfile, PositionData& pdata, BoxData& boxdata, PARTADJ *adjs, float *vols) +{ + ofstream out(outfile.c_str()); + if (!out) + { + cout << format("Unable to open %s") % outfile << endl; + exit(0); + } + out.write((char *)&pdata.np, sizeof(pid_t)); + out.write((char *)&boxdata.nvp, sizeof(pid_t)); + cout << format("nvp = %d") % boxdata.nvp << endl; + + /* Tell us where the original particles were */ + out.write((char *)boxdata.orig, sizeof(pid_t)*boxdata.nvp); + /* Volumes*/ + out.write((char *)vols,sizeof(float)*boxdata.nvp); + /* Adjacencies */ + for (pid_t i = 0; i < boxdata.nvp; i++) + { + out.write((char*)&adjs[i].nadj, sizeof(pid_t)); + if (adjs[i].nadj > 0) + out.write((char *)adjs[i].adj, adjs[i].nadj*sizeof(pid_t)); + else + (cout << "0").flush(); + } + out.close(); +} int main(int argc, char *argv[]) { PositionData pdata; @@ -270,26 +319,18 @@ int main(int argc, char *argv[]) { else boxdata.bf = 0.1; - boxdata.prepareBox(pdata); + boxdata.prepareBox(pdata, b, numdiv); pdata.destroy(); boxdata.addGuardPoints(); - adjs = new PARTADJ[np]; + adjs = new PARTADJ[boxdata.nvpall]; if (adjs == 0) { cout << "Unable to allocate adjs" << endl; return 0; } - xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF; - for (i=nvpbuf;i xmax) xmax = parts[3*i]; - if (parts[3*i+1] < ymin) ymin = parts[3*i+1]; - if (parts[3*i+1] > ymax) ymax = parts[3*i+1]; - if (parts[3*i+2] < zmin) zmin = parts[3*i+2]; - if (parts[3*i+2] > zmax) zmax = parts[3*i+2]; - } + boxdata.findBoundary(); cout << format("Added guard points to total %d points (should be %d)") % boxdata.nvpall % (boxdata.nvpbuf + 6*(NGUARD+1)*(NGUARD+1)) << endl; @@ -297,7 +338,7 @@ int main(int argc, char *argv[]) { /* Do tesselation*/ printf("File read. Tessellating ...\n"); fflush(stdout); - exitcode = delaunadj(boxdata.parts, boxdata.nvp, boxdata.nvpbuf, boxdata.nvpall, &adjs); + int exitcode = delaunadj(boxdata.parts, boxdata.nvp, boxdata.nvpbuf, boxdata.nvpall, &adjs); if (exitcode != 0) { printf("Error while tesselating. Stopping here."); fflush(stdout); @@ -306,9 +347,9 @@ int main(int argc, char *argv[]) { /* Now calculate volumes*/ printf("Now finding volumes ...\n"); fflush(stdout); - vols = new float[nvp]; + vols = new float[boxdata.nvp]; - for (pid_t i = 0; i < nvp; i++) + for (pid_t i = 0; i < boxdata.nvp; i++) { /* Just the original particles * Assign adjacency coordinate array*/ /* Volumes */ @@ -316,7 +357,7 @@ int main(int argc, char *argv[]) { { for (int d = 0; d < 3; d++) { - deladjs[3*j + d] = parts[3*adjs[i].adj[j]+d] - parts[3*i+d]; + deladjs[3*j + d] = boxdata.parts[3*adjs[i].adj[j]+d] - boxdata.parts[3*i+d]; if (deladjs[3*j+d] < -boxdata.width2[d]) deladjs[3*j+d] += boxdata.width[d]; @@ -326,51 +367,28 @@ int main(int argc, char *argv[]) { } exitcode = vorvol(deladjs, points, intpoints, adjs[i].nadj, &(vols[i])); - vols[i] *= np/V0; + vols[i] *= pdata.np/pdata.V0; if ((i % 1000) == 0) cout << format("%d: %d, %f") % i % adjs[i].nadj % vols[i] << endl; } /* Get the adjacencies back to their original values */ - for (pid_t i=0; i 0) - out.write((char *)adjs[i].adj, adjs[i].nadj*sizeof(pid_t)); - else - (cout << "0").flush(); - } - out.close(); + saveTesselation(outfile, pdata, boxdata, adjs, vols); delete[] adjs; return(0); diff --git a/c_tools/zobov2/voz1b1/voz_io.cpp b/c_tools/zobov2/voz1b1/voz_io.cpp index 6bf033f..712bd5b 100644 --- a/c_tools/zobov2/voz1b1/voz_io.cpp +++ b/c_tools/zobov2/voz1b1/voz_io.cpp @@ -1,3 +1,4 @@ +#include #include #include #include "voz_io.hpp" @@ -12,15 +13,15 @@ bool PositionData::readFrom(const string& fname) try { UnformattedRead f(fname.c_str()); - + f.beginCheckpoint(); np = f.readInt32(); - f.endCheckPoint(); + f.endCheckpoint(); for (int j = 0; j < 3; j++) { xyz[j] = new float[np]; - + f.beginCheckpoint(); for (int p = 0; p < np; p++) xyz[j][p] = f.readReal32(); @@ -39,7 +40,7 @@ bool PositionData::readFrom(const string& fname) { return false; } - + return true; } @@ -49,13 +50,13 @@ void PositionData::findExtrema() for (int i = 0; i < 3; i++) { - xyz_min[i] = BF; + xyz_min[i] = BF; xyz_max[i] = -BF; } for (int i = 0; i < 3; i++) { - for (pid_t p = 0; p < np; p++) + for (pid_t p = 0; p < np; p++) { xyz_min[p] = min(xyz_min[p], xyz[p][i]); xyz_max[p] = max(xyz_max[p], xyz[p][i]); diff --git a/c_tools/zobov2/voz1b1/voz_io.hpp b/c_tools/zobov2/voz1b1/voz_io.hpp index 01205bb..77a1e5c 100644 --- a/c_tools/zobov2/voz1b1/voz_io.hpp +++ b/c_tools/zobov2/voz1b1/voz_io.hpp @@ -8,7 +8,8 @@ struct PositionData float *xyz[3]; pid_t np; float xyz_min[3], xyz_max[3]; - + float V0; + void destroy() { for (int j = 0; j < 3; j++) From 1553601028c119afe178521c893d4436b63841be Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 11 Jun 2013 10:48:04 +0200 Subject: [PATCH 7/9] Initialize volume of the box. --- c_tools/zobov2/voz1b1/voz_io.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/c_tools/zobov2/voz1b1/voz_io.cpp b/c_tools/zobov2/voz1b1/voz_io.cpp index 712bd5b..d8dea31 100644 --- a/c_tools/zobov2/voz1b1/voz_io.cpp +++ b/c_tools/zobov2/voz1b1/voz_io.cpp @@ -54,12 +54,14 @@ void PositionData::findExtrema() xyz_max[i] = -BF; } + V0=1; for (int i = 0; i < 3; i++) { for (pid_t p = 0; p < np; p++) { - xyz_min[p] = min(xyz_min[p], xyz[p][i]); - xyz_max[p] = max(xyz_max[p], xyz[p][i]); + xyz_min[i] = min(xyz_min[i], xyz[i][p]); + xyz_max[i] = max(xyz_max[i], xyz[i][p]); } + V0 *= (xyz_max[i]-xyz_min[i]); } } From 03d7fda57fe40108d142a1bad3bb85c3b6986768 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Tue, 11 Jun 2013 11:51:07 +0200 Subject: [PATCH 8/9] Fixes. Now run through classical testcase. --- c_tools/mock/loaders/sdf_loader.cpp | 2 +- c_tools/zobov2/voz1b1/voz1b1.cpp | 43 +++++++++++++++++------------ 2 files changed, 27 insertions(+), 18 deletions(-) diff --git a/c_tools/mock/loaders/sdf_loader.cpp b/c_tools/mock/loaders/sdf_loader.cpp index f2b0c32..55cb279 100644 --- a/c_tools/mock/loaders/sdf_loader.cpp +++ b/c_tools/mock/loaders/sdf_loader.cpp @@ -166,7 +166,7 @@ public: } if (load_flags & NEED_VELOCITY) { - const char *name_template[3] = { "x", "y", "z" }; + const char *name_template[3] = { "vx", "vy", "vz" }; for (int q = 0; q < 3; q++) { d->Vel[q] = new float[numPartToLoad]; diff --git a/c_tools/zobov2/voz1b1/voz1b1.cpp b/c_tools/zobov2/voz1b1/voz1b1.cpp index 0bad2ce..d9ff28a 100644 --- a/c_tools/zobov2/voz1b1/voz1b1.cpp +++ b/c_tools/zobov2/voz1b1/voz1b1.cpp @@ -1,3 +1,4 @@ +#include #include #include #include @@ -59,10 +60,12 @@ void BoxData::checkParticle(double *xyz, bool& in_main, bool& in_buf) in_buf = in_buf && (fabs(xyz[d]) < totwidth2[d]); in_main = in_main && (fabs(xyz[d]) <= width2[d]); } + assert(!in_main || in_buf); } void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) { + double BF = std::numeric_limits::max(); guard_added = false; for (int i = 0; i < 3; i++) @@ -82,7 +85,7 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) exit(0); } g[i] = (bf / 2.)*(1. + sqrt(1 - 2.*s[i]*s[i]/(bf*bf))); - cout << format("s[%d] = %f, bf = %f, g[%d] = %f.") % s[i] % i % bf % g[i] % i << endl; + cout << format("s[%d] = %f, bf = %f, g[%d] = %f.") % i% s[i] % bf % i % g[i] << endl; c[i] = b_id[i] * width[i]; } @@ -99,18 +102,20 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] }; bool is_it_in_buf, is_it_in_main; - checkParticle(xyz, is_it_in_buf, is_it_in_main); + checkParticle(xyz, is_it_in_main, is_it_in_buf); if (is_it_in_buf) nvpbuf++; if (is_it_in_main) nvp++; + assert(nvp <= nvpbuf); } nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard points */ + nvpall = nvpbuf; - parts = new coordT[3*nvpbuf]; - orig = new pid_t[nvpbuf]; + parts = new coordT[3*nvpall]; + orig = new pid_t[nvpall]; if (parts == 0) { @@ -123,12 +128,12 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) exit(0); } - nvp = 0; nvpall = 0; /* nvp = number of particles without buffer */ + nvp = 0; /* nvp = number of particles without buffer */ for (int j = 0; j < 3; j++) { - xyz_min[j] = -std::numeric_limits::max(); - xyz_max[j] = std::numeric_limits::max(); + xyz_min[j] = BF; + xyz_max[j] = -BF; } for (pid_t i = 0; i < pdata.np; i++) { @@ -138,6 +143,7 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) checkParticle(xyz, is_it_in_main, is_it_in_buf); if (is_it_in_main) { + assert(nvp < nvpall); for (int j = 0; j < 3; j++) { parts[3*nvp+j] = xyz[j]; @@ -160,6 +166,7 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) if (is_it_in_buf && !is_it_in_main) { + assert(nvpbuf < nvpall); for (int j = 0; j < 3; j++) { parts[3*nvpbuf+j] = xyz[j]; @@ -173,7 +180,6 @@ void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs) } cout << format("nvpbuf = %d") % nvpbuf << endl; cout << format("x: %f,%f; y: %f,%f; z:%f,%f\n") % xyz_min[0] % xyz_max[0] % xyz_min[1] % xyz_max[1] % xyz_min[2] % xyz_max[2] << endl; - nvpall = nvpbuf; double predict = 1; for (int j = 0; j < 3; j++) @@ -188,6 +194,8 @@ void BoxData::addGuardPoints() int rot_si[3][3] = { {1, 0, 0}, {1,0,0}, {0,1,0} }; int rot_sj[3][3] = { {0, 1, 0}, {0,0,1}, {0,0,1} }; + pid_t nvpcur = nvpbuf; + for (int k = 0; k < 3; k++) { @@ -196,17 +204,19 @@ void BoxData::addGuardPoints() { for (int j=0; j <= NGUARD; j++) { + assert(nvpcur < nvpall); /* Bottom */ for (int a = 0; a < 3; a++) - parts[3*nvpall+a] = -width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] - rot_g[k][a] * g[a]; + parts[3*nvpcur+a] = -width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] - rot_g[k][a] * g[a]; - nvpall++; + nvpcur++; + assert(nvpcur < nvpall); /* Top */ for (int a = 0; a < 3; a++) - parts[3*nvpall+a] = (2*rot_g[k][a]-1)*width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] + rot_g[k][a] * g[a]; + parts[3*nvpcur+a] = (2*rot_g[k][a]-1)*width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] + rot_g[k][a] * g[a]; - nvpall++; + nvpcur++; } } } @@ -218,15 +228,15 @@ void BoxData::findBoundary() for (int j = 0; j < 3; j++) { - xyz_min[j] = -BF; - xyz_max[j] = BF; + xyz_min[j] = BF; + xyz_max[j] = -BF; } for (pid_t i = nvpbuf; i < nvpall; i++) { for (int j = 0; j < 3; j++) { xyz_min[j] = std::min(xyz_min[j], parts[3*i+j]); - xyz_max[j] = std::min(xyz_max[j], parts[3*i+j]); + xyz_max[j] = std::max(xyz_max[j], parts[3*i+j]); } } } @@ -271,7 +281,6 @@ int main(int argc, char *argv[]) { char *suffix, *outDir, *posfile; PARTADJ *adjs; float *vols; - pid_t *orig; double border, boxsize[3]; int b[3], numdiv[3]; double totalvol; @@ -375,7 +384,7 @@ int main(int argc, char *argv[]) { /* Get the adjacencies back to their original values */ for (pid_t i=0; i Date: Tue, 11 Jun 2013 11:58:29 +0200 Subject: [PATCH 9/9] Updated launchers to follow the new location of jozov2 --- python_tools/void_python_tools/backend/launchers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python_tools/void_python_tools/backend/launchers.py b/python_tools/void_python_tools/backend/launchers.py index a4599f6..e53b241 100644 --- a/python_tools/void_python_tools/backend/launchers.py +++ b/python_tools/void_python_tools/backend/launchers.py @@ -276,7 +276,7 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None, cmd = "./%s >> %s 2>&1" % (vozScript, logFile) os.system(cmd) - cmd = "%s/../c_tools/zobov2/jozov2 %s %s %s %s %s %g %s >> %s 2>&1" % \ + cmd = "%s/../c_tools/zobov2/jozov2/jozov2 %s %s %s %s %s %g %s >> %s 2>&1" % \ (binPath, \ zobovDir+"/adj_"+sampleName+".dat", \ zobovDir+"/vol_"+sampleName+".dat", \