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:
- I somehow tried to use
size_t
everywhere, so that almost no casting would necessary. Does it actually make a difference in performance? - The lowest value of
i
in the mainfor
loop has to be equal to zero.
Butfor (i=N_i; i>=0; i--)
resulted in wrong behaviour (probably of the unsigned type ofi
) so I change it tofor (i=N_i; (i+1)>0; i--)
- is there any better way? Is such an extra addition negligible in terms of performance? - Making
double dE, dt;
local variables made the code faster. However I tried to pass some of the other variables tooptimalPath
to aviod accessing them as globals, but apparently it didn't give any speed up. Any general advice here? - 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? - 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). - 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 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]++;
}
}