Added example

src/testlib/Makefile

@@ -0,0 +1,2 @@
+all:
+	gfortran -llapack -lblas -shared ./funclib.f90 -o funclib.so

src/testlib/__init__.py

@@ -0,0 +1,81 @@
+import ctypes as ct
+import ctypes
+import numpy as np
+
+from ctypes import CDLL, POINTER, c_int, c_double, c_bool
+
+# import the shared library
+fortlib = ct.CDLL('./funclib.so') 
+
+c_float_p = ctypes.POINTER(ctypes.c_float)
+c_integer_p = ctypes.POINTER(ctypes.c_int)
+c_bool_p = ctypes.POINTER(ctypes.c_bool)
+
+# Specify arguments and result types
+fortlib.sum2.argtypes = [ct.POINTER(ct.c_double)]
+fortlib.sum2.restype = ct.c_double
+
+# fortlib.fkpca.argtypes = [ct.POINTER(ct.c_double)]
+# fortlib.fkpca.restype = ct.c_double
+
+
+def kpca(K, n=2, centering=True):
+    """Calculates `n` first principal components for the kernel :math:`K`.
+
+    The PCA is calculated using an OpenMP parallel Fortran routine.
+
+    A square, symmetric kernel matrix is required. Centering of the kernel matrix
+    is enabled by default, although this isn't a strict requirement.
+
+    :param K: 2D kernel matrix
+    :type K: numpy array
+    :param n: Number of kernel PCAs to return (default=2)
+    :type n: integer
+    :param centering: Whether to center the kernel matrix (default=True)
+    :type centering: bool
+
+    :return: array containing the principal components
+    :rtype: numpy array
+    """
+
+    assert K.shape[0] == K.shape[1], "ERROR: Square matrix required for Kernel PCA."
+    assert np.allclose(K, K.T, atol=1e-8), "ERROR: Symmetric matrix required for Kernel PCA."
+    assert n <= K.shape[0], "ERROR: Requested more principal components than matrix size."
+
+    size = K.shape[0]
+    print("K", list(K))
+
+    result = np.zeros((size, size), dtype="double")
+    print("result", list(result))
+
+
+    fortlib.fkpca(K.ctypes.data_as(POINTER(c_double)), ctypes.byref(c_int(size)), ctypes.byref(c_bool(centering)), result.ctypes.data_as(POINTER(c_double)))
+
+
+    print("result", list(result))
+
+    # pca = fortlib.fkpca(K.ctypes.data_as(c_float_p), ctypes.c_int(size), ctypes.c_bool(centering))
+    # print(pca)
+
+    # return pca[:n]
+
+if __name__ == "__main__":
+
+    a = ct.c_double(5)
+    b = fortlib.sum2(ct.byref(a))
+    print("sum function", b == 7.0)
+
+    # 
+    x = np.array([[2, 1], [1, 2]], dtype="double")
+
+    N = 3
+    a = np.random.rand(N, N)
+    kernel = np.tril(a) + np.tril(a, -1).T
+
+    print(x)
+    pca = kpca(kernel)
+
+
+
+
+

src/testlib/funclib.f90

@@ -0,0 +1,138 @@
+module funclib
+
+use iso_c_binding, only: c_double, c_int, c_bool
+implicit none
+
+contains
+
+function sum2(a) result(b) bind(c, name='sum2')
+
+    use iso_c_binding
+    implicit none
+
+    real(c_double), intent(in)  :: a
+    real(c_double)              :: b
+
+    b = a + 2.d0
+
+end function sum2
+
+subroutine fkpca(k, n, centering, kpca) bind(c)
+
+    implicit none
+
+   real(c_double), intent(in) :: k(n, n)
+   integer(c_int), intent(in) :: n
+   logical(c_bool), intent(in) :: centering
+   real(c_double), intent(out) :: kpca(n, n)
+
+   ! Eigenvalues
+   real, dimension(n) :: eigenvals
+
+   real, allocatable, dimension(:) :: work
+
+   integer :: lwork
+   integer :: info
+
+   integer :: i
+
+   real :: inv_n
+   real, allocatable, dimension(:) :: temp
+   real :: temp_sum
+
+   write(*,*) "hello world"
+
+   kpca(:, :) = k(:, :)
+
+   write(*,*) "hello world"
+   write(*,*) kpca
+   write(*,*) k
+   write(*,*) "hello world"
+
+   ! This first part centers the matrix,
+   ! basically Kpca = K - G@K - K@G + G@K@G, with G = 1/n
+   ! It is a bit hard to follow, sry, but it is very fast
+   ! and requires very little memory overhead.
+
+   if (centering) then
+
+       write(*,*) "Should be false"
+
+      inv_n = 1.0d0/n
+
+      allocate (temp(n))
+      temp(:) = 0.0d0
+
+      !$OMP PARALLEL DO
+      do i = 1, n
+         temp(i) = sum(k(i, :))*inv_n
+      end do
+      !$OMP END PARALLEL DO
+
+      temp_sum = sum(temp(:))*inv_n
+
+      !$OMP PARALLEL DO
+      do i = 1, n
+         kpca(i, :) = kpca(i, :) + temp_sum
+      end do
+      !$OMP END PARALLEL DO
+
+      !$OMP PARALLEL DO
+      do i = 1, n
+         kpca(:, i) = kpca(:, i) - temp(i)
+      end do
+      !$OMP END PARALLEL DO
+
+      !$OMP PARALLEL DO
+      do i = 1, n
+         kpca(i, :) = kpca(i, :) - temp(i)
+      end do
+      !$OMP END PARALLEL DO
+
+      deallocate (temp)
+
+   end if
+
+   ! This 2nd part solves the eigenvalue problem with the least
+   ! memory intensive solver, namely DSYEV(). DSYEVD() is twice
+   ! as fast, but requires a lot more memory, which quickly
+   ! becomes prohibitive.
+
+   ! Dry run which returns the optimal "lwork"
+   allocate (work(1))
+   call dsyev("V", "U", n, kpca, n, eigenvals, work, -1, info)
+   lwork = nint(work(1)) + 1
+   deallocate (work)
+
+   write(*,*) "what is l work"
+   write(*,*) lwork
+
+   ! Get eigenvectors
+   allocate (work(lwork))
+   call dsyev("V", "U", n, kpca, n, eigenvals, work, lwork, info)
+   deallocate (work)
+
+   if (info < 0) then
+
+      write (*, *) "ERROR: The ", -info, "-th argument to DSYEV() had an illegal value."
+
+   else if (info > 0) then
+
+      write (*, *) "ERROR: DSYEV() failed to compute an eigenvalue."
+
+   end if
+
+   ! This 3rd part sorts the kernel PCA matrix such that the first PCA is kpca(1)
+   kpca = kpca(:, n:1:-1)
+   kpca = transpose(kpca)
+
+   !$OMP PARALLEL DO
+   do i = 1, n
+      kpca(i, :) = kpca(i, :)*sqrt(eigenvals(n - i + 1))
+   end do
+   !$OMP END PARALLEL DO
+
+
+end subroutine fkpca
+
+end module