{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Using matplotlib backend: TkAgg\n",
        "Populating the interactive namespace from numpy and matplotlib\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.special import i0, j0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy import optimize"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "g0 = lambda b: i0(b) * exp(-b)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def get_ac(ac):\n",
      "    a = ac[:ac.shape[0]/2]\n",
      "    c = ac[ac.shape[0]/2:]\n",
      "    return a, c "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def obj_f_2(params):\n",
      "    a, c = get_ac(params)\n",
      "    x = linspace(0, 20, 2000)\n",
      "    s = sum([ci * j0(ai * x) ** 2 for (ci, ai) in zip(c, a)])\n",
      "    return sqrt((g0(x ** 2) - s) ** 2).sum() * (x[1] - x[0])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 6
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def obj_f_inf(params):\n",
      "    a, c = get_ac(params)\n",
      "    x = linspace(0, 20, 2000)\n",
      "    s = sum([ci * j0(ai * x) ** 2 for (ci, ai) in zip(c, a)])\n",
      "    return ((g0(x ** 2) - s) ** 2).max()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cons = (\n",
      "        {'type': 'eq', # sum c_i = 1\n",
      "         'fun': lambda params: params[:params.shape[0] / 2].sum() - 1.},\n",
      "        {'type': 'eq', # sum (c_i * a_i) ** 2 = 2\n",
      "         'fun': lambda params: (params[:params.shape[0] / 2] * params[params.shape[0] / 2:] ** 2).sum() - 2.},\n",
      "        {'type': 'ineq', # radii > 0\n",
      "         'fun': lambda params: params[params.shape[0] / 2:].min()},\n",
      "        {'type': 'ineq', # radii <= 2.5\n",
      "         'fun': lambda params: 2.5 - params[params.shape[0] / 2:].max()}\n",
      "        )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "case_arraygen = lambda n: hstack(((1 + arange(n) * 2.5 / n), ones(n) / float(n)))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hstack(((1 + arange(4) * 2.5 / 4, ones(4) / 4.))) == case_arraygen(4)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 10,
       "text": [
        "array([ True,  True,  True,  True,  True,  True,  True,  True], dtype=bool)"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "case_arraygen(20)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "array([ 1.   ,  1.125,  1.25 ,  1.375,  1.5  ,  1.625,  1.75 ,  1.875,\n",
        "        2.   ,  2.125,  2.25 ,  2.375,  2.5  ,  2.625,  2.75 ,  2.875,\n",
        "        3.   ,  3.125,  3.25 ,  3.375,  0.05 ,  0.05 ,  0.05 ,  0.05 ,\n",
        "        0.05 ,  0.05 ,  0.05 ,  0.05 ,  0.05 ,  0.05 ,  0.05 ,  0.05 ,\n",
        "        0.05 ,  0.05 ,  0.05 ,  0.05 ,  0.05 ,  0.05 ,  0.05 ,  0.05 ])"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cases = [(case_arraygen(4), obj_f_2), \n",
      "         (case_arraygen(5), obj_f_2), \n",
      "         (case_arraygen(6), obj_f_2), \n",
      "         (case_arraygen(7), obj_f_2), \n",
      "         (case_arraygen(10), obj_f_2)]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i, (params0, obj_f) in enumerate(cases):\n",
      "    print i, params0, obj_f"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0 [ 1.     1.625  2.25   2.875  0.25   0.25   0.25   0.25 ] <function obj_f_2 at 0x04A6DF30>\n",
        "1 [ 1.   1.5  2.   2.5  3.   0.2  0.2  0.2  0.2  0.2] <function obj_f_2 at 0x04A6DF30>\n",
        "2 [ 1.          1.41666667  1.83333333  2.25        2.66666667  3.08333333\n",
        "  0.16666667  0.16666667  0.16666667  0.16666667  0.16666667  0.16666667] <function obj_f_2 at 0x04A6DF30>\n",
        "3 [ 1.          1.35714286  1.71428571  2.07142857  2.42857143  2.78571429\n",
        "  3.14285714  0.14285714  0.14285714  0.14285714  0.14285714  0.14285714\n",
        "  0.14285714  0.14285714] <function obj_f_2 at 0x04A6DF30>\n",
        "4 [ 1.    1.25  1.5   1.75  2.    2.25  2.5   2.75  3.    3.25  0.1   0.1\n",
        "  0.1   0.1   0.1   0.1   0.1   0.1   0.1   0.1 ] <function obj_f_2 at 0x04A6DF30>\n"
       ]
      }
     ],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%time\n",
      "results = []\n",
      "for i, (params0, obj_f) in enumerate(cases):\n",
      "    res = optimize.minimize(obj_f, params0, \n",
      "                            method='SLSQP',\n",
      "                            options={'maxiter': 20},\n",
      "                            constraints=cons)\n",
      "    results.append(res)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Wall time: 18.8 s\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "for i, res in enumerate(results):\n",
      "    print(\"---- %i points\" % int(len(res.x / 2)))\n",
      "    print(res)\n",
      "    for j in range(int(len(res.x)/2)):\n",
      "        print(\"lrCirclesDouble(%i, 1) = %f # coeff circle %i\" % (j + 1, res.x[j + len(res.x) / 2], j + 1))\n",
      "        print(\"lrCirclesDouble(%i, 2) = %f # radius circle %i\" % (j + 1, res.x[j], j + 1))\n",
      "    print()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "---- 8 points\n",
        "  status: 9\n",
        " success: False\n",
        "    njev: 21\n",
        "    nfev: 270\n",
        "     fun: 1.5142105724341708\n",
        "       x: array([ 0.0448529 , -0.03021077,  2.97214753,  2.42150225, -0.01274741,\n",
        "       -0.01496447,  0.7140736 , -0.01418282])\n",
        " message: 'Iteration limit exceeded'\n",
        "     jac: array([  7.80102554e+02,  -6.82055811e+02,  -9.77405781e+01,\n",
        "         2.39090502e+00,   1.07795319e+04,   1.15683884e+04,\n",
        "         4.45916596e+02,   5.31407425e+02,   0.00000000e+00])\n",
        "     nit: 21\n",
        "lrCirclesDouble(1, 1) = -0.012747 # coeff circle 1\n",
        "lrCirclesDouble(1, 2) = 0.044853 # radius circle 1\n",
        "lrCirclesDouble(2, 1) = -0.014964 # coeff circle 2\n",
        "lrCirclesDouble(2, 2) = -0.030211 # radius circle 2\n",
        "lrCirclesDouble(3, 1) = 0.714074 # coeff circle 3\n",
        "lrCirclesDouble(3, 2) = 2.972148 # radius circle 3\n",
        "lrCirclesDouble(4, 1) = -0.014183 # coeff circle 4\n",
        "lrCirclesDouble(4, 2) = 2.421502 # radius circle 4\n",
        "()\n",
        "---- 10 points\n",
        "  status: 9\n",
        " success: False\n",
        "    njev: 21\n",
        "    nfev: 319\n",
        "     fun: 893.72527154637839\n",
        "       x: array([ 0.9408399 ,  0.62182301,  1.08429362,  0.87432457,  1.69636943,\n",
        "       -0.10305518, -0.6269415 , -0.54183445,  1.14040183,  0.77232284])\n",
        " message: 'Iteration limit exceeded'\n",
        "     jac: array([  343.53822327,  4552.79451752,  1551.80874634, -4695.21598816,\n",
        "       -1046.68219757,  3775.85966492,  5272.16714478,  3365.53674316,\n",
        "        4047.48488617,  2328.04063416,     0.        ])\n",
        "     nit: 21\n",
        "lrCirclesDouble(1, 1) = -0.103055 # coeff circle 1\n",
        "lrCirclesDouble(1, 2) = 0.940840 # radius circle 1\n",
        "lrCirclesDouble(2, 1) = -0.626942 # coeff circle 2\n",
        "lrCirclesDouble(2, 2) = 0.621823 # radius circle 2\n",
        "lrCirclesDouble(3, 1) = -0.541834 # coeff circle 3\n",
        "lrCirclesDouble(3, 2) = 1.084294 # radius circle 3\n",
        "lrCirclesDouble(4, 1) = 1.140402 # coeff circle 4\n",
        "lrCirclesDouble(4, 2) = 0.874325 # radius circle 4\n",
        "lrCirclesDouble(5, 1) = 0.772323 # coeff circle 5\n",
        "lrCirclesDouble(5, 2) = 1.696369 # radius circle 5\n",
        "()\n",
        "---- 12 points\n",
        "  status: 9\n",
        " success: False\n",
        "    njev: 21\n",
        "    nfev: 364\n",
        "     fun: 146.79237066627917\n",
        "       x: array([ 1.81339974, -0.62011745,  1.54525265,  1.59164722,  2.4200528 ,\n",
        "        1.84299224,  0.2280895 ,  0.09564532, -0.50406859, -0.55498723,\n",
        "        0.44197458,  0.46131993])\n",
        " message: 'Iteration limit exceeded'\n",
        "     jac: array([ -244.97188187,   708.38697433,   796.7755909 ,   608.20923996,\n",
        "        -240.55032921,  -541.35148621,  2208.0682106 ,  5284.67614555,\n",
        "        2514.69530296,  2452.59099197,  1731.76121902,  2174.29994011,\n",
        "           0.        ])\n",
        "     nit: 21\n",
        "lrCirclesDouble(1, 1) = 0.228090 # coeff circle 1\n",
        "lrCirclesDouble(1, 2) = 1.813400 # radius circle 1\n",
        "lrCirclesDouble(2, 1) = 0.095645 # coeff circle 2\n",
        "lrCirclesDouble(2, 2) = -0.620117 # radius circle 2\n",
        "lrCirclesDouble(3, 1) = -0.504069 # coeff circle 3\n",
        "lrCirclesDouble(3, 2) = 1.545253 # radius circle 3\n",
        "lrCirclesDouble(4, 1) = -0.554987 # coeff circle 4\n",
        "lrCirclesDouble(4, 2) = 1.591647 # radius circle 4\n",
        "lrCirclesDouble(5, 1) = 0.441975 # coeff circle 5\n",
        "lrCirclesDouble(5, 2) = 2.420053 # radius circle 5\n",
        "lrCirclesDouble(6, 1) = 0.461320 # coeff circle 6\n",
        "lrCirclesDouble(6, 2) = 1.842992 # radius circle 6\n",
        "()\n",
        "---- 14 points\n",
        "  status: 9\n",
        " success: False\n",
        "    njev: 21\n",
        "    nfev: 411\n",
        "     fun: 69.090271420899128\n",
        "       x: array([ 0.29620416,  0.93421717,  1.22525633,  1.38202028,  3.11537323,\n",
        "        2.50110766,  2.07151724, -0.02545283, -0.05229518, -0.28703569,\n",
        "        0.09626451,  0.56893896, -0.26608561,  0.32169302])\n",
        " message: 'Iteration limit exceeded'\n",
        "     jac: array([ -751.64078999,  -186.46660233,  -705.96837616,   190.08371258,\n",
        "         247.71982193,  -162.75920868,   212.2142725 , -9428.33320427,\n",
        "       -3798.70682812, -3049.51406384, -2760.56082439, -1396.47885227,\n",
        "       -1679.88785172, -1969.59872818,     0.        ])\n",
        "     nit: 21\n",
        "lrCirclesDouble(1, 1) = -0.025453 # coeff circle 1\n",
        "lrCirclesDouble(1, 2) = 0.296204 # radius circle 1\n",
        "lrCirclesDouble(2, 1) = -0.052295 # coeff circle 2\n",
        "lrCirclesDouble(2, 2) = 0.934217 # radius circle 2\n",
        "lrCirclesDouble(3, 1) = -0.287036 # coeff circle 3\n",
        "lrCirclesDouble(3, 2) = 1.225256 # radius circle 3\n",
        "lrCirclesDouble(4, 1) = 0.096265 # coeff circle 4\n",
        "lrCirclesDouble(4, 2) = 1.382020 # radius circle 4\n",
        "lrCirclesDouble(5, 1) = 0.568939 # coeff circle 5\n",
        "lrCirclesDouble(5, 2) = 3.115373 # radius circle 5\n",
        "lrCirclesDouble(6, 1) = -0.266086 # coeff circle 6\n",
        "lrCirclesDouble(6, 2) = 2.501108 # radius circle 6\n",
        "lrCirclesDouble(7, 1) = 0.321693 # coeff circle 7\n",
        "lrCirclesDouble(7, 2) = 2.071517 # radius circle 7\n",
        "()\n",
        "---- 20 points\n",
        "  status: 9\n",
        " success: False\n",
        "    njev: 21\n",
        "    nfev: 529\n",
        "     fun: 1.7766701053627527\n",
        "       x: array([ 0.67584164,  1.46998201,  1.70290549,  1.56124607,  1.70592654,\n",
        "        2.12681978,  2.27909545,  2.44647647,  2.70980765,  3.0313936 ,\n",
        "       -0.09916783, -0.14045072, -0.14205628,  0.08612575, -0.12090092,\n",
        "        0.04471519, -0.01746348,  0.22541161,  0.25204293,  0.29373551])\n",
        " message: 'Iteration limit exceeded'\n",
        "     jac: array([ -459.29809998,  -176.89344704,  -187.51368636,   122.5751541 ,\n",
        "        -157.1198234 ,    35.75088869,   -11.90532859,   144.48435917,\n",
        "         101.61165623,    98.05111745, -4975.51527296, -2629.48916389,\n",
        "       -2319.29090156, -2490.59279244, -2315.33353958, -1931.3035195 ,\n",
        "       -1822.12778057, -1716.09995891, -1572.23719129, -1430.56669591,\n",
        "           0.        ])\n",
        "     nit: 21\n",
        "lrCirclesDouble(1, 1) = -0.099168 # coeff circle 1\n",
        "lrCirclesDouble(1, 2) = 0.675842 # radius circle 1\n",
        "lrCirclesDouble(2, 1) = -0.140451 # coeff circle 2\n",
        "lrCirclesDouble(2, 2) = 1.469982 # radius circle 2\n",
        "lrCirclesDouble(3, 1) = -0.142056 # coeff circle 3\n",
        "lrCirclesDouble(3, 2) = 1.702905 # radius circle 3\n",
        "lrCirclesDouble(4, 1) = 0.086126 # coeff circle 4\n",
        "lrCirclesDouble(4, 2) = 1.561246 # radius circle 4\n",
        "lrCirclesDouble(5, 1) = -0.120901 # coeff circle 5\n",
        "lrCirclesDouble(5, 2) = 1.705927 # radius circle 5\n",
        "lrCirclesDouble(6, 1) = 0.044715 # coeff circle 6\n",
        "lrCirclesDouble(6, 2) = 2.126820 # radius circle 6\n",
        "lrCirclesDouble(7, 1) = -0.017463 # coeff circle 7\n",
        "lrCirclesDouble(7, 2) = 2.279095 # radius circle 7\n",
        "lrCirclesDouble(8, 1) = 0.225412 # coeff circle 8\n",
        "lrCirclesDouble(8, 2) = 2.446476 # radius circle 8\n",
        "lrCirclesDouble(9, 1) = 0.252043 # coeff circle 9\n",
        "lrCirclesDouble(9, 2) = 2.709808 # radius circle 9\n",
        "lrCirclesDouble(10, 1) = 0.293736 # coeff circle 10\n",
        "lrCirclesDouble(10, 2) = 3.031394 # radius circle 10\n",
        "()\n"
       ]
      }
     ],
     "prompt_number": 15
    }
   ],
   "metadata": {}
  }
 ]
}