#Author: John Lawrence F. Quiñones
#Subject: CMSC 150 B-4L

#Author's note: Please set the working directory of R Studio to the directory of the extracted zip 
#               file containing "QuiñonesEx04.R" and "QuiñonesEx03.R" in order for this program to work

source('QuiñonesEx03.R')

VA= function(x1, x2, x3, x4, x5, x6, x7, x8, x9) 4 * x1 + -1 * x2 + 0 * x3 + -1 * x4 + 0 * x5 + 0 * x6 + 0 * x7 + 0 * x8 + 0 * x9 + -80;
VB= function(x1, x2, x3, x4, x5, x6, x7, x8, x9) -1 * x1 + 4 * x2 + -1 * x3 + 0 * x4 + -1 * x5 + 0 * x6 + 0 * x7 + 0 * x8 + 0 * x9 + -30;
VC= function(x1, x2, x3, x4, x5, x6, x7, x8, x9) 0 * x1 + -1 * x2 + 4 * x3 + 0 * x4 + 0 * x5 + -1 * x6 + 0 * x7 + 0 * x8 + 0 * x9 + -80;
VD= function(x1, x2, x3, x4, x5, x6, x7, x8, x9) -1 * x1 + 0 * x2 + 0 * x3 + 4 * x4 + -1 * x5 + 0 * x6 + -1 * x7 + 0 * x8 + 0 * x9 + -50;
VE= function(x1, x2, x3, x4, x5, x6, x7, x8, x9) 0 * x1 + -1 * x2 + 0 * x3 + -1 * x4 + 4 * x5 + -1 * x6 + 0 * x7 + -1 * x8 + 0 * x9 + 0;
VF= function(x1, x2, x3, x4, x5, x6, x7, x8, x9) 0 * x1 + 0 * x2 + -1 * x3 + 0 * x4 + -1 * x5 + 4 * x6 + 0 * x7 + 0 * x8 + -1 * x9 + -50;
VG= function(x1, x2, x3, x4, x5, x6, x7, x8, x9) 0 * x1 + 0 * x2 + 0 * x3 + -1 * x4 + 0 * x5 + 0 * x6 + 4 * x7 + 0 * x8 + 0 * x9 + -70;
VH= function(x1, x2, x3, x4, x5, x6, x7, x8, x9) 0 * x1 + 0 * x2 + 0 * x3 + 0 * x4 + 0 * x5 + 0 * x6 + -1 * x7 + 4 * x8 + -1 * x9 + -70;
VI= function(x1, x2, x3, x4, x5, x6, x7, x8, x9) 0 * x1 + 0 * x2 + 0 * x3 + 0 * x4 + 0 * x5 + -1 * x6 + 0 * x7 + -1 * x8 + 4 * x9 + -120;
ls = list(eq1=VA, eq2=VB, eq3=VC, eq4=VD, eq5=VE, eq6=VF, eq7=VG, eq8=VH, eq9=VI)


lbldList = AugCoeffMatrix(ls)
print(lbldList)

findPivotRow = function(n, mat){
  for (i in 1:n){
    #finding the pivot row
    pivotRow = NA
    for (j in i:n){
      if(is.na(pivotRow)){
        pivotRow = j
      }else{
        if(mat[pivotRow, 1] < mat[j, 1]){
          print("Pivot:")
          print(pivotRow)
          print("Eval:")
          print(j)
          print(mat[pivotRow, 1])
          print("is less than")
          print(mat[j, 1])
          pivotRow = j
        }
      }
    }
    
    #print("Pivot row:")
    #print(pivotRow)
  }
  return(pivotRow)
}
  
getRowArrangement = function(i, n, pivotRow, mat){
  rowArrangement = c()
  #dynamically getting the row names arrangement but with i and pivotRow swapped positions
  #print("Row Arrangement: ")
  for (k in 1:n){
    if(k == i){
      rowArrangement = c(rowArrangement, rownames(mat)[pivotRow])
    }else if(k == pivotRow){
      rowArrangement = c(rowArrangement, rownames(mat)[i])
    }else{
      rowArrangement = c(rowArrangement, rownames(mat)[k])
    }
    #print(rowArrangement)
  }
  return(rowArrangement)
}

GaussianMethod= function(ls){
  variables = ls[[2]]
  mat = ls[[1]]
  
  n = nrow(mat)
  
  for (i in 1:(n-1)){
    pivotRow = findPivotRow(n, mat)
      
    #if the first element of the pivot row is zero, return NA
    if((mat[pivotRow, i] == 0) & i != (n-1)){
      print("No solution for: ")
      print(mat)
      print("Pivot row: ")
      print(mat[pivotRow, ])
      return(NA)
    }
    
    #print("Initial Matrix:")
    #print(mat)
    
    rowArrangement = getRowArrangement(i, n, pivotRow, mat)
  
    #pivoted matrix
    mat = mat[match(rowArrangement, rownames(mat)), ]
    #print("pivoted matrix:")
    #print(mat)
    
    for (l in (i+1):n){
      pivotElement = mat[i, i]
      #print("Pivot element: ")
      #print(pivotElement)
      multiplier = (mat[l, i])/pivotElement
      
      #print("multiplier: ")
      #print(multiplier)
      
      normalizedRow = multiplier * mat[i, ]
      #print("Normalized Row:")
      #print(normalizedRow)
      
      differenceRow = mat[l, ] - normalizedRow
      #fix this part Mr. AI
      mat[l, ] <- differenceRow
      
      #print("difference:")
      #print(differenceRow)
      #print("Resulting Matrix:")
      #print(mat)
    }
  }
  
  #backwards substitution
  varValues = c()
  
  #getting the value of x[n]
  xSubN = mat[n, ncol(mat)] / mat[n,n]
  
  #initialize x vector with 0 as initial values
  x = c()
  for (i in 1:n){
    x <- c(x, 0)
  }
  
  #putting xSubN in x[n]
  x[n] <- xSubN
  
  for (i in (n-1): 1){
    sum = 0
    for (j in (i+1):n){
      prod = mat[i,j] * x[j]
      sum <- prod + sum
    }
    diff = mat[i, ncol(mat)] - sum
    quo = diff/mat[i,i]
    x[i] = quo
  }
  
  retVal = list(variables=variables, augcoeffmatrix = round(mat, digits=4), solution = x)
  
  return(retVal)
}

GaussJordanMethod = function(ls){
  variables = ls[[2]]
  mat = ls[[1]]
  
  n = nrow(mat)
  
  for (i in 1:n){
    if (i != n){
      pivotRow = findPivotRow(n, mat)
      
      #if the first element of the pivot row is zero, return NA
      if((mat[pivotRow, i] == 0) & i != (n-1)){
        print("No solution!")
        return(NA)
      }
      
      rowArrangement = c()
      #dynamically getting the row names arrangement but with i and pivotRow swapped positions
      #print("Row Arrangement: ")
      for (k in 1:n){
        if(k == i){
          rowArrangement = c(rowArrangement, rownames(mat)[pivotRow])
        }else if(k == pivotRow){
          rowArrangement = c(rowArrangement, rownames(mat)[i])
        }else{
          rowArrangement = c(rowArrangement, rownames(mat)[k])
        }
        #print(rowArrangement)
      }
      
      #pivoted matrix
      mat = mat[match(rowArrangement, rownames(mat)), ]
      #print("pivoted matrix:")
      #print(mat)
    }
    #print("Row to be normalized: ")
    #print(mat[i, ])
    #print("Divisor: ")
    #print(mat[i, i])
    
    #normalizing the pivot row
    mat[i, ] <- mat[i, ] / mat[i, i]
    #print("NORMALIZED: ")
    #print(mat)
    
    #subtracting the normalized row to mat[j, ] to create a diagonal matrix and extract the unknown variables
    for (j in 1:n){
      if (i == j){
        next
      }
      normalizedRow = mat[j, i] * mat[i, ]
      mat[j, ] <- mat[j, ] - normalizedRow
      
      #print("Resulting matrix with normalized row: ")
      #print(mat)
    }
  }
  
  #print("Final matrix: ")
  #print(mat)
  
  #saving the rhs values (solutions) to a vector
  x = c()
  for (i in 1:n){
    x[i] <- mat[i, ncol(mat)]
  }
  
  #print("Solutions: ")
  #print(x)
  
  retVal = list(variables=variables, augcoeffmatrix = round(mat), solution = x)
  
  return(retVal)
}

GaussianMethod(lbldList)
GaussJordanMethod(lbldList)