🛷 Move a slice, hold the system fixed¶
This tutorial is about a performance minded feature of Bertini 2.
Regeneration and witness-set work share a shape: most equations are fixed – the polynomial system, plus “below” linear slices that cut the dimension – and only a small part moves along the path variable. The fixed equations should be evaluated just once per step: never duplicated, never scaled by the path coefficient, never differentiated in \(t\).
straight_line_homotopy() builds exactly that. Hand it the target and start block
of rows that move and, as fixed=, the equations that do not, and it returns
keeping the fixed equations as their own evaluation blocks, and deforming only the moving parts. Those rows are the whole of the path-variable dependence: they are evaluated once per step rather than blended, and contribute exactly zero to \(dH/dt\).
It hands back the homotopy and the systems at both ends – .homotopy, .target,
.start – so neither end is yours to assemble. That matters for more than typing: the rows must
come fixed-first, and nothing would catch it if they did not.
Move one slice¶
A fixed unit circle, sliced by a line that slides from the x-axis (\(y=0\)) to the diagonal (\(y-x=0\)). At \(t=1\) the slice is the x-axis, so the start points are the circle’s intersections with it, \((\pm 1, 0)\); at \(t=0\) the slice is the diagonal:
x, y = bertini.Variable('x'), bertini.Variable('y')
fixed = bertini.System()
fixed.add_variable_group(bertini.VariableGroup([x, y]))
fixed.add_function(x*x + y*y - 1) # the fixed unit circle
start_moving = bertini.System() # the slice at t=1: the x-axis y = 0
start_moving.add_variable_group(bertini.VariableGroup([x, y]))
start_moving.add_function(y)
end_moving = bertini.System() # the slice at t=0: the diagonal y - x = 0
end_moving.add_variable_group(bertini.VariableGroup([x, y]))
end_moving.add_function(y - x)
gamma = bertini.coefficient(bertini.multiprec.complex_mp('0.6', '0.8')) # off the real axis
b = nag_algorithm.straight_line_homotopy(end_moving, start_moving, fixed=fixed, gamma=gamma)
start_points = [np.array([bertini.multiprec.complex_mp('1'), bertini.multiprec.complex_mp('0')]),
np.array([bertini.multiprec.complex_mp('-1'), bertini.multiprec.complex_mp('0')])]
solver = bertini.HomotopySolver(b, start_points) # b.target comes with it
solver.solve()
roots = sorted((round(complex(s[0]).real, 4), round(complex(s[1]).real, 4))
for s in solver.all_solutions())
r = round(1 / np.sqrt(2), 4)
assert roots == sorted([(r, r), (-r, -r)]) # circle ∩ diagonal = (±1/√2, ±1/√2)
The circle never moved; only the slice did. The two start points slid along the circle to the two diagonal intersections.
A static slice and a moving slice¶
Now in three variables, the genuinely regeneration-flavored case: a fixed unit sphere (a
surface) cut to dimension zero by two slices – a static one z = 0 and a moving one
that slides y = 0 \(\to\) y - x = 0. Two blocks are fixed (the sphere and the static
slice); only the moving slice carries \(t\):
x, y, z = bertini.Variable('x'), bertini.Variable('y'), bertini.Variable('z')
fixed = bertini.System()
fixed.add_variable_group(bertini.VariableGroup([x, y, z]))
fixed.add_function(x*x + y*y + z*z - 1) # the sphere
fixed.add_linear(np.array([[0, 0, 1]]), np.array([x, y, z])) # static slice z = 0
start_moving = bertini.System(); start_moving.add_variable_group(bertini.VariableGroup([x, y, z]))
start_moving.add_function(y) # moving slice at t=1
end_moving = bertini.System(); end_moving.add_variable_group(bertini.VariableGroup([x, y, z]))
end_moving.add_function(y - x) # moving slice at t=0
b = nag_algorithm.straight_line_homotopy(end_moving, start_moving, fixed=fixed, gamma=gamma)
start_points = [np.array([bertini.multiprec.complex_mp(str(a)), bertini.multiprec.complex_mp('0'),
bertini.multiprec.complex_mp('0')]) for a in (1, -1)]
solver = bertini.HomotopySolver(b, start_points)
solver.solve()
roots = sorted((round(complex(s[0]).real, 4), round(complex(s[1]).real, 4), round(complex(s[2]).real, 4))
for s in solver.all_solutions())
r = round(1 / np.sqrt(2), 4)
assert roots == sorted([(r, r, 0.0), (-r, -r, 0.0)])
The fixed system is left out of the homotopy’s “motion”. Ask the homotopy for \(dH/dt\): the rows of the two fixed blocks (the sphere and the static slice) are exactly zero – they are evaluated once and never differentiated as the slice moves – while only the moving row can be nonzero:
# affine here, so the rows line up with the three equations above: the builder projectivizes
# by default, which adds a homogenizing coordinate and a patch row (see below).
Haff = nag_algorithm.straight_line_homotopy(end_moving, start_moving, fixed=fixed,
gamma=gamma, projectivize=False).homotopy
pt = np.array([bertini.multiprec.complex_mp('0.3'),
bertini.multiprec.complex_mp('0.4'),
bertini.multiprec.complex_mp('0.5')])
# H evaluates at the precision of the point it is given; nothing to align first.
dHdt = Haff.eval_time_derivative(pt, bertini.multiprec.complex_mp('0.5'))
assert abs(complex(dHdt[0])) == 0.0 # sphere row: out of dH/dt
assert abs(complex(dHdt[1])) == 0.0 # static slice row: out of dH/dt
assert abs(complex(dHdt[2])) > 0.0 # only the moving slice carries t
Projective by default¶
The homotopy above is built over projective coordinates: the builder homogenizes and patches the three systems before combining them, so a path heading to infinity reaches an ordinary point instead of running off. That is what the zero-dimensional solver has always done with the systems it builds for itself, and a homotopy can only be made projective where it is built – once its blend exists, the operand systems are fixed (for performance and recording reasons).
You can see projectivity in the shape: a homogenizing coordinate per affine group, and a patch row.
assert b.homotopy.is_homogeneous()
assert b.homotopy.num_variables() == 4 # x, y, z, and the homogenizing coordinate
assert b.homotopy.num_functions() == 4 # the three rows, plus the patch
# your own systems are untouched -- the builder works on clones
assert not fixed.is_homogeneous()
The solver bridges the two for you: hand it your affine target and affine start points, as above,
and it brings both into the homotopy’s coordinates. Pass projectivize=False to build the
homotopy exactly as written, which is what you want when the affine coordinates are the point –
all-real tracking, for instance. Mixing the two, a projective system with an affine one, is
refused rather than quietly mangled.
Deform a product of linears into a polynomial¶
The regeneration “add a degree” step itself is the same construction with a different moving pair:
deform a product of linears into the actual polynomial, holding a slice fixed. Here the product
\((x-1)(x+1)\) deforms into the circle, with the slice y = 1/2 static. At \(t=1\) the
moving row is \(\gamma\,(x-1)(x+1)\), whose roots on the slice are \((\pm 1, 1/2)\):
x, y = bertini.Variable('x'), bertini.Variable('y')
fixed = bertini.System(); fixed.add_variable_group(bertini.VariableGroup([x, y]))
fixed.add_linear(np.array([[0, 1]]), np.array([x, y]), ['-1/2']) # static slice y = 1/2
start_moving = bertini.System(); start_moving.add_variable_group(bertini.VariableGroup([x, y]))
start_moving.add_products_of_linears([[[1, 0, -1], [1, 0, 1]]]) # (x-1)(x+1), a structured block
end_moving = bertini.System(); end_moving.add_variable_group(bertini.VariableGroup([x, y]))
end_moving.add_function(x*x + y*y - 1) # the polynomial
b = nag_algorithm.straight_line_homotopy(end_moving, start_moving, fixed=fixed, gamma=gamma)
start_points = [np.array([bertini.multiprec.complex_mp(str(a)), bertini.multiprec.complex_mp('0.5')])
for a in (1, -1)]
solver = bertini.HomotopySolver(b, start_points)
solver.solve()
roots = sorted((round(complex(s[0]).real, 4), round(complex(s[1]).real, 4))
for s in solver.all_solutions())
s3 = round(np.sqrt(3) / 2, 4)
assert roots == sorted([(s3, 0.5), (-s3, 0.5)]) # circle ∩ {y = 1/2}
The product-of-linears is a first-class structured block (it is never expanded into nodes), the static slice stays put, and only the one regenerating row carries \(t\). That is the unit of work a regeneration cascade repeats – and at every step the fixed equations are evaluated once.
Complete example¶
"""Move a slice, hold the system fixed -- Bertini 2 tutorial (moving_slice).
Run: python moving_slice.py
"""
import numpy as np
import bertini
from bertini import nag_algorithm
def move_one_slice(gamma):
"""A fixed unit circle sliced by a line sliding from y=0 to y-x=0."""
x, y = bertini.Variable('x'), bertini.Variable('y')
fixed = bertini.System()
fixed.add_variable_group(bertini.VariableGroup([x, y]))
fixed.add_function(x*x + y*y - 1) # the fixed unit circle
start_moving = bertini.System() # the slice at t=1: the x-axis y = 0
start_moving.add_variable_group(bertini.VariableGroup([x, y]))
start_moving.add_function(y)
end_moving = bertini.System() # the slice at t=0: the diagonal y - x = 0
end_moving.add_variable_group(bertini.VariableGroup([x, y]))
end_moving.add_function(y - x)
b = nag_algorithm.straight_line_homotopy(end_moving, start_moving, fixed=fixed, gamma=gamma)
start_points = [np.array([bertini.multiprec.complex_mp('1'), bertini.multiprec.complex_mp('0')]),
np.array([bertini.multiprec.complex_mp('-1'), bertini.multiprec.complex_mp('0')])]
solver = bertini.HomotopySolver(b, start_points) # b.target is the t=0 system, assembled for us
solver.solve()
roots = sorted((round(complex(s[0]).real, 4), round(complex(s[1]).real, 4))
for s in solver.all_solutions())
r = round(1 / np.sqrt(2), 4)
assert roots == sorted([(r, r), (-r, -r)]) # circle ∩ diagonal = (±1/√2, ±1/√2)
def static_and_moving_slice(gamma):
"""A fixed unit sphere cut by a static slice z=0 and a moving slice y=0 -> y-x=0."""
x, y, z = bertini.Variable('x'), bertini.Variable('y'), bertini.Variable('z')
fixed = bertini.System()
fixed.add_variable_group(bertini.VariableGroup([x, y, z]))
fixed.add_function(x*x + y*y + z*z - 1) # the sphere
fixed.add_linear(np.array([[0, 0, 1]]), np.array([x, y, z])) # static slice z = 0
start_moving = bertini.System(); start_moving.add_variable_group(bertini.VariableGroup([x, y, z]))
start_moving.add_function(y) # moving slice at t=1
end_moving = bertini.System(); end_moving.add_variable_group(bertini.VariableGroup([x, y, z]))
end_moving.add_function(y - x) # moving slice at t=0
b = nag_algorithm.straight_line_homotopy(end_moving, start_moving, fixed=fixed, gamma=gamma)
start_points = [np.array([bertini.multiprec.complex_mp(str(a)), bertini.multiprec.complex_mp('0'),
bertini.multiprec.complex_mp('0')]) for a in (1, -1)]
solver = bertini.HomotopySolver(b, start_points)
solver.solve()
roots = sorted((round(complex(s[0]).real, 4), round(complex(s[1]).real, 4), round(complex(s[2]).real, 4))
for s in solver.all_solutions())
r = round(1 / np.sqrt(2), 4)
assert roots == sorted([(r, r, 0.0), (-r, -r, 0.0)])
# The fixed system is left out of the motion: dH/dt is zero on the fixed blocks. Affine here,
# so the rows line up with the three equations; the builder projectivizes by default, which
# adds a homogenizing coordinate and a patch row.
Haff = nag_algorithm.straight_line_homotopy(end_moving, start_moving, fixed=fixed,
gamma=gamma, projectivize=False).homotopy
pt = np.array([bertini.multiprec.complex_mp('0.3'),
bertini.multiprec.complex_mp('0.4'),
bertini.multiprec.complex_mp('0.5')])
# H evaluates at the precision of the point it is given; nothing to align first.
dHdt = Haff.eval_time_derivative(pt, bertini.multiprec.complex_mp('0.5'))
assert abs(complex(dHdt[0])) == 0.0 # sphere row: out of dH/dt
assert abs(complex(dHdt[1])) == 0.0 # static slice row: out of dH/dt
assert abs(complex(dHdt[2])) > 0.0 # only the moving slice carries t
# projective by default: a homogenizing coordinate per affine group, plus a patch row
assert b.homotopy.is_homogeneous()
assert b.homotopy.num_variables() == 4 and b.homotopy.num_functions() == 4
assert not fixed.is_homogeneous() # your own systems are untouched
def deform_product_into_polynomial(gamma):
"""Deform a product of linears (x-1)(x+1) into the circle, slice y=1/2 static."""
x, y = bertini.Variable('x'), bertini.Variable('y')
fixed = bertini.System(); fixed.add_variable_group(bertini.VariableGroup([x, y]))
fixed.add_linear(np.array([[0, 1]]), np.array([x, y]), ['-1/2']) # static slice y = 1/2
start_moving = bertini.System(); start_moving.add_variable_group(bertini.VariableGroup([x, y]))
start_moving.add_products_of_linears([[[1, 0, -1], [1, 0, 1]]]) # (x-1)(x+1), a structured block
end_moving = bertini.System(); end_moving.add_variable_group(bertini.VariableGroup([x, y]))
end_moving.add_function(x*x + y*y - 1) # the polynomial
b = nag_algorithm.straight_line_homotopy(end_moving, start_moving, fixed=fixed, gamma=gamma)
start_points = [np.array([bertini.multiprec.complex_mp(str(a)), bertini.multiprec.complex_mp('0.5')])
for a in (1, -1)]
solver = bertini.HomotopySolver(b, start_points)
solver.solve()
roots = sorted((round(complex(s[0]).real, 4), round(complex(s[1]).real, 4))
for s in solver.all_solutions())
s3 = round(np.sqrt(3) / 2, 4)
assert roots == sorted([(s3, 0.5), (-s3, 0.5)]) # circle ∩ {y = 1/2}
def main():
gamma = bertini.coefficient(bertini.multiprec.complex_mp('0.6', '0.8')) # off the real axis
move_one_slice(gamma)
static_and_moving_slice(gamma)
deform_product_into_polynomial(gamma)
if __name__ == '__main__':
main()