<jupyter_start><jupyter_text>Table of Contents 1&nbsp;&nbsp;Short study of the Lempel-Ziv complexity1.1&nbsp;&nbsp;Short definition1.2&nbsp;&nbsp;Python implementation1.3&nbsp;&nbsp;Tests (1/2)1.4&nbsp;&nbsp;Cython implementation1.5&nbsp;&nbsp;Numba implementation1.6&nbsp;&nbsp;Tests (2/2)1.7&nbsp;&nbsp;Benchmarks1.8&nbsp;&nbsp;Complexity ?1.9&nbsp;&nbsp;Conclusion1.10&nbsp;&nbsp;(Experimental) Julia implementation1.11&nbsp;&nbsp;Ending notes  Short study of the Lempel-Ziv complexityIn this short [Jupyter notebook](https://www.Jupyter.org/) aims at defining and explaining the [Lempel-Ziv complexity](https://en.wikipedia.org/wiki/Lempel-Ziv_complexity).[I](http://perso.crans.org/besson/) will give examples, and benchmarks of different implementations.- **Reference:** <NAME> and <NAME>, *« On the Complexity of Finite Sequences »*, IEEE Trans. on Information Theory, January 1976, p. 75–81, vol. 22, n°1. ---- Short definitionThe Lempel-Ziv complexity is defined as the number of different substrings encountered as the stream is viewed from begining to the end.As an example:```python>>> s = '1001111011000010'>>> lempel_ziv_complexity(s)   1 / 0 / 01 / 1110 / 1100 / 00106```Marking in the different substrings, this sequence $s$ has complexity $\mathrm{Lempel}$-$\mathrm{Ziv}(s) = 6$ because $s = 1001111011000010 = 1 / 0 / 01 / 1110 / 1100 / 0010$.- See the page https://en.wikipedia.org/wiki/Lempel-Ziv_complexity for more details. Other examples:```python>>> lempel_ziv_complexity('1010101010101010')   1 / 0 / 103>>> lempel_ziv_complexity('1001111011000010000010')   1 / 0 / 01 / 1110 / 1100 / 0010 / 000 / 0107>>> lempel_ziv_complexity('100111101100001000001010')   1 / 0 / 01 / 1110 / 1100 / 0010 / 000 / 010 / 108``` ---- Python implementation<jupyter_code>def lempel_ziv_complexity(binary_sequence):
    """Lempel-Ziv complexity for a binary sequence, in simple Python code."""
    u, v, w = 0, 1, 1
    v_max = 1
    length = len(binary_sequence)
    complexity = 1
    while True:
        if binary_sequence[u + v - 1] == binary_sequence[w + v - 1]:
            v += 1
            if w + v >= length:
                complexity += 1
                break
        else:
            if v > v_max:
                v_max = v
            u += 1
            if u == w:
                complexity += 1
                w += v_max
                if w > length:
                    break
                else:
                    u = 0
                    v = 1
                    v_max = 1
            else:
                v = 1
    return complexity<jupyter_output><empty_output><jupyter_text>---- Tests (1/2)<jupyter_code>s = '1001111011000010'
lempel_ziv_complexity(s)  # 1 / 0 / 01 / 1110 / 1100 / 0010
%timeit lempel_ziv_complexity(s)
lempel_ziv_complexity('1010101010101010')  # 1 / 0 / 10
lempel_ziv_complexity('1001111011000010000010')  # 1 / 0 / 01 / 1110
lempel_ziv_complexity('100111101100001000001010')  # 1 / 0 / 01 / 1110 / 1100 / 0010 / 000 / 010 / 10
%timeit lempel_ziv_complexity('100111101100001000001010')<jupyter_output>19.4 µs ± 2.31 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)<jupyter_text>We can start to see that the time complexity of this function seems to grow exponentially as the complexity grows. ---- Cython implementationAs [this blog post](https://jakevdp.github.io/blog/2013/06/15/numba-vs-cython-take-2/) explains it, we can easily try to use [Cython](http://Cython.org/) in a notebook cell.> See [the Cython documentation](http://docs.cython.org/en/latest/src/quickstart/build.htmlusing-the-jupyter-notebook) for more information.<jupyter_code>%load_ext cython
%%cython
from __future__ import division
import cython

ctypedef unsigned int DTYPE_t

@cython.boundscheck(False) # turn off bounds-checking for entire function, quicker but less safe
def lempel_ziv_complexity_cython(str binary_sequence not None):
    """Lempel-Ziv complexity for a binary sequence, in simple Cython code (C extension)."""
    cdef DTYPE_t u = 0
    cdef DTYPE_t v = 1
    cdef DTYPE_t w = 1
    cdef DTYPE_t v_max = 1
    cdef DTYPE_t length = len(binary_sequence)
    cdef DTYPE_t complexity = 1
    # that was the only needed part, typing statically all the variables
    while True:
        if binary_sequence[u + v - 1] == binary_sequence[w + v - 1]:
            v += 1
            if w + v >= length:
                complexity += 1
                break
        else:
            if v > v_max:
                v_max = v
            u += 1
            if u == w:
                complexity += 1
                w += v_max
                if w > length:
                    break
                else:
                    u = 0
                    v = 1
                    v_max = 1
            else:
                v = 1
    return complexity<jupyter_output><empty_output><jupyter_text>Let try it!<jupyter_code>s = '1001111011000010'
lempel_ziv_complexity_cython(s)  # 1 / 0 / 01 / 1110 / 1100 / 0010
%timeit lempel_ziv_complexity_cython(s)
lempel_ziv_complexity_cython('1010101010101010')  # 1 / 0 / 10
lempel_ziv_complexity_cython('1001111011000010000010')  # 1 / 0 / 01 / 1110
lempel_ziv_complexity_cython('100111101100001000001010')  # 1 / 0 / 01 / 1110 / 1100 / 0010 / 000 / 010 / 10
%timeit lempel_ziv_complexity_cython('100111101100001000001010')<jupyter_output>259 ns ± 13.9 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)<jupyter_text>> $\implies$ Yay! It seems faster indeed! ---- Numba implementationAs [this blog post](https://jakevdp.github.io/blog/2013/06/15/numba-vs-cython-take-2/) explains it, we can also try to use [Numba](http://Numba.PyData.org/) in a notebook cell.<jupyter_code>from numba import jit
@jit("int32(boolean[:])")
def lempel_ziv_complexity_numba_x(binary_sequence):
    """Lempel-Ziv complexity for a binary sequence, in Python code using numba.jit() for automatic speedup (hopefully)."""
    u, v, w = 0, 1, 1
    v_max = 1
    length = len(binary_sequence)
    complexity = 1
    while True:
        if binary_sequence[u + v - 1] == binary_sequence[w + v - 1]:
            v += 1
            if w + v >= length:
                complexity += 1
                break
        else:
            if v > v_max:
                v_max = v
            u += 1
            if u == w:
                complexity += 1
                w += v_max
                if w > length:
                    break
                else:
                    u = 0
                    v = 1
                    v_max = 1
            else:
                v = 1
    return complexity

def str_to_numpy(s):
    """str to np.array of bool"""
    return np.array([int(i) for i in s], dtype=np.bool)

def lempel_ziv_complexity_numba(s):
    return lempel_ziv_complexity_numba_x(str_to_numpy(s))<jupyter_output><empty_output><jupyter_text>Let try it!<jupyter_code>str_to_numpy(s)
s = '1001111011000010'
lempel_ziv_complexity_numba(s)  # 1 / 0 / 01 / 1110 / 1100 / 0010
%timeit lempel_ziv_complexity_numba(s)
lempel_ziv_complexity_numba('1010101010101010')  # 1 / 0 / 10
lempel_ziv_complexity_numba('1001111011000010000010')  # 1 / 0 / 01 / 1110
lempel_ziv_complexity_numba('100111101100001000001010')  # 1 / 0 / 01 / 1110 / 1100 / 0010 / 000 / 010 / 10
%timeit lempel_ziv_complexity_numba('100111101100001000001010')<jupyter_output>9.35 µs ± 1.22 µs per loop (mean ± std. dev. of 7 runs, 100000 loops each)<jupyter_text>> $\implies$ Well... It doesn't seem that much faster from the naive Python code.> We specified the signature when calling [`@numba.jit`](http://numba.pydata.org/numba-doc/latest/user/jit.html), and used the more appropriate data structure (string is probably the smaller, numpy array are probably faster).> But even these tricks didn't help that much.> I tested, and without specifying the signature, the fastest approach is using string, compared to using lists or numpy arrays.> Note that the [`@jit`](http://numba.pydata.org/numba-doc/latest/user/jit.html)-powered function is compiled at runtime when first being called, so the signature used for the *first* call is determining the signature used by the compile function ---- Tests (2/2)To test more robustly, let us generate some (uniformly) random binary sequences.<jupyter_code>from numpy.random import binomial

def bernoulli(p, size=1):
    """One or more samples from a Bernoulli of probability p."""
    return binomial(1, p, size)
bernoulli(0.5, 20)<jupyter_output><empty_output><jupyter_text>That's probably not optimal, but we can generate a string with:<jupyter_code>''.join(str(i) for i in bernoulli(0.5, 20))
def random_binary_sequence(n, p=0.5):
    """Uniform random binary sequence of size n, with rate of 0/1 being p."""
    return ''.join(str(i) for i in bernoulli(p, n))
random_binary_sequence(50)
random_binary_sequence(50, p=0.1)
random_binary_sequence(50, p=0.25)
random_binary_sequence(50, p=0.5)
random_binary_sequence(50, p=0.75)
random_binary_sequence(50, p=0.9)<jupyter_output><empty_output><jupyter_text>And so, this function can test to check that the three implementations (naive, Cython-powered, Numba-powered) always give the same result.<jupyter_code>def tests_3_functions(n, p=0.5, debug=True):
    s = random_binary_sequence(n, p=p)
    c1 = lempel_ziv_complexity(s)
    if debug:
        print("Sequence s = {} ==> complexity C = {}".format(s, c1))
    c2 = lempel_ziv_complexity_cython(s)
    c3 = lempel_ziv_complexity_numba(s)
    assert c1 == c2 == c3, "Error: the sequence {} gave different values of the Lempel-Ziv complexity from 3 functions ({}, {}, {})...".format(s, c1, c2, c3)
    return c1
tests_3_functions(5)
tests_3_functions(20)
tests_3_functions(50)
tests_3_functions(500)
tests_3_functions(5000)<jupyter_output>Sequence s = 000011000001100001101000111111110101111100000101101001111110000001110110110000010111011011010111001011100100010100000110111101000010110010110001011000010000001101001100111111000101110001010000111000100101000111100111100111010111011100111100011101101001001111100100101100101001001001010010101101001111110000100111101010111100010100100001101010000001110010101010001010011010001111000011001011110100100010111111111100011110011000110100111010110011111101100111011000110010000010001101010010100000001010111111010101101111100000000010001110000000111101010011101010111111000001010111110001000000001101100110110111110011111010001100010111110110100111001001100100010110110010011001110000011001111010110000101000101111000100001110011110011000010100010100101001101011110010100110101100011000011000111100111011010001100000100100101111101010110000001011101100111001101100111111110100110011010111010100011100010000010000000011110100000000100001010110011110011011101001100001001100000100100111101110010[...]<jupyter_text>---- BenchmarksOn two example of strings (binary sequences), we can compare our three implementation.<jupyter_code>%timeit lempel_ziv_complexity('100111101100001000001010')
%timeit lempel_ziv_complexity_cython('100111101100001000001010')
%timeit lempel_ziv_complexity_numba('100111101100001000001010')
%timeit lempel_ziv_complexity('10011110110000100000101000100100101010010111111011001111111110101001010110101010')
%timeit lempel_ziv_complexity_cython('10011110110000100000101000100100101010010111111011001111111110101001010110101010')
%timeit lempel_ziv_complexity_numba('10011110110000100000101000100100101010010111111011001111111110101001010110101010')<jupyter_output>139 µs ± 4.33 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
1.94 µs ± 83.4 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
22 µs ± 364 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)<jupyter_text>Let check the time used by all the three functions, for longer and longer sequences:<jupyter_code>%timeit tests_3_functions(10, debug=False)
%timeit tests_3_functions(20, debug=False)
%timeit tests_3_functions(40, debug=False)
%timeit tests_3_functions(80, debug=False)
%timeit tests_3_functions(160, debug=False)
%timeit tests_3_functions(320, debug=False)
def test_cython(n):
    s = random_binary_sequence(n)
    c = lempel_ziv_complexity_cython(s)
    return c
%timeit test_cython(10)
%timeit test_cython(20)
%timeit test_cython(40)
%timeit test_cython(80)
%timeit test_cython(160)
%timeit test_cython(320)
%timeit test_cython(640)
%timeit test_cython(1280)
%timeit test_cython(2560)
%timeit test_cython(5120)
%timeit test_cython(10240)
%timeit test_cython(20480)<jupyter_output>52.1 ms ± 584 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
178 ms ± 6.03 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)<jupyter_text>---- Complexity ?$\implies$ The function `lempel_ziv_complexity_cython` seems to be indeed (almost) linear in $n$, the length of the binary sequence $S$.But let check more precisely, as it could also have a complexity of $\mathcal{O}(n \log n)$.<jupyter_code>import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
sns.set(context="notebook", style="darkgrid", palette="hls", font="sans-serif", font_scale=1.4)<jupyter_output><empty_output><jupyter_text>It's durty, but let us capture manually the times given by the experiments above.<jupyter_code>x = [10, 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 10240, 20480]
y = [18, 30, 55, 107, 205, 471, 977, 2270, 5970, 17300, 56600, 185000]
plt.figure()
plt.plot(x, y, 'o-')
plt.xlabel("Length $n$ of the binary sequence $S$")
plt.ylabel(r"Time in $\mu\;\mathrm{s}$")
plt.title("Time complexity of Lempel-Ziv complexity")
plt.show()
plt.figure()
plt.semilogx(x, y, 'o-')
plt.xlabel("Length $n$ of the binary sequence $S$")
plt.ylabel(r"Time in $\mu\;\mathrm{s}$")
plt.title("Time complexity of Lempel-Ziv complexity, semilogx scale")
plt.show()
plt.figure()
plt.semilogy(x, y, 'o-')
plt.xlabel("Length $n$ of the binary sequence $S$")
plt.ylabel(r"Time in $\mu\;\mathrm{s}$")
plt.title("Time complexity of Lempel-Ziv complexity, semilogy scale")
plt.show()
plt.figure()
plt.loglog(x, y, 'o-')
plt.xlabel("Length $n$ of the binary sequence $S$")
plt.ylabel(r"Time in $\mu\;\mathrm{s}$")
plt.title("Time complexity of Lempel-Ziv complexity, loglog scale")
plt.show()<jupyter_output><empty_output><jupyter_text>It is linear in $\log\log$ scale, so indeed the algorithm seems to have a linear complexity.To sum-up, for a sequence $S$ of length $n$, it takes $\mathcal{O}(n)$ basic operations to compute its Lempel-Ziv complexity $\mathrm{Lempel}-\mathrm{Ziv}(S)$. ---- Conclusion- The Lempel-Ziv complexity is not too hard to implement, and it indeed represents a certain complexity of a binary sequence, capturing the regularity and reproducibility of the sequence.- Using the [Cython](http://Cython.org/) was quite useful to have a $\simeq \times 100$ speed up on our manual naive implementation !- The algorithm is not easy to analyze, we have a trivial $\mathcal{O}(n^2)$ bound but experiments showed it is more likely to be $\mathcal{O}(n \log n)$ in the worst case, and $\mathcal{O}(n)$ in practice for "not too complicated sequences" (or in average, for random sequences). ---- (Experimental) [Julia](http://julialang.org) implementationI want to (quickly) try to see if I can use [Julia](http://julialang.org) to write a faster version of this function.See [issue 1](https://github.com/Naereen/Lempel-Ziv_Complexity/issues/1).**Disclaimer:** I am still learning Julia!<jupyter_code>%%time
%%script julia

"""Lempel-Ziv complexity for a binary sequence, in simple Julia code."""
function lempel_ziv_complexity(binary_sequence)
    u, v, w = 0, 1, 1
    v_max = 1
    size = length(binary_sequence)
    complexity = 1
    while true
        if binary_sequence[u + v] == binary_sequence[w + v]
            v += 1
            if w + v >= size
                complexity += 1
                break
            end
        else
            if v > v_max
                v_max = v
            end
            u += 1
            if u == w
                complexity += 1
                w += v_max
                if w > size
                    break
                else
                    u = 0
                    v = 1
                    v_max = 1
                end
            else
                v = 1
            end
        end
    end
    return complexity
end

s = "1001111011000010"
lempel_ziv_complexity(s)  # 1 / 0 / 01 / 1110 / 1100 / 0010

M = 100;
N = 10000;
for _ in 1:M
    s = join(rand(0:1, N));
    lempel_ziv_complexity(s);
end
lempel_ziv_complexity(s)  # 1 / 0 / 01 / 1110 / 1100 / 0010<jupyter_output>"Lempel-Ziv complexity for a binary sequence, in simple Julia code."
lempel_ziv_complexity (generic function with 1 method)
"1001111011000010"
6
100
10000
778
CPU times: user 0 ns, sys: 4 ms, total: 4 ms
Wall time: 4.86 s<jupyter_text>And to compare it fairly, let us use [Pypy](http://pypy.org) for comparison.<jupyter_code>%%time
%%pypy

def lempel_ziv_complexity(binary_sequence):
    """Lempel-Ziv complexity for a binary sequence, in simple Python code."""
    u, v, w = 0, 1, 1
    v_max = 1
    length = len(binary_sequence)
    complexity = 1
    while True:
        if binary_sequence[u + v - 1] == binary_sequence[w + v - 1]:
            v += 1
            if w + v >= length:
                complexity += 1
                break
        else:
            if v > v_max:
                v_max = v
            u += 1
            if u == w:
                complexity += 1
                w += v_max
                if w > length:
                    break
                else:
                    u = 0
                    v = 1
                    v_max = 1
            else:
                v = 1
    return complexity

s = "1001111011000010"
lempel_ziv_complexity(s)  # 1 / 0 / 01 / 1110 / 1100 / 0010

from random import random

M = 100
N = 10000
for _ in range(M):
    s = ''.join(str(int(random() < 0.5)) for _ in range(N))
    lempel_ziv_complexity(s)<jupyter_output>CPU times: user 4 ms, sys: 0 ns, total: 4 ms
Wall time: 7.89 s
{'fraction_non_alphanumeric': 0.10582453789021053, 'fraction_numerical': 0.14145306304797006, 'mean_word_length': 3.218610129564193, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 93, 'https://': 6, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 12, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '309946', 'license': 'MIT', 'n_tokens_mistral': 7896, 'n_tokens_neox': 5952, 'n_words': 1683, 'path': 'Short_study_of_the_Lempel-Ziv_complexity.ipynb', 'repo_name': 'UnofficialJuliaMirrorSnapshots/LempelZiv.jl-9cfb0768-07a8-5146-8897-ffaa113e453e', 'size': '17918'}