%pylab inline

#!/usr/bin/env python
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt

# Model to fit the data with
def hill(x, baseline, amplitude, tconstant, hillcoef):
    """
    Return the hill equation curve for the x values.
    """
    return baseline+amplitude*(x**hillcoef)/(x**hillcoef+tconstant**hillcoef)


# Bootstrap resampling for fit function
class DataFit:
    
    def __init__(self, name, function):
        self.name = name
        data = np.loadtxt(self.name, delimiter=',')
        self.xdata = data[:, 0]
        self.ydata = data[:, 1]
        self.xfit = np.arange(self.xdata[0], self.xdata.max()+0.1, 0.1)

        # set the fit function/parameters
        self.set_fitfunc(function)
        self.compute_yfit()

    def set_fitfunc(self, namefunc):
        import inspect
        self.fitfunc = namefunc
        f = inspect.getargspec(self.fitfunc)
        n_params = len(f.args)-1 # minus the x argument
        self.p0 = [1]*n_params
        self.p0_args = f.args[1:]

    def compute_yfit(self):
        from scipy.optimize import curve_fit
        self.popt, self.pcov = curve_fit(self.fitfunc, self.xdata,
                self.ydata, p0=self.p0)
        params = []
        params.extend(self.popt)
        params.insert(0, self.xfit)
        yf = apply(self.fitfunc, params)
        self.yfit = yf

    def bootstrap_fitfunc(self, n, pos=1, alpha=0.05):
        """
        Bootstrap n times the fit of the data and store the n fits in a matrix,
        as well as the parameters of interest (position = pos) we want the
        1-alpha confidence interval
        """
        from scipy.optimize import curve_fit
        # calculate residual from fit
        params = list(self.popt)
        params.insert(0, self.xdata)
        ydatafit = apply(self.fitfunc, params)
        r = self.ydata - ydatafit
        # create matrix to store result of bootstrap
        self.matrix_bootstrap = np.zeros((self.xfit.size, n))
        boot_param = np.zeros(n)
        for i in xrange(n):
            new_popt = curve_fit(self.fitfunc, self.xdata,
                    ydatafit+r[sp.random.random_integers(0, len(r)-1,
                        len(r))], p0=self.popt)[0]
            p_newfit = list(new_popt)
            p_newfit.insert(0, self.xfit)
            self.matrix_bootstrap[:, i] = apply(self.fitfunc, p_newfit)
            boot_param[i] = new_popt[pos]
        # rank rows of matrix and extract confidence band limits
        ranked_matrix = self.matrix_bootstrap.copy()
        ranked_matrix[:].sort()
        self.ci_fitup = ranked_matrix[:, round(n*(1-alpha/2))-1]
        self.ci_fitlow = ranked_matrix[:, round(n*(alpha/2))-1]
        # rank parameter array and extract confidence interval
        self.bootstrap_param = boot_param
        boot_param.sort()
        self.ci_param = list([boot_param[round(n*(alpha/2)-1)],
                boot_param[n*(1-alpha/2)-1]])

# alpha to be fixed (% confidence interval) and number of resampling
alpha = 0.05
n = 10000

# Data file (2 columns, time and fluorescence value)
adult_hk2y_0glu = "all-data-adult/a-HK2-0glu/all_hk2yfp-0glu-adult.txt"
!cat "all-data-adult/a-HK2-0glu/all_hk2yfp-0glu-adult.txt"

# Fit of the data
d2y = DataFit(adult_hk2y_0glu, hill)
d2y.p0 = [100., -40., 10., 2.]
d2y.compute_yfit()

# Bootstrap of the residual of the fit
d2y.bootstrap_fitfunc(n, pos=2, alpha=alpha)

# tweak the plots
def format_axe(ax):
    ax.spines['right'].set_color('none')
    ax.spines['top'].set_color('none')
    ax.xaxis.set_ticks_position('bottom')
    ax.yaxis.set_ticks_position('left')
    ax.tick_params(axis='both', direction='out', width=3, length=7,
                   labelsize=16, pad=8)
    ax.spines['left'].set_linewidth(3)
    ax.spines['bottom'].set_linewidth(3)

fig = plt.figure(figsize=(12, 10))
ax1 = fig.add_subplot(221) # Data and original fit
ax2 = fig.add_subplot(222) # Individual bootstrapped fits (10000)
ax3 = fig.add_subplot(223) # Histogram of the time constants and 95% confidence intervals
ax4 = fig.add_subplot(224) # Original fit + gray shaded area for 95% confidence intervals of the fit

format_axe(ax1)
format_axe(ax2)
format_axe(ax3)
format_axe(ax4)

ax1.set_xlim(-1, 46)
ax1.set_ylim(20, 120)
ax2.set_xlim(-1, 46)
ax2.set_ylim(20, 120)
ax3.set_xlim(6, 14)
ax4.set_xlim(-1, 46)
ax4.set_ylim(20, 120)

ax1.set_xlabel('Time (min)', size=20)
ax2.set_xlabel('Time (min)', size=20)
ax3.set_xlabel('t$_{1/2}$ (min)', size=20)
ax4.set_xlabel('Time (min)', size=20)

ax1.set_ylabel('fluo (a.u.)', size=20)
ax2.set_ylabel('fluo (a.u.)', size=20)
ax3.set_ylabel('#', size=20)
ax4.set_ylabel('fluo (a.u.)', size=20)

ax1.plot(d2y.xdata, d2y.ydata, 'o', markersize=11, markerfacecolor='w',
         markeredgecolor='k', markeredgewidth=2)
ax1.plot(d2y.xfit, d2y.yfit, '-r', lw=3, alpha=1)

ax2.plot(d2y.xfit, d2y.matrix_bootstrap, '-', color='0.5', lw=1)
ax2.plot(d2y.xfit, d2y.yfit, '-y', lw=3)

ax3.hist(d2y.bootstrap_param, bins=100, facecolor='0.5', edgecolor='w')
ax3.axvline(x=d2y.popt[2], linewidth=3, color='r')
ax3.axvline(d2y.ci_param[0], ymin=0, ymax=0.3, linewidth=2, color='r')
ax3.axvline(d2y.ci_param[1], ymin=0, ymax=0.3, linewidth=2, color='r')

ax4.plot(d2y.xdata, d2y.ydata, 'o', markersize=11, markerfacecolor='w',
         markeredgecolor='k', markeredgewidth=2)
ax4.plot(d2y.xfit, d2y.yfit, '-r', lw=3, alpha=1)
ax4.fill_between(d2y.xfit, d2y.ci_fitlow, d2y.ci_fitup, color='k',
        alpha=0.5, edgecolor='none') 
ax4.plot(d2y.xfit, d2y.ci_fitup, '-r', lw=2, alpha=0.5)
ax4.plot(d2y.xfit, d2y.ci_fitlow, '-r', lw=2, alpha=0.5)

plt.tight_layout()