r''' A module for functions of two component power law.
.. codeauthor:: Yoshifumi Futaana
The function form is
.. math::
J(E) = \left\{
\begin{array}{ll}
k_0 E^{\gamma_0} & (0<E<E_0)\\
k_1 E^{\gamma_1} & (E>E_0)\\
\end{array}
\right.
where
.. math::
E_0 = \left(\frac{k_1}{k_0}\right)^{\frac{1}{\gamma_0-\gamma_1}}
The parameter should satisfy
.. math::
k_0 > 0 \\
k_1 > 0 \\
\gamma_1 < 0
>>> cls = BiPower(100, -0.5, 1000, -1)
>>> func = cls.mkfunc()
>>> func(1)
100.0
>>> func(10) == 100 / math.sqrt(10)
True
>>> func(100)
10.0
>>> func(1000)
1.0
'''
import math
import logging
_logger = logging.getLogger(__name__)
import numpy as np
import scipy.optimize
import scipy.integrate
[docs]class BiPower:
'''
'''
def __init__(self, k0, r0, k1, r1):
'''
'''
if k0 <= 0:
msg = 'k0 should be positive (%g)' % k0
_logger.error(msg)
raise RuntimeError(msg)
if k1 <= 0:
msg = 'k1 should be positive (%g)' % k1
_logger.error(msg)
raise RuntimeError(msg)
self.k0 = k0
self.r0 = r0
self.k1 = k1
self.r1 = r1
self.E0 = (k1 / k0) ** (1. / (r0 - r1))
[docs] def mkfunc(self):
''' Return the function: J(E)
'''
def func(E):
if E <= 0:
raise RuntimeError("Domain error. E should be positive.")
f0 = self.k0 * E ** (self.r0)
f1 = self.k1 * E ** (self.r1)
return min([f0, f1])
return func
[docs] def get_knee(self):
''' Return the knee (*E0*).
Return the knee point.
For example,
>>> k0 = 1000
>>> r0 = -1
>>> k1 = 10000
>>> r1 = -2
will make the knee at 10.
>>> bip = BiPower(k0, r0, k1, r1)
>>> print('%.1f' % bip.get_knee())
10.0
'''
return self.E0
[docs] def getRandom(self, n=10000):
'''
'''
raise NotImplementedError('BiPower.getRandom')
[docs] def integral(self, x0, x1):
''' Return the integral over [x0, x1].
Returned is a tuple, with the first element holding the estimated
value of the integral and the second element holding
an upper bound on the error.
This method just calls ``scipy.integrate.quad``.
Use scipy.integrate.Inf and -scipy.integrate.Inf for infinity.
'''
if x0 < 0 or x1 < 0:
raise ValueError(
'Invalid integral range. x0=%g x1=%g (>=0)' % (x0, x1))
if x0 > x1:
raise ValueError(
'Invalid integral range. x0=%g x1=%g (x0 < x1)' % (x0, x1))
f = self.mkfunc()
return scipy.integrate.quad(f, x0, x1)
[docs]def mkfunc(k0, r0, k1, r1):
'''
'''
kls = BiPower(k0, r0, k1, r1)
return kls.mkfunc()
[docs]def mklnfunc(log10k0, r0, log10k1, r1):
''' Creating a function from with parameters in log space.
Parameter is log10(k0), r0, log10(k1) and r1.
Input/Output for the generated functions is in linear space.
A bipwer function, f0, is made with linear parameters.
>>> f0 = mkfunc(100, -1, 1000, -2)
>>> y0 = np.array([f0(x) for x in (0.1, 1, 10, 100)])
A bipwer function, f1, is made with log parameres.
Parameter k0 (=100) and k1 (=1000) is in log space, i.e. 2 and 3,
respectively.
>>> f1 = mklnfunc(2, -1, 3, -2)
Note that the generated function, f1, will take the argument (x) in
the linear space, thus, the argument should be the same for f0 and f1.
>>> y1 = np.array([f1(x) for x in (0.1, 1, 10, 100)])
f0 and f1 should be the same.
>>> (y0 == y1).all()
True
'''
kls = BiPower(10 ** log10k0, r0, 10 ** log10k1, r1)
return kls.mkfunc()
[docs]def get_knee(k0, r0, k1, r1):
''' Return the knee (`E0`).
See also :func:`BiPower.get_knee`.
>>> print('%.1f' % get_knee(1000, -1, 10000, -2))
10.0
'''
kls = BiPower(k0, r0, k1, r1)
return kls.get_knee()
[docs]def integral(k0, r0, k1, r1, x0, x1):
''' Return the integral over [x0, x1].
See also :func:`BiPower.integral`.
'''
kls = BiPower(k0, r0, k1, r1)
return kls.integral(x0, x1)
[docs]def getRandom(E0, r0, E1, r1, n=10000):
'''
'''
kls = BiPower(k0, r0, k1, r1)
return kls.getRandom(n=n)
[docs]def fitting(x_array, y_array, initprm=None):
r''' Return fitted parameters of k0, k1, r0, and r1.
As a sample, we try to fit the data that follows exactly the bi-linear.
.. math::
f(x) = 1000 \times x^{-1} ~~~\mathrm{for}~x < 10 \\
f(x) = 10000 \times x^{-2} ~~~\mathrm{for}~x >= 10
>>> x_array = [0.1, 1, 10, 100, 1000]
>>> y_array = [10000, 1000, 100, 1, 0.01]
>>> prm = fitting(x_array, y_array, initprm=(3, -1, 4, -2))
>>> prm = fitting(x_array, y_array)
:param x_array: Array of x.
:type x_array: ``numpy.array`` type. **Not MaskedArray**.
You have to use compressed() method to make new array
from ``MaskedArray``.
:param y_array: Array of y.
:type y_array: ``numpy.array`` type. **Not MaskedArray**.
You have to use compressed() method to make new array
from ``MaskedArray``.
:keyword initprm: You can specify the initial values for iteration.
The parameter is passed to ``scipy.optimize.leastsq`` function.
*initprm* should be a tuple of (log10k0, r0, log10k1, r1) or ``None``.
:return: (k0, r0, k1, r1), success_code
'''
x_array0 = np.array(x_array).copy()
y_array0 = np.array(y_array).copy()
# x_array and y_array should not be masked array
klsma = np.ma.masked_array([]).__class__
if klsma == x_array.__class__ or klsma == y_array.__class__:
raise ValueError(
'Type of x_array or y_array should not be masked array.')
# All y should be positive.
if y_array0.min() <= 0:
raise ValueError(
'All y_array elements should be positive. (min=%g)' %
y_array0.min())
# Number of element should be the same
nx = x_array0.size
ny = y_array0.size
if nx != ny:
raise ValueError('x_array size (%d) != y_array size (%d)' % (nx, ny))
if nx < 4:
raise ValueError('array size (%d) must be >= 4' % (nx))
_logger.debug('x_array0 = %s' % repr(x_array0))
_logger.debug('y_array0 = %s' % repr(y_array0))
ln_x = np.ma.log10(x_array0)
ln_y = np.ma.log10(y_array0)
_logger.debug('ln_x = %s' % repr(ln_x))
_logger.debug('ln_y = %s' % repr(ln_y))
# Now the part of fitting.
# Fitting should be done in log-space.
def fitfunc(p, x):
# p[0] is ln(k0), p[1] is r0, p[2] is ln(k1) and p[3] is r1
func = mklnfunc(p[0], p[1], p[2], p[3]) # This should not be used.
y = np.array([func(xi) for xi in x])
return np.log10(y)
def errfunc(p, x, y):
val = fitfunc(p, x) - np.log10(y)
return fitfunc(p, x) - np.log10(y)
if initprm is None:
# Estimate the initial values.
# initprm = (3, -1, 4, -2)
# Simplest way is using (x0, y0) and (x1, y1)
r0 = (ln_y[0] - ln_y[1]) / (ln_x[0] - ln_x[1])
r1 = (ln_y[-2] - ln_y[-1]) / (ln_x[-2] - ln_x[-1])
_logger.info('Estimated r0 = %g' % r0)
_logger.info('Estimated r1 = %g' % r1)
lnk0 = (ln_y[0] - ln_x[0] * r0)
lnk1 = (ln_y[-1] - ln_x[-1] * r1)
_logger.info('Estimated log10k0 = %g' % lnk0)
_logger.info('Estimated log10k1 = %g' % lnk1)
initprm = (lnk0, r0, lnk1, r1)
# Estimate of the initial values
p1, success = scipy.optimize.leastsq(errfunc, initprm,
args=(x_array0, y_array0))
return (10 ** p1[0], p1[1], 10 ** p1[2], p1[3]), success
import unittest
import doctest
[docs]def doctests():
return unittest.TestSuite((
doctest.DocTestSuite(),
))
if __name__ == '__main__':
# Sample program for plotting.
import numpy as np
bipower = mkfunc(100, -0.5, 1000, -1) # f(E) = 100 * E^-0.5 or 1000* E^-1
x = np.arange(-3, 6.001, 0.1)
x = 10. ** x
y = np.array([bipower(E) for E in x])
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(211)
ax.plot(x, y)
ax.set_yscale('log')
ax2 = fig.add_subplot(212)
ax2.plot(x, y)
ax2.set_xscale('log')
ax2.set_yscale('log')
fig.savefig('bipower.png')
# For unit test.
unittest.main(defaultTest='doctests')