2
\$\begingroup\$

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:

Landscape defining the pluvial basins

Color-coded basin areas

Flagged cells (assigned 'wrong' color - see text.)

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
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.