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)
def _init_LP(self):
if self._lp_init:
return
slog('Init LP')
self.lp = MixedIntegerLinearProgram(solver='GLPK', maximization=False)
#self.lp.solver_parameter(backend.glp_simplex_or_intopt, backend.glp_simplex_only) # LP relaxation
self.lp.solver_parameter("iteration_limit", LP_ITERATION_LIMIT)
# self.lp.solver_parameter("timelimit", LP_TIME_LIMIT)
# add constraints once here
# constraints
self.alpha = self.lp.new_variable(dim=2)
beta2 = self.lp.new_variable(dim=2)
beta3 = self.lp.new_variable(dim=3)
# alphas are indicator vars
for a in self.author_graph:
self.lp.add_constraint(sum(self.alpha[a][p] for p in self.parts) == 1)
# beta2 is the sum of beta3s
slog('Init LP - pair constraints')
for a, b in self.author_graph.edges():
if self.author_graph[a][b]['denom'] <= 2:
continue
self.lp.add_constraint(0.5 * sum(beta3[a][b][p] for p in self.parts) - beta2[a][b], min=0, max=0)
for p in self.parts:
self.lp.add_constraint(self.alpha[a][p] - self.alpha[b][p] - beta3[a][b][p], max=0)
self.lp.add_constraint(self.alpha[b][p] - self.alpha[a][p] - beta3[a][b][p], max=0)
# store indiv potential linear function as a dict to improve performance
self.objF_indiv_dict = {}
self.alpha_dict = {}
for a in self.author_graph:
self.alpha_dict[a] = {}
for p in self.parts:
var_id = self.alpha_dict[a][p] = self.alpha[a][p].dict().keys()[0]
self.objF_indiv_dict[var_id] = 0 # init variables coeffs to zero
# pairwise potentials
slog('Obj func - pair potentials')
objF_pair_dict = {}
s = log(1 - self.TAU) - log(self.TAU)
for a, b in self.author_graph.edges():
if self.author_graph[a][b]['denom'] <= 2:
continue
var_id = beta2[a][b].dict().keys()[0]
objF_pair_dict[var_id] = -self.author_graph[a][b]['weight'] * s
self.objF_pair = self.lp(objF_pair_dict)
self._lp_init = True
slog('Init LP Done')
def __init__(self, solver=None):
r"""
Initializes the instance
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>`.
EXAMPLES::
sage: S=SAT(solver="LP"); S
an ILP-based SAT Solver
"""
SatSolver.__init__(self)
self._LP = MixedIntegerLinearProgram()
self._vars = self._LP.new_variable(binary=True)
开发者ID:mcognetta，项目名称:sage，代码行数:22，代码来源:sat_lp.py
示例4: __init__
def __init__(self, solver=None, verbose=0):
r"""
Initializes the instance
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:`~sage.numerical.mip.MixedIntegerLinearProgram.solve` of
the class :class:`~sage.numerical.mip.MixedIntegerLinearProgram`.
- ``verbose`` -- integer (default: ``0``). Sets the level of verbosity
of the LP solver. Set to 0 by default, which means quiet.
EXAMPLES::
sage: S=SAT(solver="LP"); S
an ILP-based SAT Solver
"""
SatSolver.__init__(self)
self._LP = MixedIntegerLinearProgram(solver=solver)
self._LP_verbose = verbose
self._vars = self._LP.new_variable(binary=True)
def arc(self, s=2, solver=None, verbose=0):
r"""
Return the ``s``-arc with maximum cardinality.
A `s`-arc is a subset of points in a BIBD that intersects each block on
at most `s` points. It is one possible generalization of independent set
for graphs.
A simple counting shows that the cardinality of a `s`-arc is at most
`(s-1) * r + 1` where `r` is the number of blocks incident to any point.
A `s`-arc in a BIBD with cardinality `(s-1) * r + 1` is called maximal
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[:]
#.........这里部分代码省略.........
def binpacking(items,maximum=1,k=None):
r"""
Solves the bin packing problem.
The Bin Packing problem is the following :
Given a list of items of weights `p_i` and a real value `K`, what is
the least number of bins such that all the items can be put in the
bins, while keeping sure that each bin contains a weight of at most `K` ?
For more informations : http://en.wikipedia.org/wiki/Bin_packing_problem
Two version of this problem are solved by this algorithm :
* Is it possible to put the given items in `L` bins ?
* What is the assignment of items using the
least number of bins with the given list of items ?
INPUT:
- ``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)):
#.........这里部分代码省略.........
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:
#.........这里部分代码省略.........
class SatLP(SatSolver):
def __init__(self, solver=None):
r"""
Initializes the instance
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>`.
EXAMPLES::
sage: S=SAT(solver="LP"); S
an ILP-based SAT Solver
"""
SatSolver.__init__(self)
self._LP = MixedIntegerLinearProgram()
self._vars = self._LP.new_variable(binary=True)
def var(self):
"""
Return a *new* variable.
EXAMPLES::
sage: S=SAT(solver="LP"); S
an ILP-based SAT Solver
sage: S.var()
1
"""
nvars = n = self._LP.number_of_variables()
while nvars==self._LP.number_of_variables():
n += 1
self._vars[n] # creates the variable if needed
return n
def nvars(self):
"""
Return the number of variables.
EXAMPLES::
sage: S=SAT(solver="LP"); S
an ILP-based SAT Solver
sage: S.var()
1
sage: S.var()
2
sage: S.nvars()
2
"""
return self._LP.number_of_variables()
def add_clause(self, lits):
"""
Add a new clause to set of clauses.
INPUT:
- ``lits`` - a tuple of integers != 0
.. note::
If any element ``e`` in ``lits`` has ``abs(e)`` greater
than the number of variables generated so far, then new
variables are created automatically.
EXAMPLES::
sage: S=SAT(solver="LP"); S
an ILP-based SAT Solver
sage: for u,v in graphs.CycleGraph(6).edges(labels=False):
....: u,v = u+1,v+1
....: S.add_clause((u,v))
....: S.add_clause((-u,-v))
"""
if 0 in lits:
raise ValueError("0 should not appear in the clause: {}".format(lits))
p = self._LP
p.add_constraint(p.sum(self._vars[x] if x>0 else 1-self._vars[-x] for x in lits)
>=1)
def __call__(self):
"""
Solve this instance.
OUTPUT:
- If this instance is SAT: A tuple of length ``nvars()+1``
where the ``i``-th entry holds an assignment for the
``i``-th variables (the ``0``-th entry is always ``None``).
- If this instance is UNSAT: ``False``
EXAMPLES::
#.........这里部分代码省略.........
def solve_magic_hexagon(solver=None):
r"""
Solves the magic hexagon problem
We use the following convention for the positions::
1 2 3
4 5 6 7
8 9 10 11 12
13 14 15 16
17 18 19
INPUT:
- ``solver`` -- string (default:``None``)
EXAMPLES::
sage: from slabbe.magic_hexagon import solve_magic_hexagon
sage: solve_magic_hexagon() # long time (90s if GLPK, <1s if Gurobi)
[15, 14, 9, 13, 8, 6, 11, 10, 4, 5, 1, 18, 12, 2, 7, 17, 16, 19, 3]
"""
p = MixedIntegerLinearProgram(solver=solver)
x = p.new_variable(binary=True)
A = range(1,20)
# exactly one tile at each position pos
for pos in A:
S = p.sum(x[pos,tile] for tile in A)
name = "one tile at {}".format(pos)
p.add_constraint(S==1, name=name)
# each tile used exactly once
for tile in A:
S = p.sum(x[pos,tile] for pos in A)
name = "tile {} used once".format(tile)
p.add_constraint(S==1, name=name)
lines = [(1,2,3), (4,5,6,7), (8,9,10,11,12), (13,14,15,16), (17,18,19),
(1,4,8), (2,5,9,13), (3,6,10,14,17), (7,11,15,18), (12,16,19),
(8,13,17), (4,9,14,18), (1,5,10,15,19), (2,6,11,16), (3,7,12)]
# the sum is 38 on each line
for line in lines:
S = p.sum(tile*x[pos,tile] for tile in A for pos in line)
name = "sum of line {} must be 38".format(line)
p.add_constraint(S==38, name=name)
p.solve()
soln = p.get_values(x)
nonzero = sorted(key for key in soln if soln[key]!=0)
return [tile for (pos,tile) in nonzero]
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
开发者ID:videlec，项目名称:max_plus，代码行数:88，代码来源:rank.py
示例12: __init__
class hard_EM:
def __init__(self, author_graph, TAU=0.5001, nparts=5, init_partition=None):
self.parts = range(nparts)
self.TAU = TAU
self.author_graph = nx.convert_node_labels_to_integers(author_graph, discard_old_labels=False)
self._lp_init = False
# init hidden vars
if init_partition:
self.partition = init_partition
else:
self._rand_init_partition()
self.m_step()
def _rand_init_partition(self):
slog('Random partitioning with seed: %s' % os.getpid())
random.seed(os.getpid())
self.partition = {}
nparts = len(self.parts)
for a in self.author_graph:
self.partition[a] = randint(0, nparts - 1)
def _init_LP(self):
if self._lp_init:
return
slog('Init LP')
self.lp = MixedIntegerLinearProgram(solver='GLPK', maximization=False)
#self.lp.solver_parameter(backend.glp_simplex_or_intopt, backend.glp_simplex_only) # LP relaxation
self.lp.solver_parameter("iteration_limit", LP_ITERATION_LIMIT)
# self.lp.solver_parameter("timelimit", LP_TIME_LIMIT)
# add constraints once here
# constraints
self.alpha = self.lp.new_variable(dim=2)
beta2 = self.lp.new_variable(dim=2)
beta3 = self.lp.new_variable(dim=3)
# alphas are indicator vars
for a in self.author_graph:
self.lp.add_constraint(sum(self.alpha[a][p] for p in self.parts) == 1)
# beta2 is the sum of beta3s
slog('Init LP - pair constraints')
for a, b in self.author_graph.edges():
if self.author_graph[a][b]['denom'] <= 2:
continue
self.lp.add_constraint(0.5 * sum(beta3[a][b][p] for p in self.parts) - beta2[a][b], min=0, max=0)
for p in self.parts:
self.lp.add_constraint(self.alpha[a][p] - self.alpha[b][p] - beta3[a][b][p], max=0)
self.lp.add_constraint(self.alpha[b][p] - self.alpha[a][p] - beta3[a][b][p], max=0)
# store indiv potential linear function as a dict to improve performance
self.objF_indiv_dict = {}
self.alpha_dict = {}
for a in self.author_graph:
self.alpha_dict[a] = {}
for p in self.parts:
var_id = self.alpha_dict[a][p] = self.alpha[a][p].dict().keys()[0]
self.objF_indiv_dict[var_id] = 0 # init variables coeffs to zero
# pairwise potentials
slog('Obj func - pair potentials')
objF_pair_dict = {}
s = log(1 - self.TAU) - log(self.TAU)
for a, b in self.author_graph.edges():
if self.author_graph[a][b]['denom'] <= 2:
continue
var_id = beta2[a][b].dict().keys()[0]
objF_pair_dict[var_id] = -self.author_graph[a][b]['weight'] * s
self.objF_pair = self.lp(objF_pair_dict)
self._lp_init = True
slog('Init LP Done')
def log_phi(self, a, p):
author = self.author_graph.node[a]
th = self.theta[p]
res = th['logPr']
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')
#.........这里部分代码省略.........
def steiner_tree(self,G, weighted = False):
"""
Returns a tree of minimum weight connecting the given
set of vertices.
+
+ Definition :
+
+ Computing a minimum spanning tree in a graph can be done in `n
+ \log(n)` time (and in linear time if all weights are
+ equal). On the other hand, if one is given a large (possibly
+ weighted) graph and a subset of its vertices, it is NP-Hard to
+ find a tree of minimum weight connecting the given set of
+ vertices, which is then called a Steiner Tree.
+
+ `Wikipedia article on Steiner Trees
+ <http://en.wikipedia.org/wiki/Steiner_tree_problem>`_.
+
+ INPUT:
+
+ - ``vertices`` -- the vertices to be connected by the Steiner
+ Tree.
+
+ - ``weighted`` (boolean) -- Whether to consider the graph as
+ weighted, and use each edge's label as a weight, considering
+ ``None`` as a weight of `1`. If ``weighted=False`` (default)
+ all edges are considered to have a weight of `1`.
+
+ .. NOTE::
+
+ * This problem being defined on undirected graphs, the
+ orientation is not considered if the current graph is
+ actually a digraph.
+
+ * The graph is assumed not to have multiple edges.
+
+ ALGORITHM:
+
+ Solved through Linear Programming.
+
+ COMPLEXITY:
+
+ NP-Hard.
+
+ Note that this algorithm first checks whether the given
+ set of vertices induces a connected graph, returning one of its
+ spanning trees if ``weighted`` is set to ``False``, and thus
+ 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)
#.........这里部分代码省略.........
def OA_and_oval(q):
r"""
Return a `OA(q+1,q)` whose blocks contains `\leq 2` zeroes in the last `q`
columns.
This `OA` is build from a projective plane of order `q`, in which there
exists an oval `O` of size `q+1` (i.e. a set of `q+1` points no three of which
are [colinear/contained in a common set of the projective plane]).
Removing an element `x\in O` and all sets that contain it, we obtain a
`TD(q+1,q)` in which `O` intersects all columns except one. As `O` is an
oval, no block of the `TD` intersects it more than twice.
INPUT:
- ``q`` -- a prime power
.. NOTE::
This function is called by :func:`construction_3_6`, an
implementation of Construction 3.6 from [AC07]_.
EXAMPLES::
sage: from sage.combinat.designs.orthogonal_arrays_recursive import OA_and_oval
sage: _ = OA_and_oval
"""
from sage.rings.arith import is_prime_power
from sage.combinat.designs.block_design import projective_plane
from orthogonal_arrays import OA_relabel
assert is_prime_power(q)
B = projective_plane(q, check=False)
# We compute the oval with a linear program
from sage.numerical.mip import MixedIntegerLinearProgram
p = MixedIntegerLinearProgram()
b = p.new_variable(binary=True)
V = B.ground_set()
p.add_constraint(p.sum([b[i] for i in V]) == q+1)
for bl in B:
p.add_constraint(p.sum([b[i] for i in bl]) <= 2)
p.solve()
b = p.get_values(b)
oval = [x for x,i in b.items() if i]
assert len(oval) == q+1
# We remove one element from the oval
x = oval.pop()
oval.sort()
# We build the TD by relabelling the point set, and removing those which
# contain x.
r = {}
B = list(B)
# (this is to make sure that the first set containing x in B is the one
# which contains no other oval point)
B.sort(key=lambda b:int(any([xx in oval for xx in b])))
BB = []
for b in B:
if x in b:
for xx in b:
if xx == x:
continue
r[xx] = len(r)
else:
BB.append(b)
assert len(r) == (q+1)*q # all points except x have an image
assert len(set(r.values())) == len(r) # the images are different
# Relabelling/sorting the blocks and the oval
BB = [[r[xx] for xx in b] for b in BB]
oval = [r[xx] for xx in oval]
for b in BB:
b.sort()
oval.sort()
# Turning the TD into an OA
BB = [[xx%q for xx in b] for b in BB]
oval = [xx%q for xx in oval]
assert len(oval) == q
# We relabel the "oval" as relabelled as [0,...,0]
OA = OA_relabel(BB+([[0]+oval]),q+1,q,blocks=[[0]+oval])
OA = [[(x+1)%q for x in B] for B in OA]
OA.remove([0]*(q+1))
assert all(sum([xx == 0 for xx in b[1:]]) <= 2 for b in OA)
return OA
#.........这里部分代码省略.........
[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\"")
客服电话
APP下载
官方微信
请发表评论