test_Eigen.praat 5.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188
  1. # test_Eigen.praat
  2. # djmw 20161116, 20180829
  3. appendInfoLine: "test_Eigen.praat"
  4. @testInterface
  5. eps = 1e-7
  6. for i to 10
  7. @testDiagonal: i
  8. endfor
  9. @test2by2
  10. @test3by3
  11. @testgeneralSquare
  12. procedure testInterface
  13. for .i to 5
  14. .numberOfColumns = randomInteger (3, 12)
  15. .tableofreal = Create TableOfReal: "t", 100, .numberOfColumns
  16. Formula: ~ randomGauss (0, 1)
  17. .pca = To PCA
  18. .eigen = Extract Eigen
  19. .numberOfEigenvalues = Get number of eigenvalues
  20. assert .numberOfEigenvalues == .numberOfColumns
  21. .dimension = Get eigenvector dimension
  22. assert .dimension == .numberOfColumns
  23. for .j to .numberOfEigenvalues
  24. .eigenvalue [.j] = Get eigenvalue: .j
  25. endfor
  26. for .j to .numberOfEigenvalues
  27. .sump = Get eigenvalue: .j
  28. for .k from .j to .numberOfEigenvalues
  29. .sum = Get sum of eigenvalues: .j, .k
  30. assert .sum >= .eigenvalue [.j]
  31. endfor
  32. endfor
  33. for .j to .numberOfEigenvalues
  34. for .k from .j to .dimension
  35. .val[.j,.k] = Get eigenvector element: .j, .k
  36. endfor
  37. endfor
  38. for .j to .numberOfEigenvalues
  39. for .k to .dimension
  40. .val[.k] = Get eigenvector element: .j, .k
  41. endfor
  42. Invert eigenvector: .j
  43. for .k to .dimension
  44. .valk = Get eigenvector element: .j, .k
  45. assert .valk == - .val[.k]
  46. endfor
  47. endfor
  48. removeObject: .tableofreal, .pca, .eigen
  49. endfor
  50. endproc
  51. procedure assertApproximatelyEqual: .val1, .val2, .eps, .comment$
  52. .diff = abs (.val1 -.val2)
  53. .tekst$ = .comment$ + " " + string$ (.val1) + ", " + string$ (.val2)
  54. if .val1 == 0
  55. assert .diff < .eps; '.tekst$'
  56. else
  57. .reldif = .diff / abs(.val1)
  58. assert .reldif < .eps ; '.tekst$'
  59. endif
  60. endproc
  61. procedure test2by2
  62. .dim = 2
  63. .mat = Create simple Matrix: "2x2s", .dim, .dim, "1"
  64. Set value: 1, 1, 2
  65. Set value: 2, 2, 2
  66. .eigenvalues# = {3, 1}
  67. # 20180829 clumsy because we cannot yet do mat##={{},{}}
  68. .eigenvec1# = {1/sqrt (2), 1/sqrt (2)}
  69. .eigenvec2# = {-1/sqrt (2), 1/sqrt (2)}
  70. .eigenvectors## = zero## (.dim, .dim)
  71. for .j to .dim
  72. .eigenvectors## [1, .j] = .eigenvec1# [.j]
  73. .eigenvectors## [2, .j] = .eigenvec2# [.j]
  74. endfor
  75. .eigen = To Eigen
  76. appendInfoLine: tab$, "2x2 symmetrical"
  77. @testeigen: .eigen, .dim, .eigenvalues#, .eigenvectors##
  78. removeObject: .mat, .eigen
  79. endproc
  80. procedure testeigen: .eigen, .dim, .eigenvalues#, .eigenvectors##
  81. selectObject: .eigen
  82. .numberOfEigenvalues = Get number of eigenvalues
  83. assert .numberOfEigenvalues == .dim
  84. for .i to .dim
  85. .eval = Get eigenvalue: .i
  86. .comment$ = string$ (.dim) + " eigenvalue" + string$ (.i)
  87. @assertApproximatelyEqual: .eval, .eigenvalues# [.i], eps, .comment$
  88. for .j to .dim
  89. .evecj = Get eigenvector element: .i, .j
  90. .comment$ = "eigenvector[" + string$ (.i) + "] [" +string$ (.j) + "]"
  91. .val = .eigenvectors## [.i,.j]
  92. @assertApproximatelyEqual: .evecj, .val, eps, .comment$
  93. endfor
  94. endfor
  95. endproc
  96. procedure test3by3
  97. .dim = 3
  98. .mat = Create simple Matrix: "3x3s", .dim, .dim, "0"
  99. Set value: 1, 1, 2
  100. Set value: 2, 2, 3
  101. Set value: 2, 3, 4
  102. Set value: 3, 2, 4
  103. Set value: 3, 3, 9
  104. .eigenvalues# = {11, 2 , 1}
  105. .eigenvec1# = {0, 1/sqrt (5), 2/sqrt (5)}
  106. .eigenvec2# = {1, 0, 0}
  107. .eigenvec3# = {0, -2/sqrt (5), 1/sqrt (5)}
  108. .eigenvectors## = zero## (.dim, .dim)
  109. for .j to .dim
  110. .eigenvectors## [1, .j] = .eigenvec1# [.j]
  111. .eigenvectors## [2, .j] = .eigenvec2# [.j]
  112. .eigenvectors## [3, .j] = .eigenvec3# [.j]
  113. endfor
  114. .eigen = To Eigen
  115. appendInfoLine: tab$, "3x3 symmetrical"
  116. @testeigen: .eigen, .dim, .eigenvalues#, .eigenvectors##
  117. removeObject: .mat, .eigen
  118. endproc
  119. procedure diagonalData: .dim
  120. .name$ = string$(.dim) + "x" + string$ (.dim)
  121. .mat = Create simple Matrix: .name$, .dim, .dim, "0"
  122. .eigenvalues# = zero# (.dim)
  123. .eigenvectors## = zero## (.dim, .dim)
  124. for .i to .dim
  125. .val = .dim - .i + 1
  126. Set value: .i, .i, .val
  127. .eigenvalues# [.i] = .val
  128. .eigenvectors## [.i,.i] = 1
  129. endfor
  130. endproc
  131. procedure testDiagonal: .dim
  132. @diagonalData: .dim
  133. .matname$ = selected$ ("Matrix")
  134. .eigen = To Eigen
  135. appendInfoLine: tab$, .matname$ +" diagonal"
  136. @testeigen: .eigen, .dim, diagonalData.eigenvalues#, diagonalData.eigenvectors##
  137. removeObject: .eigen, diagonalData.mat
  138. endproc
  139. procedure testgeneralSquare
  140. .dim = 3
  141. .name$ = "3x3square"
  142. .mat = Create simple Matrix: .name$, .dim, .dim, "0"
  143. Set value: 1, 2, 1
  144. Set value: 2, 3, 1
  145. Set value: 3, 1, 1
  146. .given_re# = {1, -1/2, -1/2}
  147. .given_im# = {0, 0.5*sqrt(3), -0.5*sqrt(3)}
  148. Eigen (complex)
  149. .eigenvectors = selected ("Matrix", 1)
  150. .eigenvalues = selected ("Matrix", 2)
  151. selectObject: .eigenvalues
  152. .nrow = Get number of rows
  153. assert .nrow = 3
  154. # lite version of equality: check for occurrence
  155. # the eigenvalues of a real square matrix are not "sorted". We only know that complex conjugate eigenvalues occur
  156. # have the one with positive imaginary part first.
  157. .eval_re# = {object [.eigenvalues, 1, 1], object [.eigenvalues, 2, 1],object [.eigenvalues, 3, 1]}
  158. .eval_im# = {object [.eigenvalues, 1, 2], object [.eigenvalues, 2, 2],object [.eigenvalues, 3, 2]}
  159. if .eval_re# [1] > 1-eps and .eval_re# [1] < 1+eps
  160. assert .eval_re# [2] / .given_re# [2] > 1-eps and .eval_re# [2] / .given_re# [2] < 1+eps
  161. assert .eval_re# [3] / .given_re# [3] > 1-eps and .eval_re# [3] / .given_re# [3] < 1+eps
  162. else
  163. assert .eval_re# [1] / .given_re# [2] > 1-eps and .eval_re# [1] / .given_re# [2] < 1+eps
  164. assert .eval_re# [2] / .given_re# [3] > 1-eps and .eval_re# [2] / .given_re# [3] < 1+eps
  165. assert .eval_re# [3] > 1-eps and .eval_re# [3] < 1+eps and .eval_im# [3] == 0
  166. endif
  167. selectObject: .eigenvectors
  168. .ncol = Get number of columns
  169. assert .ncol = 6
  170. removeObject: .mat, .eigenvectors, .eigenvalues
  171. endproc
  172. appendInfoLine: "test_Eigen.praat OK"