Initialize distribution of nt seqs from file. Print (aggregate) unique amino acid...

diff --git a/reaction.cpp b/reaction.cpp
index 99f0daf..c95d419 100644 (file)
--- a/reaction.cpp
+++ b/reaction.cpp
@@ -280,12 +280,12 @@ void NTVirusReaction::fire()
         unsigned siteNT = *it;
 
         // nucleotide to codon
-        int oldNT = str2nt.at(v.nt_seq[siteNT]);  // silent conversion from enum nt
-        int newNT;
+        nt oldNT = str2nt.at(v.nt_seq[siteNT]);  // silent conversion from enum nt
+        nt newNT;
         do {    // pick rand on {0,1,2,3} until we get a new one
-            newNT = gsl_rng_uniform_int(rnd,4-1);
+            newNT = static_cast<nt>(gsl_rng_uniform_int(rnd,4-1));
         } while (newNT == oldNT);
-        v.nt_seq[siteNT] = newNT;
+        v.nt_seq[siteNT] = nt2str.at(newNT);
 
         // codon to amino acid
         unsigned siteAA = siteNT / 3;

diff --git a/ss.cpp b/ss.cpp
index c06481b..6f73581 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -80,6 +80,49 @@ void initialize_Ising(RunParameters_SS &r, Species_parray &species, Rxn_parray &
 
 }
 
+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
 
@@ -111,7 +154,7 @@ void run(RunParameters_SS &r, unsigned seed) {
     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);
 
@@ -174,12 +217,24 @@ void print_species_traj(double t, const Species_parray &print_spec, bool verbose
                 }
             }
             else if (NTVirusSpecies *V = dynamic_cast<NTVirusSpecies*> (*spec)) {
+                // aggregate and print unique amino acid sequences
+                std::map<std::string,unsigned int> aa_seqs;
                 NTVirus_map::iterator  it = V->pop.begin(),
                                       end = V->pop.end();
                 for (; it != end; ++it) {
-                    std::cout << ' ' << it->second
-                              << '(' << it->first.nt_seq << ')';
+                    //std::cout << ' ' << it->second
+                    //          << '(' << aaseq2str(it->first.aa_seq) << ')';
+                    std::string s = aaseq2str(it->first.aa_seq);
+                    if (aa_seqs.count(s)==0) aa_seqs[s]  = it->second;  // new sequence
+                    else                     aa_seqs[s] += it->second;  // add to existing pool
+                }
+                std::map<std::string,unsigned int>::iterator ss = aa_seqs.begin(),
+                                                             send = aa_seqs.end();
+                for (; ss != send; ++ss) {
+                    std::cout << ' ' << ss->second
+                              << '(' << ss->first << ')';
                 }
+
             }
         }
 
@@ -303,6 +358,42 @@ 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';
+    }
+    
+}
+
+
+
+// load epitope definitions from file
+
+void importEpitope_Potts(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, NTVirusSpecies *V, Species_parray &print_spec)
+{
+
+}
+
+
 // load epitope definitions from file
 
 void importEpitope(RunParameters_SS &r, Species_parray &species, Rxn_parray &reactions, VirusSpecies *V, Species_parray &print_spec)

diff --git a/ss.h b/ss.h
index 39539fa..b52d6f5 100644 (file)
--- a/ss.h
+++ b/ss.h
@@ -158,6 +158,7 @@ public:
     
     bool usePotts;
     std::string seq2stateInfile;
+    std::vector<std::string> initPop_NT;
 
 
     // constructor sets default parameters
@@ -244,9 +245,12 @@ 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