From e4edbe688d7cc029752b176aee3433072e6ff6b0 Mon Sep 17 00:00:00 2001
From: Dariusz Murakowski <murakdar@mit.edu>
Date: Wed, 15 May 2013 10:10:01 -0400
Subject: [PATCH] More output for non-epitope runs. Assuage compiler warnings
 of unused parameters and signed/unsigned int comparison.

---
 Makefile        |  2 +-
 hamiltonian.cpp | 10 +++++--
 population.cpp  |  2 +-
 virus.cpp       |  2 +-
 wf.cpp          | 85 ++++++++++++++++++++++++++++++++++++++++++++++-----------
 wf.h            |  2 +-
 6 files changed, 81 insertions(+), 22 deletions(-)

diff --git a/Makefile b/Makefile
index 231fde7..e7f3816 100644
--- a/Makefile
+++ b/Makefile
@@ -3,7 +3,7 @@ SRCS = wf.cpp hamiltonian.cpp population.cpp virus.cpp
 EXECNAME = wf
 
 CXX = c++
-CFLAGS = $(DBGFLAG) -Wall -Wextra -pedantic
+CFLAGS = $(DBGFLAG) -Wall -Wextra -pedantic -Wno-unused-parameter
 INCLUDEDIR = -I/usr/local/pkg/gsl/gsl-1.15/include
 # = `gsl-config --cflags`
 LIBDIR = -L/usr/local/pkg/gsl/gsl-1.15/lib
diff --git a/hamiltonian.cpp b/hamiltonian.cpp
index b5b94ba..842504a 100644
--- a/hamiltonian.cpp
+++ b/hamiltonian.cpp
@@ -52,6 +52,7 @@ Hamiltonian::Hamiltonian(std::string &FILENAME) {
                         JIndex[inputKey]=neighborJ;
                         
                 }
+                //std::cout << site << " " << neighborSite << " " << neighborJ << "\n";
                 
 			}
             
@@ -92,6 +93,11 @@ Hamiltonian::Hamiltonian(std::string &FILENAME) {
 		}
         
 	}
+
+    //for(int i=0; i<size; i++)
+        //for(int j=0; j<size; j++)
+            //std::cout << i << " " << j << " " << J[i][j] << "\n";
+            //printf("%d %d %lf\n",i,j,J[i][j]);
     
 }
 
@@ -237,7 +243,7 @@ bool EpitopeHamiltonian::escaped(const std::set<unsigned int> &mutated_sites) co
 
     bool escape=false;
     
-    for (int i=0;i<epitopeWT.size();i++) {
+    for (unsigned i=0;i<epitopeWT.size();i++) {
     
         if (mutated_sites.count(epitopeWT[i])==1) { escape=true; break; }
         
@@ -245,7 +251,7 @@ bool EpitopeHamiltonian::escaped(const std::set<unsigned int> &mutated_sites) co
     
     if (!escape) {
     
-        for (int i=0;i<epitopeMut.size();i++) {
+        for (unsigned i=0;i<epitopeMut.size();i++) {
         
             if (mutated_sites.count(epitopeMut[i])==0) { escape=true; break; }
             
diff --git a/population.cpp b/population.cpp
index a222a87..3c3cdc2 100644
--- a/population.cpp
+++ b/population.cpp
@@ -28,7 +28,7 @@ Population::Population(const Hamiltonian &H, unsigned int N, double mu) {
 
 Population::Population(const Hamiltonian &H, unsigned int N, double mu, const std::vector<std::set<unsigned int> > &initPop, const std::vector<double> &initFrac) {
 	
-    for (int i=0;i<initPop.size();i++) {
+    for (unsigned i=0;i<initPop.size();i++) {
     
         Virus V(H, mu, initPop[i]);
         pop[V] = (unsigned int) N * initFrac[i];
diff --git a/virus.cpp b/virus.cpp
index 007a5d1..fd69023 100644
--- a/virus.cpp
+++ b/virus.cpp
@@ -51,7 +51,7 @@ void Virus::mutate(const Hamiltonian &H, gsl_rng* r) {
     
 	while (n<1) n=gsl_ran_binomial(r,mu,L);
     
-    for (int i=0;i<n;i++) {
+    for (unsigned i=0;i<n;i++) {
     
         unsigned int x=gsl_rng_uniform_int(r,H.size);
         
diff --git a/wf.cpp b/wf.cpp
index 8c14ef8..c89e797 100644
--- a/wf.cpp
+++ b/wf.cpp
@@ -13,6 +13,27 @@
 #include "wf.h"
 
 
+// generate a random number generator seed from /dev/urandom
+// borrowed from SSC src/runtime/ssc_rts.c
+static unsigned sim_random_seed() {
+    unsigned random_seed;
+    FILE *rng = fopen("/dev/urandom", "r");
+    if (rng == NULL) {
+        fprintf(stderr,
+                "error: cannot read /dev/urandom randomness generator\n");
+        exit(1);
+    }
+    size_t num_random_read = fread(&random_seed, sizeof(random_seed), 1, rng);
+    if (num_random_read != 1) {
+        fprintf(stderr,
+                "error: cannot read /dev/urandom randomness generator\n");
+        exit(1);
+    }
+    fclose(rng);
+    return random_seed;
+}
+
+
 // Run the program
 
 void run(RunParameters &r) {
@@ -21,14 +42,15 @@ void run(RunParameters &r) {
     
     // Initialize RNG and set initial state, if importing from file
     
-    static gsl_rng *rnd = gsl_rng_alloc(gsl_rng_taus);	//pointer to random number generator state
+    static gsl_rng *rnd = gsl_rng_alloc(gsl_rng_taus2);	//pointer to random number generator state
     
-    srand((unsigned)time(0));
-    gsl_rng_set(rnd,rand());
+    //srand((unsigned)time(0));
+    //gsl_rng_set(rnd,rand());
+    gsl_rng_set(rnd,sim_random_seed());
     
     r.setFiles();
-    //FILE *popout=fopen(r.trajectoryOutfile.c_str(),"w");
-    FILE *supout=fopen(r.supplementaryOutfile.c_str(),"w");
+    FILE *popout=fopen(r.trajectoryOutfile.c_str(),"w");        // <outfile>.dat
+    FILE *supout=fopen(r.supplementaryOutfile.c_str(),"w");     // <outfile>.sum
     
     if (r.importState) importState(r);
     
@@ -48,7 +70,7 @@ void run(RunParameters &r) {
             for (i=0; i<r.g; ++i) {
         
                 P.next_generation(H, rnd, r.useRelative, r.useVerbose);
-                //P.write_population(popout,i);
+                P.write_population(popout,i);
             
                 if (r.runUntilEscape && P.escaped(H)) break;
             
@@ -71,22 +93,32 @@ void run(RunParameters &r) {
     // Run (w/out targeted epitope)
     
     else {
-    
-        Hamiltonian H(r.couplingsInfile);
-        H.set_temp(r.bh, r.bJ);
-        Population P(H, r.n, r.mu);
-        
-        for (unsigned int i=0; i<r.g; ++i) {
+
+        //for (unsigned int n=0;n<r.num_runs;n++) {
         
-            P.next_generation(H, rnd, r.useRelative, r.useVerbose);
-            //P.write_population(popout,i);
-                        
-        }
+            Hamiltonian H(r.couplingsInfile);
+            H.set_temp(r.bh, r.bJ);
+            Population P(H, r.n, r.mu);
+            
+            unsigned int i;
+            for (i=0; i<r.g; ++i) {
+            
+                P.next_generation(H, rnd, r.useRelative, r.useVerbose);
+                P.write_population(popout,i);
+
+                if (r.runUntilEscape && P.escaped(H)) break;
+
+            }
+            
+        //}
         
     }
     
     gsl_rng_free(rnd);  //Free up memory from random number generator
 
+    fclose(popout);     // close file handles
+    fclose(supout);
+
 }
 
 
@@ -151,6 +183,25 @@ void importState(RunParameters &r) {
 
 }
 
+void usage()
+{
+    std::cout <<
+"Command line rules:\n"
+"\n"
+" -d  (string)        working directory\n"
+" -i  (string)        input couplings file\n"
+" -o  (string)        output file stem\n"
+" -s  (string)        input state file, containing initial population fraction and targeted epitope\n"
+" -n  (int/double)    population size\n"
+" -g  (int/double)    number of generations\n"
+" -mu (double)        mutation rate\n"
+" -b  (double)        \"inverse temperature\" multiplier\n"
+" -e  (int) (int)     targeted epitope spans sites [e#1, ..., e#2], assuming all WT\n"
+" -r                  flag to use relative energy condition for survival probability\n"
+" -esc                flag to run until escape observed (or max number of generations reached)\n"
+" -v                  flag for verbose output\n";
+    std::cout.flush();
+}
 
 /*
  * Command line rules:
@@ -197,6 +248,8 @@ int main(int argc, char *argv[]) {
         else if (strcmp(argv[i],"-r")==0)   { r.useRelative=true;    }
         else if (strcmp(argv[i],"-esc")==0) { r.runUntilEscape=true; }
         else if (strcmp(argv[i],"-v")==0)   { r.useVerbose=true;     }
+
+        else if (strcmp(argv[i],"-h")==0)   { usage(); return 0;    }
         
         else printf("Unrecognized command! '%s'\n",argv[i]);
                 
diff --git a/wf.h b/wf.h
index f437a2e..7250345 100644
--- a/wf.h
+++ b/wf.h
@@ -113,7 +113,7 @@ public:
         n  = 100000;
         g  = 150;
         num_runs = 1000;
-        mu = 6.0 * pow(10.0,-5.0);
+        mu = 6.0e-5;
         bh = 1.0;
         bJ = 1.0;
         
-- 
2.7.4