Binning irregularly sampled data to a regular grid is another important operator, often found in practical applications. The adjoint of binning is interpolation, and the simplest case of interpolation is nearest-neighbor interpolation, where each data point is assigned with the nearest grid point value.
Irregular data typically come in ordinate-coordinate pairs. According to the minimal interface design, coordinates as well as the grid information belong to the initialization interface, while the corresponding ordinates go directly to the operator interface. At the initialization step, module bin1 transforms the grid and coordinate information into two data-size arrays: the integer array x with the nearest bin locations and the logical array m, which masks the data points outside of the grid. Additionally, there is a real array w with the weighting function, optionally supplied by the user.
module bin1 use adj_mod
integer, private :: nd integer, dimension (:), allocatable, private :: x real, dimension (:), allocatable, private :: w logical, dimension (:), allocatable, private :: m
subroutine bin1_init (coord, o1, d1, n1, weight) real, dimension (:), intent (in) :: coord, weight real, intent (in) :: o1, d1 integer, intent (in) :: n1 optional :: weight
nd = size (coord) if (.not. allocated (x)) allocate (x (nd), m (nd), w (nd))
x = 1.5 + (coord - o1) / d1 m = (x < 1 .or. x > n1)
if (present (weight)) then m = m .or. weight == 0. where (.not. m) w = weight else where (.not. m) w = 1. end if end subroutine bin1_init
function bin1_op (adj, add, mod, ord) result (stat) integer :: stat logical, intent (in) :: adj, add real, dimension (:) :: mod, ord
integer :: i
call adjnull (adj, add, mod, ord)
do i = 1, nd ; if (m (i)) cycle if (adj) then mod (x (i)) = mod (x (i)) + w (i) * ord ( i ) else ord ( i ) = ord ( i ) + w (i) * mod (x (i)) end if end do
stat = 0 end function bin1_op
subroutine bin1_close () deallocate (x, m, w) end subroutine bin1_close
end module bin1
The output of the interpolation operator is ``smooth'', each data point obtains a unique value from the nearest bin. In a shorthand Fortran 90 notation, we could write the forward operation in one line as where (.not. m) ord = ord + w * mod (x). However, writing the adjoint operation as where (.not. m) mod (x) = mod (x) + w * ord is not correct, because the output of the adjoint process (binning) is ``rough'': each bin may receive zero, one, or more than one input from the data.