! Prototype program for sequential
! update of variables, without blocking
! with a single total iterations option
! of a simple model, where

! alpha is a parameter of the repetitive structure "square"
! beta  is a parameter of "bird"
! x and obs are variables of "bird", "square"-pairs.     
!
! x depends stochastically on alpha and beta,
! obs depends stochastically on x.
! obs is found in data.

! alpha conforms to the uniform continuous distribution from 0 to 1
!
!

! TODO: proposal distribution subroutines
! TODO: distribution syntax must follow that of NAG

!
! reads: initial values, data, adjacency matrix
!

PROGRAM main
  USE definitions
  USE proposal
  IMPLICIT NONE 
  
  ! HARD-CODED: the variables in the model and their type (real or integer)
  TYPE(variables_real) :: alpha
  TYPE(variables_int) :: beta
  
  TYPE(variables_int) :: x
  TYPE(variables_int) :: obs
  
  ! HARD-CODED: the output file name
  CHARACTER(LEN=*), PARAMETER :: output_filename='output.txt'

  INTEGER :: i
  INTEGER, PARAMETER :: multiplier=100 
  ! this is the length multiplier 
  ! for a proposal buffer, for example if alpha has 300 values, then 
  ! the buffer for alpha would be 100 * 300 = 30 000 in size.

  ! HARD-CODED: The thinning factor and the burn-in iteration count
  INTEGER, PARAMETER :: thinning=100, burn_in=100


  CALL random_init()
  
  ! HARD-CODED: the dimension (3 squares)
  ALLOCATE(alpha % one_dim(1: 3))
  ALLOCATE(alpha % buffer(SIZE(alpha % one_dim) * multiplier))

  ! HARD-CODED: the dimension (2 birds)
  ALLOCATE(beta % one_dim(1: 2))
  ALLOCATE(beta % buffer(SIZE(beta % one_dim) * multiplier))
  
  ! HARD-CODED: the dimension (birds-horizontal, squares-vertical)
  ALLOCATE(x % two_dim(1:3, 1:2))
  ALLOCATE(x % buffer(SIZE(x % one_dim) * multiplier))

  ALLOCATE(obs % two_dim(1:3 , 1:2))
  ALLOCATE(obs % buffer(SIZE(obs % one_dim) * multiplier))
  
  ! HARD-CODED: 10 values are missing from the data (in obs)
  ! ALLOCATE(obs % two_dim_missing(1:10, 2))

  ! HARD-CODED: obs is the only variable from data, 
  ! so we don't allocate other 'missing'-arrays here.

  ! then set all updates to zero.
  ! wanted updates are not set, this time the user has requested a total
  ! update number to be done.

  alpha % one_dim % update_count = 0
  alpha % one_dim % successful_changes = 0
  alpha % buffer_index = -1

  beta  % one_dim % update_count = 0
  beta  % one_dim % successful_changes = 0
  beta  % buffer_index = -1
  
  x     % two_dim % update_count = 0
  x     % two_dim % successful_changes = 0
  x     % buffer_index = -1
  
  obs   % two_dim % update_count = 0
  obs   % two_dim % successful_changes = 0
  obs   % buffer_index = -1
  

  ! setting the initial values and reading data.
  CALL read_data(alpha, beta, x, obs)
  CALL set_initial_values(alpha, beta, x, obs)
  WRITE(*, *) 'Initial values read'
  CALL check_initial_values(alpha, beta, x, obs)
  WRITE(*, *) 'Initial values checked'

  OPEN(UNIT=15, FILE=output_filename)
  WRITE(*, *) 'Output opened'

  ! HARD-CODED: the number of iterations (here 1000)
  DO i=1, 1000
     ! HARD-CODED: the update strategy (here: sequential)
     CALL update_all(alpha, beta, x, obs)
     WRITE(*, *) 'updated the ', i, 'th time'
     IF (i > burn_in .AND. MOD(i, thinning) == 0) THEN
        CALL write_output(alpha, beta, x, obs)
     END IF
  END DO

  WRITE(*, *) 'Finished updating'
  CALL write_summary(alpha, beta, x, obs)
  WRITE(*, *) 'end of program reached'

  CLOSE(15)

CONTAINS    
  SUBROUTINE random_init()
    ! This subroutine initializes the default random number generator.
    IMPLICIT NONE
    INTEGER :: size, status
    INTEGER, DIMENSION(8) :: time
    INTEGER, DIMENSION(:), ALLOCATABLE :: seed
    CALL RANDOM_SEED(size=size)
    ALLOCATE(seed(size), stat=status)
    IF (status > 0) STOP 'Failed to allocate!'
    CALL DATE_AND_TIME(values = time)
    seed = 100*time(7) + time(8)/10
    CALL RANDOM_SEED(put=seed)
  END SUBROUTINE random_init

  SUBROUTINE read_data(alpha, beta, x, obs)
    ! This subroutine reads the data.
    IMPLICIT NONE

    TYPE(variables_real), INTENT(INOUT) :: alpha
    TYPE(variables_int), INTENT(INOUT) :: beta

    TYPE(variables_int), INTENT(INOUT) :: x
    TYPE(variables_int), INTENT(INOUT) :: obs    

    INTEGER :: i, j, temp, next
    ! INTEGER :: space_index
    TYPE(variable_int) :: current
    ! HARD-CODED: the maximum line length that occurs in the data file
    !    CHARACTER(LEN=3000) :: buffer 
    
    !   CHARACTER(LEN=3000) :: temp
    ! HARD-CODED: the data file name(s)
    ! these have been generated: the observations_1 file has zeros on the
    ! places where the original had *-signs. This makes reading it easier.
    ! the observations_missing has the places of the missing variables
    ! in a matrix of boolean values (0, 1).
    CHARACTER(LEN=*), PARAMETER :: data_filename='observations_1.txt'
    CHARACTER(LEN=*), PARAMETER :: data_missing='observations_missing.txt'
    
    OPEN(UNIT=16, FILE=data_filename)
    OPEN(UNIT=17, FILE=data_missing)
    
    next = 1
    
    ! these loop lengths come from the generator. 
    ! These are the dimensions of the data array (obs), and should match
    ! the dimensions of the matrix in the file.

    ! HARD-CODED: the loop length
    DO i = 1, 3

       READ(16, *, ERR=13, END=17) obs % two_dim(i, :) % value

       obs % two_dim(i, :) % is_data = .TRUE.

       DO j=1, 2

          READ(17, FMT='(I2)', ADVANCE='NO', ERR=13, END=17) temp
          IF (temp == 1) THEN
             obs % two_dim(i, j) % is_data = .FALSE.
             obs % two_dim_missing(next, 1) = i
             obs % two_dim_missing(next, 2) = j
             next = next + 1
          END IF
       END DO

       !       READ(16, FMT='(A)', ERR=13) buffer
       ! HARD-CODED: the loop length
       !      DO j=1, 200
       !        space_index = INDEX(buffer, ' ')
       !       temp = TRIM(buffer(1:space_index))
       !      buffer = buffer(space_index + 1:LEN_TRIM(buffer))
       !     
       !    IF (temp == '*') THEN
       !      ! HARD-CODED: we have a variable called obs and it is two-dimensional
       !     obs % two_dim(i, j) % is_data = .FALSE.
       !    obs % two_dim(i, j) % initial_value_set = .FALSE.
       !   CYCLE
       !  END IF
       !  READ (temp, *, ERR=17) obs % two_dim(i, j) % value
       !  obs % two_dim(i, j) % is_data = .TRUE.        
       
       !       END DO
    END DO
17  CLOSE(16)
    CLOSE(17)
    RETURN
13  WRITE(*, *) 'Could not read ', data_filename
    STOP
  END SUBROUTINE read_data

  SUBROUTINE set_initial_values(alpha, beta, x, obs)
    ! This subroutine reads the initial values of the parameters and sets them.
    IMPLICIT NONE

    TYPE(variables_real), INTENT(INOUT) :: alpha
    TYPE(variables_int), INTENT(INOUT) :: beta

    TYPE(variables_int), INTENT(INOUT) :: x
    TYPE(variables_int), INTENT(INOUT) :: obs

    CHARACTER(LEN=30) :: line, var_name
    ! HARD-CODED: the filename of the initial values file
    CHARACTER(LEN=*), PARAMETER :: initial_values_filename='initial_values.txt'
    INTEGER :: i, j, linenum
    INTEGER :: vertical_start, vertical_end, horizontal_start, horizontal_end
    INTEGER :: ind

    ! The format of the file is specified in the Design document 
    ! (section Simulation parameters input files)

    OPEN(UNIT=16, FILE=initial_values_filename)
    WRITE(*, *) 'Reading ', initial_values_filename

    READ (16, FMT='(A)', ERR=14, END=12) line
    linenum = 1

    IF (.NOT. line == '?? initial values') THEN
       WRITE(*, *) 'Expected initial values legend in ' , initial_values_filename
       WRITE(*, *) 'on line 1 not found. Quitting.'
       CLOSE(1)
       STOP
    END IF
    
    DO WHILE (.TRUE.)       
       READ (16, FMT='(A)', ERR=14, END=12) line
       linenum = linenum + 1

       line = TRIM(line)
   !    PRINT(*, *) line
       
       IF (LEN_TRIM(line) == 0) CYCLE
       IF (line(1:1) == '#') CYCLE
       
       ! now we have a "? alpha 1:something" at hand
       IF(line(1:1) /= '?')  THEN
          WRITE(*, *) 'Expected legend on line ', linenum, &
               ' not found, quitting.'
          STOP
       END IF
       line = line (3:LEN_TRIM(line))
       ind = INDEX(line, ' ')
       IF (ind <= 1) THEN 
          var_name = line
       ELSE
          var_name = TRIM(line(1:ind))
       END IF
!       WRITE(*, *) "variable name: ", var_name
       line = line(LEN_TRIM(var_name)+2:LEN_TRIM(line)) 
       ! trim the variable name out 
       ! and the trailing space
       
       vertical_start  = 0
       vertical_end    = 0
       horizontal_start = 0
       horizontal_end   = 0
       
       SELECT CASE (var_name)
          ! HARD-CODED: the name of the variable
       CASE ('alpha') 
          WRITE(*, *) 'Case alpha reached'
          ind = INDEX(line, ':')
          line(ind:ind) = ' '

          ! HARD-CODED: alpha is one-dimensional so we consider the format
          ! "alpha 4:6" for example
          
          WRITE(*, *) line
          READ (line, FMT=*, ERR=14, END=12) vertical_start, vertical_end
          WRITE(*, *) 'Vertical indices are: ', vertical_start, vertical_end
          DO i = vertical_start, vertical_end
             READ(16, *, ERR=14, END=12) alpha % one_dim (i) % value
             alpha % one_dim(i) % initial_value_set = .TRUE.
          END DO
          ! end of case alpha
          
          ! HARD-CODED: same things as in case alpha
       CASE ('beta') 
          ind = INDEX(line, ':')
          line(ind:ind) = ' '
          !WRITE(*, *) line
          READ (line, FMT=*, ERR=14, END=12) vertical_start, vertical_end
          !WRITE(*, *) 'Vertical indices are: ', vertical_start, vertical_end
          DO i = vertical_start, vertical_end
             READ(16, *, ERR=14, END=12) beta % one_dim (i) % value
             beta % one_dim(i) % initial_value_set = .TRUE.
          END DO
          ! end of case beta

          ! HARD-CODED: the name of the variable
       CASE ('x')

          ! HARD-CODED: x is two-dimensional so we consider the 
          ! format "x 1:20 5:400" for example
          DO i = 1, 2
             ind = INDEX(line, ':')
             line(ind:ind) = ' '
          END DO

          READ (line, FMT=*, ERR=14, END=12) vertical_start, vertical_end, horizontal_start, horizontal_end 
          DO i = vertical_start, vertical_end
             READ(16, *, ERR=14, END=12) x % two_dim (i, horizontal_start:horizontal_end) % value
             x % two_dim(i, horizontal_start:horizontal_end) % initial_value_set = .TRUE.
             
          END DO
          ! end of case x

          ! HARD-CODED: same things as in case x
       CASE ('obs')
          DO i = 1, 2
             ind = INDEX(line, ':')
             line(ind:ind) = ' '
          END DO
          
          READ (line, FMT=*, ERR=14, END=12) vertical_start, vertical_end, horizontal_start, horizontal_end 
          DO i = vertical_start, vertical_end
             READ(16, *, ERR=14, END=12) obs % two_dim (i, horizontal_start:horizontal_end) % value
             obs % two_dim(i, horizontal_start:horizontal_end) % initial_value_set = .TRUE.             
          END DO
          ! end of case obs

       CASE DEFAULT
       END SELECT

    END DO ! Ends the file-reading loop.
12  CLOSE(16)

    ! test printout
    !    DO i = vertical_start, vertical_end
    !       DO j = horizontal_start, horizontal_end
    !          WRITE (*, *) x % two_dim (i, j) % value
    !       END DO
    !    END DO

    RETURN
14  WRITE(*, *) 'Error on line ', linenum
    CLOSE(16)
    STOP
  END SUBROUTINE set_initial_values


  SUBROUTINE check_initial_values(alpha, beta, x, obs)
    IMPLICIT NONE

    TYPE(variables_real), INTENT(INOUT) :: alpha
    TYPE(variables_int), INTENT(INOUT) :: beta
    
    TYPE(variables_int), INTENT(INOUT) :: x
    TYPE(variables_int), INTENT(INOUT) :: obs    

    ! We must now ensure that all the variables not in data
    ! have an initial value.
    ! That is, for every array that holds data or variables
    ! the following condition must hold:
    ! array(:[,:]) % initial_value_set .OR. array(:[,:]) % is_data
   
    ! in this hardcoded case, every array = {alpha, beta, x, obs}.
    

    ! HARD-CODED: we have a variable group called alpha and it's one-dimensional
    IF(ANY((.NOT. alpha % one_dim(:) % initial_value_set .AND. .NOT. alpha % one_dim(:) % is_data))) THEN 
       WRITE(*, *) 'Insufficient initial values for group alpha'
       WRITE(*, *) 'Simulation cannot start.'
       STOP
    END IF

    ! HARD-CODED: we have a variable group called beta and it's one-dimensional
    IF(ANY((.NOT. beta % one_dim(:) % initial_value_set .AND. .NOT. beta % one_dim(:) % is_data))) THEN 
       WRITE(*, *) 'Insufficient initial values for variable group beta'
       WRITE(*, *) 'Simulation cannot start.'
       STOP
    END IF
    
    ! HARD-CODED: we have a variable group called x and it's two-dimensional
    IF(ANY((.NOT. x % two_dim(:, :) % initial_value_set .AND. .NOT. x % two_dim(:, :) % is_data))) THEN 
       WRITE(*, *) 'Insufficient initial values for variable group x'
       WRITE(*, *) 'Simulation cannot start.'
       STOP
    END IF

    ! HARD-CODED: we have a variable group called obs and it's two-dimensional
    IF(ANY((.NOT. obs % two_dim(:, :) % initial_value_set .AND. .NOT. obs % two_dim(:, :) % is_data))) THEN 
       WRITE(*, *) 'Insufficient initial values for variable group obs'
       WRITE(*, *) 'Simulation cannot start.'
       STOP
    END IF

  END SUBROUTINE check_initial_values
    
SUBROUTINE update_all(alpha, beta, x, obs)
    ! This subroutine updates all the parameters once.
    
    ! assumed: 
    ! initial values have been set, data read into memory 

    IMPLICIT NONE


    TYPE(variables_real), INTENT(INOUT) :: alpha
    TYPE(variables_int), INTENT(INOUT) :: beta
    
    TYPE(variables_int), INTENT(INOUT) :: x
    TYPE(variables_int), INTENT(INOUT) :: obs    

    INTEGER :: i, j, k, limit, limit2, temp
    REAL :: p_acc, new_frequency, frequency
    TYPE(variable_int) :: current_int
    TYPE(variable_real) :: current_real

    !Prototype, no need for random walk temporary variables since random walk
    !is not used yet.
    !INTEGER :: int_transition
    !REAL :: real_transition

    REAL, SAVE, DIMENSION(100000) :: prob
    INTEGER, SAVE :: prob_index = -1
    
    IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
       CALL RANDOM_NUMBER(prob)
       prob_index = 1
    END IF

    ! looping through alpha...

    ! HARD-CODED: we have 300 alphas
    limit = 3
    DO i = 1, limit
       p_acc = 1

       ! HARD-CODED: we have a variable group called alpha and it's one-dimensional
       current_real = alpha % one_dim(i)

       ! if we need to generate more proposals:

       ! HARD-CODED: the name of the variable group ("alpha")
       IF (alpha % buffer_index > SIZE(alpha % buffer) .OR. alpha % buffer_index < 1) THEN
          CALL generate('alpha', alpha % buffer)
          alpha % buffer_index = 1
       END IF

       ! this is new value generation on a fixed proposal strategy, 
       ! the generator has chosen this according to instructions

       ! HARD CODED: alpha's proposal strategy (fixed or rw, fixed chosen)
       current_real % new_value = alpha % buffer(alpha % buffer_index)

       ! If this was random walk instead:
       !current_real % new_value = current_real % value +  & 
       !     alpha % buffer(alpha % buffer_index)

       alpha % buffer_index = alpha % buffer_index + 1

       ! now alpha has a new value.
 
       ! here we have hardcoded the formula for alpha's likelihood.

       ! HARD-CODED: how to calculate alpha's accepting probability
       
       ! this is: p(alpha')
       CALL uniform_continuous_density(0.0, 1.0, current_real % new_value, new_frequency)
       ! this is: p(alpha)
       CALL uniform_continuous_density(0.0, 1.0, current_real % value, frequency)

       p_acc = p_acc * new_frequency/frequency
       ! this density function is assumed to first take the parameters 
       ! and then the point for which the density is computed.


       ! HARD-CODED: x depends on alpha and beta
       ! HARD-CODED: the loop length
       DO j = 1, 2
          ! HARD-CODED: this is how x depends on alpha and beta

          ! this is p(x | alpha', beta)
          CALL uniform_discrete_frequency(NINT(current_real % new_value), &
               beta % one_dim(j) % value, & 
               x % two_dim(i, j) % value, new_frequency)

          ! this is p(x | alpha, beta)
          CALL uniform_discrete_frequency(NINT(current_real % value), &
               beta % one_dim(j) % value, & 
               x % two_dim(i, j) % value, frequency)


          p_acc = p_acc * new_frequency / frequency

       END DO

       ! The q(alpha', alpha) / q(alpha, alpha')
       ! if fixed, this is P(alpha) / P(alpha')

       ! HARD-CODED: the proposal distribution of alpha
       CALL uniform_continuous_density(0.0, 1.0, current_real % new_value, new_frequency)
       CALL uniform_continuous_density(0.0, 1.0, current_real % value, frequency)


       p_acc = p_acc * frequency/new_frequency

       ! If random walk instead:
       ! real_transition = current_real % new_value - current_real % value
       ! p_acc = p_acc * normal_density(0, 1, -real_transition)/ &
       !     normal_density(0, 1, real_transition)
              
       ! this is the M-H acceptation...
       IF(p_acc > prob(prob_index)) THEN
          current_real % value = current_real % new_value 
          current_real % successful_changes = current_real % successful_changes + 1
       END IF
       
       
       prob_index = prob_index + 1
       
       
       IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
          CALL RANDOM_NUMBER(prob)
          prob_index = 1
       END IF
       
    END DO ! End of update all alpha.
    
    ! The variables to be updated with Gibbs will go through 
    ! a different update sequence...
    ! 
    
    ! HARD-CODED: same things as in case alpha
    WRITE(*, *) 'alphas updated'

    limit = 2 ! the number of betas
    DO i = 1, limit
       ! update all betas
       p_acc = 1
       current_int = beta % one_dim(i)

       ! if we need to generate more proposals:
       IF (beta % buffer_index > SIZE(beta % buffer) .OR. beta % buffer_index < 1) THEN
          CALL generate('beta', beta % buffer)
          WRITE(*, *) 'generate called on beta'
          beta % buffer_index = 1
       END IF

       ! this is new value generation on a fixed proposal strategy, 
       ! the generator has chosen this according to instructions
       current_int % new_value = beta % buffer(beta % buffer_index)

       ! If this was random walk instead:
       !current_int % new_value = current_int % value +  & 
       !     beta % buffer(beta % buffer_index)

       beta % buffer_index = beta % buffer_index + 1

       ! now beta has a new value.

       ! here we have hardcoded the formula for beta's likelihood.
       WRITE(*, *) 'about to call uniform disrete with values', 1, 2, current_int % new_value
       WRITE(*, *) 'about to call uniform discrete with values', 1, 2, current_int % value


       CALL uniform_discrete_frequency(1, 2, current_int % new_value, new_frequency)
       CALL uniform_discrete_frequency(1, 2, current_int % value, frequency)

       WRITE(*, *) new_frequency, frequency
       p_acc = p_acc * new_frequency / frequency 
       ! this density function is assumed to first take the parameters 
       ! and then the point for which the density is computed.
       
       DO j = 1, 3          
          CALL uniform_discrete_frequency(NINT(alpha % one_dim(j) % value), &
               current_int % new_value, & 
               x % two_dim(i, j) % value, new_frequency)

          CALL uniform_discrete_frequency(NINT(alpha % one_dim(j) % value), &
               current_int % value, & 
               x % two_dim(i, j) % value, frequency)
          

          p_acc = p_acc * new_frequency / frequency
          
       END DO

       ! The q(b', b) / q(b, b')
       ! if fixed, this is P(b) / P(b')
       WRITE(*, *) current_int % new_value

       CALL uniform_discrete_frequency(1, 2, current_int % new_value, new_frequency)
       CALL uniform_discrete_frequency(1, 2, current_int % value, frequency)

       p_acc = p_acc * frequency/new_frequency

       ! If random walk instead:
       ! real_transition = current_int % new_value - current_int % value
       ! p_acc = p_acc * normal_density(0, 1, -real_transition)/ &
       !     normal_density(0, 1, real_transition)
       
       ! this is the M-H acceptation...
       IF(p_acc > prob(prob_index)) THEN
          current_int % value = current_int % new_value 
          current_int % successful_changes = current_int % successful_changes + 1
       END IF
       
       prob_index = prob_index + 1

       ! and then check for remaining randoms.
       
       IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
          CALL RANDOM_NUMBER(prob)
          prob_index = 1
       END IF
       
    END DO

    WRITE(*, *) 'betas updated'

    limit = 3
    limit2 = 2
    DO i = 1, limit
       DO j = 1, limit2
          ! update all x
          p_acc = 1
          current_int = x % two_dim(i, j)
          
          ! we need to generate more proposals
          IF (x % buffer_index > SIZE(x % buffer) .OR. x % buffer_index < 1) THEN
             CALL generate('x', x % buffer)
             x % buffer_index = 1
             WRITE(*, *) 'called generate for x'
          END IF
          ! this is new value generation on a fixed proposal strategy, 
          ! the generator has chosen this according to instructions
          current_int % new_value = x % buffer(x % buffer_index)
          x % buffer_index = x % buffer_index + 1
          
          ! If this was random walk instead:
          !current_int % new_value = current_int % value +  & 
          !     x % buffer(x % buffer_index)
          !x % buffer_index = x % buffer_index + 1
          
          ! now x has a new value.
          
          ! here we have hardcoded the formula for x's likelihood.
          WRITE(*, *) 'calling uniform discrete with values ', & 
               NINT(alpha % one_dim(i) % value), beta % one_dim(j) % value, &
               current_int % new_value

          CALL uniform_discrete_frequency(NINT(alpha % one_dim(i) % value), &
               beta % one_dim(j) % value, current_int % new_value, new_frequency)
          
          WRITE(*, *) 'calling uniform discrete with values ', & 
               NINT(alpha % one_dim(i) % value), beta % one_dim(j) % value, &
               current_int % value

          CALL uniform_discrete_frequency(NINT(alpha % one_dim(i) % value), &
               beta % one_dim(j) % value, current_int % value, frequency)
               
          WRITE(*, *) 'about to compare x''s frequencies'
          WRITE(*, *) 'dividing ', new_frequency, ' by ', frequency
     
          p_acc = p_acc * new_frequency / frequency          
          
          ! this density subroutine is assumed to first take the parameters 
          ! and then the point for which the density is computed, and 
          ! lastly the OUT-variable.


          WRITE(*, *) 'calling binomial freq with values ', current_int % new_value, 0.5, obs % two_dim(i, j) % value

          ! prior dist. of obs
          CALL binomial_frequency(current_int % new_value, 0.5, obs % two_dim(i, j) % value, new_frequency)

          WRITE(*, *) 'calling binomial freq with values ', current_int % value, 0.5, obs % two_dim(i, j) % value

          CALL binomial_frequency(current_int % value, 0.5, obs % two_dim(i, j) % value, frequency)
          
          p_acc = p_acc * new_frequency / frequency
          
          ! proposal dist. of x
          CALL binomial_frequency(3, 0.5, current_int % value, frequency)
          CALL binomial_frequency(3, 0.5, current_int % new_value, new_frequency)
          
          ! The q(a', a) / q(a, a')
          ! if fixed, this is P(a) / P(a')
          p_acc = p_acc * new_frequency/frequency
          
          WRITE(*, *) 'considering ', p_acc, ' against ', prob(prob_index)
          
          ! this is the M-H acceptation...
          IF(p_acc > prob(prob_index)) THEN
             WRITE(*, *) 'accepted'
             current_int % value = current_int % new_value 
             current_int % successful_changes = current_int % successful_changes + 1
          END IF
          
          prob_index = prob_index + 1
          
          ! and then check for remaining randoms.
          
          IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
             CALL RANDOM_NUMBER(prob)
             prob_index = 1
          END IF
          
       END DO
    END DO

    ! HARD-CODED: obs has 10 values missing in the data.
    ! now we need to update them
    !limit=10
    !DO k = 1, limit       
    !   i=obs % two_dim_missing(k, 1)
    !   j=obs % two_dim_missing(k, 2)
       ! update obs(i,j)
    !   p_acc = 1
    !   current_int = obs % two_dim(i, j)
       
       ! we need to generate more proposals
    !   IF (obs % buffer_index > SIZE(obs % buffer)) THEN
    !      CALL generate('obs', obs % buffer)
    !      obs % buffer_index = 1
    !   END IF
       ! this is new value generation on a fixed proposal strategy, 
       ! the generator has chosen this according to instructions
    !   current_int % new_value = obs % buffer(obs % buffer_index)
    !   obs % buffer_index = obs % buffer_index + 1
       

    !   CALL binomial_frequency(x % two_dim(i, j) % value, 0.5, current_int % new_value, new_frequency)
    !   CALL binomial_frequency(x % two_dim(i, j) % value, 0.5, obs % two_dim(i, j) % value, frequency)

    !   p_acc = p_acc * new_frequency/frequency
       
       ! and then q(obs', obs)/ q(obs, obs')
    !   CALL binomial_frequency(3, 0.5, current_int % value, frequency)
    !   CALL binomial_frequency(3, 0.5, current_int % new_value, new_frequency)
       
    !   p_acc = p_acc * frequency/new_frequency

    !   IF(p_acc > prob(prob_index)) THEN
    !      current_int % value = current_int % new_value 
    !      current_int % successful_changes = current_int % successful_changes + 1
    !   END IF
       
    !   prob_index = prob_index + 1
       
       ! and then check for remaining randoms.
       
    !   IF (prob_index > SIZE(prob) .OR. prob_index < 1) THEN
    !      CALL RANDOM_NUMBER(prob)
    !      prob_index = 1
    !   END IF
    !          
    !END DO
    ! update other parameters, if the model would have more    
    
  END SUBROUTINE update_all

  SUBROUTINE write_output(alpha, beta, x, obs)
    IMPLICIT NONE


    TYPE(variables_real), INTENT(IN) :: alpha
    TYPE(variables_int), INTENT(IN) :: beta
    
    TYPE(variables_int), INTENT(IN) :: x
    TYPE(variables_int), INTENT(IN) :: obs    

    ! output the parameters which the user wants to output
    ! (the parameters to output and the file name are hard-coded here)
    INTEGER :: i
    LOGICAL, SAVE :: written = .FALSE.
    
    IF(.NOT. written) THEN
       ! HARD-CODED: the parameters to write and their order
       WRITE(15, *) '# x_1_1 x_2_1 x_3_1 x_1_2 x_2_2 x_3_2 b_1 b_2'
       written = .TRUE.
    END IF
    
    ! HARD-CODED: which parameters to write
    WRITE(15, *) x % two_dim(:, :) % value, beta % one_dim(:) % value
        
  END SUBROUTINE write_output

  SUBROUTINE write_summary(alpha, beta, x, obs)
    IMPLICIT NONE

    TYPE(variables_real), INTENT(IN) :: alpha
    TYPE(variables_int), INTENT(IN) :: beta
    
    TYPE(variables_int), INTENT(IN) :: x
    TYPE(variables_int), INTENT(IN) :: obs    

    ! output the summary wanted
    ! HARD-CODED: the file name of the summary file

  END SUBROUTINE write_summary

  SUBROUTINE uniform_continuous_density(left, right, value, prob)
    REAL, INTENT(IN) :: value, left, right
    REAL, INTENT(OUT) :: prob
    
    prob = 1/(right - left)
    
  END SUBROUTINE uniform_continuous_density

  SUBROUTINE uniform_discrete_frequency(left, right, value, prob)
    INTEGER, INTENT(IN) :: value, left, right
    REAL, INTENT(OUT) :: prob
    
    prob = REAL(1 /((REAL(right) - REAL(left)) +1))

  END SUBROUTINE uniform_discrete_frequency

  SUBROUTINE binomial_frequency(n, p, value, prob)
    REAL, INTENT(IN) :: p
    REAL, INTENT(OUT) :: prob
    INTEGER, INTENT(IN) :: n, value

    IF (value >= n) THEN 
       prob = 0
       RETURN
    END IF

    ! a little help from the fortran book...
    prob = REAL(kert(n)) / REAL((kert(value)*kert(n - value)))
    prob = prob * p**value*(1-p)**(n - value)

  END SUBROUTINE binomial_frequency

  FUNCTION kert(n) RESULT(kert_arvo)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: n
    INTEGER :: kert_arvo

    INTEGER :: temp

    kert_arvo=1
    IF (n > 1) THEN
       temp = n
       DO WHILE (temp > 1)
          kert_arvo = temp * kert_arvo
          temp = temp - 1
       END DO
    END IF
  END FUNCTION kert
  
END PROGRAM main
