wisp
 
(Arne Babenhauserheide)
2017-03-23: Refactor and save results for multiple tests

Refactor and save results for multiple tests

diff --git a/examples/benchmark.w b/examples/benchmark.w
--- a/examples/benchmark.w
+++ b/examples/benchmark.w
@@ -44,9 +44,9 @@ define : running-stddev nums
   . running-stddev-2
 
 define* : benchmark-run-single fun #:key (min-seconds 0.1)
+  ;; trigger garbage collection before stats collection to avoid polluting the data
+  gc
   let profiler : (loop-num 4)
-    ;; trigger garbage collection before stats collection to avoid polluting the data
-    gc
     let : : t : get-internal-real-time
       with-output-to-string
         lambda ()
@@ -62,15 +62,20 @@ define* : benchmark-run-single fun #:key
             /  seconds loop-num ;; this wastes less than {(4 * ((4^(i-1)) - 1)) / 4^i} fractional data but gains big in simplicity
             profiler (* 4 loop-num) ;; for fast functions I need to go up rapidly, for slow ones I need to avoid overshooting
 
-define* : benchmark-run fun #:key (max-relative-uncertainty 0.2)
-    pretty-print fun
-    let lp : (min-seconds 1.e-5) (sampling-steps 4) ;; min seconds: 10μs ~ 30k cycles, start with at least 3 sampling steps to make the approximations in stddev-unbiased-normal good enough
+;; Define targets for the data aquisition
+define max-relative-uncertainty 0.3 ;; 3 sigma from 0
+define min-aggregated-runtime-seconds 1.e-5 ;; 10μs ~ 30k cycles
+define max-absolute-uncertainty-seconds 1.e-3 ;; 1ms, required to ensure that the model uses the higher values (else they would have huge uncertainties). If you find you need more, use a smaller test case.
+
+define* : benchmark-run fun
+    ;; pretty-print fun
+    let lp : (min-seconds min-aggregated-runtime-seconds) (sampling-steps 4) ;; start with at least 3 sampling steps to make the approximations in stddev-unbiased-normal good enough
         let*
           : res : list-ec (: i sampling-steps) : benchmark-run-single fun #:min-seconds min-seconds
             std : stddev-unbiased-normal res
             mean : / (apply + res) (length res)
-           pretty-print : list mean '± std min-seconds sampling-steps
-           if : < std {mean * max-relative-uncertainty}
+           ;; pretty-print : list mean '± std min-seconds sampling-steps
+           if : and {std < {mean * max-relative-uncertainty}} {std < max-absolute-uncertainty-seconds}
               . mean
               lp (* 2 min-seconds) (* 2 sampling-steps) ;; should decrease σ by factor 2 or √2 (for slow functions)
 
@@ -109,30 +114,30 @@ define : logiota steps start stepsize
           logstep : / (- (log (+ start (* stepsize (- steps 1)))) logstart) (- steps 1)
         map inexact->exact : map round : map exp : iota steps logstart logstep 
 
-define : benchmark-list-append
-  . "Test (append a b) with lists of different lengths."
-  define : bench-append param-list
-    zip param-list 
-        map 
-         lambda (x)
-           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 10)
-    concatenate
-      list 
-        let : (param-list (zip (logiota steps 1 10000) (logiota steps 1 0)))
-               bench-append param-list
-        ;; let : (param-list (zip (logiota steps 20 0) (logiota steps 1 10000)))
-        ;;        bench-append param-list
-        ;; let : (param-list (zip (logiota steps 1 1000) (logiota steps 1 0)))
-        ;;        bench-append param-list
-        ;; let : (param-list (zip (logiota steps 1 0) (logiota steps 1 1000)))
-        ;;        bench-append param-list
-        ;; let : (param-list (zip (logiota steps 1 1000) (logiota steps 100000 0)))
-        ;;        bench-append param-list
-        ;; let : (param-list (zip (logiota steps 100000 0) (logiota steps 1 1000)))
-        ;;        bench-append param-list
+define : bench-append param-list
+  . "Test (append a b) with lists of lengths from the param-list."
+  zip param-list 
+      map 
+       lambda (x)
+         let : (N (list-ref x 0)) (m (list-ref x 1))
+             benchmark (append a b) :let ((a (iota N))(b (iota m)))
+       . param-list
+
+define : benchmark-list-append-N-1 steps
+  . "Test (append a b) with lists with lengths N and 1."
+  bench-append (zip (logiota steps 1 1000) (logiota steps 1 0))
+
+define : benchmark-list-append-N-100 steps
+  . "Test (append a b) with lists with lengths N and 1."
+  bench-append (zip (logiota steps 1 1000) (logiota steps 100 0))
+
+define : benchmark-list-append-1-m steps
+  . "Test (append a b) with lists with lengths N and 1."
+  bench-append (zip (logiota steps 1 0) (logiota steps 1 1000))
+
+define : benchmark-list-append-100-m steps
+  . "Test (append a b) with lists with lengths N and 1."
+  bench-append (zip (logiota steps 100 0) (logiota steps 1 1000))
 
 
 ;; prepare a multi-function fit
@@ -156,13 +161,13 @@ define*
        H-N-m x pos #:key all const OlogN OsqrtN ON ONlogN ON²
                                  . Ologm Osqrtm Om Omlogm Om²
                                  . OlogNm ONlogm OmlogN ONm
-                                 . ON²m Om²N OsinN/N Osinm/m
+                                 . ON²m Om²N
        . "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.
 "
        when all
-           let lp : (l '(const OlogN OsqrtN ON ONlogN ON² Ologm Osqrtm Om Omlogm Om² OlogNm ONlogm OmlogN ONm ON²m Om²N OsinN/N Osinm/m))
+           let lp : (l '(const OlogN OsqrtN ON ONlogN ON² Ologm Osqrtm Om Omlogm Om² OlogNm ONlogm OmlogN ONm ON²m Om²N))
                when : not : null? l
                       define-quoted (car l) #t
                       lp : cdr l
@@ -189,8 +194,6 @@ x are parameters to be optimized, pos is
              or0 ONm    : * (list-ref x 14) : * N m
              or0 ON²m   : * (list-ref x 15) : * (expt N 2) m
              or0 Om²N   : * (list-ref x 16) : * (expt m 2) N
-             or0 OsinN/N : * (list-ref x 17) : / (sin (/ (- N (list-ref x 18)) (list-ref x 19))) N ; sin(x)/x
-             or0 Osinm/m : * (list-ref x 20) : / (sin (/ (- m (list-ref x 21)) (list-ref x 22))) m ; sin(x)/x
 
 
 define : interleave lx lz
@@ -205,7 +208,7 @@ define : print-fit x σ
     . "Print the big-O parameters which are larger than σ (their standard deviation)."
     let : : number-format "~,1,,,,,'ee±~,1,,,,,'ee"
       let big-O
-        : names : list "" "log(N)" "sqrt(N)" "N log(N)" "N^2" "log(m)" "sqrt(m)" "m" "m log(m)" "m^2" "log(N + m)" "N log(m)" "m log(N)" "N m" "N^2 m" "m^2 N" "sin((N-x)/p)/N" "sin((N-x*)/p)/N" "sin((N-x)/p*)/N" "sin((m-x)/p)/m" "sin((m-x*)/p)/m" "sin((m-x)/p*)/m" 
+        : names : list "" "log(N)" "sqrt(N)" "N log(N)" "N^2" "log(m)" "sqrt(m)" "m" "m log(m)" "m^2" "log(N + m)" "N log(m)" "m log(N)" "N m" "N^2 m" "m^2 N"
           x x
           σ σ
         cond
@@ -223,20 +226,20 @@ define : flatten li
          append-ec (: i li) i
 
 ;; TODO: add filename and title and fix the units
-define* : plot-benchmark-result bench H
+define* : plot-benchmark-result bench H #:key filename
      let*
-        : ensemble-member-count 256
-          ensemble-member-plot-skip 4
+        : ensemble-member-count 128
+          ensemble-member-plot-skip 0
           y_0 : apply min : map car : map cdr bench
           y_m : apply max : map car : map cdr bench
           nb : apply max : interleave (map car (map car bench)) (map car (map cdr (map car bench)))
           ;; "const" "log(N)" "sqrt(N)" "N" "N^2" "N^3" "log(m)" "sqrt(m)" "m" "m^2" "m^3" "log(N + m)" "N log(m)" "m log(N)" "N m" "N^2 m" "m^2 N"
-          x^b : list y_0 (/ y_m (log nb)) (/ y_m (sqrt nb)) (/ y_m nb) (/ y_m nb nb) (/ y_m nb nb nb) (/ y_m (log nb)) (/ y_m (sqrt nb)) (/ y_m nb) (/ y_m nb nb) (/ y_m nb nb nb) (/ y_m nb nb) (/ y_m nb nb) (/ y_m nb nb nb) (/ y_m nb nb nb) (/ y_m nb nb nb nb) (/ y_m nb nb nb nb) y_0 (/ nb 100) 10000 y_0 (/ nb 100) 10000 ; inital guess: constant starting at the first result
+          x^b : list y_0 (/ y_m (log nb)) (/ y_m (sqrt nb)) (/ y_m nb) (/ y_m nb nb) (/ y_m nb nb nb) (/ y_m (log nb)) (/ y_m (sqrt nb)) (/ y_m nb) (/ y_m nb nb) (/ y_m nb nb nb) (/ y_m nb nb) (/ y_m nb nb) (/ y_m nb nb nb) (/ y_m nb nb nb) (/ y_m nb nb nb nb) (/ y_m nb nb nb nb)  ; inital guess: constant starting at the first result
           x^b-std : list-ec (: i x^b) i ; 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⁰-stds : list-ec (: i y⁰) {i * 0.2} ; enforcing 20% max std in benchmark-run
+          y⁰-stds : list-ec (: i y⁰) : min max-absolute-uncertainty-seconds {max-relative-uncertainty * i} ; enforcing 20% max std in benchmark-run
           y⁰-std : list-ref (sort y⁰ <) : round : / (length y⁰) 8 ; lower octile median
           R : make-covariance-matrix-with-offdiagonals-using-stds y⁰-stds
           optimized : EnSRF H x^b P y⁰ R y⁰-pos ensemble-member-count
@@ -268,6 +271,7 @@ define* : plot-benchmark-result bench H
         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 "ystds = [float(i) for i in '~A'[1:-1].split(' ')]\n" y⁰-stds
           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
@@ -277,10 +281,10 @@ define* : plot-benchmark-result bench H
           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], marker='H', mew=0, ms=10, linewidth=0.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 "eb=pl.errorbar(*zip(*sorted(zip(ypos1, y0))), yerr=ystds, 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], marker='h', mew=0, ms=10, linewidth=0.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"
+          format port "eb=pl.errorbar(*zip(*sorted(zip(ypos2, y0))), yerr=ystds, alpha=0.6, marker='x', mew=2, ms=10, linewidth=0, label='measurements vs. m')\neb[-1][0].set_linewidth(1)\n"
           format port "pl.plot(sorted(ypos1+ypos2), pl.log(sorted(ypos1+ypos2))*(max(y0) / pl.log(max(ypos1+ypos2))), label='log(x)')\n"
           format port "pl.plot(sorted(ypos1+ypos2), pl.sqrt(sorted(ypos1+ypos2))*(max(y0) / pl.sqrt(max(ypos1+ypos2))), label='sqrt(x)')\n"
           format port "pl.plot(sorted(ypos1+ypos2), pl.multiply(sorted(ypos1+ypos2), max(y0) / max(ypos1+ypos2)), label='x')\n"
@@ -296,19 +300,25 @@ scalarMap = mpl.cm.ScalarMappable(norm=c
                           let : (offset (/ (apply max (append y⁰ y-opt)) 2)) (spreading (/ (apply max (append y⁰ y-opt)) (- (apply max member) (apply min member))))
                               format port "pl.plot(~A, ~A, marker='.', color=scalarMap.to_rgba(~A), linewidth=0, label='', alpha=0.6, zorder=-1)\n"
                                           . (/ step 1) (+ offset (* spreading (list-ref member param-idx))) param-idx
-          format port "pl.legend(loc='upper left')\n"
+          format port "pl.legend(loc='upper left', fancybox=True, framealpha=0.5)\n"
           format port "pl.xlabel('position [arbitrary units]')\n"
           format port "pl.ylabel('value [arbitrary units]')\n"
           format port "pl.title('~A')\n" "Operation scaling behaviour"
           format port "pl.xscale('log')\n"
           ;; format port "pl.yscale('log')\n"
-          format port "pl.show()\n"
+          if filename
+              format port "pl.savefig('~A')\n" filename
+              format port "pl.show()\n"
           format port "exit()\n"
           close-pipe port
 
 
 define : main args
    let*
-      : bench : benchmark-list-append ;; benchmark results
-        H : lambda (x pos) (H-N-m x pos #:const #t #:ON #t #:ONlogN #t)
-      plot-benchmark-result bench H
+      : H : lambda (x pos) (H-N-m x pos #:const #t #:ON #t #:ONlogN #t)
+        steps 200
+        pbr plot-benchmark-result
+      pbr (benchmark-list-append-N-1 steps) H #:filename "/tmp/benchmark-N-1.pdf"
+      pbr (benchmark-list-append-N-100 steps) H #:filename "/tmp/benchmark-N-100.pdf"
+      pbr (benchmark-list-append-1-m steps) H #:filename "/tmp/benchmark-1-m.pdf"
+      pbr (benchmark-list-append-100-m steps) H #:filename "/tmp/benchmark-100-m.pdf"