From: Dariusz Murakowski <murakdar@mit.edu>
Date: Thu, 23 Apr 2015 06:18:24 +0000 (-0400)
Subject: Break normal functionality for testing! Separate beta (inverse temperature).
X-Git-Url: http://src.murakowski.org/?a=commitdiff_plain;h=5baf59a90fb800b69f104499a9cd12ea0577ada9;p=VirEvoDyn.git

Break normal functionality for testing! Separate beta (inverse temperature).
---

diff --git a/ham_ss.cpp b/ham_ss.cpp
index 9eba31a..e988411 100644
--- a/ham_ss.cpp
+++ b/ham_ss.cpp
@@ -129,9 +129,18 @@ double Hamiltonian::get_energy(const std::set<unsigned int> &mutated_sites) cons
 
 
 
+PottsHamiltonian::PottsHamiltonian(const std::string &couplingsFile, const std::string &seq2stateFile) : bh(1.0), bJ(1.0)
+{
+    getCouplings(couplingsFile);
+    getSeq2State(seq2stateFile);
+}
+
+
 // Read correlations and number of sites in from a file
+// Code adapted from John Barton.
 
-PottsHamiltonian::PottsHamiltonian(std::string &FILENAME) {
+void PottsHamiltonian::getCouplings(const std::string &FILENAME)
+{
 
     FILE *input = fopen(FILENAME.c_str(),"r");    // <paramfile>.j
     if (input == NULL) { perror((std::string("ERROR: ") + FILENAME).c_str()); exit(1); }
@@ -157,18 +166,25 @@ PottsHamiltonian::PottsHamiltonian(std::string &FILENAME) {
 	
 }
 
+void PottsHamiltonian::getSeq2State(const std::string &FILENAME)
+{
+}
+
 
 // Set energies based on an initial configuration
+// Code adapted from John Barton.
 // Note: config[i] = 0 denotes the 2nd most common state,
-// whereas config[i] = J[i].size() is the top state
-// Also, "gauge fixing" sets field and couplings between
+// whereas config[i] = J[i].size() is the top state.
+// Also, "gauge fixing" sets field at and couplings between
 // the commonest states to zero.
 
-double PottsHamiltonian::get_energy(const std::vector<unsigned> &config) const {
+double PottsHamiltonian::get_energy(const PottsState &config) const
+{
     
     // Compute energy of the input configuration
     
-    double E = 0;
+    double Efield = 0;
+    double Ecoupling = 0;
     
     for (unsigned i=0;i<config.size();i++) {
 
@@ -176,12 +192,12 @@ double PottsHamiltonian::get_energy(const std::vector<unsigned> &config) const {
         
         if (config[i]!=J[i].size()) {
         
-            E -= J[i][config[i]];
+            Efield -= J[i][config[i]];
             
             for (unsigned j=i+1;j<config.size();j++) {
             
                 if (config[j]!=J[i].size())
-                    E -= J[index(i,j,config.size())][sindex(config[i],config[j],J[i].size(),J[j].size())];
+                    Ecoupling -= J[index(i,j,config.size())][sindex(config[i],config[j],J[i].size(),J[j].size())];
                 
             }
             
@@ -189,6 +205,7 @@ double PottsHamiltonian::get_energy(const std::vector<unsigned> &config) const {
         
     }
     
+    double E = bh*Efield + bJ*Ecoupling;
     return E;
     
 }
diff --git a/ham_ss.h b/ham_ss.h
index d18306f..951bce4 100644
--- a/ham_ss.h
+++ b/ham_ss.h
@@ -14,6 +14,8 @@ typedef std::vector<std::vector<double> > Coupling;
 
 typedef std::set<unsigned int> MutatedSiteSequence;
 
+typedef std::vector<unsigned> PottsState;
+
 
 class Virus;
 
@@ -77,14 +79,17 @@ public:
     double bJ;      // couplings "inverse temperature" multiplier
     Coupling J;
     
-    PottsHamiltonian(std::string &FILENAME);
+    PottsHamiltonian(const std::string &couplingsfile, const std::string &seq2stateFile);
     PottsHamiltonian() : bh(1.0), bJ(1.0) { }
     
-    double get_energy(const std::vector<unsigned> &config) const;
+    double get_energy(const PottsState &config) const;
     
     void set_temp(double x)           { bh=x; bJ=x; }
     void set_temp(double x, double y) { bh=x; bJ=y; }
 
+    void getCouplings(const std::string &FILENAME);
+    void getSeq2State(const std::string &FILENAME);
+
 };
 
 #endif  // HAM_SS_H
diff --git a/ss.cpp b/ss.cpp
index c4fecad..8131cb7 100644
--- a/ss.cpp
+++ b/ss.cpp
@@ -488,6 +488,16 @@ void usage()
 
 int main(int argc, char *argv[]) {
 
+    std::string fname = argv[1];
+    PottsHamiltonian H(fname,"");
+    for (unsigned i=0; i<H.J.size(); i++) {
+        for (unsigned j=0; j<H.J[i].size(); j++) {
+            std::cout << H.J[i][j] << ' ';
+        }
+        std::cout << '\n';
+    }
+    return 0;
+
     RunParameters_SS r;
 
     unsigned seed = sim_random_seed();