Source code for openquake.hazardlib.tom
# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2012-2019 GEM Foundation
#
# OpenQuake is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# OpenQuake is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>.
"""
Module :mod:`openquake.hazardlib.tom` contains implementations of probability
density functions for earthquake temporal occurrence modeling.
"""
import abc
import math
import numpy
import scipy.stats
from openquake.baselib.slots import with_slots
registry = {}
[docs]@with_slots
class BaseTOM(metaclass=abc.ABCMeta):
"""
Base class for temporal occurrence model.
:param time_span:
The time interval of interest, in years.
:raises ValueError:
If ``time_span`` is not positive.
"""
_slots_ = ['time_span', 'occurrence_rate']
@classmethod
def __init_subclass__(cls):
registry[cls.__name__] = cls
def __init__(self, time_span, occurrence_rate=None):
if time_span <= 0:
raise ValueError('time_span must be positive')
self.time_span = time_span
self.occurrence_rate = occurrence_rate
[docs] @abc.abstractmethod
def get_probability_one_or_more_occurrences(self):
"""
Calculate and return the probability of event to happen one or more
times within the time range defined by constructor's ``time_span``
parameter value.
"""
[docs] @abc.abstractmethod
def get_probability_n_occurrences(self):
"""
Calculate the probability of occurrence of a number of events in the
constructor's ``time_span``.
"""
[docs] @abc.abstractmethod
def sample_number_of_occurrences(self, seeds=None):
"""
Draw a random sample from the distribution and return a number
of events to occur.
The method uses the numpy random generator, which needs a seed
in order to get reproducible results. If the seed is None, it
should be set outside of this method.
"""
[docs] @abc.abstractmethod
def get_probability_no_exceedance(self):
"""
Compute and return, for a number of ground motion levels and sites,
the probability that a rupture with annual occurrence rate given by
``occurrence_rate`` and able to cause ground motion values higher than
a given level at a site with probability ``poes``, does not cause any
exceedance in the time window specified by the ``time_span`` parameter
given in the constructor.
"""
[docs]@with_slots
class FatedTOM(BaseTOM):
def __init__(self, time_span, occurrence_rate=None):
self.time_span = time_span
self.occurrence_rate = occurrence_rate
[docs] def get_probability_one_or_more_occurrences(self, occurrence_rate):
return 1
[docs] def get_probability_n_occurrences(self, occurrence_rate, num):
if num != 1:
return 0
else:
return 1
[docs] def sample_number_of_occurrences(self, seeds=None):
return 1
[docs] def get_probability_no_exceedance(self, occurrence_rate, poes):
return 1-poes
[docs]@with_slots
class PoissonTOM(BaseTOM):
"""
Poissonian temporal occurrence model.
"""
[docs] def get_probability_one_or_more_occurrences(self, occurrence_rate):
"""
Calculates probability as ``1 - e ** (-occurrence_rate*time_span)``.
:param occurrence_rate:
The average number of events per year.
:return:
Float value between 0 and 1 inclusive.
"""
return 1 - math.exp(- occurrence_rate * self.time_span)
[docs] def get_probability_n_occurrences(self, occurrence_rate, num):
"""
Calculate the probability of occurrence of ``num`` events in the
constructor's ``time_span``.
:param occurrence_rate:
Annual rate of occurrence
:param num:
Number of events
:return:
Probability of occurrence
"""
return scipy.stats.poisson(occurrence_rate * self.time_span).pmf(num)
[docs] def sample_number_of_occurrences(self, occurrence_rate, seeds=None):
"""
Draw a random sample from the distribution and return a number
of events to occur.
The method uses the numpy random generator, which needs a seed
in order to get reproducible results. If the seed is None, it
should be set outside of this method.
:param occurrence_rate:
The average number of events per year.
:param seeds:
Random number generator seeds, one per each occurrence_rate
:return:
Sampled integer number of events to occur within model's
time span.
"""
if isinstance(seeds, numpy.ndarray): # array of seeds
assert len(seeds) == len(occurrence_rate), (
len(seeds), len(occurrence_rate))
rates = occurrence_rate * self.time_span
occ = []
for rate, seed in zip(rates, seeds):
numpy.random.seed(seed)
occ.append(numpy.random.poisson(rate))
return numpy.array(occ)
elif isinstance(seeds, int):
numpy.random.seed(seeds)
return numpy.random.poisson(occurrence_rate * self.time_span)
[docs] def get_probability_no_exceedance(self, occurrence_rate, poes):
"""
The probability is computed using the following formula ::
(1 - e ** (-occurrence_rate * time_span)) ** poes
:param occurrence_rate:
The average number of events per year.
:param poes:
2D numpy array containing conditional probabilities that the
rupture occurrence causes a ground shaking value exceeding a
ground motion level at a site. First dimension represent sites,
second dimension intensity measure levels. ``poes`` can be obtained
calling the :meth:`method
<openquake.hazardlib.gsim.base.GroundShakingIntensityModel.get_poes>`.
:return:
2D numpy array containing probabilities of no exceedance. First
dimension represents sites, second dimensions intensity measure
levels.
"""
p = self.get_probability_one_or_more_occurrences(occurrence_rate)
return (1 - p) ** poes