From 73c4d9f2e026383cd9a1ed36a2708cb535df7985 Mon Sep 17 00:00:00 2001
From: Dariusz Murakowski <murakdar@mit.edu>
Date: Fri, 24 Apr 2015 02:34:30 -0400
Subject: [PATCH] Potts state mapping+tracking complete. Fix bounds check error
 in energy calculation.

---
 ham_ss.cpp | 11 ++++++++++-
 ham_ss.h   |  4 +++-
 pop_ss.cpp | 20 ++++++++++++++++++++
 pop_ss.h   | 10 +++++++---
 ss.cpp     | 11 ++++++++---
 5 files changed, 48 insertions(+), 8 deletions(-)

diff --git a/ham_ss.cpp b/ham_ss.cpp
index 78a7834..94232d7 100644
--- a/ham_ss.cpp
+++ b/ham_ss.cpp
@@ -251,7 +251,7 @@ double PottsHamiltonian::get_energy(const PottsState &config) const
             
             for (unsigned j=i+1;j<config.size();j++) {
             
-                if (config[j]!=J[i].size())
+                if (config[j]!=J[j].size())
                     Ecoupling -= J[index(i,j,config.size())][sindex(config[i],config[j],J[i].size(),J[j].size())];
                 
             }
@@ -265,4 +265,13 @@ double PottsHamiltonian::get_energy(const PottsState &config) const
     
 }
 
+// put the Potts state representation of the amino acid sequence into the given 'config' object
+void PottsHamiltonian::coarse_grain(const std::vector<aa> &aa_seq, PottsState &config) const
+{
+    //assert(this->size == config.size() && this->size == aa_seq.size());
+    for (size_t i = 0; i < this->size; ++i) {
+        config[i] = this->aa2state[i][aa_seq[i]];
+    }
+}
+
 
diff --git a/ham_ss.h b/ham_ss.h
index bbfbbdd..ea0a946 100644
--- a/ham_ss.h
+++ b/ham_ss.h
@@ -16,7 +16,7 @@ typedef std::vector<std::vector<double> > Coupling;
 
 typedef std::set<unsigned int> MutatedSiteSequence;
 
-typedef std::vector<aa> PottsState;
+typedef std::vector<unsigned> PottsState;
 
 
 class Virus;
@@ -87,6 +87,8 @@ public:
     PottsHamiltonian() : bh(1.0), bJ(1.0) { }
     
     double get_energy(const PottsState &config) const;
+
+    void coarse_grain(const std::vector<aa> &aa_seq, PottsState &config) const;
     
     void set_temp(double x)           { bh=x; bJ=x; }
     void set_temp(double x, double y) { bh=x; bJ=y; }
diff --git a/pop_ss.cpp b/pop_ss.cpp
index fb4851b..ba8ddad 100644
--- a/pop_ss.cpp
+++ b/pop_ss.cpp
@@ -62,9 +62,29 @@ NTVirus::NTVirus(const std::string &NT)
 , aa_seq(translate_nt_seq(NT))
 , L_nt(NT.length())
 , L_aa(this->L_nt/3)
+//, config()  // default constructor? warning: order disregarded
 {
 }
 
+NTVirus::NTVirus(const PottsHamiltonian &H, const std::string &NT)
+: nt_seq(NT)
+, aa_seq(translate_nt_seq(NT))
+, L_nt(NT.length())
+, L_aa(this->L_nt/3)
+{
+    //assert(this->L_aa == H.size);
+    if (this->L_aa != H.size)
+        throw "mismatch between Hamiltonian size and sequence length.";
+    config.resize(L_aa);
+    H.coarse_grain(aa_seq,config);
+    update_energy(H);
+}
+
+void NTVirus::update_energy(const PottsHamiltonian &H)
+{
+    this->energy = H.get_energy(config);
+}
+
 bool operator<(const NTVirus& lhs, const NTVirus& rhs) {
     return lhs.nt_seq < rhs.nt_seq;
 }
diff --git a/pop_ss.h b/pop_ss.h
index 083d656..ea92489 100644
--- a/pop_ss.h
+++ b/pop_ss.h
@@ -11,9 +11,10 @@
 #include "seqTools.h"
 
 
-typedef std::vector<aa> PottsState;
+typedef std::vector<unsigned> PottsState;
 
 class Hamiltonian;
+class PottsHamiltonian;
 
 
 class Virus {
@@ -51,11 +52,14 @@ class NTVirus {
 public:
     double energy;
     std::string nt_seq;
-    PottsState aa_seq;
+    std::vector<aa> aa_seq;
+    PottsState config;
     unsigned int L_nt;
     unsigned int L_aa;
 
-    NTVirus(const std::string &);
+    NTVirus(const std::string &NT);
+    NTVirus(const PottsHamiltonian &H, const std::string &NT);
+    void update_energy(const PottsHamiltonian &H);
 };
 
 // needed for STL map
diff --git a/ss.cpp b/ss.cpp
index 6523be8..3e84214 100644
--- a/ss.cpp
+++ b/ss.cpp
@@ -501,16 +501,21 @@ int main(int argc, char *argv[]) {
     }
     for (unsigned i=0; i<H.aa2state.size(); i++) {
         for (int aa=aa_BEGIN; aa != aa_END; ++aa) {
-            std::cout << H.aa2state[i][aa] << ' ';
+            printf("%c(%d) ", int2aa.at(aa),H.aa2state[i][aa]);
+            //std::cout << H.aa2state[i][aa] << ' ';
         }
         std::cout << '\n';
     }
     //NTVirus v("TTT");
-    NTVirus v(argv[3]);
+    NTVirus v(H,argv[3]);
     //std::cout << int2aa.at(codon2aa.at(str2codon.at(v.nt_seq))) << '\n';
     std::cout << aa2str.at(v.aa_seq[0]) << '\n';
     std::cout << aaseq2str(v.aa_seq) << '\n';
-    printf("%d %d\n",v.L_nt,v.L_aa);
+    for (size_t i=0; i<v.L_aa; i++) {
+        std::cout << v.config[i] << ' ';
+    }
+    std::cout << '\n';
+    printf("Nnt=%d Naa=%d E=%f\n",v.L_nt,v.L_aa,v.energy);
     return 0;
 //*/
 
-- 
2.7.4