Update Algorithms

Author: carefree0910

NN/Basic/Layers.py

@@ -17,39 +17,47 @@ def conv_bp(n, n_filters, out_h, out_w, dx_padded,
             for j in range(out_h):
                 for k in range(out_w):
                     for h in range(dx_padded.shape[1]):
-                        dx_padded[i, h, j * sd:filter_height + j * sd, k * sd:filter_width + k * sd] += (
-                            inner_weight[f][h] * delta[i, f, j, k]
-                        )
+                        jsd, ksd = j * sd, k * sd
+                        for p in range(filter_height):
+                            for q in range(filter_width):
+                                dx_padded[i, h, jsd+p, ksd+q] += (
+                                    inner_weight[f][h][p][q] * delta[i, f, j, k]
+                                )
 
 
 @numba.jit([
     "void(int64, int64, int64, int64, float32[:,:,:,:], float32[:,:,:,:],"
-    "int64, int64, int64)"
+    "int64, int64, int64, int32[:,:,:,:,:])"
 ], nopython=True)
 def max_pool(n, n_channels, out_h, out_w, x, out,
-             pool_height, pool_width, sd):
+             pool_height, pool_width, sd, pos_cache):
     for i in range(n):
         for j in range(n_channels):
             for k in range(out_h):
                 for l in range(out_w):
-                    window = x[i, j, k * sd:pool_height + k * sd, l * sd:pool_width + l * sd]
-                    out[i, j, k, l] = np.max(window)
+                    ksd, lsd = k * sd, l * sd
+                    _max = x[i, j, ksd, lsd]
+                    pos = (0, 0)
+                    for p in range(pool_height):
+                        for q in range(pool_width):
+                            if x[i, j, ksd+p, lsd+q] > _max:
+                                _max = x[i, j, ksd+p, lsd+q]
+                                pos = (p, q)
+                    pos_cache[i, j, k, l] = pos
+                    out[i, j, k, l] = _max
 
 
 @numba.jit([
-    "void(int64, int64, int64, int64, float32[:,:,:,:],"
-    "int64, int64, int64, float32[:,:,:,:], float32[:,:,:,:])"
+    "void(int64, int64, int64, int64, int64, float32[:,:,:,:], float32[:,:,:,:], int32[:,:,:,:,:])"
 ], nopython=True)
-def max_pool_bp(n, n_channels, out_h, out_w, x_cache,
-                pool_height, pool_width, sd, dx, delta):
+def max_pool_bp(n, n_channels, out_h, out_w, sd, dx, delta, pos_cache):
     for i in range(n):
         for j in range(n_channels):
             for k in range(out_h):
                 for l in range(out_w):
-                    window = x_cache[i, j, k * sd:pool_height + k * sd, l * sd:pool_width + l * sd]
-                    dx[i, j, k * sd:pool_height + k * sd, l * sd:pool_width + l * sd] = (
-                        (window == np.max(window)) * delta[i, j, k, l]
-                    )
+                    ksd, lsd = k * sd, l * sd
+                    pos = pos_cache[i, j, k, l]
+                    dx[i, j, ksd+pos[0], lsd+pos[1]] = delta[i, j, k, l]
 
 
 # Abstract Layers
@@ -519,11 +527,13 @@ def _activate(self, x, *args):
             self._pool_cache["method"] = "reshape"
         else:
             out = np.zeros((n, n_channels, self.out_h, self.out_w), dtype=np.float32)
+            pos_cache = np.zeros((n, n_channels, self.out_h, self.out_w, 2), dtype=np.int32)
             max_pool(
                 n, n_channels, self.out_h, self.out_w, x, out,
-                pool_height, pool_width, sd
+                pool_height, pool_width, sd, pos_cache
             )
             self._pool_cache["method"] = "original"
+            self._pool_cache["pos_cache"] = pos_cache
         return out
 
     def _derivative(self, y, *args):
@@ -549,10 +559,8 @@ def _derivative(self, y, *args):
             sd = self._stride
             dx = np.zeros_like(self.x_cache)
             n, n_channels, *_ = self.x_cache.shape
-            _, pool_height, pool_width = self._shape[1]
             max_pool_bp(
-                n, n_channels, self.out_h, self.out_w, self.x_cache,
-                pool_height, pool_width, sd, dx, delta
+                n, n_channels, self.out_h, self.out_w, sd, dx, delta, self._pool_cache["pos_cache"]
             )
         else:
             raise LayerError("Undefined pooling method '{}' found".format(method))

Notebooks/numba/Basic(zh-cn).ipynb

@@ -75,7 +75,7 @@
     "+ 注意：\n",
     "    + `numba`不支持 list comprehension，详情可参见[这里](https://github.com/numba/numba/issues/504)\n",
     "    + `jit`会在某种程度上“预编译”你的代码，这意味着它会在某种程度上固定住各个变量的数据类型；所以在`jit`下定义数组时，如果想要使用的是`float`数组的话，就不能像上述`wrong_add`里那样用`[0] * len(x)`定义、而应该在`0`后面加一个小数点：`[0.] * len(x)`\n",
-    "    + `jit`能够加速的不限于`for`，但一般而言加速`for`会比较常见、效果也比较显著。我在我实现的`numpy`版本的卷积神经网络（`CNN`）中用了`jit`后、可以把代码加速 **20** 倍左右。具体代码可以参见[这里](https://github.com/carefree0910/MachineLearning/blob/master/NN/Basic/Layers.py#L9)，不过如果不想看源代码的话，可以参见[CNN(zh-cn).ipynb][1]，我在其中做了一些相应的、比较简单的实验\n",
+    "    + `jit`能够加速的不限于`for`，但一般而言加速`for`会比较常见、效果也比较显著。我在我实现的`numpy`版本的卷积神经网络（`CNN`）中用了`jit`后、可以把代码加速 **60 多倍**。具体代码可以参见[这里](https://github.com/carefree0910/MachineLearning/blob/master/NN/Basic/Layers.py#L9)，不过如果不想看源代码的话，可以参见[CNN(zh-cn).ipynb][1]，我在其中做了一些相应的、比较简单的实验\n",
     "\n",
     "[1]: https://github.com/carefree0910/MachineLearning/blob/master/Notebooks/numba/CNN(zh-cn).ipynb"
    ]

Notebooks/numba/CNN(zh-cn).ipynb

@@ -9,12 +9,12 @@
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "0.00113585\n",
-      "0.000733545\n",
-      "0.00112681\n",
-      "3.63 s ± 194 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
-      "300 ms ± 20.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
-      "8.69 ms ± 223 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
+      "0.00116325\n",
+      "0.000750613\n",
+      "0.00115522\n",
+      "3.32 s ± 115 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
+      "300 ms ± 19.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
+      "8.34 ms ± 171 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
      ]
     }
    ],
@@ -29,7 +29,6 @@
     "                window = x[i, ..., j:j+filter_height, p:p+filter_width]\n",
     "                for q in range(n_filters):\n",
     "                    rs[i, q, j, p] += np.sum(w[q] * window)\n",
-    "    return rs\n",
     "\n",
     "@nb.jit(nopython=True)\n",
     "def jit_conv_kernel(x, w, rs, n, n_channels, height, width, n_filters, filter_height, filter_width, out_h, out_w):\n",
@@ -84,18 +83,49 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "592 ms ± 25.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
-      "8.5 ms ± 150 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
+      "288 ms ± 33.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
+      "64.8 ms ± 815 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
+      "7.76 ms ± 104 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
      ]
     }
    ],
+   "source": [
+    "@nb.jit(nopython=True)\n",
+    "def jit_conv_kernel2(x, w, rs, n, n_channels, height, width, n_filters, filter_height, filter_width, out_h, out_w):\n",
+    "    for i in range(n):\n",
+    "        for j in range(out_h):\n",
+    "            for p in range(out_w):\n",
+    "                for q in range(n_filters):\n",
+    "                    for r in range(n_channels):\n",
+    "                        for s in range(filter_height):\n",
+    "                            for t in range(filter_width):\n",
+    "                                rs[i, q, j, p] += x[i, r, j+s, p+t] * w[q, r, s, t]\n",
+    "                                \n",
+    "assert np.allclose(conv(x, w, jit_conv_kernel, args), conv(x, w, jit_conv_kernel, args))\n",
+    "%timeit conv(x, w, jit_conv_kernel, args)\n",
+    "%timeit conv(x, w, jit_conv_kernel2, args)\n",
+    "%timeit cs231n_conv(x, w, args)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "+ 可以看到，使用`jit`和使用纯`numpy`进行编程的很大一点不同就是，不要畏惧用`for`；事实上一般来说，代码“长得越像 C”、速度就会越快"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
    "source": [
     "def max_pool_kernel(x, rs, *args):\n",
     "    n, n_channels, pool_height, pool_width, out_h, out_w = args\n",
@@ -115,6 +145,21 @@
     "                for q in range(out_w):\n",
     "                    window = x[i, j, p:p+pool_height, q:q+pool_width]\n",
     "                    rs[i, j, p, q] += np.max(window)\n",
+    "                    \n",
+    "@nb.jit(nopython=True)\n",
+    "def jit_max_pool_kernel2(x, rs, *args):\n",
+    "    n, n_channels, pool_height, pool_width, out_h, out_w = args\n",
+    "    for i in range(n):\n",
+    "        for j in range(n_channels):\n",
+    "            for p in range(out_h):\n",
+    "                for q in range(out_w):\n",
+    "                    _max = x[i, j, p, q]\n",
+    "                    for r in range(pool_height):\n",
+    "                        for s in range(pool_width):\n",
+    "                            _tmp = x[i, j, p+r, q+s]\n",
+    "                            if _tmp > _max:\n",
+    "                                _max = _tmp\n",
+    "                    rs[i, j, p, q] += _max\n",
     "\n",
     "def max_pool(x, kernel, args):\n",
     "    n, n_channels = args[:2]\n",
@@ -130,8 +175,20 @@
     "args = (n, n_channels, pool_height, pool_width, out_h, out_w)\n",
     "\n",
     "assert np.allclose(max_pool(x, max_pool_kernel, args), max_pool(x, jit_max_pool_kernel, args))\n",
+    "assert np.allclose(max_pool(x, jit_max_pool_kernel, args), max_pool(x, jit_max_pool_kernel2, args))\n",
     "%timeit max_pool(x, max_pool_kernel, args)\n",
-    "%timeit max_pool(x, jit_max_pool_kernel, args)"
+    "%timeit max_pool(x, jit_max_pool_kernel, args)\n",
+    "%timeit max_pool(x, jit_max_pool_kernel2, args)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%load_ext line_profiler\n",
+    "%lprun -f max_pool(x, jit_max_pool_kernel2, args)"
    ]
   }
  ],
@@ -156,5 +213,5 @@
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 1
 }

Util/Bases.py

@@ -2,10 +2,13 @@
 import cv2
 import time
 import math
+import ctypes
+import multiprocessing
 import numpy as np
 import tensorflow as tf
 import matplotlib.pyplot as plt
 from PIL import Image
+from multiprocessing import Pool
 from mpl_toolkits.mplot3d import Axes3D
 
 from Util.Util import VisUtil
@@ -214,63 +217,57 @@ def f1_score(y, y_pred):
 
     # noinspection PyUnusedLocal
     @staticmethod
-    def _multi_clf(x, clfs, task, kwargs, stack=np.vstack):
-        # n_cores = kwargs.get("n_cores", 2)
-        # n_cores = multiprocessing.cpu_count() if n_cores <= 0 else n_cores
-        # if n_cores == 1:
-        #     matrix = np.array([clf.predict(x, n_cores=1) for clf in clfs], dtype=np.float32).T
-        # else:
-        #     pool = Pool(max_workers=n_cores)
-        #     batch_size = int(len(clfs) / n_cores)
-        #     batch_clfs, cursor = [], 0
-        #     for i in range(n_cores):
-        #         if i == n_cores - 1:
-        #             batch_clfs.append(clfs[cursor:])
-        #         else:
-        #             batch_clfs.append(clfs[cursor:cursor + batch_size])
-        #         cursor += batch_size
-        #     x_size = np.prod(x.shape)  # type: int
-        #     shared_base = multiprocessing.Array(ctypes.c_float, int(x_size))
-        #     shared_matrix = np.ctypeslib.as_array(shared_base.get_obj()).reshape(x.shape)
-        #     shared_matrix[:] = x
-        #     del x, clfs, shared_base
-        #     matrix = stack(
-        #         pool.map(task, ((shared_matrix, clfs, 1) for clfs in batch_clfs))
-        #     ).T.astype(np.float32)
-        # return matrix
-        return np.array([clf.predict(x) for clf in clfs], dtype=np.float32).T
+    def _multi_clf(x, clfs, task, kwargs, stack=np.vstack, target="single"):
+        if target != "parallel":
+            return np.array([clf.predict(x) for clf in clfs], dtype=np.float32).T
+        n_cores = kwargs.get("n_cores", 2)
+        n_cores = multiprocessing.cpu_count() if n_cores <= 0 else n_cores
+        if n_cores == 1:
+            matrix = np.array([clf.predict(x, n_cores=1) for clf in clfs], dtype=np.float32).T
+        else:
+            pool = Pool(processes=n_cores)
+            batch_size = int(len(clfs) / n_cores)
+            clfs = [clfs[i*batch_size:(i+1)*batch_size] for i in range(n_cores)]
+            x_size = np.prod(x.shape)  # type: int
+            shared_base = multiprocessing.Array(ctypes.c_float, int(x_size))
+            shared_matrix = np.ctypeslib.as_array(shared_base.get_obj()).reshape(x.shape)
+            shared_matrix[:] = x
+            matrix = stack(
+                pool.map(task, ((shared_matrix, clfs, n_cores) for clfs in clfs))
+            ).T.astype(np.float32)
+        return matrix
 
     # noinspection PyUnusedLocal
-    def _multi_data(self, x, task, kwargs, stack=np.hstack):
-        # n_cores = kwargs.get("n_cores", 2)
-        # n_cores = multiprocessing.cpu_count() if n_cores <= 0 else n_cores
-        # if n_cores == 1:
-        #     matrix = task((x, self, n_cores))
-        # else:
-        #     pool = Pool(max_workers=n_cores)
-        #     batch_size = int(len(x) / n_cores)
-        #     batch_base, batch_data, cursor = [], [], 0
-        #     x_dim = x.shape[1]
-        #     for i in range(n_cores):
-        #         if i == n_cores - 1:
-        #             batch_data.append(x[cursor:])
-        #             batch_base.append(multiprocessing.Array(ctypes.c_float, (len(x) - cursor) * x_dim))
-        #         else:
-        #             batch_data.append(x[cursor:cursor + batch_size])
-        #             batch_base.append(multiprocessing.Array(ctypes.c_float, batch_size * x_dim))
-        #         cursor += batch_size
-        #     shared_arrays = [
-        #         np.ctypeslib.as_array(shared_base.get_obj()).reshape(-1, x_dim)
-        #         for shared_base in batch_base
-        #     ]
-        #     for i, data in enumerate(batch_data):
-        #         shared_arrays[i][:] = data
-        #     del x, batch_base, batch_data
-        #     matrix = stack(
-        #         pool.map(task, ((x, self, n_cores) for x in shared_arrays))
-        #     )
-        # return matrix.astype(np.float32)
-        return task((x, self, 1))
+    def _multi_data(self, x, task, kwargs, stack=np.hstack, target="single"):
+        if target != "parallel":
+            return task((x, self, 1))
+        n_cores = kwargs.get("n_cores", 2)
+        n_cores = multiprocessing.cpu_count() if n_cores <= 0 else n_cores
+        if n_cores == 1:
+            matrix = task((x, self, n_cores))
+        else:
+            pool = Pool(processes=n_cores)
+            batch_size = int(len(x) / n_cores)
+            batch_base, batch_data, cursor = [], [], 0
+            x_dim = x.shape[1]
+            for i in range(n_cores):
+                if i == n_cores - 1:
+                    batch_data.append(x[cursor:])
+                    batch_base.append(multiprocessing.Array(ctypes.c_float, (len(x) - cursor) * x_dim))
+                else:
+                    batch_data.append(x[cursor:cursor + batch_size])
+                    batch_base.append(multiprocessing.Array(ctypes.c_float, batch_size * x_dim))
+                cursor += batch_size
+            shared_arrays = [
+                np.ctypeslib.as_array(shared_base.get_obj()).reshape(-1, x_dim)
+                for shared_base in batch_base
+            ]
+            for i, data in enumerate(batch_data):
+                shared_arrays[i][:] = data
+            matrix = stack(
+                pool.map(task, ((x, self, n_cores) for x in shared_arrays))
+            )
+        return matrix.astype(np.float32)
 
     @staticmethod
     def _get_train_repeat(x, batch_size):