subroutine error_analysis(tSingPartPhase, iShiftVary, mean_ProjE_re, ProjE_Err_re, &
mean_ProjE_im, ProjE_Err_im, mean_Shift, Shift_Err, tNoProjEValue, tNoShiftValue, &
equilshift, equilproje)
implicit none
logical, intent(in) :: tSingPartPhase
integer, intent(inout) :: iShiftVary
integer, intent(in), optional :: equilshift, equilproje
integer :: corrlength_denom, corrlength_num, corrlength_shift, corrlength_imnum
integer :: equil_time_denom, equil_time_num, equil_time_shift, equil_time_imnum
real(dp), allocatable :: temp(:)
integer :: final_corrlength_proj, ierr
real(dp) :: mean_denom, error_denom, eie_denom
real(dp) :: mean_num, error_num, eie_num
real(dp) :: eie_shift
real(dp) :: mean_imnum, error_imnum, eie_imnum
real(dp) :: covariance_re, correction_re
real(dp) :: covariance_im, correction_im
logical :: tFail_num, tFail_denom, tFail_imnum, tFail_shift
logical :: tPrintIntermediateBlocking
logical, intent(out) :: tNoProjEValue, tNoShiftValue
real(dp), intent(out) :: mean_ProjE_re, ProjE_Err_re, mean_ProjE_im, ProjE_Err_im, mean_Shift, Shift_Err
character(len=*), parameter :: t_r = 'Error_analysis'
logical :: tFailRead
if (iProcIndex /= Root) return !Just do this in serial
write(stdout, "(A)")
write(stdout, "(A)") "Calculating approximate errorbar for energy estimates..."
write(stdout, "(A)")
if (present(equilshift) .or. present(equilproje)) then
tGivenEquilibrium = .true.
if (present(equilshift)) then
Shift_Start = equilshift
else
Shift_Start = 0
end if
if (present(equilproje)) then
ProjE_Start = equilproje
else
ProjE_Start = 0
end if
iShiftVary = min(ProjE_Start, Shift_Start)
else
tGivenEquilibrium = .false.
end if
tNoProjEValue = .false.
tNoShiftValue = .false.
if ((tSingPartPhase .and. (.not. tGivenEquilibrium)) .or. (SftDamp < 1.0e-7_dp)) then
write(stdout, "(A)") "Calculation has not entered variable shift phase. Error analysis therefore not performed."
write(stdout, "(A)") "Direct reblocking of instantaneous energy possible, but this would contain finite sampling bias."
tNoProjEValue = .true.
tNoShiftValue = .true.
return
else
write(stdout, "(A,I12)") "Attempting automatic reblocking analysis on data from iteration: ", iShiftVary
end if
!Read and store numerator_data, pophf_data, shift_data, and if complex, also imnumerator_data
call read_fcimcstats(iShiftVary, tFailRead)
if (tFailRead) then
tNoProjEValue = .true.
tNoShiftValue = .true.
return
end if
! STEP 2) Performs an preliminary blocking analysis, with a fixed blocklength of 2
! This is predominantly to yield the correlation length
! Files for plotting (for manual reblocking) are outputted in PLOT files
if (Errordebug > 0) then
tPrintIntermediateBlocking = .true.
else
tPrintIntermediateBlocking = .false.
end if
call automatic_reblocking_analysis(pophf_data, 2, corrlength_denom, &
'PLOT_denom', tPrintIntermediateBlocking, 1)
call automatic_reblocking_analysis(numerator_data, 2, corrlength_num, &
'PLOT_num', tPrintIntermediateBlocking, 2)
call automatic_reblocking_analysis(shift_data, 2, corrlength_shift, &
'PLOT_shift', tPrintIntermediateBlocking, 3)
if (lenof_sign == 2 .and. inum_runs == 1) then
call automatic_reblocking_analysis(imnumerator_data, 2, &
corrlength_imnum, 'PLOT_imnum', tPrintIntermediateBlocking, 4)
end if
if (.not. tGivenEquilibrium) then
! STEP 3) Attempts an automatic removal of equilibration time
! From the data that's had serial correlation removed
! (See subroutine for how this is achieved)
allocate(temp(size(pophf_data)), stat=ierr)
if (ierr /= 0) call stop_all(t_r, 'Alloc error')
temp = pophf_data
call reblock_data(temp, corrlength_denom)
if (errordebug > 0) call print_vector(temp, filename='PLOT_blocked_pophf')
call attempt_detect_equil_time(temp, equil_time_denom, tFail_denom)
if (tFail_denom) then
write(stdout, "(A)")
write(stdout, "(A)") "*** ERROR *** Failure to automatically detect equilibration time for " &
& //"projected energy denominator"
write(stdout, "(A)") "Skipping blocking analysis of projected energy, and energy estimate " &
& //"will be simple average over "
write(stdout, "(A)") "all iterations (including growth phase), which may contain correlated " &
& //"sampling bias. Use with caution."
write(stdout, "(A)") "Manual reblocking or continued running suggested for accurate " &
& //"projected energy estimate."
tNoProjEValue = .true.
end if
deallocate(temp)
allocate(temp(size(numerator_data)), stat=ierr)
if (ierr /= 0) call stop_all(t_r, 'Alloc error')
temp = numerator_data
call reblock_data(temp, corrlength_num)
if (errordebug > 0) call print_vector(temp, filename='PLOT_blocked_num')
call attempt_detect_equil_time(temp, equil_time_num, tFail_num)
if (tFail_num) then
write(stdout, "(A)")
write(stdout, "(A)") "*** ERROR *** Failure to automatically detect equilibration time for " &
& //"projected energy numerator"
write(stdout, "(A)") "Skipping blocking analysis of projected energy, and energy estimate " &
& //"will be simple average over "
write(stdout, "(A)") "all iterations (including growth phase), which may contain correlated " &
& //"sampling bias. Use with caution."
write(stdout, "(A)") "Manual reblocking or continued running suggested for accurate " &
& //"projected energy estimate."
tNoProjEValue = .true.
end if
deallocate(temp)
allocate(temp(size(shift_data)), stat=ierr)
if (ierr /= 0) call stop_all(t_r, 'Alloc error')
temp = shift_data
call reblock_data(temp, corrlength_shift)
if (errordebug > 0) call print_vector(temp, filename='PLOT_blocked_shift')
call attempt_detect_equil_time(temp, equil_time_shift, tFail_shift)
if (tFail_shift) then
write(stdout, "(A)")
write(stdout, "(A)") "*** ERROR *** Failure to automatically detect equilibration time for shift value."
write(stdout, "(A)") "Skipping blocking analysis and calculation of average shift."
write(stdout, "(A)") "Continued running suggested if accurate shift estimate required."
tNoShiftValue = .true.
end if
deallocate(temp)
if (lenof_sign == 2 .and. inum_runs == 1) then
allocate(temp(size(imnumerator_data)), stat=ierr)
if (ierr /= 0) call stop_all(t_r, 'Alloc error')
temp = imnumerator_data
call reblock_data(temp, corrlength_imnum)
if (errordebug > 0) call print_vector(temp, filename='PLOT_blocked_imnum')
call attempt_detect_equil_time(temp, equil_time_imnum, tFail_imnum)
if (tFail_imnum) then
write(stdout, "(A)")
write(stdout, "(A)") "*** ERROR *** Failure to automatically detect equilibration time for " &
& //"imaginary projected energy numerator"
write(stdout, "(A)") "Skipping blocking analysis of projected energy, and energy estimate " &
& //"will be simple average over "
write(stdout, "(A)") "all iterations (including growth phase), which may contain correlated " &
& //"sampling bias. Use with caution."
write(stdout, "(A)") "Continued running suggested for accurate projected energy estimate."
tNoProjEValue = .true.
end if
deallocate(temp)
end if
! STEP 4) Removes the equilibration time from the start of the data set
! Allow shift and projected energy to have seperate equilibration times
! If complex, the equilibration time for projected energy is the longest of all three quantities
if (.not. tNoProjEValue) then
relaxation_time_proje = max(equil_time_denom * corrlength_denom, equil_time_num * corrlength_num)
if (lenof_sign == 2 .and. inum_runs == 1) then
relaxation_time_proje = max(relaxation_time_proje, equil_time_imnum * corrlength_imnum)
end if
write(stdout, "(A,I8,A)") "Relaxation time for projected energy estimated to be ", relaxation_time_proje, &
" update cycles."
end if
if (.not. tNoShiftValue) then
relaxation_time_shift = equil_time_shift * corrlength_shift
write(stdout, "(A,I8,A)") "Relaxation time for shift estimated to be ", relaxation_time_shift, &
" update cycles."
end if
else
write(stdout, "(A,I8,A)") "Relaxation time for projected energy estimated to be ", relaxation_time_proje, &
" update cycles."
write(stdout, "(A,I8,A)") "Relaxation time for shift estimated to be ", relaxation_time_shift, " update cycles."
end if !Endif automatic relaxation time calculation
if (.not. tNoProjEValue) then
call resize(pophf_data, relaxation_time_proje)
call resize(numerator_data, relaxation_time_proje)
if (lenof_sign == 2 .and. inum_runs == 1) call resize(imnumerator_data, relaxation_time_proje)
!Now find all block lengths, and write them to file, so that blocklengths can be checked.
!Call routine here to write out a file (Blocks_proje) with the projected energy mean and error for all block sizes.
if (lenof_sign == 2 .and. inum_runs == 1) then
call print_proje_blocks(pophf_data, numerator_data, 2, 'Blocks_proje_re')
#ifdef CMPLX_
call print_proje_blocks(pophf_data, imnumerator_data, 2, 'Blocks_proje_im')
#endif
else
call print_proje_blocks(pophf_data, numerator_data, 2, 'Blocks_proje')
end if
!Now perform automatic reblocking again, to get expected blocklength
call automatic_reblocking_analysis(pophf_data, 2, corrlength_denom, 'Blocks_denom', .true., 1)
call automatic_reblocking_analysis(numerator_data, 2, corrlength_num, 'Blocks_num', .true., 2)
if (lenof_sign == 2 .and. inum_runs == 1) then
call automatic_reblocking_analysis(imnumerator_data, 2, corrlength_imnum, 'Blocks_imnum', .true., 4)
end if
end if
if (.not. tNoShiftValue) then
call resize(shift_data, relaxation_time_shift)
call automatic_reblocking_analysis(shift_data, 2, corrlength_shift, 'Blocks_shift', .true., 3)
end if
! STEP 5) Now gathers together the properly reblocked data and find statistics
if (.not. tNoProjEValue) then
final_corrlength_proj = max(corrlength_num, corrlength_denom)
if (lenof_sign == 2 .and. inum_runs == 1) then
final_corrlength_proj = max(final_corrlength_proj, corrlength_imnum)
end if
if (Errordebug > 0) write(stdout, "(A,I10)") "Final projected energy correlation length", final_corrlength_proj
call reblock_data(pophf_data, final_corrlength_proj)
call reblock_data(numerator_data, final_corrlength_proj)
call analyze_data(pophf_data, mean_denom, error_denom, eie_denom)
call analyze_data(numerator_data, mean_num, error_num, eie_num)
write(stdout, "(A,F22.10,A,G20.8,A,G22.10)") "ProjE_denominator:", mean_denom, " +/- ", error_denom, &
" Relative error: ", abs(error_denom / mean_denom)
if (lenof_sign == 2 .and. inum_runs == 1) then
call reblock_data(imnumerator_data, final_corrlength_proj)
call analyze_data(imnumerator_data, mean_imnum, error_imnum, eie_imnum)
write(stdout, "(A,F22.10,A,G20.8,A,G22.10)") "ProjE_numerator (Re):", mean_num, " +/- ", error_num, &
" Relative error: ", abs(error_num / mean_num)
write(stdout, "(A,F22.10,A,G20.8,A,G22.10)") "ProjE_numerator (Im):", mean_imnum, " +/- ", error_imnum, &
" Relative error: ", abs(error_imnum / mean_imnum)
else
write(stdout, "(A,F22.10,A,G20.8,A,G22.10)") "ProjE_numerator: ", mean_num, " +/- ", error_num, &
" Relative error: ", abs(error_num / mean_num)
end if
end if
! write(stdout,*) "OVERALL", mean2/mean1, "+-",
!abs(mean2/mean1)*((abs(error1/mean1))**2.0_dp+(abs(error2/mean2))**2.0_dp)**0.5_dp
! this is before the covariance, so don't print it now
if (.not. tNoShiftValue) then
!Now, also do this for the shift estimate
call reblock_data(shift_data, corrlength_shift)
call analyze_data(shift_data, mean_shift, shift_Err, eie_shift)
if (Errordebug > 0) then
write(stdout, "(A,F20.10,A,G20.8,A,G22.10)") "Shift: ", &
mean_shift, " +/- ", shift_Err, " Relative error: ", abs(shift_Err / mean_shift)
end if
end if
! STEP 6) Refine statistics using covariance
if (.not. tNoProjEValue) then
covariance_re = calc_covariance(pophf_data, numerator_data)
correction_re = 2.0_dp * covariance_re / ((size(pophf_data)) * mean_denom * mean_num)
if (lenof_sign == 2 .and. inum_runs == 1) then
covariance_im = calc_covariance(pophf_data, imnumerator_data)
correction_im = 2.0_dp * covariance_im / ((size(pophf_data)) * mean_denom * mean_imnum)
end if
mean_ProjE_re = mean_num / mean_denom
ProjE_Err_re = abs(mean_ProjE_re) * sqrt((abs(error_denom / mean_denom))**2.0_dp &
+ (abs(error_num / mean_num))**2.0_dp - correction_re)
if (lenof_sign / inum_runs == 2) then
mean_ProjE_im = mean_imnum / mean_denom
ProjE_Err_im = abs(mean_ProjE_im) * sqrt((abs(error_denom / mean_denom))**2.0_dp &
+ (abs(error_imnum / mean_imnum))**2.0_dp - correction_im)
if (ErrorDebug > 0) then
write(stdout, "(A,F20.8)") "Covariance correction (Re):", correction_re
write(stdout, "(A,F20.8)") "Covariance correction (Im):", correction_im
write(stdout, "(A,F20.10,A,G20.8)") "Final projected energy (Re): ", mean_ProjE_re, " +/- ", ProjE_Err_re
write(stdout, "(A,F20.10,A,G20.8)") "Final projected energy (Im): ", mean_ProjE_im, " +/- ", ProjE_Err_im
end if
else
if (ErrorDebug > 0) then
write(stdout, "(A,F20.8)") "Covariance correction:", correction_re
write(stdout, "(A,F20.10,A,G20.8)") "Final projected energy: ", mean_ProjE_re, " +/- ", ProjE_Err_re
end if
end if
end if
deallocate(pophf_data, numerator_data, shift_data)
if (lenof_sign == 2 .and. inum_runs == 1) deallocate(imnumerator_data)
end subroutine error_analysis