Code Snippets

The functions are based on the code snippets of the menus of the Jupyter notebook extension of PHCpy. Every function is self contained and illustrates one particular feature.

the blackbox solver

def solve_random_trinomials():
    """
    Illustrates the solution of random trinomials.
    """
    from phcpy.solver import random_trinomials
    f = random_trinomials()
    for pol in f: print(pol)
    from phcpy.solver import solve
    sols = solve(f)
    for sol in sols: print(sol)
    print(len(sols), "solutions found")

def solve_specific_trinomials():
    """
    Illustrates the solution of specific trinomials.
    """
    f = ['x^2*y^2 + 2*x - 1;', 'x^2*y^2 - 3*y + 1;']
    from phcpy.solver import solve
    sols = solve(f)
    for sol in sols: print(sol)

def solution_from_string_to_dictionary():
    """
    Illustrates the dictionary format of a solution.
    """
    p = ['x + y - 1;', '2*x - 3*y + 1;']
    from phcpy.solver import solve
    sols = solve(p)
    print(sols[0])
    from phcpy.solutions import strsol2dict
    dsol = strsol2dict(sols[0])
    print(dsol.keys())
    for key in dsol.keys(): print('the value for', key, 'is', dsol[key])

def verify_by_evaluation():
    """
    Illustrates the verification of solutions by evaluation.
    """
    p = ['x + y - 1;', '2*x - 3*y + 1;']
    from phcpy.solver import solve
    sols = solve(p)
    from phcpy.solutions import strsol2dict, evaluate
    dsol = strsol2dict(sols[0])
    eva = evaluate(p, dsol)
    for val in eva: print(val)

def make_a_solution():
    """
    Illustrates the making of a solution.
    """
    from phcpy.solutions import make_solution
    s0 = make_solution(['x', 'y'], [float(3.14), complex(0, 2.7)])
    print(s0)
    s1 = make_solution(['x', 'y'], [int(2), int(3)])
    print(s1)

def filter_solution_lists():
    """
    Illustrates the filtering of solution lists.
    """
    from phcpy.solutions import make_solution, is_real, filter_real
    s0 = make_solution(['x', 'y'], [float(3.14), complex(0, 2.7)])
    print(is_real(s0, 1.0e-8))
    s1 = make_solution(['x', 'y'], [int(2), int(3)])
    print(is_real(s1, 1.0e-8))
    realsols = filter_real([s0, s1], 1.0e-8, 'select')
    for sol in realsols: print(sol)

def coordinates_names_values():
    """
    Illustrates coordinates, names, values of a solution.
    """
    from phcpy.solver import solve
    p = ['x^2*y^2 + x + 1;', 'x^2*y^2 + y + 1;']
    s = solve(p)
    print(s[0])
    from phcpy.solutions import coordinates, make_solution
    (names, values) = coordinates(s[0])
    print(names)
    print(values)
    s0 = make_solution(names, values)
    print(s0)

def fix_retrieve_seed():
    """
    Illustrates fixing and retrieving the seed.
    """
    from phcpy.dimension import set_seed, get_seed
    set_seed(2024)
    print(get_seed())

def solve_with_four_threads():
    """
    Ilustrates solving with four threads.
    """
    from phcpy.solver import solve
    from phcpy.families import cyclic
    nbrt = 4 # number of tasks
    pols = cyclic(6)
    print('solving the cyclic 6-roots problem :')
    for pol in pols: print(pol)
    from datetime import datetime
    starttime = datetime.now()
    sols = solve(pols)
    stoptime = datetime.now()
    elapsed = stoptime - starttime
    print(f'elapsed time with no multithreading :')
    print(elapsed)
    starttime = datetime.now()
    sols = solve(pols, tasks=nbrt)
    stoptime = datetime.now()
    elapsed = stoptime - starttime
    print(f'elapsed time with {nbrt} threads :')
    print(elapsed)

def four_root_counts():
    """
    Illustrates four root counts.
    """
    f = ['x^3*y^2 + x*y^2 + x^2;', 'x^5 + x^2*y^3 + y^2;']
    from phcpy.starters import total_degree
    print('the total degree :', total_degree(f))
    from phcpy.starters import m_homogeneous_bezout_number as mbz
    (bz, part) = mbz(f)
    print('a multihomogeneous Bezout number :', bz)
    from phcpy.starters import linear_product_root_count as lrc
    (rc, ssrc) = lrc(f)
    print('a linear-product root count :', rc)
    from phcpy.volumes import stable_mixed_volume
    (mv, smv) = stable_mixed_volume(f)
    print('the mixed volume :', mv)
    print('the stable mixed volume :', smv)

def newton_and_deflation():
    """
    Illustrates Newton's method and deflation.
    """
    p = ['(29/16)*x^3 - 2*x*y;', 'x^2 - y;']
    from phcpy.solutions import make_solution
    s = make_solution(['x', 'y'],[float(1.0e-6), float(1.0e-6)])
    print(s)
    from phcpy.deflation import double_newton_step
    s2 = double_newton_step(p, [s])
    print(s2[0])
    s3 = double_newton_step(p, s2)
    print(s3[0])
    from phcpy.deflation import double_deflate
    sd = double_deflate(p, [s])
    print(sd[0])

def overconstrained_deflation():
    """
    Illustrates deflation of an overconstrained system.
    """
    from phcpy.solutions import make_solution
    from phcpy.deflation import double_deflate
    sol = make_solution(['x', 'y'], [float(1.0e-6), float(1.0e-6)])
    print(sol)
    pols = ['x**2;', 'x*y;', 'y**2;']
    sols = double_deflate(pols, [sol], tolrnk=1.0e-8)
    print(sols[0])
    sols = double_deflate(pols, [sol], tolrnk=1.0e-4)
    print(sols[0])

def equation_and_variable_scaling():
    """
    Illustrates equation and variable scaling.
    """
    print('solving without scaling ...')
    from phcpy.solver import solve
    p = ['0.000001*x^2 + 0.000004*y^2 - 4;', '0.000002*y^2 - 0.001*x;']
    psols = solve(p)
    for sol in psols: print(sol)
    print('solving after scaling ...')
    from phcpy.scaling import double_scale_system as scalesys
    from phcpy.scaling import double_scale_solutions as scalesols
    (q, c) = scalesys(p)
    print('the scaled polynomial system :')
    for pol in q: print(pol)
    qsols = solve(q)
    ssols = scalesols(len(q), qsols, c)
    for sol in ssols: print(sol)

path trackers

def total_degree_start_system():
    """
    The product of the degrees of the polynomials (the total degree)
    provide an upper bound on the number of solutions.
    A total degree start system is a simple system that has
    as many solutions as the product of the degrees.
    """
    from phcpy.starters import total_degree
    from phcpy.starters import total_degree_start_system
    from phcpy.trackers import double_track
    p = ['x^2 + 4*y^2 - 4;', '2*y^2 - x;']
    d = total_degree(p)
    print('the total degree :', d)
    (q, qsols) = total_degree_start_system(p)
    print('the number of start solutions :', len(qsols))
    print('the start system :', q)
    s = double_track(p, q, qsols)
    print('the number of solutions :', len(s))
    for sol in s: print(sol)

def track_one_path():
    """
    Illustrates the tracking of one solution path.
    The order in which the solutions appear at the end
    depends on the gamma constant.
    """
    from phcpy.starters import total_degree_start_system
    from phcpy.trackers import double_track
    p = ['x^2 + 4*y^2 - 4;', '2*y^2 - x;']
    (q, qsols) = total_degree_start_system(p)
    g1, s1 = double_track(p, q, [qsols[2]])
    print('first gamma :', g1)
    print(s1[0])
    g2, s2 = double_track(p, q, [qsols[2]])
    print('second gamma :', g2)
    print(s2[0])

def track_with_fixed_gamma():
    """
    Fixing the gamma in the homotopies fixes the order
    in which the solutions appear at the end of the paths.
    """
    from phcpy.starters import total_degree_start_system
    from phcpy.trackers import double_track
    p = ['x^2 + 4*y^2 - 4;', '2*y^2 - x;']
    (q, qsols) = total_degree_start_system(p)
    g3, s3 = double_track(p, q, [qsols[2]], \
        gamma=complex(0.824372806319,0.56604723848934))
    print('gamma :', g3)
    print('the solution at the end:')
    print(s3[0])

def get_next_point_on_path():
    """
    A step-by-step path tracker gives control to the user
    who can ask for the next point on a path.
    """
    from phcpy.starters import total_degree_start_system
    p = ['x**2 + 4*x**2 - 4;', '2*y**2 - x;']
    (q, s) = total_degree_start_system(p)
    from phcpy.trackers import initialize_double_tracker
    from phcpy.trackers import initialize_double_solution
    from phcpy.trackers import next_double_solution
    initialize_double_tracker(p, q)
    initialize_double_solution(len(p), s[0])
    s1 = next_double_solution()
    print('the next point on the solution path :')
    print(s1)
    print(next_double_solution())
    print(next_double_solution())
    initialize_double_solution(len(p), s[1])
    points = [next_double_solution() for i in range(11)]
    from phcpy.solutions import strsol2dict
    dicpts = [strsol2dict(sol) for sol in points]
    xvals = [sol['x'] for sol in dicpts]
    print('the x-coordinates on the path :')
    for x in xvals: print(x)

def plot_trajectories():
    """
    The step-by-step path tracker is applied to plot the
    trajectories of the solutions using matplotlib.
    """
    import matplotlib.pyplot as plt
    p = ['x^2 + y - 3;', 'x + 0.125*y^2 - 1.5;']
    print('constructing a total degree start system ...')
    from phcpy.starters import total_degree_start_system
    q, qsols = total_degree_start_system(p)
    print('number of start solutions :', len(qsols))
    from phcpy.trackers import initialize_double_tracker
    from phcpy.trackers import initialize_double_solution
    from phcpy.trackers import next_double_solution
    initialize_double_tracker(p, q, False)
    from phcpy.solutions import strsol2dict
    plt.ion()
    fig = plt.figure()
    for k in range(len(qsols)):
        if(k == 0):
           axs = fig.add_subplot(221)
        elif(k == 1):
           axs = fig.add_subplot(222)
        elif(k == 2):
            axs = fig.add_subplot(223)
        elif(k == 3):
           axs = fig.add_subplot(224)
        startsol = qsols[k]
        initialize_double_solution(len(p),startsol)
        dictsol = strsol2dict(startsol)
        xpoints =  [dictsol['x']]
        ypoints =  [dictsol['y']]
        for k in range(300):
            ns = next_double_solution()
            dictsol = strsol2dict(ns)
            xpoints.append(dictsol['x'])
            ypoints.append(dictsol['y'])
            tval = dictsol['t'].real
            if(tval >= 1.0):
                break
        print(ns)
        xre = [point.real for point in xpoints]
        yre = [point.real for point in ypoints]
        axs.set_xlim(min(xre)-0.3, max(xre)+0.3)
        axs.set_ylim(min(yre)-0.3, max(yre)+0.3)
        dots, = axs.plot(xre,yre,'r-')
        fig.canvas.draw()
    fig.canvas.draw()
    ans = input('hit return to continue')

def polyhedral_homotopies():
    """
    Polyhedral homotopies solve random coefficient systems
    tracking an optimal number of solution paths,
    that is equal to the mixed volume.
    """
    from phcpy.volumes import mixed_volume
    from phcpy.volumes import double_polyhedral_homotopies
    from phcpy.trackers import double_track
    p = ['x^3*y^2 - 3*x^3 + 7;','x*y^3 + 6*y^3 - 9;']
    print('the mixed volume :', mixed_volume(p))
    (q, qsols) = double_polyhedral_homotopies()
    print('the number of start solutions :', len(qsols))
    gamma, psols = double_track(p, q, qsols)
    print('the number of solutions at the end :', len(psols))
    for sol in psols: print(sol)

sweep homotopies

def quadratic_turning_point():
    """
    Arc length parameter continuation is applied in a real sweep
    which ends at a quadratic turning point.
    """
    from phcpy.sweepers import double_real_sweep
    from phcpy.solutions import make_solution
    circle = ['x^2 + y^2 - 1;', 'y*(1-s) + (y-2)*s;']
    first = make_solution(['x', 'y', 's'], [1.0, 0.0, 0.0])
    second = make_solution(['x', 'y', 's'], [-1.0, 0.0, 0.0])
    startsols = [first, second]
    newsols = double_real_sweep(circle, startsols)
    for sol in newsols: print(sol)

def complex_parameter_homotopy_continuation():
    """
    Sweeps the parameter space with a convex linear combination
    of the parameters.  By a random gamma constant, no singularities
    are encountered during this complex sweep.
    """
    from phcpy.sweepers import double_complex_sweep
    from phcpy.solutions import make_solution
    circle = ['x^2 + y^2 - 1;']
    first = make_solution(['x', 'y'], [1.0, 0.0])
    second = make_solution(['x', 'y'], [-1.0, 0.0])
    startsols = [first, second]
    par = ['y']
    start = [0, 0]
    target = [2, 0]
    newsols = double_complex_sweep(circle, startsols, 2, par, start, target)
    for sol in newsols: print(sol)

Schubert calculus

def lines_meeting_four_lines():
    """
    Applies Pieri homotopies to compute all lines meeting
    four given lines in 3-space.
    """
    from phcpy.schubert import pieri_root_count
    from phcpy.schubert import random_complex_matrix, run_pieri_homotopies
    (m, p, q) = (2, 2, 0)
    dim = m*p + q*(m+p)
    roco = pieri_root_count(m, p, q)
    print('the root count :', roco)
    L = [random_complex_matrix(m+p, m) for _ in range(dim)]
    (f, fsols) = run_pieri_homotopies(m, p, q, L)
    for sol in fsols: print(sol)
    print('number of solutions :', len(fsols))

def line_producing_interpolating_curves():
    """
    Applies Pieri homotopies to compute line producing curves,
    interpolating at given lines in 3-space.
    """
    from phcpy.schubert import pieri_root_count
    from phcpy.schubert import random_complex_matrix, run_pieri_homotopies
    (m, p, q) = (2, 2, 1)
    dim = m*p + q*(m+p)
    roco = pieri_root_count(m, p, q)
    print('the root count :', roco)
    L = [random_complex_matrix(m+p, m) for _ in range(dim)]
    points = random_complex_matrix(dim, 1)
    (f, fsols) = run_pieri_homotopies(m, p, q, L, 0, points)
    print('number of solutions :', len(fsols))

def resolve_some_schubert_conditions():
    """
    Resolves an example of a Schubert condition,
    applying the Littlewood-Richardson rule.
    """
    from phcpy.schubert import resolve_schubert_conditions
    brackets = [[2, 4, 6], [2, 4, 6], [2, 4, 6]]
    roco = resolve_schubert_conditions(6, 3, brackets)
    print('number of solutions :', roco)

def solve_generic_schubert_problem():
    """
    Runs the Littlewood-Richardson homotopies to solve
    a generic instance of a Schubert problem.
    """
    brackets = [[2, 4, 6], [2, 4, 6], [2, 4, 6]]
    from phcpy.schubert import double_littlewood_richardson_homotopies as lrh
    (count, flags, sys, sols) = lrh(6, 3, brackets, verbose=False)
    print('the root count :', count)
    for sol in sols: print(sol)
    print('the number of solutions :', len(sols))

power series expansions

def example4pade():
    """
    The function f(z) = ((1 + 1/2*z)/(1 + 2*z))^(1/2) is
    a solution x(s) of (1-s)*(x^2 - 1) + s*(3*x^2 - 3/2) = 0
    """
    from phcpy.solutions import make_solution
    from phcpy.series import double_newton_at_point
    from phcpy.series import double_pade_approximants
    pol = ['(x^2 - 1)*(1-s) + (3*x^2 - 3/2)*s;']
    variables = ['x', 's']
    sol1 = make_solution(variables, [1, 0])
    sol2 = make_solution(variables, [-1, 0])
    sols = [sol1, sol2]
    print('start solutions :')
    for sol in sols: print(sol)
    srs = double_newton_at_point(pol, sols, idx=2)
    print('The series :')
    for ser in srs: print(ser)
    pad = double_pade_approximants(pol, sols, idx=2)
    print('the Pade approximants :')
    for app in pad: print(app)

def viviani_expansion(vrblvl=0):
    """
    Computes the power series expansion for the Viviani curve,
    from a natural parameter perspective, in double precision.
    """
    from phcpy.series import double_newton_at_series
    pols = [ '2*t^2 - x;', \
             'x^2 + y^2 + z^2 - 4;' , \
             '(x-1)^2 + y^2 - 1;']
    lser = [ '2*t^2;', '2*t;', '2;']
    nser = double_newton_at_series(pols, lser, maxdeg=12, nbr=8)
    variables = ['x', 'y', 'z']
    for (var, pol) in zip(variables, nser): print(var, '=', pol)

def apollonius_expansions():
    """
    Compare the series expansions at two solutions
    for the problem of Apollonius.
    """
    from phcpy.series import double_newton_at_series
    pols = [ 'x1^2 + 3*x2^2 - r^2 - 2*r - 1;', \
             'x1^2 + 3*x2^2 - r^2 - 4*x1 - 2*r + 3;', \
    '3*t^2 + x1^2 - 6*t*x2 + 3*x2^2 - r^2 + 6*t - 2*x1 - 6*x2 + 2*r + 3;']
    lser1 = ['1;', '1 + 0.536*t;', '1 + 0.904*t;']
    lser2 = ['1;', '1 + 7.464*t;', '1 + 11.196*t;']
    nser1 = double_newton_at_series(pols, lser1, idx=4, nbr=7)
    nser2 = double_newton_at_series(pols, lser2, idx=4, nbr=7)
    variables = ['x', 'y', 'z']
    print('the first solution series :')
    for (var, pol) in zip(variables, nser1): print(var, '=', pol)
    print('the second solution series :')
    for (var, pol) in zip(variables, nser2): print(var, '=', pol)

positive dimensional solution sets

def twisted_cubic_witness_set():
    """
    Makes a witness set for the twisted cubic.
    """
    twisted = ['x^2 - y;', 'x^3 - z;']
    from phcpy.sets import double_embed
    embtwist = double_embed(3, 1, twisted)
    print('embedded system augmented with a random hyperplane :')
    for pol in embtwist:
        print(pol)
    from phcpy.solver import solve
    sols = solve(embtwist)
    print('number of generic points :', len(sols))
    for sol in sols: print(sol)

def homotopy_membership_test():
    """
    Homotopies to decide whether a point belongs to a positive
    dimensional solution set.
    """
    from phcpy.families import cyclic
    c4 = cyclic(4)
    from phcpy.sets import double_embed
    c4e1 = double_embed(4, 1, c4)
    from phcpy.solver import solve
    sols = solve(c4e1)
    from phcpy.solutions import filter_zero_coordinates as filter
    genpts = filter(sols, 'zz1', 1.0e-8, 'select')
    print('generic points on the cyclic 4-roots set :')
    for sol in genpts: print(sol)
    from phcpy.sets import double_membertest
    pt0 = [1, 0, -1, 0, 1, 0, -1, 0]
    ismbr = double_membertest(c4e1, sols, 1, pt0)
    print('Is', pt0, 'a member?', ismbr)
    pt1 = [1, 0, 1, 0, -1, 0, -1, 0]
    ismbr = double_membertest(c4e1, sols, 1, pt1)
    print('Is', pt1, 'a member?', ismbr)

def factor_a_cubic():
    """
    Illustrates the numerical factorization of a cubic.
    """
    cubic = '(x+1)*(x^2 + y^2 + 1);'
    from phcpy.sets import double_hypersurface_set
    (wit, pts) = double_hypersurface_set(2, cubic)
    for pol in wit:
        print(pol)
    print('number of witness points :', len(pts))
    for (idx, sol) in enumerate(pts):
        print('Solution', idx+1, ':')
        print(sol)
    from phcpy.factor import double_monodromy_breakup, write_factorization
    deco = double_monodromy_breakup(wit, pts, dim=1)
    write_factorization(deco)

def cascades_of_homotopies():
    """
    A cascade of homotopies computes generic points on all positive
    components of the solution set, for all dimensions.
    """
    pol1 = '(x^2 + y^2 + z^2 - 1)*(y - x^2)*(x - 0.5);'
    pol2 = '(x^2 + y^2 + z^2 - 1)*(z - x^3)*(y - 0.5);'
    pol3 = '(x^2 + y^2 + z^2 - 1)*(z - x*y)*(z - 0.5);'
    pols = [pol1, pol2, pol3]
    from phcpy.cascades import double_top_cascade, double_cascade_filter
    (embpols, sols0, sols1) = double_top_cascade(3, 2, pols)
    print('at dimension 2, degree :', len(sols0))
    (wp1, ws0, ws1) = double_cascade_filter(2, embpols, sols1, tol=1.0e-8)
    print('at dimension 1, candidate generic points :', len(ws0))
    (wp0, ws0, ws1) = double_cascade_filter(1, wp1, ws1, tol=1.0e-8)
    print('candidate isolated points :', len(ws0))

def numerical_irreducible_decomposition():
    """
    Runs a test example on a numerical irreducible decomposition.
    """
    pol0 = '(x1-1)*(x1-2)*(x1-3)*(x1-4);'
    pol1 = '(x1-1)*(x2-1)*(x2-2)*(x2-3);'
    pol2 = '(x1-1)*(x1-2)*(x3-1)*(x3-2);'
    pol3 = '(x1-1)*(x2-1)*(x3-1)*(x4-1);'
    pols = [pol0, pol1, pol2, pol3]
    from phcpy.decomposition import solve, write_decomposition
    deco = solve(pols, verbose=False)
    write_decomposition(deco)

def sphere_cylinder_intersection():
    """
    A sphere intersected by a cylinder defines a quartic curve.
    """
    sphere = 'X^2 + Y^2 + Z^2 - 1;'
    cylinder = 'X^2 + 1.0e-14*Y^2 + (Z - 0.5)^2 - 1;'
    from phcpy.sets import double_hypersurface_set
    (spheqs, sphpts) = double_hypersurface_set(3, sphere)
    (cyleqs, cylpts) = double_hypersurface_set(3, cylinder)
    from phcpy.diagonal import double_diagonal_solve
    quaeqs, quapts = double_diagonal_solve\
        (3, 2, spheqs, sphpts, 2, cyleqs, cylpts)
    for pol in quaeqs: print(pol)
    for sol in quapts: print(sol)