(in-package :ll-test)

(defun least-square-pow-and-coeff (x y)
  "Fits to y = a * x ^ b, returns (b a)"
  (let* ((logx (.log x))
	 (logy (.log y))
	 (Sx2 (msum (.^ logx 2)))
	 (Sx (msum logx))
	 (Sy (msum logy))
	 (Sxy (msum (.* logx logy)))
	 (N (size logx))
	 (a (./ (.- (.* N Sxy) (.* Sx Sy))
		(.- (.* N Sx2) (.^ Sx 2))))
	 (b (./ (.- Sy (.* a Sx)) N)))
    (list a (.exp b))))
    
(defun test-power-law ()
  (let* ((x (drange 10 1 10))
	 (y (.* 6.3 (.^ x 3.4))))
    (least-square-pow-and-coeff x y)))

(defun seconds-of-internal-real-time (fun)
  "A primitive timer"
  (let ((start (get-internal-real-time)))
    (funcall fun)
    (let ((end (get-internal-real-time)))
      (/ (- end start) 1.0 INTERNAL-TIME-UNITS-PER-SECOND))))

(defun get-times (fun inputs num)
  (mapcar (lambda (x num) 
	    (prog1 
		(/ (seconds-of-internal-real-time 
		    (lambda ()
		      (dotimes (i num)
			(funcall fun x))))
		   num)
	      (sb-ext:gc :full t)))
	  inputs num))

(defun calc-scaling-time (fun matlist Ns &optional (num (mapcar (constantly 1) Ns)))
  "Returns the scaling power and coefficient"
  (let* ((y (get-times fun matlist num)))
    (least-square-pow-and-coeff (dmat (list Ns))
				(dmat (list y)))))


(defclass test-scaling-time (test-base)
  ((size :initarg size :initform nil)
   (arg :initarg arg)
   (times :initarg :times) 
   (a :initarg :a)
   (b :initarg :b)))

(defmethod initialize-instance :after ((test test-scaling-time) 
				       &key (guess-b 2))
  (with-slots (times size) test
    (setf times (mapcar (lambda (N) 
			  (max 1 (floor (.^ (/ 1000 N) guess-b))))
			size))))



  
(defun test-scaling ()
  (let* ((x (loop for i from 100 below 500 by 100 collect i))
	 (Ns (mapcar (lambda (x) (max 1 (floor (.^ (/ 1000 x) 2)))) x))
	 (args (mapcar (lambda (N) (drandom N N)) x )))
    (calc-scaling-time (lambda (x) (.^ x -2)) args x Ns)))