docs: def classical(vulnerability_function, hazard_imls, hazard_poes, loss_ratios):
      """
      :param vulnerability_function:
          an instance of
          :py:class:`openquake.risklib.scientific.VulnerabilityFunction`
          representing the vulnerability function used to compute the curve.
      :param hazard_imls:
          the hazard intensity measure type and levels
      :type hazard_poes:
          the hazard curve
      :param loss_ratios:
          a tuple of C loss ratios
      :returns:
          an array of shape (2, C)
      """
      assert len(hazard_imls) == len(hazard_poes), (
          len(hazard_imls), len(hazard_poes))
      vf = vulnerability_function
      imls = vf.mean_imls()
      lrem = vf.loss_ratio_exceedance_matrix(loss_ratios)
  
      # saturate imls to hazard imls
      min_val, max_val = hazard_imls[0], hazard_imls[-1]
      numpy.putmask(imls, imls < min_val, min_val)
      numpy.putmask(imls, imls > max_val, max_val)
  
      # interpolate the hazard curve
      poes = interpolate.interp1d(hazard_imls, hazard_poes)(imls)
  
      # compute the poos
      pos = pairwise_diff(poes)
      lrem_po = numpy.empty(lrem.shape)
      for idx, po in enumerate(pos):
          lrem_po[:, idx] = lrem[:, idx] * po  # column * po
      return numpy.array([loss_ratios, lrem_po.sum(axis=1)])