code: def hessian(f, varlist, constraints=()):
      """Compute Hessian matrix for a function f wrt parameters in varlist
      which may be given as a sequence or a row/column vector. A list of
      constraints may optionally be given.
  
      Examples
      ========
  
      >>> from sympy import Function, hessian, pprint
      >>> from sympy.abc import x, y
      >>> f = Function('f')(x, y)
      >>> g1 = Function('g')(x, y)
      >>> g2 = x**2 + 3*y
      >>> pprint(hessian(f, (x, y), [g1, g2]))
      [                   d               d            ]
      [     0        0    --(g(x, y))     --(g(x, y))  ]
      [                   dx              dy           ]
      [                                                ]
      [     0        0        2*x              3       ]
      [                                                ]
      [                     2               2          ]
      [d                   d               d           ]
      [--(g(x, y))  2*x   ---(f(x, y))   -----(f(x, y))]
      [dx                   2            dy dx         ]
      [                   dx                           ]
      [                                                ]
      [                     2               2          ]
      [d                   d               d           ]
      [--(g(x, y))   3   -----(f(x, y))   ---(f(x, y)) ]
      [dy                dy dx              2          ]
      [                                   dy           ]
  
      References
      ==========
  
      .. [1] https://en.wikipedia.org/wiki/Hessian_matrix
  
      See Also
      ========
  
      sympy.matrices.matrixbase.MatrixBase.jacobian
      wronskian
      """
      # f is the expression representing a function f, return regular matrix
      if isinstance(varlist, MatrixBase):
          if 1 not in varlist.shape:
              raise ShapeError(
                  "`varlist` must be a column or row vector."
              )
          if varlist.cols == 1:
              varlist = varlist.T
          varlist = varlist.tolist()[0]
      if is_sequence(varlist):
          n = len(varlist)
          if not n:
              raise ShapeError("`len(varlist)` must not be zero.")
      else:
          raise ValueError("Improper variable list in hessian function")
      if not getattr(f, "diff"):
          # check differentiability
          raise ValueError(
              "Function `f` (%s) is not differentiable" % f
          )
      m = len(constraints)
      N = m + n
      out = zeros(N)
      for k, g in enumerate(constraints):
          if not getattr(g, "diff"):
              # check differentiability
              raise ValueError(
                  "Function `f` (%s) is not differentiable" % f
              )
          for i in range(n):
              out[k, i + m] = g.diff(varlist[i])
      for i in range(n):
          for j in range(i, n):
              out[i + m, j + m] = f.diff(varlist[i]).diff(varlist[j])
      for i in range(N):
          for j in range(i + 1, N):
              out[j, i] = out[i, j]
      return out
  
imports: import random
  from sympy.core.basic import Basic
  from sympy.core.singleton import S
  from sympy.core.symbol import Symbol
  from sympy.core.sympify import sympify
  from sympy.functions.elementary.trigonometric import cos, sin
  from sympy.utilities.decorator import doctest_depends_on
  from sympy.utilities.exceptions import sympy_deprecation_warning
  from sympy.utilities.iterables import is_sequence
  from .exceptions import ShapeError
  from .decompositions import _cholesky, _LDLdecomposition
  from .matrixbase import MatrixBase
  from .repmatrix import MutableRepMatrix, RepMatrix
  from .solvers import _lower_triangular_solve, _upper_triangular_solve
label_desc: Use image node_id 17 for calling a global function with example usage: hessian(f, varlist, constraints) and returns: out