(Arne Babenhauserheide)
2017-03-21: clean list benchmark plus fit clean list benchmark plus fit
diff --git a/examples/benchmark.w b/examples/benchmark.w
--- a/examples/benchmark.w
+++ b/examples/benchmark.w
@@ -1,16 +1,19 @@
#!/usr/bin/env sh
# -*- wisp -*-
-exec guile -L $(dirname $(dirname $(realpath "$0"))) --language=wisp -e '(@@ (examples benchmark) main)' -s "$0" "$@"
+exec guile -L $(dirname $(dirname $(realpath "$0"))) --language=wisp -e '(@@ (examples benchmark) main)' -l $(dirname $(realpath "$0"))/cholesky.w -l $(dirname $(realpath "$0"))/ensemble-estimation.w -s "$0" "$@"
; !#
define-module : examples benchmark
import : statprof
ice-9 optargs
+ ice-9 format
srfi srfi-1
+ srfi srfi-42 ; list-ec
ice-9 pretty-print
system vm program
+
define : benchmark-run fun
let profiler : : loop-num 100
statprof-start
@@ -76,20 +79,20 @@ define : benchmark-list-append
let : (N (list-ref x 0)) (m (list-ref x 1))
benchmark (append a b) :let ((a (iota N))(b (iota m)))
. param-list
- let : (steps 4)
+ let : (steps 30)
concatenate
list
let : (param-list (zip (iota steps 1 1000) (iota steps 1 0)))
bench-append param-list
- let : (param-list (zip (iota steps 1 0) (iota steps 1 100)))
+ let : (param-list (zip (iota steps 1 0) (iota steps 1 1000)))
bench-append param-list
let : (param-list (zip (iota steps 1 1000) (iota steps 1 0)))
bench-append param-list
- let : (param-list (zip (iota steps 1 0) (iota steps 1 100)))
+ let : (param-list (zip (iota steps 1 0) (iota steps 1 1000)))
bench-append param-list
let : (param-list (zip (iota steps 1 1000) (iota steps 100000 0)))
bench-append param-list
- let : (param-list (zip (iota steps 100000 0) (iota steps 1 100)))
+ let : (param-list (zip (iota steps 100000 0) (iota steps 1 1000)))
bench-append param-list
;; stddev from rosetta code: http://rosettacode.org/wiki/Standard_deviation#Scheme
@@ -106,12 +109,143 @@ define : running-stddev nums
stddev nums
. running-stddev-2
+
+;; prepare a multi-function fit
+import
+ only : examples ensemble-estimation
+ . EnSRF make-covariance-matrix-with-offdiagonals-using-stds
+ . standard-deviation-from-deviations x-deviations->y-deviations
+ . x^steps
+ only : ice-9 popen
+ . open-output-pipe close-pipe
+
+define : H x pos
+ . "Observation operator. It generates modelled observations from the input.
+
+x are parameters to be optimized, pos is another input which is not optimized. For plain functions it could be the position of the measurement on the x-axis. We currently assume absolute knowledge about the position.
+"
+ let : (N (first pos)) (m (second pos))
+ +
+ list-ref x 0 ; constant value
+ ;; pure N
+ * (list-ref x 1) : log : + 1 N ; avoid breakage at pos 0
+ ; * (list-ref x 2) : sqrt N
+ * (list-ref x 3) N
+ ; * (list-ref x 4) : expt N 2
+ ; * (list-ref x 5) : expt N 3
+ ;; pure m
+ * (list-ref x 6) : log : + 1 m ; avoid breakage at pos 0
+ ; * (list-ref x 7) : sqrt m
+ * (list-ref x 8) m
+ ; * (list-ref x 9) : expt m 2
+ ; * (list-ref x 10) : expt m 3
+ ;; mixed terms
+ * (list-ref x 11) : log : + 1 N m
+ * (list-ref x 12) : * N (log (+ 1 m))
+ * (list-ref x 13) : * m (log (+ 1 N))
+ ; * (list-ref x 14) : * N m
+ ; * (list-ref x 15) : * (expt N 2) m
+ ; * (list-ref x 16) : * (expt m 2) N
+
+
+define : interleave lx lz
+ cond
+ (null? lx) lz
+ else
+ cons : car lx
+ interleave lz : cdr lx
+
+
+define : print-fit x σ
+ let : : msg "
+~,1,,,,,'ee±~,1,,,,,'ee + ~,1,,,,,'ee±~,1,,,,,'ee log(N) + ~,1,,,,,'ee±~,1,,,,,'ee sqrt(N) + ~,1,,,,,'ee±~,1,,,,,'ee N + ~,1,,,,,'ee±~,1,,,,,'ee N^2 + ~,1,,,,,'ee±~,1,,,,,'ee N^3 +
+~,1,,,,,'ee±~,1,,,,,'ee log(m) + ~,1,,,,,'ee±~,1,,,,,'ee sqrt(m) + ~,1,,,,,'ee±~,1,,,,,'ee m + ~,1,,,,,'ee±~,1,,,,,'ee m^2 + ~,1,,,,,'ee±~,1,,,,,'ee m^3 +
+~,1,,,,,'ee±~,1,,,,,'ee log(N + m) + ~,1,,,,,'ee±~,1,,,,,'ee N log(m) + ~,1,,,,,'ee±~,1,,,,,'ee m log(N)+ ~,1,,,,,'ee±~,1,,,,,'ee N m + ~,1,,,,,'ee±~,1,,,,,'ee N^2 m + ~,1,,,,,'ee±~,1,,,,,'ee m^2 N
+"
+ apply format
+ append (list #t msg) (interleave x σ)
+
+
+define : flatten li
+ append-ec (: i li) i
+
define : main args
- map
- lambda : mismatch-fun
- write (procedure-name mismatch-fun)
- newline
- let : (mis (mismatch-fun (benchmark-list-append)))
- map : lambda (x) : pretty-print (stddev x)
- apply zip mis
- list mismatch-to-const-N-m mismatch-to-linear-N-m
+ ;; map
+ ;; lambda : mismatch-fun
+ ;; write (procedure-name mismatch-fun)
+ ;; newline
+ ;; let : (mis (mismatch-fun (benchmark-list-append)))
+ ;; map : lambda (x) : pretty-print (stddev x)
+ ;; apply zip mis
+ ;; list mismatch-to-const-N-m mismatch-to-linear-N-m
+ let*
+ : bench : benchmark-list-append ;; benchmark results
+ ; fitting to cost estimates
+ ensemble-member-count 32
+ ensemble-member-plot-skip 4
+ x^b : list-ec (: i 17) (car (cdr (car bench))) ; inital guess: constant starting at the first result
+ x^b-std : list-ec (: i 17) (car (cdr (car bench))) ; inital guess: 100% uncertainty
+ P : make-covariance-matrix-with-offdiagonals-using-stds x^b-std
+ y⁰-pos : map car bench
+ y⁰ : append-map cdr bench
+ y⁰-std : stddev y⁰
+ R : make-covariance-matrix-with-offdiagonals-using-stds : list-ec (: i (length bench)) y⁰-std
+ optimized : EnSRF H x^b P y⁰ R y⁰-pos ensemble-member-count
+ x-opt : list-ref optimized 0
+ x-deviations : list-ref optimized 1
+ x-std
+ list-ec (: i (length x-opt))
+ apply standard-deviation-from-deviations : list-ec (: j x-deviations) : list-ref j i
+ y-deviations : x-deviations->y-deviations H x-opt x-deviations y⁰-pos
+ y-stds : list-ec (: i y-deviations) : apply standard-deviation-from-deviations i
+ y-opt : map (λ (x) (H x-opt x)) y⁰-pos
+ x^b-deviations-approx
+ list-ec (: i ensemble-member-count)
+ list-ec (: j (length x^b))
+ * : random:normal
+ sqrt : list-ref (list-ref P j) j ; only for diagonal P!
+ y^b-deviations : x-deviations->y-deviations H x^b x^b-deviations-approx y⁰-pos
+ y-std
+ apply standard-deviation-from-deviations
+ flatten y-deviations
+ y-stds : list-ec (: i y-deviations) : apply standard-deviation-from-deviations i
+ y^b-stds : list-ec (: i y^b-deviations) : apply standard-deviation-from-deviations i
+
+ ;; print-fit x-std
+ print-fit x-opt x-std
+ ; now plot the result
+ let : : port : open-output-pipe "python2"
+ format port "import pylab as pl\nimport matplotlib as mpl\n"
+ format port "y0 = [float(i) for i in '~A'[1:-1].split(' ')]\n" y⁰
+ format port "yerr = ~A\n" y⁰-std
+ format port "ypos1 = [float(i) for i in '~A'[1:-1].split(' ')]\n" : list-ec (: i y⁰-pos) : first i
+ format port "ypos2 = [float(i) for i in '~A'[1:-1].split(' ')]\n" : list-ec (: i y⁰-pos) : second i
+ format port "yinit = [float(i) for i in '~A'[1:-1].split(' ')]\n" : list-ec (: i y⁰-pos) : H x^b i
+ format port "yinitstds = [float(i) for i in '~A'[1:-1].split(' ')]\n" y^b-stds
+ format port "yopt = [float(i) for i in '~A'[1:-1].split(' ')]\n" : list-ec (: i y⁰-pos) : H x-opt i
+ format port "yoptstds = [float(i) for i in '~A'[1:-1].split(' ')]\n" y-stds
+ format port "pl.errorbar(*zip(*sorted(zip(ypos1, yinit))), yerr=zip(*sorted(zip(ypos1, yinitstds)))[1], label='prior vs N')\n"
+ format port "pl.errorbar(*zip(*sorted(zip(ypos1, yopt))), yerr=zip(*sorted(zip(ypos1, yoptstds)))[1], label='optimized vs N')\n"
+ format port "eb=pl.errorbar(*zip(*sorted(zip(ypos1, y0))), yerr=yerr, alpha=0.6, marker='x', mew=2, ms=10, linewidth=0, label='measurements vs N')\neb[-1][0].set_linewidth(1)\n"
+ format port "pl.errorbar(*zip(*sorted(zip(ypos2, yinit))), yerr=zip(*sorted(zip(ypos2, yinitstds)))[1], label='prior vs. m')\n"
+ format port "pl.errorbar(*zip(*sorted(zip(ypos2, yopt))), yerr=zip(*sorted(zip(ypos2, yoptstds)))[1], label='optimized vs. m')\n"
+ format port "eb=pl.errorbar(*zip(*sorted(zip(ypos2, y0))), yerr=yerr, alpha=0.6, marker='x', mew=2, ms=10, linewidth=0, label='measurements vs. m')\neb[-1][0].set_linewidth(1)\n"
+ list-ec (: step 0 (length x^steps) 16)
+ let : : members : list-ref x^steps (- (length x^steps) step 1)
+ list-ec (: member-idx 0 (length members) ensemble-member-plot-skip) ; reversed
+ let : : member : list-ref members member-idx
+ format port "paired = pl.get_cmap('Paired')
+cNorm = mpl.colors.Normalize(vmin=~A, vmax=~A)
+scalarMap = mpl.cm.ScalarMappable(norm=cNorm, cmap=paired)\n" 0 (length member)
+ list-ec (: param-idx 0 (length member) 16) ; step = 16
+ ; plot parameter 0
+ format port "pl.plot(~A, ~A, marker='.', color=scalarMap.to_rgba(~A), linewidth=0, label='', alpha=0.6, zorder=-1)\n" (/ step 1) (+ 80 (* (/ (apply + y-opt) (length y-opt)) (list-ref member param-idx))) param-idx
+ format port "pl.legend(loc='upper right')\n"
+ format port "pl.xlabel('position [arbitrary units]')\n"
+ format port "pl.ylabel('value [arbitrary units]')\n"
+ format port "pl.title('ensemble optimization results')\n"
+ format port "pl.xscale('log')\n"
+ format port "pl.yscale('log')\n"
+ format port "pl.show()\n"
+ format port "exit()\n"
+ close-pipe port
diff --git a/examples/ensemble-estimation.w b/examples/ensemble-estimation.w
--- a/examples/ensemble-estimation.w
+++ b/examples/ensemble-estimation.w
@@ -34,7 +34,10 @@ exec guile -L $(dirname $(dirname $(real
;; x'^a = x'^b - αK·H(x'^b)
define-module : examples ensemble-estimation
- . #:export : EnSRF H standard-deviation-from-deviations make-covariance-matrix-with-offdiagonals-using-stds
+ . #:export (EnSRF H standard-deviation-from-deviations
+ make-covariance-matrix-with-offdiagonals-using-stds
+ x-deviations->y-deviations x^steps)
+
use-modules : srfi srfi-42 ; list-ec
srfi srfi-9 ; records
oop goops ; generic functions
@@ -276,8 +279,8 @@ Limitations: y is a single value. R and
if : equal? 1-AK '()
set! 1-AK : list-ec (: i K) {1 - i} ; init
set! 1-AK : list-ec (: i (length K)) : * (list-ref 1-AK i) {1 - (abs (list-ref K i))}
- display 1-AK ; TODO: What does this actually signify?
- newline
+ ;; display 1-AK ; TODO: What does this actually signify?
+ ;; newline
step
cdr observations-to-process
cdr observation-variances