diff --git a/VERSION b/VERSION index 11d421a81..e7a40bd76 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -3.0.0rc1 +3.1.0.dev1 diff --git a/core/include/bertini2/eigen_extensions.hpp b/core/include/bertini2/eigen_extensions.hpp index dd508eea7..4c0b01a7f 100644 --- a/core/include/bertini2/eigen_extensions.hpp +++ b/core/include/bertini2/eigen_extensions.hpp @@ -594,6 +594,44 @@ namespace bertini { return Vec(size).unaryExpr([](NumberType const&) { return RandomUnit(); }); } + + /** + \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 + inline + bool IsSamePoint(Eigen::MatrixBase const& a, Eigen::MatrixBase const& b, double tol) + { + if (a.size() != b.size()) + return false; + if (a.size() == 0) + return true; + return (a - b).template lpNorm() <= 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 + inline + bool IsDistinct(Eigen::MatrixBase const& a, Eigen::MatrixBase const& b, double tol) + { + return !IsSamePoint(a, b, tol); + } + } diff --git a/core/include/bertini2/nag_algorithms/zero_dim_solve.hpp b/core/include/bertini2/nag_algorithms/zero_dim_solve.hpp index 80c8350e1..db8bb1577 100644 --- a/core/include/bertini2/nag_algorithms/zero_dim_solve.hpp +++ b/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 - SolnCont> SolutionsWhere(Pred pred, bool user_coords = true) const + SolnCont> SolutionsWhere(Pred pred, bool user_coords = true, + bool merge_multiplicities = false) const { auto const& sols = user_coords ? SolutionsUserCoords() : SolutionsInternalCoords(); auto const& md = SolutionMetadata(); SolnCont> 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> FiniteSolutions(bool user_coords = true) const + SolnCont> 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> RealSolutions(bool user_coords = true) const + SolnCont> 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> NonsingularSolutions(bool user_coords = true) const + SolnCont> 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> SingularSolutions(bool user_coords = true) const + SolnCont> 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 const& point, double tol, bool user_coords = true) const + { + auto const& sols = user_coords ? SolutionsUserCoords() : SolutionsInternalCoords(); + auto const& md = SolutionMetadata(); + + std::vector 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 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 CoincidentMetadataForPoint(Vec const& point, double tol, bool user_coords = true) const + { + auto const& sols = user_coords ? SolutionsUserCoords() : SolutionsInternalCoords(); + auto const& md = SolutionMetadata(); + std::vector 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 */ diff --git a/core/include/bertini2/system/slice.hpp b/core/include/bertini2/system/slice.hpp index 618185b5c..2aed676b9 100644 --- a/core/include/bertini2/system/slice.hpp +++ b/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 gen) + static Slice Make(VariableGroup const& v, unsigned dim, bool homogeneous, bool orthogonal, bool real, std::function gen) { const unsigned num_vars = static_cast(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(dim, num_vars); + if (real) + coeffs = bertini::RandomConjugateOrthonormalMatrix(dim, num_vars) + .unaryExpr([](real_mp const& r){ return complex_mp(r); }); + else + coeffs = bertini::RandomConjugateOrthonormalMatrix(dim, num_vars); DefaultPrecision(prev_precision); } else diff --git a/core/include/bertini2/system/system.hpp b/core/include/bertini2/system/system.hpp index 6095b1e10..cd9b9ce3d 100644 --- a/core/include/bertini2/system/system.hpp +++ b/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 + Vec CoordinatesOfGroup(Vec 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(variable_groups_[affine_counter++].size()); break; + case VariableGroupType::Homogeneous: group_size = static_cast(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(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. diff --git a/core/test/classes/eigen_test.cpp b/core/test/classes/eigen_test.cpp index 3f04d0bc6..4b4270298 100644 --- a/core/test/classes/eigen_test.cpp +++ b/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 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 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 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() diff --git a/core/test/classes/seeded_randomness_test.cpp b/core/test/classes/seeded_randomness_test.cpp index c4bd5927d..31903ad40 100644 --- a/core/test/classes/seeded_randomness_test.cpp +++ b/core/test/classes/seeded_randomness_test.cpp @@ -153,6 +153,30 @@ BOOST_AUTO_TEST_CASE(per_path_streams_are_deterministic_and_domain_separated) BOOST_CHECK(RandomMp(30) != path7_draw); // never collides with a path stream } +// ---- issue #294: the friendly factory draws are seeded and real-when-asked ---- + +// The bounded-modulus factory draws (behind bertini.random_real / random_complex / random_vector) +// are continuous, seed-reproducible, and -- for the real draw -- genuinely real. This is the +// generic-direction tool the notebook wanted, distinct from the quantized orthonormal random_matrix. +BOOST_AUTO_TEST_CASE(bounded_modulus_factories_are_seeded_and_real_when_asked) +{ + using bertini::multiprecision::RandomRealBoundedModulus; + using bertini::multiprecision::RandomComplexBoundedModulus; + + SetGlobalSeed(42); + auto const r1 = RandomRealBoundedModulus(); + auto const c1 = RandomComplexBoundedModulus(); + + BOOST_CHECK_EQUAL(r1.imag(), real_mp(0)); // a real draw is actually real + + SetGlobalSeed(42); // same seed reproduces the draws, in order + BOOST_CHECK_EQUAL(RandomRealBoundedModulus(), r1); + BOOST_CHECK_EQUAL(RandomComplexBoundedModulus(), c1); + + SetGlobalSeed(43); // a different seed => a different draw + BOOST_CHECK(RandomRealBoundedModulus() != r1); // (not the seed-independent orthonormal footgun) +} + // ---- the acceptance test: same seed => digest-identical homotopies ---- BOOST_AUTO_TEST_CASE(same_seed_builds_digest_identical_homotopies) diff --git a/core/test/classes/slice_test.cpp b/core/test/classes/slice_test.cpp index c37559ccd..5a509adbe 100644 --- a/core/test/classes/slice_test.cpp +++ b/core/test/classes/slice_test.cpp @@ -85,6 +85,14 @@ BOOST_AUTO_TEST_CASE(slice_basic_real) BOOST_CHECK_EQUAL(s.Dimension(),2); BOOST_CHECK_EQUAL(s.NumVariables(),3); + + // issue #294: a real slice's coefficients must actually be real. Previously the orthogonal + // (default) path fell through to the COMPLEX conjugate-orthonormal matrix, so every "real" + // slice came out complex. Every entry of the augmented coefficient matrix has zero imaginary part. + auto const& C = s.Coefficients(); + for (int ii = 0; ii < C.rows(); ++ii) + for (int jj = 0; jj < C.cols(); ++jj) + BOOST_CHECK_EQUAL(C(ii,jj).imag(), real_mp(0)); } diff --git a/core/test/classes/system_test.cpp b/core/test/classes/system_test.cpp index 13dbb069c..88e888e88 100644 --- a/core/test/classes/system_test.cpp +++ b/core/test/classes/system_test.cpp @@ -1934,6 +1934,70 @@ BOOST_AUTO_TEST_CASE(make_homotopy_with_colliding_pathvar_throws) } +// issue #297: CopyFunctions appends another system's functions (sugar behind sys.copy_functions). +BOOST_AUTO_TEST_CASE(copy_functions_appends_another_systems_functions) +{ + Var x = Variable::Make("x"); + Var y = Variable::Make("y"); + + System a; + a.AddFunction(pow(x, 2) + y); + a.AddFunction(x - y); + a.AddVariableGroup(VariableGroup{x, y}); + BOOST_CHECK_EQUAL(a.GetNaturalFunctions().size(), 2u); + + // build a new system up from a's equations + System b; + b.CopyFunctions(a); + BOOST_CHECK_EQUAL(b.GetNaturalFunctions().size(), 2u); + + // and keep adding + b.AddFunction(x + y); + BOOST_CHECK_EQUAL(b.GetNaturalFunctions().size(), 3u); + BOOST_CHECK_EQUAL(a.GetNaturalFunctions().size(), 2u); // a is untouched +} + + +// Cluster G: CoordinatesOfGroup / FIFOIndexOfGroup -- project a user-coordinate point onto one +// variable group. A system with an affine group [x,y] (FIFO 0) and a projective group [u,v] (FIFO 1). +BOOST_AUTO_TEST_CASE(coordinates_of_group_projection) +{ + Var x = Variable::Make("x"); + Var y = Variable::Make("y"); + Var u = Variable::Make("u"); + Var v = Variable::Make("v"); + + System sys; + sys.AddFunction(x + y); + sys.AddFunction(u - v); + sys.AddVariableGroup(VariableGroup{x, y}); // FIFO group 0 (affine) + sys.AddHomVariableGroup(VariableGroup{u, v}); // FIFO group 1 (projective) + + // user-coordinate point [x, y, u, v] = [1, 2, 3, 4] + Vec pt(4); + pt << complex_mp(1), complex_mp(2), complex_mp(3), complex_mp(4); + + // by FIFO index: group 0 -> the affine coords (1,2); group 1 -> the projective coords (3,4) + auto g0 = sys.CoordinatesOfGroup(pt, 0); + BOOST_CHECK_EQUAL(g0.size(), 2); + BOOST_CHECK_EQUAL(g0(0), complex_mp(1)); + BOOST_CHECK_EQUAL(g0(1), complex_mp(2)); + + auto g1 = sys.CoordinatesOfGroup(pt, 1); + BOOST_CHECK_EQUAL(g1.size(), 2); + BOOST_CHECK_EQUAL(g1(0), complex_mp(3)); + BOOST_CHECK_EQUAL(g1(1), complex_mp(4)); + + // by group object (identity match) + BOOST_CHECK_EQUAL(sys.FIFOIndexOfGroup(VariableGroup{x, y}), 0u); + BOOST_CHECK_EQUAL(sys.FIFOIndexOfGroup(VariableGroup{u, v}), 1u); + + // out-of-range index / unknown group throw + BOOST_CHECK_THROW(sys.CoordinatesOfGroup(pt, 2), std::out_of_range); + BOOST_CHECK_THROW(sys.FIFOIndexOfGroup(VariableGroup{x, u}), std::runtime_error); +} + + BOOST_AUTO_TEST_SUITE_END() diff --git a/core/test/nag_algorithms/zero_dim.cpp b/core/test/nag_algorithms/zero_dim.cpp index f99f2f007..86e6d0b2b 100644 --- a/core/test/nag_algorithms/zero_dim.cpp +++ b/core/test/nag_algorithms/zero_dim.cpp @@ -723,6 +723,44 @@ BOOST_AUTO_TEST_CASE(multiplicity_representative_marks_one_per_cluster) BOOST_CHECK_EQUAL(reps2, 4u); } +// issue #299 / #302: merge_multiplicities collapses a multiplicity-m cluster to one point, and +// metadata_for(point) looks a solution's metadata up by point. {x^2, y^2} has one solution (0,0) +// of multiplicity 4 (Bezout 4 = four coincident endpoints). +BOOST_AUTO_TEST_CASE(merge_multiplicities_and_metadata_for_point) +{ + using namespace bertini; + + auto x = node::Variable::Make("x"); + auto y = node::Variable::Make("y"); + System sys; + sys.AddFunction(pow(x, 2)); + sys.AddFunction(pow(y, 2)); + sys.AddVariableGroup(VariableGroup{x, y}); + + auto zd = algorithm::ZeroDimSolver::Cauchy, System>(sys); + zd.DefaultSetup(); + zd.Solve(); + + // #299: the multiplicity-4 solution collapses to ONE representative with merge, stays 4 without. + BOOST_CHECK_EQUAL(zd.FiniteSolutions(true, /*merge_multiplicities=*/true ).size(), 1u); + BOOST_CHECK_EQUAL(zd.FiniteSolutions(true, /*merge_multiplicities=*/false).size(), 4u); + + // #302: look the (0,0) solution up by point -> the representative, carrying multiplicity 4. + auto const pt = zd.FiniteSolutions(true, /*merge_multiplicities=*/true)[0]; // (0,0) + auto const rep = zd.MetadataForPoint(pt, 1e-5); + BOOST_CHECK(rep.multiplicity_representative); + BOOST_CHECK_EQUAL(rep.multiplicity, 4); + + // coincident=true: every copy in the cluster (all four coincident endpoints). + BOOST_CHECK_EQUAL(zd.CoincidentMetadataForPoint(pt, 1e-5).size(), 4u); + + // a point matching nothing throws. + auto far = pt; + for (int i = 0; i < far.size(); ++i) + far(i) = far(i) + 1000.0; + BOOST_CHECK_THROW(zd.MetadataForPoint(far, 1e-5), std::runtime_error); +} + // Regression: a fixed-multiple (MultiplePrecisionTracker) zero-dim solve used to throw at the start // of tracking -- "start point ... differing precision from default (20!=16)" -- because the tracker // (built at DefaultPrecision) and the config-driven ambient/thread precision (DoublePrecision) diff --git a/python/bertini/__init__.py b/python/bertini/__init__.py index 2ae47ba3e..7de4c87d2 100644 --- a/python/bertini/__init__.py +++ b/python/bertini/__init__.py @@ -87,6 +87,7 @@ Named = symbolics.NamedExpression # Named(expr, "a"): a user-named subexpression System = system.System default_precision = multiprec.default_precision +is_distinct_up_to = multiprec.is_distinct_up_to # tolerance point-inequality, infinity norm (#304) # the multiprecision number types, hoisted to the top level (they also live in bertini.multiprec) from .multiprec import complex_mp, real_mp, int_mp, rational_mp @@ -104,6 +105,7 @@ from ._calculus import jacobian from .random import random_matrix +from .random import random_vector, random_real, random_complex # exact-coefficient coercion at the top level (was bertini.linalg.coefficient / as_coefficients) from ._coefficients import coefficient, coefficients @@ -123,14 +125,52 @@ from . import operators # `from bertini.operators import *` -> just the math ops +# --- sympy interop (#295) ---------------------------------------------------------------------- +# Make every function-tree node auto-convert to sympy (the `_sympy_` protocol), so sympy.sympify(node), +# sympy.Matrix(array_of_nodes), and sympy.det(J) work directly. sympy is an optional dependency; the +# import only happens when the method is actually called (so importing bertini never needs sympy). +def _node_to_sympy(self): + from bertini.sympy_bridge import to_sympy + return to_sympy(self) + + +symbolics.AbstractNode._sympy_ = _node_to_sympy + + +# `bertini.sympy_bridge` is exposed lazily: accessing it imports the module (which raises a helpful +# ImportError if sympy is missing) without making sympy a hard dependency of `import bertini`. +def __getattr__(name): + if name == 'sympy_bridge': + import importlib + return importlib.import_module('bertini.sympy_bridge') # import_module avoids re-entering __getattr__ + raise AttributeError("module {!r} has no attribute {!r}".format(__name__, name)) + + +# --- numpy-compatible elementwise helpers for the mp types (#298, #301) ------------------------ +# Vectorized real/imag/abs/conj/round/sum/norm/is_real that stay mp-native (numpy's ufuncs/attrs on +# the custom dtypes are unreliable -- see docs/source/known_gotchas.rst). These are attributes of the +# top-level module (bertini.abs, bertini.real, ...). The builtin-shadowing names (abs, round, sum) are +# deliberately kept OUT of __all__, so `from bertini import *` never clobbers the Python builtins. +from . import _numpy_helpers as _numpy_helpers +real = _numpy_helpers.real +imag = _numpy_helpers.imag +conj = _numpy_helpers.conj +norm = _numpy_helpers.norm +is_real = _numpy_helpers.is_real +abs = _numpy_helpers.abs # noqa: A001 (bertini.abs; not exported via *) +round = _numpy_helpers.round # noqa: A001 +sum = _numpy_helpers.sum # noqa: A001 + + # https://stackoverflow.com/questions/44834/what-does-all-mean-in-python # "a list of strings defining what symbols in a module will be exported when from import * is used on the module" __all__ = ['solve','save','load','annotate','solutions_of','provenance','recording','records_dir','runs','tracks','provenance_graph','plot_chain','Solution','SolveResult','records', 'Variable','variables','gather_variables','VariableGroup','Named','system','System', - 'jacobian','random_matrix','coefficient','coefficients', + 'jacobian','random_matrix','random_vector','random_real','random_complex','coefficient','coefficients', 'complex_mp','real_mp','int_mp','rational_mp', - 'nag_algorithm','default_precision', + 'nag_algorithm','default_precision','is_distinct_up_to', + 'real','imag','conj','norm','is_real', 'tracking','endgame','logging','symbolics','parse','multiprec','random','parallel', 'operators', # everyday classes hoisted to the top level diff --git a/python/bertini/_numpy_helpers.py b/python/bertini/_numpy_helpers.py new file mode 100644 index 000000000..6ea7b81d9 --- /dev/null +++ b/python/bertini/_numpy_helpers.py @@ -0,0 +1,162 @@ +# This file is part of Bertini 2. +# +# python/bertini/_numpy_helpers.py is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# python/bertini/_numpy_helpers.py is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this file. If not, see . +# +# Copyright(C) Bertini2 Development Team + +"""Vectorized element-wise helpers for the multiprecision types (issues #298, #301). + +NumPy's ``.real`` / ``np.abs`` / ``np.round`` / identity-seeded reductions are unreliable on the +custom ``real_mp`` / ``complex_mp`` dtypes -- that boundary cannot be patched in the bindings (see +``docs/source/known_gotchas.rst``). These helpers do the elementwise work themselves, over a scalar, +a list, or a numpy array, and return **mp-native** results (a numpy object array for array input), so +your values stay arbitrary-precision instead of collapsing to float64. + + bertini.real(pt) # real parts, as real_mp + bertini.abs(pt) # magnitudes, as real_mp + bertini.is_real(pt) # is every coordinate real (imag within tol)? +""" + +import numpy as _np +from decimal import Decimal as _Decimal, ROUND_HALF_EVEN as _ROUND_HALF_EVEN + +from bertini.multiprec import real_mp as _real_mp, complex_mp as _complex_mp +from bertini.multiprec import abs as _mp_abs, conj as _mp_conj + + +def _is_container(x): + return isinstance(x, (list, tuple)) or (isinstance(x, _np.ndarray) and x.ndim > 0) + + +def _elementwise(fn, x): + """Apply scalar ``fn`` over a scalar / list / numpy array, returning mp-native results + (a numpy object array, preserving shape, for container input).""" + if not _is_container(x): + return fn(x) + arr = _np.asarray(x, dtype=object) + out = _np.empty(arr.shape, dtype=object) + a_flat = arr.reshape(-1) + o_flat = out.reshape(-1) + for i in range(a_flat.size): + o_flat[i] = fn(a_flat[i]) + return out + + +# --- scalar operations (robust across complex_mp / real_mp / python numbers) ------------------ + +def _real_scalar(v): + if isinstance(v, _complex_mp): + return v.real + if isinstance(v, _real_mp): + return v + return getattr(v, 'real', v) + + +def _imag_scalar(v): + if isinstance(v, _complex_mp): + return v.imag + if isinstance(v, _real_mp): + return _real_mp(0) + return getattr(v, 'imag', 0.0) + + +def _abs_scalar(v): + if isinstance(v, (_complex_mp, _real_mp)): + return _mp_abs(v) + return abs(v) + + +def _conj_scalar(v): + if isinstance(v, _complex_mp): + return _mp_conj(v) + if isinstance(v, _real_mp): + return v + return v.conjugate() if hasattr(v, 'conjugate') else v + + +def _round_real_mp(r, decimals): + """Round a real_mp to ``decimals`` places, staying arbitrary-precision (via Decimal).""" + q = _Decimal(1).scaleb(-decimals) + return _real_mp(str(_Decimal(repr(r)).quantize(q, rounding=_ROUND_HALF_EVEN))) + + +def _round_scalar(v, decimals): + if isinstance(v, _complex_mp): + return _complex_mp(_round_real_mp(v.real, decimals), _round_real_mp(v.imag, decimals)) + if isinstance(v, _real_mp): + return _round_real_mp(v, decimals) + return round(v, decimals) + + +# --- the public helpers ----------------------------------------------------------------------- + +def real(x): + """Real part(s), as ``real_mp`` -- over a scalar / list / array (replaces numpy ``.real``).""" + return _elementwise(_real_scalar, x) + + +def imag(x): + """Imaginary part(s), as ``real_mp`` -- over a scalar / list / array (replaces numpy ``.imag``).""" + return _elementwise(_imag_scalar, x) + + +def abs(x): + """Magnitude(s), as ``real_mp`` -- over a scalar / list / array (replaces ``np.abs``).""" + return _elementwise(_abs_scalar, x) + + +def conj(x): + """Complex conjugate(s) -- over a scalar / list / array.""" + return _elementwise(_conj_scalar, x) + + +def round(x, decimals=0): + """Round to ``decimals`` places, staying mp-native -- over a scalar / list / array (``np.round``).""" + return _elementwise(lambda v: _round_scalar(v, decimals), x) + + +def is_real(point, tol=1e-10): + """Is every coordinate of ``point`` real -- i.e. is each ``|imag| < tol``? Returns a bool. + + The one-liner behind the notebook's real-solution filter:: + + just_real = [pt for pt in solutions if bertini.is_real(pt)] + """ + flat = _np.asarray(point, dtype=object).reshape(-1) + return all(float(_abs_scalar(_imag_scalar(c))) < tol for c in flat) + + +def sum(x): + """Sum of a 1-D collection, staying mp-native (sidesteps numpy's identity-reduction gotcha).""" + flat = list(_np.asarray(x, dtype=object).reshape(-1)) + if not flat: + return 0 + total = flat[0] + for v in flat[1:]: + total = total + v + return total + + +def norm(x): + """Euclidean (2-)norm of a 1-D collection, as ``real_mp`` -- ``sqrt(sum |x_i|^2)``.""" + flat = _np.asarray(x, dtype=object).reshape(-1) + acc = None + for v in flat: + a = _abs_scalar(v) + term = a * a + acc = term if acc is None else acc + term + if acc is None: + return _real_mp(0) + # mp square root via a half-power (stays arbitrary-precision) + return acc ** _real_mp("0.5") diff --git a/python/bertini/_slice_ops.py b/python/bertini/_slice_ops.py index 3f6ccf7fb..121a0eb11 100644 --- a/python/bertini/_slice_ops.py +++ b/python/bertini/_slice_ops.py @@ -57,10 +57,62 @@ def from_coefficients(cls, coefficients, variables, homogeneous=False): return _native['from_coefficients'](vg, M, homogeneous) +def _coerce_slice_variables(variables, method_name): + """Coerce the `variables` argument of a random Slice factory to a single VariableGroup (issue #293). + + Accepts a VariableGroup, a one-element list of groups (unwrapped), or a flat list of Variables + (wrapped). A list of *several* groups -- the common mistake of passing ``sys.variable_groups()`` + -- raises a precise message instead of the opaque converter TypeError. + """ + from bertini._pybertini.container import VariableGroup as _VariableGroup + from bertini._pybertini.function_tree.symbol import Variable as _Variable + + if isinstance(variables, _VariableGroup): + return variables + # Any other iterable: a list/tuple/ndarray of Variables, or a sequence of VariableGroups (e.g. the + # ListOfVariableGroup that system.variable_groups() returns -- NOT a plain Python list). + try: + seq = list(variables) + except TypeError: + seq = None + if seq is not None: + if len(seq) == 1 and isinstance(seq[0], _VariableGroup): + return seq[0] # a one-element sequence of groups -- unwrap + if seq and all(isinstance(g, _VariableGroup) for g in seq): + raise TypeError( + "Slice.{m} wants ONE variable group, but was given a sequence of {n} of them (this is " + "what system.variable_groups() returns). Pass a single group, e.g. " + "Slice.{m}(sys.variable_groups()[0], ...), or a flat list of Variables." + .format(m=method_name, n=len(seq))) + if seq and all(isinstance(v, _Variable) for v in seq): + return _VariableGroup(seq) # a flat list of Variables -- wrap + raise TypeError( + "Slice.{m}: `variables` must be a VariableGroup or a list of Variables (got a {t})" + .format(m=method_name, t=type(variables).__name__)) + + +def _make_random_slice_factory(native, method_name): + def factory(cls, variables, dim, homogeneous=False, orthogonal=True): + vg = _coerce_slice_variables(variables, method_name) + return native(vg, dim, homogeneous, orthogonal) + factory.__name__ = method_name + factory.__doc__ = ( + "Make a random {kind} slice of `dim` linear forms over `variables` (a VariableGroup or a flat " + "list of Variables). homogeneous=True zeroes the constant column; orthogonal=True (default) " + "orthonormalizes the coefficient block.".format( + kind='real' if method_name == 'random_real' else 'complex')) + return classmethod(factory) + + def install(Slice): - """Replace ``Slice.from_coefficients`` with the friendly, coercing classmethod (idempotent).""" + """Replace ``Slice.from_coefficients`` with the friendly, coercing classmethod (idempotent); + give the random factories clear variable-group coercion + errors (issue #293).""" if getattr(Slice, "_b2_slice_ops_installed", False): return _native['from_coefficients'] = Slice.from_coefficients # native static: (variables, mpfr_matrix) Slice.from_coefficients = classmethod(from_coefficients) + for _name in ('random_complex', 'random_real'): + if hasattr(Slice, _name): + _native[_name] = getattr(Slice, _name) + setattr(Slice, _name, _make_random_slice_factory(_native[_name], _name)) Slice._b2_slice_ops_installed = True diff --git a/python/bertini/_system_ops.py b/python/bertini/_system_ops.py index 90116688c..8840577b5 100644 --- a/python/bertini/_system_ops.py +++ b/python/bertini/_system_ops.py @@ -186,15 +186,76 @@ def add_slices_as_products(self, slices): return add_products_of_linears(self, factors) +def add_variable_group(self, *variables): + """Add an affine variable group -- accepts several forms, unambiguously (issue #293). + + All of these work (a single :class:`Variable`, a :class:`VariableGroup`, and a ``list`` are + mutually distinguishable, so there is no ambiguity):: + + sys.add_variable_group(x, y, z) # loose variables + sys.add_variable_group([x, y, z]) # a list + sys.add_variable_group(x) # one variable -> a one-variable group + sys.add_variable_group(bertini.VariableGroup([x, y, z])) # the explicit form + + Returns ``self``. + """ + from bertini._pybertini.container import VariableGroup as _VariableGroup + from bertini._pybertini.function_tree.symbol import Variable as _Variable + + if len(variables) == 1: + obj = variables[0] + if isinstance(obj, _VariableGroup): + vg = obj + elif isinstance(obj, _Variable): + vg = _VariableGroup([obj]) + elif isinstance(obj, (list, tuple, np.ndarray)): + elts = list(obj) + if not elts or not all(isinstance(v, _Variable) for v in elts): + raise TypeError("add_variable_group([...]): the list must be non-empty and hold only Variables") + vg = _VariableGroup(elts) + else: + raise TypeError( + f"add_variable_group does not know how to make a variable group from a " + f"{type(obj).__name__}; pass a VariableGroup, a list of Variables, or loose Variables") + else: + if not variables: + raise TypeError("add_variable_group() needs at least one variable") + if not all(isinstance(v, _Variable) for v in variables): + raise TypeError( + "add_variable_group(x, y, z): every positional argument must be a Variable; to add a " + "VariableGroup or a list, pass it as the single argument") + vg = _VariableGroup(list(variables)) + + _native['add_variable_group'](self, vg) + return self + + +def clone(self): + """A copy of this System, ready to extend (issue #296). + + Shares the immutable node DAG (variables, functions, subexpressions) with the original but gets its + own evaluation memory, so its variables line up with the original's (clone a set-up system, give the + clone different functions, concatenate the two). The copy is unsealed and independently mutable; + adding/removing functions on one does not affect the other. For a fully serialized deep copy use + ``copy.deepcopy``. + """ + from bertini.system import clone as _clone + return _clone(self) + + def install(System): """Attach the friendly methods to the bound ``System`` class (idempotent). - Captures the native ``randomize`` first so the friendly override can still reach it. + Captures the native ``randomize`` / ``add_variable_group`` first so the friendly overrides can + still reach them. """ if getattr(System, "_b2_system_ops_installed", False): return _native['randomize'] = System.randomize + _native['add_variable_group'] = System.add_variable_group System.add_functions = add_functions + System.add_variable_group = add_variable_group + System.clone = clone System.add_linear_forms = add_linear_forms System.add_linear = add_linear System.add_products_of_linears = add_products_of_linears diff --git a/python/bertini/nag_algorithm/__init__.py b/python/bertini/nag_algorithm/__init__.py index 7253dc939..0394781e7 100644 --- a/python/bertini/nag_algorithm/__init__.py +++ b/python/bertini/nag_algorithm/__init__.py @@ -261,7 +261,8 @@ def _witness_repr(self): ) -def _zerodim_to_dataframe(self, *, user_coords=True, omit_infinite=True, merge_multiplicities=True): +def _zerodim_to_dataframe(self, *, user_coords=True, omit_infinite=True, merge_multiplicities=True, + group=None): """The solve as a pandas DataFrame -- one row per solution, the "database of solutions". Columns are ``solution`` -- the whole solution point, kept in a single cell -- then every @@ -292,6 +293,10 @@ def _zerodim_to_dataframe(self, *, user_coords=True, omit_infinite=True, merge_m ``True`` by default. Pass ``False`` to keep every endpoint, including the ``m-1`` duplicate copies. The grouping is the solver's own (the C++ clustering that computes multiplicity), read off ``multiplicity_representative``; this does not re-cluster. + group : VariableGroup or int, optional + Project each solution onto one variable group -- the ``solution`` cell then holds only that + group's coordinates (Cluster G). Pass the ``VariableGroup`` object or its 0-based FIFO index. + ``None`` (default) keeps the whole point. Returns ------- @@ -330,7 +335,10 @@ def _zerodim_to_dataframe(self, *, user_coords=True, omit_infinite=True, merge_m # internal buffer, so storing the live vector (or its elements) and letting pandas read it # later would collapse every cell to one value. The copy keeps the native element type, so # a multiprecision solve loses no precision. - row = {'solution': points[i].copy()} + sol = points[i].copy() + if group is not None: + sol = system.coordinates_of(sol, group) # project onto one variable group (Cluster G) + row = {'solution': sol} for field in _SOLUTION_METADATA_FIELDS: row[field] = getattr(m, field) row['system'] = system # a reference, so rows from different solves stay identifiable @@ -353,6 +361,49 @@ def _attach_to_dataframe(): _attach_to_dataframe() +# --- group= : project each returned solution onto one variable group (Cluster G) ------------ +# +# The reusable primitive is C++ (System.coordinates_of); here we just surface a `group=` keyword on +# every point getter so `solver.solutions(group=xyz)` returns only that group's coordinates. The +# projection goes through the solver's target_system (whose FIFO layout the user-coord points follow). +_POINT_GETTERS_FOR_GROUP = ( + 'all_solutions', 'solutions', 'finite_solutions', 'real_solutions', + 'nonsingular_solutions', 'singular_solutions', 'infinite_solutions', 'nonsolutions', +) + + +def _make_group_getter(native): + def getter(self, *args, group=None, **kwargs): + points = native(self, *args, **kwargs) + if group is None: + return points + sysm = self.target_system() + return [sysm.coordinates_of(p, group) for p in points] + getter.__name__ = getattr(native, '__name__', 'solutions') + getter.__doc__ = (getattr(native, '__doc__', '') or '') + ( + "\n\ngroup=: project each returned point onto that " + "variable group -- return only its coordinates (Cluster G). Default None keeps the whole point.") + return getter + + +def _attach_group_kwarg(): + """Give every solver class a group= keyword on its point getters (idempotent).""" + for name in dir(_pybnalag): + if not name.startswith(('ZeroDimSolver', 'HomotopySolver')): + continue + cls = getattr(_pybnalag, name) + if not isinstance(cls, type) or getattr(cls, '_b2_has_group_kwarg', False): + continue + for getter_name in _POINT_GETTERS_FOR_GROUP: + native = getattr(cls, getter_name, None) + if native is not None: + setattr(cls, getter_name, _make_group_getter(native)) + cls._b2_has_group_kwarg = True + + +_attach_group_kwarg() + + # --- ZeroDimSolver: a friendly factory over the bound ZeroDimSolver classes --- # Each bound solver class is named ZeroDimSolver (the start system is NO diff --git a/python/bertini/random/__init__.py b/python/bertini/random/__init__.py index b1df9cde0..c1d29bd1f 100644 --- a/python/bertini/random/__init__.py +++ b/python/bertini/random/__init__.py @@ -33,6 +33,12 @@ def random_matrix(rows, cols, real=False, units=False, orthonormal=True, symboli ----- Reproducible via :func:`bertini.random.set_random_seed`. A projection is just a linear functional with zero constant term, so its gradient row is exactly ``random_matrix(1, n)``. + + For a **generic real** direction (e.g. a real projection), prefer :func:`bertini.random_vector` + with ``real=True``: a real *orthonormal* matrix is QR-factored from a matrix of real units (+/-1), + so its entries are **quantized** (a real ``random_matrix(3, 1)`` has entries +/-1/sqrt(3)) and look + seed-independent -- fine for conditioning, but not a generic direction. ``random_vector`` draws + continuous bounded-modulus reals instead (generic and seed-reproducible). """ if orthonormal: M = _np.array(_pybrand.conjugate_orthonormal_matrix(rows, cols, real)) diff --git a/python/bertini/symbolics/__init__.py b/python/bertini/symbolics/__init__.py index 83d8335a0..62186a526 100644 --- a/python/bertini/symbolics/__init__.py +++ b/python/bertini/symbolics/__init__.py @@ -66,21 +66,70 @@ def sqrt(x): VariableGroup.__str__ = lambda vg: '[{}]'.format(','.join([str(v) for v in vg])) -def variables(base, indices, fmt='{base}{index}'): - """Make a list of integer-indexed Variables. +def _variablegroup_matmul(self, coefficients): + """``vg @ coeffs`` -- the single linear-combination node sum_i vg[i]*coeffs[i]. + + ``coeffs`` is a length-``len(vg)`` 1-D iterable/array of exact values (see + :func:`bertini.coefficient`; Python floats are refused) or function-tree nodes. Handy for a + projection / linear functional: ``pi = vg @ bertini.random_vector(len(vg), real=True)``. + """ + import numpy as _np + from bertini._coefficients import coefficient as _coefficient + from bertini._pybertini.function_tree import AbstractNode as _AbstractNode + + c = _np.asarray(coefficients, dtype=object).ravel() + if c.size != len(self): + raise ValueError( + "vg @ coeffs: expected {} coefficients (one per variable) but got {}" + .format(len(self), c.size)) + terms = [var * (coef if isinstance(coef, _AbstractNode) else _coefficient(coef)) + for var, coef in zip(self, c)] + result = terms[0] + for t in terms[1:]: + result = result + t + return result + + +VariableGroup.__matmul__ = _variablegroup_matmul + + +def variables(base, indices=None, fmt='{base}{index}'): + """Make a list of Variables. + + Two forms: :: - base -- name prefix, e.g. 'x' - indices -- an int n (shorthand for range(n)) or any iterable of ints - fmt -- str.format template using {base} and {index}; - default '{base}{index}' gives x0, x1, x2, ... + # explicit names -- pass a list/tuple of names as the sole argument: + x, y, z = bertini.variables(['x', 'y', 'z']) + + # integer-indexed -- a name prefix + a count (or an iterable of indices): + v = bertini.variables('v', 3) # -> [v0, v1, v2] + + Parameters + ---------- + base : str or iterable of str + A name prefix (with ``indices``), OR -- when ``indices`` is omitted -- an iterable of the + explicit variable names. + indices : int or iterable of int, optional + An int ``n`` (shorthand for ``range(n)``) or any iterable of ints. Omit to use the + explicit-names form. + fmt : str + ``str.format`` template using ``{base}`` and ``{index}``; default ``'{base}{index}'`` + gives ``x0, x1, x2, ...``. Returns a ``list[Variable]``. Wrap in a VariableGroup if desired:: pb.VariableGroup(pb.variables('x', 5)) """ from bertini._pybertini.function_tree.symbol import Variable + # explicit-names form: variables(['x', 'y', 'z']) + if indices is None: + if isinstance(base, str): + raise TypeError( + "variables('x') needs a count or indices, e.g. variables('x', 3); to name variables " + "explicitly pass a list of names, e.g. variables(['x', 'y', 'z'])") + return [Variable(name) for name in base] if isinstance(indices, int): indices = range(indices) return [Variable(fmt.format(base=base, index=i)) for i in indices] diff --git a/python/docs/source/known_gotchas.rst b/python/docs/source/known_gotchas.rst index ce13f3acb..b0b7c2033 100644 --- a/python/docs/source/known_gotchas.rst +++ b/python/docs/source/known_gotchas.rst @@ -79,3 +79,28 @@ In short: anywhere you would reach for ``np.sum(a)`` / ``np.prod(a)`` / ``np.mea ``real_mp`` or ``complex_mp`` array, reach for ``np.add.reduce(a, initial=...)`` (or ``np.multiply.reduce(a, initial=...)``), ``np.dot`` / ``np.linalg.norm``, or a plain Python ``sum`` instead. + +The bertini helpers do this for you +----------------------------------- + +The same boundary makes ``arr.real`` / ``np.abs(arr)`` / ``np.round(arr)`` unreliable on these +dtypes. So bertini ships elementwise helpers that operate over a scalar, a list, or a numpy array and +return **mp-native** results (never a float64 collapse): + +.. code-block:: python + + import bertini + + bertini.real(pt) # real parts, as real_mp (replaces arr.real) + bertini.imag(pt) # imaginary parts, as real_mp + bertini.abs(pt) # magnitudes, as real_mp (replaces np.abs) + bertini.conj(pt) # complex conjugates + bertini.round(pt, 8) # rounded, staying mp (replaces np.round) + bertini.sum(pt) # sum, staying mp (sidesteps the reduction gotcha above) + bertini.norm(pt) # Euclidean norm, as real_mp + bertini.is_real(pt) # is every coordinate real (|imag| < tol)? -> bool + +They live at the top level (``bertini.abs``, ...); the builtin-shadowing names (``abs``, ``round``, +``sum``) are deliberately kept out of ``from bertini import *``, so a star-import never clobbers the +Python builtins. For the tolerance point comparison behind de-duplication, see +:func:`bertini.is_distinct_up_to`. diff --git a/python/test/classes/numpy_helpers_test.py b/python/test/classes/numpy_helpers_test.py new file mode 100644 index 000000000..69a4be3bc --- /dev/null +++ b/python/test/classes/numpy_helpers_test.py @@ -0,0 +1,63 @@ +"""Vectorized mp helpers (#298, #301): real/imag/abs/conj/round/sum/norm/is_real stay mp-native over +a scalar / list / numpy array, so values do not collapse to float64.""" + +import numpy as np +import pytest + +import bertini as pb +from bertini.multiprec import real_mp, complex_mp + + +@pytest.fixture +def pt(): + return np.array([complex_mp("1.5"), complex_mp(2, 3)], dtype=object) + + +def test_real_imag_stay_real_mp(pt): + re = pb.real(pt) + im = pb.imag(pt) + assert all(isinstance(v, real_mp) for v in re) + assert all(isinstance(v, real_mp) for v in im) + assert [str(v) for v in re] == ['1.5', '2'] + assert [str(v) for v in im] == ['0', '3'] + + +def test_abs_is_real_mp_magnitude(pt): + a = pb.abs(pt) + assert all(isinstance(v, real_mp) for v in a) + assert abs(float(a[1]) - (13 ** 0.5)) < 1e-12 # |2 + 3i| = sqrt(13) + + +def test_conj_negates_imaginary(pt): + c = pb.conj(pt) + assert complex(c[1]) == complex(2, -3) + + +def test_round_stays_mp(pt): + r = pb.round(real_mp("2.34567"), 2) + assert isinstance(r, real_mp) + assert str(r) == '2.35' + # complex rounds both parts + rc = pb.round(complex_mp("1.23456", "7.89123"), 2) + assert isinstance(rc, complex_mp) + + +def test_is_real_predicate(pt): + assert pb.is_real(pt) is False # has a 2+3i entry + assert pb.is_real([complex_mp("1.0"), complex_mp("2.0")]) is True + assert pb.is_real([complex_mp(0, 1e-20)], tol=1e-10) is True # tiny imaginary within tol + + +def test_sum_and_norm_stay_mp(): + s = pb.sum([real_mp(1), real_mp(2), real_mp(3)]) + assert isinstance(s, real_mp) + assert float(s) == 6.0 + n = pb.norm([real_mp(3), real_mp(4)]) + assert isinstance(n, real_mp) + assert abs(float(n) - 5.0) < 1e-40 # full-precision 5, not float noise + + +def test_helpers_work_on_scalars_too(): + assert isinstance(pb.real(complex_mp(1, 2)), real_mp) + assert isinstance(pb.abs(complex_mp(3, 4)), real_mp) + assert abs(float(pb.abs(complex_mp(3, 4))) - 5.0) < 1e-12 diff --git a/python/test/classes/sympy_bridge_test.py b/python/test/classes/sympy_bridge_test.py index 93b325b93..280375356 100644 --- a/python/test/classes/sympy_bridge_test.py +++ b/python/test/classes/sympy_bridge_test.py @@ -186,3 +186,28 @@ def test_end_to_end_solve_matches_sympy(sxy): assert len(got) == len(want) == 2 for g in got: assert min(np.linalg.norm(g - w) for w in want) < 1e-8 + + +# --- #295: nodes auto-sympify (the _sympy_ protocol) + lazy bertini.sympy_bridge --- + +def test_node_sympifies_via_protocol(): + x, y = pb.Variable('x'), pb.Variable('y') + f = x**2 + y**2 - 1 + # sympy.sympify picks up the node's _sympy_ method directly + assert sp.sympify(f) == to_sympy(f) + assert sp.sympify(f) == sp.Symbol('x')**2 + sp.Symbol('y')**2 - 1 + + +def test_sympy_matrix_and_det_over_node_array(): + x, y = pb.Variable('x'), pb.Variable('y') + g = x**2 + y + J = np.array([[g.differentiate(x), g.differentiate(y)], + [(x * y).differentiate(x), (x * y).differentiate(y)]], dtype=object) + M = sp.Matrix(J) # each node auto-sympifies + # det [[2x, 1], [y, x]] = 2x^2 - y + assert sp.simplify(M.det() - (2 * sp.Symbol('x')**2 - sp.Symbol('y'))) == 0 + + +def test_sympy_bridge_reachable_at_top_level(): + # bertini.sympy_bridge is exposed lazily (module __getattr__) + assert pb.sympy_bridge.to_sympy is to_sympy diff --git a/python/test/classes/ui_ergonomics_test.py b/python/test/classes/ui_ergonomics_test.py new file mode 100644 index 000000000..d67ae43f0 --- /dev/null +++ b/python/test/classes/ui_ergonomics_test.py @@ -0,0 +1,139 @@ +"""UI ergonomics regressions: variable/group construction (#293), System accessors (#296/#297), +and node.eval accepting a point/dict (#300).""" + +import numpy as np +import pytest + +import bertini as pb +from bertini.multiprec import complex_mp + + +# --- #293: variables([...]) list form --------------------------------------------------------- + +def test_variables_explicit_names_list(): + x, y, z = pb.variables(['x', 'y', 'z']) + assert [str(v) for v in (x, y, z)] == ['x', 'y', 'z'] + + +def test_variables_integer_indexed_still_works(): + vs = pb.variables('v', 3) + assert [str(v) for v in vs] == ['v0', 'v1', 'v2'] + + +def test_variables_prefix_without_count_is_a_clear_error(): + with pytest.raises(TypeError): + pb.variables('x') + + +# --- #293: add_variable_group accepts several forms ------------------------------------------- + +@pytest.mark.parametrize("build", [ + lambda x, y, z: (lambda s: s.add_variable_group(x, y, z)), # loose variadic + lambda x, y, z: (lambda s: s.add_variable_group([x, y, z])), # a list + lambda x, y, z: (lambda s: s.add_variable_group(pb.VariableGroup([x, y, z]))), # explicit +]) +def test_add_variable_group_forms(build): + x, y, z = pb.variables(['x', 'y', 'z']) + sys = pb.System() + build(x, y, z)(sys) + assert sys.num_variable_groups() == 1 + assert sys.num_variables() == 3 + + +def test_add_variable_group_single_variable(): + (x,) = pb.variables(['x']) + sys = pb.System() + sys.add_variable_group(x) + assert sys.num_variable_groups() == 1 + + +def test_add_variable_group_rejects_mixed_variadic(): + x, y, z = pb.variables(['x', 'y', 'z']) + sys = pb.System() + with pytest.raises(TypeError): + sys.add_variable_group(x, [y, z]) + + +# --- #293: Slice.random_* clear error on a list of groups ------------------------------------- + +def test_slice_random_unwraps_single_group_sequence(): + x, y, z = pb.variables(['x', 'y', 'z']) + sys = pb.System() + sys.add_variable_group([x, y, z]) + sys.add_functions([x * x + y * y + z * z - 1]) + s = pb.Slice.random_complex(sys.variable_groups(), 1) # a one-group sequence unwraps + assert s.dimension() == 1 + + +def test_slice_random_multi_group_sequence_raises_clear_error(): + x, y, u, v = pb.variables(['x', 'y', 'u', 'v']) + sys = pb.System() + sys.add_variable_group([x, y]) + sys.add_variable_group([u, v]) + with pytest.raises(TypeError, match="ONE variable group"): + pb.Slice.random_complex(sys.variable_groups(), 1) + + +# --- #297 / #296: System.functions(), copy_functions(), clone() ------------------------------ + +def _sphere_system(): + x, y, z = pb.variables(['x', 'y', 'z']) + sys = pb.System() + sys.add_variable_group([x, y, z]) + sys.add_functions([x * x + y * y + z * z - 1, x - y]) + return sys, (x, y, z) + + +def test_functions_returns_all_functions(): + sys, _ = _sphere_system() + fns = sys.functions() + assert len(fns) == sys.num_functions() == 2 + + +def test_copy_functions_builds_from_another_system(): + sys, (x, y, z) = _sphere_system() + other = pb.System() + other.copy_functions(sys) + assert len(other.functions()) == 2 + # add_functions(functions()) is the equivalent long form + other2 = pb.System() + other2.add_functions(sys.functions()) + assert len(other2.functions()) == 2 + + +def test_clone_returns_extendable_system(): + sys, (x, y, z) = _sphere_system() + c = sys.clone() + assert len(c.functions()) == 2 + c.add_function(x + y + z) # clone is independently extendable + assert len(c.functions()) == 3 + assert len(sys.functions()) == 2 # original untouched + + +# --- #300: node.eval accepts a point array and a dict ---------------------------------------- + +def test_eval_accepts_array_dict_and_kwargs_consistently(): + x, y, z = pb.variables(['x', 'y', 'z']) + f = x * x + y - z + pt = [complex_mp(2), complex_mp(3), complex_mp(4)] # variables() order is alphabetical: x, y, z + by_kwargs = complex(f.eval(x=pt[0], y=pt[1], z=pt[2])) + by_array = complex(f.eval(pt)) + by_nparray = complex(f.eval(np.array(pt, dtype=object))) + by_dict = complex(f.eval({x: pt[0], y: pt[1], z: pt[2]})) + by_dict_names = complex(f.eval({'x': pt[0], 'y': pt[1], 'z': pt[2]})) + assert by_kwargs == by_array == by_nparray == by_dict == by_dict_names == (2 * 2 + 3 - 4) + + +def test_eval_returns_mp_native_not_float(): + x = pb.Variable('x') + f = x * x + out = f.eval([complex_mp(3)]) + assert isinstance(out, complex_mp) # NOT a python float -- the projection through-line + assert isinstance(out.real, type(out.real)) # .real is real_mp + + +def test_eval_array_wrong_length_raises(): + x, y = pb.variables(['x', 'y']) + f = x + y + with pytest.raises(RuntimeError): + f.eval([complex_mp(1)]) # 1 value, 2 variables diff --git a/python/test/random/factories_test.py b/python/test/random/factories_test.py new file mode 100644 index 000000000..bf65d815c --- /dev/null +++ b/python/test/random/factories_test.py @@ -0,0 +1,56 @@ +"""Issue #294: friendly random factories (random_real / random_complex / random_vector) and the +``vg @ coeffs`` projection sugar. + +A ``random_real`` is a genuine ``real_mp`` (not a ``complex_mp`` with zero imaginary part), the draws +are seed-reproducible and continuous (so a real projection direction is generic, unlike the quantized +orthonormal ``random_matrix``), and ``vg @ coeffs`` builds the single linear-combination node. +""" + +import numpy as np + +import bertini as pb +from bertini.multiprec import real_mp, complex_mp + + +def test_random_real_is_real_mp(): + r = pb.random_real() + assert isinstance(r, real_mp) + + +def test_random_complex_is_complex_mp(): + c = pb.random_complex() + assert isinstance(c, complex_mp) + + +def test_random_vector_element_types(): + vr = pb.random_vector(4, real=True) + vc = pb.random_vector(4, real=False) + assert len(vr) == 4 and len(vc) == 4 + assert all(isinstance(x, real_mp) for x in vr) + assert all(isinstance(x, complex_mp) for x in vc) + + +def test_random_vector_seed_reproducible_and_seed_sensitive(): + pb.random.set_random_seed(1234) + a = [repr(x) for x in pb.random_vector(5, real=True)] + + pb.random.set_random_seed(1234) + b = [repr(x) for x in pb.random_vector(5, real=True)] + assert a == b # same seed -> same vector + + pb.random.set_random_seed(4321) + c = [repr(x) for x in pb.random_vector(5, real=True)] + assert a != c # a different seed -> a different vector + + +def test_variablegroup_matmul_is_the_linear_combination(): + x, y, z = pb.Variable('x'), pb.Variable('y'), pb.Variable('z') + vg = pb.VariableGroup([x, y, z]) + coeffs = pb.random_vector(3, real=True) + + pi = vg @ coeffs # sum_i coeffs[i] * vg[i], a single node + + pt = [complex_mp(2), complex_mp(-3), complex_mp(5)] + got = complex(pi.eval(x=pt[0], y=pt[1], z=pt[2])) + want = complex(coeffs[0]) * 2 + complex(coeffs[1]) * (-3) + complex(coeffs[2]) * 5 + assert abs(got - want) < 1e-9, (got, want) diff --git a/python/test/zero_dim/observer_temp_system_test.py b/python/test/zero_dim/observer_temp_system_test.py index a74b1c4a4..ae458d210 100644 --- a/python/test/zero_dim/observer_temp_system_test.py +++ b/python/test/zero_dim/observer_temp_system_test.py @@ -38,7 +38,10 @@ def test_observer_with_temporary_system_does_not_crash(): zd.add_observer(collector) zd.solve() # must not segfault - assert len(zd.finite_solutions()) == 4 # (1,1) [singular], (-2,4), (4,16) + # finite_solutions() merges multiplicities by default (#299): (1,1) is a singular double root, + # so it counts once -> 3 distinct solutions {(1,1), (-2,4), (4,16)}; without merging it is 4 endpoints. + assert len(zd.finite_solutions()) == 3 + assert len(zd.finite_solutions(merge_multiplicities=False)) == 4 assert len(collector.series) == 6 # every path collected # the collected per-path data is usable after the solve (touches the shared nodes again). # as_dataframe() needs the optional pandas dependency (absent in CI); fall back to the diff --git a/python/test/zero_dim/solutions_ui_test.py b/python/test/zero_dim/solutions_ui_test.py new file mode 100644 index 000000000..1f1f0125b --- /dev/null +++ b/python/test/zero_dim/solutions_ui_test.py @@ -0,0 +1,118 @@ +"""Solution-surface UI: dedup by default (#299), metadata-by-point (#302), tolerance point +comparison (#304), and group projection of solutions (Cluster G).""" + +import numpy as np +import pytest + +import bertini as pb +from bertini import ZeroDimSolver + + +def _solve(functions, groups): + sys = pb.System() + for g in groups: + sys.add_variable_group(g) + sys.add_functions(functions) + solver = ZeroDimSolver(sys) + solver.solve() + return sys, solver + + +@pytest.fixture +def double_root(): + # {x^2, y^2} -> the single solution (0,0) of multiplicity 4 + x, y = pb.variables(['x', 'y']) + sys, solver = _solve([x * x, y * y], [[x, y]]) + return sys, solver + + +@pytest.fixture +def four_simple_roots(): + # x^2 - 1, y^2 - 1 -> four distinct simple roots (+-1, +-1) + x, y = pb.variables(['x', 'y']) + sys, solver = _solve([x * x - 1, y * y - 1], [[x, y]]) + return sys, solver + + +# --- #299: dedup multiplicities by default --------------------------------------------------- + +def test_solutions_dedup_multiplicities_by_default(double_root): + _, solver = double_root + assert len(solver.finite_solutions()) == 1 # merged (the default) + assert len(solver.finite_solutions(merge_multiplicities=False)) == 4 + assert len(solver.solutions()) == 1 + assert len(solver.solutions(merge_multiplicities=False)) == 4 + + +def test_simple_roots_unaffected_by_merge(four_simple_roots): + _, solver = four_simple_roots + assert len(solver.finite_solutions()) == 4 + assert len(solver.finite_solutions(merge_multiplicities=False)) == 4 + + +# --- #302: metadata_for(point) --------------------------------------------------------------- + +def test_metadata_for_point_returns_representative(double_root): + _, solver = double_root + pt = solver.finite_solutions()[0] # (0,0) + md = solver.metadata_for(pt, tol=1e-5) + assert md.multiplicity == 4 # a single record, learns m + assert md.multiplicity_representative + + +def test_metadata_for_point_coincident_is_a_list(double_root): + _, solver = double_root + pt = solver.finite_solutions()[0] + all_md = solver.metadata_for(pt, tol=1e-5, coincident=True) + assert isinstance(all_md, list) + assert len(all_md) == 4 # all coincident copies + + +def test_metadata_for_missing_point_raises(double_root): + _, solver = double_root + far = solver.finite_solutions()[0].copy() + for i in range(len(far)): + far[i] = far[i] + 1000 + with pytest.raises(RuntimeError): + solver.metadata_for(far, tol=1e-5) + + +# --- #304: is_distinct_up_to on solution points ---------------------------------------------- + +def test_is_distinct_up_to_tells_solutions_apart(four_simple_roots): + _, solver = four_simple_roots + sols = solver.finite_solutions() + assert pb.is_distinct_up_to(sols[0], sols[1], 1e-6) # different roots + assert not pb.is_distinct_up_to(sols[0], sols[0], 1e-6) # same point + + +# --- Cluster G: project solutions onto one variable group ------------------------------------ + +def test_group_projection_of_solutions(): + # two affine groups [x] and [y]; x^2-1, y^2-1 -> (+-1, +-1) + x, y = pb.variables(['x', 'y']) + sys, solver = _solve([x * x - 1, y * y - 1], [[x], [y]]) + full = solver.finite_solutions() + assert len(full) == 4 and all(len(p) == 2 for p in full) + + # by FIFO index + g0 = solver.finite_solutions(group=0) + assert all(len(p) == 1 for p in g0) # just the x coordinate + + # by VariableGroup object + vg1 = sys.variable_groups()[1] + g1 = solver.finite_solutions(group=vg1) + assert all(len(p) == 1 for p in g1) + + # the x-coordinates (deduped over the group) are {+1, -1} + xs = sorted({round(float(p[0].real)) for p in g0}) + assert xs == [-1, 1] + + +def test_group_projection_in_to_dataframe(): + pd = pytest.importorskip('pandas') + x, y = pb.variables(['x', 'y']) + sys, solver = _solve([x * x - 1, y * y - 1], [[x], [y]]) + df = solver.to_dataframe(group=0) + assert len(df) == 4 + assert all(len(sol) == 1 for sol in df.solution) # solution cell holds only group 0 diff --git a/python_bindings/include/zero_dim_export.hpp b/python_bindings/include/zero_dim_export.hpp index 2e255c720..b80d38641 100644 --- a/python_bindings/include/zero_dim_export.hpp +++ b/python_bindings/include/zero_dim_export.hpp @@ -289,7 +289,7 @@ void ZDVisitor::visit(PyClass& cl) const "get ALL the computed solutions, one per tracked path (finite, at-infinity, and failed alike). by default they are in the coordinates of YOUR variables (dehomogenized, depatched). pass user_coords=False to decline, getting the solver's internal coordinates instead: homogenized, lying on the target system's patch -- the representation to use for continuing work. the container is computed at most once per solve; repeated calls and indexing do not recompute it. for the filtered view see solutions(); for the at-infinity ones see infinite_solutions.") .def("solutions", +[](AlgoT const& self, bool singular, bool real, bool nonreal, bool nonsingular, - bool infinite, bool nonsolution, bool user_coords){ + bool infinite, bool nonsolution, bool user_coords, bool merge_multiplicities){ boost::python::list out; // Each category flag toggles inclusion of one kind of endpoint. A genuine finite // solution is returned iff its singular-class AND its real-class are both enabled; the @@ -302,14 +302,14 @@ void ZDVisitor::visit(PyClass& cl) const bool real_ok = m.is_real ? real : nonreal; return sing_ok && real_ok; }; - for (auto const& p : self.SolutionsWhere(pred, user_coords)) out.append(p); + for (auto const& p : self.SolutionsWhere(pred, user_coords, merge_multiplicities)) out.append(p); return out; }, (boost::python::arg("self"), boost::python::arg("singular") = true, boost::python::arg("real") = true, boost::python::arg("nonreal") = true, boost::python::arg("nonsingular") = true, boost::python::arg("infinite") = false, boost::python::arg("nonsolution") = false, - boost::python::arg("user_coords") = true), + boost::python::arg("user_coords") = true, boost::python::arg("merge_multiplicities") = true), "the solutions, filtered by category (returns points, not metadata). By DEFAULT every finite " "genuine solution -- real and complex, simple and multiple -- and nothing else. Each keyword " "toggles a category: singular / nonsingular select by conditioning, real / nonreal by realness " @@ -317,41 +317,43 @@ void ZDVisitor::visit(PyClass& cl) const "enabled), infinite=True also returns the at-infinity endpoints, and nonsolution=True also " "returns the nonsolutions. E.g. solutions(real=False) -> complex finite solutions only; " "solutions(singular=False) -> nonsingular finite solutions; solutions(infinite=True) adds the " - "divergent paths. user_coords=False gives the solver's internal coordinates. See also " - "real_solutions / nonsingular_solutions / singular_solutions / infinite_solutions / nonsolutions " - "for the common single-category views, and all_solutions for the raw per-path list.") + "divergent paths. merge_multiplicities=True (the DEFAULT) collapses a multiplicity-m solution " + "to its single representative (pass False to get all m coincident copies). user_coords=False " + "gives the solver's internal coordinates. See also real_solutions / nonsingular_solutions / " + "singular_solutions / infinite_solutions / nonsolutions for the common single-category views, " + "and all_solutions for the raw per-path list.") .def("finite_solutions", - +[](AlgoT const& self, bool user_coords){ + +[](AlgoT const& self, bool user_coords, bool merge_multiplicities){ boost::python::list out; - for (auto const& p : self.FiniteSolutions(user_coords)) out.append(p); + for (auto const& p : self.FiniteSolutions(user_coords, merge_multiplicities)) out.append(p); return out; }, - (boost::python::arg("self"), boost::python::arg("user_coords") = true), - "the FINITE solutions: successful endpoints the solver calls finite (is_finite applies endpoint_finite_threshold). includes singular, nonsingular, and real solutions alike. user coordinates by default; user_coords=False for internal coordinates.") + (boost::python::arg("self"), boost::python::arg("user_coords") = true, boost::python::arg("merge_multiplicities") = true), + "the FINITE solutions: successful endpoints the solver calls finite (is_finite applies endpoint_finite_threshold). includes singular, nonsingular, and real solutions alike. merge_multiplicities=True (default) collapses each multiple solution to one representative. user coordinates by default; user_coords=False for internal coordinates.") .def("real_solutions", - +[](AlgoT const& self, bool user_coords){ + +[](AlgoT const& self, bool user_coords, bool merge_multiplicities){ boost::python::list out; - for (auto const& p : self.RealSolutions(user_coords)) out.append(p); + for (auto const& p : self.RealSolutions(user_coords, merge_multiplicities)) out.append(p); return out; }, - (boost::python::arg("self"), boost::python::arg("user_coords") = true), - "the REAL finite solutions (is_real applies the configured tolerance).") + (boost::python::arg("self"), boost::python::arg("user_coords") = true, boost::python::arg("merge_multiplicities") = true), + "the REAL finite solutions (is_real applies the configured tolerance). merge_multiplicities=True (default) collapses each multiple solution to one representative.") .def("nonsingular_solutions", - +[](AlgoT const& self, bool user_coords){ + +[](AlgoT const& self, bool user_coords, bool merge_multiplicities){ boost::python::list out; - for (auto const& p : self.NonsingularSolutions(user_coords)) out.append(p); + for (auto const& p : self.NonsingularSolutions(user_coords, merge_multiplicities)) out.append(p); return out; }, - (boost::python::arg("self"), boost::python::arg("user_coords") = true), - "the NONSINGULAR finite solutions (simple, well-conditioned roots).") + (boost::python::arg("self"), boost::python::arg("user_coords") = true, boost::python::arg("merge_multiplicities") = true), + "the NONSINGULAR finite solutions (simple, well-conditioned roots). (Nonsingular solutions are simple, so merge_multiplicities is a no-op; kept for a uniform signature.)") .def("singular_solutions", - +[](AlgoT const& self, bool user_coords){ + +[](AlgoT const& self, bool user_coords, bool merge_multiplicities){ boost::python::list out; - for (auto const& p : self.SingularSolutions(user_coords)) out.append(p); + for (auto const& p : self.SingularSolutions(user_coords, merge_multiplicities)) out.append(p); return out; }, - (boost::python::arg("self"), boost::python::arg("user_coords") = true), - "the SINGULAR finite solutions (multiple or ill-conditioned roots).") + (boost::python::arg("self"), boost::python::arg("user_coords") = true, boost::python::arg("merge_multiplicities") = true), + "the SINGULAR finite solutions (multiple or ill-conditioned roots). merge_multiplicities=True (default) collapses each multiple solution to one representative.") .def("nonsolutions", +[](AlgoT const& self, bool user_coords){ boost::python::list out; @@ -375,6 +377,28 @@ void ZDVisitor::visit(PyClass& cl) const return_internal_reference<>(), "get the prepared target system: the homogenized, auto-patched clone of the system you supplied. its patch is the one internal-coordinate solutions lie on; use its dehomogenize_point/homogenize_point/variable_ordering to move between representations.") .def("solution_metadata", &AlgoT::SolutionMetadata, return_internal_reference<>(), "get the metadata for the solutions at the target time") + .def("metadata_for", + // point taken BY VALUE (a copy from numpy) -- no writable Eigen::Ref, so the ADR-0001 + // adjacent-scalar hazard does not apply. + +[](AlgoT const& self, + Vec point, + double tol, bool coincident, bool user_coords) -> boost::python::object { + if (coincident) { + boost::python::list out; + for (auto const& m : self.CoincidentMetadataForPoint(point, tol, user_coords)) out.append(m); + return out; + } + return boost::python::object(self.MetadataForPoint(point, tol, user_coords)); + }, + (boost::python::arg("self"), boost::python::arg("point"), boost::python::arg("tol"), + boost::python::arg("coincident") = false, boost::python::arg("user_coords") = true), + "the SolutionMetaData for the solution matching `point` (issue #302). Matches by the infinity-norm " + "tolerance `tol` (see is_distinct_up_to). The return type NEVER depends on the point's multiplicity: " + "by default returns exactly ONE record -- the multiplicity-cluster representative (it carries " + ".multiplicity, so you still learn m). coincident=True instead ALWAYS returns a LIST of every " + "coincident copy's record (their per-path condition number / residual / precision), length 1 for a " + "simple root. Raises if no solution matches, or if the point matches more than one distinct cluster " + "(reduce tol). `point` is in user coordinates unless user_coords=False.") .def("endgame_boundary_solutions", &AlgoT::EndgameBoundarySolutions, return_internal_reference<>(), "get the solutions (per-path point data) at the endgame boundary, where regular tracking switches to the endgame") .def("endgame_boundary_metadata", &AlgoT::EndgameBoundaryMetadata, return_internal_reference<>(), "get the MidpathCheckReport from the path-crossing check at the endgame boundary: how many crossings were detected, which paths, how many re-track attempts were made, and whether the check ultimately passed") .def("report", &AlgoT::Report, "a concise end-of-solve diagnostic summary (a SolveReport): how every path ended up -- finite solutions, diverged, or FAILED (by named reason) -- plus singular/real counts, max condition number, the path-crossing outcome, and all_paths_resolved. print(solver.report()) for a human-readable summary; a count alone can hide a path the tracker silently lost.") diff --git a/python_bindings/src/mpfr_export.cpp b/python_bindings/src/mpfr_export.cpp index d6ed10e65..e93350c0c 100644 --- a/python_bindings/src/mpfr_export.cpp +++ b/python_bindings/src/mpfr_export.cpp @@ -38,6 +38,7 @@ #include "mpfr_export.hpp" +#include "bertini2/eigen_extensions.hpp" // bertini::IsDistinct (issue #304) #include #include @@ -317,6 +318,35 @@ namespace bertini{ def("default_precision", def_prec1, "get the default precision for variable-precision numbers. is digits, not bits."); def("default_precision", def_prec2, "set the default precision for variable-precision numbers. should be a positive number. is digits, not bits."); + + // is_distinct_up_to (issue #304): tolerance-based point inequality in the infinity norm. + // Overloaded for complex_mp and real_mp vectors (points and, e.g., projection-value vectors). + // Points are taken by const& (read-only), so the ADR-0001 writable-Ref + adjacent-scalar + // hazard does not apply. + def("is_distinct_up_to", + +[](Vec const& p, Vec const& q, double tol) -> bool { + return bertini::IsDistinct(p, q, tol); }, + (arg("p"), arg("q"), arg("tol")), + "True if points p and q differ by more than tol in the infinity norm (max_i |p_i - q_i|); " + "False if they are the same to within tol. The tolerance-based point-equality test " + "(issue #304). Accepts complex_mp or real_mp vectors; different-length points are distinct."); + def("is_distinct_up_to", + +[](Vec const& p, Vec const& q, double tol) -> bool { + return bertini::IsDistinct(p, q, tol); }, + (arg("p"), arg("q"), arg("tol")), + "True if real points p and q differ by more than tol in the infinity norm (issue #304)."); + // ...and the double-precision points (complex_dbl / real_dbl), so callers need not think + // about which numeric type they are holding. + def("is_distinct_up_to", + +[](Vec> const& p, Vec> const& q, double tol) -> bool { + return bertini::IsDistinct(p, q, tol); }, + (arg("p"), arg("q"), arg("tol")), + "True if complex-double points p and q differ by more than tol in the infinity norm (issue #304)."); + def("is_distinct_up_to", + +[](Vec const& p, Vec const& q, double tol) -> bool { + return bertini::IsDistinct(p, q, tol); }, + (arg("p"), arg("q"), arg("tol")), + "True if real-double points p and q differ by more than tol in the infinity norm (issue #304)."); } diff --git a/python_bindings/src/node_export.cpp b/python_bindings/src/node_export.cpp index bf8775f09..66521abcd 100644 --- a/python_bindings/src/node_export.cpp +++ b/python_bindings/src/node_export.cpp @@ -222,26 +222,95 @@ namespace bertini{ throw std::runtime_error("could not interpret a supplied value as a number in eval"); } - // f.eval(x=2, y=5) --- evaluate this expression at a point given as keyword - // arguments naming the variables. No System is required; the values are bound - // by variable name. Every variable of the expression must be supplied, and - // every keyword must name a variable of the expression (see EvalExpression). - static object NodeEvalRaw(tuple args, dict kwargs) + // Resolve a Python object naming a variable -- a Variable node or a name string -- to its name. + static std::string VariableNameOf(object const& key) { - if (len(args) != 1) - throw std::runtime_error("eval takes the variable values as keyword arguments, e.g. f.eval(x=2, y=5)"); + extract as_str(key); + if (as_str.check()) + return as_str(); + extract> as_var(key); + if (as_var.check()) + return as_var()->name(); + throw std::runtime_error("eval: dictionary keys must be Variables or variable-name strings"); + } + // f.eval(...) --- evaluate this expression at a point. No System is required; values bind by + // variable name and the result is a complex_mp. Three call forms (issue #300): + // (a) keyword arguments naming the variables: f.eval(x=2, y=5) + // (b) a single positional dict {Variable-or-name: value}: f.eval({x: 2, y: 5}) + // (c) a single positional 1-D array / list, mapped in order to the expression's variables() + // (which are sorted by name); override that ordering with a variables= keyword: + // f.eval(pt) / f.eval(pt, variables=[x, y, z]) + // Every variable of the expression must be supplied a value. + static object NodeEvalRaw(tuple args, dict kwargs) + { + long const nargs = len(args); + if (nargs < 1) + throw std::runtime_error("eval: missing self"); std::shared_ptr self = extract>(args[0]); std::map values; - list items = dict(kwargs).items(); - for (long i = 0; i < len(items); ++i) + + // (a) keyword form -- f.eval(x=2, y=5) + if (nargs == 1) + { + list items = dict(kwargs).items(); + for (long i = 0; i < len(items); ++i) + { + object pair = items[i]; + std::string name = extract(pair[0]); + if (name == "variables") // an ordering hint is only meaningful for the array form + continue; + values[name] = CoerceToMpfrComplex(object(pair[1])); + } + return object(bertini::EvalExpression(self, values)); + } + + if (nargs > 2) + throw std::runtime_error("eval: pass either keyword arguments or a single positional point " + "(a dict, or a 1-D array/list); got too many positional arguments"); + + object point = args[1]; + + // (b) a dict {Variable-or-name: value} + extract as_dict(point); + if (as_dict.check()) + { + dict d = as_dict(); + list items = d.items(); + for (long i = 0; i < len(items); ++i) + { + object pair = items[i]; + values[VariableNameOf(object(pair[0]))] = CoerceToMpfrComplex(object(pair[1])); + } + return object(bertini::EvalExpression(self, values)); + } + + // (c) a 1-D array/list mapped to a variable ordering (variables= override, else variables()) + VariableGroup vars; + if (kwargs.has_key("variables")) + { + object vobj = kwargs["variables"]; + extract as_vg(vobj); + if (as_vg.check()) + vars = as_vg(); + else + for (long i = 0; i < len(vobj); ++i) + vars.push_back(extract>(vobj[i])()); + } + else { - object pair = items[i]; - std::string name = extract(pair[0]); - values[name] = CoerceToMpfrComplex(object(pair[1])); + vars = bertini::node::GatherVariables(self); // the expression's own variables, sorted by name } + long const n = len(point); + if (static_cast(n) != vars.size()) + throw std::runtime_error("eval: the point has " + std::to_string(n) + " entries but the " + "expression has " + std::to_string(vars.size()) + " variables; " + "pass a variables= ordering if they differ"); + for (long i = 0; i < n; ++i) + values[vars[static_cast(i)]->name()] = CoerceToMpfrComplex(object(point[i])); + return object(bertini::EvalExpression(self, values)); } @@ -256,11 +325,12 @@ namespace bertini{ class_("AbstractNode", no_init) .def(NodeVisitor()) .def("eval", raw_function(&NodeEvalRaw, 1), - "evaluate this expression at a point given as keyword arguments naming the " - "variables, e.g. f.eval(x=2, y=5). No System is needed. Evaluation is in " - "multiple precision at the current default precision; native Python floats " - "carry only float64 of information. Every variable of the expression must be " - "supplied a value, and every keyword must name a variable of the expression.") + "evaluate this expression at a point, returning a complex_mp. No System is needed. " + "Three forms (issue #300): keyword args f.eval(x=2, y=5); a dict f.eval({x: 2, y: 5}); " + "or a 1-D array/list f.eval(pt) mapped in order to the expression's variables() (sorted " + "by name) -- override the ordering with variables=, e.g. f.eval(pt, variables=[x, y, z]). " + "Evaluation is at the current default precision; native Python floats carry only float64 " + "of information. Every variable of the expression must be supplied a value.") .def("variables", &NodeVariables, (arg("self")), "The distinct variables appearing in this expression, sorted by name.") ; diff --git a/python_bindings/src/random_export.cpp b/python_bindings/src/random_export.cpp index 6d4479714..6dbe36784 100644 --- a/python_bindings/src/random_export.cpp +++ b/python_bindings/src/random_export.cpp @@ -35,6 +35,35 @@ void ExportRandom(){ def("real_unit", +[]() { return complex_mp(bertini::RandomUnit()); }, "Make a random real number of unit modulus (i.e. +1 or -1), as a complex number with imaginary part 0, in the current default precision"); + // --- friendly factories (issue #294): make a random real / complex / vector directly --- + // Continuous bounded-modulus draws (away from 0 and infinity -- the Bertini genericity draw), + // so a random projection / linear functional is generic and changes with set_random_seed. Unlike + // the orthonormal random_matrix, these are NOT quantized: a real random_vector is the right tool + // for a generic real projection direction. + def("random_real", + +[]() -> bertini::real_mp { return bertini::multiprecision::RandomRealBoundedModulus().real(); }, + "Make a random real number (real_mp) of bounded modulus (box-uniform in [-1,1], away from 0), at the current default precision. Reproducible via set_random_seed."); + + def("random_complex", + +[]() -> complex_mp { return bertini::multiprecision::RandomComplexBoundedModulus(); }, + "Make a random complex number (complex_mp) of bounded modulus (away from 0 and infinity), at the current default precision. Reproducible via set_random_seed."); + + def("random_vector", + +[](unsigned size, bool real) -> object { + if (real) { + bertini::Vec v(size); + for (unsigned i = 0; i < size; ++i) + v(i) = bertini::multiprecision::RandomRealBoundedModulus().real(); + return object(v); + } + bertini::Vec v(size); + for (unsigned i = 0; i < size; ++i) + v(i) = bertini::multiprecision::RandomComplexBoundedModulus(); + return object(v); + }, + (arg("size"), arg("real") = false), + "A random length-`size` vector of bounded-modulus numbers -- real_mp when real=True, else complex_mp -- at the current default precision. The natural random projection / linear-functional coefficient vector (generic and seed-reproducible, unlike the quantized orthonormal random_matrix)."); + def("conjugate_orthonormal_matrix", +[](unsigned rows, unsigned cols, bool real) -> bertini::Mat { if (real) { diff --git a/python_bindings/src/system_export.cpp b/python_bindings/src/system_export.cpp index c2cb4fecd..cd2a48959 100644 --- a/python_bindings/src/system_export.cpp +++ b/python_bindings/src/system_export.cpp @@ -202,6 +202,42 @@ namespace bertini{ .def("have_path_variable", &SystemBaseT::HavePathVariable, (arg("self")), "Asks whether the System has a path variable defined") .def("function", &SystemBaseT::Function, (arg("self"), arg("index")), "Get a function with a given index. Problems ensue if out of range -- uses un-rangechecked version of underlying getter") + .def("functions", + +[](SystemBaseT const& self) { + boost::python::list out; + for (auto const& f : self.GetNaturalFunctions()) out.append(f); + return out; + }, + (arg("self")), + "The system's functions, as a list of function-tree nodes (issue #297; structured blocks are expanded). So `critpt_sys.add_functions(sys.functions())` copies them all in.") + .def("copy_functions", + +[](SystemBaseT& self, SystemBaseT const& other) -> SystemBaseT& { + self.CopyFunctions(other); + return self; + }, + return_internal_reference<>(), + (arg("self"), arg("other")), + "Append another system's functions to this one and return self (issue #297; sugar for add_functions(other.functions())).") + .def("coordinates_of", + +[](SystemBaseT const& self, Vec point, boost::python::object group) -> boost::python::object { + boost::python::extract as_idx(group); + unsigned idx = as_idx.check() ? as_idx() + : self.FIFOIndexOfGroup(boost::python::extract(group)()); + return boost::python::object(self.CoordinatesOfGroup(point, idx)); + }, + (arg("self"), arg("point"), arg("group")), + "Project a user-coordinate point onto one variable group: return just that group's coordinates. " + "`group` is either the VariableGroup object or its 0-based FIFO index. Affine groups return their " + "affine coordinates; projective groups are returned as-is (not dehomogenized). Handy for an " + "augmented system (e.g. a critical-point system) where you only care about one group.") + .def("coordinates_of", + +[](SystemBaseT const& self, Vec point, boost::python::object group) -> boost::python::object { + boost::python::extract as_idx(group); + unsigned idx = as_idx.check() ? as_idx() + : self.FIFOIndexOfGroup(boost::python::extract(group)()); + return boost::python::object(self.CoordinatesOfGroup(point, idx)); + }, + (arg("self"), arg("point"), arg("group"))) .def("symbolic_jacobian", +[](SystemBaseT const& self, bool usercoordinates) { auto J = self.SymbolicJacobian(usercoordinates);