#Welcome to the code to calculate Linking Coherency Index Scores
# All you need to to is 1) define the number of alternatives (nalts) and 2) total number of criteria (totalcrit) 
# in the supermatrix and 3) import your weighted Supermatrix and run the code

#Here is where you tell the program how many alternatives and criteria are in your model
nalts=3
totalcrit=6

#Simply provide the Supermatrix as a textfile. It should be in the form: (numbers only no labels)
#                 Criteria  Alternatives
#       Criteria   
#   Alternatives
#
#change to your file path where you pasted the supermatrix in text document
weightedsuper=read.table("C:/Users/olcooper/Desktop/weightedsuper.txt")
weightedsuper=as.matrix(weightedsuper,nalts+totalcrit,nalts+totalcrit)




#choose the alternative and criterion that you will use to anchor (really doesn't matter but this way you can choose)

altanchor=2
critanchor=2

#what power to raise Supermatrix to - must be even number

powersuper=100


##################################################################################################
#now for the functions that will be called in the main method

#define a cutoff value to deal with floats essentially if you are smaller than ... you are = to zero
cu=.0000001

###This will raise the Supermatrix to powers to provide the limit matrix
raisepowers<-function(matrix1,powers){
  limit=matrix1
  for (i in 1:powers){
    limit=limit%*%matrix1
  }
  return(limit)
}

#This function is called by the function "estimations" that will be used to normalizes the column into ratios 
# of specific entries in each column. This is part of step 1b

normalizecolumn<-function(weightedsuper,colnum,anchr,totalcrit,nalts,top){
#note the top 0,1 may seem counter intuitive here because you normalize the other half to use in the linking
#process
  
  if (top==1){
    size=nalts
    anchor=weightedsuper[totalcrit+anchr,colnum]
    place=totalcrit
  }
  else{
    colnum=colnum+totalcrit
    size=totalcrit
    anchor=weightedsuper[anchr,colnum]
    place=0
  }
  
  normalizedcol=matrix(0,size,1)

  for (i in 1:size){
    normalizedcol[i,1]=weightedsuper[place+i,colnum]/anchor
  }

  return(normalizedcol)
}

#This function is called by "estimations" to convert the normalized entries into data that is in the same units as part of step 1b which 
#will then be used to get the linking estimations in step 1c

altlinkingpin<-function(weighted,altanc,critanchor,step1,nalts,totalcrit,top){
  if(top==1){
    size=totalcrit
    linkinganchor=weighted[critanchor,totalcrit+altanc]
  }
  else{
    size=nalts
    linkinganchor=weighted[totalcrit+altanc,critanchor]
  }
  
  constants=matrix(0,size,1)
  linkedstep2=step1
  
  if(top==1){
    for (i in 1:totalcrit){
      constants[i,1]=weighted[i,altanc+totalcrit]/linkinganchor
    }
    for (i in 1:totalcrit){
      for (j in 1:nalts){
        linkedstep2[j,i]=linkedstep2[j,i]*constants[i,1]
      }
    } 
    
  }
  else{
    for (i in 1:nalts){
      constants[i,1]=weighted[totalcrit+i,critanchor]/linkinganchor
    }
    for (i in 1:nalts){
      for (j in 1:totalcrit){
        linkedstep2[j,i]=linkedstep2[j,i]*constants[i,1]
      }
    }  
  }
  return(linkedstep2)   
}


#this function is called by estimations and is used to get the data for the bottom half of the entries in the linking estimations

bottomhalffromlinked<-function(lvs,totalcrit,nalts){
  matrixtopass=matrix(0,nalts+totalcrit,nalts+totalcrit)
  for (i in 1:totalcrit){
    coltotal=0
    for (k in 1:nalts){
      coltotal=coltotal+lvs[k,i]
    }
    for (j in 1:nalts){
      matrixtopass[totalcrit+j,i]=lvs[j,i]/coltotal
    }
  }
  return(matrixtopass)
}

#This is step 1 that also calls the "normalize column, altlinkingpin, and bottom half from linked to calculate
#the linking estimations and then the other half of the linking estimations. This is step 1c

estimations<-function(c,weightedsuper,totalcrit,nalts,anchor,top){
   #step1 inside if loops because of dimensions
  linker=c
  if(top==0){
    step1=matrix(0,totalcrit,nalts)
    for (i in 1:nalts){
      step1[,i]=normalizecolumn(weightedsuper,i,linker,totalcrit,nalts,top)
    }
    
    step2=altlinkingpin(weightedsuper,anchor,linker,step1,nalts,totalcrit,top)
    #step2=altlinkingpin(weightedsuper,anchor,linker,step1,nalts,totalcrit,top)
    allsameunit=t(step2)

  }
  
  else
  {
    step1=matrix(0,nalts,totalcrit)
    for (i in 1:totalcrit){
      step1[,i]=normalizecolumn(weightedsuper,i,linker,totalcrit,nalts,top)
    }
  
    step2=altlinkingpin(weightedsuper,linker,altanchor,step1,nalts,totalcrit,top)
    allsameunit=step2
    
  }

  recoveredweightedsuper=bottomhalffromlinked(allsameunit,totalcrit,nalts)
  
  
  for (i in 1:nalts){
    rowtotal=0
    for (k in 1:totalcrit){
      rowtotal=rowtotal+allsameunit[i,k]
    }
    
    for (j in 1:totalcrit){
      recoveredweightedsuper[j,totalcrit+i]=allsameunit[i,j]/rowtotal
    }
  }
  return(recoveredweightedsuper)
}

#Here we calculate the LCIs for each comparison of the linking estimations in Step 2a.

lcimatrixpieces<-function(top,linkingsample,cu,nalts,totalcrit){
  
  if (top==0){
    startrow=1+totalcrit
    startcolumn=1
    size=totalcrit
    oppsize=nalts
  } else {
    startrow=1
    startcolumn=1+totalcrit
    size=nalts
    oppsize=totalcrit
  }
  
  
  lcimatrix=diag(1,size,size)
  
  
  for (i in 1:(size-1)){
    a=linkingsample[,,i]
    for (j in (i+1):size){
      sumd=0
      b=linkingsample[,,j]
      for (m in startrow:(startrow+oppsize-1)){
        for (n in startcolumn:(startcolumn+size-1)){
          if (b[m,n]>cu){
            b[m,n]=1/b[m,n]
          } else {
            b[m,n]=0
          }
        }
      }
      for (k in startrow:(startrow+oppsize-1)){
        for (l in startcolumn:(startcolumn+size-1)){
          
          d=a[k,l]*b[k,l]
          
          if (d<1 && d>cu){
            d=1/d
          }
          sumd=sumd+d
        }                  
      }
      tentativelci=sumd/(nalts*totalcrit)
      lcimatrix[i,j]=tentativelci
      
    }
  }
  
  for (l in 1:size){
    for (m in 2:size){
      lcimatrix[m,l]=lcimatrix[l,m]
    }
  }
  return(lcimatrix)
}


#Now we calculate LCI Scores from LCIs in Step 2b.

lciscores<-function(lcimatrix,squaresize){
  scores=c(squaresize)
  for (i in 1:squaresize){
    sumi=0
    for (j in 1:squaresize){
      sumi=lcimatrix[i,j]+sumi
    }

    scores[i]=((sumi-1)/(squaresize-1))
  }
  return(scores)
}  


#Here is the main code

#Creating some new variables

  pagesforcrit=totalcrit
  pagesforalts=nalts
  
  linkingestimationscrit=array(0,c(nalts+totalcrit,nalts+totalcrit,pagesforcrit))
  linkingestimationsalts=array(0,c(nalts+totalcrit,nalts+totalcrit,pagesforalts))
  

#Step 1 create the linking estimations

    for (s in 1:totalcrit){
      linkingestimationscrit[,,s]=estimations(s,weightedsuper,totalcrit,nalts,altanchor,0)
    }
    for (s in 1:nalts){
      linkingestimationsalts[,,s]=estimations(s,weightedsuper,totalcrit,nalts,critanchor,1)
    }
#Step 2 Calculate the LCI  
  #2a
    lciscrit=lcimatrixpieces(0,linkingestimationscrit,cu,nalts,totalcrit)
    lcisalts=lcimatrixpieces(1,linkingestimationsalts,cu,nalts,totalcrit)
  #2b
    lciscorescrit=lciscores(lciscrit,totalcrit)
    lciscoresalts=lciscores(lcisalts,nalts)

#Some formatting to make the results more readable

lciscorescrit=matrix(lciscorescrit,1,totalcrit)
colnames(lciscorescrit)=sprintf("Crit%d", 1:totalcrit)
lciscoresalts=matrix(lciscoresalts,1,nalts)
colnames(lciscoresalts)=sprintf("Alt%d", 1:nalts)

#You can call the limit matrix from the command prompt to see what it looks like 
#in each iteration or once your Supermatrix is coherent

limitmatrix=raisepowers(weightedsuper,powersuper)

#Printing results to the screen

print("lciscores alts")
print(lciscoresalts)
print("lciscores crit")
print(lciscorescrit)
print("Step 3 Choose the highest LCI Score")
print("Step 4 Revise the priority vector chosen in Step 3")
print("Step 5 Update 'weightedsuper' and repeat until the level of coherency in the Supermatrix is acceptable")
    
#Here you can write the output to file if you would like

#write.table(limit,'C:/(fill in where you want it)/limit.txt', sep='\t',row.names = FALSE,)
#write.table(weightedsuper,'C:/(fill in where you want it)/newweightedsuper.txt', sep='\t',row.names = FALSE,)