<jupyter_start><jupyter_text>NgspiceX - Monte Carlo simulation example1.0 | 22.08.2021 | jchabloz | First draft1.1 | 21.12.2021 | jchabloz | New run with NgSpice 35 Table of ContentsNgspiceX - Monte Carlo simulation exampleIntroductionSimulation preparationExpected resultsMonte Carlo simulation loopSimulation resultsComparison with expected resultsConclusion  IntroductionThis notebook intends to demonstrate how to use [ngspice-connect](https://github.com/jchabloz/ngspice-connect) to run efficiently a [Monte Carlo simulation](https://www.wikiwand.com/en/Monte_Carlo_method) with [ngspice](http://ngspice.sourceforge.net/) by using Python instead of resorting to the simulator internal scripting language.In this demo, we create a simple resistive divider circuit with the value of the resistances including some statistical variations. We distinguish a *lot (or process) variation*, that would affect the nominal value of the resistances in a similar manner for all the resistances in the circuit from a *device variation* that affects separately each resistor and that models the mismatch between them.  Simulation preparationThe demo start by importing all the required packages and creating an instance of the ``NgSpice`` class to access the ngspice shared library.<jupyter_code>from ngspicex import NgSpice
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from math import sqrt
from tqdm import tqdm
from time import perf_counter
ngx = NgSpice()<jupyter_output>******
 ** ngspice-35 shared library
 ** Creation Date: Tue Dec 21 19:51:18 UTC 2021
 ******<jupyter_text>We define a few parameters for the simulation:<jupyter_code>V1 = 1.0            #[V]    Input voltage
sigma_R_lot = 100   #[Ohms] Standard deviation (lot)
sigma_R_dev = 200   #[Ohms] Standard deviation (device)
R = 1e4             #[Ohms] Resistances nominal value
lot = 0             #[-]    0: don't use lot variations, 1: use lot variations
dev = 1             #[-]    0: don't use device variations, 1: use device variations
niter = 500         #[-]    Number of iterations<jupyter_output><empty_output><jupyter_text>The next step is to create the circuit in ngspice. If we want to use parameters, the easiest is to use the ``send_circ()`` function, especially for a small netlist. The netlist is sent line after line to ngspice.The statistical variations of the resistors value are introduced via parameters and the built-in ngspice `agauss` function.<jupyter_code>ngx.send_circ(
    "*Test Monte-Carlo of resistive divider",
    ".param __LOT__ = {:d}".format(lot),
    ".param __DEV__ = {:d}".format(dev),
    ".param r0 = {}".format(R),
    ".param dr0 = agauss(0.0, '{}*__LOT__', 1) rval='r0 + dr0'".format(sigma_R_lot),
    ".param dr1 = agauss(0.0, '{}*__DEV__', 1)".format(sigma_R_dev),
    "R1 n1 n2 'rval + dr1'",
    "R2 n2 0 'rval + dr1'",
    "V1 n1 0 DC {}".format(V1),
    ".op",
    ".end")<jupyter_output>Source Deck
 No compatibility mode selected!
 Circuit: *test monte-carlo of resistive divider<jupyter_text>Expected resultsLet us consider the case for ``lot = 0, dev = 1``.What are the expected statistics for the voltage at the middle point of the resistive divider and the current through it?Since we consider the case for which both resistances of the divider are nominally equal $R_1 = R_2 = R$, we can calculate the expected standard deviations for both the output voltage $\sigma_{V_2}$ and the current $\sigma_I$ as$$\begin{align}\sigma_{V_2} &= \frac{1}{\sqrt{2}}\frac{\sigma_R}{R}\cdot\frac{V_1}{2}\\\sigma_{I} &= \frac{1}{\sqrt{2}}\frac{\sigma_R}{R}\cdot\frac{V_1}{2R}\\\end{align}$$.<jupyter_code>sigma_V2 = 1/sqrt(2)*sigma_R_dev/R*V1/2
sigma_I = 1/sqrt(2)*sigma_R_dev/R*V1/2/R

print("Expected standard deviation on V(n2) (lot=0, dev=1) = {:.3f}mV".format(sigma_V2*1e3))
print("Expected standard deviation on I(V1) (lot=0, dev=1) = {:.3f}uA".format(sigma_I*1e6))<jupyter_output>Expected standard deviation on V(n2) (lot=0, dev=1) = 7.071mV
Expected standard deviation on I(V1) (lot=0, dev=1) = 0.707uA<jupyter_text>Monte Carlo simulation loopUsing Python, it is actually quite easy to run a simulation loop, collect the results and concatenate them into a single [pandas](https://pandas.pydata.org) DataFrame. Using [tqdm](https://tqdm.github.io/) to show a progress bar makes for a nice way to monitor the progress of the loop.Note that in order to make sure that the values of the parameters are recomputed for each iteration, the `reset` commands needs to be used.<jupyter_code>dfx=[]
t0 = perf_counter()
for i in tqdm(range(niter)):
    ngx.reset()
    ngx.run(silent=True)
    dfx.append(ngx.get_all_vectors())
df_out = pd.concat(dfx)
t1 = perf_counter()
print("Monte Carlo loop time: {:.3f}s".format(t1 - t0))<jupyter_output>100%|████████████████████████████████████████| 500/500 [00:01<00:00, 461.65it/s]<jupyter_text>Simulation resultsHaving directly the simulation results in a pandas DataFrame allows to use all the available powerful tools of pandas and associated libraries to analyze them.For example, using the `describe()` method allows to get a quick overview of the content and statistics of a DataFrame.<jupyter_code>df_out.describe().T<jupyter_output><empty_output><jupyter_text>We can also plot the distributions for the variables of interest. For example using the nice [seaborn library](https://seaborn.pydata.org/).<jupyter_code>sns.set_theme(context='notebook', style='whitegrid')
f, axs = plt.subplots(1, 2, figsize=(14,6), sharey=True)
sns.histplot(data=df_out, x="i(v1)", ax=axs[0], stat='probability', kde=True)
_ = sns.histplot(data=df_out, x="n2", ax=axs[1], stat='probability', kde=True)<jupyter_output><empty_output><jupyter_text>Comparison with expected results<jupyter_code>sigma_V2_sim = df_out.n2.std()
print("Expected standard deviation on V(n2) = {:.3f}mV".format(sigma_V2*1e3))
print("Simulated standard deviation on V(n2) = {:.3f}mV".format(sigma_V2_sim*1e3))
print("Error = {:.1f}%".format((sigma_V2_sim/sigma_V2 - 1)*100.0))
sigma_I_sim = df_out['i(v1)'].std()
print("Expected standard deviation on I(V1) = {:.3f}uA".format(sigma_I*1e6))
print("Simulated standard deviation on I(V1) = {:.3f}uA".format(sigma_I_sim*1e6))
print("Error = {:.1f}%".format((sigma_I_sim/sigma_I - 1)*100.0))<jupyter_output>Expected standard deviation on I(V1) = 0.707uA
Simulated standard deviation on I(V1) = 0.736uA
Error = 4.1%
{'fraction_non_alphanumeric': 0.12648716343143393, 'fraction_numerical': 0.03318722604884158, 'mean_word_length': 5.2392578125, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 28, 'https://': 5, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '303363', 'license': 'MIT', 'n_tokens_mistral': 2275, 'n_tokens_neox': 2120, 'n_words': 670, 'path': 'examples/montecarlo_resdiv/ngspicex_montecarlo_resdiv.ipynb', 'repo_name': 'jchabloz/ngspice-connect', 'size': '6388'}