import numpy as np # Library with mathematical computation tools
import itertools # Iterative library with combinatoric tools
import math # To use floor

''' Values of the generating functions (completion of the curve) over exceptional divisors. '''

# Function to obtain an integral number 
# Used for the number of relevant exceptional divisors and 
# the number of curves in the maximal contact completion 
# string is a string for the elements to count

def getIntegralNumber(string):
    print(string)
    int_number = input()
    while True:
        try:
            int_number = int(int_number)
            break
        except:
            print("You have to type an integer.")
            int_number = input()
    return int_number


def get_curve_data():
    # Get the number of RED
    string_number_divisors = "How many relevant exceptional divisors does the curve have?"
    number_divisors = getIntegralNumber(string_number_divisors)


    # Get the number of MCC
    string_number_functions = "How many functions do you need to generate the maximal contact completion of the curve?"
    number_functions = getIntegralNumber(string_number_functions)

    
    # Get the number of functions needed for monomials
    string_number_functions_monomial = "Consider the dual graph of the minimal embedded resolution of the curve. " + \
    "How many functions do you need to complete it " + \
    "(one for each exceptional divisor of valency one)?"
    number_functions_monomial = getIntegralNumber(string_number_functions_monomial)

    # Get the number of components
    string_number_components = \
    "How many components does the curve have?"
    number_components = getIntegralNumber(string_number_components)


    # Get the order of vanishing of curves in the completion along RED and log-discrepancies
    print("\nNow, for each divisor, introduce the orders of vanishing and its log-discrepancy." \
         "\nLet z_0 represent the transversal branch to the curve " + \
          "and the components of the curve be the last " + str(number_components) + \
         ": z_" + str(number_functions - number_components) + ", z_" + str(number_functions - 1) + \
         ". Also, let z_0, ..., z_" + str(number_functions_monomial) + \
          " be the curves needed to complete the dual graph. " + \
          "Introduce the orders of vanishing separated by single spaces and after input the log-discrepancy." + \
          " \nThe curve is represented by z_" + str(number_functions) + \
          ".\nThe divisors vary from E_1 to E_" + str(number_divisors))

    # the number of columns is number of functions + 2, 
    # since we add the order of vanishing of the whole curve and the log-discrepancy
    Matrix_orders_vanishing = [[0 for x in range(number_functions + 2)] for y in range(number_divisors)] 

    for div in range(number_divisors):
        Matrix_orders_vanishing[div][0:number_functions] = \
        [int(item) for item in input("\nEnter the orders of vanishing along E_" + str(div + 1) + "\n").split()]

        Matrix_orders_vanishing[div][number_functions] = \
        sum(Matrix_orders_vanishing[div][number_functions - number_components : number_functions])

        string_log_discrepancy = "Introduce the log-discrepancy of the divisor E_" + str(div+1)
        Matrix_orders_vanishing[div][number_functions + 1] = getIntegralNumber(string_log_discrepancy)
        

    # Print the collected data as a matrix
    print("\nThe matrix of orders of vanishing and log-discrepancies is:")
    print(np.array(Matrix_orders_vanishing))


    return (number_divisors,number_functions,number_functions_monomial,number_components,Matrix_orders_vanishing)
    


''' Compute the Milnor numbers of the components and total curve 
    and the number of jumping numbers counted with multiplicity. 
    Milnor(C) = sum(1,r) Milnor(C_i) + sum(i < j) (Ci,Cj) - number_components + 1 .
    number of jn counted with multiplicity = (Milnor(C) - number_components + 1)/2 .'''

def get_Milnor(number_divisors,number_functions,number_functions_monomial,number_components,Matrix_orders_vanishing):
    # Instructions for user
    print("\nNow introduce the number i of the divisor E_i such that it is the last relevant exceptional divisor" +\
        " of the irreducible curve C_i, generated by the germ z_" + str(number_functions - number_components - 1) + "+ i.")

    corresp_comp_red = [0 for x in range(number_components)]
    Milnor_number = [0 for x in range(number_components)]
    number_jn_multiplicity = 0

    print("Initially number_jn_multiplicity = " + str(number_jn_multiplicity))

    for comp in range(number_components):
        # Obtain the correspondence between curves and relevant exceptional divisors.
        string_corresp_comp_red = "Introduce the number chosen for the last relevant exceptional divisor" +\
        " of the irreducible curve C_" + str(comp + 1) + " = V(z_" +\
        str(number_functions - number_components + comp) + "):"

        corresp_comp_red[comp] = getIntegralNumber(string_corresp_comp_red) - 1

        # Compute the corresponding Milnor number as 
        # order of vanishing of the curve - log-discrepancy + 1
        Milnor_number[comp] = \
        Matrix_orders_vanishing[corresp_comp_red[comp]][number_functions - number_components + comp] - \
        Matrix_orders_vanishing[corresp_comp_red[comp]][number_functions + 1] + 1

        number_jn_multiplicity += Milnor_number[comp]
        print("The Milnor number of the curve C_" + str(comp + 1) + " = V(z_" +\
        str(number_functions - number_components + comp) + ") is: " + str(Milnor_number[comp]))

        for comp2 in range(comp+1,number_components):
            number_jn_multiplicity += \
            2 * Matrix_orders_vanishing[corresp_comp_red[comp]][number_functions - number_components + comp2]        


    Milnor_number_C = number_jn_multiplicity - number_components + 1
    print("The Milnor number of the curve C is: " + str(Milnor_number_C))

    number_jn_multiplicity =  int((Milnor_number_C - number_components + 1) / 2)
    print("The amount of jumping numbers counted with multiplicity is: " + str(number_jn_multiplicity))


    return(number_jn_multiplicity)



''' Generate a matrix with exponent values, 
    the vanishing of the corresponding monomial over divisors and 
    minimal vanishing. '''

def create_jumping_numbers(number_divisors,number_functions,number_functions_monomial,number_components,Matrix_orders_vanishing):
    
    # Get the maximum exponents for each MCC
    bound_exponents = [0 for x in range(number_functions_monomial)]
    
    for func in range(number_functions_monomial):
        # We could do: string = "\nIntroduce the bound exponent for the function z_" + str(func) + "if -1, then we compute it"
        # bound_exponents[j] = getIntegralNumber(string)
        bounding_cond = \
        [(Matrix_orders_vanishing[div][number_functions] - Matrix_orders_vanishing[div][number_functions + 1])/\
         Matrix_orders_vanishing[div][func] for div in range(number_divisors)]
        bound_exponents[func] = math.floor(np.amax(bounding_cond))

    print("\nThe maximum (inclusive) exponents considered for monomials in the z_i are:")
    print(bound_exponents)



    # Generate a (np.prod(bound_exponents)) x (number_functions_monomial + number_divisors + 1) matrix for which:
    # first (number_functions) columns represent a monomial in the semi-roots
    # next number_divisors, from number_functions to number_divisors-1, give the jn condition associated to the monomial
    # the last column is the jumping number associated to the considered monomial
    ''' Still has to be optimized!!! '''
    
    MonomioGen = [[0 for x in range(number_functions_monomial + number_divisors + 1)] \
                  for y in range(np.prod(bound_exponents))]

    # compute the jn conditions
    i = 0
    for mon in itertools.product(*map(range, np.array(bound_exponents))):
        # monomial
        MonomioGen[i][0:number_functions_monomial] = np.array(mon)

        for div in range(number_divisors):
            # the next represent the corresponding condition in the divisor E_(div + 1)
            # monomial times (@) order of vanishing of each function over the vanishing of the curve
            monomial_array_ld = np.append(np.array(MonomioGen[i][:number_functions_monomial]),1)

            order_van_array = \
            np.append(Matrix_orders_vanishing[div][:number_functions_monomial], \
                      Matrix_orders_vanishing[div][number_functions + 1])

            MonomioGen[i][number_functions_monomial + div] = \
            (monomial_array_ld @ order_van_array) / Matrix_orders_vanishing[div][number_functions] 

        # The jumping number is the minimal condition
        MonomioGen[i][number_functions_monomial + number_divisors] = \
        np.amin(MonomioGen[i][number_functions_monomial : number_functions_monomial + number_divisors])

        i += 1



    # Turn the matrix into an array for use of numpy functions and sort the matrix by the corresponding jumping number.
    MonomioGen = np.array(MonomioGen)
    MonomioGen = MonomioGen[MonomioGen[:,number_functions_monomial + number_divisors].argsort()]


    # Make a matrix (array) with elements of minimum vanishing < 1.

    # Number of monomials with minimum vanishing along divisors < 1
    number_dif_jn = sum(i < 1 for i in MonomioGen[:,number_functions_monomial + number_divisors])

    # Keep those elements in JumpingNumberList
    JumpingNumberList = MonomioGen[:number_dif_jn,:]

    # Number of different jumping numbers
    number_jn = len(set(JumpingNumberList[:,number_functions_monomial + number_divisors]))


    return (number_dif_jn,number_jn,JumpingNumberList)
    
    

''' Print the set of jumping numbers < 1, 
    together with the associated monomial and 
    the divisors at which they are attained. '''

def print_jumping_numbers(number_divisors,number_functions,number_functions_monomial,number_components,Matrix_orders_vanishing,\
                          number_jn_multiplicity,number_dif_jn,number_jn,JumpingNumberList):

    # The log-canonical threshold, 1st jumping number
    print("\n\nThe log-canonical threshold of the curve is ", end = '') 
    print(JumpingNumberList[0][number_functions_monomial + number_divisors], end = '')

    # Print the fractions corresponding to the lct
    divisors = np.where(JumpingNumberList[0][number_functions_monomial : number_functions_monomial + number_divisors] ==\
                        JumpingNumberList[0][number_functions_monomial + number_divisors])
    for supp in divisors[0]:
        jump_fraction = int(round(JumpingNumberList[0][number_functions_monomial + number_divisors] * \
                                  Matrix_orders_vanishing[supp][number_functions]))
        print(" = " + str(jump_fraction) + "/" + str(Matrix_orders_vanishing[supp][number_functions]), end = '')

    # Print divisors in the support of the lct
    print(" attained at divisors E_", divisors[0]+1)

    number_jn_diff_support = 1


    # Print rest of jumping numbers
    for i in range(1,number_dif_jn):
        # Support of current jumping number
        divisors = np.where(JumpingNumberList[i][number_functions_monomial : number_functions_monomial + number_divisors] == \
                            JumpingNumberList[i][number_functions_monomial + number_divisors])

        # If the jumping number is different from previous one
        if JumpingNumberList[i][number_functions_monomial + number_divisors] != \
        JumpingNumberList[i-1][number_functions_monomial + number_divisors]:
            print("\n")
            # Print the monomial
            for exp in range(number_functions_monomial):
                print("z_" + str(exp) + "^" + str(int(JumpingNumberList[i][exp])) + " ", end = '')

            # Print the jumping number
            print(" generates the jumping number " + \
                  str(JumpingNumberList[i][number_functions_monomial + number_divisors]), end = '')

            # Print the fractions corresponding to the jumping number and divisors in the support
            for supp in divisors[0]:
                jump_fraction = int(round(JumpingNumberList[i][number_functions_monomial + number_divisors] * \
                                          Matrix_orders_vanishing[supp][number_functions]))
                print(" = " + str(jump_fraction) + "/" + str(Matrix_orders_vanishing[supp][number_functions]), end = '')

            # Print the support of the jumping number (the set of divisors)
            print(" at divisors E_", divisors[0]+1)

            number_jn_diff_support += 1

        # If jumping number is equal to previous one but support is different    
        elif JumpingNumberList[i][number_functions_monomial + number_divisors] == \
        JumpingNumberList[i-1][number_functions_monomial + number_divisors] and\
        not np.array_equal(divisors,\
         np.where(JumpingNumberList[i-1][number_functions_monomial : number_functions_monomial + number_divisors] == \
                  JumpingNumberList[i-1][number_functions_monomial + number_divisors])):
            print("  ", end = '')
            # Print the monomial
            for exp in range(number_functions_monomial):
                print("z_" + str(exp) + "^" + str(int(JumpingNumberList[i][exp])) + " ", end = '')

            # Print the jumping number
            print(" also generates the jumping number " + \
                  str(JumpingNumberList[i][number_functions_monomial + number_divisors]), end = '')

            # Print the fractions corresponding to the jumping number and divisors in the support
            for supp in divisors[0]:
                jump_fraction = int(round(JumpingNumberList[i][number_functions_monomial + number_divisors] * \
                                          Matrix_orders_vanishing[supp][number_functions]))
                print(" = " + str(jump_fraction) + "/" + str(Matrix_orders_vanishing[supp][number_functions]), end = '')

            # Print the support of the jumping number (the set of divisors)
            print(" at divisors E_", divisors[0]+1, end = '')

            print(" (Different support).")
            number_jn_diff_support += 1



    print("\nThe number of monomial candidates < 1 was: " + str(number_dif_jn))
    print("The amount of different jumping numbers < 1 is: " + str(number_jn))
    print("The amount of jumping numbers (with different support) < 1 is: " + str(number_jn_diff_support))
    print("Amount of jumping numbers < 1 counted with multiplicity: " + str(number_jn_multiplicity))

    
(number_divisors,number_functions,number_functions_monomial,number_components,Matrix_orders_vanishing) = get_curve_data()
number_jn_multiplicity = get_Milnor(number_divisors,number_functions,number_functions_monomial,number_components,Matrix_orders_vanishing)
(number_dif_jn,number_jn,JumpingNumberList) = create_jumping_numbers(number_divisors,number_functions,number_functions_monomial,number_components,Matrix_orders_vanishing)
print_jumping_numbers(number_divisors,number_functions,number_functions_monomial,number_components,Matrix_orders_vanishing,\
                          number_jn_multiplicity,number_dif_jn,number_jn,JumpingNumberList)
    
    

# Example ThreeDiv
# 3, 5, 5, 2
# [[ 2 3 6 12 9 21 5] [ 4 6 15 30 18 48 13] [ 3 5 9 18 15 33 8]] 
# Example
# 4, 10, 4, 8
# [2 1 1 1 2 2 2 2 1 1 3] 
# [1 2 2 2 1 1 1 1 2 2 3] 
# [1 2 3 2 1 1 1 1 3 2 4] 
# [1 2 2 3 1 1 1 1 2 3 4]
# 3, 4, 1, 1, 1, 1, 3, 4