#!/usr/bin/env python3

import sys

# This script generates a configuration of 3-staircases in which there is a
# grid of particles in one phase, and the rest of space is filled with
# particles in another phase.

# size of grid
L1=100
L2=100

# grid of particles in phase 1
N11=12
N12=12

# position of lower left particle in phase 1
x01=30+N12
x02=30

# only print from here to there
l11=16
l12=16
l21=80
l22=88

# R2 of the model
R2=3

def main():
    # particles in phase 1
    particles1=[]
    # particles in phase 2
    particles2=[]

    fill_particles1(particles1)
    fill_particles2(particles2,particles1)

    # join particles
    particles=particles1+particles2

    # compute which particles have a Voronoi cells containing a given site
    voronoi_index=particles_in_voronoi(particles)
    # compute the Voronoi cells of each particle
    voronois=voronoi_of_particles(voronoi_index,particles)

    # check whether the particles are in the GFc
    GFc=particles_in_GFc(voronois,particles1,particles2)

    # boundary of GFc
    support=GFc_support(GFc,voronois)
    GFc_boundary=boundary(support)

    # print
    print_headers()
    print_particles(particles1,particles2,GFc)
    print_boundary(GFc_boundary)
    print_footers()

# generate grid of particles in phase 1
def fill_particles1(particles1):
    for k1 in range(N11):
        for k2 in range(N12):
            particles1.append((x01+2*k1-k2,x02+k1+3*k2,1))

# fill the rest of space with particles in phase 2
def fill_particles2(particles2, particles1):
    for x1 in range(L1):
        for x2 in range(L2):
            # index of particle
            k1=(x1+3*x2)/7
            k2=(2*x1-x2)/7

            # check the indices are integers:
            if k1 != int(k1) or k2 != int(k2):
                continue

            # check overlap with particles in phase 1
            if check_overlap_1((x1,x2),particles1):
                continue

            particles2.append((x1,x2,2))

# check overlap with particles in phase 1
def check_overlap_1(x,particles1):
    for p in particles1:
        if check_overlap(x,p):
            return True
    return False

# check overlap between two 3-staircases
def check_overlap(x,y):
    if max(abs(y[0]-x[0]),abs(y[1]-x[1]))>2:
        return False
    if abs((y[1]-x[1])+(y[0]-x[0]))<=2:
        return True
    return False


# given a site, compute the particles whose Voronoi cell the site is in
def particles_in_voronoi(particles):
    out=[]
    for x1 in range(L1):
        for x2 in range(L2):
            # minimize distance
            mindist=L1+L2+1
            candidate_list=[]
            for p in particles:
                dd=dist_support((x1,x2),p)
                if dd<mindist:
                    # reset candidate list
                    candidate_list=[p]
                    mindist=dd
                elif dd==mindist:
                    # append candidate
                    candidate_list.append(p)

            # add list to output (keep x1,x2 as index)
            out.append(((x1,x2),candidate_list))
    return out

# distance between a point and the support of a particle
def dist_support(x,p):
    # if x is in support then 0
    if abs(p[0]-x[0])+abs(p[1]-x[1])<=2 and x[0]>=p[0] and x[1]>=p[1]:
        return 0

    # base
    out=abs(x[0]-p[0])+abs(x[1]-p[1])
    # correct this depending on which site in \sigma_p is actually closest
    # lower left quadrant: p is the closest
    if x[0]-p[0]<=0 and x[1]-p[1]<=0:
        return out
    # lines next to lower left quadrant: neighbor of p closest
    if x[0]-p[0]==1 and x[1]-p[1]<=0 or x[1]-p[1]==1 and x[0]-p[0]<=0:
        return out-1
    # otherwise: second neightbor of p closest
    return out-2

# from a table associating sites to voronoi cells, get particles
def voronoi_of_particles(voronoi_index,particles):
    out=[]
    # output is indexed in the same order as particles
    for p in particles:
        vois=[]
        for v in voronoi_index:
            if p in v[1]:
                vois.append(v[0])
        out.append((p,vois))

    return out

# find which particles are in the GFc
def particles_in_GFc(voronois,particles1,particles2):
    out=[]
    # loop over particles
    for v in voronois:
        # check all other particles
        for w in voronois:
            # skip v
            if v==w:
                continue
            # check whether the distance between the cells is <=R2
            if set_dist(v[1],w[1])>R2:
                continue
            # check whether w is #-correct
            if not ((w[0][2]==1 and sharp_correct1(w[0],particles1)) or (w[0][2]==2 and sharp_correct2(w[0],particles2))):
                # in GFc
                out.append(v[0])
                break

    return out
                

# distance between two sets
def set_dist(s1,s2):
    if len(s1)==0 or len(s2)==0:
        print("error: computing the distance between two sets, one of which is empty",file=sys.stderr)
        exit(-1)
        return -1
    mind=abs(s1[0][0]-s2[0][0])+abs(s1[0][1]-s2[0][1])
    for x in s1:
        for y in s2:
            if abs(x[0]-y[0])+abs(x[1]-y[1])<mind:
                mind=abs(x[0]-y[0])+abs(x[1]-y[1])
    return mind

# check whether a particle of type 1 is #-correct:
def sharp_correct1(x,particles1):
    if (x[0]+2,x[1]+1,1) not in particles1:
        return False
    if (x[0]-2,x[1]-1,1) not in particles1:
        return False
    if (x[0]-3,x[1]+2,1) not in particles1:
        return False
    if (x[0]+3,x[1]-2,1) not in particles1:
        return False
    if (x[0]+1,x[1]-3,1) not in particles1:
        return False
    if (x[0]-1,x[1]+3,1) not in particles1:
        return False
    return True

# check whether a particle of type 2 is #-correct:
def sharp_correct2(x,particles2):
    if (x[0]+1,x[1]+2,2) not in particles2:
        return False
    if (x[0]-1,x[1]-2,2) not in particles2:
        return False
    if (x[0]+2,x[1]-3,2) not in particles2:
        return False
    if (x[0]-2,x[1]+3,2) not in particles2:
        return False
    if (x[0]-3,x[1]+1,2) not in particles2:
        return False
    if (x[0]+3,x[1]-1,2) not in particles2:
        return False
    return True


# support of a GFc
def GFc_support(GFc,voronois):
    out=[]
    for p in GFc:
        # find Voronoi cell
        for v in voronois:
            if p==v[0]:
                out=out+v[1]
                break

    # remove duplicates
    return list(set(out))

# boundary of a set
def boundary(S):
    out=[]
    for s in S:
        if (s[0]+1,s[1]) not in S:
            out.append(((s[0]+0.5,s[1]-0.5),(s[0]+0.5,s[1]+0.5)))
        if (s[0]-1,s[1]) not in S:
            out.append(((s[0]-0.5,s[1]-0.5),(s[0]-0.5,s[1]+0.5)))
        if (s[0],s[1]+1) not in S:
            out.append(((s[0]-0.5,s[1]+0.5),(s[0]+0.5,s[1]+0.5)))
        if (s[0],s[1]-1) not in S:
            out.append(((s[0]-0.5,s[1]-0.5),(s[0]+0.5,s[1]-0.5)))
    return out


# print particles
def print_headers():
    print(r'''\documentclass{standalone}

\usepackage{tikz}
\usepackage{shapes}

\begin{document}
\begin{tikzpicture}

''')
    
    # print grid
    print(r'\grid{'+str(l21-l11+1+3)+'}{'+str(l22-l12+1+3)+'}{('+str(l11-1)+','+str(l12-1)+')}')

def print_particles(particles1,particles2,GFc):
    for p in particles1:
        if p[0]>=l11 and p[0]<=l21 and p[1]>=l12 and p[1]<=l22:
            if p in GFc:
                print(r'\staircase{green}{('+str(p[0])+','+str(p[1])+')}')
            else:
                print(r'\staircase{cyan}{('+str(p[0])+','+str(p[1])+')}')
    for p in particles2:
        if p[0]>=l11 and p[0]<=l21 and p[1]>=l12 and p[1]<=l22:
            if p in GFc:
                print(r'\staircase{yellow}{('+str(p[0])+','+str(p[1])+')}')
            else:
                print(r'\staircase{red}{('+str(p[0])+','+str(p[1])+')}')

def print_boundary(boundary):
    for seg in boundary:
        if seg[0][0]>=l11 and seg[0][0]<=l21 and seg[0][1]>=l12 and seg[0][1]<=l22:
            if seg[1][0]>=l11 and seg[1][0]<=l21 and seg[1][1]>=l12 and seg[1][1]<=l22:
                print(r'\draw[line width=8pt]('+str(seg[0][0])+','+str(seg[0][1])+')--('+str(seg[1][0])+','+str(seg[1][1])+');')

def print_footers():
    print(r'''

\end{tikzpicture}
\end{document}''')

main()