test_run.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import vectfit as m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Function vectfit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.619933845662216e-11"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# construting signal\n",
    "Ns = 101\n",
    "s = np.linspace(3., 7., Ns)\n",
    "poles = [5.0+0.1j, 5.0-0.1j]\n",
    "residues = [[0.5-11.0j, 0.5+11.0j],\n",
    "            [1.5-20.0j, 1.5+20.0j]]\n",
    "f = m.evaluate(s, poles, residues)\n",
    "weight = 1.0/f\n",
    "poles_init = [3.5 + 0.035j, 3.5 - 0.035j]\n",
    "\n",
    "# run vector fitting\n",
    "# single iteration\n",
    "poles_fit, residues_fit, cf, f_fit, rms = m.vectfit(f, s, poles_init, weight)\n",
    "\n",
    "# rms error\n",
    "rms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "### evaluate at finer points\n",
    "s_test = np.linspace(3., 7., 5*Ns)\n",
    "f_test = m.evaluate(s_test, poles, residues, cf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xt8VOW99/3Pb3KcEJJAEkg4gyAK\nCqlSCqVaWs/YotCDClar9bYgtdZu9966u9vejz5u9bndu9ZWYVOt2mqr3RXqoaD1UDw8hVK0VEVU\nNCAEcoBAAoSc57r/mJk4hAlMMmfyfb9e85o1a1bW+rFC5jvXOlyXOecQERHpzpPsAkREJDUpIERE\nJCwFhIiIhKWAEBGRsBQQIiISlgJCRETCUkCIiEhYCggREQlLASEiImFlJruAaJSUlLgxY8YkuwwR\nkbTyxhtv7HHOlR5rubQOiDFjxrBhw4ZklyEiklbM7ONIltMhJhERCUsBISIiYSkgREQkrLQ+ByEi\n/Ud7eztVVVW0tLQku5S0kZuby4gRI8jKyurTzysgRCQtVFVVMXDgQMaMGYOZJbuclOeco76+nqqq\nKsaOHdundegQk0gcVFfD5z8PNTXJruT40dLSQnFxscIhQmZGcXFxVC0uBYRIHNy2pIbXX/Nx63VK\niFhSOPROtPtLASESQ14vmMfH0o9ewDegjqUryzDzzxdJNwoIkRiqrIQvnPUAzL8CvjMRLwdZOPUd\ntm5NdmUivaeAEImh8nKom7jO/yJ3Py2FuykYV0JZWXLrkti49957Ofnkk1m4cCHDhw/nRz/6EeA/\nIfzd736X8ePHM2XKFN58800APvroIyoqKsjPz09m2X2mgBCJsaoBG8npKALg3HkvUYPSIWnWroU7\n7vA/x8D999/PqlWrmDBhAjfeeCO33norAKtXr2bLli1s2bKF5cuXs3jxYgBOOOEENm7cGJNtJ4MC\nQiSGnHP4ij7km5+5hCxPFp9a+BErViS7qn5q7Vo46yz44Q/9z1GGxKJFi6isrGTu3Ln85Cc/Oey9\np556iiuuuAIzY8aMGTQ0NFBdXR3V9lKBAkIkhnYf2s2BtgNMLp3M6KLRbG3QyYekWbMG2tqgs9P/\nvGZNVKtbtmwZw4YN489//jM33njjYe/t3LmTkSNHdr0eMWIEO3fujGp7qUABIRJDH+79EIDxg8cz\nqnAU2xu3J7mifmz2bMjOhowM//Ps2XHblHPuiHnHwyW5upNaJIYq91UCMG7QOEYVjuKFj15IckX9\n2MyZ8NJL/pbD7Nn+13EyYsQIduzY0fW6qqqKYcOGxW17iaIWhEgMVR/wH3ceXjCcUQWj2HVgF+2d\n7Umuqh+bORNuuSWu4QAwd+5cfvWrX+GcY926dRQWFlJeXh7XbSaCWhAiMVR9sJoBWQPIz85nVOEo\nHI6dB3YypmhMskuTOJozZw6rVq1i/Pjx5OXl8dBDDyW7pJhQQIjEUM3BGsry/Ze1Dhs4rGueAuL4\nsG3btrDzzYz77rsvscUkgA4xicRQaEAMGTAEgNqDtcksSeIgPz+f5cuXd90o15PgjXJDhw5NUGWx\npRaESAxVH6zmlCGnADA03/+hUNdUl8ySJA5uuukmbrrppmMupxvlRKRLzcEahg7wB0NpXimggJD0\npYAQiZH2znYaWhq6Di3lZOZQmFNIbZMOMUl6UkCIxMje5r0AFHuLu+YNzR+qFoSkrbgFhJmNNLM/\nm9lmM9tkZjcE5g82sxfMbEvgeVBgvpnZvWb2oZm9ZWanxas2kXiob64HoDjvk4AYMmCIWhCStuLZ\ngugA/sk5dzIwA1hiZpOAm4GXnHMTgJcCrwEuACYEHtcCS+NYm0jM1R8KBERIC6LYW9zVspDEi/XQ\nr6HdfXfX2trK2WefTUVFBU888QTXXHMN7777LgD/8R//0bXctm3b8Hq9VFRUdM177rnnmDhxIuPH\nj+fOO+/smr9w4UIGDx7M73//+9j8A3opbgHhnKt2zr0ZmD4AbAaGAxcBjwQWewS4ODB9EfAr57cO\nKDKz9L8VUfqNcC2Iwd7BCogkuu02eP11CPTKHbVgd9+PPfbYEe/9/e9/p729nY0bN3LJJZfwwAMP\nMGnSJODwgIDDr27q7OxkyZIlrF69mnfffZff/va3XcHy2GOPMXfu3NgU3wcJOQdhZmOATwF/BYY6\n56rBHyLAkMBiw4EdIT9WFZjXfV3XmtkGM9uwe/fueJYt0ivhWhCDvYPZ17wvWSX1W14vmMHSpeDz\n+Z+jHfr1aN1919XVcfnll7Nx40YqKir46KOPmD17Nhs2bODmm2+mubmZioqKsC2P9evXM378eMaN\nG0d2djaXXnopTz31VN8LjaG4B4SZ5QNPAt9zzu0/2qJh5h3RRaJzbrlzbppzblppaWmsyhSJ2p5D\newAoySvpmjcodxBN7U20drQmq6x+qbISFiyAvDz/67w8WLiQqIZ+PVp330OGDOGBBx7gjDPOYOPG\njZxwwgld79155514vV42btwYtuWRyl2FxzUgzCwLfzg85pwLDptSGzx0FHgOXuJRBYwM+fERwK54\n1icSS/XN9eRk5JCXldc1b7B3MAD7WtSKSKTycigogJYWyM31PxcUkJJDv6ZyV+HxvIrJgAeBzc65\n/wp562ngysD0lcBTIfOvCFzNNANoDB6KEkkHDS0NDPIOOuyPOxgQX/7avpidKJXI1NbCokWwbp3/\nOVX3fyp3FR7PrjZmAd8A3jaz4L3m/wbcCfzOzL4FbAe+FnhvFTAH+BA4BFwVx9pEYq6hpYHCnMLD\n5g3yDgLgjXf3cOt1Ndy/IgW/wh6nQod6TXY/ellZWbS3t5OVlXXEe5/+9KfZsmULW7duZfjw4Tz+\n+OP85je/SUKVR4pbQDjnXif8eQWAs8Is74Al8apHJN4aWxspyi3qeu31QsvgwXAtuNxGlq4sY6n5\nD3k0NyexUEm4a6+9lilTpnDaaadx++23H/ZeZmYmP//5zznvvPPo7Ozk6quvZvLkyUmq9HDqrE8k\nRhpbGinM/aQFUVkJ3/7qbp4B8O4ljybmTd3K3c+dkrQaJTo9dfcNMHv2bGaHDGu6JmQM7Lvuuou7\n7rqrx3XMmTOHOXPmxKjK2FFXGyIx0tDScFgLorwcSsrHApDpraWFXArGlaTkiVJJnIyMDBobGw+7\nUa4nCxcu5JVXXiE3NzcBlR1JLQiRGGlsbTziHMQ+34ngjKvPfI3MYd+g2ikd0t1DDz3ET3/608Pm\nzZo1K+IBg0aOHHnYSemjCXdZbCIpIERipHsLAmDlCg+D7yoie/pofnaXwuF4cNVVV3HVVf3jGhod\nYhKJgbbONlo6WvifXxcecTnlYO9g9raouw1JPwoIkRhobGkE4OP3C7j1usMTYpB3kLrbkLSkgBCJ\nktcLQ0Y1AOBaBrF0Zdlh/f6owz5JVwoIkShVVsJ5p77pf9FSSB5NLJz6Tle/PwoISVcKCJEolZdD\n1kj/PaHZLd4jLmcdlDtIfTEdJ442HsSaNWv4y1/+0qf1btu27bC7px9++GFKS0u55ppruubdcccd\njB8/nokTJ/L8888DdPUSm52dzZ49e/q07aPRVUwiMVCXmQHAo/NfZ832Uw67nDXYgnDOpUwnbNI3\n999/P6tXr2bs2LFHvLdmzRry8/P57Gc/2+v1BgNiwYIFXfMuueQSfv7znwPw7rvv8vjjj7Np0yZ2\n7drF2WefzQcffNDVS+yYMWP6/G86GgWESAxce30D65+B6T++iq8VHX45a1FuET7n42DbQQbmDExS\nhceX7z33PTbWbDz2gr1QUVbBPeff0+P7oeNBXH311Yd1+b1t2zaWLVtGRkYGjz76KD/72c846aST\nWLRoEdu3bwfgnnvuYdasWbzyyivccMMNgL/X1ldffZWbb76ZzZs3U1FRwZVXXsmgQYMO2/ZTTz3F\npZdeSk5ODmPHjmX8+PGsX7+emTNnxnQfdKeAEImBxlb/VUzd74MIndfQ0qCASGPLli3jueee489/\n/jMlJSWHvTdmzBgWLVpEfn4+N910EwALFizgxhtv5HOf+xzbt2/nvPPOY/Pmzdx9993cd999zJo1\ni4MHD5Kbm8udd97J3XffzbPPPgv4DzGF2rlzJzNmzOh6nagxIxQQIjHQ2NKIYWEDIDQgRhaOPOJ9\n6b2jfdNPFS+++GLX0KEA+/fv58CBA8yaNYvvf//7LFy4kPnz5zNixIhjritZY0YoIERioKGlgYKc\nAjx25HUfoQEh/YfP52Pt2rV4u41zevPNN3PhhReyatUqZsyYwYsvvnjMdSVrzAhdxSQSA42th/fk\nGkoB0T8MHDiQAwcOdL0+99xzu04yA2zc6D9n8tFHH3Hqqafyr//6r0ybNo333nvviJ/tbu7cuTz+\n+OO0traydetWtmzZwvTp0+P3jwlQQIjEQLjBgoIUEP3Dl7/8ZVauXElFRQWvvfYa9957Lxs2bGDK\nlClMmjSJZcuWAf6T1aeccgpTp07F6/VywQUXMGXKFDIzM5k6dSo/+clPjlj35MmT+frXv86kSZM4\n//zzue+++8jIyIj7v0mHmERioPtgQaEG5fqvSFFApL+jjQdx4okn8tZbbx0274knnjhiuZ/97Gdh\nf/6ll17qmu5+khrgBz/4AT/4wQ8iKzRG1IIQiYHugwWFCs5XQEikvF4vq1evPuxGuXCCN8q1t7fj\n8cT+41wtCJEYaGhpYPKQ8MNEZnoyyc/OV0AcJ6IdDyISl1xyCZdccskxlwveKBcvCgiRGAg3WFCo\notwidbcRA6lwN3o6jQcR7vLY3tAhJpEoOef8h5iOERBqQUQnNzeX+vr6qD/0+gvnHPX19VENV6oW\nhEiUmtqb6HSdPZ6kBgVELIwYMYKqqip2796d7FLSRm5ubkQ34vVEASESpeBgQT2dpAZ/QOzcH/+u\nEY5nWVlZYTvJk/jRISaRKAVbBmpByPFGASESpWBHfUc9B5GjgJD0o4AQiVKkLYjG1kZ8zpeoskSi\npoAQiVIk5yAGeQd1jQkhki4UECJRiugQk/pjkjSkgBCJUqSHmEKXFUkHCgiRKDW2NJLlySI3s+cb\nkoIBsa9Zd1NL+lBAiESpoaWBotyio3YBoRaEpCMFhEiUjjZYUJACQtKRAkIkSsfqqA8UEJKeFBAi\nUQoeYjqagpyCrmVF0oUCQiRKRxssKCjTk8nA7IEKCEkrcQsIM/ulmdWZ2Tsh8/63me00s42Bx5yQ\n924xsw/N7H0zOy9edYnEWkNLA0U5R29BQKA/plYFhKSPeLYgHgbODzP/J865isBjFYCZTQIuBSYH\nfuZ+M4v/iNwiMRDJSWrw302tFoSkk7gFhHPuVWBvhItfBDzunGt1zm0FPgSmx6s2kVjp8HVwsO3g\nMU9Sg3p0lfSTjHMQ3zGztwKHoAYF5g0HdoQsUxWYJ5LS9rfuB/ytg2Mpyi3SjXKSVhIdEEuBE4AK\noBr4z8D8cHcYhR1X0MyuNbMNZrZBI0tJskXSzUaQWhCSbhIaEM65Wudcp3POB/yCTw4jVQEjQxYd\nAezqYR3LnXPTnHPTSktL41uwyDH0KiA0JoSkmYQGhJmVh7ycBwSvcHoauNTMcsxsLDABWJ/I2kT6\norctiP2t+zUmhKSNuI1JbWa/BWYDJWZWBfwYmG1mFfgPH20Dvg3gnNtkZr8D3gU6gCXOuc541SYS\nK8GAiPQktcOxv3V/RIEikmxxCwjn3GVhZj94lOVvB26PVz0i8dDbFkTwZxQQkg50J7VIFPoaECLp\nQAEhEoWGlgYMY2DOwGMuq4CQdKOAEIlCQ0sDhbmFeOzYf0rBeyWW/FMDNTXxrkwkegoIkSj05nxC\ncLnNlXu59TolhKQ+BYRIFCINCK8Xxpb7l3O5jSxdWYaZf75IqlJAiEQh0oCorITLJm4HZ5DbQB5N\nLJz6Dlu3JqBIkT5SQIhEIdKAKC+HwrFDoLWAjNw9tJBLwbgSysoSUKRIH8XtPgiR/qA35yBqrYyB\nnoF8vuLvjMrbTbVTOkhqU0CIRCHSwYIAVqyAimXFeIpKue8OhYOkPh1iEumjDl8HB9oORDRYUJB6\ndJV0ooAQ6aPgWBC96TZDASHpRAEh0keNLY0ALP1JUcQ3vmnYUUknCgiRPgp+0G95qyDiG980JoSk\nEwWESB94vXDaZ/0f9K5lcMQ3vgXHhOjwdSSgSpHoKCBE+qCyEs446S3/i5aiiG98K84rBtDY1JIW\nFBAifVBeDhnD/MOmZ7fkRnzj22DvYADqm+vjXaJI1HQfhEgf1Wf4A+LFK17g8aqiiG58K/b6WxD1\nhxQQkvoUECJ9NO+yBt551Zh1zxLOiKC7b/jkEJNaEJIOdIhJpI8aWhooyCmIaCyIILUgJJ0oIET6\nqKG192NLqwUh6UQBIdJHwdHkemNg9kAyPZlqQUhaUECI9FFvenINMjMGewezt3lvnKoSiR0FhEgf\n7Wve1+uAAP95CB1iknSggBDpo/rm+q6Tzr1RnKeAkPSggBDpA+cc9Yf6GBDeYp2DkLSggBDpg0Pt\nh2jtbO26Kqk3dIhJ0oUCQqQPgieZ+3yI6VA9zrlYlyUSUwoIkT4ItgD62oJo7WzlUPuhWJclElMK\nCJE+CJ5DCHa+1xvBn9GlrpLqFBAifdDVgujjISaAPYf2xLQmkVjrMSDM7NeB5xsSV45Ieug6B9GH\nQ0xDBgwBFBCS+o7WgjjdzEYDV5vZIDMbHPpIVIEiqSh4iKkvLYhgQNQ11cW0JpFYO1p338uA54Bx\nwBuAhbznAvNF+qX65noGZA0gJzOn1z8bDIjaptpYlyUSUz22IJxz9zrnTgZ+6Zwb55wbG/JQOEi/\nVt9c36cT1ACFOYVkebLUgpCUd8yT1M65xYkoRCSd1B+q79P5B/B32DdkwBAFhKQ8XcUk0gd7m/f2\n6fxDkAJC0kHcAsLMfmlmdWb2Tsi8wWb2gpltCTwPCsw3M7vXzD40s7fM7LR41SUSC/XNfW9BgAJC\n0kM8WxAPA+d3m3cz8JJzbgLwUuA1wAXAhMDjWmBpHOsSiVpfO+oLUkBIOohbQDjnXgW63yp6EfBI\nYPoR4OKQ+b9yfuuAIjMrj1dtItHwOR/7Wvb1+SQ1fBIQ6o9JUlmiz0EMdc5VAwSehwTmDwd2hCxX\nFZgnknL2Nu/F53xdl6v2xZABQ2juaKapvSmGlYnEVqqcpLYw88J+tTKza81sg5lt2L17d5zLEjlS\n8NBQtAERui6RVJTogKgNHjoKPAf/OqqAkSHLjQB2hVuBc265c26ac25aaWlpXIsVCSf4oT50wNA+\nr0MBIekg0QHxNHBlYPpK4KmQ+VcErmaaATQGD0WJpJrag/47oGPRgrjmhjpqamJSlkjMxfMy198C\na4GJZlZlZt8C7gTOMbMtwDmB1wCrgErgQ+AXwHXxqkskWrE8xPTux9Xcep0SQlLT0fpiiopz7rIe\n3jorzLIOWBKvWkRiqa6pDo95+nwVk9cLLe1D4YfgBtSxdGUZSw1yc6G5OcbFikQhVU5Si6SNuqY6\nSvJKyPBk9OnnKythwSlboKkYCnaSRxMLp77D1q0xLlQkSgoIkV6qO1QX1Qnq8nIoGFsCB4bjGbiD\nFnIpGFdCWVkMixSJgbgdYhI5XtU11UV1/gGg1soYlVdMfskmZs/fTbVTOkjqUQtCpJdqD9ZGHRAr\nVsDZZ45lX1Eb9z1ZxooVMSpOJIYUECK9FIsWBMDwguHUNtXS4euIQVUisaeAEOmFyh3NHGg7gNcX\nfUAMGzgMn/PpZjlJWQoIkV748U2bAFi3IjfqdQ0bOAyAnft3Rr0ukXjQSWqRCHi90NICDHMwCda8\nPh6L8t6FYEDsOhC2VxmRpFMLQiQClZWwYOo7ZA+sBCDnwKCo711QQEiqU0CIRCB470J7gb8H4bb9\nY6O+d2HogKF4zKOAkJSlgBCJUK2V8alpW8hwHr59rlFDdPcuZHgyKMsvo+pAVYwqFIktnYMQidCK\nFXDVU/upqxzG0idjM57VqMJRbG/cHpN1icSaWhAivbBz/06GD4zdYIejC0fzccPHMVufSCwpIER6\nYeeBnV0nl2NhTNEYtjdux+d8MVunSKwoIER6IR4tiHZfO9UHND6WpB4FhEiEmtqaaGxtZHhB7AJi\nTNEYAD5u1GEmST0KCJEIBS9HjWkLomg0gM5DSEpSQIhEaOcBf5cYsTwHMbowEBBqQUgKUkCIRCjY\nZ1IsDzENyB5ASV4J2xq2xWydIrGigBCJ0I79OwAYUTAipusdXThaLQhJSQoIkQhta9hGSV4J+dn5\nMV3vmKIxbN2nAakl9SggRCK0rWFb11VHsTRh8AQq91Vq4CBJOQoIkQhtbdjK2KKxMV/vhOIJtPva\n+ewF26mpifnqRfpMASESAZ/z8XHDx3FpQZxYfCIAGyrf49brlBCSOhQQIhGoOVhDa2drzFsQXi+c\nMckfEK74Q5auLMPMP18k2RQQIhEIXoYa6xZEZSVcNr4WWgpg8BbyaIp6ICKRWFFAiEQgXgFRXg6F\nY0uh/kQ8xe/RQm7UAxGJxIrGgxCJQPAy1Hicg6i1MibkjaSp7DUunr+baqd0kNSgFoRIBLY2bGXo\ngKF4s2J/cmDFClhw8RSqM+r5zyeKWLEi5psQ6RMFhEgEPqj/gAnFE+K2/smlk3E4Nu/eHLdtiPSW\nAkLkGKqrYd2H7zMqb2LctnHKkFMAeLvu7bhtQ6S3FBAix/Dv332P9uw6tr4cu076uptQPIGcjBze\nrlVASOrQSWqRHni90NICjGiAU2DtmmmYQW4uNDfHdluZnkxOLj1ZLQhJKWpBiPSgshIWTH2H7OK3\nAMitHxnXexROHXKqAkJSigJCpAfl5VAwtoS2kkrozKR130lxvUfh1CGnsuvALvY2743PBkR6SQEh\nchS1Vsa4KW8zhhIWX9xADfG7R2HK0CkA/KPmH3HbhkhvJCUgzGybmb1tZhvNbENg3mAze8HMtgSe\nByWjNpFQK1ZA7vhtTJk0nfueLIvrPQqnDzsdgA27NsRvIyK9kMwWxBeccxXOuWmB1zcDLznnJgAv\nBV6LJFVLRwvv73mfKUOmxH1bJXkljC0ay/pd6+O+LZFIpNIhpouARwLTjwAXJ7EWEQA21W2i03VS\nUVaRkO1NHz6dv+38W0K2JXIsyQoIB/zJzN4ws2sD84Y656oBAs9DklSbSJeNNRsBmFo2NSHb+/Sw\nT/Nx48fMPLtWgwdJ0iUrIGY5504DLgCWmNmZkf6gmV1rZhvMbMPu3bvjV6EI8I/af5Cfnc+4QeMS\nsr3pw6cD8Ncdf9XgQZJ0SQkI59yuwHMdsBKYDtSaWTlA4Lmuh59d7pyb5pybVlpamqiSpZ/aWLOR\nKUOn4LH4/6l4vXDmhNPAl4EbuU6DB0nSJTwgzGyAmQ0MTgPnAu8ATwNXBha7Engq0bWJhHLO8Vbt\nW0wdmpjDS5WVsGDSVjy7ToPRr2jwIEm6ZHS1MRRYaWbB7f/GOfecmf0N+J2ZfQvYDnwtCbWJdNmy\ndwuNrY2cVn5aQrYXvDHPt202zLyH5qxODR4kSZXwgHDOVQJHfCVzztUDZyW6HpGerN2xFoCZI2Ym\nbJu1VsaFw6fyx4x2vnTZc9Qc+HrCti3SXSpd5iqSUtZWraUwp5CTS09O2DZXrIDf/PeXybAMplz5\nlgYPkqRSQIiEUV0Nj65ZS0XpZxJygjpUQU4Bnx7+aV6ofCGh2xXpTgEhEsYPr99CU/7bNG04NSnb\nnzN+Dn/b+TfqmsJezCeSEAoIkRBeL5jBgxu3gjk2vHBuUi41vfDEC3E4Vm9ZndgNi4RQQIiECI4B\nkTluNXRm4d1ekZRLTT9V9inK88t5dsuzid2wSAgFhEiI4KWmHeNewXbMoLW9OCmXmpoZXzrxS6ze\n8hyf+0Kzut2QpFBAiHSzPTMTyv/OdUOKWTR/d1zHgDiaSyZfQlP7Qf5S96y63ZCkUECIdHPFj18C\n4BvX3xz3MSB64vXC2eNnw8GhuMm/U7cbkhQKCJFunvngGQZ7B3cN4JMMlZWwYMpmMjddDCc+izdn\nl7rdkIRTQIiEaOts49kPnmXuxLlkepLRE41f8FxI5z++CVkttJy6Ut1uSMIpIERCrNm2hsbWRuaf\nND/ZpVBrZSyaMZqJHSMZfNZSqhma7JKkn1FAiARUV8OVdz1JXuYAzjnhnGSXw4oVcP+T5Xzvon+j\n3ruJf/npumSXJP2MAkIk4Mff2UbN4N8xfPe55GbmJrucLpdPuZyi3CL+35fv5vOfR5e8SsIoIKTf\nC949/Yv314O3gS3PXJdSVwzlZ+ez5NNLWLV1Ja+9u1mXvErCKCCk3wvePe2peAAaR+LdOj2lrhjy\neuH2uddDRw7uzNt1yaskjAJC+r3ycvBN3I9v/Itk/P1yWt2AlLpiqLISFozfTea678CUx8gtfz2l\nAkyOXwoIEWD9kN/gcVk8P6UkqXdPh9N1yevrt0BTCS3n/jsDxxWnTIDJ8UsBIf1adTV89qx6asoe\n5orTFnDWT7+ftLunj6bWylh8YRs3u3Nh7Cusy/v/dcJa4i55dwKJpIDbltSwNuO/oP0Q//zZf052\nOT3yB1YZ7Z0P88KD77OpfRGtfziDW69z3L9CTQmJD7UgpF8KXrm09AVg+s/hrYVMHjIp5U/8FuRn\n8ca/P0KLHcBduISlK4fqhLXEjQJC+qXglUsZZ/8TeDrIfeWf0+LEb2UlLBjmyPrzD2HSk2RN/0/m\nT3qPigodbpLYU0BIv1ReDk2nfERnxW/I+MsNtO2dnFJXLvWka7yKv/wrnve+RPv5N/NG3mbW/9Wn\n+yMk5hQQ0u9UV8OsL+7nTyO/z8C2Ebw+aWzKXbl0NLVWxuJ59WQ89SvYcxIfn3U1vrKNuj9CYk4n\nqaXfuW1JDX8p/BcsZxuvXf0qM0bNYkayi+qF4Anrf6+GxfP/i6dnXIP7xjnkPrKKOcWF7Co4iZoa\nUr41JKlPLQjpN7pOTH/8DFT8GvfKD/nc6Flp+427vBzKy6bAIy9CRy4t3zyfN3Pf1+EmiRkFhPQL\n1dUwdSrM/OJy+NJi2HI+3ldvTIsT00dTa2Us/kIBWQ+/AgfL2Xbh1/Cd9kDX1U0ej05eS98pIKRf\nuG1JDetrX2bdzO9D7RRy/ue/OASDAAAOvUlEQVQRWn35aXFi+mhWrID7nizj43fG89W//gLPtjPg\ny98m4+vzGDfkPcCpNSF9poCQ41Z1NWRkBA4rbfob7vI5uH2j4bFVuM7BaXVi+ljKy6FkxAm4R5+D\nP91F58Q/Unn5WbjJj6s1IX2mgJDj1m1LavC5DorP/je49GKoPRXvQ6tZeMIePq7KTMkuNaIRvLrp\nhZNGMPLxP8LBMvjqAjxXfY4RE58l2JqorkbddEhEFBByXDms1fBSByyYS/3n7oBNX4dHXqa1eXja\nH1bqSfBw09kPLuDCyVOwX6yFp/8bX8kWqi77Mu4bZ7P0H+8wbLiPV19VWMixKSDkuHLbkhp81kLp\nmbfDdZNhzBps1T2Mf/G/ePHyZ46rw0pH429N7OWFinxO+OU67Pn/D4ZsgivOge9MhFn/h6V/MoYN\nQ2EhPTLnXLJr6LNp06a5DRs2JLsMSbLqahgxAny0w+TfwRd+BIMrofIseGYZnn1j+fa83f22U7vF\n82pY/odSMjP303byszBtOYx+HXwe+PhM2PwVeO8i2D8ScCyeV8sP7yvj0kvhiSfAObqmj8eWV39k\nZm8456Ydazm1ICRthH7DDZ2+5Ya38c28k8wbx8BXLoe2fOzRZxn/1K95cd76ftNq6EmtlbFo/m7W\nL17B5D0XwkOvkvXzN+G1f4P8GphzPXx/FHznJJhzPUs3r2XYhLqulsVtS2p4/TXfEa0MtTiOf2pB\nSEqrrv7k2+ut19Xw338YwrcvrqM5o56HP1wPpzwO414Ejw8qvwjrvgdb5uBx9OtWQ0/mz4dyq+Ha\n4X9k4e/nsal6EFmlf6d9/MvYuJdwo1+F7EP+hRtGw87p/kdNBdRN9p/4BhbPqwXo+n2EtjjKyg7/\nvakFknoibUEoICRpevoQCZ2+9boalq4cCjkHYNgGGPsSjH8eyt8Ec7BvLLy9AN5egGf3iZxQuIel\nX3mZFQ1fpNodX1cpxVq4sMjO2E/b8I14RqzDN+xNGL4eBm375Ieai2D3JNg9GfaOh4Yx/se+sXCo\nmMXz6rh/RRnXzf8kzCF8kPT0O+9pWuESOwoISbhIPvC7f/h3/xC59uJamjJr+PWGeijd7A+F4X+D\nkvf8gdCZCVUzsI/OwX14AZ5dU/GRRQ6ttJOlVkMfhQuLXFpoIQcG7CZryJu0l27BSjfhSt+D0ndh\nwJ7DV9KW5w+LA8P9LY1wj6ZSaCli8UX1AGFDpKfpaMIlmunjMZjSNiDM7Hzgp0AG8IBz7s6ellVA\nHKm3H9KR/oFE+m3/WH/kDsey5zOgsAoKdkBBFRTugMLtUPwBFL8POQc/+QcdHAo7p2M7T8Pt+gzZ\nO06nrbUEsK4PsMnl+3jsq39g+c4L1WqIgdCwmPfkQgBWfuWxI4Mj5wBZhR/QXrQLT9FWfEU78Az6\nkKyCnbTl7cENrIHM1vAbaR0IzYOgZRA0D/ZPNw/2v24pgrb8Tx6tA0NeD+CKs5rJ9uXz4IrRLLrY\nH1K9CZreTscqmCL9O0pEgKVlQJhZBvABcA5QBfwNuMw592645fsaEEc7PhqPX2AipyP5kO7tdPdD\nBj46Wb7ay+Vzt9LuaeLx13Ig+6D/MFD2AfDuhbw9kFcfeN4D3pDpzLbDfyGdWbB/BOwdj+05EfZM\nxO2ZjGfPCfgOjCSHNlrJ7gqC0A8thULiRBIcbWRzUnkj71UXkkUrrbktkF9LVv522vP34Mnbjc+7\nH/Pugdx9OG+j//9GboP//413X8+hEk7bAGj3QkduBI8wy3VmgS8LOrNDpkOeO7MD05lc/LkDeFwW\nK18ewtdmH8TjMnni+TIuP6cRj8viV38s5+oLGjAyePCpMv7Xl+vxkMHyP5SF/TuC2P199la6BsRM\n4H87584LvL4FwDl3R7jl+xoQPf2S4vELvPbiWsDx30+V8L8uqgXz8YunS/nW3Bowx4PPlnDVl/zz\nH/pjCVde6F/+kecGc/kFtTjr4LEXi7jk3FocnfxuTQFf+WIdPutk5WsDmXvmbnzWyTN/GQieTvB0\nHP2R0QaZLYFHq/85ozWCec2fhEBW87F3sjP/N8JDxXCoBA6VYof8r92BkXgay/HtH0P2/lLamkaA\ny/jkm6laB2klNDiW77yQlWvLmDczghZIyO/5sOnMVnKy6mnNbofsJrKy62nPbsOTvR9fTjOW3QjZ\nTbjsQ1j2fshswWW2YZlNuMzWT/6/hv7fzWoJ+f8ceFgCP/s6M8FlgC+j52dfL5d55zJ48xoAcnOh\nOYI/y6B0DYivAuc7564JvP4G8Bnn3HfCLd/bgPB6oaUFOOF5OP9GMB/g/M/m8/+H6WneYfPDzQvz\n8x5flHskzjpyAo9c6MwJ+WaVHZiXi6cji1znobU1n862IqxtALQOxLUVYq150Oaf9rR68bUVkdWa\nQ3vzUGgZRK5rP+oHQbgWgVoHx6+eWiCRTEccLhFO59BMa4ZBRgc5ngNd09meg7RlGHg6yMo4SHsG\nZHha6cxweDzNkNGBz+PDk3EIl9GB8/iwjBbwtOMyfJinBTw+nDnM0wbmw3kcmZ5m8rKaOeTLoMM8\nmKcdrBPncZi1g6fTvy5rD/y8Dzztn3yOWEfgC2DnJ9NvXU7ehiuYN3Urdz93Sq8ONaVrQHwNOK9b\nQEx3zl0fssy1wLUAo0aNOv3jjz+OeP3V1XDTBe/wZH0drTOXYc6By8C5DDLoJNfTSUtHNp0uO/Dl\nwnAu0z/tAtM4cJ7A/JBpfF3r6lovGYdtw5wDwvysc4FtZWHOF5ifhfkMfJk4X04gg7Lw+XL8udPD\ntM+XQ5bPR7vPC74ssn2dtPnywJdFjq+DVl8+dHjJ6fTR6rwc6w8q9JBBNm1RfeDrw1/6KppwiWY6\nVsHUm7+j3q63LxdmRBoQOOdS5gHMBJ4PeX0LcEtPy59++umutxZdXO08dLhcDjnodOBzuRxyHjrc\npPL6sO+l0/Tk8j1u4/UPuLHDmt3YYc1RT183v9qVlzt33fzqYy43b16vfx0iKW3evGP/34/l31Ff\n1tuXvztgg4vkMzmShRL1wD8EaiUwFsgG/gFM7mn5vgRET7/weP0CEzmtD2kRiUSkAZFSh5gAzGwO\ncA/+y1x/6Zy7vadldZmriEjvRXqIKTMRxfSGc24VsCrZdYiI9HfqrE9ERMJSQIiISFgKCBERCUsB\nISIiYSkgREQkrJS7zLU3zGw3EPmt1IcrAfYcc6nES9W6IHVrU129o7p653isa7RzrvRYC6V1QETD\nzDZEch1woqVqXZC6tamu3lFdvdOf69IhJhERCUsBISIiYfXngFie7AJ6kKp1QerWprp6R3X1Tr+t\nq9+egxARkaPrzy0IERE5iuM6IMws18zWm9k/zGyTmf0/YZbJMbMnzOxDM/urmY1Jkbq+aWa7zWxj\n4HFNvOsK2XaGmf3dzJ4N817C91eEdSVlf5nZNjN7O7DNI7oWNr97A/vrLTM7LUXqmm1mjSH760eJ\nqCuw7SIz+72ZvWdmmwNDDYe+n/B9FkFNSdlfZjYxZJsbzWy/mX2v2zJx218p15trjLUCX3TOHTSz\nLOB1M1vtnFsXssy3gH3OufFmdilwF3BJCtQF8ITrYbjVOLsB2AwUhHkvGfsrkrogefvrC865nq5H\nvwCYEHh8BlgaeE52XQCvOee+lKBaQv0UeM4591Uzywbyur2fjH12rJogCfvLOfc+UAH+L0jATmBl\nt8Xitr+O6xZEYGyMg4GXWYFH95MuFwGPBKZ/D5xlZpYCdSWFmY0ALgQe6GGRhO+vCOtKVRcBvwr8\nztcBRWZWnuyiksXMCoAzgQcBnHNtzrmGbosldJ9FWFMqOAv4yDnX/ebguO2v4zogoOuwxEagDnjB\nOffXbosMB3YAOOc6gEagOAXqAvhKoMn4ezMbGe+aAu4B/gXw9fB+UvZXBHVBcvaXA/5kZm+Yf7z0\n7rr2V0BVYF6y6wKYGTjMudrMJiegJoBxwG7gocDhwgfMbEC3ZRK9zyKpCZKzv0JdCvw2zPy47a/j\nPiCcc53OuQpgBDDdzE7ptki4b79x/zYfQV3PAGOcc1OAF/nkW3vcmNmXgDrn3BtHWyzMvLjurwjr\nSvj+CpjlnDsNfzN/iZmd2e39pPz/4th1vYm/u4WpwM+APySgJvAf1j4NWOqc+xTQBNzcbZlE77NI\nakrW/gIgcNhrLvA/4d4OMy8m++u4D4igQJNxDXB+t7eqgJEAZpYJFAJ7k12Xc67eOdcaePkL4PQE\nlDMLmGtm24DHgS+a2aPdlknG/jpmXUnaXzjndgWe6/AfG57ebZGu/RUwAtiV7Lqcc/uDhzkDozhm\nmVlJvOvCvz+qQlrMv8f/4dx9mUTus2PWlMT9FXQB8KZzrjbMe3HbX8d1QJhZqZkVBaa9wNnAe90W\nexq4MjD9VeBlF+ebQyKpq9sxxLn4T87GlXPuFufcCOfcGPzN2Zedc5d3Wyzh+yuSupKxv8xsgJkN\nDE4D5wLvdFvsaeCKwJUmM4BG51x1susys7LguSMzm47/s6A+nnUBOOdqgB1mNjEw6yzg3W6LJXSf\nRVJTsvZXiMsIf3gJ4ri/jvermMqBRwJn/z3A75xzz5rZrcAG59zT+E9M/drMPsT/TfjSFKnru2Y2\nF+gI1PXNBNQVVgrsr0jqSsb+GgqsDHxuZAK/cc49Z2aLAJxzy/CPrz4H+BA4BFyVInV9FVhsZh1A\nM3BpvIM+xPXAY4HDJpXAVSmwz45VU9L2l5nlAecA3w6Zl5D9pTupRUQkrOP6EJOIiPSdAkJERMJS\nQIiISFgKCBERCUsBISIiYSkgREQkLAWEiIiEpYAQibHAncx/DHTs9o6ZJao7dJGYOt7vpBZJhvOB\nXc65CwHMrDDJ9Yj0iVoQIrH3NnC2md1lZmc45xqTXZBIXyggRGLMOfcB/t5k3wbusAQO5ykSSzrE\nJBJjZjYM2Ouce9TMDpLEjhZFoqGAEIm9U4H/Y2Y+oB1YnOR6RPpEvbmKiEhYOgchIiJhKSBERCQs\nBYSIiISlgBARkbAUECIiEpYCQkREwlJAiIhIWAoIEREJ6/8CqXl3EUhV4rcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot\n",
    "plt.plot(s, f[0], 'r.', label='f[0]')\n",
    "plt.plot(s, f_fit[0], 'b*', label='f_fit[0]')\n",
    "plt.plot(s_test, f_test[0], 'g', label='f_test[0]')\n",
    "plt.xlabel('s')\n",
    "plt.ylabel('f')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Large problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "Ns = 50000\n",
    "N = 1000\n",
    "s = np.linspace(1.e-2, 5.e3, Ns)\n",
    "poles = np.linspace(1.1e-2, 4.8e+3, N/2);\n",
    "poles = poles + poles*0.01j\n",
    "poles = np.sort(np.append(poles, np.conj(poles)))\n",
    "residues = np.linspace(1e+2, 1e+6, N/2);\n",
    "residues = residues + residues*0.5j\n",
    "residues = np.sort(np.append(residues, np.conj(residues)))\n",
    "residues = residues.reshape((1, N))\n",
    "f = np.zeros((1, Ns))\n",
    "for p, r in zip(poles, residues[0]):\n",
    "    f[0] += (r/(s - p)).real\n",
    "weight = 1.0/f\n",
    "\n",
    "poles_init = np.linspace(1.2e-2, 4.7e+3, N/2);\n",
    "poles_init = poles_init + poles_init*0.01j\n",
    "poles_init = np.sort(np.append(poles_init, np.conj(poles_init)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 4min 14s, sys: 6min 42s, total: 10min 57s\n",
      "Wall time: 2min 51s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.005450713259715198"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time poles_fit, residues_fit, cf, f_fit, rms = m.vectfit(f, s, poles_init, weight)\n",
    "rms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.testing.assert_allclose(f, f_fit, 1e-2)"
   ]
  }
 ],
 "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.5.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}