# Rainfall challenge: procedural implementation

This evening I have seen the text of the Rainfall Challenge posted seven years ago by Hunter McMillen. I did not look at his solution (besides, I have no clue about perl), but I was intrigued by the problem, so I set out to solve it. In this case, it seems to me that OOP is a bit of an overkill, so I did it in a procedural way. With Fortran. I imagine that in a high-level language it should be possible to write a very concise solution. Yet, the Fortran code length is in the same ballpark as Perl's in this case.

The program generates a random landscape (more like cosmic microwave background than countryside hills, but still ...) and save a 3D plot with elevation, an unique integer index per basin, and an error code specifying whether the conditions of the problem are met. For the most part, the program seems to work quite fast and well, but:

1. I wonder whether the cycle that colors the map can be made (still) faster with some clever tweak;

2. Occasionally, there are isolated points in the map where the conditions are not met (two equivalent cells have different colors). We are talking about a handful points when the map has 160,000 cells. I imagine this glitches may occur in proximity of saddles or local extrema, where the elevation, represented just in single precision, can happen to be identical in two contiguous cells, thus breaking the assumptions of the problem. However, I have not pinned down the issue yet.

Here are some sample plots in a field of 400 x 400 cells:   I look forward to your hints and feedback. Below is the code and instructions on how to run it and plot the results. For folks not familiar with the simple syntax of modern Fortran, I suggest the following quick modern Fortran tutorial.

program Basins

implicit none

integer, parameter :: XSIZE = 400
integer, parameter :: YSIZE = 400

real    :: land (0:XSIZE+1,0:YSIZE+1) = 0.
integer :: color(0:XSIZE+1,0:YSIZE+1) = 0
integer :: error(0:XSIZE+1,0:YSIZE+1) = 0

call CreateLand( land,               XSIZE, YSIZE )
call FindBasins( land, color,        XSIZE, YSIZE )
call CheckColor( land, color, error, XSIZE, YSIZE )
call SaveData  ( land, color, error, XSIZE, YSIZE )

contains

subroutine CreateLand( land, XSIZE, YSIZE )
real   , intent(out):: land( 0:XSIZE+1, 0:YSIZE+1 )
integer, intent(in) :: XSIZE, YSIZE
real   , parameter  :: PI = 4. * atan(1.)
complex, parameter  :: ZI = (0.,1.)
integer, parameter  :: NFOURX = 4
integer, parameter  :: NFOURY = 4
integer             :: ikx, iky, ix, iy
real                :: ck_abs, ck_arg, kx, ky
complex             :: ck, vx(0:XSIZE+1,1), vy(1,0:YSIZE+1)
vx = (0.,0.)
vy = (0.,0.)
do iky = -NFOURY, NFOURY
ky = 2 * PI / real(YSIZE) * real(iky)
do iy = 1, YSIZE
vy(1,iy) = exp( ZI * ky * real(iy) )
enddo
do ikx = -NFOURX, NFOURX
kx = 2 * PI / real(XSIZE) * real(ikx)
do ix = 1, XSIZE
vx(ix,1) = exp( ZI * kx * real(ix) )
enddo
call random_number(ck_abs)
call random_number(ck_arg)
ck = ck_abs * exp( ZI * ck_arg * 2. * PI ) / sqrt(real(ikx*ikx+iky*iky+1.))
land = land + real( ck * matmul(vx,vy) )
enddo
enddo
do ix = 1, XSIZE
land( ix,    0    ) = land( ix,   1   ) + 1.
land( ix, YSIZE+1 ) = land( ix, YSIZE ) + 1.
enddo
do iy = 1, YSIZE
land(    0   , iy ) = land(   1  , iy ) + 1.
land( XSIZE+1, iy ) = land( XSIZE, iy ) + 1.
enddo
end subroutine CreateLand

subroutine FindBasins( land, color, XSIZE, YSIZE )
real   , intent(in) :: land (0:XSIZE+1,0:YSIZE+1)
integer, intent(out):: color(0:XSIZE+1,0:YSIZE+1)
integer, intent(in) :: XSIZE, YSIZE
integer :: ix, iy, ix0, iy0, idx, idy
integer :: iColor, newColor, ncolors = 0
logical :: CHANGED
integer, allocatable :: colorVec(:)

!.. Cycles until all equivalent cells have the same color
do
CHANGED = .FALSE.
do iy = 1, YSIZE
do ix = 1, XSIZE
call FindLowestNeighbor( land(ix,iy), &
land(ix+1,iy), land(ix-1,iy) , &
land(ix,iy+1), land(ix,iy-1) , &
idx, idy )
ix0 = ix + idx
iy0 = iy + idy
call ReassignCellPairColor( &
land (ix0,iy0), land (ix,iy), &
color(ix0,iy0), color(ix,iy), nColors, CHANGED )
enddo
enddo
if(.not.CHANGED)exit
enddo

call ReMapColor( color, ncolors, XSIZE, YSIZE )

end subroutine FindBasins

subroutine FindLowestNeighbor( val, px, mx, py, my, idx, idy )
real   , intent(in) :: val, px, mx, py, my
integer, intent(out):: idx, idy
real :: lowest
lowest = val
idx = 0
idy = 0
if( px < lowest )then
idx = 1
lowest = px
endif
if( mx < lowest )then
idx = -1
lowest = mx
endif
if( py < lowest )then
idy = 1
lowest = py
endif
if( my < lowest )then
idy = -1
lowest = my
endif
end subroutine FindLowestNeighbor

subroutine ReassignCellPairColor( v0, v1, c0, c1, nc, CHANGED )
real   , intent(in)    :: v0, v1
integer, intent(inout) :: c0, c1, nc
logical, intent(inout) :: CHANGED
if( v0 < v1 )then
if( c1 == 0 )then
if( c0 == 0 )then
nc = nc + 1
c0 = nc
c1 = nc
else
c1 = c0
endif
CHANGED = .TRUE.
else
if( c0 == 0 )then
c0 = c1
CHANGED = .TRUE.
else
if( c0 /= c1 )then
c1 = min( c0, c1 )
c0 = c1
CHANGED = .TRUE.
endif
endif
endif
else
if( c1 == 0 )then
nc = nc + 1
c1 = nc
CHANGED = .TRUE.
endif
endif
end subroutine ReassignCellPairColor

subroutine ReMapColor( color, ncolors, XSIZE, YSIZE )
integer, intent(inout) :: color( 0:XSIZE+1, 0:YSIZE+1 )
integer, intent(in)    :: ncolors, XSIZE, YSIZE
integer :: colorVec(nColors), newColor, iColor, ix, iy
colorVec=0
do ix = 1, XSIZE
do iy = 1, YSIZE
colorVec(color(ix,iy))=1
enddo
enddo
newColor=0
do iColor = 1, nColors
newColor = newColor + colorVec(iColor)
colorVec(iColor) = newColor
end do
do ix = 1, XSIZE
do iy = 1, YSIZE
color(ix,iy) = colorVec(color(ix,iy))
enddo
enddo
end subroutine ReMapColor

subroutine CheckColor( land, color, error, XSIZE, YSIZE )
real   , intent(in) :: land (0:XSIZE+1,0:YSIZE+1)
integer, intent(in) :: color(0:XSIZE+1,0:YSIZE+1)
integer, intent(out):: error(0:XSIZE+1,0:YSIZE+1)
integer, intent(in) :: XSIZE, YSIZE
integer :: ix, iy, ix0, iy0, idx, idy
real    :: lowest
error = 0
do ix = 1, XSIZE
do iy = 1, YSIZE
call FindLowestNeighbor( land(ix,iy), &
land(ix+1,iy), land(ix-1,iy) , &
land(ix,iy+1), land(ix,iy-1) , &
idx, idy )
ix0 = ix + idx
iy0 = iy + idy
if( abs(idx) + abs(idy) > 0 .and. color(ix0,iy0) /= color(ix,iy) ) error(ix,iy) = 1
enddo
enddo
end subroutine CheckColor

subroutine SaveData( land, color, error, XSIZE, YSIZE )
real   , intent(in) :: land (0:XSIZE+1,0:YSIZE+1)
integer, intent(in) :: color(0:XSIZE+1,0:YSIZE+1)
integer, intent(in) :: error(0:XSIZE+1,0:YSIZE+1)
integer, intent(in) :: XSIZE, YSIZE
integer :: uid, ix, iy
open(newunit = uid, &
file    ="basin.txt")
do ix = 1, XSIZE
do iy = 1, YSIZE
write(uid,"(i0,x,i0,x,e14.6,x,i0,x,i0)") &
ix, iy, land(ix,iy), color(ix,iy), error(ix,iy)
enddo
write(uid,*)
enddo
close(uid)
end subroutine SaveData

end program


The code compiles with gfortran -O3

To plot the result with gnuplot:

set pm3d
set view map
set hidden3d
set size square
f='basin.txt'
splot f u 1:2:3 w l
splot f u 1:2:4 w l
splot f u 1:2:5 w l