MatrixIssue386Test.php 34 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492
  1. <?php
  2. namespace MathPHP\Tests\LinearAlgebra\Regression;
  3. use MathPHP\LinearAlgebra\MatrixFactory;
  4. use MathPHP\LinearAlgebra\Vector;
  5. /**
  6. * Github issue #386 - https://github.com/markrogoyski/math-php/issues/386
  7. * Large 49x49 singular-ish matrix (depending on what the error tolerance is set to).
  8. * Very small determinant: 1.001788e-19
  9. * When doing Matrix->solve($b), it got a divide by zero error because it tried to solve using LU decomposition.
  10. */
  11. class MatrixIssue386Test extends \PHPUnit\Framework\TestCase
  12. {
  13. private const MATRIX = [
  14. [1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  15. [0,0.25,0.583333333,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  16. [0,0,0.166666667,0.666666667,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  17. [0,0,0,0.166666667,0.583333333,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  18. [0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  19. [0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.583333333,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  20. [0,0,0,0,0,0,0,0,0.0625,0.145833333,0.041666667,0,0,0,0,0.145833333,0.340277778,0.097222222,0,0,0,0,0.041666667,0.097222222,0.027777778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  21. [0,0,0,0,0,0,0,0,0,0.041666667,0.166666667,0.041666667,0,0,0,0,0.097222222,0.388888889,0.097222222,0,0,0,0,0.027777778,0.111111111,0.027777778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  22. [0,0,0,0,0,0,0,0,0,0,0.041666667,0.145833333,0.0625,0,0,0,0,0.097222222,0.340277778,0.145833333,0,0,0,0,0.027777778,0.097222222,0.041666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  23. [0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.583333333,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  24. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0.666666667,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  25. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.041666667,0.097222222,0.027777778,0,0,0,0,0.166666667,0.388888889,0.111111111,0,0,0,0,0.041666667,0.097222222,0.027777778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  26. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.027777778,0.111111111,0.027777778,0,0,0,0,0.111111111,0.444444444,0.111111111,0,0,0,0,0.027777778,0.111111111,0.027777778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  27. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.027777778,0.097222222,0.041666667,0,0,0,0,0.111111111,0.388888889,0.166666667,0,0,0,0,0.027777778,0.097222222,0.041666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  28. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0.666666667,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  29. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0.583333333,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0],
  30. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.041666667,0.097222222,0.027777778,0,0,0,0,0.145833333,0.340277778,0.097222222,0,0,0,0,0.0625,0.145833333,0.041666667,0,0,0,0,0,0,0,0,0,0],
  31. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.027777778,0.111111111,0.027777778,0,0,0,0,0.097222222,0.388888889,0.097222222,0,0,0,0,0.041666667,0.166666667,0.041666667,0,0,0,0,0,0,0,0,0],
  32. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.027777778,0.097222222,0.041666667,0,0,0,0,0.097222222,0.340277778,0.145833333,0,0,0,0,0.041666667,0.145833333,0.0625,0,0,0,0,0,0,0,0],
  33. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0.583333333,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0],
  34. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0],
  35. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0.583333333,0.166666667,0,0,0],
  36. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0.666666667,0.166666667,0,0],
  37. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0.583333333,0.25,0],
  38. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1],
  39. [0.125,0,0,0,0,0,0,0.59375,0,0,0,0,0,0,0.260416667,0,0,0,0,0,0,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  40. [0,0.03125,0.072916667,0.020833333,0,0,0,0,0.1484375,0.346354167,0.098958333,0,0,0,0,0.065104167,0.151909722,0.043402778,0,0,0,0,0.005208333,0.012152778,0.003472222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  41. [0,0,0.020833333,0.083333333,0.020833333,0,0,0,0,0.098958333,0.395833333,0.098958333,0,0,0,0,0.043402778,0.173611111,0.043402778,0,0,0,0,0.003472222,0.013888889,0.003472222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  42. [0,0,0,0.020833333,0.072916667,0.03125,0,0,0,0,0.098958333,0.346354167,0.1484375,0,0,0,0,0.043402778,0.151909722,0.065104167,0,0,0,0,0.003472222,0.012152778,0.005208333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  43. [0,0,0,0,0,0,0.125,0,0,0,0,0,0,0.59375,0,0,0,0,0,0,0.260416667,0,0,0,0,0,0,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  44. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0,0,0,0,0,0,0.260416667,0,0,0,0,0,0,0.59375,0,0,0,0,0,0,0.125,0,0,0,0,0,0],
  45. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.005208333,0.012152778,0.003472222,0,0,0,0,0.065104167,0.151909722,0.043402778,0,0,0,0,0.1484375,0.346354167,0.098958333,0,0,0,0,0.03125,0.072916667,0.020833333,0,0,0],
  46. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.003472222,0.013888889,0.003472222,0,0,0,0,0.043402778,0.173611111,0.043402778,0,0,0,0,0.098958333,0.395833333,0.098958333,0,0,0,0,0.020833333,0.083333333,0.020833333,0,0],
  47. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.003472222,0.012152778,0.005208333,0,0,0,0,0.043402778,0.151909722,0.065104167,0,0,0,0,0.098958333,0.346354167,0.1484375,0,0,0,0,0.020833333,0.072916667,0.03125,0],
  48. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0,0,0,0,0,0,0.260416667,0,0,0,0,0,0,0.59375,0,0,0,0,0,0,0.125],
  49. [0.125,0.59375,0.260416667,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  50. [0,0,0,0,0,0,0,0.03125,0.1484375,0.065104167,0.005208333,0,0,0,0.072916667,0.346354167,0.151909722,0.012152778,0,0,0,0.020833333,0.098958333,0.043402778,0.003472222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  51. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0.098958333,0.043402778,0.003472222,0,0,0,0.083333333,0.395833333,0.173611111,0.013888889,0,0,0,0.020833333,0.098958333,0.043402778,0.003472222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  52. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0.098958333,0.043402778,0.003472222,0,0,0,0.072916667,0.346354167,0.151909722,0.012152778,0,0,0,0.03125,0.1484375,0.065104167,0.005208333,0,0,0,0,0,0,0,0,0,0],
  53. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.125,0.59375,0.260416667,0.020833333,0,0,0],
  54. [0,0,0,0.020833333,0.260416667,0.59375,0.125,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  55. [0,0,0,0,0,0,0,0,0,0,0.005208333,0.065104167,0.1484375,0.03125,0,0,0,0.012152778,0.151909722,0.346354167,0.072916667,0,0,0,0.003472222,0.043402778,0.098958333,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  56. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.003472222,0.043402778,0.098958333,0.020833333,0,0,0,0.013888889,0.173611111,0.395833333,0.083333333,0,0,0,0.003472222,0.043402778,0.098958333,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  57. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.003472222,0.043402778,0.098958333,0.020833333,0,0,0,0.012152778,0.151909722,0.346354167,0.072916667,0,0,0,0.005208333,0.065104167,0.1484375,0.03125,0,0,0,0,0,0,0],
  58. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0.260416667,0.59375,0.125],
  59. [0.015625,0.07421875,0.032552083,0.002604167,0,0,0,0.07421875,0.352539063,0.154622396,0.012369792,0,0,0,0.032552083,0.154622396,0.06781684,0.005425347,0,0,0,0.002604167,0.012369792,0.005425347,0.000434028,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  60. [0,0,0,0.002604167,0.032552083,0.07421875,0.015625,0,0,0,0.012369792,0.154622396,0.352539063,0.07421875,0,0,0,0.005425347,0.06781684,0.154622396,0.032552083,0,0,0,0.000434028,0.005425347,0.012369792,0.002604167,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],
  61. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.002604167,0.012369792,0.005425347,0.000434028,0,0,0,0.032552083,0.154622396,0.06781684,0.005425347,0,0,0,0.07421875,0.352539063,0.154622396,0.012369792,0,0,0,0.015625,0.07421875,0.032552083,0.002604167,0,0,0],
  62. [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.000434028,0.005425347,0.012369792,0.002604167,0,0,0,0.005425347,0.06781684,0.154622396,0.032552083,0,0,0,0.012369792,0.154622396,0.352539063,0.07421875,0,0,0,0.002604167,0.032552083,0.07421875,0.015625],
  63. ];
  64. private const B = [
  65. 0,
  66. -0.100074402,
  67. -0.279235927,
  68. -0.506044043,
  69. -0.768516341,
  70. -0.955919161,
  71. -1.117129522,
  72. -1.352203887,
  73. -1.630181622,
  74. -1.939321332,
  75. -1.88849433,
  76. -2.106903819,
  77. -2.392150397,
  78. -2.714509406,
  79. -3.062767193,
  80. -2.774951587,
  81. -3.045226542,
  82. -3.37367948,
  83. -3.732547604,
  84. -4.11141885,
  85. -3.593643233,
  86. -3.909183413,
  87. -4.272821737,
  88. -4.65943525,
  89. -5.059667743,
  90. -0.479425294,
  91. -0.610468305,
  92. -0.818221399,
  93. -1.071388315,
  94. -1.358070774,
  95. -3.194062398,
  96. -3.487868198,
  97. -3.834975015,
  98. -4.208858884,
  99. -4.599605356,
  100. -0.035621956,
  101. -1.022823342,
  102. -1.985046905,
  103. -2.898794698,
  104. -3.741752572,
  105. -0.633344193,
  106. -1.781366011,
  107. -2.885886041,
  108. -3.919931388,
  109. -4.858250183,
  110. -0.530845567,
  111. -1.211057585,
  112. -3.330455609,
  113. -4.402550402,
  114. ];
  115. private const X = [
  116. 0.00000000,
  117. -0.01139036,
  118. -0.08909642,
  119. -0.27152341,
  120. -0.50022549,
  121. -0.67596776,
  122. -0.76851634,
  123. -0.31994364,
  124. -0.33853214,
  125. -0.43014886,
  126. -0.63201718,
  127. -0.87876797,
  128. -1.06563542,
  129. -1.16352499,
  130. -0.95982455,
  131. -0.99271596,
  132. -1.11151769,
  133. -1.35077506,
  134. -1.63175101,
  135. -1.83943015,
  136. -1.94724405,
  137. -1.89621359,
  138. -1.94975644,
  139. -2.10643553,
  140. -2.39593601,
  141. -2.72139070,
  142. -2.95521779,
  143. -3.07528632,
  144. -2.78628705,
  145. -2.85917348,
  146. -3.04990309,
  147. -3.38257539,
  148. -3.74455966,
  149. -3.99882324,
  150. -4.12821379,
  151. -3.33432751,
  152. -3.41892359,
  153. -3.62924964,
  154. -3.98528957,
  155. -4.36573594,
  156. -4.62948616,
  157. -4.76298568,
  158. -3.59364323,
  159. -3.68372552,
  160. -3.90295272,
  161. -4.26917767,
  162. -4.65726700,
  163. -4.92466621,
  164. -5.05966774,
  165. ];
  166. /**
  167. * @test Large 49x49 matrix with lots of zeros and small floating point values should be recognized as a singular matrix.
  168. * The singular determination uses the determinant to see if it is zero.
  169. * The determinant has to take in the floating point error tolerance, otherwise, smallNumber E-19 is not going to equal 0.
  170. *
  171. * R Shows this matrix as being singular:
  172. * > M = rbind(c(1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  173. * c(0,0.25,0.583333333,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  174. * c(0,0,0.166666667,0.666666667,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  175. * c(0,0,0,0.166666667,0.583333333,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  176. * c(0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  177. * c(0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.583333333,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  178. * c(0,0,0,0,0,0,0,0,0.0625,0.145833333,0.041666667,0,0,0,0,0.145833333,0.340277778,0.097222222,0,0,0,0,0.041666667,0.097222222,0.027777778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  179. * c(0,0,0,0,0,0,0,0,0,0.041666667,0.166666667,0.041666667,0,0,0,0,0.097222222,0.388888889,0.097222222,0,0,0,0,0.027777778,0.111111111,0.027777778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  180. * c(0,0,0,0,0,0,0,0,0,0,0.041666667,0.145833333,0.0625,0,0,0,0,0.097222222,0.340277778,0.145833333,0,0,0,0,0.027777778,0.097222222,0.041666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  181. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.583333333,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  182. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0.666666667,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  183. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.041666667,0.097222222,0.027777778,0,0,0,0,0.166666667,0.388888889,0.111111111,0,0,0,0,0.041666667,0.097222222,0.027777778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  184. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.027777778,0.111111111,0.027777778,0,0,0,0,0.111111111,0.444444444,0.111111111,0,0,0,0,0.027777778,0.111111111,0.027777778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  185. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.027777778,0.097222222,0.041666667,0,0,0,0,0.111111111,0.388888889,0.166666667,0,0,0,0,0.027777778,0.097222222,0.041666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  186. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0.666666667,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  187. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0.583333333,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0),
  188. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.041666667,0.097222222,0.027777778,0,0,0,0,0.145833333,0.340277778,0.097222222,0,0,0,0,0.0625,0.145833333,0.041666667,0,0,0,0,0,0,0,0,0,0),
  189. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.027777778,0.111111111,0.027777778,0,0,0,0,0.097222222,0.388888889,0.097222222,0,0,0,0,0.041666667,0.166666667,0.041666667,0,0,0,0,0,0,0,0,0),
  190. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.027777778,0.097222222,0.041666667,0,0,0,0,0.097222222,0.340277778,0.145833333,0,0,0,0,0.041666667,0.145833333,0.0625,0,0,0,0,0,0,0,0),
  191. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0.583333333,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0),
  192. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0),
  193. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0.583333333,0.166666667,0,0,0),
  194. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0.666666667,0.166666667,0,0),
  195. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0.583333333,0.25,0),
  196. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1),
  197. * c(0.125,0,0,0,0,0,0,0.59375,0,0,0,0,0,0,0.260416667,0,0,0,0,0,0,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  198. * c(0,0.03125,0.072916667,0.020833333,0,0,0,0,0.1484375,0.346354167,0.098958333,0,0,0,0,0.065104167,0.151909722,0.043402778,0,0,0,0,0.005208333,0.012152778,0.003472222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  199. * c(0,0,0.020833333,0.083333333,0.020833333,0,0,0,0,0.098958333,0.395833333,0.098958333,0,0,0,0,0.043402778,0.173611111,0.043402778,0,0,0,0,0.003472222,0.013888889,0.003472222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  200. * c(0,0,0,0.020833333,0.072916667,0.03125,0,0,0,0,0.098958333,0.346354167,0.1484375,0,0,0,0,0.043402778,0.151909722,0.065104167,0,0,0,0,0.003472222,0.012152778,0.005208333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  201. * c(0,0,0,0,0,0,0.125,0,0,0,0,0,0,0.59375,0,0,0,0,0,0,0.260416667,0,0,0,0,0,0,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  202. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0,0,0,0,0,0,0.260416667,0,0,0,0,0,0,0.59375,0,0,0,0,0,0,0.125,0,0,0,0,0,0),
  203. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.005208333,0.012152778,0.003472222,0,0,0,0,0.065104167,0.151909722,0.043402778,0,0,0,0,0.1484375,0.346354167,0.098958333,0,0,0,0,0.03125,0.072916667,0.020833333,0,0,0),
  204. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.003472222,0.013888889,0.003472222,0,0,0,0,0.043402778,0.173611111,0.043402778,0,0,0,0,0.098958333,0.395833333,0.098958333,0,0,0,0,0.020833333,0.083333333,0.020833333,0,0),
  205. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.003472222,0.012152778,0.005208333,0,0,0,0,0.043402778,0.151909722,0.065104167,0,0,0,0,0.098958333,0.346354167,0.1484375,0,0,0,0,0.020833333,0.072916667,0.03125,0),
  206. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0,0,0,0,0,0,0.260416667,0,0,0,0,0,0,0.59375,0,0,0,0,0,0,0.125),
  207. * c(0.125,0.59375,0.260416667,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  208. * c(0,0,0,0,0,0,0,0.03125,0.1484375,0.065104167,0.005208333,0,0,0,0.072916667,0.346354167,0.151909722,0.012152778,0,0,0,0.020833333,0.098958333,0.043402778,0.003472222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  209. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0.098958333,0.043402778,0.003472222,0,0,0,0.083333333,0.395833333,0.173611111,0.013888889,0,0,0,0.020833333,0.098958333,0.043402778,0.003472222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  210. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0.098958333,0.043402778,0.003472222,0,0,0,0.072916667,0.346354167,0.151909722,0.012152778,0,0,0,0.03125,0.1484375,0.065104167,0.005208333,0,0,0,0,0,0,0,0,0,0),
  211. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.125,0.59375,0.260416667,0.020833333,0,0,0),
  212. * c(0,0,0,0.020833333,0.260416667,0.59375,0.125,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  213. * c(0,0,0,0,0,0,0,0,0,0,0.005208333,0.065104167,0.1484375,0.03125,0,0,0,0.012152778,0.151909722,0.346354167,0.072916667,0,0,0,0.003472222,0.043402778,0.098958333,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  214. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.003472222,0.043402778,0.098958333,0.020833333,0,0,0,0.013888889,0.173611111,0.395833333,0.083333333,0,0,0,0.003472222,0.043402778,0.098958333,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  215. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.003472222,0.043402778,0.098958333,0.020833333,0,0,0,0.012152778,0.151909722,0.346354167,0.072916667,0,0,0,0.005208333,0.065104167,0.1484375,0.03125,0,0,0,0,0,0,0),
  216. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0.260416667,0.59375,0.125),
  217. * c(0.015625,0.07421875,0.032552083,0.002604167,0,0,0,0.07421875,0.352539063,0.154622396,0.012369792,0,0,0,0.032552083,0.154622396,0.06781684,0.005425347,0,0,0,0.002604167,0.012369792,0.005425347,0.000434028,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  218. * c(0,0,0,0.002604167,0.032552083,0.07421875,0.015625,0,0,0,0.012369792,0.154622396,0.352539063,0.07421875,0,0,0,0.005425347,0.06781684,0.154622396,0.032552083,0,0,0,0.000434028,0.005425347,0.012369792,0.002604167,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  219. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.002604167,0.012369792,0.005425347,0.000434028,0,0,0,0.032552083,0.154622396,0.06781684,0.005425347,0,0,0,0.07421875,0.352539063,0.154622396,0.012369792,0,0,0,0.015625,0.07421875,0.032552083,0.002604167,0,0,0),
  220. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.000434028,0.005425347,0.012369792,0.002604167,0,0,0,0.005425347,0.06781684,0.154622396,0.032552083,0,0,0,0.012369792,0.154622396,0.352539063,0.07421875,0,0,0,0.002604167,0.032552083,0.07421875,0.015625))
  221. * > is.singular.matrix(M)
  222. * [1] TRUE
  223. * @throws \Exception
  224. */
  225. public function testLargeSingularMatrixWithLotsOfFloatingPointValues()
  226. {
  227. // Given
  228. $A = MatrixFactory::create(self::MATRIX);
  229. // When
  230. $isSingular = $A->isSingular();
  231. $isNonsingular = $A->isNonsingular();
  232. // Then
  233. $this->assertTrue($isSingular);
  234. $this->assertFalse($isNonsingular);
  235. }
  236. /**
  237. * @test If the error tolerance is set to an extremely low value, then the matrix will not identify as being singular.
  238. *
  239. * R behavior:
  240. * > det(M)
  241. * [1] 1.001788e-19
  242. * > is.singular.matrix(M, tol=1e-19)
  243. * [1] FALSE
  244. * @throws \Exception
  245. */
  246. public function testLargeSingularMatrixWithLotsOfFloatingPointValuesUsingErrorTolerance()
  247. {
  248. // Given
  249. $A = MatrixFactory::create(self::MATRIX);
  250. // And
  251. $ε = 1e-19;
  252. $A->setError($ε);
  253. // When
  254. $isSingular = $A->isSingular();
  255. $isNonsingular = $A->isNonsingular();
  256. // Then
  257. $this->assertFalse($isSingular);
  258. $this->assertTrue($isNonsingular);
  259. }
  260. /**
  261. * @test Determinant of large singular matrix with lots of zeros and floating point values.
  262. *
  263. * R behavior:
  264. * > det(M)
  265. * [1] 1.001788e-19
  266. * @throws \Exception
  267. */
  268. public function testDeterminantOfSingularMatrix()
  269. {
  270. // Given
  271. $A = MatrixFactory::create(self::MATRIX);
  272. // When
  273. $det = $A->det();
  274. // Then
  275. $this->assertEqualsWithDelta(1.001788e-19, $det, 1e-25);
  276. }
  277. /**
  278. * @test Augmenting A with b and then computing the RREF solves Ax = b
  279. * The right most column of augmented Ab is x
  280. * @throws \Exception
  281. */
  282. public function testSolveRref()
  283. {
  284. // Given
  285. $A = MatrixFactory::create(self::MATRIX);
  286. $b = new Vector(self::B);
  287. // When
  288. $Ab = $A->augment($b->asColumnMatrix());
  289. $rref = $Ab->rref();
  290. $x = new Vector(\array_column($rref->getMatrix(), $rref->getN() - 1));
  291. // Then
  292. $this->assertEqualsWithDelta(self::X, $x->getVector(), 0.00000001);
  293. }
  294. /**
  295. * @test Solve (bug that was reported)
  296. * Original implementation of Matrix->solve($b) would use LU decomposition and fail by division by zero for this matrix.
  297. * @throws \Exception
  298. */
  299. public function testSolve()
  300. {
  301. // Given
  302. $A = MatrixFactory::create(self::MATRIX);
  303. $b = new Vector(self::B);
  304. // When
  305. $x = $A->solve($b);
  306. // Then
  307. $this->assertEqualsWithDelta(self::X, $x->getVector(), 0.00000001);
  308. }
  309. /*
  310. * R code to create the matrix (for verification purposes)
  311. * > A = rbind(c(1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  312. * c(0,0.25,0.583333333,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  313. * c(0,0,0.166666667,0.666666667,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  314. * c(0,0,0,0.166666667,0.583333333,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  315. * c(0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  316. * c(0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.583333333,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  317. * c(0,0,0,0,0,0,0,0,0.0625,0.145833333,0.041666667,0,0,0,0,0.145833333,0.340277778,0.097222222,0,0,0,0,0.041666667,0.097222222,0.027777778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  318. * c(0,0,0,0,0,0,0,0,0,0.041666667,0.166666667,0.041666667,0,0,0,0,0.097222222,0.388888889,0.097222222,0,0,0,0,0.027777778,0.111111111,0.027777778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  319. * c(0,0,0,0,0,0,0,0,0,0,0.041666667,0.145833333,0.0625,0,0,0,0,0.097222222,0.340277778,0.145833333,0,0,0,0,0.027777778,0.097222222,0.041666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  320. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0.583333333,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  321. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0.666666667,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  322. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.041666667,0.097222222,0.027777778,0,0,0,0,0.166666667,0.388888889,0.111111111,0,0,0,0,0.041666667,0.097222222,0.027777778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  323. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.027777778,0.111111111,0.027777778,0,0,0,0,0.111111111,0.444444444,0.111111111,0,0,0,0,0.027777778,0.111111111,0.027777778,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  324. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.027777778,0.097222222,0.041666667,0,0,0,0,0.111111111,0.388888889,0.166666667,0,0,0,0,0.027777778,0.097222222,0.041666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  325. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0.666666667,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  326. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0.583333333,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0,0,0,0,0,0,0),
  327. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.041666667,0.097222222,0.027777778,0,0,0,0,0.145833333,0.340277778,0.097222222,0,0,0,0,0.0625,0.145833333,0.041666667,0,0,0,0,0,0,0,0,0,0),
  328. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.027777778,0.111111111,0.027777778,0,0,0,0,0.097222222,0.388888889,0.097222222,0,0,0,0,0.041666667,0.166666667,0.041666667,0,0,0,0,0,0,0,0,0),
  329. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.027777778,0.097222222,0.041666667,0,0,0,0,0.097222222,0.340277778,0.145833333,0,0,0,0,0.041666667,0.145833333,0.0625,0,0,0,0,0,0,0,0),
  330. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0,0,0,0,0,0,0.583333333,0,0,0,0,0,0,0.25,0,0,0,0,0,0,0),
  331. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0),
  332. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.25,0.583333333,0.166666667,0,0,0),
  333. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0.666666667,0.166666667,0,0),
  334. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.166666667,0.583333333,0.25,0),
  335. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1),
  336. * c(0.125,0,0,0,0,0,0,0.59375,0,0,0,0,0,0,0.260416667,0,0,0,0,0,0,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  337. * c(0,0.03125,0.072916667,0.020833333,0,0,0,0,0.1484375,0.346354167,0.098958333,0,0,0,0,0.065104167,0.151909722,0.043402778,0,0,0,0,0.005208333,0.012152778,0.003472222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  338. * c(0,0,0.020833333,0.083333333,0.020833333,0,0,0,0,0.098958333,0.395833333,0.098958333,0,0,0,0,0.043402778,0.173611111,0.043402778,0,0,0,0,0.003472222,0.013888889,0.003472222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  339. * c(0,0,0,0.020833333,0.072916667,0.03125,0,0,0,0,0.098958333,0.346354167,0.1484375,0,0,0,0,0.043402778,0.151909722,0.065104167,0,0,0,0,0.003472222,0.012152778,0.005208333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  340. * c(0,0,0,0,0,0,0.125,0,0,0,0,0,0,0.59375,0,0,0,0,0,0,0.260416667,0,0,0,0,0,0,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  341. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0,0,0,0,0,0,0.260416667,0,0,0,0,0,0,0.59375,0,0,0,0,0,0,0.125,0,0,0,0,0,0),
  342. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.005208333,0.012152778,0.003472222,0,0,0,0,0.065104167,0.151909722,0.043402778,0,0,0,0,0.1484375,0.346354167,0.098958333,0,0,0,0,0.03125,0.072916667,0.020833333,0,0,0),
  343. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.003472222,0.013888889,0.003472222,0,0,0,0,0.043402778,0.173611111,0.043402778,0,0,0,0,0.098958333,0.395833333,0.098958333,0,0,0,0,0.020833333,0.083333333,0.020833333,0,0),
  344. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.003472222,0.012152778,0.005208333,0,0,0,0,0.043402778,0.151909722,0.065104167,0,0,0,0,0.098958333,0.346354167,0.1484375,0,0,0,0,0.020833333,0.072916667,0.03125,0),
  345. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0,0,0,0,0,0,0.260416667,0,0,0,0,0,0,0.59375,0,0,0,0,0,0,0.125),
  346. * c(0.125,0.59375,0.260416667,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  347. * c(0,0,0,0,0,0,0,0.03125,0.1484375,0.065104167,0.005208333,0,0,0,0.072916667,0.346354167,0.151909722,0.012152778,0,0,0,0.020833333,0.098958333,0.043402778,0.003472222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  348. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0.098958333,0.043402778,0.003472222,0,0,0,0.083333333,0.395833333,0.173611111,0.013888889,0,0,0,0.020833333,0.098958333,0.043402778,0.003472222,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  349. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0.098958333,0.043402778,0.003472222,0,0,0,0.072916667,0.346354167,0.151909722,0.012152778,0,0,0,0.03125,0.1484375,0.065104167,0.005208333,0,0,0,0,0,0,0,0,0,0),
  350. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.125,0.59375,0.260416667,0.020833333,0,0,0),
  351. * c(0,0,0,0.020833333,0.260416667,0.59375,0.125,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  352. * c(0,0,0,0,0,0,0,0,0,0,0.005208333,0.065104167,0.1484375,0.03125,0,0,0,0.012152778,0.151909722,0.346354167,0.072916667,0,0,0,0.003472222,0.043402778,0.098958333,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  353. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.003472222,0.043402778,0.098958333,0.020833333,0,0,0,0.013888889,0.173611111,0.395833333,0.083333333,0,0,0,0.003472222,0.043402778,0.098958333,0.020833333,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  354. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.003472222,0.043402778,0.098958333,0.020833333,0,0,0,0.012152778,0.151909722,0.346354167,0.072916667,0,0,0,0.005208333,0.065104167,0.1484375,0.03125,0,0,0,0,0,0,0),
  355. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.020833333,0.260416667,0.59375,0.125),
  356. * c(0.015625,0.07421875,0.032552083,0.002604167,0,0,0,0.07421875,0.352539063,0.154622396,0.012369792,0,0,0,0.032552083,0.154622396,0.06781684,0.005425347,0,0,0,0.002604167,0.012369792,0.005425347,0.000434028,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  357. * c(0,0,0,0.002604167,0.032552083,0.07421875,0.015625,0,0,0,0.012369792,0.154622396,0.352539063,0.07421875,0,0,0,0.005425347,0.06781684,0.154622396,0.032552083,0,0,0,0.000434028,0.005425347,0.012369792,0.002604167,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0),
  358. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.002604167,0.012369792,0.005425347,0.000434028,0,0,0,0.032552083,0.154622396,0.06781684,0.005425347,0,0,0,0.07421875,0.352539063,0.154622396,0.012369792,0,0,0,0.015625,0.07421875,0.032552083,0.002604167,0,0,0),
  359. * c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.000434028,0.005425347,0.012369792,0.002604167,0,0,0,0.005425347,0.06781684,0.154622396,0.032552083,0,0,0,0.012369792,0.154622396,0.352539063,0.07421875,0,0,0,0.002604167,0.032552083,0.07421875,0.015625))
  360. */
  361. /*
  362. * R code to create b vector (for verification purposes)
  363. * > b = c(0,
  364. * -0.100074402,
  365. * -0.279235927,
  366. * -0.506044043,
  367. * -0.768516341,
  368. * -0.955919161,
  369. * -1.117129522,
  370. * -1.352203887,
  371. * -1.630181622,
  372. * -1.939321332,
  373. * -1.88849433,
  374. * -2.106903819,
  375. * -2.392150397,
  376. * -2.714509406,
  377. * -3.062767193,
  378. * -2.774951587,
  379. * -3.045226542,
  380. * -3.37367948,
  381. * -3.732547604,
  382. * -4.11141885,
  383. * -3.593643233,
  384. * -3.909183413,
  385. * -4.272821737,
  386. * -4.65943525,
  387. * -5.059667743,
  388. * -0.479425294,
  389. * -0.610468305,
  390. * -0.818221399,
  391. * -1.071388315,
  392. * -1.358070774,
  393. * -3.194062398,
  394. * -3.487868198,
  395. * -3.834975015,
  396. * -4.208858884,
  397. * -4.599605356,
  398. * -0.035621956,
  399. * -1.022823342,
  400. * -1.985046905,
  401. * -2.898794698,
  402. * -3.741752572,
  403. * -0.633344193,
  404. * -1.781366011,
  405. * -2.885886041,
  406. * -3.919931388,
  407. * -4.858250183,
  408. * -0.530845567,
  409. * -1.211057585,
  410. * -3.330455609,
  411. * -4.402550402
  412. * )
  413. */
  414. /*
  415. * R code to solve Ax = b
  416. * > solve(A, b)
  417. * [,1]
  418. * [1,] 0.00000000
  419. * [2,] -0.01139036
  420. * [3,] -0.08909642
  421. * [4,] -0.27152341
  422. * [5,] -0.50022549
  423. * [6,] -0.67596776
  424. * [7,] -0.76851634
  425. * [8,] -0.31994364
  426. * [9,] -0.33853214
  427. * [10,] -0.43014886
  428. * [11,] -0.63201718
  429. * [12,] -0.87876797
  430. * [13,] -1.06563542
  431. * [14,] -1.16352499
  432. * [15,] -0.95982455
  433. * [16,] -0.99271596
  434. * [17,] -1.11151769
  435. * [18,] -1.35077506
  436. * [19,] -1.63175101
  437. * [20,] -1.83943015
  438. * [21,] -1.94724405
  439. * [22,] -1.89621359
  440. * [23,] -1.94975644
  441. * [24,] -2.10643553
  442. * [25,] -2.39593601
  443. * [26,] -2.72139070
  444. * [27,] -2.95521779
  445. * [28,] -3.07528632
  446. * [29,] -2.78628705
  447. * [30,] -2.85917348
  448. * [31,] -3.04990309
  449. * [32,] -3.38257539
  450. * [33,] -3.74455966
  451. * [34,] -3.99882324
  452. * [35,] -4.12821379
  453. * [36,] -3.33432751
  454. * [37,] -3.41892359
  455. * [38,] -3.62924964
  456. * [39,] -3.98528957
  457. * [40,] -4.36573594
  458. * [41,] -4.62948616
  459. * [42,] -4.76298568
  460. * [43,] -3.59364323
  461. * [44,] -3.68372552
  462. * [45,] -3.90295272
  463. * [46,] -4.26917767
  464. * [47,] -4.65726700
  465. * [48,] -4.92466621
  466. * [49,] -5.05966774
  467. */
  468. }