feat: interp1d linear interpolant (#141)

* feat: interp1d linear interpolant

* fix ci

* another fix to ci

* fix: tests

* clippy

Author: martinjrobins

.github/workflows/ci.yml

@@ -180,6 +180,26 @@ jobs:
     - name: Setup tmate session
       uses: mxschmitt/action-tmate@v3
       if: ${{ github.event_name == 'workflow_dispatch' && inputs.debug_enabled }}
+    - name: Workaround missing libcurses.tbd on macOS
+      if : matrix.os == 'macos-latest' && matrix.llvm != ''
+      run: |
+        SDK_PATH=$(xcrun --sdk macosx --show-sdk-path)
+        if [ ! -f "$SDK_PATH/usr/lib/libcurses.tbd" ] && [ -f "$SDK_PATH/usr/lib/libncurses.tbd" ]; then
+          # 1. Patch hardcoded libcurses.tbd paths in LLVM cmake config files
+          LLVM_PREFIX="${LLVM_PATH:-$GITHUB_WORKSPACE/llvm}"
+          LLVM_CMAKE_DIR="$LLVM_PREFIX/lib/cmake/llvm"
+          if [ -d "$LLVM_CMAKE_DIR" ]; then
+            grep -rl 'libcurses\.tbd' "$LLVM_CMAKE_DIR" 2>/dev/null | while IFS= read -r f; do
+              sed -i '' 's|/usr/lib/libcurses\.tbd|/usr/lib/libncurses.tbd|g' "$f"
+            done
+          fi
+          # 2. Provide a symlink via CMAKE_PREFIX_PATH so find_library(curses)
+          #    can resolve to libncurses when the path is resolved dynamically
+          TMP_CURSES_DIR=$(mktemp -d)
+          mkdir -p "$TMP_CURSES_DIR/lib"
+          ln -s "$SDK_PATH/usr/lib/libncurses.tbd" "$TMP_CURSES_DIR/lib/libcurses.tbd"
+          echo "CMAKE_PREFIX_PATH=$TMP_CURSES_DIR" >> "$GITHUB_ENV"
+        fi
     - name: Run tests (LLVM)
       if : matrix.tests && matrix.llvm != ''
       run: cargo test --verbose --features llvm${{ matrix.llvm[1] }} --features "${{ matrix.features }}"

book/src/functions.md

@@ -20,6 +20,7 @@ The DiffSL supports the following mathematical functions that can be used in an
 * `abs(x)` - absolute value of x
 * `sigmoid(x)` - sigmoid function of x
 * `heaviside(x)` - Heaviside step function of x
+* `interp1d(x_i, y_i, q)` - 1D piecewise linear interpolation
 
 You can use these functions as part of an expression in the DSL. For example, to define a variable `a` that is the sine of another variable `b`, you can write:
 
@@ -28,6 +29,33 @@ b { 1.0 }
 a { sin(b) }
 ```
 
+### `interp1d(x_i, y_i, q)`
+
+The `interp1d` function performs 1D piecewise linear interpolation. `x_i` and `y_i`
+must be compile-time constant 1D dense tensors of equal length. The query point
+`q` can be a scalar or a tensor. Out-of-range queries are clamped to the
+endpoints.
+
+`x_i` must be **strictly increasing**; an error is reported at compile time if
+it is not.
+
+If `x_i` is uniformly spaced (constant increment `dx` between consecutive values),
+the interpolation runs in **O(1)** per query point.  For non-uniform `x_i`, a
+binary search is used, giving **O(log n)** cost per query point.
+
+```diffsl
+xs_i { 0.0, 10.0, 20.0 }
+ys_i { 1.0,  5.0, 15.0 }
+r { interp1d(xs_i, ys_i, 5.0) }     # => 3.0
+
+# batched queries
+q_j { 0.0, 5.0, 10.0, 15.0, 20.0 }
+r_j { interp1d(xs_i, ys_i, q_j) }   # => [1.0, 3.0, 5.0, 10.0, 15.0]
+
+# the first two arguments may be any constant 1D dense expression
+r { interp1d(xs_i * 2.0, ys_i + 1.0, 5.0) }
+```
+
 ## Pre-defined variables
 
 There are two predefined variables in DiffSL. The first is the scalar `t` which is the current time, this allows the equations to be written as functions of time. For example

diffsl/src/ast/mod.rs

@@ -525,7 +525,7 @@ impl<'a> AstKind<'a> {
     }
 }
 
-#[derive(Debug, Copy, Clone)]
+#[derive(Debug, Copy, Clone, Hash, PartialEq, Eq)]
 pub struct StringSpan {
     pub pos_start: usize,
     pub pos_end: usize,
@@ -914,8 +914,14 @@ impl<'a> Ast<'a> {
                 monop.child.collect_indices(indices);
             }
             AstKind::Call(call) => {
-                for c in &call.args {
-                    c.collect_indices(indices);
+                if call.fn_name == "interp1d" {
+                    if let Some(q) = call.args.get(2) {
+                        q.collect_indices(indices);
+                    }
+                } else {
+                    for c in &call.args {
+                        c.collect_indices(indices);
+                    }
                 }
             }
             AstKind::CallArg(arg) => {

diffsl/src/continuous/builder.rs

@@ -534,13 +534,18 @@ impl<'s> ModelInfo<'s> {
             }
             AstKind::Call(call) => {
                 // check name in allowed functions
-                let functions = ["sin", "cos", "tan", "pow", "exp", "log", "sqrt", "abs"];
+                let functions = [
+                    "sin", "cos", "tan", "pow", "exp", "log", "sqrt", "abs", "interp1d",
+                ];
                 if functions.contains(&call.fn_name) {
-                    // built in functions all have 1 arg
-                    // built in functions should have no keyword args
-                    if call.args.len() != 1 {
+                    let expected_nargs = match call.fn_name {
+                        "interp1d" => 3,
+                        "pow" => 2,
+                        _ => 1,
+                    };
+                    if call.args.len() != expected_nargs {
                         self.errors.push(Output::new(
-                            format!("incorrect number of given arguments ({} instead of {}) for function {}", call.args.len(), 1, call.fn_name),
+                            format!("incorrect number of given arguments ({} instead of {}) for function {}", call.args.len(), expected_nargs, call.fn_name),
                             expr.span,
                         ));
                     }

diffsl/src/discretise/discrete_model.rs

@@ -1054,6 +1054,37 @@ impl<'s> DiscreteModel<'s> {
         self.reset.as_ref()
     }
 
+    pub fn all_tensors(&self) -> Vec<&Tensor<'_>> {
+        let mut tensors = Vec::new();
+        tensors.extend(self.constant_defns.iter());
+        tensors.extend(self.input_dep_defns.iter());
+        tensors.extend(self.time_dep_defns.iter());
+        tensors.extend(self.state_dep_defns.iter());
+        tensors.extend(self.state_dep_post_f_defns.iter());
+        tensors.extend(self.dstate_dep_defns.iter());
+        tensors.push(&self.state);
+        if let Some(sd) = &self.state_dot {
+            tensors.push(sd);
+        }
+        if let Some(inp) = &self.input {
+            tensors.push(inp);
+        }
+        if let Some(lhs) = &self.lhs {
+            tensors.push(lhs);
+        }
+        tensors.push(&self.rhs);
+        if let Some(out) = &self.out {
+            tensors.push(out);
+        }
+        if let Some(reset) = &self.reset {
+            tensors.push(reset);
+        }
+        if let Some(stop) = &self.stop {
+            tensors.push(stop);
+        }
+        tensors
+    }
+
     pub fn take_state0_input_deps(&mut self) -> Vec<(usize, usize)> {
         std::mem::take(&mut self.state0_input_deps)
     }
@@ -2026,7 +2057,7 @@ mod tests {
         );
         assert!(c.elmts()[0].has_values());
 
-        let layout = DataLayout::new(&model);
+        let layout = DataLayout::new(&model).unwrap();
         assert_eq!(layout.get_tensor_constants("C").unwrap(), &[2.0, 5.0]);
     }
 
@@ -2149,4 +2180,84 @@ mod tests {
             );
         }
     }
+
+    #[test]
+    fn test_interp1d_basic() {
+        let text = "
+            xs_i { 0.0, 10.0, 20.0 }
+            ys_i { 1.0, 5.0, 15.0 }
+            u_i { z = 1 }
+            F_i { interp1d(xs_i, ys_i, 5.0) }
+        ";
+        let model = parse_ds_string(text).unwrap();
+        let discrete = DiscreteModel::build("$name", &model).unwrap();
+        assert_eq!(discrete.constant_defns().len(), 2); // xs and ys are constants
+        assert_eq!(discrete.rhs.elmts().len(), 1);
+    }
+
+    #[test]
+    fn test_interp1d_batched() {
+        let text = "
+            xs_i { 0.0, 10.0, 20.0 }
+            ys_i { 1.0, 5.0, 15.0 }
+            q_j { 0.0, 15.0 }
+            u_j { z = 1, w = 1 }
+            F_j { interp1d(xs_i, ys_i, q_j) }
+        ";
+        let model = parse_ds_string(text).unwrap();
+        let discrete = DiscreteModel::build("$name", &model).unwrap();
+        assert_eq!(discrete.rhs.elmts().len(), 1);
+        assert_eq!(discrete.state.elmts().len(), 2);
+    }
+
+    #[test]
+    fn test_interp1d_sortedness_error() {
+        let text = "
+            xs_i { 10.0, 0.0, 20.0 }
+            ys_i { 1.0, 5.0, 15.0 }
+            u_i { z = 1 }
+            F_i { interp1d(xs_i, ys_i, 5.0) }
+        ";
+        let model = parse_ds_string(text).unwrap();
+        let discrete = DiscreteModel::build("$name", &model).unwrap();
+        let err = DataLayout::new(&discrete).unwrap_err();
+        assert!(
+            err.has_error_contains("strictly increasing"),
+            "expected sortedness error, got: {}",
+            err
+        );
+    }
+
+    #[test]
+    fn test_interp1d_nonconstant_x() {
+        let text = "
+            in { a = 1 }
+            u_i { z = 1 }
+            F_i { interp1d(u_i, u_i, 5.0) }
+        ";
+        let model = parse_ds_string(text).unwrap();
+        let err = DiscreteModel::build("$name", &model).unwrap_err();
+        assert!(
+            err.has_error_contains("compile-time constant"),
+            "expected constant error in: {}",
+            err.as_error_message(text)
+        );
+    }
+
+    #[test]
+    fn test_interp1d_shape_mismatch() {
+        let text = "
+            xs_i { 0.0, 10.0, 20.0 }
+            ys_j { 1.0, 5.0 }
+            u_i { z = 1 }
+            F_i { interp1d(xs_i, ys_j, 5.0) }
+        ";
+        let model = parse_ds_string(text).unwrap();
+        let err = DiscreteModel::build("$name", &model).unwrap_err();
+        assert!(
+            err.has_error_contains("cannot find index j"),
+            "expected index mismatch error in: {}",
+            err.as_error_message(text)
+        );
+    }
 }

diffsl/src/discretise/env.rs

@@ -524,6 +524,9 @@ impl Env {
         ast: &Ast,
         indices: &Vec<char>,
     ) -> Option<Layout> {
+        if call.fn_name == "interp1d" {
+            return self.get_layout_interp1d(call, ast, indices);
+        }
         let layouts = call
             .args
             .iter()
@@ -539,6 +542,90 @@ impl Env {
         }
     }
 
+    fn get_layout_interp1d(
+        &mut self,
+        call: &ast::Call,
+        ast: &Ast,
+        indices: &Vec<char>,
+    ) -> Option<Layout> {
+        if call.args.len() != 3 {
+            self.errs.push(ValidationError::new(
+                format!(
+                    "interp1d requires 3 arguments (x_i, y_i, q), got {}",
+                    call.args.len()
+                ),
+                ast.span,
+            ));
+            return None;
+        }
+
+        for (i, arg) in [&call.args[0], &call.args[1]].iter().enumerate() {
+            let deps = arg.get_dependents();
+            for name in &deps {
+                if self.get(name).is_none_or(|v| !v.is_constant()) {
+                    self.errs.push(ValidationError::new(
+                        format!(
+                            "interp1d argument {} depends on '{}' which is not compile-time constant",
+                            i + 1,
+                            name
+                        ),
+                        ast.span,
+                    ));
+                    return None;
+                }
+            }
+        }
+
+        let shared_index = call.args[0].get_indices().into_iter().next().unwrap_or({
+            // No index found — the expression is scalar, which will fail the rank check below.
+            // Use a dummy index so get_layout can give a meaningful error.
+            ' '
+        });
+        let ctx = vec![shared_index];
+        let x_layout = self.get_layout(&call.args[0], &ctx)?;
+        let y_layout = self.get_layout(&call.args[1], &ctx)?;
+
+        if x_layout.rank() != 1 {
+            self.errs.push(ValidationError::new(
+                format!(
+                    "interp1d first argument must be 1D, got rank {}",
+                    x_layout.rank()
+                ),
+                ast.span,
+            ));
+            return None;
+        }
+        if !x_layout.is_dense() {
+            self.errs.push(ValidationError::new(
+                "interp1d first argument must be a dense vector".to_string(),
+                ast.span,
+            ));
+            return None;
+        }
+        if !y_layout.is_dense() {
+            self.errs.push(ValidationError::new(
+                "interp1d second argument must be a dense vector".to_string(),
+                ast.span,
+            ));
+            return None;
+        }
+        if y_layout.shape() != x_layout.shape() {
+            self.errs.push(ValidationError::new(
+                format!(
+                    "interp1d: x and y must have the same shape, got {} and {}",
+                    x_layout.shape(),
+                    y_layout.shape()
+                ),
+                ast.span,
+            ));
+            return None;
+        }
+
+        let mut q_layout = self.get_layout(&call.args[2], indices)?;
+        q_layout.to_dense();
+        Some(q_layout)
+    }
+
     pub fn get_layout(&mut self, ast: &Ast, indices: &Vec<char>) -> Option<Layout> {
         let layout = match &ast.kind {
             AstKind::Assignment(a) => self.get_layout(a.expr.as_ref(), indices),

diffsl/src/discretise/error.rs

@@ -82,3 +82,5 @@ impl fmt::Display for ValidationErrors {
         Ok(())
     }
 }
+
+impl std::error::Error for ValidationErrors {}