From e57d9a731aeca13399d865ad7560e010b94ba7c0 Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Wed, 26 Aug 2009 11:52:16 -0500 Subject: [PATCH] Initial import --- .gitignore | 7 + Makefile | 41 ++ findrtop.c | 81 ++++ jozov.c | 550 +++++++++++++++++++++ mytools/Makefile | 47 ++ mytools/config.mk | 6 + mytools/loadZobov.cpp | 136 ++++++ mytools/loadZobov.hpp | 28 ++ mytools/zobovConf.c | 1068 +++++++++++++++++++++++++++++++++++++++++ mytools/zobovConf.h | 243 ++++++++++ readfiles.c | 121 +++++ readme.txt | 25 + voz.h | 8 + voz1b1.c | 334 +++++++++++++ vozinit.c | 153 ++++++ voztie.c | 182 +++++++ vozutil.c | 241 ++++++++++ 17 files changed, 3271 insertions(+) create mode 100644 .gitignore create mode 100644 Makefile create mode 100644 findrtop.c create mode 100644 jozov.c create mode 100644 mytools/Makefile create mode 100644 mytools/config.mk create mode 100644 mytools/loadZobov.cpp create mode 100644 mytools/loadZobov.hpp create mode 100644 mytools/zobovConf.c create mode 100644 mytools/zobovConf.h create mode 100644 readfiles.c create mode 100644 readme.txt create mode 100644 voz.h create mode 100644 voz1b1.c create mode 100644 vozinit.c create mode 100644 voztie.c create mode 100644 vozutil.c diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..248a728 --- /dev/null +++ b/.gitignore @@ -0,0 +1,7 @@ +*~ +*.o +*.prog +voz1b1 +jozov +vozinit +voztie diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..b7e7c94 --- /dev/null +++ b/Makefile @@ -0,0 +1,41 @@ +# Compiler choice: + +# Gcc +CC = icc +CFLAGS = -O3 -ggdb + +MLIBS = -lm + +#QLIB = -L../qhull-2003.1/src/ +#QINC = -I../qhull-2003.1/src/ +QLIB = -L../qhull2002.1/src/ +QINC = -I../qhull2002.1/src/ + +############### + +all: vozinit voz1b1 voztie jozov + +jozov: jozov.o findrtop.o Makefile + $(CC) $(CFLAGS) -o jozov jozov.o findrtop.o $(MLIBS) +jozov.o: jozov.c Makefile + $(CC) $(CFLAGS) -c jozov.c +findrtop.o: findrtop.c Makefile + $(CC) $(CFLAGS) -c findrtop.c + +voz1b1: voz1b1.o readfiles.o vozutil.o voz.h Makefile + $(CC) -o voz1b1 $(CFLAGS) voz1b1.o readfiles.o vozutil.o -L. $(QLIB) -lqhull $(MLIBS) +voz1b1.o: voz1b1.c Makefile + $(CC) $(CFLAGS) $(QINC) -c voz1b1.c +vozutil.o: vozutil.c Makefile + $(CC) $(CFLAGS) $(QINC) -c vozutil.c + +vozinit: vozinit.o readfiles.o voz.h Makefile + $(CC) -o vozinit $(CFLAGS) vozinit.o readfiles.o -L. $(MLIBS) +vozinit.o: vozinit.c Makefile + $(CC) $(CFLAGS) -c vozinit.c $(QINC) + +voztie: voztie.o readfiles.o Makefile + $(CC) -o voztie $(CFLAGS) voztie.o readfiles.o +voztie.o: voztie.c Makefile + $(CC) $(CFLAGS) -c voztie.c + diff --git a/findrtop.c b/findrtop.c new file mode 100644 index 0000000..8dbeeb3 --- /dev/null +++ b/findrtop.c @@ -0,0 +1,81 @@ +#include +#include + +/*------------------------------------------------------------------------------ + Find nb elements of Real array a having the largest value. + Returns index iord of these elements, ordered so iord[0] corresponds + to element a[iord[0]] having the largest value. + If nb > na, then the last nb-na elements of iord are undefined. + Elements of a that are equal are left in their original order. + + Courtesy Andrew Hamilton. +*/ +void findrtop(double *a, int na, int *iord, int nb) +{ +#undef ORDER +#define ORDER(ia, ja) a[ia] > a[ja] || (a[ia] == a[ja] && ia < ja) + + int i, ia, ib, it, n; + + n = (na <= nb)? na : nb; + if (n <= 0) return; + + /* heap first n elements, so smallest element is at top of heap */ + for (ib = (n >> 1); ib < n; ib++) { + iord[ib] = ib; + } + for (ia = (n >> 1) - 1; ia >= 0; ia--) { + i = ia; + for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) { + if (ib+1 < n) { + if (ORDER(iord[ib], iord[ib+1])) ib++; + } + if (ORDER(ia, iord[ib])) { + iord[i] = iord[ib]; + i = ib; + } else { + break; + } + } + iord[i] = ia; + } + + /* now compare rest of elements of array to heap */ + for (ia = n; ia < na; ia++) { + /* if new element is greater than smallest, sift it into heap */ + i = 0; + if (ORDER(ia, iord[i])) { + for (ib = (i << 1) + 1; ib < n; ib = (i << 1) + 1) { + if (ib+1 < n) { + if (ORDER(iord[ib], iord[ib+1])) ib++; + } + if (ORDER(ia, iord[ib])) { + iord[i] = iord[ib]; + i = ib; + } else { + break; + } + } + iord[i] = ia; + } + } + + /* unheap iord so largest element is at top */ + for (ia = n - 1; ia > 0; ia--) { + it = iord[ia]; + i = 0; + iord[ia] = iord[i]; + for (ib = (i << 1) + 1; ib < ia; ib = (i << 1) + 1) { + if (ib+1 < ia) { + if (ORDER(iord[ib], iord[ib+1])) ib++; + } + if (ORDER(it, iord[ib])) { + iord[i] = iord[ib]; + i = ib; + } else { + break; + } + } + iord[i] = it; + } +} diff --git a/jozov.c b/jozov.c new file mode 100644 index 0000000..f5fb4e7 --- /dev/null +++ b/jozov.c @@ -0,0 +1,550 @@ +#include +/* jovoz.c by Mark Neyrinck */ +#include +#include +#include + +#define BIGFLT 1e30 /* Biggest possible floating-point number */ +#define NLINKS 100 /* Number of possible links with the same rho_sl */ +#define FF fflush(stdout) + +typedef struct Particle { + float dens; + int nadj; + int ncnt; + int *adj; +} PARTICLE; + +typedef struct Zone { + int core; /* Identity of peak particle */ + int np; /* Number of particles in zone */ + int npjoin; /* Number of particles in the joined void */ + int nadj; /* Number of adjacent zones */ + int nhl; /* Number of zones in final joined void */ + float leak; /* Volume of last leak zone*/ + int *adj; /* Each adjacent zone, with ... */ + float *slv; /* Smallest Linking Volume */ + float denscontrast; /* density contrast */ + double vol; /* Total volume of all particles in the zone */ + double voljoin; /* Total volume of all particles in the joined void */ +} ZONE; + +typedef struct ZoneT { + int nadj; /* Number of zones on border */ + int *adj; /* Each adjacent zone, with ... */ + float *slv; /* Smallest Linking Volume */ +} ZONET; + +void findrtop(double *a, int na, int *iord, int nb); + +int main(int argc,char **argv) { + + FILE *adj, *vol, *zon, *zon2, *txt; + PARTICLE *p; + ZONE *z; + ZONET *zt; + char *adjfile, *volfile, *zonfile, *zonfile2, *txtfile; + int i, j,k,l, h, h2,hl,n,np, np2, nzones, nhl, nhlcount, nhl2; + int *jumper, *jumped, *numinh; + int *zonenum, *zonelist, *zonelist2; + int link[NLINKS], link2, nl; + float lowvol, voltol, prob; + + int q,q2; + + int za, nin; + int testpart; + char already, interior, *inyet, *inyet2, added, beaten; + int *nm, **m; + + float maxvol, minvol; + double *sorter, e1,maxdenscontrast; + int *iord; + + e1 = exp(1.)-1.; + + if (argc != 7) { + printf("Wrong number of arguments.\n"); + printf("arg1: adjacency file\n"); + printf("arg2: volume file\n"); + printf("arg3: output file containing particles in each zone\n"); + printf("arg4: output file containing zones in each void\n"); + printf("arg5: output text file\n"); + printf("arg6: Density threshold (0 for no threshold)\n\n"); + exit(0); + } + adjfile = argv[1]; + volfile = argv[2]; + zonfile = argv[3]; + zonfile2 = argv[4]; + txtfile = argv[5]; + if (sscanf(argv[6],"%f",&voltol) == 0) { + printf("Bad density threshold.\n"); + exit(0); + } + if (voltol <= 0.) { + printf("Proceeding without a density threshold.\n"); + voltol = 1e30; + } + + adj = fopen(adjfile, "r"); + if (adj == NULL) { + printf("Unable to open %s\n",adjfile); + exit(0); + } + fread(&np,1, sizeof(int),adj); + + printf("adj: %d particles\n", np); + FF; + + p = (PARTICLE *)malloc(np*sizeof(PARTICLE)); + /* Adjacencies*/ + for (i=0;i 0) + p[i].adj = (int *)malloc(p[i].nadj*sizeof(int)); + else p[i].adj = 0; + p[i].ncnt = 0; /* Temporarily, it's an adj counter */ + } + for (i=0;i 0) + for (k=0;k 1e30)) { + printf("Whacked-out volume found, of particle %d: %f\n",i,p[i].dens); + p[i].dens = 1.; + } + p[i].dens = 1./p[i].dens; /* Get density from volume */ + } + fclose(vol); + + jumped = (int *)malloc(np*sizeof(int)); + jumper = (int *)malloc(np*sizeof(int)); + numinh = (int *)malloc(np*sizeof(int)); + + /* find jumper */ + for (i = 0; i < np; i++) { + minvol = p[i].dens; jumper[i] = -1; + for (j=0; j< p[i].nadj; j++) { + if (p[p[i].adj[j]].dens < minvol) { + jumper[i] = p[i].adj[j]; + minvol = p[jumper[i]].dens; + } + } + numinh[i] = 0; + } + + printf("About to jump ...\n"); FF; + + /* Jump */ + for (i = 0; i < np; i++) { + jumped[i] = i; + while (jumper[jumped[i]] > -1) + jumped[i] = jumper[jumped[i]]; + numinh[jumped[i]]++; + } + printf("Post-jump ...\n"); FF; + + nzones = 0; + for (i = 0; i < np; i++) + if (numinh[i] > 0) nzones++; + printf("%d initial zones found\n", nzones); + + z = (ZONE *)malloc(nzones*sizeof(ZONE)); + if (z == NULL) { + printf("Unable to allocate z\n"); + exit(0); + } + zt = (ZONET *)malloc(nzones*sizeof(ZONET)); + if (zt == NULL) { + printf("Unable to allocate zt\n"); + exit(0); + } + zonenum = (int *)malloc(np*sizeof(int)); + if (zonenum == NULL) { + printf("Unable to allocate zonenum\n"); + exit(0); + } + + h = 0; + for (i = 0; i < np; i++) + if (numinh[i] > 0) { + z[h].core = i; + zonenum[i] = h; + h++; + } else { + zonenum[i] = -1; + } + + /* Count border particles */ + for (i = 0; i < np; i++) + for (j = 0; j < p[i].nadj; j++) { + testpart = p[i].adj[j]; + if (jumped[i] != jumped[testpart]) + zt[zonenum[jumped[i]]].nadj++; + } + + for (h=0;h p[i].dens) { + /* there could be a weakest link through testpart */ + already = 0; + for (za = 0; za < zt[h].nadj; za++) + if (zt[h].adj[za] == zonenum[jumped[testpart]]) { + already = 1; + if (p[testpart].dens < zt[h].slv[za]) { + zt[h].slv[za] = p[testpart].dens; + } + } + if (already == 0) { + zt[h].adj[zt[h].nadj] = zonenum[jumped[testpart]]; + zt[h].slv[zt[h].nadj] = p[testpart].dens; + zt[h].nadj++; + } + } else { /* There could be a weakest link through i */ + already = 0; + for (za = 0; za < zt[h].nadj; za++) + if (zt[h].adj[za] == zonenum[jumped[testpart]]) { + already = 1; + if (p[i].dens < zt[h].slv[za]) { + zt[h].slv[za] = p[i].dens; + } + } + if (already == 0) { + zt[h].adj[zt[h].nadj] = zonenum[jumped[testpart]]; + zt[h].slv[zt[h].nadj] = p[i].dens; + zt[h].nadj++; + } + } + } + } + } + printf("Found zone adjacencies\n"); FF; + + /* Free particle adjacencies */ + for (i=0;i maxvol) maxvol = p[i].dens; + if (p[i].dens < minvol) minvol = p[i].dens; + } + printf("Densities range from %e to %e.\n",minvol,maxvol);FF; + + zon2 = fopen(zonfile2,"w"); + if (zon2 == NULL) { + printf("Problem opening zonefile %s.\n\n",zonfile2); + exit(0); + } + fwrite(&nzones,1,4,zon2); + + for (h = 0; h voltol) { + beaten = 1; + z[h].leak = lowvol; + continue; + } + + for (l=0; l < nl; l++) + if (p[z[link[l]].core].dens < p[z[h].core].dens) + beaten = 1; + if (beaten == 1) { + z[h].leak = lowvol; + continue; + } + /* Add everything linked to the link(s) */ + nhl2 = 0; + for (l=0; l < nl; l++) { + if (inyet2[link[l]] == 0) { + zonelist2[nhl2] = link[l]; + inyet2[link[l]] = 1; + nhl2 ++; + added = 1; + while ((added == 1) && (beaten == 0)) { + added = 0; + for (hl = 0; (hl < nhl2) && (beaten == 0); hl++) { + h2 = zonelist2[hl]; + if (inyet2[h2] == 1) { + interior = 1; /* Guilty until proven innocent */ + for (za = 0; za < z[h2].nadj; za ++) { + link2 = z[h2].adj[za]; + if ((inyet[link2]+inyet2[link2]) == 0) { + interior = 0; + if (z[h2].slv[za] <= lowvol) { + if (p[z[link2].core].dens < p[z[h].core].dens) { + beaten = 1; + break; + } + zonelist2[nhl2] = link2; + inyet2[link2] = 1; + nhl2++; + added = 1; + } + } + } + if (interior == 1) inyet2[h2] = 2; + } + } + } + } + } + for (hl = 0; hl < nhl2; hl++) + inyet2[zonelist2[hl]] = 0; + + /* See if there's a beater */ + if (beaten == 1) { + z[h].leak = lowvol; + } else { + for (h2 = 0; h2 < nhl2; h2++) { + zonelist[nhl] = zonelist2[h2]; + inyet[zonelist2[h2]] = 1; + nhl++; + z[h].npjoin += z[zonelist2[h2]].np; + } + } + if (nhl/10000 > nhlcount) { + nhlcount = nhl/10000; + printf(" %d",nhl); FF; + } + } while((lowvol < BIGFLT) && (beaten == 0)); + + z[h].denscontrast = z[h].leak/p[z[h].core].dens; + if (z[h].denscontrast < 1.) z[h].denscontrast = 1.; + + /* find biggest denscontrast */ + if (z[h].denscontrast > maxdenscontrast) { + maxdenscontrast = (double)z[h].denscontrast; + } + + /* Don't sort; want the core zone to be first */ + + if (nhlcount > 0) { /* Outputs the number of zones in large voids */ + printf(" h%d:%d\n",h,nhl); + FF; + } + /* Calculate volume */ + z[h].voljoin = 0.; + for (q = 0; q .mydepends + +distclean: clean + @rm -f .mydepends + +clean: + @rm -f *.o + @rm -f $(PROGS) + +.mydepends: $(SOURCES) + @touch .mydepends + @make depend + +%.prog: + @echo "[L] $@" + @$(CXX) -o $@ $^ $(LIBS) + +zobovConf.c zobovConf.h: showZobov.ggo Makefile + @echo "[OPT] $@" + gengetopt -i $< -f zobovConf -a zobovConf_info -F zobovConf -C + +%.o: %.c + @echo "[C] $< ..." + @$(CC) -c -o $@ $< $(CFLAGS) + +%.o: %.cpp + @echo "[C++] $< ..." + @$(CXX) -c -o $@ $< $(CXXFLAGS) + +include .mydepends diff --git a/mytools/config.mk b/mytools/config.mk new file mode 100644 index 0000000..d41d32c --- /dev/null +++ b/mytools/config.mk @@ -0,0 +1,6 @@ +CC= gcc +CXX= g++ +CPPFLAGS= -I$(HOME)/Science/Software/CosmoToolbox/install/include +LDFLAGS= -L$(HOME)/Science/Software/CosmoToolbox/install/lib -lCosmoTool +CXXFLAGS= $(CPPFLAGS) -ggdb -O0 -ffast-math +CFLAGS= $(CPPFLAGS) -ggdb -O0 -ffast-math diff --git a/mytools/loadZobov.cpp b/mytools/loadZobov.cpp new file mode 100644 index 0000000..521c0a8 --- /dev/null +++ b/mytools/loadZobov.cpp @@ -0,0 +1,136 @@ +#include +#include +#include +#include +#include +#include +#include "loadZobov.hpp" + +using namespace std; + +bool loadZobov(const char *descName, const char *adjName, const char *volName, ZobovRep& z) +{ + ifstream descFile(descName); + ifstream adjFile(adjName); + ifstream volFile(volName); + int32_t numParticles, numZones, numPinZone; + int32_t totalParticles; + int32_t numVoids; + + adjFile.read((char *)&numParticles, sizeof(numParticles)); + adjFile.read((char *)&numZones, sizeof(numZones)); + if (!adjFile) + return false; + + cout << "Number of particles = " << numParticles << endl; + cout << "Number of zones = " << numZones << endl; + + totalParticles = 0; + z.allZones.resize(numZones); + for (int zone = 0; zone < numZones; zone++) + { + adjFile.read((char *)&numPinZone, sizeof(numPinZone)); + if (!adjFile) + { + cout << "Problem on the zone " << zone << " / " << numZones << endl; + return false; + } + + z.allZones[zone].pId.resize(numPinZone); + adjFile.read((char *)&z.allZones[zone].pId[0], sizeof(int)*numPinZone); + + totalParticles += numPinZone; + } + cout << "Zoned " << totalParticles << endl; + + if (totalParticles != numParticles) + { + cerr << "The numbers of particles are inconsistent ! (" << totalParticles << " vs " << numParticles << ")"<< endl; + abort(); + } + + volFile.read((char *)&numVoids, sizeof(numVoids)); + if (!volFile) + return false; + + cout << "Number of voids = " << numVoids << endl; + + z.allVoids.resize(numVoids); + for (int v = 0; v < numVoids; v++) + { + int32_t numZinV; + + volFile.read((char *)&numZinV, sizeof(numZinV)); + if (!volFile) + return false; + + z.allVoids[v].zId.resize(numZinV); + + int *zId = new int[numZinV]; + + volFile.read((char *)zId, sizeof(int)*numZinV); + for (int k = 0; k < numZinV; k++) + z.allVoids[v].zId[k] = &z.allZones[zId[k]]; + + delete[] zId; + } + + cout << "Loading description" << endl; + + string line; + getline(descFile, line); + getline(descFile, line); + getline(descFile, line); + while (!descFile.eof()) + { + istringstream lineStream(line.c_str()); + int orderId, volId, coreParticle, numParticlesInZone, numZonesInVoid, numInVoid; + float coreDensity, volumeZone, volumeVoid, densityContrast; + float probability; + + lineStream + >> orderId + >> volId + >> coreParticle + >> coreDensity + >> volumeZone + >> numParticlesInZone + >> numZonesInVoid + >> volumeVoid + >> numInVoid + >> densityContrast + >> probability; + if (!lineStream) + { + cerr << "Error in text stream" << endl; + abort(); + } + + z.allVoids[volId].proba = probability; + z.allVoids[volId].volume = volumeVoid; + z.allVoids[volId].numParticles = numInVoid; + z.allVoids[volId].coreParticle = coreParticle; + + // Sanity check + int actualNumber = 0; + for (int j = 0; j < z.allVoids[volId].zId.size(); j++) + actualNumber += z.allVoids[volId].zId[j]->pId.size(); + + if (actualNumber != numInVoid) + { + cerr << "Sanity check failed." + << " The number of particles in the description (" + << numInVoid + << ") is different from the one in the file (" + << actualNumber << ")" << endl; + } + getline(descFile, line); + } + + cout << "Done loading" << endl; + + + + + return true; +} diff --git a/mytools/loadZobov.hpp b/mytools/loadZobov.hpp new file mode 100644 index 0000000..f660558 --- /dev/null +++ b/mytools/loadZobov.hpp @@ -0,0 +1,28 @@ +#ifndef __LOAD_ZOBOV_HPP +#define __LOAD_ZOBOV_HPP + +#include + +struct ZobovZone +{ + std::vector pId; +}; + +struct ZobovVoid +{ + std::vector zId; + float proba; + int numParticles, coreParticle; + float volume; +}; + +struct ZobovRep +{ + std::vector allZones; + std::vector allVoids; +}; + +bool loadZobov(const char *descName, + const char *adjName, const char *volName, ZobovRep& z); + +#endif diff --git a/mytools/zobovConf.c b/mytools/zobovConf.c new file mode 100644 index 0000000..3019686 --- /dev/null +++ b/mytools/zobovConf.c @@ -0,0 +1,1068 @@ +/* + File autogenerated by gengetopt version 2.22 + generated with the following command: + gengetopt -i showZobov.ggo -f zobovConf -a zobovConf_info -F zobovConf -C + + The developers of gengetopt consider the fixed text that goes in all + gengetopt output files to be in the public domain: + we make no copyright claims on it. +*/ + +/* If we use autoconf. */ +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include + +#include "getopt.h" + +#include "zobovConf.h" + +const char *zobovConf_info_purpose = ""; + +const char *zobovConf_info_usage = "Usage: showZobov [OPTIONS]..."; + +const char *zobovConf_info_description = ""; + +const char *zobovConf_info_help[] = { + " -h, --help Print help and exit", + " -V, --version Print version and exit", + " -d, --desc=STRING The description file name for the voids", + " -a, --adj=STRING Adjacent file name", + " -v, --void=STRING Void/zone bind filename", + " -i, --id=INT Void id to extract", + " -s, --sProba Sort against probability (default=off)", + " -l, --sVol Sort against volume (default=on)", + " -m, --mProba=DOUBLE Minimal probability to accept (default=`0.0')", + " -b, --ramsesDir=STRING Ramses base output directory", + " -r, --ramsesId=INT Ramses output id", + " -f, --configFile=STRING Configuration file", + " -q, --quiet Quiet mode (default=off)", + " -p, --move=STRING Move the center by (x,y,z) (default=`(0,0,0)')", + " --galax Output is galax particles (default=off)", + " -o, --output=STRING Output file (default=`voidDesc')", + 0 +}; + +typedef enum {ARG_NO + , ARG_FLAG + , ARG_STRING + , ARG_INT + , ARG_DOUBLE +} zobovConf_arg_type; + +static +void clear_given (struct zobovConf_info *args_info); +static +void clear_args (struct zobovConf_info *args_info); + +static int +zobovConf_internal (int argc, char * const *argv, struct zobovConf_info *args_info, + struct zobovConf_params *params, const char *additional_error); + +static int +zobovConf_required2 (struct zobovConf_info *args_info, const char *prog_name, const char *additional_error); +struct line_list +{ + char * string_arg; + struct line_list * next; +}; + +static struct line_list *cmd_line_list = 0; +static struct line_list *cmd_line_list_tmp = 0; + +static void +free_cmd_list(void) +{ + /* free the list of a previous call */ + if (cmd_line_list) + { + while (cmd_line_list) { + cmd_line_list_tmp = cmd_line_list; + cmd_line_list = cmd_line_list->next; + free (cmd_line_list_tmp->string_arg); + free (cmd_line_list_tmp); + } + } +} + + +static char * +gengetopt_strdup (const char *s); + +static +void clear_given (struct zobovConf_info *args_info) +{ + args_info->help_given = 0 ; + args_info->version_given = 0 ; + args_info->desc_given = 0 ; + args_info->adj_given = 0 ; + args_info->void_given = 0 ; + args_info->id_given = 0 ; + args_info->sProba_given = 0 ; + args_info->sVol_given = 0 ; + args_info->mProba_given = 0 ; + args_info->ramsesDir_given = 0 ; + args_info->ramsesId_given = 0 ; + args_info->configFile_given = 0 ; + args_info->quiet_given = 0 ; + args_info->move_given = 0 ; + args_info->galax_given = 0 ; + args_info->output_given = 0 ; +} + +static +void clear_args (struct zobovConf_info *args_info) +{ + args_info->desc_arg = NULL; + args_info->desc_orig = NULL; + args_info->adj_arg = NULL; + args_info->adj_orig = NULL; + args_info->void_arg = NULL; + args_info->void_orig = NULL; + args_info->id_orig = NULL; + args_info->sProba_flag = 0; + args_info->sVol_flag = 1; + args_info->mProba_arg = 0.0; + args_info->mProba_orig = NULL; + args_info->ramsesDir_arg = NULL; + args_info->ramsesDir_orig = NULL; + args_info->ramsesId_orig = NULL; + args_info->configFile_arg = NULL; + args_info->configFile_orig = NULL; + args_info->quiet_flag = 0; + args_info->move_arg = gengetopt_strdup ("(0,0,0)"); + args_info->move_orig = NULL; + args_info->galax_flag = 0; + args_info->output_arg = gengetopt_strdup ("voidDesc"); + args_info->output_orig = NULL; + +} + +static +void init_args_info(struct zobovConf_info *args_info) +{ + + + args_info->help_help = zobovConf_info_help[0] ; + args_info->version_help = zobovConf_info_help[1] ; + args_info->desc_help = zobovConf_info_help[2] ; + args_info->adj_help = zobovConf_info_help[3] ; + args_info->void_help = zobovConf_info_help[4] ; + args_info->id_help = zobovConf_info_help[5] ; + args_info->sProba_help = zobovConf_info_help[6] ; + args_info->sVol_help = zobovConf_info_help[7] ; + args_info->mProba_help = zobovConf_info_help[8] ; + args_info->ramsesDir_help = zobovConf_info_help[9] ; + args_info->ramsesId_help = zobovConf_info_help[10] ; + args_info->configFile_help = zobovConf_info_help[11] ; + args_info->quiet_help = zobovConf_info_help[12] ; + args_info->move_help = zobovConf_info_help[13] ; + args_info->galax_help = zobovConf_info_help[14] ; + args_info->output_help = zobovConf_info_help[15] ; + +} + +void +zobovConf_print_version (void) +{ + printf ("%s %s\n", ZOBOVCONF_PACKAGE, ZOBOVCONF_VERSION); +} + +static void print_help_common(void) { + zobovConf_print_version (); + + if (strlen(zobovConf_info_purpose) > 0) + printf("\n%s\n", zobovConf_info_purpose); + + if (strlen(zobovConf_info_usage) > 0) + printf("\n%s\n", zobovConf_info_usage); + + printf("\n"); + + if (strlen(zobovConf_info_description) > 0) + printf("%s\n", zobovConf_info_description); +} + +void +zobovConf_print_help (void) +{ + int i = 0; + print_help_common(); + while (zobovConf_info_help[i]) + printf("%s\n", zobovConf_info_help[i++]); +} + +void +zobovConf_init (struct zobovConf_info *args_info) +{ + clear_given (args_info); + clear_args (args_info); + init_args_info (args_info); +} + +void +zobovConf_params_init(struct zobovConf_params *params) +{ + if (params) + { + params->override = 0; + params->initialize = 1; + params->check_required = 1; + params->check_ambiguity = 0; + params->print_errors = 1; + } +} + +struct zobovConf_params * +zobovConf_params_create(void) +{ + struct zobovConf_params *params = + (struct zobovConf_params *)malloc(sizeof(struct zobovConf_params)); + zobovConf_params_init(params); + return params; +} + +static void +free_string_field (char **s) +{ + if (*s) + { + free (*s); + *s = 0; + } +} + + +static void +zobovConf_release (struct zobovConf_info *args_info) +{ + + free_string_field (&(args_info->desc_arg)); + free_string_field (&(args_info->desc_orig)); + free_string_field (&(args_info->adj_arg)); + free_string_field (&(args_info->adj_orig)); + free_string_field (&(args_info->void_arg)); + free_string_field (&(args_info->void_orig)); + free_string_field (&(args_info->id_orig)); + free_string_field (&(args_info->mProba_orig)); + free_string_field (&(args_info->ramsesDir_arg)); + free_string_field (&(args_info->ramsesDir_orig)); + free_string_field (&(args_info->ramsesId_orig)); + free_string_field (&(args_info->configFile_arg)); + free_string_field (&(args_info->configFile_orig)); + free_string_field (&(args_info->move_arg)); + free_string_field (&(args_info->move_orig)); + free_string_field (&(args_info->output_arg)); + free_string_field (&(args_info->output_orig)); + + + + clear_given (args_info); +} + + +static void +write_into_file(FILE *outfile, const char *opt, const char *arg, char *values[]) +{ + if (arg) { + fprintf(outfile, "%s=\"%s\"\n", opt, arg); + } else { + fprintf(outfile, "%s\n", opt); + } +} + + +int +zobovConf_dump(FILE *outfile, struct zobovConf_info *args_info) +{ + int i = 0; + + if (!outfile) + { + fprintf (stderr, "%s: cannot dump options to stream\n", ZOBOVCONF_PACKAGE); + return EXIT_FAILURE; + } + + if (args_info->help_given) + write_into_file(outfile, "help", 0, 0 ); + if (args_info->version_given) + write_into_file(outfile, "version", 0, 0 ); + if (args_info->desc_given) + write_into_file(outfile, "desc", args_info->desc_orig, 0); + if (args_info->adj_given) + write_into_file(outfile, "adj", args_info->adj_orig, 0); + if (args_info->void_given) + write_into_file(outfile, "void", args_info->void_orig, 0); + if (args_info->id_given) + write_into_file(outfile, "id", args_info->id_orig, 0); + if (args_info->sProba_given) + write_into_file(outfile, "sProba", 0, 0 ); + if (args_info->sVol_given) + write_into_file(outfile, "sVol", 0, 0 ); + if (args_info->mProba_given) + write_into_file(outfile, "mProba", args_info->mProba_orig, 0); + if (args_info->ramsesDir_given) + write_into_file(outfile, "ramsesDir", args_info->ramsesDir_orig, 0); + if (args_info->ramsesId_given) + write_into_file(outfile, "ramsesId", args_info->ramsesId_orig, 0); + if (args_info->configFile_given) + write_into_file(outfile, "configFile", args_info->configFile_orig, 0); + if (args_info->quiet_given) + write_into_file(outfile, "quiet", 0, 0 ); + if (args_info->move_given) + write_into_file(outfile, "move", args_info->move_orig, 0); + if (args_info->galax_given) + write_into_file(outfile, "galax", 0, 0 ); + if (args_info->output_given) + write_into_file(outfile, "output", args_info->output_orig, 0); + + + i = EXIT_SUCCESS; + return i; +} + +int +zobovConf_file_save(const char *filename, struct zobovConf_info *args_info) +{ + FILE *outfile; + int i = 0; + + outfile = fopen(filename, "w"); + + if (!outfile) + { + fprintf (stderr, "%s: cannot open file for writing: %s\n", ZOBOVCONF_PACKAGE, filename); + return EXIT_FAILURE; + } + + i = zobovConf_dump(outfile, args_info); + fclose (outfile); + + return i; +} + +void +zobovConf_free (struct zobovConf_info *args_info) +{ + zobovConf_release (args_info); +} + +/** @brief replacement of strdup, which is not standard */ +char * +gengetopt_strdup (const char *s) +{ + char *result = NULL; + if (!s) + return result; + + result = (char*)malloc(strlen(s) + 1); + if (result == (char*)0) + return (char*)0; + strcpy(result, s); + return result; +} + +int +zobovConf (int argc, char * const *argv, struct zobovConf_info *args_info) +{ + return zobovConf2 (argc, argv, args_info, 0, 1, 1); +} + +int +zobovConf_ext (int argc, char * const *argv, struct zobovConf_info *args_info, + struct zobovConf_params *params) +{ + int result; + result = zobovConf_internal (argc, argv, args_info, params, NULL); + + if (result == EXIT_FAILURE) + { + zobovConf_free (args_info); + exit (EXIT_FAILURE); + } + + return result; +} + +int +zobovConf2 (int argc, char * const *argv, struct zobovConf_info *args_info, int override, int initialize, int check_required) +{ + int result; + struct zobovConf_params params; + + params.override = override; + params.initialize = initialize; + params.check_required = check_required; + params.check_ambiguity = 0; + params.print_errors = 1; + + result = zobovConf_internal (argc, argv, args_info, ¶ms, NULL); + + if (result == EXIT_FAILURE) + { + zobovConf_free (args_info); + exit (EXIT_FAILURE); + } + + return result; +} + +int +zobovConf_required (struct zobovConf_info *args_info, const char *prog_name) +{ + int result = EXIT_SUCCESS; + + if (zobovConf_required2(args_info, prog_name, NULL) > 0) + result = EXIT_FAILURE; + + if (result == EXIT_FAILURE) + { + zobovConf_free (args_info); + exit (EXIT_FAILURE); + } + + return result; +} + +int +zobovConf_required2 (struct zobovConf_info *args_info, const char *prog_name, const char *additional_error) +{ + int error = 0; + + /* checks for required options */ + if (! args_info->desc_given) + { + fprintf (stderr, "%s: '--desc' ('-d') option required%s\n", prog_name, (additional_error ? additional_error : "")); + error = 1; + } + + if (! args_info->adj_given) + { + fprintf (stderr, "%s: '--adj' ('-a') option required%s\n", prog_name, (additional_error ? additional_error : "")); + error = 1; + } + + if (! args_info->void_given) + { + fprintf (stderr, "%s: '--void' ('-v') option required%s\n", prog_name, (additional_error ? additional_error : "")); + error = 1; + } + + if (! args_info->ramsesDir_given) + { + fprintf (stderr, "%s: '--ramsesDir' ('-b') option required%s\n", prog_name, (additional_error ? additional_error : "")); + error = 1; + } + + if (! args_info->ramsesId_given) + { + fprintf (stderr, "%s: '--ramsesId' ('-r') option required%s\n", prog_name, (additional_error ? additional_error : "")); + error = 1; + } + + + /* checks for dependences among options */ + + return error; +} + + +static char *package_name = 0; + +/** + * @brief updates an option + * @param field the generic pointer to the field to update + * @param orig_field the pointer to the orig field + * @param field_given the pointer to the number of occurrence of this option + * @param prev_given the pointer to the number of occurrence already seen + * @param value the argument for this option (if null no arg was specified) + * @param possible_values the possible values for this option (if specified) + * @param default_value the default value (in case the option only accepts fixed values) + * @param arg_type the type of this option + * @param check_ambiguity @see zobovConf_params.check_ambiguity + * @param override @see zobovConf_params.override + * @param no_free whether to free a possible previous value + * @param multiple_option whether this is a multiple option + * @param long_opt the corresponding long option + * @param short_opt the corresponding short option (or '-' if none) + * @param additional_error possible further error specification + */ +static +int update_arg(void *field, char **orig_field, + unsigned int *field_given, unsigned int *prev_given, + char *value, char *possible_values[], const char *default_value, + zobovConf_arg_type arg_type, + int check_ambiguity, int override, + int no_free, int multiple_option, + const char *long_opt, char short_opt, + const char *additional_error) +{ + char *stop_char = 0; + const char *val = value; + int found; + char **string_field; + + stop_char = 0; + found = 0; + + if (!multiple_option && prev_given && (*prev_given || (check_ambiguity && *field_given))) + { + if (short_opt != '-') + fprintf (stderr, "%s: `--%s' (`-%c') option given more than once%s\n", + package_name, long_opt, short_opt, + (additional_error ? additional_error : "")); + else + fprintf (stderr, "%s: `--%s' option given more than once%s\n", + package_name, long_opt, + (additional_error ? additional_error : "")); + return 1; /* failure */ + } + + + if (field_given && *field_given && ! override) + return 0; + if (prev_given) + (*prev_given)++; + if (field_given) + (*field_given)++; + if (possible_values) + val = possible_values[found]; + + switch(arg_type) { + case ARG_FLAG: + *((int *)field) = !*((int *)field); + break; + case ARG_INT: + if (val) *((int *)field) = strtol (val, &stop_char, 0); + break; + case ARG_DOUBLE: + if (val) *((double *)field) = strtod (val, &stop_char); + break; + case ARG_STRING: + if (val) { + string_field = (char **)field; + if (!no_free && *string_field) + free (*string_field); /* free previous string */ + *string_field = gengetopt_strdup (val); + } + break; + default: + break; + }; + + /* check numeric conversion */ + switch(arg_type) { + case ARG_INT: + case ARG_DOUBLE: + if (val && !(stop_char && *stop_char == '\0')) { + fprintf(stderr, "%s: invalid numeric value: %s\n", package_name, val); + return 1; /* failure */ + } + break; + default: + ; + }; + + /* store the original value */ + switch(arg_type) { + case ARG_NO: + case ARG_FLAG: + break; + default: + if (value && orig_field) { + if (no_free) { + *orig_field = value; + } else { + if (*orig_field) + free (*orig_field); /* free previous string */ + *orig_field = gengetopt_strdup (value); + } + } + }; + + return 0; /* OK */ +} + + +int +zobovConf_internal (int argc, char * const *argv, struct zobovConf_info *args_info, + struct zobovConf_params *params, const char *additional_error) +{ + int c; /* Character of the parsed option. */ + + int error = 0; + struct zobovConf_info local_args_info; + + int override; + int initialize; + int check_required; + int check_ambiguity; + + package_name = argv[0]; + + override = params->override; + initialize = params->initialize; + check_required = params->check_required; + check_ambiguity = params->check_ambiguity; + + if (initialize) + zobovConf_init (args_info); + + zobovConf_init (&local_args_info); + + optarg = 0; + optind = 0; + opterr = params->print_errors; + optopt = '?'; + + while (1) + { + int option_index = 0; + + static struct option long_options[] = { + { "help", 0, NULL, 'h' }, + { "version", 0, NULL, 'V' }, + { "desc", 1, NULL, 'd' }, + { "adj", 1, NULL, 'a' }, + { "void", 1, NULL, 'v' }, + { "id", 1, NULL, 'i' }, + { "sProba", 0, NULL, 's' }, + { "sVol", 0, NULL, 'l' }, + { "mProba", 1, NULL, 'm' }, + { "ramsesDir", 1, NULL, 'b' }, + { "ramsesId", 1, NULL, 'r' }, + { "configFile", 1, NULL, 'f' }, + { "quiet", 0, NULL, 'q' }, + { "move", 1, NULL, 'p' }, + { "galax", 0, NULL, 0 }, + { "output", 1, NULL, 'o' }, + { NULL, 0, NULL, 0 } + }; + + c = getopt_long (argc, argv, "hVd:a:v:i:slm:b:r:f:qp:o:", long_options, &option_index); + + if (c == -1) break; /* Exit from `while (1)' loop. */ + + switch (c) + { + case 'h': /* Print help and exit. */ + zobovConf_print_help (); + zobovConf_free (&local_args_info); + exit (EXIT_SUCCESS); + + case 'V': /* Print version and exit. */ + zobovConf_print_version (); + zobovConf_free (&local_args_info); + exit (EXIT_SUCCESS); + + case 'd': /* The description file name for the voids. */ + + + if (update_arg( (void *)&(args_info->desc_arg), + &(args_info->desc_orig), &(args_info->desc_given), + &(local_args_info.desc_given), optarg, 0, 0, ARG_STRING, + check_ambiguity, override, 0, 0, + "desc", 'd', + additional_error)) + goto failure; + + break; + case 'a': /* Adjacent file name. */ + + + if (update_arg( (void *)&(args_info->adj_arg), + &(args_info->adj_orig), &(args_info->adj_given), + &(local_args_info.adj_given), optarg, 0, 0, ARG_STRING, + check_ambiguity, override, 0, 0, + "adj", 'a', + additional_error)) + goto failure; + + break; + case 'v': /* Void/zone bind filename. */ + + + if (update_arg( (void *)&(args_info->void_arg), + &(args_info->void_orig), &(args_info->void_given), + &(local_args_info.void_given), optarg, 0, 0, ARG_STRING, + check_ambiguity, override, 0, 0, + "void", 'v', + additional_error)) + goto failure; + + break; + case 'i': /* Void id to extract. */ + + + if (update_arg( (void *)&(args_info->id_arg), + &(args_info->id_orig), &(args_info->id_given), + &(local_args_info.id_given), optarg, 0, 0, ARG_INT, + check_ambiguity, override, 0, 0, + "id", 'i', + additional_error)) + goto failure; + + break; + case 's': /* Sort against probability. */ + + + if (update_arg((void *)&(args_info->sProba_flag), 0, &(args_info->sProba_given), + &(local_args_info.sProba_given), optarg, 0, 0, ARG_FLAG, + check_ambiguity, override, 1, 0, "sProba", 's', + additional_error)) + goto failure; + + break; + case 'l': /* Sort against volume. */ + + + if (update_arg((void *)&(args_info->sVol_flag), 0, &(args_info->sVol_given), + &(local_args_info.sVol_given), optarg, 0, 0, ARG_FLAG, + check_ambiguity, override, 1, 0, "sVol", 'l', + additional_error)) + goto failure; + + break; + case 'm': /* Minimal probability to accept. */ + + + if (update_arg( (void *)&(args_info->mProba_arg), + &(args_info->mProba_orig), &(args_info->mProba_given), + &(local_args_info.mProba_given), optarg, 0, "0.0", ARG_DOUBLE, + check_ambiguity, override, 0, 0, + "mProba", 'm', + additional_error)) + goto failure; + + break; + case 'b': /* Ramses base output directory. */ + + + if (update_arg( (void *)&(args_info->ramsesDir_arg), + &(args_info->ramsesDir_orig), &(args_info->ramsesDir_given), + &(local_args_info.ramsesDir_given), optarg, 0, 0, ARG_STRING, + check_ambiguity, override, 0, 0, + "ramsesDir", 'b', + additional_error)) + goto failure; + + break; + case 'r': /* Ramses output id. */ + + + if (update_arg( (void *)&(args_info->ramsesId_arg), + &(args_info->ramsesId_orig), &(args_info->ramsesId_given), + &(local_args_info.ramsesId_given), optarg, 0, 0, ARG_INT, + check_ambiguity, override, 0, 0, + "ramsesId", 'r', + additional_error)) + goto failure; + + break; + case 'f': /* Configuration file. */ + + + if (update_arg( (void *)&(args_info->configFile_arg), + &(args_info->configFile_orig), &(args_info->configFile_given), + &(local_args_info.configFile_given), optarg, 0, 0, ARG_STRING, + check_ambiguity, override, 0, 0, + "configFile", 'f', + additional_error)) + goto failure; + + break; + case 'q': /* Quiet mode. */ + + + if (update_arg((void *)&(args_info->quiet_flag), 0, &(args_info->quiet_given), + &(local_args_info.quiet_given), optarg, 0, 0, ARG_FLAG, + check_ambiguity, override, 1, 0, "quiet", 'q', + additional_error)) + goto failure; + + break; + case 'p': /* Move the center by (x,y,z). */ + + + if (update_arg( (void *)&(args_info->move_arg), + &(args_info->move_orig), &(args_info->move_given), + &(local_args_info.move_given), optarg, 0, "(0,0,0)", ARG_STRING, + check_ambiguity, override, 0, 0, + "move", 'p', + additional_error)) + goto failure; + + break; + case 'o': /* Output file. */ + + + if (update_arg( (void *)&(args_info->output_arg), + &(args_info->output_orig), &(args_info->output_given), + &(local_args_info.output_given), optarg, 0, "voidDesc", ARG_STRING, + check_ambiguity, override, 0, 0, + "output", 'o', + additional_error)) + goto failure; + + break; + + case 0: /* Long option with no short option */ + /* Output is galax particles. */ + if (strcmp (long_options[option_index].name, "galax") == 0) + { + + + if (update_arg((void *)&(args_info->galax_flag), 0, &(args_info->galax_given), + &(local_args_info.galax_given), optarg, 0, 0, ARG_FLAG, + check_ambiguity, override, 1, 0, "galax", '-', + additional_error)) + goto failure; + + } + + break; + case '?': /* Invalid option. */ + /* `getopt_long' already printed an error message. */ + goto failure; + + default: /* bug: option not considered. */ + fprintf (stderr, "%s: option unknown: %c%s\n", ZOBOVCONF_PACKAGE, c, (additional_error ? additional_error : "")); + abort (); + } /* switch */ + } /* while */ + + + + if (check_required) + { + error += zobovConf_required2 (args_info, argv[0], additional_error); + } + + zobovConf_release (&local_args_info); + + if ( error ) + return (EXIT_FAILURE); + + return 0; + +failure: + + zobovConf_release (&local_args_info); + return (EXIT_FAILURE); +} + +#ifndef CONFIG_FILE_LINE_SIZE +#define CONFIG_FILE_LINE_SIZE 2048 +#endif +#define ADDITIONAL_ERROR " in configuration file " + +#define CONFIG_FILE_LINE_BUFFER_SIZE (CONFIG_FILE_LINE_SIZE+3) +/* 3 is for "--" and "=" */ + +static int +_zobovConf_configfile (char * const filename, int *my_argc) +{ + FILE* file; + char my_argv[CONFIG_FILE_LINE_BUFFER_SIZE+1]; + char linebuf[CONFIG_FILE_LINE_SIZE]; + int line_num = 0; + int result = 0, equal; + char *fopt, *farg; + char *str_index; + size_t len, next_token; + char delimiter; + + if ((file = fopen(filename, "r")) == NULL) + { + fprintf (stderr, "%s: Error opening configuration file '%s'\n", + ZOBOVCONF_PACKAGE, filename); + return EXIT_FAILURE; + } + + while ((fgets(linebuf, CONFIG_FILE_LINE_SIZE, file)) != NULL) + { + ++line_num; + my_argv[0] = '\0'; + len = strlen(linebuf); + if (len > (CONFIG_FILE_LINE_BUFFER_SIZE-1)) + { + fprintf (stderr, "%s:%s:%d: Line too long in configuration file\n", + ZOBOVCONF_PACKAGE, filename, line_num); + result = EXIT_FAILURE; + break; + } + + /* find first non-whitespace character in the line */ + next_token = strspn (linebuf, " \t\r\n"); + str_index = linebuf + next_token; + + if ( str_index[0] == '\0' || str_index[0] == '#') + continue; /* empty line or comment line is skipped */ + + fopt = str_index; + + /* truncate fopt at the end of the first non-valid character */ + next_token = strcspn (fopt, " \t\r\n="); + + if (fopt[next_token] == '\0') /* the line is over */ + { + farg = NULL; + equal = 0; + goto noarg; + } + + /* remember if equal sign is present */ + equal = (fopt[next_token] == '='); + fopt[next_token++] = '\0'; + + /* advance pointers to the next token after the end of fopt */ + next_token += strspn (fopt + next_token, " \t\r\n"); + + /* check for the presence of equal sign, and if so, skip it */ + if ( !equal ) + if ((equal = (fopt[next_token] == '='))) + { + next_token++; + next_token += strspn (fopt + next_token, " \t\r\n"); + } + str_index += next_token; + + /* find argument */ + farg = str_index; + if ( farg[0] == '\"' || farg[0] == '\'' ) + { /* quoted argument */ + str_index = strchr (++farg, str_index[0] ); /* skip opening quote */ + if (! str_index) + { + fprintf + (stderr, + "%s:%s:%d: unterminated string in configuration file\n", + ZOBOVCONF_PACKAGE, filename, line_num); + result = EXIT_FAILURE; + break; + } + } + else + { /* read up the remaining part up to a delimiter */ + next_token = strcspn (farg, " \t\r\n#\'\""); + str_index += next_token; + } + + /* truncate farg at the delimiter and store it for further check */ + delimiter = *str_index, *str_index++ = '\0'; + + /* everything but comment is illegal at the end of line */ + if (delimiter != '\0' && delimiter != '#') + { + str_index += strspn(str_index, " \t\r\n"); + if (*str_index != '\0' && *str_index != '#') + { + fprintf + (stderr, + "%s:%s:%d: malformed string in configuration file\n", + ZOBOVCONF_PACKAGE, filename, line_num); + result = EXIT_FAILURE; + break; + } + } + + noarg: + if (!strcmp(fopt,"include")) { + if (farg && *farg) { + result = _zobovConf_configfile(farg, my_argc); + } else { + fprintf(stderr, "%s:%s:%d: include requires a filename argument.\n", + ZOBOVCONF_PACKAGE, filename, line_num); + } + continue; + } + len = strlen(fopt); + strcat (my_argv, len > 1 ? "--" : "-"); + strcat (my_argv, fopt); + if (len > 1 && ((farg && *farg) || equal)) + strcat (my_argv, "="); + if (farg && *farg) + strcat (my_argv, farg); + ++(*my_argc); + + cmd_line_list_tmp = (struct line_list *) malloc (sizeof (struct line_list)); + cmd_line_list_tmp->next = cmd_line_list; + cmd_line_list = cmd_line_list_tmp; + cmd_line_list->string_arg = gengetopt_strdup(my_argv); + } /* while */ + + if (file) + fclose(file); + return result; +} + +int +zobovConf_configfile (char * const filename, + struct zobovConf_info *args_info, + int override, int initialize, int check_required) +{ + struct zobovConf_params params; + + params.override = override; + params.initialize = initialize; + params.check_required = check_required; + params.check_ambiguity = 0; + params.print_errors = 1; + + return zobovConf_config_file (filename, args_info, ¶ms); +} + +int +zobovConf_config_file (char * const filename, + struct zobovConf_info *args_info, + struct zobovConf_params *params) +{ + int i, result; + int my_argc = 1; + char **my_argv_arg; + char *additional_error; + + /* store the program name */ + cmd_line_list_tmp = (struct line_list *) malloc (sizeof (struct line_list)); + cmd_line_list_tmp->next = cmd_line_list; + cmd_line_list = cmd_line_list_tmp; + cmd_line_list->string_arg = gengetopt_strdup (ZOBOVCONF_PACKAGE); + + result = _zobovConf_configfile(filename, &my_argc); + + if (result != EXIT_FAILURE) { + my_argv_arg = (char **) malloc((my_argc+1) * sizeof(char *)); + cmd_line_list_tmp = cmd_line_list; + + for (i = my_argc - 1; i >= 0; --i) { + my_argv_arg[i] = cmd_line_list_tmp->string_arg; + cmd_line_list_tmp = cmd_line_list_tmp->next; + } + + my_argv_arg[my_argc] = 0; + + additional_error = (char *)malloc(strlen(filename) + strlen(ADDITIONAL_ERROR) + 1); + strcpy (additional_error, ADDITIONAL_ERROR); + strcat (additional_error, filename); + result = + zobovConf_internal (my_argc, my_argv_arg, args_info, + params, + additional_error); + + free (additional_error); + free (my_argv_arg); + } + + free_cmd_list(); + if (result == EXIT_FAILURE) + { + zobovConf_free (args_info); + exit (EXIT_FAILURE); + } + + return result; +} diff --git a/mytools/zobovConf.h b/mytools/zobovConf.h new file mode 100644 index 0000000..82245be --- /dev/null +++ b/mytools/zobovConf.h @@ -0,0 +1,243 @@ +/** @file zobovConf.h + * @brief The header file for the command line option parser + * generated by GNU Gengetopt version 2.22 + * http://www.gnu.org/software/gengetopt. + * DO NOT modify this file, since it can be overwritten + * @author GNU Gengetopt by Lorenzo Bettini */ + +#ifndef ZOBOVCONF_H +#define ZOBOVCONF_H + +/* If we use autoconf. */ +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include /* for FILE */ + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + +#ifndef ZOBOVCONF_PACKAGE +/** @brief the program name */ +#define ZOBOVCONF_PACKAGE "showZobov" +#endif + +#ifndef ZOBOVCONF_VERSION +/** @brief the program version */ +#define ZOBOVCONF_VERSION "0" +#endif + +/** @brief Where the command line options are stored */ +struct zobovConf_info +{ + const char *help_help; /**< @brief Print help and exit help description. */ + const char *version_help; /**< @brief Print version and exit help description. */ + char * desc_arg; /**< @brief The description file name for the voids. */ + char * desc_orig; /**< @brief The description file name for the voids original value given at command line. */ + const char *desc_help; /**< @brief The description file name for the voids help description. */ + char * adj_arg; /**< @brief Adjacent file name. */ + char * adj_orig; /**< @brief Adjacent file name original value given at command line. */ + const char *adj_help; /**< @brief Adjacent file name help description. */ + char * void_arg; /**< @brief Void/zone bind filename. */ + char * void_orig; /**< @brief Void/zone bind filename original value given at command line. */ + const char *void_help; /**< @brief Void/zone bind filename help description. */ + int id_arg; /**< @brief Void id to extract. */ + char * id_orig; /**< @brief Void id to extract original value given at command line. */ + const char *id_help; /**< @brief Void id to extract help description. */ + int sProba_flag; /**< @brief Sort against probability (default=off). */ + const char *sProba_help; /**< @brief Sort against probability help description. */ + int sVol_flag; /**< @brief Sort against volume (default=on). */ + const char *sVol_help; /**< @brief Sort against volume help description. */ + double mProba_arg; /**< @brief Minimal probability to accept (default='0.0'). */ + char * mProba_orig; /**< @brief Minimal probability to accept original value given at command line. */ + const char *mProba_help; /**< @brief Minimal probability to accept help description. */ + char * ramsesDir_arg; /**< @brief Ramses base output directory. */ + char * ramsesDir_orig; /**< @brief Ramses base output directory original value given at command line. */ + const char *ramsesDir_help; /**< @brief Ramses base output directory help description. */ + int ramsesId_arg; /**< @brief Ramses output id. */ + char * ramsesId_orig; /**< @brief Ramses output id original value given at command line. */ + const char *ramsesId_help; /**< @brief Ramses output id help description. */ + char * configFile_arg; /**< @brief Configuration file. */ + char * configFile_orig; /**< @brief Configuration file original value given at command line. */ + const char *configFile_help; /**< @brief Configuration file help description. */ + int quiet_flag; /**< @brief Quiet mode (default=off). */ + const char *quiet_help; /**< @brief Quiet mode help description. */ + char * move_arg; /**< @brief Move the center by (x,y,z) (default='(0,0,0)'). */ + char * move_orig; /**< @brief Move the center by (x,y,z) original value given at command line. */ + const char *move_help; /**< @brief Move the center by (x,y,z) help description. */ + int galax_flag; /**< @brief Output is galax particles (default=off). */ + const char *galax_help; /**< @brief Output is galax particles help description. */ + char * output_arg; /**< @brief Output file (default='voidDesc'). */ + char * output_orig; /**< @brief Output file original value given at command line. */ + const char *output_help; /**< @brief Output file help description. */ + + unsigned int help_given ; /**< @brief Whether help was given. */ + unsigned int version_given ; /**< @brief Whether version was given. */ + unsigned int desc_given ; /**< @brief Whether desc was given. */ + unsigned int adj_given ; /**< @brief Whether adj was given. */ + unsigned int void_given ; /**< @brief Whether void was given. */ + unsigned int id_given ; /**< @brief Whether id was given. */ + unsigned int sProba_given ; /**< @brief Whether sProba was given. */ + unsigned int sVol_given ; /**< @brief Whether sVol was given. */ + unsigned int mProba_given ; /**< @brief Whether mProba was given. */ + unsigned int ramsesDir_given ; /**< @brief Whether ramsesDir was given. */ + unsigned int ramsesId_given ; /**< @brief Whether ramsesId was given. */ + unsigned int configFile_given ; /**< @brief Whether configFile was given. */ + unsigned int quiet_given ; /**< @brief Whether quiet was given. */ + unsigned int move_given ; /**< @brief Whether move was given. */ + unsigned int galax_given ; /**< @brief Whether galax was given. */ + unsigned int output_given ; /**< @brief Whether output was given. */ + +} ; + +/** @brief The additional parameters to pass to parser functions */ +struct zobovConf_params +{ + int override; /**< @brief whether to override possibly already present options (default 0) */ + int initialize; /**< @brief whether to initialize the option structure zobovConf_info (default 1) */ + int check_required; /**< @brief whether to check that all required options were provided (default 1) */ + int check_ambiguity; /**< @brief whether to check for options already specified in the option structure zobovConf_info (default 0) */ + int print_errors; /**< @brief whether getopt_long should print an error message for a bad option (default 1) */ +} ; + +/** @brief the purpose string of the program */ +extern const char *zobovConf_info_purpose; +/** @brief the usage string of the program */ +extern const char *zobovConf_info_usage; +/** @brief all the lines making the help output */ +extern const char *zobovConf_info_help[]; + +/** + * The command line parser + * @param argc the number of command line options + * @param argv the command line options + * @param args_info the structure where option information will be stored + * @return 0 if everything went fine, NON 0 if an error took place + */ +int zobovConf (int argc, char * const *argv, + struct zobovConf_info *args_info); + +/** + * The command line parser (version with additional parameters - deprecated) + * @param argc the number of command line options + * @param argv the command line options + * @param args_info the structure where option information will be stored + * @param override whether to override possibly already present options + * @param initialize whether to initialize the option structure my_args_info + * @param check_required whether to check that all required options were provided + * @return 0 if everything went fine, NON 0 if an error took place + * @deprecated use zobovConf_ext() instead + */ +int zobovConf2 (int argc, char * const *argv, + struct zobovConf_info *args_info, + int override, int initialize, int check_required); + +/** + * The command line parser (version with additional parameters) + * @param argc the number of command line options + * @param argv the command line options + * @param args_info the structure where option information will be stored + * @param params additional parameters for the parser + * @return 0 if everything went fine, NON 0 if an error took place + */ +int zobovConf_ext (int argc, char * const *argv, + struct zobovConf_info *args_info, + struct zobovConf_params *params); + +/** + * Save the contents of the option struct into an already open FILE stream. + * @param outfile the stream where to dump options + * @param args_info the option struct to dump + * @return 0 if everything went fine, NON 0 if an error took place + */ +int zobovConf_dump(FILE *outfile, + struct zobovConf_info *args_info); + +/** + * Save the contents of the option struct into a (text) file. + * This file can be read by the config file parser (if generated by gengetopt) + * @param filename the file where to save + * @param args_info the option struct to save + * @return 0 if everything went fine, NON 0 if an error took place + */ +int zobovConf_file_save(const char *filename, + struct zobovConf_info *args_info); + +/** + * Print the help + */ +void zobovConf_print_help(void); +/** + * Print the version + */ +void zobovConf_print_version(void); + +/** + * Initializes all the fields a zobovConf_params structure + * to their default values + * @param params the structure to initialize + */ +void zobovConf_params_init(struct zobovConf_params *params); + +/** + * Allocates dynamically a zobovConf_params structure and initializes + * all its fields to their default values + * @return the created and initialized zobovConf_params structure + */ +struct zobovConf_params *zobovConf_params_create(void); + +/** + * Initializes the passed zobovConf_info structure's fields + * (also set default values for options that have a default) + * @param args_info the structure to initialize + */ +void zobovConf_init (struct zobovConf_info *args_info); +/** + * Deallocates the string fields of the zobovConf_info structure + * (but does not deallocate the structure itself) + * @param args_info the structure to deallocate + */ +void zobovConf_free (struct zobovConf_info *args_info); + +/** + * The config file parser (deprecated version) + * @param filename the name of the config file + * @param args_info the structure where option information will be stored + * @param override whether to override possibly already present options + * @param initialize whether to initialize the option structure my_args_info + * @param check_required whether to check that all required options were provided + * @return 0 if everything went fine, NON 0 if an error took place + * @deprecated use zobovConf_config_file() instead + */ +int zobovConf_configfile (char * const filename, + struct zobovConf_info *args_info, + int override, int initialize, int check_required); + +/** + * The config file parser + * @param filename the name of the config file + * @param args_info the structure where option information will be stored + * @param params additional parameters for the parser + * @return 0 if everything went fine, NON 0 if an error took place + */ +int zobovConf_config_file (char * const filename, + struct zobovConf_info *args_info, + struct zobovConf_params *params); + +/** + * Checks that all the required options were specified + * @param args_info the structure to check + * @param prog_name the name of the program that will be used to print + * possible errors + * @return + */ +int zobovConf_required (struct zobovConf_info *args_info, + const char *prog_name); + + +#ifdef __cplusplus +} +#endif /* __cplusplus */ +#endif /* ZOBOVCONF_H */ diff --git a/readfiles.c b/readfiles.c new file mode 100644 index 0000000..d6ac62f --- /dev/null +++ b/readfiles.c @@ -0,0 +1,121 @@ +#include +#include + +#define DL for (d=0;d<3;d++) /* Dimension loop */ +#define BF 1e30 + +/* Positions */ +/* Returns number of particles read */ +int posread(char *posfile, float ***p, float fact) { + + FILE *pos; + int np,dum,d,i; + float xmin,xmax,ymin,ymax,zmin,zmax; + float *ptemp; + + pos = fopen(posfile, "r"); + if (pos == NULL) { + printf("Unable to open position file %s\n\n",posfile); + exit(0); + } + /* Fortran77 4-byte headers and footers */ + /* Delete "dum" statements if you don't need them */ + + /* Read number of particles */ + fread(&dum,1,4,pos); fread(&np,1, sizeof(int),pos); fread(&dum,1,4,pos); + + /* Allocate the arrays */ + (*p) = (float **)malloc(np*sizeof(float *)); + ptemp = (float *)malloc(np*sizeof(float)); + + printf("np = %d\n",np); + + /* Fill the arrays */ + fread(&dum,1,4,pos); + fread(ptemp,np,4,pos); + for (i=0; ixmax) xmax = (*p)[i][0]; + if ((*p)[i][1]ymax) ymax = (*p)[i][1]; + if ((*p)[i][2]zmax) zmax = (*p)[i][2]; + } + printf("np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np,xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout); + + return(np); +} + +/* Velocities */ +/* Returns number of particles read */ +/* Used in voboz, but not zobov, which doesn't use velocities */ +int velread(char *velfile, float ***v, float fact) { + + FILE *vel; + int np,dum,d,i; + float xmin,xmax,ymin,ymax,zmin,zmax; + + vel = fopen(velfile, "r"); + if (vel == NULL) { + printf("Unable to open velocity file %s\n\n",velfile); + exit(0); + } + /* Fortran77 4-byte headers and footers */ + /* Delete "dum" statements if you don't need them */ + + /* Read number of particles */ + fread(&dum,1,4,vel); fread(&np,1, sizeof(int),vel); fread(&dum,1,4,vel); + + /* Allocate the arrays */ + (*v) = (float **)malloc(3*sizeof(float*)); + for (i=0;i<3;i++) (*v)[i] = (float *)malloc(np*sizeof(float)); + + /* Fill the arrays */ + fread(&dum,1,4,vel); fread((*v)[0],np,4,vel); fread(&dum,1,4,vel); + fread(&dum,1,4,vel); fread((*v)[1],np,4,vel); fread(&dum,1,4,vel); + fread(&dum,1,4,vel); fread((*v)[2],np,4,vel); fread(&dum,1,4,vel); + + fclose(vel); + + /* Convert from code units into physical units (km/sec) */ + + for (i=0; i xmax) xmax = (*v)[0][i]; + if ((*v)[1][i] < ymin) ymin = (*v)[1][i]; if ((*v)[1][i] > ymax) ymax = (*v)[1][i]; + if ((*v)[2][i] < zmin) zmin = (*v)[2][i]; if ((*v)[2][i] > zmax) zmax = (*v)[2][i]; + } + printf("vx: %f,%f; vy: %f,%f; vz: %f,%f\n",xmin,xmax, ymin,ymax, zmin,zmax); + + return(np); +} diff --git a/readme.txt b/readme.txt new file mode 100644 index 0000000..cb91c8a --- /dev/null +++ b/readme.txt @@ -0,0 +1,25 @@ +This is the readme file for ZOBOV (ZOnes Bordering On Voidness), +version 1.0 (Jan 5, 2008) a void-finding algorithm for sets of +particles (e.g. a cosmological N-body simulation) by Mark Neyrinck. +It is an inversion of the halo-finding algorithm VOBOZ (VOronoi BOund +Zones), which was also written by Mark Neyrinck, with the guidance of +Nick Gnedin and Andrew Hamilton. + +The gz file which once contained this file should include the +following files: + +Makefile jozov.c readme.txt voz1b1.c voztie.c +findrtop.c readfiles.c voz.h vozinit.c vozutil.c + +For instructions on installing and running ZOBOV, please see +http://www.ifa.hawaii.edu/~neyrinck/zobov/. To use ZOBOV, you will +need to download the Qhull package, at http://www.qhull.org. For help +in understanding ZOBOV, please see the paper describing it +(arXiv:0712.3049), If these resources are insufficient, feel free to +email Mark.Neyrinck@colorado.edu with questions. + +This is free software. It may be freely copied, modified, and +redistributed, as long as ZOBOV and its author are acknowledged. +There is no warranty or other guarantee of fitness for ZOBOV; it is +provided "as is." I welcome bug reports (but do not guarantee their +fixing), which should be sent to Mark.Neyrinck@colorado.edu. diff --git a/voz.h b/voz.h new file mode 100644 index 0000000..3a59560 --- /dev/null +++ b/voz.h @@ -0,0 +1,8 @@ +#define MAXVERVER 2000 +#define NGUARD 42 /*Actually, the number of SPACES between guard points + in each dim */ + +typedef struct Partadj { + int nadj; + int *adj; +} PARTADJ; diff --git a/voz1b1.c b/voz1b1.c new file mode 100644 index 0000000..1d97121 --- /dev/null +++ b/voz1b1.c @@ -0,0 +1,334 @@ +#include "qhull_a.h" +#include "voz.h" + +#define DL for (d=0;d<3;d++) +#define BF 1e30 + +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); + +int main(int argc, char *argv[]) { + int exitcode; + int i, j, np; + float **r; + coordT rtemp[3], *parts; + coordT deladjs[3*MAXVERVER], points[3*MAXVERVER]; + pointT intpoints[3*MAXVERVER]; + FILE *pos, *out; + char *posfile, outfile[80], *suffix; + PARTADJ *adjs; + float *vols; + float predict, xmin,xmax,ymin,ymax,zmin,zmax; + int *orig; + + int isitinbuf; + char isitinmain, d; + int numdiv, nvp, nvpall, nvpbuf; + float width, width2, totwidth, totwidth2, bf, s, g; + float border, boxsize; + float c[3]; + int b[3]; + double totalvol; + + if (argc != 9) { + printf("Wrong number of arguments.\n"); + printf("arg1: position file\n"); + printf("arg2: border size\n"); + printf("arg3: box size\n"); + printf("arg4: suffix\n"); + printf("arg5: number of divisions\n"); + printf("arg6-8: b[0-2]\n\n"); + exit(0); + } + posfile = argv[1]; + if (sscanf(argv[2],"%f",&border) != 1) { + printf("That's no border size; try again.\n"); + exit(0); + } + if (sscanf(argv[3],"%f",&boxsize) != 1) { + printf("That's no boxsize; try again.\n"); + exit(0); + } + suffix = argv[4]; + if (sscanf(argv[5],"%d",&numdiv) != 1) { + printf("%s is no number of divisions; try again.\n",argv[5]); + exit(0); + } + if (numdiv == 1) { + printf("Only using one division; should only use for an isolated segment.\n"); + } + if (numdiv < 1) { + printf("Cannot have a number of divisions less than 1. Resetting to 1.\n"); + numdiv = 1; + } + if (sscanf(argv[6],"%d",&b[0]) != 1) { + printf("That's no b index; try again.\n"); + exit(0); + } + if (sscanf(argv[7],"%d",&b[1]) != 1) { + printf("That's no b index; try again.\n"); + exit(0); + } + if (sscanf(argv[8],"%d",&b[2]) != 1) { + printf("That's no b index; try again.\n"); + exit(0); + } + + /* Boxsize should be the range in r, yielding a range 0-1 */ + np = posread(posfile,&r,1./boxsize); + printf("%d particles\n",np);fflush(stdout); + xmin = BF; xmax = -BF; ymin = BF; ymax = -BF; zmin = BF; zmax = -BF; + for (i=0; ixmax) xmax = r[i][0]; + if (r[i][1]ymax) ymax = r[i][1]; + if (r[i][2]zmax) zmax = r[i][2]; + } + printf("np: %d, x: %f,%f; y: %f,%f; z: %f,%f\n",np,xmin,xmax, ymin,ymax, zmin,zmax); fflush(stdout); + + 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]+0.5)*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 = (int *)malloc(nvpbuf*sizeof(int)); + + 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]) 0) && + ((fabs(rtemp[0])>width2)||(fabs(rtemp[1])>width2)||(fabs(rtemp[2])>width2))) { + + /*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);*/ + parts[3*nvpbuf] = rtemp[0]; + parts[3*nvpbuf+1] = rtemp[1]; + parts[3*nvpbuf+2] = rtemp[2]; + orig[nvpbuf] = i; + + nvpbuf++; + 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("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); + + 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); + + /* 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(int),out); + else printf("0"); + } + fclose(out); + + return(0); +} diff --git a/vozinit.c b/vozinit.c new file mode 100644 index 0000000..c26d860 --- /dev/null +++ b/vozinit.c @@ -0,0 +1,153 @@ +#include "qhull_a.h" +#include "voz.h" + +#define DL for (d=0;d<3;d++) +#define BF 1e30 + +int posread(char *posfile, float ***p, float fact); + +int main(int argc, char *argv[]) { + int i, np; + float **rfloat, rtemp[3]; + FILE *pos, *scr; + char *posfile, scrfile[80], systemstr[90], *suffix; + float xmin,xmax,ymin,ymax,zmin,zmax; + + int isitinbuf; + char isitinmain, d; + int numdiv; + int nvp, nvpall, nvpbuf, nvpmin, nvpmax, nvpbufmin, nvpbufmax; /* yes, the insurance */ + float width, width2, totwidth, totwidth2, bf, s, g; + float border, boxsize; + float c[3]; + int b[3]; + + if (argc != 6) { + printf("Wrong number of arguments.\n"); + printf("arg1: position file\n"); + printf("arg2: buffer size (default 0.1)\n"); + printf("arg3: box size\n"); + printf("arg4: number of divisions (default 2)\n"); + printf("arg5: suffix describing this run\n\n"); + exit(0); + } + posfile = argv[1]; + suffix = argv[2]; + if (sscanf(suffix,"%f",&border) != 1) { + printf("That's no border size; try again.\n"); + exit(0); + } + suffix = argv[3]; + if (sscanf(suffix,"%f",&boxsize) != 1) { + printf("That's no boxsize; try again.\n"); + exit(0); + } + suffix = argv[4]; + if (sscanf(suffix,"%d",&numdiv) != 1) { + printf("That's no number of divisions; try again.\n"); + exit(0); + } + if (numdiv < 2) { + printf("Cannot have a number of divisions less than 2. Resetting to 2:\n"); + numdiv = 2; + } + + suffix = argv[5]; + + /* Read the position file */ + np = posread(posfile,&rfloat,1./boxsize); + /* Boxsize should be the range in r, yielding a range 0-1 */ + + 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("Not enough guard points for given border.\nIncrease guards to >= %f\n.", + totwidth/sqrt(0.5*bf*bf)); + printf("bf = %f\n",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); + + nvpmax = 0; nvpbufmax = 0; nvpmin = np; nvpbufmin = np; + + for (b[0] = 0; b[0] < numdiv; b[0]++) { + c[0] = ((float)b[0]+0.5)*width; + for (b[1] = 0; b[1] < numdiv; b[1]++) { + c[1] = ((float)b[1]+0.5)*width; + for (b[2] = 0; b[2] < numdiv; b[2]++) { + c[2] = ((float)b[2]+0.5)*width; + + nvp = 0; /* Number of particles excluding buffer */ + nvpbuf = 0; /* Number of particles to tesselate, including + 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] ++; + isitinbuf = isitinbuf && (fabs(rtemp[d]) < totwidth2); + isitinmain = isitinmain && (fabs(rtemp[d]) <= width2); + } + if (isitinbuf) { + nvpbuf++; + } + if (isitinmain) { + 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]; + } + } + if (nvp > nvpmax) nvpmax = nvp; + if (nvpbuf > nvpbufmax) nvpbufmax = nvpbuf; + if (nvp < nvpmin) nvpmin = nvp; + if (nvpbuf < nvpbufmin) nvpbufmin = nvpbuf; + + printf("b=(%d,%d,%d), c=(%f,%f,%f), nvp=%d, nvpbuf=%d\n", + b[0],b[1],b[2],c[0],c[1],c[2],nvp,nvpbuf); + } + } + } + printf("Nvp range: %d,%d\n",nvpmin,nvpmax); + printf("Nvpbuf range: %d,%d\n",nvpbufmin,nvpbufmax); + + /* Output script file */ + sprintf(scrfile,"scr%s",suffix); + printf("Writing script file to %s.\n",scrfile);fflush(stdout); + scr = fopen(scrfile,"w"); + if (scr == NULL) { + printf("Problem opening script file.\n"); + fflush(stdout); + exit(0); + } + fprintf(scr,"#!/bin/bash -f\n"); + for (b[0]=0;b[0] +#include +#include +#include "voz.h" + +int main(int argc, char *argv[]) { + + FILE *part, *adj, *vol; + char partfile[80], *suffix, adjfile[80], volfile[80]; + float *vols, volstemp; + + PARTADJ *adjs; + + int numdiv,np,np2,na; + + int i,j,k,p,nout; + int nvp,npnotdone,nvpmax, nvpsum, *orig; + double avgnadj, avgvol; + + if (argc != 3) { + printf("Wrong number of arguments.\n"); + printf("arg1: number of divisions (default 2)\n"); + printf("arg2: suffix describing this run\n\n"); + exit(0); + } + if (sscanf(argv[1],"%d",&numdiv) != 1) { + printf("That's no number of divisions; try again.\n"); + exit(0); + } + if (numdiv < 2) { + printf("Cannot have a number of divisions less than 2. Resetting to 2:\n"); + numdiv = 2; + } + suffix = argv[2]; + + np = -1; nvpmax = -1; nvpsum = 0; + + for (i = 0; i < numdiv; i++) { + for (j = 0; j < numdiv; j++) { + for (k = 0; k < numdiv; k++) { + sprintf(partfile,"part.%s.%02d.%02d.%02d",suffix,i,j,k); + part = fopen(partfile,"r"); + if (part == NULL) { + printf("Unable to open file %s.\n\n",partfile); + exit(0); + } + fread(&np2,1,sizeof(int),part); + fread(&nvp,1,sizeof(int),part); + if (np == -1) + np = np2; + else + if (np2 != np) { + printf("Incompatible total particle numbers: %d,%d\n\n",np,np2); + exit(0); + } + if (nvp > nvpmax) nvpmax = nvp; + fclose(part); + } + } + } + printf("We have %d particles to tie together.\n",np); fflush(stdout); + printf("The maximum number of particles in a file is %d.\n",nvpmax); + + adjs = (PARTADJ *)malloc(np*sizeof(PARTADJ)); + 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)); + if (orig == NULL) printf("Couldn't allocate orig.\n"); + if ((vols == NULL) || (orig == NULL) || (adjs == NULL)) { + printf("Not enough memory to allocate. Exiting.\n"); + exit(0); + } + for (p=0;p -1.) + if (fabs(vols[orig[p]]-volstemp)/volstemp > 1e-3) { + printf("Inconsistent volumes for p. %d: (%g,%g)!\n", + orig[p],vols[orig[p]],volstemp); + exit(0); + } + vols[orig[p]] = volstemp; + } + + for (p=0;p 0) { + adjs[orig[p]].nadj = na; + adjs[orig[p]].adj = (int *)malloc(na*sizeof(int)); + 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); + } else { + printf("0"); fflush(stdout); + } + } + fclose(part); + printf("%d ",k); + } + } + } + printf("\n"); + npnotdone = 0; avgnadj = 0.; avgvol = 0.; + for (p=0;p 0) + printf("%d particles not done!\n"); + printf("%d particles done more than once.\n",nvpsum-np); + avgnadj /= (double)np; + avgvol /= (double)np; + printf("Average # adjacencies = %lf (%f for Poisson)\n",avgnadj, + 48.*3.141593*3.141593/35.+2.); + printf("Average volume = %lf\n",avgvol); + + /* Now the output! */ + + sprintf(adjfile,"adj%s.dat",suffix); + sprintf(volfile,"vol%s.dat",suffix); + + printf("Outputting to %s, %s\n\n",adjfile,volfile); + + adj = fopen(adjfile,"w"); + if (adj == NULL) { + printf("Unable to open %s\n",adjfile); + exit(0); + } + fwrite(&np,1, sizeof(int),adj); + /* Adjacencies: first the numbers of adjacencies, + and the number we're actually going to write per particle */ + for (i=0;i 0) { + nout = 0; + for (j=0;j i) nout++; + fwrite(&nout,1,sizeof(int),adj); + for (j=0;j i) + fwrite(&(adjs[i].adj[j]),1,sizeof(int),adj); + } + + fclose(adj); + + /* Volumes */ + vol = fopen(volfile,"w"); + if (vol == NULL) { + printf("Unable to open %s\n",volfile); + exit(0); + } + fwrite(&np,1, sizeof(int),vol); + fwrite(vols,sizeof(float),np,vol); + + fclose(vol); + + return(0); +} diff --git a/vozutil.c b/vozutil.c new file mode 100644 index 0000000..1e2aee0 --- /dev/null +++ b/vozutil.c @@ -0,0 +1,241 @@ +#include "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; + + 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"); + 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, NULL, &numfacets, &numsimplicial, + &totneighbors, &numridges, &numcoplanars, &numtricoplanars); + qh_vertexneighbors(); + vertices= qh_facetvertices (qh facet_list, NULL, NULL); + 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) { + 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; +}