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

This commit is contained in:
P.M. Sutter 2013-03-28 21:19:07 -05:00
commit 34400fbaca
3 changed files with 75 additions and 92 deletions

View file

@ -101,6 +101,8 @@ int main(int argc,char **argv) {
exit(0); exit(0);
} }
fread(&np,1, sizeof(int),adj); fread(&np,1, sizeof(int),adj);
if (mockIndex < 0)
mockIndex = np;
printf("adj: %d particles\n", np); printf("adj: %d particles\n", np);
FF; FF;
@ -125,14 +127,21 @@ int main(int argc,char **argv) {
assert(i < j); assert(i < j);
if (p[i].ncnt == p[i].nadj) if (p[i].ncnt == p[i].nadj)
{ {
p[i].adj = (pid_t *)realloc(p[i].adj, (p[i].nadj+1)*sizeof(pid_t)); int q;
p[i].nadj++; printf("OVERFLOW for particle %d (pending %d). List of accepted:\n", i, j);
for (q=0;q<p[i].nadj;q++)
printf(" %d\n", p[i].adj[q]);
abort();
} }
if (p[j].ncnt == p[j].nadj) if (p[j].ncnt == p[j].nadj)
{ {
p[j].adj = (pid_t *)realloc(p[j].adj, (p[j].nadj+1)*sizeof(pid_t)); int q;
p[j].nadj++; printf("OVERFLOW for particle %d (pending %d). List of accepted:\n", j, i);
for (q=0;q<p[j].nadj;q++)
printf(" %d\n", p[j].adj[q]);
abort();
} }
p[i].adj[p[i].ncnt] = j; p[i].adj[p[i].ncnt] = j;
p[j].adj[p[j].ncnt] = i; p[j].adj[p[j].ncnt] = i;
p[i].ncnt++; p[j].ncnt++; p[i].ncnt++; p[j].ncnt++;

View file

@ -21,16 +21,16 @@ int main(int argc, char *argv[]) {
char *posfile, outfile[200], *suffix, *outDir; char *posfile, outfile[200], *suffix, *outDir;
PARTADJ *adjs; PARTADJ *adjs;
float *vols; float *vols;
float predict, xmin,xmax,ymin,ymax,zmin,zmax; realT predict, xmin,xmax,ymin,ymax,zmin,zmax;
pid_t *orig; pid_t *orig;
int isitinbuf; int isitinbuf;
char isitinmain, d; char isitinmain, d;
int numdiv; int numdiv;
pid_t nvp, nvpall, nvpbuf; pid_t nvp, nvpall, nvpbuf;
float width, width2, totwidth, totwidth2, bf, s, g; realT width, width2, totwidth, totwidth2, bf, s, g;
float border, boxsize; float border, boxsize;
float c[3]; realT c[3];
int b[3]; int b[3];
double totalvol; double totalvol;
@ -118,7 +118,7 @@ int main(int argc, char *argv[]) {
exit(0); exit(0);
} }
DL c[d] = ((float)b[d]+0.5)*width; DL c[d] = ((float)b[d])*width;
printf("c: %f,%f,%f\n",c[0],c[1],c[2]); printf("c: %f,%f,%f\n",c[0],c[1],c[2]);
/* Assign temporary array*/ /* Assign temporary array*/
nvpbuf = 0; /* Number of particles to tesselate, including nvpbuf = 0; /* Number of particles to tesselate, including
@ -159,7 +159,7 @@ int main(int argc, char *argv[]) {
for (i=0; i<np; i++) { for (i=0; i<np; i++) {
isitinmain = 1; isitinmain = 1;
DL { DL {
rtemp[d] = r[i][d] - c[d]; rtemp[d] = (realT)r[i][d] - (realT)c[d];
if (rtemp[d] > 0.5) rtemp[d] --; if (rtemp[d] > 0.5) rtemp[d] --;
if (rtemp[d] < -0.5) rtemp[d] ++; if (rtemp[d] < -0.5) rtemp[d] ++;
isitinmain = isitinmain && (fabs(rtemp[d]) <= width2); isitinmain = isitinmain && (fabs(rtemp[d]) <= width2);
@ -186,7 +186,7 @@ int main(int argc, char *argv[]) {
isitinmain = 1; isitinmain = 1;
DL { DL {
rtemp[d] = r[i][d] - c[d]; rtemp[d] = (realT)r[i][d] - (realT)c[d];
if (rtemp[d] > 0.5) rtemp[d] --; if (rtemp[d] > 0.5) rtemp[d] --;
if (rtemp[d] < -0.5) rtemp[d] ++; if (rtemp[d] < -0.5) rtemp[d] ++;
isitinbuf = isitinbuf && (fabs(rtemp[d])<totwidth2); isitinbuf = isitinbuf && (fabs(rtemp[d])<totwidth2);
@ -223,40 +223,40 @@ int main(int argc, char *argv[]) {
for (i=0; i<NGUARD+1; i++) { for (i=0; i<NGUARD+1; i++) {
for (j=0; j<NGUARD+1; j++) { for (j=0; j<NGUARD+1; j++) {
/* Bottom */ /* Bottom */
parts[3*nvpall] = -width2 + (float)i * s; parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 + (float)j * s; parts[3*nvpall+1] = -width2 + (realT)j * s;
parts[3*nvpall+2] = -width2 - g; parts[3*nvpall+2] = -width2 - g;
nvpall++; nvpall++;
/* Top */ /* Top */
parts[3*nvpall] = -width2 + (float)i * s; parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 + (float)j * s; parts[3*nvpall+1] = -width2 + (realT)j * s;
parts[3*nvpall+2] = width2 + g; parts[3*nvpall+2] = width2 + g;
nvpall++; nvpall++;
} }
} }
for (i=0; i<NGUARD+1; i++) { /* Don't want to overdo the corners*/ for (i=0; i<NGUARD+1; i++) { /* Don't want to overdo the corners*/
for (j=0; j<NGUARD+1; j++) { for (j=0; j<NGUARD+1; j++) {
parts[3*nvpall] = -width2 + (float)i * s; parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = -width2 - g; parts[3*nvpall+1] = -width2 - g;
parts[3*nvpall+2] = -width2 + (float)j * s; parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++; nvpall++;
parts[3*nvpall] = -width2 + (float)i * s; parts[3*nvpall] = -width2 + (realT)i * s;
parts[3*nvpall+1] = width2 + g; parts[3*nvpall+1] = width2 + g;
parts[3*nvpall+2] = -width2 + (float)j * s; parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++; nvpall++;
} }
} }
for (i=0; i<NGUARD+1; i++) { for (i=0; i<NGUARD+1; i++) {
for (j=0; j<NGUARD+1; j++) { for (j=0; j<NGUARD+1; j++) {
parts[3*nvpall] = -width2 - g; parts[3*nvpall] = -width2 - g;
parts[3*nvpall+1] = -width2 + (float)i * s; parts[3*nvpall+1] = -width2 + (realT)i * s;
parts[3*nvpall+2] = -width2 + (float)j * s; parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++; nvpall++;
parts[3*nvpall] = width2 + g; parts[3*nvpall] = width2 + g;
parts[3*nvpall+1] = -width2 + (float)i * s; parts[3*nvpall+1] = -width2 + (realT)i * s;
parts[3*nvpall+2] = -width2 + (float)j * s; parts[3*nvpall+2] = -width2 + (realT)j * s;
nvpall++; nvpall++;
} }
} }

View file

@ -1,3 +1,4 @@
#include <assert.h>
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
@ -96,7 +97,11 @@ int main(int argc, char *argv[]) {
exit(0); exit(0);
} }
for (p=0;p<np;p++) for (p=0;p<np;p++)
{
vols[p] = -1.; vols[p] = -1.;
adjs[p].nadj = 0;
adjs[p].adj = 0;
}
nvpsum = 0; nvpsum = 0;
for (i = 0; i < numdiv; i++) { for (i = 0; i < numdiv; i++) {
@ -118,13 +123,9 @@ int main(int argc, char *argv[]) {
for (p=0;p<nvp;p++) { for (p=0;p<nvp;p++) {
fread(&volstemp,1,sizeof(float),part); fread(&volstemp,1,sizeof(float),part);
if (vols[orig[p]] > -1.) if (vols[orig[p]] > -1.)
// PMS
if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3 && orig[p] < mockIndex) { if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3 && orig[p] < mockIndex) {
//if (fabs(vols[orig[p]]-volstemp)/volstemp > 1.5e-3) {
// END PMS
printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n", printf("Inconsistent volumes for p. %d: (%10g,%10g)!\n",
orig[p],vols[orig[p]],volstemp); orig[p],vols[orig[p]],volstemp);
// TEST
//exit(0); //exit(0);
} }
vols[orig[p]] = volstemp; vols[orig[p]] = volstemp;
@ -132,10 +133,13 @@ int main(int argc, char *argv[]) {
for (p=0;p<nvp;p++) { for (p=0;p<nvp;p++) {
fread(&na,1,sizeof(int),part); fread(&na,1,sizeof(int),part);
pid_t pid = orig[p];
if (na > 0) { if (na > 0) {
assert(adjs[pid].nadj == 0);// || adjs[pid].nadj == na);
adjs[orig[p]].nadj = na; adjs[orig[p]].nadj = na;
if (adjs[pid].adj == 0)
adjs[orig[p]].adj = (pid_t *)malloc(na*sizeof(pid_t)); adjs[orig[p]].adj = (pid_t *)malloc(na*sizeof(pid_t));
if (adjs[orig[p]].adj == NULL) { if (adjs[orig[p]].adj == 0) {
printf("Couldn't allocate adjs[orig[%d]].adj.\n",p); printf("Couldn't allocate adjs[orig[%d]].adj.\n",p);
exit(0); exit(0);
} }
@ -179,7 +183,6 @@ int main(int argc, char *argv[]) {
int numAdjSaved = 0; int numAdjSaved = 0;
for (j = 0; j < adjs[i].nadj; j++) { for (j = 0; j < adjs[i].nadj; j++) {
//if ( vols[adjs[i].adj[j]] != -99) {
if ( adjs[adjs[i].adj[j]].nadj != 0) { if ( adjs[adjs[i].adj[j]].nadj != 0) {
adjs[i].adj[numAdjSaved] = adjs[i].adj[j]; adjs[i].adj[numAdjSaved] = adjs[i].adj[j];
numAdjSaved++; numAdjSaved++;
@ -187,21 +190,8 @@ int main(int argc, char *argv[]) {
} }
adjs[i].nadj = numAdjSaved; adjs[i].nadj = numAdjSaved;
} }
/*
for (i = 0; i < mockIndex; i++) {
printf("ADJ: %d %d : ", i, adjs[i].nadj);
for (j = 0; j < adjs[i].nadj; j++) {
printf(" %d", adjs[i].adj[j]);
}
printf("\n");
}
*/
printf("Removed %d mock galaxies and %d adjacent galaxies.\n", np-mockIndex, printf("Removed %d mock galaxies and %d adjacent galaxies.\n", np-mockIndex,
numRemoved); numRemoved);
printf("There are %d galaxies remaining.\n", mockIndex-numRemoved); printf("There are %d galaxies remaining.\n", mockIndex-numRemoved);
@ -238,39 +228,23 @@ printf("\n");
printf("Unable to open %s\n",adjfile); printf("Unable to open %s\n",adjfile);
exit(0); exit(0);
} }
// PMS
fwrite(&mockIndex,1, sizeof(int),adj); fwrite(&mockIndex,1, sizeof(int),adj);
//fwrite(&np,1, sizeof(int),adj);
// END OMS
/* Adjacencies: first the numbers of adjacencies, /* Adjacencies: first the numbers of adjacencies,
and the number we're actually going to write per particle */ and the number we're actually going to write per particle */
// Recount the number of adjacencies after merge // Recount the number of adjacencies after merge
for(i=0;i<mockIndex;i++) for(i=0;i<mockIndex;i++)
cnt_adj[i] = 0; cnt_adj[i] = adjs[i].nadj;
for(i=0;i<mockIndex;i++)
{
for (j=0;j<adjs[i].nadj;j++)
{
cnt_adj[adjs[i].adj[j]]++;
cnt_adj[i]++;
}
}
// PMS
for (i=0;i<mockIndex;i++) for (i=0;i<mockIndex;i++)
//for (i=0;i<np;i++)
// END PMS
fwrite(&cnt_adj[i],1,sizeof(int),adj); fwrite(&cnt_adj[i],1,sizeof(int),adj);
/* Now the lists of adjacencies (without double counting) */ /* Now the lists of adjacencies (without double counting) */
// PMS
for (i=0;i<mockIndex;i++) { for (i=0;i<mockIndex;i++) {
//for (i=0;i<np;i++) {
//if (adjs[i].nadj > 0) {
// END PMS
nout = 0; nout = 0;
for (j=0;j<adjs[i].nadj; j++) if (adjs[i].adj[j] > i) nout++; for (j=0;j<adjs[i].nadj; j++)
if (adjs[i].adj[j] > i)
nout++;
fwrite(&nout,1,sizeof(int),adj); fwrite(&nout,1,sizeof(int),adj);
for (j=0;j<adjs[i].nadj; j++) { for (j=0;j<adjs[i].nadj; j++) {
pid_t id = adjs[i].adj[j]; pid_t id = adjs[i].adj[j];