Bayesian Optimization

(require '[cljplot.build :as b]

         '[cljplot.core :refer :all]

         '[fastmath.core :as m]

         '[fastmath.kernel :as k]

         '[fastmath.random :as r]

         '[fastmath.optimization :as o]

         '[fastmath.classification :as cl]

         '[clojure2d.color :as c]

         '[clojure.java.io :as io]

         '[clojure.string :as str])

(defn black-box-function

  ^double [^double x ^double y]

  (inc (- (- (m/sq x))

       (m/sq (dec y)))))

(def bounds [[2 4] [-3 3]])

(black-box-function 1.0 0.0)

(o/maximize :bobyqa black-box-function {:bounds bounds})

(def bo (o/bayesian-optimization black-box-function {:bounds bounds}))

(map #(select-keys % [:x :y]) (take 5 bo))

(select-keys (nth bo 25) [:x :y])

(first bo)

(def palette (reverse (:rdylbu-9 c/palette-presets)))

(def gradient (c/gradient palette))

(defn draw-2d

  [f bounds bayesian-optimizer iteration]

  (let [{:keys [util-fn gp xs]} (nth bayesian-optimizer iteration)

        cfg {:x (first bounds)

             :y (second bounds)

             :palette palette

             :contours 30

             :gradient gradient}]

  (xy-chart {:width 700 :height 700}

            (b/series [:contour-2d f (assoc cfg :position [0 1])]

                      [:scatter xs {:size 10 

                                    :color (c/darken (first palette))

                                    :margins {:x [0 0] :y [0 0]}

                                    :position [0 1]

                                    :label "2d function with guessed points."}]

                      [:function-2d util-fn (assoc cfg :position [0 0] 

                                                   :label "Utility function")]

                      [:contour-2d (fn [x y] (gp [x y]))

                       (assoc cfg :position [1 1]

                              :label "Gaussian processes - mean (interpolation)")]

                      [:contour-2d (fn [x y] (second (gp [x y] true))) 

                       (assoc cfg :position [1 0]

                              :label "Gaussian processes - std dev")])

                  (b/add-axes :bottom)

                  (b/add-axes :left))))

(defn black-box-function2

  ^double [^double x ^double y]

  (+ 100.0 (* (m/sinc (+ 2.0 x y))

           (m/sin (+ x x))

           (m/cos (+ y y))

           (m/sin (* 2.0 x y)))))

(def bounds2 [[-1.0 0.5] [-2 -0.5]])

(def draw-black-box2 (partial draw-2d black-box-function2 bounds2))

(def bo (o/bayesian-optimization black-box-function2 {:bounds bounds2}))

(save (draw-black-box2 bo 0) "/results/black-box2.jpg")

(save (draw-black-box2 bo 10) "/results/black-box2-10.jpg")

(save (draw-black-box2 bo 40) "/results/black-box2-40.jpg")

(def bo (o/bayesian-optimization black-box-function2 {:bounds bounds2

                                                      :normalize? false}))

(save (draw-black-box2 bo 10) "/results/black-box2-not-normalized.jpg")

(def bo (o/bayesian-optimization black-box-function2 {:bounds bounds2

                                                      :noise 0.01

                                                      :utility-function-type :ucb

                                                      :utility-param 0.1}))

(save (draw-black-box2 bo 15) "/results/black-box2-ucb-low.jpg")

(def bo (o/bayesian-optimization black-box-function2 {:bounds bounds2

                                                      :noise 0.01

                                                      :utility-function-type :ucb

                                                      :utility-param 8}))

(save (draw-black-box2 bo 15) "/results/black-box2-ucb-high.jpg")

(def bo (o/bayesian-optimization black-box-function2 {:bounds bounds2

                                                      :noise 0.01

                                                      :utility-function-type :poi

                                                      :utility-param 0.01}))

(save (draw-black-box2 bo 15) "/results/black-box2-poi-low.jpg")

(def bo (o/bayesian-optimization black-box-function2 {:bounds bounds2

                                                      :noise 0.01

                                                      :utility-function-type :poi

                                                      :utility-param 0.9}))

(save (draw-black-box2 bo 15) "/results/black-box2-poi-high.jpg")

(def bo (o/bayesian-optimization black-box-function2 {:bounds bounds2

                                                      :noise 0.01

                                                      :utility-function-type :ei

                                                      :utility-param 0.01}))

(save (draw-black-box2 bo 15) "/results/black-box2-ei-low.jpg")

(def bo (o/bayesian-optimization black-box-function2 {:bounds bounds2

                                                      :noise 0.01

                                                      :utility-function-type :ei

                                                      :utility-param 0.9}))

(save (draw-black-box2 bo 15) "/results/black-box2-ei-high.jpg")

(def bo (o/bayesian-optimization black-box-function2 {:bounds bounds2

                                                      :optimizer :cmaes

                                                      :utility-function-type :ei

                                                      :kernel (k/kernel :mattern-12 0.1)

                                                      :noise 0.1}))

(save (draw-black-box2 bo 15) "/results/black-box2-narrow.jpg")

(def bo (o/bayesian-optimization black-box-function2 {:bounds bounds2

                                                      :utility-function-type :ei

                                                      :kernel (k/kernel :mattern-12 2)

                                                      :noise 0.1}))

(save (draw-black-box2 bo 15) "/results/black-box2-wide.jpg")

(def dataset (with-open [data (io/reader 
sonar.csv)]

             (mapv #(-> %

                        (str/trim)

                        (str/split #",")

                        (->> (map read-string))) (line-seq data))))

(count dataset)

(def xs (map butlast dataset))

(def ys (map (comp keyword last) dataset))

(frequencies ys)

(defn ada-boost-params

  ^double [^double trees ^double nodes]

  (let [ada-boost (cl/ada-boost {:number-of-trees (int trees)

                                 :max-nodes (int nodes)} xs ys)]

    (:accuracy (cl/cv ada-boost))))

(def ada-bounds [[1 (count dataset)]

                 [2 (count dataset)]])

(ada-boost-params 10 20)

(def ada-boost-bo (o/bayesian-optimization ada-boost-params

                                           {:init-points 5

                                            :kernel (k/kernel :mattern-52 55.0)

                                            :noise 0.05

                                            :utility-param 0.1

                                            :bounds ada-bounds}))

(select-keys (first ada-boost-bo) [:x :y])

(select-keys (nth ada-boost-bo 20) [:x :y])

(defn draw-2d-2

  [bounds bayesian-optimizer iteration]

  (let [{:keys [gp xs]} (nth bayesian-optimizer iteration)

        cfg {:x (first bounds)

             :y (second bounds)

             :palette palette

             :contours 30}]

    (xy-chart {:width 700 :height 300}

              (b/series [:contour-2d (fn [x y] (gp [x y])) cfg]

                        [:contour-2d (fn [x y] (second (gp [x y] true))) 

                         (assoc cfg :position [1 0]

                                :label "Gaussian processes - std dev")]

                        [:scatter xs {:size 10 

                                      :color (c/darken (first palette))

                                      :margins {:x [0 0] :y [0 0]}

                                      :position [0 0]

                                      :label "Gaussian processes - mean (interpolation)"}])

              (b/add-axes :bottom)

              (b/add-axes :left))))

(save (draw-2d-2 ada-bounds ada-boost-bo 20) "/results/ada-boost.jpg")

(defn svm-params

  ^double [^double cp ^double cn]

  (let [cp (m/pow 10 cp)

        cn (m/pow 10 cn)

        svm (cl/svm {:c-or-cp cp :cn cn

                     :kernel (k/kernel :gaussian)} xs ys)]

    (m/log (:accuracy (cl/cv svm)))))

(def svm-bounds [[-6 6] [-6 6]])

(m/exp (svm-params 1 1))

(def svm-bo (o/bayesian-optimization svm-params 

                                     {:init-points 5

                                      :kernel (k/kernel :mattern-52 1.5)

                                      :utility-param 0.3

                                      :noise 0.05

                                      :bounds svm-bounds}))

(m/exp (:y (nth svm-bo 10)))

(save (draw-2d-2 svm-bounds svm-bo 10) "/results/svm.jpg")