My code is limited in speed by the following rebinning subroutine. Code is also here on Github: futils/futils.f90 at main · Nicholaswogan/futils · GitHub . Suppose you have values defined on some bins. This subroutine can rebin to a new set of bins.
I’m wondering if someone knows of an implementation that is faster? or if they can spot optimizations that I could make to the current version? Thanks!
subroutine rebin(old_bins, old_vals, new_bins, new_vals, ierr)
real(dp), intent(in) :: old_bins(:)
real(dp), intent(in) :: old_vals(:)
real(dp), intent(in) :: new_bins(:)
real(dp), intent(out) :: new_vals(:)
integer, optional, intent(out) :: ierr
integer :: i, j, l, n_old, n_new
real(dp) :: b1, b2
n_old = size(old_vals)
n_new = size(new_vals)
! option to check inputs.
if (present(ierr)) then
ierr = 0
! check shape
if (n_old+1 /= size(old_bins)) then
ierr = -1
return
endif
if (n_new+1 /= size(new_bins)) then
ierr = -2
return
endif
! check that new bins overlap with the old bins
if (new_bins(1) < old_bins(1) .and. new_bins(2) < old_bins(1)) then
ierr = -3
return
endif
if (new_bins(n_new) > old_bins(n_old+1) .and. new_bins(n_new+1) > old_bins(n_old+1)) then
ierr = -4
return
endif
! check that bins are all increasing
do i = 1,n_old
if (old_bins(i+1) <= old_bins(i)) then
ierr = -5
return
endif
enddo
do i = 1,n_new
if (new_bins(i+1) <= new_bins(i)) then
ierr = -6
return
endif
enddo
endif
new_vals = 0.0_dp
l = 1
do i = 1,n_new
b1 = new_bins(i+1) - new_bins(i)
do j = l,n_old
! four different cases
! ______ (old bin)
! ________ (new bin)
if (old_bins(j) <= new_bins(i) .and. &
old_bins(j+1) > new_bins(i) .and. old_bins(j+1) < new_bins(i+1)) then
b2 = old_bins(j+1) - new_bins(i)
new_vals(i) = new_vals(i) + (b2/b1)*old_vals(j)
! ____ (old bin)
! ________ (new bin)
elseif (old_bins(j) >= new_bins(i) .and. old_bins(j+1) <= new_bins(i+1)) then
b2 = old_bins(j+1) - old_bins(j)
new_vals(i) = new_vals(i) + (b2/b1)*old_vals(j)
! ____ (old bin)
! ________ (new bin)
elseif (old_bins(j) >= new_bins(i) .and. &
old_bins(j) < new_bins(i+1) .and. old_bins(j+1) > new_bins(i+1)) then
b2 = new_bins(i+1) - old_bins(j)
new_vals(i) = new_vals(i) + (b2/b1)*old_vals(j)
! ________ (old bin)
! ____ (new bin)
elseif (old_bins(j) < new_bins(i) .and. old_bins(j+1) > new_bins(i+1)) then
new_vals(i) = new_vals(i) + old_vals(j)
! _____ (old bin)
! ____ (new bin)
elseif (old_bins(j) > new_bins(i+1)) then
! time to move onto the next new_bin
l = j - 1
exit
else
! nothing
endif
enddo
enddo
end subroutine
@ivanpribec, this re-binning routine is actually a bottleneck now in my code. It takes more time than sorting. I have a little benchmark here: GitHub - Nicholaswogan/futils
Build with cmake, then run ./tests/test_futils. I get
Time to rebin (s): 2.4625E-07
Yesterday I was able to optimize it a little bit. But not enough to make a ton of difference.
In practice, the function is for down-binning data defined on bins (i.e. piecewise functions). An example is spectra of the infrared radiation emitted by Earth to space. I can calculate the emission of infrared photons within many tiny wavelength bins (R = 1000 case in figure below), and then I can re-bin these results to a lower resolution (orange line is rebinning of R = 1000 case).
Is the calculation possibly related to a moving average for smoothing (with a given window function) and so FFT could be used for convolution calculation?
You could try with just the do j = 1,n_old loop and keep a counter for the current new bin. Each old bin either fits entirely into the current new bin or there is a remainder to distributed over (perhaps several) following new bins. I think it would be sufficient to carry forward the density (value/width) and upper limit of the undistributed bin.
I don’t have time for a more detailed response right now. Perhaps later if you need more detail.
How about this one?
I’m using temporary variables.
Time to rebin (s), Original: 9.5300E-07
Time to rebin (s), Temp Var: 7.5267E-07 ←
subroutine rebin(old_bins, old_vals, new_bins, new_vals, ierr)
real(dp), intent(in) :: old_bins(:)
real(dp), intent(in) :: old_vals(:)
real(dp), intent(in) :: new_bins(:)
real(dp), intent(out) :: new_vals(:)
integer, optional, intent(out) :: ierr
integer :: i, j, l, n_old, n_new
real(dp) :: b2, b1_inv, v_new
real(dp) :: b_old0, b_old1
real(dp) :: b_new0, b_new1
real(dp) :: v_old
n_old = size(old_vals)
n_new = size(new_vals)
! option to check inputs.
if (present(ierr)) then
ierr = 0
! check shape
if (n_old+1 /= size(old_bins)) then
ierr = -1
return
endif
if (n_new+1 /= size(new_bins)) then
ierr = -2
return
endif
! check that new bins overlap with the old bins
if (new_bins(1) < old_bins(1) .and. new_bins(2) < old_bins(1)) then
ierr = -3
return
endif
if (new_bins(n_new) > old_bins(n_old+1) .and. new_bins(n_new+1) > old_bins(n_old+1)) then
ierr = -4
return
endif
! check that bins are all increasing
do i = 1,n_old
if (old_bins(i+1) <= old_bins(i)) then
ierr = -5
return
endif
enddo
do i = 1,n_new
if (new_bins(i+1) <= new_bins(i)) then
ierr = -6
return
endif
enddo
endif
new_vals = 0.0_dp
l = 1
do i = 1,n_new
b_new0 = new_bins(i)
b_new1 = new_bins(i+1)
b1_inv = 1d0/(b_new1 - b_new0)
v_new = new_vals(i)
do j = l,n_old
b_old0 = old_bins(j)
b_old1 = old_bins(j+1)
v_old = old_vals(j)
! Several different cases
! ____ (old bin)
! ________ (new bin)
if (b_old0 > b_new0 .and. b_old1 < b_new1) then
b2 = b_old1 - b_old0
v_new = v_new + (b2*b1_inv)*v_old
! ______ (old bin)
! ________ (new bin)
elseif (b_old0 <= b_new0 .and. &
b_old1 > b_new0 .and. b_old1 < b_new1) then
b2 = b_old1 - b_new0
v_new = v_new + (b2*b1_inv)*v_old
! ____ (old bin)
! ________ (new bin)
elseif (b_old0 >= b_new0 .and. &
b_old0 < b_new1 .and. b_old1 > b_new1) then
b2 = b_new1 - b_old0
v_new = v_new + (b2*b1_inv)*v_old
! need to move to the next new_bin
l = j
exit
! ________ (old bin)
! ____ (new bin)
elseif (b_old0 <= b_new0 .and. b_old1 >= b_new1) then
v_new = v_new + v_old
! need to move to the next new_bin
l = j
exit
! There are two remaining cases:
! _____ (old bin)
! ____ (new bin)
!
! and
! ____ (old bin)
! _____ (new bin)
!
! These cases can occur on either end of re-binning.
! However this is OK. The part of old_bin that does not
! overlap with new_bin on the ends, will just be ignored.
endif
enddo
new_vals(i) = v_new
enddo
end subroutine
@NaoMatch Thanks so much! This is faster on my computer too.
Why though?? Why is using temporary variables faster? Why doesn’t the compiler make this optimization for me? I suppose it uses slightly more stack memory, which might not be wanted in some special circumstances.
You have a j-loop within an i-loop. For each pair of values new_bins(i) and new_bins(i+1), you are searching linearly through the (probably large) set of values of old_bins(j) until you find the bin that matches the specified condition. Your bins are probably in monotonically increasing order. If so, using binary search would change the complexity of the search procedure from O(N) to O(lg N), where N is the number of bins.
Programming tricks such as using temporary variables, pointers, etc., can help speed up a program, but they rarely suffice to compensate for a poor choice of algorithm (such as linear search instead of binary search).
The nested loop is not as bad as it looks at first glance. The number of times the stuff in the inner
loop is evaluated (all the conditionals) is n_old. This is because the inner loop is being broken many times.
A binary search does not help either. In all my use cases, i know in advance that no “searching” is necessary, and that the inner loop will be relevant on j == 1
I’m pretty sure we are reaching maximal speed, for the operation I need to do. There is a possibility to make an approximation that is faster as suggested by @septc .
I think we need to help the compiler as much as the compiler helps us.
To use temporary variables is to help the compiler optimize.
The compiler is very smart, but not always smart. I think we should use simple temporary variables (ie: v_new) because they do not require much human effort. (Access to arrays should be kept as low as possible.)
Of course the smartest way is to improve the algorithm.
