    MODULE TOOL


    CONTAINS

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

    SUBROUTINE networksize(GG, MAX, path1, path2, N)
    IMPLICIT NONE
    CHARACTER(100), INTENT(IN) :: path1, path2
    CHARACTER(40) :: numchr
    INTEGER :: g, i, io
    INTEGER, INTENT(IN) :: GG, MAX
    REAL, DIMENSION(MAX,MAX) :: temp
    INTEGER, DIMENSION(GG), INTENT(OUT) :: N

    DO g = 1,GG
        WRITE(numchr,*) g
        !  Read rows from file "**.dat"
        OPEN (g, file = trim(adjustl(path1)) // trim(adjustl(path2)) // '/age/age' // trim(adjustl(numchr)) // '.dat', status='old')         
        N(g)=0
        DO i=1, max
            READ(g,*,IOSTAT=io)  temp(i,1)   
            IF (io < 0) THEN
                EXIT
            END IF
            N(g) = N(g) + 1
        END DO
        CLOSE(g)    
    END DO
    END SUBROUTINE networksize


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

    SUBROUTINE define_variable(GG, age, sex, race, sport, network, group, N, path1, path2, MAX)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: GG, MAX
    INTEGER, DIMENSION(GG), INTENT(IN) :: N
    INTEGER :: i, j, g
    REAL, DIMENSION(MAX,MAX,GG), INTENT(INOUT) :: age, sex, race, sport, network
    REAL, DIMENSION(MAX,19,GG), INTENT(INOUT) :: group
    CHARACTER(100), INTENT(IN) :: path1, path2
    CHARACTER(40) :: numchr


    DO g=1, GG

        WRITE(numchr,*) g
        OPEN (100, file = trim(adjustl(path1)) // trim(adjustl(path2)) // '/age/age' // trim(adjustl(numchr)) // '.dat', status='old') 
        OPEN (200, file = trim(adjustl(path1)) // trim(adjustl(path2)) // '/sex/sex' // trim(adjustl(numchr)) // '.dat', status='old') 
        OPEN (300, file = trim(adjustl(path1)) // trim(adjustl(path2)) // '/race/race' // trim(adjustl(numchr)) // '.dat', status='old') 
        OPEN (400, file = trim(adjustl(path1)) // trim(adjustl(path2)) // '/sport/sport' // trim(adjustl(numchr)) // '.dat', status='old') 
        OPEN (500, file = trim(adjustl(path1)) // trim(adjustl(path2)) // '/network/network' // trim(adjustl(numchr)) // '.dat', status='old') 
        OPEN (600, file = trim(adjustl(path1)) // trim(adjustl(path2)) // '/group/group' // trim(adjustl(numchr)) // '.dat', status='old') 

        READ(100,*) ((age(i,j,g),j=1,N(g)),i=1,N(g))
        READ(200,*) ((sex(i,j,g),j=1,N(g)),i=1,N(g))
        READ(300,*) ((race(i,j,g),j=1,N(g)),i=1,N(g))
        READ(400,*) ((sport(i,j,g),j=1,N(g)),i=1,N(g))
        READ(500,*) ((network(i,j,g),j=1,N(g)),i=1,N(g))
        READ(600,*) ((group(i,j,g),j=1,19),i=1,N(g))

        CLOSE(100)
        CLOSE(200)
        CLOSE(300)
        CLOSE(400)
        CLOSE(500)
        CLOSE(600)
    END DO

    END SUBROUTINE define_variable 


    !************************************************!
    FUNCTION  zeros(d1,d2) RESULT(val_fun)
    IMPLICIT  NONE
    INTEGER, INTENT(IN) :: d1, d2
    INTEGER :: i, j
    REAL, DIMENSION(d1,d2) :: val_fun  
    val_fun=0.0
    RETURN
    END FUNCTION  zeros


    !************************************************!
    FUNCTION  ones(d1,d2) RESULT(val_fun)
    IMPLICIT  NONE
    INTEGER, INTENT(IN) :: d1, d2
    INTEGER :: i, j
    REAL, DIMENSION(d1,d2) :: val_fun  
    val_fun=1.0
    RETURN
    END FUNCTION  ones


    !************************************************!
    FUNCTION  EYE(d) RESULT(val_fun)
    IMPLICIT  NONE
    INTEGER, INTENT(IN) :: d
    INTEGER :: i, j
    REAL, DIMENSION(d,d) :: val_fun
    DO i=1, d
        DO j=1, d
            val_fun(i,j)=0.0
        END DO
        val_fun(i,i)=1.0
    END DO
    RETURN
    END FUNCTION  EYE


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

    FUNCTION reshape_cov(s,n) RESULT(val_fun)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: n
    INTEGER :: i, j
    REAL, DIMENSION(n,n), INTENT(IN) :: s
    REAL, DIMENSION(n*(n+1)/2) :: val_fun

    IF (n==1) THEN
        val_fun(1)=s(1,1)
    ELSE
        DO j=1, n        
            DO i=1, j            
                val_fun(j*(j-1)/2+i)=s(i,j)                    
            END DO
        END DO
    END IF
    END FUNCTION reshape_cov


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

    FUNCTION trace(matrix,m) RESULT(val_fun)
    INTEGER, INTENT(IN) :: m
    REAL, DIMENSION(m,m), INTENT(IN) :: matrix
    REAL :: val_fun
    INTEGER :: i, j

    val_fun=0.0
    DO i=1, m
        DO j=1, m
            IF(j==i) val_fun=val_fun+matrix(i,j)
        END DO 
    END DO
    RETURN
    END FUNCTION trace

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

    FUNCTION mvnpdf(x,v,n) RESULT(val_fun)
    REAL, PARAMETER          :: pi=3.1415927
    INTEGER                  :: INDX(n), INFO,IPIV(n) , J  
    INTEGER, INTENT(IN)      :: n
    REAL, INTENT(IN)         :: x(n), v(n,n) 
    REAL                     :: x_t(n), inv_v(n,n)
    REAL(kind=8)             :: val_fun,det


    det=FindDet(v, n)
    x_t=x
    CALL SGESV(n,1,v,n,IPIV,x_t,n,INFO)
    val_fun= (2*pi)**(-(n)/2.0)*(det**(-0.5))*exp(-0.5*dot_product(x,x_t)) 
    RETURN
    END FUNCTION mvnpdf


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

    FUNCTION logmvnpdf(x,v,n) RESULT(val_fun)
    INTEGER, INTENT(IN)      :: n
    REAL, INTENT(IN)         :: x(n), v(n,n) 
    REAL                     :: x_t(n), inv_v(n,n)
    REAL(kind=8)             :: val_fun,det
    REAL, PARAMETER :: pi=3.1415927
    INTEGER :: INDX(n),  INFO,IPIV(n) , J  

    det=FindDet(v, n)
    x_t=x
    CALL SGESV(n,1,v,n,IPIV,x_t,n,INFO)
    val_fun= (-0.5*n)*log(2*pi)-0.5*log(det)-0.5*dot_product(x,x_t) 
    RETURN
    END FUNCTION logmvnpdf



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

    FUNCTION diag(x,n) RESULT(digx)

    INTEGER, INTENT(IN)      :: n
    REAL, DIMENSION(n,n)       :: x
    REAL, DIMENSION(n)     :: digx 
    INTEGER :: i
    digx=0.0
    DO i=1,n
        digx(i)=x(i,i)
    END DO
    RETURN
    END FUNCTION diag


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

    FUNCTION diagonalize(x,n) RESULT(digx)

    INTEGER, INTENT(IN)      :: n
    REAL, DIMENSION(n)       :: x
    REAL, DIMENSION(n,n)     :: digx 
    INTEGER :: i
    digx=0.0
    DO i=1,n
        digx(i,i)=x(i)
    END DO
    RETURN
    END FUNCTION diagonalize



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


    SUBROUTINE FINDInv(matrix, inverse, n, errorflag)
    !Subroutine to find the inverse of a square matrix
    !Author : Louisda16th a.k.a Ashwith J. Rego
    !Reference : Algorithm has been well explained in:
    !http://math.uww.edu/~mcfarlat/inverse.htm           
    !http://www.tutor.ms.unimelb.edu.au/matrix/matrix_inverse.html
    IMPLICIT NONE
    !Declarations
    INTEGER, INTENT(IN) :: n
    INTEGER, INTENT(OUT) :: errorflag  !Return error status. -1 for error, 0 for normal
    REAL, INTENT(IN), DIMENSION(n,n) :: matrix  !Input matrix
    REAL, INTENT(OUT), DIMENSION(n,n) :: inverse !Inverted matrix

    LOGICAL :: FLAG = .TRUE.
    INTEGER :: i, j, k, l
    REAL :: m
    REAL, DIMENSION(n,2*n) :: augmatrix !augmented matrix

    !Augment input matrix with an identity matrix
    DO i = 1, n
        DO j = 1, 2*n
            IF (j <= n ) THEN
                augmatrix(i,j) = matrix(i,j)
            ELSE IF ((i+n) == j) THEN
                augmatrix(i,j) = 1
            Else
                augmatrix(i,j) = 0
            ENDIF
        END DO
    END DO

    !Reduce augmented matrix to upper traingular form
    DO k =1, n-1
        IF (augmatrix(k,k) == 0) THEN
            FLAG = .FALSE.
            DO i = k+1, n
                IF (augmatrix(i,k) /= 0) THEN
                    DO j = 1,2*n
                        augmatrix(k,j) = augmatrix(k,j)+augmatrix(i,j)
                    END DO
                    FLAG = .TRUE.
                    EXIT
                ENDIF
                IF (FLAG .EQV. .FALSE.) THEN
                    PRINT*, "Matrix is non - invertible"
                    inverse = 0
                    errorflag = -1
                    return
                ENDIF
            END DO
        ENDIF
        DO j = k+1, n			
            m = augmatrix(j,k)/augmatrix(k,k)
            DO i = k, 2*n
                augmatrix(j,i) = augmatrix(j,i) - m*augmatrix(k,i)
            END DO
        END DO
    END DO

    !Test for invertibility
    DO i = 1, n
        IF (augmatrix(i,i) == 0) THEN
            PRINT*, "Matrix is non - invertible"
            inverse = 0
            errorflag = -1
            return
        ENDIF
    END DO

    !Make diagonal elements as 1
    DO i = 1 , n
        m = augmatrix(i,i)
        DO j = i , (2 * n)				
            augmatrix(i,j) = (augmatrix(i,j) / m)
        END DO
    END DO

    !Reduced right side half of augmented matrix to identity matrix
    DO k = n-1, 1, -1
        DO i =1, k
            m = augmatrix(i,k+1)
            DO j = k, (2*n)
                augmatrix(i,j) = augmatrix(i,j) -augmatrix(k+1,j) * m
            END DO
        END DO
    END DO				

    !store answer
    DO i =1, n
        DO j = 1, n
            inverse(i,j) = augmatrix(i,j+n)
        END DO
    END DO
    errorflag = 0
    END SUBROUTINE FINDinv


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


    FUNCTION FindDet(matrix, n) RESULT(Det)

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

    !Function to find the determinant of a square matrix
    !Author : Louisda16th a.k.a Ashwith J. Rego
    !Description: The subroutine is based on two key points:
    !1] A determinant is unaltered when row operations are performed: Hence, using this principle,
    !row operations (column operations would work as well) are used
    !to convert the matrix into upper traingular form
    !2]The determinant of a triangular matrix is obtained by finding the product of the diagonal elements
    !
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: n
    REAL, DIMENSION(n,n) :: matrix	
    REAL(kind=8) :: m, temp, Det
    INTEGER :: i, j, k, l
    LOGICAL :: DetExists = .TRUE.
    l = 1
    !Convert to upper triangular form
    DO k = 1, n-1
        IF (matrix(k,k) == 0) THEN
            DetExists = .FALSE.
            DO i = k+1, n
                IF (matrix(i,k) /= 0) THEN
                    DO j = 1, n
                        temp = matrix(i,j)
                        matrix(i,j)= matrix(k,j)
                        matrix(k,j) = temp
                    END DO
                    DetExists = .TRUE.
                    l=-l
                    EXIT
                ENDIF
            END DO
            IF (DetExists .EQV. .FALSE.) THEN
                Det = 0
                return
            END IF
        ENDIF
        DO j = k+1, n
            m = matrix(j,k)/matrix(k,k)
            DO i = k+1, n
                matrix(j,i) = matrix(j,i) - m*matrix(k,i)
            END DO
        END DO
    END DO

    !Calculate determinant by finding product of diagonal elements
    Det = l
    DO i = 1, n
        Det = Det * matrix(i,i)
    END DO
    RETURN
    END FUNCTION FindDet

    !**********************************!
    FUNCTION cov(s,n,t) RESULT(val_fun)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: n, t
    INTEGER :: i, j
    REAL, DIMENSION(n) :: mean_i, var_i
    REAL, DIMENSION(n,t), INTENT(IN) :: s
    REAL, DIMENSION(n,n) :: val_fun

    DO i=1, n
        mean_i(i)=sum(s(i,:))/t
        var_i(i)=0.0
        DO j=1, t
            var_i(i)=var_i(i)+(s(i,j)-mean_i(i))**2
        END DO
        var_i(i)=var_i(i)/t
    END DO

    IF (n==1) THEN
        val_fun(1,1)=var_i(1)
        !RESHAPE(val_fun, (/n/)) 
    ELSE
        DO i=1, n
            val_fun(i,i)=var_i(i)
            DO j=i+1, n
                val_fun(i,j)=dot_product(s(i,:)-mean_i(i),s(j,:)-mean_i(j))/t
                val_fun(j,i)=val_fun(i,j)
            END DO
        END DO
    END IF
    END FUNCTION cov

    !**********************************!
    FUNCTION var(s,t) RESULT(val_fun)
    IMPLICIT NONE
    INTEGER, INTENT(IN) :: t
    INTEGER :: i
    REAL :: mean
    REAL, DIMENSION(t), INTENT(IN) :: s
    REAL :: val_fun

    mean=sum(s)/t
    val_fun=0.0
    DO i=1, t  
        val_fun=val_fun+(s(i)-mean)**2
    END DO
    val_fun=val_fun/t

    END FUNCTION var


   


    END MODULE TOOL