# The GAP code in this function has been taken, essentially verbatim, from 
# an appendix to the paper
# Hulpke, Stanovsky, Vojtechovsky: Connected quandles and transitive groups

############################################################
############################################################
InstallGlobalFunction(QuandleIsomorphismRepresentatives,
function(L)
local C,R, P, Q, invariant, ExtendHomomorphismByClosingSource,
ExtendIsomorphism, IsomorphismQuandles, QuandlesUpToIsomorphism,
QuandleGenerators, T, Subquandle;

#########################
invariant:=function(QQ)
local rt, A, F, Q;
Q:=MagmaByMultiplicationTable(QQ);
rt:=RightMultiplicationGroupOfQuandle(Q);
if HasIdGroup(rt) then
rt:=IdGroup(rt);
else
rt:=false;
fi;
A:=[rt];
return A;
end;
#########################

#########################
Subquandle := function( Q, gens )
# returns the subquandle of Q generated by gens as a list of elements 
    local sub, i, j;
    sub := ShallowCopy( gens );
    repeat 
        gens := ShallowCopy( sub );
        for i in gens do for j in gens do
            AddSet( sub, Q[i][j] );
        od; od;
    until Length( sub ) = Length( gens );
    return sub;
end;
#########################

#########################
QuandleGenerators := function( Q )
# returns a small generating set of Q
    local gens, diff, n;
    gens := [ ];
    n := Length( Q );
    diff := [1..n];
    while diff <> [ ] do
        AddSet( gens, diff[ 1 ] );
        diff := Difference( diff, Subquandle( Q, gens ) );
    od;
    return Set( gens );
end;
#########################


ExtendHomomorphismByClosingSource := function( f, L, M )
# <L>, <M> are multiplication tables, <f> is a partial map from a subset of elements of <L>
# to a subset of elements of <M>.  This function attempts to extend <f> into a homomorphism
# by extending the source of <f> into (the smallest possible) subquandle of <L>.

    local oldS, newS, pairs, x, y, newNow, p, z, fz;    
    oldS := [ ];
    newS := f[ 2 ];

    repeat  
        pairs := [];
        for x in oldS do for y in newS do 
            Add( pairs, [ x, y ] ); 
            Add( pairs, [ y, x ] );
        od; od;
        for x in newS do for y in newS do
            Add( pairs, [ x, y ] );
        od; od;
        newNow := [];
        for p in pairs do
            x := p[ 1 ];
            y := p[ 2 ];
            z := L[ x ][ y ];
            fz := M[ f[ 1 ][ x ] ][ f[ 1 ][ y ] ];
            if f[ 1 ][ z ] = 0 then
                f[ 1 ][ z ] := fz; AddSet( f[ 2 ], z ); AddSet( f[ 3 ], fz );
                Add( newNow, z );
            else 
                if not f[ 1 ][ z ] = fz then return fail; fi;
            fi;
        od;
        oldS := Union( oldS, newS );
        newS := ShallowCopy( newNow );
    until IsEmpty( newS );
    return f;           
end;

ExtendIsomorphism := function( f, L, GenL, M )
# Given a partial map <f> from quandle <L> to quandle <M>, it attempts to extend
# <f> into an isomorphism betweem <L> and <M>.
# <GenL> is a small generating set of <L>.
    local x, possible_images, y, g, n;
    n := Length( L );
    f := ExtendHomomorphismByClosingSource( f, L, M );
    if f = fail or Length( f[ 2 ] ) > Length( f[ 3 ] ) then return fail; fi;
    if Length( f[ 2 ] ) = Length( L ) then return f; fi; #isomorphism found
    
    x := GenL[ 1 ];
    GenL := List( [ 2..Length( GenL ) ], i -> GenL[ i ] ); 
    possible_images := Filtered( [1..n], y -> not y in f[ 3 ] );    
    for y in possible_images do
        g := StructuralCopy( f );
        g[ 1 ][ x ] := y; AddSet( g[ 2 ], x ); AddSet( g[ 3 ], y );
        g := ExtendIsomorphism( g, L, GenL, M );
        if not g = fail then return g; fi; #isomorphism found
    od;
    return fail;    
end;

IsomorphismQuandles := function( L, M )
# returns an isomorphism of L to M if such exists, else returns fail.
    local GenL, map, iso;
    GenL := QuandleGenerators( L );
    map := 0 * [ 1.. Length( L ) ]; map[ 1 ] := 1; # THIS IS ONLY OK IN INDECOMPOSABLE QUANDLES
    iso := ExtendIsomorphism( [ map, [ 1 ], [ 1 ] ], L, GenL, M);
    if not iso = fail then return SortingPerm( iso[ 1 ] ); fi;
    return fail;
end;

QuandlesUpToIsomorphism := function( ls )
# given a list of quandles ls, returns representatives up to isomorphism
    local qs, L, D, G, with_same_D, is_new, K;
    qs := [];
    for L in ls do
        D := invariant( L );
        G := QuandleGenerators( L );
        # will be testing only quandles with the same discriminator
        with_same_D := Filtered( qs, K -> K[2] = D );
        is_new := true;
        for K in with_same_D do # K is a set of ordered tuples [quandle, discriminator]
            if not IsomorphismQuandles( L, K[1] ) = fail then
                is_new := false;
                break;
            fi;
        od;
        if is_new then Add( qs, [ L, D ] ); fi;
    od;
    # returning only quandles, not their discriminators
    return List( qs, L -> L[1] );
end;

T:=List(L,x->MultiplicationTable(Elements(x)));

return List(QuandlesUpToIsomorphism(T), x->MagmaByMultiplicationTable(x));

end);
############################################################
############################################################