Program c09ezfe

!     C09EZF Example Program Text
!     Mark 26.1 Release. NAG Copyright 2016.

!     .. Use Statements ..
      Use nag_library, Only: c09abf, c09ecf, c09edf, c09eyf, c09ezf, dnrm2,    &
                             nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: mse, thresh
      Integer                          :: cindex, denoised, i, ifail, ilev, j, &
                                          lda, ldb, ldd, lenc, m, n, nf, nwcn, &
                                          nwct, nwl
      Character (10)                   :: mode, wavnam, wtrans
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), an(:,:), b(:,:), c(:),       &
                                          d(:,:), e(:,:)
      Integer, Allocatable             :: dwtlvm(:), dwtlvn(:)
      Integer                          :: icomm(180)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, log, real, sqrt
!     .. Executable Statements ..
      Write (nout,*) 'C09EZF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
!     Read problem parameters
      Read (nin,*) m, n
      Read (nin,*) wavnam, mode
      Write (nout,99999) wavnam, mode, m, n

!     Allocate arrays to hold the original data, A, original data plus noise,
!     AN, reconstruction using denoised coefficients, B, and randomly
!     generated noise, X.
      lda = m
      ldb = m
      Allocate (a(lda,n),an(lda,n),b(ldb,n),e(m,n))

!     Read in the original data
      Do i = 1, m
        Read (nin,*) a(i,1:n)
      End Do

!     Output the original data
      Write (nout,99997)
      Do i = 1, m
        Write (nout,99998) a(i,1:n)
      End Do

!     Fill the array AN with the original data in A plus some noise
!     and return a VisuShrink denoising threshold, thresh.
      Call create_noise(a,an,lda,m,n,thresh)

!     Output the noisy data
      Write (nout,99996)
      Do i = 1, m
        Write (nout,99998) an(i,1:n)
      End Do

!     Query wavelet filter dimensions
!     For Multi-Resolution Analysis, decomposition, wtrans = 'M'
      wtrans = 'Multilevel'
      ifail = 0
      Call c09abf(wavnam,wtrans,mode,m,n,nwl,nf,nwct,nwcn,icomm,ifail)

!     Allocate arrays to hold the coefficients, C, and the dimensions
!     of the coefficients at each level, DWTLVM, DWTLVN
      lenc = nwct
      Allocate (c(lenc),dwtlvm(nwl),dwtlvn(nwl))

!     Perform a forwards multi-level transform on the noisy data
      ifail = 0
      Call c09ecf(m,n,an,lda,lenc,c,nwl,dwtlvm,dwtlvn,icomm,ifail)

!     Reconstruct without thresholding of detail coefficients
      ifail = 0
      Call c09edf(nwl,lenc,c,m,n,b,ldb,icomm,ifail)

!     Calculate the Mean Square Error of the noisy reconstruction
      e(:,:) = a(:,:) - b(:,:)
      mse = dnrm2(m*n,e,1)
      mse = mse**2
      mse = mse/real(m*n,kind=nag_wp)
      Write (nout,99995) mse

!     Now perform the denoising by extracting each of the detail
!     coefficients at each level and applying hard thresholding

!     Allocate a 2D array to hold the detail coefficients
      ldd = dwtlvm(nwl)
      Allocate (d(ldd,dwtlvn(nwl)))

      denoised = 0
!     For each level
      Do ilev = nwl, 1, -1

!       Select detail coefficients
        Do cindex = 1, 3

!         Extract coefficients into the 2D array D
          ifail = 0
          Call c09eyf(ilev,cindex,lenc,c,d,ldd,icomm,ifail)

!         Perform the hard thresholding operation
          Do j = 1, dwtlvn(nwl-ilev+1)
            Do i = 1, dwtlvm(nwl-ilev+1)
              If (abs(d(i,j))<thresh) Then
                d(i,j) = 0.0_nag_wp
                denoised = denoised + 1
              End If
            End Do
          End Do

!         Insert the denoised coefficients back into C
          ifail = 0
          Call c09ezf(ilev,cindex,lenc,c,d,ldd,icomm,ifail)

        End Do

      End Do

!     Output the number of coefficients that were set to zero
      Write (nout,99994) denoised, nwct - dwtlvm(1)*dwtlvn(1)

!     Reconstruct original data following thresholding of detail coefficients
      ifail = 0
      Call c09edf(nwl,lenc,c,m,n,b,ldb,icomm,ifail)

!     Calculate the Mean Square Error of the denoised reconstruction
      e(:,:) = a(:,:) - b(:,:)
      mse = dnrm2(m*n,e,1)
      mse = mse**2
      mse = mse/real(m*n,kind=nag_wp)
      Write (nout,99993) mse

!     Output the denoised reconstruction
      Write (nout,99992)
      Do i = 1, m
        Write (nout,99998) b(i,1:n)
      End Do

99999 Format (1X,' MLDWT :: Wavelet  : ',A,/,1X,'          End mode : ',A,/,   &
        1X,'          M        : ',I4,/,1X,'          N        : ',I4)
99998 Format (8(F8.4,1X),:)
99997 Format (/,1X,' Original data            A  : ')
99996 Format (/,1X,' Original data plus noise AN : ')
99995 Format (/,1X,' Without denoising Mean Square Error is ',F9.6)
99994 Format (/,1X,' Number of coefficients denoised is ',I3,' out of ',I3)
99993 Format (/,1X,' With denoising Mean Square Error is ',F9.6)
99992 Format (/,1X,' Reconstruction of denoised input D : ')

    Contains

!     Subroutine fills the output array AN with the data in A
!     plus some noise taken from a normal distribution, and
!     returns the VisuShrink denoising threshold, thresh.

      Subroutine create_noise(a,an,lda,m,n,thresh)

!       .. Use Statements ..
        Use nag_library, Only: g05kff, g05skf
!       .. Parameters ..
        Integer, Parameter             :: lseed = 1
!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out) :: thresh
        Integer, Intent (In)           :: lda, m, n
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In) :: a(lda,n)
        Real (Kind=nag_wp), Intent (Out) :: an(lda,n)
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: var, xmu
        Integer                        :: genid, i, ifail, lstate, subid
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: x(:,:)
        Integer                        :: seed(lseed)
        Integer, Allocatable           :: state(:)
!       .. Executable Statements ..

!       Set up call to g05skf in order to create some random noise from
!       a normal distribution to add to the original data.
!       Initial call to RNG initializer to get size of STATE array
        seed(1) = 642521
        genid = 3
        subid = 0
        lstate = 0
        Allocate (state(lstate))
        ifail = 0
        Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!       Reallocate STATE
        Deallocate (state)
        Allocate (state(lstate))

!       Initialize the generator to a repeatable sequence
        ifail = 0
        Call g05kff(genid,subid,seed,lseed,state,lstate,ifail)

!       Set the distribution parameters for the random noise.
        xmu = 0.0_nag_wp
        var = 0.1E-3_nag_wp

        Allocate (x(m,n))

!       Generate the noise variates
        ifail = 0
        Do i = 1, n
          Call g05skf(m,xmu,var,state,x(1,i),ifail)
        End Do

!       Add the noise to the original input and save in AN
        an(:,:) = a(:,:) + x(:,:)

!       Calculate the threshold based on VisuShrink denoising
        thresh = sqrt(var)*sqrt(2._nag_wp*log(real(m*n,kind=nag_wp)))

      End Subroutine create_noise

    End Program c09ezfe