System.eval accepts any point precision; mixed time/space precisions throw (closes #377)

core/include/bertini2/system/system.hpp

@@ -333,9 +333,9 @@ namespace bertini {
 		/**
 		 Evaluate the system, provided a path variable is defined for the system, in place.
 
-		 \throws std::runtime_error, if a path variable is NOT defined, and you passed it a value.  Also throws if the number of variables doesn't match.
+		 \throws std::runtime_error, if a path variable is NOT defined, and you passed it a value.  Also throws if the number of variables doesn't match, or if (for multiprecision types) the precisions of ``variable_values`` and ``path_variable_value`` differ -- both are supplied in the same call, so align them before evaluating.
 		 \tparam T the number-type for return.  Probably complex_dbl=std::complex<double>, or complex_mp=bertini::complex_mp.
-		 
+
 		 \param function_values The vector to write the function values into.
 		 \param variable_values The values of the variables, for the evaluation.
 		 \param path_variable_value The current value of the path variable.
@@ -352,6 +352,24 @@ namespace bertini {
 			if (!have_path_variable_)
 				throw std::runtime_error("trying to use a time value for evaluation of system, but no path variable defined.");
 
+			#ifndef BERTINI_DISABLE_PRECISION_CHECKS
+				// The variables and the time are BOTH caller-supplied in this one call, so a
+				// precision mismatch between them is caller incoherence, not internal-state
+				// drift -- silently promoting either one would guess intent and mask exactly
+				// the ambient-precision-drift bugs this check exists to expose.  (Unlike the
+				// system's own precision, which follows the input point; see SetVariables.)
+				if constexpr (!std::is_same<T,complex_dbl>::value) {
+					if (variable_values.size() > 0
+					    && Precision(variable_values) != Precision(path_variable_value))
+						throw std::runtime_error(
+							"precision of the variable values ("
+							+ std::to_string(Precision(variable_values))
+							+ ") differs from the precision of the path-variable value ("
+							+ std::to_string(Precision(path_variable_value))
+							+ ") in the same evaluation call; align them before evaluating");
+				}
+			#endif
+
 			SetVariables(variable_values.eval());
 			SetPathVariable(path_variable_value);
 
@@ -514,8 +532,8 @@ namespace bertini {
 
 		/**
 		 Evaluate the Jacobian of the system, provided a path variable is defined for the system, in place.
-		 
-		 \throws std::runtime_error, if a path variable is NOT defined, and you passed it a value.  Also throws if the number of variables doesn't match.
+
+		 \throws std::runtime_error, if a path variable is NOT defined, and you passed it a value.  Also throws if the number of variables doesn't match, or if (for multiprecision types) the precisions of ``variable_values`` and ``path_variable_value`` differ -- both are supplied in the same call, so align them before evaluating.
 
 		 \tparam T the number-type for return.  Probably complex_dbl=std::complex<double>, or complex_mp=bertini::complex_mp.
 
@@ -530,10 +548,26 @@ namespace bertini {
 
 			if (variable_values.size()!=static_cast<Eigen::Index>(NumVariables()))
 				throw std::runtime_error("trying to evaluate jacobian, but number of variables doesn't match.");
-			
+
 			if (!HavePathVariable())
 				throw std::runtime_error("trying to use a time value for computation of jacobian, but no path variable defined.");
-
+
+			#ifndef BERTINI_DISABLE_PRECISION_CHECKS
+				// same coherence rule as EvalInPlace: both arguments came from this one
+				// call, so a precision mismatch between them is caller incoherence, not
+				// internal-state drift (which SetVariables absorbs by following the input)
+				if constexpr (!std::is_same<T,complex_dbl>::value) {
+					if (variable_values.size() > 0
+					    && Precision(variable_values) != Precision(path_variable_value))
+						throw std::runtime_error(
+							"precision of the variable values ("
+							+ std::to_string(Precision(variable_values))
+							+ ") differs from the precision of the path-variable value ("
+							+ std::to_string(Precision(path_variable_value))
+							+ ") in the same jacobian call; align them before evaluating");
+				}
+			#endif
+
 			SetVariables(variable_values.eval());
 			SetPathVariable(path_variable_value);
 			JacobianInPlace(J);
@@ -877,13 +911,18 @@ namespace bertini {
 		 \tparam T the number-type for return.  Probably complex_dbl=std::complex<double>, or complex_mp=bertini::complex_mp.
 		 \throws std::runtime_error if the number of variables doesn't match.
 
-		 The ordering of the variables matters.  
+		 The ordering of the variables matters.
 
 		 * The AffHomUng ordering is 1) variable groups, with homogenizing variable first. 2) homogeneous variable groups. 3) ungrouped variables.
 		 * The FIFO ordering uses the order in which the variable groups were added.
 
 		 The path variable is not considered a variable for this operation.  It is set separately.
-
+
+		 A multiprecision point of any precision is accepted: the input defines the
+		 working precision, and the system's internal precision is changed to match.
+		 Callers need not pre-align a point to the system.  Under
+		 BERTINI_DISABLE_PRECISION_CHECKS no alignment (nor any check) is performed.
+
 		 \param new_values The new updated values for the variables.
 
 		 \see SetPathVariable
@@ -903,8 +942,10 @@ namespace bertini {
 				if (new_values.size() > 0)
 				{
 					if constexpr (!std::is_same<T,complex_dbl>::value) {
-						if (Precision(new_values) != this->precision())
-							throw std::runtime_error("precision of input point in SetVariables (" + std::to_string(Precision(new_values)) + ") must match the precision of the system (" + std::to_string(this->precision()) + ").");
+						const auto point_prec = Precision(new_values);
+						if (point_prec != this->precision())
+							this->precision(point_prec);   // the input defines the
+							                               // working precision (#377)
 					}
 				}
 			#endif

core/test/classes/system_test.cpp

@@ -497,6 +497,75 @@ BOOST_AUTO_TEST_CASE(system_evaluate_mpfr)
 }
 
 
+/**
+\class bertini::System
+\test \b system_evaluate_mpfr_any_point_precision Evaluate at multiprecision points whose precision differs from the system's -- the input defines the working precision, and the system follows it in both directions.
+*/
+BOOST_AUTO_TEST_CASE(system_evaluate_mpfr_any_point_precision)
+{
+	bertini::DefaultPrecision(30);
+
+	std::string str = "function f; variable_group x1, x2; y = x1*x2; f = y*y;";
+	bertini::System sys;
+	[[maybe_unused]] bool s = bertini::parsing::classic::parse(str.begin(), str.end(), sys);
+
+	{ // higher-precision point: the system follows the input up
+		bertini::DefaultPrecision(50);
+		Vec<mpfr> values(2);
+		values << mpfr(2), mpfr(3);
+		auto v = sys.Eval(values);
+		BOOST_CHECK_EQUAL(sys.precision(), 50u);
+		BOOST_CHECK_EQUAL(Precision(v), 50u);
+		BOOST_CHECK_EQUAL(v(0), mpfr(36));
+
+		auto J = sys.Jacobian(values);
+		BOOST_CHECK_EQUAL(J(0,0), mpfr(36));   // 2*x1*x2^2 at (2,3)
+		BOOST_CHECK_EQUAL(J(0,1), mpfr(24));   // 2*x1^2*x2 at (2,3)
+	}
+
+	{ // lower-precision point: the system follows the input down, too
+		bertini::DefaultPrecision(20);
+		Vec<mpfr> values(2);
+		values << mpfr(2), mpfr(3);
+		auto v = sys.Eval(values);
+		BOOST_CHECK_EQUAL(sys.precision(), 20u);
+		BOOST_CHECK_EQUAL(v(0), mpfr(36));
+	}
+}
+
+
+/**
+\class bertini::System
+\test \b system_evaluate_mixed_time_precision_throws The variables and the time value are both supplied in one evaluation call, so a precision mismatch between THEM is caller incoherence and throws; aligned inputs evaluate at their common precision.
+*/
+BOOST_AUTO_TEST_CASE(system_evaluate_mixed_time_precision_throws)
+{
+	bertini::DefaultPrecision(30);
+
+	Var x = Variable::Make("x");
+	Var t = Variable::Make("t");
+	bertini::System S;
+	S.AddUngroupedVariable(x);
+	S.AddPathVariable(t);
+	S.AddFunction((1-t)*x);
+
+	Vec<mpfr> v(1);
+	v << mpfr(2);                                  // precision 30
+
+	bertini::DefaultPrecision(50);
+	mpfr time_hi(1);                               // precision 50: incoherent with v
+	BOOST_CHECK_THROW(S.Eval(v, time_hi), std::runtime_error);
+	Mat<mpfr> J(S.NumTotalFunctions(), S.NumVariables());
+	BOOST_CHECK_THROW(S.JacobianInPlace(J, v, time_hi), std::runtime_error);
+
+	bertini::DefaultPrecision(30);
+	mpfr time_ok(1);                               // precision 30: coherent
+	auto out = S.Eval(v, time_ok);
+	BOOST_CHECK_EQUAL(Precision(out), 30u);
+	BOOST_CHECK_EQUAL(out(0), mpfr(0));            // (1-1)*2
+}
+
+
 BOOST_AUTO_TEST_CASE(system_jacobian)
 {
 	auto x = Variable::Make("x");

python/test/classes/system_test.py

@@ -388,3 +388,53 @@ def test_mult_system_node():
     assert np.abs(sysEval[0].imag / (0.42) - 1) <= tol_d
     assert np.abs(sysEval[1].real / (39.3240) - 1) <= tol_d
     assert np.abs(sysEval[1].imag / (-37.5584) - 1) <= tol_d
+
+
+def test_system_eval_any_point_precision():
+    """A multiprecision point of any precision is accepted by eval: the input
+    defines the working precision and the system follows it, in both directions.
+    No more pre-alignment footgun for library consumers."""
+    s = pb.parse.system('function f; variable_group x, y; f = x*y;')
+    s.precision(30)
+
+    pb.default_precision(50)
+    v = np.array((mpfr_complex(2), mpfr_complex(3)))
+    e = s.eval(v)
+    assert s.precision() == 50
+    assert e[0] == mpfr_complex(6)
+
+    pb.default_precision(20)
+    v = np.array((mpfr_complex(2), mpfr_complex(3)))
+    e = s.eval(v)
+    assert s.precision() == 20
+    assert e[0] == mpfr_complex(6)
+
+
+def test_system_eval_mixed_time_space_precision_raises_nicely():
+    """Space and time are both supplied in ONE eval call, so a precision mismatch
+    between them is caller incoherence and must raise -- with a message naming
+    BOTH precisions and saying they must be aligned."""
+    pb.default_precision(30)
+    x = Variable("mtx")
+    t = Variable("mtt")
+    s = System()
+    vg = pb.VariableGroup()
+    vg.append(x)
+    s.add_variable_group(vg)
+    s.add_path_variable(t)
+    s.add_function((1 - t) * x)
+
+    v = np.array((mpfr_complex(2),))            # precision 30
+    pb.default_precision(50)
+    time_hi = mpfr_complex(1)                   # precision 50: incoherent
+
+    with pytest.raises(RuntimeError) as excinfo:
+        s.eval(v, time_hi)
+    msg = str(excinfo.value)
+    assert '30' in msg and '50' in msg
+    assert 'align' in msg.lower()
+
+    pb.default_precision(30)
+    time_ok = mpfr_complex(1)
+    e = s.eval(v, time_ok)
+    assert e[0] == mpfr_complex(0)              # (1-1)*2