       (m/join (m/as-vector (:weight-gradient this)) (m/as-vector (:bias-gradient this)))))
 
 
-(defn max-pooling-forward-default!
-  "The forward pass writes to two things; output and output-indexes.  The indexes
-record which item from the input row we actually looked at.  Returns
-[output output-indexes].  Note that this does not change the number of channels
-so it needs to remember the max per-channel."
-  [input output output-indexes
-   {:keys [^long height ^long width
-           ^long padx ^long pady
-           ^long k-width ^long k-height
-           ^long stride-w ^long stride-h
-           ^long num-channels] :as conv-layer-config}]
-  ;;Each input row has k-width*num-channels*k-height in it.
-  ;;Each output index gets num-channels written to it.
-  (let [conv-sequence (convolution-operation-sequence conv-layer-config)
-        height (long height)
-        width (long width)
-        padx (long padx)
-        pady (long pady)
-        k-width (long k-width)
-        k-height (long k-height)
-        stride-w (long stride-w)
-        stride-h (long stride-h)
-        num-channels (long num-channels)
-        output-width (get-padded-strided-dimension width padx k-width stride-w)]
-    (doseq [conv-item conv-sequence]
-      (let [{:keys [^long input-x ^long input-y ^long conv-x ^long conv-y
-                    ^long output-x ^long output-y]} conv-item
-            input-x (long input-x)
-            input-y (long input-y)
-            conv-x (long conv-x)
-            conv-y (long conv-y)
-            output-x (long output-x)
-            output-y (long output-y)
-            output-offset (* num-channels (+ (* output-y output-width) output-x))]
-        (loop [conv-x conv-x]
-          (when (< conv-x k-width)
-            (let [input-offset-x (+ input-x conv-x)
-                  valid-input? (and (>= input-y 0)
-                                    (< input-y height)
-                                    (>= input-offset-x 0)
-                                    (< input-offset-x width))
-                  kernel-index (+ (* conv-y k-width) conv-x)
-                  input-offset (* num-channels (+ (* input-y width) input-offset-x))]
-              (loop [chan 0]
-                (when (< chan num-channels)
-                  (let [input-val (double (if valid-input?
-                                            (m/mget input (+ input-offset chan))
-                                            0.0))
-                        output-offset (+ output-offset chan)
-                        existing-value (double (m/mget output output-offset))]
-                    (when (or (= kernel-index 0)
-                              (> input-val existing-value))
-                      (m/mset! output output-offset input-val)
-                      (m/mset! output-indexes output-offset kernel-index)))
-                  (recur (inc chan)))))
-            (recur (inc conv-x))))))
-    [output output-indexes]))
-
-
-#?(:cljs(defn max-pooling-forward!
-          [input output output-indexes conv-config]
-          (max-pooling-forward-default! input output output-indexes conv-config))
-
-   :clj(defn max-pooling-forward!
-         [input output output-indexes conv-config]
-         (let [input-ary (mp/as-double-array input)
-               output-ary (mp/as-double-array output)
-               output-indexes-ary (mp/as-double-array output-indexes)
-               {:keys [width height k-width k-height padx pady stride-w stride-h
-                       num-channels]} conv-config]
-           (if (and input-ary output-ary output-indexes-ary)
-             (ConvOps/maxPooling width height k-width k-height padx pady
-                                 stride-w stride-h num-channels
-                                 input-ary output-ary output-indexes-ary)
-             (max-pooling-forward-default! input output output-indexes conv-config)))))
-
-
-(defn max-pooling-backward!
+(defn planar-max-pooling-forward!
+  [layer input]
+  (let [^ConvLayerConfig config (:conv-config layer)
+        output-width (output-width config)
+        output-height (output-height config)
+        n-output (* output-width output-height (.num-channels config))
+        output (or (:output layer)
+                   (b/new-array [n-output]))
+        ^ints output-indexes (or (:output-indexes layer)
+                                 (int-array n-output))
+        ^doubles input-ary (mp/as-double-array input)
+        ^doubles output-ary (mp/as-double-array output)]
+    (convolution-outer-kernel
+     config
+     (convolution-roll-unroll-inner-kernel
+      (let [input-val (double (if input-valid?
+                                (aget input-ary input-addr)
+                                0.0))
+            output-addr (+ (* out-y output-width)
+                           out-x
+                           chan-output-offset)
+            k-idx (+ (* k-y k-width) k-x)
+            output-val (aget output-ary output-addr)]
+        (when (or (= 0 k-idx)
+                  (> input-val output-val))
+          (aset output-indexes output-addr k-idx)
+          (aset output-ary output-addr output-val)))))
+    (assoc layer
+           :output output
+           :output-indexes output-indexes)))
+
+
+(defn planar-max-pooling-backward!
   "Calculates the input gradient using the inverse of the convolution step
-combined with the output gradient and the output indexes which tell you which kernel
-index the output came from."
-  [output-gradient output-indexes input-gradient
-   {:keys [^long width ^long height
-           ^long k-width ^long k-height
-           ^long padx ^long pady
-           ^long stride-w ^long stride-h
-           ^long num-channels] :as conv-layer-config}]
-  (m/mset! input-gradient 0.0)
-  (let [output-width (get-padded-strided-dimension width padx k-width stride-w)
-        output-height (get-padded-strided-dimension height pady k-height stride-h)
-        num-pixels (* output-width output-height)]
-    (loop [pixel 0]
-      (when (< pixel num-pixels)
-        (let [output-offset (* pixel num-channels)
-              output-x (rem pixel output-width)
-              output-y (quot pixel output-width)]
-          (loop [chan 0]
-            (when (< chan num-channels)
-              (let [output-offset (+ output-offset chan)
-                    kernel-index (long (m/mget output-indexes output-offset))
-                    conv-x (rem kernel-index k-width)
-                    conv-y (quot kernel-index k-width)
-                    input-x (- (+ (* output-x stride-w) conv-x) padx)
-                    input-y (- (+ (* output-y stride-h) conv-y) pady)]
-                (when (and (>= input-y 0)
-                           (< input-y height)
-                           (>= input-x 0)
-                           (< input-x width))
-                  (let [^double output-gradient-value (m/mget output-gradient output-offset)
-                        input-offset (+ (* num-channels (+ (* input-y width)
-                                                           input-x))
-                                        chan)
-                        ^double input-gradient-value (m/mget input-gradient input-offset)]
-                    (m/mset! input-gradient input-offset (+ input-gradient-value
-                                                            output-gradient-value)))))
-              (recur (inc chan)))))
-        (recur (inc pixel))))))
-
+  combined with the output gradient and the output indexes which tell you which kernel
+  index the output came from."
+  [layer input output-gradient]
+  (let [^ConvLayerConfig config (:conv-config layer)
+        ^ints output-indexes (:output-indexes layer)
+        input-gradient (or (:input-gradient layer)
+                           (b/new-array (m/shape input)))
+        ^doubles input-ary (mp/as-double-array input-gradient)
+        ^doubles output-ary (mp/as-double-array output-gradient)]
+    ;;Zero accumulator
+    (m/fill! input-gradient 0.0)
+    (convolution-outer-kernel
+     config
+     (let [output-addr (+ (* out-y output-width)
+                          out-x
+                          chan-output-offset)
+           k-idx (aget output-indexes output-addr)
+           output-val (aget output-ary output-addr)
+           k-y (quot k-idx k-width)
+           k-x (rem k-idx k-width)
+           input-x (+ input-rel-x k-x)
+           input-y (+ input-rel-y k-y)
+           addr (+ (* input-y width)
+                   input-x
+                   chan-input-offset)]
+       (when (and (> input-x 0)
+                  (<= input-x width)
+                  (> input-y 0)
+                  (<= input-y height))
+         (aset input-ary addr output-val))))
+    (assoc layer :input-gradient input-gradient)))
 
 ;;Max pooling layer.  There are other pooling layer types (average,a sochiastic)
 ;;that may be implemented later but for now we only need max pooling.

(defrecord Pooling [output output-indexes input-gradient conv-layer-config]
  cp/PModule
  (cp/calc [this input]
    (max-pooling-forward! input output output-indexes conv-layer-config))

  (cp/output [m]
    (:output m))

  cp/PNeuralTraining
  (forward [this input]
    (cp/calc this input)
    this)

  (backward [this input output-gradient]
    (max-pooling-backward! output-gradient output-indexes
                           input-gradient conv-layer-config)
    this)

  (input-gradient [this]
    (:input-gradient this)))

                     padx pady stride-w stride-h num-kernels
                     num-channels))

(defmacro convolution-outer-kernel
  [config & body]
  `(let [^ConvLayerConfig config# ~config
         ~'output-width (output-width config#)
         ~'output-height (output-height config#)
         ~'num-channels (.num-channels config#)
         ~'num-kernels (.num-kernels config#)
         ~'input-planar-stride (* (.width config#) (.height config#))
         ~'output-planar-stride (* output-width output-height ~'num-kernels)
         ~'output-channel-stride (* (.k-width config#) (.k-height config#))
         ~'output-column-stride (* ~'output-channel-stride ~'num-channels)
         ~'width (.width config#)
         ~'height (.height config#)
         ~'k-width (.k-width config#)
         ~'k-height (.k-height config#)
         ~'stride-h (.stride-h config#)
         ~'stride-w (.stride-w config#)
         ~'padx (.padx config#)
         ~'pady (.pady config#)]
     (c-for
      [~'chan 0 (< ~'chan ~'num-channels) (inc ~'chan)]
      (let [~'chan-input-offset (* ~'chan ~'input-planar-stride)
            ~'chan-output-offset (* ~'chan ~'output-planar-stride)]
       (c-for
        [~'out-y 0 (< ~'out-y ~'output-height) (inc ~'out-y)]
        (let [~'input-rel-y (- (* ~'out-y ~'stride-h) ~'pady)]
          (c-for
           [~'out-x 0 (< ~'out-x ~'output-width) (inc ~'out-x)]
           (let [~'input-rel-x (- (* ~'out-x ~'stride-w) ~'padx)]
             ~@body))))))))


(defmacro convolution-roll-unroll-inner-kernel
  [& body]
  `(let [~'chan-conv-offset (* ~'chan ~'output-channel-stride)]
    (c-for
     [~'k-y 0 (< ~'k-y ~'k-height) (inc ~'k-y)]
     (c-for
      [~'k-x 0 (< ~'k-x ~'k-width) (inc ~'k-x)]
      (let [~'input-x (+ ~'input-rel-x ~'k-x)
            ~'input-y (+ ~'input-rel-y ~'k-y)
            ~'output-conv-addr (+ (* ~'out-y ~'out-x ~'output-column-stride)
                                  ~'chan-conv-offset
                                  (* ~'k-y ~'k-width)
                                  ~'k-x)
            ~'input-addr  (+ (* ~'input-y ~'width)
                             ~'input-x
                             ~'chan-input-offset)
            ~'input-valid? (and (> ~'input-x 0)
                                (<= ~'input-x ~'width)
                                (> ~'input-y 0)
                                (<= ~'input-y ~'height))]
        ~@body)))))



(defn planar-input->convolution!
  [input output ^ConvLayerConfig config]
  (let [^doubles input-ary (mp/as-double-array input)
        ^doubles output-ary (mp/as-double-array output)]
    (convolution-outer-kernel
     config
     (convolution-roll-unroll-inner-kernel
      (let [input-val (double (if input-valid?
                                (aget input-ary input-addr)
                                0.0))]
        (aset output-ary output-conv-addr input-val)))))
  output)


(defn convolution->planar-output!
  "Sum the convolution up to the planar input."
  [conv-input-gradient input-gradient ^ConvLayerConfig config]
  ;;I am using input to mean upstream or in this case destination so that
  ;;this code can look as similar to the code above as possible
  (let [^doubles input-ary (mp/as-double-array input-gradient)
        ^doubles output-ary (mp/as-double-array conv-input-gradient)]
    ;;Zero accumulator
    (m/fill! input-gradient 0.0)
    (convolution-outer-kernel
     config
     (convolution-roll-unroll-inner-kernel
      (when input-valid?
        (let [input-val (aget input-ary input-addr)
              output-val (aget output-ary output-conv-addr)]
          (aset input-ary input-addr (+ input-val output-val)))))))
  input-gradient)


(defn create-convolution-matrix
  [^ConvLayerConfig config]
  (let [output-width (output-width config)
        output-height (output-height config)
        kernel-stride (* (.k-width config) (.k-height config))
        n-cols (* kernel-stride (.num-channels config))
        n-rows (* output-width output-height)]
    (b/new-array [n-rows n-cols])))


(defn planar-convolution-forward!
  [layer input]
  (let [weights (:weights layer)
        bias (:bias layer)
        ^ConvLayerConfig config (:conv-config layer)
        conv-matrix (or (:conv-matrix layer)
                        (create-convolution-matrix config))
        output-width (output-width config)
        output-height (output-height config)
        output-stride (* output-width (.num-kernels config))
        output (or (:output layer)
                   (b/new-array [output-height output-stride]))]
    (planar-input->convolution! input conv-matrix config)
    (m/assign! output bias)
    (if (blas/supports-blas? weights)
      (blas/gemm! output false true 1.0 weights conv-matrix 1.0)
      (m/add! output (m/inner-product weights (m/transpose conv-matrix))))
    (assoc layer :conv-matrix conv-matrix :output output)))


(defn planar-convolution-backward!
  [layer input output-gradient]
  (let [weights (:weights layer)
        bias (:bias layer)
        ^ConvLayerConfig config (:conv-config layer)
        output-width (output-width config)
        output-height (output-height config)
        n-output (* output-width output-height)
        n-kernels (.num-kernels config)
        conv-matrix (:conv-matrix layer)
        weight-gradient (or (:weight-gradient layer)
                            (b/new-array (m/shape weights)))
        bias-gradient (or (:bias-gradient layer)
                          (b/new-array (m/shape bias)))
        output-gradient-matrix (or (:output-gradient-matrix layer)
                                   (b/new-array [n-kernels n-output]))
        input-gradient (or (:input-gradient layer)
                           (b/new-array (m/shape input)))]
    (m/assign! (m/as-vector output-gradient-matrix) output-gradient)
    (m/add! bias-gradient output-gradient)
    ;;conv-matrix is assumed to hold the actual input
    (if (blas/supports-blas? weights)
      (do
        (blas/gemm! weight-gradient false false 1.0 output-gradient-matrix conv-matrix 1.0)
        (blas/gemm! conv-matrix true false 1.0 output-gradient-matrix weights 0.0))
      (do
        (m/add! weight-gradient (m/inner-product output-gradient-matrix conv-matrix))
        (m/assign! conv-matrix (m/inner-product (m/transpose output-gradient-matrix) weights))))
    (convolution->planar-output! conv-matrix input-gradient config)
    (assoc layer
           :weight-gradient weight-gradient
           :bias-gradient bias-gradient
           :output-gradient-matrix output-gradient-matrix
           :input-gradient input-gradient)))


;; Forward: Take the input which is expected to be a single vector of data
;; and create a matrix that contains a row for each convolution.  So for example if you have
;; 2x2 kernels and you have a 3x3 matrix of monotonically incrementing indexes we produce
;; a new matrix
