Faster sample_entropy

raphaelvallat · Oct 21, 2018 · b20ae72 · b20ae72
1 parent c10ab43
commit b20ae72
Show file tree

Hide file tree

Showing 12 changed files with 153 additions and 33 deletions.
diff --git a/README.rst b/README.rst
@@ -70,7 +70,7 @@ Entropy
     0.9945519071575192
     0.8482185855709181
     2.0754913760787277
-    2.1919237573930315
+    2.192416747827227
 
 Fractal dimension
 -----------------
@@ -114,21 +114,21 @@ Here are some benchmarks computed on an average PC (i7-7700HQ CPU @ 2.80 Ghz - 8
     %timeit spectral_entropy(x, 100, method='fft')
     %timeit svd_entropy(x, order=3, delay=1)
     %timeit app_entropy(x, order=2) # Slow
-    %timeit sample_entropy(x, order=2) # Slow
+    %timeit sample_entropy(x, order=2) # Numba
     # Fractal dimension
     %timeit petrosian_fd(x)
     %timeit katz_fd(x)
-    %timeit higuchi_fd(x) # Numba (fast)
+    %timeit higuchi_fd(x) # Numba
     # Other
-    %timeit detrended_fluctuation(x) # Numba (fast)
+    %timeit detrended_fluctuation(x) # Numba
 
 .. parsed-literal::
 
     127 µs ± 3.86 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
     150 µs ± 859 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
     42.4 µs ± 306 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
     4.59 ms ± 62.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
-    4.61 ms ± 163 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
+    2.03 ms ± 39.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
     16.4 µs ± 251 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
     32.4 µs ± 578 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
     17.4 µs ± 274 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

diff --git a/docs/build/doctrees/environment.pickle b/docs/build/doctrees/environment.pickle
diff --git a/docs/build/doctrees/generated/entropy.sample_entropy.doctree b/docs/build/doctrees/generated/entropy.sample_entropy.doctree
diff --git a/docs/build/doctrees/index.doctree b/docs/build/doctrees/index.doctree
diff --git a/docs/build/html/_modules/entropy/entropy.html b/docs/build/html/_modules/entropy/entropy.html
@@ -87,7 +87,8 @@
 
   <h1>Source code for entropy.entropy</h1><div class="highlight"><pre>
 <span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
-<span class="kn">from</span> <span class="nn">math</span> <span class="k">import</span> <span class="n">factorial</span>
+<span class="kn">from</span> <span class="nn">numba</span> <span class="k">import</span> <span class="n">jit</span>
+<span class="kn">from</span> <span class="nn">math</span> <span class="k">import</span> <span class="n">factorial</span><span class="p">,</span> <span class="n">log</span>
 <span class="kn">from</span> <span class="nn">sklearn.neighbors</span> <span class="k">import</span> <span class="n">KDTree</span>
 <span class="kn">from</span> <span class="nn">scipy.signal</span> <span class="k">import</span> <span class="n">periodogram</span><span class="p">,</span> <span class="n">welch</span>
 
@@ -367,6 +368,61 @@ <h1>Source code for entropy.entropy</h1><div class="highlight"><pre>
     <span class="k">return</span> <span class="n">phi</span>
 
 
+<span class="nd">@jit</span><span class="p">(</span><span class="s1">&#39;f8(f8[:], i4, f8)&#39;</span><span class="p">,</span> <span class="n">nopython</span><span class="o">=</span><span class="kc">True</span><span class="p">)</span>
+<span class="k">def</span> <span class="nf">_numba_sampen</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">mm</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> <span class="n">r</span><span class="o">=</span><span class="mf">0.2</span><span class="p">):</span>
+    <span class="sd">&quot;&quot;&quot;</span>
+<span class="sd">    Fast evaluation of the sample entropy using Numba.</span>
+<span class="sd">    &quot;&quot;&quot;</span>
+    <span class="n">n</span> <span class="o">=</span> <span class="n">x</span><span class="o">.</span><span class="n">size</span>
+    <span class="n">n1</span> <span class="o">=</span> <span class="n">n</span> <span class="o">-</span> <span class="mi">1</span>
+    <span class="n">mm</span> <span class="o">+=</span> <span class="mi">1</span>
+    <span class="n">mm_dbld</span> <span class="o">=</span> <span class="mi">2</span> <span class="o">*</span> <span class="n">mm</span>
+
+    <span class="c1"># Define threshold</span>
+    <span class="n">r</span> <span class="o">*=</span> <span class="n">x</span><span class="o">.</span><span class="n">std</span><span class="p">()</span>
+
+    <span class="c1"># initialize the lists</span>
+    <span class="n">run</span> <span class="o">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">*</span> <span class="n">n</span>
+    <span class="n">run1</span> <span class="o">=</span> <span class="n">run</span><span class="p">[:]</span>
+    <span class="n">r1</span> <span class="o">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">*</span> <span class="p">(</span><span class="n">n</span> <span class="o">*</span> <span class="n">mm_dbld</span><span class="p">)</span>
+    <span class="n">a</span> <span class="o">=</span> <span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">*</span> <span class="n">mm</span>
+    <span class="n">b</span> <span class="o">=</span> <span class="n">a</span><span class="p">[:]</span>
+    <span class="n">p</span> <span class="o">=</span> <span class="n">a</span><span class="p">[:]</span>
+
+    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n1</span><span class="p">):</span>
+        <span class="n">nj</span> <span class="o">=</span> <span class="n">n1</span> <span class="o">-</span> <span class="n">i</span>
+
+        <span class="k">for</span> <span class="n">jj</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">nj</span><span class="p">):</span>
+            <span class="n">j</span> <span class="o">=</span> <span class="n">jj</span> <span class="o">+</span> <span class="n">i</span> <span class="o">+</span> <span class="mi">1</span>
+            <span class="k">if</span> <span class="nb">abs</span><span class="p">(</span><span class="n">x</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">-</span> <span class="n">x</span><span class="p">[</span><span class="n">i</span><span class="p">])</span> <span class="o">&lt;</span> <span class="n">r</span><span class="p">:</span>
+                <span class="n">run</span><span class="p">[</span><span class="n">jj</span><span class="p">]</span> <span class="o">=</span> <span class="n">run1</span><span class="p">[</span><span class="n">jj</span><span class="p">]</span> <span class="o">+</span> <span class="mi">1</span>
+                <span class="n">m1</span> <span class="o">=</span> <span class="n">mm</span> <span class="k">if</span> <span class="n">mm</span> <span class="o">&lt;</span> <span class="n">run</span><span class="p">[</span><span class="n">jj</span><span class="p">]</span> <span class="k">else</span> <span class="n">run</span><span class="p">[</span><span class="n">jj</span><span class="p">]</span>
+                <span class="k">for</span> <span class="n">m</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">m1</span><span class="p">):</span>
+                    <span class="n">a</span><span class="p">[</span><span class="n">m</span><span class="p">]</span> <span class="o">+=</span> <span class="mi">1</span>
+                    <span class="k">if</span> <span class="n">j</span> <span class="o">&lt;</span> <span class="n">n1</span><span class="p">:</span>
+                        <span class="n">b</span><span class="p">[</span><span class="n">m</span><span class="p">]</span> <span class="o">+=</span> <span class="mi">1</span>
+            <span class="k">else</span><span class="p">:</span>
+                <span class="n">run</span><span class="p">[</span><span class="n">jj</span><span class="p">]</span> <span class="o">=</span> <span class="mi">0</span>
+        <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">mm_dbld</span><span class="p">):</span>
+            <span class="n">run1</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">run</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
+            <span class="n">r1</span><span class="p">[</span><span class="n">i</span> <span class="o">+</span> <span class="n">n</span> <span class="o">*</span> <span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">run</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
+        <span class="k">if</span> <span class="n">nj</span> <span class="o">&gt;</span> <span class="n">mm_dbld</span> <span class="o">-</span> <span class="mi">1</span><span class="p">:</span>
+            <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">mm_dbld</span><span class="p">,</span> <span class="n">nj</span><span class="p">):</span>
+                <span class="n">run1</span><span class="p">[</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">run</span><span class="p">[</span><span class="n">j</span><span class="p">]</span>
+
+    <span class="n">m</span> <span class="o">=</span> <span class="n">mm</span> <span class="o">-</span> <span class="mi">1</span>
+
+    <span class="k">while</span> <span class="n">m</span> <span class="o">&gt;</span> <span class="mi">0</span><span class="p">:</span>
+        <span class="n">b</span><span class="p">[</span><span class="n">m</span><span class="p">]</span> <span class="o">=</span> <span class="n">b</span><span class="p">[</span><span class="n">m</span> <span class="o">-</span> <span class="mi">1</span><span class="p">]</span>
+        <span class="n">m</span> <span class="o">-=</span> <span class="mi">1</span>
+
+    <span class="n">b</span><span class="p">[</span><span class="mi">0</span><span class="p">]</span> <span class="o">=</span> <span class="n">n</span> <span class="o">*</span> <span class="n">n1</span> <span class="o">/</span> <span class="mi">2</span>
+    <span class="n">a</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="nb">float</span><span class="p">(</span><span class="n">aa</span><span class="p">)</span> <span class="k">for</span> <span class="n">aa</span> <span class="ow">in</span> <span class="n">a</span><span class="p">])</span>
+    <span class="n">b</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="nb">float</span><span class="p">(</span><span class="n">bb</span><span class="p">)</span> <span class="k">for</span> <span class="n">bb</span> <span class="ow">in</span> <span class="n">b</span><span class="p">])</span>
+    <span class="n">p</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">true_divide</span><span class="p">(</span><span class="n">a</span><span class="p">,</span> <span class="n">b</span><span class="p">)</span>
+    <span class="k">return</span> <span class="o">-</span><span class="n">log</span><span class="p">(</span><span class="n">p</span><span class="p">[</span><span class="o">-</span><span class="mi">1</span><span class="p">])</span>
+
+
 <div class="viewcode-block" id="app_entropy"><a class="viewcode-back" href="../../generated/entropy.app_entropy.html#entropy.app_entropy">[docs]</a><span class="k">def</span> <span class="nf">app_entropy</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">order</span><span class="o">=</span><span class="mi">2</span><span class="p">,</span> <span class="n">metric</span><span class="o">=</span><span class="s1">&#39;chebyshev&#39;</span><span class="p">):</span>
     <span class="sd">&quot;&quot;&quot;Approximate Entropy.</span>
 
@@ -439,8 +495,6 @@ <h1>Source code for entropy.entropy</h1><div class="highlight"><pre>
 
 <span class="sd">    Notes</span>
 <span class="sd">    -----</span>
-<span class="sd">    Original code from the mne-features package.</span>
-
 <span class="sd">    Sample entropy is a modification of approximate entropy, used for assessing</span>
 <span class="sd">    the complexity of physiological time-series signals. It has two advantages</span>
 <span class="sd">    over approximate entropy: data length independence and a relatively</span>
@@ -458,8 +512,11 @@ <h1>Source code for entropy.entropy</h1><div class="highlight"><pre>
 <span class="sd">    :math:`C(m, r)` is the number of embedded vectors of length</span>
 <span class="sd">    :math:`m` having a Chebyshev distance inferior to :math:`r`.</span>
 
-<span class="sd">    Code adapted from the mne-features package by Jean-Baptiste Schiratti</span>
-<span class="sd">    and Alexandre Gramfort.</span>
+<span class="sd">    Note that if metric == &#39;chebyshev&#39; and x.size &lt; 5000 points, then the</span>
+<span class="sd">    sample entropy is computed using a fast custom Numba script. For other</span>
+<span class="sd">    metric types or longer time-series, the sample entropy is computed using</span>
+<span class="sd">    a code from the mne-features package by Jean-Baptiste Schiratti</span>
+<span class="sd">    and Alexandre Gramfort (requires sklearn).</span>
 
 <span class="sd">    References</span>
 <span class="sd">    ----------</span>
@@ -487,10 +544,12 @@ <h1>Source code for entropy.entropy</h1><div class="highlight"><pre>
 <span class="sd">        &gt;&gt;&gt; print(sample_entropy(x, order=3, metric=&#39;euclidean&#39;))</span>
 <span class="sd">            2.725</span>
 <span class="sd">    &quot;&quot;&quot;</span>
-    <span class="n">phi</span> <span class="o">=</span> <span class="n">_app_samp_entropy</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">order</span><span class="o">=</span><span class="n">order</span><span class="p">,</span> <span class="n">metric</span><span class="o">=</span><span class="n">metric</span><span class="p">,</span> <span class="n">approximate</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
-    <span class="k">if</span> <span class="n">np</span><span class="o">.</span><span class="n">allclose</span><span class="p">(</span><span class="n">phi</span><span class="p">[</span><span class="mi">0</span><span class="p">],</span> <span class="mi">0</span><span class="p">)</span> <span class="ow">or</span> <span class="n">np</span><span class="o">.</span><span class="n">allclose</span><span class="p">(</span><span class="n">phi</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="mi">0</span><span class="p">):</span>
-        <span class="k">raise</span> <span class="ne">ValueError</span><span class="p">(</span><span class="s1">&#39;Sample Entropy is not defined.&#39;</span><span class="p">)</span>
+    <span class="n">x</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">asarray</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">dtype</span><span class="o">=</span><span class="n">np</span><span class="o">.</span><span class="n">float64</span><span class="p">)</span>
+    <span class="k">if</span> <span class="n">metric</span> <span class="o">==</span> <span class="s1">&#39;chebyshev&#39;</span> <span class="ow">and</span> <span class="n">x</span><span class="o">.</span><span class="n">size</span> <span class="o">&lt;</span> <span class="mi">5000</span><span class="p">:</span>
+        <span class="k">return</span> <span class="n">_numba_sampen</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">mm</span><span class="o">=</span><span class="n">order</span><span class="p">,</span> <span class="n">r</span><span class="o">=</span><span class="mf">0.2</span><span class="p">)</span>
     <span class="k">else</span><span class="p">:</span>
+        <span class="n">phi</span> <span class="o">=</span> <span class="n">_app_samp_entropy</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">order</span><span class="o">=</span><span class="n">order</span><span class="p">,</span> <span class="n">metric</span><span class="o">=</span><span class="n">metric</span><span class="p">,</span>
+                                <span class="n">approximate</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>
         <span class="k">return</span> <span class="o">-</span><span class="n">np</span><span class="o">.</span><span class="n">log</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">divide</span><span class="p">(</span><span class="n">phi</span><span class="p">[</span><span class="mi">1</span><span class="p">],</span> <span class="n">phi</span><span class="p">[</span><span class="mi">0</span><span class="p">]))</span></div>
 </pre></div>
 

diff --git a/docs/build/html/_sources/index.rst.txt b/docs/build/html/_sources/index.rst.txt
@@ -65,7 +65,7 @@ Entropy
     0.9945519071575192
     0.8482185855709181
     2.0754913760787277
-    2.1919237573930315
+    2.192416747827227
 
 Fractal dimension
 -----------------
@@ -123,7 +123,7 @@ Here are some benchmarks computed on an average PC (i7-7700HQ CPU @ 2.80 Ghz - 8
     150 µs ± 859 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
     42.4 µs ± 306 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
     4.59 ms ± 62.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
-    4.61 ms ± 163 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
+    2.03 ms ± 39.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
     16.4 µs ± 251 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
     32.4 µs ± 578 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
     17.4 µs ± 274 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

diff --git a/docs/build/html/generated/entropy.sample_entropy.html b/docs/build/html/generated/entropy.sample_entropy.html
@@ -119,7 +119,6 @@ <h1>entropy.sample_entropy</h1>
 </tbody>
 </table>
 <p class="rubric">Notes</p>
-<p>Original code from the mne-features package.</p>
 <p>Sample entropy is a modification of approximate entropy, used for assessing
 the complexity of physiological time-series signals. It has two advantages
 over approximate entropy: data length independence and a relatively
@@ -134,8 +133,11 @@ <h1>entropy.sample_entropy</h1>
 <span class="math notranslate nohighlight">\(m + 1\)</span> having a Chebyshev distance inferior to <span class="math notranslate nohighlight">\(r\)</span> and
 <span class="math notranslate nohighlight">\(C(m, r)\)</span> is the number of embedded vectors of length
 <span class="math notranslate nohighlight">\(m\)</span> having a Chebyshev distance inferior to <span class="math notranslate nohighlight">\(r\)</span>.</p>
-<p>Code adapted from the mne-features package by Jean-Baptiste Schiratti
-and Alexandre Gramfort.</p>
+<p>Note that if metric == ‘chebyshev’ and x.size &lt; 5000 points, then the
+sample entropy is computed using a fast custom Numba script. For other
+metric types or longer time-series, the sample entropy is computed using
+a code from the mne-features package by Jean-Baptiste Schiratti
+and Alexandre Gramfort (requires sklearn).</p>
 <p class="rubric">References</p>
 <table class="docutils citation" frame="void" id="ra05fe78bc240-1" rules="none">
 <colgroup><col class="label" /><col /></colgroup>