adaptive.clj

(defn adaptive
  "Accepts two 'integrator' functions of:

  - `f`: some integrand
  - `a` and `b`: the lower and upper endpoints of integration
  - `opts`, a dictionary of configuration options

  And returns a new function of the same signature that adaptively subdivides
  the region $a, b$ into intervals if integration fails to converge."
  ([integrator] (adaptive integrator integrator))
  ([open-integrator closed-integrator]
   (fn rec
     ([f a b] (rec f a b {}))
     ([f a b {:keys [fuzz-factor]
              :or {fuzz-factor 0}
              :as opts}]

      ;; The user can wrap either `a` or `b` in markers that tag the endpoints
      ;; as explicitly open or closed; `(closed a) b`, for example, or `(open
      ;; a) (open b)`. An interval is `closed?` if both endpoints or closed and
      ;; treated as open otherwise.
      (letfn [(integrate [l r]
                (if (ui/closed? l r)
                  (closed-integrator f l r opts)
                  (open-integrator f l r opts)))]

        ;; This `(ua/kahan-sum)` is an aggregator that is a bit more accurate
        ;; for summing up floating point values; the simpler version would just
        ;; have `0.0` here and `sum` as the final return value.
        (loop [sum   (ua/kahan-sum)
               stack [[a b]]]

          ;; I COULD do this recursion explicitly, but wanted to keep an
          ;; explicit stack around in case the recursion dropped many steps.
          ;; This might be silly, since I'll descend by log_2 and unless I make
          ;; up some crazy-misbehaved function I won't blow the stack...
          (if (empty? stack)
            (first sum)
            (let [[l r]  (peek stack)
                  rest   (pop stack) 
                  {:keys [converged? result]} (integrate l r)]
              (if converged?
                (recur (ua/kahan-sum sum result) rest)

                ;; The midpoint is explicitly closed.
                (let [mid (ui/closed
                           (split-point l r fuzz-factor))]
                  (recur sum (conj rest [l mid] [mid r]))))))))))))

;; This version's a little tighter. No more explicit stack, and also no more
;; `kahan-sum` trick to hold off floating point addition errors.

(defn adaptive
  ([integrator] (adaptive integrator integrator))
  ([open-integrator closed-integrator]
   (letfn [(integrate [f l r opts]
             (if (ui/closed? l r)
               (closed-integrator f l r opts)
               (open-integrator f l r opts)))]
     (fn rec
       ([f a b] (rec f a b {}))
       ([f a b opts]
        (let [{:keys [converged? result]} (integrate f a b opts)]
          (if converged?
            result
            (let [mid (ui/closed
                       (split-point a b (:fuzz-factor opts 0)))]
              (+ (rec f a mid opts)
                 (rec f mid b opts))))))))))

full_attempt.cljc

(defn split-point
  "If `fuzz-factor` is non-zero, pick a point within $fuzz-factor \\over 2$ of the
  midpoint to split some interval.

  Else, returns the midpoint (ie `fuzz-factor` defaults to 0)."
  ([a b] (split-point a b 0))
  ([a b fuzz-factor]
   {:pre [(>= fuzz-factor 0)
          (< fuzz-factor 1)]}
   (let [width  (- b a)
         offset (if (zero? fuzz-factor)
                  0.5
                  (+ 0.5 (* fuzz-factor (dec (rand 2.0)))))]
     (+ a (* offset width)))))

;; Intervals... these no longer wrap the values themselves, but live alongside.

(def open        [::open ::open])
(def closed      [::closed ::closed])
(def open-closed [::open ::closed])
(def closed-open [::closed ::open])

(defn closed? [x] (= x closed))
(def open? (complement closed?))

(defn close-l [[_ r]] [::closed r])
(defn close-r [[l _]] [l ::closed])
(defn open-l [[_ r]] [::open r])
(defn open-r [[l _]] [l ::open])

(defn- fill-defaults [opts]
  (merge {:maxterms *integrate-n*
          :fuzz-factor *quadrature-neighborhood-width*
          :interval open}
         opts))

(defn adaptive
  "Accepts two 'integrator' functions of:

  - `f`: some integrand
  - `a` and `b`: the lower and upper endpoints of integration
  - `opts`, a dictionary of configuration options

  And returns a new function of the same signature that adaptively subdivides
  the region $a, b$ into intervals if integration fails to converge."
  ([integrator] (adaptive integrator integrator))
  ([open-integrator closed-integrator]
   (fn rec
     ([f a b] (rec f a b {}))
     ([f a b opts]
      (let [opts      (fill-defaults opts)
            integrate (fn [l r interval]
                        (if (closed? interval)
                          (closed-integrator f l r opts)
                          (open-integrator f l r opts)))]
        (loop [sum   (ua/kahan-sum)
               stack [[a b (:interval opts)]]]
          (if (empty? stack)
            (first sum)
            (let [[l r interval] (peek stack)
                  remaining      (pop stack)
                  {:keys [converged? result]} (integrate l r interval)]
              (if converged?
                (recur (ua/kahan-sum sum result) remaining)
                (let [midpoint (split-point l r (:fuzz-factor opts))]
                  (recur sum (conj remaining
                                   [l midpoint (close-r interval)]
                                   [midpoint r (close-l interval)]))))))))))))