Fixed and added example for SARE python bindings

Author: GiovaGa

pyalp/pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "alp-graphblas"
-version = "0.8.41"
+version = "0.8.42"
 description = "Python bindings for ALP GraphBLAS (minimal package layout)"
 authors = [ { name = "ALP" } ]
 readme = "README.md"

pyalp/src/matrix_wrappers.hpp

@@ -41,7 +41,7 @@ void buildMatrix(
 
     grb::RC io_rc;
     (void)io_rc;
-    io_rc = grb::buildMatrixUnique( M, data_ptr_i, data_ptr_j , data_ptr_v, nnz, grb::SEQUENTIAL );
+    io_rc = grb::buildMatrixUnique( M, data_ptr_i, data_ptr_i + nnz, data_ptr_j, data_ptr_j + nnz, data_ptr_v, data_ptr_v + nnz, grb::SEQUENTIAL );
     assert( io_rc == grb::SUCCESS );
 }
 
@@ -57,7 +57,7 @@ grb::Matrix<ScalarType> matrix_factory(
 
     // Helper for dispatch
     bool handled = false;
-    auto try_type = [&](auto dummy) {
+    auto try_type = [&](const auto dummy) {
         using IntType = decltype(dummy);
         if (py::dtype::of<IntType>().is(data1.dtype()) && py::dtype::of<IntType>().is(data2.dtype())) {
             buildMatrix<IntType, ScalarType>(

pyalp/src/pyalp/common_bindings.hpp

@@ -3,6 +3,7 @@
 
 #include <pybind11/pybind11.h>
 #include <pybind11/numpy.h>
+#include <pybind11/stl.h>
 
 #include <graphblas.hpp>
 
@@ -14,6 +15,7 @@
 
 namespace py = pybind11;
 
+
 // Register all pyalp bindings. Module-local registration can be enabled by
 // instantiating with ModuleLocal = true. When ModuleLocal==true the
 // py::module_local() policy is applied to class bindings to avoid symbol
@@ -24,7 +26,7 @@ void register_pyalp(py::module_ &m) {
     // Common bindings for all backends
     m.def("backend_name", [](){ return "backend"; });
 
-    if constexpr (ModuleLocal) {
+    if (ModuleLocal) {
         py::class_<grb::Matrix< ScalarType >>(m, "Matrix", py::module_local())
         .def(py::init([](size_t m_, size_t n_,
                 py::array data1,
@@ -51,7 +53,20 @@ void register_pyalp(py::module_ &m) {
          py::arg("m"),
          py::arg("k_array")
         )
-        .def("to_numpy", &to_numpy, "Convert to numpy array");
+        .def("to_numpy", &to_numpy< ScalarType >, "Convert to numpy array");
+
+        py::class_<grb::Vector< StateType >>(m, "State", py::module_local())
+        .def(py::init<size_t>())
+        .def(py::init([](size_t m,
+                             py::array_t<StateType> data3) {
+                grb::Vector< StateType > vec(m); // call the basic constructor
+                buildVector(vec, data3); // initialize with data
+                return vec;
+            }),
+         py::arg("m"),
+         py::arg("k_array")
+        )
+        .def("to_numpy", &to_numpy< StateType >, "Convert to numpy array");
     } else {
         py::class_<grb::Matrix< ScalarType >>(m, "Matrix")
         .def(py::init([](size_t m_, size_t n_,
@@ -77,10 +92,54 @@ void register_pyalp(py::module_ &m) {
          py::arg("m"),
          py::arg("k_array")
         )
-        .def("to_numpy", &to_numpy, "Convert to numpy array");
+        .def("to_numpy", &to_numpy< ScalarType >, "Convert to numpy array");
+
+        py::class_<grb::Vector< StateType >>(m, "State")
+        .def(py::init<size_t>())
+        .def(py::init([](size_t m,
+                             py::array_t<StateType> data3) {
+                grb::Vector< StateType > vec(m); // call the basic constructor
+                buildVector(vec, data3); // initialize with data
+                return vec;
+            }),
+         py::arg("m"),
+         py::arg("k_array")
+        )
+        .def("to_numpy", &to_numpy< StateType >, "Convert to numpy array");
+
     }
+	
+	py::class_< std::vector<grb::Vector<StateType>> >(m, "stdVectorStates")
+        .def(py::init<size_t>())
+		.def(py::init([](
+						size_t m,
+						py::array_t< StateType > arr ) {
+				(void) m;
 
-    m.def("buildVector", &buildVector, "Fill Vector from 1 NumPy array");
+				const py::buffer_info buf = arr.request();
+				
+				if (buf.ndim != 2) {
+					throw std::runtime_error("Input array must be 2-dimensional");
+				}
+				const size_t sz = buf.shape[0];
+                std::vector< grb::Vector<StateType> > vec (sz);
+				assert( static_cast<size_t>( buf.shape[1] ) <= m );
+				auto ptr = static_cast<StateType*>(buf.ptr);
+
+    			grb::RC io_rc = grb::SUCCESS;
+				for (size_t i = 0; i < sz; i++) {
+					io_rc = io_rc ? io_rc :
+						grb::buildVector( vec[i], ptr, ptr + buf.shape[1], grb::SEQUENTIAL );
+				}
+                return vec;
+            }),
+         py::arg("m"),
+         py::arg("k_array")
+		 );
+        // .def("get_vec", & );
+
+    m.def("buildVector", &buildVector<ScalarType>, "Fill Vector from 1 NumPy array");
+    m.def("buildVectorInt8", &buildVector<StateType>, "Fill Vector from 1 NumPy array");
     m.def("print_my_numpy_array", &print_my_numpy_array, "Print a numpy array as a flattened std::vector");
     m.def("conjugate_gradient", &conjugate_gradient, "Pass alp data to alp CG solver",
       py::arg("L"),
@@ -114,4 +173,5 @@ void register_pyalp(py::module_ &m) {
       py::arg("seed") = 0,
       py::arg("verbose") = 0
     );
+
 }

pyalp/src/simulated_annealing_replica_exchange.hpp

@@ -33,8 +33,9 @@ ScalarType SARE_QUBO(
     ) {
     ScalarType best_energy = std::numeric_limits<ScalarType>::max();
 
-	if( !verbose )
+	if( !verbose ){
 	    std::cout << "SARE_QUBO:  start \n";
+	}
 	// get user process ID
 	const size_t s = grb::spmd<>::pid();
 	(void)s;
@@ -75,10 +76,10 @@ ScalarType SARE_QUBO(
 
 	if( !verbose ) {
 	    // output
-	    // std::cout << " solver_iterations = " << solver_iterations << "\n";
-	    // std::cout << " tol = " << tol << "\n";
-	    // std::cout << " iterations = " << iterations << "\n";
-	    // std::cout << " residual = " << residual << "\n";
+	    std::cout << " seed = " << seed << "\n";
+	    std::cout << " #iterations = " << solver_iterations << "\n";
+	    std::cout << " best_energy = " << best_energy << "\n";
+	    std::cout << " solver time = " << single_time << " s \n";
 	}
 
 	if( !verbose ) {
@@ -99,7 +100,7 @@ ScalarType SARE_Ising(
     const size_t seed = 0,
     size_t verbose = 0
     ) {
-    ScalarType best_energy = std::numeric_limits<ScalarType>::max();
+    ScalarType best_energy = 0;
 
 	if( !verbose )
 	    std::cout << "SARE_Ising:  start \n";
@@ -144,10 +145,10 @@ ScalarType SARE_Ising(
 
 	if( !verbose ) {
 	    // output
-	    // std::cout << " solver_iterations = " << solver_iterations << "\n";
-	    // std::cout << " tol = " << tol << "\n";
-	    // std::cout << " iterations = " << iterations << "\n";
-	    // std::cout << " residual = " << residual << "\n";
+	    std::cout << " seed = " << seed << "\n";
+	    std::cout << " #iterations = " << solver_iterations << "\n";
+	    std::cout << " best_energy = " << best_energy << "\n";
+	    std::cout << " solver time = " << single_time << " s \n";
 	}
 
 	if( !verbose ) {

pyalp/src/vector_wrappers.hpp

@@ -20,46 +20,39 @@
 namespace py = pybind11;
 
 
-using BaseScalarType = double;
-#ifdef _CG_COMPLEX
- using ScalarType = std::complex< BaseScalarType >;
-#else
- using ScalarType = BaseScalarType;
-#endif
-
-void buildVector(grb::Vector< ScalarType >& V, py::array_t<ScalarType> arrv) {
+template< typename T >
+void buildVector(grb::Vector< T >& V, py::array_t<T> arrv) {
 
     // Check array is 1D
     py::buffer_info info_v = arrv.request();
     if (info_v.ndim != 1) throw std::runtime_error("Array must be 1D");
-    ScalarType* data_ptr_v = static_cast<ScalarType*>(info_v.ptr);
+    T* data_ptr_v = static_cast<T*>(info_v.ptr);
 
     grb::RC io_rc;
     (void)io_rc;
     io_rc = grb::buildVector( V, data_ptr_v, data_ptr_v + info_v.size, grb::SEQUENTIAL );
     assert( io_rc == grb::SUCCESS );
 }
 
-py::array_t<ScalarType>
-to_numpy(grb::Vector< ScalarType >& x) {
-    grb::PinnedVector< ScalarType > pinnedVector;
-    pinnedVector = grb::PinnedVector< ScalarType >( x, grb::SEQUENTIAL );
+template< typename T >
+py::array_t< T > to_numpy( const grb::Vector< T >& x ) {
+    grb::PinnedVector< T > pinnedVector;
+    pinnedVector = grb::PinnedVector< T >( x, grb::SEQUENTIAL );
+	const size_t sz = pinnedVector.size();
 
     std::cout << "create numpy array from grb::vector\n";
 
-    ScalarType* data = new ScalarType[grb::size(x)];
-    for( size_t k = 0; k < grb::size(x); ++k ) {
-	const auto &value = pinnedVector.getNonzeroValue( k );
-	data[k]=value;
-    }
-
-    // Capsule to manage memory (will delete[] when array is destroyed in Python)
-    py::capsule free_when_done(data, [](void *f) {
-	delete[] reinterpret_cast<ScalarType*>(f);
-    });
+	auto result = py::array_t<T>(sz);
+	py::buffer_info buf = result.request();
+	T* ptr = static_cast< T* >( buf.ptr );
 
-    // Create NumPy array that shares memory with C++
-    py::array_t<ScalarType> arr({grb::size(x)}, {sizeof(ScalarType)}, data, free_when_done);
-    return arr;
+    for( size_t k = 0; k < sz; ++k ) {
+		ptr[ k ] = 0;
+	}
+    for( size_t k = 0; k < pinnedVector.nonzeroes(); ++k ) {
+		const auto &value = pinnedVector.getNonzeroValue( k );
+		ptr[pinnedVector.getNonzeroIndex( k )] = value;
+    }
 
+    return result;
 }

tests/python/sare_test.py

@@ -0,0 +1,110 @@
+"""
+Test script for the pyalp_ref (GraphBLAS-like) Python module, specifically the simulated annealing Replica-Exchange (SARE) stuff.
+
+This file should also work as a simple example of usage of such features.
+
+Usage:
+    python sare_test.py
+
+Dependencies:
+    - numpy
+    - scipy.sparse
+    - pyalp_ref (should be available in the Python path)
+"""
+
+import numpy as np
+from scipy.sparse import dok_array
+from time import time
+
+import pyalp
+# select backend
+backendname = "pyalp_omp" # or "pyalp_omp" "pyalp_nonblocking"
+pyalp = pyalp.get_backend(backendname)
+
+# set here your parameters
+
+n_replicas = 32
+seed = 0
+np.random.seed(seed)
+
+# This is the data for the problem matrix
+# note that data types have been fixed at compile time!
+N = 5
+J = dok_array((N, N), dtype=np.float64)
+
+J[0,1] = 2
+J[1,3] = 4
+J[2,3] = -1
+J[4,1] = 3
+J[4,1] = -2
+
+# You can also import a matrix as follows:
+from scipy.io import mmread
+# J = mmread(matpath).astype(np.float64)
+
+# read the local fields from a file
+# h = np.loadtxt(vecpath, dtype=np.float64)
+
+# or get the vector from the diagonal of the matrix
+# if np.any(J.diagonal() != 0):
+    # h = J.diagonal()
+    # J.setdiag(0)
+N = J.shape[0]
+
+# Make sure the matrix is symmetric with zero diagonal, as that is a precondition!
+J = (J + J.T)/2
+
+# if diagonal is not zero we can call
+assert np.all(J.diagonal() == 0)
+
+
+# create coordinates and value arrays
+J = J.todok()
+idata,jdata = np.array([[x,y] for x,y in J.keys()], dtype=np.int32).T
+vdata = np.array( list(J.values()), dtype=np.float64 )
+idata, jdata, vdata = np.sort(idata), jdata[np.argsort(idata)], vdata[np.argsort(idata)]
+
+# initial states
+states_numpy = np.random.randint(0,2, (n_replicas, N), dtype=np.int8 )
+
+# Initialize temperatures
+betas_numpy = np.logspace( 1e-2, 1e+2, n_replicas, dtype=np.float64 )
+
+# Initialize energies
+energies_numpy = np.diag(states_numpy@J@states_numpy.T)/2 + np.dot(states_numpy,h )
+
+niterations = 100
+verbose = 0
+
+#########################
+# Create the pyalp_ref Matrix and Vector objects
+J_alp = pyalp.Matrix( N, N, idata, jdata, vdata )
+h_alp = pyalp.Vector( N, h )
+states =  [pyalp.State(N, state) for state in states_numpy ]
+betas = pyalp.Vector( n_replicas, betas_numpy )
+# energies will be initialized by the library
+energies   = pyalp.Vector( n_replicas, energies_numpy  )
+
+# Preallocate best state
+best_state = pyalp.State( N, np.zeros(N, dtype=np.int8) )
+
+t0 = time()
+
+# For Ising solver
+best_energy = pyalp.SARE_Ising( J_alp, h_alp, states, energies, betas, best_state, niterations, seed, verbose )
+
+# There is also a QUBO solver, with the only difference being that h is fixed 0
+# best_energy = pyalp.SARE_Ising( J_alp, states, energies, betas, best_state, niterations, seed, verbose )
+t1 = time()
+
+best_state_numpy = best_state.to_numpy()
+
+print( best_energy )
+print( best_state_numpy )
+print( f"Solver time is: {(t1-t0)*1000} ms" )
+
+# check that energy is correct:
+assert best_state_numpy@(J@best_state_numpy/2 + h )  == best_energy
+
+
+