Save mapping of Potts state to amino acid as well. Virus object has no access though.

diff --git a/ham_ss.cpp b/ham_ss.cpp
index 94232d7..df392f3 100644 (file)
--- a/ham_ss.cpp
+++ b/ham_ss.cpp
@@ -5,6 +5,7 @@
 #include <cstdlib>
 #include <map>
 #include <vector>
+#include <algorithm>    // rotate
 
 #include "ham_ss.h"
 #include "seqTools.h"
@@ -201,6 +202,17 @@ void PottsHamiltonian::getSeq2State(const std::string &FILENAME)
        
     fclose(input);
 
+    
+    for (size_t i=0; i<Schar.size(); i++) {
+        state2aa.push_back(std::vector<aa>(Schar[i].size()));
+        std::vector<aa> last = state2aa.back();
+        for (size_t j=0; j<Schar[i].size(); j++) {
+            last.push_back(str2aa.at(Schar[i][j]));
+        }
+        // move the first element to last place, because that is how states are ordered
+        std::rotate(last.begin(),++last.begin(),last.end());    // could have used std::next() but eh
+    }
+
 
     //std::vector<std::array<int,aa_END> > aa2state;
     for (unsigned i=0; i<Schar.size(); i++) {

diff --git a/ham_ss.h b/ham_ss.h
index 76d0c62..c489016 100644 (file)
--- a/ham_ss.h
+++ b/ham_ss.h
@@ -82,6 +82,7 @@ public:
     double bJ;      // couplings "inverse temperature" multiplier
     Coupling J;
     std::vector<std::array<int,aa_END> > aa2state;
+    std::vector<std::vector<aa> > state2aa;
     
     PottsHamiltonian(const std::string &couplingsfile, const std::string &seq2stateFile);
     PottsHamiltonian() : bh(1.0), bJ(1.0) { }

diff --git a/ss.cpp b/ss.cpp
index 6f73581..b3ea4ad 100644 (file)
--- a/ss.cpp
+++ b/ss.cpp
@@ -224,6 +224,7 @@ void print_species_traj(double t, const Species_parray &print_spec, bool verbose
                 for (; it != end; ++it) {
                     //std::cout << ' ' << it->second
                     //          << '(' << aaseq2str(it->first.aa_seq) << ')';
+                    // unfortunately, the virus does not know about the Hamiltonian's state2aa, so we can't print the "internal" representation in a useful manner. TODO....
                     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