diff --git a/binomials.ipynb b/binomials.ipynb new file mode 100644 index 0000000..40531e5 --- /dev/null +++ b/binomials.ipynb @@ -0,0 +1,622 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "\n", + "longitud = 50\n", + "numero = 50000\n", + "altura = 7\n", + "\n", + "vvv = [random_vector(ZZ, longitud, x=0, y=20) for i in range(numero)]\n", + "www = vvv[::-1]\n", + "zero = vector(ZZ,longitud*[0])" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def minux(v):\n", + " for i in xrange(len(v)):\n", + " v[i] *=-1" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1 loop, best of 3: 443 ms per loop\n" + ] + } + ], + "source": [ + "%%timeit\n", + "for v in vvv:\n", + " v=-v" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1 loop, best of 3: 2.01 s per loop\n" + ] + } + ], + "source": [ + "%%timeit\n", + "for v in vvv:\n", + " minux(v)" + ] + }, + { + "cell_type": "code", + "execution_count": 30, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1 loop, best of 3: 591 ms per loop\n" + ] + } + ], + "source": [ + "%%timeit\n", + "for v in vvv:\n", + " v *= -1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def pm_split(v):\n", + " \"\"\" slightly faster than pm_split2\"\"\"\n", + " l = len(v)\n", + " up = l*[0]\n", + " um = l*[0]\n", + " for i,a in enumerate(v):\n", + " if a>0:\n", + " up[i]=a\n", + " else:\n", + " um[i]=-a\n", + " return vector(ZZ,up), vector(ZZ,um)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def pm_split(v):\n", + " \"\"\" slightly faster than pm_split2\"\"\"\n", + " up = copy(zero)\n", + " um = copy(zero)\n", + " for i,a in enumerate(v):\n", + " if a>0:\n", + " up[i]=a\n", + " else:\n", + " um[i]=-a\n", + " return up, um\n" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "\n", + "def pm_split(v):\n", + " \"\"\" slightly faster than pm_split2\"\"\"\n", + " up = copy(zero)\n", + " um = copy(zero)\n", + " for i,a in enumerate(v):\n", + " if a>0:\n", + " up[i]=a\n", + " elif a<0:\n", + " um[i]=-a\n", + " return up, um\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1 loop, best of 3: 3.3 s per loop\n" + ] + } + ], + "source": [ + "%%timeit\n", + "for v in vvv:\n", + " pm_split(v)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "A = random_matrix(ZZ, 7, 25, x=1, y=15)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "A.apply_map?" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def less_or_equal(a,b):\n", + " \"True if and only if a is less or equal than b. Faster than zipping lists by a factor of 100\"\n", + " for i in xrange(len(a)):\n", + " if a[i]>b[i]:\n", + " return False\n", + " return True\n" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def less_or_equal2(a,b):\n", + " \"True if and only if a is less or equal than b. Faster than zipping lists by a factor of 100\"\n", + " for i,a in enumerate(a):\n", + " if a>b[i]:\n", + " return False\n", + " return True\n" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "10 loops, best of 3: 64.2 ms per loop\n" + ] + } + ], + "source": [ + "%%timeit\n", + "for k in xrange(numero):\n", + " less_or_equal(vvv[k],www[k])" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "10 loops, best of 3: 55.5 ms per loop\n" + ] + } + ], + "source": [ + "%%timeit\n", + "for k in xrange(numero):\n", + " less_or_equal2(vvv[k],www[k])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def supports(v):\n", + " \"\"\"\n", + " Compute the positive and negative support of v\n", + " \"\"\"\n", + " psupp = []\n", + " nsupp = []\n", + " for (k, c) in enumerate(v):\n", + " if c>0:\n", + " psupp.append(k)\n", + " elif c<0:\n", + " nsupp.append(k)\n", + " return psupp, nsupp" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1 loop, best of 3: 744 ms per loop\n" + ] + } + ], + "source": [ + "%%timeit\n", + "for i in xrange(numero):\n", + " supports(vvv[i])" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1 loop, best of 3: 1.36 s per loop\n" + ] + } + ], + "source": [ + "%%timeit\n", + "for i in xrange(numero):\n", + " copy2(vvv[i],www[i])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 109, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "c = random_vector(ZZ,longitud, x=1, y=15)\n", + "A = random_matrix(ZZ, altura, longitud, x=0, y=50)\n", + "b = random_vector(ZZ, altura, x=90, y=150)\n", + "u = random_vector(ZZ, longitud, x=1, y=14)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 110, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "load(\"binomials.py\")" + ] + }, + { + "cell_type": "code", + "execution_count": 111, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "g = BinomialBasis(A,b,c,u)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Cost functions" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# Cost function. Best is cost1\n", + "\n", + "from itertools import izip\n", + "\n", + "def cost1(v):\n", + " return c*v\n", + "\n", + "def cost2(v):\n", + " return sum(c[i]*v[i] for i in xrange(longitud))\n", + "\n", + "def cost3(v):\n", + " return sum(ci*vi for ci, vi in izip(c,v))\n", + "\n", + "def cost4(v):\n", + " return sum(ci*vi for ci, vi in zip(c,v))\n", + "\n", + "def cost5(v):\n", + " return sum(c[i]*v[i] for i in range(longitud))\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 44, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "100 loops, best of 3: 8.95 ms per loop\n" + ] + } + ], + "source": [ + "%%timeit \n", + "for v in zzz: \n", + " cost1(v)" + ] + }, + { + "cell_type": "code", + "execution_count": 45, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "10 loops, best of 3: 79.1 ms per loop\n" + ] + } + ], + "source": [ + "%%timeit \n", + "for v in zzz: \n", + " cost2(v)" + ] + }, + { + "cell_type": "code", + "execution_count": 46, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "10 loops, best of 3: 83.3 ms per loop\n" + ] + } + ], + "source": [ + "%%timeit \n", + "for v in zzz: \n", + " cost3(v)" + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "10 loops, best of 3: 98 ms per loop\n" + ] + } + ], + "source": [ + "%%timeit \n", + "for v in zzz: \n", + " cost4(v)" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "10 loops, best of 3: 93.2 ms per loop\n" + ] + } + ], + "source": [ + "%%timeit \n", + "for v in zzz: \n", + " cost5(v)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 7.5.1", + "language": "", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.13" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/binomials.py b/binomials.py new file mode 100755 index 0000000..f47e386 --- /dev/null +++ b/binomials.py @@ -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.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): + + + + diff --git a/ejemploucha/ejemplosoto.cost b/ejemploucha/ejemplosoto.cost new file mode 100644 index 0000000..d030ada --- /dev/null +++ b/ejemploucha/ejemplosoto.cost @@ -0,0 +1,2 @@ +1 11 +-252 -537 -965 -237 -724 0 0 0 0 0 0 diff --git a/ejemploucha/ejemplosoto.gro b/ejemploucha/ejemplosoto.gro new file mode 100644 index 0000000..60500ae --- /dev/null +++ b/ejemploucha/ejemplosoto.gro @@ -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 diff --git a/ejemploucha/ejemplosoto.mat b/ejemploucha/ejemplosoto.mat new file mode 100644 index 0000000..827974c --- /dev/null +++ b/ejemploucha/ejemplosoto.mat @@ -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 diff --git a/ejemploucha/ejemplosoto.zsol b/ejemploucha/ejemplosoto.zsol new file mode 100644 index 0000000..87814b9 --- /dev/null +++ b/ejemploucha/ejemplosoto.zsol @@ -0,0 +1,2 @@ +1 11 +0 0 0 0 0 1 1 1 1 1 1415 diff --git a/ipproblem.py b/ipproblem.py old mode 100644 new mode 100755 index cbb1b02..92c6fde --- a/ipproblem.py +++ b/ipproblem.py @@ -1,4 +1,5 @@ -# coding=utf-8 +#!/usr/bin/env sage -python +# encoding: utf-8 """ Algorithms for computing test sets of Integer Programming problems. @@ -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() @@ -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] @@ -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: @@ -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() + + + \ No newline at end of file diff --git a/obj3x7.sobj b/obj3x7.sobj new file mode 100644 index 0000000..936a659 Binary files /dev/null and b/obj3x7.sobj differ diff --git a/objeto.sobj b/objeto.sobj new file mode 100644 index 0000000..969946b Binary files /dev/null and b/objeto.sobj differ diff --git a/optimus_check.ipynb b/optimus_check.ipynb index b76ae78..9bfd142 100644 --- a/optimus_check.ipynb +++ b/optimus_check.ipynb @@ -11,18 +11,18 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 2, "metadata": { "collapsed": false }, "outputs": [], "source": [ - "#%load_ext sage" + "%load_ext sage" ] }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 13, "metadata": { "collapsed": false }, @@ -272,7 +272,49 @@ }, { "cell_type": "code", - "execution_count": 40, + "execution_count": 15, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "set_verbose(0)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%load_ext sage" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%run ipproblem.py" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 3, "metadata": { "collapsed": false }, @@ -284,10 +326,10 @@ "Integer programming problem\n", " max\\{c.x : A x <= b, 0<=x<=u\\}\n", "with\n", - "c=(78, 64, 58, 75, 63, 51, 79, 66, 85, 36),\n", + "c=(66, 84, 70, 41, 47, 96, 59, 26, 66, 29),\n", "A=\n", - "[24 20 21 13 10 28 26 12 14 29],\n", - "b=(81),\n", + "[28 18 16 19 18 10 16 24 10 14],\n", + "b=(73),\n", "u=(1, 1, 1, 1, 1, 1, 1, 1, 1, 1).\n" ] } @@ -304,29 +346,26 @@ }, { "cell_type": "code", - "execution_count": 41, + "execution_count": 4, "metadata": { "collapsed": false }, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - " " - ] - } - ], + "outputs": [], "source": [ - "%%prun \n", + "#%%prun \n", "#%%time\n", "\n", + "P=IPProblem(c=(80, 39, 29, 61, 66, 74, 52, 29, 40, 62), \n", + " A=[[14, 20, 27, 10, 19, 14, 23, 18, 23, 20]], \n", + " b=[66],\n", + " u=(1, 1, 1, 1, 1, 1, 1, 1, 1, 1))\n", + "\n", "L=P.minimal_test_set()\n" ] }, { "cell_type": "code", - "execution_count": 42, + "execution_count": 5, "metadata": { "collapsed": false }, @@ -334,26 +373,40 @@ { "data": { "text/plain": [ - "97" + "((1, 0, 0, 1, 1, 1, 0, 0, 0, 0), 281)" ] }, - "execution_count": 42, + "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "len(L)" + "x=P.walk_to_best(P.get_feasible())\n", + "x, P.cost(x)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": { "collapsed": false }, - "outputs": [], - "source": [] + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "((1, 0, 0, 0, 0, 0, 0, 0, 0, 0), 80)\n" + ] + } + ], + "source": [ + "set_verbose(0)\n", + "f=P.walkback(region=lambda v: add([v[i]**2 for i in range(len(v))])<2)\n", + "print f\n", + "set_verbose(0)" + ] }, { "cell_type": "code", @@ -367,9 +420,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Sage 6.6", - "language": "", - "name": "sage_6_6" + "display_name": "Python 2", + "language": "python", + "name": "python2" }, "language_info": { "codemirror_mode": { @@ -381,7 +434,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython2", - "version": "2.7.8" + "version": "2.7.9" } }, "nbformat": 4, diff --git a/optimus_prime.ipynb b/optimus_prime.ipynb index 4c04b89..d3aa9c0 100644 --- a/optimus_prime.ipynb +++ b/optimus_prime.ipynb @@ -1,747 +1,791 @@ { - "metadata": { - "name": "", - "signature": "sha256:5b324d41b9a26a91cd4e9cc161937852cf2c2de1b7f6d9516918e9c75afa2729" - }, - "nbformat": 3, - "nbformat_minor": 0, - "worksheets": [ + "cells": [ { - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Algoritmos de U-W-Z para G-S-U.\n", - "\n", - "Esta hoja contiene los resultados de la implementaci\u00f3n de los algoritmos del art\u00edculo (Urbaniak, Weismantel, Ziegler 1997). De momento no pod\u00e9is ejecutarla, pero s\u00ed pod\u00e9is leerla, para ver el tipo de cosas que ya podemos calcular. \n", - "\n", - "El n\u00facleo de lo que he hecho es un programa en SAGE de unas 200 l\u00edneas. Cuando nos veamos despu\u00e9s de pascua, os dir\u00e9 c\u00f3mo ten\u00e9is que llamarlo desde SAGE.\n", - "\n", - "\n" - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "#%load_ext sage" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 3 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "%run ipproblem" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 4 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 4 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Un problema de IP se declara de la siguiente forma:" - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 4 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P1=IPProblem(c=[1,2,1],\n", - " A=[[1,0,0],[0,0,1]], \n", - " b=[10,10] )" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 5 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "O bien, simplemente como" - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P1=IPProblem([1,2,1],\n", - " [[1,0,0],[0,0,1]], \n", - " [10,10] )" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 6 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Para comprobar los datos que hemos metido, podemos \"imprimir\" el problema" - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "print P1" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "output_type": "stream", - "stream": "stdout", - "text": [ - "Integer programming problem\n", - " max\\{c.x : A x <= b, 0<=x<=u\\}\n", - "with\n", - "c=(1, 2, 1),\n", - "A=\n", - "[1 0 0]\n", - "[0 0 1],\n", - "b=(10, 10),\n", - "u=(10, 0, 10).\n" - ] - } - ], - "prompt_number": 7 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Algunas funciones \u00fatiles:\n", - "\n", - "* `is_feasible(v)` determina si el vector v es \u201cfeasible\u201d.\n", - "* `order(v,w)` determina si los vectores est\u00e1n ordenados para el orden $\\prec_c$ del art\u00edculo\n", - "* `pm_split(v)` devuelve un par de vectores, correspondientes a $v^+$ y $v^-$ del art\u00edculo.\n", - "* `succ(v)` devuelve el vector $v^\\prec$ del art\u00edculo\n", - "\n", - "Hay m\u00e1s funciones, pero no voy a dar ejemplos de todas. " - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P1.is_feasible(vector(ZZ,[0,0,0]))" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 8, - "text": [ - "True" - ] - } - ], - "prompt_number": 8 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P1.order(vector(ZZ,[0,2,0]),vector(ZZ,[13,0,1]))" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 9, - "text": [ - "-1" - ] - } - ], - "prompt_number": 9 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Para las comprobaciones, uso el IP del Ejemplo 2.3 del art\u00edculo" - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P=IPProblem(c=[1,3,2], A=[[1,2,3]], b=[3], u=[1,1,1])\n", - "print P" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "output_type": "stream", - "stream": "stdout", - "text": [ - "Integer programming problem\n", - " max\\{c.x : A x <= b, 0<=x<=u\\}\n", - "with\n", - "c=(1, 3, 2),\n", - "A=\n", - "[1 2 3],\n", - "b=(3),\n", - "u=(1, 1, 1).\n" - ] - } - ], - "prompt_number": 10 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "C\u00e1lculo de un test set, simple y llan0; lo hago para los dos problemas que hemos definido, `P1` y `P`. `test_set` usa el algoritmo 2.1" - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "set_verbose(0)\n", - "P.test_set()" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 11, - "text": [ - "[(1, 0, 0),\n", - " (0, 1, 0),\n", - " (0, 0, 1),\n", - " (-1, 1, 0),\n", - " (-1, 0, 1),\n", - " (0, 1, -1),\n", - " (1, 1, -1)]" - ] - } - ], - "prompt_number": 11 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "T=P1.test_set()\n", - "T" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 12, - "text": [ - "[(1, 0, 0),\n", - " (0, 0, 1),\n", - " (1, 0, -1),\n", - " (-1, 0, 2),\n", - " (2, 0, -2),\n", - " (-2, 0, 3),\n", - " (-3, 0, 4),\n", - " (3, 0, -3),\n", - " (-4, 0, 5),\n", - " (4, 0, -4),\n", - " (-5, 0, 6),\n", - " (-6, 0, 7),\n", - " (5, 0, -5),\n", - " (-7, 0, 8),\n", - " (-8, 0, 9),\n", - " (6, 0, -6),\n", - " (-9, 0, 10),\n", - " (7, 0, -7),\n", - " (8, 0, -8),\n", - " (9, 0, -9),\n", - " (10, 0, -10)]" - ] - } - ], - "prompt_number": 12 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Como se puede ver, sobran vectores. Hay que definir varias funciones de reducci\u00f3n, pero no me entretengo. S\u00f3lo un ejemplo:" - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P1.can_reduce_by(v=vector(ZZ, [1,0,-1]), w=vector([10,0,-10]))" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 13, - "text": [ - "True" - ] - } - ], - "prompt_number": 13 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Ahora podemos reducir el test set que hemos encontrado." - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "set_verbose(0) # or 2\n", - "P1.inter_reduce(T)" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 14, - "text": [ - "[(1, 0, 0), (0, 0, 1)]" - ] - } - ], - "prompt_number": 14 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 14 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 14 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "La funci\u00f3n `minimal_test_set` es una implementaci\u00f3n del algoritmo 2.7. Devuelve el \u00fanico test-set minimal." - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P.minimal_test_set()" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 15, - "text": [ - "[(1, 0, 0), (0, 1, 0), (0, 0, 1), (-1, 1, 0), (-1, 0, 1), (0, 1, -1)]" - ] - } - ], - "prompt_number": 15 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P1.minimal_test_set()" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 16, - "text": [ - "[(1, 0, 0), (0, 0, 1)]" - ] - } - ], - "prompt_number": 16 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Por \u00faltimo, `reduced_test_set` primero calcula un test set, sin reducir, y luego aplica reducci\u00f3n. Debe dar lo mismo que antes; el \u00fanico test-set minimal." - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "set_verbose(0)\n", - "P.reduced_test_set()" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 17, - "text": [ - "[(1, 0, 0), (0, 1, 0), (0, 0, 1), (-1, 1, 0), (-1, 0, 1), (0, 1, -1)]" - ] - } - ], - "prompt_number": 17 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P1.reduced_test_set()" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 18, - "text": [ - "[(1, 0, 0), (0, 0, 1)]" - ] - } - ], - "prompt_number": 18 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Algoritmo 2.10. Mejor que el Grobner simple, pero no necesariamente reducido." - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "set_verbose(0)\n", - "P1.test_set_non_reducibles()" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 19, - "text": [ - "[(1, 0, 0),\n", - " (0, 0, 1),\n", - " (1, 0, -1),\n", - " (-1, 0, 2),\n", - " (2, 0, -2),\n", - " (-2, 0, 3),\n", - " (3, 0, -3),\n", - " (-3, 0, 4),\n", - " (4, 0, -4),\n", - " (-4, 0, 5),\n", - " (5, 0, -5),\n", - " (-5, 0, 6),\n", - " (6, 0, -6),\n", - " (-6, 0, 7),\n", - " (7, 0, -7),\n", - " (-7, 0, 8),\n", - " (8, 0, -8),\n", - " (-8, 0, 9),\n", - " (9, 0, -9),\n", - " (-9, 0, 10),\n", - " (10, 0, -10)]" - ] - } - ], - "prompt_number": 19 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P.test_set_non_reducibles()" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 20, - "text": [ - "[(1, 0, 0), (0, 1, 0), (0, 0, 1), (-1, 1, 0), (-1, 0, 1), (0, 1, -1)]" - ] - } - ], - "prompt_number": 20 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Anulo los c\u00f3mputos de tiempo hasta que se estabilicen los algoritmos." - ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "#%timeit -n100 P.groebner_simple()" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 21 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "#%timeit -n100 P.groebner_with_reduction()" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 22 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "#%timeit -n100 P.groebner_reduced()" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 23 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "#%timeit -n100 P.groebner_non_reducibles()" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 24 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "C\u00e1lculo del \u00f3ptimo a partir de un factible. Siempre se usa el test set minimal.\n" + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Algoritmos de U-W-Z para G-S-U.\n", + "\n", + "Esta hoja contiene los resultados de la implementación de los algoritmos del artículo (Urbaniak, Weismantel, Ziegler 1997). De momento no podéis ejecutarla, pero sí podéis leerla, para ver el tipo de cosas que ya podemos calcular. \n", + "\n", + "El núcleo de lo que he hecho es un programa en SAGE de unas 200 líneas. Cuando nos veamos después de pascua, os diré cómo tenéis que llamarlo desde SAGE.\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#%load_ext sage" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%run ipproblem" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Un problema de IP se declara de la siguiente forma:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "P1=IPProblem(c=[1,2,1],\n", + " A=[[1,0,0],[0,0,1]], \n", + " b=[10,10] )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "O bien, simplemente como" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "P1=IPProblem([1,2,1],\n", + " [[1,0,0],[0,0,1]], \n", + " [10,10] )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Para comprobar los datos que hemos metido, podemos \"imprimir\" el problema" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Integer programming problem\n", + " max\\{c.x : A x <= b, 0<=x<=u\\}\n", + "with\n", + "c=(1, 2, 1),\n", + "A=\n", + "[1 0 0]\n", + "[0 0 1],\n", + "b=(10, 10),\n", + "u=(10, 0, 10).\n" ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "%run ipproblem\n", - "P=IPProblem(c=[1,3,2], A=[[1,2,3]], b=[3], u=[1,1,1])\n", - "set_verbose(0)\n", - "P.walk_to_best([0,0,0])\n" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 25, - "text": [ - "(1, 1, 0)" - ] - } - ], - "prompt_number": 25 - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Tambi\u00e9n podemos pedir el camino completo desde el inicial al final." + } + ], + "source": [ + "print P1" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Algunas funciones útiles:\n", + "\n", + "* `is_feasible(v)` determina si el vector v es “feasible”.\n", + "* `order(v,w)` determina si los vectores están ordenados para el orden $\\prec_c$ del artículo\n", + "* `pm_split(v)` devuelve un par de vectores, correspondientes a $v^+$ y $v^-$ del artículo.\n", + "* `succ(v)` devuelve el vector $v^\\prec$ del artículo\n", + "\n", + "Hay más funciones, pero no voy a dar ejemplos de todas. " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P1.is_feasible(vector(ZZ,[0,0,0]))" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "-1" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P1.order(vector(ZZ,[0,2,0]),vector(ZZ,[13,0,1]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Para las comprobaciones, uso el IP del Ejemplo 2.3 del artículo" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Integer programming problem\n", + " max\\{c.x : A x <= b, 0<=x<=u\\}\n", + "with\n", + "c=(1, 3, 2),\n", + "A=\n", + "[1 2 3],\n", + "b=(3),\n", + "u=(1, 1, 1).\n" ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P.walk_to_best([0,0,0], path=True)" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 26, - "text": [ - "[(0, 0, 0), (0, 1, 0), (1, 1, 0)]" - ] - } - ], - "prompt_number": 26 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P.walk_to_best([0,0,1])" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 27, - "text": [ - "(1, 1, 0)" - ] - } - ], - "prompt_number": 27 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P.walk_to_best([10,10,10]) #\u00a0not feasible" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 28, - "text": [ - "False" - ] - } - ], - "prompt_number": 28 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "P1.walk_to_best([0,0,0], path=True)" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 29, - "text": [ - "[(0, 0, 0),\n", - " (1, 0, 0),\n", - " (2, 0, 0),\n", - " (3, 0, 0),\n", - " (4, 0, 0),\n", - " (5, 0, 0),\n", - " (6, 0, 0),\n", - " (7, 0, 0),\n", - " (8, 0, 0),\n", - " (9, 0, 0),\n", - " (10, 0, 0),\n", - " (10, 0, 1),\n", - " (10, 0, 2),\n", - " (10, 0, 3),\n", - " (10, 0, 4),\n", - " (10, 0, 5),\n", - " (10, 0, 6),\n", - " (10, 0, 7),\n", - " (10, 0, 8),\n", - " (10, 0, 9),\n", - " (10, 0, 10)]" - ] - } - ], - "prompt_number": 29 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 29 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 29 } ], - "metadata": {} + "source": [ + "P=IPProblem(c=[1,3,2], A=[[1,2,3]], b=[3], u=[1,1,1])\n", + "print P" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Cálculo de un test set, simple y llan0; lo hago para los dos problemas que hemos definido, `P1` y `P`. `test_set` usa el algoritmo 2.1" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0),\n", + " (0, 1, 0),\n", + " (0, 0, 1),\n", + " (-1, 1, 0),\n", + " (-1, 0, 1),\n", + " (0, 1, -1),\n", + " (1, 1, -1)]" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "set_verbose(0)\n", + "P.test_set()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0),\n", + " (0, 0, 1),\n", + " (1, 0, -1),\n", + " (-1, 0, 2),\n", + " (2, 0, -2),\n", + " (-2, 0, 3),\n", + " (-3, 0, 4),\n", + " (3, 0, -3),\n", + " (-4, 0, 5),\n", + " (4, 0, -4),\n", + " (-5, 0, 6),\n", + " (-6, 0, 7),\n", + " (5, 0, -5),\n", + " (-7, 0, 8),\n", + " (-8, 0, 9),\n", + " (6, 0, -6),\n", + " (-9, 0, 10),\n", + " (7, 0, -7),\n", + " (8, 0, -8),\n", + " (9, 0, -9),\n", + " (10, 0, -10)]" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "T=P1.test_set()\n", + "T" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Como se puede ver, sobran vectores. Hay que definir varias funciones de reducción, pero no me entretengo. Sólo un ejemplo:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P1.can_reduce_by(v=vector(ZZ, [1,0,-1]), w=vector([10,0,-10]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Ahora podemos reducir el test set que hemos encontrado." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0), (0, 0, 1)]" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "set_verbose(0) # or 2\n", + "P1.inter_reduce(T)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "La función `minimal_test_set` es una implementación del algoritmo 2.7. Devuelve el único test-set minimal." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0), (0, 1, 0), (0, 0, 1), (-1, 1, 0), (-1, 0, 1), (0, 1, -1)]" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P.minimal_test_set()" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0), (0, 0, 1)]" + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P1.minimal_test_set()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Por último, `reduced_test_set` primero calcula un test set, sin reducir, y luego aplica reducción. Debe dar lo mismo que antes; el único test-set minimal." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0), (0, 1, 0), (0, 0, 1), (-1, 1, 0), (-1, 0, 1), (0, 1, -1)]" + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "set_verbose(0)\n", + "P.reduced_test_set()" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0), (0, 0, 1)]" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P1.reduced_test_set()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Algoritmo 2.10. Mejor que el Grobner simple, pero no necesariamente reducido." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0),\n", + " (0, 0, 1),\n", + " (1, 0, -1),\n", + " (-1, 0, 2),\n", + " (2, 0, -2),\n", + " (-2, 0, 3),\n", + " (3, 0, -3),\n", + " (-3, 0, 4),\n", + " (4, 0, -4),\n", + " (-4, 0, 5),\n", + " (5, 0, -5),\n", + " (-5, 0, 6),\n", + " (6, 0, -6),\n", + " (-6, 0, 7),\n", + " (7, 0, -7),\n", + " (-7, 0, 8),\n", + " (8, 0, -8),\n", + " (-8, 0, 9),\n", + " (9, 0, -9),\n", + " (-9, 0, 10),\n", + " (10, 0, -10)]" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "set_verbose(0)\n", + "P1.test_set_non_reducibles()" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0), (0, 1, 0), (0, 0, 1), (-1, 1, 0), (-1, 0, 1), (0, 1, -1)]" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P.test_set_non_reducibles()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Anulo los cómputos de tiempo hasta que se estabilicen los algoritmos." + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#%timeit -n100 P.groebner_simple()" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#%timeit -n100 P.groebner_with_reduction()" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#%timeit -n100 P.groebner_reduced()" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#%timeit -n100 P.groebner_non_reducibles()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Cálculo del óptimo a partir de un factible. Siempre se usa el test set minimal.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(1, 1, 0)" + ] + }, + "execution_count": 25, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%run ipproblem\n", + "P=IPProblem(c=[1,3,2], A=[[1,2,3]], b=[3], u=[1,1,1])\n", + "set_verbose(0)\n", + "P.walk_to_best([0,0,0])\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "También podemos pedir el camino completo desde el inicial al final." + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(0, 0, 0), (0, 1, 0), (1, 1, 0)]" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P.walk_to_best([0,0,0], path=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(1, 1, 0)" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P.walk_to_best([0,0,1])" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "False" + ] + }, + "execution_count": 28, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P.walk_to_best([10,10,10]) # not feasible" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(0, 0, 0),\n", + " (1, 0, 0),\n", + " (2, 0, 0),\n", + " (3, 0, 0),\n", + " (4, 0, 0),\n", + " (5, 0, 0),\n", + " (6, 0, 0),\n", + " (7, 0, 0),\n", + " (8, 0, 0),\n", + " (9, 0, 0),\n", + " (10, 0, 0),\n", + " (10, 0, 1),\n", + " (10, 0, 2),\n", + " (10, 0, 3),\n", + " (10, 0, 4),\n", + " (10, 0, 5),\n", + " (10, 0, 6),\n", + " (10, 0, 7),\n", + " (10, 0, 8),\n", + " (10, 0, 9),\n", + " (10, 0, 10)]" + ] + }, + "execution_count": 29, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "P1.walk_to_best([0,0,0], path=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] } - ] -} \ No newline at end of file + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.13" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/problema_1.ipynb b/problema_1.ipynb index 5b8d688..66adc0e 100644 --- a/problema_1.ipynb +++ b/problema_1.ipynb @@ -1,191 +1,253 @@ { - "metadata": { - "name": "", - "signature": "sha256:b1be195f6c48a319e1c3362ed98a61fb12052a56bb2f8ac10f76d0ecfa77fc9a" - }, - "nbformat": 3, - "nbformat_minor": 0, - "worksheets": [ + "cells": [ + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%run ipproblem" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "c=(4, 1, 1, 3, 4, 2, 2, 2, 2, 3, 1, 3)\n", + "A=[[1, 0, 0, 3, 1, 3, 3, 2, 1, 1, 3, 3],\n", + " [0, 3, 2, 0, 3, 2, 2, 2, 3, 0, 1, 2],\n", + " [0, 1, 2, 2, 1, 2, 0, 1, 2, 0, 3, 3],\n", + " [3, 0, 1, 2, 0, 0, 1, 1, 0, 1, 3, 1]]\n", + "b=(17, 2, 13, 9)\n", + "u=(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, { - "cells": [ + "cell_type": "code", + "execution_count": 7, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "set_verbose(0)\n", + "R=IPProblem(c,A,b,u)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "%run ipproblem" - ], - "language": "python", + "data": { + "text/plain": [ + "(168, 66)" + ] + }, + "execution_count": 8, "metadata": {}, - "outputs": [], - "prompt_number": 1 - }, + "output_type": "execute_result" + } + ], + "source": [ + "H1=R.test_set_non_reducibles()\n", + "H=R.inter_reduce(H1)\n", + "len(H1), len(H)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "c=(4, 1, 1, 3, 4, 2, 2, 2, 2, 3, 1, 3)\n", - "A=[[1, 0, 0, 3, 1, 3, 3, 2, 1, 1, 3, 3],\n", - " [0, 3, 2, 0, 3, 2, 2, 2, 3, 0, 1, 2],\n", - " [0, 1, 2, 2, 1, 2, 0, 1, 2, 0, 3, 3],\n", - " [3, 0, 1, 2, 0, 0, 1, 1, 0, 1, 3, 1]]\n", - "b=(17, 2, 13, 9)\n", - "u=(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)\n", - "\n" - ], - "language": "python", + "data": { + "text/plain": [ + "66" + ] + }, + "execution_count": 9, "metadata": {}, - "outputs": [], - "prompt_number": 2 - }, + "output_type": "execute_result" + } + ], + "source": [ + "G=R.minimal_test_set()\n", + "len(G)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Comprobaciones" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [], - "language": "python", + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 10, "metadata": {}, - "outputs": [], - "prompt_number": 2 - }, + "output_type": "execute_result" + } + ], + "source": [ + "(all(R.is_improvement(v) for v in G) and \n", + " all(R.is_improvement(v) for v in H) and \n", + " len(H)==len(G) and \n", + " all(v in G for v in H) \n", + " and all(v in H1 for v in G))" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "set_verbose(0)\n", - "R=IPProblem(c,A,b,u)" - ], - "language": "python", + "data": { + "text/plain": [ + "(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)" + ] + }, + "execution_count": 11, "metadata": {}, - "outputs": [], - "prompt_number": 3 - }, + "output_type": "execute_result" + } + ], + "source": [ + "v=R.get_feasible()\n", + "v" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "H1=R.test_set_non_reducibles()\n", - "H=R.inter_reduce(H1)\n", - "len(H1), len(H)" - ], - "language": "python", + "data": { + "text/plain": [ + "[(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0),\n", + " (1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1)]" + ] + }, + "execution_count": 12, "metadata": {}, - "outputs": [] - }, + "output_type": "execute_result" + } + ], + "source": [ + "R.walk_to_best(v, path=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "G=R.minimal_test_set()\n", - "len(G)" - ], - "language": "python", + "data": { + "text/plain": [ + "[4, 7, 10, 13]" + ] + }, + "execution_count": 13, "metadata": {}, - "outputs": [] - }, + "output_type": "execute_result" + } + ], + "source": [ + "[R.cost(v) for v in _]" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Comprobaciones" + "name": "stdout", + "output_type": "stream", + "text": [ + "[(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), (0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0), (0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0), (0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0), (0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0), (0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0), (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0), (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0), (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1), (1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0), (1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0), (1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0), (1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0), (1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0), (1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0), (1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0), (1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1), (0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 0, 0), (0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0), (0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 0, 0), (0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0), (0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0), (0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0), (0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1), (0, 0, 0, 1, 0, -1, 0, 0, 0, 0, 0, 0), (0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0), (0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 0), (0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0), (0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0), (0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1), (0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0), (0, 0, 0, 0, 0, 1, 0, -1, 0, 0, 0, 0), (0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0), (0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0), (0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1), (0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0), (0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0), (0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0), (0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1), (0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0), (0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0), (0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1), (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0), (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, -1), (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1), (1, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0), (1, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0), (-1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0), (-1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0), (-1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0), (-1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1), (-1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0), (-1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0), (-1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1), (0, 0, 1, -1, 0, 0, 0, 0, 0, 1, 0, 0), (0, 0, 0, -1, 0, 1, 0, 0, 0, 1, 0, 0), (1, 0, 0, -1, 0, 1, -1, 0, 0, 0, 0, 0), (1, 0, 0, -1, 0, 1, 0, -1, 0, 0, 0, 0), (1, 0, 0, -1, 0, 1, 0, 0, 0, -1, 0, 0), (1, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, -1), (0, 0, 0, -1, 0, 0, 1, 0, 0, 1, 0, 0), (0, 0, 0, -1, 0, 0, 0, 1, 0, 1, 0, 0), (0, 0, -1, -1, 0, 0, 0, 1, 0, 1, 0, 0), (0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 1), (1, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, -1), (0, 0, 0, -1, 0, -1, 0, 0, 0, 1, 0, 1)]\n" ] - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "(all(R.is_improvement(v) for v in G) and \n", - " all(R.is_improvement(v) for v in H) and \n", - " len(H)==len(G) and \n", - " all(v in G for v in H) \n", - " and all(v in H1 for v in G))" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 6, - "text": [ - "True" - ] - } - ], - "prompt_number": 6 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "v=R.get_feasible()\n", - "v" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 7, - "text": [ - "(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)" - ] - } - ], - "prompt_number": 7 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "R.walk_to_best(v, path=True)" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 8, - "text": [ - "[(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),\n", - " (1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),\n", - " (1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0),\n", - " (1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1)]" - ] - } - ], - "prompt_number": 8 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "[R.cost(v) for v in _]" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 9, - "text": [ - "[4, 7, 10, 13]" - ] - } - ], - "prompt_number": 9 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [], - "language": "python", - "metadata": {}, - "outputs": [] } ], - "metadata": {} + "source": [ + "print G" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] } - ] -} \ No newline at end of file + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 7.5.1", + "language": "", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.13" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/problema_1b.ipynb b/problema_1b.ipynb index 2029ebf..c785235 100644 --- a/problema_1b.ipynb +++ b/problema_1b.ipynb @@ -1,388 +1,415 @@ { - "metadata": { - "name": "", - "signature": "sha256:0aaa10cffc9ae9f9651a957482fa8267e78dc329502acb33a9aafe48571a03b7" - }, - "nbformat": 3, - "nbformat_minor": 0, - "worksheets": [ + "cells": [ { - "cells": [ - { - "cell_type": "code", - "collapsed": false, - "input": [ - "%run ipproblem" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 1 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "c=(4, 1, 1, 3, 4, 2, 2, 2, 2, 3, 1, 3)\n", - "A=[[1, 0, 0, 3, 1, 3, 3, 2, 1, 1, 3, 3],\n", - " [0, 3, 2, 0, 3, 2, 2, 2, 3, 0, 1, 2],\n", - " [0, 1, 2, 2, 1, 2, 0, 1, 2, 0, 3, 3],\n", - " [3, 0, 1, 2, 0, 0, 1, 1, 0, 1, 3, 1]]\n", - "b=(17, 2, 13, 9)\n", - "u=(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)\n", - "\n" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 2 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 2 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "set_verbose(0)\n", - "R=IPProblem(c,A,b,u)" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 3 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "G1=R.test_set_non_reducibles()" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 25 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "len(G1)" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 26, - "text": [ - "82" - ] - } - ], - "prompt_number": 26 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "def check_set(R,G1):\n", - " for v in G1:\n", - " vp,vm=R.pm_split(v)\n", - " print (v,R.order(v, vector(ZZ,len(v)*[0]))>0, R.is_feasible(vp), R.is_feasible(vm))" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 27 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "G=R.minimal_test_set()" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 28 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "len(G)" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 29, - "text": [ - "66" - ] - } - ], - "prompt_number": 29 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "#check_set(R,G)" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 9 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "H=R.inter_reduce(G1)\n", - "len(H)" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 10, - "text": [ - "61" - ] - } - ], - "prompt_number": 10 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "Z=R.inter_reduce(G)\n", - "len(Z)" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 11, - "text": [ - "66" - ] - } - ], - "prompt_number": 11 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "Z==G" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 12, - "text": [ - "True" - ] - } - ], - "prompt_number": 12 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "for v in H:\n", - " if not v in G:\n", - " print v\n", - " " - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 13 - }, + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%run ipproblem" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "c=(4, 1, 1, 3, 4, 2, 2, 2, 2, 3, 1, 3)\n", + "A=[[1, 0, 0, 3, 1, 3, 3, 2, 1, 1, 3, 3],\n", + " [0, 3, 2, 0, 3, 2, 2, 2, 3, 0, 1, 2],\n", + " [0, 1, 2, 2, 1, 2, 0, 1, 2, 0, 3, 3],\n", + " [3, 0, 1, 2, 0, 0, 1, 1, 0, 1, 3, 1]]\n", + "b=(17, 2, 13, 9)\n", + "u=(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "set_verbose(0)\n", + "R=IPProblem(c,A,b,u)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "G1=R.test_set_non_reducibles()" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "l=[]\n", - "for v in G:\n", - " if not v in H:\n", - " l.append(v)\n", - "print l" - ], - "language": "python", + "data": { + "text/plain": [ + "82" + ] + }, + "execution_count": 26, "metadata": {}, - "outputs": [ - { - "output_type": "stream", - "stream": "stdout", - "text": [ - "[(1, 0, 0, -1, 0, 1, 0, 0, 0, -1, 0, 0), (1, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, -1), (0, 0, -1, -1, 0, 0, 0, 1, 0, 1, 0, 0), (1, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, -1), (0, 0, 0, -1, 0, -1, 0, 0, 0, 1, 0, 1)]\n" - ] - } - ], - "prompt_number": 14 - }, + "output_type": "execute_result" + } + ], + "source": [ + "len(G1)" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def check_set(R,G1):\n", + " for v in G1:\n", + " vp,vm=R.pm_split(v)\n", + " print (v,R.order(v, vector(ZZ,len(v)*[0]))>0, R.is_feasible(vp), R.is_feasible(vm))" + ] + }, + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "G=R.minimal_test_set()" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "check_set(R,l)" - ], - "language": "python", + "data": { + "text/plain": [ + "66" + ] + }, + "execution_count": 29, "metadata": {}, - "outputs": [ - { - "output_type": "stream", - "stream": "stdout", - "text": [ - "((1, 0, 0, -1, 0, 1, 0, 0, 0, -1, 0, 0), True, True, True)\n", - "((1, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, -1), True, True, True)\n", - "((0, 0, -1, -1, 0, 0, 0, 1, 0, 1, 0, 0), True, True, True)\n", - "((1, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, -1), True, True, True)\n", - "((0, 0, 0, -1, 0, -1, 0, 0, 0, 1, 0, 1), True, True, True)\n" - ] - } - ], - "prompt_number": 15 - }, + "output_type": "execute_result" + } + ], + "source": [ + "len(G)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "#check_set(R,G)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "R.get_feasible()" - ], - "language": "python", + "data": { + "text/plain": [ + "61" + ] + }, + "execution_count": 10, "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 16, - "text": [ - "(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)" - ] - } - ], - "prompt_number": 16 - }, + "output_type": "execute_result" + } + ], + "source": [ + "H=R.inter_reduce(G1)\n", + "len(H)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "v=R.minimal[0]" - ], - "language": "python", + "data": { + "text/plain": [ + "66" + ] + }, + "execution_count": 11, "metadata": {}, - "outputs": [], - "prompt_number": 17 - }, + "output_type": "execute_result" + } + ], + "source": [ + "Z=R.inter_reduce(G)\n", + "len(Z)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "R.walk_to_best(v, path=True)" - ], - "language": "python", + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 12, "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 18, - "text": [ - "[(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),\n", - " (1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),\n", - " (1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0),\n", - " (1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1)]" - ] - } - ], - "prompt_number": 18 - }, + "output_type": "execute_result" + } + ], + "source": [ + "Z==G" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "for v in H:\n", + " if not v in G:\n", + " print v\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "[R.cost(v) for v in _]" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "metadata": {}, - "output_type": "pyout", - "prompt_number": 19, - "text": [ - "[4, 7, 10, 13]" - ] - } - ], - "prompt_number": 19 - }, + "name": "stdout", + "output_type": "stream", + "text": [ + "[(1, 0, 0, -1, 0, 1, 0, 0, 0, -1, 0, 0), (1, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, -1), (0, 0, -1, -1, 0, 0, 0, 1, 0, 1, 0, 0), (1, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, -1), (0, 0, 0, -1, 0, -1, 0, 0, 0, 1, 0, 1)]\n" + ] + } + ], + "source": [ + "l=[]\n", + "for v in G:\n", + " if not v in H:\n", + " l.append(v)\n", + "print l" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "set_verbose(1)\n", - "#R.test_set()" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 20 - }, + "name": "stdout", + "output_type": "stream", + "text": [ + "((1, 0, 0, -1, 0, 1, 0, 0, 0, -1, 0, 0), True, True, True)\n", + "((1, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, -1), True, True, True)\n", + "((0, 0, -1, -1, 0, 0, 0, 1, 0, 1, 0, 0), True, True, True)\n", + "((1, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, -1), True, True, True)\n", + "((0, 0, 0, -1, 0, -1, 0, 0, 0, 1, 0, 1), True, True, True)\n" + ] + } + ], + "source": [ + "check_set(R,l)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [ - "set_verbose(0)" - ], - "language": "python", + "data": { + "text/plain": [ + "(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)" + ] + }, + "execution_count": 16, "metadata": {}, - "outputs": [], - "prompt_number": 21 - }, + "output_type": "execute_result" + } + ], + "source": [ + "R.get_feasible()" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "v=R.minimal[0]" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [], - "language": "python", + "data": { + "text/plain": [ + "[(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0),\n", + " (1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1)]" + ] + }, + "execution_count": 18, "metadata": {}, - "outputs": [], - "prompt_number": 21 - }, + "output_type": "execute_result" + } + ], + "source": [ + "R.walk_to_best(v, path=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "collapsed": false + }, + "outputs": [ { - "cell_type": "code", - "collapsed": false, - "input": [], - "language": "python", + "data": { + "text/plain": [ + "[4, 7, 10, 13]" + ] + }, + "execution_count": 19, "metadata": {}, - "outputs": [] + "output_type": "execute_result" } ], - "metadata": {} + "source": [ + "[R.cost(v) for v in _]" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "set_verbose(1)\n", + "#R.test_set()" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "set_verbose(0)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.10" } - ] -} \ No newline at end of file + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/resolutor.py b/resolutor.py new file mode 100644 index 0000000..fb10630 --- /dev/null +++ b/resolutor.py @@ -0,0 +1,72 @@ +import subprocess +def run_4ti2(command, data, results): + """Run 4ti2 and """ + #cleanup_42() + setup_42(data) + s = subprocess.call( + ["./%s %s"% (command, TMPF)], + shell=True, + # stdout=subprocess.PIPE, + # stderr=subprocess.PIPE, + ) + if s: + raise RuntimeError("ha fallado 4ti2") + + return [load_42_matrix(r) for r in results] + +def load_42_matrix(ext): + """docstring for load_42_matrix""" + with open(TMPF+ext) as f: + s=f.read().splitlines() + return [line_to_list(l) for l in s[1:]] + + +def write_42_file(ext, string): + """docstring for write_42""" + fname = TMPF + ext + with open(fname, 'w') as f: + f.write(string) + +def matrix_to_42(A): + """Print a representation of matrix M as 4ti2 input. M should be an integer matrix.""" + M = matrix(ZZ,A) + n, m = M.dimensions() + return "%s %s\n" % (n,m) + "\n".join([row_to_42(r) for r in M.rows()]) + "\n" + +def row_to_42(r): + """Convert a vector to 4ti2 input format. Should be an integer vector""" + return " ".join([str(i) for i in r]) + +def setup_42(data): + for (k, s) in data.iteritems(): + write_42_file(k,s) + + +def bigsolve(A,b,c,u): + + (d,n)=A.dimensions() + + lower = identity_matrix(n).augment(zero_matrix(n,d)).augment(identity_matrix(n)) + + bigA = A.augment(identity_matrix(d)).augment(zero_matrix(d,n)).stack(lower) + bigb = matrix(b).augment(matrix(c)) + bigc = matrix(c).augment(zero_matrix(1,n+d)) + + + zsol = vector(ZZ,n*[0]+list(b)+n*[1]) + #bigA, bigb, bigc, zsol + data = { + ".mat": matrix_to_42(bigA), + ".zsol": matrix_to_42(matrix(zsol)), + ".cost": matrix_to_42(-matrix(bigc)) + } + GG = matrix(ZZ,run_4ti2("../42/pruebas/groebner -t lp", data, [".gro"])[0])[:,0:n] + return GG.rows() + +def line_to_list(s): + """docstring for line_to_list""" + return [Integer(i) for i in s.strip().split()] + + +TMPF = "temporary_file__" + diff --git a/resolutor.pyc b/resolutor.pyc new file mode 100644 index 0000000..6c6675b Binary files /dev/null and b/resolutor.pyc differ diff --git a/roer.sobj b/roer.sobj new file mode 100644 index 0000000..c47f0a5 Binary files /dev/null and b/roer.sobj differ diff --git a/temporary_file__.cost b/temporary_file__.cost new file mode 100644 index 0000000..740fd72 --- /dev/null +++ b/temporary_file__.cost @@ -0,0 +1,2 @@ +1 28 +-4 -1 -1 -3 -4 -2 -2 -2 -2 -3 -1 -3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diff --git a/temporary_file__.gro b/temporary_file__.gro new file mode 100644 index 0000000..f24f7d6 --- /dev/null +++ b/temporary_file__.gro @@ -0,0 +1,67 @@ +66 28 +-1 0 0 0 0 -1 0 0 0 1 0 1 0 0 -1 1 1 0 0 0 0 1 0 0 0 -1 0 -1 +-1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 +-1 0 0 0 0 0 0 0 0 0 0 1 -2 -2 -3 2 1 0 0 0 0 0 0 0 0 0 0 -1 +-1 0 0 0 0 0 0 0 0 0 1 0 -2 -1 -3 0 1 0 0 0 0 0 0 0 0 0 -1 0 +-1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 2 1 0 0 0 0 0 0 0 0 -1 0 0 +-1 0 0 0 0 0 0 1 0 0 0 0 -1 -2 -1 2 1 0 0 0 0 0 0 -1 0 0 0 0 +-1 0 0 0 0 0 1 0 0 0 0 0 -2 -2 0 2 1 0 0 0 0 0 -1 0 0 0 0 0 +-1 0 0 0 0 1 0 0 0 0 0 0 -2 -2 -2 3 1 0 0 0 0 -1 0 0 0 0 0 0 +-1 0 0 1 0 -1 0 0 0 0 0 1 -2 0 -3 0 1 0 0 -1 0 1 0 0 0 0 0 -1 +-1 0 0 1 0 -1 0 0 0 1 0 0 0 2 0 0 1 0 0 -1 0 1 0 0 0 -1 0 0 +-1 0 0 1 0 -1 0 1 0 0 0 0 -1 0 -1 0 1 0 0 -1 0 1 0 -1 0 0 0 0 +-1 0 0 1 0 -1 1 0 0 0 0 0 -2 0 0 0 1 0 0 -1 0 1 -1 0 0 0 0 0 +-1 0 0 1 0 0 0 0 0 0 0 0 -2 0 -2 1 1 0 0 -1 0 0 0 0 0 0 0 0 +-1 0 1 0 0 0 0 0 0 0 0 0 1 -2 -2 2 1 0 -1 0 0 0 0 0 0 0 0 0 +-1 0 1 0 0 0 0 0 0 1 0 0 0 -2 -2 1 1 0 -1 0 0 0 0 0 0 -1 0 0 +-1 0 1 1 0 0 0 0 0 0 0 0 -2 -2 -4 0 1 0 -1 -1 0 0 0 0 0 0 0 0 + 0 0 -1 0 0 0 0 0 0 0 0 0 0 2 2 1 0 0 1 0 0 0 0 0 0 0 0 0 + 0 0 -1 0 0 0 0 0 0 0 1 0 -3 1 -1 -2 0 0 1 0 0 0 0 0 0 0 -1 0 + 0 0 -1 1 0 0 0 0 0 -1 0 0 -2 2 0 0 0 0 1 -1 0 0 0 0 0 1 0 0 + 0 0 0 -1 0 0 0 0 0 0 0 0 3 0 2 2 0 0 0 1 0 0 0 0 0 0 0 0 + 0 0 0 -1 0 0 0 0 0 0 0 1 0 -2 -1 1 0 0 0 1 0 0 0 0 0 0 0 -1 + 0 0 0 -1 0 0 0 0 0 0 1 0 0 -1 -1 -1 0 0 0 1 0 0 0 0 0 0 -1 0 + 0 0 0 -1 0 0 0 0 0 1 0 0 2 0 2 1 0 0 0 1 0 0 0 0 0 -1 0 0 + 0 0 0 -1 0 0 0 1 0 0 0 0 1 -2 1 1 0 0 0 1 0 0 0 -1 0 0 0 0 + 0 0 0 -1 0 0 1 0 0 0 0 0 0 -2 2 1 0 0 0 1 0 0 -1 0 0 0 0 0 + 0 0 0 -1 0 1 0 0 0 0 0 0 0 -2 0 2 0 0 0 1 0 -1 0 0 0 0 0 0 + 0 0 0 0 0 -1 0 0 0 0 0 0 3 2 2 0 0 0 0 0 0 1 0 0 0 0 0 0 + 0 0 0 0 0 -1 0 0 0 0 1 0 0 1 -1 -3 0 0 0 0 0 1 0 0 0 0 -1 0 + 0 0 0 0 0 -1 0 1 0 0 0 0 1 0 1 -1 0 0 0 0 0 1 0 -1 0 0 0 0 + 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 2 -1 0 0 0 0 0 1 -1 0 0 0 0 0 + 0 0 0 0 0 0 -1 0 0 0 0 0 3 2 0 1 0 0 0 0 0 0 1 0 0 0 0 0 + 0 0 0 0 0 0 -1 0 0 0 1 0 0 1 -3 -2 0 0 0 0 0 0 1 0 0 0 -1 0 + 0 0 0 0 0 0 -1 1 0 0 0 0 1 0 -1 0 0 0 0 0 0 0 1 -1 0 0 0 0 + 0 0 0 0 0 0 0 -1 0 0 0 0 2 2 1 1 0 0 0 0 0 0 0 1 0 0 0 0 + 0 0 0 0 0 0 0 -1 0 0 1 0 -1 1 -2 -2 0 0 0 0 0 0 0 1 0 0 -1 0 + 0 0 0 0 0 0 0 0 0 -1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 + 0 0 0 0 0 0 0 0 0 -1 0 1 -2 -2 -3 0 0 0 0 0 0 0 0 0 0 1 0 -1 + 0 0 0 0 0 0 0 0 0 -1 1 0 -2 -1 -3 -2 0 0 0 0 0 0 0 0 0 1 -1 0 + 0 0 0 0 0 0 0 0 0 0 -1 0 3 1 3 3 0 0 0 0 0 0 0 0 0 0 1 0 + 0 0 0 0 0 0 0 0 0 0 0 -1 3 2 3 1 0 0 0 0 0 0 0 0 0 0 0 1 + 0 0 0 0 0 0 0 0 0 0 1 -1 0 1 0 -2 0 0 0 0 0 0 0 0 0 0 -1 1 + 0 0 0 0 0 0 0 1 0 -1 0 0 -1 -2 -1 0 0 0 0 0 0 0 0 -1 0 1 0 0 + 0 0 0 0 0 0 0 1 0 0 0 -1 1 0 2 0 0 0 0 0 0 0 0 -1 0 0 0 1 + 0 0 0 0 0 0 1 0 0 -1 0 0 -2 -2 0 0 0 0 0 0 0 0 -1 0 0 1 0 0 + 0 0 0 0 0 0 1 0 0 0 0 -1 0 0 3 0 0 0 0 0 0 0 -1 0 0 0 0 1 + 0 0 0 0 0 1 0 0 0 -1 0 0 -2 -2 -2 1 0 0 0 0 0 -1 0 0 0 1 0 0 + 0 0 0 0 0 1 0 0 0 0 0 -1 0 0 1 1 0 0 0 0 0 -1 0 0 0 0 0 1 + 0 0 0 1 0 -1 0 0 0 -1 0 0 1 2 0 -1 0 0 0 -1 0 1 0 0 0 1 0 0 + 0 0 0 1 0 0 -1 0 0 -1 0 0 1 2 -2 0 0 0 0 -1 0 0 1 0 0 1 0 0 + 0 0 0 1 0 0 0 -1 0 -1 0 0 0 2 -1 0 0 0 0 -1 0 0 0 1 0 1 0 0 + 0 0 0 1 0 0 0 0 0 -1 0 -1 1 2 1 0 0 0 0 -1 0 0 0 0 0 1 0 1 + 0 0 0 1 0 1 0 0 0 -1 0 -1 -2 0 -1 0 0 0 0 -1 0 -1 0 0 0 1 0 1 + 0 0 1 -1 0 0 0 0 0 0 0 0 3 -2 0 1 0 0 -1 1 0 0 0 0 0 0 0 0 + 0 0 1 0 0 -1 0 0 0 0 0 0 3 0 0 -1 0 0 -1 0 0 1 0 0 0 0 0 0 + 0 0 1 0 0 0 -1 0 0 0 0 0 3 0 -2 0 0 0 -1 0 0 0 1 0 0 0 0 0 + 0 0 1 0 0 0 0 -1 0 0 0 0 2 0 -1 0 0 0 -1 0 0 0 0 1 0 0 0 0 + 0 0 1 0 0 0 0 0 0 -1 0 0 1 -2 -2 0 0 0 -1 0 0 0 0 0 0 1 0 0 + 0 0 1 0 0 0 0 0 0 0 0 -1 3 0 1 0 0 0 -1 0 0 0 0 0 0 0 0 1 + 0 0 1 1 0 0 0 -1 0 -1 0 0 0 0 -3 -1 0 0 -1 -1 0 0 0 1 0 1 0 0 + 1 0 0 -1 0 0 -1 0 0 0 0 0 5 2 2 0 -1 0 0 1 0 0 1 0 0 0 0 0 + 1 0 0 -1 0 0 0 -1 0 0 0 0 4 2 3 0 -1 0 0 1 0 0 0 1 0 0 0 0 + 1 0 0 -1 0 0 0 0 0 -1 0 0 3 0 2 0 -1 0 0 1 0 0 0 0 0 1 0 0 + 1 0 0 -1 0 0 0 0 0 0 0 -1 5 2 5 0 -1 0 0 1 0 0 0 0 0 0 0 1 + 1 0 0 0 0 0 -1 0 0 -1 0 0 3 2 0 -1 -1 0 0 0 0 0 1 0 0 1 0 0 + 1 0 0 0 0 0 0 -1 0 -1 0 0 2 2 1 -1 -1 0 0 0 0 0 0 1 0 1 0 0 + 1 0 0 0 0 0 0 0 0 -1 0 -1 3 2 3 -1 -1 0 0 0 0 0 0 0 0 1 0 1 diff --git a/temporary_file__.mat b/temporary_file__.mat new file mode 100644 index 0000000..2bc07d2 --- /dev/null +++ b/temporary_file__.mat @@ -0,0 +1,17 @@ +16 28 +1 0 0 3 1 3 3 2 1 1 3 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 3 2 0 3 2 2 2 3 0 1 2 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 1 2 2 1 2 0 1 2 0 3 3 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 +3 0 1 2 0 0 1 1 0 1 3 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 +1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 +0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 +0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 +0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 +0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 +0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 +0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 +0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 +0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 +0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 +0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 +0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 diff --git a/temporary_file__.zsol b/temporary_file__.zsol new file mode 100644 index 0000000..73bf237 --- /dev/null +++ b/temporary_file__.zsol @@ -0,0 +1,2 @@ +1 28 +0 0 0 0 0 0 0 0 0 0 0 0 17 2 13 9 1 1 1 1 1 1 1 1 1 1 1 1 diff --git a/truncated.ipynb b/truncated.ipynb new file mode 100644 index 0000000..8c5af24 --- /dev/null +++ b/truncated.ipynb @@ -0,0 +1,1294 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 28, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def random_vvv(l,n):\n", + " return [vector(ZZ,[ZZ.random_element() for i in xrange(l)]) for j in xrange(n)]\n", + "\n", + "def random_lll(l,n):\n", + " return [ [ZZ.random_element() for i in xrange(l)] for j in xrange(n)]\n", + "\n", + "longitud = 50\n", + "numero = 100000\n", + "\n", + "#xxx = random_lll(longitud, numero)\n", + "#yyy = random_lll(longitud, numero)\n", + "#zzz = [vector(ZZ,z) for z in xxx]\n", + "#ttt = [vector(ZZ,z) for z in yyy]\n", + "\n", + "#vvv = random_vvv(5000, 1000)\n", + "#yyy = random_vvv(5000, 1000)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 159, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Primer ejemplo" + ] + }, + { + "cell_type": "code", + "execution_count": 264, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "reset()\n", + "load(\"resolutor.py\")\n", + "load(\"truncated.py\")" + ] + }, + { + "cell_type": "code", + "execution_count": 265, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(\n", + "[ 3 1 11 2 3 5 3] \n", + "[ 4 5 0 1 7 4 6] \n", + "[ 5 6 1 9 2 3 3], (31, 27, 38), (23, 15, 6, 7, 1, 53, 4),\n", + "\n", + "(38, 38, 38, 38, 38, 38, 38)\n", + ")" + ] + }, + "execution_count": 265, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A1, b1, c1, u1" + ] + }, + { + "cell_type": "code", + "execution_count": 267, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "4588 criterion 1, 1421 criterion 2, 8437 criterion 3, 64 non feasible, 541 reductions, of which 374 to zero\n", + " " + ] + } + ], + "source": [ + "%%prun\n", + "\n", + "vecs, g = bubu(A1,b1,c1,vector(ZZ, 7*[38]))\n" + ] + }, + { + "cell_type": "code", + "execution_count": 183, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "174" + ] + }, + "execution_count": 183, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(vecs)" + ] + }, + { + "cell_type": "code", + "execution_count": 184, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "GG = bigsolve(A1,b1,c1,u1)" + ] + }, + { + "cell_type": "code", + "execution_count": 185, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "63" + ] + }, + "execution_count": 185, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(GG)" + ] + }, + { + "cell_type": "code", + "execution_count": 186, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "for v in GG:\n", + " if not (v in vecs or -v in vecs): \n", + " print copy(v), g.reduce(copy(v))" + ] + }, + { + "cell_type": "code", + "execution_count": 187, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "ename": "AttributeError", + "evalue": "'PartialBasis' object has no attribute 'self_reduce'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m/Applications/SageMath-7.5.1.app/Contents/Resources/sage/local/lib/python2.7/site-packages/sage/all_cmdline.pyc\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mself_reduce\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[0;31mAttributeError\u001b[0m: 'PartialBasis' object has no attribute 'self_reduce'" + ] + } + ], + "source": [ + "g.self_reduce()" + ] + }, + { + "cell_type": "code", + "execution_count": 188, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "174" + ] + }, + "execution_count": 188, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(g.vectors)" + ] + }, + { + "cell_type": "code", + "execution_count": 189, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "for v in GG:\n", + " if not (v in vecs or -v in vecs): \n", + " print copy(v), g.reduce(copy(v))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Re-test" + ] + }, + { + "cell_type": "code", + "execution_count": 238, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "74 criterion 1, 9 criterion 2, 72 criterion 3, 10 non feasible, 25 reductions, of which 10 to zero\n" + ] + } + ], + "source": [ + "load(\"resolutor.py\")\n", + "load(\"truncated.py\")\n", + "A = matrix(ZZ,[[999, 18, 625, 416, 773]])\n", + "b = vector(ZZ, [1415])\n", + "c = vector(ZZ, (252, 537, 965, 237, 724))\n", + "u = vector(ZZ,5*[1])\n", + "\n", + "vv, g = bubu(A,b,c,u)" + ] + }, + { + "cell_type": "code", + "execution_count": 239, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0, 0, 0),\n", + " (0, 1, 0, 0, 0),\n", + " (0, 0, 1, 0, 0),\n", + " (0, 0, 0, 1, 0),\n", + " (0, 0, 0, 0, 1),\n", + " (0, 0, 0, -1, 1),\n", + " (0, 0, 1, 1, -1),\n", + " (0, 1, -1, -1, 1),\n", + " (0, 1, 0, 1, -1),\n", + " (-1, 0, 0, 0, 1),\n", + " (0, 0, 1, 0, -1),\n", + " (0, -1, 0, 0, 1),\n", + " (0, 0, 1, -1, 0),\n", + " (0, -1, 1, -1, 0),\n", + " (-1, 0, 1, 0, 0),\n", + " (0, 1, 0, -1, 0),\n", + " (1, 0, 0, -1, 0),\n", + " (-1, 1, 0, 1, 0),\n", + " (0, -1, 1, 0, 0),\n", + " (-1, 1, 0, 0, 0)]" + ] + }, + "execution_count": 239, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g.vectors" + ] + }, + { + "cell_type": "code", + "execution_count": 240, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "20" + ] + }, + "execution_count": 240, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(g.vectors)" + ] + }, + { + "cell_type": "code", + "execution_count": 241, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "17" + ] + }, + "execution_count": 241, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "GG = bigsolve(A,b,c,u)\n", + "len(GG)" + ] + }, + { + "cell_type": "code", + "execution_count": 243, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(-1, 0, 0, 0, 0),\n", + " (-1, 0, 0, 1, 0),\n", + " (0, -1, 0, -1, 1),\n", + " (0, -1, 0, 0, 0),\n", + " (0, -1, 0, 1, 0),\n", + " (0, 0, -1, 0, 0),\n", + " (0, 0, -1, 0, 1),\n", + " (0, 0, -1, 1, 0),\n", + " (0, 0, 0, -1, 0),\n", + " (0, 0, 0, 0, -1),\n", + " (0, 0, 0, 1, -1),\n", + " (0, 1, -1, 0, 0),\n", + " (0, 1, -1, 1, 0),\n", + " (0, 1, 0, 0, -1),\n", + " (1, -1, 0, 0, 0),\n", + " (1, 0, -1, 0, 0),\n", + " (1, 0, 0, 0, -1)]" + ] + }, + "execution_count": 243, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + " GG\n" + ] + }, + { + "cell_type": "code", + "execution_count": 244, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "(-1, 1, 0, 1, 0)\n", + "(0, 0, 1, 1, -1)\n", + "(0, 1, -1, -1, 1)\n" + ] + } + ], + "source": [ + "for z in sorted(g.vectors):\n", + " if not (z in GG or -z in GG): print z" + ] + }, + { + "cell_type": "code", + "execution_count": 245, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "for z in GG:\n", + " if not (z in g.vectors or -z in g.vectors): print copy(z), g.reduce(z)" + ] + }, + { + "cell_type": "code", + "execution_count": 198, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "ename": "AttributeError", + "evalue": "'PartialBasis' object has no attribute 'self_reduce'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m/Applications/SageMath-7.5.1.app/Contents/Resources/sage/local/lib/python2.7/site-packages/sage/all_cmdline.pyc\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mself_reduce\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvectors\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mAttributeError\u001b[0m: 'PartialBasis' object has no attribute 'self_reduce'" + ] + } + ], + "source": [ + "g.self_reduce()\n", + "len(g.vectors)" + ] + }, + { + "cell_type": "code", + "execution_count": 199, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "for z in GG:\n", + " if not (z in g.vectors or -z in g.vectors): \n", + " print copy(z), g.reduce(z)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "load(\"resolutor.py\")\n", + "load(\"truncated.py\")" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "265 criterion 1, 87 criterion 2, 536 criterion 3, 126 non feasible, 67 reductions, of which 27 to zero\n" + ] + } + ], + "source": [ + "#%%prun\n", + "\n", + "A1 = matrix(ZZ, [[3,1,11,2,3,5,3], [4, 5, 0, 1, 7,4, 6], [5,6,1,9,2,3,3]])\n", + "c1 = vector(ZZ, [23,15,6,7,1,53,4])\n", + "b1 = vector(ZZ, [10,11,12])\n", + "u1 = vector(ZZ, 7*[1])\n", + "vects, g = bubu(A1,b1,c1,u1)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "47" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(vects)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "ename": "RuntimeError", + "evalue": "ha fallado 4ti2", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mRuntimeError\u001b[0m Traceback (most recent call last)", + "\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mgg\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbigsolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mb1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mc1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mu1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 2\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/Users/soto/Documents/Maths/Papers/optimizacion/gsu/resolutor.py\u001b[0m in \u001b[0;36mbigsolve\u001b[0;34m(A, b, c, u)\u001b[0m\n\u001b[1;32m 61\u001b[0m \u001b[0;34m\".cost\"\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mmatrix_to_42\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mmatrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbigc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 62\u001b[0m }\n\u001b[0;32m---> 63\u001b[0;31m \u001b[0mGG\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmatrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mZZ\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mrun_4ti2\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"../42/pruebas/groebner -t lp\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m\".gro\"\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 64\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mGG\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrows\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 65\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;32m/Users/soto/Documents/Maths/Papers/optimizacion/gsu/resolutor.py\u001b[0m in \u001b[0;36mrun_4ti2\u001b[0;34m(command, data, results)\u001b[0m\n\u001b[1;32m 11\u001b[0m )\n\u001b[1;32m 12\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mRuntimeError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"ha fallado 4ti2\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 14\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 15\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mload_42_matrix\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mr\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mresults\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mRuntimeError\u001b[0m: ha fallado 4ti2" + ] + } + ], + "source": [ + "gg = bigsolve(A1,b1,c1,u1)\n", + "len(gg)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "# Ejemplo duro de roer" + ] + }, + { + "cell_type": "code", + "execution_count": 47, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "\n", + "load(\"truncated.py\")" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1648 criterion 1, 833 criterion 2, 1973 criterion 3, 513 non feasible, 184 reductions, of which 94 to zero\n", + "1648 criterion 1, 833 criterion 2, 1973 criterion 3, 513 non feasible, 184 reductions, of which 94 to zero\n", + "1648 criterion 1, 833 criterion 2, 1973 criterion 3, 513 non feasible, 184 reductions, of which 94 to zero\n", + "1648 criterion 1, 833 criterion 2, 1973 criterion 3, 513 non feasible, 184 reductions, of which 94 to zero\n", + "1648 criterion 1, 833 criterion 2, 1973 criterion 3, 513 non feasible, 184 reductions, of which 94 to zero\n", + "1648 criterion 1, 833 criterion 2, 1973 criterion 3, 513 non feasible, 184 reductions, of which 94 to zero\n", + "1648 criterion 1, 833 criterion 2, 1973 criterion 3, 513 non feasible, 184 reductions, of which 94 to zero\n", + "1648 criterion 1, 833 criterion 2, 1973 criterion 3, 513 non feasible, 184 reductions, of which 94 to zero\n", + "1648 criterion 1, 833 criterion 2, 1973 criterion 3, 513 non feasible, 184 reductions, of which 94 to zero\n", + "3 loops, best of 3: 449 ms per loop\n" + ] + } + ], + "source": [ + "%%timeit -r3 -n3\n", + "c1=vector(ZZ,(4, 1, 1, 3, 4, 2, 2, 2, 2, 3, 1, 3))\n", + "A1=matrix(ZZ, [[1, 0, 0, 3, 1, 3, 3, 2, 1, 1, 3, 3],\n", + " [0, 3, 2, 0, 3, 2, 2, 2, 3, 0, 1, 2],\n", + " [0, 1, 2, 2, 1, 2, 0, 1, 2, 0, 3, 3],\n", + " [3, 0, 1, 2, 0, 0, 1, 1, 0, 1, 3, 1]])\n", + "b1=vector(ZZ,(17, 2, 13, 9))\n", + "u1=vector(ZZ,(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))\n", + "v, g = bubu(A1,b1,c1,u1)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "reset()\n", + "load(\"truncated.py\")\n", + "load(\"resolutor.py\")" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "1648 criterion 1, 833 criterion 2, 1973 criterion 3, 513 non feasible, 184 reductions, of which 94 to zero\n", + " " + ] + } + ], + "source": [ + "%%prun\n", + "c1=vector(ZZ,(4, 1, 1, 3, 4, 2, 2, 2, 2, 3, 1, 3))\n", + "A1=matrix(ZZ, [[1, 0, 0, 3, 1, 3, 3, 2, 1, 1, 3, 3],\n", + " [0, 3, 2, 0, 3, 2, 2, 2, 3, 0, 1, 2],\n", + " [0, 1, 2, 2, 1, 2, 0, 1, 2, 0, 3, 3],\n", + " [3, 0, 1, 2, 0, 0, 1, 1, 0, 1, 3, 1]])\n", + "b1=vector(ZZ,(17, 2, 13, 9))\n", + "u1=vector(ZZ,(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1))\n", + "v, g = bubu(A1,b1,c1,u1)\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 621, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, -1),\n", + " (0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1),\n", + " (0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, -1),\n", + " (0, 0, 0, -1, 0, 0, 0, -1, 0, 1, 0, 1),\n", + " (1, 0, 0, -1, 0, 0, 0, -1, 0, 0, 0, 1),\n", + " (-1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0),\n", + " (1, 0, 0, -1, 0, 1, 0, -1, 0, 0, 0, 0),\n", + " (1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1),\n", + " (0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1),\n", + " (0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, -1),\n", + " (0, 0, 0, -1, 0, 0, -1, 0, 0, 1, 0, 1),\n", + " (1, 0, 0, -1, 0, 0, -1, 0, 0, 0, 0, 1),\n", + " (-1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, -1, 0, 1, -1, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1),\n", + " (0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, -1),\n", + " (0, 0, 0, -1, 0, -1, 0, 0, 0, 1, 0, 1),\n", + " (1, 0, 0, -1, 0, -1, 0, 0, 0, 0, 0, 1),\n", + " (1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, -1),\n", + " (-1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1),\n", + " (-1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0),\n", + " (1, 0, 0, -1, 0, 1, 0, 0, 0, -1, 0, 0),\n", + " (-1, 0, 1, 1, 0, -1, 0, 0, 0, 1, 0, 0),\n", + " (1, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, -1),\n", + " (1, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, -1),\n", + " (0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1),\n", + " (0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 1),\n", + " (0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0),\n", + " (0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0),\n", + " (0, 0, 0, -1, 0, 0, 0, 1, 0, 1, 0, 0),\n", + " (0, 0, 0, 1, 0, 0, 1, -1, 0, -1, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, -1, 1, 0, 1, 0, 0),\n", + " (-1, 0, 0, 1, 0, 0, -1, 1, 0, 1, 0, 0),\n", + " (0, 0, 0, 1, 0, 1, 0, -1, 0, -1, 0, 0),\n", + " (0, 0, -1, -1, 0, 0, 0, 1, 0, 1, 0, 0),\n", + " (0, 0, 0, -1, 0, 0, 1, 0, 0, 1, 0, 0),\n", + " (0, 0, 0, 1, 0, 0, -1, 1, 0, -1, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0),\n", + " (0, 0, 0, 1, 0, 0, -1, 1, 0, 0, 0, 0),\n", + " (1, 0, 0, 0, 0, 0, -1, 1, 0, 0, 0, 0),\n", + " (0, 0, 0, 1, 0, 1, -1, 0, 0, -1, 0, 0),\n", + " (0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, -1, 0, 1, 0, 0, 0, 1, 0, 0),\n", + " (0, 0, 0, 1, 0, -1, 0, 1, 0, -1, 0, 0),\n", + " (0, 0, 0, 1, 0, -1, 1, 0, 0, -1, 0, 0),\n", + " (0, 0, 0, 1, 0, -1, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, 1, -1, 0, 0, 0, 0, 0, 1, 0, 0),\n", + " (0, 0, -1, 1, 0, 0, 0, 1, 0, -1, 0, 0),\n", + " (0, 0, 1, 0, 0, 0, 0, -1, 0, 1, 0, 0),\n", + " (0, 0, -1, 1, 0, 0, 1, 0, 0, -1, 0, 0),\n", + " (0, 0, 1, 0, 0, 0, -1, 0, 0, 1, 0, 0),\n", + " (0, 0, -1, 1, 0, 1, 0, 0, 0, -1, 0, 0),\n", + " (1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0),\n", + " (-1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1),\n", + " (-1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0),\n", + " (1, 0, 0, 0, 0, 1, 0, -1, 0, -1, 0, 0),\n", + " (-1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0),\n", + " (1, 0, 0, 0, 0, 1, -1, 0, 0, -1, 0, 0),\n", + " (1, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0),\n", + " (0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1),\n", + " (1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1),\n", + " (1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0),\n", + " (0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0),\n", + " (0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0),\n", + " (0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0),\n", + " (0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0),\n", + " (1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0),\n", + " (0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0),\n", + " (0, 0, 0, 0, 0, 1, 0, -1, 0, 0, 0, 0),\n", + " (1, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0),\n", + " (0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 0),\n", + " (0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0),\n", + " (1, 0, 1, 0, 0, 0, 0, -1, 0, 0, 0, 0),\n", + " (1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0),\n", + " (0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 0, 0),\n", + " (1, 0, 1, 0, 0, 0, -1, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0),\n", + " (0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0)]" + ] + }, + "execution_count": 621, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g.vectors" + ] + }, + { + "cell_type": "code", + "execution_count": 622, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "102" + ] + }, + "execution_count": 622, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(g.vectors)" + ] + }, + { + "cell_type": "code", + "execution_count": 623, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "g.self_reduce()" + ] + }, + { + "cell_type": "code", + "execution_count": 624, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "69" + ] + }, + "execution_count": 624, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(g.vectors)" + ] + }, + { + "cell_type": "code", + "execution_count": 625, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "66" + ] + }, + "execution_count": 625, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "GG = bigsolve(A1,b1,c1,u1)\n", + "len(GG)" + ] + }, + { + "cell_type": "code", + "execution_count": 626, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 1),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, -1),\n", + " (0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 1),\n", + " (-1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0),\n", + " (1, 0, 0, -1, 0, 1, 0, -1, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 1),\n", + " (-1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, -1, 0, 1, -1, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 1),\n", + " (0, 0, 0, -1, 0, -1, 0, 0, 0, 1, 0, 1),\n", + " (1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0, -1),\n", + " (-1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1),\n", + " (-1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0),\n", + " (1, 0, 0, -1, 0, 1, 0, 0, 0, -1, 0, 0),\n", + " (1, 0, -1, -1, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0, -1),\n", + " (0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1),\n", + " (0, 0, 0, -1, 0, 0, 0, 0, 0, 1, 0, 1),\n", + " (0, 0, -1, 0, 0, 0, 0, 0, 0, 1, 0, 0),\n", + " (0, 0, 0, 1, 0, 0, 0, 0, 0, -1, 0, 0),\n", + " (0, 0, 0, -1, 0, 0, 0, 1, 0, 1, 0, 0),\n", + " (0, 0, -1, -1, 0, 0, 0, 1, 0, 1, 0, 0),\n", + " (0, 0, 0, -1, 0, 0, 1, 0, 0, 1, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0),\n", + " (0, 0, 0, 1, 0, 0, -1, 0, 0, 0, 0, 0),\n", + " (0, 0, 0, -1, 0, 1, 0, 0, 0, 1, 0, 0),\n", + " (0, 0, 0, 1, 0, -1, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, 1, -1, 0, 0, 0, 0, 0, 1, 0, 0),\n", + " (1, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0),\n", + " (0, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0),\n", + " (-1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1),\n", + " (-1, 0, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0),\n", + " (-1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0),\n", + " (1, 0, -1, 0, 0, 0, 0, 0, 0, -1, 0, 0),\n", + " (0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 1),\n", + " (1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1),\n", + " (0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0),\n", + " (0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -1, 0),\n", + " (0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -1, 0),\n", + " (0, 0, 0, 1, 0, 0, 0, 0, 0, 0, -1, 0),\n", + " (0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, 0),\n", + " (1, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0),\n", + " (0, 0, 0, 0, 0, 0, 0, -1, 0, 1, 0, 0),\n", + " (0, 0, 0, 0, 0, -1, 0, 0, 0, 1, 0, 0),\n", + " (0, 0, 0, 0, 0, 1, 0, -1, 0, 0, 0, 0),\n", + " (0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0, 0),\n", + " (0, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0, 0),\n", + " (1, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0),\n", + " (0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 0),\n", + " (0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0),\n", + " (0, 0, -1, 0, 0, 1, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0),\n", + " (0, 0, -1, 1, 0, 0, 0, 0, 0, 0, 0, 0),\n", + " (1, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0)]" + ] + }, + "execution_count": 626, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "g.vectors" + ] + }, + { + "cell_type": "code", + "execution_count": 574, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "for v in GG:\n", + " if not (v or -v in g.vectors):\n", + " vv = copy(v)\n", + " print v, g.reduce(vv)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "## Ejemplo de Ucha" + ] + }, + { + "cell_type": "code", + "execution_count": 550, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "reset()\n", + "load(\"truncated.py\")\n", + "load(\"resolutor.py\")\n", + "\n", + "def mload_42_matrix(ff):\n", + " \"\"\"docstring for load_42_matrix\"\"\"\n", + " with open(ff) as f:\n", + " s=f.read().splitlines()\n", + " return [line_to_list(l) for l in s[1:]]\n", + "\n", + "\n", + "Ap = matrix(matrix(ZZ, mload_42_matrix(\"ejemploucha/ejemplosoto.mat\"))[0][:5])\n", + "cp = -vector(ZZ, mload_42_matrix(\"ejemploucha/ejemplosoto.cost\")[0])[:5]\n", + "zsol = vector(ZZ,mload_42_matrix(\"ejemploucha/ejemplosoto.zsol\")[0])\n", + "\n", + "b = vector(ZZ,[1415])\n", + "u = vector(ZZ, 5*[1])" + ] + }, + { + "cell_type": "code", + "execution_count": 551, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "load(\"resolutor.py\")\n", + "gg = bigsolve(Ap,b,cp,u)" + ] + }, + { + "cell_type": "code", + "execution_count": 552, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "17" + ] + }, + "execution_count": 552, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(gg)" + ] + }, + { + "cell_type": "code", + "execution_count": 553, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "load(\"truncated.py\")" + ] + }, + { + "cell_type": "code", + "execution_count": 554, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "74 criterion 1, 9 criterion 2, 72 criterion 3, 10 non feasible, 25 reductions, of which 10 to zero\n" + ] + } + ], + "source": [ + "\n", + "vecs, g = bubu(Ap,b,cp,u)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 555, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[(1, 0, 0, 0, 0),\n", + " (0, 1, 0, 0, 0),\n", + " (0, 0, 1, 0, 0),\n", + " (0, 0, 0, 1, 0),\n", + " (0, 0, 0, 0, 1),\n", + " (0, 0, 0, -1, 1),\n", + " (0, 0, 1, 1, -1),\n", + " (0, 1, -1, -1, 1),\n", + " (0, 1, 0, 1, -1),\n", + " (-1, 0, 0, 0, 1),\n", + " (0, 0, 1, 0, -1),\n", + " (0, -1, 0, 0, 1),\n", + " (0, 0, 1, -1, 0),\n", + " (0, -1, 1, -1, 0),\n", + " (-1, 0, 1, 0, 0),\n", + " (0, 1, 0, -1, 0),\n", + " (1, 0, 0, -1, 0),\n", + " (-1, 1, 0, 1, 0),\n", + " (0, -1, 1, 0, 0),\n", + " (-1, 1, 0, 0, 0)]" + ] + }, + "execution_count": 555, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "vecs" + ] + }, + { + "cell_type": "code", + "execution_count": 556, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "20" + ] + }, + "execution_count": 556, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(g.vectors)" + ] + }, + { + "cell_type": "code", + "execution_count": 557, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "6 (0, 0, 1, 1, -1)\n", + "7 (0, 1, -1, -1, 1)\n", + "17 (-1, 1, 0, 1, 0)\n" + ] + } + ], + "source": [ + "for i,v in enumerate(vecs): \n", + " if not (v in gg or -v in gg): \n", + " print i,v" + ] + }, + { + "cell_type": "code", + "execution_count": 558, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "for w in gg:\n", + " if not (w in vecs or -w in vecs): print w\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 559, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "g.self_reduce()" + ] + }, + { + "cell_type": "code", + "execution_count": 560, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "for w in gg:\n", + " if not (w in g.vectors or -w in g.vectors): print w" + ] + }, + { + "cell_type": "code", + "execution_count": 561, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "(17, 17)" + ] + }, + "execution_count": 561, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "len(g.vectors), len(gg)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "SageMath 7.5.1", + "language": "", + "name": "sagemath" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.13" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/truncated.py b/truncated.py new file mode 100755 index 0000000..3830136 --- /dev/null +++ b/truncated.py @@ -0,0 +1,561 @@ +#!/usr/bin/env python + +from collections import deque +from itertools import izip, imap +import pdb + + +def matprod(A,b): + return [sum(ai*bi for ai,bi in izip(row, b)) for row in A] + +def dotprod(u,v): + return sum(ui*vi for (ui,vi) in izip(u,v)) + + +def getz(v): + """Check if v is greater or equal than zero, component-wise.""" + # return all(vi>=0 for vi in v) #(Slower by a factor of 2) + for vi in v: + if vi<0: return False + return True + +def isz(v): + """Check if v is zero. Slightly faster than .is_zero()""" + # # # return all(vi>=0 for vi in v) #(Slower by a factor of 2) + for vi in v: + if vi!=0: return False + return True + # return v.is_zero() + + +def let(a,b): + "True if and only if a is less or equal than b. Faster than zipping lists" + for i in xrange(len(a)): + if a[i]>b[i]: + return False + return True + +def let_pp(a,b): + "True if and only if a^+ is less than or equal b^+." + for i in xrange(len(a)): + if a[i]>0 and a[i]>max(b[i], 0): + return False + return True + + +def let_pm(a,b): + "True if and only if a^+ is less or equal than b^-." + for i in xrange(len(a)): + if a[i]>0 and a[i]>max(-b[i], 0): + return False + return True + +def let_mp(a,b): + "True if and only if a^- is less or equal than b^+." + for i in xrange(len(a)): + if a[i]<0 and -a[i]>max(b[i], 0): + return False + return True + +def let_mm(a,b): + "True if and only if a^- is less or equal than b^-." + for i in xrange(len(a)): + if a[i]<0 and -a[i]>max(-b[i], 0): + return False + return True + +def disjoint_pp(v, w): + """true if v^+ and w^+ have disjoint support.""" + return not set(i for i in xrange(len(v)) if v[i]>0).intersection( + j for j in xrange(len(w)) if w[j]>0) + +def disjoint_mm(v, w): + """true if v^- and w^- have disjoint support.""" + # return not set(i for i in xrange(len(v)) if v[i]<0).intersection( + # j for j in xrange(len(w)) if w[j]<0) + return not bool(set(i for (i,vi) in enumerate(v) if vi<0).intersection( + j for (j,wi) in enumerate(w) if wi<0)) + +def max_ppp(z, u, w): + """check that z^+ ≤ u^+ ∨ v^+""" + for l in xrange(len(z)): + if z[l]>0 and z[l]>max(u[l],w[l],0): + return False + return True + # + # for (l,v) in enumerate(z): + # if v>0 and v>max(u[l],w[l],0): + # return False + # return True + +def max_mmm(z, u, w): + """check that z^- ≤ u^- ∨ v^-""" + for l in xrange(len(z)): + if z[l]<0 and -z[l]>max(-u[l],-w[l],0): + return False + return True + +def ne_max_max(z, u, w): + """check that z^+ ∨ u^+ ≠ u^+ ∨ w^+""" + for l in xrange(len(z)): + if max(z[l], u[l],0) != max(u[l], w[l],0): + return True + return False + +def ne_max_max_minus(z, u, w): + """check that z^- ∨ u^- ≠ u^- ∨ w^-""" + for l in xrange(len(z)): + if max(-z[l], -u[l],0) != max(-u[l], -w[l],0): + return True + return False + + + +def pm_split(v): + """ slightly faster than pm_split2""" + l = len(v) + up = l*[0] + um = l*[0] + for i,a in enumerate(v): + if a>0: + up[i]=a + else: + um[i]=-a + return vector(ZZ,up), vector(ZZ,um) + + +class PartialBasis(object): + """docstring for PartialBasis""" + def __init__(self, A, c, b, u): + super(PartialBasis, self).__init__() + self.d, self.n = A.dimensions() + self.A = A + self.c = c + self.b = b + self.u = u + self.vectors = [] + self.cv = [] + self.Av = [] + self.pairs = [] + # + self.vectors_p = [] + self.vectors_m = [] + self.Av_p = [] + self.Av_m = [] + # + self.zerox = vector(ZZ,self.n*[0]) + self.zerob = vector(ZZ, self.d*[0]) + for e in identity_matrix(ZZ,self.n).rows(): + self.add_element(copy(e)) + + + + def supports(self, v): + p, n = [], [] + for (i,c) in enumerate(v): + if c>0: + p.append(i) + elif c<0: + n.append(i) + return (set(p),set(n)) + + def binomial_order(self, b): + """docstring for binomial_order""" + c = b*self.c + minux = False + if c<0: + minux = True + elif c == 0: + for k in xrange(self.n): + l = self.n-k-1 + if b[l]>0: + minux = True + break + elif b[l]<0: + minux = False + break + else: + raise RuntimeError("Sorting the zero vector!! %s" % v) + if minux: + for (i,v) in enumerate(b): + b[i] = -v + return abs(c) + + + + def add_element(self, e, pos=None, pairs=True): + """ + add element e to partial basis, update critical pairs to be computed, + and precompute as much stuff as possible for reduction. + """ + c = self.binomial_order(e) + l = len(self.vectors) + Av = self.A*e + p,m = self.supports(e) + Ap,Am = self.supports(Av) + + if pos is None: + self.vectors.append(e) + self.cv.append(c) + self.Av.append(Av) + self.vectors_p.append(p) + self.vectors_m.append(m) + self.Av_p.append(Ap) + self.Av_m.append(Am) + else: + i = pos + self.vectors[i]=e + self.cv[i]=c + self.Av[i]=Av + self.vectors_p[i]=set(p) + self.vectors_m[i]=set(m) + self.Av_p[i]=Ap + self.Av_m[i]=Am + + + + if pairs and l: # true if this is not the first vector + for i in xrange(l): + self.pairs.append([i,l]) + + + def pop_element(self, i, pairs=False): + self.vectors.pop(i) + self.cv.pop(i) + self.Av.pop(i) + self.vectors_p.pop(i) + self.vectors_m.pop(i) + self.Av_m.pop(i) + self.Av_p.pop(i) + if pairs: + j = 0 + l = len(self.pairs) + while j=len(self.vectors):print 30*"-", "problem: %s at %s %s (len=%s)" % (i,j,p, len(self.vectors)) + + j += 1 + elif p[0]=len(self.vectors):print 30*"-", "problem: %s at %s %s (len=%s)" % (i,j,p, len(self.vectors)) + + j += 1 + else: + # print "nothing at %s, %s, %s" % (j,p,i) + j += 1 + + + # print len(self.vectors) + # print self.pairs + + + + + def criterion_1(self, i, j): + """True if (i,j) is skipabble: + The S-poly of the pair (i,j) is skippable if + its leading terms are disjoint. + """ + return not ( + bool(self.Av_p[i].intersection(self.Av_p[j])) or + bool(self.vectors_p[i].intersection(self.vectors_p[j])) or + bool(self.vectors_m[i].intersection(self.vectors_m[j])) + ) + + def criterion_3(self, i, j): + """Malkin's criterion.""" + return ( + bool(self.Av_m[i].intersection(self.Av_m[j])) or + bool(self.vectors_p[i].intersection(self.vectors_p[j])) or + bool(self.vectors_m[i].intersection(self.vectors_m[j])) + ) + + + def criterion_2(self, i, j): + """Gebauer and Möller""" + # it is assumed i0: + up[i]=a + elif a<0: + um[i]=-a + return up, um + + def feasible(self, v): + """docstring for feasible""" + if let_pp(v,self.u) and let_mp(v,self.u): + avp, avm = self.pm_split(v) + return let(self.A*avp, self.b) and let(self.A*avm, self.b) + return False + + +def bubu(A,b,c,u, interval=1000): + A = matrix(ZZ,A) + b = vector(ZZ, b) + c = vector(ZZ,c) + u = vector(ZZ,u) + m, n = A.dimensions() + conditions = [m == len(b), n == len(c), n == len(u) ] + + if not all(conditions): + raise ValueError("Not all conditions met %s" %conditions) + + # pdb.set_trace() + g = PartialBasis(A, c, b, u) + c1 = 0 + c2 = 0 + c3 = 0 + rx = 0 + elig = 0 + z = 0 + + + while g.pairs: # or unfinished? + # if len(g.vectors)%interval==0: + # print "cleaning (%s):" % len(g.vectors), + # g.clean() + # print "%s" % len(g.vectors) + + + if len(g.vectors)>2000: break + i, j = g.next_pair() + + if g.criterion_1(i,j): + c1 += 1 + continue + + if g.criterion_3(i,j): + c3 += 1 + continue + + # if g.criterion_2(i,j): # Gebauer and Möller criterion + # c2 += 1 + # continue + + + + # consider u[i]-u[j] + r = copy(g.zerox) + for l in xrange(g.n): + r[l] = g.vectors[j][l]-g.vectors[i][l] + g.binomial_order(r) + + + if g.feasible(r): + rx +=1 + if (rx%interval) == 0: + print "reduction: %s, basis: %s." % (rx, len(g.vectors)) + # g.clean() + # g.self_reduce() + c = g.reduce(r) + if not isz(c): + g.add_element(c) + else: + z += 1 + else: + # print r + elig += 1 + continue + + print "%s criterion 1, %s criterion 2, %s criterion 3, %s non feasible, %s reductions, of which %s to zero" % (c1, c2, c3, elig, rx, z) + return g.vectors, g + + + + +A1 = matrix(ZZ, [[3,1,11,2,3,5,3], [4, 5, 0, 1, 7,4, 6], [5,6,1,9,2,3,3]]) +c1 = vector(ZZ, [23,15,6,7,1,53,4]) +b1 = vector(ZZ, [31,27,38]) +u1 = vector(ZZ, 7*[38]) + + +B = matrix(ZZ, [ +[9, 6, 6, 8, 3, 6, 2, 4, 6, 3, 9, 4, 4, 8, 6, 9, 4, 2, 8, 8], +[8, 8, 0, 5, 4, 4, 7, 0, 2, 3, 0, 6, 7, 7, 0, 6, 7, 8, 6, 0], +[0, 9, 4, 1, 2, 0, 6, 5, 8, 5, 5, 5, 0, 0, 9, 3, 1, 8, 4, 8], +[5, 0, 1, 6, 0, 7, 7, 5, 8, 0, 2, 7, 3, 9, 0, 6, 9, 8, 7, 3], +[3, 0, 2, 7, 5, 8, 1, 2, 1, 5, 7, 8, 3, 8, 4, 3, 6, 2, 9, 3], +]) +c2 = vector(ZZ,[3, 2, 9, 1, 7, 5, 6, 2, 7, 8, 6, 9, 2, 6, 8, 7, 6, 2, 2, 1]) +b2 = vector(ZZ, 5*[20]) +u2 = vector(ZZ, 20*[1]) + diff --git a/wbtl.ipynb b/wbtl.ipynb new file mode 100644 index 0000000..3fb6d3e --- /dev/null +++ b/wbtl.ipynb @@ -0,0 +1,158 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The sage extension is already loaded. To reload it, use:\n", + " %reload_ext sage\n" + ] + } + ], + "source": [ + "%load_ext sage" + ] + }, + { + "cell_type": "code", + "execution_count": 134, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "%attach wbtl.py" + ] + }, + { + "cell_type": "code", + "execution_count": 135, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "127740" + ] + }, + "execution_count": 135, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "N = 100\n", + "L = 100000\n", + "cost = random_vector(ZZ, N, 50, 100)\n", + "def R(N): return random_vector(ZZ, N, -2, 20)\n", + "P = R(N)\n", + "T = [R(N) for t in range(L)]\n", + "Cost_so_far = 2*abs(fast_dot(cost, P))\n", + "Cost_so_far" + ] + }, + { + "cell_type": "code", + "execution_count": 136, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " " + ] + } + ], + "source": [ + "%prun W=list(worthy_children(P, T, cost, Cost_so_far))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.9" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/wbtl.py b/wbtl.py new file mode 100755 index 0000000..bd4167d --- /dev/null +++ b/wbtl.py @@ -0,0 +1,162 @@ +#!/usr/bin/env sage -python +# encoding: utf-8 + +""" +Algorithms for computing test sets of Integer Programming problems. + +# ============================================== +# = This is ongoing work by J. Soto y JM Ucha. = +# ============================================== + +This file implements the algorithm “walk back, then look”. + +""" + +# import operator # not needed anymore +# import itertools # not needed anymore +from sage.all import * + + + + +def fast_add(l,m): + """ + Fast addition of two lists. First versions used ideas in + + http://stackoverflow.com/questions/18713321/element-wise-addition-of-2-lists-in-python + + using, essentially, + + return map(operator.add, l,m) + + Slower alternatives: + + return [a+b for a,b in zip(l,m)] + return list(itertools.imap(operator.add, l,m)) + + Turns out, fastest is just to wrap Sage's addition. Can't reinvent the wheel, after all. + + >>> fast_add([],[]) + () + + >>> fast_add([1,2],[1]) + Traceback (most recent call last): + ... + TypeError: unsupported operand parent(s) for '+': 'Ambient free module of rank 2 over the principal ideal domain Integer Ring' and 'Ambient free module of rank 1 over the principal ideal domain Integer Ring' + + + Caveat: returns a `vector`; not a `list`. + + >>> fast_add(10000*[1],10000*[3]) == vector(ZZ,10000*[4]) + True + + + """ + return vector(ZZ,l)+vector(ZZ,m) + +# def isgtz(l): +# """Check whether l[i]>0 for all i.""" +# return all(li>0 for li in l) + +def isgetz(l): + """ + Given a list or vector `l`, return `True` if all its components are greater or equal than 0. Uses a generator and `all` for speed. + + >>> isgetz(1000*[1]) + True + >>> isgetz(1000*[1]+[-1]) + False + + + On the empty list it returns `True`. + >>> isgetz([]) + True + + """ + return all(li>=0 for li in l) + +def fast_dot(c,x): + """ + Quickly compute the scalar product `c*x`. + + First versions used some iterators ideas: + + return add(itertools.imap(operator.mul, c,x)) + return add(map(operator.mul, c,x)) + return add(ci*xi for ci,xi in zip(c,x)) + + but it turns out the fastest is to use Sage's vectors. + """ + return vector(ZZ,c)*vector(ZZ,x) + + +def worthy_children(P, T, cost, Oc=infinity): + """ + Find all children of point P from test set T, such that + 1. They are positive. + 2. Their cost is strictly less than the cost at the optimum so far, Oc. (Oc is the Optimum cost, not the optimum point itself.) + 3. Are feasible? + """ + for t in T: + candidate = fast_add(P,t) + if isgetz(candidate) and fast_dot(cost, candidate)<=Oc: + yield candidate + +def only_new(New, Cache): + """ + Check for each element of New whether it already belongs to cache. + Return a generator of elements not in the Cache, and update Cache. + + """ + for P in New: + if P not in Cache: + Cache.append(P) + yield P + +def add_to_V(current, T, c, best_cost_Omega, V, cached): + for v in worthy_children(current, T, c, best_cost_Omega): + if v not in cached: + cached.append(v) + V.append(v) + + +def walk_back_and_look_behind(P, T, Omega, c, A, b): + """ + Implement the walk back and look behind algorithm for the problem + min cx + st Ax=b + x in RegionOmega + with the given test set T, optimal point P for the linear problem, and non-linear region function indicator Omega. That is, x in RegionOmega == True if and only if Omega(x)==True + """ + best_cost_Omega = infinity + V = list() + S = list() + + if Omega(P): + S = [P] + else: + V = list(worthy_children(P, T, c)) + + cached = V[:] # shallow copy + while V: + current = V.pop() # 0 = first, nothing = last. + current_cost = fast_dot(c, current) + if Omega(current): + if current_cost < best_cost_Omega: + S = [current] + best_cost_Omega = current_cost + # elif current_cost == best_cost_Omega: + else: # worthy_children() only gives points with cost <= thatn best_cost_Omega + S.append(current) + enlarge = True + else: + enlarge = True + if enlarge: + # V += list(only_new(worthy_children(current, T, c, best_cost_Omega), cached)) + add_to_V(current, T, c, best_cost_Omega, V, cached) + enlarge = False + return S + +if __name__ == '__main__': + import doctest + doctest.testmod() \ No newline at end of file