From 788d1ec58ad10e5eaff4f63e81b359386db2f90a Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Sat, 25 Jul 2026 13:15:52 +0000 Subject: [PATCH 1/3] feat(nag): standalone NewtonRefine -- tracker-free Newton against a (deflated) System The sharpening primitive: refine a point against a square autonomous (or path-variable) System to a requested consecutive-step tolerance, at the precision of the inputs. Overdetermined systems are refused with squaring instructions (randomize, as the tracking layer does). Tests walk the isosingular hierarchy on hand-deflated, hand-randomized systems (small-prime R matrices, nonsingularity verified in comments): - nonsingular root: quadratic convergence to 1e-40 - double root on the plain system: provably stalls (the disease, isolated) - double cone singularity: one deflation ([f; grad f], randomized) -> 1e-45 - whitney handle at (0,0,1): deflation + pinning slice, randomized -> 1e-45 - whitney pinch point: stage-1 deflation FAILS (still singular there); stage-2 (gradient + 2x2 minors of its Jacobian, randomized) -> 1e-45 Co-Authored-By: Claude Fable 5 --- core/CMakeLists.txt | 2 + .../bertini2/nag_algorithms/newton_refine.hpp | 149 +++++++++++ core/test/nag_algorithms/newton_refine.cpp | 244 ++++++++++++++++++ 3 files changed, 395 insertions(+) create mode 100644 core/include/bertini2/nag_algorithms/newton_refine.hpp create mode 100644 core/test/nag_algorithms/newton_refine.cpp diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt index 06ed90e95..be7e68396 100644 --- a/core/CMakeLists.txt +++ b/core/CMakeLists.txt @@ -239,6 +239,7 @@ set(io_parsing_settings_headers set(nag_algorithms_headers include/bertini2/nag_algorithms/midpath_check.hpp + include/bertini2/nag_algorithms/newton_refine.hpp include/bertini2/nag_algorithms/numerical_irreducible_decomposition.hpp include/bertini2/nag_algorithms/output.hpp include/bertini2/nag_algorithms/sharpen.hpp @@ -931,6 +932,7 @@ if (ENABLE_UNIT_TESTING) test/nag_algorithms/zero_dim_records.cpp test/nag_algorithms/numerical_irreducible_decomposition.cpp test/nag_algorithms/trace.cpp + test/nag_algorithms/newton_refine.cpp ) add_executable(test_nag_algorithms ${B2_NAG_ALGORITHMS_TEST}) diff --git a/core/include/bertini2/nag_algorithms/newton_refine.hpp b/core/include/bertini2/nag_algorithms/newton_refine.hpp new file mode 100644 index 000000000..fe0970e0b --- /dev/null +++ b/core/include/bertini2/nag_algorithms/newton_refine.hpp @@ -0,0 +1,149 @@ +//This file is part of Bertini 2. +// +//bertini2/nag_algorithms/newton_refine.hpp 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. +// +//bertini2/nag_algorithms/newton_refine.hpp 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 bertini2/nag_algorithms/newton_refine.hpp. If not, see . +// +// silviana amethyst, university of wisconsin eau claire + +/** +\file bertini2/nag_algorithms/newton_refine.hpp + +\brief Standalone Newton refinement of a point against a System -- no tracker required. + +The sharpening primitive: given a square system and an approximate solution, iterate +Newton's method until consecutive approximations agree to a requested tolerance. The +system may be autonomous (no path variable -- the ordinary case for deflated critical +point systems) or carry a path variable with a caller-supplied time. + +The whole point of the standalone (versus the tracker-owned \c Refine) is refinement on +DEFLATED systems at singular points: deflation restores quadratic convergence exactly +where the tracker's own system cannot converge, and a deflated system is a plain +autonomous System that no tracker is configured around. Overdetermined systems are +refused with an instructive message -- square them by randomization first, as the +tracking layer does. + +Evaluation runs at the precision of the supplied system and point; lift both before +calling to sharpen beyond their current precision. +*/ + +#pragma once + +#include "bertini2/system/system.hpp" +#include "bertini2/trackers/config.hpp" +#include "bertini2/linalg/lu_solver.hpp" + +namespace bertini { +namespace algorithm { + +/** +\brief The outcome of a standalone Newton refinement. + +\tparam ComplexT the complex number type the refinement ran in. +*/ +template +struct NewtonRefineResult +{ + tracking::SuccessCode code; ///< Success when the step norm reached the tolerance; FailedToConverge or MatrixSolveFailure otherwise. + Vec point; ///< The refined point (the best iterate reached, even on failure). + NumErrorT achieved; ///< Infinity norm of the last Newton step -- the consecutive-approximation agreement actually achieved. + unsigned iterations; ///< Number of Newton iterations taken. +}; + +/** +\brief Newton-refine a point against a square System, without a tracker. + +Iterates full Newton steps until the infinity norm of the step falls at or below +\p tolerance, or \p max_iterations steps have been taken. For a system with a path +variable, evaluation is at time \p time; for an autonomous system \p time is ignored. + +\tparam ComplexT the complex number type to iterate in. + +\param S the system to refine against. Must be square: \c NumTotalFunctions() + (user functions plus patches) equal to \c NumVariables(). Deflated systems are + the intended customers -- square them by randomization if overdetermined. +\param start the approximate solution to refine. +\param tolerance stop when the infinity norm of a Newton step is at or below this. +\param max_iterations refuse to iterate more than this many times. +\param time the path-variable value for non-autonomous systems; ignored otherwise. + +\return a NewtonRefineResult carrying the refined point, the achieved step norm, + the iteration count, and the SuccessCode. + +\throws std::runtime_error if the system is not square, with a message saying how + to square it. +*/ +template +NewtonRefineResult NewtonRefine(System const& S, + Vec const& start, + NumErrorT tolerance, + unsigned max_iterations, + ComplexT time = ComplexT(0)) +{ + const auto n_funcs = S.NumTotalFunctions(); + const auto n_vars = S.NumVariables(); + if (n_funcs != n_vars) + { + std::stringstream ss; + ss << "NewtonRefine requires a SQUARE system, but this one has " + << n_funcs << " total functions (user functions plus patches) over " + << n_vars << " variables. Square an overdetermined system by " + << "randomization (multiply by a generic full-rank matrix) before refining."; + throw std::runtime_error(ss.str()); + } + if (start.size() != static_cast(n_vars)) + { + std::stringstream ss; + ss << "NewtonRefine: start point has " << start.size() + << " coordinates but the system has " << n_vars << " variables."; + throw std::runtime_error(ss.str()); + } + + NewtonRefineResult result{tracking::SuccessCode::FailedToConverge, + start, static_cast(-1), 0}; + + Vec f(n_funcs); + Mat J(n_funcs, n_vars); + Vec step(n_vars); + linalg::PartialPivLU lu; + + for (unsigned it = 0; it < max_iterations; ++it) + { + if (S.HavePathVariable()) + S.SetAndReset(result.point, time); + else + S.SetAndReset(result.point); + S.EvalInPlace(f); + S.JacobianInPlace(J); + + if (lu.Factor(J) != MatrixSuccessCode::Success) + { + result.code = tracking::SuccessCode::MatrixSolveFailure; + return result; + } + lu.Solve(f, step); // step = +J^{-1} f = -(Newton step) + result.point -= step; + result.iterations = it + 1; + result.achieved = static_cast( + step.template lpNorm()); + + if (result.achieved <= tolerance) + { + result.code = tracking::SuccessCode::Success; + return result; + } + } + return result; +} + +} // namespace algorithm +} // namespace bertini diff --git a/core/test/nag_algorithms/newton_refine.cpp b/core/test/nag_algorithms/newton_refine.cpp new file mode 100644 index 000000000..a0e729560 --- /dev/null +++ b/core/test/nag_algorithms/newton_refine.cpp @@ -0,0 +1,244 @@ +//This file is part of Bertini 2. +// +//test/nag_algorithms/newton_refine.cpp 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. +// +//test/nag_algorithms/newton_refine.cpp 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 test/nag_algorithms/newton_refine.cpp. If not, see . +// +// silviana amethyst, university of wisconsin eau claire + +/** +\file test/nag_algorithms/newton_refine.cpp + +Tests for the standalone NewtonRefine -- above all, refinement of SINGULAR points on +DEFLATED systems, which is the whole point: deflation restores quadratic convergence +exactly where the plain system cannot converge. The specimens walk the isosingular +hierarchy: the double cone's isolated singularity (one deflation), a non-origin point +of the whitney umbrella's handle (a smooth point of the singular curve; one deflation +plus a pinning slice), and the umbrella's origin (the pinch point -- a singular point +OF the singular embedded curve, needing the SECOND deflation stage). + +Overdetermined deflated systems are squared by randomization with hardcoded +small-prime matrices whose nonsingularity at the target point is verified by hand in +comments -- the same recipe the tracking layer uses in production. +*/ + +#include + +#include "bertini2/nag_algorithms/newton_refine.hpp" +#include "bertini2/system/system.hpp" + +using bertini::System; +using bertini::Vec; +using bertini::complex_mp; +using bertini::DefaultPrecision; +using bertini::node::Variable; +using bertini::algorithm::NewtonRefine; +using bertini::tracking::SuccessCode; + +BOOST_AUTO_TEST_SUITE(standalone_newton_refine) + + +// a nonsingular root: full quadratic convergence, achieved accuracy at the ask +BOOST_AUTO_TEST_CASE(nonsingular_root_refines_quadratically) +{ + DefaultPrecision(60); + auto x = Variable::Make("x"); + System S; + S.AddUngroupedVariable(x); + S.AddFunction(pow(x,2) - 2); + + Vec start(1); + start << complex_mp("1.4142"); // ~1.4e-5 from sqrt(2) + + auto r = NewtonRefine(S, start, 1e-40, 50); + BOOST_CHECK(r.code == SuccessCode::Success); + BOOST_CHECK(r.achieved <= 1e-40); + BOOST_CHECK(r.iterations <= 10); // quadratic: ~3 doublings needed + using mpfr_float = bertini::mpfr_float; + mpfr_float residual = abs(pow(r.point(0),2) - complex_mp(2)); + BOOST_CHECK(residual < mpfr_float("1e-38")); +} + + +// overdetermined systems are refused with instructions, not mangled +BOOST_AUTO_TEST_CASE(overdetermined_system_is_refused) +{ + DefaultPrecision(30); + auto x = Variable::Make("x"); + System S; + S.AddUngroupedVariable(x); + S.AddFunction(pow(x,2)); + S.AddFunction(x - 1); + + Vec start(1); + start << complex_mp("0.5"); + BOOST_CHECK_THROW(NewtonRefine(S, start, 1e-20, 10), std::runtime_error); +} + + +// the disease, isolated: at a double root the plain system converges only +// linearly (steps halve), so a tight tolerance is out of reach in few iterations +BOOST_AUTO_TEST_CASE(plain_newton_stalls_at_a_double_root) +{ + DefaultPrecision(60); + auto x = Variable::Make("x"); + System S; + S.AddUngroupedVariable(x); + S.AddFunction(pow(x,2)); + + Vec start(1); + start << complex_mp("1e-6"); + auto r = NewtonRefine(S, start, 1e-40, 25); + BOOST_CHECK(r.code == SuccessCode::FailedToConverge); + BOOST_CHECK(r.achieved > 1e-40); // ~1e-6/2^25 ~ 3e-14: nowhere near +} + + +// silviana specimen 1: the double cone x^2+y^2-z^2, singular at the origin. +// Deflation appends the gradient; the overdetermined [f; grad f] (4 fns, 3 vars) +// is squared by the hardcoded randomization +// R = [ 1 2 3 5 ] +// [ 7 11 13 17 ] +// [19 23 29 31 ] +// At the origin grad f = 0 and J([f;fx;fy;fz]) has rows 0,(2,0,0),(0,2,0),(0,0,-2), +// so J(R.F)(0) has rows 2*(2,3,-5), 2*(11,13,-17), 2*(23,29,-31) with +// det(base) = -70 != 0: the deflated randomized system is REGULAR at the origin. +BOOST_AUTO_TEST_CASE(double_cone_singularity_refines_on_deflated_system) +{ + DefaultPrecision(60); + auto x = Variable::Make("x"); + auto y = Variable::Make("y"); + auto z = Variable::Make("z"); + auto f = pow(x,2) + pow(y,2) - pow(z,2); + auto fx = 2*x; + auto fy = 2*y; + auto fz = -2*z; + + System S; + bertini::VariableGroup vg{x, y, z}; + S.AddVariableGroup(vg); + S.AddFunction( 1*f + 2*fx + 3*fy + 5*fz); + S.AddFunction( 7*f + 11*fx + 13*fy + 17*fz); + S.AddFunction(19*f + 23*fx + 29*fy + 31*fz); + + Vec start(3); + start << complex_mp("1e-6"), complex_mp("-2e-6"), complex_mp("5e-7"); + + auto r = NewtonRefine(S, start, 1e-45, 50); + BOOST_CHECK(r.code == SuccessCode::Success); + BOOST_CHECK(r.achieved <= 1e-45); + using mpfr_float = bertini::mpfr_float; + mpfr_float dist = max(abs(r.point(0)), max(abs(r.point(1)), abs(r.point(2)))); + BOOST_CHECK(dist < mpfr_float("1e-40")); // landed ON the singularity +} + + +// silviana specimen 2: a NON-ORIGIN point of the whitney umbrella's handle. +// f = x^2 - y^2 z is singular along the whole handle x=y=0; the point (0,0,1) is a +// SMOOTH point of that singular curve. One deflation (the gradient +// {2x, -2yz, -y^2}) plus the pinning slice z-1 gives 4 fns over 3 vars, squared by +// R = [ 1 2 3 5 ] +// [ 7 11 13 17 ] +// [19 23 29 31 ] +// At (0,0,1) the stacked Jacobian has rows (2,0,0),(0,-2,0),(0,0,0),(0,0,1), so +// J(R.F) has rows (2,-4,5),(14,-22,17),(38,-46,31) with det = 312 != 0: REGULAR. +BOOST_AUTO_TEST_CASE(whitney_handle_point_refines_on_deflated_system) +{ + DefaultPrecision(60); + auto x = Variable::Make("x"); + auto y = Variable::Make("y"); + auto z = Variable::Make("z"); + auto g1 = 2*x; // f_x + auto g2 = -2*y*z; // f_y + auto g3 = -pow(y,2); // f_z + auto g4 = z - 1; // the pinning slice through the target point + + System S; + bertini::VariableGroup vg{x, y, z}; + S.AddVariableGroup(vg); + S.AddFunction( 1*g1 + 2*g2 + 3*g3 + 5*g4); + S.AddFunction( 7*g1 + 11*g2 + 13*g3 + 17*g4); + S.AddFunction(19*g1 + 23*g2 + 29*g3 + 31*g4); + + Vec start(3); + start << complex_mp("1e-7"), complex_mp("-1e-7"), + complex_mp(1) + complex_mp("1e-7"); + + auto r = NewtonRefine(S, start, 1e-45, 50); + BOOST_CHECK(r.code == SuccessCode::Success); + BOOST_CHECK(r.achieved <= 1e-45); + using mpfr_float = bertini::mpfr_float; + mpfr_float dist = max(abs(r.point(0)), + max(abs(r.point(1)), abs(r.point(2) - complex_mp(1)))); + BOOST_CHECK(dist < mpfr_float("1e-40")); // landed ON the handle point +} + + +// silviana specimen 3: the whitney umbrella's ORIGIN -- the pinch point, a singular +// point OF the singular embedded curve. The first deflation F1 = {2x, -2yz, -y^2} +// is itself singular there (its Jacobian has rank 1 at 0), so the SECOND stage +// appends 2x2 minors of J(F1): det[[2,0],[0,-2z]] = -4z and the row-1/row-3 +// minor -4y. F2 = {2x, -2yz, -y^2, -4z, -4y} (5 fns, 3 vars), squared by +// R = [ 1 2 3 5 7 ] +// [11 13 17 19 23 ] +// [29 31 37 41 43 ] +// At the origin the stacked Jacobian has rows (2,0,0),0,0,(0,0,-4),(0,-4,0), so +// J(R.F2)(0) has rows (2,-28,-20),(22,-92,-76),(58,-172,-164); factoring 2 from +// column 1, det(base) = -2304 != 0: the SECOND-stage system is REGULAR. +BOOST_AUTO_TEST_CASE(whitney_pinch_point_needs_and_gets_second_deflation) +{ + DefaultPrecision(60); + auto x = Variable::Make("x"); + auto y = Variable::Make("y"); + auto z = Variable::Make("z"); + auto g1 = 2*x; + auto g2 = -2*y*z; + auto g3 = -pow(y,2); + auto m1 = -4*z; // second-stage minor + auto m2 = -4*y; // second-stage minor + + // control: the FIRST-stage deflation alone, square as-is, is still singular at + // the pinch point -- Newton limps and cannot reach a tight tolerance + { + System S1; + bertini::VariableGroup vg{x, y, z}; + S1.AddVariableGroup(vg); + S1.AddFunction(g1); + S1.AddFunction(g2); + S1.AddFunction(g3); + Vec start(3); + start << complex_mp("1e-7"), complex_mp("1e-7"), complex_mp("1e-7"); + auto r1 = NewtonRefine(S1, start, 1e-45, 25); + BOOST_CHECK(r1.code != SuccessCode::Success); + } + + // the second-stage deflated randomized system nails it + System S2; + bertini::VariableGroup vg{x, y, z}; + S2.AddVariableGroup(vg); + S2.AddFunction( 1*g1 + 2*g2 + 3*g3 + 5*m1 + 7*m2); + S2.AddFunction(11*g1 + 13*g2 + 17*g3 + 19*m1 + 23*m2); + S2.AddFunction(29*g1 + 31*g2 + 37*g3 + 41*m1 + 43*m2); + + Vec start(3); + start << complex_mp("1e-7"), complex_mp("1e-7"), complex_mp("1e-7"); + + auto r = NewtonRefine(S2, start, 1e-45, 50); + BOOST_CHECK(r.code == SuccessCode::Success); + BOOST_CHECK(r.achieved <= 1e-45); + using mpfr_float = bertini::mpfr_float; + mpfr_float dist = max(abs(r.point(0)), max(abs(r.point(1)), abs(r.point(2)))); + BOOST_CHECK(dist < mpfr_float("1e-40")); // landed ON the pinch point +} + + +BOOST_AUTO_TEST_SUITE_END() From a4204759f9d0d9bc2818e50773d7c5d9a7548484 Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Sat, 25 Jul 2026 13:23:56 +0000 Subject: [PATCH 2/3] fix(nag): NewtonRefine sizes its LU workspace; bindings + python wrapper for newton_refine The stateful PartialPivLU requires ChangeSize before Factor/Solve -- without it both run over a 0-dimensional system and 'succeed' without computing (caught by the double-root stall test demanding failure and getting a phantom success). All six core tests green: nonsingular quadratic convergence, overdetermined refusal, double-root stall, and the three silviana specimens -- double cone origin (one deflation), whitney handle point (deflation + pinning slice), whitney pinch point (stage-2 deflation via Jacobian minors, after stage-1 provably fails). Sub-millisecond each. Bindings: newton_refine free function (ADR-0001 shape: vector returned, time by value), registered in bertini_python. Python: bertini.newton_refine wrapper (numpydoc), deflated-system sharpening from the port next. Co-Authored-By: Claude Fable 5 --- .../bertini2/nag_algorithms/newton_refine.hpp | 13 ++-- core/test/nag_algorithms/newton_refine.cpp | 12 +-- python/bertini/_refine.py | 65 ++++++++++++++++ python_bindings/CMakeLists.txt | 2 + .../include/newton_refine_export.hpp | 38 ++++++++++ python_bindings/src/bertini_python.cpp | 3 + python_bindings/src/newton_refine_export.cpp | 75 +++++++++++++++++++ 7 files changed, 198 insertions(+), 10 deletions(-) create mode 100644 python/bertini/_refine.py create mode 100644 python_bindings/include/newton_refine_export.hpp create mode 100644 python_bindings/src/newton_refine_export.cpp diff --git a/core/include/bertini2/nag_algorithms/newton_refine.hpp b/core/include/bertini2/nag_algorithms/newton_refine.hpp index fe0970e0b..eea93e895 100644 --- a/core/include/bertini2/nag_algorithms/newton_refine.hpp +++ b/core/include/bertini2/nag_algorithms/newton_refine.hpp @@ -39,7 +39,7 @@ calling to sharpen beyond their current precision. #pragma once #include "bertini2/system/system.hpp" -#include "bertini2/trackers/config.hpp" +#include "bertini2/common/config.hpp" #include "bertini2/linalg/lu_solver.hpp" namespace bertini { @@ -53,7 +53,7 @@ namespace algorithm { template struct NewtonRefineResult { - tracking::SuccessCode code; ///< Success when the step norm reached the tolerance; FailedToConverge or MatrixSolveFailure otherwise. + SuccessCode code; ///< Success when the step norm reached the tolerance; FailedToConverge or MatrixSolveFailure otherwise. Vec point; ///< The refined point (the best iterate reached, even on failure). NumErrorT achieved; ///< Infinity norm of the last Newton step -- the consecutive-approximation agreement actually achieved. unsigned iterations; ///< Number of Newton iterations taken. @@ -108,13 +108,16 @@ NewtonRefineResult NewtonRefine(System const& S, throw std::runtime_error(ss.str()); } - NewtonRefineResult result{tracking::SuccessCode::FailedToConverge, + NewtonRefineResult result{SuccessCode::FailedToConverge, start, static_cast(-1), 0}; Vec f(n_funcs); Mat J(n_funcs, n_vars); Vec step(n_vars); linalg::PartialPivLU lu; + lu.ChangeSize(static_cast(n_vars)); // the LU workspace is stateful: + // without this, Factor/Solve run + // over a 0-dimensional system for (unsigned it = 0; it < max_iterations; ++it) { @@ -127,7 +130,7 @@ NewtonRefineResult NewtonRefine(System const& S, if (lu.Factor(J) != MatrixSuccessCode::Success) { - result.code = tracking::SuccessCode::MatrixSolveFailure; + result.code = SuccessCode::MatrixSolveFailure; return result; } lu.Solve(f, step); // step = +J^{-1} f = -(Newton step) @@ -138,7 +141,7 @@ NewtonRefineResult NewtonRefine(System const& S, if (result.achieved <= tolerance) { - result.code = tracking::SuccessCode::Success; + result.code = SuccessCode::Success; return result; } } diff --git a/core/test/nag_algorithms/newton_refine.cpp b/core/test/nag_algorithms/newton_refine.cpp index a0e729560..06d5638d4 100644 --- a/core/test/nag_algorithms/newton_refine.cpp +++ b/core/test/nag_algorithms/newton_refine.cpp @@ -32,17 +32,19 @@ comments -- the same recipe the tracking layer uses in production. */ #include +#include #include "bertini2/nag_algorithms/newton_refine.hpp" #include "bertini2/system/system.hpp" using bertini::System; using bertini::Vec; +using bertini::Mat; using bertini::complex_mp; using bertini::DefaultPrecision; using bertini::node::Variable; using bertini::algorithm::NewtonRefine; -using bertini::tracking::SuccessCode; +using bertini::SuccessCode; BOOST_AUTO_TEST_SUITE(standalone_newton_refine) @@ -63,7 +65,7 @@ BOOST_AUTO_TEST_CASE(nonsingular_root_refines_quadratically) BOOST_CHECK(r.code == SuccessCode::Success); BOOST_CHECK(r.achieved <= 1e-40); BOOST_CHECK(r.iterations <= 10); // quadratic: ~3 doublings needed - using mpfr_float = bertini::mpfr_float; + using mpfr_float = bertini::real_mp; mpfr_float residual = abs(pow(r.point(0),2) - complex_mp(2)); BOOST_CHECK(residual < mpfr_float("1e-38")); } @@ -136,7 +138,7 @@ BOOST_AUTO_TEST_CASE(double_cone_singularity_refines_on_deflated_system) auto r = NewtonRefine(S, start, 1e-45, 50); BOOST_CHECK(r.code == SuccessCode::Success); BOOST_CHECK(r.achieved <= 1e-45); - using mpfr_float = bertini::mpfr_float; + using mpfr_float = bertini::real_mp; mpfr_float dist = max(abs(r.point(0)), max(abs(r.point(1)), abs(r.point(2)))); BOOST_CHECK(dist < mpfr_float("1e-40")); // landed ON the singularity } @@ -176,7 +178,7 @@ BOOST_AUTO_TEST_CASE(whitney_handle_point_refines_on_deflated_system) auto r = NewtonRefine(S, start, 1e-45, 50); BOOST_CHECK(r.code == SuccessCode::Success); BOOST_CHECK(r.achieved <= 1e-45); - using mpfr_float = bertini::mpfr_float; + using mpfr_float = bertini::real_mp; mpfr_float dist = max(abs(r.point(0)), max(abs(r.point(1)), abs(r.point(2) - complex_mp(1)))); BOOST_CHECK(dist < mpfr_float("1e-40")); // landed ON the handle point @@ -235,7 +237,7 @@ BOOST_AUTO_TEST_CASE(whitney_pinch_point_needs_and_gets_second_deflation) auto r = NewtonRefine(S2, start, 1e-45, 50); BOOST_CHECK(r.code == SuccessCode::Success); BOOST_CHECK(r.achieved <= 1e-45); - using mpfr_float = bertini::mpfr_float; + using mpfr_float = bertini::real_mp; mpfr_float dist = max(abs(r.point(0)), max(abs(r.point(1)), abs(r.point(2)))); BOOST_CHECK(dist < mpfr_float("1e-40")); // landed ON the pinch point } diff --git a/python/bertini/_refine.py b/python/bertini/_refine.py new file mode 100644 index 000000000..bc9bea260 --- /dev/null +++ b/python/bertini/_refine.py @@ -0,0 +1,65 @@ +"""Standalone Newton refinement -- the sharpening primitive. + +Wraps the core's tracker-free ``NewtonRefine``: iterate Newton's method on a square +system until consecutive approximations agree to a requested tolerance, at the +precision of the inputs. The intended customers are DEFLATED systems at singular +points, where deflation restores the quadratic convergence the plain system loses. +""" + +import numpy as np + +import _pybertini + +from .multiprec import complex_mp + + +def newton_refine(system, start, tolerance=1e-30, max_iterations=30, time=None): + """Newton-refine a point against a square system, without a tracker. + + Iterates full Newton steps until the infinity norm of a step falls at or below + ``tolerance`` or ``max_iterations`` steps have been taken. Evaluation runs at + the precision of the supplied system and point -- lift both (and raise + ``bertini.default_precision``) to sharpen beyond their current precision. + + Refining a SINGULAR point requires the system to be its deflation (else Newton + converges only linearly and the requested tolerance is out of reach); square an + overdetermined deflated system by randomization first, exactly as the tracking + layer does. + + Parameters + ---------- + system : bertini.System + The square system to refine against (total functions, including patches, + equal to variables). + start : array_like + The approximate solution, length equal to the system's variable count. + tolerance : float, optional + Stop when the infinity norm of a Newton step is at or below this. + max_iterations : int, optional + Refuse to iterate more than this many times. + time : complex, optional + Path-variable value for non-autonomous systems; must be omitted (or None) + for autonomous systems, which is the ordinary sharpening case. + + Returns + ------- + point : numpy.ndarray + The refined point (the best iterate reached, even on failure). + code : bertini.tracking.SuccessCode + ``Success`` when the tolerance was reached; ``FailedToConverge`` or + ``MatrixSolveFailure`` otherwise. + achieved : float + Infinity norm of the last Newton step -- the agreement actually achieved. + iterations : int + Number of Newton iterations taken. + """ + arr = np.asarray([c if isinstance(c, complex_mp) else complex_mp(c) + for c in np.asarray(start).ravel()]) + if time is None: + point, code, achieved, iterations = _pybertini.newton_refine( + system, arr, float(tolerance), int(max_iterations)) + else: + point, code, achieved, iterations = _pybertini.newton_refine( + system, arr, float(tolerance), int(max_iterations), + complex_mp(time)) + return np.asarray(point), code, achieved, iterations diff --git a/python_bindings/CMakeLists.txt b/python_bindings/CMakeLists.txt index 1c74f1624..4bca24021 100644 --- a/python_bindings/CMakeLists.txt +++ b/python_bindings/CMakeLists.txt @@ -132,6 +132,7 @@ set(PYBERTINI_HEADERS include/eigenpy_interaction.hpp include/function_tree_export.hpp include/mpfr_export.hpp + include/newton_refine_export.hpp include/random_export.hpp include/node_export.hpp include/symbol_export.hpp @@ -165,6 +166,7 @@ set(PYBERTINI_SOURCES src/endgame_double_export.cpp src/endgame_mp_export.cpp src/endgame_amp_export.cpp + src/newton_refine_export.cpp src/random_export.cpp src/mpfr_export.cpp src/linalg_export.cpp diff --git a/python_bindings/include/newton_refine_export.hpp b/python_bindings/include/newton_refine_export.hpp new file mode 100644 index 000000000..b8c5505f3 --- /dev/null +++ b/python_bindings/include/newton_refine_export.hpp @@ -0,0 +1,38 @@ +//This file is part of Bertini 2. +// +//python/newton_refine_export.hpp 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/newton_refine_export.hpp 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 python/newton_refine_export.hpp. If not, see . +// +// silviana amethyst, university of wisconsin eau claire + +/** +\file python/newton_refine_export.hpp + +\brief Exposes the standalone NewtonRefine (tracker-free Newton against a System, +deflated systems above all) to Python. +*/ + +#pragma once + +#include "python_common.hpp" + +namespace bertini { +namespace python { + +/** +\brief Export the standalone newton_refine free function to the current module scope. +*/ +void ExportNewtonRefine(); + +} // namespace python +} // namespace bertini diff --git a/python_bindings/src/bertini_python.cpp b/python_bindings/src/bertini_python.cpp index d345f1e09..b8d736b4b 100644 --- a/python_bindings/src/bertini_python.cpp +++ b/python_bindings/src/bertini_python.cpp @@ -46,6 +46,7 @@ #include "eigenpy_interaction.hpp" // EnableEigenPy() #include "parallel_export.hpp" // ExportParallel() #include "records_export.hpp" // ExportRecords() +#include "newton_refine_export.hpp" // ExportNewtonRefine() #include "bertini2/fast_allocator.hpp" // InstallFastAllocator() namespace bertini { namespace python { @@ -151,6 +152,8 @@ namespace bertini ExportRecords(); + ExportNewtonRefine(); + ExportNID(); ExportInfo(); diff --git a/python_bindings/src/newton_refine_export.cpp b/python_bindings/src/newton_refine_export.cpp new file mode 100644 index 000000000..a24b983a0 --- /dev/null +++ b/python_bindings/src/newton_refine_export.cpp @@ -0,0 +1,75 @@ +//This file is part of Bertini 2. +// +//python/newton_refine_export.cpp 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/newton_refine_export.cpp 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 python/newton_refine_export.cpp. If not, see . +// +// silviana amethyst, university of wisconsin eau claire + +#include "newton_refine_export.hpp" + +#include "bertini2/nag_algorithms/newton_refine.hpp" +#include "bertini2/system/system.hpp" + +namespace bertini { +namespace python { + +namespace { + +// ADR-0001: the scalar time is taken BY VALUE (never const& adjacent to an +// eigenpy-converted vector), and the refined point is RETURNED rather than +// written through a writable Ref. +boost::python::tuple NewtonRefineAtTime(bertini::System const& sys, + Vec const& start, + double tolerance, + unsigned max_iterations, + complex_mp time) +{ + auto result = bertini::algorithm::NewtonRefine(sys, start, tolerance, + max_iterations, time); + return boost::python::make_tuple(result.point, result.code, + result.achieved, result.iterations); +} + +boost::python::tuple NewtonRefineAutonomous(bertini::System const& sys, + Vec const& start, + double tolerance, + unsigned max_iterations) +{ + return NewtonRefineAtTime(sys, start, tolerance, max_iterations, + complex_mp(0)); +} + +} // namespace + + +void ExportNewtonRefine() +{ + boost::python::def( + "newton_refine", &NewtonRefineAutonomous, + "newton_refine(system, start, tolerance, max_iterations) -> (point, code, achieved, iterations)\n\n" + "Standalone Newton refinement of a point against a SQUARE system -- no tracker.\n" + "Iterates full Newton steps until consecutive approximations agree to `tolerance`\n" + "in the infinity norm, at the precision of the supplied system and point.\n\n" + "The intended customers are DEFLATED systems at singular points: deflation\n" + "restores quadratic convergence exactly where the plain system cannot converge.\n" + "Overdetermined systems raise with instructions to square by randomization.\n\n" + "Returns (refined point, SuccessCode, achieved step norm, iterations taken)."); + + boost::python::def( + "newton_refine", &NewtonRefineAtTime, + "newton_refine(system, start, tolerance, max_iterations, time) -> (point, code, achieved, iterations)\n\n" + "As newton_refine/4, for a system with a path variable, evaluated at `time`."); +} + +} // namespace python +} // namespace bertini From ac14884c733797792398d438bcff1a5ca09cdec2 Mon Sep 17 00:00:00 2001 From: silviana amethyst Date: Sat, 25 Jul 2026 13:28:06 +0000 Subject: [PATCH 3/3] test(python): newton_refine interface tests; wire bertini.newton_refine into the package Interface-only per doctrine (correctness is gated in C++): the cone specimen end-to-end (4 iterations, 1e-6 offset -> <1e-40 from the singularity), the overdetermined refusal surfacing as a friendly RuntimeError, and the honest FailedToConverge on a plain double root. Co-Authored-By: Claude Fable 5 --- python/bertini/__init__.py | 1 + python/bertini/_refine.py | 2 +- python/test/classes/newton_refine_test.py | 59 +++++++++++++++++++++++ 3 files changed, 61 insertions(+), 1 deletion(-) create mode 100644 python/test/classes/newton_refine_test.py diff --git a/python/bertini/__init__.py b/python/bertini/__init__.py index b39230556..b834d34cb 100644 --- a/python/bertini/__init__.py +++ b/python/bertini/__init__.py @@ -106,6 +106,7 @@ from ._calculus import jacobian from ._randomize import randomize +from ._refine import newton_refine # dense linear algebra for the mp types (solve / LU), backed by eigenpy's decompositions from . import linalg from .random import random_matrix diff --git a/python/bertini/_refine.py b/python/bertini/_refine.py index bc9bea260..13c213ba8 100644 --- a/python/bertini/_refine.py +++ b/python/bertini/_refine.py @@ -8,7 +8,7 @@ import numpy as np -import _pybertini +from bertini import _pybertini from .multiprec import complex_mp diff --git a/python/test/classes/newton_refine_test.py b/python/test/classes/newton_refine_test.py new file mode 100644 index 000000000..59dedd105 --- /dev/null +++ b/python/test/classes/newton_refine_test.py @@ -0,0 +1,59 @@ +"""Interface tests for ``bertini.newton_refine`` -- the standalone sharpening +primitive. Correctness (the isosingular-hierarchy specimens) is gated in C++ +(``core/test/nag_algorithms/newton_refine.cpp``); these check the Python surface: +argument handling, return shape, and the friendly refusal of overdetermined +systems.""" + +import numpy as np +import pytest + +import bertini +from bertini.multiprec import complex_mp as C + + +def _cone_deflated_randomized(): + """The double cone's deflated, hand-randomized (small primes) sharpening system.""" + x, y, z = bertini.variables(list('xyz')) + S = bertini.System() + S.add_variable_group([x, y, z]) + f = x**2 + y**2 - z**2 + fx, fy, fz = 2 * x, 2 * y, -2 * z + S.add_functions([1 * f + 2 * fx + 3 * fy + 5 * fz, + 7 * f + 11 * fx + 13 * fy + 17 * fz, + 19 * f + 23 * fx + 29 * fy + 31 * fz]) + return S + + +@pytest.mark.parametrize("precision", [60], indirect=True) +def test_refines_cone_singularity_on_deflated_system(precision): + S = _cone_deflated_randomized() + start = [C('1e-6'), C('-2e-6'), C('5e-7')] + pt, code, achieved, its = bertini.newton_refine( + S, start, tolerance=1e-45, max_iterations=50) + assert str(code) == 'Success' + assert achieved <= 1e-45 + assert its <= 10 # quadratic, not limping + assert isinstance(pt, np.ndarray) and pt.shape == (3,) + assert max(abs(complex(c)) for c in pt) < 1e-40 # ON the singularity + + +@pytest.mark.parametrize("precision", [60], indirect=True) +def test_overdetermined_system_raises_with_instructions(precision): + x = bertini.variables(['x'])[0] + S = bertini.System() + S.add_variable_group([x]) + S.add_functions([x**2, x - 1]) + with pytest.raises(RuntimeError, match='[Ss]quare'): + bertini.newton_refine(S, [C('0.5')], tolerance=1e-20, max_iterations=10) + + +@pytest.mark.parametrize("precision", [60], indirect=True) +def test_plain_double_root_reports_failure_honestly(precision): + x = bertini.variables(['x'])[0] + S = bertini.System() + S.add_variable_group([x]) + S.add_functions([x**2]) + pt, code, achieved, its = bertini.newton_refine( + S, [C('1e-6')], tolerance=1e-40, max_iterations=25) + assert str(code) == 'FailedToConverge' + assert achieved > 1e-40 # linear halving cannot get there