Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Revert "Wrote unit test for Hilbert transform" #163

Merged
merged 1 commit into from
Nov 12, 2024
Merged

Revert "Wrote unit test for Hilbert transform" #163

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
47 changes: 1 addition & 46 deletions src/misc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ module misc
#endif

use precision, only: r128, r64, i64
use params, only: kB, twopi, pi
use params, only: kB, twopi

implicit none

Expand Down Expand Up @@ -1653,49 +1653,4 @@ subroutine subtitle(text)
string2print(75 - length + 1 : 75) = text
if(this_image() == 1) write(*,'(A75)') string2print
end subroutine subtitle

subroutine Hilbert_transform(fx, Hfx)
!! Does Hilbert tranform for a given function
!! Ref - EQ (4)3, R. Balito et. al.
!! - An algorithm for fast Hilbert transform of real
!! functions
!!
!! fx is the function
!! Hfx is the Hilbert transform of the function
!!

real(r64), allocatable, intent(out) :: Hfx(:)
real(r64), intent(in) :: fx(:)

! Local variables
integer :: n, k, nfx
real(r64) :: term2, term3, b

nfx = size(fx)
allocate(Hfx(nfx))

! Hilbert function is zero at the edges
Hfx(1) = 0.0_r64
Hfx(nfx) = 0.0_r64

! Both Hfx and fx follow 1-indexing system unlike the paper
do k = 1, nfx - 2 ! Run over the internal points
term2 = 0.0_r64 ! 2nd term in Bilato Eq. 4
term3 = 0.0_r64 ! 3rd term in Bilato Eq. 4

do n = 1, nfx - 2 - k ! Partial sum over internal points
b = log((n + 1.0_r64)/n)
term2 = term2 - (1.0_r64 - (n + 1.0_r64)*b)*fx(k + n + 1) + &
(1.0_r64 - n*b)*fx(k + n + 2)
end do

do n = 1, k - 1 ! Partial sum over internal points
b = log((n + 1.0_r64)/n)
term3 = term3 + (1.0_r64 - (n + 1.0_r64)*b)*fx(k - n + 1) - &
(1.0_r64 - n*b)*fx(k - n)
end do

Hfx(k + 1) = (-1.0_r64/pi)*(fx(k + 2) - fx(k) + term2 + term3)
end do
end subroutine Hilbert_transform
end module misc
50 changes: 48 additions & 2 deletions src/screening.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,7 @@ module screening_module
use numerics_module, only: numerics
use misc, only: linspace, mux_vector, binsearch, Fermi, print_message, &
compsimps, twonorm, write2file_rank2_real, write2file_rank1_real, &
distribute_points, sort, qdist, operator(.umklapp.), Bose, &
hilbert_transform
distribute_points, sort, qdist, operator(.umklapp.), Bose
use wannier_module, only: wannier
use delta, only: delta_fn, get_delta_fn_pointer

Expand Down Expand Up @@ -59,6 +58,53 @@ subroutine calculate_qTF(crys, el)
end if
end subroutine calculate_qTF

subroutine Hilbert_transform(f, Hf)
!! Does Hilbert tranform of spectral head of bare polarizability
!! Ref - EQ (4), R. Balito et. al.
!!
!! Hf H.f(x), the Hilbert transform
!! f(x) the function

real(r64), intent(in) :: f(:)
real(r64), allocatable, intent(out) :: Hf(:)

! Locals
real(r64) :: term2, term3, term1, b
integer(i64) :: nxs, k, n

!Number points on domain grid
nxs = size(f)

allocate(Hf(nxs))

!Assume that f vanishes at the edges, and Hf also
Hf(1) = 0.0_r64
Hf(nxs) = 0.0_r64

do k = 2, nxs - 1 !Run over the internal points
term2 = 0.0_r64 !2nd term in Bilato Eq. 4
term3 = 0.0_r64 !3rd term in Bilato Eq. 4

do n = 2, nxs - 1 - k !Partial sum over internal points
b = log((n + 1.0_r64)/n)

term2 = term2 - (1.0_r64 - (n + 1.0_r64)*b)*f(k + n) + &
(1.0_r64 - n*b)*f(k + n + 1)
end do

do n = 2, k - 2 !Partial sum over internal points
b = log((n + 1.0_r64)/n)

term3 = term3 + (1.0_r64 - (n + 1.0_r64)*b)*f(k - n) - &
(1.0_r64 - n*b)*f(k - n - 1)
end do

term1 = f(k + 1) - f(k - 1)

Hf(k) = (-1.0_r64/pi)*(term1 + term2 + term3)
end do
end subroutine Hilbert_transform

!!$ subroutine head_polarizability_real_3d_T(Reeps_T, Omegas, spec_eps_T, Hilbert_weights_T)
!!$ !! Head of the bare real polarizability of the 3d Kohn-Sham system using
!!$ !! Hilbert transform for a given set of temperature-dependent quantities.
Expand Down
86 changes: 2 additions & 84 deletions test/test_misc.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,12 +7,12 @@ program test_misc
twonorm, binsearch, mux_vector, demux_vector, interpolate, coarse_grained, &
unique, linspace, compsimps, mux_state, demux_state, demux_mesh, expm1, &
Fermi, Bose, Pade_continued, precompute_interpolation_corners_and_weights, &
interpolate_using_precomputed, operator(.umklapp.), shrink, hilbert_transform
interpolate_using_precomputed, operator(.umklapp.), shrink

implicit none

integer :: itest
integer, parameter :: num_tests = 32
integer, parameter :: num_tests = 28
type(testify) :: test_array(num_tests), tests_all
integer(i64) :: index, quotient, remainder, int_array(5), v1(3), v2(3), &
v1_muxed, v2_muxed, ik, ik1, ik2, ik3, ib1, ib2, ib3, wvmesh(3), &
Expand All @@ -23,9 +23,6 @@ program test_misc
real_array(5), result, q1(3, 4), q2(3, 4), q3(3, 4)
real(r64), allocatable :: integrand(:), domain(:), im_axis(:), real_func(:), &
widc(:, :), f_coarse(:), f_interp(:), array_of_reals(:)
real(r64), allocatable :: hfx1_even(:), hfx1_odd(:), hfx2_even(:), hfx2_odd(:), &
ind_even(:), ind_odd(:), x_even(:), x_odd(:), xmin, xmax
integer(i64) :: n_even, n_odd

print*, '<<module misc unit tests>>'

Expand Down Expand Up @@ -336,88 +333,9 @@ program test_misc
array_of_reals = [1, 2, 3, 4, 5]*1.0_r64
call shrink(array_of_reals, 2_i64)
call test_array(itest)%assert(array_of_reals, [1, 2]*1.0_r64)

! Hilbert transform tests (H)
! fx1 -> function 1, fx2 -> function 2
! - hfx1_even stores hilbert transform calculated for fx1, and for even number
! - of points
xmin = -30.0
xmax = 30.0
n_even = 4000
n_odd = 4001
! ind_even are indices to compare in case of even number of points
! ind_odd are indices to compare in case of odd number of points
allocate(ind_even(6),ind_odd(5))
ind_even = [801, 1201, 1601, 2001, 2401, 2801]
ind_odd = [889, 1333, 1777, 2221, 2665]

itest = itest + 1
test_array(itest) = testify("Hilbert transform: even function, even points")
allocate(x_even(n_even), hfx1_even(n_even))
call linspace(x_even, xmin, xmax, n_even)
call Hilbert_transform(fx1_array(x_even), hfx1_even)
call test_array(itest)%assert(hfx1_even(ind_even), hfx1_array(x_even(ind_even)), &
tol = 2e-4_r64)

itest = itest + 1
test_array(itest) = testify("Hilbert transform: odd function, even points")
allocate(hfx2_even(n_even))
call Hilbert_transform(fx2_array(x_even), hfx2_even)
call test_array(itest)%assert(hfx2_even(ind_even), hfx2_array(x_even(ind_even)), &
tol = 1e-4_r64)

itest = itest + 1
test_array(itest) = testify("Hilbert transform: even function, odd points")
allocate(x_odd(n_odd), hfx1_odd(n_odd))
call linspace(x_odd, xmin, xmax, n_odd)
call Hilbert_transform(fx1_array(x_odd), hfx1_odd)
call test_array(itest)%assert(hfx1_odd(ind_odd), hfx1_array(x_odd(ind_odd)), &
tol = 4e-4_r64)

itest = itest + 1
test_array(itest) = testify("Hilbert transform: odd function, odd points")
allocate(hfx2_odd(n_odd))
call Hilbert_transform(fx2_array(x_odd), hfx2_odd)
call test_array(itest)%assert(hfx2_odd(ind_odd), hfx2_array(x_odd(ind_odd)), &
tol = 1e-5_r64)

tests_all = testify(test_array)
call tests_all%report

if(tests_all%get_status() .eqv. .false.) error stop -1

contains
! reference functions for the Hilbert transform test
! fx1 is an even function
function fx1_array(x) result(fx1)
real(r64), intent(in) :: x(:)
real(r64), allocatable :: fx1(:)
allocate(fx1(size(x)))
fx1 = 1/(1 + x**2)
end function fx1_array

! Hfx1 is actual hilbert transform of fx1, is an odd function
function hfx1_array(x) result(hfx1)
real(r64), intent(in) :: x(:)
real(r64), allocatable :: hfx1(:)
allocate(hfx1(size(x)))
hfx1 = x/(1 + x**2)
end function hfx1_array

! fx2 is an odd function
function fx2_array(x) result(fx2)
real(r64), intent(in) :: x(:)
real(r64), allocatable :: fx2(:)
allocate(fx2(size(x)))
fx2 = sin(x)/(1 + x**2)
end function fx2_array

! Hfx2 is actual hilbert transform of fx2, is an even function
function hfx2_array(x) result(hfx2)
real(r64), intent(in) :: x(:)
real(r64), allocatable :: hfx2(:)
real(r64), parameter :: e = 2.718281
allocate(hfx2(size(x)))
hfx2 = (1/e - cos(x))/(1 + x**2)
end function hfx2_array
end program test_misc
Loading