Key reference: McKay (2003). Algorithm 42.9, Replicate Figure 42.3

# Binary Hopfield network (pp506-507) w/o gradient descent
I <- 25; # 25 neurons
N <- 4; # 4 desired letters/patterns/memories; see below
letter.d <- c( 1, 1, 1, 1,-1,
              -1, 1,-1,-1, 1,
              -1, 1,-1,-1, 1,
              -1, 1,-1,-1, 1,
              -1, 1, 1, 1,-1);
letter.j <- c(  1, 1, 1, 1, 1,
               -1,-1,-1, 1,-1,
               -1,-1,-1, 1,-1,
                1,-1,-1, 1,-1,
                1, 1, 1,-1,-1);
letter.c <- c( -1, 1, 1, 1, 1,
                1,-1,-1,-1,-1,
                1,-1,-1,-1,-1,
                1,-1,-1,-1,-1,
               -1, 1, 1, 1, 1);
letter.m <- c(  1,-1,-1,-1, 1,
                1, 1,-1, 1, 1,
                1,-1, 1,-1, 1,
                1,-1,-1,-1, 1,
                1,-1,-1,-1, 1);

# show letters:
draw.letter <- function (x) {
  for (i in 1:length(x)) {
    cat(ifelse(x[i]==1, "*", " "));    
    if(i%%5 == 0) {cat("\n")};
  }
};
draw.letter(letter.d)
## **** 
##  *  *
##  *  *
##  *  *
##  ***
draw.letter(letter.c)
##  ****
## *    
## *    
## *    
##  ****
draw.letter(letter.j)
## *****
##    * 
##    * 
## *  * 
## ***
draw.letter(letter.m)
## *   *
## ** **
## * * *
## *   *
## *   *
# start network
x <- matrix(c(letter.d, letter.j, letter.c, letter.m), nrow = N, byrow = T)
w <- t(x) %*% x; # initialize wts with Hebb rule: outer products of matrix multiplication
for (i in 1:I) { w[i,i] = 0 } #  ensure zero self-wts (no self-feed)

# error correction with sequential updates
mutate.letter <- function (letter, num.error) {
  draw.letter(letter);
  pos <- sample(1:length(letter), num.error);
  letter[pos] <- ifelse(letter[pos]==1, -1, 1);
  draw.letter(letter);
  return(letter);
}

hopfield.binary <- function (error.letter, iteration, memory = w) {
  draw.letter(error.letter);
  a <- integer(I); # activities for each neuron
  current.letter <- error.letter;
  for (n in 1:iteration) { # run  5 iterations, enough to restore 1-5 mutations
#    for (i in 1:I) { # for each neuron
#      a[i] <- 0;
#      for (j in 1:length(current.letter)) {
#        a[i] <- a[i] + w[i,j] * current.letter[j]
#      }
#    }
    a <- w %*% current.letter; 
    current.letter <- ifelse(a>0, 1, -1)
    draw.letter(current.letter)
  }
}

x <- mutate.letter(letter.d, 5)
## **** 
##  *  *
##  *  *
##  *  *
##  *** 
## **** 
##  ** *
##  *  *
## ** **
##   *
hopfield.binary(x, iteration = 10)
## **** 
##  ** *
##  *  *
## ** **
##   *  
## **** 
##  *  *
##  *  *
##  *  *
##  *** 
## **** 
##  *  *
##  *  *
##  *  *
##  *** 
## **** 
##  *  *
##  *  *
##  *  *
##  *** 
## **** 
##  *  *
##  *  *
##  *  *
##  *** 
## **** 
##  *  *
##  *  *
##  *  *
##  *** 
## **** 
##  *  *
##  *  *
##  *  *
##  *** 
## **** 
##  *  *
##  *  *
##  *  *
##  *** 
## **** 
##  *  *
##  *  *
##  *  *
##  *** 
## **** 
##  *  *
##  *  *
##  *  *
##  *** 
## **** 
##  *  *
##  *  *
##  *  *
##  ***
# Applications: Pattern recognition; Error correction; Optimization for NP-complete constrainted optimization (e.g., traveling saleman) (need gradient descent)