;;; Lisplab, level2-array-functions.lisp
;;; Level2, algbraic functions on arrays
;;; 
;;; TOOD: Make fast methods also for integers.


;;; Copyright (C) 2009 Joern Inge Vestgaarden
;;;
;;; This program is free software; you can redistribute it and/or modify
;;; it under the terms of the GNU General Public License as published by
;;; the Free Software Foundation; either version 2 of the License, or
;;; (at your option) any later version.
;;;
;;; This program is distributed in the hope that it will be useful,
;;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;;; GNU General Public License for more details.
;;;
;;; You should have received a copy of the GNU General Public License along
;;; with this program; if not, write to the Free Software Foundation, Inc.,
;;; 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

(in-package :lisplab) 

;;; Array specialization for the main algebraic operators and functions.
;;; It is obviously usefull to have .+, .-, etc. and .sin .cos, etc.
;;; to operate on arrays. But what about predicates: .=, .<, etc.  
;;; It might be an idea to do specialized these for arrays, but so far
;;; they are not.
 

(defmacro define-array-binary-operator (new old)
  ;; Creates methods for binary operators specialized on arrys.
  ;; This is not easy to do in an efficient way, which is one 
  ;; main reason to make specialized matrices (such as blas-real). 
  ;; However I have put in optimizations for double floats. It 
  ;; might be usefull to also do so for complex double-float and some 
  ;; integer types and bytes. 
  (let ((a (gensym))
	(b (gensym))
	(c (gensym))
	(i (gensym)))
    `(progn

       ;; two arrays
       (defmethod ,new ((,a array) (,b array))
	 (if (and (eql (element-type ,a) 'double-float)
		  (subtypep (type-of ,a) 'simple-array)
		  (eql (element-type ,b) 'double-float)
		  (subtypep (type-of ,b) 'simple-array))
	     (let ((,c (copy ,a)))
	       (declare (type (simple-array double-float) ,b ,c))
	       (dotimes (,i (min (size ,c) (size ,a)))
		 (setf (row-major-aref ,c ,i) 
		       (,old (row-major-aref ,c ,i) (row-major-aref ,b ,i))))
	       ,c)
	     (let ((,c (mcreate ,a 0)))
	       (dotimes (,i (min (size ,c) (size ,a)))
		 (setf (vref ,c ,i)
		       (,new (vref ,a ,i) (vref ,b ,i))))
	       ,c)))

       ;; array and number
       (defmethod ,new ((,a array) (,b number))
	 (if (and (eql (element-type ,a) 'double-float)
		  (subtypep (type-of ,a) 'simple-array)
		  (realp ,b))
	     (let ((,b (coerce ,b 'double-float))
		   (,c (copy ,a)))
	       (declare (type (simple-array double-float) ,c))
	       (dotimes (,i (size ,c))
		 (setf (row-major-aref ,c ,i) 
		       (,old (row-major-aref ,c ,i) ,b)))
	       ,c)
	     (let ((,c (mcreate ,a 0)))
	       (dotimes (,i (size ,c))
		 (setf (vref ,c ,i)
		       (,new (vref ,a ,i) ,b)))
	       ,c)))
      
       ;; number and array
       (defmethod ,new ((,a number) (,b array))
	 (if (and (eql (element-type ,b) 'double-float)
		  (subtypep (type-of ,b) 'simple-array)
		  (realp ,a))
	     (let ((,b (coerce ,a 'double-float))
		   (,c (copy ,b)))
	       (declare (type (simple-array double-float) ,c))
	       (dotimes (,i (size ,c))
		 (setf (row-major-aref ,c ,i)
		       (,old ,b (row-major-aref ,c ,i))))
	       ,c)
	     (let ((,c (mcreate ,b 0)))
	       (dotimes (,i (size ,c))
		 (setf (vref ,c ,i)
		       (,new ,b (vref ,b ,i))))
	       ,c))))))

(define-array-binary-operator .add +)
(define-array-binary-operator .sub -)
(define-array-binary-operator .mul *)
(define-array-binary-operator .div /)
(define-array-binary-operator .expt expt)

(defmacro each-array-element-df-to-df (x form) 
  "Applies a form on each element of an array. The form must 
make real output for real arguments"
  (let ((i (gensym))
	(y (gensym)))
    `(if (and (eql (element-type ,x) 'double-float)
	      (subtypep (type-of ,x) 'simple-array))
	 (let ((,y (copy ,x)))
	   (declare (type (simple-array double-float) ,y))
	   (dotimes (,i (size ,y))
	     (let ((,x (row-major-aref ,y ,i)))
	       (declare (type double-float ,x))	       
	       (setf (row-major-aref ,y ,i) ,form)))
	   ,y)
	 (let ((,y (mcreate ,x 0)))
	   (dotimes (,i (size ,y))
	     (let ((,x (vref ,x ,i)))
	       (setf (vref ,y ,i) ,form)))
	   ,y))))

(defmacro each-array-element-df-to-complex-df (x form) 
  "Applies a form on each element of an array. The form must 
make complex output for real arguments"
  (let ((i (gensym))
	(y (gensym)))
    `(if (and (eql (element-type ,x) 'double-float)
	      (subtypep (type-of ,x) 'simple-array))
	 (let ((,y (make-array (dim ,x) :element-type '(complex double-float))))
	   (declare (type (simple-array (complex double-float)) ,y))
	   (dotimes (,i (size ,y))
	     (let ((,x (row-major-aref ,x ,i)))
	       (declare (type double-float ,x))	       
	       (setf (row-major-aref ,y ,i) ,form)))
	   ,y)
	 (let ((,y (mcreate ,x 0))) ; TOOD must make sure to allow complex values
	   (dotimes (,i (size ,y))
	     (let ((,x (vref ,x ,i)))
	       (setf (vref ,y ,i) ,form)))
	   ,y))))

;;; Trignometric functions

(defmethod .sin ((x array))
  (each-array-element-df-to-df x (.sin x)))

(defmethod .cos ((x array))
  (each-array-element-df-to-df x (.cos x)))

(defmethod .tan ((x array))
  (each-array-element-df-to-df x (.tan x)))

(defmethod .asin ((x array))
  (each-array-element-df-to-df x (.asin x)))

(defmethod .acos ((x array))
  (each-array-element-df-to-df x (.acos x)))

(defmethod .atan ((x array))
  (each-array-element-df-to-df x (.atan x)))


;;; Hyperbolic functions

(defmethod .sinh ((x array))
  (each-array-element-df-to-df x (.sinh x)))

(defmethod .cosh ((x array))
  (each-array-element-df-to-df x (.cosh x)))

(defmethod .tanh ((x array))
  (each-array-element-df-to-df x (.tanh x)))

(defmethod .asinh ((x array))
  (each-array-element-df-to-df x (.asinh x)))

(defmethod .acosh ((x array))
  (each-array-element-df-to-df x (.acosh x)))

(defmethod .atanh ((x array))
  (each-array-element-df-to-df x (.atanh x)))

;;; Log and exponent

(defmethod .log ((x array) &optional base)  
  (each-array-element-df-to-df x (.log x base)))

(defmethod .exp ((x array))  
  (each-array-element-df-to-df x (.exp x)))

;;; Bessel functions

(defmethod .besj (n (x array))
  (each-array-element-df-to-df x (.besj n x)))
		       
(defmethod .besy (n (x array))
  (each-array-element-df-to-df x (.besy n x)))

(defmethod .besi (n (x array))
  (each-array-element-df-to-df x (.besi n x)))

(defmethod .besk (n (x array))
  (each-array-element-df-to-df x (.besk n x)))

;;; Hankel functions. NB! These are complex with real arguments.

(defmethod .besh1 (n (x array))
  (each-array-element-df-to-complex-df x (.besh1 n x)))

(defmethod .besh2 (n (x array))
  (each-array-element-df-to-complex-df x (.besh2 n x)))