BenchMapCoordinates.py

# -*- coding: utf-8 -*-
""" 
Benchmark for spline interpolation of a 3D wind field, using the function map_coordinates. 
The spline interpolation is about 46000 times slower than a linear interpolation.
It is also about 10000 times slower than an equivalent program, written in the programming
language Julia.
"""
import time
import numpy as np
from scipy import ndimage

class Timer(object):
    def __init__(self, verbose=False):
        self.verbose = verbose

    def __enter__(self):
        self.start = time.time()
        return self

    def __exit__(self, *args):
        self.end = time.time()
        self.secs = self.end - self.start
        self.msecs = self.secs * 1000  # millisecs
        if self.verbose:
            print 'elapsed time: %f ms' % self.msecs

def benchMapCoordinates(u, v, w, x, y, z, **kwargs):
    x_wind = ndimage.map_coordinates(u, [[x], [y], [z]], **kwargs)
    y_wind = ndimage.map_coordinates(v, [[x], [y], [z]], **kwargs)
    z_wind = ndimage.map_coordinates(w, [[x], [y], [z]], **kwargs)
    return x_wind[0], y_wind[0], z_wind[0]

if __name__ == "__main__":
    u = np.random.rand(2026, 51, 251)
    v = np.random.rand(2026, 51, 251)
    w = np.random.rand(2026, 51, 251)

    with Timer() as t1:
        benchMapCoordinates(u, v, w, 1, 2, 3, order=1)
    print "time for benchMapCoordinates order = 1 [ms]:   ", t1.secs * 1e3 

    u, v, w = [ndimage.spline_filter(item, order=3) for item in [u, v, w]] 
    with Timer() as t1:
        benchMapCoordinates(u, v, w, 1, 2, 3, order=3, prefilter=False)
    print "time for benchMapCoordinates order = 3 [ms]:   ", t1.secs * 1e3