Skip to content
New issue

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

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

Already on GitHub? Sign in to your account

improved integer vectors efficiency -Enhancement #36830

Merged
Merged
Changes from 8 commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
dc62757
added cardinality method to IntegerVectors_nk class
amanmoon Dec 6, 2023
95a984e
added an efficient cardinality function in IntegerVectors_nk and fixe…
amanmoon Dec 7, 2023
6c8c167
added rank, unrank and cardinality function in IntegerVectors, Intege…
amanmoon Dec 7, 2023
42b5620
fixed typos and made code more redable
amanmoon Dec 7, 2023
db8dbd2
fixed more typos and made code redable
amanmoon Dec 7, 2023
f536d45
changed EXAMPLE to EXAMPLES
amanmoon Dec 8, 2023
7fe22d1
made elif to if and added new EXAMPLES in IntegerVector_nk.cardinality()
amanmoon Dec 8, 2023
167acc9
formated Integer_vector.py according to sagemath guidelines
amanmoon Dec 8, 2023
cad751d
created helper function for unrank in IntegerVector
amanmoon Dec 8, 2023
48c188d
fixed the IndexError in IntegerVectors_n unrank method and made the I…
amanmoon Dec 9, 2023
ec8cb35
changed the name of unrank(helper) function and removed unnecessary e…
amanmoon Dec 12, 2023
2f1e749
added the else statement
amanmoon Dec 12, 2023
314c257
Merge branch 'sagemath:develop' into enhancement/IntegerVectors_effic…
amanmoon Dec 15, 2023
246b6d6
added examples in _unrank_helper method and corner cases test
amanmoon Dec 15, 2023
fcfb917
Merge branch 'enhancement/IntegerVectors_efficiency' of into enhancem…
amanmoon Dec 15, 2023
2fb35a6
corrected typos and indentation
amanmoon Dec 15, 2023
847b705
put the description of _unrank_helper at start of the method
amanmoon Dec 15, 2023
dedee8d
implemented alternate way of looping in the rank method
amanmoon Dec 15, 2023
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
249 changes: 230 additions & 19 deletions src/sage/combinat/integer_vector.py
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@
from sage.rings.integer import Integer


def is_gale_ryser(r,s):
def is_gale_ryser(r, s):
r"""
Tests whether the given sequences satisfy the condition
of the Gale-Ryser theorem.
Expand Down Expand Up @@ -314,20 +314,20 @@ def gale_ryser_theorem(p1, p2, algorithm="gale",
"""
from sage.matrix.constructor import matrix

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

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

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

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

Expand All @@ -340,12 +340,12 @@ def gale_ryser_theorem(p1, p2, algorithm="gale",
k = i + 1
while k < m and r[i] == r[k]:
k += 1
if t >= k - i: # == number rows of the same length
if t >= k - i: # == number rows of the same length
for j in range(i, k):
r[j] -= 1
c[j] = 1
t -= k - i
else: # Remove the t last rows of that length
else: # Remove the t last rows of that length
for j in range(k-t, k):
r[j] -= 1
c[j] = 1
Expand All @@ -366,17 +366,17 @@ def gale_ryser_theorem(p1, p2, algorithm="gale",
k1, k2 = len(p1), len(p2)
p = MixedIntegerLinearProgram(solver=solver)
b = p.new_variable(binary=True)
for (i,c) in enumerate(p1):
p.add_constraint(p.sum([b[i,j] for j in range(k2)]) == c)
for (i,c) in enumerate(p2):
p.add_constraint(p.sum([b[j,i] for j in range(k1)]) == c)
for (i, c) in enumerate(p1):
p.add_constraint(p.sum([b[i, j] for j in range(k2)]) == c)
for (i, c) in enumerate(p2):
p.add_constraint(p.sum([b[j, i] for j in range(k1)]) == c)
p.set_objective(None)
p.solve()
b = p.get_values(b, convert=ZZ, tolerance=integrality_tolerance)
M = [[0]*k2 for i in range(k1)]
for i in range(k1):
for j in range(k2):
M[i][j] = b[i,j]
M[i][j] = b[i, j]
return matrix(M)

else:
Expand Down Expand Up @@ -839,7 +839,10 @@ def __init__(self, n):
sage: TestSuite(IV).run()
"""
self.n = n
IntegerVectors.__init__(self, category=InfiniteEnumeratedSets())
if self.n == 0:
IntegerVectors.__init__(self, category=EnumeratedSets())
else:
IntegerVectors.__init__(self, category=InfiniteEnumeratedSets())

def _repr_(self):
"""
Expand Down Expand Up @@ -898,6 +901,83 @@ def __contains__(self, x):
return False
return sum(x) == self.n

def rank(self, x):
"""
Return the rank of a given element.

INPUT:

- ``x`` -- a list with ``sum(x) == n``

EXAMPLES::

sage: IntegerVectors(n=5).rank([5,0])
1
sage: IntegerVectors(n=5).rank([3,2])
3
"""
if sum(x) != self.n:
raise ValueError("argument is not a member of IntegerVectors({},{})".format(self.n, None))

n, k = self.n, len(x)

r = binomial(k+n-1, n+1)

for i in range(k - 1):
k -= 1
n -= x[i]
r += binomial(k + n - 1, n - 1)

This comment was marked as resolved.

This comment was marked as resolved.

return r

def unrank(self, x):
"""
Return the element at given rank x.

INPUT:

- ``x`` -- an integer.

EXAMPLES::

sage: IntegerVectors(n=5).unrank(2)
[4, 1]
sage: IntegerVectors(n=10).unrank(10)
[1, 9]
"""
ptr = 0
rtn = [self.n]
while self.rank(rtn) < x:
rtn.append(0)
rtn.pop()

while True:
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe this while loop can be factored out into a helper function in IntegerVectors?

In any case, we should store current_rank = self.rank(rtn), so we don't compute it twice.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

should i do the same for rank() method?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think no. If I see it correctly, it would only be 3 lines there.

if self.rank(rtn) < x:
rtn[ptr+1] = rtn[ptr]
rtn[ptr] = 0
ptr += 1
elif self.rank(rtn) > x:
rtn[ptr] -= 1
rtn[ptr-1] += 1
else:
return self._element_constructor_(rtn)

def cardinality(self):
"""
Return the cardinality of ``self``.

EXAMPLES::

sage: IntegerVectors(n=0).cardinality()
1
sage: IntegerVectors(n=10).cardinality()
+Infinity
"""
if self.n == 0:
return Integer(1)
else:
return PlusInfinity()


class IntegerVectors_k(UniqueRepresentation, IntegerVectors):
"""
Expand All @@ -912,7 +992,10 @@ def __init__(self, k):
sage: TestSuite(IV).run()
"""
self.k = k
IntegerVectors.__init__(self, category=InfiniteEnumeratedSets())
if self.k == 0:
IntegerVectors.__init__(self, category=EnumeratedSets())
else:
IntegerVectors.__init__(self, category=InfiniteEnumeratedSets())

def _repr_(self):
"""
Expand Down Expand Up @@ -968,6 +1051,86 @@ def __contains__(self, x):
return False
return len(x) == self.k

def rank(self, x):
"""
Return the rank of a given element.

INPUT:

- ``x`` -- a list with ``len(x) == k``

EXAMPLES::

sage: IntegerVectors(k=5).rank([0,0,0,0,0])
0
sage: IntegerVectors(k=5).rank([1,1,0,0,0])
7
"""
if len(x) != self.k:
raise ValueError("argument is not a member of IntegerVectors({},{})".format(None, self.k))

n, k = sum(x), self.k

r = sum(binomial(k + i - 1, k - 1) for i in range(n))
Copy link
Contributor

@jukkakohonen jukkakohonen Dec 9, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It may be useful to notice that
sum(binomial(k + i - 1, k - 1) for i in range(n)) == binomial(n+k-1, k)
avoiding summation (please check yourself before using).

Currently it seems that IntegerVectors(k=1)[n] takes time roughly quadratic to n, which is worse than the old linear version.

sage: V=IntegerVectors(k=1)
....: %time V[1000]
CPU times: user 791 ms, sys: 0 ns, total: 791 ms
Wall time: 791 ms
[1000]
sage: V=IntegerVectors(k=1)
....: %time V[2000]
CPU times: user 3.13 s, sys: 0 ns, total: 3.13 s
Wall time: 3.13 s
[2000]
sage: V=IntegerVectors(k=1)
....: %time V[4000]
CPU times: user 10.9 s, sys: 0 ns, total: 10.9 s
Wall time: 10.9 s
[4000]

How does it fare if you replace the summation with a single binomial?


for i in range(k - 1):
k -= 1
n -= x[i]
r += binomial(k + n - 1, n - 1)

return r

def unrank(self, x):
"""
Return the element at given rank x.

INPUT:

- ``x`` -- an integer such that x < len(self) ``

EXAMPLES::

sage: IntegerVectors(k=5).unrank(10)
[1, 0, 0, 0, 1]
sage: IntegerVectors(k=5).unrank(15)
[0, 0, 2, 0, 0]
"""
if self.k == 0 and x != 0:
raise IndexError(f"Index {x} is out of range for the IntegerVector.")
else:
Copy link
Collaborator

@mantepse mantepse Dec 11, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

else: is unnecessary here, because if can only raise an error. Similarly below.

n, ptr = 0, 0
rtn = [0]*self.k
while self.rank(rtn) <= x:
n += 1
rtn[ptr] = n
rtn[ptr] -= 1
while True:
if self.rank(rtn) < x:
rtn[ptr+1] = rtn[ptr]
rtn[ptr] = 0
ptr += 1
elif self.rank(rtn) > x:
rtn[ptr] -= 1
rtn[ptr-1] += 1
else:
return self._element_constructor_(rtn)

def cardinality(self):
"""
Return the cardinality of ``self``.

EXAMPLES::

sage: IntegerVectors(k=0).cardinality()
1
sage: IntegerVectors(k=10).cardinality()
+Infinity
"""
if self.k == 0:
return Integer(1)
else:
return PlusInfinity()


class IntegerVectors_nk(UniqueRepresentation, IntegerVectors):
"""
Expand Down Expand Up @@ -1010,11 +1173,11 @@ def _list_rec(self, n, k):
res = []

if k == 1:
return [ (n, ) ]
return [(n, )]

for nbar in range(n + 1):
n_diff = n - nbar
for rest in self._list_rec( nbar , k - 1):
for rest in self._list_rec(nbar, k - 1):
res.append((n_diff,) + rest)
return res

Expand Down Expand Up @@ -1160,10 +1323,58 @@ def rank(self, x):
for i in range(k - 1):
k -= 1
n -= x[i]
r += binomial(k + n - 1, k)
r += binomial(k + n - 1, n - 1)

return r

def unrank(self, x):
"""
Return the element at given rank x.

INPUT:

- ``x`` -- an integer such that ``x < len(self)``

EXAMPLES::

sage: IntegerVectors(4,5).unrank(30)
[1, 0, 1, 0, 2]
sage: IntegerVectors(2,3).unrank(5)
[0, 0, 2]
"""
if x >= len(self):
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using len(self) seems problematic. For me

sage: V=IntegerVectors(100000,6); c=V.cardinality(); c; V[c-1]

produces

83345834041685416895001
...
OverflowError: cannot fit 'int' into an index-sized integer

How about using cardinality instead?

raise IndexError(f"Index {x} is out of range for the IntegerVector.")
else:
ptr = 0
rtn = [0]*self.k
rtn[ptr] = self.n
while True:
if self.rank(rtn) < x:
rtn[ptr+1] = rtn[ptr]
rtn[ptr] = 0
ptr += 1
elif self.rank(rtn) > x:
rtn[ptr] -= 1
rtn[ptr-1] += 1
else:
return self._element_constructor_(rtn)

def cardinality(self):
"""
Return the cardinality of ``self``.

EXAMPLES::

sage: IntegerVectors(3,5).cardinality()
35
sage: IntegerVectors(99, 3).cardinality()
5050
sage: IntegerVectors(10^9 - 1, 3).cardinality()
500000000500000000
"""
n, k = self.n, self.k
return Integer(binomial(n + k - 1, n))


class IntegerVectors_nnondescents(UniqueRepresentation, IntegerVectors):
r"""
Expand Down Expand Up @@ -1320,11 +1531,11 @@ def __init__(self, n=None, k=None, **constraints):
category = FiniteEnumeratedSets()
else:
category = EnumeratedSets()
elif k is not None and 'max_part' in constraints: # n is None
elif k is not None and 'max_part' in constraints: # n is None
category = FiniteEnumeratedSets()
else:
category = EnumeratedSets()
IntegerVectors.__init__(self, category=category) # placeholder category
IntegerVectors.__init__(self, category=category) # placeholder category

def _repr_(self):
"""
Expand Down