rush-models-notebook.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from collections import defaultdict, Counter\n",
    "import copy\n",
    "import pandas as pd\n",
    "import sys\n",
    "import datetime as dt\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy\n",
    "import statsmodels.api as sm\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SimpleMultivariateHawkesProcess:\n",
    "    _ndim = None # no. of dimensions\n",
    "    _lambda = None # ndim x 1 vector\n",
    "    _alpha = None # ndim x ndim matrix\n",
    "    _beta = None # ndim x ndim matrix\n",
    "    _observation_window = None\n",
    "    _EPS = 10**-9 # numbers below this are set to 0\n",
    "\n",
    "    def __init__(self, ndim=None, lam=None, alpha=None, beta=None, observation_window=None):\n",
    "        self._ndim = ndim\n",
    "        self._lambda = lam\n",
    "        self._alpha = alpha\n",
    "        self._beta = beta\n",
    "        self._observation_window = observation_window\n",
    "    \n",
    "    def fit(self, ndim, alpha, data_epochs, data_dims, observation_window, niter=1000):\n",
    "        print 'Fitting...'\n",
    "        print '\\tPrecomputing A(i), B(i)...'\n",
    "        #A, B = self._par_precompute(ndim, alpha, data_epochs, data_dims, observation_window)\n",
    "        A, B = self._precompute(ndim, alpha, data_epochs, data_dims, observation_window)\n",
    "                \n",
    "        print '\\tMinimizing negative loglikelihood...'\n",
    "        alpha_ravel = alpha.ravel()\n",
    "        x0 = np.random.uniform(low=10**-5.0, high=alpha_ravel[0] - self._EPS, size=(ndim * (ndim + 1)))\n",
    "        bounds = [(self._EPS, None) for i in range(ndim)] +\\\n",
    "                 [(self._EPS, None) for i in range(ndim * ndim)] \n",
    "                 #[(self._EPS, None) for i in range(ndim * ndim)]\n",
    "                            \n",
    "        x, f, d = scipy.optimize.fmin_l_bfgs_b(self._seq_negloglikelihood,\n",
    "                                               x0=x0, args=(ndim, alpha, observation_window, A, B),\n",
    "                                               factr=10.0**7.0, iprint=99, pgtol=10**-5, approx_grad=False,\n",
    "                                               callback=None, disp=99,\n",
    "                                               bounds=bounds, maxiter=niter, maxfun=niter)\n",
    "\n",
    "        return x, f, d\n",
    "    \n",
    "    @staticmethod\n",
    "    def _precompute(ndim, alpha, data_epochs, data_dims, observation_window):\n",
    "        T0, T = observation_window\n",
    "        nsamples = len(data_epochs)\n",
    "\n",
    "        A = []\n",
    "        B = [np.zeros(shape=(ndim, ndim), dtype=np.float64) for k in range(nsamples)]\n",
    "        for k in range(nsamples):\n",
    "            # memoize locations of timestamps for each dimension for this sample (2s for 18K samples)\n",
    "            epochs_per_dimension = []\n",
    "            for i in range(ndim):\n",
    "                idx_i = data_dims[k] == i\n",
    "                epochs_i = data_epochs[k][idx_i]\n",
    "                epochs_per_dimension.append(epochs_i)\n",
    "            \n",
    "            A_k = [np.zeros(shape=(ndim, len(epochs_per_dimension[m])), dtype=np.float64)\n",
    "                   for m in range(ndim)]\n",
    "            #B_mn = np.zeros(shape=(ndim, ndim), dtype=np.float64)\n",
    "            for m in range(ndim):\n",
    "                epochs_m = epochs_per_dimension[m]\n",
    "                if len(epochs_m) == 0:\n",
    "                    continue # no transactions on dimension m\n",
    "                #A_m = np.zeros(shape=(ndim, len(epochs_m)), dtype=np.float64)\n",
    "                for n in range(ndim):\n",
    "                    epochs_n = epochs_per_dimension[n]\n",
    "                    if len(epochs_n) == 0: # no transactions on dimension n\n",
    "                        continue\n",
    "                    len_n = len(epochs_n)\n",
    "                    A_mn = A_k[m][n] # influence of dimension n on each timestamp of dimension m\n",
    "\n",
    "                    j = 0\n",
    "                    if (m != n): # mutual excitation    \n",
    "                        \n",
    "                        while j < len_n and epochs_n[j] < epochs_m[0]:\n",
    "                            j += 1\n",
    "                            \n",
    "                        #A_mn[0] = np.sum(np.exp(-alpha[m,n] * (epochs_m[0] - epochs_n[epochs_n < epochs_m[0]]) / 1000.0))\n",
    "                        A_mn[0] = np.sum(np.exp(-alpha[m,n] * (epochs_m[0] - epochs_n[:j]) / 1000.0))\n",
    "\n",
    "                        for i in range(1, len(epochs_m)):\n",
    "                            prev_j = j\n",
    "                            while j < len_n and epochs_n[j] < epochs_m[i]:\n",
    "                                j += 1\n",
    "                            A_mn[i] = np.exp(-alpha[m,n] * (epochs_m[i] - epochs_m[i-1]) / 1000.0) * A_mn[i-1]\n",
    "                            #A_mn[i] += np.sum(np.exp(-alpha[m,n] *\n",
    "                            #                         (epochs_m[i] - epochs_n[(epochs_n >= epochs_m[i-1]) &\n",
    "                            #                                                 (epochs_n < epochs_m[i])]) / 1000.0))\n",
    "                            A_mn[i] += np.sum(np.exp(-alpha[m,n] *\n",
    "                                                     (epochs_m[i] - epochs_n[prev_j:j]) / 1000.0))\n",
    "                        \n",
    "                        B[k][m,n] = np.sum(1 - np.exp(-alpha[m,n] * (T - epochs_n) / 1000.0))\n",
    "                    else: # dim_m = dim_n, self-excitation\n",
    "                        for i in range(1, len(epochs_m)):\n",
    "                            A_mn[i] = np.exp(-alpha[m,n] * (epochs_m[i] - epochs_m[i-1]) / 1000.0) * (1 + A_mn[i-1])\n",
    "                        B[k][m,n] = np.sum(1 - np.exp(-alpha[m,n] * (T - epochs_m) / 1000.0))\n",
    "                \n",
    "                #A_k.append(A_m)\n",
    "            \n",
    "            A.append(A_k)\n",
    "            #B[k] = B_mn\n",
    "\n",
    "        return A, B\n",
    "    \n",
    "    def _seq_negloglikelihood(self, params, ndim, alpha, observation_window, A, B):\n",
    "        nsamples = len(A)\n",
    "        nll = 0.0\n",
    "        grad = np.zeros(shape=((ndim + 1) * ndim), dtype=np.float64) # first row is grad of lambda\n",
    "        for k in range(nsamples):\n",
    "            sample_nll, sample_grad = self._negloglikelihood(params, ndim, alpha, observation_window, A[k], B[k])\n",
    "            nll += sample_nll\n",
    "            grad += sample_grad\n",
    "        nll /= nsamples\n",
    "        grad /= nsamples\n",
    "        return nll, grad\n",
    "    \n",
    "    @staticmethod\n",
    "    def _negloglikelihood(params, ndim, alpha, observation_window, A, B):\n",
    "        lam = params[:ndim] # first ndim elements\n",
    "        beta = params[ndim:].reshape((ndim, ndim)) # next ndim * ndim elements\n",
    "        T0, T = observation_window\n",
    "        T_length = (T - T0)/1000.0 # scale down\n",
    "        \n",
    "        loglikelihood = 0.0\n",
    "        grad_lambda = np.zeros(shape=(ndim), dtype=np.float64)\n",
    "        grad_beta = np.zeros(shape=(ndim, ndim), dtype=np.float64)\n",
    "\n",
    "        for m in range(ndim):\n",
    "            Z = lam[m] + np.sum(beta[m,:][:,None] * A[m], axis=0)\n",
    "            loglikelihood += np.sum(np.log(Z)) - np.sum(beta[m,:] * B[m,:] / alpha[m,:]) - lam[m] * T_length\n",
    "            grad_lambda[m] = np.sum(1.0/Z) - T_length\n",
    "            grad_beta[m] = np.sum(A[m]/Z, axis=1) - B[m,:]/alpha[m,:]\n",
    "            \n",
    "            #print np.sum(np.log(Z)) - np.sum(beta[m,:] * B[m,:] / alpha[m,:]) - lam[m] * T_length \\\n",
    "            #   - len(A[m][0]) * np.log(1000) # compare with Nan\n",
    "            #print 1000.0 * np.sum(1.0/Z) - T_length * 1000.0 # compare with Nan\n",
    "            #print 1000.0 * grad_beta[m] # compare with nan\n",
    "        \n",
    "        gradient = np.concatenate((grad_lambda, grad_beta.ravel()))\n",
    "        return -loglikelihood, -gradient\n",
    "\n",
    "    def check_gradient(self, ndim, alpha, data_epochs, data_dims, observation_window, delta=10**-9, niter=50):\n",
    "        print 'Checking gradient...'\n",
    "        print '\\tPrecomputing A, B...'\n",
    "        A, B = self._precompute(ndim, alpha, data_epochs, data_dims, observation_window)\n",
    "        \n",
    "        nsamples = len(data_epochs)\n",
    "        for i in range(niter):\n",
    "            row_idx = np.random.randint(0, nsamples)\n",
    "            a = A[row_idx]\n",
    "            b = B[row_idx]\n",
    "            x0 = np.array([np.random.uniform(10**-3, 10**3), np.random.uniform(10**-3, 10**3),  # lambda\n",
    "                           np.random.uniform(10**-3, 10**3), np.random.uniform(10**-3, 10**3),  # beta\n",
    "                           np.random.uniform(10**-3, 10**3), np.random.uniform(10**-3, 10**3)]) # beta\n",
    "\n",
    "            nll, grad = self._negloglikelihood(x0, ndim, alpha, observation_window, a, b)\n",
    "            param_idx = np.random.randint(len(x0))\n",
    "            x0[param_idx] += delta\n",
    "            new_nll, new_grad = self._negloglikelihood(x0, ndim, alpha, observation_window, a, b)\n",
    "            num_grad = (new_nll - nll)/delta            \n",
    "            if abs((grad[param_idx] - num_grad)/num_grad)  > 0.01:\n",
    "                print 'Gradient checking failed for parameter:', param_idx\n",
    "                print 'Analytical gradient:', grad[param_idx],\n",
    "                print 'Numerical gradient:', num_grad\n",
    "        print '\\tGradient checking complete.'\n",
    "        \n",
    "    def simulate(self, t0, prev_t, prev_d, horizon_duration):\n",
    "        simulated_iats = []\n",
    "        simulated_dims = []\n",
    "        prev_t_copy = copy.copy(list(prev_t))\n",
    "        prev_d_copy = copy.copy(list(prev_d))\n",
    "        prev_t_array = np.array(prev_t_copy, dtype=np.float64)\n",
    "        prev_d_array = np.array(prev_d_copy, dtype=np.int)\n",
    "        cumulative_iat = 0.0\n",
    "        while cumulative_iat < horizon_duration:\n",
    "            upper_bound = np.sum(self.intensity(t0 + cumulative_iat, (prev_t_array, prev_d_array), include_t0=True))\n",
    "            s = 1000.0 * np.random.exponential(1.0/upper_bound) # params are in 1/1000s units            \n",
    "                \n",
    "            cumulative_iat += s # move ahead in time (may accept or reject)\n",
    "            if cumulative_iat >= horizon_duration:\n",
    "                break\n",
    "\n",
    "            lambda_t_dim = self.intensity(t0 + cumulative_iat, (prev_t_array, prev_d_array))\n",
    "            lambda_t = np.sum(lambda_t_dim)\n",
    "\n",
    "            u = np.random.uniform(0.0, 1.0)\n",
    "            if u <= lambda_t/upper_bound: # accept\n",
    "                prev_t_copy.append(t0 + cumulative_iat)\n",
    "                prev_t_array = np.array(prev_t_copy, dtype=np.float64)\n",
    "                simulated_iats.append(cumulative_iat)\n",
    "                \n",
    "                cumprob = np.cumsum(lambda_t_dim/lambda_t)\n",
    "                u2 = np.random.uniform(0.0, 1.0)\n",
    "                dim = np.searchsorted(cumprob, u2)\n",
    "                \n",
    "                prev_d_copy.append(dim)\n",
    "                prev_d_array = np.array(prev_d_copy, dtype=np.int)\n",
    "                simulated_dims.append(dim)\n",
    "\n",
    "        return simulated_iats, simulated_dims\n",
    "    \n",
    "    def predict_dist(self, t0, prev_t, prev_d, niter=1000):\n",
    "        #step = 3600.0 - (t0 % 3600)\n",
    "        simulated_iats = np.zeros(shape=(niter), dtype=np.float64)\n",
    "        simulated_dims = np.zeros(shape=(niter), dtype=np.int)\n",
    "        for i in range(niter):\n",
    "            cumulative_iat = 0.0\n",
    "            while True:\n",
    "                upper_bound = np.sum(self.intensity(t0 + cumulative_iat, (prev_t, prev_d), include_t0=True))\n",
    "                s = 1000.0 * np.random.exponential(1.0/upper_bound) # params are in 1/1000s units\n",
    "\n",
    "                #print 'Candidate next IAT:', t + s\n",
    "                #if s >= step: # upper bound no longer valid\n",
    "                #    t += step\n",
    "                #    step = 3600.0 # one hour in seconds\n",
    "                #    #print '\\tLarger than step, moved to:', t\n",
    "                #    continue\n",
    "                \n",
    "                cumulative_iat += s # move ahead in time (may accept or reject)\n",
    "                #step -= s\n",
    "                lambda_t_dim = self.intensity(t0 + cumulative_iat, (prev_t, prev_d))\n",
    "                lambda_t = np.sum(lambda_t_dim)\n",
    "                u = np.random.uniform(0.0, 1.0)\n",
    "                \n",
    "                if u <= lambda_t/upper_bound: # accept\n",
    "                    simulated_iats[i] = cumulative_iat\n",
    "                    if self._ndim > 1:\n",
    "                        simulated_dims[i] = np.random.choice(range(self._ndim), p=lambda_t_dim/lambda_t)\n",
    "                    break\n",
    "    \n",
    "        return simulated_iats, simulated_dims\n",
    "\n",
    "    def intensity(self, t0, prev_e, include_t0=False):\n",
    "        prev_t, prev_d = prev_e\n",
    "        if not include_t0:\n",
    "            prev_t_ = prev_t[prev_t < t0]\n",
    "            prev_d_ = prev_d[prev_t < t0]\n",
    "        else:\n",
    "            prev_t_ = prev_t[prev_t <= t0]\n",
    "            prev_d_ = prev_d[prev_t <= t0]\n",
    "        t0_minus_ti = (t0 - prev_t_) / 1000.0\n",
    "        return self._lambda + np.sum(self._beta[:, prev_d_] * np.exp(-self._alpha[:, prev_d_] * t0_minus_ti), axis=1)\n",
    "    \n",
    "    def intensity_integral(self, start_t, end_t, data_epochs, data_dims):\n",
    "        # memoize locations of timestamps for each dimension for this sample (2s for 18K samples)\n",
    "        epochs_per_dimension = []\n",
    "        for i in range(self._ndim):\n",
    "            idx_i = data_dims == i\n",
    "            epochs_i = data_epochs[idx_i]\n",
    "            epochs_per_dimension.append(epochs_i)\n",
    "        \n",
    "        integral_value = 0.0\n",
    "        \n",
    "        for m in range(self._ndim):\n",
    "            integral_value += self._lambda[m] * (end_t - start_t) / 1000.0\n",
    "            for n in range(self._ndim):\n",
    "                beta_mn = self._beta[m,n]\n",
    "                alpha_mn = self._alpha[m,n]\n",
    "                epochs_n = epochs_per_dimension[n]\n",
    "                \n",
    "                epochs_n_ = epochs_n[epochs_n < start_t]\n",
    "                term1 = np.sum(np.exp(-alpha_mn * (start_t - epochs_n_) / 1000.0) -\n",
    "                               np.exp(-alpha_mn * (end_t - epochs_n_) / 1000.0))\n",
    "                \n",
    "                epochs_n_ = epochs_n[(epochs_n >= start_t) & (epochs_n < end_t)]\n",
    "                term2 = np.sum(1 - np.exp(-alpha_mn * (end_t - epochs_n_) / 1000.0))\n",
    "                \n",
    "                integral_value += beta_mn/alpha_mn * (term1 + term2)\n",
    "        \n",
    "        return integral_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "class TimeVaryingMultivariateHawkesProcess:\n",
    "    _ndim = None # no. of dimensions\n",
    "    _lambda = None # ndim x 1 vector\n",
    "    _alpha = None # ndim x ndim matrix\n",
    "    _beta = None # ndim x ndim matrix\n",
    "    _mu = None # ndim x nfeatures matrix\n",
    "    _observation_window = None\n",
    "    _EPS = 10**-9 # numbers below this are set to 0\n",
    "\n",
    "    def __init__(self, ndim=None, lam=None, alpha=None, beta=None, mu=None, observation_window=None):\n",
    "        self._ndim = ndim\n",
    "        self._lambda = lam\n",
    "        self._alpha = alpha\n",
    "        self._beta = beta\n",
    "        self._mu = mu\n",
    "        self._observation_window = observation_window\n",
    "    \n",
    "    def fit(self, ndim, alpha, F, D, data_epochs, data_dims, observation_window, niter=1000):\n",
    "        print 'Fitting...'\n",
    "        print '\\tPrecomputing A(i), B(i)...'\n",
    "        A, B = self._precompute(ndim, alpha, data_epochs, data_dims, observation_window)\n",
    "                \n",
    "        print '\\tMinimizing negative loglikelihood...'\n",
    "        alpha_ravel = alpha.ravel()\n",
    "        nfeatures = D.shape[0]\n",
    "        x0 = np.random.uniform(low=10**-5.0, high=alpha_ravel[0] - self._EPS, size=(ndim * (ndim + 1) + ndim*nfeatures))\n",
    "        bounds = [(self._EPS, None) for i in range(ndim)] +\\\n",
    "                 [(self._EPS, None) for i in range(ndim * ndim)] +\\\n",
    "                 [(self._EPS, None) for i in range(ndim * nfeatures)]\n",
    "                 #[(self._EPS, None) for i in range(ndim * ndim)]\n",
    "                 \n",
    "        x, f, d = scipy.optimize.fmin_l_bfgs_b(self._seq_negloglikelihood,\n",
    "                                               x0=x0, args=(ndim, alpha, data_dims, F, D, observation_window, A, B),\n",
    "                                               factr=10.0**7.0, iprint=99, pgtol=10**-5, approx_grad=False,\n",
    "                                               disp=99,\n",
    "                                               bounds=bounds, callback=self._callback, maxiter=niter, maxfun=niter)\n",
    "\n",
    "        return x, f, d\n",
    "    \n",
    "    @staticmethod\n",
    "    def _precompute(ndim, alpha, data_epochs, data_dims, observation_window):\n",
    "        T0, T = observation_window\n",
    "        nsamples = len(data_epochs)\n",
    "\n",
    "        A = []\n",
    "        B = [np.zeros(shape=(ndim, ndim), dtype=np.float64) for k in range(nsamples)]\n",
    "        for k in range(nsamples):\n",
    "            # memoize locations of timestamps for each dimension for this sample (2s for 18K samples)\n",
    "            epochs_per_dimension = []\n",
    "            for i in range(ndim):\n",
    "                idx_i = data_dims[k] == i\n",
    "                epochs_i = data_epochs[k][idx_i]\n",
    "                epochs_per_dimension.append(epochs_i)\n",
    "            \n",
    "            A_k = [np.zeros(shape=(ndim, len(epochs_per_dimension[m])), dtype=np.float64)\n",
    "                   for m in range(ndim)]\n",
    "            #B_mn = np.zeros(shape=(ndim, ndim), dtype=np.float64)\n",
    "            for m in range(ndim):\n",
    "                epochs_m = epochs_per_dimension[m]\n",
    "                if len(epochs_m) == 0:\n",
    "                    continue # no transactions on dimension m\n",
    "                #A_m = np.zeros(shape=(ndim, len(epochs_m)), dtype=np.float64)\n",
    "                for n in range(ndim):\n",
    "                    epochs_n = epochs_per_dimension[n]\n",
    "                    if len(epochs_n) == 0: # no transactions on dimension n\n",
    "                        continue\n",
    "                    len_n = len(epochs_n)\n",
    "                    A_mn = A_k[m][n] # influence of dimension n on each timestamp of dimension m\n",
    "\n",
    "                    j = 0\n",
    "                    \n",
    "                    if (m != n): # mutual excitation    \n",
    "                        \n",
    "                        while j < len_n and epochs_n[j] < epochs_m[0]:\n",
    "                            j += 1\n",
    "                            \n",
    "                        #A_mn[0] = np.sum(np.exp(-alpha[m,n] * (epochs_m[0] - epochs_n[epochs_n < epochs_m[0]]) / 1000.0))\n",
    "                        A_mn[0] = np.sum(np.exp(-alpha[m,n] * (epochs_m[0] - epochs_n[:j]) / 1000.0))\n",
    "\n",
    "                        for i in range(1, len(epochs_m)):\n",
    "                            prev_j = j\n",
    "                            while j < len_n and epochs_n[j] < epochs_m[i]:\n",
    "                                j += 1\n",
    "                            A_mn[i] = np.exp(-alpha[m,n] * (epochs_m[i] - epochs_m[i-1]) / 1000.0) * A_mn[i-1]\n",
    "                            #A_mn[i] += np.sum(np.exp(-alpha[m,n] *\n",
    "                            #                         (epochs_m[i] - epochs_n[(epochs_n >= epochs_m[i-1]) &\n",
    "                            #                                                 (epochs_n < epochs_m[i])]) / 1000.0))\n",
    "                            A_mn[i] += np.sum(np.exp(-alpha[m,n] *\n",
    "                                                     (epochs_m[i] - epochs_n[prev_j:j]) / 1000.0))\n",
    "                        \n",
    "                        B[k][m,n] = np.sum(1 - np.exp(-alpha[m,n] * (T - epochs_n) / 1000.0))\n",
    "                    else: # dim_m = dim_n, self-excitation\n",
    "                        for i in range(1, len(epochs_m)):\n",
    "                            A_mn[i] = np.exp(-alpha[m,n] * (epochs_m[i] - epochs_m[i-1]) / 1000.0) * (1 + A_mn[i-1])\n",
    "                        B[k][m,n] = np.sum(1 - np.exp(-alpha[m,n] * (T - epochs_m) / 1000.0))\n",
    "                \n",
    "                #A_k.append(A_m)\n",
    "            \n",
    "            A.append(A_k)\n",
    "            #B[k] = B_mn\n",
    "\n",
    "        return A, B\n",
    "\n",
    "    def _seq_negloglikelihood(self, params, ndim, alpha, data_dims, F, D, observation_window, A, B):\n",
    "        nsamples = len(A)\n",
    "        nfeatures = D.shape[0]\n",
    "        nll = 0.0\n",
    "        grad = np.zeros(shape=((ndim + 1) * ndim + ndim*nfeatures), dtype=np.float64) # first row is grad of lambda\n",
    "        for k in range(nsamples):\n",
    "            sample_nll, sample_grad = self._negloglikelihood(params, ndim, alpha, data_dims[k],\n",
    "                                                             F[k], D, observation_window, A[k], B[k])\n",
    "            nll += sample_nll\n",
    "            grad += sample_grad\n",
    "        nll /= nsamples\n",
    "        grad /= nsamples\n",
    "        return nll, grad\n",
    "    \n",
    "    @staticmethod\n",
    "    def _negloglikelihood(params, ndim, alpha, data_dims, F, D, observation_window, A, B):\n",
    "        nfeatures = D.shape[0]\n",
    "        lam = params[:ndim] # first ndim elements\n",
    "        beta = params[ndim:ndim+ndim*ndim].reshape((ndim, ndim)) # next ndim * ndim elements\n",
    "        mu = params[ndim+ndim*ndim:ndim+ndim*ndim + ndim*nfeatures].reshape((ndim,nfeatures))\n",
    "        T0, T = observation_window\n",
    "        T_length = (T - T0)/1000.0 # scale down\n",
    "        \n",
    "        loglikelihood = 0.0\n",
    "        grad_lambda = np.zeros(shape=(ndim), dtype=np.float64)\n",
    "        grad_beta = np.zeros(shape=(ndim, ndim), dtype=np.float64)\n",
    "        grad_mu = np.zeros(shape=(ndim, nfeatures), dtype=np.float64)\n",
    "\n",
    "        for m in range(ndim):\n",
    "            mu_m = mu[m,:]\n",
    "            dimindex = data_dims == m\n",
    "            F_m = F[:,dimindex]\n",
    "            C_m = np.dot(mu_m, F_m)\n",
    "\n",
    "            Z = lam[m] + np.sum(beta[m,:][:,None] * A[m], axis=0)\n",
    "            Z += C_m\n",
    "            loglikelihood += np.sum(np.log(Z)) - np.sum(beta[m,:] * B[m,:] / alpha[m,:]) - lam[m] * T_length\n",
    "            loglikelihood -= np.sum(mu_m * D)\n",
    "            grad_lambda[m] = np.sum(1.0/Z) - T_length\n",
    "            grad_beta[m] = np.sum(A[m]/Z, axis=1) - B[m,:]/alpha[m,:]\n",
    "            grad_mu[m] = [np.sum(F_m[j]/Z) - D[j] for j in range(nfeatures)]\n",
    "        \n",
    "        gradient = np.concatenate((grad_lambda, grad_beta.ravel(), grad_mu.ravel()))\n",
    "        return -loglikelihood, -gradient\n",
    "\n",
    "    def check_gradient(self, ndim, alpha, data_epochs, data_dims, F, D, observation_window, delta=10**-8, niter=50):\n",
    "        print 'Checking gradient...'\n",
    "        print '\\tPrecomputing A, B...'\n",
    "        A, B = self._precompute(ndim, alpha, data_epochs, data_dims, observation_window)\n",
    "        \n",
    "        nfeatures = D.shape[0]\n",
    "        nsamples = len(data_epochs)\n",
    "        for i in range(niter):\n",
    "            row_idx = np.random.randint(0, nsamples)\n",
    "            a = A[row_idx]\n",
    "            b = B[row_idx]\n",
    "            f = F[row_idx]\n",
    "            d = data_dims[row_idx]\n",
    "            alpha_ravel = alpha.ravel()\n",
    "            x0 = np.random.uniform(low=10**-5.0, high=alpha_ravel[0] - self._EPS,\n",
    "                                   size=(ndim * (ndim + 1) + ndim*nfeatures))\n",
    "            nll, grad = self._negloglikelihood(x0, ndim, alpha, d, f, D, observation_window, a, b)\n",
    "            param_idx = np.random.randint(len(x0))\n",
    "            x0[param_idx] += delta\n",
    "            new_nll, new_grad = self._negloglikelihood(x0, ndim, alpha, d, f, D, observation_window, a, b)\n",
    "            num_grad = (new_nll - nll)/delta            \n",
    "            if abs((grad[param_idx] - num_grad)/num_grad)  > 0.01:\n",
    "                print 'Gradient checking failed for parameter:', param_idx\n",
    "                print 'Analytical gradient:', grad[param_idx],\n",
    "                print 'Numerical gradient:', num_grad\n",
    "        print '\\tGradient checking complete.'\n",
    "        \n",
    "    def simulate(self, t0, prev_t, prev_d, horizon_duration):\n",
    "        simulated_iats = []\n",
    "        simulated_dims = []\n",
    "        prev_t_copy = copy.copy(list(prev_t))\n",
    "        prev_d_copy = copy.copy(list(prev_d))\n",
    "        prev_t_array = np.array(prev_t_copy, dtype=np.float64)\n",
    "        prev_d_array = np.array(prev_d_copy, dtype=np.int)\n",
    "        cumulative_iat = 0.0\n",
    "        while cumulative_iat < horizon_duration:\n",
    "            upper_bound = np.sum(self.intensity(t0 + cumulative_iat, (prev_t_array, prev_d_array), include_t0=True))\n",
    "            s = 1000.0 * np.random.exponential(1.0/upper_bound) # params are in 1/1000s units            \n",
    "                \n",
    "            cumulative_iat += s # move ahead in time (may accept or reject)\n",
    "            if cumulative_iat >= horizon_duration:\n",
    "                break\n",
    "\n",
    "            lambda_t_dim = self.intensity(t0 + cumulative_iat, (prev_t_array, prev_d_array))\n",
    "            lambda_t = np.sum(lambda_t_dim)\n",
    "\n",
    "            u = np.random.uniform(0.0, 1.0)\n",
    "            if u <= lambda_t/upper_bound: # accept\n",
    "                prev_t_copy.append(t0 + cumulative_iat)\n",
    "                prev_t_array = np.array(prev_t_copy, dtype=np.float64)\n",
    "                simulated_iats.append(cumulative_iat)\n",
    "                \n",
    "                cumprob = np.cumsum(lambda_t_dim/lambda_t)\n",
    "                u2 = np.random.uniform(0.0, 1.0)\n",
    "                dim = np.searchsorted(cumprob, u2)\n",
    "                \n",
    "                prev_d_copy.append(dim)\n",
    "                prev_d_array = np.array(prev_d_copy, dtype=np.int)\n",
    "                simulated_dims.append(dim)\n",
    "\n",
    "        return simulated_iats, simulated_dims\n",
    "    \n",
    "    def predict_dist(self, t0, prev_t, prev_d, niter=1000):\n",
    "        step = 3600.0 - (t0 % 3600)\n",
    "        simulated_iats = np.zeros(shape=(niter), dtype=np.float64)\n",
    "        simulated_dims = np.zeros(shape=(niter), dtype=np.int)\n",
    "        for i in range(niter):\n",
    "            cumulative_iat = 0.0\n",
    "            while True:\n",
    "                #print 'trial start'\n",
    "                upper_bound = np.sum(self.intensity(t0 + cumulative_iat, (prev_t, prev_d), include_t0=True))\n",
    "                s = 1000.0 * np.random.exponential(1.0/upper_bound) # params are in 1/1000s units\n",
    "\n",
    "                #print 'Candidate next IAT:', t + s\n",
    "                if s >= step: # upper bound no longer valid\n",
    "                    cumulative_iat += step\n",
    "                    step = 3600.0 # one hour in seconds\n",
    "                    #print '\\tLarger than step, moved to:', t\n",
    "                    continue\n",
    "                \n",
    "                cumulative_iat += s # move ahead in time (may accept or reject)\n",
    "                step -= s\n",
    "                lambda_t_dim = self.intensity(t0 + cumulative_iat, (prev_t, prev_d))\n",
    "                lambda_t = np.sum(lambda_t_dim)\n",
    "                u = np.random.uniform(0.0, 1.0)\n",
    "                \n",
    "                if u <= lambda_t/upper_bound: # accept\n",
    "                    simulated_iats[i] = cumulative_iat\n",
    "                    if self._ndim > 1:\n",
    "                        #print 'simulating dim'\n",
    "                        simulated_dims[i] = np.random.choice(range(self._ndim), p=lambda_t_dim/lambda_t)\n",
    "                        #print 'done simulating dim'\n",
    "                    break\n",
    "    \n",
    "        return simulated_iats, simulated_dims\n",
    "\n",
    "    def intensity(self, t0, prev_e, include_t0=False):\n",
    "        prev_t, prev_d = prev_e\n",
    "        if not include_t0:\n",
    "            prev_t_ = prev_t[prev_t < t0]\n",
    "            prev_d_ = prev_d[prev_t < t0]\n",
    "        else:\n",
    "            prev_t_ = prev_t[prev_t <= t0]\n",
    "            prev_d_ = prev_d[prev_t <= t0]\n",
    "        \n",
    "        t0_minus_ti = (t0 - prev_t_) / 1000.0\n",
    "        intensity = self._lambda + np.sum(self._beta[:, prev_d_] * np.exp(-self._alpha[:, prev_d_] * t0_minus_ti), axis=1)\n",
    "        \n",
    "        t0_ = time.gmtime(t0)\n",
    "        hour_of_day, day_of_week, day_of_month = map(int, time.strftime(\"%H %w %d\", t0_).split(' '))\n",
    "        day_of_week = (day_of_week - 1) % 7\n",
    "        \n",
    "        f = np.array([int(hour_of_day == h) for h in range(24)] +\n",
    "                     [int(day_of_week >= 0 and day_of_week <= 3)] +\n",
    "                     [int(day_of_week == 4)] +\n",
    "                     [int(day_of_week >= 5 and day_of_week <= 6)] +\n",
    "                     [int(day_of_month == 1)], dtype=np.int)\n",
    "\n",
    "        #print intensity.shape, self._mu.shape, f.shape, (self._mu * f).shape, np.sum(self._mu * f, axis=0).shape\n",
    "        intensity += np.sum(self._mu * f, axis=1)\n",
    "        \n",
    "        return intensity\n",
    "    \n",
    "    def intensity_integral(self, start_t, end_t, data_epochs, data_dims):\n",
    "        # memoize locations of timestamps for each dimension for this sample (2s for 18K samples)\n",
    "        epochs_per_dimension = []\n",
    "        for i in range(self._ndim):\n",
    "            idx_i = data_dims == i\n",
    "            epochs_i = data_epochs[idx_i]\n",
    "            epochs_per_dimension.append(epochs_i)\n",
    "        \n",
    "        integral_value = 0.0\n",
    "        \n",
    "        for m in range(self._ndim):\n",
    "            integral_value += self._lambda[m] * (end_t - start_t) / 1000.0\n",
    "            for n in range(self._ndim):\n",
    "                beta_mn = self._beta[m,n]\n",
    "                alpha_mn = self._alpha[m,n]\n",
    "                epochs_n = epochs_per_dimension[n]\n",
    "                \n",
    "                epochs_n_ = epochs_n[epochs_n < start_t]\n",
    "                term1 = np.sum(np.exp(-alpha_mn * (start_t - epochs_n_) / 1000.0) -\n",
    "                               np.exp(-alpha_mn * (end_t - epochs_n_) / 1000.0))\n",
    "                \n",
    "                epochs_n_ = epochs_n[(epochs_n >= start_t) & (epochs_n < end_t)]\n",
    "                term2 = np.sum(1 - np.exp(-alpha_mn * (end_t - epochs_n_) / 1000.0))\n",
    "                \n",
    "                integral_value += beta_mn/alpha_mn * (term1 + term2)\n",
    "        \n",
    "        # mu term integral\n",
    "        d = np.zeros(28, dtype=np.float64)\n",
    "        start_timestamp = pd.to_datetime(start_t, unit='s')\n",
    "        end_timestamp = pd.to_datetime(end_t, unit='s')\n",
    "\n",
    "        #print 't_i-1', start_timestamp, 't_i', end_timestamp\n",
    "        \n",
    "        start_ts_roundup = start_timestamp.floor(freq='H')\n",
    "        start_seconds = (start_timestamp - start_ts_roundup).total_seconds()\n",
    "        start_hour = start_timestamp.hour\n",
    "        end_ts_roundwn = end_timestamp.ceil(freq='H')\n",
    "        end_seconds = (end_ts_roundwn - end_timestamp).total_seconds()\n",
    "        end_hour = end_timestamp.hour\n",
    "\n",
    "        #print 't_i-1', start_timestamp, 't_i', end_timestamp\n",
    "        \n",
    "        #hour_periods = pd.date_range(start_timestamp, end_timestamp, freq='s', closed='left')\n",
    "        hour_periods2 = pd.date_range(start_ts_roundup, end_ts_roundwn, freq='H', closed='left')\n",
    "        #print hour_periods\n",
    "        #hour_counts = Counter(hour_periods.hour)\n",
    "        hour_counts2 = Counter(hour_periods2.hour)\n",
    "        hour_counts2[start_hour] -= start_seconds/3600.0\n",
    "        hour_counts2[end_hour] -= end_seconds/3600.0\n",
    "        #print hour_counts2\n",
    "        #print start_timestamp, start_ts_roundup, end_timestamp, end_ts_roundwn\n",
    "        for h in range(24):\n",
    "            #d[h] = hour_counts[h] / 1000.0\n",
    "            d[h] = (hour_counts2[h] * 3600.0) / 1000.0\n",
    "            #print h, hour_counts[h], hour_counts2[h] * 3600.0\n",
    "            #assert abs(d[h] - (hour_counts2[h] * 3600.0) / 1000.0) < 10**-9\n",
    "        \n",
    "        start_ts_roundup = start_timestamp.floor(freq='D')\n",
    "        start_seconds = (start_timestamp - start_ts_roundup).total_seconds()\n",
    "        start_day = start_timestamp.weekday()\n",
    "        end_ts_roundwn = end_timestamp.ceil(freq='D')\n",
    "        end_seconds = (end_ts_roundwn - end_timestamp).total_seconds()\n",
    "        end_day = end_timestamp.weekday()\n",
    "        \n",
    "        #print start_timestamp, start_ts_roundup, start_seconds/3600.0, end_timestamp, end_ts_roundwn, end_seconds/3600.0\n",
    "        #day_counts = Counter(hour_periods.weekday)\n",
    "        day_periods = pd.date_range(start_ts_roundup, end_ts_roundwn, freq='D', closed='left')\n",
    "        day_counts2 = Counter(day_periods.weekday)\n",
    "        day_counts2[start_day] -= start_seconds/86400.0\n",
    "        day_counts2[end_day] -= end_seconds/86400.0\n",
    "        #for d_i in range(7):\n",
    "        #    assert abs(day_counts[d_i]/86400.0 - day_counts2[d_i]) < 10**-9\n",
    "        d[24] = (day_counts2[0] + day_counts2[1] + day_counts2[2] + day_counts2[3]) * 86400.0 / 1000.0\n",
    "        d[25] = day_counts2[4] * 86400.0/ 1000.0\n",
    "        d[26] = (day_counts2[5] + day_counts2[6]) * 86400.0 / 1000.0\n",
    "        \n",
    "        num_month_starts = np.sum(day_periods.is_month_start)\n",
    "        #print num_month_starts\n",
    "        if start_timestamp.is_month_start:\n",
    "            num_month_starts -= start_seconds/86400.0\n",
    "        if end_timestamp.is_month_start:\n",
    "            num_month_starts -= end_seconds/86400.0\n",
    "        #print start_timestamp, end_timestamp, np.sum(hour_periods.is_month_start), num_month_starts*86400.0\n",
    "        #assert abs(np.sum(hour_periods.is_month_start) - num_month_starts*86400.0) < 10**-9\n",
    "        #d[27] = np.sum(hour_periods.is_month_start) / 1000.0\n",
    "        d[27] = num_month_starts * 86400.0 / 1000.0\n",
    "        d = np.array(d)\n",
    "        #print d\n",
    "        \n",
    "        mu_term = np.sum(self._mu * d)\n",
    "        integral_value += mu_term\n",
    "        \n",
    "        return integral_value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hawkes Process: Simulate Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating 100 data samples from a 2-D Hawkes process...\n",
      "Done.\n"
     ]
    }
   ],
   "source": [
    "simulated_epochs = []\n",
    "simulated_dims = []\n",
    "observation_window = [0.0, 600.0]\n",
    "t0 = observation_window[0]\n",
    "ndim = 2\n",
    "true_lambda = 1000.0 * np.array([0.1, 0.1])\n",
    "true_alpha = 1000.0 * np.array([[1.0, 1.0], [1.0, 1.0]])\n",
    "true_beta = 1000.0 * np.array([[0.5, 0.01], [0.9, 0.4]])\n",
    "true_params = np.concatenate((true_lambda, true_beta.ravel()))\n",
    "print \"Generating 100 data samples from a \" + str(ndim) + \"-D Hawkes process...\"\n",
    "for k in range(100):\n",
    "    iats, dims = SimpleMultivariateHawkesProcess(ndim=ndim, lam=true_lambda, alpha=true_alpha, beta=true_beta,\n",
    "                                                 observation_window=observation_window).simulate(t0,\n",
    "                                                                                                 np.array([]),\n",
    "                                                                                                 np.array([]),\n",
    "                                                                                                 observation_window[-1])\n",
    "    epochs = np.array(iats, dtype=np.float64)\n",
    "    dims = np.array(dims, dtype=np.int)\n",
    "    simulated_epochs.append(epochs)\n",
    "    simulated_dims.append(dims)\n",
    "print \"Done.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hawkes Process: Fit Model\n",
    "\n",
    "Note that alpha is not learned from the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting...\n",
      "\tPrecomputing A(i), B(i)...\n",
      "\tMinimizing negative loglikelihood...\n",
      "Negative loglikelihood: -2325.514811743731\n",
      "Estimated params: [0.10101926 0.10131509 0.52093352 0.00443929 0.88401316 0.40952672]\n",
      "True params: [0.1  0.1  0.5  0.01 0.9  0.4 ]\n",
      "Termination info:\n",
      "\tCONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL\n",
      "\tIterations: 19 Funcalls: 22\n"
     ]
    }
   ],
   "source": [
    "est_params, est_nll, d = SimpleMultivariateHawkesProcess().fit(ndim, true_alpha, simulated_epochs,\n",
    "                                                               simulated_dims, observation_window, niter=1000)\n",
    "print \"Negative loglikelihood:\", est_nll\n",
    "print \"Estimated params:\", est_params/1000.0\n",
    "print \"True params:\", true_params/1000.0\n",
    "print \"Termination info:\"\n",
    "print \"\\t\", d[\"task\"]\n",
    "print \"\\t\", \"Iterations:\", d[\"nit\"], \"Funcalls:\", d[\"funcalls\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time-Varying Hawkes Process: Simulate Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating 100 data samples from a 2-D Time-Varying Hawkes process...\n",
      "Done.\n"
     ]
    }
   ],
   "source": [
    "simulated_epochs = []\n",
    "simulated_dims = []\n",
    "observation_window = [0.0, 60.0]\n",
    "t0 = observation_window[0]\n",
    "ndim = 2\n",
    "true_lambda = 1000.0 * np.array([0.1, 0.1])\n",
    "true_alpha = 1000.0 * np.array([[1.0, 1.0], [1.0, 1.0]])\n",
    "true_beta = 1000.0 * np.array([[0.5, 0.01], [0.9, 0.4]])\n",
    "true_mu = np.array([[0.2]*6 + [2.0]*12 + [0.2]*6 + [0.1, 2.0, 1.5, 0.0],\n",
    "                    [0.2]*6 + [0.4]*12 + [2.0]*6 + [2.0, 0.1, 1.0, 0.0]])\n",
    "true_params = np.concatenate((true_lambda, true_beta.ravel(), true_mu.ravel()))\n",
    "print \"Generating 100 data samples from a \" + str(ndim) + \"-D Time-Varying Hawkes process...\"\n",
    "for k in range(100):\n",
    "    iats, dims = TimeVaryingMultivariateHawkesProcess(ndim=ndim, lam=true_lambda, alpha=true_alpha, beta=true_beta,\n",
    "                                                      mu=true_mu,\n",
    "                                                      observation_window=observation_window).simulate(t0,\n",
    "                                                                                                      np.array([]),\n",
    "                                                                                                      np.array([]),\n",
    "                                                                                                      observation_window[-1])\n",
    "    epochs = np.array(iats, dtype=np.float64)\n",
    "    dims = np.array(dims, dtype=np.int)\n",
    "    simulated_epochs.append(epochs)\n",
    "    simulated_dims.append(dims)\n",
    "print \"Done.\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time-Varying Hawkes Process: Fit Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Precompute D, F"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "est_params, est_nll, d = TimeVaryingMultivariateHawkesProcess().fit(ndim, true_alpha, simulated_epochs,\n",
    "                                                                    simulated_dims, observation_window, niter=1000)\n",
    "print \"Negative loglikelihood:\", est_nll\n",
    "print \"Estimated params:\", est_params/1000.0\n",
    "print \"True params:\", true_params/1000.0\n",
    "print \"Termination info:\"\n",
    "print \"\\t\", d[\"task\"]\n",
    "print \"\\t\", \"Iterations:\", d[\"nit\"], \"Funcalls:\", d[\"funcalls\"]"
   ]
  }
 ],
 "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}