add fortran tests

Author: scivision

.github/workflows/ci.yml

@@ -4,12 +4,26 @@ on:
   push:
     paths:
     - "**/*.py"
+    - "**/*.f90"
   pull_request:
     paths:
     - "**/*.py"
+    - "**/*.f90"
 
 jobs:
 
+  fortran:
+    runs-on: ubuntu-latest
+    steps:
+    - uses: actions/checkout@v1
+    - uses: actions/setup-python@v1
+      with:
+        python-version: '3.x'
+    - run: sudo apt install -yq --no-install-recommends gfortran ninja-build
+    - run: pip install meson
+    - run: meson setup build
+    - run: meson test -C build
+
   linux:
     runs-on: ubuntu-latest
     steps:

airtools/logmart.py

@@ -52,6 +52,8 @@ def logmart(A: np.ndarray, b: np.ndarray,
         raise ValueError('b must be all non-negative')
 
     b = b.copy()  # needed to avoid modifying outside this function!
+    # %% make sure there are no 0's in b
+    b[b <= 1e-8] = 1e-8
 # %% set defaults
     if x0 is None:  # backproject
         x = A.T @ b / A.sum()
@@ -62,8 +64,7 @@ def logmart(A: np.ndarray, b: np.ndarray,
         x = x0
     if not x.size == A.shape[1]:
         raise ValueError('x0 must be scalar or match Ncolumns of A')
-# %% make sure there are no 0's in b
-    b[b <= 1e-8] = 1e-8
+
     x[x < 1e-8] = 1e-8
     # W=sigma;
     # W=linspace(1,0,size(A,1))';

fortran/logmart.f90

@@ -1,10 +1,12 @@
 module art
 
+use iso_fortran_env, only: wp=>real64
+
 implicit none
 
 contains
 
-pure subroutine logmart(A,b,relax,x0,sigma,max_iter,x)
+pure subroutine logmart(A,b,relax,x0,sigma,max_iter, x)
 ! delta Chisquare.
 ! stopped if Chisquare increases.
 !
@@ -30,45 +32,38 @@ pure subroutine logmart(A,b,relax,x0,sigma,max_iter,x)
 !    x = [1,2,3]
 !    b = A*x
 
-use iso_fortran_env, only: wp=>real64
-
-
 ! --- parameter check
 real(wp), intent(in) :: A(:,:), b(:)
-real(wp), optional, value :: relax
-real(wp), intent(in),optional :: x0(:), sigma(:)
+real(wp), optional, value :: relax, sigma
+real(wp), intent(in),optional :: x0(:)
 integer, optional, value :: max_iter
 real(wp), intent(out) :: x(:)
 
-real(wp), dimension(size(b)) :: xA, op_sigma, W, op_b, arg,xold,c
+real(wp), dimension(size(b)) :: W(size(b)), x_prev,c, op_b
 integer :: i
-logical :: done
 real(wp) :: t,chi2,chiold
 
-op_b = b
 
 if (.not.size(A,1) == size(b)) error stop 'A and b row numbers must match'
+if (any(A<0)) error stop 'A must be non-negative'
+if (any(b<0)) error stop 'b must be non-negative'
+op_b = b
+! --- make sure there are no 0's in b
+where(op_b <= 1e-8) op_b = 1e-8_wp
 
 ! --- set defaults
-if (.not.present(relax)) relax = 1._wp
+if (.not.present(relax)) relax = 1
 if (.not.present(max_iter)) max_iter = 200
+if (.not.present(sigma))  sigma = 1
 
 if (.not.present(x0)) then
-  x  = matmul(transpose(A), b) / sum(A)
-  xA = matmul(A, x)
-  x  = x * maxval(b) / maxval(xA)
+  x  = matmul(transpose(A), op_b) / sum(A)
+  x  = x * maxval(op_b) / maxval(matmul(A, x))
 else
   x = x0
 endif
 
-if (.not.present(sigma)) then
-  op_sigma = 1._wp
-else
-  op_sigma = sigma
-endif
 
-! --- make sure there are no 0's in b
-where(op_b<=1e-8) op_b = 1e-8_wp
 
 ! W=sigma;
 ! W=linspace(1,0,size(A,1))';
@@ -77,27 +72,28 @@ pure subroutine logmart(A,b,relax,x0,sigma,max_iter,x)
 W = W / sum(W)
 
 ! --- iterate solution
-i=0
-done=.false.
-arg= ((matmul(A,x) - op_b) / op_sigma)**2
-chi2 = sqrt(sum(arg))
-
-do while (.not.done)
-  i = i+1
-  xold = x
-  xA = matmul(A,x)
-  t = minval(1/xA)
-  C = relax*t*(1-(xA/b))
+chi2 = chi_squared(A, op_b, x, sigma)
+
+do i = 1, max_iter
+  x_prev = x
+  t = minval(1/matmul(A,x))
+  C = relax*t*(1-(matmul(A,x)/op_b))
   x = x / (1-x*matmul(transpose(A),W*C))
 ! monitor solution
   chiold = chi2
-  chi2 = sqrt( sum(((xA - b)/op_sigma)**2) )
-  ! dchi2=(chi2-chiold)
-  done = ((chi2>chiold) .and. (i>2)) .or. (i==max_iter) .or. (chi2<0.7)
+  chi2 = chi_squared(A, op_b, x, sigma)
+  if (chi2 > chiold .and. i > 2) exit
 enddo
 
-x = xold
+x = x_prev
 
 end subroutine logmart
 
+
+pure real(wp) function chi_squared(A, b, x, sigma)
+real(wp), intent(in) :: A(:,:), b(:), x(:), sigma
+chi_squared = sqrt(sum(((matmul(A,x) - b) / sigma)**2))
+
+end function chi_squared
+
 end module art

fortran/random_utils.f90

@@ -8,15 +8,15 @@ module random_utils
 
 contains
 
-subroutine randn(noise)
+impure elemental subroutine randn(noise)
 ! implements Box-Muller Transform
 ! https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
 !
 ! Output:
 ! noise: Gaussian noise vector
 
-real(wp),intent(out) :: noise(:)
-real(wp),dimension(size(noise)) :: u1, u2
+real(wp),intent(out) :: noise
+real(wp) :: u1, u2
 
 call random_number(u1)
 call random_number(u2)

fortran/test_logmart.f90

@@ -8,21 +8,51 @@
 
 implicit none
 
-integer, parameter :: N=3
-real(wp) :: A(N,N)
-real(wp), parameter :: x_true(N)=[5,5,5]
-real(wp), parameter :: errtol=0.05_wp
-real(wp), dimension(N) :: x, noise, b,bias
-logical :: add_bias, add_noise
-
-add_bias = .false.
-add_noise = .true.
+integer, parameter :: N=4, M=4
+real(wp), dimension(N,M) :: A1, A2
+real(wp) :: x(M)
+logical :: ok=.true., add_noise = .true.
 
 call init_random_seed()
 
-A = reshape([1,0,0, &
-             0,1,0, &
-             0,0,1], shape(A), order=[2,1])
+A1 = reshape([5,0,0,0, &
+              0,5,0,0, &
+              0,0,5,0, &
+              0,0,0,5], shape(A1), order=[2,1])
+
+A2 = reshape([0,1,2,3, &
+              1,0,1,2, &
+              2,1,0,1, &
+              3,2,1,0], shape(A2), order=[2,1])
+
+x = [1._wp, 3._wp, 0.5_wp, 2._wp]
+
+if (.not. run_test(A1, x, 20,add_noise)) then
+  ok = .false.
+  write(stderr,*) 'failed on identity test'
+endif
+
+if (.not. run_test(A2, x, 2000, add_noise)) then
+  ok = .false.
+  write(stderr,*) 'failed on Fiedler test'
+endif
+
+if (.not. ok) error stop
+
+print *, 'OK: logmart'
+
+contains
+
+logical function run_test(A, x, max_iter, add_noise)
+
+real(wp), intent(in) :: A(:,:), x(:)
+logical, intent(in) :: add_noise
+integer, intent(in) :: max_iter
+
+real(wp), parameter :: errtol=0.05_wp
+real(wp), dimension(size(A,1)) :: noise, x_est, b
+
+run_test = .true.
 
 block
 integer :: i
@@ -33,30 +63,25 @@
 end block
 
 ! ---- noisy observation
-if (add_bias) then
-  call randn(bias)
-  bias = 0.01_wp * bias
-  print '(/,A,3F10.3)','bias',bias
-  A = A * spread(bias,2,N)
-endif
+b = matmul(A,x)
 
 if (add_noise) then
   call randn(noise)
   noise = 0.01_wp * noise
   print '(/,A,3F10.3)', 'noise',noise
-  b = matmul(A,x_true) + noise
+  b = b + noise
 endif
 
 ! ---- inversion
-call logmart(A,b,x=x)
+call logmart(A,b, max_iter=max_iter, x=x_est)
 
 ! --- check estimate
-if (any(abs(x-x_true) > errtol*maxval(x_true))) then
-  print *,x
+if (any(abs(x_est-x) > errtol*maxval(x))) then
+  print *,x_est
   write (stderr,*) 'larger than',errtol*100,' % error'
-  stop 1
+  run_test = .false.
 endif
 
-print '(/,A)','OK: logmart'
+end function run_test
 
 end program

matlab/logmart.m

@@ -26,6 +26,9 @@
 %% set defaults
 if (nargin<6), max_iter=200.; end
 if (nargin<5), sigma=1.; end
+%% make sure there are no 0's in y
+y(y<=1e-8)=1e-8;
+
 if (nargin<4) || isempty(x0)
     x=(A'*y)./sum(A(:));
     xA=A*x;
@@ -40,8 +43,7 @@
 validateattributes(relax, {'numeric'}, {'scalar', 'positive'})
 validateattributes(sigma, {'numeric'}, {'scalar', 'positive'})
 validateattributes(max_iter, {'numeric'}, {'scalar', 'positive'})
-%% make sure there are no 0's in y
-y(y<=1e-8)=1e-8;
+
 
 % W=sigma;
 % W=linspace(1,0,size(A,1))';