Fast pairwise SameQ test of an ordered list

Question

I have an ordered list of $N$ numbers, and I want to find out whether adjacent numbers are SameQ or not, returning a list of 0s and 1s of length $N-1$. One approach is something like:

pairwiseBoole[y_] := Boole[SameQ @@@ Transpose[{Rest[y],Most[y]}]]

but this is extremely slow:

data = Union @ RandomReal[E + {0, 10^-7}, 10^6];
Length[data]

res = pairwiseBoole[data]; //AbsoluteTiming
Count[res, 0]

997747

{0.640326, Null}

993254

If I naively compile this, I can increase the speed by an order of magnitude, but I think it should be possible to write a function that is still faster, since:

Rest[y] - Most[y]; //AbsoluteTiming

{0.009215, Null}

and SameQ is true as long as they only differ in their last binary digit (note that one might need to know about Internal`$SameQTolerance).

How can I perform this computation as quickly as possible?

Update

Based on the answers/comments, I should have provided a better example, one where the fact that SameQ uses a relative and not absolute tolerance plays a more important role. So, the kind of data I'm using is more similar to the following (corrected to output a packed array thanks to @Szabolcs):

data = Developer`ToPackedArray @ Union[
    RandomChoice[100 Sin[Subdivide[0, 7, 40]], 10^6] +
    RandomReal[{0, 10^-7}, 10^6]
];

And, it is vital that the output is the same as what one would get using pairwise SameQ on list elements. Only results that produce the same output as pairwiseBoole are acceptable.

My compiled version

By the way, here is the compiled version that I had come up with:

pairwiseBooleC = Compile[{{x,_Real,1}},
    Table[Boole[SameQ[x[[i]], x[[i+1]]]], {i, Length[x]-1}]
];

and a comparison with my pairwiseBoole:

r1 = pairwiseBoole[data]; //RepeatedTiming
r2 = pairwiseBooleC[data]; //RepeatedTiming

r1 === r2

{0.63, Null}

{0.042, Null}

True

Some background, I am writing a RegionMember type function that can handle multiple regions, returning which region a point belongs to. In order to get edge cases correct it is vital that I use something equivalent to SameQ when comparing adjacent elements.

Answer 1

On my computer simple LibraryFunction, that compares subsequent elements in a loop, is fastest.

pairwiseBooleJkuczm = Last@Compile[{{data, _Real, 1}},
  Table[
    Boole[Compile`GetElement[data, i] === Compile`GetElement[data, i + 1]],
    {i, Length@data - 1}
  ],
  CompilationTarget -> "C",
  RuntimeOptions -> {"CatchMachineIntegerOverflow" -> False, "CompareWithTolerance" -> True}
];

SeedRandom@0;
data = Developer`ToPackedArray@
   Union[RandomChoice[100 Sin[Subdivide[0, 7, 40]], 10^6] + 
     RandomReal[{0, 10^-7}, 10^6]];

(res1 = pairwiseBoole@data) // MaxMemoryUsed // RepeatedTiming
(res2 = pairwiseBooleC@data) // MaxMemoryUsed // RepeatedTiming
(res3 = pairwiseBooleJkuczm@data) // MaxMemoryUsed // RepeatedTiming
res1 === res2 === res3
(* {0.563, 175791672} *)
(* {0.044,  15981232} *)
(* {0.013,   7990696} *)
(* True *)

Answer 2

Here is an approach assuming that the SameQ tolerance is not changed from its default.

\$MachineEpsilon

$MachineEpsilon is the smallest number that when added to 1. produces a different number. This means that 1. and 1.+ $MachineEpsilon will differ by 1 in the last bit of the base 2 representation of the mantissa. Hence the two numbers will satisfy the SameQ predicate. Any larger numbers will not. This suggests the following algorithm for deciding the SameQ predicate:

Algorithm

Let $l$ and $u$ be two numbers to be compared, with $l<u$. Then, the SameQ predicate should be equivalent to:

$$\frac{u - l}{\left|u\right|} \le \epsilon$$

We can rewrite this as:

$$\epsilon \left|u\right| + (l-u) \ge 0$$

Finally, if we use a Heaviside theta function where 0. is mapped to 1, then we can use the following arithmetical expression:

$$\theta(\epsilon \left|u\right| + (l-u))$$

which will return 1 if $l$ and $u$ are SameQ, and 0 otherwise. I'm not sure about the correctness of this algorithm, but it passes all of the tests that I could think of.

Mathematica implementation

Here is the Mathematica implementation:

pairwiseSameQ[data_] := With[{u=Rest[data], l=Most[data]},
    UnitStep[$MachineEpsilon Abs[u] + (l-u)]
]

Tests

Let's compare this function to my slow function:

data = Developer`ToPackedArray @ Union[
    RandomChoice[100 Sin[Subdivide[0,7,40]], 10^6] + 
    RandomReal[{0, 10^-7}, 10^6]
];

r1 = pairwiseBoole[data]; //RepeatedTiming
r2 = pairwiseSameQ[data]; //RepeatedTiming

r1 === r2

{0.651, Null}

{0.0109, Null}

True

We get the same results, and the function is faster than the Compile alternatives in the other answers. What happens if we Compile this approach?

Compile

Following @jkuczm's approach, a compiled version might look like:

pairwiseSameQC = With[{e=$MachineEpsilon},
    Last @ Compile[{{d, _Real, 1}},
        Table[
            Boole[e Abs[Compile`GetElement[d,i+1]]>=(Compile`GetElement[d,i+1]-Compile`GetElement[d,i])],
            {i, Length@d-1}
        ],
        CompilationTarget->"C",
        RuntimeOptions->"Speed"
    ]
];

One final comparison:

r3 = pairwiseSameQC[data]; //RepeatedTiming

r1 === r2 === r3

{0.0016, Null}

True

Now, that's pretty fast!

Answer 3

Best so far is BoolEval[Rest[y] == Most[y]]

---

data = Union @ RandomReal[E + {0, 10^-7}, 10^6];

---

pwb1[y_] := Boole[SameQ @@@ Transpose[{Rest[y], Most[y]}]]

pwb1[data]; // RepeatedTiming
(* {0.730, Null} *)

---

pwb2[y_] := Inner[Boole @* SameQ, Rest[y], Most[y], List]

pwb2[data]; // RepeatedTiming
(* {0.9117, Null} *)

---

pwb3[y_] := Boole @ Inner[SameQ, Rest[y], Most[y], List]

pwb3[data]; // RepeatedTiming
(* {0.66, Null} *)

---

pwb4[y_] := Boole[Developer`PartitionMap[SameQ, y, 2, 1]]

pwb4[data]; // RepeatedTiming
(* {0.55, Null} *)

---

pwb5[y_] := Unitize @ Differences[y]

pwb5[data]; // RepeatedTiming
(* {0.024, Null} *)

---

pwb6[y_] := 1 - Unitize @ Differences[y]

pwb6[data]; // RepeatedTiming
(* {0.0269, Null} *)

---

pwb7[y_] := 1 - Unitize[Rest[y] - Most[y]]

pwb7[data]; // RepeatedTiming
(* {0.0203, Null} *)

---

pwb8[y_] := 1 - Unitize[Rest[y] - Most[y], 1*^-17]

pwb8[data]; // RepeatedTiming
(* {0.0233, Null} *)

---

pwb9[y_] := ListConvolve[{1, 1}, y, {-1, 1}, {}, Times, SameQ]

pwb9[data]; // RepeatedTiming
(* {0.712, Null} *)

---

pwb10[y_] := BoolEval[Rest[y] == Most[y]]

pwb10[data]; // RepeatedTiming
(* {0.00671, Null} *)

---

pwb11[y_] := Subtract[1, Unitize[Subtract[Rest[y], Most[y]]]]

pwb11[data]; // RepeatedTiming
(* {0.0165, Null} *)