code: def _smw_solver(s, A, AtA, Qi, di):
      r"""
      Returns a solver for the linear system:
  
      .. math::
  
          (sI + ABA^\prime) y = x
  
      The returned function f satisfies f(x) = y as defined above.
  
      B and its inverse matrix are block diagonal.  The upper left block
      of :math:`B^{-1}` is Qi and its lower right block is diag(di).
  
      Parameters
      ----------
      s : scalar
          See above for usage
      A : ndarray
          p x q matrix, in general q << p, may be sparse.
      AtA : square ndarray
          :math:`A^\prime  A`, a q x q matrix.
      Qi : square symmetric ndarray
          The matrix `B` is q x q, where q = r + d.  `B` consists of a r
          x r diagonal block whose inverse is `Qi`, and a d x d diagonal
          block, whose inverse is diag(di).
      di : 1d array_like
          See documentation for Qi.
  
      Returns
      -------
      A function for solving a linear system, as documented above.
  
      Notes
      -----
      Uses Sherman-Morrison-Woodbury identity:
          https://en.wikipedia.org/wiki/Woodbury_matrix_identity
      """
  
      # Use SMW identity
      qmat = AtA / s
      m = Qi.shape[0]
      qmat[0:m, 0:m] += Qi
  
      if sparse.issparse(A):
          qmat[m:, m:] += sparse.diags(di)
  
          def solver(rhs):
              ql = A.T.dot(rhs)
              # Based on profiling, the next line can be the
              # majority of the entire run time of fitting the model.
              ql = sparse.linalg.spsolve(qmat, ql)
              if ql.ndim < rhs.ndim:
                  # spsolve squeezes nx1 rhs
                  ql = ql[:, None]
              ql = A.dot(ql)
              return rhs / s - ql / s**2
  
      else:
          d = qmat.shape[0]
          qmat.flat[m * (d + 1) :: d + 1] += di
          qmati = np.linalg.solve(qmat, A.T)
  
          def solver(rhs):
              # A is tall and qmati is wide, so we want
              # A * (qmati * rhs) not (A * qmati) * rhs
              ql = np.dot(qmati, rhs)
              ql = np.dot(A, ql)
              return rhs / s - ql / s**2
  
      return solver
  
imports: import warnings
  import numpy
  import pandas
  import patsy
  from scipy import sparse
  from scipy.stats.distributions import norm
  from statsmodels.base._penalties import Penalty
  import statsmodels.base.model
  from statsmodels.tools import data
  from statsmodels.tools.decorators import cache_readonly
  from statsmodels.tools.sm_exceptions import ConvergenceWarning
label_desc: Use image node_id 5 for calling a global function with example usage: _smw_solver(s, A, AtA, Qi, di) and returns: solver, unknown, unknown