!###############################################################################
!                                                                              #
!         CCCCCC  SSSSSS  IIIII  RRRRRR   OOOO          L      IIIII           #
!         C       S         I    R    R  O    O         L        I             #
!         C       SSSSSS    I    RRRRRR  O    O  =====  L        I             #
!         C            S    I    R RR    O    O         L        I             #
!         CCCCCC  SSSSSS  IIIII  R   RR   OOOO          LLLLL  IIIII           #
!                                                                              #
!         Long Pocket Laboratories                                             #
!         120 Meiers Road                                                      #
!         Indooroopilly, QLD 4068                                              #
!         AUSTRALIA                                                            #
!         Ph. (07) 3214 2392                                                   #
!         Fx. (07) 3214 2881                                                   #
!                                                                              #
!         Objective:      Bayesian Mixture model with known N of components    #
!         Author:         Toni Reverter (Tony.Reverter-Gomez@csiro.au)         #
!         Last Revision:  January 2003                                         #
!                                                                              #
!###############################################################################

PROGRAM mixture

IMPLICIT NONE


INTEGER, PARAMETER :: minni = 5,           & ! Min No Rec in a given cluster
                      minclusters = 1,     &
                      maxclusters = 5,     &
                      chainlength = 12000, &
                      burning = 2000

REAL, PARAMETER :: minmix = 0.005, & ! Minimum mixing proportion
                   zero = 0.0, half = 0.5, one = 1.0, two = 2.0,   &
                   vsmall = TINY(1.0), vlarge = HUGE(1.0)

INTEGER :: n,          & ! No of records
           nclusters     ! from minclusters to maxclusters

INTEGER, ALLOCATABLE, DIMENSION(:) :: &
           ni,    & ! No records in each cluster
           id,    & ! Cluster No to which each record belongs
           itemp   

TYPE rec_element 
     INTEGER :: location
     REAL :: value
END TYPE rec_element
type(rec_element), ALLOCATABLE, DIMENSION(:) ::  &
        y           ! Data records

REAL, ALLOCATABLE, DIMENSION(:) :: &
        mix,      & ! Mixing proportions (add to 1.0)
        mn,       & ! Mean of each cluster
        var,      & ! Variance of each cluster
        avgmix,   & ! Avg Mixing proportions (post-burning)
        avgmn,    & ! Avg Mean of each cluster (post-burning)
        avgvar,   & ! Avg Variance of each cluster (post-burning)
        postprob    ! Posterior Probability

REAL :: cummprob,            & ! Cummulative Probability
        u,rtemp,prevvar

REAL :: loglikel, mean0, ss0, var0, df

REAL :: random_normal, &
        normal_eval,   &
        loglikel_eval, &
        random_chisq,  & 
        random_gamma,  & 
        random_gamma1, & 
        random_gamma2, & 
        random_exponential

INTEGER :: i,j,k,round,attemps,dt(8)
CHARACTER(LEN=10) :: date,time,zone

!###############################################################
OPEN(5,file='dat',form='formatted',status='old')
OPEN(50,file='postprob',form='formatted',status='unknown')
!###############################################################

CALL RANDOM_SEED()


!**********************
! Read input data
!**********************
n = 0
DO
   READ(5,*,IOSTAT=k)
   IF( k .NE. 0 ) EXIT
   n = n + 1
ENDDO
ALLOCATE(y(n),id(n))
REWIND(5)
DO i = 1, n
   y(i)%location = i
   READ(5,*)y(i)%value
ENDDO
CLOSE(5)

i=1
CALL sort_rec(y,i,n)

!Test that the sorting is OK
DO i = 1, n
   write(33,*)y(i)%location,y(i)%value
enddo
!stop

!**********************
! Start Big Loop
!**********************

DO nclusters = minclusters, maxclusters

   mean0 = SUM(y%value)/n
   ss0 = zero
   DO i = 1, n
      ss0 = ss0 + (y(i)%value-mean0)**2
   ENDDO
   var0 = ss0/(n-1)
   ss0 = one/var0

   ALLOCATE(ni(nclusters),mix(nclusters),mn(nclusters),var(nclusters), &
            avgmix(nclusters),avgmn(nclusters),avgvar(nclusters),      &
            postprob(nclusters),itemp(nclusters))

   ! randomly allocate records to each cluster
   ! then compute means and var within cluster

   ni = 0
   mix = one/nclusters
   mn = zero
   var = zero

   k = 0
   DO j = 1, nclusters-1
      DO i = 1, int(n/nclusters)
         k = k + 1
         id(k) = j
         ni(j) = ni(j) + 1
         mn(j) = mn(j) + y(k)%value
      ENDDO
   ENDDO
   DO i = k+1, n
      id(i) = nclusters
      ni(nclusters) = ni(nclusters) + 1
      mn(nclusters) = mn(nclusters) + y(i)%value
   ENDDO

   DO j = 1 , nclusters
   mn(j) = mn(j) / ni(j)
   ENDDO
   DO i = 1, n
      var(id(i)) = var(id(i)) + (y(i)%value-mn(id(i)))**2
   ENDDO
   DO j = 1 , nclusters
      var(j) = var(j) / (ni(j)-1)
   ENDDO

   ! Start Iterating

   avgmix = zero
   avgmn = zero
   avgvar = zero

   DO round = 1, chainlength

      ! Update to which cluster each record belongs to

      k = 0
      ni = 0
      DO i = 1, n
         rtemp = -one*vsmall
         DO j = 1, nclusters
            u = normal_eval(y(i)%value,mn(j),var(j))
            u = mix(j) * u
            IF( u > rtemp )THEN
               id(i) = j
               rtemp = u
            ENDIF
         ENDDO
         ni(id(i)) = ni(id(i)) + 1
      ENDDO
      DO j = 1 , nclusters
         IF( ni(j) < minni )THEN
            ni(j) = minni
            DO i = 1, ni(j)
               CALL RANDOM_NUMBER(u)
               id(INT(u*n)+1) = j
               id(i) = j
            ENDDO
         ENDIF
      ENDDO

      ! Update mean and variance within cluster

      mn = zero
      var = zero
      DO i = 1, n
         mn(id(i)) = mn(id(i)) + y(i)%value
      ENDDO
      DO j = 1 , nclusters
         mn(j) = mn(j) / ni(j)
      ENDDO
      DO i = 1, n
         var(id(i)) = var(id(i)) + (y(i)%value-mn(id(i)))**2
      ENDDO
      DO j = 1 , nclusters
         prevvar = var(j)
         var(j) = var(j) / (ni(j) - one)
         if(var(j)<0.001*prevvar)var(j)=0.1*prevvar
      ENDDO
   
      ! sample posterior for mn(j)

      DO j = 1, nclusters
         u = RANDOM_NORMAL()
         IF( var(j) > zero )THEN
            rtemp = one/(ni(j)/var(j) + nclusters)
            rtemp = rtemp * ( mn(j)*ni(j)/var(j) + nclusters*mean0 )
            mn(j) = rtemp + u*SQRT(one/(ni(j)/var(j) + nclusters))
         ELSE
            mn(j) = mn(j) + u*SQRT(one/(nclusters))
         ENDIF
      ENDDO
      
      ! sample posterior for var(j)

      DO j = 1, nclusters
         rtemp = zero
         DO i = 1, n
            IF( id(i) == j )rtemp = rtemp + (y(i)%value-mn(j))**2
         ENDDO
         rtemp = 2.0*ss0 + rtemp
         df = 2.0 + ni(j)
!        df = var0/(nclusters+2) + ni(j)
         u = RANDOM_CHISQ(df)
         var(j) = one/u * rtemp
      ENDDO

      ! sample posterior for mixing proprotions mix(j)
      DO j = 1, nclusters
        itemp(j) = 1 + ni(j)
      ENDDO
      CALL RANDOM_DIRICHLET(itemp,nclusters,mix)
      IF( MINVAL(mix) < minmix )THEN
         DO j = 1, nclusters
            IF( mix(j) < minmix ) mix(j) = minmix
         ENDDO
         DO j = 1, nclusters
            mix(j) = mix(j)/SUM(mix)
         ENDDO
      ENDIF
   
      IF( round > burning )THEN
         DO j = 1, nclusters
            avgmix(j) = avgmix(j) + mix(j)
            avgmn(j) = avgmn(j) + mn(j)
            avgvar(j) = avgvar(j) + var(j)
         ENDDO
      ENDIF

   ENDDO

   avgmix = avgmix/(chainlength-burning)
   avgmn  = avgmn/(chainlength-burning)
   avgvar = avgvar/(chainlength-burning)

   ! Compute posterior probabilities of each data point
   ! to belong to each of the clusters

   DO i = 1, n
      DO k = 1, n
         IF( y(k)%location == i )THEN
            postprob = zero
            cummprob = zero
            DO j = 1, nclusters
               u = normal_eval(y(k)%value,avgmn(j),avgvar(j))
               postprob(j) = avgmix(j) * u
               cummprob = cummprob + postprob(j)
            ENDDO
            DO j = 1, nclusters
               postprob(j) = postprob(j)/cummprob
            ENDDO
            WRITE(50,'(7f12.5)')y(k)%value,postprob
         ENDIF
      ENDDO
   ENDDO

   ! Compute Log Likelihood of this model

   loglikel = zero
   DO i = 1, n
      cummprob = zero
      DO j = 1, nclusters
         u = normal_eval(y(i)%value,avgmn(j),avgvar(j))
         cummprob = cummprob + avgmix(j)*u
      ENDDO
      loglikel = loglikel + log(cummprob)
   ENDDO

   ! Print out results for this model

   WRITE(*,*)
   WRITE(*,'(a15,i5)')'N Clusters',nclusters
   WRITE(*,'(a20)')REPEAT("=",20)
   WRITE(*,'(a15,7f12.5)')'Mixing Prop',avgmix
   WRITE(*,'(a15,7f12.5)')'Means',avgmn
   WRITE(*,'(a15,7f12.5)')'Variances',avgvar
   WRITE(*,*)
   WRITE(*,'(a15,1f15.3)')'LogL',loglikel
   WRITE(*,'(a15,1f15.3)')'AIC',-2*loglikel+2*(3*nclusters-1)
   WRITE(*,'(a15,1f15.3)')'BIC',-2*loglikel+(3*nclusters-1)*log(1.0*n)
   WRITE(*,*)
      
   DEALLOCATE(ni,mix,mn,var,avgmix,avgmn,avgvar,postprob,itemp)

ENDDO

END PROGRAM mixture
!###############################################################
!###############################################################
FUNCTION loglikel_eval(n,x,mean,var) RESULT(fn_val)
IMPLICIT NONE

INTEGER, INTENT(IN) :: n
REAL, INTENT(IN) :: x(n), mean, var
REAL :: fn_val, pi=3.141592654, ss
INTEGER :: i

ss = 0.0
DO i = 1, n
   ss = ss + (x(i)-mean)**2
ENDDO

fn_val = -n/2*log(2*pi)-1/2*log(var)-1/(2*var)*ss

END FUNCTION loglikel_eval
!###############################################################
FUNCTION normal_eval(x,mean,sigma2) RESULT(fn_val)
IMPLICIT NONE

REAL :: fn_val
REAL, INTENT(IN) :: x, mean, sigma2
REAL :: pi=3.141592654, sigma

sigma=sqrt(sigma2)

fn_val=(1/sqrt(2*pi*sigma2))*exp(-.5*(x-mean)**2/sigma2)

END FUNCTION normal_eval
!###############################################################
SUBROUTINE random_dirichlet(alpha,nsize,theta)
IMPLICIT NONE

INTEGER :: i, nsize
INTEGER,INTENT(IN)  :: alpha(nsize)
REAL,INTENT(OUT) :: theta(nsize)
REAL :: norm, random_gamma, zero = 0.0, half = 0.5, one = 1.0, two = 2.0


norm = zero
DO i = 1, nsize
   theta(i) = random_gamma(one*alpha(i))
   norm = norm +  theta(i)
ENDDO

DO i = 1, nsize
   theta(i) = theta(i)/norm
ENDDO

END SUBROUTINE random_dirichlet
!###############################################################
FUNCTION random_normal() RESULT(fn_val)
IMPLICIT NONE

! Adapted from the following Fortran 77 code
!      ALGORITHM 712, COLLECTED ALGORITHMS FROM ACM.
!      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
!      VOL. 18, NO. 4, DECEMBER, 1992, PP. 434-435.

REAL :: fn_val, half = 0.5

!     Local variables

REAL     :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472,           &
            r1 = 0.27597, r2 = 0.27846, u, v, x, y, q

!     Generate P = (u,v) uniform in rectangle enclosing acceptance region

DO
  CALL RANDOM_NUMBER(u)
  CALL RANDOM_NUMBER(v)
  v = 1.7156 * (v - half)

!     Evaluate the quadratic form
  x = u - s
  y = ABS(v) - t
  q = x**2 + y*(a*y - b*x)

!     Accept P if inside inner ellipse
  IF (q < r1) EXIT
!     Reject P if outside outer ellipse
  IF (q > r2) CYCLE
!     Reject P if outside acceptance region
  IF (v**2 < -4.0*LOG(u)*u**2) EXIT
END DO

!     Return ratio of P's coordinates as the normal deviate
fn_val = v/u
RETURN

END FUNCTION random_normal
!###############################################################
FUNCTION random_chisq(ndf) RESULT(fn_val)
IMPLICIT NONE

!     Generates a random variate from the chi-squared distribution with
!     ndf degrees of freedom


REAL, INTENT(IN) :: ndf
REAL             :: fn_val, random_gamma, half = 0.5, two = 2.0

fn_val = two * random_gamma(half*ndf)
RETURN

END FUNCTION random_chisq
!###############################################################
FUNCTION random_gamma(s) RESULT(fn_val)
IMPLICIT NONE

! Adapted from Fortran 77 code from the book:
!     Dagpunar, J. 'Principles of random variate generation'
!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

!     FUNCTION GENERATES A RANDOM GAMMA VARIATE.
!     CALLS EITHER random_gamma1 (S > 1.0)
!     OR random_exponential (S = 1.0)
!     OR random_gamma2 (S < 1.0).

!     S = SHAPE PARAMETER OF DISTRIBUTION (0 < REAL).

REAL, INTENT(IN)    :: s
REAL                :: fn_val, one = 1.0, zero = 0.0, &
                       random_gamma1, random_gamma2, random_exponential

IF (s <= zero) THEN
  WRITE(*, *) 'SHAPE PARAMETER VALUE MUST BE POSITIVE'
  STOP
END IF

IF (s > one) THEN
  fn_val = random_gamma1(s)
ELSE IF (s < one) THEN
  fn_val = random_gamma2(s)
ELSE
  fn_val = random_exponential()
END IF

RETURN
END FUNCTION random_gamma
!###############################################################
FUNCTION random_gamma1(s) RESULT(fn_val)
IMPLICIT NONE

! Uses the algorithm in
! Marsaglia, G. and Tsang, W.W. (2000) `A simple method for generating
! gamma variables', Trans. om Math. Software (TOMS), vol.26.

! Generates a random gamma deviate for shape parameter s >= 1.

REAL, INTENT(IN)    :: s
REAL                :: fn_val,random_normal
REAL                :: zero = 0.0, half = 0.5, one = 1.0

! Local variables
REAL, SAVE  :: c, d
REAL        :: u, v, x

d = s - one/3.
c = one/SQRT(9.0*d)

! Start of main loop
DO

! Generate v = (1+cx)^3 where x is random normal; repeat if v <= 0.

  DO
    x = random_normal()
    v = (one + c*x)**3
    IF (v > zero) EXIT
  END DO

! Generate uniform variable U

  CALL RANDOM_NUMBER(u)
  IF (u < one - 0.0331*x**4) THEN
    fn_val = d*v
    EXIT
  ELSE IF (LOG(u) < half*x**2 + d*(one - v + LOG(v))) THEN
    fn_val = d*v
    EXIT
  END IF
END DO

RETURN
END FUNCTION random_gamma1
!###############################################################
FUNCTION random_gamma2(s) RESULT(fn_val)
IMPLICIT NONE

! Adapted from Fortran 77 code from the book:
!     Dagpunar, J. 'Principles of random variate generation'
!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
! A GAMMA DISTRIBUTION WITH DENSITY PROPORTIONAL TO
! GAMMA2**(S-1) * EXP(-GAMMA2),
! USING A SWITCHING METHOD.

!    S = SHAPE PARAMETER OF DISTRIBUTION
!          (REAL < 1.0)

REAL, INTENT(IN)    :: s
REAL                :: fn_val

REAL :: zero = 0.0, one = 1.0, vsmall = TINY(1.0)

!     Local variables
REAL       :: r, x, w
REAL, SAVE :: a, p, c, uf, vr, d

IF (s <= zero .OR. s >= one) THEN
  WRITE(*, *) 'SHAPE PARAMETER VALUE OUTSIDE PERMITTED RANGE'
  STOP
END IF

a = one - s
p = a/(a + s*EXP(-a))
IF (s < vsmall) THEN
    WRITE(*, *) 'SHAPE PARAMETER VALUE TOO SMALL'
    STOP
END IF

c = one/s
uf = p*(vsmall/a)**s
vr = one - vsmall
d = a*LOG(a)

DO
  CALL RANDOM_NUMBER(r)
  IF (r >= vr) THEN
    CYCLE
  ELSE IF (r > p) THEN
    x = a - LOG((one - r)/(one - p))
    w = a*LOG(x)-d
  ELSE IF (r > uf) THEN
    x = a*(r/p)**c
    w = x
  ELSE
    fn_val = zero
    RETURN
  END IF

  CALL RANDOM_NUMBER(r)
  IF (one-r <= w .AND. r > zero) THEN
    IF (r*(w + one) >= one) CYCLE
    IF (-LOG(r) <= w) CYCLE
  END IF
  EXIT
END DO

fn_val = x
RETURN

END FUNCTION random_gamma2
!###############################################################
FUNCTION random_exponential() RESULT(fn_val)
IMPLICIT NONE

! Adapted from Fortran 77 code from the book:
!     Dagpunar, J. 'Principles of random variate generation'
!     Clarendon Press, Oxford, 1988.   ISBN 0-19-852202-9

! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM
! A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL
! TO EXP(-random_exponential), USING INVERSION.

REAL  :: fn_val
REAL  :: zero = 0.0

!     Local variable
REAL  :: r

DO
  CALL RANDOM_NUMBER(r)
  IF (r > zero) EXIT
END DO

fn_val = -LOG(r)
RETURN

END FUNCTION random_exponential
!***********************************************************
RECURSIVE SUBROUTINE sort_rec(array,first,last)

IMPLICIT NONE

TYPE rec_element
     INTEGER :: location
     REAL :: value
END TYPE rec_element

INTEGER, INTENT(IN) :: first,last
TYPE(rec_element), INTENT(INOUT) :: array(last)

INTEGER :: i,j, indxloc
REAL :: median, indxval

i = first
j = last

median = array(INT((first+last)/2))%value

do
   if( i > j ) exit
   do
      if( median <=  array(i)%value ) exit
      i = i+1
   enddo
   do
      if( array(j)%value <= median ) exit
      j = j-1
   enddo
   if ( i < j ) then

      indxloc = array(i)%location
      indxval = array(i)%value

      array(i)%location = array(j)%location
      array(i)%value = array(j)%value

      array(j)%location = indxloc
      array(j)%value = indxval

      i = i+1
      j = j-1
   elseif ( i == j ) then
      i=i+1
   endif
enddo

if ( first < j ) CALL sort_rec(array,first,j)
if ( i < last ) CALL sort_rec(array,i,last)

END SUBROUTINE sort_rec

!***************************************************************
!*********************************************************************