       activation)))
 
 ;; Gradient descent
-
+;; 1) forward propagate
+;;  * sum((weights * inputs) + biases) for each neuron, layer by layer
+;; 2) back propagate deltas
+;;  * After propagating forward, we need to figure out the error (often called the
+;;    delta) for each neuron.  For the outputs this is just the
+;;    (output - expected_output) * activation-fn-prime, but for the hidden units the output error has to
+;;    be propagated back across the weights to distributed the error according to
+;;    which hidden units were responsible for it.
+;; 3) compute gradients
+;;  * Then for each weight you multiply its output delta by its input activation to get
+;;    its gradient.  This will correspond to the direction and magnitude of the
+;;    error for that weight, so any update should happen in the opposite direction.
+;; 4) update weights
+;;  * multiply gradient by the learning-rate to smooth out jitter in updates.
 (defn backprop
   [{:keys [biases weights layer-sizes] :as net} input expected-output]
-  (println "backprop")
-  (println "input: " input " => " expected-output)
-  (println "\nforward -----------------------------")
+  ;(println "\n--------------- backprop ---------------------")
   (let [bias-gradients (map #(mat/zero-array (mat/shape %)) biases)
         weight-gradients (map #(mat/zero-array (mat/shape %)) weights)
         [activations zs] (reduce

                             (let [z (mat/add (mat/mmul layer-weights (last activations)) layer-biases)
                                   activation (mat-sigmoid z)]
                               [(conj activations activation) (conj zs z)]))
                           [[input] []]
                           (map vector biases weights))
        _ (println "\nbackward -----------------------------")
        output (last activations)
        output-error (mat/sub output expected-output)
        output-deltas (mat/emul output-error (mat-sigmoid-prime (last zs)))
        ; TODO: I don't like that we have to reshape here...
        ;output-deltas (mat/reshape output-deltas [1 (first (mat/shape output-deltas))])
        bias-gradients [output-deltas]
        layer-activations (mat/transpose [(last (drop-last activations))])
        ;layer-activations (mat/reshape layer-activations [1 (first (mat/shape layer-activations))])
        _ (println "output-deltas: " (mat/shape output-deltas))
        ;_ (println "activations: " (mat/shape layer-activations))
        output-weight-gradients (mat/mmul output-deltas layer-activations)
        weight-gradients [output-weight-gradients]
        layer-indices (reverse (range (- (count layer-sizes) 2)))
        [_ bias-gradients weight-gradients]
        (reduce
          (fn [[delta d-biases d-weights] i] ; layer weights, zs and activations, in top down order
            (let [lws (nth weights i)
                  lzs (nth zs i)
                  las (nth activations i)
                  sp (mat-sigmoid-prime lzs)
                  _ (println "ws: " (count weights) " i: " (inc i))
                  _ (println "weights: " (mat/shape (mat/transpose (nth weights (inc i)))))
                  _ (println "deltas: " (mat/shape delta))
                  _ (println "sp: " (mat/shape sp) sp)
                  weight-delta (mat/mmul (mat/transpose (nth weights (inc i))) delta)
                  _ (println "deltas: " (mat/shape weight-delta))
                  delta (mat/mmul weight-delta sp)
                  weight-grad (mat/mmul delta (mat/transpose (nth activations (dec i))))
                  d-biases (cons delta d-biases)
                  d-weights (cons weight-grad d-weights)]
              [delta d-biases d-weights]))
          [output-deltas bias-gradients weight-gradients]
          layer-indices)]
    [bias-gradients weight-gradients]))


(defn update-mini-batch
  [{:keys [biases weights] :as net} learning-rate data labels batch]
  (println "update-mini-batch")
  (let [bias-gradients (map #(mat/zero-array (mat/shape %)) biases)
        weight-gradients (map #(mat/zero-array (mat/shape %)) weights)
        batch-size (count batch)

(defn xor-test
  []
  (let [net (network [2 3 1])
        optim-options {:n-epochs 10000
                       :batch-size 2
                       :learning-rate 0.1}
        _ (println "initial: " (:weights net))
        trained (sgd net optim-options XOR-DATA XOR-LABELS)
        _ (println "trained: " (:weights net))
        score (evaluate trained XOR-DATA XOR-LABELS)
        label-count (count XOR-LABELS)
        score-percent (float (/ score label-count))]
