feat: UI quality-of-life batch (#293–#304) + solution group projection — 3.1.0.dev1

VERSION

@@ -1 +1 @@
-3.0.0rc1
+3.1.0.dev1

core/include/bertini2/eigen_extensions.hpp

@@ -594,6 +594,44 @@ namespace bertini {
 		return Vec<NumberType>(size).unaryExpr([](NumberType const&) { return RandomUnit<NumberType>(); });
 	}
 
+
+	/**
+	\brief True if two points are the same to within \p tol in the infinity norm.
+
+	The tolerance-based point-equality test: ``max_i |a_i - b_i| <= tol``.  Points of different
+	lengths are never the same.  This is the primitive behind :func:`bertini.is_distinct_up_to` and
+	the solver's point->metadata lookup (issue #304).
+
+	\param a One point.
+	\param b The other point.
+	\param tol The (absolute) infinity-norm tolerance.
+	*/
+	template <typename DerivedA, typename DerivedB>
+	inline
+	bool IsSamePoint(Eigen::MatrixBase<DerivedA> const& a, Eigen::MatrixBase<DerivedB> const& b, double tol)
+	{
+		if (a.size() != b.size())
+			return false;
+		if (a.size() == 0)
+			return true;
+		return (a - b).template lpNorm<Eigen::Infinity>() <= tol;
+	}
+
+	/**
+	\brief True if two points differ by more than \p tol in the infinity norm (the negation of
+	:func:`IsSamePoint`).  Points of different lengths are distinct.  See issue #304.
+
+	\param a One point.
+	\param b The other point.
+	\param tol The (absolute) infinity-norm tolerance.
+	*/
+	template <typename DerivedA, typename DerivedB>
+	inline
+	bool IsDistinct(Eigen::MatrixBase<DerivedA> const& a, Eigen::MatrixBase<DerivedB> const& b, double tol)
+	{
+		return !IsSamePoint(a, b, tol);
+	}
+
 }
 
 

core/include/bertini2/nag_algorithms/zero_dim_solve.hpp

@@ -31,6 +31,7 @@
 #pragma once
 
 #include "bertini2/num_traits.hpp"
+#include "bertini2/eigen_extensions.hpp"   // bertini::IsDistinct (issues #302, #304)
 
 #include "bertini2/detail/visitable.hpp"
 #include "bertini2/tracking.hpp"
@@ -1227,14 +1228,22 @@ run the endgame, classify the endpoints, report.  See the forward-declare doc ab
 			solver's internal coordinates), mirroring SolutionsUserCoords / SolutionsInternalCoords.
 			*/
 			template<typename Pred>
-			SolnCont<Vec<BaseComplexT>> SolutionsWhere(Pred pred, bool user_coords = true) const
+			SolnCont<Vec<BaseComplexT>> SolutionsWhere(Pred pred, bool user_coords = true,
+			                                           bool merge_multiplicities = false) const
 			{
 				auto const& sols = user_coords ? SolutionsUserCoords() : SolutionsInternalCoords();
 				auto const& md   = SolutionMetadata();
 				SolnCont<Vec<BaseComplexT>> out;
 				for (size_t i = 0; i < md.size() && i < sols.size(); ++i)
+				{
+					// merge_multiplicities: keep only the chosen representative of each multiplicity
+					// cluster, so a multiplicity-m solution appears once rather than m times.  A no-op for
+					// simple / at-infinity / failed endpoints (each is its own representative).  Issue #299.
+					if (merge_multiplicities && !md[i].multiplicity_representative)
+						continue;
 					if (pred(md[i]))
 						out.push_back(sols[i]);
+				}
 				return out;
 			}
 
@@ -1244,31 +1253,31 @@ run the endgame, classify the endpoints, report.  See the forward-declare doc ab
 			Includes singular, nonsingular, and real solutions alike.
 			\see RealSolutions, SingularSolutions, NonsingularSolutions, Nonsolutions
 			*/
-			SolnCont<Vec<BaseComplexT>> FiniteSolutions(bool user_coords = true) const
+			SolnCont<Vec<BaseComplexT>> FiniteSolutions(bool user_coords = true, bool merge_multiplicities = false) const
 			{
 				return SolutionsWhere([](auto const& m){
-					return m.endgame_success_code == SuccessCode::Success && m.is_finite && !m.is_nonsolution; }, user_coords);
+					return m.endgame_success_code == SuccessCode::Success && m.is_finite && !m.is_nonsolution; }, user_coords, merge_multiplicities);
 			}
 
 			/// \brief The real finite solutions (is_real applies the configured tolerance).
-			SolnCont<Vec<BaseComplexT>> RealSolutions(bool user_coords = true) const
+			SolnCont<Vec<BaseComplexT>> RealSolutions(bool user_coords = true, bool merge_multiplicities = false) const
 			{
 				return SolutionsWhere([](auto const& m){
-					return m.endgame_success_code == SuccessCode::Success && m.is_finite && !m.is_nonsolution && m.is_real; }, user_coords);
+					return m.endgame_success_code == SuccessCode::Success && m.is_finite && !m.is_nonsolution && m.is_real; }, user_coords, merge_multiplicities);
 			}
 
 			/// \brief The nonsingular finite solutions (simple, well-conditioned roots).
-			SolnCont<Vec<BaseComplexT>> NonsingularSolutions(bool user_coords = true) const
+			SolnCont<Vec<BaseComplexT>> NonsingularSolutions(bool user_coords = true, bool merge_multiplicities = false) const
 			{
 				return SolutionsWhere([](auto const& m){
-					return m.endgame_success_code == SuccessCode::Success && m.is_finite && !m.is_nonsolution && !m.is_singular; }, user_coords);
+					return m.endgame_success_code == SuccessCode::Success && m.is_finite && !m.is_nonsolution && !m.is_singular; }, user_coords, merge_multiplicities);
 			}
 
 			/// \brief The singular finite solutions (multiple or ill-conditioned roots).
-			SolnCont<Vec<BaseComplexT>> SingularSolutions(bool user_coords = true) const
+			SolnCont<Vec<BaseComplexT>> SingularSolutions(bool user_coords = true, bool merge_multiplicities = false) const
 			{
 				return SolutionsWhere([](auto const& m){
-					return m.endgame_success_code == SuccessCode::Success && m.is_finite && !m.is_nonsolution && m.is_singular; }, user_coords);
+					return m.endgame_success_code == SuccessCode::Success && m.is_finite && !m.is_nonsolution && m.is_singular; }, user_coords, merge_multiplicities);
 			}
 
 			/**
@@ -1307,6 +1316,72 @@ run the endgame, classify the endpoints, report.  See the forward-declare doc ab
 				return solution_final_metadata_;
 			}
 
+			/**
+			\brief The metadata of the solution matching \p point -- the cluster REPRESENTATIVE (issue #302).
+
+			Matches \p point against the solutions (user coordinates by default) with the infinity-norm
+			tolerance \p tol (see IsDistinct), and returns the representative's metadata.  The return type
+			never depends on the point's multiplicity: this ALWAYS returns exactly one record (it carries
+			``multiplicity``, so you still learn m).  Throws if no solution matches, or if \p point matches
+			more than one distinct cluster (reduce \p tol).
+
+			\param point The point to look up, in user (dehomogenized) coordinates unless \p user_coords is false.
+			\param tol The infinity-norm match tolerance.
+			\param user_coords Whether \p point is in user coordinates (else the solver's internal coordinates).
+			\see CoincidentMetadataForPoint
+			*/
+			SolutionMetaDataT MetadataForPoint(Vec<BaseComplexT> const& point, double tol, bool user_coords = true) const
+			{
+				auto const& sols = user_coords ? SolutionsUserCoords() : SolutionsInternalCoords();
+				auto const& md   = SolutionMetadata();
+
+				std::vector<size_t> matches;
+				for (size_t i = 0; i < md.size() && i < sols.size(); ++i)
+					if (!bertini::IsDistinct(sols[i], point, tol))
+						matches.push_back(i);
+
+				if (matches.empty())
+					throw std::runtime_error("metadata_for: no solution matches the given point within tol");
+
+				// More than one representative among the matches => tol grouped two distinct clusters.
+				std::vector<size_t> reps;
+				for (auto i : matches)
+					if (md[i].multiplicity_representative)
+						reps.push_back(i);
+				if (reps.size() > 1)
+					throw std::runtime_error("metadata_for: the point matches more than one distinct solution cluster; reduce tol");
+				if (reps.size() == 1)
+					return md[reps.front()];
+				// The representative itself fell just outside tol but its duplicates matched; return a
+				// matched record (it still carries the cluster's multiplicity).
+				return md[matches.front()];
+			}
+
+			/**
+			\brief All metadata records for the paths coincident with \p point (issue #302).
+
+			Like MetadataForPoint, but returns EVERY coincident copy's record -- their individual per-path
+			diagnostics (condition number, residual, precision) -- ALWAYS as a list (length 1 for a simple
+			root).  Throws if no solution matches.
+
+			\param point The point to look up.
+			\param tol The infinity-norm match tolerance.
+			\param user_coords Whether \p point is in user coordinates.
+			\see MetadataForPoint
+			*/
+			std::vector<SolutionMetaDataT> CoincidentMetadataForPoint(Vec<BaseComplexT> const& point, double tol, bool user_coords = true) const
+			{
+				auto const& sols = user_coords ? SolutionsUserCoords() : SolutionsInternalCoords();
+				auto const& md   = SolutionMetadata();
+				std::vector<SolutionMetaDataT> out;
+				for (size_t i = 0; i < md.size() && i < sols.size(); ++i)
+					if (!bertini::IsDistinct(sols[i], point, tol))
+						out.push_back(md[i]);
+				if (out.empty())
+					throw std::runtime_error("metadata_for: no solution matches the given point within tol");
+				return out;
+			}
+
 			/**
 			\brief Get the solutions as computed at the endgame boundary
 			*/

core/include/bertini2/system/slice.hpp

@@ -103,10 +103,12 @@ namespace bertini {
 		{
 			typedef void (*funtype) (complex_mp&, unsigned); // the type for number generation
 			// bounded-modulus draw (away from 0 and infinity), matching patches and the start systems;
-			// kept REAL so a real slice stays real.  (A deeper pass on slice generation -- the constant
-			// column and the orthogonal=false path -- is still TODO.)
+			// kept REAL so a real slice stays real.  The orthogonal path is real too (real=true below):
+			// it QR-factors a matrix of REAL units, yielding a real orthogonal coefficient block (issue
+			// #294 -- previously the orthogonal path hardcoded the complex orthonormal matrix, so a
+			// "real" slice came out complex).
 			funtype gen = bertini::multiprecision::RandomRealBoundedModulusAssign;
-			return Make(v, dim, homogeneous, orthogonal, gen);
+			return Make(v, dim, homogeneous, orthogonal, /*real=*/true, gen);
 		}
 
 		/**
@@ -117,15 +119,23 @@ namespace bertini {
 			typedef void (*funtype) (complex_mp&, unsigned); // the type for number generation
 			// bounded-modulus draw (away from 0 and infinity), matching patches and the start systems.
 			funtype gen = bertini::multiprecision::RandomComplexBoundedModulusAssign;
-			return Make(v, dim, homogeneous, orthogonal, gen);
+			return Make(v, dim, homogeneous, orthogonal, /*real=*/false, gen);
 		}
 
 		/**
 		\brief Factory for generating slices.  Generates the variable-coefficient block (optionally
 		orthonormalized by a QR factorization) and the constant column, then assembles the augmented
 		matrix the LinearFormsBlock holds.
+
+		\param v The variable group the slice is over.
+		\param dim The number of linear forms (the slice's dimension).
+		\param homogeneous Whether the slice is homogeneous (zero constant column).
+		\param orthogonal Whether to orthonormalize the coefficient block via a QR factorization.
+		\param real Whether the coefficients are real (a real orthonormal block on the orthogonal path,
+		            matching the real \p gen used on the non-orthogonal path and for the constant column).
+		\param gen The scalar generator used for the non-orthogonal coefficients and the constant column.
 		*/
-		static Slice Make(VariableGroup const& v, unsigned dim, bool homogeneous, bool orthogonal, std::function<void(complex_mp&, unsigned)> gen)
+		static Slice Make(VariableGroup const& v, unsigned dim, bool homogeneous, bool orthogonal, bool real, std::function<void(complex_mp&, unsigned)> gen)
 		{
 			const unsigned num_vars = static_cast<unsigned>(v.size());
 
@@ -136,9 +146,15 @@ namespace bertini {
 				// conjugate-orthonormal coefficient matrix (orthonormal linear forms), drawn the b1 way
 				// via RandomConjugateOrthonormalMatrix (ADR-0041) -- it generates square and truncates,
 				// so the old transpose dance is gone.  Built at max precision, like the rest of the slice.
+				// A real slice QR-factors a matrix of REAL units (a real orthogonal block); a complex
+				// slice uses complex units.  (Issue #294: the real path must stay real.)
 				auto prev_precision = DefaultPrecision();
 				DefaultPrecision(MaxPrecisionAllowed());
-				coeffs = bertini::RandomConjugateOrthonormalMatrix<complex_mp>(dim, num_vars);
+				if (real)
+					coeffs = bertini::RandomConjugateOrthonormalMatrix<real_mp>(dim, num_vars)
+					             .unaryExpr([](real_mp const& r){ return complex_mp(r); });
+				else
+					coeffs = bertini::RandomConjugateOrthonormalMatrix<complex_mp>(dim, num_vars);
 				DefaultPrecision(prev_precision);
 			}
 			else

core/include/bertini2/system/system.hpp

@@ -1468,6 +1468,93 @@ namespace bertini {
 		}
 
 
+		/**
+		 \brief Append another system's functions to this one (issue #297).
+
+		 Syntactic sugar for ``AddFunctions(other.GetNaturalFunctions())`` -- handy for building a new
+		 system up from an existing one's equations (e.g. a critical-point system that starts from the
+		 original functions).  The functions are shared nodes, so the two systems' variables line up.
+
+		 \param other The system whose functions to copy in.
+		*/
+		void CopyFunctions(System const& other)
+		{
+			AddFunctions(other.GetNaturalFunctions());
+		}
+
+
+		/**
+		 \brief The FIFO index (0-based) of a variable group -- its position in the canonical FIFO
+		 ordering (``time_order_of_variable_groups_``), across affine and projective groups together.
+
+		 Matches by variable-node identity, so pass a group whose variables are the system's own.  Throws
+		 if the group is not one of this system's groups.
+
+		 \param g The variable group to locate.
+		*/
+		unsigned FIFOIndexOfGroup(VariableGroup const& g) const
+		{
+			unsigned affine_counter = 0, hom_counter = 0, fifo_i = 0;
+			for (auto const& iter : time_order_of_variable_groups_)
+			{
+				bool match = false;
+				switch (iter)
+				{
+					case VariableGroupType::Affine:      match = (variable_groups_[affine_counter++] == g); break;
+					case VariableGroupType::Homogeneous: match = (hom_variable_groups_[hom_counter++] == g); break;
+					case VariableGroupType::Ungrouped:   break;
+					default: throw std::runtime_error("unacceptable VariableGroupType in FIFOIndexOfGroup");
+				}
+				if (match)
+					return fifo_i;
+				++fifo_i;
+			}
+			throw std::runtime_error("FIFOIndexOfGroup: the given variable group is not a group of this system");
+		}
+
+
+		/**
+		 \brief Project a user-coordinate point onto one variable group: return just that group's
+		 coordinates.
+
+		 The point is in user (dehomogenized) coordinates -- length ``NumNaturalVariables()``, laid out in
+		 the FIFO group order.  Each group occupies a contiguous slice; this returns the slice for the
+		 group at FIFO position \p group_fifo_index.  Affine groups return their (dehomogenized) affine
+		 coordinates; **projective groups are returned as-is -- their coordinates are not dehomogenized**.
+		 Handy for an augmented system (critical points, deflation) where only one group's coordinates
+		 matter.
+
+		 \tparam T the point's number type (complex_dbl or complex_mp).
+		 \param user_point The point, in user coordinates.
+		 \param group_fifo_index The FIFO position of the group to extract (see FIFOIndexOfGroup).
+		*/
+		template<typename T>
+		Vec<T> CoordinatesOfGroup(Vec<T> const& user_point, unsigned group_fifo_index) const
+		{
+			unsigned affine_counter = 0, hom_counter = 0, fifo_i = 0, offset = 0;
+			for (auto const& iter : time_order_of_variable_groups_)
+			{
+				unsigned group_size;
+				switch (iter)
+				{
+					case VariableGroupType::Affine:      group_size = static_cast<unsigned>(variable_groups_[affine_counter++].size()); break;
+					case VariableGroupType::Homogeneous: group_size = static_cast<unsigned>(hom_variable_groups_[hom_counter++].size()); break;
+					case VariableGroupType::Ungrouped:   group_size = 1; break;
+					default: throw std::runtime_error("unacceptable VariableGroupType in CoordinatesOfGroup");
+				}
+				if (fifo_i == group_fifo_index)
+				{
+					if (offset + group_size > static_cast<unsigned>(user_point.size()))
+						throw std::runtime_error("CoordinatesOfGroup: point is shorter than the system's variable structure implies");
+					return user_point.segment(offset, group_size);
+				}
+				offset += group_size;
+				++fifo_i;
+			}
+			throw std::out_of_range("CoordinatesOfGroup: group index out of range");
+		}
+
+
 
 		/**
 		 Get the affine variable groups in the problem.

core/test/classes/eigen_test.cpp

@@ -584,7 +584,39 @@ using bertini::KahanMatrix;
 		BOOST_CHECK((A-B).norm() < 1e-38);
 		BOOST_CHECK_EQUAL(Precision(B),100);
 
-	} 
+	}
+
+BOOST_AUTO_TEST_SUITE_END()
+
+
+// issue #304: IsDistinct / IsSamePoint -- the infinity-norm tolerance point comparison.
+BOOST_AUTO_TEST_SUITE(point_distinctness)
+
+BOOST_AUTO_TEST_CASE(is_distinct_infinity_norm)
+{
+	using namespace bertini;
+
+	Vec<complex_mp> a(2), b(2), c(2);
+	a << complex_mp(1), complex_mp(2);
+	b << complex_mp(1), complex_mp(2) + complex_mp("1e-9");   // 1e-9 away in one coordinate
+	c << complex_mp(1), complex_mp(3);                        // 1 away
+
+	BOOST_CHECK(!IsDistinct(a, b, 1e-6));    // within tol -> same point
+	BOOST_CHECK( IsSamePoint(a, b, 1e-6));
+	BOOST_CHECK( IsDistinct(a, c, 1e-6));    // far apart -> distinct
+	BOOST_CHECK( IsDistinct(a, b, 1e-12));   // tighter tol -> distinct
+
+	Vec<complex_mp> shorter(1);
+	shorter << complex_mp(1);
+	BOOST_CHECK(IsDistinct(a, shorter, 1e-6));   // different length -> distinct, never same
+
+	// the same predicate works on plain doubles (so callers need not think about the number type)
+	Vec<double> da(2), db(2);
+	da << 1.0, 2.0;
+	db << 1.0, 2.0 + 1e-9;
+	BOOST_CHECK(!IsDistinct(da, db, 1e-6));
+	BOOST_CHECK( IsDistinct(da, db, 1e-12));
+}
 
 BOOST_AUTO_TEST_SUITE_END()