nhlib performance measure tools

perfmeasure.py

import sys
from decimal import Decimal
import time
import multiprocessing

import numpy

from nhlib import const
from nhlib.site import SiteCollection
from nhlib.source import AreaSource
from nhlib.pmf import PMF
from nhlib.geo import NodalPlane
from nhlib.scalerel import PeerMSR
from nhlib.gsim.sadigh_1997 import SadighEtAl1997
from nhlib.calc import hazard_curves_poissonian as hazard_curves

from tests.acceptance import _peer_test_data as test_data


def run(area_discretization, num_sites, gmpe):
    hypocenter_probability = (
        Decimal(1) / len(test_data.SET1_CASE11_HYPOCENTERS)
    )
    hypocenter_pmf = PMF([
        (hypocenter_probability, hypocenter)
        for hypocenter in test_data.SET1_CASE11_HYPOCENTERS
    ])
    sources = [AreaSource(source_id='area', name='area',
        tectonic_region_type=const.TRT.ACTIVE_SHALLOW_CRUST,
        mfd=test_data.SET1_CASE11_MFD,
        nodal_plane_distribution=PMF([(1, NodalPlane(0.0, 90.0, 0.0))]),
        hypocenter_distribution=hypocenter_pmf,
        upper_seismogenic_depth=0.0,
        lower_seismogenic_depth=10.0,
        magnitude_scaling_relationship = PeerMSR(),
        rupture_aspect_ratio=test_data.SET1_RUPTURE_ASPECT_RATIO,
        polygon=test_data.SET1_CASE11_SOURCE_POLYGON,
        area_discretization=area_discretization,
        rupture_mesh_spacing=1.0
    )]
    sites = SiteCollection([test_data.SET1_CASE11_SITE1] * num_sites)
    gsims = {const.TRT.ACTIVE_SHALLOW_CRUST: gmpe}
    truncation_level = 1.0
    time_span = 1.0
    imts = {test_data.IMT: test_data.SET1_CASE11_IMLS}
    start = time.time()
    hazard_curves(sources, sites, imts, time_span,
                  gsims, truncation_level)
    stop = time.time()
    return stop - start


class Runner(multiprocessing.Process):
    def __init__(self, tasks_queue, res_queue):
        super(Runner, self).__init__()
        self.tasks_queue = tasks_queue
        self.res_queue = res_queue

    def run(self):
        while True:
            task = self.tasks_queue.get()
            if task is None:
                break
            discr, num_sites, num_ruptures, gmpe = task
            gmpe_name = type(gmpe).__name__
            time_spent = run(discr, num_sites, gmpe)
            print >> sys.stderr, \
                  "GMPE=%s, discretization=%d, num sites=%d, " \
                  "num ruptures=%d: %.1f sec spent" % \
                  (gmpe_name, discr, num_sites, num_ruptures, time_spent)
            self.res_queue.put((num_sites, num_ruptures, time_spent))
            self.tasks_queue.task_done()


NUM_PROCESSES = 4

def measure():
    GMPE = SadighEtAl1997()
    sbands, smax = 17, 15000
    spowers = numpy.empty(sbands)
    spowers.fill(numpy.log10(smax) / sbands)
    NUM_SITES = numpy.unique((10 ** spowers.cumsum()).astype(int))
    # mapping of area discretizations to num ruptures
    NUM_RUPTURES = {#10: 27810, 12: 19260, 14: 14130, 17: 9540, 20: 6660,
                    23: 5130, 26: 4050, 30: 2970, 35: 2070, 40: 1440}

    tasks_queue = multiprocessing.JoinableQueue()
    res_queue = multiprocessing.Queue()
    processes = [Runner(tasks_queue, res_queue) for i in xrange(NUM_PROCESSES)]
    for process in processes:
        process.start()

    num_tasks = 0
    for num_sites in sorted(NUM_SITES, reverse=True):
        for discr in sorted(NUM_RUPTURES):
            num_tasks += 1
            tasks_queue.put((discr, num_sites, NUM_RUPTURES[discr], GMPE))
    for process in processes:
        tasks_queue.put(None)

    results = []
    while len(results) < num_tasks:
        num_sites, num_ruptures, time_spent = res_queue.get()
        results.append((num_sites, num_ruptures, time_spent))
    numpy.save(sys.stdout, numpy.array(results))


if __name__ == '__main__':
    measure()

perfplot.py

import sys
import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D; Axes3D

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

data = numpy.array(sorted([tuple(row) for row in numpy.load(sys.stdin)]))
sites, ruptures, times = data.transpose()
sites = numpy.unique(sites.astype(int))
ruptures = numpy.unique(ruptures.astype(int))
times = times.reshape((len(sites), len(ruptures)))
x, y = numpy.meshgrid(numpy.log10(ruptures), numpy.log10(sites))
z = numpy.log10(times)

surf = ax.plot_wireframe(x, y, z)

ax.set_ylabel('sites (powers of 10)')
ax.set_xlabel('ruptures (powers of 10)')
ax.set_zlabel('time to compute (seconds, powers of 10)')

plt.show()