See AWS for Paper Results 4d citibike.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4d Citibike Experimental Results for Paper\n",
    "\n",
    "We've moved this to run on AWS instances"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "warning in stationary: failed to import cython module: falling back to numpy\n"
     ]
    }
   ],
   "source": [
    "import dp4gp_datasets\n",
    "import dp4gp\n",
    "import random\n",
    "import numpy as np\n",
    "import GPy\n",
    "import matplotlib.pyplot as plt\n",
    "import dp4gp_histogram\n",
    "import pandas as pd\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## House prices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def get_house_prices():\n",
    "    #Skip this if you want, and load precomputed data from the csv file below:\n",
    "    #dp4gp_datasets.prepare_preloaded_prices('test.csv', boundingbox=[480e3, 130e3, 580e3, 230e3], N=10000, col_list=['QS501EW'])\n",
    "    #Load precomputed dataset:\n",
    "    dataset = pd.read_csv('price_dataset10k.csv') #london, 10k purchases\n",
    "\n",
    "    #Reduce the size of the dataset and split into training and test data:\n",
    "    random.seed(123)\n",
    "    #dataset = dataset[dataset['years']>2005]\n",
    "    dataset = dataset.ix[random.sample(dataset.index, 200)]\n",
    "\n",
    "    #get into useful form\n",
    "    #east, north, time, education -> price\n",
    "    inputs = np.vstack([dataset['easting'].values,dataset['northing'].values,dataset['seconds'].values,(dataset['QS501EW_6']/dataset['QS501EW_0']).values]).T\n",
    "\n",
    "    #inputs = inputs[:,0:2]\n",
    "    ys = dataset['price'].values\n",
    "    \n",
    "    return inputs, ys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_citibike_data():\n",
    "    allcb = dp4gp_datasets.load_citibike(station=None)\n",
    "    subcb = allcb[allcb['usertype']=='Subscriber']\n",
    "    random.seed(123)\n",
    "\n",
    "    cb = subcb.ix[random.sample(subcb.index, 5000)]\n",
    "    inputs = np.c_[cb['start station latitude'],cb['end station latitude'],cb['start station longitude'],cb['end station longitude']]\n",
    "    ys = cb['tripduration'].values\n",
    "    \n",
    "    return inputs, ys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "def get_noDP_prediction(training_inputs, training_ys, test_inputs, sens, eps, delta,noise,modvar,kernval,kern_ls,steps):\n",
    "    rbf = GPy.kern.RBF(training_inputs.shape[1],kernvar, kern_ls,ARD=True)\n",
    "    mod = GPy.models.GPRegression(training_inputs,training_ys,rbf)\n",
    "    mod.Gaussian_noise = noise\n",
    "    dpgp = dp4gp.DPGP_cloaking(mod,sens,eps,delta)\n",
    "    preds, mu, cov = dpgp.draw_prediction_samples(test_inputs,1,1,0)\n",
    "    return mu, None\n",
    "\n",
    "\n",
    "def get_cloaking_prediction(training_inputs, training_ys, test_inputs, sens, eps, delta,noise,modvar,kernval,kern_ls,steps):\n",
    "    rbf = GPy.kern.RBF(training_inputs.shape[1],kernvar, kern_ls,ARD=True)\n",
    "    #rbf = GPy.kern.RBF(training_inputs.shape[1],modvar,[5e3,5e3],ARD=True)\n",
    "    mod = GPy.models.GPRegression(training_inputs,training_ys,rbf)\n",
    "    mod.Gaussian_noise = noise\n",
    "    dpgp = dp4gp.DPGP_cloaking(mod,sens,eps,delta)\n",
    "    preds, mu, cov = dpgp.draw_prediction_samples(test_inputs,1,1,1000)\n",
    "    return preds, cov\n",
    "    \n",
    "def get_integral_prediction(training_inputs, training_ys, test_inputs, sens, eps, delta,noise,modvar,kernval,kern_ls,steps):\n",
    "    Xtest, free_inputs, step = dp4gp.compute_Xtest(training_inputs,steps=steps)\n",
    "    print step\n",
    "    dpgp = dp4gp_histogram.DPGP_histogram(sens,eps,delta)\n",
    "    dpgp.prepare_model(Xtest,training_inputs,step,training_ys,lengthscale=kern_ls)\n",
    "    #dpgp.optimize()\n",
    "    dpgp.model.optimize(messages=True)\n",
    "    preds, cov = dpgp.draw_prediction_samples(test_inputs)\n",
    "    return preds, cov\n",
    "\n",
    "def get_standard_prediction(training_inputs, training_ys, test_inputs, sens, eps, delta,noise,modvar,kernval,kern_ls,steps):\n",
    "    rbf = GPy.kern.RBF(training_inputs.shape[1],kernvar,kern_ls,ARD=True)\n",
    "    mod = GPy.models.GPRegression(training_inputs,training_ys,rbf)\n",
    "    mod.Gaussian_noise = noise\n",
    "    dpgp = dp4gp.DPGP_normal_prior(mod,sens,eps,delta)\n",
    "    preds, mu, cov = dpgp.draw_prediction_samples(test_inputs,1)\n",
    "    return preds, cov\n",
    "\n",
    "def get_pseudo_prediction(training_inputs, training_ys, test_inputs, sens, eps, delta,noise,modvar,kernval,kern_ls,steps):\n",
    "    rbf = GPy.kern.RBF(training_inputs.shape[1],kernvar,kern_ls,ARD=True)\n",
    "    print training_ys.ndim\n",
    "    print training_inputs.shape\n",
    "    mod = GPy.models.SparseGPRegression(training_inputs,training_ys,kernel=rbf,num_inducing=40) #no idea how many inducing!\n",
    "    mod.inference_method = GPy.inference.latent_function_inference.FITC()\n",
    "    mod.set_Z(training_inputs[0:40,:]) #grab random inputs as pseudoinputs\n",
    "    mod.Gaussian_noise = noise  \n",
    "    dpgp = dp4gp.DPGP_pseudo_prior(mod,sens,eps,delta)\n",
    "    preds, mu, cov = dpgp.draw_prediction_samples(test_inputs,1)\n",
    "    return preds, cov"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(100, 4900)\n",
      "*\n",
      "Stopped before convergence\n",
      "(3.9836146099035994, 3.2552472614374586, array([[ 1.26196074]]), 0.01, 0.0)\n",
      "[ 0.12933827  0.12996228  0.10469202  0.15097224]\n",
      "(100, 4900)\n",
      "*\n",
      ". . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ."
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-5-655fc6b0da55>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m\n\u001b[0;32m     36\u001b[0m         \u001b[0mresults\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     37\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mfn\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlabel\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mzip\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfns\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mlabels\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 38\u001b[1;33m             \u001b[0mpreds\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcov\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfn\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtraining_inputs\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mtraining_ys\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mtest_inputs\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msens\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0meps\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m0.01\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m10.0\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1.0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkernvar\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkern_ls\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msteps\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     39\u001b[0m             \u001b[0mRMSE\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msqrt\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmean\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtest_ys\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mpreds\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m**\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     40\u001b[0m             \u001b[0mresults\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m{\u001b[0m\u001b[1;34m'label'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mlabel\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'preds'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mpreds\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'cov'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mcov\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'RMSE'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[0mRMSE\u001b[0m\u001b[1;33m}\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-4-98466bdd8b53>\u001b[0m in \u001b[0;36mget_cloaking_prediction\u001b[1;34m(training_inputs, training_ys, test_inputs, sens, eps, delta, noise, modvar, kernval, kern_ls, steps)\u001b[0m\n\u001b[0;32m     15\u001b[0m     \u001b[0mmod\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mGaussian_noise\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnoise\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     16\u001b[0m     \u001b[0mdpgp\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdp4gp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mDPGP_cloaking\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmod\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0msens\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0meps\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mdelta\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 17\u001b[1;33m     \u001b[0mpreds\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmu\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcov\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdpgp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdraw_prediction_samples\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtest_inputs\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1000\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     18\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mpreds\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcov\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     19\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/lionfish/Documents/Research/dp/dp4gp.pyc\u001b[0m in \u001b[0;36mdraw_prediction_samples\u001b[1;34m(self, Xtest, N, Nattempts, Nits)\u001b[0m\n\u001b[0;32m     96\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mdraw_prediction_samples\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mXtest\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mN\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mNattempts\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m7\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mNits\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1000\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     97\u001b[0m         \u001b[0mGPymean\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcovar\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmodel\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpredict\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mXtest\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 98\u001b[1;33m         \u001b[0mmean\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnoise\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcov\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdraw_noise_samples\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mXtest\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mN\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mNattempts\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mNits\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     99\u001b[0m         \u001b[1;31m#TODO: In the long run, remove DP4GP's prediction code and just use GPy's\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    100\u001b[0m         \u001b[1;31m#print GPymean-mean\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/lionfish/Documents/Research/dp/dp4gp.pyc\u001b[0m in \u001b[0;36mdraw_noise_samples\u001b[1;34m(self, Xtest, N, Nattempts, Nits)\u001b[0m\n\u001b[0;32m    389\u001b[0m             \u001b[0mcs\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mC\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    390\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 391\u001b[1;33m         \u001b[0mls\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfindLambdas_repeat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcs\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mNattempts\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mNits\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    392\u001b[0m         \u001b[0mM\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcalcM\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mls\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    393\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/lionfish/Documents/Research/dp/dp4gp.pyc\u001b[0m in \u001b[0;36mfindLambdas_repeat\u001b[1;34m(self, cs, Nattempts, Nits)\u001b[0m\n\u001b[0;32m    325\u001b[0m             \u001b[0msys\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstdout\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mflush\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    326\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 327\u001b[1;33m             \u001b[0mls\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfindLambdas_grad\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcs\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mNits\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    328\u001b[0m             \u001b[1;32mif\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mls\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m<\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m0.01\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    329\u001b[0m                 \u001b[1;32mcontinue\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/lionfish/Documents/Research/dp/dp4gp.pyc\u001b[0m in \u001b[0;36mfindLambdas_grad\u001b[1;34m(self, cs, maxit)\u001b[0m\n\u001b[0;32m    289\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mit\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mmaxit\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    290\u001b[0m             \u001b[0mlsbefore\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mls\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcopy\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 291\u001b[1;33m             \u001b[0mdelta_ls\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m-\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdL_dl\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mls\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mlr\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    292\u001b[0m             \u001b[0mls\u001b[0m \u001b[1;33m=\u001b[0m  \u001b[0mls\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mdelta_ls\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    293\u001b[0m             \u001b[0mls\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mls\u001b[0m\u001b[1;33m<\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/lionfish/Documents/Research/dp/dp4gp.pyc\u001b[0m in \u001b[0;36mdL_dl\u001b[1;34m(self, ls, cs)\u001b[0m\n\u001b[0;32m    271\u001b[0m         \u001b[0mgrads\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mls\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    272\u001b[0m         \u001b[1;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 273\u001b[1;33m             \u001b[0mgrads\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m-\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtrace\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mMinv\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mcs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mcs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    274\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgrads\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    275\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m/home/lionfish/anaconda2/lib/python2.7/site-packages/numpy/core/fromnumeric.pyc\u001b[0m in \u001b[0;36mtrace\u001b[1;34m(a, offset, axis1, axis2, dtype, out)\u001b[0m\n\u001b[0;32m   1378\u001b[0m         \u001b[1;32mreturn\u001b[0m \u001b[0masarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtrace\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0moffset\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1379\u001b[0m     \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1380\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0masanyarray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mtrace\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0moffset\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis1\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis2\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1381\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1382\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "#inputs, ys = get_house_prices()\n",
    "inputs, ys = get_citibike_data()\n",
    "\n",
    "#FOR CITIBIKE\n",
    "#squash data into 0-2000 seconds range\n",
    "ys[ys>2000] = 2000\n",
    "ys[ys<0] = 0\n",
    "sens = 2000-0\n",
    "kernvar = 10.0\n",
    "kern_ls = np.array([0.01,0.01,0.01,0.01])*5.0\n",
    "\n",
    "#FOR HOUSE PRICES\n",
    "#squash data into £10k-£500k range\n",
    "#ys[ys>5e5] = 5e5\n",
    "#ys[ys<1e4] = 1e4\n",
    "#sens = 5e5-1e4\n",
    "#kernvar = 1.0\n",
    "#kern_ls = [15e3,15e3,50*31536000,5.0]\n",
    "\n",
    "ys_mean = np.mean(ys)\n",
    "ys_std = np.std(ys)\n",
    "ys = ys - ys_mean\n",
    "ys = ys / ys_std\n",
    "sens = sens / ys_std\n",
    "\n",
    "training_inputs = inputs[0:-100,:]\n",
    "training_ys = ys[0:-100][:,None]\n",
    "test_inputs = inputs[-100:,:]\n",
    "test_ys = ys[-100:][:,None]\n",
    "\n",
    "fns = [get_noDP_prediction, get_integral_prediction, get_cloaking_prediction]#,get_pseudo_prediction,get_standard_prediction]\n",
    "labels = [\"No DP\",\"Integral\",\"Cloaking\"]   #\"Pseudo\",\"Standard\"]\n",
    "\n",
    "for steps in [1,3,6,10]:\n",
    "    for eps in [0.01,0.1,0.2,0.5,1.0]:\n",
    "        results = []\n",
    "        for fn,label in zip(fns,labels):\n",
    "            preds, cov = fn(training_inputs,training_ys,test_inputs,sens,eps,0.01,10.0,1.0, kernvar, kern_ls, steps)\n",
    "            RMSE = np.sqrt(np.mean((test_ys-preds)**2))\n",
    "            results.append({'label':label, 'preds':preds, 'cov':cov, 'RMSE':RMSE})\n",
    "        print \"EPSILON = %0.3f\" % eps\n",
    "        print \"STEPS = %0.3f\" % steps\n",
    "        for r in results:\n",
    "            print r['label'],r['RMSE']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for r in results:\n",
    "    print r['label'],r['RMSE']#*ys_std"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=[14,10])\n",
    "for i,r in enumerate(results):\n",
    "    \n",
    "    plt.subplot(2,2,i+1)\n",
    "    plt.title(r['label'])\n",
    "    plt.errorbar(test_ys,r['preds'],1.96*np.sqrt(np.diag(r['cov'])),fmt=None,color='k',alpha=0.2)\n",
    "    plt.plot(test_ys,r['preds'],'+k')\n",
    "    #plt.xlim([0,1e6])\n",
    "    #plt.ylim([0,1e6])\n",
    "    #plt.plot([0,1e6],[0,1e6],'k-')\n",
    "    plt.xlim([-2,4])\n",
    "    plt.ylim([-2,4])\n",
    "    plt.plot([-2,4],[-2,4],'k-')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def plot_cloaking_prediction(training_inputs, training_ys, test_inputs, sens, eps, delta,noise,modvar):\n",
    "    rbf = GPy.kern.RBF(training_inputs.shape[1],1.0,[5e3,5e3,50*31536000,5.0],ARD=True)\n",
    "    #rbf = GPy.kern.RBF(training_inputs.shape[1],modvar,[5e3,5e3],ARD=True)\n",
    "    mod = GPy.models.GPRegression(training_inputs,training_ys,rbf)\n",
    "    mod.Gaussian_noise = noise\n",
    "    dpgp = dp4gp.DPGP_cloaking(mod,sens,eps,delta)\n",
    "    return dpgp\n",
    "dpgp = plot_cloaking_prediction(training_inputs,training_ys[:,None],test_inputs,sens,1.0,0.01,5.0,1.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dpgp.plot(fixed_inputs=[(2,np.mean(training_inputs[:,2])),(3,np.mean(training_inputs[:,3]))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.mgrid[0:2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.mgrid[0:2,0:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "np.mgrid[0:2,0:3,0:4][0].flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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": 0
}