{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "import itertools\n",
    "\n",
    "import numpy\n",
    "from sympy import init_session"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IPython console for SymPy 1.0 (Python 2.7.12-64-bit) (ground types: python)\n",
      "\n",
      "These commands were executed:\n",
      ">>> from __future__ import division\n",
      ">>> from sympy import *\n",
      ">>> x, y, z, t = symbols('x y z t')\n",
      ">>> k, m, n = symbols('k m n', integer=True)\n",
      ">>> f, g, h = symbols('f g h', cls=Function)\n",
      ">>> init_printing()\n",
      "\n",
      "Documentation can be found at http://docs.sympy.org/1.0/\n"
     ]
    }
   ],
   "source": [
    "init_session()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "C = symbols(\"c0:4\")\n",
    "S = symbols(\"s0:4\")\n",
    "A = symbols(\"a0:4\")\n",
    "B = symbols(\"b0:4\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def dot(es1, es2):\n",
    "    return sum([e[0]*e[1] for e in zip(es1, es2)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def do_integral(f, s):\n",
    "    return Integral(f, (s, 0, 1))\n",
    "        # The bounds of 0 and 1 are specific to the Fourier FA. Make\n",
    "        # sure you normalize the inputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "da_integral = reduce(do_integral, S[1:], cos(pi * dot(S, C))).doit()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "# Copied from hiora_cartpole.fourier_fa.\n",
    "order = 3\n",
    "n_dims = 4\n",
    "c_matrix = np.array(\n",
    "               list( itertools.product(range(order+1), repeat=n_dims) ),\n",
    "               dtype=np.int32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def sum_term(integral, c, c_vec):\n",
    "    return integral.subs(zip(c, c_vec))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "sum_terms = [sum_term(da_integral, C, c_vec) for c_vec in c_matrix]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "np_sum_terms = [lambdify(S[0], sum_term, 'numpy') for sum_term in sum_terms]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "def phi(npsts, theta, s0):\n",
    "    ns0 = (s0 - -2.5) / 5.0\n",
    "    return np.dot(theta, \n",
    "                  np.array([npst(ns0) for npst in npsts]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "theta = np.load(\"theta.npy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "res = np.array([phi(np_sum_terms, theta[512:768], x) \n",
    " for x in np.arange(-2.38, 2.5, 0.5*1.19)])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  -129838.63097278,   -298998.59748568,   3087412.84846354,\n",
       "         8130942.03887799,  10758471.51583669,   8664232.14626823,\n",
       "         3636646.02919975,   -222434.83286789,   -491482.63696056])"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "other = np.array([-3748598.374407076,\n",
    " -8333255.9176837215,\n",
    " 92906846.75614552,\n",
    " 242969379.49722022,\n",
    " 320543060.70953935,\n",
    " 257463642.38676718,\n",
    " 107526913.72252564,\n",
    " -7061727.3744605975,\n",
    " -14631018.954087665])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.03463658,  0.03588017,  0.03323127,  0.03346488,  0.03356326,\n",
       "        0.03365225,  0.0338208 ,  0.03149864,  0.03359183])"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res/other"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 28.87120995,  27.87055186,  30.09213582,  29.88207004,\n",
       "        29.79447966,  29.71569067,  29.56760511,  31.74739893,  29.76914718])"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "other/res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": [
    "res2 = np.array([phi(np_sum_terms, theta[768:1024], x) \n",
    " for x in np.arange(-2.38, 2.5, 0.5*1.19)])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  -174316.21632421,   -505178.62138594,   2743738.37912048,\n",
       "         7828381.28926773,  10698118.73072955,   8891660.84160946,\n",
       "         3996424.61865612,     49410.50215361,   -409740.34118969])"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "collapsed": false,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ -5.14232838e+06,  -1.49027693e+07,   8.09402822e+07,\n",
       "         2.30937248e+08,   3.15594503e+08,   2.62303995e+08,\n",
       "         1.17894526e+08,   1.45760981e+06,  -1.20873401e+07])"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res2*29.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is strange. Why are the results I get with numeric integration approximately 30 times those of the analytic integration?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "run_control": {
     "frozen": false,
     "read_only": false
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12+"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}