;;;   mat.scm : Matlab/Octave-like matrix operations for Gauche Scheme

;;;   Copyright (c) 2007 Issac Trotts, All rights reserved.
;;;   
;;;   Redistribution and use in source and binary forms, with or without
;;;   modification, are permitted provided that the following conditions
;;;   are met:
;;;   
;;;   1. Redistributions of source code must retain the above copyright
;;;      notice, this list of conditions and the following disclaimer.
;;;  
;;;   2. Redistributions in binary form must reproduce the above copyright
;;;      notice, this list of conditions and the following disclaimer in the
;;;      documentation and/or other materials provided with the distribution.
;;;  
;;;   3. Neither the name of the authors nor the names of its contributors
;;;      may be used to endorse or promote products derived from this
;;;      software without specific prior written permission.
;;;  
;;;   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
;;;   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
;;;   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
;;;   A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
;;;   OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
;;;   SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
;;;   TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
;;;   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
;;;   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
;;;   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
;;;   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

; History
; =======
; mar 14, 2007
;    Fixed a bug that made (randn) with no args not work.
;    Got rid of the overloading for * because it caused bugs like
;     (* 0.5 0.0) => 0.5.  Note: the bug was in Gauche 0.8.8, but is
;     fixed in 0.8.9 according to Shiro.
;    Added linspace and plot.
;
; mar 13, 2007
;    Added =! syntax so we can say (=! (m n) (size A)) in some contexts.
;    Added <- function for more concise assignment to entries in 
;    matrices.
;    Fixed a bug in matrix-reduce that was causing 1x1 matrices to be
;    returned for column vectors.
;    Output of matrices now shows less precision for greater readability.
;
; earlier
;    vectorized math functions (cos A) etc.

; Todo
; ====
; - Bridge with Octave.  
; - Slicing, for example (A ': 2) to get column 2 
; - Complex matrices.  This could be done by making a pair of
;   f64arrays, one for the real part and one for the complex part.
; - Greatly increase the performance somehow.  Currently Octave is
;   much faster.  For example try (mean (randn 1000)).
; - Block matrix notation
; - Look into the possibility of automatically converting 1x1 matrices
;   to ordinary numbers.

(define *image-viewer* "eog")
(define *numsize* 7)  ; use this to adjust the length of printed numbers


(use math.mt-random) ; better than sys-random?
(use math.const)     ; for pi
(use gauche.array)
(use srfi-1)  ; iota
(use srfi-11) ; let-values
(use gauche.test)
(use gauche.collection)

(define-macro (assert e)
  `(if (not ,e)
       (error "Assertion failed: " ,(x->string e))))


;; Thanks to Alex Shinn for posting this improved version of define-values
;; that makes =! work better (though still not perfectly).
;;
(define-syntax define-values
 (syntax-rules ()
   ((_ (var  ...) expr)
    (define-values-sub () (var ...) () expr))
   ((_ . else)
    (syntax-error "malformed define-values" (define-values . else)))
   ))

(define-syntax define-values-sub
 (syntax-rules ()
   ((_ (tmp ...) (last-var) (var ...) expr)
    ;; place the last var at the start of the list
    (define-values-sub (tmp1 tmp ...) () (last-var var ...) expr))
   ((_ (last-tmp tmp ...) () (last-var var ...) expr)
    (begin (define var #f) ... ;(undefined)
           (define last-var
             (receive (tmp ... last-tmp) expr
               (set! var tmp) ...
               last-tmp))))
   ((_ (tmp ...) (v v2 ...) (var ...) expr)
    (define-values-sub (tmp ... tmp1) (v2 ...) (var ... v) expr))
   ))

(define-syntax =!
  (syntax-rules ()
    ((=! (var1 . rest) val-maker)
     (define-values (var1 . rest) val-maker))
    ((=! var val-maker)
     (define var (let-values ((vals val-maker))
                             (if (null? (cdr vals))
                                 (car vals)
                                 vals))))))

;; Now we can write (<- A 1 2  4.0) to get A(1,2)=4.0.  Compare it to
;;                  (array-set! A 1 2 4.0) 
;;
;; and (<- v 2 5.0) to get v(2)=5.0.  Compare it to
;;     (array-set! v 2 0 5.0)  
;;
(define (<- mat i . args)
  (if (= (length args) 1)
      (let-values (((m n) (size mat))
                   (val (car args)))
                  (cond
                   ((= m 1) (array-set! mat 0 i val))
                   ((= n 1) (array-set! mat i 0 val))
                   (else
                    (error "Matrix must be Mx1 or 1xN to use (=! (mat i) x)")
                    )))
      (let ((j (car args))
            (val (cadr args)))
        (array-set! mat i j val)))
  mat)


(let ()
  (=! (a b c) (values 1 2 3))
  (assert (= a 1))
  (assert (= b 2))
  (assert (= c 3)))

;; Unfortunately we cannot do 
;;
;; (let ()
;;   (=! (a b) (values 1 2))
;;   (print a)
;;   (=! (c b) (values 3 4)))
;;
;; The =!'s must always occur at the top of the block enclosing them.


(define <matrix> <f64array>)

(define (matrix? thing)
  (and (is-a? thing <matrix>)
       (= (array-rank thing) 2)))

; Now you can read the i j element of mat by saying (mat i j).
(define-method object-apply ((A <matrix>) (i <integer>) (j <integer>))
  (check-arg matrix? A)
  (array-ref A i j))

;; (A ': 1) to get column 1 of A, like A(:,1) in Octave
;; This does not work
;;(define-method object-apply ((A <matrix>) (colon <symbol>) (j <integer>))
;;  (assert (eqv? colon ':))
;;  (matrix-column A j))

;; (A 1 ':) to get row 1 of A, like A(1,:) in Octave
;; This does not work
;;(define-method object-apply ((A <matrix>) (i <integer>) (colon <symbol>))
;;  (assert (eqv? colon ':))
;;  (matrix-row A i))

(define-method object-apply ((v <matrix>) (i <integer>))
  (check-arg matrix? v)
  (let-values (((m n) (size-values v)))
              (cond
               ((= m 1)
                (array-ref v 0 i))
               ((= n 1)
                (array-ref v i 0))
               (else
                (error "Cannot use a single index for a non-vector matrix")))))




;;;; Printing 

; Make the REPL pretty-print matrices

  
(define-method write-object ((A <matrix>) port)
  (matrix-print A port))

(define (matrix-print A port)
  (=! (m n) (size A))
  (dotimes (i m)
           (dotimes (j n)
                    (let1 aij (array-ref A i j)
                          (display (number->string-sized aij *numsize*) port)
                          ;(display aij)
                          (display (string-copy " ") port)))
           (newline port)))

;; Convert the given number to a string, trying to make it exactly
;; len characters long.  It may end up being longer.
;;
(define (number->string-sized num len)
  (let* ((nstr (number->string num))
         (slen (string-length nstr))
         (space #\ ))
    
    (cond
     ((= slen len) nstr)
     ((< slen len) (string-append (make-string (- len slen) space)
                                  nstr))
     ((> slen len)
      
      (let1 pos (string-scan nstr #\.)
            (if pos
                (let1 parts (string-split nstr #/\.|e/i)
                      (if (>= (length parts) 2)
                          (let* ((wholepart (car parts))
                                 (fracpart (cadr parts))
                                 (fraclen (string-length fracpart)))
                            (string-append
                             wholepart
                             "."
                             ;; FIXME: we should round here
                             (substring fracpart 0
                                        (max 0 (- fraclen
                                                  (- slen len))))
                             (if (string-scan nstr #\e)
                                 (string-append "e" (caddr parts))
                                 "")))
                          nstr))
                nstr  ; failed to shrink. maybe we should use sci notation here
                ))))))

(use gauche.array)

;(define old-* *)

;(define-method * ((a <number>) . rest)
;  (let-values (((nums mats) (partition number? 
;  (old-* a b))

; matrix multiplication 
(define m* array-mul)

; The following define-method led to bad results (* 1.0 0.0 ) => 1.0.
; use m* instead.

; (define-method * ((a <f64array>) . rest)
;   (print "is this on?")
;   (check-arg matrix? a)  ; Have to check the dims
;   (dolist (b rest)
;           (check-arg matrix? b))
;   (check-arg matrix? b)
;   (apply array-mul a rest))


; element-wise multiplication
(define-method .* ((A <matrix>) . rest)
  (check-arg matrix? A)
  (dolist (B rest)
          (check-arg matrix? B))
  (apply array-mul-elements A rest))

(define-method .* ((x <number>) . rest)
  (dolist (B rest)
          (check-arg number? B))
  (apply * x rest))

; element-wise division
(define-method ./ ((A <matrix>) . rest)
  (check-arg matrix? A)
  (dolist (B rest)
          (check-arg matrix? B))
  (if (null? rest)
      (map / A)
      (apply array-div-elements A rest)))

(define-method ./ ((x <number>) . rest)
  (dolist (B rest)
          (check-arg number? B))
  (apply / x rest))

; element-wise exponentiation
(define-method .^ ((A <matrix>) (e <number>))
  (check-arg matrix? A)
  (let-values (((m n) (size-values A)))
              (fun->matrix m n (lambda (i j)
                                 (expt (array-ref A i j) e)))))

(define-method .^ ((x <number>) (e <number>))
  (expt x e))


(define (mscale s A)
  (=! (m n) (size A))
  (fun->matrix m n
               (lambda (i j)
                 (* s (array-ref A i j)))))


;; The following define-method code led to trouble: (* 0.5 0.0) => 0.5
;; Use mscale instead

; scale a matrix by a scalar on the left
; (define-method s* ((s <number>) (A <matrix>))
;  (check-arg matrix? A)
;  (let-values (((m n) (size-values A)))
;              (fun->matrix m n
;                           (lambda (i j)
;                             (* s (array-ref A i j))))))

(define-method / ((A <matrix>) (s <number>))
  (check-arg matrix? A)
  (* (/ s) A))
         
(define-method + ((A <matrix>) . rest)
  (check-arg matrix? A)
  (dolist (B rest)
          (check-arg matrix? B))
  (apply array-add-elements A rest))

(define-method - ((A <matrix>) . rest)
  (check-arg matrix? A)
  (dolist (B rest)
          (check-arg matrix? B))
  (if (null? rest)
      (map - A)
      (apply array-sub-elements A rest)))

(define-method ^ ((A <matrix>) (n <integer>))
  (if (= n 0)
      (eye A)
      (* A (^ A (- n 1)))))

(define-method length ((A <matrix>))
  (check-arg matrix? A)
  (let-values (((m n) (size-values A)))
              m))

(define (matrix-negate A)
  (map - A))

(define-syntax domatrix
  (syntax-rules ()
    ((domatrix (i j A . retval) . body)
     (let-values (((m n) (size-values A)))
                 (dotimes (i m)
                          (dotimes (j n)
                                   . body))
                 . retval))))

(define-method map (f (A <matrix>))
  (let1 ret (zeros A)
        (domatrix (i j A  ret)
                  (let* ((aij (array-ref A i j))
                         (faij (f aij)))
                  (mset! ret i j faij)))))


;; This is dirty.  Is there a better way?
(define-macro (vectorize f)
  (let1 oldname (string->symbol (string-append "__" (x->string f)))
        `(begin
           (define ,oldname ,f)
           (define-method ,f ((A <matrix>))
             (map ,oldname A)))))

;; These lines allow expressions such as (tan A) like in Octave.
;;
(vectorize exp)
(vectorize log)
(vectorize sqrt)
(vectorize cos)
(vectorize sin)
(vectorize tan)
(vectorize asin)
(vectorize acos)
(vectorize atan)
        
(define (uniform-matrix m n fill)
  (let1 ret (make-f64array (shape 0 m 0 n))
        (domatrix (i j ret  ret)
                  (array-set! ret i j fill))))

(define (mset! mat i j val)
  (array-set! mat i j val))

(define (eye n)
  (identity-array n))

(define-method eye ((A <matrix>))
  (=! (m n) (size A))
  (assert (= m n))
  (eye m))

(define (zeros . args)
  (let-optionals* args ((m 1)
                        (n m))
                  (uniform-matrix m n 0.0)))

; Make a matrix of zeros the same size as A.
(define-method zeros ((A <matrix>))
  (let-values (((m n) (size-values A)))
              (zeros m n)))
                  
(define (ones . args)
  (let-optionals* args ((m 1)
                        (n m))
                  (uniform-matrix m n 1.0)))

(define-method ones ((A <matrix>))
  (let-values (((m n) (size-values A)))
              (ones m n)))


; multiple return values with size info
(define (size-values mat)
  (check-arg matrix? mat)
  (let1 s (array-shape mat)
        (values (array-ref s 0 1)
                (array-ref s 1 1))))

(define size size-values)

; pair containing size info
(define (size-pair mat)
  (let-values (((m n) (size-values mat)))
              (cons m n)))

; matrix containing size info
(define (msize mat)
  (let-values (((m n) (size-values mat)))
              (list->matrix (list m n))))

(define (row->list mat i)
  (check-arg matrix? mat)
  (let1 ret '()
        (let-values (((m n) (size-values mat)))
                    (dotimes (j n)
                             (push! ret (array-ref mat i j))))
        (reverse ret)))

(define (matrix-row mat i)
  (=! (m n) (size mat))
  (let1 ret (zeros 1 n)
        (dotimes (j n)
                 (array-set! ret 0 j
                             (array-ref mat i j)))
        ret))

           
(define (matrix-col mat j)
  (=! (m n) (size mat))
  (let1 ret (zeros m 1)
        (dotimes (i m)
                 (array-set! ret i 0
                             (array-ref mat i j)))
        ret))


(define (matrix->list mat)
  (check-arg matrix? mat)
  (let1 ret '()
        (let-values (((m n) (size-values mat)))
                    (dotimes (i m)
                             (push! ret (row->list mat i))))
        (reverse ret)))

(define (list->matrix lists)
  (check-arg list? lists)
  (if (null? lists)
      (zeros 0 0)
      (if (not (list? (car lists)))
          (transpose (list->matrix (list lists)))
          (let* ((m (length lists))
                 (n (length (car lists)))
                 (ret (zeros m n)))
            (if (> (length lists) 1)
                ;; Check that all rows are the same length
                (assert (apply = (map length lists))))
            (do ((i 0 (+ i 1))
                 (lists lists (cdr lists)))
                ((null? lists) ret)
              (do ((j 0 (+ j 1))
                   (ls (car lists) (cdr ls)))
                  ((null? ls))
                (array-set! ret i j (car ls))))))))

;; (vec 1 2 3) makes a column vector [1 2 3]'
(define (vec . args)
  (dolist (x args)
          (check-arg number? x))
  (list->matrix args))

(define (vec? mat)
  (and (matrix? mat)
       (let-values (((m n) (size mat)))
                   (or (= m 1) (= n 1)))))

(define inv array-inverse)

(define det determinant)

;; I don't use array-transpose from Gauche/lib/matrix.scm here because
;; it returns a object of type <array>, which is not compatible with
;; the <f64arrays> being used in mat.scm.
;;
(define (transpose mat)
  (let-values (((m n) (size-values mat)))
              (let1 ret (zeros n m)
                    (domatrix (i j ret  ret)
                              (array-set! ret i j
                                          (array-ref mat j i))))))

(define-method T ((a <matrix>))
  (check-arg matrix? a)
  (array-transpose a))

(define (fun->matrix m n f)
  (let1 ret (zeros m n)
        (domatrix (i j ret  ret)
                  (array-set! ret i j (f i j)))))

(define (reshape m2 n2 mat)
  (check-arg integer? m2)
  (check-arg integer? n2)
  (check-arg matrix? mat)
  (let ((ret (zeros m2 n2))
        (i2 0)
        (j2 0))
    (domatrix (i j mat  ret)
              (let1 matij (array-ref mat i j)
                    (array-set! ret i2 j2 matij))
              (inc! j2)
              (when (= j2 n2)
                    (set! j2 0)
                    (inc! i2)))))

;; Same as foo(:) in Octave
(define (flatten mat)
  (let-values (((m n) (size-values mat)))
              (reshape (* m n) 1 mat)))

(define float exact->inexact)

(define (matrix-reduce f init mat)
  (let-values
   (((m n) (size-values mat)))
   (cond
    ((= m 1)
     ;; mat has only one row, so the return value is a number
     (let1 ret init 
           (dotimes (j n ret)
                    (set! ret (f ret (array-ref mat 0 j))))))

    ((= n 1)
     ;; mat is a single column, so return a number
     (let1 ret init 
           (dotimes (i m ret)
                    (set! ret (f ret (array-ref mat i 0))))))

    (else 
       ;; multiple rows, so the return value is a row vector
       (let1 ret (zeros 1 n)
             (dotimes (j n)
                      (let1 acc init
                            (dotimes (i m)
                                     (set! acc
                                           (f acc (array-ref mat i j))))
                            (array-set! ret 0 j acc)))
             ret)))))

(define (sum mat)
  (matrix-reduce + 0.0 mat))

(define (prod mat)
  (matrix-reduce * 1.0 mat))

(define-method max ((mat <matrix>))
  (matrix-reduce max (array-ref mat 0 0) mat))

(define-method min ((mat <matrix>))
  (matrix-reduce min (array-ref mat 0 0) mat))

(define __rng (make <mersenne-twister> :seed (sys-time)))

; (define (rand1)
;   (float
;    (/ (sys-random)
;       (+ RAND_MAX 1))))

; rand : a single random number 0<=x<1
(define (rand1)
   (mt-random-real __rng))

; rand : uniform random numbers 0<=x<1)
(define (rand . args)
  (let-optionals* args ((m 1)
                        (n m))
                  (if (= m n 1)
                      (rand1)
                      (fun->matrix m n (lambda (i j) (rand1))))))
                  
; irand (random integers) 0<=i<max
(define (irand max . options)
  (let-keywords* options ((m :m 1)
                          (n :n m))
                 (fun->matrix m n (lambda (i j)
                                    (modulo (sys-random) max)))))

;; randn : Gaussian (normal) random numbers with mean zero and
;; unit variance, as described on p. 312 of David MacKay's book
(define (randn . options)
  (let-optionals* options ((m 1)
                           (n m))
                  (let ((u1 (rand m n))
                        (u2 (rand m n)))
                    (.* (cos (mscale (* 2.0 pi) u1))
                        (sqrt (mscale 2.0 (log (./ u2))))))))
                    

;; Save an image of matrix A in PGM format
(define (matrix-save-pgm A filename)
  (check-arg matrix? A)
  (check-arg string? filename)
  (let ((mx (max (max A)))
        (mn (min (min A))))
    (let-values (((m n) (size-values A)))
                (with-output-to-file filename
                  (lambda ()
                    (print "P5")
                    (print n " " m)
                    (print "255")
                    (domatrix (i j A)
                              (let* ((Aij (array-ref A i j))
                                     (c (* 255 (/ (- Aij mn) (- mx mn))))
                                     (c (x->integer c)))
                                (write-byte c))))))))
                                  
(define (imagesc A)
  (check-arg matrix? A)
  (let1 pgmname (string-append (sys-tmpnam) ".pgm")
        (matrix-save-pgm A pgmname)
        (sys-system (string-append *image-viewer* " " pgmname " &"))
        #f))
        

; If X is a single row or column, return its mean.
; Otherwise, return the means of its columns.
;
(define (mean X)
  (check-arg matrix? X)
  (let-values (((m n) (size-values X)))
              (if (or (= m 1) (= n 1))
                  (/ (sum X) (* m n))
                  (/ (sum X) m))))

; variance
(define (var X)
  (- (mean (.^ X 2))
     (.^ (mean X) 2)))

; standard deviation
(define (std X)
  (sqrt (var X)))

; root mean square
(define rms std)

(define (copy X)
  (check-arg matrix? X)
  (let-values (((m n) (size-values X)))
              (fun->matrix m n
                           (lambda (i j)
                             (array-ref X i j)))))

;; Make vectors columnar.  Same as x(:) in Octave when x is a vector.
;;
(define (columnify x)
  (check-arg matrix? x)
  (let-values (((m n) (size-values x)))
              (assert (or (= m 1) (= n 1)))
              (if (= n 1)
                  x
                  (reshape n m x))))

;; optimization-friendly numerical gradient
(define (numerical-gradient f . options)
  (let-optionals*
   options ((step 1e-6))
   (lambda (x) 
     (let* ((fx (f x))
            (x (columnify x))
            (N (length x))
            (g (zeros x)))
        (dotimes (n N)
          (let ((x2 (copy x)))
            (array-set! x2 n 0 
                        (+ (array-ref x n 0) step))
            (let ((df_dxn (/ (- (f x2) fx) step)))
              (array-set! g n 0 df_dxn))))
        g))))

;; For example (define g (grad f))
;;
(define grad numerical-gradient)

(define (squared x)
  (* x x))

(define (cubed x)
  (* x x x))

(define (gradient-descent f x . options)
  (let-keywords*
   options ((g (numerical-gradient f))
            (step 1e-1)
            (ftol 1e-6)
            (gtol 1e-6)
            (max-steps 1000))
   (let ((fx (f x)))
     (if (or (= max-steps 0)
             (<= fx ftol))
         (values x max-steps)
         (let ((gx (g x)))
           (if (<= (std gx) gtol)
               (values x max-steps)
               (let ((x2 (- x (* step gx)))
                     (max-steps (- max-steps 1)))
                 (gradient-descent 
                  f x2
                  :step step
                  :ftol ftol
                  :gtol gtol
                  :max-steps max-steps))))))))
  

(define (test-gradient-descent)
  (test "gradient-descent" (vec 1 -2)
        (lambda ()
          (let1 x0 (vec 0 0)
                (gradient-descent 
                 (lambda (x)
                   (+ (squared (- (x 0) 1.0))
                      (squared (- (x 1) -2.0))))
                 x0)))
        (lambda (a b)
          (<= (std (- a b)) 1e-3))))
            

(define (gaussian mu sigma x)
  (assert (> sigma 0.0))
  (/ (exp (* 0.5 (expt (/ (- x mu) sigma) 2.0)))
     (* sigma (sqrt (* 2.0 pi)))))

;  (* 0.5 (expt (/ (- x mu) sigma) 2.0)))
       
  
(define (new-gaussian mu sigma)
  (lambda (msg . args)
    (case msg
      ((sample) (+ mu (* sigma (randn))))
      ((eval)
       (let1 x (car args)
             (gaussian mu sigma x))))))

(define (new-uniform lo hi)
  (lambda (msg . args)
    (case msg
      ((sample)
       (+ lo (* (- hi lo) (rand))))
      ((eval)
       (let1 x (car args)
             (if (or (< x lo)
                     (> x hi))
                 0.0
                 (/ (- hi lo))))))))
                     

;; constructor for a joint distribution made up of
;; independent factors
(define (new-indep-joint factor1 . rest)
  (let ((factors (cons factor1 rest)))
    (lambda (msg . args)
      (case msg
        ((sample)
         (list->matrix
          (map (lambda (factor)
                 (factor 'sample))
               factors)))
        ((case eval)
         (let1 x (car args)
               (reduce * 1.0
                       (map (lambda (factor i)
                              (factor 'eval (x i)))
                            factors
                            (iota (length factors))))))))))
       


;; 1D minimization using probabilities
;;First, let's see what distro we get for xmin when we sample
;;a bunch of as bs and cs.

(define p-a (new-gaussian 1.0 2.0))
(define p-b (new-gaussian 0.0 10.0))
(define p-c (new-gaussian 0.0 10.0))

;; A distribution over smiling quadratics 
(define p-quadratic (new-indep-joint p-a p-b p-c))

(define (rand-quadratic)
  (let1 coeffs (p-quadratic 'sample)
        (lambda (x)
          (case x
            ((coeffs) coeffs) ; extract the coeffs
            (else (+ (* (coeffs 0) x x) ; evaluate
                     (* (coeffs 1) x)
                     (coeffs 2)))))))

;; like in Octave
(define (linspace base limit N)
  (assert (> N 1))
  (let1 ret (zeros N 1)
    (dotimes (n N)
             (let1 x (+ base (* (- limit base) (/ n (- N 1))))
                   (array-set! ret n 0 x)))
    ret))


(define (plot . args)
  (define (gnuplot)
    (sys-system "echo 'plot \"-data-\" with lines' | gnuplot -persist"))
  (cond
   ((= (length args) 1)
    (let1 y (car args)
          (with-output-to-file "-data-"
            (lambda ()
              (print (columnify y))))
          (gnuplot)))
   ((= (length args) 2)
    (let ((x (car args))
          (y (cadr args)))
      (assert (vec? x))
      (assert (vec? y))
      (let ((xlen (length x))
            (ylen (length y)))
        (assert (= xlen ylen))
        (with-output-to-file "-data-"
          (lambda ()
            (dotimes (i xlen)
                     (print (x i) " " (y i)))))
        (gnuplot))))))

(define (plot-test)
  (=! x (linspace (- pi) pi 100))
  (plot x (sin x)))

; MacKay style gradient descent

; diag


; hist 

; meshgrid

; surf

; cumsum

; eig

; pca

; x->matrix : convert lists and vectors to matrices