Quadratic Regression 3D.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from deap import tools, creator, base, algorithms\n",
    "import numpy\n",
    "import random\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "%matplotlib inline\n",
    "\n",
    "points = []\n",
    "mean = 0.0\n",
    "model_coef = []\n",
    "\n",
    "#Rosenbrock z = (a-x)**2 + B(y-x**2)**2\n",
    "a = 1\n",
    "b = 100\n",
    "def rosenbrocke(x,y): \n",
    "    return (a-x)**2 + b*(y-x**2)**2\n",
    "\n",
    "def rosenbrock():\n",
    "    x = random.randint(-10,10)\n",
    "    y = random.randint(-10,10)\n",
    "    z = (a-x)**2 + b*(y-x**2)**2\n",
    "    return (x,y,z)\n",
    "\n",
    "def quadratic3d(coef):\n",
    "        x = point[0]\n",
    "        y = point[1]\n",
    "        return coef[0]*(x**2) + coef[1]*(y**2) + coef[2]*x + coef[3]*y + coef[4]\n",
    "    \n",
    "def quadratic3de(x,y):\n",
    "    return model_coef[0]*(x**2) + model_coef[1]*(y**2) + model_coef[2]*x + model_coef[3]*y + model_coef[4]\n",
    "\n",
    "def test_min3D(coords):\n",
    "    return (quadratic3de(coords[0],coords[1]),)\n",
    "\n",
    "def test_regress(coef):\n",
    "    SSres = 0 \n",
    "    SStot = 0\n",
    "\n",
    "    for point in points:        \n",
    "        zi = quadratic3d(coef)\n",
    "        SSres = SSres + (point[2] - zi)**2\n",
    "        SStot = SStot + (point[2] - mean)**2\n",
    "\n",
    "    err = SSres/SStot\n",
    "\n",
    "    return (err,)\n",
    "\n",
    "def model ():\n",
    "    creator.create(\"bestFit\", base.Fitness, weights = (-1.0,))\n",
    "    creator.create(\"Individual\",list,fitness=creator.bestFit)\n",
    "\n",
    "    toolbox = base.Toolbox()\n",
    "    toolbox.register(\"Coeff\",random.uniform,-30,30) #values picked at random, not sure how to improve\n",
    "    toolbox.register(\"individual\",tools.initRepeat,creator.Individual,toolbox.Coeff, 5) \n",
    "    toolbox.register(\"population\",tools.initRepeat,list,toolbox.individual, n = 50)\n",
    "    toolbox.register(\"evaluate\",test_regress)\n",
    "\n",
    "    toolbox.register(\"mate\",tools.cxBlend,alpha=.5)\n",
    "    toolbox.register(\"mutate\",tools.mutGaussian,mu=0,sigma=1,indpb=.5)\n",
    "    toolbox.register(\"select\", tools.selTournament,tournsize =3)\n",
    "    toolbox.register(\"map\",map)\n",
    "\n",
    "    cxprob, mutprob, ngens = 0.5,0.1,100\n",
    "    pop = toolbox.population()\n",
    "    ind = toolbox.individual()\n",
    "    hof = tools.HallOfFame(1)\n",
    "    pop = algorithms.eaMuCommaLambda(pop, toolbox,10,50, cxprob, mutprob, ngens,halloffame = hof, verbose = 0)\n",
    "    return hof[0]\n",
    "\n",
    "def findMin3D():\n",
    "    creator.create(\"bestFit\", base.Fitness, weights = (-1.0,))\n",
    "    creator.create(\"Individual\",list,fitness=creator.bestFit)\n",
    "\n",
    "    toolbox = base.Toolbox()\n",
    "    toolbox.register(\"coord\",random.uniform,-10,10) #same range as randomly generated\n",
    "    toolbox.register(\"individual\",tools.initRepeat,creator.Individual,toolbox.coord, 2) \n",
    "    toolbox.register(\"population\",tools.initRepeat,list,toolbox.individual, n = 50)\n",
    "    toolbox.register(\"evaluate\",test_min)\n",
    "\n",
    "    toolbox.register(\"mate\",tools.cxBlend,alpha=.5)\n",
    "    toolbox.register(\"mutate\",tools.mutGaussian,mu=0,sigma=1,indpb=.5)\n",
    "    toolbox.register(\"select\", tools.selTournament,tournsize =3)\n",
    "    toolbox.register(\"map\",map)\n",
    "\n",
    "    cxprob, mutprob, ngens = 0.5,0.1,100\n",
    "    pop = toolbox.population()\n",
    "    ind = toolbox.individual()\n",
    "    hof = tools.HallOfFame(1)\n",
    "    pop = algorithms.eaMuCommaLambda(pop, toolbox,10,50, cxprob, mutprob, ngens,halloffame = hof, verbose = 0)\n",
    "    return hof[0]\n",
    "\n",
    "#------------------------------\n",
    "random.seed(10)\n",
    "xi,yi,zi,mi = [],[],[],[] #empty lists for plotting\n",
    "fig = plt.figure()\n",
    "\n",
    "#getting initial data\n",
    "for i in range(10):\n",
    "    points.append(rosenbrock())\n",
    "\n",
    "#end setting up, begin main loop\n",
    "for i in range(1):    \n",
    "    #find mean for SStott,there might be a way to speed up recalculating but for now its not the speed problem\n",
    "    mean = 0\n",
    "    for point in points:\n",
    "        mean = mean + point[2]\n",
    "    mean = mean/len(points)\n",
    "\n",
    "    model_coef = model() #GA\n",
    "    #model_coef = MatrixRegress3D() \n",
    "    #model_coef = model_coef[0]\n",
    "\n",
    "    \n",
    "    newmin = findMin() #GA\n",
    "    print newmin\n",
    "\n",
    "    points.append([newmin[0],newmin[1],rosenbrocke(newmin[0],newmin[1])])\n",
    "        \n",
    "    #plotting\n",
    "    plt.clf()\n",
    "    xi,yi,zi,mi = [],[],[],[] #empty lists for plotting    \n",
    "    \n",
    "\n",
    "    for point in points:\n",
    "        xi.append(point[0])\n",
    "        yi.append(point[1]) \n",
    "        zi.append(point[2])\n",
    "        mi.append(quadratic3de(point[0],point[1]))\n",
    "\n",
    "\n",
    "    ax = fig.add_subplot(111,projection='3d')\n",
    "    ax.scatter(xi,yi,zi)\n",
    "    ax.plot_trisurf(xi,yi,mi,cmap=cm.spring)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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
}