Check number of inserted sites. Use unsigned int instead of size_t.

diff --git a/ham_ss.cpp b/ham_ss.cpp
index 573b355..0fb6dad 100644 (file)
--- a/ham_ss.cpp
+++ b/ham_ss.cpp
@@ -290,28 +290,28 @@ 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()
+unsigned PottsHamiltonian::insertReferenceStates()
 {
     // find indices of singular states
-    std::vector<size_t> Xidx;
+    std::vector<unsigned> Xidx;
     //Xidx.push_back(0);      // prepend with 0, for starting interval
-    for (size_t i = 0; i < state2aa.size(); ++i) {
+    for (unsigned 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;
+    unsigned numX = Xidx.size(); // - 2;
+    unsigned sizeOld = size - numX;
     //assert(J.size() == (sizeOld+1)*(sizeOld)/2);
 
     // bail out early if no need to insert sites
-    //if (numX == 0)
-    //    return 0;
+    if (numX == 0)
+        return 0;
 
     Coupling newJ;
-    newJ.reserve(J.size() + numX*size - numX*(numX-1)/2);
-    //newJ.reserve(size*(size+1)/2);
+    //newJ.reserve(J.size() + numX*size - numX*(numX-1)/2);
+    newJ.reserve(size*(size+1)/2);
 
     Coupling::iterator Jbegin = J.begin();
     //Coupling::iterator end = J.begin();
@@ -323,6 +323,8 @@ int PottsHamiltonian::insertReferenceStates()
     //    }
     //}
 
+    unsigned inserted = 0;
+
     unsigned srcBeg = 0;
     unsigned srcEnd = 0;
     unsigned Xdone = 0;
@@ -331,7 +333,7 @@ int PottsHamiltonian::insertReferenceStates()
     //unsigned dstBeg = 0;
     //unsigned dstEnd = 0;
     // iterate over regions in the destination vector
-    for (size_t sec = 0; sec < size; ++sec) {
+    for (unsigned sec = 0; sec < size; ++sec) {
         unsigned secLen = size - sec;   // dst
         //dstBeg = index(sec,sec+1,size) - size + sec;
         //equivalently, srcBeg = srcEnd; ... but this may prevent loop unrolling?
@@ -341,12 +343,13 @@ int PottsHamiltonian::insertReferenceStates()
         sectionBeg = srcBeg;
         copyToEnd = true;
         //dstBeg = index(sec,sec+1,size) - secLen;
-        for (size_t i = Xdone; i < numX; ++i) {
+        for (unsigned i = Xdone; i < numX; ++i) {
             if (sec <= Xidx[i]) {
                 srcEnd = sectionBeg + Xidx[i] + Xdone - i - sec;
                 // copy
                 newJ.insert(newJ.end(), Jbegin+srcBeg, Jbegin+srcEnd);
                 newJ.push_back(std::vector<double>());
+                inserted += 1;
                 // copy from begin to Xidx[i]-sec
                 // copy Xidx[i-1],Xidx[i] + i-1;
                 srcBeg = srcEnd;
@@ -354,6 +357,7 @@ int PottsHamiltonian::insertReferenceStates()
             else if (sec == Xidx[i]+1) {
                 // insert a whole bunch of nothings
                 newJ.insert(newJ.end(), secLen, std::vector<double>());
+                inserted += secLen;
                 ++Xdone;    // skip this index in subsequent blocks
                 copyToEnd = false;
                 break;
@@ -373,7 +377,10 @@ int PottsHamiltonian::insertReferenceStates()
         }
     }
 
-    int inserted = 0;
+    if (inserted != numX*size - numX*(numX-1)/2) {
+        std::cerr << "ERROR in insertReferenceStates: incorrectly calculated padding" << std::endl;
+        throw "incorrectly calculated correction....";
+    }
     //for (size_t idx = 0; idx < Xidx.size(); ++idx) {
     //    i = Xidx[idx];
     //    end += i;

diff --git a/ham_ss.h b/ham_ss.h
index 78a4208..3a8cdec 100644 (file)
--- a/ham_ss.h
+++ b/ham_ss.h
@@ -97,7 +97,7 @@ public:
     void getCouplings(const std::string &FILENAME);
     void getSeq2State(const std::string &FILENAME);
 
-    int insertReferenceStates();    // XXX
+    unsigned insertReferenceStates();    // XXX
 
 };