Skip to content

cached best; added first walk back. #4

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
anteprandium wants to merge 15 commits into
base: master
Choose a base branch
from
Open

cached best; added first walk back. #4

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
9efbda9
cached best; added first walk back.
anteprandium Oct 20, 2015
eadb76c
Bug fix
anteprandium Oct 20, 2015
5a4beb8
Added sketch of proof
anteprandium Oct 21, 2015
53fb1e2
fixed bug in ipproblem and started wack back
anteprandium Mar 23, 2016
5735f9d
added doctests
anteprandium Mar 23, 2016
c589c35
trunc1st
anteprandium Mar 8, 2017
d03a0d3
changed eligibility
anteprandium Mar 8, 2017
ca1c532
Removed all splitting stuff and improved speed by a factor of 100
anteprandium Mar 9, 2017
ee4c843
Fix to criterion 2; added self_reduce and changed back eligibility cr…
anteprandium Mar 28, 2017
ca46552
reduction recycles vectors and matrices
anteprandium Mar 31, 2017
453e55e
fixed the criterion
anteprandium Apr 3, 2017
83b24a5
Fixed binomial ordering; added p,m supports to precomputed data
anteprandium Apr 4, 2017
c99cb23
Removed order, and use new pm_split
anteprandium Apr 4, 2017
b576e5d
commit pre-unvector
anteprandium Apr 6, 2017
8ef887e
pre-basiselement
anteprandium Apr 7, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
622 changes: 622 additions & 0 deletions binomials.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

206 changes: 206 additions & 0 deletions binomials.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,206 @@
#!/usr/bin/env python


class BinomialBasis(object):
"""

Fastest way to compute cost: c*v


"""
def __init__(self, A, b,c, u):
super(BinomialBasis, self).__init__()
self.d, self.n = A.dimensions()
self.A = A
self.c = c
self.b = b
self.u = u
self.bins = [] # binomials
self.Abins = [] # A*v for v in binomials
self.bins_psupp = [] # positive support of binomial
self.bins_nsupp = [] # negative support of binomials
self.Abins_psupp = [] # positive support of Av
self.Abins_nsupp = [] # negative support of Av
self.cbins = [] # c*v for v in binomials
self.pairs = [] # pairs still to go.
self.zerox = vector(ZZ,self.n*[0]) # for fast copying.
self.zerob = vector(ZZ, self.d*[0]) # for fast copying.

#
#
for e in identity_matrix(ZZ,self.n).rows():
self.push_binomial(e)

def supports(self, v):
"""
Compute the positive and negative support of v
"""
psupp = []
nsupp = []
for (k, c) in enumerate(v):
if c>0:
psupp.append(k)
elif c<0:
nsupp.append(k)
return psupp, nsupp


def push_binomial(self, v, Av=None):
"""
Add binomial v and associated data into place.

May change sign of v and Av in place.

"""
cost = self.c*v
if Av is None:
av = self.A*v
else:
av = av

if cost<0:
# modify v and av in place.
for k in xrange(self.n):
v[k] = -v[k]
for k in xrange(self.d):
av[k] = -av[k]

if cost == 0:
raise RuntimeError("zero cost vector %s. Deal with this." %v)

p, n = self.supports(v)
ap, an = self.supports(av)
l = len(self.bins)

self.bins.append(v)
self.cbins.append(abs(cost))
self.Abins.append(av)
self.bins_psupp.append(p)
self.bins_nsupp.append(n)
self.Abins_psupp.append(ap)
self.Abins_nsupp.append(an)

if l: # true if this is not the first vector
for i in xrange(l):
self.pairs.append((i,l))

def pop_binomials(i, pairs=True):
self.bins.pop(i)
self.cbins.pop(i)
self.Abins.pop(i)
self.bins_psupp.pop(i)
self.bins_nsupp.pop(i)
self.Abins_psupp.pop(i)
self.Abins_nsupp.pop(i)
if pairs:
for (k,p) in enumerate(self.pairs):
if i in p:
self.pairs.pop(k)

def next_pair(self):
return self.pairs.pop[-1]

def order(self, i,j):
"""
Set i, j in gr-rev-lex order:

vi > vj if c*vi > c*vj or
equality and the last nonzero coordinate of vi-vj is negative.

We sort indices and not vectors so that we don't have to compute
c*v and c*v more than once.
"""
if self.cbins[i]>self.cbins[j]:
return (i,j)
elif self.cbins[i]<self.cbins[j]:
return (j,i)
else: # equal cost :(
for k in xrange(self.n):
l = self.n-k
if self.bins[l] > self.bins[l] :
return (j,i)
elif self.bins[l] < self.bins[l] :
return (i,j)
# The vectors are equal. This should never happen:
raise RuntimeError("Vectors %s and %s are equal (%s)" % (i,j,self.bins[i]))


def criterion_1(self, i,j):
if (
bool(set.intersect(self.Abins_psupp[i], self.Abins_psupp[j])) or
bool(set.intersect(self.bins_nsupp[i], self.bins_nsupp[j])) or
bool(set.intersect(self.bins_psupp[i], self.bins_psupp[j]))
):
return False
else:
return True

def feasible(self, v):
"""docstring for feasible"""
if self.less_or_equal_pp(v,self.b) and self.less_or_equal_mp(v,self.b):
avp, avm = pm_split(v)
return self.less_or_equal_pp(self.A*avp, self.b) and self.less_or_equal_pp(self.A*avm, self.b)
return False

def pm_split(v):
""" slightly faster than pm_split2"""
up = copy(self.zerox)
um = copy(self.zerox)
for i,a in enumerate(v):
if a>0:
up[i]=a
elif a<0:
um[i]=-a
return up, um




def less_or_equal(self, a,b):
"True if and only if a is less or equal than b."
for (i,ai) in enumerate(a):
if ai>b[i]:
return False
return True

def less_or_equal_pp(self, a,b):
"True if and only if a^+ is less or equal than b^+."
for (i,ai) in enumerate(a):
if ai>0 and ai>b[i]:
return False
return True

def less_or_equal_mp(self, a,b):
"True if and only if a^+ is less or equal than b^+."
for (i,ai) in enumerate(a):
if ai<0 and -ai>b[i]:
return False
return True






def bu_buch(self):
"""docstring for fname"""
skip_c1 = 0


while self.pairs:
i,j = self.next_pair()
if criterion_1(i,j):
skip_c1 +=1
continue

i,j = self.order(i,j)
# fast r = ui-uj
r = copy(g.zerox)
for l in xrange(g.n):
r[l] = self.bins[i][l]-self.bins[j][l]
# end of fast r=ui-uj
if self.feasible(v):




2 changes: 2 additions & 0 deletions ejemploucha/ejemplosoto.cost
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
1 11
-252 -537 -965 -237 -724 0 0 0 0 0 0
18 changes: 18 additions & 0 deletions ejemploucha/ejemplosoto.gro
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,18 @@
17 11
-1 0 0 0 0 1 0 0 0 0 999
-1 0 0 1 0 1 0 0 -1 0 583
0 -1 0 -1 1 0 1 0 1 -1 -339
0 -1 0 0 0 0 1 0 0 0 18
0 -1 0 1 0 0 1 0 -1 0 -398
0 0 -1 0 0 0 0 1 0 0 625
0 0 -1 0 1 0 0 1 0 -1 -148
0 0 -1 1 0 0 0 1 -1 0 209
0 0 0 -1 0 0 0 0 1 0 416
0 0 0 0 -1 0 0 0 0 1 773
0 0 0 1 -1 0 0 0 -1 1 357
0 1 -1 0 0 0 -1 1 0 0 607
0 1 -1 1 0 0 -1 1 -1 0 191
0 1 0 0 -1 0 -1 0 0 1 755
1 -1 0 0 0 -1 1 0 0 0 -981
1 0 -1 0 0 -1 0 1 0 0 -374
1 0 0 0 -1 -1 0 0 0 1 -226
7 changes: 7 additions & 0 deletions ejemploucha/ejemplosoto.mat
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,7 @@
6 11
999 18 625 416 773 0 0 0 0 0 1
1 0 0 0 0 1 0 0 0 0 0
0 1 0 0 0 0 1 0 0 0 0
0 0 1 0 0 0 0 1 0 0 0
0 0 0 1 0 0 0 0 1 0 0
0 0 0 0 1 0 0 0 0 1 0
2 changes: 2 additions & 0 deletions ejemploucha/ejemplosoto.zsol
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
1 11
0 0 0 0 0 1 1 1 1 1 1415
82 changes: 80 additions & 2 deletions ipproblem.py
100644 → 100755
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
# coding=utf-8
#!/usr/bin/env sage -python
# encoding: utf-8

"""
Algorithms for computing test sets of Integer Programming problems.
Expand Down Expand Up @@ -49,6 +50,7 @@ def __init__(self, c, A, b, u=None):
self.b = vector(ZZ,b)
self.c = vector(ZZ,c)
self.minimal=None
self.best=None
self.non_reducible=None
self.rows=self.A.nrows()
self.cols=self.A.ncols()
Expand Down Expand Up @@ -334,12 +336,14 @@ def test_set_non_reducibles(self):
def walk_to_best(self,s, path=False):
"""Given a feasible solution $s$, walk from $s$ to the optimum. Return the optimum, or, if `path` is `True`, the full path from `s` to the optimum.
"""
if self.best:
return self.best
s_0=vector(ZZ,s)
if not self.is_feasible(s_0):
return False
T=self.minimal or self.non_reducible
if not T:
self.minimal()
self.minimal_test_set()
T=self.minimal
F=[self.is_feasible(s_0+t) for t in T]
path_l=[s_0]
Expand All @@ -362,6 +366,7 @@ def walk_to_best(self,s, path=False):
cont=any(F)
else:
cont=False
self.best=s_0
if path:
return path_l
else:
Expand Down Expand Up @@ -389,3 +394,76 @@ def __repr__(self):
b=%s,
u=%s."""
return s % (self.c, self.A, self.b, self.u)


def walkback(self, region=lambda x: True):
"""
Start a walk back from the optimus of the regular problem
to an optimus which is also within the region determined by
the function `region` (Might be non-linear).

>>> P=IPProblem(c=[3, 4, 3, 6, 3],A=[[3, 2, 1, 5, 7]], b=[4], u=[1,1,1,1,1] )
>>> P.minimal_test_set()
[(1, 0, 0, 0, 0), (0, 1, 0, 0, 0), (0, 0, 1, 0, 0), (-1, 1, 0, 0, 0), (1, 0, -1, 0, 0), (0, 1, -1, 0, 0)]
>>> P.walk_to_best(P.get_feasible())
(0, 1, 1, 0, 0)
>>> P.walkback(region=lambda x: x[0]**2+x[1]**2+x[2]**2<2)
((0, 1, 0, 0, 0), 4)


!!! Nota:
Para demostrar que este algoritmo siempre calcula el óptimo con la restricción
adicional de la región, hay que argumentarlo. Creo que
una línea argumental positiva es la siguiente: Si el “esqueleto” que resulta
de considerar el test set da un grafo dirigido con fin en el óptimo, el grafo
con direcciones inversas empezando en el óptimo, con el algoritmo de
travesía de grafos en profundidad y con las aristas ordenadas en orden creciente,
recorre cada punto desde el coste mayor hacia abajo, y se queda con el máximo
cada vez que entra en la región determinada por `region`.
"""
# Sort a minimal test set by increasing cost.
T=sorted(self.minimal_test_set(), key=self.cost)
first_point=self.walk_to_best(self.get_feasible(), False)
visited=[]
current_best=False
current_cost=False
S=[first_point] # the first point is always feasible.
while S:
v=S.pop(0) # Get first in queue
verbose("|S|=%d"%len(S), level=2)
c=self.cost(v)

# Cross out from future examinations
if v in visited:
break
else:
visited.append(v)
# see if it's better than what we have so far
# if self.is_feasible(v) and region(v):
if region(v): # we can skip feasibility, since we only include feasibles.
if current_cost is None or c>current_cost:
current_best,current_cost=v,c
verbose("current best=%s, current_cost=%s"%(current_best, current_cost), level=1)
# Now, look at sons.
for d in T:
w = v-d # You'll be worsening things here.
new_c = self.cost(w)
if (self.is_feasible(w) and
not w in S and
not w in visited and
(current_cost is None or new_c > current_cost)):
verbose("appending %s"%w, level=2)
S.append(w)

if current_cost:
return current_best, current_cost
else:
return None


# if __name__ == '__main__':
# import doctest
# doctest.testmod()



Binary file added obj3x7.sobj
View file
Open in desktop
Binary file not shown.
Binary file added objeto.sobj
View file
Open in desktop
Binary file not shown.
Loading