#! /usr/bin/python

###############################
# Author : Alexandre Mouradian
# year : 2013
###############################

#input : - infile.xml -> node TA
#		 - graph -> description of the topology

#output : outfile.xml -> NTA + query.q

#graph example :
#	
#p 0 48.1540811791 78.4171121718
#p 1 22.3034450305 86.8491780822
#p 2 73.2186900176 90.7916104965
#p 3 41.5013633282 49.1759685779
#p 4 46.5326056803 78.3737165383
#0 1 2 3 4 
#1 0 4 
#2 0 4 
#3 0 4 
#4 0 1 2 3 

import sys
import os
import copy
import subprocess
import time

import numpy as np
import math
try:
	from lxml import etree as ElementTree
except ImportError, e:
	try:
		from elementtree import ElementTree
	except ImportError, e:
		print '***'
		print '*** Error: Must install either ElementTree or lxml.'
		print '***'
		raise ImportError, 'must install either ElementTree or lxml'

class Node:
	def __init__(self):
		self.neighbors=[]
		self.myCliques=[]
	ID=0
	rank=100
	coordinate=0
	mts=1
	x=0.0
	y=0.0

class Data:
	A=0
	B=0
	C=0

#global variables
V=0
graph=[]
graph2=[]
boolGraph=[]
boolGraph2=[]
link=[]
data=[]
dataBorder=[]
nMax=10
indMax=9
STEP=1.0
T=[]
probabilities=[]
T2check=[]
cg=[]
Cap=99

def parseGraph(graphFile):#retreiving the graph from the file (parse 2 times)
	ifi=open(graphFile,'r')
	l=ifi.readline()
	i=0
	while l!="":#creates nodes
		ls=l.split()
		if ls[0]!="p":
			graph.append(Node())
			graph[i].ID=int(ls[0])
			i+=1
		l=ifi.readline()
	ifi=open(graphFile,'r')
	l=ifi.readline()
	
	i=0
	while l!="":#assigns neighbors to nodes
		ls=l.split()
		if ls[0]!="p":
			for j in range(1,len(ls)):
				for k in range(len(graph)):
					if graph[k].ID==int(ls[j]):
						graph[i].neighbors.append(graph[k])
			i+=1
		l=ifi.readline()

	ifi=open(graphFile,'r')
	l=ifi.readline()

	i=0
	while l!="":#assigns locations to nodes
		ls=l.split()
		if ls[0]=="p":
			graph[i].x=float(ls[2])
			graph[i].y=float(ls[3])
			i+=1
		l=ifi.readline()


def setRank(node):#recursive algorithm to set ranks, call with the sink as parameter
	for i in range(len(node.neighbors)):
		if node.neighbors[i].rank>node.rank+1:
			node.neighbors[i].rank=node.rank+1
			setRank(node.neighbors[i])

def F(theta, cn, ci, R):
	inter=0.0
	ret=0.0
	inter=math.sqrt(pow(R,2)-pow((cn-1)*R+ci,2)*pow(math.cos(theta),2))
	ret=(((cn-1)*R+ci)/2)*(-math.cos(theta)*inter-(pow(R,2)*math.atan(((cn-1)*R+ci)*math.cos(theta)/inter))/((cn-1)*R+ci));
	return ret

#generate data for RTXP coordinates
def initData():
	for i in range(nMax+1):
		c=[]
		d=[]
		for j in range(indMax+1):
			c.append(Data())
			d.append(Data())
		data.append(c)
		dataBorder.append(d)


	xc1=0.0
	yc1=0.0
	xc2=0.0
	yc2=0.0
	i=0.0
	A=0.0
	C=0.0
	B=0.0
	pA=0.0
	pC=0.0
	pB=0.0
	R=10.0
	n=1
	
	while n<nMax:
		i=0.0
		while i<R:
			if not (n==1 and i==0.0):
				if n!=1:
					yc1=(pow((n-1)*R,2)-pow(R,2)+pow((n-1)*R+i,2))/(2*((n-1)*R+i))
					xc1=np.sqrt(pow((n-1)*R,2)-pow(yc1,2))
					#deduces the angles for A
					theta1=math.atan(yc1/xc1)
					theta2=math.atan(-yc1/xc1)+math.pi

				#computes the intersection points for zone C
				yc2=(pow(n*R,2)-pow(R,2)+pow((n-1)*R+i,2))/(2*((n-1)*R+i))
				xc2=np.sqrt(pow(n*R,2)-pow(yc2,2))

				#deduces the angles for C
				theta1p=math.atan(yc2/xc2)
				theta2p=math.atan(-yc2/xc2)+math.pi

				#comoputes the actual surface of A	and of C depending on the current i which is the "offset"	
				if n!=1:
					A=n*(n-2)*pow(R,2)*(theta2-theta1)/2+(pow((n-1)*R+i,2)*(math.sin(theta2)*math.cos(theta2)-math.sin(theta1)*math.cos(theta1))/2)+F(theta2,n,i,R)-F(theta1,n,i,R)
				C=(1-pow(n,2))*pow(R,2)*(theta2p-theta1p)/2-(pow((n-1)*R+i,2)*(math.sin(theta2p)*math.cos(theta2p)-math.sin(theta1p)*math.cos(theta1p))/2)+F(theta2p,n,i,R)-F(theta1p,n,i,R)
			
			#transforms i in an integer in order to store computed value with i and n as indexes
			ind=int(i/STEP)

			if n>1:
				#general case
				pA=A/(math.pi*pow(R,2))
				pC=C/(math.pi*pow(R,2))
				pB=1-pA-pC
				data[n][ind].A=pA
				data[n][ind].B=pB
				data[n][ind].C=pC
				#border data
				pB=1-pA;
				dataBorder[n][ind].A=pA
				dataBorder[n][ind].B=pB
				dataBorder[n][ind].C=0
			else:#n=1 => there is no surface A
				if n==1 and i==0.0:
					pB=1.0
					pC=0.0
				else:
					pC=C/(math.pi*pow(R,2))
					pB=1-pC
				data[n][ind].A=0
				data[n][ind].B=pB
				data[n][ind].C=pC

			i+=STEP
		n+=1

#assign RTXP coordinates
def setCoordinate():
	for i in range(len(graph)):
		if graph[i].rank==0:
			graph[i].coordinate=1
		else:
			up=0
			lev=0
			down=0
			for j in range(len(graph[i].neighbors)):
				if graph[i].neighbors[j].rank < graph[i].rank:
					down+=1
				elif graph[i].neighbors[j].rank == graph[i].rank:
					lev+=1
				elif graph[i].neighbors[j].rank > graph[i].rank:
					up+=1
			tot=up+lev+down
			pup=float(up)/float(tot)
			plev=float(lev)/float(tot)
			pdown=float(down)/float(tot)
			r=graph[i].rank
			mini=100.0
			res=0
			current=0.0

			for k in range(len(data[r])):
				if pup>0.0:
					current=math.sqrt(pow(((data[r][k].A)-pdown),2)+pow(((data[r][k].B)-plev),2)+pow(((data[r][k].C)-pup),2))
				else:
					#Border case
					current=math.sqrt(pow(((dataBorder[r][k].A)-pdown),2)+pow(((dataBorder[r][k].B)-plev),2))

				if current < mini:
					mini=current
					res=float(k*STEP)


			graph[i].coordinate=int((res+mini)*10.0)
			

def neighb2():#to be used on a copy of the graph: it is the 2hop connection graph
	for i in range(len(graph)):
		for j in range(len(graph[i].neighbors)):
			for k in range(len(graph)):
				if graph[k].ID==graph[i].neighbors[j].ID:
					for l in range(len(graph[k].neighbors)):
						if graph[k].neighbors[l].ID!=graph[i].ID:
							graph2[i].neighbors.append(copy.copy(graph[k].neighbors[l]))

def neighbool2(s):
	s2=[]
	for i in range(len(s)):
		s2.append(copy.copy(s[i]))
	for i in range(len(s)):
		for j in range(len(s)):
			if s[i][j]==1:
				for k in range(len(s)):
					if s[j][k]==1 and k!=i:
						s2[i][k]=1
	return s2

def boolToStruct(b):
	g=[]
	for i in range(len(b)):
		g.append(Node())
		g[i].ID=i
	for i in range(len(b)):
		for j in range(i,len(b)):
			if b[i][j]==1 and i!=j:
				g[i].neighbors.append(g[j])
				g[j].neighbors.append(g[i])
	return g


def verifRanks():#check if node ranks are up to date,if not remove the topology from "T2check"
	for i in range(len(T)):
		T2check.append(i)
	for i in range(len(T)):
		cg=[]
		cg=boolToStruct(T[i])
		cg[0].rank=0
		setRank(cg[0])
		for j in range(len(cg)):
			if cg[j].rank>graph[j].rank:
				T2check.remove(i)
				break

def prob(d):
	Gt=1.0
	Gr=1.0
	ple=3.0
	Nb=800
	L=1
	f=2.4*10**9
	c=3*10**8
	lambd=c/f
	N0=10**(-17.4)*10**-3
	B=0.5*10**6
	alpha=1
	beta=2

	#***************
	Pt=1*10**-3
	#***************

	K=Gt*Gr*lambd**2/(((4*math.pi)**2)*L)
	SNR=K*Pt*(d**-ple)/(N0*B)
	BER=alpha/(2*beta*SNR)
	pl=(1-BER)**Nb
	return pl

def distance(x1,y1,x2,y2):
	return math.sqrt((x1-x2)**2+(y1-y2)**2)

#generate topologies
def initT(mat):
	current=[]
	for i in range(len(mat)):
		current.append(mat[i][:])
	for i in range(len(graph)):
		for j in range(i,len(graph)):
			if current[i][j]==1 and i!=j:
				current[i][j]=0
				current[j][i]=0
				found=0
				for k in range(len(T)):
					if T[k]==current:
						found=1
				if found!=1:
					T.append(current)
					initT(current)
				icurrent=[]
				for q in range(len(current)):
					icurrent.append(current[q][:])
				current=icurrent
				current[i][j]=1
				current[j][i]=1
	return 0


#set link probabilities : i
def linkProb():
	lim=1
	for i in range(len(boolGraph)):
		for j in range(len(boolGraph)):
			if boolGraph[i][j]==1:
				lim=5#number of retransmissions
				d=distance(graph[i].x,graph[i].y,graph[j].x,graph[j].y)
				pr=prob(d)
				link[i][j]=1-(1-pr)**(lim+1)

#add probabilities to topologies
def trans2():
	linkProb()
	for j in range (len(T)):
		p=-1
		for k in range(len(graph)):
			for l in range(k,len(graph)):
				if k!=l:
					if boolGraph[k][l]==1:
						if T[j][k][l]==1:
							if p==-1:
								p=link[k][l]
							else:
								p=p*link[k][l]
						if T[j][k][l]==0:
							if p==-1:
								p=(1-link[k][l])
							else:
								p=p*(1-link[k][l])
		probabilities.append(p)

def setAutomata(inFileName, outFileName):
	V=len(graph)
	sysString1=""
	sysString2="system"
	docO = ElementTree.parse(inFileName)#the new templates will be added to this tree (O is for ouput)
	ntaO = docO.getroot()
	for i in range(V):
		if graph[i].mts==1:
			sysString1=sysString1+" n"+str(i)+" = "+"Node("+str(graph[i].rank)+","+str(graph[i].coordinate)+",1,"+str(i)+")"+";\n"#strings which create the system objects, case sender
		else:
			sysString1=sysString1+" n"+str(i)+" = "+"Node("+str(graph[i].rank)+","+str(graph[i].coordinate)+",0,"+str(i)+")"+";\n"#strings which create the system objects, case not sender
		if i==0:
			sysString2=sysString2+" n"+str(i) #define the object which are part of the system
		else:
			sysString2=sysString2+", n"+str(i) 
	sysString2=sysString2+";"
	system=ntaO.find("system")
	system.text=sysString1+sysString2#put the text in the xml file
	
	declaration=ntaO.find("declaration")# this part sets the general declarations
	decString="const int N="+str(V)+";\n const bool connect1[N][N]=\n{"	
	for i in range(len(boolGraph)):
		decString=decString+"{"
		for j in range(len(boolGraph)):
			decString=decString+str(boolGraph[i][j])
			if j<len(boolGraph)-1:
				decString=decString+","
		if i<len(boolGraph)-1:
			decString=decString+"},\n"
		else:
			decString=decString+"}"
	decString=decString+"};\n\n"

	decString=decString+"const bool connect2[N][N]=\n{"	
	for i in range(len(boolGraph2)):
		decString=decString+"{"
		for j in range(len(boolGraph2)):
			decString=decString+str(boolGraph2[i][j])
			if j<len(boolGraph2)-1:
				decString=decString+","
		if i<len(boolGraph2)-1:
			decString=decString+"},\n"
		else:
			decString=decString+"}"
	decString=decString+"};\n\n"
	
	declaration.text=decString+declaration.text #put the declaration string in the xml

	docO.write(outFileName)#write the tree in the output file

def findDuplicate():
	for i in range(len(graph)):
		for j in range(i,len(graph)):
			if graph[i].rank==graph[j].rank and graph[i].coordinate==graph[j].coordinate and i!=j:
				for k in range(len(graph[i].neighbors)):
					if graph[j].ID==graph[i].neighbors[k].ID:
						return [i,j]
	return [0,0]

def findDuplicate():
	for i in range(len(graph)):
		for j in range(len(graph[i].neighbors)):
			for k in range(len(graph[i].neighbors)):
				if graph[i].neighbors[j].rank==graph[i].neighbors[k].rank and graph[i].neighbors[j].coordinate==graph[i].neighbors[k].coordinate and j!=k:
					return [graph[i].neighbors[j].ID,graph[i].neighbors[k].ID]
	return [0,0]

def main():
	global boolGraph2
	global boolGraph
	global link
	args = sys.argv[1:]
	if len(args) != 3:
		print 'usage: generateCases.py infile.xml outfile.xml graphFile'
		sys.exit(-1)
	parseGraph(sys.argv[3])
	ofiQueries=open("rtxp.q",'w')

	for i in range(len(graph)):#allocate list
		boolGraph.append(len(graph)*[0])
		link.append(len(graph)*[0.0])

	NBlink=0
	for i in range(len(graph)):#produce a matrix of booleans which represent the graph
		for j in range(len(graph[i].neighbors)):
			boolGraph[graph[i].ID][graph[i].neighbors[j].ID]=1
			NBlink+=1
	
	for i in range(len(graph)):
		graph2.append(Node())
		graph2[i].ID=graph[i].ID

	for i in range(len(graph)):
		for j in range(len(graph[i].neighbors)):
			graph2[i].neighbors.append(graph2[graph[i].neighbors[j].ID])

	neighb2()

	for i in range(len(graph2)):#allocate list
		boolGraph2.append(len(graph2)*[0])

	for i in range(len(graph2)):#produce a matrix of booleans which represent the 2hop graph
		for j in range(len(graph2[i].neighbors)):
			boolGraph2[graph2[i].ID][graph2[i].neighbors[j].ID]=1

	V=len(graph)

	NBlink=NBlink/2
	NBT=2**NBlink
	T.append(copy.copy(boolGraph))
	initT(boolGraph)

	graph[0].rank=0
	setRank(graph[0])
	trans2()
	print link
	print probabilities

	sumline=0
	for j in range(len(probabilities)):
		sumline+=probabilities[j]


	verifRanks()
	
	sumMBG=0
	for i in T2check:
		sumMBG+=probabilities[i]

	maxRank=0
	maxRankID=0
	for i in range(len(graph)):#search a node with a maximum rank value
		if maxRank<graph[i].rank:
			maxRank=graph[i].rank
			maxRankID=graph[i].ID

	initData()
	setCoordinate()
	
	res=findDuplicate()
	while not (res[0]==0 and res[1]==0):
		graph[res[0]].coordinate+=1
		res=findDuplicate()
	
	buff="A[] ("
	for i in range(len(graph)-1):
		if i<len(graph)-2:
			buff+="msgR["+str(i+1)+"]==1 and "
		else:
			buff+="msgR["+str(i+1)+"]==1)"

	ofiQueries.write(buff+" imply a<=(3*(backoffD+rcvSlot+backoffD)+buzz)*DC\n")#write the query in a .q file
	
	n=0	
	os.popen("rm cases")
	for i in T2check:
		n+=1
		boolGraph=T[i]
		boolGraph2=neighbool2(T[i])
		setAutomata(args[0],"testout"+str(n)+".xml")
		os.popen("echo "+str(n)+" "+str(i)+" "+str(probabilities[i])+" >> cases")

if __name__ == '__main__':
    main()