Fazang: a reverse-mode automatic differentiation library

Announcements
21

Fazang now supports ODE sensitivity solution (based on CVODES) with respect to given parameters
without explicitly asking for Jacobian. For that the user-defined ODE must include right-hand-side (RHS) definitions for var type and real types, due to the lack of meta-programming.

User-defined RHS would be like this.

module ode_mod
  use fazang
  use, intrinsic :: iso_c_binding
  implicit none

  real(rk), parameter :: omega = 0.5d0
  real(rk), parameter :: d1 = 1.0d0
  real(rk), parameter :: d2 = 1.0d0

contains
  ! right-hand-side for data input
  subroutine eval_rhs(t, y, fy)
    implicit none
    real(c_double), intent(in) :: t, y(:)
    real(c_double), intent(inout) :: fy(size(y))
    fy(1) = y(2)
    fy(2) = sin(omega * d1 * d2 * t)
  end subroutine eval_rhs

  ! right-hand-side for var input with parameters
  ! y, p, and output fy must all be of var type
  subroutine eval_rhs_pvar(t, y, fy, p)
    implicit none
    real(c_double), intent(in) :: t
    type(var), intent(in) :: y(:), p(:)
    type(var), intent(inout) :: fy(size(y))
    fy(1) = y(2)
    fy(2) = sin(p(1) * p(2) * p(3) * t)
  end subroutine eval_rhs_pvar
end module ode_mod

Now we can solve the defined ODE

program cvodes_demo
  use ode_mod
  use fazang
  implicit none

  type(var) :: yt(2, 3)
  type(cvodes_tol) :: tol
  real(rk), parameter :: ts(3) = [1.2d0, 2.4d0, 4.8d0]
  real(rk), parameter :: y00(2) = [0.2d0, 0.8d0]
  type(var) :: param(3)
  real(rk) :: y0(2), ga(2)
  integer :: i, j

  y0 = y00                      ! init condition
  param = var([omega, d1, d2])  ! parameters
  tol = cvodes_tol(CV_BDF, 1.d-10, 1.d-10, 1000_8)

  yt = cvodes_sol(0.d0, y0, ts, param, eval_rhs,&
       & eval_rhs_pvar, tol)
! ...
end program cvodes_demo

Note that the call to cvodes_sol includes parameter argument param and two RHS functions, one for real input and one for var input.

The sensitivities are obtained the same way by calling grad and adj on the output and input, respectively.

  call yt(1, 1) % grad()
  write(*, *) "dy_1/ d_omega at time ts(1):", param(1)%adj()

See examples and user doc for details.

5 Likes
22

This seems like a REALLY interesting and useful project! For me a lot since I’m starting to learn and use automatic differentiation. I’ll give this project a try but I would also love fpm compatibility to easily add it to my projects

23

I’m hopelessly overwhelmed in past a few months and lack bandwidth to revisit this. Also I’m not familiar with fpm so any help is appreciated.

24

I’ve just realized this post was from a whole year ago! I’ll check if I can help you with fpm compatibility. I don’t have much experience with packages with this complexity (mostly the sundials bit), but this can also be a learning experience for me :slight_smile: