# Suggestions on performance

I have following piece of code that is called by a program many thousands of times as part of Monte-Carlo simulation. using gprof, I see that 38% of the time was spent on this function. Are there any obvious areas of improvement in this code?

double StudentTCopula::uniform_to_default_time_student(double u, const
std::vector<double>& times, const std::vector<double>& values)
{
if (u == 0.0)
return 99999.0;
if (u == 1.0)
return 0.0;
size_t num_points = times.size();
size_t index = 0;
for (size_t i{1};i<num_points;++i){
if (u <= values[i - 1] && u > values[i]){
index = i;
break;
}
}
double tau = 0.0;
if (index == num_points + 1) {
auto t1 = times[num_points - 1];
auto q1 = values[num_points - 1];
auto t2 = times[num_points];
auto q2 = values[num_points];
auto lam = log(q1 / q2) / (t2 - t1);
tau = t2 - log(u / q2) / lam;
} else if (index == 0){
auto t1 = times.back();
auto q1 = values.back();
auto t2 = times[index];
auto q2 = values[index];
tau = (t1 * log(q2 / u) + t2 * log(u / q1)) / log(q2 / q1);
} else {
auto t1 = times[index - 1];
auto q1 = values[index - 1];
auto t2 = times[index];
auto q2 = values[index];
tau = (t1 * log(q2 / u) + t2 * log(u / q1)) / log(q2 / q1);
}
return tau;
}


I have tried removing the intermediate variables and directly genearting a long return expression with the calculations, which did not affect gprof results.

• Please edit your title to describe what the code does, not your issues with it. See How to Ask. Jul 14 at 20:44
• Please clarify: (1) are times.size() and values.size() related in any way? (2) how large is time.size()? (3) are values sorted indeed?
– vnp
Jul 15 at 0:51
• Are you sure that times[num_points] and values[num_points] are in bounds? To me, it looks like you may be reading data past the end of the vectors. Are you compiling your debug builds with warnings and sanitizers?
– Pkkm
Jul 15 at 15:38
• Post more of the code. We don't know if your bottleneck is in this function or some where else in the program. Jul 19 at 15:34

• Comparing floating point representations for equality has its own risks - consider

    if (u <= 0 + EPS)
return ‹BIG_VALUE›;
if (1 <= u + EPS)
return 0;


For suitable values of ‹BIG_VALUE›(there should be a better name) and EPS (possibly 0).

• Most of the handling for index != num_points + 1 is common. Don't repeat yourself:

  } else {
auto t1 = (index == 0) ? times.back() : times[index - 1];
auto q1 = (index == 0) ? values.back() : values[index - 1];
auto t2 = times[index];
auto q2 = values[index];
tau = (t1 * log(q2 / u) + t2 * log(u / q1)) / log(q2 / q1);
}


Just before the end of uniform_to_default_time_student(), code indentation is slightly inconsistent.

G. Sliepen: use log(a/b)=log(a)−log(b) to avoid the division:

auto log_u = log(u),
log_q1 = log(q1),
log_q2 = log(q2);
tau = (t1 * (log_q2 - log_u) + t2 * (log_u - log_q1)) / (log_q2 - log_q1);


(implementation imperfections by greybeard)

I have no idea whether substituting two more divisions is useful:

  if (index == num_points + 1) {
auto t1 = times[num_points - 1];
auto q1 = values[num_points - 1];
auto t2 = times[num_points];
auto _q2 = 1 / values[num_points];
auto _lam = (t2 - t1) / log(q1 * _q2);
tau = t2 - log(u * _q2) * _lam;
}

• Thank you all for your suggestions. Jul 17 at 12:16

Disclaimer: Review based on assumptions

From the code snippet, I am guessing that values are in descending order. In this case you may use std::lower_bound to do a binary search instead of a for loop. Apart from this you can do a micro optimization by replacing log(q2/q1) with log(q2/u) + log(u/q1).

• If they’re not sorted, is it justified to stop at the first pair of bounds spanning u? Jul 14 at 20:50
• I dont think so. That wiy I assumed that it is sorted. But algorithms are weird and stating your assumptions upfront has served me extremely well so far. Jul 15 at 12:22

Division is usually a relatively slow operation. You're using the same divisor for a couple of division operations in a few places, so you might gain a little by converting something like this:

tau = (t1 * log(q2 / u) + t2 * log(u / q1)) / log(q2 / q1);


...to something more like this:

double q1i = 1.0 / q1;
tau = t1 * (log(q2 / u) + t2 * log(u * q1i)) / log(q2 * q1i);


May help a little (but no guarantee--your compiler may recognize the same possibility on its own).

At least in my testing, a * log(b) + c * log(d) is slower than just using pow to compute the power, so changing:

double tau1 = (t1 * log(q2 / u) + t2 * log(u / q1)) / log(q2 / q1);


...to:

double tau2 = std::log(std::pow(q2 / u, t1) * std::pow(u / q1, t2)) / log(q2 / q1);


...seems to improve speed by a factor of about 2. I'd guess there's less likelihood of the compiler recognizing the identity involved here, but I suppose it's possible (or it could turn out that with the CPU/flags you're using, the original equation is just as fast; impossible to be sure).

• You can also use the identity $\log(a/b) = \log(a) - \log(b)$ to avoid the division. Jul 16 at 9:50
• @G.Sliepen: You can, but I doubt it'll do much (if any) good. Jul 16 at 17:14
• @G.Sliepen But then you'll have to compute the logarithm twice, which likely takes longer than a floating point division Jul 17 at 21:56
• @chrysante: for a single application, you are right. But there are three logs in the tau expression anyway, and three quotients in its original form with just three values as numerators and denominators. Taking the log of each of these allows to substitute three divisions with three subtractions. Jul 18 at 16:39

I don't understand this stanza:

  if (index == num_points + 1) {

auto t1 = times[num_points - 1];
auto q1 = values[num_points - 1];
auto t2 = times[num_points];
auto q2 = values[num_points];
auto lam = log(q1 / q2) / (t2 - t1);
tau = t2 - log(u / q2) / lam;
} else ...


After initializing at zero, we have the invariant that i, and therefore index, shall always be strictly less than the number of points. As @Pkkm observed, a [num_points] de-reference Would Be Bad, if it ran.

Surely the stanza is equivalent to:

  if (0) {
/* NOTREACHED */
} else ...


Why did your toolchain not issue "C4702 unreachable code" or a similar warning? Does it lack a -Wunreachable-code switch? (Gcc, alas, removed that switch in 2010, though it does offer -Wunreachable-code-ctrl.)

Recommend deleting it, or changing that conditional.

Please assert times.size() == values.size()`. Or convert to passing around an array of structs, each holding a time and a value.

It wasn't obvious to me how to go from this definition to the OP's calculations. Within the source code, please cite a reference.