code: def calculate_BODF(sub_network, skip_pre=False):
      
  
      if not skip_pre:
          calculate_PTDF(sub_network)
  
      num_branches = sub_network.PTDF.shape[0]
  
      #build LxL version of PTDF
      branch_PTDF = sub_network.PTDF*sub_network.K
  
      denominator = csr_matrix((1/(1-np.diag(branch_PTDF)),(r_[:num_branches],r_[:num_branches])))
  
      sub_network.BODF = branch_PTDF*denominator
  
      #make sure the flow on the branch itself is zero
      np.fill_diagonal(sub_network.BODF,-1)
docs: Calculate the Branch Outage Distribution Factor (BODF) for
      sub_network.
  
      Sets sub_network.BODF as a (dense) numpy array.
  
      The BODF is a num_branch x num_branch 2d array.
  
      For the outage of branch l, the new flow on branch k is
      given in terms of the flow before the outage
  
      f_k^after = f_k^before + BODF_{kl} f_l^before
  
      Note that BODF_{ll} = -1.
  
      Parameters
      ----------
      sub_network : pypsa.SubNetwork
      skip_pre: bool, default False
          Skip the preliminary step of computing the PTDF.