Finding time-series direction reversal of certain magnitude

I'd like to find the points of a time-series that are a certain distance away (in value, not in time) from the previous maximum, which I consider a reversal.

For example, for {1, 2, 3, 5, 10, 8, 6, 3} with a threshold of 4, the reversal point would be 6, since it's 4 units away from the previous maximum of 10.

I want to find all such reversal points, but at each one I need to reset the running maximum (or split the list and start again on the new list).

Here is some code to find the first reversal. It computes the distance from each element to a rolling maximum, then looks for the first difference bigger than the threshold.

ts = {1, 2, 3, 5, 10, 8, 6, 3};
rollmax = FoldList[Max, First[ts], Rest[ts]];
delta = rollmax - ts;
sel = Map[# >= 4 &, delta];
index = Position[sel, True, 1, 1]


I can do this in a procedural way, but it's probably not the preferred way. I'm new to functional programming and I don't quite know what sort of patterns are available that would help me here.

And if there was a better way of writing the above code, please let me know.

-

You seem to be on the right track. If I understand your question I believe this will help:

f = If[#2 + 4 <= #, -\[Infinity], Max[##]] &;

FoldList[f, First@ts, Rest@ts]

Position[%, -\[Infinity], {1}]


The above assumes that you want to reset the new maximum to the value after the reversal (3).
If you want to reset it to the reversal point value itself, try this:

f = If[#2 + 4 <= Max@#, {#2}, Max[##]] &;

FoldList[f, First@ts, Rest@ts]

Position[%, {_}, {1}]


I argue the superiority of the FoldList method over Module/MapIndexed. The latter introduces a variable that it does not need to, it is longer, and it is slower.

SeedRandom[1]
ts = RandomInteger[20, 50000];

Timing[
r1 =
Module[{max = -\[Infinity]},
MapIndexed[(max = Max[max, #1];
If[max - #1 >= 4, Sow[#2]; max = -\[Infinity];]) &, ts];] //
Reap // Last // Flatten;
]


{0.1966, Null}

f = If[#2 + 4 <= #, -\[Infinity], Max[##]] &;

Timing[
r2 = Join @@ Position[FoldList[f, #, {##2}] & @@ ts, -\[Infinity], 1];
]


{0.0874, Null}

r1 === r2


True

-
BTW, if you need to get the parts indexed by Position look to Extract; but I suspect if you are already programming with FoldList you know of this. –  Mr.Wizard Feb 14 '12 at 19:21
The only func programming I knew was map/reduce/filter and the x:xs concept from Haskell. I found FoldList when searching for something suitable in the Virtual Book. BTW, my first version looked like this FoldList[Max, 0, ts][[2;;]], but then I remembered x:xs –  Adal Feb 15 '12 at 1:44
@Adal that is an interesting avenue to Mathematica; some day I should learn Haskell. I argue that FoldList as you started with is superior to the answer you presently have selected. See my updated answer for comparative timings. –  Mr.Wizard Feb 15 '12 at 9:33

Here is another way of obtaining the positions of reversals using Reap-Sow and MapIndexed. I've used a longer ts than yours to demonstrate multiple reversals.

ts = {1, 2, 3, 5, 10, 8, 6, 3, 5, 7, 4, 3, 2, 1, 6, 9, 5};
Module[{max = -Infinity},
MapIndexed[
(max = Max[max, #1];If[max - #1 == 4, Sow[#2]; max = -Infinity;]) &, ts
];
] // Reap // Last // Flatten

Out[1]= {7, 12, 17}

-
Very interesting, I wasn't aware that you could change variables in an expression. Did you use Module strictly for avoiding polluting the global namespace, or for making it faster (if it is)? –  Adal Feb 15 '12 at 1:39
@Adal It was solely to localize the variable. –  rm -rf Feb 15 '12 at 3:15
For ts = {11, 10, 9, 7, 2, 4, 6, 9, 7, 5, 8, 9, 10, 11, 6, 3, 7};, only the first two positions returned are correct. –  Ｊ. Ｍ. Feb 15 '12 at 10:13
I think it's because of the strict check that the delta is 4, not bigger than 4. Another answer had the same problem. –  Adal Feb 15 '12 at 10:16
@J.M. Indeed, if you change it to >=4, it returns the correct positions. In fact, I didn't pay close attention to Adal's code, because Mr.Wizard had already answered with a FoldList solution and I didn't want to cloud my thoughts. I just went by "a certain distance away" which I took to mean equality. –  rm -rf Feb 15 '12 at 10:24

Here is a solution based on rules and recursion (not very efficient, but IMO rather interesting):

Clear[fn];
fn[lst_, delta_, startIndex_: 0] :=
Flatten@
ReplaceList[
lst,
{Shortest[PatternSequence[x__, y_, z___]], p_, q___} /;
y - p == delta && y >= Max[x] :>
({ # + startIndex, fn[{q}, delta, # + startIndex]} &[
Length[{x, y, z, p}]]),
1]


On the same test, it gives:

ts = {1, 2, 3, 5, 10, 8, 6, 3, 5, 7, 4, 3, 2, 1, 6, 9, 5};
fn[ts, 4]

(*
==> {7, 12, 17}
*)

-
Cool, this basically describes the shape of the pattern that I'm looking for. From the manual it seems that Shortest is not needed. What I don't understand is why using PatternSequence[x__, y_, z___, p_, q___] instead of PatternSequence[x__, y_, z___], p_, q___ doesn't work, but I've got more reading on patterns to do. –  Adal Feb 15 '12 at 2:12
fn[] fails on ts = {11, 10, 9, 7, 2, 4, 6, 9, 7, 5, 8, 9, 10, 11, 6, 3, 7};. –  Ｊ. Ｍ. Feb 15 '12 at 10:20
@J.M. I took it that the first element should not be considered a maximum. If it should, one should change x__ to x___, and it will work. I also used strict equality, while another possibility is to use inequality y-p>=delta. My point was to illustrate the idea though. –  Leonid Shifrin Feb 15 '12 at 13:01
@Adal You are right about the Shortest, it seems. As for why the difference, my guess is that the set of cases tested by the pattern-matcher is constructed differently for the two placements of PatternSequence, but alas no time at the moment to dig in and see why. –  Leonid Shifrin Feb 15 '12 at 13:24

Using the same example data as R.M.,

ts = {1, 2, 3, 5, 10, 8, 6, 3, 5, 7, 4, 3, 2, 1, 6, 9, 5};


Try something like this:

Flatten@Position[
FoldList[With[{b =
Boole[#2 + 4 <= First[#1]]}, {(1 - b) Max[{First[#1], #2}] +
b #2, b}] &, {First@ts, 0}, Rest@ts], {_, 1}]


{7, 12, 17}

This approach as the advantage that you can easily modify it to get the values of the troughs.

With[{m =
FoldList[
With[{b =
Boole[#2 + 4 <= First[#1]]}, {(1 - b) Max[{First[#1], #2}] +
b #2, b}] &, {First@ts, 0}, Rest@ts]},
Pick[m[[All, 1]], m[[All, 2]], 1] ]


{6, 3, 5}

Getting the values of the maxima is a little more involved from this starting point but this works:

Select[Last /@
SplitBy[FoldList[
With[{b =
Boole[#2 + 4 <= First[#1]]}, {(1 - b) Max[{First[#1], #2}] +
b #2, b}] &, {First@ts, 0}, Rest@ts], #[[2]] & ], #[[2]] == 0 &][[All, 1]]


{10, 7, 9}

-
The first snippet returns four positions, when it should return only three, on ts = {11, 10, 9, 7, 2, 4, 6, 9, 7, 5, 8, 9, 10, 11, 6, 3, 7};. –  Ｊ. Ｍ. Feb 15 '12 at 10:18
@J.M. You're right, I didn't test the case of consecutive falls larger than the threshold. The OP doesn't actually say how they want to handle that. Do they want the first fall or the last fall in a sequence from a peak? I'll have another look at how to handle that later. –  Verbeia Feb 15 '12 at 21:21

I might be a bit late, but hopefully not too late. Here's a routine for finding the positions of "reversals" based on some finite difference trickery, which exploits the analogy between differences for discrete data and derivatives in the usual calculus:

findReversalPositions[data_?VectorQ, h_?NumericQ] :=
Module[{n = Length[data], si = Sign[Differences[data]], ch},
ch = Flatten[Position[
ListConvolve[{1, 1}, si, {-1, 1}, 0, Times, Times], -1]] + 1;
ch = Switch[First[si], 1, Append[ch, n], -1, Prepend[ch, 1], 0, Return[]];
Composition[Extract[Apply[Range, First[#]], Last[#] + 1] &,
MapAt[First, #, {2}] &] /@ DeleteCases[{#,
Flatten[Position[(First[#] - Rest[#]) &[
Take[data, #]], _?(# >= h &)]]} & /@
Partition[ch, 2], {__, {}}]]


h is the "threshold" parameter, which is 4 in the OP. I haven't quite figured out what to do for the case of the first section of the data being constant (e.g. {9, 9, 5, ...}), so I decided on a hard Return[] in that case for the time being. (I'll edit this post once I figure out what to do. Also, the current version fixes a nasty bug in the previous version noted by kguler.)

Here are a few demonstrations:

(* OP's data *)
findReversalPositions[{1, 2, 3, 5, 10, 8, 6, 3}, 4]
{7}

(* R.M.'s data *)
findReversalPositions[{1, 2, 3, 5, 10, 8, 6, 3, 5, 7, 4, 3, 2, 1, 6, 9, 5}, 4]
{7, 12, 17}

(* R.M.'s data, flipped version *)
findReversalPositions[{11, 10, 9, 7, 2, 4, 6, 9, 7, 5, 8, 9, 10, 11, 6, 3, 7}, 4]
{4, 10, 15}

(* kguler's data, with different "thresholds" *)
findReversalPositions[{1, 2, 3, 8, 7, 6, 12, 8}, 2]
{6, 8}

findReversalPositions[{1, 2, 3, 8, 7, 6, 12, 8}, 4]
{8}

findReversalPositions[{1, 2, 3, 8, 7, 6, 12, 8}, 5]
{}

-
I get an error message on inputs like {1, 2, 3, 8, 7, 6, 12, 8} where a "reversal" does not happen after some local maxima. –  kguler Feb 15 '12 at 14:19
@kguler: indeed, findReversalPositions[{1, 2, 3, 8, 7, 6, 12, 8}, 3] breaks this method. I had not considered this case; let me think about how to fix this... –  Ｊ. Ｍ. Feb 15 '12 at 14:29