1 # Copyright (c) 2024 Cisco and/or its affiliates.
2 # Licensed under the Apache License, Version 2.0 (the "License");
3 # you may not use this file except in compliance with the License.
4 # You may obtain a copy of the License at:
6 # http://www.apache.org/licenses/LICENSE-2.0
8 # Unless required by applicable law or agreed to in writing, software
9 # distributed under the License is distributed on an "AS IS" BASIS,
10 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
11 # See the License for the specific language governing permissions and
12 # limitations under the License.
14 """Module holding PLRsearch class."""
18 import multiprocessing
21 from collections import namedtuple
25 from scipy.special import erfcx, erfc
27 # TODO: Teach FD.io CSIT to use multiple dirs in PYTHONPATH,
28 # then switch to absolute imports within PLRsearch package.
29 # Current usage of relative imports is just a short term workaround.
30 from . import Integrator
31 from . import stat_trackers
32 from .log_plus import log_plus, log_minus
36 """A class to encapsulate data relevant for the search method.
38 The context is performance testing of packet processing systems.
39 The system, when being offered a steady stream of packets,
40 can process some of them successfully, other are considered "lost".
42 See docstring of the search method for algorithm description.
44 Two constants are stored as class fields for speed.
46 Method other than search (and than __init__)
47 are just internal code structure.
49 TODO: Those method names should start with underscore then.
52 xerfcx_limit = math.pow(math.acos(0), -0.5)
53 log_xerfcx_10 = math.log(xerfcx_limit - math.exp(10) * erfcx(math.exp(10)))
58 trial_duration_per_trial,
59 packet_loss_ratio_target,
60 trial_number_offset=0,
64 """Store rate measurer and additional parameters.
66 The measurer must never report negative loss count.
68 TODO: Copy AbstractMeasurer from MLRsearch.
70 :param measurer: The measurer to call when searching.
71 :param trial_duration_per_trial: Each trial has larger duration
72 than the previous trial. This is the increment, in seconds.
73 :param packet_loss_ratio_target: The algorithm tries to estimate
74 the offered load leading to this ratio on average.
75 Trial ratio is number of packets lost divided by packets offered.
76 :param trial_number_offset: The "first" trial number will be 1+this.
77 Use this to ensure first iterations have enough time to compute
78 reasonable estimates for later trials to use.
79 :param timeout: The search ends if it lasts more than this many seconds.
80 :type measurer: MLRsearch.AbstractMeasurer
81 :type trial_duration_per_trial: float
82 :type packet_loss_ratio_target: float
83 :type trial_number_offset: int
86 self.measurer = measurer
87 self.trial_duration_per_trial = float(trial_duration_per_trial)
88 self.packet_loss_ratio_target = float(packet_loss_ratio_target)
89 self.trial_number_offset = int(trial_number_offset)
90 self.timeout = float(timeout)
91 self.trace_enabled = bool(trace_enabled)
93 def search(self, min_rate, max_rate):
94 """Perform the search, return average and stdev for throughput estimate.
96 Considering measurer and packet_loss_ratio_target (see __init__),
97 find such an offered load (called critical load) that is expected
98 to hit the target loss ratio in the limit of very long trial duration.
99 As the system is probabilistic (and test duration is finite),
100 the critical ratio is only estimated.
101 Return the average and standard deviation of the estimate.
103 In principle, this algorithm performs trial measurements,
104 each with varied offered load (which is constant during the trial).
105 During each measurement, Bayesian inference is performed
106 on all the measurement results so far.
107 When timeout is up, the last estimate is returned,
108 else another trial is performed.
110 It is assumed that the system under test, even though not deterministic,
111 still follows the rule of large numbers. In another words,
112 any growing set of measurements at a particular offered load
113 will converge towards unique (for the given load) packet loss ratio.
114 This means there is a deterministic (but unknown) function
115 mapping the offered load to average loss ratio.
116 This function is called loss ratio function.
117 This also assumes the average loss ratio
118 does not depend on trial duration.
120 The actual probability distribution of loss counts, achieving
121 the average ratio on trials of various duration
122 can be complicated (and can depend on offered load), but simply assuming
123 Poisson distribution will make the algorithm converge.
124 Binomial distribution would be more precise,
125 but Poisson is more practical, as it effectively gives
126 less information content to high ratio results.
128 Even when applying other assumptions on the loss ratio function
129 (increasing function, limit zero ratio when load goes to zero,
130 global upper limit on rate of packets processed), there are still
131 too many different shapes of possible loss functions,
132 which makes full Bayesian reasoning intractable.
134 This implementation radically simplifies things by examining
135 only two shapes, each with finitely many (in this case just two)
136 parameters. In other words, two fitting functions
137 (each with two parameters and one argument).
138 When restricting model space to one of the two fitting functions,
139 the Bayesian inference becomes tractable (even though it needs
140 numerical integration from Integrator class).
142 The first measurement is done at the middle between
143 min_rate and max_rate, to help with convergence
144 if max_rate measurements give loss below target.
145 TODO: Fix overflow error and use min_rate instead of the middle.
147 The second measurement is done at max_rate, next few measurements
148 have offered load of previous load minus excess loss rate.
149 This simple rule is found to be good when offered loads
150 so far are way above the critical rate. After few measurements,
151 inference from fitting functions converges faster that this initial
152 "optimistic" procedure.
154 Offered loads close to (limiting) critical rate are the most useful,
155 as linear approximation of the fitting function
156 becomes good enough there (thus reducing the impact
157 of the overall shape of fitting function).
158 After several trials, usually one of the fitting functions
159 has better predictions than the other one, but the algorithm
160 does not track that. Simply, it uses the estimate average,
161 alternating between the functions.
162 Multiple workarounds are applied to try and avoid measurements
163 both in zero loss region and in big loss region,
164 as their results tend to make the critical load estimate worse.
166 The returned average and stdev is a combination of the two fitting
169 :param min_rate: Avoid measuring at offered loads below this,
170 in packets per second.
171 :param max_rate: Avoid measuring at offered loads above this,
172 in packets per second.
173 :type min_rate: float
174 :type max_rate: float
175 :returns: Average and stdev of critical load estimate.
176 :rtype: 2-tuple of float
178 stop_time = time.time() + self.timeout
179 min_rate = float(min_rate)
180 max_rate = float(max_rate)
182 f"Started search with min_rate {min_rate!r}, "
183 f"max_rate {max_rate!r}"
185 trial_result_list = []
186 trial_number = self.trial_number_offset
187 focus_trackers = (None, None)
188 transmit_rate = (min_rate + max_rate) / 2.0
189 lossy_loads = [max_rate]
190 zeros = 0 # How many consecutive zero loss results are happening.
193 logging.info(f"Trial {trial_number!r}")
194 results = self.measure_and_compute(
195 self.trial_duration_per_trial * trial_number,
202 measurement, average, stdev, avg1, avg2, focus_trackers = results
204 # TODO: Ratio of fill rate to drain rate seems to have
205 # exponential impact. Make it configurable, or is 4:3 good enough?
206 if measurement.loss_ratio >= self.packet_loss_ratio_target:
207 for _ in range(4 * zeros):
208 lossy_loads.append(measurement.intended_load)
209 if measurement.loss_ratio > 0.0:
212 if stop_time <= time.time():
213 return average, stdev
214 trial_result_list.append(measurement)
215 if (trial_number - self.trial_number_offset) <= 1:
217 elif (trial_number - self.trial_number_offset) <= 3:
218 next_load = measurement.relative_forwarding_rate / (
219 1.0 - self.packet_loss_ratio_target
222 next_load = (avg1 + avg2) / 2.0
224 if lossy_loads[0] > next_load:
225 diminisher = math.pow(2.0, 1 - zeros)
226 next_load = lossy_loads[0] + diminisher * next_load
227 next_load /= 1.0 + diminisher
228 # On zero measurement, we need to drain obsoleted low losses
229 # even if we did not use them to increase next_load,
230 # in order to get to usable loses at higher loads.
231 if len(lossy_loads) > 3:
232 lossy_loads = lossy_loads[3:]
234 f"Zeros {zeros!r} orig {(avg1 + avg2) / 2.0!r} "
235 f"next {next_load!r} loads {lossy_loads!r}"
237 transmit_rate = min(max_rate, max(min_rate, next_load))
240 def lfit_stretch(trace, load, mrr, spread):
241 """Stretch-based fitting function.
243 Return the logarithm of average packet loss per second
244 when the load (argument) is offered to a system with given
245 mrr and spread (parameters).
246 Stretch function is 1/(1+Exp[-x]). The average itself is definite
247 integral from zero to load, of shifted and x-scaled stretch function.
248 As the integrator is sensitive to discontinuities,
249 and it calls this function at large areas of parameter space,
250 the implementation has to avoid rounding errors, overflows,
251 and correctly approximate underflows.
253 TODO: Explain how the high-level description
254 has been converted into an implementation full of ifs.
256 :param trace: A multiprocessing-friendly logging function (closure).
257 :param load: Offered load (positive), in packets per second.
258 :param mrr: Parameter of this fitting function, equal to limiting
259 (positive) average number of packets received (as opposed to lost)
260 when offered load is many spreads more than mrr.
261 :param spread: The x-scaling parameter (positive). No nice semantics,
262 roughly corresponds to size of "tail" for loads below mrr.
263 :type trace: function (str, object) -> NoneType
267 :returns: Logarithm of average number of packets lost per second.
270 # TODO: What is the fastest way to use such values?
273 log_spread = math.log(spread)
274 # TODO: chi is from https://en.wikipedia.org/wiki/Nondimensionalization
275 chi = (load - mrr) / spread
277 trace("stretch: load", load)
279 trace("spread", spread)
284 load - mrr + (log_plus(0, -chi) - log_plus(0, chi0)) * spread
286 trace("big loss direct log_lps", log_lps)
288 two_positive = log_plus(chi, 2 * chi0 - log_2)
289 two_negative = log_plus(chi0, 2 * chi - log_2)
290 if two_positive <= two_negative:
291 log_lps = log_minus(chi, chi0) + log_spread
292 trace("small loss crude log_lps", log_lps)
294 two = log_minus(two_positive, two_negative)
295 three_positive = log_plus(two_positive, 3 * chi - log_3)
296 three_negative = log_plus(two_negative, 3 * chi0 - log_3)
297 three = log_minus(three_positive, three_negative)
299 log_lps = two + log_spread
300 trace("small loss approx log_lps", log_lps)
302 log_lps = math.log(log_plus(0, chi) - log_plus(0, chi0))
303 log_lps += log_spread
304 trace("small loss direct log_lps", log_lps)
308 def lfit_erf(trace, load, mrr, spread):
309 """Erf-based fitting function.
311 Return the logarithm of average packet loss per second
312 when the load (argument) is offered to a system with given
313 mrr and spread (parameters).
314 Erf function is Primitive function to normal distribution density.
315 The average itself is definite integral from zero to load,
316 of shifted and x-scaled erf function.
317 As the integrator is sensitive to discontinuities,
318 and it calls this function at large areas of parameter space,
319 the implementation has to avoid rounding errors, overflows,
320 and correctly approximate underflows.
322 TODO: Explain how the high-level description
323 has been converted into an implementation full of ifs.
325 :param trace: A multiprocessing-friendly logging function (closure).
326 :param load: Offered load (positive), in packets per second.
327 :param mrr: Parameter of this fitting function, equal to limiting
328 (positive) average number of packets received (as opposed to lost)
329 when offered load is many spreads more than mrr.
330 :param spread: The x-scaling parameter (positive). No nice semantics,
331 roughly corresponds to size of "tail" for loads below mrr.
332 :type trace: function (str, object) -> NoneType
336 :returns: Logarithm of average number of packets lost per second.
339 # Beware, this chi has the sign opposite to the stretch function chi.
340 # TODO: The stretch sign is just to have less minuses. Worth changing?
341 chi = (mrr - load) / spread
343 trace("Erf: load", load)
345 trace("spread", spread)
349 trace("positive, b roughly bigger than m", None)
350 if chi > math.exp(10):
351 first = PLRsearch.log_xerfcx_10 + 2 * (math.log(chi) - 10)
352 trace("approximated first", first)
354 first = math.log(PLRsearch.xerfcx_limit - chi * erfcx(chi))
355 trace("exact first", first)
357 second = math.log(PLRsearch.xerfcx_limit - chi * erfcx(chi0))
358 second -= chi0 * chi0
359 intermediate = log_minus(first, second)
360 trace("first", first)
362 trace("negative, b roughly smaller than m", None)
363 exp_first = PLRsearch.xerfcx_limit + chi * erfcx(-chi)
364 exp_first *= math.exp(-chi * chi)
366 # TODO: Why has the following line chi there (as opposed to chi0)?
367 # In general the functions would be more readable if they explicitly
368 # return math.log(func(chi) - func(chi0))
369 # for some function "func", at least for some branches.
370 second = math.log(PLRsearch.xerfcx_limit - chi * erfcx(chi0))
371 second -= chi0 * chi0
372 intermediate = math.log(exp_first - math.exp(second))
373 trace("exp_first", exp_first)
374 trace("second", second)
375 trace("intermediate", intermediate)
376 result = intermediate + math.log(spread) - math.log(erfc(-chi0))
377 trace("result", result)
381 def find_critical_rate(
382 trace, lfit_func, min_rate, max_rate, loss_ratio_target, mrr, spread
384 """Given ratio target and parameters, return the achieving offered load.
386 This is basically an inverse function to lfit_func
387 when parameters are fixed.
388 Instead of implementing effective implementation
389 of the inverse function, this implementation uses
390 brute force binary search. It is bisecting (nim_rate, max_rate) interval
391 until the critical load is found (or interval becomes degenerate).
392 This implementation assures min and max rate limits are honored.
394 TODO: Use some method with faster convergence?
396 :param trace: A multiprocessing-friendly logging function (closure).
397 :param lfit_func: Fitting function, typically lfit_spread or lfit_erf.
398 :param min_rate: Lower bound for binary search [pps].
399 :param max_rate: Upper bound for binary search [pps].
400 :param loss_ratio_target: Fitting function should return loss rate
401 giving this ratio at the returned load and parameters [1].
402 :param mrr: The mrr parameter for the fitting function [pps].
403 :param spread: The spread parameter for the fittinmg function [pps].
404 :type trace: function (str, object) -> None
405 :type lfit_func: Function from 3 floats to float.
406 :type min_rate: float
407 :type max_rate: float
408 :type loss_ratio_target: float
411 :returns: Load [pps] which achieves the target with given parameters.
414 trace("Finding critical rate for loss_ratio_target", loss_ratio_target)
418 while loss_ratio != loss_ratio_target:
419 rate = (rate_hi + rate_lo) / 2.0
420 if rate in (rate_hi, rate_lo):
422 loss_rate = math.exp(lfit_func(trace, rate, mrr, spread))
423 loss_ratio = loss_rate / rate
424 if loss_ratio > loss_ratio_target:
425 trace("halving down", rate)
427 elif loss_ratio < loss_ratio_target:
428 trace("halving up", rate)
434 def log_weight(trace, lfit_func, trial_result_list, mrr, spread):
435 """Return log of weight of trial results by the function and parameters.
437 Integrator assumes uniform distribution, but over different parameters.
438 Weight and likelihood are used interchangeably here anyway.
440 Each trial has an intended load, a sent count and a loss count
441 (probably counting unsent packets as loss, as they signal
442 the load is too high for the traffic generator).
443 The fitting function is used to compute the average loss rate.
444 Geometric distribution (with average loss per trial) is used
445 to get likelihood of one trial result, the overal likelihood
446 is a product of all trial likelihoods.
447 As likelihoods can be extremely small, logarithms are tracked instead.
449 The current implementation does not use direct loss rate
450 from the fitting function, as the input and output units may not match
451 (e.g. intended load in TCP transactions, loss in packets).
452 Instead, the expected average loss is scaled according to the number
453 of packets actually sent.
455 TODO: Copy MeasurementResult from MLRsearch.
457 :param trace: A multiprocessing-friendly logging function (closure).
458 :param lfit_func: Fitting function, typically lfit_spread or lfit_erf.
459 :param trial_result_list: List of trial measurement results.
460 :param mrr: The mrr parameter for the fitting function.
461 :param spread: The spread parameter for the fitting function.
462 :type trace: function (str, object) -> None
463 :type lfit_func: Function from 3 floats to float.
464 :type trial_result_list: list of MLRsearch.MeasurementResult
467 :returns: Logarithm of result weight for given function and parameters.
471 trace("log_weight for mrr", mrr)
472 trace("spread", spread)
473 for result in trial_result_list:
474 trace("for tr", result.intended_load)
475 trace("lc", result.loss_count)
476 trace("d", result.intended_duration)
477 # _rel_ values use units of intended_load (transactions per second).
478 log_avg_rel_loss_per_second = lfit_func(
479 trace, result.intended_load, mrr, spread
481 # _abs_ values use units of loss count (maybe packets).
482 # There can be multiple packets per transaction.
483 log_avg_abs_loss_per_trial = log_avg_rel_loss_per_second + math.log(
484 result.offered_count / result.intended_load
486 # Geometric probability computation for logarithms.
487 log_trial_likelihood = log_plus(0.0, -log_avg_abs_loss_per_trial)
488 log_trial_likelihood *= -result.loss_count
489 log_trial_likelihood -= log_plus(0.0, +log_avg_abs_loss_per_trial)
490 log_likelihood += log_trial_likelihood
491 trace("avg_loss_per_trial", math.exp(log_avg_abs_loss_per_trial))
492 trace("log_trial_likelihood", log_trial_likelihood)
493 return log_likelihood
495 def measure_and_compute(
502 focus_trackers=(None, None),
505 """Perform both measurement and computation at once.
507 High level steps: Prepare and launch computation worker processes,
508 perform the measurement, stop computation and combine results.
510 Integrator needs a specific function to process (-1, 1) parameters.
511 As our fitting functions use dimensional parameters,
512 so a transformation is performed, resulting in a specific prior
513 distribution over the dimensional parameters.
514 Maximal rate (line rate) is needed for that transformation.
516 Two fitting functions are used, computation is started
517 on temporary worker process per fitting function. After the measurement,
518 average and stdev of the critical rate (not log) of each worker
519 are combined and returned. Raw averages are also returned,
520 offered load for next iteration is chosen based on them.
521 The idea is that one fitting function might be fitting much better,
522 measurements at its avg are best for relevant results (for both),
523 but we do not know which fitting function it is.
525 Focus trackers are updated in-place. If a focus tracker in None,
526 new instance is created.
528 TODO: Define class for result object, so that fields are documented.
529 TODO: Re-use processes, instead creating on each computation?
530 TODO: As only one result is needed fresh, figure out a way
531 how to keep the other worker running. This will alow shorter
532 duration per trial. Special handling at first and last measurement
533 will be needed (to properly initialize and to properly combine results).
535 :param trial_duration: Length of the measurement in seconds.
536 :param transmit_rate: Offered load in packets per second.
537 :param trial_result_list: Results of previous measurements.
538 :param min_rate: Practical minimum of possible ofered load.
539 :param max_rate: Practical maximum of possible ofered load.
540 :param focus_trackers: Pair of trackers initialized
541 to speed up the numeric computation.
542 :param max_samples: Limit for integrator samples, for debugging.
543 :type trial_duration: float
544 :type transmit_rate: float
545 :type trial_result_list: list of MLRsearch.MeasurementResult
546 :type min_rate: float
547 :type max_rate: float
548 :type focus_trackers: 2-tuple of None or stat_trackers.VectorStatTracker
549 :type max_samples: None or int
550 :returns: Measurement and computation results.
551 :rtype: _ComputeResult
554 f"measure_and_compute started with self {self!r}, trial_duration "
555 f"{trial_duration!r}, transmit_rate {transmit_rate!r}, "
556 f"trial_result_list {trial_result_list!r}, max_rate {max_rate!r}, "
557 f"focus_trackers {focus_trackers!r}, max_samples {max_samples!r}"
561 stretch_focus_tracker, erf_focus_tracker = focus_trackers
562 if stretch_focus_tracker is None:
563 stretch_focus_tracker = stat_trackers.VectorStatTracker(dimension)
564 stretch_focus_tracker.unit_reset()
565 if erf_focus_tracker is None:
566 erf_focus_tracker = stat_trackers.VectorStatTracker(dimension)
567 erf_focus_tracker.unit_reset()
568 old_trackers = stretch_focus_tracker.copy(), erf_focus_tracker.copy()
570 def start_computing(fitting_function, focus_tracker):
571 """Just a block of code to be used for each fitting function.
573 Define function for integrator, create process and pipe ends,
574 start computation, return the boss pipe end.
576 :param fitting_function: lfit_erf or lfit_stretch.
577 :param focus_tracker: Tracker initialized to speed up the numeric
579 :type fitting_function: Function from 3 floats to float.
580 :type focus_tracker: None or stat_trackers.VectorStatTracker
581 :returns: Boss end of communication pipe.
582 :rtype: multiprocessing.Connection
585 boss_pipe_end, worker_pipe_end = multiprocessing.Pipe()
586 # Starting the worker first. Contrary to documentation
587 # https://docs.python.org/3/library/multiprocessing.html#multiprocessing.connection.Connection
588 # sending of large object without active listener on the other side
589 # results in a deadlock, not in a ValueError.
590 # See https://stackoverflow.com/questions/15137292/large-objects-and-multiprocessing-pipes-and-send
591 worker = multiprocessing.Process(
592 target=Integrator.try_estimate_nd,
593 args=(worker_pipe_end, 5.0, self.trace_enabled),
598 # Only now it is safe to send the function to compute with.
599 def value_logweight_func(trace, x_mrr, x_spread):
600 """Return log of critical rate and log of likelihood.
602 This is a closure. The ancestor function got
603 trial_result_list as a parameter, and we are accessing it.
604 As integrator has strict conditions on function signature,
605 trial_result_list cannot be an explicit argument
606 of the current function.
607 This is also why we have to define this closure
608 at each invocation of the ancestor function anew.
610 The dimensional spread parameter is the (dimensional) mrr
611 raised to the power of x_spread scaled to interval (0, 1).
612 The dimensional mrr parameter distribution has shape of
613 1/(1+x^2), but x==1 corresponds to max_rate
614 and 1.0 pps is added to avoid numerical problems in fitting
617 TODO: x^-2 (for x>1.0) might be simpler/nicer prior.
619 :param trace: Multiprocessing-safe logging function (closure).
620 :param x_mrr: The first dimensionless param
621 from (-1, 1) interval.
622 :param x_spread: The second dimensionless param
623 from (-1, 1) interval.
624 :type trace: function (str, object) -> None
626 :type x_spread: float
627 :returns: Log of critical rate [pps] and log of likelihood.
628 :rtype: 2-tuple of float
630 mrr = max_rate * (1.0 / (x_mrr + 1.0) - 0.5) + 1.0
631 spread = math.exp((x_spread + 1.0) / 2.0 * math.log(mrr))
632 logweight = self.log_weight(
633 trace, fitting_function, trial_result_list, mrr, spread
636 self.find_critical_rate(
641 self.packet_loss_ratio_target,
646 return value, logweight
648 dilled_function = dill.dumps(value_logweight_func)
650 (dimension, dilled_function, focus_tracker, max_samples)
654 erf_pipe = start_computing(self.lfit_erf, erf_focus_tracker)
655 stretch_pipe = start_computing(self.lfit_stretch, stretch_focus_tracker)
658 measurement = self.measurer.measure(trial_duration, transmit_rate)
661 def stop_computing(name, pipe):
662 """Just a block of code to be used for each worker.
664 Send stop object, poll for result, then either
665 unpack response, log messages and return, or raise traceback.
667 TODO: Define class/structure for the return value?
669 :param name: Human friendly worker identifier for logging purposes.
670 :param pipe: Boss end of connection towards worker to stop.
672 :type pipe: multiprocessing.Connection
673 :returns: Computed value tracker, actual focus tracker,
674 and number of samples used for this iteration.
675 :rtype: _PartialResult
677 # If worker encountered an exception, we get it in the recv below,
678 # but send will report a broken pipe.
679 # EAFP says we should ignore the error (instead of polling first).
680 # https://devblogs.microsoft.com/python
681 # /idiomatic-python-eafp-versus-lbyl/
684 except BrokenPipeError:
686 if not pipe.poll(10.0):
687 raise RuntimeError(f"Worker {name} did not finish!")
688 result_or_traceback = pipe.recv()
696 ) = result_or_traceback
697 except ValueError as exc:
699 f"Worker {name} failed with the following traceback:\n"
700 f"{result_or_traceback}"
702 logging.info(f"Logs from worker {name!r}:")
703 for message in debug_list:
704 logging.info(message)
705 for message in trace_list:
706 logging.debug(message)
708 f"trackers: value {value_tracker!r} focus {focus_tracker!r}"
710 return _PartialResult(value_tracker, focus_tracker, sampls)
712 stretch_result = stop_computing("stretch", stretch_pipe)
713 erf_result = stop_computing("erf", erf_pipe)
714 result = PLRsearch._get_result(measurement, stretch_result, erf_result)
716 f"measure_and_compute finished with trial result "
717 f"{result.measurement!r} avg {result.avg!r} stdev {result.stdev!r} "
718 f"stretch {result.stretch_exp_avg!r} erf {result.erf_exp_avg!r} "
719 f"new trackers {result.trackers!r} old trackers {old_trackers!r} "
720 f"stretch samples {stretch_result.samples!r} erf samples "
721 f"{erf_result.samples!r}"
726 def _get_result(measurement, stretch_result, erf_result):
727 """Process and collate results from measure_and_compute.
729 Turn logarithm based values to exponential ones,
730 combine averages and stdevs of two fitting functions into a whole.
732 :param measurement: The trial measurement obtained during computation.
733 :param stretch_result: Computation output for stretch fitting function.
734 :param erf_result: Computation output for erf fitting function.
735 :type measurement: MeasurementResult
736 :type stretch_result: _PartialResult
737 :type erf_result: _PartialResult
738 :returns: Combined results.
739 :rtype: _ComputeResult
741 stretch_avg = stretch_result.value_tracker.average
742 erf_avg = erf_result.value_tracker.average
743 stretch_var = stretch_result.value_tracker.get_pessimistic_variance()
744 erf_var = erf_result.value_tracker.get_pessimistic_variance()
745 avg_log = (stretch_avg + erf_avg) / 2.0
746 var_log = (stretch_var + erf_var) / 2.0
747 var_log += (stretch_avg - erf_avg) * (stretch_avg - erf_avg) / 4.0
748 stdev_log = math.sqrt(var_log)
749 low, upp = math.exp(avg_log - stdev_log), math.exp(avg_log + stdev_log)
750 avg = (low + upp) / 2
752 trackers = (stretch_result.focus_tracker, erf_result.focus_tracker)
753 sea = math.exp(stretch_avg)
754 eea = math.exp(erf_avg)
755 return _ComputeResult(measurement, avg, stdev, sea, eea, trackers)
758 # Named tuples, for multiple local variables to be passed as return value.
759 _PartialResult = namedtuple(
760 "_PartialResult", "value_tracker focus_tracker samples"
762 """Two stat trackers and sample counter.
764 :param value_tracker: Tracker for the value (critical load) being integrated.
765 :param focus_tracker: Tracker for focusing integration inputs (sample points).
766 :param samples: How many samples were used for the computation.
767 :type value_tracker: stat_trackers.ScalarDualStatTracker
768 :type focus_tracker: stat_trackers.VectorStatTracker
772 _ComputeResult = namedtuple(
774 "measurement avg stdev stretch_exp_avg erf_exp_avg trackers",
776 """Measurement, 4 computation result values, pair of trackers.
778 :param measurement: The trial measurement result obtained during computation.
779 :param avg: Overall average of critical rate estimate.
780 :param stdev: Overall standard deviation of critical rate estimate.
781 :param stretch_exp_avg: Stretch fitting function estimate average exponentiated.
782 :param erf_exp_avg: Erf fitting function estimate average, exponentiated.
783 :param trackers: Pair of focus trackers to start next iteration with.
784 :type measurement: MeasurementResult
787 :type stretch_exp_avg: float
788 :type erf_exp_avg: float
789 :type trackers: 2-tuple of stat_trackers.VectorStatTracker