๐Ÿงฑ 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 list(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

clone is the right tool here precisely because it shares the variables: the cloneโ€™s \(x_1\) is the baseโ€™s \(x_1\), so the blocks line up. If you instead need a fully independent serialized copy (for example to send a system to another process), use copy.deepcopy or pickle.

Complete exampleยถ

The whole tutorial as one runnable script โ€“ assemble nothing, just run it:

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 list(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()