Last active
October 16, 2020 04:28
-
-
Save sritchie/49a23138f597f86f89fb83864071fd3f to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
(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)))))))))) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
(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)])))))))))))) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment