{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "0", "metadata": { "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "import logging\n", "import warnings\n", "\n", "import cmomy\n", "import numpy as np\n", "\n", "rng = cmomy.random.default_rng(0)\n", "\n", "np.set_printoptions(precision=4)\n", "warnings.filterwarnings(\"ignore\")\n", "\n", "\n", "logger = logging.getLogger()\n", "logger.setLevel(logging.ERROR)" ] }, { "cell_type": "markdown", "id": "1", "metadata": {}, "source": [ "# Data Organization\n", "\n", "\n", "All of the extrapolation and interpolation models in the thermoextrap package expect input data to be organized in a certain fashion. To help manage the data, there are data objects that help organize it. Even the inputs to these data object, however, must be organized appropriately. Here we will use data from our ideal gas test system (from {mod}`thermoextrap.idealgas`) to demonstrate this organization, as well as the various options for what types of data that may be provided as input." ] }, { "cell_type": "code", "execution_count": 2, "id": "2", "metadata": {}, "outputs": [], "source": [ "# need imports\n", "%matplotlib inline\n", "import numpy as np\n", "import xarray as xr\n", "\n", "# import thermoextrap\n", "import thermoextrap as xtrap\n", "\n", "# Import idealgas module\n", "from thermoextrap import idealgas" ] }, { "cell_type": "code", "execution_count": 3, "id": "3", "metadata": {}, "outputs": [], "source": [ "# Define reference beta\n", "beta_ref = 5.6\n", "\n", "# And maximum order\n", "order = 2\n", "\n", "npart = 1000 # Number of particles (in single configuration)\n", "nconfig = 100_000 # Number of configurations\n", "\n", "# Generate all the data we could want\n", "xdata, udata = idealgas.generate_data((nconfig, npart), beta_ref)" ] }, { "cell_type": "markdown", "id": "4", "metadata": {}, "source": [ "Refer to {mod}`thermoextrap.data` for more information on the data classes." ] }, { "cell_type": "markdown", "id": "5", "metadata": {}, "source": [ "## Basics\n", "\n", "Rather than passing data directly to `__init__` methods for creating data class objects and simultaneously telling it which dimensions mean what (or expecting that specific dimensions mean a certain thing), {mod}`thermoextrap` uses {mod}`xarray` to label the dimensions of inputs. While this is also useful in the background, it helps to clarify what is expected of user inputs.\n", "\n", "Currently, `xdata` is of the shape (nconfig), or the number of configurations generated with each entry being the average $x$ location for the associated configuration." ] }, { "cell_type": "code", "execution_count": 4, "id": "6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "100000 (100000,)\n" ] } ], "source": [ "print(nconfig, xdata.shape)" ] }, { "cell_type": "markdown", "id": "7", "metadata": {}, "source": [ "The dimension over which independent samples vary is the \"record\" dimension, with its default name in {mod}`thermoextrap.data` being 'rec'. So when we create an {class}`xarray.DataArray` object to house the input $x$ data, we must label that dimension 'rec'. Same goes for the input potential energy data. Note that the list provided to the argument `dims` is a list of strings naming the dimensions in the array passed to {class}`xarray.DataArray`." ] }, { "cell_type": "code", "execution_count": 5, "id": "8", "metadata": {}, "outputs": [], "source": [ "xdata = xr.DataArray(xdata, dims=[\"rec\"])\n", "udata = xr.DataArray(udata, dims=[\"rec\"])" ] }, { "cell_type": "markdown", "id": "9", "metadata": {}, "source": [ "Now when we create a data object in {mod}`thermoextrap` to hold the data, we tell it that the \"record\" dimension, `rec_dim` is named 'rec', which is the default, but it could be named something different as long as you provided that name to `rec_dim`.\n", "\n", "Note that the `xv` is the argument for the observable $x$ and `uv` is the potential energy or appropriate Hamiltonian or thermodynamic conjugate variable." ] }, { "cell_type": "code", "execution_count": 6, "id": "10", "metadata": {}, "outputs": [], "source": [ "data = xtrap.DataCentralMomentsVals.from_vals(\n", " order=order, rec_dim=\"rec\", xv=xdata, uv=udata, central=True\n", ")" ] }, { "cell_type": "markdown", "id": "11", "metadata": {}, "source": [ "A couple more notes are in order about the inputs to any of the {mod}`thermoextrap.data` object variants. First, you only need to provide the order you expect to extrapolate to up front if you're using the {meth}`~thermoextrap.data.DataCentralMomentsVals.from_vals` constructor. This is because you need to specify the order of moments that will be calculated from the raw data.\n", "\n", "The next argument to be aware of is `central`. This is True by default and tells the data object to work with central moments for calculating derivatives in the background, which it turns out is much more numerically stable than non-central moments. You probably want `central` to be True, but know that you can change it if you wish." ] }, { "cell_type": "markdown", "id": "12", "metadata": {}, "source": [ "## Data Structure\n", "\n", "```{eval-rst}\n", ".. currentmodule:: thermoextrap.data\n", "```\n", "\n", "\n", "\n", "A lot of data is already computed as soon as we create our data object. The original raw data is still stored in {attr}`~DataCentralMomentsVals.xv` and {attr}`~DataCentralMomentsVals.uv`, and order is {attr}`~DataCentralMomentsVals.order`, but we can already see the central moments appearing if we look at..." ] }, { "cell_type": "code", "execution_count": 7, "id": "13", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (rec: 100000)> Size: 800kB\n",
       "array([0.181 , 0.1687, 0.1713, ..., 0.1763, 0.1761, 0.1829])\n",
       "Dimensions without coordinates: rec
" ], "text/plain": [ " Size: 800kB\n", "array([0.181 , 0.1687, 0.1713, ..., 0.1763, 0.1761, 0.1829])\n", "Dimensions without coordinates: rec" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.xv" ] }, { "cell_type": "code", "execution_count": 8, "id": "14", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray ()> Size: 8B\n",
       "array(0.1748)
" ], "text/plain": [ " Size: 8B\n", "array(0.1748)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.xave" ] }, { "cell_type": "markdown", "id": "15", "metadata": {}, "source": [ "{attr}`~DataCentralMomentsVals.xave` is the average observable value." ] }, { "cell_type": "code", "execution_count": 9, "id": "16", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (umom: 3)> Size: 24B\n",
       "array([1.0000e+00, 1.7485e+02, 3.0601e+04])\n",
       "Dimensions without coordinates: umom
" ], "text/plain": [ " Size: 24B\n", "array([1.0000e+00, 1.7485e+02, 3.0601e+04])\n", "Dimensions without coordinates: umom" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.u" ] }, { "cell_type": "markdown", "id": "17", "metadata": {}, "source": [ "{attr}`~DataCentralMomentsVals.u` are the moments of the potential energy, i.e., $\\langle U^i \\rangle$ for the $i^\\mathrm{th}$ index in the array.\n", "\n", "For the central moments of the potential energy, $\\langle (U - \\langle U \\rangle )^i \\rangle$, you can look at {attr}`~DataCentralMomentsVals.du`" ] }, { "cell_type": "code", "execution_count": 10, "id": "18", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (umom: 3)> Size: 24B\n",
       "array([ 1.    ,  0.    , 28.2438])\n",
       "Dimensions without coordinates: umom
" ], "text/plain": [ " Size: 24B\n", "array([ 1. , 0. , 28.2438])\n", "Dimensions without coordinates: umom" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.du" ] }, { "cell_type": "markdown", "id": "19", "metadata": {}, "source": [ "The other necessary component for calculating derivatives is $\\langle x U^i \\rangle$, which is in {attr}`DataCentralMomentsVals.xu`" ] }, { "cell_type": "code", "execution_count": 11, "id": "20", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (umom: 3)> Size: 24B\n",
       "array([1.7485e-01, 3.0601e+01, 5.3604e+03])\n",
       "Dimensions without coordinates: umom
" ], "text/plain": [ " Size: 24B\n", "array([1.7485e-01, 3.0601e+01, 5.3604e+03])\n", "Dimensions without coordinates: umom" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.xu" ] }, { "cell_type": "markdown", "id": "21", "metadata": {}, "source": [ "Or if working with central moments, $\\langle (x - \\langle x \\rangle) (U - \\langle U \\rangle)^i \\rangle$, is in {attr}`DataCentralMomentsVals.dxdu`" ] }, { "cell_type": "code", "execution_count": 12, "id": "22", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (umom: 3)> Size: 24B\n",
       "array([0.    , 0.0282, 0.0078])\n",
       "Dimensions without coordinates: umom
" ], "text/plain": [ " Size: 24B\n", "array([0. , 0.0282, 0.0078])\n", "Dimensions without coordinates: umom" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.dxdu" ] }, { "cell_type": "markdown", "id": "23", "metadata": {}, "source": [ "All of this information is condensed in {attr}`~DataCentralMomentsVals.values`, which takes only exactly what we need for computing derivatives." ] }, { "cell_type": "code", "execution_count": 13, "id": "24", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (xmom: 2, umom: 3)> Size: 48B\n",
       "array([[1.0000e+05, 1.7485e+02, 2.8244e+01],\n",
       "       [1.7485e-01, 2.8244e-02, 7.8139e-03]])\n",
       "Dimensions without coordinates: xmom, umom
" ], "text/plain": [ " Size: 48B\n", "array([[1.0000e+05, 1.7485e+02, 2.8244e+01],\n", " [1.7485e-01, 2.8244e-02, 7.8139e-03]])\n", "Dimensions without coordinates: xmom, umom" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.values" ] }, { "cell_type": "markdown", "id": "25", "metadata": {}, "source": [ "Understanding this internal structure will help to understand possible inputs as well. In {attr}`~DataCentralMoments.values`, the data object has stored the central moments of $\\langle x \\rangle$ and $\\langle U^i \\rangle$ and cross moments $\\langle x U^i \\rangle$. Note that the second dimension, \"umom\", short for \"U moments\", is just `order` plus 1. That makes sense if we remember that the zeroth order derivative is the observable itself and we specify that we want to use up to `order` derivatives. So the second dimension involves setting the exponent $i$ on $U$ equal to the index of that dimension, or order of that moment. The first dimension does the same thing, but with $x$. Regardless of the order, however, we only ever need $x$ raised to the zeroth or first power in the average.\n", "\n", "The first row in `values` contains all moments of just $U$, i.e., $\\langle U^0 \\rangle$, $\\langle U^1 \\rangle$, $\\langle U^2 \\rangle$, etc. The second row in `values` contains all moments of $x$ multiplied by $U$, i.e., $\\langle x U^0 \\rangle$, $\\langle x U^1 \\rangle$, etc. But note that beyond the powers of 0 and 1 for the first row, and just 0 for the second row, all values shown are central moments, e.g., $\\langle (x - \\langle x \\rangle)(U - \\langle U \\rangle)^i \\rangle$ or $\\langle (U - \\langle U \\rangle)^i \\rangle$." ] }, { "cell_type": "markdown", "id": "26", "metadata": {}, "source": [ "In other words, {attr}`~DataCentralMomentsVals.values` is a special array with structure...\n", "\n", "for i + j <= 1:
\n", "`data.values[0, 0]` = {sum of weights or count}
\n", "`data.values[1, 0]` = {ave of x} = $\\langle x \\rangle$
\n", "`data.values[0, 1]` = {ave of u} = $\\langle U \\rangle$
\n", "\n", "for i + j > 1:
\n", "`data.values[i, j]` = $\\langle (x - \\langle x \\rangle)^i (U - \\langle U \\rangle)^j \\rangle$\n", "\n", "To summarize, {attr}`~DataCentralMomentsVals.values` contains the bare bones of what is required for calculating derivatives and will be shared in some form or another across all data classes, with this information passed to the functions that compute derivative." ] }, { "cell_type": "markdown", "id": "27", "metadata": {}, "source": [ "## Input formats and resampling\n", "\n", "Since {attr}`~DataCentralMomentsVals.values` reflects the internal structure, you can just provide it directly (or something similar to it in terms of moments) if you prefer to do that. You'll just need to use a different data class {class}`DataCentralMoments` and a different constructor method {meth}`DataCentralMoments.from_vals`.\n", "\n", "While {class}`DataCentralMomentsVals` is designed to work with 'values' (i.e., individual observations), {class}`DataCentralMoments` is designed to work with moments. Both classes can be constructed from 'values', but {class}`DataCentralMomentsVals` retains the underlying values (for resampling, etc), {class}`DataCentralMoments` converts the values to moments, and goes from there. Basically, if you have pre-computed moments (e.g., from a simulation), {class}`DataCentralMoments` is probably what you want to use. Note that resampling for {class}`DataCentralMoments` is based on resampling over multiple samples the moments." ] }, { "cell_type": "markdown", "id": "28", "metadata": {}, "source": [ "For example, if we construct an {class}`DataCentralMoments` object using the {meth}`DataCentralMoments.from_data` constructor, we have:" ] }, { "cell_type": "code", "execution_count": 14, "id": "29", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (xmom: 2, umom: 3)> Size: 48B\n",
       "array([[1.0000e+05, 1.7485e+02, 2.8244e+01],\n",
       "       [1.7485e-01, 2.8244e-02, 7.8139e-03]])\n",
       "Dimensions without coordinates: xmom, umom
" ], "text/plain": [ " Size: 48B\n", "array([[1.0000e+05, 1.7485e+02, 2.8244e+01],\n", " [1.7485e-01, 2.8244e-02, 7.8139e-03]])\n", "Dimensions without coordinates: xmom, umom" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data_noboot = xtrap.DataCentralMoments.from_data(data.values)\n", "xr.testing.assert_allclose(data_noboot.values, data.values)\n", "data_noboot.values" ] }, { "cell_type": "markdown", "id": "30", "metadata": {}, "source": [ "Which is identical to `data` above. Note that the order here is inferred from the passed moments array. Likewise, we could have just created this from values using {meth}`DataCentralMoments.from_vals`" ] }, { "cell_type": "code", "execution_count": 15, "id": "31", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (xmom: 2, umom: 3)> Size: 48B\n",
       "array([[1.0000e+05, 1.7485e+02, 2.8244e+01],\n",
       "       [1.7485e-01, 2.8244e-02, 7.8139e-03]])\n",
       "Dimensions without coordinates: xmom, umom
" ], "text/plain": [ " Size: 48B\n", "array([[1.0000e+05, 1.7485e+02, 2.8244e+01],\n", " [1.7485e-01, 2.8244e-02, 7.8139e-03]])\n", "Dimensions without coordinates: xmom, umom" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data_noboot = xtrap.DataCentralMoments.from_vals(\n", " xv=xdata, uv=udata, rec_dim=\"rec\", central=True, order=order\n", ")\n", "xr.testing.assert_allclose(data_noboot.values, data.values)\n", "data_noboot.values" ] }, { "cell_type": "markdown", "id": "32", "metadata": {}, "source": [ "However, since `data_noboot` is based on just a single average, bootstrapping makes little sense. For example:" ] }, { "cell_type": "code", "execution_count": 16, "id": "33", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "caught error!\n", "not implemented for scalar\n" ] } ], "source": [ "data_noboot = xtrap.DataCentralMoments.from_raw(data.values)\n", "\n", "try:\n", " data_noboot.resample(nrep=3).values\n", "except ValueError as e:\n", " print(\"caught error!\")\n", " print(e)" ] }, { "cell_type": "markdown", "id": "34", "metadata": {}, "source": [ "versus..." ] }, { "cell_type": "code", "execution_count": 17, "id": "35", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (rep: 3, xmom: 2, umom: 3)> Size: 144B\n",
       "array([[[1.0000e+05, 1.7485e+02, 2.8276e+01],\n",
       "        [1.7485e-01, 2.8276e-02, 7.8536e-03]],\n",
       "\n",
       "       [[1.0000e+05, 1.7483e+02, 2.8195e+01],\n",
       "        [1.7483e-01, 2.8195e-02, 6.9730e-03]],\n",
       "\n",
       "       [[1.0000e+05, 1.7487e+02, 2.8451e+01],\n",
       "        [1.7487e-01, 2.8451e-02, 9.0861e-03]]])\n",
       "Dimensions without coordinates: rep, xmom, umom
" ], "text/plain": [ " Size: 144B\n", "array([[[1.0000e+05, 1.7485e+02, 2.8276e+01],\n", " [1.7485e-01, 2.8276e-02, 7.8536e-03]],\n", "\n", " [[1.0000e+05, 1.7483e+02, 2.8195e+01],\n", " [1.7483e-01, 2.8195e-02, 6.9730e-03]],\n", "\n", " [[1.0000e+05, 1.7487e+02, 2.8451e+01],\n", " [1.7487e-01, 2.8451e-02, 9.0861e-03]]])\n", "Dimensions without coordinates: rep, xmom, umom" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.resample(nrep=3).values" ] }, { "cell_type": "markdown", "id": "36", "metadata": {}, "source": [ "Note that whenever you call `.resample()` a new dimension is created for the output called 'rep', short for \"repetitions.\" This is similar to the 'rec' dimension, but helps you keep track of whether you're working with original data, or a bootstrapped sample." ] }, { "cell_type": "markdown", "id": "37", "metadata": {}, "source": [ "So clearly if you like to provide raw moments and the uncertainty quantification in `thermoextrap` to work, you will need to do this from blocks of your data or from repeated simulations. But the code will all still work if you prefer to calculate your own moments (saving them periodically from simulations rather than saving all of the configurations, energies, or observables frequently).\n", "\n", "As an example, we can work with blocks of data, adding an axis called 'block' that we will ask the constructor to average over by specifying the `dim` argument." ] }, { "cell_type": "code", "execution_count": 18, "id": "38", "metadata": {}, "outputs": [], "source": [ "# Make 100 averaged observations\n", "xx = xr.DataArray(xdata.values.reshape(100, -1), dims=[\"rec\", \"block\"])\n", "uu = xr.DataArray(udata.values.reshape(100, -1), dims=[\"rec\", \"block\"])" ] }, { "cell_type": "code", "execution_count": 19, "id": "39", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (rec: 100, xmom: 2, umom: 3)> Size: 5kB\n",
       "array([[[ 1.0000e+03,  1.7504e+02,  2.8681e+01],\n",
       "        [ 1.7504e-01,  2.8681e-02,  1.0778e-02]],\n",
       "\n",
       "       [[ 1.0000e+03,  1.7487e+02,  2.9507e+01],\n",
       "        [ 1.7487e-01,  2.9507e-02,  3.6761e-03]],\n",
       "\n",
       "       [[ 1.0000e+03,  1.7471e+02,  2.8956e+01],\n",
       "        [ 1.7471e-01,  2.8956e-02,  3.4057e-02]],\n",
       "\n",
       "       [[ 1.0000e+03,  1.7470e+02,  2.9232e+01],\n",
       "        [ 1.7470e-01,  2.9232e-02, -6.4238e-03]],\n",
       "\n",
       "       [[ 1.0000e+03,  1.7469e+02,  2.8072e+01],\n",
       "        [ 1.7469e-01,  2.8072e-02,  5.0937e-03]],\n",
       "\n",
       "       [[ 1.0000e+03,  1.7457e+02,  2.6426e+01],\n",
       "        [ 1.7457e-01,  2.6426e-02,  3.2884e-03]],\n",
       "\n",
       "       [[ 1.0000e+03,  1.7480e+02,  2.8580e+01],\n",
       "        [ 1.7480e-01,  2.8580e-02,  2.9648e-03]],\n",
       "...\n",
       "       [[ 1.0000e+03,  1.7491e+02,  2.8840e+01],\n",
       "        [ 1.7491e-01,  2.8840e-02,  1.4932e-02]],\n",
       "\n",
       "       [[ 1.0000e+03,  1.7478e+02,  2.7380e+01],\n",
       "        [ 1.7478e-01,  2.7380e-02,  2.4111e-02]],\n",
       "\n",
       "       [[ 1.0000e+03,  1.7492e+02,  2.6600e+01],\n",
       "        [ 1.7492e-01,  2.6600e-02, -3.0268e-03]],\n",
       "\n",
       "       [[ 1.0000e+03,  1.7489e+02,  2.8386e+01],\n",
       "        [ 1.7489e-01,  2.8386e-02, -8.0114e-03]],\n",
       "\n",
       "       [[ 1.0000e+03,  1.7477e+02,  2.6906e+01],\n",
       "        [ 1.7477e-01,  2.6906e-02, -2.0863e-04]],\n",
       "\n",
       "       [[ 1.0000e+03,  1.7486e+02,  2.6945e+01],\n",
       "        [ 1.7486e-01,  2.6945e-02,  2.7821e-02]],\n",
       "\n",
       "       [[ 1.0000e+03,  1.7496e+02,  2.7179e+01],\n",
       "        [ 1.7496e-01,  2.7179e-02, -1.3813e-02]]])\n",
       "Dimensions without coordinates: rec, xmom, umom
" ], "text/plain": [ " Size: 5kB\n", "array([[[ 1.0000e+03, 1.7504e+02, 2.8681e+01],\n", " [ 1.7504e-01, 2.8681e-02, 1.0778e-02]],\n", "\n", " [[ 1.0000e+03, 1.7487e+02, 2.9507e+01],\n", " [ 1.7487e-01, 2.9507e-02, 3.6761e-03]],\n", "\n", " [[ 1.0000e+03, 1.7471e+02, 2.8956e+01],\n", " [ 1.7471e-01, 2.8956e-02, 3.4057e-02]],\n", "\n", " [[ 1.0000e+03, 1.7470e+02, 2.9232e+01],\n", " [ 1.7470e-01, 2.9232e-02, -6.4238e-03]],\n", "\n", " [[ 1.0000e+03, 1.7469e+02, 2.8072e+01],\n", " [ 1.7469e-01, 2.8072e-02, 5.0937e-03]],\n", "\n", " [[ 1.0000e+03, 1.7457e+02, 2.6426e+01],\n", " [ 1.7457e-01, 2.6426e-02, 3.2884e-03]],\n", "\n", " [[ 1.0000e+03, 1.7480e+02, 2.8580e+01],\n", " [ 1.7480e-01, 2.8580e-02, 2.9648e-03]],\n", "...\n", " [[ 1.0000e+03, 1.7491e+02, 2.8840e+01],\n", " [ 1.7491e-01, 2.8840e-02, 1.4932e-02]],\n", "\n", " [[ 1.0000e+03, 1.7478e+02, 2.7380e+01],\n", " [ 1.7478e-01, 2.7380e-02, 2.4111e-02]],\n", "\n", " [[ 1.0000e+03, 1.7492e+02, 2.6600e+01],\n", " [ 1.7492e-01, 2.6600e-02, -3.0268e-03]],\n", "\n", " [[ 1.0000e+03, 1.7489e+02, 2.8386e+01],\n", " [ 1.7489e-01, 2.8386e-02, -8.0114e-03]],\n", "\n", " [[ 1.0000e+03, 1.7477e+02, 2.6906e+01],\n", " [ 1.7477e-01, 2.6906e-02, -2.0863e-04]],\n", "\n", " [[ 1.0000e+03, 1.7486e+02, 2.6945e+01],\n", " [ 1.7486e-01, 2.6945e-02, 2.7821e-02]],\n", "\n", " [[ 1.0000e+03, 1.7496e+02, 2.7179e+01],\n", " [ 1.7496e-01, 2.7179e-02, -1.3813e-02]]])\n", "Dimensions without coordinates: rec, xmom, umom" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Create directly from values of moments - notice that this is DataCentralMoments, not DataCentralMomentsVals\n", "# Effectively just means that the 'rec' dim will not be collapsed when using data for extrapolation, etc.\n", "# Behaves similarly to the 'rep' dim when resampling\n", "data_fv = xtrap.DataCentralMoments.from_vals(\n", " xv=xx, uv=uu, dim=\"block\", order=order, central=True\n", ")\n", "# So 'rec' is for each separate block average\n", "data_fv.values" ] }, { "cell_type": "markdown", "id": "40", "metadata": {}, "source": [ "Again, note that we did not use the {class}`DataCentralMomentsVals` class above, but instead the {class}`DataCentralMoments` class. The former is for processing and storing simulation \"timeseries\" of observable and potential energy values, while the latter takes pre-computed moments, including multiple replicates of precomputed moments as above. Behind the scenes, this will influence how bootstrapped confidence intervals are computed.\n", "\n", "What's functionally different, though, is that the 'rec' dim also appears in `.values`. That means that when this data is used in models for extrapolation or interpolation, that dimension will also be preserved. So prediction to a new $\\beta$ value will result in an output matching the size of the 'rec' dimension in the same way that it would match the 'rep' dimension created through resampling.\n", "\n", "If we resample over this data set, we see that we just take `nrep` random samples from it, putting those samples into a new dimension called 'rep'." ] }, { "cell_type": "code", "execution_count": 20, "id": "41", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (rep: 3, xmom: 2, umom: 3)> Size: 144B\n",
       "array([[[1.0000e+05, 1.7485e+02, 2.8257e+01],\n",
       "        [1.7485e-01, 2.8257e-02, 8.0428e-03]],\n",
       "\n",
       "       [[1.0000e+05, 1.7485e+02, 2.8462e+01],\n",
       "        [1.7485e-01, 2.8462e-02, 7.3831e-03]],\n",
       "\n",
       "       [[1.0000e+05, 1.7485e+02, 2.8340e+01],\n",
       "        [1.7485e-01, 2.8340e-02, 9.6959e-03]]])\n",
       "Dimensions without coordinates: rep, xmom, umom
" ], "text/plain": [ " Size: 144B\n", "array([[[1.0000e+05, 1.7485e+02, 2.8257e+01],\n", " [1.7485e-01, 2.8257e-02, 8.0428e-03]],\n", "\n", " [[1.0000e+05, 1.7485e+02, 2.8462e+01],\n", " [1.7485e-01, 2.8462e-02, 7.3831e-03]],\n", "\n", " [[1.0000e+05, 1.7485e+02, 2.8340e+01],\n", " [1.7485e-01, 2.8340e-02, 9.6959e-03]]])\n", "Dimensions without coordinates: rep, xmom, umom" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data_fv.resample(nrep=3).values" ] }, { "cell_type": "markdown", "id": "42", "metadata": {}, "source": [ "If we had computed the moments from the blocked data ourselves, we could also create a data object with the {meth}`DataCentralMoments.from_ave_raw` constructor (below). Many other constructors exist, including from central moments if you like. If you use those, please take a look at the documentation to make sure you are specifying or using the correct dimension naming conventions, such as 'rec', 'xmom', 'umom', etc. Remember, if you are extrapolating an observable that has an explicit dependence on the extrapolation observable, you also need to specify the 'deriv' dimension that describes the observed derivatives with respect to the extrapolation variable (see the [Temperature extrapolation case 2](./Temperature_Extrap_Case2.ipynb) notebook)." ] }, { "cell_type": "code", "execution_count": 21, "id": "43", "metadata": {}, "outputs": [], "source": [ "# Compute moments of U, i.e., averages to integer powers up to maximum order desired\n", "mom_u = xr.DataArray(np.arange(order + 1), dims=[\"umom\"])\n", "uave = (uu**mom_u).mean(\"block\")\n", "xuave = (xx * uu**mom_u).mean(\"block\")\n", "data_fa = xtrap.DataCentralMoments.from_ave_raw(\n", " u=uave, xu=xuave, central=True, w=xx.sizes[\"block\"]\n", ")\n", "\n", "xr.testing.assert_allclose(data_fv.values, data_fa.values)" ] }, { "cell_type": "markdown", "id": "44", "metadata": {}, "source": [ "The above `.values` should be identical to those from the `from_vals` constructor.\n", "\n", "At this point, we have seen how the same data objects that interface with extrapolation or interpolation models can be created from different inputs. Other features also exist, such as specifying weights with the argument `w` to a constructor to change the weights used during averaging." ] }, { "cell_type": "markdown", "id": "45", "metadata": {}, "source": [ "## Vector observables\n", "\n", "Finally, we can also have vector observables, like RDFs, for example. This is easy to accomplish with any of the above constructors or types of data input. All that is required is to add another dimension to our {class}`xarray.DataArray` input. Typically, we will call this dimension 'vals', short for \"values\", which is the default name for this dimension when using the {meth}`DataCentralMomentsVals.from_vals` constructor." ] }, { "cell_type": "code", "execution_count": 22, "id": "46", "metadata": {}, "outputs": [], "source": [ "# Extrapolate both average x and average x**2\n", "x_xsq_data = xr.DataArray(\n", " np.vstack([xdata.values, xdata.values**2]).T,\n", " dims=[\"rec\", \"vals\"],\n", " coords={\"vals\": [\"x\", \"xsq\"]},\n", ")\n", "data_vec = xtrap.DataCentralMomentsVals.from_vals(\n", " order=order, rec_dim=\"rec\", xv=x_xsq_data, uv=udata, central=True\n", ")" ] }, { "cell_type": "code", "execution_count": 23, "id": "47", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (vals: 2, xmom: 2, umom: 3)> Size: 96B\n",
       "array([[[1.0000e+05, 1.7485e+02, 2.8244e+01],\n",
       "        [1.7485e-01, 2.8244e-02, 7.8139e-03]],\n",
       "\n",
       "       [[1.0000e+05, 1.7485e+02, 2.8244e+01],\n",
       "        [3.0601e-02, 9.8847e-03, 4.3329e-03]]])\n",
       "Coordinates:\n",
       "  * vals     (vals) <U3 24B 'x' 'xsq'\n",
       "Dimensions without coordinates: xmom, umom
" ], "text/plain": [ " Size: 96B\n", "array([[[1.0000e+05, 1.7485e+02, 2.8244e+01],\n", " [1.7485e-01, 2.8244e-02, 7.8139e-03]],\n", "\n", " [[1.0000e+05, 1.7485e+02, 2.8244e+01],\n", " [3.0601e-02, 9.8847e-03, 4.3329e-03]]])\n", "Coordinates:\n", " * vals (vals) \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (rep: 3, vals: 2, xmom: 2, umom: 3)> Size: 288B\n",
       "array([[[[1.0000e+05, 1.7484e+02, 2.8261e+01],\n",
       "         [1.7484e-01, 2.8261e-02, 7.9176e-03]],\n",
       "\n",
       "        [[1.0000e+05, 1.7484e+02, 2.8261e+01],\n",
       "         [3.0598e-02, 9.8903e-03, 4.3704e-03]]],\n",
       "\n",
       "\n",
       "       [[[1.0000e+05, 1.7482e+02, 2.8255e+01],\n",
       "         [1.7482e-01, 2.8255e-02, 7.9522e-03]],\n",
       "\n",
       "        [[1.0000e+05, 1.7482e+02, 2.8255e+01],\n",
       "         [3.0590e-02, 9.8870e-03, 4.3943e-03]]],\n",
       "\n",
       "\n",
       "       [[[1.0000e+05, 1.7485e+02, 2.8214e+01],\n",
       "         [1.7485e-01, 2.8214e-02, 6.9934e-03]],\n",
       "\n",
       "        [[1.0000e+05, 1.7485e+02, 2.8214e+01],\n",
       "         [3.0601e-02, 9.8735e-03, 4.0474e-03]]]])\n",
       "Coordinates:\n",
       "  * vals     (vals) <U3 24B 'x' 'xsq'\n",
       "Dimensions without coordinates: rep, xmom, umom
" ], "text/plain": [ " Size: 288B\n", "array([[[[1.0000e+05, 1.7484e+02, 2.8261e+01],\n", " [1.7484e-01, 2.8261e-02, 7.9176e-03]],\n", "\n", " [[1.0000e+05, 1.7484e+02, 2.8261e+01],\n", " [3.0598e-02, 9.8903e-03, 4.3704e-03]]],\n", "\n", "\n", " [[[1.0000e+05, 1.7482e+02, 2.8255e+01],\n", " [1.7482e-01, 2.8255e-02, 7.9522e-03]],\n", "\n", " [[1.0000e+05, 1.7482e+02, 2.8255e+01],\n", " [3.0590e-02, 9.8870e-03, 4.3943e-03]]],\n", "\n", "\n", " [[[1.0000e+05, 1.7485e+02, 2.8214e+01],\n", " [1.7485e-01, 2.8214e-02, 6.9934e-03]],\n", "\n", " [[1.0000e+05, 1.7485e+02, 2.8214e+01],\n", " [3.0601e-02, 9.8735e-03, 4.0474e-03]]]])\n", "Coordinates:\n", " * vals (vals) \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (vals: 2, xmom: 2, umom: 3)> Size: 96B\n",
       "array([[[1.0000e+05, 1.7485e+02, 2.8244e+01],\n",
       "        [1.7485e-01, 2.8244e-02, 7.8139e-03]],\n",
       "\n",
       "       [[1.0000e+05, 1.7485e+02, 2.8244e+01],\n",
       "        [3.0601e-02, 9.8847e-03, 4.3329e-03]]])\n",
       "Coordinates:\n",
       "  * vals     (vals) <U3 24B 'x' 'xsq'\n",
       "Dimensions without coordinates: xmom, umom
" ], "text/plain": [ " Size: 96B\n", "array([[[1.0000e+05, 1.7485e+02, 2.8244e+01],\n", " [1.7485e-01, 2.8244e-02, 7.8139e-03]],\n", "\n", " [[1.0000e+05, 1.7485e+02, 2.8244e+01],\n", " [3.0601e-02, 9.8847e-03, 4.3329e-03]]])\n", "Coordinates:\n", " * vals (vals) \n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (rep: 3, vals: 2, xmom: 2, umom: 3)> Size: 288B\n",
       "array([[[[1.0000e+05, 1.7485e+02, 2.8203e+01],\n",
       "         [1.7485e-01, 2.8203e-02, 6.0452e-03]],\n",
       "\n",
       "        [[1.0000e+05, 1.7485e+02, 2.8203e+01],\n",
       "         [3.0600e-02, 9.8685e-03, 3.7024e-03]]],\n",
       "\n",
       "\n",
       "       [[[1.0000e+05, 1.7486e+02, 2.8246e+01],\n",
       "         [1.7486e-01, 2.8246e-02, 8.5305e-03]],\n",
       "\n",
       "        [[1.0000e+05, 1.7486e+02, 2.8246e+01],\n",
       "         [3.0604e-02, 9.8867e-03, 4.5961e-03]]],\n",
       "\n",
       "\n",
       "       [[[1.0000e+05, 1.7486e+02, 2.8203e+01],\n",
       "         [1.7486e-01, 2.8203e-02, 1.1006e-02]],\n",
       "\n",
       "        [[1.0000e+05, 1.7486e+02, 2.8203e+01],\n",
       "         [3.0603e-02, 9.8740e-03, 5.4435e-03]]]])\n",
       "Coordinates:\n",
       "  * vals     (vals) <U3 24B 'x' 'xsq'\n",
       "Dimensions without coordinates: rep, xmom, umom
" ], "text/plain": [ " Size: 288B\n", "array([[[[1.0000e+05, 1.7485e+02, 2.8203e+01],\n", " [1.7485e-01, 2.8203e-02, 6.0452e-03]],\n", "\n", " [[1.0000e+05, 1.7485e+02, 2.8203e+01],\n", " [3.0600e-02, 9.8685e-03, 3.7024e-03]]],\n", "\n", "\n", " [[[1.0000e+05, 1.7486e+02, 2.8246e+01],\n", " [1.7486e-01, 2.8246e-02, 8.5305e-03]],\n", "\n", " [[1.0000e+05, 1.7486e+02, 2.8246e+01],\n", " [3.0604e-02, 9.8867e-03, 4.5961e-03]]],\n", "\n", "\n", " [[[1.0000e+05, 1.7486e+02, 2.8203e+01],\n", " [1.7486e-01, 2.8203e-02, 1.1006e-02]],\n", "\n", " [[1.0000e+05, 1.7486e+02, 2.8203e+01],\n", " [3.0603e-02, 9.8740e-03, 5.4435e-03]]]])\n", "Coordinates:\n", " * vals (vals)