{
  "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": [
        "from rmellipse.uobjects import RMEMeas\nfrom rmellipse.propagators import RMEProp\nimport matplotlib.pyplot as plt"
      ]
    },
    {
      "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": [
        "N = 1000\n\nmyprop = RMEProp(montecarlo_sims=N, 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 the monte-carlo data. It is also an xarray.DataArray,\nwhere the first dimension is always called umech_id, and the first\nindex of the dimension represents the nominal value. However, the index is\ninstead a numeric index, starting at zero and counting up. Indexes starting from\n1 and up of this dimension represent samples of the probability distribution of\nthe measurement. The RMEMeas will randomly sample from this distribution to\nperform the monte-carlo analysis.\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 (1001,) (index 0 for the\nnominal and index 1-1001 representing the samples from the gaussian\ndistribution)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "use_sample_mean = False\n\nV = RMEMeas.from_dist(\n\tname='voltage', nom=2, std=0.01, samples=N, dist='gaussian', use_sample_mean=False\n)\n\nI = RMEMeas.from_dist(\n\tname='current', nom=1.5, std=0.01, samples=N, dist='gaussian', use_sample_mean=False\n)\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 will pick random samples from the\ndata stored in the`RMEMeas.mc` attribute. We defined 1000 montecarlo sims in our propagator\nso v and i will have shape (1001,) when passed through our power function (\n1000 random samples and 1 copy of the nominal value).\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\tprint(v.dims, v.shape, i.shape)\n\treturn 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": [
        "lower, upper = p.uncbounds(k=2)\ndof = p.dof()\nlower, upper = p.confint(0.95)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can plot the montecarlo distributions of our results as well.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(1, 1)\nax.hist(V.mc[1:], label=V.name)\nax.hist(I.mc[1:], label=I.name)\nax.hist(p.mc[1:], label=p.name)\nax.legend(loc='best')\nax.set_xlabel('Value')\nax.set_ylabel('Quantity')"
      ]
    }
  ],
  "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.13.3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}