{
 "metadata": {
  "name": "",
  "signature": "sha256:3d9763ad8e713bc4c7be2162464d383f293b498ef88147f61f44be6689e11eb8"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Linear algebra for multidimensional polynomial fitting"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import itertools\n",
      "# A small data set\n",
      "#T_in = np.array([20,24,35,40])+273.15\n",
      "#x_in = np.array([47,55,70,78,82])/100.0\n",
      "#rho_in = np.array([[1047,1033,1020,1000,990],[997,983,970,950,940],[947,933,920,900,895],[987,883,870,850,845]])\n",
      "#\n",
      "# large data set\n",
      "T_in   = np.array([   -45 ,    -40 ,    -35 ,    -30 ,    -25 ,    -20 ,    -15 ,    -10])+273.15 # Kelvin\n",
      "x_in   = np.array([     5 ,     10 ,     15 ,     20 ,     25 ,     30 ,     35 ])/100.0 # mass fraction\n",
      "        \n",
      "rho_in = np.array([\n",
      "          [1064.0,    1054.6,    1045.3,    1036.3,    1027.4,    1018.6,    1010.0],\n",
      "          [1061.3,    1052.1,    1043.1,    1034.3,    1025.6,    1017.0,    1008.6],\n",
      "          [1057.6,    1048.8,    1040.1,    1031.5,    1023.1,    1014.8,    1006.7],\n",
      "          [1053.1,    1044.6,    1036.2,    1028.0,    1019.9,    1012.0,    1004.1],\n",
      "          [1047.5,    1039.4,    1031.5,    1023.7,    1016.0,    1008.4,    1000.9],\n",
      "          [1040.7,    1033.2,    1025.7,    1018.4,    1011.2,    1004.0,     997.0],\n",
      "          [1032.3,    1025.3,    1018.5,    1011.7,    1005.1,     998.5,     992.0],\n",
      "          [1021.5,    1015.3,    1009.2,    1003.1,     997.1,     991.2,     985.4]]) # kg/m3"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Minimizing the squared error $\\epsilon(\\mathbf{c}) = \\sqrt{\\sum (\\mathbf{z} - \\mathbf{A} \\cdot \\mathbf{c})^2 }$ can be achieved by solving the system of orthogonal equations given by $\\mathbf{A}^\\text{T}\\mathbf{A} \\cdot \\mathbf{c} =\\mathbf{A}^\\text{T}\\mathbf{z}$. Using Python tools, we can leave this to the software and we only have to construct the Van der Monde matrix $\\mathbf{A}$ of the independent variable and equate it with the result vector of the dependent variable."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def getCoeffs1d(x,z,order):\n",
      "    if (len(x)<order+1): \n",
      "        raise ValueError(\"You have only {0} elements and try to fit {1} coefficients, please reduce the order.\".format(len(x),order+1))\n",
      "    A = np.vander(x,order+1)[:,::-1]\n",
      "    #Anew = np.dot(A.T,A)\n",
      "    #znew = np.dot(A.T,z)\n",
      "    #coeffs = np.linalg.solve(Anew, znew)\n",
      "    coeffs = np.linalg.lstsq(A, z)[0]\n",
      "    return coeffs\n",
      "\n",
      "xorder = 2\n",
      "x = T_in[0:4]#T_in#T_in[0:4]#\n",
      "z = rho_in[0:4,0]#rho_in#rho_in[0:4,0]#\n",
      "coeffs = getCoeffs1d(x,z,xorder)\n",
      "print coeffs\n",
      "zf = np.polynomial.polynomial.polyval(x,coeffs)\n",
      "print z[0]\n",
      "print zf[0]\n",
      "print z[-1]\n",
      "print zf[-1]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We can extend the whole procedure to 2D, given that we have a solution matrix $\\mathbf{Z}$ instead of a vector. Since this potentially involves a large number of coefficients, we disregard the higher order terms in order to avoid overfitting. This is done by discarding terms with a sum of exponents higher than the largest single exponent. The pair of exponents for elements $\\mathbf{x}$ and $\\mathbf{y}$, $i$ and $j$, has to satisfy $ i+j \\leq \\max(k,l)$ with $k$ and $l$ being the highest exponents in $x$ and $y$ direction respectively. The matrix of exponent pairs $\\mathbf{E}$ for $k<l$ is therefore defined as $$ \\mathbf{E}_{i,j} =\n",
      " \\begin{pmatrix}\n",
      "   (0,0) &  (0,1) & \\cdots  & \\cdots  & \\cdots  & (0,l)  \\\\\n",
      "   (1,0) & \\ddots &         &         & (1,l-1) & (0,0)  \\\\\n",
      "  \\vdots &        &         & \\udots  & \\udots  & \\vdots \\\\\n",
      "   (k,0) & \\cdots & (k,l-k) &  (0,0)  & \\cdots  & (0,0)  \n",
      " \\end{pmatrix}\n",
      " \\text{ yielding, for example, }\n",
      "  \\begin{pmatrix}\n",
      "   (0,0) & (0,1) & (0,2) & (0,3) & (0,4) \\\\\n",
      "   (1,0) & (1,1) & (1,2) & (1,3) & (0,0) \\\\\n",
      "   (2,0) & (2,1) & (2,2) & (0,0) & (0,0) \n",
      " \\end{pmatrix}\n",
      " \\text{ for $k=2$ and $l=4$. }\n",
      " $$\n",
      "Following the matrix notation would result in a 4-dimensional functional matrix from the Cartesian product of all elements of $\\mathbf{E}$ and the input vectors $\\mathbf{x}$ and $\\mathbf{y}$. However, having many linear algebra solvers available in 2D, it is more practical to manually reduce the dimensionality of $\\mathbf{A}$ to two. Each row of $\\mathbf{A}$ corresponds to an unique pair elements of the input vectors $\\mathbf{x}$ and $\\mathbf{y}$. Evaluating the resulting expression $x_n^i y_n^j$ for all non-zero entries of $\\mathbf{E}$ fills the columns of $\\mathbf{A}$ with values. Every row of $\\mathbf{A}$ can therefore be mapped to an entry in the solution matrix $\\mathbf{Z}$. A data set consisting of 10 entries in $x$ and 20 in $y$ requires $\\mathbf{Z}$ to have 200 elements and thus $\\mathbf{A}$ to have 200 rows. The example given above, $k=2$ and $l=4$, leads to 10 columns in $\\mathbf{A}$. After minimizing $\\epsilon(\\mathbf{c})$, information from $\\mathbf{E}$ can be used to convert the coefficient vector $\\mathbf{c}$ to a matrix to be used with two-dimensional polynomials."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def getCoeffs2dmatrix(x_in,y_in,z_in,x_order,y_order):\n",
      "    \n",
      "    x_order += 1\n",
      "    y_order += 1\n",
      "    \n",
      "    #To avoid overfitting, we only use the upper left triangle of the coefficient matrix\n",
      "    x_exp = range(x_order)\n",
      "    y_exp = range(y_order)\n",
      "    limit = max(x_order,y_order)\n",
      "    \n",
      "    xy_exp = []\n",
      "    \n",
      "    # Construct the upper left triangle of coefficients    \n",
      "    for i in x_exp:\n",
      "        for j in y_exp:\n",
      "            if(i+j<limit): xy_exp.append((i,j))\n",
      "                \n",
      "    x_num = len(x_in)\n",
      "    y_num = len(y_in)\n",
      "                \n",
      "    cols = len(xy_exp)\n",
      "    eqns = x_num * y_num\n",
      "    #if (eqns<cols):\n",
      "    #    raise ValueError(\"You have only {0} equations and try to fit {1} coefficients, please reduce the order.\".format(eqns,cols))   \n",
      "    if (x_num<x_order):\n",
      "        raise ValueError(\"You have only {0} x-entries and try to fit {1} x-coefficients, please reduce the x_order.\".format(x_num,x_order))\n",
      "    if (y_num<y_order):\n",
      "        raise ValueError(\"You have only {0} y-entries and try to fit {1} y-coefficients, please reduce the y_order.\".format(y_num,y_order))\n",
      "    \n",
      "    #Create functional matrix\n",
      "    A = np.zeros((x_num,y_num,x_order,y_order))\n",
      "    for i in range(x_num):\n",
      "        for j in range(y_num):\n",
      "            for (xk,yk) in xy_exp:\n",
      "                A[i][j][xk][yk] = x[i]**xk * y[j]**yk\n",
      "                \n",
      "    raise NotImplementedError(\"No 4-dimensional solver implemented\")\n",
      "    \n",
      "    \n",
      "def getCoeffs2d(x_in,y_in,z_in,x_order,y_order):\n",
      "    \n",
      "    x_order += 1\n",
      "    y_order += 1\n",
      "    \n",
      "    #To avoid overfitting, we only use the upper left triangle of the coefficient matrix\n",
      "    x_exp = range(x_order)\n",
      "    y_exp = range(y_order)\n",
      "    limit = max(x_order,y_order)\n",
      "    \n",
      "    xy_exp = []\n",
      "    \n",
      "    # Construct the upper left triangle of coefficients    \n",
      "    for i in x_exp:\n",
      "        for j in y_exp:\n",
      "            if(i+j<limit): xy_exp.append((i,j))\n",
      "    \n",
      "    # Construct input pairs   \n",
      "    xx, yy = np.meshgrid(x_in,y_in,indexing='ij')\n",
      "    xx = np.array(xx.flat)\n",
      "    yy = np.array(yy.flat)\n",
      "    zz = np.array(z_in.flat)\n",
      "    \n",
      "    x_num = len(x_in)\n",
      "    y_num = len(y_in)\n",
      "                \n",
      "    cols = len(xy_exp)\n",
      "    eqns = x_num * y_num\n",
      "    #if (eqns<cols):\n",
      "    #    raise ValueError(\"You have only {0} equations and try to fit {1} coefficients, please reduce the order.\".format(eqns,cols))   \n",
      "    if (x_num<x_order):\n",
      "        raise ValueError(\"You have only {0} x-entries and try to fit {1} x-coefficients, please reduce the x_order.\".format(x_num,x_order))\n",
      "    if (y_num<y_order):\n",
      "        raise ValueError(\"You have only {0} y-entries and try to fit {1} y-coefficients, please reduce the y_order.\".format(y_num,y_order))\n",
      "    \n",
      "    # Build the functional matrix\n",
      "    A = np.zeros((eqns,cols))\n",
      "    for i in range(eqns): # row loop\n",
      "        for j, (xj,yj) in enumerate(xy_exp): # makes columns\n",
      "            A[i][j] = xx[i]**xj * yy[i]**yj\n",
      "             \n",
      "    coeffs = np.linalg.lstsq(A, zz)[0]\n",
      "    \n",
      "    #Rearrange coefficients to a matrix shape\n",
      "    C = np.zeros((x_order,y_order))\n",
      "    for i, (xi,yi) in enumerate(xy_exp): # makes columns\n",
      "        C[xi][yi] = coeffs[i]\n",
      "        \n",
      "    return C\n",
      "\n",
      "xorder = 3\n",
      "yorder = 3\n",
      "\n",
      "x = T_in\n",
      "y = x_in\n",
      "z = rho_in\n",
      "\n",
      "coeffs = getCoeffs2d(x,y,z,xorder,yorder)\n",
      "print coeffs\n",
      "\n",
      "print z[0][0]\n",
      "print np.polynomial.polynomial.polyval2d(x[0],y[0],coeffs)\n",
      "\n",
      "print z[-1][-1]\n",
      "print np.polynomial.polynomial.polyval2d(x[-1],y[-1],coeffs)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%pylab inline\n",
      "import matplotlib\n",
      "matplotlib.rcParams['savefig.dpi'] = 2 * matplotlib.rcParams['savefig.dpi']\n",
      "from mpl_toolkits.mplot3d import Axes3D\n",
      "\n",
      "# Construct input pairs   \n",
      "xx, yy = np.meshgrid(x,y,indexing='ij')\n",
      "xx = np.array(xx.flat)\n",
      "yy = np.array(yy.flat)\n",
      "zz = np.array(z.flat)\n",
      "\n",
      "zf = np.polynomial.polynomial.polyval2d(xx,yy,coeffs)\n",
      "\n",
      "fig  = plt.figure()\n",
      "ax   = plt.axes(projection='3d')\n",
      "ax.scatter(xx, yy, zz, color='blue', label='Original data')\n",
      "ax.scatter(xx, yy, zf, color='red',  label='Fitted data')\n",
      "#ax.plot(xx[0], yy[0], zs=zz[0], label='Fitted line', zdir='z')\n",
      "ax.legend()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}