Detect Maxima Using Persistent Homology

Question

Inspired by this SO answer about peak detection.

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 code-golf, 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.

Answer 1

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!

Answer 2

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.

Answer 3

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

Answer 4

Python, 322 307 298 279 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

Answer 5

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.

Answer 6

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 += (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(key => -r(key)._1).map(i => r(i)._2).take(n)
}