oscar-system / Singular.jl

Julia package for the Singular library
Other
31 stars 33 forks source link

Access to Singular interpreter variables from the kernel? #595

Closed fingolfin closed 1 year ago

fingolfin commented 1 year ago

Currently @LukasKuehne and @bschroter are trying to use Homalg with OSCAR. For various reasons I don't want to go into right now, this requires access to singular interpreter variables.

In a nutshell, consider this trivial example:

julia> Singular.libSingular.call_interpreter("int a = 42;")
4-element Vector{Any}:
 false
      ""
      ""
      ""

Now, how can I access the value of the Singular interpreter variable a from Julia? Of course we can ask the Singular interpreter for it, but all that gives us is a string representation:

julia> Singular.libSingular.call_interpreter("a;")
4-element Vector{Any}:
 false
      "42\n"
      ""
      ""

Of course an int is a bit special as it is ring independent. In practice I think they need to access polynomials and singular matrices (over which domains?), and I guess thus also rings. So besides wanting to get those from the interpreter, we'd also have to wrap them properly into Singular.jl objects...

@thofma pointed me at lookup_library_symbol but this calls the C++ function lookup_singular_library_symbol_wo_rng and it also requires a package name, which is not relevant here.

In the SingularInterface code, we used this helper function to this end:

Obj FuncSingularValueOfVar(Obj self, Obj name)
{
    Int      len;
    Obj      tmp, tmp2;
    intvec * v;
    int      i, j, k;
    Int      rows, cols;
    /* number n;   */

    idhdl h = ggetid(reinterpret_cast<char *>(CHARS_STRING(name)));
    if (h == NULL)
        return Fail;
    switch (IDTYP(h)) {
    case INT_CMD:
        return ObjInt_Int((Int)(IDINT(h)));
    case STRING_CMD:
        return MakeString(IDSTRING(h));
    case INTVEC_CMD:
        v = IDINTVEC(h);
        len = (Int)v->length();
...
}
fingolfin commented 1 year ago

Ping @tthsqe12 @hannes14

tthsqe12 commented 1 year ago

I think there is always a package. At the top level it is called "main" or something like that. Will have an answer in 9 hours.

fingolfin commented 1 year ago

Thanks @tthsqe12. Too me some sleuthing around, but this seems to work (at least in this simple case):

julia> Singular.lookup_library_symbol("Top", "a")
42

Great! Next, I tried creating a ring in the interpreter, which also worked:

julia> Singular.libSingular.call_interpreter("ring r=0,(x,y,z),dp;")
4-element Vector{Any}:
 false
      ""
      ""
      ""

julia> Singular.lookup_library_symbol("Top", "r")
2-element Vector{Any}:
 Singular Polynomial Ring (QQ),(x,y,z),(dp(3),C)
 Dict{Any, Any}()

Then I wanted to extra a polynomial in that ring:

julia> Singular.libSingular.call_interpreter("poly f = 2x2z - 42y;")
4-element Vector{Any}:
 false
      ""
      ""
      ""

julia> Singular.libSingular.call_interpreter("f;")
4-element Vector{Any}:
 false
      "2x2z-42y\n"
      ""
      ""

julia> Singular.lookup_library_symbol("Top", "f")
ERROR: Singular symbol Top::f does not exist

Hmmm, OK... so the internal name of the C++ function in libsingular-polymake we call gave me a hint: it is called lookup_singular_library_symbol_wo_rng. And IIRC, ring-dependant variables are stored "inside the ring". So let's try that:

julia> Singular.lookup_library_symbol("r", "f")
ERROR: MethodError: no method matching base_ring(::Nothing)
Closest candidates are:
  base_ring(::Union{AbstractAlgebra.Generic.LaurentSeriesFieldElem{T}, AbstractAlgebra.Generic.LaurentSeriesRingElem{T}} where T<:RingElement) at ~/.julia/packages/AbstractAlgebra/ZmWFo/src/generic/LaurentSeries.jl:40
  base_ring(::Union{AbstractAlgebra.Generic.PuiseuxSeriesFieldElem{T}, AbstractAlgebra.Generic.PuiseuxSeriesRingElem{T}} where T<:RingElement) at ~/.julia/packages/AbstractAlgebra/ZmWFo/src/generic/PuiseuxSeries.jl:56
  base_ring(::T) where T<:Union{gfp_mat, nmod_mat} at ~/.julia/packages/Nemo/pzCmt/src/flint/nmod_mat.jl:126
  ...
Stacktrace:
 [1] (::Singular.var"#184#200")(vptr::Ptr{Nothing}, R::Nothing)
   @ Singular ~/.julia/packages/Singular/hfI7H/src/caller.jl:50
 [2] convert_normal_value(valueptr::Ptr{Nothing}, typ::Int64, R::Nothing)
   @ Singular ~/.julia/packages/Singular/hfI7H/src/caller.jl:152
 [3] convert_return(value::Vector{Any}, R::Nothing)
   @ Singular ~/.julia/packages/Singular/hfI7H/src/caller.jl:223
 [4] lookup_library_symbol(pack::String, name::String)
   @ Singular ~/.julia/packages/Singular/hfI7H/src/caller.jl:402
 [5] top-level scope
   @ REPL[41]:1

OK, it gets further, but fails because the parent ring is not known to the system. Indeed, lookup_library_symbol does this: return convert_return(res, nothing) and we'd have to sneak a different value in there instead of nothing.

But then I also discovered this:

julia> Singular.lookup_library_symbol("Top", "r")
2-element Vector{Any}:
 Singular Polynomial Ring (QQ),(x,y,z),(dp(3),C)
 Dict{Symbol, Singular.spoly{Singular.n_Q}}(:f => 2*x^2*z - 42*y)

julia> Singular.libSingular.call_interpreter("poly g = xyz3;")
4-element Vector{Any}:
 false
      ""
      ""
      ""

julia> Singular.lookup_library_symbol("Top", "r")
2-element Vector{Any}:
 Singular Polynomial Ring (QQ),(x,y,z),(dp(3),C)
 Dict{Symbol, Singular.spoly{Singular.n_Q}}(:f => 2*x^2*z - 42*y, :g => x*y*z^3)

So it seems asking for the ring also gives us all the polynomials defined in there? That could be handy (although in our case, we really only have the name of the polynomial... though I guess we could figure out the name of the parent ring, too, with some effort).

Alas, this quickly seems to run into limitations of the system: when I tried to work with an ideal, I got a crash:

julia> Singular.libSingular.call_interpreter("ideal I = f,g;")
4-element Vector{Any}:
 false
      ""
      ""
      ""

julia> Singular.lookup_library_symbol("Top", "r")
2-element Vector{Any}:
julia-1.8(96953,0x103044580) malloc: *** error for object 0x89804: pointer being freed was not allocated
julia-1.8(96953,0x103044580) malloc: *** set a breakpoint in malloc_error_break to debug

signal (6): Abort trap: 6
in expression starting at none:0
__pthread_kill at /usr/lib/system/libsystem_kernel.dylib (unknown line)
Allocations: 82218589 (Pool: 82140439; Big: 78150); GC: 32
Abort trap: 6
tthsqe12 commented 1 year ago

The crash shouldn't happen. Will investigate. I would also think that Singular.lookup_library_symbol("r", "f") should either work or fail much sooner.

hannes14 commented 1 year ago

Trying r,D=Singular.lookup_library_symbol("Top", "r") one gets the ring r (as julia object) and all the (global) objects which have r as parent in form of a Dict

tthsqe12 commented 1 year ago

This one is fun. here is a MNWE. Including this code consistently executes without crash

using Revise, Singular

Singular.libSingular.call_interpreter("ring r=0,(x,y,z),dp;")
z = Singular.lookup_library_symbol("Top", "r")

Singular.libSingular.call_interpreter("poly f = x;")
z = Singular.lookup_library_symbol("Top", "r")

Singular.libSingular.call_interpreter("ideal I = y;")
z = Singular.lookup_library_symbol("Top", "r")

nothing

But trying to look at z consistently crashes.

julia> z
2-element Vector{Any}:

signal (11): Segmentation fault
tthsqe12 commented 1 year ago

Very interesting that this serious issue has not reared its ugle head yet. First, here is what can happen

Singular.libSingular.call_interpreter("ring r=0,(x,y,z),dp;")
z = Singular.lookup_library_symbol("Top", "r")

# all good so far

Singular.libSingular.call_interpreter("poly f = x;")
z = Singular.lookup_library_symbol("Top", "r")

# At this point the poly `p` in `:f => p` is walking around in julia clothes even though
# the kernel pointer has NOT be coppied.
# Therefore, garbage collection of `p` will corrupt the poly that the Singular interpreter owns.

Singular.libSingular.call_interpreter("ideal I = y;")
z = Singular.lookup_library_symbol("Top", "r")

# The only thing that can happen here is said garbage collection occurs, and stuff crashes.

@hannes14 can you suggest an alternative to https://github.com/oscar-system/libsingular-julia/blob/e91e08ea4eec25a3fab876fbe414b7c7ab366016/caller.cpp#L252 that copies the data?

It looks like https://github.com/oscar-system/libsingular-julia/issues/39 finally needs addressing.

hannes14 commented 1 year ago

substitute jl_arrayset(current, jl_box_voidpointer(IDDATA(h)), 2); by sleftv x; x.Copy((leftv)h); jl_arrayset(current, jl_box_voidpointer(x.data), 2);

tthsqe12 commented 1 year ago

@hannes14 does the current ring to be set correct for this copy to work?

hannes14 commented 1 year ago

yes, currRing needs to be set to the base ring of these elements, for example: ring save=currRing;rChangeCurrRing(IDRING(h)); .... rChangeCurrRing(save);

tthsqe12 commented 1 year ago

@hannes14 thanks, seems to work. @fingolfin The crash seems to be fixed now. Is there anything else (besides possibly exporting lookup_library_symbol and adding its doc string to the manual) to address for this issue?

fingolfin commented 1 year ago

@tthsqe12 fixed how and where? your example still crashes for me in latest Singular.jl and I see no recent updates to it?

tthsqe12 commented 1 year ago

It is fixed on my computer :) I'll have new version of Singular.jl in Oscar.jl once some issues with the Oscar test code are resolved.

tthsqe12 commented 1 year ago

Closing this because the following session exhibits normal behaviour with Singular.jl version 0.14

julia> Singular.call_interpreter("int a = 42;")
julia> Singular.call_interpreter("a;")
julia> Singular.lookup_library_symbol("Top", "a")
julia> Singular.call_interpreter("ring r=0,(x,y,z),dp;")
julia> Singular.lookup_library_symbol("Top", "r")
julia> Singular.call_interpreter("poly f = 2x2z - 42y;")
julia> Singular.call_interpreter("f;")
julia> Singular.lookup_library_symbol("Top", "f")
julia> Singular.lookup_library_symbol("Top", "r")
julia> Singular.call_interpreter("poly g = xyz3;")
julia> Singular.lookup_library_symbol("Top", "r")
julia> Singular.call_interpreter("ideal I = f,g;")
julia> Singular.lookup_library_symbol("Top", "r")
fingolfin commented 1 year ago

@tthsqe12 thank you, works great for me, too!

@LukasKuehne @bschroter so that part is done. Have you managed to get the code by Mohamed to work which directly uses Singular.jl for computations? Then in combination with the code snippets by Dan above, you should now be in business...

LukasKuehne commented 1 year ago

Dear all, Thank you very much for your great support. It mostly works now for me too. Can I just ask two more things:

1) Your code, Dan, from above doesn't fully work for me yet: (I've updated Singular.jl to v0.14.0):

julia> Singular.call_interpreter("ring r=0,(x,y,z),dp;");
julia> Singular.lookup_library_symbol("Top", "r")
2-element Vector{Any}:
 Singular Polynomial Ring (QQ),(x,y,z),(dp(3),C)
 Dict{Any, Any}()

julia> Singular.call_interpreter("poly f = 2x2z - 42y;");
julia> Singular.call_interpreter("f;")
4-element Vector{Any}:
 false
      "2x2z-42y\n"
      ""
      ""
julia> Singular.lookup_library_symbol("Top", "f")
ERROR: Singular symbol Top::f does not exist
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] lookup_library_symbol(pack::String, name::String)
   @ Singular ~/.julia/packages/Singular/IjiQo/src/caller.jl:421
 [3] top-level scope
   @ REPL[87]:1

So it seems, I can look up rings (and get the ring elements from the dictionary that this look-up yields) but I can't look up the ring elements directly. I could live with that detour. Or am I missing something? Some other component I need to update?

2) Is there a direct way to cast the rings, polynomials and ideals I'm getting through this look up (or via the associated dictionary) to Oscar objects? (as our code lives in Oscar and mostly deals with objects in oscar so far.)

Thanks a lot again!

tthsqe12 commented 1 year ago

I haven't had time to look into the possibility of only looking up a particular symbol in a ring. Since rings themselves contain a symbol table, they behave much like a package and should allow for such lookup function. For now - at least until the next singular.jl release - you just have to copy the whole contents of the ring into a dictionary and then lookup the poly, ideal, module, etc. from the dictionary. As for the casting into an Oscar object, this is currently usualy done with the oscarring(singularpoly) syntax, but keep in mind that the Singular.jl objects do quite well on their own when confined to operations within the singular kernel.

fingolfin commented 1 year ago

To make what Dan said explicit: you can do

julia> R, d = Singular.lookup_library_symbol("Top", "r");;

julia> f = d[:f]
2*x^2*z - 42*y

julia> S, (x,y,z) = PolynomialRing(QQ,[:x,:y,:z])
(Multivariate Polynomial Ring in x, y, z over Rational Field, fmpq_mpoly[x, y, z])

julia> S(f)
2*x^2*z - 42*y

julia> typeof(ans)
fmpq_mpoly
LukasKuehne commented 1 year ago

Thank you both for your quick replies.

I haven't had time to look into the possibility of only looking up a particular symbol in a ring. Since rings themselves contain a symbol table, they behave much like a package and should allow for such lookup function. For now - at least until the next singular.jl release - you just have to copy the whole contents of the ring into a dictionary and then lookup the poly, ideal, module, etc. from the dictionary.

Sure that's no problem and quite convenient to use. I was just confused since you wrote Singular.lookup_library_symbol("Top", "f") in your code example above. That's why I thought one can also look up ring elements directly.

As for the casting into an Oscar object, this is currently usualy done with the oscarring(singularpoly) syntax, but keep in mind that the Singular.jl objects do quite well on their own when confined to operations within the singular kernel.

Okay cool, then I'll keep them as Singular.jl objects for now.

But just for curiosity/completeness: When I have a Singular polynomial ring, do I have to manually create the analogous Oscar ring (as in your example Max) with the same variables etc. or can I also somehow also copy this ring directly over to Oscar?

Thanks again!