- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I have been playing with KMEANS for a while, but I wanted to try it in Fortran.
Anyway, here is a F90 version of the f77 version on the web. Rosettacode.org.
!**********************************************************************
! 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 ITER ! maximum iterations
REAL BIG ! arbitrary large number
PARAMETER (ITER = 1000, BIG = 1E33)
! local variables
INTEGER H ! count iterations
INTEGER I ! count points
INTEGER I1 ! point marked as initial center
INTEGER J ! count dimensions
INTEGER L ! count clusters
INTEGER L0 ! present cluster ID
INTEGER L1 ! new cluster ID
REAL BEST ! shortest distance to a center
REAL D2 ! squared distance
REAL TOT ! a total
REAL 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
!***********************************************************************
! test program (extra credit #1)
!***********************************************************************
PROGRAM TPEC1
IMPLICIT NONE
Integer, PARAMETER :: N = 30, P = 2, K = 6
REal, Parameter :: TWOPI = 6.2831853
INTEGER I, L, Z(N), IFAULT
REAL X(P,N), C(P,K), R, THETA, W, WORK(N)
! Begin
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
write(*,*)"A :: ",I, R,THETA
X(1,I) = R * COS(THETA) ! Cartesian coordinates
X(2,I) = R * SIN(THETA)
!write(*,*)I,X(1,I),X(2,I)
END DO
! Call subroutine
CALL KMPP (X, P, N, K, C, Z, WORK, IFAULT)
PRINT *, 'kmpp returns with error code ', IFAULT
! Print lists of points in each cluster
DO L = 1, K
PRINT *, 'Cluster ', L, ' contains points: '
10 FORMAT (I6, $)
20 FORMAT ()
DO I = 1, N
IF (Z(I) .EQ. L) PRINT 10, I
END DO
PRINT 20
END DO
! Write CSV file with Y-coordinates in different columns by cluster
OPEN (UNIT=1, FILE='tpec1.csv', STATUS='UNKNOWN', 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
it needs more work, but it is a neat little program.
Link Copied
8 Replies
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
n = RGBnode%nn
call RGBnode%fill(1,n )I have to declare an integer n(2) and set it to pass n RGBnode%nn into the subroutine, it will not pass directly - I just do not understand why I need to do this?
Line 194 is the relevant subroutine.
module rgbimage_m
implicit none
private
public :: rgbimage
INTEGER, PARAMETER :: dp = selected_real_kind(15, 307)
REAL (KIND=dp) :: g = 9.806, pi = 3.14159265D0
type rgbimage
!! usage
!! 1) init
!! 2a) fill_image
!! or
!! 2b) set_pixel
!! 3) normalize
!! 4) write
integer Colred(3), colgre(3), colblu(3), Colyel(3)
integer, dimension(:,:,:), allocatable :: input
integer, dimension(:,:,:), allocatable :: rgb
!! pixel arrays of rgb values
!! indices (i,j,k)
!! i: position x_i
!! j: position y_j
!! k=1: red, k=2: green, k=3: blue
integer :: n(2)
!! image dimensions: [height, width]
integer :: nn(2)
!! image input dimensions: [height, width]
contains
procedure :: init => rgbimage_init ! inits image
procedure :: fill_image => rgbimage_fill_image ! fill image with constant rgb value
procedure :: get_pixel => rgbimage_get_pixel ! gets one pixel
procedure :: normalize => rgbimage_normalize ! normalizes all pixels onto range [0, 255]
procedure :: set_pixel => rgbimage_set_pixel ! sets one pixel
procedure :: write => rgbimage_write ! outputs image to file
procedure :: fill => rgbimage_fill ! outputs image to file
procedure, private :: inside => rgbimage_inside
procedure, private :: valid => rgbimage_valid
end type
contains
subroutine rgbimage_init(this, height, width)
!! initialize image.
!! sets dimensions, allocates pixels and sets colors to 0.
class(rgbimage), intent(out) :: this
integer, intent(in) :: height, width
integer astat
integer, parameter :: ihpixf = 4096, jvpixf = 4096 ! pixel size, jvpixf
this%n = [height, width]
this%nn = [ihpixf,jvpixf]
allocate (this%rgb(height,width,3), source=0, STAT = astat)
if(astat == 0) then
write(*,*)"Allocated matrix"
else
stop "Failed to allocate matrix"
end if
allocate(this%input(3,ihpixf,jvpixf),STAT = astat)
if(astat == 0) then
write(*,*)"Allocated matrix"
else
stop "Failed to allocate matrix"
end if
this%Colgre(1) = 248
this%Colgre(2) = 248
this%Colgre(3) = 255
this%Colyel(1) = 248
this%Colyel(2) = 100
this%Colyel(3) = 100
this%Colblu(1) = 50
this%Colblu(2) = 56
this%Colblu(3) = 200
this%Colred(1) = 200
this%Colred(2) = 256
this%Colred(3) = 256
return
end subroutine
logical function rgbimage_valid(this, check_rgb_vals)
!! checks if the image has valid dimensions and optionally valid rgb values.
class(rgbimage), intent(in) :: this
logical, intent(in), optional :: check_rgb_vals
!! check if rgb values are in allowed range [0, 255]?
!! default: dont check
! always check that dimensions match
rgbimage_valid = ( all(this%n > 0) .and. &
& (size(this%rgb, dim=1) == this%n(1)) .and. &
& (size(this%rgb, dim=2) == this%n(2)) .and. &
& (size(this%rgb, dim=3) == 3) )
! optionally: check if rgb values are in allowed range
if (present(check_rgb_vals)) then
if (check_rgb_vals) rgbimage_valid = ( rgbimage_valid .and. &
& (all(this%rgb >= 0)) .and. &
& (all(this%rgb <= 255)) )
end if
end function
logical function rgbimage_inside(this, x, y)
!! checks if given coordinates are inside the image
class(rgbimage), intent(in) :: this
integer, intent(in) :: x, y
rgbimage_inside = ((x > 0) .and. (x <= this%n(1)) .and. (y > 0) .and. (y <= this%n(2)))
end function
subroutine rgbimage_set_pixel(this, x, y, rgb)
class(rgbimage), intent(inout) :: this
integer, intent(in) :: x, y
!! coordinates
integer, intent(in) :: rgb(3)
!! red, green, blue values
if (this%inside(x, y)) then
! use given data at first
this%rgb(x,y,:) = rgb
! check if given data was out of bounds
where (this%rgb(x,y,:) > 255)
this%rgb(x,y,:) = 255
elsewhere (this%rgb(x,y,:) < 0)
this%rgb(x,y,:) = 0
end where
end if
end subroutine
function rgbimage_get_pixel(this, x, y) result(rgb)
class(rgbimage), intent(in) :: this
integer, intent(in) :: x, y
!! coordinates
integer :: rgb(3)
!! red, green, blue values
if (this%inside(x, y)) then
rgb = this%rgb(x,y,:)
else
rgb = 0
end if
end function
subroutine rgbimage_normalize(this)
!! normalize colors to be in range [0, 255]
class(rgbimage), intent(inout) :: this
where (this%rgb(:,:,:) > 255)
this%rgb(:,:,:) = 255
elsewhere (this%rgb(:,:,:) < 0)
this%rgb(:,:,:) = 0
end where
end subroutine
subroutine rgbimage_fill_image(this, rgb)
!! fill whole image with given rgb values.
class(rgbimage), intent(inout) :: this
integer, intent(in) :: rgb(3)
!! red, green, blue values
integer :: i
if (this%valid()) then
do i = 1, 3
this%rgb(:,:,i) = rgb(i)
end do
end if
end subroutine
subroutine rgbimage_fill(this, i, n)
implicit none
class(rgbimage), intent(out) :: this
integer, intent(in) :: i, n(2)
integer j
do 100 j = 1, n(1)
write(*,*)"Here :: ",i,j
100 end do
return
end subroutine
subroutine rgbimage_write(this, fname)
class(rgbimage), intent(in) :: this
character(*), intent(in) :: fname
character*100 :: error
!! file path, e.g. "tmp/out.ppm"
integer :: iounit, ios, i,j,k
open (newunit=iounit, file=fname, iostat=ios, action='WRITE')
if (ios /= 0) then
write(error,100)fname
100 format( "Error opening file: ", A20)
stop 'Error opening file'
endif
! write header
write (iounit, '(A)') 'P6'
write (iounit, '(I0, A, I0)') this%n(1), " ", this%n(2)
write (iounit, '(A)') '255'
do i = 1, this%n(1)
do j = 1, this%n(2)
write (iounit, '(3A1)', advance='no') [(achar(this%rgb(i,j,k)), k=1,3)]
end do
end do
close (unit=iounit, iostat=ios)
if (ios /= 0) error stop "Error closing file"
end subroutine
! --------------------------------------------
!
! Notes
! o With a parameter ipixout set at 1, 2 or others,
! this subroutine will generate PPM-P6(binary), PPM-P3(text),
! or BMP(24bit depth without color table).
!
! o Some parts follow DEC-FORTRAN convention that had been defacto-standard long ago.
! Some compilers today may not accept if "ipixout" is set other than 2.
!
! o g77 (ver. 3.3.3) works for all three choices.
! o Intel compiler (ver. 9 or later) works for all three choices.
!
! --------------------------------------------
!
subroutine pixout(rgb,nframe)
implicit none
! interface arg.
integer ihpixf, jvpixf
parameter(ihpixf = 4096, jvpixf = 4096) ! pixel size, eacg must be multiple of 4, if BMP is chosen as output format.
character*1 rgb(3,ihpixf,jvpixf) ! RGB data array
integer nframe
! local
character*12 fnameout
integer i, j, k
integer itmp, icnt
character*14 frmtstr
character*54 headmsw
character*4 byt4
character*2 byt2
! choices
integer ipixout
parameter(ipixout = 3) ! 1 / 2 / other= PPM6, PPM3, BMP(24bit)
if (ipixout .EQ. 1) then
! PPM P6
write(fnameout,'(''smpl'',i3.3,''.ppm'')') nframe ! name of PPM file
open(unit=2,file=fnameout,status='unknown')
write(*,*) 'Now writing PPM (P6) file : ', fnameout
! header
write(2,'(''P6'', 2(1x,i4),'' 255 '',$)')ihpixf, jvpixf
! image data
itmp = ihpixf ! jvpixf ! 3
write(frmtstr,'(''('',i8.8,''A,$)'')') itmp ! make output "format"
write(2,fmt=frmtstr) (((rgb(k,i,j),k=1,3),i=1,ihpixf),j=jvpixf,1,-1)
close(2)
else if (ipixout .EQ. 2) then
! PPM P3 ! rather "safer" choice for many Fortran compiler(s).
write(fnameout,'(''smpl'',i3.3,''.ppm'')') nframe ! name of PPM file
open(unit=2,file=fnameout,status='unknown')
write(*,*) 'Now writing PPM (P3) file : ', fnameout
! header
write(2,'(A)') 'P3'
write(2,'(2(1x,i4),'' 255 '')') ihpixf, jvpixf
icnt = 0
! image data
do j = jvpixf, 1, -1 ! here, j (vertical address) runs from top to bottom.
do i = 1, ihpixf, 1
do k = 1, 3
itmp = ichar(rgb(k,i,j))
icnt = icnt + 4
if (icnt .LT. 60) then
write(2,fmt='(1x,i3,$)') itmp ! "$" is not standard.
else
write(2,fmt='(1x,i3)') itmp
icnt = 0
endif
enddo
enddo
enddo
write(2,'(A)') ' '
close(2)
else
! BMP (24bit depth)... this part works only when width is multiple of 4.
itmp = mod(ihpixf, 4)
if (itmp .NE. 0) then
write(*,*) 'width must be multiple of 4'
stop
endif
write(fnameout,'(''smpl'',i3.3,''.bmp'')') nframe ! name of BMP file
open(unit=2,file=fnameout,status='unknown')
write(*,*) 'Now writing BMP(24bit) file : ', fnameout
! header 1 (file header ; 1--14 byte)
headmsw( 1: 2) = 'BM' ! declaring this is BMP file
itmp = 54 + ihpixf * jvpixf * 3 ! total file size = header + data
call num2bit4(itmp,byt4)
headmsw( 3: 6) = byt4(1:4)
itmp = 0 ! may be 0
call num2bit2(itmp,byt2)
headmsw( 7:
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I haven't tried your code. This seems like a regression error (bug) where the UDT%var as argument is seen as UDT alone.
Try: call RGBnode%fill(1,(RGBnode%nn))
The extra parenthesis will construct a copy of the member variable (unfortunate as it is, but removes the necessity for a user defined temp.
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
John, et. al.
I must caution that that fix only applies to arguments that are intended to be changed/defined by the called procedure.
Jim Dempsey
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
In pseudo-code
Module AA
TYPE A
INTEGER N(2)
End TYPE
FUNCTION START(N)
set N(2)
end function
FUNCTION USE(N)
USE N for work
end function
END Module AA
MAIN PROGRAM
TYPE A == set new
CALL A%Start(A%N)
!! A%N has now two integer values
CALL %Use((A%N))
!! Do not change it -- it was set in Start.
End main
This is how it is working at the moment, is that the intended use of the (())
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
It has been fun to get the KMeans working in Intel Fortran and putting a drawing program with it, so I get straight plots.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Nice. Where I can get more information about these colors in rgb.f90? I want to make a color scale, from blue to red, 6~8 values in between. Thanks.
case(1) !Medium Violet Red
this%colors(i,1) = 199
this%colors(i,2) = 21
this%colors(i,3) = 133
this%list(i) = "MVP"
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
In the ZIP file are three solutions in Fortran. They work on VS 2019 - but PSXE 2020.4. I had to go back to get a working Fortran with the latest VS.
The KMEANS solution works on a simple set of FFT data. The data is enclosed.
RGB.f90 is enclosed in the KM++ solution folder. It has the 140 colors of the WEB Safe W3 coded into the routine. You call a color by number. ignore the RGB in the Kmeans it is an old version.
I enclose the web safe PDF that lists all the colors. I have the subroutine set up for 140 colors, but only 70 have the full RGB codes in place. If you see Pink it is a place holder. I have found it is much easier to stick to the color names in this list.
I include the code so you have samples of the calling routines.
This is not production code, it is for experimenting with data.
Reply
Topic Options
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page