!**************************************************
program f90_KMEANPP
!**************************************************
    implicit none
    character::strcom*100         !Comment
    integer::mcol                 !Number of variables
    integer::ndata                !Number of data sets
    integer::kk                   !Number of cluster
    integer::mds                  !minimum data size of each cluster
    real(8),allocatable::xd(:,:)  !Array for input data
    integer,allocatable::cl(:)    !Index of cluster
    real(8),allocatable::dis(:,:) !Distance from the center of a group
    real(8),allocatable::w1(:)
    real(8),allocatable::wxd(:,:,:)
    real(8),allocatable::wxm(:,:)
    real(8),allocatable::wva(:,:)
    real(8),allocatable::wvo(:,:)
    real(8),allocatable::wwvo(:,:,:)
    integer,allocatable::num(:)
    real(8),allocatable::wx(:)
    real(8),allocatable::wm(:)
    integer,parameter::itmax=100
    real(8),parameter::eps=1D-10
    real(8)::dismin
    real(8),allocatable::dist(:) ! Sum of distance of each group
    real(8)::dists0,dists1       ! Total sum of distance
    integer::nzero               ! Index of no element
    integer::ranseedflag=0

    integer::iii,i,j,k,j1,j2,kw
    integer::imethod
    character::fnameR*50,fnameW*50
    character::dummy*100,linebuf*1000

    !-------------------------------------------
    ! Input of cammand line arguments
    !-------------------------------------------
    call getarg(1,dummy);read(dummy,*) imethod
    call getarg(2,dummy);read(dummy,*) kk
    call getarg(3,dummy);read(dummy,*) mds
    call getarg(4,fnameR)
    call getarg(5,fnameW)
    !-------------------------------------------
    ! Data input
    !-------------------------------------------
    open(11,file=fnameR,status='old')
    read(11,'(a)') strcom
    read(11,*) mcol,ndata
    if(mds<=mcol)mds=mcol+1
    allocate(xd(1:ndata,1:mcol))
    allocate(cl(1:ndata))
    allocate(dis(1:ndata,1:kk))
    allocate(w1(1:ndata))
    allocate(wxd(1:kk,1:ndata,1:mcol))
    allocate(wxm(1:kk,1:mcol))
    allocate(wvo(1:mcol,1:mcol))
    allocate(wwvo(1:kk,1:mcol,1:mcol))
    allocate(num(1:kk))
    allocate(wx(1:mcol))
    allocate(wm(1:mcol))
    allocate(dist(1:kk))

    do k=1,ndata
        read(11,*) (xd(k,j),j=1,mcol)
    end do
    close(11)
    !-------------------------------------------
    ! Normalization
    !-------------------------------------------
!    call STAND(ndata,mcol,xd)
    !-------------------------------------------
    ! Initial value of cluster (k-means++)
    !-------------------------------------------
    call CLINI(ndata,mcol,kk,imethod,xd,cl,ranseedflag)
    !-------------------------------------------
    ! Start of Iteration work
    !-------------------------------------------
    dists0=1.0D30
    iii=0
    do while(iii<=itmax*10)
        iii=iii+1
        call SELLDAT(ndata,mcol,kk,xd,cl,wxd,num)
        do k=1,kk
            do j=1,mcol
                wxm(k,j)=0.0D0
                do i=1,num(k)
                    wxm(k,j)=wxm(k,j)+wxd(k,i,j)
                end do
                wxm(k,j)=wxm(k,j)/dble(num(k))
            end do
        end do
        ! Preparation for calculation of mahalabinos distance
        if(imethod==0)then
            do k=1,kk
                allocate(wva(1:num(k),1:mcol))
                do j1=1,num(k)
                    do j2=1,mcol
                        wva(j1,j2)=wxd(k,j1,j2)
                    end do
                end do
                call VAR(num(k),mcol,wva,wvo)
                deallocate(wva)
                ! ----------------------------------------------------
                ! It is not certain                                   
                ! whether all the time inverse matrix is obtained.    
                ! ----------------------------------------------------
                call MATINV(mcol,wvo,mcol)
                do i=1,mcol
                    do j=1,mcol
                        wwvo(k,i,j)=wvo(i,j)
                    end do
                end do
            end do
        end if
        ! Calculation of distance
        do i=1,ndata
            do j=1,mcol
                wx(j)=xd(i,j)
            end do
            do k=1,kk
                do j=1,mcol
                    wm(j)=wxm(k,j)
                end do
                !------------------------------------
                ! imethod=0: Mahalanobis distance
                ! imethod=1: Cutyblock distance
                ! imethod=2: Euclidean distance
                !------------------------------------
                if(imethod==0)then
                    do j1=1,mcol
                        do j2=1,mcol
                            ! Inverse matrix of V-CV matrix
                            wvo(j1,j2)=wwvo(k,j1,j2)
                        end do
                    end do
                    dis(i,k)=MAHALANOBIS(mcol,wx,wm,wvo)
                else
                    dis(i,k)=MINKOWSKI(mcol,wx,wm,imethod)
                end if
            end do
            ! Classification
            dismin=dis(i,1)
            do k=2,kk
                if(dis(i,k)<dismin)dismin=dis(i,k)
            end do
            do k=1,kk
                if(dismin==dis(i,k))cl(i)=k
            end do
        end do
        ! Judgement of convergence
        do k=1,kk
            dist(k)=0.0D0
        end do
        do i=1,ndata
            do k=1,kk
                if(cl(i)==k)dist(k)=dist(k)+dis(i,k)
            end do
        end do
        dists1=0.0D0
        do k=1,kk
            dists1=dists1+dist(k)
        end do
write(6,'(i5,4(f10.3))') iii,dist(1),dist(2),dist(3),dists1
        if(abs(dists1-dists0)<eps)then
            nzero=1
            do k=1,kk
                if(num(k)<mds)nzero=0
            end do
            if(1<=nzero)then
                exit
            else
                write(6,*) "num:",(num(k),k=1,kk)
                if(itmax<=iii)exit
                call CLINI(ndata,mcol,kk,imethod,xd,cl,ranseedflag)
                call SELLDAT(ndata,mcol,kk,xd,cl,wxd,num)
                dists1=1.0D30
            end if
        end if
        dists0=dists1
    end do
    !-------------------------------------------
    ! End of Iteration work
    !-------------------------------------------

write(6,*) "***"
write(6,*) mcol,ndata,kk
write(6,*) num(1),num(2),num(3)
!do i=1,ndata
!    write(6,'(i5,3(f10.1),i5)') i,dis(i,1),dis(i,2),dis(i,3),cl(i)
!end do


    open(12, file=fnameW, status='replace')
    write(12,'(a)'),"v1 v2 v3 v4 cl"
    do i=1,ndata
        write(12,'(4(f10.1),i5)') (xd(i,j),j=1,mcol),cl(i)
    end do
    close(12)

stop

contains

    !**************************************************
    ! Normalization
    !**************************************************
    subroutine STAND(ndata,mcol,xd)
        integer,intent(in)::ndata,mcol
        real(8),intent(inout)::xd(1:ndata,1:mcol)

        integer::i,k
        real(8)::xm(1:mcol),sd(1:mcol)

        !Average
        do i=1,mcol
            xm(i)=0.0D0
            do k=1,ndata
                xm(i)=xm(i)+xd(k,i)
            end do
            xm(i)=xm(i)/dble(ndata)
        end do
        !Standard deviation
        do i=1,mcol
            sd(i)=0.0D0
            do k=1,ndata
                sd(i)=sd(i)+(xd(k,i)-xm(i))*(xd(k,i)-xm(i))
            end do
            sd(i)=sqrt(sd(i)/dble(ndata-1))
        end do
        !Normalization
        do i=1,mcol
            do k=1,ndata
                xd(k,i)=(xd(k,i)-xm(i))/sd(i)
            end do
        end do
    end subroutine STAND
    !**************************************************
    ! Sellect Data
    !**************************************************
    subroutine SELLDAT(ndata,mcol,kk,xd,cl,wxd,num)
        integer,intent(in)::ndata,mcol,kk
        integer,intent(in)::cl(1:ndata)
        real(8),intent(in)::xd(1:ndata,1:mcol)
        real(8),intent(out)::wxd(1:kk,1:ndata,1:mcol)
        integer,intent(out)::num(1:kk)

        integer::i,j,k

        do k=1,kk
            num(k)=0
            do i=1,ndata
                do j=1,mcol
                    wxd(k,i,j)=0.0D0
                end do
            end do
        end do
        do i=1,ndata
            do k=1,kk
                if(cl(i)==k)then
                    num(k)=num(k)+1
                    do j=1,mcol
                        wxd(k,num(k),j)=xd(i,j)
                    end do
                end if
            end do
        end do
    end subroutine SELLDAT
    !**************************************************
    ! Variance-covariance matrix
    !**************************************************
    subroutine VAR(nrow,ncol,xd,varm)
        integer,intent(in)::nrow,ncol
        real(8),intent(in)::xd(1:nrow,1:ncol)
        real(8),intent(out)::varm(1:ncol,1:ncol)
        integer::i,j,k,nn,mm
        real(8)::xm(1:ncol),wmat(1:nrow,1:ncol)

        nn=nrow
        mm=ncol
        do i=1,nn
            do j=1,mm
                wmat(i,j)=xd(i,j)
            end do
        end do
        do j=1,mm
            xm(j)=0.0D0
            do i=1,nn
                xm(j)=xm(j)+wmat(i,j)
            end do
            xm(j)=xm(j)/dble(nn)
        end do
        do i=1,mm
            do j=1,mm
                varm(i,j)=0.0D0
                do k=1,nn
                    varm(i,j)=varm(i,j)+(wmat(k,i)-xm(i))*(wmat(k,j)-xm(j))
                end do
                varm(i,j)=varm(i,j)/dble(nn-1)
            end do
        end do
    end subroutine VAR
    !**************************************************
    ! Minkowski distance
    !**************************************************
    real(8) function MINKOWSKI(mm,xx,xm,ip)
        integer,intent(in)::mm,ip
        real(8),intent(in)::xx(1:mm),xm(1:mm)
        integer::i
        real(8)::sum

        sum=0.0D0
        do i=1,mm
            sum=sum+abs(xx(i)-xm(i))**dble(ip)
        end do
        MINKOWSKI=(sum)**(1.0D0/dble(ip))
    end function MINKOWSKI
    !**************************************************
    ! Mahalanobis distance
    !**************************************************
    real(8) function MAHALANOBIS(mm,xx,xm,va)
        integer,intent(in)::mm
        real(8),intent(in)::xx(1:mm),xm(1:mm),va(1:mm,1:mm)
        integer::i,j,k
        real(8)::wmat(1:mm,1:mm),wv1(1:mm),wv2(1:mm),ws

        do i=1,mm
            do j=1,mm
                wmat(i,j)=va(i,j)
            end do
        end do
        do j=1,mm
            wv1(j)=xx(j)-xm(j)
        end do
        do i=1,mm
            wv2(i)=0.0D0
            do j=1,mm
                wv2(i)=wv2(i)+wmat(i,j)*wv1(j)
            end do
        end do
        ws=0.0D0
        do i=1,mm
            ws=ws+wv1(i)*wv2(i)
        end do
        MAHALANOBIS=sqrt(ws)
    end function MAHALANOBIS
    !**************************************************
    ! Inverse matrix
    !**************************************************
    subroutine MATINV(n,ary,nt)
        integer,intent(in)::n,nt
        real(8),intent(inout)::ary(1:nt,1:nt)

        integer::i,j,k,l,l1,irow,icol,jrow,jcol
        real(8)::aamax,s,t
        real(8)::factor=1.0D-30
        integer::ipivot(n)
        real(8)::pivot(n)
        integer::nindex(n,2)

        do j=1,n
            ipivot(j)=0
        end do
        do i=1,n
            aamax=0.0D0
            do j=1,n
                if(ipivot(j)-1/=0)then
                    do k=1,n
                        if(0<ipivot(k)-1)return
                        if(ipivot(k)-1<0)then
                            if(abs(aamax)-abs(ary(j,k))<0)then
                                irow=j
                                icol=k
                                aamax=ary(j,k)
                            end if
                        end if
                    end do
                end if
            end do
            ipivot(icol)=ipivot(icol)+1
            if(irow-icol/=0)then
                do l=1,n
                    s=ary(irow,l)
                    ary(irow,l)=ary(icol,l)
                    ary(icol,l)=s
                end do
            end if
            nindex(i,1)=irow
            nindex(i,2)=icol
            pivot(i)=ary(icol,icol)
            if(abs(pivot(i))<factor)then
                write(6,*) '#Singular'
                stop
            end if
            ary(icol,icol)=1.0D0
            do l=1,n
                ary(icol,l)=ary(icol,l)/pivot(i)
            end do
            do l1=1,n
                if(l1-icol/=0)then
                    t=ary(l1,icol)
                    ary(l1,icol)=0.0D0
                    do l=1,n
                        ary(l1,l)=ary(l1,l)-ary(icol,l)*t
                    end do
                end if
            end do
        end do
        do i=1,n
            l=n+1-i
            if(nindex(l,1)-nindex(l,2)/=0)then
                jrow=nindex(l,1)
                jcol=nindex(l,2)
                do k=1,n
                    s=ary(k,jrow)
                    ary(k,jrow)=ary(k,jcol)
                    ary(k,jcol)=s
                end do
            end if
        end do
    end subroutine MATINV
    !**************************************************
    ! Initial values (k-mean++)
    !**************************************************
    subroutine CLINI(ndata,mcol,kk,imethod,xd,cl,ranseedflag)
        integer,intent(in)::ndata,mcol,kk,imethod
        real(8),intent(in)::xd(1:ndata,1:mcol)
        integer,intent(out)::cl(1:ndata)
        integer,intent(inout)::ranseedflag

        real(8)::wmi(1:kk,1:mcol)
        real(8)::wx(1:mcol),wm(1:mcol)
        real(8)::dis1(1:ndata),dis2(1:ndata),disr(1:ndata)
        real(8)::dis(1:ndata,1:kk),dismin,dismax
        integer::i,j,k,ir,ip
        real(8)::work
        integer :: clock             ! for random_seed()
        integer :: nrand             ! for random_seed()
        integer,allocatable::seed(:) ! for random_seed()

        ip=imethod
        if(imethod==0)ip=2

        if(ranseedflag==0)then
            call random_seed(size=nrand)
            allocate(seed(nrand))
            call system_clock(count=clock)
            seed = clock
            call random_seed(put = seed)
            ranseedflag=1
        end if

        k=1
        call random_number(work)
        ir=1+int(work*dble(ndata))
        write(6,*) "ir=",ir
        do j=1,mcol
            wmi(k,j)=xd(ir,j)
        end do
        do i=1,ndata
            disr(i)=1.0D30
        end do
        do k=1,kk-1
            do i=1,ndata
                do j=1,mcol
                    wx(j)=xd(i,j)
                    wm(j)=wmi(k,j)
                end do
                dis1(i)=disr(i)
                dis2(i)=MINKOWSKI(mcol,wx,wm,ip)
            end do
            do i=1,ndata
                disr(i)=min(dis1(i),dis2(i))
            end do
            dismax=disr(1)
            do i=2,ndata
                if(dismax<disr(i))then
                    dismax=disr(i)
                    ir=i
                end if
            end do
            do j=1,mcol
                wmi(k+1,j)=xd(ir,j)
            end do
        end do

        do i=1,ndata
            do j=1,mcol
                wx(j)=xd(i,j)
            end do
            do k=1,kk
                do j=1,mcol
                    wm(j)=wmi(k,j)
                end do
                dis(i,k)=MINKOWSKI(mcol,wx,wm,ip)
            end do
            dismin=dis(i,1)
            do k=2,kk
                if(dis(i,k)<dismin)dismin=dis(i,k)
            end do
            do k=1,kk
                if(dismin==dis(i,k))cl(i)=k
            end do
        end do
    end subroutine CLINI
!**************************************************
end program f90_KMEANPP
!**************************************************