# Tracking Paths Step by Step

One main benefit of phcpy is the interactivity and the inversion of control. Instead of the path tracker deciding the pace of computing points on the path, the user can ask for the next point on the path, which is convenient to plot the points on the path.

## plotting solutions paths

In the example below, the real parts of the solution paths are plotted, in Fig. 15.

    from phcpy.dimension import set_seed, get_seed
    from phcpy.solutions import strsol2dict
    from phcpy.starters import total_degree_start_system
    from phcpy.trackers import initialize_double_tracker
    from phcpy.trackers import initialize_double_solution
    from phcpy.trackers import next_double_solution

The construction of the homotopy depends on the generation of random numbers. To obtain consistently the same plots, the seed of the random number generator is fixed.

    set_seed(12871)
    print('the seed :', get_seed())

What is printed is `the seed : 12871`.

The system that will be solved is defined as follows:

    p = ['x^2 + y - 3;', 'x + 0.125*y^2 - 1.5;']

For this intersection of two quadrics, we construct a total degree start system and compute four start solutions:

    q, qsols = total_degree_start_system(p)
    for pol in q:
        print(pol)

Then here is the start system:

    x^2 - 1;
    y^2 - 1;

    print('number of start solutions :', len(qsols))

and `4` is printed.

    initialize_double_tracker(p, q, False)

    plt.ion()

The code

    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-')
        dots, = axs.plot(xre,yre,'ro')
        fig.canvas.draw()
    fig.canvas.draw()

prints the solutions at the end of the paths

    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x :  1.00000071115204E+00  -4.02767598012276E-06
     y :  1.99999857771171E+00   8.05535854116567E-06
    = err :  2.614E-06 = rco :  1.000E+00 = res :  9.361E-13 =
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x :  9.99994847664516E-01  -1.64638311049799E-06
     y :  2.00001030465627E+00   3.29275198635159E-06
    == err :  7.006E-06 = rco :  1.000E+00 = res :  1.187E-11 =
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x :  9.99999462998807E-01   3.01282709881038E-06
     y :  2.00000107400770E+00  -6.02565521719758E-06
    == err :  6.245E-06 = rco :  1.000E+00 = res :  7.729E-12 =
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x : -3.00000000000000E+00   1.39180902310149E-15
     y : -6.00000000000001E+00   5.55716009286787E-15
    == err :  4.185E-07 = rco :  1.000E+00 = res :  2.071E-14 =

and shows then the plot

Fig. 15 The real parts of four solution paths.

Typically, with an adaptive step size control, the points are closer to each other at the start and end of the paths, and where the paths turn. The `trackers` module exports the original path trackers, which use aposteriori step size control. An aposteriori step size control algorithm determines the step size based on the performance of the corrector.

## plotting paths and poles

The apriori step size control determines the step size based on the location of the nearest pole and the curvature. In addition to the (real parts of the paths), the location of the nearest poles is plotted.

    import matplotlib.pyplot as plt

    from phcpy.dimension import set_seed, get_seed
    from phcpy.solutions import strsol2dict
    from phcpy.starters import total_degree_start_system
    from phcpy.curves import set_default_parameters, write_parameters
    from phcpy.curves import initialize_double_artificial_homotopy
    from phcpy.curves import set_double_solution, get_double_solution
    from phcpy.curves import double_predict_correct
    from phcpy.curves import double_t_value, double_closest_pole

The seed is fixed to obtain the same plots in each run.

    set_seed(12871)

The system that will be solved is defined as

    p = ['x^2 + y - 3;', 'x + 0.125*y^2 - 1.5;']

and a start system based on the total degree is constructed:

    x^2 - 1;
    y^2 - 1;

as the output of the code

    q, qsols = total_degree_start_system(p)
    for pol in q:
        print(pol)

The list qsols has the four solutions of the start system q.

Before lauching the path trackers, the parameters must be set. Default values are used

    set_default_parameters()
    write_parameters()

and shown below:

    Values of the HOMOTOPY CONTINUATION PARAMETERS :
     1. gamma : (-0.8063005962200716-0.5915060004219376j)
     2. degree of numerator of Pade approximant    : 5
     3. degree of denominator of Pade approximant  : 1
     4. maximum step size                          : 0.1
     5. minimum step size                          : 1e-06
     6. multiplication factor for the pole radius  : 0.5
     7. multiplication factor for the curvature    : 0.005
     8. tolerance on the residual of the predictor : 0.001
     9. tolerance on the residual of the corrector : 1e-08
    10. tolerance on zero series coefficients      : 1e-12
    11. maximum number of corrector steps          : 4
    12. maximum steps on a path                    : 1000

Then the homotopy is constructed:

    initialize_double_artificial_homotopy(p, q, False)

Then the code below

    plt.ion()
    fig1 = plt.figure()
    allpoles = []
    for k in range(len(qsols)):
        if(k == 0):
            axs = fig1.add_subplot(221)
        elif(k == 1):
            axs = fig1.add_subplot(222)
        elif(k == 2):
            axs = fig1.add_subplot(223)
        elif(k == 3):
            axs = fig1.add_subplot(224)
        startsol = qsols[k]
        set_double_solution(len(p), startsol)
        dictsol = strsol2dict(startsol)
        xpoints =  [dictsol['x']]
        ypoints =  [dictsol['y']]
        poles = []
        for k in range(100):
            ns = get_double_solution()
            dictsol = strsol2dict(ns)
            xpoints.append(dictsol['x'])
            ypoints.append(dictsol['y'])
            tval = dictsol['t'].real
            if(tval == 1.0):
                break
            double_predict_correct()
            pole = double_closest_pole()
            tval = double_t_value()
            locp = (tval+pole[0], pole[1])
            poles.append(locp)
        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,'b-')
        dots, = axs.plot(xre,yre,'bo')
        fig1.canvas.draw()
        allpoles.append(poles)
    fig1.canvas.draw()

prints the solutions at the end of the path:

    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x :  1.00000458683462E+00   7.95321001653953E-06
     y :  1.99999082635867E+00  -1.59064683422953E-05
    == err :  1.458E-05 = rco :  3.155E-12 = res :  1.529E-12 =
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x :  9.99999434968695E-01   9.78667863138325E-07
     y :  2.00000113004483E+00  -1.95736651949408E-06
    == err :  1.658E-05 = rco :  6.664E-12 = res :  1.976E-12 =
    t :  1.00000000000000E+00   0.00000000000000E+00\
    m : 1
    the solution for t :
     x :  1.00000600449478E+00   2.32224203692682E-07
     y :  1.99998799104108E+00  -4.64448075210226E-07
    == err :  1.671E-05 = rco :  9.377E-13 = res :  3.579E-12 =
    t :  1.00000000000000E+00   0.00000000000000E+00
    m : 1
    the solution for t :
     x : -3.00000000000000E+00   0.00000000000000E+00
     y : -6.00000000000000E+00   0.00000000000000E+00
    == err :  5.551E-16 = rco :  1.965E-01 = res :  0.000E+00 =

and produces the plot in Fig. 16.

Fig. 16 The real parts of four solution paths.

The poles are shown in Fig. 17 plotted with the code in

    fig2 = plt.figure()
    for k in range(len(qsols)):
        if(k == 0):
            axs = fig2.add_subplot(221)
        elif(k == 1):
            axs = fig2.add_subplot(222)
        elif(k == 2):
            axs = fig2.add_subplot(223)
        elif(k == 3):
            axs = fig2.add_subplot(224)
        poles = allpoles[k]
        pl0 = [pole[0] for pole in poles]
        pl1 = [pole[1] for pole in poles]
        axs.set_xlim(-0.2, 1.2)
        axs.set_ylim(-0.5, 0.5)
        dots, = axs.plot(pl0,pl1,'r+')
        fig2.canvas.draw()
    fig2.canvas.draw()

Fig. 17 The nearest poles of four solution paths.

Observe that for this example, more poles are located closer to the middle and the end of the paths.