Modelling and predicting extreme behavior in critical real-time systems with advanced statistics

Modelling and Predicting Extreme Behavior in Critical Real-Time Systems with Advanced Statistics Sergi Vilardell Moreno A dissertation submitted in partial satisfaction of the requirements for the degree Doctor of Philosophy in Computer Architecture Supervisors: Dra. Isabel Serra & Dr. Jaume Abella Tutor: Dr. Francisco J. Cazorla January 25, 2023 Acknowledgements I wish to thank in general all people involved, directly or not, in my journey to complete my PhD. studies. While I am the main driver of this thesis, it has been a collaborative effort with a strong network of technical and emotional support. First and foremost, thanks to my first co-advisor Dr. Jaume Abella, and my tutor Dr. Francisco J. Cazorla, for the opportunity to develop my thesis under their excellent guidance. I would also like to thank Dr. Enrico Mezzetti for his contributions throughout these years. Finally, my most sincere gratitude towards my second co-advisor Dr. Isabel Serra for being the responsible of my development as a scientist and putting her trust in me from the very beginning by recommending me to pursue my thesis in Barcelona Supercomputing Center. I want to also thank Laurent Rioux for the opportunity to make an internship at Thales. This thesis would not be possible without the great work of the students and engineers of the Computer Architecture / Operating System group by paving the way towards incorporating statistical analysis in the Computer Architecture field. On the personal side, I wish to thank my family, Joan, Puri, and Enric, for raising me up to what I am today and always pushing me towards my real capabilities with tenderness and support. To my dear friends, the Brotherhood, for all those wonderful crucial years for our growth as people and scientists. In particular, my deep gratitude to Alejandro for being always there despite time and distance. Finally, to my dear love and main support during these years, Olaya, who demonstrated that not only technical ability is necessary to carry on a thesis, but that her love and emotional support has been the key to this journey. This work has been partially supported by Barcelona SuperComputing Center and the Spanish Ministry of Economy and Competitiveness (MINECO) under grant PID2019-110854RB-I00 / AEI / 10.13039/501100011033 and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 772773). i ii Abstract CriticalReal-Time Embedded Systems(CRTES)areusedindomainsliketransportation(e.g. avionics, automotive, space, and railway), healthcare, and industrial machinery. This subset of embedded systems requires undergoing a stringent Validation & Verification (V&V) process before they are allowed to enter in operation since any misbehavior can result in harm of humans or even fatalities. Software timing behavior is a key element to cover in the V&V process, providing with evidence that software runs timely. Software timing analysis, in turn, requires deriving bounds to each application task’s execution time. These bounds are referred to as Worst-Case Execution Time (WCET) estimates. As CRTES implement more complex safety-related functionalities in every new product, more complex software and consequently more performant computing hardware is used to satisfy the highperformance requirements. The side effect, however, of using more complex hardware and software is challenging state-of-the-art software timing analysis techniques. Measurement-Based Probabilistic Timing Analysis (MBPTA) techniques have been proposed to handle such hardware and software complexity providing tight and trustworthy WCET bounds (estimates). Specifically, Extreme Value Theory (EVT) has been used to provide with models for the most extreme occurrences in form of a probabilistic distribution. The output of the timing model is referred as a probabilistic WCET (pWCET). However trustworthy, EVT models can be cumbersome to apply and they sometimes can be exceedingly pessimistic which adds extra cost into timing budgets. This thesis investigates MBPTA techniques and develops novel methodologies within this framework in three distinct fronts. Firstly, by improving the tightness of pWCET models on sky-high quantiles with two models. A first one that combines risk analysis with EVT for a safe and accurate pWCET. And a second one that introduces Markov’s Inequality to the pWCET estimation problem, which provides with trustworthy guarantees with less requirements for its correct application. Secondly, in order to boost the use of data coming from performance monitoring counters - increasingly used by MBPTA techniques to tighten estimates-, this thesis shows two mathematically-based ways of merging multiple disjointed readings based on order statistics and copula models. Finally, this thesis proposes a model for the contention of competing tasks, when the timing profile obtained is limited, that allows to provide with more extreme WCET scenarios based on the dependencies between tasks. Summarizing, this thesis pushes the state-of-the-art forward in the V&V methodologies for CRTES in the framework of MBPTA in terms of WCET estimation and data gathering. iii iv Contents Acknowledgements i Abstract iii Contents v List of Figures ix List of Tables xi List of Acronyms xiii I Introduction, Background, and Experimental Methodology 1 1 Introduction 3 1.1 Certification with Higher Performance Requirements . . . . . . . . . . . . . . . . . . . 3 1.2 Measurement-Based Timing Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 1.3 HardwareEventMonitors................................... 4 1.4 MBPTA and Extreme Value Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.5 EVTforpWCETAnalysis ................................... 6 1.6 MBPTAwithLimitedData................................... 7 1.7 ThesisContributions ...................................... 8 1.7.1 Sky-high Quantile Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 1.7.2 Hardware Event Monitor Data Merging . . . . . . . . . . . . . . . . . . . . . . . 9 1.7.3 ContentionModelling ................................. 9 1.8 StructureoftheThesis ..................................... 9 1.9 ListofPublications ....................................... 10 2 Background 11 2.1 Measurement-Based Probabilistic Timing Analysis . . . . . . . . . . . . . . . . . . . . . 11 2.1.1 MBPTARequirements................................. 12 2.1.2 State-of-the-art on MBPTA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2 ExtremeValueTheory ..................................... 13 2.2.1 Law of Large Numbers: the Central Limit Theorem . . . . . . . . . . . . . . . . 13 2.2.2 Law of Extremes: Block Maxima . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.3 Threshold Law: Peaks Over Threshold . . . . . . . . . . . . . . . . . . . . . . . 15 2.2.4 State-of-the-art for EVT in MBPTA . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3 EVTAlternatives ........................................ 17 2.3.1 Chebyshev’s Inequality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.3.2 Markov’sInequality .................................. 17 2.3.3 State-of-the-art for EVT Alternatives for MBPTA . . . . . . . . . . . . . . . . . . 18 2.4 Hardware Event Monitors in MBTA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.4.1 Disproportion Between the Number of HEMs and PMCs . . . . . . . . . . . . . 18 2.4.2 State-of-the-art on HEM Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 2.5 Timing Validation in Automotive Systems . . . . . . . . . . . . . . . . . . . . . . . . . . 20 v Contents 2.5.1 State-of-the-art on MBPTA with Dependency . . . . . . . . . . . . . . . . . . . . 21 3 Experimental Methodology 23 3.1 BenchmarksandPlatforms................................... 23 3.1.1 Railway Case-Study on a LEON3+ Platform . . . . . . . . . . . . . . . . . . . . 24 3.1.2 Microbenchmarks on an NXP T2080 Platform . . . . . . . . . . . . . . . . . . . 24 3.1.3 Apollo Autonomous Driving Framework on an NVIDIA Jetson Platform . . . . 25 3.2 AnalyticalDistributions .................................... 27 3.3 Statistical Tests and Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.3.1 HypothesisTesting................................... 29 3.3.2 ConfidenceIntervals.................................. 30 3.3.3 Testing for Independence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.3.4 Test for Identical Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.4 Estimating the Extreme Value Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.4.1 HillEstimator...................................... 31 3.4.2 CoefficientofVariation................................. 32 II Sky-High Quantile Estimation for CRTES 35 4 Sky-high Quantile Estimation with Weibull Tails 37 4.1 Introduction........................................... 37 4.2 EVTlimits ............................................ 38 4.3 On the Use of Light Tails and Risk Analysis for WCET Estimation . . . . . . . . . . . . 39 4.3.1 Risk and Survivability Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 4.3.2 IHR Distributions in Survivability Analysis . . . . . . . . . . . . . . . . . . . . . 40 4.4 Equivalence Between IHR and Non-heavy Tails . . . . . . . . . . . . . . . . . . . . . . . 40 4.5 Weibull Tails (TailW) for pWCET Estimation . . . . . . . . . . . . . . . . . . . . . . . . 41 4.5.1 Formal Definition of tailW ............................... 42 4.6 FittingProtocol ......................................... 42 4.7 Evaluation ............................................ 43 4.7.1 Assessing Model Hypotheses: H-Convexity and Light Tails . . . . . . . . . . . 44 4.7.2 Assessment with Large Data Sets . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.7.3 Comparing exp,tailW and logc Models........................ 46 4.8 Summary............................................. 48 5 Sky-high Quantile Estimation with Markov’s Inequality 49 5.1 Introduction........................................... 49 5.2 Chebyshev and Markov Inequalities for pWCET Estimation . . . . . . . . . . . . . . . . 50 5.2.1 Markov’s Inequality on Low Probabilities . . . . . . . . . . . . . . . . . . . . . . 50 5.3 Power-of-kfunctions for Markov’s Inequality . . . . . . . . . . . . . . . . . . . . . . . . 51 5.3.1 Tightness of MIK for Increasing Values of k..................... 52 5.4 Handling Markov Sampling Uncertainty . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.4.1 Sample Moment Estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 5.4.2 Understanding the Behavior of max𝑘......................... 55 5.4.3 Deriving max𝑘from Unknown Distributions . . . . . . . . . . . . . . . . . . . . 56 5.5 RESTK and EVT PWCET Estimates on Distributions . . . . . . . . . . . . . . . . . . . . 57 5.6 RailwayUseCase........................................ 59 5.7 Summary............................................. 60 III Merging Hardware Event Monitor Data for CRTES 61 6 Merging Hardware Event Monitor Data for Complex MPSoCs Using Order Statistics 63 6.1 Introduction........................................... 63 6.2 Motivation............................................ 64 6.2.1 HEMVariability .................................... 64 vi Contents 6.2.2 Distribution....................................... 66 6.2.3 Reasons Behind the Observed Variability . . . . . . . . . . . . . . . . . . . . . . 66 6.3 ProblemFormalization..................................... 67 6.4 HRM: a Technique to Merge HEMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 6.4.1 Approach ........................................ 68 6.4.2 Procedure ........................................ 69 6.4.3 QuantileEstimation .................................. 70 6.4.4 CorrelationBoundary ................................. 71 6.4.5 Matrix Completion Techniques . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.5 ExperimentalEvaluation.................................... 73 6.5.1 Validation Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 6.5.2 Independence and Identical Distribution . . . . . . . . . . . . . . . . . . . . . . 73 6.5.3 Correlation Between HEMs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 6.5.4 Overheads........................................ 75 6.6 Summary............................................. 76 7 Merging Hardware Event Monitor Data for Complex MPSoCs Using Copulas 77 7.1 Introduction........................................... 77 7.2 MUCH: Multi-Correlation HEM Reading and Merging . . . . . . . . . . . . . . . . . . 78 7.2.1 HRMAnalysis ..................................... 78 7.2.2 MathematicalApproach................................ 79 7.2.3 Procedure ........................................ 81 7.3 Evaluation ............................................ 82 7.3.1 ValidationApproach.................................. 82 7.3.2 Results.......................................... 83 7.4 Summary............................................. 85 IV Contention Modelling for CRTES 87 8 Clean Execution Times for MBPTA 89 8.1 Introduction........................................... 89 8.2 Timing Analysis Validation on Apollo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 8.2.1 Prediction Module Timing Behavior . . . . . . . . . . . . . . . . . . . . . . . . . 90 8.2.2 Aggregation of Overlaps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 8.3 CleanET ............................................. 92 8.3.1 SettingtheObjective.................................. 92 8.3.2 ModellingApproach.................................. 92 8.3.3 Mathematical Development of CleanET . . . . . . . . . . . . . . . . . . . . . . . 93 8.3.4 Obtaining r....................................... 94 8.3.5 IntrinsicPessimism................................... 95 8.3.6 CleanET for Multiple Overlappings . . . . . . . . . . . . . . . . . . . . . . . . . 96 8.4 Evaluation ............................................ 97 8.4.1 Dilation Factors (𝑟𝑖)................................... 97 8.4.2 Validation of the Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 8.4.3 UsingCleanETResults.................................100 8.5 Summary.............................................101 V Conclusions and Future Work 103 9 Conclusions and Future Work 105 9.1 Conclusions ...........................................105 9.2 FutureWork...........................................106 Bibliography 109 vii List of Acronyms r.v. random variable. 11 STA Static Timing Analysis. 4 V&V Validation and Verification. 3 WCET Worst-Case Execution Time. 3 xiv Part I Introduction, Background, and Experimental Methodology 1 Chapter 1 Introduction Historically, computer systems have been used mostly for automation and entertainment. This trend has changed in recent years with computers being increasingly used in a number of critical domains of our life including transportation and health, among others. In these domains, computers are embedded into devices like a car, a train, an airplane, and medical machinery. Unlike other Embedded Systems (ES) in which a malfunctioning of the computer causes mostly discomfort, a malfunction in the so-called Critical Embedded Systems (CES) can cause serious harm to the environment, or people, even resulting in fatalities. In order to ensure that these systems are safe to use in critical scenarios, it is imperative that their design satisfies functional correctness [61] tested according to a well-defined Validation and Verification (V&V) process. Functional correctness ensures that the system functions as intended without anomalies. On the other hand, timing correctness ensures that the performed tasks complete their execution within the assigned time budget or deadline. In this thesis, we focus on CES whose correctness also depends on their timely response, which are categorized as Critical Real-Time Embedded Systems (CRTES). In the domain of CRTES, the functional and timing correctness have equally important weight through the V&V process, given that not meeting the pre-assigned time budget, or deadline, can lead to catastrophic consequences. Thus, in timing analysis deriving a conservative estimation of the Worst-Case Execution Time (WCET) is essential. The interest in keeping up with the challenges of timing analysis is due to the increasing growth in CRTES in services like healthcare and industries like avionics, railway, and automotive vehicles [72]. Furthermore, the systems are growing exponentially in complexity: the car manufacturer Volvo stated that in 2020 their cars needed around 100 million lines of code to implement the key functionalities of a car such as electronic fuel injection, transmission control, navigation systems, etc. Their measurements show that since 2000, the automotive software increases an order of magnitude every decade [9]. Their functionalities are also way more complex too. In the case of autonomous cars, the system must handle data coming from radars and LIDARS and process it using state-of-the-art Neural Networks algorithms in real-time, requiring very high levels of performance. To satisfy those needs, CRTES need to rely on high-performance hardware such as multi-level cache memories design and multi-core architectures. While powerful, this kind of design tends to be harder to analyze which makes them more challenging to work with towards an accurate WCET estimation. 1.1 Certification with Higher Performance Requirements In past years, the prevailing architecture paradigm for the design of CRTES was the federated architecture. These architectures comprise multiple subsystems, each with critical functionalities, which are physically separated from each other. The increasing complexity of new functionalities in recent years challenge the scalability of federated systems [96]. The higher performance needs of today’s functionalities would rise the cost, weight, and size of CRTES. In that regard, the CRTES industry shifted towards integrated architecture paradigms that integrate several functionalities in the same computing hardware. Integrated architectures can also naturally leverage the parallelism of multicore systems and the capabilities of high-performance features like multi-level cache. To that end, software solutions like resource partitioning are used to contain the impact that software partitions can generate on each other. While the adoption of integrated architectures has aided the cost-weight-size 3 Chapter 1. Introduction problem and provided with satisfactory performance, the certification on multicore systems becomes notably more challenging. The loss of guarantees in terms of independence between subsystems within integrated architectures affects their timing predictability. Depending on the degree of integration and partitioning the dependency, and therefore predictability, will be affected differently. This usually results in a trade-off between performance and degree of integration. In terms of certification for CRTES, new designs trend towards favoring timing predictability [140] instead of increased average performance. 1.2 Measurement-Based Timing Analysis Increasingly complex software functionalities in CRTES require high-performance multicore processors (computing hardware). However, more complex processors have drawbacks in relation to the V&V process. Complex CRTES should comply to the same safety standards as the simpler architectures and given the interconnected nature of integrated architecture and the complexity of multicores, the V&V process for timing is very challenging in terms of effort and accuracy [5]. A highquality WCET estimation in this domain must be affordable in terms of costs and trustworthiness, so that the produced execution time estimate tightly upper-bounds the real execution time. Timing analysis techniques can be broadly categorized into two groups: Static Timing Analysis (STA) and Measurement-Based Timing Analysis (MBTA) [172]. STA does not need to execute code on real hardware and instead it relies on creating abstract models of the architecture and analyzing possible execution paths for the task under analysis. This approach faces limitations due to the limited information about the hardware provided by the hardware manufacturer and its accuracy [1]. In fact, some processor manufacturers do not disclose full information about their systems’ microarchitecture in order to protect intellectual property. Even assuming the validity and completeness of the information disclosed, the technical manual of a modern Multi-Processor System on Chip (MPSoC) like the Zynq UltraScale+ is more than a thousand pages long [179]. The methods that worked for simpler single core processor architectures do not scale well with the complexity that CRTES have nowadays. In less complex architectures, the approaches for V&V were usually guided by previous engineering experience [41]. MBTA considers a different approach, in which timing measurements of a given task in a given hardware guide the analysis. The CRTES system engineer is responsible for designing experiments with stressful initial conditions for the task running in real hardware or simulator. Those experiments should reflect the timing behavior at operation so that the bounds derived from it are reliable. However, measurements in themselves are not treated as bounds [173], given the general impossibility of reproducing and forcing those execution conditions that lead to the worst-case scenario. Industry standard solutions to derive bounds from measurements rely on adding safety margins. For instance, on top of the high-water mark, a safety margin (e.g. 20%) can be added. The logic behind it is to compensate the unobserved effects that can impact observed values. However, for complex hardware deriving a safety margin is significantly more complex. 1.3 Hardware Event Monitors Theexecutiontime ofaprogramina complexcomputingplatformis theproductof manyinteractions happening within the different hardware blocks. In complex CRTES these interactions are hard to control and model. For instance, the functionality of multi-level cache memories is deeply connected to the variability of the program’s execution time. Smaller level 1 caches directly connected to the core have low latency and their data can be quickly accessed. Due to cost and convenience, caches at a higher-level (e.g. L2) are bigger in size, more distant from the core, and the time required to access them increases. A program that tends to perform accesses to the higher-level caches, or the off-chip memory, will suffer higher latency to fetch data and instructions than the one which mostly operates at low-level caches. Many of these events describing the usage of hardware computing resources made by a task can be tracked via Hardware Event Monitors (HEMs) which can be hence used for software timing analysis purposes. In fact, several MBTA techniques build on HEMs to model and control contention for 4 1.4. MBPTA and Extreme Value Theory instance by limiting the maximum number of accesses a task is allowed to do to a given resource [32, 125]. HEMs are programmed and accessed by the Performance Monitoring Unit (PMU) in the processor, which is also referred to as statistics unit. The PMU offers a set of visible registers called Performance Monitoring Counters (PMCs) which can be configured to read the number of occurrences of a given HEM. 1.4 MBPTA and Extreme Value Theory In an effort towards relieving some of the burden on the application of MBTA to complex systems, Measurement-Based Probabilistic Timing Analysis (MBPTA) approaches have been gaining traction [3, 21, 23, 44] in CRTES over more than a decade ago. The variety of execution conditions complex CRTES can produce causes huge differences between runs of the same task and great variance in the distribution of execution times. Hence, instead of aiming for a fixed and likely pessimistic WCET, a better alternative is to provide a probabilistic view into timing analysis [36]. Instead of fixing a single value that upper-bounds all executions, MBPTA provides a distribution that upper-bounds the probability of exceeding any specific execution time, although we are interested in high values. For each arbitrarily rare occurrence, an MBPTA tool provides a probability of exceeding it, thus a probabilistic Worst-Case Execution Time (pWCET) estimate. That is, on MBPTA the results are given in the form of a distribution. When dealing with the distribution of arbitrarily rare occurrences, Extreme Value Theory (EVT) is the theoretical best approach [40]. In particular, the appeal of extreme value analysis is to estimate the probabilities of events which are rarer than the ones observed. With EVT one can assess for instance, the likelihood of an earthquake bigger than the ones we have observed so far. Given the generality of its results, EVT has been applied in multiple fields. EVT is still present in fields which have very little in common between them like traffic safety [79], finance [27], earthquake analysis [31], biophysical science [120] and even musical analysis [146]. The different nature of each application field requires to adjust the approach to each problem instance, but the underlying laws are still the same. Of course, for instance, in case of earthquakes one cannot generate more data so it makes sense that EVT is used in this domain. Nonetheless, EVT is also applied in fields where data is readily available, but the cost or even the chance to obtain it is barely possible or feasible. In the case of WCET estimation, the execution time of a program may be relatively short, but obtaining extreme execution times may take weeks or even months. Therefore, one must resort to construct a pWCET model with smaller yet representative set of data. And it is in the word representative where a big part of the challenges arises in MBPTA. Obtaining stressful benchmarks for the CRTES at analysis that generates the most extreme execution conditions in the system is a very challenging task. In my view, there are two sides to the MBPTA problem. The first is the technical aspect of computer science that is concerned about designing hardware that is compliant with MBPTA experimentation, and working towards obtaining representative extreme conditions [18, 117]. The second one, the focus of this thesis, is the modelling aspect of the estimation of unseen extreme execution times under the upper-bounding condition. There have been several mentions to extreme values without proper definition. In EVT there are two classical ways of thinking about extremes. Since this thesis is concerned about WCET estimation, we consider extremes as the high-end values of our measurements, although EVT also work for the lowest values. •The first way is considering that the extremes are only the maximum value on your sample. Imagine we have a sample of size 𝑛and we divide it into separate groups of size 𝑆. Then, for each group, we obtain the maximum value on it. This is called the Block Maxima (BM) approach. •Thesecond wayis consideringextremes asthegroupof values aboveacertain threshold. Thisis the Peaks over Threshold (PoT) approach. To illustrate them, Figure 1.1a and 1.1b represent the set of values deemed as extremes using the BM and PoT approach respectively. In this example, the sample of data is the same, yet the set of values deemed as extreme is different. There are values classified as extremes in Figure 1.1b which are classified as belonging to the body of the 5 Chapter 1. Introduction Body 1 3 (a) Block Maxima 1 3 Body (b) Peaks Over Threshold Figure 1.1: Both classical EVT ways of considering extremes. distribution in Figure 1.1a. Not only that, but one can also see how changing the size of the block 𝑆or the location of the threshold will change the group of extreme values. This is not inconsistentwithinEVTand both approachesyieldgoodmodelsfor extremes. Yet, theselection of the block size or the threshold is a well-known problem with many different approaches as we will see in Chapter 3. The key finding, and the core of EVT, is that these extreme values have a distribution of their own. The resulting extremes from the BM approach tend to a family of distributions called the Generalized Extreme Value (GEV) distribution [65], while the extremes from PoT converge to another family of distributions, named the Generalized Pareto Distribution (GPD) [15]. These results, under certain conditions, are general to any random variable; which allows us to model rainfall annual data [91] and execution times from CRTES under the same theoretical umbrella. These theoretical family of distributions are named so because there are different distributions within them depending on the Extreme Value Index (EVI) parameter represented as 𝜉. The EVI is what characterizes the shape of the most extreme values of a distribution. For the PoT model the EVI determines the shape of the tail of the distribution in three ways. We refer to the tail as the probability mass at the upper end point of the distribution. Figure 1.2 shows these shapes for different distributions and values of the EVI. If the EVI is 0, the shape of the tail will be exponential, meaning that the probability of greater values decays as an exponential function. If the EVI is negative, the probability of the extremes will decay faster than the exponential. This kind of tail is named a light tail. Otherwise, positive values of the EVI will make the probability decay slower than the exponential, which makes for a heavy tail. The methods to estimate the EVI used in this thesis will be explained in Chapter 2. 1.5 EVT for pWCET Analysis The integration of EVT in MBPTA has received increasing attention in the last decade [44, 82]. The research community working on MBPTA has been debating on which kind of tails are more appropriate for pWCET. It has been shown by argument and empirically that exponential tails are a safe upper-bound in general [3, 149]. Nonetheless, exponential tails can be too pessimistic to be usable, especially for very high quantiles. The pessimism is in part due to the very nature of the distribution. The exponential tail assumes that the execution time may grow to infinity. However, in CRTES this not the case. In the design process of CRTES, each task has an assigned deadline. Therefore, the estimation of a WCET is of utmost importance to ensure the task is executed within its time budget. Then, the theoretically best bound would be one that takes into account the finiteness 6 1.6. MBPTA with Limited Data ξ Weibull ξ Fréchet ξ Gumbel ξ Beta ξ Gamma Weibull Gamma Beta Gumbel Fréchet Figure 1.2: Shape of distribution tails with different EVI 𝜉. of execution times in CRTES programs. Light tails in EVT fulfill this purpose as their domain is bounded with a maximum value and decay faster than the exponential. This is however, both a blessing and a curse. When using a light tail model to fit our execution time samples, it will take the maximum of those samples as the reference for the bounding. If in our samples we do not get an execution time close to the WCET, the light tail model will underestimate eventually. In reality, while light tails models tend to be very accurate in practice, in order to overestimate the WCET, they need to be fed with values in the sample that are very close to the real WCET, which is the very same problem we are trying to solve. In summary, while classical EVT and its two fundamental theorems offer the best theoretical solutions to the pWCET problem, neither the GEV or the GPD cannot prevent the uncertainty in estimating sky-high quantiles which are outside of the observed sample, thus in general estimating a tight and trustworthy upper-bound for the WCET with EVT is challenging. In this thesis, we identify some of these limitations and implement solutions that combine EVT with other statistical tools, and also find new tools outside EVT to estimate extreme values. 1.6 MBPTA with Limited Data Frequently, CRTESengineerscannotaffordanexhaustivesetofbenchmarksortestvectorstostressthe system towards the WCET scenario. That is the case of timing analysis for Autonomous Driving (AD) systems, where the demonstrations via driving test serve as the benchmarks for providing guarantees over the system behavior at operation. However, designing benchmarks that are representative of real-life driving cases is a very challenging task [80]. As mentioned before, the difficulty in assessing WCET scenarios is the complexity of the hardware and software. In the case of AD the software complexity is at the core of its functionality, with Machine Learning algorithms deciphering what the perceptive modules receive in real-time. In these systems, there can be different tasks accessing concurrently the same hardware resources on the system, which causes timing interference or contention impact. In this thesis, we take another point of view in the TA of these systems. Instead of trying to provide with a representative set of benchmarks, we model the intrinsic dependency that occurs within these systems. We aim at modelling the time dilation that occurs when similar tasks are running at the same time. A model of this kind can give us an estimation of a WCET scenario that may occur in very stressful situations with extreme competition among tasks. Ultimately, this is expected to lessen the burden of the benchmark selection for the V&V. 7 Chapter 1. Introduction 1.7 Thesis Contributions The contributions of this thesis are towards the design of new and improved methods for the timing analysis of software running on complex CRTES. It is a great challenge to improve on the classical fundamental theorems of EVT because they provide with the theoretical best solution for predicting extreme values. However, MBPTA’s characterizing trait is the upper-bounding constraint, which is difficult to achieve both in theory and practice. Also, in CRTES there is a lot of data coming from the monitoring support in the hardware, which would allow creating more sophisticated MB(P)TA techniques accounting for varying contention scenarios. Such data is typically not used owing to inherent hardware limitations, constraining to sampling only a restricted number of variables at once. We provide solutions that allow merging monitoring data from multiple separate experiments as if they were extracted from the same execution conditions. Last but not least, apart from (p)WCET estimation for system design, there is a need for validation (testing) campaigns to assess that the system is properly designed and WCET bounds are not exceeded. However, the lack of controllability of how tasks overlap in multicores precludes from having sufficiently exhaustive test campaigns. In summary, we can divide the contributions in three blocks: i) improving the current MBPTA methods in tightness and safety for the pWCET estimation problem, ii) creating new methods to analyze the monitoring data of CRTES by putting it in an orderly way that can be used by MBPTA tools, and iii) devising new methods to produce execution time observations for tasks overlapping in arbitrary (and worst-case) ways that test (validation) campaigns fail to produce. More specifically, the contributions are divided in three parts: 1. Sky-high quantile estimation techniques. 2. Hardware Event Monitor data merging techniques. 3. Contention modelling techniques. 1.7.1 Sky-high Quantile Estimation The estimation of very rare occurrences for the WCET is still an ongoing problem. Those extreme values, a.k.a. sky-high quantiles, are typically in the exceedance probability range of 𝑝={1−10−6,1− 10−12}. In order to fulfill MBPTA requirements, the pWCET estimation should seek to upper-bound which is an added difficulty on the estimation of extreme values. In that regard, in this thesis we maketwocontributionstowardsimprovingtheWCETestimationproblemusingmeasurement-based techniques which aim at: i) Improving the accuracy of the estimation of sky-high quantiles, and ii) keep providing with the same safe guarantees of the estimation while needing less requirements. 1. On our first contribution towards sky-high quantile estimation, we develop a new parametric model with analytical expression which can produce safe and accurate estimates while overcoming some limitations of the usual EVT approaches. This technique was devised by making use of both EVT and Risk Analysis, which helps us to make safe guarantees. The proposed parametric model for tail estimation is the tail of a Weibull distribution, or as we named it tailW. The solution is tested against real case-study hardware data from a CRTES and improves on the current best EVT solution. This methodology is publicly available in an R package on the official repository (CRAN). 2. Our second contribution towards the pWCET estimaion problem is making use of Markov’s Inequality to upper-bound sky-high quantiles. Unlike EVT, and even tailW, Markov’s Inequality does not try to fit the data but rather upper-bounds data using simple probabilistic tools. We modify the simplest version of Markov’s Inequality with the power-of-kfunction, which provides with much tighter upper-bounds. This works theoretically because one needs information about the distribution of the sample at analysis. The challenge of this problem is in working out a practical solution that allows to take an unknown distribution and still being able to provide trustworthy and tight bounds. We validate this method against reference theoretical distributions and real case-study hardware data from a CRTES, which improves on the current best EVT solution. 8 1.8. Structure of the Thesis 1.7.2 Hardware Event Monitor Data Merging Multicore processors do not allow to read all the HEMs at once. HEMs which provide information about the system at operation, can be read in small groups. Multicore processors have only few (4-6) PMCs while the number of HEMs is in the order of hundreds [67, 68, 69]. Therefore, not much information from the hardware can be extracted at once at operation under the same execution conditions. To make things worse, conditions greatly vary from run to run, making the comparison of runs a relevant concern. This presents a challenge for the use of HEMS in MBTA since, while information from the system at operation can be obtained, only a small subset of the HEMs can be measured at the same time. In this thesis we explore methods to capture as much data as possible that can be used to the advantage of more robust V&V techniques. •The first contribution makes use of statistics of order to merge HEMs. We do so by selecting a common HEM that will be present in all collections of HEMs to read on a given run. With the use of this common HEM, called the anchor HEM, we can proceed to perform a merge based on the order statistics of the anchor. We evaluate if the HEMs keep the relationship they have with one another after the merge. •The second contribution makes use of Copula Theory in order to preserve the relationship between HEMs after the merge. This methodology uses a more data expensive procedure to merge HEMs because one needs to measure the relationship between all HEMs before reading. Then, modelling said relationships with a Multivariate Gaussian copula, we can perform the merge and obtain a final dataset were all values are sorted on similar execution conditions. In both methods we use the T2080 platform for validation purposes. 1.7.3 Contention Modelling High-performance processors used in CRTES provide the required computing capacity, but this comes at the cost of additional challenges to the timing analysis of the systems. The fact that multiple cores access shared resources at the same time can produce delay in the execution time of tasks. Therefore, when multiple tasks are executed at the same time, they tend to create contention on each other. In this final contribution, we design a model for inter-task contention that occurs in multicore systems. At system operation we have the timeline of the execution times of different tasks, although they are rarely executed in isolation from other competing tasks. In this work, we show how the contention model we constructed allows us to estimate isolated (solo) execution time and WCET scenarios in systems were this information cannot be obtained. Here we work on a mathematical model of contention to estimate the general impact of contending tasks in CRTES. This model allows us to estimate less pessimistic WCET conditions when working with limited information. We validated this model using execution times from a demo of the Apollo autonomous driving software. 1.8 Structure of the Thesis •In Chapter 2 we introduce some background to understand each of the contributions made on this thesis. First, a deeper understanding of the limitations of EVT for the pWCET estimation problembased ontheoretical argumentsispresented. Then, we explainindetailthe monitoring data coming from CRTES that cannot be extracted at once. Finally, we explain in detail the contention suffered in MPSoCs and the need to model it for the pWCET problem. •Chapter 3 shows all the practical and theoretical tools that were used in this thesis which are divided in several blocks. The first block is the statistical one. The final one shows the benchmarks used to validate our solutions. •Chapter 4 is devoted to the first solution to the WCET problem, the usage of Weibull tails to estimate sky-high quantiles. In this chapter we combine EVT with Risk analysis to provide a tight yet trustworthy solution to the WCET problem. •In Chapter 5 another solution for the WCET problem is shown. We explain how Markov’s Inequality, which upper-bounds by construction, can be the foundation for tight and trustworthy 9 Chapter 2. Background Peaks Over Threshold (PoT). Given a r.v. 𝑋with a CDF, F, and a threshold, 𝑢>0, such that 𝑦=𝑥−𝑢, the excess r.v., 𝑋𝑢, defined as (𝑋−𝑢|𝑋>𝑢)is given by the CDF, 𝐹𝑢, defined as: 𝐹𝑢(𝑦)=𝑃(𝑋−𝑢≤𝑦|𝑋>𝑢)=𝐹(𝑢+𝑦)−𝐹(𝑢) 1−𝐹(𝑢), 𝑦 ≥0.(2.12) Then 𝐹𝑢is the excess distribution function, which is the distribution of the excess value over a threshold 𝑢. A PoT model is a semi-parametric model where the law of 𝑋for 𝑥<𝑢is described by the empirical distribution, and for 𝑥>𝑢is defined by a parametric model for the tail. It is semi-parametric because, if we want to apply EVT models, we are concerned with the tail of the distribution. Therefore, for the values of the r.v. higher than the threshold 𝑋>𝑢we use a parametric model 𝑆𝑢(𝑥𝑢), with 𝑥𝑢=𝑥−𝑢, and for the cumulative probability up to the threshold 𝑢we use the empirical probability mass given by the sample. More formally, 𝑃(𝑋>𝑥)=𝑆(𝑥)=𝑆(𝑢)𝑆𝑢(𝑥𝑢),(2.13) where 𝑆(𝑥)=1−𝐹(𝑥)and 𝑆𝑢(𝑥𝑢)=1−𝐹𝑢(𝑥𝑢)are the CCDFs. Theorem 3 (Pickands–Balkema–de Haan theorem [15]).Let 𝐹be a distribution function such that 𝐹(𝑥)<1with 𝐹𝑢being its conditional excess distribution function. Then, 𝐹𝑢converges in probability to the Generalized Pareto Distribution (GPD) for large 𝑢. That is, 𝐹𝑢ℒ −→ 𝐺(𝑦;𝜎,𝜉)as 𝑢→ ∞, where: 𝐺(𝑦;𝜎,𝜉)=(1−1+𝜉𝑦−𝜇 𝜎−1 𝜉if 𝜉≠0, 1−exp−𝑦 𝜎if 𝜉=0.(2.14) With 𝜎 > 0, and 𝑦≥0when 𝜉≥0and 𝜇≤𝑦≤𝜇−𝜎/𝜉when 𝜉 < 0. This result is crucial for tail estimation. If the conditions for the theorem are met, the GPD family of functions results in accurate estimates for the most extreme values on the tail. However, the conditions are met when the threshold 𝑢→ ∞, which brings uncertainty in the implementation of estimates with the GPD since finite values of 𝑢need to be used. 2.2.4 State-of-the-art for EVT in MBPTA Extreme Value Theory, a consolidated approach for modelling and predicting the occurrence of rare events, has emerged as the preferred option for modelling the WCET of software programs in CRTES. EVT has been considered particularly fit for probabilistic modelling of WCET as the latter is normally considered a rare event in the program’s timing behavior. EVT is at the foundation of several MBPTA approaches [3, 20, 35, 82, 149, 165], which have been already positively assessed in some industrial use cases [170]. Probabilistic approaches building on sample and population sizes have also been built for overlapping concerns across task scheduling and timing analysis [127], hence being orthogonal to WCET estimation. Several works assess the necessary conditions for a correct application of EVT to the timing problem since an inattentive application of the EVT statistical tools can severely affect both trustworthiness and quality of the derived pWCET bounds [71, 105, 118]. Several studies focus on EVT applicability preconditions on the (timing) observations being i.i.d. r.v.s [35, 44]. EVT has also been shown to be applicable also to stationary data preserving extremal independence [142]. Different statistical tests have been assessed for that purpose in the real-time literature [3, 35, 44, 137]. Also, platform randomization [97] or data sample randomization [104] have been used to meet the i.i.d. requirement. However, even in case such statistical preconditions are met, the reliability of the obtained pWCET bounds is still affected by the choice in the EVT inputs (i.e. selection of samples belonging to the tail) and parameters of the fit distribution. Several methods have been proposed for selecting and assessing the quality of EVT parameters and sample selection and, in turn, their impact on the trustworthiness of the computed bounds [3, 12, 138]. Several authors have considered general EVT distributions, without restricting them to the exponential law [76, 104, 105], whereas others build on the particular characteristics of the problem 16 2.3. EVT Alternatives modelled – finiteness of the execution time of critical real-time programs – to fit exponential tails only [3, 44, 59, 82], which have been shown to provide usable bounds and confidence intervals [149]. Probabilistic and statistical approaches [35, 48] have been increasingly considered as a promising solution to cope with the rise in complexity of hardware and software systems, as determined by unprecedented increasing computing performance requirements [54]. 2.3 EVT Alternatives In this thesis we propose alternative models for estimating very rare occurrences in the context of MBPTA. Our aim is to propose new methodologies that maintain safety guarantees while needing less requirements of application. 2.3.1 Chebyshev’s Inequality Additionally, the other front is to mitigate EVT’s model uncertainty. Given that EVT has uncertainty in the modelling process that cannot be avoided by construction due to the threshold estimation, we want to resort to methods which can produce safe upper-bounds which do not have uncertainties embedded into them. In that regard, there is a probabilistic inequality called Chebyshev’s Inequality which upper-bounds the cumulative probability of exceeding a given value of the r.v. at study. This inequality upper-bounds by construction any positive r.v. which makes for a suitable candidate for starting a pWCET methodology. Let us introduce this probabilistic tool. Theorem 4 (Chebyshev’s Inequality [156]).Let 𝑋be a non-negative r.v., 𝑏>0, and 𝑓a non-negative and increasing function. Chebyshev’s inequality states that: 𝑃(𝑋≥𝑏) ≤ 𝐸(𝑓(𝑋)) 𝑓(𝑏).(2.15) For WCET estimation, 𝑋corresponds to the execution time distribution to be bounded, and 𝑏to an execution time for which we want to find its upper-bound probability. Regarding 𝑓, it needs to be defined to realize the general Chebyshev’s inequality into a specific upper-bound function. The function 𝑓can be any non-negative function so that for a particular domain D, Property 2.16 holds. ∀𝑥∈𝐷, 𝑓 (𝑥)>0,(2.16) 𝑓must also be an increasing function so that for a given interval 𝐼, Property 2.17 holds. 𝑎, 𝑏 ∈𝐼|𝑎<𝑏, ⇒𝑓(𝑎) ≤ 𝑓(𝑏).(2.17) Interestingly, Chebyshev’s inequality does not require determining where the actual tail distribution starts but instead works with the entire distribution, hence removing EVT’s model uncertainty for tail selection. In fact, Chebyshev’s inequality carries no model uncertainty. Observation 1. Chebyshev’s inequality is a model uncertainty free general model for pWCET estimation. Chebyshev’s inequality also applies to continuous and discrete distributions regardless of their characteristics (e.g. shape, variance, kurtosis, etc.). Also, it is non-parametric, i.e. it makes no assumption on parameters for the studied distribution. 2.3.2 Markov’s Inequality Markov’s inequality is a specific instantiation of Chebyshev’s inequality. Corollary 1 (Markov’s Inequality [111]).Let 𝑋>0and let the function 𝑓be the identity function 𝑓(𝑋)=𝑋. Hence, Markov’s inequality yields: 𝑃(𝑋≥𝑏) ≤ 𝐸(𝑋) 𝑏.(2.18) 17 Chapter 2. Background As for the baseline Chebyshev’s inequality, Markov’s inequality holds for any real-valued r.v. with a finite expected value and positive value 𝑏. Also, it (i) is a trustworthy upper-bound, by construction, of the underlying distribution; (ii) has no model uncertainty; (iii) is non-parametric; and (iv) can be applied to discrete and continuous distributions. In Chapter 5, we will assess Chebyshev’s Inequality viability for an MBPTA methodology. 2.3.3 State-of-the-art for EVT Alternatives for MBPTA Markov’s and Chebyshev’s inequalities have been historically applied in a wide variety of fields such as Engineering [160], Big Data sampling [143], and radiative transfer [163]. In the context of real-time computing, moment-based bounds on tail probabilities have also been considered in the scope of probabilistic schedulability analysis [37, 38, 56, 167]. Chernoff bounds [37, 38] and generalizations thereof [167] are exploited to compute the cumulative distribution of the interference caused by higher-priority tasks on a task response time, ultimately delivering a probability for deadline misses. While building on application of Markov’s inequality to the timing dimension, these approaches do not address the problem of computing probabilistic bounds to tasks’ execution time, which are instead assumed to be available, but offer a scalable alternative to the computational complexity of convolutions. Some works [106, 135] consider the use of Chebyshev’s inequality for WCET and/or cache hit and miss rates estimation. However, those works consider Chebyshev’s inequality only to estimate the impact of the variance on those metrics, and focus on the analysis of statistical uncertainties. Other works consider using higher moments to improve concentration inequalities similar to Markov [101]. Finally, authors in [185] consider a similar approach to those works for WCET estimation, but discard Chebyshev’s inequality altogether given the pessimism expected for high quantiles. 2.4 Hardware Event Monitors in MBTA Powerfulmeasurement-based timing analysissolutionshavebeen proposedformulticores forcritical applications building on hardware and software profiling [24, 57, 144]. In this context, the PMU in MPSoCs offers information relevant for timing analysis [73, 115], enabling V&V and software time budgeting of time-critical applications. For instance, the usage of some shared resources can be budgeted, monitored, and enforced using event quotas building on HEMs reached through the PMU [46, 126, 131, 183]. HEMs are used to monitor when tasks exceed their usage quotas which are suspended. Furthermore, HEM data has been shown to aid sampling representativeness [33]. Indeed, HEM information isalready usedas apillar tocertify criticalavionics systems[134], so PMUs, and the HEMs they allow monitoring, become the basis of industrial-quality multicore interference mitigation and estimation techniques. However, HEM information is not straightforward to use. Due to the impossibility of controlling execution conditions in complex CRTES. A program run several times and performing the same number of instructions, thus maintaining functionality, will yield different HEM results. In this section we analyze the variability of HEM data within the context of the NXP T2080 Reference Board. We show that several of the 262 HEMs in the T2080 present significant variation (Section 6.2.1) and follow different distributions (Section 6.2.2). We also dig down into some of the reasons behind the observed variation (Section 6.2.3). 2.4.1 Disproportion Between the Number of HEMs and PMCs Despite the increase in number and specialization of HEMs, their observability and accessibility in reference CRTES platforms is typically constrained by the availability of a relatively (but consistently) lower number of PMCs. The latter represent, in fact, the most natural way of making HEM information available to the user. The number of PMCs available in modern MPSoCs typically ranges between 4 and 8 per core, which inherently clashes with the number of HEMs an analysis would need to track. In the case of ARM A53/A57/A71 cores, the number of PMCs is 6, which also matches the number of PMCs in NXP cores e500mc/e6500. 18 2.4. Hardware Event Monitors in MBTA Overall, the clear imbalance between the large number of HEMs and available PMCs calls for approaches that reliably merge HEM readings across different experiments. It is also worth mentioning that MPSoCs include an increasing number of complex shared resources. This will naturally result into more HEMs tracked by timing analysis techniques to capture the effect of contention. Hence, several runs are required to read all the HEMs of interest, which are later ‘merged’ off-line to analyze the program behavior and reason about contention. For instance, to decide whether some tasks can be scheduled concurrently, we need to budget how much each one is expected to access each shared resource, which requires consistent reads of a large number of HEMs. To make things worse, several runs of the same experiment in an MPSoC can result in inevitable variations in the timing behavior of the program, though its functional behavior is the same. This is due to the impossibility to control the entire hardware and software initial state in each run. In practical terms, this translates into variability in HEM readings (ashigh as 59% for processor cycles in our target system for relevant HEMs), with no variability observed in instruction count (as analyzed in Section 6.2.3). The engineer is confronted with a set of values (readings) for each HEM, that need to be merged to allow reasoning about multicore contention. Unfortunately, since HEM values from different runs to be merged can be subject to different (large) noise, it is challenging to merge them consistently so that merged HEM vectors – those where all HEMs of interest are included – resemble the values that would have been obtained if they could have been read all of them simultaneously in the same run. 2.4.2 State-of-the-art on HEM Analysis In the CRTES domain, several works build on HEMs for the estimation of bounds to software timing. Paulisch et al. [126] create an analysis and runtime monitoring solution for limiting task contention in multicores by tracking and controlling HEM. In the same vein, Diaz et al. [55] build on HEM to produce an ILP-based contention model for an AURIX automotive microcontroller. Likewise, Santinelli et al. [77], build on the HEM of a multicore system to derive probabilistic WCET estimates. Griffin et al. [75] derive a method to select the HEM with highest contribution to software timing and predict execution time under unseen configurations. More recently, information fromHEMs hasbeen exploitedas the cornerstone ofindustrial-qualityap- proaches [134, 162] for providing the necessary evidence for supporting the certification of multicore CRTES, in conformance with the requirements from domain-specific certification authorities [62]. Several works in the mainstream (high-performance) domain reason on the sources of variability in HEM values when executing several times the same piece of software. This covers from the operating system noise [119], application variability [7, 121] and the particular HEM-Reading library, to the complexity of the hardware [184]. For instance, [119] focuses on the cycle count HEM and shows that its variability is often related to the executable layout and operating system issues. Also, at software level, [184] assesses the accuracy of various high-level counter APIs with focus on cycle count and total retired instruction HEMs. In our work, we use no operating system and access directly, with no library, the HEMs (via the PMCs) so they are not subject to software-induced variability. In [121], authors focus on task-parallel programs in high-performance environments with highly-dynamic execution conditions, including dynamic task scheduling, that cause tasks to execute in different orders and in different cores across executions. Authors propose techniques to determine which HEM readings belong to each task and hence, combine them to derive all HEMs for a task. Interestingly, the reading of each group of HEMs is performed once, so authors do not assess the impact of variability in HEM readings due to hardware and software related variability. We, instead, focus in much more predictable environments, as needed for CRTES and consider the variability of HEM readings. HEM sampling or multiplexing consists in time-sharing the PMCs over a set of HEMs: at each interval boundary, whose duration is a configuration parameter, the PMCs are reprogrammed to read a different set ofHEMs. HEM sampling is, for instance, adopted by the Linux kernel’s perf event subsystem. The potential inaccuracies introduced by the interpolation made by sampling techniques have been studied elsewhere [103, 180]. Other works rearrange samples derived from HEMs from long executions to improve performance analysis [147]. At hardware level, other works [171] focus on specific HEMs (e.g. retired instructions, branches, 19 Chapter 2. Background loads/stores) and develop low-level hardware hypotheses on the reasons behind some of these HEMs suffering from various forms of under and over count. Authors do several recommendations on the hardware design to reduce the observed variability. Our goal, instead, is managing at software level HEM variability and limitations to read HEMs on existing boards (e.g. NXP T2080). Incomplete sets of data have been considered with Matrix Completion (MC) methods [83, 136]. There are fundamental differences between HEM merging and MC. 1. MC requires input data be a random matrix, where values in all rows and columns belong to one measure with its own distribution for each value, which does not hold for the problem at hand (e.g. the outcome of a HEM measure is a column, thus with its own distribution). 2. MC aims at producing synthetic data to complete missing data, thus bringing risks due to inferring the distribution of real data to produce new data, which may match some characteristics of real input data but miss others. In fact, since input data (HEM values) do not match requirement (1), and the fraction of missing values is large, MC populates the data array with values whose mean and standard deviation differs drastically from those for real (observed) data. Overall, MC does not fit the needs of the problem at hand by construction since its prerequisites are not met. Several solutions have been proposed to master multicore interference in industrial contexts. Some approaches attempt to reduce or fully remove interference across tasks running concurrently in multicoresbysegregatingaccessestodifferenthardwareblocks. Theseapproacheshavebeendevised for on-chip and off-chip memories, including banks of shared caches, as well as banks and ranks of DDR memories [94, 107, 110, 128, 155, 182]. Other approaches perform segregation over time, rather than over space, by splitting execution of tasks into memory and computation phases, thus letting schedulers guarantee that memory phases from different tasks do not occur simultaneously [25, 51, 129]. However, segregation over time is not always doable due to the characteristics of the hardware or the application itself (e.g. application semantics cannot be changed due to overwhelming V&V costs). Therefore, if time segregation is not feasible, even if space segregation is used, multicore interference can still occur in hardware resources not visible at software level, such as interconnects, as well as buffers and queues internal to shared caches for instance [161]. In those cases, solutions are needed to master interference in multicores and account for it during timing analysis. 2.5 Timing Validation in Automotive Systems The need for deriving time budgets for safety-related software components in automotive emanates from safety requirements, as described in ISO 26262 automotive functional safety standard [89]. Safety requirements determine the maximum response time affordable for a given functionality and the fault tolerant time interval, which stands for the time allowed since a fault occurs until an action is taken to recover or to bring the system to a safe state. Both determine the end-to-end time for a given functionality to be performed (including sensing and actuation time). Once discounted the time devoted to interfacing with physical components, the remaining time is the time budget allocated for the software component to complete its execution. Automotive systems are designed and verified following appropriate practices to achieve safety compliance. For instance, WCET estimation tools may be used, together with appropriate response timeanalysismethods for taskscheduling. Nevertheless, a validationstepis neededbeforedeploying the system to detect system faults due to, for instance, the violation of some assumptions caused during integration. Whether an application adheres to its timing constraints during the validation process is normally assessed in hardware-in-the-loop testing environments, where the complete application can be run. Appropriate test equipment is used to collect information on the execution of the application under analysis in general, and each of its tasks in particular, thus with a much higher observability than in the system during operation. However, due to the high coupling between the different tasks, the operating system, and the input/output interface of the application, individual tasks cannot be run in isolation or at externally-controlled time instants. Hence, although start and end times can be obtained for the different jobs of each task (observability), it is not possible to 20 2.5. Timing Validation in Automotive Systems enforce specific release times at will (controllability). Overall, tasks overlap with each other in an intricate manner depending on release times dictated by input data and timers, as well as by their duration. Execution time measurements, therefore, reflect arbitrary-looking overlaps across tasks, representative for the input data used during tests. Whether other overlaps are possible and whether those can lead to higher execution times is, in general, unknown and hard to control in practice. Hence, other than considering an engineering margin on top of the high watermark execution time to account for scenarios not triggered by the tests used, end users lack tools to use those execution time measurements in a more informed manner for validation purposes. While engineering margins set based on experience have worked for hardware platforms with low execution time variability, the increasing use of Graphics Processing Units (GPUs) and other high-performance hardware to match the performance needs of Autonomous Driving frameworks brings much higher performance variability due to contention in the use of shared resources (e.g. shared caches and main memory bandwidth). Hence, end users lack means to quantify whether unobserved scenarios could lead to timing violations. One of the focus of this thesis is to aim the timing validation methodologies by providing with an analysis tool for the contention derived from overlapping tasks. Continuing the trend of resorting to measurement-based analysis our contention modelling builds on a state-of-the-art autonomous driving framework. 2.5.1 State-of-the-art on MBPTA with Dependency The impact of dependencies in WCET estimation has been mostly considered for the particular case where they exist across jobs of a given task. In particular, probabilistic WCET estimation was originally formulated on top of the assumption of i.i.d. execution time observations for a task [59]. Therefore, solutions for pWCET estimation have often built on top of Extreme Value Theory for i.i.d. processes [65]. However, some researchers noted that execution time dependencies may exist across jobs of a given task, thus jeopardizing i.i.d. properties, so methods considering those dependencies would be more convenient [40, 100]. Amongst those works, Bernat et al. [22] pioneered in the area of execution time dependencies and analyzed them at the scope of program components. Santinelli et al. [142] considered a particular type of dependencies across jobs – stationary processes – andconcludedthattheycouldleadtoWCETunderestimationifnotaccountedforcarefully. Similarly, Melani et al. [114] showed that appropriate statistical tests (mostly correlation and independence tests) can be used to account for those dependencies satisfactorily. Lima and Bate [104] proposed a solution to mitigate the impact of dependencies across jobs to facilitate WCET estimation. Finally, Abella et al. [2] have recently shown that the source of those dependencies across jobs imposes different constraints on task scheduling. 21 Chapter 2. Background 22 Chapter 3 Experimental Methodology The workflow used in this thesis is carefully designed in order to ensure the results obtained can be trusted; from the generation of the data until the end goal. The diagram in Figure 3.1 shows the general workflow of all experimental methodologies followed in this thesis. In Section 3.1 we present the Data Generation process which in this thesis has been performed in two ways. •The first option is to design a set of programs that stresses the hardware (board) or simulator used to obtain execution times or HEM data for the analysis. •The second option, described in Section 3.2, is to use analytical parametric distributions as benchmarks for pWCET methodologies’ validation. Gathering a set of parametric distributions with similar tail profile than the most extreme conditions of CRTES at operation is useful to validate pWCET estimation methodologies. Once we have chosen the benchmarks for data generation, we perform the Data Gathering using sampling techniques. In Section 3.3 the Data Validation step is detailed. After sampling, the samples coming from real hardware or simulator should pass a validation process with an i.i.d. test to ensure that the methodologies can be trustingly applied. In contrast, the data sampled from parametric distributions is obtained using standard statistical libraries from R [133], in which the i.i.d. condition is given by construction. Finally, Section 3.4 shows the Modelling part of the experimentation. Once the data is ready and validated, we apply the models and methods we propose in this thesis and also state-of-the-art methodologies like EVT. The application of EVT models will be further explained in this chapter. After the confidence intervals have been computed, we obtain our End Result, be it a pWCET or other results from our methodologies. Even if the methodologies we investigated in this thesis are different in nature, the core of the experimental methodology remains the same to ensure trustworthy results and methods. Benchmark Programs Benchmark Distributions Hardware Sample i.i.d. Test EVT Model or other Methodologies Confidence Intervals pWCET Estimation or other Models Data Generation Data ValidationData Gathering Modelling End result Simulator Figure 3.1: Diagram of the experimental workflow with the tools used in this thesis. 3.1 Benchmarks and Platforms In this section we describe the benchmarks used in this thesis to generate the data that will be used to validate and test our methodologies. We distinguish between a real-case study on real hardware, 23 Chapter 3. Experimental Methodology Table 3.1: Workloads on the T2080 for validation purposes. Core0 Core1 Core2 Core3 W1 IMUL_UL2 FADD_MEM FMUL_MEM LMUL_MEM W2 LADD_DL1 LMUL_UL2 FADD_UL2 FMUL_UL2 W3 LMUL_MEM LMUL_UL2 FADD_UL2 DMUL_UL2 W4 LADD_UL2 FADD_MEM FMUL_MEM LADD_MEM W5 FADD_MEM FMUL_MEM LADD_MEM LMUL_MEM W6 FADD_UL2 FMUL_UL2 LADD_UL2 LDIV_UL2 W7 FADD_UL1 FMUL_UL1 LADD_UL1 LMUL_MEM W8 FMUL_UL1 FADD_MEM LMUL_L1 LMUL_MEM W9 FMUL_MEM FADD_DL1 FADD_MEM LMUL_DL1 W10 LADD_MEM LADD_DL1 LADD_UL2 LADD_MEM W11 FADD_DL1 FADD_UL2 FADD_MEM FMUL_UL1 W12 FADD_UL2 LADD_UL2 LDIV_DL1 LMUL_UL2 W13 LADD_UL2 FADD_MEM FMUL_MEM LADD_MEM W14 LADD_UL2 LMUL_UL2 FADD_MEM FMUL_MEM W15 FADD_MEM FMUL_MEM LADD_DL1 LDIV_DL1 W16 LADD_DL1 LDIV_DL1 FADD_MEM FMUL_MEM benchmarks that are run on real hardware boards to collect HEM data, a representative autonomous driving software framework, and a set of reference analytical distributions which are representative of tail profiles in CRTES. 3.1.1 Railway Case-Study on a LEON3+ Platform The European Train Control System (ETCS) is a safety-critical application (SIL 4) responsible of signaling and control in the European Rail Traffic Management System (ERTMS) framework. ETCS protects the train motion by constantly monitoring traveled distance and speed, and is programmed to activate the emergency brake system whenever unauthorized speed values are detected. The ETCS subsystem comprises three main tasks that are executed sequentially to provide the required safety function: the Odometry module, estimating a set of parameters based on the inputs collected from the train environment (e.g., estimated position); the Service module, managing the Service braking system; and the Emergency module, actually controlling the Emergency braking system. This function has strict real-time requirements and needs to be certified at the highest integrity level in the corresponding railway safety standards. While all three tasks do exhibit strict real-time requirements, we focus our evaluation on the Emergency module, the core of the ETCS CRTES module. The ETCS validation suite, which is made available with the application, includes 10 different input vectors (TEST0 to TEST9) corresponding to the operating conditions for functional and timing validation regarded as relevant by the application owner. These benchmarks are used for the techniques proposed in Chapters 4 and 5. Platform. The platform for this experiment is an Field-Programmable Gate Array (FPGA) implementation of a LEON3+ architecture comprising first level instruction and data caches, and a unified L2 cache, where the sources of execution time variation have been conveniently controlled by hardware means to guarantee the representativeness of the measurements collected at analysis with respect to the timing behavior during operation [85, 159]. In particular, this processor builds upon the concepts of time upper-bounding and time randomization to enforce representativeness. 3.1.2 Microbenchmarks on an NXP T2080 Platform In general, programs can have built-in sources of non-determinism (e.g. time- or input-dependent values). Also, they may easily be subject to variability due to minimal variations in the Operating System [171]. In order to reduce these sources of variability, we construct specific test-cases, which 24 3.1. Benchmarks and Platforms e6500 e6500 e6500 e6500 d$ i$d$ i$ d$ i$d$ i$ d$ i$d$ i$ d$ i$d$ i$ L2 Cache Bank Bank Bank Bank CoreNet DDR Memory DMA CPC Other IO · · · Figure 3.2: Block diagram of the T2080. also aim at triggering a wide set of HEMs. To that end, we have created different benchmarks comprising different core and cache (memory) patterns. At core level, we create 3 benchmarks using intensively the integer pipeline, using integer (I) and long (L) operands, and the floating point (F) pipeline. We use short latency addition (ADD) operations and long-latency multiplication (MUL) operations. At the cache hierarchy level, benchmarks operate on a vector whose size we vary so most load/store operations hit in the data cache (DL1), the L2 (UL2) or memory (MEM). From these 18 benchmarks (𝐼, 𝐿, 𝐹) · (ADD,MUL) · (DL1,UL2,MEM), we have generated 16 workloads, as shown in Table 3.1. These benchmarks are used for the techniques presented in Chapters 6 and 7. In both chapters, we target a NXP T2080 Reference Design Board [68, 69] increasingly considered in the avionics domain, with Airbus already achieving multi-core certification on this board [113], and with Rockwell Collins pursuing such certification [134]. The T2080 equips 4 e6500 cores (see Figure 3.2), each comprising private instruction and data cache as well as a private MMU. A second level cache is shared between all the cores. A “CoreNet” coherence fabric provides access to the memory controller as well as other peripherals present in the board. Some features are deactivated in our setup for predictability reasons, such as SMT (Hyperthreading in Intel terminology) and the CoreNet Platform Cache. We have run our tests in a bare-metal setup, using the software development kit provided by the boardmanufacturer (NXP)to configure the platform and loadimagestoit through a JTAG debugging interface. In the bare-metal setup, we access PMCs directly without the use of a specific library, e.g. PAPI, to minimize the impact of readings. In each experiment, we run one benchmark per core. The task in core0 is the reference task on which we perform the analysis, for the tasks in the other cores would be performed analogously. In each run of every experiment, we collect measurements when the task in core0 finishes its execution. We consider single-path benchmarks to isolate platform-level variability, so that in all runs the number of instructions executed (INSTRUCTIONS_COMPLETED) in core0 is exactly the same. Across any two runs of an experiment, we reset the state of caches, TLBs, and Branch Target Buffer. To that end, we execute a micro-benchmark that generates a massive number of misses in all those stateful blocks. While ISA-specific solutions exist that allow obtaining the same effect with specific instructions, we considered the micro-benchmark solution to be more platform agnostic. 3.1.3 Apollo Autonomous Driving Framework on an NVIDIA Jetson Platform In Chapter 8 we present a methodology to assess contention in CRTES. There, we target Apollo [11] AD framework, as the largest existing AD project with more than 120 partners including top-tier AI and tech companies, and car manufacturers. 25 Chapter 3. Experimental Methodology Light Heavy Exponential Figure 3.4: Example of a CV plot for a sample from Gaussian distribution. 3.4.2 Coefficient of Variation The Coefficient of Variation (CV) is a powerful test to distinguish between polynomial (light and heavy) and exponential tails [90]. It uses the residual coefficient of variation. The theoretical residual CV from the excess random variable 𝑋𝑢from a threshold 𝑢is given by CV(𝑢) ≡ CV(𝑋𝑢)=p𝑉(𝑢)/𝑀(𝑢),(3.10) where 𝑀(𝑢)=𝐸(𝑋𝑢)is the residual mean and 𝑉(𝑢)=Var(𝑋𝑢)is the residual variance. For a given sample we can compute the empirical residual CV with the standard deviation and mean of the sample. It was shown in [132] and [15] that the residual CV of a random variable, for a sufficiently large threshold, is almost constant and tends to the residual CV of a certain GPD. In the case of the GPD, the residual CV is constant and independent of the scale parameter, 𝐶𝜉=p1/(1−2𝜉).(3.11) Therefore estimating the CV one can estimate 𝜉: 𝜉=1 2 1−1 𝐶2 𝜉!.(3.12) This serves as a distinguisher for the kind of tail one is working with. Consider the case when 𝐶𝜉=1, then 𝜉=0, therefore our sample, 𝑋, is in the domain of attraction of a Gumbel distribution (exponential tails). Similarly, if 𝐶𝜉>1we are dealing with heavy tails, and if 𝐶𝜉<1we are dealing with light tails. In the same way as with the Hill plot, we can draw a plot that help us distinguish which tail we are dealing with. In Figure 3.4 we show an example of the CV plot (with the residual CV being the black line) for a sample of size 𝑛=105coming from a Gaussian distribution. First of all, we only use those order statistics that are closest to the tail. In Figure 3.4 only the last thousand order statistics are shown, which are the ones closest to the tail in this case. There are three zones to remark here. The red line represent the statistical test developed in [90] to reject the hypothesis of exponential tail. In the green zone in between the red lines, there is the region where we cannot 32 3.4. Estimating the Extreme Value Index reject the hypothesis of exponential tails. Our example data shows that on the order statistic between 99200 and 99250 the exponential tail of our distribution begins. The red regions are the ones where we can reject the exponential hypothesis because we are either dealing with heavy tails (upper zone) or light tails (bottom zone). This technique has been used throughout the thesis as a first assessment on the kind of tails we are working with. The tools to implement this concept are in the CRAN package "ercv" [53]. 33 Chapter 3. Experimental Methodology 34 Part II Sky-High Quantile Estimation for CRTES 35 Chapter 4 Sky-high Quantile Estimation with Weibull Tails 4.1 Introduction The increasing automation of CRTES, such as those in cars and planes, leads to more complex and performance-demanding on-board software and the subsequent adoption of multicores and accelerators. This causes software’s execution time dispersion to increase due to variable-latency resources such as caches, networks on chip, advanced memory controllers and the like. CRTES undergo strict V&V processes against their corresponding functional safety standards (i.e. ISO26262 in automotive [89] and DO178B/C[139] in avionics). Timing V&V processes require collecting evidence supporting that software will execute correctly and timely. In particular, those processes provide evidence that the risk of violating deadlines for critical real-time software is residual, since functional safety standards acknowledge that risk cannot be completely avoided. Thus, they impose thorough V&Vprocessesthatallow deemingsuchriskas“sufficientlylow”sothatitcanbeneglected. In this line, pWCET estimates allow to quantify such residual risk. EVT [38], appropriate for risk analysis, is used to model the right/upper tail of execution time distributions. In the context of pWCET estimation, exponential tails delivered by EVT, see Figure 1.2, have been shown by argument [7] and empirically [156] to provide a reliable tail model for pWCET estimation. However, tightness of those exponential tails is limited. Tails lighter than exponential ones (so with 𝜉 < 0) can deliver tighter bounds, as discussed in [3] and later in Section 4.2. Yet, in the context of EVT, either GEV or GPD, distributions with 𝜉 < 0have a compact support, i.e. they have an absolute maximum value that cannot be exceeded. Hence, light tails in the case of EVT have an intrinsic risk of delivering optimistic tail distributions. This has some key implications in the fitting process, since a sufficiently large sample is needed to guarantee that light tail fitting is reliable for arbitrarily low exceedance probabilities. The target of this chapter is overcoming this limitation of the data and delivering a practical solution to obtain pWCET estimates tighter than those of exponential tails while preserving reliability. We do so by complementing EVT with survivability analysis as the theoretical ground for our hypothesis. Analysis. We formally show that tail modelling in the context of survivability analysis [43] targets analogous questions to those of risk analysis. Then, we show that log-concave1distributions [112], inherited from survivability analysis, deliver the tightest distribution models, but existing fitting methods fail to model tails [14, 58], needed for pWCET estimation. Proposal. We propose the use of a subset of log-concave distributions: Weibull tail distributions with increasing hazard rate, i.e. with shape 𝛽≥1, neither Reverse Weibull as in EVT nor full Weibull as in survivability analysis. Our approach provides analogous accuracy to that of logconcave distributions, without limitations to extend them to arbitrarily low exceedance probabilities, as needed for pWCET estimation. 1A function 𝑓is logarithmically concave log-concave for short, if the function log(𝑓(𝑥)) is concave. 37 Chapter 4. Sky-high Quantile Estimation with Weibull Tails Assessment. We compare our method and EVT alternatives with a bootstrap analysis on large data samples(107measurements) froma railwaycase studyas groundtruth. Our results provideevidence on the reliability of our approach and significant pWCET reductions with respect to exponential tails obtained with EVT, thus allowing to trustingly increase system utilization noticeably. 4.2 EVT limits Threshold Figure 4.1: EVT estimation on the mixture distribution (three Gaussian with parameters 𝜇={5,50,100},𝜎={10,10,10}and weights 𝑤={0.60,0.39,0.01}, respectively.) CCDF Figure 4.2: CCDF for GPD (𝜉 < 0) and exponential tails from EVT for TEST8. Both models are fitted with a sample of 𝑛=1000 observations out of all 𝑛=107observations made. CRTES programs have a finite WCET as discussed in Section 2.1. A distribution with a finite maximum value is said to have compact support. In that situation the most appropriate tail profile to model the input data with will be a light one, i.e. with extreme value index 𝜉 < 0. In [44] it has been argued that heavy-tails, i.e. 𝜉 > 0, are too pessimistic to model distributions with compact support. As an upper limit for the model, exponential tails, 𝜉=0can be proposed. As an example, we show the results of this reasoning. In Figure 4.2 we show the fit of an execution time profile from a program running in CRTES. The execution time profile contains 𝑛=107samples. Although for modelling purposes, we tried to fit a GPD with an exponential tail and a light tail using only a sample of 𝑛=103realizations. In the figure we can appreciate how the exponential model in this case is 38 4.3. On the Use of Light Tails and Risk Analysis for WCET Estimation overly pessimistic for probabilities as low as 𝑝=10−6. On the other hand, the GPD model with light tails is closer to the real distribution, however the model underestimates, thus being unsafe in the context of MBPTA. The reason the GPD with light tails underestimates, is because its compact support nature. When trying to fit the small sample of 𝑛=103with a light tail, the maximum value on that sample will serve as the reference for the GPD. Therefore, if the maximum value for the sample is smaller than the maximum value of the whole sample, the GPD will underestimate. Nonetheless, the GPD with light tails is a good model in a vacuum for these kinds of distributions. However, in modelling the context always matters. Therefore, while in general the GPD with light tails is a good model, in the context of MBPTA it is a risky model. We have just discussed how the exponential tail can be used as an upper-bound model for execution time profiles of CRTES. Recall that we discussed in Section 2.2.3 how the second fundamental theorem of EVT, Theorem 3, was fulfilled when the threshold 𝑢→ ∞ and that it could bring uncertainty to the estimation. Figure 4.1 shows a mixture of three Gaussians with 𝜇={5,50,100}, 𝜎={10,10,10}and weights 𝑤={0.6,0.39,0.01}. Each mode around Time {5,50,100}represents a Gaussian in the mixture. For each different tail threshold 𝑢such that 𝑢>60, we fit an exponential tail, given that Gaussian distributions have exponential tails. In Figure 4.1 we see three scenarios that we exemplify with approximate ranges of 𝑢. For values of 𝑢around [60,80](purple lines), we see how the exceedance probability of the mixture is underestimated in the Time range [80,120] (probabilities 10−3−10−5). For 𝑢in the range [120,140](blue-green lines), GPD over estimates the reference distribution. For values of 𝑢above 140 (green lines), the estimate becomes increasingly tight. Depending on the threshold selected, even if the exponential is the appropriate model to fit the tail of a mixture of Gaussians, it can lead to optimistic pWCET. In that regard, one of the goals of this thesis is to overcome the limits of EVT when working in the context of MBPTA. More specifically, in this work our goal is to devise a more flexible model for the tails than the GPD with light tails and the exponential by making use of Risk Analysis theory. 4.3 On the Use of Light Tails and Risk Analysis for WCET Estimation 4.3.1 Risk and Survivability Analysis EVT is used in risk analysis to predict extreme (rare) events with the objective of proving that risk is below specific thresholds (e.g. financial risk). Survivability analysis, while it also focuses on predicting extreme events, it has opposite goals: proving that survivability is above specific thresholds (e.g. human life duration). Hence, both analyses target the modelling of extreme events, but with different objectives. The type of distributions used to model those extreme events (i.e. tail distributions) differs across both analyses. •EVT (either GPD or GEV) is used in the context of risk analysis, and it has been used so far in pWCET estimation building on the idea that exceeding a specific execution time bound is a risk. •For survivability analysis tail distributions can be split into two types: Decreasing Hazard Rate (DHR) and Increasing Hazard Rate (IHR) distributions. The boundary between those two categories corresponds to exponential distributions, which can be regarded as part of both. In the context of pWCET estimation, we focus on IHR distributions, since they include those distributions that, as values get higher, the probability of realization increases. In our context this means that, as the program runs, there is an increasing probability of finishing the execution, which is the case of real-time programs that need to have a finite execution time to meet their deadline. Formally stated, a random variable 𝑋is IHR if the hazard rate function is increasing, where the hazard rate function is defined as: ℎ(𝑥)=𝑓(𝑥) 1−𝐹(𝑥), 𝑥 ∈support(𝑋),(4.1) where 𝑓and 𝐹stand for the Probability Density Function (PDF) and the CDF of 𝑋, respectively. In Equation 4.1, support(𝑋)is a function representing the subset of the domain in which the random 39 Chapter 4. Sky-high Quantile Estimation with Weibull Tails variable (𝑋) probability is defined (i.e. it is not zero). In fact, the hazard rate function which assesses the IHR property is equivalent to the convexity 𝐻function, where 𝐻(𝑥)=−log(1−𝐹(𝑥)), called cumulative hazard rate. 4.3.2 IHR Distributions in Survivability Analysis Log-concave distributions, as they can have an arbitrarily large number of parameters, are one of the best tools to fit data. On the other hand, they are generally defined only for the range of data observed. Hence, they are unable to model the distribution beyond that range, which is not useful for pWCET estimation. Some methods smooth the distribution by means of convolutions with other laws (e.g. Gaussian) to better fit the mode of the distribution [14, 58, 112]. While such an approach delivers a distribution that spans beyond the range of the data observed, it is tuned to model central behavior (not the tails) and so inherits the original problem we aim to tackle: an appropriate law needs to be identified for tail modelling. Weibull distributions, among others, have often been used to model IHR distributions. However, they are unable to fit all tail distributions, so they are used only in specific contexts where the problem at hand matches the shape of those distributions [112, 174]. Therefore, while survivability analysis opens new opportunities to model pWCET, this is yet unexplored and existing distributions, in general, may not fit the needs of pWCET estimation. In this chapter we address this challenge by contextualizing the needs of pWCET estimation and defining distribution families able to model pWCET distributions reliably, tightly and without incurring on the limitations imposed by logconcave distributions. 4.4 Equivalence Between IHR and Non-heavy Tails In this section we introduce the connection between risk and survivability analysis in the context of pWCET estimation, which will lay out the ground for our hypothesis in further sections. pWCET estimates should be obtained theoretically under the assumption of the law of extreme events characterized by quicker decay in the tail than an exponential law, or equal, in the limit case. Exponential tails have been regarded as the appropriate (limit) model in practice [3, 149]. An exponential decay is a memoryless process where the probability of the process to complete is constant regardless of how long the process has been progressing. In the context of pWCET estimation, this corresponds to having a constant probability for the program to finish its execution despite the time elapsed since the program started running. Instead, the theoretical solution for pWCET estimation indicates that, the longer the program has been running, the higher the probability of finishing. That is, if 𝑋 corresponds to execution time as a non-negative random variable, this assumption can be formally stated as follows if 𝑠<𝑡and 𝑥>0: 𝑃(𝑋>𝑡+𝑥|𝑋>𝑡) ≤ 𝑃(𝑋>𝑠+𝑥|𝑋>𝑠),(4.2) In the context of extreme events, where 𝑢is the threshold upon which values belong to the tail of the distribution, 𝑠, 𝑡 have to be large enough so that 𝑠, 𝑡 >𝑢, for some 𝑢>02. In the general case of EVT (e.g. GPD), the tail decay can be described as follows for heavy, exponential and light tails respectively: if 𝜉 > 0,then 𝑃(𝑋>𝑡+𝑥|𝑋>𝑡) ≥ 𝑃(𝑋>𝑥)for 𝑥>𝑡. if 𝜉=0,then 𝑃(𝑋>𝑡+𝑥|𝑋>𝑡)=𝑃(𝑋>𝑥)for 𝑥>𝑡. if 𝜉 < 0,then 𝑃(𝑋>𝑡+𝑥|𝑋>𝑡) ≤ 𝑃(𝑋>𝑥)for 𝑥>𝑡. Note that the assumption described by Equation 4.2 for extreme events implies 𝜉≤0. Hence, the assumption for pWCET estimation matches the formulation above for light and exponential tails from GPD, inherited from risk analysis since, in the context of pWCET estimation, the longer the 2Note that the threshold 𝑢is not larger than the theoretical threshold corresponding to the asymptotic behavior of the tail from the second fundamental theorem in EVT [15, 132]. 40 4.5. Weibull Tails (TailW) for pWCET Estimation program has been running, the higher the probability of finishing. This matches the known concept in reliability modelling of IHR, which is therefore appropriate for pWCET estimation. Building on Equation 4.2, and given a fixed 𝑥, we have the following equivalent assumption: 𝑃(𝑋<𝑡+𝑥|𝑋>𝑡) ≥ 𝑃(𝑋<𝑠+𝑥|𝑋>𝑠),if 𝑠<𝑡. Given that 𝑠<𝑡, then 𝑃(𝑋>𝑡)<𝑃(𝑋>𝑠), and hence, we can elaborate the equation above as follows: 𝑃(𝑋<𝑡+𝑥)−𝑃(𝑋<𝑡) 𝑃(𝑋>𝑡)≥𝑃(𝑋<𝑠+𝑥)−𝑃(𝑋<𝑠) 𝑃(𝑋>𝑠). If we use the CDF expressions instead, we have the excess distribution: 𝐹(𝑡+𝑥)−𝐹(𝑡) 1−𝐹(𝑡)≥𝐹(𝑠+𝑥)−𝐹(𝑠) 1−𝐹(𝑠),if 𝑠<𝑡. (4.3) If 𝑥tends to 0, then the cumulative probability ranges, between 𝑡and 𝑡+𝑥and between 𝑠and 𝑠+𝑥, reduce to the particular probabilities at 𝑡and 𝑠respectively: ℎ(𝑠)=𝑓(𝑠) 1−𝐹(𝑠)≤𝑓(𝑡) 1−𝐹(𝑡)=ℎ(𝑡),if 𝑠<𝑡. (4.4) As shown, Equation 4.4 – which we derive from Equation 4.2 – builds upon the hazard rate function shown in Equation 4.1. Hence, it models IHR distributions for survivability analysis, analogously to the GPD formulation for risk analysis. Note that in Equation 4.4 the equality case corresponds to a constant decay rate, hence a constant hazard rate function. Log Concavity: In order to use IHR distributions for pWCET estimation, we build upon the following theorem proven in [43] and [84]: Theorem 5. Given a non-negative random variable 𝑋, with 𝑓and 𝐹the PDF and CDF, respectively (where 𝐻(𝑥)=−log(1−𝐹(𝑥)), 𝑥 ∈support(𝑋)), log(𝑓)concave ⇒𝑋IHR ⇔𝐻convex.(4.5) Note that 𝑋is IHR in the tail, i.e. (𝑋|𝑋>𝑢)is IHR for some threshold 𝑢>0, if and only if Equation 4.2 holds for all 𝑠, 𝑡 >𝑢and, therefore, 𝑋is log-concave. Thus, by using log-concave distributions, IHR holds by construction. A non-negative function is log-concave if its domain is a convex set, and if it satisfies the inequality 𝑓(𝜃𝑥+(1−𝜃)𝑦)≥𝑓(𝑥)𝜃𝑓(𝑦)1−𝜃, for all 𝑥, 𝑦 in the domain of fand 0< 𝜃 < 1. In order to test IHR one could make use of the log-concavity of the probability density function, which would give a convex Hfunction and hence, IHR. Given an appropriate threshold 𝑢so that (𝑋|𝑋>𝑢)is IHR (and log-concave), we can fit a logconcave densityfunction tothe tailby usingthe maximumlikelihoodapproach, as detailedin [14, 58]. Regardless of whether we fit the best log-concave distribution or a distribution function family preserving log-concavity but with much fewer parameters, as we do in this work, the exceedance threshold (𝑢above) must be estimated to use the appropriate set of tail values from the sample for fitting. In particular, we build upon the work by Hazelton [84] that provides a procedure for testing whether we can reject the hypothesis of log-concavity for a given threshold 𝑢. 4.5 Weibull Tails (TailW) for pWCET Estimation In Section 4.2 we have concluded that light tails with compact support are likely optimistic (thus unreliable) for pWCET estimation. On the other hand, exponential tails are the limit distribution for appropriate pWCET distribution models, hence being reliable but likely pessimistic. In order to 41 Chapter 4. Sky-high Quantile Estimation with Weibull Tails 4.8 Summary The need for performing a proper analysis is crucial for the V&V of CRTES. While EVT has been proposed and applied successfully to derive WCET estimates, its results can be improved upon. The first alternative we proposed in this thesis to compute pWCET is the Tail of a Weibull (tailW), which combines characteristics from EVT and Risk Analysis. As shown in the results, tailW is less pessimistic than the exponential tail while not being optimistic like GPD with light tails model, thus producing tight pWCET models. 48 Chapter 5 Sky-high Quantile Estimation with Markov’s Inequality 5.1 Introduction DerivingWCET estimatesforsoftwareprogramswith probabilisticmeans (a.k.a. pWCET estimation) has received significant attention during the last years as a way to deal with the increased complexity of the processors used in real-time systems. Many works build on EVT that is fed with a sample of the collected data (execution times). In its application, EVT carries two sources of uncertainty [26]: model uncertainty that is intrinsic to the EVT model and relates to determining the subset of the sample that belongs to the (upper) tail, and hence, is actually used by EVT for prediction; and statistical uncertainty that is induced by the sampling process and hence is inherent to all samplebased methods. •Statistical uncertainty encompasses as the first aspect the testing conditions under which the experiments are performed in reference to those that can arise during system operation. Testing conditions that are representative or worse than operation conditions are the basis to attain representativeness of the sample data (execution time) [4, 12, 78, 118] so that the pWCET estimate holds during system operation. A second aspect of statistical uncertainty relates to the natural uncertainty of a sampling process that, in general, reduces as the sample size increases, andthatishandledwithconfidenceintervals. Samplinguncertaintyimpactssummarystatistics (e.g. mean) and tail fitting methods, whose goodness – either of their hypotheses or outcome – is assessed with specific methods [12, 137]. •Model uncertainty, instead, relates to uncertainties intrinsic to the mathematical model used for tail prediction. In the case of EVT, model uncertainty relates to determining the threshold from which the upper tail starts. This threshold plays a key role on the trustworthiness (safeness) of EVT results since only samples above it (i.e. the maxima data set) are fed into EVT for pWCET estimation. There is not an exact mathematical method to derive this threshold. Instead, current methods estimate the tail of a distribution [30] based on plot inference [50, 52, 86] and regression analysis [29]. As discussed in Section 2.2.3, the second fundamental theorem of EVT is satisfied when, for the excess distribution function 𝐹𝑢in Equation 2.12, the threshold 𝑢→ ∞. It is crucial to make a good estimation of this threshold because it will affect the pWCET as seen in Figure 4.1. Given that the computation of the threshold is not exact and there does not exist a closed form to compute it, one must rely on estimations to assess the tail of a distribution. This phenomenon is an example of model uncertainty. The model itself, in this case EVT, is the one bringing uncertainty in the calculations before dealing with the sampling process. To tackle this issue, we introduce Markov’s Inequality, a probabilistic tool which produces provably safe upper-bounds by construction, which fits the needs for pWCET estimation. In this chapter, we show that Markov’s Inequality is a provably safe probabilistic tool which upper- 49 Chapter 5. Sky-high Quantile Estimation with Markov’s Inequality (a) Gaussian distribution (b) Weibull distribution (c) Beta distribution (d) Gaussian Mixture Figure 5.1: Markov’s Inequality bound for the reference distributions. bounds by construction since it does not have any source of model uncertainty. The modification of Markov with the power-of-k function leads to tight upper-bounds even for very low probabilities. Finally, we show how to apply this probabilistic tool for any collection of samples coming from execution times of Real-Time systems. We carry out the development of our methodology using analytical reference distributions which are representative of the tail profiles from execution time distributions of Real-Time systems. This reference distributions are described in Section 3.2 of the Experimental chapter. Furthermore, we assess the validity of our methodology with a railway case-study data coming from the European Train Control System (ETCS), which consists on 10 different test from the Emergency Module. 5.2 Chebyshev and Markov Inequalities for pWCET Estimation This section analyzes the applicability of Markov’s inequality as an alternative model to EVT for the problem of trustworthy pWCET estimation and shows that it is not subject to any model uncertainty, thus resulting in provably safe bounds for the analyzed distribution. Recalling the discussion from Chapter 2 in Section 2.1; one upper-bounds the execution time profile of the task under analysis, when the resulting execution time at a given probability 𝑝for the pWCET ,etbound(𝑝), is higher (or equal) than the one of the task under analysis, etreal(𝑝). This can be expressed as tightness(p)=etbound(𝑝) etreal(𝑝). We introduced in Section 2.3.1 Chebyshev’s Inequality, and one of its realizations – Markov’s Inequality–, which are both probabilistic tools that satisfy MBPTA requirements for the pWCET problem because produce safe upper-bounds by construction. This chapter will be dedicated to study the viability and applicability of these inequalities to estimate sky-high quantiles. 5.2.1 Markov’s Inequality on Low Probabilities Besides trustworthiness, pWCET estimates are also required to be reasonably tight, especially for the range of relevant probabilities usually considered for pWCET estimation, e.g. [10−6,10−15]. In this line, our analysis shows that Markov’s inequality tends to be hardly useful for pWCET estimation. This is better illustrated in Figure 5.1 which shows for all considered distributions the probability bounds given by Markov’s inequality. We can observe that estimates are very pessimistic, orders of magnitude higher than the real probability. This includes the range of probabilities of interest for pWCET estimation. In fact, we seethat Markov’s inequalitynever goesbelow10−2forall distributions for the execution time value range plotted. Observation 2. Markov’s inequality in its original form is too pessimistic to be usable in practice for pWCET estimation. 50 5.3. Power-of-kfunctions for Markov’s Inequality (a) Gaussian distribution (b) Weibull distribution (c) Beta distribution (d) Gaussian Mixture Figure 5.2: MIK bounds for the reference distributions. 5.3 Power-of-kfunctions for Markov’s Inequality One of the main insights of this work is that the key reason for Markov’s inequality resulting in loose pWCET bounds lies on the fact that it builds on the identity function, 𝑓(𝑋)=𝑋, of a random variable. In this section, we show how a different function can lead to increased tightness on the produced pWCET estimates while preserving trustworthiness. In particular, we contend that the power function for any 𝑘∈R>0={𝑘∈R|𝑘>0}, i.e. 𝑓(𝑋)=𝑋𝑘, or power-of-𝑘function for short, can be safely used instead of the identity function to obtain tighter and trustworthy pWCET bounds. Definition 4. Let 𝑋be a discrete random variable and 𝑘a positive real value. The expected value of 𝑋𝑘, also defined as the 𝑘th (theoretical) moment, is: 𝐸(𝑋𝑘)=Õ 𝑥 𝑥𝑘𝑃(𝑋=𝑥).(5.1) Corollary 2 (Markov’s Inequality to the power-of-k).Let 𝑋>0and let the function 𝑓be the power-of-𝑘 function 𝑓(𝑋)=𝑋𝑘. Markov’s inequality to the power-of-𝑘yields: 𝑃(𝑋≥𝑏) ≤ 𝐸(𝑋𝑘) 𝑏𝑘.(5.2) Hence, the probability that 𝑋takes a value greater or equal to 𝑏is bounded by 𝐸(𝑋𝑘)/𝑏𝑘. This makes Markov’s inequality with 𝑓(𝑥)=𝑥𝑘(MIK for short) a safe pWCET estimate when 𝑋represents the execution time of a program. Proof. Theorem 4 holds true when Property 2.16 and Property 2.17 are fulfilled. The power-of-k function does not fulfill those properties in general. However, when the 𝑥domain is restricted to the positive real numbers R>0={𝑥∈R|𝑥>0}, which in fact includes the domain of execution time profiles, the power-of-kfunction does fulfill Properties 2.16 and 2.17 since 𝑥is positive, so 𝑥𝑘is also positive and an increasing function. □ Overall, for this application scenario (𝑥∈R>0), Equation 5.2 is an upper-bound when using the power-of-𝑘function onto the reference distribution for any value of 𝑘∈R>0. Hence, it can be leveraged for pWCET estimation. It is worth noting that other functions can exist that fulfill Properties 2.16 and 2.17. While exploring them is part of our future work, as shown in Section 5.3.1, MIK (i.e. 𝑓(𝑋)=𝑋𝑘) achieves very tight pWCET estimates which leaves small room for improvement. Observation 3. For every value of 𝑘∈R>0, MIK (i.e. Markov’s inequality with 𝑓(𝑋)=𝑋𝑘) is a safe pWCET estimate when 𝑋represents the execution time of a program. 51 Chapter 5. Sky-high Quantile Estimation with Markov’s Inequality (a) Gaussian distribution (b) Weibull distribution (c) Beta distribution (d) Gaussian Mixture Figure 5.3: Evolution of MIK bounds with the value of 𝑘. Note that there is no theoretical constraint on the maximum value of 𝑘, which can be any positive real number 𝑘. 5.3.1 Tightness of MIK for Increasing Values of k Once we have established the safe use of MIK for pWCET estimation, we illustrate the impact of varying 𝑘on tightness. Figure 5.2 shows for several values of 𝑘,{10,15,25,50}, and the reference distributions presented before, that MIK dramatically increases the tightness provided by Markov’s inequality (Figure 5.1), while remaining a trustworthy upper-bound for every value of 𝑘. We also observe that MIK tightness remains for high exceedance probabilities, which hence makes it a promising model to provide pWCET estimates. Observation 4. Markov’s inequality with 𝑓(𝑋)=𝑋𝑘heavily reduces the pessimism of Markov’s inequality. Intuitively, from Figure 5.2, higher values of 𝑘result in tighter estimates, i.e. minimizing the distance between the reference distribution and the upper-bound distribution. However, this is not always the case. For instance, if we take a closer look at the Beta distribution (Figure 5.2(c)), we see that at cut-off probabilities 10−3and 10−6the tightest MIK estimates are not obtained for the highest value of 𝑘evaluated (50). Observation 5. For a given threshold probability, higher values of 𝑘do not necessarily result in a tighter MIK bound. This is better illustrated with the examples in Figure 5.3 that shows quantitatively the evolution of the MIK bound obtained for varying values of 𝑘. In this experiment, for the value of the target distribution at each probability, we evaluate MIK for different values of 𝑘. As it can be seen, for every threshold probability and distribution the value of 𝑘resulting in the tightest estimation is different. For instance, for the Gaussian distribution (Figure 5.3(a)) and target probability 10−6,𝑘=71 produces the tightest estimate, while for 10−9and 10−12 the best 𝑘is 97 and 121, respectively. As a general trend, we see that the lower the target exceedance probability, the higher the value of the best 𝑘is. Yet, the highest value of 𝑘evaluated for each target probability does not produce the tightest bound. Overall, for each probability there exists a value of 𝑘producing the tightest upper-bound, with the optimal value of 𝑘depending on the actual reference distribution. Observation 6. Increasing the tightness of MIK for each probability is an optimization problem on 𝑘which only increases accuracy and does not affect trustworthiness. In order to address this optimization problem, we propose MEMIK (Minimum Envelope for MIK, i.e. Markov’s Inequality to the power-of-K). MEMIK combines the results of MIK bounds obtained for any value of 𝑘by keeping, for each point in an interval, the value of 𝑘producing the tightest estimate. This set of points form an envelope that is a provable trustworthy and tight tail bound by construction for any exceedance probability. Therefore, MEMIK improves the pessimistic upperbounds of Markov’s inequality (see Figure 5.1), with a much tighter envelope that is usable for 52 5.3. Power-of-kfunctions for Markov’s Inequality Algorithm 1 Compute an envelope using power-of-𝑘 1: function MEMIK(𝑡range,𝑡step,𝑝all,max𝑘(𝑝all),𝑘step) 2: for 𝑡∈𝑡range, 𝑡step do 3: for 𝑝∈𝑝all do 4: mikbest(𝑝) ← ∞ 5: for 𝑘∈ ([1,max𝑘(𝑝)], 𝑘step)do 6: vpred ←EVAL(𝑡, 𝑘, 𝑝) 7: if vpred <mikbest(𝑝)then 8: mikbest(𝑝)=vpred 9: 𝑘best(𝑝)=𝑘 10: end if 11: end for 12: envelope(𝑡, 𝑝) ←<mikbest(𝑝), 𝑘best(𝑝)> 13: end for 14: end for 15: return envelope 16: end function k (a) Gaussian distribution (b) Weibull distribution (c) Beta distribution (d) Gaussian Mixture Figure 5.4: MEMIK bound (envelope) on the reference distributions. pWCET estimation. Formally, the MEMIK bound is defined as follows: 𝑃(𝑋≥𝑏) ≤ min 𝑘 𝐸(𝑋𝑘) 𝑏𝑘for 𝑘>0.(5.3) MEMIK, see Algorithm 1, which uses point-wise power-of-𝑘Markov’s inequality values, performs a simple complete exploration of MIK values over a given time range 𝑡range and over a configurable range of 𝑘, determined by the maximum value max𝑘to be explored for each probability 𝑝in the set of probabilities of interest 𝑝all. The granularity of MEMIK exploration over 𝑡and 𝑘is determined by the 𝑡step and 𝑘step parameters respectively. For each probability 𝑝in the interval of interest 𝑝all (line 3) and 𝑘in the range determined by max𝑘(𝑝)(line 5), the algorithm estimates the value of the target distribution, EVAL(t, k, p) (line 6). To that end, we evaluate Equation 5.2 with 𝐸(𝑋𝑘), which corresponds to the theoretical 𝑘th moment of the target distribution from Equation 5.1, obtaining vpred that we compare to the best MIK value so far for all considered 𝑘(line 7). The minimum MIK value produced for a given 𝑡and across all 𝑘values (mikbest(𝑝)) is stored, together with the corresponding 𝑘best(𝑝), in the data structure envelope (line 12). Eventually, after iterating over the whole-time interval, the algorithm returns the envelope data structure (line 15) which holds the point-wise definition of the approximation envelope. Note that, if for any value 𝑡, the value of 𝑘best(𝑝)matches max𝑘(𝑝), then max𝑘(𝑝)can be increased to find tighter bounds. We applied MEMIK to our reference distributions, for which we can derive the theoretical moments. For this experiment, we varied the value of 𝑘up to 150 with 𝑘step =1. We obtained the envelopes depicted in Figure 5.4 which provides evidence that MEMIK produces tight and trustworthy estimates for all distributions, with an observed error of around 5%. 53 Chapter 5. Sky-high Quantile Estimation with Markov’s Inequality Figure 5.5: MEMIK with sample moments (𝑛=1000 and 𝑛sims =100) on the reference Gaussian distribution with (RESTK) and without (NO RESTK) restricting 𝑘. Also, MEMIK evaluated with theoretical moments. Maxk Figure 5.6: Minimum max𝑘values for the reference distributions used in this work. Overall, this section provides the key result that our proposal, Markov’s inequality to the power-of-k (MIK), unlike EVT, suffers no model uncertainty at the theoretical level and hence, provides correctby-construction pWCET estimates that are much tighter than those provided by the default Markov’s inequality. This leaves sampling uncertainty as the problem to address. 5.4 Handling Markov Sampling Uncertainty So far, we have been reasoning on examples for which we could compute the theoretical 𝑘th moments for each distribution. This was possible since the distributions considered were known and, hence, we could compute exactly the value of each moment (i.e. 𝐸(𝑋𝑘)for each value of 𝑘) using its analytical closed form. However, we need to consider the scenario in which only samples are available. Hence, as for any other sample-based method, we need to deal with the underlying sampling uncertainty. A commonality of sample-related methods like EVT [35], and something that we also assume, is that, input samples are independent and identically distributed [40] (i.i.d.) or at least exhibit extremal independence [142]. The i.i.d. property can be pursued with platform randomization [97] or data (time measurements) sample randomization [104]. For the case of the Markov’s inequality, this translates into deriving the sample moments, referred to as ˆ 𝐸(𝑋𝑘)(for each value of 𝑘). In particular, we need an estimator for high-order moments that can produce good estimates for any distribution. 5.4.1 Sample Moment Estimation The 𝑘th moment of a random variable 𝑋can be estimated as: ˆ 𝐸(𝑋𝑘)=1 𝑛 𝑛 Õ 𝑖=1 𝑋𝑘 𝑖.(5.4) In general, this estimator is the best one to deal with high-order sample moments [81], as it is asymptotically unbiased. Given that it asymptotically tends to a Gaussian distribution [16], the properties of the Central Limit Theorem’s apply to it. However, the estimator is asymptotically unbiased [64] only when using large amounts of data. For instance, for a sample of a Gaussian distribution with 𝜇=100,𝜎=10 and 𝑛=103, the difference between the 3rd exact moment, and the sample moment using Equation 5.4 is about 0.02%, it is between 1% and 3% for the 50th moment, and can be up to 160% for the 100th moment. 54 5.4. Handling Markov Sampling Uncertainty When the sample moment, i.e. the estimate of the 𝑘th moment, is higher than the theoretical moment, there is a risk ofunderestimatingtheupper tailof the distributionby assigningtoa certainprobability a smaller quantile than it has in reality. The approach we propose to limit that risk consists in setting a maximum value of 𝑘allowed for each different probability, which we refer to as max𝑘. In order to illustrate the effect of not controlling max𝑘, Figure 5.5 assesses the tightness of MEMIK over a particular set of nine probabilities (from 10−7to 10−15). The red triangles represent the theoretical bound of MEMIK using Equation 5.3 with 𝐸(𝑋𝑘)being the theoretical moments, while the orange dots (NO RESTK) represent the application over multiple simulations of Equation 5.3 using the sample moment from Equation 5.4. On both applications we set up a high value for max𝑘,150. We can see how the loss of consistency of the sample moment estimator on Equation 5.4 results in bad estimates. An intuitive way to control the gap between the 𝑘th sample and the theoretical moments is to vastly increase the sample size, which in our domain would require an unaffordable number of runs of the task under analysis. Alternatively, one can control the range of values of kexplored in Equation 5.4. By doing so, we trade some tightness for trustworthiness. That is, if we explore values of 𝑘until a low max𝑘limit, we can see in Figure 5.3 that the theoretical bound is not optimal in terms of accuracy. On the other hand, small values of max𝑘also limit the inaccuracy of the sample moment estimator. Note that it is not possible to identify a general optimal value of max𝑘for any kind of data under analysis. The appropriate max𝑘value changes across distribution types, across the same distribution type with different parameters, and even across probabilities for a given distribution. For this reason, we propose the restricted k method (RESTK) that builds on the information gathered directly from the samples to derive max𝑘so as to produce trustworthy and tight results. 5.4.2 Understanding the Behavior of max𝑘 We gain insight on the behavior of max𝑘along 3 axes. We analyze i) whether for a given distribution there exists a pattern for max𝑘that can provide tight and safe results using the sample moment estimator; ii) whether this pattern can be predicted using only the information from the sample; and iii) whether the pattern can be generalized for any distribution. We focus on the same example distributions used in previous sections. We fix the interval [1,150] as exploration range for 𝑘. In order to account for sample uncertainty, we perform 𝑛sims =103 Monte-Carlo simulations, each one considering a random sample of size 𝑛=103. We first compute Markov to the power-of-𝑘(MIK) using Equation 5.2 with the sample moment estimator in Equation 5.4, for all selected 𝑘. In each simulation, we increase values of 𝑘and find the first (smallest) value of 𝑘that produces underestimation. This is computed by comparing the estimation with the actual value of the distribution. That is, we take the value of 𝑘for which the estimated quantile is smaller than the theoretical quantile. Then, we set 𝑘−1as our max𝑘. For each Monte-Carlo simulation, we compute the max𝑘for all target probabilities. When all simulations are performed, we keep the smallest max𝑘for each probability. As a result, for each experiment of 𝑛sims Monte-Carlo simulations, we obtain one value of max𝑘for each probability under study. We plot those values in Figure 5.6 from which we derive two main conclusions. We observe that the values in Figure 5.6 follow a linear distribution. For each distribution we fit minima max𝑘values to a linear model and we derive the resulting correlation coefficient. The correlation coefficient quantifies the strength of the linear relation between two variables. It ranges between -1 and 1, with 1 or -1 indicating perfect correlation (all points would lie along a straight line). The distributions we use in this work include 4 types of unimodal distributions, 2 multi-modal distributions and several parameters thereof (see Table 3.2 in Section 5.5). For all those distributions, Table 5.1 shows that the correlation coefficient is very high and steadily stays above 0.99. Even the empirical distributions derived from the case study analyzed in Section 5.6 result in a high coefficient of correlation (0.98 on average), despite they tend to produce more variability in the estimations. Hence, empirical evidence in support of linearity for the minimum observed value of max𝑘is strong for the distribution tails and range of probabilities representative for the WCET. Besides, the application of RESTK includes a method to assess whether the linearity property holds, building on the observed data. It is also noted that, similar empirical reasoning is used to support 55 Chapter 5. Sky-high Quantile Estimation with Markov’s Inequality statistical arguments whether phenomena adhere to specific distributions builds on empirical tests. For instance, in the case of EVT, previous work uses QQ-plots to assess, based on observation, whether some tails can be considered exponential [105]. Table 5.1: Correlation Coefficient for all the distributions used in this work. Gauss1 Gauss2 Weib1 Weib2 Beta1 Beta2 Gam1 Gam2 Mix1 Mix2 Mix3 Mix4 .999 .999 .999 .999 .999 .998 .999 .999 .998 .998 .997 .998 Overall, by restricting max𝑘, one can avoid under-estimating the upper tail of the modeled distribution. This is exemplified in Figure 5.5, where the green squares (RESTK) represent the estimates obtained for the Gaussian distribution when applying the max𝑘values in Figure 5.6. By restricting max𝑘, we address the lack of trustworthiness in Figure 5.5 (NO RESTK) and produce tight and trustworthy bounds. Analogous results are obtained for the other distributions. 5.4.3 Deriving max𝑘from Unknown Distributions When deriving Figure 5.6, we built on the values of the theoretical quantiles so as to determine the value of 𝑘for which the sample moment starts underestimating. Given a sample of size 10𝑝, we can estimate confidently quantiles from exceedance probabilities bigger than 1/10(𝑝−1). In this case, based on the law of large numbers, it is very likely to see around 10 realizations whose probability is of the order of 10(𝑝−1)[151]. That is, on a sample of size 𝑛=1000 we will see around 10 realizations whose probability is 0.01 (1%). Therefore, for a sample of size 104, quantiles corresponding to exceedance probabilities 10−3and 10−2can be estimated easily with the usual quantile estimation functions from statistical packages [88]. The computation of confidence intervals for quantile estimation can be done using distribution-free methods like Kaigh and Lagenbruch or bootstrap [152]; and in any case the accuracy of the estimation can be increased using a bootstrap technique to correct variability. RESTK estimates at least three max𝑘points to construct its model and assess linearity. The latter is assessed by deriving the correlation coefficient for these three points. If such coefficient is above a threshold th =0.95, we deem max𝑘boundary to be linear and vice-versa (in which case RESTK cannot be applied). For instance, the quantiles corresponding to exceedance probabilities 10−3,10−4, and 10−5can be estimated very accurately with a sample of size 𝑛=106. These reference points allow us to assess when RESTK underestimates, and hence generate three max𝑘points, one for each probability. With those points, we can generate the regression line that projects max𝑘for any probability of interest (e.g., those in Figure 5.6) and assess that the correlation coefficient is above the desired threshold. Algorithm 2 generalizes the RESTK process starting from a main sample of the distribution under analysis (size 10𝑝), selecting the range of 𝑘to explore and the number of 𝑛sims to run. First, RESTK estimates the quantiles at given probabilities 𝑝test from the main sample (Line 3), e.g. 10−(𝑝−1),10−(𝑝−2), and10−(𝑝−3). For each simulation, 𝑛boot bootstrap samples of size 10𝑝−3are generated from the main sample (Line 5). For each of these samples, we compute the maximum 𝑘and maximum tightness for alltheprobabilitiestotest. Thepredictedvaluevpred isobtainedbycomputingMIKfromEquation5.2 with the sample moment estimator, ˆ 𝐸(𝑋𝑘)in Equation 5.4, (Line 10). The tightness of the predicted value is computed (Line 11) by using as reference value vref the estimated quantiles obtained before (Line3). The algorithm finds the values of 𝑘inthe considered range thatproduce the tightestestimate (Lines 15-16) and terminates its exploration as soon as a 𝑘that underestimates (tightness <1) is found (Line 12). After exploring all selected probabilities for all 𝑘and all simulations, the algorithm returns the smallest max𝑘across all simulations (Lines 20-21). Once the max𝑘(e.g. for 10−(𝑝−1),10−(𝑝−2), and 10−(𝑝−3)) are obtained, RESTK builds a linear model for all possible probabilities (Line 24). The final check (Line 25) will ensure that the correlation of max𝑘is above th =0.95, otherwise RESTK provides no pWCET estimate. As we show in Table 5.1, the correlation should always be close to 1 for max𝑘. The threshold th =0.95 is a standard and stringent threshold for confidence in statistics, and used as a way to discard estimates that do not meet the safety criteria of finding a proper max𝑘. RESTK enables the computation of a value for max𝑘that reduces the risk of underestimation for any probability of interest. This can be directly exploited by running MEMIK (Algorithm 1) on a predetermined range for 𝑘and for each probability 𝑝under study by using max𝑘(𝑝)as the upper bound for 56 5.5. RESTK and EVT PWCET Estimates on Distributions Algorithm 2 Computing the boundary necessary to apply RESTK approach. 1: function RESTKBoundary(sample,𝑘range,𝑘step,𝑛boot,𝑛sims,𝑝test,𝑝all,th) 2: for 𝑝∈𝑝test do 3: 𝑞est ←estimateQuantiles(sample, 𝑝) 4: for sim ∈𝑛sims do 5: sampleboot ←bootstrap(sample, 𝑛boot) 6: max𝑘(𝑝)←∞ 7: tightnessbest(𝑝)←∞ 8: for 𝑘∈𝑘range, 𝑘step do 9: vref ←𝑞est 10: vpred ←MIK(𝑘, sampleboot, 𝑝) 11: tightness ←vpred/vref 12: if tightness <1then 13: break 14: end if 15: if tightness <tightnessbest(𝑝)then 16: current𝑘(𝑝)=𝑘 17: end if 18: end for 19: if current𝑘(𝑝)<max𝑘(𝑝)then 20: max𝑘(𝑝)=current𝑘(𝑝) 21: end if 22: end for 23: end for 24: max𝑘(𝑝all) ← buildLinearModel(max𝑘(𝑝test), 𝑝all) 25: if corr(max𝑘(𝑝all)) <th then 26: return no pWCET estimate 27: end if 28: return max𝑘(𝑝all) 29: end function 𝑘, instead of considering an arbitrary range. Also, note that with RESTK, we use Equation 5.2 with sampled moments ( ˆ 𝐸(𝑋𝑘)) as EVAL(𝑡, 𝑘, 𝑝)function in MEMIK. The rest of the MEMIK algorithm remains unchanged when using the RESTK method. 5.5 RESTK and EVT PWCET Estimates on Distributions Our implementation of RESTK and MEMIK is programmed in 𝑅[133]. We run experiments on an Intel Core i5-7600K CPU clocked at 3.8GHz. The maximum execution time required per experiment was very short, 50 milliseconds or lower. We analyze values of 𝑘in the range 𝑘∈ [1,150]with 𝑘step =1to estimate max𝑘for all reference distributions, which, as shown in Figure 5.6 is a wide enough range to find the best max𝑘across distributions and probabilities. For all methods compared in this section, we use a sample size of 𝑛=106. For RESTK, we set the number of bootstrap simulations to 𝑛boot =2000. For EVT, we use the PoT methodology to fit tails and the CV Plot [52] to find the appropriate threshold for the PoT model. We use two EVT models fitted for pWCET estimation, namely exponential and light tails models. For each specific model, a different threshold using the CV plot will be found to ensure the best possible fit. •Exponential: with an exponential model, the shape of the GPD is fixed to 𝜉=0, which only leaves us to estimate the threshold 𝑢and the scale 𝜎. The threshold is estimated with the CV Plot fixing 𝜉=0, which finds where the exponential tail begins. Once we find the tail, we separate it from the rest of the sample and estimate the scale 𝜎with it. •GPD light tails: for the light tails model, we need the value of 𝜉, with 𝜉 < 0, best fitting the data. Using the CV Plot we find the threshold 𝑢where the light tail begins. Then, we separate 57 Chapter 6. Merging Hardware Event Monitor Data for Complex MPSoCs Using Order Statistics 0 20 40 60 80 100 120 140 160 CYCLES_LMQ_LOSES_DLINK_… STALL_FOR_WDB_CYCLES STALL_FOR_RLDB_CYCLES STALL_FOR_RLT_CYCLES CYCLES_SFX1_IDLE CYCLES_SFX0_IDLE STALL_FOR_DTQ_CYCLES CYCLES_CB_FULL_OR_CLOSE_… PROCESSOR_CYCLES CYCLES_0_INSTRUCTIONS_CO… CYCLES_DECODE_STALLED CYCLES_IB_FULL_OR_CLOSE_… CYCLES_CFX_SCHEDULE_STAL… BIU_MASTER_REQUESTS BIU_MASTER_DATA_SIDE_RE… BIU_MASTER_GLOBAL_REQU… L2_SNOOP_HITS L2_SNOOPS_CAUSING_MINT L2_RELOADS_FROM_CORENET L2_SNOOP_PUSHES L2_STORE_ALLOCATES CYCLES_A_PRE_SYNC_SERIALI… FPU_INPUT_DATA_STALL L2_DATA_MISSES CYCLES_SGB_LOSES_DLINK_A… CYCLES_VCFX_IDLE BLINK_REQUEST L2_MISSES L2_ACCESSES L2_DEMAND_ACCESSES L2_DATA_ACCESSES CYCLES_BU_IDLE CYCLES_FPU_IDLE CYCLES_VPU_IDLE CYCLES_VFPU_IDLE CYCLES_VSFX_IDLE 0 1 2 5 7 7 4 7 7 7 7 7 7 5 5 5 5 5 5 5 5 7 7 5 3 7 4 5 6 6 6 7 7 7 7 7 Variability (%) Figure 6.1: Observed variability for several HEMs in the T2080. values subject to different noise. HRM identifies an anchor HEM, and defines groups of HEMs, each group with PMCs - 1 HEMs plus the anchor. HRM performs several runs for each group of HEMs, ranks HEM values in each group using order statistics on the anchor HEM, and merges those HEMs with the same rank in different groups. Order statistics are non-parametric and hence, can handle the different distributions of the HEM values observed. Since noise-free HEM values cannot be obtained in general in complex MPSoCs, we evaluate HRM comparing the correlation across HEMs merged by HRM against their correlation when those HEMs are measured in the same run, thus under identical noise. Our results show that HRM captures accurately the correlation between HEMs, as opposed to blindly merging HEMs read in different runs. 6.2 Motivation We motivate this work by showing how several of the 262 HEMs in the T2080 present significant variation (Section 6.2.1) and follow different distributions (Section 6.2.2). We also dig down into some of the reasons behind the observed variation (Section 6.2.3). 6.2.1 HEM Variability On the NXP T2080 [69] we run a four-task workload with each task pinned to one of its e6500 cores. Each task performs integer and floating operations at the core level over several large vectors so that data operated is fetched from main memory, causing frequent misses in all cache levels, and thus exercising several HEMs. For this experiment, as well as the remaining ones throughout this chapter, we run on baremetal to remove potential interference coming from the operating system (the specifics of our experimental framework are described in Section 3.1.2). We divide the experiment into several sub-experiments, in each of which we read 6 HEMs (the total number of PMC available in the T2080). Hence, reading all 262 HEMs requires 44 sub-experiments, each of which we repeat 100 times to capture the impact of noise on HEM readings. In all runs we focus on the HEMs for core 0. Each run finishes when the task in core 0 finishes. Figure 6.1 shows the maximum relative variability observed for several HEMs, i.e. var =(max − min)/min, with bars in the figure sorted from higher to lower. Each bar is tagged (see bottom part of the figure) with the order of magnitude 𝑚of the value of each HEM in the experiment. For instance, for 𝑚=3,103≤hem𝑖<104. This information allows assessing the potential impact of the variability on execution time, whose magnitude for this experiment is tens of millions (107) of cycles. So is that for the number of committed instructions. In the NXP T2080, the maximum duration of an event triggered by an instruction can be in the order of hundred cycles (102). Hence, HEMs below 104arguably have low impact on performance. This, of course, is related to events not involving the execution of system software, e.g. a TLB miss, whose impact is not covered by multicore contention timing analysis but instead captured by the system-level timing analysis. We differentiate some cases for our experiments. 64 6.2. Motivation ECDF ECDF ECDFECDF ECDFECDF ECDF (a) (c) (b) (d) (e) Figure 6.2: Histogram and empirical CDF (ECDF) types: (a) Normal, (b) Concave, (c) Convex, (d) Clustered, (e) hard-to-fit. •Relevant high variability. Some HEMs present high variability while their magnitude is relevant, 104-107. These HEMs are the focus of our study as they can significantly impact the timing of the application and hence, the bounds that can be derived to it. In this category we find PROCESSOR_CYCLES with a variability of 45% from 3.6·107to 5.2·107(in other experiments the variability of this HEM reached 59%). 37 HEMs fall in this category if we set 1% as threshold for low-variability. •Irrelevant or low variability. Other HEMs have low variability in absolute terms, thus having little impact on performance. There are 70 HEMs in this category, including the three on the left of Figure 6.1 whose variability is over 180% but their value is below 300, hence, insignificant with respect to the cycle count. Other HEMs, 5 in total for this experiment, while exhibiting values above 104, incurred less than 1% variability, with limited impact on performance. •Not exercised. Finally, other HEMs, 150 in our case, were not exercised by the program under analysis making both the minimum and maximum value be zero. As the set of HEMs exercised can change across different experiments, the particular non-exercised HEMs will like vary. In fact, this is the motivation behind having different benchmarks in the experimental evaluation. The observed HEM variability does not depend on the particular subset of HEMs that are enable/disabled when collecting observations. Interestingly, the ‘low variability’ category comprises HEMs presenting no variability. Those are related to the functional execution of the program capturing the number of completed instructions including SFX, CFX, store, load, stores, taken and non-taken branches. For instance, the total number of executed instructions (INSTRUCTIONS_CMPLTD) of the first task is exactly the same in all runs (15,646,749). This leads us to conclude that the observed variability does not come from the software that always performs the same function (e.g. it traverses always the same execution path), and instead the variability is induced by the hardware. 65 Chapter 6. Merging Hardware Event Monitor Data for Complex MPSoCs Using Order Statistics 6.2.2 Distribution Focusing on the relevant high variability HEMs, we identified their variable behavior falls into five main classes of distribution. These types are depicted in Figure 6.2 that shows the histogram (bars) and the cumulative distribution function or CDF (line) of observed values for one HEM in each category for illustrative purposes. The x-axis shows HEM value, the left y-axis the frequency of occurrence for the histogram (for a 500 observations sample), and the right y-axis the fraction of observations for the CDF. a. Normal. HEMs in this category show a symmetric behavior that resembles a normal distribution. b. Concave. The distribution resembles a uniform with leaning towards the smallest values, which gives a concave cumulative distribution function. c. Convex. Distribution with the probability mass concentrated on the highest values of the distribution, giving a convex cumulative distribution function. d. Clustered. HEMs in this category show a clustered behavior around two values or more values. Distributions with more than one clear mode also fall into this category. e. Hard to fit. Finally, some distributions follow no obvious distribution (hard to fit category) apparently characterized by having two modes and a long tail. Out of the 37 relevant HEMs in our experiment, their distribution is as follows: 2 (5.4%) Normal, 13 (35.1%) concave, 19 (51.4%) convex, 1 (2.7%) clustered, and 2 (5.4%) hard-to-fit. The case of the HEM PROCESSOR_CYCLES is particularly relevant for a two-fold reason. It is the main HEM used in timing analysis for making predictions, and it presents a hard to fit distribution (see Figure 6.2(e)). This HEM presents 59% variability from around 2.0·107to 3.2·107. 6.2.3 Reasons Behind the Observed Variability The T2080 implements a complex architecture with an aggressive core (the e6500), so some form of hardware-induced HEM variability is therefore expected. We have observed that the HEMs with relevant high variability capture the activity in a wide range of hardware units, from the (on-core) integer issue queue to the internal queues of the L2 cache. High variability can be due to the complex nature of the T2080 and its sources of multicore interference: specific hardware scheduling choices in the multiple shared queues and buffers in the core-to-L2 interconnect, internal to the L2, the CoreNet Coherence Fabric, and the memory controller, may lead to variable latencies for specific requests. Specific and controlled execution scenarios allow narrowing down the sources of execution time variability. As an example, we have performed some bare-metal experiments where all cores hit L2 cache sustainedly, with a task 𝜏overlapping its full execution with the others. The intent is that interference occurs solely in the L2 cache. Variability observed across executions (up to more than 40%) could be attributed to minor initial processor state differences causing slight time shifts between L2 accesses across runs, and leading to cascade effects in L2 queues. In general, however, the limited information about the internal functioning of some of these resources, e.g. CCF, simply prevents identifying some of the reasons behind the observed variability. Also, as the programs used in this study typically perform the same activities repeatedly, contention for requests of a given core can stay repeatedly low or repeatedly high, leading to cumulatively high variability. Further, such systematic patterns may inadvertently switch from low to high contention scenarios (or vice versa) due to several reasons, such as the effects of loop control instructions in a program, which might alter systematic behavior inside the loop, as well as the impact of DRAM refresh operations, just to name some examples. Authors in [171] perform a hardware analysis of several Intel architectures and formulate several hypotheses on the reasons behind various forms of under and over counting affecting some HEMs (retired instructions, branches, load/stores, floating point, etc). Extending this to modern MPSoCs confronts with the inclusion of large hardware IP blocks with limited description and the increasing number of HEMs monitoring events highly sensitive to such variability. Also, note that knowing the reasons behind such variability would help assessing whether the device allows some configurations under which the variability reduces. However, if the behavior causing the variability is intrinsic to 66 6.3. Problem Formalization Table 6.1: Main terms used in this work. Term Definition nh,np,nb,nr number of HEMs, PMCs, sub-experiments, and runs ℎ𝑖HEM with id 𝑖 𝑏𝑗𝑗th sub-experiment 𝑟𝑗,𝑘 run 𝑘of a given sub-experiment 𝑏𝑗 𝑚𝑖 𝑗,𝑘 measured value for ℎ𝑖in run 𝑟𝑗,𝑘 𝑀𝑖 𝑗set of measured values for ℎ𝑖in 𝑏𝑗 𝑣𝑖 𝑗,𝑘 range of variation of measured value 𝑚𝑖 𝑗,𝑘 ℎ𝑎anchor HEM sr𝑗,𝑙 Run of 𝑏𝑗with the 𝑙th lowest value of ℎ𝑎 SR𝑙Concatenation of sr𝑗,𝑙 for all {𝑏𝑗}𝑗=1,···,nb SM𝑙HEM readings in sr𝑗,𝑙 for all {𝑏𝑗}𝑗=1,···,nb the complexity/functioning of the device, a solution like HRM is still needed – whether or not the root can be explained. 6.3 Problem Formalization Combining the discussion Section 2.4.1 about the disproportion between HEMs and PMCs with the analysis of high variability between measurements of the previous section, there lies a challenge to merge the information from different HEM readings into a single dataset as if they were measured together. In this section we will formalize the HEM merging problem. We are interested in collecting the values of a set of relevant HEMs, ℋ ∋ {ℎ1, ℎ2,··· , ℎnh}, whilst a given program executes on the target platform in response to a given input (the main terms used in this chapter are listed in Table 6.1). Inanidealscenario, allnh HEMsarecollectedatonceonasingleprogramexecutionwithoutincurring theuncontrolled(platformorsystemlevel)jitterorvariabilitythatmayariseacrossexecutions. Under such favorable conditions, we obtain a set of measurements (values) for each HEM in ℎ𝑖∈ {ℋ} that cumulatively capture the activity performed by the program. This is referred to as Scenario 1 in Figure 6.3, in which the row shows the single execution, columns the HEMs and the cell their respective values. In a more realistic scenario, program executions on the target platform are subject to noise so that in each execution the measured values for a given HEM ℎ𝑖can potentially vary. Note that we use the term ‘noise’ to generically refer to the varying execution conditions across experiments, either due to different initial hardware and system software state (in this respect, our experiments are executed baremetal reducing the variability due to system software). We are not after quantifying such noise, but we just recognize that it is in general uncontrollable, beyond the measures we take in order to reduce it as shown in Section 3.1.2. To capture the impact of noise, several runs of each experiment need to be performed. The noise of the different runs is represented as different levels of grey in Figure 6.3 (Scenario 2). In this scenario, noise can occur but at least all HEMs can be read at once, so all HEMs in each run are exposed to the same noise. This makes possible to reason about their relationships for statistical inference. In general, however, it is not possible to read all nh HEMs at once in a single execution as the number of HEMs that can be tracked simultaneously is determined by the number of available PMCs. Assuming our platform support np configurable PMCs1typically comparatively small with respect to the number of supported HEMs, with np ≪nh. For this reason, HEMs are necessarily collected in groups of at most np elements. Hence, to measure all HEMs for a given program we 1Without lack of generality, we assume there are no constraints on which specific HEM can be read from each PMC. Some processors exhibit such constraints, due to hierarchy of multiplexors to route HEM readings to a specific PMC. This scenario would just restrict which HEMs can be read in the same run, but would not affect HRM, as some HEMs (e.g. PROCESSOR_CYCLES) can be read along with any other group of HEMs. 67 Chapter 6. Merging Hardware Event Monitor Data for Complex MPSoCs Using Order Statistics Scenario 3: With Noise.  Only 2 HEMs can be read at once Scenario 1: No Noise.  All HEMs can be read at once run1 run1 Scenario 2: With Noise.  All HEMs can be read at once run2 runj-1 runj h1h2hnh-1 hnh ... ... ... ... ... ... ... ... run1 run2 runk-1 runk h1h2hnh-1 hnh ... ... ... ... ... ... ... ... ... Naïve merging (run1 and runk) Figure 6.3: Scenarios in HEM reading. must perform a set of at least nb ≥ ⌈nh/np⌉sub-experiment (𝑏1,··· , 𝑏𝑗,··· , 𝑏nb), each capturing the values of at most np distinct HEMs and cumulatively covering all HEMs. To capture the variability in measured values, several runs of the same sub-experiment 𝑏𝑗are carried out, see Scenario 3 in Figure 6.3, with crosses showing the HEMs not read in a given run. In this case, we assume only 2 HEMs can be read in each run. Also, as shown at the bottom of Scenario 3, naively merging HEMs (from the first run and kth run in this case) results in merging HEMs values obtained under different noise levels, potentially resulting in inconsistent values that cannot be reliably used. In this chapter, we address the challenge of merging the readings (measurements) for all ℎ𝑖∈ ℋ, each one measured several times in a different sub-experiment (Scenario 3) to obtain noise-consistent measurements for all HEMS (Scenario 2), preserving their relationships with execution time to favor timing analysis. Note that, noise-free HEM values (Scenario 1) are arguably hard to achieve, if at all possible, in MPSoCs. In particular, we aim at obtaining vectors with values for all HEMs under similar noise, as if all of them could have been read simultaneously in every single run. 6.4 HRM: a Technique to Merge HEMs Table 6.2 introduces an example with the main inputs and outputs to be generated by any HEM merging approach. In particular, it shows the measurements made when the number of HEMs is nh =12 and the number of PMCs is np =3, hence being required nb =4sub-experiments. In the example, nr =5runs areperformed persub-experiment. Onthe left, itis reportedthe sub-experiment and run id. In the top part, the HEM id. We use 𝑟𝑗,𝑘 to refer to the run 𝑘of sub-experiment 𝑏𝑗. We refer to measured value of HEM ℎ𝑖in 𝑏𝑗and run 𝑘as 𝑚𝑖 𝑗,𝑘. In terms of outputs, a HEM merging mechanism must aim at producing a list of nr all-HEM readings (vectors) where each vector includes all HEMs. Each of the nr measurements of each HEM is placed exactly in one of those nr vectors. This is illustrated at the bottom of Table 6.2, where each run of each sub-experiment is merged with another run from each other sub-experiment so that each run is represented exactly once in the merged result. For instance, in the example, the 1st run of the first sub-experiment (𝑚1 1,1,𝑚2 1,1,𝑚3 1,1) is merged with the 𝑥′th run of the 4th sub-experiment, and with one run of each other sub-experiment represented as the 𝑥th run of the 𝑗th sub-experiment. 6.4.1 Approach Our approach, HRM, builds on non-parametric order statistics, which allows relating random variables based on the order of the sampled values of the variables, regardless of their distributions. In particular, HRM aims at merging the HEM measurements from different sub-experiments in such a way that the noise experienced by the different measurements is as similar as possible. HRM must also allow merging HEMs regardless of the distribution of the data to be merged. Non-parametric 68 6.4. HRM: a Technique to Merge HEMs Table 6.2: Example with nh =12,np =3,nb =4, a nr =5. ℎ1ℎ2ℎ3··· ℎ𝑖··· ℎ10 ℎ11 ℎ12 𝑟1,1𝑚1 1,1𝑚2 1,1𝑚3 1,1 . . .. . .. . .. . . 𝑏1𝑟1,𝑘 𝑚1 1,𝑘 𝑚2 1,𝑘 𝑚3 1,𝑘 . . .. . .. . .. . . 𝑟1,5𝑚1 1,5𝑚2 1,5𝑚3 1,5 ... . . .. . . 𝑏𝑗𝑟𝑗,𝑘 ···𝑚𝑖 𝑗,𝑘 ··· . . .. . . ... 𝑟4,1𝑚10 4,1𝑚11 4,1𝑚12 4,1 . . .. . .. . .. . . 𝑏4𝑟4,𝑘 𝑚10 4,𝑘 𝑚11 4,𝑘 𝑚12 4,𝑘 . . .. . .. . .. . . 𝑟4,5𝑚10 4,5𝑚11 4,5𝑚12 4,5 ⇓ SM1𝑚1 1,1𝑚2 1,1𝑚3 1,1··· 𝑚𝑖 𝑗,𝑥 ···𝑚10 4,𝑥′𝑚11 4,𝑥′𝑚12 4,𝑥′ . . .. . .. . . SM𝑖𝑚1 1,𝑘 𝑚2 1,𝑘 𝑚3 1,𝑘 ··· 𝑚𝑖 𝑗,𝑦 ···𝑚10 4,𝑦′𝑚11 4,𝑥′𝑚12 4,𝑦′ . . .. . .. . . SMnr 𝑚1 1,5𝑚2 1,5𝑚3 1,5··· 𝑚𝑖 𝑗,𝑧 ···𝑚10 4,𝑧′𝑚11 4,𝑥′𝑚12 4,𝑧′ order statistics, which resort to the order of data regardless of their distribution, allow relating runs across measurements through the use of an ‘anchor’ HEM, referred to as ℎ𝑎, measured in all sub-experiments. HRM derives the relation between HEMs in different sub-experiments via their relation to ℎ𝑎. This is illustrated in the left side of Figure 6.4 that shows how the individual readings of ℎ𝑖and ℎ𝑗(𝑚𝑖and 𝑚𝑗respectively) from different sub-experiments are related to those of the ℎ𝑎in each sub-experiment referred to as. HRM provides the following properties. First, it preserves the distribution of each individual HEM. It also preserves the joint distribution between each HEM, ℎ𝑖, and the anchor HEM, ℎ𝑎(recall that the joint distribution between the HEMs read in the same sub-experiment is maintained). Finally, HRM estimates the most reliable joint distribution across HEMs in different sub-experiments. Next, we detail the procedure followed by HRM to provide those properties, followed by the mathematical foundation of the approach. 6.4.2 Procedure The application process of HRM includes four main steps. STEP ➀. HRM starts by selecting the anchor HEM, ℎ𝑎, that will be read in all sub-experiments. In each sub-experiment np −1PMCs are used to read different HEM. That is, from all available np PMC, HRM uses one of them in each sub-experiment 𝑏𝑗for the anchor, and the other np −1PMCs for other HEM. HRM approximates unobserved HEM relationships via their individual (observed) relationship with ℎ𝑎. Thus, the selection of ℎ𝑎is critically important as it determines how effective is HRM to merge HEMs for the problem under study. As the problem at hand relates to timing analysis, we chose ℎ𝑎to be as relevant as possible to timing. In the case of the T2080, execution time is measured via the HEM PROCESSOR_CYCLES, and hence ℎ𝑎=PROCESSOR_CYCLES. STEP ➁. After performing nr runs of each sub-experiment, HRM sorts the runs of each subexperiment by ℎ𝑎, from lowest to highest. As a result, each element in the sorted list for each 69 Chapter 6. Merging Hardware Event Monitor Data for Complex MPSoCs Using Order Statistics sub-experiment x sub-experiment y mx,1...nr imx,1...nr a my,1...nr a my,1...nr j HRM order statistics 1 2 3 ... nr-1 nr mimamami ... ... ... ... sub-experiment x y Figure 6.4: Introduction to the HRM approach. sub-experiment will indicate an order statistic, with the 𝑘th order statistic of a sample being its 𝑘th-lowest value. Each sub-experiment is characterized by a small fixed set of HEMs, limited by the number of PMCs available, np. Each sub-experiment in {𝑏𝑗}𝑗=1,···,nb is represented by a set of nr runs of dimension np, 𝑟𝑗,𝑘 :𝑚𝑎 𝑗,𝑘 , 𝑚(np−1)(𝑗−1)+1 𝑗,𝑘 ,··· , 𝑚(np−1)(𝑗−1)+(np−1) 𝑗,𝑘 . The selection of nr should be based on prior knowledge of ℎ𝑎random variable behavior. Without such knowledge, one must resort to nr ≥30, as this is the minimum size to estimate the main properties of a distribution through the central limit theorem. Runs in each sub-experiment and across them, should be designed to ensure that they are independent and identically distributed to enable the probabilistic reasoning on which HRM builds. To achieve this property, we empty the processor state between runs (see Section 3.1.2). We assess it by performing statistical independence and identical distribution tests (see Section 6.5.2). STEP ➂. Once all sub-experiments are sorted based on ℎ𝑎, we merge the different sub-experiments so that the 𝑘th measurement in the list for ℎ𝑖in a given sub-experiment is merged with the 𝑘th measurement of ℎ𝑗in another sub-experiment. Naturally, HEM measurements in the same run of the same sub-experiment remain at exactly the same position in the sorted list, so they remain together upon merging. Let sr𝑗,𝑙 be the run of sub-experiment 𝑏𝑗with the 𝑙th lowest value of ℎ𝑎.sr𝑗,𝑙 is defined as sr𝑗,𝑙 =𝑟𝑗,𝑘 where 𝑚𝑎 𝑗,𝑘 is the 𝑙-lowest value in the set {𝑚𝑎 𝑗,𝑘 }𝑘=1,···,nr. Finally, the concatenation produces a vector with completed representation of HEMs SR𝑙=(sr1,𝑙 ,··· ,sr𝑗,𝑙 ,··· ,srnb,𝑙)for each 𝑙=1,··· ,nr. Since the 𝑙-lowest value of 𝑚𝑎 𝑗,𝑘 is well-defined, we can assume that for each 𝑗=1,··· ,nb,𝑚𝑎 𝑗,𝑘 ≤𝑚𝑎 𝑗,𝑘+1 for all 𝑘=1,··· ,(nr −1). Therefore, 𝑚𝑎 𝑗,𝑘 is the (𝑘:nr)-order statistic of the sample of size nr, {𝑚𝑎 𝑗,𝑘 }𝑘=1,···,nr. HRM merges the values read in the same ordered run across all sub-experiments, referred to as SM𝑙. For each SM𝑙, HRM produces one reading for each HEM and nb readings for ℎ𝑎. STEP ➃. After merging, we compute the summarized order statistics for ℎ𝑎. In particular, we compute the quantiles of the distribution of all values of ℎ𝑎across all sub-experiments so that we obtain exactly nr quantiles, i.e. one for each row of our merged list of HEM values. The resulting array yields ˆ 𝑚𝑎=quantile(0,··· , 𝑘/(nr −1),··· ,1), where 𝑘=0,··· ,(nr −1). HRM estimates the nr equal spaced quantiles of ℎ𝑎using the sample of size (nr ·nb)obtained from joining all ℎ𝑎values. 6.4.3 Quantile Estimation Several methods for quantile estimation can be considered. Let {𝑥(𝑘)}𝑘=1,···,𝑛 be an ordered sample of size 𝑛. In general, a method for quantile estimation corresponds to weighted averages of consecutive order statistics. Given fixed values for a function 𝛾and a constant 𝑚, the 𝑝-quantile is defined by 𝑞(𝑝)=(1−𝛾(𝑗, 𝑚))𝑥(𝑗)+𝛾(𝑗, 𝑚)𝑥(𝑗+1), where (𝑗−𝑚)/𝑛≤𝑝<(𝑗−𝑚+1)/𝑛,𝑥(𝑗)is the (𝑗:𝑛)−order statistic. We consider a continuous representation of quantile estimation with 𝛾(𝑗, 𝑚)=𝑝·𝑛+𝑚−𝑗 and 𝑚=1−𝑝, which is equivalent to do linear interpolation between the points {(𝑝𝑘, 𝑥(𝑘))} where 𝑝𝑘attempts to estimate the mode of 𝐹(𝑥(𝑘)). Then 𝑞(𝑝)is a continuous function of 𝑝and 𝑝(𝑘)= (𝑘−1)/(𝑛−1). We refer the interested reader to [88] for a review of programming quantile estimation. 70 6.4. HRM: a Technique to Merge HEMs 6.4.4 Correlation Boundary HRM produces a solution that preserves the observed information and reliably builds unobserved information by preserving joint distributions. That is, HRM preserves the correlation across HEMs. In particular, HRM describes the relationship between the expected value of a target HEM, the anchor ℎ𝑎, and the values observed for a different HEM ℎ𝑖. The relationship across different HEMs, namely ℎ𝑖and ℎ𝑗, not observed together, is built therefore through ℎ𝑎. HRM builds such relationship by estimating the covariance matrix across all HEMs. The covariance matrix can be used because each ℎ𝑖is a random variable with at least 4 finite moments. In our case, each HEM has infinite finite moments since all HEM values are bounded, i.e. they count finite events per cycle during a finite interval, since the measurement starts until the measured value is collected. Therefore, each HEM value is a bounded number, thus guaranteeing the existence of infinite finite moments. By being random variables with finite moments, we can describe the relationship across expected values of HEMs through a multivariate normal distribution based on the central limit theorem, which ultimately ensures the existence of the covariance matrix that characterizes the relationship between the expected values of HEMs asymptotically. In particular, correlation across HEMs is based on Pearson correlation coefficient [66]. It can be obtained via a Least-Squares fit, where a value 1 represents a perfect positive relationship, −1a perfect negative relationship, and 0the absence of any apparent relationship across variables. Let 𝑋and 𝑌be random variables, and denote by cor(𝑋, 𝑌)the Pearson correlation, obtained as cov(𝑋,𝑌) 𝜎𝑥𝜎𝑦, where 𝜎and cov describe the variances and covariance, respectively. Lemma. Let (𝑌, 𝑋1, 𝑋2)be a random vector with multivariate standardized normal distribution. Then, the correlation between 𝑋1and 𝑋2is in the interval 𝜌1𝜌2±q1−𝜌2 1q1−𝜌2 2, where 𝜌𝑖=cor(𝑌, 𝑋𝑖)for 𝑖=1,2. Proof. Let Σbe the covariance matrix the joint distribution between the random variables, 𝑌,𝑋1and 𝑋2.Σcan be described as the correlation matrix ©« 1𝜌1𝜌2 𝜌11𝜌 𝜌2𝜌1ª®¬.(6.1) where 𝜌cannot be arbitrarily set between [−1,1], since the matrix must be positive semidefinite. A Hermitian matrix is positive semidefinite if and only if all principal minors are non-negative. Building on Silvester’s criterium, only the minors defined by submatrices starting from the upper left corner need being checked. The 2-by-2 submatrix, 1𝜌1 𝜌11, is trivial, and the 3-by-3 matrix produces the result to prove. Note that, by operating the result of the multivariate standardized normal distribution with the corresponding 𝜇and 𝜎of the random variables of the HEMs whose joint distribution we are studying, we can directly obtain the result for the multivariate non-standardized normal distribution. Moreover, since the random variables studied (the HEMs) have 4 finite moments (in fact they have infinite moments), and based on the central limit theorem, the Lemma guarantees that observed values converge asymptotically to the expected values. Building on the Lemma, we can prove that HRM guarantees the three properties described in Section 6.4.1. Theorem Let 𝐻be the joint distribution of all HEMs, and assume the set of ordered sub-experiments as shown in the example in Table 6.2. Let be {ˆ 𝑚𝑎 𝑘}𝑘=1,···,nr the set of nr equal spaced quantiles of ℎ𝑎 from the sample {𝑚𝑎 𝑗,𝑘 }𝑗,𝑘 of size (nb ·nr). Consider the complete (merged) vector with all HEMs defined as: 71 Chapter 6. Merging Hardware Event Monitor Data for Complex MPSoCs Using Order Statistics Table 6.3: HEMs with observed relevant variability. Name ID CYCLES_LSU_SCHE_STALLED 1 CYCLES_LSU_ISSUE_STALLED 2 BLINK_REQUEST 3 L2_MISSES 4 L2_DEMAND_ACCESSES 5 L2_ACCESSES 6 L2_STORE_ALLOCATES 7 L2_DATA_MISSES 8 L2_RELOADS_FROM_CORENET 9 L2_SNOOP_HITS 10 L2_SNOOP_PUSHES 11 STALL_FOR_RLT_CYCLES 12 STALL_FOR_WDB_CYCLES 13 BIU_MASTER_REQUESTS 14 BIU_GLOBAL_REQUESTS 15 ˆ 𝑚𝑎 𝑘, 𝑚1 1,𝑘 ,··· , 𝑚(np−1) 1,𝑘 ,··· 𝑚(np−1)(𝑗−1)+1 𝑗,𝑘 ,··· , 𝑚(np−1)(𝑗−1)+(np−1) 𝑗,𝑘 ,··· 𝑚(np−1)(nb−1)+1 nb,𝑘 ,··· , 𝑚(np−1)(nb−1)+(np−1) nb,𝑘 , for each 𝑘=1,··· ,nr. Then, the empirical joint distribution described by the complete vectors complies the (formalized) properties of HRM: •Property 1. Preserves the marginal distribution of 𝐻for all HEMs. •Property 2. Preserves the joint distribution across HEMs in the same sub-experiment. •Property 3. Estimates, with minimum error on correlation, the joint distribution between HEMs in different sub-experiments. Proof.Property 1: The marginal distribution of all HEMs but ℎ𝑎is preserved, since no modifications are produced in the observed values of those HEMs – they are just sorted. Only ℎ𝑎is modified, since it is replaced by the order statistics. As its distribution has infinite moments, replacing ℎ𝑎by its order statistics of a larger sample, leads to a higher amount of information (i.e. nb ·nr values instead of nr), improving the sample and hence, preserving the marginal distribution of ℎ𝑎. Property 2: The re-ordering procedure preserves measurements of different HEMs in the same subexperiment together. Therefore, their joint distribution is preserved identical. Property 3: Regarding the joint distribution of HEMs in different sub-experiments, HRM estimates such jointdistributionforeach pairof HEMs. Notethat, sincethose pairsof HEMsare never observed in the same sub-experiment, data collected provides no information about their joint distribution. The only relation across sub-experiments is had through ℎ𝑎, which is observed in all of them, so the joint distribution to be estimated needs to preserve this common relation. Based on the Lemma, such relation is preserved if the estimated correlation for describing the real joint distribution of two HEMs in different sub-experiments is in the interval 𝜌1𝜌2±q1−𝜌2 1q1−𝜌2 2, where 𝜌1and 𝜌2are the correlation between each of those two HEMs and ℎ𝑎. Since there is no additional information about the actual correlation between those two HEMs, any value in the interval is equally probable. Thus, the correlation value proposed by HRM is 𝜌1𝜌2, since this is the value that minimizes the absolute error with respect to the real value. 72 6.5. Experimental Evaluation 6.4.5 Matrix Completion Techniques HRM aims at merging actual observations rather than filling missing values with synthetic data. The latter, which may be realized with Matrix Completion methods [83, 136], as discussed before in Section 2.4.2, is not appropriate in our case. This is so because Matrix Completion requires that values in each row and column belong to a different distribution, which is not our case, since each column is a different HEM with its own distribution. As a consequence, the use of Matrix Completion methods for our problem leads to inadequate value distributions where, for instance, the mean and standard distribution of the synthetic data for all HEMs is extremely different from those for actual observations. For instance, in our experiments, the mean for synthetic data is ≈20𝑥smaller than that of real data, whereas the standard deviation is between 0.36𝑥and 5𝑥that of real data. 6.5 Experimental Evaluation 6.5.1 Validation Methodology For the experimental evaluation of this work we resort to the microbenchmarks run in the NXP T2080 Reference Board described in Section 3.1.2 and represented in Figure 3.2. The benchmarks we use to generate validation data are the microbenchmarks described in Section 3.1.2, which are designed to stress different cache levels to conform a representative set of cases. In order to create a realistic scenario for the validation, we use the set of 16 workloads described in Table 3.1, where we run one benchmark per core and the readings are performed from core0. The validation of any HEM merging methodology is complex on real hardware as we do not have the noise-free value for each HEM (ℎ𝑖), as explained in Section 6.3. This prevents us from directly comparing the estimated value for each HEM with its corresponding noise-free value. Thus, we can only evaluate HRM comparing the correlation for the HEM merged with HRM against the real – measured – correlation. In order to evaluate estimated and real correlations, first, for each workload, we perform 100 runs for each of the 53 =261/5sub-experiments. Hence, we collect readings for all 262 HEMs with 5 HEMs plus ℎ𝑎read in each sub-experiment2, except the last group (sub-experiment) that only includes 1 HEM and the ℎ𝑎. We validate HRM for 15 HEMs having high relevant variability, see Table 6.3. To that end we on purpose place those HEMs in different groups so that their mutual correlation is not observed in the data used for HRM. For each of the 120 pairs3of HEMs we estimate their correlation ˆ 𝜌𝑖,𝑗, after merging them with HRM. We also collect 100 runs for a set of experiments in which those 15 HEMs and the anchor are observed in the same group. Thus, for each pair of HEMs (ℎ𝑖and ℎ𝑗), as well as ℎ𝑎, we obtain their actual correlation 𝜌𝑖,𝑗 from those measurements. This allows us comparing their real correlation 𝜌𝑖,𝑗 with the estimated correlation after merging with HRM ˆ 𝜌𝑖,𝑗. In particular we measure the absolute distance (difference) between |𝜌𝑖,𝑗 −ˆ 𝜌𝑖,𝑗 |, so that the maximum difference obtained for a pair of HEMs is 2. This happens when the estimated correlation is 1(or −1) and the real one is −1(or 1). 6.5.2 Independence and Identical Distribution HRM builds on these statistical properties for ℎ𝑎,PROCESSOR_CYCLES, to apply order statistics. In practice, this holds since all values of ℎ𝑎have been collected from the repeated execution of the same workload, withthesameinputs,andenforcingthesamehardwareandsoftwarestateasmuchasitcan be controlled. We have further evaluated these properties quantitatively. We performed an ANOVA test [63] to assess identical distribution of PROCESSOR_CYCLES across sub-experiments. The result of the test is a p-value 𝑝=0.57, so the test is not rejected comparing the law on the expected value of PROCESSOR_CYCLES, and tells us that the noise is identically distributed across sub-experiments. We assess independence within each sub-experiment with a Ljung-Box test [108] with lag =10. With a significance level 𝛼=0.05, independence is not rejected in 96% of the sub-experiments, so 2For the sake of convenience, we refer to the HEM read in the same sub-experiment as being in the same (HEM) group. 3All possible pairs with the 16 HEMs analyzed (the 15 relevant and ℎ𝑎). 73 Chapter 7. Merging Hardware Event Monitor Data for Complex MPSoCs Using Copulas ℎ1ℎ2ℎ3ℎ4ℎ5ℎ6ℎ7 𝑟1 𝑟2 𝑟3 𝑟4 𝑟5 𝑟6 𝑟7 𝑟8 𝑟9 𝑟10 𝑟11 𝑟12 𝑟13 𝑟14 𝑟15 𝑟16 𝑟17 𝑟18 𝑟19 𝑟20 𝑟21 𝑟22 𝑟23 𝑟24 𝑟25 𝑟26 𝑟27 𝑟28 Figure 7.3: Runs collected with MUCH when np <nh. In particular, np =3and nh =7in the example. •For each ℎ𝑖, the measurements collected in each sub-experiment allow obtaining a grouped sample, that is {ℎ𝑖}, which includes all measurements of ℎ𝑖across all sub-experiments. •For each ℎ𝑖, the values of ˆ 𝜇𝑖and ˆ 𝜎𝑖obtained from {ℎ𝑖}. •For each pair of HEMs ℎ𝑖and ℎ𝑗, we obtain ˆ 𝜌ij, which is the empirical counterpart of 𝜌ij.ˆ 𝜌ij is obtained from the sub-experiment where both ℎ𝑖and ℎ𝑗are observed together. Therefore, since we have ˆ 𝜎𝑖,ˆ 𝜎𝑗and ˆ 𝜌ij, we can also obtain their empirical covariance ˆ 𝜎ij. Then, we can produce the empirical correlation matrix ˆ 𝑆∼𝑆, and the corresponding (empirical) covariance matrix ˆ Σ. Once we have constructed this information, we have a model for how each HEM in the experiment relates to all other HEMs in terms of expected values. Since we have ˆ 𝜇and ˆ Σfor all HEMs, we can describe the corresponding MVGD as follows: 𝑋∼ 𝒩nh ˆ 𝜇,ˆ Σ. However, now we need to generate actual merged HEM vectors using measured data in accordance with the MVGD model. From a mathematical point of view, the challenge is to reorder the grouped sample {ℎ𝑖}for each HEM ℎ𝑖, such that for each pair of HEMs ℎ𝑖and ℎ𝑗the empirical correlation of the grouped samples is close to 𝜌ij. Equivalently, the challenge is to reorder the grouped samples such that the new corresponding (empirical) covariance matrix is close to ˆ Σ. This corresponds to an optimization problem with an existing solution, since the set of potential orders (all permutations) is large but finite. However, it is not a trivial problem from a computational point of view. We construct the solution with a probabilistic approach using the preliminary results of copula theory [122]. Application of copula theory. For each ℎ𝑖, the empirical distribution function 𝐹𝑖 emp lets us transform the grouped sample into a uniform sample. A uniform sample can be transformed into a Gaussian sample by applying the inverse function of the cumulative distribution function of a standard Gaussian distribution, Φ. Therefore, the grouped sample is transformed one-to-one into a standard Gaussian distribution: Φ−1𝐹𝑖 emp {ℎ𝑖}.(7.1) 80 7.2. MUCH: Multi-Correlation HEM Reading and Merging For instance, if {ℎ𝑖}is {1,2,5,8}, this would map to a uniform sample {0.2,0.4,0.6,0.8}1. Then, those values would map to a standard Gaussian distribution as {−0.842,−0.253,0.253,0.842}2. We apply this process to all HEMs, thus obtaining for each measurement of each HEM itscounterpart value forthestandardGaussian distribution. Then, we refer toas ˆ Σ0tothecovariance matrix obtained from the transformed data to differentiate it from the one obtained from the measured data (ˆ Σ). Finally, we generate a sample sampMVGD of 𝑛MVGD runs (e.g. 𝑛MVGD =10000) from the MVGD, which we define using ˆ Σ0: 𝑋∼ 𝒩nh 0,ˆ Σ0,(7.2) which produces a joint sample with marginal standard Gaussian distribution (i.e. sampled values follow such distribution). At this point, sampMVGD provides a matrix with as many columns as HEMs (nh), as many rows as HEM vectors we want to generate (e.g. as many as total measurements per HEM), and preserving the correlations across all pairs of HEMs (𝜌ij) simultaneously. However, values in the matrix correspond to a standard Gaussian distribution instead of being HEM values read. Thus, we use the actual sampMVGD to produce the indexes for order statistics. In other words, if for a given HEM ℎ𝑖,sampMVGD has a particular order of values (e.g. 𝑘th lowest first, 𝑙th lowest second, 𝑚th lowest third, and so on and so forth), we set the actual observed values for ℎ𝑖in the very same order to generate the merged HEM vectors. The easiest way to do this is setting 𝑛MVGD to the actual number of measured values per HEM. For instance, if for a given HEM ℎ𝑖the sampMVGD has produced the values {1.121,−0.870,−0.172,0.343}, and the actual values observed are 9,10,12,17, they would be sorted as follows: {17,9,10,12}, thus preserving the same ordering, but this time using the actual values read. By following their sampMVGD ordering for all HEMs to organize the actual values measured, we generate as many merged HEM vectors as actual values have been observed for each HEM. For instance, recalling the example in Figure 7.3, where we have 12 values per HEM, we could set 𝑛MVGD =12, and would sort the values for each of the 7 HEMs in the same order as their synthetic values in sampMVGD. In fact, once this process is complete, we could assess ˆ 𝜌′ ij for all pairs of HEMs in the merged vectors and compare them with the original values ˆ 𝜌ij obtained from pairwise HEM measurements. Some (small) discrepancy is expected due to statistical reasons (i.e. sampling processes can always produce inaccuracies). Such discrepancy could be reduced with an iterative process where 𝑋in Equation 7.2 is obtained as many times as needed and measured HEM values sorted accordingly to obtain new merged vectors where ˆ 𝜌′′ ij is compared to ˆ 𝜌ij. This process could be repeated a fixed number of times or until a specific criterion is fulfilled. However, this step is purely optional. Recalling the example for HRM limitations before, MUCH would successfully preserve all correlations acrossthe three HEMssimultaneously, overcoming the limitation ofHRM. In particular, MUCH does not favor any particular random variable (HEM) when merging and, instead, all correlations are preserved as accurately as possible at the same time. Instead, HRM strategy is a supervised one where correlations with a particular variable (anchor HEM) are perfectly preserved at the expense of causing large inaccuracies for other correlations if they are weakly correlated with the anchor. 7.2.3 Procedure For the sake of completion, we provide the application process of MUCH, which includes five main steps. The procedure can be followed visually on Figure 7.2. •STEP ➀. The HEM selection for each sub-experiment does not play a role in MUCH. Each HEM will be measured with every other HEM at least once, to capture the relation between them, i.e. to compute the empirical correlation matrix ˆ 𝑆. While generating the minimum number of sub-experiments allowing to capture all pairwise HEM correlations is convenient, it is not strictly mandatory for the application of MUCH, so combinations can be generated with greedy 1Given a HEM ℎ𝑖for which we have 𝑛values, the uniform sample probability space is split into 𝑛+1identical parts. Out of the 𝑛+2boundary values, we exclude 0and 1since they cannot be used later for the Gaussian distribution as their counterpart values would be −∞and +∞respectively. 2As for the uniform distribution, those values distribute the probability space into 𝑛+1parts with identical accumulated probability. 81 Chapter 7. Merging Hardware Event Monitor Data for Complex MPSoCs Using Copulas Table 7.1: HEMs with observed relevant variability. Name HRM id PROCESSOR_CYCLES ℎ𝑎 CYCLES_LSU_SCHE_STALLED 1 CYCLES_LSU_ISSUE_STALLED 1 BLINK_REQUEST 1 L2_MISSES 1 L2_DEMAND_ACCESSES 1 L2_ACCESSES 2 L2_STORE_ALLOCATES 2 L2_DATA_MISSES 2 L2_RELOADS_FROM_CORENET 2 L2_SNOOP_HITS 2 L2_SNOOP_PUSHES 3 STALL_FOR_RLT_CYCLES 3 STALL_FOR_WDB_CYCLES 3 BIU_MASTER_REQUESTS 3 BIU_GLOBAL_REQUESTS 3 algorithms if needed. For each sub-experiment at least 30 runs (nr ≥30) are needed to allow the use of the Central Limit Theory [87]. In general, the higher nr, the more accurate ˆ 𝑆will be. As a matter of fact, in this work we set nr =50. •STEP ➁. Once the values are gathered, map them to a standard MVGD, and compute the covariance matrix ˆ Σ0. •STEP ➂. Compute the MVGD using ˆ Σ0as shown in Equation 7.2, and generate a sample equal in size (𝑛MVGD) to the number of collected values for each HEM in a sub-experiment or in all sub-experiments. Note that the method could also be applied with larger 𝑛MVGD values. •STEP ➃. (OPTIONAL) As an optimization step, we can compute the correlation matrix of the generated sample ˆ 𝑆′from the MVGD and compare it to the measured correlation matrix ˆ 𝑆, for instance, obtaining the Minimum Square Error (MSE). Then, we can repeat step ➂and keep the ˆ 𝑆′with lowest MSE compared to ˆ 𝑆. Such process can be repeated as many times as wanted as a way to further increase accuracy without requiring additional runs on the target platform. •STEP ➄. Copy the order statistics of the sample of the MVGD into the experimental data. Now, the experimental data is finally merged in accordance with the MVGD. 7.3 Evaluation This section presents the experimental framework, the validation approach followed to evaluate MUCH and compare it with HRM, and the results of the evaluation. 7.3.1 Validation Approach In this work the experimental validation is performed on the same platform, the NXP T2080, as with HRM in Chapter 6 with the same microbenchmarks and workloads described in Section 3.1.2. The reference against which to compare MUCH is the actual correlation of each pair of HEMs when measured together in the platform, so that we can validate whether correlations in the merged HEM vector are accurate with respect to real correlations. For completeness, we compare MUCH against HRM in terms of both, accuracy with respect to the real correlations and number of runs required. For the sake of comparison, we focus on the same 16 HEMs regarded as relevant in HRM [166], which we list in Table 7.1 for completeness. Note that, while all HEMs are treated homogeneously by MUCH – thus meaning that all pairwise HEM correlations are measured and then processed together in an identical basis – the same does not apply to HRM. In particular, HRM needs a HEM to be the anchor. Then, given that the T2080 MPSoC has 6 PMCs and one of them is used by the anchor, the remaining 15 HEMs need to be distributed across 3 sub-experiments. The sub-experiment where each HEM is read for HRM is shown in Table 7.1 in the HRM id column. 82 7.3. Evaluation Figure 7.4: Correlation Difference for each HEMs pair for workload 1 (top-left), workload 2 (top-right), workload 3 (bottom-left), and workload 5 (bottom-right). Note that, by using 16 HEMs, there are 120 different pairs of HEMs for which we evaluate the actual correlation obtained for the merged HEM vectors with MUCH and HRM, and compare them against their real empirical correlation, ˆ 𝜌𝑖,𝑗. In particular, we compute the accuracy for both methods as 𝜌𝑖,𝑗 much −ˆ 𝜌𝑖,𝑗and 𝜌𝑖,𝑗 hrm −ˆ 𝜌𝑖,𝑗. 7.3.2 Results We evaluate the correlation accuracy obtained, for both MUCH and HRM, against the real empirical correlation. For this first comparison, we set nr =50. The number of sub-experiments is only 3 for HRM (nb =3). For MUCH, while theoretically we could observe 120 pairs of HEMs with nb =8 sub-experiments (15 pairs per sub-experiment with 6 PMCs), we needed nb =10 just following a greedy process to create sub-experiments where we iterate over HEMs from 1 to 16, and for each one we create sub-experiments adding the lowest order HEM not yet observed with any of already selected HEMs in the sub-experiment. Therefore, 𝑛=150 for HRM and 𝑛=1000 for MUCH. We show detailed results for the 120 pairs of HEMs in Figure 7.4 for workloads W1, W2, W3, and W5. In particular, W2 corresponds to the case where the improvement of MUCH with respect to HRM is only moderate, W3 to an extreme case with huge improvement, and W1 and W5 to two cases with typical high improvement. The HEM pair values are sorted from lowest difference to highest difference with respect to the real correlation for each technique. As shown, MUCH provides higher accuracy since its differences with respect to the real correlation are much lower than those of HRM. Moreover, the difference in the worst case for MUCH is up to 0.25 for very few pairs of HEMs, whereas HRM reaches values above 0.5 for a non-negligible number of pairs, and even above 0.75 in some cases. Note that the maximum theoretical difference is 2.0, which would occur when the 83 Chapter 7. Merging Hardware Event Monitor Data for Complex MPSoCs Using Copulas MUCHMUCH n n n n Figure 7.5: Mean square error for HRM and MUCH as a function of the number of runs. Workload 1 (top-left), workload 2 (top-right), workload 3 (bottom-left), and workload 5 (bottom-right). estimated correlation is 1.0 (or -1.0), and the real correlation is -1.0 (or 1.0). As a second comparison, we study the dependence of each method on the total number of runs 𝑛=nr ·nb, which illustrates the trade-off between cost (in terms of number of runs) and accuracy for both methods. Again, we consider the same 4 workloads, where we vary nr (values 50, 100, 200 and 400), and obtain for each workload, the MSE for their difference with respect to the real correlation across the 120 pairs of HEMs. In particular, to produce a statistically significant comparison, we bootstrapped 50 samples for each of the methods and each nr value. For instance, for nr =100 this implies that we generate a random sample of 100 runs for each sub-experiment and apply the corresponding method on that sample. Then, we repeat the process 50 times, thus obtaining 50 estimates for each method and nr value. Figure 7.5 shows those results, where dots indicate individual measurements and the line corresponds to the mean across them. Note that both axes are in logarithmic scale. First, we observe that HRM obtains negligible gains from increasing 𝑛. Those gains are only noticeable for W2, where pairwise correlations with the anchor HEM are indeed high, allowing HRM to be almost as accurate as MUCH. In any case, HRM apparently plateaus at 𝑛=1200 (nr =400). For MUCH, we observe significant gains in all cases but W3, where increasing 𝑛produces limited improvements in accuracy. However, in the other 3 cases we observe improved accuracy as we increase 𝑛, and, apparently, such improvement does not plateau even with nr =400, thus offering opportunities to further increase accuracy if the number of runs is increased beyond that number. When comparing MUCH and HRM, we note that in general, MUCH allows reaching much more accurate merged HEM vectors. It is of prominent importance the case where 𝑛≈1000, because 84 7.4. Summary Table 7.2: MSE for the 16 workloads under iso-runs (𝑛=2100). Workload MUCH HRM W1 0.003 0.135 W2 0.004 0.024 W3 0.004 0.295 W4 0.006 0.272 W5 0.005 0.110 W6 0.007 0.083 W7 0.006 0.021 W8 0.006 0.068 W9 0.020 0.228 W10 0.007 0.168 W11 0.017 0.207 W12 0.005 0.158 W13 0.008 0.282 W14 0.005 0.106 W15 0.007 0.068 W16 0.017 0.250 it allows performing an iso-cost comparison across both methods, i.e. with the same number of total runs. In this case, where 𝑛much =1000 (nrmuch =50) and 𝑛hrm =1200 (nrhrm =400), thus with a slight advantage for HRM, we observe that MUCH is significantly better than HRM in 3 out of 4 workloads, and in the remaining one, where the MSE is already pretty low for both methods, MUCH is slightly better despite its slightly lower number of runs. Finally, note that in both methods, increasing 𝑛reduces dispersion of the bootstrap, thus making merged HEM vectors more stable in terms of accuracy with respect to the real correlations. For completion, we have evaluated all workloads with 𝑛=2100, so nrhrm =700 and nrmuch =210, again with a bootstrap with 50 samples. The mean MSE across the 50 samples is shown in Table 7.2. As expected, MUCH is systematically more accurate than HRM for all workloads, and the difference across both methods is tiny (e.g. W2 and W7) only when HRM is highly accurate, since MUCH is always highly accurate. In fact, the worst MSE for MUCH (0.020 for W9) is indeed better than the best MSE for HRM (0.021 for W7). So far we have shown that MUCH outperforms HRM under iso-cost (identical number of runs 𝑛), and naturally, under identical number of runs per sub-experiment (iso-nr), where MUCH has a higher 𝑛value. Moreover, we have shown that in all cases MUCH is highly accurate. However, the number of sub-experiments needed by MUCH is much more dependent on nh and np than that of HRM. For instance, if nh is high or np is very low, MUCH may need many sub-experiments whereas HRM only needs 𝑛hrm =nr ·nh−1 np−1. In particular, for MUCH we need to observe nh 2=nh·(nh−1) 2pairs of HEMs, and each sub-experiment allows observing up to np 2=np·(np−1) 2pairs. Assuming that sub-experiments for MUCH can be optimized to generate always unobserved pairs of HEMs only, the number of sub-experiments would be ratio between both values, and hence, the total number of runs would be 𝑛much =nr ·nh·(nh−1) np·(np−1). Figure 7.6 shows 𝑛much and 𝑛hrm for the case np =6, as in the T2080, and nr =100, when varying nh. As shown, 𝑛much grows much faster than 𝑛hrm as we increase nh. Therefore, there may be cases where 𝑛much might not be affordable and the only affordable solution is HRM. In those cases, despite the limitations of HRM, such solution has been shown to be systematically better than any other alternative (except MUCH) [166], and thus, it would be the best choice. 7.4 Summary In the previous chapter, we proposed the HRM approach to merge HEMs from different runs accurately preserving their correlation with respect to one anchor HEM (i.e. processor cycles) building on order statistics. However, HRM does not always preserve the correlation between other pairs of 85 Chapter 7. Merging Hardware Event Monitor Data for Complex MPSoCs Using Copulas Figure 7.6: Number of total runs as a function of the number of HEMs to arrange. The parameters on 𝑛much and 𝑛hrm are np =6,nr =100. HEMs that might be lost to a large extent. This chapter copes with HRM limitations by proposing the MUlti-Correlation HEM reading and merging approach (MUCH). MUCH builds on multivariate Gaussian distributions to merge HEMs from different runs while preserving pairwise correlations across each individual pair of HEMs simultaneously. Our results on an NXP T2080 MPSoC used for avionics systems show that MUCH largely outperforms HRM for an identical number of input runs. However, MUCH requires significantly more data to achieve such performance. Therefore, both HRM and MUCH have their use depending on the needs and resources of the engineer. 86 Part IV Contention Modelling for CRTES 87 Chapter 8 Clean Execution Times for MBPTA 8.1 Introduction Timing validation for CRTES (e.g. automotive systems) occurs in late integration stages when it is hard to control how the instances of software tasks overlap in time. To make things worse, in complex software systems, like those for autonomous driving, tasks schedule has a strong eventdriven nature, which further complicates relating those task-overlapping scenarios (TOS), which produce contention, captured during the software timing budgeting and those observed during validation phases. Multiple approaches exist to account for contention due to concurrency on the access to shared hardware resources, as part of the design and verification process of CRTES. While it is not the purpose of this thesis surveying on this topic, we identify the main families of techniques. Some techniques control the impact of contention by building upon some hardware and/or software support [93, 126, 183]. Other techniques use such support to completely avoid any impact due to contention [8, 17, 25, 129, 130]. Finally, some other approaches, rather than controlling or mitigating contention, aim at upper-bounding it [47, 145]. While each of those techniques has its assumptions and requirements and has been proven appropriate for WCET estimation as part of the verification process, a validation step is still needed during system integration phases, in accordance with safety-related systems development processes. In this context, performing a reliable timing validation process building on measurements with arbitrary overlap across tasks is a difficult challenge. Tasks overlapping heavily impacts individual tasks’ execution time and is an aspect of multicorebased systems for which no satisfactory solution exists yet. Given a reference analysis task, the number of other tasks (and the particular tasks) that overlap with it changes its execution time. Likewise the degree of overlap, i.e. how long tasks overlap, affects also its timing behavior. Task overlap varies in non-obvious manners across individual measurements and also during operation. Therefore, measurements are subject to hard-to-control contention, which brings uncertainty to the validation process. While the user might have means to observe the variability among tasks in observed TOS during tests, it is hard to enforce specific TOS. In this chapter we propose CleanET, that stands for clean execution times, an approach to distill both, contention-free measurements as well as measurements subject to specific contention levels, to enable a reliable timing validation process that captures any TOS that might occur during operation. CleanET builds upon statistical dependence analysis to derive the dilation factor 𝑟affecting execution times due to simultaneous task execution. CleanET resamples the execution times measured during testing to derive tasks execution time under (1) no overlapping scenarios (i.e. in isolation), (2) full single overlapping (i.e. when the task overlaps 100% of the time with exactly one other task), and (3) data with any specific overlap representing the worst overlap of interest (e.g. full overlap with 𝑛other tasks). CleanET estimates provide additional means of validation for derived time budgets and alleviates the pressure on end users to produce tests cases with different TOS. 89 Chapter 8. Clean Execution Times for MBPTA Proof. Start by laying out the expression for cov(𝑋, 𝑋𝑃): 2cov(𝑋, 𝑃𝑋)=var(𝑋)+var(𝑃𝑋)−var(𝑋−𝑃𝑋).(8.23) Determine var(𝑃𝑋)and var(𝑋−𝑃𝑋). var(𝑃𝑋)=var(𝑃)var(𝑋)+var(𝑃)𝐸(𝑋)2+var(𝑋)𝐸(𝑃)2,(8.24) var(𝑋−𝑃𝑋)=var(1−𝑃)var(𝑋)+var(1−𝑃)𝐸(𝑋)2+var(𝑋)𝐸(1−𝑃)2.(8.25) By subtracting those two equations we obtain the right-hand side of Equation 8.23: var(𝑋)𝐸(𝑃)2−𝐸(1−𝑃)2+1=2cov(𝑋, 𝑃𝑋),(8.26) var(𝑋)𝐸(𝑃)2−(1−𝐸(𝑃))2=2cov(𝑋, 𝑃𝑋),(8.27) var(𝑋)𝐸(𝑃)=cov(𝑋, 𝑃𝑋)>0.(8.28) Which holds in our case because 𝑃∈ [0,1]and 𝑋∈R+.□ Given that 𝑋and 𝑃are not independent, the relationship between 𝑟and 𝛽yields: 𝛽=𝑟−1 𝑟+cov(𝑋, 𝑉) var(𝑉)>𝑟−1 𝑟,(8.29) since cov(𝑋, 𝑉)>0. Note that the covariance between 𝑋and 𝑉is generally unknown. Thus, by ignoring it, we overestimate the factor 𝑟−1 𝑟, which, given that 𝑟>1as indicated before, implies that we overestimate 𝑟. The consequence of overestimating 𝑟is that, by considering scenarios where the overlapping is higher than measured, the overestimated 𝑟leads to higher predicted execution times than those in the real system. In our case, this implies that, if those predicted measurements respect timing budgets, then real system measurements would necessarily also respect the budgets. Also, if the system can overrun the timing budget, predicted measurements obtained with CleanET will indicate even higher overruns, so faults will be reliably detected. However, in some cases, time budgets may be respected and CleanET report that they are violated. Hence, while this may cost some tightness by imposing the allocation of larger time budgets than needed, our approach does not challenge the safety of the system tested. 8.3.6 CleanET for Multiple Overlappings In the previous section we focused the case where one specific job either does not overlap or partially overlap with another job. The model can be easily extended to account for several jobs overlapping at the same time (from 1to 𝑚). This is done by defining Equation 8.2 for each number of overlaps, so having 𝑊1 𝑡,··· , 𝑊𝑚 𝑡, where the property in each case is whether there are exactly 𝑖overlaps. For instance, for 𝑖=2,𝑊2 𝑡=1if and only if the job overlaps with exactly 2 other jobs in time 𝑡. For any other number of overlapping jobs (e.g. 0,1,3, etc), 𝑊2 𝑡=0. Hence, 𝑉𝑖stands for the total time where property 𝑊𝑖 𝑡holds. We can, therefore, extend our original definition of 𝑌(𝑌=𝑈+𝑉=(1−𝑞)𝑌+𝑞𝑌) decomposing 𝑉 and 𝑞across the different types of overlapping (from 1to 𝑚overlapping jobs) as follows: 𝑌=𝑈+ 𝑚 Õ 𝑖=1 𝑉𝑖= 1− 𝑚 Õ 𝑖=1 𝑞𝑖!𝑌+ 𝑚 Õ 𝑖=1 𝑞𝑖𝑌. (8.30) This allows us to reformulate Equation 8.7 as follows: 𝑌= 1− 𝑚 Õ 𝑖=1 𝑝𝑖!+ 𝑚 Õ 𝑖=1 𝑟𝑖𝑝𝑖!𝑋, (8.31) 96 8.4. Evaluation whichcan betransformed intothe followingequation where linear regression canbe directly applied: 𝑌= 𝑚 Õ 𝑖=1𝑟𝑖−1 𝑟𝑖 𝑉𝑖+𝑋, (8.32) where𝑌and all𝑉𝑖are known, and 𝑋and all 𝑟𝑖dilation factors are obtained through linear regression. 8.4 Evaluation For the evaluation of CleanET, we first estimate 𝑟𝑖dilation factors for the different task overlaps. Then, we validate how dilation factor estimation matches expectations. Finally, we show how those dilation factors can be used to estimate execution time bounds for different overlapping scenarios, thus allowing to validate whether time budgets are violated. 8.4.1 Dilation Factors (𝑟𝑖) Single Dilation Factor We have applied CleanET to obtain the dilation factors, 𝑟𝑖, for the different numbers of overlapped tasks. Results are shown in Table 8.1 considering a single dilation factor (whether overlap exists or not), as per Section 8.3.4. As shown, execution time with no overlap at all is around 16.32ms, with a standard deviation of 1.64ms. However, the impact of overlap is huge since 𝑟=6.85 with a standard deviation of 0.78. Hence, on average, we could expect the execution time will full overlap to be around 112ms (16.32ms ·6.85). For the sake of completeness, the table also provides the number of execution time observations used (289), which are those for class 0 before, and the degree of correlation measured between non-overlapped and overlapped execution time measurements, which is, as expected, very high (0.9). Table 8.1: CleanET applied to filtered data with Overlap as co-variate. Overlap (𝑟)6.85 ±0.78 Constant (𝑝=0)16.32 ±1.64 Observations 289 Adjusted R-squared (correlation): 0.90 All Dilation Factors As a second step, we have applied CleanET to obtain individual dilation factors for 1, 2 and 3 jobs overlapping with the jobs of the task under analysis, as per Section 8.3.6. Results are shown below in Table 8.2. We observe that results for 1 or 2 jobs overlapping are very similar and ranges (e.g. 𝜇±𝜎) overlap almost completely. Statistically, we cannot prove that they are distinguishable and, therefore, we conclude that overlapping with 1 or 2 jobs must be considered together, thus defining that the property in Equation 8.2 holds if the task overlaps with exactly 1 or 2 jobs. In the case of 3 jobs overlapping, we observe that (1) the value of 𝑟3is far lower than that for 𝑟1and 𝑟2, which would mean that overlapping with 3 jobs leads to a lower dilation factor (so a lower execution time increase) than overlapping with 1 or 2 jobs, which is against intuition. However, the real problem with 3 overlaps relates to the fact that, out of the 289 class 0 measurements, only 9 have 3 jobs overlapping with the task under analysis for some time, which is, in practice, too scarce data to raise any reliable prediction. In fact, the (very high) standard deviation for this case already indicates this behavior indirectly. The number of measurements including data for each overlap case is provided in Table 8.3 for completeness. Finally, note that, while the case with no overlap (“Constant” in the tables) has no meaningful change with respect to the case where we fit only 𝑟, it is not absolutely identical. This relates to the fact that all parameters in Equation 8.32 are fit together, thus dealing to minor variations when varying the parameters to fit. 97 Chapter 8. Clean Execution Times for MBPTA Table 8.2: Linear model applied on filtered data with Overlap time of one, two and three processes at a time. Overlap 1 (𝑟1)6.97 ±0.81 Overlap 2 (𝑟2)6.55 ±0.81 Overlap 3 (𝑟3)3.23 ±2.42 Constant (𝑝=0)16.32 ±1.65 Observations 289 Adjusted R-squared (correlation): 0.90 Table 8.3: Number of class 0 measurements with part of its execution with 0, 1, 2, 3 jobs overlapping. Number of Measurements Percentage overlapping jobs (w.r.t. 289) Overlap 0 137 47.4% Overlap 1 288 99.7% Overlap 2 218 75.4% Overlap 3 9 3.1% Statistically Significant Dilation Factors After concluding that 1 and 2 overlapping measurements are not statistically distinguishable, we apply CleanET again considering only two dilation factors: 𝑟1−2for 1 or 2 overlaps, and 𝑟3for 3 overlaps, as shown in Table 8.4 below. We note that 𝑟1−2is not distinguishable in practice with 𝑟in Table 8.1 when considering just one dilation factor since the amount of data for 3 overlaps is too little to cause meaningful differences. We also note that, by considering 1 and 2 overlaps together instead of separated, the case for 3 overlaps varies drastically, which reflects the weakness of the fit for 3 overlaps due to the too little data available for this case. Table 8.4: Linear model applied on filtered data with Overlap time of one and two processes added, and three processes. Overlap 1-2 (𝑟1−2)6.88 ±0.78 Overlap 3 (𝑟3)2.75 ±1.68 Constant (𝑝=0)16.31 ±1.64 Observations 289 Adjusted R-squared (correlation): 0.90 Overall, our results show that scenarios with 0, 1 or 2 overlaps can be reliably resampled from the data available. If the case with 3 overlaps is regarded as relevant for operation conditions, then additional input data would be required including many more measurements corresponding to that scenario, so that CleanET could deliver a reliable dilation factor for this case. 8.4.2 Validation of the Method In general, these solutions cannot be validated due to the lack of reference data, which in this case would correspond to measurements with the complete execution time with a constant number of jobs overlapping (e.g. 100% of the execution time overlapping with exactly 1 job). However, in our case, we have measurements with exactly 1 job overlapping during their complete execution. In fact, since 1 and 2 overlappings have been shown to be indistinguishable, and 3 overlappings occur seldom (too occasionally to be statistically significant), we consider as reference measurements those in which some overlapping (with either 1, 2 or 3 jobs) occurs during the whole execution. Hence, we performed the following: 98 8.4. Evaluation ●●●● ●●●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ●● ● ● ●●●● ● ● Figure 8.4: QQ-plot of full overlap resampled data with CleanET with respect to empirical data with full overlap. ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ●●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ●● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ECCDF Figure 8.5: Execution time bounds for no overlapping (𝑝=0) and full overlapping (𝑝=1). 1. Split the set of measurements into 2 categories: one with those measurements with 100% overlap (OVL1 group, 152 values), and those with non-full overlap (OVLmix, the remaining 137 values). 2. Apply CleanET on OVLmix data. 3. Compare the results obtained with CleanET resampling data in OVLmix to the case with full overlap, with the actual data with those overlapping characteristics (OVL1 data) with a QQ-plot (see Figure 8.4). In the figure, we would like to have a linear relation between the empirical data and the data resampled with CleanET (straight line). However, the fact that 𝑟is overestimated leads to non-linearity. Results show that CleanET, using OVLmix data, produces higher values than those observed empirically (OVL1), thus corroborating our expectation of having overestimated execution times for high overlaps due to the overestimation of 𝑟. This supports empirically our expectations with a reference data set different to that used by CleanET by partitioning the data into OVL1 and OVLmix groups. 99 Chapter 8. Clean Execution Times for MBPTA ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●●●●●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ●●●● ●●● ●●●● ●●●●● ●●● ●● ●●●●●●● ●● ●●●●● ●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●● ● ●●●●●●● ●●●●●●●●● ●●●●●●● ●●●●●● ●● ●●●●●●●●●●●●● ● ●●●●●●●● ● ●●●● ●● ●●●● ● ● ● ● ● ● Figure 8.6: Impact of the dilation factor (𝑟) on the execution time bounds as we vary the degree of overlapping (𝑝). 8.4.3 Using CleanET Results We have used CleanET to generate, from the actual measurements obtained from the system with arbitrary overlaps, measurements corresponding to overlapping scenarios of interest. In particular, since the purpose of timing validation is assessing that no overrun occurs, we have considered the case with full overlap. For the sake of illustration, we also show results for the case of no-overlap. Figure 8.5 plots the empirical complementary cumulative density function (ECCDF) for the three scenarios: actual data measured (empirical), data obtained with CleanET with no overlap (basal), and data obtained with CleanET with full overlap. Together with the empirical data, we fit an exponential tail to the highest values of the extreme cases (basal and full overlap), which has been recommended as a suitable approach to predict high execution times [3, 44]. For that purpose, we build on peaks-over-threshold methods selecting the threshold as indicated in [90]. Measurements resampled with CleanET allow considering cases that could not be explicitly tested in the system under test. As shown, as expected, the full overlap case leads to execution times higher than those of the empirical (measured) case, whereas the basal case, instead, leads to the lowest execution times one would expect in the real system. Note that, as discussed in Section 8.3.5, 𝑟values are overestimated. Hence, the full overlap case is an overestimation of what would be in practice the behavior of the real system with full overlap, whose timing behavior would be somewhere in the region between the empirical measurements for arbitrary overlaps and the estimated full-overlap measurements obtained with CleanET. Note also that overestimating 𝑟makes that, whenever we consider lower overlaps than in practice (e.g. basal case), execution times obtained are in practice lowerbounds of the true basal case. In any case, the basal case is irrelevant to validate whether time budgets allocated suffice. For the sake of completeness, we provide in Figure 8.6 how execution times expected vary as we vary the degree of overlap (𝑝). We use peak-over-threshold again, and we collect the values obtained at different exceedance levels (lvl), from 10−3to 10−12. As expected, a linear increase of the degree of overlap, leads to nearly-linear increase of the execution times expected, with small disturbances due to the uncertainty associated to threshold selection in the fitting process of a distribution (an exponential distribution in our case). For instance, based on the results obtained, if we consider an exceedance level of 10−6, we would conclude that execution times could be slightly below 500ms (assuming full overlap being the worst-case during operation). Thus, validation tests would be passed if such value is lower than the deadline for the task under analysis. Overall, CleanET allows exploring any arbitrary degree of overlap of interest, thus providing end users with means to validate their system against scenarios that cannot be triggered during system testing, thus relieving end users from having to create additional test cases with the hope of those test cases to produce the execution scenarios that need being considered. 100 8.5. Summary 8.5 Summary The complexity of current MPSoCs is mainly due to the intricate interactions between their multiple hardware blocks in the execution of complex software like the one governing Autonomous Driving. The competition for resources in these hardware blocks impacts greatly the execution time of each task. In this work we propose cleanET, a mathematical model that expresses how tasks that overlap in time produce a time dilation on each other. CleanET not only can model the dilation factor of a given set of time data, it can also use that factor to estimate a worst-case scenario where the contention is at its peak. We used as case-study data from a real autonomous driving framework, Apollo, to assess the dilation factor of overlapping tasks. CleanET allows generating, apart from the worst-case full overlap scenario, a no overlap scenario so that we can observe a model for the execution times in isolation with only the data during system operation. 101 Chapter 8. Clean Execution Times for MBPTA 102 Part V Conclusions and Future Work 103 Chapter 9 Conclusions and Future Work 9.1 Conclusions The growth in the complexity of CRTES continues to challenge the real-time community to provide trustworthy timing analysis solutions. The execution of increasingly-complex mixed-critically software, like machine learning, on complex multicore hardware makes timing analysis, and providing evidence of correctness, difficult tasks. In order to address the complexity rise in CRTES, several measurement-based techniques have been proposed. In particular, measurement probabilistic timing analysis techniques based on Extreme Value Theory have been shown to deliver provably trustworthy WCET estimates on execution scenarios that have not been observed. In order for MBPTA techniques to be applied, the execution time occurrences under validation should be representative of all possible execution conditions and scenarios that may arise during system operation, hardware- and software-wise. The particularity of MBPTA resides in providing WCET bounds via measurements in the form of a probability distribution. This thesis pushes the state-of-the-art in the MBPTA domain: We have proposed new tools to i) provide tighter and trustworthy pWCET estimates, ii) merge disjoint hardware event monitor data, and iii) assess contention from competing tasks. More specifically, we achieved the following: •Sky-high Quantile Estimation: Light tails are the theoretical best fit for pWCET estimation. However, the lack of reliable fitting methods (for light tails as part of EVT) that deliver reliable pWCET estimates for arbitrarily low exceedance probabilities imposes the use of exponential tails, which although proven reliable, are increasingly pessimistic for decreasing exceedance probabilities. We proposed two different solutions for this problem. –We proved that risk analysis, where EVT is used, and survivability analysis tackle the same fundamental problem by predicting the occurrence of rare events. Then, building on distributions from survivability analysis, we proposed the use of Weibull tails (tailW), which are proven to be as reliable as reference log-concave distributions, but enabling the tight modelling of arbitrarily low exceedance probabilities. –We presented for the first time a method based on Markov’s Inequality for pWCET estimation that represents a solid alternative approach to EVT. In particular, we showed that MIK (Markov Inequalityto the power-of-k) has no model uncertaintyand proposed amethod to handle sampling uncertainty(RESTK)that consistently providesmore trustworthy, tighter and stable results than EVT in different scenarios, including a railway case study. These promising results suggest that RESTK can be effectively used as a standalone method for pWCET estimation, or even as an alternative approach to validate EVT results in those cases where EVT is already consolidated. In this line, the fact that MIK (RESTK) and EVT build on completely different mathematical foundations provides stronger evidence on the trustworthiness of the obtained pWCET estimates. •HEM Data Merging: Measurement-based timing analysis methods increasingly build on HEMs to measure and estimate the timing behavior of time-critical applications running on 105 Bibliography [54] Deloitte. Semiconductors – the Next Wave Opportunities and winning strategies for semiconductor companies, 2019. [55] E. Díaz, E. Mezzetti, L. Kosmidis, J. Abella, and F. J. Cazorla. Modelling multicore contention on the aurixtm tc27x. In Design Automation Conference (DAC), 2018. [56] J. L. Diaz, D. F. Garcia, K. Kim, C. Lee, L. Lo Bello, J. M. Lopez, S. Min, and O. Mirabella. Stochastic analysis of periodic real-time systems. In Real-Time Systems Symposium (RTSS), pages 289–300, 2002. [57] B. Dreyer, C. Hochberger, A.Lange, S.Wegener, and A. Weiss. Continuousnon-intrusive hybrid wcet estimation using waypoint graphs. In International Workshop on Worst-Case Execution Time Analysis (WCET), 2016. [58] L. Duembgen, A. Huesler, and K. Rufibach. Active set and EM algorithms for log-concave densities based on complete and censored data. Technical report, IMSV, Univ. of Bern, 2010. [59] S. Edgar and A. Burns. Statistical analysis of wcet for scheduling. In Real-Time Systems Symposium (RTSS), pages 215–224, 2001. [60] B. Efron. Bootstrap methods: Another look at the jackknife. The Annals of Statistics, 7(1):1–26, 1979. [61] R. Ernst. Codesign of embedded systems: status and trends. Design Test of Computers, 15(2):45– 54, 1998. [62] Federal Aviation Administration, Certification Authorities Software Team (CAST). CAST-32A Multi-core Processors, 2016. [63] R. A. Fisher. Xv.—the correlation between relatives on the supposition of mendelian inheritance. Transactions of the Royal Society of Edinburgh, 52(2):399–433, 1919. [64] R. A. Fisher. Moments and product moments of sampling distributions. Proceedings of the London Mathematical Society, s2-30(1):199–238, 1930. [65] R. A. Fisher and L. H. C. Tippett. Limiting forms of the frequency distribution of the largest or smallest member of a sample. Mathematical Proceedings of the Cambridge Philosophical Society, 24(2):180–190, 1928. [66] D. Freedman, R. Pisani, and R. Purves. Statistics: Fourth International Student Edition. W.W. Norton & Company, 2007. [67] Freescalesemicondutor. ARM®architecturereferencemanual.armv8, forarmv8-aarchitecture profile. ARM DDI 0487F.b (ID040120). [68] Freescale semicondutor. e6500 Core Reference Manual. https://www.nxp.com/docs/en/ reference-manual/E6500RM.pdf, 2014. E6500RM. Rev 0. 06/2014. [69] Freescalesemicondutor. QorIQT2080ReferenceManual, 2016. AlsosupportsT2081.Document Number: T2080RM. Rev. 3, 11/2016. [70] M. Fusi, F. Mazzocchetti, A. Farres, L. Kosmidis, R. Canal, F. J. Cazorla, and J. Abella. On the use of probabilistic worst-case execution time estimation for parallel applications in high performance systems. Mathematics, 8(3), 2020. [71] S. Jiménez Gil, I. Bate, G. Lima, L. Santinelli, A. Gogonel, and L. Cucu-Grosjean. Open challenges for probabilistic measurement-based worst-case execution time. Embedded Systems Letters, 9(3):69–72, 2017. [72] S. Girbal, M. Moretó, A. Grasset, J. Abella, E. Quiñones, F. J. Cazorla, and S. Yehia. On the convergence of mainstream and mission-critical markets. In Design Automation Conference (DAC), New York, NY, USA, 2013. ACM. [73] T. Grass, A. Rico, M. Casas, M. Moreto, and A. Ramirez. Evaluating execution time predictability of task-based programs on multi-core processors. In Parallel Processing Workshops, pages 218–229. Springer International Publishing, 2014. 112 Bibliography [74] D. Griffin and A. Burns. Realism in Statistical Analysis of Worst Case Execution Times. In Björn Lisper, editor, International Workshop on Worst-Case Execution Time Analysis (WCET), volume 15, pages 44–53, Dagstuhl, Germany, 2010. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik. The printed version of the WCET’10 proceedings are published by OCG (www.ocg.at) - ISBN 978-3-85403-268-7. [75] D. Griffin, B. Lesage, I. Bate, F. Soboczenski, and R. I. Davis. Forecast-based interference: Modelling multicore interference from observable factors. In International Conference on Real- Time Networks and Systems (RTNS), page 198–207, New York, NY, USA, 2017. ACM. [76] F. Guet, L. Santinelli, and J. Morio. On the Reliability of the Probabilistic Worst-Case Execution Time Estimates. In Embedded Real-time Software and Systems (ERTS2) Conference, 2016. [77] F.Guet, L.Santinelli,andJ.Morio. Probabilisticanalysisofcachememoriesandcachememories impacts on multi-core embedded systems. In Symposium on Industrial Embedded Systems (SIES), pages 1–10, 2016. [78] F. Guet, L. Santinelli, and J. Morio. On the representativity of execution time measurements: Studying dependence and multi-mode tasks. In International Workshop on Worst-Case Execution Time Analysis (WCET), pages 3:1–3:13, 2017. [79] Y. Guo, T. Sayed, L. Zheng, and M. Essa. An extreme value theory based approach for calibration of microsimulation models for safety analysis. Simulation Modelling Practice and Theory, 106:102172, 2021. [80] M. Gyllenhammar, R. Johansson, F. Warg, D. Chen, H. Heyn, M. Sanfridson, J. Söderberg, A. Thorsén, and S. Ursing. Towards an operational design domain that supports the safety argumentation of an automated driving system. In European Congress on Embedded Real Time Systems (ERTS), 2020. [81] P. R. Halmos. The Theory of Unbiased Estimation. The Annals of Mathematical Statistics, 17(1):34 – 43, 1946. [82] J. Hansen, S. Hissam, and G. A. Moreno. Statistical-Based WCET Estimation and Validation. In Niklas Holsti, editor, International Workshop on Worst-Case Execution Time Analysis (WCET), volume 10, pages 1–11, Dagstuhl, Germany, 2009. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik. also published in print by Austrian Computer Society (OCG) with ISBN 978-3- 85403-252-6. [83] T. Hastie, R. Mazumder, J. D. Lee, and R. Zadeh. Matrix completion and low-rank svd via fast alternating least squares. J. Mach. Learn. Res., 16(1):3367–3402, 2015. [84] M. L. Hazelton. Assessing log-concavity of multivariate densities. Statistics and Probability Letters, 81(1):121–125, 2011. [85] C. Hernández, J. Abella, F. J. Cazorla, A. Bardizbanyan, J. Andersson, F. Cros, and F. Wartel. Design and Implementation of a Time Predictable Processor: Evaluation With a Space Case Study. In Marko Bertogna, editor, Euromicro Conference on Real-Time Systems (ECRTS), volume 76, pages 16:1–16:23, Dagstuhl, Germany, 2017. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik. [86] B. M. Hill. A Simple General Approach to Inference About the Tail of a Distribution. The Annals of Statistics, 3(5):1163 – 1174, 1975. [87] R. V. Hogg and E. A. Tanis. Probability and statistical inference. Prentice Hall, 1997. [88] R. Hyndman and Y. Fan. Sample quantiles in statistical packages. The American Statistician, 50(4):361–365, 1996. [89] InternationalOrganizationforStandardization. ISO/DIS 26262. Road Vehicles – Functional Safety, 2009. [90] J. Daoudi J. del Castillo and R. Lockhart. Methods to distinguish between polynomial and exponential tails. Scandinavian Journal of Statistics, 41(2):382–393, 2014. 113 Bibliography [91] A. F. Jenkinson. The frequency distribution of the annual maximum (or minimum) values of meteorological elements. Quarterly Journal of the Royal Meteorological Society, 81(348):158–171, 1955. [92] N. L. Johnson. Continuous univariate distributions. Vol. 2. Wiley and Sons, New York, 2nd ed. / norman l. johnson, samuel kotz, n. balakrishnan. edition, 1994. [93] H. Kim, D. de Niz, B. Andersson, M. Klein, O. Mutlu, and R. Rajkumar. Bounding memory interference delay in cots-based multi-core systems. In Real-Time and Embedded Technology and Applications Symposium (RTAS), pages 145–154. [94] T. Kloda, M. Solieri, R. Mancuso, N. Capodieci, P. Valente, and M. Bertogna. Deterministic memory hierarchy and virtualization for modern multi-core embedded systems. In Real-Time and Embedded Technology and Applications Symposium (RTAS), pages 1–14, 2019. [95] M. Kolmogorov. Sulla determinazione empírica di uma legge di distribuzione. Giornale dell’Istituto Italiano degli Attuari, 1933. [96] H. Kopetz. The complexity challenge in embedded system design. In International Symposium on Object and Component-Oriented Real-Time Distributed Computing (ISORC), pages 3–12, 2008. [97] L. Kosmidis, J. Abella, E. Quiñones, and F. J. Cazorla. A cache design for probabilistically analysable real-time systems. In Enrico Macii, editor, Design, Automation and Test in Europe (DATE), pages 513–518. EDA Consortium San Jose, CA, USA / ACM DL, 2013. [98] L. Kosmidis, E. Quiñones, J. Abella, G. Farrall, F. Wartel, and F. J. Cazorla. Containing timingrelated certification cost in automotive systems deploying complex hardware. In Design Automation Conference (DAC), page 1–6, New York, NY, USA, 2014. ACM. [99] L. Kosmidis, E. Quiñones, J. Abella, T. Vardanega, C. Hernandez, A. Gianarro, I. Broster, and F. J. Cazorla. Fitting processor architectures for measurement-based probabilistic timing analysis. Microprocessors and Microsystems, 47:287–302, 2016. [100] S. Kotz and S. Nadarajah. Extreme value distributions: theory and applications. World Scientific, 2000. [101] B. Light. Concentration inequalities using higher moments information. arXiv, 2020. [102] H. W. Lilliefors. On the kolmogorov-smirnov test for normality with mean and variance unknown. Journal of the American Statistical Association, 62(318):399–402, 1967. [103] R. V. Lim. Computationally efficient multiplexing of events on hardware counters. In Linux Symposium, 2014. [104] G. Lima and I. Bate. Valid Application of EVT in Timing Analysis by Randomising Execution Time Measurements. In Real-Time and Embedded Technology and Applications Symposium (RTAS), pages 187–198, 2017. [105] G. Lima, D. Dias, and E. Barros. Extreme value theory for estimating task execution time bounds: A careful look. In Euromicro Conference on Real-Time Systems (ECRTS), pages 200–211, 2016. [106] J. S. Liu and S. M. Shahrier. On predictability of caches for real-time applications. In International Workshop on Modeling, Analysis and Simulation of Computer and Telecommunication Systems (MASCOT), pages 52–56, 1994. [107] L. Liu, Z. Cui, M. Xing, Y. Bao, M. Chen, and C. Wu. A software memory partition approach for eliminating bank-level interference in multicore systems. In Pact, pages 367–376. ACM, 2012. [108] G. M. Ljung and G. E. P. Box. On a measure of lack of fit in time series models. Biometrika, 65(2):297–303, 1978. [109] Y. Lu, T. Nolte, I. Bate, and L. Cucu-Grosjean. A new way about using statistical analysis of worst-case execution times. SIGBED Rev., 8(3):11–14, 2011. 114 Bibliography [110] R. Mancuso, R. Dudko, E. Betti, M. Cesati, M. Caccamo, and R. Pellizzoni. Real-time cache management framework for multi-core architectures. In Real-Time and Embedded Technology and Applications Symposium (RTAS), pages 45–54. IEEE, 2013. [111] A. Markov. On certain applications of algebraic continued fractions. Ph.D. thesis, St. Petersburg, 1884. [112] A.W. Marshall and I. Olkin. Life distributions. Structure of Nonparametric, Semiparametric, and Parametric Families. Springer, 2007. [113] Maspatechnologies. Maspatechnologies Helps Airbus Pass Multicore Certification on NXP T2080. https://maspatechnologies.com/wp-content/uploads/2022/08/A3R_ MulticoreCertification_ADSMPT.pdf, 2022. [114] A. Melani, E. Noulard, and L. Santinelli. Learning from probabilities: Dependences within real-time systems. In Conference on Emerging Technologies Factory Automation (ETFA), pages 1–8, 2013. [115] E. Mezzetti, L. Kosmidis, J. Abella, and F. J. Cazorla. High-integrity performance monitoring units in automotive chips for reliable timing v v. Micro, 38(1):56–65, 2018. [116] T. Mikosch. Regular variation, subexponentiality and their applications in probability theory. International Journal of Production Economics, 1999. [117] S. Milutinovic. On the limits of probabilistic timing analysis. PhD thesis, UPC, Departament d’Arquitectura de Computadors, 2019. [118] S. Milutinovic, E. Mezzetti, J. Abella, T. Vardanega, and F. J. Cazorla. On uses of Extreme Value Theory fit for industrial-quality WCET analysis. In International Symposium on Industrial Embedded Systems (SIES), pages 1–6. IEEE, 2017. [119] T. Mytkowicz, A. Diwan, M. Hauswirth, and P. Sweeney. Producing wrong data without doing anything obviously wrong! volume 44, pages 265–276, 2009. [120] T. Naoi, Y. Kagawa, K. Nagino, S. Niwa, and K. Hayashi. Extreme value analysis of the velocity of axonal transport by kinesin and dynein. Biophysical Journal, 121(3, Supplement 1):400a, 2022. [121] R. Neill, A.Drebes, and A. Pop. Fuse: Accurate multiplexingof hardware performance counters across executions. Transactions on Architecture and Code Optimization, 14(4), 2017. [122] R. B. Nelsen. An Introduction to Copulas. Springer Publishing Company, Incorporated, 2010. [123] J. Neyman, E. Pearson, and K. Pearson. Ix. on the problem of the most efficient tests of statistical hypotheses. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 231(694-706):289–337, 1933. [124] J. Neyman and E. S. Pearson. On the problem of the most efficient tests of statistical hypotheses. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 231:289–337, 1933. [125] J. Nowotsch, M. Paulitsch, D. Bühler, H. Theiling, S. Wegener, and M. Schmidt. Multi-core interference-sensitive wcet analysis leveraging runtime resource capacity enforcement. In Euromicro Conference on Real-Time Systems (ECRTS), pages 109–118, 2014. [126] J. Nowotsch, M. Paulitsch, D. Bühler, H. Theiling, S. Wegener, and M. Schmidt. Multi-core interference-sensitive WCET analysis leveraging runtime resource capacity enforcement. In Euromicro Conference on Real-Time Systems (ECRTS), pages 109–118, 2014. [127] S. Osborne and J. H. Anderson. Simultaneous multithreading and hard real time: Can it be safe? In Marcus Völp, editor, Euromicro Conference on Real-Time Systems (ECRTS), volume 165, pages 14:1–14:25. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2020. [128] X. Pan and F. Mueller. Controller-aware memory coloring for multicore real-time systems. In Symposium on Applied Computing (SAC), pages 584–592. ACM, 2018. 115 Bibliography [129] R. Pellizzoni, E. Betti, S. Bak, G. Yao, J. Criswell, M. Caccamo, and R. Kegley. A predictable execution model for COTS-based embedded systems. In Real-Time and Embedded Technology and Applications Symposium (RTAS), pages 269–279, 2011. [130] R. Pellizzoni, B. D. Bui, M. Caccamo, and L. Sha. Coscheduling of CPU and I/O Transactions in COTS-Based Embedded Systems. In Real-Time Systems Symposium (RTSS), pages 221–231, 2008. [131] R. Pellizzoni, A. Schranzhofer, J. Chen, M. Caccamo, and L. Thiele. Worst case delay analysis for memory interference in multicore systems. In Design, Automation Test in Europe (DATE), pages 741–746, 2010. [132] J. Pickands. Statistical inference using extreme order statistics. The Annals of Statistics, 3(1):119– 131, 1975. [133] R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2021. [134] D. Radack, H. G. Tiedeman Jr, and P. J. Parkinson. Civil Certification of Multi-core Processing Systems in Commercial Avionics, 2018. [135] A. Ravindar and Y. N. Srikant. Estimation of probabilistic bounds on phase CPI and relevance in WCET analysis. In International Conference on Embedded Software, page 165–174, New York, NY, USA, 2012. ACM. [136] B. Recht. A simpler approach to matrix completion. J. Mach. Learn. Res., 12(null):3413–3430, 2011. [137] F. Reghenzani, G. Massari, and W. Fornaciari. Probabilistic-WCET reliability: Statistical testing of EVT hypotheses. Microprocess. Microsystems, 77:103–135, 2020. [138] F. Reghenzani, L. Santinelli, and W. Fornaciari. Dealing with uncertainty in pWCET estimations. Transactions on Embedded Computing Systems, 19(5):33:1–33:23, 2020. [139] RTCA and EUROCAE. DO-178C / ED-12C, Software Considerations in Airborne Systems and Equipment Certification, 2011. [140] C. El Salloum, M. Elshuber, O. Höftberger, H. Isakovic, and A. Wasicek. The across mpsoc – a new generation of multi-core processors designed for safety-critical embedded systems. In Euromicro Conference on Digital System Design, pages 105–113, 2012. [141] L. Santinelli, F. Guet, and J. Morio. Revising measurement-based probabilistic timing analysis. In Real-Time and Embedded Technology and Applications Symposium (RTAS), pages 199–208, 2017. [142] L. Santinelli, J. Morio, G. Dufour, and D. Jacquemart. On the Sustainability of the Extreme Value Theory for WCET Estimation. In Heiko Falk, editor, International Workshop on Worst-Case Execution Time Analysis (WCET), volume 39, pages 21–30, Dagstuhl, Germany, 2014. Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik. [143] A. Satyanarayana. Intelligent sampling for big data using bootstrap sampling and Chebyshev inequality. In Canadian Conference on Electrical and Computer Engineering (CCECE), pages 1–6, 2014. [144] K. Schmidt, D.Marx, J. Harnisch, A.Mayer, U. Dannebaum, andH. Christlbauer. Non-intrusive tracing at first instruction. In SAE Technical Paper. SAE International, 2015. [145] A. Schranzhofer, R. Pellizzoni, J. J. Chen, L. Thiele, and M. Caccamo. Worst-case response time analysis of resource access models in multi-core systems. In Design Automation Conference (DAC), pages 332–337, 2010. [146] J. Serrà, A. Corral, M. Boguñá, M. Haro, and J. Ll. Arcos. Measuring the evolution of contemporary western popular music. Scientific reports, 2(1):1–6, 2012. [147] H. Servat, G. Llort, J. Giménez, K. Huck, and J. Labarta. Folding: detailed analysis with coarse sampling. In Tools for High Performance Computing 2011, pages 105–118. Springer, 2012. 116 Bibliography [148] D. Shapiro. Introducing xavier, the nvidia ai supercomputer for the future of autonomous transportation. NVIDIA blog, 2016. [149] K. P. Silva, L. F. Arcaro, and R. Silva De Oliveira. On Using GEV or Gumbel Models When Applying EVT for Probabilistic WCET Estimation. In Real-Time Systems Symposium (RTSS), pages 220–230, 2017. [150] N. Smirnov. On the estimation of the discrepancy between empirical curves of distribution for two independent samples. Moscow University Mathematics Bulletin, 1939. [151] D. Sornette. Critical Phenomena in Natural Sciences: Chaos, Fractals, Selforganization and Disorder: Concepts and Tools. Springer, 2006. [152] S. M. Steinberg and C. E. Davis. Distribution-free confidence intervals for quantiles in small samples. Communications in Statistics - Theory and Methods, 14(4):979–990, 1985. [153] Z. Stephenson, J. Abella, and T. Vardanega. Supporting industrial use of probabilistic timing analysis with explicit argumentation. In IEEE International Conference on Industrial Informatics (INDIN), pages 734–740, 2013. [154] Student. The probable error of a mean. Biometrika, 6(1):1–25, 1908. [155] N. Suzuki, H. Kim, D. de Niz, B. Andersson, L. Wrage, M. H. Klein, and R. Rajkumar. Coordinated bank and cache coloring for temporal protection of memory accesses. In CSE, pages 685–692. IEEE Computer Society, 2013. [156] P. Tchebichef. Des valeurs moyennes. Journal de mathématiques pures et appliquées, 12(2):177–184, 1867. [157] D. Trilla. Non-functional considerations of time-randomized processor architectures. PhD thesis, UPC, Departament d’Arquitectura de Computadors, 2020. [158] Udacity. An Open Source Self-Driving Car. https://github.com/udacity/self-driving- car/, 2017. [159] http://www.gaisler.com/index.php/products/processors/leon3.Leon3 Processor. [160] V. Utkin. Calculating the reliability of machine parts on the basis of the Chebyshev inequality. Russian Engineering Research, 32, 2012. [161] P. Kumar Valsan, H. Yun, and F. Farshchi. Taming non-blocking caches to improve isolation in multicore real-time systems. In Real-Time and Embedded Technology and Applications Symposium (RTAS), pages 1–12, 2016. [162] S. H. VanderLeest and C. Evripidou. An approach to verification of interference concerns for multicore systems (CAST-32A). In SAE Technical Paper. SAE International, 2020. [163] G. V.G. Baranoski, J. G. Rokne, and G. Xu. Applying the exponential Chebyshev inequality to the nondeterministic computation of form factors. Journal of Quantitative Spectroscopy and Radiative Transfer, 69(4):447–467, 2001. [164] S. Vilardell and À. Pineda. distTails: A Collection of Full Defined Distribution Tails, 2019. R package version 0.1.2. [165] S. Vilardell, I. Serra, J. Abella, J. del Castillo, and F. J. Cazorla. Software timing analysis for complex hardware with survivability and risk analysis. In International Conference on Computer Design (ICCD), pages 227–236, 2019. [166] S. Vilardell, I. Serra, R. Santalla, E. Mezzetti, J. Abella, and F. J. Cazorla. Hrm: merging hardware event monitors for improved timing analysis of complex mpsocs. Transactions on Computer-Aided Design of Integrated Circuits and Systems (TCAD), 2020. [167] G. von der Brüggen, N. Piatkowski, K. Chen, J. Chen, and K. Morik. Efficiently approximating the probability of deadline misses in real-time systems. In Euromicro Conference on Real-Time Systems (ECRTS), 2018. [168] B. von Mises. Fundamentalsätze der wahrscheinlichkeitsrechnung. Mathematische Zeitschrift, 4:1–97, 1919. 117 Bibliography [169] H. Wang and K. J. Jonas. The likelihood of severe covid-19 outcomes among plhiv with various comorbidities: a comparative frequentist and bayesian meta-analysis approach. Journal of the International AIDS Society, 24(11):e25841, 2021. [170] F. Wartel, L. Kosmidis, A. Gogonel, A. Baldovin, Z. R. Stephenson, B. Triquet, E. Quiñones, C. Lo, E.Mezzetti, I. Broster, J. Abella, L.Cucu-Grosjean, T. Vardanega, and F. J.Cazorla. Timing analysis of an avionics case study on complex hardware/software platforms. In Wolfgang Nebel and David Atienza, editors, Design, Automation and Test in Europe (DATE), pages 397– 402. ACM, 2015. [171] V. M.Weaver, D. Terpstra, and S.Moore. Non-determinismand overcount onmodern hardware performance counter implementations. In International Symposium on Performance Analysis of Systems and Software (ISPASS), pages 215–224, 2013. [172] R. Wilhelm, J. Engblom, A. Ermedahl, N. Holsti, S. Thesing, D. Whalley, G. Bernat, C. Ferdinand, R. Heckmann, T. Mitra, F. Mueller, I. Puaut, P. Puschner, J. Staschulat, and P. Stenström. The worst-case execution-time problem–overview of methods and survey of tools. Transactions on Embedded Computing Systems, 7(3):36:1–36:53, 2008. [173] R. Wilhelm, J. Engblom, A. Ermedahl, N. Holsti, S. Thesing, D. Whalley, G. Bernat, C. Ferdinand, R. Heckmann, T. Mitra, F. Mueller, I. Puaut, P. Puschner, J. Staschulat, and P. Stenström. Theworst-caseexecution-timeproblem—overviewofmethodsandsurveyoftools. Transactions on Embedded Computing Systems, 7(3), 2008. [174] D.J. Wilkins. The bathtub curve and product failure behavior, part one: The bathtub curve, infant mortality and burn-in. Reliability Hotwire, (21), 2002. [175] S. S. Wilks. The large-sample distribution of the likelihood ratio for testing composite hypotheses. The Annals of Mathematical Statistics, 9(1):60–62, 1938. [176] S. N. Wood. Core Statistics. Cambridge University Press, 2015. [177] J. Worms and S. Touati. Parametric and Non-Parametric Statistics for Program Performance Analysis and Comparison. [Research Report] RR-8875, INRIA Sophia Antipolis - I3S; Université NiceSophia Antipolis; Université Versailles Saint Quentin en Yvelines; Laboratoire de mathématiques deVersailles, 2017. [178] J. Worms and S. Touati. Modelling program’s performance with gaussian mixtures for parametric statistics. Transactions on Multi-Scale Computing Systems (TMSCS), 4(3):383–395, 2018. [179] Xilinx. Zynq ultrascale+ device technical reference manual. https://docs.xilinx.com/v/u/ en-US/ug1085-zynq-ultrascale-trm, 2020. [180] H. Xu, Q. Wang, S. Song, L. K. John, and X. Liu. Can we trust profiling results? understanding and fixing the inaccuracy in modern profilers. In International Conference on Supercomputing, page 284–295, New York, NY, USA, 2019. ACM. [181] X. Yan and X. G. Su. Linear Regression Analysis: Theory and Computing, chapter 2. WorldScientific Publishing Co., Inc., River Edge, NJ, USA, 2009. [182] H. Yun, R. Mancuso, Z. P. Wu, and R. Pellizzoni. PALLOC: DRAM bank-aware memory allocator for performance isolation on multicore platforms. In Real-Time and Embedded Technology and Applications Symposium (RTAS), pages 155–166. IEEE, 2014. [183] H. Yun, G. Yao, R. Pellizzoni, M. Caccamo, and L. Sha. Memguard: Memory bandwidth reservation system for efficient performance isolation in multi-core platforms. In Real-Time and Embedded Technology and Applications Symposium (RTAS), pages 55–64, 2013. [184] D. Zaparanuks, M. Jovic, and M. Hauswirth. Accuracy of performance counter measurements. In International Symposium on Performance Analysis of Systems and Software (ISPASS), pages 23–32, 2009. [185] P. G. Zaykov and J. Kubalčik. Worst-case measurement-based statistical tool. In Aerospace Conference, pages 1–10, 2019. 118 Bibliography [186] M. Ziccardi, E. Mezzetti, T. Vardanega, J. Abella, and F. J. Cazorla. Epc: Extended path coverage for measurement-based probabilistic timing analysis. In IEEE Real-Time Systems Symposium (RTSS), pages 338–349, 2015. [187] D. Ziegenbein and A. Hamann. Timing-aware control software design for automotive systems. In Design Automation Conference (DAC), pages 56:1–56:6, New York, NY, USA, 2015. ACM. 119