# Stochasticke simulacne metody
# Radoslav Harman, KAMS FMFI UK
# Cvicenie 8 (Isingov model) - strucny zaznam prikazov pre R

ising <- function(N, tau, Bx1, n, Graph=0) {
  # Isingov model simulovany diskretnym Metropolisovym algoritmom
  # N ... rozmer spinovej mriezky
  # tau ... konstanta umerna teplote systemu
  # Bx1 ... startovaci stav (matica NxN hodnot -1 a 1)
  # n ... pozadovany pocet realizacii
  # Graph ... ci sa ma zobrazovat ilustracny graf

  Bx <- Bx1              # momentalny stav

  for(i in 1:n){

    # Zobrazenie momentalneho stavu
    if((Graph > 0) & (i %% Graph == 0)) image(x=1:N, y=1:N, z=Bx)

    # Novy kandidat na stav By, ktory sa lisi od Bx v jednom bode preklopenia iu,ju
    iu <- ceiling(N * runif(1))
    ju <- ceiling(N * runif(1))
 
    # Vypocet pomeru pravd. kandidatskeho a momentalneho stavu ratio=p[By]/p[Bx]
    s <- 0
    if(iu > 1) s <- s + Bx[iu - 1, ju]
    if(ju > 1) s <- s + Bx[iu, ju - 1]
    if(iu < N) s <- s + Bx[iu + 1, ju]
    if(ju < N) s <- s + Bx[iu, ju + 1]
    ratio <- exp(-(2 / tau) * Bx[iu, ju] * s)

    # Potencialna zmena stavu za kandidata
    alpha <- min(c(1, ratio))
    if(runif(1) < alpha) Bx[iu, ju] <- (-Bx[iu, ju])

  }

  return(Bx)
}

ising(200, 1, matrix(sample(c(-1, 1), 40000, replace=TRUE), ncol=200), 5000000, 100000)
ising(200, 2, matrix(sample(c(-1, 1), 40000, replace=TRUE), ncol=200), 5000000, 100000)
ising(200, 3, matrix(sample(c(-1, 1), 40000, replace=TRUE), ncol=200), 5000000, 100000)
ising(200, 4, matrix(sample(c(-1, 1), 40000, replace=TRUE), ncol=200), 5000000, 100000)