{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Equilibrium of a system of weights connected by strings/springs\n",
    "\n",
    "In this notebook we show how to solve the following problem: Find the equlibrium of a system of masses connected by a system of strings, with some masses being assigned fixed coordinates (attached to the wall, say). See the next picture.\n",
    "\n",
    "![](basic.png)\n",
    "\n",
    "Suppose we have $n$ masses with weights $w_1,\\ldots,w_n$, and the length of the string between $i$ and $j$ is $\\ell_{ij}$ for some set $L$ of pairs of indices $(i,j)$ (we assume $\\ell_{ij}$ is not defined if there is no connection). The strings themselves have no mass. We also have a set $F$ of indices such that the $i$-th point is fixed to have coordinates $f_i$ if $i\\in F$. The equilibrium of the system is a configuration which minimizes potential energy. With this setup we can write our problem as:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{array}{ll}\n",
    "minimize & g\\cdot \\sum_i w_ix_i^{(2)} \\\\\n",
    "s.t.     & \\|x_i-x_j\\|\\leq \\ell_{ij},\\ ij\\in L \\\\\n",
    "         & x_i = f_i,\\ i\\in F\n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "\n",
    "where $x\\in (\\mathbf{R}^n)^2$, $x_i^{(2)}$ denotes the second (vertical) coordinate of $x_i$ and $g$ is the gravitational constant.\n",
    "\n",
    "Here is a sample problem description."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "w = [0.0, 1.1, 2.2, 0.0, 2.1, 2.2, 0.2]\n",
    "l = {(0,1): 1.0, (1,2): 1.0, (2,3): 1.0, (1,4): 1.0, (4,5): 0.3, (5,2): 1.0, (5,6): 0.5, (1,3): 8.0}\n",
    "f = {0: (0.0,1.0), 3: (2.0,1.0)}\n",
    "g = 9.81"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can formulate the problem using Mosek Fusion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mosek.fusion import *\n",
    "\n",
    "# w - masses of points\n",
    "# l - lengths of strings\n",
    "# f - coordinates of fixed points\n",
    "# g - gravitational constant\n",
    "def stringModel(w, l, f, g):\n",
    "    n, m = len(w), len(l)\n",
    "    starts = [ lKey[0] for lKey in l.keys() ]\n",
    "    ends = [ lKey[1] for lKey in l.keys() ]\n",
    "\n",
    "    M = Model(\"strings\")\n",
    "\n",
    "    # Coordinates of points\n",
    "    x = M.variable(\"x\", [n, 2])\n",
    "\n",
    "    # A is the signed incidence matrix of points and strings\n",
    "    A = Matrix.sparse(m, n, list(range(m))+list(range(m)), starts+ends, [1.0]*m+[-1.0]*m)\n",
    "\n",
    "    # ||x_i-x_j|| <= l_{i,j}\n",
    "    c = M.constraint(\"c\", Expr.hstack(Expr.constTerm(list(l.values())), Expr.mul(A, x)), \n",
    "        Domain.inQCone() )\n",
    "\n",
    "    # x_i = f_i for fixed points\n",
    "    for i in f:\n",
    "        M.constraint(x.slice([i,0], [i+1,2]), Domain.equalsTo(list(f[i])).withShape([1,2]))\n",
    "\n",
    "    # sum (g w_i x_i_2)\n",
    "    M.objective(ObjectiveSense.Minimize, \n",
    "        Expr.mul(g, Expr.dot(w, x.slice([0,1], [n,2]))))\n",
    "\n",
    "    # Solve\n",
    "    M.solve()\n",
    "    if M.getProblemStatus(SolutionType.Interior) == ProblemStatus.PrimalAndDualFeasible:\n",
    "        return x.level().reshape([n,2]), c.dual().reshape([m,3])\n",
    "    else:\n",
    "        return None, None\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is a quick description of how we use vectorization to deal with all the conic constraints in one go. The matrix $A$ is the incidence matrix between the masses and the strings, with coefficients $+1, -1$ for the two endpoints of each string. It is chosen so that the product $Ax$ has rows of the form\n",
    "\n",
    "$$\n",
    "(x_i^{(1)} - x_j^{(1)}, x_i^{(2)} - x_j^{(2)})\n",
    "$$\n",
    "\n",
    "for all pairs $i,j$ for which $\\ell_{ij}$ is bounded. Stacking the values of $\\ell$ in the left column produces a matrix with each row of the form\n",
    "\n",
    "$$\n",
    "(\\ell_{ij}, x_i^{(1)} - x_j^{(1)}, x_i^{(2)} - x_j^{(2)})\n",
    "$$\n",
    "\n",
    "and a conic constraint is imposed on all the rows, as required.\n",
    "\n",
    "The objective and linear constraints show examples of slicing the variable $x$.\n",
    "\n",
    "The function returns the coordinates of the masses and the values of the dual conic variables. A zero dual value indicates that a particular string is hanging loose, and a nonzero value means it is fully stretched. \n",
    "\n",
    "All we need now is to define a display function and we can look at some plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "# x - coordinates of the points\n",
    "# c - dual values of string length constraints\n",
    "# d - pairs of points to connect\n",
    "def display(x, c, d):\n",
    "    import matplotlib.pyplot as plt\n",
    "    fig, ax = plt.subplots()\n",
    "    # Plot points\n",
    "    ax.scatter(x[:,0], x[:,1], color=\"r\")\n",
    "    # Plot fully stretched strings (nonzero dual value) as solid lines, else dotted lines\n",
    "    for i in range(len(c)):\n",
    "        col = \"b\" if c[i][0] > 1e-4 else \"b--\"\n",
    "        ax.plot([x[d[i][0]][0], x[d[i][1]][0]], [x[d[i][0]][1], x[d[i][1]][1]], col)\n",
    "    ax.axis(\"equal\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xt4VNW5x/Hvi9wMIlqxVRQIKBdRqVqO9yoVe2r1KCqgeKISL02tUuuxNy0FqxVrr9pWa8ULARqVijdUlKqI12INLQqoKCBBLireEAyCkHX+eCcSISGTzGTvmdm/z/PMMzN7trNfJuO796z1rrUshICIiCRLq7gDEBGR6Cn5i4gkkJK/iEgCKfmLiCSQkr+ISAIp+YuIJJCSv4hIAin5i4gkkJK/iEgCtY47gIZ07tw5FBcXxx2GiEhemT179nshhF0b2y9nk39xcTGVlZVxhyEiklfMrCqd/dTsIyKSQEr+IiIJpOQvIpJASv4iIgmk5C8ikkBK/iIiCZSV5G9mt5vZu2Y2r4HXzcz+ZGYLzexlMzsoG8eVNFRUQHExtGrl9xUVcUckspm+n7HJ1pV/OXDcNl7/NtArdSsDbsrScWVbKiqgrAyqqiAEvy8r0/9gkhv0/YxVVpJ/COFp4INt7DIYmBjcLGAnM9s9G8eWbRg1CqqreZKB/JsDfVt1tW8XiVvq+zmJEqro5tv0/YxMVG3+ewBv1Xm+LLXtC8yszMwqzaxy1apVEYVWwJYuZT1tOZuJlFLOBtp8vl0kbq9Vbc8xPMEIJnA1P9/8gr6fkYgq+Vs928JWG0IYF0IYEEIYsOuujU5NIY3p1o12bOAvXMhc+vMrLv98u0hcPvoI/u//YH9e5jmOwAj8jLGbd9D3MxJRJf9lQNc6z/cEVkR07OQaOxaKijiRhziDOxjLKOa1H+DbRWLw7LPQqxf88Y8w7LDl1GCUcQs9SE1HU1Sk72dEokr+U4GzU1U/hwKrQwgrIzp2cpWUwLhx0L07f+QSOrVaw7ldHmHj6SVxRyYJs3at3/frB4cfDrNnQ1G/YrZrsx2j9pgAZtC9u39fS/T9jEK2Sj3vBP4J9DGzZWZ2npldYGYXpHaZBiwGFgK3ABdm47iShpISWLKEXcO7/LliF15c3Jnrr487KEmKN9+EIUPg6KOhpga+9CV44AHo2BHKy+G739uOPZfN8heXLFHij1BWpnQOIZzRyOsBuCgbx5LmO/10uPNOGD0aBg/2n98iLWHtWvjVr+D3v4fttoPLL4eNG6FtW3/9qqv88eWXxxtnkmmEb4KYwU03Qbt2cP75frElkm2vvAK9e8M118Bpp8Hrr8PPf7458b/2mpfyX3QR7LZbvLEmmZJ/wnTp4ldjTz8NN98cdzRSSD76yO/33hsGDoRZs2DiRNhji6LuK6+E7beHn/wk8hClDiX/BDr3XDj2WP+fTyXVkqlly+DMM2G//by5p21buOMOOOSQrfedNw8mT4aLLwZVc8dLyT+BzLyooqYGvvtdH1kv0lTr1sEvfwl9+sCUKVBa6lP0bMsVV3hn749+FEmIsg1K/gnVo4d3yD36KEyaFHc0km9WroS+fWHMGDj+eHj1Vbj6ai/Tb8h//gP33usDvL70pehilfop+SfYyJFec33JJfD223FHI/ngvff8frfd4MQTYeZMuPtuv5hozJgxsPPOnvwlfkr+CdaqFdx2m8+lNXJk3NFILnvnHa8Q69HD2/jN4IYbvH4/Hf/6Fzz0kDf3dOrUsrFKepT8E65vX2+Hvecev4nUtX49/Pa3PiZkwgTvI+rYsenvM2YM7LILfP/72Y9RmkfJX/jRj+DAA73u+oNtTcwtiVJdDf37e1XY0UfD/Pnwu981/cr9uedg+nT46U+bd+KQlqHkL7RpA7ffDu+/r/ZY8c5c8M7bM8/0ooAHH/SBW80xejR85St+cSG5Q8lfADjgAL8ymzgRHnkk7mgkDu+/730/3bvDiy/6ttGj4Vvfav57Pvmk3y6/fNuVQBI9JX/53OjRsM8+3q778cdxRyNR+ewz+POfvV3/ppvgO99Jr3qnMSF4W3+XLv6dktyi5C+fa9fOq3+WLYPLLos7GolCTQ0ceaSPuD3oIJgzB268ETp3zvy9H3vM5+8fNQrat8/8/SS7lPzlCw47DH7wA78CfOqpuKORllK7ZnqrVnDeeXD//Z6s998/O+8fgv+S7NbN319yj5K/bOXqq6FnT6/rrq6OOxrJptWr4cc/9iae2tLesjKf4tvqW2y1mR5+2Gv7R4/2X5SSe5T8ZSsdOsAtt8DChT4GQPLfpk3epNe7t8/qetZZ3tzTEmrb+nv2hBEjWuYYkjklf6nXMcd4x98f/rC58kPy16mn+i+5Xr3873nbbS03l/599/k8Pldc4WXEkpss5OiUjgMGDAiVlZVxh5Foq1fDvvv6fCyzZ29ejEPyQ1WVJ/h27WDqVG/CO/307DbvbKmmBr76VdiwwQeFtc7KWoHSFGY2O4QwoLH9dOUvDerUCf76V5+D/Zpr4o5G0vXJJ97W3revl3ACnHQSDB/esokffJK3efPgF79Q4s91Sv6yTf/zP/C//wtjx8LcuXFHI9sSgi+P2KePd9qfcopf6Udl0yZP+vvuG+1xpXmU/KVRf/yjN/2ce64vwi256Xvf8+kYdt/d59O54w7o2jW6499xh6/Pe+WVjS/qIvHTn0ga1bmzNx9UVsJ118UdjdS1YsXmyfhKS2H8eHjhBV+nIUqffeZJ/4AD/BeH5D4lf0nLaad5LfiYMfD663FHI59+6v0wvXt7+z7AoYemt5RiS5g0CRYtgquu0lV/vtCfSdJiBn/5i1eOnH++V3VI9ELwwVn77OPTJvz3f8Oll8Yb04YNnvT/67+8j0jyg5K/pK1LF6/7f+YZrwKS6P3ylzB0qM+L/8QTvibuXnvFG9Ptt3tZ6VVXtXw1kWSP6vylSULwKX7/+U8v6evePe6ICt+qVV6j3707LFniU25/5zu5UUr56aew994e27PPKvnnAtX5S4swg3Hj/CTw3e/6vbSMDRu8g71Xr80LoRQXe1VPLiR+8O/C8uX+i0SJP78o+UuTFRfDr37lS/NNnBh3NIVp2jSfYfPSS32m1d/9Lu6ItlZd7Z3OAwf6dCCSX5T8pVkuugiOOMKXfXz77bijKSy33gonnOCPH37Ym3n69o03pvrcdBO884639Uv+UfKXZmnVyicHq67W2qzZ8OGH8Mor/vi003xg3dy5cPzx8cbVkLVr4dpr4ZvfhK9/Pe5opDmU/KXZ+vTx4fz33gtTpsQdTX7auNErp3r18rl3QoAdd/SVtXJ5Ir0//xnee8/b+iU/KflLRn70I1/+76KLfAFwSd+MGf7Zfe973r4/aVJ+dJquXg2//a03TR1ySNzRSHMp+UtGWrf2Ou8PPvD2f0nPtGkwaBCsWeO/mmbM8KmQ88H113szldr681tWkr+ZHWdmC8xsoZlttfS3mZWa2Sozm5O6nZ+N40pu+OpXfcH3SZM8qUn91qzxpQ3BR+beeKO38w8Zkh9X/OBJ/w9/gJNP9l8tkr8yTv5mth1wI/BtoB9whpn1q2fXySGEA1K3WzM9ruSWn/8c+vXz2v+PP447mtxSUwPl5T4Pz4kn+sCo1q3hwgth++3jjq5pfv97//teeWXckUimsnHlfzCwMISwOISwAbgLGJyF95U80q6dV/8sXw4//Wnc0eSO55/3dvFzzvHxEQ8+CO3bxx1V87z3nlchnXYa9O8fdzSSqWwk/z2At+o8X5batqUhZvaymU0xs3pnGTezMjOrNLPKVatWZSE0idKhh8Ill3j1ysyZcUcTv5de8rEQK1bA3/7mc+wffHDcUTXfb37jpb2/+EXckUg2ZCP519daueWg/weB4hBCf+BxYEJ9bxRCGBdCGBBCGLDrrrtmITSJ2tVXQ8+ePvNndXXc0USvutonXAPvCykvhwULoKQkv6c6fvttuOEGX9Vtn33ijkayIRtfx2VA3Sv5PYEVdXcIIbwfQlifenoL8LUsHFdyUFGRj1BdtMjn/k+KEOCuu3wk7gkn+MhXgBEjYIcd4o0tG6691ucaStLftNBlI/m/CPQysx5m1hYYDkytu4OZ7V7n6UnAq1k4ruSob3wDysp8UrLa6pZCNnu2j3I94wxf9ewf/4CvfCXuqLJn+XJvyjv7bB+MJoUh4+QfQtgIjASm40n97yGE+WZ2lZmdlNrtYjObb2YvARcDpZkeV3Lbb37ja8meey6sX9/4/vlq5UqfeO311+GWW+DFF+Goo+KOKruuucYXZ69dMUwKg+bzlxbz8MO+stOYMYVVGrh+PTz6qC9rCXD//f5rp1OneONqCVVVfrV/7rlawCdfaD5/id0JJ3hH5zXXwMsvxx1N5kKABx6Afff1QU7z5vn2k08uzMQP3oFv5ktGSmFR8pcWdf31sPPOfuW4cWPc0TTfvHk+g+XJJ/uYhunTYb/94o6qZS1aBOPH+8C9rvUWZ0s+U/KXFtW5s5cIzp7t0wLko3XrfMGSf/8b/vQnmDPHp2codL/8JbRpA5dfHnck0hKU/KXFDRvmV8xXXOEdo/ngs8/gjjt8aobtt4e//x3eeAO+/31PiIVuwQKfq+nCC73jXgqPkr+0ODP4y198WoPzzvOEmsumT/cBWiUlXrYJvkzhLrvEG1eUrrzST3qaqqNwKflLJHbf3Zt9nn3Wl//LRa+/7hOvHXecX/k/8AB861txRxW9efN8wNr3vw9f/nLc0UhLUamnRCYET6zPP+8Jpnv3uCParKbGZyVdscLr2S++2Dt2k2joUP/F8+abyfq1UyhU6ik5xwxuvtlPAmVlfh+nTZtgwgSfj6dVK5987Y034Mc/Tm7inzMH7rnHJ+hT4i9sSv4SqeJinyfmH//wxBuXp56CAQOgtNQ7dsGfF9K0DM1xxRWw005w6aVxRyItTclfInfhhXDkkb7s48qV0R57yRKvPho40JeenDzZO6HFp6aYOhV++EM/AUhhU/KXyLVq5TN/rlvnC79H2fxTVubTTlx1Fbz2mi9Mki9LKLa0MWO8qecHP4g7EomCkr/Eok8fT8D33ecLmLeUmhpvy6/9hXHDDV7DPnp0/i2h2JKef97nK/rJT6Bjx7ijkSgo+UtsLr0UvvY1GDkS3n8/++//wgtw+OFw1lkwbpxv691bUxXUZ/RoL+u86KK4I5GoKPlLbFq3httv97b3Sy7J3vsuX+5zzx96qM9KWV6u6Yi3ZeZMmDEDLrsMOnSIOxqJipK/xKp/f/jZz7xpZtq07LznFVd4R+7ll/vArREj8nsJxZYUgrf1d+kCF1wQdzQSJQ3yktht2AAHHQSrV8P8+bDjjk3770Pw2vTevf1k8vbbXrvfs2fLxFtIHnvMJ6m74QY1+RQKDfKSvNG2rTf/rFjhHY5NMWeOl20OG+YJDGC33ZT40xGCN4d17Qrnnx93NBI1JX/JCQcf7HX/N9/sbdCNefddL9s86CD/tXDTTT55nKRv2jTvFB89OrkjmpNMzT6SM6qrfTbNEHzlr6Kihve96iqfb37kSG+z3nnn6OIsBCH4iOYPP/TS1yRMU50UavaRvFNU5IO/Fi3aujonBHjoIXjiCX/+wx/C3Llw3XVK/M3xwAO+OM2YMUr8SaXkLznl6KO96uT662p4YfeToVUrXt1jEN/+6gpOPNGTPXhJYt++8caar2pqPOn37g1nnhl3NBIXJX/JOb8+aDJ7sJySt39NWbiJ/VdMZ9bcIq47s5L77os7uvw3ZYr/arriCh9rIcmkNn/JPcXFTK76L4ZzNxBoy3p2ZRU7t/mEjgP6suOOPgVB3fv6tm35Wrt2msdn0ybYf3//HF5+GbbbLu6IJNvSbfPXeV9yz9KlDGMpF/ABG2lNKeWsYUfWfNaRjzv05aOPYOlSWLMGPv7Y79O5hmnTpnknjfru27Zt+Y+hJdx5J7z6qq9JrMSfbEr+knu6daNVVRV/4Iecy3jO4C4O55++9NdjQ7bavabGK4VqTwR17+vbVve1996DxYs3b1u7Nr0Q27Vr3kljy20dO0bU4VpRwcafjeHKpY/Sv81Ghnz6b6AkggNLrlLyl9wzdiyUlTG0egojuYFySjm86CXfXo9WrWCHHfyWqZoaPwGke/Kou+2dd2Dhws2vVVend8z27bPza6Rjxwau5isqoKyMSdWnsZBe3P/ZYFpd8Lj3+JXoBJBUavOX3FRRAaNGMaLqSu63U1h56zSKzh0ed1RNsnHjF08kzTmh1N6vW5feMYuK6jkxvPAY7de9x70MYU+WsZi9MPBfUkuWtOAnIHFQm7/kt5ISKCmh9EmYeAzc3344/xt3TE3UurWviJWNVbE++8xPAg2dNLZ18li6rjNLOZCNtOUQ/sXnfd5Ll2YemOQtXflLTqupgb32gl69fN1faYbiYgZX/ZF/cTBL6UobNvl2XfkXJI3wlYLQqpVPyfz44/DWW3FHk5/e+cnveZgTOJuJmxN/UVGDfSiSDEr+kvPOPttLOSdNijuS/FTx6RA20ZrSLv/wAv/u3X1pM3X2JpqafSQvDBzoUz4vWKCBWk0Rgq9x0KEDzJoVdzQSBTX7SEEpLYU33vCFxiV9//kPzJvnn59IXVlJ/mZ2nJktMLOFZnZZPa+3M7PJqddfMLPibBxXkmPoUL96LS+PO5L8Ul7uA9JOPz3uSCTXZJz8zWw74Ebg20A/4Awz67fFbucBH4YQ9gauA36d6XElWXbYwU8AkyenP3gq6dav9+ESp5yiaa9la9m48j8YWBhCWBxC2ADcBQzeYp/BwITU4ynAIDO13ErTlJZ67bpm9kzPQw/BBx+oyUfql43kvwdQtwhvWWpbvfuEEDYCq4FdtnwjMyszs0ozq1y1alUWQpNCctRRUFyspp90lZdDly5w7LFxRyK5KBvJv74r+C1LiNLZhxDCuBDCgBDCgF133TULoUkhqa35f+IJDU5tzNtvwyOPeJmsZu+U+mQj+S8DutZ5viewoqF9zKw10An4IAvHloQZMUI1/+moqPC5+9XkIw3JRvJ/EehlZj3MrC0wHJi6xT5TgRGpx0OBGSFXBxhITuvRw2v+y8vTm8M/iULwz+eww6BPn7ijkVyVcfJPteGPBKYDrwJ/DyHMN7OrzOyk1G63AbuY2ULgUmCrclCRdJWW+tTJqvmv3+zZqu2XxmmEr+SdtWtht91g+HC49da4o8k9I0fCbbfBypXZmVFU8otG+ErB2mEHGDbMlyL85JO4o8kt69fDHXd4bb8Sv2yLkr/kJdX81+/BB+HDD9XkI41T8pe89PWve+evav6/aPx42HNPGDQo7kgk1yn5S16qrfmfMQOqquKOJjesXAmPPqrafkmPkr/kLdX8f9Hf/uYrn40Y0fi+Ikr+kreKi+Eb31DNP2yu7T/8cOjdO+5oJB8o+UteKy2FRYvguefijiRelZXwyivq6JX0KflLXhsyxEs/k97xO348bL89nHZa3JFIvlDyl7zWoYNq/j/9FO68E049FTp1ijsayRdK/pL3amv+77037kjiMXUqfPSRmnykaZT8Je8deST07Jncpp/ycuja1Tu/RdKl5C95L8k1/ytWwPTpqu2XplPyl4JQW9s+cWK8cURt0iSv7VeTjzSVkr8UhO7d4ZhjklXzX1vbf+SRsPfecUcj+UbJXwpGaSksXgzPPht3JNH417/gtdd01S/No+QvBePUU5NV819e7rX9w4bFHYnkIyV/KRgdOvggpyTU/NfW9g8ZAjvuGHc0ko+U/KWglJb6Sl+FXvN///2wejWcc07ckUi+UvKXglJb8z9+fNyRtKzycujWzRezF2kOJX8pKGZ+9f/kk7BkSdzRtIzly+Gxx7y8tZX+D5Zm0ldHCs7ZZ/tJoFBr/mtr+zVvv2RCyV8KTt2a/5qauKPJrtra/qOOgr32ijsayWdK/lKQSkvhzTcLr+Z/1ixYsEC1/ZI5JX8pSKecAh07Fl7Nf3k5FBXB0KFxRyL5TslfClLdmv+1a+OOJjvWrYO77vLE37Fj3NFIvlPyl4JVWuqDvQql5v/+++Hjj9XkI9mh5C8F64gjvFO0UJp+xo/3ReuPPjruSKQQKPlLwSqkmv+33oLHH1dtv2SPvkZS0Gpr/idMiDuSzEya5GWeZ58ddyRSKJT8paB16waDBnnyz9ea/9ra/qOP9qkrRLJByV8KXm3N/zPPxB1J8/zzn/DGG+rolexS8peCl+81/+PHe+mqavslm5T8peAVFcHpp8Pdd+dfzX91NUye7Au27LBD3NFIIcko+ZvZl8zsMTN7I3W/cwP7bTKzOanb1EyOKdIctTX/99wTdyRNc999sGaNmnwk+zK98r8MeCKE0At4IvW8PutCCAekbidleEyRJjv8cF/kPN+afsrLoUcP+PrX445ECk2myX8wUFtENwE4OcP3E2kRtTX/M2d6528+WLoUnnhCtf3SMjL9Sn0lhLASIHX/5Qb2a29mlWY2y8waPEGYWVlqv8pVq1ZlGJrIF511Vn7N8z9xopd5at5+aQkWQtj2DmaPA7vV89IoYEIIYac6+34YQtiq3d/MuoQQVphZT2AGMCiEsGhbxx0wYECorKxM598gkrZvfhMWLoRFi3L7ajoE6NXLxynMmBF3NJJPzGx2CGFAY/s1+vUPIRwbQtivntsDwDtmtnvqgLsD7zbwHitS94uBmcCBTfi3iGTNOef4VA9PPx13JNv23HN+glJHr7SUTK99pgK1P0pHAA9suYOZ7Wxm7VKPOwNHAK9keFyRZjn5ZNhxx9zv+C0v99LOIUPijkQKVabJ/1rgm2b2BvDN1HPMbICZ3ZraZx+g0sxeAp4Erg0hKPlLLGpr/qdMyd2a/08+8XUIhg3zwV0iLSGj5B9CeD+EMCiE0Ct1/0Fqe2UI4fzU4+dDCPuHEL6aur8tG4GLNFdtzf+UKXFHUr977/Xa/nPOiTsSKWQ53OUl0jIOO8w7U3O16ae83CdwO/LIuCORQqbkL4lTW/P/1FOweHHc0XxRVZVX95SWepwiLUXJXxIpV2v+J070uDRvv7Q0JX9JpK5d4dhjc2ue/9p5+485Brp3jzsaKXRK/pJYtTX/Tz0VdyTumWe8GUq1/RIFJX9JrFyr+S8v93UHTjkl7kgkCZT8JbG23x6GD/eSzzVr4o1l7Vqv7T/tNNX2SzSU/CXRSkt9wZS4a/7vvdfHHqjJR6Ki5C+Jduih0Lt3/E0/48f7egNHHBFvHJIcSv6SaLU1/08/7ROpxeHNN32dAdX2S5SU/CXx4q75r63tP+useI4vyaTkL4m3554+z38cNf81NX7cQYN87n6RqCj5i+A1/1VV0df8P/OMN/uoo1eipuQvAgweDJ06Rd/xO368jzVQbb9ETclfhFTN/9feYMqkatbYjlBcDBUVLXrMtWu9xPT0032dAZEoKfmLAFRUUPrcd6gORdzNUG8DKitr0RPAlCmq7Zf4NLqAe1y0gLtEqriYUFVFDxbTlg1UcCb7Mp+i7l/2CYBawMCBsGIFLFigEk/JnnQXcG8dRTAiOW/pUgzoyWKeZBAH8yJGDb2q3qD/MOjfH/bf3++Li6FVhr+ZFy/2zuWxY5X4JR5K/iLgdZZVVdzNMBaxF8voysv0Z27RIcyZ04d77vEpl8EXVq89EdS932mn9A83YYLm7Zd4qdlHBLxtv6zMJ/qpVVQE48ZBSQlr18L8+TB3Lrz88ubbhx9u3r1bt80ng9pb797QuvUXj1Pzs5/Tc+mT9Gm/lOm3vgUlJZH9M6XwqdlHpClqE/CoUbB0qWfysWM/377DDnDIIX6rFQIsX771CWH6dNi40fdp2xb69UudDDbOpv89d/Hx+gOpophffXo5lE394vFFIqIrf5Es27ABXnvtiyeEuXO9c3ezGqrpwPZ86st2tVCnsiSPrvxFYtK27eZmn7res12Zy36cwR1sorUnfvBfGiIRU/IXiUjn7h34RtVM+rLgiy9oUh+JgQZ5iURl7Nith/IWFfl2kYjpyl8kKrWduue1h/Wptv46ncoiUVLyF4lSSQnckno8c0mckUjCqdlHRCSBlPxFRBJIyV9EJIGU/EVEEkjJX0QkgZT8RUQSKKPkb2bDzGy+mdWYWYNzSZjZcWa2wMwWmtllmRxTREQyl+mV/zzgVODphnYws+2AG4FvA/2AM8ysX4bHFRGRDGQ0yCuE8CqAbXspooOBhSGExal97wIGA69kcmwREWm+KNr89wDeqvN8WWrbVsyszMwqzaxy1apVEYQmIpJMjV75m9njwG71vDQqhPBAGseo72dBvYsIhBDGAePA5/NP471FRKQZGk3+IYRjMzzGMqBrned7Aisa2FdERCIQRbPPi0AvM+thZm2B4cDUCI4rIiINyLTU8xQzWwYcBjxsZtNT27uY2TSAEMJGYCQwHXgV+HsIYX5mYYuISCYyrfa5D7ivnu0rgOPrPJ8GTMvkWCIikj0a4SsikkBK/iIiCaTkLyKSQEr+IiIJpOQvIpJASv4iIgmk5C8ikkBK/iIiCaTkLyKSQEr+IiIJpOQvIpJASv4iIgmk5C8ikkBK/iIiCaTkLyKSQEr+IiIJpOQvIpJASv4iIgmk5C8ikkBK/iJRqqiAWbPgqZlQXOzPRWKg5C8SlYoKKCuD9Z/686oqf64TgMRAyV8kKqNGQXX1F7dVV/t2kYi1jjsAkcRYuhSAA5hT73aRKCn5i0SlWzeoquJ6/m/r7SIRU7OPSFTGjoWioi9uKyry7SIRU/IXiUpJCYwbB927g5nfjxvn20UipmYfkSiVlCjZS07Qlb+ISAIp+YuIJJCSv4hIAin5i4gkkJK/iEgCKfmLiCSQhRDijqFeZrYKqIo7jjR0Bt6LO4g0KM7sUpzZpTizp3sIYdfGdsrZ5J8vzKwyhDAg7jgaozizS3Fml+KMnpp9REQSSMlfRCSBlPwzNy7uANKkOLNLcWaX4oyY2vxFRBJIV/4iIgmk5N9EZvYlM3vMzN5I3e/cwH6bzGxO6jY1wviOM7MFZrbQzC6r5/V2ZjY59foLZlYcVWxbxNFYnKVmtqrOZ3h+DDHebmbvmtm8Bl43M/tT6t/wspl2VLEYAAADgElEQVQdFHWMqTgai3Ogma2u81mOiTrGVBxdzexJM3vVzOab2Q/q2Sf2zzTNOHPiM81ICEG3JtyA3wCXpR5fBvy6gf3WxhDbdsAioCfQFngJ6LfFPhcCf009Hg5MztE4S4EbYv5bHwUcBMxr4PXjgUcAAw4FXsjROAcCD8X5Wabi2B04KPW4I/B6PX/32D/TNOPMic80k5uu/JtuMDAh9XgCcHKMsWzpYGBhCGFxCGEDcBceb111458CDDIzizBGSC/O2IUQngY+2MYug4GJwc0CdjKz3aOJbrM04swJIYSVIYR/px6vAV4F9thit9g/0zTjzHtK/k33lRDCSvAvCfDlBvZrb2aVZjbLzKI6QewBvFXn+TK2/tJ+vk8IYSOwGtglkujqiSGlvjgBhqR++k8xs67RhNYk6f47csFhZvaSmT1iZvvGHUyqufFA4IUtXsqpz3QbcUKOfaZNpZW86mFmjwO71fPSqCa8TbcQwgoz6wnMMLO5IYRF2YmwQfVdwW9ZzpXOPi0tnRgeBO4MIaw3swvwXyvHtHhkTZMLn2U6/o0P+V9rZscD9wO94grGzHYA7gEuCSF8vOXL9fwnsXymjcSZU59pc+jKvx4hhGNDCPvVc3sAeKf2Z2jq/t0G3mNF6n4xMBO/emhpy4C6V8h7Aisa2sfMWgOdiL7JoNE4QwjvhxDWp57eAnwtotiaIp3PO3YhhI9DCGtTj6cBbcyscxyxmFkbPKFWhBDurWeXnPhMG4szlz7T5lLyb7qpwIjU4xHAA1vuYGY7m1m71OPOwBHAKxHE9iLQy8x6mFlbvEN3y0qjuvEPBWaEVA9WhBqNc4t23pPwdtdcMxU4O1WhciiwurZJMJeY2W61/TpmdjD+//37McRhwG3AqyGEPzSwW+yfaTpx5spnmgk1+zTdtcDfzew8YCkwDMDMBgAXhBDOB/YBbjazGvxLcW0IocWTfwhho5mNBKbjFTW3hxDmm9lVQGUIYSr+pZ5kZgvxK/7hLR1XM+O82MxOAjam4iyNOk4zuxOv6uhsZsuAK4A2qX/DX4FpeHXKQqAaOCfqGNOMcyjwPTPbCKwDhsdwwge/CDoLmGtmc1LbfgZ0qxNrLnym6cSZK59ps2mEr4hIAqnZR0QkgZT8RUQSSMlfRCSBlPxFRBJIyV9EJIGU/EVEEkjJX0QkgZT8RUQS6P8BtCToTxvdspMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0e0c912828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x,c = stringModel(w, l, f, g)\n",
    "\n",
    "if x is not None:\n",
    "    display(x, c, list(l.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How about we find a discrete approximation to the [catenary](#https://en.wikipedia.org/wiki/Catenary):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8lNW9x/HPyUYSCJuASlhVRBOXWuK+N/G63Nal2l5tVGz1qmBre+1yq/beKq231i5Wq2CtpWhJa6ltLbUqGlxRUIKCGhRFFIEoskNIQrZz/5gn8zzPZEgmMOsz3/frlZfJPL9MTkbyzck55znHWGsREZFgyUl1A0REJP4U7iIiAaRwFxEJIIW7iEgAKdxFRAJI4S4iEkAKdxGRAFK4i4gEkMJdRCSA8lL1hYcNG2bHjRuXqi8vIpKRlixZstFaO7y3upSF+7hx46irq0vVlxcRyUjGmNWx1GlYRkQkgBTuIiIBpHAXEQkghbuISAAp3EVEAkjhLiISQAp3EZEAUriLiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgAKdxFRAJI4S4iEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgGkcBcRCSCFu4hIACncRUQCSOEuIhJACncRkQBSuIuIBJDCXUQkgBTuIiIBpHAXEQkghbuISAAp3EVEAkjhLiISQAp3EZEAUriLiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgA9RruxpiZxphPjTFv7ea6McbcbYxZaYx5wxjz2fg3czdqamDcOMjJCf23piZpX1pEpFcpzKhYeu6zgLN6uH42MMF5uxqYsffNikFNDVx+OaxeDdaG/nv55Qp4EUkPNTVwxRX+jLriiqRlVK/hbq19AdjcQ8l5wEM2ZBEw2Bizf7wauFvXXAOdnXzIWP7BF0KPdXbCZZcl/EuLiPRq8mRob+dvXMBaSkOPtbfDtdcm5cvHY8y9FFjj+Xit81g3xpirjTF1xpi6DRs27N1X3bkTgEqe5nwe5T3GhR63Fqqq9u65RUT2RlUVdHSwkAou5K+cw7/ca42NSWlCPMLdRHnMRiu01t5vra2w1lYMHz48Dl8ajuI1IIeTecl9cP78uDy3iMgecTLobOYBhtN5JulNiEe4rwVGez4eBTTE4Xl7lhNq+p+4GOhkPfuzmRL3unrvIpIKTvasZAzbGAJ0chc3uNdNtP5w/MUj3OcClzurZo4DtllrP47D8/bsmmsAyAeO4HXAcCyvutfVexeRVHCy5xReAgynRfba02XM3RjzJ2AhMNEYs9YYc6Ux5lpjTFcLHwdWASuB3wJTE9Zar+nTw+8u5HjAspKJtHlrpianKSIiQDhzNjGQjykFOnmKM/w1nuxKJGNt1OHxhKuoqLB1dXV79yRTp8KM0MrLsazkIw6kgoUs5gS3JkXfn4hkIWfI5WDqeY8yDmcpb3CUe33KlL0Od2PMEmttRW91mX2HqudFWkwFYKnjWNq9NRp7F5FkcLKmmTze41DA8gpH+2uS1GuHTA93gMpKAEawlX35GMjhbB53r2vsXUSSwcmaE1kIGA7gPYq8XU0nq5Il88O9tjb87gJOACy1nOmv0di7iCSS02tvB17ns4RGESb5azxZlQyZH+4Q/o14EKsZykYgh/N5xL0+Izk7IohIlnJ67afyPJBDKR8xBM/NSknutUNQwt3zG3EeZwCWf3C+v0a9dxFJBCdb2oCXOQmwLOBEf02Se+0QlHAHGDAAgAqWMZAtQC6X8JB7Xb13EUkEJ1vOZB6Qwwg+YRzr3OtTpqSkWcEJ9/vuC7/7OGcDlof5ir9GvXcRiSdPpjxLFWB5IbLXnsQVMl7BCffqaigsBOBEXqWEbUAuk3nArVHvXUTiycmUf+MxIIdhfMpEPnCv9++fmnYRpHAHeMAN8sc4B7A8xBX+GvXeRSQePPfQPO2MFizw3kAJ8JvfJLdNHsEK9+pqKCsD4BQWUsJWIJfLmOnWqPcuIvHgrJA5k8dxe+2r3OtlZaFMSpFghTtAfX343a6x99lc7q9R711E9oZnhcxTnEnUFTKeLEqF4IU7hMfeT+IVSpyVM1/ij+519d5FZG84GVJFLV0rZCbyvns9BevaIwUz3D1j7086v1Uf4cv+GvXeRWRPeHrtL3A6YHkxDda1RwpmuHtWzpxAHYPZBORyDo+6Neq9i8iecLLjRF4CctifdRycJitkvIIZ7uDrvb/g3DX2BF/w7xip3ruI9IWn177YOUfiVSJ2303hChmv4Ia7Z+XM4axgOJ8AOXyOZ90a9d5FpC+czPgsSwHDWFYxivXu9RSvkPEKbriDb7Z6MZMIjY2d6j+tSfu9i0gsnKzYTiFvcQRgeYPD/DUpXiHjFexwh/Cs9Vg+ZiwfAIajWOpe137vIhILJysOpx4wHMYyBtLiXk+DFTJewQ93z6z16xwJWOo5gk0MdGvKy5PfLhHJHE5GvMdYPmI80Mlr3uPzIC1WyHgFP9wh/Bt1CI0c7oyVHckb7vXly6GmJjVtE5H0VlMTygjgOF4FDKfyLPnemjTrtUOmH5DdF87BtW1AAR2AoZ5DKOPd0PXcXGhv3+2ni0iWys+H9nae4SQqeQHoxJLnr0lijmbHAdl94eypnA9U8SRgnHWqjo4O9d5FxK+mJtzp+zyPA4ZqZvtrUrRfe2+yp+cO4d47gKEdyOGfnM3nmede7+xMbptEJH3l5UFHBw9wOf/JLAztdFLgr0lyhqrnHo3nN+x1/BqAi/ibe91a3dgkIiFTp4b+ogemcD8AP+YH/po07bVDtvXcIfybGCCPXXSQz4+5kZv5qVuTotdERNKI85f+9fySX/MtCmhmF/3911Pwl7567rvz4IPhd2dwDQD/w4/9Neq9i2Q3Twb8musBeIQL/TV/+EMyW9Rn2ddzh/DsN0B/ttHEQM7nL/zdu3Okeu8i2cvptZ/MMyzgdIawgc2McK+ncHWdeu49mTUr/O4zzpadj3Khf1sC9d5FspOzzcA2ilnAaYBlceTmYJ4RgHSVneHu2RL4WF5jFKuBHA7nLbdGm4qJZJ+amvA2A+W8DRgOpZ4D+citKShIm83BepKd4Q6+LYHrKQMsKyhjBQe4NaWlyW+XiKTO5MkAvMTRrGM00MkyDvfXzJzZ/fPSUPaGu2dL4IE0U8k8wHA0nnmAhgYNz4hkC8/SxzN4FjBcwmz/NgNptKVvb7JzQtXLc2NTDm1YcpnB1VyL27PX5KpIFnCy4BZ+wK1MI5c22unnr0mDLNCEaqw8G/78iJsAuI7p/hr13kWCzXOuw63cAsD9XOmvSeMblqKJqedujDkLuAvIBR6w1t4ecX0sMBMYDmwGLrXWru3pOdOm5w6QkxP+jVzMdpop4fP8nX/yRbcmDX5ji0gC1NTApZcCcAIvsJCTGcQmtjLMrUmjrUni1nM3xuQC9wJnA2XAJcaYsoiynwMPWWuPAKYBP+l7k1PIczPC85wCWB7jPLZT5NZoz3eRYLriCgDWsi8LnfOWu+3VnuY3LEUTy7DMMcBKa+0qa20r8DBwXkRNGdB1pNGzUa6nN8/SyKNZygTeBnKYyAq3Rnu+iwSPZ9fHrhOWjmERB7DGrcmQpY+RYgn3UvB+p6x1HvNaBuF7cy8ASowx++x985LItzSyHOjkE0bxD852ay67LPntEpHEcZY+/oavspWhQAevcIK/JkOWPkaKJdxNlMciB6C/A5xqjHkdOBVYB3S7N9cYc7Uxps4YU7dhw4Y+NzahPEsj84Gv8ysALuJRt8ZaHagtEhRVVeGlj1PDuz7e5K+prMzIXjvEFu5rgdGej0cBDd4Ca22DtfaL1tqjgJudx7ZFPpG19n5rbYW1tmL48OF70ewEqa8PL4f6Nd+mH820k8+FPOzWzJ+v4RmRTOe5E/U0aukkjwFs52bucGuMSbtzUfsilnBfDEwwxow3xhQAFwNzvQXGmGHGmK7nupHQypnM5Jk4eZZTAfgbX2Yjg9yar30t2a0SkXi66ioAVjGa5/kcYHmdI/01GTiJ6tVruFtr24GvA/OAt4E51tp6Y8w0Y8y5TtlpwApjzLvAvsBtCWpv4lVXh/Z8B46njkN5i9D+Eu+4Na2t6r2LZKqaGmhpAeBIlgGG41jAQax2a3JzM3Y4povuUI3Gs+41dKB26Ei+GVzJtfzerdPad5HM49zX8kNuZho/wtBBp3+TAZg9O23DXXeo7o3q6vCdq/nANP4XcCddwjS5KpJZysvBWtqBadwKwH3OoT1hGTyJ6qWee088+84MYhPbGcpRLOY1jnFr1HsXyQyev8gPYgXvczAj+Jj1jPTXpfnPtHru8eDZS2IFEwhNulTwsnfjfm0LLJIZnDXtczmb952f55XeLb4h4/aP6YnCvSfTp8PI0G/1/djMhcwB4DRecGsaGjQ8I5LuPGvaL3DuXZnCrymhxa0ZOTL0Mx8QCvferFsXfvcRLiaPFtoo5Dz+6tZo7btI+vKsaT+ZZ+mkgCIamc43/XWen/UgULjHwvOn2oucAsBcLmC1d6xOWxOIpCfnZ7OOI1nAqYSGVz/jr5k9O/ntSjCFeyymTw+tewWOo47jWAAYyrxr363Vvu8i6aaqKjxBeiKLAPh35jKRVW5NANa0R6Nwj5XntPOFnEIurTQxgMneE5t0qLZIenGGY87hUVopJJ8WHuN8f43nZztIFO6x8qx9B3iSfwPgIb7GGvZz67R6RiQ9OD+Lb3IIT3AuYFnkXcYMgVnTHo3CvS9qa8PDM1U8zyReBQyH8K5bo9UzIqlXVRX6WQQqeA0wVDKPz/KWWzNyZEZvDNYbhXtfef6Eq+M4cmijiQFUM8ut0eoZkdTxrI45k8ec4Zhmar1nM0DgVsdEUrj3VcTwzOPOP5g/cjkrGevWafWMSGo4P3uv8Fme4hwAFnKcvyZANyvtjsJ9T3iGZ85kPsfzImA4jOVujQ72EEm+0tLw6piTeBmAzzOXSbzh1hgTqJuVdkfhvqc8wzMvcyr5tLCLIs7gCbdGwzMiyTN1qmecfRHt9KOQJv4ZuTomw/dpj5XCfU9FDM8sZhIAtZzJ05zu1jn7WYhIgjlLkf/ERSzhGMDyNgf7awK8OiaSwn1veIZnjmQ5X+KPAJzJPLemo0PDMyKJ5vyMtQFfcY7FvI67GOc9ETTgq2MiKdz3lmd4Zg6XMoBtWPI5iBVujYZnRBLHszpmNGuAXPbhU+7hv/x1AV8dE0nhvrcihmc+ZCzQyftM4Cd8161z9pEWkThzVsf8F3ewnlKggzWM8tcEcO+Y3ijc46G2Nnzu6j5s5w6+A8BN/JT1DHXrdPeqSHw5q2PeYyy/cn7uZjGZItrdmoKCrBln91K4x8usWeF3v8udHOIcrD3Oe+iu7l4ViR/PXajlvAMYjmUhk4kYAp05M/ltSwMK93iJGJ55myPIYxct9OdknnXrNP4usvemTg2Psx/O67RRSD+aWcSJ/rosWh0TSeEeT7W14ZObAJZyBAALOJWZeO5Y1fi7yJ6rqQkve/wJN/AWRwKdrPLeIQ5ZtzomksI93tatCx+sXc67fJNfAnAls9hGsVun8XeRPeNMoK5jODfxcwB+xI2MZIO/LstWx0RSuCeC5w64X/EdxrAKyGF/PnFrGhqgvDz5bRPJZJ7tBcazmtC2H8v4AXf467JwdUwkhXsiRIy/r+YgcmmlmRKOcfa7AGD5ck2wisTKM4F6CG/SRhEFNPEmR/nrsnic3UvhnigR4+/LOQSwLOY4/s9ZsgVoglUkFp4blW7gDlZQDlg+YJy/rqwsq8fZvRTuieQZfz+YD7iFmwC4mZ+yggPcOk2wivTMGWdfyCTudDpHM7jKP84+ciTU16eidWlJ4Z5onvH3H3K7c3pTDmW8TZu3ThOsItE54+xtuIdcVzGPa4lYv57lE6iRFO6JVl3tOxigjuPoz3Y6KaDUu6mRbnAS6a68PDzOPoxNWPIYygaejjxVSROo3Sjck2H6dN8E66cMAzrYwH5U8aRbN39+6OYMEQl1dpaHDsA5moVsZwg5tPEJ+/rrNIEalcI9WTwTrMW08SInADCff+Mn3gnWGTM0wSrimUD9NrdTx7EAvEE5+d66LL9RqScK92TyTLCexKvcxK0A3MQdvMpn3DpNsEq2cyZQn+I0fsn3ALiT6ynnPX+dxtl3S+GebJ4J1tu4lZN4HoBjqWMHhW6dJlglWw0ZAtaynqGcSaj3/gX+zre4x1+ncfYeKdyTLeIGpxc5nRF8DOQynI1une5glWxUWgpbt9IGjKIByGEs7zOXC/11U6ZonL0XMYW7MeYsY8wKY8xKY8z3o1wfY4x51hjzujHmDWPMOfFvaoDU1oZutnCsp5R8mtlFf0r5yK1bvlwBL9nDszJmXz51Drhu5EMm+OsqK0OLFKRHvYa7MSYXuBc4GygDLjHGlEWU/QCYY609CrgY0Cvfm/p63x2s6xkBdNDAKE7iObdOAS/ZoLw8vDLmcF5nC8MwtLOFIf46TaDGLJae+zHASmvtKmttK/AwcF5EjQUGOu8PAu8CbtktzwTrEBpZ4KygeYlTuIr73DrtQSNB5lnyeBFznC18LfUcSqH3RKXBgzWB2gexhHspsMbz8VrnMa9bgEuNMWuBx4FvRHsiY8zVxpg6Y0zdhg0bopVkH88E64m8yr1cDcDvuJpfel9GrYGXIPIcunETt/JXLgLgz1zIoaz0127ZkuzWZbRYwt1EecxGfHwJMMtaOwo4B/iDMabbc1tr77fWVlhrK4YPH9731gZRxB2sU3mA6/g1AN/mLp7gDLdWa+AlSDyHbvyBi/kJ/wPA//JDvsyj/lqtjOmzWMJ9LTDa8/Eoug+7XAnMAbDWLgQKgWHxaGBWiLiD9R6+SRVPAHAOT/IWB7u1l16qgJfMV1MTvp/jJY7mcv4IwJf4E7fyI3+t7kDdI7GE+2JggjFmvDGmgNCE6dyImo+ASgBjzKGEwl3jLn1RW+sL+Kf5d8p5EzAcznI+9v6uVMBLJvME+woO4CRnM7BjeYk5RIR4ZaUmUPeQsTZyhCVKUWhp46+AXGCmtfY2Y8w0oM5aO9dZPfNbYAChIZvvWWuf6uk5KyoqbF1d3V5/A4FTVRUegwQoZTUNjCGHNhoppsg7wRTD/zuRtJOTA9aykUHOvR15HMgKVnKIv07BHpUxZom1tqLXuljCPREU7j3wLAsDGMwmtjGUPFpoosjdW2PwYE0ySWYZMgS2bqWZPEpoooN8RvAx6xnprysr097suxFruOsO1XQUsQZ+I/vQj520U8hAdrr7wG/dGvphEckExcXhu08HsoMO8hnI1u7BrkM34kLhnq7WrQv1zIE8YAcDyKeFFooZzFa3TgEvmaC4GJqbaQf6Ox2VIhrZGHmTktayx43CPZ1t2RIO+HxgKwPJZRdNDKIEz3CMAl7S2ZAh0NxMG1BCI20Uk08z2yjxb9+rYca4UrinO0/AF9PGFgaTQyuNDGYgm906Bbyko4ihmBb6U0ALOylWsCeYwj0TbNkS3qaghBYn4NvYwRAFvKQvz1DMQBppYQB5tLCNAQr2JFC4ZwrPNgUDaWYzg8mllR0MoZjt/knWnBytg5fUqakJdUacoZgB7KSF/uTRQiP9KaTDrS0qUrAniMI9U1RX+27BHkQTWxhELrtopsQZy3RYqxudJDWmTg3foNREPsW0sItiCmimiSL60enWFhVBU1OKGhp8CvdMUl0dCm7PEM0OBpDn7AXfnxaavH/wKuAlmaZODe8Vs5GBlLCTdvpRRCON0cbYFewJpXDPRJ4hmiLaaaKYIhppox8D2Mka9nNrL71Uu0lK4nmCfQUHMpzNdJJPCVvYGbkqRkMxSaFwz0QRQzT5QBMlDGITlnzGsJYlHO7Wz5ih/eAlcaqqwsFeyykcwrtADiNoYDtD/dvKaigmaRTumapriMZZJgmwlWGU8iGQQwVL+SvnuvXz5+tEJ4m/8vLwXki/ZTJn8BxgOJh3WB957IOGYpJK4Z7pPOvgAdYyniNZAhgu4lF+yE1u7fLloQOIReKhtDS8B9I3+TlX83sATuI5VhBxEufIkRqKSTKFexBEBPxSjuZc/g7ANH7MF0Nb7Yc0NITWH4vsqa6ljs5h1pXM425uAGAyM3mRz/nrR47UlgIpoHAPii1bfJuN/YMLuYlbAfg7F3EIb7q1zc0KeNkznqWOAGNYxTPOaWG/4FvM4ip/vYI9ZRTuQbJuXWirVMdt3MrfOBewrOAwStjiroVvbtbNTtI3nonTHRRSzE7WMB6wPMMp3MDd/vqyMgV7Cincg6a+3nei0wU8xkrGk0sbjQymH2281rWSputmJ62kkd6UloYnTp/nBAbSSDPF5NFCA/tyOgv89ZWV2rY3xRTuQVRb6zt0+0A+ooUCBrMJSy6TWMZP+S+3fv58TbRKdBHj6//NbZzGAiCHfWmgjSL2Z6P/c6ZM0QlKaUDhHlTTp/vWwucBWxjGcbwMwPf5BSfxvFvf0KBhGvGrqvKNrx9FHXdwIwBn8CSfRC51NCb0b2769GS2UnZD4R5kXWvhi4rCDy3kJP6HWwB4iZMpYau7ZUHXMI3Ww4tnGGYTAylmB0uZBFh+wfU8xTn++pEjobMz9G9O0oLCPRs0NfmWSk5jGnV8hhzaaWQQ/WlmDhe49cuXqxefrSKGYe7nqwxjC80MIJddvMMEbuAe/+do4jQtKdyzRcRSyUm8wS4K2J81QA7/wV+p4gm3XpOt2SdiGOZ4FnANvwMM43mPdgqZyCr/52jiNG0p3LPJunW+lTR5QANjuIrfADCfMylmB+sZ6n6OJluzg2cY5gNGUchOFnEiAN/hDlZxsL++a3xdE6dpS+GebWprfROtAL9lCos5ilxaaWYA+7GB7/Fjt0CTrcEVMQxzLXdzAB+xyznn9B0O4md83/85Gl/PCAr3bNQ10eq54amCZbRTSBlvAIafcRMjaKCZvFBB1zBNUZFCPghqaiAvLzwMs50ihrCR3/B1ACaxiFaKow/DaHw9Iyjcs1l9fbdefD1H8gBXAJ1sYH+K2cVtfMctaGnRWHym6xpb7wgdd/ddbmMQjWxlH6CDOXyROo73f46GYTKOsdam5AtXVFTYurq6lHxtiaK0NPynOUA7MJ4PWMtYAIaygVWMZxCeLVtzc+HBB/XneaaoqYHLLgv9FQasZygTeY9thA5VP4gVvM2hXX+rubQ/TFoxxiyx1lb0Vqeeu4REmWxdw3ju4jqgk82MYDA7uJzfuZ/T0aF18ZmivDz0/8oJ9guZw35sZBtDMXQwi2reixbsGobJWAp3cUWZbL2eGbSSx8EsBwx/4Kv0o4nnnZUUQGhdvDE6zi8dTZ0a+n/j7Ls+l7PJo4W/8SUAjmQJneQzmT/5P0/DMBlP4S5+XZOtnjXx+cAKynmKz5HPLlop5DReZCyr3AlXCO0YqJBPD12h7uziuJ0i9mcN5/EvOiigH028zLEsJcpf92VlWg0TAAp3iW7dOt/mYwBn8BytFHE5vwcsHzGeYlo5i7n+z50xI7QSQ6tqki8i1AFOZT6D2MknjAIs1/ErWujP8Sz2f25ubqi3rpuSAkHhLrs3fXqoF+8Ziwd4kCtpIZ8DeReAeXwBQzvf4Xa3qGs8Pj9fIZ8MNTWhexE8oX4N92Jo5wXnZKQyltFGLvc4pyb5VFZCe7t66wGicJfedY3F5+aGH+pHJyuZyDLKKWELkMMv+B557OKXfMP93Pb2UMhruCYxunrqnsnSW/kBObRyP1OAHIawgfcZQz2f6T5hmpensfWAUrhLbKqrQ0EdMeF6BG+znaHM4YsU0EIH+XybuyighTudG2LCNCYfP1GGX37E98ljF7cwDUsehTTxL85iMyM4gLX+z++aMG1rU289qKy1vb4BZwErgJXA96NcvxNY6ry9C2zt7TknTZpkJYOVlVkb6iv63u7hP20+LRY6LXTaXFrsjUyLWmunTEn1d5F5Kiu7vY5f506bQ2v4Ne9Hs32Qr0R/zSH0HJKxgDobS273WgC5wPvAAUABsAwo66H+G8DM3p5X4R4As2dbW1AQNUDuZKrNpzkcOIY2ez5/sR3RwsYYBX1Ppkzp9pq1ga3icQvtzkOdth877e+4bPehXlaW6u9E4iDWcI9lWOYYYKW1dpW1thV4GDivh/pLIHLRrARSdTXs2tVtPB7gW0ynlSIe4isU0Ygll0e5iFw6OJB3WMVot9had8hG2xq4qqq6Db28zQTG8j75dFLL2UAOA9jK3ziXFvrzNf7Q/XkKC7UKJgvFEu6lwBrPx2udx7oxxowFxgPP7H3TJGN4x+ON8V26jIdpYiCvMol9WQcYVjGRA1lNEY3cyK3+55o/P/QcxkBJSXattKmpgQED3O/f2YIX4Jv8gn40UcYKPuIAAEpZTT0T2cEQLuCx7s/XNVna3Kxx9SwUS7ibKI/tbkOai4FHrLUdUZ/ImKuNMXXGmLoNGzbE2kbJFNXVoZtfItbHAxzN63zCKJoo4GSewdBBC8Xczv9i6GA0H/AcJ/g/qbHRXWkT5B59Vw/90kth587ww49zBvuxFkMHd3MDrRRi6OBM/kUrOaxlHGW81/35ukJdk6VZLZZwXwvev6EZBTTspvZiehiSsdbeb62tsNZWDB8+PPZWSmbpWh8/ezYUFPguFdHOC1TSST6/YzKD2AwY1jKO01mAoY1DeYMlHO5/Tm+PPicns1fcdK10idJDf4HjOIi3MbTz78xjPaWAYQgb+DMX0kk+T/L5rlNv/bqGXxTqAjFNqOYBqwgNt3RNqJZHqZsIfIiz02Rvb5pQzSI9TLxasK1gv8RsW0CTMwFrnf+22zF8YB/iP3Y/Sdj1NmBA6Oukm9mzre3fv8e238dX7f585EyOut9/P3bayTxgWzP1e5eEIF6rZULPxTmElji+D9zsPDYNONdTcwtweyzPZxXu2SvKUj7v23YK7Zk85llpY8NhV8QOexpP2+Uc2HvYpyL0YghyC/Z1yu2JPG8Laez2PfajyZ7HI7aJvN6/Py1pzEpxDfdEvCncs1wMQdgK9gp+YweyyUJHRK++wxaywx5BnZ1JdWxhH+2tsDC2XwBRliPG8tYJ9j6utIfxuifM/d/HYDbaKfyq9x66eunGGwS7AAAHiElEQVRiYw93HdYhqTd1qm+53+4spILruZu3OIIWigjN9XfN91ugkyKaGUkDJ/M8U7mXo1mWwIb7vUwF9/INFnE8DezvtDEnoo2WQpqYxBLu5ut8lrdie/IpU0JzGZL1Yj2sQ+Eu6SXGoAeo52Bu5v9YwMlsYQid4Z1TvAu8QoFq6KQfrRSxk4FsZ18+YX8+ZgxrGMsqSmlgPz5hHzbTjxaaKWIj+/Ap+/IxI1nFAaxlNOsYyafsxw5K2El/2sjHhgM88utCDu0MYxMn8zy38z0O4qPYXgdj4NprFejSjcJdMl8fgr7LGvbl5/w3tXyOtYyhiWLayccN3mgre/dU189OJ3m0U8xOxrGaSp7me/yU/djct6dToEsMFO4SPFVVvmWDfdUGLORY5lPFO5SxmjFsZxCN9GcXRXSSw07600EuuXRQzE5y6KSYZopoYihbGM1qDuFtzuApjqdu73feq6zUjozSJ7GGe7cdQEXSVmQI9rFnnw+cwiucwivxbVesCgvhgQe0Bl2SQlv+SubqulnK+xbl7tiU6LqhyNs2bQMgSaRwl2CJFviRb7NnQ//+e/41BgzoHtyRbwpySTENy0j2qa5W8ErgqecuIhJACncRkQBSuIuIBJDCXUQkgBTuIiIBpHAXEQkghbuISAAp3EVEAkjhLiISQAp3EZEAUriLiASQwl1EJIAU7iIiAaRwFxEJIIW7iEgAKdxFRAJI4S4iEkAKdxGRAFK4i4gEkMJdRCSAFO4iIgGkcBcRCSCFu4hIACncRUQCSOEuIhJAMYW7MeYsY8wKY8xKY8z3d1PzZWPMcmNMvTHmj/FtpoiI9EVebwXGmFzgXuAMYC2w2Bgz11q73FMzAbgRONFau8UYMyJRDRYRkd7F0nM/BlhprV1lrW0FHgbOi6j5T+Bea+0WAGvtp/FtpoiI9EUs4V4KrPF8vNZ5zOtg4GBjzEvGmEXGmLPi1UAREem7XodlABPlMRvleSYApwGjgBeNMYdZa7f6nsiYq4GrAcaMGdPnxoqISGxi6bmvBUZ7Ph4FNESp+Ye1ts1a+wGwglDY+1hr77fWVlhrK4YPH76nbRYRkV7EEu6LgQnGmPHGmALgYmBuRM2jwOkAxphhhIZpVsWzoSIiErtew91a2w58HZgHvA3MsdbWG2OmGWPOdcrmAZuMMcuBZ4HvWms3JarRIiLSM2Nt5PB5clRUVNi6urqUfG0RkUxljFlira3orU53qIqIBJDCXUQkgFI2LGOM2QCsjtPTDQM2xum5EilT2gmZ09ZMaSdkTlvVzviLZ1vHWmt7XW6YsnCPJ2NMXSxjUKmWKe2EzGlrprQTMqetamf8paKtGpYREQkghbuISAAFJdzvT3UDYpQp7YTMaWumtBMyp61qZ/wlva2BGHMXERG/oPTcRUTEI6PCvbcToYwx/Ywxf3auv2KMGZf8VsbUzlOMMa8ZY9qNMReloo2etvTW1hucE7beMMbMN8aMTdN2XmuMedMYs9QYs8AYU5aO7fTUXWSMscaYlK32iOE1vcIYs8F5TZcaY65Kx3Y6NWlxElwMr+mdntfzXWPM1mjPExfW2ox4A3KB94EDgAJgGVAWUTMVuM95/2Lgz2naznHAEcBDwEVp/pqeDhQ7709J49d0oOf9c4En07GdTl0J8AKwCKhI4//3VwD3pKJ9fWznBOB1YIjz8Yh0bWtE/TeAmYlqTyb13GM5Eeo84EHn/UeASmNMtP3oE6nXdlprP7TWvgF0JrltkWJp67PW2ibnw0WEtnxOtljaud3zYX+6nzmQDLH8GwX4EXAH0JLMxkWIta2plkknwfX1Nb0E+FOiGpNJ4R7LiVDhGhvazXIbsE9SWhelDY5o7UwXfW3rlcATCW1RdDG10xhznTHmfULBeX2S2ubVazuNMUcBo621jyWzYVHE+v/+QmdI7hFjzOgo1xMtk06Ci/nnyRneHA88k6jGZFK4x3IiVCw1iZYObYhVzG01xlwKVAA/S2iLooupndbae621BwL/Dfwg4a3qrsd2GmNygDuBbyetRbsXy2v6T2CctfYIoBb3r+Jk6utJcJcADxhjBie4XdH05Wf/YuARa21HohqTSeEe64lQowGMMXnAIGBzUloXpQ2OaO1MFzG11RhTBdwMnGut3ZWktnn19TV9GDg/oS2Krrd2lgCHAc8ZYz4EjgPmpmhStdfX1Fq7yfP/+7fApCS1zStuJ8ElQV/+nV5MAodkgIyaUM0jdLrTeNzJivKImuvwT6jOScd2empnkdoJ1Vhe06MITRJNSPN2TvC8/wWgLh3bGVH/HKmbUI3lNd3f8/4FwKI0bedZwIPO+8MIDY3sk45tdeomAh/i3GeUsPak4h/WXrx45wDvOmFzs/PYNEI9SoBC4C/ASuBV4IA0befRhH7L7wQ2AfVp/JrWAuuBpc7b3DRt511AvdPGZ3sK1VS2M6I2ZeEe42v6E+c1Xea8poekaTsN8EtgOfAmcHG6vqbOx7cAtye6LbpDVUQkgDJpzF1ERGKkcBcRCSCFu4hIACncRUQCSOEuIhJACncRkQBSuIuIBJDCXUQkgP4fR6bVmpXv/3sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0de01852e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 1000\n",
    "w = [1.0]*n\n",
    "l = {(i,i+1): 1.0/n for i in range(n-1)}\n",
    "f = {0: (0.0,1.0), n-1: (0.7,1.0)}\n",
    "g = 9.81\n",
    "\n",
    "x,c = stringModel(w, l, f, g)\n",
    "if x is not None:\n",
    "    display(x, c, list(l.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also have more suspension points and more complicated shapes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl4VNX9x/H3l7AJgopERCAJKHGvIJFWK66oiFXchUYRraaKouIKv4iigtalblVUROpChIJVxH1BEGsFCagoWhaVAIIQQEGMbMn5/XEmEkNIJmRm7iyf1/PMMzN3bu58Lxk+uXPuueeYcw4REUku9YIuQEREIk/hLiKShBTuIiJJSOEuIpKEFO4iIklI4S4ikoQU7iIiSUjhLiKShBTuIiJJqH5Qb9yyZUuXlZUV1NuLiCSkWbNmrXLOpde0XmDhnpWVRWFhYVBvLyKSkMysKJz11CwjIpKEFO4iIklI4S4ikoQU7iIiSUjhLiKShBTuIiJJSOEuIpKEFO4iIklI4S4ikoRqDHczG21mK83si+28vp+ZfWRmG83s+siXKCICFBRAVhbUq+fvCwqCriiuhXPk/jTQo5rX1wBXAfdFoiARkW0UFEBeHhQVgXP+Pi9PAV+NGsPdOTcNH+Dbe32lc24msDmShYmI/Co/H0pK6M+jHMkHrKcplJT45VKlmLa5m1memRWaWWFxcXEs31pEEtnixQC8wcl8yJF8xf6/WS7bimm4O+dGOudynHM56ek1jlgpIuJlZLCaFiyhHS0p5jAKf10uVQtsyF8RkbANH87f+hVTuiWN9nzrlzVpAsOHB1tXHFO4i0jc++6YXB6xLbSqt4omZSWQmemDPTc36NLiVo3hbmZjgWOAlma2FLgVaADgnHvczPYECoHmQJmZXQMc4JxbF7WqRSSlDBsGpdQn67B0aHwMTF0UdElxr8Zwd871qeH174G2EatIRKSCr7+GUaPgr3+FL6q82kaqoitURSSuDR0KDRqo12NtKdxFJG598YW/Tumqq6B166CrSSwKdxGJW0OGQLNmcOONQVeSeBTuIhKXZsyAiRPhhhugRYugq0k8CncRiUv5+ZCeDldfHXQliUn93EUk7kye7G8PPOCbZaT2dOQuInHFOX/U3q4dXHZZ0NUkLh25i0hceeUV394+ahQ0bhx0NYlLR+4iEjfKyvxRe8eOcOGFQVeT2HTkLiJxY9w437d93Dior3SqEx25i0hc2LwZbrkFDjkEzjkn6GoSn/42ikhcGD3ajyPz6qt+mlSpG/0TikjgfvkFbr8djjgCevYMuprkoCN3EQnciBGwbBmMHQtmQVeTHHTkLiKBWrcO7roLTjoJjjoq6GqSh8JdRIJRUABZWTywy1BWr4bhR74RdEVJReEuIrFXUAB5eRQXrefvXMtZvECXu872yyUiFO4iEnv5+VBSwp94jfXszK0MhZISzcgRQQp3EYm9xYsppAszOYx9mcfBzP11+XYVFMD06fD+VMjK0lF+DRTuIhJzm9rtzcWMpjXL+YjDt76QkVH1D4Sacdi4wT8vKvLPFfDbVWO4m9loM1tpZlVOTWvew2a20MzmmNmhkS9TRJLJ3w57gc/5HY9xObuy1i9s0gSGD6/6B0LNOAvZh4Xs45epGada4Ry5Pw30qOb1k4GOoVse8Fjdy0pioR4C1Kunr5aSkr74AoZNOoQ+hy/itMw5vmN7ZiaMHAm5uVX/0Paaa6prxklxNV7E5JybZmZZ1azSC3jWOeeA6Wa2q5m1ds4tj1CNyaP8q2VJCQ6w8q+WsP0PtUgS2bIFLr4YdtkFHno5C9IXhfeDGRlQVMQ+LNx2uVQpEm3ubYAlFZ4vDS2TyvLzWV3SmJ68xguc7Zfpq6WkkAcfhJkz4ZFH/BR6YRs+3DfbVFRdM45EJNyruljYVbmiWZ6ZFZpZYXFxcQTeOsEsXsyu/Mh3tOFG7mEDjX5dLpLs5s+HIUOgVy8499xa/nBurm+2aRSavaOmZhyJSLgvBdpVeN4WWFbVis65kc65HOdcTnqt/mwniYwM0ijjAQayiPY8yDW/LhdJZmVlcMkl0KiRH0dmh8aPyc2FP/wBjj4GFi1SsNcgEuE+Cegb6jXzB2Ct2tu3I/TV8jim0IuJDCef73dqr6+WkvQefxw++MBPeL3XXkFXkxpqPKFqZmOBY4CWZrYUuBVoAOCcexx4HegJLARKgIuiVWzCKz/SyM/n3qIbOZAvuLnr24zK3SfYukSiqKgIbroJTjwR+vULuprUYb6TS+zl5OS4wsLCQN47Xlx3nT+SmTULOncOuhqRyHMOevSADz/0XSCzsuq2vWOO8fdTp9axsARmZrOcczk1racrVAM0ZAi0aAEDB/r/BCLJ5umn4e234e676x7sUjsK9wDtuivccQe8/z5MnBh0NSKRtXw5XHstdOsGl18edDWpR+EesEsvhQMPhOuvh40bg65GJDKc84G+YQOMGqU5UYOgf/KA1a8P998P33wDDz8cdDUikTF+PLz8sp8XNTs76GpSk8I9Dpx4IpxyCgwbBitXBl2NSN2sWgUDBkBOjj+fJMFQuMeJ++7zIxHcckvQlYjUzdVXw48/wujR/pupBEPhHif22w/694cnn4Q5c4KuRmTHvPIKPP+8Hy7p4IODria1KdzjyK23+tHyrr1WXSMl8fz4I1x2mQ/1wYODrkYU7nGkRQu47TaYPBlefTXoakTCUGF+guvbjeP75WWMHg0NGwZdmCjc48xll/kmmuuug02bgq5GpBrl8xMUFfGuO46n1vfm+rQHyZmnCWjigcI9zjRoAH//OyxYAI8+GnQ1ItUITX23gnTOZTwdmc/QLfmanyBOKNzj0Mknw0kn+SaaVauCrkZkOxYvppR69OR1fmA3ruFBdmKD5ieIEwr3OGTmj97Xr4ehQ4OuRmQ7MjK4lvuZTQ75DKN/+fTJmp8gLijc49SBB/r298cfh7lzg65GZFsPdnuBh7magdzPMEIXaGjqu7ihcI9jQ4fCzjura6TEn5degmsLcjjzsMXcl/EP/3VTU9/FFYV7HGvZ0vd9f/tteGPPfn70paws30tBJCAzZvj87toVnpuaQb2ib/08epr6Lq4o3OPcFS3G0tEWcN3Km9js0vy0Nnl5CngJxDffwKmnQuvWMGmSb4WJmYICmD4d3p+qg5wwKNzjXMNbB3Ofu5b/sT8X85RfWFKi7mYSc2vWQM+esGULvP467LFHDN+8vE/9xg3+uQ5yaqRwj3eLF3MKr7EbaxjLn1nI3r8uF4mVjRvhjDPg22/9UL777hvjAkJ96n9DBznVUrjHu4wM0nBMoxs7s57zGcMW0tTdTGKmrAwuugimTfPT5nXrFkARoYOZTnxKJz7dZrlsS+Ee74YPhyZNOIgveZzLmMEfGN5gqLqbScwMGQJjx8Kdd0KfPgEVETqYeZCBPMjAbZbLtsIKdzPrYWbzzGyhmQ2q4vVMM5tsZnPMbKqZtY18qSkqN9d3L8vMpLeNJ7fpS9xR+n9M31u9EiT6Ro3yoX7ppTBom//5MRQ6yPkN9amvnnOu2huQBnwNdAAaAp8BB1RaZwJwYejxccBzNW23S5cuTmrvxx+dy8hwbu+9nfvpp6CrkcCNGeNcZqZzZv5+zJiIbfqtt5xLS3PupJOc27QpYpvdcVHc10QCFLoa8tU5F9aRe1dgoXPuG+fcJmAc0KvSOgcAk0OPp1TxukTILrvAc8/5LmmawizFhXqQ/FC0lnWuaUR7kHz2GZx9tr9Sevx4P6Bd4HJzfV969akPSzjh3gZYUuH50tCyij4Dzgo9PgNoZma71708qcpRR8FNN/mvzC+9FHQ1EphQD5IBPEwrVnILt7GypGmde5B8952f07d5c3jtNX8viSeccLcqllW+GP564Ggz+wQ4GvgO2LLNhszyzKzQzAqLi4trXaxsddttcOihvi10+fKgq5FAhHqKrGQPGrOBYdxMJkX0L7qJhQt3bJPr1vlgX7fOB3tbnT1LWOGE+1KgXYXnbYFlFVdwzi1zzp3pnOsM5IeWra28IefcSOdcjnMuJz09vQ5lS8OG/tt3SYnvpqaxZ1JQqKfIOnahKx/zFftzAc/xFBez777w3nu129zmzXDuufDFFzBhAhxySBRqlpgJJ9xnAh3NrL2ZNQR6A5MqrmBmLc2sfFuDgdGRLVOqst9+cN998NZbmtgjJYV6kKyhBS1Yw77MZ2STgRQ98ipDh8KRR/rVXnjBX1Fa3QGAc3DFFf6z9Nhjfj4BSXDhnHUFegLz8b1m8kPLbgdOCz0+G1gQWmcU0Kimbaq3TGSUlTl38snONW7s3Ny5QVcjMTdmjGtRb43rz6Pb7UFy+OHOgXMHHeTcM884t3Hj1p8t731y165/c+Dc4MGxLF52BGH2lgkr3KNxU7hHzvLlzrVs6VynThX+40pKKC31PQOHDNn+Ops2Offssz7cwbk2bZx74ar3nWvSxDlwYznPgXN90v7lSp9Nze6FiSTccNcVqklgzz3hqafg00/91YSSOtau9U0qLVpsf50GDeCCC2DOHN88k50NTcaNhpISHqE/fXmWbkzjn6UXUG+IxmpJFgr3JHHaab6L8733wtSpQVcjsbJ6tb+vLtzLmfn5ed97D04ufpZnuIABPEIDNjOR02nEJo3VkkQU7knk/vthn32gb1/48cegq5FYWLPG3+9ei6tKNm6EK5o+TT+epROf8i7H04If/IsaqyVpKNyTSNOmMGYMLFvmez5I8isP93CO3MEfmB91FIxY35fr6z/Ix3TlcGb4FzVWS1JRuCeZrl391HzPP+9vktxq0yzz1lv+wrevvoJ//xvufTqdBpltNP9pklK4J6HBg+GII6B/fzWhJrtwmmXKyuD22317e+vWUFgIZ56JxmpJcgr3JFS/vh9crLTUt7+XlgZdkURLebjvumvVr69e7YcTuPVWOP98PwVpdnbs6pPgKNyTVIcO8I9/wPvvw99b3gX16mlS4SS0erUfKbR+/W1fmznTN8O89x48/jg884w/LyOpQeGexC6sX8BZaS9x84/XMdsdokmFk9CaNds2yTgHTzyxdfiB//wH/vpX37QuqUPhnsTs5nyeKL2EZqyjG/9hBemaVDjJrFnz25OpJSVw4YVw2WVw3HEwezYcdlhw9UlwFO7JbPFidmcN13I/JTThdF5mA410ljWJrF69NdwXLIA//MF3h73tNj9kb236v0tyUbgns9AFKfncxbP0ZTqH04exbGnXPuDCJFLKm2VefBFycvxEG2+8Abfc4k+zSOrSrz+ZVZhU+ALG8BBXMZEzyMt6W+O/J4nVq32/9bPO8kNAf/KJhusVT+GezHJz/YUpmZlgxlWZk7jljDn8c9re3HRT0MVJXS1dCj/84AeM698fpk3T6AGyVRUdqCSp5Ob+5uKUoQ5WXekHGGvZEm68McDaZIdNmwbnnOMfn3++JmuRbenIPcWY+f7vvXtvnWRbEodzfvat447b2mf9xBODrUnik8I9BdWr5y9oOekk3//5xReDrkjCsXatb1u/4QY4/fStf5jDHTRMUovCPUU1bOgHj+raFfr0qf1kyhJbc+b43jCTJvmhnSdM8EP3gro7StUU7imsaVPfF7pjR+jVyw8oJfHnued8//Wff4YpU2DgQLDnC1hzwdUAtDjrWF11LNtQuKe4Fi38ULAtW/pRA//3v6ArknIbNvgrTfv2hd//3l9t2q0bPsjz8li92vdnbbHscw0rIdtQuAtt2sDbb/u2+BNPhCVLgq5IFi3yQf7EE/7E9zvv+LlyAT98REkJs+kMwK78oGElZBthhbuZ9TCzeWa20MwGVfF6hplNMbNPzGyOmfWMfKkSTR07wptv+pN2J54Iq1YFXVHqeuMN6NIF5s+Hl16Cv/2t0qiPoeEjPuFQjDLqU/ab5SIQRribWRrwKHAycADQx8wOqLTazcB451xnoDcwItKFSvR17gyvvOKPGnv2hJ9+Crqi1FJa6sddP+UUaNsWZs3yvWK2EbpSqSPz2YeF2ywXgfCO3LsCC51z3zjnNgHjgF6V1nFA89DjXYBlkStRYumoo+Bf//Ltu2eeubVHhkTXqlX+D+rtt/s29o8+8pOdVyk0rMRKWrFX+X81zX8qlYQT7m2Aiq2wS0PLKhoKnG9mS4HXgQFVbcjM8sys0MwKi4uLd6BciYXTToPRo+Hdd/3Vj5rJKbpmzPCTakyd6keL+Oc/fx0SqGqhYSVW1G9DK1Zo/lOpUjjhXtUQ/5WHneoDPO2cawv0BJ4zs2227Zwb6ZzLcc7lpKen175aiZm+fX1/6hde8OOWaKCxyHMORozwJ07T0uC//4VLLw1zUo3cXFY03ZtWA87T/KdSpXDGllkKtKvwvC3bNrv8BegB4Jz7yMwaAy2BlZEoUoIxcKBvLrjzTt9VUt/6I+fnn/3VwQUFvjnmuedqd6Xphg3+5PevPWhEKgnnyH0m0NHM2ptZQ/wJ00mV1lkMHA9gZvsDjQG1uySBYcN8F+o774T7W9yhuVgjYN4832/9+efhjjv8SezaDiGwYoW/b9Uq8vVJcqjxyN05t8XMrgTeAtKA0c65uWZ2O1DonJsEXAc8aWYD8U02/ZzTF/lkYAYjjnyeNU815rofhrCFddxYdJ9PfFBzQC298AJcdBE0buwvHjvhhB3bjsJdamJBZXBOTo4r1PXuiSEri41Fy8mhkC84iP6M4B8MoF5mhm/vlRpt3uwvRnrgAX/UPmECtGtX889tzyuv+BPfH3+sOVJTjZnNcs7l1LSerlCVmi1eTCM28QY9OIC5jOAKTmMSq4vWB11ZQli2DI491gf7gAF+LPa6BDtsPXJXm7tsj8Jdaha6OKYty/iCg3mEK3iHE+iUNof//jfg2uLc1Kn+4rBPPvFt7A8/7EfkrKvvv/f3e+xR921JclK4S80qzMVqwBWM4L+Nj6Ph7s046ii45x4oKwu2xHjjHNx9Nxx/POy2m28+6dMncttfsQJ23RUaNYrcNiW5KNylZpXmYiUzky6j+jN7fjPOPNO3JZ96qsajKffjj3DGGTBoEJx9NsycCQceGNn3WLFCJ1Olegp3CU9urj95Wlb260Uzu+zihyp49FF/NWvnzvDhh0EXGqzPPvOTarz2Gjz4IIwbB82aRf59VqxQe7tUT+EudWLmr2D96CPfRHD00b45IhWbaZ5+2k+q8csvvq396qvDvNp0B3z/vY7cpXoKd4mIQw/1oxieeaZvjvjTn1KnmWbDBt/t/6KL4PDD/cnTP/4xuu+pZhmpicJdIqa8mWbECJg8GTp1gv/8J+iqouvbb32QP/kkDB7sJz2Jdg+W8qEHFO5SHYW7RJQZXH65b6Zp3BiOOcZPNpGMzTSvveYn1fj6az9x9Z13VppUI0pWhkZsUpu7VEfhLlFx6KF+TPizzvJHtKecAskyynNpKQwZ4pueMjN9c9Spp8bu/cv7uOvIXaqjcJeoad7c9xYZMQLee8/3pkn0ZpriYujRww+odvHFfpjevfeObQ0aV0bCoXCXqCpvppk+HXbaKbGbaaZP999IPvgARo2Cp57y+xRrCncJh8JdYqJzZ998cfbZiddM4xz84x9+CsIGDfz5hL/8Jbh6FO4SDoW7xEzz5jB2LDz2GEyZ4nvTfDDkbT8+fJyOE79+Pfz5z3DVVb45ZtYs/4cqSN9/r6EHpGYKd4kpM7jsMt/E0aR0HccOO46bii6j1AFFRb7DeJwE/FdfQdeuMH687wkzcaIfJyZo6uMu4VC4SyA6dYJZDY/gDF7kHgaRzioepT/rStIgPz/o8hg/3o+TvmqV77s+eLD/chG4ggJWTJpBq3nvx+U3HYkf8fBxlRTVfOmXjOF8+vIMrVnOlTzKXizj8qJBzJkTTE2bNsE118B558HvfuevNj3++GBq2UZBAeTlsWLjLuzJ93H3TUfii8JdgpORQSM28wz9mMtBfMxhnMMEnrZ+HHIIdOvm2+g3boxNOd995yfVeOghPy7M1KnQpk1s3jss+flQUsICOvIdocJKSuLim47EH4W7BKfCOPEAh1HIP5tcydIRr3Dfff7E4Z//7Gct+r//8weq0TJ5sj9R+tlnfgiFBx+MzKQaEbV4MWtpRhlptGPxb5aLVKZwl+BUMU48I0ey+2XncN11MG+en0T6iCP8SJMdOvh5Q998M3L95MvK4K674MQToWVLP/b6uedGZtsRl5HBEvysWKfz8m+Wi1QWVribWQ8zm2dmC81sUBWvP2Bmn4Zu883sx8iXKkmpinHiy9Wr50N34kQ/QNfgwTBjBpx8MnTsCPfdB6tX7/hb//AD9OrlvxWce66fLWn//eu8R9EzfDiLGu0LQBaL/LImTfw3IJHKnHPV3oA04GugA9AQ+Aw4oJr1BwCja9puly5dnEhtbdzo3NixznXr5hw416iRc337Ojd9unNlZeFvZ/Zs59q3d65BA+cefrh2PxukRy782IFzy9nTucxM58aMCbokiTGg0NWQr865sI7cuwILnXPfOOc2AeOAXtWs3wcYu8N/bUSq0bAh9O4N06bB55/7K0VffNFPkpGT44cEKCmpfhtPPeXHXd+0Cd5/HwYMiN6kGpG2KP0wGjeGVmXLt/mmI1JROOHeBlhS4fnS0LJtmFkm0B54r+6liVTvoIP8FH/LlvnByTZtgksu8T1cBg70bfYUFPx6BewvGfvyl2MWcsklvifOJ5/4kE8kixZtPUUhUp1wwr2qj5Hbzrq9gRecc6VVbsgsz8wKzaywOFEGFpG416yZH5xszhx/RN+jhw/9/faDE/q25sWiQ5nv9uaIJeMY/f4+3Nzrc958E9LTg6689hYt8n+rRGoSTrgvBdpVeN4WWLaddXtTTZOMc26kcy7HOZeTnoj/sySumW3tG79kCQzf9R7mle3DWbzIvsxnAR15lVO449NTSUsLutodU1Tkj9xFahJOuM8EOppZezNriA/wSZVXMrN9gd2AjyJbokjttWoF/7d2EN/QgVu5BTBuZSin8HrC9gv/+Wc/kqaO3CUcNYa7c24LcCXwFvAVMN45N9fMbjez0yqs2gcYFzqbKxK8jAzqU8r13A/AFhr8ujwRlV/EpXCXcIQ146Nz7nXg9UrLbqn0fGjkyhKJgOHDIS+PnUt+Zi++Yz7ZCd0vfNEif69wl3DoClVJXhWugM1mPvMbHeyfJ2j3wfJwV5u7hEPhLsktdAVsdt6xzNs5J2GDHXyzTMOGsOeeQVciiUDhLikhO9sPVVCX4QqCVt7HPS7GlZe4p4+JpITsbH+/YEGwddSF+rhLbSjcJSWUh/v8+cHWURcKd6kNhbukhA4dIC0tccP9l19g5UqdTJXwKdwlJTRo4AM+UcNdfdylthTukjKysxM33NXHXWpL4S4pozzcIzWLUywp3KW2FO6SMrKzfdv1d98FXUntFRX5pqXWrYOuRBKFwl1SRiL3mFm0yA+Joz7uEi59VCRlJHq4q0lGakPhLimjTRs/blhChXtoJqlF05eTNWOcfy4ShrBGhRRJBmYJ1mOmoADy8thQUsr3tCZz/VzI88MXJ/IYORIbOnKXlJJQ4Z6fDyUlzOFgAFqw2s/+nZ8fcGGSCBTuklKys+Hbb/1k2nEvNGPU25wIgCufzjhBZ5KS2FK4S0rJzobSUvjmm6ArCUNoxqjNoRmkzmHCb5aLVEfhLikloXrMDB8OTZowh0PIZh6tKE7omaQktnRCVVJKQoV76KTp7AsP44jSD/yoYcOH62SqhEXhLillt90gPT1Bwh1YdVIui0thwL194Po+QZcjCUTNMpJaCgrIXjeT+U9O9VcFxXm/8U8+8feHHhpsHZJ4wgp3M+thZvPMbKGZDdrOOuea2ZdmNtfMno9smSIREOo3nr3xc+aT7QdsycuL64CfPdvfd+4cbB2SeGoMdzNLAx4FTgYOAPqY2QGV1ukIDAb+6Jw7ELgmCrWK1E2o33gLVrOcvVhFi7jvNz57NrRv75uTRGojnCP3rsBC59w3zrlNwDigV6V1LgUedc79AOCcWxnZMkUiINQ/PA0/5u/L5R/jOO43Pnu2mmRkx4QT7m2AJRWeLw0tqygbyDazD81supn1iFSBIhET6h9+DQ9Qj1IWk/mb5fFm7VpYuFDhLjsmnHC3Kpa5Ss/rAx2BY4A+wCgz23WbDZnlmVmhmRUWFxfXtlaRugn1G2/NCg5jJu/SPa77jetkqtRFOOG+FGhX4XlbYFkV67zsnNvsnPsWmIcP+99wzo10zuU453LS09N3tGaRHZObCyNHQmYm3ZnMDH7PugdHx22/cZ1MlboIJ9xnAh3NrL2ZNQR6A5MqrTMROBbAzFrim2kS4QJvSTW5ubBoEd2n5FNKfd7f87ygK9qu2bP9MMWtWgVdiSSiGsPdObcFuBJ4C/gKGO+cm2tmt5vZaaHV3gJWm9mXwBTgBufc6mgVLVJXhx8OO+0E77wTdCXbN3s2dOkSdBWSqMK6QtU59zrweqVlt1R47IBrQzeRuNeoERx1FLz7btCVVO3nn+F//4Pz4veLhcQ5XaEqKat7d/jqq/icMPuzz8A5nUyVHadwl5R1wgn+fvLkYOuoSvnJVIW77CiFu6Ssgw/2g4jFY9PM7Nmwxx6w115BVyKJSuEuKatePTj+eB/urvKVGwGbNcsftVtVV5mIhEHhLimte3dYvty3vceLDRtg7lw1yUjdKNwlpXXv7u/jqUvk55/7qQAV7lIXCndJaZmZsM8+8dXurpOpEgkKd0l53bvD1KmweXPQlXizZ/shfrOygq5EEpnCXVLeCSfA+vXw8cdBV+KVD/Ork6lSFwp3SXnHHuuDNPCmmYICNmfuw5zCjRw68/G4niFK4p/CXVLebrtBTk7A4R6aAnDW4t3ZRCMOXTc17qcAlPimcBfBt7tPnw4//RRQAaEpAB9hAADNWBf3UwBKfFO4i+DDfcsWmDYtoAIWL2YLaUzlaDqwkB68+etykR2hcBcBjjgCGjcOsL97RgYTOIfvaMf9XEda+WRncToFoMQ/hbsIPti7dQuu3d0NG87dNpj9+IpTecUvjOMpACX+KdxFQk44wV/2v3x57N/77fRcPnO/44bdR1PP8FdXjRzKkRTIAAAIsklEQVQZt1MASvxTuIuEdN/k56OZvNcF/gqiGPZUuftuPwJk7nf3QlkZLFqkYJc6UbiLABQUcMjwc9mdVbzL8VBUFLOuiDNnwpQpMHCgnyFKJBIU7iIA+fnU++VnOvEJ4zmXn2gas66I99wDu+zi/5aIRIrCXQR+7XK4F8v5hZ34PR/zNR2i3hVxwQL497+hf39o3jyqbyUpRuEuAr92OXyGC3mHE1hBK3Io5K3086P6tvfdBw0bwtVXR/VtJAWFFe5m1sPM5pnZQjMbVMXr/cys2Mw+Dd0uiXypIlE0fDg0aYIB3ZlMITlk2FJ6Fj/N3XdHZ6am77+HZ56Bfv2gVavIb19SW43hbmZpwKPAycABQB8zO6CKVf/lnOsUuo2KcJ0i0ZWb67seZmaCGe0zHf99ci7nnFuPQYOgd2/4+efIvuVDD/lhhq+/PrLbFQGoH8Y6XYGFzrlvAMxsHNAL+DKahYnEXG7ub7ofNgXGXuyH3x082E/FN3EidOhQ97datw4eewzOOstPFiISaeE0y7QBllR4vjS0rLKzzGyOmb1gZu2q2pCZ5ZlZoZkVFhcX70C5IrFlBjfeCK+/DkuW+NEjIzFEwRNPwNq1cNNNdd+WSFXCCfeqpgyo3AL5CpDlnPsd8C7wTFUbcs6NdM7lOOdy0tPTa1epSIBOOgkKC6FNG+jRw58I3dF2+I0b4YEH4PjjoUuXyNYpUi6ccF8KVDwSbwssq7iCc261c25j6OmTgD6yknT23hs++gjOPBNuuMG34JSU1H47Y8b4IQ501C7RFE64zwQ6mll7M2sI9AYmVVzBzFpXeHoa8FXkShSJHzvvDOPHw513wrhx8Mc/+pECwlVWBvfeC507+2GGRaKlxnB3zm0BrgTewof2eOfcXDO73cxOC612lZnNNbPPgKuAftEqWCRoZv4E62uv+WDPyYHJk8P72Zdfhnnz/FG75kiVaDIXjQ68YcjJyXGFhYWBvLdIpCxcCKef7nvS3HuvHx9me6HtHBx+OKxcCfPnQ/1w+qqJVGJms5xzOTWtpytURepgn318O/zpp8N118EFF2y/HX7aNJgxw/drV7BLtCncReqoWTOYMAGGDYPnn4cjj/SDSlZ2zz2Qng4XXRT7GiX1KNxFIqBePT+A5CuvwNdf+3b4KVO2vv75576v/FVXwU47BVenpA6Fu0gEnXIKfPwxtGzpZ3Z66IJCXGYW9/zuOZraz/RPnxB0iZIiFO4iEbbvvr5t/U+dlnDNmBzOWvx3xtKHPPcELa7tF9MZniR1KdxFoqB5c3ix+CiGcisvcRal1ONyRsRsAhARnbMXiZJ6S4q4ldtZT1NmkUNHvvYvRHkCEBFQuItET0YGFBVxLzdtu1wkytQsIxItoQlAfqNJE79cJMoU7iLRUmkCEDIz/fMKY8aLRIuaZUSiqdIEICKxoiN3EZEkpHAXEUlCCncRkSSkcBcRSUIKdxGRJKRwFxFJQgp3EZEkpHAXEUlCCncRkSQU2ATZZlYMVJyMrCWwKpBigpFK+6t9TV6ptL/xsq+Zzrn0mlYKLNwrM7PCcGb0ThaptL/a1+SVSvubaPuqZhkRkSSkcBcRSULxFO4jgy4gxlJpf7WvySuV9jeh9jVu2txFRCRy4unIXUREIiSwcDezFmb2jpktCN3vVsU6nczsIzOba2ZzzOy8IGqtCzPrYWbzzGyhmQ2q4vVGZvav0OszzCwr9lVGRhj7eq2ZfRn6XU42s8wg6oyEmva1wnpnm5kzs4TpZVFZOPtqZueGfrdzzez5WNcYSWF8jjPMbIqZfRL6LPcMos4aOecCuQH3AINCjwcBd1exTjbQMfR4L2A5sGtQNe/APqYBXwMdgIbAZ8ABldbpDzweetwb+FfQdUdxX48FmoQeX57M+xparxkwDZgO5ARddxR/rx2BT4DdQs/3CLruKO/vSODy0OMDgEVB113VLchmmV7AM6HHzwCnV17BOTffObcg9HgZsBKosfN+HOkKLHTOfeOc2wSMw+93RRX/HV4Ajjczi2GNkVLjvjrnpjjnSkJPpwNtY1xjpITzewW4A38QsyGWxUVYOPt6KfCoc+4HAOfcyhjXGEnh7K8Dmoce7wIsi2F9YQsy3Fs555YDhO73qG5lM+uK/0v6dQxqi5Q2wJIKz5eGllW5jnNuC7AW2D0m1UVWOPta0V+AN6JaUfTUuK9m1hlo55x7NZaFRUE4v9dsINvMPjSz6WbWI2bVRV44+zsUON/MlgKvAwNiU1rtRHWCbDN7F9izipfya7md1sBzwIXOubJI1BYjVR2BV+6eFM46iSDs/TCz84Ec4OioVhQ91e6rmdUDHgD6xaqgKArn91of3zRzDP7b2AdmdpBz7sco1xYN4exvH+Bp59zfzexw4LnQ/sZVNkU13J1z3bf3mpmtMLPWzrnlofCu8qucmTUHXgNuds5Nj1Kp0bIUaFfheVu2/QpXvs5SM6uP/5q3JjblRVQ4+4qZdcf/cT/aObcxRrVFWk372gw4CJgaamHbE5hkZqc55wpjVmVkhPsZnu6c2wx8a2bz8GE/MzYlRlQ4+/sXoAeAc+4jM2uMH3cmrpqjgmyWmQRcGHp8IfBy5RXMrCHwEvCsc25CDGuLlJlARzNrH9qX3vj9rqjiv8PZwHsudKYmwdS4r6GmiieA0xK8XbbafXXOrXXOtXTOZTnnsvDnFxIx2CG8z/BE/MlyzKwlvpnmm5hWGTnh7O9i4HgAM9sfaAwUx7TKcAR4Vnp3YDKwIHTfIrQ8BxgVenw+sBn4tMKtU9BnoWu5nz2B+fhzBfmhZbfj/7OD/2BMABYCHwMdgq45ivv6LrCiwu9yUtA1R2tfK607lQTtLRPm79WA+4Evgc+B3kHXHOX9PQD4EN+T5lPgxKBrruqmK1RFRJKQrlAVEUlCCncRkSSkcBcRSUIKdxGRJKRwFxFJQgp3EZEkpHAXEUlCCncRkST0/0DHejllrrA0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0e0c01d160>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 20\n",
    "w = [1.0]*n\n",
    "l = {(i,i+1): 0.09 for i in range(n-1)}\n",
    "l.update({(5,14): 0.3})\n",
    "f = {0: (0.0,1.0), 13: (0.5,0.9), 17: (0.7,1.1)}\n",
    "g = 9.81\n",
    "\n",
    "x,c = stringModel(w, l, f, g)\n",
    "if x is not None:\n",
    "    display(x, c, list(l.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Duality and feasibility\n",
    "\n",
    "The dual problem is as follows:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{array}{ll}\n",
    "maximize & -\\sum_{ij\\in L}\\ell_{ij}y_{ij} - \\sum_{i\\in F}f_i\\circ z_i\\\\\n",
    "s.t.     & y_{ij}\\geq \\|v_{ij}\\|,\\ ij\\in L \\\\\n",
    "         & \\sum_{j~:~ij\\in L} v_{ij}\\mathrm{sgn}_{ij} + \\left(\\begin{array}{c}0\\\\ gw_i\\end{array}\\right) +z_i = 0, \\ i=1,\\ldots,n\n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\mathrm{sgn}_{ij}=+1$ if $i>j$ and $-1$ otherwise and $\\circ$ is the dot product. The variables are $(y_{ij},v_{ij})\\in \\mathbf{R}\\times\\mathbf{R}^2$ for $ij\\in L$ and $z_i\\in\\mathbf{R}^2$ for $i\\in F$ (we assume $z_i=0$ for $i\\not\\in F$).\n",
    "\n",
    "Obviously (!) the linear constraints describe the equilibrium of forces at every mass. The ingredients are: the vectors of forces applied through adjacent strings ($v_{ij}$), gravity, and the attaching force holding a fixed point in its position. By proper use of vectorization this is much easier to express in Fusion than it looks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dualStringModel(w, l, f, g):\n",
    "    n, m = len(w), len(l)\n",
    "    starts = [ lKey[0] for lKey in l.keys() ]\n",
    "    ends = [ lKey[1] for lKey in l.keys() ]\n",
    "\n",
    "    M = Model(\"dual strings\")\n",
    "\n",
    "    x = M.variable(Domain.inQCone(m,3))       #(y,v)\n",
    "    y = x.slice([0,0],[m,1])\n",
    "    v = x.slice([0,1],[m,3])\n",
    "    z = M.variable([n,2])\n",
    "\n",
    "    # z_i = 0 if i is not fixed\n",
    "    for i in range(n):\n",
    "        if i not in f:\n",
    "            M.constraint(z.slice([i,0], [i+1,2]), Domain.equalsTo(0.0))\n",
    "\n",
    "    B = Matrix.sparse(m, n, list(range(m))+list(range(m)), starts+ends, [1.0]*m+[-1.0]*m).transpose()\n",
    "    w2 = Matrix.sparse(n, 2, range(n), [1]*n, [-wT*g for wT in w])\n",
    "\n",
    "    # sum(v_ij *sgn(ij)) + z_i = -(0, gw_i) for all vertices i\n",
    "    M.constraint(Expr.add( Expr.mul(B, v), z ), Domain.equalsTo(w2))\n",
    "\n",
    "    # Objective -l*y -fM*z\n",
    "    fM = Matrix.sparse(n, 2, list(f.keys())+list(f.keys()), [0]*len(f)+[1]*len(f), \n",
    "                       [pt[0] for pt in f.values()] + [pt[1] for pt in f.values()])\n",
    "    \n",
    "    M.objective(ObjectiveSense.Maximize, Expr.neg(Expr.add(Expr.dot(list(l.values()), y),Expr.dot(fM, z))))\n",
    "    M.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us quickly discuss the possible situations regarding feasibility:\n",
    "\n",
    "* The system has an equilibrium --- the problem is **primal feasible** and **dual feasible**.\n",
    "* The strings are too short and it is impossible to stretch the required distance between fixed points --- the problem is **primal infeasible**.\n",
    "* The system has a component that is not connected to any fixed point, hence some masses can keep falling down indefinitely, causing the problem **primal unbounded**. Clearly the forces within such component cannot be balanced, so the problem is **dual infeasible**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Springs\n",
    "\n",
    "We can extend this to consider infinitely strechable springs instead of fixed-length strings connecting the masses. The next model appears in [Applications of SOCP](#http://stanford.edu/~boyd/papers/pdf/socp.pdf) by Lobo, Boyd, Vandenberghe, Lebret. We will now interpret $\\ell_{ij}$ as the base length of the spring and assume that the elastic potential energy stored in the spring at length $x$ is \n",
    "\n",
    "$$\n",
    "E_{ij}=\\left\\{\\begin{array}{ll}0 & x\\leq \\ell_{ij}\\\\ \\frac{k}{2}(x-\\ell_{ij})^2 & x>\\ell_{ij}\\end{array}\\right.\n",
    "$$\n",
    "\n",
    "That leads us to consider the following second order cone program minimizing the total potential energy:\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{array}{ll}\n",
    "minimize & g\\cdot \\sum_i w_ix_i^{(2)} + \\frac{k}{2}\\sum_{ij\\in L} t_{ij}^2 \\\\\n",
    "s.t.     & \\|x_i-x_j\\|\\leq \\ell_{ij}+t_{ij},\\ ij\\in L \\\\\n",
    "         & 0\\leq t_{ij},\\ ij\\in L \\\\\n",
    "         & x_i = f_i,\\ i\\in F\n",
    "\\end{array}\n",
    "\\end{equation}\n",
    "\n",
    "If $t$ denotes the vector of $t_{ij}$ then using a rotated quadratic cone for $(1,T,t)$:\n",
    "\n",
    "$$\n",
    "2\\cdot 1\\cdot T\\geq \\|t\\|^2\n",
    "$$\n",
    "\n",
    "will place a bound on $\\frac12\\sum t_{ij}^2$. We now have a simple extension of the first model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# w - masses of points\n",
    "# l - lengths of strings\n",
    "# f - coordinates of fixed points\n",
    "# g - gravitational constant\n",
    "# k - stiffness coefficient\n",
    "def elasticModel(w, l, f, g, k):\n",
    "    n, m = len(w), len(l)\n",
    "    starts = [ lKey[0] for lKey in l.keys() ]\n",
    "    ends = [ lKey[1] for lKey in l.keys() ]\n",
    "\n",
    "    M = Model(\"strings\")\n",
    "    x = M.variable(\"x\", [n, 2])                 # Coordinates\n",
    "    t = M.variable(m, Domain.greaterThan(0.0))  # Streching\n",
    "\n",
    "    T = M.variable(1)                           # Upper bound\n",
    "    M.constraint(Expr.vstack(T, Expr.constTerm(1.0), t), Domain.inRotatedQCone())\n",
    "\n",
    "    # A is the signed incidence matrix of points and strings\n",
    "    A = Matrix.sparse(m, n, list(range(m))+list(range(m)), starts+ends, [1.0]*m+[-1.0]*m)\n",
    "\n",
    "    # ||x_i-x_j|| <= l_{i,j} + t_{i,j}\n",
    "    c = M.constraint(\"c\", Expr.hstack(Expr.add(t, Expr.constTerm(list(l.values()))), Expr.mul(A, x)), \n",
    "        Domain.inQCone() )\n",
    "\n",
    "    # x_i = f_i for fixed points\n",
    "    for i in f:\n",
    "        M.constraint(x.slice([i,0], [i+1,2]), Domain.equalsTo(list(f[i])).withShape([1,2]))\n",
    "\n",
    "    # sum (g w_i x_i_2) + k*T\n",
    "    M.objective(ObjectiveSense.Minimize, \n",
    "        Expr.add(Expr.mul(k,T), Expr.mul(g, Expr.dot(w, x.slice([0,1], [n,2])))))\n",
    "\n",
    "    # Solve\n",
    "    M.solve()\n",
    "    if M.getProblemStatus(SolutionType.Interior) == ProblemStatus.PrimalAndDualFeasible:\n",
    "        return x.level().reshape([n,2]), c.dual().reshape([m,3])\n",
    "    else:\n",
    "        return None, None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3Xl8VNXdx/HPj7CEAIIQqCyGoAKWqlVL1frUpUUUN2iBWmi01VJRsY/W3TaVIja2olZrH0Hjrk3dilbqVtxwaYslVlzYBJEERGVfZAiQcJ4/zgRCmCST5M7cycz3/XrlNTN3bu78yCt8c+65555jzjlERCS9tAq7ABERCZ7CXUQkDSncRUTSkMJdRCQNKdxFRNKQwl1EJA0p3EVE0pDCXUQkDSncRUTSUOuwPjg3N9fl5+eH9fEiIi3SO++8s8Y5172h/UIL9/z8fEpLS8P6eBGRFsnMyuLZT90yIiJpSOEuIpKGFO4iImlI4S4ikoYU7iIiaUjhLiKShhTuIiJpSOEuIpKGFO4iEpySEsjPh1at/GNJSdgVZazQ7lAVkTRTUgLjx0Mk4l+XlfnXAAUF4dWVodRyF5FgFBZCJMIonuRHRFvskYjfLkmnlruIBKO8HAe8yKm0Yuce2yX51HIXkWDk5fEmxxGhA735dI/tknxquYtIMIqKmPqTtmRV7aAHq/y2nBwoKgq3rgylcBeRQHw+pIDp7KRn1he0qtoJffv6YNfF1FAo3EUkEPfdB5VVreh1VE9o3xNmLQu7pIymcBeRZqushLvvhpNOgh07wq5GQBdURSQAzz0Hy5fDhAlhVyLVFO4i0mxTp0Lv3nDmmWFXItUU7iLSLIsXw8yZcMEF0FodvSlD4S4izXLXXT7Uf/azsCuRmhTuItJkW7fCAw/A978PPXuGXY3UpHAXkSZ7/HFYv14XUlORwl1EmmzqVPjqV+GEE8KuRGpTuItIk8yZ478mTACzsKuR2hTuItIk06ZBhw5wzjlhVyKxKNxFpNHWrYNHH4Wzz4bOncOuRmJRuItIoz30EFRUwEUXhV2J1EXhLiKNsnOn75I59lj4+tfDrkbqonAXkUZ55RV/V6qGP6Y2hbuINMrUqZCbC6NHh12J1KfBcDez+81slZl9WMf7ZmZ3mNkSM3vfzI4MvkyJS0kJ5OdDq1b+saQk7IokzSxfDjNm+KkG2rULuxqpTzwt9weBYfW8fyrQP/o1HpjW/LKk0UpKYPx4KCsD5/zj+PEKeAlUcbH/9brggrArkYY0GO7OuTeAdfXsMgJ42HmzgS5mplkmkq2wEBeJMJuj+Rff8tsiESgsDLcuSRvbt8O998Jpp/kTQ0ltQUzQ2RtYXuP1iui2z2rvaGbj8a178rQierDKyzHgPB6gJ5/xKkN2bRcJwt/+Bp9/rgupLUUQF1Rj3XjsYu3onCt2zg12zg3u3r17AB8tu0T/WI5iOq9zAqvJ3WO7SHNNnQr9+sEpp9SzU0kJzJ4Nr8/SdZ+QBRHuK4D9a7zuA6wM4LjSGEVFkJPDKKazkyyeYQTk5PjtIs00bx68/jpceCFkZdWxU/V1n20V/rWu+4QqiHCfAfw4OmrmGGCjc26vLhlJsIICKC7m8Lz15PMJ07PP9le/CgrCrkxasugIrGmH3Ek7tvHTff5a976FhRCJsISDWMJBfpuu+4SmwT53M3sUOBHINbMVwG+ANgDOubuA54HTgCVABDgvUcVKAwoKsIICRl0Jd9zRjw2nQ5ewa5KWK9oS/zJiPMw5nMXj5F5xEXTaFrvRUNf1HV33CUU8o2XGOud6OufaOOf6OOfuc87dFQ12oqNkLnbOHeicO9Q5V5r4sqU+o0bBjh3w7LNhVyItWrQlfjU3sZl9uIC76m+JR6/vHBRtu9feLsmlO1TT0NFHQ69eMH162JVIi1ZeToT2PMS5tCfC0by9a3tM0es+e9B1n9Ao3NNQq1YwciS8+CJ8+WXY1UiLlZfHFK4mQgeeZDSt2blre0zR6z60y/av+/bVdZ8QKdzT1KhRfkrWF14IuxJpqcouu52buIYxPMrpRH+RGmqJFxTAMcfACSfCsmUK9hAp3NPUccdB9+7qmpGmu/Kt79GqXRum9L7Dr6OnlniLEsQdqpKCsrLge9/zq+VUVEB2dtgVSUvy2mvw17/CDTe0Zv9f/zvscqQJ1HJPY6NG+T73mTPDrkRakspKuOQSf4PpFVeEXY00lcI9jX3nO9Cli7pmpHHuvhs+/BD+8Ado3z7saqSpFO5prG1bGD7cz7+9fXvY1UhLsHYtXHcdDBniu/Wk5VK4p7mRI2HDBt+HKtKQ666DTZvgj3/011Cl5VK4p7mTT4YOHeCpp8KuRFLde+/5LpmLL4avfS3saqS5FO5prn17OP10Pxd3VVXY1Uiqcs5fRN13X5g0KexqJAgK9wwwahSsWgVvvRV2JZKqnngC3ngDbrzRB7y0fAr3DHDaaX6cu0bNSCyRCFx1FRxxBIwbF3Y1EhSFewbo2NGvnvPUU7BzZ9jVSKq56SZYvtxfRK1zIQ5pcRTuGWLUKPj0U/jPf8KuROISXSSDVq0SulzdsmUwZQqMHeunrJD0oXDPEGeeCW3aqGumRaherq6szF/pTOBydVde6f9+TJkS+KElZAr3DNGli78xZfp0nxeSwqKLZEzlIqZxgd+WgOXqXn3V/z788pfQp0+gh5YUoHDPIKNGwSefwNy5YVci9YouhvFLfsev+N1e24NQPX9Mv36+9S7pR+GeQUaM8Kfg6ppJcbsWw3C4mNubb9o0mDcPbr014BlDS0pg9mx4fVZCrxVIwxTuGaR7dzj+eIV7ytu1XJ3hqv+LBrhc3Zo1MHEinHRSwPPHVF8r2FbhXyfwWoE0TOGeYUaNgoULYf78sCuROkWXq3PWiiqyAl8k47rrYPPmBMwfE71WsIcEXCuQ+CjcM8z3v+8fNddMiisoYGdOJ6qyOwa6XN3cubvnjxk0KJBD7ha9JnA4czmcuXttl+RSuGeY3r3hW99S10xLUFXlL3wGoqQE1zefS454g262lklfezKgA9cQvSZwO5dxO5fttV2SS+GegUaN8i24pUvDrkTqUx3uzR66Gu0Lf7T8W7zJ8RTt/CX7XnZu8H3hu64V1BDgtQJpHIV7Bho50j9O/8aNCb8DUpquehbPiopmHqiwkE8i3TmXhziQxYzjvsT0hUevFdC3rxbUTgFaIDsD9ftXCUfaIKZv+C5XUeMOSNB/xBRRWbl7HqBIpHnL3a0uizCMN6gii0u4gyyiB05EX3hBgX6HUoRa7pmosJDj3Cze5hj+y+F+m0Y1pJRt23Y/37Kl6cdZswaGtHmdcvoyk6Fcwv/tflN94WlN4Z6JyssZxosA/B//u8d2SQ01u2Jqjy6M19q1fsqJxfTn7+1+wBBqrLWovvC0p3DPRHl5DGMmJ/ESLzGUSrJ2bZfU0NyW+7p1/ialRYtgxnOtOem+seoLzzAK90wUHdVwMXeygv15ljPUkksxNVvujQ336mBfsACeeQaGDsUH+bJlviM/wHHzkroU7pkoOqrhjLwP6MNypmVfrpZcimlqt8z69X5R9Hnz/Lq5p5wSfG3SMijcM1VBAa3LPuaCG/ZnZsXxLD5KwZ5KmtJy37DBB/sHH8DTT8OwYYmpTVoGhXuG+9nPoHVruOuusCuRmhrbct+40Qf7e+/5u49POy1xtUnLoHDPcPvt5+9YfeCBpo/KkOA1puW+caPvfpk71wf7GWcktjZpGRTuwoQJvq/28cfDrkSqxRvumzb57pd33oEnn/TLKYqAwl3wCyN/7WswdWrYlUi1eLplNm+GU0+F0lJ44gm/GItINYW7YAYXXeRDYs6csKsR2B3ubdrEbrlXB/vbb8Njj+2eylmkWlzhbmbDzGyRmS0xs2tjvJ9nZq+Z2btm9r6Z6XJOC3POOdChg1rvqaI63HNy9m65f/klnH66X83u0Uf9NROR2hoMdzPLAu4ETgUGAWPNrPY0/78GnnDOHQGMARQRLcw++/iAf+wxf9u6hKs63Dt02LPlvmWLD/Z//Qv+8hf4wQ/CqU9SXzwt96OAJc65pc657cBjQO3ePQfsE33eGVgZXImSLBdd5EPlwQfDrkRitdy3bPEjYd56C/78ZzjrrPDqk9QXT7j3BpbXeL0iuq2mScDZZrYCeB5qzka1m5mNN7NSMytdvXp1E8qVRDrsMPj2t2HatN3TzUo4qsO9Y0cf6pGIHwnzxhvwyCMwZky49UnqiyfcYy2hW3ttmLHAg865PsBpwCNmttexnXPFzrnBzrnB3bt3b3y1knATJsDHH8NLL4VdSWbb1S2zYA6b/v4aw7u9xeuzdvLww/CjH4Vbm7QM8YT7CmD/Gq/7sHe3yzjgCQDn3L+BbCA3iAIluUaOhB49dGE1bBX/nU82W8natoV3OYJXK47lwTbjKUArZkl84gn3OUB/M+tnZm3xF0xn1NqnHBgCYGZfxYe7+l1aoHbt/JQEzz6r6d3DVPHG22RTwX84ik105l7Gcc72+7SgisStwXB3zlUCPwf+ASzAj4qZZ2aTzWx4dLcrgPPN7D3gUeBc55q9rK+EpHrFveLicOvIZBWbK2lFFRXkcBjv8VMe9G/oL67EKa41VJ1zz+MvlNbcNrHG8/nA/wRbmoSlb18/KuOee2DiRGjbNuyKMk9FTle2RDrSgy94i+N2v6EFVSROukNVYpowAVatgqeeCruSzLS093FsI5vJTKQTX/qNWlBFGkHhLjENHQoHHqgLq2HYsQPmruxBuzZVjMt7WUvjSZPE1S0jmadVK7jwQrjqKr/4w6GHhl1R5rj7bj+2feDALFov/DjscqSFUstd6nTeeX70zLRpYVeSOTZuhOt/VUFn28h+i16D/Hwo0fBHaTyFu9SpWzd/J+Qjj/h5wyXxfn/OPNZszqaPK6c9FVBW5ocvKeClkRTuUq8JE/wshH/+c9iVpL/ycrjt7wdxDg+zjm6sJnoXdySi8e3SaAp3qdc3vwnf+Ia/sKo7FxKrsBCMnfyWX7OOrqyh2+43Nb5dGknhLvUy8633efPgzZ5n+Sut6gcO3Dvv+LOjy/a5nzyW053VfJfXdu+g8e3SSAp3adAYHqML65n6xUjffFc/cKCcgyuvhNxcuOaW7pCTQ4QccojO9avx7dIECndpUM7kazmXB3iS0XzEQX6j+oED89xzMGsWTJoEnc8/C4qLiVgHctiq8e3SZAp3aVh5OUN5iZ20Zjz37LFdmqey0t9LMGDA7jl9do4toMJlk/Obq2DZMgW7NInCXRqWl8dpvMjp/J1/cSyLq1vv6gdutnvvhYULYcoUvxg2wNat/jEnJ7y6pOVTuEvDioogJ4d7OZ92bONKblE/cAA2bYLf/AaOPx6GD9+9PVKjq12kqRTu0rCCAiguZr++2RRyIzMYwcuXzFB3QTNNmeInZ7vlFj8qqZrCXYKgcJf4FBTAsmX8Yuvv6NcPLnt2CJWVYRfVcq1YAX/4A4wd6+8lqEnhLkFQuEujZGfDzTfDhx/6/mJpmuuug6oquPHGvd9Tn7sEQeEujTZyJJxwgg+oDRvCrqblmTsXHnoILr3U3w9Wm1ruEgSFuzSaGdx+O6xdC5Mnh11Ny1J9w9K++8KvfhV7H4W7BEHhLk1y+OEwbhz86U/w0UdhV9NyvPgivPKKHyXTpUvsfRTuEgSFuzTZb38L7dvDFVeEXUnLUH3D0kEH+YVQ6qJwlyAo3KXJvvIV+PWv4dlnYebMsKtJfQ8+6Cdg+/3v6190XOEuQVC4S7Nceqlfa/Xyy9HQyHp8+aW/AH3ssf6CdH0U7hIEhbs0S7t2/iacefP8/FYS2y23wOefw6237nnDUiwKdwmCwl2abcQI+M53YOJEWL8+7GpSz8qV/t6As86CY45peP9IxP8BaNcu8bVJ+lK4S7NVD41cvx6uvz7salLPxImwYwf87nfx7R+J+FZ7Qy18kfoo3CUQhx0GP/sZ3Hmnn+VQvA8+gPvvh5//HA44IL7vqQ53keZQuEtgbrjBh5KGRu529dXQubMfVRQvhbsEQeEugenRw48Ief55f7NOpps50/8crrsOunaN//sU7hIEhbsE6pJL/E06l1/u+5kzVVWVv2GpXz+4+OLGfa/CXYKgcJdAtW3rh/stWAB33x12NeF5+GF4/31/w1JjR71s3apwl+ZTuEvgzjwThgzx86esWxd2Ncm3ZYvvYz/6aPjBDxr//Wq5SxAU7hI4M7jtNj8d8KRJYVeTfLfd5se2115hKV6RiJ+zR6Q5FO6SEIceCuPHw9SpvosmU3z+ue+KGTkSvv3tph1DLXcJgsJdEmbyZOjYbjuXHzkLWrXyK1OUlIRdVkJNmgTbtvmAbyqFuwRB4S4J031mCRN3TOTFihN51J0FZWW+OZ+mAT9/PtxzD0yYAP37N/04CncJgsJdEqewkJ/v+AP7sIFzeYhP6emTq7Aw7MoS4uqroVMnP669yUpKiGzYRs4dv8+IMx1JHIW7JE55OW3ZwUQm4zDO4Fk2sg+Ul4ddWeBefRWee84vnZeb28SDlJRQef5FbKcdOWxJ+zMdSay4wt3MhpnZIjNbYmbX1rHPWWY238zmmdlfgi1TWqS8PACu4DaeYQQfcijDmcHWPs3os0hBO3f6dVHz8vxNXE1WWMi6rX5Q/Kf08tvS+ExHEqvBcDezLOBO4FRgEDDWzAbV2qc/8Evgf5xzXwN+kYBapaUpKtrVeXwqL/II5/Amx3FW7itpdfdqSQm8+66f9TE7uxkHKi9nFT0A2EKHPbaLNFY8LfejgCXOuaXOue3AY8CIWvucD9zpnFsP4JxbFWyZ0iIVFPgVPPr2BTPG9J3N1PNKefbdPpx3nm/xtnRbt/qumMGDYcyYZh4sL4+2+L96pzBzj+0ijRVPuPcGltd4vSK6raYBwAAz+6eZzTazYbEOZGbjzazUzEpXr17dtIqlZSkogGXLfJIvW8aF9x/FjTf61u6ll4JzYRfYPLffDitW+BuWWjX3ClZREZuzuwOwD5v8tpwcfwYk0kit49gn1j12tf9Ltgb6AycCfYA3zewQ59yGPb7JuWKgGGDw4MEt/L+1NNW11/ppCW65xc+W2FIX+Fi1ynfFDB8OJ5wQwAELCtg0/ytwI3TiS3/GU1Tk/0CKNFI84b4C2L/G6z7Ayhj7zHbO7QA+MbNF+LCfE0iVklbMYMoUv3LT5Mmw777wixZ4leb66/31zptuCu6Ym48+CYBOc16FwcEdVzJPPCeSc4D+ZtbPzNoCY4AZtfb5G/AdADPLxXfTLA2yUEkvZn7WyFGj4LLL4KGHwq6ocRYt8vVfcAEcfHBwx9282T/us09wx5TM1GDL3TlXaWY/B/4BZAH3O+fmmdlkoNQ5NyP63slmNh+oAq5yzq1NZOHS8mVl+b73jRth3Di/YtH3vhd2VfG55hrfHf6b3wR73E3RrvZOnYI9rmSeeLplcM49Dzxfa9vEGs8dcHn0SyRu7drB00/D0KHwwx/CCy/Ad78bdlX1e/11eOYZuPFGv/pUkKpb7gp3aS7doSqh69jR3905YACMGAFzUvhKTfUNS336JOY6webNvsuqQ4eG9xWpj8JdUkLXrn7N0R494NRT/SRcqeixx6C01A9iScSc65s3+1Z7U+aBF6lJ4S4po2dPeOklv1TfySf74fGppKLC37B0+OFw9tmJ+YxNm9QlI8FQuEtKOeAA34KPRHw//BdfhF3Rbn/6k5/L69ZbA7hhqQ7VLXeR5lK4S8o55BB4/nm/VN0pp/jl+sK2Zo3vijn99MRe8FW4S1AU7pKSjjkG/vY33/d+xhm+JR+mG27wwTtlSmI/Z/NmjXGXYCjcJWUNHQp/+Qv8+98w+thP2d63fyjL9S1e7NeCPf98GDSo4f2bQ33uEhSFu6S00aOh+KezeeG93vy4/AaqnCV9EYtrr/Xj8SdNSvxnqVtGgqJwl5Q37qUx3MyVPM4YhvECFbRN2iIW//wnPPWUvyN1v/0S/nEKdwmMwl1SX3k5V3IrZ/IML3MyPfmca/kdS8uyEvqxzsEVV0CvXnB5Eu69dk597hIchbukvuhiFU/zfX7P1ZzILG7hSg7kY045xU9fkIiVnZ58Et5+G3772+TcMbptm/93qOUuQVC4S+qLLteXheMabuZpRlKWfTDXj3qP+fNh5Eg/9fnEicGtSLdtm+9rP+ww+PGPgzlmQzSvjARJ4S6pr9ZyffTtS+97r2fiX7/OJ5/AjBlwxBG+hd2vH5x5pp+rpqqq6R95553wySdw881+9spk0HS/EiRzIa1zNnjwYFdaWhrKZ0t6WrYM7rkH7rvP39mal+eHL44b56c2iNe6dXDggXD00fDiiwkrdy9z5/o/UtOn+7MRkVjM7B3nXINLuajlLmkjP9/34Cxf7vvLBwyA667zIT96tJ+3Jp5FuYuK/Hjzm29OeMl7ULeMBEnhLmmnTZvdYf7RR35q3lmz/GRkAwb4u0zrWp996VI/h8x558Ghhya1bIW7BErhLmmtf3/fAl+xwt/z1KuXH7Pepw/86Ed+4Q3n8G/m5/PLAx+nTWWEyYc/lfRa1ecuQVK4S0bIzvZh/sYbMG8eXHihX/XpxBNhUJ+N/PG8ufyjbABP8EOuclPodc05SZ3iALTEngRL4S4ZZ9Ag+OMf4dNP4YEHoPOapfxix80M4x+0ZRtXckvS7oDdpaSEzdfcAECn/zks6X9YJP0o3CVj5eTAuefC7B3fYDZH0Ybt9OJTOrLF7xDUoPmGlJTA+PFsXl8JQMfl85M6d46kJ4W7SF4eRzOHb/Bf+rNkj+1JUVgIkQjv8XXasJ3WVCX/zEHSjsJdJHoHbDfWsoZcvy0nx29PhugZwgIOpoqsvbaLNEXrsAsQCV1BAQC5F1Tw/pZcfydsUdGu7QmXlwdlZRzMIlpTued2kSZSuIsAFBSQ+y6snUbyV+YuKoLx41kTyaU7a/y2ZJ45SFpSt4xIVLduvqt769Ykf3B07pw1rXuSyxp/5lBcnLwzB0lLarmLROVGu9vXrvU3OSVVQQFrfgG5PzgIpp6V5A+XdKSWu0hUt27+cc2a5H92VZWfsKz6D4xIcyncRaKqgzWMcN+wwU9qpnCXoCjcRaJqdsskW/UfFIW7BEXhLhIVZreMwl2CpnAXiera1T+GGe7Vf2BEmkvhLhLVpg106aJuGUkPCneRGrp1U7eMpAeFu0gNubnhtdyzs/2NqSJBULiL1BBmyz03F8yS/9mSnhTuIjXk5oYb7iJBUbiL1BBWt8zatQp3CVZc4W5mw8xskZktMbNr69lvtJk5MxscXIkiydOtG2zZAhUVyf1ctdwlaA2Gu5llAXcCpwKDgLFmNijGfp2AS4C3gy5SJFnCuktV4S5Bi6flfhSwxDm31Dm3HXgMGBFjvxuAKUCS2zwiwQnjLtXKSli/XuEuwYon3HsDy2u8XhHdtouZHQHs75x7NsDaRJIujMnD1q3b87NFghBPuMcanOV2vWnWCrgNuKLBA5mNN7NSMytdvXp1/FWKJEkY3TK6gUkSIZ5wXwHsX+N1H2BljdedgEOAWWa2DDgGmBHroqpzrtg5N9g5N7h79+5Nr1okQbrNmg7Amh9eDPn5UFKS8M9UuEsixBPuc4D+ZtbPzNoCY4AZ1W865zY653Kdc/nOuXxgNjDcOVeakIpFEqWkhG5X/RSAtXSFsjIYPz7hAa9wl0RoMNydc5XAz4F/AAuAJ5xz88xsspkNT3SBIklTWEibrZvIJsLbHOW3RSJQWJjQj1W4SyLENc7dOfe8c26Ac+5A51xRdNtE59yMGPueqFa7tEjl5QBU0ZoPOHSv7QlRUsKaa6YA0O3YgUnpBpLMoDtURarl5QEwlJl0ZtNe2wNXUgLjx7Nww1doT4Ts8o+S0g0kmUHhLlKtqAhychjEAhbTn52Yn6axqCgxn1dYCJEIb3Ec22nrtyWhG0gyQ+uwCxBJGQUFAAy85H0q1rWnvPex5N900a7tgYt29+SwhSN5Z6/tIs2hlrtITQUFDHj6JgAW3fdW4oIdIC8PByyjH99i9h7bRZpL4S5Sy8CB/nHRogR/UFERK7L7s4WOHMxCvy2R3UCSURTuIrX06AGdOych3AsKWHjpNAAOZhH07QvFxYk9W5CMoT53kVrMfOs94eEOLOw9BICDV74KPRP/eZI51HIXiSFp4b7QnyXst1/iP0syi8JdJIaBA2HFCr9wRyItXAgHH6y1UyV4CneRGKovqn70UWI/Z8ECH+4iQVO4i8SQjBEzGzfCZ58p3CUxFO4iMRx0kO8qSWS4Vx9b4S6JoHAXiaF9ez8yMZHhvjA6tP2rX03cZ0jmUriL1CHRI2YWLoTWreGAAxL3GZK5FO4idRg40F9Qda7hfZtiwQLf/dOmTWKOL5lN4S5Sh4ED4csvYeXKhvdtiuphkCKJoHAXqUMiR8zs2AFLlqi/XRJH4S5Sh0SG+9KlUFmplrskjsJdpA69e0OHDokJ9wUL/KPCXRJF4S5SBzMYMCAx4V49DLL67EAkaAp3kbqUlDBw0QwWvbgU8vMDXdt04ULo2dNPGiaSCAp3kViii1fnR+bxCfmsLdsU6OLVCxfqYqokluZzF4klunj1NtoCrchlHZ0jG+h13mp6PQC9eu351bPn7sfs7HqOW1KC+1UhC8vnUtBxBpRkaXEOSQiFu0gs0UWqx/IYEXLowSo2sC8rd/RiZaQ/b77px79v3773t3btunf49+oFvRa/Tq+7isna1o2NdOHgL+fA+Pv9NyngJWAKd5FY8vKgrIxvUso3Kd29vW9f+NdowN+5um6dD/m6vubP9zM/VlUBnAC8vutQlWRBJOLPEhTuEjCFu0gsRUW+jz0S2b2t1uLVZtCtm/869NC6D1VVBWvWwMr9jmQlPXmXw3mZoQzlZb9D9CxBJEjmEjVxRgMGDx7sSktLG95RJCwlJb5VXV7uW/JFRc1rYefnQ1nZ3tv79oVly5p+XMkoZvaOc25wQ/tptIxIXQoKfOju3Okfm9t1UlTkW/811TobEAmKwl0kWQoKoLjYt9TN/GNxsfrbJSHU5y6STAUFCnNJCrXcRUTSkMJdRCQNKdxFRNIzt9J4AAAEvUlEQVSQwl1EJA0p3EVE0pDCXUQkDSncRUTSkMJdRCQNhTa3jJmtBmJMtJFScoE1YRfRBKo7uVR3cmV63X2dc90b2im0cG8JzKw0ngl6Uo3qTi7VnVyqOz7qlhERSUMKdxGRNKRwr19x2AU0kepOLtWdXKo7DupzFxFJQ2q5i4ikIYV7DWbW1cxeMrPF0cd9Y+xzuJn928zmmdn7ZvbDMGqN1jLMzBaZ2RIzuzbG++3M7PHo+2+bWX7yq9xbHHVfbmbzoz/fV8ysbxh11tZQ3TX2G21mzsxSYkRHPHWb2VnRn/k8M/tLsmuMJY7fkzwze83M3o3+rpwWRp21arrfzFaZ2Yd1vG9mdkf03/S+mR2ZsGKcc/qKfgFTgGujz68FboqxzwCgf/R5L+AzoEsItWYBHwMHAG2B94BBtfaZANwVfT4GeDwFfsbx1P0dICf6/KKWUnd0v07AG8BsYHBLqBvoD7wL7Bt93aOF1F0MXBR9PghYlgJ1Hw8cCXxYx/unAS8ABhwDvJ2oWtRy39MI4KHo84eA79XewTn3kXNucfT5SmAV0OANBQlwFLDEObfUObcdeAxff001/z1/BYaYmSWxxlgarNs595pzLhJ9ORvok+QaY4nn5w1wA76RUJHM4uoRT93nA3c659YDOOdWJbnGWOKp2wH7RJ93BlYmsb6YnHNvAOvq2WUE8LDzZgNdzKxnImpRuO/pK865zwCijz3q29nMjsK3Kj5OQm219QaW13i9Irot5j7OuUpgI9AtKdXVLZ66axqHb+mErcG6zewIYH/n3LPJLKwB8fy8BwADzOyfZjbbzIYlrbq6xVP3JOBsM1sBPA/8b3JKa5bG/v43WcatoWpmLwP7xXirsJHH6Qk8AvzEObcziNoaKVYLvPbQp3j2Sba4azKzs4HBwAkJrSg+9dZtZq2A24Bzk1VQnOL5ebfGd82ciD9LetPMDnHObUhwbfWJp+6xwIPOuVvN7FvAI9G6w/j/GK+k/Z/MuHB3zp1U13tm9oWZ9XTOfRYN75inp2a2D/Ac8OvoqVUYVgD713jdh71PS6v3WWFmrfGnrvWdMiZDPHVjZifh/+Ce4JzblqTa6tNQ3Z2AQ4BZ0Z6v/YAZZjbcOVeatCr3Fu/vyWzn3A7gEzNbhA/7OckpMaZ46h4HDANwzv3bzLLx87ekQrdSXeL6/Q+CumX2NAP4SfT5T4Bnau9gZm2Bp/H9Zk8msbba5gD9zaxftKYx+PprqvnvGQ286qJXdULUYN3R7o27geEp0v8LDdTtnNvonMt1zuU75/Lx1wrCDnaI7/fkb/iL2JhZLr6bZmlSq9xbPHWXA0MAzOyrQDawOqlVNt4M4MfRUTPHABuru4IDF/bV5VT6wvdHvwIsjj52jW4fDNwbfX42sAOYW+Pr8JDqPQ34CN/nXxjdNhkfKuB/2Z8ElgD/AQ4I+2ccZ90vA1/U+PnOCLvmeOqute8sUmC0TJw/bwP+AMwHPgDGhF1znHUPAv6JH0kzFzg5BWp+FD+Cbge+lT4OuBC4sMbP+s7ov+mDRP6O6A5VEZE0pG4ZEZE0pHAXEUlDCncRkTSkcBcRSUMKdxGRNKRwFxFJQwp3EZE0pHAXEUlD/w/zwcqrgqfKaQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f0dd049e1d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 20\n",
    "w = [1.0]*n\n",
    "l = {(i,i+1): 0.09 for i in range(n-1)}\n",
    "l.update({(5,14): 0.3})\n",
    "f = {0: (0.0,1.0), 13: (0.5,0.9), 17: (0.7,1.1)}\n",
    "g = 9.81\n",
    "k = 800\n",
    "\n",
    "x, c = elasticModel(w, l, f, g, k)\n",
    "if x is not None:\n",
    "    display(x, c, list(l.keys()))"
   ]
  },
  {
   "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). "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}