Fix #446: reject a substitution that drops the symbol's assumptions

include/numsim_cas/core/substitute.h

@@ -8,8 +8,15 @@
 
 namespace numsim::cas::detail {
 
+// Domains that attach assumptions to symbols overload this for their holder
+// type (found by ADL); this fallback covers the rest.
+template <class TargetBase>
+inline void validate_substitution(expression_holder<TargetBase> const &,
+                                  expression_holder<TargetBase> const &) {}
+
 struct substitute_fn {
-  // explicit typed call
+  // explicit typed call — skips the assumption check, so the substitution
+  // visitors use it to recurse into children
   template <class ExprBase, class TargetBase>
   constexpr auto operator()(std::type_identity<ExprBase>,
                             std::type_identity<TargetBase>,
@@ -34,17 +41,17 @@ struct substitute_fn {
                       std::type_identity<TargetBase>{}, expr, old_val, new_val);
   }
 
-  // ergonomic call: substitute(expr, old, new)
+  // ergonomic call: substitute(expr, old, new) — the public entry, which
+  // rejects a replacement that does not carry the assumptions asserted on
+  // the symbol it replaces
   template <class ExprBase, class TargetBase>
-  constexpr auto operator()(expression_holder<ExprBase> const &expr,
-                            expression_holder<TargetBase> const &old_val,
-                            expression_holder<TargetBase> const &new_val) const
-      noexcept(noexcept((*this)(std::type_identity<ExprBase>{},
-                                std::type_identity<TargetBase>{}, expr, old_val,
-                                new_val)))
-          -> decltype((*this)(std::type_identity<ExprBase>{},
-                              std::type_identity<TargetBase>{}, expr, old_val,
-                              new_val)) {
+  auto operator()(expression_holder<ExprBase> const &expr,
+                  expression_holder<TargetBase> const &old_val,
+                  expression_holder<TargetBase> const &new_val) const
+      -> decltype((*this)(std::type_identity<ExprBase>{},
+                          std::type_identity<TargetBase>{}, expr, old_val,
+                          new_val)) {
+    validate_substitution(old_val, new_val);
     return (*this)(std::type_identity<ExprBase>{},
                    std::type_identity<TargetBase>{}, expr, old_val, new_val);
   }

include/numsim_cas/scalar/visitors/scalar_substitution.h

@@ -3,6 +3,7 @@
 
 #include <numsim_cas/core/substitute.h>
 #include <numsim_cas/scalar/visitors/scalar_rebuild_visitor.h>
+#include <numsim_cas/substitution_guard.h>
 
 namespace numsim::cas {
 

include/numsim_cas/substitution_guard.h

@@ -0,0 +1,153 @@
+#ifndef SUBSTITUTION_GUARD_H
+#define SUBSTITUTION_GUARD_H
+
+#include <string>
+
+#include <numsim_cas/core/cas_error.h>
+#include <numsim_cas/core/substitute.h>
+#include <numsim_cas/scalar/scalar_assume.h>
+#include <numsim_cas/tensor/tensor_assume.h>
+#include <numsim_cas/tensor_to_scalar/tensor_to_scalar_scalar_wrapper.h>
+
+namespace numsim::cas {
+
+// A fold may have consumed the assumptions asserted on a symbol and erased
+// the operation that used them (skew(S) -> 0 for symmetric S), so a
+// replacement that does not provably carry the same assumptions would leave
+// an unjustified result behind. Unprovable counts as not carried.
+
+namespace detail {
+
+[[noreturn]] inline void reject_substitution(std::string const &fact) {
+  throw invalid_expression_error(
+      "substitute: replacement does not carry the assumption '" + fact +
+      "' asserted on the symbol it replaces");
+}
+
+inline char const *numeric_assumption_name(numeric_assumption const &a) {
+  static_assert(std::variant_size_v<numeric_assumption> == 13);
+  static constexpr char const *names[] = {
+      "positive", "negative", "nonzero", "nonnegative", "nonpositive",
+      "integer",  "even",     "odd",     "rational",    "irrational",
+      "real",     "complex",  "prime"};
+  return names[a.index()];
+}
+
+inline char const *
+tensor_algebra_assumption_name(tensor_algebra_assumption const &a) {
+  static_assert(std::variant_size_v<tensor_algebra_assumption> == 5);
+  static constexpr char const *names[] = {
+      "orthogonal", "positive_definite", "positive_semidefinite",
+      "proper_rotation", "improper_rotation"};
+  return names[a.index()];
+}
+
+} // namespace detail
+
+inline void
+validate_substitution(expression_holder<scalar_expression> const &old_val,
+                      expression_holder<scalar_expression> const &new_val) {
+  if (!old_val.is_valid() || !new_val.is_valid() || !old_val.get().is_symbol())
+    return;
+  auto const &facts = old_val.get().assumptions().data();
+  if (facts.empty())
+    return;
+  infer_assumptions(new_val);
+  for (auto const &a : facts) {
+    if (!new_val.get().assumptions().contains(a))
+      detail::reject_substitution(detail::numeric_assumption_name(a));
+  }
+}
+
+inline void
+validate_substitution(expression_holder<tensor_expression> const &old_val,
+                      expression_holder<tensor_expression> const &new_val) {
+  if (!old_val.is_valid() || !new_val.is_valid() || !old_val.get().is_symbol())
+    return;
+
+  if (auto const &sp = old_val.get().space(); sp) {
+    bool perm_ok = true;
+    if (std::holds_alternative<Symmetric>(sp->perm))
+      perm_ok = is_symmetric(new_val);
+    else if (std::holds_alternative<Skew>(sp->perm))
+      perm_ok = is_skew(new_val);
+    else if (std::holds_alternative<Minor>(sp->perm))
+      perm_ok = is_minor(new_val);
+    else if (std::holds_alternative<Major>(sp->perm))
+      perm_ok = is_major(new_val);
+    else if (std::holds_alternative<MinorMajor>(sp->perm))
+      perm_ok = is_minor_major(new_val);
+    else if (!std::holds_alternative<General>(sp->perm)) {
+      // payload-carrying tags (Young) have no predicate: require the same tag
+      auto const &np = new_val.get().space();
+      perm_ok = np && np->perm == sp->perm;
+    }
+    if (!perm_ok)
+      detail::reject_substitution("permutation symmetry");
+
+    if (!std::holds_alternative<AnyTraceTag>(sp->trace)) {
+      bool trace_ok = false;
+      if (std::holds_alternative<VolumetricTag>(sp->trace))
+        trace_ok = is_volumetric(new_val);
+      else if (std::holds_alternative<DeviatoricTag>(sp->trace))
+        trace_ok = is_deviatoric(new_val);
+      else {
+        auto const &np = new_val.get().space();
+        trace_ok = np && np->trace == sp->trace;
+      }
+      if (!trace_ok)
+        detail::reject_substitution("trace constraint");
+    }
+  }
+
+  for (auto const &a : old_val.get().tensor_algebra_assumptions().data()) {
+    bool ok = false;
+    if (std::holds_alternative<orthogonal>(a))
+      ok = is_orthogonal(new_val);
+    else if (std::holds_alternative<proper_rotation>(a))
+      ok = is_proper_rotation(new_val);
+    else if (std::holds_alternative<improper_rotation>(a))
+      ok = is_improper_rotation(new_val);
+    else if (std::holds_alternative<positive_definite>(a))
+      ok = is_positive_definite(new_val);
+    else if (std::holds_alternative<positive_semidefinite>(a))
+      ok = is_positive_semidefinite(new_val);
+    if (!ok)
+      detail::reject_substitution(detail::tensor_algebra_assumption_name(a));
+  }
+}
+
+inline void validate_substitution(
+    expression_holder<tensor_to_scalar_expression> const &old_val,
+    expression_holder<tensor_to_scalar_expression> const &new_val) {
+  if (!old_val.is_valid() || !new_val.is_valid() || !old_val.get().is_symbol())
+    return;
+  // A t2s symbol is a wrapped scalar, and the facts sit on the wrapped
+  // expression; a non-wrapper replacement must carry them itself.
+  auto unwrap = [](expression_holder<tensor_to_scalar_expression> const &h) {
+    expression_holder<scalar_expression> inner;
+    if (is_same<tensor_to_scalar_scalar_wrapper>(h))
+      inner = h.get<tensor_to_scalar_scalar_wrapper>().expr();
+    return inner;
+  };
+  auto old_inner = unwrap(old_val);
+  if (!old_inner.is_valid()) {
+    for (auto const &a : old_val.get().assumptions().data()) {
+      if (!new_val.get().assumptions().contains(a))
+        detail::reject_substitution(detail::numeric_assumption_name(a));
+    }
+    return;
+  }
+  if (auto new_inner = unwrap(new_val); new_inner.is_valid()) {
+    validate_substitution(old_inner, new_inner);
+    return;
+  }
+  for (auto const &a : old_inner.get().assumptions().data()) {
+    if (!new_val.get().assumptions().contains(a))
+      detail::reject_substitution(detail::numeric_assumption_name(a));
+  }
+}
+
+} // namespace numsim::cas
+
+#endif // SUBSTITUTION_GUARD_H

include/numsim_cas/tensor/visitors/tensor_substitution.h

@@ -3,6 +3,7 @@
 
 #include <numsim_cas/core/substitute.h>
 #include <numsim_cas/scalar/visitors/scalar_substitution.h>
+#include <numsim_cas/substitution_guard.h>
 #include <numsim_cas/tensor/visitors/tensor_rebuild_visitor.h>
 
 namespace numsim::cas {
@@ -35,14 +36,16 @@ class tensor_substitution final : public tensor_rebuild_visitor {
 
   scalar_holder_t apply_scalar(scalar_holder_t const &expr) override {
     if constexpr (std::is_same_v<TargetBase, scalar_expression>) {
-      return substitute(expr, m_old, m_new);
+      return substitute(std::type_identity<scalar_expression>{},
+                        std::type_identity<TargetBase>{}, expr, m_old, m_new);
     } else {
       return expr;
     }
   }
 
   t2s_holder_t apply_t2s(t2s_holder_t const &expr) override {
-    return substitute(expr, m_old, m_new);
+    return substitute(std::type_identity<tensor_to_scalar_expression>{},
+                      std::type_identity<TargetBase>{}, expr, m_old, m_new);
   }
 
 private:

include/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_substitution.h

@@ -3,6 +3,7 @@
 
 #include <numsim_cas/core/substitute.h>
 #include <numsim_cas/scalar/visitors/scalar_substitution.h>
+#include <numsim_cas/substitution_guard.h>
 #include <numsim_cas/tensor/visitors/tensor_substitution.h>
 #include <numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_rebuild_visitor.h>
 
@@ -38,7 +39,8 @@ class tensor_to_scalar_substitution final
 
   scalar_holder_t apply_scalar(scalar_holder_t const &expr) override {
     if constexpr (std::is_same_v<TargetBase, scalar_expression>) {
-      return substitute(expr, m_old, m_new);
+      return substitute(std::type_identity<scalar_expression>{},
+                        std::type_identity<TargetBase>{}, expr, m_old, m_new);
     } else {
       return expr;
     }
@@ -47,7 +49,11 @@ class tensor_to_scalar_substitution final
   // Unconditional: a tensor subtree can carry scalar and t2s children, so a
   // needle of any domain may hide inside it.
   tensor_holder_t apply_tensor(tensor_holder_t const &expr) override {
-    return substitute(expr, m_old, m_new);
+    // Descend unconditionally: a t2s needle can sit inside a tensor child.
+    // The typed call skips re-validating what the public entry point
+    // already checked.
+    return substitute(std::type_identity<tensor_expression>{},
+                      std::type_identity<TargetBase>{}, expr, m_old, m_new);
   }
 
 private:

src/numsim_cas/tensor/tensor_solve.cpp

@@ -65,7 +65,11 @@ std::vector<tensor_solver::expr_holder_t> tensor_solver::solve() const {
   // 5. Get constant term by substituting X = 0
   auto zero_tensor =
       make_expression<tensor_zero>(m_x.get().dim(), m_x.get().rank());
-  auto b = substitute(m_expr, m_x, zero_tensor);
+  // X = 0 reads off the constant term; it is a probe, not a claim that zero
+  // satisfies X's assumptions, so it uses the unchecked typed call.
+  auto b = substitute(std::type_identity<tensor_expression>{},
+                      std::type_identity<tensor_expression>{}, m_expr, m_x,
+                      zero_tensor);
 
   // 6. Build solution
   auto neg_b = -b;

tests/CoreBugFixTest.h

@@ -6,6 +6,7 @@
 #include <cmath>
 #include <numsim_cas/core/substitute.h>
 #include <numsim_cas/tensor/visitors/tensor_substitution.h>
+#include <numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_substitution.h>
 
 namespace numsim::cas {
 
@@ -2153,6 +2154,90 @@ TEST(SubstitutionSpace, OperatorDerivedTagSurvives) {
   EXPECT_EQ(to_string(sym(h)), "0{2}");
 }
 
+// A fold can consume a symbol's assumption and erase the operation, so a
+// replacement that does not provably carry the same assumption is rejected.
+TEST(SubstitutionAssumptions, TensorReplacementMustCarryTheAssumption) {
+  auto [S, G] =
+      make_tensor_variable(std::tuple{"S", std::size_t{3}, std::size_t{2}},
+                           std::tuple{"G", std::size_t{3}, std::size_t{2}});
+  S.assumption(Symmetric{});
+  ASSERT_EQ(to_string(skew(S)), "0{2}");
+  EXPECT_THROW((void)substitute(skew(S), S, G), invalid_expression_error);
+  EXPECT_THROW((void)substitute(sym(S), S, G), invalid_expression_error);
+
+  // a provably symmetric replacement keeps the fold justified
+  auto [X] =
+      make_tensor_variable(std::tuple{"X", std::size_t{3}, std::size_t{2}});
+  EXPECT_EQ(to_string(substitute(skew(S), S, sym(X))), "0{2}");
+  auto [H] =
+      make_tensor_variable(std::tuple{"H", std::size_t{3}, std::size_t{2}});
+  H.assumption(Symmetric{});
+  EXPECT_EQ(to_string(substitute(skew(S), S, H)), "0{2}");
+
+  // an unannotated symbol is unprovable, not merely unequal: still rejected
+  EXPECT_THROW((void)substitute(S + G, S, G), invalid_expression_error);
+  // a skew replacement violates the assumption outright
+  EXPECT_THROW((void)substitute(skew(S), S, skew(X)), invalid_expression_error);
+  // substituting into an unannotated symbol is unaffected
+  EXPECT_NO_THROW((void)substitute(skew(G), G, X));
+}
+
+TEST(SubstitutionAssumptions, AlgebraAssumptionsAreChecked) {
+  auto [P, Q] =
+      make_tensor_variable(std::tuple{"P", std::size_t{3}, std::size_t{2}},
+                           std::tuple{"Q", std::size_t{3}, std::size_t{2}});
+  P.assumption(positive_definite{});
+  EXPECT_THROW((void)substitute(det(P), P, Q), invalid_expression_error);
+  Q.assumption(positive_definite{});
+  EXPECT_NO_THROW((void)substitute(det(P), P, Q));
+
+  auto [R, T] =
+      make_tensor_variable(std::tuple{"R", std::size_t{3}, std::size_t{2}},
+                           std::tuple{"T", std::size_t{3}, std::size_t{2}});
+  R.assumption(orthogonal{});
+  EXPECT_THROW((void)substitute(inv(R), R, T), invalid_expression_error);
+  T.assumption(orthogonal{});
+  EXPECT_NO_THROW((void)substitute(inv(R), R, T));
+}
+
+TEST(SubstitutionAssumptions, ScalarReplacementMustCarryTheAssumption) {
+  auto [p, q] = make_scalar_variable("p", "q");
+  p.assumption(positive{});
+  ASSERT_EQ(to_string(abs(p)), "p");
+  EXPECT_THROW((void)substitute(abs(p), p, q), invalid_expression_error);
+  q.assumption(positive{});
+  EXPECT_EQ(to_string(substitute(abs(p), p, q)), "q");
+
+  auto [n] = make_scalar_variable("n");
+  n.assumption(negative{});
+  EXPECT_THROW((void)substitute(abs(p), p, n), invalid_expression_error);
+}
+
+TEST(SubstitutionAssumptions, T2sWrapperForwardsToTheWrappedScalar) {
+  auto [p, q] = make_scalar_variable("p", "q");
+  p.assumption(positive{});
+  auto wp = make_expression<tensor_to_scalar_scalar_wrapper>(p);
+  auto wq = make_expression<tensor_to_scalar_scalar_wrapper>(q);
+  auto [A] =
+      make_tensor_variable(std::tuple{"A", std::size_t{3}, std::size_t{2}});
+  EXPECT_THROW((void)substitute(wp * trace(A), wp, wq),
+               invalid_expression_error);
+  q.assumption(positive{});
+  EXPECT_NO_THROW((void)substitute(wp * trace(A), wp, wq));
+}
+
+// The solver substitutes X = 0 to read off the constant term; that is an
+// internal probe, not a claim that the replacement satisfies X's assumptions.
+TEST(SubstitutionAssumptions, SolverProbeIsNotValidated) {
+  auto [X, B] =
+      make_tensor_variable(std::tuple{"X", std::size_t{3}, std::size_t{2}},
+                           std::tuple{"B", std::size_t{3}, std::size_t{2}});
+  X.assumption(positive_definite{});
+  EXPECT_NO_THROW((void)solve(X - B, X));
+  auto [c] = make_scalar_variable("c");
+  EXPECT_NO_THROW((void)solve(c * X - B, X));
+}
+
 // Substituting a dim-2 argument into a dim-3 projector contraction is a
 // shape error: rejected when the rebuilt node is constructed.
 TEST(RoundTwoReview, DimChangingSubstitutionThrows) {