python c-extension for get_system_E().

get_system_E() is added to sqaod.py.
shinmorino · Nov 11, 2018 · 94206d2 · 94206d2
1 parent fde7bf8
commit 94206d2
Show file tree

Hide file tree

Showing 6 changed files with 84 additions and 0 deletions.
diff --git a/sqaodc/pyglue/annealer.inc b/sqaodc/pyglue/annealer.inc
@@ -820,6 +820,34 @@ PyObject *annealer_make_solution(PyObject *module, PyObject *args) {
     return Py_None;    
 }
 
+template<class real>
+PyObject *internal_get_system_E(PyObject *objExt, PyObject *objG, PyObject *objBeta) {
+    typedef NpConstScalarType<real> NpConstScalar;
+    NpConstScalar G(objG), beta(objBeta);
+    if (G.err || beta.err)
+        return NULL;
+    real E = pyobjToCppObj<real>(objExt)->getSystemE(G, beta);
+    PyObject *objE = newScalarObj(E);
+    // Py_INCREF(objE);
+    return objE;
+}
+
+extern "C"
+PyObject *annealer_get_system_E(PyObject *module, PyObject *args) {
+    PyObject *objExt, *objG, *objKT, *dtype;
+    if (!PyArg_ParseTuple(args, "OOOO", &objExt, &objG, &objKT, &dtype))
+        return NULL;
+
+    ASSERT_DTYPE(dtype);
+
+    TRY {
+        if (isFloat64(dtype))
+            return internal_get_system_E<double>(objExt, objG, objKT);
+        else // if (isFloat32(dtype))
+            return internal_get_system_E<float>(objExt, objG, objKT);
+    } CATCH_ERROR_AND_RETURN;
+    return NULL;
+}
 
 template<class real>
 PyObject *internal_anneal_one_step(PyObject *objExt, PyObject *objG, PyObject *objBeta) {
@@ -876,6 +904,7 @@ PyMethodDef annealer_methods[] = {
 	{"calculate_E", annealer_calculate_E, METH_VARARGS},
 	{"prepare", annealer_prepare, METH_VARARGS},
 	{"make_solution", annealer_make_solution, METH_VARARGS},
+	{"get_system_E", annealer_get_system_E, METH_VARARGS},
 	{"anneal_one_step", annealer_anneal_one_step, METH_VARARGS},
 	{NULL},
 };

diff --git a/sqaodpy/sqaod/common/bipartite_graph_annealer_base.py b/sqaodpy/sqaod/common/bipartite_graph_annealer_base.py
@@ -98,6 +98,10 @@ def prepare(self) :
 
     def make_solution(self) :
         self._cext.make_solution(self._cobj, self.dtype)
+
+    def get_system_E(self, G, beta) :
+        G, beta = self.dtype(G), self.dtype(beta)
+        return self._cext.get_system_E(self._cobj, G, beta, self.dtype)
 
     def anneal_one_step(self, G, beta) :
         G, beta = self.dtype(G), self.dtype(beta)

diff --git a/sqaodpy/sqaod/common/dense_graph_annealer_base.py b/sqaodpy/sqaod/common/dense_graph_annealer_base.py
@@ -92,6 +92,10 @@ def prepare(self) :
     def make_solution(self) :
         self._cext.make_solution(self._cobj, self.dtype)
 
+    def get_system_E(self, G, beta) :
+        G, beta = self.dtype(G), self.dtype(beta)
+        return self._cext.get_system_E(self._cobj, G, beta, self.dtype)
+
     def anneal_one_step(self, G, beta) :
         G, beta = self.dtype(G), self.dtype(beta)
         self._cext.anneal_one_step(self._cobj, G, beta, self.dtype)
diff --git a/sqaodpy/sqaod/common/preference.py b/sqaodpy/sqaod/common/preference.py
@@ -8,6 +8,16 @@ class Algorithm :
     sa_default = 'sa_default'
     sa_naive = 'sa_naive'
     sa_coloring = 'sa_coloring'
+
+    @staticmethod
+    def is_sqa(algo) :
+        if algo == Algorithm.default or \
+           algo == Algorithm.naive or \
+           algo == Algorithm.coloring :
+            return True
+
+        return False
+
 
 algorithm = Algorithm()
 """

diff --git a/sqaodpy/sqaod/py/bipartite_graph_annealer.py b/sqaodpy/sqaod/py/bipartite_graph_annealer.py
@@ -258,6 +258,26 @@ def make_solution(self) :
         """
         pass
 
+    def get_system_E(self, G, beta) :
+        # average energy
+        E = np.mean(self.get_E())
+
+        m = self._m
+        algo = self._get_algorithm()
+        if sqaod.algorithm.is_sqa(algo) :
+            q0, q1 = self._q0, self._q1
+            spinDotSum = 0.
+            for im in range(m) :
+                q00 = np.asarray(q0[im], np.float64)
+                q01 = np.asarray(q0[(im + 1) % m], np.float64)
+                spinDotSum += q00.dot(q01)
+                q10 = np.asarray(q1[im], np.float64)
+                q11 = np.asarray(q1[(im + 1) % m], np.float64)
+                spinDotSum += q10.dot(q11)
+            E += - 1. / (2 * beta) * np.log(np.tanh(G * beta / m)) * spinDotSum
+
+        return E
+
     def anneal_one_step(self, G, beta) :
         """ Run annealing one step.
 

diff --git a/sqaodpy/sqaod/py/dense_graph_annealer.py b/sqaodpy/sqaod/py/dense_graph_annealer.py
@@ -261,6 +261,23 @@ def make_solution(self) :
         """
         pass
 
+    def get_system_E(self, G, beta) :
+        # average energy
+        E = np.mean(self.get_E())
+
+        m = self._m
+        algo = self._get_algorithm()
+        if sqaod.algorithm.is_sqa(algo) :
+            q = self._q
+            spinDotSum = 0.
+            for im in range(m) :
+                q0 = np.asarray(q[im], np.float64)
+                q1 = np.asarray(q[(im + 1) % m], np.float64)
+                spinDotSum += q0.dot(q1)
+            E -= 0.5 / beta * np.log(np.tanh(G * beta / m)) * spinDotSum
+
+        return E
+
     def anneal_one_step(self, G, beta) :
         """ Run annealing one step.