Add differential evolution.

Update dependencies.
Tweak stats mean-col etc.

project.clj

@@ -7,7 +7,7 @@
 ;; the terms of this license.
 ;; You must not remove this notice, or any other, from this software.
 
-(defproject matlib "0.1.8"
+(defproject matlib "0.1.9"
   :description "A Clojure library of optimisation and control theory tools and convenience functions based on Neanderthal."
   :license {:name "EPL-2.0 OR GPL-2.0-or-later WITH Classpath-exception-2.0"
             :url "https://www.eclipse.org/legal/epl-2.0/"}
@@ -17,10 +17,11 @@
   ; jvm-opts required by Neanderthal on JDK > 8.
   :jvm-opts ^:replace ["--add-opens=java.base/jdk.internal.ref=ALL-UNNAMED"]
   :dependencies [[org.clojure/clojure "1.10.1"]
-                 [uncomplicate/neanderthal "0.34.0"]
+                 [uncomplicate/neanderthal "0.36.0"]
+                 ;[org.bytedeco/mkl-platform-redist "2020.1-1.5.3"]
                  ; No need to specify slf4j-api, as it’s required by logback.
                  ; These loggers are here mainly to keep Neanderthal quiet.
-                 [org.clojure/tools.logging "0.6.0"]
+                 [org.clojure/tools.logging "1.1.0"]
                  [ch.qos.logback/logback-classic "1.2.3"]]
   ; lein with-env-vars repl
   ; Is this working? Specified in $PATH in .bashrc.

src/matlib/de.clj

@@ -0,0 +1,72 @@
+(ns matlib.de
+  "Differential evolution heuristic gradient-free optimisation.
+  See, for example, https://en.wikipedia.org/wiki/Differential_evolution ."
+  (:require
+      [matlib.core :as ml-core]
+      [matlib.linalg :as ml-linalg]
+      [uncomplicate.neanderthal.real :refer [entry entry!]]
+      [uncomplicate.neanderthal.native :as native]
+      [uncomplicate.neanderthal.linalg :as n-linalg]
+      [uncomplicate.neanderthal.core :as n-core]
+      [uncomplicate.neanderthal.vect-math :as vect-math]
+      [uncomplicate.neanderthal.random :as random]))  
+
+(defn perturb
+  "Perturb a state vector `x`."
+  [x]
+  (let [x0 (native/dv x)]
+    (while (= 0.0 (reduce * x0)) (n-core/axpy! (random/rand-normal! (native/dv x)) x0))
+    (vect-math/mul x0 (random/rand-normal! (native/dv x)))))
+
+(defn population
+  "Initialise a population based on an example `x`.
+  `x` cannot have zero entries.
+  `NP` must be at least 5."
+  ([x]
+   (let [d (n-core/dim x)
+         NP (cond (< d 10) (* 5 d)
+                  (> d 100) (int (/ d 2))
+                  :else 50)]
+     (population x NP)))
+  ([x NP]
+   (for [n (range NP)] (perturb x))))
+
+(defn update-individual
+  "Returns updated `x`."
+  [x a b c CR F]
+  (let [R (rand-int (n-core/dim x))
+        y (n-core/copy x)]
+    (doseq
+      [i (range (n-core/dim x))]
+      (if (or (< (rand) CR) (= i R))
+        (entry! y i (+ (entry a i) (* F (- (entry b i) (entry c i)))))))
+    y))
+
+(defn solve
+  "Find the minimum of `f` given a population (vec) of `x`s.
+  Constraints should be handled in `f`.
+  If evaluations of `f` are expensive, consider memoizing `f` with lru."
+  ([f xs CR F maxiter scores n]
+   (let [x (first xs)
+         ys (shuffle (rest xs))
+         a (first ys)
+         b (nth ys 2)
+         c (nth ys 3)]
+     ;(print n scores "\r")
+     (if (and (> (- (reduce max scores) (reduce min scores)) ml-core/sq-eps) (< n maxiter))
+       ; rotate xs with new one at end
+       (let [y (update-individual x a b c CR F)
+             fx (f x)
+             fy (f y)
+             new-xs (conj (vec (rest xs)) (if (> fy fx) x y))
+             new-scores (conj (vec (rest scores)) (min fx fy))]
+         (recur f new-xs CR F maxiter new-scores (inc n)))
+       {:sol x
+        :f (f x)
+        :iterations n
+        :maxiter maxiter
+        :CR CR
+        :F F
+        :NP (count xs)
+        :success (< n maxiter)}))))
+

src/matlib/optim.clj

@@ -1,6 +1,8 @@
 (ns matlib.optim
   "Various optimisation algorithms.
   L-BFGS and gradient descent are implemented.
+  The gradient is estimated by central finite difference.
+  Differential evolution (gradient-free heuristic optimisation) is implemented.
 
   [NW-06]  
   'Numerical Optimization (second Ed.)'  
@@ -19,12 +21,14 @@
   (:require
     [matlib.core :refer :all]
     [matlib.linalg :refer :all]
+    [matlib.de]
     [uncomplicate.neanderthal.real :refer [entry entry!]]
     [uncomplicate.neanderthal.native :refer :all :exclude [sv]]
     [uncomplicate.neanderthal.linalg :refer :all]
     [uncomplicate.neanderthal.core :refer :all :exclude [entry entry!]]
     [uncomplicate.neanderthal.vect-math :as vect-math]
-    [uncomplicate.neanderthal.random :as random]))
+    [uncomplicate.neanderthal.random :as random]
+    [clojure.core.memoize :as memo]))
 
 (def ^:private ip (/ (- (Math/sqrt 5.0) 1.0) 2.0))
 (def ^:private ip2 (/ (- 3.0 (Math/sqrt 5.0)) 2.0))
@@ -296,6 +300,29 @@
        (print "\t\ta: steplength\ts: step\tq: df/dx\tp: search dirn\n"))
      (merge (alg-7-5 f (copy x) S Y X 0 linesearch tol maxiter history output) {:lsmethod lsmethod}))))
 
+(defn de
+  "Differential evolution (a heuristic gradient-free optimisation).
+  Find the minimum of `f: ℝⁿ -> ℝ` given a population of `x`s.
+  `f`   function to be solved, `f: ℝⁿ -> ℝ`.  
+  `x`   initial solution guess (as a vector) from which the population is generated.  
+  options (default):  
+  `:CR`       combination rate (0.9)  
+  `:F`        differential weight (0.8)
+  `:maxiter`  maximum iterations (10000)  
+  note:  
+  `f` is memoized with a lru cache.  
+  Constraints should be handled in `f`.
+  "
+  ([f x & params]
+   (let [xs (if (:NP params) (matlib.de/population x (:NP params)) (matlib.de/population x))
+         {:keys [CR F maxiter scores n] :or {CR 0.9
+                                             F 0.8
+                                             maxiter 10000
+                                             scores (map f xs)
+                                             n 0}} params
+         mem-f (memo/lru f {} :lru/threshold (* 2 (count xs)))]
+     (matlib.de/solve f xs CR F maxiter scores n))))
+
 (defn- booth
   "Booth function f: ℝ² -> ℝ, minimum at f(1, 3) = 0."
   [v]

src/matlib/stats.clj

@@ -13,12 +13,12 @@
 (defn row-mean
   "Mean across rows."
   ([M]
-   (scal! (/ 1.0 (mrows M)) (dge (map sum (cols M))))))
+   (scal! (/ 1.0 (mrows M)) (dv (map sum (cols M))))))
 
 (defn col-mean
   "Mean across columns."
   ([M]
-   (scal! (/ 1.0 (ncols M)) (dge (map sum (rows M))))))
+   (scal! (/ 1.0 (ncols M)) (dv (map sum (rows M))))))
 
 (defn subtract-col-mean
   "Subtract the mean over columns from each row."