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

Some bugfixes and improvements #3672

Merged

Conversation

HechtiDerLachs
Copy link
Collaborator

This is work in progress for all the small fixes necessary to make some old scripts of @simonbrandhorst and myself run again.

@HechtiDerLachs HechtiDerLachs marked this pull request as draft May 2, 2024 11:33
Comment on lines 1119 to 1156
# Try to eliminate variables first. This will often speed up computations significantly.
if simplify_ring
Q, pr = quo(R, I)
W, id, id_inv = simplify(Q)
@assert domain(id) === codomain(id_inv) === Q
@assert codomain(id) === domain(id_inv) === W
res_simp = minimal_primes(modulus(W); algorithm, factor_generators, simplify_ring=false)
result = [I + ideal(R, lift.(id_inv.(W.(gens(j))))) for j in res_simp]
for p in result
end
return result
end

# This will in many cases lead to an easy simplification of the problem
if factor_generators
J = typeof(I)[ideal(R, elem_type(R)[])]
for g in gens(I)
K = typeof(I)[]
is_zero(g) && continue
for (b, k) in factor(g)
for j in J
push!(K, j + ideal(R, b))
end
end
J = K
end
result = unique!(filter!(!is_one, vcat([minimal_primes(j; algorithm, factor_generators=false) for j in J]...)))
# The list might not consist of minimal primes only. We have to discard the embedded ones
final_list = typeof(I)[]
for p in result
any((q !== p && is_subset(q, p)) for q in result) && continue
push!(final_list, p)
end
for p in final_list
set_attribute!(p, :is_prime=>true)
end
return final_list
end
Copy link
Collaborator Author

Choose a reason for hiding this comment

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

@afkafkafk13 , @hannes14 : Could you let me knew what you think of these? I tried to run a script again from about 5 months ago and it broke down on some primary decompositions (which used to go though at some point). Now I'm trying to get back there with these tweaks. And I thought that they might be useful in general.

@hannes14 : This was also used to prepare the example I sent you. Interestingly it is solved rather quickly with algorithm=:charSets, but never terminated for me with :GTZ.

Comment on lines 1127 to 1128
for p in result
end
Copy link
Collaborator

Choose a reason for hiding this comment

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

This looks a bit weird: either some line was lost or some leftovers were not deleted

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

Indeed. I will clean that up.

@afkafkafk13
Copy link
Collaborator

General question: Do you have any idea, what increased the cost of these older scripts? I find this a bit scary!
To me this sounds like some optimization (at some other and maybe completely unrelated) point having very negative (and up to now undetected) side effects.

The improvements make sense in any case. I have used such tweaks many times for large settings in Singular, where primary decomposition as a black box did not go through. But you need to be aware that you are deliberately circumventing the heuristics of primdec.lib and this should only be necessary in nasty examples. Is your offending example really sufficiently nasty?

@HechtiDerLachs
Copy link
Collaborator Author

Thanks for your comments.

No, in my opinion the offending example is not nasty. I sent it to @hannes14 via email, but I can put it here again:

julia> minimal_primes(list[2])
R = Multivariate polynomial ring in 4 variables over QQ
gens(R) = QQMPolyRingElem[θ_2, θ_1, (s0//s1), (w//x)]
g = -θ_1^2*(s0//s1) + θ_1^2 + (s0//s1)*(w//x)^4
g = θ_2^8 + 1
g = -θ_2^6 + θ_2^4 - θ_2^2 + θ_1^2

Be aware that s0//s1 etc. are not fractions, but only symbols for variable names. So this is really a polynomial computation and the three gs are the generators of the ideal in question (printed in a for-loop with g as local variable).

To me this did not look scary at all. But minimal_primes gives up on this with GTZ.

@HechtiDerLachs
Copy link
Collaborator Author

Does anyone have an explanation for the failing tests? I suppose the same things are tested in other threads, aren't they? Why do they fail on ubuntu 1.10 long?

@lgoettgens ?

@lgoettgens
Copy link
Member

Does anyone have an explanation for the failing tests? I suppose the same things are tested in other threads, aren't they? Why do they fail on ubuntu 1.10 long?

@lgoettgens ?

ubuntu-1.10-long spent 2.5 hours in /home/runner/work/Oscar.jl/Oscar.jl/test/AlgebraicGeometry/Schemes/elliptic_surface.jl and timed out

withoutexperimental spent about 2 hours in /home/runner/work/Oscar.jl/Oscar.jl/test/AlgebraicGeometry/Schemes/elliptic_surface.jl and timed out

ubuntu-1.10-book failed due to some change in output in decker-schmitt-invariant-theory: https://github.com/oscar-system/Oscar.jl/actions/runs/8927137202/job/24519841523?pr=3672#step:8:3632

@afkafkafk13
Copy link
Collaborator

afkafkafk13 commented May 2, 2024

Thanks for your comments.

No, in my opinion the offending example is not nasty. I sent it to @hannes14 via email, but I can put it here again:

julia> minimal_primes(list[2])
R = Multivariate polynomial ring in 4 variables over QQ
gens(R) = QQMPolyRingElem[θ_2, θ_1, (s0//s1), (w//x)]
g = -θ_1^2*(s0//s1) + θ_1^2 + (s0//s1)*(w//x)^4
g = θ_2^8 + 1
g = -θ_2^6 + θ_2^4 - θ_2^2 + θ_1^2

Be aware that s0//s1 etc. are not fractions, but only symbols for variable names. So this is really a polynomial computation and the three gs are the generators of the ideal in question (printed in a for-loop with g as local variable).

To me this did not look scary at all. But minimal_primes gives up on this with GTZ.

I entered the computation in Singular and found the following observations:
ring r=0,(a,b,c,d),dp; ideal I=b^2*c+b^2+c*d^4, a^8+1, -a^6+a^4-a^2+b^2;
internally uses a ring (0,c),(a,b,d),(lp(3),C) and stalls

ring r=0,(d,c,b,a),dp; ideal I=b^2*c+b^2+c*d^4, a^8+1, -a^6+a^4-a^2+b^2;
internally uses a ring (0,d),(c,b,a),(lp(3),C) and returns the result in a fraction of a second.

So the difference is in the internal choice of the "maximal independent set" in this primary decomposition algorithm and of course in the difficulty of the lp-groebner basis to be computed internally. This choice is not done by an elaborate heuristic; hence the situation is brittle, whenever the different choices have significantly different runtimes. So probalby some change further up in the hierarchy changed the sequence of the variables, which in turn triggers a different choice of a maximal independent set.
Actually, the most solid (but currently rather unrealistic) variant would be to run the different choices in parallel and let the first one win. (Maybe you want to at least make sure that the variables from the tower of algebraic field extensions are added in a suitable order to avoid unnecessary spoly among the minimal polynomials w.r.t. lex ordering.)

@HechtiDerLachs
Copy link
Collaborator Author

I changed the default algorithm for the place where primary decomposition is called in the elliptic surfaces. Maybe that already helps to have the tests run more reliably (#3676).

@HechtiDerLachs HechtiDerLachs marked this pull request as ready for review May 3, 2024 13:27
Copy link
Collaborator

@afkafkafk13 afkafkafk13 left a comment

Choose a reason for hiding this comment

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

Thank you.

Concerning the long-running tests: It makes sense to keep them, but disable them. We should, however, remember to run them locally for all PR concerning optimizations in the basic functionality for IdealSheaf, ..Divisor etc to detect regression.


L, loc_map = localization(W, U)

I = ideal(L, [y^2, 0]) # Dealing with zero generators lead to a bug once, so this is a real test.
Copy link
Collaborator

Choose a reason for hiding this comment

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

Suggested change
I = ideal(L, [y^2, 0]) # Dealing with zero generators lead to a bug once, so this is a real test.
I = ideal(L, [y^2, 0]) # Dealing with zero generators led to a bug once, so this is a real test.

@joschmitt
Copy link
Member

The change to minimal_primes changes the ordering in one of the book tests. How do we handle this again?
Otherwise I would like to have this in, so that the CI stops running into #3676 (I hope).

@benlorenz
Copy link
Member

The change to minimal_primes changes the ordering in one of the book tests. How do we handle this again? Otherwise I would like to have this in, so that the CI stops running into #3676 (I hope).

You can change the output line in the jlcon file (test/book/specialized/decker-schmitt-invariant-theory/H3.jlcon).

Comment on lines +879 to +884
function _maximal_associated_points(
I::AbsIdealSheaf;
covering=default_covering(scheme(I)),
use_decomposition_info::Bool=true,
algorithm::Symbol=:GTZ
)
Copy link
Collaborator

Choose a reason for hiding this comment

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

Quoting the singular manual:

https://www.singular.uni-kl.de/Manual/latest/sing_1426.htm#SEC1507
"In small characteristic and for algebraic extensions, the procedures via Gianni, Trager, Zacharias may not terminate. "

Copy link
Collaborator

Choose a reason for hiding this comment

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

Copy link
Collaborator

Choose a reason for hiding this comment

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

The algebraic extensions only refer to those implemented as extended basefield in Singular. If you encode the algebraic extension by hand, the problem does not occur (hence the encoding by hand). This is the reason why these particular calls to Singular should not use the algebraic extensions in Singular, but do the encoding themselves (this was subject of discussion in another PR).
The restriction on the characteristic results from the sufficiently general linear coordinate change, which is used (the open set needs to be non-empty, which might be a problem in very small characteristics). This should be stated as a warning somewhere in the docstring, but allowing characteristics like 101 or 32003 is vital for many practical applications where coefficients are known to explode in characteristic zero.

@joschmitt
Copy link
Member

I adjusted the output of the failing book test.

Copy link

codecov bot commented May 7, 2024

Codecov Report

Attention: Patch coverage is 69.11765% with 21 lines in your changes are missing coverage. Please review.

Project coverage is 81.31%. Comparing base (789bd83) to head (70da52f).
Report is 25 commits behind head on master.

Additional details and impacted files
@@            Coverage Diff             @@
##           master    #3672      +/-   ##
==========================================
- Coverage   82.89%   81.31%   -1.59%     
==========================================
  Files         577      577              
  Lines       78308    78568     +260     
==========================================
- Hits        64916    63886    -1030     
- Misses      13392    14682    +1290     
Files Coverage Δ
experimental/Schemes/critical_locus.jl 0.00% <ø> (ø)
src/Rings/mpoly-ideals.jl 93.02% <100.00%> (+0.36%) ⬆️
src/Rings/mpolyquo-localizations.jl 71.11% <100.00%> (-0.44%) ⬇️
src/Rings/MPolyMap/MPolyQuo.jl 88.88% <33.33%> (-2.78%) ⬇️
experimental/Schemes/elliptic_surface.jl 4.44% <0.00%> (-70.64%) ⬇️
experimental/Schemes/IdealSheaves.jl 72.21% <36.36%> (-5.85%) ⬇️
...erimental/Schemes/MorphismFromRationalFunctions.jl 39.93% <20.00%> (-7.53%) ⬇️

... and 67 files with indirect coverage changes

Comment on lines +24 to +25
# We need to avoid creating a polynomial ring with zero variables
dim(X) == 0 && vector_space_dimension(OO(X)) == 0 && return SimplifiedAffineScheme(X, X, identity_map(X), identity_map(X), check=false)
Copy link
Collaborator

Choose a reason for hiding this comment

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

May I suggest moving the point to 0? There is the function rational_point_coordinates.
It will push some of the complexity into the gluings ... but things like saturation should be easier.
Not sure though how exactly simplify is hooked into everything.

Copy link
Collaborator Author

Choose a reason for hiding this comment

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

That's actually not that easy to program. I don't think it pays off since these schemes are already very simple and stuff like saturation will not be applied here.

@simonbrandhorst
Copy link
Collaborator

I would suggest to merge this and relegate further improvements to a new PR.

@HechtiDerLachs
Copy link
Collaborator Author

Good to go from my side.

@afkafkafk13 afkafkafk13 merged commit 671c689 into oscar-system:master May 8, 2024
26 checks passed
joschmitt added a commit to joschmitt/Oscar.jl that referenced this pull request May 8, 2024
simonbrandhorst pushed a commit that referenced this pull request May 10, 2024
This reverts commit 85594b1.

Commit is made obsolete by PR #3672
(85594b1)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

6 participants