Multivariate power-law fit to the smoothing width based on number of iterations and nearest-neighbor weighting coefficient

fit_smoothing_width.py

#!/usr/bin/env python

import numpy
import matplotlib.pyplot as plt
from scipy.ndimage.filters import convolve1d
from pandas import DataFrame
from sklearn import linear_model

nx = 2001
x = numpy.linspace(-10., 10., nx)
dx = x[1] - x[0]

f = numpy.zeros((nx,))
half = (nx-1)//2
f[half] = 1.

outers = numpy.linspace(-2., -0.5, 40)
iterations = 1000
iters = numpy.arange(1, iterations + 1)

Outers, Iters = numpy.meshgrid(outers, iters)

maxes = numpy.zeros(Outers.shape)
widths = numpy.zeros(Outers.shape)

for outerIndex, outer in enumerate(outers):
    weights = numpy.array([10**outer, 1., 10**outer])
    weights = weights/weights.sum()

    fsmooth = f
    for iter in range(iterations):
        fsmooth = convolve1d(fsmooth, weights, mode='wrap')
        fmax = numpy.amax(fsmooth)
        width = 2.*numpy.interp(0.5*fmax, fsmooth[-1:half:-1], x[-1:half:-1])
        maxes[iter, outerIndex] = fmax
        widths[iter, outerIndex] = width/dx

mask = widths > 3.0

data_dict = {'weight': Outers[mask],
             'iteration': numpy.log10(Iters[mask]),
             'width': numpy.log10(widths[mask])}

df = DataFrame(data_dict, columns=['weight', 'iteration', 'width'])

X = df[['weight', 'iteration']]
Y = df['width']

# with sklearn
regr = linear_model.LinearRegression()
regr.fit(X, Y)

print('Intercept: \n', regr.intercept_, 10**regr.intercept_)
print('Coefficients: \n', regr.coef_)

#A = 10**regr.intercept_
# p_weight = regr.coef_[0]
# p_iter = regr.coef_[1]
# width_fit = A*10**(Outers*p_weight)*Iters**p_iter

p_weight = 0.45
p_iter = 0.5
intercept = numpy.mean(data_dict['width'] - (p_weight*data_dict['weight'] +
                                             p_iter*data_dict['iteration']))
A = 10**intercept
width_fit = A*10**(Outers*p_weight)*Iters**p_iter

print('A: {}'.format(A))

# A = 2.65
# p_iter = 0.5
# p_width = 2.25
# weight = (width/(A*sqrt(iterations)))**p_width

vmin = numpy.amin(widths)
vmax = numpy.amax(widths)

plt.figure(1)
plt.pcolormesh(Outers, Iters, widths, vmin=vmin, vmax=vmax)
plt.colorbar()

plt.figure(2)
plt.pcolormesh(Outers, Iters, width_fit, vmin=vmin, vmax=vmax)
plt.colorbar()


plt.figure(3)
plt.pcolormesh(Outers, Iters, widths - width_fit)
plt.colorbar()

plt.figure(4)
plt.pcolormesh(Outers, Iters, mask)

plt.show()