From 30f0f1c005cd9ec0dba4385e55fd88fc3b31f27f Mon Sep 17 00:00:00 2001 From: Dariusz Murakowski Date: Fri, 24 Apr 2015 04:28:52 -0400 Subject: [PATCH] Start runtime framework for executing Potts system. Break out initialization routine. --- ss.cpp | 88 +++++++++++++++++++++++++++++++++++++++++------------------------- ss.h | 8 ++++++ 2 files changed, 63 insertions(+), 33 deletions(-) diff --git a/ss.cpp b/ss.cpp index 35396f5..c06481b 100644 --- a/ss.cpp +++ b/ss.cpp @@ -39,37 +39,12 @@ static unsigned sim_random_seed() { } -// 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); @@ -87,7 +62,6 @@ void run(RunParameters_SS &r, unsigned seed) { 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); // .ep @@ -104,6 +78,42 @@ void run(RunParameters_SS &r, unsigned seed) { 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); @@ -155,8 +165,7 @@ void print_species_traj(double t, const Species_parray &print_spec, bool verbose std::cout << '\t' << (*spec)->count; if (verbose) { - VirusSpecies *V = dynamic_cast (*spec); - if (V) { + if (VirusSpecies *V = dynamic_cast (*spec)) { virus_map::iterator it = V->pop.begin(), end = V->pop.end(); for (; it != end; ++it) { @@ -164,6 +173,14 @@ void print_species_traj(double t, const Species_parray &print_spec, bool verbose << '(' << it->first.mutated_sites << ')'; } } + else if (NTVirusSpecies *V = dynamic_cast (*spec)) { + NTVirus_map::iterator it = V->pop.begin(), + end = V->pop.end(); + for (; it != end; ++it) { + std::cout << ' ' << it->second + << '(' << it->first.nt_seq << ')'; + } + } } } @@ -451,10 +468,11 @@ void usage() "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" @@ -465,6 +483,7 @@ void usage() " -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(); @@ -499,6 +518,7 @@ int main(int argc, char *argv[]) { 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); } } @@ -518,6 +538,8 @@ int main(int argc, char *argv[]) { 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]); } diff --git a/ss.h b/ss.h index 5738ecd..39539fa 100644 --- a/ss.h +++ b/ss.h @@ -155,6 +155,9 @@ public: double rate_aprime; double rate_g; double rate_h; + + bool usePotts; + std::string seq2stateInfile; // constructor sets default parameters @@ -192,6 +195,8 @@ public: rate_g = 0.03; rate_h = 0.03; + usePotts = false; + } @@ -204,6 +209,7 @@ public: 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"; @@ -218,6 +224,7 @@ public: 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"; @@ -239,6 +246,7 @@ void run(RunParameters_SS &r); 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 -- 2.7.4