Linear Algebra: Method selection

[fieker]

as an example: determinant of a matrix. Available algorithms: (at least a subset)

• for fields via rref or similar (n^3)
• for rings fraction free evaluation (n^4)
• for euclidean rings: HNF (n^3)
• for univariate polynomials: interpolation/ evaluation
• special implementations for fmpz_mat, fmpq_mat, ...

GIven a matrix aver a ring, how should the "correct" method be selected? There is no type-match for euclidean rings - nor can there be: We've just implemented localisations of (maximal) orders of number fields. In some cases, this is a euclidean ring Similar: quadratic fields: in some cases (important for crypto), the ring of integers is effective euclidean (Z[i]) but cannot distinguished by type.

What is a generic method for dealing with this?

[wbhart]

I don't think you want me to check if you have a Euclidean ring, do you? And it's not a field. It's not a polynomial ring. So I only see one method that can/will be applied.

Do you have a better algorithm in mind?

[fieker]

On Mon, Jun 17, 2019 at 09:31:51AM -0700, wbhart wrote:

I don't think you want me to check if you have a Euclidean ring, do you? And it's not a field. It's not a polynomial ring. So I only see one method that can/will be applied.

Do you have a better algorithm in mind?

I do not know how to generically select the "best" algorithm other than testing is isUsefulEuclidean is set. I do know your aversion to checks, but a determinant, say 100x100 is MUCH more expensive than isUsefulEuclidean. Otherwise, generic algorithms will be slow as they will only be able to select

• special implementation
• field
• ring

As I said, I do not have a good solution, but I have many good problems. How do you envision to call a function that does an hnf for a matrix over Z[i], secure in the knowledge that Z[i] is a euclidean domain?

Or is the intend to always use the generic division free det alog with n^4 complexity unless

• there is a special type (fmpz_mat, ...)
• I actually call a different function? (det_euc which assumes the ring is euclidean without checking anything)

I do not have a good, nice solution. This would be "traits". At this point, in all examples I currently have, they are cheap to set on creation.

I don't know if it is a good trade-off to say: bit testing and the if statement at runtime is too expensive, rather call a method with the wrong complexity - or force the user to have more knowledge. On the other hand, testing this everywhere is error-prone.

Summary: there is a problem, thus an opportunity, but no solution so far. Especially not a good one.

Any ideas? That's why I created the issue...

-- You are receiving this because you authored the thread. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/339#issuecomment-502757156

[wbhart]

Do you have a way of implementing iseuclidean_useful for the rings you have in mind?

I have no problem calling such a function and using the hnf. But I don't see how it will solve your problem (since I don't understand precisely your problem).

[wbhart]

And you'll need a gcd and divrem with canonical representatives, of course. Do you have those, even for Z[i]?

[wbhart]

By the way, I discussed with Mohamed the issue of iseuclidean. What he seems to need is an isinterface_euclidean, which will tell him if a ring supports the entire Euclidean interface. What matters to him is composability, and whether or not functionality composes depends on whether the entire interface is implemented.

So we have the following:

• iseuclidean : mathematically tries to check if it is Euclidean or not, can fail if we don't have the mathematics

• isinterface_euclidean : the entire Euclidean interface is implemented for this ring (there is a sketch of a Euclidean Interface in the docs, but no doubt it will evolve; @kalmarek 's conformance testing will be useful here to make sure we deliver what we promise)

• iseuclidean_type : the cheapest test as it is compile time, but only applies for rings like Z, K[x] (for K a field) and so on, that we can determine are Euclidean just from the type

• iseuclidean_fast : a fast runtime check for Euclidean that will not slow the system down like the mathematical check

The problem I foresee is that we need to decide which one(s) to use in library code. The mathematical one is not useful, since it can take arbitrarily long. However, there is always someone who is going to want their pet ring to behave like a Euclidean ring because they spent a week running the computation to determine that it is one.

The only way I see to work around that is to have a flag iseuclidean in each ring parent. This will get set by the iseuclidean test, and we'll have to make iseuclidean_fast return true if that flag is set. It can of course also set the flag if it itself determines that a ring is Euclidean, but this shouldn't actually make a big difference to runtime since it is designed to be fast anyway.

The next problem to resolve is how to deal with Euclidean rings vs Euclidean domains. Code can make different assumptions in the one case vs the other. Or, we can just have isdomain, with a similar list of predicates.

As for dispatch on type, theoretically we could introduce a EuclideanRing type to wrap any existing ring with a type that belongs to a Euclidean class. Then we could just write routines for Euclidean rings/domains that accept anything that belongs to this abstract type. This would be better that the current union of fields and residue rings with a crash when it isn't actually Euclidean. But it's worse in that the user would always be required to wrap their ring elements if they wanted to run Euclidean functionality.

The only place this would be helpful is when trying to get code to dispatch to other functions that it wouldn't otherwise accept, because we now know the ring is Euclidean, e.g. the functions we currently have that only accept fields or residue elements. If some other ring comes along and we determine it is euclidean, we can't call that code on it because it only accepts fields and residue rings. But if we make that code also accept anything that belongs to EuclideanDomain, we only have to wrap it at the point we do the dispatch and now the code will work on it.

For example, det(M) for general rings would first test if the ring is Euclidean. If so, there is currently no way to pass it to the Field/ResElem code. But we can do:

if iseuclidean_fast(base_ring(M))
   M_wrapped = change_base_ring(EuclideanDomain(base_ring(M)), M)
   d = det(M_wrapped)
   return d.data
end

This would be the easiest way of not duplicating oodles of functions, and or ending up with a rabbit warren of if statements everywhere which will be brittle. The EuclideanDomain can provide only the functions promised by the Euclidean interface. This way we ensure that if the code doesn't work, someone relied on functionality that isn't necessarily there.

[kalmarek]

@fieker @wbhart the julia way is to do it via traits, and (since julia falls into "no-automatic-conversion" language group) ask user to specifically compute iseuclidean, memoize it, and dispatch (=choose the optimal method) via trait; The same happens everywhere (e.g. DiffEqs):

struct EuclideanDomain end
struct NonEuclideanDomain end

det(::Type{EuclideanDomain}, M) = ... # optimized
det(::Type{NonEuclideanDomain}, M) = ... # generic

det(M) = ismemoized(iseuclidean, base_ring(M)) ? det(iseuclidean(base_ring(M)), M) : det(NonEuclideanDomain(), M)

Here ismemoized checks if the property has been actually computed. Of course dispatch will happen at run-time (and it's based on value = slow), but since it'll be contained only to dispatch the cost will be negligible compared e.g. to 100x100 determinant. In my measurings the cost of fetching such value is of order 50µs -- nothing we want do run in a hot loop, but generally acceptable. And this could be adapted to store only those values which can not be computed from type. e.g. by adding

@inline ismemoized(::Type{typeof{iseuclidean}}, ::Field) = true
@inline iseuclidean(::Field) = EuclideanDomain()

we will make the dispatch for fields a compile-time problem.

[wbhart]

Sounds interesting. Is memoisation thread safe?

I thought traits were compile time only. So how can a trait be set at runtime?

[wbhart]

Ok, I see roughly what the answer to my second question is. There's just a main function which checks the trait and dispatches accordingly.

This doesn't quite solve the same problem, though we should use that pattern for similar cases as it will be familiar to other Julia coders.

The problem I want to solve is how to make our existing functions accept types they didn't accept before. Otherwise we just end up duplicating the functions. I guess what you are suggesting is pretty much what I was suggesting, except for the memoisation thing, which I didn't know about.

[wbhart]

Can you show how I would write a single function in that paradigm that would accept one of two traits? Let's say either EuclideanRing or IntegralDomain?

[fieker]

On Mon, Jun 17, 2019 at 04:07:13PM -0700, wbhart wrote:

Can you show how I would write a single function in that paradigm that would accept one of two traits? Let's say either EuclideanRing or IntegralDomain?

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/339#issuecomment-502883170

Marek's traits are the same as "my" isuseful_euclidean, the new approach is to NOT store this information on the ring - which would require the ring to be changed to contain this, but to have a global outside storage for them. I'd say, Marek's approach is perfect for "outside" rings (where we cannot store additional info in the ring), but for ours I'd argue to make space for such traits, or bit-flags. In particular in view of serialisation....

Julia also does this, e.g. for iterators, where algorithm selection in collect is based on a number of "traits": IteratorSize IteratorEltype and others

The problem is not so much for say solve(A, b) to solve Ax = b in the "best" way, here I could call solve_euc, but for e.g. modules over (unknown) rings: Mohamed taught us how to extend linear algebra to module operations generically (hence the only possible approach for unknown situations). This depends on the linear algebra to automatically find a good algorithm. Thus traits (if thats what you want to call them).

It needs

• an interface to
  • test if known
  • test if can be computed
  • set/ compute
• to exist for useful properties (euc, field)
• for mathematical properties (useless in practise, but important to mathematicians)

• be able to be extended easily

  eg: if isknown_field(R) call field stuff if isknown_euc(R) call euc stuff if isknown_domain(R) ...

The order of operations cannot be fixed in advance, it will depend on the importance for the method. I agree with Bill that this is not good as the decision tree is complicated and prone to errors. But I don't know any alternatives. Maybe have _solve with just type based dispatch and solve with traits calling _solve appropriately?

It needs to be complemented by ways of setting those traits. Some, like isknown_field(R) might be non-trivial (isprime, isirreducible) so should not be always done automatically. However, I need to have a way "prepare_for_linalg(R)" to set those traits that apply? or isfield(R) to test, and set (cache)?

(I know its tricky: intelligent users will do R = () isfield(R) ... and get fast results if R is suitable while other will not do the isfield() and suffer in runtime )

Those traits need 2 bits at least tested or not present or not (three values: has property, does not, not known)

How do things like size(R), characteristic(R) fit into this? Or rand(R) vs. rand(R, 1:10)? hasrand(R)? hasrand_size(R)? hassmall_rand(R)?

Other: in Z[i] I know how to do the euclidean interface, similar for other small quadratics (done this in Magma)

canonical_unit(a) is not canonical nor unique in general. It is mainly a hint/ way to improve the element modulo units. we have it for fields: cu(a) = a, so xinv(cu(x)) = 1 Z cu(a) = sign(a) k[x] cu(a) = lead(a) then it gets more complicated. In Z/mZ cu(a) is defined s.th. xinv(cu(x)) = gcd(x, m) (and similar for all euclidean quotients) (this thus relies on the gcd to be normalised)

in Z_k/id: no idea (well not quite true: if I allow factoring of the ideal, I can make it unique)

localisations: for Z: OK Z_k: inhertited from Z_k/id in one case, not there in the 2nd

Result: we cannot gurarantee uniqueness - unless s.o. comes up with a good new algorithm.

[wbhart]

I am not sure I see the usefulness of isknown_euclidean. What will we immediately do if it is not known? We will call iseuclidean_fast. We may as well just always call it. The field itself can store three values as a nullable boolean.

julia> mutable struct bill
         iseuclidean::Union{Bool, Nothing}

         function bill()
            return new(nothing)
         end
       end

julia> a = bill()
bill(nothing)

julia> a.iseuclidean

julia> a.iseuclidean = true
true

julia> a.iseuclidean
true

julia> a.iseuclidean = false
false

[wbhart]

I agree we need to implement our traits per ring, not per ring type. So memoisation as it is implemented in Julia may not be workable (unless objects can be memoised).

But we can still implement the same memoisation interface for our parent objects, just as we have mimicked lots of other functionality for our parents. The question is, will that paradigm actually work for us.

We need to be able to select by type when traitA and traitB hold and we need to be able to select by type when traitA or traitB hold, and every combination thereof.

It would also be better if we had a function ismemoized say, that returns true iff the property has been computed and it is true. We need something else (ismemoized_false, say) that can determine if a property has been computed and it is false. This is a much rarer case, and I'd like to not have to check both whether something has been computed AND whether it is true in each and every case it is computed. Otherwise, getting the modus ponens and modus tollens right everywhere is a nightmare.

A third predicated ismemoised_known which determines just whether the property has been computed could also be useful in very rare cases.

[wbhart]

Where do I read about Julia memoization. I can only find this:

https://github.com/JuliaCollections/Memoize.jl

which just adds a cache to a function, which is not very flexible.

[thofma]

There is no Julia memoization.

[kalmarek]

Sounds interesting. Is memoisation thread safe?

i wouldn't bother at the moment; I assume it's not (there are no thread-safe data structures in julia, at the moment), and I'd just wait for partr to stabilize and get wider acceptance/use. In the worst case (i.e. in one-two years) You can roll your own solution with collision-free dictionary or else.

The problem I want to solve is how to make our existing functions accept types they didn't accept before. Otherwise we just end up duplicating the functions.

Didn't understand this part.

Can you show how I would write a single function in that paradigm that would accept one of two traits? Let's say either EuclideanRing or IntegralDomain?

use Unions: ::Type{<:Union{EuclideanRing, IntegralDomain}}; you can also add an abstract type DomainTypeTrait and inherit

@fieker doing id predicate do something; elseif different_predicate something else; ... is not really extendable way of doing this. If I implement my hypothetical ring of DiagonalMatrices i can't really hook into if... elseif... "dispatch. However I can implement

ismemoized(..., ::DiagonalMatrix) = true;
iseuclidean(..., ::DiagonalMatrix) = DIAGONAL();
det(::Type{DIAGONAL), M) = # super optimized det for diagonal matrices

to have it going super fast. We may think if iseuclidean predicate is the best one for choosing the algorithm for det (no it is not), but you should get the idea; trait based dispatch is easily extendable to user constructs.

BitVectors could be used to implement bitflags, but we'd need two of them to encode "not known", as @fieker wrote.

(I know its tricky: intelligent users will do R = () isfield(R) ... and get fast results if R is suitable while other will not do the isfield() and suffer in runtime )

we should not guard users against their lack of knowledge -- if You know more about the domain, share that knowledge with us. If You don't -- suffer the complexity. Or at least I don't think we should do it at this stage/level. What I believe is we're trying to get an easily extendable low-level solution here.

@wbhart memoize caches by (type of function, args), I'm not sure why it is not flexible?

[kalmarek]

There is no Julia memoization.

i.e. it's not in Base or stdlib

[fieker]

On Tue, Jun 18, 2019 at 01:50:20AM -0700, wbhart wrote:

I am not sure I see the usefulness of isknown_euclidean. What will we immediately do if it is not known? We will call iseuclidean_fast. We may as well just always call it. The field itself can store three values as a nullable boolean. For euclidean there is no point. For is field there is: Z/nZ is a field iff n is prime. once I checked this I need to store this information. I now know if I have a domain (field) or not, I never need to test again, regardless of the result This is about the interface, not the names directly... I still feel there needs to be the slow, complete version...

julia> mutable struct bill
         iseuclidean::Union{Bool, Nothing}

         function bill()
            return new(nothing)
         end
       end

julia> a = bill()
bill(nothing)

julia> a.iseuclidean

julia> a.iseuclidean = true
true

julia> a.iseuclidean
true

julia> a.iseuclidean = false
false

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/339#issuecomment-503010881

[fieker]

On Tue, Jun 18, 2019 at 03:23:00AM -0700, kalmarek wrote:

Sounds interesting. Is memoisation thread safe?

i wouldn't bother at the moment; I assume it's not (there are no thread-safe data structures in julia, at the moment), and I'd just wait for partr to stabilize and get wider acceptance/use. In the worst case (i.e. in one-two years) You can roll your own solution with collision-free dictionary or else.

The problem I want to solve is how to make our existing functions accept types they didn't accept before. Otherwise we just end up duplicating the functions.

Didn't understand this part.

Can you show how I would write a single function in that paradigm that would accept one of two traits? Let's say either EuclideanRing or IntegralDomain?

use Unions: ::Type{<:Union{EuclideanRing, IntegralDomain}}; you can also add an abstract type DomainTypeTrait and inherit

@fieker doing id predicate do something; elseif different_predicate something else; ... is not really extendable way of doing this. If I implement my hypothetical ring of DiagonalMatrices i can't really hook into if... elseif... "dispatch. However I can implement

ismemoized(..., ::DiagonalMatrix) = true;
iseuclidean(..., ::DiagonalMatrix) = DIAGONAL();
det(::Type{DIAGONAL), M) = # super optimized det for diagonal matrices

I don't get this. If I want to call a special function for my matices, I do det_diag(M). What I want to do is det(M) always and this needs to select the "proper" det function. So far, we replaced det_for_euc by det(M, Type{Euc})

function det(M::Matrix) isdiag_type(typeof(M)) && return det(M, Type{Diag})

R = base_ring(M)

iseuc(R) && return det(M, Type{Euc}) isfield(R) && return det(M, Type{Field}) isdomain(R) && return det(M, Type{Domain})

isgreen(R) && iscommutative(R) && return det(M, Type{Strange}) (isgreen(R) || iscommutative(R)) && return det(M, Type{Strange})

return det_df(M) end

This looks like it should be

function trait(R::Ring) iseuc(R) return Type{Euc} isfield(R) return Type{Field} end

det(M) = det(M, trait(base_ring(M)))

or s.th. similar...

to have it going super fast. We may think if iseuclidean predicate is the best one for choosing the algorithm for det (no it is not), but you should get the idea; trait based dispatch is easily extendable to user constructs.

BitVectors could be used to implement bitflags, but we'd need two of them to encode "not known", as @fieker wrote.

(I know its tricky: intelligent users will do R = () isfield(R) ... and get fast results if R is suitable while other will not do the isfield() and suffer in runtime )

we should not guard users against their lack of knowledge -- if You know more about the domain, share that knowledge with us. If You don't -- suffer the complexity. Or at least I don't think we should do it at this stage/level. What I believe is we're trying to get an easily extendable low-level solution here.

@wbhart memoize caches by (type of function, args), I'm not sure why it is not flexible?

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/339#issuecomment-503044306

[wbhart]

@fieker and I had a long discussion about this, and here is what we came up with.

• We want a set of mathematical predicates, iseuclidean(R), isdomain(R), etc., which store a flag in the ring/parent when the predicate is computed. This function should take an optional argument fast which should default to false. Thus, when the user asks, iseuclidean(R), it can take as long as it wants to, but will return either true/false/maybe. When fast is set to true it takes at most one millisecond. We would never call it with fast=false in generic code.

• We need a function hastrait(R, iseuclidean) which first checks (compile time only, depends only on types) whether the Euclidean interface is implemented, then calls the above function with fast = true to determine it has the mathematical property of being Euclidean. This function must return only true or false, no other value.

• We only want to dispatch on the return value of hastrait, (never purely on the value of iseuclidean, for which we would use a simple if statement), which we can do as follows:

function myfun(M::Matrix)
   R = base_ring(M)

   if hastrait(R, iseuclidean, isdomain)
      return myfun(M, Trait(Euclidean, Domain))
   elseif hastrait(R, isfield)
      return myfun(M, Trait(Field)) 
   end

   # generic fallback
end                   

The Trait function will return a boolean tuple type in a fixed order with the Euclidean and Domain fields set to true, meaning both hold.

In the first case above, the corresponding myfun function would look like:

function myfun(M::Matrix, ::Trait(Euclidean, Domain))

The Trait function here would return the appropriate type that needs to match.

Of course we also need to set it up so we can specify && and || for all traits, but I vaguely think this is possible with Unions and various combinations of trait fields in the tuples set to true.

I don't think we can compress the two lines of each if statement into one, since we want to both test if the traits are satisfied and to dispatch to whatever function we need to based on that, and still go on to other tests or fallback code if the traits are not satisfied.

Do you see any potential problems with this @kalmarek ?

[wbhart]

And refresh your browser. I edited it.

[kalmarek]

@fieker

function det(M::Matrix) 
  isdiag_type(typeof(M)) && return det(M, Type{Diag}) 
  R = base_ring(M)
  iseuc(R) && return det(M, Type{Euc}) 
  isfield(R) && return det(M, Type{Field}) 
  isdomain(R) && return det(M, Type{Domain})
  isgreen(R) && iscommutative(R) && return det(M, Type{Strange}) 
  (isgreen(R) || iscommutative(R)) && return det(M, Type{Strange})
  return det_df(M)
end

the problem with this code is that it will never run my specialized version of det; and we definitely don't want to see det_diag -- this is very much against multiple dispatch;

@wbhart fast=true arg sounds very sensible; however if based myfun again is not user extendable. if we have abstract type DomainAbstractTrait, users should be able subtype it with MyStrangeTrait <:DomainAbstractTrait and to extend myfunc(M, MyStrangeTrait) to get things running; fallback should not be done on the base of if ... else, but rather myfunc(M, ::Type{<:DomainAbstractTrait}).

This way we could define produce_traits(R) = (DomainAbstractTrait,) function and call myfunc(M) = myfunc(M, prodce_traits(base_ring(M))...).

@wbhart, for standard traits we need to agree on the order of traits, I agree. and packaging of traits is also something to think about (I packaged is as a tuple above, but I didn't think seriously about it).

[kalmarek]

and maybe there should be a place for StandardTraits and UserDefinedTraits?

I'd say it's impossible to get this right on first try, before we try actual code ;)

[fieker]

On Tue, Jun 18, 2019 at 06:47:06AM -0700, kalmarek wrote:

@fieker

function det(M::Matrix) 
isdiag_type(typeof(M)) && return det(M, Type{Diag}) 
R = base_ring(M)
iseuc(R) && return det(M, Type{Euc}) 
isfield(R) && return det(M, Type{Field}) 
isdomain(R) && return det(M, Type{Domain})
isgreen(R) && iscommutative(R) && return det(M, Type{Strange}) 
(isgreen(R) || iscommutative(R)) && return det(M, Type{Strange})
return det_df(M)
end

the problem with this code is that it will never run my specialized version of det; and we definitely don't want to see det_diag -- this is very much against multiple dispatch;

you would write function det(M::diag_type) ... end

this is dispatch on traits of the ring, not the actual object. The object traits are a completely different problem...

@wbhart fast=true arg sounds very sensible; however if based myfun again is not user extendable. if we have abstract type DomainAbstractTrait, users should be able subtype it with MyStrangeTrait <:DomainAbstractTrait and to extend myfunc(M, MyStrangeTrait) to get things running; fallback should not be done on the base of if ... else, but rather myfunc(M, ::Type{<:DomainAbstractTrait}).

This way we could define produce_traits(R) = (DomainAbstractTrait,) function and call myfunc(M) = myfunc(M, prodce_traits(base_ring(M)...).

@wbhart, for standard traits we need to agree on the order of traits, I agree. and packaging of traits is also something to think about (I packaged is as a tuple above, but I didn't think seriously about it).

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/339#issuecomment-503142722

[kalmarek]

@fieker If the actual method will be always deducible from type we woudn't be discussing this ;-)

[wbhart]

@kalmarek I don't have any problem with the fallback function dispatching on some general trait.

I am not entirely sure if you are right about not using if statements too.

I think we can use the covariance of tuples in Julia to handle the trait functions. The functions will accept (unions of) Tuples of Bool and Type{true}, I think. They will be passed Tuples of Type{true} and Type{false}.

(@fieker BTW, this doesn't handle tuples of booleans that are shorter than the full list of traits. But this is only a problem if people hand code it rather than use the Trait functions designed to get the order and length right for them.)

@kalmarek How would I determine which traits to check without if statements? You certainly don't want it to just compute all known traits for every call. You only want it to check the ones you specify, and you may only want some traits to be checked if others failed, etc. I don't see how you manage that without if statements in the main function.

[wbhart]

By the way, I certainly don't care about user defined traits at this point. We can put that in later.

Typically the user will want to define their own trait system for their own module, i.e. this whole thing assumes there will be a linear algebra traits system, a ring trait system, a module trait system and so on.

I don't think having a list of 300 traits to every function is workable.

[wbhart]

The main problem I have with implementing this is Tuples are the only covariant type in Julia and we need a function that takes a tuple of Type{true}'s and Type{false}'s and does a "logical and" on them. We'll have to give that a funky name so we don't get into type piracy.

[wbhart]

Actually what I said about Tuples being the only covariant type in Julia is not true. Julia now supports making type paramters covariant, contravariant or invariant, for all types. So we can define a custom Trait abstract type if we want.

[fieker]

On Tue, Jun 18, 2019 at 07:07:50AM -0700, kalmarek wrote:

@fieker If the actual method will be always deducible from type we woudn't be discussing this ;-)

It certainly will not always be the case, I agree. But for this proposal, its about ring traits, not matrix traits...

I agree, a more general model might be useful at times, however, I have a feeling we're quickly going down to over-engineering for all eventualities. It would be cool to have an extensible fast mechanism, but this is probably impossible. Gap does this with their predicate dispatch/ method selection Polymake even more general with some optimisation algorithm based on graphs (if my memory of Michael's explanation is correct)

Neither is fast. neither is (for beginners/ users) predictable. nor extensible nor perfect. Both operate under different assumptions: Gap is interpreted and thus can never be that fast. It is intended for longer computations to make informed decisions at the beginning. Polymake is even worse: this seems to select a strategy for the next couple of days/ weeks/ months. Thus even an hour for the selection is fine (please correct me on anything)

Of course, nothing can improve on 2x2 matrices over small ints.

We're trying to come up with s.th. feasible that solves our (current and forseeable future) problems. So far, this is about algorithms selected by ring properties. It is also not about mathematical properties, but actual useful ones, thus in computer algebra, in the context of linear algebra, a fairly small set of properties exact/ inexact domain field euc I actually don't see many more generic ones... maybe very long term Dedekind, but that's a large interface. Also remember, this is on top of type based selection, thus matrices over fields do not need to check the field trait and can directly do s.th. special.

Possibly Mohamed would add localisation (as he can reduce lin. alg to the ring) but this is not really orthogonal for our types: localisation of euc is euc, the only other case is some localisations of Dedeking and thats very special (so far).

Matrix traits/ types is different, there I can see even user defined types: banded matrices, banded plus some other, nilpotent, diagonal, triagonal, square, invertable ... Ideally some of those have different types, others, in the context of det, solve and nullspace should be trivial to handle by the ring based selection (rref trivial by obvious pivoting). But maybe, matrices should have properties/ traits/ more information.

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/339#issuecomment-503151773

[wbhart]

@kalmarek I agree it would be nice if users could subtype our AbstractTrait type, but this just doesn't support functions which must have all of a certain collection of traits (plural) or all of another collection of traits (plural).

Think (isdomain && iseuclidean) || isartinian

How would you do that with subtyping? I really don't think we need something that complex. @fieker is right. Let's solve the problem we need urgently solved, rather than overengineer something for every eventuality. We certainly need && and || in the functions that are being dispatched to so that we can keep the traits as orthogonal as possible. Otherwise we will have 2^n functions for everything in the limit.

[kalmarek]

@wbhart You know more about covariance of tuples, I trust you on this;

If computing traits is cheap (adds a dozen or two of miliseconds), I'm fine computing all of them. We don't have to add all of them to each function, args... would do., but obviously passing only those that matter is much more elegant... hmm...

If we don't allow trait system to be extendable to unknown traits then it becomes much easier, at the cost of being not possible to add them later. I wouldn't discard user defined types and traits (after all this is why we define interfaces, right?) but that's my personal choice. and here it's a design decision we have to make.

But maybe @fieker is right that part of the problem could be fixed through myfunc(::MyType, some, standard, traits...). At the moment I think it'd be best to actually start prototyping to see what kind of problems we can encounter in the long-run.

[wbhart]

@kalmarek I don't think we can afford to reimplement Gap here. We need something simple the compiler can obviously cope with.

Computing all of 100 traits for a few milliseconds each is obviously a bad idea. It has to work on specific traits we specify.

Moreover, I still don't see how your proposal would implement && and || on the functions being dispatched to. I can use Union{bunch of traits} to get ||, but I can't do && this way, and that is essential to limit the number of traits to 100 instead of 2^100.

I think we have a workable model. We should implement it and see how far we can push it. Of course we won't get it right, but lets hope we have many developers who can help with lifting it to a better solution once we have a better idea.

[wbhart]

Ok, we can support variable length Tuples after all, using Vararg tuples. So that will allow for future extensibility if we use it right in the implementation.

And no need to use Type{true}, we can just define abstract types for all our traits, which can in fact be in trees.

abstract type TraitAbstract end

abstract type TraitEuclidean <: TraitAbstract end

abstract type TraitField <: TraitEuclidean end

abstract type TraitDomain <: TraitAbstract end

abstract type TraitSpecial <: TraitAbstract end

function myfun(M::Int, ::Type{<:Tuple{TraitEuclidean, TraitDomain, Vararg{TraitAbstract}}})
   return 1
end

function myfun(M::Int, ::Type{<:Tuple{Vararg{TraitAbstract}}})
   return 2
end

function myfun(M::Int)
   if M == 1
      return myfun(M, Tuple{TraitEuclidean, TraitDomain})
   elseif M < 100
      return myfun(M, Tuple{})
    else
      return myfun(M, Tuple{TraitEuclidean, TraitDomain, TraitSpecial})
   end
end

julia> myfun(1)
1

julia> myfun(5)
2

julia> myfun(101)
1

[kalmarek]

@wbhart I wasn't thinking about 100s of them ;-)

• for a few traits this could work:

  
  abstract type AbstractTrait end
  const AT = AbstractTrait;
  struct unknown <: AbstractTrait end
  struct euclidean <:AbstractTrait end
  struct domain <:AbstractTrait end

myfunc(x, ::Tuple{domain, <:AT <:AT}) = 3x myfunc(x, ::Tuple{domain, <:AT, euclidean}) = 2x myfunc(2, (domain(), unknown(), euclidean())) # returns 4 myfunc(2, (domain(), unknown(), unknown())) # returns 6

* for a dozen or so we could define
```julia
struct Traits{A, B, C, ...} where A<:AT, B<:AT, C<:AT
...
end

parametrized by bunch of these and define
myfunc(x, t::Traits) = myfunc(x, (Traits[domain()], Traits[euclidean()])...)
myfunc(x, ::domain, ::euclidean) = 2x
myfunc(x, ::domain, ::unknown) = 3x

(indexing to Traits can be resolved at compile time).

of course for more than a dozen there is not point of computing all of the traits, we agree on this

EDIT: My first proposal is essentially yours without Vararg

[wbhart]

Your first solution is a bit simpler than mine; it doesn't have the ability to inherit from traits, and you are actually instantiating the types, whereas I do not, so everything remains compile time in my case.

Your code in the second case won't work as written, even if you overload getindex. You have to have a fixed number of traits if you do that. You can't put "..." in the type.

It's also brittle in that if you want to insert a trait in the dispatched-to function's type signature, you need to add a new field, then modify a bunch of other code to get those right too. I'm not so keen on that. This will lead to errors that only show up at runtime, which I'm keen to avoid where possible.

My approach is possibly also extensible by the user.

Also, why overload [] and not just write a function Traits(...) from the get-go. I would suggest that the function also sort out the order for the developer, rather than require them to specify Traits[...] over and over.

Also, you will end up with massively long function signatures in your second suggestion with any number of unknowns for traits that don't matter. Imagine dozens of them.

My approach is heavier on the compiler, but doesn't suffer from these problems, I think.

[wbhart]

Oh, actually you can do

julia> struct cas{T, S, Vararg}
       end

julia> a = cas{Int, Bool, Int}()
cas{Int64,Bool,Int64}()

which I was not aware of.

[kalmarek]

@wbhart the first one was me being lazy and not writing ::Type{<:} everywhere ;)

It's also brittle in that if you want to insert a trait in the dispatched-to function's type signature, you need to add a new field, then modify a bunch of other code to get those right too. I'm not so keen on that. This will lead to errors that only show up at runtime, which I'm keen to avoid where possible.

I didn't understand that. There are no fields (needed) in Traits, it could be an empty struct, the whole "predicate" information is in the type. Any changes to Traits will be opaque to the rest of the code, as long as getindex returns what it should. And overloading getindex will allow You to add some logic there (should we need it). Moreover the ordering of traits (e.g. in tuple approach) is contained in the Traits "class" and does not spread to the rest of the code. I'd say it's a pretty good functional separation, but You may have different opinion.

Finally it is very easy to test ;) both Traits indexing and myfun (without referring to Traits). Not to mention that it produces very clean signatures (without unused arguments and so on).

anyway, i guess it's high time we produced list of functions that may use those traits, see how the corresponding methods should look like and see what we actually need.

and it's time to go home:)

[wbhart]

I was talking about your second example

myfunc(x, ::domain, ::euclidean) = 2x
myfunc(x, ::domain, ::unknown) = 3x

These are brittle. If I add new traits, I'll need to modify these. The order of traits matters too because of how dispatch works, and it is not possible to always get 100% orthogonal traits, so we may end up having to insert traits in between existing ones.

But again I ask, what happens if there are dozens of traits. Can you imagine what these functions are going to look like!?

[kalmarek]

If we add traits... at the end of myfunc signatures, then no, we don't need to modify those; If we make changes to Traits (by adding a new trait), all the changes will be contained in Traits (namely its indexing). This is what I meant by functional separation. All the rest is done by myfunc(x, t::Traits) = myfunc(x, (selection of traits here)). And of course if You add new algorithm to myfunc which is faster in the presence of ::Green trait, then well yes, you have to write

• getindex(t::Traits, ::Type{Green}) (all of what I wrote with empty type instances could/should be done with Types).
• the actual implementation of the new algorithm with the signature: myfunc(x, ::Type{<:AT}, ::Type{<:AT}, ::Type{Green}, traits...).

Note also, that traits lists at the end will be separated as well by the type of x, as @fieker suggested. Hence I don't think there will actually be methods that would benefit from passing dozens of traits. If there are, please name a few.

[wbhart]

One possible solution to the IsEuclidean issue was something I started thinking about at our Saarbruecken meeting.

We could drop all abstract types from our concrete types entirely. Then we write a single top level function (let's take det as an example), which figures out the properties of the object(s) passed to it, i.e. whether it is a Ring, Field, Domain, EuclideanDomain, etc., then wraps the objects with a new type, e.g. Mat{EuclideanDomain}, and then just write our det functions for whichever kinds of inputs the algorithm expects, e.g. det(M::Mat{EuclideanDomain}).

This would be quite some work, but it would cleanly separate the logic of determining properties from the logic for "dispatch on algorithm", all the while preserving performance (all of the wrapper types would be immutable, for performance reasons).

This would be a lot of work to implement, but we could try it out in a branch on a small set of functionality and see if it works.

It would also make it a lot easier to check that our functions don't require functionality not provided by a given interface, but that is just an added bonus of this approach.

I've no idea if this is a good idea or not, but I do have the feeling that we are going to have to move to something like this to keep Oscar moving forwards.

The major disadvantage is that it makes @which much less useful. But so long as there is only ever one top level function where all the logic for deciding properties of objects is contained, we won't end up with a sprawling mess like certain other unmentionable CAS's.

[wbhart]

In particular, it will still be possible by looking at the type signature of a function to determine for what kinds of inputs it is designed to work.

[fieker]

What is the difference between this and just calling det_euc directly? It would be less cool, but easier to understand...

On Wed, Sep 25, 2019 at 06:43:24AM -0700, wbhart wrote:

One possible solution to the IsEuclidean issue was something I started thinking about at our Saarbruecken meeting.

We could drop all abstract types from our concrete types entirely. Then we write a single top level function (let's take det as an example), which figures out the properties of the object(s) passed to it, i.e. whether it is a Ring, Field, Domain, EuclideanDomain, etc., then wraps the objects with a new type, e.g. Mat{EuclideanDomain}, and then just write our det functions for whichever kinds of inputs the algorithm expects, e.g. det(M::Mat{EuclideanDomain}).

This would be quite some work, but it would cleanly separate the logic of determining properties from the logic for "dispatch on algorithm", all the while preserving performance (all of the wrapper types would be immutable, for performance reasons).

This would be a lot of work to implement, but we could try it out in a branch on a small set of functionality and see if it works.

It would also make it a lot easier to check that our functions don't require functionality not provided by a given interface, but that is just an added bonus of this approach.

I've no idea if this is a good idea or not, but I do have the feeling that we are going to have to move to something like this to keep Oscar moving forwards.

The major disadvantage is that it makes @which much less useful. But so long as there is only ever one top level function where all the logic for deciding properties of objects is contained, we won't end up with a sprawling mess like certain other unmentionable CAS's.

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/339#issuecomment-535029534

[wbhart]

Types are more powerful than names. Consider unions, default args, etc.

Moreover, you probably want to name the implementation after the algorithm, if anything, not the kinds of objects they should operate on, which is why it's nice to have that documented in the types.

It's also way less typing when actually calling the function (if you ever have the need to call it directly), though certainly it adds the extra overhead of having to wrap all the objects before you do so.

[wbhart]

Another reason: types are enforced by a type checker. Names are not. Type checkers do not, however, care about algorithms. They only care about what kind of data you have.

[fieker]

On Wed, Sep 25, 2019 at 07:06:34AM -0700, wbhart wrote:

Types are more powerful than names. Consider unions, default args, etc.

Moreover, you probably want to name the implementation after the algorithm, if anything, not the kinds of objects they should operate on, which is why it's nice to have that documented in the types. Why? I'd as happily call them by type then by algo - in fact prior to using Flint/Nemo I've never called a function by algo name...

It's also way less typing when actually calling the function (if you ever have the need to call it directly), though certainly it adds the extra overhead of having to wrap all the objects before you do so. I don't know det(change_type_to_euclidean(M)) vs det_euclidean(M)

and there is the issue of readability... yet another dispatched call does not help to follow and predict the code path. Ifs are easy to read, types are harder. Mind you, it all getting used to it. I'm not sure if this is an attempt to use the types for everything, just because we can, or if they are actually helpful here.

Basically: we have to implement the decisino tree and either at the end of it call a function with the type in the name or change the type. I guess functionally and work-wise equivalent.

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/339#issuecomment-535040234

[wbhart]

Alright, lets not do this then. Types are more flexible and I can call all functions that do the same thing by the same name if I use types. I can also make one argument accept one type and another argument accept another type, etc.

But I'm not going to bother with trying to explain why types exist. Let's just drop the idea. Obviously not worthwhile.

But it still means we have no really workable solution for a Euclidean interface.

[fieker]

On Wed, Sep 25, 2019 at 09:33:17AM -0700, wbhart wrote:

Alright, lets not do this then. Types are more flexible and I can call all functions that do the same thing by the same name if I use types. Yep. I can also make one argument accept one type and another argument accept another type, etc. Yep.

But I'm not going to bother with trying to explain why types exist. Let's just drop the idea. Obviously not worthwhile. Types exists for a reason, Untyped stuff exists for a reason. Neither solve all problems. I'm happy to discuss and experiment. Maybe we should try this for 2 functions either way and see what looks nicer, is easier to use and easier to teach ...

But it still means we have no really workable solution for a Euclidean interface. yep...

-- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: https://github.com/Nemocas/AbstractAlgebra.jl/issues/339#issuecomment-535105264