;;
;; Tensor Arithmetics
;;
(+ 1 [| 1 2 3 |])
;=>[|2 3 4|]

(+ [| 1 2 3 |] 1)
;=>[|2 3 4|]

(+ [| 1 2 3 |]_i [| 1 2 3 |]_i)
;=>[|2 4 6|]_i

(+ [| 10 20 30 |] [| 1 2 3 |])
;=>[| [| 11 12 13 |] [| 21 22 23 |] [| 31 32 33 |] |]

(+ [| 100 200 300 |]_i
   [|[| 1 2 3 |]
     [| 10 20 30 |]|]_j_i)
;=>[| [| 101 110 |] [| 202 220 |] [| 303 330 |] |]_i_j

(+ [|[| 11 12 |]
     [| 21 22 |]
     [| 31 32 |]|]_i_j
   [| 100 200 300 |]_i)
;=>[| [| 111 112 |] [| 221 222 |] [| 331 332 |] |]_i_j

(+ [| 100 200 300 |]_i
   [|[| 11 12 |]
     [| 21 22 |]
     [| 31 32 |]|]_i_j)
;=>[| [| 111 112 |] [| 221 222 |] [| 331 332 |] |]_i_j

;;
;; Derivative
;;
(∂/∂ (f x y z) x)
;=>(f_1 x y z)

(∂/∂ [| (f x) (g x) |] x)
;=>[| (f_1 x) (g_1 x) |]

(∂/∂ (f x y z) [| x y z |])
;=>[| (f_1 x y z) (f_2 x y z) (f_3 x y z) |]

([| (∂/∂ $ x) (∂/∂ $ y) |] (f x y))
;=>[| (f_1 x y) (f_2 x y) |]

([| (∂/∂ $ x) (∂/∂ $ y) |] [| (f x y) (g x y) |])
;=>[| [| (f_1 x y) (g_1 x y) |] [| (f_2 x y) (g_2 x y) |] |]

;;
;; Nabla
;;
(define $∇ ∂/∂)

(∇ (f x y) [| x y |])
;=>[| (f_1 x y) (f_2 x y) |]

(∇ [| (f x y) (g x y) |] [| x y |])
;=>[| [| (f_1 x y) (f_2 x y) |] [| (g_1 x y) (g_2 x y) |] |]

;;
;; Contraction
;;
(contract + (* [|1 2 3|]~i [|10 20 30|]_i))
;=>
140

(define $trace (lambda [%t] (with-symbols {i} (contract + t~i_i))))

(trace [|[|10 20 30|] [|20 40 60|] [|30 60 90|]|])
;=>
140

;;
;; Divergence
;;
(define $div (compose ∇ (trace $)))

(div [| (f x y z) (g x y z) (h x y z) |] [| x y z |])
;=>(+ (f_1 x y z) (g_2 x y z) (h_3 x y z))

;;
;; Taylor Expansion
;;
(define $taylor-expansion
  (lambda [%f %xs %as]
    (with-symbols {h}
      (let {[$hs (generate-tensor 1#h_%1 (tensor-size xs))]}
        (map2 *
              (map 1#(/ 1 (fact %1)) nats0)
              (map (compose (V.substitute xs as $)
                            (V.substitute hs (with-symbols {i} (- xs_i as_i)) $))
                   (iterate (compose (∇ $ xs) (V.* hs $)) f)))))))

(take 3 (taylor-expansion (f x) x 0))
;=>
;{(f 0)
; (* x (f_1 0))
; (/ (* x^2 (f_1_1 0))
;    2)}

(take 3 (taylor-expansion (f x y) [| x y |] [| 0 0 |]))
;=>
;{(f 0 0)
; (+ (* x (f_1 0 0))
;    (* y (f_2 0 0)))
; (/ (+ (* x^2 (f_1_1 0 0))
;       (* x y (f_2_1 0 0))
;       (* y x (f_1_2 0 0))
;       (* y^2 (f_2_2 0 0)))
;    2)}