Lab8_Ising-Model.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "plt.rcParams[\"figure.figsize\"] = (10, 6)  # increase default figure size\n",
    "\n",
    "np.random.seed(2016)  # for reproducibility"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Ising Model\n",
    "\n",
    "\n",
    "Στην άσκηση αυτή θα χρησιμοποιήσουμε τον αλγόριθμο **Metropolis-Hastings** για να\n",
    "κάνουμε δειγματοληψία από την αναλλοίωτη κατανομή $\\pi$ του μοντέλου [Ising](https://en.wikipedia.org/wiki/Ising_model). \n",
    "Στο επόμενο κελί κώδικά θεωρούμε πλέγμα μήκους `L = 32`, \n",
    "οπότε ο χώρος $\\mathbb{X}$ των δυνατών διαμορφώσεων έχει πληθάριθμο $2^{1024} > 10^{307}$. \n",
    "Καταλαβαίνει κανείς ότι η δειγματοληψία από έναν τόσο μεγάλο χώρο είναι\n",
    "πρακτικά αδύνατη με συμβατικές μεθόδους. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "L = 32  # length of lattice\n",
    "N = L * L  # number of particles\n",
    "nsteps = N * 800  # number of updates, we will update each spin 800 times on average \n",
    "Temperature = 30.  # in units where the Boltzmann factor is 1\n",
    "\n",
    "# Initialize Variables\n",
    "spins = 2 * (np.random.rand(L, L) > 0.5) - 1  # random initial state (table of {1, -1})\n",
    "kernel = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]])\n",
    "\n",
    "for n in range(nsteps):\n",
    "    r, c = np.random.randint(0, L, 2)  # select a vertex (row, column) at random\n",
    "    # indices of a box around r, c\n",
    "    C = [(c-1) % L, c % L, (c + 1) % L]  # periodic boundary conditions, i.e. (0,x) is a neighbour of (L-1,x) \n",
    "    R = [[(r-1) % L], [r % L], [(r+1) % L]]\n",
    "    DH = spins[r, c] * np.sum(kernel * spins[R, C]) # check the textbook to see why this expression is ΔΗ\n",
    "    if np.random.rand() < np.exp(- DH / Temperature):  # think why this realizes Metropolis-Hastings (discussed in class)\n",
    "        spins[r, c] *= -1  # flip spin \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ο κώδικας χρησιμοποιεί τον αλγόριθμο Metropolis-Hastings \n",
    "για να προσομοιώσει μια μαρκοβιανή αλυσίδα στον $\\mathbb{X}$ \n",
    "που ξεκινά από μια τυχαία αρχική διαμόρφωση (`spins`) και έχει ανναλοίωτη κατανομή το μέτρο Gibbs που αντιστοιχεί στην ενέργεια του μοντέλου Ising, δηλαδή\n",
    "\n",
    "$$H(\\sigma)=-\\sum_{x,y: x~y} \\sigma(x)\\sigma(y)$$\n",
    "και\n",
    "$$\\pi(\\sigma)=\\frac{e^{-\\beta H(\\sigma)}}{Z(\\beta)}$$\n",
    "όπου η συνάρτηση επιμερισμού $Z(\\beta)$ δίνεται από την\n",
    "$$Z(\\beta)=\\sum_{\\sigma\\in\\mathbb{X}} e^{-\\beta H(\\sigma)}.$$\n",
    "\n",
    "\n",
    "Σε κάθε επανάληψη (`nsteps`)\n",
    "επιλέγει τυχαία ένα σημείο του πλέγματος (`r, c`) και υπολογίζει τη διαφορά στην ενέργεια \n",
    "του συστήματος (`DH`) αν το spin σε εκείνο το σημείο αλλάξει πρόσημο. \n",
    "Στη συνέχεια δέχεται η απορρίπτει την αλλαγή spin με πιθανότητα $e^{-\\beta (DH)^+}$, όπως ορίζει ο αλγόριθμος Metropolis-Hastings.\n",
    "\n",
    "Στο επόμενο κελί κώδικα οπτικοποιούμε την τελική κατάσταση του συστήματος"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(spins, cmap='binary', vmin=-1, vmax=1, interpolation='nearest')\n",
    "plt.title(\"Temperature: %.1f, Grid Size: %d\" % (Temperature, L))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Δοκιμάσετε να ξανατρέξετε τα 2 προηγούμενα κελιά και δείτε πώς αλλάζει (αν αλλάζει) η τελική διάταξη του πλέγματος"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Ising Parameters\n",
    "\n",
    "Για να αποκτήσουμε μια καλύτερη διαίσθηση για το πως λειτουργεί το σύστημα,\n",
    "θα μελετήσουμε την συμπεριφορά του αν πειράξουμε τις διάφορες παραμέτρους.\n",
    "\n",
    "Πρώτα θα δούμε πώς η ενέργεια επηρεάζει την πιθανότητα αλλαγής σπιν.\n",
    "Στο παρακάτω κελί κώδικα πλοτάρουμε την πιθανότητα μετάβασης για όλες τις\n",
    "πιθανές τιμές τις ενέργειας."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.arange(-4, 5)  # all possible DH after a spin flip\n",
    "Temperature = 3.0\n",
    "y = np.exp(-np.maximum(x,0) / Temperature)  \n",
    "plt.plot(x, y)\n",
    "plt.xlabel(\"Energy Difference\")\n",
    "plt.ylabel(\"Transition Probability\")\n",
    "plt.grid(True)\n",
    "plt.axvline(0, color='red', linestyle='dashed')\n",
    "#plt.axhline(1, color='red', linestyle='dashed')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Παρατηρήστε ότι όταν η ενέργεια του συστήματος μειώνεται (`DH < 0`),\n",
    "ο αλγόριθμος δέχεται αυτομάτως την προτεινόμενη αλλαγή. Αντιθέτως όταν η ενέργεια του συστήματος αυξάνεται,\n",
    "ο αλγόριθμος δεν απορρίπτει αμέσως την αλλαγή αλλά την δέχεται με φθίνουσα πιθανότητα.\n",
    "\n",
    "Στο επόμενο κελί κώδικα δοκιμάζουμε διάφορες θερμοκρασίες και βλέπουμε πως αυτές επηρεάζουν\n",
    "την πιθανότητα αλλαγής σπιν."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temperatures = [1, 2, 5, 10, 20, 50]\n",
    "cmap = plt.cm.get_cmap('Dark2')\n",
    "colors = [cmap(x) for x in np.linspace(0, 1, 6)]\n",
    "\n",
    "for col, temp in zip(colors, temperatures):\n",
    "    y = np.exp(-np.maximum(x,0) / temp)\n",
    "    plt.semilogy(x, y, label='T=%d' % temp)\n",
    "\n",
    "plt.xlabel(\"Energy Difference\")\n",
    "plt.ylabel(\"Transition Probability\")\n",
    "plt.grid(True)\n",
    "plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Παρατηρήστε πώς η θερμοκρασία (`Τ`) \"ελέγχει\" την κλίση του λογαρίθμου της πιθανότητας μετάβασης.\n",
    "Μικρή θερμοκρασία (και άρα μεγάλη κλίση) σημαίνει ότι η αλύσιδα είναι πολύ ευαίσθητη\n",
    "σε αλλαγές της ενέργειας, επομένως οι μόνες αλλαγές που δέχεται είναι αυτές που\n",
    "μειώνουν την ενέργεια του συστήματος ή την αυξάνουν ελάχιστα. Αντίθετα, μεγάλη θερμοκρασία\n",
    "(και άρα μικρή κλίση) σημαίνει ότι η αλυσίδα αδιαφορεί ουσιαστικά για\n",
    "την ενέργεια του συστήματος και δέχεται κάθε προτεινόμενη αλλαγή καθώς όπως βλέπουμε\n",
    "η πιθανότητα αλλαγής σπιν είναι κοντά στη μονάδα.\n",
    "\n",
    "### Ising Animation\n",
    "\n",
    "Για το επόμενο κελί κώδικα χρειάζεται να έχετε στον ίδιο φάκελο με τον notebook το\n",
    "script `animate_ising.py`. Αυτό το script οπτικοποιεί σε \"βίντεο\" την εξέλιξη του\n",
    "πλέγματος των σπιν.\n",
    "\n",
    "Δέχεται 3 ορίσματα\n",
    "\n",
    "1. `-L`: το μέγεθος του πλέγατος\n",
    "2. `-Temp`: τη θερμοκρασία\n",
    "3. `-up`: το ποσοστό των σπιν που είναι αρχικά θετικά\n",
    "\n",
    "Δοκιμάστε να τρέξετε το παρακάτω κελί αφού σβήσετε το `#` από τις πρώτες 2 γραμμές\n",
    "ανάλογα με το λειτουργικό σας σύστημα. \n",
    "\n",
    "Αν δεν τρέξει, ανοίξτε ένα παράθυρο εντολών (Terminal) και από τον φάκελο που έχετε αποθηκεύσει το αρχείο `animate_ising.py` δώστε την εντολή \n",
    "\n",
    "```\n",
    "python3 animate_ising.py -T 0.1 -L 32 -up 0.5`\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if you run Windows/Linux uncomment the next line\n",
    "#%matplotlib qt \n",
    "# if you run Mac OSX uncomment the next line\n",
    "#%matplotlib macosx \n",
    "%run animate_ising.py -T 0.1 -L 32 -up 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Το βίντεο συνεχίζει σε αέναο βρόχο, οπότε για να το σταματήσετε απλά κλείστε το παράθυρο.\n",
    "Αφήστε όμως το πρόγραμμα να τρέξει για αρκετό χρόνο ώστε να δείτε το σύστημα να καταλήγει \n",
    "σε μια κατάσταση με τα σπιν οργανωμένα (είτε σχεδόν όλα +1, είτε σχεδόν όλα -1).\n",
    "Όταν το κλείσετε εκτελέστε το επόμενο κελί κώδικα αν θέλετε τα επόμενα γραφήματα αυτού του\n",
    "notebook να εμφανίζονται μέσα στο notebook (και όχι pop-out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Στο παρακάτω κελί κώδικα σας δίνουμε την συνάρτηση `Ising` αν θελήσετε να την χρησιμοποιείσετε για τις\n",
    "ασκήσεις σας. Για να την τρέξετε χρειάζεται να δώσετε:\n",
    "\n",
    "1. Την αρχική κατάσταση του πλέγματος (ένα numpy array με -1 και 1)\n",
    "2. Τη θερμοκρασία\n",
    "3. (προαιρετικά) τον αριθμό των βημάτων"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def Ising(spins_init, Temp, nsteps=None):\n",
    "    \"\"\"\n",
    "    Ising Model.\n",
    "        spins_init: initial configuration\n",
    "        Temp: the Temperature\n",
    "        updates: the number of updates to perform\n",
    "    \"\"\"\n",
    "    spins = spins_init.astype(np.int8)  # copy and save some space :)\n",
    "    L = spins.shape[0]\n",
    "    N = L * L\n",
    "    kernel = np.array([[0, 1, 0], [1, 0, 1], [0, 1, 0]], dtype=np.int8)\n",
    "    Temp = Temp\n",
    "    if nsteps is None:\n",
    "        nsteps = L * L * 800  # number of updates\n",
    "\n",
    "    for n in range(nsteps):\n",
    "        r, c = np.random.randint(0, L, 2)  # select a vertex at random\n",
    "        C = [(c-1) % L, c % L, (c + 1) % L]\n",
    "        R = [[(r-1) % L], [r % L], [(r+1) % L]]\n",
    "        DH = -spins[r, c] * np.sum(kernel * spins[R, C]) / Temp\n",
    "        if np.log(np.random.rand()) < DH:\n",
    "            spins[r, c] *= -1  # switch sign\n",
    "    \n",
    "    return spins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Παραδοτέα\n",
    "\n",
    "α) Δοκιμάστε να τρέξετε τον κώδικα με όλα τα spin αρχικά στο +1. \n",
    "Τρέξτε τον μερικές φορές. Αλλάζουν τα ποιοτικά χαρακτηριστικά \n",
    "της διαμόρφωσης που επιστρέφει ο αλγόριθμος;\n",
    "\n",
    "β) Αρχίστε τώρα να κατεβάζετε τη θερμοκρασία, θέτοντας διαδοχικά \n",
    "T = 30, 10, 5, 2, 1.5, 1., 0.5, 0.1, 0.01. Δείξτε για κάθε θερμοκρασία την εικόνα που πήρατε.\n",
    "\n",
    "Η μαγνήτιση του συστήματος ορίζεται ως ο μέσος όρος των σπιν του πλέγματος Λ, δηλαδή\n",
    "\n",
    "$$ m(\\sigma) = \\frac{1}{|\\Lambda|} \\sum_{k \\in \\Lambda} \\sigma(k) $$\n",
    "\n",
    "γ) Εξετάστε πώς συμπεριφέρεται η $m(\\sigma)$ για τις καταστάσεις $\\sigma$ \n",
    "που δίνει ο αλγόριθμος καθώς κατεβάζετε τη θερμοκρασία, όπως στο ερώτημα (β).\n",
    "Δείξτε μια γραφική παράσταση της $m(\\sigma)$ ως προς τη θερμοκρασία (για τη διαμόρφωση $\\sigma$ που προσομοιώσατε).\n",
    "\n",
    "δ) Για τη θερμοκρασία `T = 0.1` αλλάξτε την τιμή της παραμέτρου nsteps. \n",
    "Πώς μοιάζουν τα ενδιάμεσα στάδια από τα οποία περνάει η αλυσίδα μέχρι να καταλήξουμε στην εικόνα που πήραμε για `nsteps = 800 × L × L`; \n",
    "Δείξτε τις εικόνες που πήρατε για `nsteps=10,20,50,100,200,400,800  × L × L`"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}