##/////////////////////////////////////////////////////////////
##Copyright 2017 Michele Zito and Marco Bertamini
##
##Occupancy.py is free software; you can redistribute it and/or modify
##it under the terms of the GNU General Public License as published by the
##Free Software Foundation.
##
##Occupancy.py is distributed in the hope that it will be useful,
##but WITHOUT ANY WARRANTY; without even the implied warranty of
##MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
##GNU General Public License for more details.
##
##For copy of the GNU General Public License see <http://www.gnu.org/licenses/>.
##
##Please retain this note in the code. If you use the program for academic research
##The most relevant published citation is the following:
##Bertamini, M. Zito, M. Scott-Samuel, N.E. & Hulleman, J. (2016). Spatial Clustering and Its Effect on Perceived Clustering, 
##Numerosity and Dispersion Attention, Perception & Psychophysics, 78(5), 1460-1471.
##/////////////////////////////////////////////////////////////////////

##/////////////////////////////////////////////////////////////////////
##README
##
##The script is called occupancy.py
##
##You do not need to edit the script. Prepare a text file in which
##the first line is the radius of the area of influence.
##
##After the first line, each line has one pair of x,y coordinates separated by a comma. 
##An example file is provided called coords.txt
##
##The units are the same for the radius and for the coordinates.
##For example pixels or mm. It does not matter how the coordinates have been 
##generated.
##
##The occupancy area computed is the area of all regions of influence
##taking into account the overlap. This value is written in the output and
##also saved in a file called result.txt. In the file the first column is 
##a timestamp, the second is the occupancy area. 
##
##The program has been tested on a Mac, not on Windows.
##
##Known issues:
##Remove any blank lines at the end of your text file.
##For number of elements above 50 the program is slow and 
##exponentially more so as the number increases.
##
##For copyright and reference see the top of the script.
##
##For information contact:
##Marco Bertamini <m.bertamini@liv.ac.uk>
##Michele Zito <michele@liverpool.ac.uk>

from numpy import * #many different maths functions
import os #handy system and path functions
from math import *
from itertools import *
import sys, copy, datetime


#### EXTRA
#euclidian distance
# usage d = distance([0,0],[10,10])
def d2(p1,p2):
    return sqrt((p2[0] - p1[0]) ** 2 +(p2[1] - p1[1]) ** 2)

# area of any convex polygon, starting from the vertices (see it.wikipedia.org/wiki/Poligono)
# usage area = gaussarea(p) where p is a list of lists: [[0,0],[2,3],[30,23]]
# 
# VERY IMPORTANT!! The vertices of the polygon must be given in the order in which they are met
# during a walk around the polygon edges (either clockwise or anticlockwise). Any other ordering 
# is bound to produce meaningless results. [http://geomalgorithms.com/a01-_area.html]
# 
def gaussarea(poly):

    temp1 = 0
    temp2 = 0
    poly2 = copy.deepcopy(poly)
    poly1 = copy.deepcopy(poly)
    poly1.append(poly[0])
    
    for i in range(len(poly2)):
        temp1 = temp1 + poly1[i][0]*poly1[i+1][1]

    for i in range(len(poly2)):
        temp2 = temp2 + poly1[i+1][0]*poly1[i][1]
       
    return (0.5 * abs(temp1 - temp2))
    

# return the xy coordinate of the centoid of a set of points
# usage area = analysecentroid(p) where p is a list of lists: [[0,0],[2,3],[30,23]]
def analysecentroid(matrix):
    sxy =0.0
    sx =0.0
    sy =0.0
    svalue=0
    centroid=[0,0]
           
    for i in range(len(matrix)):
        sx = sx + matrix[i][0]
        sy = sy + matrix[i][1]
        
    svalue = len(matrix)
    
    centroid[0] = sx/svalue
    centroid[1] = sy/svalue
    
    return centroid



#### CONVEX HULL

# This code returns the vertices of a convex hull for a set of arbitrary vertices
# the vertices of the convex hull are returned ordered
# 
## [MZ, 2014/10/15] If we know that points are indeed the points in the convex hull and they just
## need to be reallanged then a simpler computation suffices.
##
TURN_LEFT, TURN_RIGHT, TURN_NONE = (1, -1, 0)
 
def turn(p, q, r):
    """Returns -1, 0, 1 if p,q,r forms a right, straight, or left turn."""
    return cmp((q[0] - p[0])*(r[1] - p[1]) - (r[0] - p[0])*(q[1] - p[1]), 0)
 


def clockOrder(points):
    if len(points) > 2:
        p = 0
        q = 1
        hull = [points[p],points[q]]
        prevTurn = -2
        repeat = True
        while repeat:
            for k in range(2,len(points)+1):
                kk = k if k<len(points) else 0
                if kk!=2 and prevTurn != turn(points[p],points[q],points[kk]):
                    # swap q and points[kk]
                    t = points[q]
                    points[q] =  points[kk]
                    points[kk] = t
                    # start again
                    p = 0
                    q = 1
                    hull = [points[p],points[q]]
                    repeat = True
                    break;
                else:
                    if prevTurn == -2:
                        prevTurn = turn(points[p],points[q],points[kk])
                    # carry on. update hull and shift the points along
                    if k<len(points):
                        hull.append(points[k])
                        p += 1
                        q += 1
                if k == len(points):
                    repeat = False
    else:
        hull = points
    return hull





# BETTER BOUND ON CLIQUE SIZE: returnableDeg is the smallest positive degree for 
# which there are at least returnableDeg+1 vertices of degree at least returnableDeg
# (ie, there's enough vertices around to build a clique)
#
# a = list of dots
# r = radius
#
def refinedBound(a, r):
    #initializations
    neighbors =[[] for i in range(len(a))]  # adjacency lists
    degrees =[0 for i in range(len(a))]  # degree sequence: degrees[k] is the number of vertices of degree k
    returnableDeg = 0

    
    ##     FILL THE GRAPH
    # fill the adjacenct lists
    for index1, i in enumerate(a):
        # here  index  is the order number of i in a
        for index2, j in enumerate(a):
            # here i and j are points in a
            if d2(i, j) <2*r and index1 != index2:
                neighbors[index1].append(j)
        kk = len(neighbors[index1])
        degrees[kk]+=1

    ##     COMPUTE THE BOUND
    sum = 0
    for i in range(len(degrees)-1,-1,-1):
        sum += degrees[i]
        if sum>i and degrees[i]>0:
            ##     PICKS THE FIRST i FOR WHICH
            ##     degree[i]>0   there is a vertex of degree i
            ##     AND
            ##     sum>i         there is enough vertices of degree at least i to build a K_{i+1}
            ##    
            ##     Function returns i+1 as the first NON-FEASIBLE VALUE.
            return neighbors, i+1

    ##     DEFAULT RETURN (returnableDeg is always zero) 
    return neighbors, returnableDeg
#end refined Bound


# compute the coordinates that intersect two circles with centres c1 and c2 and radius r (fixed)
def cIntersect2(centre0, radius0, centre1, radius1):
#    This function checks for the intersection of two circles.
#    If one circle is wholly contained within the other c2cintersect is -1
#    If there is no intersection of the two circles c2cintersect is 0 
#    If the circles intersect c2cintersect is 1 and the function return xi1, yi1, xi2, yi2


    x0 = centre0[0]
    y0 = centre0[1]
    x1 =centre1[0]
    y1 =centre1[1]
    r0 = radius0
    r1 = radius1

    #dx and dy are the vertical and horizontal distances between the circle centers.
    dx = x1 - x0
    dy = y1 - y0

    #Determine the straight-Line distance between the centers.
    d = sqrt((dy*dy) + (dx*dx))
 
    #Check for solvability.
    if (d > (r0 + r1)):          #no solution. circles do Not intersect
        c2cIntersect = 0
        xi1 = yi1 = xi2 = yi2 =[]
    elif (d < abs(r0 - r1)):   #  no solution. one circle is contained in the other
        c2cIntersect = -1
        xi1 = yi1 = xi2 = yi2 =[]
    else:
    #    Determine the distance from point 0 To point 2. point 2' is the point where the Line through the circle
    #    intersection points crosses the Line between the circle centers.
        a = ((r0*r0) - (r1*r1) + (d*d)) / (2.0 * d)

    #    ' Determine the coordinates of point 2.
        x2 = x0 + (dx * a/d)
        y2 = y0 + (dy * a/d)

    #    Determine the distance from point 2 To either of the intersection points.
        h = sqrt((r0*r0) - (a*a))

         #     Now determine the offsets of the intersection points from point 2.
        rx = (0-dy) * (h/d)
        ry = dx * (h/d)

    # Determine the absolute intersection points.
        xi1 = x2 + rx
        xi2 = x2 - rx
        yi1 = y2 + ry
        yi2 = y2 - ry
        
        c2cIntersect = 1

    return xi1, yi1, xi2, yi2


def overlap(c, r):
    # Compute the set P of points P1 = (x1,y1),...,Pk = (xk,yk), belonging to the intersection of two circles and to all others. 
    p =[]

    q1 =[[0,0],[0,0]]    
    q2=[[0,0],[0,0]]
    q3=[[0,0],[0,0]]
    q4=[[0,0],[0,0]]
    q5=[[0,0],[0,0]]
    q6=[[0,0],[0,0]]

    k = len(c)

    # Compute the set of points belonging to the intersection of two circles and to all others
    for i in list(combinations(range(0, k), 2)):
        q1=[[0,0],[0,0]]    
        q2=[[0,0],[0,0]]
        # Find Q1 and Q2 the intersection between Ci0 and Ci1
        q1[0], q1[1], q2[0], q2[1] =cIntersect2(c[i[0]], r, c[i[1]], r)
        if q1 == [[],[]] or q2 == [[],[]]:
            return 0
        
        # variable q stores the point in Ci0 cap Ci1 that belongs to all other Cj's
        q =[]
        for h in range(0, k):    # Low dimension cases 
            if h != i[0] and h != i[1]:
                dq1c = d2(q1,c[h])
                dq2c = d2(q2,c[h])
                if (dq1c<r) and (dq2c <r):   # both inside
                    newc = [x for x in c if x != c[h]]
                    return overlap(newc, r)
                elif (dq1c >= r) and (dq2c >= r):
                    p = []
                    break
                elif q==[]:
                    q=q1 if (dq1c<r) else q2
                else:
                    qq=q1 if (dq1c<r) else q2
                    
                    if (q!=qq):
                        q = []
                        break
        
        if (k==2):
            p.append(q1)
            p.append(q2)
        elif (q!=[]):
            p.append(q)
    
    # Compute the area T of the polygon on P1, . . . , Pk using Gauss formula: 
    if k>2 and len(p)>2:
       p = clockOrder(p)
       T = gaussarea(p)
    else:
       T=0


    if (k>2) and len(p)>2:
        centroid =analysecentroid(p)
        aa = array(p)   # same as p but as a numpy array
        aaa = aa - centroid     # subtract the centroid
        ppp = aaa[argsort(arctan2(aaa[:,1], aaa[:,0]))]     #sort CW
    else:
        ppp = p;
        
    #Compute the areas S1, . . . , Sk of the k circular segments on P1,P2, P2,P3, ..., and Pk,P1  
    esse = 0
    for index in range(len(ppp)):
        first = index
        second = index +1
        if index ==len(ppp) -1:
            second =0
        a1 = d2(ppp[first], ppp[second])
        tetha1 = 2. * asin(a1/(2.*r))
        s1 = r*r * (tetha1 - sin(tetha1))/2 
        esse = esse + s1
        
    return T + esse

def occupancy(a,r):
#    print "Occupancy(", a, ",", r,")"
    area1 = pi* r * r    # the area of a full circle
    diam = 2*r    # diameter
    area = area1 * len(a)  # total area without
    
#    print "Starting Area: \t\t", area
    overarea = 0
    parity = -1
    space = " "

    # refinedBound based on groups of vcs of max deg.
    # This time DD is the largest size of a maximum clique as for any larger value there is
    # at most d vertices of degree d (and a complete graph of degree d is K_{d+1}
    neighbors, DD = refinedBound(a,r)
#    print "Maximum number of intersecting circles: %2d\n      " % (DD),

    # [MZ, 27/11/2014]
    # FIRST OPTIMIZATION: the main cycle looks at intersections between k circles (k-intersections).
    # If a vertex is not part of any (k-1)-intersection it definitely cannot be part of a k-intersection
    # (in fact it needs to be part of at least k-1 (k-1)-intersections, we are not checking this though
    # and it is questionable if this woudl actually give any big saving (we would need counters for
    # all elements that appear in a (k-1)-intersection: only those appearing in at least k-1 different
    # intersections will end up in groundset next time around)!)
    #
    # 
    # In each iteration groundset contains reference to the set of circles that belong to at least one
    # (k-1)-intersection.
    # 
    groundset = set(range(0,len(a)))
        
    for k in range(2, DD+1):
        for i in range(1,k):
            space = space + "="
        loop = 0

        # needed to update groundset
        newGroundset = set()
        
        cliques = 0
        for i in combinations(groundset, k):
            Good = True
            j = []

            # [MZ, 27/11/2014]
            # FIRST FILTERING: if the degree of each of the k-tuple elements is not at least
            # k-1 there is no change they will all intersect so the joint interesection must
            # be empty (tested in O(k) steps)
            for jj in range(0,k):
                if len(neighbors[i[jj]]) < k-1:
                    Good = False
                    break
                else:
                    j.append(a[i[jj]])
            if Good:
                # [MZ, 27/11/2014]
                # SECOND FILTERING: clique checking, if k vertices have passed the "degree test"
                # we need to check that they form a k-clique, or else their **joint** intersection
                # will be empty (tested in O(k^2) steps)
                for jj in range(0,k):
                    for kk in range(0,jj):
                        if d2(j[jj],j[kk])>=2*r:
                            Good = False
                            break
                if Good:
                    temp = overlap(j, r)
                    if temp>0:
                        cliques+=1
                        
                        for jj in range(0,k): newGroundset.add(i[jj])

                        overarea = overarea + parity * temp
#                        if (parity == 1):
#                            print space, '  +', j, '\t', temp
#                        else:
#                            print space, '  -', j, '\t', temp
        groundset = newGroundset
        parity = parity * -1

        # [MZ, 27/11/2014]
        # SECOND OPTIMIZATION: if there is less than k-1 cliques leading to non-empty
        # intersection we can break the whole computation, without the need to go up to
        # the maximum degree DD.

        if cliques < k+1: break
#    print "\b\b\b\b\b\b",
        
    finalarea = area + overarea
    return finalarea

# end OCCUPANCY



#### TEST CASES:

#a = [ [ -66.15480776, 282.8991993],  [-209.06454318, -10.18274172],  [  60.66724238, 236.43864335]]
#d = 100
#print occupancy(a,d)
#
#a = [ [-104.72204258, 207.51844101],  [ 247.23932665,  -50.34252937],  [ 171.52499851, -115.89956405]] 
#print occupancy(a,d)
# 81990.3695666
#
#a = [[32.962280490376813, 263.28537883335321], [146.1166459578493, 46.855807912736168], [60.819644919131711, 274.76062738498462]]
#print occupancy(a,d)
# 68834.643707
#
#a = [[12,240],[0,0],[-1,-244]]
#print occupancy(a,d)
# 94247.7796077
#
#a = [[-109.76498786290672, -109.5252985846228], [78.211475775509754, -191.19216664600384], [-1.4888916703626593, -14.787963349172513], [141.55342151677559, 124.95656006682363], [83.921273534137043, -122.00046857637021]]
#print occupancy(a,d)
# 127463.288765
#
#a = [[103, 2], [-65, -85], [53, -113], [58, -36], [-44, 133], [7, 35], [-58, 100], [95, -72], [36, 142], [-129, -80], [-108, 78], [-110, -28], [135, -24], [111, 73], [87, 88], [-87, 25], [-64, -117], [-84, 94], [-22, 136], [-23, 86], [-74, -141], [-104, 8], [-79, -99], [-35, -63], [69, 58], [64, -144], [151, -37], [-11, -117]]
#print occupancy(a,d)
# 126325.471429

#from here the next few lines simply read a string and turn it into a list of lists
fin = open('coords.txt')
fout = open('result.txt', "a")

a =[]
with open('coords.txt', 'r') as f:
    regionInfluence = int(f.readline())
    for line in f:
        geta = "".join(line.split())
        c = map(int, geta.split(","))
        a.append(c)

occupancyResult =occupancy(a, regionInfluence)
print occupancyResult

fout.write(str(datetime.datetime.utcnow()))
fout.write('\t')
fout.write(str(occupancyResult))
fout.write('\n')
    
fin.close()
fout.close()