{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# SOL Calibration in Python\n\nThis example  a performs a short,  open, and load (SOL) calibration with a linear\nsensitivity analysis.\n\nIn this scenario, we have some raw S-parameter measurements stored in\ntab separated text files of our calibration standards\nand of our DUT - which is a load. We also have some\ndata-defined models of our calibration standards stored in the HDF5 format. These models contain some\nlinear uncertainty mechanisms which carry some category information describing their origin.\nThis tutorial will walk through defining file reader functions, analysis functions,\nand performing the analysis with and without uncertainties.\n\nThe data for this example is stored in th 'sol_demo_data' folder of the source code repository hosted\non github. The data set is intended to act as an example use case of the software library only.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from rmellipse.uobjects import RMEMeas\nfrom rmellipse.propagators import RMEProp\nfrom pathlib import Path\nimport h5py\nimport xarray as xr\nimport numpy as np\nimport matplotlib.pyplot as plt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Demo Data\n\nThe sample data are store inside the `rmellipse.sol_demo` submodule, and the paths\nto the files can be imported as well.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# text files paths as Path objects\n# replace these with the correct data paths\n# after downloading the sample data\nshort_raw_path = sample_data_dir / 'Short.s1p'\nopen_raw_path = sample_data_dir / 'Open.s1p'\nload_raw_path = sample_data_dir / 'Load.s1p'\ndut_raw_path = sample_data_dir / 'Dut.s1p'\n\n# HDF5 file containing the standard models\ndefs_path = sample_data_dir / 'definitions.hdf5'"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Writing Functions\n\nThe first thing we need to do is write some functions. These include\nthe text file reader and analysis functions\nto calculate our error box and correct our raw data.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Text File Reader\n\nFirst we define a function that can read our text files into a DataArray. We assign\na dimension called Frequency (GHz), where we store the frequencies corresponding to each\nS-parameter. We also store our data as a complex-valued array.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def from_s1p(path) -> xr.DataArray:\n\tarr = np.loadtxt(path, float, delimiter='\\t')\n\t# define a coordinate set\n\tcoords = {'Frequency (GHz)': arr[:, 0]}\n\t# convert to 1 port complex data\n\tvalues = arr[:, 1] + 1.0j * arr[:, 2]\n\t# create an xarray data set\n\tout = xr.DataArray(values, coords=coords, dims=('Frequency (GHz)'))\n\treturn out"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lets confirm this works by  reading in our raw measurements\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "raw_dut = from_s1p(dut_raw_path)\nraw_short = from_s1p(short_raw_path)\nraw_open = from_s1p(open_raw_path)\nraw_load = from_s1p(load_raw_path)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Calculating the Error Box\n\nWe want two analysis functions. The first takes in our definitions and raw measurements, and\ncalculates the error-box terms of the 1-port calibration. Importantly, all the\nimport objects that might be ``RMEMeas`` objects are exposed as input arguments.\nThe ``**`` packs any keyword arguments ino a dictionary. This function matches any\ninputs by matching device_names, and assuming keywords are formatted\nby 'def_<device_name>' and 'meas_<device_name>' only. Also, note that the function\nsignature says that the function is expecting xr.DataArray objects. This is what will be\npassed in when the function is wrapped in a propagator.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def SOL_cal(**stds: xr.DataArray) -> xr.DataArray:\n\t# get list of standards and definitions, ordered to correspond\n\tdefs = []\n\tms = []\n\tfor k, v in stds.items():\n\t\tif 'def' in k:\n\t\t\tdef_key = k\n\t\t\tm_key = def_key.replace('def_', 'raw_')\n\t\t\tdefs.append(stds[def_key])\n\t\t\tms.append(stds[m_key])\n\tN = len(defs)\n\tfrequencies = ms[-1]['Frequency (GHz)']\n\n\t# output has the same shape as the inputs except an\n\t# additional dimensions to hold the error terms\n\toutput_shape = list(ms[0].shape) + [3]\n\toutput_dims = list(ms[0].dims) + ['errterm']\n\toutput_coords = dict(ms[0].coords)\n\toutput_coords.update({'errterm': ['e00', 'e11', 'delta']})\n\n\t# pre allocate output xarray\n\t# 3 output has 3\n\tresult = np.zeros(output_shape, complex)\n\tresult = xr.DataArray(result, dims=output_dims, coords=output_coords)\n\n\t# pre allocated temporary arrays\n\t# that will be used to solve the set of equations\n\tmshape = list(result.shape)[:-1] + [N, 3]\n\tM = np.zeros(mshape, complex)\n\n\tyshape = list(result.shape)[:-1] + [N, 1]\n\ty = np.zeros(yshape, complex)\n\n\t# I am going to work with the underlying numpy arrays\n\t# here because it is convenient for linear algebra\n\t# for each device add a row to the regressor matrix\n\tfor i in range(N):\n\t\tS11_meas = ms[i].values\n\t\tS11_def = defs[i].values\n\t\tM[..., i, 0] = 1\n\t\tM[..., i, 1] = S11_meas * S11_def\n\t\tM[..., i, 2] = -S11_def\n\t\ty[..., i, 0] = S11_meas\n\n\t# this transposes M along last 2 dimensions\n\tn_dims = len(M.shape)\n\ttranspose_dims = np.arange(n_dims)\n\ttranspose_dims[[n_dims - 1, n_dims - 2]] = transpose_dims[[n_dims - 2, n_dims - 1]]\n\tMt = np.transpose(M, transpose_dims)\n\n\t# do least squares\n\tcoeff = np.linalg.inv(Mt @ M) @ Mt @ y\n\n\t# reassign values to output\n\tresult.loc[..., 'e00'] = coeff[..., 0, 0]\n\tresult.loc[..., 'e11'] = coeff[..., 1, 0]\n\tresult.loc[..., 'delta'] = coeff[..., 2, 0]\n\n\treturn result"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lets use the function without any uncertainty propagation\nto verify it works. I will do this by reading in the definitions\nand grabbing just a view into the underlying nominal DataArray.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "with h5py.File(defs_path, 'r') as f:\n\tdef_short = RMEMeas.from_h5(f['Short']).nom\n\tdef_open = RMEMeas.from_h5(f['Open']).nom\n\tdef_load = RMEMeas.from_h5(f['Load']).nom\n\nflist = raw_short['Frequency (GHz)']\nerrbox = SOL_cal(\n\tdef_short=def_short,\n\tdef_open=def_open,\n\tdef_load=def_load,\n\traw_short=raw_short,\n\traw_open=raw_open,\n\traw_load=raw_load,\n)\n\nprint(errbox)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Correcting Raw Data\n\nFinally, lets write a function that takes in our error box\nand corrects a raw measurement.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def SOL_correct(errorbox: xr.DataArray, device: xr.DataArray) -> xr.DataArray:\n\tS11 = device\n\te00 = errorbox.sel(errterm='e00')\n\te11 = errorbox.sel(errterm='e11')\n\tdelta = errorbox.sel(errterm='delta')\n\n\tcorrected = (-e00 + S11) / (-delta + e11 * S11)\n\n\t# return the corrected result\n\tr = xr.zeros_like(device)\n\tr.loc[...] = corrected\n\treturn r"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "And lets use it to correct our raw device.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dut = SOL_correct(errbox, raw_dut)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can plot this to check that it matches our expectations. This DUT is a load\nso we expect a magnitude close to 0.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.close('all')\nfig, ax = plt.subplots(2, 1)\nfig.suptitle('DUT Corrected Data')\nax[0].plot(dut['Frequency (GHz)'], np.abs(dut))\nax[0].set_ylabel('Linear Magnitude')\n\nax[1].plot(dut['Frequency (GHz)'], np.angle(dut, deg=True))\nax[1].set_xlabel('Frequency (GHz)')\nax[1].set_ylabel('Phase (deg)')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Linear Uncertainty Analysis\nIn this step, we will add the ability to propagate uncertainties associated with\nthe calibration standard definitions.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Import and Wrap Functions\n\nThe first step is to define a propagator and wrap our analysis\nfunctions in it. We are going to turn on the sensitivity analysis only.\nWhen we wrap our functions in the propagator, any RMEMeas objects passed\nthrough the function will be run through the linear sensitivity analysis.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "myprop = RMEProp(sensitivity=True)\n\nSOL_cal = myprop.propagate(SOL_cal)\nSOL_correct = myprop.propagate(SOL_correct)\n\n# turn xarray's into uncertainty objects with RMEMeas\nraw_dut = RMEMeas.from_nom('DUT', raw_dut)\nraw_short = RMEMeas.from_nom('DUT', raw_short)\nraw_open = RMEMeas.from_nom('DUT', raw_open)\nraw_load = RMEMeas.from_nom('DUT', raw_load)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Lets also grab our definitions with the full uncertainty information\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "with h5py.File(defs_path, 'r') as f:\n\tdef_short = RMEMeas.from_h5(f['Short'])\n\tdef_open = RMEMeas.from_h5(f['Open'])\n\tdef_load = RMEMeas.from_h5(f['Load'])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Propagate Functions\n\nNow we can use our analysis functions that were wrapped\nin the propagator to do our analysis.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "errbox = SOL_cal(\n\tdef_short=def_short,\n\tdef_open=def_open,\n\tdef_load=def_load,\n\traw_short=raw_short,\n\traw_open=raw_open,\n\traw_load=raw_load,\n)\n\ndut = SOL_correct(errbox, raw_dut)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Plot Results\n\nLets plot our corrected device like last time, and\nlets instead inspect the uncertainty measurements.\nWe would like to look at the magnitude and phase,\nso lets define some functions to propagate those conversions.\n\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Nominal Values\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "@myprop.propagate\ndef calc_mag(arr):\n\tout = xr.zeros_like(arr, dtype=float)\n\tout.values = np.abs(arr)\n\treturn out\n\n\n@myprop.propagate\ndef calc_phase(arr):\n\tout = xr.zeros_like(arr, dtype=float)\n\tout.values = np.angle(arr, deg=True)\n\treturn out\n\n\nmag = calc_mag(dut)\nphase = calc_phase(dut)\n\n\nk = 2\nfig, ax = plt.subplots(2, 1)\nmag_lower = mag.uncbounds(k=k).cov\nmag_upper = mag.uncbounds(k=-k).cov\nphs_lower = phase.uncbounds(k=k, deg=True).cov\nphs_upper = phase.uncbounds(k=-k, deg=True).cov\n\nax[0].fill_between(\n\tdut.nom['Frequency (GHz)'],\n\ty1=mag_lower,\n\ty2=mag_upper,\n\tcolor='k',\n\talpha=0.5,\n\tlabel=f'k = {k} Uncertainty',\n)\n\nax[1].fill_between(\n\tdut.nom['Frequency (GHz)'],\n\ty1=phs_lower,\n\ty2=phs_upper,\n\talpha=0.5,\n\tcolor='k',\n\tlabel=f'k = {k} Uncertainty',\n)\n\nax[0].plot(mag.nom['Frequency (GHz)'], mag.nom, label='nominal')\nax[0].set_ylabel('Linear Magnitude')\n\nax[1].plot(phase.nom['Frequency (GHz)'], phase.nom, label='nominal')\nax[1].set_xlabel('Frequency (GHz)')\nax[1].set_ylabel('Phase (deg)')\n\nhandles, labels = ax[0].get_legend_handles_labels()\nfig.legend(handles, labels, loc='upper center', ncols=3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Uncertainties Only\n\nLets look at the k=2 uncertainties\nonly. Notice how the phase uncertainty increases\nas the magnitude approaches zero. This is natural, as\nphase is not well defined if your magnitude is on the origin.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(2, 1)\nax[0].plot(\n\tmag.nom['Frequency (GHz)'],\n\tmag.stdunc(k=k).cov,\n)\nax[0].set_ylabel(f'Lin Magn k={k} Uncertainty')\nax[1].plot(phase.nom['Frequency (GHz)'], phase.stdunc(k=k).cov, label='nominal')\nax[1].set_xlabel('Frequency (GHz)')\nax[1].set_ylabel(f'Phase k={k} Uncertainty (deg)')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Categorizing Uncertainties\n\nThe definitions we use have some string metadata attached\nto each uncertainty mechanisms. This metadata is stored in\n`dut.covcats`. Each linear uncertainty mechanism has been\ngiven an `\"Origin\"` designation we can use to group each\nmechanisms that belongs to a common source. By calling\n`RMEMeas` each mechanism that belongs to the same Origin\nis collected into a single larger, independent mechanism.\nBy doing this, we see our uncertainties are dominated\nby the VNA drift, cable bend variations, connection\nrepeatability, and the physical dimensions of our\nthe primary standards used to define the calibration\nkit.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "grouped_mag = mag.categorize_by('Origin')\ngrouped_phase = phase.categorize_by('Origin')\ntotal_mag_variance = grouped_mag.stdunc().cov ** 2\ntotal_phase_variance = grouped_phase.stdunc().cov ** 2\nmag_variance = []\nphase_variance = []\nfig, axs = plt.subplots(2, 1)\nfor um in grouped_mag.umech_id:\n\tmag_var_i = grouped_mag.usel(umech_id=[um]).stdunc().cov ** 2\n\tphase_var_i = grouped_phase.usel(umech_id=[um]).stdunc().cov ** 2\n\tmag_variance.append(mag_var_i / total_mag_variance * 100)\n\tphase_variance.append(phase_var_i / total_phase_variance * 100)\n\naxs[0].stackplot(dut.nom['Frequency (GHz)'], mag_variance, labels=grouped_mag.umech_id)\naxs[1].stackplot(\n\tdut.nom['Frequency (GHz)'], mag_variance, labels=grouped_phase.umech_id\n)\n\nhandles, labels = axs[0].get_legend_handles_labels()\nfig.legend(handles, labels, loc='upper center', ncols=3)\n\naxs[1].set_xlabel('Frequency (GHz)')\naxs[0].set_ylabel('Mag (% of Variation)')\naxs[1].set_ylabel('Phase (% of Variation)')"
      ]
    }
  ],
  "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
}