Excerpted from the second link given below:
A magic square is a square matrix where all row and column sums and also the
diagonal sums all have the same value.
This value or the characteristic sum for a magic square of order n is sum(1:n^2)/n.
as per the wikipedia entry: this equation should form the "magic constant"
and yet using this in a simple test that works on any magic sq could be problematic:
In [36]:
my_mat = matrix(c(8, 3, 4, 1, 5, 9, 6, 7, 2), ncol = 3) # by column by default
my_mat # this is a magic square
Out[36]:
In [38]:
# altering a copy of my_mat so it is slightly off and will fail
my_mat2 <- my_mat
my_mat2[1,1] <- 50
my_mat2
Out[38]:
In [39]:
# another test ... more similar to original but should still fail
my_mat3 <- my_mat
my_mat3[2,3] <-10
my_mat3
Out[39]:
In [41]:
# test with a bigger matrix
# source of this matrix:
# http://www.sfcc.edu.hk/academic_subjects/mathematics/web/1999_2000_projects/Tracy%20Chan/magic%20square%20S.htm
my_larger_matrix <- matrix(c(1, 98,97,96,85,84,15,14,13, 2,
95,19,20,80,79,78,77,25,26, 6,
94,27,28,72,71,70,69,33,34, 7,
92,66,65,37,38,39,40,60,59, 9,
83,58,57,45,46,47,48,52,51,18,
12,50,49,53,54,55,56,44,43,89,
11,42,41,61,62,63,64,36,35,90,
10,67,68,32,31,30,29,73,74,91,
8 ,75,76,24,23,22,21,81,82,93,
99, 3, 4, 5,16,17,86,87,88,100),
ncol=10, byrow = TRUE)
my_larger_matrix
Out[41]:
In [42]:
# modified the larger matrix just a little so it should fail ...
my_larger_matrix2 <- matrix(c(1, 95,94,92,83,12,11,10, 8,99,
98,19,27,66,58,50,42,67,75,3,
97,20,28,65,57,49,41,68,76,4,
96,80,72,37,45,53,61,32,24,5,
85,79,71,38,46,54,62,31,23,16,
84,78,70,39,47,55,63,30,22,17,
15,77,69,40,48,56,64,29,21,86,
14,25,33,60,52,55,36,73,81,87,
13,26,34,59,51,43,35,74,82,88,
2, 6, 7, 9,18,89,90,91,93,100), ncol=10)
my_larger_matrix2
Out[42]:
This version of the code exits the magic square test as early as possible as soon as a total proves it is not a magic square.
In [12]:
# This approach:
# * Attempts to walk the matrix calculating as it goes .. it exits if any sums saved is not equal to other sums ...
# * this approach terminates on "first discovery" of error instead of having to wait
# until later in the loops to find it ... (at least in theory)
#
# This approach was written before learning of built in approaches that might work faster
is_magic_matrix <- function(a_matrix, suppressMsgs=FALSE) {
rws = nrow(a_matrix)
cls = ncol(a_matrix)
if (rws != cls) {
if (suppressMsgs == FALSE) {
print("Not a Magic Square:")
print(paste0(c("row/col dimensions not equal:", rws, cls), collapse=" "))
}
return(-1)
}
# 3x3 is the smallest magic square according to the web
# so this version of the code does not accept 2x2 as valid
if (rws < 3) {
if (suppressMsgs == FALSE) {
print("Not a Magic Square:")
print(paste0("row not > 3, rows = ", rws))
}
return(-1)
}
# future enhancement note: collapse to single gridline total?
rowTot <- 0
colTot <- 0
diagAccum <- 0
diagAccum2 <- 0
d2_offset <- 1
msg1 <- "Not a Magic Square - Failed this test:"
for(i in 1:rws) {
for(j in 1:cls) {
# this approach should only perform calculations when needed:
# Note: if(i==j==1) threw an error (unexpected "==")
if(i == 1 & j == 1) {
###############################################################
# we are in the first cell: a_matrix[1, 1]
# initialize these and do a quick test:
rowTot <- sum(a_matrix[i,1:cls])
colTot <- sum(a_matrix[1:rws, j])
if (rowTot != colTot) {
# failed row/col equiv test, get out w/ -1
if (suppressMsgs == FALSE) {
print(msg1)
print(paste0("a_matrix[1,1] - row and col sums:", rowTot, colTot,
collapse=" "))
}
return(-1)
}
# intialize the left to right diagonal w/ first value
diagAccum <- a_matrix[i, j]
# intialize right to left diagonal w/ first value
diagAccum2 <- a_matrix[i, cls]
###############################################################
} else {
# i, j, and 1 are not all equivelent, therefore this condition
# starts with any cell that is not a_matrix[1, 1]
if (i == 1) {
#############################################################
# we are on first row but not in a_matrix[1][1]
# here we can calculate and check rest of the columns
#
if (colTot == sum(a_matrix[1:rws, j])) {
# passed test for this column
} else {
# not all columns equal - get out:
if (suppressMsgs == FALSE) {
print(msg1)
print(paste0("a_matrix[1, ", j, "] - col sums not equal: ",
colTot, " ", sum(a_matrix[1:rws, j])))
}
return(-1)
}
#############################################################
}
if (j == 1) {
#############################################################
# we on first column but not in a_matrix[1][1]
# here we cna calculate and check the rest of the rows
#
if (rowTot == sum(a_matrix[i,1:cls])) {
# passed test for this column
} else {
# not all columns equal - get out:
if (suppressMsgs == FALSE) {
print(msg1)
print(paste0("a_matrix[1, ", j, "] - row sums not equal: ",
rowTot, " ", sum(a_matrix[i,1:cls])))
}
return(-1)
}
#############################################################
}
if (i == j) {
# starting with a_matrix[2, 2], then [3, 3] etc.
# this will build up the left to right diagonal
diagAccum <- diagAccum + a_matrix[i, j]
diagAccum2 <- diagAccum2 + a_matrix[i, (cls - d2_offset)]
# increment the offset after each use
d2_offset <- d2_offset + 1
}
} ###### (complex conditional ends here) ########################
} ##### (inner loop ends here) ####################################
} ##### (outer loop ends here) ######################################
if (diagAccum == diagAccum2 &
rowTot == colTot &
diagAccum == rowTot) {
# passed all tests give us back the magic number as return value:
return(rowTot)
} else {
# tell us what went wrong and return -1
if (suppressMsgs == FALSE) {
print("Almost made it but not a magic Square:")
print(paste0("row total : ", rowTot ))
print(paste0("col total : ", colTot ))
print(paste0("l-2-r diag: ", diagAccum))
print(paste0("r-2-l diag: ", diagAccum2))
}
return(-1)
}
# code ends here ...
} # final closing bracket ...
In [25]:
test_magic_squareFunction <- function(mtrx, suppressMsgs=FALSE) {
rtnVal = is_magic_matrix(mtrx, suppressMsgs)
if (suppressMsgs == FALSE) {
if (rtnVal == -1) {
print(paste0("Previous messages show how this matrix failed the magic ",
"square test."))
} else {
print("my_mat is a magic square w/ a magic number of:")
print(rtnVal)
}
}
}
print(my_mat)
test_magic_squareFunction(my_mat)
In [26]:
# altering a copy of my_mat so it is slightly off and will fail
my_mat2
test_magic_squareFunction(my_mat2)
Out[26]:
In [27]:
# another test ... more similar to original but should still fail
my_mat3
test_magic_squareFunction(my_mat3)
Out[27]:
In [28]:
# test with a larger magic square (this one should pass)
test_magic_squareFunction(my_larger_matrix)
In [29]:
# modified the larger matrix just a little so it should fail ...
test_magic_squareFunction(my_larger_matrix2)
In [33]:
# not a magic square = -1
# magic square returns the magic number
# this is how to turn off all output messages and just pass the return value
is_magic_matrix(my_mat, TRUE)
is_magic_matrix(my_mat2, TRUE)
Out[33]:
Out[33]:
In [34]:
# time performance metrics ... runs two magic square tests 1000 times
system.time({
for (i in 1:1000) {
is_magic_matrix(my_larger_matrix, TRUE)
is_magic_matrix(my_larger_matrix2, TRUE)
}
})
Out[34]:
In [58]:
is_magic <- function(mat, suppressMsgs=FALSE){
if (dim(mat)[1] != dim(mat)[2]) {
if (suppressMsgs==FALSE) {
print("Dimensions uneven - this is not even a square!")
}
return(-1)
}
diag1 = sum(diag(mat))
diag2 = sum(diag(mat[nrow(mat):1,]))
if (diag1 != diag2) {
if (suppressMsgs==FALSE) {
print("Diagonals don't match!")
}
return(-1)
} else {
rwsCalc = rowSums(mat)
for(i in 1:length(rwsCalc)) {
if (rwsCalc[i] != diag2) {
if (suppressMsgs==FALSE) {
print("Sum of a row did not match diagonals!")
}
return(-1)
}
}
colsCalc = colSums(mat)
for(i in 1:length(colsCalc)) {
if (colsCalc[i] != diag2) {
if (suppressMsgs==FALSE) {
print("Sum of a column did not match a diagonal!")
}
return(-1)
}
}
if (suppressMsgs==FALSE) {
print("This is a Magic square!")
}
return(diag2)
}
}
In [59]:
test_magic_squareFunction2 <- function(mtrx, suppressMsgs=FALSE) {
rtnVal = is_magic(mtrx, suppressMsgs)
if (suppressMsgs == FALSE) {
if (rtnVal == -1) {
print(paste0("Previous messages show how this matrix failed the magic ",
"square test."))
} else {
print(rtnVal)
}
}
}
print(my_mat)
test_magic_squareFunction2(my_mat)
In [60]:
# altering a copy of my_mat so it is slightly off and will fail
my_mat2
test_magic_squareFunction2(my_mat2)
Out[60]:
In [61]:
my_mat3
test_magic_squareFunction2(my_mat3)
Out[61]:
In [62]:
# test with a larger magic square (this one should pass)
test_magic_squareFunction2(my_larger_matrix)
In [64]:
# test with a larger magic square (this one should pass)
test_magic_squareFunction2(my_larger_matrix2)
In [65]:
# not a magic square = -1
# magic square returns the magic number
# this is how to turn off all output messages and just pass the return value
is_magic(my_mat, TRUE)
is_magic(my_mat2, TRUE)
Out[65]:
Out[65]:
In [66]:
# time performance metrics ... runs two magic square tests 1000 times
system.time({
for (i in 1:1000) {
is_magic(my_larger_matrix, TRUE)
is_magic(my_larger_matrix2, TRUE)
}
})
Out[66]:
In [ ]: