basic_examples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RMWSPy\n",
    "\n",
    "This is a compilation of basic examples for the RMWSPy package for conditional spatial simulation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# open an interactive console that one can use to check variables \n",
    "# or perform calculations interactively\n",
    "%qtconsole \n",
    "# import modules\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import scipy.stats as st\n",
    "import scipy.spatial as sp\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from RMWSPy.rmwspy import *\n",
    "\n",
    "# also let's define a random seed to ensure repeatability\n",
    "np.random.seed(123)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unconditional random fields\n",
    "Let's start with the most basic case of spatial random field simulation - unconditional simulation. We just need to define a domainsize for the field as well as a spatial covariance model. We can then use FFTMA to simulate an unconditional spatial random field with the given parameters. Note that the resulting field will have a standard normal marginal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define field size\n",
    "domainsize = (100, 100)\n",
    "\n",
    "# spatial covariance which we assume to know from our observations\n",
    "# the covariance model is defined as: Sill Model(Range) \n",
    "# and it can be a linear combination of multiple models\n",
    "# available covariance models are: Nug (Nugget), Exp (Exponential), \n",
    "# Sph (Spherical), Gau (Gaussian)\n",
    "covariancemodel = '0.01 Nug(0.0) + 0.99 Exp(7.2)' \n",
    "\n",
    "fftma_sim = FFTMA(domainsize=domainsize, covmod=covariancemodel)\n",
    "field = fftma_sim.simnew()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now let's plot that field\n",
    "plt.figure(figsize=(5, 5))  \n",
    "ax = plt.gca()    \n",
    "im = ax.imshow(field,\n",
    "            origin='lower', interpolation='nearest', \n",
    "            cmap='jet', vmin=-4, vmax=4)\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if you want another realization\n",
    "# you can just run the above code block again\n",
    "# or you call simnew() again\n",
    "field = fftma_sim.simnew()\n",
    "# and plot it\n",
    "plt.figure(figsize=(5, 5))  \n",
    "ax = plt.gca()    \n",
    "im = ax.imshow(field,\n",
    "            origin='lower', interpolation='nearest', \n",
    "            cmap='jet', vmin=-4, vmax=4)\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's simulate a field with a much bigger variogram range\n",
    "# therefore we have to re-initialise FFTMA with a\n",
    "# updated covariancemodel\n",
    "domainsize = (100, 100)\n",
    "covariancemodel = '0.01 Nug(0.0) + 0.99 Exp(14)' \n",
    "\n",
    "fftma_sim = FFTMA(domainsize=domainsize, covmod=covariancemodel)\n",
    "field = fftma_sim.simnew()\n",
    "\n",
    "# and plot it again\n",
    "plt.figure(figsize=(5, 5))  \n",
    "ax = plt.gca()    \n",
    "im = ax.imshow(field,\n",
    "            origin='lower', interpolation='nearest', \n",
    "            cmap='jet', vmin=-4, vmax=4)\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Unconditional fields with anisotropy\n",
    "We can also simulate field with anisotropy. Therefore, in addition to the previously defined parameters, we need to define a scaling and rotation variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "domainsize = (120, 80)\n",
    "covariancemodel = '0.01 Nug(0.0) + 0.99 Exp(14)' \n",
    "anisotropy = (0.2, 1, 25) # x dimension has the above specified range, y dimension has 0.2*range, 25 degrees rotation around the x axis\n",
    "\n",
    "fftma_sim = FFTMA(domainsize=domainsize, covmod=covariancemodel, anisotropy=anisotropy)\n",
    "field = fftma_sim.simnew()\n",
    "\n",
    "# and plot it again\n",
    "plt.figure(figsize=(5, 5))  \n",
    "ax = plt.gca()    \n",
    "im = ax.imshow(field,\n",
    "            origin='lower', interpolation='nearest', \n",
    "            cmap='jet', vmin=-4, vmax=4)\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can do the same for a 3D field\n",
    "domainsize = (80, 60, 50)\n",
    "covmod = '1.0 Sph(14.5)'\n",
    "\n",
    "fftma_sim = FFTMA(domainsize=domainsize, covmod=covmod, anisotropy=(0.2, 0.5, 1, 45, 25)) # variogram range will be smaller in z and y directions (and then rotated)\n",
    "field = fftma_sim.simnew()\n",
    "\n",
    "# and plot again\n",
    "xyz = np.mgrid[[slice(0,n,1) for n in domainsize]].reshape(3,-1).T\n",
    "plt.figure()\n",
    "axfield = plt.subplot2grid((1,1), (0,0), projection='3d', rowspan=1)\n",
    "axfield.scatter(xyz[:,0],xyz[:,1],xyz[:,2], c=field.flatten(), cmap='jet', vmin=-3.5, vmax=3.5, marker='o',linewidth=0)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Equality constraints\n",
    "Now we want to simulate spatial random fields conditioned on a covariance model and on linear (equality '=') observations. Therefore, we can use RMWS. As RMWS works in standard normal space we (for now) assume that these obervations are sampled from a standard normal distribution. Further, we assume that we know the covariance model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define field size\n",
    "domainsize = (100, 100)\n",
    "# define how many conditional realizations we want to simulate\n",
    "nfields = 20 \n",
    "# spatial covariance which we assume to know from our observations\n",
    "# the covariance model is defined as: Sill Model(Range) \n",
    "# and it can be a linear combination of multiple models\n",
    "# available covariance models are: Nug (Nugget), Exp (Exponential), \n",
    "# Sph (Spherical), Gau (Gaussian)\n",
    "covariancemodel = '0.01 Nug(0.0) + 0.99 Exp(7.2)' \n",
    "\n",
    "# define equality observations\n",
    "# observation point coordinates\n",
    "cp = np.array([[23, 37], [26, 34], [11, 89], [69, 51], [88, 39]]) \n",
    "# observation values in standard normal space \n",
    "cv = np.array([-1.1, -1.8, 2.8, 0.2, 1.9]) \n",
    "\n",
    "# if only conditional simulation is required\n",
    "# you don't need to define a nonlineartemplate\n",
    "\n",
    "# use random mixing for the conditional simulation\n",
    "# initialize RMWS\n",
    "rm = RMWS(domainsize = domainsize,\n",
    "          covmod = covariancemodel,\n",
    "          nFields = nfields,\n",
    "          cp = cp,\n",
    "          cv = cv)\n",
    "\n",
    "# run simulations\n",
    "rm()\n",
    "\n",
    "# these are the simulated conditional fields\n",
    "condfields = rm.finalFields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's have a look at four randomly selected fields \n",
    "rn = np.random.randint(0,nfields,4)\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "for i in range(4):\n",
    "    plt.subplot(2, 2, i+1)\n",
    "    ax = plt.gca()    \n",
    "    im = ax.imshow(condfields[rn[i]],\n",
    "               origin='lower', interpolation='nearest', \n",
    "               cmap='jet', vmin=-4, vmax=4)\n",
    "    ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "    divider = make_axes_locatable(ax)\n",
    "    cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "    plt.colorbar(im, cax=cax)\n",
    "plt.show()\n",
    "\n",
    "# plot the mean field (average over all the 20 simulated fields,\n",
    "# if one simulates enough fields the mean will \n",
    "# be similar to an interpolation)\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.subplot(1, 2, 1)\n",
    "ax = plt.gca()\n",
    "im = ax.imshow(np.mean(condfields,axis=0),\n",
    "           origin='lower', interpolation='nearest', \n",
    "           cmap='jet', vmin=-4, vmax=4)\n",
    "ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "plt.title('mean field')\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "\n",
    "# plot the variance field (variance over all 20)\n",
    "plt.subplot(1, 2, 2)\n",
    "ax = plt.gca()\n",
    "im = ax.imshow(np.var(condfields,axis=0),\n",
    "           origin='lower', interpolation='nearest', \n",
    "           cmap='jet')\n",
    "ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "plt.title('variance field')\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# but one can of course also check the actual simulated values \n",
    "# at the cp locations across these four fields;\n",
    "# NOTE that minor deviations from the actual values are due \n",
    "# to numerical rounding issues\n",
    "for i, c in enumerate(cp):\n",
    "    print('should be:', cv[i])\n",
    "    print('is:', condfields[rn, c[0], c[1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inequality constraints\n",
    "Sweet, that seems to work. Now let's try to additionally include some inequality (less/equal '<=') point constraints. These can, for example, arise from measurements below detection limits. Apart from the inequalities we keep the same setup as in the first example. Note that RMWS can also handle greater equal ('>=') constraints (keywords: ge_cp and ge_cv)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define inequality (<=) observations\n",
    "# <= observation coordinates\n",
    "cp_leq = np.array([[29,39], [81,46]]) \n",
    "# <= observation values (upper bounds) in standard normal space\n",
    "cv_leq = np.array([-2.1, 2.2])  \n",
    "    \n",
    "# add the additional variable for the inequalities\n",
    "rm = RMWS(domainsize = domainsize,\n",
    "          covmod = covariancemodel,\n",
    "          nFields = nfields,\n",
    "          cp = cp,\n",
    "          cv = cv,\n",
    "          le_cp = cp_leq,\n",
    "          le_cv = cv_leq)\n",
    "# run simulations\n",
    "rm()\n",
    "\n",
    "# these are the conditional fields\n",
    "condfields = rm.finalFields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's again have a look at four fields\n",
    "rn = np.random.randint(0,nfields,4)\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "for i in range(4):\n",
    "    plt.subplot(2, 2, i+1)\n",
    "    ax = plt.gca()    \n",
    "    im = ax.imshow(condfields[rn[i]],\n",
    "               origin='lower', interpolation='nearest', \n",
    "               cmap='jet', vmin=-4, vmax=4)\n",
    "    ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "    ax.plot(cp_leq[:,1], cp_leq[:,0], '.', c='black', mew=2, ms=8)\n",
    "    divider = make_axes_locatable(ax)\n",
    "    cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "    plt.colorbar(im, cax=cax)\n",
    "plt.show()\n",
    "\n",
    "# plot the mean field \n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.subplot(1,2,1)\n",
    "ax = plt.gca()\n",
    "im = ax.imshow(np.mean(condfields,axis=0),\n",
    "           origin='lower', interpolation='nearest', \n",
    "           cmap='jet', vmin=-4, vmax=4)\n",
    "ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "ax.plot(cp_leq[:,1], cp_leq[:,0], '.', c='black', mew=2, ms=8)\n",
    "plt.title('mean field')\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "\n",
    "# plot the variance field\n",
    "plt.subplot(1,2,2)\n",
    "ax = plt.gca()\n",
    "im = ax.imshow(np.var(condfields,axis=0),\n",
    "           origin='lower', interpolation='nearest', \n",
    "           cmap='jet')\n",
    "ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "ax.plot(cp_leq[:,1], cp_leq[:,0], '.', c='black', mew=2, ms=8)\n",
    "plt.title('variance field')\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the conditioning point values again\n",
    "for i, c in enumerate(cp):\n",
    "    print('should be:', cv[i])\n",
    "    print('is:', condfields[rn, c[0], c[1]])\n",
    "    \n",
    "# also check the values at the <= constraints\n",
    "print('\\n')\n",
    "for i, c in enumerate(cp_leq):\n",
    "    print('should be less or equal to:', cv_leq[i])\n",
    "    print('is:', condfields[rn, c[0], c[1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulations in actual data space\n",
    "Until now we have assumed that all our observations are in standard normal. However, this is rather seldom the case in real world data. Now we assume that we have observed data in a different space and that we have 'fitted' a marginal distribution to this data. For simplicity we don't actually fit a distribution here but we assume that our observations follow an exponential distribution with $\\lambda=1$. As RMWS works in standard normal space though we need to transform these observations from their original data space into standard normal space. After the simulation we then back-transform from standard normal to data space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define field size\n",
    "domainsize = (100, 100)\n",
    "# how many conditional realizations do we want to simulate\n",
    "nfields = 20 \n",
    "# spatial covariance which we assume to know from our observations\n",
    "covariancemodel = '0.01 Nug(0.0) + 0.99 Exp(17.2)' \n",
    "\n",
    "# define equality observations\n",
    "# observation point coordinates\n",
    "cp = np.array([[23,37], [26, 34], [11,89], [69,51], [88,39]]) \n",
    "# observation values in actual data space\n",
    "cv = np.array([0.2, 0.9, 7.1, 5.5, 2.3]) \n",
    "# transform these values into standard normal space for the simulation\n",
    "# we assume that we have 'fitted' a marginal distribution to our observations\n",
    "# here it is an exponential distribution with lambda=1\n",
    "cvn = st.norm.ppf(st.expon.cdf(cv))\n",
    "    \n",
    "# use random mixing for the conditional simulation\n",
    "# initialize RMWS\n",
    "rm = RMWS(domainsize = domainsize,\n",
    "          covmod = covariancemodel,\n",
    "          nFields = nfields,\n",
    "          cp = cp,\n",
    "          cv = cvn)\n",
    "\n",
    "# run simulations\n",
    "rm()\n",
    "\n",
    "# these are the simulated conditional fields (in standard normal)\n",
    "condfields = rm.finalFields\n",
    "# thus we need to back-transform to their original data space \n",
    "# using the inverse distribution function\n",
    "condfields_ds = st.expon.ppf(st.norm.cdf(rm.finalFields))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now we can have a look at some of the fields in standard normal space \n",
    "# as well as in data space;\n",
    "# in order to plot the same ones we determine the random numbers first\n",
    "rn = np.random.randint(0,nfields,4)\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "for i in range(4):\n",
    "    plt.subplot(2, 2, i+1)\n",
    "    ax = plt.gca()    \n",
    "    im = ax.imshow(condfields[rn[i]],\n",
    "               origin='lower', interpolation='nearest', \n",
    "               cmap='jet', vmin=-4, vmax=4)\n",
    "    ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "    ax.plot(cp_leq[:,1], cp_leq[:,0], '.', c='black', mew=2, ms=8)\n",
    "    divider = make_axes_locatable(ax)\n",
    "    cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "    plt.colorbar(im, cax=cax)\n",
    "plt.show()\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "for i in range(4):\n",
    "    plt.subplot(2, 2, i+1)\n",
    "    ax = plt.gca()    \n",
    "    im = ax.imshow(condfields_ds[rn[i]],\n",
    "               origin='lower', interpolation='nearest', \n",
    "               cmap='jet', vmin=0, vmax=9)\n",
    "    ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "    ax.plot(cp_leq[:,1], cp_leq[:,0], '.', c='black', mew=2, ms=8)\n",
    "    divider = make_axes_locatable(ax)\n",
    "    cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "    plt.colorbar(im, cax=cax)\n",
    "plt.show()\n",
    "\n",
    "# plot the mean field \n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.subplot(1,2,1)\n",
    "ax = plt.gca()\n",
    "im = ax.imshow(np.mean(condfields,axis=0),\n",
    "           origin='lower', interpolation='nearest', \n",
    "           cmap='jet', vmin=-4, vmax=4)\n",
    "ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "ax.plot(cp_leq[:,1], cp_leq[:,0], '.', c='black', mew=2, ms=8)\n",
    "plt.title('mean field')\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "\n",
    "# plot the variance field\n",
    "plt.subplot(1,2,2)\n",
    "ax = plt.gca()\n",
    "im = ax.imshow(np.var(condfields,axis=0),\n",
    "           origin='lower', interpolation='nearest', \n",
    "           cmap='jet')\n",
    "ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "ax.plot(cp_leq[:,1], cp_leq[:,0], '.', c='black', mew=2, ms=8)\n",
    "plt.title('variance field')\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "plt.show()\n",
    "\n",
    "# plot the mean field in data space\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.subplot(1,2,1)\n",
    "ax = plt.gca()\n",
    "im = ax.imshow(np.mean(condfields_ds,axis=0),\n",
    "           origin='lower', interpolation='nearest', \n",
    "           cmap='jet', vmin=0, vmax=9)\n",
    "ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "ax.plot(cp_leq[:,1], cp_leq[:,0], '.', c='black', mew=2, ms=8)\n",
    "plt.title('mean field')\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "\n",
    "# plot the variance field in data space\n",
    "plt.subplot(1,2,2)\n",
    "ax = plt.gca()\n",
    "im = ax.imshow(np.var(condfields_ds,axis=0),\n",
    "           origin='lower', interpolation='nearest', \n",
    "           cmap='jet')\n",
    "ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "ax.plot(cp_leq[:,1], cp_leq[:,0], '.', c='black', mew=2, ms=8)\n",
    "plt.title('variance field')\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the conditioning point values again but this time in data space\n",
    "for i, c in enumerate(cp):\n",
    "    print('should be:', cv[i])\n",
    "    print('is:', condfields_ds[rn, c[0], c[1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Equality constraints and anisotropy\n",
    "Let's go back to the basic case where we want to simulate spatial random fields conditioned on a covariance model and on linear (equality '=') observations. As RMWS works in standard normal space we (for now) assume that these obervations are sampled from a standard normal distribution. Further, we assume that we know the covariance model and this time the covariance exhibits anisotropy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define field size\n",
    "domainsize = (150, 100)\n",
    "# define how many conditional realizations we want to simulate\n",
    "nfields = 20 \n",
    "# spatial covariance which we assume to know from our observations\n",
    "# the covariance model is defined as: Sill Model(Range) \n",
    "# and it can be a linear combination of multiple models\n",
    "# available covariance models are: Nug (Nugget), Exp (Exponential), \n",
    "# Sph (Spherical), Gau (Gaussian)\n",
    "covariancemodel = '0.01 Nug(0.0) + 0.99 Exp(15.2)' \n",
    "\n",
    "# define equality observations\n",
    "# observation point coordinates\n",
    "cp = np.array([[23, 37], [26, 34], [11, 89], [69, 51], [88, 39]]) \n",
    "# observation values in standard normal space \n",
    "cv = np.array([-1.1, -1.8, 2.8, 0.2, 1.9]) \n",
    "\n",
    "# define the anisotropy\n",
    "anisotropy = (0.2, 1, 35) \n",
    "\n",
    "# use random mixing for the conditional simulation\n",
    "# initialize RMWS\n",
    "rm = RMWS(domainsize = domainsize,\n",
    "          covmod = covariancemodel,\n",
    "          nFields = nfields,\n",
    "          cp = cp,\n",
    "          cv = cv,\n",
    "          anisotropy = anisotropy)\n",
    "\n",
    "# run simulations\n",
    "rm()\n",
    "\n",
    "# these are the simulated conditional fields\n",
    "condfields = rm.finalFields"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's have a look at four randomly selected fields \n",
    "rn = np.random.randint(0,nfields,4)\n",
    "\n",
    "plt.figure(figsize=(10, 10))\n",
    "for i in range(4):\n",
    "    plt.subplot(2, 2, i+1)\n",
    "    ax = plt.gca()    \n",
    "    im = ax.imshow(condfields[rn[i]],\n",
    "               origin='lower', interpolation='nearest', \n",
    "               cmap='jet', vmin=-4, vmax=4)\n",
    "    ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "    divider = make_axes_locatable(ax)\n",
    "    cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "    plt.colorbar(im, cax=cax)\n",
    "plt.show()\n",
    "\n",
    "# plot the mean field (average over all the 20 simulated fields,\n",
    "# if one simulates enough fields the mean will \n",
    "# be similar to an interpolation)\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.subplot(1, 2, 1)\n",
    "ax = plt.gca()\n",
    "im = ax.imshow(np.mean(condfields,axis=0),\n",
    "           origin='lower', interpolation='nearest', \n",
    "           cmap='jet', vmin=-4, vmax=4)\n",
    "ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "plt.title('mean field')\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "\n",
    "# plot the variance field (variance over all 20)\n",
    "plt.subplot(1, 2, 2)\n",
    "ax = plt.gca()\n",
    "im = ax.imshow(np.var(condfields,axis=0),\n",
    "           origin='lower', interpolation='nearest', \n",
    "           cmap='jet')\n",
    "ax.plot(cp[:,1], cp[:,0], 'x', c='black', mew=2, ms=8)\n",
    "plt.title('variance field')\n",
    "divider = make_axes_locatable(ax)\n",
    "cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "plt.colorbar(im, cax=cax)\n",
    "plt.show()"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}