Interpolation (WIP)

[adamkewley]

This is still a WIP, but hopefully will clean up + accelerate the implementation. Changes:

• Moved the evaluations into the .osim file
  • The osim file now contains the evaluations, rather than some external file
  • This simplifies the finalizeFromProperties and finalizeConnections phases of the model, because they used to rely on loading external files, which can fail if (e.g.) the model is copied in-memory (and, therefore, has no backing file from which the parent directory can be computed) or if the model file is opened by a GUI running in a different working directory (which requires careful handling of the current working directory during the finalization phase, which causes other bugs)
  • This also simplifies the command-line tool, because it now has a 1:1 relationship between input and output (osim goes in, osim goes out - no other files)
  • This does mean that FBP-based model files can become very large, but (imho) this is outweighed by the advantages of having everything in one file
  • The data looks something like:

<FunctionBasedPath name="functionbasedpath">
    <!-- other props -->
    <FunctionBasedPathDiscretizationSet>
        <objects>
            <FunctionBasedPathDiscretization name="functionbasedpathdiscretization">
                <!--The absolute path, in the model, to the OpenSim::Coordinate that this discretization was produced from-->
                <coordinate_abspath>/jointset/r_shoulder/r_shoulder_elev</coordinate_abspath>
                <!--The lowest OpenSim::Coordinate value that was used for discretization-->
                <x_begin>-1.6231562043547265</x_begin>
                <!--The highest OpenSim:::Coordinate value that was used for discretization-->
                <x_end>1.6755160819145563</x_end>
                <!--The number of evenly-spaced OpenSim::Coordinate values between [x_begin, x_end] (inclusive) that were used for discretization of the OpenSim::Coordinate. E.g. [x_begin, 1*(x_begin+((x_end-x_begin)/3)), 2*(x_begin+((x_end-x_begin)/3)), x_end]-->
                <num_points>64</num_points>
            </FunctionBasedPathDiscretization>
        </objects>
    </FunctionBasedPathDiscretizationSet>
    <Evaluations>
        0.27785555660465677 0.27672518770854371 0.2755762685216484 <!--- a lot more numbers -->
    </Evaluations>
</FunctionBasedPath>

• Added a variety of comments that try to explain each step of the process
  • This includes restructuring the implementation file to use standalone functions over simple datatypes
• Methods that used to warn on each call (e.g. generateDecorations) now warn once per application boot
  • This is to prevent the situation where a GUI calling generateDecorations repeatably would have its log spammed with the same warning message
• Externalized all path fitting parameters into a caller-provided FunctionBasedPath::FittingParams struct
  • This is so that external users can tweak with the fitting params without having to recompile the codebase. It enables performing experiments on different gridsizes, probing steps, maximum coordinates, etc. with scripts.
  • The main utility of this is that it lets us performance-tweak the PBP-to-FBP step
  • Here's the struct:

  struct FittingParams final {

      // maximum coords that can affect the given PointBasedPath
      //
      // if this is higher, more paths may be eligible for
      // PointBasedPath --> FunctionBasedPath conversion, because some paths may be
      // affected by more coordinates than other paths. However, be careful. Increasing
      // this also *significantly* increases the memory usage of the function-based fit
      //
      // must be 0 < v <= 16, or -1 to mean "use a sensible default"
      int maxCoordsThatCanAffectPath;

      // number of discretization steps to use for each coordinate during the "probing
      // phase"
      //
      // in the "probing phase", each coordinate is set to this number of evenly-spaced
      // values in the range [getRangeMin()..getRangeMax()] (inclusive) to see if changing
      // that coordinate has any affect on the path. The higher this value is, the longer
      // the probing phase takes, but the higher chance it has of spotting a pertubation
      //
      // must be >0, or -1 to mean "use a sensible default"
      int numProbingDiscretizations;

      // minimum amount that the moment arm of the path must change by during the "probing phase"
      // for the coorinate to be classified as affecting the path
      //
      // must be >0, or <0 to mean "use a sensible default"
      double minProbingMomentArmChange;

      // the number of discretization steps for each coordinate that passes the "probing phase" and,
      // therefore, is deemed to affect the input (point-based) path
      //
      // this is effectively "grid granulaity". More discretizations == better fit, but it can increase
      // the memory usage of the fit significantly. Assume the path is parameterized as an n-dimensional
      // surface. E.g. if you discretize 10 points over 10 dimensions then you may end up with
      // 10^10 datapoints (ouch).
      //
      // must be >0, or -1 to mean "use a sensible default"
      int numDiscretizationStepsPerDimension;

      FittingParams();
  };

• Changed signature of FunctionBasedPath::fromPointBasedPath to return a std::unique_ptr<FunctionBasedPath> rather than a FunctionBasedPath (and also, so that it takes the FittingParams, above)

  • This allows the function to fail to fit the path and return nullptr, which the caller can handle as "can't fit"
  • This allows the user to use FittingParams to limit how "far" the fitting implementation should go. E.g. if it's known that some muscles are affected by many coordinates, it's probably best to leave it as a PointBasedPath (this is likely to be the case for larger models)

• Refactored the algorithm to maximize performance + readability:

  • Functions renamed to reflect intent (e.g. getInterpDer --> Impl_GetPathLengthDerivative, getInterp --> Impl_GetPathLength)
  • Added a variety of comments that try to explain each step of the algorithm, with appropriate links
  • Renamed a variety of variables to reflect intent (e.g. dS --> discretizations, n --> closestDiscretizationSteps)
  • Removed all runtime memory allocations and replaced them with static arrays that have an upper bound (g_MaxCoordsThatCanBeInterpolated)
  • Removed most of the auxiliary datastructures that were held in Interpolation. Many of them can be computed on-the-fly by the CPU (e.g. `gridsize == (discretization.end-discretization.start)/(discretization.steps-1)). Most of these datastructures are now created + populated at runtime in stack-based arrays that have no allocation overhead
  • Replaced runtime data lookups with closed-form mathematical equivalents. E.g. instead of finding the closest discretization by iterating over a collection of discretization steps, the implementation computes the whole number part of v = begin + n*stepSize (solved for n)
  • Replaced multiply-add accumulations (e.g. acc += x*y) inside loops with std::fma (e.g. acc = std::fma(x, y, acc);). This reduces error propagation when accumulating small values onto larger values in a loop (e.g. z += beta*v)
  • Micro-optimized the implementation based on feedback from perf and valgrind. This mostly involved doing things like switching pow calls for in-source equivalents (e.g. pow(frac, 2) --> frac*frac). Some compilers do not optimize pow (yet)

[adamkewley]

This is still work-in-progress. I need to see if it's actually as correct as the original implementation (the algorithm is the same, but it has been majorly refactored). However, anecdotally, it does seem to be a little bit (~30 %) faster:

Before (from here):

Input Model (point-based paths): 
    avg. time (ms) = 439 
    integration steps attempted = 764 
    integration steps taken = 531)
Output Model (function-based paths): 
    avg. time (ms) = 194 
    integration steps attempted = 729 
    integration steps taken = 520)

After:

Input Model (point-based paths): 
    avg. time (ms) = 415 
    integration steps attempted = 764 
    integration steps taken = 531)
Output Model (function-based paths): 
    avg. time (ms) = 134 
    integration steps attempted = 746 
    integration steps taken = 512)

Again, though, I need to ensure that it's actually correct.

[adamkewley]

I fixed a few edge-case bugs that pop up in real simulations (e.g. ToyDropLanding), where coordinate values can fall onto the edges of the interpolation range.

What can happen is that the real coordinate value falls out of [min, max] and then the implementations (both this new closed-form one and the existing std::find one) would assign the nearest discretization point as 0 or npoints-1 (or something). The problem with this is that the subsequent steps are testing -1, 0, 1, 2 relative to that point. If you are on 0 and test -1 then you'll perform an out-of-bounds memory access, same goes for npoints-1+2.

The other thing that was changed was the default number of discretization points. 64 (80+3 in the original) is faaaar too high for any realistic model. Full-sized models will contain paths affected by (e.g.) 4 coordinates. That's 80^4 (~40 M) evaluations, which will blow up most machines.

In addition, the implementation has been hard-coded to only handling paths with up to 8 coordinates that affect the path. This limit can be reduced by changing the fitting parameters, but cannot be increased. Later versions can increase this limit by changing a single global. Muscles that are affected by >4 coordinates will produce much larger datasets and this particular algorithm will probably explode under those conditions. E.g. even with only 8 datapoints per dimension, it's going require 8^8 (~16 M) evaluations.

The FunctionBasedPathTool now enables setting fitting parameters externally. This is so that downstream code (test suites, benchmarking suties, etc.) can modify the fitting parameters at runtime (e.g. as part of testing how accuracy is affected by the fitting parameters).

The test suite now has a new test that collects all the datapoints from a PBP and and FBP into two std::vectors for comparison. This enables performing FBP-to-PBP comparisons in C++, which should make testing the relative quality of the fit easier (rather than relying on external MATLAB scripts, or whatever).

I noticed the derivative implementation isn't continuous. Rather, it takes the value two nearby points and performs discrete derivatization. The reason I mention this is because that implementation doubles the amount of computation (the derivative is asked for a lot at runtime). It also causes weird runtime bugs if h is invalid. I had to make h smaller for the Arm26 example to resemble the point-based model, for example. This kind of tweaking isn't going to work at scale, and the final implementation might need to have an actual derivative implementation that doesn't rely on knowing h.

[adamkewley]

I've committed the code I developed as-is, @joris997 , so you'll need to remove the hardcoded paths.

The script I'm using to test accuracy is:

import pandas
import matplotlib
import matplotlib.pyplot as plt

df = pandas.read_csv("datatable.tsv", sep="\t", index_col=0)
df.plot()
plt.show()

Which requires pip packages matplotlib and pandas. I run it from my dekstop with python3 plot.py