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_modinteger, private :: nd integer, dimension (:), allocatable, private :: x real, dimension (:), allocatable, private :: w logical, dimension (:), allocatable, private :: m

contains

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.

11/11/1997