PROGRAM TPEC1 IMPLICIT NONE REAL, Parameter :: TWOPI=6.2831853 !Integer, Parameter :: N=30000, P=2, K=6 !INTEGER :: I, L, Z(N), IFAULT !REAL :: X(P,N), C(P,K), R, THETA, W, WORK(N) Integer :: N, P, K, I, L, IFAULT Real :: R, THETA, W Integer, Allocatable :: Z(:) Real, Allocatable :: X(:,:), C(:,:), Work(:) Integer :: j !------------------------------------------------------------------------- K=6 open(unit = 10, file = "Dr.txt", Status="Old") read (10, *) N, P Allocate (X(P,N), C(P,K), WORK(N), Z(N)) Do i=1, N read (10,*) (X(j,i), j=1,P) Enddo !CALL RANDOM_SEED() !DO I = 1, N ! random points over unit circle ! CALL RANDOM_NUMBER (W) ! R = SQRT(W) ! radius ! CALL RANDOM_NUMBER (W) ! THETA = W * TWOPI ! angle ! X(1,I) = R * COS(THETA) ! Cartesian coordinates ! X(2,I) = R * SIN(THETA) !END DO CALL KMPP (X, P, N, K, C, Z, WORK, IFAULT) PRINT *, "kmpp returns with error code ", IFAULT DO L = 1, K PRINT *, "Cluster ", L, " contains points: " DO I = 1, N IF (Z(I) .EQ. L) PRINT "(I6)", I END DO PRINT* END DO !------------------------------------------------------------------------- OPEN (UNIT=1, FILE="tpec1.csv", STATUS="NEW", IOSTAT=IFAULT) IF (IFAULT .NE. 0) PRINT *, "tpec1: trouble opening file" 30 FORMAT (F8.4, $) 40 FORMAT (',', $) 50 FORMAT (F8.4) DO I = 1, N WRITE (UNIT=1, FMT=30, IOSTAT=IFAULT) X(1,I) IF (IFAULT .NE. 0) PRINT *, "tpec1: trouble writing X-coord" DO L = 1, Z(I) ! one comma per cluster ID WRITE (UNIT=1, FMT=40, IOSTAT=IFAULT) IF (IFAULT .NE. 0) PRINT *, "tpec1: trouble writing comma" END DO WRITE (UNIT=1, FMT=50, IOSTAT=IFAULT) X(2,I) IF (IFAULT .NE. 0) PRINT *, "tpec1: trouble writing Y-coord" END DO ! Write the centroids in the far column DO L = 1, K WRITE (UNIT=1, FMT=30, IOSTAT=IFAULT) C(1,L) IF (IFAULT .NE. 0) PRINT *, "tpec1: trouble writing X-coord" DO I = 1, K+1 WRITE (UNIT=1, FMT=40, IOSTAT=IFAULT) IF (IFAULT .NE. 0) PRINT *, "tpec1: trouble writing comma" END DO WRITE (UNIT=1, FMT=50, IOSTAT=IFAULT) C(2,L) IF (IFAULT .NE. 0) PRINT *, "tpec1: trouble writing Y-coord" END DO CLOSE (UNIT=1) END ! of test program !*********************************************************************** !* KMPP - K-Means++ - Traditional data clustering with a special initialization !* Public Domain - This program may be used by any person for any purpose. !* !* Origin: !* Hugo Steinhaus, 1956 !* !* Refer to: !* "kmeans++: the advantages of careful seeding" !* David Arthur and Sergei Vassilvitskii !* Proceedings of the eighteenth annual ACM-SIAM symposium !* on Discrete algorithms, 2007 !* !*____Variable_______I/O_______Description___________________Type_______ !* X(P,N) In Data points Real !* P In Dimension of the data Integer !* N In Number of points Integer !* K In # clusters Integer !* C(P,K) Out Center points of clusters Real !* Z(N) Out What cluster a point is in Integer !* WORK(N) Neither Real !* IFAULT Out Error code Integer !************************************************************************ SUBROUTINE KMPP (X, P, N, K, C, Z, WORK, IFAULT) IMPLICIT NONE INTEGER :: P, N, K, Z, IFAULT REAL :: X, C, WORK DIMENSION :: X(P,N), C(P,K), Z(N), WORK(N) !* constants INTEGER, PARAMETER :: ITER = 1000 ! maximum iterations REAL, PARAMETER :: BIG = 1E33 ! arbitrary large number !* local variables Integer :: H, & ! count iterations I, & ! count points I1, & ! point marked as initial center J, & ! count dimensions L, & ! count clusters L0, & ! present cluster ID L1 ! new cluster ID REAL :: BEST, & ! shortest distance to a center D2, & ! squared distance TOT, & ! a total W ! temp scalar LOGICAL :: CHANGE ! whether any points have been reassigned !************************************************************************ !* Begin. !************************************************************************ IFAULT = 0 IF (K < 1 .OR. K > N) THEN ! K out of bounds IFAULT = 3 RETURN END IF DO I = 1, N ! clear Z Z(I) = 0 END DO !************************************************************************ !* initial centers !************************************************************************ DO I = 1, N WORK(I) = BIG END DO CALL RANDOM_NUMBER (W) I1 = MIN(INT(W * FLOAT(N)) + 1, N) ! choose first center at random DO J = 1, P C(J,1) = X(J,I1) END DO DO L = 2, K ! initialize other centers TOT = 0. DO I = 1, N ! measure from each point BEST = WORK(I) D2 = 0. ! to prior center DO J = 1, P D2 = D2 + (X(J,I) - C(J,L-1)) **2 ! Squared Euclidean distance IF (D2 .GE. BEST) GO TO 10 ! needless to add to D2 END DO ! next J IF (D2 < BEST) BEST = D2 ! shortest squared distance WORK(I) = BEST 10 TOT = TOT + BEST ! cumulative squared distance END DO ! next data point !************************************************************************ !* Choose center with probability proportional to its squared distance !* from existing centers. !************************************************************************ CALL RANDOM_NUMBER (W) W = W * TOT ! uniform at random over cumulative distance TOT = 0. DO I = 1, N I1 = I TOT = TOT + WORK(I) IF (TOT > W) GO TO 20 END DO ! next I 20 CONTINUE DO J = 1, P ! assign center C(J,L) = X(J,I1) END DO END DO ! next center to initialize !************************************************************************ !* main loop !************************************************************************ DO H = 1, ITER CHANGE = .FALSE. !* find nearest center for each point DO I = 1, N L0 = Z(I) L1 = 0 BEST = BIG DO L = 1, K D2 = 0. DO J = 1, P D2 = D2 + (X(J,I) - C(J,L)) **2 IF (D2 .GE. BEST) GO TO 30 END DO 30 CONTINUE IF (D2 < BEST) THEN ! new nearest center BEST = D2 L1 = L END IF END DO ! next L IF (L0 .NE. L1) THEN Z(I) = L1 ! reassign point CHANGE = .TRUE. END IF END DO ! next I IF (.NOT. CHANGE) RETURN ! success !************************************************************************ !* find cluster centers !************************************************************************ DO L = 1, K ! zero population WORK(L) = 0. END DO DO L = 1, K ! zero centers DO J = 1, P C(J,L) = 0. END DO END DO DO I = 1, N L = Z(I) WORK(L) = WORK(L) + 1. ! count DO J = 1, P C(J,L) = C(J,L) + X(J,I) ! add END DO END DO DO L = 1, K IF (WORK(L) < 0.5) THEN ! empty cluster check IFAULT = 1 ! fatal error RETURN END IF W = 1. / WORK(L) DO J = 1, P C(J,L) = C(J,L) * W ! multiplication is faster than division END DO END DO END DO ! next H IFAULT = 2 ! too many iterations RETURN END ! of KMPP !************************************************************************