Start runtime framework for executing Potts system. Break out initialization routine.
authorDariusz Murakowski <murakdar@mit.edu>
Fri, 24 Apr 2015 08:28:52 +0000 (04:28 -0400)
committerDariusz Murakowski <murakdar@mit.edu>
Fri, 24 Apr 2015 08:33:36 +0000 (04:33 -0400)
ss.cpp
ss.h

diff --git a/ss.cpp b/ss.cpp
index 35396f5..c06481b 100644 (file)
--- 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);    // <infile>.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<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) {
@@ -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<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 << ')';
+                }
+            }
         }
 
     }
@@ -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 (file)
--- 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