From: Dariusz Murakowski <murakdar@mit.edu>
Date: Tue, 28 Apr 2015 21:29:40 +0000 (-0400)
Subject: Merge branch 'master' into ss_send
X-Git-Url: http://src.murakowski.org/?a=commitdiff_plain;h=9b086d97a9d2022f2688a304b789599067007da7;p=VirEvoDyn.git

Merge branch 'master' into ss_send
---

9b086d97a9d2022f2688a304b789599067007da7
diff --cc ss.cpp
index c4fecad,56f26bd..a1dc318
--- a/ss.cpp
+++ b/ss.cpp
@@@ -103,6 -78,85 +77,43 @@@ void initialize_Ising(RunParameters_SS 
      r2->V = s1;
      reactions.push_back(r2);
  
+ }
+ 
 -void initialize_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, Species_parray &print_spec)
 -{
 -
 -    std::cout << "initializing potts\n";
 -
 -    PottsHamiltonian H(r.couplingsInfile,r.seq2stateInfile);
 -    H.set_temp(r.bh, r.bJ);
 -
 -    // initialize virus
 -    NTVirusSpecies *s1 = new NTVirusSpecies("I");
 -    species.push_back(s1);
 -    if (r.importState) {
 -        importState_Potts(r);      // <infile>.st
 -        for (size_t i=0; i<r.initPop_NT.size(); i++) {
 -            NTVirus V(H, r.initPop_NT[i]);
 -            unsigned int N = (unsigned int) r.n * r.initFrac[i];
 -            s1->pop[V] = N;
 -            s1->count += N;
 -        }
 -    }
 -    else {
 -        //s1->count = r.n;
 -        //s1->pop[Virus(H)] = s1->count;   // default mu = 6.0e-5
 -    }
 -
 -    print_spec.push_back(s1);
 -
 -    if (r.useEpitope)  importEpitope_Potts(r,species,reactions,s1,print_spec);    // <infile>.ep
 -
 -    // V -> V + V'
 -    NTVirusReaction* r1 = new NTVirusReaction(r.rate_s);  // s_k
 -    r1->mu = r.mu;
 -    r1->H = H;
 -    r1->V = s1;
 -    reactions.push_back(r1);
 -    // V -> 0
 -    NTVirusDeathReaction* r2 = new NTVirusDeathReaction(r.rate_u);  // u
 -    r2->H = H;
 -    r2->V = s1;
 -    reactions.push_back(r2);
 -
 -}
+ 
+ 
+ // Run the program
+ 
+ void run(RunParameters_SS &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
+     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);
+ 
+     if (r.t_end == HUGE_VAL && r.max_steps == 0 && r.sample_interval == HUGE_VAL) {
+         std::cerr << "ERROR: must specify one of -e or -t or -E; see -h for usage" << std::endl;
+         exit(1);
+     }
+     if (r.sample_interval == HUGE_VAL)
+         r.sample_interval = r.t_end;  // ensure last time will be printed
+ 
+     r.setFiles();
+ 
+     // master arrays of (pointers to) species and reactions
+     Species_parray species;
+     Rxn_parray reactions;
+     Species_parray print_spec;
+ 
+     for (unsigned n = 0; n < r.num_runs; n++) {
+ 
+     if (r.usePotts)
 -        initialize_Potts(r,species,reactions,print_spec);
++        ;
+     else
+         initialize_Ising(r,species,reactions,print_spec);
  
      simulate(reactions, species, r.t_end, r.max_steps, r.sample_interval, print_spec, r.useVerbose);
  
@@@ -285,6 -359,41 +295,15 @@@ void importState(RunParameters_SS &r) 
  }
  
  
 -// Import initial state from a state file
 -
 -void importState_Potts(RunParameters_SS &r)
 -{
 -    std::ifstream input(r.stateInfile.c_str());   // <infile>.st
 -    if (input == NULL) { perror((std::string("ERROR in importState_Potts: ") + r.stateInfile).c_str()); exit(1); }
 -    
 -    double frac;
 -    std::string nt_str_seq;
 -
 -    std::string line;
 -    while (std::getline(input,line)) {
 -        std::stringstream readStrStrm(line);
 -        readStrStrm >> frac;
 -        readStrStrm >> nt_str_seq;
 -
 -        r.initFrac.push_back(frac);
 -        r.initPop_NT.push_back(nt_str_seq);
 -    }
 -
 -    for (size_t i=0; i<r.initFrac.size(); i++) {    // XXX
 -        std::cout << r.initFrac[i] << ' ' << r.initPop_NT[i] << '\n';
 -    }
 -    
 -}
 -
+ 
+ // helper function adapted from /usr/include/c++/4.4.7/bits/basic_string.h
+   inline std::string
+   to_string_uint(unsigned __val)
+   { return __gnu_cxx::__to_xstring<std::string>(&std::vsnprintf,
+                        4 * sizeof(unsigned),
+                        "%u", __val); }
+ 
+ 
  // load epitope definitions from file
  
  void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V, Species_parray &print_spec)
diff --cc ss.h
index 5738ecd,d971f19..ffc6fad
--- a/ss.h
+++ b/ss.h
@@@ -237,8 -267,12 +267,9 @@@ public
  
  void run(RunParameters_SS &r);
  void importState(RunParameters_SS &r);
 -void importState_Potts(RunParameters_SS &r);
  void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V, Species_parray &print_spec);
 -void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, NTVirusSpecies *V, Species_parray &print_spec);
  void simulate(Rxn_parray &reactions, Species_parray &species, double t_end, long max_steps, double sample_interval, const Species_parray &print_spec, bool verbose);
+ void initialize_Ising(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, Species_parray &print_spec);
 -void initialize_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, Species_parray &print_spec);
  
  
  #endif