From 475ff9161054e02f6696e65392bed00a5e44817a Mon Sep 17 00:00:00 2001
From: Dariusz Murakowski <murakdar@mit.edu>
Date: Tue, 24 Feb 2015 12:14:52 -0500
Subject: [PATCH] Code to read Potts coupling parameters and calculate energy,
 from John Barton 20150218.

---
 snip.cpp | 68 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 snip.h   | 48 +++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 116 insertions(+)
 create mode 100644 snip.cpp
 create mode 100644 snip.h

diff --git a/snip.cpp b/snip.cpp
new file mode 100644
index 0000000..dde3cfd
--- /dev/null
+++ b/snip.cpp
@@ -0,0 +1,68 @@
+#include <vector>
+
+#include "snip.h"
+
+
+// Read correlations and number of sites in from a file
+
+void getCorrelations(FILE *input, Vector &p) {
+
+    double c;
+    char o;
+    
+    while (fscanf(input,"%le",&c)==1) {
+    
+        p.push_back(std::vector<double>());
+        (p.back()).push_back(c);
+        
+        while (fscanf(input,"%c",&o)==1) {
+    
+            if (o=='\n' || o=='\r') break;
+            
+            fscanf(input,"%le",&c);
+            (p.back()).push_back(c);
+            
+        }
+        
+    }
+	
+}
+
+
+// Read couplings in from a file
+
+void getCouplings(FILE *input, Vector &p) {
+
+    getCorrelations(input, p);
+	
+}
+
+
+// Set energies based on an initial configuration
+
+double computeEnergy(const Vector &J, const std::vector<int> &config) {
+    
+    // Compute energy of the input configuration
+    
+    double E = 0;
+    
+    for (int i=0;i<config.size();i++) {
+        
+        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/snip.h b/snip.h
new file mode 100644
index 0000000..cd1cd1f
--- /dev/null
+++ b/snip.h
@@ -0,0 +1,48 @@
+#ifndef SNIP_H
+#define SNIP_H
+
+
+#include <vector>
+#include <math.h>
+#include <stdio.h>
+
+
+// Typedefs
+typedef std::vector<unsigned long> Key;
+typedef std::vector<std::vector<double> > Vector;
+typedef std::vector<std::vector<int> > IntVector;
+
+
+void getCouplings(FILE *, Vector &);
+void getCorrelations(FILE *, Vector &);
+double computeEnergy(const Vector &, const std::vector<int> &);
+
+
+// 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);
+    
+}
+
+
+#endif
-- 
2.7.4