import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
from scipy.signal import find_peaks # ピーク値を探す
import pandas as pd

import scienceplots # 科学論文スタイルのグラフ作成用
plt.style.use(['science', 'notebook']) # 論文向けのスタイルを有効化
#plt.rcParams['font.family'] = 'Times New Roman' # font familyの設定
plt.rcParams['font.family'] = 'serif' # font familyの設定
plt.rcParams['mathtext.fontset'] = 'cm' # math fontの設定

np.set_printoptions(precision=5)

def calf(func, x0, dx, nx): # 与えられた関数計算の出力
    print(f"{'x':>8s} {'f(x)':>15s}")
    for i in range(nx):
        x = x0 + dx*i
        y = func(x)
        print(f"{x:>8.3f} {y:>15.11f}")
    return

def func4a(x):
    y = -(x+1.5)*x*(x-1.5) 
    return y

print("f(x)=-(x+1.5)*x*(x-1.5)")
calf(func4a, 0.0, 0.2, 21)

f(x)=-(x+1.5)*x*(x-1.5)
       x            f(x)
   0.000   0.00000000000
   0.200   0.44200000000
   0.400   0.83600000000
   0.600   1.13400000000
   0.800   1.28800000000
   1.000   1.25000000000
   1.200   0.97200000000
   1.400   0.40600000000
   1.600  -0.49600000000
   1.800  -1.78200000000
   2.000  -3.50000000000
   2.200  -5.69800000000
   2.400  -8.42400000000
   2.600 -11.72600000000
   2.800 -15.65200000000
   3.000 -20.25000000000
   3.200 -25.56800000000
   3.400 -31.65400000000
   3.600 -38.55600000000
   3.800 -46.32200000000
   4.000 -55.00000000000

xx = np.linspace(-2.0, 2.0, 51)
yy = func4a(xx)
#xx,yy

df = pd.DataFrame({
    '$x$':xx,
    'func4a':yy,
})
df.style.background_gradient(cmap='rainbow')

locs, _ = find_peaks(yy) # ピーク値
locs, _, xx[locs],yy[locs]

(array([36]), {}, array([0.88]), array([1.29853]))

find_peaks?

Signature:
find_peaks(
    x,
    height=None,
    threshold=None,
    distance=None,
    prominence=None,
    width=None,
    wlen=None,
    rel_height=0.5,
    plateau_size=None,
)
Docstring:
Find peaks inside a signal based on peak properties.

This function takes a 1-D array and finds all local maxima by
simple comparison of neighboring values. Optionally, a subset of these
peaks can be selected by specifying conditions for a peak's properties.

Parameters
----------
x : sequence
    A signal with peaks.
height : number or ndarray or sequence, optional
    Required height of peaks. Either a number, ``None``, an array matching
    `x` or a 2-element sequence of the former. The first element is
    always interpreted as the  minimal and the second, if supplied, as the
    maximal required height.
threshold : number or ndarray or sequence, optional
    Required threshold of peaks, the vertical distance to its neighboring
    samples. Either a number, ``None``, an array matching `x` or a
    2-element sequence of the former. The first element is always
    interpreted as the  minimal and the second, if supplied, as the maximal
    required threshold.
distance : number, optional
    Required minimal horizontal distance (>= 1) in samples between
    neighbouring peaks. Smaller peaks are removed first until the condition
    is fulfilled for all remaining peaks.
prominence : number or ndarray or sequence, optional
    Required prominence of peaks. Either a number, ``None``, an array
    matching `x` or a 2-element sequence of the former. The first
    element is always interpreted as the  minimal and the second, if
    supplied, as the maximal required prominence.
width : number or ndarray or sequence, optional
    Required width of peaks in samples. Either a number, ``None``, an array
    matching `x` or a 2-element sequence of the former. The first
    element is always interpreted as the  minimal and the second, if
    supplied, as the maximal required width.
wlen : int, optional
    Used for calculation of the peaks prominences, thus it is only used if
    one of the arguments `prominence` or `width` is given. See argument
    `wlen` in `peak_prominences` for a full description of its effects.
rel_height : float, optional
    Used for calculation of the peaks width, thus it is only used if `width`
    is given. See argument  `rel_height` in `peak_widths` for a full
    description of its effects.
plateau_size : number or ndarray or sequence, optional
    Required size of the flat top of peaks in samples. Either a number,
    ``None``, an array matching `x` or a 2-element sequence of the former.
    The first element is always interpreted as the minimal and the second,
    if supplied as the maximal required plateau size.

    .. versionadded:: 1.2.0

Returns
-------
peaks : ndarray
    Indices of peaks in `x` that satisfy all given conditions.
properties : dict
    A dictionary containing properties of the returned peaks which were
    calculated as intermediate results during evaluation of the specified
    conditions:

    * 'peak_heights'
          If `height` is given, the height of each peak in `x`.
    * 'left_thresholds', 'right_thresholds'
          If `threshold` is given, these keys contain a peaks vertical
          distance to its neighbouring samples.
    * 'prominences', 'right_bases', 'left_bases'
          If `prominence` is given, these keys are accessible. See
          `peak_prominences` for a description of their content.
    * 'width_heights', 'left_ips', 'right_ips'
          If `width` is given, these keys are accessible. See `peak_widths`
          for a description of their content.
    * 'plateau_sizes', left_edges', 'right_edges'
          If `plateau_size` is given, these keys are accessible and contain
          the indices of a peak's edges (edges are still part of the
          plateau) and the calculated plateau sizes.

          .. versionadded:: 1.2.0

    To calculate and return properties without excluding peaks, provide the
    open interval ``(None, None)`` as a value to the appropriate argument
    (excluding `distance`).

Warns
-----
PeakPropertyWarning
    Raised if a peak's properties have unexpected values (see
    `peak_prominences` and `peak_widths`).

Warnings
--------
This function may return unexpected results for data containing NaNs. To
avoid this, NaNs should either be removed or replaced.

See Also
--------
find_peaks_cwt
    Find peaks using the wavelet transformation.
peak_prominences
    Directly calculate the prominence of peaks.
peak_widths
    Directly calculate the width of peaks.

Notes
-----
In the context of this function, a peak or local maximum is defined as any
sample whose two direct neighbours have a smaller amplitude. For flat peaks
(more than one sample of equal amplitude wide) the index of the middle
sample is returned (rounded down in case the number of samples is even).
For noisy signals the peak locations can be off because the noise might
change the position of local maxima. In those cases consider smoothing the
signal before searching for peaks or use other peak finding and fitting
methods (like `find_peaks_cwt`).

Some additional comments on specifying conditions:

* Almost all conditions (excluding `distance`) can be given as half-open or
  closed intervals, e.g., ``1`` or ``(1, None)`` defines the half-open
  interval :math:`[1, \infty]` while ``(None, 1)`` defines the interval
  :math:`[-\infty, 1]`. The open interval ``(None, None)`` can be specified
  as well, which returns the matching properties without exclusion of peaks.
* The border is always included in the interval used to select valid peaks.
* For several conditions the interval borders can be specified with
  arrays matching `x` in shape which enables dynamic constrains based on
  the sample position.
* The conditions are evaluated in the following order: `plateau_size`,
  `height`, `threshold`, `distance`, `prominence`, `width`. In most cases
  this order is the fastest one because faster operations are applied first
  to reduce the number of peaks that need to be evaluated later.
* While indices in `peaks` are guaranteed to be at least `distance` samples
  apart, edges of flat peaks may be closer than the allowed `distance`.
* Use `wlen` to reduce the time it takes to evaluate the conditions for
  `prominence` or `width` if `x` is large or has many local maxima
  (see `peak_prominences`).

.. versionadded:: 1.1.0

Examples
--------
To demonstrate this function's usage we use a signal `x` supplied with
SciPy (see `scipy.datasets.electrocardiogram`). Let's find all peaks (local
maxima) in `x` whose amplitude lies above 0.

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> from scipy.datasets import electrocardiogram
>>> from scipy.signal import find_peaks
>>> x = electrocardiogram()[2000:4000]
>>> peaks, _ = find_peaks(x, height=0)
>>> plt.plot(x)
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.plot(np.zeros_like(x), "--", color="gray")
>>> plt.show()

We can select peaks below 0 with ``height=(None, 0)`` or use arrays matching
`x` in size to reflect a changing condition for different parts of the
signal.

>>> border = np.sin(np.linspace(0, 3 * np.pi, x.size))
>>> peaks, _ = find_peaks(x, height=(-border, border))
>>> plt.plot(x)
>>> plt.plot(-border, "--", color="gray")
>>> plt.plot(border, ":", color="gray")
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.show()

Another useful condition for periodic signals can be given with the
`distance` argument. In this case, we can easily select the positions of
QRS complexes within the electrocardiogram (ECG) by demanding a distance of
at least 150 samples.

>>> peaks, _ = find_peaks(x, distance=150)
>>> np.diff(peaks)
array([186, 180, 177, 171, 177, 169, 167, 164, 158, 162, 172])
>>> plt.plot(x)
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.show()

Especially for noisy signals peaks can be easily grouped by their
prominence (see `peak_prominences`). E.g., we can select all peaks except
for the mentioned QRS complexes by limiting the allowed prominence to 0.6.

>>> peaks, properties = find_peaks(x, prominence=(None, 0.6))
>>> properties["prominences"].max()
0.5049999999999999
>>> plt.plot(x)
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.show()

And, finally, let's examine a different section of the ECG which contains
beat forms of different shape. To select only the atypical heart beats, we
combine two conditions: a minimal prominence of 1 and width of at least 20
samples.

>>> x = electrocardiogram()[17000:18000]
>>> peaks, properties = find_peaks(x, prominence=1, width=20)
>>> properties["prominences"], properties["widths"]
(array([1.495, 2.3  ]), array([36.93773946, 39.32723577]))
>>> plt.plot(x)
>>> plt.plot(peaks, x[peaks], "x")
>>> plt.vlines(x=peaks, ymin=x[peaks] - properties["prominences"],
...            ymax = x[peaks], color = "C1")
>>> plt.hlines(y=properties["width_heights"], xmin=properties["left_ips"],
...            xmax=properties["right_ips"], color = "C1")
>>> plt.show()
File:      ~/anaconda3/lib/python3.11/site-packages/scipy/signal/_peak_finding.py
Type:      function

fig = plt.figure() # グラフ領域の作成
plt.plot(xx,yy)
plt.plot([-10,10],[0,0],'--',color='gray',alpha=0.5)
plt.plot(xx[locs],yy[locs],'o') # ピーク点の表示
plt.xlim(-2.0, 2.0)
#plt.ylim(-4.0, 4.0)
plt.xlabel(r"$x$")
plt.ylabel(r"$f(x)$")
fig.savefig('p2_ex04_1.pdf')

def bisect(func, eps, x01, x02): # 二分法
    x1 = x01
    x2 = x02
    y = func(x1)
    while np.abs(x2-x1)>eps:
        xm = (x1+x2)*0.5
        if y*func(xm)>0.0:
            x1 = xm
        else:
            x2 = xm
    return xm

eps = 1.0e-15
x1, x2 = 0.1, 4.0
xzero = bisect(func4a, eps, x1, x2)
print ('zero of function: ',xzero)

zero of function:  1.4999999999999991

xzero2 = optimize.bisect(func4a, x1, x2)
print ('optimize.bisect : ',xzero2)

optimize.bisect :  1.4999999999991358

xzero3 = optimize.brentq(func4a, x1, x2)
print ('optimize.brentq : ',xzero3,"\n")

optimize.brentq :  1.5

#xx

uu = xx==0
uu

array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False,  True, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False])

def func4b(x):
    u = x==0
    t = np.where(u, 1.0, np.sin(x)/x)
    y = t + np.cos(x)
    return y

x=np.array([0,1])
u = x==0
t = np.where(u, 1.0, np.sin(x)/x)
u,t

/var/folders/69/4jdqr8157dl57hwmzlx1rd640000gn/T/ipykernel_75258/2911575268.py:3: RuntimeWarning: invalid value encountered in divide
  t = np.where(u, 1.0, np.sin(x)/x)

(array([ True, False]), array([1.     , 0.84147]))

def func4c(x):
    y = np.sinc(x/np.pi) + np.cos(x)
    return y

print("f(x)=sin(x)/x + cos(x)")
calf(func4b, -2.0, 0.2, 21)

f(x)=sin(x)/x + cos(x)
       x            f(x)
  -2.000   0.03850187687
  -1.800   0.31382436691
  -1.600   0.59553397960
  -1.400   0.87385980718
  -1.200   1.13905699278
  -1.000   1.38177329068
  -0.800   1.59340182297
  -0.600   1.76640640390
  -0.400   1.89460684977
  -0.200   1.97341323182
   0.000   2.00000000000
   0.200   1.97341323182
   0.400   1.89460684977
   0.600   1.76640640390
   0.800   1.59340182297
   1.000   1.38177329068
   1.200   1.13905699278
   1.400   0.87385980718
   1.600   0.59553397960
   1.800   0.31382436691
   2.000   0.03850187687

/var/folders/69/4jdqr8157dl57hwmzlx1rd640000gn/T/ipykernel_75258/1732970055.py:3: RuntimeWarning: invalid value encountered in scalar divide
  t = np.where(u, 1.0, np.sin(x)/x)

yy = func4b(xx)
df = pd.DataFrame({
    '$x$':xx,
    'func4b':yy,
})
df.style.background_gradient(cmap='coolwarm')

/var/folders/69/4jdqr8157dl57hwmzlx1rd640000gn/T/ipykernel_75258/1732970055.py:3: RuntimeWarning: invalid value encountered in divide
  t = np.where(u, 1.0, np.sin(x)/x)

xx2 = np.linspace(-4.0, 4.0, 81)
yy2 = func4b(xx2)
yy3 = func4c(xx2)
fig = plt.figure() # グラフ領域の作成
plt.plot(xx2,yy2)
plt.plot(xx2,yy3)
fig.savefig('p2_ex04_2.pdf')

/var/folders/69/4jdqr8157dl57hwmzlx1rd640000gn/T/ipykernel_75258/1732970055.py:3: RuntimeWarning: invalid value encountered in divide
  t = np.where(u, 1.0, np.sin(x)/x)

x1, x2 = 0.1, 3.0
xzero = bisect(func4b, eps, x1, x2)
print ('zero of function: ',xzero)
xzero2 = optimize.bisect(func4b, x1, x2)
print ('optimize.bisect : ',xzero2)
xzero3 = optimize.brentq(func4b, x1, x2)
print ('optimize.brentq : ',xzero3)

zero of function:  2.028757838110434
optimize.bisect :  2.028757838110187
optimize.brentq :  2.0287578381103835

	$x$	func4a
0	-2.000000	3.500000
1	-1.920000	2.757888
2	-1.840000	2.089504
3	-1.760000	1.491776
4	-1.680000	0.961632
5	-1.600000	0.496000
6	-1.520000	0.091808
7	-1.440000	-0.254016
8	-1.360000	-0.544544
9	-1.280000	-0.782848
10	-1.200000	-0.972000
11	-1.120000	-1.115072
12	-1.040000	-1.215136
13	-0.960000	-1.275264
14	-0.880000	-1.298528
15	-0.800000	-1.288000
16	-0.720000	-1.246752
17	-0.640000	-1.177856
18	-0.560000	-1.084384
19	-0.480000	-0.969408
20	-0.400000	-0.836000
21	-0.320000	-0.687232
22	-0.240000	-0.526176
23	-0.160000	-0.355904
24	-0.080000	-0.179488
25	0.000000	0.000000
26	0.080000	0.179488
27	0.160000	0.355904
28	0.240000	0.526176
29	0.320000	0.687232
30	0.400000	0.836000
31	0.480000	0.969408
32	0.560000	1.084384
33	0.640000	1.177856
34	0.720000	1.246752
35	0.800000	1.288000
36	0.880000	1.298528
37	0.960000	1.275264
38	1.040000	1.215136
39	1.120000	1.115072
40	1.200000	0.972000
41	1.280000	0.782848
42	1.360000	0.544544
43	1.440000	0.254016
44	1.520000	-0.091808
45	1.600000	-0.496000
46	1.680000	-0.961632
47	1.760000	-1.491776
48	1.840000	-2.089504
49	1.920000	-2.757888
50	2.000000	-3.500000

	$x$	func4b
0	-2.000000	0.038502
1	-1.920000	0.147249
2	-1.840000	0.257940
3	-1.760000	0.369965
4	-1.680000	0.482706
5	-1.600000	0.595534
6	-1.520000	0.707821
7	-1.440000	0.818936
8	-1.360000	0.928257
9	-1.280000	1.035165
10	-1.200000	1.139057
11	-1.120000	1.239344
12	-1.040000	1.335455
13	-0.960000	1.426845
14	-0.880000	1.512991
15	-0.800000	1.593402
16	-0.720000	1.667618
17	-0.640000	1.735214
18	-0.560000	1.795802
19	-0.480000	1.849035
20	-0.400000	1.894607
21	-0.320000	1.932256
22	-0.240000	1.961766
23	-0.160000	1.982966
24	-0.080000	1.995735
25	0.000000	2.000000
26	0.080000	1.995735
27	0.160000	1.982966
28	0.240000	1.961766
29	0.320000	1.932256
30	0.400000	1.894607
31	0.480000	1.849035
32	0.560000	1.795802
33	0.640000	1.735214
34	0.720000	1.667618
35	0.800000	1.593402
36	0.880000	1.512991
37	0.960000	1.426845
38	1.040000	1.335455
39	1.120000	1.239344
40	1.200000	1.139057
41	1.280000	1.035165
42	1.360000	0.928257
43	1.440000	0.818936
44	1.520000	0.707821
45	1.600000	0.595534
46	1.680000	0.482706
47	1.760000	0.369965
48	1.840000	0.257940
49	1.920000	0.147249
50	2.000000	0.038502

関数のピーク値，零点¶

ピーク値¶

二分法（bisection method）¶