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..eea93e895
--- /dev/null
+++ b/core/include/bertini2/nag_algorithms/newton_refine.hpp
@@ -0,0 +1,152 @@
+//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/common/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
+{
+ 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{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)
+ {
+ 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 = 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 = 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..06d5638d4
--- /dev/null
+++ b/core/test/nag_algorithms/newton_refine.cpp
@@ -0,0 +1,246 @@
+//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
+
+#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::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::real_mp;
+ 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::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
+}
+
+
+// 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::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
+}
+
+
+// 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::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
+}
+
+
+BOOST_AUTO_TEST_SUITE_END()
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
new file mode 100644
index 000000000..13c213ba8
--- /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
+
+from bertini 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/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
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