OCT.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "creative-graphic",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "D:\\Programming\\Python\\Python3.6\\python.exe\n",
      "3.6.4 (v3.6.4:d48eceb, Dec 19 2017, 06:54:40) [MSC v.1900 64 bit (AMD64)]\n",
      "sys.version_info(major=3, minor=6, micro=4, releaselevel='final', serial=0)\n"
     ]
    }
   ],
   "source": [
    "from pyscipopt import Model\n",
    "from pyscipopt import quicksum\n",
    "import sys\n",
    "import math\n",
    "import numpy as np\n",
    "from collections import Counter\n",
    "print(sys.executable)\n",
    "print(sys.version)\n",
    "print(sys.version_info)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "standard-immigration",
   "metadata": {},
   "outputs": [],
   "source": [
    "model = Model(\"OCT Example\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "gorgeous-black",
   "metadata": {},
   "outputs": [],
   "source": [
    "X = [[0., 1.],\n",
    "     [1., 0.],\n",
    "     [0., 0.]]\n",
    "\n",
    "y = [0, 1, 2]\n",
    "\n",
    "# max_depth\n",
    "D = 2\n",
    "# num of nodes\n",
    "T = int(math.pow(2, D + 1) - 1)\n",
    "# num of branch nodes\n",
    "TB = int(T/2)\n",
    "# num of features\n",
    "P = len(X[0])\n",
    "# num of training points\n",
    "n = len(X)\n",
    "# number of labels\n",
    "y_set = set(y)\n",
    "K = len(y_set)\n",
    "# hyper parameter-minimum leaf size\n",
    "Nmin = 1\n",
    "# complexity parameter \n",
    "alpha = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "foreign-partner",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_parents(t):\n",
    "    a = t\n",
    "    parents = []\n",
    "    while a != 1:\n",
    "        a = int(a/2)\n",
    "        parents.append(a)\n",
    "    return parents\n",
    "\n",
    "def get_AL(t):\n",
    "    a = t\n",
    "    ALs = []\n",
    "    while a != 1:\n",
    "        temp = int(a/2)\n",
    "        if temp == a/2:\n",
    "            ALs.append(int(temp))   \n",
    "        a = int(a/2)\n",
    "    return ALs\n",
    "\n",
    "def get_AR(t):\n",
    "    a = t\n",
    "    ARs = []\n",
    "    while a != 1:\n",
    "        temp = int(a/2)\n",
    "        if temp != a/2:\n",
    "            ARs.append(int(temp))   \n",
    "        a = int(a/2)\n",
    "    return ARs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "creative-facility",
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = np.zeros((n, K))\n",
    "for i in range(1, n + 1):\n",
    "    for k in range(1, K + 1):\n",
    "        Y[i - 1, k - 1] = 1 if y[i - 1] == k - 1 else -1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "center-scheduling",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.0, 1.0]\n"
     ]
    }
   ],
   "source": [
    "# p * 1\n",
    "eps = [0] * P\n",
    "X_np = np.array(X)\n",
    "\n",
    "for j in range(1, P + 1):\n",
    "    sorted_features = np.sort(X_np[:, j - 1])\n",
    "    min_diff = float(\"inf\")\n",
    "    for i in range(2, len(sorted_features)):\n",
    "        if sorted_features[i] != sorted_features[i-1]:\n",
    "            min_diff = min(min_diff, sorted_features[i] - sorted_features[i-1])\n",
    "    eps[j - 1] = min_diff\n",
    "print(eps)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "serious-fence",
   "metadata": {},
   "outputs": [],
   "source": [
    "counter = Counter(y)\n",
    "L_hat = n - counter.most_common(1)[0][1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "final-screening",
   "metadata": {},
   "outputs": [],
   "source": [
    "a = {}\n",
    "for j in range(1, P + 1):\n",
    "    for t in range(1, TB + 1):\n",
    "        a[j,t] = model.addVar(vtype=\"B\", name=\"a(%s,%s)\" % (j, t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "hidden-barbados",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{1: b(1), 2: b(2), 3: b(3)}\n"
     ]
    }
   ],
   "source": [
    "# b = np.zeros(TB)\n",
    "b = {}\n",
    "for t in range(1, TB + 1):\n",
    "    b[t] = model.addVar(vtype=\"C\", lb=0.0, name=\"b(%s)\" % (t))\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "alive-distribution",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{1: d(1), 2: d(2), 3: d(3)}\n"
     ]
    }
   ],
   "source": [
    "d = {}\n",
    "for t in range(1, TB + 1):\n",
    "    d[t] = model.addVar(vtype=\"B\", name=\"d(%s)\" % (t))\n",
    "print(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "lucky-reynolds",
   "metadata": {},
   "outputs": [],
   "source": [
    "z = {}\n",
    "for i in range(1, n + 1):\n",
    "    for t in range(TB + 1, T + 1):\n",
    "        z[i,t] = model.addVar(vtype=\"B\", name=\"z(%s,%s)\" % (i,t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "normal-probe",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{4: l(4), 5: l(5), 6: l(6), 7: l(7)}\n"
     ]
    }
   ],
   "source": [
    "l = {}\n",
    "for t in range(TB + 1, T + 1):\n",
    "    l[t] = model.addVar(vtype=\"B\", name=\"l(%s)\" % (t))\n",
    "print(l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "devoted-invite",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = {}\n",
    "for k in range(1, K + 1):\n",
    "    for t in range(TB + 1, T + 1):\n",
    "        c[k,t] = model.addVar(vtype=\"B\", name=\"c(%s,%s)\" % (k,t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "sufficient-protein",
   "metadata": {},
   "outputs": [],
   "source": [
    "Nkt = {}\n",
    "for k in range(1, K + 1):\n",
    "    for t in range(TB + 1, T + 1):\n",
    "        Nkt[k, t] = model.addVar(vtype=\"I\", name=\"Nkt(%s,%s)\" % (k, t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "south-liechtenstein",
   "metadata": {},
   "outputs": [],
   "source": [
    "Nt = {}\n",
    "for t in range(TB + 1, T + 1):\n",
    "    Nt[t] = model.addVar(vtype=\"I\", name=\"Nt(%s)\" % (t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "dominican-syndrome",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = {}\n",
    "for t in range(TB + 1, T + 1):\n",
    "    L[t] = model.addVar(vtype=\"I\", lb=0, name=\"L(%s)\" % (t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "whole-determination",
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in range(2, TB + 1):\n",
    "    parent_index = get_parents(t)[0]\n",
    "#     print(t, d_parent, parent_index)\n",
    "    model.addCons(d[t] <= d[parent_index], name=\"DCons(%s)\" % t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "equal-amsterdam",
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in range(1, TB + 1):\n",
    "    model.addCons(b[t] <= d[t], name=\"BCons(%s)\" % t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "million-olympus",
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in range(1, TB + 1):\n",
    "    model.addCons(quicksum(a[j,t] for j in range(1, P + 1)) == d[t], name=\"aCons(%s)\" % t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "explicit-collective",
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in range(TB + 1, T + 1):\n",
    "    model.addCons(quicksum(z[i,t] for i in range(1, n + 1)) >= Nmin * l[t] , name=\"zCons(%s)\" % t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "actual-dairy",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(1, n + 1):\n",
    "    for t in range(TB + 1, T + 1):\n",
    "        model.addCons(z[i, t] <= l[t], name=\"zCons2(%s, %s)\" % (i, t))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "iraqi-address",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(1, n + 1):\n",
    "    model.addCons(quicksum(z[i,t] for t in range(TB + 1, T + 1)) == 1, name=\"zCons3(%s)\" % i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "moving-thomas",
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in range(TB + 1, T + 1):\n",
    "    model.addCons(quicksum(c[k,t] for k in range(1, K + 1)) == l[t], name=\"cCons(%s)\" % t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "visible-indiana",
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in range(TB + 1, T + 1):\n",
    "    Nt[t] = quicksum(z[i, t] for i in range(1, n + 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "naughty-nudist",
   "metadata": {},
   "outputs": [],
   "source": [
    "for k in range(1, K + 1):\n",
    "    for t in range(TB + 1, T + 1):\n",
    "        Nkt[k, t] = 0\n",
    "        for i in range(1, n + 1):\n",
    "             Nkt[k, t] += (1 + Y[i - 1, k - 1]) * z[i, t]\n",
    "        Nkt[k, t] /= 2 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "nasty-raise",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(1, n + 1):\n",
    "    for t in range(TB + 1, T + 1):\n",
    "        AL = get_AL(t)\n",
    "        for m in AL:\n",
    "            left = 0\n",
    "            for j in range(1, P + 1):\n",
    "                left +=  X_np[i-1, j-1] * a[j, m]\n",
    "            left += min(eps)\n",
    "            right = b[m] + (1 + max(eps))*(1 - z[i, t])\n",
    "            model.addCons(left <= right, name=\"AlCons(%s, %s, %s)\" % (i, t, m))        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "spoken-biology",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(1, n + 1):\n",
    "    for t in range(TB + 1, T + 1):\n",
    "        AR = get_AR(t)\n",
    "        for m in AR:\n",
    "            left = 0\n",
    "            for j in range(1, P + 1):\n",
    "                left += X_np[i-1, j-1] * a[j, m]\n",
    "            right = b[m] - (1 - z[i, t])\n",
    "            model.addCons(left >= right, name=\"ArCons(%s, %s, %s)\" % (i, t, m))   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "suburban-tender",
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in range(TB + 1, T + 1):\n",
    "    for k in range(1, K + 1):\n",
    "        model.addCons(L[t] <= Nt[t] - Nkt[k, t] + n * c[k, t], name=\"LCons1(%s, %s)\" % (t, k))        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "developed-enough",
   "metadata": {},
   "outputs": [],
   "source": [
    "for t in range(TB + 1, T + 1):\n",
    "    for k in range(1, K + 1):\n",
    "        model.addCons(L[t] >= Nt[t] - Nkt[k, t] - n * (1 - c[k, t]), name=\"LCons2(%s, %s)\" % (t, k))        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "identified-inspector",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.setObjective((1/L_hat) * quicksum(L[t] for t in range(TB + 1, T + 1)) + alpha * quicksum(d[t] for t in range(1, TB + 1)), \"minimize\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "direct-correlation",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.optimize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "wired-neighborhood",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal value: 0.0\n",
      "optimal\n"
     ]
    }
   ],
   "source": [
    "if model.getStatus() == \"optimal\":\n",
    "    print(\"Optimal value:\", model.getObjVal())\n",
    "print(model.getStatus())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "meaning-scientist",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Depth  1\n",
      "Branch 1, Split Attribute 2, Split Threshold 1.000\n",
      "\n",
      "Depth  2\n",
      "Branch 2, Split Attribute 1, Split Threshold 1.000\n",
      "Branch 3, Split Attribute 2, Split Threshold 0.000\n",
      "\n",
      "Depth  3\n",
      "Leaf 4, Class Label 3, Point Count 1, Loss 0\n",
      "Leaf 5, Class Label 2, Point Count 1, Loss 0\n",
      "Leaf 6, Class Label 1, Point Count 0, Loss 0\n",
      "Leaf 7, Class Label 1, Point Count 1, Loss 0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for i in range(0, D + 1):\n",
    "    print(\"Depth \", i + 1)\n",
    "    for t in range(2**i, 2**(i+1)):\n",
    "        if i != D: # Branch\n",
    "            for j in range(1, P+1):\n",
    "                if model.getVal(a[j, t]) == 1:\n",
    "                    print(\"Branch %d, Split Attribute %d, Split Threshold %.3f\" %(t, j, model.getVal(b[t])))\n",
    "                    break\n",
    "        else:   # Leaf\n",
    "            valuemax = -float('inf')\n",
    "            argmax = -10\n",
    "            for k in range(1, K + 1):\n",
    "                if model.getVal(Nkt[k, t]) > valuemax:\n",
    "                    valuemax = model.getVal(Nkt[k, t])\n",
    "                    argmax = k\n",
    "            print(\"Leaf %d, Class Label %d, Point Count %d, Loss %d\" %(t, argmax, model.getVal(Nt[t]), model.getVal(L[t])))\n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "earlier-crash",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv_new",
   "language": "python",
   "name": "venv_new"
  },
  "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}