;    Copyright (C) 2006  Will M. Farr <farr@mit.edu>
;
;    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.


(module manifolds "generics.ss"
  (require "derivatives.ss" "tuples.ss" "operators.ss"
           (all-except (lib "43.ss" "srfi") vector-fill! vector->list))
  (require-for-syntax (all-except (lib "43.ss" "srfi") vector-fill! vector->list)
                      (only (lib "1.ss" "srfi") any))
  
  (provide <chart> chart? <manifold> manifold? <vector-field> vector-field? <one-form-field> one-form-field?
           chart->coordinate-vector-fields chart->coordinate-one-form-fields
           procedure->vector-field procedure->one-form-field
           component-field->vector-field component-field->one-form-field
           vector-field->component-field one-form-field->component-field
           define-named-coordinates
           lie-bracket)
  
  (defclass <chart> ()
    chi chiinv dim
    :auto #t
    :printer #t)
  
  (define (chart->coordinate-vector-fields chart)
    (with-slots chart (dim)
      (down-unfold 
       (lambda (i)
         (component-field->vector-field 
          (lambda (m)
            (up-unfold (lambda (j) (if (= i j) 1 0)) dim))
          chart))
       dim)))
  
  (define (chart->coordinate-one-form-fields chart)
    (with-slots chart (dim)
      (up-unfold 
       (lambda (i)
         (component-field->one-form-field
          (lambda (m)
            (down-unfold (lambda (j) (if (= i j) 1 0)) dim))
          chart))
       dim)))
  
  (defclass <manifold> () 
    standard-chart
    :auto #t
    :printer #t)
  
  (defentityclass <vector-field> (<linear-operator>)
    :auto #t
    :printer #t)
  
  (define (procedure->vector-field proc)
    (let ((result (make <vector-field>)))
      (set-instance-proc! result
                          proc)
      result)) 
  
  (define (component-field->vector-field comp chart)
    (with-slots chart (chi chiinv)
      (define ((proc f) m)
        (let ((comp (comp m)))
          (* comp
             ((D (compose f chiinv))
              (chi m)))))
      (procedure->vector-field proc)))
  
  (define ((vector-field->component-field v chart) m)
    (with-slots chart (chi chiinv)
      ((v chi) m)))
      
  (defmethod (mul (f <function>) (v <vector-field>))
    (define ((proc g) m)
      (* (f m)
         ((v g) m)))
    (procedure->vector-field proc))
  (defmethod (mul (x <number>) (v <vector-field>))
    (mul (lambda (m) x) v))
  (defmethod (mul (v <vector-field>) (f <function>))
    (mul f v))
  (defmethod (mul (v <vector-field>) (x <number>))
    (mul x v))

  (defmethod (div (v <vector-field>) (f <function>))
    (define ((proc g) m)
      (/ ((v g) m)
         (f m)))
    (procedure->vector-field proc))
  (defmethod (div (v <vector-field>) (x <number>))
    (div v (lambda (m) x)))  
  
  ;; Want to define addition/subtraction of vector fields differently from
  ;; linear operators.  Linear operator addition/subtraction of vector fields
  ;; involves addition and subtraction of functions because the vector fields 
  ;; are operators from functions -> functions.  I don't want to define
  ;; addition and subtraction of functions generically because it's not appropriate 
  ;; for all functions.  So, I define addition and subtraction of vector fields 
  ;; here as a special case.
  (defmethod (add (v1 <vector-field>) (v2 <vector-field>))
    (define ((proc f) m)
      (+ ((v1 f) m)
         ((v2 f) m)))
    (procedure->vector-field proc))
  (defmethod (minus (v <vector-field>))
    (define ((proc f) m)
      (- ((v f) m)))
    (procedure->vector-field proc))
  (defmethod (sub (v1 <vector-field>) (v2 <vector-field>))
    (define ((proc f) m)
      (- ((v1 f) m)
         ((v2 f) m)))
    (procedure->vector-field proc))
  
  (defentityclass <one-form-field> (<linear-operator>)
    :auto #t
    :printer #t)
  
  (define (procedure->one-form-field proc)
    (let ((result (make <one-form-field>)))
      (set-instance-proc! result proc)
      result))
  
  (define (component-field->one-form-field comp chart)
    (with-slots chart (chi chiinv)
      (define ((proc v) m)
        (* (comp m)
           ((v chi) m)))
      (procedure->one-form-field proc)))
  
  (define (one-form-field->component-field w chart)
    (let ((coord-vectors (chart->coordinate-vector-fields chart)))
      (with-slots chart (dim)
        (lambda (m)
          (down-unfold 
           (lambda (i)
             ((w (ref coord-vectors i)) m))
           dim)))))
    

  (defmethod (times (w <one-form-field>)) w)
  (defmethod (mul (w <one-form-field>) (f <function>))
    (mul f w))
  (defmethod (mul (w <one-form-field>) (x <number>))
    (mul x w))
  (defmethod (mul (f <function>) (w <one-form-field>))
    (define ((proc v) m)
      (* (f m)
         ((w v) m)))
    (procedure->one-form-field proc))
  (defmethod (mul (x <number>) (w <one-form-field>))
    (mul (lambda (m) x) w))
  (defmethod (div (w <one-form-field>) (f <function>))
    (define ((proc v) m)
      (/ ((w v) m)
         (f m)))
    (procedure->one-form-field proc))
  (defmethod (div (w <one-form-field>) (x <number>))
    (div w (lambda (m) x)))
  
  (define-for-syntax (mapi fn . lists)
    (let loop ((i 0) (result '()) (lists lists))
      (if (any null? lists)
          (reverse result)
          (loop (+ i 1) 
                (cons (apply fn i (map car lists))
                      result)
                (map cdr lists)))))
  
  (define-syntax define-named-coordinates
    (lambda (stx)
      (syntax-case stx ()
        ((define-named-coordinates (x y ...) chart)
         (let ((vectors (map (lambda (coord)
                               (datum->syntax-object
                                #'x
                                (symbol-append 'd/d coord)
                                #'x))
                             (syntax-object->datum #'(x y ...))))
               (forms (map (lambda (coord)
                             (datum->syntax-object
                              #'x
                              (symbol-append 'd coord)
                              #'x))
                           (syntax-object->datum #'(x y ...)))))
           (let ((nstx (datum->syntax-object #'x (length vectors) #'x)))
             #`(define-values (x y ... #,@vectors #,@forms)
                 (let ((c chart))
                   (with-slots c (chi chiinv)
                     (apply 
                      values
                      #,@(mapi (lambda (i dummy)
                                 (let ((istx (datum->syntax-object #'x i #'x)))
                                   #`(lambda (m)
                                       (ref (chi m) #,i))))
                               vectors)
                      (append 
                       (vector->list 
                        (tuple->vector 
                         (chart->coordinate-vector-fields c)))
                       (vector->list
                        (tuple->vector
                         (chart->coordinate-one-form-fields c))))))))))))))
  
  (define (lie-bracket v1 v2)
    (define ((proc f) m)
      (- ((v1 (v2 f)) m)
         ((v2 (v1 f)) m)))
    (procedure->vector-field proc)))