examples/UsingCalculationResults.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploring `calculate` and `equilibrium` xarray Datasets\n",
    "\n",
    "## xarray Datasets\n",
    "\n",
    "Results returned from calling `calculate` or `equilibrium` in pycalphad are [xarray](http://xarray.pydata.org/en/stable/) Datasets. An xarray Dataset is a data structure that represents N-dimensional tabular data. It is an N-dimensional analog to the Pandas DataFrame.\n",
    "\n",
    "This notebook will walk through the structure of xarray Datasets in pycalphad and some basics of using them. For more in-depth tutorials and documentation on using xarray Datasets and DataArray's fully, see the [xarray documentation](http://xarray.pydata.org/en/stable/index.html).\n",
    "\n",
    "## Dataset structure\n",
    "\n",
    "Each Dataset stores the conditions that properties are calculated at and the values of the properties as a function of the different conditions. There are three key terms:\n",
    "\n",
    "* `Dimensions`: these are the conditions that are calculated over, e.g. pressure (P) and temperature (T). They are essentially labels.\n",
    "* `Coordinates`: these are the actual *values* that are taken on by the dimensions. \n",
    "* `Data variables`: these are the properties calculated by pycalphad, such as the Gibbs energy, mixing energy, composition, etc.\n",
    "\n",
    "## `calculate()` results\n",
    "\n",
    "Calculate is used to sample properties of a single phase. There are five dimensions/coordinates:\n",
    "\n",
    "* `P`: pressures (in Pa).\n",
    "* `T`: temperatures (in K).\n",
    "* `component`: the string names of the components in the system\n",
    "* `internal_dof`: The internal_dof (internal degrees of freedom) is the index of the site in any phase's site fraction array. Below the FCC_A1 phase has the sublattice model (AL, ZN) and thus the internal_dof are integers 0 and 1 referring to the AL site (index 0) and the ZN site (index 1).\n",
    "* `points`: By default, the calculate function samples points over all of the internal degrees of freedom. Each coordinate point simply represents the index is a list of all configurations of the internal_dof sampled. There is no underlying physical meaning or order.\n",
    "\n",
    "There are also at least four Data variables:\n",
    "\n",
    "* `Phase`: The string name of the phase. For `calculate`, this will always be the phase name passed.\n",
    "* `X`: The composition of each component in mole fraction as a function of the temperature, pressure, and the index of the points (there is one composition for each point).\n",
    "* `Y`: The site fraction of each index in the internal_dof array for the given temperature, pressure and point.\n",
    "* `output`: \"output\" is always whatever property is calculated by the output keyword passed to `calculate`. The default is the molar Gibbs energy, GM. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.Dataset>\n",
      "Dimensions:    (N: 1, P: 1, T: 2, component: 2, internal_dof: 2, points: 4001)\n",
      "Coordinates:\n",
      "  * component  (component) <U2 'AL' 'ZN'\n",
      "  * N          (N) float64 1.0\n",
      "  * P          (P) float64 1.013e+05\n",
      "  * T          (T) float64 500.0 1e+03\n",
      "Dimensions without coordinates: internal_dof, points\n",
      "Data variables:\n",
      "    X          (N, P, T, points, component) float64 1.0 1e-15 ... 0.7439 0.2561\n",
      "    Phase      (N, P, T, points) <U6 'FCC_A1' 'FCC_A1' ... 'FCC_A1' 'FCC_A1'\n",
      "    Y          (N, P, T, points, internal_dof) float64 1.0 1e-15 ... 0.2561\n",
      "    GM         (N, P, T, points) float64 -1.559e+04 -2.01e+04 ... -4.81e+04\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "from pycalphad import Database, calculate, equilibrium, variables as v\n",
    "\n",
    "dbf = Database('alzn_mey.tdb')\n",
    "comps = ['AL', 'ZN', 'VA']\n",
    "calc_result = calculate(dbf, comps, 'FCC_A1', P=101325, T=[500, 1000])\n",
    "print(calc_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can manipulate this by selecting data by value (of a coordinate) using `sel` or index (of a coordinate) using `isel` similar to a Pandas array. Below we get the site fraction of ZN (internal_dof index of 1 selected by index) at 1000K (selected by value) for the 50th point (selected by index). \n",
    "\n",
    "The results of selecting over Data variables gives an xarray DataArray which is useful for plotting or performing computations on (see [DataArrays vs Datasets](http://xarray.pydata.org/en/stable/data-structures.html))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.DataArray 'Y' (N: 1, P: 1)>\n",
      "array([[0.97648824]])\n",
      "Coordinates:\n",
      "  * N        (N) float64 1.0\n",
      "  * P        (P) float64 1.013e+05\n",
      "    T        float64 1e+03\n"
     ]
    }
   ],
   "source": [
    "print(calc_result.Y.isel(internal_dof=1, points=49).sel(T=1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "accessing the `values` attribute on any on any DataArray returns the multidimensional NumPy array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[[[1.00000000e+00 1.00000000e-15]\n",
      "    [1.00000000e-15 1.00000000e+00]\n",
      "    [1.00000000e-15 1.00000000e+00]\n",
      "    ...\n",
      "    [1.56995650e-01 8.43004350e-01]\n",
      "    [1.12072782e-01 8.87927218e-01]\n",
      "    [7.43933641e-01 2.56066359e-01]]\n",
      "\n",
      "   [[1.00000000e+00 1.00000000e-15]\n",
      "    [1.00000000e-15 1.00000000e+00]\n",
      "    [1.00000000e-15 1.00000000e+00]\n",
      "    ...\n",
      "    [1.56995650e-01 8.43004350e-01]\n",
      "    [1.12072782e-01 8.87927218e-01]\n",
      "    [7.43933641e-01 2.56066359e-01]]]]]\n"
     ]
    }
   ],
   "source": [
    "print(calc_result.X.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `equilibrium()` results\n",
    "\n",
    "The Datasets returned by equilibrium are very similar to calculate, however there are several key differences worth discussing. In equilibrium Datasets, there are six dimensions/coordinates:\n",
    "\n",
    "* `P`: pressures (in Pa).\n",
    "* `T`: temperatures (in K).\n",
    "* `component`: (Same as calculate) The string names of the components in the system.\n",
    "* `internal_dof`: (Same as calculate, except it will be the longest possible internal_dof for all phases) The internal_dof (internal degrees of freedom) is the index of the site in any phase's site fraction array. Below the FCC_A1 phase has the sublattice model (AL, ZN) and thus the internal_dof are integers 0 and 1 referring to the AL site (index 0) and the ZN site (index 1).\n",
    "* `X_ZN`: This is the composition of the species that was passed into the conditions array. Since we passed `v.X('ZN')` to the conditions dictionary, this is `X_ZN`. \n",
    "* `vertex`: The vertex is the index of the phase in equilibrium. The vertex has no inherent physical meaning. There will automatically be enough to describe the number of phases present in any equilibria calculated, implying that vertex can never be large enough to invalidate Gibbs phase rule.\n",
    "\n",
    "There are also at least six Data variables:\n",
    "\n",
    "* `Phase`: The string name of the phase in equilibrium at the conditions. There are as many as `len(vertex)` phases. Any time there are fewer phases in equilibrium than the indices described by `vertex`, the values of phase are paded by `''`, e.g. for a single phase region for FCC_A1, the values of Phase will be `['FCC_A1', '']`. When more than one phase is present, it is important to note that they are not necessarily sorted.\n",
    "* `NP`: Phase fraction of each phase in equilibrium. When there is no other equilibrium phase (e.g. single phase `['FCC_A1', '']`) then the value of `NP` will be `nan` for the absence of a phase, rather than 0.\n",
    "* `MU`: The chemical potentials of each component for the conditions calculated.\n",
    "* `X`: The equilibrium composition of each element in each phase for the calculated conditions.\n",
    "* `Y`: The equilibrium site fraction of each site in each phase for the calculated conditions.\n",
    "* `GM`: Same as `output` for `calculate`. It is always reported no matter the value of `output`.\n",
    "* `output`: (optional) \"output\" is always whatever equilibrium property is calculated by the output keyword passed to `equilibrium`. Unlike `calculate`, this will be in addition to the `GM` because `GM` is always reported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.Dataset>\n",
      "Dimensions:    (N: 1, P: 1, T: 5, X_ZN: 20, component: 2, internal_dof: 2, vertex: 3)\n",
      "Coordinates:\n",
      "  * N          (N) float64 1.0\n",
      "  * P          (P) float64 1.013e+05\n",
      "  * T          (T) float64 500.0 600.0 700.0 800.0 900.0\n",
      "  * X_ZN       (X_ZN) float64 1e-12 0.05 0.1 0.15 0.2 ... 0.75 0.8 0.85 0.9 0.95\n",
      "  * vertex     (vertex) int64 0 1 2\n",
      "  * component  (component) <U2 'AL' 'ZN'\n",
      "Dimensions without coordinates: internal_dof\n",
      "Data variables:\n",
      "    NP         (N, P, T, X_ZN, vertex) float64 1.0 nan nan 1.0 ... 1.0 nan nan\n",
      "    GM         (N, P, T, X_ZN) float64 -1.559e+04 -1.615e+04 ... -5.068e+04\n",
      "    MU         (N, P, T, X_ZN, component) float64 -1.559e+04 ... -5.065e+04\n",
      "    X          (N, P, T, X_ZN, vertex, component) float64 1.0 1e-12 ... nan nan\n",
      "    Y          (N, P, T, X_ZN, vertex, internal_dof) float64 1.0 1e-12 ... nan\n",
      "    Phase      (N, P, T, X_ZN, vertex) <U6 'FCC_A1' '' '' ... 'LIQUID' '' ''\n",
      "    HM         (N, P, T, X_ZN) float64 5.194e+03 5.859e+03 ... 2.528e+04\n",
      "Attributes:\n",
      "    engine:   pycalphad 0.8.3+10.gfd19517e\n",
      "    created:  2020-10-27T14:30:03.243487\n"
     ]
    }
   ],
   "source": [
    "phases = ['LIQUID', 'FCC_A1', 'HCP_A3']\n",
    "eq_result = equilibrium(dbf, comps , phases, {v.X('ZN'):(0,1,0.05), v.T: (500, 1000, 100), v.P:101325}, output='HM')\n",
    "print(eq_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A common operation might be to find the phase fractions of the HCP_A3 phase as a function of composition for T=800.\n",
    "\n",
    "However, the only way we can access the values of the phase fraction is by either the indices or values of the coordinates, we would have to know which index the HCP_A3 phase is in before hand to use the `sel` or `isel` commands. \n",
    "\n",
    "Since we do not know this, we can do what is called [masking](http://xarray.pydata.org/en/stable/indexing.html#masking-with-where) to find the data values that match a condition (the Phase is FCC_A1):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<xarray.DataArray 'NP' (N: 1, X_ZN: 20, vertex: 3)>\n",
      "array([[[1.        ,        nan,        nan],\n",
      "        [1.        ,        nan,        nan],\n",
      "        [1.        ,        nan,        nan],\n",
      "        [1.        ,        nan,        nan],\n",
      "        [0.89739918,        nan,        nan],\n",
      "        [0.71825009,        nan,        nan],\n",
      "        [0.53910097,        nan,        nan],\n",
      "        [0.35995186,        nan,        nan],\n",
      "        [0.18080276,        nan,        nan],\n",
      "        [0.00165369,        nan,        nan],\n",
      "        [       nan,        nan,        nan],\n",
      "        [       nan,        nan,        nan],\n",
      "        [       nan,        nan,        nan],\n",
      "        [       nan,        nan,        nan],\n",
      "        [       nan,        nan,        nan],\n",
      "        [       nan,        nan,        nan],\n",
      "        [       nan,        nan,        nan],\n",
      "        [       nan,        nan,        nan],\n",
      "        [       nan,        nan,        nan],\n",
      "        [       nan,        nan,        nan]]])\n",
      "Coordinates:\n",
      "  * N        (N) float64 1.0\n",
      "    P        float64 1.013e+05\n",
      "    T        float64 800.0\n",
      "  * X_ZN     (X_ZN) float64 1e-12 0.05 0.1 0.15 0.2 ... 0.75 0.8 0.85 0.9 0.95\n",
      "  * vertex   (vertex) int64 0 1 2\n"
     ]
    }
   ],
   "source": [
    "print(eq_result.NP.where(eq_result.Phase=='FCC_A1').sel(P=101325, T=800))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:calphad-dev-2]",
   "language": "python",
   "name": "conda-env-calphad-dev-2-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}