# Detect Maxima Using Persistent Homology

## Background

Persistent homology is a fancy term for a fancy math concept, but in this case all we care about is how it can be used for peak detection. In the case of a 1-dimensional array we can describe the algorithm as follows:

We are given an array as input. We begin by sorting it, while maintaining references to the original indices for each value. We also create some kind of data structure to group points together into peaks. We then loop over each point in the sorted array. For each point we have three cases:

• If it is not adjacent to any points in existing groups, a new group is created with a "position" of the point's original index and given a "birth time" equal to that point's value.
• If it is adjacent to one point that is in a group, it is added to said group.
• If it is adjacent to two points in different groups, the younger group is given a "death time" equal to that point's value, and all its points are given to the older group.

At the end, the single group which remains is given a death time of the final value. Now we have a list of groups, each of which have a position that corresponds to a local maximum, and a persistence which is equal to their "birth time" minus their "death time".

## The challenge

Given a 1-dimensional sequence of inputs which support comparison and subtraction, and an integer $$\n \geq 1\$$, your task is to output the indices of the $$\n\$$ most persisent peaks, sorted by their persistence.

Your program does not need to follow the algorithm exactly as I've described it, nor does it need to be as efficient. You only need to support one class of inputs for the sequence (ie. integers, floats, etc.) and any reasonable format for I/O is allowed so long as it is consistent. You may assume that $$\n\$$ is less than or equal to the number of local maxima.

In the case that two local maxima to be merged have the same "birth time", the behavior for which one persists is undefined, and either is considered valid for this challenge. In the case that two or more maxima have the same persistence, they may appear in any order in the output.

This is , so the shortest answer in bytes wins.

## Worked example

Let's look at [1, 2, 1, 2, 3, 2, 1]. First we'll enumerate and sort [0:1, 1:2, 2:1, 3:2, 4:3, 5:2, 6:1] to [4:3, 1:2, 3:2, 5:2, 0:1, 2:1, 6:1] then we iterate over this. I'll show for each iteration what the aforementioned "group" data structure might look like, in this case simply a list of point indices, birth time and death time

[[[4], 3, None]]
[[[4], 3, None], [[1], 2, None]]
[[[4, 3], 3, None], [[1], 2, None]]
[[[4, 3, 5], 3, None], [[1], 2, None]]
[[[4, 3, 5], 3, None], [[1, 0], 2, None]]
[[[4, 3, 5, 1, 0, 2], 3, None], [[1, 0], 2, 1]]
[[[4, 3, 5, 1, 0, 2, 6], 3, 1], [[1, 0], 2, 1]]


this gives us the a peak at 4 with persistence 3-1=2 and a peak at 1 with persistence 2-1=1 thus our sorted list of peaks is [4,1]

## Test cases

Outputs here are 0-indexed.

[1, 2, 1, 2, 3, 2, 1, 2, 3, 4, 3, 2, 1]
2
->
[9, 4]

[1, 2, 3, 4, 5, 6, 7, 8, 7, 6, 5, 6, 7, 6, 5, 4, 3, 4, 5, 6, 5, 4, 5, 4, 3, 2, 1]
3
->
[7, 19, 12]

[93, 92, 89, 88, 89, 88, 87, 87, 88, 88, 90, 89, 88, 87, 86, 90, 89, 89, 88, 88, 86, 82, 76, 69, 61, 52, 44, 38, 33, 30, 27, 24, 22, 21, 18, 15, 10, 6, 2, 0, 0, 2, 2, 1, 2, 5, 8, 10, 13, 16]
4
->
[0, 49, 15, 10]


More information about persistent homology in general can be found on wikipedia and this more detailed write-up by the author of the SO answer.

• This needs a worked example or 2. Jun 27, 2023 at 21:50
• Could you provide the output for the last test case with all peak values? (i.e. with $n=7$, I think) Jun 28, 2023 at 10:00
• (It would have been simpler to just output all peaks sorted by their persistence, BTW. The argument $n$ seems a bit pointless.) Jun 28, 2023 at 10:05
• We begin by sorting it -> for consistent results, I think it's worth mentioning that this sort must be stable. (Or said otherwise, sorted by decreasing value first and increasing index second.) Jun 28, 2023 at 10:21
• In order, I added a worked example, hopefully it helps. One possible output for $n=7$ for the last test case would be [0, 49, 15, 10, 4, 41, 8], but note that swapping the 4 and 41 would be valid too. The point about excluding $n$ is fair, though I think changing it now would be a bit of a rug pull. Finally, I'm okay with output not being consistent in order to keep the requirements from being too tight, just like specifying behavior for ties, I think specifying sorting behavior would be too restrictive. Jun 28, 2023 at 14:21

# Python3, 350 bytes:

S=sorted
def f(s,n):
r,R={},0
for i,v in S(enumerate(s),key=lambda x:-x[1]):
R+=1
if(K:=[r.pop(I)for I in[*r]if i+1in r[I][2]or i-1in r[I][2]]):
*a,b=S(K,key=lambda x:x[0])
if a:r[R]=[a[0][0]-v,a[0][1],[]];b[2]+=a[0][2]+[i];R+=1;r[R]=b
else:b[2]+=[i];r[R]=b
else:r[R]=[v,i,[i]]
return[r[i][1]for i in S(r,key=lambda x:-r[x][0])][:n]


Try it online!

• instead of reversing the sorted lists, you can use negative values for the keys, i.e: key=lambda x:-x[1], also you don't need whitespace after the numbers on line 6
– c--
Jun 28, 2023 at 2:37
• @c-- Thanks, updated. Jun 28, 2023 at 2:45
• Nice answer, removed some conditionals, -71 bytes: Try it online! Jun 28, 2023 at 18:13
• @AnttiP Nice golf. I submitted another Python answer (with different algorithm), but I guess it's longer than your further-golfed version now (was shorter than the original 350 bytes). Jun 28, 2023 at 20:20
• Got a sudden urge to golf a few more bytes, -2 bytes: Try it online! Jun 30, 2023 at 15:00

# JavaScript (ES10*), 215 bytes

* ES10 guarantees a stable sort, but the OP has since clarified that this is not a strict requirement.

(or 214 bytes in ES12, but ATO doesn't seem to be working right now)

Expects (array)(n).

a=>n=>(o=[]).map(a=>a[2],a.map((v,i)=>[v,i]).sort(g=([a],[b])=>b-a).map(([v,i])=>([p,q]=o.sort(g).filter(a=>a[1].some(j=>(j-=i)*j<2)),q=q||[,[]],p)?q[p[1].push(...q[1],i),q[0]-=v,1]=[]:o.push([v,[i],i]))).slice(0,n)


Try it online!

### Commented

a =>                        // outer function taking a[]
n =>                        // inner function taking n
(o = [])                    // o[] = array of [ value, group, index ] tuples
.map(a =>                   // once o[] is finalized ...
a[2],                     //   ... keep only the original indices
a.map((v, i) => [ v, i ]) //   we first turn a[] into an array of
//   [ value, index ] pairs
.sort(g = ([a], [b]) =>   //   and sort them ...
b - a                   //     ... from highest to lowest value
)                         //
.map(([v, i]) =>          //   for each pair [ value, index ]:
( [p, q] = o.sort(g)    //     sort o[] by decreasing 'birth time'
.filter(a =>          //     and keep only the tuples ...
a[1].some(j =>      //       ... whose groups contain a value
(j -= i) * j < 2  //       adjacent to the current value
)                   //
),                    //     save at most 2 tuples in [ p, q ]
q = q || [, []],      //     force a default dummy value for q
p                     //     test p
) ?                     //     if p is defined:
q[                    //
p[1].push(          //       append to the group of p:
...q[1],          //         the values from the group of q
i                 //         followed by i
),                  //
q[0] -= v,          //       subtract the 'dead time' from the
1                   //       'birth time' of q
] = []                //       and clear the group of q
:                       //     else:
o.push([v, [i], i])   //       create a new tuple
)                         //   end of inner map()
)                           // end of outer map()
.slice(0, n)                // return only the n first elements

• Seems like this was upvoted less than 500ms after I hit the submit button. Either someone is really, really fast or it was a bot. Jun 28, 2023 at 10:14
• I was just staring the spec down, and trying to piece something together out of that, anything's impressive... Jun 28, 2023 at 10:16
• @UnrelatedString TBH, I had to stare at the spec for a while before figuring out what's going on here. I was a little more confused by the 3 links in the question which do provide a bit of context but also a lot of useless information. Jun 28, 2023 at 13:13

# Python, 322307298279 272 bytes

def P(l,n,*o):
L=len(l);*u,S=1,*[0]*L,1,sorted
for i in S(range(L),key=lambda x:-l[x]):*r,d=0,0,1;exec('k,c=i,sum(u[i-(d<0)::d])<2\nwhile~-u[k+1]:r[c]=max(r[c],l[i]-l[k]);k+=d\nd=-1;'*2);o+=((-r[0]or-r[1],i),)*(1>i%~-L*sum(u[i:i+3]));u[i+1]=1
return[*zip(*S(o))][1][:n]


Using different algorithm than the one in the question:

1. Sort descending
2. For each value, if it's a peak (immediate left and right have not been visited) iterate left and iterate right, record max diff until encountering a previously visited value or end of array. Give priority to encountering previously visited value. Store them.
3. Sort and return

Try it online!

-9 bytes thanks to UnrelatedString
-19 bytes thanks to AnttiP and Kevin
-7 bytes thanks to AnttiP

• 298 with some splat abuse Jun 29, 2023 at 1:39
• @UnrelatedString Ooh, interesting! Yea, I haven't used splat abuse in a while I forgot about it. And putting the o into the arguments, nice! Jun 29, 2023 at 5:48
• Shaved 8 bytes, might return to this later today: Try it online! Jun 29, 2023 at 15:29
• I'm playing both sides, so that I always come out on top Jun 29, 2023 at 18:56
• @AnttiP The \n'*2 can be ;'*2 for -1 byte: try it online. Jun 30, 2023 at 7:24

# Python, 153 bytes

lambda a,n:sorted(r:=range(N:=len(a)),key=lambda j:max((g:=lambda s,k:g(s,k+s)if k in r and a[k]<=a[j]else min(a[j:k%-N:s]or[0]))(1,j),g(-1,j))-a[j])[:n]


Attempt This Online!

Expects a sequence of positive integers (other positive numbers should work, too).

How?

I'm using the following algorithm which I think is equivalent to OP's:

For each position p find the nearest position with larger value and take the smallest value in between. Do this on the left and the right and keep the larger of the two minima minus the value at p. Sort positions by this quantity.

# Charcoal, 46 bytes

ＵＭθ⟦⁻⌈Ｅ⟦⮌…θ⊕κ✂θκ⟧∨⌊…λ⌕Ｅλ›νι¹¦⁰ικ⟧Ｆη⊞υ⌊⁻θυＩＥυ⊟ι


Try it online! Link is to verbose version of code. Explanation: Port of @Albert.Lang's Python answer, because it's the only post here I could understand.

ＵＭθ⟦⁻⌈Ｅ⟦⮌…θ⊕κ✂θκ⟧∨⌊…λ⌕Ｅλ›νι¹¦⁰ικ⟧


For each element of the input array, get the reversed prefix and suffix, truncate each at the first element greater than the current element, then take the minimum, or zero if the list is empty, subtract the current element from the larger of the two, and pair it with the current index.

Ｆη⊞υ⌊⁻θυ


Get the persistence and index of the n most persistent peaks.

ＩＥυ⊟ι


Output the indices.

# Scala, 516 511 bytes

Port of @Ajax1234's Python answer in Scala.

Saved many bytes thanks to the comment of @Aiden Chow

Golfed version. Try it online!

def f(s:Seq[Int],n:Int)={val r=mutable.Map[Int,(Int,Int,Seq[Int])]();var R=0;for((v,i)<-s.zipWithIndex.sortBy(vi=> -vi._1)){val K=r.keys.filter(I=>r(I)._3.contains(i+1)||r(I)._3.contains(i-1)).map(I=>r.remove(I).get).toSeq.sortBy(_._1);if(K.nonEmpty){val(a,b)=(K.init,K.last);R+=1;if(a.nonEmpty){r+=(R->(a.head._1-v,a.head._2,Nil));R+=1;r+=(R->(b._1,b._2,b._3++a.head._3++Seq(i)))}else{r+=(R->(b._1,b._2,b._3++Seq(i)))}}else{R+=1;r+=(R->(v,i,Seq(i)))}};r.keys.toSeq.sortBy(k=> -r(k)._1).map(i=>r(i)._2).take(n)}


Ungolfed version. Try it online!

def f(s:Seq[Int],n:Int)={
val r = mutable.Map[Int, (Int, Int, Seq[Int])]()
var R = 0
for ((v, i) <- (s.zipWithIndex.sortBy { case (v, i) => -v})) {
val K = r.keys
.filter(I => r(I)._3.contains(i + 1) || r(I)._3.contains(i - 1))
.map(I => r.remove(I).get)
.toSeq
.sortBy(_._1)
if (K.nonEmpty) {
val (a, b) = (K.init, K.last)
R += 1
if (a.nonEmpty) {
R += 1
r += (R -> (b._1, b._2, b._3 ++ a.head._3 ++ Seq(i)))
} else {
r += (R ->  (b._1, b._2, b._3 ++ Seq(i)))
}
} else {
R += 1
r += (R -> (v, i, Seq(i)))
}
}
r.keys.toSeq.sortBy(key => -r(key)._1).map(i => r(i)._2).take(n)
}

• Is this even fully golfed? I see a lot of whitespace here; I doubt that they are all required for the code to run correctly. Jun 28, 2023 at 1:33
• @Aiden Chow Thanks. I have updated. Jun 29, 2023 at 1:22
• That looks a lot better, but you should golf your code before posting it next time. Also in vi => -vi._1 the space after vi can be removed, so vi=> -vi._1. Jun 29, 2023 at 1:28
• There's also one part that has key=> -r(key)._1, I think that could just be k=> -r(k)._1. Jun 29, 2023 at 1:31