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);
-
-}
+ }
+
- initialize_Potts(r,species,reactions,print_spec);
+
+
+ // 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)
++ ;
+ else
+ initialize_Ising(r,species,reactions,print_spec);
simulate(reactions, species, r.t_end, r.max_steps, r.sample_interval, print_spec, r.useVerbose);
}
-// 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)
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_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, Species_parray &print_spec);
+ void initialize_Ising(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, Species_parray &print_spec);
#endif