1# 2# This file is part of CasADi. 3# 4# CasADi -- A symbolic framework for dynamic optimization. 5# Copyright (C) 2010-2014 Joel Andersson, Joris Gillis, Moritz Diehl, 6# K.U. Leuven. All rights reserved. 7# Copyright (C) 2011-2014 Greg Horn 8# 9# CasADi is free software; you can redistribute it and/or 10# modify it under the terms of the GNU Lesser General Public 11# License as published by the Free Software Foundation; either 12# version 3 of the License, or (at your option) any later version. 13# 14# CasADi is distributed in the hope that it will be useful, 15# but WITHOUT ANY WARRANTY; without even the implied warranty of 16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17# Lesser General Public License for more details. 18# 19# You should have received a copy of the GNU Lesser General Public 20# License along with CasADi; if not, write to the Free Software 21# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 22# 23# 24from casadi import * 25import casadi as c 26import numpy 27from numpy import eye, linalg, arange, matrix 28import unittest 29from types import * 30from helpers import * 31import numpy 32from itertools import * 33 34class Matrixtests(casadiTestCase): 35 def test_constructorlol(self): 36 self.message("List of list constructor") 37 a=DM(array([[1,2,3],[4,5,6],[7,8,9]])) 38 b=DM([[1,2,3],[4,5,6],[7,8,9]]) 39 self.checkarray(a,b,"List of list constructor") 40 41 def test_sum(self): 42 self.message("sum") 43 D=DM([[1,2,3],[4,5,6],[7,8,9]]) 44 self.checkarray(c.sum1(D),array([[12,15,18]]),'sum()') 45 self.checkarray(c.sum2(D),array([[6,15,24]]).T,'sum()') 46 47 48 def test_inv(self): 49 self.message("Matrix inverse") 50 a = DM([[1,2],[1,3]]) 51 self.checkarray(mtimes(c.inv(a),a),eye(2),"DM inverse") 52 53 def test_trans(self): 54 self.message("trans") 55 a = DM(0,1) 56 b = a.T 57 self.assertEqual(b.size1(),1) 58 self.assertEqual(b.size2(),0) 59 60 def test_vertcat(self): 61 self.message("vertcat") 62 A = DM.ones(2,3) 63 B = DM(4,3) 64 C = vertcat(A,B) 65 66 self.checkarray(C.shape,(6,3),"vertcat shape") 67 self.assertEqual(C.nnz(),A.nnz(),"vertcat size") 68 69 self.assertRaises(RuntimeError,lambda : horzcat(A,B)) 70 71 def test_horzcat(self): 72 self.message("horcat") 73 A = DM.ones(3,2) 74 B = DM(3,4) 75 C = horzcat(A,B) 76 77 self.checkarray(C.shape,(3,6),"horzcat shape") 78 self.assertEqual(C.nnz(),A.nnz(),"vertcat size") 79 80 self.assertRaises(RuntimeError,lambda : vertcat(A,B)) 81 82 def test_diagcat(self): 83 84 x = MX.sym("x",2,2) 85 y = MX.sym("y",Sparsity.lower(3)) 86 z = MX.sym("z",4,2) 87 88 L = [x,y,z] 89 90 fMX = Function("fMX", L,[diagcat(*L)]) 91 92 LSX = [ SX.sym("",i.sparsity()) for i in L ] 93 fSX = Function("fSX", LSX,[diagcat(*LSX)]) 94 95 for f in [fMX,fSX]: 96 for i in range(3): 97 f_in[i]=list(range(f.nnz_in(i))) 98 99 self.checkfunction(fMX,fSX) 100 101 def test_veccat(self): 102 self.message("vecccat") 103 A = DM(2,3) 104 A[0,1] = 2 105 A[1,0] = 1 106 A[1,2] = 3 107 B = DM(3,1) 108 B[0,0] = 4 109 B[1,0] = 5 110 B[2,0] = 6 111 C = veccat(*[A,B]) 112 113 self.checkarray(C.shape,(9,1),"veccat shape") 114 self.assertEqual(C.nnz(),A.nnz()+B.nnz(),"veccat size") 115 116 self.checkarray(tuple(C.nonzeros()),tuple(arange(1,7)),"numbers shape") 117 118 def test_slicestepnegative(self): 119 self.message("Slice step negative") 120 a1 = [1,2,3,4,5] 121 a2 = DM(a1) 122 123 self.checkarray(a2[0:4:-1,0],DM(a1[0:4:-1])) # gives empty set 124 self.checkarray(a2[4:0:-1,0],DM(a1[4:0:-1])) # gives [5, 4, 3, 2] 125 self.checkarray(a2[0:4:-2,0],DM(a1[0:4:-2])) # gives empty set 126 self.checkarray(a2[4:0:-2,0],DM(a1[4:0:-2])) # gives [5, 4, 3, 2] 127 self.checkarray(a2[1:4:-2,0],DM(a1[1:4:-2])) # gives empty set 128 self.checkarray(a2[4:1:-2,0],DM(a1[4:1:-2])) # gives [5, 4, 3, 2] 129 self.checkarray(a2[0:3:-2,0],DM(a1[0:3:-2])) # gives empty set 130 self.checkarray(a2[3:0:-2,0],DM(a1[3:0:-2])) # gives [5, 4, 3, 2] 131 self.checkarray(a2[::-1,0],DM(a1[::-1])) # gives [5, 4, 3, 2, 1] 132 self.checkarray(a2[::1,0],DM(a1[::1])) # gives [1,2,3,4,5] 133 self.checkarray(a2[2::-1,0],DM(a1[2::-1])) # gives [3,2,1] 134 self.checkarray(a2[:2:-1,0],DM(a1[:2:-1])) # gives [5,4] 135 self.checkarray(a2[-1::-1,0],DM(a1[-1::-1])) # gives [5,4,3,2,1] 136 self.checkarray(a2[:-1:-1,0],DM(a1[:-1:-1])) # gives [] 137 138 def test_indexingOutOfBounds(self): 139 self.message("Indexing out of bounds") 140 y = DM.zeros(4, 5) 141 self.assertRaises(RuntimeError,lambda : y[12,0] ) 142 self.assertRaises(RuntimeError,lambda : y[12,12] ) 143 self.assertRaises(RuntimeError,lambda : y[0,12] ) 144 self.assertRaises(RuntimeError,lambda : y[12,:] ) 145 self.assertRaises(RuntimeError,lambda : y[12:15,0] ) 146 self.assertRaises(RuntimeError,lambda : y[:,12] ) 147 self.assertRaises(RuntimeError,lambda : y[0,12:15] ) 148 y[-1,2] 149 self.assertRaises(RuntimeError,lambda : y[-12,2] ) 150 y[-3:-1,2] 151 self.assertRaises(RuntimeError,lambda : y[-12:-9,2] ) 152 153 def test(): 154 y[12,0] = 0 155 self.assertRaises(RuntimeError,test) 156 def test(): 157 y[12,12] = 0 158 self.assertRaises(RuntimeError,test) 159 def test(): 160 y[0,12] = 0 161 self.assertRaises(RuntimeError,test) 162 def test(): 163 y[12,:] = 0 164 self.assertRaises(RuntimeError,test) 165 def test(): 166 y[12:15,0] = 0 167 self.assertRaises(RuntimeError,test) 168 def test(): 169 y[:,12] = 0 170 self.assertRaises(RuntimeError,test) 171 def test(): 172 y[0,12:15] = 0 173 self.assertRaises(RuntimeError,test) 174 y[-1,2] = 0 175 def test(): 176 y[-12,2] = 0 177 self.assertRaises(RuntimeError,test) 178 y[-3:-1,2] = 0 179 def test(): 180 y[-12:-9,2] = 0 181 self.assertRaises(RuntimeError,test) 182 183 def huge_slice(self): 184 self.message("huge slice") 185 a = SX.sym("a",Sparsity.diag(50000)) 186 187 a[:,:] 188 189 def test_nonmonotonous_indexing(self): 190 self.message("non-monotonous indexing") 191 # Regression test for #354 192 A = DM([[1,2,3],[4,5,6],[7,8,9]]) 193 B = A[[0,2,1],0] 194 self.checkarray(DM([1,7,4]),B,"non-monotonous") 195 196 B = A[0,[0,2,1]] 197 self.checkarray(DM([1,3,2]).T,B,"non-monotonous") 198 199 def test_IM_indexing(self): 200 self.message("IM") 201 A = DM(2,2) 202 A[0,0] = 1 203 A[1,1] = 3 204 A[0,1] = 2 205 A[1,0] = 4 206 207 208 B = DM([1,2,3,4,5]) 209 210 B_ = B[A] 211 212 self.checkarray(B_,DM([[2,3],[5,4]]),"matrix indexing") 213 214 B[A] = DM([[1,2],[3,4]]) 215 216 self.checkarray(B,DM([1,1,2,4,3]),"matrix indexing assignement") 217 218 #B[A].set(DM([[10,20],[30,40]])) 219 220 #self.checkarray(B,DM([1,10,20,40,30]),"Imatrix indexing setting") 221 222 B = DM([1,2,3,4,5]) 223 224 B_ = B[A] 225 226 self.checkarray(array(B_),DM([[2,3],[5,4]]),"Imatrix indexing") 227 228 B[A] = DM([[1,2],[3,4]]) 229 230 self.checkarray(array(B),DM([1,1,2,4,3]),"Imatrix indexing assignement") 231 232 B = DM(5,1) 233 234 self.assertRaises(Exception, lambda : B.nz[A]) 235 236 def test_sparsity_indexing(self): 237 self.message("sparsity") 238 239 B = DM([[1,2,3,4,5],[6,7,8,9,10]]) 240 241 A = DM([[1,1,0,0,0],[0,0,1,0,0]]) 242 A = sparsify(A) 243 sp = A.sparsity() 244 245 246 self.checkarray(B[sp],DM([[1,2,0,0,0],[0,0,8,0,0]]),"sparsity indexing") 247 248 B[sp] = -4 249 250 self.checkarray(B,DM([[-4,-4,3,4,5],[6,7,-4,9,10]]),"sparsity indexing assignement") 251 252 B = DM([[1,2,3,4,5],[6,7,8,9,10]]) 253 254 B[sp] = 2*B 255 256 self.checkarray(B,DM([[2,4,3,4,5],[6,7,16,9,10]]),"Imatrix indexing assignement") 257 258 self.assertRaises(Exception, lambda : B[Sparsity.dense(4,4)]) 259 260 def test_index_setting(self): 261 self.message("index setting") 262 B = DM([1,2,3,4,5]) 263 264 B[0] = 8 265 self.checkarray(B,DM([8,2,3,4,5]),"index setting") 266 B[1,0] = 4 267 self.checkarray(B,DM([8,4,3,4,5]),"index setting") 268 B[:,0] = 7 269 self.checkarray(B,DM([7,7,7,7,7]),"index setting") 270 #B[0].set(3) 271 #self.checkarray(B,DM([3,7,7,7,7]),"index setting") 272 #B[0].setAll(4) 273 #self.checkarray(B,DM([4,4,4,4,4]),"index setting") 274 275 276 def test_issue298(self): 277 self.message("Issue #298") 278 a = DM(4,1) 279 b = c.reshape(a,2,2) 280 self.assertEqual(type(a),type(b)) 281 282 def test_det(self): 283 self.message("Determinant") 284 npy_det = numpy.linalg.det 285 286 a = DM(1,1) 287 a[0,0] = 5 288 self.checkarray(det(a),npy_det(a),"det()") 289 290 a = DM(5,5) 291 for i in range(5): 292 a[i,i] = i+1 293 294 self.checkarray(det(a),npy_det(a),"det()") 295 296 a = DM(5,5) 297 for i in range(4): 298 a[i,i] = i+1 299 a[0,4] = 3 300 a[4,0] = 7 301 302 self.checkarray(det(a),npy_det(a),"det()") 303 304 a = DM(5,5) 305 for i in range(5): 306 for j in range(5): 307 a[i,j] = i+j 308 309 self.checkarray(det(a),npy_det(a),"det()") 310 311 a = DM(5,5) 312 for i in range(4): 313 for j in range(5): 314 a[i,j] = i+j 315 316 self.checkarray(det(a),npy_det(a),"det()") 317 318 a = DM(5,5) 319 for i in range(5): 320 for j in range(4): 321 a[i,j] = i+j 322 323 self.checkarray(det(a),npy_det(a),"det()") 324 325 a = DM(5,5) 326 for i in range(4): 327 for j in range(5): 328 a[i,j] = i+j 329 a[4,1] = 12 330 331 self.checkarray(det(a),npy_det(a),"det()") 332 333 a = DM(5,5) 334 for i in range(5): 335 for j in range(4): 336 a[i,j] = i+j 337 a[1,4] = 12 338 339 self.checkarray(det(a),npy_det(a),"det()") 340 341 a = DM(5,5) 342 for i in range(4): 343 for j in range(5): 344 a[i,j] = i+j 345 a[4,2] = 12 346 347 self.checkarray(det(a),npy_det(a),"det()") 348 349 a = DM(5,5) 350 for i in range(5): 351 for j in range(4): 352 a[i,j] = i+j 353 a[2,4] = 12 354 355 self.checkarray(det(a),npy_det(a),"det()") 356 357 a = DM(50,50) 358 for i in range(50): 359 a[i,i] = i+1 360 361 self.checkarray(det(a)/npy_det(a),1,"det()") 362 363 @skip(not GlobalOptions.getSimplificationOnTheFly()) 364 def test_inv_sparsity(self): 365 self.message("sparsity pattern of inverse") 366 367 n = 8 368 369 sp = Sparsity.lower(n) 370 371 x = SX(sp,vertcat(*[SX.sym("a%d" % i) for i in range(sp.nnz())])) 372 373 374 x_ = DM.ones(x.sparsity()) 375 376 I_ = DM.ones(inv(x).sparsity()) 377 378 # For a reducible matrix, struct(A^(-1)) = struct(A) 379 self.checkarray(x_,I_,"inv") 380 381 sp = Sparsity.lower(n) 382 383 x = SX.sym("a", sp) 384 x[0,n-1] = 1 385 386 387 I_ = DM.ones(inv(x).sparsity()) 388 389 # An irreducible matrix has a dense inverse in general 390 self.checkarray(DM.ones(n,n),I_,"inv") 391 392 x = SX.sym("a", sp) 393 x[0,int(n/2)] = 1 394 395 s_ = DM.ones(sp) 396 s_[:,:int(n/2)+1] = 1 397 398 I_ = DM.ones(inv_minor(x).sparsity()) 399 400 s_ = densify(s_) 401 T_ = densify(I_) 402 # An irreducible matrix does not have to be dense per se 403 self.checkarray(s_,I_,"inv") 404 405 def test_mtimes(self): 406 A = DM.ones((4,3)) 407 B = DM.ones((3,8)) 408 C = DM.ones((8,7)) 409 410 self.assertRaises(RuntimeError,lambda : mtimes([])) 411 412 D = mtimes([A]) 413 414 self.assertEqual(D.shape[0],4) 415 self.assertEqual(D.shape[1],3) 416 417 D = mtimes([A,B]) 418 419 self.assertEqual(D.shape[0],4) 420 self.assertEqual(D.shape[1],8) 421 422 D = mtimes([A,B,C]) 423 424 self.assertEqual(D.shape[0],4) 425 self.assertEqual(D.shape[1],7) 426 427 def test_remove(self): 428 self.message("remove") 429 B = DM([[1,2,3,4],[5,6,7,8],[9,10,11,12],[13,14,15,16],[17,18,19,20]]) 430 431 A = DM(B) 432 A.remove([],[]) 433 self.checkarray(A, B,"remove nothing") 434 435 A = DM(B) 436 A.remove([],[1]) 437 self.checkarray(A, DM([[1,3,4],[5,7,8],[9,11,12],[13,15,16],[17,19,20]]),"remove a column") 438 439 A = DM(B) 440 A.remove([0,3],[1]) 441 self.checkarray(A, DM([[5,7,8],[9,11,12],[17,19,20]]),"remove a column and two rows ") 442 443 def test_comparisons(self): 444 m = DM 445 A = m([[5,4],[2,1]]) 446 447 for c in [6,6.0,DM([6]),np.array([6]),matrix(6)]: 448 self.checkarray(A<=c,m([[1,1],[1,1]]),"<=") 449 self.checkarray(A<c,m([[1,1],[1,1]]),"<") 450 self.checkarray(A>c,m([[0,0],[0,0]]),">") 451 self.checkarray(A>=c,m([[0,0],[0,0]]),">=") 452 self.checkarray(A==c,m([[0,0],[0,0]]),"==") 453 self.checkarray(A!=c,m([[1,1],[1,1]]),"!=") 454 455 self.checkarray(c>=A,m([[1,1],[1,1]]),"<=") 456 self.checkarray(c>A,m([[1,1],[1,1]]),"<") 457 self.checkarray(c<A,m([[0,0],[0,0]]),">") 458 self.checkarray(c<=A,m([[0,0],[0,0]]),">=") 459 if args.known_bugs or not isinstance(c,matrix): 460 self.checkarray(c==A,m([[0,0],[0,0]]),"==") 461 self.checkarray(c!=A,m([[1,1],[1,1]]),"!=") 462 463 for c in [5,5.0,DM([5]),np.array([5])]+([matrix(5)] if check_matrix else []): 464 self.checkarray(A<=c,m([[1,1],[1,1]]),"<=") 465 self.checkarray(A<c,m([[0,1],[1,1]]),"<") 466 self.checkarray(A>c,m([[0,0],[0,0]]),">") 467 self.checkarray(A>=c,m([[1,0],[0,0]]),">=") 468 self.checkarray(A==c,m([[1,0],[0,0]]),"==") 469 self.checkarray(A!=c,m([[0,1],[1,1]]),"!=") 470 471 self.checkarray(c>=A,m([[1,1],[1,1]]),"<=") 472 self.checkarray(c>A,m([[0,1],[1,1]]),"<") 473 self.checkarray(c<A,m([[0,0],[0,0]]),">") 474 self.checkarray(c<=A,m([[1,0],[0,0]]),">=") 475 if args.known_bugs or not isinstance(c,matrix): 476 self.checkarray(c==A,m([[1,0],[0,0]]),"==") 477 self.checkarray(c!=A,m([[0,1],[1,1]]),"!=") 478 479 for c in [4,4.0,DM([4]),np.array([4]),matrix(4)]: 480 self.checkarray(A<=c,m([[0,1],[1,1]]),"<=") 481 self.checkarray(A<c,m([[0,0],[1,1]]),"<") 482 self.checkarray(A>c,m([[1,0],[0,0]]),">") 483 self.checkarray(A>=c,m([[1,1],[0,0]]),">=") 484 if args.known_bugs or not isinstance(c,matrix): 485 self.checkarray(A==c,m([[0,1],[0,0]]),"==") 486 self.checkarray(A!=c,m([[1,0],[1,1]]),"!=") 487 488 self.checkarray(c>=A,m([[0,1],[1,1]]),"<=") 489 self.checkarray(c>A,m([[0,0],[1,1]]),"<") 490 self.checkarray(c<A,m([[1,0],[0,0]]),">") 491 self.checkarray(c<=A,m([[1,1],[0,0]]),">=") 492 if args.known_bugs or not isinstance(c,matrix): 493 self.checkarray(c==A,m([[0,1],[0,0]]),"==") 494 self.checkarray(c!=A,m([[1,0],[1,1]]),"!=") 495 496 for c in [1,1.0,DM([1]),np.array([1]),matrix(1)]: 497 self.checkarray(A<=c,m([[0,0],[0,1]]),"<=") 498 self.checkarray(A<c,m([[0,0],[0,0]]),"<") 499 self.checkarray(A>c,m([[1,1],[1,0]]),">") 500 self.checkarray(A>=c,m([[1,1],[1,1]]),">=") 501 self.checkarray(A==c,m([[0,0],[0,1]]),"==") 502 self.checkarray(A!=c,m([[1,1],[1,0]]),"!=") 503 504 self.checkarray(c>=A,m([[0,0],[0,1]]),"<=") 505 self.checkarray(c>A,m([[0,0],[0,0]]),"<") 506 self.checkarray(c<A,m([[1,1],[1,0]]),">") 507 self.checkarray(c<=A,m([[1,1],[1,1]]),">=") 508 if args.known_bugs or not isinstance(c,matrix): 509 self.checkarray(c==A,m([[0,0],[0,1]]),"==") 510 self.checkarray(c!=A,m([[1,1],[1,0]]),"!=") 511 512 for c in [0,DM([0]),np.array([0]),matrix(0)]: 513 self.checkarray(A<=c,m([[0,0],[0,0]]),"<=") 514 self.checkarray(A<c,m([[0,0],[0,0]]),"<") 515 self.checkarray(A>c,m([[1,1],[1,1]]),">") 516 self.checkarray(A>=c,m([[1,1],[1,1]]),">=") 517 self.checkarray(A==c,m([[0,0],[0,0]]),"==") 518 self.checkarray(A!=c,m([[1,1],[1,1]]),"!=") 519 520 self.checkarray(c>=A,m([[0,0],[0,0]]),"<=") 521 self.checkarray(c>A,m([[0,0],[0,0]]),"<") 522 self.checkarray(c<A,m([[1,1],[1,1]]),">") 523 self.checkarray(c<=A,m([[1,1],[1,1]]),">=") 524 if args.known_bugs or not isinstance(c,matrix): 525 self.checkarray(c==A,m([[0,0],[0,0]]),"==") 526 self.checkarray(c!=A,m([[1,1],[1,1]]),"!=") 527 528 def test_truth(self): 529 self.assertTrue(bool(DM([1]))) 530 self.assertFalse(bool(DM([0]))) 531 self.assertTrue(bool(DM([0.2]))) 532 self.assertTrue(bool(DM([-0.2]))) 533 self.assertRaises(Exception, lambda : bool(DM([2.0,3]))) 534 self.assertRaises(Exception, lambda : bool(DM())) 535 536 def test_listslice(self): 537 def check(d,rowbase,colbase): 538 for col in permutations(colbase): 539 for row in permutations(rowbase): 540 r = DM.zeros(len(row),len(col)) 541 for i,ii in enumerate(row): 542 for j,jj in enumerate(col): 543 r[i,j] = d[ii,jj] 544 self.checkarray(d[row,col],r,"%s[%s,%s]" % (repr(d),str(row),str(col))) 545 546 547 # get1 548 check(DM(Sparsity.dense(3,3),list(range(3*3))),[0,1,2],[0,1,2]) 549 check(DM(Sparsity.dense(4,4),list(range(4*4))),[0,1,3],[0,2,3]) 550 check(DM(Sparsity.dense(3,3),list(range(3*3))),[0,0,1],[0,0,1]) 551 check(DM(Sparsity.dense(3,3),list(range(3*3))),[0,0,2],[0,0,2]) 552 check(DM(Sparsity.dense(3,3),list(range(3*3))),[1,1,2],[1,1,2]) 553 554 sp = Sparsity.lower(4) 555 d = DM(sp,list(range(sp.nnz()))) 556 check(d,[0,1,3],[0,2,3]) 557 check(d.T,[0,1,3],[0,2,3]) 558 559 sp = Sparsity.rowcol([0,1,2],[0,1],4,4) 560 d = DM(sp,list(range(sp.nnz()))) 561 check(d,[0,3],[0,2]) 562 563 # get2 564 check(DM(Sparsity.dense(2,2),list(range(2*2))),[0,0,0],[0,0,0]) 565 check(DM(Sparsity.dense(2,2),list(range(2*2))),[0,0,1],[0,0,1]) 566 check(DM(Sparsity.dense(2,2),list(range(2*2))),[1,1,0],[1,1,0]) 567 check(DM(Sparsity.dense(2,2),list(range(2*2))),[1,1,1],[1,1,1]) 568 569 sp = Sparsity.lower(3) 570 d = DM(sp,list(range(sp.nnz()))) 571 check(d,[0,1,2],[0,1,2]) 572 check(d.T,[0,1,2],[0,1,2]) 573 574 sp = Sparsity.rowcol([0,2],[0,1],4,4) 575 d = DM(sp,list(range(sp.nnz()))) 576 check(d,[0,1,3],[0,2,3]) 577 578 def test_sparsesym(self): 579 # feature removed in 73f407e 580 return 581 self.message("sparsesym") 582 D = DM([[1,2,-3],[2,-1,0],[-3,0,5]]) 583 D = sparsify(D) 584 585 i = D.getSym() 586 #self.checkarray(list(i),[1,2,-1,-3,5]) 587 A = 2*D 588 A.setSym(i) 589 self.checkarray(A,D) 590 591 def test_diagcat(self): 592 self.message("diagcat") 593 C = diagcat(*[DM([[-1.4,-3.2],[-3.2,-28]]),DM([[15,-12,2.1],[-12,16,-3.8],[2.1,-3.8,15]]),1.8,-4.0]) 594 r = DM([[-1.4,-3.2,0,0,0,0,0],[-3.2,-28,0,0,0,0,0],[0,0,15,-12,2.1,0,0],[0,0,-12,16,-3.8,0,0],[0,0,2.1,-3.8,15,0,0],[0,0,0,0,0,1.8,0],[0,0,0,0,0,0,-4]]) 595 r = sparsify(r) 596 self.checkarray(C,r) 597 598 def test_diag_sparse(self): 599 self.message("diag sparse") 600 601 for n in [[0,1,0,0,2,3,4,5,6,0],[1,2,3,0],[0,1,2,3]]: 602 d = DM(n) 603 D = DM(n) 604 d = sparsify(d) 605 m = c.diag(d) 606 M = sparsify(c.diag(D)) 607 608 self.checkarray(m.sparsity().colind(),M.sparsity().colind()) 609 self.checkarray(m.sparsity().row(),M.sparsity().row()) 610 611 def test_sprank(self): 612 self.message("sprank") 613 614 a = DM([[1,0,0],[0,1,0],[0,0,1]]) 615 a = sparsify(a) 616 self.assertEqual(sprank(a),3) 617 618 a = DM([[1,0,0],[0,0,0],[0,0,1]]) 619 a = sparsify(a) 620 self.assertEqual(sprank(a),2) 621 622 a = DM([[0,0,0],[0,0,0],[0,0,1]]) 623 a = sparsify(a) 624 self.assertEqual(sprank(a),1) 625 626 a = DM([[0,0,0],[0,0,0],[0,0,0]]) 627 a = sparsify(a) 628 self.assertEqual(sprank(a),0) 629 630 self.assertEqual(sprank(DM.ones(1,3)),1) 631 self.assertEqual(sprank(DM.ones(3,1)),1) 632 self.assertEqual(sprank(DM.ones(2,3)),2) 633 self.assertEqual(sprank(DM.ones(3,2)),2) 634 self.assertEqual(sprank(DM.ones(3,3)),3) 635 self.assertEqual(sprank(DM.ones(3,3)),3) 636 637 A = DM(6,4) 638 A[0,0] = 1 639 A[1,2] = 1 640 A[2,2] = 1 641 A[5,3] = 1 642 643 self.assertEqual(sprank(A),3) 644 645 def test_cross(self): 646 self.message("cross products") 647 648 crossc = c.cross 649 650 self.checkarray(crossc(DM([1,0,0]),DM([0,1,0])),DM([0,0,1])) 651 652 self.checkarray(crossc(DM([1.1,1.3,1.7]),DM([2,3,13])),DM([11.8,-10.9,0.7])) 653 self.checkarray(crossc(DM([1.1,1.3,1.7]).T,DM([2,3,13]).T),DM([11.8,-10.9,0.7]).T) 654 655 self.checkarray(crossc(DM([[1.1,1.3,1.7],[1,0,0],[0,0,1],[4,5,6]]),DM([[2,3,13],[0,1,0],[0,0,1],[1,0,1]])),DM([[11.8,-10.9,0.7],[0,0,1],[0,0,0],[5,2,-5]])) 656 self.checkarray(crossc(DM([[1.1,1.3,1.7],[1,0,0],[0,0,1],[4,5,6]]).T,DM([[2,3,13],[0,1,0],[0,0,1],[1,0,1]]).T),DM([[11.8,-10.9,0.7],[0,0,1],[0,0,0],[5,2,-5]]).T) 657 658 self.checkarray(crossc(DM([[1.1,1.3,1.7],[1,0,0],[0,0,1],[4,5,6]]),DM([[2,3,13],[0,1,0],[0,0,1],[1,0,1]]),2),DM([[11.8,-10.9,0.7],[0,0,1],[0,0,0],[5,2,-5]])) 659 660 self.checkarray(crossc(DM([[1.1,1.3,1.7],[1,0,0],[0,0,1],[4,5,6]]).T,DM([[2,3,13],[0,1,0],[0,0,1],[1,0,1]]).T,1),DM([[11.8,-10.9,0.7],[0,0,1],[0,0,0],[5,2,-5]]).T) 661 662 def test_is_regular(self): 663 self.assertTrue(DM([1,2]).is_regular()) 664 self.assertFalse(DM([1,inf]).is_regular()) 665 self.assertFalse(DM.nan(2).is_regular()) 666 667 def test_sizes(self): 668 self.assertEqual(Sparsity.diag(10).nnz_diag(),10) 669 self.assertEqual(Sparsity.diag(10).nnz_upper(),10) 670 self.assertEqual(Sparsity.diag(10).nnz_lower(),10) 671 self.assertEqual(Sparsity.dense(10,10).nnz_lower(),10*11/2) 672 self.assertEqual(Sparsity.dense(10,10).nnz_upper(),10*11/2) 673 self.assertEqual(Sparsity.dense(10,10).nnz_diag(),10) 674 675 self.assertEqual(sparsify(DM([[1,1,0],[1,0,1],[0,0,0]])).nnz_diag(),1) 676 self.assertEqual(sparsify(DM([[1,1,0],[1,0,1],[0,0,0]])).nnz_lower(),2) 677 self.assertEqual(sparsify(DM([[1,1,0],[1,0,1],[0,0,0]])).nnz_upper(),3) 678 679 def test_tril2symm(self): 680 a = DM(Sparsity.upper(3),list(range(Sparsity.upper(3).nnz()))).T 681 s = tril2symm(a) 682 self.checkarray(s,DM([[0,1,3],[1,2,4],[3,4,5]])) 683 684 with self.assertRaises(Exception): 685 tril2symm(DM.ones(5,3)) 686 687 print(DM.ones(5,5).nnz_upper()-DM.ones(5,5).nnz_diag()) 688 689 with self.assertRaises(Exception): 690 tril2symm(DM.ones(5,5)) 691 692 def test_not_null(self): 693 x = MX.sym('x',3,1) 694 sp = Sparsity.upper(2) 695 MX(sp,x) 696 697 def test_segfault(self): 698 x = MX.sym('x',10,1) 699 sp = Sparsity.upper(2) 700 y = triu2symm(MX(sp,x[1:4])) 701 f = Function("f", [x],[y]) 702 703 def test_append_empty(self): 704 a = vertcat(DM(0,0),DM(0,2)) 705 706 self.assertEqual(a.size1(),0) 707 self.assertEqual(a.size2(),2) 708 709 a = vertcat(DM(0,0),DM(2,0),DM(3,0)) 710 711 self.assertEqual(a.size1(),5) 712 self.assertEqual(a.size2(),0) 713 714 def test_vertcat_empty(self): 715 a = DM(0,2) 716 v = vertcat(a,a) 717 718 self.assertEqual(v.size1(),0) 719 self.assertEqual(v.size2(),2) 720 721 a = DM(2,0) 722 v = vertcat(a,a) 723 724 self.assertEqual(v.size1(),4) 725 self.assertEqual(v.size2(),0) 726 727 def test_vertsplit(self): 728 a = DM(Sparsity.upper(5),list(range(int(5*6/2)))).T 729 v = vertsplit(a,[0,2,4,5]) 730 731 self.assertEqual(len(v),3) 732 self.checkarray(v[0],DM([[0,0,0,0,0],[1,2,0,0,0]])) 733 self.checkarray(v[1],DM([[3,4,5,0,0],[6,7,8,9,0]])) 734 self.checkarray(v[2],DM([[10,11,12,13,14]])) 735 736 v = vertsplit(a) 737 self.assertEqual(len(v),a.size1()) 738 self.checkarray(v[0],DM([[0,0,0,0,0]])) 739 self.checkarray(v[1],DM([[1,2,0,0,0]])) 740 self.checkarray(v[2],DM([[3,4,5,0,0]])) 741 self.checkarray(v[3],DM([[6,7,8,9,0]])) 742 self.checkarray(v[4],DM([[10,11,12,13,14]])) 743 744 v = vertsplit(a,[0,2,4,5]) 745 self.assertEqual(len(v),3) 746 self.checkarray(v[0],DM([[0,0,0,0,0],[1,2,0,0,0]])) 747 self.checkarray(v[1],DM([[3,4,5,0,0],[6,7,8,9,0]])) 748 self.checkarray(v[2],DM([[10,11,12,13,14]])) 749 750 v = vertsplit(a,[0,0,3,5]) 751 self.assertEqual(len(v),3) 752 self.assertEqual(v[0].size1(),0) 753 self.assertEqual(v[0].size2(),5) 754 self.checkarray(v[1],DM([[0,0,0,0,0],[1,2,0,0,0],[3,4,5,0,0]])) 755 self.checkarray(v[2],DM([[6,7,8,9,0],[10,11,12,13,14]])) 756 757 def test_horzsplit(self): 758 a = DM(Sparsity.upper(5),list(range(int(5*6/2)))).T 759 v = horzsplit(a,[0,2,4,5]) 760 761 self.assertEqual(len(v),3) 762 self.checkarray(v[0],DM([[0,0],[1,2],[3,4],[6,7],[10,11]])) 763 self.checkarray(v[1],DM([[0,0],[0,0],[5,0],[8,9],[12,13]])) 764 self.checkarray(v[2],DM([[0],[0],[0],[0],[14]])) 765 766 v = horzsplit(a) 767 self.assertEqual(len(v),a.size1()) 768 self.checkarray(v[0],DM([0,1,3,6,10])) 769 self.checkarray(v[1],DM([0,2,4,7,11])) 770 self.checkarray(v[2],DM([0,0,5,8,12])) 771 self.checkarray(v[3],DM([0,0,0,9,13])) 772 self.checkarray(v[4],DM([0,0,0,0,14])) 773 774 v = horzsplit(a,[0,2,4,5]) 775 self.assertEqual(len(v),3) 776 self.checkarray(v[0],DM([[0,0],[1,2],[3,4],[6,7],[10,11]])) 777 self.checkarray(v[1],DM([[0,0],[0,0],[5,0],[8,9],[12,13]])) 778 self.checkarray(v[2],DM([[0],[0],[0],[0],[14]])) 779 780 v = horzsplit(a,[0,0,3,5]) 781 self.assertEqual(len(v),3) 782 self.assertEqual(v[0].size1(),5) 783 self.assertEqual(v[0].size2(),0) 784 self.checkarray(v[1],DM([[0,0,0],[1,2,0],[3,4,5],[6,7,8],[10,11,12]])) 785 self.checkarray(v[2],DM([[0,0],[0,0],[0,0],[9,0],[13,14]])) 786 787 def test_blocksplit(self): 788 a = DM(Sparsity.upper(5),list(range(int(5*6/2)))).T 789 v = blocksplit(a,[0,2,4,5],[0,1,3,5]) 790 791 self.checkarray(v[0][0],DM([0,1])) 792 self.checkarray(v[0][1],DM([[0,0],[2,0]])) 793 self.checkarray(v[1][0],DM([3,6])) 794 self.checkarray(blockcat(v),a) 795 796 def test_solve(self): 797 import random 798 799 spA = [ 800 Sparsity.dense(1,1) 801 ] 802 803 for n in range(2,5): 804 spA+= [ 805 Sparsity.diag(n), 806 Sparsity.dense(n,n), 807 Sparsity.lower(n), 808 Sparsity.lower(n).T, 809 Sparsity.banded(n,1), 810 diagcat(*[Sparsity.diag(n),Sparsity.dense(n,n)]), 811 diagcat(*[Sparsity.diag(n),Sparsity.lower(n)]), 812 diagcat(*[Sparsity.diag(n),Sparsity.lower(n).T]), 813 diagcat(*[Sparsity.lower(n),Sparsity.lower(n).T]), 814 Sparsity.diag(n)+Sparsity.rowcol([0],[n-1],n,n), 815 Sparsity.diag(n)+Sparsity.rowcol([0,n-1],[n-1,0],n,n), 816 Sparsity.diag(n)+Sparsity.triplet(n,n,[0],[n-1]), 817 Sparsity.diag(n)+Sparsity.triplet(n,n,[0,n-1],[n-1,0]), 818 ] 819 820 for sA in spA: 821 822 random.seed(1) 823 a = DM(sA,[random.random() for i in range(sA.nnz())]) 824 A = SX.sym("a",a.sparsity()) 825 for sB in [ Sparsity.dense(a.size1(),1), vertcat(Sparsity.dense(1,1),Sparsity(a.size1()-1,1)),Sparsity.lower(a.size1()),Sparsity.lower(a.size1()).T]: 826 827 b = DM(sB,[random.random() for i in range(sB.nnz())]) 828 B = SX.sym("B",b.sparsity()) 829 C = solve(A,B) 830 831 f = Function("f", [A,B],[C]) 832 833 834 c = f(a,b) 835 836 c_ref = DM(linalg.solve(a,b)) 837 c_ref = sparsify(c_ref) 838 839 print(sA.dim(), sB.dim()) 840 841 842 843 #try: 844 print("foo",c,c_ref) 845 self.checkarray(c,c_ref) 846 #self.assertTrue(min((IM.ones(c_ref.sparsity())-IM.ones(c.sparsity())).nonzeros())==0) 847 #except Exception as e: 848 # c.print_dense() 849 # print("sol:") 850 # c.sparsity().spy() 851 # print("ref:") 852 # c_ref.sparsity().spy() 853 # c_ref.print_dense() 854 # a.print_dense() 855 # raise e 856 857 def test_kron(self): 858 a = sparsify(DM([[1,0,6],[2,7,0]])) 859 b = sparsify(DM([[1,0,0],[2,3,7],[0,0,9],[1,12,13]])) 860 861 c_ = c.kron(a,b) 862 863 self.assertEqual(c_.size1(),a.size1()*b.size1()) 864 self.assertEqual(c_.size2(),a.size2()*b.size2()) 865 self.assertEqual(c_.nnz(),a.nnz()*b.nnz()) 866 867 self.checkarray(c_,numpy.kron(a,b)) 868 869 def test_vec_kron(self): 870 A = SX.sym("A",2,3) 871 B = SX.sym("B",4,5) 872 P = SX.sym("P",A.size2(),B.size1()) 873 874 f = Function("f", [vec(P.T),A,B],[vec(mtimes([A,P,B]).T)]) 875 876 J = f.jacobian_old(0, 0) 877 J_in = [] 878 J_in.append(numpy.random.rand(*vec(P.T).shape)) 879 J_in.append(numpy.random.rand(*A.shape)) 880 J_in.append(numpy.random.rand(*B.shape)) 881 882 res, _ = J(*J_in) 883 884 ref = kron(J_in[1],J_in[2].T) 885 886 self.checkarray(res,ref) 887 888 def test_repmat(self): 889 a = DM([[1,2],[3,4],[5,6]]) 890 self.checkarray(repmat(a,2,3),kron(DM.ones(2,3),a)) 891 892 def test_triu(self): 893 a = DM([[1,2],[3,4]]) 894 b = triu(a) 895 self.checkarray(b, DM([[1,2],[0,4]]) ) 896 897 898 def test_tril(self): 899 a = DM([[1,2],[3,4]]) 900 b = tril(a) 901 self.checkarray(b, DM([[1,0],[3,4]]) ) 902 903 def test_nz(self): 904 a = sparsify(DM([[1,2],[0,0],[3,4]])) 905 self.checkarray(a.nz[:], DM([1,3,2,4]) ) 906 self.checkarray(len(a.nz), 4 ) 907 self.checkarray(a.nz[:-1], DM([1,3,2]) ) 908 self.checkarray(a.nz[0], DM([1]) ) 909 910 def test_norm_inf_mul(self): 911 numpy.random.seed(0) 912 913 A = numpy.random.random((10,2)) 914 B = numpy.random.random((2,8)) 915 916 self.checkarray(norm_inf_mul(A,B),norm_inf(mtimes(A,B))) 917 self.checkarray(DM(norm_0_mul(A,B)),mtimes(A,B).nnz()) 918 919 # Sparse 920 for i in range(5): 921 A[numpy.random.randint(A.shape[0]),numpy.random.randint(A.shape[1])] = 0 922 B[numpy.random.randint(B.shape[0]),numpy.random.randint(B.shape[1])] = 0 923 924 A = sparsify(A) 925 B = sparsify(B) 926 927 self.checkarray(norm_inf_mul(A,B),norm_inf(mtimes(A,B))) 928 self.checkarray(DM(norm_0_mul(A,B)),mtimes(A,B).nnz()) 929 930 931 A = numpy.random.random((8,2)) 932 B = numpy.random.random((2,10)) 933 934 self.checkarray(norm_inf_mul(A,B),norm_inf(mtimes(A,B))) 935 self.checkarray(DM(norm_0_mul(A,B)),mtimes(A,B).nnz()) 936 937 # Sparse 938 for i in range(5): 939 A[numpy.random.randint(A.shape[0]),numpy.random.randint(A.shape[1])] = 0 940 B[numpy.random.randint(B.shape[0]),numpy.random.randint(B.shape[1])] = 0 941 942 A = sparsify(A) 943 B = sparsify(B) 944 945 self.checkarray(norm_inf_mul(A,B),norm_inf(mtimes(A,B))) 946 self.checkarray(DM(norm_0_mul(A,B)),mtimes(A,B).nnz()) 947 948 def test_mul3_issue_1465(self): 949 with self.assertRaises(Exception): 950 w = SX.sym("w",2,1) 951 Q = np.eye(2) 952 mtimes(w.T,Q,w) 953 954 def test_chol(self): 955 numpy.random.seed(0) 956 957 for i in range(4): 958 A = numpy.random.random((3,3)) 959 H = mtimes(A,A.T) 960 961 R = chol(H) 962 963 assert R.is_triu() 964 self.checkarray(mtimes(R.T,R),H) 965 def test_skew(self): 966 x = DM([1,7,13]) 967 self.checkarray(inv_skew(skew(x)),x) 968 y = DM([0.2,0.9,0.4]) 969 self.checkarray(mtimes(skew(x),y),cross(x,y)) 970 971 def test_nz_overflow(self): 972 d = DM([2,3]) 973 r = d.nz[:] 974 self.checkarray(r,d) 975 976 def test_DMcrash(self): 977 with self.assertRaises(Exception): 978 DM([DM([1,2]),DM([1,2])]) 979 a = DM([DM([1]),DM([2])]) 980 self.checkarray(a,DM([1,2])) 981 982 def test_sparsity_operation(self): 983 L = [DM(1), DM(Sparsity(1,1),1), DM(Sparsity(2,1),1), DM(Sparsity.dense(2,1),1)] 984 985 for a in L: 986 for b in L: 987 c = a*b 988 989 if a.nnz()==0 or b.nnz()==0: 990 self.assertTrue(c.nnz()==0) 991 else: 992 self.assertTrue(c.nnz()>0) 993 994 self.assertTrue(sum2(DM(Sparsity(1,1),1)).nnz()==0) 995 996 def test_matlab_operations(self): 997 998 data = [ np.array([[1,3],[11,17]]) , np.array([[1,3]]) ,np.array([[1],[3]]), np.array([[3]])] 999 1000 for A in data: 1001 B = reshape(DM(A),A.shape) 1002 #self.checkarray(np.cumsum(A),cumsum(B)) 1003 self.checkarray(np.cumsum(A,0),cumsum(B,0)) 1004 self.checkarray(np.cumsum(A,1),cumsum(B,1)) 1005 1006 Bs = MX.sym("B",B.shape) 1007 A0 = cumsum(B,0) 1008 A1 = cumsum(B,1) 1009 1010 self.checkarray(np.cumsum(A,0),evalf(cumsum(MX(B),0))) 1011 self.checkarray(np.cumsum(A,1),evalf(cumsum(B,1))) 1012 1013 f = Function("f",[Bs],[A0,A1]) 1014 self.checkfunction(f,f.expand(),inputs = [B]) 1015 self.check_codegen(f,inputs=[B],check_serialize=True) 1016 1017 #self.checkarray(np.diff(A),diff(B)) 1018 1019 #self.checkarray(np.diff(A,1),diff(B,1)) 1020 #self.checkarray(np.diff(A,2),diff(B,2)) 1021 1022 self.checkarray(np.diff(A,1,0),diff(B,1,0)) 1023 #self.checkarray(np.diff(A,1,1),diff(B,1,1)) 1024 1025 1026 def test_singular_repmat(self): 1027 for X in [DM, SX, MX, Sparsity]: 1028 for n_b in [0,2]: 1029 for m_b in [0,2]: 1030 b = X(n_b, m_b) 1031 1032 for n in [0,3]: 1033 for m in [0,3]: 1034 self.assertEqual(repmat(b, n, m).shape,(n_b * n, m_b * m)) 1035 1036 1037 def test_serialize(self): 1038 for a in [DM(), DM(2), DM([1,2]),DM([[1,2],[3,4],[5,6]])]: 1039 b = DM.deserialize(a.serialize()) 1040 self.checkarray(a,b) 1041 1042 def test_iterable(self): 1043 a = DM([1,2,3]) 1044 b = list(iter(a.nz)) 1045 self.checkarray(a,DM(b)) 1046 1047 with self.assertInException("CasADi matrices are not iterable"): 1048 iter(a) 1049 1050 def test_from_file(self): 1051 a = sparsify(DM([[1,0,-6,5,0],[4,0,-4e-301,9.3e-18,0]])) 1052 1053 a.print_dense() 1054 1055 a[0,4] = 0 1056 a[1,1] = 0 1057 a.to_file("test.mtx") 1058 1059 with self.assertInException("foobar.mtx"): 1060 DM.from_file("foobar.mtx") 1061 b = DM.from_file("test.mtx") 1062 1063 self.assertTrue(a.sparsity()==b.sparsity()) 1064 self.checkarray(a,b) 1065 1066 a.to_file("test.txt") 1067 b = DM.from_file("test.txt") 1068 self.checkarray(a,b) 1069 1070 for s in ["1 00 -6 5 0\n4 0 -4e-301 9.3e-18 00\n", 1071 "1 00 -6 5 0\n4 0 -4e-301 9.3e-18 00", 1072 "1 00 -6 5 0 \n4 0 -4e-301 9.3e-18 00 ", 1073 "1\t00\t-6\t5\t0\t\n4 0 -4e-301 9.3e-18 00\n", 1074 "1 00 -6 5 0\n%foobar\n4 0 -4e-301 9.3e-18 00\n"]: 1075 1076 with open("test2.txt","w") as f: 1077 f.write(s) 1078 b = DM.from_file("test2.txt") 1079 self.assertTrue(a.sparsity()==b.sparsity()) 1080 self.assertTrue(float(norm_1(a-b))==0) 1081 1082 with open("test.txt","w") as f: 1083 f.write("1 aa 6\n4 8 9e-18\n") 1084 1085 1086 with self.assertInException("Parsing"): 1087 DM.from_file("test.txt") 1088 1089 a = DM([3, inf, -inf, np.nan]) 1090 a.to_file("test.txt") 1091 b = DM.from_file("test.txt") 1092 self.checkarray(a,b) 1093 def test_norm_2(self): 1094 a = np.array([1,2,3]) 1095 r = np.linalg.norm(a) 1096 for M in [DM,MX,SX]: 1097 self.checkarray(evalf(norm_2(M(a))),r) 1098 self.checkarray(evalf(norm_2(M(a).T)),r) 1099 1100 def test_ldl(self): 1101 H = diagcat(DM.rand(5,5),DM.rand(5,5)) 1102 H = H+H.T+2*DM.eye(10) 1103 1104 p = np.random.permutation(10) 1105 H = H[p,p] 1106 1107 1108 [D,Lt,p] = ldl(H) 1109 P = DM.eye(10)[:,p] 1110 1111 F = mtimes(mtimes(sqrt(diag(D)),DM.eye(10)+Lt),P.T) 1112 print(H-mtimes(F.T,F)) 1113 self.assertTrue(norm_fro(H-mtimes(F.T,F))<=1e-14) 1114 1115 1116 def test_im_bugs(self): 1117 a = vertcat(1,2) 1118 self.assertTrue(isinstance(a,DM)) 1119 self.checkarray(c.linspace(1,3,10),c.linspace(1.0,3.0,10)) 1120 1121if __name__ == '__main__': 1122 unittest.main() 1123