Copy J matrix in blocks, padding with empty entries.

diff --git a/ham_ss.cpp b/ham_ss.cpp
index 5bb951a..b0793fb 100644 (file)
--- a/ham_ss.cpp
+++ b/ham_ss.cpp
@@ -135,8 +135,8 @@ PottsHamiltonian::PottsHamiltonian(const std::string &couplingsFile, const std::
 {
     getCouplings(couplingsFile);
     getSeq2State(seq2stateFile);
-    //insertReferenceStates();  // XXX
     this->size = aa2state.size();
+    insertReferenceStates();  // XXX
 }
 
 
@@ -310,9 +310,10 @@ int PottsHamiltonian::insertReferenceStates()
     //    return 0;
 
     Coupling newJ;
+    newJ.reserve(J.size() + numX*size);
 
-    Coupling::iterator begin = J.begin();
-    Coupling::iterator end = J.begin();
+    Coupling::iterator Jbegin = J.begin();
+    //Coupling::iterator end = J.begin();
 
     //// iterate over 
     //for (size_t i = 0; i < numX; ++i) {
@@ -321,21 +322,45 @@ int PottsHamiltonian::insertReferenceStates()
     //    }
     //}
 
+    unsigned srcBeg = 0;
+    unsigned srcEnd = 0;
+    unsigned Xdone = 0;
+    //unsigned dstBeg = 0;
+    //unsigned dstEnd = 0;
     // 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) {
+        unsigned secLen = size - sec;   // dst
+        //dstBeg = index(sec,sec+1,size) - size + sec;
+        //equivalently, srcBeg = srcEnd; ... but this may prevent loop unrolling?
+        srcBeg = index(sec,sec+1,sizeOld) - sizeOld + sec;
+        srcEnd = srcBeg;
+        //dstBeg = index(sec,sec+1,size) - secLen;
+        for (size_t i = Xdone; i < numX; ++i) {
+            if (sec <= Xidx[i]) {
+                srcEnd += Xidx[i] - i - sec;
+                // copy
+                newJ.insert(newJ.end(), Jbegin+srcBeg, Jbegin+srcEnd);
+                newJ.push_back(std::vector<double>());
                 // copy from begin to Xidx[i]-sec
                 // copy Xidx[i-1],Xidx[i] + i-1;
+                srcBeg = srcEnd;
             }
-            else if (Xidx[i] == sec+1) {
+            else if (sec == Xidx[i]+1) {
                 // insert a whole bunch of nothings
+                newJ.insert(newJ.end(), secLen, std::vector<double>());
+                ++Xdone;    // skip this index in subsequent blocks
+                //for (size_t k = 0; k < secLen; ++k)
+                //    newJ.push_back(std::vector<double>());
             }
-            begin = end;
+            //dstBeg = dstEnd;
+            //srcBeg = srcEnd;
+            // TODO: consider using a deque and Xidx.pop_front()
         }
-        end = J.begin() + (size-sec)*(size-sec-1)/2;  // or something
+        //end = J.begin() + (size-sec)*(size-sec-1)/2;  // or something
+        //dstEnd = index(sec,sec+1,size);
+        srcEnd = index(sec,sec+1,sizeOld);
         // copy from begin to end of section
+        newJ.insert(newJ.end(), Jbegin+srcBeg, Jbegin+srcEnd);
     }
 
     int inserted = 0;