Using brute force to find an energy optimal path for a given amount of time

Question

The code below implements a brute force method to find an energy optimal path for a given amount of time t_aim using recursion (the recursion function is optimalPath).

It seems the code works as expected. The transition energy and time is provided by transEnergy and transTime. v_tmp holds indices to the found transition states. If a better path is found v_tmp is copied to v which is the return value (v_tmp and v are vectors). Overall complexity is O(N_v^N_x) where N_x is the path length, i.e. number of segments with the current segment k and N_v is the number of possible states at each segment with the current state i (or branch-offs if you want so).

Most time is spent running the optimalPath function up to the first condition:

if ((E_tmp + dE) <= E_min && (t_tmp + dt) < t_aim){...

Here are my questions:

1. I somehow tried to use size_t everywhere, so that almost no casting would necessary. Does it actually make a difference in performance?
2. The lowest value of i in the main for loop has to be equal to zero.
   But for (i=N_i; i>=0; i--) resulted in wrong behaviour (probably of the unsigned type of i) so I change it to for (i=N_i; (i+1)>0; i--) - is there any better way? Is such an extra addition negligible in terms of performance?
3. Making double dE, dt; local variables made the code faster. However I tried to pass some of the other variables to optimalPath to aviod accessing them as globals, but apparently it didn't give any speed up. Any general advice here?
4. Indexing with k + N_x * (v_tmp[k] + i * N_v) is done two times. Is it better to first compute the index, caching it and the use it?
5. Does this code have any bottleneck or is time consumption evenly distributed?
6. Would it be better to avoid recursion and solve the problem using a different approach or is recursion not that bad, depending on the conditions? In the given case the maximum recursion depth is the length of the path N_x (I think).
7. What should be heeded, when it comes to readability? For example how to group variable declarations? Can I put all doubles in a line or is it just a matter of taste?
8. Any other advice or hints?

Recursion function:

void optimalPath(double E_tmp, double t_tmp, size_t k) {

    double dE, dt;
    size_t i, j; // i - next velocity index (main for loop)

    for (i=N_i; (i+1)>0; i--) {
        // Indexing works as follows
        // M x N x P // i, j, k --> mat[ i + M * (j + k * N)]
        dE = transEnergy[k + N_x * (v_tmp[k] +  i * N_v)];
        dt = transTime[k + N_x * (v_tmp[k] +  i * N_v)];

        // checking for energy lower than of current optimal path
        // and time less than desired time t_aim
        if ((E_tmp + dE) <= E_min && (t_tmp + dt) < t_aim){
            v_tmp[k+1] = i;

            // if path complete
            if (k == N_k){

                if (t_fin == -1){
                    t_fin = t_tmp;
                }

                // if total time of current path is closer to desired time 
                // or less than time of current optimal path
                if ((t_aim - (t_tmp + dt)) < (t_aim - t_fin) || (t_tmp + dt) < t_fin){
                    // save the current path as new optimal path
                    E_min = E_tmp + dE;
                    t_fin = t_tmp + dt;

                    for (j=0; j<=N_x; j++) {
                        v[j] = (uint16_T)v_tmp[j];
                    }
                    break;
                }
                break;
            }
            else {
                // proceeding to next segment k + 1
                optimalPath(E_tmp + dE, t_tmp + dt, k + 1);
            }
        }

    }
}

Rest of the file (recursion function is in the same file actually):

#include "mex.h"
#include "matrix.h"
#include <stdlib.h>
/*
 * optimalPath.c
 *
 * Finds the optimal path using
 * 
 * - minimum total energy given according to the transition energy
 *   given by first input matrix transEnergy (N_x, N_v, N_v)
 *
 * - maximum of available time t_aim according to the transition time
 *   given by second input matrix transTime (N_x, N_v, N_v)
 *
 * The calling syntax is:
 *
 *      v = optimalPath(transEnergy, transTime, t_aim)
 *
 * v is the vector of optimal velocities
 *
 * This function will take a bearable amount of time up to N_x ~ 13, e. g.
 * N_x:13 N_v:36/18 - Time: 15.36/0.3 s
 * 
 * No loop abortion condition in this version!
*/
double t_fin, E_min, inf, t_aim, *transEnergy, *transTime;
size_t *v_tmp, N_x, N_v, N_i, N_k;
uint16_T *v;

/* The gateway function */
void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
    /* VARIABLE DECLARATIONS */
    const mwSize *dim_array;
    size_t j;

    // init globals
    t_fin = -1;
    E_min = 1e15; 
    inf = mxGetInf();

    /* ARGUMENT VERIFICATION */
    if(nrhs!=3) {
        mexErrMsgIdAndTxt("MyToolbox:optimalPath:nrhs",
                          "Two inputs required.");
    }

    if(nlhs!=1) {
        mexErrMsgIdAndTxt("MyToolbox:optimalPath:nlhs",
                          "One output required.");
    }

    /* make sure the first two input arguments are matrices */
    if( !mxIsDouble(prhs[0]) || !mxIsDouble(prhs[1]) ||
         mxIsComplex(prhs[0]) || mxIsComplex(prhs[1])) {
        mexErrMsgIdAndTxt("MyToolbox:optimalPath:notDouble",
            "Input matrices must be type double.");
    }

    /* make sure the thrid input argument is scalar */
    if( !mxIsDouble(prhs[2]) || 
         mxIsComplex(prhs[2]) ||
         mxGetNumberOfElements(prhs[2])!=1 ) {
        mexErrMsgIdAndTxt("MyToolbox:optimalPath:notScalar",
                          "Input multiplier must be a scalar.");
    }

    /* INPUT ARGUMENTS */

    /* create a pointer to the real data in the input matrix  */
    transEnergy = mxGetPr(prhs[0]);
    transTime = mxGetData(prhs[1]);

    /* get dimensions of the input matrix */
    /* N_x corresponds to N_x-1 in the original function */    
    dim_array = mxGetDimensions(prhs[1]);

    N_x = (size_t)dim_array[0];
    N_v = (size_t)dim_array[1];
    N_i = N_v - 1;
    N_k = N_x - 1;

    /* get the value of the scalar input  */
    t_aim = mxGetScalar(prhs[2]);
    v_tmp = (size_t *) mxCalloc(N_x + 1, sizeof(size_t *));
    v_tmp[0] = 0;

    /* OUTPUT ARGUMENTS */

    /* create the output (single row) matrix, i. e. vector*/
    plhs[0] = mxCreateNumericMatrix(1, (int)N_x + 1, mxUINT16_CLASS, mxREAL);
    /* get a pointer to the real data in the output matrix */
    v = (uint16_T *)mxGetPr(plhs[0]);

    printf("N_x:%d N_v:%d t_aim:%.0f\n", dim_array[0], N_v, t_aim);
    /* call the computational routine */
    optimalPath(0, 0, 0);

    // increase each item of result with +1
    for (j=0; j<=N_x; j++) {
        v[j]++;
    }
}

Answer 1

Some comments on this (not necessarily in order):

I somehow tried to use size_t everywhere, so that almost no casting would neccessary. Does it acctually make a difference in performance?

Not in general, no. However, it can substantially reduce the number of warnings about losing siginificance and casting signed to unsigned messages. It generally makes sense to use the type that's expected.

The lowest value of i in the main for loop has to be equal to zero. But for (i=N_i; i>=0; i--) resulted in wrong behaviour (probably of the unsigned type of i) so I change it to for (i=N_i; (i+1)>0; i--) - is there any better way? Is such an extra addition negligible in terms of performance?

Your other option is to check with the max value for i to see if it has wrapped (effectlively testing i != -1 but as i is unsigned, that's not really recommended). I doubt this would affect performance much, though it would have an impact on readability and maintainability. Your other option is to do something like

int iteration;
size_t i;
...
i = N + 1;
for (iteration = 0; iteration != N; ++iteration)
{
    --i;
   ...
}

Making double dE, dt; local variables made the code faster. However I tried to pass some of the other variables to optimalPath to aviod accessing them as globals, but apparently it didn't give any speed up. Any general advice here?

You should always give variables the minimum possible scope (and you can declare variables inside blocks) - it makes maintenance easier. Only do otherwise if you can prove it makes an appreciable difference to performance.

Indexing with k + N_x * (v_tmp[k] + i * N_v) is done two times. Is it better to first compute the index, caching it and the use it?

Do a performance analysis. That's a complex expression, so from the point of view of maintainer, calculating once has a lot of benefits. But you can't really tell what effect it'll have on performance without measurement.

Does this code have any bottleneck or is time consumption evenly distributed? Would it be better to avoid recursion and solve the problem using a different approach or is recursion not that bad, depending on the conditions? In the given case the maximum recursion depth is the length of the path N_x (I think).

You have no idea what your compiler will do with the code. TO see how different algorithms will affect the performance, use tool to measure the performance.

In terms of recursion, per se it's not massively expensive unless you do it a lot. But it is generally more expensive than iteration as it involves copying arguments and calling and returning. However, it can be more difficult to understand what is going on with iteration, depending on what changes get made to the local variables. If you have a good compiler, it's possible that the recursion (esp if it is recognisable as tail recursion, which this one should be) will be turned into an iteration.

What should be heeded, when it comes to readability? For example how to group variable declarations? Can I put all doubles in a line or is it just a matter of taste? Any other advices or hints?

OK, this is pretty much personal, but

• Kernighan and Ritchies used run on braces MERELY to reduce the number of pages in the book, not as a recommended coding style. Put braces on separate lines
• Declare variables as late as possible and minimise their scope. It helps the reader (and the compiler) to know where the variables aren't used. For instance, dE and dT in optimalPath can be declared inside the loop and immediately initialised. As a corollary of this, you should avoid global variables/
• Avoid uninitialised variables where possible. This is trickier in C than C++ but it might make sense to declare a new block purely to limit the scope of a bunch of variables.
• Declare each item on a line by itself in a single statement. It's easier to check differences when doing maintenance and costs nothing.
• Choose an efficient algorithm (irrespective of iteration or recursion, an algorithm that takes O(n) iterations is going to be better than one that takes O(nLogn) iterations. Then measure performance and start optimising.
• Try and arrange recursive functions so that they can benefit from tail recursion optimisations.

Answer 2

A few things that jumped out at me that had not already been covered by another answer.

1. The following line looks wrong to me ...

   if ((t_aim - (t_tmp + dt)) < (t_aim - t_fin) || (t_tmp + dt) < t_fin)
   

   If I am following the logic correctly, it will always evaluate to TRUE unless

   (t_tmp + dt) == t_fin
   

   Here's why. Consider the following inequality...

   (t_aim - (t_tmp + dt)) < (t_aim - t_fin)
   

   Now subtract t_aim from both sides...

   -(t_tmp + dt) < -t_fin
   

   Simplify by multiplying both sides by -1...

   (t_tmp + dt) > t_fin
   

   What we wind up with is the following logic...

   if (((t_tmp + dt) > t_fin) || ((t_tmp + dt) < t_fin))
   

   Which is the same as...

   if ((t_tmp + dt) != t_fin)
   

   Again, I doubt this is the intent of the statement.

2. There is an implicit conversion between double and int in the following line:

   if (t_fin == -1)
   

   My guess is that you want...

   if (t_fin == -1.0)
   

   However, the calculation of t_fin might not hit exactly -1.0, so you might want it within a some range close to -1.0 instead.

   if ((t_fin > (-1.0 - EPSILON)) && (t_fin < (-1.0 + EPSILON)))
   

   where EPSILON is a small positive double.

3. If the routine optimalPath() is only used in this file, you may wish to try applying the static keyword. This can often help the compiler perform additional optimizations.