{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sizing of clock meshes\n",
    "\n",
    "Original by Lieven Vanderberghe.<br>\n",
    "Adapted to CVX by Argyris Zymnis, 12/4/2005.<br>\n",
    "Modified by Michael Grant, 3/8/2006.<br>\n",
    "Adapted to CVXPY, with cosmetic modifications, by Judson Wilson, 5/26/2014.<br>\n",
    "\n",
    "Topic References:\n",
    "\n",
    "* Section 4, L. Vandenberghe, S. Boyd, and A. El Gamal<br>\n",
    "  \"Optimal Wire and Transistor Sizing for Circuits with Non-Tree Topology\"\n",
    "\n",
    "## Introduction\n",
    "\n",
    "We consider the problem of sizing a clock mesh, so as to minimize the\n",
    "total dissipated power under a constraint on the dominant time constant.\n",
    "The numbers of nodes in the mesh is $N$ per row or column (thus $n=(N+1)^2$\n",
    "in total). We divide the wire into m segments of width $x_i$, $i = 1,\\dots,m$\n",
    "which is constrained as $0 \\le x_i \\le W_{\\mbox{max}}$. We use a pi-model of each wire\n",
    "segment, with capacitance $\\beta_i x_i$ and conductance $\\alpha_i x_i$.\n",
    "Defining $C(x) = C_0+x_1 C_1 + x_2 C_ 2 + \\cdots + x_m C_m$ we have that the dissipated\n",
    "power is equal to $\\mathbf{1}^T C(x) \\mathbf{1}$. Thus to minimize the\n",
    "dissipated power subject to a constraint in the widths and a constraint\n",
    "in the dominant time constant, we solve the SDP\n",
    "    \\begin{array}{ll}\n",
    "    \\mbox{minimize}   & \\mathbf{1}^T C(x) \\mathbf{1}   \\\\\n",
    "    \\mbox{subject to} &  T_{\\mbox{max}} G(x) - C(x) \\succeq 0 \\\\\n",
    "                      & 0 \\le x_i \\le W_{\\mbox{max}}.\n",
    "    \\end{array}\n",
    "\n",
    "## Import and setup packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cvxpy as cp\n",
    "import numpy as np\n",
    "import scipy as scipy\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Show plots inline in ipython.\n",
    "%matplotlib inline\n",
    "\n",
    "# Plot properties.\n",
    "plt.rc('text', usetex=True)\n",
    "plt.rc('font', family='serif')\n",
    "font = {'weight' : 'normal',\n",
    "        'size'   : 16}\n",
    "plt.rc('font', **font)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computes the step response of a linear system.\n",
    "def simple_step(A, B, DT, N):\n",
    "    n  = A.shape[0]\n",
    "    Ad = scipy.linalg.expm((A * DT))\n",
    "    Bd = (Ad - np.eye(n)).dot(B)\n",
    "    Bd = np.linalg.solve(A, Bd)\n",
    "    X  = np.zeros((n, N))\n",
    "    for k in range(1, N):\n",
    "        X[:, k] = Ad.dot(X[:, k-1]) + Bd;\n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate problem data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# Circuit parameters.\n",
    "#\n",
    "\n",
    "dim=4            # Grid is dimxdim (assume dim is even).\n",
    "n=(dim+1)**2     # Number of nodes.\n",
    "m=2*dim*(dim+1)  # Number of wires.\n",
    "                 #   0 ... dim(dim+1)-1 are horizontal segments\n",
    "                 #     (numbered rowwise);\n",
    "                 #   dim(dim+1) ... 2*dim(dim+1)-1 are vertical\n",
    "                 #     (numbered columnwise)\n",
    "beta = 0.5       # Capacitance per segment is twice beta times xi.\n",
    "alpha = 1        # Conductance per segment is alpha times xi.\n",
    "G0 = 1           # Source conductance.\n",
    "C0 = np.array([ (      10,     2,     7,     5,     3),\n",
    "                (       8,     3,     9,     5,     5),\n",
    "                (       1,     8,     4,     9,     3),\n",
    "                (       7,     3,     6,     8,     2),\n",
    "                (       5,     2,     1,     9,    10) ])\n",
    "wmax = 1        # Upper bound on x.\n",
    "\n",
    "#\n",
    "# Build capacitance and conductance matrices.\n",
    "#\n",
    "\n",
    "CC = np.zeros((dim+1, dim+1, dim+1, dim+1, m+1))\n",
    "GG = np.zeros((dim+1, dim+1, dim+1, dim+1, m+1))\n",
    "\n",
    "# Constant terms.\n",
    "# - Reshape order='F' is fortran order to match original\n",
    "#   version in MATLAB code.\n",
    "CC[:, :, :, :, 0] = np.diag(C0.flatten(order='F')).reshape(dim+1, dim+1,\n",
    "                                                           dim+1, dim+1, order='F').copy()\n",
    "zo13 = np.zeros((2, 1, 2, 1))\n",
    "zo13[:,0,:,0] = np.array([(1, 0), (0, 1)])\n",
    "zo24 = np.zeros((1, 2, 1, 2))\n",
    "zo24[0,:,0,:] = zo13[:, 0, :, 0]\n",
    "pn13 = np.zeros((2, 1, 2, 1))\n",
    "pn13[:,0,:,0] = np.array([[1, -1], [-1, 1]])\n",
    "pn24 = np.zeros((1, 2, 1, 2))\n",
    "pn24[0, :, 0, :] = pn13[:, 0, :, 0]\n",
    "\n",
    "for i in range(dim+1):\n",
    "    # Source conductance.\n",
    "    # First driver in the middle of row 1.\n",
    "    GG[int(dim/2), i, int(dim/2), i, 0] = G0\n",
    "    for j in range(dim):\n",
    "        # Horizontal segments.\n",
    "        node = 1 + j + i * dim\n",
    "        CC[j:j+2, i, j:j+2, i, node] = beta * zo13[:, 0, :, 0]\n",
    "        GG[j:j+2, i, j:j+2, i, node] = alpha * pn13[:, 0, :, 0]\n",
    "        # Vertical segments.\n",
    "        node = node + dim * ( dim + 1 )\n",
    "        CC[i, j:j+2, i, j:j+2, node] = beta * zo24[0, :, 0, :]\n",
    "        GG[i, j:j+2, i, j:j+2, node] = alpha * pn24[0, :, 0, :]\n",
    "\n",
    "# Reshape for ease of use.\n",
    "CC = CC.reshape((n*n, m+1), order='F').copy()\n",
    "GG = GG.reshape((n*n, m+1), order='F').copy()\n",
    "\n",
    "#\n",
    "# Compute points the tradeoff curve, and the three sample points.\n",
    "#\n",
    "\n",
    "npts    = 50\n",
    "delays  = np.linspace(50, 150, npts)\n",
    "xdelays = [50, 100]\n",
    "xnpts   = len(xdelays)\n",
    "areas   = np.zeros(npts)\n",
    "xareas  = dict()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve problem and display results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Point 1 of 50 on the tradeoff curve (Tmax = 50.0)\n",
      "CVXOPT failed, trying robust KKT\n",
      "Point 2 of 50 on the tradeoff curve (Tmax = 52.04081632653061)\n",
      "Point 3 of 50 on the tradeoff curve (Tmax = 54.08163265306123)\n",
      "Point 4 of 50 on the tradeoff curve (Tmax = 56.12244897959184)\n",
      "Point 5 of 50 on the tradeoff curve (Tmax = 58.16326530612245)\n",
      "Point 6 of 50 on the tradeoff curve (Tmax = 60.20408163265306)\n",
      "Point 7 of 50 on the tradeoff curve (Tmax = 62.244897959183675)\n",
      "Point 8 of 50 on the tradeoff curve (Tmax = 64.28571428571429)\n",
      "Point 9 of 50 on the tradeoff curve (Tmax = 66.3265306122449)\n",
      "Point 10 of 50 on the tradeoff curve (Tmax = 68.36734693877551)\n",
      "Point 11 of 50 on the tradeoff curve (Tmax = 70.40816326530611)\n",
      "Point 12 of 50 on the tradeoff curve (Tmax = 72.44897959183673)\n",
      "Point 13 of 50 on the tradeoff curve (Tmax = 74.48979591836735)\n",
      "Point 14 of 50 on the tradeoff curve (Tmax = 76.53061224489795)\n",
      "Point 15 of 50 on the tradeoff curve (Tmax = 78.57142857142857)\n",
      "Point 16 of 50 on the tradeoff curve (Tmax = 80.61224489795919)\n",
      "Point 17 of 50 on the tradeoff curve (Tmax = 82.65306122448979)\n",
      "Point 18 of 50 on the tradeoff curve (Tmax = 84.6938775510204)\n",
      "Point 19 of 50 on the tradeoff curve (Tmax = 86.73469387755102)\n",
      "Point 20 of 50 on the tradeoff curve (Tmax = 88.77551020408163)\n",
      "Point 21 of 50 on the tradeoff curve (Tmax = 90.81632653061224)\n",
      "Point 22 of 50 on the tradeoff curve (Tmax = 92.85714285714286)\n",
      "Point 23 of 50 on the tradeoff curve (Tmax = 94.89795918367346)\n",
      "Point 24 of 50 on the tradeoff curve (Tmax = 96.93877551020408)\n",
      "Point 25 of 50 on the tradeoff curve (Tmax = 98.9795918367347)\n",
      "Point 26 of 50 on the tradeoff curve (Tmax = 101.0204081632653)\n",
      "Point 27 of 50 on the tradeoff curve (Tmax = 103.06122448979592)\n",
      "Point 28 of 50 on the tradeoff curve (Tmax = 105.10204081632654)\n",
      "Point 29 of 50 on the tradeoff curve (Tmax = 107.14285714285714)\n",
      "Point 30 of 50 on the tradeoff curve (Tmax = 109.18367346938776)\n",
      "Point 31 of 50 on the tradeoff curve (Tmax = 111.22448979591837)\n",
      "Point 32 of 50 on the tradeoff curve (Tmax = 113.26530612244898)\n",
      "Point 33 of 50 on the tradeoff curve (Tmax = 115.3061224489796)\n",
      "Point 34 of 50 on the tradeoff curve (Tmax = 117.34693877551021)\n",
      "Point 35 of 50 on the tradeoff curve (Tmax = 119.38775510204081)\n",
      "Point 36 of 50 on the tradeoff curve (Tmax = 121.42857142857143)\n",
      "Point 37 of 50 on the tradeoff curve (Tmax = 123.46938775510205)\n",
      "Point 38 of 50 on the tradeoff curve (Tmax = 125.51020408163265)\n",
      "Point 39 of 50 on the tradeoff curve (Tmax = 127.55102040816327)\n",
      "Point 40 of 50 on the tradeoff curve (Tmax = 129.59183673469389)\n",
      "Point 41 of 50 on the tradeoff curve (Tmax = 131.6326530612245)\n",
      "Point 42 of 50 on the tradeoff curve (Tmax = 133.67346938775512)\n",
      "Point 43 of 50 on the tradeoff curve (Tmax = 135.71428571428572)\n",
      "Point 44 of 50 on the tradeoff curve (Tmax = 137.75510204081633)\n",
      "Point 45 of 50 on the tradeoff curve (Tmax = 139.79591836734693)\n",
      "Point 46 of 50 on the tradeoff curve (Tmax = 141.83673469387756)\n",
      "Point 47 of 50 on the tradeoff curve (Tmax = 143.87755102040816)\n",
      "Point 48 of 50 on the tradeoff curve (Tmax = 145.9183673469388)\n",
      "Point 49 of 50 on the tradeoff curve (Tmax = 147.9591836734694)\n",
      "Point 50 of 50 on the tradeoff curve (Tmax = 150.0)\n",
      "Particular solution 1 of 2 (Tmax = 50)\n",
      "CVXOPT failed, trying robust KKT\n",
      "Solution 1:\n",
      "Vertical segments:\n",
      "[[0.65284441 0.4391586  0.52378143 0.47092764 0.2363529 ]\n",
      " [0.99999993 0.85353862 0.99999992 0.93601078 0.56994586]\n",
      " [0.92325575 0.29557654 0.80041338 0.99999998 0.99999997]\n",
      " [0.41300012 0.13553757 0.26699524 0.67049218 0.88916807]]\n",
      "Horizontal segments:\n",
      "[[1.96487539e-01 1.40591789e-01 9.70591442e-08 7.79376843e-08\n",
      "  5.27429285e-08]\n",
      " [7.07446433e-02 6.38430105e-02 1.02136471e-07 8.59913722e-08\n",
      "  6.28906472e-08]\n",
      " [6.05807467e-09 1.16285450e-08 3.91561390e-08 9.48052913e-02\n",
      "  1.58096913e-01]\n",
      " [3.82528741e-07 4.85708568e-07 5.75578696e-07 8.39862772e-02\n",
      "  5.38639181e-02]]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 29648 characters]\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Particular solution 2 of 2 (Tmax = 100)\n",
      "Solution 2:\n",
      "Vertical segments:\n",
      "[[0.2687881  0.04368684 0.17122095 0.133796   0.07360396]\n",
      " [0.41346231 0.08016135 0.30642705 0.2224136  0.1484946 ]\n",
      " [0.25755998 0.08016077 0.11200259 0.38352317 0.28159768]\n",
      " [0.13439419 0.04368697 0.02445701 0.24083502 0.24534599]]\n",
      "Horizontal segments:\n",
      "[[ 1.53896782e-09 -5.18600578e-10 -9.75218556e-10 -5.19196383e-10\n",
      "   1.57176577e-09]\n",
      " [ 9.30752726e-10 -9.56673760e-10 -1.35065528e-09 -9.96797753e-10\n",
      "   1.03852376e-09]\n",
      " [ 9.35404466e-10 -9.12219313e-10 -2.22938358e-10 -7.91865186e-10\n",
      "   1.51304362e-09]\n",
      " [ 1.31975762e-09 -8.50790152e-10 -1.39421076e-09 -8.33247519e-10\n",
      "   1.27680128e-09]]\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 32508 characters]\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 18472 characters]\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Iterate over all points, and revisit specific points\n",
    "for i in range(npts + xnpts):\n",
    "    # First pass, only gather minimal data from all cases.\n",
    "    if i < npts:\n",
    "        delay = delays[i]\n",
    "        print( ('Point {} of {} on the tradeoff curve ' \\\n",
    "              + '(Tmax = {})').format(i+1, npts, delay))\n",
    "    # Second pass, gather more data for specific cases,\n",
    "    # and make plots (later).\n",
    "    else:\n",
    "        xi = i - npts\n",
    "        delay = xdelays[xi]\n",
    "        print( ('Particular solution {} of {} ' \\\n",
    "              + '(Tmax = {})').format(xi+1, xnpts, delay))\n",
    "\n",
    "    #\n",
    "    # Construct and solve the convex model.\n",
    "    #\n",
    "\n",
    "    # Variables.\n",
    "    xt = cp.Variable(shape=(m+1)) # Element 1 of xt == 1 below.\n",
    "    G = cp.Variable((n,n), symmetric=True)  # Symmetric constraint below.\n",
    "    C = cp.Variable((n,n), symmetric=True)  # Symmetric constraint below.\n",
    "    \n",
    "    # Objective.\n",
    "    obj = cp.Minimize(cp.sum(C))\n",
    "\n",
    "    # Constraints.\n",
    "    constraints = [ xt[0] == 1,\n",
    "                    G == G.T,\n",
    "                    C == C.T,\n",
    "                    G == cp.reshape(GG*xt, (n,n)),\n",
    "                    C == cp.reshape(CC*xt, (n,n)),\n",
    "                    delay * G - C == cp.Variable(shape=(n,n), PSD=True),\n",
    "                    0 <= xt[1:],\n",
    "                    xt[1:] <= wmax,\n",
    "                  ]\n",
    "\n",
    "    # Solve problem (use CVXOPT instead of SCS to match original results;\n",
    "    # cvxopt produces lower objective values as well, but is much slower)\n",
    "    prob = cp.Problem(obj, constraints)\n",
    "    try:\n",
    "        prob.solve(solver=cp.CVXOPT)\n",
    "    except cp.SolverError:\n",
    "        print(\"CVXOPT failed, trying robust KKT\")\n",
    "        prob.solve(solver=cp.CVXOPT, kktsolver='robust')\n",
    "            \n",
    "    if prob.status not in [cp.OPTIMAL, cp.OPTIMAL_INACCURATE]:\n",
    "        raise Exception('CVXPY Error')\n",
    "    \n",
    "    # Chop off the first element of x, which is \n",
    "    # constrainted to be 1\n",
    "    x = xt.value[1:] \n",
    "\n",
    "    # First pass, only gather minimal data from all cases.\n",
    "    if i < npts:\n",
    "        areas[i] = sum(x)\n",
    "    # Second pass, gather more data for specific cases,\n",
    "    # and make plots.\n",
    "    else:\n",
    "        xareas[xi] = sum(x)\n",
    "\n",
    "        #\n",
    "        # Print display sizes.\n",
    "        #\n",
    "\n",
    "        print('Solution {}:'.format(xi+1))\n",
    "        print('Vertical segments:')\n",
    "        print(x[0:dim*(dim+1)].reshape(dim, dim+1, order='F').copy())\n",
    "        print('Horizontal segments:')\n",
    "        print(x[dim*(dim+1):].reshape(dim, dim+1, order='F').copy())\n",
    "\n",
    "        #\n",
    "        # Determine and plot the step responses.\n",
    "        #\n",
    "\n",
    "        A = -np.linalg.inv(C.value).dot(G.value)\n",
    "        B = -A.dot(np.ones(n))\n",
    "        T = np.linspace(0, 500, 2000)\n",
    "        Y = simple_step(A, B, T[1], len(T))\n",
    "        indmax = -1\n",
    "        indmin = np.inf\n",
    "        for j in range(Y.shape[0]):\n",
    "            inds = np.amin(np.nonzero(Y[j, :] >= 0.5)[0])\n",
    "            if ( inds > indmax ):\n",
    "                indmax = inds\n",
    "                jmax = j\n",
    "            if ( inds < indmin ):\n",
    "                indmin = inds\n",
    "                jmin = j\n",
    "\n",
    "        tthres = T[indmax]\n",
    "        GinvC  = np.linalg.solve(G.value, C.value)\n",
    "        tdom   = max(np.linalg.eig(GinvC)[0])\n",
    "        elmore = np.amax(np.sum(GinvC.T, 0))\n",
    "        plt.figure(figsize=(8, 8))\n",
    "        plt.plot( T, np.asarray(Y[jmax,:]).flatten(), '-',\n",
    "                  T, np.asarray(Y[jmin,:]).flatten() )\n",
    "        plt.plot( tdom   * np.array([1, 1]), [0, 1], '--',\n",
    "                  elmore * np.array([1, 1]), [0, 1], '--',\n",
    "                  tthres * np.array([1, 1]), [0, 1], '--' )\n",
    "        plt.xlim([0, 500])\n",
    "        plt.ylim([0, 1])\n",
    "        plt.text(tdom, 0.92, 'd')\n",
    "        plt.text(elmore, 0.88, 'e')\n",
    "        plt.text(tthres, 0.96, 't')\n",
    "        plt.text( T[600], Y[jmax, 600], 'v{}'.format(jmax+1))\n",
    "        plt.text( T[600], Y[jmin, 600], 'v{}'.format(jmin+1))\n",
    "        plt.title(('Solution {} (Tmax={}), fastest ' \\\n",
    "                    + 'and slowest step responses').format(xi+1, delay), fontsize=16)\n",
    "        plt.show()\n",
    "\n",
    "#\n",
    "# Plot the tradeoff curve.\n",
    "#\n",
    "\n",
    "plt.figure(figsize=(8, 8))\n",
    "ind = np.isfinite(areas)\n",
    "plt.plot(areas[ind], delays[ind])\n",
    "plt.xlabel('Area')\n",
    "plt.ylabel('Tdom')\n",
    "plt.title('Area-delay tradeoff curve')\n",
    "# Label the specific cases.\n",
    "for k in range(xnpts):\n",
    "    plt.text(xareas[k], xdelays[k], '({})'.format(k+1))\n",
    "plt.show()"
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
