Build system and test suite update

[bnavigator]

• setup.py install is deprecated. Use this opportunity to migrate to PEP517 style packaging: fixes Building wheel for dev version fails #758
• Update test matrix to test Python 3.10
• Don't cancel other CI jobs if one fails
• Build Slycot from source and test with PyPI packages instead of Conda
• Remove custom xvfb setup. Run pytest-xvfb with QtAgg backend for the zoom toolbar in the extra full conda test run

[murrayrm]

I saw something similar occur when I rebased #754. The CI tests worked on my local branch but failed in python-control (even though the code was identical) with the DISPLAY error. I re-ran the tests in python-control and they completed, so figured it was a one-time glitch.

I also noticed there was a warning about mambaforge, but I'm not able to find it. Something about how something was an experimental implementation and resolution of some dependencies may not be repeatable.

[bnavigator]

I also saw it in my local branch several times. It fails, but which of the many jobs fail is random (see #756)

[bnavigator]

Switching back from mamba to conda gets rid of the "Mamba support is still experimental and can result in differently solved environments!" warning you were mentioning @murrayrm, but even that didn't help:

https://github.com/bnavigator/python-control/commits/update-CI

[coveralls]

Coverage remained the same at 94.738% when pulling 69fb34a on bnavigator:setuptools-update into 87a9c57 on python-control:main.

[bnavigator]

I give up.

I think the random "invalid DISPLAY" and segmentation faults are not caused by us but rather either by the conda environment or the Github Actions runners. This might even resolve itself through the next release cycle of either of these.

I would say we just ignore the occassional random CI job fail and look at the remaining successful ones for determining whether a PR is good. We can always re-run a job and see if the failure is reproducible.

[bnavigator]

The X Display errors seem to be gone. New failures with the balred routines using conda Slycot producing possibly valid but differently ordered systems. Needs to be checked.

[murrayrm]

Odd that Slycot from source is OK, but the other Slycots aren't...

[bnavigator]

It's conda slycot using conda openblas vs. ubuntu netlib lapack/blas

[bnavigator]

Whack-a-mole!

Now this cvxopt test is flaky: passivity_test.py::test_ispassive_edge_cases :

https://github.com/python-control/python-control/runs/7930805491?check_suite_focus=true#step:6:3450

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

[bnavigator]

@murrayrm I think this is ready. I guess, we will have to hit "re-run failed jobs" on the GHA page until #761 is resolved.

[murrayrm]

This looks good to go, but a few small questions as I looked through the code.