code: def batch_normalization_test(
      inputs,
      gamma,
      beta,
      mean,
      var,
      axes="per-activation",
      epsilon=1e-4,
  ):
      """
      Performs batch normalization of the given inputs, using the given mean and
      variance.
  
      Parameters
      ----------
      axes : 'per-activation', 'spatial' or a tuple of ints
          The axes along which the input should be normalized. ``'per-activation'``
          normalizes per activation and is equal to ``axes=(0,)``.
          ``'spatial'`` shares normalization factors across spatial dimensions
          (i.e., all dimensions past the second), which for 4D inputs would be
          equal to ``axes=(0, 2, 3)``.
      gamma : tensor
          Scale factors. The shape must match the shape of `inputs`,
          except for the axes in `axes`. These axes should be set to 1 or be
          skipped altogether (such that `gamma.ndim == inputs.ndim - len(axes)`).
      beta : tensor
          Biases. Must match the tensor layout of `gamma`.
      mean : tensor
          Means. Usually these are running averages computed during training.
          Must match the tensor layout of `gamma`.
      var : tensor
          Variances. Usually these are running averages computed during training.
          Must match the tensor layout of `gamma`.
      epsilon : float
          Epsilon value used in the batch normalization formula. Minimum allowed
          value is 1e-5 (imposed by cuDNN).
  
      Returns
      -------
      out : tensor
          Batch-normalized inputs.
  
      Notes
      -----
      If per-activation or spatial normalization is selected, this operation
      will use the cuDNN implementation. (This requires cuDNN 5 or newer.)
  
      The returned value is equivalent to:
  
      .. code-block:: python
  
          # for per-activation normalization
          axes = (0,)
          # for spatial normalization
          axes = (0,) + tuple(range(2, inputs.ndim))
          gamma, beta, mean, var = (T.addbroadcast(t, *axes)
                                    for t in (gamma, beta, mean, var))
          out = (inputs - mean) * gamma / T.sqrt(var + epsilon) + beta
      """
      ndim = inputs.ndim
      axes, non_bc_axes = _prepare_batch_normalization_axes(axes, ndim)
  
      # have the parameter tensors been broadcasted yet?
      if gamma.ndim == ndim:
          params_ndim = ndim
      else:
          params_ndim = len(non_bc_axes)
          params_dimshuffle_pattern = ["x"] * ndim
          for i, axis in enumerate(non_bc_axes):
              params_dimshuffle_pattern[axis] = i
  
      if gamma.ndim != params_ndim or beta.ndim != params_ndim:
          raise ValueError(
              "gamma and beta dimensionality must match the "
              "number of non-normalized axes, or have the "
              "same number of dimensions as the inputs; "
              "got %d and %d instead of %d"
              % (gamma.ndim, beta.ndim, params_ndim)
          )
      if mean.ndim != params_ndim or var.ndim != params_ndim:
          raise ValueError(
              "mean and var must be of the same dimensionality "
              "as gamma and beta; got %d and %d instead of %d"
              % (mean.ndim, var.ndim, params_ndim)
          )
  
      # epsilon will be converted to floatX later. we need to check
      # for rounding errors now, since numpy.float32(1e-5) < 1e-5.
      epsilon = np.cast[theano.config.floatX](epsilon)
      if epsilon < 1e-5:
          raise ValueError(
              "epsilon must be at least 1e-5, got %s" % str(epsilon)
          )
  
      gamma = as_tensor_variable(gamma)
      beta = as_tensor_variable(beta)
      mean = as_tensor_variable(mean)
      var = as_tensor_variable(var)
  
      if params_ndim != ndim:
          gamma = gamma.dimshuffle(params_dimshuffle_pattern)
          beta = beta.dimshuffle(params_dimshuffle_pattern)
          mean = mean.dimshuffle(params_dimshuffle_pattern)
          var = var.dimshuffle(params_dimshuffle_pattern)
      else:
          gamma = T.addbroadcast(gamma, *axes)
          beta = T.addbroadcast(beta, *axes)
          mean = T.addbroadcast(mean, *axes)
          var = T.addbroadcast(var, *axes)
  
      batchnorm_op = AbstractBatchNormInference(axes=axes)
      return batchnorm_op(
          inputs, gamma, beta, mean, var, epsilon=epsilon
      )
  
imports: from __future__ import absolute_import, print_function, division
  import numpy
  import theano
  from theano import Apply, Op
  from theano.gof import local_optimizer
  from theano.gof.opt import copy_stack_trace
  from theano.tensor import as_tensor_variable, TensorType
  from theano.tensor import basic
  from theano.tensor.opt import register_specialize_device
  from theano.scalar import Composite, as_common_dtype
  from theano.scalar import add, sub, true_div, mul
label_desc: Use image node_id 4 for calling a global function with example usage: batch_normalization_test(inputs, gamma, beta, mean, var, axes, epsilon) and returns: batchnorm_op