Impulse Response Function

In [4]:

    

#Import packages for basic calculations and plotting
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gammainc, gamma
%matplotlib inline

# First start with the input parameters

t = np.arange(0,3000.,1)
A = 1.
a = 1./10**2.8
n = 1.5

# Create the impulse response curve
X2 = A * ((a ** n) * (t ** (n-1)) * np.exp(-a*t)) / gamma(n)

# Create the step response curve
Fs = A * gammainc(n,t/a)

# Test of response function of Peterson et al. (2014)
A = 1.0
a = 10**2.8
X1 = A*t**(n-1)*np.exp(-t/a)
Fs = np.cumsum(A*t**(n-1)*np.exp(-t/a))

# Create the block response curve
Fb = Fs[1:] - Fs[0:-1]
Fb = np.append(0, Fb) #This is only done for drawing the graph as you normally
# Convolute 

plt.figure(figsize=(10,8))

plt.subplot(321)
plt.bar([0],[1],width=0.1, color='gray')
plt.xlim(-1,8)
plt.ylim(0,1.2)
plt.yticks([1],['$\infty$'])
plt.ylabel('Recharge [L/T]')
plt.subplot(322)
plt.plot(t,X1, 'k')
plt.plot(t,X2, 'r')
plt.ylabel('Response [-]')
plt.xlim(0,3000)


plt.subplot(323)
plt.bar([0],[1],width=8, color='gray')
plt.ylim(0,1.2)
plt.xlim(-1,8)
plt.xticks([0, 8],[0, '$\infty$'])
plt.ylabel('Recharge [L/T]')
plt.subplot(324)
plt.plot(np.arange(len(Fs)),Fs, 'k')
plt.ylabel('Response [-]')
plt.xlim(0,3000)


plt.subplot(325)
plt.bar([0],[1],width=1, color='gray')
plt.ylim(0,1.2)
plt.xlim(-1,8)
plt.ylabel('Recharge [L/T]')
plt.xlabel(' Time [T]')
plt.subplot(326)
plt.plot(np.arange(len(Fb)),Fb, 'k')
plt.ylabel('Response [-]')
plt.xlim(0,3000)
plt.xlabel(' Time [T]')


    

    




---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-a1efa79ba059> in <module>()
      1 #Import packages for basic calculations and plotting
      2 import numpy as np
----> 3 import matplotlib.pyplot as plt
      4 from scipy.special import gammainc, gamma
      5 get_ipython().magic(u'matplotlib inline')

/Applications/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py in <module>()
   1129 
   1130 # this is the instance used by the matplotlib classes
-> 1131 rcParams = rc_params()
   1132 
   1133 if rcParams['examples.directory']:

/Applications/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py in rc_params(fail_on_error)
    973         return ret
    974 
--> 975     return rc_params_from_file(fname, fail_on_error)
    976 
    977 

/Applications/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py in rc_params_from_file(fname, fail_on_error, use_default_template)
   1098         parameters specified in the file. (Useful for updating dicts.)
   1099     """
-> 1100     config_from_file = _rc_params_in_file(fname, fail_on_error)
   1101 
   1102     if not use_default_template:

/Applications/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py in _rc_params_in_file(fname, fail_on_error)
   1016     cnt = 0
   1017     rc_temp = {}
-> 1018     with _open_file_or_url(fname) as fd:
   1019         try:
   1020             for line in fd:

/Applications/anaconda/lib/python2.7/contextlib.pyc in __enter__(self)
     15     def __enter__(self):
     16         try:
---> 17             return self.gen.next()
     18         except StopIteration:
     19             raise RuntimeError("generator didn't yield")

/Applications/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py in _open_file_or_url(fname)
    998     else:
    999         fname = os.path.expanduser(fname)
-> 1000         encoding = locale.getdefaultlocale()[1]
   1001         if encoding is None:
   1002             encoding = "utf-8"

/Applications/anaconda/lib/python2.7/locale.pyc in getdefaultlocale(envvars)
    541     else:
    542         localename = 'C'
--> 543     return _parse_localename(localename)
    544 
    545 

/Applications/anaconda/lib/python2.7/locale.pyc in _parse_localename(localename)
    473     elif code == 'C':
    474         return None, None
--> 475     raise ValueError, 'unknown locale: %s' % localename
    476 
    477 def _build_localename(localetuple):

ValueError: unknown locale: UTF-8

In [3]:

    

# Define a function that calculates and plots the impulse response function
def impulseresponse_plot(A=1,a=0.01,n=1.5):
    t = np.linspace(0,5000)
    Fs = A * gammainc(n,t/a) #Step Response Function
    Fb = Fs[1:] - Fs[0:-1] # Block Response Function
    Fb = np.append(0,Fb) #This is done for the visualisation, not for actual modelling!!
    plt.plot(t,Fb,'r')
    plt.xlabel('Time [days]')
    plt.ylabel('')
    plt.legend('IRF', loc = 1)
    return Fb

interact(impulseresponse_plot, 
        A = widgets.FloatSliderWidget(min = 0.01, max = 400, step =1, description = 'A', value = 10),
        a = widgets.FloatSliderWidget(min = 1, max = 1000, step =0.005, value = 0.01),
        n = widgets.FloatSliderWidget(min = 0.01, max = 10, step =0.10, value = 1.5))


    

    




---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-3-3a77bc3f755a> in <module>()
     12 
     13 interact(impulseresponse_plot, 
---> 14         A = widgets.FloatSliderWidget(min = 0.01, max = 400, step =1, description = 'A', value = 10),
     15         a = widgets.FloatSliderWidget(min = 1, max = 1000, step =0.005, value = 0.01),
     16         n = widgets.FloatSliderWidget(min = 0.01, max = 10, step =0.10, value = 1.5))

/Applications/anaconda/lib/python2.7/site-packages/IPython/utils/shimmodule.pyc in __getattr__(self, key)
     90             return import_item(name)
     91         except ImportError:
---> 92             raise AttributeError(key)

AttributeError: FloatSliderWidget

In [21]:

    

t = np.arange(0,3000.,1)
A = 0.3
a = 10**2.8 # 1/a to change for Peterson
n = 1.5
b = 100.

def f(A,a,n,t):
    return A*(t-b)**(n-1)*np.exp(-(t-b)/a) #change /a to * a to comply with Peterson

t_peak = (n-1)*a #time of peak
print t_peak
value_tp = f(A,a,n,t_peak)
print value_tp

def f(A,a,n,t):
    Fs = A * gammainc(n, t/a) # Step response function based on pearsonIII
    Fb = np.append(0,Fs[1:] - Fs[0:-1]) #block reponse function
    #return Fb
    #return (A*t**(n-1)*np.exp(-t/a))/(((n-1)*a)**(n-1)*np.exp(1-n))
    return A*(t-b)**(n-1)*np.exp(-(t-b)/a)

plt.plot(f(A,a,n,t))


    

    




315.47867224
3.12973232862






    
Out[21]:





[<matplotlib.lines.Line2D at 0x10a901790>]






    

In [5]:

    

#Create some random shit
from lmfit import minimize, Parameters, Parameter, report_fit
from scipy import random
x = random.rand(len(t))*0.1

f_test = f(A,a,n,t) + x

plt.plot(f_test)


    

    
Out[5]:





[<matplotlib.lines.Line2D at 0x1095f5310>]






    

In [6]:

    

# Define an objective function
def residuals(X0):
    A = X0['A'].value
    a = X0['a'].value
    n = X0['n'].value
    return (f(A,a,n,t)-f_test)

X0 = Parameters()
#           (Name,  Value,  Vary,   Min,  Max,  Expr)
X0.add_many(('A',   2.0,    True,   None, None,  None),
           ('a',    10**2.8,    True,   None, None,  None),
           ('n',    1.5,    True,   None, None,  None),
)

result = minimize(residuals, X0, method='leastsq')
report_fit(result)

#print residuals(X0)
param = result.params
print param


    

    




[[Fit Statistics]]
    # function evals   = 23
    # data points      = 3000
    # variables        = 3
    chi-square         = 2.946
    reduced chi-square = 0.001
[[Variables]]
    A:   1.03716516 +/- 0.001300 (0.13%) (init= 2)
    a:   720.078790 +/- 2.261413 (0.31%) (init= 630.9573)
    n:   1.43336742 +/- 0.002498 (0.17%) (init= 1.5)
[[Correlations]] (unreported correlations are <  0.100)
    C(a, n)                      = -0.879 
    C(A, a)                      = -0.408 
    C(A, n)                      =  0.144 
Parameters([('A', <Parameter 'A', value=1.0371651694278736 +/- 0.0013, bounds=[-inf:inf]>), ('a', <Parameter 'a', value=720.07879040758053 +/- 2.26, bounds=[-inf:inf]>), ('n', <Parameter 'n', value=1.4333674224554911 +/- 0.0025, bounds=[-inf:inf]>)])

In [7]:

    

x1,x2,x3 = param['A'].value, param['a'].value, param['n'].value
plt.plot(f_test)
plt.plot(f(x1,x2,x3,t), 'r')


    

    
Out[7]:





[<matplotlib.lines.Line2D at 0x10959eb50>]






    

In [8]:

    

result.covar


    

    
Out[8]:





array([[  1.69144846e-06,  -1.20002270e-03,   4.67774435e-07],
       [ -1.20002270e-03,   5.11398882e+00,  -4.96446578e-03],
       [  4.67774435e-07,  -4.96446578e-03,   6.24122484e-06]])

In [ ]:

    

result.covar


    

In [ ]: