From 9f723425a1b82c6dd6efc9dc9a58c15f6a6e2d46 Mon Sep 17 00:00:00 2001
From: Dariusz Murakowski <murakdar@mit.edu>
Date: Tue, 23 Jun 2015 04:49:19 -0400
Subject: [PATCH] WIP attempt to insert dummy sites into J for
 degenerate/singular Potts states representing all amino acids ("X").

---
 ham_ss.cpp | 67 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 ham_ss.h   |  2 ++
 2 files changed, 69 insertions(+)

diff --git a/ham_ss.cpp b/ham_ss.cpp
index 61d64c0..5bb951a 100644
--- a/ham_ss.cpp
+++ b/ham_ss.cpp
@@ -135,6 +135,7 @@ PottsHamiltonian::PottsHamiltonian(const std::string &couplingsFile, const std::
 {
     getCouplings(couplingsFile);
     getSeq2State(seq2stateFile);
+    //insertReferenceStates();  // XXX
     this->size = aa2state.size();
 }
 
@@ -287,4 +288,70 @@ void PottsHamiltonian::coarse_grain(const std::vector<aa> &aa_seq, PottsState &c
     }
 }
 
+// some sites are 
+// Returns the number of elements inserted into J.
+int PottsHamiltonian::insertReferenceStates()
+{
+    // find indices of singular states
+    std::vector<size_t> Xidx;
+    //Xidx.push_back(0);      // prepend with 0, for starting interval
+    for (size_t i = 0; i < state2aa.size(); ++i) {
+        if (state2aa[i].size() == 1) {  // state2aa[i] == {aa_X}
+            Xidx.push_back(i);
+        }
+    }
+    //Xidx.push_back(size);   // append with size, for ending interval
+    size_t numX = Xidx.size(); // - 2;
+    size_t sizeOld = size - numX;
+    //assert(J.size() == (sizeOld+1)*(sizeOld)/2);
+
+    // bail out early if no need to insert sites
+    //if (numX == 0)
+    //    return 0;
+
+    Coupling newJ;
+
+    Coupling::iterator begin = J.begin();
+    Coupling::iterator end = J.begin();
+
+    //// iterate over 
+    //for (size_t i = 0; i < numX; ++i) {
+    //    for (size_t j = i+1; j < numX; ++j) {
+
+    //    }
+    //}
+
+    // iterate over regions in the destination vector
+    for (size_t sec = 0; sec < size; ++sec) {
+        size_t secLen = size - sec;
+        for (size_t i = 0; i < numX; ++i) {
+            if (Xidx[i] < sec) {
+                // copy from begin to Xidx[i]-sec
+                // copy Xidx[i-1],Xidx[i] + i-1;
+            }
+            else if (Xidx[i] == sec+1) {
+                // insert a whole bunch of nothings
+            }
+            begin = end;
+        }
+        end = J.begin() + (size-sec)*(size-sec-1)/2;  // or something
+        // copy from begin to end of section
+    }
+
+    int inserted = 0;
+    //for (size_t idx = 0; idx < Xidx.size(); ++idx) {
+    //    i = Xidx[idx];
+    //    end += i;
+
+    //        end += i;
+    //        newJ.insert(newJ.end(),begin,end);  // copy range [begin,end)
+    //        newJ.push_back(std::vector<double>());  // add dummy entry
+    //        begin += i;     // next copy starts here
+    //        inserted++;
+
+    J.assign(newJ.begin(),newJ.end());
+
+    return inserted;
+}
+
 
diff --git a/ham_ss.h b/ham_ss.h
index c489016..78a4208 100644
--- a/ham_ss.h
+++ b/ham_ss.h
@@ -97,6 +97,8 @@ public:
     void getCouplings(const std::string &FILENAME);
     void getSeq2State(const std::string &FILENAME);
 
+    int insertReferenceStates();    // XXX
+
 };
 
 #endif  // HAM_SS_H
-- 
2.7.4