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 random, array, linalg, matrix, zeros, ones, ndarray, eye 28import unittest 29from types import * 30from helpers import * 31from copy import deepcopy 32 33import sys 34 35if sys.version_info >= (3, 0): 36 TupleType = tuple 37 38 39warnings.filterwarnings("ignore",category=DeprecationWarning) 40 41scipy_available = True 42try: 43 from scipy.sparse import csr_matrix 44except: 45 scipy_available = False 46 47def checkarray(self,zr,zt,name): 48 if len(zr.shape)==1 and (zt.shape[0]==1 or zt.shape[1]==1) and zr.shape[0]==zt.shape[1]*zt.shape[0]: 49 zr=reshape(zr,(zt.shape)); 50 self.assertEqual(zt.shape[0],zr.shape[0],"%s dimension error. Got %s, expected %s" % (name,str(zt.shape),str(zr.shape))) 51 self.assertEqual(len(zt.shape),len(zr.shape),"%s dimension error. Got %s, expected %s" % (name,str(zt.shape),str(zr.shape))) 52 self.assertEqual(zt.shape[1],zr.shape[1],"%s dimension error. Got %s, expected %s" % (name,str(zt.shape),str(zr.shape))) 53 for i in range(zr.shape[0]): 54 for j in range(zr.shape[1]): 55 self.assertAlmostEqual(zt[i,j],zr[i,j],10,"%s evaluation error. %s <-> %s" % (name, str(zt),str(zr))) 56 57def checkMXoperations(self,ztf,zrf,name): 58 x = MX.sym("x",1,3) 59 z=vertcat(*[x*(i+1) for i in range(8)]) 60 f = Function("f", [x],[ztf(z)]) 61 L=[1,2,3] 62 f_out = f(L) 63 zt = f_out.full() 64 zr = array([[L[0]*(i+1),L[1]*(i+1),L[2]*(i+1)] for i in range(8)]) 65 checkarray(self,zrf(zr),zt,name) 66 return (zt,zrf(zr)) 67 68def checkMXoperations2(self,ztf,zrf,name): 69 x = MX.sym("x",3,1) 70 z = horzcat(*[x*i for i in range(8)]) 71 f = Function("f", [x],[ztf(z)]) 72 L=[1,2,3] 73 f_out = f(L) 74 zt = f_out.full() 75 zr = array([[L[0]*i,L[1]*i,L[2]*i] for i in range(8)]).T 76 checkarray(self,zrf(zr),zt,name) 77 return zt 78 79def checkMXoperations3(self,ztf,zrf,name): 80 x = MX.sym("x",3,1) 81 p = horzcat(*[x[0,0],x[1,0],x[2,0]]) 82 z = vertcat(*[p*i for i in range(8)]) 83 f = Function("f", [x],[ztf(z)]) 84 L=[1,2,3] 85 f_out = f(L) 86 zt = f_out.full() 87 zr = array([[L[0]*i,L[1]*i,L[2]*i] for i in range(8)]) 88 checkarray(self,zrf(zr),zt,name) 89 return (zt,zrf(zr)) 90 91class MXtests(casadiTestCase): 92 93 def setUp(self): 94 self.pool=FunctionPool() 95 self.pool.append(lambda x: sqrt(x[0]),sqrt,"sqrt") 96 self.pool.append(lambda x: sin(x[0]),sin,"sin") 97 self.pool.append(lambda x: cos(x[0]),cos,"cos") 98 self.pool.append(lambda x: tan(x[0]),tan,"tan") 99 self.pool.append(lambda x: arctan(x[0]),arctan,"arctan") 100 self.pool.append(lambda x: arcsin(x[0]),arcsin,"arcsin") 101 self.pool.append(lambda x: arccos(x[0]),arccos,"arccos") 102 self.pool.append(lambda x: exp(x[0]),exp,"exp") 103 self.pool.append(lambda x: log(x[0]),log,"log",flags={'nozero'}) 104 self.pool.append(lambda x: x[0]**0,lambda x : x**0,"x^0",flags={'nozero'}) 105 self.pool.append(lambda x: x[0]**1,lambda x : x**1,"^1") 106 self.pool.append(lambda x: x[0]**(-2),lambda x : x**(-2),"^-2",flags={'nozero'}) 107 self.pool.append(lambda x: x[0]**(0.3),lambda x : x**(0.3),"^0.3") 108 self.pool.append(lambda x: floor(x[0]),floor,"floor") 109 self.pool.append(lambda x: ceil(x[0]),ceil,"ceil") 110 self.Jpool=FunctionPool() 111 self.Jpool.append(lambda x: sqrt(x[0]),lambda x:diag(1/(2.0*sqrt(x))),"sqrt") 112 self.Jpool.append(lambda x: sin(x[0]),lambda x:diag(cos(x)),"sin") 113 self.Jpool.append(lambda x: cos(x[0]),lambda x:diag(-sin(x)),"cos") 114 self.Jpool.append(lambda x: tan(x[0]),lambda x:diag(1.0/cos(x)**2),"tan") 115 self.Jpool.append(lambda x: arctan(x[0]),lambda x:diag( 1.0/(x**2+1)),"arctan") 116 self.Jpool.append(lambda x: arcsin(x[0]),lambda x:diag( 1.0/sqrt(1-x**2)),"arcsin") 117 self.Jpool.append(lambda x: arccos(x[0]),lambda x: diag(-1.0/sqrt(1-x**2)),"arccos") 118 self.Jpool.append(lambda x: exp(x[0]),lambda x: diag(exp(x)),"exp") 119 self.Jpool.append(lambda x: log(x[0]),lambda x: diag(1.0/x),"log") 120 self.Jpool.append(lambda x: x[0]**0,lambda x :diag(zeros(x.shape)),"x^0") 121 self.Jpool.append(lambda x: x[0]**1,lambda x : diag(ones(x.shape)),"^1") 122 self.Jpool.append(lambda x: x[0]**(-2),lambda x : diag(-2.0/x**3),"^-2") 123 self.Jpool.append(lambda x: x[0]**(0.3),lambda x :diag( 0.3/x**0.7),"^0.3") 124 self.matrixpool=FunctionPool() 125 #self.matrixpool.append(lambda x: norm_2(x[0]),linalg.norm,"norm_2") 126 #self.matrixpool.append(lambda x: norm_1(x[0]),lambda x: sum(sum(abs(x))),"norm_1") 127 #self.matrixpool.append(lambda x: norm_inf(x[0]),lambda x: abs(matrix(x)).max(),"norm_inf") 128 self.matrixbinarypool=FunctionPool() 129 self.matrixbinarypool.append(lambda a: a[0]+a[1],lambda a: a[0]+a[1],"Matrix+Matrix") 130 self.matrixbinarypool.append(lambda a: a[0]-a[1],lambda a: a[0]-a[1],"Matrix-Matrix") 131 self.matrixbinarypool.append(lambda a: a[0]*a[1],lambda a: a[0]*a[1],"Matrix*Matrix") 132 self.matrixbinarypool.append(lambda a: fmax(a[0],a[1]),lambda a: fmax(a[0],a[1]),"fmin") 133 134 self.matrixbinarypool.append(lambda a: fmin(a[0],a[1]),lambda a: fmin(a[0],a[1]),"fmax") 135 self.matrixbinarypool.append(lambda a: mtimes(a[0],a[1].T),lambda a: numpy.dot(a[0],a[1].T),"mtimes(Matrix,Matrix.T)") 136 self.matrixbinarypool.append(lambda a: arctan2(a[0],a[1]),lambda a: arctan2(a[0],a[1]),"arctan2") 137 #self.matrixbinarypool.append(lambda a: inner_mul(a[0],trans(a[1])),lambda a: c.dot(a[0].T,a[1]),name="inner_mul(Matrix,Matrix)") 138 self.matrixbinarypool.append(lambda a: mtimes(a[0],a[1].T),lambda a: numpy.dot(a[0],a[1].T),"mtimes(Matrix,Matrix.T)") 139 140 def test_MX1(self): 141 self.message("MX constructor") 142 x = MX.sym("x",2,3) 143 self.assertEqual(x.size1(),2,"MX fails to indicate its size1") 144 self.assertEqual(x.size2(),3,"MX fails to indicate its size2") 145 146 def test_MXvertcat(self): 147 self.message("MX vertcat") 148 x = MX.sym("x",1,3) 149 y = MX.sym("y",1,3) 150 z=vertcat(*(x,y)) 151 self.assertEqual(z.size1(),2,"MX fails to indicate its size1") 152 self.assertEqual(z.size2(),3,"MX fails to indicate its size2") 153 154 def test_MX_fun1(self): 155 self.message("MXFunction single input, single output") 156 # check if x->2*x 157 # evaluates correctly for x=3 158 x = MX.sym("x") 159 y = 2*x 160 f = Function("f", [x],[y]) 161 self.assertEqual(f.n_in(),1,"Function fails to indicate correct number of inputs") 162 self.assertEqual(f.n_out(),1,"Function fails to indicate correct number of outputs") 163 f_out = f(3) 164 yt = tuple(f_out.nonzeros()) 165 self.assertEqual(type(yt),TupleType,"Output of Function is expected to be tuple of floats") 166 self.assertEqual(len(yt),1,"Output of Function was tuple of floats, as expected, but length is incorrect.") 167 y=yt[0] 168 self.assertEqual(type(y),float,"Output of Function is expected to be tuple of floats") 169 self.assertAlmostEqual(y, 2*3,10) 170 171 def test_MXfunction2(self): 172 self.message("Function multi input, multi output") 173 # check if [x,y]->[y+x,y*x] 174 # evaluates correctly for x=3,y=7 175 x = MX.sym("x") 176 y = MX.sym("y") 177 f = Function("f", [x,y],[x+y,y*x]) 178 self.assertEqual(f.n_in(),2,"Function fails to indicate correct number of inputs") 179 self.assertEqual(f.n_out(),2,"Function fails to indicate correct number of outputs") 180 f_out = f(3, 7) 181 zt1 = tuple(f_out[0].nonzeros()) 182 zt2 = tuple(f_out[1].nonzeros()) 183 self.assertEqual(type(zt1),TupleType,"Output of Function is expected to be tuple of floats") 184 self.assertEqual(type(zt2),TupleType,"Output of Function is expected to be tuple of floats") 185 self.assertEqual(len(zt1),1,"Output of Function was tuple of floats, as expected, but length is incorrect.") 186 self.assertEqual(len(zt2),1,"Output of Function was tuple of floats, as expected, but length is incorrect.") 187 z1=zt1[0] 188 z2=zt2[0] 189 self.assertEqual(type(z1),float,"Output of Function is expected to be tuple of floats") 190 self.assertEqual(type(z2),float,"Output of Function is expected to be tuple of floats") 191 self.assertAlmostEqual(z2, 21,10) 192 self.assertAlmostEqual(z1, 10,10) 193 194 195 196 def test_MXfunction3(self): 197 self.message("Function single input, multi output (1)") 198 # check if [x,y]->[y+x,y*x] 199 # evaluates correctly for x=3,y=7 200 # now with single input, multi output 201 xy = MX.sym("xy",2) 202 f = Function("f", [xy],[xy[0]+xy[1],xy[0]*xy[1]]) 203 self.assertEqual(f.n_in(),1,"Function fails to indicate correct number of inputs") 204 self.assertEqual(f.n_out(),2,"Function fails to indicate correct number of outputs") 205 f_out = f([3,7]) 206 zt1 = tuple(f_out[0].nonzeros()) 207 zt2 = tuple(f_out[1].nonzeros()) 208 self.assertEqual(type(zt1),TupleType,"Output of Function is expected to be tuple of floats") 209 self.assertEqual(type(zt2),TupleType,"Output of Function is expected to be tuple of floats") 210 self.assertEqual(len(zt1),1,"Output of Function was tuple of floats, as expected, but length is incorrect.") 211 self.assertEqual(len(zt2),1,"Output of Function was tuple of floats, as expected, but length is incorrect.") 212 z1=zt1[0] 213 z2=zt2[0] 214 self.assertEqual(type(z1),float,"Output of Function is expected to be tuple of floats") 215 self.assertEqual(type(z2),float,"Output of Function is expected to be tuple of floats") 216 self.assertAlmostEqual(z2, 21,10) 217 self.assertAlmostEqual(z1, 10,10) 218 219 def test_MXfunction3b(self): 220 self.message("Function single input, multi output (2)") 221 # check if [x,y]->[y+x,y*x] 222 # evaluates correctly for x=3,y=7 223 # now with single input, multi output 224 xy = MX.sym("xy",1,2) 225 f = Function("f", [xy],[xy[0,0]+xy[0,1],xy[0,0]*xy[0,1]]) 226 self.assertEqual(f.n_in(),1,"Function fails to indicate correct number of inputs") 227 self.assertEqual(f.n_out(),2,"Function fails to indicate correct number of outputs") 228 f_out = f([3,7]) 229 zt1 = f_out[0].full() 230 zt2 = f_out[1].full() 231 232 self.assertEqual(type(zt1),ndarray,"Output of Function is expected to be numpy.ndarray") 233 self.assertEqual(zt1.shape[0],1,"Output of Function is of wrong shape.") 234 self.assertEqual(zt1.shape[1],1,"Output of Function is of wrong shape.") 235 236 self.assertEqual(type(zt2),ndarray,"Output of Function is expected to be numpy.ndarray") 237 self.assertEqual(zt2.shape[0],1,"Output of Function is of wrong shape.") 238 self.assertEqual(zt2.shape[1],1,"Output of Function is of wrong shape.") 239 240 z1=zt1[0,0] 241 z2=zt2[0,0] 242 self.assertEqual(type(z1),numpy.float64,"Output of Function is expected to be numpy.ndarray of floats") 243 self.assertEqual(type(z2),numpy.float64,"Output of Function is expected to be numpy.ndarray of floats") 244 self.assertAlmostEqual(z2, 21,10) 245 self.assertAlmostEqual(z1, 10,10) 246 247 def test_MXfunction4(self): 248 self.message("Function single input, single output , using vertcat") 249 # check if [x,y]->[y+x,y*x] 250 # evaluates correctly for x=3,y=7 251 # now with single input, single output 252 xy = MX.sym("xy",2) 253 z=vertcat(*[xy[0]+xy[1],xy[0]*xy[1]]) 254 f = Function("f", [xy],[z]) 255 self.assertEqual(f.n_in(),1,"Function fails to indicate correct number of inputs") 256 self.assertEqual(f.n_out(),1,"Function fails to indicate correct number of outputs") 257 f_out = f([3,7]) 258 zt=f_out.full() 259 self.assertEqual(type(zt),ndarray,"Output of Function is expected to be numpy.ndarray") 260 self.assertEqual(zt.shape[0],2,"Output of Function is of wrong shape.") 261 self.assertEqual(zt.shape[1],1,"Output of Function is of wrong shape.") 262 z1=zt[0,0] 263 z2=zt[1,0] 264 self.assertEqual(type(z1),numpy.float64,"Output of Function is expected to be numpy.ndarray of floats") 265 self.assertEqual(type(z2),numpy.float64,"Output of Function is expected to be numpy.ndarray of floats") 266 self.assertAlmostEqual(z2, 21,10) 267 self.assertAlmostEqual(z1, 10,10) 268 269 def test_MXfunction5(self): 270 self.message("Function single input, single output , using horzcat") 271 # check if [x,y]->[y+x,y*x] 272 # evaluates correctly for x=3,y=7 273 # now with single input, single output 274 xy = MX.sym("xy",2) 275 z=horzcat(*[xy[0]+xy[1],xy[0]*xy[1]]) 276 f = Function("f", [xy],[z]) 277 self.assertEqual(f.n_in(),1,"Function fails to indicate correct number of inputs") 278 self.assertEqual(f.n_out(),1,"Function fails to indicate correct number of outputs") 279 f_out = f([3,7]) 280 zt = f_out.full() 281 self.assertEqual(type(zt),ndarray,"Output of Function is expected to be numpy.ndarray") 282 self.assertEqual(zt.shape[0],1,"Output of Function is of wrong shape.") 283 self.assertEqual(zt.shape[1],2,"Output of Function is of wrong shape.") 284 z1=zt[0,0] 285 z2=zt[0,1] 286 self.assertEqual(type(z1),numpy.float64,"Output of Function is expected to be numpy.ndarray of floats") 287 self.assertEqual(type(z2),numpy.float64,"Output of Function is expected to be numpy.ndarray of floats") 288 self.assertAlmostEqual(z2, 21,10) 289 self.assertAlmostEqual(z1, 10,10) 290 291 def test_which_depends_empty(self): 292 for X in [SX,MX]: 293 x=X.sym("x") 294 295 for tr in [True,False]: 296 for i in [0,1,2]: 297 self.assertEqual(which_depends(x,X(0,1),i,tr),[False]*(1 if tr else 0)) 298 299 self.assertEqual(which_depends(X(0,1),x,i,tr),[False]*(0 if tr else 1)) 300 self.assertTrue(len(which_depends(X(0,1),X(0,1),i,tr))==0) 301 302 def test_issue83(self): 303 x=MX.sym("x") 304 y=MX.sym("y") 305 306 z = x + y 307 308 f = Function("f", [x,y],[z]) 309 310 fc = f(MX(3),y) 311 312 g = Function("g", [y],[fc]) 313 g_in = [7] 314 g_out = g(g_in) 315 316 self.assertAlmostEqual(g_out[0],10,10,"issue #83") 317 318 fc = f(x,MX(7)) 319 320 g = Function("g", [x],[fc]) 321 g_in = [3] 322 g_out = g(g_in) 323 324 self.assertAlmostEqual(g_out[0],10,10,"issue #83") 325 326 def test_identitySX(self): 327 self.message("identity SXFunction") 328 x = SX.sym("x") 329 f = Function("f", [x],[x]) 330 f_in = [3] 331 f_out = f(f_in) 332 self.assertAlmostEqual(f_out[0,0], 3,10) 333 334 def test_identityMX(self): 335 self.message("identity Function") 336 x = MX.sym("x") 337 f = Function("f", [x],[x]) 338 f_in = [3] 339 f_out = f(f_in) 340 self.assertAlmostEqual(f_out[0,0], 3,10) 341 342 def test_MXorder(self): 343 self.message("Function order of non-zero elements") 344 x = MX.sym("x",2,3) 345 f = Function("f", [x],[x+x]) 346 347 self.assertEqual(f.n_in(),1,"Function fails to indicate correct number of inputs") 348 self.assertEqual(f.n_out(),1,"Function fails to indicate correct number of outputs") 349 L=[1,2,3,4,5,6] 350 f_in = DM(f.sparsity_in(0),L) 351 f_out = f(f_in) 352 zt = f_out.full() 353 self.assertEqual(zt.shape[0],2,"Output of Function is of wrong shape.") 354 self.assertEqual(zt.shape[1],3,"Output of Function is of wrong shape.") 355 356 Lr=numpy.reshape(L,(2,3),'F') 357 for i in range(2): 358 for j in range(3): 359 self.assertAlmostEqual(Lr[i,j]*2, zt[i,j],10) 360 361 def test_trans(self): 362 self.message("trans") 363 a = MX(0,1) 364 b = a.T 365 self.assertEqual(b.size1(),1) 366 self.assertEqual(b.size2(),0) 367 368 def test_MXtrans(self): 369 self.message("trans(MX)") 370 x = MX.sym("x",2,3) 371 z=x.T 372 self.assertEqual(z.size1(),3,"Vec returns MX of wrong dimension") 373 self.assertEqual(z.size2(),2,"Vec returns MX of wrong dimension") 374 f = Function("f", [x],[z]) 375 self.assertEqual(f.n_in(),1,"Function fails to indicate correct number of inputs") 376 self.assertEqual(f.n_out(),1,"Function fails to indicate correct number of outputs") 377 L=[1,2,3,4,5,6] 378 f_in = DM(f.sparsity_in(0),L) 379 f_out = f(f_in) 380 zt = f_out.full() 381 382 ztr=numpy.reshape(zt,(3,2)) 383 Lr=numpy.reshape(L,(2,3),'F') 384 for i in range(2): 385 for j in range(3): 386 self.assertAlmostEqual(Lr[i,j], ztr[j,i],10) 387 388 def test_MXvec(self): 389 390 u = DM([[10*j+i for i in range(3)] for j in range(4) ]) 391 392 U = MX.sym("u",u.shape) 393 394 f = Function("f", [U],[vec(U)]) 395 f_out = f(u) 396 397 self.checkarray(vec(u),f_out,"vec") 398 399 def test_MXreshape(self): 400 self.message("reshape(MX)") 401 x = MX.sym("x",2,3) 402 z=c.reshape(x,(1,6)) 403 self.assertEqual(z.size1(),1,"Vec returns MX of wrong dimension") 404 self.assertEqual(z.size2(),6,"Vec returns MX of wrong dimension") 405 f = Function("f", [x],[z]) 406 self.assertEqual(f.n_in(),1,"Function fails to indicate correct number of inputs") 407 self.assertEqual(f.n_out(),1,"Function fails to indicate correct number of outputs") 408 L=[1,2,3,4,5,6] 409 f_in = DM(f.sparsity_in(0),L) 410 f_out = f(f_in) 411 zt = f_out.full() 412 for i in range(len(L)): 413 self.assertAlmostEqual(L[i], zt[0,i],10) 414 415 def test_MXcompose(self): 416 self.message("compositions of vec, trans, reshape with vertcat") 417 checkMXoperations(self,lambda x: x,lambda x: x,'vertcat') 418 checkMXoperations(self,lambda x: x.T,lambda x: x.T,'trans(vertcat)') 419 checkMXoperations(self,lambda x: x.T.T,lambda x: x,'trans(trans(vertcat))') 420 checkMXoperations(self,lambda x: vec(x.T),lambda x: numpy.reshape(x,(numpy.prod(x.shape),1)),'vec(trans(vertcat))') 421 checkMXoperations(self,lambda x: vec(x).T,lambda x: numpy.reshape(x.T,(numpy.prod(x.shape),1)).T,'vec(trans(vertcat))') 422 checkMXoperations(self,lambda x: c.reshape(x.T,(6,4)).T,lambda x: numpy.reshape(x,(4,6)),'reshape(vertcat)') 423 checkMXoperations(self,lambda x: c.reshape(x,(6,4)).T,lambda x: numpy.reshape(x.T,(4,6)),'reshape(trans(vertcat))') 424 checkMXoperations(self,lambda x: c.reshape(x.T,(6,4)),lambda x: numpy.reshape(x,(4,6)).T,'trans(reshape(vertcat))') 425 426 def test_MXcompose2(self): 427 self.message("compositions of vec, trans, reshape with horzcat") 428 checkMXoperations2(self,lambda x: x,lambda x: x,'horzcat') 429 checkMXoperations2(self,lambda x: x.T,lambda x: x.T,'trans(horzcat)') 430 checkMXoperations2(self,lambda x: x.T.T,lambda x: x,'trans(trans(horzcat))') 431 checkMXoperations2(self,lambda x: vec(x.T),lambda x: numpy.reshape(x,(numpy.prod(x.shape),1)),'vec(trans(horzcat))') 432 checkMXoperations2(self,lambda x: vec(x).T,lambda x: numpy.reshape(x.T,(numpy.prod(x.shape),1)).T,'vec(trans(horzcat))') 433 checkMXoperations2(self,lambda x: c.reshape(x.T,(6,4)).T,lambda x: numpy.reshape(x,(4,6)),'reshape(horzcat)') 434 checkMXoperations2(self,lambda x: c.reshape(x,(6,4)).T,lambda x: numpy.reshape(x.T,(4,6)),'reshape(trans(horzcat))') 435 checkMXoperations2(self,lambda x: c.reshape(x.T,(6,4)),lambda x: numpy.reshape(x,(4,6)).T,'trans(reshape(horzcat))') 436 437 def test_MXcompose3(self): 438 self.message("compositions of vec, trans, reshape with vertcat") 439 checkMXoperations3(self,lambda x: x,lambda x: x,'snippet') 440 checkMXoperations3(self,lambda x: x.T,lambda x: x.T,'trans(snippet)') 441 checkMXoperations3(self,lambda x: x.T.T,lambda x: x,'trans(trans(snippet))') 442 checkMXoperations3(self,lambda x: vec(x.T),lambda x: numpy.reshape(x,(numpy.prod(x.shape),1)),'vec(trans(snippet))') 443 checkMXoperations3(self,lambda x: vec(x).T,lambda x: numpy.reshape(x.T,(numpy.prod(x.shape),1)).T,'vec(trans(snippet))') 444 checkMXoperations3(self,lambda x: c.reshape(x.T,(6,4)).T,lambda x: numpy.reshape(x,(4,6)),'reshape(snippet)') 445 checkMXoperations3(self,lambda x: c.reshape(x,(6,4)).T,lambda x: numpy.reshape(x.T,(4,6)),'reshape(trans(snippet))') 446 checkMXoperations3(self,lambda x: c.reshape(x.T,(6,4)),lambda x: numpy.reshape(x,(4,6)).T,'trans(reshape(snippet))') 447 448 def test_MXcompose4(self): 449 self.message("compositions of horzcat + vertcat") 450 checkMXoperations(self,lambda x: vertcat(*[x]),lambda x: x,'vertcat(*vertcat)') 451 checkMXoperations(self,lambda x: vertcat(*[x,x*2]),lambda x: numpy.vstack((x,x*2)),'vertcat(*vertcat,vertcat)') 452 checkMXoperations(self,lambda x: horzcat(*[x]),lambda x: x,'horzcat(*vertcat)') 453 checkMXoperations(self,lambda x: horzcat(*[x,x*2]),lambda x: numpy.hstack((x,x*2)),'horzcat(*vertcat,vertcat)') 454 455 checkMXoperations2(self,lambda x: vertcat(*[x]),lambda x: x,'vertcat(*horzcat)') 456 checkMXoperations2(self,lambda x: vertcat(*[x,x*2]),lambda x: numpy.vstack((x,x*2)),'vertcat(*horzcat,horzcat)') 457 checkMXoperations2(self,lambda x: horzcat(*[x]),lambda x: x,'horzcat(*horzcat)') 458 checkMXoperations2(self,lambda x: horzcat(*[x,x*2]),lambda x: numpy.hstack((x,x*2)),'horzcat(*horzcat,horzcat)') 459 460 checkMXoperations3(self,lambda x: vertcat(*[x]),lambda x: x,'vertcat(*snippet)') 461 checkMXoperations3(self,lambda x: vertcat(*[x,x*2]),lambda x: numpy.vstack((x,x*2)),'vertcat(*snippet,snippet)') 462 checkMXoperations3(self,lambda x: horzcat(*[x]),lambda x: x,'horzcat(*snippet)') 463 checkMXoperations3(self,lambda x: horzcat(*[x,x*2]),lambda x: numpy.hstack((x,x*2)),'horzcat(*snippet,snippet)') 464 465 @known_bug() # Test refactoring, cf. #1436 466 def test_MXslicingnew(self): 467 self.message("MX slicing new") 468 469 self.message(":dense") 470 x = MX.sym("x",3,2) 471 x0=array([[0.738,0.2],[ 0.1,0.39 ],[0.99,0.999999]]) 472 self.numpyEvaluationCheck(lambda x:x[0][1,0],lambda x: x[1,0],[x],x0,'x[1,0]') 473 self.numpyEvaluationCheck(lambda x:x[0][0,1],lambda x: x[0,1],[x],x0,'x[0,1]') 474 self.numpyEvaluationCheck(lambda x:x[0][0,-1],lambda x: x[0,-1],[x],x0,'x[0,-1]') 475 self.numpyEvaluationCheck(lambda x: x[0][:,0], lambda x: matrix(x)[:,0],[x],x0,name="x[:,0]") 476 self.numpyEvaluationCheck(lambda x: x[0][:,1], lambda x: matrix(x)[:,1],[x],x0,name="x[:,1]") 477 self.numpyEvaluationCheck(lambda x: x[0][1,:], lambda x: matrix(x)[1,:],[x],x0,name="x[1,:]") 478 self.numpyEvaluationCheck(lambda x: x[0][0,:], lambda x: matrix(x)[0,:],[x],x0,name="x[0,:]") 479 self.numpyEvaluationCheck(lambda x: x[0][-1,:], lambda x: matrix(x)[-1,:],[x],x0,name="x[-1,:]") 480 self.numpyEvaluationCheck(lambda x: x[0][:,-2], lambda x: matrix(x)[:,-2],[x],x0,name="x[:,-2]") 481 self.numpyEvaluationCheck(lambda x: x[0][0:-2,0:-1], lambda x: matrix(x)[0:-2,0:-1],[x],x0,name="x[0:-2,0:-1]") 482 self.numpyEvaluationCheck(lambda x: x[0][0:2,0:2], lambda x: matrix(x)[0:2,0:2],[x],x0,name="x[0:2,0:2]") 483 self.numpyEvaluationCheck(lambda x: x[0][[0,1],0:2], lambda x: matrix(x)[[0,1],0:2],[x],x0,name="x[[0,1],0:2]") 484 self.numpyEvaluationCheck(lambda x: x[0].nz[[0,2,3]], lambda x: matrix([x[0,0],x[2,0],x[0,1]]).T,[x],x0,name="x[[0,2,3]]") 485 486 self.numpyEvaluationCheck(lambda x: x[0].nz[1], lambda x: matrix(x.T.ravel()[1]).T,[x],x0,name="x[1] on dense matrix") 487 self.numpyEvaluationCheck(lambda x: x[0].nz[-1], lambda x: matrix(x.ravel()[-1]).T,[x],x0,name="x[-1] on dense matrix") 488 self.numpyEvaluationCheck(lambda x: x[0][[0,1],0:1],lambda x: x[[0,1],0:1],[x],x0,name='x[:,0:1]') 489 self.numpyEvaluationCheck(lambda x: x[0][0:1,[0,1]],lambda x: x[0:1,[0,1]],[x],x0,name='x[0:1,:]') 490 491 self.message(":sparse") 492 493 sp=Sparsity(4,3,[0,2,2,3],[1,2,1]) 494 x=MX.sym("X",sp) 495 sx0=[0.738,0.39,0.99] 496 x0=DM(Sparsity(4,3,[0,2,2,3],[1,2,1]),[0.738,0.39,0.99]).full() 497 self.numpyEvaluationCheck(lambda x: x[0][0,0], lambda x: matrix(x)[0,0],[x],x0,name="x[0,0]",setx0=[sx0]) 498 self.numpyEvaluationCheck(lambda x: x[0][1,0], lambda x: matrix(x)[1,0],[x],x0,name="x[1,0]",setx0=[sx0]) 499 self.numpyEvaluationCheck(lambda x: x[0][0,1], lambda x: matrix(x)[0,1],[x],x0,name="x[1,0]",setx0=[sx0]) 500 self.numpyEvaluationCheck(lambda x: x[0][0,-1], lambda x: matrix(x)[0,-1],[x],x0,name="x[0,-1]",setx0=[sx0]) 501 self.numpyEvaluationCheck(lambda x: x[0][:,0], lambda x: matrix(x)[:,0],[x],x0,name="x[:,0]",setx0=[sx0]) 502 self.numpyEvaluationCheck(lambda x: x[0][:,1], lambda x: matrix(x)[:,1],[x],x0,name="x[:,1]",setx0=[sx0]) 503 self.numpyEvaluationCheck(lambda x: x[0][1,:], lambda x: matrix(x)[1,:],[x],x0,name="x[1,:]",setx0=[sx0]) 504 self.numpyEvaluationCheck(lambda x: x[0][0,:], lambda x: matrix(x)[0,:],[x],x0,name="x[0,:]",setx0=[sx0]) 505 self.numpyEvaluationCheck(lambda x: x[0][-1,:], lambda x: matrix(x)[-1,:],[x],x0,name="x[-1,:]",setx0=[sx0]) 506 self.numpyEvaluationCheck(lambda x: x[0][:,-2], lambda x: matrix(x)[:,-2],[x],x0,name="x[:,-2]",setx0=[sx0]) 507 self.numpyEvaluationCheck(lambda x: x[0][0:-2,0:-1], lambda x: matrix(x)[0:-2,0:-1],[x],x0,name="x[0:-2,0:-1]",setx0=[sx0]) 508 self.numpyEvaluationCheck(lambda x: x[0][0:2,0:2], lambda x: matrix(x)[0:2,0:2],[x],x0,name="x[0:2,0:2]",setx0=[sx0]) 509 self.numpyEvaluationCheck(lambda x: x[0][[0,1],0:2], lambda x: matrix(x)[[0,1],0:2],[x],x0,name="x[[0,1],0:2]",setx0=[sx0]) 510 self.numpyEvaluationCheck(lambda x: x[0].nz[[2,1]], lambda x: matrix([x[1,2],x[2,0]]).T,[x],x0,name="x[[2,1]]") 511 self.numpyEvaluationCheck(lambda x: x[0].nz[0:2], lambda x: matrix(sx0[0:2]).T,[x],x0,name="x[0:2] on dense matrix") 512 self.numpyEvaluationCheck(lambda x: x[0].nz[1], lambda x: matrix(sx0[1]).T,[x],x0,name="x[1]",setx0=[sx0]) 513 self.numpyEvaluationCheck(lambda x: x[0].nz[-1], lambda x: matrix(sx0[-1]).T,[x],x0,name="x[-1]",setx0=[sx0]) 514 515 def test_mx_in(self): 516 self.message("mx_out/mx_in") 517 x=MX.sym("x",2,3) 518 f = Function("f", [x],[3*x]) 519 x_in = f.mx_in() 520 x_out = f.call(x_in) 521 g = Function("g", [x_in[0]],[6*x_out[0]]) 522 n=[1,2,3,4,5,6] 523 f_in=DM(f.sparsity_in(0),n) 524 f_out = f(f_in) 525 g_in=DM(g.sparsity_in(0),n) 526 g_out = g(g_in) 527 checkarray(self,6*f_out.full(),g_out.full(),"slicing(trans)") 528 529 def test_scalarMX(self): 530 x=MX.sym("x") 531 x0=0.738 532 self.numpyEvaluationCheckPool(self.pool,[x],x0,name="scalarMX") 533 534 self.numpyEvaluationCheckPool(self.matrixpool,[x],x0,name="scalarMX") 535 536 def test_MXJacobian(self): 537 self.message("MX(1,1) unary operation, jacobian") 538 self.Jpool=FunctionPool() 539 self.message("SX(1,1) unary operation, jacobian") 540 x=MX.sym("x") 541 x0=array([[0.738]]) 542 543 def fmod(f,x): 544 J=f.jacobian_old(0, 0) 545 return J 546 547 self.numpyEvaluationCheckPool(self.Jpool,[x],x0,name="MX unary operations, jacobian",fmod=fmod) 548 549 def test_MXJacobians(self): 550 self.message("MX(3,1) unary operation, jacobian") 551 x=MX.sym("x",3,1) 552 553 x0=array([0.738,0.9,0.3]) 554 555 def fmod(f,x): 556 J=f.jacobian_old(0, 0) 557 return J 558 559 self.numpyEvaluationCheckPool(self.Jpool,[x],x0,name="MX unary operations, jacobian",fmod=fmod) 560 561 def test_MXbinary(self): 562 self.message("MX binary operations") 563 x=MX.sym("x",3,2) 564 y=MX.sym("x",3,2) 565 x0=array([[0.738,0.2],[ 0.1,0.39 ],[0.99,0.999999]]) 566 y0=array([[1.738,0.6],[ 0.7,12 ],[0,-6]]) 567 self.numpyEvaluationCheckPool(self.matrixbinarypool,[x,y],[x0,y0],name="MX") 568 569 def test_MXSparse(self): 570 self.message("MX unary operations, sparse") 571 sp=Sparsity(4,3,[0,2,2,3],[1,2,1]) 572 573 x=MX.sym("x",sp) 574 if scipy_available: 575 x0=DM(Sparsity(4,3,[0,2,2,3],[1,2,1]),[0.738,0.1,0.99]).sparse() 576 577 self.numpyEvaluationCheckPool(self.pool,[x],array(x0.todense()),name="MX",setx0=x0,excludeflags={'nozero'}) 578 self.numpyEvaluationCheckPool(self.matrixpool,[x],array(x0.todense()),name="MX",setx0=x0) 579 else: 580 x0=DM(Sparsity(4,3,[0,2,2,3],[1,2,1]),[0.738,0.1,0.99]).full() 581 582 self.numpyEvaluationCheckPool(self.pool,[x],x0,name="MX",setx0=x0,excludeflags={'nozero'}) 583 self.numpyEvaluationCheckPool(self.matrixpool,[x],x0,name="MX",setx0=x0) 584 585 def test_MXbinarySparse(self): 586 self.message("SX binary operations") 587 spx=Sparsity(4,3,[0,2,2,3],[1,2,1]) 588 spy=Sparsity(4,3,[0,2,2,3],[0,2,3]) 589 xx=MX.sym("x",spx) 590 yy=MX.sym("y",spy) 591 if scipy_available: 592 x0=DM(Sparsity(4,3,[0,2,2,3],[1,2,1]),[0.738,0.1,0.99]).sparse() 593 y0=DM(Sparsity(4,3,[0,2,2,3],[0,2,3]),[1.738,0.7,-6]).sparse() 594 595 self.numpyEvaluationCheckPool(self.matrixbinarypool,[xx,yy],[array(x0.todense()),array(y0.todense())],name="MX",setx0=[x0,y0]) 596 else: 597 x0=DM(Sparsity(4,3,[0,2,2,3],[1,2,1]),[0.738,0.1,0.99]).full() 598 y0=DM(Sparsity(4,3,[0,2,2,3],[0,2,3]),[1.738,0.7,-6]).full() 599 600 self.numpyEvaluationCheckPool(self.matrixbinarypool,[xx,yy],[x0,y0],name="MX",setx0=[x0,y0]) 601 602 def test_symbolcheck(self): 603 self.message("Check if non-symbolic inputs are caught") 604 self.assertRaises(RuntimeError, lambda : Function("f", [MX(0)],[MX.sym("x")])) 605 606 def test_unite(self): 607 self.message("unite operation") 608 import numpy 609 numpy.random.seed(42) 610 xn = numpy.random.random((3,4)) 611 x=MX(3,4) 612 y=MX.sym("x",3,4) 613 z=unite(x,y) 614 f = Function("f", [y],[z]) 615 f_in = [0]*f.n_in();f_in[0]=xn 616 f_out = f(*f_in) 617 self.checkarray(f_out,xn,"unite dense") 618 619 spx=Sparsity(4,3,[0,2,2,3],[1,2,1]) 620 spy=Sparsity(4,3,[0,1,2,3],[0,2,2]) 621 622 nx=DM.zeros(spx) 623 for k in range(nx.nnz()): 624 nx.nz[k]= numpy.random.rand() 625 ny=DM.zeros(spy) 626 for k in range(nx.nnz()): 627 ny.nz[k]= numpy.random.rand() 628 629 nxn = nx.full() 630 nyn = ny.full() 631 x=MX.sym("x",spx) 632 y=MX.sym("y",spy) 633 z=unite(x,y) 634 635 f = Function("f", [x,y],[z]) 636 f_in = [0]*f.n_in();f_in[0]=nx 637 f_in[1]=ny 638 f_out = f(*f_in) 639 self.checkarray(f_out,nxn+nyn,"unite sparse") 640 641 def test_imatrix_index(self): 642 self.message("IM indexing") 643 X = MX.sym("x",2,2) 644 Y = X.nz[np.array([[0,2],[1,1],[3,3]])] 645 646 f = Function("f", [X],[Y]) 647 f_in = [0]*f.n_in();f_in[0]=DM(f.sparsity_in(0),[1,2,3,4]) 648 f_out = f(*f_in) 649 650 self.checkarray(f_out,array([[1,3],[2,2],[4,4]]),"IM indexing") 651 652 Y = X[:,:] 653 Y.nz[np.array([[0,2]])] = DM([[9,8]]) 654 655 f = Function("f", [X],[Y]) 656 f_in = DM(f.sparsity_in(0),[1,2,3,4]) 657 f_out = f(f_in) 658 659 self.checkarray(f_out,array([[9,8],[2,4]]),"IM indexing assignment") 660 661 def test_subsass(self): 662 self.message("Check subscripted assignment") 663 664 X = MX.sym("x",2,2) 665 X[0,0]=MX(5) 666 X[0,0]=5 667 X[:,0]=8 668 669 x=MX.sym("X",3,4) 670 import numpy 671 numpy.random.seed(42) 672 xn = numpy.random.random((3,4)) 673 r = numpy.zeros((7,8)) 674 y=MX.zeros(7,8) 675 y[1:4,[2,4,6,7]]=x 676 r[1:4,[2,4,6,7]]=xn 677 fy = Function("fy", [x],[y]) 678 fy_in = [0]*fy.n_in();fy_in[0]=xn 679 fy_out = fy(*fy_in) 680 681 self.checkarray(fy_out,r,"subscripted assigment") 682 683 y=MX(7,8) 684 y[1:4,[2,4,6,7]]=x 685 r[1:4,[2,4,6,7]]=xn 686 fy = Function("fy", [x],[y]) 687 fy_out = fy(xn) 688 self.checkarray(fy_out,r,"subscripted assigment") 689 690 kl=[2,4,5,8] 691 692 s=y.sparsity() 693 for k in kl: 694 r[s.row()[k],s.get_col()[k]]=1.0 695 696 y.nz[kl]=MX(1) 697 fy = Function("fy", [x],[y]) 698 fy_out = fy(xn) 699 self.checkarray(fy_out,r,"subscripted assigment") 700 701 y.nz[kl]=x.nz[[0,1,2,3]] 702 s=y.sparsity() 703 sx=x.sparsity() 704 cnt=0 705 for k in kl: 706 r[s.row()[k],s.get_col()[k]]=xn[sx.row()[cnt],sx.get_col()[cnt]] 707 cnt+=1 708 fy = Function("fy", [x],[y]) 709 fy_out = fy(xn) 710 self.checkarray(fy_out,r,"subscripted assigment") 711 712 def test_erase(self): 713 self.message("Erase function") 714 self.message(":dense") 715 y=MX.sym("Y",7,8) 716 import numpy 717 r=2*numpy.ones((7,8)) 718 r[1:4,[2,4,6,7]]=numpy.zeros((3,4)) 719 z = y *2 720 z.erase([1,2,3],[2,4,6,7]) 721 f = Function("f", [y],[z]) 722 f_in = [0]*f.n_in();f_in[0]=DM(f.sparsity_in(0),[1]*56) 723 f_out = f(0) 724 e = f_out 725 self.checkarray(f_out,e,"erase") # fishy 726 self.message(":sparse") 727 728 def test_MXalgebraDense(self): 729 self.message("Test some dense algebraic properties of matrices") 730 # issue 96 731 n = 3 732 m = 4 733 import numpy 734 numpy.random.seed(42) 735 A_ = numpy.random.random((m,n)) 736 A = MX.sym("A",m,n) 737 b_ = numpy.random.random((m,1)) 738 b = MX.sym("b",m,1) 739 C_ = numpy.random.random((m,m)) 740 C = MX.sym("C",m,m) 741 D_ = numpy.random.random((m,n)) 742 D = MX.sym("D",m,n) 743 e_ = numpy.random.random((m,1)) 744 e = MX.sym("e",m,1) 745 x_ = numpy.random.random((n,1)) 746 x = MX.sym("x",n,1) 747 748 Axb = mtimes(A,x)+b 749 Dxe = mtimes(D,x)+e 750 a = mtimes(mtimes(Axb.T,C),Dxe) 751 752 f = Function("f", [x,A,b,C,D,e],[a]) 753 f_out = f(x_, A_, b_, C_, D_, e_) 754 755 f_ = numpy.dot(numpy.dot((numpy.dot(A_,x_)+b_).T,C_),(numpy.dot(D_,x_)+e_)) 756 757 self.checkarray(f_out,f_,"evaluation") 758 759 760 J_ = numpy.dot(numpy.dot((numpy.dot(D_,x_)+e_).T,C_.T),A_) + numpy.dot(numpy.dot((numpy.dot(A_,x_)+b_).T,C_),D_) 761 762 for w in [0, 1]: 763 f = Function("f", [x,A,b,C,D,e], [a], {"ad_weight":w, "ad_weight_sp":w}) 764 J = f.jacobian_old(0, 0) 765 J_in = [0]*J.n_in();J_in[0]=x_ 766 J_in[1]=A_ 767 J_in[2]=b_ 768 J_in[3]=C_ 769 J_in[4]=D_ 770 J_in[5]=e_ 771 J_out = J.call(J_in) 772 773 self.checkarray(J_out[0],J_,"evaluation") 774 775 def test_MXalgebraSparse(self): 776 self.message("Test some sparse algebraic properties of matrices") 777 if not(scipy_available): 778 return 779 # issue 97 780 n = 3 781 m = 4 782 import numpy 783 numpy.random.seed(42) 784 785 def randsparsity(m,n): 786 sp = Sparsity(m,n) 787 for i in range(int((n*m)/2)): 788 sp.add_nz(numpy.random.randint(m),numpy.random.randint(n)) 789 return sp 790 791 def gentest(m,n): 792 As = randsparsity(m,n) 793 A_ = DM.zeros(As) 794 for k in range(As.nnz()): 795 A_.nz[k]= numpy.random.rand() 796 A = MX.sym("A",As) 797 return (A_.sparse(),A) 798 799 (A_,A)=gentest(m,n) 800 (b_,b)=gentest(m,1) 801 (C_,C)=gentest(m,m) 802 (D_,D)=gentest(m,n) 803 (e_,e)=gentest(m,1) 804 x_ = numpy.random.random((n,1)) 805 x = MX.sym("x",n,1) 806 807 Axb = mtimes(A,x)+b 808 Dxe = mtimes(D,x)+e 809 a = mtimes(mtimes(Axb.T,C),Dxe) 810 811 f = Function("f", [x,A,b,C,D,e],[a]) 812 f_in = [0]*f.n_in();f_in[0]=x_ 813 f_in[1]=A_ 814 f_in[2]=b_ 815 f_in[3]=C_ 816 f_in[4]=D_ 817 f_in[5]=e_ 818 f_out = f(*f_in) 819 820 821 Axb_ = A_*x_+b_ 822 Dxe_ = D_*x_+e_ 823 824 f_ = Axb_.T*C_*Dxe_ 825 826 self.checkarray(f_out,f_,"evaluation") 827 828 829 J_ = (D_*x_+e_).T*C_.T*A_ + (A_*x_+b_).T*C_*D_ 830 831 for w in [0, 1]: 832 f = Function("f", [x,A,b,C,D,e], [a], {"ad_weight":w, "ad_weight_sp":w}) 833 834 J = f.jacobian_old(0, 0) 835 J_in = [0]*J.n_in();J_in[0]=x_ 836 J_in[1]=A_ 837 J_in[2]=b_ 838 J_in[3]=C_ 839 J_in[4]=D_ 840 J_in[5]=e_ 841 J_out = J.call(J_in) 842 843 self.checkarray(J_out[0],J_,"evaluation") 844 845 #@unittest.skipIf(not(scipy_available)) 846 def test_MXalgebraSparseSparse(self): 847 if not(scipy_available): 848 return 849 self.message("Test some sparse algebraic properties of matrices") 850 n = 8 851 m = 10 852 import numpy 853 numpy.random.seed(42) 854 855 def randsparsity(m,n): 856 sp = Sparsity(m,n) 857 for k in range(int((n*m)/2)): 858 i = numpy.random.randint(m) 859 j = numpy.random.randint(n) 860 if not(i == int(m/2)): 861 if n==1 or not(j == int(n/2)): 862 sp.add_nz(i,j) 863 return sp 864 865 def gentest(m,n): 866 As = randsparsity(m,n) 867 A_ = DM.zeros(As) 868 for k in range(As.nnz()): 869 A_.nz[k]= numpy.random.rand() 870 A = MX.sym("A",As) 871 return (A_.sparse(),A) 872 873 (A_,A)=gentest(m,n) 874 (b_,b)=gentest(m,1) 875 (C_,C)=gentest(m,m) 876 (D_,D)=gentest(m,n) 877 (e_,e)=gentest(m,1) 878 x_ = numpy.random.random((n,1)) 879 x = MX.sym("x",n,1) 880 881 Axb = mtimes(A,x)+b 882 Dxe = mtimes(D,x)+e 883 a = mtimes(mtimes(Axb.T,C),Dxe) 884 885 f = Function("f", [x,A,b,C,D,e],[a]) 886 f_in = [0]*f.n_in();f_in[0]=x_ 887 f_in[1]=A_ 888 f_in[2]=b_ 889 f_in[3]=C_ 890 f_in[4]=D_ 891 f_in[5]=e_ 892 f_out = f(*f_in) 893 894 895 Axb_ = A_*x_+b_ 896 Dxe_ = D_*x_+e_ 897 898 f_ = Axb_.T*C_*Dxe_ 899 900 self.checkarray(f_out,f_,"evaluation") 901 902 903 J_ = (D_*x_+e_).T*C_.T*A_ + (A_*x_+b_).T*C_*D_ 904 905 for w in [0, 1]: 906 f = Function("f", [x,A,b,C,D,e], [a], {"ad_weight":w, "ad_weight_sp":w}) 907 J = f.jacobian_old(0, 0) 908 J_in = [0]*J.n_in();J_in[0]=x_ 909 J_in[1]=A_ 910 J_in[2]=b_ 911 J_in[3]=C_ 912 J_in[4]=D_ 913 J_in[5]=e_ 914 J_out = J.call(J_in) 915 916 self.checkarray(J_out[0],J_,"evaluation") 917 918 919 def test_chaining(self): 920 self.message("Chaining SX and MX together") 921 x=SX.sym("x") 922 y=x**3 923 f=Function("f", [x],[y]) 924 J=f.jacobian_old(0, 0) 925 926 X=MX.sym("X") 927 F=Function("F", [X], J(X)) 928 929 x_=1.7 930 F_out = F(x_) 931 self.checkarray(F_out[0],3*x_**2,"Chaining eval") 932 933 def test_issue107(self): 934 self.message("Regression test for issue 107: +=") 935 x=MX.sym("x") 936 y=MX.sym("y") 937 938 z=x 939 z+=y 940 941 self.assertTrue(x.is_symbolic()) 942 self.assertFalse(z.is_symbolic()) 943 944 def test_MXd_trivial(self): 945 self.message("symbolic variables and constants jac") 946 X = MX.sym("X",10) 947 V = MX.sym("V") 948 J = jacobian(X,X) 949 self.assertTrue(isinstance(J,MX)) 950 self.assertEqual(J.nnz(),10) 951 self.assertEqual(J.size1(),10) 952 self.assertEqual(J.size2(),10) 953 954 g = Function("g", [],[J]) 955 [g_out] = g.call([]) 956 self.checkarray(g_out,eye(10),"unit matrix") 957 g = Function("g", [],[jacobian(MX.eye(3),X)]) 958 [g_out] = g.call([]) 959 self.checkarray(g_out,zeros((9,10)),"zero matrix") 960 g = Function("g", [],[jacobian(X,V)]) 961 [g_out] = g.call([]) 962 self.checkarray(g_out,zeros((10,1)),"zero matrix") 963 964 g = Function("g", [],[jacobian(MX.eye(3),V)]) 965 [g_out] = g.call([]) 966 self.checkarray(g_out,zeros((9,1)),"zero matrix") 967 968 def test_MXd_substractionl(self): 969 self.message("substraction jac") 970 V = MX.sym("V") 971 X = MX.sym("X") 972 g = Function("g", [],[jacobian(X-V, X)]) 973 [g_out] = g.call([]) 974 self.checkarray(g_out,ones((1,1)), "one") 975 976 g = Function("g", [],[jacobian(X-V, V)]) 977 [g_out] = g.call([]) 978 self.checkarray(g_out,-ones((1,1)), "one") 979 980 g = Function("g", [],[jacobian(V-X, X)]) 981 [g_out] = g.call([]) 982 self.checkarray(g_out,-ones((1,1)), "one") 983 984 g = Function("g", [],[jacobian(V-X, V)]) 985 [g_out] = g.call([]) 986 self.checkarray(g_out,ones((1,1)),"one") 987 988 def test_MXd_mapping(self): 989 self.message("mapping jac") 990 X = MX.sym("X",3) 991 Y = MX.sym("Y",2) 992 J = jacobian(vertcat(X,Y),X) 993 JJ = DM.ones(J.sparsity()) 994 self.checkarray(JJ,numpy.vstack((eye(3),zeros((2,3)))),"diag") 995 J = jacobian(vertcat(X,Y),Y) 996 JJ = DM.ones(J.sparsity()) 997 self.checkarray(JJ,numpy.vstack((zeros((3,2)),eye(2))),"diag") 998 999 def test_null(self): 1000 self.message("Function null") 1001 x = MX.sym("x") 1002 1003 f = Function("f", [x],[x**2,MX()]) 1004 1005 self.assertEqual(f.size1_out(1),0) 1006 self.assertEqual(f.size2_out(1),0) 1007 f_out = f(0) 1008 1009 f = Function("f", [x,MX()],[x**2,MX()]) 1010 1011 self.assertEqual(f.size1_out(1),0) 1012 self.assertEqual(f.size2_out(1),0) 1013 f_out = f(0,0) 1014 1015 r = f(x,MX()) 1016 self.assertTrue(r[1].is_empty(True)) 1017 1018 r = f(MX(),MX()) 1019 self.assertTrue(r[1].is_empty(True)) 1020 1021 #self.assertRaises(Exception,lambda : f([x,x],True)) 1022 #self.assertRaises(Exception,lambda : f([[],[]],True)) 1023 1024 def test_issue184(self): 1025 self.message("Regression test issue #184") 1026 x = MX.sym("x", 3) 1027 y = x[0:0] 1028 self.assertEqual(y.nnz(),0) 1029 1030 def test_indexingOutOfBounds(self): 1031 self.message("Indexing out of bounds") 1032 y = DM.zeros(4, 5) 1033 self.assertRaises(RuntimeError,lambda : y[12,0] ) 1034 self.assertRaises(RuntimeError,lambda : y[12,12] ) 1035 self.assertRaises(RuntimeError,lambda : y[0,12] ) 1036 self.assertRaises(RuntimeError,lambda : y[12,:] ) 1037 self.assertRaises(RuntimeError,lambda : y[12:15,0] ) 1038 self.assertRaises(RuntimeError,lambda : y[:,12] ) 1039 self.assertRaises(RuntimeError,lambda : y[0,12:15] ) 1040 y[-1,2] 1041 self.assertRaises(RuntimeError,lambda : y[-12,2] ) 1042 y[-3:-1,2] 1043 self.assertRaises(RuntimeError,lambda : y[-12:-9,2] ) 1044 1045 def test(): 1046 y[12,0] = 0 1047 self.assertRaises(RuntimeError,test) 1048 def test(): 1049 y[12,12] = 0 1050 self.assertRaises(RuntimeError,test) 1051 def test(): 1052 y[0,12] = 0 1053 self.assertRaises(RuntimeError,test) 1054 def test(): 1055 y[12,:] = 0 1056 self.assertRaises(RuntimeError,test) 1057 def test(): 1058 y[12:15,0] = 0 1059 self.assertRaises(RuntimeError,test) 1060 def test(): 1061 y[:,12] = 0 1062 self.assertRaises(RuntimeError,test) 1063 def test(): 1064 y[0,12:15] = 0 1065 self.assertRaises(RuntimeError,test) 1066 y[-1,2] = 0 1067 def test(): 1068 y[-12,2] = 0 1069 self.assertRaises(RuntimeError,test) 1070 y[-3:-1,2] = 0 1071 def test(): 1072 y[-12:-9,2] = 0 1073 self.assertRaises(RuntimeError,test) 1074 1075 def test_indexinglimits(self): 1076 self.message("Limits of indexing") 1077 y = casadi.MX.sym("y", 3) 1078 self.assertRaises(RuntimeError,lambda : y[[0, 5]] ) 1079 try: 1080 y[[0, 5]] = MX.sym("a") 1081 self.assertTrue(False) 1082 except RuntimeError: 1083 pass 1084 y[[0, 2]] 1085 y[[0, 2]] = MX.sym("a") 1086 1087 def test_mtimes(self): 1088 A = MX(DM.ones((4,3))) 1089 B = MX(DM.ones((3,8))) 1090 C = MX(DM.ones((8,7))) 1091 1092 self.assertRaises(RuntimeError,lambda : mtimes([])) 1093 1094 D = mtimes([A]) 1095 1096 self.assertEqual(D.shape[0],4) 1097 self.assertEqual(D.shape[1],3) 1098 1099 D = mtimes([A,B]) 1100 1101 self.assertEqual(D.shape[0],4) 1102 self.assertEqual(D.shape[1],8) 1103 1104 D = mtimes([A,B,C]) 1105 1106 self.assertEqual(D.shape[0],4) 1107 self.assertEqual(D.shape[1],7) 1108 1109 def test_truth(self): 1110 self.message("Truth values") 1111 self.assertRaises(Exception, lambda : bool(MX.sym("x"))) 1112 #self.assertRaises(Exception, lambda : bool(MX.sym("x")>0)) 1113 self.assertTrue(bool(MX(1))) 1114 self.assertFalse(bool(MX(0))) 1115 self.assertTrue(bool(MX(0.2))) 1116 self.assertTrue(bool(MX(-0.2))) 1117 self.assertRaises(Exception, lambda : bool(MX(DM([2.0,3])))) 1118 self.assertRaises(Exception, lambda : bool(MX())) 1119 1120 1121 def test_MXbool(self): 1122 self.message("bool") 1123 1124 xy = MX.sym("x",2) 1125 x = xy[0] 1126 y = xy[1] 1127 1128 f = Function("f", [xy],[vertcat(*[logic_and(x,y),logic_or(x,y),logic_not(x)])]) 1129 1130 1131 for t1 in [0,1]: 1132 for t2 in [0,1]: 1133 T1 = t1!=0 1134 T2 = t2!=0 1135 f_out = f([t1,t2]) 1136 self.checkarray(f_out,DM([T1 and T2,T1 or T2,not T1]),"bool(%d,%d): %s" % (t1,t2,str(f_out))) 1137 1138 def test_MXineq(self): 1139 self.message("SX ineq") 1140 1141 xy = MX.sym("x",2) 1142 x = xy[0] 1143 y = xy[1] 1144 1145 1146 f = Function("f", [xy],[vertcat(*[x<y,x<=y,x>=y,x==y,x!=y])]) 1147 1148 for t1 in [-10,0.1,0,1,10]: 1149 for t2 in [-10,0.1,0,1,10]: 1150 T1 = t1 1151 T2 = t2 1152 f_out = f([t1,t2]) 1153 self.checkarray(f_out,DM([T1 < T2,T1 <= T2, T1 >= T2, T1 == T2, T1 != T2]),"ineq(%d,%d)" % (t1,t2)) 1154 1155 def test_if_else_zero(self): 1156 x = MX.sym("x") 1157 y = if_else(x,5,0) 1158 f = Function("f", [x],[y]) 1159 f_in = 1 1160 f_out = f(f_in) 1161 self.assertTrue(f_out==5,"if_else_zero %s " % str(f_out)) 1162 f_in = 0 1163 f_out = f(f_in) 1164 self.assertTrue(f_out==0,"if_else_zero") 1165 1166 1167 def test_if_else(self): 1168 x = MX.sym("x") 1169 y = if_else(x,1,2) 1170 f = Function("f", [x],[y]) 1171 f_in = 1 1172 f_out = f(f_in) 1173 self.assertTrue(f_out==1,"if_else") 1174 f_in = 0 1175 f_out = f(f_in) 1176 self.assertTrue(f_out==2,"if_else") 1177 1178 def test_regression491(self): 1179 self.message("regression #491") 1180 u = SX.sym("u") 1181 x = SX.sym("x") 1182 1183 F = Function("F", [u,x],[u+1/x]) 1184 1185 U = MX.sym("U") 1186 1187 X = F(U,U) 1188 G = F(U,X) 1189 1190 for kk in range(2): 1191 gfcn = 0 1192 if kk==0: 1193 gfcn = Function("gfcn", [U], [G]).expand("e_gfcn", {"ad_weight":1}) 1194 else: 1195 gfcn = Function("gfcn", [U],[G], {"ad_weight":1}) 1196 J = gfcn.jacobian_old(0, 0) 1197 J_in = [0]*J.n_in();J_in[0]=1 1198 J_out = J.call(J_in) 1199 self.assertAlmostEqual(J_out[0],1,9) 1200 1201 def test_ticket(self): 1202 J = [] + MX.sym("x") 1203 J = MX.sym("x") + [] 1204 1205 def test_jacobian_tools(self): 1206 self.message("jacobian") 1207 1208 X = MX.sym("X") 1209 1210 Y = jacobian(X**2,X) 1211 1212 f = Function("f", [X], [Y]) 1213 1214 f_in=2.3 1215 f_out = f(f_in) 1216 1217 self.assertAlmostEqual(f_out,4.6) 1218 1219 def test_reshape(self): 1220 self.message("reshape") 1221 X = MX.sym("X",10) 1222 1223 i = DM(Sparsity.lower(3),list(range(6))) 1224 1225 i.print_dense() 1226 print(i.T.nz[:]) 1227 1228 T = X.nz[i] 1229 1230 q = T.T.nz[:]**2 1231 f = Function("f", [X], [q]) 1232 f_out = f(list(range(10))) 1233 1234 self.checkarray(DM([0,1,9,4,16,25]),f_out) 1235 1236 Y = MX.sym("Y",10) 1237 1238 ff = Function("ff", [Y],f.call([Y],True)) 1239 ff_out = ff(list(range(10))) 1240 1241 self.checkarray(DM([0,1,9,4,16,25]),ff_out) 1242 1243 J = Function("J", [X],[jacobian(q, X)]) 1244 J_out = J(list(range(10))) 1245 1246 i = horzcat(*[diag([0,2,4,6,8,10]),DM.zeros(6,4)]) 1247 i[[2,3],:] = i[[3,2],:] 1248 1249 self.checkarray(i,J_out) 1250 1251 q = (T.T).nz[:]**2 1252 J = Function("J", [X],[jacobian(q,X)]) 1253 J_in = [0]*J.n_in();J_in[0]=list(range(10)) 1254 J_out = J(*J_in) 1255 1256 i = horzcat(*[diag([0,2,4,6,8,10]),DM.zeros(6,4)]) 1257 i[[2,3],:] = i[[3,2],:] 1258 1259 self.checkarray(i,J_out) 1260 1261 def test_vertcat(self): 1262 self.message("vertcat") 1263 X = MX.sym("X",10) 1264 1265 T = vertcat(*[X[4],X[2]]) 1266 q = T**2 1267 f = Function("f", [X],[q]) 1268 f_in = [0]*f.n_in();f_in[0]=list(range(10)) 1269 f_out = f.call(f_in) 1270 1271 self.checkarray(DM([16,4]),f_out[0]) 1272 1273 Y = MX.sym("Y",10) 1274 1275 ff = Function("ff", [Y],f.call([Y],True)) 1276 ff_in = [0]*ff.n_in();ff_in[0]=list(range(10)) 1277 ff_out = ff(*ff_in) 1278 1279 self.checkarray(DM([16,4]),ff_out) 1280 J = Function("J", [X],[jacobian(q,X)]) 1281 J_in = [0]*J.n_in();J_in[0]=list(range(10)) 1282 J_out = J(*J_in) 1283 1284 i = DM.zeros(2,10) 1285 i[0,4] = 8 1286 i[1,2] = 4 1287 1288 self.checkarray(i,J_out) 1289 q = T**2 1290 J = Function("J", [X],[jacobian(q, X)]) 1291 J_in = [0]*J.n_in();J_in[0]=list(range(10)) 1292 J_out = J(*J_in) 1293 1294 self.checkarray(i,J_out) 1295 1296 def test_blockcat(self): 1297 x = MX.sym("x") 1298 1299 y = blockcat([[x,2*x],[3*x,4*x]]) 1300 f = Function("f", [x],[y]) 1301 f_out = f(3) 1302 self.checkarray(f_out,DM([[3,6],[9,12]])) 1303 1304 1305 def test_veccats(self): 1306 x= MX.sym("x",2) 1307 self.assertTrue(hash(vec(x))==hash(x)) 1308 1309 def test_constmxmtimes(self): 1310 0.1*MX.ones(2) 1311 1312 def test_is_regular(self): 1313 self.assertTrue(MX(DM([0,1])).is_regular()) 1314 self.assertFalse(MX(DM([0,inf])).is_regular()) 1315 with self.assertRaises(Exception): 1316 self.assertFalse(MX.sym("x",2).is_regular()) 1317 1318 def test_diagcat(self): 1319 C = diagcat(*[MX(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]) 1320 self.assertTrue(isinstance(C,MX)) 1321 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]]) 1322 r = sparsify(r) 1323 f = Function("f", [],[C]) 1324 f_out = f.call([]) 1325 1326 self.checkarray(f_out[0],r) 1327 1328 def test_tril2symm(self): 1329 x = MX.sym("x",Sparsity.lower(3)) 1330 f = Function("f", [x],[tril2symm(x)]) 1331 f_in = [0]*f.n_in();f_in[0]=DM(f.sparsity_in(0),list(range(6))) 1332 f_out = f(*f_in) 1333 self.checkarray(f_out,DM([[0,1,2],[1,3,4],[2,4,5]])) 1334 1335 def test_sparsity_indexing(self): 1336 self.message("sparsity") 1337 1338 B_ = DM([[1,2,3,4,5],[6,7,8,9,10]]) 1339 B = MX.sym("B",2,5) 1340 1341 A = DM([[1,1,0,0,0],[0,0,1,0,0]]) 1342 A = sparsify(A) 1343 sp = A.sparsity() 1344 import copy 1345 1346 def meval(m): 1347 f = Function("f", [B],[m]) 1348 f_in = [0]*f.n_in();f_in[0]=B_ 1349 f_out = f(*f_in) 1350 return f_out 1351 1352 self.checkarray(meval(B[sp]),DM([[1,2,0,0,0],[0,0,8,0,0]]),"sparsity indexing") 1353 1354 Bmod = copy.copy(B) 1355 Bmod[sp] = -4 1356 1357 self.checkarray(meval(Bmod),DM([[-4,-4,3,4,5],[6,7,-4,9,10]]),"sparsity indexing assignement") 1358 1359 Bmod = copy.copy(B) 1360 Bmod[sp] = 2*B 1361 1362 self.checkarray(meval(Bmod),DM([[2,4,3,4,5],[6,7,16,9,10]]),"Imatrix indexing assignement") 1363 1364 self.assertRaises(Exception, lambda : B[Sparsity.dense(4,4)]) 1365 1366 def test_symvar(self): 1367 a = MX.sym("a") 1368 b = MX.sym("b") 1369 c = MX.sym("c") 1370 e = cos(a*b) + c 1371 w = symvar(e) 1372 self.assertEqual(len(w),3) 1373 if GlobalOptions.getSimplificationOnTheFly(): 1374 self.assertTrue(is_equal(w[0],a)) 1375 self.assertTrue(is_equal(w[1],b)) 1376 self.assertTrue(is_equal(w[2],c)) 1377 1378 @known_bug() 1379 def test_vertcat_empty(self): 1380 a = MX(DM(0,2)) 1381 v = vertcat(*[a,a]) 1382 1383 self.assertEqual(v.size1(),0) 1384 self.assertEqual(v.size2(),2) 1385 1386 a = MX(DM(2,0)) 1387 v = vertcat(*[a,a]) 1388 1389 self.assertEqual(v.size1(),4) 1390 self.assertEqual(v.size2(),0) 1391 1392 def test_jacobian_empty(self): 1393 x = MX.sym("x",3) 1394 1395 s = jacobian(DM(0,0),x).shape 1396 self.assertEqual(s[0],0) 1397 self.assertEqual(s[1],3) 1398 1399 s = jacobian(x,MX.sym("x",0,4)).shape 1400 self.assertEqual(s[0],3) 1401 self.assertEqual(s[1],0) 1402 1403 def test_mul_sparsity(self): 1404 1405 N = 10 1406 x = MX.sym("x",N,N) 1407 y = MX.sym("y",N,N) 1408 1409 x_ = self.randDM(N,N) 1410 y_ = self.randDM(N,N) 1411 1412 filt = Sparsity.diag(N)+Sparsity.triplet(N,N,[1],[3]) 1413 1414 f = Function("f", [x,y],[mtimes(x,y)]) 1415 f_in = (x_, y_) 1416 g = Function("g", [x,y],[mac(x,y,MX.zeros(filt))]) 1417 g_in = (x_, y_) 1418 1419 f_out = f(*f_in) 1420 g_out = g(*g_in) 1421 1422 self.checkarray(DM.ones(filt),DM.ones(g.sparsity_out(0))) 1423 1424 self.checkarray(f_out[filt],g_out) 1425 1426 def test_mul_zero_wrong(self): 1427 with self.assertRaises(RuntimeError): 1428 mtimes(MX.sym("X",4,5),MX.zeros(3,2)) 1429 1430 def test_vertsplit(self): 1431 a = MX.sym("X",Sparsity.lower(5)) 1432 v = vertsplit(a,[0,2,4,5]) 1433 1434 Nr = int(5*6/2) 1435 1436 f = Function("f", [a],v) 1437 f_in = [0]*f.n_in();f_in[0]=DM(f.sparsity_in(0),list(range(Nr))) 1438 1439 f_out = f.call(f_in) 1440 v = [f_out[i] for i in range(len(v))] 1441 1442 self.assertEqual(len(v),3) 1443 self.checkarray(v[0],DM([[0,0,0,0,0],[1,5,0,0,0]])) 1444 self.checkarray(v[1],DM([[2,6,9,0,0],[3,7,10,12,0]])) 1445 self.checkarray(v[2],DM([[4,8,11,13,14]])) 1446 1447 v = vertsplit(a) 1448 1449 f = Function("f", [a],v) 1450 f_in = [0]*f.n_in();f_in[0]=DM(f.sparsity_in(0),list(range(Nr))) 1451 1452 f_out = f.call(f_in) 1453 v = [f_out[i] for i in range(len(v))] 1454 1455 self.assertEqual(len(v),a.size1()) 1456 self.checkarray(v[0],DM([[0,0,0,0,0]])) 1457 self.checkarray(v[1],DM([[1,5,0,0,0]])) 1458 self.checkarray(v[2],DM([[2,6,9,0,0]])) 1459 self.checkarray(v[3],DM([[3,7,10,12,0]])) 1460 self.checkarray(v[4],DM([[4,8,11,13,14]])) 1461 1462 v = vertsplit(a,2) 1463 1464 f = Function("f", [a],v) 1465 f_in = [0]*f.n_in();f_in[0]=DM(f.sparsity_in(0),list(range(Nr))) 1466 1467 f_out = f.call(f_in) 1468 v = [f_out[i] for i in range(len(v))] 1469 1470 self.assertEqual(len(v),3) 1471 self.checkarray(v[0],DM([[0,0,0,0,0],[1,5,0,0,0]])) 1472 self.checkarray(v[1],DM([[2,6,9,0,0],[3,7,10,12,0]])) 1473 self.checkarray(v[2],DM([[4,8,11,13,14]])) 1474 1475 v = vertsplit(a,[0,0,3,a.size1()]) 1476 1477 f = Function("f", [a],v) 1478 f_in = [0]*f.n_in();f_in[0]=DM(f.sparsity_in(0),list(range(Nr))) 1479 1480 f_out = f(*f_in) 1481 V = [f_out[i] for i in range(len(v))] 1482 1483 self.assertEqual(len(v),3) 1484 self.assertEqual(v[0].size1(),0) 1485 self.assertEqual(v[0].size2(),5) # why not 5? 1486 self.checkarray(V[1],DM([[0,0,0,0,0],[1,5,0,0,0],[2,6,9,0,0]])) 1487 self.checkarray(V[2],DM([[3,7,10,12,0],[4,8,11,13,14]])) 1488 1489 def test_horzsplit(self): 1490 a = MX.sym("X",Sparsity.lower(5)) 1491 v = horzsplit(a,[0,2,4,5]) 1492 1493 Nr = int(5*6/2) 1494 f = Function("f", [a],v) 1495 f_in = DM(f.sparsity_in(0),list(range(Nr))) 1496 1497 f_out = f(f_in) 1498 v = [f_out[i] for i in range(len(v))] 1499 self.assertEqual(len(v),3) 1500 self.checkarray(v[0],DM([[0,0],[1,5],[2,6],[3,7],[4,8]])) 1501 self.checkarray(v[1],DM([[0,0],[0,0],[9,0],[10,12],[11,13]])) 1502 self.checkarray(v[2],DM([[0],[0],[0],[0],[14]])) 1503 1504 v = horzsplit(a) 1505 1506 f = Function("f", [a],v) 1507 f_in = DM(f.sparsity_in(0),list(range(Nr))) 1508 f_out = f(f_in) 1509 v = [f_out[i] for i in range(len(v))] 1510 self.assertEqual(len(v),a.size1()) 1511 self.checkarray(v[0],DM([0,1,2,3,4])) 1512 self.checkarray(v[1],DM([0,5,6,7,8])) 1513 self.checkarray(v[2],DM([0,0,9,10,11])) 1514 self.checkarray(v[3],DM([0,0,0,12,13])) 1515 self.checkarray(v[4],DM([0,0,0,0,14])) 1516 1517 v = horzsplit(a,2) 1518 1519 f = Function("f", [a],v) 1520 f_in = DM(f.sparsity_in(0),list(range(Nr))) 1521 f_out = f(f_in) 1522 v = [f_out[i] for i in range(len(v))] 1523 1524 self.assertEqual(len(v),3) 1525 self.checkarray(v[0],DM([[0,0],[1,5],[2,6],[3,7],[4,8]])) 1526 self.checkarray(v[1],DM([[0,0],[0,0],[9,0],[10,12],[11,13]])) 1527 self.checkarray(v[2],DM([[0],[0],[0],[0],[14]])) 1528 1529 v = horzsplit(a,[0,0,3,a.size2()]) 1530 f = Function("f", [a],v) 1531 f_in = DM(f.sparsity_in(0),list(range(Nr))) 1532 f_out = f(f_in) 1533 V = [f_out[i] for i in range(len(v))] 1534 1535 self.assertEqual(len(v),3) 1536 self.assertEqual(v[0].size1(),5) 1537 self.assertEqual(v[0].size2(),0) 1538 self.checkarray(V[1],DM([[0,0,0],[1,5,0],[2,6,9],[3,7,10],[4,8,11]])) 1539 self.checkarray(V[2],DM([[0,0],[0,0],[0,0],[12,0],[13,14]])) 1540 1541 def test_blocksplit(self): 1542 a = MX.sym("X",Sparsity.lower(5)) 1543 v = blocksplit(a,[0,2,4,5],[0,1,3,5]) 1544 1545 Nr = int(5*6/2) 1546 fs = [Function("fs", [a],vr) for vr in v] 1547 v = [fs[i](DM(fs[i].sparsity_in(0),list(range(Nr)))) for i in range(3)] 1548 1549 self.checkarray(v[0][0],DM([0,1])) 1550 self.checkarray(v[0][1],DM([[0,0],[5,0]])) 1551 self.checkarray(v[1][0],DM([2,3])) 1552 self.checkarray(blockcat(v),DM(fs[0].sparsity_in(0),list(range(Nr)))) 1553 1554 def test_mxnulloutput(self): 1555 a = MX(5,0) 1556 b = MX.sym("x",2) 1557 1558 f = Function("f", [b],[a]) 1559 c = f(b) 1560 1561 self.assertEqual(c.size1(),5) 1562 self.assertEqual(c.size2(),0) 1563 1564 c = f.call([b],True)[0] 1565 1566 self.assertEqual(c.size1(),5) 1567 self.assertEqual(c.size2(),0) 1568 1569 a = MX(0,0) 1570 b = MX.sym("x",2) 1571 1572 f = Function("f", [b],[a]) 1573 c = f(b) 1574 1575 self.assertEqual(c.size1(),0) 1576 self.assertEqual(c.size2(),0) 1577 1578 c = f.call([b],True)[0] 1579 1580 self.assertEqual(c.size1(),0) 1581 self.assertEqual(c.size2(),0) 1582 1583 def test_mxnull(self): 1584 a = MX(5,0) 1585 b = MX(0,3) 1586 1587 c = mtimes(a,b) 1588 1589 self.assertEqual(c.nnz(),0) 1590 1591 a = MX(5,3) 1592 b = MX(3,4) 1593 1594 c = mtimes(a,b) 1595 1596 self.assertEqual(c.nnz(),0) 1597 1598 def test_mxnullop(self): 1599 c = MX(0,0) 1600 x = MX.sym("x",2,3) 1601 1602 # https://github.com/casadi/casadi/issues/2628 1603 if swig4: 1604 with self.assertRaises(TypeError): 1605 d = x + c 1606 else: 1607 with self.assertRaises(RuntimeError): 1608 d = x + c 1609 1610 if swig4: 1611 with self.assertRaises(TypeError): 1612 d = x / c 1613 else: 1614 with self.assertRaises(RuntimeError): 1615 d = x / c 1616 1617 @slow() 1618 @memory_heavy() 1619 def test_MX_shapes(self): 1620 self.message("MX unary operations") 1621 1622 #self.checkarray(DM(Sparsity.lower(4),1),DM(Sparsity.dense(4,4),1)) 1623 1624 for sp in [Sparsity.dense(0,0),Sparsity.dense(0,2),Sparsity.dense(2,0),Sparsity.dense(1,1),Sparsity.dense(2,2), Sparsity(4,3,[0,2,2,3],[1,2,1])]: 1625 for v in [0,1,0.2]: 1626 x_ = DM(sp,v) 1627 1628 xx = MX.sym("x",sp.size1(),sp.size2()) 1629 x=xx[sp] 1630 1631 for (casadiop, numpyop,name, flags) in self.pool.zip(): 1632 if 'nozero' in flags and v==0: continue 1633 r = casadiop([x]) 1634 f = Function("f", [xx],[r]) 1635 rr = f(v) 1636 self.checkarray(rr,numpyop(x_)) 1637 1638 a = DM(f.sparsity_out(0),1) 1639 b = DM.ones(DM(numpyop(x_)).sparsity()) 1640 1641 c = b-a 1642 if c.nnz()>0: 1643 # At least as sparse as DM calculus 1644 self.assertTrue(min(c.nonzeros())>=0,str([sp,v,name])) 1645 1646 for sp in [Sparsity(1,1),Sparsity.dense(1,1),Sparsity(3,4),Sparsity.dense(3,4), Sparsity(4,3,[0,2,2,3],[1,2,1]).T]: 1647 for v1 in [0,1,0.2,-0.2]: 1648 x1_ = DM(sp,v1) 1649 xx1 = MX.sym("x",sp.size1(),sp.size2()) 1650 x1=xx1[sp] 1651 xx1s = SX.sym("x",sp.size1(),sp.size2()) 1652 x1s=xx1s[sp] 1653 for sp2 in [Sparsity(1,1),Sparsity.dense(1,1),Sparsity(3,4),Sparsity.dense(3,4), Sparsity(4,3,[0,2,2,3],[1,2,1]).T]: 1654 for v2 in [0,1,0.2,-0.2]: 1655 x2_ = DM(sp2,v2) 1656 xx2 = MX.sym("x",sp2.size1(),sp2.size2()) 1657 x2=xx2[sp2] 1658 xx2s = SX.sym("x",sp2.size1(),sp2.size2()) 1659 x2s=xx2s[sp2] 1660 for (casadiop, numpyop,name, flags) in self.matrixbinarypool.zip(): 1661 if ("mul" in name or "mtimes" in name) and (sp.numel()==1 or sp2.numel()==1): continue 1662 r = casadiop([x1,x2]) 1663 f = Function("f", [xx1,xx2],[r]) 1664 f_in = [v1, v2] 1665 f_out = f(*f_in) 1666 g = Function("g", [xx1,xx2],[r]) 1667 g_out = g(v1, v2) 1668 1669 self.checkarray(f_out,numpyop([x1_,x2_]),str([sp,sp2,v1,v2,x1_,x2_,name])) 1670 1671 1672 if "mul" not in name: 1673 a = DM.ones(f.sparsity_out(0)) 1674 b = DM.ones(g.sparsity_out(0)) 1675 1676 c = b-a 1677 if c.nnz()>0: 1678 # At least as sparse as DM calculus 1679 self.assertTrue(min(c.nonzeros())>=0,str([sp,sp2,v1,v2,name,a,b])) 1680 1681 if sp.nnz()>0 and sp2.nnz()>0 and v1!=0 and v2!=0: 1682 self.checkfunction(f,g,inputs=f_in,hessian=False,failmessage=str([sp,sp2,v1,v2,x1_,x2_,name])) 1683 1684 @memory_heavy() 1685 def test_MXConstant(self): 1686 self.message("MX unary operations, constant") 1687 1688 #self.checkarray(DM(Sparsity.lower(4),1),DM(Sparsity.dense(4,4),1)) 1689 1690 for sp in [Sparsity.dense(0,0),Sparsity.dense(0,2),Sparsity.dense(2,0),Sparsity.dense(1,1),Sparsity.dense(2,2), Sparsity(4,3,[0,2,2,3],[1,2,1])]: 1691 for v in [0,1,0.2]: 1692 x_ = DM(sp,v) 1693 1694 x=MX(sp,v) 1695 1696 for (casadiop, numpyop,name, flags) in self.pool.zip(): 1697 if 'nozero' in flags and (v==0 or not sp.is_dense()): continue 1698 r = casadiop([x]) 1699 print(r) 1700 self.assertTrue(r.is_constant()) 1701 1702 self.checkarray(r.to_DM(),numpyop(x_),str([x_,name])) 1703 1704 a = DM.ones(r.to_DM().sparsity()) 1705 b = DM.ones(DM(numpyop(x_)).sparsity()) 1706 1707 c = b-a 1708 if c.nnz()>0: 1709 # At least as sparse as DM calculus 1710 self.assertTrue(min(c.nonzeros())>=0,str([sp,v,name])) 1711 1712 for sp in [Sparsity.dense(1,1),Sparsity(1,1),Sparsity(3,4),Sparsity.dense(3,4), Sparsity(4,3,[0,2,2,3],[1,2,1]).T]: 1713 for v1 in [0,1,0.2,-0.2]: 1714 x1_ = DM(sp,v1) 1715 x1=MX(sp,v1) 1716 for sp2 in [Sparsity.dense(1,1),Sparsity(1,1),Sparsity(3,4),Sparsity.dense(3,4), Sparsity(4,3,[0,2,2,3],[1,2,1]).T]: 1717 for v2 in [0,1,0.2,-0.2]: 1718 x2_ = DM(sp2,v2) 1719 x2=MX(sp2,v2) 1720 for (casadiop, numpyop,name, flags) in self.matrixbinarypool.zip(): 1721 if ("mul" in name or "mtimes" in name) and (sp.numel()==1 or sp2.numel()==1): continue 1722 r = casadiop([x1,x2]) 1723 f = Function("f", [],[r]) # Should not be needed -> constant folding 1724 f_out = f.call([]) 1725 1726 1727 self.checkarray(f_out[0],numpyop([x1_,x2_]),str([sp,sp2,v1,v2,name])) 1728 if "mul" not in name and "mtimes" not in name: 1729 a = DM.ones(f.sparsity_out(0)) 1730 b = DM.ones(DM(numpyop([x1_,x2_])).sparsity()) 1731 1732 c = b-a 1733 if c.nnz()>0: 1734 # At least as sparse as DM calculus 1735 self.assertTrue(min(c.nonzeros())>=0,str([sp,sp2,v1,v2,name])) 1736 1737 def test_graph_substitute(self): 1738 x=MX.sym("X",4,4) 1739 y=MX.sym("Y",4,4) 1740 b=MX.sym("B",4,4) 1741 1742 c = x*y 1743 d = b*y 1744 f = c+d 1745 1746 1747 C = MX.sym("C",4,4) 1748 f = graph_substitute(f,[c],[C]) 1749 1750 F = Function("F", [y,b,C],[f]) 1751 F_out = F(1, 2, 3) 1752 1753 self.checkarray(F_out,5*DM.ones(4,4)) 1754 1755 D = MX.sym("D",4,4) 1756 f = graph_substitute(f,[d],[D]) 1757 1758 F = Function("F", [D,C],[f]) 1759 F_out = F(4, 5) 1760 1761 self.checkarray(F_out,9*DM.ones(4,4)) 1762 1763 1764 def test_matrix_expand(self): 1765 n = 2 1766 a = MX.sym("a",n,n) 1767 b = MX.sym("b",n,n) 1768 c = MX.sym("c",n,n) 1769 1770 d = a+b 1771 e = d*c 1772 1773 self.assertEqual(n_nodes(e),6) 1774 1775 t0 = matrix_expand(e) 1776 1777 self.assertEqual(n_nodes(t0),5) 1778 1779 t1 = matrix_expand(e,[d]) 1780 self.assertEqual(n_nodes(t1),6) 1781 1782 print(e,t0,t1) 1783 1784 1785 outs = [] 1786 for x in [e,t0,t1]: 1787 f = Function("f", [a,b,c],[x]) 1788 1789 f_in = [1.1, 2.2, 3.3] 1790 f_out = f(*f_in) 1791 1792 outs.append(f_out) 1793 if len(outs)>1: 1794 self.checkarray(outs[0],outs[-1]) 1795 1796 print(outs) 1797 1798 def test_kron(self): 1799 a = sparsify(DM([[1,0,6],[2,7,0]])) 1800 b = sparsify(DM([[1,0,0],[2,3,7],[0,0,9],[1,12,13]])) 1801 1802 A = MX.sym("A",a.sparsity()) 1803 B = MX.sym("B",b.sparsity()) 1804 C = c.kron(A,B) 1805 1806 f = Function("f", [A,B],[C]) 1807 f_in = [a, b] 1808 f_out = f(*f_in) 1809 1810 c_ = f_out 1811 1812 self.assertEqual(c_.size1(),a.size1()*b.size1()) 1813 self.assertEqual(c_.size2(),a.size2()*b.size2()) 1814 self.assertEqual(c_.nnz(),a.nnz()*b.nnz()) 1815 1816 self.checkarray(c_,numpy.kron(a,b)) 1817 1818 def test_project(self): 1819 x = MX.sym("x",Sparsity.lower(3)) 1820 y = project(x, Sparsity.lower(3).T) 1821 1822 f = Function("f", [x],[y]) 1823 f_in=DM(f.sparsity_in(0),list(range(1,int(4*3/2+1)))) 1824 f_out = f(f_in) 1825 1826 self.checkarray(f_out,DM([[1,0,0],[0,4,0],[0,0,6]])) 1827 self.checkarray(DM.ones(f.sparsity_out(0)),DM.ones(Sparsity.lower(3).T)) 1828 1829 self.check_codegen(f,inputs=[DM([[1,0,0],[0,4,0],[0,0,6]])]) 1830 self.check_serialize(f,inputs=[DM([[1,0,0],[0,4,0],[0,0,6]])]) 1831 1832 def test_repmat(self): 1833 a = DM([[1,2],[3,4],[5,6]]) 1834 self.checkarray(repmat(a,2,3),kron(DM.ones(2,3),a)) 1835 1836 def test_transpose_bug(self): 1837 1838 A = [[-26.9091,00,00,1,00,00,00,00,00,00,00,00,00,00,00], 1839 [00,-26.9091,00,00,1,00,00,00,00,00,00,00,00,00,00], 1840 [00,00,-26.9091,00,00,1,00,00,00,00,00,00,00,00,00], 1841 [00,00,00,-14.1393,-0.348109,-0.382033,3.6491,9.40422,-5.05449,3.03347,-0.00126949,-0.000414104,00,0.00498232,0.00495617], 1842 [00,00,00,-0.348109,-13.6315,-0.194202,-6.41868,-0.0189287,3.34725,1.86987,-0.000645326,-0.000210504,00,0.0025327,0.0025194], 1843 [00,00,00,-0.382033,-0.194202,-13.6677,10.1342,-13.0101,-0.890078,-8.99621,-0.000708215,-0.000231018,00,0.00277951,0.00276493], 1844 [00,00,00,00,00,00,-27.0487,0.36775,-0.277186,0.564995,0.424389,0.0235157,0.227774,00,00], 1845 [00,00,00,00,00,00,0.113776,-27.3241,0.159254,-0.520988,-0.0235157,0.424389,0.132128,00,00], 1846 [00,00,00,00,00,00,0.227472,-0.0735535,-26.9135,0.206825,-0.227774,-0.132128,0.424389,00,00], 1847 [00,00,00,00,00,00,0.332188,-1.02565,-0.0471482,-28.3499,0.132128,-0.227774,0.0235157,00,00], 1848 [00,00,00,-0.00126949,-0.000645326,-0.000708215,0.00763923,0.0016513,-0.0044564,-0.00284703,-0.119342,-0.000941301,-0.00160961,-1.62792e-05,0.00156863], 1849 [00,00,00,-0.000414104,-0.000210504,-0.000231018,0.0024919,0.000538651,-0.00145367,-0.000928698,0.000941301,-0.119493,0.000742542,-0.00155303,-7.50989e-06], 1850 [00,00,00,00,00,00,00,00,00,00,2.1684e-19,-2.1684e-19,-0.171654,-0.000868032,0.000868032], 1851 [00,00,00,0.00181259,0.000921407,0.0010112,-0.0109074,-0.00235775,0.00636292,0.00406504,0,-0.000567731,-0.000868032,-0.00087529,00], 1852 [00,00,00,0.00166876,0.000848291,0.000930959,-0.0100419,-0.00217066,0.005858,0.00374247,0.00052268,0,-0.000868032,00,-0.000874062] 1853 ] 1854 1855 A = sparsify(DM(A)) 1856 1857 As = MX.sym("As",A.sparsity()) 1858 1859 f = Function("f", [As],[densify(As.T),densify(As).T,As.T,As,densify(As)]) 1860 1861 f_in = [0]*f.n_in() 1862 f_in[0]=A 1863 f_out = f(*f_in) 1864 1865 self.checkarray(f_out[0],A.T) 1866 self.checkarray(f_out[1],A.T) 1867 self.checkarray(f_out[2],A.T) 1868 self.checkarray(f_out[3],A) 1869 self.checkarray(f_out[4],A) 1870 1871 @requiresPlugin(Linsol,"csparse") 1872 def test_bizarre_bug(self): 1873 1874 A = [[-26.9091,00,00,1,00,00,00,00,00,00,00,00,00,00,00], 1875 [00,-26.9091,00,00,1,00,00,00,00,00,00,00,00,00,00], 1876 [00,00,-26.9091,00,00,1,00,00,00,00,00,00,00,00,00], 1877 [00,00,00,-14.1393,-0.348109,-0.382033,3.6491,9.40422,-5.05449,3.03347,-0.00126949,-0.000414104,00,0.00498232,0.00495617], 1878 [00,00,00,-0.348109,-13.6315,-0.194202,-6.41868,-0.0189287,3.34725,1.86987,-0.000645326,-0.000210504,00,0.0025327,0.0025194], 1879 [00,00,00,-0.382033,-0.194202,-13.6677,10.1342,-13.0101,-0.890078,-8.99621,-0.000708215,-0.000231018,00,0.00277951,0.00276493], 1880 [00,00,00,00,00,00,-27.0487,0.36775,-0.277186,0.564995,0.424389,0.0235157,0.227774,00,00], 1881 [00,00,00,00,00,00,0.113776,-27.3241,0.159254,-0.520988,-0.0235157,0.424389,0.132128,00,00], 1882 [00,00,00,00,00,00,0.227472,-0.0735535,-26.9135,0.206825,-0.227774,-0.132128,0.424389,00,00], 1883 [00,00,00,00,00,00,0.332188,-1.02565,-0.0471482,-28.3499,0.132128,-0.227774,0.0235157,00,00], 1884 [00,00,00,-0.00126949,-0.000645326,-0.000708215,0.00763923,0.0016513,-0.0044564,-0.00284703,-0.119342,-0.000941301,-0.00160961,-1.62792e-05,0.00156863], 1885 [00,00,00,-0.000414104,-0.000210504,-0.000231018,0.0024919,0.000538651,-0.00145367,-0.000928698,0.000941301,-0.119493,0.000742542,-0.00155303,-7.50989e-06], 1886 [00,00,00,00,00,00,00,00,00,00,2.1684e-19,-2.1684e-19,-0.171654,-0.000868032,0.000868032], 1887 [00,00,00,0.00181259,0.000921407,0.0010112,-0.0109074,-0.00235775,0.00636292,0.00406504,0,-0.000567731,-0.000868032,-0.00087529,00], 1888 [00,00,00,0.00166876,0.000848291,0.000930959,-0.0100419,-0.00217066,0.005858,0.00374247,0.00052268,0,-0.000868032,00,-0.000874062] 1889 ] 1890 1891 A = sparsify(DM(A)) 1892 1893 b = DM(list(range(15))) 1894 H = 5 1895 1896 Bs = MX.sym("Bs",b.sparsity()) 1897 As = MX.sym("As",A.sparsity()) 1898 1899 Ast = As.T 1900 1901 r= Function("r", [As,Bs],[solve(Ast,Bs,"csparse")]) 1902 R= Function("R", [As,Bs],[solve(densify(Ast),Bs,"csparse")]) 1903 1904 r_in = (A, b) 1905 r_out = r(*r_in) 1906 R_out = R(*r_in) 1907 1908 r.sparsity_out(0).spy() 1909 R.sparsity_out(0).spy() 1910 1911 self.checkarray(R_out,numpy.linalg.solve(A.T,b)) 1912 self.checkarray(r_out,R_out) 1913 1914 def test_depends_on(self): 1915 a = MX.sym("a") 1916 b = MX.sym("b") 1917 1918 self.assertTrue(depends_on(a**2,a)) 1919 self.assertTrue(depends_on(a,a)) 1920 self.assertFalse(depends_on(0,a)) 1921 1922 ab = vertcat(*(a,b)) 1923 self.assertTrue(depends_on(a**2,ab)) 1924 self.assertTrue(depends_on(a,ab)) 1925 self.assertFalse(depends_on(0,ab)) 1926 self.assertTrue(depends_on(b**2,ab)) 1927 self.assertTrue(depends_on(b,ab)) 1928 self.assertTrue(depends_on(a**2+b**2,ab)) 1929 self.assertTrue(depends_on(a+b,ab)) 1930 self.assertTrue(depends_on(vertcat(*[0,a]),a)) 1931 self.assertTrue(depends_on(vertcat(*[a,0]),a)) 1932 self.assertFalse(depends_on(vertcat(*[0,b]),a)) 1933 self.assertFalse(depends_on(vertcat(*[b,0]),a)) 1934 self.assertTrue(depends_on(vertcat(*[a**2,b**2]),ab)) 1935 self.assertTrue(depends_on(vertcat(*[a,0]),ab)) 1936 self.assertTrue(depends_on(vertcat(*[0,b]),ab)) 1937 self.assertTrue(depends_on(vertcat(*[b,0]),ab)) 1938 self.assertFalse(depends_on(vertcat(*[0,0]),ab)) 1939 1940 def test_vertcat_simp(self): 1941 x = MX.sym("x",10) 1942 y = MX.sym("y") 1943 z = MX.sym("z") 1944 x_ = DM(list(range(10))) 1945 y_ = DM([20]) 1946 z_ = DM([30]) 1947 1948 def evalvertcat(*a): 1949 f = Function("f", [x,y,z],[vertcat(*a)]) 1950 f_in = [0]*f.n_in();f_in[0]=x_ 1951 f_in[1]=y_ 1952 f_in[2]=z_ 1953 return f(*f_in) 1954 1955 self.checkarray(evalvertcat(*vertsplit(x)),x_) 1956 self.checkarray(evalvertcat(*vertsplit(x)+[y]),vertcat(*[x_,y_])) 1957 self.checkarray(evalvertcat(*[z]+vertsplit(x)+[y] + vertsplit(x)+[z]),vertcat(*[z_,x_,y_,x_,z_])) 1958 self.checkarray(evalvertcat(*vertsplit(x)[:-1]),x_[:-1]) 1959 self.checkarray(evalvertcat(*vertsplit(x)[:-1]+[y]),vertcat(*[x_[:-1],y_])) 1960 self.checkarray(evalvertcat(*[z]+vertsplit(x)[:-1]+[y] + vertsplit(x)[:-1]+[z]),vertcat(*[z_,x_[:-1],y_,x_[:-1],z_])) 1961 self.checkarray(evalvertcat(*vertsplit(x)[1:]),x_[1:]) 1962 self.checkarray(evalvertcat(*vertsplit(x)[1:]+[y]),vertcat(*[x_[1:],y_])) 1963 self.checkarray(evalvertcat(*[z]+vertsplit(x)[1:]+[y] + vertsplit(x)[1:]+[z]),vertcat(*[z_,x_[1:],y_,x_[1:],z_])) 1964 g = vertsplit(x)[5:]+vertsplit(x)[:5] 1965 self.checkarray(evalvertcat(*g),vertcat(*[x_[5:],x_[:5]])) 1966 self.checkarray(evalvertcat(*g+[y]),vertcat(*[x_[5:],x_[:5],y_])) 1967 self.checkarray(evalvertcat(*[z]+g+[y] + g+[z]),vertcat(*[z_,x_[5:],x_[:5],y_,x_[5:],x_[:5],z_])) 1968 1969 1970 1971 w = vertsplit(x,2) 1972 r = builtins.sum([vertsplit(i) for i in w],[]) 1973 1974 self.checkarray(evalvertcat(*r),x_) 1975 1976 w = vertsplit(x,2) 1977 r = builtins.sum([vertsplit(i)+[y] for i in w],[]) 1978 print("vertcat:", r) 1979 print("result:", vertcat(*r)) 1980 1981 w = vertsplit(x,2) 1982 r = builtins.sum([vertsplit(i) for i in w],[]) 1983 print("vertcat:", r) 1984 print("result:", vertcat(*r+[y])) 1985 1986 self.assertTrue(is_equal(vertcat(*vertsplit(x)),x)) 1987 1988 def test_horzcat_simp(self): 1989 x = MX.sym("x",1,10) 1990 y = MX.sym("y") 1991 z = MX.sym("z") 1992 x_ = DM(list(range(10))).T 1993 y_ = DM([20]) 1994 z_ = DM([30]) 1995 1996 def evalhorzcat(*a): 1997 f = Function("f", [x,y,z],[horzcat(*a)]) 1998 return f(x_, y_, z_) 1999 2000 self.checkarray(evalhorzcat(*horzsplit(x)),x_) 2001 self.checkarray(evalhorzcat(*horzsplit(x)+[y]),horzcat(*[x_,y_])) 2002 self.checkarray(evalhorzcat(*[z]+horzsplit(x)+[y] + horzsplit(x)+[z]),horzcat(*[z_,x_,y_,x_,z_])) 2003 self.checkarray(evalhorzcat(*horzsplit(x)[:-1]),x_[0,:-1]) 2004 self.checkarray(evalhorzcat(*horzsplit(x)[:-1]+[y]),horzcat(*[x_[0,:-1],y_])) 2005 self.checkarray(evalhorzcat(*[z]+horzsplit(x)[:-1]+[y] + horzsplit(x)[:-1]+[z]),horzcat(*[z_,x_[0,:-1],y_,x_[0,:-1],z_])) 2006 self.checkarray(evalhorzcat(*horzsplit(x)[1:]),x_[0,1:]) 2007 self.checkarray(evalhorzcat(*horzsplit(x)[1:]+[y]),horzcat(*[x_[0,1:],y_])) 2008 self.checkarray(evalhorzcat(*[z]+horzsplit(x)[1:]+[y] + horzsplit(x)[1:]+[z]),horzcat(*[z_,x_[0,1:],y_,x_[0,1:],z_])) 2009 g = horzsplit(x)[5:]+horzsplit(x)[:5] 2010 self.checkarray(evalhorzcat(*g),horzcat(*[x_[0,5:],x_[0,:5]])) 2011 self.checkarray(evalhorzcat(*g+[y]),horzcat(*[x_[0,5:],x_[0,:5],y_])) 2012 self.checkarray(evalhorzcat(*[z]+g+[y] + g+[z]),horzcat(*[z_,x_[0,5:],x_[0,:5],y_,x_[0,5:],x_[0,:5],z_])) 2013 2014 2015 w = horzsplit(x,2) 2016 r = builtins.sum([horzsplit(i) for i in w],[]) 2017 2018 self.checkarray(evalhorzcat(*r),x_) 2019 2020 w = horzsplit(x,2) 2021 r = builtins.sum([horzsplit(i)+[y] for i in w],[]) 2022 print("vertcat:", r) 2023 print("result:", horzcat(*r)) 2024 2025 w = horzsplit(x,2) 2026 r = builtins.sum([horzsplit(i) for i in w],[]) 2027 print("vertcat:", r) 2028 print("result:", horzcat(*r+[y])) 2029 2030 self.assertTrue(is_equal(horzcat(*horzsplit(x)),x)) 2031 2032 def test_vertsplit_simp(self): 2033 2034 dvars = [MX.sym("abcdefghijklm"[i]) for i in range(5) ] 2035 dvars_ = list(range(5)) 2036 2037 zz = MX.sym("zz",2) 2038 zz_ = DM([11,12]) 2039 y = MX.sym("y") 2040 z = MX.sym("z") 2041 y_ = DM([20]) 2042 z_ = DM([30]) 2043 2044 aa = MX.sym("aa",5) 2045 aa_ = list(range(100,105)) 2046 2047 def evalvertsplit(a,*args): 2048 print(vertsplit(a,*args)) 2049 f = Function("f", dvars+[y,z,zz,aa],vertsplit(a,*args)) 2050 f_in = [0]*f.n_in() 2051 for i in range(5): 2052 f_in[i]=dvars_[i] 2053 f_in[5+0]=y_ 2054 f_in[5+1]=z_ 2055 f_in[5+2]=zz_ 2056 f_in[5+3]=aa_ 2057 f_out = f.call(f_in) 2058 return [f_out[i] for i in range(f.n_out())] 2059 2060 s= evalvertsplit(vertcat(*[y]+dvars+[z])) 2061 self.checkarray(s[0],y_) 2062 for i in range(5): 2063 self.checkarray(s[1+i],dvars_[i]) 2064 self.checkarray(s[6],z_) 2065 2066 s= evalvertsplit(vertcat(*[y]+dvars+[z]),2) 2067 2068 self.checkarray(s[0],vertcat(*[y_,dvars_[0]])) 2069 self.checkarray(s[1],vertcat(*[dvars_[1],dvars_[2]])) 2070 self.checkarray(s[2],vertcat(*[dvars_[3],dvars_[4]])) 2071 self.checkarray(s[3],vertcat(*[z_])) 2072 2073 s= evalvertsplit(vertcat(*[y,zz,z,zz]),2) 2074 2075 self.checkarray(s[0],vertcat(*[y_,zz_[0]])) 2076 self.checkarray(s[1],vertcat(*[zz_[1],z_])) 2077 self.checkarray(s[2],zz_) 2078 2079 s= evalvertsplit(vertcat(*[y,zz,z,zz]),3) 2080 2081 self.checkarray(s[0],vertcat(*[y_,zz_[0],zz_[1]])) 2082 self.checkarray(s[1],vertcat(*[z_,zz_[0],zz_[1]])) 2083 2084 s= evalvertsplit(vertcat(*[zz,zz]),2) 2085 self.checkarray(s[0],zz_) 2086 self.checkarray(s[1],zz_) 2087 2088 s= evalvertsplit(vertcat(*[zz]+dvars)) 2089 self.checkarray(s[0],zz_[0]) 2090 self.checkarray(s[1],zz_[1]) 2091 2092 for i in range(5): 2093 self.checkarray(s[2+i],dvars_[i]) 2094 2095 s= evalvertsplit(vertcat(*dvars+[aa]),5) 2096 self.checkarray(s[0],DM(dvars_)) 2097 self.checkarray(s[1],DM(aa_)) 2098 2099 s= evalvertsplit(vertcat(*dvars+[aa]),4) 2100 self.checkarray(s[0],DM(dvars_[:4])) 2101 self.checkarray(s[1],DM([dvars_[-1]]+aa_[:3])) 2102 self.checkarray(s[2],DM(aa_[3:])) 2103 2104 s= evalvertsplit(vertcat(*dvars+[aa]),6) 2105 self.checkarray(s[0],DM(dvars_+[aa_[0]])) 2106 self.checkarray(s[1],DM(aa_[1:])) 2107 2108 for i in range(5): 2109 self.assertTrue(is_equal(vertsplit(vertcat(*dvars))[i],dvars[i])) 2110 2111 def test_horzsplit_simp(self): 2112 2113 dvars = [MX.sym("abcdefghijklm"[i]) for i in range(5) ] 2114 dvars_ = list(range(5)) 2115 2116 zz = MX.sym("zz",1,2) 2117 zz_ = DM([11,12]).T 2118 y = MX.sym("y") 2119 z = MX.sym("z") 2120 y_ = DM([20]) 2121 z_ = DM([30]) 2122 2123 aa = MX.sym("aa",1,5) 2124 aa_ = list(range(100,105)) 2125 2126 def evalhorzsplit(a,*args): 2127 print(horzsplit(a,*args)) 2128 f = Function("f", dvars+[y,z,zz,aa],horzsplit(a,*args)) 2129 f_in = [0]*f.n_in() 2130 for i in range(5): 2131 f_in[i]=dvars_[i] 2132 f_in[5+0]=y_ 2133 f_in[5+1]=z_ 2134 f_in[5+2]=zz_ 2135 f_in[5+3]=aa_ 2136 f_out = f(*f_in) 2137 return [f_out[i] for i in range(f.n_out())] 2138 2139 s= evalhorzsplit(horzcat(*[y]+dvars+[z])) 2140 self.checkarray(s[0],y_) 2141 for i in range(5): 2142 self.checkarray(s[1+i],dvars_[i]) 2143 self.checkarray(s[6],z_) 2144 2145 s= evalhorzsplit(horzcat(*[y]+dvars+[z]),2) 2146 2147 self.checkarray(s[0],vertcat(*[y_,dvars_[0]]).T) 2148 self.checkarray(s[1],vertcat(*[dvars_[1],dvars_[2]]).T) 2149 self.checkarray(s[2],vertcat(*[dvars_[3],dvars_[4]]).T) 2150 self.checkarray(s[3],vertcat(*[z_]).T) 2151 2152 s= evalhorzsplit(horzcat(*[y,zz,z,zz]),2) 2153 2154 self.checkarray(s[0],vertcat(*[y_,zz_[0,0]]).T) 2155 self.checkarray(s[1],vertcat(*[zz_[0,1],z_]).T) 2156 self.checkarray(s[2],zz_) 2157 2158 s= evalhorzsplit(horzcat(*[y,zz,z,zz]),3) 2159 2160 self.checkarray(s[0],vertcat(*[y_,zz_[0,0],zz_[0,1]]).T) 2161 self.checkarray(s[1],vertcat(*[z_,zz_[0,0],zz_[0,1]]).T) 2162 2163 s= evalhorzsplit(horzcat(*[zz,zz]),2) 2164 self.checkarray(s[0],zz_) 2165 self.checkarray(s[1],zz_) 2166 2167 s= evalhorzsplit(horzcat(*[zz]+dvars)) 2168 self.checkarray(s[0],zz_[0,0]) 2169 self.checkarray(s[1],zz_[0,1]) 2170 2171 for i in range(5): 2172 self.checkarray(s[2+i],dvars_[i]) 2173 2174 s= evalhorzsplit(horzcat(*dvars+[aa]),5) 2175 self.checkarray(s[0],DM(dvars_).T) 2176 self.checkarray(s[1],DM(aa_).T) 2177 2178 s= evalhorzsplit(horzcat(*dvars+[aa]),4) 2179 self.checkarray(s[0],DM(dvars_[:4]).T) 2180 self.checkarray(s[1],DM([dvars_[-1]]+aa_[:3]).T) 2181 self.checkarray(s[2],DM(aa_[3:]).T) 2182 2183 s= evalhorzsplit(horzcat(*dvars+[aa]),6) 2184 self.checkarray(s[0],DM(dvars_+[aa_[0]]).T) 2185 self.checkarray(s[1],DM(aa_[1:]).T) 2186 2187 for i in range(5): 2188 self.assertTrue(is_equal(horzsplit(horzcat(*dvars))[i],dvars[i])) 2189 2190 def test_vertsplit_derivative(self): 2191 m = MX.sym("X",10) 2192 2193 f = Function("f", [m],[vertsplit(m)[0]]) 2194 2195 f.reverse(1) 2196 2197 def test_MX_const_sp(self): 2198 x = MX.sym("x",4,1) 2199 2200 sp = Sparsity.triplet(3,3,[0,1,2,2],[0,0,1,2]) 2201 2202 f = Function("f", [x],[x.nz[DM(sp,list(range(sp.nnz())))]]) 2203 2204 g = Function("g", [x],[MX(sp,x)]) 2205 2206 f_in = [0]*f.n_in();f_in[0]=DM(list(range(1,5))) 2207 2208 self.checkfunction(f,g,inputs=f_in) 2209 2210 def test_reshape_sp(self): 2211 x = MX.sym("x",4,1) 2212 2213 f = Function("f", [x],[x.reshape((2,2))]) 2214 2215 sx = SX.sym("x",4,1) 2216 2217 g = Function("g", [sx],[sx.reshape((2,2))]) 2218 2219 f_in = [0]*f.n_in();f_in[0]=DM(list(range(1,5))) 2220 2221 self.checkfunction(f,g,inputs=f_in) 2222 2223 def test_issue1041(self): 2224 x = MX.sym("x",2) 2225 2226 y = vertsplit(x,[0,1,2])[1] 2227 2228 f = Function("f", [x],[y]) 2229 2230 H = f.hessian_old(0, 0) 2231 2232 def test_bug_1042(self): 2233 2234 x = MX.sym('x',2,1) 2235 2236 mf = Function("mf", [x],[x*x[0,0]]) 2237 2238 mfunction = mf.expand('expand_'+mf.name()) 2239 2240 mfg = mf.reverse(1) 2241 2242 mfunctiong = mfunction.reverse(1) 2243 2244 f_in = [0, 5, DM([1,2])] 2245 2246 self.checkfunction(mfg,mfunctiong,inputs=f_in) 2247 2248 def test_bug_1042bis(self): 2249 x = MX.sym('x',2,1) 2250 a = MX.sym("ax",2,1) 2251 i1 = x[0,0] 2252 z = i1*x 2253 i3 = i1*a 2254 i3= c.dot(x,a) 2255 d = Function("d", [x,a],[z,i3]) 2256 2257 dx = d.expand('expand_'+d.name()) 2258 dx_in = [0]*dx.n_in();dx_in[0]=DM([1,2]) 2259 dx_in[1]=DM([3,4]) 2260 2261 self.checkfunction(d,dx,inputs=dx_in) 2262 2263 def test_bug_1042tris(self): 2264 x = MX.sym('x',2,1) 2265 a = MX.sym("ax",2,1) 2266 d = Function("d", [x,a],[c.dot(x,a)]) 2267 2268 dx = d.expand('expand_'+d.name()) 2269 dx_in = [0]*dx.n_in();dx_in[0]=DM([1,2]) 2270 dx_in[1]=DM([3,4]) 2271 2272 self.checkfunction(d,dx,inputs=dx_in) 2273 2274 def test_bug_1046(self): 2275 x = MX.sym('x',1,1) 2276 y = MX.sym('y',1,1) 2277 z = jacobian(x,y) 2278 2279 self.assertTrue(z.nnz()==0) 2280 2281 def test_expand_free(self): 2282 x = MX.sym("x") 2283 y = MX.sym("y") 2284 f = Function('f',[x],[x+y]) 2285 with self.assertInException("free"): 2286 f.expand() 2287 2288 g = Function('g',[x,y],[f(x)]) 2289 ge = g.expand() 2290 self.checkfunction(g,g,inputs=[1,2]) 2291 2292 def test_singularcat(self): 2293 2294 for c in [MX,SX,DM]: 2295 x0 = c.zeros(10,0) 2296 2297 x1s = vertsplit(x0, [0,5,10]) 2298 2299 for x in x1s: 2300 self.checkarray(x.shape,(5,0)) 2301 2302 2303 x2 = vertcat(*x1s) 2304 self.checkarray(x2.shape,(10,0)) 2305 2306 x2 = vertcat(*[c.zeros(0,0)] + x1s + [c.zeros(0,0)]) 2307 self.checkarray(x2.shape,(10,0)) 2308 2309 x0 = c.zeros(0,1) 2310 x2 = vertcat(x0,c.zeros(0,0),x0) 2311 self.checkarray(x2.shape,(0,1)) 2312 2313 for c in [MX,SX,DM]: 2314 x0 = c.zeros(0,10) 2315 2316 x1s = horzsplit(x0, [0,5,10]) 2317 2318 for x in x1s: 2319 self.checkarray(x.shape,(0,5)) 2320 2321 x2 = horzcat(*x1s) 2322 self.checkarray(x2.shape,(0,10)) 2323 2324 x2 = horzcat(*[c.zeros(0,0)] + x1s + [c.zeros(0,0)]) 2325 self.checkarray(x2.shape,(0,10)) 2326 2327 x0 = c.zeros(1,0) 2328 x2 = horzcat(x0,c.zeros(0,0),x0) 2329 self.checkarray(x2.shape,(1,0)) 2330 2331 for c in [MX,SX,DM]: 2332 x0 = c.zeros(10,0) 2333 2334 x1s = vertsplit(x0, [0,5,10]) 2335 2336 x0 = c.zeros(0,10) 2337 x1st = horzsplit(x0, [0,5,10]) 2338 2339 x2 = diagcat(*x1s) 2340 self.checkarray(x2.shape,(10,0)) 2341 2342 x2 = diagcat(*[c.zeros(0,0)] + x1s + [c.zeros(0,0)]) 2343 self.checkarray(x2.shape,(10,0)) 2344 2345 x2 = diagcat(*x1st) 2346 self.checkarray(x2.shape,(0,10)) 2347 2348 x2 = diagcat(*[c.zeros(0,0)] + x1st + [c.zeros(0,0)]) 2349 self.checkarray(x2.shape,(0,10)) 2350 2351 x2 = diagcat(*x1s+x1st) 2352 self.checkarray(x2.shape,(10,10)) 2353 2354 x2 = diagcat(*[c.zeros(0,0)] + x1s+x1st + [c.zeros(0,0)]) 2355 self.checkarray(x2.shape,(10,10)) 2356 2357 def test_exprjac(self): 2358 2359 2360 for fun in [lambda x: x[0]*x[1],lambda x: x[0]*sin(x[1])]: 2361 def hessian1(f, x): return hessian(f, x)[0] 2362 for op in [c.gradient, jacobian, hessian1]: 2363 print(fun, op) 2364 x = MX.sym("x",2) 2365 print(fun(x)) 2366 print(op(fun(x),x)) 2367 f = Function("f", [x],[op(fun(x),x)]) 2368 2369 x = SX.sym("x",2) 2370 fr = Function("fr", [x],[op(fun(x),x)]) 2371 2372 self.checkfunction(f,fr,inputs=[0]) 2373 2374 @memory_heavy() 2375 def test_einstein(self): 2376 dim_b = [3, 3] 2377 dim_a = [3, 3] 2378 2379 2380 2381 2382 def free(a): 2383 return set([i for i in a if i<0]) 2384 2385 2386 def dim_map(dim_a, dim_b, dim_c, ind_a, ind_b, ind_c): 2387 dim_map = {} 2388 for dim, ind in [(dim_a,ind_a), (dim_b,ind_b)]: 2389 for i,j in enumerate(ind): 2390 if j<0: 2391 if j in dim_map: 2392 assert dim_map[j]==dim[i] 2393 else: 2394 dim_map[j]=dim[i] 2395 return dim_map 2396 2397 def sub2ind(dims, sub): 2398 ret = [] 2399 for i in range(len(dims)): 2400 ret.append(sub % dims[i]) 2401 sub/= dims[i] 2402 return ret 2403 2404 def ind2sub(dims, ind): 2405 ret=0 2406 cumprod = 1 2407 for i in range(len(dims)): 2408 ret+= ind[i]*cumprod 2409 cumprod*= dims[i] 2410 return ret 2411 2412 def combine(dim_a, dim_b, ind_a, ind_b, ind_c): 2413 dim_map = {} 2414 for dim, ind in [(dim_a,ind_a), (dim_b,ind_b)]: 2415 for i,j in enumerate(ind): 2416 if j<0: 2417 if j in dim_map: 2418 assert dim_map[j]==dim[i] 2419 else: 2420 dim_map[j]=dim[i] 2421 return [dim_map[i] for i in ind_c] 2422 2423 2424 def subs(e,old,n): 2425 return [ n[old.index(i)] if i in old else i for i in e ] 2426 2427 2428 def einstein_tests(dim_a, dim_b, dim_c, ind_a, ind_b, ind_c): 2429 2430 2431 A = MX.sym("A", np.product(dim_a)) 2432 B = MX.sym("B", np.product(dim_b)) 2433 C = MX.sym("C", int(np.product(dim_c)),1) 2434 2435 A_sx = SX.sym("A", np.product(dim_a)) 2436 B_sx = SX.sym("B", np.product(dim_b)) 2437 C_sx = SX.sym("C", int(np.product(dim_c)),1) 2438 2439 np.random.seed(0) 2440 2441 A_ = np.random.random(A.shape) 2442 B_ = np.random.random(B.shape) 2443 C_ = np.random.random(C.shape) 2444 2445 2446 2447 out = casadi.einstein(A, B, C, dim_a, dim_b, dim_c, ind_a, ind_b, ind_c) 2448 out_sx = casadi.einstein(A_sx, B_sx, C_sx, dim_a, dim_b, dim_c, ind_a, ind_b, ind_c) 2449 f = Function('f',[A,B,C],[out],{"ad_weight_sp": 0}) 2450 frev = Function('f',[A,B,C],[out],{"ad_weight_sp": 1}) 2451 fr = Function('fr',[A,B,C],[my_einstein(A, B, C, dim_a, dim_b, dim_c, ind_a, ind_b, ind_c)]) 2452 f_sx = Function('f',[A_sx,B_sx,C_sx],[out_sx]) 2453 2454 fsx = f.expand() 2455 self.checkfunction(f, fr, inputs=[A_,B_,C_]) 2456 self.checkfunction(fsx, fr, inputs=[A_,B_,C_]) 2457 self.check_codegen(f, inputs=[A_,B_,C_]) 2458 self.checkfunction(fsx, f_sx, inputs=[A_,B_,C_]) 2459 2460 for i in range(3): 2461 self.check_sparsity(f.sparsity_jac(i, 0), fsx.sparsity_jac(i, 0)) 2462 self.check_sparsity(frev.sparsity_jac(i, 0), fsx.sparsity_jac(i, 0)) 2463 2464 2465 def my_einstein(A, B, C, dim_a, dim_b, dim_c, ind_a, ind_b, ind_c): 2466 try: 2467 R = DM(C) 2468 except: 2469 R = MX(C) 2470 2471 d = dim_map(dim_a, dim_b, dim_c, ind_a, ind_b, ind_c) 2472 2473 dk = list(d.keys()) 2474 2475 for val in itertools.product(*[range(d[k]) for k in dk]): 2476 ind_a_ = subs(ind_a, dk, val) 2477 ind_b_ = subs(ind_b, dk, val) 2478 ind_c_ = subs(ind_c, dk, val) 2479 2480 2481 ai = ind2sub(dim_a, ind_a_) 2482 bi = ind2sub(dim_b, ind_b_) 2483 ci = ind2sub(dim_c, ind_c_) 2484 2485 R[ci]+= A[ai]*B[bi] 2486 2487 return R 2488 2489 ind = [ [0, -1], [1, -1], [0, 1], [-1, -2], [-2, -1] ] 2490 2491 for ind_a in ind: 2492 for ind_b in ind: 2493 for ind_c in itertools.permutations(list(free(ind_a) ^ free(ind_b))): 2494 dim_c = combine(dim_a, dim_b, ind_a, ind_b, ind_c) 2495 einstein_tests(dim_a, dim_b, dim_c, ind_a, ind_b, ind_c) 2496 2497 einstein_tests([2,4,3], [2,5,3], [5, 4], [-1, -2, -3], [-1, -4, -3], [-4, -2]) 2498 2499 def test_sparsity_operation(self): 2500 L = [MX(Sparsity(1,1)),MX(Sparsity(2,1)), MX.sym("x",1,1), MX.sym("x", Sparsity(1,1)), DM(1), DM(Sparsity(1,1),1), DM(Sparsity(2,1),1), DM(Sparsity.dense(2,1),1)] 2501 2502 for a in L: 2503 for b in L: 2504 c = a*b 2505 2506 if a.nnz()==0 or b.nnz()==0: 2507 self.assertTrue(c.nnz()==0) 2508 else: 2509 self.assertTrue(c.nnz()>0) 2510 2511 @requires_linsol("lapackqr") 2512 def test_solve(self): 2513 print(123) 2514 N = 3 2515 nrhs = 50 2516 np.random.seed(0) 2517 A = np.random.random((3,3)) 2518 B = np.random.random((3,50)) 2519 2520 2521 C = solve(A, B, "lapackqr", {"max_nrhs": 50}) 2522 C1 = solve(A, B, "lapackqr", {"max_nrhs": 20}) 2523 C2 = solve(A, B, "lapackqr") 2524 2525 self.checkarray(C, C1) 2526 self.checkarray(C1, C2) 2527 def test_sparse_lt(self): 2528 x = MX.sym("x",Sparsity.lower(5)) 2529 self.assertEqual((x>0).nnz(),5*6/2) 2530 self.assertEqual((x>=0).nnz(),5*5) 2531 def test_inv(self): 2532 np.random.seed(0) 2533 2534 for X in [SX, MX]: 2535 A = X.sym("x",3,3) 2536 Av = np.random.random((3,3)) 2537 f = Function('f',[A],[inv(A),inv(DM(Av)),A.__mpower__(-1), DM(Av).__mpower__(-1)]) 2538 out = f(Av) 2539 for o in out: 2540 self.checkarray(o, np.linalg.inv(np.array(Av))) 2541 2542 2543 def test_interp1d(self): 2544 v = [7,3,4,-3] 2545 vs = MX.sym("x",4,1) 2546 2547 xq = [10,0,1,2,4,8,7,5,1.5] 2548 x = [1,2,4,8] 2549 F = Function("f",[vs],[interp1d(x,vs,xq)]) 2550 2551 self.checkarray(F(v),np.interp(xq,x,v)) 2552 2553 F = Function("f",[vs],[interp1d(x,vs,xq,"floor")]) 2554 2555 self.checkarray(F(v),DM([-3,7,7,3,4,-3,4,4,7])) 2556 2557 F = Function("f",[vs],[interp1d(x,vs,xq,"ceil")]) 2558 2559 self.checkarray(F(v),DM([-3,7,7,3,4,-3,-3,-3,3])) 2560 2561 v = [7,3,4,-3] 2562 vs = MX.sym("x",4,1) 2563 2564 xq = [10,0,1,2,3,4,3.5,3.25,1.5] 2565 x = [1,2,3,4] 2566 F = Function("f",[vs],[interp1d(x,vs,xq,"linear",True)]) 2567 2568 self.checkarray(F(v),np.interp(xq,x,v)) 2569 2570 def test_bilin_etc(self): 2571 x = MX.sym("x",3,3) 2572 y = MX.sym("y",3,1) 2573 z = MX.sym("z",3,1) 2574 2575 import numpy 2576 numpy.random.seed(42) 2577 x0 = numpy.random.random((3,3)) 2578 y0 = numpy.random.random((3,1)) 2579 z0 = numpy.random.random((3,1)) 2580 2581 for e in [(bilin(x,y,z),mtimes(mtimes(y.T,x),z)),(rank1(x,0.3,y,z),x+0.3*mtimes(y,z.T))]: 2582 f = Function('f',[x,y,z],[e[0]]) 2583 fref = Function('fref',[x,y,z],[e[1]]) 2584 f_sx = f.expand() 2585 self.checkfunction(f,fref,inputs=[x0,y0,z0]) 2586 self.checkfunction(f_sx,fref,inputs=[x0,y0,z0]) 2587 self.check_codegen(f,inputs=[x0,y0,z0]) 2588 2589 def test_det_shape(self): 2590 X = MX.sym("x",2,3) 2591 with self.assertRaises(RuntimeError): 2592 det(X) 2593 X = MX.sym("x",3,3) 2594 det(X) 2595 2596 2597 @known_bug() 2598 def test_det(self): 2599 X = MX.sym("x",3,3) 2600 x = SX.sym("x",3,3) 2601 2602 import numpy 2603 numpy.random.seed(42) 2604 x0 = numpy.random.random((3,3)) 2605 2606 f = Function('f',[x],[det(x)]) 2607 F = Function('F',[X],[det(X)]) 2608 self.checkfunction(f,F,inputs=[x0]) 2609 2610 def test_mtimes_mismatch_segfault(self): 2611 with self.assertInException("incompatible dimensions"): 2612 mtimes(DM(Sparsity.lower(5)),MX.sym('x',100)) 2613 2614 def test_monitor(self): 2615 x = MX.sym("x") 2616 y = sqrt(x.monitor("hey")) 2617 2618 f = Function('f',[x],[y]) 2619 with capture_stdout() as out: 2620 f(3) 2621 self.assertTrue(out[0]=="hey:\n[3]\n") 2622 2623 with capture_stdout() as out2: 2624 self.check_codegen(f,inputs=[3]) 2625 if args.run_slow: 2626 self.assertTrue("hey:\n[3]\n" in out2[0]) 2627 2628 def test_codegen_specials(self): 2629 x = MX.sym("x") 2630 y = MX.sym("y") 2631 2632 for z in [ x**2, if_else(y>0,2*x,x*y), fmin(x,y), fmax(x,y), sign(x*y)]: 2633 f = Function('f',[x,y],[z]) 2634 self.check_codegen(f,inputs=[1,2]) 2635 self.check_codegen(f,inputs=[1,-2]) 2636 2637 2638 @memory_heavy() 2639 def test_getsetnonzeros(self): 2640 import numpy 2641 numpy.random.seed(42) 2642 for S in [Sparsity.lower(5),Sparsity.dense(5,5)]: 2643 M = MX.sym("X",S.nnz()) 2644 2645 m = sin(MX(S,M)) 2646 for i in [0,2]: 2647 for ind in [(i,i),(1,i),(i,1),(slice(None),i),(i,slice(None)),(slice(i,i+2),slice(i,i+2)),(slice(i,i+2),slice(None)),([i,i],[0,0]),([],[])]: 2648 E = m.__getitem__(ind) 2649 e = cos(E) 2650 e = dot(e,e) 2651 f = Function('f',[M],[e]) 2652 self.checkfunction(f,f.expand(),inputs=[ numpy.random.random((S.nnz(),1))]) 2653 2654 mc = m+0 2655 2656 Y = MX.sym("y",E.nnz()) 2657 y = MX(E.sparsity(),Y) 2658 mc.__setitem__(ind,y) 2659 e = cos(mc) 2660 e = dot(e,e) 2661 2662 f = Function('f',[M,Y],[e]) 2663 self.checkfunction(f,f.expand(),inputs=[ numpy.random.random((S.nnz(),1)), numpy.random.random((E.nnz(),1))]) 2664 2665 def test_evalf(self): 2666 x = MX.sym("x") 2667 2668 p = MX.sym("p") 2669 f = Function('f',[x],[sin(x)]) 2670 y = f.call([p],False,True)[0] 2671 y = substitute(y,p,3) 2672 2673 with self.assertInException("not defined"): 2674 y.to_DM() 2675 self.checkarray(evalf(y),sin(3)) 2676 with self.assertInException("since variables [x] are free"): 2677 evalf(x) 2678 2679 2680 def test_mmin(self): 2681 2682 def mx_eval(f,X,x): 2683 y = X.sym("y",x.sparsity()) 2684 return Function('f',[y],f.call([y],X is MX,False)) 2685 2686 for X in [SX,MX]: 2687 for mod in [lambda f:f, lambda f:f.expand(), lambda f: mx_eval(f,X,x)]: 2688 x = X.sym("x",2,2) 2689 2690 f = mod(Function('f',[x],[mmin(x),mmax(x)])) 2691 2692 2693 [mmin_res,mmax_res] = f(DM([[-2,3],[-3,5]])) 2694 2695 self.checkarray(mmin_res, -3) 2696 self.checkarray(mmax_res, 5) 2697 2698 x = X.sym("x",Sparsity.lower(2)) 2699 2700 f = mod(Function('f',[x],[mmin(x),mmax(x)])) 2701 2702 [mmin_res,mmax_res] = f(sparsify(DM([[-2,0],[-3,-5]]))) 2703 2704 self.checkarray(mmin_res, -5) 2705 self.checkarray(mmax_res, 0) 2706 2707 [mmin_res,mmax_res] = f(sparsify(DM([[2,0],[3,5]]))) 2708 2709 self.checkarray(mmin_res, 0) 2710 self.checkarray(mmax_res, 5) 2711 2712 x = X.sym("x",Sparsity(2,2)) 2713 f = mod(Function('f',[x],[mmin(x),mmax(x)])) 2714 [mmin_res,mmax_res] = f(DM(2,2)) 2715 2716 self.checkarray(mmin_res, 0) 2717 self.checkarray(mmax_res, 0) 2718 2719 x = X.sym("x",Sparsity(2,0)) 2720 f = mod(Function('f',[x],[mmin(x),mmax(x)])) 2721 [mmin_res,mmax_res] = f(DM(2,0)) 2722 2723 self.assertTrue(mmin_res.is_empty()) 2724 self.assertTrue(mmax_res.is_empty()) 2725 2726 x = X.sym("x",Sparsity(0,0)) 2727 f = mod(Function('f',[x],[mmin(x),mmax(x)])) 2728 [mmin_res,mmax_res] = f(DM(0,0)) 2729 2730 self.assertTrue(mmin_res.is_empty()) 2731 self.assertTrue(mmax_res.is_empty()) 2732 2733 def test_doc_expression_tools(self): 2734 self.assertTrue("Given a repeated matrix, computes the sum of repeated parts." in repsum.__doc__) 2735 2736 def test_densify(self): 2737 I = sparsify(DM([[0,1,0,1],[1,1,0,1],[0,1,1,0]])).sparsity() 2738 2739 x = MX.sym("x",I) 2740 2741 f = Function("f",[x],[densify(x)]) 2742 2743 a = DM([[0,1,0,2],[3,4,0,5],[0,6,7,0]]) 2744 2745 y = f(sparsify(a)) 2746 self.assertTrue(y.sparsity()==a.sparsity()) 2747 self.checkarray(y,a) 2748 2749 self.check_codegen(f,inputs=[a]) 2750 2751 2752 x = MX.sym("x",3,4) 2753 2754 f = Function("f",[x],[x[I]]) 2755 2756 a = DM([[0,1,0,2],[3,4,0,5],[0,6,7,0]]) 2757 b = sparsify(a) 2758 a = DM([[9,1,9,2],[3,4,9,5],[9,6,7,9]]) 2759 2760 y = f(a) 2761 self.assertTrue(y.sparsity()==b.sparsity()) 2762 self.checkarray(y,b) 2763 2764 self.check_codegen(f,inputs=[a]) 2765 self.check_serialize(f,inputs=[a]) 2766 2767 def test_constant_from_file(self): 2768 2769 open("test.txt", "w").write("5.6 1e-5") 2770 with self.assertInException("Failed to read a double from 'test.txt'. Expected 3 doubles."): 2771 x = MX(Sparsity.lower(2), "test.txt") 2772 2773 open("test.txt", "w").write("5.6 abc 1e-5") 2774 with self.assertInException("Failed to read a double from 'test.txt'. Expected 3 doubles."): 2775 x = MX(Sparsity.lower(2), "test.txt") 2776 2777 with self.assertInException("Cannot open file 'testQWHAL567p.txt'."): 2778 x = MX(Sparsity.lower(2), "testQWHAL567p.txt") 2779 2780 open("test.txt", "w").write("5.6 1e-5 -12") 2781 x = MX(Sparsity.lower(2), "test.txt") 2782 2783 a = MX.sym("a") 2784 2785 f = Function("f",[a],[a*x]) 2786 self.checkarray(f(2),DM(Sparsity.lower(2),[2*5.6,2*1e-5,-2*12])) 2787 2788 self.check_codegen(f,inputs=[2]) 2789 2790 with self.assertInException("Not defined for ConstantFile"): 2791 x.to_DM() 2792 2793 def test_nonzeros_param(self): 2794 x = MX.sym("x",2) 2795 y = MX.sym("y") 2796 2797 for ad_weight_sp in [0,1]: 2798 opts = {"helper_options":{"ad_weight_sp":ad_weight_sp}} 2799 z = x.nz[y] 2800 J = jacobian(z,x,opts) 2801 self.assertTrue(J.size()==(1,2)) 2802 self.assertTrue(J.nnz()==2) 2803 J = jacobian(z,y,opts) 2804 self.assertTrue(J.size()==(1,1)) 2805 self.assertTrue(J.nnz()==0) 2806 2807 q = MX.sym("q") 2808 z = MX(x) 2809 z.nz[y] = q 2810 2811 J = jacobian(z,x,opts) 2812 self.assertTrue(J.size()==(2,2)) 2813 self.assertTrue(J.nnz()==4) 2814 J = jacobian(z,y,opts) 2815 self.assertTrue(J.size()==(2,1)) 2816 self.assertTrue(J.nnz()==0) 2817 J = jacobian(z,q,opts) 2818 self.assertTrue(J.size()==(2,1)) 2819 self.assertTrue(J.nnz()==2) 2820 2821 def test_low(self): 2822 v = MX.sym("v",5) 2823 p0 = [-1,0,0.2,1,1.3,2,2.5,3,3.5,4,4.5] 2824 p = MX.sym("p",len(p0)) 2825 2826 2827 f = Function("f",[v,p],[low(v, p)]) 2828 2829 inputs = [list(range(5)),p0] 2830 self.check_codegen(f,inputs=inputs) 2831 self.check_serialize(f,inputs=inputs) 2832 2833 self.checkarray(f(*inputs),DM([0, 0, 0, 1, 1, 2, 2, 3, 3, 3, 3])) 2834 2835 def test_linear_interpn(self): 2836 N = 4 2837 n_dim=1 2838 nq = 3 2839 n_out = 2 2840 grid = [] 2841 x0 = [] 2842 for i in range(n_dim): 2843 g = sorted(list(np.random.random(N))) 2844 grid.append(g) 2845 x0.append(np.linspace(g[0],g[-1],nq+2)[1:-1]) 2846 2847 x0r = list(zip(*x0)) 2848 2849 D = np.random.random([n_out]+[N]*n_dim) 2850 2851 d = D.ravel(order='F') 2852 2853 strides = [n_out] 2854 for i in range(n_dim): 2855 strides.append(strides[-1]*N) 2856 2857 F = interpolant('F', 'linear', grid, d) 2858 2859 ref = vcat([F(x).T for x in x0r]) 2860 2861 x = [MX.sym("x%d" % i,1,N) for i in range(n_dim)] 2862 v = MX.sym("v",N**n_dim*n_out) 2863 xq = [MX.sym("xq%d" % i,nq) for i in range(n_dim)] 2864 2865 e = MX.interpn_linear(x,v,xq) 2866 2867 f = Function('f',[hcat(x),v,vcat(xq)],[e]) 2868 inputs = [hcat([hcat(e) for e in grid]),d,vcat(x0)] 2869 r = f(*inputs) 2870 self.checkarray(r, ref) 2871 self.check_codegen(f,inputs=inputs,std="c99") 2872 self.check_serialize(f,inputs=inputs) 2873 2874 def test_paramslice(self): 2875 N = 4 2876 M = 6 2877 A = DM.rand(N,M) 2878 As = MX.sym("A", A.shape) 2879 ASs = MX.sym("AS",tril(As.sparsity())) 2880 print(ASs.sparsity().spy()) 2881 vr=list(range(N)) 2882 vc=list(range(M)) 2883 vk=list(range(M*N)) 2884 for r in [1, slice(1,3), slice(0,N,2), slice(None)]: 2885 rs = MX.sym("r",A[r,:].shape[0]) 2886 for c in [1, slice(1,3), slice(0,M,2), slice(None)]: 2887 cs = MX.sym("c",A[:,c].shape[1]) 2888 for rr in [r, rs]: 2889 for cc in [c, cs]: 2890 f = Function('f',[As,rs,cs],[As[rr,cc]]) 2891 self.checkarray(f(A,vr[r],vc[c]),A[r,c]) 2892 if isinstance(rr,MX) or isinstance(cc,MX): 2893 with self.assertInException("Parametric slicing only supported for dense matrices"): 2894 ASs[rr,cc] 2895 for k in [1, slice(1,3), slice(0,N,2), slice(0,2*N,2), slice(None)]: 2896 ks = MX.sym("r",A[k].shape[0]) 2897 for kk in [k, ks]: 2898 f = Function('f',[As,ks],[As[kk]]) 2899 self.checkarray(f(A,vk[k]),A[k]) 2900 with self.assertInException("Parametric slicing only supported for dense matrices"): 2901 ASs[ks] 2902 2903 def test_mapping(self): 2904 A = MX.sym("A",4,4) 2905 i = DM([[0,3],[1,2]]) 2906 self.checkarray(i,A[i].mapping()) 2907 2908 2909 def test_convexify(self): 2910 A = diagcat(1,2,-1,blockcat([[1.2,1.3],[1.3,4]]),sparsify(blockcat([[0,1,0],[1,4,7],[0,7,9]])),DM(2,2)) 2911 2912 margin = 1e-7 2913 2914 np.random.seed(0) 2915 p = np.random.permutation(A.shape[0]) 2916 2917 2918 [D,V] = np.linalg.eig(np.array(A)) 2919 Dr= fmax(abs(D),1e-7) 2920 Dc= fmax(D,1e-7) 2921 2922 Ar_ref = mtimes(mtimes(V,diag(Dr)),V.T) 2923 Ac_ref = mtimes(mtimes(V,diag(Dc)),V.T) 2924 As = MX.sym("As",A.sparsity()) 2925 2926 for opts,ref in [({"strategy":"eigen-reflect"},Ar_ref), ({"strategy":"eigen-clip"},Ac_ref), ({"strategy":"regularize"},A+(4+margin)*DM.eye(A.shape[0]))]: 2927 2928 2929 ops = [lambda e: e, lambda e: triu(e), lambda e: tril(e), 2930 lambda e: e[p,p], lambda e: triu(e[p,p]), lambda e: tril(e[p,p])] 2931 if "regularize" in str(opts): ops = [lambda e: e, lambda e: e[p,p]] 2932 for op in ops: 2933 2934 f= Function("f",[As],[convexify(op(A),opts)]) 2935 2936 self.checkarray(f(A),op(ref),digits=8) 2937 2938 self.check_serialize(f,inputs=[A]) 2939 2940 self.check_codegen(f,inputs=[A]) 2941 2942 2943if __name__ == '__main__': 2944 unittest.main() 2945