######### variables ###########

kantennummer_max=20
emax=30
fmax=30

###############################



###### create databases #######

database=matrix(list(),emax,fmax)
database_out=matrix(list(),emax,fmax)

###############################



## create starting point and ##
## put it into the database ###

tetrahedron <- matrix(c(
0,1,1,1,
1,0,1,1,
1,1,0,1,
1,1,1,0
),nrow=4, byrow=TRUE)

database[4,4][[1]]=list(tetrahedron)

###############################



## functions to test whether the poly- ##
## tope already exists in the database ##

sums=function(polytop)
{
c(sort(colSums(polytop)),sort(rowSums(polytop)))
}

test_sums=function(polytop)
{
liste=database[[nrow(polytop),ncol(polytop)]]
w=1
if(length(liste)>0){
	for(i in 1:length(liste)){
		if(all(sums(polytop)==sums(liste[[i]]))){
			w0=0
		}else{
			w0=1}
	w=w*w0}
}else{
	w=2}
w
}

########################################



######## function to write a ########
######## polytop into the db ########

write_to_db=function(polytop){
en=ncol(polytop)
fn=nrow(polytop)
if(test_sums(polytop)==1){
	#print("matrix isn't in the DB yet")
	database[fn,en][[1]]=list(database[fn,en][[1]][[1]], polytop)
	database
}else if(test_sums(polytop)==2){
	#print("matrix isn't in the DB entry yet and DB is empty")
	database[fn,en][[1]]=list(polytop)
	database
}else{
	#print("matrix already exists in the DB")
	database
}
}

#####################################



###### process functions to #########
###### transform polytopes ##########

prozess_p=function(e, f, n, ecke, db) # cut off a vertex
{
polytop=db[f,e][[1]][[n]]
eckenrang=sum(polytop[,ecke])
polytop_backup=polytop
polytop1=polytop[polytop[,ecke]==1,]
zusatz1=matrix(1,eckenrang,eckenrang)
zusatz1=zusatz1-diag(eckenrang)
polytop1=cbind(polytop1,zusatz1)
polytop0=polytop[polytop[,ecke]==0,]
polytop=matrix(0,(nrow(polytop_backup)+1),(ncol(polytop_backup)+eckenrang))
polytop[0:eckenrang,0:(ncol(polytop_backup)+eckenrang)]=polytop1
polytop[(eckenrang+1):nrow(polytop_backup),0:ncol(polytop_backup)]=polytop0
polytop[nrow(polytop_backup)+1,(ncol(polytop_backup)+1):(ncol(polytop_backup)+eckenrang)]=1
polytop=polytop[,-ecke]
en=nrow(polytop)
fn=ncol(polytop)
db=write_to_db(polytop)
}



prozess_q=function(e, f, n, flaeche, db) # put a pyramid onto a face
{
polytop=db[f,e][[1]][[n]]
flaechenrang=sum(polytop[flaeche,])
polytop_backup=polytop
polytop1=polytop[,polytop[flaeche,]==1]
zusatz1=matrix(1,flaechenrang,flaechenrang)
zusatz1=zusatz1-diag(flaechenrang)
polytop1=rbind(polytop1,zusatz1)
polytop0=polytop[,polytop[flaeche,]==0]
polytop=matrix(0,(nrow(polytop_backup)+flaechenrang),(ncol(polytop_backup)+1))
polytop[0:(nrow(polytop_backup)+flaechenrang),0:flaechenrang]=polytop1
polytop[0:nrow(polytop_backup),(flaechenrang+1):ncol(polytop_backup)]=polytop0
polytop[(nrow(polytop_backup)+1):(nrow(polytop_backup)+flaechenrang),ncol(polytop_backup)+1]=1
polytop=polytop[-flaeche,]
en=nrow(polytop)
fn=ncol(polytop)
db=write_to_db(polytop)
}



prozess_r=function(e, f, n, ecke1, ecke2, db) # cut off a vertex (cut goes through another vertex
{
polytop=db[f,e][[1]][[n]]
eckenrang=sum(polytop[,ecke1])
polytop_backup=polytop
polytop1=polytop[polytop[,ecke1]==1,]
zus1=matrix(0,eckenrang,(eckenrang-1))
zus2=matrix(0,eckenrang,(eckenrang-1))
zus1[0:(eckenrang-1),0:(eckenrang-1)]=diag(eckenrang-1)
zus2[2:(eckenrang),0:(eckenrang-1)]=diag(eckenrang-1)
zusatz1=zus1+zus2
polytop1=cbind(polytop1,zusatz1)
polytop0=polytop[polytop[,ecke1]==0,]
polytop=matrix(0,(nrow(polytop_backup)+1),(ncol(polytop_backup)+eckenrang-1))
polytop[0:eckenrang,0:(ncol(polytop_backup)+eckenrang-1)]=polytop1
polytop[(eckenrang+1):nrow(polytop_backup),0:ncol(polytop_backup)]=polytop0
polytop[(nrow(polytop_backup)+1),ecke2]=1
polytop[(nrow(polytop_backup)+1),(ncol(polytop_backup)+1):(ncol(polytop_backup)+eckenrang-1)]=1
polytop=polytop[,-ecke1]
en=nrow(polytop)
fn=ncol(polytop)
db=write_to_db(polytop)
}

####################################



########### main program ###########

for(k in 6:kantennummer_max){
	eckennummer=k-2
	for(eckenposition in 1:(eckennummer-3)){
		e=eckennummer-eckenposition+1
		f=k-e+2
		if(length(database[e,f][[1]])>0){
			for(n in length(database[e,f][[1]])){
				for(ecke in 1:e){
					#print(c("p",e,f,n,ecke))
					database=prozess_p(e,f,n,ecke,database)
					#print("done")
				}
				for(flaeche in 1:f){
					#print(c("q",e,f,n,flaeche))
					database=prozess_q(e,f,n,flaeche,database)
					#print("done")
				}
				#database=prozess_r(e,f,n,1,1,database) #TODO
			}	
		}
	}
}

###################################



#### get numbers out of the db ####

for(e in 1:ncol(database)){
	for(f in 1:nrow(database)){
		database_out[e,f]=length(database[e,f][[1]])
	}
}
database_out

###################################