         weight-gradients [(mat/outer-product output-delta layer-activations)]
 
         ;; Now compute deltas and gradients for the hidden layers
-        layer-indices (reverse (range 1 (dec (count layer-sizes))))
+        layer-indices (reverse (range 1 (dec n-layers)))
         [_ bias-gradients weight-gradients]
         (reduce
           (fn [[delta bias-grads weight-grads] i]

;       (v = velocity, m = momentum, n=learning rate)

; TODO: switch to passing layer-specs that have size and activation fn
(defn network
  [layer-sizes & opts]
  {:n-layers (count layer-sizes)
   :layer-sizes layer-sizes
   :biases (map rand-vector (next layer-sizes))
   :weights (map weight-matrix (rest layer-sizes) (drop-last layer-sizes))
   :cost-fn (get opts :cost-fn (QuadraticCost.))})

(defn feed-forward
  [{:keys [biases weights] :as net} input]

;; 4) update weights
;;  * multiply gradient by the learning-rate to smooth out jitter in updates.
(defn backprop
  [{:keys [biases weights cost-fn n-layers] :as net} input expected-output]
  (let [bias-gradients (map #(mat/zero-array (mat/shape %)) biases)
        weight-gradients (map #(mat/zero-array (mat/shape %)) weights)
        [activations zs] (reduce
