-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsmacof.f90
131 lines (103 loc) · 3.37 KB
/
smacof.f90
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
subroutine smacof_embed(numvertices, target_dimension, weights, lengths, Vplus, maximum_iterations, X)
implicit none
! input
! dimension
integer, intent(in) :: numvertices
! iterations cap
integer, intent(in) :: maximum_iterations
! target dimension
integer, intent(in) :: target_dimension
! the edge lengths
real, dimension (numvertices, numvertices), intent(in) :: lengths
! the edge weights
real, dimension (numvertices, numvertices), intent(in) :: weights
! storage during calculation and output
integer step
! V_+ is th emoore-penrose inverse of the V matrix (the 'hessian') that holds the transformed edges, weights
! calculated (once) using numpy for convenience
real, dimension (numvertices, numvertices) :: Vplus
real :: stress
! previous step
real, dimension (numvertices, target_dimension) :: Z
! B = B(Z)
real, dimension (numvertices, numvertices) :: B
! output
! array of resulting coordinates
real, dimension (numvertices, target_dimension), intent(out) :: X
call RANDOM_NUMBER(X)
do step=1,maximum_iterations
Z = X
call BTransform(numvertices, target_dimension, weights, lengths, Z, B)
X = matmul(Vplus, matmul(B, Z))
call FindStress(numvertices, target_dimension, weights, lengths, X, stress)
print*, "step: ",step,"stress: ",stress
end do
end subroutine smacof_embed
subroutine FindStress(numvertices, target_dimension, weights, lengths, Z, stress)
implicit none
! input
integer, intent(in) :: target_dimension
integer, intent(in) :: numvertices
real, intent(in), dimension (numvertices, target_dimension) :: Z
real, dimension (numvertices, numvertices), intent(in) :: lengths
real, dimension (numvertices, numvertices), intent(in) :: weights
! intermediate
real :: distance
integer i, j, k
! output
real, intent(out) :: stress
stress = 0.0
! i != j
do j=2,numvertices
do i=1,j-1
distance = 0.0
do k=1,target_dimension
distance = distance + (Z(i, k) - Z(j,k))**2
end do
distance = sqrt(distance)
stress = stress + (distance - lengths(i, j))**2 * weights(i, j)
end do
end do
end subroutine FindStress
subroutine BTransform(numvertices, target_dimension, weights, lengths, Z, B)
implicit none
! input
integer, intent(in) :: target_dimension
integer, intent(in) :: numvertices
real, intent(in), dimension (numvertices, target_dimension) :: Z
real, dimension (numvertices, numvertices), intent(in) :: lengths
real, dimension (numvertices, numvertices), intent(in) :: weights
! intermediate
real :: distance
integer i, j, k
! output
real, dimension (numvertices, numvertices), intent(out) :: B
! initialize
B = 0.0
! i != j
do j=2,numvertices
do i=1,j-1
distance = 0.0
do k=1,target_dimension
distance = distance + (Z(i, k) - Z(j,k))**2
end do
distance = sqrt(distance)
if (distance > 0) then
B(i, j) = - weights(i, j) * lengths(i, j) / distance
B(j, i) = - weights(j, i) * lengths(j, i) / distance
else
B(i, j) = 0
B(j, i) = 0
end if
end do
end do
! i = j
do i=1, numvertices
do j=1, i-1
B(i, i) = B(i, i) - B(i, j)
end do
do j=i+1, numvertices
B(i, i) = B(i, i) - B(i, j)
end do
end do
end subroutine BTransform