Online.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coordinate Descent"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Lower Bound, Take 4\n",
    "\n",
    "Ensure feasibility of \"Lower Bound, Take 2\" by adjusting alpha as necessary."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Lower Bound, Take 3\n",
    "\n",
    "Ensure feasibility by allowing a stochastic mixture with the MLE.  Doesn't work (not DCP)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Assume $r_{\\min} = 0$ for simplicity.  Idea for online solving $$\n",
    "\\begin{aligned}\n",
    "&\\!\\min_{Q \\succeq 0} &\\qquad& \\sum_{(w,r)} w r Q_{w, r}, \\\\\n",
    "&\\text{subject to} &  & \\sum_{(w,r)} w Q_{w,r} = 1, \\\\\n",
    "&                  &  & \\sum_{(w,r)} Q_{w,r} = 1, \\\\\n",
    "&                  &  & \\sum_n \\log(Q_{w_n, r_n}) \\geq \\phi\n",
    "\\end{aligned}\n",
    "$$ where $\\phi = -\\frac{1}{2} \\chi^{2,\\alpha}_{(1)} + \\sum_n \\log(Q^{\\text{(mle)}}_{w_n, r_n})$.  Because the support of $Q$ is at most the empirical support plus $(w_{\\min}, 0)$ and $(w_{\\max}, 0)$ we will maintain two variables $q_{\\min}$ and $q_{\\max}$ corresponding to $w_{\\min}$ and $w_{\\max}$ respectively.  Otherwise we need one primal variable for each data point.  However we will use two primal variables corresponding to $r = 0$ and $r = 1$ respectively.  We will split each datapoint into two points whose fractional counts are proportional to $r$ and $(1 - r)$ respectively.\n",
    "\n",
    "At time $t$ we receive $(w_t, r_t)$ and we want to determine $q_t$.  We are allowed to choose the new $q_{\\min}$ and $q_{\\max}$ arbitrarily.  For points $q_{<t}$ we are allowed to scale them by $\\psi_0$ and $\\psi_1$, corresponding to points with $r = 0$ and $r = 1$ respectively.  We are also allowed to stochastically mix in the maximum likelihood solution according to $\\psi_{\\text{mle}}$.  We assume $(q_{0,<t}, q_{1,<t}, q_{<t,\\min}, q_{<t,\\max})$ is feasible before receiving $(w_t, r_t)$.  Then \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\!\\min_{q_{0,t}, q_{1,t}, q_{\\min}, q_{\\max}, \\psi_0, \\psi_1 \\succeq 0, \\psi_{\\text{mle}} \\in [0, 1]} &\\qquad& q_{1,t} w_t + \\psi_1 v_{<t} + \\psi_{\\text{mle}} v_{\\text{mle}} \\\\\n",
    "&\\text{subject to} &  & w_t q_{0,t} + w_t q_{1,t} + w_{\\min} q_{\\min} + w_{\\max} q_{\\max} + \\psi_0 w_{0,<t} q_{0,<t} + \\psi_1 w_{1,<t} q_{1,<t} = 1, \\\\\n",
    "&                  &  & q_{0,t} + q_{1,t} + q_{\\min} + q_{\\max} + \\psi_0 q_{0,<t} + \\psi_1 q_{1,<t} = 1, \\\\\n",
    "&                  &  & (1 - r_t) \\log(q_{0,t}) + r_t \\log(q_{1,t}) + (t - r_{<t}) \\log(\\psi_0) + r_{<t} \\log(\\psi_1) \\geq \\phi - \\mathcal{L}_{<t} = -\\frac{1}{2} \\chi^{2,\\alpha}_{(1)} + \\mathcal{L}^{\\text{(mle)}}_t - \\mathcal{L}_{<t}\n",
    "\\end{aligned}\n",
    "$$ where $v_{<t}$ is the previous lower bound, $\\mathcal{L}^{\\text{(mle)}}_t$ is the mle likelihood of the observed data including point $t$, and $\\mathcal{L}_{<t}$ is the previously obtained likelihood of the observed data for the lower bound."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Lower Bound, Take 2\n",
    "\n",
    "Better, but still has feasibility issues as the stream progresses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Assume $r_{\\min} = 0$ for simplicity.  Idea for online solving $$\n",
    "\\begin{aligned}\n",
    "&\\!\\min_{Q \\succeq 0} &\\qquad& \\sum_{(w,r)} w r Q_{w, r}, \\\\\n",
    "&\\text{subject to} &  & \\sum_{(w,r)} w Q_{w,r} = 1, \\\\\n",
    "&                  &  & \\sum_{(w,r)} Q_{w,r} = 1, \\\\\n",
    "&                  &  & \\sum_n \\log(Q_{w_n, r_n}) \\geq \\phi\n",
    "\\end{aligned}\n",
    "$$ where $\\phi = -\\frac{1}{2} \\chi^{2,\\alpha}_{(1)} + \\sum_n \\log(Q^{\\text{(mle)}}_{w_n, r_n})$.  Because the support of $Q$ is at most the empirical support plus $(w_{\\min}, 0)$ and $(w_{\\max}, 0)$ we will maintain two variables $q_{\\min}$ and $q_{\\max}$ corresponding to $w_{\\min}$ and $w_{\\max}$ respectively.  Otherwise we need one primal variable for each data point.  \n",
    "\n",
    "At time $t$ we receive $(w_t, r_t)$ and we want to determine $q_t$.  We are allowed to choose the new $q_{\\min}$ and $q_{\\max}$ arbitrarily.  For points $q_{<t}$ we are allowed to scale them by $\\psi_0$ and $\\psi_1$, corresponding to points with $r = 0$ and $r = 1$ respectively.  We assume $(q_{0,<t}, q_{1,<t}, q_{<t,\\min}, q_{<t,\\max})$ is feasible before receiving $(w_t, r_t)$.  Then \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\!\\min_{q_t, q_{\\min}, q_{\\max},\\psi_0, \\psi_1 \\succeq 0} &\\qquad& q_t w_t r_t + \\psi_1 v_{<t} \\\\\n",
    "&\\text{subject to} &  & w_t q_t + w_{\\min} q_{\\min} + w_{\\max} q_{\\max} + \\psi_0 w_{0,<t} q_{0,<t} + \\psi_1 w_{1,<t} q_{1,<t} = 1, \\\\\n",
    "&                  &  & q_t + q_{\\min} + q_{\\max} + \\psi_0 q_{0,<t} + \\psi_1 q_{1,<t} = 1, \\\\\n",
    "&                  &  & \\log(q_t) + (t - r_{<t}) \\log(\\psi_0) + r_{<t} \\log(\\psi_1) \\geq \\phi - \\mathcal{L}_{<t} = -\\frac{1}{2} \\chi^{2,\\alpha}_{(1)} + \\mathcal{L}^{\\text{(mle)}}_t - \\mathcal{L}_{<t}\n",
    "\\end{aligned}\n",
    "$$ where $v_{<t}$ is the previous lower bound, $\\mathcal{L}^{\\text{(mle)}}_t$ is the mle likelihood of the observed data including point $t$, and $\\mathcal{L}_{<t}$ is the previously obtained likelihood of the observed data for the lower bound.\n",
    "\n",
    "After computing $q_t$ we increment $q_{0,<t}$ and $q_{1,<t}$ proportional to $1 - r_t$ and $r_t$ respectively; and $w_{0,<t}$ and $w_{1,<t}$ by $w_t (1 - r_t)$ and $w_t r_t$ respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Lower Bound\n",
    "\n",
    "Starts out good, but then runs into infeasibility."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Assume $r_{\\min} = 0$ for simplicity.  Idea for online solving $$\n",
    "\\begin{aligned}\n",
    "&\\!\\min_{Q \\succeq 0} &\\qquad& \\sum_{(w,r)} w r Q_{w, r}, \\\\\n",
    "&\\text{subject to} &  & \\sum_{(w,r)} w Q_{w,r} = 1, \\\\\n",
    "&                  &  & \\sum_{(w,r)} Q_{w,r} = 1, \\\\\n",
    "&                  &  & \\sum_n \\log(Q_{w_n, r_n}) \\geq \\phi\n",
    "\\end{aligned}\n",
    "$$ where $\\phi = -\\frac{1}{2} \\chi^{2,\\alpha}_{(1)} + \\sum_n \\log(Q^{\\text{(mle)}}_{w_n, r_n})$.  Because the support of $Q$ is at most the empirical support plus $(w_{\\min}, 0)$ and $(w_{\\max}, 0)$ we will maintain two variables $q_{\\min}$ and $q_{\\max}$ corresponding to $w_{\\min}$ and $w_{\\max}$ respectively.  Otherwise we need one primal variable for each data point.\n",
    "\n",
    "At time $t$ we receive $(w_t, r_t)$ and we want to determine $q_t$.  We are allowed to choose the new $q_{\\min}$ and $q_{\\max}$ arbitrarily.  For points $q_{<t}$ we only allowed to scale them by $\\psi$.  We assume $(q_{<t}, q_{<t,\\min}, q_{<t,\\max})$ is feasible before receiving $(w_t, r_t)$.  Then \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\!\\min_{q_t, q_{\\min}, q_{\\max}, \\psi \\succeq 0} &\\qquad& q_t w_t r_t + \\psi v_{<t} \\\\\n",
    "&\\text{subject to} &  & w_t q_t + w_{\\min} q_{\\min} + w_{\\max} q_{\\max} + \\psi w_{<t} q_{<t} = 1, \\\\\n",
    "&                  &  & q_t + q_{\\min} + q_{\\max} + \\psi q_{<t} = 1, \\\\\n",
    "&                  &  & \\log(q_t) + t \\log(\\psi) \\geq \\phi - \\mathcal{L}_{<t} = -\\frac{1}{2} \\chi^{2,\\alpha}_{(1)} + \\mathcal{L}^{\\text{(mle)}}_t - \\mathcal{L}_{<t}\n",
    "\\end{aligned}\n",
    "$$ where $v_{<t}$ is the previous lower bound, $\\mathcal{L}^{\\text{(mle)}}_t$ is the mle likelihood of the observed data including point $t$, and $\\mathcal{L}_{<t}$ is the previously obtained likelihood of the observed data for the lower bound. Substituting $q_{<t} + q_{<t,\\min} + q_{<t,\\max} = 1$ and $w_{<t} q_{<t} + w_{\\min} q_{<t,\\min} + w_{\\max} q_{<t,\\max} = 1$ yields \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\!\\min_{q_t, q_{\\min}, q_{\\max} \\succeq 0, \\psi \\in [0, 1]} &\\qquad& q_t w_t r_t + \\psi v_{<t} \\\\\n",
    "&\\text{subject to} &  & w_t q_t + w_{\\min} q_{\\min} + w_{\\max} q_{\\max} + \\psi (1 - w_{\\min} q_{<t,\\min} - w_{\\max} q_{<t,\\max}) = 1, \\\\\n",
    "&                  &  & q_t + q_{\\min} + q_{\\max} + \\psi (1 - q_{<t,\\max} - q_{<t,\\min}) = 1 \\\\\n",
    "&                  &  & -\\log(q_t) - t \\log(\\psi) - \\frac{1}{2} \\chi^{2,\\alpha}_{(1)} + \\mathcal{L}^{\\text{(mle)}}_t  - \\mathcal{L}_{<t} \\leq 0\n",
    "\\end{aligned}\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### MLE"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "Idea for online solving $$\n",
    "\\begin{aligned}\n",
    "&\\!\\max_{Q \\succeq 0} &\\qquad& \\sum_n \\log(Q_{w_n, r_n}), \\\\\n",
    "&\\text{subject to} &  & \\sum_{(w,r)} w Q_{w,r} = 1, \\\\\n",
    "&                  &  & \\sum_{(w,r)} Q_{w,r} = 1. \n",
    "\\end{aligned}\n",
    "$$\n",
    "Because the support of $Q$ is at most the empirical support plus $w_{\\min}$ and $w_{\\max}$ we will maintain two variables $q_{\\min}$ and $q_{\\max}$ corresponding to $w_{\\min}$ and $w_{\\max}$ respectively.  Otherwise we need one primal variable for each data point.\n",
    "\n",
    "At time $t$ we receive $(w_t, r_t)$ and we want to determine $q_t$.  We are allowed to choose the new $q_{\\min}$ and $q_{\\max}$ arbitrarily.  For points $q_{<t}$ we only allowed to scale them by $\\psi$.  We assume $(q_{<t}, q_{<t,\\min}, q_{<t,\\max})$ is feasible before receiving $(w_t, r_t)$.  Then \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\!\\max_{q_t, q_{\\min}, q_{\\max}, \\psi > 0} &\\qquad& t \\log(\\psi) + \\log(q_t), \\\\\n",
    "&\\text{subject to} &  & w_t q_t + w_{\\min} q_{\\min} + w_{\\max} q_{\\max} + \\psi w_{<t} q_{<t} = 1, \\\\\n",
    "&                  &  & q_t + q_{\\min} + q_{\\max} + \\psi q_{<t} = 1\n",
    "\\end{aligned}\n",
    "$$ \n",
    "\n",
    "Substituting $q_{<t} + q_{<t,\\min} + q_{<t,\\max} = 1$ and $w_{<t} q_{<t} + w_{\\min} q_{<t,\\min} + w_{\\max} q_{<t,\\max} = 1$ yields \n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\!\\max_{q_t, q_{\\min}, q_{\\max}, \\psi > 0} &\\qquad& t \\log(\\psi) + \\log(q_t), \\\\\n",
    "&\\text{subject to} &  & w_t q_t + w_{\\min} q_{\\min} + w_{\\max} q_{\\max} + \\psi (1 - w_{\\min} q_{<t,\\min} - w_{\\max} q_{<t,\\max}) = 1, \\\\\n",
    "&                  &  & q_t + q_{\\min} + q_{\\max} + \\psi (1 - q_{<t,\\max} - q_{<t,\\min}) = 1 \n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "At the beginning of time we can initialize with $$\n",
    "\\begin{aligned}\n",
    "q_{0,\\min} &= \\frac{1 - w_{\\min}}{w_{\\max} - w_{\\min}} \\\\\n",
    "q_{0,\\max} &= \\frac{w_{\\max} - 1}{w_{\\max} - w_{\\min}}\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "code_folding": [
     0,
     1,
     26,
     79,
     177,
     247,
     315,
     432,
     517,
     660
    ],
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1, 2, 1, 1, '0.0733', '0.073', {'vlb': '0.1465', 'vmle': '1.0', 'alpha': 0.05}]\n",
      "[1, 2, 0, 2, '0.481', '0.962', {'vlb': '0.038075', 'vmle': '0.5', 'alpha': 0.05}]\n",
      "[1, 0, 0, 5, '0.247', '1.236', {'vlb': '0.021139', 'vmle': '0.3', 'alpha': 0.05}]\n",
      "[1, 2, 1, 17, '0.0223', '0.379', {'vlb': '0.22307', 'vmle': '0.5593', 'alpha': 0.05}]\n",
      "[1, 2, 1, 65, '0.0102', '0.662', {'vlb': '0.42774', 'vmle': '0.63908', 'alpha': 0.05}]\n",
      "[1, 2, 0, 257, '0.00427', '1.097', {'vlb': '0.6275', 'vmle': '0.74588', 'alpha': 0.05}]\n",
      "[1, 2, 1, 1025, '0.000895', '0.917', {'vlb': '0.66016', 'vmle': '0.72023', 'alpha': 0.05}]\n",
      "[1, 1000, 0, 3805, '0.000129', '0.489', {'vlb': '0.67922', 'vmle': '0.71246', 'alpha': 0.05}]\n",
      "[1, 0, 0, 4097, '0.00025', '1.025', {'vlb': '0.67433', 'vmle': '0.70771', 'alpha': 0.05}]\n",
      "[1, 2, 1, 16385, '5.97e-05', '0.978', {'vlb': '0.71508', 'vmle': '0.73192', 'alpha': 0.05}]\n",
      "[1, 1000, 1, 23861, '5.34e-06', '0.127', {'vlb': '0.73212', 'vmle': '0.81397', 'alpha': 0.05}]\n",
      "[1, 1000, 1, 30958, '7e-06', '0.217', {'vlb': '0.7432', 'vmle': '0.81855', 'alpha': 0.05}]\n",
      "[1, 1000, 0, 32227, '8.16e-05', '2.630', {'vlb': '0.73097', 'vmle': '0.81724', 'alpha': 0.05}]\n",
      "[1, 0, 0, 65537, '1.54e-05', '1.008', {'vlb': '0.72676', 'vmle': '0.7734', 'alpha': 0.05}]\n"
     ]
    }
   ],
   "source": [
    "class OnlineCoordinateDescentMLE:\n",
    "    def __init__(self, wmin, wmax):\n",
    "        from cvxopt import matrix\n",
    "\n",
    "        assert wmax > 1\n",
    "        assert wmin >= 0\n",
    "        assert wmin < wmax\n",
    "        \n",
    "        self.wmin = wmin\n",
    "        self.wmax = wmax\n",
    "        self.qmin = (wmax - 1) / (wmax - wmin)\n",
    "        self.qmax = (1 - wmin) / (wmax - wmin)\n",
    "        self.obj = 0\n",
    "        self.vmin = 0\n",
    "        self.lastphi = 0\n",
    "        \n",
    "        self.G = matrix([ [ -1, 0, 0, 0 ],\n",
    "                          [ 0, -1, 0, 0 ],\n",
    "                          [ 0, 0, -1, 0 ],\n",
    "                          [ 0, 0, 0, -1 ],\n",
    "                        ],\n",
    "                        tc='d').T\n",
    "        self.h = matrix([ 0, 0, 0, 0 ], tc='d')\n",
    "        self.b = matrix([ 1 / wmax, 1 ], tc='d')        \n",
    "        self.t = 0\n",
    "        \n",
    "    def update(self, c, w, r):\n",
    "        from cvxopt import matrix, solvers\n",
    "        \n",
    "        assert c > 0\n",
    "    \n",
    "        safet = max(self.t, 1)\n",
    "        x0 = matrix([ c / (c + safet), \n",
    "                      self.qmin * safet / (c + safet), \n",
    "                      self.qmax * safet / (c + safet), \n",
    "                      safet / (c + safet) ], \n",
    "                    tc='d')\n",
    "\n",
    "        def F(x=None, z=None):\n",
    "            import math\n",
    "            \n",
    "            if x is None: return 0, x0\n",
    "            \n",
    "            if x[0] <= 0 or x[3] <= 0:\n",
    "                return None\n",
    "            \n",
    "            f = -c * math.log(x[0]) / safet - self.t * math.log(x[3]) / safet\n",
    "            jf = matrix([ -c / (safet * x[0]), 0, 0, -self.t / (safet * x[3]) ], tc='d').T\n",
    "            if z is None: return f, jf\n",
    "            hf = z[0] * matrix([ [ (c / safet) * 1/x[0]**2, 0, 0, 0 ], \n",
    "                                 [ 0, 0, 0, 0 ],\n",
    "                                 [ 0, 0, 0, 0 ],\n",
    "                                 [ 0, 0, 0, (self.t / safet) * 1/x[3]**2 ] \n",
    "                               ], tc='d')\n",
    "            return f, jf, hf\n",
    "          \n",
    "        A = matrix([ \n",
    "                     [ float(w) / self.wmax, \n",
    "                       self.wmin / self.wmax, \n",
    "                      1, \n",
    "                      (1 / self.wmax - (self.wmin / self.wmax) * self.qmin - self.qmax) ],\n",
    "                     [ 1, 1, 1, (1 - self.qmin - self.qmax) ]\n",
    "                   ],\n",
    "                   tc='d')\n",
    "                \n",
    "        soln = solvers.cp(F=F, G=self.G, h=self.h, A=A.T, b=self.b, options={'show_progress': False})\n",
    "        from pprint import pformat\n",
    "        assert soln['status'] == 'optimal', pformat([ soln, self.t ])\n",
    "        \n",
    "        self.obj -= safet * soln['primal objective']\n",
    "        self.lastq = soln['x'][0]\n",
    "        self.qmin = soln['x'][1]\n",
    "        self.qmax = soln['x'][2]\n",
    "        self.lastphi = soln['x'][3]\n",
    "        self.vmin = soln['x'][0] * w * r + soln['x'][3] * self.vmin\n",
    "        self.t += c\n",
    "        \n",
    "        return self.lastq\n",
    "\n",
    "class OnlineCoordinateDescentLB: \n",
    "    class Flass:\n",
    "        def __init__(self):\n",
    "            pass\n",
    "        \n",
    "    def __init__(self, wmin, wmax, alpha): \n",
    "        assert wmax > 1\n",
    "        assert wmin >= 0\n",
    "        assert wmin < wmax\n",
    "        \n",
    "        self.wmin = wmin\n",
    "        self.wmax = wmax\n",
    "        self.qmin = (wmax - 1) / (wmax - wmin)\n",
    "        self.qmax = (1 - wmin) / (wmax - wmin)\n",
    "        \n",
    "        self.vlb = 0\n",
    "        self.wq0 = 0\n",
    "        self.wq1 = 0\n",
    "        self.q0t = 0\n",
    "        self.q1t = 0\n",
    "        self.t = 0\n",
    "        self.rt = 0\n",
    "        self.llb = 0\n",
    "        \n",
    "        self.alpha = alpha\n",
    "        self.mle = OnlineCoordinateDescentMLE(wmin=wmin, wmax=wmax)\n",
    "        \n",
    "        from scipy.stats import chi2\n",
    "        import cvxpy as cp\n",
    "        \n",
    "        self.vars = OnlineCoordinateDescentLB.Flass()\n",
    "        self.vars.qt = cp.Variable(nonneg=True)\n",
    "        self.vars.qmin = cp.Variable(nonneg=True)\n",
    "        self.vars.qmax = cp.Variable(nonneg=True)\n",
    "        self.vars.psi0 = cp.Variable(nonneg=True)\n",
    "        self.vars.psi1 = cp.Variable(nonneg=True)\n",
    "        \n",
    "        self.params = OnlineCoordinateDescentLB.Flass()\n",
    "        self.params.w = cp.Parameter(nonneg=True)\n",
    "        self.params.wr = cp.Parameter(nonneg=True)\n",
    "        self.params.vlb = cp.Parameter(nonneg=True)\n",
    "        self.params.wq0 = cp.Parameter(nonneg=True)\n",
    "        self.params.wq1 = cp.Parameter(nonneg=True)\n",
    "        self.params.q0t = cp.Parameter(nonneg=True)\n",
    "        self.params.q1t = cp.Parameter(nonneg=True)\n",
    "        self.params.c = cp.Parameter(nonneg=True)\n",
    "        self.params.tminusrt = cp.Parameter(nonneg=True)\n",
    "        self.params.rt = cp.Parameter(nonneg=True)\n",
    "        self.params.constraintrhs = cp.Parameter()\n",
    "        \n",
    "        self.prob = cp.Problem(cp.Minimize(self.params.wr * self.vars.qt + self.params.vlb * self.vars.psi1), [\n",
    "              self.params.w * self.vars.qt\n",
    "            + (self.wmin / self.wmax) * self.vars.qmin \n",
    "            + self.vars.qmax \n",
    "            + self.params.wq0 * self.vars.psi0 \n",
    "            + self.params.wq1 * self.vars.psi1\n",
    "            == 1 / self.wmax,  \n",
    "              self.vars.qt \n",
    "            + self.vars.qmin \n",
    "            + self.vars.qmax \n",
    "            + self.params.q0t * self.vars.psi0\n",
    "            + self.params.q1t * self.vars.psi1 == 1,\n",
    "              self.params.c * cp.log(self.vars.qt) \n",
    "            + self.params.tminusrt * cp.log(self.vars.psi0)\n",
    "            + self.params.rt * cp.log(self.vars.psi1) \n",
    "            >= self.params.constraintrhs\n",
    "        ])\n",
    "        \n",
    "    def innersolve(self, c, w, r, alpha):\n",
    "        from scipy.stats import chi2\n",
    "\n",
    "        safet = max(self.t, 1)\n",
    "        halfchisq = 0.5 * chi2.isf(q=alpha, df=1)\n",
    "        \n",
    "        self.params.w.value = w / self.wmax\n",
    "        self.params.wr.value = w * r\n",
    "        self.params.vlb.value = self.vlb\n",
    "        self.params.wq0.value = self.wq0 / self.wmax\n",
    "        self.params.wq1.value = self.wq1 / self.wmax\n",
    "        self.params.q0t.value = self.q0t\n",
    "        self.params.q1t.value = self.q1t\n",
    "        self.params.c.value = c / safet\n",
    "        self.params.tminusrt.value = (self.t - self.rt) / safet\n",
    "        self.params.rt.value = self.rt / safet\n",
    "        self.params.constraintrhs.value = (-halfchisq + self.mle.obj - self.llb) / safet\n",
    "        \n",
    "        self.prob.solve(verbose=False)\n",
    "        \n",
    "        return (self.prob.value, \n",
    "                self.prob.status, \n",
    "                (1 - r) * self.vars.qt.value if self.vars.qt.value is not None else None, \n",
    "                r * self.vars.qt.value if self.vars.qt.value is not None else None, \n",
    "                self.vars.qmin.value, \n",
    "                self.vars.qmax.value, \n",
    "                self.vars.psi0.value, \n",
    "                self.vars.psi1.value\n",
    "               )\n",
    "        \n",
    "    def updatev3(self, c, w, r):\n",
    "        import cvxpy as cp\n",
    "        import math\n",
    "        from scipy.special import xlogy\n",
    "        \n",
    "        qmle = self.mle.update(c, w, r)\n",
    "        \n",
    "        q0t = cp.Variable(nonneg=True)\n",
    "        q1t = cp.Variable(nonneg=True)\n",
    "        qmin = cp.Variable(nonneg=True)\n",
    "        qmax = cp.Variable(nonneg=True)\n",
    "        psi0 = cp.Variable(nonneg=True)\n",
    "        psi1 = cp.Variable(nonneg=True)\n",
    "        psimle = cp.Variable(nonneg=True)\n",
    "        \n",
    "        safet = max(self.t, 1)\n",
    "\n",
    "#        prob = cp.Problem(cp.Minimize((float(w) * q1t + self.vlb * psi1) * (1 - psimle) + self.mle.vhat * psimle), [\n",
    "        prob = cp.Problem(cp.Minimize(float(w) * q1t + self.vlb * psi1 + self.mle.vhat * psimle), [\n",
    "              float(w / self.wmax) * q0t\n",
    "            + float(w / self.wmax) * q1t\n",
    "            + (self.wmin / self.wmax) * qmin \n",
    "            + qmax \n",
    "            + psi0 * (self.wq0 / self.wmax)\n",
    "            + psi1 * (self.wq1 / self.wmax) \n",
    "            == 1 / self.wmax,  \n",
    "            q0t + q1t + qmin + qmax + psi0 * self.q0t + psi1 * self.q1t == 1,\n",
    "              float(c * (1 - r) / safet) * cp.log(q0t) \n",
    "            + float(c * r / safet) * cp.log(q1t) \n",
    "            + float((self.t - self.rt) / safet) * cp.log(psi0)\n",
    "            + float(self.rt / safet) * cp.log(psi1) \n",
    "#             >= ((-self.halfchisq + self.mle.obj) / safet) * cp.inv_pos(1 - psimle) - (self.llb / safet),\n",
    "            >= ((-self.halfchisq + self.mle.obj) / safet) * (1 + psimle) - (self.llb / safet),\n",
    "            psimle <= 1,\n",
    "            psimle >= 0.9\n",
    "        ])\n",
    "        prob.solve(verbose=False)\n",
    "        assert prob.status[:7] == 'optimal', prob.solve(verbose=True)\n",
    "        \n",
    "        self.vlb = (w * q1t.value + self.vlb * psi1.value) * (1 - psimle.value) + psimle.value * self.mle.vhat\n",
    "        self.qmin = qmin.value\n",
    "        self.qmax = qmax.value\n",
    "        self.llb += (  xlogy(c * (1 - r), q0t.value)\n",
    "                     + xlogy(c * r, q1t.value)\n",
    "                     + xlogy(self.t - self.rt, psi0.value)\n",
    "                     + xlogy(self.rt, psi1.value)\n",
    "                    )\n",
    "        self.llb = (1 - psimle.value) * self.llb + psimle.value * self.mle.obj\n",
    "     \n",
    "        self.t += c\n",
    "        self.rt += c * r\n",
    "        self.q0t = q0t.value + psi0.value * self.q0t\n",
    "        self.wq0 = w * q0t.value + psi0.value * self.wq0\n",
    "        self.q1t = q1t.value + psi1.value * self.q1t\n",
    "        self.wq1 = w * q1t.value + psi1.value * self.wq1\n",
    "        \n",
    "        myq = (1 - r) * q0t.value + r * q1t.value\n",
    "        myq = (1 - psimle.value) * myq + psimle.value * qmle\n",
    "        \n",
    "        return myq, { 'self.q0t': self.q0t, \n",
    "                      'self.q1t': self.q1t,\n",
    "                      'qmin': qmin.value,\n",
    "                      'qmax': qmax.value,\n",
    "                      'vlb': self.vlb,\n",
    "                      'llb': self.llb,\n",
    "                      'vmle': self.mle.vhat,\n",
    "                      'self.rt': self.rt,\n",
    "                      'psimle': psimle.value,\n",
    "                    }\n",
    "\n",
    "    def __initturg__(self, wmin, wmax, alpha): \n",
    "        assert wmax > 1\n",
    "        assert wmin >= 0\n",
    "        assert wmin < wmax\n",
    "        \n",
    "        self.wmin = wmin\n",
    "        self.wmax = wmax\n",
    "        self.qmin = (wmax - 1) / (wmax - wmin)\n",
    "        self.qmax = (1 - wmin) / (wmax - wmin)\n",
    "        \n",
    "        self.vlb = 0\n",
    "        self.wq0 = 0\n",
    "        self.wq1 = 0\n",
    "        self.q0t = 0\n",
    "        self.q1t = 0\n",
    "        self.t = 0\n",
    "        self.rt = 0\n",
    "        self.llb = 0\n",
    "        \n",
    "        self.alpha = alpha\n",
    "        self.mle = OnlineCoordinateDescentMLE(wmin=wmin, wmax=wmax)\n",
    "        \n",
    "        from scipy.stats import chi2\n",
    "        import cvxpy as cp\n",
    "        \n",
    "        self.vars = OnlineCoordinateDescentLB.Flass()\n",
    "        self.vars.q0t = cp.Variable(nonneg=True)\n",
    "        self.vars.q1t = cp.Variable(nonneg=True)\n",
    "        self.vars.qmin = cp.Variable(nonneg=True)\n",
    "        self.vars.qmax = cp.Variable(nonneg=True)\n",
    "        self.vars.psi0 = cp.Variable(nonneg=True)\n",
    "        self.vars.psi1 = cp.Variable(nonneg=True)\n",
    "        \n",
    "        self.params = OnlineCoordinateDescentLB.Flass()\n",
    "        self.params.w = cp.Parameter(nonneg=True)\n",
    "        self.params.wcost = cp.Parameter(nonneg=True)\n",
    "        self.params.vlb = cp.Parameter(nonneg=True)\n",
    "        self.params.wq0 = cp.Parameter(nonneg=True)\n",
    "        self.params.wq1 = cp.Parameter(nonneg=True)\n",
    "        self.params.q0t = cp.Parameter(nonneg=True)\n",
    "        self.params.q1t = cp.Parameter(nonneg=True)\n",
    "        self.params.coneminusr = cp.Parameter(nonneg=True)\n",
    "        self.params.cr = cp.Parameter(nonneg=True)\n",
    "        self.params.tminusrt = cp.Parameter(nonneg=True)\n",
    "        self.params.rt = cp.Parameter(nonneg=True)\n",
    "        self.params.constraintrhs = cp.Parameter()\n",
    "        \n",
    "        self.prob = cp.Problem(cp.Minimize(self.params.wcost * self.vars.q1t + self.params.vlb * self.vars.psi1), [\n",
    "              self.params.w * self.vars.q0t\n",
    "            + self.params.w * self.vars.q1t\n",
    "            + (self.wmin / self.wmax) * self.vars.qmin \n",
    "            + self.vars.qmax \n",
    "            + self.params.wq0 * self.vars.psi0 \n",
    "            + self.params.wq1 * self.vars.psi1\n",
    "            == 1 / self.wmax,  \n",
    "              self.vars.q0t \n",
    "            + self.vars.q1t \n",
    "            + self.vars.qmin \n",
    "            + self.vars.qmax \n",
    "            + self.params.q0t * self.vars.psi0\n",
    "            + self.params.q1t * self.vars.psi1 == 1,\n",
    "              self.params.coneminusr * cp.log(self.vars.q0t) \n",
    "            + self.params.cr * cp.log(self.vars.q1t) \n",
    "            + self.params.tminusrt * cp.log(self.vars.psi0)\n",
    "            + self.params.rt * cp.log(self.vars.psi1) \n",
    "            >= self.params.constraintrhs\n",
    "        ])\n",
    "\n",
    "    def innersolveflass(self, c, w, r, alpha):   \n",
    "        # doesn't work, not sure why (?)\n",
    "\n",
    "        from cvxopt import matrix, spdiag, solvers\n",
    "        import numpy as np\n",
    "        from scipy.stats import chi2\n",
    "        \n",
    "        assert 0 < c\n",
    "        assert 0 <= r\n",
    "        assert r <= 1\n",
    "\n",
    "        safet = max(self.t, 1)\n",
    "        halfchisq = 0.5 * chi2.isf(q=alpha, df=1)\n",
    "         \n",
    "        print([\n",
    "            w / self.wmax,\n",
    "            w,\n",
    "            self.vlb,\n",
    "            self.wq0 / self.wmax,\n",
    "            self.wq1 / self.wmax,\n",
    "            self.q0t,\n",
    "            self.q1t,\n",
    "            c * (1 - r) / safet,\n",
    "            c * r / safet,\n",
    "            (self.t - self.rt) / safet,\n",
    "            self.rt / safet,\n",
    "            (-halfchisq + self.mle.obj - self.llb) / safet\n",
    "        ])\n",
    "        \n",
    "        \n",
    "        G = matrix(-np.eye(6), tc='d')\n",
    "        h = matrix(0, size=(6,1), tc='d')\n",
    "        \n",
    "        A = matrix([ [ w / self.wmax,\n",
    "                       w / self.wmax,\n",
    "                       self.wmin / self.wmax, \n",
    "                       1, \n",
    "                       self.wq0 / self.wmax,\n",
    "                       self.wq1 / self.wmax \n",
    "                     ],\n",
    "                     [ 1, 1, 1, 1, self.q0t, self.q1t ]\n",
    "                   ],\n",
    "                   tc='d')\n",
    "        b = matrix([ 1 / self.wmax, 1 ], tc='d')\n",
    "        \n",
    "        cost = matrix(0, size=(6,1), tc='d')\n",
    "        cost[1] = float(w) \n",
    "        cost[5] = self.vlb \n",
    "        \n",
    "        x0 = matrix([ c * (1 - r) / (c + safet),\n",
    "                      c * r / (c + safet),\n",
    "                      self.qmin * safet / (c + safet), \n",
    "                      self.qmax * safet / (c + safet), \n",
    "                      safet / (c + safet),\n",
    "                      safet / (c + safet)\n",
    "                    ], tc='d')\n",
    "\n",
    "        def F(x=None, z=None):\n",
    "            from scipy.special import xlogy\n",
    "\n",
    "            if x is None: return 1, x0\n",
    "            \n",
    "            if any(z < 0 for z in x):\n",
    "                 return None\n",
    "                \n",
    "            f = (  xlogy(c * (1 - r), x[0])\n",
    "                 + xlogy(c * r, x[1])\n",
    "                 + xlogy(self.t - self.rt, x[3])\n",
    "                 + xlogy(self.rt, x[4])\n",
    "                 + halfchisq\n",
    "                 - self.mle.obj\n",
    "                 + self.llb)\n",
    "            f *= -1 / safet\n",
    "            \n",
    "            jf = matrix(0, size=(1, 6), tc='d')\n",
    "            jf[0] = c * (1 - r) / x[0] if c * (1 - r) > 0 else 0 \n",
    "            jf[1] = c * r / x[1] if c * r > 0 else 0 \n",
    "            jf[4] = (self.t - self.rt) / x[4] if self.t > self.rt else 0\n",
    "            jf[5] = self.rt / x[5] if self.rt > 0 else 0\n",
    "            jf *= -1 / safet\n",
    "            \n",
    "            if z is None: return f, jf\n",
    "            \n",
    "            hf = spdiag([\n",
    "                -c * (1 - r) / x[0]**2 if c * (1 - r) > 0 else 0,\n",
    "                 -c * r / x[1]**2 if c * r > 0 else 0,               \n",
    "                0,\n",
    "                0,\n",
    "                -(self.t - self.rt) / x[4]**2 if self.t > self.rt else 0,\n",
    "                -self.rt / x[5]**2 if self.rt > 0 else 0 \n",
    "            ]) \n",
    "\n",
    "            hf *= -z[0] / safet\n",
    "\n",
    "            return f, jf, hf\n",
    "\n",
    "        soln = solvers.cpl(c=cost, F=F, G=G, h=h, A=A.T, b=b)\n",
    "\n",
    "        from pprint import pformat\n",
    "        import numpy\n",
    "        assert soln['status'][:7] == 'optimal', pformat({ 'soln': soln, \n",
    "                                                  'solnx': [ z for z in soln['x'] ], \n",
    "                                                  'datum': (c, w, r),\n",
    "                                                  'F(x=x0)': F(x=x0),\n",
    "                                                  'A': numpy.matrix(A),\n",
    "                                                  'b': [ z for z in b ],\n",
    "                                                  'A.x0 - b': [ z for z in A.T*x0 - b ],\n",
    "                                                  'G.x0 - h': [ z for z in G*x0 - h ],\n",
    "                                                  'F(x=soln)': F(x=soln['x']),\n",
    "                                                  'A.x - b': [ z for z in A.T*soln['x'] - b ],\n",
    "                                                  'G.x - h': [ z for z in G*soln['x'] - h ],\n",
    "                                                })\n",
    "        return (soln['primal objective'], soln['status'], \n",
    "                (1 - r) * soln['x'][0],\n",
    "                r * soln['x'][0]\n",
    "               ) + tuple(soln['x'][1:])\n",
    "\n",
    "    def innersolveturg(self, c, w, r, alpha):\n",
    "        from scipy.stats import chi2\n",
    "\n",
    "        safet = max(self.t, 1)\n",
    "        halfchisq = 0.5 * chi2.isf(q=alpha, df=1)\n",
    "        \n",
    "        self.params.w.value = w / self.wmax\n",
    "        self.params.wcost.value = w\n",
    "        self.params.vlb.value = self.vlb\n",
    "        self.params.wq0.value = self.wq0 / self.wmax\n",
    "        self.params.wq1.value = self.wq1 / self.wmax\n",
    "        self.params.q0t.value = self.q0t\n",
    "        self.params.q1t.value = self.q1t\n",
    "        self.params.coneminusr.value = c * (1 - r) / safet\n",
    "        self.params.cr.value = c * r / safet\n",
    "        self.params.tminusrt.value = (self.t - self.rt) / safet\n",
    "        self.params.rt.value = self.rt / safet\n",
    "        self.params.constraintrhs.value = (-halfchisq + self.mle.obj - self.llb) / safet\n",
    "        \n",
    "        self.prob.solve(verbose=False)\n",
    "        \n",
    "        return (self.prob.value, \n",
    "                self.prob.status, \n",
    "                self.vars.q0t.value, \n",
    "                self.vars.q1t.value, \n",
    "                self.vars.qmin.value, \n",
    "                self.vars.qmax.value, \n",
    "                self.vars.psi0.value, \n",
    "                self.vars.psi1.value\n",
    "               )\n",
    "                \n",
    "    def update(self, c, w, r):\n",
    "        import math\n",
    "        from scipy.special import xlogy\n",
    "        \n",
    "        self.mle.update(c, w, r)\n",
    "        \n",
    "        alpha = self.alpha\n",
    "        \n",
    "        (pvalue, pstatus, q0t, q1t, qmin, qmax, psi0, psi1) = self.innersolve(c, w, r, alpha)\n",
    "        \n",
    "        if pstatus[:7] != 'optimal':\n",
    "            alphalb = 0\n",
    "            alphaub = alpha\n",
    "            \n",
    "            while alphaub - alphalb >= 1e-3:\n",
    "                alphatest = 0.5 * (alphalb + alphaub)\n",
    "                (pvalue, pstatus, q0t, q1t, qmin, qmax, psi0, psi1) = self.innersolve(c, w, r, alphatest)\n",
    "            \n",
    "                if pstatus[:7] == 'optimal':\n",
    "                    alphalb = alphatest\n",
    "                else:\n",
    "                    alphaub = alphatest\n",
    "                    \n",
    "            alpha = alphalb\n",
    "            (pvalue, pstatus, q0t, q1t, qmin, qmax, psi0, psi1) = self.innersolve(c, w, r, alpha)\n",
    "            \n",
    "        assert pstatus[:7] == 'optimal', { 'alpha': alpha, 'pstatus': pstatus }\n",
    "        \n",
    "        self.vlb = pvalue\n",
    "        self.qmin = qmin\n",
    "        self.qmax = qmax\n",
    "        self.llb += (  xlogy(c * (1 - r), q0t)\n",
    "                     + xlogy(c * r, q1t)\n",
    "                     + xlogy(self.t - self.rt, psi0)\n",
    "                     + xlogy(self.rt, psi1)\n",
    "                    )\n",
    "     \n",
    "        self.t += c\n",
    "        self.rt += c * r\n",
    "        self.q0t = q0t + psi0 * self.q0t\n",
    "        self.wq0 = w * q0t + psi0 * self.wq0\n",
    "        self.q1t = q1t + psi1 * self.q1t\n",
    "        self.wq1 = w * q1t + psi1 * self.wq1\n",
    "        \n",
    "        return (1 - r) * q0t + r * q1t, { # 'self.q0t': self.q0t, \n",
    "#                                           'self.q1t': self.q1t,\n",
    "#                                           'qmin': qmin.value,\n",
    "#                                           'qmax': qmax.value,\n",
    "                                            'vlb': self.vlb,\n",
    "                                            'vmle': self.mle.vmin,\n",
    "#                                           'self.rt': self.rt,\n",
    "                                            'alpha': alpha,\n",
    "                                        }\n",
    "\n",
    "    def __initv1__(self, wmin, wmax, alpha):\n",
    "        from scipy.stats import chi2\n",
    "        from cvxopt import matrix\n",
    "\n",
    "        assert wmax > 1\n",
    "        assert wmin >= 0\n",
    "        assert wmin < wmax\n",
    "        \n",
    "        self.wmin = wmin\n",
    "        self.wmax = wmax\n",
    "        self.qmin = (wmax - 1) / (wmax - wmin)\n",
    "        self.qmax = (1 - wmin) / (wmax - wmin)\n",
    "        \n",
    "        self.G = matrix([ [ -1, 0, 0, 0 ],\n",
    "                          [ 0, -1, 0, 0 ],\n",
    "                          [ 0, 0, -1, 0 ],\n",
    "                          [ 0, 0, 0, -1 ],\n",
    "                        ],\n",
    "                        tc='d').T\n",
    "        self.h = matrix([ 0, 0, 0, 0 ], tc='d')\n",
    "        self.b = matrix([ 1 / wmax, 1 ], tc='d')        \n",
    "        self.t = 0\n",
    "\n",
    "        self.halfchisq = 0.5 * chi2.isf(q=alpha, df=1)\n",
    "        self.llb = 0\n",
    "        self.vlb = 0\n",
    "        \n",
    "        self.mle = OnlineCoordinateDescentMLE(wmin=wmin, wmax=wmax)\n",
    "    \n",
    "    def updatev1(self, c, w, r):\n",
    "        import cvxpy as cp\n",
    "        import math\n",
    "        \n",
    "        self.mle.update(c, w, r)\n",
    "                \n",
    "        qt = cp.Variable(nonneg=True)\n",
    "        qmin = cp.Variable(nonneg=True)\n",
    "        qmax = cp.Variable(nonneg=True)\n",
    "        psi = cp.Variable(nonneg=True)\n",
    "        \n",
    "        safet = max(self.t, 1)\n",
    "\n",
    "        prob = cp.Problem(cp.Minimize(float(w * r) * qt + self.vlb * psi), [\n",
    "              float(w / self.wmax) * qt \n",
    "            + (self.wmin / self.wmax) * qmin \n",
    "            + qmax \n",
    "            + (1 / self.wmax - (self.wmin / self.wmax) * self.qmin - self.qmax) * psi == 1 / self.wmax,\n",
    "            qt + qmin + qmax + (1 - self.qmin - self.qmax) * psi == 1,\n",
    "            float(c / safet) * cp.log(qt) + (self.t / safet) * cp.log(psi) >= (-self.halfchisq + self.mle.obj - self.llb) / safet\n",
    "        ])\n",
    "        prob.solve(verbose=False)\n",
    "        \n",
    "        if prob.status[:7] != 'optimal':\n",
    "            # just maximize likelihood to recover (?)\n",
    "            pass\n",
    "        \n",
    "        from pprint import pformat\n",
    "        assert prob.status[:7] == 'optimal', pformat({ 'datum': [ c, w, r ], \n",
    "                                                   'x': [ z.value for z in [ qt, qmin, qmax, psi ] ],\n",
    "                                                   'prob': prob.status,\n",
    "                                                   't': self.t, \n",
    "                                                   'lmle': self.mle.obj, \n",
    "                                                   'llb': self.llb,# + c * math.log(qt.value) + self.t * math.log(psi.value),\n",
    "                                                   'halfchisq': self.halfchisq,\n",
    "                                                 })\n",
    "        \n",
    "        self.vlb = prob.value\n",
    "        self.qmin = qmin.value\n",
    "        self.qmax = qmax.value\n",
    "        self.lastphi = psi.value\n",
    "        self.llb += c * math.log(qt.value) + self.t * math.log(psi.value)\n",
    "        self.t += c\n",
    "\n",
    "        return { 'soln': [ z.value for z in [ qt, qmin, qmax, psi ] ], 'lmle': self.mle.obj, \n",
    "                 'llb': self.llb, 'halfchisq': self.halfchisq, 'vlb': self.vlb, 'vhat': self.mle.vhat }\n",
    "       \n",
    "        \n",
    "        if False:\n",
    "            from cvxopt import matrix, solvers\n",
    "            import math\n",
    "\n",
    "            assert c > 0 \n",
    "\n",
    "            self.mle.update(c, w, r)\n",
    "            lmle = self.mle.obj\n",
    "\n",
    "            safet = max(self.t, 1)\n",
    "            x0 = matrix([ c / (c + safet), \n",
    "                          self.qmin * safet / (c + safet), \n",
    "                          self.qmax * safet / (c + safet), \n",
    "                          safet / (c + safet) ], tc='d')\n",
    "\n",
    "            def F(x=None, z=None):\n",
    "                if x is None: return 1, x0\n",
    "\n",
    "                if x[0] <= 0 or x[3] <= 0:\n",
    "                     return None\n",
    "\n",
    "                f = -c * math.log(x[0]) - self.t * math.log(x[3]) - self.halfchisq + lmle - self.llb\n",
    "                f /= safet\n",
    "                jf = matrix([ -c / (safet * x[0]), 0, 0, -self.t / (safet * x[3]) ], tc='d').T\n",
    "                if z is None: return f, jf\n",
    "                hf = z[0] * matrix([ [ (c / safet) * 1/x[0]**2, 0, 0, 0 ], \n",
    "                                     [ 0, 0, 0, 0 ],\n",
    "                                     [ 0, 0, 0, 0 ],\n",
    "                                     [ 0, 0, 0, (self.t / safet) * 1/x[3]**2 ] \n",
    "                                   ], tc='d')\n",
    "                return f, jf, hf\n",
    "\n",
    "            cost = matrix([ float(w * r), 0, 0, self.vlb ], tc='d')\n",
    "            A = matrix([ \n",
    "                         [ float(w) / self.wmax, \n",
    "                           self.wmin / self.wmax, \n",
    "                          1, \n",
    "                          (1 / self.wmax - (self.wmin / self.wmax) * self.qmin - self.qmax) ],\n",
    "                         [ 1, 1, 1, (1 - self.qmin - self.qmax) ]\n",
    "                       ],\n",
    "                       tc='d')\n",
    "\n",
    "            soln = solvers.cpl(c=cost, F=F, G=self.G, h=self.h, A=A.T, b=self.b, options={'show_progress': True, \n",
    "                                                                                          'maxiters': 100 })\n",
    "            from pprint import pformat\n",
    "            assert soln['status'] == 'optimal', pformat({ 'datum': [ c, w, r ], \n",
    "                                                          'x': [ z for z in soln['x'] ], \n",
    "                                                          'soln': soln,\n",
    "                                                          't': self.t, \n",
    "                                                          'lmle': lmle, \n",
    "                                                          'llb': self.llb, \n",
    "                                                          'F(x0)': F(x=x0),\n",
    "                                                          'F(x)': F(soln['x'])\n",
    "                                                        })\n",
    "            print(pformat((F(x=x0), F(x=soln['x']))))\n",
    "\n",
    "            self.vlb = soln['primal objective']\n",
    "            self.qmin = soln['x'][1]\n",
    "            self.qmax = soln['x'][2]\n",
    "            self.lastphi = soln['x'][3]\n",
    "            self.llb += c * math.log(soln['x'][0]) + self.t * math.log(soln['x'][3])\n",
    "            self.t += c\n",
    "\n",
    "            return { 'soln': [ z for z in soln['x'] ], 'lmle': lmle, \n",
    "                     'llb': self.llb, 'halfchisq': self.halfchisq, 'vlb': self.vlb, 'vhat': self.mle.vhat }\n",
    "\n",
    "class Test:\n",
    "    def flass():\n",
    "        import environments.ControlledRangeVariance\n",
    "        import MLE.MLE\n",
    "        \n",
    "        def batchtoonline(samples, seed=45):\n",
    "            import numpy as np\n",
    "\n",
    "            state = np.random.RandomState(seed)\n",
    "\n",
    "            n = sum(c for c, w, r in samples)\n",
    "            while n > 0:\n",
    "                p = np.array([ c for c, w, r in samples ], dtype='float64') / n\n",
    "                what = state.choice(len(samples), p=p)\n",
    "                c = min(samples[what][0], 1)\n",
    "                yield (c, samples[what][1], samples[what][2])\n",
    "                samples[what] = (samples[what][0] - c, samples[what][1], samples[what][2]) \n",
    "                n -= c\n",
    "                \n",
    "        env = environments.ControlledRangeVariance.ControlledRangeVariance(seed=45, wsupport=[0,2,1000], expwsq=100)\n",
    "        _, samples = env.sample(1 + (1 << 16))\n",
    "        ocd = OnlineCoordinateDescentLB(wmin=0, wmax=1000, alpha=0.05)\n",
    "\n",
    "        from pprint import pformat\n",
    "        \n",
    "        t = 0\n",
    "        for n, (c, w, r) in enumerate(batchtoonline(samples)):\n",
    "            qt, qex = ocd.update(c, w, r)\n",
    "            t += c\n",
    "            if (n & (n - 1) == 0 and n & 0xAAAAAAAA == 0) or w == 1000:\n",
    "                vlb = '{:.5}'.format(qex['vlb'])\n",
    "                vmle = '{:.5}'.format(qex['vmle'])\n",
    "                print([ c, w, r, t, '{:.3g}'.format(qt), '{:.3f}'.format(t * qt), { 'vlb': vlb, 'vmle': vmle, 'alpha': qex['alpha'] } ], \n",
    "                      flush=True)\n",
    "\n",
    "Test.flass()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# Histogram Based"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 353,
   "metadata": {
    "hidden": true,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(35819, 0, 0), (5658, 2, 0), (24056, 2, 1), (2, 1000, 0), (2, 1000, 1)]\n",
      "[1, 2, 1, 1, 'n/a', 'n/a', {'vlb': '0.000', 'vmle': '1.000', 'alpha': 0.05}]\n",
      "[1, 2, 0, 2, 'n/a', 'n/a', {'vlb': '-0.000', 'vmle': '0.500', 'alpha': 0.05}]\n",
      "[1, 0, 0, 5, 'n/a', 'n/a', {'vlb': '0.003', 'vmle': '0.400', 'alpha': 0.05}]\n",
      "[1, 2, 1, 17, 'n/a', 'n/a', {'vlb': '0.208', 'vmle': '0.589', 'alpha': 0.05}]\n",
      "[1, 2, 1, 65, 'n/a', 'n/a', {'vlb': '0.425', 'vmle': '0.647', 'alpha': 0.05}]\n",
      "[1, 2, 0, 257, 'n/a', 'n/a', {'vlb': '0.625', 'vmle': '0.748', 'alpha': 0.05}]\n",
      "[1, 2, 1, 1025, 'n/a', 'n/a', {'vlb': '0.659', 'vmle': '0.721', 'alpha': 0.05}]\n",
      "[1, 1000, 0, 3805, 'n/a', 'n/a', {'vlb': '0.678', 'vmle': '0.710', 'alpha': 0.05}]\n",
      "[1, 0, 0, 4097, 'n/a', 'n/a', {'vlb': '0.675', 'vmle': '0.706', 'alpha': 0.05}]\n",
      "[1, 2, 1, 16385, 'n/a', 'n/a', {'vlb': '0.716', 'vmle': '0.731', 'alpha': 0.05}]\n",
      "[1, 1000, 1, 23861, 'n/a', 'n/a', {'vlb': '0.733', 'vmle': '0.782', 'alpha': 0.05}]\n",
      "[1, 1000, 1, 30958, 'n/a', 'n/a', {'vlb': '0.750', 'vmle': '0.799', 'alpha': 0.05}]\n",
      "[1, 1000, 0, 32227, 'n/a', 'n/a', {'vlb': '0.746', 'vmle': '0.784', 'alpha': 0.05}]\n",
      "[1, 0, 0, 65537, 'n/a', 'n/a', {'vlb': '0.743', 'vmle': '0.781', 'alpha': 0.05}]\n"
     ]
    }
   ],
   "source": [
    "from importlib import reload\n",
    "import environments.ControlledRangeVariance\n",
    "import MLE.MLE \n",
    "\n",
    "reload(environments.ControlledRangeVariance)\n",
    "reload(MLE.MLE)\n",
    "\n",
    "def batchtoonline(samples, seed=45):\n",
    "    import numpy as np\n",
    "    \n",
    "    state = np.random.RandomState(seed)\n",
    "    \n",
    "    n = sum(c for c, w, r in samples)\n",
    "    while n > 0:\n",
    "        p = np.array([ c for c, w, r in samples ], dtype='float64') / n\n",
    "        what = state.choice(len(samples), p=p)\n",
    "        c = min(samples[what][0], 1)\n",
    "        yield (c, samples[what][1], samples[what][2])\n",
    "        samples[what] = (samples[what][0] - c, samples[what][1], samples[what][2]) \n",
    "        n -= c\n",
    "    \n",
    "env = environments.ControlledRangeVariance.ControlledRangeVariance(seed=45, wsupport=[0,2,1000], expwsq=100)\n",
    "happrox = MLE.MLE.Online.HistApprox(wmin=0, wmax=1000, numbuckets=10)\n",
    "onlineci = MLE.MLE.Online.CI(wmin=0, wmax=1000, rmin=0, rmax=1, alpha=0.05)\n",
    "onlinemle = MLE.MLE.Online.MLE(wmin=0, wmax=1000, rmin=0, rmax=1)\n",
    "_, samples = env.sample(1 + (1 << 16))\n",
    "from pprint import pformat\n",
    "print(pformat(samples), flush=True)\n",
    "\n",
    "t = 0\n",
    "for n, (c, w, r) in enumerate(batchtoonline(samples)):\n",
    "    happrox.update(c, w, r)\n",
    "    onlineci.update(happrox.iterator)\n",
    "    onlinemle.update(happrox.iterator)\n",
    "    t += c\n",
    "    \n",
    "    if (n & (n - 1) == 0 and n & 0xAAAAAAAA == 0) or w == 1000:\n",
    "        vmle = MLE.MLE.estimate(happrox.iterator, wmin=0, wmax=1000)[1]['vmin']\n",
    "        vlb = MLE.MLE.asymptoticconfidenceinterval(happrox.iterator, wmin=0, wmax=1000)[0][0]\n",
    "        print([ c, w, r, t, 'n/a', 'n/a', { 'vlb': '{:.3f}'.format(vlb), 'vmle': '{:.3f}'.format(vmle), 'alpha': 0.05 } ], \n",
    "                              flush=True)\n",
    "        \n",
    "#         from pprint import pformat\n",
    "#         print(pformat(\n",
    "#                 {\n",
    "#                     'n': n,\n",
    "#                     'onlineci': onlineci.getqfunc(),\n",
    "#                     'onlinemle': onlinemle.getqfunc(),\n",
    "#                     'batchmle': MLE.MLE.estimate(happrox.iterator, wmin=0, wmax=1000), \n",
    "#                     'batchci': MLE.MLE.asymptoticconfidenceinterval(happrox.iterator, wmin=0, wmax=1000)\n",
    "#                 }), \n",
    "#               flush=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Forest Covertype"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Full Information Online Learning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "code_folding": [
     0,
     1
    ],
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "**** lr = 1 ****\n",
      "n       \temp loss\tsince last\n",
      "32      \t0.562   \t0.562     \n",
      "64      \t0.594   \t0.625     \n",
      "128     \t0.594   \t0.594     \n",
      "256     \t0.551   \t0.508     \n",
      "512     \t0.496   \t0.441     \n",
      "1024    \t0.470   \t0.443     \n",
      "2048    \t0.431   \t0.393     \n",
      "4096    \t0.395   \t0.359     \n",
      "8192    \t0.380   \t0.365     \n",
      "16384   \t0.376   \t0.371     \n",
      "32768   \t0.355   \t0.335     \n",
      "65536   \t0.354   \t0.353     \n",
      "116160  \t0.352   \t0.350     \n",
      "test accuracy: [0.66119363 0.66225697 0.66277654]\n"
     ]
    }
   ],
   "source": [
    "class OnlineDRO:\n",
    "    class EasyAcc:\n",
    "        def __init__(self):\n",
    "            self.n = 0\n",
    "            self.sum = 0\n",
    "            \n",
    "        def __iadd__(self, other):\n",
    "            self.n += 1\n",
    "            self.sum += other\n",
    "            return self\n",
    "            \n",
    "        def mean(self):\n",
    "            return self.sum / max(self.n, 1)\n",
    " \n",
    "    def flass():\n",
    "        from sklearn.datasets import fetch_covtype\n",
    "        from sklearn.decomposition import PCA\n",
    "        from sklearn.linear_model import SGDClassifier\n",
    "        from sklearn.metrics import accuracy_score\n",
    "        from math import ceil\n",
    "        import numpy as np\n",
    "  \n",
    "        cov = fetch_covtype()\n",
    "        cov.data = PCA(whiten=True).fit_transform(cov.data)\n",
    "        classes = np.unique(cov.target - 1)\n",
    "        ndata = len(cov.target)\n",
    "        order = np.random.RandomState(seed=42).permutation(ndata)\n",
    "        ntrain = ceil(0.2 * ndata)\n",
    "        Object = lambda **kwargs: type(\"Object\", (), kwargs)()\n",
    "        train = Object(data = cov.data[order[:ntrain]], target = cov.target[order[:ntrain]] - 1)\n",
    "        test = Object(data = cov.data[order[ntrain:]], target = cov.target[order[ntrain:]] - 1)\n",
    "        \n",
    "        for lr in (1, ):\n",
    "            print(\"**** lr = {} ****\".format(lr))\n",
    "            print('{:8.8s}\\t{:8.8s}\\t{:10.10s}'.format('n', 'emp loss', 'since last'))\n",
    "            \n",
    "            classweights = { k: lr for k, _ in Counter(train.target).items() }\n",
    "            cls = SGDClassifier(loss='log', class_weight=classweights, shuffle=False)\n",
    "            loss = OnlineDRO.EasyAcc()\n",
    "            sincelast = OnlineDRO.EasyAcc()\n",
    "            blocksize = 32\n",
    "\n",
    "            for pno in range(1):\n",
    "                order = np.random.RandomState(seed=42+pno).permutation(len(train.data))\n",
    "                for n, ind in enumerate(zip(*(iter(order),)*blocksize)):\n",
    "                    v = np.array([ np.outer(t, np.append(t, [1])).ravel() for z in ind for t in ( train.data[z], ) ])\n",
    "                    actual = [ train.target[z] for z in ind ]\n",
    "                    if n > 0:\n",
    "                        pred = cls.predict(v)\n",
    "                        for p, a in zip(pred, actual):\n",
    "                            loss += 0 if p == a else 1\n",
    "                            sincelast += 0 if p == a else 1\n",
    "                        if (n & (n - 1) == 0): # and n & 0xAAAAAAAA == 0):\n",
    "                            print('{:<8d}\\t{:<8.3f}\\t{:<10.3f}'.format(loss.n, loss.mean(), sincelast.mean()), flush=True)\n",
    "                            sincelast = OnlineDRO.EasyAcc()\n",
    "\n",
    "                    cls.partial_fit(v, actual, classes=classes)\n",
    "\n",
    "                print('{:<8d}\\t{:<8.3f}\\t{:<10.3f}'.format(loss.n, loss.mean(), sincelast.mean()), flush=True)\n",
    "                sincelast = OnlineDRO.EasyAcc()\n",
    "\n",
    "                preds = cls.predict(np.array([np.outer(d, np.append(d, [1])).ravel() for d in test.data]))\n",
    "                ascores = []\n",
    "                for b in range(16):\n",
    "                    bootie = np.random.RandomState(90210+b).choice(len(test.target), replace=True, size=len(test.target))\n",
    "                    ascores.append(accuracy_score(y_true=test.target[bootie], y_pred=preds[bootie]))\n",
    "\n",
    "                print(\"test accuracy: {}\".format(np.quantile(ascores, [0.05, 0.5, 0.95])))\n",
    "    \n",
    "OnlineDRO.flass()      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Partial Information Online Learning, Softmax Logging Policy\n",
    "\n",
    "Uniform $(\\tau = 0)$ and softmax $(\\tau = 4)$ are pretty similar for off-policy learning but uniform has larger regret."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "code_folding": [
     0,
     1,
     54
    ],
    "hidden": true,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** lr = 0.0031622776601683794 tau = 0 ***\n",
      "n       \temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.469   \t0.469     \t0.143   \t0.143     \n",
      "64      \t0.484   \t0.500     \t0.143   \t0.143     \n",
      "128     \t0.492   \t0.500     \t0.143   \t0.143     \n",
      "256     \t0.508   \t0.523     \t0.143   \t0.143     \n",
      "512     \t0.490   \t0.473     \t0.143   \t0.143     \n",
      "1024    \t0.461   \t0.432     \t0.143   \t0.143     \n",
      "2048    \t0.448   \t0.435     \t0.143   \t0.143     \n",
      "4096    \t0.431   \t0.414     \t0.143   \t0.143     \n",
      "8192    \t0.412   \t0.394     \t0.143   \t0.143     \n",
      "16384   \t0.387   \t0.362     \t0.143   \t0.143     \n",
      "32768   \t0.370   \t0.353     \t0.143   \t0.143     \n",
      "65536   \t0.348   \t0.326     \t0.143   \t0.143     \n",
      "116160  \t0.331   \t0.309     \t0.143   \t0.143     \n",
      "test accuracy: [0.69911189 0.69984876 0.70084917]\n",
      "*** lr = 0.0031622776601683794 tau = 4 ***\n",
      "n       \temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t0.500   \t0.531     \t0.239   \t0.237     \n",
      "128     \t0.484   \t0.469     \t0.275   \t0.311     \n",
      "256     \t0.527   \t0.570     \t0.299   \t0.323     \n",
      "512     \t0.510   \t0.492     \t0.354   \t0.409     \n",
      "1024    \t0.510   \t0.510     \t0.385   \t0.415     \n",
      "2048    \t0.499   \t0.487     \t0.412   \t0.440     \n",
      "4096    \t0.479   \t0.460     \t0.436   \t0.459     \n",
      "8192    \t0.460   \t0.440     \t0.455   \t0.474     \n",
      "16384   \t0.426   \t0.392     \t0.483   \t0.511     \n",
      "32768   \t0.403   \t0.381     \t0.495   \t0.507     \n",
      "65536   \t0.373   \t0.342     \t0.506   \t0.518     \n",
      "116160  \t0.350   \t0.320     \t0.511   \t0.516     \n",
      "test accuracy: [0.69698575 0.69779522 0.69884297]\n"
     ]
    }
   ],
   "source": [
    "class OnlineDRO:\n",
    "    class EasyAcc:\n",
    "        def __init__(self):\n",
    "            self.n = 0\n",
    "            self.sum = 0\n",
    "            \n",
    "        def __iadd__(self, other):\n",
    "            self.n += 1\n",
    "            self.sum += other\n",
    "            return self\n",
    "            \n",
    "        def mean(self):\n",
    "            return self.sum / max(self.n, 1)\n",
    " \n",
    "    def flass():\n",
    "        from scipy.special import softmax\n",
    "        from sklearn.datasets import fetch_covtype\n",
    "        from sklearn.decomposition import PCA\n",
    "        from sklearn.linear_model import SGDClassifier\n",
    "        from sklearn.metrics import accuracy_score\n",
    "        from math import ceil\n",
    "        import numpy as np\n",
    "        \n",
    "        cov = fetch_covtype()\n",
    "        cov.data = PCA(whiten=True).fit_transform(cov.data)\n",
    "        classes = np.unique(cov.target - 1)\n",
    "        ndata = len(cov.target)\n",
    "        order = np.random.RandomState(seed=42).permutation(ndata)\n",
    "        ntrain = ceil(0.2 * ndata)\n",
    "        Object = lambda **kwargs: type(\"Object\", (), kwargs)()\n",
    "        train = Object(data = cov.data[order[:ntrain]], target = cov.target[order[:ntrain]] - 1)\n",
    "        test = Object(data = cov.data[order[ntrain:]], target = cov.target[order[ntrain:]] - 1)\n",
    "        \n",
    "        blocksize = 32\n",
    "        for lr, tau in ( (x, y) for x in np.logspace(-2.5, -2, 1) for y in (0, 4, ) ):\n",
    "            print(\"*** lr = {} tau = {} ***\".format(lr, tau), flush=True)\n",
    "            print('{:8.8s}\\t{:8.8s}\\t{:10.10s}\\t{:8.8s}\\t{:10.10s}'.format(\n",
    "                'n', 'emp loss', 'since last', 'log pv', 'since last')\n",
    "            )\n",
    "            \n",
    "            cls = SGDClassifier(loss='log', shuffle=False)\n",
    "            loss = OnlineDRO.EasyAcc()\n",
    "            sincelast = OnlineDRO.EasyAcc()\n",
    "            logpv = OnlineDRO.EasyAcc()\n",
    "            logpvsl = OnlineDRO.EasyAcc()\n",
    "            \n",
    "            loggerrand = np.random.RandomState(seed=2112)\n",
    "            logchoices = [None]*len(train.data)\n",
    "            pchoices = [None]*len(train.data)\n",
    " \n",
    "            for pno in range(1):\n",
    "                order = np.random.RandomState(seed=42+pno).permutation(len(train.data))\n",
    "                for n, ind in enumerate(zip(*(iter(order),)*blocksize)):\n",
    "                    v = np.array([ np.outer(t, np.append(t, [1])).ravel() for z in ind for t in ( train.data[z], ) ]) \n",
    "                    if n == 0 and pno == 0:\n",
    "                        for i, z in enumerate(ind):\n",
    "                            if logchoices[z] is None:\n",
    "                                choice = loggerrand.choice(a=classes, size=1)\n",
    "                                logchoices[z] = choice[0]\n",
    "                                pchoices[z] = 1.0 / len(classes)\n",
    "                    else:\n",
    "                        predlogp = cls.predict_proba(v)\n",
    "                        soft = softmax(tau * predlogp, axis=1)\n",
    "\n",
    "                        for i, z in enumerate(ind):\n",
    "                            if logchoices[z] is None:\n",
    "                                choice = loggerrand.choice(a=classes, p=soft[i,:], size=1)\n",
    "                                logchoices[z] = choice[0]\n",
    "                                pchoices[z] = soft[i, choice[0]]\n",
    "\n",
    "                        pred = cls.predict(v)\n",
    "                        actual = [ train.target[z] for z in ind ]\n",
    "                        for i, (p, a) in enumerate(zip(pred, actual)):\n",
    "                            loss += 0 if p == a else 1\n",
    "                            sincelast += 0 if p == a else 1\n",
    "                            logpv += soft[i, a]\n",
    "                            logpvsl += soft[i, a]\n",
    "                            \n",
    "                        if (n & (n - 1) == 0): # and n & 0xAAAAAAAA == 0):\n",
    "                            print('{:<8d}\\t{:<8.3f}\\t{:<10.3f}\\t{:<8.3f}\\t{:<10.3f}'.format(\n",
    "                                        loss.n, loss.mean(), sincelast.mean(), logpv.mean(), logpvsl.mean()),\n",
    "                                  flush=True)\n",
    "\n",
    "                            sincelast = OnlineDRO.EasyAcc()\n",
    "                            logpvsl = OnlineDRO.EasyAcc()\n",
    "\n",
    "                    x = np.array([ v[i] for i, z in enumerate(ind) if logchoices[z] == train.target[z] ])\n",
    "                    y = np.array([ logchoices[z] for i, z in enumerate(ind) if logchoices[z] == train.target[z] ])\n",
    "                    w = np.array([ (lr / len(classes)) * (1 / pchoices[z]) \n",
    "                                   for i, z in enumerate(ind) if logchoices[z] == train.target[z] ])\n",
    "                    if np.any(x):\n",
    "                        cls.partial_fit(x, y, classes=classes, sample_weight=w)\n",
    "\n",
    "                print('{:<8d}\\t{:<8.3f}\\t{:<10.3f}\\t{:<8.3f}\\t{:<10.3f}'.format(\n",
    "                             loss.n, loss.mean(), sincelast.mean(), logpv.mean(), logpvsl.mean()),\n",
    "                       flush=True)                \n",
    "                sincelast = OnlineDRO.EasyAcc()\n",
    "                logpvsl = OnlineDRO.EasyAcc()\n",
    "\n",
    "                preds = cls.predict(np.array([np.outer(d, np.append(d, [1])).ravel() for d in test.data]))\n",
    "                ascores = []\n",
    "                for b in range(16):\n",
    "                    bootie = np.random.RandomState(90210+b).choice(len(test.target), replace=True, size=len(test.target))\n",
    "                    ascores.append(accuracy_score(y_true=test.target[bootie], y_pred=preds[bootie]))\n",
    "\n",
    "                print(\"test accuracy: {}\".format(np.quantile(ascores, [0.05, 0.5, 0.95])))\n",
    "    \n",
    "OnlineDRO.flass()      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Bound Online Learning\n",
    "\n",
    "Either MLE, lower bound, or upper bound.  Some count decay seems better than no count decay.  Upper bound seems to want less count decay than the MLE or lower bound.  All forms of bound learning have lower regret than IPS learning.\n",
    "\n",
    "\n",
    "TODO: \"delayed batch\" online learning."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "code_folding": [
     1,
     14,
     17,
     34,
     52,
     59,
     66,
     74,
     82,
     90,
     141,
     150
    ],
    "hidden": true,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** lr = 0.0031622776601683794 tau = 4 what = mle gamma = 1 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.094   \t0.094     \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t1.024   \t1.954     \t0.469   \t0.469     \t0.305   \t0.369     \n",
      "128     \t1.023   \t1.022     \t0.477   \t0.484     \t0.353   \t0.401     \n",
      "256     \t0.930   \t0.836     \t0.488   \t0.500     \t0.391   \t0.430     \n",
      "512     \t0.823   \t0.716     \t0.473   \t0.457     \t0.429   \t0.467     \n",
      "1024    \t0.811   \t0.800     \t0.450   \t0.428     \t0.469   \t0.510     \n",
      "2048    \t0.775   \t0.739     \t0.422   \t0.395     \t0.505   \t0.541     \n",
      "4096    \t0.742   \t0.709     \t0.414   \t0.405     \t0.520   \t0.535     \n",
      "8192    \t0.738   \t0.734     \t0.412   \t0.411     \t0.523   \t0.525     \n",
      "16384   \t0.743   \t0.748     \t0.394   \t0.375     \t0.539   \t0.556     \n",
      "32768   \t0.760   \t0.777     \t0.376   \t0.357     \t0.555   \t0.570     \n",
      "65536   \t0.797   \t0.834     \t0.354   \t0.333     \t0.570   \t0.585     \n",
      "116160  \t0.828   \t0.869     \t0.339   \t0.320     \t0.578   \t0.590     \n",
      "test accuracy: [0.68822947 0.68919707 0.69008829]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = mle gamma = 0.9999 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.094   \t0.094     \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t1.022   \t1.951     \t0.469   \t0.469     \t0.305   \t0.369     \n",
      "128     \t1.022   \t1.022     \t0.477   \t0.484     \t0.353   \t0.401     \n",
      "256     \t0.929   \t0.836     \t0.488   \t0.500     \t0.391   \t0.430     \n",
      "512     \t0.822   \t0.714     \t0.473   \t0.457     \t0.429   \t0.467     \n",
      "1024    \t0.871   \t0.919     \t0.456   \t0.439     \t0.464   \t0.500     \n",
      "2048    \t0.803   \t0.736     \t0.419   \t0.382     \t0.504   \t0.544     \n",
      "4096    \t0.784   \t0.766     \t0.405   \t0.391     \t0.523   \t0.542     \n",
      "8192    \t0.756   \t0.729     \t0.394   \t0.383     \t0.539   \t0.554     \n",
      "16384   \t0.763   \t0.770     \t0.382   \t0.370     \t0.549   \t0.560     \n",
      "32768   \t0.770   \t0.778     \t0.371   \t0.361     \t0.558   \t0.566     \n",
      "65536   \t0.804   \t0.838     \t0.361   \t0.351     \t0.562   \t0.567     \n",
      "116160  \t0.829   \t0.862     \t0.351   \t0.337     \t0.568   \t0.574     \n",
      "test accuracy: [0.66472626 0.66579283 0.66703474]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = mle gamma = 0.999 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.094   \t0.094     \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t1.009   \t1.924     \t0.469   \t0.469     \t0.305   \t0.368     \n",
      "128     \t1.011   \t1.014     \t0.477   \t0.484     \t0.352   \t0.400     \n",
      "256     \t0.926   \t0.841     \t0.492   \t0.508     \t0.391   \t0.429     \n",
      "512     \t0.816   \t0.705     \t0.471   \t0.449     \t0.437   \t0.484     \n",
      "1024    \t0.781   \t0.745     \t0.455   \t0.439     \t0.467   \t0.497     \n",
      "2048    \t0.775   \t0.770     \t0.448   \t0.441     \t0.484   \t0.500     \n",
      "4096    \t0.765   \t0.754     \t0.420   \t0.393     \t0.513   \t0.542     \n",
      "8192    \t0.753   \t0.741     \t0.407   \t0.393     \t0.529   \t0.545     \n",
      "16384   \t0.758   \t0.763     \t0.393   \t0.379     \t0.542   \t0.555     \n",
      "32768   \t0.759   \t0.761     \t0.387   \t0.380     \t0.545   \t0.549     \n",
      "65536   \t0.801   \t0.843     \t0.363   \t0.339     \t0.562   \t0.579     \n",
      "116160  \t0.845   \t0.903     \t0.336   \t0.302     \t0.581   \t0.604     \n",
      "test accuracy: [0.69481497 0.69572663 0.69657913]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = mle gamma = 0.99 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.094   \t0.094     \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t0.891   \t1.688     \t0.453   \t0.438     \t0.301   \t0.361     \n",
      "128     \t0.906   \t0.921     \t0.477   \t0.500     \t0.349   \t0.397     \n",
      "256     \t0.865   \t0.824     \t0.484   \t0.492     \t0.389   \t0.430     \n",
      "512     \t0.773   \t0.681     \t0.473   \t0.461     \t0.434   \t0.479     \n",
      "1024    \t0.737   \t0.702     \t0.456   \t0.439     \t0.466   \t0.499     \n",
      "2048    \t0.731   \t0.724     \t0.427   \t0.397     \t0.502   \t0.538     \n",
      "4096    \t0.726   \t0.722     \t0.408   \t0.389     \t0.524   \t0.545     \n",
      "8192    \t0.738   \t0.749     \t0.396   \t0.385     \t0.536   \t0.548     \n",
      "16384   \t0.738   \t0.739     \t0.386   \t0.375     \t0.547   \t0.558     \n",
      "32768   \t0.748   \t0.757     \t0.378   \t0.371     \t0.552   \t0.556     \n",
      "65536   \t0.779   \t0.811     \t0.364   \t0.350     \t0.560   \t0.569     \n",
      "116160  \t0.824   \t0.883     \t0.340   \t0.308     \t0.577   \t0.598     \n",
      "test accuracy: [0.70416343 0.70524237 0.70603624]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = ub gamma = 1 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.094   \t0.094     \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t3.228   \t6.362     \t0.547   \t0.625     \t0.279   \t0.317     \n",
      "128     \t1.981   \t0.733     \t0.492   \t0.438     \t0.347   \t0.415     \n",
      "256     \t1.404   \t0.828     \t0.539   \t0.586     \t0.359   \t0.371     \n",
      "512     \t1.138   \t0.872     \t0.562   \t0.586     \t0.367   \t0.374     \n",
      "1024    \t0.981   \t0.824     \t0.553   \t0.543     \t0.390   \t0.413     \n",
      "2048    \t0.863   \t0.746     \t0.505   \t0.458     \t0.429   \t0.469     \n",
      "4096    \t0.809   \t0.755     \t0.476   \t0.447     \t0.463   \t0.496     \n",
      "8192    \t0.789   \t0.769     \t0.455   \t0.434     \t0.485   \t0.507     \n",
      "16384   \t0.792   \t0.795     \t0.418   \t0.382     \t0.518   \t0.552     \n",
      "32768   \t0.806   \t0.819     \t0.385   \t0.353     \t0.547   \t0.576     \n",
      "65536   \t0.836   \t0.867     \t0.356   \t0.326     \t0.569   \t0.592     \n",
      "116160  \t0.863   \t0.898     \t0.334   \t0.307     \t0.583   \t0.600     \n",
      "test accuracy: [0.69934747 0.70008326 0.70105032]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = ub gamma = 0.9999 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.094   \t0.094     \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t3.224   \t6.355     \t0.547   \t0.625     \t0.279   \t0.317     \n",
      "128     \t1.979   \t0.733     \t0.492   \t0.438     \t0.347   \t0.415     \n",
      "256     \t1.407   \t0.835     \t0.543   \t0.594     \t0.359   \t0.371     \n",
      "512     \t1.110   \t0.813     \t0.561   \t0.578     \t0.370   \t0.381     \n",
      "1024    \t1.023   \t0.936     \t0.540   \t0.520     \t0.403   \t0.436     \n",
      "2048    \t0.857   \t0.691     \t0.500   \t0.460     \t0.433   \t0.463     \n",
      "4096    \t0.820   \t0.783     \t0.473   \t0.446     \t0.465   \t0.496     \n",
      "8192    \t0.802   \t0.784     \t0.446   \t0.418     \t0.492   \t0.518     \n",
      "16384   \t0.808   \t0.813     \t0.408   \t0.371     \t0.527   \t0.563     \n",
      "32768   \t0.831   \t0.854     \t0.371   \t0.333     \t0.558   \t0.589     \n",
      "65536   \t0.874   \t0.917     \t0.343   \t0.316     \t0.578   \t0.597     \n",
      "116160  \t0.896   \t0.924     \t0.330   \t0.312     \t0.584   \t0.592     \n",
      "test accuracy: [0.68947406 0.68995437 0.69088324]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = ub gamma = 0.999 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.094   \t0.094     \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t3.192   \t6.290     \t0.547   \t0.625     \t0.279   \t0.317     \n",
      "128     \t1.942   \t0.691     \t0.492   \t0.438     \t0.346   \t0.413     \n",
      "256     \t1.446   \t0.950     \t0.531   \t0.570     \t0.371   \t0.397     \n",
      "512     \t1.109   \t0.773     \t0.545   \t0.559     \t0.384   \t0.396     \n",
      "1024    \t1.054   \t0.998     \t0.534   \t0.523     \t0.404   \t0.425     \n",
      "2048    \t0.918   \t0.783     \t0.499   \t0.464     \t0.432   \t0.459     \n",
      "4096    \t0.928   \t0.938     \t0.468   \t0.436     \t0.467   \t0.503     \n",
      "8192    \t0.931   \t0.934     \t0.454   \t0.440     \t0.482   \t0.497     \n",
      "16384   \t0.938   \t0.945     \t0.424   \t0.394     \t0.511   \t0.540     \n",
      "32768   \t0.953   \t0.969     \t0.400   \t0.377     \t0.532   \t0.552     \n",
      "65536   \t0.967   \t0.980     \t0.390   \t0.379     \t0.541   \t0.551     \n",
      "116160  \t0.986   \t1.011     \t0.373   \t0.351     \t0.554   \t0.569     \n",
      "test accuracy: [0.66320628 0.66409966 0.66516354]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = ub gamma = 0.99 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.094   \t0.094     \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t2.905   \t5.717     \t0.531   \t0.594     \t0.279   \t0.316     \n",
      "128     \t1.788   \t0.671     \t0.500   \t0.469     \t0.339   \t0.400     \n",
      "256     \t1.448   \t1.108     \t0.520   \t0.539     \t0.376   \t0.414     \n",
      "512     \t1.486   \t1.525     \t0.531   \t0.543     \t0.395   \t0.414     \n",
      "1024    \t1.509   \t1.532     \t0.524   \t0.518     \t0.413   \t0.430     \n",
      "2048    \t1.279   \t1.049     \t0.494   \t0.463     \t0.427   \t0.442     \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4096    \t1.151   \t1.023     \t0.477   \t0.460     \t0.451   \t0.474     \n",
      "8192    \t1.100   \t1.050     \t0.467   \t0.456     \t0.466   \t0.482     \n",
      "16384   \t1.123   \t1.145     \t0.450   \t0.434     \t0.486   \t0.507     \n",
      "32768   \t1.110   \t1.097     \t0.428   \t0.407     \t0.505   \t0.524     \n",
      "65536   \t1.122   \t1.133     \t0.417   \t0.406     \t0.516   \t0.527     \n",
      "116160  \t1.132   \t1.144     \t0.405   \t0.390     \t0.527   \t0.540     \n",
      "test accuracy: [0.62157682 0.6224546  0.62311993]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = lb gamma = 1 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.094   \t0.094     \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t0.154   \t0.215     \t0.453   \t0.438     \t0.238   \t0.235     \n",
      "128     \t0.156   \t0.158     \t0.484   \t0.516     \t0.217   \t0.196     \n",
      "256     \t0.221   \t0.286     \t0.535   \t0.586     \t0.231   \t0.244     \n",
      "512     \t0.511   \t0.801     \t0.500   \t0.465     \t0.334   \t0.438     \n",
      "1024    \t0.616   \t0.722     \t0.501   \t0.502     \t0.385   \t0.436     \n",
      "2048    \t0.672   \t0.728     \t0.451   \t0.400     \t0.461   \t0.536     \n",
      "4096    \t0.710   \t0.748     \t0.420   \t0.390     \t0.501   \t0.541     \n",
      "8192    \t0.720   \t0.730     \t0.409   \t0.399     \t0.519   \t0.537     \n",
      "16384   \t0.747   \t0.774     \t0.381   \t0.353     \t0.549   \t0.579     \n",
      "32768   \t0.770   \t0.793     \t0.367   \t0.353     \t0.560   \t0.572     \n",
      "65536   \t0.803   \t0.835     \t0.353   \t0.339     \t0.569   \t0.577     \n",
      "116160  \t0.832   \t0.870     \t0.342   \t0.328     \t0.574   \t0.580     \n",
      "test accuracy: [0.68575963 0.68657664 0.68735814]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = lb gamma = 0.9999 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.094   \t0.094     \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t0.154   \t0.215     \t0.453   \t0.438     \t0.238   \t0.235     \n",
      "128     \t0.156   \t0.158     \t0.484   \t0.516     \t0.217   \t0.196     \n",
      "256     \t0.222   \t0.288     \t0.535   \t0.586     \t0.231   \t0.244     \n",
      "512     \t0.513   \t0.804     \t0.494   \t0.453     \t0.343   \t0.456     \n",
      "1024    \t0.610   \t0.706     \t0.494   \t0.494     \t0.395   \t0.447     \n",
      "2048    \t0.680   \t0.751     \t0.440   \t0.386     \t0.466   \t0.537     \n",
      "4096    \t0.732   \t0.783     \t0.396   \t0.351     \t0.520   \t0.573     \n",
      "8192    \t0.741   \t0.751     \t0.381   \t0.366     \t0.542   \t0.564     \n",
      "16384   \t0.766   \t0.791     \t0.379   \t0.376     \t0.547   \t0.551     \n",
      "32768   \t0.772   \t0.779     \t0.368   \t0.357     \t0.557   \t0.567     \n",
      "65536   \t0.795   \t0.818     \t0.354   \t0.339     \t0.568   \t0.580     \n",
      "116160  \t0.830   \t0.876     \t0.334   \t0.308     \t0.582   \t0.600     \n",
      "test accuracy: [0.69846216 0.69920548 0.70053398]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = lb gamma = 0.999 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.094   \t0.094     \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t0.153   \t0.213     \t0.453   \t0.438     \t0.238   \t0.235     \n",
      "128     \t0.154   \t0.155     \t0.484   \t0.516     \t0.217   \t0.196     \n",
      "256     \t0.218   \t0.282     \t0.531   \t0.578     \t0.234   \t0.251     \n",
      "512     \t0.522   \t0.827     \t0.463   \t0.395     \t0.366   \t0.498     \n",
      "1024    \t0.623   \t0.724     \t0.446   \t0.430     \t0.436   \t0.505     \n",
      "2048    \t0.657   \t0.692     \t0.425   \t0.403     \t0.485   \t0.535     \n",
      "4096    \t0.684   \t0.711     \t0.403   \t0.382     \t0.518   \t0.550     \n",
      "8192    \t0.695   \t0.706     \t0.398   \t0.393     \t0.530   \t0.543     \n",
      "16384   \t0.705   \t0.714     \t0.392   \t0.386     \t0.539   \t0.548     \n",
      "32768   \t0.722   \t0.739     \t0.384   \t0.376     \t0.546   \t0.554     \n",
      "65536   \t0.746   \t0.771     \t0.370   \t0.357     \t0.557   \t0.567     \n",
      "116160  \t0.770   \t0.802     \t0.356   \t0.338     \t0.566   \t0.579     \n",
      "test accuracy: [0.67352289 0.6747535  0.67572648]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = lb gamma = 0.99 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.094   \t0.094     \t0.469   \t0.469     \t0.241   \t0.241     \n",
      "64      \t0.144   \t0.195     \t0.453   \t0.438     \t0.238   \t0.235     \n",
      "128     \t0.140   \t0.136     \t0.484   \t0.516     \t0.217   \t0.196     \n",
      "256     \t0.228   \t0.315     \t0.527   \t0.570     \t0.239   \t0.261     \n",
      "512     \t0.480   \t0.732     \t0.484   \t0.441     \t0.359   \t0.479     \n",
      "1024    \t0.558   \t0.637     \t0.464   \t0.443     \t0.426   \t0.493     \n",
      "2048    \t0.608   \t0.657     \t0.433   \t0.402     \t0.476   \t0.527     \n",
      "4096    \t0.620   \t0.633     \t0.417   \t0.402     \t0.506   \t0.536     \n",
      "8192    \t0.633   \t0.646     \t0.415   \t0.412     \t0.517   \t0.528     \n",
      "16384   \t0.637   \t0.642     \t0.412   \t0.410     \t0.523   \t0.530     \n",
      "32768   \t0.643   \t0.649     \t0.411   \t0.409     \t0.527   \t0.530     \n",
      "65536   \t0.652   \t0.661     \t0.408   \t0.406     \t0.530   \t0.533     \n",
      "116160  \t0.661   \t0.672     \t0.400   \t0.388     \t0.538   \t0.549     \n",
      "test accuracy: [0.6130104  0.6143147  0.61559318]\n"
     ]
    }
   ],
   "source": [
    "class OnlineDRO:\n",
    "    class EasyAcc:\n",
    "        def __init__(self):\n",
    "            self.n = 0\n",
    "            self.sum = 0\n",
    "            \n",
    "        def __iadd__(self, other):\n",
    "            self.n += 1\n",
    "            self.sum += other\n",
    "            return self\n",
    "            \n",
    "        def mean(self):\n",
    "            return self.sum / max(self.n, 1)\n",
    "\n",
    "    class OnlineCressieReadLB:\n",
    "        from math import inf\n",
    "        \n",
    "        def __init__(self, alpha, gamma=1, wmin=0, wmax=inf):\n",
    "            import numpy as np\n",
    "            \n",
    "            self.alpha = alpha\n",
    "            self.gamma = gamma\n",
    "            self.n = 0\n",
    "            self.sumw = 0\n",
    "            self.sumwsq = 0\n",
    "            self.sumwr = 0\n",
    "            self.sumwsqr = 0\n",
    "            self.sumwsqrsq = 0\n",
    "            self.wmin = wmin\n",
    "            self.wmax = wmax\n",
    "            \n",
    "            self.duals = None\n",
    "            self.mleduals = None\n",
    "            \n",
    "        def update(self, c, w, r):\n",
    "            if c > 0:\n",
    "                assert w + 1e-6 >= self.wmin and w <= self.wmax + 1e-6, 'w = {} < {} < {}'.format(self.wmin, w, self.wmax)\n",
    "                assert r >= 0 and r <= 1, 'r = {}'.format(r)\n",
    "                \n",
    "                decay = self.gamma ** c\n",
    "                self.n = decay * self.n + c\n",
    "                self.sumw = decay * self.sumw + c * w\n",
    "                self.sumwsq = decay * self.sumwsq + c * w**2\n",
    "                self.sumwr = decay * self.sumwr + c * w * r\n",
    "                self.sumwsqr = decay * self.sumwsqr + c * (w**2) * r\n",
    "                self.sumwsqrsq = decay * self.sumwsqrsq + c * (w**2) * (r**2)\n",
    "                    \n",
    "                self.duals = None\n",
    "                self.mleduals = None\n",
    "                \n",
    "            return self\n",
    "        \n",
    "        def recomputeduals(self):\n",
    "            from MLE.MLE import CrMinusTwo as CrMinusTwo\n",
    "            \n",
    "            self.duals = CrMinusTwo.intervalimpl(self.n, self.sumw, self.sumwsq, \n",
    "                                                 self.sumwr, self.sumwsqr, self.sumwsqrsq,\n",
    "                                                 self.wmin, self.wmax, self.alpha, raiseonerr=True)\n",
    "            \n",
    "        def recomputedualsmle(self):\n",
    "            from MLE.MLE import CrMinusTwo as CrMinusTwo\n",
    "            \n",
    "            self.mleduals = CrMinusTwo.estimateimpl(self.n, self.sumw, self.sumwsq, \n",
    "                                                    self.sumwr, self.sumwsqr, None, None,\n",
    "                                                    self.wmin, self.wmax, raiseonerr=True)\n",
    "        \n",
    "        def qlb(self, c, w, r):\n",
    "            if self.duals is None:\n",
    "                self.recomputeduals()\n",
    "                \n",
    "                assert self.duals is not None\n",
    "                \n",
    "            return self.duals[1][0]['qfunc'](c, w, r) if self.duals[1][0] is not None else 1\n",
    "        \n",
    "        def qub(self, c, w, r):\n",
    "            if self.duals is None:\n",
    "                self.recomputeduals()\n",
    "                \n",
    "                assert self.duals is not None\n",
    "                \n",
    "            return self.duals[1][1]['qfunc'](c, w, r) if self.duals[1][1] is not None else 1\n",
    "        \n",
    "        def qmle(self, c, w, r):\n",
    "            if self.mleduals is None:\n",
    "                self.recomputedualsmle()\n",
    "                \n",
    "                assert self.mleduals is not None\n",
    "                \n",
    "            return self.mleduals[1]['qfunc'](c, w, r) if self.mleduals[1] is not None else 1\n",
    "           \n",
    "    def flass():\n",
    "        from scipy.special import softmax\n",
    "        from sklearn.datasets import fetch_covtype\n",
    "        from sklearn.decomposition import PCA\n",
    "        from sklearn.linear_model import SGDClassifier\n",
    "        from sklearn.metrics import accuracy_score\n",
    "        from math import ceil, exp\n",
    "        import numpy as np\n",
    "        \n",
    "        cov = fetch_covtype()\n",
    "        cov.data = PCA(whiten=True).fit_transform(cov.data)\n",
    "        classes = np.unique(cov.target - 1)\n",
    "        ndata = len(cov.target)\n",
    "        order = np.random.RandomState(seed=42).permutation(ndata)\n",
    "        ntrain = ceil(0.2 * ndata)\n",
    "        Object = lambda **kwargs: type(\"Object\", (), kwargs)()\n",
    "        train = Object(data = cov.data[order[:ntrain]], target = cov.target[order[:ntrain]] - 1)\n",
    "        test = Object(data = cov.data[order[ntrain:]], target = cov.target[order[ntrain:]] - 1)\n",
    "        \n",
    "        blocksize = 32\n",
    "        for lr, tau, what, gamma in ( (x, y, z, g) for x in np.logspace(-2.5, -2, 1) for y in (4, ) \n",
    "                                      for z in ('mle', 'ub', 'lb') for g in (1, 0.9999, 0.999, 0.99, ) ):\n",
    "            print(\"*** lr = {} tau = {} what = {} gamma = {} ***\".format(lr, tau, what, gamma), flush=True)\n",
    "            print('{:8.8s}\\t{:8.8s}\\t{:10.10s}\\t{:8.8s}\\t{:10.10s}\\t{:8.8s}\\t{:10.10s}'.format(\n",
    "                'n', 'eff n', 'since last', 'emp loss', 'since last', 'log pv', 'since last')\n",
    "            )\n",
    "            \n",
    "            cls = SGDClassifier(loss='log', shuffle=False)\n",
    "            loss = OnlineDRO.EasyAcc()\n",
    "            sincelast = OnlineDRO.EasyAcc()\n",
    "            logpv = OnlineDRO.EasyAcc()\n",
    "            logpvsl = OnlineDRO.EasyAcc()\n",
    "            effn = OnlineDRO.EasyAcc()\n",
    "            effnsl = OnlineDRO.EasyAcc()\n",
    "            \n",
    "            loggerrand = np.random.RandomState(seed=2112)\n",
    "            logchoices = [None]*len(train.data)\n",
    "            pchoices = [None]*len(train.data)\n",
    "                        \n",
    "            ocrl = OnlineDRO.OnlineCressieReadLB(alpha=0.05, \n",
    "                                                 gamma=gamma,\n",
    "                                                 wmin=0,\n",
    "                                                 wmax=exp(tau) + len(classes) - 1\n",
    "                                                )\n",
    "            qfunc = ocrl.qmle if what == 'mle' else ocrl.qlb if what == 'lb' else ocrl.qub\n",
    " \n",
    "            for pno in range(1):\n",
    "                order = np.random.RandomState(seed=42+pno).permutation(len(train.data))\n",
    "                for n, ind in enumerate(zip(*(iter(order),)*blocksize)):\n",
    "                    v = np.array([ np.outer(t, np.append(t, [1])).ravel() for z in ind for t in ( train.data[z], ) ]) \n",
    "                    if n == 0 and pno == 0:\n",
    "                        for i, z in enumerate(ind):\n",
    "                            if logchoices[z] is None:\n",
    "                                choice = loggerrand.choice(a=classes, size=1)\n",
    "                                logchoices[z] = choice[0]\n",
    "                                pchoices[z] = 1.0 / len(classes)\n",
    "                    else:\n",
    "                        predlogp = cls.predict_proba(v)\n",
    "                        soft = softmax(tau * predlogp, axis=1)\n",
    "\n",
    "                        for i, z in enumerate(ind):\n",
    "                            if logchoices[z] is None:\n",
    "                                choice = loggerrand.choice(a=classes, p=soft[i,:], size=1)\n",
    "                                logchoices[z] = choice[0]\n",
    "                                pchoices[z] = soft[i, choice[0]]\n",
    "\n",
    "                        pred = cls.predict(v)\n",
    "                        actual = [ train.target[z] for z in ind ]\n",
    "                        for i, (p, a) in enumerate(zip(pred, actual)):\n",
    "                            loss += 0 if p == a else 1\n",
    "                            sincelast += 0 if p == a else 1\n",
    "                            logpv += soft[i, a]\n",
    "                            logpvsl += soft[i, a]\n",
    "                            \n",
    "                        if (n & (n - 1) == 0): # and n & 0xAAAAAAAA == 0):\n",
    "                            print('{:<8d}\\t{:<8.3f}\\t{:<10.3f}\\t{:<8.3f}\\t{:<10.3f}\\t{:<8.3f}\\t{:<10.3f}'.format(\n",
    "                                         loss.n, effn.mean(), effnsl.mean(), loss.mean(), sincelast.mean(), logpv.mean(), logpvsl.mean()),\n",
    "                                   flush=True)   \n",
    "\n",
    "                            sincelast = OnlineDRO.EasyAcc()\n",
    "                            logpvsl = OnlineDRO.EasyAcc()\n",
    "                            effnsl = OnlineDRO.EasyAcc()\n",
    "                            \n",
    "                    for i, z in enumerate(ind):\n",
    "                        r = 1 if logchoices[z] == train.target[z] else 0\n",
    "                        w = 1 / pchoices[z] \n",
    "                        ocrl.update(1, w, r) \n",
    "                        \n",
    "                    if n == 0 and pno == 0:\n",
    "                        sampweight = np.array([ lr for i, z in enumerate(ind) if logchoices[z] == train.target[z] ])\n",
    "                    else:   \n",
    "                        sampweight = np.array([ lr * w * ocrl.n \n",
    "                                                   * max(0, qfunc(1, w, 1))\n",
    "                                                for i, z in enumerate(ind) \n",
    "                                                if logchoices[z] == train.target[z]\n",
    "                                                for w in (1 / pchoices[z],)# if logchoices[z] == pred[i] else 0,)\n",
    "                                              ])\n",
    "                        \n",
    "                    effn += sampweight.sum() / (lr * blocksize)\n",
    "                    effnsl += sampweight.sum() / (lr * blocksize)\n",
    "                                            \n",
    "                    x = np.array([ v[i] for i, z in enumerate(ind) if logchoices[z] == train.target[z] ])\n",
    "                    y = np.array([ logchoices[z] for i, z in enumerate(ind) if logchoices[z] == train.target[z] ])\n",
    "                    \n",
    "                    if np.any(x):\n",
    "                        cls.partial_fit(x, y, classes=classes, sample_weight=sampweight)\n",
    "\n",
    "                print('{:<8d}\\t{:<8.3f}\\t{:<10.3f}\\t{:<8.3f}\\t{:<10.3f}\\t{:<8.3f}\\t{:<10.3f}'.format(\n",
    "                             loss.n, effn.mean(), effnsl.mean(), loss.mean(), sincelast.mean(), logpv.mean(), logpvsl.mean()),\n",
    "                       flush=True)\n",
    "\n",
    "                sincelast = OnlineDRO.EasyAcc()\n",
    "                logpvsl = OnlineDRO.EasyAcc()\n",
    "                effnsl = OnlineDRO.EasyAcc()\n",
    "\n",
    "                preds = cls.predict(np.array([np.outer(d, np.append(d, [1])).ravel() for d in test.data]))\n",
    "                ascores = []\n",
    "                for b in range(16):\n",
    "                    bootie = np.random.RandomState(90210+b).choice(len(test.target), replace=True, size=len(test.target))\n",
    "                    ascores.append(accuracy_score(y_true=test.target[bootie], y_pred=preds[bootie]))\n",
    "\n",
    "                print(\"test accuracy: {}\".format(np.quantile(ascores, [0.05, 0.5, 0.95])))\n",
    "    \n",
    "OnlineDRO.flass()      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "### Delayed Batch Bound Online Learning\n",
    "\n",
    "Process larger batches as sets of smaller batches to emulate time delay in policy updates without changing optimization properties (i.e., SGD batch size).\n",
    "\n",
    "Everything still seems to work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "code_folding": [
     1,
     2,
     6,
     11,
     14,
     17,
     34,
     52,
     59,
     66,
     74,
     82,
     90,
     142
    ],
    "hidden": true,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** lr = 0.0031622776601683794 tau = 4 what = mle gamma = 1 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "256     \t0.141   \t0.141     \t0.422   \t0.422     \t0.452   \t0.452     \n",
      "512     \t0.481   \t0.822     \t0.451   \t0.480     \t0.458   \t0.463     \n",
      "1024    \t0.541   \t0.601     \t0.444   \t0.438     \t0.475   \t0.493     \n",
      "2048    \t0.596   \t0.651     \t0.424   \t0.404     \t0.506   \t0.537     \n",
      "4096    \t0.641   \t0.685     \t0.410   \t0.396     \t0.521   \t0.536     \n",
      "8192    \t0.645   \t0.649     \t0.395   \t0.381     \t0.537   \t0.553     \n",
      "16384   \t0.667   \t0.690     \t0.386   \t0.377     \t0.544   \t0.552     \n",
      "32768   \t0.683   \t0.698     \t0.367   \t0.347     \t0.562   \t0.579     \n",
      "65536   \t0.700   \t0.717     \t0.353   \t0.340     \t0.571   \t0.579     \n",
      "115712  \t0.725   \t0.758     \t0.344   \t0.331     \t0.574   \t0.579     \n",
      "test accuracy: [0.67735833 0.67855076 0.67900847]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = mle gamma = 0.9999 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "256     \t0.141   \t0.141     \t0.422   \t0.422     \t0.452   \t0.452     \n",
      "512     \t0.481   \t0.822     \t0.449   \t0.477     \t0.458   \t0.463     \n",
      "1024    \t0.541   \t0.601     \t0.443   \t0.438     \t0.474   \t0.491     \n",
      "2048    \t0.587   \t0.634     \t0.420   \t0.396     \t0.505   \t0.536     \n",
      "4096    \t0.623   \t0.658     \t0.415   \t0.410     \t0.516   \t0.526     \n",
      "8192    \t0.633   \t0.643     \t0.412   \t0.409     \t0.521   \t0.527     \n",
      "16384   \t0.654   \t0.675     \t0.396   \t0.381     \t0.536   \t0.550     \n",
      "32768   \t0.656   \t0.659     \t0.393   \t0.391     \t0.537   \t0.539     \n",
      "65536   \t0.685   \t0.713     \t0.375   \t0.356     \t0.551   \t0.565     \n",
      "115712  \t0.715   \t0.755     \t0.355   \t0.329     \t0.565   \t0.582     \n",
      "test accuracy: [0.68849302 0.68939823 0.6900614 ]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = mle gamma = 0.999 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "256     \t0.141   \t0.141     \t0.422   \t0.422     \t0.452   \t0.452     \n",
      "512     \t0.480   \t0.820     \t0.447   \t0.473     \t0.456   \t0.459     \n",
      "1024    \t0.535   \t0.589     \t0.438   \t0.428     \t0.480   \t0.503     \n",
      "2048    \t0.593   \t0.651     \t0.406   \t0.375     \t0.518   \t0.556     \n",
      "4096    \t0.622   \t0.651     \t0.403   \t0.400     \t0.525   \t0.533     \n",
      "8192    \t0.643   \t0.665     \t0.409   \t0.414     \t0.523   \t0.520     \n",
      "16384   \t0.644   \t0.645     \t0.401   \t0.392     \t0.533   \t0.543     \n",
      "32768   \t0.669   \t0.694     \t0.378   \t0.356     \t0.550   \t0.568     \n",
      "65536   \t0.697   \t0.725     \t0.360   \t0.341     \t0.564   \t0.578     \n",
      "115712  \t0.728   \t0.769     \t0.344   \t0.325     \t0.574   \t0.586     \n",
      "test accuracy: [0.69617682 0.6970788  0.69786084]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = mle gamma = 0.99 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "256     \t0.141   \t0.141     \t0.422   \t0.422     \t0.452   \t0.452     \n",
      "512     \t0.483   \t0.825     \t0.436   \t0.449     \t0.467   \t0.481     \n",
      "1024    \t0.554   \t0.626     \t0.411   \t0.387     \t0.502   \t0.538     \n",
      "2048    \t0.604   \t0.654     \t0.393   \t0.375     \t0.531   \t0.560     \n",
      "4096    \t0.625   \t0.647     \t0.406   \t0.418     \t0.526   \t0.522     \n",
      "8192    \t0.631   \t0.636     \t0.406   \t0.406     \t0.528   \t0.530     \n",
      "16384   \t0.642   \t0.652     \t0.406   \t0.406     \t0.529   \t0.530     \n",
      "32768   \t0.655   \t0.669     \t0.392   \t0.377     \t0.540   \t0.550     \n",
      "65536   \t0.677   \t0.700     \t0.373   \t0.354     \t0.553   \t0.567     \n",
      "115712  \t0.707   \t0.745     \t0.357   \t0.336     \t0.563   \t0.577     \n",
      "test accuracy: [0.68596886 0.68669927 0.68783414]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = ub gamma = 1 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "256     \t0.141   \t0.141     \t0.422   \t0.422     \t0.452   \t0.452     \n",
      "512     \t0.581   \t1.022     \t0.434   \t0.445     \t0.471   \t0.490     \n",
      "1024    \t0.641   \t0.700     \t0.426   \t0.418     \t0.496   \t0.522     \n",
      "2048    \t0.648   \t0.656     \t0.411   \t0.396     \t0.520   \t0.543     \n",
      "4096    \t0.694   \t0.739     \t0.424   \t0.437     \t0.510   \t0.501     \n",
      "8192    \t0.691   \t0.688     \t0.431   \t0.438     \t0.504   \t0.497     \n",
      "16384   \t0.694   \t0.697     \t0.405   \t0.379     \t0.529   \t0.554     \n",
      "32768   \t0.701   \t0.709     \t0.377   \t0.349     \t0.552   \t0.575     \n",
      "65536   \t0.715   \t0.728     \t0.357   \t0.337     \t0.567   \t0.582     \n",
      "115712  \t0.741   \t0.775     \t0.342   \t0.322     \t0.576   \t0.587     \n",
      "test accuracy: [0.68314404 0.68403151 0.68514487]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = ub gamma = 0.9999 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "256     \t0.141   \t0.141     \t0.422   \t0.422     \t0.452   \t0.452     \n",
      "512     \t0.583   \t1.025     \t0.434   \t0.445     \t0.471   \t0.490     \n",
      "1024    \t0.646   \t0.709     \t0.424   \t0.414     \t0.496   \t0.521     \n",
      "2048    \t0.653   \t0.660     \t0.410   \t0.396     \t0.518   \t0.539     \n",
      "4096    \t0.694   \t0.735     \t0.428   \t0.446     \t0.505   \t0.492     \n",
      "8192    \t0.672   \t0.649     \t0.441   \t0.455     \t0.497   \t0.489     \n",
      "16384   \t0.678   \t0.683     \t0.421   \t0.401     \t0.515   \t0.533     \n",
      "32768   \t0.703   \t0.729     \t0.394   \t0.366     \t0.537   \t0.558     \n",
      "65536   \t0.718   \t0.733     \t0.378   \t0.362     \t0.548   \t0.559     \n",
      "115712  \t0.750   \t0.792     \t0.356   \t0.327     \t0.563   \t0.583     \n",
      "test accuracy: [0.68395083 0.6848996  0.6861399 ]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = ub gamma = 0.999 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "256     \t0.141   \t0.141     \t0.422   \t0.422     \t0.452   \t0.452     \n",
      "512     \t0.598   \t1.056     \t0.436   \t0.449     \t0.472   \t0.492     \n",
      "1024    \t0.698   \t0.798     \t0.427   \t0.418     \t0.484   \t0.497     \n",
      "2048    \t0.735   \t0.772     \t0.418   \t0.410     \t0.504   \t0.525     \n",
      "4096    \t0.797   \t0.858     \t0.440   \t0.461     \t0.485   \t0.466     \n",
      "8192    \t0.786   \t0.776     \t0.455   \t0.470     \t0.477   \t0.468     \n",
      "16384   \t0.791   \t0.796     \t0.452   \t0.448     \t0.484   \t0.492     \n",
      "32768   \t0.806   \t0.821     \t0.435   \t0.419     \t0.501   \t0.518     \n",
      "65536   \t0.815   \t0.825     \t0.425   \t0.414     \t0.512   \t0.524     \n",
      "115712  \t0.837   \t0.866     \t0.406   \t0.381     \t0.529   \t0.550     \n",
      "test accuracy: [0.61968464 0.62079263 0.621828  ]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = ub gamma = 0.99 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "256     \t0.141   \t0.141     \t0.422   \t0.422     \t0.452   \t0.452     \n",
      "512     \t0.571   \t1.001     \t0.436   \t0.449     \t0.471   \t0.490     \n",
      "1024    \t0.915   \t1.259     \t0.454   \t0.473     \t0.461   \t0.450     \n",
      "2048    \t0.833   \t0.751     \t0.457   \t0.460     \t0.469   \t0.478     \n",
      "4096    \t0.923   \t1.013     \t0.480   \t0.503     \t0.448   \t0.428     \n",
      "8192    \t0.951   \t0.978     \t0.487   \t0.494     \t0.444   \t0.439     \n",
      "16384   \t0.944   \t0.937     \t0.471   \t0.456     \t0.462   \t0.480     \n",
      "32768   \t0.955   \t0.966     \t0.462   \t0.453     \t0.473   \t0.484     \n",
      "65536   \t0.975   \t0.995     \t0.444   \t0.427     \t0.493   \t0.512     \n",
      "115712  \t0.980   \t0.988     \t0.424   \t0.398     \t0.511   \t0.535     \n",
      "test accuracy: [0.6197153  0.62073454 0.62161608]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = lb gamma = 1 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "256     \t0.141   \t0.141     \t0.422   \t0.422     \t0.452   \t0.452     \n",
      "512     \t0.427   \t0.713     \t0.461   \t0.500     \t0.446   \t0.439     \n",
      "1024    \t0.513   \t0.598     \t0.434   \t0.406     \t0.477   \t0.508     \n",
      "2048    \t0.587   \t0.662     \t0.415   \t0.396     \t0.506   \t0.535     \n",
      "4096    \t0.614   \t0.641     \t0.423   \t0.432     \t0.505   \t0.504     \n",
      "8192    \t0.628   \t0.643     \t0.408   \t0.393     \t0.523   \t0.541     \n",
      "16384   \t0.663   \t0.697     \t0.384   \t0.361     \t0.544   \t0.564     \n",
      "32768   \t0.676   \t0.690     \t0.363   \t0.341     \t0.563   \t0.583     \n",
      "65536   \t0.691   \t0.705     \t0.351   \t0.339     \t0.572   \t0.581     \n",
      "115712  \t0.714   \t0.744     \t0.341   \t0.328     \t0.577   \t0.584     \n",
      "test accuracy: [0.68785082 0.68856132 0.69022921]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = lb gamma = 0.9999 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "256     \t0.141   \t0.141     \t0.422   \t0.422     \t0.452   \t0.452     \n",
      "512     \t0.426   \t0.712     \t0.461   \t0.500     \t0.446   \t0.439     \n",
      "1024    \t0.511   \t0.596     \t0.435   \t0.408     \t0.477   \t0.508     \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2048    \t0.585   \t0.658     \t0.416   \t0.396     \t0.506   \t0.535     \n",
      "4096    \t0.607   \t0.630     \t0.415   \t0.415     \t0.513   \t0.520     \n",
      "8192    \t0.638   \t0.668     \t0.397   \t0.379     \t0.534   \t0.554     \n",
      "16384   \t0.658   \t0.678     \t0.386   \t0.375     \t0.545   \t0.557     \n",
      "32768   \t0.678   \t0.698     \t0.363   \t0.341     \t0.563   \t0.582     \n",
      "65536   \t0.677   \t0.677     \t0.365   \t0.366     \t0.563   \t0.562     \n",
      "115712  \t0.693   \t0.714     \t0.355   \t0.342     \t0.570   \t0.579     \n",
      "test accuracy: [0.66936365 0.67011073 0.6713677 ]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = lb gamma = 0.999 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "256     \t0.141   \t0.141     \t0.422   \t0.422     \t0.452   \t0.452     \n",
      "512     \t0.419   \t0.697     \t0.461   \t0.500     \t0.445   \t0.437     \n",
      "1024    \t0.501   \t0.584     \t0.436   \t0.410     \t0.477   \t0.510     \n",
      "2048    \t0.578   \t0.655     \t0.407   \t0.378     \t0.514   \t0.551     \n",
      "4096    \t0.603   \t0.628     \t0.411   \t0.415     \t0.518   \t0.522     \n",
      "8192    \t0.622   \t0.640     \t0.402   \t0.393     \t0.530   \t0.541     \n",
      "16384   \t0.646   \t0.671     \t0.377   \t0.352     \t0.553   \t0.576     \n",
      "32768   \t0.665   \t0.684     \t0.361   \t0.344     \t0.566   \t0.578     \n",
      "65536   \t0.683   \t0.701     \t0.349   \t0.337     \t0.573   \t0.581     \n",
      "115712  \t0.703   \t0.728     \t0.336   \t0.320     \t0.581   \t0.592     \n",
      "test accuracy: [0.68891039 0.69026848 0.69160128]\n",
      "*** lr = 0.0031622776601683794 tau = 4 what = lb gamma = 0.99 ***\n",
      "n       \teff n   \tsince last\temp loss\tsince last\tlog pv  \tsince last\n",
      "256     \t0.141   \t0.141     \t0.422   \t0.422     \t0.452   \t0.452     \n",
      "512     \t0.360   \t0.579     \t0.449   \t0.477     \t0.450   \t0.449     \n",
      "1024    \t0.434   \t0.507     \t0.433   \t0.416     \t0.479   \t0.508     \n",
      "2048    \t0.526   \t0.619     \t0.404   \t0.375     \t0.516   \t0.553     \n",
      "4096    \t0.571   \t0.615     \t0.397   \t0.391     \t0.528   \t0.541     \n",
      "8192    \t0.587   \t0.603     \t0.388   \t0.379     \t0.540   \t0.551     \n",
      "16384   \t0.608   \t0.630     \t0.379   \t0.370     \t0.548   \t0.557     \n",
      "32768   \t0.622   \t0.636     \t0.364   \t0.349     \t0.561   \t0.574     \n",
      "65536   \t0.642   \t0.663     \t0.347   \t0.331     \t0.574   \t0.588     \n",
      "115712  \t0.658   \t0.677     \t0.333   \t0.315     \t0.585   \t0.599     \n",
      "test accuracy: [0.68671917 0.68775777 0.68868503]\n"
     ]
    }
   ],
   "source": [
    "class OnlineDRO:\n",
    "    class EasyAcc:\n",
    "        def __init__(self):\n",
    "            self.n = 0\n",
    "            self.sum = 0\n",
    "            \n",
    "        def __iadd__(self, other):\n",
    "            self.n += 1\n",
    "            self.sum += other\n",
    "            return self\n",
    "            \n",
    "        def mean(self):\n",
    "            return self.sum / max(self.n, 1)\n",
    "\n",
    "    class OnlineCressieReadLB:\n",
    "        from math import inf\n",
    "        \n",
    "        def __init__(self, alpha, gamma=1, wmin=0, wmax=inf):\n",
    "            import numpy as np\n",
    "            \n",
    "            self.alpha = alpha\n",
    "            self.gamma = gamma\n",
    "            self.n = 0\n",
    "            self.sumw = 0\n",
    "            self.sumwsq = 0\n",
    "            self.sumwr = 0\n",
    "            self.sumwsqr = 0\n",
    "            self.sumwsqrsq = 0\n",
    "            self.wmin = wmin\n",
    "            self.wmax = wmax\n",
    "            \n",
    "            self.duals = None\n",
    "            self.mleduals = None\n",
    "            \n",
    "        def update(self, c, w, r):\n",
    "            if c > 0:\n",
    "                assert w + 1e-6 >= self.wmin and w <= self.wmax + 1e-6, 'w = {} < {} < {}'.format(self.wmin, w, self.wmax)\n",
    "                assert r >= 0 and r <= 1, 'r = {}'.format(r)\n",
    "                \n",
    "                decay = self.gamma ** c\n",
    "                self.n = decay * self.n + c\n",
    "                self.sumw = decay * self.sumw + c * w\n",
    "                self.sumwsq = decay * self.sumwsq + c * w**2\n",
    "                self.sumwr = decay * self.sumwr + c * w * r\n",
    "                self.sumwsqr = decay * self.sumwsqr + c * (w**2) * r\n",
    "                self.sumwsqrsq = decay * self.sumwsqrsq + c * (w**2) * (r**2)\n",
    "                    \n",
    "                self.duals = None\n",
    "                self.mleduals = None\n",
    "                \n",
    "            return self\n",
    "        \n",
    "        def recomputeduals(self):\n",
    "            from MLE.MLE import CrMinusTwo as CrMinusTwo\n",
    "            \n",
    "            self.duals = CrMinusTwo.intervalimpl(self.n, self.sumw, self.sumwsq, \n",
    "                                                 self.sumwr, self.sumwsqr, self.sumwsqrsq,\n",
    "                                                 self.wmin, self.wmax, self.alpha, raiseonerr=True)\n",
    "            \n",
    "        def recomputedualsmle(self):\n",
    "            from MLE.MLE import CrMinusTwo as CrMinusTwo\n",
    "            \n",
    "            self.mleduals = CrMinusTwo.estimateimpl(self.n, self.sumw, self.sumwsq, \n",
    "                                                    self.sumwr, self.sumwsqr, None, None,\n",
    "                                                    self.wmin, self.wmax, raiseonerr=True)\n",
    "        \n",
    "        def qlb(self, c, w, r):\n",
    "            if self.duals is None:\n",
    "                self.recomputeduals()\n",
    "                \n",
    "                assert self.duals is not None\n",
    "                \n",
    "            return self.duals[1][0]['qfunc'](c, w, r) if self.duals[1][0] is not None else 1\n",
    "        \n",
    "        def qub(self, c, w, r):\n",
    "            if self.duals is None:\n",
    "                self.recomputeduals()\n",
    "                \n",
    "                assert self.duals is not None\n",
    "                \n",
    "            return self.duals[1][1]['qfunc'](c, w, r) if self.duals[1][1] is not None else 1\n",
    "        \n",
    "        def qmle(self, c, w, r):\n",
    "            if self.mleduals is None:\n",
    "                self.recomputedualsmle()\n",
    "                \n",
    "                assert self.mleduals is not None\n",
    "                \n",
    "            return self.mleduals[1]['qfunc'](c, w, r) if self.mleduals[1] is not None else 1\n",
    "           \n",
    "    def flass():\n",
    "        from scipy.special import softmax\n",
    "        from sklearn.datasets import fetch_covtype\n",
    "        from sklearn.decomposition import PCA\n",
    "        from sklearn.linear_model import SGDClassifier\n",
    "        from sklearn.metrics import accuracy_score\n",
    "        from math import ceil, exp\n",
    "        import numpy as np\n",
    "        \n",
    "        cov = fetch_covtype()\n",
    "        cov.data = PCA(whiten=True).fit_transform(cov.data)\n",
    "        classes = np.unique(cov.target - 1)\n",
    "        ndata = len(cov.target)\n",
    "        order = np.random.RandomState(seed=42).permutation(ndata)\n",
    "        ntrain = ceil(0.2 * ndata)\n",
    "        Object = lambda **kwargs: type(\"Object\", (), kwargs)()\n",
    "        train = Object(data = cov.data[order[:ntrain]], target = cov.target[order[:ntrain]] - 1)\n",
    "        test = Object(data = cov.data[order[ntrain:]], target = cov.target[order[ntrain:]] - 1)\n",
    "        \n",
    "        subblocksize = 32\n",
    "        delay = 8\n",
    "        blocksize = delay * subblocksize\n",
    "        for lr, tau, what, gamma in ( (x, y, z, g) for x in np.logspace(-2.5, -2, 1) for y in (4, ) \n",
    "                                      for z in ('mle', 'ub', 'lb') for g in (1, 0.9999, 0.999, 0.99,) ):\n",
    "            print(\"*** lr = {} tau = {} what = {} gamma = {} ***\".format(lr, tau, what, gamma), flush=True)\n",
    "            print('{:8.8s}\\t{:8.8s}\\t{:10.10s}\\t{:8.8s}\\t{:10.10s}\\t{:8.8s}\\t{:10.10s}'.format(\n",
    "                'n', 'eff n', 'since last', 'emp loss', 'since last', 'log pv', 'since last')\n",
    "            )\n",
    "            \n",
    "            cls = SGDClassifier(loss='log', shuffle=False)\n",
    "            loss = OnlineDRO.EasyAcc()\n",
    "            sincelast = OnlineDRO.EasyAcc()\n",
    "            logpv = OnlineDRO.EasyAcc()\n",
    "            logpvsl = OnlineDRO.EasyAcc()\n",
    "            effn = OnlineDRO.EasyAcc()\n",
    "            effnsl = OnlineDRO.EasyAcc()\n",
    "            \n",
    "            loggerrand = np.random.RandomState(seed=2112)\n",
    "            logchoices = [None]*len(train.data)\n",
    "            pchoices = [None]*len(train.data)\n",
    "                        \n",
    "            ocrl = OnlineDRO.OnlineCressieReadLB(alpha=0.05, \n",
    "                                                 gamma=gamma,\n",
    "                                                 wmin=0,\n",
    "                                                 wmax=exp(tau) + len(classes) - 1\n",
    "                                                )\n",
    "            qfunc = ocrl.qmle if what == 'mle' else ocrl.qlb if what == 'lb' else ocrl.qub\n",
    " \n",
    "            for pno in range(1):\n",
    "                order = np.random.RandomState(seed=42+pno).permutation(len(train.data))\n",
    "                for n, ind in enumerate(zip(*(iter(order),)*blocksize)):\n",
    "                    v = np.array([ np.outer(t, np.append(t, [1])).ravel() for z in ind for t in ( train.data[z], ) ]) \n",
    "                    if n == 0 and pno == 0:\n",
    "                        for i, z in enumerate(ind):\n",
    "                            if logchoices[z] is None:\n",
    "                                choice = loggerrand.choice(a=classes, size=1)\n",
    "                                logchoices[z] = choice[0]\n",
    "                                pchoices[z] = 1.0 / len(classes)\n",
    "                    else:\n",
    "                        predlogp = cls.predict_proba(v)\n",
    "                        soft = softmax(tau * predlogp, axis=1)\n",
    "\n",
    "                        for i, z in enumerate(ind):\n",
    "                            if logchoices[z] is None:\n",
    "                                choice = loggerrand.choice(a=classes, p=soft[i,:], size=1)\n",
    "                                logchoices[z] = choice[0]\n",
    "                                pchoices[z] = soft[i, choice[0]]\n",
    "\n",
    "                        pred = cls.predict(v)\n",
    "                        actual = [ train.target[z] for z in ind ]\n",
    "                        for i, (p, a) in enumerate(zip(pred, actual)):\n",
    "                            loss += 0 if p == a else 1\n",
    "                            sincelast += 0 if p == a else 1\n",
    "                            logpv += soft[i, a]\n",
    "                            logpvsl += soft[i, a]\n",
    "                            \n",
    "                        if (n & (n - 1) == 0): # and n & 0xAAAAAAAA == 0):\n",
    "                            print('{:<8d}\\t{:<8.3f}\\t{:<10.3f}\\t{:<8.3f}\\t{:<10.3f}\\t{:<8.3f}\\t{:<10.3f}'.format(\n",
    "                                         loss.n, effn.mean(), effnsl.mean(), loss.mean(), sincelast.mean(), logpv.mean(), logpvsl.mean()),\n",
    "                                   flush=True)   \n",
    "\n",
    "                            sincelast = OnlineDRO.EasyAcc()\n",
    "                            logpvsl = OnlineDRO.EasyAcc()\n",
    "                            effnsl = OnlineDRO.EasyAcc()\n",
    "                            \n",
    "                    for i, z in enumerate(ind):\n",
    "                        r = 1 if logchoices[z] == train.target[z] else 0\n",
    "                        w = 1 / pchoices[z] \n",
    "                        ocrl.update(1, w, r) \n",
    "                                                   \n",
    "                    for d in range(delay):\n",
    "                        x = np.array([ v[i] \n",
    "                                       for i, z in enumerate(ind)\n",
    "                                       if (d-1)*subblocksize <= i and i < d*subblocksize\n",
    "                                       if logchoices[z] == train.target[z] ])\n",
    "                        y = np.array([ logchoices[z] \n",
    "                                       for i, z in enumerate(ind) \n",
    "                                       if (d-1)*subblocksize <= i and i < d*subblocksize\n",
    "                                       if logchoices[z] == train.target[z] ])\n",
    "                        \n",
    "                        if n == 0 and pno == 0:\n",
    "                            sampweight = np.array([ lr \n",
    "                                                    for i, z in enumerate(ind) \n",
    "                                                    if (d-1)*subblocksize <= i and i < d*subblocksize\n",
    "                                                    if logchoices[z] == train.target[z] ])\n",
    "                        else:\n",
    "                            sampweight = np.array([ lr * w * ocrl.n \n",
    "                                                       * max(0, qfunc(1, w, 1))\n",
    "                                                    for i, z in enumerate(ind) \n",
    "                                                    if (d-1)*subblocksize <= i and i < d*subblocksize\n",
    "                                                    if logchoices[z] == train.target[z]\n",
    "                                                    for w in (1 / pchoices[z],)\n",
    "                                                  ])\n",
    "                        \n",
    "                        effn += sampweight.sum() / (lr * subblocksize)\n",
    "                        effnsl += sampweight.sum() / (lr * subblocksize)\n",
    "                                             \n",
    "                        if np.any(x):\n",
    "                            cls.partial_fit(x, y, classes=classes, sample_weight=sampweight)\n",
    "\n",
    "                print('{:<8d}\\t{:<8.3f}\\t{:<10.3f}\\t{:<8.3f}\\t{:<10.3f}\\t{:<8.3f}\\t{:<10.3f}'.format(\n",
    "                             loss.n, effn.mean(), effnsl.mean(), loss.mean(), sincelast.mean(), logpv.mean(), logpvsl.mean()),\n",
    "                       flush=True)\n",
    "                \n",
    "#                 from pprint import pformat\n",
    "#                 print(pformat(ocrl.__dict__))\n",
    "\n",
    "                sincelast = OnlineDRO.EasyAcc()\n",
    "                logpvsl = OnlineDRO.EasyAcc()\n",
    "                effnsl = OnlineDRO.EasyAcc()\n",
    "\n",
    "                preds = cls.predict(np.array([np.outer(d, np.append(d, [1])).ravel() for d in test.data]))\n",
    "                ascores = []\n",
    "                for b in range(16):\n",
    "                    bootie = np.random.RandomState(90210+b).choice(len(test.target), replace=True, size=len(test.target))\n",
    "                    ascores.append(accuracy_score(y_true=test.target[bootie], y_pred=preds[bootie]))\n",
    "\n",
    "                print(\"test accuracy: {}\".format(np.quantile(ascores, [0.05, 0.5, 0.95])))\n",
    "    \n",
    "OnlineDRO.flass()      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Auto-Temperature + Delayed Batch Bound Online Learning\n",
    "\n",
    "Best regret achieved with lower bound optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "code_folding": [
     1,
     90,
     118,
     197,
     203,
     217,
     257
    ],
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.025 what = mle gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t1.025    [1.910     ]\t7.307    [6.613     ]\t0.432    [0.441     ]\t0.548    [0.552     ]\n",
      "1024    \t0.764    [0.502     ]\t6.931    [6.556     ]\t0.464    [0.496     ]\t0.525    [0.502     ]\n",
      "2048    \t0.763    [0.763     ]\t6.492    [6.053     ]\t0.474    [0.484     ]\t0.518    [0.512     ]\n",
      "4096    \t1.145    [1.527     ]\t6.465    [6.437     ]\t0.478    [0.481     ]\t0.518    [0.517     ]\n",
      "8192    \t1.044    [0.944     ]\t6.264    [6.063     ]\t0.469    [0.459     ]\t0.526    [0.533     ]\n",
      "16384   \t0.990    [0.936     ]\t6.390    [6.515     ]\t0.461    [0.453     ]\t0.526    [0.526     ]\n",
      "32768   \t0.779    [0.567     ]\t6.319    [6.249     ]\t0.455    [0.449     ]\t0.533    [0.540     ]\n",
      "65536   \t0.714    [0.649     ]\t6.579    [6.838     ]\t0.462    [0.469     ]\t0.524    [0.516     ]\n",
      "115712  \t0.825    [0.970     ]\t6.676    [6.803     ]\t0.467    [0.475     ]\t0.520    [0.515     ]\n",
      "test accuracy: [0.54671919 0.54791323 0.5490438 ]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.05 what = mle gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t1.025    [1.910     ]\t6.796    [5.593     ]\t0.432    [0.441     ]\t0.544    [0.543     ]\n",
      "1024    \t0.941    [0.857     ]\t6.406    [6.015     ]\t0.476    [0.520     ]\t0.496    [0.448     ]\n",
      "2048    \t0.884    [0.826     ]\t6.031    [5.657     ]\t0.478    [0.480     ]\t0.494    [0.492     ]\n",
      "4096    \t0.877    [0.870     ]\t5.730    [5.429     ]\t0.479    [0.480     ]\t0.498    [0.502     ]\n",
      "8192    \t0.820    [0.763     ]\t5.545    [5.359     ]\t0.465    [0.451     ]\t0.513    [0.528     ]\n",
      "16384   \t0.819    [0.818     ]\t5.510    [5.474     ]\t0.469    [0.474     ]\t0.507    [0.500     ]\n",
      "32768   \t0.724    [0.629     ]\t5.500    [5.490     ]\t0.462    [0.454     ]\t0.516    [0.525     ]\n",
      "65536   \t0.636    [0.548     ]\t5.364    [5.228     ]\t0.436    [0.409     ]\t0.542    [0.569     ]\n",
      "115712  \t0.717    [0.821     ]\t5.372    [5.382     ]\t0.440    [0.445     ]\t0.538    [0.533     ]\n",
      "test accuracy: [0.55968527 0.56140264 0.56268704]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.1 what = mle gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t1.025    [1.910     ]\t6.314    [4.628     ]\t0.432    [0.441     ]\t0.534    [0.523     ]\n",
      "1024    \t1.128    [1.232     ]\t5.873    [5.431     ]\t0.501    [0.570     ]\t0.453    [0.373     ]\n",
      "2048    \t0.988    [0.848     ]\t5.418    [4.964     ]\t0.510    [0.519     ]\t0.440    [0.427     ]\n",
      "4096    \t0.921    [0.854     ]\t4.986    [4.553     ]\t0.502    [0.494     ]\t0.456    [0.471     ]\n",
      "8192    \t0.937    [0.953     ]\t4.724    [4.463     ]\t0.484    [0.466     ]\t0.472    [0.489     ]\n",
      "16384   \t0.923    [0.909     ]\t4.642    [4.560     ]\t0.468    [0.452     ]\t0.487    [0.502     ]\n",
      "32768   \t0.865    [0.806     ]\t4.522    [4.403     ]\t0.453    [0.439     ]\t0.504    [0.520     ]\n",
      "65536   \t0.834    [0.803     ]\t4.497    [4.473     ]\t0.448    [0.443     ]\t0.510    [0.517     ]\n",
      "115712  \t0.802    [0.762     ]\t4.485    [4.468     ]\t0.421    [0.386     ]\t0.535    [0.566     ]\n",
      "test accuracy: [0.63792655 0.63885058 0.64035227]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.2 what = mle gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t5.345    [5.345     ]\t0.422    [0.422     ]\t0.503    [0.503     ]\n",
      "512     \t0.542    [0.943     ]\t4.516    [3.687     ]\t0.414    [0.406     ]\t0.504    [0.506     ]\n",
      "1024    \t0.765    [0.987     ]\t4.126    [3.735     ]\t0.430    [0.445     ]\t0.484    [0.464     ]\n",
      "2048    \t0.791    [0.817     ]\t3.904    [3.682     ]\t0.457    [0.484     ]\t0.459    [0.433     ]\n",
      "4096    \t0.814    [0.837     ]\t3.787    [3.671     ]\t0.465    [0.474     ]\t0.454    [0.450     ]\n",
      "8192    \t0.888    [0.962     ]\t3.745    [3.702     ]\t0.467    [0.469     ]\t0.451    [0.447     ]\n",
      "16384   \t0.930    [0.972     ]\t3.692    [3.638     ]\t0.455    [0.443     ]\t0.462    [0.473     ]\n",
      "32768   \t0.891    [0.853     ]\t3.634    [3.577     ]\t0.431    [0.407     ]\t0.483    [0.505     ]\n",
      "65536   \t0.884    [0.876     ]\t3.606    [3.578     ]\t0.419    [0.406     ]\t0.495    [0.506     ]\n",
      "115712  \t0.873    [0.858     ]\t3.607    [3.609     ]\t0.403    [0.383     ]\t0.508    [0.525     ]\n",
      "test accuracy: [0.65399497 0.65475496 0.65563974]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.025 what = mle gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t1.025    [1.910     ]\t7.307    [6.613     ]\t0.432    [0.441     ]\t0.548    [0.552     ]\n",
      "1024    \t0.764    [0.502     ]\t6.931    [6.556     ]\t0.464    [0.496     ]\t0.525    [0.502     ]\n",
      "2048    \t0.763    [0.763     ]\t6.492    [6.053     ]\t0.474    [0.484     ]\t0.518    [0.512     ]\n",
      "4096    \t1.145    [1.527     ]\t6.465    [6.437     ]\t0.478    [0.481     ]\t0.518    [0.517     ]\n",
      "8192    \t1.384    [1.623     ]\t6.316    [6.168     ]\t0.475    [0.473     ]\t0.518    [0.518     ]\n",
      "16384   \t1.177    [0.969     ]\t6.443    [6.570     ]\t0.474    [0.473     ]\t0.510    [0.503     ]\n",
      "32768   \t1.022    [0.868     ]\t6.472    [6.501     ]\t0.467    [0.460     ]\t0.514    [0.517     ]\n",
      "65536   \t0.965    [0.907     ]\t6.458    [6.443     ]\t0.467    [0.467     ]\t0.517    [0.519     ]\n",
      "115712  \t0.901    [0.819     ]\t6.348    [6.205     ]\t0.466    [0.464     ]\t0.518    [0.521     ]\n",
      "test accuracy: [0.53321956 0.5346992  0.53558666]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.05 what = mle gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t1.025    [1.910     ]\t6.796    [5.593     ]\t0.432    [0.441     ]\t0.544    [0.543     ]\n",
      "1024    \t0.941    [0.857     ]\t6.406    [6.015     ]\t0.476    [0.520     ]\t0.496    [0.448     ]\n",
      "2048    \t0.884    [0.826     ]\t6.031    [5.657     ]\t0.478    [0.480     ]\t0.494    [0.492     ]\n",
      "4096    \t0.877    [0.870     ]\t5.730    [5.429     ]\t0.479    [0.480     ]\t0.498    [0.502     ]\n",
      "8192    \t0.794    [0.712     ]\t5.545    [5.359     ]\t0.467    [0.455     ]\t0.511    [0.524     ]\n",
      "16384   \t0.822    [0.849     ]\t5.522    [5.500     ]\t0.482    [0.496     ]\t0.496    [0.481     ]\n",
      "32768   \t0.706    [0.590     ]\t5.423    [5.323     ]\t0.462    [0.443     ]\t0.514    [0.531     ]\n",
      "65536   \t0.714    [0.722     ]\t5.345    [5.267     ]\t0.452    [0.441     ]\t0.525    [0.537     ]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.1 what = mle gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/pmineiro/miniconda3/envs/elfcb/lib/python3.7/site-packages/sklearn/linear_model/base.py:309: RuntimeWarning: invalid value encountered in true_divide\n",
      "  prob /= prob.sum(axis=1).reshape((prob.shape[0], -1))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "512     \t1.025    [1.910     ]\t6.314    [4.628     ]\t0.432    [0.441     ]\t0.534    [0.523     ]\n",
      "1024    \t1.128    [1.232     ]\t5.873    [5.431     ]\t0.501    [0.570     ]\t0.453    [0.373     ]\n",
      "2048    \t0.988    [0.848     ]\t5.418    [4.964     ]\t0.510    [0.519     ]\t0.440    [0.427     ]\n",
      "4096    \t0.931    [0.873     ]\t4.986    [4.554     ]\t0.500    [0.490     ]\t0.456    [0.472     ]\n",
      "8192    \t0.907    [0.884     ]\t4.703    [4.421     ]\t0.485    [0.470     ]\t0.473    [0.490     ]\n",
      "16384   \t0.784    [0.660     ]\t4.638    [4.573     ]\t0.483    [0.481     ]\t0.476    [0.479     ]\n",
      "32768   \t0.766    [0.748     ]\t4.550    [4.462     ]\t0.447    [0.411     ]\t0.509    [0.543     ]\n",
      "65536   \t0.733    [0.701     ]\t4.513    [4.476     ]\t0.436    [0.426     ]\t0.521    [0.532     ]\n",
      "115712  \t0.744    [0.758     ]\t4.507    [4.499     ]\t0.426    [0.411     ]\t0.531    [0.544     ]\n",
      "test accuracy: [0.63185577 0.63279218 0.63372536]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.2 what = mle gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t5.345    [5.345     ]\t0.422    [0.422     ]\t0.503    [0.503     ]\n",
      "512     \t0.542    [0.943     ]\t4.516    [3.687     ]\t0.414    [0.406     ]\t0.504    [0.506     ]\n",
      "1024    \t0.765    [0.987     ]\t4.126    [3.735     ]\t0.430    [0.445     ]\t0.484    [0.464     ]\n",
      "2048    \t0.791    [0.817     ]\t3.904    [3.682     ]\t0.457    [0.484     ]\t0.459    [0.433     ]\n",
      "4096    \t0.814    [0.837     ]\t3.787    [3.671     ]\t0.465    [0.474     ]\t0.454    [0.450     ]\n",
      "8192    \t0.860    [0.906     ]\t3.729    [3.670     ]\t0.462    [0.458     ]\t0.456    [0.458     ]\n",
      "16384   \t0.916    [0.972     ]\t3.672    [3.615     ]\t0.452    [0.443     ]\t0.464    [0.472     ]\n",
      "32768   \t0.836    [0.755     ]\t3.646    [3.620     ]\t0.442    [0.432     ]\t0.472    [0.480     ]\n",
      "65536   \t0.825    [0.815     ]\t3.627    [3.608     ]\t0.410    [0.379     ]\t0.499    [0.526     ]\n",
      "115712  \t0.823    [0.820     ]\t3.633    [3.640     ]\t0.392    [0.369     ]\t0.515    [0.535     ]\n",
      "test accuracy: [0.66067352 0.66141792 0.66278783]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.025 what = ub gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t1.053    [1.966     ]\t7.305    [6.610     ]\t0.432    [0.441     ]\t0.548    [0.552     ]\n",
      "1024    \t0.785    [0.516     ]\t6.928    [6.552     ]\t0.465    [0.498     ]\t0.525    [0.501     ]\n",
      "2048    \t0.788    [0.792     ]\t6.470    [6.012     ]\t0.474    [0.483     ]\t0.518    [0.512     ]\n",
      "4096    \t1.099    [1.409     ]\t6.427    [6.385     ]\t0.473    [0.471     ]\t0.521    [0.524     ]\n",
      "8192    \t1.160    [1.222     ]\t6.597    [6.767     ]\t0.478    [0.482     ]\t0.507    [0.493     ]\n",
      "16384   \t1.018    [0.875     ]\t6.502    [6.406     ]\t0.466    [0.455     ]\t0.517    [0.528     ]\n",
      "32768   \t0.872    [0.726     ]\t6.507    [6.512     ]\t0.473    [0.479     ]\t0.509    [0.500     ]\n",
      "65536   \t1.018    [1.164     ]\t6.502    [6.498     ]\t0.482    [0.492     ]\t0.497    [0.486     ]\n",
      "115712  \t1.918    [3.089     ]\t6.735    [7.040     ]\t0.496    [0.513     ]\t0.484    [0.466     ]\n",
      "test accuracy: [0.45444419 0.45516546 0.45614919]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.05 what = ub gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t1.053    [1.966     ]\t6.795    [5.590     ]\t0.432    [0.441     ]\t0.544    [0.543     ]\n",
      "1024    \t0.896    [0.740     ]\t6.404    [6.012     ]\t0.484    [0.537     ]\t0.490    [0.437     ]\n",
      "2048    \t0.868    [0.839     ]\t6.014    [5.625     ]\t0.496    [0.508     ]\t0.482    [0.474     ]\n",
      "4096    \t0.861    [0.853     ]\t5.850    [5.685     ]\t0.476    [0.456     ]\t0.501    [0.519     ]\n",
      "8192    \t0.820    [0.779     ]\t5.793    [5.736     ]\t0.473    [0.469     ]\t0.501    [0.502     ]\n",
      "16384   \t0.884    [0.949     ]\t5.586    [5.379     ]\t0.475    [0.478     ]\t0.500    [0.500     ]\n",
      "32768   \t0.962    [1.040     ]\t5.571    [5.557     ]\t0.480    [0.484     ]\t0.497    [0.494     ]\n",
      "65536   \t1.157    [1.351     ]\t5.508    [5.444     ]\t0.478    [0.475     ]\t0.499    [0.502     ]\n",
      "115712  \t1.129    [1.093     ]\t5.597    [5.714     ]\t0.500    [0.530     ]\t0.474    [0.441     ]\n",
      "test accuracy: [0.472635   0.47394844 0.47486871]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.1 what = ub gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t1.053    [1.966     ]\t6.313    [4.626     ]\t0.432    [0.441     ]\t0.534    [0.523     ]\n",
      "1024    \t1.159    [1.265     ]\t5.870    [5.428     ]\t0.501    [0.570     ]\t0.454    [0.373     ]\n",
      "2048    \t1.017    [0.874     ]\t5.421    [4.971     ]\t0.510    [0.520     ]\t0.440    [0.427     ]\n",
      "4096    \t0.973    [0.929     ]\t4.980    [4.540     ]\t0.503    [0.496     ]\t0.454    [0.468     ]\n",
      "8192    \t0.948    [0.923     ]\t4.746    [4.512     ]\t0.482    [0.461     ]\t0.473    [0.492     ]\n",
      "16384   \t0.916    [0.885     ]\t4.605    [4.465     ]\t0.457    [0.432     ]\t0.496    [0.519     ]\n",
      "32768   \t0.994    [1.071     ]\t4.599    [4.593     ]\t0.444    [0.431     ]\t0.508    [0.520     ]\n",
      "65536   \t0.990    [0.986     ]\t4.648    [4.697     ]\t0.448    [0.451     ]\t0.506    [0.504     ]\n",
      "115712  \t0.967    [0.937     ]\t4.579    [4.489     ]\t0.429    [0.404     ]\t0.525    [0.551     ]\n",
      "test accuracy: [0.58680017 0.58770699 0.58861274]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.2 what = ub gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t5.345    [5.345     ]\t0.422    [0.422     ]\t0.503    [0.503     ]\n",
      "512     \t0.555    [0.970     ]\t4.513    [3.681     ]\t0.414    [0.406     ]\t0.504    [0.506     ]\n",
      "1024    \t0.783    [1.010     ]\t4.120    [3.726     ]\t0.428    [0.441     ]\t0.485    [0.465     ]\n",
      "2048    \t0.771    [0.760     ]\t3.893    [3.666     ]\t0.454    [0.480     ]\t0.461    [0.437     ]\n",
      "4096    \t0.797    [0.824     ]\t3.772    [3.651     ]\t0.469    [0.483     ]\t0.452    [0.442     ]\n",
      "8192    \t0.867    [0.937     ]\t3.748    [3.725     ]\t0.466    [0.464     ]\t0.451    [0.450     ]\n",
      "16384   \t0.933    [0.998     ]\t3.696    [3.645     ]\t0.456    [0.445     ]\t0.460    [0.469     ]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.025 what = ub gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t1.058    [1.975     ]\t7.305    [6.609     ]\t0.432    [0.441     ]\t0.548    [0.552     ]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/pmineiro/miniconda3/envs/elfcb/lib/python3.7/site-packages/sklearn/linear_model/base.py:309: RuntimeWarning: invalid value encountered in true_divide\n",
      "  prob /= prob.sum(axis=1).reshape((prob.shape[0], -1))\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1024    \t0.789    [0.521     ]\t6.926    [6.548     ]\t0.465    [0.498     ]\t0.524    [0.501     ]\n",
      "2048    \t0.798    [0.808     ]\t6.475    [6.024     ]\t0.474    [0.482     ]\t0.519    [0.513     ]\n",
      "4096    \t1.123    [1.448     ]\t6.396    [6.317     ]\t0.475    [0.477     ]\t0.520    [0.520     ]\n",
      "8192    \t1.118    [1.112     ]\t6.549    [6.701     ]\t0.488    [0.500     ]\t0.499    [0.478     ]\n",
      "16384   \t1.268    [1.417     ]\t6.615    [6.682     ]\t0.495    [0.503     ]\t0.485    [0.472     ]\n",
      "32768   \t2.960    [4.652     ]\t6.644    [6.673     ]\t0.507    [0.519     ]\t0.471    [0.457     ]\n",
      "65536   \t8.751    [14.541    ]\t7.068    [7.493     ]\t0.525    [0.543     ]\t0.444    [0.417     ]\n",
      "115712  \t10.477   [12.719    ]\t7.204    [7.382     ]\t0.517    [0.506     ]\t0.459    [0.478     ]\n",
      "test accuracy: [0.50706957 0.50841421 0.50952918]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.05 what = ub gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t1.058    [1.975     ]\t6.795    [5.590     ]\t0.432    [0.441     ]\t0.544    [0.543     ]\n",
      "1024    \t0.902    [0.747     ]\t6.403    [6.010     ]\t0.485    [0.539     ]\t0.490    [0.437     ]\n",
      "2048    \t0.880    [0.859     ]\t6.011    [5.619     ]\t0.498    [0.510     ]\t0.482    [0.474     ]\n",
      "4096    \t0.945    [1.010     ]\t5.774    [5.538     ]\t0.478    [0.459     ]\t0.501    [0.521     ]\n",
      "8192    \t1.347    [1.748     ]\t5.475    [5.175     ]\t0.498    [0.518     ]\t0.482    [0.463     ]\n",
      "16384   \t2.018    [2.690     ]\t5.679    [5.883     ]\t0.533    [0.569     ]\t0.441    [0.401     ]\n",
      "32768   \t2.267    [2.516     ]\t5.817    [5.955     ]\t0.562    [0.590     ]\t0.413    [0.385     ]\n",
      "65536   \t5.685    [9.103     ]\t6.268    [6.718     ]\t0.560    [0.558     ]\t0.400    [0.386     ]\n",
      "115712  \t4.417    [2.770     ]\t6.269    [6.272     ]\t0.557    [0.553     ]\t0.400    [0.399     ]\n",
      "test accuracy: [0.51070547 0.51174353 0.512766  ]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.1 what = ub gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t1.058    [1.975     ]\t6.313    [4.626     ]\t0.432    [0.441     ]\t0.534    [0.523     ]\n",
      "1024    \t1.167    [1.276     ]\t5.870    [5.426     ]\t0.501    [0.570     ]\t0.454    [0.373     ]\n",
      "2048    \t1.055    [0.943     ]\t5.391    [4.913     ]\t0.509    [0.517     ]\t0.441    [0.429     ]\n",
      "4096    \t0.975    [0.895     ]\t4.947    [4.502     ]\t0.504    [0.500     ]\t0.453    [0.464     ]\n",
      "8192    \t1.065    [1.154     ]\t4.684    [4.421     ]\t0.488    [0.472     ]\t0.470    [0.488     ]\n",
      "16384   \t1.152    [1.238     ]\t4.717    [4.750     ]\t0.480    [0.471     ]\t0.473    [0.477     ]\n",
      "32768   \t1.579    [2.006     ]\t4.839    [4.961     ]\t0.478    [0.476     ]\t0.468    [0.463     ]\n",
      "65536   \t1.523    [1.467     ]\t4.775    [4.710     ]\t0.473    [0.468     ]\t0.477    [0.485     ]\n",
      "115712  \t1.517    [1.510     ]\t4.671    [4.535     ]\t0.455    [0.432     ]\t0.495    [0.520     ]\n",
      "test accuracy: [0.60193004 0.60324133 0.60412718]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.2 what = ub gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t5.345    [5.345     ]\t0.422    [0.422     ]\t0.503    [0.503     ]\n",
      "512     \t0.558    [0.974     ]\t4.513    [3.681     ]\t0.414    [0.406     ]\t0.504    [0.506     ]\n",
      "1024    \t0.787    [1.017     ]\t4.119    [3.724     ]\t0.429    [0.443     ]\t0.484    [0.465     ]\n",
      "2048    \t0.773    [0.759     ]\t3.888    [3.657     ]\t0.456    [0.482     ]\t0.460    [0.436     ]\n",
      "4096    \t0.854    [0.934     ]\t3.760    [3.633     ]\t0.477    [0.498     ]\t0.447    [0.434     ]\n",
      "8192    \t1.055    [1.255     ]\t3.723    [3.685     ]\t0.469    [0.461     ]\t0.450    [0.454     ]\n",
      "16384   \t1.066    [1.078     ]\t3.722    [3.722     ]\t0.472    [0.475     ]\t0.445    [0.440     ]\n",
      "32768   \t1.116    [1.166     ]\t3.666    [3.609     ]\t0.458    [0.444     ]\t0.458    [0.470     ]\n",
      "65536   \t1.145    [1.175     ]\t3.623    [3.581     ]\t0.437    [0.416     ]\t0.476    [0.495     ]\n",
      "115712  \t1.115    [1.075     ]\t3.603    [3.575     ]\t0.415    [0.386     ]\t0.495    [0.520     ]\n",
      "test accuracy: [0.63125391 0.63224895 0.63314501]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.025 what = lb gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t0.315    [0.490     ]\t8.000    [8.000     ]\t0.443    [0.465     ]\t0.530    [0.516     ]\n",
      "1024    \t0.391    [0.466     ]\t7.688    [7.376     ]\t0.431    [0.418     ]\t0.551    [0.572     ]\n",
      "2048    \t0.465    [0.540     ]\t7.118    [6.548     ]\t0.419    [0.408     ]\t0.566    [0.581     ]\n",
      "4096    \t0.497    [0.529     ]\t6.763    [6.409     ]\t0.409    [0.399     ]\t0.578    [0.591     ]\n",
      "8192    \t0.521    [0.546     ]\t6.610    [6.457     ]\t0.391    [0.373     ]\t0.597    [0.617     ]\n",
      "16384   \t0.545    [0.569     ]\t6.505    [6.400     ]\t0.373    [0.355     ]\t0.615    [0.632     ]\n",
      "32768   \t0.563    [0.581     ]\t6.426    [6.348     ]\t0.369    [0.365     ]\t0.619    [0.623     ]\n",
      "65536   \t0.521    [0.478     ]\t6.328    [6.229     ]\t0.446    [0.523     ]\t0.546    [0.472     ]\n",
      "115712  \t0.490    [0.450     ]\t6.269    [6.191     ]\t0.467    [0.495     ]\t0.525    [0.498     ]\n",
      "test accuracy: [0.51810206 0.51917239 0.52021529]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.05 what = lb gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t0.315    [0.490     ]\t7.337    [6.675     ]\t0.443    [0.465     ]\t0.526    [0.507     ]\n",
      "1024    \t0.391    [0.466     ]\t6.647    [5.956     ]\t0.438    [0.432     ]\t0.536    [0.547     ]\n",
      "2048    \t0.458    [0.526     ]\t6.102    [5.557     ]\t0.434    [0.430     ]\t0.545    [0.553     ]\n",
      "4096    \t0.485    [0.512     ]\t5.803    [5.503     ]\t0.432    [0.431     ]\t0.546    [0.548     ]\n",
      "8192    \t0.500    [0.514     ]\t5.610    [5.418     ]\t0.420    [0.408     ]\t0.558    [0.569     ]\n",
      "16384   \t0.529    [0.558     ]\t5.571    [5.532     ]\t0.400    [0.379     ]\t0.578    [0.599     ]\n",
      "32768   \t0.557    [0.584     ]\t5.545    [5.520     ]\t0.380    [0.359     ]\t0.598    [0.617     ]\n",
      "65536   \t0.572    [0.587     ]\t5.506    [5.468     ]\t0.361    [0.343     ]\t0.615    [0.632     ]\n",
      "115712  \t0.588    [0.609     ]\t5.523    [5.544     ]\t0.345    [0.324     ]\t0.630    [0.649     ]\n",
      "test accuracy: [0.68372331 0.68474255 0.68549985]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.1 what = lb gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t0.315    [0.490     ]\t6.642    [5.284     ]\t0.443    [0.465     ]\t0.516    [0.487     ]\n",
      "1024    \t0.383    [0.451     ]\t5.742    [4.841     ]\t0.438    [0.434     ]\t0.518    [0.520     ]\n",
      "2048    \t0.456    [0.528     ]\t5.216    [4.691     ]\t0.423    [0.408     ]\t0.533    [0.548     ]\n",
      "4096    \t0.495    [0.534     ]\t4.897    [4.579     ]\t0.414    [0.404     ]\t0.542    [0.552     ]\n",
      "8192    \t0.524    [0.554     ]\t4.698    [4.499     ]\t0.390    [0.365     ]\t0.564    [0.586     ]\n",
      "16384   \t0.550    [0.575     ]\t4.627    [4.557     ]\t0.372    [0.354     ]\t0.580    [0.596     ]\n",
      "32768   \t0.573    [0.596     ]\t4.601    [4.574     ]\t0.355    [0.338     ]\t0.595    [0.610     ]\n",
      "65536   \t0.586    [0.599     ]\t4.609    [4.618     ]\t0.353    [0.350     ]\t0.598    [0.602     ]\n",
      "115712  \t0.592    [0.600     ]\t4.602    [4.593     ]\t0.348    [0.341     ]\t0.603    [0.609     ]\n",
      "test accuracy: [0.66968529 0.67036245 0.6718437 ]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.2 what = lb gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t5.345    [5.345     ]\t0.422    [0.422     ]\t0.503    [0.503     ]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "512     \t0.320    [0.499     ]\t4.681    [4.016     ]\t0.449    [0.477     ]\t0.477    [0.451     ]\n",
      "1024    \t0.398    [0.477     ]\t4.215    [3.750     ]\t0.440    [0.432     ]\t0.481    [0.485     ]\n",
      "2048    \t0.465    [0.532     ]\t3.979    [3.743     ]\t0.417    [0.393     ]\t0.500    [0.520     ]\n",
      "4096    \t0.512    [0.559     ]\t3.798    [3.617     ]\t0.398    [0.379     ]\t0.512    [0.524     ]\n",
      "8192    \t0.533    [0.553     ]\t3.697    [3.596     ]\t0.388    [0.378     ]\t0.520    [0.528     ]\n",
      "16384   \t0.564    [0.595     ]\t3.669    [3.640     ]\t0.369    [0.350     ]\t0.535    [0.549     ]\n",
      "32768   \t0.591    [0.617     ]\t3.680    [3.692     ]\t0.355    [0.341     ]\t0.545    [0.556     ]\n",
      "65536   \t0.618    [0.645     ]\t3.718    [3.756     ]\t0.342    [0.328     ]\t0.555    [0.565     ]\n",
      "115712  \t0.645    [0.680     ]\t3.780    [3.861     ]\t0.334    [0.323     ]\t0.561    [0.568     ]\n",
      "test accuracy: [0.69895161 0.70030916 0.70169414]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.025 what = lb gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t0.313    [0.486     ]\t8.000    [8.000     ]\t0.443    [0.465     ]\t0.531    [0.517     ]\n",
      "1024    \t0.384    [0.455     ]\t7.711    [7.422     ]\t0.434    [0.424     ]\t0.549    [0.567     ]\n",
      "2048    \t0.450    [0.515     ]\t7.116    [6.520     ]\t0.425    [0.417     ]\t0.561    [0.573     ]\n",
      "4096    \t0.471    [0.492     ]\t6.682    [6.248     ]\t0.421    [0.417     ]\t0.567    [0.573     ]\n",
      "8192    \t0.493    [0.516     ]\t6.415    [6.148     ]\t0.402    [0.383     ]\t0.587    [0.606     ]\n",
      "16384   \t0.511    [0.528     ]\t6.329    [6.243     ]\t0.388    [0.374     ]\t0.601    [0.615     ]\n",
      "32768   \t0.532    [0.553     ]\t6.245    [6.160     ]\t0.421    [0.454     ]\t0.568    [0.536     ]\n",
      "65536   \t0.496    [0.459     ]\t6.379    [6.513     ]\t0.453    [0.484     ]\t0.536    [0.504     ]\n",
      "115712  \t0.487    [0.475     ]\t6.365    [6.346     ]\t0.447    [0.439     ]\t0.542    [0.550     ]\n",
      "test accuracy: [0.58256456 0.58331809 0.58431313]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.05 what = lb gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t0.313    [0.486     ]\t7.349    [6.697     ]\t0.443    [0.465     ]\t0.526    [0.508     ]\n",
      "1024    \t0.388    [0.463     ]\t6.663    [5.978     ]\t0.438    [0.434     ]\t0.537    [0.548     ]\n",
      "2048    \t0.447    [0.505     ]\t6.154    [5.644     ]\t0.433    [0.427     ]\t0.546    [0.555     ]\n",
      "4096    \t0.464    [0.482     ]\t5.807    [5.459     ]\t0.432    [0.431     ]\t0.547    [0.549     ]\n",
      "8192    \t0.472    [0.479     ]\t5.593    [5.379     ]\t0.429    [0.426     ]\t0.550    [0.554     ]\n",
      "16384   \t0.484    [0.497     ]\t5.524    [5.455     ]\t0.426    [0.422     ]\t0.553    [0.555     ]\n",
      "32768   \t0.499    [0.514     ]\t5.503    [5.483     ]\t0.409    [0.392     ]\t0.570    [0.586     ]\n",
      "65536   \t0.514    [0.529     ]\t5.454    [5.405     ]\t0.397    [0.385     ]\t0.581    [0.591     ]\n",
      "115712  \t0.509    [0.503     ]\t5.444    [5.431     ]\t0.404    [0.414     ]\t0.573    [0.563     ]\n",
      "test accuracy: [0.58322827 0.5844207  0.58574113]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.1 what = lb gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t8.000    [8.000     ]\t0.422    [0.422     ]\t0.544    [0.544     ]\n",
      "512     \t0.313    [0.486     ]\t6.648    [5.297     ]\t0.443    [0.465     ]\t0.516    [0.488     ]\n",
      "1024    \t0.378    [0.444     ]\t5.751    [4.853     ]\t0.438    [0.432     ]\t0.519    [0.522     ]\n",
      "2048    \t0.438    [0.498     ]\t5.219    [4.688     ]\t0.428    [0.418     ]\t0.530    [0.541     ]\n",
      "4096    \t0.468    [0.497     ]\t4.897    [4.574     ]\t0.426    [0.425     ]\t0.532    [0.534     ]\n",
      "8192    \t0.477    [0.486     ]\t4.720    [4.543     ]\t0.424    [0.422     ]\t0.534    [0.535     ]\n",
      "16384   \t0.506    [0.535     ]\t4.636    [4.553     ]\t0.394    [0.364     ]\t0.561    [0.589     ]\n",
      "32768   \t0.522    [0.539     ]\t4.584    [4.532     ]\t0.375    [0.357     ]\t0.577    [0.594     ]\n",
      "65536   \t0.531    [0.539     ]\t4.597    [4.611     ]\t0.372    [0.369     ]\t0.582    [0.586     ]\n",
      "115712  \t0.532    [0.535     ]\t4.572    [4.539     ]\t0.371    [0.370     ]\t0.582    [0.583     ]\n",
      "test accuracy: [0.63460637 0.63570413 0.63726391]\n",
      "*** lr = 0.0031622776601683794 taumax = 8 target = 0.2 what = lb gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.141    [0.141     ]\t5.345    [5.345     ]\t0.422    [0.422     ]\t0.503    [0.503     ]\n",
      "512     \t0.314    [0.487     ]\t4.696    [4.047     ]\t0.449    [0.477     ]\t0.478    [0.454     ]\n",
      "1024    \t0.386    [0.458     ]\t4.229    [3.762     ]\t0.438    [0.428     ]\t0.481    [0.484     ]\n",
      "2048    \t0.437    [0.488     ]\t3.973    [3.717     ]\t0.417    [0.396     ]\t0.496    [0.512     ]\n",
      "4096    \t0.465    [0.493     ]\t3.799    [3.626     ]\t0.418    [0.418     ]\t0.496    [0.495     ]\n",
      "8192    \t0.491    [0.517     ]\t3.708    [3.616     ]\t0.401    [0.385     ]\t0.509    [0.522     ]\n",
      "16384   \t0.526    [0.561     ]\t3.674    [3.640     ]\t0.370    [0.338     ]\t0.534    [0.560     ]\n",
      "32768   \t0.546    [0.566     ]\t3.661    [3.648     ]\t0.349    [0.328     ]\t0.550    [0.566     ]\n",
      "65536   \t0.554    [0.562     ]\t3.647    [3.634     ]\t0.342    [0.336     ]\t0.554    [0.559     ]\n",
      "115712  \t0.559    [0.566     ]\t3.637    [3.624     ]\t0.337    [0.331     ]\t0.558    [0.562     ]\n",
      "test accuracy: [0.67261499 0.67346588 0.67449963]\n"
     ]
    }
   ],
   "source": [
    "class OnlineDRO:\n",
    "    class EasyAcc:\n",
    "        def __init__(self):\n",
    "            self.n = 0\n",
    "            self.sum = 0\n",
    "            \n",
    "        def __iadd__(self, other):\n",
    "            self.n += 1\n",
    "            self.sum += other\n",
    "            return self\n",
    "            \n",
    "        def mean(self):\n",
    "            return self.sum / max(self.n, 1)\n",
    "\n",
    "    class OnlineCressieReadLB:\n",
    "        from math import inf\n",
    "        \n",
    "        def __init__(self, alpha, gamma=1, wmin=0, wmax=inf):\n",
    "            import numpy as np\n",
    "            \n",
    "            self.alpha = alpha\n",
    "            self.gamma = gamma\n",
    "            self.n = 0\n",
    "            self.sumw = 0\n",
    "            self.sumwsq = 0\n",
    "            self.sumwr = 0\n",
    "            self.sumwsqr = 0\n",
    "            self.sumwsqrsq = 0\n",
    "            self.wmin = wmin\n",
    "            self.wmax = wmax\n",
    "            \n",
    "            self.duals = None\n",
    "            self.mleduals = None\n",
    "            \n",
    "        def update(self, c, w, r):\n",
    "            if c > 0:\n",
    "                assert w + 1e-6 >= self.wmin and w <= self.wmax + 1e-6, 'w = {} < {} < {}'.format(self.wmin, w, self.wmax)\n",
    "                assert r >= 0 and r <= 1, 'r = {}'.format(r)\n",
    "                \n",
    "                decay = self.gamma ** c\n",
    "                self.n = decay * self.n + c\n",
    "                self.sumw = decay * self.sumw + c * w\n",
    "                self.sumwsq = decay * self.sumwsq + c * w**2\n",
    "                self.sumwr = decay * self.sumwr + c * w * r\n",
    "                self.sumwsqr = decay * self.sumwsqr + c * (w**2) * r\n",
    "                self.sumwsqrsq = decay * self.sumwsqrsq + c * (w**2) * (r**2)\n",
    "                    \n",
    "                self.duals = None\n",
    "                self.mleduals = None\n",
    "                \n",
    "            return self\n",
    "        \n",
    "        def recomputeduals(self):\n",
    "            from MLE.MLE import CrMinusTwo as CrMinusTwo\n",
    "            \n",
    "            self.duals = CrMinusTwo.intervalimpl(self.n, self.sumw, self.sumwsq, \n",
    "                                                 self.sumwr, self.sumwsqr, self.sumwsqrsq,\n",
    "                                                 self.wmin, self.wmax, self.alpha, raiseonerr=True)\n",
    "            \n",
    "        def recomputedualsmle(self):\n",
    "            from MLE.MLE import CrMinusTwo as CrMinusTwo\n",
    "            \n",
    "            self.mleduals = CrMinusTwo.estimateimpl(self.n, self.sumw, self.sumwsq, \n",
    "                                                    self.sumwr, self.sumwsqr, None, None,\n",
    "                                                    self.wmin, self.wmax, raiseonerr=True)\n",
    "        \n",
    "        def qlb(self, c, w, r):\n",
    "            if self.duals is None:\n",
    "                self.recomputeduals()\n",
    "                \n",
    "                assert self.duals is not None\n",
    "                \n",
    "            return self.duals[1][0]['qfunc'](c, w, r) if self.duals[1][0] is not None else 1\n",
    "        \n",
    "        def qub(self, c, w, r):\n",
    "            if self.duals is None:\n",
    "                self.recomputeduals()\n",
    "                \n",
    "                assert self.duals is not None\n",
    "                \n",
    "            return self.duals[1][1]['qfunc'](c, w, r) if self.duals[1][1] is not None else 1\n",
    "        \n",
    "        def qmle(self, c, w, r):\n",
    "            if self.mleduals is None:\n",
    "                self.recomputedualsmle()\n",
    "                \n",
    "                assert self.mleduals is not None\n",
    "                \n",
    "            return self.mleduals[1]['qfunc'](c, w, r) if self.mleduals[1] is not None else 1\n",
    "           \n",
    "    def autotune(pre, target, taumax):\n",
    "        from scipy.optimize import root_scalar\n",
    "        \n",
    "        def f(tau):\n",
    "            from scipy.special import softmax\n",
    "            import numpy as np\n",
    "\n",
    "            soft = softmax(tau * pre, axis=1)\n",
    "            minsoft = np.min(soft, axis=1)\n",
    "            \n",
    "            return np.mean(minsoft) - target\n",
    "            \n",
    "        fmax = f(taumax)\n",
    "        \n",
    "        if fmax >= 0:\n",
    "            return taumax\n",
    "        \n",
    "        taumin = 0\n",
    "        fmin = f(taumin)\n",
    "        \n",
    "        assert fmin > 0, { 'fmin': fmin, 'target': target }\n",
    "        \n",
    "        root = root_scalar(f, bracket=(taumin, taumax))\n",
    "        \n",
    "        assert root.converged, root\n",
    "                \n",
    "        return root.root\n",
    "            \n",
    "    def flass():\n",
    "        from scipy.special import softmax\n",
    "        from sklearn.datasets import fetch_covtype\n",
    "        from sklearn.decomposition import PCA\n",
    "        from sklearn.linear_model import SGDClassifier\n",
    "        from sklearn.metrics import accuracy_score\n",
    "        from math import ceil, exp\n",
    "        import numpy as np\n",
    "        \n",
    "        cov = fetch_covtype()\n",
    "        cov.data = PCA(whiten=True).fit_transform(cov.data)\n",
    "        classes = np.unique(cov.target - 1)\n",
    "        ndata = len(cov.target)\n",
    "        order = np.random.RandomState(seed=42).permutation(ndata)\n",
    "        ntrain = ceil(0.2 * ndata)\n",
    "        Object = lambda **kwargs: type(\"Object\", (), kwargs)()\n",
    "        train = Object(data = cov.data[order[:ntrain]], target = cov.target[order[:ntrain]] - 1)\n",
    "        test = Object(data = cov.data[order[ntrain:]], target = cov.target[order[ntrain:]] - 1)\n",
    "        \n",
    "        subblocksize = 32\n",
    "        delay = 8\n",
    "        blocksize = delay * subblocksize\n",
    "        for lr, taumax, target, what, gamma in ( (x, 8, y, z, g) for x in np.logspace(-2.5, -2, 1) \n",
    "                                         for z in ('mle', 'ub', 'lb') for g in (0.9999, 0.999,) for y in (0.025, 0.05, 0.1, 0.2, ) ):\n",
    "            print(\"*** lr = {} taumax = {} target = {} what = {} gamma = {} ***\".format(lr, taumax, target, what, gamma), flush=True)\n",
    "            print('{:8.8s}\\t{:8.8s} [{:10.10s}]\\t{:8.8s} [{:10.10s}]\\t{:8.8s} [{:10.10s}]\\t{:8.8s} [{:10.10s}]'.format(\n",
    "                'n', 'eff n', 'since last', 'av tau', 'since last', 'emp loss', 'since last', 'log pv', 'since last')\n",
    "            )\n",
    "            \n",
    "            try:\n",
    "                target /= len(classes)\n",
    "\n",
    "                cls = SGDClassifier(loss='log', shuffle=False)\n",
    "                loss = OnlineDRO.EasyAcc()\n",
    "                sincelast = OnlineDRO.EasyAcc()\n",
    "                logpv = OnlineDRO.EasyAcc()\n",
    "                logpvsl = OnlineDRO.EasyAcc()\n",
    "                effn = OnlineDRO.EasyAcc()\n",
    "                effnsl = OnlineDRO.EasyAcc()\n",
    "                avtau = OnlineDRO.EasyAcc()\n",
    "                avtausl = OnlineDRO.EasyAcc()\n",
    "\n",
    "                loggerrand = np.random.RandomState(seed=2112)\n",
    "                logchoices = [None]*len(train.data)\n",
    "                pchoices = [None]*len(train.data)\n",
    "\n",
    "                ocrl = OnlineDRO.OnlineCressieReadLB(alpha=0.05, \n",
    "                                                     gamma=gamma,\n",
    "                                                     wmin=0,\n",
    "#                                                      wmax=exp(taumax) + len(classes) - 1\n",
    "                                                    )\n",
    "                qfunc = ocrl.qmle if what == 'mle' else ocrl.qlb if what == 'lb' else ocrl.qub\n",
    "\n",
    "                for pno in range(1):\n",
    "                    order = np.random.RandomState(seed=42+pno).permutation(len(train.data))\n",
    "                    for n, ind in enumerate(zip(*(iter(order),)*blocksize)):\n",
    "                        v = np.array([ np.outer(t, np.append(t, [1])).ravel() for z in ind for t in ( train.data[z], ) ]) \n",
    "                        if n == 0 and pno == 0:\n",
    "                            pred = np.zeros(blocksize)\n",
    "                            for i, z in enumerate(ind):\n",
    "                                if logchoices[z] is None:\n",
    "                                    choice = loggerrand.choice(a=classes, size=1)\n",
    "                                    logchoices[z] = choice[0]\n",
    "                                    pchoices[z] = 1.0 / len(classes)\n",
    "                        else:\n",
    "                            predlogp = cls.predict_proba(v)\n",
    "                            tau = OnlineDRO.autotune(predlogp, target, taumax)\n",
    "                            avtau += tau\n",
    "                            avtausl += tau\n",
    "                            soft = softmax(tau * predlogp, axis=1)\n",
    "\n",
    "                            for i, z in enumerate(ind):\n",
    "                                if logchoices[z] is None:\n",
    "                                    choice = loggerrand.choice(a=classes, p=soft[i,:], size=1)\n",
    "                                    logchoices[z] = choice[0]\n",
    "                                    pchoices[z] = soft[i, choice[0]]\n",
    "\n",
    "                            pred = cls.predict(v)\n",
    "                            actual = [ train.target[z] for z in ind ]\n",
    "                            for i, (p, a) in enumerate(zip(pred, actual)):\n",
    "                                loss += 0 if p == a else 1\n",
    "                                sincelast += 0 if p == a else 1\n",
    "                                logpv += soft[i, a]\n",
    "                                logpvsl += soft[i, a]\n",
    "\n",
    "                            if (n & (n - 1) == 0): # and n & 0xAAAAAAAA == 0):\n",
    "                                print('{:<8d}\\t{:<8.3f} [{:<10.3f}]\\t{:<8.3f} [{:<10.3f}]\\t{:<8.3f} [{:<10.3f}]\\t{:<8.3f} [{:<10.3f}]'.format(\n",
    "                                             loss.n, effn.mean(), effnsl.mean(), avtau.mean(), avtausl.mean(), loss.mean(), sincelast.mean(), logpv.mean(), logpvsl.mean()),\n",
    "                                       flush=True)   \n",
    "                                sincelast = OnlineDRO.EasyAcc()\n",
    "                                logpvsl = OnlineDRO.EasyAcc()\n",
    "                                effnsl = OnlineDRO.EasyAcc()\n",
    "                                avtausl = OnlineDRO.EasyAcc()\n",
    "\n",
    "                        for i, z in enumerate(ind):\n",
    "                            r = 1 if logchoices[z] == train.target[z] else 0\n",
    "                            w = 1 / pchoices[z] if pred[i] == logchoices[z] else 0\n",
    "                            ocrl.update(1, w, r) \n",
    "\n",
    "                        for d in range(delay):\n",
    "                            x = np.array([ v[i] \n",
    "                                           for i, z in enumerate(ind)\n",
    "                                           if (d-1)*subblocksize <= i and i < d*subblocksize\n",
    "                                           if logchoices[z] == train.target[z] ])\n",
    "                            y = np.array([ logchoices[z] \n",
    "                                           for i, z in enumerate(ind) \n",
    "                                           if (d-1)*subblocksize <= i and i < d*subblocksize\n",
    "                                           if logchoices[z] == train.target[z] ])\n",
    "\n",
    "                            if n == 0 and pno == 0:\n",
    "                                sampweight = np.array([ lr \n",
    "                                                        for i, z in enumerate(ind) \n",
    "                                                        if (d-1)*subblocksize <= i and i < d*subblocksize\n",
    "                                                        if logchoices[z] == train.target[z] ])\n",
    "                            else:\n",
    "                                sampweight = np.array([ lr * w * ocrl.n \n",
    "                                                           * max(0, qfunc(1, w, 1))\n",
    "                                                        for i, z in enumerate(ind) \n",
    "                                                        if (d-1)*subblocksize <= i and i < d*subblocksize\n",
    "                                                        if logchoices[z] == train.target[z]\n",
    "                                                        for w in (1 / pchoices[z],)\n",
    "                                                      ])\n",
    "\n",
    "                            effn += sampweight.sum() / (lr * subblocksize)\n",
    "                            effnsl += sampweight.sum() / (lr * subblocksize)\n",
    "\n",
    "                            if np.any(x):\n",
    "                                cls.partial_fit(x, y, classes=classes, sample_weight=sampweight)\n",
    "\n",
    "                    print('{:<8d}\\t{:<8.3f} [{:<10.3f}]\\t{:<8.3f} [{:<10.3f}]\\t{:<8.3f} [{:<10.3f}]\\t{:<8.3f} [{:<10.3f}]'.format(\n",
    "                                     loss.n, effn.mean(), effnsl.mean(), avtau.mean(), avtausl.mean(), loss.mean(), sincelast.mean(), logpv.mean(), logpvsl.mean()),\n",
    "                               flush=True)   \n",
    "                    sincelast = OnlineDRO.EasyAcc()\n",
    "                    logpvsl = OnlineDRO.EasyAcc()\n",
    "                    effnsl = OnlineDRO.EasyAcc()\n",
    "                    avtausl = OnlineDRO.EasyAcc()\n",
    "\n",
    "                    preds = cls.predict(np.array([np.outer(d, np.append(d, [1])).ravel() for d in test.data]))\n",
    "                    ascores = []\n",
    "                    for b in range(16):\n",
    "                        bootie = np.random.RandomState(90210+b).choice(len(test.target), replace=True, size=len(test.target))\n",
    "                        ascores.append(accuracy_score(y_true=test.target[bootie], y_pred=preds[bootie]))\n",
    "\n",
    "                    print(\"test accuracy: {}\".format(np.quantile(ascores, [0.05, 0.5, 0.95])))\n",
    "            except KeyboardInterrupt:\n",
    "                raise\n",
    "            except:\n",
    "                pass\n",
    "    \n",
    "OnlineDRO.flass()      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "heading_collapsed": true
   },
   "source": [
    "# 20 Newsgroups\n",
    "\n",
    "Really hard ... 20 actions and only 22K examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### Full Information Online Learning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "code_folding": [
     1,
     14
    ],
    "hidden": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "**** ngram = 2 norm = l2 lr = 1 nfeat = 20 ****\n",
      "n       \temp loss\tsince last\n",
      "32      \t0.844   \t0.844     \n",
      "64      \t0.891   \t0.938     \n",
      "128     \t0.898   \t0.906     \n",
      "256     \t0.898   \t0.898     \n",
      "512     \t0.842   \t0.785     \n",
      "1024    \t0.814   \t0.787     \n",
      "2048    \t0.746   \t0.678     \n",
      "4096    \t0.661   \t0.577     \n",
      "8192    \t0.583   \t0.505     \n",
      "11264   \t0.548   \t0.456     \n",
      "test accuracy: [0.48064923 0.48771907 0.49654806]\n"
     ]
    }
   ],
   "source": [
    "class OnlineDRO:\n",
    "    class EasyAcc:\n",
    "        def __init__(self):\n",
    "            self.n = 0\n",
    "            self.sum = 0\n",
    "            \n",
    "        def __iadd__(self, other):\n",
    "            self.n += 1\n",
    "            self.sum += other\n",
    "            return self\n",
    "            \n",
    "        def mean(self):\n",
    "            return self.sum / max(self.n, 1)\n",
    " \n",
    "    def flass():\n",
    "        from scipy.sparse import vstack\n",
    "        from sklearn.datasets import fetch_20newsgroups\n",
    "        from sklearn.feature_extraction.text import HashingVectorizer\n",
    "        from sklearn.linear_model import SGDClassifier\n",
    "        from sklearn.metrics import accuracy_score\n",
    "        import numpy as np\n",
    "        \n",
    "        train = fetch_20newsgroups(subset='train', remove=('headers', 'footers', 'quotes'))\n",
    "        classes = np.unique(train.target)\n",
    "        test = fetch_20newsgroups(subset='test', remove=('headers', 'footers', 'quotes'))\n",
    "        \n",
    "        for ngram, norm, lr, nfeat in ( (2, 'l2', 1, 20), ):\n",
    "            print(\"**** ngram = {} norm = {} lr = {} nfeat = {} ****\".format(ngram, norm, lr, nfeat))\n",
    "            print('{:8.8s}\\t{:8.8s}\\t{:10.10s}'.format('n', 'emp loss', 'since last'))\n",
    "\n",
    "            vectorizer = HashingVectorizer(n_features = 1 << nfeat, norm=norm, ngram_range=(1, ngram), alternate_sign=True)\n",
    "            docs = vectorizer.transform(train.data)\n",
    "            testdocs = vectorizer.transform(test.data)\n",
    "            \n",
    "            classweights = { k: lr for k in classes }\n",
    "            cls = SGDClassifier(loss='log', class_weight=classweights, shuffle=False)\n",
    "            \n",
    "            loss = OnlineDRO.EasyAcc()\n",
    "            sincelast = OnlineDRO.EasyAcc()\n",
    "            blocksize = 32\n",
    "\n",
    "            for pno in range(1):\n",
    "                order = np.random.RandomState(seed=42+pno).permutation(len(train.data))\n",
    "                for n, ind in enumerate(zip(*(iter(order),)*blocksize)):\n",
    "                    v = vstack([ docs[z] for z in ind ])\n",
    "                    actual = [ train.target[z] for z in ind ]\n",
    "                    if n > 0:\n",
    "                        pred = cls.predict(v)\n",
    "                        for p, a in zip(pred, actual):\n",
    "                            loss += 0 if p == a else 1\n",
    "                            sincelast += 0 if p == a else 1\n",
    "                        if (n & (n - 1) == 0): # and n & 0xAAAAAAAA == 0):\n",
    "                            print('{:<8d}\\t{:<8.3f}\\t{:<10.3f}'.format(loss.n, loss.mean(), sincelast.mean()), flush=True)\n",
    "                            sincelast = OnlineDRO.EasyAcc()\n",
    "\n",
    "                    cls.partial_fit(v, actual, classes=classes)\n",
    "\n",
    "                print('{:<8d}\\t{:<8.3f}\\t{:<10.3f}'.format(loss.n, loss.mean(), sincelast.mean()), flush=True)\n",
    "                sincelast = OnlineDRO.EasyAcc()\n",
    "\n",
    "                preds = cls.predict(testdocs)\n",
    "                ascores = []\n",
    "                for b in range(16):\n",
    "                    bootie = np.random.RandomState(90210+b).choice(len(test.target), replace=True, size=len(test.target))\n",
    "                    ascores.append(accuracy_score(y_true=test.target[bootie], y_pred=preds[bootie]))\n",
    "\n",
    "                print(\"test accuracy: {}\".format(np.quantile(ascores, [0.05, 0.5, 0.95])))\n",
    "    \n",
    "OnlineDRO.flass()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### Partial Information Online Learning, Softmax Logging Policy\n",
    "\n",
    "Uniform $(\\tau = 0)$ and softmax $(\\tau = 4)$ are pretty similar for off-policy learning but uniform has larger regret."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "code_folding": [
     1,
     14
    ],
    "hidden": true,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "**** ngram = 2 norm = l2 lr = 0.5 nfeat = 20 tau = 0 ****\n",
      "n       \temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.969   \t0.969     \t0.050   \t0.050     \n",
      "64      \t0.906   \t0.844     \t0.050   \t0.050     \n",
      "128     \t0.938   \t0.969     \t0.050   \t0.050     \n",
      "256     \t0.945   \t0.953     \t0.050   \t0.050     \n",
      "512     \t0.930   \t0.914     \t0.050   \t0.050     \n",
      "1024    \t0.927   \t0.924     \t0.050   \t0.050     \n",
      "2048    \t0.927   \t0.928     \t0.050   \t0.050     \n",
      "4096    \t0.921   \t0.915     \t0.050   \t0.050     \n",
      "8192    \t0.897   \t0.873     \t0.050   \t0.050     \n",
      "11264   \t0.881   \t0.838     \t0.050   \t0.050     \n",
      "test accuracy: [0.15676447 0.16370154 0.17113648]\n",
      "**** ngram = 2 norm = l2 lr = 0.5 nfeat = 20 tau = 4 ****\n",
      "n       \temp loss\tsince last\tlog pv  \tsince last\n",
      "32      \t0.969   \t0.969     \t0.050   \t0.050     \n",
      "64      \t0.906   \t0.844     \t0.050   \t0.050     \n",
      "128     \t0.938   \t0.969     \t0.050   \t0.050     \n",
      "256     \t0.945   \t0.953     \t0.050   \t0.050     \n",
      "512     \t0.926   \t0.906     \t0.050   \t0.050     \n",
      "1024    \t0.930   \t0.934     \t0.050   \t0.051     \n",
      "2048    \t0.923   \t0.917     \t0.051   \t0.051     \n",
      "4096    \t0.916   \t0.908     \t0.051   \t0.052     \n",
      "8192    \t0.894   \t0.873     \t0.053   \t0.054     \n",
      "11264   \t0.879   \t0.839     \t0.054   \t0.056     \n",
      "test accuracy: [0.17455523 0.18102762 0.19095194]\n"
     ]
    }
   ],
   "source": [
    "class OnlineDRO:\n",
    "    class EasyAcc:\n",
    "        def __init__(self):\n",
    "            self.n = 0\n",
    "            self.sum = 0\n",
    "            \n",
    "        def __iadd__(self, other):\n",
    "            self.n += 1\n",
    "            self.sum += other\n",
    "            return self\n",
    "            \n",
    "        def mean(self):\n",
    "            return self.sum / max(self.n, 1)\n",
    " \n",
    "    def flass():\n",
    "        from scipy.sparse import vstack\n",
    "        from scipy.special import softmax\n",
    "        from sklearn.datasets import fetch_20newsgroups\n",
    "        from sklearn.feature_extraction.text import HashingVectorizer\n",
    "        from sklearn.linear_model import SGDClassifier\n",
    "        from sklearn.metrics import accuracy_score\n",
    "        import numpy as np\n",
    "        \n",
    "        train = fetch_20newsgroups(subset='train', remove=('headers', 'footers', 'quotes'))\n",
    "        classes = np.unique(train.target)\n",
    "        test = fetch_20newsgroups(subset='test', remove=('headers', 'footers', 'quotes'))\n",
    "        \n",
    "        blocksize = 32\n",
    "        for ngram, norm, lr, nfeat, tau in ( (2, 'l2', z, 20, x) for x in (0, 4,) for z in (0.5, ) ):\n",
    "            print(\"**** ngram = {} norm = {} lr = {} nfeat = {} tau = {} ****\".format(ngram, norm, lr, nfeat, tau))\n",
    "            print('{:8.8s}\\t{:8.8s}\\t{:10.10s}\\t{:8.8s}\\t{:10.10s}'.format(\n",
    "                'n', 'emp loss', 'since last', 'log pv', 'since last')\n",
    "            )\n",
    "            \n",
    "            vectorizer = HashingVectorizer(n_features = 1 << nfeat, norm=norm, ngram_range=(1, ngram), alternate_sign=True)\n",
    "            docs = vectorizer.transform(train.data)\n",
    "            testdocs = vectorizer.transform(test.data)\n",
    "            \n",
    "            classweights = { k: lr for k in classes }\n",
    "            cls = SGDClassifier(loss='log', class_weight=classweights, shuffle=False)\n",
    "            \n",
    "            loss = OnlineDRO.EasyAcc()\n",
    "            sincelast = OnlineDRO.EasyAcc()\n",
    "            logpv = OnlineDRO.EasyAcc()\n",
    "            logpvsl = OnlineDRO.EasyAcc()\n",
    "            \n",
    "            loggerrand = np.random.RandomState(seed=2112)\n",
    "            logchoices = [None]*len(train.data)\n",
    "            pchoices = [None]*len(train.data)\n",
    " \n",
    "            for pno in range(1):\n",
    "                order = np.random.RandomState(seed=42+pno).permutation(len(train.data))\n",
    "                for n, ind in enumerate(zip(*(iter(order),)*blocksize)):\n",
    "                    v = vstack([ docs[z] for z in ind ])\n",
    "                    if n == 0 and pno == 0:\n",
    "                        for i, z in enumerate(ind):\n",
    "                            if logchoices[z] is None:\n",
    "                                choice = loggerrand.choice(a=classes, size=1)\n",
    "                                logchoices[z] = choice[0]\n",
    "                                pchoices[z] = 1.0 / len(classes)\n",
    "                    else:\n",
    "                        predlogp = cls.predict_proba(v)\n",
    "                        soft = softmax(tau * predlogp, axis=1)\n",
    "\n",
    "                        for i, z in enumerate(ind):\n",
    "                            if logchoices[z] is None:\n",
    "                                choice = loggerrand.choice(a=classes, p=soft[i,:], size=1)\n",
    "                                logchoices[z] = choice[0]\n",
    "                                pchoices[z] = soft[i, choice[0]]\n",
    "\n",
    "                        pred = cls.predict(v)\n",
    "                        actual = [ train.target[z] for z in ind ]\n",
    "                        for i, (p, a) in enumerate(zip(pred, actual)):\n",
    "                            loss += 0 if p == a else 1\n",
    "                            sincelast += 0 if p == a else 1\n",
    "                            logpv += soft[i, a]\n",
    "                            logpvsl += soft[i, a]\n",
    "                            \n",
    "                        if (n & (n - 1) == 0): # and n & 0xAAAAAAAA == 0):\n",
    "                            print('{:<8d}\\t{:<8.3f}\\t{:<10.3f}\\t{:<8.3f}\\t{:<10.3f}'.format(\n",
    "                                        loss.n, loss.mean(), sincelast.mean(), logpv.mean(), logpvsl.mean()),\n",
    "                                  flush=True)\n",
    "\n",
    "                            sincelast = OnlineDRO.EasyAcc()\n",
    "                            logpvsl = OnlineDRO.EasyAcc()\n",
    "\n",
    "#                     y = np.array([ logchoices[z] for i, z in enumerate(ind) ])\n",
    "#                     w = np.array([ (lr / len(classes)) * (1 / pchoices[z]) \n",
    "#                                     if logchoices[z] == train.target[z]\n",
    "#                                     else -(lr / len(classes)**2) * (1/pchoices[z])\n",
    "#                                    for i, z in enumerate(ind) ])\n",
    "                    \n",
    "                    y = np.array([ logchoices[z] for i, z in enumerate(ind) if logchoices[z] == train.target[z] ])\n",
    "                    w = np.array([ (lr / len(classes)) * (1 / pchoices[z]) \n",
    "                                   for i, z in enumerate(ind) if logchoices[z] == train.target[z] ])\n",
    "                    if np.any(y):\n",
    "                        x = vstack([ v[i] for i, z in enumerate(ind) if logchoices[z] == train.target[z] ])\n",
    "#                         x = vstack([ v[i] for i, z in enumerate(ind) ])\n",
    "                        cls.partial_fit(x, y, classes=classes, sample_weight=w)\n",
    "\n",
    "                print('{:<8d}\\t{:<8.3f}\\t{:<10.3f}\\t{:<8.3f}\\t{:<10.3f}'.format(\n",
    "                             loss.n, loss.mean(), sincelast.mean(), logpv.mean(), logpvsl.mean()),\n",
    "                       flush=True)                \n",
    "                sincelast = OnlineDRO.EasyAcc()\n",
    "                logpvsl = OnlineDRO.EasyAcc()\n",
    "\n",
    "                preds = cls.predict(testdocs)\n",
    "                ascores = []\n",
    "                for b in range(16):\n",
    "                    bootie = np.random.RandomState(90210+b).choice(len(test.target), replace=True, size=len(test.target))\n",
    "                    ascores.append(accuracy_score(y_true=test.target[bootie], y_pred=preds[bootie]))\n",
    "\n",
    "                print(\"test accuracy: {}\".format(np.quantile(ascores, [0.05, 0.5, 0.95])))\n",
    "    \n",
    "OnlineDRO.flass()      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hidden": true
   },
   "source": [
    "### Auto-Temperature and Bound Optimization\n",
    "\n",
    "Nothing working great, but regret is better than other approaches. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": [
     0,
     1,
     14,
     90
    ],
    "hidden": true,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "*** lr = 0.5 taumax = 8 target = 0.025 what = mle gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.586    [1.136     ]\t8.000    [8.000     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t1.075    [1.565     ]\t8.000    [8.000     ]\t0.938    [0.939     ]\t0.060    [0.062     ]\n",
      "2048    \t0.618    [0.161     ]\t8.000    [8.000     ]\t0.941    [0.945     ]\t0.056    [0.053     ]\n",
      "4096    \t0.508    [0.398     ]\t8.000    [8.000     ]\t0.942    [0.943     ]\t0.057    [0.059     ]\n",
      "8192    \t0.296    [0.084     ]\t8.000    [8.000     ]\t0.942    [0.941     ]\t0.058    [0.058     ]\n",
      "11008   \t0.318    [0.378     ]\t8.000    [8.000     ]\t0.942    [0.942     ]\t0.058    [0.059     ]\n",
      "test accuracy: [0.05961232 0.064923   0.06880643]\n",
      "*** lr = 0.5 taumax = 8 target = 0.05 what = mle gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.586    [1.136     ]\t8.000    [8.000     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t1.075    [1.565     ]\t8.000    [8.000     ]\t0.938    [0.939     ]\t0.060    [0.062     ]\n",
      "2048    \t0.616    [0.157     ]\t7.965    [7.930     ]\t0.942    [0.947     ]\t0.056    [0.053     ]\n",
      "4096    \t0.718    [0.820     ]\t7.863    [7.761     ]\t0.945    [0.947     ]\t0.056    [0.055     ]\n",
      "8192    \t0.556    [0.394     ]\t7.780    [7.697     ]\t0.948    [0.951     ]\t0.053    [0.050     ]\n",
      "11008   \t0.548    [0.528     ]\t7.802    [7.867     ]\t0.946    [0.939     ]\t0.054    [0.058     ]\n",
      "test accuracy: [0.04855948 0.05031864 0.05420207]\n",
      "*** lr = 0.5 taumax = 8 target = 0.1 what = mle gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.586    [1.136     ]\t8.000    [8.000     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t1.075    [1.565     ]\t7.743    [7.486     ]\t0.938    [0.939     ]\t0.060    [0.063     ]\n",
      "2048    \t0.599    [0.123     ]\t7.118    [6.493     ]\t0.943    [0.949     ]\t0.054    [0.049     ]\n",
      "4096    \t0.668    [0.738     ]\t7.047    [6.975     ]\t0.941    [0.939     ]\t0.056    [0.058     ]\n",
      "8192    \t0.654    [0.640     ]\t6.711    [6.375     ]\t0.942    [0.942     ]\t0.057    [0.059     ]\n",
      "11008   \t0.677    [0.738     ]\t6.746    [6.849     ]\t0.936    [0.919     ]\t0.063    [0.077     ]\n",
      "test accuracy: [0.05187865 0.0568906  0.06143787]\n",
      "*** lr = 0.5 taumax = 8 target = 0.2 what = mle gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.586    [1.136     ]\t7.118    [6.236     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t0.799    [1.013     ]\t6.870    [6.622     ]\t0.947    [0.959     ]\t0.050    [0.044     ]\n",
      "2048    \t0.836    [0.873     ]\t6.247    [5.624     ]\t0.952    [0.957     ]\t0.050    [0.051     ]\n",
      "4096    \t0.948    [1.059     ]\t6.261    [6.275     ]\t0.949    [0.945     ]\t0.052    [0.054     ]\n",
      "8192    \t0.865    [0.781     ]\t5.838    [5.415     ]\t0.942    [0.935     ]\t0.056    [0.061     ]\n",
      "11008   \t0.850    [0.813     ]\t5.749    [5.492     ]\t0.937    [0.925     ]\t0.060    [0.070     ]\n",
      "test accuracy: [0.0673128  0.06950345 0.07527881]\n",
      "*** lr = 0.5 taumax = 8 target = 0.025 what = mle gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.586    [1.136     ]\t8.000    [8.000     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t0.427    [0.269     ]\t8.000    [8.000     ]\t0.938    [0.939     ]\t0.059    [0.061     ]\n",
      "2048    \t0.251    [0.075     ]\t8.000    [8.000     ]\t0.943    [0.949     ]\t0.055    [0.051     ]\n",
      "4096    \t0.207    [0.162     ]\t8.000    [8.000     ]\t0.947    [0.950     ]\t0.052    [0.049     ]\n",
      "8192    \t0.173    [0.139     ]\t8.000    [8.000     ]\t0.946    [0.945     ]\t0.053    [0.055     ]\n",
      "11008   \t0.150    [0.088     ]\t8.000    [8.000     ]\t0.945    [0.944     ]\t0.053    [0.054     ]\n",
      "test accuracy: [0.0545008  0.05662507 0.05835104]\n",
      "*** lr = 0.5 taumax = 8 target = 0.05 what = mle gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.586    [1.136     ]\t8.000    [8.000     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t0.427    [0.269     ]\t8.000    [8.000     ]\t0.938    [0.939     ]\t0.059    [0.061     ]\n",
      "2048    \t0.251    [0.075     ]\t8.000    [8.000     ]\t0.943    [0.949     ]\t0.055    [0.051     ]\n",
      "4096    \t0.207    [0.162     ]\t8.000    [8.000     ]\t0.947    [0.950     ]\t0.052    [0.049     ]\n",
      "8192    \t0.228    [0.250     ]\t7.846    [7.691     ]\t0.946    [0.945     ]\t0.053    [0.054     ]\n",
      "11008   \t0.438    [0.997     ]\t7.809    [7.702     ]\t0.945    [0.943     ]\t0.054    [0.056     ]\n",
      "test accuracy: [0.05018587 0.05250929 0.05639272]\n",
      "*** lr = 0.5 taumax = 8 target = 0.1 what = mle gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.586    [1.136     ]\t8.000    [8.000     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t0.427    [0.269     ]\t8.000    [8.000     ]\t0.938    [0.939     ]\t0.059    [0.061     ]\n",
      "2048    \t0.251    [0.075     ]\t7.930    [7.861     ]\t0.943    [0.949     ]\t0.055    [0.051     ]\n",
      "4096    \t0.186    [0.121     ]\t7.907    [7.884     ]\t0.949    [0.955     ]\t0.051    [0.047     ]\n",
      "8192    \t0.172    [0.157     ]\t7.907    [7.907     ]\t0.941    [0.933     ]\t0.057    [0.063     ]\n",
      "11008   \t0.293    [0.615     ]\t7.815    [7.548     ]\t0.941    [0.940     ]\t0.058    [0.061     ]\n",
      "test accuracy: [0.06668216 0.07249071 0.0765401 ]\n",
      "*** lr = 0.5 taumax = 8 target = 0.2 what = mle gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.586    [1.136     ]\t7.118    [6.236     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t0.799    [1.013     ]\t6.877    [6.635     ]\t0.947    [0.959     ]\t0.050    [0.044     ]\n",
      "2048    \t1.039    [1.279     ]\t6.112    [5.346     ]\t0.951    [0.955     ]\t0.051    [0.052     ]\n",
      "4096    \t0.751    [0.462     ]\t6.240    [6.368     ]\t0.943    [0.936     ]\t0.057    [0.063     ]\n",
      "8192    \t0.740    [0.728     ]\t5.849    [5.457     ]\t0.946    [0.948     ]\t0.055    [0.052     ]\n",
      "11008   \t0.740    [0.742     ]\t5.688    [5.222     ]\t0.943    [0.934     ]\t0.057    [0.062     ]\n",
      "test accuracy: [0.04139007 0.04454328 0.047265  ]\n",
      "*** lr = 0.5 taumax = 8 target = 0.025 what = ub gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.597    [1.159     ]\t8.000    [8.000     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t0.484    [0.370     ]\t8.000    [8.000     ]\t0.932    [0.928     ]\t0.059    [0.062     ]\n",
      "2048    \t0.373    [0.262     ]\t8.000    [8.000     ]\t0.938    [0.944     ]\t0.058    [0.056     ]\n",
      "4096    \t0.342    [0.312     ]\t8.000    [8.000     ]\t0.939    [0.939     ]\t0.059    [0.061     ]\n",
      "8192    \t0.441    [0.540     ]\t8.000    [8.000     ]\t0.936    [0.933     ]\t0.064    [0.068     ]\n",
      "11008   \t0.433    [0.410     ]\t8.000    [8.000     ]\t0.939    [0.949     ]\t0.062    [0.056     ]\n",
      "test accuracy: [0.05065056 0.05337228 0.05625996]\n",
      "*** lr = 0.5 taumax = 8 target = 0.05 what = ub gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.597    [1.159     ]\t8.000    [8.000     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t0.484    [0.370     ]\t8.000    [8.000     ]\t0.932    [0.928     ]\t0.059    [0.062     ]\n",
      "2048    \t0.373    [0.262     ]\t8.000    [8.000     ]\t0.938    [0.944     ]\t0.058    [0.056     ]\n",
      "4096    \t0.342    [0.312     ]\t8.000    [8.000     ]\t0.939    [0.939     ]\t0.059    [0.061     ]\n",
      "8192    \t0.441    [0.540     ]\t8.000    [8.000     ]\t0.936    [0.933     ]\t0.064    [0.068     ]\n",
      "11008   \t0.477    [0.575     ]\t7.857    [7.441     ]\t0.938    [0.945     ]\t0.062    [0.057     ]\n",
      "test accuracy: [0.04779607 0.05131439 0.05592804]\n",
      "*** lr = 0.5 taumax = 8 target = 0.1 what = ub gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.597    [1.159     ]\t8.000    [8.000     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t0.484    [0.370     ]\t8.000    [8.000     ]\t0.932    [0.928     ]\t0.059    [0.062     ]\n",
      "2048    \t0.463    [0.442     ]\t7.430    [6.861     ]\t0.940    [0.948     ]\t0.056    [0.052     ]\n",
      "4096    \t1.107    [1.752     ]\t7.044    [6.658     ]\t0.944    [0.947     ]\t0.054    [0.053     ]\n",
      "8192    \t0.874    [0.640     ]\t6.865    [6.685     ]\t0.945    [0.947     ]\t0.054    [0.054     ]\n",
      "11008   \t0.761    [0.460     ]\t6.779    [6.529     ]\t0.945    [0.945     ]\t0.054    [0.053     ]\n",
      "test accuracy: [0.0520778  0.05622677 0.06037573]\n",
      "*** lr = 0.5 taumax = 8 target = 0.2 what = ub gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.597    [1.159     ]\t7.112    [6.224     ]\t0.936    [0.938     ]\t0.056    [0.062     ]\n",
      "1024    \t0.816    [1.036     ]\t6.855    [6.597     ]\t0.947    [0.959     ]\t0.050    [0.044     ]\n",
      "2048    \t0.862    [0.908     ]\t6.141    [5.428     ]\t0.951    [0.954     ]\t0.051    [0.052     ]\n",
      "4096    \t0.933    [1.004     ]\t6.443    [6.745     ]\t0.946    [0.942     ]\t0.054    [0.058     ]\n",
      "8192    \t0.905    [0.877     ]\t5.934    [5.424     ]\t0.941    [0.935     ]\t0.059    [0.064     ]\n",
      "11008   \t0.895    [0.867     ]\t5.740    [5.176     ]\t0.939    [0.936     ]\t0.060    [0.062     ]\n",
      "test accuracy: [0.06040892 0.064923   0.06960303]\n",
      "*** lr = 0.5 taumax = 8 target = 0.025 what = ub gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.599    [1.163     ]\t8.000    [8.000     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t0.486    [0.373     ]\t8.000    [8.000     ]\t0.932    [0.928     ]\t0.059    [0.062     ]\n",
      "2048    \t0.456    [0.427     ]\t8.000    [8.000     ]\t0.928    [0.925     ]\t0.069    [0.078     ]\n",
      "4096    \t0.407    [0.357     ]\t8.000    [8.000     ]\t0.928    [0.928     ]\t0.068    [0.068     ]\n",
      "8192    \t0.329    [0.250     ]\t8.000    [8.000     ]\t0.929    [0.930     ]\t0.068    [0.069     ]\n",
      "11008   \t0.347    [0.397     ]\t8.000    [8.000     ]\t0.929    [0.929     ]\t0.069    [0.069     ]\n",
      "test accuracy: [0.08576739 0.09087892 0.09323553]\n",
      "*** lr = 0.5 taumax = 8 target = 0.05 what = ub gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.599    [1.163     ]\t8.000    [8.000     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t0.486    [0.373     ]\t8.000    [8.000     ]\t0.932    [0.928     ]\t0.059    [0.062     ]\n",
      "2048    \t0.456    [0.427     ]\t8.000    [8.000     ]\t0.928    [0.925     ]\t0.069    [0.078     ]\n",
      "4096    \t0.407    [0.357     ]\t8.000    [8.000     ]\t0.928    [0.928     ]\t0.068    [0.068     ]\n",
      "8192    \t0.385    [0.364     ]\t7.988    [7.975     ]\t0.930    [0.933     ]\t0.067    [0.066     ]\n",
      "11008   \t0.550    [0.990     ]\t7.946    [7.826     ]\t0.934    [0.944     ]\t0.065    [0.058     ]\n",
      "test accuracy: [0.05592804 0.05821827 0.06346256]\n",
      "*** lr = 0.5 taumax = 8 target = 0.1 what = ub gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.599    [1.163     ]\t8.000    [8.000     ]\t0.936    [0.938     ]\t0.057    [0.062     ]\n",
      "1024    \t0.486    [0.373     ]\t8.000    [8.000     ]\t0.932    [0.928     ]\t0.059    [0.062     ]\n",
      "2048    \t0.457    [0.427     ]\t7.999    [7.998     ]\t0.928    [0.925     ]\t0.069    [0.078     ]\n",
      "4096    \t0.693    [0.930     ]\t7.489    [6.980     ]\t0.930    [0.932     ]\t0.067    [0.066     ]\n",
      "8192    \t0.623    [0.553     ]\t7.183    [6.876     ]\t0.935    [0.941     ]\t0.064    [0.060     ]\n",
      "11008   \t0.752    [1.097     ]\t7.174    [7.149     ]\t0.938    [0.944     ]\t0.062    [0.055     ]\n",
      "test accuracy: [0.05476633 0.05828465 0.06220127]\n",
      "*** lr = 0.5 taumax = 8 target = 0.2 what = ub gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.599    [1.163     ]\t7.111    [6.223     ]\t0.936    [0.938     ]\t0.056    [0.062     ]\n",
      "1024    \t0.817    [1.036     ]\t6.853    [6.595     ]\t0.947    [0.959     ]\t0.050    [0.044     ]\n",
      "2048    \t0.678    [0.538     ]\t6.218    [5.584     ]\t0.947    [0.947     ]\t0.052    [0.054     ]\n",
      "4096    \t0.727    [0.777     ]\t5.803    [5.388     ]\t0.946    [0.945     ]\t0.053    [0.055     ]\n",
      "8192    \t0.789    [0.851     ]\t5.542    [5.281     ]\t0.937    [0.927     ]\t0.061    [0.068     ]\n",
      "11008   \t0.811    [0.870     ]\t5.512    [5.422     ]\t0.930    [0.911     ]\t0.065    [0.079     ]\n",
      "test accuracy: [0.07152815 0.07534519 0.08078864]\n",
      "*** lr = 0.5 taumax = 8 target = 0.025 what = lb gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.018    [0.000     ]\t8.000    [8.000     ]\t0.930    [0.926     ]\t0.051    [0.051     ]\n",
      "1024    \t0.009    [0.000     ]\t8.000    [8.000     ]\t0.926    [0.922     ]\t0.051    [0.051     ]\n",
      "2048    \t0.023    [0.037     ]\t8.000    [8.000     ]\t0.938    [0.950     ]\t0.052    [0.053     ]\n",
      "4096    \t0.031    [0.040     ]\t8.000    [8.000     ]\t0.944    [0.950     ]\t0.054    [0.057     ]\n",
      "8192    \t0.036    [0.042     ]\t8.000    [8.000     ]\t0.946    [0.947     ]\t0.056    [0.057     ]\n",
      "11008   \t0.038    [0.042     ]\t8.000    [8.000     ]\t0.946    [0.946     ]\t0.056    [0.056     ]\n",
      "test accuracy: [0.05028545 0.05250929 0.05755443]\n",
      "*** lr = 0.5 taumax = 8 target = 0.05 what = lb gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.018    [0.000     ]\t8.000    [8.000     ]\t0.930    [0.926     ]\t0.051    [0.051     ]\n",
      "1024    \t0.009    [0.000     ]\t8.000    [8.000     ]\t0.926    [0.922     ]\t0.051    [0.051     ]\n",
      "2048    \t0.023    [0.037     ]\t8.000    [8.000     ]\t0.938    [0.950     ]\t0.052    [0.053     ]\n",
      "4096    \t0.031    [0.040     ]\t8.000    [8.000     ]\t0.944    [0.950     ]\t0.054    [0.057     ]\n",
      "8192    \t0.036    [0.042     ]\t8.000    [8.000     ]\t0.946    [0.947     ]\t0.056    [0.057     ]\n",
      "11008   \t0.038    [0.042     ]\t8.000    [8.000     ]\t0.946    [0.946     ]\t0.056    [0.056     ]\n",
      "test accuracy: [0.05028545 0.05250929 0.05755443]\n",
      "*** lr = 0.5 taumax = 8 target = 0.1 what = lb gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.018    [0.000     ]\t8.000    [8.000     ]\t0.930    [0.926     ]\t0.051    [0.051     ]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1024    \t0.009    [0.000     ]\t8.000    [8.000     ]\t0.926    [0.922     ]\t0.051    [0.051     ]\n",
      "2048    \t0.023    [0.037     ]\t8.000    [8.000     ]\t0.938    [0.950     ]\t0.052    [0.053     ]\n",
      "4096    \t0.031    [0.040     ]\t8.000    [8.000     ]\t0.944    [0.950     ]\t0.054    [0.057     ]\n",
      "8192    \t0.036    [0.042     ]\t8.000    [8.000     ]\t0.946    [0.947     ]\t0.056    [0.057     ]\n",
      "11008   \t0.038    [0.042     ]\t8.000    [8.000     ]\t0.946    [0.946     ]\t0.056    [0.056     ]\n",
      "test accuracy: [0.05028545 0.05250929 0.05755443]\n",
      "*** lr = 0.5 taumax = 8 target = 0.2 what = lb gamma = 0.9999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.018    [0.000     ]\t8.000    [8.000     ]\t0.930    [0.926     ]\t0.051    [0.051     ]\n",
      "1024    \t0.009    [0.000     ]\t8.000    [8.000     ]\t0.926    [0.922     ]\t0.051    [0.051     ]\n",
      "2048    \t0.023    [0.037     ]\t8.000    [8.000     ]\t0.938    [0.950     ]\t0.052    [0.053     ]\n",
      "4096    \t0.031    [0.040     ]\t8.000    [8.000     ]\t0.944    [0.950     ]\t0.054    [0.057     ]\n",
      "8192    \t0.036    [0.042     ]\t7.814    [7.628     ]\t0.946    [0.947     ]\t0.056    [0.057     ]\n",
      "11008   \t0.038    [0.043     ]\t7.760    [7.603     ]\t0.946    [0.946     ]\t0.056    [0.057     ]\n",
      "test accuracy: [0.05028545 0.05250929 0.05755443]\n",
      "*** lr = 0.5 taumax = 8 target = 0.025 what = lb gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.018    [0.000     ]\t8.000    [8.000     ]\t0.930    [0.926     ]\t0.051    [0.051     ]\n",
      "1024    \t0.009    [0.000     ]\t8.000    [8.000     ]\t0.926    [0.922     ]\t0.051    [0.051     ]\n",
      "2048    \t0.017    [0.025     ]\t8.000    [8.000     ]\t0.938    [0.950     ]\t0.052    [0.053     ]\n",
      "4096    \t0.026    [0.034     ]\t8.000    [8.000     ]\t0.944    [0.950     ]\t0.055    [0.057     ]\n",
      "8192    \t0.030    [0.034     ]\t8.000    [8.000     ]\t0.946    [0.947     ]\t0.056    [0.057     ]\n",
      "11008   \t0.032    [0.037     ]\t8.000    [8.000     ]\t0.946    [0.946     ]\t0.056    [0.056     ]\n",
      "test accuracy: [0.05028545 0.05250929 0.05755443]\n",
      "*** lr = 0.5 taumax = 8 target = 0.05 what = lb gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.018    [0.000     ]\t8.000    [8.000     ]\t0.930    [0.926     ]\t0.051    [0.051     ]\n",
      "1024    \t0.009    [0.000     ]\t8.000    [8.000     ]\t0.926    [0.922     ]\t0.051    [0.051     ]\n",
      "2048    \t0.017    [0.025     ]\t8.000    [8.000     ]\t0.938    [0.950     ]\t0.052    [0.053     ]\n",
      "4096    \t0.026    [0.034     ]\t8.000    [8.000     ]\t0.944    [0.950     ]\t0.055    [0.057     ]\n",
      "8192    \t0.030    [0.034     ]\t8.000    [8.000     ]\t0.946    [0.947     ]\t0.056    [0.057     ]\n",
      "11008   \t0.032    [0.037     ]\t8.000    [8.000     ]\t0.946    [0.946     ]\t0.056    [0.056     ]\n",
      "test accuracy: [0.05028545 0.05250929 0.05755443]\n",
      "*** lr = 0.5 taumax = 8 target = 0.1 what = lb gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.018    [0.000     ]\t8.000    [8.000     ]\t0.930    [0.926     ]\t0.051    [0.051     ]\n",
      "1024    \t0.009    [0.000     ]\t8.000    [8.000     ]\t0.926    [0.922     ]\t0.051    [0.051     ]\n",
      "2048    \t0.017    [0.025     ]\t8.000    [8.000     ]\t0.938    [0.950     ]\t0.052    [0.053     ]\n",
      "4096    \t0.026    [0.034     ]\t8.000    [8.000     ]\t0.944    [0.950     ]\t0.055    [0.057     ]\n",
      "8192    \t0.030    [0.034     ]\t8.000    [8.000     ]\t0.946    [0.947     ]\t0.056    [0.057     ]\n",
      "11008   \t0.032    [0.037     ]\t8.000    [8.000     ]\t0.946    [0.946     ]\t0.056    [0.056     ]\n",
      "test accuracy: [0.05028545 0.05250929 0.05755443]\n",
      "*** lr = 0.5 taumax = 8 target = 0.2 what = lb gamma = 0.999 ***\n",
      "n       \teff n    [since last]\tav tau   [since last]\temp loss [since last]\tlog pv   [since last]\n",
      "256     \t0.035    [0.035     ]\t8.000    [8.000     ]\t0.934    [0.934     ]\t0.051    [0.051     ]\n",
      "512     \t0.018    [0.000     ]\t8.000    [8.000     ]\t0.930    [0.926     ]\t0.051    [0.051     ]\n",
      "1024    \t0.009    [0.000     ]\t8.000    [8.000     ]\t0.926    [0.922     ]\t0.051    [0.051     ]\n",
      "2048    \t0.017    [0.025     ]\t8.000    [8.000     ]\t0.938    [0.950     ]\t0.052    [0.053     ]\n",
      "4096    \t0.026    [0.034     ]\t8.000    [8.000     ]\t0.944    [0.950     ]\t0.055    [0.057     ]\n",
      "8192    \t0.029    [0.033     ]\t7.961    [7.923     ]\t0.946    [0.947     ]\t0.056    [0.057     ]\n",
      "11008   \t0.032    [0.037     ]\t7.949    [7.913     ]\t0.946    [0.946     ]\t0.056    [0.056     ]\n",
      "test accuracy: [0.05028545 0.05250929 0.05755443]\n"
     ]
    }
   ],
   "source": [
    "class OnlineDRO:\n",
    "    class EasyAcc:\n",
    "        def __init__(self):\n",
    "            self.n = 0\n",
    "            self.sum = 0\n",
    "            \n",
    "        def __iadd__(self, other):\n",
    "            self.n += 1\n",
    "            self.sum += other\n",
    "            return self\n",
    "            \n",
    "        def mean(self):\n",
    "            return self.sum / max(self.n, 1)\n",
    "\n",
    "    class OnlineCressieReadLB:\n",
    "        from math import inf\n",
    "        \n",
    "        def __init__(self, alpha, gamma=1, wmin=0, wmax=inf):\n",
    "            import numpy as np\n",
    "            \n",
    "            self.alpha = alpha\n",
    "            self.gamma = gamma\n",
    "            self.n = 0\n",
    "            self.sumw = 0\n",
    "            self.sumwsq = 0\n",
    "            self.sumwr = 0\n",
    "            self.sumwsqr = 0\n",
    "            self.sumwsqrsq = 0\n",
    "            self.wmin = wmin\n",
    "            self.wmax = wmax\n",
    "            \n",
    "            self.duals = None\n",
    "            self.mleduals = None\n",
    "            \n",
    "        def update(self, c, w, r):\n",
    "            if c > 0:\n",
    "                assert w + 1e-6 >= self.wmin and w <= self.wmax + 1e-6, 'w = {} < {} < {}'.format(self.wmin, w, self.wmax)\n",
    "                assert r >= 0 and r <= 1, 'r = {}'.format(r)\n",
    "                \n",
    "                decay = self.gamma ** c\n",
    "                self.n = decay * self.n + c\n",
    "                self.sumw = decay * self.sumw + c * w\n",
    "                self.sumwsq = decay * self.sumwsq + c * w**2\n",
    "                self.sumwr = decay * self.sumwr + c * w * r\n",
    "                self.sumwsqr = decay * self.sumwsqr + c * (w**2) * r\n",
    "                self.sumwsqrsq = decay * self.sumwsqrsq + c * (w**2) * (r**2)\n",
    "                    \n",
    "                self.duals = None\n",
    "                self.mleduals = None\n",
    "                \n",
    "            return self\n",
    "        \n",
    "        def recomputeduals(self):\n",
    "            from MLE.MLE import CrMinusTwo as CrMinusTwo\n",
    "            \n",
    "            self.duals = CrMinusTwo.intervalimpl(self.n, self.sumw, self.sumwsq, \n",
    "                                                 self.sumwr, self.sumwsqr, self.sumwsqrsq,\n",
    "                                                 self.wmin, self.wmax, self.alpha, raiseonerr=True)\n",
    "            \n",
    "        def recomputedualsmle(self):\n",
    "            from MLE.MLE import CrMinusTwo as CrMinusTwo\n",
    "            \n",
    "            self.mleduals = CrMinusTwo.estimateimpl(self.n, self.sumw, self.sumwsq, \n",
    "                                                    self.sumwr, self.sumwsqr, None, None,\n",
    "                                                    self.wmin, self.wmax, raiseonerr=True)\n",
    "        \n",
    "        def qlb(self, c, w, r):\n",
    "            if self.duals is None:\n",
    "                self.recomputeduals()\n",
    "                \n",
    "                assert self.duals is not None\n",
    "                \n",
    "            return self.duals[1][0]['qfunc'](c, w, r) if self.duals[1][0] is not None else 1\n",
    "        \n",
    "        def qub(self, c, w, r):\n",
    "            if self.duals is None:\n",
    "                self.recomputeduals()\n",
    "                \n",
    "                assert self.duals is not None\n",
    "                \n",
    "            return self.duals[1][1]['qfunc'](c, w, r) if self.duals[1][1] is not None else 1\n",
    "        \n",
    "        def qmle(self, c, w, r):\n",
    "            if self.mleduals is None:\n",
    "                self.recomputedualsmle()\n",
    "                \n",
    "                assert self.mleduals is not None\n",
    "                \n",
    "            return self.mleduals[1]['qfunc'](c, w, r) if self.mleduals[1] is not None else 1\n",
    "           \n",
    "    def autotune(pre, target, taumax):\n",
    "        from scipy.optimize import root_scalar\n",
    "        \n",
    "        def f(tau):\n",
    "            from scipy.special import softmax\n",
    "            import numpy as np\n",
    "\n",
    "            soft = softmax(tau * pre, axis=1)\n",
    "            minsoft = np.min(soft, axis=1)\n",
    "            \n",
    "            return np.mean(minsoft) - target\n",
    "            \n",
    "        fmax = f(taumax)\n",
    "        \n",
    "        if fmax >= 0:\n",
    "            return taumax\n",
    "        \n",
    "        taumin = 0\n",
    "        fmin = f(taumin)\n",
    "        \n",
    "        assert fmin > 0, { 'fmin': fmin, 'target': target }\n",
    "        \n",
    "        root = root_scalar(f, bracket=(taumin, taumax))\n",
    "        \n",
    "        assert root.converged, root\n",
    "                \n",
    "        return root.root\n",
    "            \n",
    "    def flass():\n",
    "        from scipy.sparse import vstack\n",
    "        from scipy.special import softmax\n",
    "        from sklearn.datasets import fetch_20newsgroups\n",
    "        from sklearn.feature_extraction.text import HashingVectorizer\n",
    "        from sklearn.linear_model import SGDClassifier\n",
    "        from sklearn.metrics import accuracy_score\n",
    "        from math import exp\n",
    "        import numpy as np\n",
    "        \n",
    "        train = fetch_20newsgroups(subset='train', remove=('headers', 'footers', 'quotes'))\n",
    "        classes = np.unique(train.target)\n",
    "        test = fetch_20newsgroups(subset='test', remove=('headers', 'footers', 'quotes'))\n",
    "                \n",
    "        subblocksize = 32\n",
    "        delay = 8\n",
    "        blocksize = delay * subblocksize\n",
    "        for lr, taumax, target, what, gamma in ( (x, 8, y, z, g) for x in (0.5,)\n",
    "                                                                 for z in ('mle', 'ub', 'lb') \n",
    "                                                                 for g in (0.9999, 0.999,) \n",
    "                                                                 for y in (0.025, 0.05, 0.1, 0.2, ) ):\n",
    "            \n",
    "            ngram = 2\n",
    "            norm = 'l2'\n",
    "            nfeat = 20\n",
    "            \n",
    "            print(\"*** lr = {} taumax = {} target = {} what = {} gamma = {} ***\".format(lr, taumax, target, what, gamma), flush=True)\n",
    "            print('{:8.8s}\\t{:8.8s} [{:10.10s}]\\t{:8.8s} [{:10.10s}]\\t{:8.8s} [{:10.10s}]\\t{:8.8s} [{:10.10s}]'.format(\n",
    "                'n', 'eff n', 'since last', 'av tau', 'since last', 'emp loss', 'since last', 'log pv', 'since last')\n",
    "            )\n",
    "            \n",
    "            try:\n",
    "                target /= len(classes)\n",
    "\n",
    "                vectorizer = HashingVectorizer(n_features = 1 << nfeat, norm=norm, ngram_range=(1, ngram), alternate_sign=True)\n",
    "                docs = vectorizer.transform(train.data)\n",
    "                testdocs = vectorizer.transform(test.data)\n",
    "            \n",
    "                classweights = { k: lr for k in classes }\n",
    "                cls = SGDClassifier(loss='log', class_weight=classweights, shuffle=False)\n",
    "                \n",
    "                loss = OnlineDRO.EasyAcc()\n",
    "                sincelast = OnlineDRO.EasyAcc()\n",
    "                logpv = OnlineDRO.EasyAcc()\n",
    "                logpvsl = OnlineDRO.EasyAcc()\n",
    "                effn = OnlineDRO.EasyAcc()\n",
    "                effnsl = OnlineDRO.EasyAcc()\n",
    "                avtau = OnlineDRO.EasyAcc()\n",
    "                avtausl = OnlineDRO.EasyAcc()\n",
    "\n",
    "                loggerrand = np.random.RandomState(seed=2112)\n",
    "                logchoices = [None]*len(train.data)\n",
    "                pchoices = [None]*len(train.data)\n",
    "\n",
    "                ocrl = OnlineDRO.OnlineCressieReadLB(alpha=0.05, \n",
    "                                                     gamma=gamma,\n",
    "                                                     wmin=0,\n",
    "                                                     wmax=exp(taumax) + len(classes) - 1\n",
    "                                                    )\n",
    "                qfunc = ocrl.qmle if what == 'mle' else ocrl.qlb if what == 'lb' else ocrl.qub\n",
    "\n",
    "                for pno in range(1):\n",
    "                    order = np.random.RandomState(seed=42+pno).permutation(len(train.data))\n",
    "                    for n, ind in enumerate(zip(*(iter(order),)*blocksize)):\n",
    "                        v = vstack([ docs[z] for z in ind ])\n",
    "                        if n == 0 and pno == 0:\n",
    "                            pred = np.zeros(blocksize)\n",
    "                            for i, z in enumerate(ind):\n",
    "                                if logchoices[z] is None:\n",
    "                                    choice = loggerrand.choice(a=classes, size=1)\n",
    "                                    logchoices[z] = choice[0]\n",
    "                                    pchoices[z] = 1.0 / len(classes)\n",
    "                        else:\n",
    "                            predlogp = cls.predict_proba(v)\n",
    "                            tau = OnlineDRO.autotune(predlogp, target, taumax)\n",
    "                            avtau += tau\n",
    "                            avtausl += tau\n",
    "                            soft = softmax(tau * predlogp, axis=1)\n",
    "\n",
    "                            for i, z in enumerate(ind):\n",
    "                                if logchoices[z] is None:\n",
    "                                    choice = loggerrand.choice(a=classes, p=soft[i,:], size=1)\n",
    "                                    logchoices[z] = choice[0]\n",
    "                                    pchoices[z] = soft[i, choice[0]]\n",
    "\n",
    "                            pred = cls.predict(v)\n",
    "                            actual = [ train.target[z] for z in ind ]\n",
    "                            for i, (p, a) in enumerate(zip(pred, actual)):\n",
    "                                loss += 0 if p == a else 1\n",
    "                                sincelast += 0 if p == a else 1\n",
    "                                logpv += soft[i, a]\n",
    "                                logpvsl += soft[i, a]\n",
    "\n",
    "                            if (n & (n - 1) == 0): # and n & 0xAAAAAAAA == 0):\n",
    "                                print('{:<8d}\\t{:<8.3f} [{:<10.3f}]\\t{:<8.3f} [{:<10.3f}]\\t{:<8.3f} [{:<10.3f}]\\t{:<8.3f} [{:<10.3f}]'.format(\n",
    "                                             loss.n, effn.mean(), effnsl.mean(), avtau.mean(), avtausl.mean(), loss.mean(), sincelast.mean(), logpv.mean(), logpvsl.mean()),\n",
    "                                       flush=True)   \n",
    "                                sincelast = OnlineDRO.EasyAcc()\n",
    "                                logpvsl = OnlineDRO.EasyAcc()\n",
    "                                effnsl = OnlineDRO.EasyAcc()\n",
    "                                avtausl = OnlineDRO.EasyAcc()\n",
    "\n",
    "                        for i, z in enumerate(ind):\n",
    "                            r = 1 if logchoices[z] == train.target[z] else 0\n",
    "                            w = 1 / pchoices[z] if pred[i] == logchoices[z] else 0\n",
    "                            ocrl.update(1, w, r) \n",
    "\n",
    "                        for d in range(delay):\n",
    "                            y = np.array([ logchoices[z] \n",
    "                                           for i, z in enumerate(ind) \n",
    "                                           if (d-1)*subblocksize <= i and i < d*subblocksize\n",
    "                                           if logchoices[z] == train.target[z] ])\n",
    "\n",
    "                            if n == 0 and pno == 0:\n",
    "                                sampweight = np.array([ lr \n",
    "                                                        for i, z in enumerate(ind) \n",
    "                                                        if (d-1)*subblocksize <= i and i < d*subblocksize\n",
    "                                                        if logchoices[z] == train.target[z] ])\n",
    "                            else:\n",
    "                                sampweight = np.array([ lr * w * ocrl.n \n",
    "                                                           * max(0, qfunc(1, w, 1))\n",
    "                                                        for i, z in enumerate(ind) \n",
    "                                                        if (d-1)*subblocksize <= i and i < d*subblocksize\n",
    "                                                        if logchoices[z] == train.target[z]\n",
    "                                                        for w in (1 / pchoices[z],)\n",
    "                                                      ])\n",
    "\n",
    "                            effn += sampweight.sum() / (lr * subblocksize)\n",
    "                            effnsl += sampweight.sum() / (lr * subblocksize)\n",
    "\n",
    "                            if np.any(y):\n",
    "                                x = vstack([ v[i] \n",
    "                                             for i, z in enumerate(ind)\n",
    "                                             if (d-1)*subblocksize <= i and i < d*subblocksize\n",
    "                                             if logchoices[z] == train.target[z] ])\n",
    "                                cls.partial_fit(x, y, classes=classes, sample_weight=sampweight)\n",
    "\n",
    "                    print('{:<8d}\\t{:<8.3f} [{:<10.3f}]\\t{:<8.3f} [{:<10.3f}]\\t{:<8.3f} [{:<10.3f}]\\t{:<8.3f} [{:<10.3f}]'.format(\n",
    "                                     loss.n, effn.mean(), effnsl.mean(), avtau.mean(), avtausl.mean(), loss.mean(), sincelast.mean(), logpv.mean(), logpvsl.mean()),\n",
    "                               flush=True)   \n",
    "                    sincelast = OnlineDRO.EasyAcc()\n",
    "                    logpvsl = OnlineDRO.EasyAcc()\n",
    "                    effnsl = OnlineDRO.EasyAcc()\n",
    "                    avtausl = OnlineDRO.EasyAcc()\n",
    "\n",
    "                    preds = cls.predict(testdocs)\n",
    "                    ascores = []\n",
    "                    for b in range(16):\n",
    "                        bootie = np.random.RandomState(90210+b).choice(len(test.target), replace=True, size=len(test.target))\n",
    "                        ascores.append(accuracy_score(y_true=test.target[bootie], y_pred=preds[bootie]))\n",
    "\n",
    "                    print(\"test accuracy: {}\".format(np.quantile(ascores, [0.05, 0.5, 0.95])))\n",
    "            except KeyboardInterrupt:\n",
    "                raise\n",
    "            except:\n",
    "                pass\n",
    "    \n",
    "OnlineDRO.flass()      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hidden": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}