C     PROGRAM FOR COMPUTING THE BAYESIAN FILTERING ESTIMATES 
C     PX0(N): VECTOR OF INITIAL DIST OF X_0
C     PM(NM): PRIOR OF MU,THEN POSTERIOR
C     PS(NS): PRIOR OF SIGMA,THEN POSTERIOR
C     
      PROGRAM DONDDIST
      include 'par12t2.f'
      INTEGER K,L,M,W
      REAL*8 PX0(N),PM(NM),PS(NS),PR(NR),P(N,NM,NS,NR),MU(NM),SG(NS),
     &RO(NR),X(N),PX(N)
c     ,PMSR(NM,NS,NR)
c     ASSIGN INITIAL DISTRIBUTION OF X_0 ---------------
      DO 10 K=1,N
         PX0(K)=0.0d0
 10   CONTINUE
c      WRITE(*,*)"END LOOP 10"
      PX0(I0)=1.0d0
c      WRITE(*,*) "PX0",PX0
c      pause "px0"
C     --------------------------------------------------
C     ASSIGN PRIOR FOR MU, PM(NM) ----------------------
      DO 20 L=1,NM
         PM(L)=1.0d0/dble(NM)
 20   CONTINUE
c      WRITE(*,*) "END LOOP 20"
c      write(*,*) "PM",PM
c      pause "pm"
C     --------------------------------------------------
C     ASSIGN PRIOR FOR SIGMA, PS(NS) -------------------
      DO 30 M=1,NS
         PS(M)=1.0d0/dble(NS)
 30   CONTINUE
c      WRITE(*,*) "END LOOP 30"
c      write(*,*) "ps", ps
c      pause "ps"
C     --------------------------------------------------
C     ASSIGN PRIOR FOR RHO, PR(NR)
      DO 40 W=1,NR
         PR(W)=1.0D0/DBLE(NR)
 40   CONTINUE
c     COMPUTE INITIAL DIST OF P(X_0,MU,SIGMA | Z_0) ----
      DO 600 K=1,N
         DO 400 L=1,NM
            DO 200 M=1,NS
            DO 190 W=1,NR
                  P(K,L,M,W)=PX0(K)*PM(L)*PS(M)*PR(W)
 190           CONTINUE
 200           CONTINUE
 400           CONTINUE
 600           CONTINUE
C     ------------------------------------------------
c      WRITE(*,*) "P INITIAL", P
c      pause "initial p"
C     ------------------------------------------------
      CALL DONDDIST3(P)
C     COMPUTE THE POSTERIOR OF MU-----------------------
      OPEN(UNIT=17,FILE="12MU.LST",STATUS="UNKNOWN")
      DO 900 L=1,NM
        PM(L)=0d0
        MU(L)=MUMIN+(dble(L)-1.000d0)*MUL
        DO 930 K=1,N
           DO 940 W=1,NR
           DO 950 M=1,NS
           PM(L)=PM(L)+P(K,L,M,W)
 950       CONTINUE
 940       CONTINUE
 930    CONTINUE
         WRITE(17,'(g16.6E4,G16.6E4)') MU(L),PM(L)
 900  CONTINUE
      CLOSE(17)
c      WRITE(*,*) "END LOOP 900"
C     -----------------------------------------------------
C     COMPUTE THE POSTERIOR OF SIGMA-----------------------
      OPEN(UNIT=27,FILE="12SG.LST",STATUS="UNKNOWN")
      DO 1000 M=1,NS
         PS(M)=0d0
         SG(M)=SGMIN+(dble(M)-1.00d0)*SGL
         DO 1030 K=1,N
            DO 1040 W=1,NR
            DO 1050 L=1,NM
               PS(M)=PS(M)+P(K,L,M,W)
 1050       CONTINUE
 1040       CONTINUE
 1030    CONTINUE
         WRITE(27,'(G16.6E4,G16.6E4)') SG(M),PS(M)
 1000 CONTINUE
      CLOSE(27)
C     COMPUTE THE POSTERIOR OF RHO-----------------------
      OPEN(UNIT=37,FILE="12RO.LST",STATUS="UNKNOWN")
      DO 1100 W=1,NR
         PR(W)=0d0
         RO(W)=ROMIN+(dble(W)-1.00d0)*ROL
         DO 1130 K=1,N
            DO 1140 M=1,NS
            DO 1150 L=1,NM
               PR(W)=PR(W)+P(K,L,M,W)
 1150       CONTINUE
 1140       CONTINUE
 1130    CONTINUE
         WRITE(37,'(G16.6E4,G16.6E4)') RO(W),PR(W)
 1100 CONTINUE
      CLOSE(37)
C     COMPUTE THE POSTERIOR OF X-----------------------
      OPEN(UNIT=47,FILE="12X.LST",STATUS="UNKNOWN")
      DO 1200 K=1,N
         PX(K)=0d0
         X(K)=XMIN-LL+(dble(K)-1.00d0)*0.031250d0
         DO 1230 W=1,NR
            DO 1240 M=1,NS
            DO 1250 L=1,NM
               PX(K)=PX(K)+P(K,L,M,W)
 1250       CONTINUE
 1240       CONTINUE
 1230    CONTINUE
         WRITE(47,'(G16.6E4,G16.6E4)') X(K), PX(K)
 1200 CONTINUE
      CLOSE(47)
C     -----------------------------------------------------
      END

      SUBROUTINE DONDDIST3(P)
      EXTERNAL LD1,LD2,LD3
      INCLUDE 'par12t2.f'
      INTEGER I,K,L,M,NTI,J,W
      REAL*8 P(N,NM,NS,NR),A(N),B(N),X2,DT
      REAL*8 MU,SG,XB,XC,LD1,LD2,LD3
      REAL*8 PI1(N)
c     store data in X(I) and DT
c      OPEN(UNIT=32,FILE='midout12.flt',STATUS="UNKNOWN")
      OPEN(UNIT=10, FILE='MUp.dat', STATUS="UNKNOWN")
      OPEN(UNIT=20, FILE='SGp.dat', STATUS="UNKNOWN")
      OPEN(UNIT=30, FILE='ROp.dat', STATUS="UNKNOWN")
      OPEN(UNIT=45, FILE='Xb.dat', STATUS="UNKNOWN")
      OPEN(UNIT=15, FILE='MUb.dat', STATUS="UNKNOWN")
      OPEN(UNIT=25, FILE='SGb.dat', STATUS="UNKNOWN")
      OPEN(UNIT=35, FILE='ROb.dat', STATUS="UNKNOWN")
      OPEN(UNIT=1,FILE='t2df.dat',      
     &STATUS='OLD')
      DO 6 I=1,NDATA
         READ(1,*) X2,DT
C         write(*,*) X2,dt
c     loop 6: updating the cond prob from time t_1 to t_ndata
c     loop 100: approXimate p(K,L,M,W,t_i-)
         DO 100 M=1,NS
            SG=SGMIN+(dble(M)-1.00d0)*SGL
            DO 110 L=1,NM
               MU=MUMIN+(dble(L)-1.00000d0)*MUL
c               write(*,*) mu,SG
C     ASSIGN A(K), B(K) ACCORDINGLY------------------------------------
            IF (MU .GE. 0D0) THEN
c               pause "mu ge 0"
c     "*8" should be changed to "*16" when 1/16
               DO 120 K=1,N
         A(K)=0.50D0*(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
     & *(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
     & + MU*(DBLE(K)-1.0D0+(XMIN-LL)*32.0D0)
         B(K)=0.50D0*(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
     & *(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
 120           CONTINUE
            ELSE 
               DO 125 K=1,N
         A(K)=0.50D0*(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
     & *(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
         B(K)=0.50D0*(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
     & *(SG*(DBLE(K)-1.0D0+(XMIN-LL)*32.0d0))
     & - MU*(DBLE(K)-1.0D0+(XMIN-LL)*32.0D0)
 125           CONTINUE
            END IF
c            write(*,*) "end if of A B",A,B
c     -----------------------------------------------------
c     refine the partition based on the interarrival time 
c     -----------------------------------------------------
      DO 133 W=1,NR
C         RO=ROMIN+(DBLE(W)-1.00D0)*ROL
         NTI=IDNINT(DT/LT)+1
c         print*,lt,nti
         DO 55 J=1,NTI
         DO 50 K=1,N
            PI1(K)=P(K,L,M,W)
 50      CONTINUE
                DO 130 K=1,N
            IF (K .EQ. 1) THEN
c               pause "k=1"
               P(K,L,M,W)=PI1(K)+(-(A(K)+B(K))*PI1(K)
     &+B(K+1)*PI1(K+1))*DT/DBLE(NTI)
            ELSE IF (K .EQ. N) THEN
c               pause "k=n"
               P(K,L,M,W)=PI1(K)+(A(K-1)*PI1(K-1)-(A(K)+B(K))
     &*PI1(K))*DT/DBLE(NTI)
            ELSE
               P(K,L,M,W)=PI1(K)+(A(K-1)*PI1(K-1)-(A(K)+B(K))
     &*PI1(K)+B(K+1)*PI1(K+1))*DT/DBLE(NTI)
            END IF
c            if (pi2(k) .gt. 0) then
c               print*, "i,k,l,m,w",i,k,l,m,w,pi1(k),pi1(k-1)
c               print*,pi1(k+1),a(k-1),a(k),b(k),b(k+1),pi2(k)
c               pause
c               end if
            if (p(k,l,m,w) .lt. 0d0) then
          write(*,*) "sg",sg,"mu",mu
               write(*,*) "i,k,l,m,W",i,k,l,m,W,"DT",DT
               write(*,*) "a_k,b_k,a_k-1,b_k+1",a(k),b(k),a(k-1),b(k+1)
               write(*,*) "PI1(K-1),PI1(K),PI1(K+1)",
     &PI1(K-1),PI1(K),PI1(K+1)
               pause "negative PROBABILITY"
            end if
 130            CONTINUE
 55         CONTINUE
 133        CONTINUE
c          write (*,*) 'end loop 133'
 110     CONTINUE
 100     CONTINUE
c         if (i .eq. 6) then
c         pause "i=6,end loop 100"
c         end if
         IF (I .EQ. 1) THEN
            DO 450 L=1,NM
               DO 452 M=1,NS
                  DO 454 W=1,NR
                  DO 456 K=1,N
c                    IF (P(K,L,M,W) .GT. 0D0) THEN
c                       WRITE(32,*) K,L,M,W,P(K,L,M,W)
c                     END IF
C                        XA=(DBLE(K)+(XMIN-LL)*8D0-1D0)/8D0 
C                 WRITE(*,*) "XA",XA
C                 PAUSE
C                 IF (X2 .EQ. XA) THEN
C                    WRITE(*,*) "XA",XA,X2,p(k,l,m,W),K,L,M,I
C                    PAUSE
C                 END IF
 456              CONTINUE
 454              CONTINUE
 452           CONTINUE
 450        CONTINUE
         END IF
c         pause
c     ----------------------------------------------------------------
C     "IF STRUCTURE": ASSIGN FOR DIFFERENT INTENSITY FUNCTIONS TO OBTAIN
C     P(T_I)
       XB=X2*2D0
       XC=X2*4D0
c       write (*,*) 'Xb',xb,XC
C     CASE 1: X2 IS AN INTEGER OR HALF
         IF (DMOD(XB,1d0) .EQ. 0) THEN
C       write (*,*) "inside if"
            CALL DEND(LD1,P,X2)
C     CASE 2: X2 IS ODD QUARTERS
         ELSE IF (DMOD(XC,1d0) .EQ. 0) THEN
C            write(*,*) "INSIDE else if"
            CALL DEND(LD2,P,X2)
C     CASE 3: X2 IS ODD EIGHTH
         ELSE
C            write (*,*) "INSIDE else"
            CALL DEND(LD3,P,X2)
         END IF
            CALL WOUTPUT(P,I,X2,DT)
 6    continue
      close(1)
      CLOSE(10)
      CLOSE(15)
      CLOSE(20)
      CLOSE(25) 
      CLOSE(30)
      CLOSE(35)
      CLOSE(45)      
         RETURN
         END

      SUBROUTINE DEND(FUNC,P,X)
      include 'par12t2.f'
      INTEGER K,L,M,W
      REAL*8 FUNC,P(N,NM,NS,NR),X,XA,TL
C     SUM UP THE TOTEL TL
c         write(*,*) "func", x, func(x, 78d0)
         TL=0D0
         DO 180 L=1,NM
            DO 182 M=1,NS
               DO 184 W=1,NR
               DO 186 K=1,N
C                 XA=(DBLE(K)-1D0+(XMIN-LL)*8D0)/8D0  (FOR 1/8)
C                 WRITE(*,*) P1(K,L,M),FUNC(X,XA)
C      FOR 1/32
                  XA=(DBLE(K)-1D0+(XMIN-LL)*32D0)/32D0 
                 TL=TL+P(K,L,M,W)*FUNC(X,XA,W)
c                     IF (P1(K,L,M,W) .GT. 0D0) THEN
c                        WRITE(32,*) K,L,x,xa,P1(K,L,M,W)
c                      write(32,*) func(x,xa,w),func(x,xa,w)*p1(k,l,m,w)
c                      END IF
 186            CONTINUE
 184            CONTINUE
 182         CONTINUE
 180      CONTINUE
          IF (TL .EQ. 0D0) THEN 
             PAUSE "TL=0"
          END IF
c         write(32,*) "END LOOP 180", "tl",tl
         DO 280 L=1,NM
            DO 282 M=1,NS
               DO 284 W=1,NR
               DO 286 K=1,N
c                 XA=(DBLE(K)-1D0+(XMIN-LL)*8D0)/8D0  (for 1/8)
                  XA=(DBLE(K)-1D0+(XMIN-LL)*32D0)/32D0 
                 P(K,L,M,W)=P(K,L,M,W)*FUNC(X,XA,W)/TL
c                     IF (P(K,L,M,W) .GT. 0D0) THEN
c                        WRITE(32,*) K,L,xa,P(K,L,M,W)
c                     END IF
 286            CONTINUE
 284            CONTINUE
 282         CONTINUE
 280      CONTINUE
c         write(*,*) "END LOOP 280"
      RETURN
      END

      
      SUBROUTINE WOUTPUT(P,I,X2,DT)
      INCLUDE 'par12t2.f'
      INTEGER I,K,L,M,W
      REAL*8 PM(NM),PS(NS),PR(NR),P(N,NM,NS,NR),MU(NM),SG(NS),
     &RO(NR),X2,DT,EMU,EMU2,SMU,ESG,ESG2,SSG,ERO,ERO2,SRO,EX,EX2,
     &SX,PX(N),X(N)
      EMU=0D0
      EMU2=0D0
      SMU=0D0
      ESG=0D0
      ESG2=0D0
      SSG=0D0
      ERO=0D0
      ERO2=0D0
      SRO=0D0
      EX=0D0
      EX2=0D0
      SX=0D0
      DO 901 L=1,NM
        PM(L)=0d0
        MU(L)=MUMIN+(dble(L)-1.000d0)*MUL
        DO 931 K=1,N
           DO 941 W=1,NR
           DO 951 M=1,NS
           PM(L)=PM(L)+P(K,L,M,W)
 951    CONTINUE
 941  CONTINUE
 931  CONTINUE
         if ((i .eq. 250) .or. (i .eq. 500) .or. (i .eq. 750) .or. 
     & (i .eq. 1000) .or. (i .eq. 1500) .or. 
     & (i .eq. 3000) .or. (i .eq. 1) .or.
     & (i .eq. 5000) .or. (i .eq. 10000) .or. (i .eq.15000)
     & .or. (i .eq. 20000) .or. (i .eq. 30000)
     & )  then 
         WRITE(10,'(G16.6E4)') PM(L)
         end if
         EMU=EMU+MU(L)*PM(L)
         EMU2=EMU2+MU(L)*MU(L)*PM(L)
 901  CONTINUE
        SMU=DSQRT(EMU2-EMU*EMU)
      WRITE(15,'(G16.6E4,g16.6E4)') EMU,SMU
      DO 1001 M=1,NS
         PS(M)=0d0
         SG(M)=SGMIN+(dble(M)-1.00d0)*SGL
         DO 1031 K=1,N
            DO 1041 W=1,NR
            DO 1051 L=1,NM
               PS(M)=PS(M)+P(K,L,M,W)
 1051       CONTINUE
 1041    CONTINUE
 1031 CONTINUE
         if ((i .eq. 250) .or. (i .eq. 500) .or. (i .eq. 750) .or. 
     & (i .eq. 1000) .or. (i .eq. 1500) .or. 
     & (i .eq. 3000) .or. (i .eq. 1) .or.
     & (i .eq. 5000) .or. (i .eq. 10000) .or. (i .eq.15000)
     & .or. (i .eq. 20000) .or. (i .eq. 30000)
     & )  then 
         WRITE(20,'(G16.6E4)') PS(M)
         end if
        ESG=ESG+SG(M)*PS(M)
         ESG2=ESG2+SG(M)*SG(M)*PS(M)
 1001 CONTINUE
        SSG=DSQRT(ESG2-ESG*ESG)
         WRITE(25,'(G16.6E4,g16.6E4)') ESG,SSG
      DO 1101 W=1,NR
         PR(W)=0d0
         RO(W)=ROMIN+(dble(W)-1.00d0)*ROL
         DO 1131 K=1,N
            DO 1141 M=1,NS
            DO 1151 L=1,NM
               PR(W)=PR(W)+P(K,L,M,W)
 1151       CONTINUE
 1141    CONTINUE
 1131 CONTINUE
         if ((i .eq. 250) .or. (i .eq. 500) .or. (i .eq. 750) .or. 
     & (i .eq. 1000) .or. (i .eq. 1500) .or. 
     & (i .eq. 3000) .or. (i .eq. 1) .or.
     & (i .eq. 5000) .or. (i .eq. 10000) .or. (i .eq.15000)
     & .or. (i .eq. 20000) .or. (i .eq. 30000)
     & )  then 
         WRITE(30,'(G16.6E4)') PR(W)
         end if
        ERO=ERO+RO(W)*PR(W)
         ERO2=ERO2+RO(W)*RO(W)*PR(W)
 1101 CONTINUE
        SRO=DSQRT(ERO2-ERO*ERO)
      WRITE(35,'(G16.6E4,g16.6E4)') ERO,SRO
        EX=0D0
         EX2=0D0
      DO 1201 K=1,N
         PX(K)=0d0
         X(K)=XMIN-LL+(dble(K)-1.00d0)*0.031250d0
         DO 1231 W=1,NR
            DO 1241 M=1,NS
            DO 1251 L=1,NM
               PX(K)=PX(K)+P(K,L,M,W)
 1251       CONTINUE
 1241    CONTINUE
 1231 CONTINUE
        EX=EX+X(K)*PX(K)
         EX2=EX2+X(K)*X(K)*PX(K)
 1201 CONTINUE
        SX=DSQRT(EX2-EX*EX)
      WRITE(45,'(G16.6E4,g16.6E4)') EX,SX
      RETURN
      END

      FUNCTION LD1(X,Y,W)
      INCLUDE 'par12t2.f'
      REAL*8 X,LD1,Y,RHO
      INTEGER XA,W
c     add 0.000001 for protection of roundoff error
      XA=IDINT(8.0D0*X)-IDNINT(8.0D0*Y+0.000000001)
      RHO=ROMIN+(DBLE(W)-1.0D0)*ROL
C      WRITE(*,*) "INSIDE LD1"
      IF (XA .EQ. 0) THEN
         LD1=(1.0D0-RHO)*(1.0D0+BETA*RHO)
      ELSE IF (IABS(XA) .eq. 1) Then
         LD1=0.50d0*(1.0d0-rho)*rho*(1.0d0+rho*beta)
     &+(1.0d0-rho)*beta
      ELSE
         LD1=0.50d0*(1.0D0-RHO)*(RHO**(IABS(XA)))*
     & (1.0D0+BETA*(RHO+1.0d0/RHO))
      END IF
C      write(*,*) X,Y,W,XA,RHO, ROMIN,ROL,BETA,LD1
C      PAUSE
      RETURN
      END

      FUNCTION LD2(X,Y,W)
      INCLUDE 'par12t2.f'
      REAL*8 X,LD2,Y,RHO
      INTEGER XA,W
      XA=IDINT(8.0D0*X)-IDNINT(8.0D0*Y+0.000000001)
      RHO=ROMIN+(DBLE(W)-1.0D0)*ROL
      IF (XA .EQ. 0) THEN
         LD2=(1.0D0-RHO)*(1.0D0+ALPHA*RHO)
      ELSE IF (IABS(XA) .eq. 1) Then
         LD2=0.50d0*(1.0d0-rho)*rho*(1.0d0+rho*alpha)
     &+(1.0d0-rho)*alpha
      ELSE
         LD2=0.50D0*(1.0D0-RHO)*(RHO**(IABS(XA)))*
     & (1.0D0+ALPHA*(RHO+1/RHO))
      END IF
      RETURN
      END
      
      FUNCTION LD3(X,Y,W)
      INCLUDE 'par12t2.f'
      REAL*8 X,Y,LD3,RHO
      INTEGER XA,W
      XA=IDINT(8.0D0*X)-IDNINT(8.0D0*Y+0.000000001)
      RHO=ROMIN+(DBLE(W)-1.0D0)*ROL
      IF (XA .EQ. 0) THEN
        LD3=(1.0D0-RHO)*(1.0D0-ALPHA-BETA)
      ELSE
        LD3=0.50D0*(1.0D0-ALPHA-BETA)*(1.0D0-RHO)*(RHO**(IABS(XA))) 
      END IF 
      RETURN
      END