From: Dariusz Murakowski <murakdar@mit.edu>
Date: Thu, 13 Mar 2014 23:07:53 +0000 (-0400)
Subject: Start of Monte Carlo (static epitope calculation), copy WF code.
X-Git-Url: http://src.murakowski.org/?a=commitdiff_plain;h=6794bd0860cb18f6917e7c1e5e126b0076f8ce32;p=VirEvoDyn.git

Start of Monte Carlo (static epitope calculation), copy WF code.
---

diff --git a/Makefile b/Makefile
index 63367d5..2d81891 100644
--- a/Makefile
+++ b/Makefile
@@ -2,6 +2,10 @@
 SRCS = wf.cpp hamiltonian.cpp population.cpp virus.cpp
 EXECNAME = wf
 
+SRCS_MC = mc.cpp hamiltonian.cpp population.cpp virus.cpp
+EXECNAME_MC = mc
+
+
 CXX = c++
 CFLAGS = $(DBGFLAG) -Wall -Wextra -pedantic -Wno-unused-parameter
 INCLUDEDIR = -I/usr/local/pkg/gsl/gsl-1.15/include
@@ -21,12 +25,15 @@ endif
 
 .PHONY: clean
 
-all: $(EXECNAME)
+all: $(EXECNAME) $(EXECNAME_MC)
 #	@echo done making $(EXECNAME)
 
 $(EXECNAME): $(SRCS)
 	$(CXX) $(SRCS) $(CFLAGS) $(LIBDIR) $(INCLUDEDIR) -o $(EXECNAME) $(LIBS)
 
+$(EXECNAME_MC): $(SRCS_MC)
+	$(CXX) $(SRCS_MC) $(CFLAGS) $(LIBDIR) $(INCLUDEDIR) -o $(EXECNAME_MC) $(LIBS)
+
 # concatenate all the source files before compiling
 #	./concat.sh $(EXECNAME) $(SRCS)
 #	$(CXX) -x c++ $(EXECNAME).combined $(CFLAGS) $(LIBDIR) $(INCLUDEDIR) -o $(EXECNAME) $(LIBS)
@@ -39,6 +46,6 @@ $(EXECNAME): $(SRCS)
 #	$(CXX) -c $(CFLAGS) $(INCLUDES) $< -o $@
 
 clean:
-	$(RM) *.o $(EXECNAME)
+	$(RM) *.o $(EXECNAME) $(EXECNAME_MC)
 #	$(RM) *.o *.combined $(EXECNAME)
 
diff --git a/mc.cpp b/mc.cpp
new file mode 100644
index 0000000..d0ddcf0
--- /dev/null
+++ b/mc.cpp
@@ -0,0 +1,403 @@
+#include <iostream>
+#include <string>
+#include <fstream>
+#include <cstdio>
+#include <cstdlib>
+#include <cmath>
+
+#include <gsl/gsl_rng.h>
+
+#include "hamiltonian.h"
+#include "population.h"
+#include "virus.h"
+#include "mc.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, unsigned seed) {
+    
+    // Initialize RNG and set initial state, if importing from file
+    
+    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());
+    gsl_rng_set(rnd,seed);
+    
+    r.setFiles();
+    FILE *popout=fopen(r.trajectoryOutfile.c_str(),"w");        // <outfile>.dat
+    FILE *supout=r.useTwoSite ? NULL : fopen(r.supplementaryOutfile.c_str(),"w");     // <outfile>.sum
+    if (popout == NULL) { perror("File error"); exit(1); }
+    if (r.useTwoSite && supout == NULL) { perror("File error"); exit(1); }
+    
+    if (r.importState) importState(r);      // <infile>.st
+    if (r.useEpitope)  importEpitope(r);    // <infile>.ep
+    fflush(stdout);
+    
+
+    if(r.useTwoSite) {
+
+            fprintf(popout,"gen\tV00\tV10\tV01\tV11\n");
+
+            for (unsigned int n=0;n<r.num_runs;n++) {
+
+                //fprintf(popout,"gen\tV00\tV10\tV01\tV11\n");
+            
+                TwoSiteHamiltonian H(r.h1,r.h2,r.J12);
+                H.set_temp(r.bh, r.bJ);
+                Population P(H, r.n, r.mu, r.initPop, r.initFrac);
+                
+                P.write_two_site_population(popout,H,0);
+
+                unsigned int i;
+                for (i=0; i<r.g; ++i) {
+                
+                    //printf("%d\t%d\n",n,i);
+                    P.next_generation(H, rnd, r.useRelative, r.useVerbose);
+
+                    if ((i+1) % r.write_mod == 0) {
+                        P.write_two_site_population(popout,H,i+1);
+                    }
+
+                    //if (r.runUntilEscape && P.escaped(H)) break;
+
+                }
+                
+            }
+
+    }
+
+
+    else {
+        
+        // Run (w/ targeted epitope)
+        
+        if (r.useEpitope) {
+        
+            for (unsigned int n=0;n<r.num_runs;n++) {
+            
+                EpitopeHamiltonian H(r.couplingsInfile);
+                for (unsigned ep=0; ep<r.numEpitopes; ++ep)
+                    H.set_epitope(r.eWT[ep], r.eMut[ep], r.penalty);
+                H.set_temp(r.bh, r.bJ);
+                Population P(H, r.n, r.mu, r.initPop, r.initFrac);
+
+                // print epitopes
+                /*
+                std::cout << "-----------\n";
+                for (unsigned ep=0; ep<H.penalty.size(); ++ep) {
+                    for (unsigned i=0; i<H.epitopeWT[ep].size(); ++i) {
+                        std::cout << H.epitopeWT[ep][i] << " ";
+                    }
+                    std::cout << "| ";
+                    for (unsigned i=0; i<H.epitopeMut[ep].size(); ++i) {
+                        std::cout << H.epitopeMut[ep][i] << " ";
+                    }
+                    std::cout << "\n";
+                }
+                fflush(stdout);
+                */
+                
+                unsigned int i;
+                for (i=0; i<r.g; ++i) {
+            
+                    P.next_generation(H, rnd, r.useRelative, r.useVerbose);
+
+                    if ((i+1) % r.write_mod == 0) {
+                        P.write_population(popout,i+1);
+                    }
+                
+                    if (r.runUntilEscape && P.escaped(H)) break;
+                
+                }
+                
+                std::set<unsigned int> esc_var;
+                unsigned int esc_num=P.escape_variant(H, esc_var);
+                
+                fprintf(supout,"%d\t%d",i,esc_num);
+                for (std::set<unsigned int>::iterator iter=esc_var.begin(); iter!=esc_var.end(); ++iter) fprintf(supout,"\t%d",*iter);
+                fprintf(supout,"\n");
+                
+                fflush(supout);
+            
+            }
+            
+        }
+        
+        
+        // Run (w/out targeted epitope)
+        
+        else {
+
+            //for (unsigned int n=0;n<r.num_runs;n++) {
+            
+                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
+    if(!r.useTwoSite) fclose(supout);
+
+}
+
+
+// Import initial state from a state file
+
+void importState(RunParameters &r) {
+
+    FILE *input=fopen(r.stateInfile.c_str(),"r");   // <infile>.st
+    if (input == NULL) { perror((std::string("ERROR in importState: ") + r.stateInfile).c_str()); exit(1); }
+    
+    char o;
+    double frac;
+	unsigned int site;
+    
+    // Import initial state of the population
+    
+	while (fscanf(input,"%le",&frac)==1) {
+
+        std::cout << frac; // XXX
+        
+        r.initFrac.push_back(frac);
+        
+        std::set<unsigned int> tempSet;
+        
+		while (fscanf(input,"%c",&o)==1) {
+
+            std::cout << o; // XXX
+			
+			if (o=='\n' || o==';') break;
+			
+			else {
+				
+				fscanf(input,"%u",&site);
+                std::cout << site; // XXX
+                tempSet.insert(site);
+                
+            }
+			
+		}
+        
+        r.initPop.push_back(tempSet);
+        
+        if (o==';') break;
+        
+	}
+
+    std::cout << "\n";  // XXX
+    
+}
+
+
+// load epitope definitions from file
+
+void importEpitope(RunParameters &r) {
+
+    std::ifstream input(r.epitopeInfile.c_str());   // <infile>.ep
+    if (!input) { perror((std::string("ERROR in importEpitope: ") + r.epitopeInfile).c_str()); exit(1); }
+
+    std::string readStr;
+    while (std::getline(input,readStr)) {
+
+        std::string word;
+        unsigned int site;
+        std::vector<unsigned int> tmpEp;
+        size_t pos0 = 0;
+        size_t posBar = readStr.find('|',pos0);
+        size_t posEnd = std::string::npos; //readStr.size();
+
+        // could use std::noskipws to keep tab ('\t') characters?
+
+        std::stringstream readStrStrm(std::string(readStr,pos0,posBar-pos0));
+        while (readStrStrm >> word) {
+            std::cout << word << " "; // XXX
+            std::istringstream i(word);
+            if (i >> site)
+                tmpEp.push_back(site);
+            else    // must have encountered '|'
+                break;
+        }
+        r.eWT.push_back(tmpEp);
+
+        tmpEp.clear();  // reset temp epitope
+        readStrStrm.str("");  readStrStrm.clear();  // reset stream
+
+        std::cout << "| "; // XXX
+
+        readStrStrm.str(std::string(readStr,posBar+1,posEnd/*-posBar-1*/));
+        while (readStrStrm >> word) {
+            std::cout << word << " "; // XXX
+            std::istringstream i(word);
+            if (i >> site)
+                tmpEp.push_back(site);
+            else
+                break;
+        }
+        r.eMut.push_back(tmpEp);
+
+        r.numEpitopes++;
+
+        std::cout << "\n"; // XXX
+
+    }
+
+}
+
+
+// insert into initial state of population
+
+void add_to_two_site_pop(RunParameters &r, bool s1, bool s2, double frac) {
+
+    r.initFrac.push_back(frac);
+
+    std::set<unsigned int> tempSet;
+    if (s1) tempSet.insert(0);
+    if (s2) tempSet.insert(1);
+    r.initPop.push_back(tempSet);
+
+}
+
+
+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\n"
+" -e  (string)        input file containing 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"
+" -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"
+"\n"
+" -2site              flag for two-site two-allele system\n"
+" -h1  (double)       value of field at site 1\n"
+" -h2  (double)       value of field at site 2\n"
+" -J12 (double)       value of coupling between sites 1 and 2\n"
+"                     (note the sign convention on these is P ~ exp(-H),\n"
+"                      H = sum_i h_i s_i + sum_{i<j} J_{ij} s_i s_j\n"
+"                      although internal representation is opposite)\n"
+" -write_mod  (int)   write out state every 'write_mod' generations\n"
+;   std::cout.flush();
+}
+
+/*
+ * Command line rules:
+ *
+ *  -d  (string)        working directory
+ *  -i  (string)        input couplings file
+ *  -o  (string)        output file stem
+ *  -s  (string)        input state file, containing initial population fraction and targeted epitope
+ *  -e  (string)        input file containing targeted epitope
+ *  -n  (int/double)    population size
+ *  -g  (int/double)    number of generations
+ *  -mu (double)        mutation rate
+ *  -b  (double)        "inverse temperature" multiplier
+ *  -r                  flag to use relative energy condition for survival probability
+ *  -esc                flag to run until escape observed (or max number of generations reached)
+ *  -v                  flag for verbose output
+ *  -2site              flag for two-site two-allele system
+ */
+
+int main(int argc, char *argv[]) {
+    
+    RunParameters r;
+
+    unsigned seed = sim_random_seed();
+    
+    // Process command line input
+    
+    for (int i=1;i<argc;i++) {
+        
+        if      (strcmp(argv[i],"-d")==0) { if (++i==argc) break; else r.directory=argv[i];                         }
+        else if (strcmp(argv[i],"-i")==0) { if (++i==argc) break; else r.infile=argv[i];                            }
+        else if (strcmp(argv[i],"-o")==0) { if (++i==argc) break; else r.outfile=argv[i];                           }
+        else if (strcmp(argv[i],"-s")==0) { if (++i==argc) break; else { r.statefile=argv[i]; r.importState=true; } }
+        else if (strcmp(argv[i],"-e")==0) { if (++i==argc) break; else { r.epitopefile=argv[i]; r.useEpitope=true; } }
+        
+        else if (strcmp(argv[i],"-n")==0)  { if (++i==argc) break; else r.n=(unsigned int)strtodouble(argv[i]); }
+        else if (strcmp(argv[i],"-g")==0)  { if (++i==argc) break; else r.g=(unsigned int)strtodouble(argv[i]); }
+        else if (strcmp(argv[i],"-mu")==0) { if (++i==argc) break; else r.mu=strtodouble(argv[i]);              }
+        else if (strcmp(argv[i],"-b")==0)  { if (++i==argc) break; else { r.bh=strtodouble(argv[i]); r.bJ=r.bh; } }
+        else if (strcmp(argv[i],"-bh")==0) { if (++i==argc) break; else r.bh=strtodouble(argv[i]);              }
+        else if (strcmp(argv[i],"-bJ")==0) { if (++i==argc) break; else r.bJ=strtodouble(argv[i]);              }
+        
+        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],"-numruns")==0)  { if (++i==argc) break; else r.num_runs=(unsigned int)strtodouble(argv[i]); }
+
+        else if (strcmp(argv[i],"-2site")==0)   { r.useTwoSite=true;     }
+        else if (strcmp(argv[i],"-h1")==0)  { if (++i==argc) break; else r.h1=-strtodouble(argv[i]);  }
+        else if (strcmp(argv[i],"-h2")==0)  { if (++i==argc) break; else r.h2=-strtodouble(argv[i]);  }
+        else if (strcmp(argv[i],"-J12")==0) { if (++i==argc) break; else r.J12=-strtodouble(argv[i]); }
+
+        else if (strcmp(argv[i],"-p00")==0) { if (++i==argc) break; else add_to_two_site_pop(r,false,false,strtodouble(argv[i])); }
+        else if (strcmp(argv[i],"-p10")==0) { if (++i==argc) break; else add_to_two_site_pop(r,true ,false,strtodouble(argv[i])); }
+        else if (strcmp(argv[i],"-p01")==0) { if (++i==argc) break; else add_to_two_site_pop(r,false,true ,strtodouble(argv[i])); }
+        else if (strcmp(argv[i],"-p11")==0) { if (++i==argc) break; else add_to_two_site_pop(r,true ,true ,strtodouble(argv[i])); }
+
+        else if (strcmp(argv[i],"-seed")==0) { if (++i==argc) break; else seed=(unsigned)strtodouble(argv[i]); }
+
+        else if (strcmp(argv[i],"-write_mod")==0) { if (++i==argc) break; else r.write_mod=(unsigned int)strtodouble(argv[i]); }
+
+        else if (strcmp(argv[i],"-h")==0 || strcmp(argv[i],"--help")==0)   { usage(); return 0;    }
+        
+        else printf("Unrecognized command! '%s'\n",argv[i]);
+                
+    }
+    
+    std::cout << "seed = " << seed << "\n";
+    run(r,seed);
+    
+    return 0;
+    
+}	
+	
diff --git a/mc.h b/mc.h
new file mode 100644
index 0000000..724da81
--- /dev/null
+++ b/mc.h
@@ -0,0 +1,185 @@
+#ifndef MC_H
+#define MC_H
+
+#include <vector>
+#include <set>
+#include <math.h>
+#include <string>
+#include <sstream>
+#include <cstring>
+#include <iostream>
+#include <stdio.h>
+
+
+// STRING MANIPULATION
+
+// Converts generic to string
+
+template <class T>
+
+inline std::string tostring (const T& t) {
+    
+    std::stringstream ss;
+    ss << t;
+    return ss.str();
+    
+}
+
+
+// Converts a string to an integer
+
+inline int strtoint(const std::string &s) {
+    
+    std::istringstream i(s);
+    int x;
+    
+    if (!(i >> x)) printf("String to integer conversion failed!");
+    return x;
+    
+}
+
+
+// Converts a string to an unsigned integer
+
+inline int strtouint(const std::string &s) {
+    
+    std::istringstream i(s);
+    unsigned int x;
+    
+    if (!(i >> x)) printf("String to unsigned integer conversion failed!");
+    return x;
+    
+}
+
+
+// Converts a string to a double
+
+inline double strtodouble(const std::string &s) {
+    
+    std::istringstream i(s);
+    double x;
+    
+    if (!(i >> x)) printf("String to double conversion failed!");
+    return x;
+
+}
+
+
+// PROGRAM SETTINGS
+
+// This class holds the parameters needed for running the simulation
+
+class RunParameters {
+    
+public:
+    
+    std::string directory;      // Path to the directory where the inut file is located
+                                // Output will also be sent to this directory
+    std::string infile;         // Input file from which couplings are to be read
+    std::string outfile;        // Output file (prefix) where data is to be written
+    std::string statefile;      // Input state file which contains initial population fractions
+    std::string epitopefile;    // Input file name for targeted epitopes
+    
+    unsigned int n;             // Population size
+    unsigned int g;             // Number of generations
+    unsigned int num_runs;      // Number of times to run the simulation for gathering overall statistics
+    double mu;                  // Mutation rate
+    double bh;                  // Field "inverse temperature" multiplier
+    double bJ;                  // Coupling "inverse temperature" multiplier
+    int estart;                 // Epitope start position
+    int eend;                   // Epitope end position
+    bool runUntilEscape;        // If true, run the simulation until the population escapes
+    bool importState;           // If true, import state from a state file
+    bool useEpitope;            // Include selective pressure on an epitope 
+    bool useRelative;           // If true, use relative energy for survival probability
+    bool useVerbose;            // If true, print extra information while program is running
+
+    bool useTwoSite;            // If true, use two-site two-allele model
+    double h1;
+    double h2;
+    double J12;
+
+    unsigned int write_mod;     // Write state of the population every __ generations
+    
+    std::vector<std::set<unsigned int> > initPop;   // Initial population sequences
+    std::vector<double>                  initFrac;  // Initial population fraction
+    std::vector<std::vector<unsigned int> >           eWT;       // Sites that are consensus (WT) in the targeted epitope
+    std::vector<std::vector<unsigned int> >           eMut;      // Sites that are mutant in the targeted epitope
+    double                               penalty;   // Energy penalty if sequence contains the targeted epitope
+    unsigned int numEpitopes;   // how many epitopes are targeted
+    
+    std::string couplingsInfile;
+    std::string trajectoryOutfile;
+    std::string supplementaryOutfile;
+    std::string stateInfile;
+    std::string epitopeInfile;
+    
+    
+    RunParameters() {
+    
+        directory="";
+        
+        n  = 100000;
+        g  = 150;
+        num_runs = 1000;
+        mu = 6.0e-5;
+        bh = 1.0;
+        bJ = 1.0;
+        
+        estart = 0;
+        eend   = 0;
+        
+        penalty = 10.0;
+        
+        runUntilEscape=false;
+        importState=false;
+        useEpitope=false;
+        useRelative=false;
+        useVerbose=false;
+
+        useTwoSite = false;
+        h1  = 0.0;
+        h2  = 0.0;
+        J12 = 0.0;
+
+        write_mod = 1;
+
+        numEpitopes = 0;
+        
+    }
+    void setFiles() {
+        
+        if (strcmp(directory.c_str(),"")!=0) {
+        
+            couplingsInfile=directory+"/"+infile+".j";
+            trajectoryOutfile=directory+"/"+outfile+".dat";
+            supplementaryOutfile=directory+"/"+outfile+".sum";
+            stateInfile=directory+"/"+statefile+".st";
+            epitopeInfile=directory+"/"+epitopefile+".ep";
+            
+        }
+        
+        else {
+        
+            couplingsInfile=infile+".j";
+            trajectoryOutfile=outfile+".dat";
+            supplementaryOutfile=outfile+".sum";
+            stateInfile=statefile+".st";
+            epitopeInfile=epitopefile+".ep";
+        
+        }
+        
+        //printf("%s\n%s\n%s\n",couplingsInfile.c_str(),trajectoryOutfile.c_str(),stateInfile.c_str());
+    
+    }
+    ~RunParameters() {}
+    
+};
+
+
+void run(RunParameters &r);
+void importState(RunParameters &r);
+void importEpitope(RunParameters &r);
+
+
+#endif // MC_H