from __future__ import print_function
import sys
from itertools import chain
from random import random

class PoissonLieGroup(SageObject):

    def __init__(self, cartan_type, gamma=dict(), gammaprime=None, **kwargs):
        self.cartan_type = lambda : CartanType(cartan_type)
        self.rank = lambda : self.cartan_type().rank()
        self.Q = self.cartan_type().root_system().root_lattice()
        self.P = self.cartan_type().root_system().weight_space()
        self.P.print_options(latex_prefix='\\omega', prefix='omega') 
        self.roots = self.Q.positive_roots()
        
        # Store gamma
        self.gamma = gamma
        self.gamma_inv = dict( (j,i) for (i,j) in gamma.iteritems() )

        # Set up gamma as a morphism on the weight lattice. Here gamma needs to have the same indexing set of the root system
        self.gamma_morphism = self.Q.module_morphism(codomain=self.Q, on_basis=lambda i: self.Q.simple_root(gamma[i]) if i in gamma else self.Q(0))

        # Store 
        # Helper function to find roots related by <.
        def related(a):
            b = self.gamma_morphism(a)
            # Apparently < is defined only on roots of the same height
            while a.height() == b.height():
                yield (a,b)
                b = self.gamma_morphism(b)
        
        self._related_roots = list(chain(*[related(a) for a in self.roots]))

        self._simple_related_roots = [ (a,self.gamma_morphism(a)) for a in self.Q.simple_roots() if a.height() == self.gamma_morphism(a).height() ]
        
        # Type specific init
        getattr(self, '_init_type_'+self.cartan_type().type())(self.rank(), **kwargs)

        # Setup the variables, these are counted by 1!
        # TODO: for the moment we do not take into account the relations defining the group and work in SR. 
        # This is faster (?)  but not precise.
        # We might eventually want to work in an appropriate quotient ring of the polynomial ring generated by the x's.
        self.X = matrix(self._n, lambda i,j: var('x_%d_%d'%(i+1,j+1), latex_name='x_{%d,%d}'%(i+1,j+1)))
        
        # compute the cartan part of the r matrix
        A = self.cartan_type().cartan_matrix()
        scr = dict(enumerate(self.Q.simple_coroots()))
        # Half of t0
        r0 = Matrix(self.rank(), lambda i,j: self.form(self.h[scr[i]],self.h[scr[j]])).inverse()/2
        # Skew-symmetric variables
        r0 += Matrix(self.rank(), lambda i,j: sgn(j-i)*var('l_%d_%d'%tuple(sorted((i,j)))))
        # Equation given by gamma
        eqns = []
        idx = dict(map(reversed, enumerate(self.cartan_type().index_set())))
        # Do we want more equations?
        if gammaprime:
            if any(gamma[i] != gammaprime[i] for i in gamma):
                raise ValueError("gammaprime should be a map extending gamma.")
            gamma = gammaprime
        for i in gamma:
            ga_otime_one = diagonal_matrix(A.column(idx[gamma[i]]))*r0
            one_otime_a = r0*diagonal_matrix(A.column(idx[i]))
            eqns += map(lambda (x,y): sum(list(x)+list(y)), zip(ga_otime_one.columns(),one_otime_a.rows()))
        # solve
        if self.rank() > 1:
            variables = r0.variables() 
            sol = solve(eqns, variables,solution_dict=True)[0]
            for x in sol: 
                if not sol[x].is_constant():
                    sol[x] = sol[x].subs([ y == 0 for y in  sol[x].variables()])
            r0 = r0.subs(sol)
        # remember r0 just for easy comparison with old code
        self._r0_matrix = r0
        # Finally make dictionary
        self._r0 = dict(((scr[i],scr[j]),r0[i,j]) for i in scr for j in scr)

        # Compute the substitutions needed for vector fields computations
        # Cheater! we do not compute the exponential of our elements but just the first two terms
        t = var('t')
        idxs = list(chain(*[self.e.iteritems(),self.f.iteritems(),self.h.iteritems()]))
        self._left_vf_subs = dict( (r,dict(zip(self.X.list(),(self.X*(1+t*M)).list()))) for (r,M) in idxs )
        self._right_vf_subs = dict( (r,dict(zip(self.X.list(),((1+t*M)*self.X).list()))) for (r,M) in idxs )

        # setup Lie algebra
        self.lie = LieAlgebra(SR, cartan_type=cartan_type)
        self.lie.print_options(prefix='X', latex_prefix='\\xi')

        # Lie algebra version of gamma
        g = matrix(SR, self.rank(), lambda i,j: SR.var('t_%d_%d'%(i,j)))
        I = identity_matrix(self.rank())
        for k in self.gamma:
            g.set_column(k-1,I[self.gamma[k]-1])
        # the following is the correct thing to do but it is really expensive
        # the workaround gives the same matrix up to a scalar
        #K = matrix(self.rank(), lambda i,j: self.lie.killing_form( self.lie.basis()[self.Q.simple_coroot(i+1)], self.lie.basis()[self.Q.simple_coroot(j+1)] ) )
        K = self.cartan_type().cartan_matrix()/diagonal_matrix(self.cartan_type().symmetrizer())
        V = K.matrix_from_columns([i-1 for i in self.gamma]).kernel()
        sol = solve(flatten([list(g*v) for v in V.basis()]), g.variables())[0]
        #sol = solve(union(*[g*v for v in V.basis()]), g.variables())[0]
        g = g.subs(sol)
        @cached_function
        def lie_gamma(key):
            if key in self.Q:
                if not all(i in self.gamma.keys() for i in key.support()):
                    return self.lie(0)
                else:
                    return self.lie.basis()[key.map_support(lambda k: self.gamma[k])]
            else:
                return sum( c*self.lie.basis()[self.Q.simple_coroot(i)] for (i,c) in  enumerate(g*key.to_vector(),1) )
        self.lie.gamma = self.lie.module_morphism(on_basis=lie_gamma, codomain=self.lie)

        # Lie algebra version of gamma_star
        gs = K.inverse()*g.transpose()*K
        @cached_function
        def lie_gamma_star(key):
            if key in self.Q:
                if not all(i in self.gamma_inv.keys() for i in key.support()):
                    return self.lie(0)
                else:
                    return self.lie.basis()[key.map_support(lambda k: self.gamma_inv[k])]
            else:
                return sum( c*self.lie.basis()[self.Q.simple_coroot(i)] for (i,c) in  enumerate(gs*key.to_vector(),1) )
        self.lie.gamma_star = self.lie.module_morphism(on_basis=lie_gamma_star, codomain=self.lie)

        # Lie algebra triangular decomposition
        self.lie.pi_plus = self.lie.module_morphism(on_basis=lambda a: (1 if a in self.e.keys() else 0)*self.lie.basis()[a], codomain=self.lie)
        self.lie.pi_minus = self.lie.module_morphism(on_basis=lambda a: (1 if a in self.f.keys() else 0)*self.lie.basis()[a], codomain=self.lie)
        self.lie.pi_zero = self.lie.module_morphism(on_basis=lambda a: (1 if a in self.h.keys() else 0)*self.lie.basis()[a], codomain=self.lie)

        # Lie algebra R matrix
        lg = self.lie.gamma.matrix()
        lgs = self.lie.gamma_star.matrix()
        self.lie.R = self.lie.module_morphism(matrix= (1+lg)/(1-lg)*self.lie.pi_plus.matrix() - (1+lgs)/(1-lgs)*self.lie.pi_minus.matrix() + (1/(1-lg)-1/(1-lgs))*self.lie.pi_zero.matrix(), codomain=self.lie)

        # setup numerical evaluation
        self.reset_numerical()

        # setup BFZ word
        self.set_word()

        # matrix of form in fundamental weight basis
        self._form_on_weights = Matrix(self.rank(), lambda i,j: self.form(self.h[scr[i]],self.h[scr[j]])).inverse()

        # r0 asymmetry
        self._r0_asymmetry = self._r0_matrix - self._form_on_weights/2

        # setup the torus grading
        preserved = dict(enumerate([i for i in range(1,self.rank()+1) if i not in self.gamma.values()]))
        self.T = CombinatorialFreeModule(QQ, basis_keys=preserved.values(), prefix='rho', latex_prefix='\\rho' )
        pi = matrix(len(preserved), self.rank(), lambda i,j: var('p_%d_%d'%(i,j)))
        IP = identity_matrix(SR, self.rank())
        IT = identity_matrix(SR, len(preserved))
        eqns = flatten([(pi*IP.column(preserved[i]-1) - IT.column(i)).list() for i in preserved])
        A = self.cartan_type().cartan_matrix()
        eqns += flatten([ (pi*(A.column(i-1)-A.column(j-1))).list() for (i,j) in self.gamma.iteritems() ])
        sol = solve(eqns, pi.variables())[0]
        self._PtoT = self.P.module_morphism(matrix=pi.subs(sol) , codomain=self.T)

    def _init_type_A(self, rk, **kwargs):
        self._n = rk+1
        self.h = dict((self.Q.simple_coroot(i+1),self._e(i,i) - self._e(i+1,i+1)) for i in range(rk))
        self.e = dict((r,self._e(min(dict(r))-1,max(dict(r)))) for r in self.roots)
        self.f = dict((-r,M.transpose()) for (r,M) in self.e.iteritems())

    # TODO: check that this is doing what it should 
    def _init_type_C(self, rk, **kwargs):
        self._n = rk*2
        n = self._n
        self.h = dict((self.Q.simple_coroot(i+1),self._e(i,i) - self._e(i+1,i+1) + self._e(n-2-i,n-2-i) - self._e(n-1-i,n-1-i)) for i in range(rk-1))
        self.h[self.Q.simple_coroot(rk)] = self._e(rk-1,rk-1) - self._e(rk,rk)
        def to_pos(a):
            da = dict(a)
            uda = dict( [(i,1) for i in da]+[(2*rk-i,1) for i in da if da[i] == 2])
            return (min(uda)-1,max(uda))
        self.e = dict((r,self._e(i,j)+bool(i+j!=2*rk-1)*int((-1)**(i-j+1))*self._e(2*rk-1-j,2*rk-1-i)) for r in self.roots for (i,j) in [to_pos(r)])
        # we only scale f to avoid square roots
        self.f = dict((-r,M.transpose()/self.form(M,M.transpose())) for (r,M) in self.e.iteritems()) 
       # This is the old code to compute the right quotient ring
       # self.R = PolynomialRing(QQ, [var('x_%d_%d'%(i,j), latex_name='x_{%d,%d}') for i in range(n) for j in range(n) ])
       # self.X = matrix(n, self.R.gens())
       # JJ = matrix(n, lambda i,j: (1 if mod(i,2)==0 else -1 )* (1 if i+j == n-1 else 0))
       # relations = (self.X.transpose()*JJ*self.X-JJ).list()
       # self.I = self.R.ideal(relations).groebner_basis()

    def __repr__(self):
        return "Poisson-Lie group of type " + str(self.cartan_type()) + " with BD-triple " + str(self.gamma)

    def reset_numerical(self, min=1, max=2):
        # The following line does not necessarily produce something in the group
        # R = matrix(self._n, lambda i,j: randint(min,max))
        # The following line produces someting in the group but values are higher
        # Note: we multiply by the denominator of M because the matrices in self.f might have denominators
        # Note: we should be multiplying by the cartan part as well. Since we do not we do not have denominators but we lose on randomness
        R = prod( 1+randint(min,max)*M.denominator()*M for M in sorted(self.e.values()+self.f.values(), key=lambda k: random()))
        self._n_subs = dict(zip(self.X.list(),R.list()))
        self._n_right_vf_subs = dict( (r,dict((x, v.subs(self._n_subs)) for (x,v) in S.iteritems())) for (r,S) in self._right_vf_subs.iteritems())
        self._n_left_vf_subs = dict( (r,dict((x, v.subs(self._n_subs)) for (x,v) in S.iteritems())) for (r,S) in self._left_vf_subs.iteritems())
        self.n_poisson_bracket.clear_cache()
        self.n_left_gradient.clear_cache()
        self.n_right_gradient.clear_cache()
    
    def _e(self,i,j):
        M = matrix(self._n)
        M[i,j] = 1
        return M

    def form(self, X, Y):
        return (X*Y).trace()

    def form_on_weights(self, la, mu):
        return vector(la)*self._form_on_weights*vector(mu)

    @cached_method
    def left_derivative(self, f, a):
        return derivative(f.subs(self._left_vf_subs[a]),t).subs(t=0)

    @cached_method
    def right_derivative(self, f, a):
        return derivative(f.subs(self._right_vf_subs[a]),t).subs(t=0)

    @cached_method
    def right_gradient(self, f):
        #return dict( (r,derivative(f.subs(subs),t).subs(t=0)) for (r,subs) in self._right_vf_subs.iteritems() )
        return dict( (r,self.right_derivative(f,r)) for (r,subs) in self._right_vf_subs.iteritems() )

    @cached_method
    def right_algebra_gradient(self, f):
        w = vector( self.right_derivative(f,a) for a in self.Q.simple_coroots() )
        dfh = sum( d*self.lie.basis()[a] for (d,a) in zip(self._form_on_weights*w, self.Q.simple_coroots()))
        return sum(self.right_derivative(f,a)*self.lie.basis()[-a] for a in chain(*[self.e.keys(),self.f.keys()])) + dfh

    @cached_method
    def left_gradient(self, f):
        #return dict( (r,derivative(f.subs(subs),t).subs(t=0)) for (r,subs) in self._left_vf_subs.iteritems() )
        return dict( (r,self.left_derivative(f,r)) for (r,subs) in self._left_vf_subs.iteritems() )

    @cached_method
    def left_algebra_gradient(self, f):
        w = vector( self.left_derivative(f,a) for a in self.Q.simple_coroots() )
        dfh = sum( d*self.lie.basis()[a] for (d,a) in zip(self._form_on_weights*w, self.Q.simple_coroots()))
        return sum(self.left_derivative(f,a)*self.lie.basis()[-a] for a in chain(*[self.e.keys(),self.f.keys()])) + dfh

    @cached_method
    def algebra_poisson_bracket(self, f, g):
        left = self.lie.killing_form(self.lie.R(self.left_algebra_gradient(f)),self.left_algebra_gradient(g))
        right = self.lie.killing_form(self.lie.R(self.right_algebra_gradient(f)),self.right_algebra_gradient(g))
        out = (left-right)/(2*gcd(self.lie.killing_form_matrix().coefficients()))
        self.algebra_poisson_bracket.set_cache(-out,g,f)
        return out

    @cached_method
    def poisson_bracket(self, f, g):
        # This uses Chari-Pressely conventions
        ldf = self.left_gradient(f)
        rdf = self.right_gradient(f)
        ldg = self.left_gradient(g)
        rdg = self.right_gradient(g)
        out = sum( self._r0[a,b]*(ldf[a]*ldg[b]-rdf[a]*rdg[b]) for a in self.Q.simple_coroots() for b in self.Q.simple_coroots() ) 
        out += sum( ldf[-r]*ldg[r]-rdf[-r]*rdg[r] for r in self.roots ) 
        out += sum( ldf[-a]*ldg[b]-rdf[-a]*rdg[b]-ldf[b]*ldg[-a]+rdf[b]*rdg[-a] for (a,b) in self._related_roots)
        self.poisson_bracket.set_cache(-out,g,f)
        return out

    def omega(self, fns):
        return Matrix(SR, len(fns), lambda i,j: self.poisson_bracket(fns[i],fns[j])/(fns[i]*fns[j]))

    @cached_method
    def n_right_gradient(self, f):
        return dict( (r,derivative(f.subs(subs),t).subs(t=0)) for (r,subs) in self._n_right_vf_subs.iteritems() )

    @cached_method
    def n_left_gradient(self, f):
        return dict( (r,derivative(f.subs(subs),t).subs(t=0)) for (r,subs) in self._n_left_vf_subs.iteritems() )

    @cached_method
    def n_poisson_bracket(self, f, g):
        # This uses Chari-Pressely conventions
        ldf = self.n_left_gradient(f)
        rdf = self.n_right_gradient(f)
        ldg = self.n_left_gradient(g)
        rdg = self.n_right_gradient(g)
        out = sum( self._r0[a,b]*(ldf[a]*ldg[b]-rdf[a]*rdg[b]) for a in self.Q.simple_coroots() for b in self.Q.simple_coroots() ) 
        out += sum( ldf[-r]*ldg[r]-rdf[-r]*rdg[r] for r in self.roots ) 
        out += sum( ldf[-a]*ldg[b]-rdf[-a]*rdg[b]-ldf[b]*ldg[-a]+rdf[b]*rdg[-a] for (a,b) in self._related_roots)
        self.n_poisson_bracket.set_cache(-out,g,f)
        return out
    
    def n_omega(self, fns):
        n_fns = map(lambda f: f.subs(self._n_subs), fns)
        if 0 in n_fns:
            raise ValueError("At least one of the input functions is 0 on the current random element.")
        return Matrix(len(fns), lambda i,j: self.n_poisson_bracket(fns[i],fns[j])/(n_fns[i]*n_fns[j]))

    def is_log_canonical_set(self, fns):
        k = 2
        while True:
            try:
                O = self.n_omega(fns)
            except ValueError:
                self.reset_numerical(1,k)
                k+=1
            else:
                break
        while True:
            self.reset_numerical(1,k)
            try:
                return not O - self.n_omega(fns)
            except ValueError:
                k +=1
            else:
                break

    def set_word(self, word=None, shuffle=False, random_word=True, sign=-1):
        r"""
        
        Set the BFZ word to ``word``.
        If ``word`` is not given computes a (maxiamlly if not shuffle) compatible one.
        If ``shuffle`` do not force trailing minors.
        """
        if not word:
            ct = self.cartan_type()
            W = CoxeterGroup(ct)
            s = W.simple_reflections()

            def random_reduced_word_iter(w):
                while True:
                    ds = w.descents()
                    if not ds:
                        return
                    i = ds[randint(0,len(ds)-1)]
                    w = w.apply_simple_reflection(i, 'right')
                    yield i
    
            if random_word:
                reduced_word = lambda x: list(reversed(list(random_reduced_word_iter(x))))
            else:
                reduced_word = lambda x: x.reduced_word()
           
            if self.gamma:
                a_idx = self.gamma.values() if sign==-1 else self.gamma.keys()
                ct_a = ct.subtype(a_idx)
                ct_a = ct_a.relabel(dict(zip(ct_a.index_set(),a_idx)))
                d_idx = self.gamma.keys() if sign==-1 else self.gamma.values()
                ct_d = ct.subtype(d_idx)
                ct_d = ct_d.relabel(dict(zip(ct_d.index_set(),d_idx)))
                a = reduced_word(CoxeterGroup(ct_a).w0)
                d = reduced_word(CoxeterGroup(ct_d).w0)
            else:
                a = []
                d = []

            b = reduced_word(W.prod(s[i] for i in a) * W.w0)
            c = reduced_word(W.w0 * W.prod(s[i] for i in d))
            first = [ sign*i for i in a+b ]
            second = [ -sign*i for i in c+d ]
            if shuffle:
                word = []
                while first and second:
                    word.append((first if randint(0,1) else second).pop())
                word = first + second + list(reversed(word))
            else: 
                word = first + second
        self.word = word
    
    def bfz_weight_pairs(self):
        la = self.P.fundamental_weights()
        s = self.P.simple_reflections()
        word = self.word

        def act(sequence, weight):
            return reduce(lambda x,f: f(x), sequence, weight)
        
        to_map = [ ([], [s[i] for i in reversed(word) if i > 0], la[k]) for k in range(self.rank(),0,-1) ]
        to_map += [ (reversed([s[-i] for i in word[:k+1] if i < 0]), [s[i] for i in word[k+1:] if i > 0], la[abs(word[k])]) for k in range(len(word))]
        return [ (act(u,l),act(v,l)) for (u,v,l) in to_map ]

    def core_monomial(self, (la,mu), return_side=False):
        out = [(la,mu)]
        first_side = None
        # WARNING this code works for simpliy-laced cases
        # BEWARE of adjusting to use roots and coroots in the future
        bad_la = [ j for j in self.gamma.values() if la[j] > 0 ]
        bad_mu = [ j for j in self.gamma.keys() if mu[j] > 0 ]
        if len(bad_la)+len(bad_mu) > 1:
            raise ValueError("Something is wrong here: the pair (%s, %s) is bad both in rows and columns or in more than on place"%(la,mu))
        elif len(bad_la)+len(bad_mu) == 1:
            w = self.P.fundamental_weights()
            side = Zmod(2)( 0 if bad_la else 1 )
            first_side = side
            i = (bad_la+bad_mu)[0]
            while i:
                if side == 1:
                    new_mu = -w[self.gamma[i]]
                    new_la = new_mu.to_dominant_chamber()
                    i = ([ j for j in self.gamma.values() if new_la[j] > 0 ]+[None])[0]
                else:
                    new_la = -w[self.gamma_inv[i]]
                    new_mu = new_la.to_dominant_chamber()
                    i = ([ j for j in self.gamma.keys() if new_mu[j] > 0 ]+[None])[0]
                side +=1
                out.append((new_la,new_mu))
        if return_side:
            return (first_side, out)
        else: 
            return out

    def core_monomials(self):
        return [ self.core_monomial(wts) for wts in self.bfz_weight_pairs() ]
   
    def core_monomial_2(self, (la,mu)):
        leafs = [(la,mu,0)]
        dg = DiGraph(1)
        wts = {0:(la,mu)}
        w = self.P.fundamental_weights()
        while leafs:
            new_leafs = []
            for (la,mu,ord) in leafs:
                i_la = [ i for i in self.gamma.values() if la[i] > 0 ]
                i_mu = [ i for i in self.gamma.keys() if mu[i] > 0 ]
                if i_la and i_mu:
                    raise ValueError("The pair (%s, %s) is bad both in rows and columns."%(la,mu))
                for i in i_la:
                    new_la = -w[self.gamma_inv[i]]
                    new_mu = new_la.to_dominant_chamber()
                    new_leafs.append((new_la,new_mu,dg.order()))
                    wts[dg.order()] = (new_la, new_mu)
                    dg.add_edge(ord,dg.order(),'rows')
                for i in i_mu:
                    new_mu = -w[self.gamma[i]]
                    new_la = new_mu.to_dominant_chamber()
                    wts[dg.order()] = (new_la, new_mu)
                    new_leafs.append((new_la,new_mu,dg.order()))
                    dg.add_edge(ord,dg.order(),'cols')
            leafs = new_leafs
        return dg,wts

    def canonical_function_2(self, (la,mu)):
        dg,wts = self.core_monomial_2((la,mu))
        mon = dict((i,weights_to_function(wt)) for (i,wt) in wts.iteritems())
        out = [(1,mon)]
        last = out
        while last:
            new = []
            for (sgn,mon) in last:
                for e in dg.edges():
                    der = self.left_derivative if e[2] == 'cols' else self.right_derivative
                    for (a,b) in self._simple_related_roots:
                        (s,t) = (-a,b) if e[2] == 'cols' else  (b,-a)
                        x = der(mon[e[0]], s).expand()
                        y = der(mon[e[1]], t).expand()
                        if x!=0 and y!=0:
                            new_mon = copy(mon)
                            new_mon[e[0]] = x
                            new_mon[e[1]] = y
                            new.append((-1*sgn,new_mon))
            out += new
            last = new
        out = map(lambda (sgn,mon): (sgn,tuple(mon.values())), out)
        return sum(map(lambda (sgn,mon): sgn*prod(mon), uniq(out)))

    def core_monomial_3(self, (la,mu)):
        dg = DiGraph()
        powers = {(la,mu):1}
        w = self.P.fundamental_weights()
        last = [(la,mu)]
        while last:
            new = []
            for (la,mu) in last:
                i_la = [ i for i in self.gamma.values() if la[i] > 0 ]
                i_mu = [ i for i in self.gamma.keys() if mu[i] > 0 ]
                if i_la and i_mu:
                    raise ValueError("The pair (%s, %s) is bad both in rows and columns."%(la,mu))
                for i in i_la:
                    new_la = -w[self.gamma_inv[i]]
                    new_mu = new_la.to_dominant_chamber()
                    dg.add_edge((la,mu),(new_la,new_mu),'rows')
                    if (new_la,new_mu) in powers:
                        powers[(new_la,new_mu)] += 1
                    else:
                        new.append((new_la,new_mu))
                        powers[(new_la,new_mu)] = 1
                for i in i_mu:
                    new_mu = -w[self.gamma[i]]
                    new_la = new_mu.to_dominant_chamber()
                    dg.add_edge((la,mu),(new_la,new_mu),'cols')
                    if (new_la,new_mu) in powers:
                        powers[(new_la,new_mu)] += 1
                    else:
                        new.append((new_la,new_mu))
                        powers[(new_la,new_mu)] = 1
            last = new
        return dg,powers

    def canonical_function_3(self, (la,mu)):
        dg,powers = self.core_monomial_3((la,mu))
        mon = dict((wt,weights_to_function(wt)**powers[wt]) for wt in powers)
        out = [(1,mon)]
        last = out
        while last:
            new = []
            for (sgn,mon) in last:
                for e in dg.edges():
                    der = self.left_derivative if e[2] == 'cols' else self.right_derivative
                    for (a,b) in self._simple_related_roots:
                        (s,t) = (-a,b) if e[2] == 'cols' else  (b,-a)
                        x = der(mon[e[0]], s).expand()
                        y = der(mon[e[1]], t).expand()
                        if x!=0 and y!=0:
                            new_mon = copy(mon)
                            new_mon[e[0]] = x
                            new_mon[e[1]] = y
                            new.append((-1*sgn,new_mon))
            out += new
            last = new
        out = map(lambda (sgn,mon): (sgn,tuple(mon.values())), out)
        return sum(map(lambda (sgn,mon): sgn*prod(mon), uniq(out)))



    def canonical_function(self, (la,mu)):
        side,mon = self.core_monomial((la,mu), return_side=True)
        fns = map(weights_to_function, mon)
        return self._canonicalize_tuple(fns, side)

    def _canonicalize_tuple(self, mon, side):
        out = [(1,tuple(mon))]
        for i in range(len(mon)-1):
            last = out
            out = []
            der = self.left_derivative if is_odd(side+i) else self.right_derivative
            while last:
                new = []
                for (sgn,mon) in last:
                    for (a,b) in self._simple_related_roots:
                        (s,t) = (-a,b) if is_odd(side+i) else (b,-a)
                        x = der(mon[i], s).expand()
                        y = der(mon[i+1], t).expand()
                        if x!=0 and y!=0:
                            new.append((-1*sgn,mon[:i]+(x,y)+mon[i+2:]))
                out += last
                last = new
        return sum(map(lambda (sgn,mon): sgn*prod(mon), uniq(out)))

    def cluster(self):
        return map(self.canonical_function, self.bfz_weight_pairs())

    def good_bfz_weight_pairs(self, sign=-1):
        idx = Zmod(2)( sign > 0 )
        v = vector
        gamma1 = map(lambda (a,_): vector(a.associated_coroot()), self._related_roots)
        gamma2 = map(lambda (_,b): vector(b.associated_coroot()), self._related_roots)
        return [ wt for wt in self.bfz_weight_pairs() if all(sign*a*v(wt[idx+1]) >= 0 for a in gamma1) and all(sign*b*v(wt[idx]) >= 0 for b in gamma2) ]

    def bad_bfz_weight_pairs(self, sign=-1):
        good = self.good_bfz_weight_pairs(sign=sign)
        return [wt for wt in self.bfz_weight_pairs() if wt not in good ]

    def bracket_from_weights(self, i, j):
        wts = self.bfz_weight_pairs()
        if i > j:
            sgn = 1
            i,j = j,i
        else:
            sgn = -1
        la,mu = wts[i]
        xi,nu = wts[j]
        A = self.form_on_weights
        B = self._r0_asymmetry
        return sgn * (( A(la,xi) - A(mu,nu) ) / 2 + vector(la)*B*vector(xi) - vector(mu)*B*vector(nu))

    def frozen_indices(self):
        word = range(-self.rank(),0) + self.word
        possible = range(self.rank())+sorted(dict(map(reversed, enumerate(map(abs,word)))).values())
        monomials = [(i,m) for (i,m) in enumerate(self.core_monomials()) if i in possible]
        tails = map(lambda (i,m): m[1:], monomials)
        actual = [ i for (i,m) in monomials if m not in tails]
        return actual

    def gamma_orbits(self):
        starts = [ i for i in range(1,self.rank()+1) if i not in self.gamma.values() ]
        def path(x):
            while x:
                yield x
                x = self.gamma.get(x,None)
        return dict( (i,list(path(i))) for i in starts)

    def torus_degree(self, m):
        return map(self._PtoT, map(sum,zip(*m)))

    def b_matrix(self):
        # WARNING: this code may make assumptions valid only in type A, refactor this once we start dealing with other types
        cluster = self.cluster()
        n = len(cluster)
        k = 2
        while True:
            try:
                O = matrix(QQ,self.n_omega(cluster))
            except ValueError:
                self.reset_numerical(1,k)
                k += randint(0,1)
            else:
                break
        #if O.is_invertible():
        #    return O.inverse()
        B = matrix(n, lambda i,j: sgn(i-j)*var('b_%d_%d'%tuple(sorted([i,j]))))
        BO = B*O
        # the following line assumes type A (the entries on the diagonal might otherwise be different from 1)
        eqns = [ BO[i,j] == (1 if i==j else 0) for i in range(n) for j in range(n) if i not in self.frozen_indices() ]
        sol = solve(eqns, B.variables())
        if len(sol) != 1:
            print("WARNING wrong number of solutions in first step")
        B = B.subs(sol[0])

        degrees =map(lambda x: map(vector, self.torus_degree(x)), self.core_monomials())
        eqns = flatten([flatten(map(list,map(sum,zip(*map(lambda (c,(v,w)): (c*v,c*w), zip(B[i],degrees)))))) for i in range(n) if i not in self.frozen_indices() ])
        sol = solve(eqns, B.variables())
        if len(sol) != 1:
            print("WARNING wrong number of solutions (%d) in second step"%len(sol))
        B = B.subs(sol[0])
        # remove arrows between frozen
        return B.subs(dict((v,0) for v in B.variables()))

    def b_matrix_fast(self):
        # WARNING: this code makes assumptions valid only in type A, refactor this once we start dealing with other types
        cluster = self.cluster()
        k = 2
        while True:
            try:
                O = matrix(QQ,self.n_omega(cluster))
            except ValueError:
                self.reset_numerical(1,k)
                k += randint(0,1)
            else:
                break
        # try to save some time
        if O.is_invertible():
            return O.inverse()
        m = len(cluster)
        frozen = self.frozen_indices()
        n = m-len(frozen)
        L = identity_matrix(n).columns()
        for i in frozen:
            L.insert(i, vector(QQ,n))
        L = matrix(L).transpose()
        rows = O.solve_left(L).rows()
        degrees = zip(*map(lambda x: map(vector, self.torus_degree(x)), self.core_monomials()))
        def val(w):
            return [sum(map(prod,zip(w,d))) for d in degrees]
        K = O.left_kernel()
        v = tuple(var('t%d'%i) for i in range(K.rank()))
        fix = vector(map(sum,zip(*map(lambda (i,w): v[i]*w , enumerate(O.left_kernel().basis())))))
        wtfix = map(sum,zip(*map(lambda (i,w): [v[i]*u for u in val(w)] , enumerate(O.left_kernel().basis()))))
        B = []
        for r in rows:
            wtr = val(r)
            eqns = flatten(map(tuple,map(sum,zip(wtr,wtfix))))
            sol = solve(eqns, v)
            B.append(r+fix.subs(sol[0]))
        for i in frozen:
            B.insert(i, vector(QQ,m))
        B = matrix(B)
        for i in frozen:
            B[i] = - B.column(i)
        return B

def GSV_matrices(rk, gamma=dict(), with_Y=False, from_one=True):
    r"""

    Returns the matrices used to compute the functions in the GSV seed associated to gamma 

    INPUT:

    - ``rk`` -- integer: the rank of the group

    - ``gamma`` -- a dictionary representing the BD data

    - ``from_one`` -- bool (default: True) whether we label the Dynkin diagram from 1
    
    - ``with_Y`` -- bool (default: False) wheter to return functions on the double

    NOTES:

    - If from_one is true then the resulting matrices have variables indexed from 1 

    LIMITATIONS:

    - Currently only works for gamma oriented with no cycles
    """
    # work in GL_{rk+1}
    n = rk+1

    # rescale gamma to use old code
    if from_one:
        gamma = dict((i-1,j-1) for (i,j) in gamma.iteritems())
        offset = 1
    else:
        offset = 0

    # reverse map
    ammag = dict( (j,i) for (i,j) in gamma.iteritems() )

    # subdivide variable matrices according to gamma
    X = matrix(n, lambda i,j: var('x_%d_%d'%(i+offset,j+offset), latex_name='x_{%d,%d}'%(i+offset,j+offset)))
    if with_Y:
        Y = matrix(n, lambda i,j: var('y_%d_%d'%(i+offset,j+offset), latex_name='y_{%d,%d}'%(i+offset,j+offset)))
    else:
        Y = copy(X)
    
    # principal minors are always good
    L = [X]

    # run ends
    X_runs = [ range(i+1,j+1) for lst in ([k for k in [-1]+range(n) if k not in gamma ],) for (i,j) in zip(lst,lst[1:]) ]
    Y_runs = [ range(i+1,j+1) for lst in ([k for k in [-1]+range(n) if k not in ammag ],) for (i,j) in zip(lst,lst[1:]) ]

    # helper: open path starting at (i,Z)
    def path(i,Z):
        p = [(i,Z)]
        while True:
            j = rk-1-i
            f = (gamma if Z == 'X' else ammag)
            if j in f:
                Z = 'Y' if Z == 'X' else 'X'
                i = f[j]
                p.append((i,Z))
            else:
                break
        return p

    # helper: blocks
    def X_block(i):
        k = [max(lst) for lst in X_runs if i in lst][0]
        l = [min(lst) for lst in X_runs if rk-i in lst][0]
        return X[l:,:k+1]

    def Y_block(i):
        k = [max(lst) for lst in Y_runs if i in lst][0]
        l = [min(lst) for lst in Y_runs if rk-i in lst][0]
        return Y[:k+1,l:]

    # helper: glue blocks
    def glue_along_path(p):
        p = copy(p)
        (i,Z) = p.pop()
        M = (X_block if Z == 'X' else Y_block)(i)
        while p:
            (i,Z) = p.pop()
            if Z == 'X':
                N = X_block(i)
                m = M.ncols()
                n = N.ncols()
                N = block_matrix(2,1,[matrix(m+i+1-N.nrows(),n),N])
                M = block_matrix(2,1,[M,matrix(m+i+1-M.nrows(),m)])
                M = block_matrix(1,2,[M,N])
            else:
                N = Y_block(i)
                m = M.nrows()
                n = N.nrows()
                N = block_matrix(1,2,[matrix(n,m+i+1-N.ncols()),N])
                M = block_matrix(1,2,[M,matrix(m,m+i+1-M.ncols())])
                M = block_matrix(2,1,[M,N])
        # clear subdivision
        M.subdivide()
        return M

    L += [ glue_along_path(path(i,'X')) for i in range(rk) if i not in gamma ]
    L += [ glue_along_path(path(i,'Y')) for i in range(rk) if i not in ammag ]

    return L


def GSV_functions(rk, gamma=dict(), with_Y=False, from_one=True):
    r"""

    Returns a dictionary of the functions in the GSV seed associated to gamma labeled by their position in the grid

    INPUT:

    - ``rk`` -- integer: the rank of the group

    - ``gamma`` -- a dictionary representing the BD data

    - ``from_one`` -- bool (default: True) whether we label the Dynkin diagram from 1
    
    - ``with_Y`` -- bool (default: False) wheter to return functions on the double

    NOTES:

    - If from_one is true then the resulting polynomials have variables indexed from 1 and the positions on the grid start from 1

    LIMITATIONS:

    - Currently only works for gamma oriented with no cycles
    """
    L = GSV_matrices(rk, gamma=gamma, with_Y=with_Y, from_one=from_one)
    F = {}
    for M in L:
        for k in range(M.nrows()):
            f = M[k:,k:].det()
            (i,j) = map(int, str(M[k,k]).split('_')[1:])
            F[(i,j)] = f
    return F

def BFZ_matrix(ct, word, reindex=False, quiver=False):
    ct = CartanType(ct)
    A = ct.cartan_matrix()
    lenw = len(word)
    word = dict([(-i,-i) for i in range(ct.rank(),0,-1)] + list(enumerate(word,1)))
    def nxt(k):
        return min( [l for l in word if l > k and abs(word[l]) == abs(word[k])] + [lenw+1] )
    cols_idx = dict(enumerate([ i for i in range(1,lenw+1) if nxt(i) <= lenw ]))
    if reindex or quiver:
        rows_idx = dict( list(cols_idx.iteritems()) + list(enumerate( (i for i in word if i not in cols_idx.values()), len(cols_idx))))
    else:
        rows_idx = dict((i+ct.rank() if i<0 else i+ct.rank()-1,i) for i in word)
    def bkl(k,l):
        p = max(k,l)
        q = min(nxt(k),nxt(l))
        if p == q:
            return -sign(k-l)*sign(word[p])
        elif p<q and sign(word[p])*sign(word[q])*(k-l)*(nxt(k)-nxt(l)) > 0:
            return -sign(k-l)*sign(word[p])*A[abs(word[p])-1,abs(word[q])-1]
        else:
            return 0
    B = matrix(len(rows_idx),len(cols_idx), lambda i,j: bkl(rows_idx[i],cols_idx[j]))
    if quiver:
        Q = ClusterQuiver(B)
        Q._digraph.layout(pos={i:(10*rows_idx[i],10*abs(word[rows_idx[i]])) for i in rows_idx},save_pos=True)
        return Q
    if reindex:
        return (B,rows_idx)
    else:
        return B

def BFZ_seed(ct, word):
    ct = CartanType(ct)
    fns = dict(zip(range(-ct.rank(),0)+range(1,len(word)+1), BFZ_minors_weight_pairs(ct, word)))
    (B,idx) = BFZ_matrix(ct, word, reindex=True)
    return (B, [fns[idx[i]] for i in range(len(fns))])

def BFZ_minors_weight_pairs(ct, word):
    r"""

    Returns the minors in the BFZ seed attached corresponding to the word ii
    
    INPUT: 

    - ``ct`` -- integer: a cartan type (has to be finite)

    - ``word`` -- list: a doubled reduced word on the letters [1 .. ct.rank()] + [-cr.rank() .. -1]

    NOTES:
    
    - This function returns a list of pairs. In each pair the first entry labels the rows of the minor, i.e. the bottom index.
    """
    ct = CartanType(ct)
    R = ct.root_system()
    P = R.weight_lattice()
    P.print_options(latex_prefix='\\omega')
    la = P.fundamental_weights()
    s = P.simple_reflections()

    wt_pairs = [ (reduce(lambda f,g: g(f), [s[i] for i in reversed(word) if i > 0], la[k]), la[k]) for k in range(ct.rank(),0,-1) ]
    wt_pairs += [
            (reduce(lambda f,g: g(f), [s[i] for i in word[k+1:] if i > 0], la[abs(word[k])]),
                reduce(lambda f,g: g(f), reversed([s[-i] for i in word[:k+1] if i < 0]), la[abs(word[k])]))
            for k in range(len(word))
            ]
    return wt_pairs

def BFZ_minors_type_A_indices(rk, word, from_one=True):
    r"""

    Returns the minors of the BFZ seed in type A expressing weights as matrix indices

    INPUT:

    - ``rk`` -- integer: the rank of the group

    - ``word`` -- list: a doubled reduced word on the letters [1 .. rk] + [-rk .. -1]

    - ``from_one`` -- bool (default ``True``) whether to start counting from 1 (we need to count from 0 internally)

    NOTES:
    
    - This function returns a list of pairs. In each pair the first entry labels the rows of the minor.
    """
    wt_pairs = BFZ_minors_weight_pairs(['A',rk], word)
    return map( lambda (la,mu): (_to_indices(la, from_one=from_one), _to_indices(mu, from_one=from_one)), wt_pairs)

def BFZ_minors_type_A_functions(rk, word, from_one=True):
    r"""

    Returns the minors of the BFZ seed in type A expressing weights as polynomials

    INPUT:

    - ``rk`` -- integer: the rank of the group

    - ``word`` -- list: a doubled reduced word on the letters [1 .. rk] + [-rk .. -1]

    - ``from_one`` -- bool (default ``True``) whether to start counting from 1

    NOTES:
    
    - This function returns a list of pairs. In each pair the first entry labels the rows of the minor.
    """
    offset = 1 if from_one else 0
    X = matrix(rk+1, lambda i,j: var('x_%d_%d'%(i+offset,j+offset), latex_name='x_{%d,%d}'%(i+offset,j+offset)))
    wt_pairs = BFZ_minors_type_A_indices(rk, word, from_one=False)
    return tuple(X.matrix_from_rows_and_columns(ro,co).det().expand() for (ro,co) in wt_pairs)

#def find_bad_BFZ_weight_pairs(ct, gamma, word=None):
#    r"""
#
#    Return the minors in the BFZ seed attached to word that are not log canonical with respect to gamma.
#
#    INPUT:
#
#    - ``ct`` -- a Cartan Type
#
#    - ``gamma`` -- a non-empty BD funtrion
#
#    - ``word`` -- a double reduced word. If none is provided we use a maxiamlly compatible one
#
#    WARNING:
#
#    If no word is provided the result may not be the expected one because the word used may be different than the one ou are thinking of!!!!
#    """
#    if not word:
#        word = compatible_word(ct, gamma)
#    wt_pairs = BFZ_minors_weight_pairs(ct, word)
#    related_roots =  PoissonLieGroup(ct, gamma=gamma)._related_roots
#    column_roots = set([ a for (a,b) in related_roots ])
#    row_roots = set([ b for (a,b) in related_roots ])
#    return [(w1,w2) for (w1,w2) in wt_pairs if any(vector(a)*vector(w1)>0 for a in row_roots) or any(vector(a)*vector(w2)>0 for a in column_roots)]
#
#def find_good_BFZ_weight_pairs(ct, gamma, word=None):
#    r"""
#
#    Return the minors in the BFZ seed attached to word that are log canonical with respect to gamma.
#
#    INPUT:
#
#    - ``ct`` -- a Cartan Type
#
#    - ``gamma`` -- a non-empty BD funtrion
#
#    - ``word`` -- a double reduced word. If none is provided we use a maxiamlly compatible one
#
#    WARNING:
#
#    If no word is provided the result may not be the expected one because the word used may be different than the one ou are thinking of!!!!
#    """
#    if not word:
#        word = compatible_word(ct, gamma)
#    wt_pairs = BFZ_minors_weight_pairs(ct, word)
#    related_roots =  PoissonLieGroup(ct, gamma=gamma)._related_roots
#    column_roots = set([ a for (a,b) in related_roots ])
#    row_roots = set([ b for (a,b) in related_roots ])
#    return [(w1,w2) for (w1,w2) in wt_pairs if all(vector(a)*vector(w1)<=0 for a in row_roots) and all(vector(a)*vector(w2)<=0 for a in column_roots)]

def _to_indices(wt, from_one=False):
    r"""

    Returns the weight wt as a list of indices

    INPUT:

    - ``wt`` -- a weight in a fundamental representation of SL

    - ``from_one`` -- bool (default ``True``) whether to start counting from 1 (we need to count from 0 internally)
    """
    v = vector(wt)
    A = [ range(k+1) for k in range(len(v)) if v[k] == v.coefficients()[-1]]
    B = [ range(k+1) for k in range(len(v)) if v[k] == -v.coefficients()[-1]]
    idx = flatten(A)
    for i in flatten(B):
        idx.remove(i)
    if v.coefficients()[-1] == -1:
        idx = [ i for i in range(len(v)+1) if i not in idx ]
    if from_one:
        return [i+1 for i in idx]
    else:
        return idx

def weights_to_function((rows,cols)):
    r"""
    Return the polynomial version of the minor with the given rows and columns

    This is only a Type-A function

    We count from 1
    """

    rk = rows.parent().rank()
    X = matrix(rk+1, lambda i,j: var('x_%d_%d'%(i+1,j+1), latex_name='x_{%d,%d}'%(i+1,j+1)))
    return X.matrix_from_rows_and_columns(_to_indices(rows),_to_indices(cols)).det().expand()
