bertiniteam / b2

Bertini 2.0: The redevelopment of Bertini in C++.
92 stars 34 forks source link

Precision error with 2 variables and 2 functions and more than 5 runs through the loop #155

Closed mclower11 closed 6 years ago

mclower11 commented 6 years ago

Below is the code for running through bertini with 2 functions, 2 variables. They are random polynomials, some having singularities and some not. If the loop length is greater than 5, it won't run and gives an error message about the precision.

import numpy as np
import random
import pybertini as pyb
import time

average1=[]
average2=[]
x=pyb.function_tree.symbol.Variable("x")
y=pyb.function_tree.symbol.Variable("y")
function1=[]
function2=[]
for pp in range(1,4):
    a=random.randint(-10,10)
    b=random.randint(0,2)
    c=random.randint(-10,10)
    d=random.randint(-10,10)
    e=random.randint(-10,10)
    func1=((x+a)**b)*(y+c)
    func2=(x+d)*(y+e)
    function1.append(func1)
    function2.append(func2)

    start1=time.time()
    sys = pyb.System()
    sys.add_function(func1)
    sys.add_function(func2)
    varg=pyb.VariableGroup()
    varg.append(x)
    varg.append(y)
    sys.add_variable_group(varg)

    #then make a start system
    t=pyb.Variable('t')
    td=pyb.system.start_system.TotalDegree(sys)
    gamma=pyb.function_tree.symbol.Rational.rand()
    hom=(1-t)*sys + t*gamma*td
    hom.add_path_variable(t)
    newton=pyb.tracking.config.NewtonConfig()
    newton.max_num_newton_iterations=2
    newton.min_num_newton_iterations=2
#   then we have to track the homotopy
    tr=pyb.tracking.AMPTracker(hom)
    tr.set_newton(newton)
    start_time=pyb.multiprec.Complex("1")
    eg_boundary=pyb.multiprec.Complex("0.1")

    midpath_points = [None]*td.num_start_points()
    for ii in range(td.num_start_points()):
            midpath_points[ii] = pyb.multiprec.Vector()
            code = tr.track_path(result=midpath_points[ii], start_time=start_time, end_time=eg_boundary, start_point=td.start_point_mp(ii))

    #then we make the endgame
    eg=pyb.endgame.AMPCauchyEG(tr)

    #we then need to actually run the endgame
    assert(eg.cycle_number()==0)
    assert(eg.final_approximation()==pyb.minieigen.VectorXmp())

    #Need to make an empty list where the solutions will be added
    final_points = []

    target_time = pyb.multiprec.Complex(0)
    codes = []
    for ii in range(td.num_start_points()):
            eg_boundary.precision(midpath_points[ii][0].precision())
            target_time.precision(midpath_points[ii][0].precision())
            codes.append(eg.run(start_time=eg_boundary, target_time=target_time, start_point=midpath_points[ii]))

    end1=time.time()
    timee1=end1-start1
    average1.append(timee1)

#this makes the observer for looking at our tracker along the way
#g=pyb.tracking.observers.amp.GoryDetailLogger()

    start2=time.time()
    sys = pyb.System()
    sys.add_function(func1)
    sys.add_function(func2)

    varg=pyb.VariableGroup()
    varg.append(x)
    varg.append(y)
    sys.add_variable_group(varg)

    #then make a start system
    t=pyb.Variable('t')
    td=pyb.system.start_system.TotalDegree(sys)
    gamma=pyb.function_tree.symbol.Rational.rand()
    hom=(1-t)*sys + t*gamma*td
    hom.add_path_variable(t)
    newton=pyb.tracking.config.NewtonConfig()
    newton.max_num_newton_iterations=3
    newton.min_num_newton_iterations=3
    #then we have to track the homotopy
    tr=pyb.tracking.AMPTracker(hom)
    tr.set_newton(newton)
    start_time=pyb.multiprec.Complex("1")
    eg_boundary=pyb.multiprec.Complex("0.1")

    midpath_points = [None]*td.num_start_points()
    for ii in range(td.num_start_points()):
            midpath_points[ii] = pyb.multiprec.Vector()
            code = tr.track_path(result=midpath_points[ii], start_time=start_time, end_time=eg_boundary, start_point=td.start_point_mp(ii))

    #then we make the endgame
    eg=pyb.endgame.AMPCauchyEG(tr)

    #we then need to actually run the endgame
    assert(eg.cycle_number()==0)
    assert(eg.final_approximation()==pyb.minieigen.VectorXmp())

    #Need to make an empty list where the solutions will be added
    final_points = []

    target_time = pyb.multiprec.Complex(0)
    codes = []
    for ii in range(td.num_start_points()):
            eg_boundary.precision(midpath_points[ii][0].precision())
            target_time.precision(midpath_points[ii][0].precision())
            codes.append(eg.run(start_time=eg_boundary, target_time=target_time, start_point=midpath_points[ii]))

    end2=time.time()
    timee2=end2-start2
    average2.append(timee2)

averagetime1=sum(average1)/float(len(average1))
averagetime2=sum(average2)/float(len(average2))

print('this is the average time for 2 steps')
print(averagetime1)
print('this is the average time for 3 steps')
print(averagetime2)

This is the error message given when I try to do a loop of length greater than 5.

python3: /usr/local/include/bertini2/endgames/cauchy.hpp:389: bertini::SuccessCode bertini::endgame::CauchyEndgame<PrecT>::CircleTrack(const CT&) [with CT = bertini::complex; PrecT = bertini::endgame::AMPEndgame]: Assertion `Precision(current_time)==Precision(current_sample) && "current time and sample for circle track must be of same precision"' failed.
Aborted (core dumped)
jbcolli2 commented 6 years ago

Is this the same problem we saw with precision for Issue #154, where the precision is drifting, but the components are not keeping up?

ofloveandhate commented 6 years ago

@mclower11 when did you build the b2 you are using? this might just be due to an old commit for minieigen.

from the root of the b2 repo:

 cd python/minieigen/
git show --oneline -s

should currently show

5e4d9fb

if not, then you need to update minieigen. move to the minieigen directory and run

git pull origin master

if you do have the correct commit, then we have something to do on our hands!

mclower11 commented 6 years ago

It looks like I have the same thing.

5e4d9fb (HEAD, origin/master, origin/HEAD, master) changed get_item to return a reference to the original value, instead of a copy
ofloveandhate commented 6 years ago

i just ran your code once.

this is the average time for 2 steps
8.114491939544678
this is the average time for 3 steps
0.6285444100697836

what OS are you on? what compiler did you use?

ofloveandhate commented 6 years ago

got the error on my fourth try.

ofloveandhate commented 6 years ago

this is a problem for me to solve! i'll try to do so asap.

ofloveandhate commented 6 years ago

this produces the error almost instantaneously:

# Below is the code for running through bertini with 2 functions, 2 variables. They are random polynomials, some having singularities and some not. If the loop length is greater than 5, it won't run and gives an error message about the precision.

import numpy as np
import random
import pybertini as pyb
import time

n = 1;

random.seed(128232541238)

average1=[]
average2=[]
x=pyb.function_tree.symbol.Variable("x")
y=pyb.function_tree.symbol.Variable("y")
function1=[]
function2=[]
for pp in range(n):
    a=random.randint(-10,10)
    b=random.randint(0,2)
    c=random.randint(-10,10)
    d=random.randint(-10,10)
    e=random.randint(-10,10)
    func1=((x+a)**b)*(y+c)
    func2=(x+d)*(y+e)
    function1.append(func1)
    function2.append(func2)

    start1=time.time()
    sys = pyb.System()
    sys.add_function(func1)
    sys.add_function(func2)
    varg=pyb.VariableGroup()
    varg.append(x)
    varg.append(y)
    sys.add_variable_group(varg)

    #then make a start system
    t=pyb.Variable('t')
    td=pyb.system.start_system.TotalDegree(sys)
    gamma=pyb.function_tree.symbol.Rational.rand()
    hom=(1-t)*sys + t*gamma*td
    hom.add_path_variable(t)
    newton=pyb.tracking.config.NewtonConfig()
    newton.max_num_newton_iterations=2
    newton.min_num_newton_iterations=2
#   then we have to track the homotopy
    tr=pyb.tracking.AMPTracker(hom)
    tr.set_newton(newton)
    start_time=pyb.multiprec.Complex("1")
    eg_boundary=pyb.multiprec.Complex("0.1")

    midpath_points = [None]*td.num_start_points()
    for ii in range(td.num_start_points()):
            midpath_points[ii] = pyb.multiprec.Vector()
            code = tr.track_path(result=midpath_points[ii], start_time=start_time, end_time=eg_boundary, start_point=td.start_point_mp(ii))

    #then we make the endgame
    eg=pyb.endgame.AMPCauchyEG(tr)

    #we then need to actually run the endgame
    assert(eg.cycle_number()==0)
    assert(eg.final_approximation()==pyb.minieigen.VectorXmp())

    #Need to make an empty list where the solutions will be added
    final_points = []

    target_time = pyb.multiprec.Complex(0)
    codes = []
    for ii in range(td.num_start_points()):
            eg_boundary.precision(midpath_points[ii][0].precision())
            target_time.precision(midpath_points[ii][0].precision())
            codes.append(eg.run(start_time=eg_boundary, target_time=target_time, start_point=midpath_points[ii]))

    end1=time.time()
    timee1=end1-start1
    average1.append(timee1)

#this makes the observer for looking at our tracker along the way
#g=pyb.tracking.observers.amp.GoryDetailLogger()

    start2=time.time()
    sys = pyb.System()
    sys.add_function(func1)
    sys.add_function(func2)

    varg=pyb.VariableGroup()
    varg.append(x)
    varg.append(y)
    sys.add_variable_group(varg)

    #then make a start system
    t=pyb.Variable('t')
    td=pyb.system.start_system.TotalDegree(sys)
    gamma=pyb.function_tree.symbol.Rational.rand()
    hom=(1-t)*sys + t*gamma*td
    hom.add_path_variable(t)
    newton=pyb.tracking.config.NewtonConfig()
    newton.max_num_newton_iterations=3
    newton.min_num_newton_iterations=3
    #then we have to track the homotopy
    tr=pyb.tracking.AMPTracker(hom)
    tr.set_newton(newton)
    start_time=pyb.multiprec.Complex("1")
    eg_boundary=pyb.multiprec.Complex("0.1")

    midpath_points = [None]*td.num_start_points()
    for ii in range(td.num_start_points()):
            midpath_points[ii] = pyb.multiprec.Vector()
            code = tr.track_path(result=midpath_points[ii], start_time=start_time, end_time=eg_boundary, start_point=td.start_point_mp(ii))

    #then we make the endgame
    eg=pyb.endgame.AMPCauchyEG(tr)

    #we then need to actually run the endgame
    assert(eg.cycle_number()==0)
    assert(eg.final_approximation()==pyb.minieigen.VectorXmp())

    #Need to make an empty list where the solutions will be added
    final_points = []

    target_time = pyb.multiprec.Complex(0)
    codes = []
    for ii in range(td.num_start_points()):
            eg_boundary.precision(midpath_points[ii][0].precision())
            target_time.precision(midpath_points[ii][0].precision())
            codes.append(eg.run(start_time=eg_boundary, target_time=target_time, start_point=midpath_points[ii]))

    end2=time.time()
    timee2=end2-start2
    average2.append(timee2)

averagetime1=sum(average1)/float(len(average1))
averagetime2=sum(average2)/float(len(average2))

print('this is the average time for 2 steps')
print(averagetime1)
print('this is the average time for 3 steps')
print(averagetime2)
ofloveandhate commented 6 years ago

i looked into it a bit today but ran out of time. will take a deeper look tomorrow. sorry you are having this trouble. thanks for raising the issue.

ofloveandhate commented 6 years ago

i got it figured out. the tracker was giving non-success while advancing time, so the endgame should have terminated before getting to the assert. yay asserts! use them! i have a fix incoming. testing using your example now.

the reason the tracker was giving non-success is that the path is going to infinity. your timings will change, probably for the better, if you work projectively since going to infinity shouldn't happen.

ofloveandhate commented 6 years ago

@mclower11 @jbcolli2 take a pull to get the fix.

then, make install to the core, and then recompile and install the c++ part of the python bindings, so just make -j whatever install. no need to re-install through pip since the pure python part of pybertini is unchanged.

please comment and close if this solves the problem for you.

mclower11 commented 6 years ago

It is working now! Thank you

ofloveandhate commented 6 years ago

thankyou!