Проблема осадков: процедурная реализация

Question

Сегодня вечером я видел текст «Вызов дождя», опубликованный семь лет назад Хантером Макмилленом. Я не смотрел на его решение (кроме того, я понятия не имею о perl), но я был заинтригован проблемой, поэтому я решил ее решить. В этом случае мне кажется, что ООП — это немного излишне, поэтому я сделал это процедурным способом. С Фортраном. Я полагаю, что на языке высокого уровня должно быть возможно написать очень краткое решение. Тем не менее, длина кода Fortran в этом случае находится на том же уровне, что и у Perl.

Программа генерирует случайный ландшафт (больше похожий на космический микроволновый фон, чем на холмы сельской местности, но все же …) и сохраняет трехмерный график с высотой, уникальным целочисленным индексом для каждого бассейна и кодом ошибки, указывающим, выполняются ли условия проблемы. . По большей части программа вроде работает довольно быстро и хорошо, но:

Интересно, можно ли сделать (еще) быстрее цикл, который окрашивает карту, с помощью некоторой умной настройки;
Иногда на карте встречаются отдельные точки, где условия не выполняются (две эквивалентные ячейки имеют разные цвета). Мы говорим о нескольких точках, когда карта имеет 160 000 ячеек. Я предполагаю, что эти сбои могут возникать в непосредственной близости от седел или локальных экстремумов, где высота, представленная только с одинарной точностью, может оказаться одинаковой в двух смежных ячейках, тем самым нарушая предположения проблемы. Однако я еще не решил проблему.

Вот несколько примеров графиков в поле размером 400 x 400 ячеек:

Жду ваших советов и отзывов. Ниже приведен код и инструкции о том, как его запустить и построить график результатов. Для людей, не знакомых с простым синтаксисом современного Фортрана, я предлагаю следующее быстрое современное руководство по Fortran.

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

Код компилируется с gfortran -O3

Чтобы отобразить результат с помощью 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

0

Похожие записи:

Добавить комментарий Отменить ответ