test_decomp_update.py 67 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700
  1. import itertools
  2. import numpy as np
  3. from numpy.testing import assert_, assert_allclose, assert_equal
  4. from pytest import raises as assert_raises
  5. from scipy import linalg
  6. import scipy.linalg._decomp_update as _decomp_update
  7. from scipy.linalg._decomp_update import qr_delete, qr_update, qr_insert
  8. def assert_unitary(a, rtol=None, atol=None, assert_sqr=True):
  9. if rtol is None:
  10. rtol = 10.0 ** -(np.finfo(a.dtype).precision-2)
  11. if atol is None:
  12. atol = 10*np.finfo(a.dtype).eps
  13. if assert_sqr:
  14. assert_(a.shape[0] == a.shape[1], 'unitary matrices must be square')
  15. aTa = np.dot(a.T.conj(), a)
  16. assert_allclose(aTa, np.eye(a.shape[1]), rtol=rtol, atol=atol)
  17. def assert_upper_tri(a, rtol=None, atol=None):
  18. if rtol is None:
  19. rtol = 10.0 ** -(np.finfo(a.dtype).precision-2)
  20. if atol is None:
  21. atol = 2*np.finfo(a.dtype).eps
  22. mask = np.tri(a.shape[0], a.shape[1], -1, np.bool_)
  23. assert_allclose(a[mask], 0.0, rtol=rtol, atol=atol)
  24. def check_qr(q, r, a, rtol, atol, assert_sqr=True):
  25. assert_unitary(q, rtol, atol, assert_sqr)
  26. assert_upper_tri(r, rtol, atol)
  27. assert_allclose(q.dot(r), a, rtol=rtol, atol=atol)
  28. def make_strided(arrs):
  29. strides = [(3, 7), (2, 2), (3, 4), (4, 2), (5, 4), (2, 3), (2, 1), (4, 5)]
  30. kmax = len(strides)
  31. k = 0
  32. ret = []
  33. for a in arrs:
  34. if a.ndim == 1:
  35. s = strides[k % kmax]
  36. k += 1
  37. base = np.zeros(s[0]*a.shape[0]+s[1], a.dtype)
  38. view = base[s[1]::s[0]]
  39. view[...] = a
  40. elif a.ndim == 2:
  41. s = strides[k % kmax]
  42. t = strides[(k+1) % kmax]
  43. k += 2
  44. base = np.zeros((s[0]*a.shape[0]+s[1], t[0]*a.shape[1]+t[1]),
  45. a.dtype)
  46. view = base[s[1]::s[0], t[1]::t[0]]
  47. view[...] = a
  48. else:
  49. raise ValueError('make_strided only works for ndim = 1 or'
  50. ' 2 arrays')
  51. ret.append(view)
  52. return ret
  53. def negate_strides(arrs):
  54. ret = []
  55. for a in arrs:
  56. b = np.zeros_like(a)
  57. if b.ndim == 2:
  58. b = b[::-1, ::-1]
  59. elif b.ndim == 1:
  60. b = b[::-1]
  61. else:
  62. raise ValueError('negate_strides only works for ndim = 1 or'
  63. ' 2 arrays')
  64. b[...] = a
  65. ret.append(b)
  66. return ret
  67. def nonitemsize_strides(arrs):
  68. out = []
  69. for a in arrs:
  70. a_dtype = a.dtype
  71. b = np.zeros(a.shape, [('a', a_dtype), ('junk', 'S1')])
  72. c = b.getfield(a_dtype)
  73. c[...] = a
  74. out.append(c)
  75. return out
  76. def make_nonnative(arrs):
  77. return [a.astype(a.dtype.newbyteorder()) for a in arrs]
  78. class BaseQRdeltas:
  79. def setup_method(self):
  80. self.rtol = 10.0 ** -(np.finfo(self.dtype).precision-2)
  81. self.atol = 10 * np.finfo(self.dtype).eps
  82. def generate(self, type, mode='full'):
  83. np.random.seed(29382)
  84. shape = {'sqr': (8, 8), 'tall': (12, 7), 'fat': (7, 12),
  85. 'Mx1': (8, 1), '1xN': (1, 8), '1x1': (1, 1)}[type]
  86. a = np.random.random(shape)
  87. if np.iscomplexobj(self.dtype.type(1)):
  88. b = np.random.random(shape)
  89. a = a + 1j * b
  90. a = a.astype(self.dtype)
  91. q, r = linalg.qr(a, mode=mode)
  92. return a, q, r
  93. class BaseQRdelete(BaseQRdeltas):
  94. def test_sqr_1_row(self):
  95. a, q, r = self.generate('sqr')
  96. for row in range(r.shape[0]):
  97. q1, r1 = qr_delete(q, r, row, overwrite_qr=False)
  98. a1 = np.delete(a, row, 0)
  99. check_qr(q1, r1, a1, self.rtol, self.atol)
  100. def test_sqr_p_row(self):
  101. a, q, r = self.generate('sqr')
  102. for ndel in range(2, 6):
  103. for row in range(a.shape[0]-ndel):
  104. q1, r1 = qr_delete(q, r, row, ndel, overwrite_qr=False)
  105. a1 = np.delete(a, slice(row, row+ndel), 0)
  106. check_qr(q1, r1, a1, self.rtol, self.atol)
  107. def test_sqr_1_col(self):
  108. a, q, r = self.generate('sqr')
  109. for col in range(r.shape[1]):
  110. q1, r1 = qr_delete(q, r, col, which='col', overwrite_qr=False)
  111. a1 = np.delete(a, col, 1)
  112. check_qr(q1, r1, a1, self.rtol, self.atol)
  113. def test_sqr_p_col(self):
  114. a, q, r = self.generate('sqr')
  115. for ndel in range(2, 6):
  116. for col in range(r.shape[1]-ndel):
  117. q1, r1 = qr_delete(q, r, col, ndel, which='col',
  118. overwrite_qr=False)
  119. a1 = np.delete(a, slice(col, col+ndel), 1)
  120. check_qr(q1, r1, a1, self.rtol, self.atol)
  121. def test_tall_1_row(self):
  122. a, q, r = self.generate('tall')
  123. for row in range(r.shape[0]):
  124. q1, r1 = qr_delete(q, r, row, overwrite_qr=False)
  125. a1 = np.delete(a, row, 0)
  126. check_qr(q1, r1, a1, self.rtol, self.atol)
  127. def test_tall_p_row(self):
  128. a, q, r = self.generate('tall')
  129. for ndel in range(2, 6):
  130. for row in range(a.shape[0]-ndel):
  131. q1, r1 = qr_delete(q, r, row, ndel, overwrite_qr=False)
  132. a1 = np.delete(a, slice(row, row+ndel), 0)
  133. check_qr(q1, r1, a1, self.rtol, self.atol)
  134. def test_tall_1_col(self):
  135. a, q, r = self.generate('tall')
  136. for col in range(r.shape[1]):
  137. q1, r1 = qr_delete(q, r, col, which='col', overwrite_qr=False)
  138. a1 = np.delete(a, col, 1)
  139. check_qr(q1, r1, a1, self.rtol, self.atol)
  140. def test_tall_p_col(self):
  141. a, q, r = self.generate('tall')
  142. for ndel in range(2, 6):
  143. for col in range(r.shape[1]-ndel):
  144. q1, r1 = qr_delete(q, r, col, ndel, which='col',
  145. overwrite_qr=False)
  146. a1 = np.delete(a, slice(col, col+ndel), 1)
  147. check_qr(q1, r1, a1, self.rtol, self.atol)
  148. def test_fat_1_row(self):
  149. a, q, r = self.generate('fat')
  150. for row in range(r.shape[0]):
  151. q1, r1 = qr_delete(q, r, row, overwrite_qr=False)
  152. a1 = np.delete(a, row, 0)
  153. check_qr(q1, r1, a1, self.rtol, self.atol)
  154. def test_fat_p_row(self):
  155. a, q, r = self.generate('fat')
  156. for ndel in range(2, 6):
  157. for row in range(a.shape[0]-ndel):
  158. q1, r1 = qr_delete(q, r, row, ndel, overwrite_qr=False)
  159. a1 = np.delete(a, slice(row, row+ndel), 0)
  160. check_qr(q1, r1, a1, self.rtol, self.atol)
  161. def test_fat_1_col(self):
  162. a, q, r = self.generate('fat')
  163. for col in range(r.shape[1]):
  164. q1, r1 = qr_delete(q, r, col, which='col', overwrite_qr=False)
  165. a1 = np.delete(a, col, 1)
  166. check_qr(q1, r1, a1, self.rtol, self.atol)
  167. def test_fat_p_col(self):
  168. a, q, r = self.generate('fat')
  169. for ndel in range(2, 6):
  170. for col in range(r.shape[1]-ndel):
  171. q1, r1 = qr_delete(q, r, col, ndel, which='col',
  172. overwrite_qr=False)
  173. a1 = np.delete(a, slice(col, col+ndel), 1)
  174. check_qr(q1, r1, a1, self.rtol, self.atol)
  175. def test_economic_1_row(self):
  176. # this test always starts and ends with an economic decomp.
  177. a, q, r = self.generate('tall', 'economic')
  178. for row in range(r.shape[0]):
  179. q1, r1 = qr_delete(q, r, row, overwrite_qr=False)
  180. a1 = np.delete(a, row, 0)
  181. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  182. # for economic row deletes
  183. # eco - prow = eco
  184. # eco - prow = sqr
  185. # eco - prow = fat
  186. def base_economic_p_row_xxx(self, ndel):
  187. a, q, r = self.generate('tall', 'economic')
  188. for row in range(a.shape[0]-ndel):
  189. q1, r1 = qr_delete(q, r, row, ndel, overwrite_qr=False)
  190. a1 = np.delete(a, slice(row, row+ndel), 0)
  191. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  192. def test_economic_p_row_economic(self):
  193. # (12, 7) - (3, 7) = (9,7) --> stays economic
  194. self.base_economic_p_row_xxx(3)
  195. def test_economic_p_row_sqr(self):
  196. # (12, 7) - (5, 7) = (7, 7) --> becomes square
  197. self.base_economic_p_row_xxx(5)
  198. def test_economic_p_row_fat(self):
  199. # (12, 7) - (7,7) = (5, 7) --> becomes fat
  200. self.base_economic_p_row_xxx(7)
  201. def test_economic_1_col(self):
  202. a, q, r = self.generate('tall', 'economic')
  203. for col in range(r.shape[1]):
  204. q1, r1 = qr_delete(q, r, col, which='col', overwrite_qr=False)
  205. a1 = np.delete(a, col, 1)
  206. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  207. def test_economic_p_col(self):
  208. a, q, r = self.generate('tall', 'economic')
  209. for ndel in range(2, 6):
  210. for col in range(r.shape[1]-ndel):
  211. q1, r1 = qr_delete(q, r, col, ndel, which='col',
  212. overwrite_qr=False)
  213. a1 = np.delete(a, slice(col, col+ndel), 1)
  214. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  215. def test_Mx1_1_row(self):
  216. a, q, r = self.generate('Mx1')
  217. for row in range(r.shape[0]):
  218. q1, r1 = qr_delete(q, r, row, overwrite_qr=False)
  219. a1 = np.delete(a, row, 0)
  220. check_qr(q1, r1, a1, self.rtol, self.atol)
  221. def test_Mx1_p_row(self):
  222. a, q, r = self.generate('Mx1')
  223. for ndel in range(2, 6):
  224. for row in range(a.shape[0]-ndel):
  225. q1, r1 = qr_delete(q, r, row, ndel, overwrite_qr=False)
  226. a1 = np.delete(a, slice(row, row+ndel), 0)
  227. check_qr(q1, r1, a1, self.rtol, self.atol)
  228. def test_1xN_1_col(self):
  229. a, q, r = self.generate('1xN')
  230. for col in range(r.shape[1]):
  231. q1, r1 = qr_delete(q, r, col, which='col', overwrite_qr=False)
  232. a1 = np.delete(a, col, 1)
  233. check_qr(q1, r1, a1, self.rtol, self.atol)
  234. def test_1xN_p_col(self):
  235. a, q, r = self.generate('1xN')
  236. for ndel in range(2, 6):
  237. for col in range(r.shape[1]-ndel):
  238. q1, r1 = qr_delete(q, r, col, ndel, which='col',
  239. overwrite_qr=False)
  240. a1 = np.delete(a, slice(col, col+ndel), 1)
  241. check_qr(q1, r1, a1, self.rtol, self.atol)
  242. def test_Mx1_economic_1_row(self):
  243. a, q, r = self.generate('Mx1', 'economic')
  244. for row in range(r.shape[0]):
  245. q1, r1 = qr_delete(q, r, row, overwrite_qr=False)
  246. a1 = np.delete(a, row, 0)
  247. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  248. def test_Mx1_economic_p_row(self):
  249. a, q, r = self.generate('Mx1', 'economic')
  250. for ndel in range(2, 6):
  251. for row in range(a.shape[0]-ndel):
  252. q1, r1 = qr_delete(q, r, row, ndel, overwrite_qr=False)
  253. a1 = np.delete(a, slice(row, row+ndel), 0)
  254. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  255. def test_delete_last_1_row(self):
  256. # full and eco are the same for 1xN
  257. a, q, r = self.generate('1xN')
  258. q1, r1 = qr_delete(q, r, 0, 1, 'row')
  259. assert_equal(q1, np.ndarray(shape=(0, 0), dtype=q.dtype))
  260. assert_equal(r1, np.ndarray(shape=(0, r.shape[1]), dtype=r.dtype))
  261. def test_delete_last_p_row(self):
  262. a, q, r = self.generate('tall', 'full')
  263. q1, r1 = qr_delete(q, r, 0, a.shape[0], 'row')
  264. assert_equal(q1, np.ndarray(shape=(0, 0), dtype=q.dtype))
  265. assert_equal(r1, np.ndarray(shape=(0, r.shape[1]), dtype=r.dtype))
  266. a, q, r = self.generate('tall', 'economic')
  267. q1, r1 = qr_delete(q, r, 0, a.shape[0], 'row')
  268. assert_equal(q1, np.ndarray(shape=(0, 0), dtype=q.dtype))
  269. assert_equal(r1, np.ndarray(shape=(0, r.shape[1]), dtype=r.dtype))
  270. def test_delete_last_1_col(self):
  271. a, q, r = self.generate('Mx1', 'economic')
  272. q1, r1 = qr_delete(q, r, 0, 1, 'col')
  273. assert_equal(q1, np.ndarray(shape=(q.shape[0], 0), dtype=q.dtype))
  274. assert_equal(r1, np.ndarray(shape=(0, 0), dtype=r.dtype))
  275. a, q, r = self.generate('Mx1', 'full')
  276. q1, r1 = qr_delete(q, r, 0, 1, 'col')
  277. assert_unitary(q1)
  278. assert_(q1.dtype == q.dtype)
  279. assert_(q1.shape == q.shape)
  280. assert_equal(r1, np.ndarray(shape=(r.shape[0], 0), dtype=r.dtype))
  281. def test_delete_last_p_col(self):
  282. a, q, r = self.generate('tall', 'full')
  283. q1, r1 = qr_delete(q, r, 0, a.shape[1], 'col')
  284. assert_unitary(q1)
  285. assert_(q1.dtype == q.dtype)
  286. assert_(q1.shape == q.shape)
  287. assert_equal(r1, np.ndarray(shape=(r.shape[0], 0), dtype=r.dtype))
  288. a, q, r = self.generate('tall', 'economic')
  289. q1, r1 = qr_delete(q, r, 0, a.shape[1], 'col')
  290. assert_equal(q1, np.ndarray(shape=(q.shape[0], 0), dtype=q.dtype))
  291. assert_equal(r1, np.ndarray(shape=(0, 0), dtype=r.dtype))
  292. def test_delete_1x1_row_col(self):
  293. a, q, r = self.generate('1x1')
  294. q1, r1 = qr_delete(q, r, 0, 1, 'row')
  295. assert_equal(q1, np.ndarray(shape=(0, 0), dtype=q.dtype))
  296. assert_equal(r1, np.ndarray(shape=(0, r.shape[1]), dtype=r.dtype))
  297. a, q, r = self.generate('1x1')
  298. q1, r1 = qr_delete(q, r, 0, 1, 'col')
  299. assert_unitary(q1)
  300. assert_(q1.dtype == q.dtype)
  301. assert_(q1.shape == q.shape)
  302. assert_equal(r1, np.ndarray(shape=(r.shape[0], 0), dtype=r.dtype))
  303. # all full qr, row deletes and single column deletes should be able to
  304. # handle any non negative strides. (only row and column vector
  305. # operations are used.) p column delete require fortran ordered
  306. # Q and R and will make a copy as necessary. Economic qr row deletes
  307. # requre a contigous q.
  308. def base_non_simple_strides(self, adjust_strides, ks, p, which,
  309. overwriteable):
  310. if which == 'row':
  311. qind = (slice(p,None), slice(p,None))
  312. rind = (slice(p,None), slice(None))
  313. else:
  314. qind = (slice(None), slice(None))
  315. rind = (slice(None), slice(None,-p))
  316. for type, k in itertools.product(['sqr', 'tall', 'fat'], ks):
  317. a, q0, r0, = self.generate(type)
  318. qs, rs = adjust_strides((q0, r0))
  319. if p == 1:
  320. a1 = np.delete(a, k, 0 if which == 'row' else 1)
  321. else:
  322. s = slice(k,k+p)
  323. if k < 0:
  324. s = slice(k, k + p +
  325. (a.shape[0] if which == 'row' else a.shape[1]))
  326. a1 = np.delete(a, s, 0 if which == 'row' else 1)
  327. # for each variable, q, r we try with it strided and
  328. # overwrite=False. Then we try with overwrite=True, and make
  329. # sure that q and r are still overwritten.
  330. q = q0.copy('F')
  331. r = r0.copy('F')
  332. q1, r1 = qr_delete(qs, r, k, p, which, False)
  333. check_qr(q1, r1, a1, self.rtol, self.atol)
  334. q1o, r1o = qr_delete(qs, r, k, p, which, True)
  335. check_qr(q1o, r1o, a1, self.rtol, self.atol)
  336. if overwriteable:
  337. assert_allclose(q1o, qs[qind], rtol=self.rtol, atol=self.atol)
  338. assert_allclose(r1o, r[rind], rtol=self.rtol, atol=self.atol)
  339. q = q0.copy('F')
  340. r = r0.copy('F')
  341. q2, r2 = qr_delete(q, rs, k, p, which, False)
  342. check_qr(q2, r2, a1, self.rtol, self.atol)
  343. q2o, r2o = qr_delete(q, rs, k, p, which, True)
  344. check_qr(q2o, r2o, a1, self.rtol, self.atol)
  345. if overwriteable:
  346. assert_allclose(q2o, q[qind], rtol=self.rtol, atol=self.atol)
  347. assert_allclose(r2o, rs[rind], rtol=self.rtol, atol=self.atol)
  348. q = q0.copy('F')
  349. r = r0.copy('F')
  350. # since some of these were consumed above
  351. qs, rs = adjust_strides((q, r))
  352. q3, r3 = qr_delete(qs, rs, k, p, which, False)
  353. check_qr(q3, r3, a1, self.rtol, self.atol)
  354. q3o, r3o = qr_delete(qs, rs, k, p, which, True)
  355. check_qr(q3o, r3o, a1, self.rtol, self.atol)
  356. if overwriteable:
  357. assert_allclose(q2o, qs[qind], rtol=self.rtol, atol=self.atol)
  358. assert_allclose(r3o, rs[rind], rtol=self.rtol, atol=self.atol)
  359. def test_non_unit_strides_1_row(self):
  360. self.base_non_simple_strides(make_strided, [0], 1, 'row', True)
  361. def test_non_unit_strides_p_row(self):
  362. self.base_non_simple_strides(make_strided, [0], 3, 'row', True)
  363. def test_non_unit_strides_1_col(self):
  364. self.base_non_simple_strides(make_strided, [0], 1, 'col', True)
  365. def test_non_unit_strides_p_col(self):
  366. self.base_non_simple_strides(make_strided, [0], 3, 'col', False)
  367. def test_neg_strides_1_row(self):
  368. self.base_non_simple_strides(negate_strides, [0], 1, 'row', False)
  369. def test_neg_strides_p_row(self):
  370. self.base_non_simple_strides(negate_strides, [0], 3, 'row', False)
  371. def test_neg_strides_1_col(self):
  372. self.base_non_simple_strides(negate_strides, [0], 1, 'col', False)
  373. def test_neg_strides_p_col(self):
  374. self.base_non_simple_strides(negate_strides, [0], 3, 'col', False)
  375. def test_non_itemize_strides_1_row(self):
  376. self.base_non_simple_strides(nonitemsize_strides, [0], 1, 'row', False)
  377. def test_non_itemize_strides_p_row(self):
  378. self.base_non_simple_strides(nonitemsize_strides, [0], 3, 'row', False)
  379. def test_non_itemize_strides_1_col(self):
  380. self.base_non_simple_strides(nonitemsize_strides, [0], 1, 'col', False)
  381. def test_non_itemize_strides_p_col(self):
  382. self.base_non_simple_strides(nonitemsize_strides, [0], 3, 'col', False)
  383. def test_non_native_byte_order_1_row(self):
  384. self.base_non_simple_strides(make_nonnative, [0], 1, 'row', False)
  385. def test_non_native_byte_order_p_row(self):
  386. self.base_non_simple_strides(make_nonnative, [0], 3, 'row', False)
  387. def test_non_native_byte_order_1_col(self):
  388. self.base_non_simple_strides(make_nonnative, [0], 1, 'col', False)
  389. def test_non_native_byte_order_p_col(self):
  390. self.base_non_simple_strides(make_nonnative, [0], 3, 'col', False)
  391. def test_neg_k(self):
  392. a, q, r = self.generate('sqr')
  393. for k, p, w in itertools.product([-3, -7], [1, 3], ['row', 'col']):
  394. q1, r1 = qr_delete(q, r, k, p, w, overwrite_qr=False)
  395. if w == 'row':
  396. a1 = np.delete(a, slice(k+a.shape[0], k+p+a.shape[0]), 0)
  397. else:
  398. a1 = np.delete(a, slice(k+a.shape[0], k+p+a.shape[1]), 1)
  399. check_qr(q1, r1, a1, self.rtol, self.atol)
  400. def base_overwrite_qr(self, which, p, test_C, test_F, mode='full'):
  401. assert_sqr = True if mode == 'full' else False
  402. if which == 'row':
  403. qind = (slice(p,None), slice(p,None))
  404. rind = (slice(p,None), slice(None))
  405. else:
  406. qind = (slice(None), slice(None))
  407. rind = (slice(None), slice(None,-p))
  408. a, q0, r0 = self.generate('sqr', mode)
  409. if p == 1:
  410. a1 = np.delete(a, 3, 0 if which == 'row' else 1)
  411. else:
  412. a1 = np.delete(a, slice(3, 3+p), 0 if which == 'row' else 1)
  413. # don't overwrite
  414. q = q0.copy('F')
  415. r = r0.copy('F')
  416. q1, r1 = qr_delete(q, r, 3, p, which, False)
  417. check_qr(q1, r1, a1, self.rtol, self.atol, assert_sqr)
  418. check_qr(q, r, a, self.rtol, self.atol, assert_sqr)
  419. if test_F:
  420. q = q0.copy('F')
  421. r = r0.copy('F')
  422. q2, r2 = qr_delete(q, r, 3, p, which, True)
  423. check_qr(q2, r2, a1, self.rtol, self.atol, assert_sqr)
  424. # verify the overwriting
  425. assert_allclose(q2, q[qind], rtol=self.rtol, atol=self.atol)
  426. assert_allclose(r2, r[rind], rtol=self.rtol, atol=self.atol)
  427. if test_C:
  428. q = q0.copy('C')
  429. r = r0.copy('C')
  430. q3, r3 = qr_delete(q, r, 3, p, which, True)
  431. check_qr(q3, r3, a1, self.rtol, self.atol, assert_sqr)
  432. assert_allclose(q3, q[qind], rtol=self.rtol, atol=self.atol)
  433. assert_allclose(r3, r[rind], rtol=self.rtol, atol=self.atol)
  434. def test_overwrite_qr_1_row(self):
  435. # any positively strided q and r.
  436. self.base_overwrite_qr('row', 1, True, True)
  437. def test_overwrite_economic_qr_1_row(self):
  438. # Any contiguous q and positively strided r.
  439. self.base_overwrite_qr('row', 1, True, True, 'economic')
  440. def test_overwrite_qr_1_col(self):
  441. # any positively strided q and r.
  442. # full and eco share code paths
  443. self.base_overwrite_qr('col', 1, True, True)
  444. def test_overwrite_qr_p_row(self):
  445. # any positively strided q and r.
  446. self.base_overwrite_qr('row', 3, True, True)
  447. def test_overwrite_economic_qr_p_row(self):
  448. # any contiguous q and positively strided r
  449. self.base_overwrite_qr('row', 3, True, True, 'economic')
  450. def test_overwrite_qr_p_col(self):
  451. # only F orderd q and r can be overwritten for cols
  452. # full and eco share code paths
  453. self.base_overwrite_qr('col', 3, False, True)
  454. def test_bad_which(self):
  455. a, q, r = self.generate('sqr')
  456. assert_raises(ValueError, qr_delete, q, r, 0, which='foo')
  457. def test_bad_k(self):
  458. a, q, r = self.generate('tall')
  459. assert_raises(ValueError, qr_delete, q, r, q.shape[0], 1)
  460. assert_raises(ValueError, qr_delete, q, r, -q.shape[0]-1, 1)
  461. assert_raises(ValueError, qr_delete, q, r, r.shape[0], 1, 'col')
  462. assert_raises(ValueError, qr_delete, q, r, -r.shape[0]-1, 1, 'col')
  463. def test_bad_p(self):
  464. a, q, r = self.generate('tall')
  465. # p must be positive
  466. assert_raises(ValueError, qr_delete, q, r, 0, -1)
  467. assert_raises(ValueError, qr_delete, q, r, 0, -1, 'col')
  468. # and nonzero
  469. assert_raises(ValueError, qr_delete, q, r, 0, 0)
  470. assert_raises(ValueError, qr_delete, q, r, 0, 0, 'col')
  471. # must have at least k+p rows or cols, depending.
  472. assert_raises(ValueError, qr_delete, q, r, 3, q.shape[0]-2)
  473. assert_raises(ValueError, qr_delete, q, r, 3, r.shape[1]-2, 'col')
  474. def test_empty_q(self):
  475. a, q, r = self.generate('tall')
  476. # same code path for 'row' and 'col'
  477. assert_raises(ValueError, qr_delete, np.array([]), r, 0, 1)
  478. def test_empty_r(self):
  479. a, q, r = self.generate('tall')
  480. # same code path for 'row' and 'col'
  481. assert_raises(ValueError, qr_delete, q, np.array([]), 0, 1)
  482. def test_mismatched_q_and_r(self):
  483. a, q, r = self.generate('tall')
  484. r = r[1:]
  485. assert_raises(ValueError, qr_delete, q, r, 0, 1)
  486. def test_unsupported_dtypes(self):
  487. dts = ['int8', 'int16', 'int32', 'int64',
  488. 'uint8', 'uint16', 'uint32', 'uint64',
  489. 'float16', 'longdouble', 'longcomplex',
  490. 'bool']
  491. a, q0, r0 = self.generate('tall')
  492. for dtype in dts:
  493. q = q0.real.astype(dtype)
  494. with np.errstate(invalid="ignore"):
  495. r = r0.real.astype(dtype)
  496. assert_raises(ValueError, qr_delete, q, r0, 0, 1, 'row')
  497. assert_raises(ValueError, qr_delete, q, r0, 0, 2, 'row')
  498. assert_raises(ValueError, qr_delete, q, r0, 0, 1, 'col')
  499. assert_raises(ValueError, qr_delete, q, r0, 0, 2, 'col')
  500. assert_raises(ValueError, qr_delete, q0, r, 0, 1, 'row')
  501. assert_raises(ValueError, qr_delete, q0, r, 0, 2, 'row')
  502. assert_raises(ValueError, qr_delete, q0, r, 0, 1, 'col')
  503. assert_raises(ValueError, qr_delete, q0, r, 0, 2, 'col')
  504. def test_check_finite(self):
  505. a0, q0, r0 = self.generate('tall')
  506. q = q0.copy('F')
  507. q[1,1] = np.nan
  508. assert_raises(ValueError, qr_delete, q, r0, 0, 1, 'row')
  509. assert_raises(ValueError, qr_delete, q, r0, 0, 3, 'row')
  510. assert_raises(ValueError, qr_delete, q, r0, 0, 1, 'col')
  511. assert_raises(ValueError, qr_delete, q, r0, 0, 3, 'col')
  512. r = r0.copy('F')
  513. r[1,1] = np.nan
  514. assert_raises(ValueError, qr_delete, q0, r, 0, 1, 'row')
  515. assert_raises(ValueError, qr_delete, q0, r, 0, 3, 'row')
  516. assert_raises(ValueError, qr_delete, q0, r, 0, 1, 'col')
  517. assert_raises(ValueError, qr_delete, q0, r, 0, 3, 'col')
  518. def test_qr_scalar(self):
  519. a, q, r = self.generate('1x1')
  520. assert_raises(ValueError, qr_delete, q[0, 0], r, 0, 1, 'row')
  521. assert_raises(ValueError, qr_delete, q, r[0, 0], 0, 1, 'row')
  522. assert_raises(ValueError, qr_delete, q[0, 0], r, 0, 1, 'col')
  523. assert_raises(ValueError, qr_delete, q, r[0, 0], 0, 1, 'col')
  524. class TestQRdelete_f(BaseQRdelete):
  525. dtype = np.dtype('f')
  526. class TestQRdelete_F(BaseQRdelete):
  527. dtype = np.dtype('F')
  528. class TestQRdelete_d(BaseQRdelete):
  529. dtype = np.dtype('d')
  530. class TestQRdelete_D(BaseQRdelete):
  531. dtype = np.dtype('D')
  532. class BaseQRinsert(BaseQRdeltas):
  533. def generate(self, type, mode='full', which='row', p=1):
  534. a, q, r = super().generate(type, mode)
  535. assert_(p > 0)
  536. # super call set the seed...
  537. if which == 'row':
  538. if p == 1:
  539. u = np.random.random(a.shape[1])
  540. else:
  541. u = np.random.random((p, a.shape[1]))
  542. elif which == 'col':
  543. if p == 1:
  544. u = np.random.random(a.shape[0])
  545. else:
  546. u = np.random.random((a.shape[0], p))
  547. else:
  548. ValueError('which should be either "row" or "col"')
  549. if np.iscomplexobj(self.dtype.type(1)):
  550. b = np.random.random(u.shape)
  551. u = u + 1j * b
  552. u = u.astype(self.dtype)
  553. return a, q, r, u
  554. def test_sqr_1_row(self):
  555. a, q, r, u = self.generate('sqr', which='row')
  556. for row in range(r.shape[0] + 1):
  557. q1, r1 = qr_insert(q, r, u, row)
  558. a1 = np.insert(a, row, u, 0)
  559. check_qr(q1, r1, a1, self.rtol, self.atol)
  560. def test_sqr_p_row(self):
  561. # sqr + rows --> fat always
  562. a, q, r, u = self.generate('sqr', which='row', p=3)
  563. for row in range(r.shape[0] + 1):
  564. q1, r1 = qr_insert(q, r, u, row)
  565. a1 = np.insert(a, np.full(3, row, np.intp), u, 0)
  566. check_qr(q1, r1, a1, self.rtol, self.atol)
  567. def test_sqr_1_col(self):
  568. a, q, r, u = self.generate('sqr', which='col')
  569. for col in range(r.shape[1] + 1):
  570. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  571. a1 = np.insert(a, col, u, 1)
  572. check_qr(q1, r1, a1, self.rtol, self.atol)
  573. def test_sqr_p_col(self):
  574. # sqr + cols --> fat always
  575. a, q, r, u = self.generate('sqr', which='col', p=3)
  576. for col in range(r.shape[1] + 1):
  577. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  578. a1 = np.insert(a, np.full(3, col, np.intp), u, 1)
  579. check_qr(q1, r1, a1, self.rtol, self.atol)
  580. def test_tall_1_row(self):
  581. a, q, r, u = self.generate('tall', which='row')
  582. for row in range(r.shape[0] + 1):
  583. q1, r1 = qr_insert(q, r, u, row)
  584. a1 = np.insert(a, row, u, 0)
  585. check_qr(q1, r1, a1, self.rtol, self.atol)
  586. def test_tall_p_row(self):
  587. # tall + rows --> tall always
  588. a, q, r, u = self.generate('tall', which='row', p=3)
  589. for row in range(r.shape[0] + 1):
  590. q1, r1 = qr_insert(q, r, u, row)
  591. a1 = np.insert(a, np.full(3, row, np.intp), u, 0)
  592. check_qr(q1, r1, a1, self.rtol, self.atol)
  593. def test_tall_1_col(self):
  594. a, q, r, u = self.generate('tall', which='col')
  595. for col in range(r.shape[1] + 1):
  596. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  597. a1 = np.insert(a, col, u, 1)
  598. check_qr(q1, r1, a1, self.rtol, self.atol)
  599. # for column adds to tall matrices there are three cases to test
  600. # tall + pcol --> tall
  601. # tall + pcol --> sqr
  602. # tall + pcol --> fat
  603. def base_tall_p_col_xxx(self, p):
  604. a, q, r, u = self.generate('tall', which='col', p=p)
  605. for col in range(r.shape[1] + 1):
  606. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  607. a1 = np.insert(a, np.full(p, col, np.intp), u, 1)
  608. check_qr(q1, r1, a1, self.rtol, self.atol)
  609. def test_tall_p_col_tall(self):
  610. # 12x7 + 12x3 = 12x10 --> stays tall
  611. self.base_tall_p_col_xxx(3)
  612. def test_tall_p_col_sqr(self):
  613. # 12x7 + 12x5 = 12x12 --> becomes sqr
  614. self.base_tall_p_col_xxx(5)
  615. def test_tall_p_col_fat(self):
  616. # 12x7 + 12x7 = 12x14 --> becomes fat
  617. self.base_tall_p_col_xxx(7)
  618. def test_fat_1_row(self):
  619. a, q, r, u = self.generate('fat', which='row')
  620. for row in range(r.shape[0] + 1):
  621. q1, r1 = qr_insert(q, r, u, row)
  622. a1 = np.insert(a, row, u, 0)
  623. check_qr(q1, r1, a1, self.rtol, self.atol)
  624. # for row adds to fat matrices there are three cases to test
  625. # fat + prow --> fat
  626. # fat + prow --> sqr
  627. # fat + prow --> tall
  628. def base_fat_p_row_xxx(self, p):
  629. a, q, r, u = self.generate('fat', which='row', p=p)
  630. for row in range(r.shape[0] + 1):
  631. q1, r1 = qr_insert(q, r, u, row)
  632. a1 = np.insert(a, np.full(p, row, np.intp), u, 0)
  633. check_qr(q1, r1, a1, self.rtol, self.atol)
  634. def test_fat_p_row_fat(self):
  635. # 7x12 + 3x12 = 10x12 --> stays fat
  636. self.base_fat_p_row_xxx(3)
  637. def test_fat_p_row_sqr(self):
  638. # 7x12 + 5x12 = 12x12 --> becomes sqr
  639. self.base_fat_p_row_xxx(5)
  640. def test_fat_p_row_tall(self):
  641. # 7x12 + 7x12 = 14x12 --> becomes tall
  642. self.base_fat_p_row_xxx(7)
  643. def test_fat_1_col(self):
  644. a, q, r, u = self.generate('fat', which='col')
  645. for col in range(r.shape[1] + 1):
  646. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  647. a1 = np.insert(a, col, u, 1)
  648. check_qr(q1, r1, a1, self.rtol, self.atol)
  649. def test_fat_p_col(self):
  650. # fat + cols --> fat always
  651. a, q, r, u = self.generate('fat', which='col', p=3)
  652. for col in range(r.shape[1] + 1):
  653. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  654. a1 = np.insert(a, np.full(3, col, np.intp), u, 1)
  655. check_qr(q1, r1, a1, self.rtol, self.atol)
  656. def test_economic_1_row(self):
  657. a, q, r, u = self.generate('tall', 'economic', 'row')
  658. for row in range(r.shape[0] + 1):
  659. q1, r1 = qr_insert(q, r, u, row, overwrite_qru=False)
  660. a1 = np.insert(a, row, u, 0)
  661. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  662. def test_economic_p_row(self):
  663. # tall + rows --> tall always
  664. a, q, r, u = self.generate('tall', 'economic', 'row', 3)
  665. for row in range(r.shape[0] + 1):
  666. q1, r1 = qr_insert(q, r, u, row, overwrite_qru=False)
  667. a1 = np.insert(a, np.full(3, row, np.intp), u, 0)
  668. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  669. def test_economic_1_col(self):
  670. a, q, r, u = self.generate('tall', 'economic', which='col')
  671. for col in range(r.shape[1] + 1):
  672. q1, r1 = qr_insert(q, r, u.copy(), col, 'col', overwrite_qru=False)
  673. a1 = np.insert(a, col, u, 1)
  674. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  675. def test_economic_1_col_bad_update(self):
  676. # When the column to be added lies in the span of Q, the update is
  677. # not meaningful. This is detected, and a LinAlgError is issued.
  678. q = np.eye(5, 3, dtype=self.dtype)
  679. r = np.eye(3, dtype=self.dtype)
  680. u = np.array([1, 0, 0, 0, 0], self.dtype)
  681. assert_raises(linalg.LinAlgError, qr_insert, q, r, u, 0, 'col')
  682. # for column adds to economic matrices there are three cases to test
  683. # eco + pcol --> eco
  684. # eco + pcol --> sqr
  685. # eco + pcol --> fat
  686. def base_economic_p_col_xxx(self, p):
  687. a, q, r, u = self.generate('tall', 'economic', which='col', p=p)
  688. for col in range(r.shape[1] + 1):
  689. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  690. a1 = np.insert(a, np.full(p, col, np.intp), u, 1)
  691. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  692. def test_economic_p_col_eco(self):
  693. # 12x7 + 12x3 = 12x10 --> stays eco
  694. self.base_economic_p_col_xxx(3)
  695. def test_economic_p_col_sqr(self):
  696. # 12x7 + 12x5 = 12x12 --> becomes sqr
  697. self.base_economic_p_col_xxx(5)
  698. def test_economic_p_col_fat(self):
  699. # 12x7 + 12x7 = 12x14 --> becomes fat
  700. self.base_economic_p_col_xxx(7)
  701. def test_Mx1_1_row(self):
  702. a, q, r, u = self.generate('Mx1', which='row')
  703. for row in range(r.shape[0] + 1):
  704. q1, r1 = qr_insert(q, r, u, row)
  705. a1 = np.insert(a, row, u, 0)
  706. check_qr(q1, r1, a1, self.rtol, self.atol)
  707. def test_Mx1_p_row(self):
  708. a, q, r, u = self.generate('Mx1', which='row', p=3)
  709. for row in range(r.shape[0] + 1):
  710. q1, r1 = qr_insert(q, r, u, row)
  711. a1 = np.insert(a, np.full(3, row, np.intp), u, 0)
  712. check_qr(q1, r1, a1, self.rtol, self.atol)
  713. def test_Mx1_1_col(self):
  714. a, q, r, u = self.generate('Mx1', which='col')
  715. for col in range(r.shape[1] + 1):
  716. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  717. a1 = np.insert(a, col, u, 1)
  718. check_qr(q1, r1, a1, self.rtol, self.atol)
  719. def test_Mx1_p_col(self):
  720. a, q, r, u = self.generate('Mx1', which='col', p=3)
  721. for col in range(r.shape[1] + 1):
  722. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  723. a1 = np.insert(a, np.full(3, col, np.intp), u, 1)
  724. check_qr(q1, r1, a1, self.rtol, self.atol)
  725. def test_Mx1_economic_1_row(self):
  726. a, q, r, u = self.generate('Mx1', 'economic', 'row')
  727. for row in range(r.shape[0] + 1):
  728. q1, r1 = qr_insert(q, r, u, row)
  729. a1 = np.insert(a, row, u, 0)
  730. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  731. def test_Mx1_economic_p_row(self):
  732. a, q, r, u = self.generate('Mx1', 'economic', 'row', 3)
  733. for row in range(r.shape[0] + 1):
  734. q1, r1 = qr_insert(q, r, u, row)
  735. a1 = np.insert(a, np.full(3, row, np.intp), u, 0)
  736. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  737. def test_Mx1_economic_1_col(self):
  738. a, q, r, u = self.generate('Mx1', 'economic', 'col')
  739. for col in range(r.shape[1] + 1):
  740. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  741. a1 = np.insert(a, col, u, 1)
  742. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  743. def test_Mx1_economic_p_col(self):
  744. a, q, r, u = self.generate('Mx1', 'economic', 'col', 3)
  745. for col in range(r.shape[1] + 1):
  746. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  747. a1 = np.insert(a, np.full(3, col, np.intp), u, 1)
  748. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  749. def test_1xN_1_row(self):
  750. a, q, r, u = self.generate('1xN', which='row')
  751. for row in range(r.shape[0] + 1):
  752. q1, r1 = qr_insert(q, r, u, row)
  753. a1 = np.insert(a, row, u, 0)
  754. check_qr(q1, r1, a1, self.rtol, self.atol)
  755. def test_1xN_p_row(self):
  756. a, q, r, u = self.generate('1xN', which='row', p=3)
  757. for row in range(r.shape[0] + 1):
  758. q1, r1 = qr_insert(q, r, u, row)
  759. a1 = np.insert(a, np.full(3, row, np.intp), u, 0)
  760. check_qr(q1, r1, a1, self.rtol, self.atol)
  761. def test_1xN_1_col(self):
  762. a, q, r, u = self.generate('1xN', which='col')
  763. for col in range(r.shape[1] + 1):
  764. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  765. a1 = np.insert(a, col, u, 1)
  766. check_qr(q1, r1, a1, self.rtol, self.atol)
  767. def test_1xN_p_col(self):
  768. a, q, r, u = self.generate('1xN', which='col', p=3)
  769. for col in range(r.shape[1] + 1):
  770. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  771. a1 = np.insert(a, np.full(3, col, np.intp), u, 1)
  772. check_qr(q1, r1, a1, self.rtol, self.atol)
  773. def test_1x1_1_row(self):
  774. a, q, r, u = self.generate('1x1', which='row')
  775. for row in range(r.shape[0] + 1):
  776. q1, r1 = qr_insert(q, r, u, row)
  777. a1 = np.insert(a, row, u, 0)
  778. check_qr(q1, r1, a1, self.rtol, self.atol)
  779. def test_1x1_p_row(self):
  780. a, q, r, u = self.generate('1x1', which='row', p=3)
  781. for row in range(r.shape[0] + 1):
  782. q1, r1 = qr_insert(q, r, u, row)
  783. a1 = np.insert(a, np.full(3, row, np.intp), u, 0)
  784. check_qr(q1, r1, a1, self.rtol, self.atol)
  785. def test_1x1_1_col(self):
  786. a, q, r, u = self.generate('1x1', which='col')
  787. for col in range(r.shape[1] + 1):
  788. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  789. a1 = np.insert(a, col, u, 1)
  790. check_qr(q1, r1, a1, self.rtol, self.atol)
  791. def test_1x1_p_col(self):
  792. a, q, r, u = self.generate('1x1', which='col', p=3)
  793. for col in range(r.shape[1] + 1):
  794. q1, r1 = qr_insert(q, r, u, col, 'col', overwrite_qru=False)
  795. a1 = np.insert(a, np.full(3, col, np.intp), u, 1)
  796. check_qr(q1, r1, a1, self.rtol, self.atol)
  797. def test_1x1_1_scalar(self):
  798. a, q, r, u = self.generate('1x1', which='row')
  799. assert_raises(ValueError, qr_insert, q[0, 0], r, u, 0, 'row')
  800. assert_raises(ValueError, qr_insert, q, r[0, 0], u, 0, 'row')
  801. assert_raises(ValueError, qr_insert, q, r, u[0], 0, 'row')
  802. assert_raises(ValueError, qr_insert, q[0, 0], r, u, 0, 'col')
  803. assert_raises(ValueError, qr_insert, q, r[0, 0], u, 0, 'col')
  804. assert_raises(ValueError, qr_insert, q, r, u[0], 0, 'col')
  805. def base_non_simple_strides(self, adjust_strides, k, p, which):
  806. for type in ['sqr', 'tall', 'fat']:
  807. a, q0, r0, u0 = self.generate(type, which=which, p=p)
  808. qs, rs, us = adjust_strides((q0, r0, u0))
  809. if p == 1:
  810. ai = np.insert(a, k, u0, 0 if which == 'row' else 1)
  811. else:
  812. ai = np.insert(a, np.full(p, k, np.intp),
  813. u0 if which == 'row' else u0,
  814. 0 if which == 'row' else 1)
  815. # for each variable, q, r, u we try with it strided and
  816. # overwrite=False. Then we try with overwrite=True. Nothing
  817. # is checked to see if it can be overwritten, since only
  818. # F ordered Q can be overwritten when adding columns.
  819. q = q0.copy('F')
  820. r = r0.copy('F')
  821. u = u0.copy('F')
  822. q1, r1 = qr_insert(qs, r, u, k, which, overwrite_qru=False)
  823. check_qr(q1, r1, ai, self.rtol, self.atol)
  824. q1o, r1o = qr_insert(qs, r, u, k, which, overwrite_qru=True)
  825. check_qr(q1o, r1o, ai, self.rtol, self.atol)
  826. q = q0.copy('F')
  827. r = r0.copy('F')
  828. u = u0.copy('F')
  829. q2, r2 = qr_insert(q, rs, u, k, which, overwrite_qru=False)
  830. check_qr(q2, r2, ai, self.rtol, self.atol)
  831. q2o, r2o = qr_insert(q, rs, u, k, which, overwrite_qru=True)
  832. check_qr(q2o, r2o, ai, self.rtol, self.atol)
  833. q = q0.copy('F')
  834. r = r0.copy('F')
  835. u = u0.copy('F')
  836. q3, r3 = qr_insert(q, r, us, k, which, overwrite_qru=False)
  837. check_qr(q3, r3, ai, self.rtol, self.atol)
  838. q3o, r3o = qr_insert(q, r, us, k, which, overwrite_qru=True)
  839. check_qr(q3o, r3o, ai, self.rtol, self.atol)
  840. q = q0.copy('F')
  841. r = r0.copy('F')
  842. u = u0.copy('F')
  843. # since some of these were consumed above
  844. qs, rs, us = adjust_strides((q, r, u))
  845. q5, r5 = qr_insert(qs, rs, us, k, which, overwrite_qru=False)
  846. check_qr(q5, r5, ai, self.rtol, self.atol)
  847. q5o, r5o = qr_insert(qs, rs, us, k, which, overwrite_qru=True)
  848. check_qr(q5o, r5o, ai, self.rtol, self.atol)
  849. def test_non_unit_strides_1_row(self):
  850. self.base_non_simple_strides(make_strided, 0, 1, 'row')
  851. def test_non_unit_strides_p_row(self):
  852. self.base_non_simple_strides(make_strided, 0, 3, 'row')
  853. def test_non_unit_strides_1_col(self):
  854. self.base_non_simple_strides(make_strided, 0, 1, 'col')
  855. def test_non_unit_strides_p_col(self):
  856. self.base_non_simple_strides(make_strided, 0, 3, 'col')
  857. def test_neg_strides_1_row(self):
  858. self.base_non_simple_strides(negate_strides, 0, 1, 'row')
  859. def test_neg_strides_p_row(self):
  860. self.base_non_simple_strides(negate_strides, 0, 3, 'row')
  861. def test_neg_strides_1_col(self):
  862. self.base_non_simple_strides(negate_strides, 0, 1, 'col')
  863. def test_neg_strides_p_col(self):
  864. self.base_non_simple_strides(negate_strides, 0, 3, 'col')
  865. def test_non_itemsize_strides_1_row(self):
  866. self.base_non_simple_strides(nonitemsize_strides, 0, 1, 'row')
  867. def test_non_itemsize_strides_p_row(self):
  868. self.base_non_simple_strides(nonitemsize_strides, 0, 3, 'row')
  869. def test_non_itemsize_strides_1_col(self):
  870. self.base_non_simple_strides(nonitemsize_strides, 0, 1, 'col')
  871. def test_non_itemsize_strides_p_col(self):
  872. self.base_non_simple_strides(nonitemsize_strides, 0, 3, 'col')
  873. def test_non_native_byte_order_1_row(self):
  874. self.base_non_simple_strides(make_nonnative, 0, 1, 'row')
  875. def test_non_native_byte_order_p_row(self):
  876. self.base_non_simple_strides(make_nonnative, 0, 3, 'row')
  877. def test_non_native_byte_order_1_col(self):
  878. self.base_non_simple_strides(make_nonnative, 0, 1, 'col')
  879. def test_non_native_byte_order_p_col(self):
  880. self.base_non_simple_strides(make_nonnative, 0, 3, 'col')
  881. def test_overwrite_qu_rank_1(self):
  882. # when inserting rows, the size of both Q and R change, so only
  883. # column inserts can overwrite q. Only complex column inserts
  884. # with C ordered Q overwrite u. Any contiguous Q is overwritten
  885. # when inserting 1 column
  886. a, q0, r, u, = self.generate('sqr', which='col', p=1)
  887. q = q0.copy('C')
  888. u0 = u.copy()
  889. # don't overwrite
  890. q1, r1 = qr_insert(q, r, u, 0, 'col', overwrite_qru=False)
  891. a1 = np.insert(a, 0, u0, 1)
  892. check_qr(q1, r1, a1, self.rtol, self.atol)
  893. check_qr(q, r, a, self.rtol, self.atol)
  894. # try overwriting
  895. q2, r2 = qr_insert(q, r, u, 0, 'col', overwrite_qru=True)
  896. check_qr(q2, r2, a1, self.rtol, self.atol)
  897. # verify the overwriting
  898. assert_allclose(q2, q, rtol=self.rtol, atol=self.atol)
  899. assert_allclose(u, u0.conj(), self.rtol, self.atol)
  900. # now try with a fortran ordered Q
  901. qF = q0.copy('F')
  902. u1 = u0.copy()
  903. q3, r3 = qr_insert(qF, r, u1, 0, 'col', overwrite_qru=False)
  904. check_qr(q3, r3, a1, self.rtol, self.atol)
  905. check_qr(qF, r, a, self.rtol, self.atol)
  906. # try overwriting
  907. q4, r4 = qr_insert(qF, r, u1, 0, 'col', overwrite_qru=True)
  908. check_qr(q4, r4, a1, self.rtol, self.atol)
  909. assert_allclose(q4, qF, rtol=self.rtol, atol=self.atol)
  910. def test_overwrite_qu_rank_p(self):
  911. # when inserting rows, the size of both Q and R change, so only
  912. # column inserts can potentially overwrite Q. In practice, only
  913. # F ordered Q are overwritten with a rank p update.
  914. a, q0, r, u, = self.generate('sqr', which='col', p=3)
  915. q = q0.copy('F')
  916. a1 = np.insert(a, np.zeros(3, np.intp), u, 1)
  917. # don't overwrite
  918. q1, r1 = qr_insert(q, r, u, 0, 'col', overwrite_qru=False)
  919. check_qr(q1, r1, a1, self.rtol, self.atol)
  920. check_qr(q, r, a, self.rtol, self.atol)
  921. # try overwriting
  922. q2, r2 = qr_insert(q, r, u, 0, 'col', overwrite_qru=True)
  923. check_qr(q2, r2, a1, self.rtol, self.atol)
  924. assert_allclose(q2, q, rtol=self.rtol, atol=self.atol)
  925. def test_empty_inputs(self):
  926. a, q, r, u = self.generate('sqr', which='row')
  927. assert_raises(ValueError, qr_insert, np.array([]), r, u, 0, 'row')
  928. assert_raises(ValueError, qr_insert, q, np.array([]), u, 0, 'row')
  929. assert_raises(ValueError, qr_insert, q, r, np.array([]), 0, 'row')
  930. assert_raises(ValueError, qr_insert, np.array([]), r, u, 0, 'col')
  931. assert_raises(ValueError, qr_insert, q, np.array([]), u, 0, 'col')
  932. assert_raises(ValueError, qr_insert, q, r, np.array([]), 0, 'col')
  933. def test_mismatched_shapes(self):
  934. a, q, r, u = self.generate('tall', which='row')
  935. assert_raises(ValueError, qr_insert, q, r[1:], u, 0, 'row')
  936. assert_raises(ValueError, qr_insert, q[:-2], r, u, 0, 'row')
  937. assert_raises(ValueError, qr_insert, q, r, u[1:], 0, 'row')
  938. assert_raises(ValueError, qr_insert, q, r[1:], u, 0, 'col')
  939. assert_raises(ValueError, qr_insert, q[:-2], r, u, 0, 'col')
  940. assert_raises(ValueError, qr_insert, q, r, u[1:], 0, 'col')
  941. def test_unsupported_dtypes(self):
  942. dts = ['int8', 'int16', 'int32', 'int64',
  943. 'uint8', 'uint16', 'uint32', 'uint64',
  944. 'float16', 'longdouble', 'longcomplex',
  945. 'bool']
  946. a, q0, r0, u0 = self.generate('sqr', which='row')
  947. for dtype in dts:
  948. q = q0.real.astype(dtype)
  949. with np.errstate(invalid="ignore"):
  950. r = r0.real.astype(dtype)
  951. u = u0.real.astype(dtype)
  952. assert_raises(ValueError, qr_insert, q, r0, u0, 0, 'row')
  953. assert_raises(ValueError, qr_insert, q, r0, u0, 0, 'col')
  954. assert_raises(ValueError, qr_insert, q0, r, u0, 0, 'row')
  955. assert_raises(ValueError, qr_insert, q0, r, u0, 0, 'col')
  956. assert_raises(ValueError, qr_insert, q0, r0, u, 0, 'row')
  957. assert_raises(ValueError, qr_insert, q0, r0, u, 0, 'col')
  958. def test_check_finite(self):
  959. a0, q0, r0, u0 = self.generate('sqr', which='row', p=3)
  960. q = q0.copy('F')
  961. q[1,1] = np.nan
  962. assert_raises(ValueError, qr_insert, q, r0, u0[:,0], 0, 'row')
  963. assert_raises(ValueError, qr_insert, q, r0, u0, 0, 'row')
  964. assert_raises(ValueError, qr_insert, q, r0, u0[:,0], 0, 'col')
  965. assert_raises(ValueError, qr_insert, q, r0, u0, 0, 'col')
  966. r = r0.copy('F')
  967. r[1,1] = np.nan
  968. assert_raises(ValueError, qr_insert, q0, r, u0[:,0], 0, 'row')
  969. assert_raises(ValueError, qr_insert, q0, r, u0, 0, 'row')
  970. assert_raises(ValueError, qr_insert, q0, r, u0[:,0], 0, 'col')
  971. assert_raises(ValueError, qr_insert, q0, r, u0, 0, 'col')
  972. u = u0.copy('F')
  973. u[0,0] = np.nan
  974. assert_raises(ValueError, qr_insert, q0, r0, u[:,0], 0, 'row')
  975. assert_raises(ValueError, qr_insert, q0, r0, u, 0, 'row')
  976. assert_raises(ValueError, qr_insert, q0, r0, u[:,0], 0, 'col')
  977. assert_raises(ValueError, qr_insert, q0, r0, u, 0, 'col')
  978. class TestQRinsert_f(BaseQRinsert):
  979. dtype = np.dtype('f')
  980. class TestQRinsert_F(BaseQRinsert):
  981. dtype = np.dtype('F')
  982. class TestQRinsert_d(BaseQRinsert):
  983. dtype = np.dtype('d')
  984. class TestQRinsert_D(BaseQRinsert):
  985. dtype = np.dtype('D')
  986. class BaseQRupdate(BaseQRdeltas):
  987. def generate(self, type, mode='full', p=1):
  988. a, q, r = super().generate(type, mode)
  989. # super call set the seed...
  990. if p == 1:
  991. u = np.random.random(q.shape[0])
  992. v = np.random.random(r.shape[1])
  993. else:
  994. u = np.random.random((q.shape[0], p))
  995. v = np.random.random((r.shape[1], p))
  996. if np.iscomplexobj(self.dtype.type(1)):
  997. b = np.random.random(u.shape)
  998. u = u + 1j * b
  999. c = np.random.random(v.shape)
  1000. v = v + 1j * c
  1001. u = u.astype(self.dtype)
  1002. v = v.astype(self.dtype)
  1003. return a, q, r, u, v
  1004. def test_sqr_rank_1(self):
  1005. a, q, r, u, v = self.generate('sqr')
  1006. q1, r1 = qr_update(q, r, u, v, False)
  1007. a1 = a + np.outer(u, v.conj())
  1008. check_qr(q1, r1, a1, self.rtol, self.atol)
  1009. def test_sqr_rank_p(self):
  1010. # test ndim = 2, rank 1 updates here too
  1011. for p in [1, 2, 3, 5]:
  1012. a, q, r, u, v = self.generate('sqr', p=p)
  1013. if p == 1:
  1014. u = u.reshape(u.size, 1)
  1015. v = v.reshape(v.size, 1)
  1016. q1, r1 = qr_update(q, r, u, v, False)
  1017. a1 = a + np.dot(u, v.T.conj())
  1018. check_qr(q1, r1, a1, self.rtol, self.atol)
  1019. def test_tall_rank_1(self):
  1020. a, q, r, u, v = self.generate('tall')
  1021. q1, r1 = qr_update(q, r, u, v, False)
  1022. a1 = a + np.outer(u, v.conj())
  1023. check_qr(q1, r1, a1, self.rtol, self.atol)
  1024. def test_tall_rank_p(self):
  1025. for p in [1, 2, 3, 5]:
  1026. a, q, r, u, v = self.generate('tall', p=p)
  1027. if p == 1:
  1028. u = u.reshape(u.size, 1)
  1029. v = v.reshape(v.size, 1)
  1030. q1, r1 = qr_update(q, r, u, v, False)
  1031. a1 = a + np.dot(u, v.T.conj())
  1032. check_qr(q1, r1, a1, self.rtol, self.atol)
  1033. def test_fat_rank_1(self):
  1034. a, q, r, u, v = self.generate('fat')
  1035. q1, r1 = qr_update(q, r, u, v, False)
  1036. a1 = a + np.outer(u, v.conj())
  1037. check_qr(q1, r1, a1, self.rtol, self.atol)
  1038. def test_fat_rank_p(self):
  1039. for p in [1, 2, 3, 5]:
  1040. a, q, r, u, v = self.generate('fat', p=p)
  1041. if p == 1:
  1042. u = u.reshape(u.size, 1)
  1043. v = v.reshape(v.size, 1)
  1044. q1, r1 = qr_update(q, r, u, v, False)
  1045. a1 = a + np.dot(u, v.T.conj())
  1046. check_qr(q1, r1, a1, self.rtol, self.atol)
  1047. def test_economic_rank_1(self):
  1048. a, q, r, u, v = self.generate('tall', 'economic')
  1049. q1, r1 = qr_update(q, r, u, v, False)
  1050. a1 = a + np.outer(u, v.conj())
  1051. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  1052. def test_economic_rank_p(self):
  1053. for p in [1, 2, 3, 5]:
  1054. a, q, r, u, v = self.generate('tall', 'economic', p)
  1055. if p == 1:
  1056. u = u.reshape(u.size, 1)
  1057. v = v.reshape(v.size, 1)
  1058. q1, r1 = qr_update(q, r, u, v, False)
  1059. a1 = a + np.dot(u, v.T.conj())
  1060. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  1061. def test_Mx1_rank_1(self):
  1062. a, q, r, u, v = self.generate('Mx1')
  1063. q1, r1 = qr_update(q, r, u, v, False)
  1064. a1 = a + np.outer(u, v.conj())
  1065. check_qr(q1, r1, a1, self.rtol, self.atol)
  1066. def test_Mx1_rank_p(self):
  1067. # when M or N == 1, only a rank 1 update is allowed. This isn't
  1068. # fundamental limitation, but the code does not support it.
  1069. a, q, r, u, v = self.generate('Mx1', p=1)
  1070. u = u.reshape(u.size, 1)
  1071. v = v.reshape(v.size, 1)
  1072. q1, r1 = qr_update(q, r, u, v, False)
  1073. a1 = a + np.dot(u, v.T.conj())
  1074. check_qr(q1, r1, a1, self.rtol, self.atol)
  1075. def test_Mx1_economic_rank_1(self):
  1076. a, q, r, u, v = self.generate('Mx1', 'economic')
  1077. q1, r1 = qr_update(q, r, u, v, False)
  1078. a1 = a + np.outer(u, v.conj())
  1079. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  1080. def test_Mx1_economic_rank_p(self):
  1081. # when M or N == 1, only a rank 1 update is allowed. This isn't
  1082. # fundamental limitation, but the code does not support it.
  1083. a, q, r, u, v = self.generate('Mx1', 'economic', p=1)
  1084. u = u.reshape(u.size, 1)
  1085. v = v.reshape(v.size, 1)
  1086. q1, r1 = qr_update(q, r, u, v, False)
  1087. a1 = a + np.dot(u, v.T.conj())
  1088. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  1089. def test_1xN_rank_1(self):
  1090. a, q, r, u, v = self.generate('1xN')
  1091. q1, r1 = qr_update(q, r, u, v, False)
  1092. a1 = a + np.outer(u, v.conj())
  1093. check_qr(q1, r1, a1, self.rtol, self.atol)
  1094. def test_1xN_rank_p(self):
  1095. # when M or N == 1, only a rank 1 update is allowed. This isn't
  1096. # fundamental limitation, but the code does not support it.
  1097. a, q, r, u, v = self.generate('1xN', p=1)
  1098. u = u.reshape(u.size, 1)
  1099. v = v.reshape(v.size, 1)
  1100. q1, r1 = qr_update(q, r, u, v, False)
  1101. a1 = a + np.dot(u, v.T.conj())
  1102. check_qr(q1, r1, a1, self.rtol, self.atol)
  1103. def test_1x1_rank_1(self):
  1104. a, q, r, u, v = self.generate('1x1')
  1105. q1, r1 = qr_update(q, r, u, v, False)
  1106. a1 = a + np.outer(u, v.conj())
  1107. check_qr(q1, r1, a1, self.rtol, self.atol)
  1108. def test_1x1_rank_p(self):
  1109. # when M or N == 1, only a rank 1 update is allowed. This isn't
  1110. # fundamental limitation, but the code does not support it.
  1111. a, q, r, u, v = self.generate('1x1', p=1)
  1112. u = u.reshape(u.size, 1)
  1113. v = v.reshape(v.size, 1)
  1114. q1, r1 = qr_update(q, r, u, v, False)
  1115. a1 = a + np.dot(u, v.T.conj())
  1116. check_qr(q1, r1, a1, self.rtol, self.atol)
  1117. def test_1x1_rank_1_scalar(self):
  1118. a, q, r, u, v = self.generate('1x1')
  1119. assert_raises(ValueError, qr_update, q[0, 0], r, u, v)
  1120. assert_raises(ValueError, qr_update, q, r[0, 0], u, v)
  1121. assert_raises(ValueError, qr_update, q, r, u[0], v)
  1122. assert_raises(ValueError, qr_update, q, r, u, v[0])
  1123. def base_non_simple_strides(self, adjust_strides, mode, p, overwriteable):
  1124. assert_sqr = False if mode == 'economic' else True
  1125. for type in ['sqr', 'tall', 'fat']:
  1126. a, q0, r0, u0, v0 = self.generate(type, mode, p)
  1127. qs, rs, us, vs = adjust_strides((q0, r0, u0, v0))
  1128. if p == 1:
  1129. aup = a + np.outer(u0, v0.conj())
  1130. else:
  1131. aup = a + np.dot(u0, v0.T.conj())
  1132. # for each variable, q, r, u, v we try with it strided and
  1133. # overwrite=False. Then we try with overwrite=True, and make
  1134. # sure that if p == 1, r and v are still overwritten.
  1135. # a strided q and u must always be copied.
  1136. q = q0.copy('F')
  1137. r = r0.copy('F')
  1138. u = u0.copy('F')
  1139. v = v0.copy('C')
  1140. q1, r1 = qr_update(qs, r, u, v, False)
  1141. check_qr(q1, r1, aup, self.rtol, self.atol, assert_sqr)
  1142. q1o, r1o = qr_update(qs, r, u, v, True)
  1143. check_qr(q1o, r1o, aup, self.rtol, self.atol, assert_sqr)
  1144. if overwriteable:
  1145. assert_allclose(r1o, r, rtol=self.rtol, atol=self.atol)
  1146. assert_allclose(v, v0.conj(), rtol=self.rtol, atol=self.atol)
  1147. q = q0.copy('F')
  1148. r = r0.copy('F')
  1149. u = u0.copy('F')
  1150. v = v0.copy('C')
  1151. q2, r2 = qr_update(q, rs, u, v, False)
  1152. check_qr(q2, r2, aup, self.rtol, self.atol, assert_sqr)
  1153. q2o, r2o = qr_update(q, rs, u, v, True)
  1154. check_qr(q2o, r2o, aup, self.rtol, self.atol, assert_sqr)
  1155. if overwriteable:
  1156. assert_allclose(r2o, rs, rtol=self.rtol, atol=self.atol)
  1157. assert_allclose(v, v0.conj(), rtol=self.rtol, atol=self.atol)
  1158. q = q0.copy('F')
  1159. r = r0.copy('F')
  1160. u = u0.copy('F')
  1161. v = v0.copy('C')
  1162. q3, r3 = qr_update(q, r, us, v, False)
  1163. check_qr(q3, r3, aup, self.rtol, self.atol, assert_sqr)
  1164. q3o, r3o = qr_update(q, r, us, v, True)
  1165. check_qr(q3o, r3o, aup, self.rtol, self.atol, assert_sqr)
  1166. if overwriteable:
  1167. assert_allclose(r3o, r, rtol=self.rtol, atol=self.atol)
  1168. assert_allclose(v, v0.conj(), rtol=self.rtol, atol=self.atol)
  1169. q = q0.copy('F')
  1170. r = r0.copy('F')
  1171. u = u0.copy('F')
  1172. v = v0.copy('C')
  1173. q4, r4 = qr_update(q, r, u, vs, False)
  1174. check_qr(q4, r4, aup, self.rtol, self.atol, assert_sqr)
  1175. q4o, r4o = qr_update(q, r, u, vs, True)
  1176. check_qr(q4o, r4o, aup, self.rtol, self.atol, assert_sqr)
  1177. if overwriteable:
  1178. assert_allclose(r4o, r, rtol=self.rtol, atol=self.atol)
  1179. assert_allclose(vs, v0.conj(), rtol=self.rtol, atol=self.atol)
  1180. q = q0.copy('F')
  1181. r = r0.copy('F')
  1182. u = u0.copy('F')
  1183. v = v0.copy('C')
  1184. # since some of these were consumed above
  1185. qs, rs, us, vs = adjust_strides((q, r, u, v))
  1186. q5, r5 = qr_update(qs, rs, us, vs, False)
  1187. check_qr(q5, r5, aup, self.rtol, self.atol, assert_sqr)
  1188. q5o, r5o = qr_update(qs, rs, us, vs, True)
  1189. check_qr(q5o, r5o, aup, self.rtol, self.atol, assert_sqr)
  1190. if overwriteable:
  1191. assert_allclose(r5o, rs, rtol=self.rtol, atol=self.atol)
  1192. assert_allclose(vs, v0.conj(), rtol=self.rtol, atol=self.atol)
  1193. def test_non_unit_strides_rank_1(self):
  1194. self.base_non_simple_strides(make_strided, 'full', 1, True)
  1195. def test_non_unit_strides_economic_rank_1(self):
  1196. self.base_non_simple_strides(make_strided, 'economic', 1, True)
  1197. def test_non_unit_strides_rank_p(self):
  1198. self.base_non_simple_strides(make_strided, 'full', 3, False)
  1199. def test_non_unit_strides_economic_rank_p(self):
  1200. self.base_non_simple_strides(make_strided, 'economic', 3, False)
  1201. def test_neg_strides_rank_1(self):
  1202. self.base_non_simple_strides(negate_strides, 'full', 1, False)
  1203. def test_neg_strides_economic_rank_1(self):
  1204. self.base_non_simple_strides(negate_strides, 'economic', 1, False)
  1205. def test_neg_strides_rank_p(self):
  1206. self.base_non_simple_strides(negate_strides, 'full', 3, False)
  1207. def test_neg_strides_economic_rank_p(self):
  1208. self.base_non_simple_strides(negate_strides, 'economic', 3, False)
  1209. def test_non_itemsize_strides_rank_1(self):
  1210. self.base_non_simple_strides(nonitemsize_strides, 'full', 1, False)
  1211. def test_non_itemsize_strides_economic_rank_1(self):
  1212. self.base_non_simple_strides(nonitemsize_strides, 'economic', 1, False)
  1213. def test_non_itemsize_strides_rank_p(self):
  1214. self.base_non_simple_strides(nonitemsize_strides, 'full', 3, False)
  1215. def test_non_itemsize_strides_economic_rank_p(self):
  1216. self.base_non_simple_strides(nonitemsize_strides, 'economic', 3, False)
  1217. def test_non_native_byte_order_rank_1(self):
  1218. self.base_non_simple_strides(make_nonnative, 'full', 1, False)
  1219. def test_non_native_byte_order_economic_rank_1(self):
  1220. self.base_non_simple_strides(make_nonnative, 'economic', 1, False)
  1221. def test_non_native_byte_order_rank_p(self):
  1222. self.base_non_simple_strides(make_nonnative, 'full', 3, False)
  1223. def test_non_native_byte_order_economic_rank_p(self):
  1224. self.base_non_simple_strides(make_nonnative, 'economic', 3, False)
  1225. def test_overwrite_qruv_rank_1(self):
  1226. # Any positive strided q, r, u, and v can be overwritten for a rank 1
  1227. # update, only checking C and F contiguous.
  1228. a, q0, r0, u0, v0 = self.generate('sqr')
  1229. a1 = a + np.outer(u0, v0.conj())
  1230. q = q0.copy('F')
  1231. r = r0.copy('F')
  1232. u = u0.copy('F')
  1233. v = v0.copy('F')
  1234. # don't overwrite
  1235. q1, r1 = qr_update(q, r, u, v, False)
  1236. check_qr(q1, r1, a1, self.rtol, self.atol)
  1237. check_qr(q, r, a, self.rtol, self.atol)
  1238. q2, r2 = qr_update(q, r, u, v, True)
  1239. check_qr(q2, r2, a1, self.rtol, self.atol)
  1240. # verify the overwriting, no good way to check u and v.
  1241. assert_allclose(q2, q, rtol=self.rtol, atol=self.atol)
  1242. assert_allclose(r2, r, rtol=self.rtol, atol=self.atol)
  1243. q = q0.copy('C')
  1244. r = r0.copy('C')
  1245. u = u0.copy('C')
  1246. v = v0.copy('C')
  1247. q3, r3 = qr_update(q, r, u, v, True)
  1248. check_qr(q3, r3, a1, self.rtol, self.atol)
  1249. assert_allclose(q3, q, rtol=self.rtol, atol=self.atol)
  1250. assert_allclose(r3, r, rtol=self.rtol, atol=self.atol)
  1251. def test_overwrite_qruv_rank_1_economic(self):
  1252. # updating economic decompositions can overwrite any contigous r,
  1253. # and positively strided r and u. V is only ever read.
  1254. # only checking C and F contiguous.
  1255. a, q0, r0, u0, v0 = self.generate('tall', 'economic')
  1256. a1 = a + np.outer(u0, v0.conj())
  1257. q = q0.copy('F')
  1258. r = r0.copy('F')
  1259. u = u0.copy('F')
  1260. v = v0.copy('F')
  1261. # don't overwrite
  1262. q1, r1 = qr_update(q, r, u, v, False)
  1263. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  1264. check_qr(q, r, a, self.rtol, self.atol, False)
  1265. q2, r2 = qr_update(q, r, u, v, True)
  1266. check_qr(q2, r2, a1, self.rtol, self.atol, False)
  1267. # verify the overwriting, no good way to check u and v.
  1268. assert_allclose(q2, q, rtol=self.rtol, atol=self.atol)
  1269. assert_allclose(r2, r, rtol=self.rtol, atol=self.atol)
  1270. q = q0.copy('C')
  1271. r = r0.copy('C')
  1272. u = u0.copy('C')
  1273. v = v0.copy('C')
  1274. q3, r3 = qr_update(q, r, u, v, True)
  1275. check_qr(q3, r3, a1, self.rtol, self.atol, False)
  1276. assert_allclose(q3, q, rtol=self.rtol, atol=self.atol)
  1277. assert_allclose(r3, r, rtol=self.rtol, atol=self.atol)
  1278. def test_overwrite_qruv_rank_p(self):
  1279. # for rank p updates, q r must be F contiguous, v must be C (v.T --> F)
  1280. # and u can be C or F, but is only overwritten if Q is C and complex
  1281. a, q0, r0, u0, v0 = self.generate('sqr', p=3)
  1282. a1 = a + np.dot(u0, v0.T.conj())
  1283. q = q0.copy('F')
  1284. r = r0.copy('F')
  1285. u = u0.copy('F')
  1286. v = v0.copy('C')
  1287. # don't overwrite
  1288. q1, r1 = qr_update(q, r, u, v, False)
  1289. check_qr(q1, r1, a1, self.rtol, self.atol)
  1290. check_qr(q, r, a, self.rtol, self.atol)
  1291. q2, r2 = qr_update(q, r, u, v, True)
  1292. check_qr(q2, r2, a1, self.rtol, self.atol)
  1293. # verify the overwriting, no good way to check u and v.
  1294. assert_allclose(q2, q, rtol=self.rtol, atol=self.atol)
  1295. assert_allclose(r2, r, rtol=self.rtol, atol=self.atol)
  1296. def test_empty_inputs(self):
  1297. a, q, r, u, v = self.generate('tall')
  1298. assert_raises(ValueError, qr_update, np.array([]), r, u, v)
  1299. assert_raises(ValueError, qr_update, q, np.array([]), u, v)
  1300. assert_raises(ValueError, qr_update, q, r, np.array([]), v)
  1301. assert_raises(ValueError, qr_update, q, r, u, np.array([]))
  1302. def test_mismatched_shapes(self):
  1303. a, q, r, u, v = self.generate('tall')
  1304. assert_raises(ValueError, qr_update, q, r[1:], u, v)
  1305. assert_raises(ValueError, qr_update, q[:-2], r, u, v)
  1306. assert_raises(ValueError, qr_update, q, r, u[1:], v)
  1307. assert_raises(ValueError, qr_update, q, r, u, v[1:])
  1308. def test_unsupported_dtypes(self):
  1309. dts = ['int8', 'int16', 'int32', 'int64',
  1310. 'uint8', 'uint16', 'uint32', 'uint64',
  1311. 'float16', 'longdouble', 'longcomplex',
  1312. 'bool']
  1313. a, q0, r0, u0, v0 = self.generate('tall')
  1314. for dtype in dts:
  1315. q = q0.real.astype(dtype)
  1316. with np.errstate(invalid="ignore"):
  1317. r = r0.real.astype(dtype)
  1318. u = u0.real.astype(dtype)
  1319. v = v0.real.astype(dtype)
  1320. assert_raises(ValueError, qr_update, q, r0, u0, v0)
  1321. assert_raises(ValueError, qr_update, q0, r, u0, v0)
  1322. assert_raises(ValueError, qr_update, q0, r0, u, v0)
  1323. assert_raises(ValueError, qr_update, q0, r0, u0, v)
  1324. def test_integer_input(self):
  1325. q = np.arange(16).reshape(4, 4)
  1326. r = q.copy() # doesn't matter
  1327. u = q[:, 0].copy()
  1328. v = r[0, :].copy()
  1329. assert_raises(ValueError, qr_update, q, r, u, v)
  1330. def test_check_finite(self):
  1331. a0, q0, r0, u0, v0 = self.generate('tall', p=3)
  1332. q = q0.copy('F')
  1333. q[1,1] = np.nan
  1334. assert_raises(ValueError, qr_update, q, r0, u0[:,0], v0[:,0])
  1335. assert_raises(ValueError, qr_update, q, r0, u0, v0)
  1336. r = r0.copy('F')
  1337. r[1,1] = np.nan
  1338. assert_raises(ValueError, qr_update, q0, r, u0[:,0], v0[:,0])
  1339. assert_raises(ValueError, qr_update, q0, r, u0, v0)
  1340. u = u0.copy('F')
  1341. u[0,0] = np.nan
  1342. assert_raises(ValueError, qr_update, q0, r0, u[:,0], v0[:,0])
  1343. assert_raises(ValueError, qr_update, q0, r0, u, v0)
  1344. v = v0.copy('F')
  1345. v[0,0] = np.nan
  1346. assert_raises(ValueError, qr_update, q0, r0, u[:,0], v[:,0])
  1347. assert_raises(ValueError, qr_update, q0, r0, u, v)
  1348. def test_economic_check_finite(self):
  1349. a0, q0, r0, u0, v0 = self.generate('tall', mode='economic', p=3)
  1350. q = q0.copy('F')
  1351. q[1,1] = np.nan
  1352. assert_raises(ValueError, qr_update, q, r0, u0[:,0], v0[:,0])
  1353. assert_raises(ValueError, qr_update, q, r0, u0, v0)
  1354. r = r0.copy('F')
  1355. r[1,1] = np.nan
  1356. assert_raises(ValueError, qr_update, q0, r, u0[:,0], v0[:,0])
  1357. assert_raises(ValueError, qr_update, q0, r, u0, v0)
  1358. u = u0.copy('F')
  1359. u[0,0] = np.nan
  1360. assert_raises(ValueError, qr_update, q0, r0, u[:,0], v0[:,0])
  1361. assert_raises(ValueError, qr_update, q0, r0, u, v0)
  1362. v = v0.copy('F')
  1363. v[0,0] = np.nan
  1364. assert_raises(ValueError, qr_update, q0, r0, u[:,0], v[:,0])
  1365. assert_raises(ValueError, qr_update, q0, r0, u, v)
  1366. def test_u_exactly_in_span_q(self):
  1367. q = np.array([[0, 0], [0, 0], [1, 0], [0, 1]], self.dtype)
  1368. r = np.array([[1, 0], [0, 1]], self.dtype)
  1369. u = np.array([0, 0, 0, -1], self.dtype)
  1370. v = np.array([1, 2], self.dtype)
  1371. q1, r1 = qr_update(q, r, u, v)
  1372. a1 = np.dot(q, r) + np.outer(u, v.conj())
  1373. check_qr(q1, r1, a1, self.rtol, self.atol, False)
  1374. class TestQRupdate_f(BaseQRupdate):
  1375. dtype = np.dtype('f')
  1376. class TestQRupdate_F(BaseQRupdate):
  1377. dtype = np.dtype('F')
  1378. class TestQRupdate_d(BaseQRupdate):
  1379. dtype = np.dtype('d')
  1380. class TestQRupdate_D(BaseQRupdate):
  1381. dtype = np.dtype('D')
  1382. def test_form_qTu():
  1383. # We want to ensure that all of the code paths through this function are
  1384. # tested. Most of them should be hit with the rest of test suite, but
  1385. # explicit tests make clear precisely what is being tested.
  1386. #
  1387. # This function expects that Q is either C or F contiguous and square.
  1388. # Economic mode decompositions (Q is (M, N), M != N) do not go through this
  1389. # function. U may have any positive strides.
  1390. #
  1391. # Some of these test are duplicates, since contiguous 1d arrays are both C
  1392. # and F.
  1393. q_order = ['F', 'C']
  1394. q_shape = [(8, 8), ]
  1395. u_order = ['F', 'C', 'A'] # here A means is not F not C
  1396. u_shape = [1, 3]
  1397. dtype = ['f', 'd', 'F', 'D']
  1398. for qo, qs, uo, us, d in \
  1399. itertools.product(q_order, q_shape, u_order, u_shape, dtype):
  1400. if us == 1:
  1401. check_form_qTu(qo, qs, uo, us, 1, d)
  1402. check_form_qTu(qo, qs, uo, us, 2, d)
  1403. else:
  1404. check_form_qTu(qo, qs, uo, us, 2, d)
  1405. def check_form_qTu(q_order, q_shape, u_order, u_shape, u_ndim, dtype):
  1406. np.random.seed(47)
  1407. if u_shape == 1 and u_ndim == 1:
  1408. u_shape = (q_shape[0],)
  1409. else:
  1410. u_shape = (q_shape[0], u_shape)
  1411. dtype = np.dtype(dtype)
  1412. if dtype.char in 'fd':
  1413. q = np.random.random(q_shape)
  1414. u = np.random.random(u_shape)
  1415. elif dtype.char in 'FD':
  1416. q = np.random.random(q_shape) + 1j*np.random.random(q_shape)
  1417. u = np.random.random(u_shape) + 1j*np.random.random(u_shape)
  1418. else:
  1419. ValueError("form_qTu doesn't support this dtype")
  1420. q = np.require(q, dtype, q_order)
  1421. if u_order != 'A':
  1422. u = np.require(u, dtype, u_order)
  1423. else:
  1424. u, = make_strided((u.astype(dtype),))
  1425. rtol = 10.0 ** -(np.finfo(dtype).precision-2)
  1426. atol = 2*np.finfo(dtype).eps
  1427. expected = np.dot(q.T.conj(), u)
  1428. res = _decomp_update._form_qTu(q, u)
  1429. assert_allclose(res, expected, rtol=rtol, atol=atol)