I’ve been writing some functions used to find paths between two types of triangles – alphas and betas. Alphas are triangles that have been in a zone we consider important, have an “interesting” value above a given threshold, and are “active”. Betas are essentially anything that isn’t an Alpha.

The position of the zone and the geometry of the model can change between invocations.

This is written in C++03, compiled into a MEX file (.mexw64) to be executed by MATLAB R2016B on a Linux machine. Those are all hard limits.

This code uses a good deal of functions and data from an external libraries and objects. However, most of the methods used are very simple array lookups, nothing performance-hindering.

Everything works correctly so far in testing, but performance has become a significant problem.

**The code:**

`// Doxygen block exists here // Various includes go here // Only needed because ultimately MATLAB needs an error code, not a C++ // exception #define SUCCESS 0 #define DYN_ALLOC_ERR 1 #define PTHREAD_ERR 2 /* * Design notes: I considered having a modified version of the geometry checking * section that just pushed _new_ alphas onto the global vector of alphas * and set a flag to cull old ones, but that actually seemed significantly less * efficient than just making a new one each time since all the same checks need * to be made regardless however, if validAlphas gets very large, this could * be highly inefficient due to push_back, so the idea of pushing on * just new ones and culling may actually be the best solution. * * Also: Could maintain a sum of validAlphas and use that to re-size it upon * re-creation to save some overhead from resizing, maybe? */ // The indices of two triangles who have a valid alpha and a beta that can // be seen from it. 120k paths * 8 bytes per pair = ~1GB on the heap. // No, ushort won't be enough here. struct ABPair { unsigned int alphaTriIndex; unsigned int betaTriIndex; }; // Useful for multithreading, stolen and modified from craytracer.h struct ThreadData { CRayTracer* rt; pthread_t threadId; // ID returned by pthread_create unsigned uThreadID; // Index std::vector<ABPair*> validPathsThread; // valid pairs that this thread // found unsigned int numTris; // Input argument, the number of // triangles in the mesh double distThreshold; // Input argument, the maximum // distance between triangles }; // Exceptions for experimentation class PThreadException: public std::exception { virtual const char* what() const throw() { return "Exception occured in a pthread_attr_init or pthread_create\n"; } }; class DynAllocationException: public std::exception { virtual const char* what() const throw() { return "Exception occured when attempting to malloc or calloc\n"; } }; // Note: Globals must exist here so that when the MEX file exits and goes // Back to MATLAB, this information is maintained. // (AFAIK, mexMakeMemoryPersistant() wouldn't make sense here. I can link // to discussions on that topic) // An indicator for every tri to tell if it has been removed (neccesary for // maintaining previous state to check if we need to call findPaths()) static bool* triActive = NULL; // A map from a given face to the element it resides in static unsigned int* faceToElementMap = NULL; // All valid paths. Must be maintained, because we don't know if // findPaths() will be called. It may not be if geometry hasnt changed. static std::vector<ABPair*> validPaths; // The previous state of what alphas were considered valid. Neccesary to see // if a change has occured that isn't geometry-based. static unsigned int* validAlphaIndex; // Not neccesary as a global if it's just getting re-generated each time. // However, if we do decide to append and cull instead of regenerating, this // needs to stay. static std::vector<unsigned int> validAlphas; // Needed so we can accurately determine alphas. // I removed this in the past, thinking it wasn't needed. As it turns out, // it's absolutely needed and very helpful. static bool* hasBeenInZone; // Useful everywhere CRayTracerClass* rayTracer = NULL; NanoRTWrapperClass nanoRTWrapper = NULL; // Function declarations // Not required, but prevents warnings depending on how functions are ordered // and call each other // (Including the mexFunction here would be redundant, as it exists in mex.h) void exitFcn(); bool isTriInZoneRadius(const unsigned int itri); bool checkForModelChanges(const unsigned int numTris, const float* iValues, const double iThreshold ); void initialize(const float* elemFace, const unsigned int numElems, const unsigned int facePerElMax, unsigned int* numTri, unsigned int* numFace ); void* findPathsThread(void *data); void findPathsThreadSpooler(const unsigned int numTris, const double distThreshold ); void mapFacesToElements(const float* elemFace, const unsigned int numElems, const unsigned int facePerElMax ); bool checkPairValid(const unsigned int i, const unsigned int j, const double distThreshold ); bool isTriAlpha(const unsigned int itri, const float* iValues, const double iThreshold ); void findPaths(const unsigned int numTris, const double distThreshold ); //mainfunc declaration goes here /** * @brief exitFcn - Cleans up malloc'd or calloc'd memory if someone in the * MATLAB script calls "clear mexfilename" or "clear all". */ void exitFcn() { //mexPrintf("exitFcn() called\n"); if(triActive) { free(triActive); } if(faceToElementMap) { free(faceToElementMap); } if(validAlphaIndex) { free(validAlphaIndex); } if(hasBeenInZone) { free(hasBeenInZone); } for(unsigned int i = 0; i < validPaths.size(); i++) { free(validPaths[i]); } } /** * @brief Checks if a given tri is currently in the zone's external radius. * Implementation stolen from CRayTracerClass::trace_inverseray_tri * Not sure if we need to raytrace, so I omitted it * @param itri - The index of the triangle to check * @return True if in the radius, false if not */ bool isTriInZoneRadius(const unsigned int itri) { //Omitted } /** * @brief Checks if the model has changed (either in terms of alphas or * geometry) and re-generates the vector of alphas * @param numTris - The number of triangles in the finite mesh * @param iValues - The ivalue at each node * @param iThreshold - The interesting value threshold beyond which an alpha * is interesting enough to be valid * @return True if the list of alphas or the geometry has changed, false if * neither have */ bool checkForModelChanges(const unsigned int numTris, const float* iValues, const double iThreshold ) { bool modelChanged = false; bool isAlpha; bool currentlyActive; // Two checks need to happen - geometry changes and for the list of valid // alphas to change // Also regenerates the vector of valid alphas from scratch as it goes for(unsigned int i = 0; i < numTris; i++) { // Active means it has 1 exposed face, not 2 (internal) or 0 (gone) currentlyActive = nanoRTWrapper->getTriActive(i); // Has the geometry changed? if(currentlyActive != triActive[i]) { modelChanged = true; triActive[i] = currentlyActive; } // Get whether this triangle is an alpha: isAlpha = isTriAlpha(i, iValues, iThreshold); // Triangle is a valid alpha now, but wasn't before if((isAlpha == true) && (validAlphaIndex[i] == false)) { validAlphaIndex[i] = true; modelChanged = true; } // Was valid before, is no longer valid now else if((isAlpha == false) && (validAlphaIndex[i] == true)) { validAlphaIndex[i] = false; modelChanged = true; //cullalphasFlag = true; } // Generating the set of all valid alphas if(isAlpha) { validAlphas.push_back(i); } } return modelChanged; } /** * @brief Initializes this MEX file for its first run * @param rt - A pointer to the raytracer object * @param numTris - The total number of triangles in the finite mesh * @param numFaces - The total number of faces in the finite mesh * @param elemFace - The map of elements to the faces that they have * @param numElems - The number of elements in the finite mesh * @param facePerElMax - The maximum number of faces per element */ void initialize(const float* elemFace, const unsigned int numElems, const unsigned int facePerElMax, unsigned int* numTri, unsigned int* numFace ) { DynAllocationException e; // Fetch number of tris and faces // Must be done every time, since we're storing locally and not globally // However: // They're never modified // They never change between calls from the MATLAB script // They're used frequently in many functions // I think that's a strong candidate for being a global unsigned int numTris = nanoRTWrapper->getTriCount(); *numTri = numTris; unsigned int numFaces = nanoRTWrapper->getFaceCount(); *numFace = numFaces; /* * Allocate some space for things we need to be persistent between runs of * this MEX file. And check that the allocation succeeded, of course. */ if(triActive == NULL) { if(NULL == (triActive = static_cast<bool*>(calloc(numTris, sizeof(bool)))) ) { throw e; } } if(hasBeenInZone == NULL) { if(NULL == (hasBeenInZone = static_cast<bool*>(calloc(numTris, sizeof(bool)))) ) { throw e; } } if(validAlphaIndex == NULL) { if(NULL == (validAlphaIndex = static_cast<unsigned int*>(calloc(numTris, sizeof(unsigned int)))) ) { throw e; } } if(faceToElementMap == NULL) { if(NULL == (faceToElementMap = static_cast<unsigned int*>(calloc(numFaces, sizeof(unsigned int)))) ) { throw e; } mapFacesToElements(elemFace, numElems, facePerElMax); } return; } /** * @brief Is something that can be used by pthread_create(). Threads will skip * over some of the work, and do isValidPair on others. Thus...multithreading. * @param data - The data structure that will hold the results and arguments */ void* findPathsThread(void *data) { struct ThreadData* thisThreadsData = static_cast<struct ThreadData*>(data); const unsigned uThreadID = thisThreadsData->uThreadID; const unsigned uNumThreads = rayTracer->m_uNumThreads; const double distThreshold = thisThreadsData->distThreshold; const unsigned int numTris = thisThreadsData->numTris; std::vector<ABPair*>& validPathsThread = thisThreadsData->validPathsThread; // Loop over all valid alphas for(unsigned int i = 0; i < validAlphas.size(); i++) { if ((i % uNumThreads) == uThreadID) { // Loop over all triangles (potential betas) for(unsigned int j = 0; j < numTris; j++) { if(checkPairValid(i, j, distThreshold)) { ABPair* temp = static_cast<ABPair*>(malloc(sizeof(ABPair))); temp->alphaTriIndex = validAlphas[i]; temp->betaTriIndex = j; validPathsThread.push_back(temp); } } } } return NULL; } /** * @brief Creates as many threads as the system has available, and then uses * pthread_create() to dish out the work of findPaths() * @param numTris - The number of triangles in the finite mesh * @param distThreshold - The maximum distance an alpha and beta can be * apart */ void findPathsThreadSpooler(const unsigned int numTris, const double distThreshold ) { std::vector<ThreadData> threadData(rayTracer->m_nProc); pthread_attr_t attr; int rc; PThreadException e; // I think this is checking to make sure something doesn't already exist, // not sure what though if((rc = pthread_attr_init(&attr))) { throw e; } // We know how many threads the system supports // So all this does is walk through an array of them and start them up for(unsigned uThread = 0; uThread < rayTracer->m_uNumThreads; uThread++) { ThreadData& data = threadData[uThread]; data.rt = rayTracer; data.uThreadID = uThread; data.numTris = numTris; data.distThreshold = distThreshold; if(rayTracer->m_uNumThreads > 1) { if((rc = pthread_create(&data.threadId, &attr, &findPathsThread, &data))) { throw e; } } else { findPathsThread(&data); } } // Join all threads for(unsigned uThread = 0; uThread < rayTracer->m_uNumThreads; uThread++) { std::vector<ABPair*>& validPathsThread = threadData[uThread].validPathsThread; if(rayTracer->m_uNumThreads > 1) { void* res; if((rc = pthread_join(threadData[uThread].threadId, &res))) { throw e; } } // validPathsThread is the set of ABPairs that this thread found // while validPaths is the globally maintained set of valid paths // Take each thread's results and merge it into the overall results validPaths.insert(validPaths.end(), validPathsThread.begin(), validPathsThread.end()); } return; } /* void cullAlphas() { for(unsigned int i = 0; i < validAlphas.size(); i++) { if(!isValidAlpha(validAlphas[i])) { validAlphas.erase(i); } } } */ /** * @brief Determines the elements that each face belongs to * @details the MATLAB script maintains a map of all faces per element. * This is the opposite of what we want. Accessing it linearly * walks by column, not by row. Also, MATLAB stores everything 1-indexed. * Finally, the MATLAB script left them stored as the default, which are singles. * @param elemFace - A MATLAB facePerElMax by numElems array, storing which * faces belong to each element (elements being the row number) * @param numElems - The total number of elements (rows) in the array * @param facePerElMax - The max number of faces per element (the number of * columns) */ void mapFacesToElements(const float* elemFace, const unsigned int numElems, const unsigned int facePerElMax ) { unsigned int i; // elemFace[0] = 1. We don't know how elemFace will be structured precisely, // so we need to keep going until we find a face in it that equals our number // of faces, since it's 1-indexed. for(i = 0; i < (numElems * facePerElMax); i++) { faceToElementMap[static_cast<unsigned int>(elemFace[i]) - 1] = (i % numElems); // Is the next face for that element a NaN? If so, we can skip it. Keep // skipping until the next element WON'T be NaN. // Don't cast here, as NaN only exists for floating point numbers, // not integers. while(isnan(elemFace[i + 1]) && ((i + 1) < (numElems * facePerElMax))) { i++; } } } /** * @brief checkPairValid - Checks if a pair of an alpha index (of validAlphas), * beta index form a valid path * @param i - Index into validAlphas * @param j - Index into all tris (potential beta) * @param distThreshold - The max distance the tri's centers can be apart * @return Whether the pair forms a valid path */ bool checkPairValid(const unsigned int i, const unsigned int j, const double distThreshold ) { double path_dist_sqrd; double path_dist; double alphaCoords[3]; double betaCoords[3]; nanort::Ray<double> ray; // If they're not an alpha currently, they must be a potential beta, // must also be alive if(!validAlphaIndex[j] && triActive[j]) { alphaCoords[0] = rayTracer->m_vecTriFixedInfo[validAlphas[i]].center.x(); alphaCoords[1] = rayTracer->m_vecTriFixedInfo[validAlphas[i]].center.y(); alphaCoords[2] = rayTracer->m_vecTriFixedInfo[validAlphas[i]].center.z(); betaCoords[0] = rayTracer->m_vecTriFixedInfo[j].center.x(); betaCoords[1] = rayTracer->m_vecTriFixedInfo[j].center.y(); betaCoords[2] = rayTracer->m_vecTriFixedInfo[j].center.z(); // Determine distance squared between alpha and beta // (x2-x1)^2 + (y2-y1)^2 +(z2-z1)^2 path_dist_sqrd = pow((betaCoords[0] - alphaCoords[0]), 2) + pow((betaCoords[1] - alphaCoords[1]), 2) + pow((betaCoords[2] - alphaCoords[2]), 2); // Doing this instead of doing the sqrt to save doing the sqrt when not // needed for performance if(path_dist_sqrd <= pow(distThreshold, 2)) { path_dist = sqrt(path_dist_sqrd); // Set up a nanort::Ray's origin, direction, and max distance ray.org[0] = alphaCoords[0]; // x ray.org[1] = alphaCoords[1]; // y ray.org[2] = alphaCoords[2]; // z ray.dir[0] = (betaCoords[0] - alphaCoords[0]) / path_dist; ray.dir[1] = (betaCoords[1] - alphaCoords[1]) / path_dist; ray.dir[2] = (betaCoords[2] - alphaCoords[2]) / path_dist; // TODO: Subtract some EPSILON here so it doesn't report a hit because it // hit the beta itself (assuming that's how it works) ray.max_t = path_dist; // Call ShootRay() to check if there is a path (calls nanoRT) if(!(nanoRTWrapper->shootRay(ray))) { return true; } else { // There's no path return false; } } else { // The distance is too far between alpha and beta return false; } } else { // The beta is either dead or currently an alpha return false; } } /** * @brief Determines if a given triangle is a valid alpha. * @param itri - The triangle index to check * @return True if it is an alpha, false if it is not */ bool isTriAlpha(const unsigned int itri, const float* iValues, const double iThreshold ) { double tri_avg_interesting; const unsigned int* tri_nodes; // Do the simple checks first, as it's more performant to do so // alternate consideration (acccuracy, wouldn't affect performance) //if(triActive[itri] && (hasBeenAlpha[itri] || isTriInZoneRadius(itri))) if(triActive[itri] && (hasBeenInZone[itri] || isTriInZoneRadius(itri))) { // Retrieve the average iValue of this triangle tri_nodes = nanoRTWrapper->getTriNodes(itri); tri_avg_interesting = (iValues[tri_nodes[0]] + iValues[tri_nodes[1]] + iValues[tri_nodes[2]]) / 3; if(tri_avg_interesting > iThreshold) { return true; } } return false; } /** * @brief Uses the raytracer object's current state as well as arguments to * generate pairs of unobstructed paths between alphas and betas. * @param numTris - The number of triangles in the finite element mesh * @param distThreshold - The max distance an alpha and beta pair can be * apart before not being considered in calculations */ void findPaths(const unsigned int numTris, const double distThreshold ) { // This function once held more importance, but yes, it could be omitted // at this point. // Spool up some threads to take care of the work try { findPathsThreadSpooler(numTris, distThreshold); } catch(DynAllocationException& e) { throw e; } return; } // Doxygen header (omitted) int mainFunc(args) { // Initialize the program if we're on a first run try { initialize(elemFace, numElems, facePerElMax, &numTris, &numFaces); } catch(DynAllocationException& e) { return DYN_ALLOC_ERR; } // Need to see if we need to call findPaths if(checkForModelChanges(numTris, iValues, iThreshold)) { // Remove old list of valid paths for(unsigned int i = 0; i < validPaths.size(); i++) { free(validPaths[i]); } validPaths.clear(); try { findPaths(numTris, distThreshold); } catch(PThreadException& e) { return PTHREAD_ERR; } } //mexPrintf("Number of valid paths: %d\n", validPaths.size()); /* * Walk over all paths and do some calculations */ return SUCCESS; } // Doxygen block goes here, omitted void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) { // Register exit function mexAtExit(exitFcn); // Prep for writing out results // Checking to make sure # of arguments was right from MATLAB // Input argument handling to convert from mxArrays to double*, float*, etc //*errcode = mainFunc(some args) // retrieve execution time in clock cycles, convert to seconds, print // Put the outputs in plhs } `

**Callgraph(?):**

This isn’t exactly a callgraph, but it might be useful to get an idea of the flow of the program.

**The Problem: Performance**

For medium-size models (104k tris, 204k faces, 51k elems) it can take up to a couple seconds for this to complete, even though the worst of it is multi-threaded on a powerful 4C/8T machine. (roughly 100*104k size loop)

For any models where the number of alphas is very large (50K) it can take up to *three minutes* for a single execution to complete because of how large that double-nested for loop must become. (50k^2 size loop)

**Possible optimizations:**

An optimization that may be worth considering is creating a sphere around all alphas, and use its center and radius to cull the remaining triangles down to a smaller size based on the threshold distance. However, the benefit of this is extremely variable and may actually be zero for smaller meshes. And to create the sphere, we need to find the two triangle centers that are the furthest apart, an O(alphas^2) operation…very slow if there are many. And there’s no easy way to tell in advance if using this technique would be beneficial for the given state of the mesh, so it’s not easy to toggle on and off as needed.

As well, it’s possible something can be done with nanoRT’s BVHs, but I’m not very familiar with BVH’s or what they’d let me do in this

**Note: How it’s being used:**

The MATLAB script will likely call this many times. In small models, it may finish its own loop within tenths of a second and then call ours again. In larger ones, there may be half a second between calls. In total, this may be called hundreds of times.

**Note: How it’s being built:**

This isn’t built using the `MEX`

command in MATLAB, nor by using Visual Studio. Instead, g++ is used to create an object file (.o) and then g++ is used again to create the .mexw64 file in a method I’m not entirely familiar with. (This is also a hard limit I cannot touch)

I occasionally compiled with very aggressive warnings enabled to catch things like sign conversion, promotions, bad casts, etc.

**Profiling:**

I would love to be able to profile this code. However, it seems impossible. MEX files built using `MEX`

command in MATLAB can be done. MEX files compiled in Visual Studio can be profiled. But we’re not doing either of those, and so when I try to profile with either MATLAB or Visual Studio, it just doesn’t work.

Even if I could, I don’t think it would reveal anything surprising. The numbers we’re working with are large, so the double-nested loop at the core of it grows very large.

If need be, I could stick a ton of clock() calls around the start/end of every method to get some more info, though that’s not precise. (Line-by-line)

**Final Notes:**

I’m fresh out of college and this is my first real production code. I’m more familiar with C than C++, so I’m sure that bled into the code. Suggestions for style are always welcome, though performance is the most major concern here.