# Given a file containing a point set and an edge-coloring of the 
# induced straight-line drawing of K_n, this script computes the
# minimum alpha that can be obtained over all matchings and details
# and derives the asymptotic crossing density in the doubling construction.
# Note that a matching in the input is ignored.

import sys
import networkx as nx

LEFT = 0
RIGHT = 1

# Functions from Appendix C
def a_alpha(x): return (x**2/24 - 3*x/28)
def b_0_0_alpha(x): return (x**2/16 - x/24)
def b_0_1_alpha(x): return (x**2/16 - x/48 + 1/336)
def b_1_0_alpha(x): return b_0_1_alpha(x)
def b_1_1_alpha(x): return (x**2/16 + 1/112)
def b_2_1_alpha(x): return (x**2/16 + x/48 + 3/112)
def c_0_0_alpha(x): return (x/12)
def c_0_1_alpha(x): return (x/12 + 1/168)
def c_1_0_alpha(x): return c_0_1_alpha(x)
def c_1_1_alpha(x): return (x/12 + 1/84)
def c_2_1_alpha(x): return (x/12 + 1/56)


# Returns whether vertices with indices i j k in this order form a counterclockwise triangle
def ccw(i,j,k):
    [ix,iy] = coordinates[i]
    [jx,jy] = coordinates[j]
    [kx, ky] = coordinates[k]
    alpha = (ky-iy) * (jx-ix) - (jy-iy) * (kx-ix)
    if alpha == 0:
        print("Collinear points", i, j, k)
        sys.exit()
    return alpha > 0

def edge_to_index(i,j):
    if i == j:
        print("There are no loops in K_n")
        sys.exit()
    a = min(i,j)
    b = max(i,j)
    return (N*(N-1) - (N-a-1)*(N-a)) // 2 + b - a - 1

# Computes the contribution to alpha from a single vertex p
# in the case where p_1 p_2 is not a matching edge of either p_1 nor p_2 and 
# the color of p_1 p_2 (c_bar) equals the color of the matching edge of p
def compute_alpha_one_equal_colors(c_bar, color_sides):
    left_bar = color_sides[c_bar][LEFT] 
    right_bar = color_sides[c_bar][RIGHT]
    value = 0
    for c in range(K):
        if c == c_bar:
            v_1, v_2 = 0, 0
            # We consider the offsets z^l_1, z^l_2, z^r_1, z_l^2 for the left and right side at the same time.
            # z^l_2 = z^r_2 = 1 independent of the choice of matching edge for p_2
            # Depending on the matching edge of p_1, either z^l_1 = 0 and z^r_1 = 2
            v_1 += 4 * (b_0_1_alpha(left_bar) + b_2_1_alpha(right_bar))
            v_1 += 2 * (c_0_1_alpha(left_bar) + c_2_1_alpha(right_bar))
            # or z^r_1 = 0 and z^l_1 = 2
            v_2 += 4 * (b_2_1_alpha(left_bar) + b_0_1_alpha(right_bar))
            v_2 += 2 * (c_2_1_alpha(left_bar) + c_0_1_alpha(right_bar))
            # We choose the better option
            value += min(v_1, v_2)
        else: # z_1, z_2 = 0 for all colors that are not c_bar
            value += 4 * b_0_0_alpha(color_sides[c][LEFT])
            value += 4 * b_0_0_alpha(color_sides[c][RIGHT])
    return value

# Computes the contribution to alpha from a single vertex p
# in the case where the new edge p_1 p_2 is not a matching edge of either p_1 nor p_2 and 
# the color of p_1 p_2 (c_bar) does not equal the color of the matching edge of p (c_prime)
def compute_alpha_one_different_colors(c_bar, c_prime, color_sides):
    left_bar = color_sides[c_bar][LEFT] 
    right_bar = color_sides[c_bar][RIGHT]
    left_prime = color_sides[c_prime][LEFT] 
    right_prime = color_sides[c_prime][RIGHT]
    value = 0

    for c in range(K):
        if c == c_bar or c == c_prime:
            pass
        else: # z_1, z_2 = 0 for all colors that are neither c_bar or c_prime
            value += 4 * b_0_0_alpha(color_sides[c][LEFT])
            value += 4 * b_0_0_alpha(color_sides[c][RIGHT])

    v_1, v_2, v_3, v_4 = 0, 0, 0, 0 # Try all four possible matchings for p_1 and p_2
    # Case p_1 and p_2 are matched to m_0(p)_r
    v_1 += 4 * (b_1_0_alpha(left_prime) + b_0_1_alpha(right_prime) + b_1_1_alpha(left_bar) + b_0_0_alpha(right_bar))
    v_1 += 2 * (c_1_1_alpha(left_bar)   + c_0_0_alpha(right_bar))
    # Case p_1 is matched to m_0(p)_r and p_2 is matched to m_0(p)_l
    v_2 += 4 * (b_1_1_alpha(left_prime) + b_0_0_alpha(right_prime) + b_1_0_alpha(left_bar) + b_0_1_alpha(right_bar))
    v_2 += 2 * (c_1_0_alpha(left_bar)   + c_0_1_alpha(right_bar))
    # Case p_1 is matched to m_0(p)_l and p_2 is matched to m_0(p)_r
    v_3 += 4 * (b_0_0_alpha(left_prime) + b_1_1_alpha(right_prime) + b_0_1_alpha(left_bar) + b_1_0_alpha(right_bar))
    v_3 += 2 * (c_0_1_alpha(left_bar)   + c_1_0_alpha(right_bar))
    # Case p_1 and p_2 are matched to m_0(p)_l
    v_4 += 4 * (b_0_1_alpha(left_prime) + b_1_0_alpha(right_prime) + b_0_0_alpha(left_bar) + b_1_1_alpha(right_bar))
    v_4 += 2 * (c_0_0_alpha(left_bar)   + c_1_1_alpha(right_bar))
    value += min(v_1, v_2, v_3, v_4)
    return value

# Computes the contribution to alpha from a single vertex p
# in the case where the new edge p_1 p_2 is the matching edge of p_1, so 
# the color of p_1 p_2 (c_bar) is equal to the color of the matching edge of p
def compute_alpha_two(c_bar, color_sides):
    left_bar = color_sides[c_bar][LEFT] 
    right_bar = color_sides[c_bar][RIGHT]
    value = 0
    for c in range(K):
        if c == c_bar:
            # The choice of matching edge of p_2 does not impact the z_1 and z_2.
            value += 4 * (b_1_1_alpha(left_bar) + b_1_1_alpha(right_bar))
            value += 2 * (c_1_1_alpha(left_bar) + c_1_1_alpha(right_bar))
        else: # z_1, z_2 = 0 for all colors that are neither c_bar or c_prime
            value += 4 * b_0_0_alpha(color_sides[c][LEFT])
            value += 4 * b_0_0_alpha(color_sides[c][RIGHT])
    return value


# Computes for each vertex q != p the contribution to alpha of p if p was matched to
# Output as list with entries (q, alpha).
def possible_halving_edges_and_values(p):
    candidates = []
    evaluation = []
    for q in range(N): # Compute the optimal alpha if the matching edge of p was edge pq
        if p == q: continue
        c_bar = coloring[p][q]
        color_sides = [[0,0] for _ in range(K)]
        for r in range(N): # Count colored edges on either side of pq
            if r == p or r == q: continue
            if ccw(p,q,r):
                color_sides[coloring[p][r]][LEFT] += 1
            else:
                color_sides[coloring[p][r]][RIGHT] += 1
        total_colors = [c[LEFT] + c[RIGHT] for c in color_sides]
        total_colors[c_bar] += 1

        # Compute the minimum contribution to alpha over all possible choices of details at p

        value = compute_alpha_two(c_bar, color_sides) # Case: p_1 is matched to p_2
        for c_prime in range(K): # Case: p_1 p_2 gets color c and p_1,p_2 are matched to q_1, q_2
            if c_prime == c_bar: # Subcase: the color of pq coincides with the color of p_1 ip_2
                value = min(value, compute_alpha_one_equal_colors(c_bar, color_sides))
            else:
                value = min(value, compute_alpha_one_different_colors(c_bar, c_prime, color_sides))

        candidates += [q]
        evaluation += [value]
    return zip(candidates, evaluation)

    
assert len(sys.argv) >= 2, "Requires file path of a drawing"

# Extract elements of filename
initial_crossings = int(sys.argv[1].split(".")[-2].split("_")[-1])
file_extension = sys.argv[1].split(".")[-1]
K = int(file_extension[2:])

# Read input file
with open(sys.argv[1], "r") as input:
    N = int(input.readline().split(" ")[0])
    M = int(N * (N-1)/2)

    assert file_extension.startswith("co"), "File does not contain point set"

    coordinates = [list(map(int, input.readline().split(" "))) for i in range(N)]
    coloring_string = input.readline()

# Extract colors of edges and
coloring = [[-1 for _ in range(N)] for _ in range (N)]
counter = 0
for p in range(N):
    for q in range(p+1, N):
        col = int(coloring_string[counter])
        coloring[p][q] = col
        coloring[q][p] = col
        counter = counter + 1

# Solve matching problem using external library
bip_graph = nx.Graph()
bip_graph.add_nodes_from(range(N), bipartite = 0)
bip_graph.add_nodes_from(range(N, M+N), bipartite = 1)
for p in range(N):
    # Set edge weights to minimal local contributions to alpha
    bip_graph.add_weighted_edges_from([(p, N + edge_to_index(p,q), x) for (q,x) in possible_halving_edges_and_values(p)])
bipartite_matching = nx.algorithms.bipartite.minimum_weight_full_matching(bip_graph)

# Extract contribution of the vertices to alpha when using the optimal matching
vertex_contributions = sum([bip_graph.edges[p, bipartite_matching[p]]["weight"] for p in range(N)])

total_alpha = initial_crossings + a_alpha(N) + vertex_contributions
asymptotic_density = total_alpha * 24/(N**4)
print("Upper bound on geometric " + str(K) + "-colored crossing constant:", asymptotic_density)