{ "metadata": { "name": "" }, "nbformat": 3, "nbformat_minor": 0, "worksheets": [ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Scaling the Coefficients\n", "\n", "This notebook aims to scale the influence coefficients and run a large Cahn-Hilliard simulation using an MKS model. The notebook will\n", "\n", " - create some training data\n", " - demonstrate the `resize_coeff` method\n", " - examine how the mean square error changes when the coefficients are resized\n", " - run a large Cahn-Hilliard simulation" ] }, { "cell_type": "code", "collapsed": false, "input": [ "%matplotlib inline\n", "%load_ext autoreload\n", "%autoreload 2\n", "\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from pymks import MKSRegressionModel\n", "from pymks import FiPyCHModel\n", "from sklearn import metrics\n", "\n", "mse = metrics.mean_squared_error" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 1 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create the training data on a 21$\\times$21 grid" ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.random.seed(101)\n", "X_train = np.random.random((200, 21, 21))\n", "fipymodel = FiPyCHModel()\n", "y_train = fipymodel.predict(X_train)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 2 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## How to scale the coefficents?\n", "\n", "The `resized_coeff` method scales the coefficients. For example" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The `resize_coeff` method pads the coefficients with zeros for an arbitrary dimensional array. This assumes that the coefficients become small away from the origin." ] }, { "cell_type": "code", "collapsed": false, "input": [ "model = MKSRegressionModel()\n", "coeff = np.arange(16).reshape((4, 4, 1))\n", "model.Fcoeff = np.fft.fftn(coeff, axes=(0, 1))\n", "print coeff[:,:,0]" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0 1 2 3]\n", " [ 4 5 6 7]\n", " [ 8 9 10 11]\n", " [12 13 14 15]]\n" ] } ], "prompt_number": 3 }, { "cell_type": "markdown", "metadata": {}, "source": [ "becomes" ] }, { "cell_type": "code", "collapsed": false, "input": [ "model.resize_coeff((6, 6))\n", "print np.fft.ifftn(model.Fcoeff, axes=(0, 1))[:, :, 0].astype(int)" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "[[ 0 1 0 0 2 3]\n", " [ 4 5 0 0 6 7]\n", " [ 0 0 0 0 0 0]\n", " [ 0 0 0 0 0 0]\n", " [ 8 9 0 0 10 11]\n", " [12 13 0 0 14 15]]\n" ] }, { "output_type": "stream", "stream": "stderr", "text": [ "-c:2: ComplexWarning: Casting complex values to real discards the imaginary part\n" ] } ], "prompt_number": 4 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 05-0\n", "\n", "Without viewing the `resize_coeff` method try and take the \n", "\n", "```python\n", "[[ 0 1 2 3],\n", " [ 4 5 6 ],\n", " [ 8 9 10 11],\n", " [12 13 14 15]]\n", "```\n", "\n", "array and reproduce the centrally padded `(6, 6)` array without using loops. Use `np.array_split` and `np.concatenate`.\n", " " ] }, { "cell_type": "code", "collapsed": false, "input": [ "??model.resize_coeff" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 5 }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Evaluate the mean square error while varying `Nspace`" ] }, { "cell_type": "code", "collapsed": false, "input": [ "Ns = [21, 42, 84, 168]\n", "errors = []\n", "np.random.seed(102)\n", "\n", "for Nspace in Ns:\n", " print Nspace,\n", " X_test = np.random.random((2, Nspace, Nspace))\n", " fipymodel = FiPyCHModel()\n", " y_test = fipymodel.predict(X_test)\n", " model = MKSRegressionModel(Nbin=5)\n", " model.fit(X_train, y_train)\n", " model.resize_coeff(X_test.shape[1:])\n", " y_pred = model.predict(X_test)\n", " errors.append(mse(y_test, y_pred))" ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "21 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "42 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "84 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "168\n" ] } ], "prompt_number": 6 }, { "cell_type": "code", "collapsed": false, "input": [ "plt.semilogx(np.array(Ns)**2, errors)\n", "plt.title(r'Mean Square Error (MSE) v Domain Size ($S$)')\n", "plt.ylabel(r'MSE')\n", "a = plt.xlabel(r'$S$')" ], "language": "python", "metadata": {}, "outputs": [ { "metadata": {}, "output_type": "display_data", "png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAETCAYAAADQ97psAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X1wI/V5B/DvmiOhIZwl+UJDCMdZdnkpaWrLuktawmCQ\nfZSESQHbZyaQpGnRyUnbafuHX65NB1NajpOm7UyGDlhipnfkBc6SSRmG0jtrwZSElCLvHiEQXs6r\nwAEhIZYlCD0I3G3/2NlFkvVqr7wr7fcz47GkXe0+uz9pH+3vZVdQVVUFERE5UpvVARARkXWYBIiI\nHIxJgIjIwZgEiIgcjEmAiMjBmARMls1m0dbWhrGxsVXTRkZG0Na2MbtckiQMDg6ira0NHo8Hu3bt\nQi6X25B116OtrW3V3ymnnII333yz4esOh8OQZbmuMqu2X8ttTy6Xw9zcHGKx2LpiLl623++HLMvr\nWmY9JEmC3++v+z3l9tlalleJXqbVmFEWLUMlU62srKiCIKgej2fVNJfLpba1tW1IHC6XS52amlLT\n6bSazWbVUCik9vX1bci66yEIgirLsprL5Qr+Gm1lZcXYH/WUWbX9Wm17urq61hV3/vLT6bQaDodV\nQRBURVHWtdxaZbNZNZlM1vWeSvtsLcsrJ79M8yUSCTWZTBp/iURCVdX1l0WraMkksLS0pCYSCTUa\njarZbHZD160fUHbu3KlKkmS8Pj8/r46MjKiCIBS85vV6VUEQ1MHBwYJYZ2ZmVLfbrQqCoPb19Rlf\n8qWlJdXn86nhcFh1u91qV1dXwXr0efLXo9u5c2fB8l0ul+rxeIx1qaqqLi4uFnyR5ufn1cHBwZri\nGhgYUCcmJoz3V9o+nSAIajqdLrkvi5epKMqqdcTjcWMdIyMjxjpKxZNv3759aiwWU1W19jKrZb9W\n2h5VVdXJyUnjIJRvYGBAjUajxvOJiQk1FAqtmq/U8kOhUMG8lfaJz+dTQ6GQ6na71cHBQXVxcdGY\nNxwOG8soV875nw8zPov5y4vH46ogCAV/oijW9DlS1cIyzX8tPwlHo1FVlmVVVcuXhdM0ZRKYmJgo\neK5nev1LNDk5qeZyOVVRlIIP9kbQDyixWKzgixkKhdREImF8IVZWVlS3260+/PDDai6XU0OhkDoy\nMlKwjCNHjqjZbFYdGRlRJycnVVX94EsViUSM5eYfpHVdXV3q4OBgyV9Zi4uLqtvtVmVZVrPZrOrz\n+YxfwZWSQC1xjY2Nqel0uuL25av0K7Z4meWei6JoxKPv8+J5i/X19Rmv11pm1fZrte1RVVVNJpMl\n90M0Gi14va+vTxVFseTyi7cnmUwWHJir7ZO5uTk1m82qXq9Xdbvdai6XU5PJZMFns1w5FycBMz6L\npZJ0PB5Xu7u71Ww2W9PnSN9nxfumOJ78JFWuLJym6ZLAzMxMwWnc4uKikc2j0agqSZIaDodVRVHU\n+fl548O7UfQvkP7h1emP9S/azMzMqg9g/vz5v3Z2795d8EXOn0+SpLLVPNFoVB0cHDR+9ekHp4mJ\nCXVqasqYL5lM1nwmUCmu/INlte3TFf/yEwRB7e7uLrnM4uf79u1Tx8bGjOeKohjrKPcLNH+9ulrL\nTFduv1bbHj2uUvtBj0F/XGoeffnFBzo9qdeyT/KXGwqFCj4HgiAYv5rLlXNxEljvZ7FUEtC3P51O\n1/w50uMvtd6uri51cnJyVXIuVxZOs8nqNol67d69G4lEwng+OzuLnTt3AgC8Xi+SySTGx8chiiJy\nuRx27NhhSZzt7e1Go93y8jK2b99eMH1paQmJRAIej8d4TRAE4/Gtt94KURSN6V6v15iW/x61wlU/\ngsEggsEgACAWi6GrqwtLS0tIp9PGPgOAzs7OmrerUlz5j6ttX75kMll2+/JfL36eyWQKnnd2diKb\nzZZ9ry6bzcLlcq16vVqZ6UrtV0VRsG3btpq2Jz9Gncvlgs/ngyiKWFpawujoaMl1l5K/H6rtk/xY\nXC5XwfN8lco533o/i6UEAgGEw2Fs27at5s9RuTINBoMYHR3F/Pw8BgcHEY/H0dvba2xTqbJwGtv1\nDpJlGel02niuH8zLyWazBR+Q5eVlo8dHJpPBtdde29B4KxkZGcG9996LRCKBUChUMK27uxvDw8PI\nZDLGXzKZBAAkEgmIooiHH34Yhw8fxvDwcF3rTSQSBQd5QPsy+Hw+SJIEj8eDo0ePGtMURSm7rPwv\nST1xVdq+Yl6vF9u2bTP+Nm/eXNN2dnR0FBxIyh0I6lGpzKrtV91at2d0dBSHDx8uue5K4vG4kbDM\n2Cfr/fwVL6vcPpNledUBfXJyElu2bMGNN94IoL7PUTG9l9DmzZsxNDSEmZkZpFKpNW9Ly7L6VKSU\naDSqKopitOYXy6+eCIVCRj1fMpnc8OqfYsWn9S6XS+3q6jJOs4tP+ZPJpLqysqLu3r3baCzTT531\n+fr6+oxT4qWlpVXVYcWn09ls1mjoW1paUpeWltR9+/apbW1tRv2v2+1WJUlSV1ZWCtoE9GoURVGM\naWuJS69ayd++UvXFgiAYceT/lVpm8XNFUVRBEIx1DA8PG1UhxfOWWu9ayqzSfq22PXpc5aogFEUx\nGusrxa0vP3/9ehVRPftkcnKyoM1Mrw6qVM7F1UHr/Szmv2d+fl51u90FVVGlvielPkd6/PmKjwXh\ncLigkZjVQRpbJgFV1Q7u+b0l8uV/CCYnJ41EEY/HN7whuNjKykpBl8K+vr6C3iP505LJpNrV1WX0\nTMn/gOr1p36/3zhoi6KoLi0tFdQxLy4uqn6/f1UciqIYyxAEQfX7/QUNjXqPDo/HYzzWhUIhoy47\nkUgUxF9rXNW2T1eqDl1v2CxeZql1JBIJYx27du0y1lFq3nz5PV7qKbNq+7XS9qiqdqDbtWtX2bi6\nuroK6vSr7S+/32/0dql3n0xOThqNuvp26vOWK2dJkozPmxmfxfz3DA8Pq21tbQXbF4vFavocqWph\nmUqSpEajUTUajRo9BYvbUqqVhVM0LAlUOxhXmq6fAehnBMXyk4Be2Poyi78QVF2lhshWFQ6Hy/7I\naKSJiQl1bm5uw9frBPWWKctC05AkUNyjpJ7pkiQVVAEV9/WPx+Oq2+0u6A8cjUYLuohSfZyYBLLZ\nrCWD5zhAqXHqLVOWhaYhvYPK9QKpZbrecq/TexTohoeHVzVU6fMEAoF6wqQ81cqs1bS3t2N0dBSy\nLK/6zDXK3NwcpqamNmRdTlRPmbIsPiCoqrl3FtMLYOfOnTh8+HDd04mIaOOY3kU0k8msazoREW0c\nU5OALMsVq2SqTScioo1lapuAoihQFAXLy8vIZDJG1Y8+YKXc9HxOq5smIjLLWmr3TT0TGBoawtDQ\nEARBQC6XMw7oAwMDFacXU7VeS5b+3XTTTbZYXq3vq2W+avOUm17P62bvNzuUX6PLzozyq3davWXd\nrGVnl/LbiO/eWjXkshHBYBAvvvgienp6AGDVUO3i6XbU399vi+XV+r5a5qs2T7np9b5uB2bG1uiy\nq3XeSvPUO80pZbee5ZlZfnb+7pneO2i9BEFYV1Yja01PT2N6etrqMGiNWH7Na63HTttdQI6am51/\nYVJ1LD/n4ZkAEVEL4JkAERHVjUmAiMjBmASIiByMSYCIyMGYBIiIHIxJgIjIwZgEiIgcjEmAiMjB\nmASIiByMSYCIyMGYBIiIHIxJgIjIwZgEiIgcjEmAiMjBmASIiByMSYCIyMGYBIiIHIxJgIjIwZgE\niIgcjEmAiMjBmASIiByMSYCIyMGYBIiIHIxJgIjIwZgEiIgcrGFJIBKJlHw9FoshFothamqqUasm\nIqIaNSQJJJNJzM/Pr3pdFEUMDAwgGAxCURSIotiI1RMRUY02NWKhgiCUfF1RFCiKgmAwCK/XC0VR\nEAgEGhECtbh33gGuvhr4rd8Ctm4Fzj1X+9Mff+xjQJmPIRHlMT0JyLKMQCCAffv2rZoWDAaNx5Ik\n4brrrjN79eQQDzwAvPUWEAwCL70EvPwy8NhjHzx+++0PEkLx/3PPBc4+G/jQh6zeCiLrmZ4EMplM\n1XkkSUJfXx96enrMXj05xP79wNe/DgwNlZ7+618Dx45pSUFPDPPz2v+XXgJ+/nPgzDNLn0Xo/zdv\n3tBNIrKEqUlAPwuoRhRF7N2718xVk4O89hrw+OPA7Gz5eT76UeDCC7W/Ut5/X1uOniBeegn48Y+1\nMww9cZx6avkEsXUr8PGPA23sX0dNztQkoNf5Ly8vI5PJQJZl9Pb2IpvNwuVyAQCi0SjGx8cBaMmg\nVNKYnp42Hvf396O/v9/MMKnJfec72hnA6aevfRmbNmkH8q1bS09XVSCT+SBB6P+feOKDx7kc8MlP\nlq5u2roVOOcc4LTT1h4jUSULCwtYWFhY93IEVVXV9YdTKBaLIRwOIx6Po6enB36/H6lUCslkErt2\n7YLH40Emk0EikcDll19eGJAgoAEhUYtQVeCii4BYDLj4YmtjOX5cq3LSk0J+snj5ZeCVVwC3u/LZ\nhNvNBmwyx1qPnQ1JAuvBJECV/O//AtdfD7zwgv0PnidOAK+/vvpsIv/xyZPlG6+3bgU+8QnglFOs\n3hJqBkwC5Ajf+IbWs+fv/s7qSMyRy5VPEC+/DPzqV8BZZ5VOEPr/j3zE6q0gO2ASoJb3zjtaApDl\n8nX5rebdd7VqpeKzCf3/sWNaI3hnJ3DjjcCf/qnW3kHOwyRALe/gQeCuu7SunqQ5eRJ44w3gmWeA\nW24BfvELYN8+4Kqr7F9dRuZiEqCWd+WVwA03aG0CtJqqAg89BExMAFu2AOEwsGOH1VHRRmESoJb2\n6qvA7/2eVjXCOvDK3n8fOHAAuOkm4HOfA/7pn4CuLqujokZb67GTQ12oKXznO8DwMBNALTZtAv7s\nz4Dnn9cS52c+A/z1X2uNzETFmATI9lRVu0zEn/yJ1ZE0l9NP13pRPfusdnZwwQXAbbdp4xuIdEwC\nZHtPPKH1uf+DP7A6kuZ05pnA7bdrl9pIpYDzz9eqi06csDoysgMmAbI9/SyAvV3W57zzgERC62UV\njQI+H3DokNVRkdXYMEy2dvy4dn2ep57S/pM5VBX4j/8Apqa0QWfhMMCL+jY3NgxTS7r/fsDvZwIw\nmyAA11wD/OQn2v8rrwS+8hVtABo5C5MA2RobhBvr1FO1+zK88AKwbZtWRTQxAWSzVkdGG4VJgGzr\n1Ve1C8ZdfbXVkbS+M84A/uEfgKef1hLAeecB//qv2mUrqLUxCZBt3X03MDKi3UeYNsYnPqE1Gi8s\nAI88onUrvece7fIU1JrYMEy2pKraAWj/fnYNtdKjjwLj41oSiESAyy6zOiIqhw3D1FL+53+0xsvP\nftbqSJzt0ku1cRrj49oo5Kuu0i5WR62DSYBsiWMD7EMQgNFR4Kc/BQYHgcsv1y5b/dprVkdGZmAS\nINs5fhyIx4Evf9nqSCjfhz8M/NVfadck2rJFuy7R3/898OabVkdG68EkQLbz/e9rl0A++2yrI6FS\nXC7tGkSyrN3U5rzzgH/7N+C996yOjNaCSYBsh2MDmsPWrVpZHToEPPAAcNFFwNyc1qhPzYO9g8hW\njh3TLl/w6qvAaadZHQ3VY35ea0D+yEe0nkQXX2x1RM7C3kHUEr79bWDXLiaAZjQ4CEiSNgL5S18C\nrr1Waz8ge2MSINvgfQOaX1ub1qD//PNa997PfQ748z/X7n1M9sQkQLbx+OPAKafwvrit4LTTtGsQ\nPfec1qvooouAW24B3n7b6sioGJMA2QbHBrSejg7gX/5FuwbUT3+q9SSKxbQ7nZE9sGGYbOH//k+7\nXPRPfqJdv4Za05NPao3Hb7wB7NsHfOELTPpmYcMwNbXvf1+rQ2YCaG3bt2sXpguHgclJ7VpETz5p\ndVTOxiRAtsAGYecQBO0M4KmntEbkq68GrrsOUBSrI3OmhiWBSCRS8vW5uTmIoohYLNaoVVOTefll\nrWvhF79odSS0kTZt0i5K98ILwKc+pXUI+Ju/AZaXrY7MWRqSBJLJJObn51e9LkkSACAQCAAAZFlu\nxOqpydx9t3aBMo4NcKbTTwe++U3g2We1S09ccIHWXnD8uNWROUNDkoBQpqVndnYWbrcbAOD1epFM\nJhuxemoiHBtAujPPBG6/HfjhD7XeROefr/1AOHHC6sham+lJQJZl45d+sWw2C4/HYzxf5nmf4/3w\nh1o/8u3brY6E7OK887RrEN17LzAzA/T1AYcPWx1V6zI9CWQymYrT2f2T8nFsAJXzh38I/OAHwE03\nAX/xF8DOncCRI1ZH1XpMTQKVzgIAwOVyGUliZWUFHR0dZq6emszbb2u/+G64wepIyK4EAbjmGu1u\nZldfDVx5JfDVr2qdCcgcm8xcmKIoUBQFy8vLyGQykGUZvb29yGazcLlcGB0dRSqVQiAQQDqdxuDg\nYMnlTE9PG4/7+/vR399vZphkE/fdp/3aO+ssqyMhuzv1VOAb39C6lEYiQG+vdnezPXu0+xs40cLC\nAhYWFta9nIaMGI7FYgiHw4jH4+jp6YHf70cqlTKmeb1eKIqCYDC4OiCOGHaMQAAYGwNGRqyOhJrN\na69p1UT33w/87d9qVy798Ietjspaaz128rIRZImXXtIa/F55hV1Dae2eeQaYmtL+33qrdhnyNocO\ngWUSoKZyyy3a5YVvv93qSKgVLCxo1yQCtOoiJ9YgMwlQ01BVoLsbOHgQ8PutjoZaxcmTwOysVj10\n223aWYGTMAlQ03jsMa0O9+mn2TWUzPfuu9p/p7URrPXYaWrvIKJa/Pu/c2wANY7TDv7rxTMB2lC/\n/jVwzjnaDUY+/nGroyFqHbyfADWF++7T7jvLBEBkD0wCtKF4sTgie2F1EG2Yn/1Mu1DcK6+w3pbI\nbKwOIts7cEC7gxQTAJF98EyANsTJk9rYgHhcGylMRObimQDZ2mOPaXeQ8vmsjoSI8jEJ0IbYvx/4\n2tc4NoDIblgdRA2njw147jngt3/b6miIWhOrg8i25uaASy5hAiCyIyYBajj9MhFEZD+sDqKGUhTg\nM58BXn0V+NCHrI6GqHWxOohs6e67gS99iQmAyK54JkANc/Ik0NWlXS+ot9fqaIhaG88EyHb++7+B\nzZuBnh6rIyGicpgEqGF43wAi+2N1EDXEW29pYwNeeAE480yroyFqfawOIltJJIBLL2UCILI7JgFq\nCP0yEURkb2tKArlczuw4qIUsLWm3j/z8562OhIiqqZgErrjiCuPx17/+deNxIBBoXETU9A4c4NgA\nomaxqdLE5eVl4/GTTz7Z8GCo+Z08qSWB+++3OhIiqgXbBMhUCwuA282xAUTNgkmATMUbyRM1l4rj\nBNra2uD1egEAiqIUPD558mRjAuI4gab15pvA1q3Aiy8CH/uY1dEQOctaj50V2wQymUzdC0wkEnC7\n3YjH47jzzjtXTZ+bm4PL5YKiKAgGg3Uvn+wrkQAuu4wJgKiZVKwOcrlcZf9KEUURoigiEAhAURQc\nOXKkYLosy/B6vQgEAvB6vZBl2bwtIcvxvgFEzadiEpBlGX6/H2+++SZkWYbH40F3dzfuu+++kvMH\nAgHccccdALSziJ4SrYOTk5MAtCqlXl5asmUcPQo8/zzHBhA1m4pJIBgMIh6PY/PmzZicnIQoijh6\n9ChuvfXWsu/J5XKIRCLYs2fPqmm9vb3o7OyEx+OBx+NZf/RkGwcOANdfD5x6qtWREFE9qvYO6uzs\nBFD7L/f29naMj49jZmYG6XS6YFo2m0V3dzdisRiCweCq6dSc9LEBvEwEUfOp2DCsE0URAwMDVeeT\nJAmCIKC3txc+nw+JRALj4+PG9FgshlAohM2bN8Plcq2arpuenjYe9/f3o7+/v5YwySKPPAJs2QJ8\n+tNWR0LkHAsLC1hYWFj3ciomgV27dqG7uxuZTAaiKCKdTiMUCmF0dLTk/KIowufzAdB+9e/YscN4\nrDcmb968GQCMxuNS8pMA2R8bhIk2XvEP5JtvvnlNy6l6PwG9R097ezvS6TQkScLQ0FDJeXO5HGZn\nZwFo1Ud79+4FAPj9fqRSKQBAJBKB1+tFJpMp2UWU4wSaSy4HnHuu1jC8ZYvV0RA511qPnRWTwNjY\nWMkFC4Jg9AIyG5NAc7nrLuA//1O7jzARWachg8UOHjyIjo4ODA8PY3BwcM3BUevavx+YmLA6CiJa\nq6rVQZIkYWZmBouLixgYGMDY2Bi2bdvWuIB4JtA0XnwRuOQS4Ngxdg0lslpDqoOKiaKImZkZyLKM\nF198se6V1RQQk0DT+OY3gePHgX/+Z6sjIaKGVAflE0UR8XgcS0tL2L17d90rotZy4oQ2NuDBB62O\nhIjWo2ISkGW5oCooFAqVvCgcOc/DD2s3kefYAKLmVtOlpPW+/8abBAEHDx5sTECsDmoK118PfPaz\nwF/+pdWREBHQoOqgw4cPGwsHYKxAf07OlMtp1UDf+pbVkRDRetXVMLwReCZgf7EYcOiQdv8AIrKH\ntR47eXtJqhsvE0HUOngmQHV5/nng0ks5NoDIbngmQBviwAHghhuYAIhaBc8EqGYnTmgXi/uv/wI+\n9SmroyGifDwToIYTReCss5gAiFoJkwDVjA3CRK2H1UFUk2xWqwpSFKCjw+poiKgYq4OooQ4eBHbu\nZAIgajVMAlST/ftZFUTUipgEqKrnngN+9jPgiiusjoSIzMYkQFUdOAB8+cvAppovPE5EzYINw1TR\niRPA1q3A4cPARRdZHQ0RlcOGYWqI+Xng7LOZAIhaFZMAVcQGYaLWxuogKmtlBejs1MYGeDxWR0NE\nlbA6iEx38KDWI4gJgKh1MQlQWbxMBFHrYxKgkp59VrtnwOCg1ZEQUSMxCVBJHBtA5AxsGKZV3n9f\nGxsgisCFF1odDRHVwjYNw4lEAqIoYmxsrOR0SZIwNzeHWCxm9qrJJPPzWhJgAiBqfaYmAVEUIYoi\nAoEAFEXBkSNHVs1z2223YWhoCNlsFrIsm7l6MgnHBhA5R8Oqg/x+P1KpVMFriUQC6XQa4+Pj5QNi\ndZClMhnA6wXSacDttjoaIqqVbaqDcrkcIpEI9uzZs2paKpXC8vIyZFlGJBIxe9VkgnvvBf7oj5gA\niJzC9CTQ3t6O8fFxzMzMIJ1Or5q+ZcsW9Pb2AgDm5ubMXj2tE6uCiJzF1A6AkiRBEAT09vbC5/Mh\nkUgUVP10dHSgs7MTAOByufDkk09iaGho1XKmp6eNx/39/ejv7zczTCrjmWeAV1/l2ACiZrCwsICF\nhYV1L8fUJCCKInw+HwAgm81ix44dxmOXy4Xh4WEkEolV04vlJwHaOAcOAF/5CnDKKVZHQkTVFP9A\nvvnmm9e0HFMbhnO5HGZnZwEAiqJg7969AAobiWOxGDweD1KplDG9ICA2DFvi/feBc84BHnkEuOAC\nq6Mhonqt9djJwWIEAHjwQeAf/xH40Y+sjoSI1sI2vYOoObFBmMiZeCZAWF4Gurq0m8m7XFZHQ0Rr\nwTMBWrN77wU+/3kmACInYhIg3jeAyMGYBBzu6aeB118HAgGrIyEiKzAJOBzHBhA5GxuGHey997Sx\nAY8+Cpx/vtXRENF6sGGY6nbokHbFUCYAIudiEnCw/fuBr33N6iiIyEqsDnKoX/0K6O4GXnoJaG+3\nOhoiWi9WB1Fd7rkH+MIXmACInI5JwKF4mQgiApgEHOnHPwZ++Uvg8sutjoSIrMYk4EAHDgBf/SrH\nBhARG4Yd5733gE9+EvjBD4Df+R2royEis7BhmGry0EPawZ8JgIgAJgHHYYMwEeVjdZCDvPGGdgbA\nsQFErYfVQVTVPfcAV13FBEBEH2AScBBeJoKIijEJOMRTT2m3kbzsMqsjISI7YRJwiP37tfsGtLHE\niSgPG4Yd4De/0cYGPP64dtE4Imo9bBimsh56SLtnABMAERVjEnAANggTUTmsDmpxv/wlcN55wLFj\nwBlnWB0NETUKq4OopO99D/jiF5kAiKg0JoEWx8tEEFElpieBRCIBURQxNjZWcb5IJGL2qqnIkSPA\nygrQ3291JERkV6YmAVEUIYoiAoEAFEXBkSNHSs6XTCYxPz9v5qqphP37tfsGcGwAEZWzycyFBQIB\nBAIBAEAmk0FPT0/J+QRBMHO1VMJvfqO1B/zoR1ZHQkR2ZvpvxFwuh0gkgj179pScLsuykSiocR58\nELjwQqCry+pIiMjOTE8C7e3tGB8fx8zMDNLp9KrpmUzG7FVSCWwQJqJamJoEJEmCLMsAAJ/Ph0Qi\nUTCdZwEb4xe/AB59FBgetjoSIrI7U9sERFGEz+cDAGSzWezYscN47HK5oCgKFEXB8vIyMpkMZFlG\nb2/vquVMT08bj/v7+9HP7i11+d73gD/+Y44NIGplCwsLWFhYWPdyTB0xnMvlMDs7CwBQFAV79+4F\nAPj9fqRSKWO+WCyGcDiMeDy+qvGYI4bXR1WB3/994FvfYtdQIidZ67GTl41oMZIEDA0BS0vsGkrk\nJLxsBAHg2AAiqg/PBFrIu+9q9w144gnA67U6GiLaSDwTIDz4IHDRRUwARFQ7JoEWwrEBRFQvVge1\niNdf10YIHzsGfPSjVkdDRBuN1UEO993vAldfzQRARPVhEmgBqsqqICJaGyaBFiBJwNtvA5dcYnUk\nRNRsmARaAMcGENFasWG4yb37LnD22UAqBWzbZnU0RGQVNgw71AMPAJ/+NBMAEa0Nk0CTY4MwEa0H\nq4Oa2M9/Dvzu7wKvvAKcfrrV0RCRlVgd5EDf/S5wzTVMAES0dkwCTYpjA4jIDEwCTSqVAo4f59gA\nIlofJoEmpZ8FCILVkRBRM2PDcBN65x1tbIAkAeeea3U0RGQHbBh2kAceAHp6mACIaP2YBJoQG4SJ\nyCysDmoyHBtARKWwOsghvv1tYGiICYCIzMEk0EQ4NoCIzMYk0ESefRZ47z3g4outjoSIWgXbBJrM\nW28BZ5xhdRREZDdrPXYyCRARtQA2DBMRUd2YBIiIHIxJgIjIwTaZvcBEIgG32414PI4777xz1fRY\nLAYAWFpawm233Wb26omIqA6mngmIoghRFBEIBKAoCo4cObJq+sDAAILBIBRFgSiKZq6ebGBhYcHq\nEGgdWH75zfUoAAAELElEQVTOY2oSCAQCuOOOOwAAmUwGPT09BdMVRUEymQQAeL1eKIpi5urJBngQ\naW4sP+cxvU0gl8shEolgz549q6YFg0EEg0EAgCRJ2L59u9mrN43ZX4a1Lq/W99UyX7V5yk2v93U7\nMDO2RpddrfNWmqfeaU4pu/Usz8zys/N3z/Qk0N7ejvHxcczMzCCdTpecR5Ik9PX1rTpTsJNm+yAy\nCRRiEqg8zSllt57lOSUJmDpYTJIkCIKA3t5eTE1NoaOjA+Pj46vmi0QiJV8HtAEPRERUv7Uczk3t\nHSSKInw+HwAgm81ix44dxmOXywUAiEajRgLQG5HzcbQwEdHGMbU6aPfu3VAUBbFYDG63G9deey0A\nYGBgAACQTCYxNTWF7u5ueDyeqr/60+k0RFFEJBIpW7VE9pXL5Yzyy+VyVodDazA2NmZ1CLQGiqJg\n165duOuuu6rOa2oSaG9vNxp/9+7da7yeSqUAaMkgk8ng6NGjyGQyOHToUMH75+bmIIqiMZZAkiT4\n/X4MDAwgkUiYGSqZYHJysuB5cfmlUin4/X64XC72BLOZamUHaN8/sqdq5ScIAmKxGG688caqy7Js\nxHA0GsXc3JzxXP/A6dVDsixjaGgI7e3tSCaTGBkZsSROKq2W8tMfZ7NZ9Pb2bnyQVFItZZdOp+F2\nu41qXLKPWsqvs7MTmUwGsVis6lm4ZUlg9+7d8Hq9xvPZ2Vm43W4A2hgCfTxBMpnEwMAAtm3bZkWY\nVEYt5ReLxdDe3g6fz4dIJGJVqFSklrJTFAXZbBaKorAq1mZqKb+5uTl0dnbC7/cjGo1WXJ5trh2U\nzWbh8XiM58vLyxBFEeFwGDMzMwWZj+ynVPlt374doihCURSeydlYqbILBALo7OxENptljz2bK1V+\nPp8PoigilUpV/e6Zfu2g9SjuGRQIBFb1HiL7Ki4/fRwIy9D+SvXKc7lcOHz4sAXRUL2Ky6+zsxOd\nnZ01ffdscybgcrmQyWQAACsrK+jo6LA4IqoHy695seya23rLzzZJYHR01OhBkk6nMTg4aHFEVA+W\nX/Ni2TW39ZafZUkgkUgglUoZ/Vj13iOiKMLlctn6khLE8mtmLLvmZnb52e4ew0REtHFsUx1EREQb\nj0mAiMjBmASIiByMSYCIyMFsNViMyK4kScLi4iI8Hg8ymQy8Xi8HwVFL4JkAUQ1mZ2cRDAYxNDQE\nAPD7/RZHRGQOJgGiKtLpNLLZrPHc7/ejvb3dwoiIzMNxAkQ16O7uhtfrRSgUMs4GiFoB2wSIanD0\n6FHMzc0ZN/NgIqBWwSRAVEU6nUZnZyeGhoZ4lzRqOWwTIKpAUZSC2yxKksQLrFFLYZsAUQX6zYz0\nS/V2dXXh8ssvtzIkIlMxCRARORirg4iIHIxJgIjIwZgEiIgcjEmAiMjBmASIiByMSYCIyMH+H9oC\nVRleG+aIAAAAAElFTkSuQmCC\n", "text": [ "" ] } ], "prompt_number": 7 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The MSE does not vary substantially as $N$ is scaled." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It would be nice if the above with cross-validation and much more data (also using Sklearn) (do this in parallel in the next notebook). Ask how this is done?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A full Cahn-Hilliard simulation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Fit the data again." ] }, { "cell_type": "code", "collapsed": false, "input": [ "model = MKSRegressionModel(Nbin=5)\n", "model.fit(X_train, y_train)" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 8 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Create a random sample to run a simulation." ] }, { "cell_type": "code", "collapsed": false, "input": [ "np.random.seed(11)\n", "X0 = np.random.random((1, 256, 256))" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 9 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Resize the coefficients to match the sample" ] }, { "cell_type": "code", "collapsed": false, "input": [ "model.resize_coeff(X0.shape[1:])" ], "language": "python", "metadata": {}, "outputs": [], "prompt_number": 10 }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run a simulation using the MKS!" ] }, { "cell_type": "code", "collapsed": false, "input": [ "X1 = X0.copy()\n", "data = []\n", "for step in range(1001):\n", " y = model.predict(X1)\n", " X1 = X1 + fipymodel.dt * y\n", " if step % 100 == 0:\n", " data.append(X1[0])\n", " print step," ], "language": "python", "metadata": {}, "outputs": [ { "output_type": "stream", "stream": "stdout", "text": [ "0 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "100 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "200 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "300 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "400 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "500 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "600 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "700 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "800 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "900 " ] }, { "output_type": "stream", "stream": "stdout", "text": [ "1000\n" ] } ], "prompt_number": 11 }, { "cell_type": "markdown", "metadata": {}, "source": [ "The above is too slow. It needs to be optimized. About 10000 steps are required for a full simulation." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Profiling in IPython (todo)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Profiling the `MKSRegressionClass`\n", "\n", "The `MKSRegerssionClass` is quite slow. The FFT is not even taking the majority of the time in `predict`. The next notebook will show how to speed up the `MKSRegressionClass` using Numexpr and parallel IPython." ] }, { "cell_type": "code", "collapsed": false, "input": [ "%load_ext line_profiler\n", "\n", "X1 = X0.copy()\n", "%lprun -f model.predict model.predict(X1)" ], "language": "python", "metadata": {}, "outputs": [] } ], "metadata": {} } ] }