Plotting_workshop_2021.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "Plotting_workshop_2021.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/koysean/plotting_workshop/blob/main/Plotting_workshop_2021.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7yDof6fhi7q2"
      },
      "source": [
        "# Plotting Course 2021\n",
        "Authored by Jaron Ma and Sean Koyama\n",
        "Thi is a Google Colab Notebook, which is similar to a Jupyter Notebook in that it combines coding blocks with text/image annotations. It runs Python code on Google's server's and feeds the output to your browser."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "vxfr0bwujRPv"
      },
      "source": [
        "### Basic Plotting\n",
        "As a brief refresher, numpy and matplotlib are designed to imitate MATLAB's advantages in creating, manipulating, and plotting matrices of data, which are sometimes considered vectors/arrays/tables. If you've used these packages before, the script below should look pretty familiar:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1nFehn8Fjn1T"
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "\n",
        "x = np.linspace(-0.3,1.3,40)\n",
        "y = x**2 * (1-x)**2\n",
        "plt.plot(x,y)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-7Uk4YTHj0xi"
      },
      "source": [
        "The np.plot command expects as a bare minimum an input that can be interpreted as a list of data. If there are two inputs of equal size, the first is interpreted as the x-coordinates of the data, and the second is interpreted as the y-cooridnates of the data, with the datapoints connected by a line in the order that they appear."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "DU8knDLxkLzH"
      },
      "source": [
        "###Beyond the Default\n",
        "\n",
        "The default output from np.plot is not bad on its own: The blue color of the plot used is pleasant, and the axis limits are scaled to the data. This makes it easy to get in the habit of falling back on the defaults, but our first tip for plotting is to be aware of how easy it is not to use the default settings. An equivalent command as above would have been:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "LjujWvcZkaWR"
      },
      "source": [
        "plt.plot(x,y,linestyle='-',color='#1f77b4',linewidth=1.5,marker='')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ElgMRlwrlaQP"
      },
      "source": [
        "Here, we are using **name-value pairs**. The properies named linestyle, color, etc., can be listed followed by the values you want them to take on. The values given above correspond to the default values chosen by MATLAB (as of R2021b), which makes it clearer how to change the plot's apperance. For example, to mix things up:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "44wkgw2z5QGI"
      },
      "source": [
        "plt.plot(x,y,linestyle=':',color=[0.3, 0.7470 ,0.0410], linewidth=3, marker='D',markersize=8)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mYlJ9tuz5rj7"
      },
      "source": [
        "###Getting and Setting\n",
        "The next tip we have will further help customize your plots: consider taking the time to retain your object handles so that you can **get and set** properties in the code itself.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "7Uh2B4IB537n"
      },
      "source": [
        "fig_1 = plt.figure(); \n",
        "ax1_1 = fig_1.add_subplot()\n",
        "plot1, = ax1_1.plot(x,y)\n",
        "plot1.set(linestyle=':',color=[0.3, 0.7470 ,0.0410], linewidth=3, marker='D',markersize=8)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2SSAzUI8C5lr"
      },
      "source": [
        "The code above does a few things: It initializes a figure object and stores the handle as fig_1, it places an axes object inside the figure handle identified by fig_1, and it places a line plot on the axes handle identified by ax_1. By storing the object handles as variables, you can go back and use the set command to change the properties after the fact. You can use the get command to print out a list of all of the available properties for a given object:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Aw5D843nC8NW"
      },
      "source": [
        "print(plot1.pickable())"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "K3QRA1ZlEZVq"
      },
      "source": [
        "That leads us to some basic tips:\n",
        "* Make your fonts bigger\n",
        "* Consider choosing colors with a tool like ColorBrewer (hexcodes are useful here)\n",
        "* Consider rotating your y-axis labels\n",
        "* Export as a vector file when you can (e.g., .pdf and .svg)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "UO0yQDOUEf0q"
      },
      "source": [
        "plt.rcParams.update({'font.size': 14})\n",
        "plot1.set_color('#8da0cb')\n",
        "ax1_1.set_ylabel('Y-axis', rotation=0)\n",
        "\n",
        "fig_1"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ogssaazqF9j1"
      },
      "source": [
        "### Subplots"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "f6laHcFSGItz"
      },
      "source": [
        "As a reminder, you can make subplots in matplotlib using the subplots command. The usual format is subplots(m,n), where m and n specify the dimesions (rows, columns) of the subplot grid."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "lcINfVcBGHXL"
      },
      "source": [
        "fig_2, (ax_2_1, ax_2_2) = plt.subplots(1, 2)\n",
        "\n",
        "line_2_2 = ax_2_1.plot(x,(x-0.5)**2)\n",
        "line_2_3 = ax_2_2.plot(x,(x-0.5)**3)\n",
        "\n",
        "ax_2_1.set_title('parabola')\n",
        "ax_2_2.set_title('whatever this shape is')\n",
        "fig_2.suptitle('In MATLAB this would be sgtitle')\n",
        "\n",
        "fig_2.subplots_adjust(top=0.77)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hmk8OI0G_SIh"
      },
      "source": [
        "## Example: Group 14 Solids\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "6S1cwVxbpfei"
      },
      "source": [
        "import matplotlib as mpl\n",
        "import matplotlib.cm as cm\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import numpy.polynomial.polynomial as poly\n",
        "from scipy.optimize import curve_fit"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "3oZE6N1rOHRB"
      },
      "source": [
        "First import data from github, then parse the CSV files into numpy arrays. With this code the data will be arranged into 2 x *n* arrays, with the first row containing volumes (Å<sup>3</sup>) and the second row containing energy (eV).\n",
        "\n",
        "*Note: raw energy values should not be directly compared as they do not share a well-defined reference point*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "JCzdCu_x3iQK"
      },
      "source": [
        "!wget https://raw.githubusercontent.com/koysean/plotting_workshop/main/c_e_v.csv\n",
        "!wget https://raw.githubusercontent.com/koysean/plotting_workshop/main/si_e_v.csv\n",
        "!wget https://raw.githubusercontent.com/koysean/plotting_workshop/main/ge_e_v.csv\n",
        "\n",
        "\n",
        "c_e_v = [[], []]\n",
        "si_e_v = [[], []]\n",
        "ge_e_v = [[], []]\n",
        "\n",
        "# Apparently numpy has a builtin function, np.readfile(),\n",
        "# that does all of this automatically :)\n",
        "with open('c_e_v.csv', 'r') as file:\n",
        "  for line in file.readlines():\n",
        "    split_line = line.split(',')\n",
        "    c_e_v[0].append(float(split_line[0]))\n",
        "    c_e_v[1].append(float(split_line[1]))\n",
        "  c_e_v = np.asarray(c_e_v)\n",
        "\n",
        "with open('si_e_v.csv', 'r') as file:\n",
        "  for line in file.readlines():\n",
        "    split_line = line.split(',')\n",
        "    si_e_v[0].append(float(split_line[0]))\n",
        "    si_e_v[1].append(float(split_line[1]))\n",
        "  si_e_v = np.asarray(si_e_v)\n",
        "\n",
        "with open('ge_e_v.csv', 'r') as file:\n",
        "  for line in file.readlines():\n",
        "    split_line = line.split(',')\n",
        "    ge_e_v[0].append(float(split_line[0]))\n",
        "    ge_e_v[1].append(float(split_line[1]))\n",
        "  \n",
        "  ge_e_v = np.asarray(ge_e_v)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TOkQZwabPT4R"
      },
      "source": [
        "We can find the equilibrium volume and energy by looking at where the energy is at a minimum."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "kmtf4yELU2AK"
      },
      "source": [
        "min_volumes = {\n",
        "    'c':  c_e_v[0][np.argmin(c_e_v[1])],\n",
        "    'si':  si_e_v[0][np.argmin(si_e_v[1])],\n",
        "    'ge':  ge_e_v[0][np.argmin(ge_e_v[1])]\n",
        "}\n",
        "\n",
        "min_energies = {\n",
        "    'c':  np.min(c_e_v[1]),\n",
        "    'si':  np.min(si_e_v[1]),\n",
        "    'ge':  np.min(ge_e_v[1])\n",
        "}\n",
        "\n",
        "print('(Roughly) Equilibrium values of (V, E):')\n",
        "print('C:\\t{:.3f}, {:.3f}'.format(min_volumes['c'], min_energies['c']))\n",
        "print('Si:\\t{:.3f}, {:.3f}'.format(min_volumes['si'], min_energies['si']))\n",
        "print('Ge:\\t{:.3f}, {:.3f}'.format(min_volumes['ge'], min_energies['ge']))"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U4n89z2vQHyU"
      },
      "source": [
        "With this data, we can plot the energy-volume curve!"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OU-COrxkUdR4"
      },
      "source": [
        "fig, ax = plt.subplots()\n",
        "\n",
        "sc_c = ax.scatter(c_e_v[0], c_e_v[1])\n",
        "sc_si = ax.scatter(si_e_v[0], si_e_v[1])\n",
        "sc_ge = ax.scatter(ge_e_v[0], ge_e_v[1])\n",
        "\n",
        "ax.legend([sc_c, sc_si, sc_ge], ['C', 'Si', 'Ge'])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "G-u6BV8g7t4l"
      },
      "source": [
        "We can try a rudimentary second-order polynomial fit to the data, excluding the points where the curves are obviously non-parabolic."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Yan8K46t7t4r"
      },
      "source": [
        "poly_c = poly.polyfit(c_e_v[0][4:16], c_e_v[1][4:16], 2)\n",
        "poly_si = poly.polyfit(si_e_v[0][4:16], si_e_v[1][4:16], 2)\n",
        "poly_ge = poly.polyfit(ge_e_v[0][4:16], ge_e_v[1][4:16], 2)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "t3DzQ_lD7t4x"
      },
      "source": [
        "We can plot these fits along with the real data."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "PnKZa-Au7t41"
      },
      "source": [
        "fig, ax = plt.subplots()\n",
        "\n",
        "sc_c = ax.scatter(c_e_v[0], c_e_v[1])\n",
        "ax.plot(c_e_v[0], poly.polyval(c_e_v[0], poly_c))\n",
        "\n",
        "sc_si = ax.scatter(si_e_v[0], si_e_v[1])\n",
        "ax.plot(si_e_v[0], poly.polyval(si_e_v[0], poly_si))\n",
        "\n",
        "sc_ge = ax.scatter(ge_e_v[0], ge_e_v[1])\n",
        "ax.plot(ge_e_v[0], poly.polyval(ge_e_v[0], poly_ge))\n",
        "\n",
        "ax.legend([sc_c, sc_si, sc_ge], ['C', 'Si', 'Ge'])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2vzHk_R4QNGG"
      },
      "source": [
        "Not too good. We can do better and also extract physical parameters from this data by fitting to an equation of state. In this example, we fit to the Birch–Murnaghan equation of state:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gcSUIedUdC7w"
      },
      "source": [
        "$E(V) = E_0 + \\dfrac{9V_0B_0}{16} \\left\\{ \\left[   \\left(  \\dfrac{V_0}{V}\\right)^{\\frac{2}{3}} - 1 \\right]^3 B_0'  + \\left[  \\left(  \\dfrac{V_0}{V}\\right)^{\\frac{2}{3}}-1 \\right]^2 \\left[ 6 - 4 \\left(  \\dfrac{V_0}{V}\\right)^{\\frac{2}{3}}    \\right]    \\right\\}$"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "qKFQmENh9Vqy"
      },
      "source": [
        "def bm_eos(v, v0, b0, bp0, e0):\n",
        "  ''' Birch-Murnaghan EOS. \n",
        "  Arguments:\n",
        "    v: volume input\n",
        "    v0: equilibrium volume\n",
        "    b0: bulk modulus at equilibrium\n",
        "    bp0: derivative of bulk modulus wrt pressure at equilibrium\n",
        "    e0: equilibrium energy\n",
        "  \n",
        "  Outputs:\n",
        "     Internal energy of the system, E(V)\n",
        "  '''\n",
        "\n",
        "  gamma = (v0/v)**(2/3)\n",
        "\n",
        "  return e0 + 9/16 * v0 * b0 * (\n",
        "      (gamma - 1)**3 * bp0\n",
        "      + (gamma - 1)**2 * (6 - 4*gamma)\n",
        "  )"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MijyUlMPhrvZ"
      },
      "source": [
        "Using scipy's optimize library, we can fit our data to the B-M equation of state. The fit coefficients are given in the order of the parameters as defined in the function above, i.e. (v0, b0, bp0, e0):"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Ry3ThQF9_OQW"
      },
      "source": [
        "fit_c, conv_c = curve_fit(bm_eos, c_e_v[0], c_e_v[1], p0 = [min_volumes['c'], 2, 2, min_energies['c']])\n",
        "\n",
        "fit_si, conv_si = curve_fit(bm_eos, si_e_v[0], si_e_v[1], p0 = [min_volumes['si'], 2, 2, min_energies['si']])\n",
        "\n",
        "fit_ge, conv_ge = curve_fit(bm_eos, ge_e_v[0], ge_e_v[1], p0 = [min_volumes['ge'], 2, 2, min_energies['ge']])\n",
        "\n",
        "print('Fit parameters for C:', fit_c)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "JHn0VoqXtoS6"
      },
      "source": [
        "Again, plotting against the real data:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "sIWMknywBw2_"
      },
      "source": [
        "fig, ax = plt.subplots()\n",
        "\n",
        "sc_c = ax.scatter(c_e_v[0], c_e_v[1])\n",
        "ax.plot(c_e_v[0], bm_eos(c_e_v[0], *fit_c))\n",
        "\n",
        "sc_si = ax.scatter(si_e_v[0], si_e_v[1])\n",
        "ax.plot(si_e_v[0], bm_eos(si_e_v[0], *fit_si))\n",
        "\n",
        "sc_ge = ax.scatter(ge_e_v[0], ge_e_v[1])\n",
        "ax.plot(ge_e_v[0], bm_eos(ge_e_v[0], *fit_ge))\n",
        "\n",
        "ax.legend([sc_c, sc_si, sc_ge], ['C', 'Si', 'Ge'])"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5bEDZDTxuvV3"
      },
      "source": [
        "And it returns a pretty good fit!"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ebom08YZux-E"
      },
      "source": [
        "So far we've just created the most basic plots we can with this data and the fits. See if you can make these plots any nicer! Think about what you can communicate with this data and then try to make a clear and `a e s t h e t i c` plot.\n",
        "\n",
        "Remember, the data you have to work with for each material (C, Si, Ge) is:\n",
        "\n",
        "*   Raw DFT energy-volume data (i.e. c_e_v)\n",
        "*   The Birch-Murnaghan equation of state and its fitted parameters, including bulk modulus (i.e. bm_eos(v, *fit_c), fit_c)\n",
        "*   The quadratic fitting coefficients (i.e. poly_c)\n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "\n",
        "\n",
        "Some questions to get you thinking about what or how to make a good plot from this data:\n",
        "\n",
        "*   What are the quantities of interest?\n",
        "*   What do the fits show that the raw data doesn't?\n",
        "*   What do you want your audience to learn from your plot?\n",
        "*   How can you minimize the amount of thinking required on behalf of your audience in interpreting your plot; in other words, how can you make your plot efficient?\n",
        "\n",
        "Once you have a visualization, ask yourself those questions again to see if you can find further ways to improve your plot. Here are a few more questions you can ask yourself:\n",
        "\n",
        "*   Is the feature or quantity you're trying to communicate the most obvious part of the plot?\n",
        "*   Are the colors visually distinct? What about to someone with colorblindness (red-green colorblindness is the most common type)?\n",
        "*   Does understanding your plot require frequent visual lookup (i.e. looking back and forth at the plot and a legend to make sense of the plot)?\n"
      ]
    }
  ]
}