1// Copyright 2016 Ismael Jimenez Martinez. All rights reserved.
2//
3// Licensed under the Apache License, Version 2.0 (the "License");
4// you may not use this file except in compliance with the License.
5// You may obtain a copy of the License at
6//
7//     http://www.apache.org/licenses/LICENSE-2.0
8//
9// Unless required by applicable law or agreed to in writing, software
10// distributed under the License is distributed on an "AS IS" BASIS,
11// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12// See the License for the specific language governing permissions and
13// limitations under the License.
14
15// Source project : https://github.com/ismaelJimenez/cpp.leastsq
16// Adapted to be used with google benchmark
17
18#include "benchmark/benchmark.h"
19
20#include <algorithm>
21#include <cmath>
22#include "check.h"
23#include "complexity.h"
24
25namespace benchmark {
26
27// Internal function to calculate the different scalability forms
28BigOFunc* FittingCurve(BigO complexity) {
29  static const double kLog2E = 1.44269504088896340736;
30  switch (complexity) {
31    case oN:
32      return [](int64_t n) -> double { return static_cast<double>(n); };
33    case oNSquared:
34      return [](int64_t n) -> double { return std::pow(n, 2); };
35    case oNCubed:
36      return [](int64_t n) -> double { return std::pow(n, 3); };
37    case oLogN:
38      /* Note: can't use log2 because Android's GNU STL lacks it */
39      return [](int64_t n) { return kLog2E * log(static_cast<double>(n)); };
40    case oNLogN:
41      /* Note: can't use log2 because Android's GNU STL lacks it */
42      return [](int64_t n) { return kLog2E * n * log(static_cast<double>(n)); };
43    case o1:
44    default:
45      return [](int64_t) { return 1.0; };
46  }
47}
48
49// Function to return an string for the calculated complexity
50std::string GetBigOString(BigO complexity) {
51  switch (complexity) {
52    case oN:
53      return "N";
54    case oNSquared:
55      return "N^2";
56    case oNCubed:
57      return "N^3";
58    case oLogN:
59      return "lgN";
60    case oNLogN:
61      return "NlgN";
62    case o1:
63      return "(1)";
64    default:
65      return "f(N)";
66  }
67}
68
69// Find the coefficient for the high-order term in the running time, by
70// minimizing the sum of squares of relative error, for the fitting curve
71// given by the lambda expression.
72//   - n             : Vector containing the size of the benchmark tests.
73//   - time          : Vector containing the times for the benchmark tests.
74//   - fitting_curve : lambda expression (e.g. [](int64_t n) {return n; };).
75
76// For a deeper explanation on the algorithm logic, please refer to
77// https://en.wikipedia.org/wiki/Least_squares#Least_squares,_regression_analysis_and_statistics
78
79LeastSq MinimalLeastSq(const std::vector<int64_t>& n,
80                       const std::vector<double>& time,
81                       BigOFunc* fitting_curve) {
82  double sigma_gn = 0.0;
83  double sigma_gn_squared = 0.0;
84  double sigma_time = 0.0;
85  double sigma_time_gn = 0.0;
86
87  // Calculate least square fitting parameter
88  for (size_t i = 0; i < n.size(); ++i) {
89    double gn_i = fitting_curve(n[i]);
90    sigma_gn += gn_i;
91    sigma_gn_squared += gn_i * gn_i;
92    sigma_time += time[i];
93    sigma_time_gn += time[i] * gn_i;
94  }
95
96  LeastSq result;
97  result.complexity = oLambda;
98
99  // Calculate complexity.
100  result.coef = sigma_time_gn / sigma_gn_squared;
101
102  // Calculate RMS
103  double rms = 0.0;
104  for (size_t i = 0; i < n.size(); ++i) {
105    double fit = result.coef * fitting_curve(n[i]);
106    rms += pow((time[i] - fit), 2);
107  }
108
109  // Normalized RMS by the mean of the observed values
110  double mean = sigma_time / n.size();
111  result.rms = sqrt(rms / n.size()) / mean;
112
113  return result;
114}
115
116// Find the coefficient for the high-order term in the running time, by
117// minimizing the sum of squares of relative error.
118//   - n          : Vector containing the size of the benchmark tests.
119//   - time       : Vector containing the times for the benchmark tests.
120//   - complexity : If different than oAuto, the fitting curve will stick to
121//                  this one. If it is oAuto, it will be calculated the best
122//                  fitting curve.
123LeastSq MinimalLeastSq(const std::vector<int64_t>& n,
124                       const std::vector<double>& time, const BigO complexity) {
125  CHECK_EQ(n.size(), time.size());
126  CHECK_GE(n.size(), 2);  // Do not compute fitting curve is less than two
127                          // benchmark runs are given
128  CHECK_NE(complexity, oNone);
129
130  LeastSq best_fit;
131
132  if (complexity == oAuto) {
133    std::vector<BigO> fit_curves = {oLogN, oN, oNLogN, oNSquared, oNCubed};
134
135    // Take o1 as default best fitting curve
136    best_fit = MinimalLeastSq(n, time, FittingCurve(o1));
137    best_fit.complexity = o1;
138
139    // Compute all possible fitting curves and stick to the best one
140    for (const auto& fit : fit_curves) {
141      LeastSq current_fit = MinimalLeastSq(n, time, FittingCurve(fit));
142      if (current_fit.rms < best_fit.rms) {
143        best_fit = current_fit;
144        best_fit.complexity = fit;
145      }
146    }
147  } else {
148    best_fit = MinimalLeastSq(n, time, FittingCurve(complexity));
149    best_fit.complexity = complexity;
150  }
151
152  return best_fit;
153}
154
155std::vector<BenchmarkReporter::Run> ComputeBigO(
156    const std::vector<BenchmarkReporter::Run>& reports) {
157  typedef BenchmarkReporter::Run Run;
158  std::vector<Run> results;
159
160  if (reports.size() < 2) return results;
161
162  // Accumulators.
163  std::vector<int64_t> n;
164  std::vector<double> real_time;
165  std::vector<double> cpu_time;
166
167  // Populate the accumulators.
168  for (const Run& run : reports) {
169    CHECK_GT(run.complexity_n, 0) << "Did you forget to call SetComplexityN?";
170    n.push_back(run.complexity_n);
171    real_time.push_back(run.real_accumulated_time / run.iterations);
172    cpu_time.push_back(run.cpu_accumulated_time / run.iterations);
173  }
174
175  LeastSq result_cpu;
176  LeastSq result_real;
177
178  if (reports[0].complexity == oLambda) {
179    result_cpu = MinimalLeastSq(n, cpu_time, reports[0].complexity_lambda);
180    result_real = MinimalLeastSq(n, real_time, reports[0].complexity_lambda);
181  } else {
182    result_cpu = MinimalLeastSq(n, cpu_time, reports[0].complexity);
183    result_real = MinimalLeastSq(n, real_time, result_cpu.complexity);
184  }
185
186  std::string run_name = reports[0].benchmark_name().substr(
187      0, reports[0].benchmark_name().find('/'));
188
189  // Get the data from the accumulator to BenchmarkReporter::Run's.
190  Run big_o;
191  big_o.run_name = run_name;
192  big_o.run_type = BenchmarkReporter::Run::RT_Aggregate;
193  big_o.aggregate_name = "BigO";
194  big_o.iterations = 0;
195  big_o.real_accumulated_time = result_real.coef;
196  big_o.cpu_accumulated_time = result_cpu.coef;
197  big_o.report_big_o = true;
198  big_o.complexity = result_cpu.complexity;
199
200  // All the time results are reported after being multiplied by the
201  // time unit multiplier. But since RMS is a relative quantity it
202  // should not be multiplied at all. So, here, we _divide_ it by the
203  // multiplier so that when it is multiplied later the result is the
204  // correct one.
205  double multiplier = GetTimeUnitMultiplier(reports[0].time_unit);
206
207  // Only add label to mean/stddev if it is same for all runs
208  Run rms;
209  rms.run_name = run_name;
210  big_o.report_label = reports[0].report_label;
211  rms.run_type = BenchmarkReporter::Run::RT_Aggregate;
212  rms.aggregate_name = "RMS";
213  rms.report_label = big_o.report_label;
214  rms.iterations = 0;
215  rms.real_accumulated_time = result_real.rms / multiplier;
216  rms.cpu_accumulated_time = result_cpu.rms / multiplier;
217  rms.report_rms = true;
218  rms.complexity = result_cpu.complexity;
219  // don't forget to keep the time unit, or we won't be able to
220  // recover the correct value.
221  rms.time_unit = reports[0].time_unit;
222
223  results.push_back(big_o);
224  results.push_back(rms);
225  return results;
226}
227
228}  // end namespace benchmark
229