Skip to content

Commit 119c7b1

Browse files
jeciscjecisc
committed
Cut dependency to Math-Core-Process
1 parent d89ee57 commit 119c7b1

File tree

2 files changed

+215
-134
lines changed

2 files changed

+215
-134
lines changed

‎src/Math-Matrix/PMJacobiTransformation.class.st‎

Lines changed: 146 additions & 115 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
"
2-
[[[
2+
```
33
| m jacobi eigenvalues eigenvectors |
44
m := PMSymmetricMatrix rows: #((84 -79 58 55)
55
(-79 84 -55 -58)
@@ -8,12 +8,16 @@ m := PMSymmetricMatrix rows: #((84 -79 58 55)
88
jacobi := PMJacobiTransformation matrix: m.
99
eigenvalues := jacobi evaluate.
1010
eigenvectors := jacobi transform columnsCollect: [ :each | each].
11-
]]]
11+
```
1212
"
1313
Class {
1414
#name : #PMJacobiTransformation,
15-
#superclass : #PMIterativeProcess,
15+
#superclass : #Object,
1616
#instVars : [
17+
'precision',
18+
'desiredPrecision',
19+
'maximumIterations',
20+
'result',
1721
'lowerRows',
1822
'transform'
1923
],
@@ -32,28 +36,50 @@ PMJacobiTransformation class >> new [
3236
^ self error: 'Illegal creation message for this class'
3337
]
3438

39+
{ #category : #accessing }
40+
PMJacobiTransformation >> desiredPrecision: aNumber [
41+
"Defines the desired precision for the result."
42+
43+
aNumber > 0 ifFalse: [ ^ self error: 'Illegal precision: ' , aNumber printString ].
44+
desiredPrecision := aNumber
45+
]
46+
47+
{ #category : #operation }
48+
PMJacobiTransformation >> evaluate [
49+
"Perform the iteration until either the desired precision is attained or the number of iterations exceeds the maximum."
50+
51+
| iterations |
52+
iterations := 0.
53+
[
54+
iterations := iterations + 1.
55+
precision := self evaluateIteration.
56+
self hasConverged or: [ iterations >= maximumIterations ] ] whileFalse.
57+
self finalizeIterations.
58+
^ result
59+
]
60+
3561
{ #category : #operation }
3662
PMJacobiTransformation >> evaluateIteration [
63+
3764
| indices |
3865
indices := self largestOffDiagonalIndices.
39-
self transformAt: ( indices at: 1) and: ( indices at: 2).
40-
^precision
66+
self transformAt: (indices at: 1) and: (indices at: 2).
67+
^ precision
4168
]
4269

4370
{ #category : #transformation }
4471
PMJacobiTransformation >> exchangeAt: anInteger [
45-
"Private"
72+
"Private"
73+
4674
| temp n |
4775
n := anInteger + 1.
4876
temp := result at: n.
49-
result at: n put: ( result at: anInteger).
77+
result at: n put: (result at: anInteger).
5078
result at: anInteger put: temp.
51-
transform do:
52-
[ :each |
53-
temp := each at: n.
54-
each at: n put: ( each at: anInteger).
55-
each at: anInteger put: temp.
56-
].
79+
transform do: [ :each |
80+
temp := each at: n.
81+
each at: n put: (each at: anInteger).
82+
each at: anInteger put: temp ]
5783
]
5884

5985
{ #category : #operation }
@@ -63,131 +89,136 @@ PMJacobiTransformation >> finalizeIterations [
6389

6490
| n |
6591
n := 0.
66-
result := lowerRows collect:
67-
[:each |
68-
n := n + 1.
69-
each at: n].
92+
result := lowerRows collect: [ :each |
93+
n := n + 1.
94+
each at: n ].
7095
self sortEigenValues
7196
]
7297

98+
{ #category : #testing }
99+
PMJacobiTransformation >> hasConverged [
100+
101+
^ precision <= desiredPrecision
102+
]
103+
104+
{ #category : #initialization }
105+
PMJacobiTransformation >> initialize [
106+
107+
super initialize.
108+
desiredPrecision := Float machineEpsilon.
109+
maximumIterations := 50
110+
]
111+
73112
{ #category : #initialization }
74113
PMJacobiTransformation >> initialize: aSymmetricMatrix [
75-
"Private"
76-
| n |
77-
n := aSymmetricMatrix numberOfRows.
78-
lowerRows := Array new: n.
79-
transform := Array new: n.
80-
1 to: n do:
81-
[ :k |
82-
lowerRows at: k put: ( ( aSymmetricMatrix rowAt: k) copyFrom: 1 to: k).
83-
transform at: k put: ( ( Array new: n) atAllPut: 0; at: k put: 1; yourself).
84-
]
114+
"Private"
115+
116+
| numberOfRows |
117+
numberOfRows := aSymmetricMatrix numberOfRows.
118+
lowerRows := Array new: numberOfRows.
119+
transform := Array new: numberOfRows.
120+
1 to: numberOfRows do: [ :k |
121+
lowerRows at: k put: ((aSymmetricMatrix rowAt: k) copyFrom: 1 to: k).
122+
transform at: k put: ((Array new: numberOfRows)
123+
atAllPut: 0;
124+
at: k put: 1;
125+
yourself) ]
85126
]
86127

87-
{ #category : #information }
128+
{ #category : #accessing }
88129
PMJacobiTransformation >> largestOffDiagonalIndices [
89-
"Private"
130+
"Private"
131+
90132
| n m abs |
91133
n := 2.
92134
m := 1.
93-
precision := ( ( lowerRows at: n) at: m) abs.
94-
1 to: lowerRows size do:
95-
[ :i |
96-
1 to: ( i - 1) do:
97-
[ :j |
98-
abs := ( ( lowerRows at: i) at: j) abs.
99-
abs > precision
100-
ifTrue: [ n := i.
101-
m := j.
102-
precision := abs.
103-
].
104-
].
105-
].
106-
^Array with: m with: n
107-
]
108-
109-
{ #category : #display }
135+
precision := ((lowerRows at: n) at: m) abs.
136+
1 to: lowerRows size do: [ :i |
137+
1 to: i - 1 do: [ :j |
138+
abs := ((lowerRows at: i) at: j) abs.
139+
abs > precision ifTrue: [
140+
n := i.
141+
m := j.
142+
precision := abs ] ] ].
143+
^ Array with: m with: n
144+
]
145+
146+
{ #category : #accessing }
147+
PMJacobiTransformation >> maximumIterations: anInteger [
148+
"Defines the maximum number of iterations."
149+
150+
(anInteger isInteger and: [ anInteger > 1 ]) ifFalse: [ ^ self error: 'Invalid maximum number of iteration: ' , anInteger printString ].
151+
maximumIterations := anInteger
152+
]
153+
154+
{ #category : #printing }
110155
PMJacobiTransformation >> printOn: aStream [
111-
"Append to the argument aStream, a sequence of characters that describes the receiver."
112-
| first |
113-
first := true.
114-
lowerRows do:
115-
[ :each |
116-
first ifTrue: [ first := false]
117-
ifFalse:[ aStream cr].
118-
each printOn: aStream.
119-
].
156+
"Append to the argument aStream, a sequence of characters that describes the receiver."
157+
158+
lowerRows do: [ :each | each printOn: aStream ] separatedBy: [ aStream cr ]
120159
]
121160

122161
{ #category : #transformation }
123162
PMJacobiTransformation >> sortEigenValues [
124-
"Private - Use a bubble sort."
125-
| n bound m |
126-
n := lowerRows size.
127-
bound := n.
128-
[ bound = 0 ]
129-
whileFalse: [ m := 0.
130-
1 to: bound - 1 do:
131-
[ :j |
132-
( result at: j) abs > ( result at: j + 1) abs
133-
ifFalse:[ self exchangeAt: j.
134-
m := j.
135-
].
136-
].
137-
bound := m.
138-
].
139-
]
140-
141-
{ #category : #information }
163+
"Private - Use a bubble sort."
164+
165+
| numberOfRows bound |
166+
numberOfRows := lowerRows size.
167+
bound := numberOfRows.
168+
[ bound = 0 ] whileFalse: [
169+
| m |
170+
m := 0.
171+
1 to: bound - 1 do: [ :j |
172+
(result at: j) abs > (result at: j + 1) abs ifFalse: [
173+
self exchangeAt: j.
174+
m := j ] ].
175+
bound := m ]
176+
]
177+
178+
{ #category : #transformation }
142179
PMJacobiTransformation >> transform [
143-
^PMMatrix rows: transform
180+
181+
^ PMMatrix rows: transform
144182
]
145183

146184
{ #category : #transformation }
147185
PMJacobiTransformation >> transformAt: anInteger1 and: anInteger2 [
148-
"Private"
186+
"Private"
187+
149188
| d t s c tau apq app aqq arp arq |
150-
apq := ( lowerRows at: anInteger2) at: anInteger1.
151-
apq = 0
152-
ifTrue: [ ^nil].
153-
app := ( lowerRows at: anInteger1) at: anInteger1.
154-
aqq := ( lowerRows at: anInteger2) at: anInteger2.
189+
apq := (lowerRows at: anInteger2) at: anInteger1.
190+
apq = 0 ifTrue: [ ^ nil ].
191+
app := (lowerRows at: anInteger1) at: anInteger1.
192+
aqq := (lowerRows at: anInteger2) at: anInteger2.
155193
d := aqq - app.
156194
arp := d * 0.5 / apq.
157-
t := arp > 0 ifTrue: [ 1 / ( ( arp squared + 1) sqrt + arp)]
158-
ifFalse:[ 1 / ( arp - ( arp squared + 1) sqrt)].
159-
c := 1 / ( t squared + 1) sqrt.
195+
t := arp > 0
196+
ifTrue: [ 1 / ((arp squared + 1) sqrt + arp) ]
197+
ifFalse: [ 1 / (arp - (arp squared + 1) sqrt) ].
198+
c := 1 / (t squared + 1) sqrt.
160199
s := t * c.
161-
tau := s / ( 1 + c).
162-
1 to: ( anInteger1 - 1)
163-
do: [ :r |
164-
arp := ( lowerRows at: anInteger1) at: r.
165-
arq := ( lowerRows at: anInteger2) at: r.
166-
( lowerRows at: anInteger1) at: r put: ( arp - ( s * (tau * arp + arq))).
167-
( lowerRows at: anInteger2) at: r put: ( arq + ( s * (arp - (tau * arq)))).
168-
].
169-
( anInteger1 + 1) to: ( anInteger2 - 1)
170-
do: [ :r |
171-
arp := ( lowerRows at: r) at: anInteger1.
172-
arq := ( lowerRows at: anInteger2) at: r.
173-
( lowerRows at: r) at: anInteger1 put: ( arp - ( s * (tau * arp + arq))).
174-
( lowerRows at: anInteger2) at: r put: ( arq + ( s * (arp - (tau * arq)))).
175-
].
176-
( anInteger2 + 1) to: lowerRows size
177-
do: [ :r |
178-
arp := ( lowerRows at: r) at: anInteger1.
179-
arq := ( lowerRows at: r) at: anInteger2.
180-
( lowerRows at: r) at: anInteger1 put: ( arp - ( s * (tau * arp + arq))).
181-
( lowerRows at: r) at: anInteger2 put: ( arq + ( s * (arp - (tau * arq)))).
182-
].
183-
1 to: lowerRows size
184-
do: [ :r |
185-
arp := ( transform at: r) at: anInteger1.
186-
arq := ( transform at: r) at: anInteger2.
187-
( transform at: r) at: anInteger1 put: ( arp - ( s * (tau * arp + arq))).
188-
( transform at: r) at: anInteger2 put: ( arq + ( s * (arp - (tau * arq)))).
189-
].
190-
( lowerRows at: anInteger1) at: anInteger1 put: ( app - (t * apq)).
191-
( lowerRows at: anInteger2) at: anInteger2 put: ( aqq + (t * apq)).
192-
( lowerRows at: anInteger2) at: anInteger1 put: 0.
200+
tau := s / (1 + c).
201+
1 to: anInteger1 - 1 do: [ :r |
202+
arp := (lowerRows at: anInteger1) at: r.
203+
arq := (lowerRows at: anInteger2) at: r.
204+
(lowerRows at: anInteger1) at: r put: arp - (s * (tau * arp + arq)).
205+
(lowerRows at: anInteger2) at: r put: arq + (s * (arp - (tau * arq))) ].
206+
anInteger1 + 1 to: anInteger2 - 1 do: [ :r |
207+
arp := (lowerRows at: r) at: anInteger1.
208+
arq := (lowerRows at: anInteger2) at: r.
209+
(lowerRows at: r) at: anInteger1 put: arp - (s * (tau * arp + arq)).
210+
(lowerRows at: anInteger2) at: r put: arq + (s * (arp - (tau * arq))) ].
211+
anInteger2 + 1 to: lowerRows size do: [ :r |
212+
arp := (lowerRows at: r) at: anInteger1.
213+
arq := (lowerRows at: r) at: anInteger2.
214+
(lowerRows at: r) at: anInteger1 put: arp - (s * (tau * arp + arq)).
215+
(lowerRows at: r) at: anInteger2 put: arq + (s * (arp - (tau * arq))) ].
216+
1 to: lowerRows size do: [ :r |
217+
arp := (transform at: r) at: anInteger1.
218+
arq := (transform at: r) at: anInteger2.
219+
(transform at: r) at: anInteger1 put: arp - (s * (tau * arp + arq)).
220+
(transform at: r) at: anInteger2 put: arq + (s * (arp - (tau * arq))) ].
221+
(lowerRows at: anInteger1) at: anInteger1 put: app - (t * apq).
222+
(lowerRows at: anInteger2) at: anInteger2 put: aqq + (t * apq).
223+
(lowerRows at: anInteger2) at: anInteger1 put: 0
193224
]

0 commit comments

Comments
 (0)