}
-// 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;
+void initialize_Ising(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, Species_parray &print_spec)
+{
Hamiltonian H(r.couplingsInfile);
H.set_temp(r.bh, r.bJ);
- for (unsigned n = 0; n < r.num_runs; n++) {
-
// initialize virus
VirusSpecies *s1 = new VirusSpecies("I");
species.push_back(s1);
s1->pop[Virus(H)] = s1->count; // default mu = 6.0e-5
}
- Species_parray print_spec;
print_spec.push_back(s1);
if (r.useEpitope) importEpitope(r,species,reactions,s1,print_spec); // <infile>.ep
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)
+ ;
+ else
+ initialize_Ising(r,species,reactions,print_spec);
simulate(reactions, species, r.t_end, r.max_steps, r.sample_interval, print_spec, r.useVerbose);
std::cout << '\t' << (*spec)->count;
if (verbose) {
- VirusSpecies *V = dynamic_cast<VirusSpecies*> (*spec);
- if (V) {
+ if (VirusSpecies *V = dynamic_cast<VirusSpecies*> (*spec)) {
virus_map::iterator it = V->pop.begin(),
end = V->pop.end();
for (; it != end; ++it) {
<< '(' << it->first.mutated_sites << ')';
}
}
+ else if (NTVirusSpecies *V = dynamic_cast<NTVirusSpecies*> (*spec)) {
+ NTVirus_map::iterator it = V->pop.begin(),
+ end = V->pop.end();
+ for (; it != end; ++it) {
+ std::cout << ' ' << it->second
+ << '(' << it->first.nt_seq << ')';
+ }
+ }
}
}
"Command line rules:\n"
"\n"
" -d (string) working directory\n"
-" -i (string) input couplings file\n"
+" -i (string) input couplings file (*.j)\n"
+" -I (string) Potts state <-> amino acid mapping file (*-seq2state.dat)\n"
" -o (string) output file stem\n"
-" -s (string) input state file, containing initial population fraction\n"
-" -ep (string) input file containing targeted epitope\n"
+" -s (string) input state file, containing initial population fraction (*.st)\n"
+" -ep (string) input file containing targeted epitope (*.ep)\n"
" -epN0 (int) initial number of naive CTL (corresp. to latest epitope)\n"
" -epM0 (int) initial number of memory CTL (corresp. to latest epitope)\n"
" -epT0 (int) initial number of effector CTL (corresp. to latest epitope)\n"
" -mu (double) mutation rate\n"
" -b (double) \"inverse temperature\" multiplier\n"
" -v flag for verbose output\n"
+" -potts flag for nucleotide/amino-acid/Potts states\n"
" -numruns (int) number of trajectories to simulate\n"
" -rate_X (double) X in [s,u,p,a,r,d,dprime,b,e,w,aprime,g,h]\n"
; std::cout.flush();
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],"-I")==0) { if (++i==argc) break; else r.seq2stateInfile=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],"-ep")==0) { if (++i==argc) break; else { r.epitopefiles.push_back(argv[i]); r.useEpitope=true; ++r.num_CTL_clones; r.init_CTL_numN.push_back(0); r.init_CTL_numM.push_back(0); r.init_CTL_numT.push_back(0); } }
else if (strcmp(argv[i],"-v")==0) { r.useVerbose=true; }
+ else if (strcmp(argv[i],"-potts")==0) { r.usePotts=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],"-seed")==0) { if (++i==argc) break; else seed=(unsigned)strtodouble(argv[i]); }
double rate_aprime;
double rate_g;
double rate_h;
+
+ bool usePotts;
+ std::string seq2stateInfile;
// constructor sets default parameters
rate_g = 0.03;
rate_h = 0.03;
+ usePotts = false;
+
}
supplementaryOutfile=directory+"/"+outfile+".sum";
stateInfile=directory+"/"+statefile+".st";
epitopeInfile=directory+"/"+epitopefile+".ep";
+ seq2stateInfile=directory+"/"+seq2stateInfile+"-seq2state.dat";
for (size_t i = 0; i != epitopefiles.size(); ++i) {
epitopefiles[i] = directory+"/"+epitopefiles[i]+".ep";
supplementaryOutfile=outfile+".sum";
stateInfile=statefile+".st";
epitopeInfile=epitopefile+".ep";
+ seq2stateInfile=seq2stateInfile+"-seq2state.dat";
for (size_t i = 0; i != epitopefiles.size(); ++i) {
epitopefiles[i] = epitopefiles[i]+".ep";
void importState(RunParameters_SS &r);
void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *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);
#endif