code: def median_squared_scaled_error(
      y_true,
      y_pred,
      sp=1,
      horizon_weight=None,
      multioutput="uniform_average",
      square_root=False,
      **kwargs,
  ):
      """Median squared scaled error (MdSSE) or root median squared scaled error (RMdSSE).
  
      If `square_root` is False then calculates MdSSE, otherwise calculates RMdSSE if
      `square_root` is True. Both MdSSE and RMdSSE output is non-negative floating
      point. The best value is 0.0.
  
      This is a squared variant of the MdASE loss metric. Like MASE and other
      scaled performance metrics this scale-free metric can be used to compare
      forecast methods on a single series or between series.
  
      This metric is also suited for intermittent-demand series because it
      will not give infinite or undefined values unless the training data
      is a flat timeseries. In this case the function returns a large value
      instead of inf.
  
      Works with multioutput (multivariate) timeseries data
      with homogeneous seasonal periodicity.
  
      Parameters
      ----------
      y_true : pd.Series, pd.DataFrame or np.array of shape (fh,) or (fh, n_outputs) \
               where fh is the forecasting horizon
          Ground truth (correct) target values.
      y_pred : pd.Series, pd.DataFrame or np.array of shape (fh,) or (fh, n_outputs) \
               where fh is the forecasting horizon
          Forecasted values.
      y_train : pd.Series, pd.DataFrame or np.array of shape (n_timepoints,) or \
               (n_timepoints, n_outputs), default = None
          Observed training values.
      sp : int
          Seasonal periodicity of training data.
      horizon_weight : array-like of shape (fh,), default=None
          Forecast horizon weights.
      multioutput : {'raw_values', 'uniform_average'}  or array-like of shape \
              (n_outputs,), default='uniform_average'
          Defines how to aggregate metric for multivariate (multioutput) data.
          If array-like, values used as weights to average the errors.
          If 'raw_values', returns a full set of errors in case of multioutput input.
          If 'uniform_average', errors of all outputs are averaged with uniform weight.
  
      Returns
      -------
      loss : float
          RMdSSE loss.
          If multioutput is 'raw_values', then RMdSSE is returned for each
          output separately.
          If multioutput is 'uniform_average' or an ndarray of weights, then the
          weighted average RMdSSE of all output errors is returned.
  
      See Also
      --------
      mean_absolute_scaled_error
      median_absolute_scaled_error
      mean_squared_scaled_error
  
      References
      ----------
      M5 Competition Guidelines.
      https://mofc.unic.ac.cy/wp-content/uploads/2020/03/M5-Competitors-Guide-Final-10-March-2020.docx
  
      Hyndman, R. J and Koehler, A. B. (2006). "Another look at measures of
      forecast accuracy", International Journal of Forecasting, Volume 22, Issue 4.
  
      Examples
      --------
      >>> from sktime.performance_metrics.forecasting import median_squared_scaled_error
      >>> y_train = np.array([5, 0.5, 4, 6, 3, 5, 2])
      >>> y_true = np.array([3, -0.5, 2, 7, 2])
      >>> y_pred = np.array([2.5, 0.0, 2, 8, 1.25])
      >>> median_squared_scaled_error(y_true, y_pred, y_train=y_train, square_root=True)
      0.16666666666666666
      >>> y_train = np.array([[0.5, 1], [-1, 1], [7, -6]])
      >>> y_true = np.array([[0.5, 1], [-1, 1], [7, -6]])
      >>> y_pred = np.array([[0, 2], [-1, 2], [8, -5]])
      >>> median_squared_scaled_error(y_true, y_pred, y_train=y_train, square_root=True)
      0.1472819539849714
      >>> median_squared_scaled_error(y_true, y_pred, y_train=y_train, \
      multioutput='raw_values', square_root=True)
      array([0.08687445, 0.20203051])
      >>> median_squared_scaled_error(y_true, y_pred, y_train=y_train, \
      multioutput=[0.3, 0.7], square_root=True)
      0.16914781383660782
      """
      y_train = _get_kwarg(
          "y_train", metric_name="median_squared_scaled_error", **kwargs
      )
  
      # Other input checks
      _, y_true, y_pred, multioutput = _check_reg_targets(
          y_true, y_pred, multioutput
      )
      if horizon_weight is not None:
          check_consistent_length(y_true, horizon_weight)
  
      if y_train.ndim == 1:
          y_train = np.expand_dims(y_train, 1)
  
      # Check test and train have same dimensions
      if y_true.ndim != y_train.ndim:
          raise ValueError(
              "Equal dimension required for y_true and y_train"
          )
  
      if (y_true.ndim > 1) and (y_true.shape[1] != y_train.shape[1]):
          raise ValueError(
              "Equal number of columns required for y_true and y_train"
          )
  
      # naive seasonal prediction
      y_train = np.asarray(y_train)
      y_pred_naive = y_train[:-sp]
  
      # median squared error of naive seasonal prediction
      mdse_naive = median_squared_error(
          y_train[sp:], y_pred_naive, multioutput=multioutput
      )
  
      mdse = median_squared_error(
          y_true,
          y_pred,
          horizon_weight=horizon_weight,
          multioutput=multioutput,
      )
  
      if square_root:
          loss = np.sqrt(mdse / np.maximum(mdse_naive, EPS))
      else:
          loss = mdse / np.maximum(mdse_naive, EPS)
      return loss
  
imports: import numpy
  from scipy.stats import gmean
  from sklearn.metrics import mean_absolute_error
  from sklearn.metrics import mean_squared_error
  from sklearn.metrics import median_absolute_error
  from sklearn.metrics._regression import _check_reg_targets
  from sklearn.utils.stats import _weighted_percentile
  from sklearn.utils.validation import check_consistent_length
  from sktime.utils.stats import _weighted_geometric_mean
label_desc: Use image node_id 7 for calling a global function with example usage: median_squared_scaled_error(y_true, y_pred, sp, horizon_weight, multioutput, square_root) and returns: loss