diff.clj

;; Symbolic differentiation
;;
;; translated from aprolog/examples/diff.apl
;; http://homepages.inf.ed.ac.uk/jcheney/programs/aprolog/examples/diff.apl
;;

(ns nominal.diff
  (:refer-clojure :exclude [==])
  (:use [clojure.core.logic :exclude [is] :as l]
        [clojure.core.logic.nominal :exclude [fresh hash] :as nom])
  (:use [clojure.test])
  (:require [clojure.pprint :as pp]))

;; TODO(namin): Should this goal be part of core.logic.nominal?
(defn swapo [a b t out]
  (fn [s]
    (let [[t s] (swap-noms t [a b] s)]
      (bind s (== t out)))))

(defn numbero [x]
  (predc x number? `number?))

(defn nomo [x]
  (predc x nom? `nomo?))

(defn diffo [x e ep]
  (conde
    [(numbero e) (== 0 ep)]
    [(nomo e) (== x e) (== 1 ep)]
    [(nomo e) (!= x e) (== e ep)]
    [(fresh [e1 e2 e1p e2p]
       (== ['+ e1 e2] e)
       (== ['+ e1p e2p] ep)
       (diffo x e1 e1p)
       (diffo x e2 e2p))]
    [(fresh [e1 e2 e1p e2p]
       (== ['* e1 e2] e)
       (== ['+ ['* e1p e2] ['* e1 e2p]] ep)
       (diffo x e1 e1p)
       (diffo x e2 e2p))]
    [(fresh [e1 e1p]
       (== ['sin e1] e)
       (== ['* ['cos e1] e1p] ep)
       (diffo x e1 e1p))]
    [(fresh [e1 e1p]
       (== ['cos e1] e)
       (== ['* ['- ['sin e1]] e1p] ep)
       (diffo x e1 e1p))]
    [(fresh [e1 e1p]
       (== ['neg e1] e)
       (== ['neg e1p] ep)
       (diffo x e1 e1p))]
    ;; fundamental theorem of calculus, special case
    [(fresh [e1 low high]
       (nom/fresh [y]
         (== ['integ (nom/tie y e1) low high] e)
         (numbero low)
         (== x high)
         (conde
           [(== x y) (== e1 ep)]
           [(!= x y)
            (nom/hash x e1)
            (swapo x y e1 ep)])))]))

(is (= (run* [q]
         (nom/fresh [x]
           (diffo x ['* x x] q)))
      '([+ [* 1 a_0] [* a_0 1]])))

(is (= (run* [q]
         (nom/fresh [x]
           (diffo x ['integ (nom/tie x x) 0 x] q)))
      '(a_0)))

(is (= (run* [q]
         (nom/fresh [x y]
           (diffo x ['integ (nom/tie y y) 0 x] q)))
      '(a_0)))

;; This is actually 2x according to Wolfram Alpha, but is not covered by our simplified rules.
;; http://www.wolframalpha.com/input/?i=d%2Fdx+%28int+x+dy%2C+y%3D0..x%29
(is (= (run* [q]
         (nom/fresh [x y]
           (diffo x ['integ (nom/tie y x) 0 x] q)))
      '()))