diff --git a/include/numsim_cas/core/substitute.h b/include/numsim_cas/core/substitute.h index b66bd833..65476da2 100644 --- a/include/numsim_cas/core/substitute.h +++ b/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 +inline void validate_substitution(expression_holder const &, + expression_holder 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 constexpr auto operator()(std::type_identity, std::type_identity, @@ -34,17 +41,17 @@ struct substitute_fn { std::type_identity{}, 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 - constexpr auto operator()(expression_holder const &expr, - expression_holder const &old_val, - expression_holder const &new_val) const - noexcept(noexcept((*this)(std::type_identity{}, - std::type_identity{}, expr, old_val, - new_val))) - -> decltype((*this)(std::type_identity{}, - std::type_identity{}, expr, old_val, - new_val)) { + auto operator()(expression_holder const &expr, + expression_holder const &old_val, + expression_holder const &new_val) const + -> decltype((*this)(std::type_identity{}, + std::type_identity{}, expr, old_val, + new_val)) { + validate_substitution(old_val, new_val); return (*this)(std::type_identity{}, std::type_identity{}, expr, old_val, new_val); } diff --git a/include/numsim_cas/scalar/visitors/scalar_substitution.h b/include/numsim_cas/scalar/visitors/scalar_substitution.h index a4c11af0..30266887 100644 --- a/include/numsim_cas/scalar/visitors/scalar_substitution.h +++ b/include/numsim_cas/scalar/visitors/scalar_substitution.h @@ -3,6 +3,7 @@ #include #include +#include namespace numsim::cas { diff --git a/include/numsim_cas/substitution_guard.h b/include/numsim_cas/substitution_guard.h new file mode 100644 index 00000000..986ed1f7 --- /dev/null +++ b/include/numsim_cas/substitution_guard.h @@ -0,0 +1,153 @@ +#ifndef SUBSTITUTION_GUARD_H +#define SUBSTITUTION_GUARD_H + +#include + +#include +#include +#include +#include +#include + +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 == 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 == 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 const &old_val, + expression_holder 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 const &old_val, + expression_holder 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(sp->perm)) + perm_ok = is_symmetric(new_val); + else if (std::holds_alternative(sp->perm)) + perm_ok = is_skew(new_val); + else if (std::holds_alternative(sp->perm)) + perm_ok = is_minor(new_val); + else if (std::holds_alternative(sp->perm)) + perm_ok = is_major(new_val); + else if (std::holds_alternative(sp->perm)) + perm_ok = is_minor_major(new_val); + else if (!std::holds_alternative(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(sp->trace)) { + bool trace_ok = false; + if (std::holds_alternative(sp->trace)) + trace_ok = is_volumetric(new_val); + else if (std::holds_alternative(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(a)) + ok = is_orthogonal(new_val); + else if (std::holds_alternative(a)) + ok = is_proper_rotation(new_val); + else if (std::holds_alternative(a)) + ok = is_improper_rotation(new_val); + else if (std::holds_alternative(a)) + ok = is_positive_definite(new_val); + else if (std::holds_alternative(a)) + ok = is_positive_semidefinite(new_val); + if (!ok) + detail::reject_substitution(detail::tensor_algebra_assumption_name(a)); + } +} + +inline void validate_substitution( + expression_holder const &old_val, + expression_holder 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 const &h) { + expression_holder inner; + if (is_same(h)) + inner = h.get().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 diff --git a/include/numsim_cas/tensor/visitors/tensor_substitution.h b/include/numsim_cas/tensor/visitors/tensor_substitution.h index 9cd71e33..932c15d2 100644 --- a/include/numsim_cas/tensor/visitors/tensor_substitution.h +++ b/include/numsim_cas/tensor/visitors/tensor_substitution.h @@ -3,6 +3,7 @@ #include #include +#include #include 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) { - return substitute(expr, m_old, m_new); + return substitute(std::type_identity{}, + std::type_identity{}, 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{}, + std::type_identity{}, expr, m_old, m_new); } private: diff --git a/include/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_substitution.h b/include/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_substitution.h index b79deb73..09969f43 100644 --- a/include/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_substitution.h +++ b/include/numsim_cas/tensor_to_scalar/visitors/tensor_to_scalar_substitution.h @@ -3,6 +3,7 @@ #include #include +#include #include #include @@ -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) { - return substitute(expr, m_old, m_new); + return substitute(std::type_identity{}, + std::type_identity{}, 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{}, + std::type_identity{}, expr, m_old, m_new); } private: diff --git a/src/numsim_cas/tensor/tensor_solve.cpp b/src/numsim_cas/tensor/tensor_solve.cpp index df482e22..0a2c9600 100644 --- a/src/numsim_cas/tensor/tensor_solve.cpp +++ b/src/numsim_cas/tensor/tensor_solve.cpp @@ -65,7 +65,11 @@ std::vector tensor_solver::solve() const { // 5. Get constant term by substituting X = 0 auto zero_tensor = make_expression(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{}, + std::type_identity{}, m_expr, m_x, + zero_tensor); // 6. Build solution auto neg_b = -b; diff --git a/tests/CoreBugFixTest.h b/tests/CoreBugFixTest.h index 2b96298b..f36b2a38 100644 --- a/tests/CoreBugFixTest.h +++ b/tests/CoreBugFixTest.h @@ -6,6 +6,7 @@ #include #include #include +#include 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(p); + auto wq = make_expression(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) {