Note
Function that returns the best subsurface from the final surface.
We are here because a non periodic resulting surface is required. To do that, a wider periodic surface is created, and it matches the required moments and acf.
However, sub-sampling the surface into a smaller surface that matches the required size will result in degraded moments and acf. Hence, several locations are tested to find the best subsurface.
Note that the right moments can always be obtained by substitution, respecting the order of heights. However, the acf will be slightly impacted.
subroutine sub_surf() !================================================================================================ !<@note Function that returns the best subsurface from the final surface. !< !< We are here because a non periodic resulting surface is required. To do that, a wider !< periodic surface is created, and it matches the required moments and acf. !< !< However, sub-sampling the surface into a smaller surface that matches the required size !< will result in degraded moments and acf. Hence, several locations are tested to find the !< best subsurface. !< !< Note that the right moments can always be obtained by substitution, respecting the order of heights. !< However, the acf will be slightly impacted. !< !<@endnote !------------------------------------------------------------------------------------------------ implicit none integer(kind=I4) :: i, w, h, l, we, he, dw, dh, ndw, ndh, idw, jdh, ib, ie, jb, je, best_idw, best_jdh integer(kind=I4) :: n_seek, inc_dw, inc_dh, nn_res, res_ratio real(kind=R8) :: best_acf, res_acf, size_ratio character(len=100) :: text type(MOMENT_STAT) :: m_res integer(kind=I4), dimension(1:2) :: best_ind integer(kind=I4), allocatable, dimension(:) :: order_tmp real(kind=R8), allocatable, dimension(:,:) :: sav_surf, surf_tmp, acf_tmp, tab_res_acf real(kind=R8), allocatable, dimension(:) :: tab_tmp read(JOB,*) n_seek ; LINE_READ = LINE_READ +1 ; write(SPY,*) "line: ", LINE_READ, "n_seek ", n_seek ! n_seek is the number of subsurfaces explored if ( PARAM%periodic ) return ! reset the FFTW configuration because from now on, the acf will be calculated on smaller surfaces call end_fftw3() ! extended surface size we = PARAM%width he = PARAM%height allocate( sav_surf( 1:we, 1:he ) ) ! save the extended surface sav_surf( 1:we, 1:he ) = PARAM%surf( 1:we, 1:he ) ! reset previous arrays deallocate( PARAM%surf ) deallocate( PARAM%imp_acf ) deallocate( PARAM%acf_surf ) call fftw_plan_with_nthreads( nthreads = 1 ) NB_THREADS_FFT = omp_get_max_threads() ! the new image size is the one of the subsurface, because it is the size defined by the user w = PARAM%sub_width h = PARAM%sub_height l = PARAM%sub_npts call tab_init_fftw3_real( long = w, larg = h, plan_flag = FFTW_MEASURE ) PARAM%width = w PARAM%height = h PARAM%npts = l allocate( PARAM%surf( 1:w, 1:h ) ) allocate( PARAM%acf_surf( 1:w, 1:h ) ) ! build a smaller set of heights that match the required moments call build_heights( vec_out = PARAM%vect_h(1:l), & ! use_fct_expo = ( PARAM%m_end%ku < 1.10 * PARAM%m_end%sk**2 + 1. ), & ! stats_in = PARAM%m_end, & ! lg = l ) ! allocate( PARAM%imp_acf ( 1:w, 1:h ) ) ! recalculate the theoretical acf, for the new size surface call calc_imp_acf( long = w, & ! IN larg = h, & ! IN apod = PARAM%apod, & ! IN tau1 = PARAM%l_acf1, & ! IN tau2 = PARAM%l_acf2, & ! IN alpha = log(PARAM%acf__z), & ! IN ang = PARAM%a_acf * PI_R8/180, & ! IN tab_acf = PARAM%imp_acf(1:w, 1:h)) ! OUT ! difference between old and new size: hence, there will be dw*dh possible locations ! for the subsurface dw = we - w dh = he - h ! initialize the best acf result (mean absolute difference) and subsurface locations best_acf = 1.e6_R8 best_idw = 0 best_jdh = 0 size_ratio = dw / dh ! Number of locations along x and y. It respects the total number to be explored, as set ! by the user and the size ratio of the surface ndw = nint( sqrt( n_seek / size_ratio ) ) ndh = nint( real(n_seek, kind=R8) / ndw ) ! ndw and ndh are modified to be mumtiples of the number of threads, but the product ! should not be too far from n_seek if ( size_ratio > 1. ) then ndw = NB_THREADS_FFT * ( int(ndw / NB_THREADS_FFT) + 1 ) ndh = NB_THREADS_FFT * ( int(ndh / NB_THREADS_FFT) ) else ndw = NB_THREADS_FFT * ( int(ndw / NB_THREADS_FFT) ) ndh = NB_THREADS_FFT * ( int(ndh / NB_THREADS_FFT) + 1 ) endif ! don't exceed the maximums ndw = min( ndw, dw ) ndh = min( ndh, dh ) ! increments for looping inc_dw = dw / ndw inc_dh = dh / ndh ! result storage allocate( tab_res_acf(1:ndw, 1:ndh) ) tab_res_acf(1:ndw, 1:ndh) = -1 allocate( surf_tmp( 1:w, 1:h ) ) allocate( acf_tmp( 1:w, 1:h ) ) allocate( order_tmp(1:l) ) allocate( tab_tmp(1:l) ) call progress_bar_terminal(val = 0, max_val = ndw * ndh, init = .true.) !$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(NB_THREADS_FFT) !$OMP DO SCHEDULE (STATIC, ndw/NB_THREADS_FFT) PRIVATE(ib, ie, jb, je, jdh, surf_tmp, tab_tmp, order_tmp, acf_tmp, i, m_res, res_acf, nn_res, res_ratio, text) do idw = 1, ndw ib = idw * inc_dw ie = ib + w - 1 do jdh = 1, ndh jb = jdh * inc_dh je = jb + h - 1 surf_tmp( 1:w, 1:h ) = sav_surf( ib:ie, jb:je ) !....................................................................................... tab_tmp(1:l) = reshape( surf_tmp(1:w, 1:h), [l] ) ! store the subsurface height order call init_order( order = order_tmp(1:l), & ! OUT n = l ) ! IN call sort_array2( tab_inout = tab_tmp(1:l), & ! INOUT tab0 = order_tmp(1:l), & ! INOUT n = l ) ! IN ! replace old heights with new ones that matches required moments do i = 1, l tab_tmp( order_tmp(i) ) = PARAM%vect_h(i) enddo surf_tmp(1:w, 1:h) = reshape( tab_tmp(1:l), [w, h] ) call calc_moments( tab = tab_tmp(1:l), & ! IN mx = m_res, & ! OUT nb_mom = 4 ) ! IN surf_tmp(1:w, 1:h) = ( surf_tmp(1:w, 1:h) - m_res%mu ) / m_res%si !....................................................................................... ! calculate the acf call acf_wiener( tab_in = surf_tmp(1:w, 1:h), & ! IN tab_out = acf_tmp(1:w, 1:h), & ! OUT w = w, & ! IN h = h, & ! IN multi_fft = .true. ) ! IN ! calculate the result (mean absolute difference between theoretical acf and calculated acf) call calc_res_acf( acf_surf = acf_tmp(1:w, 1:h), & ! IN imp_acf = PARAM%imp_acf(1:w, 1:h), & ! IN acf__z = PARAM%acf__z, & ! IN crit_acf = res_acf, & ! OUT w = w, & ! IN h = h ) ! IN tab_res_acf( idw, jdh ) = res_acf ! update progressbar !$OMP CRITICAL nn_res = count( tab_res_acf > 0 ) call progress_bar_terminal(val = nn_res, max_val = ndw * ndh, init = .false.) !$OMP END CRITICAL enddo enddo !$OMP END DO !$OMP END PARALLEL ! find the best subsurface best_ind(1:2) = minloc( tab_res_acf ) idw = best_ind(1) jdh = best_ind(2) ib = idw * ( dw / ndw ) ie = ib + w - 1 jb = jdh * ( dh / ndh ) je = jb + h - 1 PARAM%surf( 1:w, 1:h ) = sav_surf( ib:ie, jb:je ) ! redo calculations on the best subsurface !....................................................................................... tab_tmp(1:l) = reshape( PARAM%surf(1:w, 1:h), [l] ) call init_order( order = PARAM%order(1:l), & ! n = l ) ! call sort_array2( tab_inout = tab_tmp(1:l), & ! tab0 = PARAM%order(1:l), & ! n = l ) ! do i = 1, l tab_tmp( PARAM%order(i) ) = PARAM%vect_h(i) enddo PARAM%surf(1:w, 1:h) = reshape( tab_tmp(1:l), [w, h] ) call calc_moments( tab = tab_tmp(1:l), & ! mx = m_res, & ! nb_mom = 4 ) ! PARAM%surf(1:w, 1:h) = ( PARAM%surf(1:w, 1:h) - m_res%mu ) / m_res%si !....................................................................................... call tab_end_fftw3_real() deallocate( sav_surf ) deallocate( tab_tmp ) deallocate( surf_tmp ) deallocate( acf_tmp ) deallocate( order_tmp ) deallocate( tab_res_acf ) return endsubroutine sub_surf