--- /dev/null
+# Copyright (c) 2019 Cisco and/or its affiliates.
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at:
+#
+# http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+"""Module holding PLRsearch class."""
+
+import logging
+import math
+import multiprocessing
+import time
+
+import dill
+# TODO: Inform pylint about scipy (of correct version) being available.
+from scipy.special import erfcx, erfc
+
+# TODO: The preferred way to consume this code is via a pip package.
+# If your project copies code instead, make sure your pylint check does not
+# require these imports to be absolute and descending from your superpackage.
+import Integrator
+from log_plus import log_plus, log_minus
+
+
+class PLRsearch(object):
+ """A class to encapsulate data relevant for search method.
+
+ The context is performance testing of packet processing systems.
+ The system, when being offered a steady stream of packets,
+ can process some of them successfully, other are considered "lost".
+
+ See docstring of the search method for algorithm description.
+
+ Two constants are stored as class fields for speed.
+
+ Method othed than search (and than __init__)
+ are just internal code structure.
+ TODO: Those method names should start with underscore then.
+
+ TODO: Figure out how to replace #print with logging
+ without slowing down too much.
+ """
+
+ xerfcx_limit = math.pow(math.acos(0), -0.5)
+ log_xerfcx_10 = math.log(xerfcx_limit - math.exp(10) * erfcx(math.exp(10)))
+
+ def __init__(
+ self, measurer, trial_duration_per_trial, packet_loss_ratio_target,
+ trial_number_offset=0, timeout=60.0):
+ """Store rate measurer and additional parameters.
+
+ Also declare packet_loss_per_second_target field (float),
+ to be initialized later.
+
+ TODO: Copy AbstractMeasurer from MLRsearch.
+
+ :param measurer: The measurer to call when searching.
+ :param trial_duration_per_trial: Each trial has larger duration
+ than the previous trial. This is the increment, in seconds.
+ :param packet_loss_ratio_target: The algorithm tries to estimate
+ the offered load leading to this ratio on average.
+ Trial ratio is number of packets lost divided by packets offered.
+ :param trial_number_offset: The "first" trial number will be 1+this.
+ Use this to ensure first iterations have enough time to compute
+ reasonable estimates for later trials to use.
+ :param timeout: The search ends if it lasts more than this many seconds.
+ :type measurer: MLRsearch.AbstractMeasurer
+ :type trial_duration_per_trial: float
+ :type packet_loss_ratio_target: float
+ :type trial_number_offset: int
+ :type timeout: float
+ """
+ self.measurer = measurer
+ self.trial_duration_per_trial = trial_duration_per_trial
+ self.packet_loss_ratio_target = packet_loss_ratio_target
+ self.trial_number_offset = trial_number_offset
+ self.timeout = timeout
+
+ def search(self, min_rate, max_rate):
+ """Perform the search, return average and stdev for throughput estimate.
+
+ Considering measurer and packet_loss_ratio_target (see __init__),
+ find such an offered load (called critical load) that is expected
+ to hit the target loss ratio in the limit of very long trial duration.
+ As the system is probabilistic (and test duration is finite),
+ the critical ratio is only estimated.
+ Return the average and standard deviation of the estimate.
+
+ In principle, this algorithm performs trial measurements,
+ each with varied offered load (which is constant during the trial).
+ During each measurement, Bayesian inference is performed
+ on all the measurement results so far.
+ When timeout is up, the last estimate is returned,
+ else another trial is performed.
+
+ It is assumed that the system under test, even though not deterministic,
+ still follows the rule of large numbers. In another words,
+ any growing set of measurements at a particular offered load
+ will converge towards unique (for the given load) packet loss ratio.
+ This means there is a deterministic (but unknown) function
+ mapping the offered load to average loss ratio.
+ This function is called loss function.
+ This also assumes the average loss ratio
+ does not depend on trial duration.
+
+ The actual probability distribution of loss counts, achieving
+ the average ratio on trials of various duration
+ can be complicated (and can depend on offered load), but simply assuming
+ Poisson distribution will make the algorithm converge.
+ Binomial distribution would be more precise,
+ but Poisson is more practical, as it effectively gives
+ less information content to high ratio results.
+
+ Even when applying other assumptions on the loss function
+ (increasing function, limit zero ratio when load goes to zero,
+ global upper limit on rate of packets processed), there are still
+ too many different shapes of possible loss functions,
+ which makes full Bayesian reasoning intractable.
+
+ This implementation radically simplifies things by examining
+ only two shapes, each with finitely many (in this case just two)
+ parameters. In other words, two fitting functions
+ (each with two parameters and one argument).
+ When restricting model space to one of the two fitting functions,
+ the Bayesian inference becomes tractable (even though it needs
+ numerical integration from Integrator class).
+
+ The first measurement is done at min_rate to help with convergence
+ if max_rate measurements give loss below target.
+ FIXME: Fix overflow error and really use min_rate.
+ The second measurement is done at max_rate, next few measurements
+ have offered load of previous load minus excess loss ratio.
+ This simple rule is found to be good when offered loads
+ so far are way above the critical rate. After few measurements,
+ inference from fitting functions converges faster that the initial
+ "optimistic" procedure.
+
+ Offered loads close to (limiting) critical rate are the most useful,
+ as linear approximation of the fitting function
+ becomes good enough there (thus reducing the impact
+ of the overall shape of fitting function).
+ After several trials, usually one of the fitting functions
+ has better predictions than the other one, but the algorithm
+ does not track that. Simply, it uses the estimate average,
+ alternating between the functions.
+
+ The returned average and stdev is a combination of the two fitting
+ estimates.
+
+ TODO: If measurement at max_rate is already below target loss rate,
+ the algorithm always measures at max_rate, and likelihood
+ no longer has single maximum, leading to "random" estimates.
+ Find a way to avoid that.
+
+ :param min_rate: Avoid measuring at offered loads below this,
+ in packets per second.
+ :param max_rate: Avoid measuring at offered loads above this,
+ in packets per second.
+ :type min_rate: float
+ :type max_rate: float
+ :returns: Average and stdev of critical load estimate.
+ :rtype: 2-tuple of floats
+ """
+ stop_time = time.time() + self.timeout
+ min_rate = float(min_rate)
+ max_rate = float(max_rate)
+ trial_result_list = list()
+ trial_number = self.trial_number_offset
+ integrator_data = (None, None, None, None)
+ message = "Trial %(number)r computed avg %(avg)r stdev %(stdev)r"
+ message += " stretch %(a1)r erf %(a2)r difference %(diff)r"
+ transmit_rate = (min_rate + max_rate) / 2.0
+ while 1:
+ trial_number += 1
+ trial_duration = trial_number * self.trial_duration_per_trial
+ results = self.measure_and_compute(
+ trial_duration, transmit_rate, trial_result_list, max_rate,
+ integrator_data)
+ measurement, average, stdev, avg1, avg2, integrator_data = results
+ logging.info(message, {
+ "number": trial_number, "avg": average, "stdev": stdev,
+ "a1": avg1, "a2": avg2, "diff": avg2 - avg1})
+ if stop_time <= time.time():
+ return average, stdev
+ trial_result_list.append(measurement)
+ if (trial_number - self.trial_number_offset) <= 1:
+ next_load = max_rate
+ elif (trial_number - self.trial_number_offset) <= 3:
+ next_load = (measurement.receive_rate
+ / (1.0 - self.packet_loss_ratio_target))
+ else:
+ next_load = avg1 if trial_number % 2 else avg2
+ transmit_rate = min(max_rate, max(min_rate, next_load))
+
+ @staticmethod
+ def lfit_stretch(load, mrr, spread):
+ """Stretch-based fitting function.
+
+ Return the logarithm of average packet loss per second
+ when the load (argument) is offered to a system with given
+ mrr and spread (parameters).
+ Stretch function is 1/(1+Exp[-x]). The average itself is definite
+ integral from zero to load, of shifted and x-scaled stretch function.
+ As the integrator is sensitive to discontinuities,
+ and it calls this function at large areas of parameter space,
+ the implementation has to avoid rounding errors, overflows,
+ and correctly approximate underflows.
+
+ TODO: Explain how the high-level description
+ has been converted into an implementation full of ifs.
+
+ :param load: Offered load (positive), in packets per second.
+ :param mrr: Parameter of this fitting function, equal to limiting
+ (positive) average number of packets received (as opposed to lost)
+ when offered load is many spreads more than mrr.
+ :param spread: The x-scaling parameter (positive). No nice semantics,
+ roughly corresponds to size of "tail" for loads below mrr.
+ :type load: float
+ :type mrr: float
+ :type spread: float
+ :returns: Logarithm of average number of packets lost per second.
+ :rtype: float
+ """
+ # TODO: chi is from https://en.wikipedia.org/wiki/Nondimensionalization
+ chi = (load - mrr) / spread
+ chi0 = -mrr / spread
+# print "load", load, "mrr", mrr, "spread", spread, "chi", chi
+ if chi > 0:
+ log_lps = math.log(
+ load - mrr + (log_plus(0, -chi) - log_plus(0, chi0)) * spread)
+# print "big loss direct log_lps", log_lps
+ else:
+ approx = (math.exp(chi) - math.exp(2 * chi) / 2) * spread
+ if approx == 0.0:
+ log_lps = chi
+# print "small loss crude log_lps", log_lps
+ return log_lps
+ third = math.exp(3 * chi) / 3 * spread
+ if approx + third != approx + 2 * third:
+ log_lps = math.log(
+ (log_plus(0, chi) - log_plus(0, chi0)) * spread)
+# print "small loss direct log_lps", log_lps
+ else:
+ log_lps = math.log(
+ approx - (math.exp(chi0) - math.exp(2 * chi0)) * spread)
+# print "small loss approx log_lps", log_lps
+ return log_lps
+
+ @staticmethod
+ def lfit_erf(load, mrr, spread):
+ """Erf-based fitting function.
+
+ Return the logarithm of average packet loss per second
+ when the load (argument) is offered to a system with given
+ mrr and spread (parameters).
+ Erf function is Primitive function to normal distribution density.
+ The average itself is definite integral from zero to load,
+ of shifted and x-scaled erf function.
+ As the integrator is sensitive to discontinuities,
+ and it calls this function at large areas of parameter space,
+ the implementation has to avoid rounding errors, overflows,
+ and correctly approximate underflows.
+
+ TODO: Explain how the high-level description
+ has been converted into an implementation full of ifs.
+
+ :param load: Offered load (positive), in packets per second.
+ :param mrr: Parameter of this fitting function, equal to limiting
+ (positive) average number of packets received (as opposed to lost)
+ when offered load is many spreads more than mrr.
+ :param spread: The x-scaling parameter (positive). No nice semantics,
+ roughly corresponds to size of "tail" for loads below mrr.
+ :type load: float
+ :type mrr: float
+ :type spread: float
+ :returns: Logarithm of average number of packets lost per second.
+ :rtype: float
+ """
+ # Beware, this chi has the sign opposite to the stretch function chi.
+ # TODO: The stretch sign is just to have less minuses. Worth changing?
+ chi = (mrr - load) / spread
+ chi0 = mrr / spread
+# print "load", load, "mrr", mrr, "spread", spread,
+# print "chi", chi, "chi0", chi0
+ if chi >= -1.0:
+# print "positive, b ~> m"
+ if chi > math.exp(10):
+ first = PLRsearch.log_xerfcx_10 + 2 * (math.log(chi) - 10)
+# print "approximated"
+ else:
+ first = math.log(PLRsearch.xerfcx_limit - chi * erfcx(chi))
+# print "exact"
+ first -= chi * chi
+ second = math.log(PLRsearch.xerfcx_limit - chi * erfcx(chi0))
+ second -= chi0 * chi0
+ intermediate = log_minus(first, second)
+# print "first", first, "second", second,
+# print "intermediate", intermediate
+ else:
+# print "negative, b ~< m"
+ exp_first = PLRsearch.xerfcx_limit + chi * erfcx(-chi)
+ exp_first *= math.exp(-chi * chi)
+ exp_first -= 2 * chi
+ # TODO: Why has the following line chi there (as opposed to chi0)?
+ # In general the functions would be more readable if they explicitly
+ # return math.log(func(chi) - func(chi0))
+ # for some function "func", at least for some branches.
+ second = math.log(PLRsearch.xerfcx_limit - chi * erfcx(chi0))
+ second -= chi0 * chi0
+ intermediate = math.log(exp_first - math.exp(second))
+# print "exp_first", exp_first, "second", second,
+# print "intermediate", intermediate
+ result = intermediate + math.log(spread) - math.log(erfc(-chi0))
+# print "lfit erf result", result
+ return result
+
+ @staticmethod
+ def find_critical_rate(lfit_func, log_lps_target, mrr, spread):
+ """Given lps target and parameters, return the achieving offered load.
+
+ This is basically an inverse function to lfit_func
+ when parameters are fixed.
+ Instead of implementing effective implementation
+ of the inverse function, this implementation uses
+ brute force binary search.
+ The search starts at an arbitrary offered load,
+ uses exponential search to find the other bound
+ and then bisecting until the load is found (or interval
+ becoming degenerate).
+
+ TODO: Use at least some method with faster convergence.
+
+ :param lfit_func: Fitting function, typically lfit_spread or lfit_erf.
+ :param log_lps_target: Fitting function should return this
+ at the returned load and parameters.
+ :param mrr: The mrr parameter for the fitting function.
+ :param spread: The spread parameter for the fittinmg function.
+ :type lfit_func: Function from 3 floats to float.
+ :type log_lps_target: float
+ :type mrr: float
+ :type spread: float
+ :returns: Load [pps] which achieves the target with given parameters.
+ :rtype: float
+ """
+ # TODO: Should we make the initial rate configurable?
+ rate = 10000000.0
+ log_loss = lfit_func(rate, mrr, spread)
+ if log_loss == log_lps_target:
+ return rate
+ # Exponential search part.
+ if log_loss > log_lps_target:
+ rate_hi = rate
+ while 1:
+ rate_lo = rate_hi / 2.0
+ log_loss = lfit_func(rate_lo, mrr, spread)
+ if log_loss > log_lps_target:
+ rate_hi = rate_lo
+ continue
+ if log_loss == log_lps_target:
+ return rate_lo
+ break
+ else:
+ rate_lo = rate
+ while 1:
+ rate_hi = rate_lo * 2.0
+ log_loss = lfit_func(rate_hi, mrr, spread)
+ if log_loss < log_lps_target:
+ rate_lo = rate_hi
+ continue
+ if log_loss == log_lps_target:
+ return rate_hi
+ break
+ # Binary search part.
+ while rate_hi != rate_lo:
+ rate = (rate_hi + rate_lo) / 2.0
+ log_loss = lfit_func(rate, mrr, spread)
+ if rate == rate_hi or rate == rate_lo or log_loss == log_lps_target:
+# print "found", rate
+ return rate
+ if log_loss > log_lps_target:
+ rate_hi = rate
+ else:
+ rate_lo = rate
+
+ @staticmethod
+ def log_weight(lfit_func, trial_result_list, mrr, spread):
+ """Return log of weight of trial results by the function and parameters.
+
+ Integrator assumes uniform distribution, but over different parameters.
+ Weight and likelihood are used interchangeably here anyway.
+
+ Each trial has an offered load, a duration and a loss count.
+ Fitting function is used to compute the average loss per second.
+ Poisson distribution (with average loss per trial) is used
+ to get likelihood of one trial result, the overal likelihood is product.
+ As likelihoods can be extremely small, logarithms are tracked instead.
+
+ TODO: Copy ReceiveRateMeasurement from MLRsearch.
+
+ :param lfit_func: Fitting function, typically lfit_spread or lfit_erf.
+ :param result_list: List of trial measurement results.
+ :param mrr: The mrr parameter for the fitting function.
+ :param spread: The spread parameter for the fittinmg function.
+ :type lfit_func: Function from 3 floats to float.
+ :type result_list: list of MLRsearch.ReceiveRateMeasurement
+ :type mrr: float
+ :type spread: float
+ :returns: Logarithm of result weight for given function and parameters.
+ :rtype: float
+ """
+ log_likelihood = 0.0
+ for result in trial_result_list:
+# print "DEBUG for tr", result.target_tr,
+# print "lc", result.loss_count, "d", result.duration
+ log_avg_loss_per_second = lfit_func(result.target_tr, mrr, spread)
+ log_avg_loss_per_trial = (
+ log_avg_loss_per_second + math.log(result.duration))
+ # Poisson probability computation works nice for logarithms.
+ log_trial_likelihood = (
+ result.loss_count * log_avg_loss_per_trial
+ - math.exp(log_avg_loss_per_trial))
+ log_trial_likelihood -= math.lgamma(1 + result.loss_count)
+ log_likelihood += log_trial_likelihood
+# print "log_avg_loss_per_second", log_avg_loss_per_second
+# print "log_avg_loss_per_trial", log_avg_loss_per_trial
+# print "avg_loss_per_trial", math.exp(log_avg_loss_per_trial)
+# print "log_trial_likelihood", log_trial_likelihood
+# print "log_likelihood", log_likelihood
+# print "returning log_likelihood", log_likelihood
+ return log_likelihood
+
+ # FIXME: Refactor (somehow) so pylint stops complaining about
+ # too many local variables.
+ # Some work already done in subsequent changes.
+ def measure_and_compute(
+ self, trial_duration, transmit_rate,
+ trial_result_list, max_rate, integrator_data):
+ """Perform both measurement and computation at once.
+
+ High level steps: Prepare and launch computation worker processes,
+ perform the measurement, stop computation and combine results.
+
+ Integrator needs a specific function to process (-1, 1) parameters.
+ As our fitting functions use dimensional parameters,
+ so a transformation is performed, resulting in a specific prior
+ distribution over the dimensional parameters.
+ Maximal rate (line rate) is needed for that transformation.
+
+ Two fitting functions are used, computation is started
+ on temporary worker process per fitting function. After the measurement,
+ Average and stdev of the critical rate (not log) of each worker
+ are combined and returned. Raw averages are also returned,
+ offered load for next iteration is chosen from them.
+ The idea is that one fitting function might be fitting much better,
+ measurements at its avg are best for relevant results (for both),
+ but we do not know which fitting function it is.
+
+ TODO: Define class for integrator data, so that fields are documented.
+ TODO: Define class for result object, so that fields are documented.
+ TODO: More processes to the slower fitting function?
+ TODO: Re-use processes, instead creating on each computation.
+ TODO: As only one result is needed fresh, figure out a way
+ how to keep the other worker running. This will alow shorter
+ duration per trial. Special handling at first and last measurement
+ will be needed (to properly initialize and to properly combine results).
+
+ :param trial_duration: Length of the measurement in seconds.
+ :param transmit_rate: Offered load in packets per second.
+ :param trial_result_list: Results of previous measurements.
+ :param max_rate: Theoretic maximum of possible ofered load.
+ :param integrator_data: Hints to speed up the numeric computation.
+ :type trial_duration: float
+ :type transmit_rate: float
+ :type trial_result_list: list of MLRsearch.ReceiveRateMeasurement
+ :type max_rate: float
+ :type integrator_data: 4-tuple of gaussian positions and covariances
+ :returns: Measurement and computation results.
+ :rtype: 6-tuple: ReceiveRateMeasurement, floats, integrator data.
+ """
+ # Preparation phase.
+ dimension = 2
+ stretch_bias_avg, erf_bias_avg, stretch_bias_cov, erf_bias_cov = (
+ integrator_data)
+ packet_loss_per_second_target = (
+ transmit_rate * self.packet_loss_ratio_target)
+ def start_computing(fitting_function, bias_avg, bias_cov):
+ """Just a block of code to be used for each fitting function.
+
+ Define function for integrator, create process and pipe ends,
+ start computation, return the boss pipe end.
+
+ :param fitting_function: lfit_erf or lfit_stretch.
+ :param bias_avg: Tuple of floats to start searching around.
+ :param bias_cov: Covariance matrix defining initial focus shape.
+ :type fitting_function: Function from 3 floats to float.
+ :type bias_avg: 2-tuple of floats
+ :type bias_cov: 2-tuple of 2-tuples of floats
+ :returns: Boss end of communication pipe.
+ :rtype: multiprocessing.Connection
+ """
+ def value_logweight_func(x_mrr, x_spread):
+ """Return log of critical rate and log of likelihood.
+
+ This is a closure. The ancestor function got
+ trial_result_list as a parameter, and we are accessing it.
+ As integrator has strict conditions on function signature,
+ trial_result_list cannot be an explicit argument
+ of the current function.
+ This is also why we have to define this closure
+ at each invocation of the ancestor function anew.
+
+ The dimensional spread parameter is the (dimensional) mrr
+ raised to the power of x_spread scaled to interval (0, 1).
+ The dimensional mrr parameter distribution has shape of
+ 1/(1+x^2), but x==1 corresponds to max_rate
+ and 1.0 pps is added to avoid numerical problems in fitting
+ functions.
+
+ TODO: x^-2 (for x>1.0) might be simpler/nicer prior.
+
+ :param x_mrr: The first dimensionless param
+ from (-1, 1) interval.
+ :param x_spread: The second dimensionless param
+ from (-1, 1) interval.
+ :returns: Log of critical rate [pps] and log of likelihood.
+ :rtype: 2-tuple of float
+ """
+ mrr = max_rate * (1.0 / (x_mrr + 1.0) - 0.5) + 1.0
+ spread = math.exp((x_spread + 1.0) / 2.0 * math.log(mrr))
+# print "mrr", mrr, "spread", spread
+ logweight = self.log_weight(
+ fitting_function, trial_result_list, mrr, spread)
+ value = math.log(self.find_critical_rate(
+ fitting_function,
+ math.log(packet_loss_per_second_target), mrr, spread))
+ return value, logweight
+ dilled_function = dill.dumps(value_logweight_func)
+ boss_pipe_end, worker_pipe_end = multiprocessing.Pipe()
+ boss_pipe_end.send(
+ (dimension, dilled_function, bias_avg, bias_cov))
+ worker = multiprocessing.Process(
+ target=Integrator.try_estimate_nd, args=(worker_pipe_end,))
+ worker.daemon = True
+ worker.start()
+ return boss_pipe_end
+ erf_pipe = start_computing(
+ self.lfit_erf, erf_bias_avg, erf_bias_cov)
+ stretch_pipe = start_computing(
+ self.lfit_stretch, stretch_bias_avg, stretch_bias_cov)
+ # Measurement phase.
+ measurement = self.measurer.measure(trial_duration, transmit_rate)
+ # Processing phase.
+ erf_pipe.send(None)
+ stretch_pipe.send(None)
+ if not stretch_pipe.poll(1.0):
+ raise RuntimeError("Stretch worker did not finish!")
+ result_or_traceback = stretch_pipe.recv()
+ try:
+ (stretch_avg, stretch_stdev, stretch_bias_avg,
+ stretch_bias_cov, debug_list, _) = result_or_traceback
+ except ValueError:
+ raise RuntimeError(
+ "Stretch worker failed with the following traceback:\n{tr}"
+ .format(tr=result_or_traceback))
+ logging.info("Logs from stretch worker:")
+ for message in debug_list:
+ logging.debug(message)
+ if not erf_pipe.poll(1.0):
+ raise RuntimeError("Erf worker did not finish!")
+ result_or_traceback = erf_pipe.recv()
+ try:
+ (erf_avg, erf_stdev, erf_bias_avg,
+ erf_bias_cov, debug_list, _) = result_or_traceback
+ except ValueError:
+ raise RuntimeError(
+ "Erf worker failed with the following traceback:\n{tr}"
+ .format(tr=result_or_traceback))
+ logging.info("Logs from erf worker:")
+ for message in debug_list:
+ logging.debug(message)
+ avg = math.exp((stretch_avg + erf_avg) / 2.0)
+ var = (stretch_stdev * stretch_stdev + erf_stdev * erf_stdev) / 2.0
+ var += (stretch_avg - erf_avg) * (stretch_avg - erf_avg) / 4.0
+ stdev = avg * math.sqrt(var)
+ integrator_data = (
+ stretch_bias_avg, erf_bias_avg, stretch_bias_cov, erf_bias_cov)
+ return (
+ measurement, avg, stdev, math.exp(stretch_avg),
+ math.exp(erf_avg), integrator_data)
Code Review / csit.git / blobdiff