80
$\begingroup$

I want to find :

  • all local maxima in range
  • all local minima in range

From those points I can interpolate and combine functions upper and lower boundary. What I am really interested in, is the mean function of those boundaries.

enter image description here

Data model for this plot:

GetRLine3[MMStdata_, IO_: 1][x_: x] := ListInterpolation[#, InterpolationOrder -> IO, Method -> "Spline"][x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata); data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4], Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@ Range[-4, 4, 0.08]}]; xLimits = {Min@#1, Max@#1} & @@ Transpose[data]; f = D[GetRLine3[{data}, 3][x], x];

Edit: As my effort: 

minimums = DeleteDuplicates[Round[x /. Last[FindMinimum[f, {x, #}]] & /@ Transpose[data][[1]], 0.0001]] minimumvalues = (f /. x -> #)[[1]] & /@ minimums; minimumData := Transpose[{minimums, minimumvalues}]; maximums = DeleteDuplicates[Round[x /. Last[FindMaximum[f, {x, #}]] & /@ Transpose[data][[1]], 0.0001]]; maximumsvalues = (f /. x -> #)[[1]] & /@ maximums; maximumsData := Transpose[{maximums, maximumsvalues}]; maxf = Max[{GetRLine3[{maximumsData}, 3][x], f}] minf = Min[{GetRLine3[{minimumData}, 3][x], f}] mf = Mean[{maxf, minf}]

This was what I was trying to make: enter image description here

I still get quite few warnings and I'm sure it's not the best solution. I don't like the DeleteDuplicates@Round@ part, but it was necessarily to get the interpolation function working.

Improve this question
$\endgroup$
1
  • $\begingroup$ There's always a maximum in between two minima. So if you have found all minima (I don't know how to be sure of that, though), you can just look for a maximum in each of the intervals between subsequent minima (and for possibly two more on the left/right of the first/last minimum). $\endgroup$ Commented May 15, 2012 at 10:50

5 Answers 5

Reset to default
84
$\begingroup$

This can be done using event location within NDSolve. I start off as below (note f is slightly modified from what you have, mostly to rescale it).

GetRLine3[MMStdata_, IO_: 1][x_: x] := ListInterpolation[#, InterpolationOrder -> IO, Method -> "Spline"][ x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata); data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4], Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@ Range[-4, 4, 0.08]}]; xLimits = {Min@#1, Max@#1} & @@ Transpose[data]; f = First[100*D[GetRLine3[{data}, 3][x], x]];

We'll recapture f using NDSolve, and locate the points where the derivative vanishes in the process.

vals = Reap[ soln = y[x] /. First[NDSolve[{y'[x] == Evaluate[D[f, x]], y[-9.9] == (f /. x -> -9.9)}, y[x], {x, -9.9, 30}, Method -> {"EventLocator", "Event" -> y'[x], "EventAction" :> Sow[{x, y[x]}]}]]][[2, 1]];

Visual check:

Plot[f, {x, -9.9, 30}, Epilog -> {PointSize[Medium], Red, Point[vals]}]

enter image description here

Improve this answer
$\endgroup$
3
  • 2
    $\begingroup$ Unbeatable simplicity (+1). $\endgroup$ Commented May 15, 2012 at 16:21
  • 2
    $\begingroup$ I assume you forgot {Take[val, {2, -1, 2}],Drop[val, {2, -1, 2}]} and to differentiate maximum list from minimum list, we just need to compare first element. Not bad at all actually +1. (Not sure however, if this would work, if function has a segment with a constant growth of 0) $\endgroup$ Commented May 15, 2012 at 19:48
  • 2
    $\begingroup$ Very very cool! $\endgroup$ Commented Mar 19, 2015 at 13:37
38
$\begingroup$

Sounds like a job for Ted's RootSearch package.

Clear[f]; SeedRandom[1]; GetRLine3[MMStdata_, IO_: 1][x_: x] := ListInterpolation[#, InterpolationOrder -> IO, Method -> "Spline"][x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata); data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4], Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@ Range[-4, 4, 0.08]}]; xLimits = {Min@#1, Max@#1} & @@ Transpose[data]; f[x_] = First[D[GetRLine3[{data}, 3][x], x]];

Note that I changed the way yo defined f slightly.

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]; Needs["Ersek`RootSearch`"] {{a, b}} = InterpolatingFunctionDomain[Head[f[x]]]; criticalPoints = x /. RootSearch[f'[x] == 0, {x, a, b}]; mins = Select[criticalPoints, f''[#] > 0 &]; maxs = Select[criticalPoints, f''[#] < 0 &];

Well, I get an error so I'm not certain how well it worked, but it does look good.

Plot[f[x], {x, -9, 30}, Epilog ->{PointSize[Medium], Blue, Point[{#, f[#]} & /@ mins], Darker[Red], Point[{#, f[#]} & /@ maxs] }]]

enter image description here

Some comments

  • I didn't actually notice the importance of separating the maxima from the minima right away and have edited accordingly
  • Of course, it's feasible that $f'(x)=0$ without a max or min occurring at $x$, although I think that the random nature of your function makes the probability of this occurring zero.
  • Not sure we can really prove we've got all the extremes.

Edit in Response to Artes

If you use FindMaximum from two different starting values you will likely get two slightly different approximations to the same root. Here's an example with exactly one max.

g[x_] = Exp[-(Cos[x] - x)^2]; x1 = x /. Last[FindMaximum[g[x], {x, 0.8}]]; x2 = x /. Last[FindMaximum[g[x], {x, 0.7}]]; x1 - x2

8.35447*10^-9

Improve this answer
$\endgroup$
3
  • $\begingroup$ Evaluating e.g. FindMaximum[f, {x, #}] & /@ {-5, -4.5} we got 2 maxima, where the difference between them : Subtract @@ (#2[[1, 2]] & @@@ %)is about 1.33227*10^-14. Does RootSearch treat them as two different ones or as the only one local maximum ? $\endgroup$ Commented May 15, 2012 at 12:08
  • $\begingroup$ @Artes Due to the random nature of the function, I can't reproduce this behavior. I think it is more likely, though that your FindMaximum code is finding two slightly different approximations to the same max. See my edit for an example. $\endgroup$ Commented May 15, 2012 at 12:15
  • $\begingroup$ Indeed, it seems to be the same maximum. In your example we can set appropriately WorkingPrecision etc. to ensure the issue, however for some experimental or interpolated data it could be difficult do decide. $\endgroup$ Commented May 15, 2012 at 16:37
26
$\begingroup$

Since you're fitting your data with a cubic spline, and defining $f(x)$ as the derivative of the interpolation function, then if we solve for $f'(x)=0$ for the critical points, the resulting function is a piecewise linear function. We can pull the interpolation x-grid data and use it to construct linear functions which are easy to solve for the roots without resorting to the FindMinimum, etc...

Clear[f]; SeedRandom[1]; GetRLine3[MMStdata_, IO_: 1][x_: x] := ListInterpolation[#, InterpolationOrder -> IO,Method -> "Spline"] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@MMStdata); data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4], Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@Range[-4, 4, 0.08]}]; xLimits = {Min@#1, Max@#1} & @@ Transpose[data]; f = First[GetRLine3[{data}, 3][x]] Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"]; grid = First@InterpolatingFunctionCoordinates[f]; Show[{Plot[f''[x], {x, -9, 30}, PlotRange -> All],ListPlot[Thread[{grid, f''[grid]}], PlotStyle -> Red], Graphics[Map[{Dashed, Line[#]} &,Partition[Thread [{grid, f''[grid]}], 2, 1]]]}]

enter image description here

The dashed curve is the reconstruction from the interpolation grid points shown in red. Next, $\it{eqn}$ defines the linear parameterizations where we have an intersection with the $x-$axis if $0<t<1$ when $y=0$. Since the interpolation function is a polynomial, the minima and maxima alternate and can be pulled from the zeros list in order.

NB: The dashed curve is not exactly $f''(x)$ here, as the red interpolation sample points don't lie exactly on the critical points in the first figure. According to the documentation, taking the derivative of an InterpolatingFunction returns a new InterpolatingFunction, however, it seems the grid values are the same for both. I guess you can always find $f'''(x)=0$ to get those.

pts = Thread[{grid, f''[grid]}]; eqn = Map[#[[1]] (1 - t) + #[[2]] (t) &, Partition[pts, 2, 1]]; zeros = Select[ Flatten[Table[{x, t} /. Solve[eqn[[j]] == {x, 0}, {x, t}], {j, Length[eqn]}], 1], 0 <= #[[2]] <= 1 &][[All, 1]]; zerosM = zeros[[1 ;; -1 ;; 2]]; zerosN = zeros[[2 ;; -1 ;; 2]]; Plot[f'[x], {x, -9, 30}, Epilog -> {PointSize[Medium], Red, Point[{#, f'[#]} & /@ zerosM], Blue, Point[{#, f'[#]} & /@ zerosN]}, PlotRange -> All]

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ What a great answer. At first, it did not work with my real data and I assumed your code was improperly written - it was not the case. The bug was that I attached [x] at the end of a interpolation function and it messed up further computations. I can only imagine how much time you saved me, it's pity that I can +1 only once. $\endgroup$ Commented May 16, 2012 at 12:49
22
$\begingroup$

I just come to sell my new MeshFunctions powered "golden hammer" :)

GetRLine3[MMStdata_, IO_: 1][x_: x] := ListInterpolation[#, InterpolationOrder -> IO, Method -> "Spline"][ x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata); data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4], Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@ Range[-4, 4, 0.08]}]; xLimits = {Min@#1, Max@#1} & @@ Transpose[data]; f = Function[x, Evaluate[D[GetRLine3[{data}, 3][x], x][[1]]]]; funcline = Plot[f[x], {x, ##}, PlotStyle -> GrayLevel[.8]] & @@ xLimits; maximapts = Plot[f[x], {x, ##}, PlotStyle -> None, PlotRange -> All, MeshFunctions -> Function[{x, y}, f'[x]], Mesh -> {{0}}, MeshStyle -> Directive[AbsolutePointSize[4], Red], RegionFunction -> Function[{x, y}, f''[x] < 0] ] & @@ xLimits; minimapts = Plot[f[x], {x, ##}, PlotStyle -> None, PlotRange -> All, MeshFunctions -> Function[{x, y}, f'[x]], Mesh -> {{0}}, MeshStyle -> Directive[AbsolutePointSize[4], Blue], RegionFunction -> Function[{x, y}, f''[x] > 0] ] & @@ xLimits; Show[{funcline, maximapts, minimapts}]

styled minima and maxima on line

Update

To extract the found extreme points, we can do it like this:

Cases[maximapts, GraphicsComplex[ pts_, body__ ]:>( pts[[#]]& /@ Cases[{body}, Point[idxes_]:>idxes, ∞] ), ∞] 
{{{{1.43947,0.041135},{4.75921,0.130593},<<30>>,{21.4528,-0.000886152}}}}
Improve this answer
$\endgroup$
5
  • $\begingroup$ your method works very well for me. Can you please tell me how to obtain a list of the coordinates of the maxima and minima? $\endgroup$ Commented Nov 27, 2014 at 17:56
  • $\begingroup$ @Pincopallino Thanks! To extract the coordinates, you might want to check the InputForm of the graphics for clues. (e.g. Cases[maximapts,GraphicsComplex[pts_,body__]:>Column[Short/@{pts,Cases[{body},Point[idxes_]:>idxes,∞]},Frame->All],∞][[1]] -- you might want to re-enter the Point because of SE's formatting on linebreaks in comment. ) $\endgroup$ Commented Nov 27, 2014 at 18:14
  • $\begingroup$ @Pincopallino Please see my update in the answer. $\endgroup$ Commented Nov 27, 2014 at 18:32
  • $\begingroup$ Thank you very much for your prompt response! $\endgroup$ Commented Nov 29, 2014 at 17:23
  • $\begingroup$ @Pincopallino No problem! But I should add a warning that this method will work only where the points are differentiable. $\endgroup$ Commented Nov 29, 2014 at 18:14
10
$\begingroup$

Just for the record, if you are willing to (re)sample your function to have {x, y} pairs instead of a regular function or InterpolatingFunction, there are multiple methods available (see for example here and here):

Original function: 

getRLine3[MMStdata_, IO_: 1][x_: x] := ListInterpolation[#, InterpolationOrder -> IO, Method -> "Spline"][ x] & /@ (({{#[[1]]}, #[[2]]}) & /@ # & /@ MMStdata); data = Transpose[{# + RandomReal[]*0.1 & /@ Range[-10, 30, 0.4], Tanh[#] + (Sech[2 x - 0.5]/1.5 + 1.5) /. x -> # & /@ Range[-4, 4, 0.08]}]; f = Head@First@D[getRLine3[{data}, 3][x], x]; (* sample function at regular times *) times = Range[0, 30, .01]; vals = f[times];

Using FindPeaks:

{maxP, max} = Transpose@FindPeaks[vals]; {minP, min} = Transpose@FindPeaks[-vals]; Plot[f[x], {x, 0, 30}, Epilog -> {AbsolutePointSize[5], Green, Point@Transpose@{times[[maxP]], max}, Red, Point@Transpose@{times[[minP]], -min}}]

Mathematica graphics

Using MaxDetect and MinDetect:

{minD, minD} = {MinDetect@vals, MaxDetect@vals}; Plot[f[x], {x, 0, 30}, Epilog -> {AbsolutePointSize[5], Green, Point@Pick[Transpose@{times, vals}, maxD, 1], Red, Point@Pick[Transpose@{times, vals}, minD, 1] }]

Mathematica graphics

MaxDetect+MinDetect is at least an order of magnitude faster than FindPeaks but of course the latter offers smoothing capabilities and further methods to finetune detection.

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.