2022-08-20T10:05:12.4606678Z =================================== FAILURES ===================================
2022-08-20T10:05:12.4607612Z _________________ test_ispassive_edge_cases[test_input2-True] __________________
2022-08-20T10:05:12.4607862Z 
2022-08-20T10:05:12.4608255Z test_input = (array([[-3.e+12,  0.e+00],
2022-08-20T10:05:12.4608627Z        [ 0.e+00, -2.e+12]]), array([[0],
2022-08-20T10:05:12.4609085Z        [1]]), array([[-1,  2]]), array([[0.]]))
2022-08-20T10:05:12.4609312Z expected = True
2022-08-20T10:05:12.4609448Z 
2022-08-20T10:05:12.4609567Z     @pytest.mark.parametrize(
2022-08-20T10:05:12.4609948Z         "test_input,expected",
2022-08-20T10:05:12.4610188Z         [((A, B, C, D*0.0), True),
2022-08-20T10:05:12.4610415Z          ((A_d, B, C, D), True),
2022-08-20T10:05:12.4610640Z          ((A*1e12, B, C, D*0), True),
2022-08-20T10:05:12.4610854Z          ((A, B*0, C*0, D), True),
2022-08-20T10:05:12.4611073Z          ((A*0, B, C, D), True),
2022-08-20T10:05:12.4611293Z          ((A*0, B*0, C*0, D*0), True)])
2022-08-20T10:05:12.4611683Z     def test_ispassive_edge_cases(test_input, expected):
2022-08-20T10:05:12.4611958Z         A = test_input[0]
2022-08-20T10:05:12.4612178Z         B = test_input[1]
2022-08-20T10:05:12.4612610Z         C = test_input[2]
2022-08-20T10:05:12.4612970Z         D = test_input[3]
2022-08-20T10:05:12.4613469Z         sys = ss(A, B, C, D)
2022-08-20T10:05:12.4613753Z >       assert(passivity.ispassive(sys) == expected)
2022-08-20T10:05:12.4613952Z 
2022-08-20T10:05:12.4614077Z control/tests/passivity_test.py:115: 
2022-08-20T10:05:12.4614373Z _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
2022-08-20T10:05:12.4614742Z control/passivity.py:292: in ispassive
2022-08-20T10:05:12.4615096Z     return solve_passivity_LMI(sys, rho=ofp_index, nu=ifp_index) is not None
2022-08-20T10:05:12.4615456Z control/passivity.py:177: in solve_passivity_LMI
2022-08-20T10:05:12.4615750Z     sol = cvx.solvers.sdp(c, Gs=Gs, hs=hs)
2022-08-20T10:05:12.4616316Z /usr/share/miniconda3/envs/test-env/lib/python3.10/site-packages/cvxopt/coneprog.py:4126: in sdp
2022-08-20T10:05:12.4616810Z     sol = conelp(c, G, h, dims, A = A, b = b, primalstart = ps, dualstart = ds, kktsolver = kktsolver, options = options)
2022-08-20T10:05:12.4617643Z /usr/share/miniconda3/envs/test-env/lib/python3.10/site-packages/cvxopt/coneprog.py:1395: in conelp
2022-08-20T10:05:12.4618038Z     misc.update_scaling(W, lmbda, ds, dz)
2022-08-20T10:05:12.4618322Z _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
2022-08-20T10:05:12.4618485Z 
2022-08-20T10:05:12.4618810Z W = {'beta': [], 'd': <0x1 matrix, tc='d'>, 'di': <0x1 matrix, tc='d'>, 'dnl': <0x1 matrix, tc='d'>, ...}
2022-08-20T10:05:12.4619335Z lmbda = <6x1 matrix, tc='d'>, s = <13x1 matrix, tc='d'>
2022-08-20T10:05:12.4619636Z z = <13x1 matrix, tc='d'>
2022-08-20T10:05:12.4619773Z 
2022-08-20T10:05:12.4619885Z     def update_scaling(W, lmbda, s, z):
2022-08-20T10:05:12.4620221Z         """
2022-08-20T10:05:12.4620607Z         Updates the Nesterov-Todd scaling matrix W and the scaled variable
2022-08-20T10:05:12.4621042Z         lmbda so that on exit
2022-08-20T10:05:12.4621247Z     
2022-08-20T10:05:12.4621547Z               W * zt = W^{-T} * st = lmbda.
2022-08-20T10:05:12.4621773Z     
2022-08-20T10:05:12.4622145Z         On entry, the nonlinear, 'l' and 'q' components of the arguments s
2022-08-20T10:05:12.4622629Z         and z contain W^{-T}*st and W*zt, i.e, the new iterates in the current
2022-08-20T10:05:12.4622920Z         scaling.
2022-08-20T10:05:12.4623119Z     
2022-08-20T10:05:12.4623489Z         The 's' components contain the factors Ls, Lz in a factorization of
2022-08-20T10:05:12.4624030Z         the new iterates in the current scaling, W^{-T}*st = Ls*Ls',
2022-08-20T10:05:12.4624341Z         W*zt = Lz*Lz'.
2022-08-20T10:05:12.4624645Z         """
2022-08-20T10:05:12.4624832Z     
2022-08-20T10:05:12.4625008Z     
2022-08-20T10:05:12.4625390Z         # Nonlinear and 'l' blocks
2022-08-20T10:05:12.4625606Z         #
2022-08-20T10:05:12.4626008Z         #    d :=  d .* sqrt( s ./ z )
2022-08-20T10:05:12.4626266Z         #    lmbda := lmbda .* sqrt(s) .* sqrt(z)
2022-08-20T10:05:12.4626501Z     
2022-08-20T10:05:12.4626752Z         if 'dnl' in W:
2022-08-20T10:05:12.4627026Z             mnl = len(W['dnl'])
2022-08-20T10:05:12.4627256Z         else:
2022-08-20T10:05:12.4627460Z             mnl = 0
2022-08-20T10:05:12.4627713Z         ml = len(W['d'])
2022-08-20T10:05:12.4627934Z         m = mnl + ml
2022-08-20T10:05:12.4628165Z         s[:m] = base.sqrt( s[:m] )
2022-08-20T10:05:12.4628397Z         z[:m] = base.sqrt( z[:m] )
2022-08-20T10:05:12.4629372Z     
2022-08-20T10:05:12.4629630Z         # d := d .* s .* z
2022-08-20T10:05:12.4629940Z         if 'dnl' in W:
2022-08-20T10:05:12.4630299Z             blas.tbmv(s, W['dnl'], n = mnl, k = 0, ldA = 1)
2022-08-20T10:05:12.4630693Z             blas.tbsv(z, W['dnl'], n = mnl, k = 0, ldA = 1)
2022-08-20T10:05:12.4631031Z             W['dnli'][:] = W['dnl'][:] ** -1
2022-08-20T10:05:12.4631420Z         blas.tbmv(s, W['d'], n = ml, k = 0, ldA = 1, offsetA = mnl)
2022-08-20T10:05:12.4631840Z         blas.tbsv(z, W['d'], n = ml, k = 0, ldA = 1, offsetA = mnl)
2022-08-20T10:05:12.4632171Z         W['di'][:] = W['d'][:] ** -1
2022-08-20T10:05:12.4632391Z     
2022-08-20T10:05:12.4632808Z         # lmbda := s .* z
2022-08-20T10:05:12.4633059Z         blas.copy(s, lmbda, n = m)
2022-08-20T10:05:12.4633323Z         blas.tbmv(z, lmbda, n = m, k = 0, ldA = 1)
2022-08-20T10:05:12.4633622Z     
2022-08-20T10:05:12.4633806Z     
2022-08-20T10:05:12.4634057Z         # 'q' blocks.
2022-08-20T10:05:12.4634267Z         #
2022-08-20T10:05:12.4634539Z         # Let st and zt be the new variables in the old scaling:
2022-08-20T10:05:12.4634792Z         #
2022-08-20T10:05:12.4635010Z         #     st = s_k,   zt = z_k
2022-08-20T10:05:12.4635228Z         #
2022-08-20T10:05:12.4635560Z         # and a = sqrt(st' * J * st),  b = sqrt(zt' * J * zt).
2022-08-20T10:05:12.4635811Z         #
2022-08-20T10:05:12.4636209Z         # 1. Compute the hyperbolic Householder transformation 2*q*q' - J
2022-08-20T10:05:12.4636519Z         #    that maps st/a to zt/b.
2022-08-20T10:05:12.4636747Z         #
2022-08-20T10:05:12.4637058Z         #        c = sqrt( (1 + st'*zt/(a*b)) / 2 )
2022-08-20T10:05:12.4637420Z         #        q = (st/a + J*zt/b) / (2*c).
2022-08-20T10:05:12.4637667Z         #
2022-08-20T10:05:12.4637894Z         #    The new scaling point is
2022-08-20T10:05:12.4638110Z         #
2022-08-20T10:05:12.4638578Z         #        wk := betak * sqrt(a/b) * (2*v[k]*v[k]' - J) * q
2022-08-20T10:05:12.4638818Z         #
2022-08-20T10:05:12.4639082Z         #    with betak = W['beta'][k].
2022-08-20T10:05:12.4639302Z         #
2022-08-20T10:05:12.4639636Z         # 3. The scaled variable:
2022-08-20T10:05:12.4639831Z         #
2022-08-20T10:05:12.4640044Z         #        lambda_k0 = sqrt(a*b) * c
2022-08-20T10:05:12.4640423Z         #        lambda_k1 = sqrt(a*b) * ( (2vk*vk' - J) * (-d*q + u/2) )_1
2022-08-20T10:05:12.4640655Z         #
2022-08-20T10:05:12.4640839Z         #    where
2022-08-20T10:05:12.4641026Z         #
2022-08-20T10:05:12.4641267Z         #        u = st/a - J*zt/b
2022-08-20T10:05:12.4641620Z         #        d = ( vk0 * (vk'*u) + u0/2 ) / (2*vk0 *(vk'*q) - q0 + 1).
2022-08-20T10:05:12.4641864Z         #
2022-08-20T10:05:12.4642182Z         # 4. Update scaling
2022-08-20T10:05:12.4642379Z         #
2022-08-20T10:05:12.4642577Z         #        v[k] := wk^1/2
2022-08-20T10:05:12.4642824Z         #              = 1 / sqrt(2*(wk0 + 1)) * (wk + e).
2022-08-20T10:05:12.4643060Z         #        beta[k] *=  sqrt(a/b)
2022-08-20T10:05:12.4643271Z     
2022-08-20T10:05:12.4643455Z         ind = m
2022-08-20T10:05:12.4643731Z         for k in range(len(W['v'])):
2022-08-20T10:05:12.4643952Z     
2022-08-20T10:05:12.4644314Z             v = W['v'][k]
2022-08-20T10:05:12.4644524Z             m = len(v)
2022-08-20T10:05:12.4644728Z     
2022-08-20T10:05:12.4645040Z             # ln = sqrt( lambda_k' * J * lambda_k )
2022-08-20T10:05:12.4645428Z             ln = jnrm2(lmbda, n = m, offset = ind)
2022-08-20T10:05:12.4645658Z     
2022-08-20T10:05:12.4645978Z             # a = sqrt( sk' * J * sk ) = sqrt( st' * J * st )
2022-08-20T10:05:12.4646340Z             # s := s / a = st / a
2022-08-20T10:05:12.4646607Z             aa = jnrm2(s, offset = ind, n = m)
2022-08-20T10:05:12.4646898Z             blas.scal(1.0/aa, s, offset = ind, n = m)
2022-08-20T10:05:12.4647123Z     
2022-08-20T10:05:12.4647455Z             # b = sqrt( zk' * J * zk ) = sqrt( zt' * J * zt )
2022-08-20T10:05:12.4647818Z             # z := z / a = zt / b
2022-08-20T10:05:12.4648054Z             bb = jnrm2(z, offset = ind, n = m)
2022-08-20T10:05:12.4648334Z             blas.scal(1.0/bb, z, offset = ind, n = m)
2022-08-20T10:05:12.4648565Z     
2022-08-20T10:05:12.4648857Z             # c = sqrt( ( 1 + (st'*zt) / (a*b) ) / 2 )
2022-08-20T10:05:12.4649290Z             cc = math.sqrt( ( 1.0 + blas.dot(s, z, offsetx = ind, offsety =
2022-08-20T10:05:12.4649582Z                 ind, n = m) ) / 2.0 )
2022-08-20T10:05:12.4649796Z     
2022-08-20T10:05:12.4650041Z             # vs = v' * st / a
2022-08-20T10:05:12.4650311Z             vs = blas.dot(v, s, offsety = ind, n = m)
2022-08-20T10:05:12.4650549Z     
2022-08-20T10:05:12.4651253Z             # vz = v' * J *zt / b
2022-08-20T10:05:12.4651724Z             vz = jdot(v, z, offsety = ind, n = m)
2022-08-20T10:05:12.4651960Z     
2022-08-20T10:05:12.4652334Z             # vq = v' * q where q = (st/a + J * zt/b) / (2 * c)
2022-08-20T10:05:12.4652810Z             vq = (vs + vz ) / 2.0 / cc
2022-08-20T10:05:12.4653032Z     
2022-08-20T10:05:12.4653355Z             # vu = v' * u  where u =  st/a - J * zt/b
2022-08-20T10:05:12.4653655Z             vu = vs - vz
2022-08-20T10:05:12.4653865Z     
2022-08-20T10:05:12.4654053Z             # lambda_k0 = c
2022-08-20T10:05:12.4654285Z             lmbda[ind] = cc
2022-08-20T10:05:12.4654602Z     
2022-08-20T10:05:12.4654873Z             # wk0 = 2 * vk0 * (vk' * q) - q0
2022-08-20T10:05:12.4655466Z             wk0 = 2 * v[0] * vq - ( s[ind] + z[ind] ) / 2.0 / cc
2022-08-20T10:05:12.4655699Z     
2022-08-20T10:05:12.4655989Z             # d = (v[0] * (vk' * u) - u0/2) / (wk0 + 1)
2022-08-20T10:05:12.4656492Z             dd = (v[0] * vu - s[ind]/2.0 + z[ind]/2.0) / (wk0 + 1.0)
2022-08-20T10:05:12.4656755Z     
2022-08-20T10:05:12.4657096Z             # lambda_k1 = 2 * v_k1 * vk' * (-d*q + u/2) - d*q1 + u1/2
2022-08-20T10:05:12.4657651Z             blas.copy(v, lmbda, offsetx = 1, offsety = ind+1, n = m-1)
2022-08-20T10:05:12.4658095Z             blas.scal(2.0 * (-dd * vq + 0.5 * vu), lmbda, offset = ind+1,
2022-08-20T10:05:12.4658422Z                n = m-1)
2022-08-20T10:05:12.4658904Z             blas.axpy(s, lmbda, 0.5 * (1.0 - dd/cc), offsetx = ind+1, offsety
2022-08-20T10:05:12.4659242Z                = ind+1, n = m-1)
2022-08-20T10:05:12.4659529Z             blas.axpy(z, lmbda, 0.5 * (1.0 + dd/cc), offsetx = ind+1, offsety
2022-08-20T10:05:12.4659850Z                = ind+1, n = m-1)
2022-08-20T10:05:12.4660056Z     
2022-08-20T10:05:12.4660408Z             # Scale so that sqrt(lambda_k' * J * lambda_k) = sqrt(aa*bb).
2022-08-20T10:05:12.4660855Z             blas.scal(math.sqrt(aa*bb), lmbda, offset = ind, n = m)
2022-08-20T10:05:12.4661118Z     
2022-08-20T10:05:12.4661396Z             # v := (2*v*v' - J) * q
2022-08-20T10:05:12.4661843Z             #    = 2 * (v'*q) * v' - (J* st/a + zt/b) / (2*c)
2022-08-20T10:05:12.4662220Z             blas.scal(2.0 * vq, v)
2022-08-20T10:05:12.4662532Z             v[0] -= s[ind] / 2.0 / cc
2022-08-20T10:05:12.4662939Z             blas.axpy(s, v,  0.5/cc, offsetx = ind+1, offsety = 1, n = m-1)
2022-08-20T10:05:12.4663341Z             blas.axpy(z, v, -0.5/cc, offsetx = ind, n = m)
2022-08-20T10:05:12.4663589Z     
2022-08-20T10:05:12.4663824Z             # v := v^{1/2} = 1/sqrt(2 * (v0 + 1)) * (v + e)
2022-08-20T10:05:12.4664050Z             v[0] += 1.0
2022-08-20T10:05:12.4664310Z             blas.scal(1.0 / math.sqrt(2.0 * v[0]), v)
2022-08-20T10:05:12.4664548Z     
2022-08-20T10:05:12.4664751Z             # beta[k] *= ( aa / bb )**1/2
2022-08-20T10:05:12.4665097Z             W['beta'][k] *= math.sqrt( aa / bb )
2022-08-20T10:05:12.4665331Z     
2022-08-20T10:05:12.4665510Z             ind += m
2022-08-20T10:05:12.4665708Z     
2022-08-20T10:05:12.4665895Z     
2022-08-20T10:05:12.4666158Z         # 's' blocks
2022-08-20T10:05:12.4666360Z         #
2022-08-20T10:05:12.4666633Z         # Let st, zt be the updated variables in the old scaling:
2022-08-20T10:05:12.4666886Z         #
2022-08-20T10:05:12.4667184Z         #     st = Ls * Ls', zt = Lz * Lz'.
2022-08-20T10:05:12.4667411Z         #
2022-08-20T10:05:12.4667735Z         # where Ls and Lz are the 's' components of s, z.
2022-08-20T10:05:12.4667992Z         #
2022-08-20T10:05:12.4668302Z         # 1.  SVD Lz'*Ls = Uk * lambda_k^+ * Vk'.
2022-08-20T10:05:12.4668526Z         #
2022-08-20T10:05:12.4668739Z         # 2.  New scaling is
2022-08-20T10:05:12.4668957Z         #
2022-08-20T10:05:12.4669289Z         #         r[k] := r[k] * Ls * Vk * diag(lambda_k^+)^{-1/2}
2022-08-20T10:05:12.4669698Z         #         rti[k] := r[k] * Lz * Uk * diag(lambda_k^+)^{-1/2}.
2022-08-20T10:05:12.4669951Z         #
2022-08-20T10:05:12.4670124Z     
2022-08-20T10:05:12.4670493Z         work = matrix(0.0, (max( [0] + [r.size[0] for r in W['r']])**2, 1))
2022-08-20T10:05:12.4671033Z         ind = mnl + ml + sum([ len(v) for v in W['v'] ])
2022-08-20T10:05:12.4671299Z         ind2, ind3 = ind, 0
2022-08-20T10:05:12.4671594Z         for k in range(len(W['r'])):
2022-08-20T10:05:12.4671924Z             r, rti = W['r'][k], W['rti'][k]
2022-08-20T10:05:12.4672164Z             m = r.size[0]
2022-08-20T10:05:12.4672356Z     
2022-08-20T10:05:12.4672564Z             # r := r*sk = r*Ls
2022-08-20T10:05:12.4672858Z             blas.gemm(r, s, work, m = m, n = m, k = m, ldB = m, ldC = m,
2022-08-20T10:05:12.4673126Z                 offsetB = ind2)
2022-08-20T10:05:12.4673545Z             blas.copy(work, r, n = m**2)
2022-08-20T10:05:12.4673885Z     
2022-08-20T10:05:12.4674070Z             # rti := rti*zk = rti*Lz
2022-08-20T10:05:12.4674352Z             blas.gemm(rti, z, work, m = m, n = m, k = m, ldB = m, ldC = m,
2022-08-20T10:05:12.4674620Z                 offsetB = ind2)
2022-08-20T10:05:12.4674939Z             blas.copy(work, rti, n = m**2)
2022-08-20T10:05:12.4675176Z     
2022-08-20T10:05:12.4696294Z             # SVD Lz'*Ls = U * lmbds^+ * V'; store U in sk and V' in zk.
2022-08-20T10:05:12.4696880Z             blas.gemm(z, s, work, transA = 'T', m = m, n = m, k = m, ldA = m,
2022-08-20T10:05:12.4697207Z                 ldB = m, ldC = m, offsetA = ind2, offsetB = ind2)
2022-08-20T10:05:12.4697652Z             lapack.gesvd(work, lmbda, jobu = 'A', jobvt = 'A', m = m, n = m,
2022-08-20T10:05:12.4698006Z                 ldA = m, U = s, Vt = z, ldU = m, ldVt = m, offsetS = ind,
2022-08-20T10:05:12.4698319Z                 offsetU = ind2, offsetVt = ind2)
2022-08-20T10:05:12.4698543Z     
2022-08-20T10:05:12.4698741Z             # r := r*V
2022-08-20T10:05:12.4699127Z             blas.gemm(r, z, work, transB = 'T', m = m, n = m, k = m, ldB = m,
2022-08-20T10:05:12.4699410Z                 ldC = m, offsetB = ind2)
2022-08-20T10:05:12.4699675Z             blas.copy(work, r, n = m**2)
2022-08-20T10:05:12.4700021Z     
2022-08-20T10:05:12.4700204Z             # rti := rti*U
2022-08-20T10:05:12.4700602Z             blas.gemm(rti, s, work, n = m, m = m, k = m, ldB = m, ldC = m,
2022-08-20T10:05:12.4700867Z                 offsetB = ind2)
2022-08-20T10:05:12.4701098Z             blas.copy(work, rti, n = m**2)
2022-08-20T10:05:12.4701318Z     
2022-08-20T10:05:12.4701645Z             # r := r*lambda^{-1/2}; rti := rti*lambda^{-1/2}
2022-08-20T10:05:12.4702010Z             for i in range(m):
2022-08-20T10:05:12.4702263Z >               a = 1.0 / math.sqrt(lmbda[ind+i])
2022-08-20T10:05:12.4702566Z E               ZeroDivisionError: float division by zero
2022-08-20T10:05:12.4702754Z 
2022-08-20T10:05:12.4703110Z /usr/share/miniconda3/envs/test-env/lib/python3.10/site-packages/cvxopt/misc.py:628: ZeroDivisionError