.. _cbenchmarkbranchingrst:

=====================================
Measures branching in C++ from python
=====================================


.. only:: html

    **Links:** :download:`notebook <cbenchmark_branching.ipynb>`, :downloadlink:`html <cbenchmark_branching2html.html>`, :download:`PDF <cbenchmark_branching.pdf>`, :download:`python <cbenchmark_branching.py>`, :downloadlink:`slides <cbenchmark_branching.slides.html>`, :githublink:`GitHub|_doc/notebooks/cbenchmark_branching.ipynb|*`


This notebook looks into a couple of ways to write code, which one is
efficient, which one is not when it comes to write fast and short loops.
Both experiments are around
`branching <https://en.wikipedia.org/wiki/Branch_predictor>`__. The
notebook relies on C++ code implemented in
`cbenchmark.cpp <https://github.com/sdpython/cpyquickhelper/blob/master/cpyquickhelper/numbers/cbenchmark.cpp>`__
and
`repeat_fct.h <https://github.com/sdpython/cpyquickhelper/blob/master/cpyquickhelper/numbers/repeat_fct.h>`__.

.. code:: ipython3

    from jyquickhelper import add_notebook_menu
    add_notebook_menu()






.. contents::
    :local:





.. code:: ipython3

    %matplotlib inline

`numpy <http://www.numpy.org/>`__ is multithreaded. For an accurate
comparison, this needs to be disabled. This can be done as follows or by
setting environment variable ``MKL_NUM_THREADS=1``.

.. code:: ipython3

    try:
        import mkl
        mkl.set_num_threads(1)
    except ModuleNotFoundError as e:
        print('mkl not found', e)
        import os
        os.environ['MKL_NUM_THREADS']='1'

First experiment: comparison C++ syntax
---------------------------------------

This all started with article `Why is it faster to process a sorted
array than an unsorted
array? <https://stackoverflow.com/questions/11227809/why-is-it-faster-to-process-a-sorted-array-than-an-unsorted-array/11227902#11227902>`__.
It compares different implementation fo the following function for which
we try different implementations for the third line in next cell. The
last option is taken `Checking whether a number is positive or negative
using bitwise
operators <https://stackoverflow.com/questions/3779202/checking-whether-a-number-is-positive-or-negative-using-bitwise-operators>`__
which avoids
`branching <https://en.wikipedia.org/wiki/Branch_predictor>`__.

.. code:: ipython3

    # int nb = 0;
    # for(auto it = values.begin(); it != values.end(); ++it)
    #     if (*it >= th) nb++; // this line changes
    #     if (*it >= th) nb++; // and is repeated 10 times inside the loop.
    #     // ... 10 times
    # return nb;

The third line is also repeated 10 times to avoid the loop being too
significant.

.. code:: ipython3

    from cpyquickhelper.numbers.cbenchmark_dot import measure_scenario_A, measure_scenario_B
    from cpyquickhelper.numbers.cbenchmark_dot import measure_scenario_C, measure_scenario_D
    from cpyquickhelper.numbers.cbenchmark_dot import measure_scenario_E, measure_scenario_F
    from cpyquickhelper.numbers.cbenchmark_dot import measure_scenario_G, measure_scenario_H
    from cpyquickhelper.numbers.cbenchmark_dot import measure_scenario_I, measure_scenario_J

.. code:: ipython3

    import pandas
    
    def test_benchmark(label, values, th, repeat=10, number=20):
        funcs = [(k, v) for k, v in globals().copy().items() if k.startswith("measure_scenario")]
        rows = []
        for k, v in funcs:
            exe = v(values, th, repeat, number)
            d = exe.todict()
            d['doc'] = v.__doc__.split('``')[1]
            d['label'] = label
            d['name'] = k
            rows.append(d)        
        df = pandas.DataFrame(rows)
        return df
    
    test_benchmark("sorted", list(range(10)), 5)






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>doc</th>
          <th>label</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>name</th>
          <th>number</th>
          <th>repeat</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>6.518500e-08</td>
          <td>2.870060e-07</td>
          <td>if (values[i] &gt;= th) ++nb;</td>
          <td>sorted</td>
          <td>1.581000e-06</td>
          <td>0.000001</td>
          <td>measure_scenario_A</td>
          <td>20.0</td>
          <td>10.0</td>
        </tr>
        <tr>
          <th>1</th>
          <td>5.926000e-08</td>
          <td>2.583084e-07</td>
          <td>if (*it &gt;= th) ++nb;</td>
          <td>sorted</td>
          <td>1.186000e-06</td>
          <td>0.000001</td>
          <td>measure_scenario_B</td>
          <td>20.0</td>
          <td>10.0</td>
        </tr>
        <tr>
          <th>2</th>
          <td>6.123500e-08</td>
          <td>2.682277e-07</td>
          <td>if (*it &gt;= th) nb++;</td>
          <td>sorted</td>
          <td>1.580000e-06</td>
          <td>0.000001</td>
          <td>measure_scenario_C</td>
          <td>20.0</td>
          <td>10.0</td>
        </tr>
        <tr>
          <th>3</th>
          <td>1.185000e-08</td>
          <td>6.738195e-08</td>
          <td>nb += *it &gt;= th ? 1 : 0;</td>
          <td>sorted</td>
          <td>3.950000e-07</td>
          <td>0.000000</td>
          <td>measure_scenario_D</td>
          <td>20.0</td>
          <td>10.0</td>
        </tr>
        <tr>
          <th>4</th>
          <td>6.123500e-08</td>
          <td>2.682277e-07</td>
          <td>if (*it &gt;= th) nb += 1;</td>
          <td>sorted</td>
          <td>1.580000e-06</td>
          <td>0.000001</td>
          <td>measure_scenario_E</td>
          <td>20.0</td>
          <td>10.0</td>
        </tr>
        <tr>
          <th>5</th>
          <td>7.900000e-09</td>
          <td>5.530000e-08</td>
          <td>nb += (*it - th) &gt;= 0 ? 1 : 0;</td>
          <td>sorted</td>
          <td>3.950000e-07</td>
          <td>0.000000</td>
          <td>measure_scenario_F</td>
          <td>20.0</td>
          <td>10.0</td>
        </tr>
        <tr>
          <th>6</th>
          <td>1.580500e-08</td>
          <td>7.742840e-08</td>
          <td>nb += (*it - th) &lt; 0 ? 1 : 0;</td>
          <td>sorted</td>
          <td>3.960000e-07</td>
          <td>0.000000</td>
          <td>measure_scenario_G</td>
          <td>20.0</td>
          <td>10.0</td>
        </tr>
        <tr>
          <th>7</th>
          <td>1.383000e-08</td>
          <td>7.261929e-08</td>
          <td>nb += *it &lt; th ? 1 : 0;</td>
          <td>sorted</td>
          <td>3.960000e-07</td>
          <td>0.000000</td>
          <td>measure_scenario_H</td>
          <td>20.0</td>
          <td>10.0</td>
        </tr>
        <tr>
          <th>8</th>
          <td>1.777500e-08</td>
          <td>8.188513e-08</td>
          <td>nb += 1 ^ ((unsigned int)(*it) &gt;&gt; (sizeof(int)...</td>
          <td>sorted</td>
          <td>3.950000e-07</td>
          <td>0.000000</td>
          <td>measure_scenario_I</td>
          <td>20.0</td>
          <td>10.0</td>
        </tr>
        <tr>
          <th>9</th>
          <td>1.383000e-08</td>
          <td>7.261929e-08</td>
          <td>nb += values[i] &gt;= th ? 1 : 0;</td>
          <td>sorted</td>
          <td>3.960000e-07</td>
          <td>0.000000</td>
          <td>measure_scenario_J</td>
          <td>20.0</td>
          <td>10.0</td>
        </tr>
      </tbody>
    </table>
    </div>



Times are not very conclusive on such small lists.

.. code:: ipython3

    values = list(range(100000))
    df_sorted = test_benchmark("sorted", values, len(values)//2, repeat=200)
    df_sorted






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>doc</th>
          <th>label</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>name</th>
          <th>number</th>
          <th>repeat</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.001116</td>
          <td>0.004947</td>
          <td>if (values[i] &gt;= th) ++nb;</td>
          <td>sorted</td>
          <td>0.047184</td>
          <td>0.014096</td>
          <td>measure_scenario_A</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>1</th>
          <td>0.001136</td>
          <td>0.005057</td>
          <td>if (*it &gt;= th) ++nb;</td>
          <td>sorted</td>
          <td>0.038252</td>
          <td>0.014688</td>
          <td>measure_scenario_B</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>2</th>
          <td>0.001133</td>
          <td>0.005100</td>
          <td>if (*it &gt;= th) nb++;</td>
          <td>sorted</td>
          <td>0.044802</td>
          <td>0.013577</td>
          <td>measure_scenario_C</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0.000118</td>
          <td>0.000528</td>
          <td>nb += *it &gt;= th ? 1 : 0;</td>
          <td>sorted</td>
          <td>0.004973</td>
          <td>0.001436</td>
          <td>measure_scenario_D</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.000831</td>
          <td>0.003742</td>
          <td>if (*it &gt;= th) nb += 1;</td>
          <td>sorted</td>
          <td>0.048095</td>
          <td>0.013466</td>
          <td>measure_scenario_E</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>5</th>
          <td>0.000118</td>
          <td>0.000564</td>
          <td>nb += (*it - th) &gt;= 0 ? 1 : 0;</td>
          <td>sorted</td>
          <td>0.013677</td>
          <td>0.001338</td>
          <td>measure_scenario_F</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>6</th>
          <td>0.000170</td>
          <td>0.000786</td>
          <td>nb += (*it - th) &lt; 0 ? 1 : 0;</td>
          <td>sorted</td>
          <td>0.007175</td>
          <td>0.001582</td>
          <td>measure_scenario_G</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>7</th>
          <td>0.000090</td>
          <td>0.000414</td>
          <td>nb += *it &lt; th ? 1 : 0;</td>
          <td>sorted</td>
          <td>0.004574</td>
          <td>0.001296</td>
          <td>measure_scenario_H</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>8</th>
          <td>0.000119</td>
          <td>0.000518</td>
          <td>nb += 1 ^ ((unsigned int)(*it) &gt;&gt; (sizeof(int)...</td>
          <td>sorted</td>
          <td>0.002834</td>
          <td>0.002243</td>
          <td>measure_scenario_I</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>9</th>
          <td>0.000081</td>
          <td>0.000360</td>
          <td>nb += values[i] &gt;= th ? 1 : 0;</td>
          <td>sorted</td>
          <td>0.004324</td>
          <td>0.001459</td>
          <td>measure_scenario_J</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
      </tbody>
    </table>
    </div>



The article some implementations will be slower if the values are not
sorted.

.. code:: ipython3

    import random
    random.shuffle(values)
    values = values.copy()
    values[:10]




.. parsed-literal::
    [73680, 56372, 84502, 10263, 93712, 12350, 98785, 54243, 7299, 30309]



.. code:: ipython3

    df_shuffled = test_benchmark("shuffled", values, len(values)//2, repeat=200)
    df_shuffled






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>doc</th>
          <th>label</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>name</th>
          <th>number</th>
          <th>repeat</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.000704</td>
          <td>0.003076</td>
          <td>if (values[i] &gt;= th) ++nb;</td>
          <td>shuffled</td>
          <td>0.019643</td>
          <td>0.013091</td>
          <td>measure_scenario_A</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>1</th>
          <td>0.000697</td>
          <td>0.003042</td>
          <td>if (*it &gt;= th) ++nb;</td>
          <td>shuffled</td>
          <td>0.016656</td>
          <td>0.013203</td>
          <td>measure_scenario_B</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>2</th>
          <td>0.000761</td>
          <td>0.003440</td>
          <td>if (*it &gt;= th) nb++;</td>
          <td>shuffled</td>
          <td>0.027062</td>
          <td>0.013112</td>
          <td>measure_scenario_C</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0.000128</td>
          <td>0.000563</td>
          <td>nb += *it &gt;= th ? 1 : 0;</td>
          <td>shuffled</td>
          <td>0.004075</td>
          <td>0.001579</td>
          <td>measure_scenario_D</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.001281</td>
          <td>0.005590</td>
          <td>if (*it &gt;= th) nb += 1;</td>
          <td>shuffled</td>
          <td>0.028929</td>
          <td>0.018232</td>
          <td>measure_scenario_E</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>5</th>
          <td>0.000126</td>
          <td>0.000554</td>
          <td>nb += (*it - th) &gt;= 0 ? 1 : 0;</td>
          <td>shuffled</td>
          <td>0.003189</td>
          <td>0.001542</td>
          <td>measure_scenario_F</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>6</th>
          <td>0.000144</td>
          <td>0.000634</td>
          <td>nb += (*it - th) &lt; 0 ? 1 : 0;</td>
          <td>shuffled</td>
          <td>0.005403</td>
          <td>0.001844</td>
          <td>measure_scenario_G</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>7</th>
          <td>0.000128</td>
          <td>0.000560</td>
          <td>nb += *it &lt; th ? 1 : 0;</td>
          <td>shuffled</td>
          <td>0.003976</td>
          <td>0.001567</td>
          <td>measure_scenario_H</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>8</th>
          <td>0.000223</td>
          <td>0.000974</td>
          <td>nb += 1 ^ ((unsigned int)(*it) &gt;&gt; (sizeof(int)...</td>
          <td>shuffled</td>
          <td>0.007231</td>
          <td>0.002795</td>
          <td>measure_scenario_I</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
        <tr>
          <th>9</th>
          <td>0.000140</td>
          <td>0.000687</td>
          <td>nb += values[i] &gt;= th ? 1 : 0;</td>
          <td>shuffled</td>
          <td>0.016393</td>
          <td>0.001547</td>
          <td>measure_scenario_J</td>
          <td>20.0</td>
          <td>200.0</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    df = pandas.concat([df_sorted, df_shuffled])
    dfg = df[["doc", "label", "average"]].pivot("doc", "label", "average")

.. code:: ipython3

    ax = dfg.plot.bar(rot=30)
    labels = [l.get_text() for l in ax.get_xticklabels()]
    ax.set_xticklabels(labels, ha='right')
    ax.set_title("Comparison of all implementations");



.. image:: cbenchmark_branching_16_0.png


It seems that inline tests (``cond ? value1 : value2``) do not stop the
branching and it should be used whenever possible.

.. code:: ipython3

    sdf = df[["doc", "label", "average"]]
    dfg2 = sdf[sdf.doc.str.contains('[?^]')].pivot("doc", "label", "average")
    ax = dfg2.plot.bar(rot=30)
    labels = [l.get_text() for l in ax.get_xticklabels()]
    ax.set_xticklabels(labels, ha='right')
    ax.set_title("Comparison of implementations using ? :");



.. image:: cbenchmark_branching_18_0.png


.. code:: ipython3

    sdf = df[["doc", "label", "average"]]
    dfg2 = sdf[sdf.doc.str.contains('if')].pivot("doc", "label", "average")
    ax = dfg2.plot.bar(rot=30)
    labels = [l.get_text() for l in ax.get_xticklabels()]
    ax.set_xticklabels(labels, ha='right')
    ax.set_ylim([0.0004, 0.0020])
    ax.set_title("Comparison of implementations using tests");



.. image:: cbenchmark_branching_19_0.png


*sorted*, *not sorted* does not seem to have a real impact in this case.
It shows *branching* really slows down the execution of a program.
Branching happens whenever the program meets a loop condition or a test.
Iterator ``*it`` are faster than accessing an array with notation
``[i]`` which adds a cost due to an extra addition.

Second experiment: dot product
------------------------------

The goal is to compare the dot product from
`numpy.dot <https://docs.scipy.org/doc/numpy-1.15.0/reference/generated/numpy.dot.html>`__
and a couple of implementation in C++ which look like this:

.. code:: ipython3

    # float vector_dot_product_pointer(const float *p1, const float *p2, size_t size)
    # {
    #     float sum = 0;
    #     const float * end1 = p1 + size;
    #     for(; p1 != end1; ++p1, ++p2)
    #         sum += *p1 * *p2;
    #     return sum;
    # }
    # 
    # 
    # float vector_dot_product(py::array_t<float> v1, py::array_t<float> v2)
    # {
    #     if (v1.ndim() != v2.ndim())
    #         throw std::runtime_error("Vector v1 and v2 must have the same dimension.");
    #     if (v1.ndim() != 1)
    #         throw std::runtime_error("Vector v1 and v2 must be vectors.");
    #     return vector_dot_product_pointer(v1.data(0), v2.data(0), v1.shape(0));
    # }

numpy vs C++
~~~~~~~~~~~~

numpy.dot
^^^^^^^^^

.. code:: ipython3

    %matplotlib inline

.. code:: ipython3

    import numpy
    
    def simple_dot(values):
        return numpy.dot(values, values)
    
    values = list(range(10000000))
    values = numpy.array(values, dtype=numpy.float32)
    vect = values / numpy.max(values)
    simple_dot(vect)




.. parsed-literal::
    3333333.2



.. code:: ipython3

    vect.dtype




.. parsed-literal::
    dtype('float32')



.. code:: ipython3

    from timeit import Timer
    
    def measure_time(stmt, context, repeat=10, number=50):
        tim = Timer(stmt, globals=context)
        res = numpy.array(tim.repeat(repeat=repeat, number=number))
        mean = numpy.mean(res)
        dev = numpy.mean(res ** 2)
        dev = (dev - mean**2) ** 0.5
        return dict(average=mean, deviation=dev, min_exec=numpy.min(res),
                    max_exec=numpy.max(res), repeat=repeat, number=number,
                    size=context['values'].shape[0])
    
    measure_time("simple_dot(values)", context=dict(simple_dot=simple_dot, values=vect))




.. parsed-literal::
    {'average': 0.17493714570000024,
     'deviation': 0.01055679328637683,
     'min_exec': 0.16323108099999928,
     'max_exec': 0.1946009330000038,
     'repeat': 10,
     'number': 50,
     'size': 10000000}



.. code:: ipython3

    res = []
    for i in range(10, 200000, 2500):
        t = measure_time("simple_dot(values)", repeat=100,
                         context=dict(simple_dot=simple_dot, values=vect[:i].copy()))
        res.append(t)
    
    import pandas
    dot = pandas.DataFrame(res)
    dot.tail()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>number</th>
          <th>repeat</th>
          <th>size</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>75</th>
          <td>0.001811</td>
          <td>0.000932</td>
          <td>0.009380</td>
          <td>0.001230</td>
          <td>50</td>
          <td>100</td>
          <td>187510</td>
        </tr>
        <tr>
          <th>76</th>
          <td>0.001646</td>
          <td>0.000415</td>
          <td>0.003583</td>
          <td>0.001249</td>
          <td>50</td>
          <td>100</td>
          <td>190010</td>
        </tr>
        <tr>
          <th>77</th>
          <td>0.001981</td>
          <td>0.000676</td>
          <td>0.004695</td>
          <td>0.001267</td>
          <td>50</td>
          <td>100</td>
          <td>192510</td>
        </tr>
        <tr>
          <th>78</th>
          <td>0.002436</td>
          <td>0.000855</td>
          <td>0.004890</td>
          <td>0.001313</td>
          <td>50</td>
          <td>100</td>
          <td>195010</td>
        </tr>
        <tr>
          <th>79</th>
          <td>0.001886</td>
          <td>0.000353</td>
          <td>0.002882</td>
          <td>0.001327</td>
          <td>50</td>
          <td>100</td>
          <td>197510</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    res = []
    for i in range(100000, 10000000, 1000000):
        t = measure_time("simple_dot(values)", repeat=10,
                         context=dict(simple_dot=simple_dot, values=vect[:i].copy()))
        res.append(t)
        
    huge_dot = pandas.DataFrame(res)
    huge_dot.head()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>number</th>
          <th>repeat</th>
          <th>size</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>0</th>
          <td>0.000914</td>
          <td>0.000144</td>
          <td>0.001219</td>
          <td>0.000820</td>
          <td>50</td>
          <td>10</td>
          <td>100000</td>
        </tr>
        <tr>
          <th>1</th>
          <td>0.015862</td>
          <td>0.004123</td>
          <td>0.027136</td>
          <td>0.012628</td>
          <td>50</td>
          <td>10</td>
          <td>1100000</td>
        </tr>
        <tr>
          <th>2</th>
          <td>0.035084</td>
          <td>0.003601</td>
          <td>0.042250</td>
          <td>0.030178</td>
          <td>50</td>
          <td>10</td>
          <td>2100000</td>
        </tr>
        <tr>
          <th>3</th>
          <td>0.048146</td>
          <td>0.001690</td>
          <td>0.050411</td>
          <td>0.045873</td>
          <td>50</td>
          <td>10</td>
          <td>3100000</td>
        </tr>
        <tr>
          <th>4</th>
          <td>0.070299</td>
          <td>0.005897</td>
          <td>0.086010</td>
          <td>0.066763</td>
          <td>50</td>
          <td>10</td>
          <td>4100000</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    import matplotlib.pyplot as plt
    fig, ax = plt.subplots(1, 2, figsize=(14,4))
    dot.plot(x='size', y="average", ax=ax[0])
    huge_dot.plot(x='size', y="average", ax=ax[1], logy=True)
    ax[0].set_title("numpy dot product execution time");
    ax[1].set_title("numpy dot product execution time");



.. image:: cbenchmark_branching_31_0.png


numpy.einsum
^^^^^^^^^^^^

.. code:: ipython3

    def simple_dot_einsum(values):
        return numpy.einsum('i,i->', values, values)
    
    values = list(range(10000000))
    values = numpy.array(values, dtype=numpy.float32)
    vect = values / numpy.max(values)
    simple_dot_einsum(vect)




.. parsed-literal::
    3333333.8



.. code:: ipython3

    measure_time("simple_dot_einsum(values)",
                 context=dict(simple_dot_einsum=simple_dot_einsum, values=vect))




.. parsed-literal::
    {'average': 0.2862720898999996,
     'deviation': 0.019172035812954895,
     'min_exec': 0.26404798899999093,
     'max_exec': 0.3281182589999929,
     'repeat': 10,
     'number': 50,
     'size': 10000000}



.. code:: ipython3

    res = []
    for i in range(10, 200000, 2500):
        t = measure_time("simple_dot_einsum(values)", repeat=100,
                         context=dict(simple_dot_einsum=simple_dot_einsum, values=vect[:i].copy()))
        res.append(t)
    
    import pandas
    einsum_dot = pandas.DataFrame(res)
    einsum_dot.tail()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>number</th>
          <th>repeat</th>
          <th>size</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>75</th>
          <td>0.004400</td>
          <td>0.001858</td>
          <td>0.009241</td>
          <td>0.003290</td>
          <td>50</td>
          <td>100</td>
          <td>187510</td>
        </tr>
        <tr>
          <th>76</th>
          <td>0.003400</td>
          <td>0.000065</td>
          <td>0.003683</td>
          <td>0.003320</td>
          <td>50</td>
          <td>100</td>
          <td>190010</td>
        </tr>
        <tr>
          <th>77</th>
          <td>0.003717</td>
          <td>0.000703</td>
          <td>0.007098</td>
          <td>0.003467</td>
          <td>50</td>
          <td>100</td>
          <td>192510</td>
        </tr>
        <tr>
          <th>78</th>
          <td>0.003628</td>
          <td>0.000256</td>
          <td>0.004555</td>
          <td>0.003404</td>
          <td>50</td>
          <td>100</td>
          <td>195010</td>
        </tr>
        <tr>
          <th>79</th>
          <td>0.003660</td>
          <td>0.000224</td>
          <td>0.004962</td>
          <td>0.003448</td>
          <td>50</td>
          <td>100</td>
          <td>197510</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    import matplotlib.pyplot as plt
    fig, ax = plt.subplots(1, 1, figsize=(7,4))
    dot.plot(x='size', y="average", ax=ax, label="numpy.dot", logy=True)
    einsum_dot.plot(x='size', y="average", ax=ax, logy=True,label="numpy.einsum")
    ax.set_title("numpy einsum / dot dot product execution time");



.. image:: cbenchmark_branching_36_0.png


The function
`einsum <https://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html>`__
is slower (see `Einsum - Einstein summation in deep
learning <https://rockt.github.io/2018/04/30/einsum>`__ appears to be
slower but it is usually faster when it comes to chain operations as it
reduces the number of intermediate allocations to do.

pybind11
^^^^^^^^

Now the custom implementation. We start with an empty function to get a
sense of the cost due to to pybind11.

.. code:: ipython3

    from cpyquickhelper.numbers.cbenchmark_dot import empty_vector_dot_product
    empty_vector_dot_product(vect, vect)




.. parsed-literal::
    0.0



.. code:: ipython3

    def empty_c11_dot(vect):
        return empty_vector_dot_product(vect, vect)
    
    measure_time("empty_c11_dot(values)", 
                 context=dict(empty_c11_dot=empty_c11_dot, values=vect), repeat=10)




.. parsed-literal::
    {'average': 0.0001423012999993034,
     'deviation': 4.8738723576776925e-05,
     'min_exec': 0.00010429600000350092,
     'max_exec': 0.00025876499999810676,
     'repeat': 10,
     'number': 50,
     'size': 10000000}



Very small. It should not pollute our experiments.

.. code:: ipython3

    from cpyquickhelper.numbers.cbenchmark_dot import vector_dot_product
    vector_dot_product(vect, vect)




.. parsed-literal::
    3334629.0



.. code:: ipython3

    def c11_dot(vect):
        return vector_dot_product(vect, vect)
    
    measure_time("c11_dot(values)", 
                 context=dict(c11_dot=c11_dot, values=vect), repeat=10)




.. parsed-literal::
    {'average': 0.8533396122000013,
     'deviation': 0.24492425073976992,
     'min_exec': 0.6820375839999997,
     'max_exec': 1.2789236320000015,
     'repeat': 10,
     'number': 50,
     'size': 10000000}



.. code:: ipython3

    res = []
    for i in range(10, 200000, 2500):
        t = measure_time("c11_dot(values)", repeat=10,
                         context=dict(c11_dot=c11_dot, values=vect[:i].copy()))
        res.append(t)
    
    import pandas
    cus_dot = pandas.DataFrame(res)
    cus_dot.tail()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>number</th>
          <th>repeat</th>
          <th>size</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>75</th>
          <td>0.023365</td>
          <td>0.001148</td>
          <td>0.025472</td>
          <td>0.021720</td>
          <td>50</td>
          <td>10</td>
          <td>187510</td>
        </tr>
        <tr>
          <th>76</th>
          <td>0.024021</td>
          <td>0.000643</td>
          <td>0.025093</td>
          <td>0.022678</td>
          <td>50</td>
          <td>10</td>
          <td>190010</td>
        </tr>
        <tr>
          <th>77</th>
          <td>0.023822</td>
          <td>0.001232</td>
          <td>0.025876</td>
          <td>0.021724</td>
          <td>50</td>
          <td>10</td>
          <td>192510</td>
        </tr>
        <tr>
          <th>78</th>
          <td>0.023852</td>
          <td>0.000867</td>
          <td>0.025552</td>
          <td>0.022637</td>
          <td>50</td>
          <td>10</td>
          <td>195010</td>
        </tr>
        <tr>
          <th>79</th>
          <td>0.024870</td>
          <td>0.000563</td>
          <td>0.025481</td>
          <td>0.023793</td>
          <td>50</td>
          <td>10</td>
          <td>197510</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    fig, ax = plt.subplots(1, 2, figsize=(14,4))
    dot.plot(x='size', y="average", ax=ax[0], label="numpy")
    cus_dot.plot(x='size', y="average", ax=ax[0], label="pybind11")
    dot.plot(x='size', y="average", ax=ax[1], label="numpy", logy=True)
    cus_dot.plot(x='size', y="average", ax=ax[1], label="pybind11")
    ax[0].set_title("numpy and custom dot product execution time");
    ax[1].set_title("numpy and custom dot product execution time");



.. image:: cbenchmark_branching_46_0.png


Pretty slow. Let’s see what it does to compute dot product 16 by 16.

BLAS
~~~~

Internally, `numpy <http://www.numpy.org/>`__ is using
`BLAS <http://www.netlib.org/blas/>`__. A direct call to it should give
the same results.

.. code:: ipython3

    from cpyquickhelper.numbers.direct_blas_lapack import cblas_sdot
    
    def blas_dot(vect):
        return cblas_sdot(vect, vect)
    
    measure_time("blas_dot(values)", context=dict(blas_dot=blas_dot, values=vect), repeat=10)




.. parsed-literal::
    {'average': 0.18376368689999795,
     'deviation': 0.013271289356658413,
     'min_exec': 0.16785606499999517,
     'max_exec': 0.20954452600000195,
     'repeat': 10,
     'number': 50,
     'size': 10000000}



.. code:: ipython3

    res = []
    for i in range(10, 200000, 2500):
        t = measure_time("blas_dot(values)", repeat=10,
                         context=dict(blas_dot=blas_dot, values=vect[:i].copy()))
        res.append(t)
    
    import pandas
    blas_dot = pandas.DataFrame(res)
    blas_dot.tail()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>number</th>
          <th>repeat</th>
          <th>size</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>75</th>
          <td>0.002430</td>
          <td>0.000223</td>
          <td>0.002573</td>
          <td>0.001950</td>
          <td>50</td>
          <td>10</td>
          <td>187510</td>
        </tr>
        <tr>
          <th>76</th>
          <td>0.002326</td>
          <td>0.000258</td>
          <td>0.002600</td>
          <td>0.001931</td>
          <td>50</td>
          <td>10</td>
          <td>190010</td>
        </tr>
        <tr>
          <th>77</th>
          <td>0.002530</td>
          <td>0.000197</td>
          <td>0.002635</td>
          <td>0.001977</td>
          <td>50</td>
          <td>10</td>
          <td>192510</td>
        </tr>
        <tr>
          <th>78</th>
          <td>0.002071</td>
          <td>0.000440</td>
          <td>0.002698</td>
          <td>0.001519</td>
          <td>50</td>
          <td>10</td>
          <td>195010</td>
        </tr>
        <tr>
          <th>79</th>
          <td>0.002458</td>
          <td>0.000400</td>
          <td>0.003242</td>
          <td>0.001869</td>
          <td>50</td>
          <td>10</td>
          <td>197510</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    fig, ax = plt.subplots(1, 2, figsize=(14,4))
    dot.plot(x='size', y="average", ax=ax[0], label="numpy")
    cus_dot.plot(x='size', y="average", ax=ax[0], label="pybind11")
    blas_dot.plot(x='size', y="average", ax=ax[0], label="blas")
    dot.plot(x='size', y="average", ax=ax[1], label="numpy", logy=True)
    cus_dot.plot(x='size', y="average", ax=ax[1], label="pybind11")
    blas_dot.plot(x='size', y="average", ax=ax[1], label="blas")
    ax[0].set_title("numpy and custom dot product execution time");
    ax[1].set_title("numpy and custom dot product execution time");



.. image:: cbenchmark_branching_51_0.png


Use of branching: 16 multplications in one row
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The code looks like what follows. If there is more than 16
multiplications left, we use function *vector_dot_product_pointer16*,
otherwise, there are done one by one like the previous function.

.. code:: ipython3

    # float vector_dot_product_pointer16(const float *p1, const float *p2)
    # {
    #     float sum = 0;
    #     
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    # 
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    #     sum += *(p1++) * *(p2++);
    # 
    #     return sum;
    # }
    # 
    # #define BYN 16
    # 
    # float vector_dot_product_pointer16(const float *p1, const float *p2, size_t size)
    # {
    #     float sum = 0;
    #     size_t i = 0;
    #     if (size >= BYN) {
    #         size_t size_ = size - BYN;
    #         for(; i < size_; i += BYN, p1 += BYN, p2 += BYN)
    #             sum += vector_dot_product_pointer16(p1, p2);
    #     }
    #     for(; i < size; ++p1, ++p2, ++i)
    #         sum += *p1 * *p2;
    #     return sum;
    # }
    # 
    # float vector_dot_product16(py::array_t<float> v1, py::array_t<float> v2)
    # {
    #     if (v1.ndim() != v2.ndim())
    #         throw std::runtime_error("Vector v1 and v2 must have the same dimension.");
    #     if (v1.ndim() != 1)
    #         throw std::runtime_error("Vector v1 and v2 must be vectors.");
    #     return vector_dot_product_pointer16(v1.data(0), v2.data(0), v1.shape(0));
    # }

.. code:: ipython3

    from cpyquickhelper.numbers.cbenchmark_dot import vector_dot_product16
    vector_dot_product16(vect, vect)




.. parsed-literal::
    3333331.75



.. code:: ipython3

    def c11_dot16(vect):
        return vector_dot_product16(vect, vect)
    
    measure_time("c11_dot16(values)", 
                 context=dict(c11_dot16=c11_dot16, values=vect), repeat=10)




.. parsed-literal::
    {'average': 0.5628858367999967,
     'deviation': 0.006558496803961198,
     'min_exec': 0.5514294899999896,
     'max_exec': 0.5766948529999922,
     'repeat': 10,
     'number': 50,
     'size': 10000000}



.. code:: ipython3

    res = []
    for i in range(10, 200000, 2500):
        t = measure_time("c11_dot16(values)", repeat=10,
                         context=dict(c11_dot16=c11_dot16, values=vect[:i].copy()))
        res.append(t)
    
    cus_dot16 = pandas.DataFrame(res)
    cus_dot16.tail()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>number</th>
          <th>repeat</th>
          <th>size</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>75</th>
          <td>0.010757</td>
          <td>0.000345</td>
          <td>0.011495</td>
          <td>0.010465</td>
          <td>50</td>
          <td>10</td>
          <td>187510</td>
        </tr>
        <tr>
          <th>76</th>
          <td>0.010169</td>
          <td>0.000474</td>
          <td>0.011362</td>
          <td>0.009392</td>
          <td>50</td>
          <td>10</td>
          <td>190010</td>
        </tr>
        <tr>
          <th>77</th>
          <td>0.010440</td>
          <td>0.000449</td>
          <td>0.011532</td>
          <td>0.010153</td>
          <td>50</td>
          <td>10</td>
          <td>192510</td>
        </tr>
        <tr>
          <th>78</th>
          <td>0.011015</td>
          <td>0.000679</td>
          <td>0.011867</td>
          <td>0.009231</td>
          <td>50</td>
          <td>10</td>
          <td>195010</td>
        </tr>
        <tr>
          <th>79</th>
          <td>0.011091</td>
          <td>0.000540</td>
          <td>0.012168</td>
          <td>0.010415</td>
          <td>50</td>
          <td>10</td>
          <td>197510</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    fig, ax = plt.subplots(1, 2, figsize=(14,4))
    dot.plot(x='size', y="average", ax=ax[0], label="numpy")
    cus_dot.plot(x='size', y="average", ax=ax[0], label="pybind11")
    cus_dot16.plot(x='size', y="average", ax=ax[0], label="pybind11x16")
    dot.plot(x='size', y="average", ax=ax[1], label="numpy", logy=True)
    cus_dot.plot(x='size', y="average", ax=ax[1], label="pybind11")
    cus_dot16.plot(x='size', y="average", ax=ax[1], label="pybind11x16")
    ax[0].set_title("numpy and custom dot product execution time");
    ax[1].set_title("numpy and custom dot product execution time");



.. image:: cbenchmark_branching_57_0.png


We are far from *numpy* but the branching has clearly a huge impact and
the fact the loop condition is evaluated only every 16 iterations does
not explain this gain. Next experiment with
`SSE <https://fr.wikipedia.org/wiki/Streaming_SIMD_Extensions>`__
instructions.

Optimized to remove function call
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

We remove the function call to get the following version.

.. code:: ipython3

    # float vector_dot_product_pointer16_nofcall(const float *p1, const float *p2, size_t size)
    # {
    #     float sum = 0;    
    #     const float * end = p1 + size;
    #     if (size >= BYN) {
    #         #if(BYN != 16)
    #             #error "BYN must be equal to 16";
    #         #endif
    #         unsigned int size_ = (unsigned int) size;
    #         size_ = size_ >> 4;  // division by 16=2^4
    #         size_ = size_ << 4;  // multiplication by 16=2^4
    #         const float * end_ = p1 + size_;
    #         for(; p1 != end_;)
    #         {
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             sum += *p1 * *p2; ++p1, ++p2;
    # 
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             sum += *p1 * *p2; ++p1, ++p2;
    #             sum += *p1 * *p2; ++p1, ++p2;
    #         }
    #     }
    #     for(; p1 != end; ++p1, ++p2)
    #         sum += *p1 * *p2;
    #     return sum;
    # }
    # 
    # float vector_dot_product16_nofcall(py::array_t<float> v1, py::array_t<float> v2)
    # {
    #     if (v1.ndim() != v2.ndim())
    #         throw std::runtime_error("Vector v1 and v2 must have the same dimension.");
    #     if (v1.ndim() != 1)
    #         throw std::runtime_error("Vector v1 and v2 must be vectors.");
    #     return vector_dot_product_pointer16_nofcall(v1.data(0), v2.data(0), v1.shape(0));
    # }

.. code:: ipython3

    from cpyquickhelper.numbers.cbenchmark_dot import vector_dot_product16_nofcall
    vector_dot_product16_nofcall(vect, vect)




.. parsed-literal::
    3334629.0



.. code:: ipython3

    def c11_dot16_nofcall(vect):
        return vector_dot_product16_nofcall(vect, vect)
    
    measure_time("c11_dot16_nofcall(values)",
                 context=dict(c11_dot16_nofcall=c11_dot16_nofcall, values=vect), repeat=10)




.. parsed-literal::
    {'average': 1.2945859600000034,
     'deviation': 0.029453146426095998,
     'min_exec': 1.270528181000003,
     'max_exec': 1.3768783189999851,
     'repeat': 10,
     'number': 50,
     'size': 10000000}



.. code:: ipython3

    res = []
    for i in range(10, 200000, 2500):
        t = measure_time("c11_dot16_nofcall(values)", repeat=10,
                         context=dict(c11_dot16_nofcall=c11_dot16_nofcall, values=vect[:i].copy()))
        res.append(t)
    
    cus_dot16_nofcall = pandas.DataFrame(res)
    cus_dot16_nofcall.tail()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>number</th>
          <th>repeat</th>
          <th>size</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>75</th>
          <td>0.023715</td>
          <td>0.001017</td>
          <td>0.025011</td>
          <td>0.021424</td>
          <td>50</td>
          <td>10</td>
          <td>187510</td>
        </tr>
        <tr>
          <th>76</th>
          <td>0.024223</td>
          <td>0.000463</td>
          <td>0.024711</td>
          <td>0.023164</td>
          <td>50</td>
          <td>10</td>
          <td>190010</td>
        </tr>
        <tr>
          <th>77</th>
          <td>0.024831</td>
          <td>0.000362</td>
          <td>0.025402</td>
          <td>0.024273</td>
          <td>50</td>
          <td>10</td>
          <td>192510</td>
        </tr>
        <tr>
          <th>78</th>
          <td>0.025039</td>
          <td>0.000294</td>
          <td>0.025515</td>
          <td>0.024486</td>
          <td>50</td>
          <td>10</td>
          <td>195010</td>
        </tr>
        <tr>
          <th>79</th>
          <td>0.026592</td>
          <td>0.001819</td>
          <td>0.031072</td>
          <td>0.024430</td>
          <td>50</td>
          <td>10</td>
          <td>197510</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    fig, ax = plt.subplots(1, 2, figsize=(14,4))
    dot.plot(x='size', y="average", ax=ax[0], label="numpy")
    cus_dot.plot(x='size', y="average", ax=ax[0], label="pybind11")
    cus_dot16.plot(x='size', y="average", ax=ax[0], label="pybind11x16")
    cus_dot16_nofcall.plot(x='size', y="average", ax=ax[0], label="pybind11x16_nofcall")
    dot.plot(x='size', y="average", ax=ax[1], label="numpy", logy=True)
    cus_dot.plot(x='size', y="average", ax=ax[1], label="pybind11")
    cus_dot16.plot(x='size', y="average", ax=ax[1], label="pybind11x16")
    cus_dot16_nofcall.plot(x='size', y="average", ax=ax[1], label="pybind11x16_nofcall")
    ax[0].set_title("numpy and custom dot product execution time");
    ax[1].set_title("numpy and custom dot product execution time");



.. image:: cbenchmark_branching_64_0.png


Weird, branching did not happen when the code is not inside a separate
function.

SSE instructions
~~~~~~~~~~~~~~~~

We replace one function in the previous implementation.

.. code:: ipython3

    # #include <xmmintrin.h>
    # 
    # float vector_dot_product_pointer16_sse(const float *p1, const float *p2)
    # {
    #     __m128 c1 = _mm_load_ps(p1);
    #     __m128 c2 = _mm_load_ps(p2);
    #     __m128 r1 = _mm_mul_ps(c1, c2);
    #     
    #     p1 += 4;
    #     p2 += 4;
    #     
    #     c1 = _mm_load_ps(p1);
    #     c2 = _mm_load_ps(p2);
    #     r1 = _mm_add_ps(r1, _mm_mul_ps(c1, c2));
    #     
    #     p1 += 4;
    #     p2 += 4;
    #     
    #     c1 = _mm_load_ps(p1);
    #     c2 = _mm_load_ps(p2);
    #     r1 = _mm_add_ps(r1, _mm_mul_ps(c1, c2));
    #     
    #     p1 += 4;
    #     p2 += 4;
    #     
    #     c1 = _mm_load_ps(p1);
    #     c2 = _mm_load_ps(p2);
    #     r1 = _mm_add_ps(r1, _mm_mul_ps(c1, c2));
    # 
    #     float r[4];
    #     _mm_store_ps(r, r1);
    # 
    #     return r[0] + r[1] + r[2] + r[3];
    # }

.. code:: ipython3

    from cpyquickhelper.numbers.cbenchmark_dot import vector_dot_product16_sse
    vector_dot_product16_sse(vect, vect)




.. parsed-literal::
    3333332.0



.. code:: ipython3

    def c11_dot16_sse(vect):
        return vector_dot_product16_sse(vect, vect)
    
    measure_time("c11_dot16_sse(values)", 
                 context=dict(c11_dot16_sse=c11_dot16_sse, values=vect), repeat=10)




.. parsed-literal::
    {'average': 0.22386672670000393,
     'deviation': 0.004340312827217225,
     'min_exec': 0.21708269199999108,
     'max_exec': 0.23220287799998118,
     'repeat': 10,
     'number': 50,
     'size': 10000000}



.. code:: ipython3

    res = []
    for i in range(10, 200000, 2500):
        t = measure_time("c11_dot16_sse(values)", repeat=10,
                         context=dict(c11_dot16_sse=c11_dot16_sse, values=vect[:i].copy()))
        res.append(t)
    
    cus_dot16_sse = pandas.DataFrame(res)
    cus_dot16_sse.tail()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>number</th>
          <th>repeat</th>
          <th>size</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>75</th>
          <td>0.003629</td>
          <td>0.000200</td>
          <td>0.003825</td>
          <td>0.003051</td>
          <td>50</td>
          <td>10</td>
          <td>187510</td>
        </tr>
        <tr>
          <th>76</th>
          <td>0.003571</td>
          <td>0.000230</td>
          <td>0.003793</td>
          <td>0.003063</td>
          <td>50</td>
          <td>10</td>
          <td>190010</td>
        </tr>
        <tr>
          <th>77</th>
          <td>0.003642</td>
          <td>0.000267</td>
          <td>0.003821</td>
          <td>0.003110</td>
          <td>50</td>
          <td>10</td>
          <td>192510</td>
        </tr>
        <tr>
          <th>78</th>
          <td>0.003583</td>
          <td>0.000285</td>
          <td>0.003846</td>
          <td>0.003146</td>
          <td>50</td>
          <td>10</td>
          <td>195010</td>
        </tr>
        <tr>
          <th>79</th>
          <td>0.004467</td>
          <td>0.000770</td>
          <td>0.005784</td>
          <td>0.003754</td>
          <td>50</td>
          <td>10</td>
          <td>197510</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    fig, ax = plt.subplots(1, 2, figsize=(14,4))
    dot.plot(x='size', y="average", ax=ax[0], label="numpy")
    cus_dot16_sse.plot(x='size', y="average", ax=ax[0], label="pybind11x16_sse")
    dot.plot(x='size', y="average", ax=ax[1], label="numpy", logy=True)
    cus_dot16_sse.plot(x='size', y="average", ax=ax[1], label="pybind11x16_sse")
    cus_dot.plot(x='size', y="average", ax=ax[1], label="pybind11")
    cus_dot16.plot(x='size', y="average", ax=ax[1], label="pybind11x16")
    ax[0].set_title("numpy and custom dot product execution time");
    ax[1].set_title("numpy and custom dot product execution time");



.. image:: cbenchmark_branching_71_0.png


Better even though it is still slower than *numpy*. It is closer. Maybe
the compilation option are not optimized, *numpy* was also compiled with
the Intel compiler. To be accurate, multi-threading must be disabled on
*numpy* side. That’s the purpose of the first two lines.

AVX 512
~~~~~~~

Last experiment with `AVX
512 <https://fr.wikipedia.org/wiki/Advanced_Vector_Extensions>`__
instructions but it does not work on all processor. I could not test it
on my laptop as these instructions do not seem to be available. More can
be found on wikipedia `CPUs with
AVX-512 <https://en.wikipedia.org/wiki/AVX-512#CPUs_with_AVX-512>`__.

.. code:: ipython3

    import platform
    platform.processor()




.. parsed-literal::
    'Intel64 Family 6 Model 78 Stepping 3, GenuineIntel'



.. code:: ipython3

    import numpy
    values = numpy.array(list(range(10000000)), dtype=numpy.float32)
    vect = values / numpy.max(values)

.. code:: ipython3

    from cpyquickhelper.numbers.cbenchmark_dot import vector_dot_product16_avx512
    vector_dot_product16_avx512(vect, vect)




.. parsed-literal::
    3333332.0



.. code:: ipython3

    def c11_dot16_avx512(vect):
        return vector_dot_product16_avx512(vect, vect)
    
    measure_time("c11_dot16_avx512(values)",
                 context=dict(c11_dot16_avx512=c11_dot16_avx512, values=vect), repeat=10)




.. parsed-literal::
    {'average': 0.24621207209999624,
     'deviation': 0.030678971301695272,
     'min_exec': 0.21640002699999172,
     'max_exec': 0.3037627230000055,
     'repeat': 10,
     'number': 50,
     'size': 10000000}



.. code:: ipython3

    res = []
    for i in range(10, 200000, 2500):
        t = measure_time("c11_dot16_avx512(values)", repeat=10,
                         context=dict(c11_dot16_avx512=c11_dot16_avx512, values=vect[:i].copy()))
        res.append(t)
    
    cus_dot16_avx512 = pandas.DataFrame(res)
    cus_dot16_avx512.tail()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>number</th>
          <th>repeat</th>
          <th>size</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>75</th>
          <td>0.003631</td>
          <td>0.001094</td>
          <td>0.005245</td>
          <td>0.001953</td>
          <td>50</td>
          <td>10</td>
          <td>187510</td>
        </tr>
        <tr>
          <th>76</th>
          <td>0.004108</td>
          <td>0.001309</td>
          <td>0.006254</td>
          <td>0.002541</td>
          <td>50</td>
          <td>10</td>
          <td>190010</td>
        </tr>
        <tr>
          <th>77</th>
          <td>0.003714</td>
          <td>0.000933</td>
          <td>0.005624</td>
          <td>0.002511</td>
          <td>50</td>
          <td>10</td>
          <td>192510</td>
        </tr>
        <tr>
          <th>78</th>
          <td>0.005389</td>
          <td>0.001072</td>
          <td>0.007044</td>
          <td>0.003420</td>
          <td>50</td>
          <td>10</td>
          <td>195010</td>
        </tr>
        <tr>
          <th>79</th>
          <td>0.005102</td>
          <td>0.001096</td>
          <td>0.008251</td>
          <td>0.004557</td>
          <td>50</td>
          <td>10</td>
          <td>197510</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    fig, ax = plt.subplots(1, 2, figsize=(14,4))
    dot.plot(x='size', y="average", ax=ax[0], label="numpy")
    cus_dot16.plot(x='size', y="average", ax=ax[0], label="pybind11x16")
    cus_dot16_sse.plot(x='size', y="average", ax=ax[0], label="pybind11x16_sse")
    cus_dot16_avx512.plot(x='size', y="average", ax=ax[0], label="pybind11x16_avx512")
    dot.plot(x='size', y="average", ax=ax[1], label="numpy", logy=True)
    cus_dot16.plot(x='size', y="average", ax=ax[1], label="pybind11x16")
    cus_dot16_sse.plot(x='size', y="average", ax=ax[1], label="pybind11x16_sse")
    cus_dot16_avx512.plot(x='size', y="average", ax=ax[1], label="pybind11x16_avx512")
    cus_dot.plot(x='size', y="average", ax=ax[1], label="pybind11")
    ax[0].set_title("numpy and custom dot product execution time");
    ax[1].set_title("numpy and custom dot product execution time");



.. image:: cbenchmark_branching_79_0.png


If the time is the same, it means that options AVX512 are not available.

.. code:: ipython3

    from cpyquickhelper.numbers.cbenchmark import get_simd_available_option
    get_simd_available_option()




.. parsed-literal::
    'Available options:  __SSE__ __SSE2__ __SSE3__ __SSE4_1__'



Last call with OpenMP
---------------------

.. code:: ipython3

    from cpyquickhelper.numbers.cbenchmark_dot import vector_dot_product_openmp
    vector_dot_product_openmp(vect, vect, 2)




.. parsed-literal::
    3348343.25



.. code:: ipython3

    vector_dot_product_openmp(vect, vect, 4)




.. parsed-literal::
    3343080.75



.. code:: ipython3

    def c11_dot_openmp2(vect):
        return vector_dot_product_openmp(vect, vect, nthreads=2)
    
    def c11_dot_openmp4(vect):
        return vector_dot_product_openmp(vect, vect, nthreads=4)
    
    measure_time("c11_dot_openmp2(values)",
                 context=dict(c11_dot_openmp2=c11_dot_openmp2, values=vect), repeat=10)




.. parsed-literal::
    {'average': 0.6101010486999939,
     'deviation': 0.023996900833560427,
     'min_exec': 0.5678723409999975,
     'max_exec': 0.6511117819999299,
     'repeat': 10,
     'number': 50,
     'size': 10000000}



.. code:: ipython3

    measure_time("c11_dot_openmp4(values)",
                 context=dict(c11_dot_openmp4=c11_dot_openmp4, values=vect), repeat=10)




.. parsed-literal::
    {'average': 0.329937042999984,
     'deviation': 0.016537939361682446,
     'min_exec': 0.30559936400004517,
     'max_exec': 0.35731013800000255,
     'repeat': 10,
     'number': 50,
     'size': 10000000}



.. code:: ipython3

    res = []
    for i in range(10, 200000, 2500):
        t = measure_time("c11_dot_openmp2(values)", repeat=10,
                         context=dict(c11_dot_openmp2=c11_dot_openmp2, values=vect[:i].copy()))
        res.append(t)
    
    cus_dot_openmp2 = pandas.DataFrame(res)
    cus_dot_openmp2.tail()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>number</th>
          <th>repeat</th>
          <th>size</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>75</th>
          <td>0.011458</td>
          <td>0.000959</td>
          <td>0.012970</td>
          <td>0.010097</td>
          <td>50</td>
          <td>10</td>
          <td>187510</td>
        </tr>
        <tr>
          <th>76</th>
          <td>0.011582</td>
          <td>0.000723</td>
          <td>0.012886</td>
          <td>0.010232</td>
          <td>50</td>
          <td>10</td>
          <td>190010</td>
        </tr>
        <tr>
          <th>77</th>
          <td>0.012416</td>
          <td>0.000774</td>
          <td>0.013324</td>
          <td>0.011078</td>
          <td>50</td>
          <td>10</td>
          <td>192510</td>
        </tr>
        <tr>
          <th>78</th>
          <td>0.013305</td>
          <td>0.000465</td>
          <td>0.013685</td>
          <td>0.012125</td>
          <td>50</td>
          <td>10</td>
          <td>195010</td>
        </tr>
        <tr>
          <th>79</th>
          <td>0.009506</td>
          <td>0.001906</td>
          <td>0.012279</td>
          <td>0.007230</td>
          <td>50</td>
          <td>10</td>
          <td>197510</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    res = []
    for i in range(10, 200000, 2500):
        t = measure_time("c11_dot_openmp4(values)", repeat=10,
                         context=dict(c11_dot_openmp4=c11_dot_openmp4, values=vect[:i].copy()))
        res.append(t)
    
    cus_dot_openmp4 = pandas.DataFrame(res)
    cus_dot_openmp4.tail()






.. raw:: html

    <div>
    <style scoped>
        .dataframe tbody tr th:only-of-type {
            vertical-align: middle;
        }

        .dataframe tbody tr th {
            vertical-align: top;
        }

        .dataframe thead th {
            text-align: right;
        }
    </style>
    <table border="1" class="dataframe">
      <thead>
        <tr style="text-align: right;">
          <th></th>
          <th>average</th>
          <th>deviation</th>
          <th>max_exec</th>
          <th>min_exec</th>
          <th>number</th>
          <th>repeat</th>
          <th>size</th>
        </tr>
      </thead>
      <tbody>
        <tr>
          <th>75</th>
          <td>0.006757</td>
          <td>0.000376</td>
          <td>0.007589</td>
          <td>0.006341</td>
          <td>50</td>
          <td>10</td>
          <td>187510</td>
        </tr>
        <tr>
          <th>76</th>
          <td>0.007416</td>
          <td>0.001347</td>
          <td>0.011127</td>
          <td>0.006291</td>
          <td>50</td>
          <td>10</td>
          <td>190010</td>
        </tr>
        <tr>
          <th>77</th>
          <td>0.006946</td>
          <td>0.000783</td>
          <td>0.008544</td>
          <td>0.005562</td>
          <td>50</td>
          <td>10</td>
          <td>192510</td>
        </tr>
        <tr>
          <th>78</th>
          <td>0.006466</td>
          <td>0.000636</td>
          <td>0.007366</td>
          <td>0.005439</td>
          <td>50</td>
          <td>10</td>
          <td>195010</td>
        </tr>
        <tr>
          <th>79</th>
          <td>0.006986</td>
          <td>0.000389</td>
          <td>0.007906</td>
          <td>0.006533</td>
          <td>50</td>
          <td>10</td>
          <td>197510</td>
        </tr>
      </tbody>
    </table>
    </div>



.. code:: ipython3

    fig, ax = plt.subplots(1, 2, figsize=(14,4))
    dot.plot(x='size', y="average", ax=ax[0], label="numpy")
    cus_dot16.plot(x='size', y="average", ax=ax[0], label="pybind11x16")
    cus_dot16_sse.plot(x='size', y="average", ax=ax[0], label="pybind11x16_sse")
    cus_dot_openmp2.plot(x='size', y="average", ax=ax[0], label="cus_dot_openmp2")
    cus_dot_openmp4.plot(x='size', y="average", ax=ax[0], label="cus_dot_openmp4")
    dot.plot(x='size', y="average", ax=ax[1], label="numpy", logy=True)
    cus_dot16.plot(x='size', y="average", ax=ax[1], label="pybind11x16")
    cus_dot16_sse.plot(x='size', y="average", ax=ax[1], label="pybind11x16_sse")
    cus_dot_openmp2.plot(x='size', y="average", ax=ax[1], label="cus_dot_openmp2")
    cus_dot_openmp4.plot(x='size', y="average", ax=ax[1], label="cus_dot_openmp4")
    cus_dot.plot(x='size', y="average", ax=ax[1], label="pybind11")
    ax[0].set_title("numpy and custom dot product execution time");
    ax[1].set_title("numpy and custom dot product execution time");



.. image:: cbenchmark_branching_89_0.png


Parallelization does not solve everything, efficient is important.

Back to numpy
~~~~~~~~~~~~~

This article `Why is matrix multiplication faster with numpy than with
ctypes in
Python? <https://stackoverflow.com/questions/10442365/why-is-matrix-multiplication-faster-with-numpy-than-with-ctypes-in-python>`__
gives some kints on why *numpy* is still faster. By looking at the code
of the dot product in *numpy*:
`arraytypes.c.src <https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/arraytypes.c.src>`__,
it seems that *numpy* does a simple dot product without using branching
or uses the library `BLAS <http://www.netlib.org/blas/>`__ which is the
case in this benchmark (code for dot product:
`sdot.c <http://www.netlib.org/clapack/cblas/sdot.c>`__). And it does
use *branching*. See also function
`blas_stride <https://github.com/numpy/numpy/blob/master/numpy/core/src/multiarray/common.h#L298>`__.
These libraries then play with `compilation
options <https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html>`__
and optimize for speed. This benchmark does not look into
`cython-blis <https://github.com/explosion/cython-blis>`__ which
implements some `BLAS <http://www.netlib.org/blas/>`__ functions with an
assembly language and has different implementations depending on the
platform it is used. A little bit more on C++ optimization `How to
optimize C and C++ code in
2018 <https://medium.com/@aka.rider/how-to-optimize-c-and-c-code-in-2018-bd4f90a72c2b>`__.