Finding extreme point by parabolic approximation¶
- PyAstronomy.pyasl.quadExtreme(x, y, mode='max', dp=(1, 1), exInd=None, fullOutput=False, fullPoint=False)¶
Find the extreme (minimum or maximum) by means of a parabolic fit.
This function searches for the maximum (or minimum) in the given ordinate values, fits a second-order polynomial to the surroundings of that point, and, finally, determines the thus approximated abscissa value of the extreme point.
- Parameters
- xarray
Abscissa values.
- yarray
Ordinate values.
- modestring, {min, max}, optional
Determines whether a minimum or a maximum is searched for. The default is a maximum.
- dptuple with two integers, optional
Determines the width around the extreme point, which will be used in the fitting. The default is one point to the right and left, i.e., dp = (1,1).
- extIndinteger, optional
If given, the function will assume that this is the index of the extreme point and not search for it.
- fullOutputboolean, optional
If True, the output will also cover the points used in the fitting and the resulting polynomial.
- fullPointboolean, optional
If True, the return value epos will be a tuple holding the abscissa and ordinate values of the extreme point. The default is False.
- Returns
- eposfloat or tuple
Position of the extreme point. If fullPoint is True, a tuple with the abscissa and ordinate values of the extreme point.
- miint
The index of the extreme point (maximum or minimum).
- xbarray, optional
Only returned if fullOutput is True. The abscissa values used in the polynomial fit. Note the the x-value of the extreme point has been subtracted.
- ybarray, optional
Only returned if fullOutput is True. The ordinate values used in the polynomial fit.
- pnumpy polynomial, optional
Only returned if fullOutput is True. The best-fit polynomial. Note that the fit refers to the xb axis, where the x-value of the extreme point has been subtracted.
Example¶
from __future__ import print_function, division
import numpy as np
from PyAstronomy import pyasl
import matplotlib.pylab as plt
# Create some data (a Gaussian)
x = np.arange(100.0)
y = np.exp(-(x-50.2714)**2/(2.*5.**2))
# Find the maximum
epos, mi = pyasl.quadExtreme(x, y, mode="max")
print("Maximum found at index: ", mi, ", value at maximum: ", y[mi])
print("Maximum found by parabolic fit: ", epos)
print()
# Find the maximum, use a wider range for the
# parabolic fit.
print("Using 5 points to each side of the maximum")
epos, mi = pyasl.quadExtreme(x, y, mode="max", dp=(5, 5))
print("Maximum found at index: ", mi, ", value at maximum: ", y[mi])
print("Maximum found by parabolic fit: ", epos)
print()
# Do as above, but get the full output
print("Using 2 points to each side of the maximum")
epos, mi, xb, yb, p = pyasl.quadExtreme(x, y, mode="max", dp=(2, 2), fullOutput=True)
# Evaluate polynomial at a number of points.
# Note that, internally, the x-value of the extreme point has
# been subtracted before the fit. Therefore, we need to re-shift
# it in the plot.
newx = np.linspace(min(xb), max(xb), 100)
model = np.polyval(p, newx)
# Plot the "data"
plt.plot(x, y, 'bp')
# Mark the points used in the fitting (shifted, because xb is shifted)
plt.plot(xb+x[mi], yb, 'rp')
# Overplot the model (shifted, because xb is shifted)
plt.plot(newx+x[mi], model, 'r--')
plt.show()
