scikit-image PR #1792

generate_cached_times.py

# This was run from a notebook (so I could choose which cells to run).
# I have also uploaded the entire notebook as skimage_fftconvolve.ipynb (via a dropbox link)

# Cell 1
from skimage.restoration import richardson_lucy
import numpy as np
from scipy.signal import convolve, fftconvolve
import timeit
%matplotlib inline

# Cell 2
times = {}
for shape in [2, 3]:
    times[shape] = {}
    for N in Ns:
        times[shape][N] = {}
        Ks = np.logspace(0.5, np.log10(N), num=10, dtype=int)
        Ks = np.unique(Ks)

        print('## N = ' + str(N))
        for k in Ns:
            if k > N: continue
            print('k = ' + str(k))
            x = np.random.rand(*((N,)*shape))
            h = np.random.rand(*((k,)*shape))

            t_fft = time('fftconvolve', x, h)
            t_direct = time('convolve', x, h)
            times[shape][N][k] = {'t_fft':t_fft, 't_direct':t_direct}

# Cell 3   
import pickle
pickle.dump(times, open( "times.p", "wb"))

times.py

import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import pickle
import numpy as np
from scipy.optimize import curve_fit
sns.set_context('talk')
def process_shape(t_fft, t_direct, ndims=2):
    t_fft = pd.DataFrame.from_dict(t_fft)
    t_direct = pd.DataFrame.from_dict(t_direct)

    # showing the heatmap of times
    sns.heatmap(np.log2(t_fft / t_direct), center=0)
    plt.title('time_fft / time_direct for {}D array'.format(ndims))
    plt.xlabel('N \n(image square/cubic)')
    plt.ylabel('k \n(kernel square/cubic')

    # the size of the image/kernel
    N, k = np.meshgrid(1.0*t_fft.columns, 1.0*t_direct.index)

    # the variables in our equation
    A = (ndims*(N*np.log(N) + k*np.log(k)) / (k**ndims * N**ndims)).flat[:]
    b = (np.asarray(t_fft) / np.asarray(t_direct)).flat[:]

    # it's nan where we didn't measure (the lower-triangle of the matrix)
    i = np.isnan(b)
    b, A = b[~i], A[~i]

    # we assume solutions of this form; it's finding the constant out front
    def predictor(x, a):
        return a*x

    p = curve_fit(predictor, A, b, p0=1)
    print(p)

    plt.figure()
    plt.semilogy(predictor(A, p[0]), 'o-', label='Predicted ratio', basey=2)
    plt.semilogy(b, 'o-', label='Actual ratio', basey=2)
    plt.semilogy(b / b, label='Decision boundary', basey=2)
    plt.legend(loc='best')
    plt.title('time_fft / time_direct ratio for {}D array'.format(ndims))
    plt.xlabel('instance')
    plt.ylabel('ratio')
    plt.show()

# times saved as a dictionary of dictionary of the form
# times[N][k] = {'t_fft':1.0, 't_direct':10}
times = pickle.load(open("times.p", "rb"))
N_points = 7
Ns = np.logspace(0.5, 2.000, num=N_points, dtype=int)

from pprint import pprint
pprint(times)

for shape in [2, 3]:
    if shape==3:
        Ns = Ns[:5]
    t_fft = {N: {k: times[shape][N][k]['t_fft'] for k in Ns
                                                if k in times[shape][N].keys()}
                                                for N in Ns}
    t_direct = {N: {k: times[shape][N][k]['t_direct'] for k in Ns
                                                if k in times[shape][N].keys()}
                                                for N in Ns}

    process_shape(t_fft, t_direct, ndims=shape)

I expanded the scripts a bit to work with 1D as well, and made them work without IPython -- want a PR?