# coding=utf-8 """This module contains all the response functions available in Pastas. """ import numpy as np from pandas import DataFrame from scipy.integrate import quad from scipy.special import gammainc, gammaincinv, k0, exp1, erfc, lambertw, \ erfcinv __all__ = ["Gamma", "Exponential", "Hantush", "Polder", "FourParam", "DoubleExponential", "One", "Edelman", "HantushWellModel"] [docs]class RfuncBase: _name = "RfuncBase" [docs] def __init__(self, up, meanstress, cutoff): self.up = up # Completely arbitrary number to prevent division by zero if 1e-8 > meanstress > 0: meanstress = 1e-8 elif meanstress < 0 and up is True: meanstress = meanstress * -1 self.meanstress = meanstress self.cutoff = cutoff [docs] def get_init_parameters(self, name): """Get initial parameters and bounds. It is called by the stressmodel. Parameters ---------- name : str Name of the stressmodel Returns ------- parameters : pandas DataFrame The initial parameters and parameter bounds used by the solver """ pass [docs] def get_tmax(self, p, cutoff=None): """Method to get the response time for a certain cutoff Parameters ---------- p: array_like array_like object with the values as floats representing the model parameters. cutoff: float, optional float between 0 and 1. Returns ------- tmax: float Number of days when 99.9% of the response has effectuated, when the cutoff is chosen at 0.999. """ pass [docs] def step(self, p, dt=1, cutoff=None, maxtmax=None): """Method to return the step function. Parameters ---------- p: array_like array_like object with the values as floats representing the model parameters. dt: float timestep as a multiple of of day. cutoff: float, optional float between 0 and 1. maxtmax: int, optional Maximum timestep to compute the block response for. Returns ------- s: numpy.array Array with the step response. """ pass [docs] def block(self, p, dt=1, cutoff=None, maxtmax=None): """Method to return the block funtion. Parameters ---------- p: array_like array_like object with the values as floats representing the model parameters. dt: float timestep as a multiple of of day. cutoff: float, optional float between 0 and 1. maxtmax: int, optional Maximum timestep to compute the block response for. Returns ------- s: numpy.array Array with the block response. """ s = self.step(p, dt, cutoff, maxtmax) return np.append(s[0], np.subtract(s[1:], s[:-1])) [docs] def get_t(self, p, dt, cutoff, maxtmax=None): """Internal method to detemine the times at which to evaluate the step- response, from t=0 Parameters ---------- p: array_like array_like object with the values as floats representing the model parameters. dt: float timestep as a multiple of of day. cutoff: float float between 0 and 1, that determines which part of the step- response is taken into account. maxtmax: float, optional The maximum time of the response, usually set to the simulation length. Returns ------- t: numpy.array Array with the times """ if isinstance(dt, np.ndarray): return dt else: tmax = self.get_tmax(p, cutoff) if maxtmax is not None: tmax = min(tmax, maxtmax) tmax = max(tmax, 3 * dt) return np.arange(dt, tmax, dt) [docs]class Gamma(RfuncBase): """Gamma response function with 3 parameters A, a, and n. Parameters ---------- up: bool or None, optional indicates whether a positive stress will cause the head to go up (True, default) or down (False), if None the head can go both ways. meanstress: float mean value of the stress, used to set the initial value such that the final step times the mean stress equals 1 cutoff: float proportion after which the step function is cut off. default is 0.999. Notes ----- The impulse response function may be written as: .. math:: \\theta(t) = At^{n-1} e^{-t/a} """ _name = "Gamma" [docs] def __init__(self, up=True, meanstress=1, cutoff=0.999): RfuncBase.__init__(self, up, meanstress, cutoff) self.nparam = 3 [docs] def get_init_parameters(self, name): parameters = DataFrame( columns=['initial', 'pmin', 'pmax', 'vary', 'name']) if self.up: parameters.loc[name + '_A'] = (1 / self.meanstress, 1e-5, 100 / self.meanstress, True, name) elif self.up is False: parameters.loc[name + '_A'] = (-1 / self.meanstress, -100 / self.meanstress, -1e-5, True, name) else: parameters.loc[name + '_A'] = (1 / self.meanstress, np.nan, np.nan, True, name) # if n is too small, the length of response function is close to zero parameters.loc[name + '_n'] = (1, 0.1, 100, True, name) parameters.loc[name + '_a'] = (10, 0.01, 1e4, True, name) return parameters [docs] def get_tmax(self, p, cutoff=None): if cutoff is None: cutoff = self.cutoff return gammaincinv(p[1], cutoff) * p[2] [docs] def gain(self, p): return p[0] [docs] def step(self, p, dt=1, cutoff=None, maxtmax=None): t = self.get_t(p, dt, cutoff, maxtmax) s = p[0] * gammainc(p[1], t / p[2]) return s [docs]class Exponential(RfuncBase): """Exponential response function with 2 parameters: A and a. Parameters ---------- up: bool or None, optional indicates whether a positive stress will cause the head to go up (True, default) or down (False), if None the head can go both ways. meanstress: float mean value of the stress, used to set the initial value such that the final step times the mean stress equals 1 cutoff: float proportion after which the step function is cut off. default is 0.999. Notes ----- The impulse response function may be written as: .. math:: \\theta(t) = A e^{-t/a} """ _name = "Exponential" [docs] def __init__(self, up=True, meanstress=1, cutoff=0.999): RfuncBase.__init__(self, up, meanstress, cutoff) self.nparam = 2 [docs] def get_init_parameters(self, name): parameters = DataFrame( columns=['initial', 'pmin', 'pmax', 'vary', 'name']) if self.up: parameters.loc[name + '_A'] = (1 / self.meanstress, 1e-5, 100 / self.meanstress, True, name) elif self.up is False: parameters.loc[name + '_A'] = (-1 / self.meanstress, -100 / self.meanstress, -1e-5, True, name) else: parameters.loc[name + '_A'] = (1 / self.meanstress, np.nan, np.nan, True, name) parameters.loc[name + '_a'] = (10, 0.01, 1000, True, name) return parameters [docs] def get_tmax(self, p, cutoff=None): if cutoff is None: cutoff = self.cutoff return -p[1] * np.log(1 - cutoff) [docs] def gain(self, p): return p[0] [docs] def step(self, p, dt=1, cutoff=None, maxtmax=None): t = self.get_t(p, dt, cutoff, maxtmax) s = p[0] * (1.0 - np.exp(-t / p[1])) return s [docs]class HantushWellModel(RfuncBase): """ A special implementation of the Hantush well function for multiple wells. Parameters ---------- up: bool, optional indicates whether a positive stress will cause the head to go up (True, default) or down (False) meanstress: float mean value of the stress, used to set the initial value such that the final step times the mean stress equals 1 cutoff: float proportion after which the step function is cut off. Default is 0.999. Notes ----- The Hantush well function is explained in [hantush_1955]_, [veling_2010]_ and [asmuth_2008]_. The impulse response function may be written as: .. math:: \\theta(t) = \\frac{A}{t} \\exp(-t/a -b/t) .. math:: p[0] = A # TBD \\frac{1}{4 \\pi kD} .. math:: p[1] = a = cS .. math:: p[2] = b = 1^2 / (4 \\lambda^2) .. math:: p[3] = r \\text{(not optimized)} where :math:`\\lambda = \\sqrt{kDc}` The parameter r (distance from the well to the observation point) is passed as a known value, and is used to scale the response function. The optimized parameters are slightly different from the original Hantush implementation: - A: the parameter is the same as the original Hantush, except that the distance (r) is set to 1.0 - a = cS: stays the same - b = 1 / (4 * lambda): r is used internally to scale with distance - r: distance, not optimized but used to scale A and b """ _name = "HantushWellModel" [docs] def __init__(self, up=False, meanstress=1, cutoff=0.999): RfuncBase.__init__(self, up, meanstress, cutoff) self.nparam = 3 [docs] def get_init_parameters(self, name): parameters = DataFrame( columns=['initial', 'pmin', 'pmax', 'vary', 'name']) if self.up: parameters.loc[name + '_A'] = (1 / self.meanstress, 0, 100 / self.meanstress, True, name) elif self.up is False: parameters.loc[name + '_A'] = (-1 / self.meanstress, -100 / self.meanstress, 0, True, name) else: parameters.loc[name + '_A'] = (1 / self.meanstress, np.nan, np.nan, True, name) parameters.loc[name + '_a'] = (100, 1e-3, 1e4, True, name) parameters.loc[name + '_b'] = (1, 1e-4, 25, True, name) return parameters [docs] def get_tmax(self, p, cutoff=None): r = 1.0 if len(p) == 3 else p[3] # approximate formula for tmax if cutoff is None: cutoff = self.cutoff cS = p[1] rho = np.sqrt(4 * r ** 2 * p[2]) k0rho = k0(rho) return lambertw(1 / ((1 - cutoff) * k0rho)).real * cS [docs] @staticmethod def gain(p): r = 1.0 if len(p) == 3 else p[3] rho = np.sqrt(4 * r ** 2 * p[2]) return p[0] * k0(rho) [docs] def step(self, p, dt=1, cutoff=None, maxtmax=None): r = 1.0 if len(p) == 3 else p[3] cS = p[1] rho = np.sqrt(4 * r ** 2 * p[2]) k0rho = k0(rho) t = self.get_t(p, dt, cutoff, maxtmax) tau = t / cS tau1 = tau[tau < rho / 2] tau2 = tau[tau >= rho / 2] w = (exp1(rho) - k0rho) / (exp1(rho) - exp1(rho / 2)) F = np.zeros_like(tau) F[tau < rho / 2] = w * exp1(rho ** 2 / (4 * tau1)) - (w - 1) * exp1( tau1 + rho ** 2 / (4 * tau1)) F[tau >= rho / 2] = 2 * k0rho - w * exp1(tau2) + (w - 1) * exp1( tau2 + rho ** 2 / (4 * tau2)) return p[0] * F / 2 [docs]class Hantush(RfuncBase): """ The Hantush well function, using the standard A, a, b parameters Parameters ---------- up: bool or None, optional indicates whether a positive stress will cause the head to go up (True, default) or down (False), if None the head can go both ways. meanstress: float mean value of the stress, used to set the initial value such that the final step times the mean stress equals 1 cutoff: float proportion after which the step function is cut off. default is 0.999. Notes ----- The Hantush well function is explained in [hantush_1955]_, [veling_2010]_ and [asmuth_2008]_. The impulse response function may be written as: .. math:: \\theta(t) = K_0(\\sqrt(4b)) \\frac{A}{t} \\exp(-t/a - ab/t) .. math:: p[0] = A = \\frac{1}{2 \\pi kD} .. math:: p[1] = a = cS .. math:: p[2] = b = r^2 / (4 \\lambda^2) where :math:`\\lambda = \\sqrt{kDc}` References ---------- .. [hantush_1955] Hantush, M. S., & Jacob, C. E. (1955). Nonâsteady radial flow in an infinite leaky aquifer. Eos, Transactions American Geophysical Union, 36(1), 95-100. .. [veling_2010] Veling, E. J. M., & Maas, C. (2010). Hantush well function revisited. Journal of hydrology, 393(3), 381-388. .. [asmuth_2008] Von Asmuth, J. R., Maas, K., Bakker, M., & Petersen, J. (2008). Modeling time series of ground water head fluctuations subjected to multiple stresses. Ground Water, 46(1), 30-40. """ _name = "Hantush" [docs] def __init__(self, up=False, meanstress=1, cutoff=0.999): RfuncBase.__init__(self, up, meanstress, cutoff) self.nparam = 3 [docs] def get_init_parameters(self, name): parameters = DataFrame( columns=['initial', 'pmin', 'pmax', 'vary', 'name']) if self.up: parameters.loc[name + '_A'] = (1 / self.meanstress, 0, 100 / self.meanstress, True, name) elif self.up is False: parameters.loc[name + '_A'] = (-1 / self.meanstress, -100 / self.meanstress, 0, True, name) else: parameters.loc[name + '_A'] = (1 / self.meanstress, np.nan, np.nan, True, name) parameters.loc[name + '_a'] = (100, 1e-3, 1e4, True, name) parameters.loc[name + '_b'] = (1, 1e-6, 25, True, name) return parameters [docs] def get_tmax(self, p, cutoff=None): # approximate formula for tmax if cutoff is None: cutoff = self.cutoff cS = p[1] rho = np.sqrt(4 * p[2]) k0rho = k0(rho) return lambertw(1 / ((1 - cutoff) * k0rho)).real * cS [docs] def gain(self, p): return p[0] * k0(np.sqrt(4 * p[2])) [docs] def step(self, p, dt=1, cutoff=None, maxtmax=None): cS = p[1] rho = np.sqrt(4 * p[2]) k0rho = k0(rho) t = self.get_t(p, dt, cutoff, maxtmax) tau = t / cS tau1 = tau[tau < rho / 2] tau2 = tau[tau >= rho / 2] w = (exp1(rho) - k0rho) / (exp1(rho) - exp1(rho / 2)) F = np.zeros_like(tau) F[tau < rho / 2] = w * exp1(rho ** 2 / (4 * tau1)) - (w - 1) * exp1( tau1 + rho ** 2 / (4 * tau1)) F[tau >= rho / 2] = 2 * k0rho - w * exp1(tau2) + (w - 1) * exp1( tau2 + rho ** 2 / (4 * tau2)) return p[0] * F / 2 [docs]class Polder(RfuncBase): """The Polder function, using the standard A, a, b parameters Notes ----- The Polder function is explained in [polder]_. The impulse response function may be written as: .. math:: \\theta(t) = \\exp(-\\sqrt(4b)) \\frac{A}{t^{-3/2}} \\exp(-t/a -b/t) .. math:: p[0] = A = \\exp(-x/\\lambda) .. math:: p[1] = a = \\sqrt{\\frac{1}{cS}} .. math:: p[2] = b = x^2 / (4 \\lambda^2) where :math:`\\lambda = \\sqrt{kDc}` References ---------- .. [polder] G.A. Bruggeman (1999). Analytical solutions of geohydrological problems. Elsevier Science. Amsterdam, Eq. 123.32 """ _name = "Polder" [docs] def __init__(self, up=True, meanstress=1, cutoff=0.999): RfuncBase.__init__(self, up, meanstress, cutoff) self.nparam = 3 [docs] def get_init_parameters(self, name): parameters = DataFrame( columns=['initial', 'pmin', 'pmax', 'vary', 'name']) parameters.loc[name + '_A'] = (1, 0, 2, True, name) parameters.loc[name + '_a'] = (10, 0.01, 1000, True, name) parameters.loc[name + '_b'] = (1, 1e-6, 25, True, name) return parameters [docs] def get_tmax(self, p, cutoff=None): if cutoff is None: cutoff = self.cutoff _, a, b = p b = a * b x = np.sqrt(b / a) inverfc = erfcinv(2 * cutoff) y = (-inverfc + np.sqrt(inverfc ** 2 + 4 * x)) / 2 tmax = a * y ** 2 return tmax [docs] def gain(self, p): # the steady state solution of Mazure g = p[0] * np.exp(-np.sqrt(4 * p[2])) if not self.up: g = -g return g [docs] def step(self, p, dt=1, cutoff=None, maxtmax=None): t = self.get_t(p, dt, cutoff, maxtmax) A, a, b = p s = A * self.polder_function(np.sqrt(b), np.sqrt(t / a)) # / np.exp(-2 * np.sqrt(b)) if not self.up: s = -s return s [docs] @staticmethod def polder_function(x, y): s = 0.5 * np.exp(2 * x) * erfc(x / y + y) + \ 0.5 * np.exp(-2 * x) * erfc(x / y - y) return s [docs]class One(RfuncBase): """Instant response with no lag and one parameter d. Parameters ---------- up: bool or None, optional indicates whether a positive stress will cause the head to go up (True) or down (False), if None (default) the head can go both ways. meanstress: float mean value of the stress, used to set the initial value such that the final step times the mean stress equals 1 cutoff: float proportion after which the step function is cut off. default is 0.999. """ _name = "One" [docs] def __init__(self, up=None, meanstress=1, cutoff=0.999): RfuncBase.__init__(self, up, meanstress, cutoff) self.nparam = 1 [docs] def get_init_parameters(self, name): parameters = DataFrame( columns=['initial', 'pmin', 'pmax', 'vary', 'name']) if self.up: parameters.loc[name + '_d'] = ( self.meanstress, 0, np.nan, True, name) elif self.up is False: parameters.loc[name + '_d'] = ( -self.meanstress, np.nan, 0, True, name) else: parameters.loc[name + '_d'] = ( self.meanstress, np.nan, np.nan, True, name) return parameters [docs] def gain(self, p): return p[0] [docs] def step(self, p, dt=1, cutoff=None, maxtmax=None): if isinstance(dt, np.ndarray): return p[0] * np.ones(len(dt)) else: return p[0] * np.ones(1) [docs] def block(self, p, dt=1, cutoff=None, maxtmax=None): return p[0] * np.ones(1) [docs]class FourParam(RfuncBase): """Four Parameter response function with 4 parameters A, a, b, and n. Parameters ---------- up: bool or None, optional indicates whether a positive stress will cause the head to go up (True, default) or down (False), if None the head can go both ways. meanstress: float mean value of the stress, used to set the initial value such that the final step times the mean stress equals 1 cutoff: float proportion after which the step function is cut off. default is 0.999. Notes ----- The impulse response function may be written as: .. math:: \\theta(t) = At^{n-1} e^{-t/a -ab/t} If Fourparam.quad is set to True, this response function uses np.quad to integrate the Four Parameter response function, which requires more calculation time. """ _name = "FourParam" [docs] def __init__(self, up=True, meanstress=1, cutoff=0.999): RfuncBase.__init__(self, up, meanstress, cutoff) self.nparam = 4 self.quad = False [docs] def get_init_parameters(self, name): parameters = DataFrame( columns=['initial', 'pmin', 'pmax', 'vary', 'name']) if self.up: parameters.loc[name + '_A'] = (1 / self.meanstress, 0, 100 / self.meanstress, True, name) elif self.up is False: parameters.loc[name + '_A'] = (-1 / self.meanstress, -100 / self.meanstress, 0, True, name) else: parameters.loc[name + '_A'] = (1 / self.meanstress, np.nan, np.nan, True, name) parameters.loc[name + '_n'] = (1, -10, 10, True, name) parameters.loc[name + '_a'] = (10, 0.01, 5000, True, name) parameters.loc[name + '_b'] = (10, 1e-6, 25, True, name) return parameters [docs] @staticmethod def function(t, p): return (t ** (p[1] - 1)) * np.exp(-t / p[2] - p[2] * p[3] / t) [docs] def get_tmax(self, p, cutoff=None): if cutoff is None: cutoff = self.cutoff if self.quad: x = np.arange(1, 10000, 1) y = np.zeros_like(x) func = self.function(x, p) func_half = self.function(x[:-1] + 1 / 2, p) y[1:] = y[0] + np.cumsum(1 / 6 * (func[:-1] + 4 * func_half + func[1:])) y = y / quad(self.function, 0, np.inf, args=p)[0] return np.searchsorted(y, cutoff) else: t1 = -np.sqrt(3 / 5) t2 = 0 t3 = np.sqrt(3 / 5) w1 = 5 / 9 w2 = 8 / 9 w3 = 5 / 9 x = np.arange(1, 10000, 1) y = np.zeros_like(x) func = self.function(x, p) func_half = self.function(x[:-1] + 1 / 2, p) y[0] = 0.5 * (w1 * self.function(0.5 * t1 + 0.5, p) + w2 * self.function(0.5 * t2 + 0.5, p) + w3 * self.function(0.5 * t3 + 0.5, p)) y[1:] = y[0] + np.cumsum(1 / 6 * (func[:-1] + 4 * func_half + func[1:])) y = y / quad(self.function, 0, np.inf, args=p)[0] return np.searchsorted(y, cutoff) [docs] @staticmethod def gain(p): return p[0] [docs] def step(self, p, dt=1, cutoff=None, maxtmax=None): if self.quad: t = self.get_t(p, dt, cutoff, maxtmax) s = np.zeros_like(t) s[0] = quad(self.function, 0, dt, args=p)[0] for i in range(1, len(t)): s[i] = s[i - 1] + quad(self.function, t[i - 1], t[i], args=p)[ 0] s = s * (p[0] / (quad(self.function, 0, np.inf, args=p))[0]) return s else: t1 = -np.sqrt(3 / 5) t2 = 0 t3 = np.sqrt(3 / 5) w1 = 5 / 9 w2 = 8 / 9 w3 = 5 / 9 if dt > 0.1: step = 0.1 # step size for numerical integration tmax = max(self.get_tmax(p, cutoff), 3 * dt) t = np.arange(step, tmax, step) s = np.zeros_like(t) # for interval [0,dt] : s[0] = (step / 2) * \ (w1 * self.function((step / 2) * t1 + (step / 2), p) + w2 * self.function((step / 2) * t2 + (step / 2), p) + w3 * self.function((step / 2) * t3 + (step / 2), p)) # for interval [dt,tmax]: func = self.function(t, p) func_half = self.function(t[:-1] + step / 2, p) s[1:] = s[0] + np.cumsum(step / 6 * (func[:-1] + 4 * func_half + func[ 1:])) s = s * (p[0] / quad(self.function, 0, np.inf, args=p)[0]) return s[int(dt / step - 1)::int(dt / step)] else: t = self.get_t(p, dt, cutoff, maxtmax) s = np.zeros_like(t) # for interval [0,dt] Gaussian quadrate: s[0] = (dt / 2) * \ (w1 * self.function((dt / 2) * t1 + (dt / 2), p) + w2 * self.function((dt / 2) * t2 + (dt / 2), p) + w3 * self.function((dt / 2) * t3 + (dt / 2), p)) # for interval [dt,tmax] Simpson integration: func = self.function(t, p) func_half = self.function(t[:-1] + dt / 2, p) s[1:] = s[0] + np.cumsum(dt / 6 * (func[:-1] + 4 * func_half + func[ 1:])) s = s * (p[0] / quad(self.function, 0, np.inf, args=p)[0]) return s [docs]class DoubleExponential(RfuncBase): """Gamma response function with 3 parameters A, a, and n. Parameters ---------- up: bool or None, optional indicates whether a positive stress will cause the head to go up (True, default) or down (False), if None the head can go both ways. meanstress: float mean value of the stress, used to set the initial value such that the final step times the mean stress equals 1 cutoff: float proportion after which the step function is cut off. default is 0.999. Notes ----- The impulse response function may be written as: .. math:: \\theta(t) = A (1 - \\alpha) e^{-t/a_1} + A \\alpha e^{-t/a_2} """ _name = "DoubleExponential" [docs] def __init__(self, up=True, meanstress=1, cutoff=0.999): RfuncBase.__init__(self, up, meanstress, cutoff) self.nparam = 4 [docs] def get_init_parameters(self, name): parameters = DataFrame( columns=['initial', 'pmin', 'pmax', 'vary', 'name']) if self.up: parameters.loc[name + '_A'] = (1 / self.meanstress, 0, 100 / self.meanstress, True, name) elif self.up is False: parameters.loc[name + '_A'] = (-1 / self.meanstress, -100 / self.meanstress, 0, True, name) else: parameters.loc[name + '_A'] = (1 / self.meanstress, np.nan, np.nan, True, name) parameters.loc[name + '_alpha'] = (0.1, 0.01, 0.99, True, name) parameters.loc[name + '_a1'] = (10, 0.01, 5000, True, name) parameters.loc[name + '_a2'] = (10, 0.01, 5000, True, name) return parameters [docs] def get_tmax(self, p, cutoff=None): if cutoff is None: cutoff = self.cutoff if p[2] > p[3]: # a1 > a2 return -p[2] * np.log(1 - cutoff) else: # a1 < a2 return -p[3] * np.log(1 - cutoff) [docs] def gain(self, p): return p[0] [docs] def step(self, p, dt=1, cutoff=None, maxtmax=None): t = self.get_t(p, dt, cutoff, maxtmax) s = p[0] * (1 - ((1 - p[1]) * np.exp(-t / p[2]) + p[1] * np.exp(-t / p[3]))) return s [docs]class Edelman(RfuncBase): """The function of Edelman, describing the propagation of an instantaneous water level change into an adjacent half-infinite aquifer. Parameters ---------- up: bool or None, optional indicates whether a positive stress will cause the head to go up (True, default) or down (False), if None the head can go both ways. meanstress: float mean value of the stress, used to set the initial value such that the final step times the mean stress equals 1 cutoff: float proportion after which the step function is cut off. default is 0.999. Notes ----- The Edelman function is emplained in [5]_. The impulse response function may be written as: .. math:: \\text{unknown} It's parameters are: .. math:: p[0] = \\beta = \\frac{\\sqrt{\\frac{4kD}{S}}}{x} References ---------- .. [5] http://grondwaterformules.nl/index.php/formules/waterloop/peilverandering """ _name = "Edelman" [docs] def __init__(self, up=True, meanstress=1, cutoff=0.999): RfuncBase.__init__(self, up, meanstress, cutoff) self.nparam = 1 [docs] def get_init_parameters(self, name): parameters = DataFrame( columns=['initial', 'pmin', 'pmax', 'vary', 'name']) beta_init = 1.0 parameters.loc[name + '_beta'] = (beta_init, 0, 1000, True, name) return parameters [docs] def get_tmax(self, p, cutoff=None): if cutoff is None: cutoff = self.cutoff return 1. / (p[0] * erfcinv(cutoff * erfc(0))) ** 2 [docs] @staticmethod def gain(p): return 1. [docs] def step(self, p, dt=1, cutoff=None, maxtmax=None): t = self.get_t(p, dt, cutoff, maxtmax) s = erfc(1 / (p[0] * np.sqrt(t))) return s