docs: def calc_accounts(S, L, Y, nr_sectors):
      """ Calculate sector specific cba and pba based accounts, imp and exp accounts
  
      The total industry output x for the calculation
      is recalculated from L and y
  
      Parameters
      ----------
      L : pandas.DataFrame
          Leontief input output table L
      S : pandas.DataFrame
          Direct impact coefficients
      Y : pandas.DataFrame
          Final demand: aggregated across categories or just one category, one
          column per country
      nr_sectors : int
          Number of sectors in the MRIO
  
  
      Returns
      -------
      Tuple
          (D_cba, D_pba, D_imp, D_exp)
  
          Format: D_row x L_col (=nr_countries*nr_sectors)
  
          - D_cba        Footprint per sector and country
          - D_pba      Total factur use per sector and country
          - D_imp       Total global factor use to satisfy total final demand in
                        the country per sector
          - D_exp       Total factor use in one country to satisfy final demand
                        in all other countries (per sector)
      """
      # diagonalize each sector block per country
      # this results in a disaggregated y with final demand per country per
      # sector in one column
      Y_diag = ioutil.diagonalize_blocks(Y.values, blocksize=nr_sectors)
      x_diag = L.dot(Y_diag)
      x_tot = x_diag.values.sum(1)
      del Y_diag
  
      D_cba = pd.DataFrame(S.values.dot(x_diag),
                           index=S.index,
                           columns=S.columns)
      # D_pba = S.dot(np.diagflat(x_tot))
      # faster broadcasted calculation:
      D_pba = pd.DataFrame(S.values*x_tot.reshape((1, -1)),
                           index=S.index,
                           columns=S.columns)
  
      # for the traded accounts set the domestic industry output to zero
      dom_block = np.zeros((nr_sectors, nr_sectors))
      x_trade = ioutil.set_block(x_diag.values, dom_block)
      D_imp = pd.DataFrame(S.values.dot(x_trade),
                           index=S.index,
                           columns=S.columns)
  
      x_exp = x_trade.sum(1)
      # D_exp = S.dot(np.diagflat(x_exp))
      # faster broadcasted version:
      D_exp = pd.DataFrame(S.values * x_exp.reshape((1, -1)),
                           index=S.index,
                           columns=S.columns)
  
      return (D_cba, D_pba, D_imp, D_exp)