docs: def classical_damage(
          fragility_functions, hazard_imls, hazard_poes,
          investigation_time, risk_investigation_time):
      """
      :param fragility_functions:
          a list of fragility functions for each damage state
      :param hazard_imls:
          Intensity Measure Levels
      :param hazard_poes:
          hazard curve
      :param investigation_time:
          hazard investigation time
      :param risk_investigation_time:
          risk investigation time
      :returns:
          an array of M probabilities of occurrence where M is the numbers
          of damage states.
      """
      spi = fragility_functions.steps_per_interval
      if spi and spi > 1:  # interpolate
          imls = numpy.array(fragility_functions.interp_imls)
          min_val, max_val = hazard_imls[0], hazard_imls[-1]
          assert min_val > 0, hazard_imls  # sanity check
          numpy.putmask(imls, imls < min_val, min_val)
          numpy.putmask(imls, imls > max_val, max_val)
          poes = interpolate.interp1d(hazard_imls, hazard_poes)(imls)
      else:
          imls = (hazard_imls if fragility_functions.format == 'continuous'
                  else fragility_functions.imls)
          poes = numpy.array(hazard_poes)
      afe = annual_frequency_of_exceedence(poes, investigation_time)
      annual_frequency_of_occurrence = pairwise_diff(
          pairwise_mean([afe[0]] + list(afe) + [afe[-1]]))
      poes_per_damage_state = []
      for ff in fragility_functions:
          frequency_of_exceedence_per_damage_state = numpy.dot(
              annual_frequency_of_occurrence, list(map(ff, imls)))
          poe_per_damage_state = 1. - numpy.exp(
              - frequency_of_exceedence_per_damage_state *
              risk_investigation_time)
          poes_per_damage_state.append(poe_per_damage_state)
      poos = pairwise_diff([1] + poes_per_damage_state + [0])
      return poos