report: further edits of methodology throughput sections

Change-Id: I571d1a47743eb31ee10caf3f3336ac7437daf878
Signed-off-by: Maciek Konstantynowicz <mkonstan@cisco.com>

5 files changed, 905 insertions, 0 deletions

diff --git a/docs/report/introduction/methodology_data_plane_throughput/index.rst b/docs/report/introduction/methodology_data_plane_throughput/index.rst
new file mode 100644
index 0000000000..f87283810f
--- /dev/null
+++ b/docs/report/introduction/methodology_data_plane_throughput/index.rst
@@ -0,0 +1,10 @@
+Data Plane Throughput
+=====================
+
+.. toctree::
+
+    methodology_data_plane_throughput
+    methodology_mlrsearch_tests
+    methodology_mrr_throughput
+    methodology_plrsearch
+
diff --git a/docs/report/introduction/methodology_data_plane_throughput/methodology_data_plane_throughput.rst b/docs/report/introduction/methodology_data_plane_throughput/methodology_data_plane_throughput.rst
new file mode 100644
index 0000000000..762a7c2675
--- /dev/null
+++ b/docs/report/introduction/methodology_data_plane_throughput/methodology_data_plane_throughput.rst
@@ -0,0 +1,137 @@
+Data Plane Throughput Tests
+---------------------------
+
+Network data plane throughput is measured using multiple test methods in
+order to obtain representative and repeatable results across the large
+set of performance test cases implemented and executed within CSIT.
+
+Following throughput test methods are used:
+
+- MLRsearch - Multiple Loss Ratio search
+- MRR - Maximum Receive Rate
+- PLRsearch - Probabilistic Loss Ratio search
+
+Description of each test method is followed by generic test properties
+shared by all methods.
+
+MLRsearch Tests
+^^^^^^^^^^^^^^^
+
+Description
+~~~~~~~~~~~
+
+Multiple Loss Ratio search (MLRsearch) tests discover multiple packet
+throughput rates in a single search, reducing the overall test execution
+time compared to a binary search. Each rate associated with a distinct
+Packet Loss Ratio (PLR) criteria. In FD.io CSIT two throughput rates are
+discovered: Non-Drop Rate (NDR, with zero packet loss, PLR=0) and
+Partial Drop Rate (PDR, with PLR<0.5%). MLRsearch is compliant with
+:rfc:`2544`.
+
+Usage
+~~~~~
+
+MLRsearch tests are run to discover NDR and PDR rates for each VPP and
+DPDK release covered by CSIT report. Results for small frame sizes
+(64b/78B, IMIX) are presented in packet throughput graphs
+(Box-and-Whisker Plots) with NDR and PDR rates plotted against the test
+cases covering popular VPP packet paths.
+
+Each test is executed at least 10 times to verify measurements
+repeatability and results are compared between releases and test
+environments. NDR and PDR packet and bandwidth throughput results for
+all frame sizes and for all tests are presented in detailed results
+tables.
+
+Details
+~~~~~~~
+
+See :ref:`mlrsearch_algorithm` section for more detail. MLRsearch is
+being standardized in IETF in `draft-vpolak-mkonstan-mlrsearch
+<https://tools.ietf.org/html/draft-vpolak-mkonstan-bmwg-mlrsearch>`_.
+
+MRR Tests
+^^^^^^^^^
+
+Description
+~~~~~~~~~~~
+
+Maximum Receive Rate (MRR) tests are complementary to MLRsearch tests,
+as they provide a maximum “raw” throughput benchmark for development and
+testing community.
+
+MRR tests measure the packet forwarding rate under the maximum load
+offered by traffic generator (dependent on link type and NIC model) over
+a set trial duration, regardless of packet loss. Maximum load for
+specified Ethernet frame size is set to the bi-directional link rate.
+
+MRR tests are much faster than MLRsearch as they rely on a single trial
+or a small set of trials with very short duration. It is this property
+that makes them suitable for continuous execution in daily performance
+trending jobs enabling detection of performance anomalies (regressions,
+progressions) resulting from data plane code changes.
+
+MRR tests are also used for VPP per patch performance jobs verifying
+patch performance vs. parent. CSIT reports include MRR throughput
+comparisons between releases and test environments. Small frame sizes
+only (64b/78B, IMIX).
+
+Details
+~~~~~~~
+
+See :ref:`mrr_throughput` section for more detail about MRR tests
+configuration.
+
+FD.io CSIT performance dashboard includes complete description of
+`daily performance trending tests
+<https://docs.fd.io/csit/master/trending/methodology/performance_tests.html>`_
+and `VPP per patch tests
+<https://docs.fd.io/csit/master/trending/methodology/perpatch_performance_tests.html>`_.
+
+PLRsearch Tests
+^^^^^^^^^^^^^^^
+
+Description
+~~~~~~~~~~~
+
+Probabilistic Loss Ratio search (PLRsearch) tests discovers a packet
+throughput rate associated with configured Packet Loss Ratio (PLR)
+criteria for tests run over an extended period of time a.k.a. soak
+testing. PLRsearch assumes that system under test is probabilistic in
+nature, and not deterministic.
+
+Usage
+~~~~~
+
+PLRsearch are run to discover a sustained throughput for PLR=10^-7
+(close to NDR) for VPP release covered by CSIT report. Results for small
+frame sizes (64b/78B) are presented in packet throughput graphs (Box
+Plots) for a small subset of baseline tests.
+
+Each soak test lasts 2hrs and is executed at least twice. Results are
+compared against NDR and PDR rates discovered with MLRsearch.
+
+Details
+~~~~~~~
+
+See :ref:`plrsearch_algorithm` section for more detail. PLRsearch is
+being standardized in IETF in `draft-vpolak-bmwg-plrsearch
+<https://tools.ietf.org/html/draft-vpolak-bmwg-plrsearch>`_.
+
+Generic Test Properties
+^^^^^^^^^^^^^^^^^^^^^^^
+
+All data plane throughput test methodologies share following generic
+properties:
+
+- Tested L2 frame sizes (untagged Ethernet):
+
+  - IPv4 payload: 64B, IMIX (28x64B, 16x570B, 4x1518B), 1518B, 9000B.
+  - IPv6 payload: 78B, IMIX (28x78B, 16x570B, 4x1518B), 1518B, 9000B.
+  - All quoted sizes include frame CRC, but exclude per frame
+    transmission overhead of 20B (preamble, inter frame gap).
+
+- Offered packet load is always bi-directional and symmetric.
+- All measured and reported packet and bandwidth rates are aggregate
+  bi-directional rates reported from external Traffic Generator
+  perspective.
diff --git a/docs/report/introduction/methodology_data_plane_throughput/methodology_mlrsearch_tests.rst b/docs/report/introduction/methodology_data_plane_throughput/methodology_mlrsearch_tests.rst
new file mode 100644
index 0000000000..6c2adfcd06
--- /dev/null
+++ b/docs/report/introduction/methodology_data_plane_throughput/methodology_mlrsearch_tests.rst
@@ -0,0 +1,291 @@
+.. _mlrsearch_algorithm:
+
+MLRsearch Tests
+---------------
+
+Multiple Loss Rate search (MLRsearch) tests use new search algorithm
+implemented in FD.io CSIT project. MLRsearch discovers multiple packet
+throughput rates in a single search, with each rate associated with a
+distinct Packet Loss Ratio (PLR) criteria. MLRsearch is being
+standardized in IETF with `draft-vpolak-mkonstan-mlrsearch-XX
+<https://tools.ietf.org/html/draft-vpolak-mkonstan-mlrsearch-00>`_.
+
+Two throughput measurements used in FD.io CSIT are Non-Drop Rate (NDR,
+with zero packet loss, PLR=0) and Partial Drop Rate (PDR, with packet
+loss rate not greater than the configured non-zero PLR). MLRsearch
+discovers NDR and PDR in a single pass reducing required execution time
+compared to separate binary searches for NDR and PDR. MLRsearch reduces
+execution time even further by relying on shorter trial durations
+of intermediate steps, with only the final measurements
+conducted at the specified final trial duration.
+This results in the shorter overall search
+execution time when compared to a standard NDR/PDR binary search,
+while guaranteeing the same or similar results.
+
+If needed, MLRsearch can be easily adopted to discover more throughput rates
+with different pre-defined PLRs.
+
+.. Note:: All throughput rates are *always* bi-directional
+   aggregates of two equal (symmetric) uni-directional packet rates
+   received and reported by an external traffic generator.
+
+Overview
+~~~~~~~~
+
+The main properties of MLRsearch:
+
+- MLRsearch is a duration aware multi-phase multi-rate search algorithm.
+
+  - Initial phase determines promising starting interval for the search.
+  - Intermediate phases progress towards defined final search criteria.
+  - Final phase executes measurements according to the final search
+    criteria.
+
+- *Initial phase*:
+
+  - Uses link rate as a starting transmit rate and discovers the Maximum
+    Receive Rate (MRR) used as an input to the first intermediate phase.
+
+- *Intermediate phases*:
+
+  - Start with initial trial duration (in the first phase) and converge
+    geometrically towards the final trial duration (in the final phase).
+  - Track two values for NDR and two for PDR.
+
+    - The values are called (NDR or PDR) lower_bound and upper_bound.
+    - Each value comes from a specific trial measurement
+      (most recent for that transmit rate),
+      and as such the value is associated with that measurement's duration and
+      loss.
+    - A bound can be invalid, for example if NDR lower_bound
+      has been measured with nonzero loss.
+    - Invalid bounds are not real boundaries for the searched value,
+      but are needed to track interval widths.
+    - Valid bounds are real boundaries for the searched value.
+    - Each non-initial phase ends with all bounds valid.
+
+  - Start with a large (lower_bound, upper_bound) interval width and
+    geometrically converge towards the width goal (measurement resolution)
+    of the phase. Each phase halves the previous width goal.
+  - Use internal and external searches:
+
+    - External search - measures at transmit rates outside the (lower_bound,
+      upper_bound) interval. Activated when a bound is invalid,
+      to search for a new valid bound by doubling the interval width.
+      It is a variant of `exponential search`_.
+    - Internal search - `binary search`_, measures at transmit rates within the
+      (lower_bound, upper_bound) valid interval, halving the interval width.
+
+- *Final phase* is executed with the final test trial duration, and the final
+  width goal that determines resolution of the overall search.
+  Intermediate phases together with the final phase are called non-initial
+  phases.
+
+The main benefits of MLRsearch vs. binary search include:
+
+- In general MLRsearch is likely to execute more search trials overall, but
+  less trials at a set final duration.
+- In well behaving cases it greatly reduces (>50%) the overall duration
+  compared to a single PDR (or NDR) binary search duration,
+  while finding multiple drop rates.
+- In all cases MLRsearch yields the same or similar results to binary search.
+- Note: both binary search and MLRsearch are susceptible to reporting
+  non-repeatable results across multiple runs for very bad behaving
+  cases.
+
+Caveats:
+
+- Worst case MLRsearch can take longer than a binary search e.g. in case of
+  drastic changes in behaviour for trials at varying durations.
+
+Search Implementation
+~~~~~~~~~~~~~~~~~~~~~
+
+Following is a brief description of the current MLRsearch
+implementation in FD.io CSIT.
+
+Input Parameters
+````````````````
+
+#. *maximum_transmit_rate* - maximum packet transmit rate to be used by
+   external traffic generator, limited by either the actual Ethernet
+   link rate or traffic generator NIC model capabilities. Sample
+   defaults: 2 * 14.88 Mpps for 64B 10GE link rate,
+   2 * 18.75 Mpps for 64B 40GE NIC maximum rate.
+#. *minimum_transmit_rate* - minimum packet transmit rate to be used for
+   measurements. MLRsearch fails if lower transmit rate needs to be
+   used to meet search criteria. Default: 2 * 10 kpps (could be higher).
+#. *final_trial_duration* - required trial duration for final rate
+   measurements. Default: 30 sec.
+#. *initial_trial_duration* - trial duration for initial MLRsearch phase.
+   Default: 1 sec.
+#. *final_relative_width* - required measurement resolution expressed as
+   (lower_bound, upper_bound) interval width relative to upper_bound.
+   Default: 0.5%.
+#. *packet_loss_ratio* - maximum acceptable PLR search criteria for
+   PDR measurements. Default: 0.5%.
+#. *number_of_intermediate_phases* - number of phases between the initial
+   phase and the final phase. Impacts the overall MLRsearch duration.
+   Less phases are required for well behaving cases, more phases
+   may be needed to reduce the overall search duration for worse behaving
+   cases.
+   Default (2). (Value chosen based on limited experimentation to date.
+   More experimentation needed to arrive to clearer guidelines.)
+
+Initial Phase
+`````````````
+
+1. First trial measures at maximum rate and discovers MRR.
+
+   a. *in*: trial_duration = initial_trial_duration.
+   b. *in*: offered_transmit_rate = maximum_transmit_rate.
+   c. *do*: single trial.
+   d. *out*: measured loss ratio.
+   e. *out*: mrr = measured receive rate.
+
+2. Second trial measures at MRR and discovers MRR2.
+
+   a. *in*: trial_duration = initial_trial_duration.
+   b. *in*: offered_transmit_rate = MRR.
+   c. *do*: single trial.
+   d. *out*: measured loss ratio.
+   e. *out*: mrr2 = measured receive rate.
+
+3. Third trial measures at MRR2.
+
+   a. *in*: trial_duration = initial_trial_duration.
+   b. *in*: offered_transmit_rate = MRR2.
+   c. *do*: single trial.
+   d. *out*: measured loss ratio.
+
+Non-initial Phases
+``````````````````
+
+1. Main loop:
+
+   a. *in*: trial_duration for the current phase.
+      Set to initial_trial_duration for the first intermediate phase;
+      to final_trial_duration for the final phase;
+      or to the element of interpolating geometric sequence
+      for other intermediate phases.
+      For example with two intermediate phases, trial_duration
+      of the second intermediate phase is the geometric average
+      of initial_strial_duration and final_trial_duration.
+   b. *in*: relative_width_goal for the current phase.
+      Set to final_relative_width for the final phase;
+      doubled for each preceding phase.
+      For example with two intermediate phases,
+      the first intermediate phase uses quadruple of final_relative_width
+      and the second intermediate phase uses double of final_relative_width.
+   c. *in*: ndr_interval, pdr_interval from the previous main loop iteration
+      or the previous phase.
+      If the previous phase is the initial phase, both intervals have
+      lower_bound = MRR2, uper_bound = MRR.
+      Note that the initial phase is likely to create intervals with invalid
+      bounds.
+   d. *do*: According to the procedure described in point 2,
+      either exit the phase (by jumping to 1.g.),
+      or prepare new transmit rate to measure with.
+   e. *do*: Perform the trial measurement at the new transmit rate
+      and trial_duration, compute its loss ratio.
+   f. *do*: Update the bounds of both intervals, based on the new measurement.
+      The actual update rules are numerous, as NDR external search
+      can affect PDR interval and vice versa, but the result
+      agrees with rules of both internal and external search.
+      For example, any new measurement below an invalid lower_bound
+      becomes the new lower_bound, while the old measurement
+      (previously acting as the invalid lower_bound)
+      becomes a new and valid upper_bound.
+      Go to next iteration (1.c.), taking the updated intervals as new input.
+   g. *out*: current ndr_interval and pdr_interval.
+      In the final phase this is also considered
+      to be the result of the whole search.
+      For other phases, the next phase loop is started
+      with the current results as an input.
+
+2. New transmit rate (or exit) calculation (for 1.d.):
+
+   - If there is an invalid bound then prepare for external search:
+
+     - *If* the most recent measurement at NDR lower_bound transmit rate
+       had the loss higher than zero, then
+       the new transmit rate is NDR lower_bound
+       decreased by two NDR interval widths.
+     - Else, *if* the most recent measurement at PDR lower_bound
+       transmit rate had the loss higher than PLR, then
+       the new transmit rate is PDR lower_bound
+       decreased by two PDR interval widths.
+     - Else, *if* the most recent measurement at NDR upper_bound
+       transmit rate had no loss, then
+       the new transmit rate is NDR upper_bound
+       increased by two NDR interval widths.
+     - Else, *if* the most recent measurement at PDR upper_bound
+       transmit rate had the loss lower or equal to PLR, then
+       the new transmit rate is PDR upper_bound
+       increased by two PDR interval widths.
+   - If interval width is higher than the current phase goal:
+
+     - Else, *if* NDR interval does not meet the current phase width goal,
+       prepare for internal search. The new transmit rate is
+       (NDR lower bound + NDR upper bound) / 2.
+     - Else, *if* PDR interval does not meet the current phase width goal,
+       prepare for internal search. The new transmit rate is
+       (PDR lower bound + PDR upper bound) / 2.
+   - Else, *if* some bound has still only been measured at a lower duration,
+     prepare to re-measure at the current duration (and the same transmit
+     rate). The order of priorities is:
+
+     - NDR lower_bound,
+     - PDR lower_bound,
+     - NDR upper_bound,
+     - PDR upper_bound.
+   - *Else*, do not prepare any new rate, to exit the phase.
+     This ensures that at the end of each non-initial phase
+     all intervals are valid, narrow enough, and measured
+     at current phase trial duration.
+
+Implementation Deviations
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+This document so far has been describing a simplified version of MLRsearch
+algorithm. The full algorithm as implemented contains additional logic,
+which makes some of the details (but not general ideas) above incorrect.
+Here is a short description of the additional logic as a list of principles,
+explaining their main differences from (or additions to) the simplified
+description,but without detailing their mutual interaction.
+
+1. *Logarithmic transmit rate.*
+   In order to better fit the relative width goal,
+   the interval doubling and halving is done differently.
+   For example, the middle of 2 and 8 is 4, not 5.
+2. *Optimistic maximum rate.*
+   The increased rate is never higher than the maximum rate.
+   Upper bound at that rate is always considered valid.
+3. *Pessimistic minimum rate.*
+   The decreased rate is never lower than the minimum rate.
+   If a lower bound at that rate is invalid,
+   a phase stops refining the interval further (until it gets re-measured).
+4. *Conservative interval updates.*
+   Measurements above current upper bound never update a valid upper bound,
+   even if drop ratio is low.
+   Measurements below current lower bound always update any lower bound
+   if drop ratio is high.
+5. *Ensure sufficient interval width.*
+   Narrow intervals make external search take more time to find a valid bound.
+   If the new transmit increased or decreased rate would result in width
+   less than the current goal, increase/decrease more.
+   This can happen if the measurement for the other interval
+   makes the current interval too narrow.
+   Similarly, take care the measurements in the initial phase
+   create wide enough interval.
+6. *Timeout for bad cases.*
+   The worst case for MLRsearch is when each phase converges to intervals
+   way different than the results of the previous phase.
+   Rather than suffer total search time several times larger
+   than pure binary search, the implemented tests fail themselves
+   when the search takes too long (given by argument *timeout*).
+
+.. _binary search: https://en.wikipedia.org/wiki/Binary_search
+.. _exponential search: https://en.wikipedia.org/wiki/Exponential_search
+.. _estimation of standard deviation: https://en.wikipedia.org/wiki/Unbiased_estimation_of_standard_deviation
+.. _simplified error propagation formula: https://en.wikipedia.org/wiki/Propagation_of_uncertainty#Simplification
diff --git a/docs/report/introduction/methodology_data_plane_throughput/methodology_mrr_throughput.rst b/docs/report/introduction/methodology_data_plane_throughput/methodology_mrr_throughput.rst
new file mode 100644
index 0000000000..fd4baca2f3
--- /dev/null
+++ b/docs/report/introduction/methodology_data_plane_throughput/methodology_mrr_throughput.rst
@@ -0,0 +1,52 @@
+.. _mrr_throughput:
+
+MRR Throughput
+--------------
+
+Maximum Receive Rate (MRR) tests are complementary to MLRsearch tests,
+as they provide a maximum "raw" throughput benchmark for development and
+testing community. MRR tests measure the packet forwarding rate under
+the maximum load offered by traffic generator over a set trial duration,
+regardless of packet loss.
+
+MRR tests are currently used for following test jobs:
+
+- Report performance comparison: 64B, IMIX for vhost, memif.
+- Daily performance trending: 64B, IMIX for vhost, memif.
+- Per-patch performance verification: 64B.
+- PLRsearch soaking tests: 64B.
+
+Maximum offered load for specific L2 Ethernet frame size is set to
+either the maximum bi-directional link rate or tested NIC model
+capacity, as follows:
+
+- For 10GE NICs the maximum packet rate load is 2x14.88 Mpps for 64B, a
+  10GE bi-directional link rate.
+- For 25GE NICs the maximum packet rate load is 2x18.75 Mpps for 64B, a
+  25GE bi-directional link sub-rate limited by 25GE NIC used on TRex TG,
+  XXV710.
+- For 40GE NICs the maximum packet rate load is 2x18.75 Mpps for 64B, a
+  40GE bi-directional link sub-rate limited by 40GE NIC used on TRex
+  TG,XL710. Packet rate for other tested frame sizes is limited by
+  PCIeGen3 x8 bandwidth limitation of ~50Gbps.
+
+MRR test code implements multiple bursts of offered packet load and has
+two configurable burst parameters: individual trial duration and number
+of trials in a single burst. This enables more precise performance
+trending by providing more results data for analysis.
+
+Burst parameter settings vary between different tests using MRR:
+
+- MRR individual trial duration:
+
+  - Report performance comparison: 1 sec.
+  - Daily performance trending: 1 sec.
+  - Per-patch performance verification: 10 sec.
+  - PLRsearch soaking tests: 5.2 sec.
+
+- Number of MRR trials per burst:
+
+  - Report performance comparison: 10.
+  - Daily performance trending: 10.
+  - Per-patch performance verification: 5.
+  - PLRsearch soaking tests: 1.
\ No newline at end of file
diff --git a/docs/report/introduction/methodology_data_plane_throughput/methodology_plrsearch.rst b/docs/report/introduction/methodology_data_plane_throughput/methodology_plrsearch.rst
new file mode 100644
index 0000000000..fdd03472c3
--- /dev/null
+++ b/docs/report/introduction/methodology_data_plane_throughput/methodology_plrsearch.rst
@@ -0,0 +1,415 @@
+.. _plrsearch_algorithm:
+
+PLRsearch
+^^^^^^^^^
+
+Abstract algorithm
+~~~~~~~~~~~~~~~~~~
+
+.. TODO: Refer to packet forwarding terminology, such as "offered load" and
+   "loss ratio".
+
+Eventually, a better description of the abstract search algorithm
+will appear at this IETF standard: `plrsearch draft`_.
+
+Motivation
+``````````
+
+Network providers are interested in throughput a device can sustain.
+
+`RFC 2544`_ assumes loss ratio is given by a deterministic function of
+offered load. But NFV software devices are not deterministic (enough).
+This leads for deterministic algorithms (such as MLRsearch with single
+trial) to return results, which when repeated show relatively high
+standard deviation, thus making it harder to tell what "the throughput"
+actually is.
+
+We need another algorithm, which takes this indeterminism into account.
+
+Model
+`````
+
+Each algorithm searches for an answer to a precisely formulated
+question. When the question involves indeterministic systems, it has to
+specify probabilities (or prior distributions) which are tied to a
+specific probabilistic model. Different models will have different
+number (and meaning) of parameters. Complicated (but more realistic)
+models have many parameters, and the math involved can be very
+convoluted. It is better to start with simpler probabilistic model, and
+only change it when the output of the simpler algorithm is not stable or
+useful enough.
+
+This document is focused on algorithms related to packet loss count
+only. No latency (or other information) is taken into account. For
+simplicity, only one type of measurement is considered: dynamically
+computed offered load, constant within trial measurement of
+predetermined trial duration.
+
+The main idea of the search apgorithm is to iterate trial measurements,
+using `Bayesian inference`_ to compute both the current estimate
+of "the throughput" and the next offered load to measure at.
+The computations are done in parallel with the trial measurements.
+
+The following algorithm makes an assumption that packet traffic
+generator detects duplicate packets on receive detection, and reports
+this as an error.
+
+Poisson distribution
+````````````````````
+
+For given offered load, number of packets lost during trial measurement
+is assumed to come from `Poisson distribution`_,
+each trial is assumed to be independent, and the (unknown) Poisson parameter
+(average number of packets lost per second) is assumed to be
+constant across trials.
+
+When comparing different offered loads, the average loss per second is
+assumed to increase, but the (deterministic) function from offered load
+into average loss rate is otherwise unknown. This is called "loss function".
+
+Given a target loss ratio (configurable), there is an unknown offered load
+when the average is exactly that. We call that the "critical load".
+If critical load seems higher than maximum offerable load, we should use
+the maximum offerable load to make search output more conservative.
+
+Side note: `Binomial distribution`_ is a better fit compared to Poisson
+distribution (acknowledging that the number of packets lost cannot be
+higher than the number of packets offered), but the difference tends to
+be relevant in loads far above the critical region, so using Poisson
+distribution helps the algorithm focus on critical region better.
+
+Of course, there are great many increasing functions (as candidates
+for loss function). The offered load has to be chosen for each trial,
+and the computed posterior distribution of critical load
+changes with each trial result.
+
+To make the space of possible functions more tractable, some other
+simplifying assumptions are needed. As the algorithm will be examining
+(also) loads close to the critical load, linear approximation to the
+loss function in the critical region is important.
+But as the search algorithm needs to evaluate the function also far
+away from the critical region, the approximate function has to be
+well-behaved for every positive offered load, specifically it cannot predict
+non-positive packet loss rate.
+
+Within this document, "fitting function" is the name for such a well-behaved
+function, which approximates the unknown loss function in the critical region.
+
+Results from trials far from the critical region are likely to affect
+the critical rate estimate negatively, as the fitting function does not
+need to be a good approximation there. Discarding some results,
+or "suppressing" their impact with ad-hoc methods (other than
+using Poisson distribution instead of binomial) is not used, as such
+methods tend to make the overall search unstable. We rely on most of
+measurements being done (eventually) within the critical region, and
+overweighting far-off measurements (eventually) for well-behaved fitting
+functions.
+
+Speaking about new trials, each next trial will be done at offered load
+equal to the current average of the critical load.
+Alternative methods for selecting offered load might be used,
+in an attempt to speed up convergence, but such methods tend to be
+scpecific for a particular system under tests.
+
+Fitting function coefficients distribution
+``````````````````````````````````````````
+
+To accomodate systems with different behaviours, the fitting function is
+expected to have few numeric parameters affecting its shape (mainly
+affecting the linear approximation in the critical region).
+
+The general search algorithm can use whatever increasing fitting
+function, some specific functions can described later.
+
+It is up to implementer to chose a fitting function and prior
+distribution of its parameters. The rest of this document assumes each
+parameter is independently and uniformly distributed over a common
+interval. Implementers are to add non-linear transformations into their
+fitting functions if their prior is different.
+
+Exit condition for the search is either critical load stdev
+becoming small enough, or overal search time becoming long enough.
+
+The algorithm should report both avg and stdev for critical load. If the
+reported averages follow a trend (without reaching equilibrium), avg and
+stdev should refer to the equilibrium estimates based on the trend, not
+to immediate posterior values.
+
+Integration
+```````````
+
+The posterior distributions for fitting function parameters will not be
+integrable in general.
+
+The search algorithm utilises the fact that trial measurement takes some
+time, so this time can be used for numeric integration (using suitable
+method, such as Monte Carlo) to achieve sufficient precision.
+
+Optimizations
+`````````````
+
+After enough trials, the posterior distribution will be concentrated in
+a narrow area of parameter space. The integration method should take
+advantage of that.
+
+Even in the concentrated area, the likelihood can be quite small, so the
+integration algorithm should track the logarithm of the likelihood, and
+also avoid underflow errors by other means.
+
+FD.io CSIT Implementation Specifics
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The search receives min_rate and max_rate values, to avoid measurements
+at offered loads not supporeted by the traffic generator.
+
+The implemented tests cases use bidirectional traffic.
+The algorithm stores each rate as bidirectional rate (internally,
+the algorithm is agnostic to flows and directions,
+it only cares about overall counts of packets sent and packets lost),
+but debug output from traffic generator lists unidirectional values.
+
+Measurement delay
+`````````````````
+
+In a sample implemenation in FD.io CSIT project, there is roughly 0.5
+second delay between trials due to restrictons imposed by packet traffic
+generator in use (T-Rex).
+
+As measurements results come in, posterior distribution computation takes
+more time (per sample), although there is a considerable constant part
+(mostly for inverting the fitting functions).
+
+Also, the integrator needs a fair amount of samples to reach the region
+the posterior distribution is concentrated at.
+
+And of course, speed of the integrator depends on computing power
+of the CPU the algorithm is able to use.
+
+All those timing related effects are addressed by arithmetically increasing
+trial durations with configurable coefficients
+(currently 10.2 seconds for the first trial,
+each subsequent trial being 0.2 second longer).
+
+Rounding errors and underflows
+``````````````````````````````
+
+In order to avoid them, the current implementation tracks natural logarithm
+(instead of the original quantity) for any quantity which is never negative.
+Logarithm of zero is minus infinity (not supported by Python),
+so special value "None" is used instead.
+Specific functions for frequent operations
+(such as "logarithm of sum of exponentials")
+are defined to handle None correctly.
+
+Fitting functions
+`````````````````
+
+Current implementation uses two fitting functions.
+In general, their estimates for critical rate differ,
+which adds a simple source of systematic error,
+on top of randomness error reported by integrator.
+Otherwise the reported stdev of critical rate estimate
+is unrealistically low.
+
+Both functions are not only increasing, but convex
+(meaning the rate of increase is also increasing).
+
+As `primitive function`_ to any positive function is an increasing function,
+and primitive function to any increasing function is convex function;
+both fitting functions were constructed as double primitive function
+to a positive function (even though the intermediate increasing function
+is easier to describe).
+
+As not any function is integrable, some more realistic functions
+(especially with respect to behavior at very small offered loads)
+are not easily available.
+
+Both fitting functions have a "central point" and a "spread",
+varied by simply shifting and scaling (in x-axis, the offered load direction)
+the function to be doubly integrated.
+Scaling in y-axis (the loss rate direction) is fixed by the requirement of
+transfer rate staying nearly constant in very high offered loads.
+
+In both fitting functions (as they are a double primitive function
+to a symmetric function), the "central point" turns out
+to be equal to the aforementioned limiting transfer rate,
+so the fitting function parameter is named "mrr",
+the same quantity our Maximum Receive Rate tests are designed to measure.
+
+Both fitting functions return logarithm of loss rate,
+to avoid rounding errors and underflows.
+Parameters and offered load are not given as logarithms,
+as they are not expected to be extreme,
+and the formulas are simpler that way.
+
+Both fitting functions have several mathematically equivalent formulas,
+each can lead to an overflow or underflow in different places.
+Overflows can be eliminated by using different exact formulas
+for different argument ranges.
+Underflows can be avoided by using approximate formulas
+in affected argument ranges, such ranges have their own formulas to compute.
+At the end, both fitting function implementations
+contain multiple "if" branches, discontinuities are a possibility
+at range boundaries.
+
+Offered load for next trial measurement is the average of critical rate estimate.
+During each measurement, two estimates are computed,
+even though only one (in alternating order) is used for next offered load.
+
+Stretch function
+----------------
+
+The original function (before applying logarithm) is primitive function
+to `logistic function`_.
+The name "stretch" is used for related function
+in context of neural networks with sigmoid activation function.
+
+Erf function
+------------
+
+The original function is double primitive function to `Gaussian function`_.
+The name "erf" comes from error function, the first primitive to Gaussian.
+
+Prior distributions
+```````````````````
+
+The numeric integrator expects all the parameters to be distributed
+(independently and) uniformly on an interval (-1, 1).
+
+As both "mrr" and "spread" parameters are positive and not not dimensionless,
+a transformation is needed. Dimentionality is inherited from max_rate value.
+
+The "mrr" parameter follows a `Lomax distribution`_
+with alpha equal to one, but shifted so that mrr is always greater than 1
+packet per second.
+
+The "stretch" parameter is generated simply as the "mrr" value
+raised to a random power between zero and one;
+thus it follows a `reciprocal distribution`_.
+
+Integrator
+``````````
+
+After few measurements, the posterior distribution of fitting function
+arguments gets quite concentrated into a small area.
+The integrator is using `Monte Carlo`_ with `importance sampling`_
+where the biased distribution is `bivariate Gaussian`_ distribution,
+with deliberately larger variance.
+If the generated sample falls outside (-1, 1) interval,
+another sample is generated.
+
+The center and the variance for the biased distribution has three sources.
+First is a prior information. After enough samples are generated,
+the biased distribution is constructed from a mixture of two sources.
+Top 12 most weight samples, and all samples (the mix ratio is computed
+from the relative weights of the two populations).
+When integration (run along a particular measurement) is finished,
+the mixture bias distribution is used as the prior information
+for the next integration.
+
+This combination showed the best behavior, as the integrator usually follows
+two phases. First phase (where the top 12 samples are dominating)
+is mainly important for locating the new area the posterior distribution
+is concentrated at. The second phase (dominated by whole sample population)
+is actually relevant for the critical rate estimation.
+
+Caveats
+```````
+
+Current implementation does not constrict the critical rate
+(as computed for every sample) to the min_rate, max_rate interval.
+
+Earlier implementations were targeting loss rate (as opposed to loss ratio).
+The chosen fitting functions do allow arbitrarily low loss ratios,
+but may suffer from rounding errors in corresponding parameter regions.
+Internal loss rate target is computed from given loss ratio
+using the current trial offered load, which increases search instability,
+especially if measurements with surprisingly high loss count appear.
+
+As high loss count measurements add many bits of information,
+they need a large amount of small loss count measurements to balance them,
+making the algorithm converge quite slowly.
+
+Some systems evidently do not follow the assumption of repeated measurements
+having the same average loss rate (when offered load is the same).
+The idea of estimating the trend is not implemented at all,
+as the observed trends have varied characteristics.
+
+Probably, using a more realistic fitting functions
+will give better estimates than trend analysis.
+
+Graphical examples
+``````````````````
+
+The following pictures show the upper and lower bound (one sigma)
+on estimated critical rate, as computed by PLRsearch, after each trial measurement
+within the 30 minute duration of a test run.
+
+Both graphs are focusing on later estimates. Estimates computed from
+few initial measurements are wildly off the y-axis range shown.
+
+L2 patch
+--------
+
+This test case shows quite narrow critical region. Both fitting functions
+give similar estimates, the graph shows the randomness of measurements,
+and a trend. Both fitting functions seem to be somewhat overestimating
+the critical rate. The final estimated interval is too narrow,
+a longer run would report estimates somewhat bellow the current lower bound.
+
+.. only:: latex
+
+    .. raw:: latex
+
+        \begin{figure}[H]
+            \centering
+                \graphicspath{{../_tmp/src/introduction/}}
+                \includegraphics[width=0.90\textwidth]{PLR_patch}
+                \label{fig:PLR_patch}
+        \end{figure}
+
+.. only:: html
+
+    .. figure:: PLR_patch.svg
+        :alt: PLR_patch
+        :align: center
+
+Vhost
+-----
+
+This test case shows quite broad critical region. Fitting functions give
+fairly differing estimates. One overestimates, the other underestimates.
+The graph mostly shows later measurements slowly bringing the estimates
+towards each other. The final estimated interval is too broad,
+a longer run would return a smaller interval within the current one.
+
+.. only:: latex
+
+    .. raw:: latex
+
+        \begin{figure}[H]
+            \centering
+                \graphicspath{{../_tmp/src/introduction/}}
+                \includegraphics[width=0.90\textwidth]{PLR_vhost}
+                \label{fig:PLR_vhost}
+        \end{figure}
+
+.. only:: html
+
+    .. figure:: PLR_vhost.svg
+        :alt: PLR_vhost
+        :align: center
+
+.. _plrsearch draft: https://tools.ietf.org/html/draft-vpolak-bmwg-plrsearch-00
+.. _RFC 2544: https://tools.ietf.org/html/rfc2544
+.. _Bayesian inference: https://en.wikipedia.org/wiki/Bayesian_statistics
+.. _Poisson distribution: https://en.wikipedia.org/wiki/Poisson_distribution
+.. _Binomial distribution: https://en.wikipedia.org/wiki/Binomial_distribution
+.. _primitive function: https://en.wikipedia.org/wiki/Antiderivative
+.. _logistic function: https://en.wikipedia.org/wiki/Logistic_function
+.. _Gaussian function: https://en.wikipedia.org/wiki/Gaussian_function
+.. _Lomax distribution: https://en.wikipedia.org/wiki/Lomax_distribution
+.. _reciprocal distribution: https://en.wikipedia.org/wiki/Reciprocal_distribution
+.. _Monte Carlo: https://en.wikipedia.org/wiki/Monte_Carlo_integration
+.. _importance sampling: https://en.wikipedia.org/wiki/Importance_sampling
+.. _bivariate Gaussian: https://en.wikipedia.org/wiki/Multivariate_normal_distribution