attach("sage.lib/lie_derivative.sage")
attach("sage.lib/linear_problem.sage")
attach("sage.lib/polynomial.sage")

# ============================================================================================================================= #

# ------------------------- #
# Multiple Hopf Bifurcation #
# ------------------------- #

# for a given vector field 'H' in variables 'V' with parameters 'P', a multiplicity 'm'
# find a sufficient condition for that 
#                   hopf bifurcations of multiplicit 'm' occur for the vector field 'H'
# where 'H' satifies 'H(0) = 0' that is 'H' vanish at zero for any parameter value
#        and 'H' is a list consisting of homogeneous polynomial vectors such that 
#            'H = [H1, H2, H3, ..., H(2m+1)]', 
#        and  'Hj' is a homogeneous polynomial vector of degree '2m+1'

def hopf(m, V, P, H, PR, PF, VR, log_opt = 0, qe_opt = 'none', gb_opt = 'none', annpol_opt = 0, pol_opt = [], output_opt = 0, ker_method_opt = 0, ker_opt = 0):
    if pol_opt == []:
      # Preparation
      LenV = len(V)
      LenH = len(H)
      if LenH < 2*m + 1:
         H = [h for h in H]
         H = H+[Matrix([[VR(0)] for j in range(LenV)]) for i in range(2*m+1-LenH)]
      if log_opt >= 1:
         TIME0 = cputime(subprocesses=True)
      # Liu Criterion
      # 1. generate the representing matrix 'DLin' of 
      #    the linear part 'H1' of 'H' with dummy parameter 'DP' for a curve
      DP = var(['crv%d' % i for i in range(len(P)+1)])
      DR = PolynomialRing(QQ, P+DP)
      DF = Frac(DR)
      Lin = Matrix([[PF(h[0].coefficient({v:1})) for v in V] for h in H[0]])
      DLin = Matrix([[ll.subs({P[i]:DF(P[i])+DF(DP[0]*DP[i+1]) for i in range(len(P))}) for ll in l] for l in Lin])
      # 3. generate the Routh-Hurwitz deteminants 'RH' of 'DLin'
      DMin = DLin.minpoly()
      RH = routh_hurwitz_determinant(DMin)
      # 4. construct the condition 
      #    'Liu = [[polynomials equal to 0], 
      #            [polynomials greater than 0], 
      #            [polynomials not equal to 0]]'
      Liu = [[DR(RH[1][-1]).subs({DR(DP[0]):0})],
             [DR(RH[0]).subs({DR(DP[0]):0})]+[DR(r).subs({DR(DP[0]):0}) for r in RH[1][0:-2]],
	     [DR(diff(RH[1][-1], DP[0]).subs({DR(DP[0]):0})).subs({DR(DP[0]):0})]]
      # 5. simplify the condition over the complex number field
      if gb_opt == 'none':
         Liu0 = ((DR.ideal(Liu[0]).saturation(DR.ideal(Liu[1]))[0]).radical()).groebner_basis()
      if gb_opt == 'default':
         Liu0 = ((DR.ideal(Liu[0]).saturation(DR.ideal(Liu[1]))[0]).radical()).groebner_basis()
      elif gb_opt == 'magma':
         Liu0 = ((DR.ideal(Liu[0]).saturation(DR.ideal(Liu[1]))[0]).radical()).groebner_basis(algorithm='magma:GroebnerBasis')
      Liu0 = [PR(f) for f in Liu0]
      Liu1 = [PR(f) for f in Liu[1]]
      Liu2 = [Liu[2][0].reduce(Liu0)]
      Liu  = [Liu0, Liu1, Liu2];
      LIU0 = Liu[0]
      if log_opt >= 1:
         TIME1 = cputime(subprocesses=True)
         print('Liu Condition:', TIME1-TIME0, 'seconds')
         if log_opt == 2:
            print('Liu Condition:', Liu)
      # 6. simplify the representation of 'H'
      H = [Matrix([liu_simplification(h[0], Liu0, VR)] for h in Hi) for Hi in H]
      # 7. simplify the Liu condition by eliminating quantifiers
      if qe_opt == 'none':  # not eliminate quantifiers
         Liu0_ = [str(l) + '=' + str(0) for l in Liu[0]]
         Liu1_ = [str(l) + '>' + str(0) for l in Liu[1]]
         Liu2_ = [str(l) + '<>' + str(0) for l in Liu[2]]
         Liu_ = Liu0_ + Liu1_ + Liu2_
         T_ = [l.variables() for l in Liu[0]] + [l.variables() for l in Liu[1]] + [l.variables() for l in Liu[2]]
         T_ = [t1 for t in T_ for t1 in t]
         DP_ = [d for d in DP if d in T_]
         psi_ = 'Ex(' + str(DP_) + ', And(' + ','.join(Liu_) + '))'    
         if log_opt >= 1:
            TIME2 = cputime(subprocesses=True)
            print('QE Liu Condition:', TIME2-TIME1, 'seconds')
            if log_opt == 2:
               print('QE Liu Condition:', psi_)
      elif qe_opt == 'qep': # eliminate quantifiers by QEPCAD
         T_ = [l.variables() for l in Liu[0]] + [l.variables() for l in Liu[1]] + [l.variables() for l in Liu[2]]
         T_ = [t1 for t in T_ for t1 in t]
         DP_ = [d for d in DP if d in T_]
         Liu_ = [[SR(l) for l in Liu[0]], [SR(l) for l in Liu[1]], [SR(l) for l in Liu[2]]]
         qf = qepcad_formula
         psi_ = qf.and_([l == 0 for l in Liu_[0]] + [l > 0 for l in Liu_[1]] + [l != 0 for l in Liu_[2]])
         psi_ = qf.exists(DP_, psi_)
         psi_ = qepcad(psi_)
         if log_opt >= 1:
            TIME2 = cputime(subprocesses=True)
            print('QE Liu Condition:', TIME2-TIME1, 'seconds')
            if log_opt == 2:
               print('QE Liu Condition:', psi_)
      elif qe_opt == 'math': # eliminate quantifiers by Reduce on Mathematica
         T_ = [l.variables() for l in Liu[0]] + [l.variables() for l in Liu[1]] + [l.variables() for l in Liu[2]]
         T_ = [t1 for t in T_ for t1 in t]
         DP_ = mathematica([SR(d) for d in DP if d in T_])
         Liu0_ = [str(l) + '==' + str(0) for l in Liu[0]]
         Liu1_ = [str(l) + '>' + str(0) for l in Liu[1]]
         Liu2_ = [str(l) + '!=' + str(0) for l in Liu[2]]
         Liu_ = Liu0_ + Liu1_ + Liu2_
         psi_ = 'Exists[' + str(DP_) + ', And[' + ','.join(Liu_) + ']]'    
         psi_ = mathematica('Reduce[' + psi_ + ']')
         if log_opt >= 1:
            TIME2 = cputime(subprocesses=True)
            print('QE Liu Condition:', TIME2-TIME1, 'seconds')
            if log_opt == 2:
               print('QE Liu Condition:', psi_)
      elif qe_opt == 'cgs+syn': # eliminate quantifiers by SyNRAC and CGSQE on Maple
         T_ = [l.variables() for l in Liu[0]] + [l.variables() for l in Liu[1]] + [l.variables() for l in Liu[2]]
         T_ = [t1 for t in T_ for t1 in t]
         DP_ = [d for d in DP if d in T_]
         Liu0_ = [str(l) + '=' + str(0) for l in Liu[0]]
         Liu1_ = [str(l) + '>' + str(0) for l in Liu[1]]
         Liu2_ = [str(l) + '<>' + str(0) for l in Liu[2]]
         Liu_ = Liu0_ + Liu1_ + Liu2_
         psi_ = 'Ex(' + str(DP_) + ', And(' + ','.join(Liu_) + '))'    
         psi_ = maple('qe(' + psi_ + ', use_cgsqe = true)')
         if log_opt >= 1:
            TIME2 = cputime(subprocesses=True)
            print('QE Liu Condition:', TIME2-TIME1, 'seconds')
            if log_opt == 2:
               print('QE Liu Condition:', psi_)
    else:
        if log_opt >= 1:
            TIME0 = cputime(subprocesses=True)
        PPP = pol_opt;
        Liu = [[PF(l) for l in liu] for liu in PPP];
        LIU0 = Liu[0]
        Liu0_ = [str(l) + '=' + str(0) for l in Liu[0]]
        Liu1_ = [str(l) + '>' + str(0) for l in Liu[1]]
        Liu2_ = [str(l) + '<>' + str(0) for l in Liu[2]]
        Liu_ = Liu0_ + Liu1_ + Liu2_
        psi_ = 'And(' + ','.join(Liu_) + ')'
    # Lyapunov Quantity
    LIU0_D = [PR(l) for l in LIU0]
    LIU0_D = [DEC.gens() for DEC in PR.ideal(LIU0_D).radical().primary_decomposition()]
    #LIU0_D = list((LIU0[0].numerator()).factor())
    #LIU0_D = [l[0] for l in LIU0_D]; 
    LYPss  = []
    LIUs   = []
    KER2s = []
    PHIss = []
    for ll in LIU0_D:
        # 0. simplify the representation of 'H'
        if log_opt >= 1:
           TIME3 = cputime(subprocesses=True)
        Liu = [PR.ideal(ll).groebner_basis(),Liu[1],Liu[2]]
        LIUs.append(Liu)
        H = [Matrix([liu_simplification(h[0], Liu[0], VR)] for h in Hi) for Hi in H]
        # 1. find a generator 'Ker2' of the kernek of the Lie derivative in 'H1' of 
        #                                        homogeneous polynomials of degree '2'    
        REP2 = linear_lie_derivative_representation(2, P, V, H[0])
        if log_opt >= 1:
           T_TIME4 = cputime(subprocesses=True)
           print('2 -th Lie Derivative:', T_TIME4-TIME3, 'seconds')
           if log_opt == 2:
              print('2 -th Lie Derivative:')
              print(REP2[0])
        BSS2 = REP2[1]; 
        REP2 = REP2[0]
        if ker_opt == 0:
          KER2 = linear_lie_derivative_kernel(P, V, REP2, Liu[0], opt = ker_method_opt).basis()[0]
        else: KER2 = ker_opt
        KER2V = sum([KER2[i]*BSS2[i] for i in range(len(BSS2))])
        KER2s.append(KER2V)
        KER2V_ = 0
        for i in range(len(V)):
          SUB2 = [1 if i == j else 0 for j in range(len(V))]
          KER2V_ = KER2V.subs({V[j]:SUB2[j] for j in range(len(V))})
          if KER2V_ != 0:
            break
        if KER2V_ == 0:
          for i in range(len(V)):
            for j in range(len(V)):
              SUB2 = [1 if i == k or j == k else 0 for k in range(len(V))]
              KER2V_ = KER2V.subs({V[j]:SUB2[j] for j in range(len(V))})
              if KER2V_ != 0:
                break
        if KER2V_ == 0:
           print('Theoretical Error: KER2.subs(V[i]) == 0 for i in range(len(V))')
           return 0    
        if log_opt >= 1:
           TIME4 = cputime(subprocesses=True)
           print('Quadratic Kernel Generator:', TIME4-T_TIME4, 'seconds')
           if log_opt == 2:
              print('Quadratic Kernel Generator:', KER2V)
        # 2. find an annihilating polynomial 'Ann1' of the linear Lie derivative of 
        #                                      homogeneous polynomials of degree '1'
        REP1 = linear_lie_derivative_representation(1, P, V, H[0])[0]
        ANN1 = REP1.minpoly()
        ANNi = ANN1
        # 3. for 'i = 2, ..., 2*m+1'
        PHIs = [KER2V]
        LYPs = []
        Liu0_ = [l for l in Liu[0]]; 
        x = parent(ANNi).gen()
        for i in range(1, 2*m+2):
            # 3-1. find annihilating polynomials of the linear Lie derivatives of 
            #                              homogeneous polynomials of degree '2,...,2*m+1'
            if log_opt >= 1:
               TIME5 = cputime(subprocesses=True)
            if i == 1:
               REPi = REP2
               BSSi = BSS2
            else:
               REPi = linear_lie_derivative_representation(i+1, P, V, H[0])
               BSSi = REPi[1]
               REPi = REPi[0].transpose()
            if log_opt >= 1:
               T_TIME6 = cputime(subprocesses=True)
               print(i+1, '-th Lie Derivative:', T_TIME6-TIME5, 'seconds')
               if log_opt == 2:
                  print(i+1, '-th Lie Derivative:')
                  print(REPi)
            if annpol_opt == 0: # use 'minpoly',
               ANNi = REPi.minpoly()
            else: # use 'lib/lie_derivative'
               ANNi = linear_lie_derivative_annihilating_polynomial(i+1, P, V, ANN1, ANNi, PF, lib_log_opt = i)
               COEi = ANNi.coefficients(sparse=False)
               COEi = [PF(c) for c in COEi] 
               ANNi = sum([COEi[i]*x**i for i in range(len(COEi))])
            if log_opt >= 1:
               TIME6 = cputime(subprocesses=True)
               print(i+1, '-th Annihilating Polynomial:', TIME6-T_TIME6, 'seconds')
               if log_opt == 2:
                  print(i+1, '-th Annihilating Polynomial:', ANNi)
            # 3-2. solve linear problems
            if i > 1:
               Ui = -sum([lie_derivative(V, H[(i-1)-j], PHIs[j]) for j in range(i-1)])
               Ui = Matrix([[Ui.coefficient(t)] for t in BSSi])
               if i != 2*m+1:
                    if log_opt >= 1:
                       VKi = linear_problem_solver(P, REPi, Ui, ANNi, vector_opt=0, lin_log_opt=i)
                    else:
                       VKi = linear_problem_solver(P, REPi, Ui, ANNi, vector_opt=0, lin_log_opt=0)
                    Vi = sum([VKi[0][j]*BSSi[j] for j in range(len(BSSi))]); 
                    PHIs.append(Vi)
                    if i%2 != 0:
                       if ker_opt == 0:
                         KERi = sum([VKi[1][j]*BSSi[j] for j in range(len(BSSi))])
                       else:
                         KERi = sum([VKi[0][j]*BSSi[j] for j in range(len(BSSi))])
                       LPYi = (KERi.subs({V[j]:SUB2[j] for j in range(len(V))})/((KER2V_/2)**(i//2-1)))
                       LPYiN = (LPYi.numerator()).reduce(Liu0_)
                       LPYiD = LPYi.denominator()
                       LYPs.append(LPYiN/LPYiD)
                       if log_opt >= 1:
                          TTTIME7 = cputime(subprocesses=True)
                          print(">>>>>>>", i+1-(3+(i-3)), '-th K:', TTTIME7-TIME6, 'seconds')
                          if log_opt >= 2:
                            print(">>>>>>>", i+1-(3+(i-3)), '-th K:', LPYiN/LPYiD)
 
                       Liu0_.append(LPYiN); 
                       if gb_opt == 'none':
                           Liu0_ = Liu0_
                       elif gb_opt == 'default':
                           Liu0_ = (PR.ideal(Liu0_).radical()).groebner_basis()
                       elif gb_opt == 'magma':
                          Liu0_ = (PR.ideal(Liu0_).radical()).groebner_basis(algorithm='magma:GroebnerBasis')
                       Liu0_ = [f for f in Liu0_]
                       H = [Matrix([liu_simplification(h[0], Liu0_, VR)] for h in Hi) for Hi in H]
               else:
                    if log_opt >= 1:
                       VKi = linear_problem_solver(P, REPi, Ui, ANNi, vector_opt=1, lin_log_opt=i)
                    else:
                       VKi = linear_problem_solver(P, REPi, Ui, ANNi, vector_opt=1, lin_log_opt=0)
                    if ker_opt == 0:
                         KERi = sum([VKi[1][j]*BSSi[j] for j in range(len(BSSi))])
                    else:
                         KERi = sum([VKi[0][j]*BSSi[j] for j in range(len(BSSi))])
                    LPYi = KERi.subs({V[j]:SUB2[j] for j in range(len(V))})/((KER2V_/2)**(i//2-1))
                    LPYiN = (LPYi.numerator()).reduce(Liu0_)
                    LPYiD = LPYi.denominator()
                    LYPs.append(LPYiN/LPYiD)
                    if log_opt >= 1:
                      TTTIME7 = cputime(subprocesses=True)
                      print(">>>>>>>", i+1-(3+(i-4)), '-th constant Q:', TTTIME7-TIME6, 'seconds')
                      if log_opt >= 2:
                        print(">>>>>>>", i+1-(3+(i-4)), '-th constant Q:', LPYiN/LPYiD)
            if log_opt >= 1:
               TIME7 = cputime(subprocesses=True)
               print(i+1, '-th Linear Algebra Problem:', TIME6-TIME5, 'seconds')
               print(i+1, '-th For Loop:', TIME7-TIME5, 'seconds')
        LYPss.append(LYPs)
        PHIss.append(PHIs)
    if log_opt >= 1:
       TIME8 = cputime(subprocesses=True)
       print('All:', TIME8-TIME0, 'seconds')
    if output_opt == 0:
       if pol_opt == []:
         print('a condition for a Hopf point: ', psi_)
       else: print('')
       for j in range(len(LYPss)):
         for i in range(len(LYPss[j])):
             print('a Lyapunov quantity (', i+1, '): ', LYPss[j][i])
       return(0)
    else:
       return([psi_, LYPss, LIUs, DP_, KER2s, PHIss])

# ----------------------------------------------------------------------------------------------------------------------------- #

def liu_simplification(h, Liu0, VR):
    T = VR(h).monomials()
    C = [h.monomial_coefficient(t) for t in T]
    D = [c.denominator() for c in C]
    N = [c.numerator() for c in C]
    N = [n.reduce(Liu0) for n in N]
    C = [N[i]/D[i] for i in range(len(C))]
    return sum([C[i]*T[i] for i in range(len(C))])

# ============================================================================================================================= #