!*************************************************************************
!*                                                                       *
!*     Subroutines for nonsmooth optimization based clusterwise linear   *
!*     regression software LMBM (last modified 11.9.2019).              *
!*                                                                       *
!*************************************************************************
!*
!*     Modules included:
!*
!*     clrmod       !
!*

MODULE clrmod       ! Subroutines for clusterwise linear regression software

    USE r_precision, ONLY : prec      ! Precision for reals.
    USE clrfunctionmod, ONLY : hauxfunc, hfunc, lisa
    IMPLICIT NONE

    ! MODULE clusteringmod includes the following subroutines.

    PUBLIC :: &
        distribution, &    ! Distribution over clusters.
        checkfinal, &      ! Final distribution over clusters.
        step1, &           ! Initialization of clusters.
        cleaning, &        ! Cleaning arrays.
        scaling, &
        rescaling
!    PRIVATE :: &

CONTAINS

    !==================================================================
    !     distribution over clusters
    !===================================================================

    SUBROUTINE distribution(x,f)

        USE param, ONLY : zero, large
        USE initclr, ONLY : &
            maxdim, &
            a4, &                           ! Auxiliary data matrix a4(nrecord,nft).
            dminim, &                       !
            nft, &                          ! Number of features in data, from user.
            nrecord, &                      ! Number of instances in data, from user.
            nrecord1, &                     ! Number of instances used in computions.
            nc, &                           ! Current number of clusters, loops from 1 to nclust.
            nel, &
            nob

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::      &
            f,                 &
            x(maxdim)

        ! Local variables
        REAL(KIND=prec) ::      &
            d1,d2,d3
        INTEGER ::  &
            i,i1,j,k

        f = zero
        DO i=1,nc
            nel(i)=0
        END DO

        DO k=1,nrecord1
            d1 = large !1.0E+12_prec ! pitäisikö olla large?
            DO i=1,nc
                CALL hauxfunc(x,i,k,d3)
                IF (d3 < d1) THEN
                    d1=d3
                    i1=i
                END IF
            END DO
            f=f+d1
            nel(i1)=nel(i1)+1
            nob(i1,nel(i1))=k
            dminim(k)=d1
            END DO

        RETURN

    END SUBROUTINE distribution


    !===================================================================
    !  final distribution over clusters
    !===================================================================
    SUBROUTINE checkfinal(x,f,f1,rmse1,rmse2,ce1,ce2,rarvo1,rarvo2)

        USE param, ONLY : zero, one, large
        USE initclr, ONLY : &
            maxdim, &                       ! Maximum number of variables for optimization.
            a, &                            ! Data matrix a(nrecord,nft), from input file.
            nft, &                          ! Number of features in data, from user.
            nrecord, &                      ! Number of instances in data, from user.
            nrecord1, &                     ! Number of instances used in computions.
            nrecord2, &                     ! Number of instances used in computions.
            objfunc, &                      ! Objective function used
            objc, &                         ! Regularization parameter C.
            nc, &                           ! Current number of clusters, loops from 1 to nclust.
            objkvad, &
            nel, &
            nob, &
            nl, &
            nlist, &
            nlisttrain

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::      &
            f,f1,                 &
            fapu,                &
            x(maxdim), &
            rmse1,rmse2, &
            ce1,ce2

        ! Local variables
        REAL(KIND=prec) ::      &
            pred(nrecord), &
            dist(nrecord), &
            rarvo1,rarvo2, &
            r1,r2, &
            fe1,fe2, &
            fextra, &
            e0_for1, e0_for2, osak, &
            rsumo1, rsumo2, rsumn1, rsumn2, &
            d2,d3,d4,d6, apu, &
            b, b1
        INTEGER ::  &
            i,i1,i2,i3,i4,j,k,ntotal,ncount

        ! Intrinsic Functions
        INTRINSIC SQRT,REAL

        rmse1 = zero
        rmse2 = zero
        ce1 = zero
        ce2 = zero
		rarvo1 = zero
		rarvo2 = zero
		

        r1 = zero
        loop1: DO i=1,nrecord           ! sum of observed values in training data
            DO i1=1,nl
                IF(i == nlist(i1)) CYCLE loop1
            END DO
            r1=r1+a(i,nft)
        END DO loop1

        r2= zero                        ! sum of observed values in test data 
        DO i=1,nl
            i1=nlist(i)
            r2=r2+a(i1,nft)
        END DO

        r1=r1/REAL(nrecord1,prec)                      ! mean of observed values in training data
        IF (nrecord2 > 0) r2 = r2/REAL(nrecord2,prec)  ! mean of observed values in test data

        fe1 = zero
        loop2: DO i=1,nrecord          ! sum_i (e_i-e_0)^2   training set
            DO i1=1,nl
                IF(i == nlist(i1)) CYCLE loop2
            END DO
            fe1=fe1+(a(i,nft)-r1)**2
        END DO loop2

        fe2 = zero
        DO i=1,nl   	! sum_i (e_i-e_0)^2   test set
            i1=nlist(i)
            fe2=fe2+(a(i1,nft)-r2)**2
        END DO

        f = zero
        fapu = zero
		e0_for1 = zero
		
        ! training set
        loop3: DO i=1,nrecord
            DO i1=1,nl
                IF(i == nlist(i1)) CYCLE loop3 ! skip the test set
            END DO
            
            d4 = large
            DO k=1,nc   ! the regression function to which a point belongs to
                CALL hfunc(x,k,i,d3)
                IF (d3 < d4) d4=d3  ! error to nearest regression function
            END DO
            f=f+d4
            
            d4 = large       ! needed to calculate RMSE and other measures used in cross-validation
            DO k=1,nc
                d2=x(k*nft)
                DO i1=1,nft-1
                    d2=d2+a(i,i1)*x(i1+(k-1)*nft)
                END DO
                d3 = (d2-a(i,nft))**2
                !d3 = ABS(d2-a(i,nft))
                IF (d3 < d4) THEN 
				    d4 = d3    ! error to the nearest regression function
					osak = d2  ! forecast value		
				END IF 
            END DO
            fapu=fapu+d4    
            e0_for1 = e0_for1 + osak    ! the sum of forecast value
			
        END DO loop3
		
		e0_for1 = e0_for1/REAL(nrecord1,prec)            ! the mean of forecast values

		rsumn1=0.0_prec
		rsumo1=0.0_prec		
        ! training set
        loop33: DO i=1,nrecord
            DO i1=1,nl
                IF(i == nlist(i1)) CYCLE loop33 ! skip the test set
            END DO
            
            d4 = large       ! needed to calculate r measure used in cross-validation
            DO k=1,nc
                d2=x(k*nft)
                DO i1=1,nft-1
                    d2=d2+a(i,i1)*x(i1+(k-1)*nft)
                END DO
                d3 = (d2-a(i,nft))**2
				
                IF (d3 < d4) THEN 
				    d4 = d3    ! error to the nearest regression function
					osak = (d2-e0_for1) ! bar{e}_i-bar{e}_0		
				END IF 
            END DO
               
            rsumo1 = rsumo1 + osak * (a(i,nft)-r1)    
			rsumn1 = rsumn1 + osak**2
			
        END DO loop33		

		CALL lisa(x,f)
        		             
        rmse1=SQRT(fapu/REAL(nrecord1,prec))  ! RMSE for training set
        ce1=one-fapu/fe1                      ! CE for test set
		rarvo1 = rsumo1/(SQRT(rsumn1*fe1))    ! r for test set

        f1 = zero
		fextra= zero
        ncount=0
		e0_for2 = zero
        ! test set
        DO i=1,nl 
            i1=nlist(i)
            d6 = large            ! pisteen i1 etäisyys training-joukkon pisteisiin
            DO k=1,nc
                DO i2=1,nel(k)    ! nel(k) tells the number of elements in cluster k
                    i3=nob(k,i2)  ! index of element i2 in cluster k
                    i4=nlisttrain(i3)
                    dist(i4)=zero
                    DO j=1,nft-1
                        dist(i4)=dist(i4)+(a(i1,j)-a(i4,j))**2 ! distance between points in training and test set
                    END DO
                    IF (dist(i4) < d6) d6=dist(i4) ! the smallest distance
                END DO
            END DO
			
            d6=1.5_prec*d6
            ntotal=0
            b=zero
            DO k=1,nc                ! keskiarvo training-joukon pisteiden output-arvoista, jotka on riittävän lähellä pistettä i1
                DO i2=1,nel(k)
                    i3=nob(k,i2)
                    i4=nlisttrain(i3)
                    IF (dist(i4) <= d6) THEN
                        b=b+a(i4,nft)
                        ntotal=ntotal+1
                    END IF
                END DO
            END DO
            b1=b/REAL(ntotal,prec)   ! pisteen i1 output-arvon arvio

            d4=large
            ncount=ncount+1
            DO k=1,nc
                d2=x(k*nft)
                DO j=1,nft-1
                    d2=d2+a(i1,j)*x(j+(k-1)*nft)
                END DO
                d3=(d2-a(i1,nft))**2
				
                !d3=(d2-b1)**2        !b1 on jonkinlainen keskiarvo pisteen i1 output-arvolle: keskiarvo training-joukon pisteiden outputeistä, jotka on lähellä pistettä i1
                !d3=ABS(d2-b1)
                IF (d4 >= d3) THEN
                    d4=d3
                    pred(ncount)=d2
					osak = d2
                END IF
            END DO
            f1=f1+d4
            e0_for2 = e0_for2 + osak
			
        END DO
		
        e0_for2 = e0_for2/REAL(nrecord2,prec)            ! the mean of forecast values	

		rsumn2=0.0_prec
		rsumo2=0.0_prec
		! test set
        loop33test: DO i=1,nl
            i1=nlist(i)
           
            d4 = large       ! needed to calculate RMSE and other measures used in cross-validation
            DO k=1,nc
                d2=x(k*nft)
                DO j=1,nft-1
                    d2=d2+a(i1,j)*x(j+(k-1)*nft)
                END DO
                d3=(d2-a(i1,nft))**2
				
                IF (d3 < d4) THEN 
				    d4 = d3             ! error to the nearest regression function
					osak = (d2-e0_for2) ! bar{e}_i-bar{e}_0		
				END IF 
            END DO   
            rsumo2 = rsumo2 + osak * (a(i1,nft)-r2)    ! the sum of forecast value
			rsumn2 = rsumn2 + osak**2
			
        END DO loop33test	
		
        IF (nrecord2 > 0) THEN
            ce2=one-f1/fe2                      ! CE for test set
            rmse2=SQRT(f1/REAL(nrecord2,prec))  ! RMSE for test set
		    rarvo2 = rsumo2/(SQRT(rsumn2*fe2))  ! r for test set			
        END IF
			
		
        WRITE(41,*)
        WRITE(41,52) nc
52      FORMAT(i10)
        WRITE(41,*)
        DO i=1,ncount
            WRITE(41,51) pred(i)
        END DO
51      FORMAT(f12.4)
        RETURN

    END SUBROUTINE checkfinal

    !===================================================================
    !  Initialize clusters
    !===================================================================

    SUBROUTINE step1(x,z,jpoints)

        USE param, ONLY : zero
        USE initclr, ONLY : &
            maxdim, &                       ! Maximum number of variables for optimization.
            maxinit, &                      ! Maximum number of candidate points of data set.
            a4, &                           ! Auxiliary data matrix a4(nrecord,nft).
            gamma1, &                       ! Parameter affecting the size of the set of starting point,
            dminim, &                       !
            nclust, &                       ! Maximum number of clusters, from user.
            nft, &                          ! Number of features in data, from user.
            nrecord, &                      ! Number of instances in data, from user.
            nrecord1, &                     ! Number of instances used in computions.
            nc, &                           ! Current number of clusters, loops from 1 to nclust.
            nel, &
            nob

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::      &
            x(maxdim), &
            z(maxinit)
        INTEGER :: &
            jpoints

        ! Local variables
        REAL(KIND=prec) ::      &
            y, &
            d0,d2,d3,d4,d6,d7, &
            d1(nclust,nrecord), &
            d5(nrecord)
        INTEGER ::  &
            l1(nrecord), &
            i,i1,j,k,k1


        DO i=1,nc-1
            DO k=1,nel(i)
                k1=nob(i,k)
                l1(k1)=i
            END DO
        END DO

        DO i=1,nc-1
            DO k=1,nrecord1
                d0=zero
                DO j=1,nft-1
                    d0=d0+x(j+(i-1)*nft)*a4(k,j)
                END DO
                d1(i,k)=d0
            END DO
        END DO

        d6=zero
        DO i=1,nrecord1
            i1=l1(i)

            y=a4(i,nft)-d1(i1,i)
            d5(i)=zero
            DO k=1,nrecord1
                d2=d1(i1,k)+y
                d3=(d2-a4(k,nft))**2
                d4=dmin1(zero,d3-dminim(k))
                d5(i)=d5(i)+d4
            END DO
            IF (d5(i) <= d6) THEN
                d6=d5(i)
            END IF
        END DO

        jpoints=0
        d7=gamma1*d6
        DO k=1,nrecord1
            i1=l1(k)

            IF (d5(k) <= d7) THEN
                jpoints=jpoints+1
                DO j=1,nft-1
                    z(j+(jpoints-1)*nft)=x(j+(i1-1)*nft)
                END DO
                z(jpoints*nft)=a4(k,nft)-d1(i1,k)

            END IF

        END DO

        CALL cleaning(jpoints,z)

        RETURN

    END SUBROUTINE step1


    !===================================================================
    !  Cleaning arrays
    !===================================================================

    SUBROUTINE cleaning(jpoints,z)

        USE param, ONLY : zero
        USE initclr, ONLY : &
            maxinit, &                      ! Maximum number of candidate points of data set.
            toler, &                        !
            nft                             ! Number of features in data, from user.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::      &
            z(maxinit)
        INTEGER :: &
            jpoints

        ! Local variables
        REAL(KIND=prec) ::      &
            z1(nft), &
            z2(maxinit), &
            d2
        INTEGER ::  &
            i,j,jpoints0,k

        ! Intrinsic Functions
        INTRINSIC ABS

        jpoints0=0
        loop1: DO i=1,jpoints
            DO j=1,nft
                z1(j)=z(j+(i-1)*nft)
            END DO
            IF (i == 1) THEN
                jpoints0=1
                DO j=1,nft
                    z2(j)=z1(j)
                END DO

            ELSE
                DO k=1,jpoints0
                    d2=zero
                    DO j=1,nft
                        d2=d2+ABS(z1(j)-z2(j+(k-1)*nft))
                    END DO
                    IF(d2 <= toler) CYCLE loop1
                END DO
                jpoints0=jpoints0+1
                DO j=1,nft
                    z2(j+(jpoints0-1)*nft)=z1(j)
                END DO
            END IF
        END DO loop1

        DO k=1,jpoints0
            DO j=1,nft
                z(j+(k-1)*nft)=z2(j+(k-1)*nft)
            END DO
        END DO
        jpoints=jpoints0

        RETURN

    END SUBROUTINE cleaning

    !===================================================================
    !  Scaling
    !===================================================================

    SUBROUTINE scaling

        USE param, ONLY : zero, one
        USE initclr, ONLY : &
            a4, &                           ! Auxiliary data matrix a4(nrecord,nft).
            kvad, &                         ! Scaling vector to the norm of x
            nft, &                          ! Number of features in data, from user.
            nrecord, &                      ! Number of instances in data, from user.
            clin, &
            dlin

        IMPLICIT NONE

        ! Local variables
        REAL(KIND=prec) ::      &
            dm,var
        INTEGER ::  &
            i,j

        ! Intrinsic Functions
        INTRINSIC ABS,SQRT

        DO i=1,nft
            dm=zero
            DO j=1,nrecord
                dm=dm+a4(j,i)
            END DO
            dm=dm/REAL(nrecord,prec)
            var=zero
            DO j=1,nrecord
                var=var+(a4(j,i)-dm)**2
            END DO
            var=var/REAL(nrecord-1,prec)
            var=SQRT(var)

            IF (var >= 1.0E-05_prec) THEN
                clin(i)=one/var
                dlin(i)=-dm/var
                DO j=1,nrecord
                    a4(j,i)=clin(i)*a4(j,i)+dlin(i)
                END DO
            ELSE
                IF (ABS(dm) <= one) THEN
                    clin(i)=one
                    dlin(i)=zero
                ELSE
                    clin(i)=one/dm
                    dlin(i)=zero
                END IF
                DO j=1,nrecord
                    a4(j,i)=clin(i)*a4(j,i)+dlin(i)
                END DO
            END IF
        END DO
		
		DO i =1,nft-1
		   kvad(i) = clin(i)**2/clin(nft)
		END DO
		
		!WRITE(*,*) "1/var", clin
		!WRITE(*,*) "mu/var", dlin

        RETURN

    END SUBROUTINE scaling


    !===================================================================
    !  Rescaling
    !===================================================================

    SUBROUTINE rescaling(x,x4)

        USE param, ONLY : zero, one, large
        USE initclr, ONLY : &
            maxdim, &                       ! Maximum number of variables for optimization.
            !maxinit, &                      ! Maximum number of candidate points of data set.
            !a, &                            ! Data matrix a(nrecord,nft), from input file.
            !a4, &                           ! Auxiliary data matrix a4(nrecord,nft).
            !gamma1, &                       ! Parameter affecting the size of the set of starting point,
            !       gamma2, &                       ! Parameter affecting the size of the set of starting point,
            !       gamma3, &                       ! Parameter affecting the size of the set of starting point,
            !       gamma4, &                       ! Parameter affecting the size of the set of starting point,
            !dminim, &                       !
            !toler, &                        !
            !       tlimit, &                       ! Time limit, from user.
            !       noutcom, &                      ! Number of the output column, form user.
            !       itest, &                        ! Selection of test set usage:
                                            !   1: no test set,
                                            !   2: one training set and one test set,
                                            !   3: cross-validation.
            !        ntrain, &                       ! Number of points in training set.
            !        nfoldmax, &                     ! Number of cross validations.
            !        iscale, &                       ! Selection of scaling:
                                            !   1: scaling,
                                            !   2: no scaling.
            !        nupdate, &                      ! Number of updates.
            !        naux, &                         ! Number of auxiliary function evaluations.
            !        nfit, &                         ! Number of function evaluations.
            !nclust, &                       ! Maximum number of clusters, from user.
            nft, &                          ! Number of features in data, from user.
        !nrecord, &                      ! Number of instances in data, from user.
            clin, &
            dlin, &
        !nrecord1, &                     ! Number of instances used in computions.
        !nrecord2, &                     ! Number of instances used in computions.
        nc                           ! Current number of clusters, loops from 1 to nclust.
        !        ns, &                           ! Switch for auxiliary and real clustering problem.
        !nel, &
        !nob, &
            !nl, &
            !nlist, &
            !nlisttrain
        !        init_clrpar, &                  ! Furher initialization of parameters.
        !        def_clrpar                      ! Default values of clustering parameters.

        IMPLICIT NONE

        ! Subroutine arguments
        REAL(KIND=prec) ::      &
            x(maxdim),x4(maxdim)

        ! Local variables
        REAL(KIND=prec) ::      &
             d0,rabs
        INTEGER ::  &
            i,j

        ! Intrinsic Functions
        INTRINSIC ABS

        rabs=ABS(clin(nft))
        if (rabs > 1.0E-08_prec) THEN
            DO i=1,nc
                d0=zero
                DO j=1,nft-1
                    x4(j+(i-1)*nft)=x(j+(i-1)*nft)*clin(j)/rabs
                    d0=d0+x(j+(i-1)*nft)*dlin(j)
                END DO
                x4(i*nft)=(x(i*nft)+d0-dlin(nft))/rabs
            END DO
        END IF

        RETURN

    END SUBROUTINE rescaling



END MODULE clrmod
