docs: def onecons_qcqp(z, f, tol=1e-6):
      """ Solves a nonconvex problem
        minimize ||x-z||_2^2
        subject to f(x) = x^T P x + q^T x + r ~ 0
        where the relation ~ is given by f.relop (either <= or ==)
      """
  
      # if constraint is ineq and z is feasible: z is the solution
      if f.relop == '<=' and f.eval(z) <= 0:
          return z
  
      if f.eigh is None:
          Psymm = (f.P + f.P.T)/2.
          f.eigh = LA.eigh(np.asarray(Psymm.todense()))
  
      lmb, Q = f.eigh
      zhat = Q.T.dot(z)
      qhat = Q.T.dot(f.qarray)
  
      # now solve a transformed problem
      # minimize ||xhat - zhat||_2^2
      # subject to sum(lmb_i xhat_i^2) + qhat^T xhat + r = 0
      # constraint is now equality from
      # complementary slackness
      xhat = lambda nu: -np.divide(nu*qhat-2*zhat, 2*(1+nu*lmb))
      phi = lambda xhat: lmb.dot(np.power(xhat, 2)) + qhat.dot(xhat) + f.r
  
      s = -np.inf
      e = np.inf
      for l in lmb:
          if l > 0: s = max(s, -1./l)
          if l < 0: e = min(e, -1./l)
      if s == -np.inf:
          s = -1.
          while phi(xhat(s)) <= 0: s *= 2.
      if e == np.inf:
          e = 1.
          while phi(xhat(e)) >= 0: e *= 2.
      while e-s > tol:
          m = (s+e)/2.
          p = phi(xhat(m))
          if p > 0: s = m
          elif p < 0: e = m
          else:
              s = e = m
              break
      nu = (s+e)/2.
      return Q.dot(xhat(nu))