!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
! reblock6.f90                                                      31/12/2004
!
! In-situ binary reblocking algorithm for error estimation (Flyvbjerg)
! demonstrated for mock Markov chain Monte Carlo simulation.
! USAGE:
!   gmake Reblock3; Reblock3
!   x energy.out
!   x -settype xydy error.out
!
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
! Mathematical constants and functions which might be useful elsewhere
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
MODULE maths
  IMPLICIT NONE
  INTEGER, PARAMETER :: DP=KIND(1.0d0)
  REAL(DP), PARAMETER :: PI=3.14159265358979323846_DP
END MODULE maths

!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
! Binary ReBlocking "class"
! As data comes in, we arrange it along a row ... leaf nodes.
! Whenever there are two root nodes at the same level l, 
! we merge them into a root node at level (l+1).
! Example: Consider N=11.  Suppose x(n)=n.  The nodes are labelled with their 
! weights.  There are three separate heaps:
!
! l=3               28
!                  / \         
!                 /   \
!                /     \
!               /       \
!              /         \
! l=2         6          22
!            / \         / \
!           /   \       /   \
! l=1      1     5     9    13      17
!         / \   / \   / \   / \     / \
! l=0    0   1 2   3 4   5 6   7   8   9   10
! 
! When node 11 is added, N changes from 11 to 12, i.e., from 1011b to 1100b.
! We then compute 10+11=21 and 17+21=38 and add them:
!
! l=3               28
!                  / \         
!                 /   \
!                /     \
!               /       \
!              /         \
! l=2         6          22             38
!            / \         / \           /  \
!           /   \       /   \         /    \
! l=1      1     5     9    13      17      21
!         / \   / \   / \   / \     / \    /  \
! l=0    0   1 2   3 4   5 6   7   8   9  10  11
! 
! For the binary reblocking application it is not necessary to maintain all the nodes.
! Instead, all one requires is:
! (a) the weight of the root node of each distinct heap.
!     In the case N=11, we have a({0,1,2,3}) = {10, 17, 0, 28}.
! (b) the sums of the squares of the weights of all nodes generated so far
!     at each level l.
!     In the case N=11, we have
!       b(0) = 0**2 + 1**2 + ... + 10**2
!       b(1) = 1**2 + 5**2 + ... + 17**2
!       b(2) = 6**2 + 22**2
!       b(3) = 28**2
!
!------------------------------------------------------------------------------
! In the MODULE brb, Lmax is the maximum height of the heap (that we plan for).
! Let 
!   N = number of data points 
!   L = block level (must be <= Lmax)
!   S = block size = 2**l.
!   M = number of blocks=FLOOR(N/S), 
!   a(l) = if there is a separate l-l heap, a_l is its root weight;
!         otherwise undefined (say zero)
!   c(l) = sum of weights of all l-l nodes (not used)
!   b(l) = sum of squares of weights of all l-l nodes
!
! Data collection (heap-building) stage: recursive version of algorithm
!   MAIN PROGRAM: add(x,n,l=0)
!   SUBROUTINE add(x,n,l=0)
!     b(l) += x**2
!     IF (n is odd)  add(a(l)+x, n/2, l+1);  ELSE;   a(l)=x; ENDIF
!   END SUBROUTINE
!
! Output stage:
!   SoB = sum of all complete blocks of size S 
!       = sum of weights of all existing heaps of level l or higher
!   MoBM = mean of all complete block means of size S 
!        = SoB/S/M
!   MoSoBM = Mean of Squares of Block Means
!        = buf%b(:,l)/S**2/M
!   VoBM = Variance of Block Means 
!        = (MoSoBM - (MoBM)**2)     *M/(M-1)     (note correction factor)
!   VoMoBM = Variance of Mean of Block Means = VoBM/M
!   EoMoBM = Error of Mean of Block Means = SQRT(VoMoBM)
!   EoEoMoBM = Error of Error of Mean of Block Means = EoMoBM/SQRT(2(M-1))
!
! Example: N=11=1101b       , L=1, S=2, M=5
!   a(0) = x12
!   a(1) = x8+x9
!   a(2) = ?
!   a(3) = x0+x1+x2+x3+x4+x5+x6+x7
!   a(4) = undef 
!        etc.
!   b(0) = x0**2 + x1**2 + ... + x12**2
!   b(1) = (x0+x1)**2 + ... + (x10+x11)**2
!   b(2) = (x0+x1+x2+x3)**2 + (x4+x5+x6+7)**2
!   b(3) = (x0+x1+x2+x3+x4+x5+x6+x7)**2
!
!   SoB(3) = a(3) = (x0+x1+x2+x3+x4+x5+x6+x7)
!   SoB(2) = a(3) = (x0+x1+x2+x3)+(x4+x5+x6+x7)
!   SoB(1) = a(3)+a(1) = (x0+x1)+(x2+x3)+(x4+x5)+(x6+x7)+(x8+x9)
!   SoB(0) = a(3)+a(1)+a(0) = x0+x1+x2+x3+x4+x5+x6+x7+x8+x9+x10
!
!             ( x0+x1+x2+x3         x4+x5+x6+x7     )
!   VoBM(2) = (-------------**2  + -------------**2 )
!             (      4                   4          )
!            ------------------------------------------
!                                2
!
!                (  x0+x1+x2+x3     x4+x5+x6+x7   )
!                ( ------------- + -------------  )
!                (       4               4        )
!            -   (------------------------------- )**2
!                (                2               )
!
!        MULTIPLIED BY THE CORRECTION FACTOR 2/(2-1)
!
! Ultimately we are interested in the mean of all values, MoB(l=0), and
! the error in this mean, EoMoBM, at a suitable value of l.
! The error in the error is only used to determine the autocorrelation "time.
! 
! The variance at level l is calculated using complete blocks only;
! for l<lmax-3, at least 7/8 of the data points fall within complete blocks,
! so let's just skate over the details here.
!
! Note: It is common for N to exceed 2**16=65536, so N**2 exceeds 2**32.
! Hence we cast integers to floating-points before squaring.
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
MODULE brb
  USE maths
  IMPLICIT NONE
  TYPE brb_buffer
     INTEGER :: Lmax,N,datasize
     REAL(DP), POINTER :: a(:,:), b(:,:)
  END TYPE brb_buffer
CONTAINS
  !======================================================================
  ! Allocate a buffer which is Lmax levels deep.
  !
  !======================================================================
  SUBROUTINE brb_allocate(buf, Lmax, datasize)
    TYPE(brb_buffer), INTENT(OUT) :: buf
    INTEGER, INTENT(IN) :: Lmax,datasize
    buf%datasize = datasize
    buf%Lmax = Lmax
    buf%N = 0
    ALLOCATE (buf%a(datasize, 0:Lmax))  !column-major:
    ALLOCATE (buf%b(datasize, 0:Lmax))  !first index changes most rapidly
    buf%b = 0.0_DP                      !a does not need to be initialised
  END SUBROUTINE brb_allocate
  !===========================================================================
  ! Append data vector (xx) to buffer, and update multi-heap
  !===========================================================================
  SUBROUTINE brb_append(buf, xx)
    TYPE(brb_buffer), INTENT(INOUT) :: buf
    REAL(DP), INTENT(IN) :: xx(buf%datasize)
    REAL(DP) :: yy(buf%datasize)
    INTEGER :: l,ivar

    ivar = buf%N   !size of buffer before incrementing
    yy = xx
    l = 0
    DO
       buf%b(:,l) = buf%b(:,l) + yy**2
       IF (IAND(ivar, 1)==0) THEN  !bit was 0
          buf%a(:,l) = yy
          EXIT !exit do loop
       ELSE !bit was 1: carry over to next level.  Assume Lmax is large enough
          yy = yy + buf%a(:,l)
          ivar = ISHFT(ivar, -1)
          l = l+1
       END IF
    END DO
    buf%N = buf%N+1
  END SUBROUTINE brb_append
  !===========================================================================
  ! Report contents of buffer.
  !===========================================================================
  SUBROUTINE brb_report(buf, xxmean, xxerr)
    TYPE(brb_buffer), INTENT(IN) :: buf
    REAL(DP), DIMENSION(buf%datasize), INTENT(OUT) :: xxmean,xxerr
    INTEGER :: l,S,N,M,Lmax_actual,key
    REAL(DP) :: fractchg
    REAL(DP), DIMENSION(buf%datasize) :: SoB,VoBM,VoMoBM
    REAL(DP), DIMENSION(buf%datasize,0:buf%Lmax) :: aEoMoBM,aEoEoMoBM

    N = buf%N
    IF (N<=1) THEN
       xxmean = -999.0_DP
       xxerr = -999.0_DP
       PRINT *,"Empty brb buffer"
       RETURN
    END IF

    !----- Find the level of the largest heap, Lmax_actual.
    DO l=0,buf%Lmax
       IF (IBITS(N, l, 1) /=0) Lmax_actual=l
    END DO

    !----- Accumulate sum from largest heap to smallest.
    ! Compute variances at each level.
    SoB = 0.0_DP
    DO l=Lmax_actual, 0, -1
       IF (IBITS(N, l, 1) /= 0)  SoB=SoB+buf%a(:,l)
       S=2**l; M=N/S; IF(M<=1) CYCLE  ! cannot do stats with only 1 block

       VoBM = (buf%b(:,l) - SoB**2/M) / (REAL(S,DP)**2 * (M-1))  !refactored
       VoMoBM = VoBM / M
       aEoMoBM(:,l) = SQRT(VoMoBM)
       aEoEoMoBM(:,l) = aEoMoBM(:,l) / SQRT(REAL(2*(M-1),DP))
       WRITE (20,"(I4,2F24.16)") l, aEoMoBM(1,l), aEoEoMoBM(1,l)
       WRITE (21,"(I4,2F24.16)") l, aEoMoBM(2,l), aEoEoMoBM(2,l)
    END DO

    !----- Now for the AI: Computer guesses true error
    ! As the block size is increased,
    ! compute the CHANGE in the error as a fraction of the ERROR in the error
    ! of the block means.  When this is negligible (less than 0.2, say),
    ! it means that the error estimate is unaffected by
    ! further increases of the block size; i.e., blocks are bigger than
    ! the correlation length.  Then we can return the error ...
    WRITE (*,"('Number of data values N = ',I24)") N
    DO key=1,buf%datasize
       DO l=0,Lmax_actual-2 !note this carefully
          fractchg = (aEoMoBM(key,l+1) - aEoMoBM(key,l)) / aEoEoMoBM(key,l)
          IF (fractchg < 0.25_DP) EXIT
       END DO
       xxmean(key) = SoB(key)/N
       xxerr(key) = aEoMoBM(key,l)
       WRITE (*,"(A,I3)")  "Optimal block size =", l
       WRITE (*,"(A,F24.16)")  "Mean =", xxmean(key)
       WRITE (*,"(A,F24.16)")  "Error =", xxerr(key)
       IF (l==Lmax_actual-2) &
            PRINT *,"Warning: Inter-block correlation remains!",&
            "Longer MC run may be required"
    END DO

  END SUBROUTINE brb_report
END MODULE brb






!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
! Main program: Example of using the brb package
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
PROGRAM hello
  USE maths
  USE brb
  USE nag_f77_g_chapter
  IMPLICIT NONE

  !%%%%%%%%%%%%%%%%%%%% 1. Declarations %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  INTEGER :: Lmax,Nmax,datasize,l,n
  REAL(DP) :: theta
  REAL(DP) :: costheta,sintheta,corrlen
  REAL(DP) :: costheta2,sintheta2,corrlen2
  REAL(DP), ALLOCATABLE :: xx(:),xxmean(:),xxerr(:)
  TYPE(brb_buffer) :: buf

  !%%%%%%%%%%%%%%%%%%%% 2. Initialisation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  OPEN (UNIT=10,STATUS="REPLACE",FILE="energy.out")
  OPEN (UNIT=20,STATUS="REPLACE",FILE="error.out")
  OPEN (UNIT=21,STATUS="REPLACE",FILE="error2.out")
  WRITE (*,"(A)") "==================================================="
  WRITE (*,"(A)") "In-situ reblocking for Monte Carlo error estimation"
  CALL G05CCF	         ! seed random number generator

  Lmax = 16   ! choose a safely big value
  Nmax = 65536  ! 2**Lmax

  datasize = 2
  ALLOCATE (xx(datasize))
  ALLOCATE (xxmean(datasize))
  ALLOCATE (xxerr(datasize))
  CALL brb_allocate (buf, Lmax, datasize)

  !%%%%%%%%%%%%%%%%%%%% 3. Mock simulation and brb %%%%%%%%%%%%%%%%%%%%%%%%%
  ! theta=0 (completely correlated) --> pi/2 (independent)
  ! correlation between x_n and x_{n+m} is costheta**m = exp(-m/xi)  roghly
  ! hence correlation length is xi=-1/LOG(costheta),
  ! costheta=EXP(-1/xi)
  ! NOTE: correlation length is at which correlation decays to 37 percent ...
  ! YES, IT's CORRECT!"
  !theta = PI*0.35_DP
  !costheta = COS(theta)
  !sintheta = SIN(theta)
  corrlen = 1.0_DP
  costheta = EXP(-1/corrlen)
  sintheta = SQRT(1-costheta**2)
  corrlen2 = 300.0_DP
  costheta2 = EXP(-1/corrlen2)
  sintheta2 = SQRT(1-costheta2**2)

  xx(1) = G05DDF(0d0, 1d0)
  xx(2) = G05DDF(0d0, 1d0) + 10.0_DP
  DO n=0,Nmax-1
     xx(1) = costheta*xx(1) + sintheta*G05DDF(0d0, 1d0)
     xx(2) = costheta2*(xx(2) - 10.0_DP) + sintheta2*G05DDF(0d0, 1d0) + 10.0_DP
     WRITE (10,"(I10,2F24.16)") n, xx(1), xx(2)
     CALL brb_append (buf, xx)
  END DO

  !%%%%%%%%%%%%%%%%%%%% 4. Report %%%%%%%%%%%%%%%%%%%%%%%%%
  CALL brb_report (buf, xxmean, xxerr)
  PRINT *
  PRINT *,xxmean, xxerr

END PROGRAM hello


!%%%% I think you can plot the output with %%%%
! xmgrace -settype xydy energy.out