code: def find_MAP(
      start=None,
      vars: Optional[Sequence[Variable]] = None,
      method="L-BFGS-B",
      return_raw=False,
      include_transformed=True,
      progressbar=True,
      maxeval=5000,
      model=None,
      *args,
      seed: Optional[int] = None,
      **kwargs,
  ):
      """Finds the local maximum a posteriori point given a model.
  
      `find_MAP` should not be used to initialize the NUTS sampler. Simply call
      ``pymc.sample()`` and it will automatically initialize NUTS in a better
      way.
  
      Parameters
      ----------
      start: `dict` of parameter values (Defaults to `model.initial_point`)
          These values will be fixed and used for any free RandomVariables that are
          not being optimized.
      vars: list of TensorVariable
          List of free RandomVariables to optimize the posterior with respect to.
          Defaults to all continuous RVs in a model. The respective value variables
          may also be passed instead.
      method: string or callable, optional
          Optimization algorithm. Defaults to 'L-BFGS-B' unless discrete variables are
          specified in `vars`, then `Powell` which will perform better. For instructions
          on use of a callable, refer to SciPy's documentation of `optimize.minimize`.
      return_raw: bool, optional defaults to False
          Whether to return the full output of scipy.optimize.minimize
      include_transformed: bool, optional defaults to True
          Flag for reporting automatically unconstrained transformed values in addition
          to the constrained values
      progressbar: bool, optional defaults to True
          Whether to display a progress bar in the command line.
      maxeval: int, optional, defaults to 5000
          The maximum number of times the posterior distribution is evaluated.
      model: Model (optional if in `with` context)
      *args, **kwargs
          Extra args passed to scipy.optimize.minimize
  
      Notes
      -----
      Older code examples used `find_MAP` to initialize the NUTS sampler,
      but this is not an effective way of choosing starting values for sampling.
      As a result, we have greatly enhanced the initialization of NUTS and
      wrapped it inside ``pymc.sample()`` and you should thus avoid this method.
      """
      model = modelcontext(model)
  
      if vars is None:
          vars = model.continuous_value_vars
          if not vars:
              raise ValueError(
                  "Model has no unobserved continuous variables."
              )
      else:
          try:
              vars = get_value_vars_from_user_vars(vars, model)
          except ValueError as exc:
              # Accommodate case where user passed non-pure RV nodes
              vars = pm.inputvars(model.replace_rvs_by_values(vars))
              if vars:
                  warnings.warn(
                      "Intermediate variables (such as Deterministic or Potential) were passed. "
                      "find_MAP will optimize the underlying free_RVs instead.",
                      UserWarning,
                  )
              else:
                  raise exc
  
      disc_vars = list(typefilter(vars, discrete_types))
      ipfn = make_initial_point_fn(
          model=model,
          jitter_rvs=set(),
          return_transformed=True,
          overrides=start,
      )
      start = ipfn(seed)
      model.check_start_vals(start)
  
      vars_dict = {var.name: var for var in vars}
      x0 = DictToArrayBijection.map(
          {
              var_name: value
              for var_name, value in start.items()
              if var_name in vars_dict
          }
      )
  
      # TODO: If the mapping is fixed, we can simply create graphs for the
      # mapping and avoid all this bijection overhead
      compiled_logp_func = DictToArrayBijection.mapf(
          model.compile_logp(jacobian=False), start
      )
      logp_func = lambda x: compiled_logp_func(
          RaveledVars(x, x0.point_map_info)
      )  # noqa E731
  
      rvs = [
          model.values_to_rvs[vars_dict[name]]
          for name, _, _ in x0.point_map_info
      ]
      try:
          # This might be needed for calls to `dlogp_func`
          # start_map_info = tuple((v.name, v.shape, v.dtype) for v in vars)
          compiled_dlogp_func = DictToArrayBijection.mapf(
              model.compile_dlogp(rvs, jacobian=False), start
          )
          dlogp_func = lambda x: compiled_dlogp_func(
              RaveledVars(x, x0.point_map_info)
          )  # noqa E731
          compute_gradient = True
      except (
          AttributeError,
          NotImplementedError,
          tg.NullTypeGradError,
      ):
          compute_gradient = False
  
      if disc_vars or not compute_gradient:
          pm._log.warning(
              "Warning: gradient not available."
              + "(E.g. vars contains discrete variables). MAP "
              + "estimates may not be accurate for the default "
              + "parameters. Defaulting to non-gradient minimization "
              + "'Powell'."
          )
          method = "Powell"
  
      if compute_gradient and method != "Powell":
          cost_func = CostFuncWrapper(
              maxeval, progressbar, logp_func, dlogp_func
          )
      else:
          cost_func = CostFuncWrapper(maxeval, progressbar, logp_func)
          compute_gradient = False
  
      try:
          opt_result = minimize(
              cost_func,
              x0.data,
              method=method,
              jac=compute_gradient,
              *args,
              **kwargs,
          )
          mx0 = opt_result["x"]  # r -> opt_result
      except (KeyboardInterrupt, StopIteration) as e:
          mx0, opt_result = cost_func.previous_x, None
          if isinstance(e, StopIteration):
              pm._log.info(e)
      finally:
          last_v = cost_func.n_eval
          if progressbar:
              assert isinstance(cost_func.progress, ProgressBar)
              cost_func.progress.total = last_v
              cost_func.progress.update(last_v)
              print(file=sys.stdout)
  
      mx0 = RaveledVars(mx0, x0.point_map_info)
      unobserved_vars = get_default_varnames(
          model.unobserved_value_vars, include_transformed
      )
      unobserved_vars_values = model.compile_fn(unobserved_vars)(
          DictToArrayBijection.rmap(mx0, start)
      )
      mx = {
          var.name: value
          for var, value in zip(unobserved_vars, unobserved_vars_values)
      }
  
      if return_raw:
          return mx, opt_result
      else:
          return mx
  
imports: import sys
  import warnings
  from typing import Optional, Sequence
  import numpy
  import pytensor.gradient
  from fastprogress.fastprogress import ProgressBar, progress_bar
  from numpy import isfinite
  from pytensor import Variable
  from scipy.optimize import minimize
  import pymc
  from pymc.blocking import DictToArrayBijection, RaveledVars
  from pymc.initial_point import make_initial_point_fn
  from pymc.model import modelcontext
  from pymc.util import get_default_varnames, get_value_vars_from_user_vars
  from pymc.vartypes import discrete_types, typefilter
label_desc: Use image node_id 1 for calling a global function with example usage: find_MAP(start, vars, method, return_raw, include_transformed, progressbar, maxeval, model) and returns: mx,  opt_result, mx