{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Simulating a Control-Bounded ADC\n\nThis example shows how to simulate the interactions between an analog system\nand a digital control while the former is excited by an analog signal.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport cbadc\nimport numpy as np"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## The Analog System\n\n<img src=\"file://images/chainOfIntegratorsGeneral.svg\" width=\"500\" align=\"center\" alt=\"The chain of integrators ADC.\">\n\nFirst, we have to decide on an analog system. For this tutorial, we will\ncommit to a chain-of-integrators ADC,\nsee :py:class:`cbadc.analog_system.ChainOfIntegrators`, as our analog\nsystem.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# We fix the number of analog states.\nN = 6\n# Set the amplification factor.\nbeta = 6250.\nrho = - 1e-2\nkappa = - 1.0\n# In this example, each nodes amplification and local feedback will be set\n# identically.\nbetaVec = beta * np.ones(N)\nrhoVec = betaVec * rho\nkappaVec = kappa * beta * np.eye(N)\n\n# Instantiate a chain-of-integrators analog system.\nanalog_system = cbadc.analog_system.ChainOfIntegrators(betaVec, rhoVec, kappaVec)\n# print the analog system such that we can very it being correctly initalized.\nprint(analog_system)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## The Digital Control\n\nIn addition to the analog system, our simulation will require us to specify a\ndigital control. For this tutorial, we will use\n:py:class:`cbadc.digital_control.DigitalControl`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Set the time period which determines how often the digital control updates.\nT = 1.0/(2 * beta)\n# Set the number of digital controls to be same as analog states.\nM = N\n# Initialize the digital control.\ndigital_control = cbadc.digital_control.DigitalControl(T, M)\n# print the digital control to verify proper initialization.\nprint(digital_control)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## The Analog Signal\n\nThe final and third component of the simulation is an analog signal.\nFor this tutorial, we will choose a\n:py:class:`cbadc.analog_signal.Sinusodial`. Again, this is one of several\npossible choices.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Set the peak amplitude.\namplitude = 0.5\n# Choose the sinusodial frequency via an oversampling ratio (OSR).\nOSR = 1 << 9\nfrequency = 1.0 / (T * OSR)\n\n# We also specify a phase an offset these are hovewer optional.\nphase = np.pi / 3\noffset = 0.0\n\n# Instantiate the analog signal\nanalog_signal = cbadc.analog_signal.Sinusodial(amplitude, frequency, phase, offset)\n# print to ensure correct parametrization.\nprint(analog_signal)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Simulating\n\nNext, we set up the simulator. Here we use the\n:py:class:`cbadc.simulator.StateSpaceSimulator` for simulating the\ninvolved differential equations as outlined in\n:py:class:`cbadc.analog_system.AnalogSystem`.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Simulate for 2^18 control cycles.\nend_time = T * (1 << 18)\n\n# Instantiate the simulator.\nsimulator = cbadc.simulator.StateSpaceSimulator(analog_system, digital_control, [\n                                analog_signal], t_stop=end_time)\n# Depending on your analog system the step above might take some time to\n# compute as it involves precomputing solutions to initial value problems.\n\n# Let's print the first 20 control decisions.\nindex = 0\nfor s in simulator:\n    if (index > 19):\n        break\n    print(f\"step:{index} -> s:{np.array(s)}\")\n    index += 1\n\n# To verify the simulation parametrization we can\nprint(simulator)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Tracking the Analog State Vector\n\nClearly the output type of the generator simulator above is the sequence of\ncontrol signals s[k]. Sometimes we are interested in also monitoring the\ninternal states of analog system during simulation.\n\nTo this end we use the\n:func:`cbadc.simulator.StateSpaceSimulator.state_vector` and an\n:func:`cbadc.simulator.extended_simulation_result`.\n\nNote that the :func:`cbadc.simulator.extended_simulation_result` is\ndefined like this\n\n.. code-block:: python\n\n  def extended_simulation_result(simulator):\n      for control_signal in simulator:\n          analog_state = simulator.state_vector()\n          yield {\n              'control_signal': np.array(control_signal),\n              'analog_state': np.array(analog_state)\n          }\n\nSo, in essence, we are creating a new generator from the old with an extended\noutput.\n\n<div class=\"alert alert-info\"><h4>Note</h4><p>The convenience function extended_simulation_result is one of many\n          such convenience functions found in the\n          :py:mod:`cbadc.simulator` module.</p></div>\n\nWe can achieve this by appending yet another generator to the control signal\nstream as:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Repeating the steps above we now get for the following\n# ten control cycles.\n\next_simulator = cbadc.simulator.extended_simulation_result(simulator)\nfor res in ext_simulator:\n    if (index > 29):\n        break\n    print(\n        f\"step:{index} -> s:{res['control_signal']}, x:{res['analog_state']}\")\n    index += 1"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n### Saving to File\n\nIn general, simulating the analog system and digital control interaction\nis a computationally much more intense procedure compared to the digital\nestimation step. This is one reason, and there are more, why\nyou would want to store the intermediate control signal sequence to a file.\n\nFor this purpose use the\n:func:`cbadc.utilities.control_signal_2_byte_stream` and\n:func:`cbadc.utilities.write_byte_stream_to_file` functions.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Instantiate a new simulator and control.\nsimulator = cbadc.simulator.StateSpaceSimulator(analog_system, digital_control, [\n                                analog_signal], t_stop=end_time)\n\n# Construct byte stream.\nbyte_stream = cbadc.utilities.control_signal_2_byte_stream(simulator, M)\n\n\ndef print_next_10_bytes(stream):\n    global index\n    for byte in stream:\n        if (index < 40):\n            print(f\"{index} -> {byte}\")\n            index += 1\n        yield byte\n\n\ncbadc.utilities.write_byte_stream_to_file(\"sinusodial_simulation.adcs\",\n                          print_next_10_bytes(byte_stream))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Evaluating the Analog State Vector in Between Control Signal Samples\n\nIf we wish to simulate the analog state vector trajectory between\ncontrol updates, this can be achieved using the Ts parameter of the\n:py:class:`cbadc.simulator.StateSpaceSimulator`. Technically you can scale\n$T_s = T / \\alpha$ for any positive number $\\alpha$. For such a\nscaling, the simulator will generate $\\alpha$ more control signals per\nunit of time. However, digital control is still restricted to only update\nthe control signals at multiples of $T$.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Set sampling time three orders of magnitude smaller than the control period\nTs = T / 1000.0\n\n# Simulate for 10000 control cycles.\nsize = 15000\nend_time = size * Ts\n\n# Initialize a new digital control.\nnew_digital_control = cbadc.digital_control.DigitalControl(T, M)\n\n# Instantiate a new simulator with a sampling time.\nsimulator = cbadc.simulator.StateSpaceSimulator(analog_system, new_digital_control, [\n                                analog_signal], t_stop=end_time, Ts=Ts)\n\n# Create data containers to hold the resulting data.\ntime_vector = np.arange(size) * Ts / T\nstates = np.zeros((N, size))\ncontrol_signals = np.zeros((M, size), dtype=np.int8)\n\n# Iterate through and store states and control_signals.\nfor index, res in enumerate(cbadc.simulator.extended_simulation_result(simulator)):\n    states[:, index] = res['analog_state']\n    control_signals[:, index] = res['control_signal']\n\n# Plot all analog state evolutions.\nplt.figure()\nplt.title(\"Analog state vectors\")\nfor index in range(N):\n    plt.plot(time_vector, states[index, :], label=f\"$x_{index + 1}(t)$\")\nplt.grid(b=True, which='major', color='gray', alpha=0.6, lw=1.5)\nplt.xlabel('$t/T$')\nplt.xlim((0, 10))\nplt.legend()\n\n# reset figure size and plot individual results.\nplt.rcParams['figure.figsize'] = [6.40, 6.40 * 2]\nfig, ax = plt.subplots(N, 2)\nfor index in range(N):\n    color = next(ax[0, 0]._get_lines.prop_cycler)['color']\n    ax[index, 0].grid(b=True, which='major', color='gray', alpha=0.6, lw=1.5)\n    ax[index, 1].grid(b=True, which='major', color='gray', alpha=0.6, lw=1.5)\n    ax[index, 0].plot(time_vector, states[index, :], color=color)\n    ax[index, 1].plot(time_vector, control_signals[index, :],\n                      '--', color=color)\n    ax[index, 0].set_ylabel(f\"$x_{index + 1}(t)$\")\n    ax[index, 1].set_ylabel(f\"$s_{index + 1}(t)$\")\n    ax[index, 0].set_xlim((0, 15))\n    ax[index, 1].set_xlim((0, 15))\n    ax[index, 0].set_ylim((-1, 1))\nfig.suptitle(\"Analog state and control contribution evolution\")\nax[-1, 0].set_xlabel(\"$t / T$\")\nax[-1, 1].set_xlabel(\"$t / T$\")\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Analog State Statistics\n\nAs in the previous section, visualizing the analog state trajectory is a\ngood way of identifying problems and possible errors. Another way of making\nsure that the analog states remain bounded is to estimate their\ncorresponding densities (assuming i.i.d samples).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Compute L_2 norm of analog state vector.\nL_2_norm = np.linalg.norm(states, ord=2, axis=0)\n# Similarly, compute L_infty (largest absolute value) of the analog state\n# vector.\nL_infty_norm = np.linalg.norm(states, ord=np.inf, axis=0)\n\n# Estimate and plot densities using matplotlib tools.\nbins = 150\nplt.rcParams['figure.figsize'] = [6.40, 4.80]\nfig, ax = plt.subplots(2, sharex=True)\nax[0].grid(b=True, which='major', color='gray', alpha=0.6, lw=1.5)\nax[1].grid(b=True, which='major', color='gray', alpha=0.6, lw=1.5)\nax[0].hist(L_2_norm, bins=bins, density=True)\nax[1].hist(L_infty_norm, bins=bins, density=True, color=\"orange\")\nplt.suptitle(\"Estimated probability densities\")\nax[0].set_xlabel(\"$\\|\\mathbf{x}(t)\\|_2$\")\nax[1].set_xlabel(\"$\\|\\mathbf{x}(t)\\|_\\infty$\")\nax[0].set_ylabel(\"$p ( \\| \\mathbf{x}(t) \\|_2 ) $\")\nax[1].set_ylabel(\"$p ( \\| \\mathbf{x}(t) \\|_\\infty )$\")\nfig.tight_layout()"
      ]
    }
  ],
  "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.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}