Add ExpressionBasedFunction

Bindings/OpenSimHeaders_common.h

@@ -22,6 +22,7 @@
 #include <OpenSim/Common/Event.h>
 #include <OpenSim/Common/Exception.h>
 #include <OpenSim/Common/ExperimentalSensor.h>
+#include <OpenSim/Common/ExpressionBasedFunction.h>
 #include <OpenSim/Common/FileAdapter.h>
 #include <OpenSim/Common/Function.h>
 #include <OpenSim/Common/FunctionSet.h>

Bindings/common.i

@@ -60,6 +60,7 @@
 %include <OpenSim/Common/Sine.h>
 %include <OpenSim/Common/PolynomialFunction.h>
 %include <OpenSim/Common/MultivariatePolynomialFunction.h>
+%include <OpenSim/Common/ExpressionBasedFunction.h>
 
 %include <OpenSim/Common/SmoothSegmentedFunctionFactory.h>
 %include <OpenSim/Common/SmoothSegmentedFunction.h>

CHANGELOG.md

@@ -18,6 +18,7 @@ v4.6
 - Added `Output`s to `ExpressionBasedCoordinateForce`, `ExpressionBasedPointToPointForce`, and `ExpressionBasedBushingForce` for accessing force values. (#3872)
 - `PointForceDirection` no longer has a virtual destructor, is `final`, and its `scale` functionality
   has been marked as `[[deprecated]]` (#3890)
+- Added `ExpressionBasedFunction` for creating `Function`s based on user-defined mathematical expressions. (#3892)
 
 v4.5.1
 ======

OpenSim/Common/CMakeLists.txt

@@ -16,7 +16,7 @@ OpenSimAddLibrary(
     KIT Common
     AUTHORS "Clay_Anderson-Ayman_Habib-Peter_Loan"
     # Clients of osimCommon need not link to ezc3d.
-    LINKLIBS PUBLIC ${Simbody_LIBRARIES} spdlog::spdlog
+    LINKLIBS PUBLIC ${Simbody_LIBRARIES} spdlog::spdlog osimLepton
              PRIVATE ${ezc3d_LIBRARY}
     INCLUDES ${INCLUDES}
     SOURCES ${SOURCES}

OpenSim/Common/ExpressionBasedFunction.cpp

@@ -0,0 +1,127 @@
+/* -------------------------------------------------------------------------- *
+ *                   OpenSim:  ExpressionBasedFunction.cpp                    *
+ * -------------------------------------------------------------------------- *
+ * The OpenSim API is a toolkit for musculoskeletal modeling and simulation.  *
+ * See http://opensim.stanford.edu and the NOTICE file for more information.  *
+ * OpenSim is developed at Stanford University and supported by the US        *
+ * National Institutes of Health (U54 GM072970, R24 HD065690) and by DARPA    *
+ * through the Warrior Web program.                                           *
+ *                                                                            *
+ * Copyright (c) 2005-2024 Stanford University and the Authors                *
+ * Author(s): Nicholas Bianco                                                 *
+ *                                                                            *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
+ * not use this file except in compliance with the License. You may obtain a  *
+ * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
+ *                                                                            *
+ * Unless required by applicable law or agreed to in writing, software        *
+ * distributed under the License is distributed on an "AS IS" BASIS,          *
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
+ * See the License for the specific language governing permissions and        *
+ * limitations under the License.                                             *
+ * -------------------------------------------------------------------------- */
+
+#include "ExpressionBasedFunction.h"
+
+#include <lepton/ExpressionProgram.h>
+#include <lepton/ParsedExpression.h>
+#include <lepton/Parser.h>
+#include <lepton/Exception.h>
+
+using namespace OpenSim;
+
+class SimTKExpressionBasedFunction : public SimTK::Function {
+public:
+    SimTKExpressionBasedFunction(const std::string& expression, 
+            const std::vector<std::string>& variables) :
+        m_expression(expression), m_variables(variables) {
+
+        // Check that the variable names are unique.
+        std::set<std::string> uniqueVariables;
+        for (const auto& variable : m_variables) {
+            if (!uniqueVariables.insert(variable).second) {
+                OPENSIM_THROW(Exception, 
+                        fmt::format("Variable '{}' is defined more than once.", 
+                        variable));
+            }
+        }
+
+        // Create the expression programs for the value and its derivatives.
+        Lepton::ParsedExpression parsedExpression = 
+                Lepton::Parser::parse(m_expression).optimize();
+        m_valueProgram = parsedExpression.createProgram();
+
+        for (int i = 0; i < static_cast<int>(m_variables.size()); ++i) {
+            Lepton::ParsedExpression diffExpression = 
+                    parsedExpression.differentiate(m_variables[i]).optimize();
+            m_derivativePrograms.push_back(diffExpression.createProgram());
+        }
+
+        try {
+            std::map<std::string, double> vars;
+            for (int i = 0; i < static_cast<int>(m_variables.size()); ++i) {
+                vars[m_variables[i]] = 0;
+            }
+            m_valueProgram.evaluate(vars);
+
+            for (int i = 0; i < static_cast<int>(m_variables.size()); ++i) {
+                m_derivativePrograms[i].evaluate(vars);
+            }
+        } catch (Lepton::Exception& ex) {
+            std::string msg = ex.what();
+            std::string undefinedVar = msg.substr(32, msg.size() - 32);
+            OPENSIM_THROW(Exception, 
+                    fmt::format("Variable '{}' is not defined. Use "
+                    "setVariables() to explicitly define this variable. Or, "
+                    "remove it from the expression.", undefinedVar));
+        }
+    }
+
+    SimTK::Real calcValue(const SimTK::Vector& x) const override {
+        OPENSIM_ASSERT(x.size() == static_cast<int>(m_variables.size()));
+        std::map<std::string, double> vars;
+        for (int i = 0; i < static_cast<int>(m_variables.size()); ++i) {
+            vars[m_variables[i]] = x[i];
+        }
+        return m_valueProgram.evaluate(vars);
+    }
+
+    SimTK::Real calcDerivative(const SimTK::Array_<int>& derivComponents, 
+            const SimTK::Vector& x) const override {
+        OPENSIM_ASSERT(x.size() == static_cast<int>(m_variables.size()));
+        OPENSIM_ASSERT(derivComponents.size() == 1);
+        if (derivComponents[0] < static_cast<int>(m_variables.size())) {
+            std::map<std::string, double> vars;
+            for (int i = 0; i < static_cast<int>(m_variables.size()); ++i) {
+                vars[m_variables[i]] = x[i];
+            }
+            return m_derivativePrograms[derivComponents[0]].evaluate(vars);
+        }
+        return 0.0;
+    }
+
+    int getArgumentSize() const override { 
+        return static_cast<int>(m_variables.size()); 
+    }
+    int getMaxDerivativeOrder() const override { return 1; }
+    SimTKExpressionBasedFunction* clone() const override {
+        return new SimTKExpressionBasedFunction(*this);
+    }
+
+private:
+    std::string m_expression;
+    std::vector<std::string> m_variables;
+    Lepton::ExpressionProgram m_valueProgram;
+    std::vector<Lepton::ExpressionProgram> m_derivativePrograms;
+};
+
+SimTK::Function* ExpressionBasedFunction::createSimTKFunction() const {
+    OPENSIM_THROW_IF_FRMOBJ(get_expression().empty(), Exception,
+            "The expression has not been set. Use setExpression().")
+
+    std::vector<std::string> variables;
+    for (int i = 0; i < getProperty_variables().size(); ++i) {
+        variables.push_back(get_variables(i));
+    }
+    return new SimTKExpressionBasedFunction(get_expression(), variables);
+}

OpenSim/Common/ExpressionBasedFunction.h

@@ -0,0 +1,134 @@
+#ifndef OPENSIM_EXPRESSION_BASED_FUNCTION_H_
+#define OPENSIM_EXPRESSION_BASED_FUNCTION_H_
+/* -------------------------------------------------------------------------- *
+ *                    OpenSim:  ExpressionBasedFunction.h                     *
+ * -------------------------------------------------------------------------- *
+ * The OpenSim API is a toolkit for musculoskeletal modeling and simulation.  *
+ * See http://opensim.stanford.edu and the NOTICE file for more information.  *
+ * OpenSim is developed at Stanford University and supported by the US        *
+ * National Institutes of Health (U54 GM072970, R24 HD065690) and by DARPA    *
+ * through the Warrior Web program.                                           *
+ *                                                                            *
+ * Copyright (c) 2005-2024 Stanford University and the Authors                *
+ * Author(s): Nicholas Bianco                                                 *
+ *                                                                            *
+ * Licensed under the Apache License, Version 2.0 (the "License"); you may    *
+ * not use this file except in compliance with the License. You may obtain a  *
+ * copy of the License at http://www.apache.org/licenses/LICENSE-2.0.         *
+ *                                                                            *
+ * Unless required by applicable law or agreed to in writing, software        *
+ * distributed under the License is distributed on an "AS IS" BASIS,          *
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   *
+ * See the License for the specific language governing permissions and        *
+ * limitations under the License.                                             *
+ * -------------------------------------------------------------------------- */
+
+#include "osimCommonDLL.h"
+#include "Function.h"
+namespace OpenSim {
+
+/** 
+ * A function based on a user-defined mathematical expression.
+ * 
+ * This class allows users to define a function based on a mathematical
+ * expression (e.g., "x*sqrt(y-8)"). The expression can be a function of any 
+ * number of independent variables. The expression is parsed and evaluated using 
+ * the Lepton library.
+ * 
+ * Set the expression using setExpression(). Any variables used in the 
+ * expression must be explicitly defined using setVariables(). This 
+ * implementation allows computation of first-order derivatives only.
+ * 
+ * # Creating Expressions
+ * 
+ * Expressions can contain variables, constants, operations, parentheses, commas, 
+ * spaces, and scientific "e" notation. The full list of supported operations is: 
+ * sqrt, exp, log, sin, cos, sec, csc, tan, cot, asin, acos, atan, sinh, cosh, 
+ * tanh, erf, erfc, step, delta, square, cube, recip, min, max, abs, +, -, *, /, 
+ * and ^. 
+ */
+class OSIMCOMMON_API ExpressionBasedFunction : public Function {
+    OpenSim_DECLARE_CONCRETE_OBJECT(ExpressionBasedFunction, Function);
+
+public:
+//==============================================================================
+// PROPERTIES
+//==============================================================================
+    OpenSim_DECLARE_PROPERTY(expression, std::string, 
+            "The mathematical expression defining this Function.");
+    OpenSim_DECLARE_LIST_PROPERTY(variables, std::string, 
+            "The independent variables used by this Function's expression. "
+            "In XML, variable names should be space-separated.");
+
+//==============================================================================
+// METHODS
+//==============================================================================
+
+    /** Default constructor. */
+    ExpressionBasedFunction() { constructProperties(); }
+
+    /** Convenience constructor.
+     *  
+     * @param expression The expression that defines this Function.
+     * @param variables The independent variable names of this expression.
+     */
+    ExpressionBasedFunction(std::string expression, 
+            const std::vector<std::string>& variables) {
+        constructProperties();
+        set_expression(std::move(expression));
+        setVariables(variables);
+    }
+
+    /**
+     * The mathematical expression that defines this Function. The expression 
+     * should be a function of the variables defined via setVariables().
+     * 
+     * @note The expression cannot contain any whitespace characters.
+     */
+    void setExpression(std::string expression) {
+        set_expression(std::move(expression));
+    }
+    /// @copydoc setExpression()
+    const std::string& getExpression() const {
+        return get_expression();
+    }
+
+    /**
+     * The independent variable names of this expression. The variables names 
+     * should be unique and should be comprised of alphabetic characters or any 
+     * characters not reserved by Lepton (i.e., +, -, *, /, and ^). Variable
+     * names can contain numbers as long they do not come first in the name
+     * (e.g., "var0"). The input vector passed to calcValue() and 
+     * calcDerivative() should be in the same order as the variables defined
+     * here.
+     */
+    void setVariables(const std::vector<std::string>& variables) {
+        for (const auto& var : variables) {
+            append_variables(var);
+        }
+    }
+    /// @copydoc setVariables()
+    std::vector<std::string> getVariables() const {
+        std::vector<std::string> variables;
+        for (int i = 0; i < getProperty_variables().size(); ++i) {
+            variables.push_back(get_variables(i));
+        }
+        return variables;
+    }
+
+    /**
+     * Return a pointer to a SimTK::Function object that implements this
+     * function.
+     */
+    SimTK::Function* createSimTKFunction() const override;
+
+private:
+    void constructProperties() {
+        constructProperty_expression("");
+        constructProperty_variables();
+    }
+};
+
+} // namespace OpenSim
+
+#endif  // OPENSIM_EXPRESSION_BASED_FUNCTION_H_

OpenSim/Common/RegisterTypes_osimCommon.cpp

@@ -41,6 +41,7 @@
 #include "MultiplierFunction.h"
 #include "PolynomialFunction.h"
 #include "MultivariatePolynomialFunction.h"
+#include "ExpressionBasedFunction.h"
 
 #include "SignalGenerator.h"
 
@@ -88,6 +89,7 @@ OSIMCOMMON_API void RegisterTypes_osimCommon()
     Object::registerType( MultiplierFunction() );
     Object::registerType( PolynomialFunction() );
     Object::registerType( MultivariatePolynomialFunction() );
+    Object::registerType( ExpressionBasedFunction() );
 
     Object::registerType( SignalGenerator() );