Feature/random forest

experiments/random_forest/scratch.ipynb

@@ -0,0 +1,245 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "id": "97224de0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 36,
+   "id": "50446c72",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([[2, 0, 0, 2],\n",
+       "       [2, 0, 1, 1],\n",
+       "       [3, 3, 2, 3]])"
+      ]
+     },
+     "execution_count": 36,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "X = [[1,2],[2,2],[3,1],[4,3]]\n",
+    "n = 3\n",
+    "\n",
+    "rng = np.random.default_rng()\n",
+    "\n",
+    "samples = np.array([rng.choice(4, size = len(X), replace=True) for i in range(n)])\n",
+    "\n",
+    "samples"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 44,
+   "id": "a436cbea",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[[[3, 1], [1, 2], [1, 2], [3, 1]],\n",
+       " [[3, 1], [1, 2], [2, 2], [2, 2]],\n",
+       " [[4, 3], [4, 3], [3, 1], [4, 3]]]"
+      ]
+     },
+     "execution_count": 44,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "[[X[i] for i in s] for s in samples]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 41,
+   "id": "06238b06",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[[[3, 1], [1, 2], [1, 2], [3, 1]],\n",
+       " [[3, 1], [1, 2], [2, 2], [2, 2]],\n",
+       " [[4, 3], [4, 3], [3, 1], [4, 3]]]"
+      ]
+     },
+     "execution_count": 41,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "res = []\n",
+    "for s in samples:\n",
+    "    res.append([X[i] for i in s])\n",
+    "\n",
+    "res"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "7235fb9a",
+   "metadata": {},
+   "outputs": [
+    {
+     "ename": "NameError",
+     "evalue": "name 'np' is not defined",
+     "output_type": "error",
+     "traceback": [
+      "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
+      "\u001b[31mNameError\u001b[39m                                 Traceback (most recent call last)",
+      "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[1]\u001b[39m\u001b[32m, line 1\u001b[39m\n\u001b[32m----> \u001b[39m\u001b[32m1\u001b[39m ~\u001b[43mnp\u001b[49m.isin([[\u001b[32m0\u001b[39m],[\u001b[32m1\u001b[39m],[\u001b[32m2\u001b[39m],[\u001b[32m3\u001b[39m]], [[\u001b[32m0\u001b[39m],[\u001b[32m2\u001b[39m],[\u001b[32m2\u001b[39m],[\u001b[32m3\u001b[39m]])\n",
+      "\u001b[31mNameError\u001b[39m: name 'np' is not defined"
+     ]
+    }
+   ],
+   "source": [
+    "~np.isin([[0],[1],[2],[3]], [[0],[2],[2],[3]])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "id": "0292374c",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "array([0, 0, 1])"
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "from glassboxml.models.random_forest import RandomForest\n",
+    "\n",
+    "m = RandomForest(10)\n",
+    "\n",
+    "X = [[0], [1], [-1], [-2],  [2], [3], [4], [5]]\n",
+    "y = [0, 0, 0, 0, 1, 1, 1, 1]\n",
+    "\n",
+    "m.fit(X,y)\n",
+    "\n",
+    "m.accuracy\n",
+    "\n",
+    "m.predict([[-5], [-3], [2]])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 38,
+   "id": "66e886fb",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "np.float64(0.5873015873015873)"
+      ]
+     },
+     "execution_count": 38,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "X = [[1,9], [0,8], [-1, 5], [-15,4],\n",
+    "     [1,10], [1,19], [-150,150], [0,0],\n",
+    "     [2,19], [12,15], [17,10], [10,-150],\n",
+    "     [20,21], [29,28], [350,300], [1000,21]]\n",
+    "y = [0, 0, 0, 0,\n",
+    "     1,1,1,1,\n",
+    "     2,2,2,2,\n",
+    "     3,3,3,3]\n",
+    "\n",
+    "model = RandomForest(10, max_features=2)\n",
+    "model.fit(X, y)\n",
+    "\n",
+    "model.accuracy"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 13,
+   "id": "0126cb68",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from glassboxml.models.decision_tree import DecisionTree\n",
+    "\n",
+    "m = DecisionTree(10)\n",
+    "\n",
+    "m.fit(X,y)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 39,
+   "id": "d976ae48",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "{0, 1, 2, 3}"
+      ]
+     },
+     "execution_count": 39,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "set([0,1,1,2,3])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9f00bded",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": ".venv (3.11.7.final.0)",
+   "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.11.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

glassboxml/models/decision_tree.py

@@ -17,14 +17,20 @@ def __init__(self, best_feature, best_threshold, left_child, right_child):
         self.right_child = right_child
 
 class DecisionTree:
-    def __init__(self, max_depth):
+    def __init__(self, max_depth, max_features=None):
         self.max_depth = max_depth
+        # None = consider every feature at every split (current/unchanged behavior).
+        # Otherwise, the number of features to randomly sample per split (e.g. sqrt(n_features)).
+        self.max_features = max_features
         self.root = None
 
     def fit(self, X, y):
         X = np.array(X)
         y = np.array(y)
 
+        if self.max_features == None:
+            self.max_features = X.shape[1]
+
         self.root = self.build_tree(X,y, 0)
 
     def build_tree(self, X, y, depth):
@@ -73,6 +79,11 @@ def calculate_gini(self, y):
         probs = counts / counts.sum()
         return 1 - np.sum(probs ** 2)
 
+    def _sample_feature_indices(self, n_features):
+        rng = np.random.default_rng()
+
+        return rng.choice(n_features, size = self.max_features, replace=False)
+
     def find_best_split(self, X, y):
         n_samples, n_features = X.shape
         best_IG = 0
@@ -83,7 +94,9 @@ def find_best_split(self, X, y):
 
         gini_parent = self.calculate_gini(y)
 
-        for f in range(n_features):
+        current_features = self._sample_feature_indices(n_features)
+
+        for f in current_features:
             values = np.unique(X[:, f])
             thresholds = (values[:-1] + values[1:])/2
             for t in thresholds:

glassboxml/models/random_forest.py

@@ -0,0 +1,66 @@
+import numpy as np
+from concurrent.futures import ProcessPoolExecutor
+from itertools import repeat
+from scipy.stats import mode
+from  glassboxml.models.decision_tree import DecisionTree
+from glassboxml.metrics.classification import accuracy
+
+def _fit_one_tree(train_X, train_y, X_oob, y_oob, max_features):
+    # Module-level (not a method) so it can be pickled and sent to a
+    # worker process - each tree only needs its own bootstrap sample and
+    # OOB set, so this can run fully independently of the other trees
+    model = DecisionTree(100, max_features)
+    model.fit(train_X, train_y)
+    tree_accuracy = accuracy(y_oob, model.predict(X_oob))
+    return model, tree_accuracy
+
+class RandomForest:
+    def __init__(self, n, max_features=None):
+        self.n = n
+        self.max_features = max_features
+
+    def fit(self, X, y):
+        self.X_train = np.array(X)
+        self.y_train = np.array(y)
+
+        if self.max_features == None:
+            self.max_features = np.int16(np.floor(np.sqrt(len(X[0]))))
+
+        #Bootstrapping
+        rng = np.random.default_rng()
+        # Create n samples for each tree where each sample is of size - len(X). 
+        # Each point is selected with replacement.
+        indice_samples = rng.choice(len(X), size = (self.n, len(X)), replace=True)
+
+        # Create train dataset from sampled indices
+        train_X = self.X_train[indice_samples]
+        train_y = self.y_train[indice_samples]
+
+        #Out-of-Bag(OOB) test computations
+        orig_indices = np.arange(len(X))
+        test_OOB = [np.setdiff1d(orig_indices, s) for s in indice_samples]
+
+        # Training n trees on n samples - each tree is fit independently of
+        # the others, so hand the n fits out across worker processes instead
+        # of fitting them one at a time on a single core
+        X_oob = [self.X_train[oob] for oob in test_OOB]
+        y_oob = [self.y_train[oob] for oob in test_OOB]
+
+        with ProcessPoolExecutor() as executor:
+            results = list(executor.map(
+                _fit_one_tree, train_X, train_y, X_oob, y_oob, repeat(self.max_features, self.n)
+            ))
+
+        self.models = [model for model, _ in results]
+        self.accuracy = np.mean([tree_accuracy for _, tree_accuracy in results])
+
+    def predict(self, X):
+        # Each tree's own vectorized predict() returns every sample's
+        # prediction in one call - stacking these gives a (n_trees, n_samples)
+        # matrix where column j holds every tree's vote for sample j
+        t_pred = np.array([t.predict(X) for t in self.models])
+
+        # Majority vote per sample = the mode down each column (axis=0)
+        y_pred, _ = mode(t_pred, axis=0)
+
+        return y_pred

pyproject.toml

@@ -3,7 +3,8 @@ name = "glassboxml"
 version = "0.1.0"
 dependencies = [
     "numpy",
-    "matplotlib"
+    "matplotlib",
+    "scipy"
 ]
 
 [project.optional-dependencies]