Code to read Potts coupling parameters and calculate energy, from John Barton 20150218.

diff --git a/snip.cpp b/snip.cpp
new file mode 100644 (file)
index 0000000..dde3cfd
--- /dev/null
+++ b/snip.cpp
@@ -0,0 +1,68 @@
+#include <vector>\r
+\r
+#include "snip.h"\r
+\r
+\r
+// Read correlations and number of sites in from a file\r
+\r
+void getCorrelations(FILE *input, Vector &p) {\r
+\r
+    double c;\r
+    char o;\r
+    \r
+    while (fscanf(input,"%le",&c)==1) {\r
+    \r
+        p.push_back(std::vector<double>());\r
+        (p.back()).push_back(c);\r
+        \r
+        while (fscanf(input,"%c",&o)==1) {\r
+    \r
+            if (o=='\n' || o=='\r') break;\r
+            \r
+            fscanf(input,"%le",&c);\r
+            (p.back()).push_back(c);\r
+            \r
+        }\r
+        \r
+    }\r
+       \r
+}\r
+\r
+\r
+// Read couplings in from a file\r
+\r
+void getCouplings(FILE *input, Vector &p) {\r
+\r
+    getCorrelations(input, p);\r
+       \r
+}\r
+\r
+\r
+// Set energies based on an initial configuration\r
+\r
+double computeEnergy(const Vector &J, const std::vector<int> &config) {\r
+    \r
+    // Compute energy of the input configuration\r
+    \r
+    double E = 0;\r
+    \r
+    for (int i=0;i<config.size();i++) {\r
+        \r
+        if (config[i]!=J[i].size()) {\r
+        \r
+            E -= J[i][config[i]];\r
+            \r
+            for (int j=i+1;j<config.size();j++) {\r
+            \r
+                if (config[j]!=J[i].size()) E -= J[index(i,j,config.size())][sindex(config[i],config[j],J[i].size(),J[j].size())];\r
+                \r
+            }\r
+            \r
+        }\r
+        \r
+    }\r
+    \r
+    return E;\r
+    \r
+}\r
+\r

diff --git a/snip.h b/snip.h
new file mode 100644 (file)
index 0000000..cd1cd1f
--- /dev/null
+++ b/snip.h
@@ -0,0 +1,48 @@
+#ifndef SNIP_H\r
+#define SNIP_H\r
+\r
+\r
+#include <vector>\r
+#include <math.h>\r
+#include <stdio.h>\r
+\r
+\r
+// Typedefs\r
+typedef std::vector<unsigned long> Key;\r
+typedef std::vector<std::vector<double> > Vector;\r
+typedef std::vector<std::vector<int> > IntVector;\r
+\r
+\r
+void getCouplings(FILE *, Vector &);\r
+void getCorrelations(FILE *, Vector &);\r
+double computeEnergy(const Vector &, const std::vector<int> &);\r
+\r
+\r
+// Given the size of a coupling or correlation vector, return the corresponding number of spins\r
+\r
+inline int sizetolength(size_t size) {\r
+    \r
+    return (int) ((sqrt(1 + 8 * size) - 1) / 2);\r
+    \r
+}\r
+\r
+\r
+// Return the location of a pair of states in canonical order, {{0,0}, {0,1}, ..., {0,qj}, {1,0}, ...}\r
+\r
+inline int sindex(int i, int j, size_t qi, size_t qj) {\r
+    \r
+    return (int) (i * qj + j);\r
+    \r
+}\r
+\r
+\r
+// 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)\r
+\r
+inline int index(int i, int j, size_t length) {\r
+    \r
+    return (int)(length + i * (length - 2) - (i * (i - 1)) / 2 - 1 + j);\r
+    \r
+}\r
+\r
+\r
+#endif\r