Revision history [back]

Somehow looks like the numerical precision is not allowing the process to go through. The condition to stop the infinite loop is that the sum of the labels in the matching adds up to 1, in your case:
0.9999999999999999 + 3.3306690738754696e-16+ 1.4802973661668753e-16

However this evaluates to a number that is slightly above 1, this matching is produced over and over ever after. Below is some code to allow eps tolerance to the comparison to 1 so that the loop can be broken. This code is a tweaked copy of the code for fractional_chromatic_index that ships with SageMath. In the meantime this will be reported as a bug in the sage-devel mailing list, thanks for reporting!

def fract_chromatic_index(h,solver = None, verbose_constraints = 0, verbose = 0, eps=0.00000001):
    from sage.numerical.mip import MixedIntegerLinearProgram

    g = copy(h)
    p = MixedIntegerLinearProgram(solver=solver, constraint_generation = True)

    # One variable per edge
    r = p.new_variable(nonnegative=True)
    R = lambda x,y : r[x,y] if x<y else r[y,x]

    # We want to maximize the sum of weights on the edges
    p.set_objective( p.sum( R(u,v) for u,v in g.edges(labels = False)))

    # Each edge being by itself a matching, its weight can not be more than
    # 1

    for u,v in g.edges(labels = False):
        p.add_constraint( R(u,v), max = 1)

    obj = p.solve(log = verbose)
    while True:
        # Updating the value on the edges of g
        for u,v in g.edges(labels = False):
            g.set_edge_label(u,v,p.get_values(R(u,v)))

        # Computing a matching of maximum weight...
        matching = g.matching(use_edge_labels=True)

        # If the maximum matching has weight at most 1, we are done !
        if sum((x[2] for x in matching)) <= 1.0+eps:
            break

        #print "sum((x[2] for x in matching)): {0}; type: {1}".format(sum((x[2] for x in matching)),type(sum((x[2] for x in matching))))
        #print "sum((x[2] for x in matching)) <= 1.0 ? {0}".format(sum((x[2] for x in matching))<=1.0)
        #print matching
        #sleep(1)
        # Otherwise, we add a new constraint
        if verbose_constraints:
            print("Adding a constraint on matching : {}".format(matching))

        p.add_constraint( p.sum( R(u,v) for u,v,_ in matching), max = 1)
        # And solve again
        obj = p.solve(log = verbose)

    # Accomplished !
    return obj

This allows for

h=graphs.PetersenGraph()
fractional_chromatic_index(h,solver="Glpk")

to output 2.9999999999999996.

def fract_chromatic_index(h,solver = None, verbose_constraints = 0, verbose = 0, eps=0.00000001):
    from sage.numerical.mip import MixedIntegerLinearProgram

    g = copy(h)
    p = MixedIntegerLinearProgram(solver=solver, constraint_generation = True)

    # One variable per edge
    r = p.new_variable(nonnegative=True)
    R = lambda x,y : r[x,y] if x<y else r[y,x]

    # We want to maximize the sum of weights on the edges
    p.set_objective( p.sum( R(u,v) for u,v in g.edges(labels = False)))

    # Each edge being by itself a matching, its weight can not be more than
    # 1

    for u,v in g.edges(labels = False):
        p.add_constraint( R(u,v), max = 1)

    obj = p.solve(log = verbose)
    while True:
        # Updating the value on the edges of g
        for u,v in g.edges(labels = False):
            g.set_edge_label(u,v,p.get_values(R(u,v)))

        # Computing a matching of maximum weight...
        matching = g.matching(use_edge_labels=True)

        # If the maximum matching has weight at most 1, we are done !
        if sum((x[2] for x in matching)) <= 1.0+eps:
            break

        #print "sum((x[2] for x in matching)): {0}; type: {1}".format(sum((x[2] for x in matching)),type(sum((x[2] for x in matching))))
        #print "sum((x[2] for x in matching)) <= 1.0 ? {0}".format(sum((x[2] for x in matching))<=1.0)
        #print matching
        #sleep(1)
        # Otherwise, we add a new constraint
        if verbose_constraints:
            print("Adding a constraint on matching : {}".format(matching))

        p.add_constraint( p.sum( R(u,v) for u,v,_ in matching), max = 1)
        # And solve again
        obj = p.solve(log = verbose)

    # Accomplished !
    return obj

This allows for

h=graphs.PetersenGraph()
fractional_chromatic_index(h,solver="Glpk")

to output 2.9999999999999996.