code: def _cost_matrix(b_m, hdmr):
      """The cost matrix stores information about hierarchical orthogonality.
      It is structured in a way that ensures orthogonality between component
      functions that are hierarchically related.
  
      Parameters
      ----------
      b_m : numpy.array
          Basis matrix
      hdmr : namedtuple
          Core parameters of hdmr expansion
  
      Returns
      -------
      cost : numpy.array
          Cost matrix
      """
      cost = np.zeros((hdmr.a_tnt, hdmr.a_tnt))
  
      range_2nd_1 = lambda x: range(
          hdmr.tnt1 + (x) * hdmr.nt2, hdmr.tnt1 + (x + 1) * hdmr.nt2
      )
      range_2nd_2 = lambda x: range(
          hdmr.tnt1 + (x) * hdmr.nt2,
          hdmr.tnt1 + (x) * hdmr.nt2 + hdmr.p_o * 2,
      )
      range_3rd_1 = lambda x: range(
          hdmr.tnt1 + hdmr.tnt2 + (x) * hdmr.nt3,
          hdmr.tnt1 + hdmr.tnt2 + (x + 1) * hdmr.nt3,
      )
      range_3rd_2 = lambda x: range(
          hdmr.tnt1 + hdmr.tnt2 + (x) * hdmr.nt3,
          hdmr.tnt1
          + hdmr.tnt2
          + (x) * hdmr.nt3
          + 3 * hdmr.p_o
          + 3 * hdmr.p_o**2,
      )
  
      if hdmr.max_order > 1:
          sr_i = np.mean(b_m, axis=0, keepdims=True)
          sr_ij = np.zeros((2 * hdmr.p_o + 1, hdmr.nt2))
          ct = 0
          for _ in _prod(range(0, hdmr.d - 1), range(1, hdmr.d)):
              sr_ij[0, :] = sr_i[0, range_2nd_1(ct)]
              sr_ij[1:, :] = (
                  b_m[:, range_2nd_2(ct)].T @ b_m[:, range_2nd_1(ct)]
              ) / hdmr.subset
              cost[np.ix_(range_2nd_1(ct), range_2nd_1(ct))] = (
                  sr_ij.T @ sr_ij
              )
              ct += 1
      if hdmr.max_order == 3:
          sr_ijk = np.zeros(
              (3 * hdmr.p_o + 3 * hdmr.p_o**2 + 1, hdmr.nt3)
          )
          ct = 0
          for _ in _prod(
              range(0, hdmr.d - 2),
              range(1, hdmr.d - 1),
              range(2, hdmr.d),
          ):
              sr_ijk[0, :] = sr_i[0, range_3rd_1(ct)]
              sr_ijk[1:, :] = (
                  b_m[:, range_3rd_2(ct)].T
                  @ b_m[:, range_3rd_1(ct)]
                  / hdmr.subset
              )
              cost[np.ix_(range_3rd_1(ct), range_3rd_1(ct))] = (
                  sr_ijk.T @ sr_ijk
              )
              ct += 1
  
      return cost
  
imports: import math
  import warnings
  from typing import Dict, Tuple
  from types import MethodType
  from itertools import combinations, product
  from collections import defaultdict, namedtuple
  import numpy
  from pandas import DataFrame
  from numpy.linalg import det, pinv, matrix_rank
  from scipy.linalg import svd, LinAlgError, solve
  from scipy import stats, special
  from .None import common_args
  from ..util import read_param_file, ResultDict
label_desc: Use image node_id 8 for calling a global function with example usage: _cost_matrix(b_m, hdmr) and returns: cost