{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Basics of RME Propagation\n\nThis is a basic example demonstrating how to propagate RMEMeas objects\nthrough functions.\n\nThis example generates a hypothetical voltage and current measurement each\nwith a Monte Carlo distribution and single linear uncertainty mechanism.\n\nWe define a propagator and a propagation function to calculate power and\nperform the uncertainty analysis.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Importing packages\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import rmellipse as rme\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport xarray as xr"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Creating a propagator\n\nThe first step is to define a propagator. This is the object that inspects\nfunctions for uncertainty objects, propagating them according to its settings.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "myprop = rme.RMEProp(montecarlo_sims=1000, sensitivity=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Creating RMEMeas objects\n\nRMEMeas objects (Rocky Mountain Ellipse Measurements) carry\nalong with them 2 key attributes `.cov` and `RMEMeas.mc`, along attributes for\npropagating metadata - `RMEMeas.covcats` and `RMEMeas.covdofs`.\n\nThe`RMEMeas.cov`is an xarray.DataArray\nobject that stores both the nominal value of the measurement, and copies of the nominal\nvalue perturbed by 1 standard deviation for each uncertainty mechanism. The\nperturbed copies are stored along the first dimension of the`RMEMeas.cov`attribute,\nwhich is always called 'umech_id'. Along that dimension\nthe first index is always labeled 'nominal' and holds the nominal value, the remaining\nindexes store the names of the linear uncertainty mechanisms. These are labelled\nto be unique. Importantly, if two RMEMeas objects being propagated share an\nuncertainty mechanism with the same umech_id (i.e. they both have a perturbed copy of the nominal\nwith the same index label), then they are assume to be fully correlated.\n\nThe `RMEMeas.mc` attribute stores Monte Carlo data in an xarray.DataArray.\nIts first dimension is called `sample_id`; every row is a stochastic sample,\nwith integer labels starting at zero. The independent nominal value remains in\nthe first row of `RMEMeas.cov`.\n\nIn this case, we will create a voltage and current measurement, both drawing\nfrom a gaussian distribution. Because we are creating a measurement using\na float, the resulting`RMEMeas.cov`attribute will be shape (2,) (index 0 for the nominal\nvalue and index 1 representing the uncertainty associated with the Gaussian\nwe sampled from). The `RMEMeas.mc` attribute will have shape (1000,), with\nevery index representing a sample from the Gaussian distribution.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def make_measurement(d1_coords, std=0.1):\n    \"\"\"Make a sample measurement with length 4 coordinate set.\"\"\"\n    nom = xr.DataArray(\n        np.ones((4, 2)),\n        dims=('d1', 'd2'),\n        coords={'d1': d1_coords, 'd2': np.arange(2)},\n    )\n\n    meas = rme.RMEMeas.from_nom(name='meas', nom=nom)\n\n    meas.add_umech(\n        name='mymechanisms',\n        value=meas.nom + np.ones(meas.nom.shape) * std,\n        dof=np.inf,\n        category={'Type': 'B', 'Origin': 'Data Sheet'},\n        add_uid=True,\n    )\n\n    for i in range(1000):\n        meas.add_mc_sample(meas.nom + np.random.normal(loc=0, scale=std, size=None))\n    return meas\n\n\nV = make_measurement([0, 1, 2, 3])\nI = make_measurement([0, 1.1, 2, 2.9])\n\nprint(V)\nprint('            dimensions of V.cov :', V.cov.dims, V.cov.shape)\nprint('             dimensions of V.mc :', V.mc.dims, V.mc.shape)\nprint('V linear uncertainty mechanisms :', V.umech_id)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Propagating Through a Function\n\nWriting a function to be propagated requires a little knowledge about how\nthe propagator works, but in trade enables you to take advantage of vectorized\noperations for significantly faster computation.\n\nFirst we define a function and wrap it the myprop.propagator decorator. This\nis equivalent to calling power = myprop.propagate(power) after defining\nthe function power. This will tell the propagator to inspect the arguments of\npower when it is called and propagate any RMEMeas objects that are\npositional or keyword arguments.\n\nThe propagator will call the function twice. Once for the linear sensitivity\nanalysis where it passes in the`RMEMeas.cov`attribute, and once for the monte-carlo\nanalysis where it passes in the`RMEMeas.mc` attribute.\n\nIn the case of the sensitivity analysis, the propagator will align the\nlinear uncertainty mechanisms of each RMEMeas object so that they can be\nused in math operations. In this example this means that the shape of the v\nand i arguments will have shape (3,) each, with the first axis aligning the\nuncertainty mechanisms that existed in both v and i. Uncertainty mechanisms\nthat don't exist in a variable are filled with a copy of the nominal in the\nnewly aligned dimension.\n\nFor the Monte Carlo analysis, the propagator uses samples from the data stored\nin the `RMEMeas.mc` attribute. We defined 1000 Monte Carlo simulations in our\npropagator, so v and i will have shape (1000,) when passed through our power\nfunction.\n\nNote that while the propagator is design to encourage vectorization for\nefficiency, if that is not possible simply you can turn off vectorization\nwith the vectorize setting. See the example on non-vectorized propagation\nfor more details.\n\nThe propagator will name your RMEMeas object using the name of the function\nyou used to propagate it.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "@myprop.propagate\ndef power(v, i):\n    print(v.dims, v.shape, i.shape)\n    return v * i\n\n\n# when we call this, note that the power function is called twice, and the shapes\n# of the v and i arguments are changed from the original definitions for the\n# `RMEMeas.cov`attribute, now (3,) after the uncertainty mechanisms were aligned.\np = power(V, I)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Evaluating Uncertainty\n\nWe can calculate the standard uncertainty of a RMEMeas object easily!\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "covunc, mcunc = p.stdunc(k=1)\n\n\nprint('sensitivity analysis uncertainty:', covunc)\nprint('monte-carlo uncertainty:', mcunc)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can easily calculate the lower/upper uncertainty bounds for a given\nexpansion factor on the uncertainty, degrees of freedom based on the linear\nsensitivity analysis, and confidence intervals based on the linear sensitivity\nanalysis.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# k = 2 uncertainty bounds estimated by\n# linear uncertainty analysis and monte carlo analysis\ncov_k2, mc_k2 = p.uncbounds(k=2)\n\n# degrees of freedom\n# on the linear uncertaint analysis\ndof = p.dof()\n\n# 95% confidence interval\n# of the linear uncertainty analysis\nlower_ci, upper_ci = p.confint(0.95)"
      ]
    }
  ],
  "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.12.4"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}