Progress Gaussian Process forecast.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Database connection ok\n"
     ]
    }
   ],
   "source": [
    "import pyodbc\n",
    "from meyerDB import cable_connection\n",
    "import PyQt5\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib qt\n",
    "import numpy as np\n",
    "from IPython.display import HTML, display, clear_output\n",
    "\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.model_selection import LeavePOut\n",
    "from sklearn.metrics import mean_absolute_percentage_error as mape\n",
    "from sklearn.metrics import mean_absolute_error as mae, mean_squared_error as mse\n",
    "from sklearn import linear_model\n",
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "from sklearn.dummy import DummyRegressor\n",
    "import sklearn.preprocessing as pp\n",
    "from scipy.optimize import curve_fit\n",
    "from sklearn.gaussian_process import GaussianProcessRegressor\n",
    "from sklearn.gaussian_process.kernels import RBF, RBF2, WhiteKernel, ConstantKernel, RationalQuadratic\n",
    "from scipy.stats import norm\n",
    "from copy import deepcopy\n",
    "\n",
    "import math\n",
    "\n",
    "# init db connection\n",
    "conn = pyodbc.connect(cable_connection)\n",
    "cursor = conn.cursor()\n",
    "print('Database connection ok')\n",
    "\n",
    "def display_table(data):\n",
    "    html = \"<table>\"\n",
    "    for row in data:\n",
    "        html += \"<tr>\"\n",
    "        for field in row:\n",
    "            try:\n",
    "                value = str(round(100*field, 1)).replace('.', ',') + '%'\n",
    "            except:\n",
    "                value = field\n",
    "            html += \"<td><h4>%s</h4><td>\"%(value)\n",
    "        html += \"</tr>\"\n",
    "    html += \"</table>\"\n",
    "    display(HTML(html))\n",
    "\n",
    "def logifunc(x,x0,k,l, A):\n",
    "    #l = 2300\n",
    "    return l / (1 + A*np.exp(-k*(x - x0)))\n",
    "\n",
    "plt.rcParams['font.size'] = 14\n",
    "plt.rcParams['font.weight'] = 'bold'\n",
    "plt.rcParams[\"legend.loc\"] = 'upper right'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get ensemble quantity\n",
    "\n",
    "def get_ensemble_quantity(linear_pred, wk, ship):\n",
    "\n",
    "\n",
    "    pred_ends = np.load('pred_ends.npy', allow_pickle=True)\n",
    "    progresses = np.load('progress.npy', allow_pickle=True)\n",
    "\n",
    "\n",
    "    data = pred_ends.item().get(ship)\n",
    "    pred_end = data[np.argwhere(data[:, 0]==wk)][0][0][1]\n",
    "\n",
    "\n",
    "    data = progresses.item().get(ship)\n",
    "    cables = data[np.argwhere(data[:, 0]==wk)][0][0][1]\n",
    "\n",
    "    estimated_completeness = cables/linear_pred\n",
    "    if estimated_completeness < 0.002: return linear_pred\n",
    "\n",
    "    pred = cables + pred_end\n",
    "\n",
    "    estimated_completeness = np.square(cables/linear_pred)\n",
    "    if estimated_completeness > 1.0: estimated_completeness = 1.0\n",
    "    w1, w2 = 1-estimated_completeness, estimated_completeness\n",
    "    avpred = (w1*linear_pred + w2*pred)/(w1+w2)\n",
    "    return avpred\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\sainmi53\\AppData\\Local\\Programs\\Python\\Python39\\lib\\site-packages\\scipy\\optimize\\minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "C:\\Users\\sainmi53\\AppData\\Local\\Programs\\Python\\Python39\\lib\\site-packages\\scipy\\optimize\\minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "C:\\Users\\sainmi53\\AppData\\Local\\Programs\\Python\\Python39\\lib\\site-packages\\scipy\\optimize\\minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "C:\\Users\\sainmi53\\AppData\\Local\\Programs\\Python\\Python39\\lib\\site-packages\\scipy\\optimize\\minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n",
      "C:\\Users\\sainmi53\\AppData\\Local\\Programs\\Python\\Python39\\lib\\site-packages\\scipy\\optimize\\minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated\n",
      "  warnings.warn('Covariance of the parameters could not be estimated',\n"
     ]
    }
   ],
   "source": [
    "# Design progress\n",
    "# Predictors:\n",
    "# Time series, GT\n",
    "# Algorithms:\n",
    "# Gaussian process regression\n",
    "\n",
    "def moving_average(a, n=3) :\n",
    "    ret = np.cumsum(a, dtype=float)\n",
    "    ret[n:] = ret[n:] - ret[:-n]\n",
    "    return ret[n - 1:] / n\n",
    "def meanfunc(x):\n",
    "    params, _ = curve_fit(logifunc, fx, fy, bounds=([-100, 1e-3, 0.9999, 1e-4], [-10, 2, 1.0001, 100]))\n",
    "    return logifunc(x, *params)\n",
    "\n",
    "\n",
    "def calc_error():\n",
    "    res = y_mean[xidxs_inv] - y[xidxs_inv]\n",
    "    err = np.sqrt(np.square(res).mean())/y[-1]\n",
    "    return err\n",
    "\n",
    "\n",
    "def pplot():\n",
    "    plt.plot(x[xidxs], y[xidxs], 'g')\n",
    "    plt.plot(x[xidxs_inv], y[xidxs_inv], 'r')\n",
    "    plt.plot(x[xidxs_inv], y_mean[xidxs_inv], 'k--', alpha = 0.7)\n",
    "    plt.fill_between(x[xidxs_inv], (y_std[xidxs_inv] + y_mean[xidxs_inv]), (y_mean[xidxs_inv] - y_std[xidxs_inv]), alpha=0.4, color='k')\n",
    "    #plt.ylim([0, 4000000])\n",
    "    plt.xlim([-100, 0])\n",
    "    plt.pause(0.37)\n",
    "    plt.cla()\n",
    "\n",
    "# Get the data\n",
    "cursor.execute(\"SELECT r.project_id, gross_tonnage, sum(amount) FROM routed as r\"\n",
    "    \" LEFT JOIN projects as p ON p.project_id=r.project_id\"\n",
    "    \" GROUP BY r.project_id, gross_tonnage\")\n",
    "data = np.array(cursor.fetchall())\n",
    "ships = data[:, 0].astype('int32')\n",
    "gts = data[:, 1]\n",
    "cables = data[:, -1]\n",
    "# init lpo split\n",
    "p = 2 #ships\n",
    "lpo = LeavePOut(p)\n",
    "\n",
    "# init linear quantity\n",
    "model = linear_model.LinearRegression()\n",
    "\n",
    "wk = 'pw'\n",
    "cursor.execute(\"UPDATE progress SET cables=0\")\n",
    "cursor.execute(\n",
    "\n",
    "    \"UPDATE progress SET progress.cables=t1.cables FROM\"\n",
    "    \" (SELECT project_id, {}, sum(amount) as cables FROM ship_readiness GROUP BY project_id, {}) t1\"\n",
    "    \" WHERE t1.project_id=progress.project_id AND t1.{}=progress.wk\".format(wk, wk, wk)\n",
    ")\n",
    "cursor.commit()\n",
    "# get progress data\n",
    "ship_data2 = {}\n",
    "for i in range(ships.shape[0]):\n",
    "    ship = ships[i]\n",
    "    cursor.execute(\"SELECT wk, cables FROM progress WHERE project_id={} ORDER BY wk\".format(ship))\n",
    "    data = cursor.fetchall()\n",
    "    ship_data2[ship] = np.array(data)\n",
    "\n",
    "\n",
    "wk = 'rw'\n",
    "cursor.execute(\"UPDATE progress SET cables=0\")\n",
    "cursor.execute(\n",
    "\n",
    "    \"UPDATE progress SET progress.cables=t1.cables FROM\"\n",
    "    \" (SELECT project_id, {}, sum(amount) as cables FROM ship_readiness GROUP BY project_id, {}) t1\"\n",
    "    \" WHERE t1.project_id=progress.project_id AND t1.{}=progress.wk\".format(wk, wk, wk)\n",
    ")\n",
    "cursor.commit()\n",
    "# get progress data\n",
    "ship_data = {}\n",
    "for i in range(ships.shape[0]):\n",
    "    ship = ships[i]\n",
    "    cursor.execute(\"SELECT wk, cables FROM progress WHERE project_id={} AND cables > 100 ORDER BY wk\".format(ship))\n",
    "    data = cursor.fetchall()\n",
    "    ship_data[ship] = np.array(data)\n",
    "\n",
    "# init gaussian process regressor\n",
    "kernel = ConstantKernel(constant_value=0.1, constant_value_bounds='fixed') * RBF(\n",
    "    length_scale=0.0001, length_scale_bounds='fixed')\n",
    "gp = GaussianProcessRegressor(kernel=kernel,\n",
    "                            alpha=0.5)\n",
    "\n",
    "\n",
    "cv_errors = []\n",
    "all_errs = []\n",
    "completeness = []\n",
    "for train_i, test_i in lpo.split(ships):\n",
    "    model.fit(gts[train_i].reshape(-1, 1), cables[train_i])\n",
    "    fx = []\n",
    "    fy = []\n",
    "    for i in range(11):\n",
    "        if i in test_i: continue\n",
    "        ship = ships[i]\n",
    "        fx.append(ship_data[ship][:, 0])\n",
    "        fy.append(ship_data[ship][:, 1].cumsum()/cables[i])\n",
    "    fy = np.concatenate(fy, axis=0)\n",
    "    fx = np.concatenate(fx, axis=0)\n",
    "    fmean_y = []\n",
    "    fmean_x = []\n",
    "\n",
    "    n = 30\n",
    "    p = np.argsort(fx)\n",
    "    fx = fx[p]\n",
    "    fy = fy[p]\n",
    "    fx2 = deepcopy(fx)\n",
    "    fy2 = deepcopy(fy)\n",
    "    fmean_y = moving_average(fy, n)\n",
    "    fmean_x = moving_average(fx, n)\n",
    "\n",
    "    # search hyperparamters\n",
    "    lengths = np.linspace(0.01, 20, 30)\n",
    "    constants = np.linspace(0.001, 2.0, 30)\n",
    "    noises = np.linspace(0.01, 2.0, 30)\n",
    "    best_error = None\n",
    "    best_params = None\n",
    "    errors = []\n",
    "    for iii in range(0):\n",
    "        # init gaussian process regressor\n",
    "        l = np.random.choice(lengths)\n",
    "        c = np.random.choice(constants)\n",
    "        noise = np.random.choice(noises)\n",
    "        kernel = ConstantKernel(constant_value=c, constant_value_bounds='fixed') * RBF(\n",
    "            length_scale=l, length_scale_bounds='fixed')\n",
    "        gp = GaussianProcessRegressor(kernel=kernel,\n",
    "                                    alpha=noise)\n",
    "        for i in train_i:\n",
    "            ship = ships[i]\n",
    "            gt = gts[i]\n",
    "            x = ship_data[ship][:, 0]\n",
    "            for wk in np.arange(-100, 1):\n",
    "                wks = wk\n",
    "                y = ship_data[ship][:, 1].cumsum()\n",
    "                try:\n",
    "                    xidxs = np.argwhere(x <= wk)[:, 0]\n",
    "                    if xidxs.shape[0] == 0: continue\n",
    "                    xidxs_inv = np.argwhere(x > wk)[:, 0]\n",
    "                    if xidxs_inv.shape[0] == 0: continue\n",
    "                except IndexError:\n",
    "                    continue\n",
    "                real_completeness = y[xidxs][-1]/y[-1]\n",
    "                if real_completeness < 0.001: continue\n",
    "                if real_completeness > 0.999: continue\n",
    "                \n",
    "                X = x[xidxs].reshape(-1, 1)\n",
    "                mf = meanfunc(x)\n",
    "                Y = y[xidxs] - mf[xidxs]\n",
    "                gp.fit(X, Y)\n",
    "                y_mean, y_cov = gp.predict(x.reshape(-1, 1), return_cov=True)\n",
    "                y_mean = y_mean + mf\n",
    "\n",
    "                #clear_output(wait=True)\n",
    "                err = calc_error()\n",
    "                errors.append(err)\n",
    "        errors = np.array(errors)\n",
    "        err = errors.mean()\n",
    "        if best_error is None or best_error > err:\n",
    "            best_error = err\n",
    "            best_params = (l, c, noise)\n",
    "        errors = []\n",
    "    (l, c, noise) = (10.0, 0.01, 0.001)#best_params #(10.0, 1.0, 0.1)\n",
    "    #print(best_params)\n",
    "    kernel = ConstantKernel(constant_value=c, constant_value_bounds='fixed') * RBF(\n",
    "        length_scale=l, length_scale_bounds='fixed')\n",
    "    gp = GaussianProcessRegressor(kernel=kernel,\n",
    "                                alpha=noise)\n",
    "\n",
    "    # test model\n",
    "    for i in test_i:\n",
    "        ship = ships[i]\n",
    "        gt = gts[np.argwhere(ships==ship)][0]\n",
    "        linear_pred = model.predict(gt.reshape(1, -1))[0]\n",
    "        x = ship_data[ship][:, 0]\n",
    "        x2 = ship_data2[ship][:, 0]\n",
    "        for wk in np.arange(-100, 1):\n",
    "            pred = get_ensemble_quantity(linear_pred, wk, ship)\n",
    "            #pred = cables[i]\n",
    "            y = ship_data[ship][:, 1].cumsum()/pred\n",
    "            y2 = ship_data2[ship][:, 1].cumsum()/pred\n",
    "            try:\n",
    "                xidxs = np.argwhere(x <= wk)[:, 0]\n",
    "                xidxs2 = np.argwhere(x2 <= wk)[:, 0]\n",
    "                if xidxs.shape[0] == 0: continue\n",
    "                xidxs_inv = np.argwhere(x > wk)[:, 0]\n",
    "                if xidxs_inv.shape[0] == 0: continue\n",
    "            except IndexError:\n",
    "                continue\n",
    "            real_completeness = y[xidxs][-1]/y[-1]\n",
    "            if real_completeness < 0.001: continue\n",
    "            if real_completeness > 0.999: continue\n",
    "            \n",
    "            X = x[xidxs].reshape(-1, 1)\n",
    "            if True: # use design process as reference\n",
    "                fx = x2[xidxs2]\n",
    "                fy = y2[xidxs2]\n",
    "            mf = meanfunc(x)\n",
    "            Y = y[xidxs] - mf[xidxs]\n",
    "            gp.fit(X, Y)\n",
    "            y_mean, y_std = gp.predict(x.reshape(-1, 1), return_std=True)\n",
    "            y_std = (1.96 * (y_std + noise*np.ones_like(y_std))) * pred\n",
    "            \n",
    "            y_mean = (y_mean + mf) * pred\n",
    "             # limit\n",
    "            idxs = np.argwhere(y_mean > mf*pred)\n",
    "            y_mean[idxs] = mf[idxs]*pred\n",
    "            y = y*pred\n",
    "\n",
    "            #clear_output(wait=True)\n",
    "            err = calc_error()\n",
    "            errors.append(err)\n",
    "            all_errs.append(err)\n",
    "            completeness.append(real_completeness)\n",
    "            pplot()\n",
    "\n",
    "\n",
    "    errors = np.array(errors)\n",
    "    cv_errors.append(errors.mean())\n",
    "          \n",
    "cv_errors = np.array(cv_errors)            \n",
    "print(cv_errors.mean(), cv_errors.std(), cv_errors.max())\n",
    "\n",
    "\n",
    "all_errs = np.array(all_errs)\n",
    "completeness = np.array(completeness)\n",
    "p = np.argsort(completeness)\n",
    "completeness = completeness[p]\n",
    "all_errs = all_errs[p]\n",
    "print(all_errs.shape)\n",
    "n = 170\n",
    "all_errs = moving_average(all_errs, n)\n",
    "\n",
    "\n",
    "plt.plot(100*completeness[n-1:], 100*all_errs, 'r')\n",
    "plt.legend(['Error - Moving average (100)'])\n",
    "#plt.title(\"S-Curve fitting with ensemble quantity prediction\")\n",
    "plt.xlabel(\"Completeness of cabling process (%)\")\n",
    "plt.ylabel('Error (weeks)')\n",
    "plt.ylim(-10, 100)\n",
    "plt.grid(which='major')\n",
    "#plt.savefig('wkdesigngpensemble.png')\n",
    "\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "c0b19f2120ef3124c13ef286ebe522ac06f45b032087a286eca843495c1a56b0"
  },
  "kernelspec": {
   "display_name": "Python 3.9.4 64-bit (system)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.4"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}