From: Dariusz Murakowski <murakdar@mit.edu>
Date: Thu, 23 Apr 2015 05:04:26 +0000 (-0400)
Subject: Copy John's Potts code from snip.{cpp,h}.
X-Git-Url: http://src.murakowski.org/?a=commitdiff_plain;h=6eb9e4dbb09be87a971f95e08254fffbf55df62d;p=VirEvoDyn.git

Copy John's Potts code from snip.{cpp,h}.
---

diff --git a/ham_ss.cpp b/ham_ss.cpp
index 817dca2..323617d 100644
--- a/ham_ss.cpp
+++ b/ham_ss.cpp
@@ -103,7 +103,6 @@ Hamiltonian::Hamiltonian(std::string &FILENAME) {
 }
 
 
-
 // Return the energy for a given set of mutated sites
 
 double Hamiltonian::get_energy(const std::set<unsigned int> &mutated_sites) const {
@@ -128,3 +127,69 @@ double Hamiltonian::get_energy(const std::set<unsigned int> &mutated_sites) cons
     
 }
 
+
+
+// Read correlations and number of sites in from a file
+
+PottsHamiltonian::PottsHamiltonian(std::string &FILENAME) {
+
+    FILE *input = fopen(FILENAME.c_str(),"r");    // <paramfile>.j
+    if (input == NULL) { perror((std::string("ERROR: ") + FILENAME).c_str()); exit(1); }
+
+    double c;
+    char o;
+    
+    while (fscanf(input,"%le",&c)==1) {
+    
+        J.push_back(std::vector<double>());
+        (J.back()).push_back(c);
+        
+        while (fscanf(input,"%c",&o)==1) {
+    
+            if (o=='\n' || o=='\r') break;
+            
+            fscanf(input,"%le",&c);
+            (J.back()).push_back(c);
+            
+        }
+        
+    }
+	
+}
+
+
+// Set energies based on an initial configuration
+// 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
+// the commonest states to zero.
+
+double PottsHamiltonian::get_energy(const std::vector<int> &config) const {
+    
+    // Compute energy of the input configuration
+    
+    double E = 0;
+    
+    for (int i=0;i<config.size();i++) {
+
+        //assert(config[i] <= J[i].size());
+        
+        if (config[i]!=J[i].size()) {
+        
+            E -= J[i][config[i]];
+            
+            for (int 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())];
+                
+            }
+            
+        }
+        
+    }
+    
+    return E;
+    
+}
+
+
diff --git a/ham_ss.h b/ham_ss.h
index 2c87010..dc949b2 100644
--- a/ham_ss.h
+++ b/ham_ss.h
@@ -4,6 +4,7 @@
 #include <string>
 #include <vector>
 #include <set>
+#include <math.h>   // sqrt()
 
 #include "pop_ss.h"
 
@@ -41,4 +42,49 @@ public:
 
 };
 
+
+// Given the size of a coupling or correlation vector, return the corresponding number of spins
+
+inline int sizetolength(size_t size) {
+    
+    return (int) ((sqrt(1 + 8 * size) - 1) / 2);
+    
+}
+
+
+// Return the location of a pair of states in canonical order, {{0,0}, {0,1}, ..., {0,qj}, {1,0}, ...}
+
+inline int sindex(int i, int j, size_t qi, size_t qj) {
+    
+    return (int) (i * qj + j);
+    
+}
+
+
+// For j>i, the pair {i,j} in the list {{0},{1},...,{length-1},{0,1}, {0,2}, ... } is located at index(i,j,length)
+
+inline int index(int i, int j, size_t length) {
+    
+    return (int)(length + i * (length - 2) - (i * (i - 1)) / 2 - 1 + j);
+    
+}
+
+
+class PottsHamiltonian {
+
+public:
+    double bh;      // fields "inverse temperature" multiplier
+    double bJ;      // couplings "inverse temperature" multiplier
+    Coupling J;
+    
+    PottsHamiltonian(std::string &FILENAME);
+    PottsHamiltonian() : bh(1.0), bJ(1.0) { }
+    
+    double get_energy(const std::vector<int> &config) const;
+    
+    void set_temp(double x)           { bh=x; bJ=x; }
+    void set_temp(double x, double y) { bh=x; bJ=y; }
+
+};
+
 #endif  // HAM_SS_H