#!/usr/bin/env python
"""
Module defines common functions that are used in curve_fit or fmin
parameter estimations.
For all fit functions, it defines the functions in two forms (ex. of 3
params):
`func(x, p1, p2, p3)`
`func_p(x, p)` with `p[0:3]`
The first form can be used, for example, with
`scipy.optimize.curve_fit` (ex. function f1x=a+b/x):
`p, cov = scipy.optimize.curve_fit(functions.f1x, x, y, p0=[p0,p1])`
It also defines two cost functions along with the fit functions, one
with the absolute sum, one with the squared sum of the deviations:
`cost_func` `sum(abs(obs-func))`
`cost2_func` `sum((obs-func)**2)`
These cost functions can be used, for example, with
`scipy.optimize.minimize`:
`p = scipy.optimize.minimize(jams.functions.cost_f1x, np.array([p1,p2]), args=(x,y), method='Nelder-Mead', options={'disp':False})`
Note the different argument orders:
`curvefit` needs `f(x,*args)` with the independent variable as the
first argument and the parameters to fit as separate remaining
arguments.
`minimize` is a general minimiser with respect to the first argument,
i.e. `func(p,*args)`.
The module provides also two common cost functions (absolute and
squared deviations) where any function in the form `func(x, p)` can be
used as second argument:
`cost_abs(p, func, x, y)`
`cost_square(p, func, x, y)`
This means, for example `cost_f1x(p, x, y)` is the same as
`cost_abs(p, functions.f1x_p, x, y)`.
For example:
`p = scipy.optimize.minimize(jams.functions.cost_abs, np.array([p1,p2]), args=(functions.f1x_p,x,y), method='Nelder-Mead', options={'disp':False})`
The current functions are (the functions have the name in the first
column. The seond form has a '\_p' appended to the name. The cost
functions, which have 'cost\_' and 'cost2\_' prepended to the name.):
arrhenius 1 param: Arrhenius temperature dependence of biochemical rates: `exp((T-TC25)*E/(T25*R*(T+T0)))`, parameter: E
f1x 2 params: General 1/x function: `a + b/x`
fexp 3 params: General exponential function: `a + b * exp(c*x)`
gauss 2 params: Gauss function: `1/(sig*sqrt(2*pi)) *exp(-(x-mu)**2/(2*sig**2))`, parameter: mu, sig
lasslop 6 params: Lasslop et al. (2010) a rectangular, hyperbolic light-response GPP with Lloyd & Taylor (1994) respiration and the maximum canopy uptake rate at light saturation decreases exponentially with VPD as in Koerner (1995)
line0 1 params: Straight line: `a*x`
line 2 params: Straight line: `a + b*x`
lloyd_fix 2 params: Lloyd & Taylor (1994) Arrhenius type with T0=-46.02 degC and Tref=10 degC
lloyd_only_rref 1 param: Lloyd & Taylor (1994) Arrhenius type with fixed exponential term
logistic 3 params: Logistic function: `a/(1+exp(-b(x-c)))`
logistic_offset 4 params: Logistic function with offset: `a/(1+exp(-b(x-c))) + d`
logistic2_offset 7 params: Double logistic function with offset `L1/(1+exp(-k1(x-x01))) - L2/(1+exp(-k2(x-x02))) + a`
poly n params: General polynomial: `c0 + c1*x + c2*x**2 + ... + cn*x**n`
sabx 2 params: sqrt(f1x), i.e. general sqrt(1/x) function: `sqrt(a + b/x)`
see 3 params: Sequential Elementary Effects fitting function: `a*(x-b)**c`
This module was written by Matthias Cuntz while at Department of
Computational Hydrosystems, Helmholtz Centre for Environmental
Research - UFZ, Leipzig, Germany, and continued while at Institut
National de Recherche pour l'Agriculture, l'Alimentation et
l'Environnement (INRAE), Nancy, France.
Copyright (c) 2012-2020 Matthias Cuntz - mc (at) macu (dot) de
Released under the MIT License; see LICENSE file for details.
* Written Dec 2012 by Matthias Cuntz (mc (at) macu (dot) de)
* Ported to Python 3, Feb 2013, Matthias Cuntz
* Added general cost functions cost_abs and cost_square, May 2013, Matthias Cuntz
* Added line0, Feb 2014, Matthias Cuntz
* Removed multiline_p, May 2020, Matthias Cuntz
* Changed to Sphinx docstring and numpydoc, May 2020, Matthias Cuntz
.. moduleauthor:: Matthias Cuntz
The following functions are provided:
.. autosummary::
cost_abs
cost_square
arrhenius
arrhenius_p
cost_arrhenius
cost2_arrhenius
f1x
f1x_p
cost_f1x
cost2_f1x
fexp
fexp_p
cost_fexp
cost2_fexp
gauss
gauss_p
cost_gauss
cost2_gauss
lasslop
lasslop_p
cost_lasslop
cost2_lasslop
line
line_p
cost_line
cost2_line
line0
line0_p
cost_line0
cost2_line0
lloyd_fix
lloyd_fix_p
cost_lloyd_fix
cost2_lloyd_fix
lloyd_only_rref
lloyd_only_rref_p
cost_lloyd_only_rref
cost2_lloyd_only_rref
sabx
sabx_p
cost_sabx
cost2_sabx
poly
poly_p
cost_poly
cost2_poly
cost_logistic
cost2_logistic
cost_logistic_offset
cost2_logistic_offset
cost_logistic2_offset
cost2_logistic2_offset
see
see_p
cost_see
cost2_see
"""
from __future__ import division, absolute_import, print_function
import numpy as np
import scipy.special as sp
try: # import package
from .logistic_function import logistic_p, logistic_offset_p, logistic2_offset_p
from ..const import T0, T25, R
except:
try: # e.g. python nee2gpp.py
from functions.logistic_function import logistic_p, logistic_offset_p, logistic2_offset_p
from const import T0, T25, R
except: # python fit_functions.py
from logistic_function import logistic_p, logistic_offset_p, logistic2_offset_p
T0 = 273.15 # Celcius <-> Kelvin [K]
T25 = 298.15 # Standard ambient temperature [K]
R = 8.3144621 # Ideal gas constant [J K^-1 mol^-1]
__all__ = ['cost_abs', 'cost_square',
'arrhenius', 'arrhenius_p', 'cost_arrhenius', 'cost2_arrhenius',
'f1x', 'f1x_p', 'cost_f1x', 'cost2_f1x',
'fexp', 'fexp_p', 'cost_fexp', 'cost2_fexp',
'gauss', 'gauss_p', 'cost_gauss', 'cost2_gauss',
'lasslop', 'lasslop_p', 'cost_lasslop', 'cost2_lasslop',
'line', 'line_p', 'cost_line', 'cost2_line',
'line0', 'line0_p', 'cost_line0', 'cost2_line0',
'lloyd_fix', 'lloyd_fix_p', 'cost_lloyd_fix', 'cost2_lloyd_fix',
'lloyd_only_rref', 'lloyd_only_rref_p', 'cost_lloyd_only_rref', 'cost2_lloyd_only_rref',
'sabx', 'sabx_p', 'cost_sabx', 'cost2_sabx',
'poly', 'poly_p', 'cost_poly', 'cost2_poly',
'cost_logistic', 'cost2_logistic',
'cost_logistic_offset', 'cost2_logistic_offset',
'cost_logistic2_offset', 'cost2_logistic2_offset',
'see', 'see_p', 'cost_see', 'cost2_see']
# -----------------------------------------------------------
# general cost functions
[docs]def cost_abs(p, func, x, y):
"""
General cost function for robust optimising `func(x,p)` vs. `y` with sum of absolute deviations.
Parameters
----------
p : iterable of floats
parameters
func : callable
`fun(x,p) -> float`
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(y-func(x,p)))
[docs]def cost_square(p, func, x, y):
"""
General cost function for optimising `func(x,p)` vs. `y` with sum of square deviations.
Parameters
----------
p : iterable of floats
parameters
func : callable
`fun(x,p) -> float`
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((y-func(x,p))**2)
# -----------------------------------------------------------
# arrhenius
[docs]def arrhenius(T, E):
"""
Arrhenius temperature dependence of rates.
Parameters
----------
T : float or array_like of floats
temperature [degC]
E : float
activation energy [J]
Returns
-------
float
function value(s)
"""
return np.exp((T-(T25-T0))*E/(T25*R*(T+T0)))
[docs]def arrhenius_p(T, p):
"""
Arrhenius temperature dependence of rates.
Parameters
----------
T : float or array_like of floats
temperature [degC]
p : iterable
`p[0]` is activation energy [J]
Returns
-------
float
function value(s)
"""
return np.exp((T-(T25-T0))*p[0]/(T25*R*(T+T0)))
[docs]def cost_arrhenius(p, T, rate):
"""
Sum of absolute deviations of obs and arrhenius function.
Parameters
----------
p : iterable of floats
`p[0]` is activation energy [J]
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(rate-arrhenius_p(T,p)))
[docs]def cost2_arrhenius(p, T, rate):
"""
Sum of squared deviations of obs and arrhenius.
Parameters
----------
p : iterable of floats
`p[0]` is activation energy [J]
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((rate-arrhenius_p(T,p))**2)
# -----------------------------------------------------------
# a+b/x
[docs]def f1x(x,a,b):
"""
General 1/x function: a + b/x
Parameters
----------
x : float or array_like of floats
independent variable
a : float
first parameter
b : float
second parameter
Returns
-------
float
function value(s)
"""
return a+b/x
[docs]def f1x_p(x,p):
"""
General 1/x function: a + b/x
Parameters
----------
x : float or array_like of floats
independent variable
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` a
`p[1]` b
Returns
-------
float
function value(s)
"""
return p[0]+p[1]/x
[docs]def cost_f1x(p,x,y):
"""
Sum of absolute deviations of obs and general 1/x function: a + b/x
Parameters
----------
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` a
`p[1]` b
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(y-f1x_p(x,p)))
[docs]def cost2_f1x(p,x,y):
"""
Sum of squared deviations of obs and general 1/x function: a + b/x
Parameters
----------
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` a
`p[1]` b
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((y-f1x_p(x,p))**2)
# -----------------------------------------------------------
# a+b*exp(c*x)
[docs]def fexp(x,a,b,c):
"""
General exponential function: a + b * exp(c*x)
Parameters
----------
x : float or array_like of floats
independent variable
a : float
first parameter
b : float
second parameter
c : float
third parameter
Returns
-------
float
function value(s)
"""
return a+b*np.exp(c*x)
[docs]def fexp_p(x,p):
"""
General exponential function: a + b * exp(c*x)
Parameters
----------
x : float or array_like of floats
independent variable
p : iterable of floats
parameters (`len(p)=3`)
`p[0]` a
`p[1]` b
`p[2]` c
Returns
-------
float
function value(s)
"""
return p[0]+p[1]*np.exp(p[2]*x)
[docs]def cost_fexp(p,x,y):
"""
Sum of absolute deviations of obs and general exponential function: a + b * exp(c*x)
Parameters
----------
p : iterable of floats
parameters (`len(p)=3`)
`p[0]` a
`p[1]` b
`p[2]` c
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(y-fexp_p(x,p)))
[docs]def cost2_fexp(p,x,y):
"""
Sum of squared deviations of obs and general exponential function: a + b * exp(c*x)
Parameters
----------
p : iterable of floats
parameters (`len(p)=3`)
`p[0]` a
`p[1]` b
`p[2]` c
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((y-fexp_p(x,p))**2)
# -----------------------------------------------------------
# Gauss: 1/(sig*sqrt(2*pi)) *exp(-(x-mu)**2/(2*sig**2))
[docs]def gauss(x,mu,sig):
"""
Gauss function: 1 / (sqrt(2*pi)*sig) * exp( -(x-mu)**2 / (2*sig**2) )
Parameters
----------
x : float or array_like of floats
independent variable
mu : float
mean
sig : float
width
Returns
-------
float
function value(s)
"""
return np.exp(-(x-mu)**2/(2.*sig**2))/(sig*np.sqrt(2.*np.pi))
[docs]def gauss_p(x,p):
"""
Gauss function: 1 / (sqrt(2*pi)*sig) * exp( -(x-mu)**2 / (2*sig**2) )
Parameters
----------
x : float or array_like of floats
independent variable
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` mean
`p[1]` width
Returns
-------
float
function value(s)
"""
return np.exp(-(x-p[0])**2/(2.*p[1]**2))/(p[1]*np.sqrt(2.*np.pi))
[docs]def cost_gauss(p,x,y):
"""
Sum of absolute deviations of obs and Gauss function: 1 / (sqrt(2*pi)*sig) * exp( -(x-mu)**2 / (2*sig**2) )
Parameters
----------
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` mean
`p[1]` width
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(y-gauss_p(x,p)))
[docs]def cost2_gauss(p,x,y):
"""
Sum of squared deviations of obs and Gauss function: 1 / (sqrt(2*pi)*sig) * exp( -(x-mu)**2 / (2*sig**2) )
Parameters
----------
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` mean
`p[1]` width
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((y-gauss_p(x,p))**2)
# -----------------------------------------------------------
# lasslop
[docs]def lasslop(Rg, et, VPD, alpha, beta0, k, Rref):
"""
Lasslop et al. (2010) is the rectangular, hyperbolic light-response
of NEE as by Falge et al. (2001), where the respiration is calculated
with Lloyd & Taylor (1994), and the maximum canopy uptake rate at
light saturation decreases exponentially with VPD as in Koerner (1995).
Parameters
----------
Rg : float or array_like of floats
Global radiation [W m-2]
et : float or array_like of floats
Exponential in Lloyd & Taylor: np.exp(E0*(1./(Tref-T0)-1./(T-T0))) []
VPD : float or array_like of floats
Vapour Pressure Deficit [Pa]
alpha : float
Light use efficiency, i.e. initial slope of light response curve [umol(C) J-1]
beta0 : float
Maximum CO2 uptake rate at VPD0=10 hPa [umol(C) m-2 s-1]
k : float
e-folding of exponential decrease of maximum CO2 uptake with VPD increase [Pa-1]
Rref : float
Respiration at Tref (10 degC) [umol(C) m-2 s-1]
Returns
-------
float
net ecosystem exchange [umol(CO2) m-2 s-1]
"""
# Lloyd & Taylor (1994)
gamma = Rref*et
# Koerner (1995)
VPD0 = 1000. # 10 hPa
kk = np.clip(-k*(VPD-VPD0), -600., 600.)
beta = np.where(VPD > VPD0, beta0*np.exp(kk), beta0)
return -alpha*beta*Rg/(alpha*Rg+beta) + gamma
[docs]def lasslop_p(Rg, et, VPD, p):
"""
Lasslop et al. (2010) is the rectangular, hyperbolic light-response
of NEE as by Falge et al. (2001), where the respiration is calculated
with Lloyd & Taylor (1994), and the maximum canopy uptake rate at
light saturation decreases exponentially with VPD as in Koerner (1995).
Parameters
----------
Rg : float or array_like of floats
Global radiation [W m-2]
et : float or array_like of floats
Exponential in Lloyd & Taylor: np.exp(E0*(1./(Tref-T0)-1./(T-T0))) []
VPD : float or array_like of floats
Vapour Pressure Deficit [Pa]
p : iterable of floats
parameters (`len(p)=4`)
`p[0]` Light use efficiency, i.e. initial slope of light response curve [umol(C) J-1]
`p[1]` Maximum CO2 uptake rate at VPD0=10 hPa [umol(C) m-2 s-1]
`p[2]` e-folding of exponential decrease of maximum CO2 uptake with VPD increase [Pa-1]
`p[3]` Respiration at Tref (10 degC) [umol(C) m-2 s-1]
Returns
-------
float
net ecosystem exchange [umol(CO2) m-2 s-1]
"""
# Lloyd & Taylor (1994)
gamma = p[3]*et
# Koerner (1995)
VPD0 = 1000. # 10 hPa
kk = np.clip(-p[2]*(VPD-VPD0), -600., 600.)
beta = np.where(VPD > VPD0, p[1]*np.exp(kk), p[1])
return -p[0]*beta*Rg/(p[0]*Rg+beta) + gamma
[docs]def cost_lasslop(p, Rg, et, VPD, NEE):
"""
Sum of absolute deviations of obs and Lasslop.
Parameters
----------
p : iterable of floats
parameters (`len(p)=4`)
`p[0]` Light use efficiency, i.e. initial slope of light response curve [umol(C) J-1]
`p[1]` Maximum CO2 uptake rate at VPD0=10 hPa [umol(C) m-2 s-1]
`p[2]` e-folding of exponential decrease of maximum CO2 uptake with VPD increase [Pa-1]
`p[3]` Respiration at Tref (10 degC) [umol(C) m-2 s-1]
Rg : float or array_like of floats
Global radiation [W m-2]
et : float or array_like of floats
Exponential in Lloyd & Taylor: np.exp(E0*(1./(Tref-T0)-1./(T-T0))) []
VPD : float or array_like of floats
Vapour Pressure Deficit [Pa]
NEE : float or array_like of floats
Observed net ecosystem exchange [umol(CO2) m-2 s-1]
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(NEE-lasslop(Rg, et, VPD, p[0], p[1], p[2], p[3])))
[docs]def cost2_lasslop(p, Rg, et, VPD, NEE):
"""
Sum of squared deviations of obs and Lasslop.
Parameters
----------
p : iterable of floats
parameters (`len(p)=4`)
`p[0]` Light use efficiency, i.e. initial slope of light response curve [umol(C) J-1]
`p[1]` Maximum CO2 uptake rate at VPD0=10 hPa [umol(C) m-2 s-1]
`p[2]` e-folding of exponential decrease of maximum CO2 uptake with VPD increase [Pa-1]
`p[3]` Respiration at Tref (10 degC) [umol(C) m-2 s-1]
Rg : float or array_like of floats
Global radiation [W m-2]
et : float or array_like of floats
Exponential in Lloyd & Taylor: np.exp(E0*(1./(Tref-T0)-1./(T-T0))) []
VPD : float or array_like of floats
Vapour Pressure Deficit [Pa]
NEE : float or array_like of floats
Observed net ecosystem exchange [umol(CO2) m-2 s-1]
Returns
-------
float
sum of squared deviations
"""
return np.sum((NEE-lasslop(Rg, et, VPD, p[0], p[1], p[2], p[3]))**2)
# -----------------------------------------------------------
# a+b*x
[docs]def line(x,a,b):
"""
Straight line: a + b*x
Parameters
----------
x : float or array_like of floats
independent variable
a : float
first parameter
b : float
second parameter
Returns
-------
float
function value(s)
"""
return a+b*x
[docs]def line_p(x,p):
"""
Straight line: a + b*x
Parameters
----------
x : float or array_like of floats
independent variable
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` a
`p[1]` b
Returns
-------
float
function value(s)
"""
return p[0]+p[1]*x
[docs]def cost_line(p,x,y):
"""
Sum of absolute deviations of obs and straight line: a + b*x
Parameters
----------
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` a
`p[1]` b
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(y-line_p(x,p)))
[docs]def cost2_line(p,x,y):
"""
Sum of squared deviations of obs and straight line: a + b*x
Parameters
----------
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` a
`p[1]` b
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((y-line_p(x,p))**2)
# -----------------------------------------------------------
# b*x
[docs]def line0(x,a):
"""
Straight line through origin: a*x
Parameters
----------
x : float or array_like of floats
independent variable
a : float
first parameter
Returns
-------
float
function value(s)
"""
return a*x
[docs]def line0_p(x,p):
"""
Straight line through origin: a*x
Parameters
----------
x : float or array_like of floats
independent variable
p : iterable of floats
`p[0]` is a
Returns
-------
float
function value(s)
"""
return p*x
[docs]def cost_line0(p,x,y):
"""
Sum of absolute deviations of obs and straight line through origin: a*x
Parameters
----------
p : iterable of floats
`p[0]` is a
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(y-line0_p(x,p)))
[docs]def cost2_line0(p,x,y):
"""
Sum of squared deviations of obs and straight line through origin: a*x
Parameters
----------
p : iterable of floats
`p[0]` is a
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((y-line0_p(x,p))**2)
# -----------------------------------------------------------
# lloyd_fix
[docs]def lloyd_fix(T, Rref, E0):
"""
Lloyd & Taylor (1994) Arrhenius type with T0=-46.02 degC and Tref=10 degC
Parameters
----------
T : float or array_like of floats
Temperature [K]
Rref : float
Respiration at Tref=10 degC [umol(C) m-2 s-1]
E0 : float
Activation energy [K]
Returns
-------
float
Respiration [umol(C) m-2 s-1]
"""
Tref = 283.15 # 10 [degC]
T0 = 227.13 # -46.02 [degC]
return Rref*np.exp(E0*(1./(Tref-T0)-1./(T-T0)))
[docs]def lloyd_fix_p(T, p):
"""
Lloyd & Taylor (1994) Arrhenius type with T0=-46.02 degC and Tref=10 degC
Parameters
----------
T : float or array_like of floats
Temperature [K]
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` Respiration at Tref=10 degC [umol(C) m-2 s-1]
`p[1]` Activation energy [K]
Returns
-------
float
Respiration [umol(C) m-2 s-1]
"""
Tref = 283.15 # 10 [degC]
T0 = 227.13 # -46.02 [degC]
return p[0]*np.exp(p[1]*(1./(Tref-T0)-1./(T-T0)))
[docs]def cost_lloyd_fix(p, T, resp):
"""
Sum of absolute deviations of obs and Lloyd & Taylor (1994) Arrhenius type.
Parameters
----------
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` Respiration at Tref=10 degC [umol(C) m-2 s-1]
`p[1]` Activation energy [K]
T : float or array_like of floats
Temperature [K]
resp : float or array_like of floats
Observed respiration [umol(C) m-2 s-1]
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(resp-lloyd_fix_p(T,p)))
[docs]def cost2_lloyd_fix(p, T, resp):
"""
Sum of squared deviations of obs and Lloyd & Taylor (1994) Arrhenius type.
Parameters
----------
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` Respiration at Tref=10 degC [umol(C) m-2 s-1]
`p[1]` Activation energy [K]
T : float or array_like of floats
Temperature [K]
resp : float or array_like of floats
Observed respiration [umol(C) m-2 s-1]
Returns
-------
float
sum of squared deviations
"""
return np.sum((resp-lloyd_fix_p(T,p))**2)
# -----------------------------------------------------------
# lloyd_only_rref
[docs]def lloyd_only_rref(et, Rref):
"""
If E0 is know in Lloyd & Taylor (1994) then one can calc
the exponential term outside the routine and the fitting
becomes linear. One could also use functions.line0.
Parameters
----------
et : float or array_like of floats
exp-term in Lloyd & Taylor
Rref : float
Respiration at Tref=10 degC [umol(C) m-2 s-1]
Returns
-------
float
Respiration [umol(C) m-2 s-1]
"""
return Rref*et
[docs]def lloyd_only_rref_p(et, p):
"""
If E0 is know in Lloyd & Taylor (1994) then one can calc
the exponential term outside the routine and the fitting
becomes linear. One could also use functions.line0.
Parameters
----------
et : float or array_like of floats
exp-term in Lloyd & Taylor
p : iterable of floats
`p[0]` is respiration at Tref=10 degC [umol(C) m-2 s-1]
Returns
-------
float
Respiration [umol(C) m-2 s-1]
"""
return p[0]*et
[docs]def cost_lloyd_only_rref(p, et, resp):
"""
Sum of absolute deviations of obs and Lloyd & Taylor with known exponential term.
Parameters
----------
p : iterable of floats
`p[0]` is respiration at Tref=10 degC [umol(C) m-2 s-1]
et : float or array_like of floats
exp-term in Lloyd & Taylor
resp : float or array_like of floats
Observed respiration [umol(C) m-2 s-1]
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(resp-lloyd_only_rref_p(et,p)))
[docs]def cost2_lloyd_only_rref(p, et, resp):
"""
Sum of squared deviations of obs and Lloyd & Taylor with known exponential term.
Parameters
----------
p : iterable of floats
`p[0]` is respiration at Tref=10 degC [umol(C) m-2 s-1]
et : float or array_like of floats
exp-term in Lloyd & Taylor
resp : float or array_like of floats
Observed respiration [umol(C) m-2 s-1]
Returns
-------
float
sum of squared deviations
"""
return np.sum((resp-lloyd_only_rref_p(et,p))**2)
# -----------------------------------------------------------
# sqrt(a + b/x) - theoretical form of Jackknife-after-bootstrap
[docs]def sabx(x, a, b):
"""
Square root of general 1/x function: sqrt(a + b/x)
Parameters
----------
x : float or array_like of floats
independent variable
a : float
first parameter
b : float
second parameter
Returns
-------
float
function value(s)
"""
return np.sqrt(a+b/x)
[docs]def sabx_p(x, p):
"""
Square root of general 1/x function: sqrt(a + b/x)
Parameters
----------
x : float or array_like of floats
independent variable
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` a
`p[1]` b
Returns
-------
float
function value(s)
"""
return np.sqrt(p[0]+p[1]/x)
[docs]def cost_sabx(p,x,y):
"""
Sum of absolute deviations of obs and square root of general 1/x function: sqrt(a + b/x)
Parameters
----------
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` a
`p[1]` b
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(y-sabx_p(x,p)))
[docs]def cost2_sabx(p,x,y):
"""
Sum of squared deviations of obs and square root of general 1/x function: sqrt(a + b/x)
Parameters
----------
p : iterable of floats
parameters (`len(p)=2`)
`p[0]` a
`p[1]` b
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((y-sabx_p(x,p))**2)
# -----------------------------------------------------------
# c0 + c1*x + c2*x**2 + ... + cn*x**n
[docs]def poly(x,*args):
"""
General polynomial: c0 + c1*x + c2*x**2 + ... + cn*x**n
Parameters
----------
x : float or array_like of floats
independent variable
*args : float
parameters `len(args)=n+1`
Returns
-------
float
function value(s)
"""
return np.polynomial.polynomial.polyval(x, list(args))
[docs]def poly_p(x,p):
"""
General polynomial: c0 + c1*x + c2*x**2 + ... + cn*x**n
Parameters
----------
x : float or array_like of floats
independent variable
p : iterable of floats
parameters (`len(p)=n+1`)
Returns
-------
float
function value(s)
"""
return np.polynomial.polynomial.polyval(x, p)
[docs]def cost_poly(p,x,y):
"""
Sum of absolute deviations of obs and general polynomial: c0 + c1*x + c2*x**2 + ... + cn*x**n
Parameters
----------
p : iterable of floats
parameters (`len(p)=n+1`)
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(y-poly_p(x,p)))
[docs]def cost2_poly(p,x,y):
"""
Sum of squared deviations of obs and general polynomial: c0 + c1*x + c2*x**2 + ... + cn*x**n
Parameters
----------
p : iterable of floats
parameters (`len(p)=n+1`)
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((y-poly_p(x,p))**2)
# -----------------------------------------------------------
# a/(1+exp(-b(x-c))) - logistic function
[docs]def cost_logistic(p, x, y):
"""
Sum of absolute deviations of obs and logistic function L/(1+exp(-k(x-x0)))
Parameters
----------
p : iterable of floats
parameters (`len(p)=3`)
`p[0]` L - Maximum of logistic function
`p[1]` k - Steepness of logistic function
`p[2]` x0 - Inflection point of logistic function
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(y-logistic_p(x,p)))
[docs]def cost2_logistic(p,x,y):
"""
Sum of squared deviations of obs and logistic function L/(1+exp(-k(x-x0)))
Parameters
----------
p : iterable of floats
parameters (`len(p)=3`)
`p[0]` L - Maximum of logistic function
`p[1]` k - Steepness of logistic function
`p[2]` x0 - Inflection point of logistic function
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((y-logistic_p(x,p))**2)
# -----------------------------------------------------------
# a/(1+exp(-b(x-c))) + d - logistic function with offset
[docs]def cost_logistic_offset(p, x, y):
"""
Sum of absolute deviations of obs and logistic function 1/x function: L/(1+exp(-k(x-x0))) + a
Parameters
----------
p : iterable of floats
parameters (`len(p)=4`)
`p[0]` L - Maximum of logistic function
`p[1]` k - Steepness of logistic function
`p[2]` x0 - Inflection point of logistic function
`p[3]` a - Offset of logistic function
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(y-logistic_offset_p(x,p)))
[docs]def cost2_logistic_offset(p,x,y):
"""
Sum of squared deviations of obs and logistic function 1/x function: L/(1+exp(-k(x-x0))) + a
Parameters
----------
p : iterable of floats
parameters (`len(p)=4`)
`p[0]` L - Maximum of logistic function
`p[1]` k - Steepness of logistic function
`p[2]` x0 - Inflection point of logistic function
`p[3]` a - Offset of logistic function
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((y-logistic_offset_p(x,p))**2)
# -----------------------------------------------------------
# L1/(1+exp(-k1(x-x01))) - L2/(1+exp(-k2(x-x02))) + a2 - double logistic function with offset
[docs]def cost_logistic2_offset(p, x, y):
"""
Sum of absolute deviations of obs and double logistic function with offset:
L1/(1+exp(-k1(x-x01))) - L2/(1+exp(-k2(x-x02))) + a
Parameters
----------
p : iterable of floats
parameters (`len(p)=7`)
`p[0]` L1 - Maximum of first logistic function
`p[1]` k1 - Steepness of first logistic function
`p[2]` x01 - Inflection point of first logistic function
`p[3]` L2 - Maximum of second logistic function
`p[4]` k2 - Steepness of second logistic function
`p[5]` x02 - Inflection point of second logistic function
`p[6]` a - Offset of double logistic function
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(y-logistic2_offset_p(x,p)))
[docs]def cost2_logistic2_offset(p,x,y):
"""
Sum of squared deviations of obs and double logistic function with offset:
L1/(1+exp(-k1(x-x01))) - L2/(1+exp(-k2(x-x02))) + a
Parameters
----------
p : iterable of floats
parameters (`len(p)=7`)
`p[0]` L1 - Maximum of first logistic function
`p[1]` k1 - Steepness of first logistic function
`p[2]` x01 - Inflection point of first logistic function
`p[3]` L2 - Maximum of second logistic function
`p[4]` k2 - Steepness of second logistic function
`p[5]` x02 - Inflection point of second logistic function
`p[6]` a - Offset of double logistic function
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((y-logistic2_offset_p(x,p))**2)
# -----------------------------------------------------------
# a*(x-b)**c - Sequential Elementary Effects fitting function
[docs]def see(x, a, b, c):
"""
Fit function of Sequential Elementary Effects: a * (x-b)**c
Parameters
----------
x : float or array_like of floats
independent variable
a : float
first parameter
b : float
second parameter
c : float
third parameter
Returns
-------
float
function value(s)
"""
return np.where((x-b)<0., 0., a*(x-b)**c)
[docs]def see_p(x, p):
"""
Fit function of Sequential Elementary Effects: a * (x-b)**c
Parameters
----------
x : float or array_like of floats
independent variable
p : iterable of floats
parameters (`len(p)=3`)
`p[0]` a
`p[1]` b
`p[2]` c
Returns
-------
float
function value(s)
"""
return np.where((x-p[1]) < 0., 0., p[0] * (x-p[1])**p[2])
[docs]def cost_see(p, x, y):
"""
Sum of absolute deviations of obs and fit function of Sequential Elementary Effects: a * (x-b)**c
Parameters
----------
p : iterable of floats
parameters (`len(p)=3`)
`p[0]` a
`p[1]` b
`p[2]` c
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of absolute deviations
"""
return np.sum(np.abs(y-see_p(x,p)))
[docs]def cost2_see(p,x,y):
"""
Sum of squared deviations of obs and fit function of Sequential Elementary Effects: a * (x-b)**c
Parameters
----------
p : iterable of floats
parameters (`len(p)=3`)
`p[0]` a
`p[1]` b
`p[2]` c
x : float or array_like of floats
independent variable
y : float or array_like of floats
dependent variable, observations
Returns
-------
float
sum of squared deviations
"""
return np.sum((y-see_p(x,p))**2)
# -----------------------------------------------------------
if __name__ == '__main__':
import doctest
doctest.testmod(optionflags=doctest.NORMALIZE_WHITESPACE)
# Rref = 1.0
# E0 = 126.
# T = 293.15
# resp = 2.0
# print(lloyd_fix(T, Rref, E0))
# #1.40590910521
# print(lloyd_fix_p(T, [Rref, E0]))
# #1.40590910521
# print(cost_lloyd_fix([Rref, E0], T, resp))
# #0.59409089479
# print(cost2_lloyd_fix([Rref, E0], T, resp))
# #0.352943991272
# print(poly(T,2,1))
# #295.15
# print(poly_p(T,[2,1]))
# #295.15