Change finite-support mju_sigmoid from a two-quadratic construction to a single quintic with zero 1st and 2nd derivatives at 0 and 1, the boundary of the smooth-support segment.

Note: this function is currently unused in the MuJoCo codebase.
PiperOrigin-RevId: 529969804
Change-Id: I2bddc72916eff67625b806510935bc7f3ce5c770

doc/APIreference/functions.rst

@@ -3160,7 +3160,7 @@ mju_sigmoid
 
 .. mujoco-include:: mju_sigmoid
 
-Sigmoid function over 0<=x<=1 constructed from half-quadratics.
+Sigmoid function over 0<=x<=1 using quintic polynomial.
 
 .. _Derivatives-api:
 

doc/changelog.rst

@@ -20,13 +20,23 @@ General
 ^^^^^^^
 
 - Added :ref:`mjd_inverseFD` for finite-differenced inverse-dynamics derivatives.
-- Added functions for operations on banded-then-dense "arrowhead" matrices. Such matrices are
+- Added functions for operations on banded-then-dense "arrowhead" matrices. Such matrices are common when doing direct
-  common when doing direct trajectory optimization.
+  trajectory optimization. See :ref:`mju_cholFactorBand` documentation for details.
-  See :ref:`mju_cholFactorBand` documentation for details.
 - Added :ref:`mj_multiRay` function for intersecting multiple rays emanating from a single point.
   This is significantly faster than calling :ref:`mj_ray` multiple times.
 - Increased ``mjMAXUIITEM`` (maximum number of UI elements per section in Simulate) to 100.
 - Added :ref:`documentation<exProvider>` for resource providers.
+- Changed the formula for :ref:`mju_sigmoid`, a finite-support sigmoid :math:`s \colon \mathbf R \rightarrow [0, 1]`.
+  Previously, the smooth part consisted of two stitched quadratics, once continuously differentiable.
+  It is now a single quintic, twice continuously differentiable:
+
+  .. math::
+     s(x) =
+     \begin{cases}
+        0,                    &       & x \le 0  \\
+        6x^5 - 15x^4 + 10x^3, & 0 \lt & x \lt 1  \\
+        1,                    & 1 \le & x \qquad
+     \end{cases}
 
 Version 2.3.5 (April 25, 2023)
 ------------------------------

doc/programming/extension.rst

@@ -301,6 +301,7 @@ Resource providers work via callbacks:
 - :ref:`mjfGetResourceDir<mjfGetResourceDir>`: This callback is optional and is used to extract the directory from a
   resource name.  For example, the resource name ``http://www.example.com/myasset.obj`` would have
   ``http://www.example.com/`` as its directory.
+
 .. _exProviderUsage:
 
 Usage

include/mujoco/mujoco.h

@@ -1216,7 +1216,7 @@ MJAPI mjtNum mju_Halton(int index, int base);
 // Call strncpy, then set dst[n-1] = 0.
 MJAPI char* mju_strncpy(char *dst, const char *src, int n);
 
-// Sigmoid function over 0<=x<=1 constructed from half-quadratics.
+// Sigmoid function over 0<=x<=1 using quintic polynomial.
 MJAPI mjtNum mju_sigmoid(mjtNum x);
 
 

introspect/functions.py

@@ -7752,7 +7752,7 @@ FUNCTIONS: Mapping[str, FunctionDecl] = dict([
                  type=ValueType(name='mjtNum'),
              ),
          ),
-         doc='Sigmoid function over 0<=x<=1 constructed from half-quadratics.',
+         doc='Sigmoid function over 0<=x<=1 using quintic polynomial.',
      )),
     ('mjd_transitionFD',
      FunctionDecl(

src/engine/engine_util_misc.c

@@ -1156,7 +1156,7 @@ char* mju_strncpy(char *dst, const char *src, int n) {
 
 
 
-// Sigmoid function over 0<=x<=1 constructed from half-quadratics.
+// sigmoid function over 0<=x<=1 using quintic polynomial
 mjtNum mju_sigmoid(mjtNum x) {
   // fast return
   if (x<=0) {
@@ -1165,15 +1165,8 @@ mjtNum mju_sigmoid(mjtNum x) {
   if (x>=1) {
     return 1;
   }
-  if (x==0.5) {
-    return 0.5;
-  }
 
-  // lower part
+  // sigmoid: f(x) = 6*x^5 - 15*x^4 + 10*x^3
-  if (x<0.5) {
+  // solution of f(0) = f'(0) = f''(0) = 0, f(1) = 1, f'(1) = f''(1) = 0
-    return 2*x*x;
+  return x*x*x * (3*x * (2*x - 5) + 10);
-  }
-
-  // higher part
-  return (1 - 2*(1-x)*(1-x));
 }

src/engine/engine_util_misc.h

@@ -123,10 +123,10 @@ MJAPI void mju_insertionSortInt(int* list, int n);
 // Halton sequence
 MJAPI mjtNum mju_Halton(int index, int base);
 
-// Call strncpy, then set dst[n-1] = 0.
+// call strncpy, then set dst[n-1] = 0
 MJAPI char* mju_strncpy(char *dst, const char *src, int n);
 
-// Sigmoid function over 0<=x<=1 constructed from half-quadratics.
+// sigmoid function over 0<=x<=1 using quintic polynomial
 MJAPI mjtNum mju_sigmoid(mjtNum x);
 
 #ifdef __cplusplus

test/engine/engine_util_misc_test.cc

@@ -18,12 +18,14 @@
 #include <gtest/gtest.h>
 #include <mujoco/mjdata.h>
 #include <mujoco/mujoco.h>
+#include "src/engine/engine_util_misc.h"
 #include "test/fixture.h"
 
 namespace mujoco {
 namespace {
 
 using ::testing::HasSubstr;
+using ::testing::DoubleNear;
 
 TEST_F(MujocoTest, PrintsMemoryWarning) {
   EXPECT_THAT(mju_warningText(mjWARN_CNSTRFULL, pow(2, 10)),
@@ -42,5 +44,31 @@ TEST_F(MujocoTest, PrintsMemoryWarning) {
               HasSubstr("1073741825 bytes"));
 }
 
+TEST_F(MujocoTest, Sigmoid) {
+  // function values
+  EXPECT_EQ(mju_sigmoid(-1),  0);
+  EXPECT_EQ(mju_sigmoid(0),   0);
+  EXPECT_EQ(mju_sigmoid(0.5), 0.5);
+  EXPECT_EQ(mju_sigmoid(1),   1);
+  EXPECT_EQ(mju_sigmoid(2),   1);
+
+  // epsilon for finite-differencing
+  const mjtNum dx = 1e-7;
+
+  // derivative at 0
+  mjtNum dy_dx_0 = (mju_sigmoid(0 + dx) - mju_sigmoid(0)) / dx;
+  EXPECT_THAT(dy_dx_0, DoubleNear(0, dx));
+
+  // derivative at 1
+  mjtNum dy_dx_1 = (mju_sigmoid(1) - mju_sigmoid(1 - dx)) / dx;
+  EXPECT_THAT(dy_dx_1, DoubleNear(0, dx));
+
+  // derivative at 0.5
+  const mjtNum x = 0.5;
+  mjtNum dy_dx_0p5 = (mju_sigmoid(x + dx) - mju_sigmoid(x - dx)) / (2*dx);
+  mjtNum expected = 30*x*x*x*x - 60*x*x*x + 30*x*x;
+  EXPECT_THAT(dy_dx_0p5, DoubleNear(expected, dx));
+}
+
 }  // namespace
 }  // namespace mujoco