IRC_Exercise_3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate as integrate\n",
    "from scipy import linalg\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def weight(x):\n",
    "    '''\n",
    "    Define weight function\n",
    "    '''\n",
    "    ans = np.sin(np.log(x + 1))\n",
    "    return ans\n",
    "\n",
    "def integrand(x, i, j):\n",
    "    '''\n",
    "    Define integrand for the entries\n",
    "    of the correlation matrix\n",
    "    '''\n",
    "    p = i + j\n",
    "    ans = x**p*weight(x)\n",
    "    return ans\n",
    "\n",
    "def M_elem(i, j):\n",
    "    '''\n",
    "    Define the (i,j) element of the \n",
    "    correlation matrix M\n",
    "    '''\n",
    "    ans = integrate.quad(integrand, 0, 1, args = (i, j))[0]\n",
    "    return ans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def create_M(ord):\n",
    "    '''\n",
    "    Create correlation matrix M \n",
    "    '''\n",
    "    ans = np.array([M_elem(i, j) for i in range(ord + 1) for j in range(ord + 1)]).reshape(ord + 1, ord + 1)\n",
    "    return ans"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.36972237 0.23722455 0.17394447 0.13710581 0.1130647  0.09616204]\n",
      " [0.23722455 0.17394447 0.13710581 0.1130647  0.09616204 0.0836387 ]\n",
      " [0.17394447 0.13710581 0.1130647  0.09616204 0.0836387  0.07399217]\n",
      " [0.13710581 0.1130647  0.09616204 0.0836387  0.07399217 0.06633541]\n",
      " [0.1130647  0.09616204 0.0836387  0.07399217 0.06633541 0.06011148]\n",
      " [0.09616204 0.0836387  0.07399217 0.06633541 0.06011148 0.05495326]]\n"
     ]
    }
   ],
   "source": [
    "ord_pol = 5\n",
    "M = create_M(ord_pol)\n",
    "print M"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Until this point all seems fine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.64460699e+00 -4.35221622e+00  7.78361074e+00 -1.18172980e+01\n",
      "   1.63748194e+01 -2.14019465e+01]\n",
      " [ 0.00000000e+00  6.78307428e+00 -3.23507993e+01  9.17517877e+01\n",
      "  -2.02659391e+02  3.85031445e+02]\n",
      " [ 0.00000000e+00  0.00000000e+00  2.75756308e+01 -1.87731150e+02\n",
      "   7.24942769e+02 -2.09607563e+03]\n",
      " [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  1.11287328e+02\n",
      "  -9.82376045e+02  4.79138903e+03]\n",
      " [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   4.47631499e+02 -4.85198631e+03]\n",
      " [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  1.79733849e+03]]\n"
     ]
    }
   ],
   "source": [
    "# Now we use Cholesky decomp. \n",
    "L = linalg.cholesky(M)\n",
    "L_inv = linalg.inv(L)\n",
    "print L_inv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0000000000000018"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def pols(x):\n",
    "    '''\n",
    "    Create array of polynomials whose \n",
    "    coefficients are given by the elements of\n",
    "    L_inv\n",
    "    '''\n",
    "    x_vec = np.array([x**i for i in range(6)])\n",
    "    ans = L_inv.T.dot(x_vec)\n",
    "    return ans\n",
    "\n",
    "# Here I choose any two polynomials of the array (0 and 1, in this case) and \n",
    "# check orthogonality with the weight. It is not working. \n",
    "def check_pol_ort(x):\n",
    "    ans = pols(x)[1]*pols(x)[1]*weight(x)\n",
    "    return ans\n",
    "\n",
    "integrate.quad(check_pol_ort, 0, 1)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1]\n",
      " [2]\n",
      " [3]\n",
      " [4]]\n",
      "[[1 2 3 4]]\n"
     ]
    }
   ],
   "source": [
    "arr = np.array([[1],[2],[3],[4]])\n",
    "print arr\n",
    "print arr.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/Manu/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:1: DeprecationWarning: object of type <type 'float'> cannot be safely interpreted as an integer.\n",
      "  if __name__ == '__main__':\n"
     ]
    },
    {
     "ename": "NameError",
     "evalue": "name 'pol' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-30-8a1e259d9c64>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0mdom\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinspace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1e3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdom\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpol\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdom\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdom\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtest_pol\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdom\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'r--'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'pol' is not defined"
     ]
    }
   ],
   "source": [
    "dom = np.linspace(-100, 100, 1e3)\n",
    "plt.plot(dom, pol(dom)[0])\n",
    "plt.plot(dom, test_pol(dom), 'r--',)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def test_pol(x):\n",
    "    ans = 1.64 - 4.35*x + 7.78*x**2 - 11.8*x**3 + 16.3*x**4 - 21.4*x**5\n",
    "    return ans\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import integrate\n",
    "\n",
    "def intmono(n,m):\n",
    "    '''\n",
    "    Calculate the integral of \\int_0^1 x^n x^m sin(ln(x+10)) dx\n",
    "    Method is straight sum (not optimal but fast)\n",
    "    \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    n: int\n",
    "        exponent of x (first)\n",
    "    m: int\n",
    "        exponent of x (second)\n",
    "        \n",
    "    Returns\n",
    "    ----------\n",
    "    float\n",
    "        result of integration\n",
    "        \n",
    "    '''\n",
    "    \n",
    "#   Number of points to perform the integration over\n",
    "    points = 10000\n",
    "    x = np.linspace(0,1,points)\n",
    "#   Create integrand\n",
    "    wgt = np.sin(np.log(x+1))\n",
    "    p1 = x**n\n",
    "    p2 = x**m\n",
    "    integrand = wgt*p1*p2/points\n",
    "    \n",
    "#   Integral is just sum\n",
    "    return sum(integrand)\n",
    "\n",
    "def intpoly(L,n,m):\n",
    "    '''\n",
    "    Produce the integral of two polynomials from 0 to 1 with a \n",
    "    specific weight function:\n",
    "    \n",
    "    \\int_0^1 Sin(Ln(x+1)) P_n(x) P_m(x)    \n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    L: array, float\n",
    "        Coefficents for set of polynomials\n",
    "    n: int\n",
    "        order of first polynomial\n",
    "    m: int\n",
    "        order of second polynomial\n",
    "    Returns\n",
    "    ----------\n",
    "    float\n",
    "        result of integration\n",
    "    \n",
    "    '''\n",
    "    \n",
    "#   Number of points to perform the integration over\n",
    "    points = 10000\n",
    "    x = np.linspace(0,1,points)\n",
    "    \n",
    "#   Create integrand using polynomial generator from numpy\n",
    "    wgt = np.sin(np.log(x+1))\n",
    "    p1 = np.poly1d(L[n,::-1])\n",
    "    p2 = np.poly1d(L[m,::-1])\n",
    "    integrand = wgt*p1(x)*p2(x)/points\n",
    "    \n",
    "#   Integral is just sum\n",
    "    return sum(integrand)\n",
    "\n",
    "# Vectorise just converts a scalar function to one that accepts vectors\n",
    "# Needed for numpy.fromfunction\n",
    "tmp = np.vectorize(intmono)\n",
    "\n",
    "# Creates an array with entries calculated from the function which takes array indexes as input\n",
    "M = np.fromfunction(tmp,(6,6))\n",
    "\n",
    "# Cholesky decomposition. This is equivalent to Gram-Schmidtt\n",
    "# orthonormalisation but more stable\n",
    "L = np.linalg.cholesky(M)\n",
    "\n",
    "# Calculate inverse\n",
    "L = np.linalg.inv(L)\n",
    "\n",
    "# Check polynomials are now orthonormal\n",
    "I = np.zeros((6,6), dtype='float') \n",
    "for i in range(6):\n",
    "    for j in range(6):\n",
    "        I[i,j] = intpoly(L,i,j)\n",
    "\n",
    "# This set global precision level for print statements    \n",
    "np.set_printoptions(precision=3)\n",
    "print(M)\n",
    "print(I)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Exercise 2\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "% matplotlib inline\n",
    "\n",
    "ang = np.linspace(0, np.pi)\n",
    "\n",
    "plt.plot(np.cos(ang), np.sin(ang), linewidth = 10, color = 'b')\n",
    "plt.plot(0.96*np.cos(ang), 0.96*np.sin(ang), linewidth = 10, color = 'g')\n",
    "plt.plot(0.92*np.cos(ang), 0.92*np.sin(ang), linewidth = 10, color = 'r')\n",
    "plt.plot(0.88*np.cos(ang), 0.88*np.sin(ang), linewidth = 10, color = 'c')\n",
    "plt.plot(0.84*np.cos(ang), 0.84*np.sin(ang), linewidth = 10, color = 'm')\n",
    "plt.plot(0.80*np.cos(ang), 0.80*np.sin(ang), linewidth = 10, color = 'y')\n",
    "plt.plot(0.76*np.cos(ang), 0.76*np.sin(ang), linewidth = 10, color = 'k')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "p38workshop",
   "language": "python",
   "name": "p38workshop"
  },
  "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": 0
}