breakpoints

breakpoints

from scipy import interpolate

import numpy
import pylab
from pylab import *

from scipy.optimize import curve_fit


def func(x, a, b, c):
    return a * np.exp(-b * x) + c


i_dataset = numpy.asarray([0, 1.00, 86648.03, 18, 1.75, 101912.00, 24, 2.00, 114349.81, 42, 2.75, 143267.10, 48, 3.00,
                           147048.69, 66, 3.75, 400687.43, 72, 4.00, 524333.43])

i_dataset = i_dataset.reshape((7, 3))

i_dataset[:, 0] /= i_dataset[:, 0].max()
i_dataset[:, 2] /= i_dataset[:, 2].max()

RESAMPLING_STEP = 0.05

dataset = numpy.zeros((1 / RESAMPLING_STEP, 3))

f = interpolate.interp1d(i_dataset[:, 0], i_dataset[:, 2], kind='slinear')
dataset[:, 0] = np.arange(0, 1, RESAMPLING_STEP)
dataset[:, 2] = f(dataset[:, 0])

print dataset[:, 0]
print dataset[:, 2]

window_sizes = numpy.arange(6, 14)

print window_sizes

linear_rms = zeros((window_sizes.shape[0], dataset.shape[0] - window_sizes.min() + 1))
exponential_rms = zeros((window_sizes.shape[0], dataset.shape[0] - window_sizes.min() + 1))

linear_rms -= Inf
exponential_rms -= Inf

for i in range(0, window_sizes.shape[0]):
    for j in range(0, dataset.shape[0] - window_sizes[i] + 1):
        position = numpy.arange(j, j + window_sizes[i])

        linear_fit = polyfit(dataset[position, 0], dataset[position, 2], 1)
        linear_curve = polyval(linear_fit, dataset[position, 0])
        linear_rms[i, j] = sqrt(mean((dataset[position, 2] - linear_curve) ** 2))

        try:
            exponential_fit, cov = curve_fit(func, dataset[position, 0], dataset[position, 2], maxfev=10000)
            exponential_curve = func(dataset[position, 0], exponential_fit[0], exponential_fit[1], exponential_fit[2])
            exponential_rms[i, j] = sqrt(mean((dataset[position, 2] - exponential_curve) ** 2))
            # f = pylab.figure()
            # f.add_subplot(2, 2, 1)
            # pylab.plot(dataset[position, 0], dataset[position, 2], '-ob')
            # pylab.ylim([0, 1])
            # pylab.title('Dataset')
            # f.add_subplot(2, 2, 2)
            # pylab.plot(dataset[position, 0], linear_curve, '-ob')
            # pylab.ylim([0, 1])
            # pylab.title('Linear fit')
            # f.add_subplot(2, 2, 3)
            # pylab.plot(dataset[position, 0], exponential_curve, '-ob')
            # pylab.ylim([0, 1])
            # pylab.title('Exponential fit')
            # pylab.show()
        except RuntimeError:
            print 'err'

# TODO: make use of RMS values
print linear_rms
print exponential_rms

v_max = linear_rms.max() if linear_rms.max() > exponential_rms.max() else exponential_rms.max()

f = pylab.figure()

f.add_subplot(2, 3, 1)
pylab.plot(i_dataset[:, 0], i_dataset[:, 2], '-ob')
pylab.ylim([0, 1])
pylab.title('Initial Dataset')

f.add_subplot(2, 3, 2)
pylab.plot(dataset[:, 0], dataset[:, 2], '-ob')
pylab.ylim([0, 1])
pylab.title('Resampled Dataset')

ax = f.add_subplot(2, 3, 3)
ax.set_ylabel('Window Sizes')
ax.yaxis.set_ticks(numpy.arange(window_sizes.max(), window_sizes.min(), -1))
ax.set_xlabel('Positions')
ax.xaxis.set_ticks(numpy.arange(0, dataset.shape[0] - window_sizes.min() + 1))

pylab.imshow(linear_rms, cmap=cm.autumn, interpolation='none',
             extent=[0, dataset.shape[0] - window_sizes.min() + 1, window_sizes.max(),
                     window_sizes.min()], vmin=0, vmax=v_max)
pylab.title('Linear RMS')

ax = f.add_subplot(2, 3, 4)
ax.set_ylabel('Window Sizes')
ax.yaxis.set_ticks(numpy.arange(window_sizes.max(), window_sizes.min(), -1))
ax.set_xlabel('Positions')
ax.xaxis.set_ticks(numpy.arange(0, dataset.shape[0] - window_sizes.min() + 1))

im = pylab.imshow(exponential_rms, cmap=cm.autumn, interpolation='none',
                  extent=[0, dataset.shape[0] - window_sizes.min() + 1, window_sizes.max(), window_sizes.min()],
                  vmin=0, vmax=v_max)
pylab.title('Exponential RMS')

ax = f.add_subplot(2, 3, 5)
ax.set_ylabel('Window Sizes')
ax.yaxis.set_ticks(numpy.arange(window_sizes.max(), window_sizes.min(), -1))
ax.set_xlabel('Positions')
ax.xaxis.set_ticks(numpy.arange(0, dataset.shape[0] - window_sizes.min() + 1))

pylab.imshow(linear_rms - exponential_rms, cmap=cm.jet, interpolation='none',
             extent=[0, dataset.shape[0] - window_sizes.min() + 1, window_sizes.max(), window_sizes.min()])
pylab.title('Linear RMS - Exponential RMS')
pylab.colorbar()

f.subplots_adjust(right=0.8)
cbar_ax = f.add_axes([0.85, 0.15, 0.05, 0.7])
f.colorbar(im, cax=cbar_ax)
plt.show()
pylab.show()