! 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, obs and f are variables of "bird", "square"-pairs.     
!
! f is a functional parameter f= SUM(& x) * alpha * beta
! x depends stochastically on f,

! 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 input
  USE output
  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_real) :: f
  TYPE(variables_int) :: x
  TYPE(variables_int) :: obs


  ! the spatial structure
  
  INTEGER, DIMENSION(300, 5) :: spatial
  ! HARD-CODED:
  ! first dimension: the number of spatial units (squares)
  ! second dimension: the number of maximum neighbours + 1
  
  ! the spatial(i, 1) defines how many neighbours square i has and 
  ! spatial(i, 2) ... spatial(i, n) define their indices.


  
  ! 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

  
  ! HARD-CODED: If random update strategy, we need these
  ! HARD-CODED: The number of parameters to update
  INTEGER :: parameters = 200 + 300 + 200*300 + 10
  INTEGER :: parameters_achieved = 0


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

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

  ALLOCATE(f % two_dim(1:300, 1:200))
  ! HARD-CODED: we don't need the buffer for a functional parameter
  
  ALLOCATE(x % two_dim(1:300, 1:200))
  ALLOCATE(x % buffer(SIZE(x % one_dim) * multiplier))  

  ALLOCATE(obs % two_dim(1:300 , 1:200))
  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 % successful_changes = 0
  alpha % buffer_index = -1
  
  beta  % one_dim % successful_changes = 0
  beta  % buffer_index = -1
  
  x     % two_dim % successful_changes = 0
  x     % buffer_index = -1
  
  obs   % two_dim % successful_changes = 0
  obs   % buffer_index = -1
  
  ! HARD-CODED: if we have random update strategy...
  alpha % one_dim % update_count = 0
  beta  % one_dim % update_count = 0
  x     % two_dim % update_count = 0
  obs   % two_dim % update_count = 0

  ! HARD-CODED: if we have random update strategy...
  ! here are the updates wanted. they are achieved from the input
  alpha % onde_dim % updates_wanted = 100
  alpha % onde_dim % updates_wanted = 100
  x     % two_dim % updates_wanted = 200
  obs     % two_dim % updates_wanted = 3
  ! if the user wants to be more specific, it's all hard-coded here
  ! NOTE: these values must be > 0 (if not, there will be an infinite loop)

  ! setting the initial values and reading data.
  CALL read_data(alpha, beta, x, obs)
  CALL set_initial_values(alpha, beta, x, obs)
  
  CALL set_spatial(spatial)

  CALL set_functional_parameters(alpha, beta, f, x, obs)

  CALL check_initial_values(alpha, beta, x, obs)


  OPEN(UNIT=15, FILE=output_filename)


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


  ! HARD-CODED: if random update, we would use this

 ! DO WHILE (parameters_achieved < parameters)
    
 !     CALL update_one(alpha, beta, f, x, obs, parameters_achieved)
 !     IF (i > burn_in .AND. MOD(i, thinning) == 0) THEN
 !         CALL write_output(alpha, beta, f, x, obs)
 !     END IF
 ! END DO
  
  CALL write_summary(alpha, beta, x, obs)


  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 set_functional_parameters(alpha, beta, f, x, obs, spatial)
    IMPLICIT NONE
    
    TYPE(variables_real), INTENT(IN) :: alpha
    TYPE(variables_int), INTENT(IN) :: beta
    TYPE(variables_real), INTENT(INOUT) :: f
    TYPE(variables_int), INTENT(IN) :: x
    TYPE(variables_int), INTENT(IN) :: obs
    
    INTEGER, DIMENSION(300, 5), INTENT(IN) :: spatial

    INTEGER :: i, j, k
    INTEGER :: sum_x

    ! HARD-CODED: the loop lengths
    DO i=1, 300 ! i= current square
        DO j=1, 200
            ! HARD-CODED: the equation which defines the functional dependency
            ! how to do this: loop through neighbours... ???
            f % two_dim(i, j) % value = alpha % one_dim(i) % value * beta % one_dim(j) % value
            
            ! if we would like to loop through neighbours
            sum_x=0
            DO k=2, (spatial(i, 1)+1) ! all adjacent squares
                sum_x = sum_x + x % two_dim(spatial(i, k), j) % value
            END DO

            ! if we would like to divide by the number of the neighbours...
            !sum_x/spatial(i, 1)

        END DO
    END DO
  END SUBROUTINE set_functional_parameters
 
SUBROUTINE update_all(alpha, beta, f, 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



    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 = 300 
    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, 200
          ! 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

          ! HARD-CODED: there is a functional parameter which needs to be updated
          DO j=1, 200
              f % two_dim(i, j) = alpha % one_dim(i) % value * beta % one_dim(j) % value
          END DO

       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


    limit = 200 ! 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)

          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.




       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 * 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, 300          
          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')


       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


    limit = 300
    limit2 = 200
    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
          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.

          ! prior dist. of x
          CALL uniform_discrete_frequency(NINT(alpha % one_dim(i) % value), &
               beta % one_dim(j) % value, current_int % new_value, new_frequency)
          
          CALL uniform_discrete_frequency(NINT(alpha % one_dim(i) % value), &
               beta % one_dim(j) % value, current_int % value, 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.



          ! prior dist. of obs
          CALL binomial_frequency(current_int % new_value, 0.5, obs % two_dim(i, j) % value, new_frequency)
          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



          ! 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
    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 update_one(alpha, beta, x, obs, parameters_achieved)
    ! This subroutine updates all the parameters once.
    
    ! assumed: 
    ! initial values have been set, data read into memory 

    ! note: now we use the variables' updates_wanted counter

    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, INTENT(INOUT) :: parameters_achieved

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

    INTEGER :: int_transition
    REAL :: real_transition

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

    CALL RANDOM_NUMBER(prob_temp)
    ! HARD-CODED: we have 300 alphas, 200 betas and 200*300 x's and 10 obs's 
    which_to_update = prob_temp * 60510 +1

    ! now which_to_update should be an integer between 1 and 60510

    ! HARD-CODED: we have 300 alphas
    IF (which_to_update <= 300) THEN
        name='alpha'
        current_real=alpha % one_dim(which_to_update)
    ! HARD-CODED: we have 200 betas
    ELSE IF (which_to_update <= 500) THEN
        name='beta'
        current_int=beta % one_dim(which_to_update)
    ELSE IF (which_to_update <= 60500) THEN
        name='x'
        i = MOD(which_to_update-500, 300)
        j = (which_to_update-500)/300
        current_int = x % two_dim(i, j)
    ELSE
        name='obs'
        ! note: here some obs come from the data! so we can't take the correct obs so easily
        current_int = obs % two_dim(obs % two_dim_missing(which_to_update-60500, 1), &
            obs % two_dim_missing(which_to_update-60500, 2))
    END IF

    ! NOTE: this is similar to the update_all routine, except:
    ! 1) the current_int or current_real is already set
    ! 2) if a parameter achieves its updates_wanted-limit, we update the parameters_achieved

    SELECT CASE (name)

    CASE ('alpha')
       
       ! NOTE
       current_real % update_count = current_real % update_count +1
       IF (current_real % update_count == current_real % updates_wanted) THEN
          parameters_achieved = parameters_achieved +1
       END IF
       
       
       p_acc = 1

       ! 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, 200
          ! 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
       
    
    ! The variables to be updated with Gibbs will go through 
    ! a different update sequence...
    ! 
    

   

    ! HARD-CODED: same things as in case alpha

    CASE ('beta')

      ! NOTE
      current_int % update_count = current_int % update_count +1
      IF (current_int % update_count == current_int % updates_wanted) THEN
          parameters_achieved = parameters_achieved +1
      END IF
       
       p_acc = 1

       ! 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)

          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.




       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 * 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, 300          
          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')


       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
       
    CASE ('x')
        
        current_int % update_count = current_int % update_count +1
        IF (current_int % update_count == current_int % updates_wanted) THEN
            parameters_achieved = parameters_achieved +1
        END IF



        p_acc = 1
        
          
        ! 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
        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.

        ! prior dist. of x
        CALL uniform_discrete_frequency(NINT(alpha % one_dim(i) % value), &
             beta % one_dim(j) % value, current_int % new_value, new_frequency)
          
        CALL uniform_discrete_frequency(NINT(alpha % one_dim(i) % value), &
             beta % one_dim(j) % value, current_int % value, 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.



        ! prior dist. of obs
        CALL binomial_frequency(current_int % new_value, 0.5, obs % two_dim(i, j) % value, new_frequency)
        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



        ! 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
          

    
    CASE ('obs')     
        
        current_int % update_count = current_int % update_count +1
        IF (current_int % update_count == current_int % updates_wanted) THEN
            parameters_achieved = parameters_achieved +1
        END IF
       
        p_acc = 1
        ! 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 SELECT
              
  END SUBROUTINE update_one

  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
