library(MASS)

# vertices are extended by constant 1 
dataVertices=read.table("stand-imset4-vertices.txt",header=FALSE,sep=",")
# 185 x 12 matrix
V=data.matrix(dataVertices)

# vertices are extended by constant 1 
dataB=read.table("stand-imset4-neigbourhood-matrix.txt",header=FALSE,sep=",")
# 185 x 185 matrix - B[i1,i2]=1 if i1 and i2 are neighbours
B=data.matrix(dataB)
# a copy of matrix B
C=B

source("imset.R")

# This cache will be used for storing values of current structural imset
cacheImset1=new.env(hash=TRUE)
sets=list(c(1,2),c(1,3),c(2,3),c(1,4),c(2,4),c(3,4),c(1,2,3),c(1,2,4),c(1,3,4),c(2,3,4),c(1,2,3,4))
# for all imsets corresponding to vertices of polytop do
for (v in 1:nrow(V)){
  # for all subsets of N of imset v
  for (i in 1:(ncol(V)-1)) setIms(sets[[11]],sets[[i]],V[v,i],cacheImset1);
  for (w in 1:ncol(B))
    # if the second imset is a geometric neighbour compute differential imset
    if (B[v,w]==1){
      for (j in 1:(ncol(V)-1)) setIms(sets[[11]],sets[[j]],-V[w,j],cacheImset1);
      # if the differential imset is elementary replace "1" by "2" in matrix C
      if (isImsElementary(cacheImset1,sets[[11]])) C[v,w]=2;
      for (j in 1:(ncol(V)-1)) setIms(sets[[11]],sets[[j]],V[w,j],cacheImset1);
    }      
  for (i in 1:(ncol(V)-1)) setIms(sets[[11]],sets[[i]],-V[v,i],cacheImset1);
}

write.matrix(C, file = "stand-imset4-neigbourhood-matrix-marked-incl-nb.txt", sep = ",")