ShortColdAsIce.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Program Submission for the course PHY 312: Numerical Methods, IISER Bhopal.\n",
    "\n",
    "Github Link:\n",
    "https://github.com/AKnightWing/ColdAsIce\n",
    "\n",
    "Date: 9th June, 2020.\n",
    "\n",
    "Team Name: Team ℏ\n",
    "\n",
    "Project Name: ColdAsIce\n",
    "\n",
    "Team Members:\n",
    "Tanmay Bhore\n",
    "Siddharth Chaini\n",
    "Siddharth Bachoti\n",
    "Vipin Khade\n",
    "Vinay Patil\n",
    "\n",
    "This Python3 file contains all the necessary functions to reproduce the results in the accompanying report.\n",
    "\n",
    "'''\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt \n",
    "from matplotlib import colors as clr\n",
    "from matplotlib.lines import Line2D\n",
    "import copy\n",
    "import random\n",
    "import time\n",
    "\n",
    "\n",
    "#SET EPS AND KB VALUE\n",
    "eps = 0.6\n",
    "kb = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialise_state(N): #N is the grid dimension (in the above example, N=4)\n",
    "    '''\n",
    "    A function to initialise a possible state of NxN atoms. Assigns 1 to all bonds.\n",
    "    '''\n",
    "    grid = np.ones((N,N,2),dtype=int)\n",
    "    \n",
    "    return np.array(grid)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_vector(p1,p2,num):\n",
    "    '''\n",
    "    A function to plot any two vectors on matplotlib via arrows\n",
    "    '''\n",
    "    if num>4:\n",
    "        w = 0.12/num\n",
    "        h=0.9/num\n",
    "    else:\n",
    "        w = 0.03\n",
    "        h=0.2\n",
    "    p1 = np.array(p1)\n",
    "    p2 = np.array(p2)\n",
    "    dp = p2-p1\n",
    "    dx = dp[0]\n",
    "    dy = dp[1]\n",
    "\n",
    "    plt.arrow(p1[0], p1[1],\n",
    "                 dx / 2, dy / 2,\n",
    "                 width=w, head_width=h,\n",
    "                 facecolor=\"k\",\n",
    "                 zorder=0)\n",
    "\n",
    "    plt.arrow(p1[0] + dx / 2,\n",
    "                 p1[1] + dy / 2,\n",
    "                 dx / 2, dy / 2,\n",
    "                 width=w, head_width=0,\n",
    "                 facecolor=\"k\",\n",
    "                 zorder=0)\n",
    "\n",
    "def get_coord_list(arr):\n",
    "    '''\n",
    "    A function to get the coordinates of all atoms. Used for visualising our model.\n",
    "    '''\n",
    "    coord_list=[]\n",
    "    num = len(arr)\n",
    "    for i in range(num):\n",
    "        temp_coord = []\n",
    "        for j in range(num):\n",
    "            current_elems = arr[i][j]\n",
    "            xpt = (num-1)-i\n",
    "            ypt = j\n",
    "            temp_coord.append((xpt,ypt))\n",
    "        coord_list.append(temp_coord)\n",
    "    return coord_list\n",
    "\n",
    "def visualise_2d_model(arr,savefig=False,savename=\".temp\",show=True):\n",
    "    '''\n",
    "    Wrapper function to visualise a 2D ice model when input an array on dim (n,n,2)\n",
    "    '''\n",
    "    num = len(arr)\n",
    "    if num<4:\n",
    "        mksize = 35/num \n",
    "    else:\n",
    "        mksize=(35/num)+0.5\n",
    "    plt.axes().set_aspect('equal')\n",
    "    coord_list = get_coord_list(arr)\n",
    "    for i in range(num):\n",
    "        for j in range(num):\n",
    "            current_up_state = arr[i][j][0]\n",
    "            current_right_state = arr[i][j][1]\n",
    "            x_current = coord_list[i][j][1]\n",
    "            y_current = coord_list[i][j][0]\n",
    "\n",
    "            lower_neighbour_up_state = arr[(i+1)%num][j][0]\n",
    "            x_up = coord_list[(i+1)%num][j][1]\n",
    "            y_up = coord_list[(i+1)%num][j][0]\n",
    "\n",
    "            left_neighbour_right_state = arr[i][j-1][1]\n",
    "            x_left = coord_list[i][j-1][1]\n",
    "            y_left = coord_list[i][j-1][0]\n",
    "\n",
    "\n",
    "            current_down_state = -(lower_neighbour_up_state)\n",
    "            current_left_state = -(left_neighbour_right_state)\n",
    "            \n",
    "            plt.plot(x_current,y_current,\n",
    "                        marker=\"o\", markersize=mksize, markeredgecolor=\"k\",\n",
    "                        markerfacecolor=\"yellowgreen\",\n",
    "                        zorder=1)\n",
    "\n",
    "            if current_up_state == 1:\n",
    "                plot_vector([x_current,y_current],[x_current,y_current+1],num)\n",
    "            elif current_up_state == -1:\n",
    "                plot_vector([x_current,y_current+1],[x_current,y_current],num)\n",
    "\n",
    "            if current_right_state == 1:\n",
    "                plot_vector([x_current,y_current],[x_current+1,y_current],num)\n",
    "            elif current_right_state == -1:\n",
    "                plot_vector([x_current+1,y_current],[x_current,y_current],num)\n",
    "\n",
    "            if current_down_state == 1:\n",
    "                plot_vector([x_current,y_current],[x_current,y_current-1],num)\n",
    "            elif current_down_state == -1:\n",
    "                plot_vector([x_current,y_current-1],[x_current,y_current],num)\n",
    "\n",
    "            if current_left_state == 1:\n",
    "                plot_vector([x_current,y_current],[x_current-1,y_current],num)\n",
    "            elif current_left_state == -1:\n",
    "                plot_vector([x_current-1,y_current],[x_current,y_current],num)\n",
    "    plt.xlim(-1,num)\n",
    "    plt.ylim(-1,num)\n",
    "    plt.axis('off')\n",
    "    if savefig:\n",
    "        plt.savefig(f\"{savename}.png\",dpi=300)\n",
    "    if show:\n",
    "        plt.show()\n",
    "    plt.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_6_vertex_configs():\n",
    "    '''\n",
    "    A function to create the 6 images for the types of single vertex configs in the report.\n",
    "    '''\n",
    "    configlist = [(1,1,-1,-1),\n",
    "            (-1,-1,1,1),\n",
    "            (1,-1,-1,1),\n",
    "            (-1,1,1,-1),\n",
    "            (-1,1,-1,1),\n",
    "            (1,-1,1,-1)]\n",
    "    \n",
    "    for i in range(len(configlist)):\n",
    "        config=configlist[i]\n",
    "        plt.axes().set_aspect('equal')\n",
    "        if config[0]==1:\n",
    "            plot_vector([0,0],[0,1],5) #Up\n",
    "        else:\n",
    "            plot_vector([0,1],[0,0],5) #Up\n",
    "        if config[1]==1:\n",
    "            plot_vector([0,0],[1,0],5) #Right\n",
    "        else:\n",
    "            plot_vector([1,0],[0,0],5) #Right\n",
    "        if config[2]==-1:\n",
    "            plot_vector([0,-1],[0,0],5) #Down\n",
    "        else:\n",
    "            plot_vector([0,0],[0,-1],5)\n",
    "        if config[3]==-1:\n",
    "            plot_vector([-1,0],[0,0],5) #Left\n",
    "        else:\n",
    "            plot_vector([0,0],[-1,0],5)\n",
    "        plt.plot(0,0,\n",
    "                    marker=\"o\", markersize=35, markeredgecolor=\"k\",\n",
    "                    markerfacecolor=\"yellowgreen\",\n",
    "                    zorder=1)\n",
    "        plt.xlim(-1,1)\n",
    "        plt.ylim(-1,1)\n",
    "        plt.axis('off')\n",
    "        plt.show()\n",
    "        plt.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def check_config(arr):\n",
    "    '''\n",
    "    Function to check if given state follows the ice rules.\n",
    "    '''\n",
    "    flag=True\n",
    "    N=len(arr)\n",
    "    for i in range(len(arr)):\n",
    "        for j in range(len(arr)):\n",
    "            \n",
    "            current_up_state = arr[i][j][0]\n",
    "            current_right_state = arr[i][j][1]\n",
    "\n",
    "            lower_neighbour_up_state = arr[(i+1)%N][j][0]\n",
    "            \n",
    "            left_neighbour_right_state = arr[i][j-1][1]\n",
    "\n",
    "\n",
    "            current_left_state = -(left_neighbour_right_state)\n",
    "            current_down_state = -(lower_neighbour_up_state)\n",
    "            \n",
    "            if (current_up_state + current_right_state + current_left_state + current_down_state) != 0:\n",
    "                flag=False\n",
    "                break\n",
    "    return flag\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def long_loop(arr2, verbose=False):\n",
    "    '''\n",
    "    Implementation of the long loop algorithm\n",
    "    '''\n",
    "    arr = copy.deepcopy(arr2)\n",
    "    N=len(arr)\n",
    "    iters=0\n",
    "    \n",
    "    n1 = np.random.randint(low=0, high=N)\n",
    "    n2 = np.random.randint(low=0, high=N)\n",
    "    inital_pt =(n1,n2)\n",
    "    prev_choice=None\n",
    "    \n",
    "    while True:\n",
    "        iters+=1\n",
    "        if n1==inital_pt[0] and n2==inital_pt[1] and iters!=1:\n",
    "            if verbose:\n",
    "                print(f\"Completed in {iters} iterations.\")\n",
    "             #assert(check_config(arr))\n",
    "            break\n",
    "        current_up_state = arr[n1][n2][0]\n",
    "        current_right_state = arr[n1][n2][1]\n",
    "\n",
    "        lower_neighbour_up_state = arr[(n1+1)%N][n2][0]\n",
    "        left_neighbour_right_state = arr[n1][n2-1][1]\n",
    "\n",
    "        current_down_state = -(lower_neighbour_up_state)    \n",
    "        current_left_state = -(left_neighbour_right_state)\n",
    "\n",
    "        current_states_dict = {\"up\":current_up_state,\"right\":current_right_state,\"down\":current_down_state,\"left\":current_left_state}\n",
    "        outgoing_state_dict={}\n",
    "        incoming_state_dict={}\n",
    "\n",
    "        for key in current_states_dict.keys():\n",
    "            if current_states_dict[key]==1:  #current state is outgoing\n",
    "                outgoing_state_dict[key]=current_states_dict[key]\n",
    "            else:\n",
    "                incoming_state_dict[key]=current_states_dict[key]\n",
    "\n",
    "        if prev_choice ==\"right\":\n",
    "            forbidden_choice=\"left\"\n",
    "        elif prev_choice ==\"up\":\n",
    "            forbidden_choice=\"down\"\n",
    "        elif prev_choice ==\"left\":\n",
    "            forbidden_choice=\"right\"\n",
    "        elif prev_choice ==\"down\":\n",
    "            forbidden_choice=\"up\"\n",
    "        else:\n",
    "            forbidden_choice=None\n",
    "\n",
    "\n",
    "        while True:\n",
    "            out_choice = np.random.choice(list(outgoing_state_dict.keys()))\n",
    "            if out_choice !=forbidden_choice:\n",
    "                break\n",
    "\n",
    "        prev_choice=out_choice\n",
    "\n",
    "        if out_choice == \"up\":\n",
    "            arr[n1][n2][0]= - (arr[n1][n2][0])\n",
    "            n1=(n1-1)%N\n",
    "            n2=n2\n",
    "            continue\n",
    "\n",
    "        if out_choice == \"right\":\n",
    "            arr[n1][n2][1]= - (arr[n1][n2][1])\n",
    "            n1=n1\n",
    "            n2=(n2+1)%N\n",
    "            continue\n",
    "\n",
    "        if out_choice == \"down\":\n",
    "            arr[(n1+1)%N][n2][0]= - (arr[(n1+1)%N][n2][0])\n",
    "            n1=(n1+1)%N\n",
    "            n2=n2\n",
    "            continue\n",
    "\n",
    "        if out_choice == \"left\":\n",
    "            arr[n1][(n2-1)%N][1]= - (arr[n1][(n2-1)%N][1])\n",
    "            n1=n1\n",
    "            n2=(n2-1)%N\n",
    "            continue\n",
    "\n",
    "\n",
    "\n",
    "    return arr\n",
    "\n",
    "def short_loop(arr2, verbose=False):\n",
    "    '''\n",
    "    Implementation of the short loop algorithm\n",
    "    '''\n",
    "    arr = copy.deepcopy(arr2)\n",
    "    N=len(arr)\n",
    "    iters=0\n",
    "    \n",
    "    n1 = np.random.randint(low=0, high=N)\n",
    "    n2 = np.random.randint(low=0, high=N)\n",
    "    inital_pts =[]\n",
    "    prev_choice=None\n",
    "    forbidden_choice = \"None\"\n",
    "    while True:\n",
    "        iters+=1\n",
    "        if iters !=1 and (n1,n2) in inital_pts:\n",
    "            idx=inital_pts.index((n1,n2))\n",
    "            reset_pts=inital_pts[:idx]\n",
    "            for thispt in reset_pts:\n",
    "                arr[thispt[0]][thispt[1]][0]=arr2[thispt[0]][thispt[1]][0]\n",
    "                arr[thispt[0]][thispt[1]][1]=arr2[thispt[0]][thispt[1]][1]\n",
    "            if verbose:\n",
    "                print(f\"Completed in {iters} iterations.\")\n",
    "            break\n",
    "        \n",
    "        current_up_state = arr[n1][n2][0]\n",
    "        current_right_state = arr[n1][n2][1]\n",
    "\n",
    "        lower_neighbour_up_state = arr[(n1+1)%N][n2][0]\n",
    "        left_neighbour_right_state = arr[n1][n2-1][1]\n",
    "\n",
    "        current_down_state = -(lower_neighbour_up_state)    \n",
    "        current_left_state = -(left_neighbour_right_state)\n",
    "\n",
    "        current_states_dict = {\"up\":current_up_state,\"right\":current_right_state,\"down\":current_down_state,\"left\":current_left_state}\n",
    "        outgoing_state_dict={}\n",
    "        incoming_state_dict={}\n",
    "\n",
    "        for key in current_states_dict.keys():\n",
    "            if current_states_dict[key]==1:  #current state is outgoing\n",
    "                outgoing_state_dict[key]=current_states_dict[key]\n",
    "            else:\n",
    "                incoming_state_dict[key]=current_states_dict[key]\n",
    "        if prev_choice ==\"right\":\n",
    "            forbidden_choice=\"left\"\n",
    "        elif prev_choice ==\"up\":\n",
    "            forbidden_choice=\"down\"\n",
    "        elif prev_choice ==\"left\":\n",
    "            forbidden_choice=\"right\"\n",
    "        elif prev_choice ==\"down\":\n",
    "            forbidden_choice=\"up\"\n",
    "        else:\n",
    "            forbidden_choice=None\n",
    "        while True:\n",
    "            out_choice = np.random.choice(list(outgoing_state_dict.keys()))\n",
    "            if out_choice !=forbidden_choice:\n",
    "                break\n",
    "        prev_choice=out_choice\n",
    "        \n",
    "        inital_pts.append((n1,n2))     #append old n1,n2\n",
    "        \n",
    "        if out_choice == \"up\":\n",
    "            arr[n1][n2][0]= - (arr[n1][n2][0])\n",
    "            n1=(n1-1)%N\n",
    "            n2=n2\n",
    "            continue\n",
    "\n",
    "        if out_choice == \"right\":\n",
    "            arr[n1][n2][1]= - (arr[n1][n2][1])\n",
    "            n1=n1\n",
    "            n2=(n2+1)%N\n",
    "            continue\n",
    "\n",
    "        if out_choice == \"down\":\n",
    "            arr[(n1+1)%N][n2][0]= - (arr[(n1+1)%N][n2][0])\n",
    "            n1=(n1+1)%N\n",
    "            n2=n2\n",
    "            continue\n",
    "\n",
    "        if out_choice == \"left\":\n",
    "            arr[n1][(n2-1)%N][1]= - (arr[n1][(n2-1)%N][1])\n",
    "            n1=n1\n",
    "            n2=(n2-1)%N\n",
    "            continue\n",
    "        \n",
    "    return arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def count_states(num,error_threshold,return_dict = False,verbose=False):\n",
    "    '''\n",
    "    Function to count total possible states for num x num atoms.\n",
    "    '''\n",
    "    if not (error_threshold<=100 and error_threshold>0):\n",
    "        print(\"Error! Please input error_threshold as a value between 0 and 100\")\n",
    "        assert (error_threshold<=100 and error_threshold>0)\n",
    "    state_dict={}\n",
    "    oldarr = short_loop(initialise_state(num), verbose=False)\n",
    "    good_iterations = 0           #Iterations that gave us a new state, so good.\n",
    "    bad_iterations = 0               #Iterations that gave us an already found state,so a waste and hence bad.\n",
    "    while True:\n",
    "        newarr = short_loop(oldarr,verbose=False)\n",
    "        name =arr_to_string(newarr)\n",
    "        if name not in state_dict:\n",
    "            count_repetitions=0\n",
    "            state_dict[name]=1\n",
    "            good_iterations+=1\n",
    "        else:\n",
    "            bad_iterations+=1\n",
    "            count_repetitions+=1\n",
    "            state_dict[name]+=1\n",
    "        percent_approx_err=good_iterations*100/(good_iterations+bad_iterations)\n",
    "        if verbose:\n",
    "            print(f\"Good iterations = {good_iterations} and bad iterations = {bad_iterations} and Error % = {percent_approx_err}\", end=\"\\r\",flush=True)\n",
    "        if percent_approx_err < error_threshold:\n",
    "            break\n",
    "        oldarr=newarr\n",
    "    if return_dict:\n",
    "        return len(state_dict),state_dict\n",
    "    else:\n",
    "        return len(state_dict)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def state2to4(arr):\n",
    "    '''\n",
    "    Convert an (n,n,2) array to an (n,n,4) array.\n",
    "    '''\n",
    "    fourstatearr=np.zeros((arr.shape[0],arr.shape[1],4))\n",
    "    N=len(arr)\n",
    "    for i in range(len(arr)):\n",
    "        for j in range(len(arr)):\n",
    "            current_up_state = arr[i][j][0]\n",
    "            current_right_state = arr[i][j][1]\n",
    "            lower_neighbour_up_state = arr[(i+1)%N][j][0]\n",
    "            left_neighbour_right_state = arr[i][j-1][1]\n",
    "            current_left_state = -(left_neighbour_right_state)\n",
    "            current_down_state = -(lower_neighbour_up_state)\n",
    "            fourstatearr[i][j][0] = current_up_state\n",
    "            fourstatearr[i][j][1] = current_right_state\n",
    "            fourstatearr[i][j][2] = current_down_state\n",
    "            fourstatearr[i][j][3] = current_left_state\n",
    "\n",
    "    return fourstatearr\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rot90_anticlock(arr2):\n",
    "    '''\n",
    "    Function to generate 90deg anticlockwise rotation\n",
    "    Up becomes left, left becomes down, down becomes right, right becomes up.\n",
    "    '''\n",
    "    fourstatearr = state2to4(arr2)\n",
    "    fourstatearr = np.rot90(fourstatearr,1)\n",
    "    arr=np.zeros((fourstatearr.shape[0],fourstatearr.shape[1],2))\n",
    "    N=len(arr)\n",
    "    for i in range(len(arr)):\n",
    "        for j in range(len(arr)):\n",
    "            current_up_state = fourstatearr[i][j][0]\n",
    "            current_right_state = fourstatearr[i][j][1]\n",
    "            current_down_state = fourstatearr[i][j][2]\n",
    "            current_left_state = fourstatearr[i][j][3]\n",
    "            \n",
    "            new_up_state = current_right_state\n",
    "            new_right_state = current_down_state\n",
    "            \n",
    "            arr[i][j][0]=new_up_state\n",
    "            arr[i][j][1]=new_right_state\n",
    "    \n",
    "    return arr.astype(int)\n",
    "\n",
    "def rot180_anticlock(arr2):\n",
    "    '''\n",
    "    Author: Siddharth Chaini\n",
    "    Function to generate 180deg anticlockwise rotation\n",
    "    Up becomes down, left becomes right, down becomes up, right becomes left.\n",
    "    '''\n",
    "    fourstatearr = state2to4(arr2)\n",
    "    fourstatearr = np.rot90(fourstatearr,2)\n",
    "    arr=np.zeros((fourstatearr.shape[0],fourstatearr.shape[1],2))\n",
    "    N=len(arr)\n",
    "    for i in range(len(arr)):\n",
    "        for j in range(len(arr)):\n",
    "            current_up_state = fourstatearr[i][j][0]\n",
    "            current_right_state = fourstatearr[i][j][1]\n",
    "            current_down_state = fourstatearr[i][j][2]\n",
    "            current_left_state = fourstatearr[i][j][3]\n",
    "            \n",
    "            new_up_state = current_down_state\n",
    "            new_right_state = current_left_state\n",
    "            \n",
    "            arr[i][j][0]=new_up_state\n",
    "            arr[i][j][1]=new_right_state\n",
    "    \n",
    "    return arr.astype(int)\n",
    "\n",
    "def rot270_anticlock(arr2):\n",
    "    '''\n",
    "    Function to generate 270deg anticlockwise rotation\n",
    "\tUp becomes right, left becomes up, down becomes left, right becomes down\n",
    "    '''\n",
    "    fourstatearr = state2to4(arr2)\n",
    "    fourstatearr = np.rot90(fourstatearr,3)\n",
    "    arr=np.zeros((fourstatearr.shape[0],fourstatearr.shape[1],2))\n",
    "    N=len(arr)\n",
    "    for i in range(len(arr)):\n",
    "        for j in range(len(arr)):\n",
    "            current_up_state = fourstatearr[i][j][0]\n",
    "            current_right_state = fourstatearr[i][j][1]\n",
    "            current_down_state = fourstatearr[i][j][2]\n",
    "            current_left_state = fourstatearr[i][j][3]\n",
    "            \n",
    "            new_up_state = current_left_state\n",
    "            new_right_state = current_up_state\n",
    "            \n",
    "            arr[i][j][0]=new_up_state\n",
    "            arr[i][j][1]=new_right_state\n",
    "    \n",
    "    return arr.astype(int)\n",
    "\n",
    "def hor_flip(arr2):\n",
    "    '''\n",
    "    Function to generate a horizontally flipped configuration\n",
    "    '''\n",
    "    arr = np.flip(arr2,1)\n",
    "    proper_arr=np.zeros_like(arr2)\n",
    "    num = len(arr)\n",
    "    for i in range(num):\n",
    "        for j in range(num):\n",
    "            current_up_state = arr[i][j][0]\n",
    "            current_left_state = arr[i][j][1]\n",
    "            \n",
    "            right_neighbour_left_state = arr[i][(j+1)%num][1]\n",
    "            current_right_state = - (right_neighbour_left_state)\n",
    "\n",
    "            proper_arr[i][j][0]=current_up_state\n",
    "            proper_arr[i][j][1]=current_right_state\n",
    "    return proper_arr.astype(int)\n",
    "\n",
    "def ver_flip(arr2):\n",
    "    '''\n",
    "    Function to generate a vertically flipped configuration\n",
    "    '''\n",
    "    arr = np.flip(arr2,0)\n",
    "    proper_arr=np.zeros_like(arr2)\n",
    "    num = len(arr)\n",
    "    for i in range(num):\n",
    "        for j in range(num):\n",
    "            current_down_state = arr[i][j][0]\n",
    "            current_right_state = arr[i][j][1]\n",
    "            \n",
    "            upper_neighbour_down_state = arr[i-1][j][0]\n",
    "            current_up_state = - (upper_neighbour_down_state)\n",
    "            \n",
    "            proper_arr[i][j][0]=current_up_state\n",
    "            proper_arr[i][j][1]=current_right_state\n",
    "    return proper_arr.astype(int)\n",
    "\n",
    "def flip_secondary_diag(arr2):\n",
    "    '''\n",
    "    Function to generate a secondary diagonal flipped configuration.\n",
    "    '''\n",
    "    arr = copy.deepcopy(arr2)\n",
    "    N = len(arr)\n",
    "    for i in range(N):\n",
    "        for j in range(N):\n",
    "            if (i+j)<=N-1:\n",
    "                dist = N-(i+j+1)\n",
    "                arr[i][j][0], arr[i+dist][j+dist][0], arr[i][j][1], arr[i+dist][j+dist][1] = arr[i+dist][j+dist][1], arr[i][j][1], arr[i+dist][j+dist][0], arr[i][j][0]                \n",
    "    return arr.astype(int)\n",
    "\n",
    "def flip_primary_diag(arr2):\n",
    "    '''\n",
    "    Function to generate a primary diagonal flipped configuration.\n",
    "    '''\n",
    "    arr = copy.deepcopy(arr2)\n",
    "    N = len(arr)\n",
    "    arr = rot90_anticlock(flip_secondary_diag(rot270_anticlock(arr)))\n",
    "    return arr.astype(int)\n",
    "\n",
    "def get_all_column_translations(arr):\n",
    "    '''\n",
    "    Return a list of arrays containing all possible column translations.\n",
    "    '''\n",
    "    result_arr_list=[]\n",
    "    N=len(arr) \n",
    "    for i in range(1,N):\n",
    "        a1 = arr[:,0:i].reshape(N,-1,2)\n",
    "        a2 = arr[:,i:].reshape(N,-1,2)\n",
    "        res = np.hstack([a2,a1])\n",
    "        result_arr_list.append(res)\n",
    "    return result_arr_list\n",
    "\n",
    "def get_all_row_translations(arr):\n",
    "    '''\n",
    "    Return a list of arrays containing all possible column translations.\n",
    "    '''\n",
    "    result_arr_list=[]\n",
    "    N=len(arr) \n",
    "    for i in range(1,N):\n",
    "        a1 = arr[0:i,:].reshape(-1,N,2)\n",
    "        a2 = arr[i:,:].reshape(-1,N,2)\n",
    "        res = np.vstack([a2,a1])\n",
    "        result_arr_list.append(res)\n",
    "    return result_arr_list\n",
    "\n",
    "def arr_to_string(arr2):\n",
    "    '''\n",
    "=   Convert a state array to a unique string.\n",
    "    '''\n",
    "    arr = copy.deepcopy(arr2)\n",
    "    name = ' '.join(map(str, arr.flatten())).replace(' ','')\n",
    "    return name\n",
    "\n",
    "def string_to_arr(s):\n",
    "    '''\n",
    "    Convert the above string back to an array.\n",
    "    '''\n",
    "    replaced_str = s.replace(\"-1\",\"0\")\n",
    "    arr=[]\n",
    "    for i in replaced_str:\n",
    "        if i=='1':\n",
    "            arr.append(1)\n",
    "        elif i==\"0\":\n",
    "            arr.append(-1)\n",
    "        else:\n",
    "            print(\"ERROR\")\n",
    "            assert(1==0)\n",
    "    arr = np.array(arr)\n",
    "    arr = arr.reshape(int(np.sqrt(len(arr)/2)),int(np.sqrt(len(arr)/2)),2)\n",
    "    return arr\n",
    "\n",
    "def remove_symmetries(all_names, verbose=False):\n",
    "    '''\n",
    "    Function to remove all symmetries for a given list of state strings (names)\n",
    "    '''\n",
    "    assert type(all_names)==list\n",
    "    for i,given_name in enumerate(all_names):\n",
    "        if verbose and i%100 ==0:\n",
    "            print(f\"Loading... {i}/{len(all_names)} done. Percent Completed = {100*i/len(all_names)}.\", end=\"\\r\",flush=True)\n",
    "        arr = string_to_arr(given_name)\n",
    "        #Column Translation symmetries\n",
    "        templist=get_all_column_translations(arr)\n",
    "        for newarr in templist:\n",
    "            name = arr_to_string(newarr)\n",
    "            if name in all_names[i+1:]:\n",
    "                idx = all_names[i+1:].index(name) + i+1\n",
    "                del all_names[idx]\n",
    "                i = all_names.index(given_name)\n",
    "        \n",
    "        #Row Translation symmetries\n",
    "        templist=get_all_row_translations(arr)\n",
    "        for newarr in templist:\n",
    "            name = arr_to_string(newarr)\n",
    "            if name in all_names[i+1:]:\n",
    "                idx = all_names[i+1:].index(name) + i+1\n",
    "                del all_names[idx]\n",
    "                i = all_names.index(given_name)\n",
    "        \n",
    "        #Check 90 degree rotation symmetry\n",
    "        name = arr_to_string(rot90_anticlock(arr))\n",
    "        if name in all_names[i+1:]:\n",
    "            idx = all_names[i+1:].index(name) + i+1\n",
    "            del all_names[idx]\n",
    "            i = all_names.index(given_name)\n",
    "\n",
    "        #Check 180 degree rotation symmetry\n",
    "        name = arr_to_string(rot180_anticlock(arr))\n",
    "        if name in all_names[i+1:]:\n",
    "            idx = all_names[i+1:].index(name) + i+1\n",
    "            del all_names[idx]\n",
    "            i = all_names.index(given_name)\n",
    "\n",
    "        #Check 270 degree rotation symmetry\n",
    "        name = arr_to_string(rot270_anticlock(arr))\n",
    "        if name in all_names[i+1:]:\n",
    "            idx = all_names[i+1:].index(name) + i+1\n",
    "            del all_names[idx]\n",
    "            i = all_names.index(given_name)\n",
    "\n",
    "        #Check horizontal flip symmetry\n",
    "        name = arr_to_string(hor_flip(arr))\n",
    "        if name in all_names[i+1:]:\n",
    "            idx = all_names[i+1:].index(name) + i+1\n",
    "            del all_names[idx]\n",
    "            i = all_names.index(given_name)\n",
    "\n",
    "        #Check vertical flip symmetry\n",
    "        name = arr_to_string(ver_flip(arr))\n",
    "        if name in all_names[i+1:]:\n",
    "            idx = all_names[i+1:].index(name) + i+1\n",
    "            del all_names[idx]\n",
    "            i = all_names.index(given_name)\n",
    "\n",
    "        #Check secondary diagonal flip symmetry\n",
    "        name = arr_to_string(flip_secondary_diag(arr))\n",
    "        if name in all_names[i+1:]:\n",
    "            idx = all_names[i+1:].index(name) + i+1\n",
    "            del all_names[idx]\n",
    "            i = all_names.index(given_name)\n",
    "\n",
    "        #Check primary diagonal flip symmetry\n",
    "        name = arr_to_string(flip_primary_diag(arr))\n",
    "        if name in all_names[i+1:]:\n",
    "            idx = all_names[i+1:].index(name) + i+1\n",
    "            del all_names[idx]\n",
    "            i = all_names.index(given_name)\n",
    "    return all_names\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_energy(arr,eps):\n",
    "    '''\n",
    "    A function to calculate the energy for a given state array and epsilon.\n",
    "    '''\n",
    "    N = len(arr)\n",
    "    s = 0\n",
    "    for i in range(len(arr)):\n",
    "        for j in range(len(arr)):  \n",
    "            current_up_state = arr[i][j][0]\n",
    "            current_right_state = arr[i][j][1]\n",
    "            lower_neighbour_up_state = arr[(i+1)%N][j][0] \n",
    "            left_neighbour_right_state = arr[i][j-1][1]\n",
    "            if  current_up_state ==1 and lower_neighbour_up_state == -1 : \n",
    "                s+= 1\n",
    "            elif current_right_state == 1 and left_neighbour_right_state == -1:\n",
    "                s+= 1\n",
    "    return (-s*eps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_atom_config(arr):\n",
    "    '''\n",
    "    A function to represent all atoms as a combination of the 6 possible vertex configurations.\n",
    "    '''\n",
    "    tupdict = {\n",
    "        (1,1,-1,-1):1,\n",
    "        (-1,-1,1,1):2,\n",
    "        (1,-1,-1,1):3,\n",
    "        (-1,1,1,-1):4,\n",
    "        (-1,1,-1,1):5,\n",
    "        (1,-1,1,-1):6\n",
    "    }\n",
    "    num=len(arr)\n",
    "    typearr = np.zeros((arr.shape[0],arr.shape[1]))\n",
    "    for i in range(len(arr)):\n",
    "        for j in range(len(arr)):\n",
    "            current_up_state = arr[i][j][0]\n",
    "            current_right_state = arr[i][j][1]\n",
    "\n",
    "            lower_neighbour_up_state = arr[(i+1)%num][j][0]\n",
    "\n",
    "            left_neighbour_right_state = arr[i][j-1][1]\n",
    "\n",
    "            current_down_state = -(lower_neighbour_up_state)\n",
    "            current_left_state = -(left_neighbour_right_state)\n",
    "\n",
    "            tup = (current_up_state,current_right_state,current_down_state,current_left_state)\n",
    "            atomtype = tupdict[tup]\n",
    "            typearr[i][j] = atomtype\n",
    "    return typearr.astype(int)\n",
    "\n",
    "def visualise_atom_config_with_bonds(arr,savefig=False,savename=\".temp\",show=True):\n",
    "    '''\n",
    "    A function to visualise the state array as a combination of the 6 possible vertex configurations. \n",
    "    '''\n",
    "    num = len(arr)\n",
    "    plt.axes().set_aspect('equal')\n",
    "    coord_list = get_coord_list(arr)\n",
    "    atom_config_arr = calculate_atom_config(arr)\n",
    "    for i in range(num):\n",
    "        for j in range(num):\n",
    "            current_up_state = arr[i][j][0]\n",
    "            current_right_state = arr[i][j][1]\n",
    "            x_current = coord_list[i][j][1]\n",
    "            y_current = coord_list[i][j][0]\n",
    "\n",
    "            lower_neighbour_up_state = arr[(i+1)%num][j][0]\n",
    "            x_up = coord_list[(i+1)%num][j][1]\n",
    "            y_up = coord_list[(i+1)%num][j][0]\n",
    "\n",
    "            left_neighbour_right_state = arr[i][j-1][1]\n",
    "            x_left = coord_list[i][j-1][1]\n",
    "            y_left = coord_list[i][j-1][0]\n",
    "\n",
    "\n",
    "            current_down_state = -(lower_neighbour_up_state)\n",
    "            current_left_state = -(left_neighbour_right_state)\n",
    "\n",
    "            if atom_config_arr[i][j]==5:\n",
    "                plt.plot(x_current,y_current,\n",
    "                            marker=\"o\", markersize=9, markeredgecolor=\"k\",\n",
    "                            markerfacecolor=\"cyan\",\n",
    "                            zorder=1)\n",
    "            elif atom_config_arr[i][j]==6:\n",
    "                plt.plot(x_current,y_current,\n",
    "                            marker=\"o\", markersize=9, markeredgecolor=\"k\",\n",
    "                            markerfacecolor=\"magenta\",\n",
    "                            zorder=1)\n",
    "            else:\n",
    "                plt.plot(x_current,y_current,\n",
    "                            marker=\"o\", markersize=9, markeredgecolor=\"k\",\n",
    "                            markerfacecolor=\"gray\",\n",
    "                            zorder=1)                \n",
    "            if current_up_state == 1:\n",
    "                plot_vector([x_current,y_current],[x_current,y_current+1],num)\n",
    "            elif current_up_state == -1:\n",
    "                plot_vector([x_current,y_current+1],[x_current,y_current],num)\n",
    "\n",
    "            if current_right_state == 1:\n",
    "                plot_vector([x_current,y_current],[x_current+1,y_current],num)\n",
    "            elif current_right_state == -1:\n",
    "                plot_vector([x_current+1,y_current],[x_current,y_current],num)\n",
    "\n",
    "            if current_down_state == 1:\n",
    "                plot_vector([x_current,y_current],[x_current,y_current-1],num)\n",
    "            elif current_down_state == -1:\n",
    "                plot_vector([x_current,y_current-1],[x_current,y_current],num)\n",
    "\n",
    "            if current_left_state == 1:\n",
    "                plot_vector([x_current,y_current],[x_current-1,y_current],num)\n",
    "            elif current_left_state == -1:\n",
    "                plot_vector([x_current-1,y_current],[x_current,y_current],num)\n",
    "    plt.xlim(-1,num)\n",
    "    plt.ylim(-1,num)\n",
    "    plt.axis('off')\n",
    "    if savefig:\n",
    "        plt.savefig(f\"{savename}.png\",dpi=300)\n",
    "    if show:\n",
    "        plt.show()\n",
    "    plt.close()\n",
    "\n",
    "def visualise_final_state(arr,title=\"\",savefig=False,savename=\".temp\",show=True):\n",
    "    '''\n",
    "    A function to generate a checkerboard plot (as given in the report) for a given state array.\n",
    "    '''\n",
    "\n",
    "    atom_config_arr = calculate_atom_config(arr)\n",
    "    colordict = {1:0,2:0,3:0,4:0,5:-1,6:1}\n",
    "    atom_config_arr_copy = np.copy(atom_config_arr)\n",
    "    for key, val in colordict.items(): atom_config_arr_copy[atom_config_arr==key] = val\n",
    "    flipped_arr = np.flip(atom_config_arr_copy,0)\n",
    "    \n",
    "    fig = plt.figure()\n",
    "    ax = plt.subplot(111)\n",
    "    \n",
    "    ax.set_aspect('equal')\n",
    "    cmap = plt.cm.coolwarm\n",
    "    norm = clr.Normalize(vmin=-1, vmax=1)\n",
    "    hexlist = [clr.to_hex(cmap(norm(-1))),clr.to_hex(cmap(norm(0))),clr.to_hex(cmap(norm(1)))]\n",
    "\n",
    "    legend_elements = [\n",
    "                    Line2D([0], [0], marker='o', color='w', label='Type 1/2/3/4',markerfacecolor=hexlist[1], markersize=10),\n",
    "                    Line2D([0], [0], marker='o', color='w', label='Type 5',markerfacecolor=hexlist[0], markersize=10),\n",
    "                    Line2D([0], [0], marker='o', color='w', label='Type 6',markerfacecolor=hexlist[2], markersize=10)\n",
    "                        ]\n",
    "    ax.pcolormesh(flipped_arr,cmap=cmap,norm=norm)\n",
    "\n",
    "    box = ax.get_position()\n",
    "    ax.set_position([box.x0, box.y0, box.width * 1, box.height])\n",
    "\n",
    "    ax.legend(handles=legend_elements,loc='center left', bbox_to_anchor=(1, 0.5))\n",
    "\n",
    "    plt.title(title)\n",
    "    if savefig:\n",
    "        plt.savefig(f\"{savename}.png\",dpi=300)\n",
    "    if show:\n",
    "        plt.show()\n",
    "    plt.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def metropolis_move(arr,temp,verbose = False):\n",
    "    '''\n",
    "    A function to make the metropolis move as part of the F model.\n",
    "    Dependencies: Uses energy(arr) function\n",
    "    and short_loop(arr) function.\n",
    "    \n",
    "    '''\n",
    "    new = short_loop(arr, verbose = verbose)\n",
    "    old_energy = calculate_energy(arr,eps)\n",
    "    new_energy = calculate_energy(new,eps)\n",
    "    delta_E = new_energy - old_energy\n",
    "\n",
    "    if temp >0:\n",
    "        beta = 1/(temp*kb)\n",
    "        if delta_E <0:\n",
    "            arr = new\n",
    "        elif np.random.uniform() < np.exp(-beta*delta_E):\n",
    "            arr = new\n",
    "    elif temp==0:\n",
    "        if delta_E <0:\n",
    "            arr = new\n",
    "    return arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_vector_polarization(arr):\n",
    "    '''\n",
    "    A function to calculate the polarization vector for a given state array.\n",
    "    '''\n",
    "    N = len(arr)**2\n",
    "    vert = arr[:,:,0].sum()\n",
    "    hor = arr[:,:,-1].sum()\n",
    "    polarization= ((1/((2**0.5)*N))*hor,(1/((2**0.5)*N))*vert)\n",
    "    return polarization[0],polarization[1]\n",
    "\n",
    "def calculate_polarization(arr):\n",
    "    '''\n",
    "    A function to calculate the polarization value for a given state array.\n",
    "    '''\n",
    "    polarization=calculate_vector_polarization(arr)\n",
    "    polval = np.sqrt(polarization[0]**2 + polarization[1]**2)\n",
    "    return polval\n",
    "\n",
    "\n",
    "def calculate_old_polarization(arr):\n",
    "    '''\n",
    "    DEPRECATED: 12th June 2020\n",
    "    A function to calculate the polarization for a given state array.\n",
    "    '''\n",
    "    vert = 0\n",
    "    hor = 0\n",
    "    N = len(arr)**2\n",
    "    for i in arr:\n",
    "        for j in i:\n",
    "            vert+=j[0]\n",
    "            hor+=j[1]\n",
    "\n",
    "    polarization= ((1/((2**0.5)*N))*hor,(1/((2**0.5)*N))*vert)\n",
    "    \n",
    "    polval = np.sqrt(polarization[0]**2 + polarization[1]**2)\n",
    "    \n",
    "    return polval\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def equilibrate(n,temp,breaking_iters=3000,return_polarlist=False,verbose=False):\n",
    "    '''\n",
    "    A function to equilibrate any array of size n x n, i.e use the metropolis move till it has plateaued.\n",
    "    '''\n",
    "    assert temp>=0\n",
    "    arr = initialise_state(n)\n",
    "    polar_list=[]\n",
    "    iter_list=[]\n",
    "    polval = calculate_polarization(arr)\n",
    "    for i in range(breaking_iters):\n",
    "        if verbose:\n",
    "            print(f\"Pol = {polval} and iter = {i}\",end=\"\\r\",flush=True)\n",
    "        arr=metropolis_move(arr,temp)\n",
    "        newpolval = calculate_polarization(arr)\n",
    "        polar_list.append(polval)\n",
    "        polval=newpolval\n",
    "    if return_polarlist:\n",
    "        return polar_list,arr\n",
    "    else:\n",
    "        return arr\n",
    "\n",
    "def equilibrate_existing(arr2,temp,breaking_iters=3000,return_polarlist=False,verbose=False):\n",
    "    '''\n",
    "    A function to equilibrate a given state array, i.e use the metropolis move till it has plateaued.\n",
    "    '''\n",
    "    assert temp>=0\n",
    "    arr = copy.deepcopy(arr2)\n",
    "    n=len(arr)\n",
    "    polar_list=[]\n",
    "    iter_list=[]\n",
    "    polval = calculate_polarization(arr)\n",
    "    for i in range(breaking_iters):\n",
    "        if verbose:\n",
    "            print(f\"Pol = {polval} and iter = {i}\",end=\"\\r\",flush=True)\n",
    "        arr=metropolis_move(arr,temp)\n",
    "        newpolval = calculate_polarization(arr)\n",
    "        polar_list.append(polval)\n",
    "        polval=newpolval\n",
    "    if return_polarlist:\n",
    "        return polar_list,arr\n",
    "    else:\n",
    "        return arr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_sp_heat_and_energy(n,temp,breaking_iters=3000):\n",
    "    '''\n",
    "    A function to calculate the specific heat, energy and polar.\n",
    "    '''\n",
    "    arr = equilibrate(n,temp,breaking_iters=breaking_iters)\n",
    "    #polar_list=[]\n",
    "    energy_list=[]\n",
    "    for i in range(100):\n",
    "        arr=metropolis_move(arr,temp)\n",
    "        en=calculate_energy(arr,eps)\n",
    "        #pol=calculate_polarization(arr)\n",
    "        energy_list.append(en)\n",
    "        #polar_list.append(pol)\n",
    "    energy_list = np.array(energy_list)\n",
    "    #polar_list = np.array(polar_list)\n",
    "    energy = energy_list.mean()\n",
    "    #pol = polar_list.mean()\n",
    "    beta = 1/(temp*kb)\n",
    "    sp_heat = beta**2 *kb* ((energy_list**2).mean() - (energy_list.mean())**2)/n**2\n",
    "    return sp_heat, energy\n",
    "    #return sp_heat, energy, pol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "def save_equilibration_snapshots(n,temp_min=0.001,temp_max=15,n_points=50,verbose=False):\n",
    "\t'''\n",
    "\tA function to save the equilibration snapshots as pngs.\n",
    "\t'''\n",
    "\ttemp_arr=np.linspace(temp_min,temp_max,n_points)\n",
    "\tif verbose:\n",
    "\t\tprint(f\" TEMP = {temp} *******\")\n",
    "\tfor temp in temp_arr:\n",
    "\t\tarr2 = equilibrate(n,temp=temp)        \n",
    "\t\ttitle=f\"T = {temp} K\"\n",
    "\t\tvisualise_final_state(arr2,title=title,savefig=True,savename=title)\n",
    "\n",
    "def generate_gif(input_dir,output_dir,framerate):\n",
    "    '''\n",
    "    A function that creates a GIF showing the equilibriation of a lattice as the temperature is decreased.\n",
    "    It uses the images generated by the function 'save_equilibriation_snapshots()'.\n",
    "    '''\n",
    "    frames = []\n",
    "    dirname = input_dir\n",
    "    temperatures = []\n",
    "    for i in os.listdir(dirname):\n",
    "        t = (i.split('T = '))[1].split(' K.png')[0]\n",
    "        temperatures.append(float(t))\n",
    "    temperatures.sort(reverse=True)#Make reverse flase in order to show the animation from 0->15K\n",
    "\n",
    "    filenames = []\n",
    "    for j in temperatures:\n",
    "        filenames.append(dirname + '/' + 'T = ' + str(j) + ' K.png')\n",
    "\n",
    "    for i in filenames:\n",
    "        frames.append(Image.open(i))\n",
    "    gifname = output_dir + 'Movie.gif'\n",
    "    duration = int(1000/framerate)\n",
    "    frames[0].save(gifname, format='GIF',append_images=frames[1:],save_all=True,duration=duration, loop=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Example of an ice model.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADnCAYAAADl9EEgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO2deZgddZnvP+85vZymGwUJnThiiNxRId644DYqKFdJhOtCjDoqEFSQa7x6x1EY0DHOOOOMXhIQuQLGDRcIcSME11yujA8GcWRTgyIqTyC4JRgQJEkvZ3nvH1WVrj5dVae6q04tp9/P85yn+1RVV337PfWt+p5fVf1+oqoYhlE8KnkLMAwjGDOnYRQUM6dhFBQzp2EUFDOnYRSUnjCniLxdRBoicmfeWuIgIg+LyKSInJi3ljiIyC4R2ZW3jjiIyD63ti/OW0tS+vIWkBKHAAIcmreQmIwA+3F0l4GFeQuYBQcBj1Ce2obSE2dOw+hFzJyGUVDMnIZRUMychlFQzJyGUVDMnIZRUMychlFQzJxGJCJycN4a5itmTqMTr8xbwHzFzGl04q3eLyIynKeQ+YaZ0wjFjbT+e1RfkZeW+YiZ04jilcCk7/2ZeQmZj5g5jSjeCvgbhE6waJsdZk4jkIBIC85Z1KJtRpg5jTDaIy04Z1GLthlh5jTCaI+0HieIyEjWYuYjZk5jBiGR1sOibUaYOY0ggiKtx8H4rn0a3cPMaQQRFmk9LNpmQGn7EBKRVcBh7tvnuz9HRORs32LXqOpD2SqbiYgsAZb7JnkHxRNFxOv36F5V/V6WuoJwTfffOi2GE22/0n1FHYSI/BfgpQGzXi4iXt9H96jq9zOUlQpSxrFS3GttfwEmAO8fOMj9ud/92Q+8S1U/nbG8GYjIx4B3AXV3UhUYxNE+hmPWfaq6IB+FU4hIFbiYqQPfMvelwCbfov+sqvdkLG8GInIZcDZTtfX2A6+2VeDPqvr4HOQlopRnTlXdJyK3AH8TMNv7cCaAb2enKpJrgbcxMyoKjt4W03f83FDVJvB33nsR+SCOOXeq6mm5CQtnM3A68Ji26V5tm8CWrEWlQZm/c14B7IuY/2tV/X1WYjrwQ5ydJIx9wJUZaek1bsQxYhj7gasy0pIqZTbntYSf+cdwzFsIVLUFfJVwgyrwg+wU9Q6q2sDZF8K+nzWAH2WnKD1Ka05V3QP8JGR2BfhahnLicBVT34f9tHAarqLOrEY0XwIeDZjeBL7iHhxLR2nN6RIWbYsUaT3Coq1F2uSERdvSRloovzmDom2hIq1HRLS1SJuQiGhb2kgLJTdnSLQtYqT1CIq2FmnTISjaljbSQsnN6fI5pkfbIkZaj6Bo+6U8hPQgQdG21F8XesGcW5gebQsXaT1Cou22nOT0FL5o6+c/89CSFqU3Z0C0LWqk9ZgWbS3Spsq0FFLmSAs9YE6Xz3m/FDjSenS6IcGYOzfmLSBNesWcpbk9yxdtjZRxo23P0BPmdKNtmfhw3gKM4tMT5iwbqvq7vDUYxcfMaRgFxcxpGAUlc3OKyBoR+a2IfFREni4iUY/7xF3n0koVBmpC/2BlnYgsTUnroSJSF5E9InKJ78n6RFr7ByvrBmpCpeq8T0Oru+5rRKQhIttEZGUatQUOd2pbOTyN2orIU0TkSyLyqIg8koI+/7oP7AddqO35IrJTRP41zfVGoqqZvoBLce6BnAT2Ar8DPgo8Hbdnhlmsq1YbrmweGqnuX376qK5eu1hXrB6dHBqp7K8NVzYDtYRal7ha/a89wMeBhXPVuuKM0frqtYv1xNMO17S0utu4u01rA+cmh5Vzr22lnrS2wFOYur1uWj1T2qdm7AddqO1GV/MEzh1pO4F/BZam7RHvlXk3JSJyKfDOtsl1HLM+jHPL1SbgTu0gbmikuvnJx46ctOaCo4b6B6dCQH2ixYbzd4z95o69W8f2Nlcl0LoEuDdikQdxbir4qKruzlOrq/du4Kkhs5s4N4GvV9VvdFpXUr0i8hRgLfAaILQzMFVNfHbPqLYbgVPbJk/iHAD3AF8EvqyqdyXZzrRtRu3/InIyTl88afJ24L9HzG836tXAz9uNKiJLh0aqt114/bJpH8iBlUy0OGf59sb4vtZncUw0Fw5h5oEkjD04R9cZRs1IK8Aapvr+iSLSqHH0nrti+9jY3tazVfWXvr+LZcg2/j3mcmEsqI1Uzrro+qf3RWidHNvbeg9OSpsr5xDely90waidzDmGcxpPE2Fmfy+zplKFl71plNe/54jQZb528e+4YdMDtJqhT8nHIfGRvUxaYdZ6k5IoulWqSEytE60m40k2RXR3oXMmLD106uBrAKilLyc5ff3CoiXR0hYdWaOvX5hsJo9OSSiTVpi13qSbS/T/zkLr4GRTB5NsK2siW2tVtaqqkuYLuCwN4Y26suu+6APh7vvH660W6xNofVJZtLp6f5WV3j/eN0ajnry9Iun+1Gqxfvf945G37aVU26sT/7OzpAjXORs4rbZ/wmkFfR5Q6fihNHnaTVseHKtPBD94UJ9ose3aPY3GpH6+S7ofAT4NLCmB1hZwG04XkqnU9ofXPTjWarK0ww59LPB1nN4pukJjUr+wbfOD9RxqW8dptd0FrMf5XzvWNqROgeRlzhZThrwUp8fuhar6HlW9VbVzE7Kq3qWqWzecv2PGTlSfaLHhvB1jqmxVX4NFCvgNeYiqvl1VdxZUq9+Qfar6XFXdmGVtVfUnqvp6VT2ILhk149q2mDLkx4Hjgb9S1fPc/zXVSx95XEp5EfBanMsltyX5h0SkVhuuXA2ccvxrFlQWHVlj9/3j9W3X7mmosnV8X+tUVZ1zI4CIDAI3AXcAH4ljxE5aRTjpuJULhhYdWWPXzvHWTVv2TKSh1d3G3wGrcXacq4taWxF5FvCPwGNU9eVz1diu1V/bP943xg+ve3Asxdq+DHg5zn7707SNGLjNrM3ZDUTkxkqVF/f1i7ZaXNiY1M+nfBZKDRE5plLlrr5+odng1mZD31xUrQAicnOlygvKVttGXXFjdyG1xqFXzPkd4GSgqaqFH2JCRLyiX6Sq5+YqpgMi8j3gZUBDVdO+5p06vtqmcoNDnhShQcgwjADMnIZRUMychlFQzJxGzyAiT8xbQ5qYOY1e4gN5C0gTM6fRS7zB/0ZE4jyhU1jMnEZP4EbaQ9omvz8PLWlh5jR6haBIe0bmKlLEzGn0Cm8ImHa4iByeuZKUMHMapSck0nq8L0staWLmNHqBqFba1ZmpSBkzp9ELBEVaj9JGWzOnUWo6RFqPUkZbM6dRduLceFDKaFtKc4rIYW5P7Oo+InSyO6vqTXNf/5KnTg8Ruc6vyzfrHN/0vbkJ9CEij3d7jfe0vsyd1ddW23/MU6eP58RY5tCuq+gChX/2MYSHgWqM5R7ttpCYxNFRlAPlHuLVNtWhFBLwEuAo3/vt7k8FnuH+XrYhIoESP2wtItuBZR0Wq6lq2v3uzhoRORro9ET+V1U1qmEjMzr0HA/Ojt+vqoUbods92ytOtzJr89aThKIcrefCJzrMv6cIxgRQ1btxOoaK4t+y0BKTyzvMv7uIxuw1ymzOK4juLfxTWQmJyZaIeWOqemdmSjrzyQ7zU+l72IimtOZ0j9w/j1ik05k1a6LOjN/MTEUMVLVOeAfVCmzIUM68pbTmdAkzYGEirUeHaFukSOsRFm0t0mZE2c0ZFm2LFmk9rg2Ytr9gkdYjLNoWLZH0LKU2Z0S0LeoOFDTc3bcyVxGDkGirOD3eGxlQanO6tBuxcJHWIyTaFjHSerRHW4u0GdIL5myPtkWNtB7Tom1BI61He7QtaiLpSUpvzoBoW/QdKOlIzpkREG0t0mZI6c3pcsCQRY20HjFvSCgSB6KtRdps6RVzXpG3gFnSrXE4u0GnGxKMLtET5izbEV1V/1feGuLiRlsjB3rCnIbRi5g5DaOgZGpOEVkgIm8WkQUpr3dppQoDNaF/sLJORJamsM6qiLxDRP4pDY2+9S7tH6ysG6gJlarzPq31ishnReSozkvPbr1p19Zdb1VELhGR14pIKuNo+rRKtY+Xpqh1oYi8RUQ6dYeSLqra1RewADgb+E9gApgE3pvSumu14crmoZHq/uWnj+rqtYt1xerRyaGRyv7acGUzzvOcs1lfFXgH8AughXP9VNPWuuKM0frqtYv1xNMO17lqdde5FLgK2OtpxXkutHC1DVj/Mp/mBnAz8FrcZ4yTal1++mgzYW0XuvvC7cC4u9+enUZt47668rC1e2Z8DXAW8CygDgy7syeB96vqx5JuZ2ikuvnJx46ctOaCo4b6B6dCQH2ixYbzd4z95o69W8f2Nld10FoF/gfwLuAYYMZRXFMYITkNra7eo4G1wEqmaurna6r6t0XRG4aILGOq1wI/TeAW4CJgs8bYQVOs7UJgFfA24Gk4Bw2vxvuBv1fVz3RaT1qkFmvdyHq2iPwI+D1wMfB8YIDgnSjp9paKyIwPBKB/sMKaC44aEuEkETkm4G+rIrJGRH6Bc+C4HOcs1JVhypNodf/+aBG5yu1n6JfAaXShpmnpTUgVeAHwdaAuIjeLyKqw6JtCbRe6X19uA3YC64FjgUG6WOM4JOpDSEQOBt4InInzD/nPkAMhf1YB1opIossJlSqHHnfKYTM+EI/+wQrHr1rQ/x9f/tNbgfNcvW92fw88Q4YhIokuJ1SqVI5beVhlllqfAFxA+BkyjNdlofe4lQuGbtj0wM9FpDXHzcSpv2fUa4CmiNwCfEhVr/cW6BuQtxy/6rD+Wdb2scCpOMnuvzL9DDkYoeWjXejY7CFVfXbQjEhzisgunOwdlzBDtm/zUBL2iNbXLyx6Ui1ymYVPrPVVKvyDiPxDkm2R8CDW1y8sWpKZViEDvYuOrNHXL5XJpmbVqOgZ9f/6T6IDNWHh4sS1DTNk+zKDQNrDCoYO+JtX73uTJHwao9nQFbt2jr+QiGj+x/vGaNRT+U59Q5I/btT1r3fdN76YiLNFilohA727do5ro673A/fMcTMHA8+b498eoFFXdt03HrlMirWdIP17ox8Om5Fqg5Db7b3XEPQMpscFj1QahERk6dBI9bYLr18WGG3rEy3OXbF9bGxv69mqOqPnuzgNQR5JG4SSavXWgdOB8krgoIjNJW4QSkNvjG2ENQh5tIAf06FhKKXaLsJpKT4TpyGoycwal7dBCEBV/6Sqn1bV5+Ocrt+L0/I2ydTN3qk0uqjqXaq6dcP5O8bqE9O/9tQnWmw4b8eYKlvDPhBVbarqJ1X1aUA/8D9xGltSb75OqtW3jtNUdRhnB7oaZ4dJnTT0zpEWziW31wF9qvpCVb0mqsU2pdruUtXL3O9+S3C+m96Bc6b0atyVxsJIsrheAxyOc5b6Mc4H8NaU1lurDVc214YrzeWnj+rqDyS/FsfUtc67XK3NNLUOjVT2e1qTXotz17sU2Ihz8FPgsm7pTfk651Gu3ibwI+Z4jbPLtV0EvBPnWmcLeF0atY29/Sw35v7Dj+3COm+sVNGBmrT6BmQdcExK660Cj09Z6zGuVq32yS1paXXX/cQu1PaA3jRr6+mdqyE715abU9aa+n7b6VXaHt/9iMh3cMZLaapq4YeY8I2XcpGqnpurmBh4ejWFmzG6jfX4bhhG1zFzGkZBMXMakbh3Khk5YOY0OlGUcTjnHWZOoxNv9H5xn9owMsLMaYTiRtrH+SbZWTRDzJxGFO1mPC0XFfMUM6cRxRvb3h9m0TY7zJxGIAGR1sOibUaYOY0wwkxo0TYjzJxGGO2R1uMwEXl8pkrmKWZOYwYRkdbDom0GmDmNIDqZ702ZqJjnmDmNIMIirYdF2wwo/ONVYYjIJYC3gzzL/VkVka/6FvsndYbcyxUReSmwJmDWq0Rksfv7nar64QxlBeKaLirSerwPeHeX5XRERFbg9DPbzmtF5Cnu7z9V1Y9kKCsVSvk8p9tX0QMxFr1SVc/otp5OiMjtOF2HRqGqmfVkF4rbt9IPAW/ogScy1Z+OfyDds1V1W5baghCRO3G6t4yiparVLPSkSe47w1xQ1T8Bj8ZY9OJua4lJnE6hft11FTFQp2+lv1HVo1X1aJwuRAAa3jT3lbsxXeKMzfqLrqvoAqU0p8s1HeaPqepPMlHSmc/QueOwyzvMN4K5NMYyn+i8SPEoszk79R/6nUxUxECdwX3v6rCYjSA9B9QZ3Deq71ylfCOfAyU2p6reQ3S0Tbvz36REHeF/pTaCdBKiUsfPtWQjn3uU1pwuYdG2SJHWIyraWqRNRtSBr5SRFspvzrCzY2EirUeHaGuRNgER0ba0kRZKbs6IaFu0SOsRdIS3SJsOQemjtJEWSm5Ol6+3vS9ipPUIirZxWhuNzgTV8ZLMVaRIL5izfbSywkVaj5Bo+6k8tPQaIdG2tJEWesCcqrqD6dG2qJHWY9oR3iJtqkyLtlrG2998lN6cLgeibYEjrUecGxKMudFTXxF6xZyJBuLNkpg3JBhzoNdSSE+Y0422ZWJ13gKM4tMT5iwbJYjeRgEwcxpGQTFzGkZBydycIvIuEdklIh8TkeeISOIBWUVkaaUKAzWhf7CyTkSWpqT1MBFpiMjDIvIpEXliGlr7ByvrBmpCpeq8T0Oru+5viEhTRG4VkVPTqC0w7Na2mnJtn+DW9iERuTyN0cz8+0EXavsBEfmDiFwgIs9MqbbRZD2UNk5ztwJ1YC+wG/gY8BxmOQQ5UKsNVzYPjVT3Lz99VFevXawrVo9ODo1U9teGK5uBWkKtS1yt/tfDwAZmOcS7X+uKM0brq9cu1hNPO1zT0upu4+42rU3gVuDUuda2Nlxpdqm2ywJq+xBwGfCEpPtBF2q70dU46e63fwAuAJ4529rGfWXeTYmIXAq8s21yA5gA9uEU4Wrgdu0gbmikuvnJx46ctOaCo4b6B6dCQH2ixYbzd4z95o69W8f2Nlcl0LoEuDdikUeALwP/rqq/zVOrq/du4Kkhs1vAHTi9Q2wqQG2XAdsjFvkzsAln+Pjf56nV1bsR5yDnp45j1r8AV7p6f9aptrG3GbUeETkF6E9jQz7eBrw8Yn67UTcCd7T/wyKydGiketuF1y+b9oF41CdanLN8e2N8X+tKnCPyXDgEOCvmsg8DXyHAqBlpBTgTODTGcp5RLwK+klNtR4l/SSnUqHG0nrt8++TYvtZ5QKTJO/D3wIsi5rcb9WpgexKjdjLnuLvBNBFgJOlKKlV42ZtGef17jghd5msX/44bNj1AK+fnEsqkFcqldxZaJ1tNJhJsKpX9NghVDfz+2qlrTGGq57W0SKUXtL5+YdGSWuQyi46s0dcvTDbzvVuuTFqhXHpnobV/spmoB77MG08jN6iqg6ral+YL5wt/GHWcOPsH4H8DzwAqqirtr1aL9bvvH29E6d99/3i91WJ90N/HeQFPilnHB4H/AyzKS6ur91dR23Bp4nR9+Rryre3T45WWPwEfB0YTar0w4X67KWITk8B+4Lc4t5I+bZa1CKQI1zn9hrwIeCFwhKq+X1VDM3tjUr+wbfOD9fpEK3ilEy22Xbun0ZjUz3dJt9+QC1T13aq6u6Ba/YbsV9XjVHVLgWvrN+Soqr5Hne5Qi6TVb8h1wHNVdbGqflBVU7l3Oi9ztpgy5MeIaUg/qnqXqm7dcP6OsfYPpj7RYsN5O8ZU2aqqv0xRd2xDFkBrE7iZmIYsgN7YhsxRa5MpQ64Hnpe2If3kcSnl+cBKnJhwZ5ydJWJdtdpw5WrglONfs6Cy6Mgau+8fr2+7dk9Dla3j+1qnqup4gvUPAjcAt+O0FHY0YietIpx03MoFQ4uOrLFr53jrpi17JtLQ6m7jHTjjZ14IXJdGbf16U67tEPA94Bac2nY04my0dqG2JwAn4lyGyqST6lIOx9COiNxYqfLivn7RVosLG5P6+ZSP6qkhIsdUqtzV1y80G9zabOibi6oVputttVhfntrqj5oNziqq1jiUdiCjNva1mjDZ1Jaqnpe3mChU9ZciB1o5f1D0ncevt0S1VeA/il7bThShQcgwjADMnEYkmdzgbQRi5jQ68Ya8BcxXzJxGJ97r/WJn0WwxcxqhuGZ8tm/S6/LSMh8xcxpRvIHp+8h7wxY00sfMaUTRbsbnWbTNDjOnEUhApAVnf7FomxFmTiOM9kjrYdE2I8ycRhjnhEy3aJsRZk5jBq75jg2ZbdE2I8ycRhBhkdbDom0GmDmNIMIirYdF2wwwcxrT6BBpPSzaZkApzenriV1FRIGT3VlVb5r7KsRAuiLybb8u36xzfNP35yZwOkcQb794YbeFxEFEbgiorQAf8E3fm6fGuVLW5zkfIt4OtKfbQmLyYIxlCvHUu6r+VkT+Chh0J21kyoj+Ds92ZiosnDi1De5gqOCU8szpdr/xsxiLFmWk4ziD+27puoqYqOofVfU+Vb0PGHMnN7xp7qsQBxPgwzGW+XrnRYpHKc3p8okO83+jBRnpWFV/jTO+RhSFiOBlQ1XvZOoAEkZpRj73U2Zzfp7oKPjJrITE5NqIefu70XvbPOJbEfMe1fKNfA6U2Jwxom1RIq1H1NG7MJG2pERF21JGWiixOV3Com1hIq1Hh2hrkTYBHaJtKSMtlN+cYdG2aJHW45qAaRZp0yEo2pY20kLJzRkRbYsWaT0+EjDNIm06BEXb0kZaKLk5XdqjbeEirUdItLVImwIh0ba0kRZ6w5zt0baokdZjWrS1SJsq06JtmSMt9IA5A6JtUSOtR1C0NdIhzg0JpaH05nQ5EG2LGmk9Yt6QYMwBN9r2DL1izm6NE9ktPp23AKP49IQ5C3SfZyxUtdPzkobRG+Y0jF7EzGkYBSVzc4rIUSLyFhF5XIrrXFqpwkBN6B+srBORpSmttyoil7t6U+mWQ0SW9g9W1g3UhErVeZ/Get11nyginxaRJWmtExh2a1vtQm0/JSKr06yttx90obZHu1oPSWudHVHVTF84lzomgXHgJuBM4HFzXFetNlzZPDRS3b/89FFdvXaxrlg9Ojk0UtlfG65sBmoJtS7BuYaqOA/s/gR4C+6I4HPVuuKM0frqtYv1xNMO17S0utu426f3EeCzwJIkta0NV5pdqu0yn9YmcDuwOmltPa1dqO1GYMJ9fR84Azgk6XqjXpkPOy8ilwLv9E3aC/QDtwFXAFtU9aE46xoaqW5+8rEjJ6254Kih/sGpEFCfaLHh/B1jv7lj79axvc1VCbQuAe4NmOVdW70E+KLGKGK3tbp67waeGjDrL8DXgH9T5wHqjmRQ22XA9oBZLeCnwMeBqwpU243Aqb5Je4EB4GacqwXfUNWHk2yjnSJ85xzB6RLjRTg7+x9E5CYROTMq+orIUhGZ8YEA9A9WWHPBUUMinCQix3RBswDPxPlQmiLyk6jom7NWgMcAZwH3isgjIvLZqOibs94KTgdjXwIaInJ7VPTNUesIjjlPAC4DdovI90XkjLSib2QfQiKywxWQJlHCR9yfLwKeAVwuIt4Z9WpVHfcW7BuQtxy/6rD+9g/Eo3+wwnErFwzdsOmBX4jIXPuQifNdyG/UK0TkZ8CHVXVzxloBqjGW8Yx6loh4Z9R3q+q+jPXGqa3fqF9wa/shVf3GrLSesmDohi8/8GP3/50rUW0k3n57AvAc4DMi4p1RN+kcb4zp1MHXEcT7wLuB36gvAj7nP3gO1ISFi2uRK1h0ZI2+fpHJpmb1P3hGvaYEWmG6UQ9MLKjeCvAs4LpZa11So69fRiabOtxdicB0o54AfLFTe5eqBi7QyZzvw/k+mCanAC9IupJGXdl133jkMrt2jmujrvcCv5rjZg4CXjLHvz1ARlrB0XpQgr8HMtN7MHDcHP/2ADG1tpoN3QZ8N8Gm3ohz4M2MIjQIeTyKE6G/jxMHvu2PWgHrWTo0Ur3twuuXzfiuAU5jwLkrto+N7W09W1V/OUetSwhuEPJo4jQIrFfVb+ap1d1OWIOQxwPAlcBHVTW0S8mMahvWIOTRAH4IXKCqoabKsLbtDUIej+KcwP4f8AXgu6raqcOxWOTdIPQoTtP0VuBtwGGqerKqfjXKmOA8aqWqWzecv2OsPjH9a099osWG83aMqbI1yQcSQhPnEtApqtqnqi+OMmbOWsEx5EXAAlVdqKrnRhkzZ70N4AfAK1W1X1VPiDJmjlofxbkU+E2cS2uPU9VXq+rmtIwJ5HKd859xrnN+F/hbYDjBupzrW8OV+vLTR3X1B1K/FrcQ57JJA9gGvDqx1pHKfk/r8tWjjZSvxd3s6t0NXIhzsEtNb5euIdeBG4FXFLy2F+OcSL4BrAKGkq6z4za7vYGAf1KAg1Je54WVKq2BmuzrG5B1wDEprvvwlLUe42qdrPaxJWWt1SSGjNCrAzXRktRWndrK5pS1VrIwpP+V+XfObiAi5+M8xPxHVT0ibz2dEJEGsB84U1UL38+NNwaJhrQqFglX6yPAGeq75FJG8v7OaRhGCGZOwygoZk4jEhF5bN4a5itmTqMTr85bwHzFzGl04kzvFxF5TJ5C5htmTiMUN9L6b7V8ZV5a5iNmTiOKV+PcJOBxZtiCRvqYOY0ozmTqKQuA4yzaZoeZ0wgkINKCc9ulRduMMHMaYbRHWnAe87JomxFmTiOM9kjrYdE2I8ycxgzcSPvCkNmTwKsylDNvMXMaQbwax4RBHAy8NUMt8xYzpxFEWKT1sGibAZ36ECosInIGsMB9+2L358Ei8l7fYhtVdXe2ymYiIk8FXuGb5B0UXyUii93f7ynCI06u6Y7vsFgFp9X26u4risbt8vLkgFkrReSv3d9/rarfClim0JTyeU4RGcZ5Zq/F1KjWXheeXhwT4J2q+pmM5c1ARC4C3o3TxQk4D0V7vdZN4uzsj6jqgoA/zxQRqQAfYqoryOcCz8Op8+W+Rder6s5s1c3E7ZNqDVO19Xfl6tV2j6o+Pmtticnyye6Un3rfxlR3/kGvcWBR3jpdrS/A6XcmTGsD+FTeOkO0f9DVeG/eWkL0vQSnR/uo2l6St865vMr8nfMKokeIvktVd2UlpgM/JryBBZxeEa7KSEuvcRNOggpjjALE77lQZnNuIbxP3f3A5zLUEomqtoBNTEWvdpo43UAas0RVm8DXCTfoOHBLdorSo2P70YcAAAFVSURBVLTmVNU/A7eGzK4C12QoJw4bcY7i7TSBr7oGNubGlUBQV6pNnGE8ytewQonN6RIWbYsUaT3Coq1F2uSERdvSRloovzmDom2hIq1HRLS1SJuQiGhb2kgLJTdnSLQtYqT1CIq2FmnTISjaljbSQsnN6fI5pkfbIkZajx/j9Bru58o8hPQgQdF2Yx5C0qIXzHkd06Nt4SKth3uG/DLTo+3NOcnpKXzR1k9Yg2EpKL05A6JtUSOtx7Roa5E2VaalkDJHWugBc7ocOFsWONJ6BEVbIx1uyltAmvSKOa/LW0BcfNHWSBk32vYMPWFON9qWiQ/nLcAoPj1hzrKhBXiMzSg+Zk7DKChmTsMoKGZOwygoZk7DKChmTsMoKGZOwygoZk6jE3WcfniiulkxukApe98zjPmAnTkNo6CYOQ2joJg5DaOgmDkNo6CYOQ2joJg5DaOg/H/3YqU34p6oiAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "if __name__==\"__main__\":\n",
    "\t#EXAMPLE OF A STATE ARRAY\n",
    "\tprint(\"Example of an ice model.\")\n",
    "\tarr = [[[-1,1],[-1,1],[1,-1],[-1,1]],\n",
    "\t         [[-1,-1],[-1,-1],[-1,1],[1,-1]],\n",
    "\t         [[-1,1],[-1,1],[1,1],[-1,1]],\n",
    "\t         [[-1,1],[-1,1],[1,1],[-1,1]]]\n",
    "\tarr = np.array(arr)\n",
    "\n",
    "\tvisualise_2d_model(arr)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "plist,arr = equilibrate(3,temp=0.1,return_polarlist=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(1000):\n",
    "    arr = initialise_state(3)\n",
    "    arr2 = short_loop(arr)\n",
    "    if not check_config(arr2):\n",
    "        print(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "check_config(arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "(1, 1, 1, -1)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-27-190a90ff61ab>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mvisualise_final_state\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0marr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m<ipython-input-11-7b44644f2693>\u001b[0m in \u001b[0;36mvisualise_final_state\u001b[1;34m(arr, title, savefig, savename, show)\u001b[0m\n\u001b[0;32m    105\u001b[0m     '''\n\u001b[0;32m    106\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 107\u001b[1;33m     \u001b[0matom_config_arr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcalculate_atom_config\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0marr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    108\u001b[0m     \u001b[0mcolordict\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m{\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m5\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m6\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m}\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    109\u001b[0m     \u001b[0matom_config_arr_copy\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0matom_config_arr\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-11-7b44644f2693>\u001b[0m in \u001b[0;36mcalculate_atom_config\u001b[1;34m(arr)\u001b[0m\n\u001b[0;32m     26\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     27\u001b[0m             \u001b[0mtup\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mcurrent_up_state\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcurrent_right_state\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcurrent_down_state\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcurrent_left_state\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 28\u001b[1;33m             \u001b[0matomtype\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtupdict\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mtup\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     29\u001b[0m             \u001b[0mtypearr\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0matomtype\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     30\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mtypearr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mint\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyError\u001b[0m: (1, 1, 1, -1)"
     ]
    }
   ],
   "source": [
    "visualise_final_state(arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "tupdict = {\n",
    "    (1,1,-1,-1):1,\n",
    "    (-1,-1,1,1):2,\n",
    "    (1,-1,-1,1):3,\n",
    "    (-1,1,1,-1):4,\n",
    "    (-1,1,-1,1):5,\n",
    "    (1,-1,1,-1):6\n",
    "}\n",
    "num=len(arr)\n",
    "typearr = np.zeros((arr.shape[0],arr.shape[1]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "ename": "KeyError",
     "evalue": "(1, 1, 1, -1)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyError\u001b[0m                                  Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-32-601a9b65ebee>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     12\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     13\u001b[0m         \u001b[0mtup\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mcurrent_up_state\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcurrent_right_state\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcurrent_down_state\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcurrent_left_state\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 14\u001b[1;33m         \u001b[0matomtype\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mtupdict\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mtup\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     15\u001b[0m         \u001b[0mtypearr\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0matomtype\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyError\u001b[0m: (1, 1, 1, -1)"
     ]
    }
   ],
   "source": [
    "for i in range(len(arr)):\n",
    "    for j in range(len(arr)):\n",
    "        current_up_state = arr[i][j][0]\n",
    "        current_right_state = arr[i][j][1]\n",
    "\n",
    "        lower_neighbour_up_state = arr[(i+1)%num][j][0]\n",
    "\n",
    "        left_neighbour_right_state = arr[i][j-1][1]\n",
    "\n",
    "        current_down_state = -(lower_neighbour_up_state)\n",
    "        current_left_state = -(left_neighbour_right_state)\n",
    "\n",
    "        tup = (current_up_state,current_right_state,current_down_state,current_left_state)\n",
    "        atomtype = tupdict[tup]\n",
    "        typearr[i][j] = atomtype\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "check_config(arr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}