filterdesign/filterdesign.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimization problem over nonnegative trigonometric polynomials. #\n",
    "\n",
    "We consider the trigonometric polynomial\n",
    " \n",
    "$$H(\\omega) = x_0 + 2 \\sum_{k=1}^n [ \\Re(x_k) \\cos(\\omega k) + \\Im(x_k) \\sin(\\omega k) ],$$\n",
    "\n",
    "where $H(\\omega)$ is a real valued function paramtrized by the complex vector $x\\in {\\bf C}^{n+1}$, and where $\\Im(x_0)=0$.\n",
    "\n",
    "The example shows how to construct a *non-negative* polynomial $H(\\omega)\\geq 0, \\: \\forall \\omega$ that satisfies,\n",
    "  \n",
    "$$1 - \\delta \\leq  H(\\omega) \\leq 1 + \\delta, \\quad  \\forall \\omega \\in [0, \\omega_p]$$\n",
    "\n",
    "while minimizing $\\sup_{\\omega\\in [\\omega_s,\\pi]} H(\\omega)$ over the variables $x$.\n",
    "\n",
    "In the signal processing literature such a trigonometric polynomial is known as (the squared amplitude response of) a Chebyshev lowpass filter. \n",
    "\n",
    "A squared amplitude response $H(\\omega)$ is always symmetric around $0$, so $\\Im(x_k)=0$, and we consider only\n",
    "\n",
    "$$H(\\omega) = x_0 + 2 \\sum_{k=1}^n x_k \\cos(\\omega k) $$\n",
    "\n",
    "over the real vector $x\\in {\\bf R}^{n+1}$. However, the concepts in the example are readily applied to the case with $x\\in {\\bf C}^{n+1}$.\n",
    "\n",
    "References:\n",
    "  1. \"Squared Functional Systems and Optimization Problems\",  Y. Nesterov, 2004.\n",
    "  \n",
    "  2. \"Convex Optimization of Non-negative Polynomials: Structured Algorithms and Applications\", Ph.D thesis, Y. Hachez, 2003.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import mosek\n",
    "from   mosek.fusion import *\n",
    "import mosek.fusion.pythonic\n",
    "from   math import cos, pi\n",
    "import numpy as np\n",
    "import sys"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nonnegativity everywhere ###\n",
    "\n",
    "Nesterov proved in [1] that $H(\\omega) \\geq 0, \\: \\forall \\omega$ if and only if \n",
    "$$x_i = \\langle T_i^{n+1}, X \\rangle, \\quad X \\in {\\mathcal H}^{n+1}_+,$$\n",
    "where ${\\mathcal H}_+$ is the cone of Hermitian semidefinite matrices and $T_i$ is a Toeplitz matrix\n",
    "$$[T_i]_{kl} = \\left \\{ \\begin{array}{ll}1, & k-l=i\\\\0 & \\text{otherwise}.\\end{array} \\right .$$\n",
    "For example, for $n=2$ we have\n",
    "$$\n",
    "   T_0 = \\left[\\begin{array}{ccc}\n",
    "   1 & 0 & 0\\\\0 & 1 & 0\\\\0 & 0 & 1\n",
    "   \\end{array}\n",
    "   \\right], \\quad\n",
    "   T_1 = \\left[\\begin{array}{ccc}\n",
    "   0 & 0 & 0\\\\1 & 0 & 0\\\\0 & 1 & 0\n",
    "   \\end{array}\n",
    "   \\right],\n",
    "   \\quad\n",
    "   T_2 = \\left[\\begin{array}{ccc}\n",
    "   0 & 0 & 0\\\\0 & 0 & 0\\\\1 & 0 & 0\n",
    "   \\end{array}\n",
    "   \\right].\n",
    "$$\n",
    "In our case we have $\\Im(x_i)=0$, i.e., $X\\in {\\mathcal S}^{n+1}_+$ is a real symmetric semidefinite matrix.\n",
    "\n",
    "We define the *cone on nonnegative trigonometric polynomials* as\n",
    "$$\n",
    "   K^n_{0,\\pi} = \\{ x\\in \\mathbf{R} \\times \\mathbf{C}^n \\mid x_i = \\langle X, T_i\\rangle, \\: i=0,\\dots,n, \\: X\\in\\mathcal{H}_+^{n+1} \\}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def T_dot_X(n, i, X, a=1.0):\n",
    "    if i>=n or i<=-n: return Expr.constTerm(0.0)\n",
    "    else: return Expr.dot(Matrix.diag(n, a, -i), X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trigpoly_0_pi(M, x):\n",
    "    '''x[i] == <T(n+1,i),X>'''\n",
    "    n = int(x.getSize()-1)\n",
    "    X = M.variable(\"X\", Domain.inPSDCone(n+1))\n",
    "    \n",
    "    for i in range(n+1):\n",
    "        M.constraint(T_dot_X(n+1,i,X) == x[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we have dropped the imaginary part of $X$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nonnegativity on $[0,a]$ ###\n",
    "\n",
    "Similarly, $H(\\omega)$ is nonnegative on $[0,a]$ if and only if\n",
    "\n",
    "$$x_i =\n",
    "\\langle X_1, T_i^{n+1} \\rangle + \n",
    "\\langle X_2, T_{i+1}^n \\rangle +\n",
    "\\langle X_2, T_{i-1}^n \\rangle -\n",
    "2 \\cos(a)\\langle X_2, T_{i}^n \\rangle, \\quad \n",
    "X_1 \\in {\\mathcal H}^{n+1}_+, \\:\n",
    "X_2 \\in {\\mathcal H}^n_+,\n",
    "$$\n",
    "or equivalently\n",
    "$$\n",
    "  K^n_{0,a} = \\{ x\\in \\mathbf{R}\\times \\mathbf{C}^n \\mid\n",
    "   x_i = \\langle X_1, T_i^{n+1} \\rangle +\n",
    "   \\langle X_2 , T_{i+1}^n \\rangle +\n",
    "   \\langle X_2 , T_{i-1}^n \\rangle -\n",
    "   2\\cos(a)\\langle X_2 , T_i^n\\rangle, \\: X_1\\in \\mathcal{H}_+^{n+1}, X_2\\in\\mathcal{H}_+^n \\}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trigpoly_0_a(M, x, a):\n",
    "    '''x[i] == <T(n+1,i),X1> + <T(n,i+1),X2> + <T(n,i-1),X2> - 2*cos(a)*<T(n,i),X2>'''\n",
    "    n = int(x.getSize()-1)\n",
    "    X1 = M.variable(Domain.inPSDCone(n+1))\n",
    "    X2 = M.variable(Domain.inPSDCone(n))\n",
    "\n",
    "    for i in range(n+1):\n",
    "        M.constraint(T_dot_X(n+1,i,X1) + T_dot_X(n,i+1,X2) + \\\n",
    "                     T_dot_X(n,i-1,X2) + T_dot_X(n,i,X2,-2*cos(a)) == x[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we have dropped the imaginary part of $X_1$ and $X_2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Nonnegativity on $[a,\\pi]$ ###\n",
    "\n",
    "Finally, $H(\\omega)$ is nonnegative on $[a,\\pi]$ if and only if\n",
    "\n",
    "$$x_i = \n",
    "\\langle X_1, T_i^{n+1} \\rangle -\n",
    "\\langle X_2, T_{i+1}^n \\rangle -\n",
    "\\langle X_2, T_{i-1}^n \\rangle +\n",
    "2 \\cos(a)\\langle X_2, T_{i}^n \\rangle, \\quad \n",
    "X_1 \\in {\\mathcal S}^{n+1}_+, \\:\n",
    "X_2 \\in {\\mathcal S}^n_+,\n",
    "$$\n",
    "or equivalently\n",
    "$$\n",
    "  K^n_{a,\\pi} = \\{ x\\in \\mathbf{R}\\times \\mathbf{C}^n \\mid\n",
    "   x_i = \\langle X_1, T_i^{n+1} \\rangle -\n",
    "   \\langle X_2 , T_{i+1}^n \\rangle -\n",
    "   \\langle X_2 , T_{i-1}^n \\rangle +\n",
    "   2\\cos(a)\\langle X_2 , T_i^n\\rangle, \\: X_1\\in \\mathcal{H}_+^{n+1}, X_2\\in\\mathcal{H}_+^n \\}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trigpoly_a_pi(M, x, a):\n",
    "    '''x[i] == <T(n+1,i),X1> - <T(n,i+1),X2> - <T(n,i-1),X2> + 2*cos(a)*<T(n,i),X2>'''\n",
    "    n = int(x.getSize()-1)\n",
    "    X1 = M.variable(Domain.inPSDCone(n+1))\n",
    "    X2 = M.variable(Domain.inPSDCone(n))\n",
    "\n",
    "    for i in range(n+1):\n",
    "        M.constraint(T_dot_X(n+1,i,X1) + T_dot_X(n,i+1,X2,-1.0) + \\\n",
    "                     T_dot_X(n,i-1,X2,-1.0) + T_dot_X(n,i,X2,2*cos(a)) == x[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we have dropped the imaginary part of $X_1$ and $X_2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An epigraph formulation ##\n",
    "The *epigraph* $H(\\omega) \\leq t$ can now be characterized simply as \n",
    "$-u \\in K^n_{[a,b]}, \\: u=(x_0-t, \\, x_{1:n}).$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def epigraph(M, x, t, a, b):\n",
    "    '''Models 0 <= H(w) <= t, for all w in [a, b], where\n",
    "         H(w) = x0 + 2*x1*cos(w) + 2*x2*cos(2*w) + ... + 2*xn*cos(n*w)'''\n",
    "    n  = int(x.getSize()-1)    \n",
    "    u = M.variable(n+1, Domain.unbounded())\n",
    "    M.constraint(t == x[0] + u[0]) \n",
    "    M.constraint(x[1:] + u[1:] == 0)\n",
    "    \n",
    "    if a==0.0 and b==pi:\n",
    "        trigpoly_0_pi(M, u)\n",
    "    elif a==0.0 and b<pi:\n",
    "        trigpoly_0_a(M, u, b)\n",
    "    elif a<pi and b==pi:\n",
    "        trigpoly_a_pi(M, u, a)\n",
    "    else:\n",
    "        raise ValueError(\"invalid interval.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A hypograph formulation ##\n",
    "Similarly, the *hypograph* $H(\\omega) \\geq t$ can be characterized as \n",
    "$u \\in K^n_{[a,b]}, \\: u=(x_0-t, \\, x_{1:n}).$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hypograph(M, x, t, a, b):\n",
    "    '''Models 0 <= t <= H(w), for all w in [a, b], where\n",
    "         H(w) = x0 + 2*x1*cos(w) + 2*x2*cos(2*w) + ... + 2*xn*cos(n*w)'''\n",
    "\n",
    "    n  = int(x.getSize()-1)    \n",
    "    u0 = M.variable(1, Domain.unbounded())    \n",
    "    M.constraint(t == x[0] - u0)\n",
    "    u = Var.vstack( u0, x[1:] )\n",
    "\n",
    "    if a==0.0 and b==pi:\n",
    "        trigpoly_0_pi(M, u)\n",
    "    elif a==0.0 and b<pi:\n",
    "        trigpoly_0_a(M, u,  b)\n",
    "    elif a<pi and b==pi:\n",
    "        trigpoly_a_pi(M, u, a)\n",
    "    else:\n",
    "        raise ValueError(\"invalid interval.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Putting it together ##"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 10\n",
    "M = Model(\"trigpoly\")\n",
    "x = M.variable(\"x\", n+1, Domain.unbounded())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Global nonnegativity ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# H(w) >= 0\n",
    "trigpoly_0_pi(M, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Passband specifications ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "wp = pi/4\n",
    "delta = 0.05\n",
    "\n",
    "# H(w) <= (1+delta), w in [0, wp]    \n",
    "epigraph(M, x, 1.0+delta, 0.0, wp)\n",
    "\n",
    "# (1-delta) <= H(w), w in [0, wp]\n",
    "hypograph(M, x, 1.0-delta, 0.0, wp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stopband specifications ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "ws = wp + pi/8\n",
    "\n",
    "# H(w) < t, w in [ws, pi]\n",
    "t = M.variable(\"t\", 1, Domain.greaterThan(0.0))\n",
    "epigraph(M, x, t, ws, pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting the objective ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "M.objective(ObjectiveSense.Minimize, t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "M.setLogHandler(sys.stdout)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Problem\n",
      "  Name                   : trigpoly        \n",
      "  Objective sense        : min             \n",
      "  Type                   : CONIC (conic optimization problem)\n",
      "  Constraints            : 67              \n",
      "  Cones                  : 0               \n",
      "  Scalar variables       : 35              \n",
      "  Matrix variables       : 7               \n",
      "  Integer variables      : 0               \n",
      "\n",
      "Optimizer started.\n",
      "Presolve started.\n",
      "Linear dependency checker started.\n",
      "Linear dependency checker terminated.\n",
      "Eliminator started.\n",
      "Freed constraints in eliminator : 23\n",
      "Eliminator terminated.\n",
      "Eliminator started.\n",
      "Freed constraints in eliminator : 0\n",
      "Eliminator terminated.\n",
      "Eliminator - tries                  : 2                 time                   : 0.00            \n",
      "Lin. dep.  - tries                  : 1                 time                   : 0.00            \n",
      "Lin. dep.  - number                 : 0               \n",
      "Presolve terminated. Time: 0.01    \n",
      "Problem\n",
      "  Name                   : trigpoly        \n",
      "  Objective sense        : min             \n",
      "  Type                   : CONIC (conic optimization problem)\n",
      "  Constraints            : 67              \n",
      "  Cones                  : 0               \n",
      "  Scalar variables       : 35              \n",
      "  Matrix variables       : 7               \n",
      "  Integer variables      : 0               \n",
      "\n",
      "Optimizer  - threads                : 20              \n",
      "Optimizer  - solved problem         : the primal      \n",
      "Optimizer  - Constraints            : 44\n",
      "Optimizer  - Cones                  : 1\n",
      "Optimizer  - Scalar variables       : 13                conic                  : 12              \n",
      "Optimizer  - Semi-definite variables: 7                 scalarized             : 429             \n",
      "Factor     - setup time             : 0.00              dense det. time        : 0.00            \n",
      "Factor     - ML order time          : 0.00              GP order time          : 0.00            \n",
      "Factor     - nonzeros before factor : 990               after factor           : 990             \n",
      "Factor     - dense dim.             : 0                 flops                  : 1.18e+05        \n",
      "ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  \n",
      "0   1.9e+01  1.0e+00  1.0e+00  0.00e+00   0.000000000e+00   0.000000000e+00   1.0e+00  0.03  \n",
      "1   6.0e+00  3.2e-01  3.2e-01  1.56e+00   -6.168509433e-03  -1.704940542e-02  3.2e-01  0.04  \n",
      "2   1.5e+00  8.0e-02  8.0e-02  1.59e+00   -3.629686898e-02  -5.729238420e-02  8.0e-02  0.05  \n",
      "3   3.1e-01  1.7e-02  1.7e-02  1.64e+00   -4.704273187e-03  -7.156318767e-03  1.7e-02  0.05  \n",
      "4   1.0e-01  5.6e-03  5.6e-03  1.91e-01   1.247879338e-03   3.220303293e-03   5.6e-03  0.05  \n",
      "5   5.8e-02  3.1e-03  3.1e-03  4.28e-01   4.186115053e-02   3.996134158e-02   3.1e-03  0.05  \n",
      "6   1.1e-02  5.9e-04  5.9e-04  4.33e-01   6.583025991e-02   6.538915159e-02   5.9e-04  0.05  \n",
      "7   2.7e-03  1.4e-04  1.4e-04  9.75e-01   7.102893174e-02   7.092631673e-02   1.4e-04  0.05  \n",
      "8   4.6e-04  2.4e-05  2.4e-05  1.01e+00   7.246509927e-02   7.245345168e-02   2.4e-05  0.06  \n",
      "9   3.3e-05  1.8e-06  1.8e-06  1.00e+00   7.271267688e-02   7.271184530e-02   1.8e-06  0.06  \n",
      "10  3.6e-06  1.9e-07  1.9e-07  1.00e+00   7.272951727e-02   7.272942593e-02   1.9e-07  0.06  \n",
      "11  5.8e-07  3.1e-08  3.1e-08  1.00e+00   7.273126748e-02   7.273125303e-02   3.1e-08  0.06  \n",
      "12  7.0e-08  3.9e-09  3.7e-09  1.00e+00   7.273156244e-02   7.273156074e-02   3.7e-09  0.06  \n",
      "13  6.3e-09  1.9e-09  3.4e-10  1.00e+00   7.273159968e-02   7.273159953e-02   3.4e-10  0.07  \n",
      "Optimizer terminated. Time: 0.07    \n",
      "\n",
      "\n",
      "Interior-point solution summary\n",
      "  Problem status  : PRIMAL_AND_DUAL_FEASIBLE\n",
      "  Solution status : OPTIMAL\n",
      "  Primal.  obj: 7.2731599679e-02    nrm: 1e+00    Viol.  con: 6e-09    var: 0e+00    barvar: 0e+00  \n",
      "  Dual.    obj: 7.2731599528e-02    nrm: 2e+00    Viol.  con: 0e+00    var: 3e-10    barvar: 2e-09  \n"
     ]
    }
   ],
   "source": [
    "M.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.33580937,  0.25254317,  0.1388759 ,  0.02044169, -0.04797007,\n",
       "       -0.05156956, -0.01627207,  0.01837092,  0.02948826,  0.02271378,\n",
       "       -0.00952669])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x.level()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.0727316])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t.level()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting the amplitude response ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XtUVXXeP/D3h4viBVETBW+hqKnk5cBJpAywNPGSZk1Ol3Gq1TxUs2wyny5TM9rU5DyreSpb489Kp8dsctSc+lWOqWC/sryiIHhPRfOClyDwgooinM/vD5AHEWQD+5y9zznv11qsdS5f9n7vTr7P5rvPPltUFURE5FsCrA5ARETmY7kTEfkgljsRkQ9iuRMR+SCWOxGRD2K5ExH5IJY7EZEPYrkTEfkgljsRkQ8KsmrFHTp00KioKKtWT0TklbKysn5W1fD6xllW7lFRUcjMzLRq9UREXklEDhsZx2kZIiIfxHInIvJBLHciIh/Ecici8kEsdyIiH8RyJyLyQSx3IiIfZNnn3Btr/fr1+Oabb3DjjTciKioKvXv3RmRkpNWx3Or48eNYsWIFAgMDER0djejoaHTu3BkiYnU0IrIpryv3DRs2YMaMGVc9dtttt2Hy5MmYNGkS2rVrZ1Eyc5WWluL999/H4sWLsWnTpmueT0xMxEcffQSe5UtEtfG6aZnnn38eJSUl2Lt3L9LS0jBz5kwUFRXhySefRGRkJF5++WUUFxdbHbNJcnJyMGTIEDzzzDO4dOkSXn/9dezcuRP79+9HWloa3njjDWRnZ2PgwIH48MMPwYucE9E1VNWSn7i4ODWLy+XSrKwsffjhhxWARkRE6Pz587W8vNy0dXhCWVmZvvLKKxoUFKQRERH6xRdf1Dn2xx9/1KSkJAWgTz75pAdTEpGVAGSqgY71iXKvbuPGjRofH68AdMyYMVpQUOCW9ZjtwoULeu+99yoA/dWvfqWFhYX1/k55eblOmzZNAeiHH37o/pBEZDm/LXfVitKbPXu2NmvWTDt37qxr1qxx27rMUFhYqMOGDVMR0XfeeadBv1tWVqbDhw/XFi1a6I4dO9yUkIjswmi5e92cuxEBAQGYMmUKNm3ahFatWuGOO+7AW2+9Zcu56WPHjmHYsGHYvHkzlixZgmeeeaZBvx8YGIhFixahTZs2uP/++3Hu3Dk3JSUib1JvuYvIfBHJF5GddTwvIvI3EckVke0iEmt+zMZxOBzIysrCxIkT8dxzzyE1NRWlpaVWx6py7NgxDB8+HHl5eUhLS8OkSZMatZyIiAgsXrwY+/bta/CbAxH5JiN77gsApFzn+dEAelf+pAJ4r+mxzBMaGoqlS5fiD3/4Az744AOkpKSgqKjI6ljIy8tDcnIyTp48ibS0NCQnJzdpecOHD8ezzz6LBQsW4IcffjAnJBF5rXrLXVW/B3C9NpwA4B+V00GbALQVEVudVRQQEIDXX38d//jHP7B+/XoMGzYMhw8b+r57tzh69CiSk5Px008/IS0tDQkJCaYs98UXX0RISAhef/11U5ZHRN7LjDn3LgCOVrufV/mY7UyePBnp6ek4ceIEhg4diuzsbI9nOHLkCJKTk1FQUID09HTTih0AwsPDMWXKFCxevJh770R+zoxyr+0c+FqPXIpIqohkikhmQUGBCatuuKSkJKxbtw7BwcFITEzEqlWrPLbuw4cPIzk5GT///DPS09MxdOhQ09fx3HPPce+diEwp9zwA3ard7wrgeG0DVXWeqjpV1RkeXu/1Xd0mJiYGmzZtQq9evTB27FjMmTPH7es8cOAAkpKScOrUKXz99deIj493y3q4905EgDnlvgzArys/NTMUwBlVPWHCct2qc+fOWLt2LcaOHYspU6bgd7/7HcrKytyyri1btuDWW29FcXExvv76a9xyyy1uWc8VV/beZ86c6db1EJF9Gfko5GIAGwHcJCJ5IvK4iDwpIk9WDlkB4CCAXAB/B/Bbt6U1WevWrfH5559j2rRpmD17NkaNGoWTJ0+auo6vvvoKycnJaNmyJdavX4+4uDhTl1+b8PBw/OY3v8HSpUtt8ckgIrKAkTOd3PHjzjNUG2P+/PnaokUL7dixo6anpzd5eeXl5frXv/5VAwMDNTY2Vk+cOGFCSuOys7MVgM6ZM8ej6yUi94I/n6HaGI899hi2bNmCDh06YNSoUZg6dSpOnz7dqGUdOXIEd955J1544QVMmDABa9asQUREhMmJr2/w4MEYOHAgPvroI4+ul4jsgeVeTUxMDLZs2YInnngCf/vb39C7d2+8//77hufiz507h1mzZmHgwIHIzMzE/Pnz8emnnyI0NNTNyWv3yCOPYPPmzTywSuSHWO41tGzZEu+99x6ysrLQv39/PPXUU+jZsydeeOEFZGdnX/P9NGVlZdi5cyemT5+O7t27Y9q0aYiLi0NOTg4ee+wxS6+W9PDDDyMwMJB770R+SGqWlac4nU7NzMy0ZN1GqSq+/PJL/P3vf0d6ejrKysrQrl07dOrUCR07dkRJSQl27NiBixcvQkRwzz334MUXX3TbxxwbY9y4ccjJycHhw4cRGBhodRwiaiIRyVJVZ73jWO7GFBYW4rPPPsO2bduQn5+P/Px8BAUFYfDgwXA4HLjtttvQo0cPq2Ne41//+hcmTZqE9PR0jBw50uo4RNRELHcCAFy8eBGRkZEYO3YsFi5caHUcImoio+XOOXcfFxISgkmTJuGLL77AxYsXrY5DRB7CcvcD99xzD86fP4/vvvvO6ihE5CEsdz8wfPhwtGzZEsuXL7c6ChF5CMvdD4SEhGDEiBFYvny5LS81SETmY7n7iXHjxuHQoUPYvXu31VGIyANY7n5izJgxAMCpGSI/wXL3E126dEFsbCzLnchPsNz9yLhx47BhwwYUFhZaHYWI3Izl7kfGjRsHl8vl0UsLEpE1WO5+JC4uDp06deLUDJEfYLn7kYCAAIwZMwarVq1CeXm51XGIyI1Y7n5m5MiROH36NHJycqyOQkRuxHL3M8nJyQCANWvWWJqDiNyL5e5nIiMjcdNNN+Hbb7+1OgoRuRHL3Q8lJydj7dq1hi8fSETeh+Xuh5KTk3H27FnOuxP5MJa7H0pKSgLAeXciX8Zy90OcdyfyfSx3PzV8+HDOuxP5MJa7n0pOTkZxcTGys7OtjkJEbsBy91OcdyfybSx3PxUREYG+ffty3p3IR7Hc/Rjn3Yl8l6FyF5EUEdkrIrki8vtanu8uIt+KSLaIbBeRMeZHJbPdfvvtOHfuHHbs2GF1FCIyWb3lLiKBAOYAGA2gP4AHRaR/jWF/BLBUVR0AHgDwrtlByXxDhw4FAGzcuNHiJERkNiN77kMA5KrqQVUtBbAEwIQaYxRAm8rbYQCOmxeR3CUqKgqdOnXCpk2brI5CRCYzUu5dABytdj+v8rHq/gTgVyKSB2AFgKdrW5CIpIpIpohkFhQUNCIumUlEkJCQwD13Ih9kpNyllse0xv0HASxQ1a4AxgD4WESuWbaqzlNVp6o6w8PDG56WTDd06FDk5uaCb7ZEvsVIuecB6FbtfldcO+3yOIClAKCqGwGEAOhgRkByr4SEBABARkaGxUmIyExGyn0LgN4i0kNEmqHigOmyGmOOALgTAESkHyrKnbuCXsDpdCIoKIhTM0Q+pt5yV9UyAFMApAHYg4pPxewSkddEZHzlsP8E8B8isg3AYgCPqmrNqRuyoZYtW2LQoEEsdyIfE2RkkKquQMWB0uqPzah2ezeA28yNRp4ydOhQLFiwAGVlZQgKMvS/BBHZHM9QJSQkJOD8+fPYtWuX1VGIyCQsd6o6qMqpGSLfwXIn9OjRA+Hh4Sx3Ih/Ccqeqk5l4piqR72C5E4CKqZl9+/ahsLDQ6ihEZAKWOwEA4uPjAQCZmZkWJyEiM7DcCQAQGxsLgOVO5CtY7gQACAsLQ58+fVjuRD6C5U5VnE4ny53IR7DcqYrT6UReXh5OnjxpdRQiaiKWO1VxOp0AgKysLIuTEFFTsdypisPhgIhwaobIB7DcqUrr1q3Rr18/ljuRD2C501WuHFTlNzYTeTeWO13F6XTi5MmTOH6c1zgn8mYsd7rKlYOqnJoh8m4sd7rKoEGDEBgYyHIn8nIsd7pKy5YtERMTw3In8nIsd7oGD6oSeT+WO13D6XTi559/xuHDh62OQkSNxHKna8TFxQEAtm7danESImosljtdY8CAAQgMDER2drbVUYiokVjudI0WLVqgX79+3HMn8mIsd6qVw+HgnjuRF2O5U60cDgdOnDjBr/8l8lIsd6rVlcvuce+dyDux3KlWgwcPBsByJ/JWLHeqVVhYGKKjo1nuRF7KULmLSIqI7BWRXBH5fR1jJonIbhHZJSKLzI1JVnA4HPzEDJGXqrfcRSQQwBwAowH0B/CgiPSvMaY3gJcA3KaqMQCmuiEreZjD4cDBgwdx5swZq6MQUQMZ2XMfAiBXVQ+qaimAJQAm1BjzHwDmqOopAFDVfHNjkhWuHFTNycmxOAkRNZSRcu8C4Gi1+3mVj1XXB0AfEVkvIptEJMWsgGQdh8MBgF9DQOSNggyMkVoeq/l1gUEAegNIBtAVwFoRuVlVT1+1IJFUAKkA0L179waHJc/q1KkTIiMjeVCVyAsZ2XPPA9Ct2v2uAGpegy0PwJeqellVfwSwFxVlfxVVnaeqTlV1hoeHNzYzeVBsbCzLncgLGSn3LQB6i0gPEWkG4AEAy2qM+QLAcAAQkQ6omKY5aGZQsobD4cCePXtQUlJidRQiaoB6y11VywBMAZAGYA+Apaq6S0ReE5HxlcPSABSKyG4A3wJ4XlUL3RWaPCc2Nhbl5eXYsWOH1VGIqAGMzLlDVVcAWFHjsRnVbiuAaZU/5EOuHFTNzs7GkCFDLE5DREbxDFW6rhtvvBFhYWH8OCSRl2G503WJCAYPHsxyJ/IyLHeql8PhwPbt21FeXm51FCIyiOVO9Ro8eDAuXLiA/fv3Wx2FiAxiuVO9rnz9L6dmiLwHy53q1a9fPzRr1ownMxF5EZY71atZs2aIiYnhnjuRF2G5kyFXLphdcUoDEdkdy50MGTx4MAoKCnDixAmroxCRASx3MoQHVYm8C8udDBk0aBAAXjCbyFuw3MmQNm3aIDo6mnvuRF6C5U6G8WsIiLwHy50MczgcyM3NxdmzZ62OQkT1YLmTYVcOqm7fvt3iJERUH5Y7GVb9u92JyN5Y7mRYZGQkwsPDWe5EXoDlToaJSNWZqkRkbyx3ahCHw4Fdu3ahtLTU6ihEdB0sd2oQh8OBy5cvY9euXVZHIaLrYLlTg/CgKpF3YLlTg/Tq1QutW7dmuRPZHMudGiQgIACDBg1iuRPZHMudGszhcGDbtm1wuVxWRyGiOrDcqcEcDgfOnTuH3Nxcq6MQUR1Y7tRgVw6q8kvEiOyL5U4NFhMTg+DgYM67E9kYy50arFmzZujfvz/LncjGWO7UKLxgNpG9GSp3EUkRkb0ikisiv7/OuF+IiIqI07yIZEcOhwP5+fm8YDaRTdVb7iISCGAOgNEA+gN4UET61zIuFMDvAGSYHZLsh2eqEtmbkT33IQByVfWgqpYCWAJgQi3j/gzgrwAumpiPbOrKhTu2bt1qcRIiqo2Rcu8C4Gi1+3mVj1UREQeAbqq63MRsZGOhoaHo06cPsrKyrI5CRLUwUu5Sy2NVR9FEJADALAD/We+CRFJFJFNEMgsKCoynJFtyOp0sdyKbMlLueQC6VbvfFcDxavdDAdwMYI2IHAIwFMCy2g6qquo8VXWqqjM8PLzxqckW4uLikJeXh/z8fKujEFENRsp9C4DeItJDRJoBeADAsitPquoZVe2gqlGqGgVgE4DxqprplsRkG3FxcQDAvXciG6q33FW1DMAUAGkA9gBYqqq7ROQ1ERnv7oBkX1c+McNyJ7KfICODVHUFgBU1HptRx9jkpscib9CmTRseVCWyKZ6hSk0SFxfHcieyIZY7NUlcXByOHj3Kg6pENsNypybhQVUie2K5U5PwoCqRPbHcqUnCwsLQu3dvljuRzbDcqcl4UJXIflju1GRXDqryKyWI7IPlTk3Gg6pE9sNypyaLjY0FAGRm8hsniOyC5U5NFhYWhr59+2Lz5s1WRyGiSix3MkV8fDwyMjJ4TVUim2C5kyni4+ORn5+PQ4cOWR2FiMByJ5PEx8cDADIyeAldIjtguZMpBgwYgJCQEJY7kU2w3MkUwcHBcDqd2LRpk9VRiAgsdzJRfHw8srOzUVpaanUUIr/HcifTxMfH49KlS9i2bZvVUYj8HsudTMODqkT2wXIn03Tr1g2RkZGcdyeyAZY7mUZEqk5mIiJrsdzJVPHx8cjNzUVhYaHVUYj8GsudTMV5dyJ7YLmTqZxOJwICAjjvTmQxljuZKjQ0FAMHDsT69eutjkLk11juZLrExERs3LiRJzMRWYjlTqZLTExESUkJr8xEZCGWO5nu9ttvBwB8//33Fich8l8sdzJdx44d0bdvX6xdu9bqKER+i+VObpGYmIh169ahvLzc6ihEfslQuYtIiojsFZFcEfl9Lc9PE5HdIrJdRP6fiNxoflTyJomJiThz5gx27NhhdRQiv1RvuYtIIIA5AEYD6A/gQRHpX2NYNgCnqg4E8CmAv5odlLxLYmIiAM67E1nFyJ77EAC5qnpQVUsBLAEwofoAVf1WVS9U3t0EoKu5McnbdOvWDVFRUSx3IosYKfcuAI5Wu59X+VhdHgewsrYnRCRVRDJFJLOgoMB4SvJKiYmJ+P7776GqVkch8jtGyl1qeazWf60i8isATgD/XdvzqjpPVZ2q6gwPDzeekrxSYmIiCgoKsHfvXqujEPkdI+WeB6BbtftdARyvOUhERgD4A4DxqnrJnHjkzTjvTmQdI+W+BUBvEekhIs0APABgWfUBIuIAMBcVxZ5vfkzyRr169UJERATWrFljdRQiv1NvuatqGYApANIA7AGwVFV3ichrIjK+cth/A2gN4F8ikiMiy+pYHPkREcHIkSOxevVquFwuq+MQ+ZUgI4NUdQWAFTUem1Ht9giTc5GPSElJwccff4zMzEwMGTLE6jhEfoNnqJJb3XXXXRARrFq1yuooRH6F5U5u1aFDB9xyyy1YubLWT8cSkZuw3MntRo8ejc2bN/O6qkQexHInt0tJSYHL5cLq1autjkLkN1ju5Ha33HIL2rdvz3l3Ig9iuZPbBQYG4q677sKqVav4kUgiD2G5k0ekpKTgp59+wrZt26yOQuQXWO7kEaNGjQIAfmqGyENY7uQRERERiI2NxbJlPHmZyBNY7uQx999/PzIyMnDw4EGroxD5PJY7ecwDDzwAAFiyZInFSYh8H8udPCYqKgq33norFi9ebHUUIp/HciePevDBB7Fz505eOJvIzVju5FH3338/AgMDufdO5GYsd/KoTp064c4778SSJUt4bVUiN2K5k8c99NBD+PHHH5GRkWF1FCKfxXInj5s4cSKaN2+ORYsWWR2FyGcZuhITkZnatGmD8ePHY+HChfjLX/6C1q1bW5rH5XJhz5492LhxI06dOoUWLVqgRYsW6NmzJ2699VY0b97c0nxmO378eNVXMJ86dQolJSWIiopCnz59cNNNN6Ft27ZWRzSdqmLHjh04cOAA8vPzkZ+fj+bNm6NHjx6IiopC//790apVK6tjmktVLfmJi4tT8l8bN25UAPrOO+9YlmHDhg163333aVhYmAKo9adVq1Y6fvx4XbBggZaWllqWtal27dqlL7zwgg4cOLDObQWgIqJJSUk6d+5cLSwstDp2k5SXl+vKlSs1NTVVu3Tpct3tDgkJ0YkTJ+o///lPPXv2rNXRrwtAphroWJY7Web222/X7t27e7w0V69erYmJiQpA27dvr6mpqbpgwQLdt2+fFhcXa35+vh46dEiXLVumTz31lEZFRSkAvfHGG/X999/XixcvejRvU6xfv17vvvtuBaDBwcE6fPhwfeONNzQjI0MPHz6sxcXFeunSJf3hhx902bJl+sorr+hNN91UNf7JJ5/UvLw8qzejQUpLS/Xjjz/WmJgYBaChoaF633336fz583Xr1q167NgxvXTpkp4+fVq3bdumn3/+uT799NPauXNnBaBhYWH6yiuv6KlTp6zelFqx3Mn2li9frgD0448/9sj6ioqKdPLkyQpAu3btqrNmzdLi4uJ6f8/lcuny5cs1Pj5eAWhUVJSuWrXKA4kbb//+/Tp69GgFoDfccIP+6U9/0oKCAkO/63K5NCsrS5944gkNCgrSkJAQnTZtmhYVFbk5ddN9+eWX2rNnTwWgN998sy5cuFAvXbpk6HfLy8t13bp1eu+99yoAbdOmjb766qtaUlLi5tQNw3In2ysvL9eYmBgdMGCAulwut67rq6++0s6dO2tgYKD+8Y9/bNTet8vl0rS0NO3bt68C0MmTJxsuTE+5cOGCzpgxQ5s3b66hoaH65ptv6rlz5xq9vIMHD+ojjzyiAQEB2qlTJ/3kk0/c/lo1xoEDB3TcuHEKQGNiYnTZsmVaXl7e6OXl5OToxIkTFYD26tVLV69ebWLapmG5k1dYsGCBAtCVK1e6ZfllZWX60ksvVe3JZWZmNnmZJSUlOn36dA0ODtbw8HBdvny5CUmbLiMjo2pK5aGHHtLjx4+btuytW7dqXFycAtBx48bp0aNHTVt2U5SXl+vbb7+tISEh2rp1a33rrbdMneZbvXq19urVSwHoww8/bIvjECx38gqXLl3SLl266LBhw0zfIywqKtKUlBQFoKmpqabPlW/fvr3qAOXTTz9t2Z/vpaWlOn36dA0MDNSuXbtqenq6W9Zz+fJlfeutt7Rly5batm1bXbRokVvWY9SBAweqjp3cfffdbjs2cOXNPCgoSCMjI3XFihVuWY9RLHfyGvPmzVMA+t5775m2zJ07d2p0dLQGBwfr3LlzTVtuTSUlJTp16tSqvwx27NjhtnXVZvfu3VV71L/+9a89chBw//79mpCQoAD0l7/8pcf3Zl0ul86dO1dbtWqlbdq00QULFnhkqigrK6vqIG1qaqpln6phuZPXcLlcOnLkSG3VqpUeOHCgycv797//raGhoRoREaEbNmwwIWH9VqxYoR07dtSQkBCdM2eO28vG5XLp7NmzNSQkRG+44Qb97LPP3Lq+mi5fvqwzZ87UoKAgjYiI8NjU1PHjx6sOFN9555165MgRj6z3iosXL+rzzz+vIqI9evTQ7777zqPrV2W5k5c5cuSItmnTRpOSkhp9IMzlcukbb7yhIqJxcXEenxc+efJkVfHcfffdps55V3fw4EEdOXKkAtDRo0e7bT1GbN26VQcMGKAA9LHHHtPTp0+7ZT0ul0sXLFig7du31xYtWujs2bObdMC0qdatW6fR0dEqIjp16tQmHbRuKJY7eZ358+c3+sSmwsJCnTBhQtVUwfnz592QsH4ul0vfeecdDQkJ0bZt2+oHH3xg2l58WVmZvv3229qyZUtt3bq1vvvuu7b45MrFixf1pZde0oCAAI2IiNCFCxeammvfvn16xx13KABNSEjQH374wbRlN8W5c+f0t7/9rQLQbt266eeff+6R18PUcgeQAmAvgFwAv6/l+eYAPql8PgNAVH3LZLlTTS6XS8eNG6dBQUENmidfv369duvWTYODg3XWrFm2KLx9+/ZpUlKSAtDk5OQmfUrH5XLpypUr1eFwKAAdO3asx6cjjNi8ebM6nU4FoImJibp169YmLS8/P1+fe+45bd68uYaFhel7771n6d56XdatW1f118vYsWM1JyfHreszrdwBBAI4AKAngGYAtgHoX2PMbwG8X3n7AQCf1LdcljvV5vTp01WfcJk6daqWlZXVOfbo0aOampqqgYGB2rNnT92yZYsHk9avvLxc582bp+3atauaqsnKymrQ769evVqHDRtWdfLUkiVLbPHmVZeysjKdO3eutm/fXgHoqFGj9JtvvmlQ5ry8PH355Ze1VatWGhAQoJMnT7Z06smI0tJSffPNNzU0NFQB6Pjx4zUjI8Mt6zKz3BMApFW7/xKAl2qMSQOQUHk7CMDPAOR6y2W5U10uX76szzzzjALQpKQk/fDDD/XYsWOqWvHxxm+//VafffZZbd68uQYHB+vTTz/ttrleM5w5c0b//Oc/a9u2bRWADhgwQKdPn65btmy56gxZl8ulJ0+e1DVr1uizzz5bdTp8ZGSkzpkzx/CZlnZQVFSkM2fO1I4dOyoA7du3r06bNk1Xr16txcXFV5V9cXGx5uTk6Jw5czQxMVFFpGp6bffu3RZuRcMVFRXpq6++WvWG3qdPH33++ed17dq1+vPPP5vyxmxmuf8CwAfV7k8G8H9qjNkJoGu1+wcAdLjeclnuVJ+5c+dqp06dqr7cKTw8vOp2QECAPvroo/rjjz9aHdOw06dP66xZszQxMVEDAgKu+nKyqKgobdWqVdVjwcHBOn78eF20aJFeuHDB6uiNVlJSonPnztURI0Zos2bNqravefPmGhkZqREREVd9gVe/fv30tdde03379lkdvUnOnj2r7777ro4cOVKDgoKu+oKy6OhoXbx4caOXbbTcpWJs3UTkfgCjVPU3lfcnAxiiqk9XG7Orckxe5f0DlWMKaywrFUAqAHTv3j3u8OHD1103kapi+/btSE9Px+7du9GvXz8MGjQIDocDHTt2tDpeoxUUFODrr7/G0aNHcfLkSeTn5yM8PBzR0dHo2bMnEhIS0K5dO6tjmur8+fNYs2YNdu7ciaKiIhQWFkJV0atXL/Tq1Qs333wz+vbtCxGxOqqpTp8+jTVr1uDQoUPIy8vDsWPH8Pjjj2PEiBGNWp6IZKmqs95xBso9AcCfVHVU5f2XAEBV/6vamLTKMRtFJAjASQDhep2FO51OzczMNLQxRERUwWi5G7kS0xYAvUWkh4g0Q8UB02U1xiwD8Ejl7V8A+OZ6xU5ERO5V75WYVLVMRKag4qBpIID5qrpLRF5DxdzPMgD/A+BjEckFUISKNwAiIrKIocvsqeoKACtqPDaj2u2LAO43NxoRETUWL5BNROSDWO5ERD6I5U5E5INY7kREPojlTkTkg+o9icltKxYpANDYU1Q7oOL7a7yZt2+Dt+cHvH8bmN96VmzCAcGGAAADpklEQVTDjaoaXt8gy8q9KUQk08gZWnbm7dvg7fkB798G5reenbeB0zJERD6I5U5E5IO8tdznWR3ABN6+Dd6eH/D+bWB+69l2G7xyzp2IiK7PW/fciYjoOmxd7iKSIiJ7RSRXRH5fy/PNReSTyuczRCTK8ynrZiD/oyJSICI5lT+/sSJnXURkvojki8jOOp4XEflb5fZtF5FYT2esj4FtSBaRM9Vegxm1jbOKiHQTkW9FZI+I7BKRZ2oZY9vXwWB+u78GISKyWUS2VW7Dq7WMsV8XGblckxU/cNOFuW2W/1HUuGShnX4AJAKIBbCzjufHAFgJQAAMBZBhdeZGbEMygOVW57xO/kgAsZW3QwHsq+X/I9u+Dgbz2/01EACtK28HA8gAMLTGGNt1kZ333IcAyFXVg6paCmAJgAk1xkwA8FHl7U8B3Cn2uUaXkfy2pqrfo+L7+esyAcA/tMImAG1FJNIz6YwxsA22pqonVHVr5e1iAHsAdKkxzLavg8H8tlb53/Vc5d3gyp+aBytt10V2LvcuAI5Wu5+Ha/+nqBqjqmUAzgC4wSPp6mckPwDcV/mn9Kci0s0z0UxjdBvtLqHyT+6VIhJjdZi6VP6p70DFnmN1XvE6XCc/YPPXQEQCRSQHQD6A1apa52tgly6yc7nX9q5X893SyBirGMn2bwBRqjoQwNf433d+b2Hn//5GbUXF6dyDAMwG8IXFeWolIq0BfAZgqqqerfl0Lb9iq9ehnvy2fw1UtVxVBwPoCmCIiNxcY4jtXgM7l3segOp7sl0BHK9rTOWFucNgnz/B682vqoWqeqny7t8BxHkom1mMvEa2pqpnr/zJrRVXHAsWkQ4Wx7qKiASjohj/qar/t5Yhtn4d6svvDa/BFap6GsAaACk1nrJdF9m53L39wtz15q8xLzoeFfOR3mQZgF9XflpjKIAzqnrC6lANISIRV+ZGRWQIKv5NFFqb6n9VZvsfAHtU9e06htn2dTCS3wteg3ARaVt5uwWAEQB+qDHMdl1k6BqqVlAvvzC3wfy/E5HxAMpQkf9RywLXQkQWo+KTDB1EJA/AK6g4mARVfR8V19UdAyAXwAUAj1mTtG4GtuEXAJ4SkTIAJQAesPofZQ23AZgMYEflnC8AvAygO+AVr4OR/HZ/DSIBfCQigah441mqqsvt3kU8Q5WIyAfZeVqGiIgaieVOROSDWO5ERD6I5U5E5INY7kREPojlTkTkg1juREQ+iOVOROSD/j975Qt3AUQKeQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f9099fa16d8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xx = x.level()\n",
    "def H(w): return xx[0] + 2*sum([xx[i]*cos(i*w) for i in range(1,len(xx))])\n",
    "w  = np.linspace(0, pi, 100)\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(w, [H(wi) for wi in w], 'k')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\"><img alt=\"Creative Commons License\" style=\"border-width:0\" src=\"https://i.creativecommons.org/l/by/4.0/80x15.png\" /></a><br />This work is licensed under a <a rel=\"license\" href=\"http://creativecommons.org/licenses/by/4.0/\">Creative Commons Attribution 4.0 International License</a>. The **MOSEK** logo and name are trademarks of <a href=\"http://mosek.com\">Mosek ApS</a>. The code is provided as-is. Compatibility with future release of **MOSEK** or the `Fusion API` are not guaranteed. For more information contact our [support](mailto:support@mosek.com). \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}