1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
Subroutine Orderdata(mat, n, k, rowscore)
C23456789012345678901234567890123456789012345678901234567890123456789012
C d = a distance matrix
C n = number of objects (rows of the data files)
C
C Based on Pierre Legendre Fortran code and Table 9.5 and 9.10
C We compute the principal coordinate of the first axis only.
C Set the precision level for eigenvalue estimation
Integer mat(n,k)
double precision sumrow(n), sumtot
double precision rowscore(n),colscore(n),toler,epsilon
epsilon=0.000001
toler= 0.000001
if(n.gt.1000) then
epsilon=0.00001
toler= 0.00001
endif
if(n.gt.10000) then
epsilon=0.0001
toler= 0.0001
endif
C Centre the distance matrix (Gower's method)
call Centre(mat, n, k, sumrow, sumtot)
call TWWS(mat,n,k,sumrow,sumtot,rowscore,colscore,toler,epsilon)
end
Subroutine Centre(mat, n, k, sumrow, sumtot)
Integer mat(n,k)
double precision d
double precision sumrow(n),sumtot
do i=1,n
sumrow(i)=0.0
enddo
do i=1,(n-1)
do j=(i+1),n
call SM(mat, n, k, i, j, d)
d = -0.5*d**2
sumrow(i)=sumrow(i) + d
sumrow(j)=sumrow(j) + d
enddo
enddo
C
fln=1.0/float(n)
sumtot=0.
do i=1,n
sumtot=sumtot+sumrow(i)
sumrow(i)=sumrow(i)*fln
enddo
sumtot=sumtot/float(n**2)
C do 16 i=1,n
C do 16 j=1,n
C 16 d(i,j)=d(i,j)-sumrow(i)-sumrow(j)+sumtot
end
Subroutine SM(mat, n, k, i, j, d)
C Compute a simple matching coefficient from a table of K-means results (integers).
C The 'n' rows are the objects; the 'k' columns are the partitions.
Integer mat(n,k)
double precision d
C
a = 0.0
do kk=1,k
if(mat(i,kk).eq.mat(j,kk)) a=a+1
enddo
d = 1.0 - (a/float(k))
end
Subroutine TWWS(mat,n,k,sumrow,sumtot,rowscore,colscore,
+ toler,epsilon)
Integer n, niter
Integer mat(n,k)
double precision sumrow(n),sumtot,d
double precision rowscore(n),colscore(n),epsilon,oldS,newS,toler,
+ oldrowsc(n)
niter=1000
C Step 2: Take the column order as arbitrary initial site scores
do 4 i=1,n
4 rowscore(i)=dfloat(i)
oldS=0.
C Iterations starting
do 20 it=1,niter
C Step 3: Calculate new variable scores (equation 5.8, p. 119)
do 6 i=1,n
6 colscore(i)=rowscore(i)
C Step 4: Calculate new site scores (equation 5.9, p. 122)
do 8 i=1,n
rowscore(i)=0.
do 8 j=1,n
call SM(mat, n, k, i, j, d)
d = -0.5*d**2
d = d-sumrow(i)-sumrow(j)+sumtot
8 rowscore(i)=rowscore(i)+d*colscore(j)
C Step 6: Normalize the site scores
call NormTWWS(rowscore,n,newS)
if(newS.lt.epsilon) then
C write(*,103) 0
goto 52
endif
C When convergence has been attained, check sign of eigenvalue.
C If ALL rowscores have changed sign during the last iteration,
C this is a negative eigenvalue.
C write(*,*) 'oldS-newS', oldS-newS
if(dabs(oldS-newS).le.toler) then
C write(*,105) toler, it
goto 22
endif
do 18 i=1,n
18 oldrowsc(i)=rowscore(i)
oldS=newS
20 continue
C End of iterations for estimating eigenvalues and eigenvectors
C write(*,101) k
22 continue
C Save Eigenvalues, % variance, cumulative % variance
C write(*,*) 'Eigenvalue =', newS
C End of main loop on axes
52 continue
C
C Normalize the principal coordinates to variance = eigenvalue
do 60 i=1,n
60 rowscore(i)=rowscore(i)*dsqrt(newS)
C write(*,*) rowscore
C 101 format(' Convergence not reached for axis:',i3/
C + ' Increase NITER or lower TOLER')
C 102 format(' N. iterations to reach convergence for axis',i3,' =',i4)
C 103 format(' There are',i4,' eigenvalues different from 0')
C 104 format(' Eigenvector',i3,' is complex [multiply values*Sqrt(-1)]')
C 105 format(" Tolerance is: ", F12.8, " NIter is: ", i4)
C 150 format(1x,1hR,2i3,6f10.5)
C 151 format(1x,1hC,2i3,10x,6f10.5)
end
Subroutine NormTWWS(rowscore,n,newS)
Integer n
double precision rowscore(n),s2,newS
C Normalization for two-way weighted summation algorithm for PCA
C (ter Braak 1987: 123)
C On output, vector 'rowscore' has length = 1.
C S = eigenvalue*(p-1) = contraction of vector rowscore in final
C iteration
C
s2=0.0
do 10 i=1,n
10 s2=s2+rowscore(i)**2
newS=dsqrt(s2)
do 20 i=1,n
20 rowscore(i)=rowscore(i)/newS
return
end
|