ASKSAGE: Sage Q&A Forum - RSS feedhttps://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Wed, 02 Sep 2020 18:10:13 +0200solving a matrix equation in the integers.https://ask.sagemath.org/question/53262/solving-a-matrix-equation-in-the-integers/ I'm trying to use sage to solve an expression of the form $xS = y$ for a fixed integer matrix $S$ and a fixed integer vector $y$. I need my solution vector $x$ to also have integer coefficients.
How I'm currently doing this is by using the mixed integer linear programming tool and my code looks something like this.
p = MixedIntegerLinearProgram(maximization = True, solver = "GLPK")
x = p.new_variable(integer = True, nonnegative = True)
p.add_constraint( x*S == y) #here the vector can be anything
p.add_constraint( x[25] <= 2020 )
p.set_objective( x[25] )
p.solve()
p.get_values(x)
The issue I have here is that the objective and the constraint are arbitrary and were just picked so that the space that the MixedIntegerLinearProgram function is trying to optimize over has a "feasable solution" (sage's words not mine). Is there a better way to find a solution to $xS = y$ in the integers (i.e a different function or package supported by sage) and if not, is there a way to have my constarint be to make the norm of the vector as small as possible (so that I'm not just artificially picking an entry to optimize my search over)?
JRHalesWed, 02 Sep 2020 18:10:13 +0200https://ask.sagemath.org/question/53262/Sage stalls during computationhttps://ask.sagemath.org/question/48270/sage-stalls-during-computation/I was using SageMath 8.8's [feedback_edge_set](http://doc.sagemath.org/html/en/reference/graphs/sage/graphs/digraph.html#sage.graphs.digraph.DiGraph.feedback_edge_set) function on a large dataset. The calculation would take a couple of days to finish.
Unfortunately, Sage would keep stalling every few hours. The CPU usage of the process would simply drop to zero, and it would stop doing anything. I had to kill the process (a simple `kill` worked, no `kill -9` necessary).
Has anyone else experienced this? Is it a known problem? Is there a workaround? It is very frustrating to have to keep manually killing and restarting the process. This problem also prevents me from running a calculation overnight: in the morning I might find that it stopped merely an hour after I started it (if I'm unlucky).
----
More information:
- I was using Sage on Linux
- This function uses [MixedIntegerLinearProgram](http://doc.sagemath.org/html/en/reference/numerical/sage/numerical/mip.html#sage.numerical.mip.MixedIntegerLinearProgram). I used the default GLPK (didn't change settings)
- I used `alarm()` to set a time limit on each computation.
- I always believed GLPK [not to be interruptible](https://lists.gnu.org/archive/html/help-glpk/2006-04/msg00059.html), so I am not even sure which interruption (either manual or alarm) work on Sage with this function. This may be relevant.
- The memory usage of the process stayed low. It did not stall because the memory got filled and the machine started swapping.
----
Minimal example:
Put this in a source file:
import datetime
i=1
while True:
alarm(60)
try:
dg=digraphs.RandomDirectedGNM(300,750);
fes = dg.feedback_edge_set()
print(len(fes))
except AlarmInterrupt:
print("timeout")
cancel_alarm()
print(i)
print(str(datetime.datetime.now()))
print("\n")
i = i+1
Run it as `sage sourcefile.sage`.
In my last test, it took 291 iterations and ~4-5 hours before it stalled.
It appears to stall only when GLPK is used as the MIP solver (i.e. `feedback_edge_set(solver='GLPK')`, which is also the default). If I use `solver='Coin'` (which needs to be installed first using `sage -i cbc`) then it does not stall. However, the calculation is several times slower with Coin than with GLPK.SzabolcsFri, 11 Oct 2019 09:27:37 +0200https://ask.sagemath.org/question/48270/Submitting code to Sagehttps://ask.sagemath.org/question/47462/submitting-code-to-sage/ I have written a user-friendly front-end to MixedIntegerLinearProgram (aided by a suggestion of rburing). The code is listed below. I'd be happy to contribute this to the sage code base, but (not being a developer) I'm not sure how.
I have tried to follow the conventions that I've found, but I'd welcome corrections and suggestions of where to go from here.
Thanks,
Mike
r"""
Print the solution to a mixed integer linear program.
Variables are assumed real unless specified as integer,
and all variables are assumed to be nonnegative.
EXAMPLES::
var('x1 x2')
maximize(20*x1 + 10*x2, {3*x1 + x2 <= 1300, x1 + 2*x2 <= 600, x2 <= 250})
# Z = 9000.0, x1 = 400.0, x2 = 100.0
var('s h a u')
maximize(0.05*s + 0.08*h + 0.10*a + 0.13*u, {a <= 0.30*(h+a),
s <= 300, u <= s, u <= 0.20*(h+a+u), s+h+a+u <= 2000})
# Z = 174.4, s = 300.0, u = 300.0, h = 980.0, a = 420.0
var('a b c d', domain='integer')
maximize(5*a + 7*b + 2*c + 10*d,
{2*a + 4*b + 7*c + 10*d <= 15, a <= 1, b <= 1, c <= 1, d <=1 })
# Z = 17.0, a = 0, b = 1, c = 0, d = 1 (Note that integer variables <=1 are binary.)
var('x y')
minimize(x + y, {x + 2*y >= 7, 2*x + y >= 6})
# Z = 4.33333333333, x = 1.66666666667, y = 2.66666666667
var('x')
var('y', domain='integer')
minimize(x + y, {x + 2*y >= 7, 2*x + y >= 6})
# Z = 4.5, x = 1.5, y = 3
var('x y', domain='integer')
minimize(x + y, {x + 2*y >= 7, 2*x + y >= 6})
# Z = 5.0, x = 2, y = 3
var('x y', domain='integer')
maximize(x + y, {x + 2*y >= 7, 2*x + y >= 6})
# GLPK: Objective is unbounded
var('x y', domain='integer')
maximize(x + y, {x + 2*y >= 7, 2*x + y <= 3})
# GLPK: Problem has no feasible solution
AUTHORS:
- Michael Miller (2019-Aug-11): initial version
"""
# ****************************************************************************
# Copyright (C) 2019 Michael Miller
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 2 of the License, or
# (at your option) any later version.
# https://www.gnu.org/licenses/
# ****************************************************************************
from sage.numerical.mip import MIPSolverException
def maximize(objective, constraints):
maxmin(objective, constraints, true)
def minimize(objective, constraints):
maxmin(objective, constraints, false)
def maxmin(objective, constraints, flag):
# Create a set of the original variables
variables = set(objective.variables())
for c in constraints:
variables.update(c.variables())
integer_variables = [v for v in variables if v.is_integer()]
real_variables = [v for v in variables if not v.is_integer()]
# Create the MILP variables
p = MixedIntegerLinearProgram(maximization=flag)
MILP_integer_variables = p.new_variable(integer=True, nonnegative=True)
MILP_real_variables = p.new_variable(real=True, nonnegative=True)
# Substitute the MILP variables for the original variables
# (Inconveniently, the built-in subs fails with a TypeError)
def Subs(expr):
const = RDF(expr.subs({v:0 for v in variables})) # the constant term
sum_integer = sum(expr.coefficient(v) * MILP_integer_variables[v] for v in integer_variables)
sum_real = sum(expr.coefficient(v) * MILP_real_variables[v] for v in real_variables)
return sum_real + sum_integer + const
objective = Subs(objective)
constraints = [c.operator()(Subs(c.lhs()), Subs(c.rhs())) for c in constraints]
# Set up the MILP problem
p.set_objective(objective)
for c in constraints:
p.add_constraint(c)
# Solve the MILP problem and print the results
try:
Z = p.solve()
print "Z =", Z
for v in integer_variables:
print v, "=", int(p.get_values(MILP_integer_variables[v]))
for v in real_variables:
print v, "=", p.get_values(MILP_real_variables[v])
print
except MIPSolverException as msg:
if str(msg)=="GLPK: The LP (relaxation) problem has no dual feasible solution":
print "GLPK: Objective is unbounded"
print
else:
print str(msg)
print
millermjSun, 11 Aug 2019 21:21:28 +0200https://ask.sagemath.org/question/47462/Trouble with mixed integer linear programminghttps://ask.sagemath.org/question/41458/trouble-with-mixed-integer-linear-programming/I'm having trouble with sage's mixed linear programming tool.
I entered the following code trying to maximize a variable with the given constraint that multiplication by a matrix would result in the zero vector.
[picture of inputs is here](https://drive.google.com/file/d/1guE79ATQHZRol0D1PCP9sHpROLi1zTpY/view?usp=sharing)
sage: p = MixedIntegerLinearProgram(maximization = True, solver = "GLPK")
sage: x = p.new_variable(integer = True, nonnegative = True)
sage: p.add_constraint(x*spoof_matrix == 0)
sage: p.set_objective(x[25])
sage: p.solve()
Where spoof_matrix is a 56x34 matrix that has the specific constraints I'm working with.
When I use the p.solve() it gives the error in the picture.
I know that there is at least one non-trivial solution since the test vector seen in the picture solves the constraints.
I don't know why sage isn't solving the optimization problem, and any help would be really great!
Also, if it's possible to not only get the maximum value but also obtain the vector that gives the maximum value, that would be really awesome. Thanks!!
Edit: People have asked for the code for how spoof matrix is created, using these functions the spoof matrix would be given by the following line.
create_spoof_matrix(reduce_n_times(25,generate_possible_factors_list(3),primes_first_n(168),168),168)
def spoof(p,a):
if p*a == 0:
return 1
if p!=1 and a != 1:
value = ((1/p)^(a+1)-1)/(1/p-1)
elif p ==1:
value = (a+1)
elif a ==1:
value = (1+1/p)
return value
def reduce_possible_primes(possible_factors,number_of_primes):
new_prime_list = []
for x in possible_factors:
for p in primes_first_n(number_of_primes):
if x % p == 0:
if p not in new_prime_list:
new_prime_list.append(p)
for x in xrange(0,len(new_prime_list)): #this is just ordering the primes from least to greatest.
for y in xrange(0,len(new_prime_list)):
if new_prime_list[x] < new_prime_list[y]:
cache = new_prime_list[x]
new_prime_list[x] = new_prime_list[y]
new_prime_list[y] = cache
return new_prime_list
def generate_possible_factors_list(size_of_factors):
possible_factors = []
for x in xrange(-1*(10^size_of_factors)+1,10^size_of_factors,2):
if x != -1:
numerator_factor = list(factor(numerator(spoof(x,2))))
if numerator_factor[-1][0] < 10^size_of_factors:
possible_factors.append(x)
return possible_factors
def reduce_possible_factors(possible_factors,possible_primes):
for x in possible_factors:
numerator_factor_bank = list(factor(numerator(spoof(x,2))))
for y in xrange(0,len(numerator_factor_bank)):
if numerator_factor_bank[y][0] not in possible_primes:
if x in possible_factors:
possible_factors.remove(x)
return possible_factors
def reduce_n_times(n,possible_factors, possible_primes, number_of_primes):
x = 0
while x < n:
possible_primes = reduce_possible_primes(possible_factors,number_of_primes)
possible_factors = reduce_possible_factors(possible_factors,possible_primes)
x = x+1
return possible_factors
def create_spoof_matrix(possible_factors,n):
possible_primes = reduce_possible_primes(possible_factors,n)
vector = zero_vector(len(possible_primes))
matrix_list = []
for x in possible_factors:
copy_vector = vector[:]
denominator_list = list(factor(x))
numerator_list = list(factor(numerator(spoof(x,2))))
for y in denominator_list:
copy_vector[possible_primes.index(y[0])] = -2*y[1]
for y in numerator_list:
copy_vector[possible_primes.index(y[0])] = y[1]
matrix_list.append(copy_vector)
return Matrix(matrix_list)JRHalesSat, 10 Mar 2018 01:53:02 +0100https://ask.sagemath.org/question/41458/How to get step by step table in Simplex method in sagemath?https://ask.sagemath.org/question/34649/how-to-get-step-by-step-table-in-simplex-method-in-sagemath/I have used the following commands to solve a linear programing problem in sagemath.
>p = MixedIntegerLinearProgram(maximization=True, solver='GLPK')
>x = p.new_variable(nonnegative=True)
> p.set_objective(12*x[0] + 3*x[1]+x[2])
> p.add_constraint(10*x[0] + 2*x[1]+x[2]<= 100)
>p.add_constraint(7*x[0] +3*x[1]+2*x[2] <= 77)
> p.add_constraint(2*x[0] + 4*x[1]+x[2] <= 80)
> p.show()
> print 'The optimal value:', QQ(p.solve())
It solved and gave the answer. But I want to generate a step by step table to explain to the students.
Is there any possible way to get the step by step table? or should we have to modify the program to do so?daviddglTue, 30 Aug 2016 08:15:37 +0200https://ask.sagemath.org/question/34649/Evaluation of a LinearFunctionhttps://ask.sagemath.org/question/33536/evaluation-of-a-linearfunction/ I would like to substitute values for the variables in a LinearFunction object. For example:
sage: lp.<x> = MixedIntegerLinearProgram()
sage: f = 2*x[0] + x[1]
sage: f
2*x_0 + x_1
sage: type(f)
<type 'sage.numerical.linear_functions.LinearFunction'>
The object f is not callable, and also the following does nothing (even though the Sage online help suggests that it should do something):
sage: f.subs({x[0]:1.0, x[1]:2.0})
2*x_0 + x_1
I can do what I want to do by writing some code to work with the output of f.dict(), but does Sage have a cleaner way to evaluate a linear function?EvanTue, 24 May 2016 15:25:53 +0200https://ask.sagemath.org/question/33536/Get index of a set constrainthttps://ask.sagemath.org/question/32119/get-index-of-a-set-constraint/ Using SageMath I'm trying to implement a method for the solution of a problem where I in an iterative way add a set of constraints, try to solve the problem to check if exist a solution, if not i remove the constraints and restart.
try:
p.solve()
except RuntimeError:
p.remove_constraint(index)
Given a set of constraint that I add, how to get their indexes so that I can remove them?tom123Sat, 09 Jan 2016 13:04:36 +0100https://ask.sagemath.org/question/32119/Binary variable in mixed integer linear programhttps://ask.sagemath.org/question/29626/binary-variable-in-mixed-integer-linear-program/Hi asksage!
I'm using the `MixedIntegerLinearProgram` class to model a mixed integer linear program. I am using binary variables created e.g. by `x = mip.new_variable(binary=True, name="x")`. Now, I am using this variable e.g. by `x[0]` or `x[3,2]`. When I create a constraint using sums over binary variables, this sum is in the integers, as it should be. However, I am wondering if there is another way to easily add such variables together over `GF(2)` such that e.g. `x[0] + x[0]` disappears from the equation?
I ask because my model is too complicated to write by hand, so I am using some auxiliary methods to construct constraints and variables, etc.
I hope you can help me out!mmlauridsenMon, 28 Sep 2015 08:58:51 +0200https://ask.sagemath.org/question/29626/How to stop a MILP before it reaches the optimal solution (using different backends)?https://ask.sagemath.org/question/25909/how-to-stop-a-milp-before-it-reaches-the-optimal-solution-using-different-backends/Sometimes you don't need an optimal solution. A feasible solution is a must, but a decent value of the objective function will do. At least that's better than waiting forever...
Some backends (most?, all?) for the MixedIntegerLinearProgram support different stop criteria. But it took me quite some time to find them. How to tell MixedIntegerLinearProgram not to continue until the optimal solution is found?
PD: I found the answer (at least part of it) before I submitted this question, but I think it's good to ask the question for reference:pangTue, 24 Feb 2015 10:19:05 +0100https://ask.sagemath.org/question/25909/MixedIntegerLinearProgram method .get_values() returns "IndexError: list index out of range"https://ask.sagemath.org/question/11053/mixedintegerlinearprogram-method-get_values-returns-indexerror-list-index-out-of-range/I'm writing an MLP to solve a scheduling problem for a set of jobs. I have several variable types, including family 's_ij', which is a boolean indicating if time i occurs before or after time j. My program works, and returns feasible solutions. I'm trying to see the assigned values for each variable, and while I can see the start time for each job, and some other families of booleans, when I try to display the values for the 's_ij' variables, I get an "IndexError". This is weird, since I'm letting the MILP package do all the work here. Here are relevant sections of the code, as well as the output:
from itertools import product
...
J = range(len(JOBS))
I = range(2*len(JOBS))
IXJ_copy = product(I, J)
...
for pair in IXJ_copy:
i = pair[0]
j = pair[1]
if i == len(I) - 1:
scheduling_mlp.add_constraint(s[(i, j)] <= 0)
else:
scheduling_mlp.add_constraint(M*(1-s[(i, j)]) + b[i+1] - t[j] - EPSILON >= 0)
scheduling_mlp.add_constraint(M*s[(i, j)] + t[j] - b[i+1] >= 0)
scheduling_mlp.add_constraint(M*e[(i, j)] + b[i] - t[j] - T(j) >= 0)
scheduling_mlp.add_constraint(M*(1-e[(i, j)]) + t[j] + T(j) - b[i] - EPSILON >= 0)
scheduling_mlp.add_constraint(M*(1-u[(i, j)]) + s[(i, j)] + e[(i, j)] >= 2)
scheduling_mlp.add_constraint(-M*(u[(i, j)]) + s[(i, j)] + e[(i, j)] <= 1)
...
print("s:")
start = scheduling_mlp.get_values(s)
print(start)
Returns
Traceback (most recent call last):
File "lp_solver.py", line 186, in <module>
start = scheduling_mlp.get_values(s)
File "mip.pyx", line 1161, in sage.numerical.mip.MixedIntegerLinearProgram.get_values (sage/numerical/mip.c:7098)
IndexError: list index out of range
As I said, when I run the code, I get feasible solutions back. I have the case handling in place while defining my constraints to prevent b[i+1] from going out of range. I'm able to get_values and print all other variables. Any ideas what's going wrong here?zbowenWed, 19 Feb 2014 06:45:52 +0100https://ask.sagemath.org/question/11053/Use of constraint_generation in MixedIntegerLinearProgramhttps://ask.sagemath.org/question/10915/use-of-constraint_generation-in-mixedintegerlinearprogram/I would like to solve an LP problem using row generation (aka constraint generation), where I wish to be able to choose how the constraints are generated (this will be done using SAT, or WPMaxSAT through a binary MIP).
Now, in MixedIntegerLinearProgram, there is a boolean option `constraint_generation`. I do not find an explanation of its use and therefore whether it can be used for my stated purpose. So, can I, and if so, how?
Information found, but which did not make th option clear to me:
From the [documentation](http://www.sagemath.org/doc/reference/numerical/sage/numerical/mip.html):
> constraint_generation whether to require the returned solver to support constraint generation (excludes Coin). False by default.
From the [source code documentation](http://https://github.com/sagemath/sage/blob/master/src/sage/numerical/backends/generic_backend.pyx#L983):
> - ``constraint_generation`` (boolean) -- whether the solver returned is to be used for constraint/variable generation. As the interface with Coin does not support constraint/variable generation, setting ``constraint_generation`` to ``False`` ensures that the backend to Coin is not returned when ``solver = None``. This is set to ``False`` by default.
(N.B.: Shouldn't that be: backend to Coin is not returned -> Coin backend may be returned?)Erik QuaeghebeurWed, 15 Jan 2014 09:51:10 +0100https://ask.sagemath.org/question/10915/Output results of sage MixedIntegerLinearProgram solution are not as expectedhttps://ask.sagemath.org/question/10827/output-results-of-sage-mixedintegerlinearprogram-solution-are-not-as-expected/For the problem below
------------------------
<pre><code>#!/usr/bin/sage -python
from sage.all import *
import sys
problem = MixedIntegerLinearProgram(maximization=False, solver = "GLPK")
y = problem.new_variable(real=True, name='y')
constraints_1 = [(0.0, 6769.0),
(8260.0, 6812.0),
(12390.0, 7105.0),
(16520.0, 7448.0),
(20650.0, 7924.0),
(24780.0, 8289.0),
(28910.0, 8630.0),
(33040.0, 9000.0)]
sos_var_1 = problem.new_variable(name='sos2_var')
eq1_1 = 0.0
eq1_2 = 0.0
eq1_3 = 0.0
for i in range(len(constraints_1)):
eq1_1 += sos_var_1[i]*constraints_1[i][0]
eq1_2 += sos_var_1[i]*constraints_1[i][1]
eq1_3 += sos_var_1[i]
problem.add_constraint(eq1_1 == 23450)
problem.add_constraint(eq1_2 == y[0])
problem.add_constraint(eq1_3 == 1)
problem.set_objective(y[0])
print problem.solve()
</code></pre>
--------------------------------
I am expecting the result to be 8136 but I am getting <pre><code>8099.05084746</code></pre>
Can anyone help or explain why this is happening?subbuThu, 12 Dec 2013 16:01:29 +0100https://ask.sagemath.org/question/10827/Why do I get the "unable to find a common ring for all elements" error message?https://ask.sagemath.org/question/10746/why-do-i-get-the-unable-to-find-a-common-ring-for-all-elements-error-message/I was wondering if you can take a look at this cell:
http://sagecell.sagemath.org/?q=ypqhgw
and let me know why I get the error message.
ThanksBehzadSun, 17 Nov 2013 20:46:09 +0100https://ask.sagemath.org/question/10746/