Repository: qclib/qclib (language: python)

Issue:
Testing



Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index a2f181b3..16dd9d4b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -130,3 +130,6 @@ dmypy.json
 
 # Pycharm Stuff
 .idea/
+
+# vscode
+*.json
diff --git a/example/isometry/Cnots and depth.ipynb b/example/isometry/Cnots and depth.ipynb
new file mode 100644
index 00000000..9a090584
--- /dev/null
+++ b/example/isometry/Cnots and depth.ipynb
@@ -0,0 +1,150 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "canadian-reporter",
+   "metadata": {},
+   "source": [
+    "# Qiskit config"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "stylish-integer",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit import IBMQ, QuantumCircuit, ClassicalRegister, QuantumRegister, execute, Aer, transpile"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "another-decision",
+   "metadata": {},
+   "source": [
+    "# Experiment"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "painted-clone",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\n",
+      "m = 2\n",
+      "\n",
+      "n = 3\n",
+      "\t knill\t\t CNOTs = 63\t depth = 117\n",
+      "\t ccd\t\t CNOTs = 24\t depth = 47\n",
+      "\t qiskit\t\t CNOTs = 24\t depth = 47\n",
+      "\n",
+      "n = 4\n",
+      "\t knill\t\t CNOTs = 192\t depth = 289\n",
+      "\t ccd\t\t CNOTs = 57\t depth = 111\n",
+      "\t qiskit\t\t CNOTs = 57\t depth = 111\n",
+      "\n",
+      "n = 5\n",
+      "\t knill\t\t CNOTs = 685\t depth = 1204\n",
+      "\t ccd\t\t CNOTs = 122\t depth = 239\n",
+      "\t qiskit\t\t CNOTs = 122\t depth = 239\n",
+      "\n",
+      "n = 6\n",
+      "\t knill\t\t CNOTs = 1170\t depth = 1567\n",
+      "\t ccd\t\t CNOTs = 251\t depth = 493\n",
+      "\t qiskit\t\t CNOTs = 251\t depth = 493\n",
+      "\n",
+      "n = 7\n",
+      "\t knill\t\t CNOTs = 3380\t depth = 5551\n",
+      "\t ccd\t\t CNOTs = 508\t depth = 1001\n",
+      "\t qiskit\t\t CNOTs = 508\t depth = 1001\n",
+      "\n",
+      "n = 8\n",
+      "\t knill\t\t CNOTs = 5276\t depth = 6713\n",
+      "\t ccd\t\t CNOTs = 1021\t depth = 2019\n",
+      "\t qiskit\t\t CNOTs = 1021\t depth = 2019\n"
+     ]
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "import sys\n",
+    "from scipy.stats import unitary_group\n",
+    "sys.path.append('../../')\n",
+    "from qclib.isometry import decompose\n",
+    "\n",
+    "m = 2\n",
+    "print('\\nm =', m)\n",
+    "\n",
+    "result = {}\n",
+    "for n in range(3, 9):\n",
+    "    print('\\nn =', n)\n",
+    "        \n",
+    "    V = unitary_group.rvs(2**n)[:,:2**m]\n",
+    "\n",
+    "    result[n] = []\n",
+    "    \n",
+    "    circuit1 = decompose(V, scheme='knill')\n",
+    "    circuit2 = decompose(V, scheme='ccd')\n",
+    "    \n",
+    "    A = QuantumRegister(n)\n",
+    "    circuit3 = QuantumCircuit(A)\n",
+    "    circuit3.iso(V, q_input=A[0:m], q_ancillas_for_output=A[m:])\n",
+    "        \n",
+    "    transpiled_circuit1 = transpile(circuit1, basis_gates=['u1','u2','u3', 'cx'], optimization_level=3)\n",
+    "    transpiled_circuit2 = transpile(circuit2, basis_gates=['u1','u2','u3', 'cx'], optimization_level=3)\n",
+    "    transpiled_circuit3 = transpile(circuit3, basis_gates=['u1','u2','u3', 'cx'], optimization_level=3)\n",
+    "    \n",
+    "    n_cx1 = transpiled_circuit1.count_ops()['cx']\n",
+    "    n_dp1 = transpiled_circuit1.depth()\n",
+    "    n_cx2 = transpiled_circuit2.count_ops()['cx']\n",
+    "    n_dp2 = transpiled_circuit2.depth()\n",
+    "    n_cx3 = transpiled_circuit3.count_ops()['cx']\n",
+    "    n_dp3 = transpiled_circuit3.depth()\n",
+    "    \n",
+    "    result[n].append(['knill', n_cx1, n_dp1])\n",
+    "    result[n].append(['ccd', n_cx2, n_dp2])\n",
+    "    result[n].append(['qiskit', n_cx3, n_dp3])\n",
+    "    \n",
+    "    print('\\t {0}\\t\\t CNOTs = {1}\\t depth = {2}'.format('knill', n_cx1, n_dp1))\n",
+    "    print('\\t {0}\\t\\t CNOTs = {1}\\t depth = {2}'.format('ccd', n_cx2, n_dp2))\n",
+    "    print('\\t {0}\\t\\t CNOTs = {1}\\t depth = {2}'.format('qiskit', n_cx3, n_dp3))\n",
+    "        "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "separated-pound",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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.7.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/example/isometry/experiments.ipynb b/example/isometry/experiments.ipynb
new file mode 100644
index 00000000..f03e7092
--- /dev/null
+++ b/example/isometry/experiments.ipynb
@@ -0,0 +1,507 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "canadian-reporter",
+   "metadata": {},
+   "source": [
+    "# Qiskit config"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "stylish-integer",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit import IBMQ, QuantumCircuit, ClassicalRegister, QuantumRegister, execute, Aer, transpile\n",
+    "from qiskit.test.mock import FakeMontreal, FakeMumbai\n",
+    "\n",
+    "IBMQ.load_account()\n",
+    "provider = IBMQ.get_provider(hub='ibm-q-research') # 'ibm-q'\n",
+    "\n",
+    "backend_0  = Aer.get_backend('qasm_simulator') \n",
+    "backend_1  = provider.get_backend('ibmq_qasm_simulator')\n",
+    "backend_2  = provider.get_backend('ibmq_santiago')\n",
+    "backend_3  = provider.get_backend('ibmq_manila')\n",
+    "backend_4  = provider.get_backend('ibmq_casablanca')\n",
+    "backend_5  = provider.get_backend('ibmq_jakarta')\n",
+    "backend_6  = FakeMontreal()\n",
+    "backend_7  = FakeMumbai()\n",
+    "\n",
+    "shots      = 8192"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "id": "moved-heritage",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"\n",
+    "    Select the backends that will be compared.\n",
+    "\"\"\"\n",
+    "#backends = [ backend_0 ]\n",
+    "#backends = [ backend_0, backend_1 ]\n",
+    "#backends = [ backend_2, backend_3 ]\n",
+    "#backends = [ backend_4, backend_5 ]\n",
+    "backends = [ backend_6, backend_7 ]\n",
+    "#backends = [ backend_2, backend_3, backend_4, backend_5 ]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "downtown-analysis",
+   "metadata": {},
+   "source": [
+    "# Experiment procedures"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "fleet-albania",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import sys\n",
+    "sys.path.append('../../')\n",
+    "from qclib.state_preparation.schmidt import initialize\n",
+    "\n",
+    "def measurement(circuit, c):\n",
+    "    n = len(c)\n",
+    "    circuit.measure(list(range(n)), c)\n",
+    "\n",
+    "    job = execute(circuit, backend, shots=shots, optimization_level=3)\n",
+    "    \n",
+    "    counts = job.result().get_counts(circuit)\n",
+    "    v = sum(counts.values())\n",
+    "    \n",
+    "    counts2 = {}\n",
+    "    for m in range(2**n):\n",
+    "        pattern = '{:0{}b}'.format(m, n)\n",
+    "        if pattern in counts:\n",
+    "            counts2[pattern] = counts[pattern]\n",
+    "        else:\n",
+    "            counts2[pattern] = 0.0\n",
+    "            \n",
+    "    return { key : value/v for (key, value) in counts2.items() }\n",
+    "    \n",
+    "def experiment(state, low_rank=0, isometry_scheme='knill'):\n",
+    "    circuit = initialize(state, low_rank=low_rank, isometry_scheme=isometry_scheme)\n",
+    "\n",
+    "    n = int(np.log2(len(state)))\n",
+    "    c = ClassicalRegister(n)\n",
+    "    circuit.barrier()\n",
+    "    circuit.add_register(c)\n",
+    "\n",
+    "    return measurement(circuit, c)\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "another-decision",
+   "metadata": {},
+   "source": [
+    "# Experiment"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "id": "painted-clone",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "n = 3\n",
+      "\n",
+      "Input vector probabilities (|x|^2):\n",
+      "{'000': 0.03, '001': 0.06, '010': 0.15000000000000002, '011': 0.049999999999999996, '100': 0.1, '101': 0.29999999999999993, '110': 0.19999999999999998, '111': 0.10999999999999999}\n",
+      "\n",
+      "m = 0\n",
+      "\n",
+      "\tExperiments using fake_montreal backend, with 27 qubits available.\n",
+      "\tIter: 0\n",
+      "\t {'000': 0.0474853515625, '001': 0.1016845703125, '010': 0.0975341796875, '011': 0.048095703125, '100': 0.1083984375, '101': 0.2457275390625, '110': 0.2293701171875, '111': 0.1217041015625}\n",
+      "\tIter: 1\n",
+      "\t {'000': 0.04541015625, '001': 0.104248046875, '010': 0.0914306640625, '011': 0.0518798828125, '100': 0.1190185546875, '101': 0.2445068359375, '110': 0.226806640625, '111': 0.11669921875}\n",
+      "\tIter: 2\n",
+      "\t {'000': 0.044921875, '001': 0.09716796875, '010': 0.0958251953125, '011': 0.053955078125, '100': 0.1068115234375, '101': 0.2503662109375, '110': 0.230224609375, '111': 0.1207275390625}\n",
+      "\tIter: 3\n",
+      "\t {'000': 0.0465087890625, '001': 0.10400390625, '010': 0.09765625, '011': 0.0467529296875, '100': 0.106689453125, '101': 0.2412109375, '110': 0.235107421875, '111': 0.1220703125}\n",
+      "\tIter: 4\n",
+      "\t {'000': 0.04443359375, '001': 0.1048583984375, '010': 0.093994140625, '011': 0.048095703125, '100': 0.111328125, '101': 0.2403564453125, '110': 0.2296142578125, '111': 0.1273193359375}\n",
+      "\tIter: 5\n",
+      "\t {'000': 0.0465087890625, '001': 0.1044921875, '010': 0.09228515625, '011': 0.0499267578125, '100': 0.1136474609375, '101': 0.2376708984375, '110': 0.2288818359375, '111': 0.1265869140625}\n",
+      "\tIter: 6\n",
+      "\t {'000': 0.045166015625, '001': 0.1033935546875, '010': 0.097900390625, '011': 0.0501708984375, '100': 0.112060546875, '101': 0.2374267578125, '110': 0.2333984375, '111': 0.1204833984375}\n",
+      "\tIter: 7\n",
+      "\t {'000': 0.0462646484375, '001': 0.0938720703125, '010': 0.0975341796875, '011': 0.04833984375, '100': 0.1109619140625, '101': 0.2447509765625, '110': 0.2337646484375, '111': 0.12451171875}\n",
+      "\tIter: 8\n",
+      "\t {'000': 0.04296875, '001': 0.1033935546875, '010': 0.103759765625, '011': 0.0499267578125, '100': 0.112548828125, '101': 0.240966796875, '110': 0.23291015625, '111': 0.113525390625}\n",
+      "\tIter: 9\n",
+      "\t {'000': 0.0479736328125, '001': 0.097900390625, '010': 0.0931396484375, '011': 0.052978515625, '100': 0.11376953125, '101': 0.245361328125, '110': 0.22802734375, '111': 0.120849609375}\n",
+      "\n",
+      "\tExperiments using fake_mumbai backend, with 27 qubits available.\n",
+      "\tIter: 0\n",
+      "\t {'000': 0.0418701171875, '001': 0.101806640625, '010': 0.102294921875, '011': 0.0543212890625, '100': 0.1065673828125, '101': 0.2423095703125, '110': 0.2342529296875, '111': 0.1165771484375}\n",
+      "\tIter: 1\n",
+      "\t {'000': 0.0428466796875, '001': 0.099365234375, '010': 0.09765625, '011': 0.0499267578125, '100': 0.11181640625, '101': 0.2359619140625, '110': 0.2451171875, '111': 0.1173095703125}\n",
+      "\tIter: 2\n",
+      "\t {'000': 0.044189453125, '001': 0.1058349609375, '010': 0.0970458984375, '011': 0.05126953125, '100': 0.1046142578125, '101': 0.249755859375, '110': 0.2337646484375, '111': 0.113525390625}\n",
+      "\tIter: 3\n",
+      "\t {'000': 0.046630859375, '001': 0.104736328125, '010': 0.099609375, '011': 0.0484619140625, '100': 0.10400390625, '101': 0.2352294921875, '110': 0.2392578125, '111': 0.1220703125}\n",
+      "\tIter: 4\n",
+      "\t {'000': 0.043701171875, '001': 0.1024169921875, '010': 0.093994140625, '011': 0.0462646484375, '100': 0.107421875, '101': 0.2518310546875, '110': 0.2322998046875, '111': 0.1220703125}\n",
+      "\tIter: 5\n",
+      "\t {'000': 0.04150390625, '001': 0.10009765625, '010': 0.1016845703125, '011': 0.0489501953125, '100': 0.11376953125, '101': 0.244873046875, '110': 0.227294921875, '111': 0.121826171875}\n",
+      "\tIter: 6\n",
+      "\t {'000': 0.0430908203125, '001': 0.103515625, '010': 0.0997314453125, '011': 0.04833984375, '100': 0.10693359375, '101': 0.2447509765625, '110': 0.2340087890625, '111': 0.11962890625}\n",
+      "\tIter: 7\n",
+      "\t {'000': 0.046142578125, '001': 0.0985107421875, '010': 0.099609375, '011': 0.050048828125, '100': 0.1065673828125, '101': 0.240966796875, '110': 0.237548828125, '111': 0.12060546875}\n",
+      "\tIter: 8\n",
+      "\t {'000': 0.04248046875, '001': 0.1024169921875, '010': 0.1002197265625, '011': 0.0487060546875, '100': 0.1068115234375, '101': 0.2415771484375, '110': 0.232666015625, '111': 0.1251220703125}\n",
+      "\tIter: 9\n",
+      "\t {'000': 0.0484619140625, '001': 0.100830078125, '010': 0.0975341796875, '011': 0.0465087890625, '100': 0.103515625, '101': 0.2471923828125, '110': 0.2314453125, '111': 0.12451171875}\n",
+      "\n",
+      "m = 1\n",
+      "\n",
+      "\tExperiments using fake_montreal backend, with 27 qubits available.\n",
+      "\tIter: 0\n",
+      "\t {'000': 0.0401611328125, '001': 0.0675048828125, '010': 0.1412353515625, '011': 0.0623779296875, '100': 0.116455078125, '101': 0.27099609375, '110': 0.1781005859375, '111': 0.1231689453125}\n",
+      "\tIter: 1\n",
+      "\t {'000': 0.0438232421875, '001': 0.0782470703125, '010': 0.1505126953125, '011': 0.0517578125, '100': 0.1033935546875, '101': 0.2640380859375, '110': 0.1923828125, '111': 0.1158447265625}\n",
+      "\tIter: 2\n",
+      "\t {'000': 0.0396728515625, '001': 0.0687255859375, '010': 0.14599609375, '011': 0.0567626953125, '100': 0.107421875, '101': 0.2784423828125, '110': 0.1927490234375, '111': 0.1102294921875}\n",
+      "\tIter: 3\n",
+      "\t {'000': 0.0357666015625, '001': 0.070068359375, '010': 0.1544189453125, '011': 0.0478515625, '100': 0.1072998046875, '101': 0.290771484375, '110': 0.1905517578125, '111': 0.103271484375}\n",
+      "\tIter: 4\n",
+      "\t {'000': 0.0360107421875, '001': 0.0750732421875, '010': 0.1448974609375, '011': 0.0560302734375, '100': 0.1090087890625, '101': 0.28173828125, '110': 0.1900634765625, '111': 0.107177734375}\n",
+      "\tIter: 5\n",
+      "\t {'000': 0.0440673828125, '001': 0.0765380859375, '010': 0.1566162109375, '011': 0.055908203125, '100': 0.1021728515625, '101': 0.271240234375, '110': 0.18505859375, '111': 0.1083984375}\n",
+      "\tIter: 6\n",
+      "\t {'000': 0.03271484375, '001': 0.0625, '010': 0.145751953125, '011': 0.05029296875, '100': 0.116455078125, '101': 0.279296875, '110': 0.1990966796875, '111': 0.1138916015625}\n",
+      "\tIter: 7\n",
+      "\t {'000': 0.040771484375, '001': 0.0654296875, '010': 0.1376953125, '011': 0.05029296875, '100': 0.1085205078125, '101': 0.2808837890625, '110': 0.1923828125, '111': 0.1240234375}\n",
+      "\tIter: 8\n",
+      "\t {'000': 0.035400390625, '001': 0.0609130859375, '010': 0.1453857421875, '011': 0.052001953125, '100': 0.107177734375, '101': 0.2916259765625, '110': 0.201416015625, '111': 0.1060791015625}\n",
+      "\tIter: 9\n",
+      "\t {'000': 0.034912109375, '001': 0.064697265625, '010': 0.1466064453125, '011': 0.0556640625, '100': 0.1162109375, '101': 0.272216796875, '110': 0.2015380859375, '111': 0.108154296875}\n",
+      "\n",
+      "\tExperiments using fake_mumbai backend, with 27 qubits available.\n",
+      "\tIter: 0\n",
+      "\t {'000': 0.0369873046875, '001': 0.0672607421875, '010': 0.1500244140625, '011': 0.059814453125, '100': 0.112548828125, '101': 0.264892578125, '110': 0.1954345703125, '111': 0.113037109375}\n",
+      "\tIter: 1\n",
+      "\t {'000': 0.0447998046875, '001': 0.0726318359375, '010': 0.147216796875, '011': 0.055908203125, '100': 0.1099853515625, '101': 0.2716064453125, '110': 0.1895751953125, '111': 0.1082763671875}\n",
+      "\tIter: 2\n",
+      "\t {'000': 0.03955078125, '001': 0.0645751953125, '010': 0.1536865234375, '011': 0.0484619140625, '100': 0.1058349609375, '101': 0.2861328125, '110': 0.19091796875, '111': 0.11083984375}\n",
+      "\tIter: 3\n",
+      "\t {'000': 0.04052734375, '001': 0.0712890625, '010': 0.1422119140625, '011': 0.0555419921875, '100': 0.1136474609375, '101': 0.27197265625, '110': 0.18115234375, '111': 0.1236572265625}\n",
+      "\tIter: 4\n",
+      "\t {'000': 0.0361328125, '001': 0.06884765625, '010': 0.141845703125, '011': 0.0574951171875, '100': 0.1092529296875, '101': 0.287109375, '110': 0.1923828125, '111': 0.10693359375}\n",
+      "\tIter: 5\n",
+      "\t {'000': 0.0372314453125, '001': 0.0667724609375, '010': 0.1331787109375, '011': 0.0584716796875, '100': 0.1124267578125, '101': 0.2801513671875, '110': 0.1922607421875, '111': 0.1195068359375}\n",
+      "\tIter: 6\n",
+      "\t {'000': 0.0399169921875, '001': 0.06494140625, '010': 0.1492919921875, '011': 0.0562744140625, '100': 0.1021728515625, '101': 0.2763671875, '110': 0.199951171875, '111': 0.111083984375}\n",
+      "\tIter: 7\n",
+      "\t {'000': 0.0418701171875, '001': 0.0657958984375, '010': 0.151611328125, '011': 0.054443359375, '100': 0.1046142578125, '101': 0.27685546875, '110': 0.190673828125, '111': 0.1141357421875}\n",
+      "\tIter: 8\n",
+      "\t {'000': 0.0401611328125, '001': 0.0653076171875, '010': 0.13525390625, '011': 0.0582275390625, '100': 0.11962890625, '101': 0.275390625, '110': 0.1917724609375, '111': 0.1142578125}\n",
+      "\tIter: 9\n",
+      "\t {'000': 0.034912109375, '001': 0.0654296875, '010': 0.148193359375, '011': 0.053466796875, '100': 0.1077880859375, '101': 0.2816162109375, '110': 0.201904296875, '111': 0.106689453125}\n"
+     ]
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "\"\"\"\n",
+    "    Random complex input vector.\n",
+    "\"\"\"\n",
+    "#n = 3\n",
+    "#N = 2**n\n",
+    "#rng = np.random.default_rng(42)\n",
+    "#input_state = ((1.0 - 0.01) * rng.random(N) + 0.01) + ((1.0 - 0.01) * rng.random(N) + 0.01) * 1j\n",
+    "#input_state = input_state/np.linalg.norm(input_state)\n",
+    "\"\"\"\"\"\"\n",
+    "\"\"\"\n",
+    "    Fixed real 8- and 16-dim input vector.\n",
+    "\"\"\"\n",
+    "input_state = [np.sqrt(0.03), np.sqrt(0.06), np.sqrt(0.15), np.sqrt(0.05), np.sqrt(0.1), np.sqrt(0.3), np.sqrt(0.2), np.sqrt(0.11) ]\n",
+    "#input_state = [np.sqrt(0.01), np.sqrt(0.02), np.sqrt(0.04), np.sqrt(0.02), np.sqrt(0.07), np.sqrt(0.08), np.sqrt(0.04), np.sqrt(0.01), np.sqrt(0.08), np.sqrt(0.02), np.sqrt(0.21), np.sqrt(0.09), np.sqrt(0.12), np.sqrt(0.08), np.sqrt(0.05), np.sqrt(0.06) ]\n",
+    "N = len(input_state)\n",
+    "n = int(np.log2(N))\n",
+    "\"\"\"\"\"\"\n",
+    "\n",
+    "print('n =', n)\n",
+    "\n",
+    "ideal = np.power(np.abs(input_state),2)\n",
+    "print('\\nInput vector probabilities (|x|^2):')\n",
+    "print({'{:0{}b}'.format(p, n) : a for p, a in enumerate(ideal)})\n",
+    "\n",
+    "maes  = {}\n",
+    "probs = {}\n",
+    "for m in range(0, 2):\n",
+    "    print('\\nm =', m)\n",
+    "    backend_probs = {}\n",
+    "    backend_maes = {}\n",
+    "    for j, backend in enumerate(backends):\n",
+    "        backend_name = backend.name()\n",
+    "        backend_config = backend.configuration()\n",
+    "        backend_qubits = backend_config.n_qubits\n",
+    "\n",
+    "        backend_probs[backend_name] = []\n",
+    "\n",
+    "        print('\\n\\tExperiments using {0} backend, with {1} qubits available.'.format(backend_name, backend_qubits))\n",
+    "\n",
+    "        for k in range(10):\n",
+    "            print('\\tIter:', k)\n",
+    "            prob = experiment(input_state, low_rank=2**m)\n",
+    "            backend_probs[backend_name].append([val for key, val in prob.items()])\n",
+    "            print('\\t', dict(sorted(prob.items())) )\n",
+    "\n",
+    "        backend_maes[backend_name] = np.sum(np.abs(np.mean(backend_probs[backend_name],axis=0)-ideal))/len(ideal)\n",
+    "        \n",
+    "    probs[m] = backend_probs\n",
+    "    maes[m]  = backend_maes\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "damaged-statistics",
+   "metadata": {},
+   "source": [
+    "# Mean Absolute Error"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "comic-circuit",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "s = 1\n",
+      "\tfake_montreal : MAE = 0.04395446777343749\n",
+      "\tfake_mumbai : MAE = 0.0486846923828125\n",
+      "s = 2\n",
+      "\tfake_montreal : MAE = 0.026739501953124998\n",
+      "\tfake_mumbai : MAE = 0.039910888671875\n",
+      "s = 3\n",
+      "\tfake_montreal : MAE = 0.011114501953124996\n",
+      "\tfake_mumbai : MAE = 0.024185180664062497\n"
+     ]
+    }
+   ],
+   "source": [
+    "for s, s_mae in maes.items():\n",
+    "    print('s = {0}'.format(s))\n",
+    "    for backend, mae in dict(sorted(s_mae.items(), key=lambda item: item[1])).items():\n",
+    "        print('\\t{0} : MAE = {1}'.format(backend, mae))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "western-forest",
+   "metadata": {},
+   "source": [
+    "# Chart Data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 12,
+   "id": "victorian-sleeve",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "s = 1\n",
+      "\tfake_montreal\n",
+      "\t[0.08218994\t0.08010254\t0.11359863\t0.07648926\t0.15465088\t0.20549316\n",
+      "\t0.15509033\t0.13238525]\n",
+      "\tfake_mumbai\n",
+      "\t[0.07937012\t0.07935791\t0.11169434\t0.08822021\t0.1640625\t\t0.18762207\n",
+      "\t0.15594482\t0.13372803]\n",
+      "s = 2\n",
+      "\tfake_montreal\n",
+      "\t[0.06756592\t0.08161621\t0.14246826\t0.06879883\t0.12402344\t0.23200684\n",
+      "\t0.16856689\t0.11495361]\n",
+      "\tfake_mumbai\n",
+      "\t[0.08748779\t0.08448486\t0.13027344\t0.0769165\t\t0.13448486\t0.20407715\n",
+      "\t0.15600586\t0.12626953]\n",
+      "s = 3\n",
+      "\tfake_montreal\n",
+      "\t[0.04702148\t0.07434082\t0.14597168\t0.05401611\t0.10860596\t0.27253418\n",
+      "\t0.18703613\t0.11047363]\n",
+      "\tfake_mumbai\n",
+      "\t[0.06831055\t0.08980713\t0.14331055\t0.06899414\t0.10761719\t0.23793945\n",
+      "\t0.17200928\t0.11201172]\n"
+     ]
+    }
+   ],
+   "source": [
+    "for s, s_probs in probs.items():\n",
+    "    print('s = {0}'.format(s))\n",
+    "    for backend, backend_probs in s_probs.items():\n",
+    "        print('\\t{0}'.format(backend))\n",
+    "        print('\\t{0}'.format(str(np.mean(backend_probs, axis=0)).replace(' ', '\\t')))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "nuclear-bobby",
+   "metadata": {},
+   "source": [
+    "# Plot Charts"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "auburn-token",
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "s = 1\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAI4CAYAAAB3HEhGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA+5UlEQVR4nO3deZhcVZ3/8fcXgiRkAUIgSEhIs4pBIhJRx0GaKFFHWRwkLBkUVBgQBhlGBxSBAokguKECiqOIiCCL4sIixhBmlPkFok4cIrskRDIsAYGELUC+vz+qEpumk3R1TlWlOu/X89RD1a1Tp791vGnvp8+950ZmIkmSJElafeu0ugBJkiRJ6i8MWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSpkQKsLaLYRI0bk2LFjW12GJEmSpDb2u9/9bmFmbtp9+1oXsMaOHcusWbNaXYYkSZKkNhYR83ra7imCkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKmStW+RCkiRJ6q2nn36aRx99lBdffLHVpajJ1ltvPTbbbDOGDRtW1+cMWJIkSVIPnn76aR555BFGjRrFoEGDiIhWl6QmyUyee+45HnroIYC6QpanCEqSJEk9ePTRRxk1ahQbbLCB4WotExFssMEGjBo1ikcffbSuzxqwJEmSpB68+OKLDBo0qNVlqIUGDRpU9+mhBixJkiRpBZy5Wrv15X9/A5YkSZIkFWLAkiRJktrQuHHjmDFjxqu2z5gxgy233LLIz4gI7rvvviJ9rS1cRVCSJEnqpd2mTmto/7ed/K5et50zZ04DK1FfOYMlSZIkSYUYsCRJkqQ2NHbsWKZNm8Zzzz3HYYcdxsYbb8zrX/96br/99le0W7BgAfvvvz+bbropHR0dfO1rX1v+3m233cbb3vY2NtpoI1772tdy7LHHsmTJkmZ/lX7FgCVJkiS1sdNPP53777+f+++/n1/+8pdccskly99bunQpe++9N+PHj+ehhx7i17/+NV/96lf55S9/CcC6667LV77yFRYuXMh///d/8+tf/5oLLrigVV+lXzBgSZIkSW3syiuv5OSTT2b48OGMHj2a4447bvl7t99+O4899hinnnoqr3nNa9h666054ogjuOKKKwDYddddeetb38qAAQMYO3Ys//zP/8wtt9zSqq/SL7jIhSRJktTGFixYwOjRo5e/3mqrrZY/nzdvHgsWLGCjjTZavu3ll19m9913B+Cee+7hhBNOYNasWTz77LO89NJL7Lrrrk2rvT9yBkuSJElqY6997WuZP3/+8tcPPvjg8uejR4+mo6ODJ598cvlj0aJFXH/99QAcffTRvO51r+Pee+/l6aef5vOf/zyZ2fTv0J8YsCRJkqQ2NnnyZM466yz++te/8pe//IWvf/3ry9/bbbfdGDZsGF/4whd47rnnePnll7njjjuWL4SxaNEihg0bxpAhQ7jrrru48MILW/U1+g1PEZQkSZJ6qZ77VDXLaaedxlFHHUVHRwdbbLEFhx9+OOeddx5QXcTi5z//Of/2b/9GR0cHL7zwAjvssANnnnkmAF/84hc58sgjOeecc9hll1048MADmT59eiu/TtuLZk8BRsRw4DvAJGAh8OnM/GEP7Q4CTgc2B14AbgD+JTOfrqef7iZMmJCzZs0q9G0kSZLUX915553suOOOrS5DLbai/SAifpeZE7pvb8UpgucDS4CRwBTgwogY10O73wJvz8wNga2pzrad2Yd+JEmSJKkpmhqwImIwsD9wSmYuzszfAD8DDu3eNjPnZ+bCLpteBrattx9JkiRJapZmX4O1PfByZt7TZdtsYI+eGkfE3wPXAcOAZ4EP9LGfI4EjAUaOHMmMGTNW4ytIkqRWuuvhRQB8+YxPc++dc4r1u92O47jogm8U60/tb8MNN2TRokWtLkMt9vzzz9eVH5odsIYAT3Xb9hQwtKfGtZmpDSNiFHAEMLeP/VwEXATVa7A6Ozv7ULokSVoT/PvUaQBsNOWrvLlw3x4jqKs777yToUN7PLzUWmTgwIHssssuvW7f7GuwFlOdjepqGLDSPw1k5kPAjcAVq9OPJEmSJDVSswPWPcCAiNiuy7bxQG/m9wcA2xToR5IkSZIaoqkBKzOfAX4MnBERgyPi7cC+wKXd20bElIgYE1VbAVOBX9fbjyRJkiQ1SyuWaf84MAh4FLgcODoz59TC1OKIGFNr93rgVqqnA/4WuJvqdVgr7adJ30GSJEmSXqXZi1yQmU8A+/Ww/UGqi1cse30ycHK9/UiSJElSqzQ9YEmSJElt64s7NLb/T97d66Z33303Bx10EPfddx9Tp07luOOO67Hd3Llz6ejo4MUXX2TAAA//u2rE2DjCkiRJUhs655xz6Ozs5A9/+EOrS2mYSqXCfffdxw9+8INWl9JrrbgGS5IkSdJqmjdvHuPGjWt1GS2VmSxdurTVZbyCAUuSJElqMxMnTuTmm2/m2GOPZciQIZx33nnssssuDBs2jNGjR1OpVFb42WuuuYaxY8dyxx13sHTpUs4++2y22WYbNtlkEyZPnswTTzyx0p89d+5cIoKLL76Y0aNHs/HGG/PNb36T22+/nZ133pmNNtqIY489dnn7pUuXcuaZZ7LVVlux2Wab8aEPfYinnnrqFX1dcskljBkzhhEjRjB16lQAbrzxRj7/+c/zox/9iCFDhjB+/HigekPwk08+mbe//e1ssMEG/PnPf+auu+5ir732Yvjw4eywww5ceeWVy3/+dddd1+uxKcGAJUmSJLWZ6dOns/vuu/ONb3yDxYsXM378eL7//e/z5JNPct1113HhhRdy7bXXvupzF198MSeeeCLTpk1jp5124mtf+xrXXnstt9xyCwsWLGDjjTfmmGOO6VUNM2fO5N577+VHP/oRxx9/PFOnTmXatGnMmTOHK6+8kltuuQWA733ve3zve9/j5ptv5s9//jOLFy9+RQAD+M1vfsPdd9/Nr3/9a8444wzuvPNO3vOe9/CZz3yGAw88kMWLFzN79uzl7S+99FIuuugiFi1axKabbspee+3FIYccwqOPPsrll1/Oxz/+cebMqS4wPnjw4F6NTSkGLEmSJKnNdXZ28oY3vIF11lmHnXfemYMPPnh5wFnmq1/9Kueeey4zZsxg2223BeBb3/oWU6dOZcstt2T99denUqlw9dVX89JLL63yZ55yyikMHDiQSZMmMXjwYA4++GA222wzRo0axe6777782rDLLruME044ga233pohQ4Zw1llnccUVV7ziZ5x22mkMGjSI8ePHM378+FeEqZ4cdthhjBs3jgEDBnDjjTcyduxYDj/8cAYMGMCb3vQm9t9/f66++upej01JLnIhSZIktbmZM2dy0kkncccdd7BkyRJeeOEFDjjggFe0Offcczn11FPZcsstl2+bN28eH/jAB1hnnb/Nu6y77ro88sgjjBo1aqU/c+TIkcufDxo06FWvFy9eDMCCBQvYaqutlr+31VZb8dJLL/HII48s37b55psvf77BBhss/+yKjB49+hXfYebMmWy00UbLt7300ksceuihQO/GpiRnsCRJkqQ2d8ghh7DPPvswf/58nnrqKY466igy8xVtbrrpJs4880yuueaa5dtGjx7NDTfcwJNPPrn88fzzz68yXNVjiy22YN68ectfP/jggwwYMOAVgWxFImKV20ePHs0ee+zxiu+wePFiLrzwQqB3Y1OSM1iSJElSb9Vxn6pmWrRoEcOHD2fgwIHcdttt/PCHP2TSpEmvaDNu3DhuvPFG3v3ud7Peeuuxzz77cNRRR3HyySdzySWXsNVWW/HYY49x6623su+++xar7eCDD+YLX/gC733ve9l0002XX1fVm/tOjRw5kl/96lcsXbr0FbNsXb3//e/npJNO4tJLL+Wggw4C4H/+538YMmQIO+64Y6/GpiRnsCRJkqQ2d8EFF3DqqacydOhQzjjjDCZPntxju/Hjx/OLX/yCI444ghtuuIFPfOIT7LPPPkyaNImhQ4fy1re+lZkzZxat7SMf+QiHHnoo73jHO+jo6GDgwIF8/etf79Vnl53Kt8kmm/CmN72pxzZDhw7lpptu4oorrmCLLbZg880358QTT+SFF14Aej82pUQjp8fWRBMmTMhZs2a1ugxJktRHu02d1rC+bzv5XQ3rW+3nzjvvZMcdd2x1GWqxFe0HEfG7zJzQfbszWJIkSZJUiAFLkiRJ0itcdtllDBky5FWPcePGtbq0NZ6LXEiSJEl6hSlTpjBlypRWl9GWnMGSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEVQQlSZKkXpp45cSG9j998vRet7377rs56KCDuO+++5g6dSrHHXdcj+3mzp1LR0cHL774IgMGtP/hf0Rw7733su2229b92XHjxnH++efT2dlZvrCa9h9hSZIkaS10zjnn0NnZyR/+8IdWl9I25syZ0/Cf4SmCkiRJUhuaN2+eN/5dAxmwJEmSpDYzceJEbr75Zo499liGDBnCeeedxy677MKwYcMYPXo0lUplhZ+95pprGDt2LHfccQdLly7l7LPPZptttmGTTTZh8uTJPPHEEyv92XPnziUiuPjiixk9ejQbb7wx3/zmN7n99tvZeeed2WijjTj22GOXt69UKvzTP/3Tqz7/0ksvAdDZ2clnP/tZ/u7v/o4hQ4aw99578/jjjzNlyhSGDRvGm9/8ZubOnfuKGq6//nq23nprRowYwac+9SmWLl0KwP3338/EiRPZZJNNGDFiBFOmTOHJJ59c/rmxY8cybdq0Xo5y3xiwJEmSpDYzffp0dt99d77xjW+wePFixo8fz/e//32efPJJrrvuOi688EKuvfbaV33u4osv5sQTT2TatGnstNNOfO1rX+Paa6/llltuYcGCBWy88cYcc8wxvaph5syZ3HvvvfzoRz/i+OOPZ+rUqUybNo05c+Zw5ZVXcsstt/T6+1xxxRVceumlPPTQQ9x///287W1v4/DDD+eJJ55gxx135PTTT39F+5/85CfMmjWL3//+9/z0pz/lu9/9LgCZyac//WkWLFjAnXfeyfz581caNhvBgCVJkiS1uc7OTt7whjewzjrrsPPOO3PwwQe/KuB89atf5dxzz2XGjBnLF4j41re+xdSpU9lyyy1Zf/31qVQqXH311ctnl1bmlFNOYeDAgUyaNInBgwdz8MEHs9lmmzFq1Ch23333uq4NO/zww9lmm23YcMMNee9738s222zDu971LgYMGMABBxzwqr5OPPFEhg8fzpgxYzj++OO5/PLLAdh2223Za6+9WH/99dl000054YQT6gp6JbjIhSRJktTmZs6cyUknncQdd9zBkiVLeOGFFzjggANe0ebcc8/l1FNPZcstt1y+bd68eXzgAx9gnXX+Nu+y7rrr8sgjjzBq1KiV/syRI0cufz5o0KBXvV68eHGv66+3r9GjRy9/vtVWW7FgwQIAHn30UY477jj+67/+i0WLFrF06VI23njjXtdRgjNYkiRJUps75JBD2GeffZg/fz5PPfUURx11FJn5ijY33XQTZ555Jtdcc83ybaNHj+aGG27gySefXP54/vnnVxmu6jF48GCeffbZ5a8ffvjh1e5z/vz5y58/+OCDbLHFFgB8+tOfJiL44x//yNNPP80PfvCDV41DozmDJUmSJPVSPfepaqZFixYxfPhwBg4cyG233cYPf/hDJk2a9Io248aN48Ybb+Td73436623Hvvssw9HHXUUJ598MpdccglbbbUVjz32GLfeeiv77rtvsdre+MY38oUvfIEHH3yQDTfckLPOOmu1+zz33HN5y1vewuLFiznvvPM44YQTgOo4bLjhhmy00UY89NBDnHvuuav9s+rlDJYkSZLU5i644AJOPfVUhg4dyhlnnMHkyZN7bDd+/Hh+8YtfcMQRR3DDDTfwiU98gn322YdJkyYxdOhQ3vrWtzJz5syite21114ceOCB7Lzzzuy66668//3vX+0+9913X3bddVfe+MY38r73vY+PfvSjAJx22mn8/ve/Z8MNN+R973sf//iP/7jaP6te0ewps1abMGFCzpo1q9VlSJKkPtptauOWWL7t5Hc1rG+1nzvvvJMdd9yx1WWoxVa0H0TE7zJzQvftzmBJkiRJUiEGLEmSJEmvcNlllzFkyJBXPcaNG9fq0tZ4LnIhSZIk6RWmTJnClClTWl1GW3IGS5IkSVqBpUuXtroEtVBf/vc3YEmSJEk9GDx4MA899BBLlixp+r2U1FqZyZIlS3jooYcYPHhwXZ/1FEFJkiSpB1tuuSULFy5k3rx5vPTSS60uR002YMAANtxwQ0aMGFHf5xpUjyRJktTW1llnHTbbbDM222yzVpeiNuIpgpIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqpOkBKyKGR8RPIuKZiJgXEYesoN2HI+J3EfF0RPwlIs6JiAFd3p8REc9HxOLa4+7mfQtJkiRJerVWzGCdDywBRgJTgAsjYlwP7TYAjgdGAG8B3gl8slubYzNzSO2xQ+NKliRJkqRVG7DqJuVExGBgf2CnzFwM/CYifgYcCpzUtW1mXtjl5UMRcRmwZ9OKlSRJkqQ6NTVgAdsDL2fmPV22zQb26MVn3wHM6bbtrIg4G7gbODkzZ/T0wYg4EjgSYOTIkcyY0WMzSZLUBj6yzXMN69tjBEmrq9kBawjwVLdtTwFDV/ahiDgcmAB8rMvmE4E/UT3d8CDg5xHxxsy8v/vnM/Mi4CKACRMmZGdnZ1/rlyRJLfbvU6c1rO/bDupsWN+S1g7NvgZrMTCs27ZhwKIVfSAi9gPOBt6bmQuXbc/MmZm5KDNfyMxLgN8C/1C+ZEmSJEnqnWYHrHuAARGxXZdt43n1qX8ARMR7gG8De2fm/66i7wSiSJWSJEmS1AdNDViZ+QzwY+CMiBgcEW8H9gUu7d42IiYClwH7Z+Zt3d7bKCLeHREDI2JAREyheo3WLxv/LSRJkiSpZ61Ypv3jwCDgUeBy4OjMnBMRY2r3sxpTa3cKsCFwfZd7Xd1Qe2894EzgMWAh8C/AfpnpvbAkSZIktUyzF7kgM58A9uth+4NUF8FY9nqFS7Jn5mPAmxtRnyRJkiT1VStmsCRJkiSpXzJgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJHVTqVSIiGKPSqXS6q8kqUkGtLoASZKkNU2lUlllKOrs7ARgxowZDa9HUvtwBkuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiDcaliRJWuaLO/S+7fz59X3mk3fXX4+ktuMMliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIU0PWBExPCJ+EhHPRMS8iDhkBe0+HBG/i4inI+IvEXFORAyotx9JkiRJapZWzGCdDywBRgJTgAsjYlwP7TYAjgdGAG8B3gl8sg/9SJIkSVJTNDVgRcRgYH/glMxcnJm/AX4GHNq9bWZemJn/lZlLMvMh4DLg7fX2I0mSJEnN0uwbDW8PvJyZ93TZNhvYoxeffQcwpy/9RMSRwJEAI0eOZMaMGXWWLUmS1hQf2ea5hvU9I47rddsn1/9q9TOje/mZNj7+uOvhRXz5jE9z751zVt24l7bbcRwnnHoWr9t8aLE+pTVBZGbzfljE7sBVmbl5l21HAFMys3Mlnzsc+Bzwxsxc2Nd+ACZMmJCzZs1avS8iSZJaZrep0xrW923rHwNA5aaFnP6rJ4r1e9pew6nc9Hix/pqtoWN+8rsa1rfUSBHxu8yc0H17s2ewFgPDum0bBixa0QciYj/gbOBdmbmwr/1IkiT1VmXSCCqTRrS6DEltqNmLXNwDDIiI7bpsG8/fTv17hYh4D/BtYO/M/N++9iNJkiRJzdDUgJWZzwA/Bs6IiMER8XZgX+DS7m0jYiLVhS32z8zb+tqPJEmSJDVLK5Zp/zgwCHgUuBw4OjPnRMSYiFgcEWNq7U4BNgSur21fHBE3rKqf5n0NSZIkSXqlZl+DRWY+AezXw/YHgSFdXu/Zl34kSZIkqVVaMYMlSZIkSf2SAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJElqokqlQkQUe1QqlVZ/JUldNP1Gw5IkSWuzSqWyylDU2dkJwIwZMxpej6SynMGSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIS5yIUmS1AQTr5zY67azH5td12emT57ep5oklecMliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEZdolSZKa6IGrHmDu1XN71fbmA29eZZuxHxwLk1evJknlGLAkSZKaqOOADjoO6Gh1GZIaxFMEJUmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKqStgRcROjSpEkiRJktpdvTNYf4yI2yPi6IjYqBEFSZIkSVK7qjdgvRP4E3AOsCAiLo+IvSIiypcmSZIkSe2lroCVmTdn5oeBzYFjgVHAL4F5EfG5iNimATVKkiRJUlvo0yIXmflMZn43M98BbA/MBT4D3BMRt0TEBwrWKEmSJEltoc+rCEbE2IioADcBbwOuB44EHgF+FBFfKVKhJEmSJLWJelcR3CAiPhQRNwP3AVOAbwNjMnPvzPxOZk4G/hn4aPlyJUmSJGnNNaDO9g8D6wI/Bt6VmTNW0O524PHVqEuSJEmS2k69Aesk4LLMfGpljTLzDqCjz1VJkiRJUhuq9xqsEcDgnt6IiNdGxKmrX5IkSZIktad6A9ZpwJYreG+L2vuSJEmStFaqN2AFkCt4b0vgr6tXjiRJkiS1r1VegxURHwY+XHuZwIUR8XS3ZgOBN1Bdsl2SJEmS1kq9WeTiWf62ImAATwFPdGuzBLgBuKBcaZIkSZLUXlYZsDLzKuAqgIi4GPhcZv650YVJkiRJUrup6xqszDzccCVJUv9RqVSIiGKPSqXS6q8kSS3Vm2uwzgG+lpl/qT1fmczME8uUJkmSGq1SqawyFHV2dgIwY8aMhtcjSe2uN9dgHQBcBvyl9nxlEjBgSZIkSVor9eYarI6enkuSJEmSXqne+2BJkiRJklagN9dg/UM9HWbm9X0vR5IkSZLaV2+uwfoF1WurohdtE1h3tSqSJEmSpDbVm4DldVeSJEmS1Au9WeRiXjMKkSRJkqR215trsDbIzGeXPV9V+2VtJUmSJGlt05tTBBdFxNsy8zZgMdXrrFbGa7AkSZIkrZV6E7A+Atzf5fmqApYkSWoDE6+c2Kt2sx+bXVd7gOmTp/epJklqd725BuuSLs+/19BqJEmSJKmN9WYG61UiYiNgJ+C1wP8Bd2Tmk+XKkiRJkqT2U1fAiogBwFTgGKDrghfPRsQFwMmZ+WLB+iRJkiSpbaxTZ/svA58APg+8HhhR++9ZwHHAl4pWJ0mSJK2GSqVCRBR7VCqVVn8lreHqPUXwUOAzmfnlLtueAKZGxPPAZ6kGLUmSJKnlKpXKKkNRZ2cnADNmzGh4Per/6p3BWgrMWcF7d+AKg5IkSZLWYvUGrEuBj63gvSOAH6xeOZIkSVrbeBqf+pNVBqyI+PiyBzAXeFtEzImIsyLiX2v//RPwFv52vyxJklSIB5/q7yqVCpm50scee+zBHnvsscp2mek+rpbqzTVY3+hh2xbAjj1s/zJw3mpVJEmSXsFrSCSpffTmRsP1nkYoSZLaxANXPcDcq+f2qu3NB968yjZjPziWjgM6VrMqSWpffbrRsCRJ6h86DugwEElSQX0KWBGxJbA9MLD7e5l5/eoWJUmSJEntqK6AFRFDgSuBScs21f7bdXn2dQvUJUmSJK2WiVdO7FW72Y/Nrqs9wPTJ0/tUk/q/emewzgLGALsDvwE+APwV+CdgInBw0eokSZLUv31xh961mz+/vvYAY0bVX4+0muoNWP8AfBaYWXu9IDNvB/4zIr4EfAqYXLA+SZLWXvUcSHrwKUlrhHpXCBwJzM/Ml4FngOFd3ruev506KEmSJElrnXoD1nxgRO35vcD7u7z3FuD5EkVJkiRJUjuq9xTBXwHvAn4CfAW4JCJ2BV4A3gF8qWx5kiRJktQ+6g1YJwIbAGTmpRGxGPggMAg4FvhW2fIkSZIkqX3UdYpgZj6bmQu7vP5JZk7JzH/MzAszc+mq+oiI4RHxk4h4JiLmRcQhK2i3U0T8MiIWRkT28P6MiHg+IhbXHnfX810kSZIkqbS+3mh4B+DNwGuBBcCszOxtwDkfWEJ1wYw3AtdFxOzMnNOt3YtU77l1AXDtCvo6NjP/o77qJUmSJKkx6r3R8DDg28D+VGe/FgNDgKUR8WPgY5n59Eo+P7j22Z0yczHwm4j4GXAocFLXtrXAdndEbFtPjZIkSWovlZsWcvqvnuhV2/jUPatsc9pew6lMqq7L9sBVDzD36rm96vvmA29eZZuxHxxLxwEdvepPa6fIfNXZdytuHPED4H3AMcCPM/P5iBhINTR9A7guM/9pJZ/fBbg1Mwd12fZJYI/M3HsFn9kWuDczo9v2GcA4IIC7gZMzc8YK+jgSOBJg5MiRu15xxRW9+8KSVttdDy8C4MtnfJp77+w+Ud132+04josu+Eax/qQ10iPVfzPfu+o6LrnmhmLdfnj/93LYAe/jntesV6zP7rbfePuG9b3s90ojvC4ebFjf7Tre4Jj3pNFjrjXfnnvu+bvMnNB9e70BaxHwrz2dlhcRRwBfzsyhK/n87sBVmbl5t89NyczOFXxmRQHrLcCfqJ5ueBDVgPfGzLx/Zd9hwoQJOWvWrJU1kVTQblOnNazv205+V8P6ltYI9dw0uA8mNvBGw9MnT29Y3w39vbL+MQ3ru13HGxzznjR6zLXmi4geA1a998FaDPzfCt5bQPXmw6v6/LBu24YBdf9ZJDNnZuaizHwhMy8Bfgv8Q739SJIkSVIp9Qas84FPRsSgrhsjYgPgk1QXpFiZe4ABEbFdl23jgRLnDSXV0wUlSZIkqSVWuchFRJzTbdN2wPyI+BXwKLAZsBfwHLDSc+8y85naYhhnRMTHqK4iuC/wdz383ADWB15Tez2w2kW+EBEbAW8BbgFeAg6keqPj41f1fSRJkiSpUXqziuAB3V6/WHu8tcu2Zaf47Q98ahX9fRz4LtVw9jhwdGbOiYgxVK+pen1mPghsBTzQ5XPPAfOAscB6wJnA64CXgbuA/epYKl6SJEmSiltlwMrMoutQZuYTwH49bH+Q6pLvy17PZQWn/GXmY1TvwyVJkiRJa4x6r8GSJEmSJK1A3QErIraOiAsj4n8j4qHafy+IiK0bUaAkSZIktYveXIO1XETsCtwMPA/8AngEGEn12qspEbFnZv6+eJWSJEmS1AbqCljAF4E/AO/NzGeXbawt03597f2J5cqTJEmSpPZR7ymCuwHndA1XALXXX6S6dLokSZIkrZXqDVjPAZus4L3hVE8dlCRJkqS1Ur0B6zrg7Ij4+64ba6/PAn5eqjBJkiRJajf1XoN1AvBT4JaIeIzqIheb1R63Av9WtjxJkiRJah91BazMfBz4+4h4D9Ub/b4W+D9gZmbe1ID6JEmSJKlt9DpgRcRA4I/AcZl5I3Bjw6qSJEmSpDbU62uwMvN5YCNgacOqkSRJktTWKpUKEVHsUalUWv2V6lLvIheXAYc3ohBJkiRJ7a9SqZCZK33sscce7LHHHqtsl5n9PmA9CLwjImZFxBkRcUxEfLzL4+hGFClJkiSpnLV9lqmR6l1F8Eu1/74WeFMP7ydw4WpVJEmSJKmhKpXKKkNRZ2cnADNmzGh4Pf1JvasI1jvjJUmSJElrjXpnsIiI1wCHAbvRZZl24JLMXFK0OkmSJElqI3XNSEXEjsC9wPnATsDLtf+eD9wXEa8vXqEkSZIktYl6Z7AuAp4Cds/MB5dtjIgxwHXAN4F3lCtPkiRJktpHvddUTQBO7RquAGqvTwXeXKowSZIkSWo39QasucDAFbw3kOoy7pIkSZK0Vqr3FMGTgC9FxAOZOXPZxoh4K3AG8KmSxUmSJElqHxOvnNirdrMfm11Xe4Dpk6f3qaZmqzdgfRYYBtwaEY8CjwKb1R6PA5+JiM8sa5yZu5UqVJIkSZLWdPUGrDtqD0mSJEnt5os79L7t/Pn1f2bMqPrq6YfqvdHw4Y0qRJIkSZLaXb2LXEiSJEmSVsCAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUSL33wZIkSZKkFXrgqgeYe/XcXrW9+cCbV9lm7AfH0nFAx2pW1TwGLEmSJEnFdBzQ0VaBqDQDliRJkrSWqdy0kNN/9USv2san7lllm9P2Gk5l0ojVLatfMGBJkiRJa5nKpBEGogZxkQtJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSVqJSqVCRBR7VCqVVn8lSZLUQANaXYAkrckqlcoqQ1FnZycAM2bMaHg9kiRpzeYMliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLUttzIQpJkrSmcJELSW3PhSgkSdKawhksSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFeIiF5La1xd36H3b+fPr+8wn766/HkmStNZzBkuSJEmSCjFgSZIkSVIhniIoST2YeOXEXred/djsuj4zffL0PtUkSZLWfM5gSZIkSVIhBixJkiRJKsSAJUmSJK1Bdps6jd2mTmNox3giothjaMf4Vn+1tYLXYElqe5WbFnL6r57oVdv41D2rbHPaXsPhY6NWtyxJklbLjh/7UgN6PaYBfaorA5aktleZNILKpBFF++z9EheSJEl/4ymCkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEK80bAkrcQDVz3A3Kvn9qrtzQfevMo2Yz84FiavXk2SJGnNZcCSpJXoOKCDjgM6Wl2GJElqE54iKEmSJEmFND1gRcTwiPhJRDwTEfMi4pAVtNspIn4ZEQsjIvvajyRJkiQ1SytmsM4HlgAjgSnAhRExrod2LwJXAh9dzX4kSZIkqSmaGrAiYjCwP3BKZi7OzN8APwMO7d42M+/OzO8Ac1anH0mSJElqlsh81dl3jfthEbsAt2bmoC7bPgnskZl7r+Az2wL3Zmb0tZ+IOBI4EmDkyJG7XnHFFaW+kqRVuOvhRQ3r+3XxYMP6vuc16zWs7+033r5hfav5lu3jXz7j09x756v+Jthn2+04jotOO7pYfz1p1/3c3yuv1ujfK475q7mP96xdx7wv9txzz99l5oTu25sdsHYHrsrMzbtsOwKYkpmdK/hMTwGr7n6WmTBhQs6aNWu1voek3ttt6rSG9X3b+sc0rO+JY0Y1rO/pk6c3rG81X7vu49C++3m7jnm7jjc45j1xH+9Zu455X0REjwGr2ddgLQaGdds2DKg3ppfqR5IkSZKKaXbAugcYEBHbddk2nh6us2pSP5IkSZJUTFMDVmY+A/wYOCMiBkfE24F9gUu7t42qgcBraq8HRsT69fYjSZIkSc3SimXaPw4MAh4FLgeOzsw5ETEmIhZHxJhau62A5/jbrNRzwN2r6qcZX0CSJEmSejKg2T8wM58A9uth+4PAkC6v5wLRvd2q+pEkSZKkVmnFDJYkSZIk9UsGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkrVEqlQoRUexRqVRa/ZUkSWuRpi/TLknSylQqlVWGos7OTgBmzJjR8HokSaqHM1iSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYW4iqAkaY0w8cqJvW47+7HZdX1m+uTpfapJkqR6OYMlSZIkSYUYsCRJkiSpEANWP1epVIiIYo9V3fxTkiRJWpt5DVY/V6lUVhmKOjs7AZgxY0bD65EkSZL6M2ewpMKcNZQkSVp7OYMlFeasoSRJ0trLgCVJWqM8cNUDzL16bq/a3nzgzatsM/aDY2Hy6tUkSVJvGbAkSWuUjgM66Digo9VlSJLUJ16DpdXi9UaSJEnS3ziD1Z99cYfetZs/v772AJ+8G/B6o2arVCqcfvrpxfo77bTTDLWSJEkFGbCkNmKglSRJWrN5iqAkSZIkFWLAkiRJkqRCPEVQKqkJ171Jrea1gJIkrZgBS30y8cqJvW47+7HZdX1m+uTpfaqpv+vt+NU73uCYqz5eCyhJ0op5iqAkSZIkFeIMVj9XuWkhp//qiV61jU/ds8o2p+01nMqkEatbliRJktQvGbD6ucqkEQYiSZIkqUkMWFIbeeCqB5h79dxetb35wJtX2WbsB8fScUDHalYlSZKkZQxYUhvpOKDDQKTWauRKmWNG1V+PJElrGAOWVksjZlSYvHo1SZIkSa1iwNJqcUZFkiRJ+hsDllSYKzdKkiStvQxYUmGu3ChJkrT2MmBJkuriLK0kSStmwJIk1cVZWkmSVmydVhcgSZIkSf2FAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxICltcpuU6ex29RpDO0YT0QUewztGM9uU6e1+utJkiSpxQa0ugCpFXb82JdaXYIkSZL6IWewJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQpoesCJieET8JCKeiYh5EXHIStr+a0Q8HBFPRcR3I2L9Lu/NiIjnI2Jx7XF3c76BJEmSJPWsFTNY5wNLgJHAFODCiBjXvVFEvBs4CXgnMBbYGji9W7NjM3NI7bFDQ6uWJEmSpFVoasCKiMHA/sApmbk4M38D/Aw4tIfmHwa+k5lzMvOvwOeAw5pWrCRJkiTVqdkzWNsDL2fmPV22zQZeNYNV2za7W7uREbFJl21nRcTCiPhtRHSWLlaSJEmS6hGZ2bwfFrE7cFVmbt5l2xHAlMzs7Nb2fuCYzLyx9no9qqcWdmTm3Ih4C/Cn2raDgG8Ab8zM+3v4uUcCRwKMHDly1yuuuKIRX2/N88ichnV9z2vWa1jf22+8fcP6vuvhRQ3rG+B18WDD+nbMX83x7tldDy/iy2d8mnvvLPc7YLsdx3HCqWc55j1o130cHPOeuI/3zDF/NffxnrXrmPfFnnvu+bvMnNB9e7MD1i7AbzNzgy7b/g3ozMy9u7WdDUzNzCtrrzcBFgIjMvPxHvq+EbguM7++shomTJiQs2bNWv0v0w6+2LjL0iaOGdWwvqdPnt6wvnebOq1hfQPctv4xDevbMX81x7tnjvmruY/3zDF/Nffxnjnmr+Y+3rN2HfO+iIgeA1azTxG8BxgQEdt12TYe6OnPrHNq73Vt90hP4aomgShSpSRJkiT1QVMDVmY+A/wYOCMiBkfE24F9gUt7aP594KMR8fqI2Bj4LPA9gIjYKCLeHREDI2JAREwB3gH8silfRJIkSZJ60Ipl2j8ODAIeBS4Hjs7MORExpnY/qzEAtWuvzgFuBubVHqfV+lgPOBN4jOppg/8C7JeZ3gtLkiRJUssMaPYPzMwngP162P4gMKTbti8DX+6h7WPAmxtUoiRJkiT1SStmsCRJkiSpXzJgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSA1WK7TZ3G0I7xRESxx9CO8ew2dVqrv5okSZK01hnQ6gIEO37sS60uQZIkSVIBzmBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKmQpgesiBgeET+JiGciYl5EHLKStv8aEQ9HxFMR8d2IWL8v/UiSJElSM7RiBut8YAkwEpgCXBgR47o3ioh3AycB7wTGAlsDp9fbjyRJkiQ1S1MDVkQMBvYHTsnMxZn5G+BnwKE9NP8w8J3MnJOZfwU+BxzWh34kSZIkqSkiM5v3wyJ2AW7NzEFdtn0S2CMz9+7Wdjbw+cz8Ue31COAxYAQwprf91N47Ejiy9nIH4O6iX6x/GAEsbHURaxnHvPkc8+ZyvJvPMW8+x7y5HO/mc8xXbKvM3LT7xgFNLmII8FS3bU8BQ3vRdtnzoXX2Q2ZeBFxUb7Frk4iYlZkTWl3H2sQxbz7HvLkc7+ZzzJvPMW8ux7v5HPP6NfsarMXAsG7bhgGLetF22fNFdfYjSZIkSU3R7IB1DzAgIrbrsm08MKeHtnNq73Vt90hmPl5nP5IkSZLUFE0NWJn5DPBj4IyIGBwRbwf2BS7tofn3gY9GxOsjYmPgs8D3+tCPesdTKJvPMW8+x7y5HO/mc8ybzzFvLse7+RzzOjV1kQuo3r8K+C6wF/A4cFJm/jAixgB/Al6fmQ/W2p4AnAgMAq4BjsrMF1bWT1O/jCRJkiR10fSAJUmSJEn9VStuNCxJkiRJ/ZIBS5IkSZIKMWCpRxERra5BaiT3cUmS1AgGLPUoaxfneRDaHI5z87mPN59j3RqOu/o793GtaVzkQq8QEf8ODAceAH6bmXdExDqZubTFpfVLEfGmzPx97Xmk/yAbzn28+dzPmy8ihgGLuvwhwX28wSJiQGa+1Oo61hbu463hft47zmBpuYj4KfBBYBjwTuDnEfGezFwaEe4rhUXED4FZEXEYVGdU/CtcY7mPN5/7efNFxH8A1wLXRMTnANzHGyci/i0iRmTmSxGxbqvrWRu4jzef+3l93BEFQESMB7bMzN0y8+PAx6neWO66iHhf7ReXB0WFRMQxwE7AqcBnI+Ij4MFnI7mPN5/7efNFxIXAG4DPALcCfx8RvwUPQBshIi4GzqX6x5pNM/NlDz4by328+dzP6+dOqGX+CgyOiD0AMnNhZp4FnARcHBFv8bSeMiLiNcCDwOeAC6j+0jrVg8+Gcx9vIvfz5ouIQcDWwGcz8/8BXwE+VH0rboPlB6COewERsRfQQXVW/D7gFx58Npb7ePO5n/eNAUvLvATcD+zW7R/M+cDlwP4Rsa6/tFZfZi4BfgNcn5lPAJcBZ1E9+PxorU1GxKYtLLM/ep7qPv5m9/HGq+3nM4Dr3M+bIzOfAxYCfxcR62Xmy5k5H9gPeD4ivl1r5x8SypgLfB74FXAK1es6ezz49PdKGe7jLfFn3M/rZsBai0XEhIh4U0QMyswFVE+XOg04bFmbzHwWmAe8rvaLzF9afdRtvP+amc8AZOZi4AfA2cApEbFf7QD0BxEx0F9YfVebRQEgMx8FvkN1H/9wl+3u4wV1G/NFtfF1P2+ePwL/AGyzbENt3z8X6IiIka0qrL/JzHuB/6zt53OBk6kefF7X5eCzs9bW3yvlzAbei/t4U2Tm/cAt7uf1GdDqAtQaEfETYFvgBWDTiDg6M38aER8DvhcRg4EbM/Oe2kcyIgZm5vOtqrmddRvvTSLiaOC/M/MpgMx8JiJ+ADwO/JjqbMs7HO++q52nPycivlP7qyeZ+eOorjz1zdo+/kv38XJ6GvOu3M/Lq13n9iKwNDP/IzO/EBFvo3rx//uBubWDnpnACGAo8EjrKm5vXcb75cz8Ttd9NzPvj4jPAmcCV0XEb4CPRXUVzQUtKrntRcSWwEPLDt4z85zaqd5XRcQ+uI8X18OYv7DsPffz3jFgrYUi4liq59PuBiTVv+h/PiLGZuYFEbGE6nUph0bE48BbgT09COqbFYz3VODbEXFl7fQpMnNxRLwReBr4u8z8U4tKbnsR8TWq5+U/ATwXET/sErK+FxHPAf8GfMh9vIyVjXlX7uflRMTPga2A/wb+sRaoTsrM/SLiF8A1wDkRcTswEdgAeKplBbe5Hsb7PVRnTW7vciB6X+06w9nAscA7Pejsu4i4mupM1Yepzs4CkJnvi4ibgKuAL7qPl7OiMe/K/XzVDFhrp2HArV0Ofj4dEQ8Dh0XEwsy8MiLuBDYHNgOOzswHWlVsP7DC8aZ6LvnVtXOYxwAfAN7tQWff1f7ythXVZdhfR3WRhXUi4gddQtaPIuL3wJa4j6+23ox5l7YdwD/ifr5aImICsEVm7lx7/e/AFcDZEXFqZr4/Ir5E9ffMp6n+//0hmflYq2puZysZ708B51G9rnaZT1A9QN05M+9odq39RUScQfX/F58CvhoRx2fmHyOq99LLzEnu42WtaMxX0Nz9fCW80fBaKCIOoHodxD6ZOafL9lOBo6j+Y1nYqvr6m1WM95HA+Mx8PCLWAwZl5tMtKrXfiIhtgKcyc2HtOp8KcAbww2XXvqms3o557Rqtge7nq6d2zcNlwPbLxjciNgYuBsjM/WrbNgPWB57NzMdbUWt/sIrxfhH4UGY+FxFDqIauX2Tm7S0qt1+onY75dGZeWps9HAb8Sy1kLb+pcESMAAbhPr7aVjbm3dptQPVMp5+7n/fMgLWWiIh/BjYFZgHPUP0L8rrA1zLzvi7t/hu4ODMvakmh/YTj3Xzdxvyu2sW4y95bfsCfmd+OiMOB/5eZd7ai1v6izjE/DJjpmJcR1Xv9TANuoTrGWdu+IdXTei7OzErrKuxfejHe/5GZn6ttWzczX25Zsf1IRAxfdhp9lwP+TwCzXVChMVY15stmEN3PV86AtRaIiJ8Co4E/Ub1/xJ+BpcBzVK+D+O6yg56I+DFwTWZe1qJy257j3XxdxnwOsB0wH7g0M3/Wpc3HgBOBO4H3A2/oOqOo+jjmzRcRHwL+KzMfqM14/zPV605+BlzS5aD/X4Fxmfmx1lXb/hzv5us65l22DcjMl2rPf0F1EYsDqf7vMSgzT2pJsf1EH8Z8cGb+e0uKbSNeg9XPRcTbqa5et3NWl9KcAOxN9aL+p4AlwLdq/4AGAntSnfZVHzjezdfDmL+J6g0R/732F7afAGTmf0R1xam3A7t4oN93jnnzRXX1xb2BH0XEF7K6ktflwFiq4XUT4Eu15lsA63c9jUr1cbybr9uYn52Zf67Nlry0bLakdm3hFVRnDTcG/r6lRbc5x7xxvA9W//c01QP7EbVf/rOA7wP/CQTwP1Qv1H0n1b9C75l/W7Za9XO8m6/7mP+e6j3dpgEfjYhdASLiSKoHRhMzc3bLqu0fHPMmiogPUl106HPAesCJEbFd7XqTs6meovnuiLg3Ir5D9Vracz3Y7xvHu/l6GPOTImLrLqekvRwRyyYFZtXavDEzZ7ao5LbnmDeWpwj2cxExCrid6i//r3TZvjPwGar3Abq4dvrDy/4fxOpxvJtvJWM+juo1QNfXxvzNwPOZ+b+tqbT/cMybKyK2oDoL/jPgPcBBVE85Picz761dcD6E6ik8i6muWnp3q+ptd453861gzJ8Fus6qZFTv8fafwNtqf8BUHznmjWXAWgtExL7ApcAxtZVhlv2jOQsYD7zfA/1yHO/mW8WY7wzsu+x8cpXhmDdXRKyXmS/Wnu8DTKZ60P+FrN6T5vXAfZm5pJV19heOd/OtZMzPqh3w75SZd0TEyMz0RsIFOOaN4zVYa4efUb3Q/CsRMTgzv1nbvgj4P6qnrqkcx7v5VjbmD1O9wbPKcsybKDNfXBZiM/NnERHAAcDRUb2P3h7AXlTvrafV5Hg330rG/JjamE+MiE4P9MtxzBvHGay1RO0fyoHAt4Fbqf6F4h1UrwH6Qytr648c7+ZzzJvPMW++ZQdDtedvA74HvJbqdW6evlOY4918jnnzOeblGbDWMhGxPdVzbl8D3JKZ97a4pH7N8W4+x7z5HPPm6nI65vFUV7Ibn5l3tLisfsvxbj7HvPkc87I8RXAtk9UV61y1rkkc7+ZzzJvPMW+u2kHQEOANwG4eBDWW4918jnnzOeZlOYMlSVIb6nqBuhrP8W4+x7z5HPMyDFiSJEmSVIg3GpYkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqZD/D/dkxGjUZ40MAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 864x576 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "s = 2\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAI4CAYAAAB3HEhGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA/EklEQVR4nO3deZhcVZ3/8fcXgiRkAUIgSkhIWEQMEpGIOg7SRok6yuIgYcmgosKAMMigCIiQAokg4IIKKI4iIrIIrmxihGZGmV8g6uAQ2SUhkmEJCCRsAfL9/VGV2CSdpKtzqirV/X49Tz1U3Tp1+luHS3M/fe49NzITSZIkSdKaW6fVBUiSJElSX2HAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEDWl1As40YMSLHjh3b6jIkSZIktbHf//73CzJz0+W397uANXbsWGbNmtXqMiRJkiS1sYiY2912TxGUJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIf1ukQtJkiSpp55++mkeffRRXnzxxVaXoiZbb7312GyzzRg2bFhdnzNgSZIkSd14+umneeSRRxg1ahSDBg0iIlpdkpokM3nuued46KGHAOoKWZ4iKEmSJHXj0UcfZdSoUWywwQaGq34mIthggw0YNWoUjz76aF2fNWBJkiRJ3XjxxRcZNGhQq8tQCw0aNKju00MNWJIkSdJKOHPVv/Xm378BS5IkSZIKMWBJkiRJbWj8+PF0dnausL2zs5MtttiiyM+ICO67774iffUXriIoSZIk9dAu02c0tP9bT3x3j9vOnj27gZWot5zBkiRJkqRCDFiSJElSGxo7diwzZszgueee46Mf/Sgbb7wxr3/967ntttte0W7+/Pnss88+bLrppowbN46vf/3ry9679dZbedvb3sZGG23Ea17zGo488kgWL17c7K/SpxiwJEmSpDZ2yimncP/993P//ffzq1/9iosuumjZe0uWLGGPPfZgwoQJPPTQQ/zmN7/ha1/7Gr/61a8AWHfddfnqV7/KggUL+O///m9+85vfcN5557Xqq/QJBixJkiSpjV1xxRWceOKJDB8+nNGjR3PUUUcte++2227jscce4+STT+ZVr3oVW221FYcccgiXXXYZADvvvDNvfetbGTBgAGPHjuVf//Vfufnmm1v1VfoEF7mQJEmS2tj8+fMZPXr0stdbbrnlsudz585l/vz5bLTRRsu2vfzyy+y6664A3HPPPRxzzDHMmjWLZ599lpdeeomdd965abX3Rc5gSZIkSW3sNa95DfPmzVv2+sEHH1z2fPTo0YwbN44nn3xy2WPhwoVce+21ABx++OG87nWv49577+Xpp5/mi1/8IpnZ9O/QlxiwJEmSpDY2ZcoUTj/9dP72t7/x17/+lW984xvL3ttll10YNmwYX/rSl3juued4+eWXueOOO5YthLFw4UKGDRvGkCFDuOuuuzj//PNb9TX6DE8RlCRJknqonvtUNcu0adM47LDDGDduHJtvvjkHH3ww55xzDlBdxOKXv/wln/70pxk3bhwvvPAC2223HaeddhoAZ599NoceeihnnnkmO+20E/vttx833nhjK79O24tmTwFGxHDgu8BkYAFwQmb+qJt2+wOnAK8GXgCuA/4tM5+up5/lTZw4MWfNmlXo20iSJKmvuvPOO9l+++1bXYZabGX7QUT8PjMnLr+9FacIngssBkYCU4HzI2J8N+1+B7w9MzcEtqI623ZaL/qRJEmSpKZoasCKiMHAPsBJmbkoM38L/AI4aPm2mTkvMxd02fQysE29/UiSJElSszT7GqzXAi9n5j1dtt0O7NZd44j4R+AaYBjwLPDBXvZzKHAowMiRI+ns7FyDryBJklrprocXAvCVU0/g3jtnF+t32+3Hc8F53yzWn9rfhhtuyMKFC1tdhlrs+eefrys/NDtgDQGeWm7bU8DQ7hrXZqY2jIhRwCHAnF72cwFwAVSvwero6OhF6ZIkaW3w2ekzANho6td4c+G+PUZQV3feeSdDh3Z7eKl+ZODAgey00049bt/sa7AWUZ2N6moYsMo/DWTmQ8D1wGVr0o8kSZIkNVKzA9Y9wICI2LbLtglAT+b3BwBbF+hHkiRJkhqiqQErM58BfgKcGhGDI+LtwF7Axcu3jYipETEmqrYEpgO/qbcfSZIkSWqWVizT/klgEPAocClweGbOroWpRRExptbu9cAtVE8H/B1wN9XrsFbZT5O+gyRJkiStoNmLXJCZTwB7d7P9QaqLVyx9fSJwYr39SJIkSVKrND1gSZIkSW3r7O0a2/9n7u5x07vvvpv999+f++67j+nTp3PUUUd1227OnDmMGzeOF198kQEDPPzvqhFj4whLkiRJbejMM8+ko6ODP/7xj60upWEqlQr33XcfP/zhD1tdSo+14hosSZIkSWto7ty5jB8/vtVltFRmsmTJklaX8QoGLEmSJKnNTJo0iZtuuokjjzySIUOGcM4557DTTjsxbNgwRo8eTaVSWelnr7rqKsaOHcsdd9zBkiVLOOOMM9h6663ZZJNNmDJlCk888cQqf/acOXOICC688EJGjx7NxhtvzLe+9S1uu+02dtxxRzbaaCOOPPLIZe2XLFnCaaedxpZbbslmm23Ghz/8YZ566qlX9HXRRRcxZswYRowYwfTp0wG4/vrr+eIXv8jll1/OkCFDmDBhAlC9IfiJJ57I29/+djbYYAP+8pe/cNddd7H77rszfPhwtttuO6644oplP/+aa67p8diUYMCSJEmS2syNN97Irrvuyje/+U0WLVrEhAkT+MEPfsCTTz7JNddcw/nnn8/PfvazFT534YUXctxxxzFjxgx22GEHvv71r/Ozn/2Mm2++mfnz57PxxhtzxBFH9KiGmTNncu+993L55Zdz9NFHM336dGbMmMHs2bO54ooruPnmmwH4/ve/z/e//31uuukm/vKXv7Bo0aJXBDCA3/72t9x999385je/4dRTT+XOO+/kve99L5/73OfYb7/9WLRoEbfffvuy9hdffDEXXHABCxcuZNNNN2X33XfnwAMP5NFHH+XSSy/lk5/8JLNnVxcYHzx4cI/GphQDliRJktTmOjo6eMMb3sA666zDjjvuyAEHHLAs4Cz1ta99jbPOOovOzk622WYbAL797W8zffp0tthiC9Zff30qlQpXXnklL7300mp/5kknncTAgQOZPHkygwcP5oADDmCzzTZj1KhR7LrrrsuuDbvkkks45phj2GqrrRgyZAinn346l1122St+xrRp0xg0aBATJkxgwoQJrwhT3fnoRz/K+PHjGTBgANdffz1jx47l4IMPZsCAAbzpTW9in3324corr+zx2JTkIheSJElSm5s5cybHH388d9xxB4sXL+aFF15g3333fUWbs846i5NPPpkttthi2ba5c+fywQ9+kHXW+fu8y7rrrssjjzzCqFGjVvkzR44cuez5oEGDVni9aNEiAObPn8+WW2657L0tt9ySl156iUceeWTZtle/+tXLnm+wwQbLPrsyo0ePfsV3mDlzJhtttNGybS+99BIHHXQQ0LOxKckZLEmSJKnNHXjggey5557MmzePp556isMOO4zMfEWbG264gdNOO42rrrpq2bbRo0dz3XXX8eSTTy57PP/886sNV/XYfPPNmTt37rLXDz74IAMGDHhFIFuZiFjt9tGjR7Pbbru94jssWrSI888/H+jZ2JTkDJYkSZLUU3Xcp6qZFi5cyPDhwxk4cCC33norP/rRj5g8efIr2owfP57rr7+e97znPay33nrsueeeHHbYYZx44olcdNFFbLnlljz22GPccsst7LXXXsVqO+CAA/jSl77E+973PjbddNNl11X15L5TI0eO5Ne//jVLlix5xSxbVx/4wAc4/vjjufjii9l///0B+J//+R+GDBnC9ttv36OxKckZLEmSJKnNnXfeeZx88skMHTqUU089lSlTpnTbbsKECVx99dUccsghXHfddXzqU59izz33ZPLkyQwdOpS3vvWtzJw5s2htH/vYxzjooIN4xzvewbhx4xg4cCDf+MY3evTZpafybbLJJrzpTW/qts3QoUO54YYbuOyyy9h888159atfzXHHHccLL7wA9HxsSolGTo+tjSZOnJizZs1qdRmSJKmXdpk+o2F933riuxvWt9rPnXfeyfbbb9/qMtRiK9sPIuL3mTlx+e3OYEmSJElSIQYsSZIkSa9wySWXMGTIkBUe48ePb3Vpaz0XuZAkSZL0ClOnTmXq1KmtLqMtOYMlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIqwhKkiRJPTTpikkN7f/GKTf2uO3dd9/N/vvvz3333cf06dM56qijum03Z84cxo0bx4svvsiAAe1/+B8R3HvvvWyzzTZ1f3b8+PGce+65dHR0lC+spv1HWJIkSeqHzjzzTDo6OvjjH//Y6lLaxuzZsxv+MzxFUJIkSWpDc+fO9ca/ayEDliRJktRmJk2axE033cSRRx7JkCFDOOecc9hpp50YNmwYo0ePplKprPSzV111FWPHjuWOO+5gyZIlnHHGGWy99dZssskmTJkyhSeeeGKVP3vOnDlEBBdeeCGjR49m44035lvf+ha33XYbO+64IxtttBFHHnnksvaVSoV/+Zd/WeHzL730EgAdHR18/vOf5x/+4R8YMmQIe+yxB48//jhTp05l2LBhvPnNb2bOnDmvqOHaa69lq622YsSIERx77LEsWbIEgPvvv59JkyaxySabMGLECKZOncqTTz657HNjx45lxowZPRzl3jFgSZIkSW3mxhtvZNddd+Wb3/wmixYtYsKECfzgBz/gySef5JprruH888/nZz/72Qqfu/DCCznuuOOYMWMGO+ywA1//+tf52c9+xs0338z8+fPZeOONOeKII3pUw8yZM7n33nu5/PLLOfroo5k+fTozZsxg9uzZXHHFFdx88809/j6XXXYZF198MQ899BD3338/b3vb2zj44IN54okn2H777TnllFNe0f6nP/0ps2bN4g9/+AM///nP+d73vgdAZnLCCScwf/587rzzTubNm7fKsNkIBixJkiSpzXV0dPCGN7yBddZZhx133JEDDjhghYDzta99jbPOOovOzs5lC0R8+9vfZvr06WyxxRasv/76VCoVrrzyymWzS6ty0kknMXDgQCZPnszgwYM54IAD2GyzzRg1ahS77rprXdeGHXzwwWy99dZsuOGGvO9972Prrbfm3e9+NwMGDGDfffddoa/jjjuO4cOHM2bMGI4++mguvfRSALbZZht233131l9/fTbddFOOOeaYuoJeCS5yIUmSJLW5mTNncvzxx3PHHXewePFiXnjhBfbdd99XtDnrrLM4+eST2WKLLZZtmzt3Lh/84AdZZ52/z7usu+66PPLII4waNWqVP3PkyJHLng8aNGiF14sWLepx/fX2NXr06GXPt9xyS+bPnw/Ao48+ylFHHcV//dd/sXDhQpYsWcLGG2/c4zpKcAZLkiRJanMHHngge+65J/PmzeOpp57isMMOIzNf0eaGG27gtNNO46qrrlq2bfTo0Vx33XU8+eSTyx7PP//8asNVPQYPHsyzzz677PXDDz+8xn3Omzdv2fMHH3yQzTffHIATTjiBiOBPf/oTTz/9ND/84Q9XGIdGcwZLkiRJ6qF67lPVTAsXLmT48OEMHDiQW2+9lR/96EdMnjz5FW3Gjx/P9ddfz3ve8x7WW2899txzTw477DBOPPFELrroIrbccksee+wxbrnlFvbaa69itb3xjW/kS1/6Eg8++CAbbrghp59++hr3edZZZ/GWt7yFRYsWcc4553DMMccA1XHYcMMN2WijjXjooYc466yz1vhn1csZLEmSJKnNnXfeeZx88skMHTqUU089lSlTpnTbbsKECVx99dUccsghXHfddXzqU59izz33ZPLkyQwdOpS3vvWtzJw5s2htu+++O/vttx877rgjO++8Mx/4wAfWuM+99tqLnXfemTe+8Y28//3v5+Mf/zgA06ZN4w9/+AMbbrgh73//+/nnf/7nNf5Z9YpmT5m12sSJE3PWrFmtLkOSJPXSLtMbt8TyrSe+u2F9q/3ceeedbL/99q0uQy22sv0gIn6fmROX3+4MliRJkiQVYsCSJEmS9AqXXHIJQ4YMWeExfvz4Vpe21nORC0mSJEmvMHXqVKZOndrqMtqSM1iSJEnSSixZsqTVJaiFevPv34AlSZIkdWPw4ME89NBDLF68uOn3UlJrZSaLFy/moYceYvDgwXV91lMEJUmSpG5sscUWLFiwgLlz5/LSSy+1uhw12YABA9hwww0ZMWJEfZ9rUD2SJElSW1tnnXXYbLPN2GyzzVpditqIpwhKkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqZCmB6yIGB4RP42IZyJibkQcuJJ2H4mI30fE0xHx14g4MyIGdHm/MyKej4hFtcfdzfsWkiRJkrSiVsxgnQssBkYCU4HzI2J8N+02AI4GRgBvAd4FfGa5Nkdm5pDaY7vGlSxJkiRJqzdg9U3KiYjBwD7ADpm5CPhtRPwCOAg4vmvbzDy/y8uHIuIS4J1NK1aSJEmS6tTUgAW8Fng5M+/psu12YLcefPYdwOzltp0eEWcAdwMnZmZndx+MiEOBQwFGjhxJZ2e3zSRJUhv42NbPNaxvjxEkralmB6whwFPLbXsKGLqqD0XEwcBE4BNdNh8H/Jnq6Yb7A7+MiDdm5v3Lfz4zLwAuAJg4cWJ2dHT0tn5JktRin50+o2F937p/R8P6ltQ/NPsarEXAsOW2DQMWruwDEbE3cAbwvsxcsHR7Zs7MzIWZ+UJmXgT8Dvin8iVLkiRJUs80O2DdAwyIiG27bJvAiqf+ARAR7wW+A+yRmf+7mr4TiCJVSpIkSVIvNDVgZeYzwE+AUyNicES8HdgLuHj5thExCbgE2Cczb13uvY0i4j0RMTAiBkTEVKrXaP2q8d9CkiRJkrrXimXaPwkMAh4FLgUOz8zZETGmdj+rMbV2JwEbAtd2udfVdbX31gNOAx4DFgD/Buydmd4LS5IkSVLLNHuRCzLzCWDvbrY/SHURjKWvV7oke2Y+Bry5EfVJkiRJUm+1YgZLkiRJkvokA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJGk5lUqFiCj2qFQqrf5KkppkQKsLkCRJWttUKpXVhqKOjg4AOjs7G16PpPbhDJYkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBUyoNUFSJIkrTXO3q7nbefNq+8zn7m7/noktR1nsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqZABrS5AkiRpbVO5YQGn/PqJHrWNY+9ZbZtpuw+n8pk1rUpSOzBgSZIkLacyeQSVySNaXYakNuQpgpIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiSpiSqVChFR7FGpVFr9lSR1MaDVBUiSJPUnlUpltaGoo6MDgM7OzobXI6ksZ7AkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrE+2BJkiQ1waQrJvW47e2P3V7XZ26ccmOvapJUnjNYkiRJklSIAUuSJEmSCjFgSZIkSVIhTQ9YETE8In4aEc9ExNyIOHAl7T4SEb+PiKcj4q8RcWZEDKi3H0mSJElqllbMYJ0LLAZGAlOB8yNifDftNgCOBkYAbwHeBXymF/1IkiRJUlM0NWBFxGBgH+CkzFyUmb8FfgEctHzbzDw/M/8rMxdn5kPAJcDb6+1HkiRJkpql2cu0vxZ4OTPv6bLtdmC3Hnz2HcDs3vQTEYcChwKMHDmSzs7OOsuWJElri49t/VzD+u6MoxrW9/7rrdfjtg/FQ7XP7N+j9o0+trnr4YV85dQTuPfO2atv3EPbbj+eY04+nde9emixPqW1QbMD1hDgqeW2PQWs8r+siDgYmAh8ojf9ZOYFwAUAEydOzI6OjrqKliRJa4/PTp/RsL5vXf/rDev71DGjetz20XwUgMtevKxH7W/saOx9sD47fQYbTf0aby7c7/fuh1v37yjcq9RazQ5Yi4Bhy20bBixc2QciYm/gDODdmbmgt/1IkiRJUqM1O2DdAwyIiG0z897atgn8/dS/V4iI9wLfAd6fmf/b234kSZLWFg/8+AHmXDmnR21v2u+m1bYZ+6GxMGXNapJUTlMDVmY+ExE/AU6NiE8AbwT2Av5h+bYRMYnqwhYfzMxbe9uPJEnS2mTcvuMYt++4VpchqUFasUz7J4FBwKPApcDhmTk7IsZExKKIGFNrdxKwIXBtbfuiiLhudf0072tIkiRJ0is1+xRBMvMJYO9utj9IdfGKpa/f2Zt+JEmSJKlVWjGDJUmSJEl9kgFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIXUFrIjYoVGFSJIkSVK7q3cG608RcVtEHB4RGzWiIEmSJElqV/UGrHcBfwbOBOZHxKURsXtERPnSJEmSJKm91BWwMvOmzPwI8GrgSGAU8CtgbkR8ISK2bkCNkiRJktQWerXIRWY+k5nfy8x3AK8F5gCfA+6JiJsj4oMFa5QkSZKkttDrVQQjYmxEVIAbgLcB1wKHAo8Al0fEV4tUKEmSJEltot5VBDeIiA9HxE3AfcBU4DvAmMzcIzO/m5lTgH8FPl6+XEmSJElaew2os/3DwLrAT4B3Z2bnStrdBjy+BnVJkiRJUtupN2AdD1ySmU+tqlFm3gGM63VVkiRJktSG6r0GawQwuLs3IuI1EXHympckSZIkSe2p3oA1DdhiJe9tXntfkiRJkvqlegNWALmS97YA/rZm5UiSJElS+1rtNVgR8RHgI7WXCZwfEU8v12wg8AaqS7ZLkiRJUr/Uk0UunuXvKwIG8BTwxHJtFgPXAeeVK02SJEmS2stqTxHMzB9n5r6ZuS9wEXDI0tddHlMz8wuZ6dLskiRJqkulUiEiij0qlUqrv5L6sbqWac/MgxtViCRJkvqnSqWy2lDU0dEBQGdnZ8PrkdZET67BOhP4emb+tfZ8VTIzjytTmiRJkiS1l57MYO0LXAL8tfZ8VRIwYEmSJEnql1YbsDJzXHfPJUmSpDV29nY9azdvXn3tAT5zd/31SGuo3vtgSZIkSZJWoifXYP1TPR1m5rW9L0eSJEmS2ldPrsG6muq1VdGDtgmsu0YVSZIkSVKb6knA8rorSZIkNUzlhgWc8usnetQ2jr1ntW2m7T6cyuQRa1qW1Cs9WeRibjMKkSRJUv9UmTzCQKQ+oyfXYG2Qmc8ufb669kvbSpIkSVJ/05NTBBdGxNsy81ZgEdXrrFbFa7AkSZIk9Us9CVgfA+7v8nx1AUuSJEmS+qWeXIN1UZfn329oNZIkSZLUxnoyg7WCiNgI2AF4DfB/wB2Z+WS5siRJkiSp/axTT+OIGBARXwL+CvwncHntn3+NiDMjYr0G1ChJkiT1SqVSISKKPSqVSqu/ktZy9c5gfQU4FDgV+AnwKLAZsA9wEjAQOKpkgZIkSVJvVSqV1Yaijo4OADo7Oxtej/q+egPWQcDnMvMrXbY9AUyPiOeBz2PAkiRJktRP1XWKILAEmL2S9+7AFQYlSZIk9WP1BqyLgU+s5L1DgB+uWTmSJEmS1L5We4pgRHyyy8s5wIciYjbwC/5+DdZewFDg7AbUKEmSJEltoSfXYH2zm22bA9t3s/0rwDlrVJEkSZIktame3Gi43tMIJUmSJKlfMjxJkiRJUiH1LtMOQERsAbyW6n2vXiEzr13ToiRJkiSpHdUVsCJiKHAFMHnppto/uy7Pvm6BuiRJkqQ1MumKST1qd/tjt9fVHuDGKTf2qib1ffWeIng6MAbYlWq4+iDQAXwXeAB4a8niJEmSJKmd1Buw/gmYDsysvZ6fmf+ZmYcCPweOLVmcJEmSJLWTegPWSGBeZr4MPAMM7/Letfz91EFJkiRJ6nfqDVjzgBG15/cCH+jy3luA50sUJUmSJEntqN5VBH8NvBv4KfBV4KKI2Bl4AXgH8OWy5UmSJElS+6g3YB0HbACQmRdHxCLgQ8Ag4Ejg22XLkyRJkqT2Udcpgpn5bGYu6PL6p5k5NTP/OTPPz8wlq+sjIoZHxE8j4pmImBsRB66k3Q4R8auIWBAR2c37nRHxfEQsqj3urue7SJIkSVJp9V6DBUBEbBcR/xIRx0bE1IjYro6PnwssprpgxlTg/IgY3027F6nec+vjq+jryMwcUnvUU4OkPqRSqRARxR6VSqXVX0mSJLWpem80PAz4DrAP1XC2CBgCLImInwCfyMynV/H5wbXP7pCZi4DfRsQvgIOA47u2zcy7gbsjYpt6apTU/1QqldWGoo6ODgA6OzsbXo8kSeq/6r0G6zyqS7F/GPhJZj4fEQOphqZv1t7/l1V8/rXAy5l5T5dttwO71VnHUqdHxBnA3cCJmdnZXaOIOBQ4FGDkyJEeYElNdNfDCwH4yqkncO+ds4v1u+3247ngvG/2uP2TTz4JGLCkvuBjWz/XsL4746iG9b3/eus1rO9G/25r5zH/5WW/5Oorru5R+5v2u2m1bT4w5QPssf8e/v9EKxWZK1zetPLGEQuBf8/M/+jmvUOAr2Tm0FV8flfgx5n56uU+NzUzO1bymW2AezMzltv+FuDPVE833J9qwHtjZt6/qu8wceLEnDVr1qqaSCpol+kzGtb3rSe+u8dtncGS+o6G/l5Z/4iG9T1pzKiG9X3jlBsb1jc45t1p9Jhr7RcRv8/Mictvr/carEXA/63kvflUbz68us8PW27bMGBhnXWQmTMzc2FmvpCZFwG/A/6p3n4kSZIkqZR6A9a5wGciYlDXjRGxAfAZqqcIrso9wICI2LbLtglAifOGEojVtpIkSZKkBlntNVgRceZym7YF5kXEr4FHgc2A3YHngFWee5eZz9QWwzg1Ij4BvBHYC/iHbn5uAOsDr6q9HljtIl+IiI2AtwA3Ay8B+1G90fHRq/s+kiRJktQoPVnkYt/lXr9Ye7y1y7alp/jtAxy7mv4+CXyPajh7HDg8M2dHxBiq11S9PjMfBLYEHujyueeAucBYYD3gNOB1wMvAXcDetZUHJUmSJKklVhuwMnNcyR+YmU8Ae3ez/UGqS74vfT2HlZzyl5mPAW8uWZckSZIkrale3WhYkiRJkrSiugNWRGwVEedHxP9GxEO1f54XEVs1okBJkiRJahd13Wg4InYGbgKeB64GHgFGUr32ampEvDMz/1C8SkmSJElqA3UFLOBs4I/A+zLz2aUba8u0X1t7f1K58iRJkiSpfdR7iuAuwJldwxVA7fXZVJdOlyRJkqR+qd4ZrOeATVby3nCqpw5KUnOcvV3P286bV99nPuNdH9Q/VCoVTjnllGL9TZs2jUqlUqw/SWo39Qasa4AzIuIvmfnbpRsj4h+B04FflixOkiQ1VqVSWW0g6ujoAKCzs7Ph9UhSu6s3YB0D/By4OSIeo7rIxWa1xy3Ap8uWJ0mSJEnto66AlZmPA/8YEe+leqPf1wD/B8zMzBsaUJ8kSZIktY0eB6yIGAj8CTgqM68Hrm9YVZIkSZLUhnq8imBmPg9sBCxpWDWStJapVCpERLGHF/9LktS31btM+yXAwY0oRJLWRpVKhcxc5WO33XZjt912W227zDRgSZLUx9UbsB4E3hERsyLi1Ig4IiI+2eVxeCOKlCSpP3MmVZLaR72rCH659s/XAG/q5v0Ezl+jiiRJ0iu4lLoktY96VxGsd8ZLkiRJkvqNemewiIhXAR8FdqHLMu3ARZm5uGh1kiRJktRG6pqRiojtgXuBc4EdgJdr/zwXuC8iXl+8QkmSJElqE/XOYF0APAXsmpkPLt0YEWOAa4BvAe8oV54kSZIktY96r6maCJzcNVwB1F6fDLy5VGGSJEmS1G7qncGaAwxcyXsDqS7jLkmS2sCkKyb1qN3tj91eV3uAG6fc2KuaJKnd1Ruwjge+HBEPZObMpRsj4q3AqcCxJYuTpJ6o3LCAU379RI/axrH3rLbNtN2HU/nMmlYlSVL/VKlUOOWUU4r1N23atLa6f1+9AevzwDDgloh4FHgU2Kz2eBz4XER8bmnjzNylVKGStDKVySOoTB7R6jIkSRLeu6/egHVH7SFJkiRJWk69Nxo+uFGFSJIkSVK7q3cVQUmSJEnSShiwJEmSJKmQeq/BkqR+oZ7lqOtdwtrlqyVJ6rucwZIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsRVBCVJ6sce+PEDzLlyTo/a3rTfTattM/ZDYxm377g1rEpSu+rpirr1rsAL7bMKrwFLkqR+bNy+4wxEklSQpwhKkiRJUiEGLEmSJEkqxFMEJUlaW529Xc/bzptX/2fGjKqvHknSajmDJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEJc5EKSVqERN2FlyprVJEmS1l4GLElaBW/CKkmS6uEpgpIkSVI/U6lUiIhij0ql0uqvtNZwBkuSJEnqZyqVympDUUdHBwCdnZ119d2I0+vb6WwSA5YkSZKkYvr76fWeIihJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIirCEqSJEn9xdnb9bztvHn1f2bMqPrq6YOcwZIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRDvgyVJ0lqucsMCTvn1Ez1qG8fes9o203YfTmXyiDUtS1KD7DJ9BgB3/senWTTnT8X6HTJ2RxYeUX3u75XGMWBJkrSWq0we4YGL1A9t/4kvN6DXasLy90rjeIqgJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgppesCKiOER8dOIeCYi5kbEgStpt0NE/CoiFkRE9rYfSZIkSWqWVsxgnQssBkYCU4HzI2J8N+1eBK4APr6G/UiSJElSUzQ1YEXEYGAf4KTMXJSZvwV+ARy0fNvMvDszvwvMXpN+JEmSJKlZInOFs+8a98MidgJuycxBXbZ9BtgtM/dYyWe2Ae7NzOhtPxFxKHAowMiRI3e+7LLLSn0lSatx18MLG9b36+LBhvV9z6vWa1jfr934tQ3rW823dB//yqkncO+dK/xNsNe23X48F0w7vFh/3WnX/dzfKytq9O8Vx3xF7uPda9cx7413vvOdv8/Mictvb3bA2hX4cWa+usu2Q4Cpmdmxks90F7Dq7mepiRMn5qxZs9boe0jquV2mz2hY37euf0TD+p40ZlTD+r5xyo0N61vN1677OLTvft6uY96u4w2OeXfcx7vXrmPeGxHRbcBq9jVYi4Bhy20bBtQb00v1I0mSJEnFNDtg3QMMiIhtu2ybQDfXWTWpH0mSJEkqpqkBKzOfAX4CnBoRgyPi7cBewMXLt42qgcCraq8HRsT69fYjSZIkSc3SimXaPwkMAh4FLgUOz8zZETEmIhZFxJhauy2B5/j7rNRzwN2r66cZX0CSJEmSujOg2T8wM58A9u5m+4PAkC6v5wCxfLvV9SNJkiRJrdKKGSxJkiRJ6pMMWJIkSZJUiAFLkiRJkgoxYGmNVCoVIqLYo1KptPorSZIkSb3W9EUu1LdUKpXVhqKOjg4AOjs7G16PJEmS1ErOYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5ZUmCsrSpIk9V+uIigV5sqKkiRJ/ZczWFIbcXZMkiRp7eYMltRGnB2TJElauxmw1CuTrpjU47a3P3Z7XZ+5ccqNvapJkiRJajVPEZQkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhLnIhtYmeLhJS76Ii4MIikiRJpTiDJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEJc5EIq6eztetZu3rz62gOMGVV/PVIbqlQqnHLKKcX6mzZtGpVKpVh/kiStigFLa+SBHz/AnCvn9KjtTfvdtNo2Yz80FqasWU2S2lulUlltIOro6ACgs7Oz4fVIklQPA5bWyLh9xzFu33GtLkOSJElaKxiwpDbSiBlDA7IkSVI5BiypjThjKEmStHZzFUFJkiRJKsSAJUmSJEmFGLAkSZIkqRCvwZIkrRUmXTGpx21vf+z2uj5z45Qbe1WTJEn1cgZLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQVxGUCqvcsIBTfv1Ej9rGsfests203YdTmTxiTcuSiqlUKpxyyinF+ps2bRqVSqVYf5IktZIBSyqsMnmEgUh9WqVSWW0g6ujoAKCzs7Pu/h/48QPMuXJOj9retN9Nq20z9kNjYUrdZUiS1CsGrD7OvzRLajfj9h3HuH3HtboMSZJ6xYDVxzX6L82SJEmS/s5FLiRJkiSpEAOWJEmSJBViwJIkSZKkQrwGS5LUc2dv17N28+bV1x5gzKj665EkaS3jDJYkSZIkFeIMVl/mX5olSZKkpnIGS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEFcRlCTVpXLDAk759RM9ahvH3rPaNtN2H05l8og1LUuSpLWCAUuSVJfK5BEGIkmSVsKA1cf5l2ZJkiSpeQxYfZx/aZYkSZKax0UuJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFjqV3aZPoNdps9g6LgJRESxx9BxE9hl+oxWfz1JkiS1mItcqF/a/hNfbnUJkiRJ6oOcwZIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgppesCKiOER8dOIeCYi5kbEgato++8R8XBEPBUR34uI9bu81xkRz0fEotrj7uZ8A0mSJEnqXitmsM4FFgMjganA+RExfvlGEfEe4HjgXcBYYCvglOWaHZmZQ2qP7RpatSRJkiStRlMDVkQMBvYBTsrMRZn5W+AXwEHdNP8I8N3MnJ2ZfwO+AHy0acVKkiRJUp2aPYP1WuDlzLyny7bbgRVmsGrbbl+u3ciI2KTLttMjYkFE/C4iOkoXK0mSJEn1iMxs3g+L2BX4cWa+usu2Q4CpmdmxXNv7gSMy8/ra6/Wonlo4LjPnRMRbgD/Xtu0PfBN4Y2be383PPRQ4FGDkyJE7X3bZZY34emufR2Y3rOt7XrVew/p+7cavbVjfdz28sGF9A7wuHmxY3475ihzv7t318EK+cuoJ3Htnud8B224/nmNOPt0x70a77uPgmHfHfbx7jvmK3Me7165j3hvvfOc7f5+ZE5ff3uyAtRPwu8zcoMu2TwMdmbnHcm1vB6Zn5hW115sAC4ARmfl4N31fD1yTmd9YVQ0TJ07MWbNmrfmXaQdnN+6ytEljRjWs7xun3NiwvneZPqNhfQPcuv4RDevbMV+R4909x3xF7uPdc8xX5D7ePcd8Re7j3WvXMe+NiOg2YDX7FMF7gAERsW2XbROA7v7MOrv2Xtd2j3QXrmoSiCJVSpIkSVIvNDVgZeYzwE+AUyNicES8HdgLuLib5j8APh4Rr4+IjYHPA98HiIiNIuI9ETEwIgZExFTgHcCvmvJFJEmSJKkbrVim/ZPAIOBR4FLg8MycHRFjavezGgNQu/bqTOAmYG7tMa3Wx3rAacBjVE8b/Ddg78z0XliSJEmSWmZAs39gZj4B7N3N9geBIctt+wrwlW7aPga8uUElSpIkSVKvtGIGS5IkSZL6JAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYLbbL9BkMHTeBiCj2GDpuArtMn9HqryZJkiT1OwNaXYBg+098udUlSJIkSSrAGSxJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBViwJIkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqRADliRJkiQVYsCSJEmSpEIMWJIkSZJUiAFLkiRJkgoxYEmSJElSIQYsSZIkSSrEgCVJkiRJhRiwJEmSJKkQA5YkSZIkFWLAkiRJkqRCDFiSJEmSVIgBS5IkSZIKMWBJkiRJUiEGLEmSJEkqxIAlSZIkSYUYsCRJkiSpEAOWJEmSJBXS9IAVEcMj4qcR8UxEzI2IA1fR9t8j4uGIeCoivhcR6/emH0mSJElqhlbMYJ0LLAZGAlOB8yNi/PKNIuI9wPHAu4CxwFbAKfX2I0mSJEnN0tSAFRGDgX2AkzJzUWb+FvgFcFA3zT8CfDczZ2fm34AvAB/tRT+SJEmS1BSRmc37YRE7Abdk5qAu2z4D7JaZeyzX9nbgi5l5ee31COAxYAQwpqf91N47FDi09nI74O6iX6xvGAEsaHUR/Yxj3nyOeXM53s3nmDefY95cjnfzOeYrt2Vmbrr8xgFNLmII8NRy254Chvag7dLnQ+vsh8y8ALig3mL7k4iYlZkTW11Hf+KYN59j3lyOd/M55s3nmDeX4918jnn9mn0N1iJg2HLbhgELe9B26fOFdfYjSZIkSU3R7IB1DzAgIrbtsm0CMLubtrNr73Vt90hmPl5nP5IkSZLUFE0NWJn5DPAT4NSIGBwRbwf2Ai7upvkPgI9HxOsjYmPg88D3e9GPesZTKJvPMW8+x7y5HO/mc8ybzzFvLse7+RzzOjV1kQuo3r8K+B6wO/A4cHxm/igixgB/Bl6fmQ/W2h4DHAcMAq4CDsvMF1bVT1O/jCRJkiR10fSAJUmSJEl9VStuNCxJkiRJfZIBS5IkSZIKMWCpWxERra5BaiT3cUmS1AgGLHUraxfneRDaHI5z87mPN59j3RqOu/o693GtbVzkQq8QEZ8FhgMPAL/LzDsiYp3MXNLi0vqkiHhTZv6h9jzS/yAbzn28+dzPmy8ihgELu/whwX28wSJiQGa+1Oo6+gv38dZwP+8ZZ7C0TET8HPgQMAx4F/DLiHhvZi6JCPeVwiLiR8CsiPgoVGdU/CtcY7mPN5/7efNFxH8APwOuiogvALiPN05EfDoiRmTmSxGxbqvr6Q/cx5vP/bw+7ogCICImAFtk5i6Z+Ungk1RvLHdNRLy/9ovLg6JCIuIIYAfgZODzEfEx8OCzkdzHm8/9vPki4nzgDcDngFuAf4yI34EHoI0QERcCZ1H9Y82mmfmyB5+N5T7efO7n9XMn1FJ/AwZHxG4AmbkgM08HjgcujIi3eFpPGRHxKuBB4AvAeVR/aZ3swWfDuY83kft580XEIGAr4POZ+f+ArwIfrr4Vt8KyA1DHvYCI2B0YR3VW/D7gag8+G8t9vPncz3vHgKWlXgLuB3ZZ7j+Yc4FLgX0iYl1/aa25zFwM/Ba4NjOfAC4BTqd68PnxWpuMiE1bWGZf9DzVffzN7uONV9vPO4Fr3M+bIzOfAxYA/xAR62Xmy5k5D9gbeD4ivlNr5x8SypgDfBH4NXAS1es6uz349PdKGe7jLfEX3M/rZsDqxyJiYkS8KSIGZeZ8qqdLTQM+urRNZj4LzAVeV/tF5i+tXlpuvP+Wmc8AZOYi4IfAGcBJEbF37QD0hxEx0F9YvVebRQEgMx8Fvkt1H/9Il+3u4wUtN+YLa+Prft48fwL+Cdh66Ybavn8WMC4iRraqsL4mM+8F/rO2n88BTqR68HlNl4PPjlpbf6+UczvwPtzHmyIz7wdudj+vz4BWF6DWiIifAtsALwCbRsThmfnziPgE8P2IGAxcn5n31D6SETEwM59vVc3tbLnx3iQiDgf+OzOfAsjMZyLih8DjwE+ozra8w/Huvdp5+rMj4ru1v3qSmT+J6spT36rt479yHy+nuzHvyv28vNp1bi8CSzLzPzLzSxHxNqoX/38AmFM76JkJjACGAo+0ruL21mW8X87M73bddzPz/oj4PHAa8OOI+C3wiaiuojm/RSW3vYjYAnho6cF7Zp5ZO9X7xxGxJ+7jxXUz5i8sfc/9vGcMWP1QRBxJ9XzaXYCk+hf9L0bE2Mw8LyIWU70u5aCIeBx4K/BOD4J6ZyXjPR34TkRcUTt9isxcFBFvBJ4G/iEz/9yiktteRHyd6nn5TwDPRcSPuoSs70fEc8CngQ+7j5exqjHvyv28nIj4JbAl8N/AP9cC1fGZuXdEXA1cBZwZEbcBk4ANgKdaVnCb62a830t11uS2Lgei99WuM7wdOBJ4lwedvRcRV1KdqfoI1dlZADLz/RFxA/Bj4Gz38XJWNuZduZ+vngGrfxoG3NLl4OeEiHgY+GhELMjMKyLiTuDVwGbA4Zn5QKuK7QNWOt5UzyW/snYO8xjgg8B7POjsvdpf3rakugz766gusrBORPywS8i6PCL+AGyB+/ga68mYd2k7Dvhn3M/XSERMBDbPzB1rrz8LXAacEREnZ+YHIuLLVH/PnED1//cHZuZjraq5na1ivI8FzqF6Xe1Sn6J6gLpjZt7R7Fr7iog4ler/F58CvhYRR2fmnyKq99LLzMnu42WtbMxX0tz9fBW80XA/FBH7Ur0OYs/MnN1l+8nAYVT/Y1nQqvr6mtWM96HAhMx8PCLWAwZl5tMtKrXPiIitgacyc0HtOp8KcCrwo6XXvqmsno557Rqtge7na6Z2zcMlwGuXjm9EbAxcCJCZe9e2bQasDzybmY+3ota+YDXj/SLw4cx8LiKGUA1dV2fmbS0qt0+onY75dGZeXJs9HAb8Wy1kLbupcESMAAbhPr7GVjXmy7XbgOqZTr90P++eAaufiIh/BTYFZgHPUP0L8rrA1zPzvi7t/hu4MDMvaEmhfYTj3XzLjfldtYtxl7637IA/M78TEQcD/y8z72xFrX1FnWP+UWCmY15GVO/1MwO4meoYZ237hlRP67kwMyutq7Bv6cF4/0dmfqG2bd3MfLllxfYhETF86Wn0XQ74PwXc7oIKjbG6MV86g+h+vmoGrH4gIn4OjAb+TPX+EX8BlgDPUb0O4ntLD3oi4ifAVZl5SYvKbXuOd/N1GfPZwLbAPODizPxFlzafAI4D7gQ+ALyh64yi6uOYN19EfBj4r8x8oDbj/a9Urzv5BXBRl4P+fwfGZ+YnWldt+3O8m6/rmHfZNiAzX6o9v5rqIhb7Uf33MSgzj29JsX1EL8Z8cGZ+tiXFthGvwerjIuLtVFev2zGrS2lOBPagelH/U8Bi4Nu1/4AGAu+kOu2rXnC8m6+bMX8T1Rsifrb2F7afAmTmf0R1xam3Azt5oN97jnnzRXX1xT2AyyPiS1ldyetSYCzV8LoJ8OVa882B9bueRqX6ON7Nt9yYn5GZf6nNlry0dLakdm3hZVRnDTcG/rGlRbc5x7xxvA9W3/c01QP7EbVf/rOAHwD/CQTwP1Qv1H0X1b9CvzP/vmy16ud4N9/yY/4Hqvd0mwF8PCJ2BoiIQ6keGE3KzNtbVm3f4Jg3UUR8iOqiQ18A1gOOi4hta9ebnEH1FM33RMS9EfFdqtfSnuXBfu843s3XzZgfHxFbdTkl7eWIWDopMKvW5o2ZObNFJbc9x7yxPEWwj4uIUcBtVH/5f7XL9h2Bz1G9D9CFtdMfXvZ/EGvG8W6+VYz5eKrXAF1bG/M3A89n5v+2ptK+wzFvrojYnOos+C+A9wL7Uz3l+MzMvLd2wfkQqqfwLKK6aundraq33TnezbeSMX8W6DqrklG9x9t/Am+r/QFTveSYN5YBqx+IiL2Ai4EjaivDLP2P5nRgAvABD/TLcbybbzVjviOw19LzyVWGY95cEbFeZr5Ye74nMIXqQf+XsnpPmtcD92Xm4lbW2Vc43s23ijE/vXbAv0Nm3hERIzPTGwkX4Jg3jtdg9Q+/oHqh+VcjYnBmfqu2fSHwf1RPXVM5jnfzrWrMH6Z6g2eV5Zg3UWa+uDTEZuYvIiKAfYHDo3ofvd2A3aneW09ryPFuvlWM+RG1MZ8UER0e6JfjmDeOM1j9RO0/lP2A7wC3UP0LxTuoXgP0x1bW1hc53s3nmDefY958Sw+Gas/fBnwfeA3V69w8facwx7v5HPPmc8zLM2D1MxHxWqrn3L4KuDkz721xSX2a4918jnnzOebN1eV0zKOprmQ3ITPvaHFZfZbj3XyOefM55mV5imA/k9UV61y1rkkc7+ZzzJvPMW+u2kHQEOANwC4eBDWW4918jnnzOeZlOYMlSVIb6nqBuhrP8W4+x7z5HPMyDFiSJEmSVIg3GpYkSZKkQgxYkiRJklSIAUuSJEmSCjFgSZIkSVIhBixJkiRJKsSAJUmSJEmFGLAkSZIkqZD/D1/qqWFhU+uPAAAAAElFTkSuQmCC\n",
+      "text/plain": [
+       "<Figure size 864x576 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "s = 3\n"
+     ]
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAI4CAYAAAB3HEhGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA/yElEQVR4nO3deZxcVZ3//9cHGknIQgiBKCFLsxqDRCQiDoNpogQdZXGQsERUVBgQRhlcABFSIBEEVEABxVFERHbEhUWMSTNfZX6BqBOHyC5ZJCMQEEjYIZ/fH1WJnaaT7urcqkqlX8/Hox5U3Tp1+lPHm/a++55zb2QmkiRJkqS1t0GjC5AkSZKk9YUBS5IkSZIKYsCSJEmSpIIYsCRJkiSpIAYsSZIkSSpIS6MLqLdhw4blmDFjGl2GJEmSpCb2+9//fklmbtF5e58LWGPGjGHOnDmNLkOSJElSE4uIBV1td4qgJEmSJBXEgCVJkiRJBTFgSZIkSVJBDFiSJEmSVJA+d5ELSZIkqaeeffZZHn/8cV555ZVGl6I622ijjdhyyy0ZPHhwVZ8zYEmSJEldePbZZ3nssccYMWIE/fv3JyIaXZLqJDN54YUXePTRRwGqCllOEZQkSZK68PjjjzNixAg22WQTw1UfExFssskmjBgxgscff7yqzxqwJEmSpC688sor9O/fv9FlqIH69+9f9fRQA5YkSZK0Gp656tt687+/AUuSJEmSCmLAkiRJkprQuHHjaG9vf9329vZ2tt5660J+RkTw0EMPFdJXX+FVBCVJkqQe2m36jJr2f9cp7+1x23nz5tWwEvWWZ7AkSZIkqSAGLEmSJKkJjRkzhhkzZvDCCy/w8Y9/nM0224y3vOUt3H333au0W7x4MQceeCBbbLEFra2tXHjhhSvfu+uuu3jXu97FkCFDeNOb3sRxxx3Hyy+/XO+vsl4xYEmSJElN7PTTT+fhhx/m4Ycf5le/+hWXX375yveWL1/Ovvvuy/jx43n00Uf5zW9+w/nnn8+vfvUrADbccEO++c1vsmTJEv77v/+b3/zmN1x88cWN+irrBQOWJEmS1MSuvfZaTjnlFIYOHcrIkSP5zGc+s/K9u+++myeeeILTTjuNN7zhDWyzzTYceeSRXH311QDsuuuu7L777rS0tDBmzBj+7d/+jTvuuKNRX2W94EUuJEmSpCa2ePFiRo4cufL16NGjVz5fsGABixcvZsiQISu3vfbaa+y5554APPDAA5xwwgnMmTOH559/nldffZVdd921brWvjzyDJUmSJDWxN73pTSxatGjl64ULF658PnLkSFpbW3n66adXPpYuXcott9wCwDHHHMOb3/xmHnzwQZ599lm++tWvkpl1/w7rEwOWJEmS1MSmTJnCWWedxd///nf++te/8q1vfWvle7vtthuDBw/ma1/7Gi+88AKvvfYa99xzz8oLYSxdupTBgwczcOBA7rvvPi655JJGfY31hlMEJUmSpB6q5j5V9TJt2jSOPvpoWltb2WqrrTjiiCO44IILgPJFLH7xi1/wuc99jtbWVl566SV23HFHzjzzTADOO+88jjrqKM455xx22WUXDj74YGbOnNnIr9P0ot6nACNiKPB9YDKwBDg5M3/SRbtDgNOBNwIvAbcC/56Zz1bTT2cTJkzIOXPmFPRtJEmStL669957GTt2bKPLUIOtbj+IiN9n5oTO2xsxRfAi4GVgODAVuCQixnXR7nfAHpm5KbAN5bNtZ/aiH0mSJEmqi7oGrIgYABwInJqZyzLzt8DPgcM7t83MRZm5pMOm14Dtqu1HkiRJkuql3muwdgBey8wHOmybC0zsqnFE/DNwMzAYeB74UC/7OQo4CmD48OG0t7evxVeQJEmNdN/flgLwjTNO5sF75xXW7/Zjx3Hpxd8urD81v0033ZSlS5c2ugw12IsvvlhVfqh3wBoIPNNp2zPAoK4aV85MbRoRI4Ajgfm97OdS4FIor8Fqa2vrRemSJGld8MXpMwAYMvV83lFw3x4jqKN7772XQYO6PLxUH9KvXz922WWXHrev9xqsZZTPRnU0GFjjnwYy81HgNuDqtelHkiRJkmqp3gHrAaAlIrbvsG080JPz+y3AtgX0I0mSJEk1UdeAlZnPATcCZ0TEgIjYA9gfuKJz24iYGhGjomw0MB34TbX9SJIkSVK9NOIy7Z8G+gOPA1cBx2TmvEqYWhYRoyrt3gLcSXk64O+A+ymvw1pjP3X6DpIkSZL0OvW+yAWZ+RRwQBfbF1K+eMWK16cAp1TbjyRJkiQ1St0DliRJktS0ztuxtv1//v4eN73//vs55JBDeOihh5g+fTqf+cxnumw3f/58WltbeeWVV2hp8fC/o1qMjSMsSZIkNaFzzjmHtrY2/vjHPza6lJoplUo89NBD/PjHP250KT3WiDVYkiRJktbSggULGDduXKPLaKjMZPny5Y0uYxUGLEmSJKnJTJo0iVmzZnHccccxcOBALrjgAnbZZRcGDx7MyJEjKZVKq/3sDTfcwJgxY7jnnntYvnw5Z599Nttuuy2bb745U6ZM4amnnlrjz54/fz4RwWWXXcbIkSPZbLPN+M53vsPdd9/NzjvvzJAhQzjuuONWtl++fDlnnnkmo0ePZsstt+SjH/0ozzzzzCp9XX755YwaNYphw4Yxffp0AG677Ta++tWvcs011zBw4EDGjx8PlG8Ifsopp7DHHnuwySab8Je//IX77ruPvffem6FDh7Ljjjty7bXXrvz5N998c4/HpggGLEmSJKnJzJw5kz333JNvf/vbLFu2jPHjx/OjH/2Ip59+mptvvplLLrmEm2666XWfu+yyyzjxxBOZMWMGO+20ExdeeCE33XQTd9xxB4sXL2azzTbj2GOP7VENs2fP5sEHH+Saa67h+OOPZ/r06cyYMYN58+Zx7bXXcscddwDwwx/+kB/+8IfMmjWLv/zlLyxbtmyVAAbw29/+lvvvv5/f/OY3nHHGGdx77728733v40tf+hIHH3wwy5YtY+7cuSvbX3HFFVx66aUsXbqULbbYgr333pvDDjuMxx9/nKuuuopPf/rTzJtXvsD4gAEDejQ2RTFgSZIkSU2ura2Nt771rWywwQbsvPPOHHrooSsDzgrnn38+5557Lu3t7Wy33XYAfPe732X69OlsvfXWbLzxxpRKJa6//npeffXVbn/mqaeeSr9+/Zg8eTIDBgzg0EMPZcstt2TEiBHsueeeK9eGXXnllZxwwglss802DBw4kLPOOourr756lZ8xbdo0+vfvz/jx4xk/fvwqYaorH//4xxk3bhwtLS3cdtttjBkzhiOOOIKWlhbe/va3c+CBB3L99df3eGyK5EUuJEmSpCY3e/ZsTjrpJO655x5efvllXnrpJQ466KBV2px77rmcdtppbL311iu3LViwgA996ENssME/zrtsuOGGPPbYY4wYMWKNP3P48OErn/fv3/91r5ctWwbA4sWLGT169Mr3Ro8ezauvvspjjz22ctsb3/jGlc832WSTlZ9dnZEjR67yHWbPns2QIUNWbnv11Vc5/PDDgZ6NTZE8gyVJkiQ1ucMOO4z99tuPRYsW8cwzz3D00UeTmau0uf322znzzDO54YYbVm4bOXIkt956K08//fTKx4svvthtuKrGVlttxYIFC1a+XrhwIS0tLasEstWJiG63jxw5kokTJ67yHZYtW8Yll1wC9GxsiuQZLEmSJKmnqrhPVT0tXbqUoUOH0q9fP+666y5+8pOfMHny5FXajBs3jttuu4199tmHjTbaiP3224+jjz6aU045hcsvv5zRo0fzxBNPcOedd7L//vsXVtuhhx7K1772Nd7//vezxRZbrFxX1ZP7Tg0fPpxf//rXLF++fJWzbB198IMf5KSTTuKKK67gkEMOAeB//ud/GDhwIGPHju3R2BTJM1iSJElSk7v44os57bTTGDRoEGeccQZTpkzpst348eP55S9/yZFHHsmtt97KZz/7Wfbbbz8mT57MoEGD2H333Zk9e3ahtX3iE5/g8MMP593vfjetra3069ePb33rWz367IqpfJtvvjlvf/vbu2wzaNAgbr/9dq6++mq22mor3vjGN3LiiSfy0ksvAT0fm6JELU+PrYsmTJiQc+bMaXQZkiSpl3abPqNmfd91yntr1reaz7333svYsWMbXYYabHX7QUT8PjMndN7uGSxJkiRJKogBS5IkSdIqrrzySgYOHPi6x7hx4xpd2jrPi1xIkiRJWsXUqVOZOnVqo8toSp7BkiRJkqSCGLAkSZIkqSAGLEmSJEkqiAFLkiRJkgpiwJIkSZKkgngVQUmSJKmHJl07qab9z5wys8dt77//fg455BAeeughpk+fzmc+85ku282fP5/W1lZeeeUVWlqa//A/InjwwQfZbrvtqv7suHHjuOiii2hrayu+sIrmH2FJkiSpDzrnnHNoa2vjj3/8Y6NLaRrz5s2r+c9wiqAkSZLUhBYsWOCNf9dBBixJkiSpyUyaNIlZs2Zx3HHHMXDgQC644AJ22WUXBg8ezMiRIymVSqv97A033MCYMWO45557WL58OWeffTbbbrstm2++OVOmTOGpp55a48+eP38+EcFll13GyJEj2WyzzfjOd77D3Xffzc4778yQIUM47rjjVrYvlUp85CMfed3nX331VQDa2tr48pe/zD/90z8xcOBA9t13X5588kmmTp3K4MGDecc73sH8+fNXqeGWW25hm222YdiwYXzhC19g+fLlADz88MNMmjSJzTffnGHDhjF16lSefvrplZ8bM2YMM2bM6OEo944BS5IkSWoyM2fOZM899+Tb3/42y5YtY/z48fzoRz/i6aef5uabb+aSSy7hpptuet3nLrvsMk488URmzJjBTjvtxIUXXshNN93EHXfcweLFi9lss8049thje1TD7NmzefDBB7nmmms4/vjjmT59OjNmzGDevHlce+213HHHHT3+PldffTVXXHEFjz76KA8//DDvete7OOKII3jqqacYO3Ysp59++irtf/rTnzJnzhz+8Ic/8LOf/Ywf/OAHAGQmJ598MosXL+bee+9l0aJFawybtWDAkiRJkppcW1sbb33rW9lggw3YeeedOfTQQ18XcM4//3zOPfdc2tvbV14g4rvf/S7Tp09n6623ZuONN6ZUKnH99devPLu0Jqeeeir9+vVj8uTJDBgwgEMPPZQtt9ySESNGsOeee1a1NuyII45g2223ZdNNN+X9738/2267Le9973tpaWnhoIMOel1fJ554IkOHDmXUqFEcf/zxXHXVVQBst9127L333my88cZsscUWnHDCCVUFvSIYsCRJkjoplUpERGGPev8FXX3P7Nmz2Wuvvdhiiy3YdNNN+c53vsOSJUtWaXPuuedy7LHHsvXWW6/ctmDBAj70oQ8xZMgQhgwZwtixY9lwww157LHHuv2Zw4cPX/m8f//+r3u9bNmyHtdfbV8jR45c+Xz06NEsXrwYgMcff5xDDjmEESNGMHjwYD7ykY+8bhxqzYAlSZLUSalUIjPX+Jg4cSITJ07stl1mGrBUc4cddhj77bcfixYt4plnnuHoo48mM1dpc/vtt3PmmWdyww03rNw2cuRIbr31Vp5++umVjxdffJERI0YUVtuAAQN4/vnnV77+29/+ttZ9Llq0aOXzhQsXstVWWwFw8sknExH86U9/4tlnn+XHP/7x68ah1rxMuyRJktRD1dynqp6WLl3K0KFD6devH3fddRc/+clPmDx58iptxo0bx2233cY+++zDRhttxH777cfRRx/NKaecwuWXX87o0aN54oknuPPOO9l///0Lq+1tb3sbX/va11i4cCGbbropZ5111lr3ee655/LOd76TZcuWccEFF3DCCScA5XHYdNNNGTJkCI8++ijnnnvuWv+sankGS5IkSWpyF198MaeddhqDBg3ijDPOYMqUKV22Gz9+PL/85S858sgjufXWW/nsZz/Lfvvtx+TJkxk0aBC77747s2fPLrS2vffem4MPPpidd96ZXXfdlQ9+8INr3ef+++/Prrvuytve9jY+8IEP8MlPfhKAadOm8Yc//IFNN92UD3zgA/zrv/7rWv+sakW9T5k12oQJE3LOnDmNLkOSJPXSbtNrd4nlu055b4/btrW1AdDe3l6bYtRw9957L2PHjm10GWqw1e0HEfH7zJzQebtnsCRJkiSpIAYsSZIkSau48sorGThw4Ose48aNa3Rp6zwvciFJkiRpFVOnTmXq1KmNLqMpeQZLkiRJWo3ly5c3ugQ1UG/+9zdgSZIkSV0YMGAAjz76KC+//HLd76WkxspMXn75ZR599FEGDBhQ1WedIihJkiR1Yeutt2bJkiUsWLCAV199tdHlqM5aWlrYdNNNGTZsWHWfq1E9kiRJUlPbYIMN2HLLLdlyyy0bXYqaiFMEJUmSJKkgBixJkiRJKohTBCVJklY4b8eet120qLrPfP7+6uuR1HQ8gyVJkiRJBTFgSZIkSVJBDFiSJEmSVBADliRJkiQVxIAlSZIkSQUxYEmSJElSQQxYkiRJklQQA5YkSZIkFcSAJUmSJEkFMWBJkiRJUkEMWJIkSZJUEAOWJEmSJBWkpdEFSJIkrWtKty/h9F8/1aO28YUHum0zbe+hlD6/tlVJagYGLEmSpE5Kk4dRmjys0WVIakJOEZQkSZKkghiwJEmSJKkgBixJkiRJKogBS5IkSZIKYsCSJEmSpIIYsCRJkiSpIAYsSZIkSSqIAUuSJEmSCmLAkiRJkqSCGLAkSZIkqSAGLEmSJEkqiAFLkiRJkgpiwJIkSZKkghiwJEmSJKkgBixJkiRJKogBS5IkSZIKYsCSJEmSpIIYsCRJkiSpIAYsSZIkSSpI3QNWRAyNiJ9GxHMRsSAiDltNu49FxO8j4tmI+GtEnBMRLR3eb4+IFyNiWeVxf/2+hSRJkiS9XiPOYF0EvAwMB6YCl0TEuC7abQIcDwwD3gm8B/h8pzbHZebAymPH2pUsSZIkSd1r6b5JcSJiAHAgsFNmLgN+GxE/Bw4HTurYNjMv6fDy0Yi4EtirbsVKkiRJUpXqGrCAHYDXMvOBDtvmAhN78Nl3A/M6bTsrIs4G7gdOycz2rj4YEUcBRwEMHz6c9vYum0mSpCbwiW1fqFnf7fGZmvWNxx9Sn1DvgDUQeKbTtmeAQWv6UEQcAUwAPtVh84nAnylPNzwE+EVEvC0zH+78+cy8FLgUYMKECdnW1tbb+iVJUoN9cfqMmvV918YX1qxvDna5uNQX1HsN1jJgcKdtg4Glq/tARBwAnA28PzOXrNiembMzc2lmvpSZlwO/A/6l+JIlSZIkqWfqHbAeAFoiYvsO28bz+ql/AETE+4DvAftm5v9203cCUUiVkiRJktQLdQ1YmfkccCNwRkQMiIg9gP2BKzq3jYhJwJXAgZl5V6f3hkTEPhHRLyJaImIq5TVav6r9t5AkSeq9UqlERBT2KJVKjf5KkjpoxGXaPw30Bx4HrgKOycx5ETGqcj+rUZV2pwKbArd0uNfVrZX3NgLOBJ4AlgD/DhyQmU5uliRJ67RSqURmrvExceJEJk6c2G27zDRgSeuYel/kgsx8Cjigi+0LKV8EY8Xr1V6SPTOfAN5Ri/okSZIkqbcacQZLkiRJktZLBixJkiRJKogBS5IkSZIKYsCSJEmSpIIYsCRJkiSpIAYsSZIkSSqIAUuSJEmSCmLAkiRJkqSCGLAkSZIkqSAGLEmSJEkqiAFLkiRJkgpiwJIkSZKkghiwJEmSJKkgLY0uQJIkqS+YdO2kHred+8Tcqj4zc8rMXtUkqXiewZIkSZKkghiwJEmSJKkgBixJkiRJKogBS5IkSZIKYsCSJEmSpIIYsCRJkiSpIAYsSZIkSSqIAUuSJEmSCmLAkiRJkqSCGLAkSZIkqSAGLEmSJEkqiAFLkiRJkgpiwJIkSZKkghiwJEmSJKkgBixJkiRJKkhLowuQJEnqSx657hHmXz+/R21nHTyr2zZjPjwGpqxdTZKKY8CSJEmqo9aDWmk9qLXRZUiqEacISpIkSVJBDFiSJEmSVBADliRJkiQVxIAlSZIkSQUxYEmSJElSQQxYkiRJklQQA5YkSZIkFcSAJUmSJEkFMWBJkiRJUkEMWJIkSZJUEAOWJEmSJBXEgCVJkiRJBTFgSZIkSVJBDFiSJEmSVBADliRJkiQVxIAlSZIkSQUxYEmSJElSQQxYkiRJklQQA5YkSZIkFcSAJUmSJEkFMWBJkiRJUkEMWJIkSZJUEAOWJEmSJBXEgCVJkiRJBTFgSZIkSVJBDFiSJEmSVBADliRJkiQVxIAlSZIkSQUxYEmSJElSQQxYkiRJklQQA5YkSZIkFcSAJUmSJEkFMWBJkiRJUkEMWJIkSZJUEAOWJEmSJBXEgCVJkiRJBTFgSZIkSVJBDFiSJEmSVBADliRJkiQVxIAlSZIkSQUxYEmSJElSQQxYkiRJklSQugesiBgaET+NiOciYkFEHLaadh+LiN9HxLMR8deIOCciWqrtR5IkSZLqpRFnsC4CXgaGA1OBSyJiXBftNgGOB4YB7wTeA3y+F/1IkiRJUl3UNWBFxADgQODUzFyWmb8Ffg4c3rltZl6Smf8vM1/OzEeBK4E9qu1HkiRJkuqlpfsmhdoBeC0zH+iwbS4wsQeffTcwrzf9RMRRwFEAw4cPp729vcqyJUnSuuIT275Qs77b4zM16/uQjTaqWd+1Pra5729L+cYZJ/PgvfO6b9xD248dxwmnncWb3ziosD6ldUG9A9ZA4JlO254B1vgvKyKOACYAn+pNP5l5KXApwIQJE7Ktra2qoiVJ0rrji9Nn1Kzvuza+sGZ9nzFqRM36ntk2s2Z9Q3nMh0w9n3cU3O8PHoa7DmkruFepseodsJYBgzttGwwsXd0HIuIA4GzgvZm5pLf9SJIkSVKt1fsiFw8ALRGxfYdt4/nH1L9VRMT7gO8B+2bm//a2H0mSJEmqh7oGrMx8DrgROCMiBkTEHsD+wBWd20bEJMoXtjgwM+/qbT+SJEmSVC+NuEz7p4H+wOPAVcAxmTkvIkZFxLKIGFVpdyqwKXBLZfuyiLi1u37q9zUkSZIkaVX1XoNFZj4FHNDF9oWUL16x4vVevelHkiRJkhqlEWewJEmSJGm9ZMCSJEmSpIIYsCRJkiSpIAYsSZIkSSqIAUuSJEmSCmLAkiRJkqSCGLAkSZIkqSAGLEmSJEkqiAFLkiRJkgpiwJIkSZKkghiwJEmS1FClUomIKOxRKpUa/ZXUh7U0ugBJkiT1baVSqdtQ1NbWBkB7e3vN65HWhmewJEmSJKkgBixJkiRJKohTBCVJktQ45+3Ys3aLFlXXHuDz91dfj7SWPIMlSZIkSQUxYEmSJElSQQxYkiRJklQQA5YkSZIkFcSLXEiSJKmhSrcv4fRfP9WjtvGFB7ptM23voZQmD1vbsqReMWBJkiSpoUqThxmItN5wiqAkSZIkFcSAJUmSJEkFMWBJkiRJUkEMWJIkSZJUEAOWJEmSJBXEgCVJkiRJBakqYEXETrUqRJIkSZKaXbVnsP4UEXdHxDERMaQWBUmSJElFKZVKRERhj1Kp1OivpHVctQHrPcCfgXOAxRFxVUTsHRFRfGmSJEnS2imVSmTmGh8TJ05k4sSJ3bbLTAOWulVVwMrMWZn5MeCNwHHACOBXwIKI+EpEbFuDGiVJkiSpKfTqIheZ+Vxm/iAz3w3sAMwHvgQ8EBF3RMSHCqxRkiRJkppCr68iGBFjIqIE3A68C7gFOAp4DLgmIr5ZSIWSJEmS1CSqvYrgJhHx0YiYBTwETAW+B4zKzH0z8/uZOQX4N+CTxZcrSZIkSeuulirb/w3YELgReG9mtq+m3d3Ak2tRlyRJkiQ1nWoD1knAlZn5zJoaZeY9QGuvq5IkSZKkJlTtGqxhwICu3oiIN0XEaWtfkiRJkiQ1p2oD1jRg69W8t1XlfUmSJEnqk6oNWAHkat7bGvj72pUjSZIkSc2r2zVYEfEx4GOVlwlcEhHPdmrWD3gr5Uu2S5IkSQ036dpJPWo394m5VbUHmDllZq9q0vqvJxe5eJ5/XBEwgGeApzq1eRm4Fbi4uNIkSZIkqbl0G7Ay8zrgOoCIuAz4Smb+pdaFSZIkSVKzqeoy7Zl5RK0KkSRJkqRm15M1WOcAF2bmXyvP1yQz88RiSpMkSZKk5tKTM1gHAVcCf608X5MEDFiSJEmS+qSerMFq7eq5JEmSJGlV1d4HS5IkSZK0Gj1Zg/Uv1XSYmbf0vhxJkiRJal49WYP1S8prq6IHbRPYcK0qkiRJkqQm1ZOA5borSZIkSeqBnlzkYkE9CpEkSZKkZteTNVibZObzK553135FW0mSJKnRHrnuEeZfP79HbWcdPKvbNmM+PIbWg5zgpdXryRTBpRHxrsy8C1hGeZ3VmrgGS5IkSeuE1oNaDUSqq54ErE8AD3d43l3AkiRJkqQ+qSdrsC7v8PyHNa1GkiRJkppYT85gvU5EDAF2At4E/B9wT2Y+XVxZkiRJktR8qgpYEdECTAeOBTpe8OL5iLgYOCUzXymwPkmSJElqGtWewfoGcBRwBnAj8DiwJXAgcCrQD/hMkQVKkiRJUrOoNmAdDnwpM7/RYdtTwPSIeBH4MgYsSZIkSX3UBlW2Xw7MW8179+AVBiVJkiT1YdUGrCuAT63mvSOBH69dOZIkSZLUvLqdIhgRn+7wcj7w4YiYB/ycf6zB2h8YBJxXgxolSZIkqSn0ZA3Wt7vYthUwtovt3wAuWKuKJEmSJKlJ9eRGw9VOI5QkSZKkPsnwJEmSJEkFqfYy7QBExNbADpTve7WKzLxlbYuSJEmSpGZUVcCKiEHAtcDkFZsq/+14efYNC6hLkiRJkppOtVMEzwJGAXtSDlcfAtqA7wOPALsXWZwkSZIkNZNqA9a/ANOB2ZXXizPzvzLzKOBnwBeKLE6SJEmSmkm1AWs4sCgzXwOeA4Z2eO8W/jF1UJIkSZL6nGoD1iJgWOX5g8AHO7z3TuDFIoqSJEmSpGZU7VUEfw28F/gp8E3g8ojYFXgJeDfw9WLLkyRJkqTmUW3AOhHYBCAzr4iIZcCHgf7AccB3iy1PkiRJkppHVVMEM/P5zFzS4fVPM3NqZv5rZl6Smcu76yMihkbETyPiuYhYEBGHrabdThHxq4hYEhHZxfvtEfFiRCyrPO6v5rtIkiRJUtF6e6PhHYF3AG8CFgNzMrOnAeci4GXKF8x4G3BzRMzNzHmd2r1C+Z5bFwM3raav4zLzP6urXpIkSZJqo9obDQ8GvgccSPns1zJgILA8Im4EPpWZz67h8wMqn90pM5cBv42InwOHAyd1bFsJbPdHxHbV1ChJkiRJjRKZr5t9t/rGET8GPgAcC9yYmS9GRD/KoenbwM2Z+ZE1fH4X4M7M7N9h2+eBiZm572o+sx3wYGZGp+3twDjKNzy+HzglM9tX08dRwFEAw4cP3/Xqq6/u2ReWtNbu+9tSAL5xxsk8eG/nE9W9t/3YcVx68bcL609S81jxe6UW3hwLa9b3A2/YqGZ977DZDjXrGxzzrtR6zLXu22uvvX6fmRM6b682YC0F/qOraXkRcSTwjcwctIbP7wlcl5lv7PS5qZnZtprPrC5gvRP4M+XphodQDnhvy8yH1/QdJkyYkHPmzFlTE0kF2m36jJr1fdcp7wWgVCpx+umnF9bvtGnTKJVKhfUnqVg1/b2y8bE163vSqBE163vmlJk16xsc867Uesy17ouILgNWtffBWgb832reW0z55sPdfX5wp22Dgar/LJKZszNzaWa+lJmXA78D/qXafiQ1v1KpRGau8TFx4kQmTpzYbbvMNFxJkqReqzZgXQR8PiL6d9wYEZsAn6d8QYo1eQBoiYjtO2wbDxQxbygpTxeUJEmSpIboNmBFxDkrHpTPNm0PLIqIqyLigoi4ClgIbAesdnogQGY+B9wInBERAyJiD2B/4Ioufm5U1ne9ofK6X0RsXHk+JCL2qWxriYiplG90/KsqvrskdatUKhERhT08OyZJ0vqtJ1cRPKjT61cqj907bFsxxe9A4Avd9Pdp4AfA48CTwDGZOS8iRlFeU/WWzFwIjAYe6fC5F4AFwBhgI+BM4M3Aa8B9wAFVXCpeknqkVCp1G4ra2toAaG9vr3k9kiRp3dZtwMrM1iJ/YGY+BRzQxfaFlC/5vuL1fFYz5S8zn6B8Hy5JkiRJWmf06kbDkrROOG/HnrddtKi6z3zeE+KSJKl61V7kgojYJiIuiYj/jYhHK/+9OCK2qUWBkiRJktQsqjqDFRG7ArOAF4FfAo8BwymvvZoaEXtl5h8Kr1KSJEmSmkC1UwTPA/4IvD8zn1+xsXKZ9lsq708qrjxJkiRJah7VBqzdgCkdwxVAZj4fEecB1xRWmSQ10KRre/63orlPzK3qMzOnzOxVTZIkad1XbcB6Adh8Ne8NpTx1UJLqqnT7Ek7/9VM9ahtfeKDbNtP2HgqfGrG2ZUmSpD6o2oB1M3B2RPwlM3+7YmNE/DNwFvCLIouTpJ4oTR5GafKwQvt0rrMkSeqNagPWCcDPgDsi4gnKF7nYsvK4E/hcseVJkiRJUvOoKmBl5pPAP0fE+yjf6PdNwP8BszPz9hrUJ0mSJElNo8cBKyL6AX8CPpOZtwG31awqSZIkSWpCPQ5YmfliRAwBlteuHElatzxy3SPMv35+j9rOOnhWt23GfHgMTFm7miRJ0rqr2jVYVwJHAE4HlNQntB7USutBrY0uQ5IkNYlqA9ZCYEpEzKF8Y+HHgOzwfmbmJUUVJ0mSJEnNpNqA9fXKf98EvL2L9xMwYEmSJEnqk6q9iuAGtSpEkiRJkppdtWewiIg3AB8HdqPDZdqByzPz5UKrkyRJkqQmUtUZqYgYCzwIXATsBLxW+e9FwEMR8ZbCK5QkSZKkJlHtGaxLgWeAPTNz4YqNETEKuBn4DvDu4sqTJEmSpOZR7ZqqCcBpHcMVQOX1acA7iipMkiRJkppNtQFrPtBvNe/1o3wZd0mSJEnqk6oNWCcBZ0bEOztujIjdgTOAE4sqTJIkSZKaTbVrsL4MDAbujIjHgceBLSuPJ4EvRcSXVjTOzN2KKlSSJEmS1nXVBqx7Kg9JkiRJUifV3mj4iFoVIkmSJEnNrto1WJIkSZKk1TBgSZIkSVJBDFiSJEmSVBADliRJkiQVxIAlSZIkSQUxYEmSJElSQQxYkiRJklQQA5YkSZIkFcSAJUmSJEkFMWBJkiRJUkEMWJIkSZJUEAOWJEmSJBXEgCVJkiRJBTFgSZIkSVJBDFiSJEmSVBADliRJkiQVxIAlSZIkSQUxYEmSJElSQQxYkiRJklQQA5YkSZIkFcSAJUmSJEkFMWBJkiRJUkEMWJIkSZJUEAOWJEmSJBXEgCVJkiRJBTFgSZIkSVJBDFiSJEmSVBADliRJkqTClEolIqKwR6lUavRXqkpLowuQJEmStP4olUrdhqK2tjYA2tvba15PvXkGS5IkSZIKYsCSJEmSpIIYsCRJ6sP6+loJqa/y337tuAZLkqR1XKlU4vTTTy+sv2nTpq08GOrrayWkvqpW//YnXTupR+3mPjG3qvYAM6fM7HHbRjJgSZK0jjMESVLzMGBJkiRJfcV5O/a87aJF1X9m1Ijq6lkPuQZLkiRJWofsNn0Gu02fwaDW8YWukxrUOr7RX61P8AyWJEmStA4a+6mv16DXY2vQ56oeue4R5l8/v0dtZx08q9s2Yz48htaDWteyqvoxYEmStK6q9VSez99fXT2S1AOtB7U2VSAqmgFLkiRJ6mNKty/h9F8/1aO28YUHum0zbe+hlCYPW9uy1gsGLEmS+igvpyz1XaXJwwxENeJFLiRJkiSpIAYsSZIkSSqIAUuSJEmSCuIaLEmS1nEuRpek5mHAkiRpHedidElqHk4RlCRJkqSCGLAkSZIkqSAGLEmSJEkqiGuwJEnqwx657hHmXz+/R21nHTyr2zZjPjyG1oNa17IqSWpedQ9YETEU+D4wGVgCnJyZP+mi3U7A14Fdgc0zM3rTjyRJWr3Wg1oNRJJUoEZMEbwIeBkYDkwFLomIcV20ewW4FvjkWvYjSZIkSXVR14AVEQOAA4FTM3NZZv4W+DlweOe2mXl/Zn4fmLc2/UiSJElSvURm1u+HRewC3JmZ/Tts+zwwMTP3Xc1ntgMe7DhFsNp+IuIo4CiA4cOH73r11VcX9ZUkdeO+vy2tWd9vjoU16/uBN2xUs7532GyHmvWt+luxj3/jjJN58N7X/U2w17YfO45Lpx1TWH9dadb93N8rr1fr3yuO+eu5j3etWce8N/baa6/fZ+aEztvrHbD2BK7LzDd22HYkMDUz21bzma4CVtX9rDBhwoScM2fOWn0PST232/QZNev7ro2PrVnfk0aNqFnfM6fMrFnfqr9m3ceheffzZh3zZh1vcMy74j7etWYd896IiC4DVr3XYC0DBnfaNhioNqYX1Y8kSZIkFabeAesBoCUitu+wbTxdrLOqUz+SJEmSVJi6BqzMfA64ETgjIgZExB7A/sAVndtGWT/gDZXX/SJi42r7kSRJkqR6acRl2j8N9AceB64CjsnMeRExKiKWRcSoSrvRwAv846zUC8D93fVTjy8gSZIkSV2pe8DKzKcy84DMHJCZo1bcHDgzF2bmwMxcWHk9PzOj02NMd/1I67NSqUREFPYolUqN/kqSJEnrlZZGFyCp50qlUrehqK2tDYD29vaa1yNJkqRVNWKKoCRJkiStlwxYkiRJklQQA5YkSZIkFcSAJUmSJEkFMWBJkiRJUkEMWJIkSZJUEAOWJEmSJBXEgCVJkiRJBTFgSZIkSVJBWhpdgKSemXTtpB61m/vE3KraA8ycMrNXNUmSJGlVnsGSJEmSpIIYsCRJkiSpIAYsSZIkSSqIAUuSJEmSCmLAkiRJkqSCGLAkSZIkqSAGLEmSJEkqiAFLkiRJkgpiwJIkSZKkgrQ0ugBJPffIdY8w//r5PWo76+BZ3bYZ8+ExtB7UupZVSZIkaQUDltREWg9qNRBJkiStw5wiKEmSJEkFMWBJkiRJUkEMWJIkSZJUEAOWJEmSJBXEgCUVrFQqERGFPUqlUqO/kiRJknrIgKW1Yph4vVKpRGau8TFx4kQmTpzYbbvMXC/GRJIkqa/wMu1aK6VSqdsA0NbWBkB7e3vN65EkSZIayTNYkiRJklQQA5YkSZIkFcSAJUmSJEkFMWBJkiRJUkEMWJIkSZJUEAOWJEmSJBXEy7RLRTpvx561W7SouvYAo0ZUX48kSZLqyjNYkiRJklQQA5YkSZIkFcSAJUlap5RKJSKisEepVGr0V5Ik9SGuwZIkrVNKpVK3oaitrQ2A9vb2mtcjSVI1PIMlSZIkSQUxYEmSJElSQQxYkiRJklQQA5YkSZIkFcSLXEiS1gmTrp3U47Zzn5hb1WdmTpnZq5okSaqWAUu94oHQ6pVuX8Lpv36qR23jCw9022ba3kMpTR62tmVJkiSpDgxYUsFKk4cZiCRJkvoo12BJkiRJUkEMWJIkSZJUEAOWJEmSJBXENViSpHXKI9c9wvzr5/eo7ayDZ3XbZsyHx8CUtatJkqSeMmBJktYprQe10npQa6PLkCSpV5wiKEmSJEkFMWBJkiRJUkEMWJIkSZJUEAOWJEmSJBXEgCVJkiRJBTFgSZIkSVJBDFiSJEmSVBADlvqU3abPYLfpMxjUOp6IKOwxqHU8u02f0eivJ0mSpAbzRsPqk8Z+6uuNLkGSJEnrIc9gSZIkSVJBDFiSJEmSVBCnCGqtPHLdI8y/fn6P2s46eFa3bcZ8eAxMWbuaJEmSpEYxYGmttB7USutBrY0uQ5IkSVonOEVQkiRJkgpiwJIkVaVUKhV6m4NSqdToryRJUmGcIihJqkqpVOo2FLW1tQHQ3t5e83okSVqXeAZLkiRJkgriGSxJUs+dt2PP2i1aVF17gFEjqq9HkqR1jGewJEmSJKkgBixJkiRJKkjdpwhGxFDg+8BkYAlwcmb+ZDVt/wM4EegP3AAck5kvVd5rB3YHXq00fzQzq5iL0jeUSiVOP/30wvqbNm2aV/yS+rjS7Us4/ddP9ahtfOGBbttM23sopcnD1rYsSZLWCY1Yg3UR8DIwHHgbcHNEzM3MeR0bRcQ+wEnAJGAx8FPg9Mq2FY7LzP+sR9HNyqt9SSpaafIwA5EkSatR1ymCETEAOBA4NTOXZeZvgZ8Dh3fR/GPA9zNzXmb+HfgK8PG6FStJkiRJVar3GqwdgNcys+OckbnAuC7ajqu817Hd8IjYvMO2syJiSUT8LiLaii5WkiRJkqoRmVm/HxaxJ3BdZr6xw7YjgamZ2dap7cPAsZl5W+X1RpSnFrZm5vyIeCfw58q2Q4BvA2/LzIe7+LlHAUcBDB8+fNerr766Fl+vaR1//PEAnH/++T3+zAN/735dRW/tsNkONev7vr8trVnfAG+OhTXr+4E3bFSzvpt1zB3vrt33t6V844yTefDeed037qHtx47jhNPOcsy70Kz7ODjmXXEf75pj/nru411r1jHvjb322uv3mTmh8/Z6B6xdgN9l5iYdtn0OaMvMfTu1nQtMz8xrK683p3xRjGGZ+WQXfd8G3JyZ31pTDRMmTMg5c+as/ZdZj/RmDdakayfVphhg5pSZNet7t+kzatY3wF0bH1uzvifV8B5BzTrmjnfXHPPXcx/vmmP+eu7jXXPMX899vGvNOua9ERFdBqx6TxF8AGiJiO07bBsPdPVn1nmV9zq2e6yrcFWRQBRSpSRJkiT1Ql0DVmY+B9wInBERAyJiD2B/4Ioumv8I+GREvCUiNgO+DPwQICKGRMQ+EdEvIloiYirwbuBXdfkikiRJktSFRtxo+NOU72v1OHAV5XtbzYuIURGxLCJGAVTWXp0DzAIWVB7TKn1sBJwJPEF52uC/Awdk5v11/SaSJEmS1EHd74OVmU8BB3SxfSEwsNO2bwDf6KLtE8A7alTi+uO8Ht53edGi6toD1HB+rSRJktSsGnEGS5IkSZLWSwYsSZIkSSqIAUuSJEmSCmLAkiRJkqSCGLAkSZIkqSAGLEmSJEkqiAFLkiRJkgpS9/tgqb5Kty/h9F8/1aO28YUHum0zbe+hlCYPW9uyJEmSpPWSAWs9V5o8zEAkSZIk1YlTBCVJkiSpIAYsSZIkSSqIAUuSJEmSCmLAkiRJkqSCGLAkSZIkqSAGLEmSJEkqiAFLkiRJkgpiwJIkSZKkghiwJEmSJKkgBixJkiRJKogBS5IkSZIKYsCSJEmSpIIYsCRJkiSpIAYsSZIkSSqIAUuSJEmSCmLAkiRJkqSCGLAkSZIkqSAGLEmSJEkqiAGrwXabPoNBreOJiMIeg1rHs9v0GY3+apIkSVKf09LoAgRjP/X1RpcgSZIkqQCewZIkSZKkghiwJEmSJKkgBixJkiRJKogBS5IkSZIKYsCSJEmSpIIYsCRJkiSpIAYsSZIkSSqIAUuSJEmSCmLAkiRJkqSCGLAkSZIkqSAGLEmSJEkqiAFLkiRJkgpiwJIkSZKkghiwJEmSJKkgBixJkiRJKogBS5IkSZIKYsCSJEmSpIIYsCRJkiSpIAYsSZIkSSqIAUuSJEmSCmLAkiRJkqSCGLAkSZIkqSAGLEmSJEkqiAFLkiRJkgpiwJIkSZKkghiwJEmSJKkgBixJkiRJKogBS5IkSZIKYsCSJEmSpIIYsCRJkiSpIAYsSZIkSSqIAUuSJEmSCmLAkiRJkqSCGLAkSZIkqSAGLEmSJEkqiAFLkiRJkgpiwJIkSZKkghiwJEmSJKkgBixJkiRJKogBS5IkSZIKYsCSJEmSpIIYsCRJkiSpIAYsSZIkSSqIAUuSJEmSClL3gBURQyPipxHxXEQsiIjD1tD2PyLibxHxTET8ICI27k0/kiRJklQPjTiDdRHwMjAcmApcEhHjOjeKiH2Ak4D3AGOAbYDTq+1HkiRJkuqlrgErIgYABwKnZuayzPwt8HPg8C6afwz4fmbOy8y/A18BPt6LfiRJkiSpLiIz6/fDInYB7szM/h22fR6YmJn7dmo7F/hqZl5TeT0MeAIYBozqaT+V944Cjqq83BG4v9Avtn4YBixpdBF9jGNef455fTne9eeY159jXl+Od/055qs3OjO36Lyxpc5FDASe6bTtGWBQD9queD6oyn7IzEuBS6stti+JiDmZOaHRdfQljnn9Oeb15XjXn2Nef455fTne9eeYV6/ea7CWAYM7bRsMLO1B2xXPl1bZjyRJkiTVRb0D1gNAS0Rs32HbeGBeF23nVd7r2O6xzHyyyn4kSZIkqS7qGrAy8zngRuCMiBgQEXsA+wNXdNH8R8AnI+ItEbEZ8GXgh73oRz3jFMr6c8zrzzGvL8e7/hzz+nPM68vxrj/HvEp1vcgFlO9fBfwA2Bt4EjgpM38SEaOAPwNvycyFlbYnACcC/YEbgKMz86U19VPXLyNJkiRJHdQ9YEmSJEnS+qoRNxqWJEmSpPWSAUuSJEmSCmLAUpciIhpdg1RL7uOSJKkWDFjqUlYW53kQWh+Oc/25j9efY90YjrvWd+7jWtd4kQutIiK+CAwFHgF+l5n3RMQGmbm8waWtlyLi7Zn5h8rzSP9B1pz7eP25n9dfRAwGlnb4Q4L7eI1FREtmvtroOvoK9/HGcD/vGc9gaaWI+BnwYWAw8B7gFxHxvsxcHhHuKwWLiJ8AcyLi41A+o+Jf4WrLfbz+3M/rLyL+E7gJuCEivgLgPl47EfG5iBiWma9GxIaNrqcvcB+vP/fz6rgjCoCIGA9snZm7ZeangU9TvrHczRHxgcovLg+KChIRxwI7AacBX46IT4AHn7XkPl5/7uf1FxGXAG8FvgTcCfxzRPwOPACthYi4DDiX8h9rtsjM1zz4rC338fpzP6+eO6FW+DswICImAmTmksw8CzgJuCwi3um0nmJExBuAhcBXgIsp/9I6zYPPmnMfryP38/qLiP7ANsCXM/P/A74JfLT8VtwFKw9AHfcCRMTeQCvls+IPAb/04LO23Mfrz/28dwxYWuFV4GFgt07/YC4CrgIOjIgN/aW19jLzZeC3wC2Z+RRwJXAW5YPPT1baZERs0cAy10cvUt7H3+E+XnuV/bwduNn9vD4y8wVgCfBPEbFRZr6WmYuAA4AXI+J7lXb+IaEY84GvAr8GTqW8rrPLg09/rxTDfbwh/oL7edUMWH1YREyIiLdHRP/MXEx5utQ04OMr2mTm88AC4M2VX2T+0uqlTuP998x8DiAzlwE/Bs4GTo2IAyoHoD+OiH7+wuq9ylkUADLzceD7lPfxj3XY7j5eoE5jvrQyvu7n9fMn4F+AbVdsqOz75wKtETG8UYWtbzLzQeC/Kvv5fOAUygefN3c4+GyrtPX3SnHmAu/HfbwuMvNh4A738+q0NLoANUZE/BTYDngJ2CIijsnMn0XEp4AfRsQA4LbMfKDykYyIfpn5YqNqbmadxnvziDgG+O/MfAYgM5+LiB8DTwI3Uj7b8m7Hu/cq8/TnRcT3K3/1JDNvjPKVp75T2cd/5T5enK7GvCP38+JV1rm9AizPzP/MzK9FxLsoL/7/IDC/ctAzGxgGDAIea1zFza3DeL+Wmd/vuO9m5sMR8WXgTOC6iPgt8KkoX0VzcYNKbnoRsTXw6IqD98w8pzLV+7qI2A/38cJ1MeYvrXjP/bxnDFh9UEQcR3k+7W5AUv6L/lcjYkxmXhwRL1Nel3J4RDwJ7A7s5UFQ76xmvKcD34uIayvTp8jMZRHxNuBZ4J8y888NKrnpRcSFlOflPwW8EBE/6RCyfhgRLwCfAz7qPl6MNY15R+7nxYmIXwCjgf8G/rUSqE7KzAMi4pfADcA5EXE3MAnYBHimYQU3uS7G+32Uz5rc3eFA9KHKOsO5wHHAezzo7L2IuJ7ymaqPUT47C0BmfiAibgeuA85zHy/O6sa8I/fz7hmw+qbBwJ0dDn5Ojoi/AR+PiCWZeW1E3Au8EdgSOCYzH2lUseuB1Y435bnk11fmMI8CPgTs40Fn71X+8jaa8mXY30z5IgsbRMSPO4SsayLiD8DWuI+vtZ6MeYe2rcC/4n6+ViJiArBVZu5cef1F4Grg7Ig4LTM/GBFfp/x75mTK/39/WGY+0aiam9kaxvsLwAWU19Wu8FnKB6g7Z+Y99a51fRERZ1D+/8VngPMj4vjM/FNE+V56mTnZfbxYqxvz1TR3P18DbzTcB0XEQZTXQeyXmfM6bD8NOJryP5YljapvfdPNeB8FjM/MJyNiI6B/Zj7boFLXGxGxLfBMZi6prPMpAWcAP1mx9k3F6umYV9Zo9XM/XzuVNQ9XAjusGN+I2Ay4DCAzD6hs2xLYGHg+M59sRK3rg27G+xXgo5n5QkQMpBy6fpmZdzeo3PVCZTrms5l5ReXs4WDg3ysha+VNhSNiGNAf9/G1tqYx79RuE8oznX7hft41A1YfERH/BmwBzAGeo/wX5A2BCzPzoQ7t/hu4LDMvbUih6wnHu/46jfl9lcW4K95becCfmd+LiCOA/y8z721EreuLKsf848Bsx7wYUb7XzwzgDspjnJXtm1Ke1nNZZpYaV+H6pQfj/Z+Z+ZXKtg0z87WGFbseiYihK6bRdzjg/yww1wsq1EZ3Y77iDKL7+ZoZsPqAiPgZMBL4M+X7R/wFWA68QHkdxA9WHPRExI3ADZl5ZYPKbXqOd/11GPN5wPbAIuCKzPx5hzafAk4E7gU+CLy14xlFVccxr7+I+Cjw/zLzkcoZ73+jvO7k58DlHQ76/wMYl5mfaly1zc/xrr+OY95hW0tmvlp5/kvKF7E4mPL/Hv0z86SGFLue6MWYD8jMLzak2CbiGqz1XETsQfnqdTtn+VKaE4B9KS/qfwZ4Gfhu5R9QP2Avyqd91QuOd/11MeZvp3xDxC9W/sL2U4DM/M8oX3FqD2AXD/R7zzGvvyhffXFf4JqI+FqWr+R1FTCGcnjdHPh6pflWwMYdp1GpOo53/XUa87Mz8y+VsyWvrjhbUllbeDXls4abAf/c0KKbnGNeO94Ha/33LOUD+2GVX/5zgB8B/wUE8D+UF+q+h/JfoffKf1y2WtVzvOuv85j/gfI93WYAn4yIXQEi4ijKB0aTMnNuw6pdPzjmdRQRH6Z80aGvABsBJ0bE9pX1JmdTnqK5T0Q8GBHfp7yW9lwP9nvH8a6/Lsb8pIjYpsOUtNciYsVJgTmVNm/LzNkNKrnpOea15RTB9VxEjADupvzL/5sdtu8MfInyfYAuq0x/eM3/g1g7jnf9rWHMx1FeA3RLZczfAbyYmf/bmErXH455fUXEVpTPgv8ceB9wCOUpx+dk5oOVBecDKU/hWUb5qqX3N6reZud4199qxvx5oONZlYzyPd7+C3hX5Q+Y6iXHvLYMWH1AROwPXAEcW7kyzIp/NGcB44EPeqBfHMe7/roZ852B/VfMJ1cxHPP6ioiNMvOVyvP9gCmUD/q/luV70rwFeCgzX25knesLx7v+1jDmZ1UO+HfKzHsiYnhmeiPhAjjmteMarL7h55QXmn8zIgZk5ncq25cC/0d56pqK43jX35rG/G+Ub/CsYjnmdZSZr6wIsZn584gI4CDgmCjfR28isDfle+tpLTne9beGMT+2MuaTIqLNA/3iOOa14xmsPqLyD+Vg4HvAnZT/QvFuymuA/tjI2tZHjnf9Oeb155jX34qDocrzdwE/BN5EeZ2b03cK5njXn2Nef4558QxYfUxE7EB5zu0bgDsy88EGl7Rec7zrzzGvP8e8vjpMxzye8pXsxmfmPQ0ua73leNefY15/jnmxnCLYx2T5inVeta5OHO/6c8zrzzGvr8pB0EDgrcBuHgTVluNdf455/TnmxfIMliRJTajjAnXVnuNdf455/TnmxTBgSZIkSVJBvNGwJEmSJBXEgCVJkiRJBTFgSZIkSVJBDFiSJEmSVBADliRJkiQVxIAlSZIkSQUxYEmSJElSQf5/OaaICagzZNkAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 864x576 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "    \n",
+    "def plot(labels, averages, stds, xticks):\n",
+    "    n_bars = len(xticks)\n",
+    "    n_labels = len(labels)\n",
+    "    \n",
+    "    b_width = (1/n_labels)*0.75\n",
+    "    fig = plt.subplots(figsize = (12, 8))\n",
+    "\n",
+    "    x_pos = np.arange(n_bars)\n",
+    "    for i, avg in enumerate(averages):\n",
+    "        plt.bar(x_pos, avg, yerr = stds[i], width = b_width, label = labels[i],\n",
+    "                align='center', alpha=0.9, ecolor='black', capsize=b_width*40)\n",
+    "        \n",
+    "        x_pos = [x + b_width for x in x_pos]\n",
+    "    \n",
+    "    plt.ylabel('probability', fontsize = 15)\n",
+    "    plt.xticks([x + b_width*(n_labels-1)/2 for x in range(n_bars)], xticks, fontsize=12, rotation=45)\n",
+    "    plt.grid(b=True, axis='y')\n",
+    "    plt.yticks(fontsize=12)\n",
+    "    plt.legend(fontsize=12)\n",
+    "    plt.tight_layout()\n",
+    "    #plt.savefig('bar_plot_with_error_bars.png')\n",
+    "    plt.show()\n",
+    "\n",
+    "def plot_result(ideal, probs):\n",
+    "    labels   = [ key for key, val in probs.items() ]\n",
+    "    averages = [ np.mean(prob, axis=0) for key, prob in probs.items() ]\n",
+    "    stds     = [ np.std(prob , axis=0) for key, prob in probs.items() ]\n",
+    "    xticks   = ['{:0{}b}'.format(p, n) for p, a in enumerate(input_state)]\n",
+    "\n",
+    "    labels.insert(0, 'ideal')\n",
+    "    averages.insert(0, ideal)\n",
+    "    stds.insert(0, [0]*len(ideal))\n",
+    "\n",
+    "    plot(labels, averages, stds, xticks)\n",
+    "\n",
+    "for s in range(1, n+1):\n",
+    "    print('s =', s)\n",
+    "    plot_result(ideal, probs[s])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "formal-mozambique",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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.7.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/example/schmidt/Cnots and depth.ipynb b/example/schmidt/Cnots and depth.ipynb
index 8b04c795..14fb6df0 100644
--- a/example/schmidt/Cnots and depth.ipynb
+++ b/example/schmidt/Cnots and depth.ipynb
@@ -15,14 +15,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from qiskit import IBMQ, QuantumCircuit, ClassicalRegister, QuantumRegister, execute, Aer, transpile\n",
-    "\n",
-    "backend = Aer.get_backend('qasm_simulator')\n",
-    "\n",
-    "shots   = 8192\n",
-    "backend_name = backend.name()\n",
-    "backend_config = backend.configuration()\n",
-    "backend_qubits = backend_config.n_qubits"
+    "from qiskit import IBMQ, QuantumCircuit, ClassicalRegister, QuantumRegister, execute, Aer, transpile"
    ]
   },
   {
@@ -35,7 +28,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 2,
    "id": "painted-clone",
    "metadata": {},
    "outputs": [
@@ -44,47 +37,37 @@
      "output_type": "stream",
      "text": [
       "\n",
-      "low rank approx. = 1\n",
+      "low rank approx. = 2\n",
       "\n",
       "n = 3\n",
-      "\t isometry\t CNOTs = 1\t depth = 3\n",
-      "\t unitary\t CNOTs = 4\t depth = 9\n",
+      "\t knill\t\t CNOTs = 4\t depth = 9\n",
+      "\t ccd\t\t CNOTs = 4\t depth = 9\n",
+      "\t full rank\t CNOTs = 4\t depth = 9\n",
       "\n",
       "n = 4\n",
-      "\t isometry\t CNOTs = 2\t depth = 3\n",
-      "\t unitary\t CNOTs = 9\t depth = 11\n",
+      "\t knill\t\t CNOTs = 7\t depth = 9\n",
+      "\t ccd\t\t CNOTs = 7\t depth = 9\n",
+      "\t full rank\t CNOTs = 9\t depth = 11\n",
       "\n",
       "n = 5\n",
-      "\t isometry\t CNOTs = 5\t depth = 9\n",
-      "\t unitary\t CNOTs = 30\t depth = 53\n",
+      "\t knill\t\t CNOTs = 43\t depth = 75\n",
+      "\t ccd\t\t CNOTs = 14\t depth = 23\n",
+      "\t full rank\t CNOTs = 30\t depth = 53\n",
       "\n",
       "n = 6\n",
-      "\t isometry\t CNOTs = 8\t depth = 9\n",
-      "\t unitary\t CNOTs = 54\t depth = 57\n",
+      "\t knill\t\t CNOTs = 55\t depth = 53\n",
+      "\t ccd\t\t CNOTs = 21\t depth = 23\n",
+      "\t full rank\t CNOTs = 54\t depth = 57\n",
       "\n",
       "n = 7\n",
-      "\t isometry\t CNOTs = 13\t depth = 11\n",
-      "\t unitary\t CNOTs = 150\t depth = 243\n",
+      "\t knill\t\t CNOTs = 140\t depth = 152\n",
+      "\t ccd\t\t CNOTs = 36\t depth = 51\n",
+      "\t full rank\t CNOTs = 150\t depth = 243\n",
       "\n",
       "n = 8\n",
-      "\t isometry\t CNOTs = 18\t depth = 11\n",
-      "\t unitary\t CNOTs = 253\t depth = 247\n",
-      "\n",
-      "n = 9\n",
-      "\t isometry\t CNOTs = 78\t depth = 121\n",
-      "\t unitary\t CNOTs = 661\t depth = 1039\n",
-      "\n",
-      "n = 10\n",
-      "\t isometry\t CNOTs = 174\t depth = 165\n",
-      "\t unitary\t CNOTs = 1090\t depth = 1080\n",
-      "\n",
-      "n = 11\n",
-      "\t isometry\t CNOTs = 147\t depth = 165\n",
-      "\t unitary\t CNOTs = 2770\t depth = 4345\n",
-      "\n",
-      "n = 12\n",
-      "\t isometry\t CNOTs = 108\t depth = 57\n",
-      "\t unitary\t CNOTs = 4474\t depth = 4346\n"
+      "\t knill\t\t CNOTs = 168\t depth = 150\n",
+      "\t ccd\t\t CNOTs = 51\t depth = 51\n",
+      "\t full rank\t CNOTs = 253\t depth = 247\n"
      ]
     }
    ],
@@ -95,11 +78,11 @@
     "from qclib.state_preparation.schmidt import initialize\n",
     "from qclib.state_preparation.schmidt import initialize_original\n",
     "\n",
-    "r = 1\n",
+    "r = 2\n",
     "print('\\nlow rank approx. =', r)\n",
     "\n",
     "result = {}\n",
-    "for n in range(3, 13):\n",
+    "for n in range(3, 9):\n",
     "    N = 2**n\n",
     "    print('\\nn =', n)\n",
     "        \n",
@@ -108,21 +91,28 @@
     "\n",
     "    result[n] = []\n",
     "    \n",
-    "    circuit1 = initialize(input_state, low_rank = r)\n",
-    "    circuit2 = initialize_original(input_state)\n",
+    "    circuit1 = initialize(input_state, low_rank = r, isometry_scheme='knill')\n",
+    "    circuit2 = initialize(input_state, low_rank = r, isometry_scheme='ccd')\n",
+    "    circuit3 = initialize(input_state)\n",
     "    \n",
     "    transpiled_circuit1 = transpile(circuit1, basis_gates=['u1','u2','u3', 'cx'], optimization_level=3)\n",
     "    transpiled_circuit2 = transpile(circuit2, basis_gates=['u1','u2','u3', 'cx'], optimization_level=3)\n",
+    "    transpiled_circuit3 = transpile(circuit3, basis_gates=['u1','u2','u3', 'cx'], optimization_level=3)\n",
+    "    \n",
     "    n_cx1 = transpiled_circuit1.count_ops()['cx']\n",
     "    n_dp1 = transpiled_circuit1.depth()\n",
     "    n_cx2 = transpiled_circuit2.count_ops()['cx']\n",
     "    n_dp2 = transpiled_circuit2.depth()\n",
+    "    n_cx3 = transpiled_circuit3.count_ops()['cx']\n",
+    "    n_dp3 = transpiled_circuit3.depth()\n",
     "\n",
-    "    result[n].append(['isometry', n_cx1, n_dp1])\n",
-    "    result[n].append(['unitary', n_cx2, n_dp2])\n",
+    "    result[n].append(['knill', n_cx1, n_dp1])\n",
+    "    result[n].append(['ccd', n_cx2, n_dp2])\n",
+    "    result[n].append(['fullrank', n_cx3, n_dp3])\n",
     "\n",
-    "    print('\\t {0}\\t CNOTs = {1}\\t depth = {2}'.format('isometry', n_cx1, n_dp1))\n",
-    "    print('\\t {0}\\t CNOTs = {1}\\t depth = {2}'.format('unitary', n_cx2, n_dp2))\n",
+    "    print('\\t {0}\\t\\t CNOTs = {1}\\t depth = {2}'.format('knill', n_cx1, n_dp1))\n",
+    "    print('\\t {0}\\t\\t CNOTs = {1}\\t depth = {2}'.format('ccd', n_cx2, n_dp2))\n",
+    "    print('\\t {0}\\t CNOTs = {1}\\t depth = {2}'.format('full rank', n_cx3, n_dp3))\n",
     "        "
    ]
   },
diff --git a/qclib/isometry.py b/qclib/isometry.py
index b875f86b..04659f43 100644
--- a/qclib/isometry.py
+++ b/qclib/isometry.py
@@ -12,102 +12,272 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-import numpy as np
-import scipy
-from qiskit import QuantumCircuit, QuantumRegister
-from qclib.unitary import unitary
-from qclib.state_preparation.bdsp import initialize
-
-def decompose(V, scheme='knill'):
-    v = np.copy(V)
-    # change a row vector to a column vector (in the case of state preparation).
-    if len(v.shape) == 1:
-        v = v.reshape(v.shape[0], 1)
-
-    N = v.shape[0]
-    M = v.shape[1]
-    n = np.log2(N)
-    m = np.log2(M)
-    
-    _check_isometry(v, n, m)
-
-    n = int(n)
-    m = int(m)
-    
-    if (N == M): # The isometry v is already unitary.
-        U = v
-    else:        # The isometry v is extended to a unitary maximizing the numbers of eigenvalues with complex argument equal to zero.
-        
-        # The complex conjugate of the null space is the transformation that generates the state |v> from |0>.
-        # V=UI_{2^n,2^m} => U^-1 V = I_{2^n,2^m}
-        w = np.conj(scipy.linalg.null_space(v.T)) # The transposition was removed because it would be nullified in U[:, M:] = w.T .
-        U = np.zeros((N,N), dtype=complex)
-        U[:, :M] = v
-        U[:, M:] = w                              # The transposition was removed because it nullified the transposition of the conjugate complex above (which was removed).
-    
-    if (scheme == 'csd'):
-        pass
-    elif (scheme == 'cbc'):
-        pass
 
-    return _knill(U, n)    
+"""
+Implements the decomposition of isometries using the methods
+defined at https://arxiv.org/abs/1501.06911.
+"""
 
 
+import numpy as np
+import scipy
+import qiskit.quantum_info as qi
+from qiskit import QuantumCircuit, QuantumRegister
+from qiskit.extensions.quantum_initializer.uc import UCGate
 
+# pylint: disable=maybe-no-member
 
-def _knill(U, n):
-    if n < 2:
-        raise ValueError(
-            "Knill decomposition does not work on a 1 qubit isometry (N=2)."
-        )
+def decompose(isometry, scheme='ccd'):
+    """
+    Decompose an isometry from m to n qubits.
+    In particular, it decomposes unitaries on n qubits (m=n) or prepare a
+    quantum state on n qubits (m=0).
+    `https://arxiv.org/abs/1501.06911`_.
+    Args:
+        isometry (list): an isometry from m to n qubits (n>=2 and m<=n), i.e., a
+            complex 2^n×2^m array with orthonormal columns.
+        scheme (str): method to decompose the isometry ('knill', 'ccd', 'csd').
+            Default is scheme='ccd'.
+    Returns:
+        QuantumCircuit: a quantum circuit with the isometry attached.
+    Raises:
+        ValueError: if isometry is not valid.
+    """
 
-    from qclib.state_preparation.schmidt import initialize_original as schmidt
+    iso = isometry.astype(complex)
+    if len(iso.shape) == 1:
+        iso = iso.reshape(iso.shape[0], 1) # change a row vector to a column
+                                           # vector (in the case of state
+                                           # preparation).
+    lines = iso.shape[0]
+    cols = iso.shape[1]
+    log_lines = np.log2(lines)
+    log_cols = np.log2(cols)
 
-    A = QuantumRegister(n)
-    circuit = QuantumCircuit(A)
-        
-    eigval, eigvec = np.linalg.eig(U)
-    arg = np.angle(eigval)
-        
-    for i in range(2**n): # The eigenvalues are not necessarily ordered.
-        if (np.abs(arg[i]) > 10**-15):
-            state = eigvec[:,i]
-            
-            circuit.compose( schmidt(state).inverse(), A, inplace=True )
+    _check_isometry(iso, log_lines, log_cols)
 
-            circuit.x(list(range(n)))
-            circuit.mcp(arg[i], list(range(n-1)), n-1)
-            circuit.x(list(range(n)))
+    log_lines = int(log_lines)
+    log_cols = int(log_cols)
 
-            circuit.compose( schmidt(state)          , A, inplace=True )
-    
-    return circuit
+    if scheme == 'csd':
+        pass
+    if scheme == 'ccd':
+        return _ccd(iso, log_lines, log_cols)
 
+    return _knill(iso, log_lines, log_cols)
 
 
 
+#   General
 
 
 
-def _check_isometry(V, n, m):
-    if not n.is_integer() or n < 0:
+def _check_isometry(iso, log_lines, log_cols):
+    if not log_lines.is_integer() or log_lines < 0:
         raise ValueError(
             "The number of rows of the isometry is not a non negative power of 2."
         )
-    if not m.is_integer() or m < 0:
+    if not log_cols.is_integer() or log_cols < 0:
         raise ValueError(
             "The number of columns of the isometry is not a non negative power of 2."
         )
-    if m > n:
+    if log_cols > log_lines:
+        raise ValueError(
+            "The input matrix has more columns than rows."
+        )
+    if not _is_isometry(iso, log_cols):
         raise ValueError(
-            "The input matrix has more columns than rows and hence it can't be an isometry."
+            "The input matrix has non orthonormal columns."
         )
-    if not is_isometry(V, m):
+
+def _is_isometry(iso, log_cols):
+    identity = np.conj(iso.T).dot(iso)
+    return np.allclose(identity, np.eye(int(2**log_cols)))
+
+
+
+#   Knill
+
+
+
+def _knill(iso, log_lines, log_cols):
+    if log_lines < 2:
         raise ValueError(
-            "The input matrix has non orthonormal columns and hence it is not an isometry."
+            "Knill decomposition does not work on a 1 qubit isometry (N=2)."
         )
 
-def is_isometry(V, m):
-    I = np.conj(V.T).dot(V)
-    return np.allclose(I, np.eye(int(2**m)))
+    if log_lines == log_cols:      # The isometry v is already unitary.
+        unitary = iso
+    else:                          # The isometry v is extended to a unitary maximizing
+                                   # the numbers of eigenvalues with complex argument equal
+                                   # to zero.
+        null_space = np.conj(scipy.linalg.null_space(iso.T)) # The complex conjugate of the null
+                                                             # space is the transformation that
+                                                             # generates the state |v> from |0>.
+                                                             # V=UI_{2^n,2^m} => U^-1 V=I_{2^n,2^m}
+                                                             # The transposition was removed
+                                                             # because it would be nullified in
+                                                             # U[:, M:] = W.T .
+        unitary = np.zeros((2**log_lines, 2**log_lines), dtype=complex)
+        unitary[:, :2**log_cols] = iso
+        unitary[:, 2**log_cols:] = null_space                # The transposition was removed
+                                                             # because it nullified the
+                                                             # transposition of the conjugate
+                                                             # complex above (which was removed).
+
+    eigval, eigvec = np.linalg.eig(unitary)
+    arg = np.angle(eigval)
+
+    reg = QuantumRegister(log_lines)
+    circuit = QuantumCircuit(reg)
+
+    from qclib.state_preparation.schmidt import initialize as schmidt # pylint: disable=import-outside-toplevel
+
+    for i in range(2**log_lines):                            # The eigenvalues are not necessarily
+                                                             # ordered.
+        if np.abs(arg[i]) > 10**-15:
+            state = eigvec[:,i]
+
+            circuit.compose( schmidt(state).inverse(), reg, inplace=True )
+
+            circuit.x(list(range(log_lines)))
+            circuit.mcp(arg[i], list(range(log_lines-1)), log_lines-1)
+            circuit.x(list(range(log_lines)))
+
+            circuit.compose( schmidt(state)          , reg, inplace=True )
+
+    return circuit
+
+
+
+#   Column-by-column
+
+
+
+def _ccd(iso, log_lines, log_cols):
+    reg = QuantumRegister(log_lines)
+    circuit = QuantumCircuit(reg)
+
+    for k in range(2**log_cols): # iteration through columns (starting with the second column k=1).
+        g_k = _g_k(iso, log_lines, k)
+        g_k.name = 'G'+str(k)
+
+        circuit.append( g_k, reg )
+
+    if log_cols > 0:             # It clears the phases, as explained in the last sentence of the
+                                 # first column on page 17.
+        phases = np.angle( np.diagonal( iso[:2**log_cols,:2**log_cols] ) )
+        diag = np.exp( -1j*phases )
+        circuit.diagonal( diag.tolist(), list(range(log_cols)) )
+
+    return circuit.inverse()
+
+def _g_k(iso, log_lines, col_index):
+    g_k = QuantumCircuit(log_lines)                   # Gate G columns index k, to be created.
+    k_bin = '{:0{}b}'.format(col_index, log_lines)    # Binary representation of column index k.
+    for i in range(log_lines):                        # G_k's subgate bit index i (s in the paper).
+        target = log_lines - i - 1
+        control = list(range(target))
+        ancilla = list(range(target+1, log_lines))
+
+        if _k_s(col_index, i) == 0 and _b(col_index, i+1) != 0:      # Condition defined in the
+                                                                     # first paragraph of the
+                                                                     # second column on page 16.
+            unitary = _mc_unitary(iso, col_index, i)                 # Generates single-qubit gate
+                                                                     # matrix for the
+                                                                     # multicontrolled operation.
+            mcg = _mc_gate(unitary, log_lines, control+ancilla, target, k_bin)
+            mcg = mcg.reverse_bits()                                 # Qiskit little-endian.
+            _update_isometry(iso, mcg)                               # Updates the isometry matrix.
+            g_k.compose( mcg, list(range(log_lines)), inplace=True ) # Append "mcg" to the
+                                                                     # operator "G_k".
+
+        unitaries = _uc_unitaries(iso, log_lines, col_index, i)      # Generates single-qubit gates
+                                                                     # matrices for the uniformly
+                                                                     # controlled operation.
+        ucg = _uc_gate(unitaries, log_lines, control, target)
+        ucg = ucg.reverse_bits()                                     # Qiskit little-endian.
+        _update_isometry(iso, ucg)                                   # Updates the isometry matrix.
+        g_k.compose( ucg, list(range(log_lines)), inplace=True )     # Append "ucg" to the
+                                                                     # operator "G_k".
+    return g_k
+
+def _mc_gate(unitary, n_qubits, control, target, k_bin):
+    gate = QuantumCircuit(n_qubits)
+
+    controls = []
+    for i in control:
+        if k_bin[i] == '1':
+            controls.append(i)
+
+    unitaries = [np.identity(2) for i in range(2**len(controls))]
+    unitaries[-1] = unitary
+    ucg = UCGate(unitaries, up_to_diagonal=True)
+    gate.append( ucg, [target] + controls[::-1])
+
+    return gate
+
+def _uc_gate(unitaries, n_qubits, control, target):
+    gate = QuantumCircuit(n_qubits)
+
+    if len(control) > 0 :
+        gate.uc(unitaries, control[::-1], target, up_to_diagonal=True) # "control" is reversed due
+                                                                       # to UCGate implementation.
+    else:
+        gate.unitary(unitaries[0], target)                             # UCGate does not work with
+                                                                       # target only.
+    return gate
+
+def _update_isometry(iso, gate):
+    gate_matrix = qi.Operator(gate).data # "gate" matrix representation.
+    iso[:,:] = gate_matrix @ iso         # Updates isometry.
+
+def _mc_unitary(iso, col_index, bit_index):
+    col = iso[:,col_index]
+    idx1 =  2*_a(col_index, bit_index+1)    * 2**bit_index + _b(col_index, bit_index+1)
+    idx2 = (2*_a(col_index, bit_index+1)+1) * 2**bit_index + _b(col_index, bit_index+1)
+
+    return _unitary([ [ col[idx1] ], [ col[idx2] ] ], basis=0)
+
+def _uc_unitaries(iso, n_qubits, col_index, bit_index):
+    start = _a(col_index, bit_index+1) + 1
+    if _b(col_index, bit_index+1) == 0:
+        start = _a(col_index, bit_index+1)
+
+    gates = []
+    for i in range(start):
+        gates.append( np.identity( 2 ) )
+
+    col = iso[:, col_index]
+    for i in range(start, 2**(n_qubits-bit_index-1)):
+        idx1 =  2*i    * 2**bit_index + _b(col_index, bit_index)
+        idx2 = (2*i+1) * 2**bit_index + _b(col_index, bit_index)
+
+        gates.append( _unitary([ [ col[idx1] ], [ col[idx2] ] ], basis=_k_s(col_index, bit_index)) )
+
+    return gates
+
+def _unitary(iso, basis=0): #  Lemma2 of https://arxiv.org/abs/1501.06911
+    iden = np.identity(2)
+
+    iso_norm = np.linalg.norm(iso, axis=0)[0]
+    psi = iso / iso_norm
+
+    psi_dagger = np.conj(psi.T)
+
+    phi = -psi_dagger.dot( iden[:,[1]] ) * iden[:,[0]] + psi_dagger.dot( iden[:,[0]] ) * iden[:,[1]]
+    phi_dagger = np.conj(phi.T)
+
+    unitary = np.kron( iden[:,[basis]], psi_dagger ) + np.kron( iden[:,[-(basis+1)]], phi_dagger)
+
+    return unitary
+
+def _a(col_index, bit_index):
+    return col_index // 2**bit_index
+
+def _b(col_index, bit_index):
+    return col_index - (_a(col_index, bit_index) * 2**bit_index)
 
+def _k_s(col_index, bit_index):
+    return int((col_index & 2**bit_index) / 2**bit_index) # Returns the bit value at bit_index
+                                                          # of col_index (k in the paper).
diff --git a/qclib/state_preparation/schmidt.py b/qclib/state_preparation/schmidt.py
index a5c35f2a..b4d04cf9 100755
--- a/qclib/state_preparation/schmidt.py
+++ b/qclib/state_preparation/schmidt.py
@@ -12,17 +12,21 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-""" State preparation with Schmidt decomposition """
+"""
+Implements the state preparation
+defined at https://arxiv.org/abs/1003.5760.
+"""
 
 import numpy as np
+import deprecation
 from qiskit import QuantumCircuit, QuantumRegister
 from qclib.state_preparation.mottonen import initialize as mottonen
-from qclib.unitary import unitary
-from qclib.isometry import decompose
+from qclib.unitary import unitary as decompose_unitary
+from qclib.isometry import decompose as decompose_isometry
 
 # pylint: disable=maybe-no-member
 
-def initialize(state, low_rank=0):
+def initialize(state_vector, low_rank=0, isometry_scheme='ccd', unitary_scheme='qsd'):
     """ State preparation using Schmidt decomposition arXiv:1003.5760
 
     For instance, to initialize the state a|0> + b|1>
@@ -31,14 +35,24 @@ def initialize(state, low_rank=0):
 
     Parameters
     ----------
-    state: list of int
+    state_vector: list of int
         A unit vector representing a quantum state.
         Values are amplitudes.
 
     low_rank: int
         ``state`` low-rank approximation (1 <= ``low_rank`` < 2**(n_qubits//2)).
-        If ``low_rank`` is not in the valid range, it will be ignored and the effective
-        ``rank`` will be used.
+        If ``low_rank`` is not in the valid range, it will be ignored.
+
+    isometry_scheme: string
+        Scheme used to decompose isometries.
+        Possible values are ``'knill'`` and ``'ccd'`` (column-by-column decomposition).
+        Default is ``isometry_scheme='ccd'``.
+
+    unitary_scheme: string
+        Scheme used to decompose unitaries.
+        Possible values are ``'csd'`` (cosine-sine decomposition) and ``'qsd'`` (quantum
+        Shannon decomposition).
+        Default is ``unitary_scheme='qsd'``.
 
     Returns
     -------
@@ -46,101 +60,82 @@ def initialize(state, low_rank=0):
         QuantumCircuit to initialize the state.
     """
 
-    state_copy = np.copy(state)
-
-    hsa, hsb, n_a, n_b = _init_params(state_copy)
-
-    state_copy.shape = (n_a, n_b)
-
-    u_matrix, diag, v_matrix = np.linalg.svd(state_copy)
-
-    if 1 <= rank <= n_a:
-        # This code is equivalent to the commented one below.
-        # It is more compact, but the one below is easier to read.
-        diag = diag[:rank]
-        u_matrix = u_matrix[:, :rank]
-        v_matrix = v_matrix[:rank, :]
-        if not np.log2(rank).is_integer(): # To use isometries, the rank needs to be a power of 2.
-            state_copy = u_matrix @ np.diag(diag) @ v_matrix
-
-            # s = np.zeros((n, m), dtype=complex)
-            # for i in range(rank):
-            #     s += d[i] * np.outer(u.T[i], v[i])
-
-    r = n_qubits % 2
-
-    n = int(2**(n_qubits//2))
-    m = int(2**(n_qubits//2 + r))
-    s.shape = (n, m)
-
-    u, d, v = np.linalg.svd(s)
-    
-    
-    rank = sum(j > 10**-16 for j in d)     # Effective rank ( rank \in {1, 2, ..., n} ).
-    if (low_rank > 0 and low_rank < rank):
-        rank = low_rank                    # Low-rank approximation
-    
-    k = int(2**np.ceil(np.log2(rank)))     # To use isometries, the rank needs to be a power of 2.
-    
-    u = u[:,:k]                            # Matrix u can be a unitary (k=n) or isometry (k<n).
-    v = v[:k,:]                            # If n<m, v.T is always an isometry of log2(k) to log2(m) (k<=n).
-                                           # If n=m, v.T can be a unitary (k=n) or isometry (k<n).
-    d = d[:k]                              # The length of the state vector needs to be a power of 2.
-    d[rank:] = np.zeros(k-rank)            # If k>rank, zeroes out the additional elements.
-    
-    if (rank == 1):
-        d = np.concatenate((d, [0]))       # The length of the state vector needs to be a power of 2 and >1.
-    
-
-    d = d / np.linalg.norm(d)
-    
-    A = QuantumRegister(n_qubits//2 + r)
-    B = QuantumRegister(n_qubits//2)
-
-    if len(diag) > 2:
-        circ = initialize(diag)
+    state = np.copy(state_vector)
+
+    size = len(state)
+    n_qubits = int(np.log2(size))
+
+    odd = n_qubits % 2
+
+    lines = int(2**(n_qubits//2))
+    cols = int(2**(n_qubits//2 + odd))
+    state.shape = (lines, cols)
+
+    svd_u, singular_values, svd_v = np.linalg.svd(state)
+
+    rank = lines
+    if low_rank > 0 and low_rank < rank:
+        e_rank = sum(j > 10**-16 for j in singular_values) # Effective rank.
+        if low_rank < e_rank:
+            e_rank = low_rank                              # Low-rank approximation
+
+        rank = int(2**np.ceil(np.log2(e_rank)))         # To use isometries, the rank needs to be
+                                                        # a power of 2.
+
+        svd_u = svd_u[:,:rank]                          # Matrix u can be a unitary (k=n) or
+                                                        # isometry (k<n).
+        svd_v = svd_v[:rank,:]                          # If n<m, v.T is always an isometry of
+                                                        # log2(k) to log2(m) (rank<=n).
+                                                        # If n=m, v.T can be a unitary (k=n) or
+                                                        # isometry (k<n).
+        singular_values = singular_values[:rank]        # The length of the state vector needs to
+                                                        # be a power of 2.
+        singular_values[e_rank:] = np.zeros(rank-e_rank) # If k>rank, zeroes out
+                                                                         # the additional elements.
+        if e_rank == 1:
+            singular_values = np.concatenate((singular_values, [0]))     # The length of the state
+                                                                         # vector needs to be a
+                                                                         # power of 2 and >1.
+    singular_values = singular_values / np.linalg.norm(singular_values)
+
+    reg_a = QuantumRegister(n_qubits//2 + odd)
+    reg_b = QuantumRegister(n_qubits//2)
+
+    circuit = QuantumCircuit(reg_a, reg_b)
+
+    if len(singular_values) > 2:
+        circ = initialize(singular_values)
     else:
-        circ = mottonen(d)
-        
-    circuit.compose(circ, B[:int( np.log2( len(d) ) )], inplace=True)
+        circ = mottonen(singular_values)
 
-    for j in range(int( np.log2( k ) )):           # Entangles only the necessary qubits, according to rank.
-        circuit.cx(B[j], A[j])
+    circuit.compose(circ, reg_b[:int( np.log2( len(singular_values) ) )], inplace=True)
 
-    def encode(U, reg):                            # Encodes the data using the most appropriate method:
-        if (U.shape[1] == 1):                      #   State preparation.
-            if (U.shape[0] > 2):
-                gate_u = initialize(U[:,0])
+    for j in range(int( np.log2( rank ) )):                # Entangles only the necessary qubits,
+        circuit.cx(reg_b[j], reg_a[j])                     # according to rank.
+
+    def _encode(unitary, circuit, reg): # Encodes the data using the most appropriate method:
+        if unitary.shape[1] == 1:                   # State preparation.
+            if unitary.shape[0] > 2:
+                gate_u = initialize(unitary[:,0])
             else:
-                gate_u = mottonen(U[:,0])
-        elif (U.shape[0] > U.shape[1]):            #   Isometry decomposition.
-            gate_u = decompose(U, scheme='knill') 
-        else:                                      #   Unitary decomposition.
-            gate_u = unitary(U, 'qsd')
-        
-        circuit.compose(gate_u, reg, inplace=True) # Apply gate U to the register.
+                gate_u = mottonen(unitary[:,0])
+        elif unitary.shape[0] > unitary.shape[1]:   # Isometry decomposition.
+            gate_u = decompose_isometry(unitary, scheme=isometry_scheme)
+        else:                                       # Unitary decomposition.
+            gate_u = decompose_unitary(unitary, decomposition=unitary_scheme)
+
+        circuit.compose(gate_u, reg, inplace=True)  # Apply gate U to the register.
 
-    encode(u_matrix, hsb)
-    encode(v_matrix.T, hsa)
+    _encode(svd_u, circuit, reg_b)
+    _encode(svd_v.T, circuit, reg_a)
 
     return circuit
 
 
-def _init_params(state_copy):
-    size = len(state_copy)
-    n_qubits = int(np.log2(size))
-    is_odd = n_qubits % 2
-    n_a = int(2 ** (n_qubits // 2))
-    n_b = int(2 ** (n_qubits // 2 + is_odd))
-    hsa = QuantumRegister(n_qubits // 2 + is_odd)
-    hsb = QuantumRegister(n_qubits // 2)
-    return hsa, hsb, n_a, n_b
-
-
-#import deprecation
-#@deprecation.deprecated(deprecated_in="0.0.7",
-#                        details="Use the initialize function instead")
-def initialize_original(state):
+
+@deprecation.deprecated(deprecated_in="0.0.7",
+                        details="Use the initialize function instead")
+def initialize_original(state_vector):
     """ State preparation using Schmidt decomposition arXiv:1003.5760.
         This function implements the original algorithm as defined in arXiv:1003.5760.
         It is kept here for didactic reasons.
@@ -153,7 +148,7 @@ def initialize_original(state):
 
     Parameters
     ----------
-    state: list of int
+    state_vector: list of int
         A unit vector representing a quantum state.
         Values are amplitudes.
 
@@ -163,37 +158,37 @@ def initialize_original(state):
         QuantumCircuit to initialize the state.
     """
 
-    state_copy = np.copy(state)
+    state = np.copy(state_vector)
 
-    size = len(state_copy)
+    size = len(state)
     n_qubits = int(np.log2(size))
 
-    is_odd = n_qubits % 2
+    odd = n_qubits % 2
 
-    state_copy.shape = (int(2**(n_qubits//2)), int(2**(n_qubits//2 + is_odd)))
+    state.shape = (int(2**(n_qubits//2)), int(2**(n_qubits//2 + odd)))
 
-    u_matrix, diag, v_matrix = np.linalg.svd(state_copy)
+    svd_u, singular_values, svd_v = np.linalg.svd(state)
 
-    diag = diag / np.linalg.norm(diag)
+    singular_values = singular_values / np.linalg.norm(singular_values)
 
-    hsa = QuantumRegister(n_qubits//2 + is_odd)
-    hsb = QuantumRegister(n_qubits//2)
+    reg_a = QuantumRegister(n_qubits//2 + odd)
+    reg_b = QuantumRegister(n_qubits//2)
 
-    circuit = QuantumCircuit(hsa, hsb)
+    circuit = QuantumCircuit(reg_a, reg_b)
 
-    if len(diag) > 2:
-        circ = initialize_original(diag)
+    if len(singular_values) > 2:
+        circ = initialize_original(singular_values)
     else:
-        circ = mottonen(diag)
-    circuit.compose(circ, hsb, inplace=True)
+        circ = mottonen(singular_values)
+    circuit.compose(circ, reg_b, inplace=True)
 
     for k in range(int(n_qubits//2)):
-        circuit.cx(B[k], A[k])
-    
-    gate_u = unitary(u, 'qsd')
-    gate_v = unitary(v.T, 'qsd')
+        circuit.cx(reg_b[k], reg_a[k])
+
+    gate_u = decompose_unitary(svd_u, 'qsd')
+    gate_v = decompose_unitary(svd_v.T, 'qsd')
 
-    circuit.compose(gate_u, hsb, inplace=True) # apply gate U to the first register
-    circuit.compose(gate_v, hsa, inplace=True) # apply gate V to the second register
+    circuit.compose(gate_u, reg_b, inplace=True) # apply gate U to the first register
+    circuit.compose(gate_v, reg_a, inplace=True) # apply gate V to the second register
 
     return circuit
diff --git a/qclib/state_preparation/ventura.py b/qclib/state_preparation/ventura.py
index d864aac0..21d56ad7 100644
--- a/qclib/state_preparation/ventura.py
+++ b/qclib/state_preparation/ventura.py
@@ -12,11 +12,16 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+"""
+https://arxiv.org/abs/quant-ph/9807054
+"""
+
 import numpy as np
 from qiskit import QuantumCircuit, QuantumRegister
 
+# pylint: disable=maybe-no-member
 
-def initialize(state, n_qubits, N):
+def initialize(state, n_qubits, n_output_values):
     """ State preparation using Ventura and Martinez algorithm quant-ph/9807054
     Algorithm that requires a polynomial number of elementary operations for
     initializing a quantum system to represent only the m known points of a
@@ -46,8 +51,8 @@ def initialize(state, n_qubits, N):
     n_qubits: int
         Number of bits to represent the binary values.
 
-    N: int
-        Number of possible output values (Ex.: N=2 for a binary function).
+    n_output_values: int
+        Number of possible output values N (Ex.: n_output_values=2 for a binary function).
 
     Returns
     -------
@@ -55,82 +60,80 @@ def initialize(state, n_qubits, N):
         QuantumCircuit to initialize the state.
     """
 
-    x = QuantumRegister(n_qubits, 'x')
-    g = QuantumRegister(n_qubits - 1, 'g')
-    classical_reg = QuantumRegister(2, 'c')
-
-    circuit = QuantumCircuit(x, g, classical_reg)
+    reg_x = QuantumRegister(n_qubits, 'x')
+    reg_g = QuantumRegister(n_qubits - 1, 'g')
+    reg_c = QuantumRegister(2, 'c')
 
-    x = x[::-1] # qiskit reverse (qiskit little-endian)
+    circuit = QuantumCircuit(reg_x, reg_g, reg_c)
 
-    z0 = [int(k) for k in '{:0{}b}'.format(0, n_qubits)]
-    for p, (i, s) in list(enumerate(state.items()))[::-1]:
-        z = [int(k) for k in '{:0{}b}'.format(i, n_qubits)]
+    reg_x = reg_x[::-1] # qiskit reverse (qiskit little-endian)
 
-        circuit.x(classical_reg[1])
-        for j,_ in enumerate(z):
-            if z0[j] != z[j]:
-                circuit.cx(classical_reg[1], x[j])
+    bits_z0 = [int(k) for k in '{:0{}b}'.format(0, n_qubits)]
+    for idx_p, (input_z, output_s) in list(enumerate(state.items()))[::-1]:
+        bits_z = [int(k) for k in '{:0{}b}'.format(input_z, n_qubits)]
 
-        z0 = z
-        circuit.cx(classical_reg[1], classical_reg[0])
-        circuit.x(classical_reg[1])
+        circuit.x(reg_c[1])
+        for j,_ in enumerate(bits_z):
+            if bits_z0[j] != bits_z[j]:
+                circuit.cx(reg_c[1], reg_x[j])
 
-        # This sign is here for the smaller values of "s"
-        # to be represented by negative amplitudes and
-        # the larger ones by positive amplitudes.
-        theta = -2 * np.arccos( np.sqrt(p / (p+1)) )
+        bits_z0 = bits_z
+        circuit.cx(reg_c[1], reg_c[0])
+        circuit.x(reg_c[1])
 
-        # In the paper this negative sign is missing.
-        # Without it the matrix S is not unitary.
-        lamb = -s*2*np.pi/N
+        theta = -2*np.arccos(np.sqrt(idx_p/(idx_p+1))) # This sign is here for the smaller values
+                                                       # of "s" to be represented by negative
+                                                       # amplitudes and the larger ones by positive
+                                                       # amplitudes.
+        lamb = -output_s*2*np.pi/n_output_values # In the paper this negative sign is missing.
+                                                 # Without it the matrix S is not unitary.
         phi = -lamb
-        circuit.cu(theta, phi, lamb, 0, classical_reg[0], classical_reg[1])
+        circuit.cu(theta, phi, lamb, 0, reg_c[0], reg_c[1])
 
-        if z[0] == 0:
-            circuit.x(x[0])
-        if z[1] == 0:
-            circuit.x(x[1])
+        if bits_z[0] == 0:
+            circuit.x(reg_x[0])
+        if bits_z[1] == 0:
+            circuit.x(reg_x[1])
 
-        circuit.ccx(x[0], x[1], g[0])
+        circuit.ccx(reg_x[0], reg_x[1], reg_g[0])
 
-        if z[0] == 0:
-            circuit.x(x[0])
-        if z[1] == 0:
-            circuit.x(x[1])
+        if bits_z[0] == 0:
+            circuit.x(reg_x[0])
+        if bits_z[1] == 0:
+            circuit.x(reg_x[1])
 
         for k in range(2, n_qubits):
-            if z[k] == 0:
-                circuit.x(x[k])
+            if bits_z[k] == 0:
+                circuit.x(reg_x[k])
 
-            circuit.ccx(x[k], g[k-2], g[k-1])
+            circuit.ccx(reg_x[k], reg_g[k-2], reg_g[k-1])
 
-            if z[k] == 0:
-                circuit.x(x[k])
+            if bits_z[k] == 0:
+                circuit.x(reg_x[k])
 
-        circuit.cx(g[n_qubits - 2], classical_reg[0])
+        circuit.cx(reg_g[n_qubits - 2], reg_c[0])
 
         for k in range(n_qubits - 1, 1, -1):
-            if z[k] == 0:
-                circuit.x(x[k])
+            if bits_z[k] == 0:
+                circuit.x(reg_x[k])
 
-            circuit.ccx(x[k], g[k-2], g[k-1])
+            circuit.ccx(reg_x[k], reg_g[k-2], reg_g[k-1])
 
-            if z[k] == 0:
-                circuit.x(x[k])
+            if bits_z[k] == 0:
+                circuit.x(reg_x[k])
 
-        if z[0] == 0:
-            circuit.x(x[0])
-        if z[1] == 0:
-            circuit.x(x[1])
+        if bits_z[0] == 0:
+            circuit.x(reg_x[0])
+        if bits_z[1] == 0:
+            circuit.x(reg_x[1])
 
-        circuit.ccx(x[0], x[1], g[0])
+        circuit.ccx(reg_x[0], reg_x[1], reg_g[0])
 
-        if z[0] == 0:
-            circuit.x(x[0])
-        if z[1] == 0:
-            circuit.x(x[1])
+        if bits_z[0] == 0:
+            circuit.x(reg_x[0])
+        if bits_z[1] == 0:
+            circuit.x(reg_x[1])
 
-    circuit.x(classical_reg[1])
+    circuit.x(reg_c[1])
 
     return circuit
diff --git a/requirements.txt b/requirements.txt
index 7c6475fa..d559e824 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,3 +1,4 @@
 qiskit
 graphviz
+deprecation
 pytest-pylint


Test patch:
diff --git a/test/test_isometry.py b/test/test_isometry.py
index d9d9ab21..8a9aa40e 100644
--- a/test/test_isometry.py
+++ b/test/test_isometry.py
@@ -12,7 +12,9 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-""" Test isometry """
+"""
+Tests for the isometry.py module.
+"""
 
 from unittest import TestCase
 import numpy as np
@@ -23,137 +25,106 @@
 from qclib.unitary import unitary
 from qclib.state_preparation.mottonen import initialize
 
-class TestInitialize(TestCase):
-    """ Testing isometry """
-    def test_state_preparation_knill(self):
-        """ Testing state preparation with isometry Knill scheme """
-        vector = np.random.rand(32) + np.random.rand(32) * 1j
-        vector = vector / np.linalg.norm(vector)
+# pylint: disable=missing-function-docstring
+# pylint: disable=missing-class-docstring
 
-        circuit = decompose(vector, scheme='knill')
+class TestIsometry(TestCase):
 
-        state = get_state(circuit)
-
-        self.assertTrue(np.allclose(vector, state))
-
-    def test_unitary_knill(self):
-        """ Testing isometry Knill scheme with unitary matrix"""
-        unitary_matrix = unitary_group.rvs(32)
+    # scheme='knill' isometry tests
 
-        circuit = decompose(unitary_matrix, scheme='knill')
-
-        state = get_state(circuit)
+    def test_state_preparation_real_knill(self):
+        self._test_state_preparation_real('knill')
 
-        self.assertTrue(np.allclose(unitary_matrix[:, 0], state))
+    def test_state_preparation_complex_knill(self):
+        self._test_state_preparation_complex('knill')
 
     def test_isometry_knill(self):
-        """ Testing isometry Knill scheme with isometry matrix"""
-        isometry_matrix = unitary_group.rvs(32)[:, :4]
+        self._test_isometry('knill')
 
-        gate = decompose(isometry_matrix, scheme='knill')
-
-        state = get_state(gate)
-
-        self.assertTrue(np.allclose(isometry_matrix[:, 0], state))
-
-        circuit = QuantumCircuit(5)
-        circuit.x(0)
-        circuit.append(gate, circuit.qubits)
-        state = get_state(circuit)
-        self.assertTrue(np.allclose(isometry_matrix[:, 1], state))
+    def test_fixed_isometry_knill(self):
+        self._test_fixed_isometry('knill')
 
-        circuit = QuantumCircuit(5)
-        circuit.x(1)
-        circuit.append(gate, circuit.qubits)
-        state = get_state(circuit)
-        self.assertTrue(np.allclose(isometry_matrix[:, 2], state))
+    # scheme='ccd' (column-by-column decomposition) isometry tests
 
-        circuit = QuantumCircuit(5)
-        circuit.x([0,1])
-        circuit.append(gate, circuit.qubits)
-        state = get_state(circuit)
-        self.assertTrue(np.allclose(isometry_matrix[:, 3], state))
+    def test_state_preparation_real_ccd(self):
+        self._test_state_preparation_real('ccd')
 
-    def test_compare_mottonen_isometry_knill(self):
-        """ Compare Mottonen and Isometry state preparation"""
+    def test_state_preparation_complex_ccd(self):
+        self._test_state_preparation_complex('ccd')
 
-        vector = np.random.rand(32) + np.random.rand(32) * 1j
-        vector = vector / np.linalg.norm(vector)
+    def test_isometry_ccd(self):
+        self._test_isometry('ccd')
 
-        circuit1 = decompose(vector, scheme='knill')
-        circuit2 = initialize(vector)
+    def test_fixed_isometry_ccd(self):
+        self._test_fixed_isometry('ccd')
 
-        state1 = get_state(circuit1)
-        state2 = get_state(circuit2)
+    # general functions for testing the schemes (knill and ccd).
 
-        self.assertTrue(np.allclose(state1, state2))
+    def _test_state_preparation_real(self, scheme):
+        state_vector = np.random.rand(32) + np.random.rand(32) * 1j
+        state_vector = state_vector / np.linalg.norm(state_vector)
 
-    def test_compare_unitary_qsd_isometry_knill(self):
-        """ Compare unitary and isometry"""
-        unitary_matrix = unitary_group.rvs(32)
+        circuit = decompose(state_vector, scheme=scheme)
 
-        circuit1 = decompose(unitary_matrix, scheme='knill')
-        circuit2 = unitary(unitary_matrix, decomposition='qsd')
+        state = get_state(circuit)
 
-        state1 = get_state(circuit1)
-        state2 = get_state(circuit2)
+        self.assertTrue(np.allclose(state_vector, state))
 
-        self.assertTrue(np.allclose(state1, state2))
+    def _test_state_preparation_complex(self, scheme):
+        state_vector = np.random.rand(32) + np.random.rand(32) * 1j
+        state_vector = state_vector / np.linalg.norm(state_vector)
 
-    def test_compare_unitary_csd_isometry_knill(self):
-        """ Compare unitary and csd"""
-        unitary_matrix = unitary_group.rvs(32)
+        circuit = decompose(state_vector, scheme=scheme)
 
-        circuit1 = decompose(unitary_matrix, scheme='knill')
-        circuit2 = unitary(unitary_matrix, decomposition='csd')
+        state = get_state(circuit)
 
-        state1 = get_state(circuit1)
-        state2 = get_state(circuit2)
+        self.assertTrue(np.allclose(state_vector, state))
 
-        self.assertTrue(np.allclose(state1, state2))
+    def _test_isometry(self, scheme):
+        log_lines = 4
+        log_cols = 2
+        isometry = unitary_group.rvs(2**log_lines)[:,:2**log_cols]
 
-    def test_compare_unitary_qiskit_isometry_knill(self):
-        """ Compare qiskit isometry and Knill isometry"""
-        unitary_matrix = unitary_group.rvs(32)
+        gate = decompose(isometry, scheme=scheme)
 
-        circuit1 = decompose(unitary_matrix, scheme='knill')
-        circuit2 = QuantumCircuit(5)
+        for j in range(2**log_cols):
+            circuit = QuantumCircuit(log_lines)
 
-        # pylint: disable=maybe-no-member
-        circuit2.unitary(unitary_matrix, list(range(5)))
+            for i, bit in enumerate('{:0{}b}'.format(j, 1)[::-1]):
+                if bit == '1':
+                    circuit.x(i)
 
-        state1 = get_state(circuit1)
-        state2 = get_state(circuit2)
+            circuit.append(gate, circuit.qubits)
+            state = get_state(circuit)
 
-        self.assertTrue(np.allclose(state1, state2))
+            self.assertTrue(np.allclose(isometry[:, j], state))
 
-    def test_schmidt_isometry_knill(self):
-        """  Testing decompose """
-        unitary_gate = np.copy([[-0.5391073,  -0.12662419, -0.73739705, -0.38674956],
-                    [0.15705405,  0.20566939,  0.32663193, -0.9090356],
-                    [-0.77065035, -0.23739918,  0.59084039,  0.02544217],
-                    [-0.30132273,  0.9409081,  -0.02155946,  0.15307435]])
+    def _test_fixed_isometry(self, scheme):
+        isometry = np.copy([[-0.5391073,  -0.12662419, -0.73739705, -0.38674956],
+                            [ 0.15705405,  0.20566939,  0.32663193, -0.9090356 ],
+                            [-0.77065035, -0.23739918,  0.59084039,  0.02544217],
+                            [-0.30132273,  0.9409081,  -0.02155946,  0.15307435]])
 
-        gate = decompose(unitary_gate, scheme='knill')
+        gate = decompose(isometry, scheme=scheme)
 
         state = get_state(gate)
 
-        self.assertTrue(np.allclose(unitary_gate[:, 0], state))
+        self.assertTrue(np.allclose(isometry[:, 0], state))
 
         circuit = QuantumCircuit(2)
         circuit.x(0)
         circuit.append(gate, circuit.qubits)
         state = get_state(circuit)
-        self.assertTrue(np.allclose(unitary_gate[:, 1], state))
+        self.assertTrue(np.allclose(isometry[:, 1], state))
 
         circuit = QuantumCircuit(2)
         circuit.x(1)
         circuit.append(gate, circuit.qubits)
         state = get_state(circuit)
-        self.assertTrue(np.allclose(unitary_gate[:, 2], state))
+        self.assertTrue(np.allclose(isometry[:, 2], state))
 
         circuit = QuantumCircuit(2)
         circuit.x([0,1])
         circuit.append(gate, circuit.qubits)
         state = get_state(circuit)
-        self.assertTrue(np.allclose(unitary_gate[:, 3], state))
+        self.assertTrue(np.allclose(isometry[:, 3], state))
diff --git a/test/test_mottonen.py b/test/test_mottonen.py
index 28f7cb7e..8b44d5e8 100644
--- a/test/test_mottonen.py
+++ b/test/test_mottonen.py
@@ -12,62 +12,77 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-""" Test state preparation """
+"""
+Tests for the mottonen.py module.
+"""
 
 from unittest import TestCase
 import numpy as np
-from qiskit import execute, Aer
+from qiskit import ClassicalRegister, execute, Aer
 from qclib.state_preparation.mottonen import initialize
 from qclib.util import get_state
 
+# pylint: disable=missing-function-docstring
+# pylint: disable=missing-class-docstring
+
 backend = Aer.get_backend('qasm_simulator')
 SHOTS = 8192
 
-
 class TestInitialize(TestCase):
-    """ Testing qclib.state_preparation.mottonen """
     @staticmethod
-    def measurement(circuit):
-        """ get state preparation counts"""
+    def measurement(circuit, n_qubits, classical_reg):
+        circuit.measure(list(range(n_qubits)), classical_reg)
 
-        circuit.measure_all()
-        job = execute(circuit, backend, shots=SHOTS, optimization_level=0)
+        job = execute(circuit, backend, shots=SHOTS, optimization_level=3)
 
         counts = job.result().get_counts(circuit)
+        sum_values = sum(counts.values())
 
         counts2 = {}
-        for j in range(2**circuit.num_qubits):
-            pattern = '{:0{}b}'.format(j, circuit.num_qubits)
+        for i in range(2**n_qubits):
+            pattern = '{:0{}b}'.format(i, n_qubits)
             if pattern in counts:
                 counts2[pattern] = counts[pattern]
             else:
                 counts2[pattern] = 0.0
 
-        return [value/SHOTS for (key, value) in counts2.items()]
+        return [ value/sum_values for (key, value) in counts2.items() ]
 
     @staticmethod
     def mottonen_experiment(state):
-        """ Creates circuit with Mottonen state preparation algorithm"""
         circuit = initialize(state)
 
-        return TestInitialize.measurement(circuit)
+        n_qubits = int(np.log2(len(state)))
+        classical_reg = ClassicalRegister(n_qubits)
+        circuit.add_register(classical_reg)
+
+        return TestInitialize.measurement(circuit, n_qubits, classical_reg)
+
+    def test_mottonen_state_real(self):
+        state_vector = np.random.rand(32)
+        state_vector = state_vector / np.linalg.norm(state_vector)
+
+        circuit = initialize(state_vector)
+
+        state = get_state(circuit)
+
+        self.assertTrue(np.allclose(state_vector, state))
 
-    def test_mottonen_state(self):
-        """ Testiong Mottonen state preparation """
-        vector = np.random.rand(32) + np.random.rand(32) * 1j
-        vector = vector / np.linalg.norm(vector)
+    def test_mottonen_state_complex(self):
+        state_vector = np.random.rand(32) + np.random.rand(32) * 1j
+        state_vector = state_vector / np.linalg.norm(state_vector)
 
-        circuit = initialize(vector)
+        circuit = initialize(state_vector)
 
         state = get_state(circuit)
 
-        self.assertTrue(np.allclose(vector, state))
+        self.assertTrue(np.allclose(state_vector, state))
 
     def test_mottonen_measure(self):
-        """ Testiong Mottonen state preparation with measurements"""
-        vector = np.random.rand(32) + np.random.rand(32) * 1j
-        vector = vector / np.linalg.norm(vector)
+        state_vector = np.random.rand(32) + np.random.rand(32) * 1j
+        state_vector = state_vector / np.linalg.norm(state_vector)
 
-        state = TestInitialize.mottonen_experiment(vector)
+        state = TestInitialize.mottonen_experiment(state_vector)
 
-        self.assertTrue(np.allclose( np.power(np.abs(vector), 2), state, rtol=1e-01, atol=0.005))
+        self.assertTrue(np.allclose( np.power(np.abs(state_vector),2), state,
+                        rtol=1e-01, atol=0.005))
diff --git a/test/test_schmidt.py b/test/test_schmidt.py
index 8942ff2c..7e3c6afc 100644
--- a/test/test_schmidt.py
+++ b/test/test_schmidt.py
@@ -12,86 +12,64 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-""" Test state preparation """
+"""
+Tests for the schmidt.py module.
+"""
 
 from unittest import TestCase
 import numpy as np
 from qclib.state_preparation.schmidt import initialize
 from qclib.util import get_state
 
+# pylint: disable=missing-function-docstring
+# pylint: disable=missing-class-docstring
 
 class TestInitialize(TestCase):
-    """ Testing qclib.state_preparation.schmidt """
     @staticmethod
     def mae(state, ideal):
         """
-         Mean Absolute Error
+        Mean Absolute Error
         """
         return np.sum(np.abs(state-ideal))/len(ideal)
 
-    def test_initialize(self):
-        """ Test state vector initialization """
-        vector = np.random.rand(32) + np.random.rand(32) * 1j
-        vector = vector / np.linalg.norm(vector)
+    def _test_initialize_rank(self, rank, max_mae=10**-17):
+        state_vector = np.random.rand(32) + np.random.rand(32) * 1j
+        state_vector = state_vector / np.linalg.norm(state_vector)
 
-        circuit = initialize(vector)
+        circuit = initialize(state_vector, low_rank=rank)
 
         state = get_state(circuit)
 
-        self.assertTrue(np.allclose(vector, state))
+        self.assertTrue(TestInitialize.mae(state,state_vector) < max_mae)
 
-    def test_initialize_rank_5(self):
-        """ Test rank 5 state vector initialization """
-        vector = np.random.rand(32) + np.random.rand(32) * 1j
-        vector = vector / np.linalg.norm(vector)
+    def test_initialize_full_rank(self):
+        state_vector = np.random.rand(32) + np.random.rand(32) * 1j
+        state_vector = state_vector / np.linalg.norm(state_vector)
 
-        circuit = initialize(a, low_rank=5)
+        circuit = initialize(state_vector)
 
         state = get_state(circuit)
 
-        self.assertTrue(np.allclose(vector, state)) # same as unitary.
+        self.assertTrue(np.allclose(state_vector, state))
 
-    def test_initialize_rank_4(self):
-        """ Test rank 4 state vector initialization """
-        vector = np.random.rand(32) + np.random.rand(32) * 1j
-        vector = vector / np.linalg.norm(vector)
+    def test_initialize_rank_5(self):
+        state_vector = np.random.rand(32) + np.random.rand(32) * 1j
+        state_vector = state_vector / np.linalg.norm(state_vector)
 
-        circuit = initialize(a, low_rank=4)
+        circuit = initialize(state_vector, low_rank=5)
 
         state = get_state(circuit)
-        
-        self.assertTrue(TestInitialize.mae(state,a) < 10**-14)
-
-    def test_initialize_rank_3(self):
-        """ Test rank 3 state vector initialization """
-        vector = np.random.rand(32) + np.random.rand(32) * 1j
-        vector = vector / np.linalg.norm(vector)
 
-        circuit = initialize(a, low_rank=3)
+        self.assertTrue(np.allclose(state_vector, state)) # same as unitary.
 
-        state = get_state(circuit)
+    def test_initialize_rank_4(self):
+        self._test_initialize_rank(4, max_mae=10**-14)
 
-        self.assertTrue(TestInitialize.mae(state,vector) < 0.04)
+    def test_initialize_rank_3(self):
+        self._test_initialize_rank(3, max_mae=0.04)
 
     def test_initialize_rank_2(self):
-        """ Test rank 2 state vector initialization """
-        vector = np.random.rand(32) + np.random.rand(32) * 1j
-        vector = vector / np.linalg.norm(vector)
-
-        circuit = initialize(a, low_rank=2)
-
-        state = get_state(circuit)
-
-        self.assertTrue(TestInitialize.mae(state,vector) < 0.055)
+        self._test_initialize_rank(2, max_mae=0.06)
 
     def test_initialize_rank_1(self):
-        """ Test rank 1 state vector initialization """
-        vector = np.random.rand(32) + np.random.rand(32) * 1j
-        vector = vector / np.linalg.norm(vector)
-
-        circuit = initialize(a, low_rank=1)
-
-        state = get_state(circuit)
-        
-        self.assertTrue(TestInitialize.mae(state,a) < 0.0825)
-
+        self._test_initialize_rank(1, max_mae=0.0825)
diff --git a/test/test_ventura.py b/test/test_ventura.py
index 4039c0d0..3439f87d 100644
--- a/test/test_ventura.py
+++ b/test/test_ventura.py
@@ -32,7 +32,7 @@ def test_binary_function(self):
         max_denominator = max(function_io.values())-1
         n_qubits = 4
 
-        circuit = initialize(function_io, n_qubits=n_qubits, N=max_denominator)
+        circuit = initialize(function_io, n_qubits=n_qubits, n_output_values=max_denominator)
         state = get_state(circuit)
 
         desired_state = []