Python MixedIntegerLinearProgram.set_objective方法代码示例

# 需要导入模块: from sage.numerical.mip import MixedIntegerLinearProgram [as 别名]
# 或者: from sage.numerical.mip.MixedIntegerLinearProgram import set_objective [as 别名]
def _delsarte_LP_building(n, d, d_star, q, isinteger,  solver, maxc = 0):
    """
    LP builder - common for the two functions; not exported.

    EXAMPLES::

        sage: from sage.coding.delsarte_bounds import _delsarte_LP_building
        sage: _,p=_delsarte_LP_building(7, 3, 0, 2, False, "PPL")
        sage: p.show()
        Maximization:
          x_0 + x_1 + x_2 + x_3 + x_4 + x_5 + x_6 + x_7
        Constraints:
          constraint_0: 1 <= x_0 <= 1
          constraint_1: 0 <= x_1 <= 0
          constraint_2: 0 <= x_2 <= 0
          constraint_3: -7 x_0 - 5 x_1 - 3 x_2 - x_3 + x_4 + 3 x_5 + 5 x_6 + 7 x_7 <= 0
          constraint_4: -7 x_0 - 5 x_1 - 3 x_2 - x_3 + x_4 + 3 x_5 + 5 x_6 + 7 x_7 <= 0
          constraint_5: -21 x_0 - 9 x_1 - x_2 + 3 x_3 + 3 x_4 - x_5 - 9 x_6 - 21 x_7 <= 0
          constraint_6: -21 x_0 - 9 x_1 - x_2 + 3 x_3 + 3 x_4 - x_5 - 9 x_6 - 21 x_7 <= 0
          constraint_7: -35 x_0 - 5 x_1 + 5 x_2 + 3 x_3 - 3 x_4 - 5 x_5 + 5 x_6 + 35 x_7 <= 0
          constraint_8: -35 x_0 - 5 x_1 + 5 x_2 + 3 x_3 - 3 x_4 - 5 x_5 + 5 x_6 + 35 x_7 <= 0
          constraint_9: -35 x_0 + 5 x_1 + 5 x_2 - 3 x_3 - 3 x_4 + 5 x_5 + 5 x_6 - 35 x_7 <= 0
          constraint_10: -35 x_0 + 5 x_1 + 5 x_2 - 3 x_3 - 3 x_4 + 5 x_5 + 5 x_6 - 35 x_7 <= 0
          constraint_11: -21 x_0 + 9 x_1 - x_2 - 3 x_3 + 3 x_4 + x_5 - 9 x_6 + 21 x_7 <= 0
          constraint_12: -21 x_0 + 9 x_1 - x_2 - 3 x_3 + 3 x_4 + x_5 - 9 x_6 + 21 x_7 <= 0
          constraint_13: -7 x_0 + 5 x_1 - 3 x_2 + x_3 + x_4 - 3 x_5 + 5 x_6 - 7 x_7 <= 0
          constraint_14: -7 x_0 + 5 x_1 - 3 x_2 + x_3 + x_4 - 3 x_5 + 5 x_6 - 7 x_7 <= 0
          constraint_15: - x_0 + x_1 - x_2 + x_3 - x_4 + x_5 - x_6 + x_7 <= 0
          constraint_16: - x_0 + x_1 - x_2 + x_3 - x_4 + x_5 - x_6 + x_7 <= 0
        Variables:
          x_0 is a continuous variable (min=0, max=+oo)
          x_1 is a continuous variable (min=0, max=+oo)
          x_2 is a continuous variable (min=0, max=+oo)
          x_3 is a continuous variable (min=0, max=+oo)
          x_4 is a continuous variable (min=0, max=+oo)
          x_5 is a continuous variable (min=0, max=+oo)
          x_6 is a continuous variable (min=0, max=+oo)
          x_7 is a continuous variable (min=0, max=+oo)

    """
    from sage.numerical.mip import MixedIntegerLinearProgram

    p = MixedIntegerLinearProgram(maximization=True, solver=solver)
    A = p.new_variable(integer=isinteger, nonnegative=not isinteger) # A>=0 is assumed
    p.set_objective(sum([A[r] for r in xrange(n+1)]))
    p.add_constraint(A[0]==1)
    for i in xrange(1,d):
        p.add_constraint(A[i]==0)
    for j in xrange(1,n+1):
        rhs = sum([Krawtchouk(n,q,j,r,check=False)*A[r] for r in xrange(n+1)])
        p.add_constraint(0*A[0] <= rhs)
        if j >= d_star:
          p.add_constraint(0*A[0] <= rhs)
        else: # rhs is proportional to j-th weight of the dual code
          p.add_constraint(0*A[0] == rhs)

    if maxc > 0:
        p.add_constraint(sum([A[r] for r in xrange(n+1)]), max=maxc)
    return A, p

# 需要导入模块: from sage.numerical.mip import MixedIntegerLinearProgram [as 别名]
# 或者: from sage.numerical.mip.MixedIntegerLinearProgram import set_objective [as 别名]
    def packing(self, solver=None, verbose=0):
        r"""
        Return a maximum packing

        A maximum packing in a hypergraph is collection of disjoint sets/blocks
        of maximal cardinality. This problem is NP-complete in general, and in
        particular on 3-uniform hypergraphs. It is solved here with an Integer
        Linear Program.

        For more information, see the :wikipedia:`Packing_in_a_hypergraph`.

        INPUT:

        - ``solver`` -- (default: ``None``) Specify a Linear Program (LP)
          solver to be used. If set to ``None``, the default one is used. For
          more information on LP solvers and which default solver is used, see
          the method
          :meth:`solve <sage.numerical.mip.MixedIntegerLinearProgram.solve>`
          of the class
          :class:`MixedIntegerLinearProgram <sage.numerical.mip.MixedIntegerLinearProgram>`.

        - ``verbose`` -- integer (default: ``0``). Sets the level of
          verbosity. Set to 0 by default, which means quiet.
          Only useful when ``algorithm == "LP"``.

        EXAMPLE::

            sage; IncidenceStructure([[1,2],[3,"A"],[2,3]]).packing()
            [[1, 2], [3, 'A']]
            sage: len(designs.steiner_triple_system(9).packing())
            3
        """
        from sage.numerical.mip import MixedIntegerLinearProgram

        # List of blocks containing a given point x
        d = [[] for x in self._points]
        for i,B in enumerate(self._blocks):
            for x in B:
                d[x].append(i)

        p = MixedIntegerLinearProgram(solver=solver)
        b = p.new_variable(binary=True)
        for x,L in enumerate(d): # Set of disjoint blocks
            p.add_constraint(p.sum([b[i] for i in L]) <= 1)

        # Maximum number of blocks
        p.set_objective(p.sum([b[i] for i in range(self.num_blocks())]))

        p.solve(log=verbose)

        return [[self._points[x] for x in self._blocks[i]]
                for i,v in p.get_values(b).iteritems() if v]

# 需要导入模块: from sage.numerical.mip import MixedIntegerLinearProgram [as 别名]
# 或者: from sage.numerical.mip.MixedIntegerLinearProgram import set_objective [as 别名]
    def _solve_linear_program(self, target):
        r"""
        Return the coefficients of a linear combination to write ``target`` in
        terms of the generators of this semigroup.

        Return ``None`` if no such combination exists.

        EXAMPLES::

            sage: from sage.rings.valuation.value_group import DiscreteValueSemigroup
            sage: D = DiscreteValueSemigroup([2,3,5])
            sage: D._solve_linear_program(12)
            {0: 1, 1: 0, 2: 2}
            sage: 1*2 + 0*3 + 2*5
            12

        """
        if len(self._generators) == 0:
            if target == 0:
                return {}
            else:
                return None

        if len(self._generators) == 1:
            from sage.rings.all import NN
            exp = target / self._generators[0]
            if exp not in NN:
                return None
            return {0 : exp}

        if len(self._generators) == 2 and self._generators[0] == - self._generators[1]:
            from sage.rings.all import ZZ
            exp = target / self._generators[0]
            if exp not in ZZ:
                return None
            return {0: exp, 1: 0}

        from sage.numerical.mip import MixedIntegerLinearProgram, MIPSolverException
        P = MixedIntegerLinearProgram(maximization=False, solver="ppl")
        x = P.new_variable(integer=True, nonnegative=True)
        constraint = sum([g*x[i] for i,g in enumerate(self._generators)]) == target
        P.add_constraint(constraint)
        P.set_objective(None)
        try:
            P.solve()
        except MIPSolverException:
            return None
        return P.get_values(x)

# 需要导入模块: from sage.numerical.mip import MixedIntegerLinearProgram [as 别名]
# 或者: from sage.numerical.mip.MixedIntegerLinearProgram import set_objective [as 别名]
def knapsack(seq, binary=True, max=1, value_only=False, solver=None, verbose=0):
    r"""
    Solves the knapsack problem

    For more information on the knapsack problem, see the documentation of the
    :mod:`knapsack module <sage.numerical.knapsack>` or the
    :wikipedia:`Knapsack_problem`.

    INPUT:

    - ``seq`` -- Two different possible types:

      - A sequence of tuples ``(weight, value, something1, something2,
        ...)``. Note that only the first two coordinates (``weight`` and
        ``values``) will be taken into account. The rest (if any) will be
        ignored. This can be useful if you need to attach some information to
        the items.

      - A sequence of reals (a value of 1 is assumed).

    - ``binary`` -- When set to ``True``, an item can be taken 0 or 1 time.
      When set to ``False``, an item can be taken any amount of times (while
      staying integer and positive).

    - ``max`` -- Maximum admissible weight.

    - ``value_only`` -- When set to ``True``, only the maximum useful value is
      returned. When set to ``False``, both the maximum useful value and an
      assignment are returned.

    - ``solver`` -- (default: ``None``) Specify a Linear Program (LP) solver to
      be used. If set to ``None``, the default one is used. For more information
      on LP solvers and which default solver is used, see the documentation of
      class :class:`MixedIntegerLinearProgram
      <sage.numerical.mip.MixedIntegerLinearProgram>`.

    - ``verbose`` -- integer (default: ``0``). Sets the level of verbosity. Set
      to 0 by default, which means quiet.

    OUTPUT:

    If ``value_only`` is set to ``True``, only the maximum useful value is
    returned. Else (the default), the function returns a pair ``[value,list]``,
    where ``list`` can be of two types according to the type of ``seq``:

    - The list of tuples `(w_i, u_i, ...)` occurring in the solution.

    - A list of reals where each real is repeated the number of times it is
      taken into the solution.

    EXAMPLES:

    If your knapsack problem is composed of three items ``(weight, value)``
    defined by ``(1,2), (1.5,1), (0.5,3)``, and a bag of maximum weight `2`, you
    can easily solve it this way::

        sage: from sage.numerical.knapsack import knapsack
        sage: knapsack( [(1,2), (1.5,1), (0.5,3)], max=2)
        [5.0, [(1, 2), (0.500000000000000, 3)]]

        sage: knapsack( [(1,2), (1.5,1), (0.5,3)], max=2, value_only=True)
        5.0

    Besides weight and value, you may attach any data to the items::

        sage: from sage.numerical.knapsack import knapsack
        sage: knapsack( [(1, 2, 'spam'), (0.5, 3, 'a', 'lot')])
        [3.0, [(0.500000000000000, 3, 'a', 'lot')]]

    In the case where all the values (usefulness) of the items are equal to one,
    you do not need embarrass yourself with the second values, and you can just
    type for items `(1,1), (1.5,1), (0.5,1)` the command::

        sage: from sage.numerical.knapsack import knapsack
        sage: knapsack([1,1.5,0.5], max=2, value_only=True)
        2.0
    """
    reals = not isinstance(seq[0], tuple)
    if reals:
        seq = [(x,1) for x in seq]

    from sage.numerical.mip import MixedIntegerLinearProgram
    p = MixedIntegerLinearProgram(solver=solver, maximization=True)

    if binary:
        present = p.new_variable(binary = True)
    else:
        present = p.new_variable(integer = True)

    p.set_objective(p.sum([present[i] * seq[i][1] for i in range(len(seq))]))
    p.add_constraint(p.sum([present[i] * seq[i][0] for i in range(len(seq))]), max=max)

    if value_only:
        return p.solve(objective_only=True, log=verbose)

    else:
        objective = p.solve(log=verbose)
        present = p.get_values(present)

        val = []
#.........这里部分代码省略.........

# 需要导入模块: from sage.numerical.mip import MixedIntegerLinearProgram [as 别名]
# 或者: from sage.numerical.mip.MixedIntegerLinearProgram import set_objective [as 别名]

#.........这里部分代码省略.........
            [0 0 0 0]
            [0 0 0 0]
            sage: gale_ryser_theorem([0,0,0],[0,0,0,0], algorithm="ryser")
            [0 0 0 0]
            [0 0 0 0]
            [0 0 0 0]

        Check that :trac:`16638` is fixed::

            sage: tests = [([4, 3, 3, 2, 1, 1, 1, 1, 0], [6, 5, 1, 1, 1, 1, 1]),
            ....:          ([4, 4, 3, 3, 1, 1, 0], [5, 5, 2, 2, 1, 1]),
            ....:          ([4, 4, 3, 2, 1, 1], [5, 5, 1, 1, 1, 1, 1, 0, 0]),
            ....:          ([3, 3, 3, 3, 2, 1, 1, 1, 0], [7, 6, 2, 1, 1, 0]),
            ....:          ([3, 3, 3, 1, 1, 0], [4, 4, 1, 1, 1])]
            sage: for s1, s2 in tests:
            ....:     m = gale_ryser_theorem(s1, s2, algorithm="ryser")
            ....:     ss1 = sorted(map(lambda x : sum(x) , m.rows()), reverse = True)
            ....:     ss2 = sorted(map(lambda x : sum(x) , m.columns()), reverse = True)
            ....:     if ((ss1 != s1) or (ss2 != s2)):
            ....:         print("Error in Ryser algorithm")
            ....:         print(s1, s2)

        REFERENCES:

        ..  [Ryser63] \H. J. Ryser, Combinatorial Mathematics,
            Carus Monographs, MAA, 1963.
        ..  [Gale57] \D. Gale, A theorem on flows in networks, Pacific J. Math.
            7(1957)1073-1082.
        """
        from sage.combinat.partition import Partition
        from sage.matrix.constructor import matrix

        if not(is_gale_ryser(p1,p2)):
            return False

        if algorithm=="ryser": # ryser's algorithm
            from sage.combinat.permutation import Permutation

            # Sorts the sequences if they are not, and remembers the permutation
            # applied
            tmp = sorted(enumerate(p1), reverse=True, key=lambda x:x[1])
            r = [x[1] for x in tmp]
            r_permutation = [x-1 for x in Permutation([x[0]+1 for x in tmp]).inverse()]
            m = len(r)

            tmp = sorted(enumerate(p2), reverse=True, key=lambda x:x[1])
            s = [x[1] for x in tmp]
            s_permutation = [x-1 for x in Permutation([x[0]+1 for x in tmp]).inverse()]
            n = len(s)

            # This is the partition equivalent to the sliding algorithm
            cols = []
            for t in reversed(s):
                c = [0] * m
                i = 0
                while t:
                    k = i + 1
                    while k < m and r[i] == r[k]:
                        k += 1
                    if t >= k - i: # == number rows of the same length
                        for j in range(i, k):
                            r[j] -= 1
                            c[j] = 1
                        t -= k - i
                    else: # Remove the t last rows of that length
                        for j in range(k-t, k):
                            r[j] -= 1
                            c[j] = 1
                        t = 0
                    i = k
                cols.append(c)

            # We added columns to the back instead of the front
            A0 = matrix(list(reversed(cols))).transpose()

            # Applying the permutations to get a matrix satisfying the
            # order given by the input
            A0 = A0.matrix_from_rows_and_columns(r_permutation, s_permutation)
            return A0

        elif algorithm == "gale":
          from sage.numerical.mip import MixedIntegerLinearProgram
          k1, k2=len(p1), len(p2)
          p = MixedIntegerLinearProgram()
          b = p.new_variable(binary = True)
          for (i,c) in enumerate(p1):
              p.add_constraint(p.sum([b[i,j] for j in xrange(k2)]) ==c)
          for (i,c) in enumerate(p2):
              p.add_constraint(p.sum([b[j,i] for j in xrange(k1)]) ==c)
          p.set_objective(None)
          p.solve()
          b = p.get_values(b)
          M = [[0]*k2 for i in xrange(k1)]
          for i in xrange(k1):
              for j in xrange(k2):
                  M[i][j] = int(b[i,j])
          return matrix(M)

        else:
            raise ValueError("The only two algorithms available are \"gale\" and \"ryser\"")

# 需要导入模块: from sage.numerical.mip import MixedIntegerLinearProgram [as 别名]
# 或者: from sage.numerical.mip.MixedIntegerLinearProgram import set_objective [as 别名]

#.........这里部分代码省略.........
        and is characterized by the following property: it is not empty and each
        block either contains `0` or `s` points of this arc. Equivalently, the
        trace of the BIBD on these points is again a BIBD (with block size `s`).

        For more informations, see :wikipedia:`Arc_(projective_geometry)`.

        INPUT:

        - ``s`` - (default to ``2``) the maximum number of points from the arc
          in each block

        - ``solver`` -- (default: ``None``) Specify a Linear Program (LP)
          solver to be used. If set to ``None``, the default one is used. For
          more information on LP solvers and which default solver is used, see
          the method
          :meth:`solve <sage.numerical.mip.MixedIntegerLinearProgram.solve>`
          of the class
          :class:`MixedIntegerLinearProgram <sage.numerical.mip.MixedIntegerLinearProgram>`.

        - ``verbose`` -- integer (default: ``0``). Sets the level of
          verbosity. Set to 0 by default, which means quiet.

        EXAMPLES::

            sage: B = designs.balanced_incomplete_block_design(21, 5)
            sage: a2 = B.arc()
            sage: a2 # random
            [5, 9, 10, 12, 15, 20]
            sage: len(a2)
            6
            sage: a4 = B.arc(4)
            sage: a4 # random
            [0, 1, 2, 5, 6, 8, 9, 10, 11, 12, 13, 14, 15, 16, 18, 20]
            sage: len(a4)
            16

        The `2`-arc and `4`-arc above are maximal. One can check that they
        intersect the blocks in either 0 or `s` points. Or equivalently that the
        traces are again BIBD::

            sage: r = (21-1)/(5-1)
            sage: 1 + r*1
            6
            sage: 1 + r*3
            16

            sage: B.trace(a2).is_t_design(2, return_parameters=True)
            (True, (2, 6, 2, 1))
            sage: B.trace(a4).is_t_design(2, return_parameters=True)
            (True, (2, 16, 4, 1))

        Some other examples which are not maximal::

            sage: B = designs.balanced_incomplete_block_design(25, 4)
            sage: a2 = B.arc(2)
            sage: r = (25-1)/(4-1)
            sage: print len(a2), 1 + r
            8 9
            sage: sa2 = set(a2)
            sage: set(len(sa2.intersection(b)) for b in B.blocks())
            {0, 1, 2}
            sage: B.trace(a2).is_t_design(2)
            False

            sage: a3 = B.arc(3)
            sage: print len(a3), 1 + 2*r
            15 17
            sage: sa3 = set(a3)
            sage: set(len(sa3.intersection(b)) for b in B.blocks()) == set([0,3])
            False
            sage: B.trace(a3).is_t_design(3)
            False

        TESTS:

        Test consistency with relabeling::

            sage: b = designs.balanced_incomplete_block_design(7,3)
            sage: b.relabel(list("abcdefg"))
            sage: set(b.arc()).issubset(b.ground_set())
            True
        """
        s = int(s)

        # trivial cases
        if s <= 0:
            return []
        elif s >= max(self.block_sizes()):
            return self._points[:]

        # linear program
        from sage.numerical.mip import MixedIntegerLinearProgram

        p = MixedIntegerLinearProgram(solver=solver)
        b = p.new_variable(binary=True)
        p.set_objective(p.sum(b[i] for i in range(len(self._points))))
        for i in self._blocks:
            p.add_constraint(p.sum(b[k] for k in i) <= s)
        p.solve(log=verbose)
        return [self._points[i] for (i,j) in p.get_values(b).items() if j == 1]

# 需要导入模块: from sage.numerical.mip import MixedIntegerLinearProgram [as 别名]
# 或者: from sage.numerical.mip.MixedIntegerLinearProgram import set_objective [as 别名]

#.........这里部分代码省略.........

    - ``items`` -- A list of real values (the items' weight)

    - ``maximum``   -- The maximal size of a bin

    - ``k``     -- Number of bins

      - When set to an integer value, the function returns a partition
        of the items into `k` bins if possible, and raises an
        exception otherwise.

      - When set to ``None``, the function returns a partition of the items
        using the least number possible of bins.

    OUTPUT:

    A list of lists, each member corresponding to a box and containing
    the list of the weights inside it. If there is no solution, an
    exception is raised (this can only happen when ``k`` is specified
    or if ``maximum`` is less that the size of one item).

    EXAMPLES:

    Trying to find the minimum amount of boxes for 5 items of weights
    `1/5, 1/4, 2/3, 3/4, 5/7`::

        sage: from sage.numerical.optimize import binpacking
        sage: values = [1/5, 1/3, 2/3, 3/4, 5/7]
        sage: bins = binpacking(values)
        sage: len(bins)
        3

    Checking the bins are of correct size ::

        sage: all([ sum(b)<= 1 for b in bins ])
        True

    Checking every item is in a bin ::

        sage: b1, b2, b3 = bins
        sage: all([ (v in b1 or v in b2 or v in b3) for v in values ])
        True

    One way to use only three boxes (which is best possible) is to put
    `1/5 + 3/4` together in a box, `1/3+2/3` in another, and `5/7`
    by itself in the third one.

    Of course, we can also check that there is no solution using only two boxes ::

        sage: from sage.numerical.optimize import binpacking
        sage: binpacking([0.2,0.3,0.8,0.9], k=2)
        Traceback (most recent call last):
        ...
        ValueError: This problem has no solution !
    """

    if max(items) > maximum:
        raise ValueError("This problem has no solution !")

    if k==None:
        from sage.functions.other import ceil
        k=ceil(sum(items)/maximum)
        while True:
            from sage.numerical.mip import MIPSolverException
            try:
                return binpacking(items,k=k,maximum=maximum)
            except MIPSolverException:
                k = k + 1

    from sage.numerical.mip import MixedIntegerLinearProgram, MIPSolverException
    p=MixedIntegerLinearProgram()

    # Boolean variable indicating whether
    # the i th element belongs to box b
    box=p.new_variable(dim=2)

    # Each bin contains at most max
    for b in range(k):
        p.add_constraint(p.sum([items[i]*box[i][b] for i in range(len(items))]),max=maximum)

    # Each item is assigned exactly one bin
    for i in range(len(items)):
        p.add_constraint(p.sum([box[i][b] for b in range(k)]),min=1,max=1)

    p.set_objective(None)
    p.set_binary(box)

    try:
        p.solve()
    except MIPSolverException:
        raise ValueError("This problem has no solution !")

    box=p.get_values(box)

    boxes=[[] for i in range(k)]

    for b in range(k):
        boxes[b].extend([items[i] for i in range(len(items)) if round(box[i][b])==1])

    return boxes

# 需要导入模块: from sage.numerical.mip import MixedIntegerLinearProgram [as 别名]
# 或者: from sage.numerical.mip.MixedIntegerLinearProgram import set_objective [as 别名]
def knapsack(seq, binary=True, max=1, value_only=False):
    r"""
    Solves the knapsack problem

    Knapsack problems:

    You have already had a knapsack problem, so you should know,
    but in case you do not, a knapsack problem is what happens
    when you have hundred of items to put into a bag which is
    too small for all of them.

    When you formally write it, here is your problem:

    * Your bag can contain a weight of at most `W`.
    * Each item `i` you have has a weight `w_i`.
    * Each item `i` has a usefulness of `u_i`.

    You then want to maximize the usefulness of the items you
    will store into your bag, while keeping sure the weight of
    the bag will not go over `W`.

    As a linear program, this problem can be represented this way
    (if you define `b_i` as the binary variable indicating whether
    the item `i` is to be included in your bag):

    .. MATH::

        \mbox{Maximize: }\sum_i b_i u_i \\
        \mbox{Such that: }
        \sum_i b_i w_i \leq W \\
        \forall i, b_i \mbox{ binary variable} \\

    (For more information,
    cf. http://en.wikipedia.org/wiki/Knapsack_problem.)

    EXAMPLE:

    If your knapsack problem is composed of three
    items (weight, value) defined by (1,2), (1.5,1), (0.5,3),
    and a bag of maximum weight 2, you can easily solve it this way::

        sage: from sage.numerical.knapsack import knapsack
        sage: knapsack( [(1,2), (1.5,1), (0.5,3)], max=2)
        [5.0, [(1, 2), (0.500000000000000, 3)]]

        sage: knapsack( [(1,2), (1.5,1), (0.5,3)], max=2, value_only=True)
        5.0

    In the case where all the values (usefulness) of the items
    are equal to one, you do not need embarrass yourself with
    the second values, and you can just type for items
    `(1,1), (1.5,1), (0.5,1)` the command::

        sage: from sage.numerical.knapsack import knapsack
        sage: knapsack([1,1.5,0.5], max=2, value_only=True)
        2.0

    INPUT:

    - ``seq`` -- Two different possible types:

      - A sequence of pairs (weight, value).
      - A sequence of reals (a value of 1 is assumed).

    - ``binary`` -- When set to True, an item can be taken 0 or 1 time.
      When set to False, an item can be taken any amount of
      times (while staying integer and positive).

    - ``max`` -- Maximum admissible weight.

    - ``value_only`` -- When set to True, only the maximum useful
      value is returned. When set to False, both the maximum useful
      value and an assignment are returned.

    OUTPUT:

    If ``value_only`` is set to True, only the maximum useful value
    is returned. Else (the default), the function returns a pair
    ``[value,list]``, where ``list`` can be of two types according
    to the type of ``seq``:

    - A list of pairs `(w_i, u_i)` for each object `i` occurring
      in the solution.
    - A list of reals where each real is repeated the number
      of times it is taken into the solution.
    """
    reals = not isinstance(seq[0], tuple)
    if reals:
        seq = [(x, 1) for x in seq]

    from sage.numerical.mip import MixedIntegerLinearProgram

    p = MixedIntegerLinearProgram(maximization=True)
    present = p.new_variable()
    p.set_objective(p.sum([present[i] * seq[i][1] for i in range(len(seq))]))
    p.add_constraint(p.sum([present[i] * seq[i][0] for i in range(len(seq))]), max=max)

    if binary:
        p.set_binary(present)
    else:
#.........这里部分代码省略.........

# 需要导入模块: from sage.numerical.mip import MixedIntegerLinearProgram [as 别名]
# 或者: from sage.numerical.mip.MixedIntegerLinearProgram import set_objective [as 别名]
def min_cover(npts, sets, solver='sage'):
    r"""
    EXAMPLES::

        sage: from max_plus.rank import min_cover
        sage: min_cover(5, [[0,1,2],[1,2,3],[2,4]], solver='sage')
        3
        sage: min_cover(5, [[0,1,2],[1,2,3],[2,4]], solver='lp_solve')   # optional -- lp_solve
        3
    """
    # check if the problem is solvable
    covered = [False for i in range(npts)]
    for set in sets:
        for ndx in set:
            covered[ndx] = True
    for c in covered:
        if not c:
            return False

    # Write the Free MPS format integer programming problem
    fh = open("tmp-rank-fmps", "w")
    fh.write("NAME min cover\n")
    fh.write("ROWS\n") # constraints
    for i in range(npts):
        fh.write(" G POINT" + str(i) + "\n")
    fh.write(" N NUMSETS\n")
    fh.write("COLUMNS\n") # variables
    for i in range(len(sets)):
        fh.write("  SET%d NUMSETS 1\n" % i)
        for point in sets[i]:
            fh.write("  SET%d POINT%d 1\n" % (i, point))
    fh.write("RHS\n") # right hand side to constraints
    for i in range(npts):
        fh.write("  COVER POINT%d 1\n" % i)
    fh.write("BOUNDS\n") # bounds on variables
    for i in range(len(sets)):
        fh.write(" BV A SET%d\n" % i)
    fh.write("ENDATA\n")
    fh.close()

    # Run the solver
    if solver == 'GLPK' or solver == 'glpk':
        # GLPK solver
        os.system("glpsol -w tmp-rank-glpsol tmp-rank-fmps > /dev/null")
        fh = open("tmp-rank-glpsol")
        fh.readline()
        line = fh.readline()
        min = int(line.split()[1])
        fh.close()
    elif solver == 'lp_solve':
        # lpsolve solver
        p = Popen(["lp_solve", "-fmps", "tmp-rank-fmps"], stdout=PIPE)
        fh = p.stdout
        fh.readline() # blank
        line = fh.readline()
        start = "Value of objective function: "
        if line[0:len(start)] != start:
            stderr.write("Unexpected output from lp_solve\n")
            min = 0
        else:
            min = int(line[len(start):])

        # read to the end of the output without storing anything
        while line:
            line = fh.readline()
        p.communicate()
        fh.close()

    elif solver == 'sage':
        from sage.numerical.mip import MixedIntegerLinearProgram
        M = MixedIntegerLinearProgram(maximization=False)
        x = M.new_variable(binary=True)

        nsets = len(sets)
        dual_sets = [[] for _ in range(npts)]

        for i,s in enumerate(sets):
            for k in s:
                dual_sets[k].append(i)

        for k in range(npts):
            M.add_constraint(M.sum(x[i] for i in dual_sets[k]) >= 1)

        M.set_objective(M.sum(x[i] for i in range(nsets)))

        min = int(M.solve())

    return min

# 需要导入模块: from sage.numerical.mip import MixedIntegerLinearProgram [as 别名]
# 或者: from sage.numerical.mip.MixedIntegerLinearProgram import set_objective [as 别名]

#.........这里部分代码省略.........
        if author['hlpful_fav_unfav']:
            res += th['logPrH']
        else:
            res += th['log1-PrH']
        if author['isRealName']:
            res += th['logPrR']
        else:
            res += th['log1-PrR']
        res += -((author['revLen'] - th['muL']) ** 2) / (2 * th['sigma2L'] + EPS) - (log_2pi + log(th['sigma2L'])) / 2.0
        return res

    def log_likelihood(self):
        ll = sum(self.log_phi(a, self.partition[a]) for a in self.author_graph.nodes())
        log_TAU, log_1_TAU = log(self.TAU), log(1 - self.TAU)
        for a, b in self.author_graph.edges():
            if self.partition[a] == self.partition[b]:
                ll += log_TAU * self.author_graph[a][b]['weight']
            else:
                ll += log_1_TAU * self.author_graph[a][b]['weight']
        return ll

    def e_step(self):
        slog('E-Step')
        if not self._lp_init:
            self._init_LP()

        slog('Obj func - indiv potentials')
        # individual potentials
        for a in self.author_graph:
            for p in self.parts:
                self.objF_indiv_dict[self.alpha_dict[a][p]] = -self.log_phi(a, p)

        objF_indiv = self.lp(self.objF_indiv_dict)
        self.lp.set_objective(self.lp.sum([objF_indiv, self.objF_pair]))

        # solve the LP
        slog('Solving the LP')
        self.lp.solve(log=3)
        slog('Solving the LP Done')

        # hard partitions for nodes (authors)
        self.partition = {}
        for a in self.author_graph:
            membship = self.lp.get_values(self.alpha[a])
            self.partition[a] = max(membship, key=membship.get)
        slog('E-Step Done')

    def m_step(self):
        slog('M-Step')
        stat = {p: [0.0] * len(self.parts) for p in ['freq', 'hlpful', 'realNm', 'muL', 'M2']}
        for a in self.author_graph:
            p = self.partition[a]
            author = self.author_graph.node[a]
            stat['freq'][p] += 1
            if author['hlpful_fav_unfav']: stat['hlpful'][p] += 1
            if author['isRealName']: stat['realNm'][p] += 1
            delta = author['revLen'] - stat['muL'][p]
            stat['muL'][p] += delta / stat['freq'][p]
            stat['M2'][p] += delta * (author['revLen'] - stat['muL'][p])

        self.theta = [{p: 0.0 for p in ['logPr', 'logPrH', 'log1-PrH', 'logPrR', 'log1-PrR', 'sigma2L', 'muL']}
                      for p in self.parts]
        sum_freq = sum(stat['freq'][p] for p in self.parts)

        for p in self.parts:
            self.theta[p]['logPr'] = log(stat['freq'][p] / (sum_freq + EPS) + EPS)

# 需要导入模块: from sage.numerical.mip import MixedIntegerLinearProgram [as 别名]
# 或者: from sage.numerical.mip.MixedIntegerLinearProgram import set_objective [as 别名]

#.........这里部分代码省略.........
	+        answering very quickly in some cases
	+        
	+        EXAMPLES:
	+
	+        The Steiner Tree of the first 5 vertices in a random graph is,
	+        of course, always a tree ::
	+
	+            sage: g = graphs.RandomGNP(30,.5)
	+            sage: st = g.steiner_tree(g.vertices()[:5])              # optional - requires GLPK, CBC or CPLEX
	+            sage: st.is_tree()                                       # optional - requires GLPK, CBC or CPLEX
	+            True
	+
	+        And all the 5 vertices are contained in this tree ::
	+

	+            sage: all([v in st for v in g.vertices()[:5] ])          # optional - requires GLPK, CBC or CPLEX
	+            True
	+
	+        An exception is raised when the problem is impossible, i.e.
	+        if the given vertices are not all included in the same
	+        connected component ::
	+
	+            sage: g = 2 * graphs.PetersenGraph()
	+            sage: st = g.steiner_tree([5,15])
	+            Traceback (most recent call last):
	+            ...
	+            ValueError: The given vertices do not all belong to the same connected component. This problem has no solution !
	+
	"""
	if G.is_directed():
		g = nx.Graph(G)
	else:
		g = G
	vertices=G.nodes()
	if g.has_multiple_edges():
		raise ValueError("The graph is expected not to have multiple edges.")
	# Can the problem be solved ? Are all the vertices in the same
	# connected component ?
	cc = g.connected_component_containing_vertex(vertices[0])
	if not all([v in cc for v in vertices]):
		raise ValueError("The given vertices do not all belong to the same connected component. This problem has no solution !")
	# Can it be solved using the min spanning tree algorithm ?
	if not weighted:
		gg = g.subgraph(vertices)
	if gg.is_connected():
		st = g.subgraph(edges = gg.min_spanning_tree())
		st.delete_vertices([v for v in g if st.degree(v) == 0])
		return st
	# Then, LP formulation
	p = MixedIntegerLinearProgram(maximization = False)
	
	# Reorder an edge
	R = lambda (x,y) : (x,y) if x<y else (y,x)
	
	
	# edges used in the Steiner Tree
	edges = p.new_variable()
	
	# relaxed edges to test for acyclicity
	r_edges = p.new_variable()

	# Whether a vertex is in the Steiner Tree
	vertex = p.new_variable()
	for v in g:
		for e in g.edges_incident(v, labels=False):
			p.add_constraint(vertex[v] - edges[R(e)], min = 0)
	
	# We must have the given vertices in our tree
	for v in vertices:
		p.add_constraint(sum([edges[R(e)] for e in g.edges_incident(v,labels=False)]), min=1)

	# The number of edges is equal to the number of vertices in our tree minus 1
	p.add_constraint(sum([vertex[v] for v in g]) - sum([edges[R(e)] for e in g.edges(labels=None)]), max = 1, min = 1)
	
	# There are no cycles in our graph
	
	for u,v in g.edges(labels = False):
		p.add_constraint( r_edges[(u,v)]+ r_edges[(v,u)] - edges[R((u,v))] , min = 0 )
		
	eps = 1/(5*Integer(g.order()))
	for v in g:
		p.add_constraint(sum([r_edges[(u,v)] for u in g.neighbors(v)]), max = 1-eps)
		
		
	# Objective
	if weighted:
		w = lambda (x,y) : g.edge_label(x,y) if g.edge_label(x,y) is not None else 1
	else:
		w = lambda (x,y) : 1
		
	p.set_objective(sum([w(e)*edges[R(e)] for e in g.edges(labels = False)]))
	
	p.set_binary(edges)     
	p.solve()
	
	edges = p.get_values(edges)
	
	st =  g.subgraph(edges=[e for e in g.edges(labels = False) if edges[R(e)] == 1])
	st.delete_vertices([v for v in g if st.degree(v) == 0])
	return st

# 需要导入模块: from sage.numerical.mip import MixedIntegerLinearProgram [as 别名]
# 或者: from sage.numerical.mip.MixedIntegerLinearProgram import set_objective [as 别名]

#.........这里部分代码省略.........
        ``gale_ryser_theorem`` is then called on these sequences, and the output
        checked for correction.::

            sage: def test_algorithm(algorithm, low = 10, high = 50):
            ...      n,m = randint(low,high), randint(low,high)
            ...      g = graphs.RandomBipartite(n, m, .3)
            ...      s1 = sorted(g.degree([(0,i) for i in range(n)]), reverse = True)
            ...      s2 = sorted(g.degree([(1,i) for i in range(m)]), reverse = True)
            ...      m = gale_ryser_theorem(s1, s2, algorithm = algorithm)
            ...      ss1 = sorted(map(lambda x : sum(x) , m.rows()), reverse = True)
            ...      ss2 = sorted(map(lambda x : sum(x) , m.columns()), reverse = True)
            ...      if ((ss1 == s1) and (ss2 == s2)): 
            ...          return True
            ...      return False

            sage: for algorithm in ["gale", "ryser"]:                            # long time
            ...      for i in range(50):                                         # long time
            ...         if not test_algorithm(algorithm, 3, 10):                 # long time
            ...             print "Something wrong with algorithm ", algorithm   # long time
            ...             break                                                # long time

        Null matrix::

            sage: gale_ryser_theorem([0,0,0],[0,0,0,0], algorithm="gale")
            [0 0 0 0]
            [0 0 0 0]
            [0 0 0 0]
            sage: gale_ryser_theorem([0,0,0],[0,0,0,0], algorithm="ryser")
            [0 0 0 0]
            [0 0 0 0]
            [0 0 0 0]

        REFERENCES:

        ..  [Ryser63] H. J. Ryser, Combinatorial Mathematics, 
                Carus Monographs, MAA, 1963. 
        ..  [Gale57] D. Gale, A theorem on flows in networks, Pacific J. Math.
                7(1957)1073-1082.
        """
        from sage.combinat.partition import Partition
        from sage.matrix.constructor import matrix

        if not(is_gale_ryser(p1,p2)):
            return False

        if algorithm=="ryser": # ryser's algorithm
            from sage.combinat.permutation import Permutation

            # Sorts the sequences if they are not, and remembers the permutation
            # applied 
            tmp = sorted(enumerate(p1), reverse=True, key=lambda x:x[1])
            r = [x[1] for x in tmp if x[1]>0]
            r_permutation = [x-1 for x in Permutation([x[0]+1 for x in tmp]).inverse()]
            m = len(r)

            tmp = sorted(enumerate(p2), reverse=True, key=lambda x:x[1])
            s = [x[1] for x in tmp if x[1]>0]
            s_permutation = [x-1 for x in Permutation([x[0]+1 for x in tmp]).inverse()]
            n = len(s)

            A0 = matrix([[1]*r[j]+[0]*(n-r[j]) for j in range(m)])

            for k in range(1,n+1):
                goodcols = [i for i in range(n) if s[i]==sum(A0.column(i))]
                if sum(A0.column(n-k))<>s[n-k]:
                    A0 = _slider01(A0,s[n-k],n-k, p1, p2, goodcols)

            # If we need to add empty rows/columns
            if len(p1)!=m:
                A0 = A0.stack(matrix([[0]*n]*(len(p1)-m)))

            if len(p2)!=n:
                A0 = A0.transpose().stack(matrix([[0]*len(p1)]*(len(p2)-n))).transpose()
                
            # Applying the permutations to get a matrix satisfying the
            # order given by the input
            A0 = A0.matrix_from_rows_and_columns(r_permutation, s_permutation)
            return A0    

        elif algorithm == "gale":
          from sage.numerical.mip import MixedIntegerLinearProgram, Sum
          k1, k2=len(p1), len(p2)
          p = MixedIntegerLinearProgram()
          b = p.new_variable(dim=2)
          for (i,c) in enumerate(p1):
              p.add_constraint(Sum([b[i][j] for j in xrange(k2)]),min=c,max=c)
          for (i,c) in enumerate(p2):
              p.add_constraint(Sum([b[j][i] for j in xrange(k1)]),min=c,max=c)
          p.set_objective(None)
          p.set_binary(b)
          p.solve()
          b = p.get_values(b)
          M = [[0]*k2 for i in xrange(k1)]
          for i in xrange(k1):
              for j in xrange(k2):
                  M[i][j] = int(b[i][j])
          return matrix(M)

        else:
            raise ValueError("The only two algorithms available are \"gale\" and \"ryser\"")