Add bond generation code in C.

job.sh

@@ -23,15 +23,18 @@ git checkout $commit
 NUM_REPLICAS=$num_temps ./setup.sh
 make
 
-# Generate samples
-seeds=$(python bonds.py $((2**LOG_N)) -z $Z --sigma $sigma --ns $num_cores)
+if [ "$GENERATE_SAMPLE" -eq 0 ]; then
+    seeds=$(python bonds.py $((2**LOG_N)) -z $Z --sigma $sigma --ns $num_cores)
+else
+    seeds=$(python seeds.py $num_cores)
+fi
 
 # Run
 # Don't read or write to /home from here
 time mpirun ./isg $seeds
 
 # Move output to home directory, clean up
-[ "$KEEP_SAMPLES" -eq 1 ] || rm samp_*.txt
+[ "$KEEP_SAMPLES" -eq 1 ] || rm -f samp_*.txt
 bzip2 -9 *.txt
 mv *.bz2 $output_dir
 mv *.h $output_dir

seeds.py

@@ -0,0 +1,11 @@
+from numpy.random import random_integers
+
+if __name__=='__main__':
+    import sys
+    num_seeds = int(sys.argv[1]) if len(sys.argv) > 1 else 1
+    seeds = random_integers(0, 2**24, num_seeds)
+
+    for seed in seeds:
+        print '{:06x}'.format(seed)
+
+

setup.sh

@@ -17,6 +17,8 @@ function setp {
 getp srcdir $srcdir
 
 getp LOG_N 7
+getp GENERATE_SAMPLE 0
+getp SIGMA "0.6"
 getp DILUTE 1
 getp Z 6
 getp Z_MAX 24

src/Makefile.in

@@ -14,7 +14,9 @@ LDLIBS = -lm -lgsl -lgslcblas
 MODULES = mod_meas.o mod_hist.o mod_corr.o
 endif
 
-isg: output.o state.o exchange.o replica.o sample.o $(MODULES)
+isg: bonds.o bond_dist.o output.o state.o exchange.o replica.o sample.o $(MODULES)
+
+genbonds: bond_dist.o bonds.o
 
 clean:
 	$(RM) *.o

src/bisect.c

@@ -0,0 +1,13 @@
+static int bisect(double a[], double x, int l, int r)
+{
+    int mid;
+
+    while (r-l > 1)
+    {
+        mid = (l+r)/2;
+        if (a[mid] < x) l = mid;
+        else r = mid;
+    }
+
+    return (x < a[l]) ? l : r;
+}

src/bond_dist.c

@@ -0,0 +1,60 @@
+#include "bond_dist.h"
+#include <math.h>
+
+#include "bisect.c"
+
+#define R(n) fabs(N/M_PI*sin(M_PI*n/N))
+
+void bond_dist_dilute(bonds *b, double sigma, rng *rand)
+{
+    int n, num_bonds = 0;
+    double *f, sum = 0., p = -2.*sigma;
+
+    f = malloc(N * sizeof(double));
+
+    /* compute probability distribution */
+
+    f[0] = 0.;
+
+    for (n = 1; n < N; n++)
+    {
+        double r, w;
+
+        r = R(n);
+        w = pow(r, p);
+        f[n] = f[n-1] + w;
+        sum += w;
+    }
+
+    for (n = 1; n < N; n++) f[n] /= sum;
+
+    /* generate bonds */
+
+    while (num_bonds < NUM_BONDS)
+    {
+        int i, j;
+        double r;
+
+        i = (int) (N * RNG_UNIFORM(rand));
+        r = RNG_UNIFORM(rand);
+        n = bisect(f, r, 1, N);
+        j = i + n;
+
+        if (j >= N)
+        {
+            int _i = i;
+            i = j - N;
+            j = _i;
+        }
+
+        if (! bonds_have(b, i, j))
+        {
+            double v = RNG_GAUSS(rand, 1.);
+            bonds_add(b, i, j, v);
+            num_bonds++;
+        }
+    }
+
+    free(f);
+}
+

src/bond_dist.h

@@ -0,0 +1,4 @@
+#include "bonds.h"
+
+void bond_dist_dilute(bonds *b, double sigma, rng *rand);
+

src/bonds.c

@@ -0,0 +1,69 @@
+#include "bonds.h"
+
+void bonds_init(bonds *b)
+{
+    int i;
+
+    for (i = 0; i < N; i++)
+        b->head[i] = NULL; 
+}
+
+void bonds_add(bonds *b, int i, int j, double v)
+{
+    node *new = malloc(sizeof(node));
+    new->j = j;
+    new->v = v;
+    new->next = b->head[i];
+    b->head[i] = new;
+}
+
+int bonds_have(bonds *b, int i, int j)
+{
+    node *next = b->head[i];
+
+    while (next)
+    {
+        if (next->j == j)
+            return 1;
+
+        next = next->next;
+    }
+
+    return 0;
+}
+
+void bonds_to_list(bonds *b, int *bond_site, double *bond_value)
+{
+    int i;
+
+    for (i = 0; i < N; i++)
+    {
+        node *next = b->head[i];
+
+        while (next)
+        {
+            *(bond_site++) = i;
+            *(bond_site++) = next->j;
+            *(bond_value++) = next->v;
+            next = next->next;
+        }
+    }
+}
+
+void bonds_free(bonds *b)
+{
+    int i;
+
+    for (i = 0; i < N; i++)
+    {
+        node *this = b->head[i];
+
+        while (this)
+        {
+            node *next = this->next;
+            free(this);
+            this = next;
+        }
+    }
+}
+

src/bonds.h

@@ -0,0 +1,26 @@
+#ifndef BONDS_H
+#define BONDS_H
+
+#include "rng.h"
+#include "sample.h"
+
+typedef struct
+{
+    int j;
+    double v;
+    void *next;
+} node;
+
+typedef struct { node *head[N]; } bonds;
+
+void bonds_init(bonds *b);
+
+void bonds_add(bonds *b, int i, int j, double v);
+
+int bonds_have(bonds *b, int i, int j);
+
+void bonds_to_list(bonds *b, int *bond_site, double *bond_value);
+
+void bonds_free(bonds *b);
+
+#endif

src/exchange.h.in

@@ -1,10 +1,10 @@
+#ifndef EXCHANGE_H
+#define EXCHANGE_H
+
 #include "sample.h"
 #include "replica.h"
 #include "rng.h"
 
-#ifndef EXCHANGE_H
-#define EXCHANGE_H
-
 #define NUM_REPLICAS        @NUM_REPLICAS@
 #define REPLICA_EXCHANGE    @REPLICA_EXCHANGE@
 

src/isg.c

@@ -1,5 +1,8 @@
 #include "isg.h"
+#include "bonds.h"
+#include "bond_dist.h"
 #include "modules.h"
+#include "rng.h"
 #include <assert.h>
 #include <ctype.h>
 #include <math.h>
@@ -112,11 +115,39 @@ int main(int argc, char *argv[])
         exit(EXIT_FAILURE);
     }
 
-    /* load sample */
     seed = strtol(seeds[rank], NULL, 16);
+
+#if GENERATE_SAMPLE
+    /* generate a new sample */
+    {
+        bonds *b = malloc(sizeof(bonds));
+        rng *rand = RNG_ALLOC();
+        int *bond_site;
+        double *bond_value;
+
+        bonds_init(b);
+        RNG_SEED(rand, seed);
+        bond_dist_dilute(b, SIGMA, rand);
+        RNG_FREE(rand);
+
+        bond_site =  malloc(2*NUM_BONDS*sizeof(int));
+        bond_value = malloc(  NUM_BONDS*sizeof(double));
+
+        bonds_to_list(b, bond_site, bond_value);
+        bonds_free(b);
+        free(b);
+
+        sample_init(&s->sample, bond_site, bond_value);
+
+        free(bond_site);
+        free(bond_value);
+    }
+#else
+    /* load sample from text file */
     OPEN_FILE("samp" SUFFIX, seed, buf, fp_in, "r");
     sample_init_read(&s->sample, fp_in);
     fclose(fp_in);
+#endif
 
     /* open output files */
 #define OPEN_OUTPUT(name, dummy) \

src/isg.h.in

@@ -4,6 +4,10 @@
 #define LOG_UPDATES_PER_MEAS @LOG_UPDATES_PER_MEAS@
 #define UPDATES_PER_MEAS @UPDATES_PER_MEAS@
 
+#define GENERATE_SAMPLE @GENERATE_SAMPLE@   /* if true, generate samples on the fly,
+                                               otherwise read from text files */
+#define SIGMA @SIGMA@
+
 #define TEMP_FILE "temps.txt"
 #define SUFFIX "_%06x.txt"
 

src/rng.h.in

@@ -22,6 +22,12 @@ typedef dsfmt_t rng;
 /* use an RNG from the GSL... */
 
 #include "gsl_rng.h"
+
+/* XXX Hack because gsl_randist uses `N` as variable... */
+#undef N
+#include "gsl_randist.h"
+#define N @N@
+
 typedef gsl_rng rng;
 #if   RNG == RNG_MT19937
 #define RNG_TYPE gsl_rng_mt19937
@@ -34,6 +40,7 @@ typedef gsl_rng rng;
 #define RNG_ALLOC()         gsl_rng_alloc(RNG_TYPE)
 #define RNG_SEED(r, seed)   gsl_rng_set(r, seed)
 #define RNG_UNIFORM(r)      gsl_rng_uniform(r)
+#define RNG_GAUSS(r, sigma) gsl_ran_gaussian(r, sigma)
 #define RNG_FREE(r)         gsl_rng_free(r)
 
 #endif  /* GSL RNGs*/

src/sample.c

@@ -22,7 +22,7 @@ void init_sample(int const *bond_site, double const *bond_value,
     int i, b;
 
     /* 1st pass: calculate vertex degrees and offsets */
-    memset(z, 0, N*sizeof(int));
+    memset(z, 0, N * sizeof(int));
 
     for (i = 0; i < NZ; i++)
         z[bond_site[i]]++;
@@ -33,8 +33,8 @@ void init_sample(int const *bond_site, double const *bond_value,
         offset[i] = offset[i-1] + z[i-1];
 
     /* 2nd pass: fill in values */
-    memset(z, 0, N*sizeof(int));
-    memset(h2m, 0., N*sizeof(double));
+    memset(z, 0, N * sizeof(int));
+    memset(h2m, 0., N * sizeof(double));
     *um = 0.;
 
     for (i = 0; i < NUM_BONDS; i++)
@@ -63,8 +63,8 @@ void sample_init_read(sample_data *s, FILE *fp)
     int *bond_site, i;
     double *bond_value;
 
-    bond_site  = malloc(NZ*sizeof(int));
-    bond_value = malloc(NUM_BONDS*sizeof(double));
+    bond_site  = malloc(NZ * sizeof(int));
+    bond_value = malloc(NUM_BONDS * sizeof(double));
 
     for (i = 0; i < NUM_BONDS; i++)
     {

src/sample.h.in

@@ -20,6 +20,7 @@ typedef struct
     int offset[N];
     double h2m[N];
     double um;
+
 } sample_data;
 
 void sample_init(sample_data *s,

submit.sh

@@ -29,22 +29,19 @@ askp jobscript "$root/job.sh" "job script"
 askp commit $(cd $root; git rev-parse HEAD) "commit"
 
 askp LOG_N 7 "log_2 N"
-askp sigma "0.600" "sigma"
-askp DILUTE 1 "dilute model"
-
-if [ "$DILUTE" -eq 1 ]; then
-    askp Z 6 "average z"
-    askp Z_MAX 24 "maximum z"
-else setp Z 0
-fi
+askp SIGMA "0.600" "sigma"
+askp Z 6 "average z"
 
 askp tempset "./temps.txt" "temperature set"
 askp LOG_WARMUP_UPDATES 3 "number of warmup decades"
 askp LOG_NUM_UPDATES 16 "number of decades"
 askp log_min_meas_per_bin 10 "log minimum measurements per histogram bin"
 askp REPLICA_EXCHANGE 0 "enable replica exchange"
 askp FULL_OUTPUT 0 "full output [O(N) space]"
-askp KEEP_SAMPLES 0 "keep samples"
+askp GENERATE_SAMPLE 0 "generate samples at run-time?"
+if [ $GENERATE_SAMPLE -eq 0 ]; then
+    askp KEEP_SAMPLES 0 "keep samples?"
+fi
 askp num_cores 1 "number of cores per job"
 askp num_jobs 100 "number of jobs"
 askp h_rss "4G" "reserve memory"