compare scipy.interpolate.griddata with Intergrid wrapper for scipy.ndimage.map_coordinates

griddata-intergrid-div5.log

# from: griddata-intergrid.py
# run: 29 Jun 2013 19:09  in ~bz/py/interpol/git/intergrid/test  i486  mac 10.8.3

versions: numpy 1.7.1  scipy 0.12.0  python 2.7.4
--------------------------------------------------------------------------------
griddata vs. intergrid:  Nx 36 Ny 36 Nz 84  subsample 2  order 1  prefilter 0
func min av max: -1 0.0025 1
griddata: 121.1 sec
Intergrid: 2.82 msec  13608 points in a (36, 36, 84) grid  3 maps  order 1
av |f - griddata|: 0.44
av |f - intergrid|: 0.44
123.618u 0.371s 2:03.99 99.9%

griddata-intergrid.py

# griddata-intergrid.py 29jun
# from unutbu https://gist.github.com/anonymous/5879738
# http://stackoverflow.com/questions/15748767/interpolation-subsampling-of-3d-data-in-python-without-vtk
# div 5: griddata 120 sec, Intergrid 2.8 msec, same avdiff ??

from __future__ import division
import sys
from time import time
import numpy as np
from scipy import interpolate, __version__
    # http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.griddata.html
# import matplotlib.pyplot as plt
import intergrid  # http://denis-bz.github.com/docs/intergrid.html


print "versions: numpy %s  scipy %s  python %s" % (
    np.__version__, __version__, sys.version.split()[0] )


def func( X, ncycle=10 ):
    r = np.sqrt( np.sum( X**2, axis=-1 ))
    return np.sin( 2 * np.pi * ncycle * r )

#...............................................................................
# generators for box grids, uniform / nonuniform --

def intgrid( *dims, **kw ):
    """ intgrid( 2, 3, 4 ) -> 24 x 3 [[0 0 0] [0 0 1] ... [1 2 3]] """
    dtype = kw.get( "dtype", int )
    return np.indices( dims, dtype=dtype ) .reshape(( len(dims), -1 )) .T

def cubegrid( *dims, **kw ):
    """ cubegrid( 2, 5 ) -> 2 x 5 points in the unit cube
        box centres: [[.25 .1] [.25 .3] [.25 .5] ... [.75 .9]]
    """
    dtype = kw.get( "dtype", np.float64 )
    return (intgrid( *dims, dtype=dtype ) + .5) / np.array( dims )

def productgrid( *xyz ):
    """ productgrid( 1d arrays x, y, z )
        -> [[x0 y0 z0] [x0 y0 z1] ... [xlast ylast zlast]]
    """
        # meshgrid mgrid ogrid ? mttiw
    from itertools import product
    return np.array( list( product( *xyz )))

def func_grid( func, points, **kw ):
    return np.array([ func( p, **kw ) for p in points ])  # faster way ?
    # unify fromfunction, func( ogrid ) ... ?

def avdiff( x, y ):
    return np.fabs( x - y ).mean()

#...............................................................................
Nx, Ny, Nz = 181, 181, 421
div = 5
subsample = 2
order = 1
prefilter = 0  # 0: smoothing B-spline, 1/3: M-N, 1: interpolating 
plot = 0
seed = 0

    # run this.py a=1 b=None 'c = expr' ...  from sh or ipython
exec( "\n".join( sys.argv[1:] ))
np.set_printoptions( 2, threshold=100, edgeitems=10, suppress=True )
np.random.seed(seed)

Nx //= div
Ny //= div
Nz //= div

title = "griddata vs. intergrid:  Nx %d Ny %d Nz %d  subsample %d  order %d  prefilter %g" % (
    Nx, Ny, Nz, subsample, order, prefilter )
print 80 * "-"
print title

    # Define a random box grid of shape (Nx * Ny * Nz, 3)
    # (0 1 so griddata inside convex hull, else NaNs)
x = np.r_[ 0, np.sort( np.random.random( Nx - 2 )), 1 ]
y = np.r_[ 0, np.sort( np.random.random( Ny - 2 )), 1 ]
z = np.r_[ 0, np.sort( np.random.random( Nz - 2 )), 1 ]
randomgrid = productgrid( x, y, z )

    # sample grid: uniform (Mx * My * Mz, 3)
Mx, My, Mz = Nx // subsample, Ny // subsample, Nz // subsample
samplegrid = cubegrid( Mx, My, Mz )

#...............................................................................
D = func_grid( func, randomgrid )  # 1d, Nx * Ny * Nz
print "func min av max: %.2g %.2g %.2g" % ( D.min(), D.mean(), D.max() )

t0 = time()
D_interpolated = interpolate.griddata( randomgrid, D, samplegrid ) .reshape(Mx, My, Mz)
    # docs.scipy griddata says xi (M, ndim) but the example uses mgrid ?
print "griddata: %.1f sec" % (time() - t0)

interfunc = intergrid.Intergrid(
    D.reshape( Nx, Ny, Nz ),
    lo=0, hi=1,
    maps=[x, y, z],
    order=order, prefilter=prefilter, verbose=1 )
intergridded = interfunc( samplegrid ) .reshape(Mx, My, Mz)

func_exact = func_grid( func, samplegrid ) .reshape(Mx, My, Mz)
print "av |f - griddata|: %.2g" % avdiff( func_exact, D_interpolated )
print "av |f - intergrid|: %.2g" % avdiff( func_exact, intergridded )