🧱 Build a system in pieces with clone and concatenate¶

Sometimes a polynomial system arrives in pieces. One block of equations describes some geometry; another block adds constraints. It is natural to build each block as its own System and then stack them into one square system to solve.

concatenate() does exactly that: it appends the second system’s functions onto a copy of the first, producing a new system. The two blocks share the same unknowns, so the result is one system over those unknowns with all the functions together. This tutorial stacks two four-variable, two-function blocks into a square \(4 \times 4\) system and solves it – and shows how clone() lets you declare the variables once and reuse them for every block.

The pieces¶

Four unknowns \(x_1, x_2, x_3, x_4\). The first block is some geometry – two equations relating the unknowns – and the second block adds two constraints. Each block is a perfectly good System on its own; neither is square (two equations, four unknowns), so neither can be solved alone.

x1, x2, x3, x4 = (bertini.Variable(n) for n in ('x1', 'x2', 'x3', 'x4'))
group = bertini.VariableGroup([x1, x2, x3, x4])

geometry = bertini.System()
geometry.add_variable_group(group)
geometry.add_function(x1*x1 - x2)        # x2 = x1²
geometry.add_function(x2*x2 - x3)        # x3 = x2²

constraints = bertini.System()
constraints.add_variable_group(group)
constraints.add_function(x3 - x4)        # x4 = x3
constraints.add_function(x4 - x1)        # x1 = x4

assert geometry.num_functions() == 2 and geometry.num_variables() == 4
assert constraints.num_functions() == 2 and constraints.num_variables() == 4

Both blocks are written over the same variables. That is the only requirement for concatenation: the two systems must share their variable ordering. Reusing the same group (and the same Variable objects) guarantees it – and even rebuilding the variables by name would, since variables in Bertini are canonical by name.

Stack them and solve¶

concatenate() returns a new square system; the two blocks are left untouched.

system = bertini.system.concatenate(geometry, constraints)

assert system.num_functions() == 4 and system.num_variables() == 4
assert geometry.num_functions() == 2        # inputs unchanged
assert constraints.num_functions() == 2
assert system.degrees() == [2, 2, 1, 1]

The combined system is square, so a total-degree start system can solve it. Its path count is the product of the degrees, \(2 \times 2 \times 1 \times 1 = 4\).

bertini.random.set_random_seed(1)           # reproducible start system + gamma
zd = bertini.ZeroDimSolver(system, endgame='cauchy', mptype='adaptive', startsystem='binomial')
zd.solve()

solutions = [np.array([complex(c) for c in s]) for s in zd.all_solutions()]
assert len(solutions) == 4

# every computed point satisfies the concatenated system
assert all(max(abs(complex(v)) for v in system.eval(s)) < 1e-7 for s in solutions)

By hand: the four equations force \(x_4 = x_1\), \(x_3 = x_4 = x_1\), \(x_2 = x_1^2\), and \(x_2^2 = x_3 = x_1\), so \(x_1^4 = x_1\). That gives \(x_1 \in \{0, 1, \omega, \omega^2\}\) (the cube roots of unity, and zero) – four solutions, of which two are real.

real = [s for s in solutions if max(abs(c.imag) for c in s) < 1e-7]
found = sorted(tuple(round(c.real, 4) for c in s) for s in real)
assert found == [(0.0, 0.0, 0.0, 0.0), (1.0, 1.0, 1.0, 1.0)]

Declare the variables once, with clone¶

Above, each block re-declared the variable group. When you are assembling several blocks, it is tidier to set up the variables once and clone() that setup for each block.

A clone shares the original’s variables (and gets its own evaluation memory), so blocks cloned from a common base are guaranteed to share variable ordering – exactly what concatenation needs. Adding functions to one clone does not touch the others.

base = bertini.System()                     # variable setup only -- no functions yet
base.add_variable_group(bertini.VariableGroup([x1, x2, x3, x4]))

geometry = bertini.system.clone(base)
geometry.add_function(x1*x1 - x2)
geometry.add_function(x2*x2 - x3)

constraints = bertini.system.clone(base)
constraints.add_function(x3 - x4)
constraints.add_function(x4 - x1)

system = bertini.system.concatenate(geometry, constraints)
assert system.num_functions() == 4

bertini.random.set_random_seed(1)
zd = bertini.ZeroDimSolver(system, endgame='cauchy', mptype='adaptive', startsystem='binomial')
zd.solve()
assert len(zd.all_solutions()) == 4

Same four solutions, assembled from cloned blocks. The pattern – set up variables once, clone per block, concatenate, solve – scales to as many blocks as your problem comes in.

Note

What makes the blocks line up is that they are over the same variables, and a variable is its name: the clone’s \(x_1\) is the base’s \(x_1\). The same holds for a system copied with copy.deepcopy or sent through pickle (to another process, say), so those combine with the original too. clone is the direct way, since it has no archive to go through.

Two systems are combined only when they are over the same variables, in the same groups and the same order. A point is a list of numbers, and the order of the variables is what says which number is which, so systems over \((x, y)\) and \((y, x)\) are refused, not quietly matched up by position:

u, v = bertini.Variable('u'), bertini.Variable('v')

one_way = bertini.System()
one_way.add_variable_group(bertini.VariableGroup([u, v]))
one_way.add_function(u - 2)

other_way = bertini.System()
other_way.add_variable_group(bertini.VariableGroup([v, u]))
other_way.add_function(v - 3)

try:
    bertini.system.concatenate(one_way, other_way)
except RuntimeError as refusal:
    print('same variables in a different order' in str(refusal))

True

A system that declares no variables at all takes the variables of the one it is combined with.

Complete example¶

concatenate_systems.py¶

"""Build a system in pieces with clone and concatenate (Bertini 2 tutorial).

Stacks two four-variable, two-function blocks into a square 4x4 system and solves it.
Run:  python concatenate_systems.py
"""

import numpy as np
import bertini


def build_pieces():
    """Build two non-square blocks (geometry + constraints) over shared variables."""
    x1, x2, x3, x4 = (bertini.Variable(n) for n in ('x1', 'x2', 'x3', 'x4'))
    group = bertini.VariableGroup([x1, x2, x3, x4])

    geometry = bertini.System()
    geometry.add_variable_group(group)
    geometry.add_function(x1*x1 - x2)        # x2 = x1²
    geometry.add_function(x2*x2 - x3)        # x3 = x2²

    constraints = bertini.System()
    constraints.add_variable_group(group)
    constraints.add_function(x3 - x4)        # x4 = x3
    constraints.add_function(x4 - x1)        # x1 = x4

    assert geometry.num_functions() == 2 and geometry.num_variables() == 4
    assert constraints.num_functions() == 2 and constraints.num_variables() == 4

    return (x1, x2, x3, x4), geometry, constraints


def stack(geometry, constraints):
    """Concatenate the two blocks into a new square system; inputs unchanged."""
    system = bertini.system.concatenate(geometry, constraints)

    assert system.num_functions() == 4 and system.num_variables() == 4
    assert geometry.num_functions() == 2        # inputs unchanged
    assert constraints.num_functions() == 2
    assert system.degrees() == [2, 2, 1, 1]

    return system


def solve_and_check(system):
    """Solve the square system and verify the four solutions."""
    bertini.random.set_random_seed(1)           # reproducible start system + gamma
    zd = bertini.ZeroDimSolver(system, endgame='cauchy', mptype='adaptive', startsystem='binomial')
    zd.solve()

    solutions = [np.array([complex(c) for c in s]) for s in zd.all_solutions()]
    assert len(solutions) == 4

    # every computed point satisfies the concatenated system
    assert all(max(abs(complex(v)) for v in system.eval(s)) < 1e-7 for s in solutions)

    # the two real solutions are the zero point and the all-ones point
    real = [s for s in solutions if max(abs(c.imag) for c in s) < 1e-7]
    found = sorted(tuple(round(c.real, 4) for c in s) for s in real)
    assert found == [(0.0, 0.0, 0.0, 0.0), (1.0, 1.0, 1.0, 1.0)]

    return solutions


def build_with_clone(variables):
    """Declare variables once, clone per block, concatenate, and solve."""
    x1, x2, x3, x4 = variables

    base = bertini.System()                     # variable setup only -- no functions yet
    base.add_variable_group(bertini.VariableGroup([x1, x2, x3, x4]))

    geometry = bertini.system.clone(base)
    geometry.add_function(x1*x1 - x2)
    geometry.add_function(x2*x2 - x3)

    constraints = bertini.system.clone(base)
    constraints.add_function(x3 - x4)
    constraints.add_function(x4 - x1)

    system = bertini.system.concatenate(geometry, constraints)
    assert system.num_functions() == 4

    bertini.random.set_random_seed(1)
    zd = bertini.ZeroDimSolver(system, endgame='cauchy', mptype='adaptive', startsystem='binomial')
    zd.solve()
    assert len(zd.all_solutions()) == 4


def main():
    variables, geometry, constraints = build_pieces()
    system = stack(geometry, constraints)
    solve_and_check(system)
    build_with_clone(variables)


if __name__ == '__main__':
    main()