Intervals

This dynamic variable holds the default "roundoff cutoff" used by all integrators to decide when to return an estimate based on a single slice, vs attempting to converge a sequence of progressively finer estimates. When this condition is satisfied:

|b - a| / |a| + |b| <= cutoff

An integrator will estimate a single slice directly. Else, it will attempt to converge a sequence.

NOTE - we don't have an interface yet to bind this dynamic variable. bind it manually to modify the cutoff for a specific call to some integrator:

(binding [roundoff-cutoff 1e-6] (integrate f a b))

(def ^:dynamic *roundoff-cutoff* 1e-14)

1e-14

Intervals

Implementations of the various intervals used by the adaptive integral interface. By default, integration endpoints are considered open.

(def open        [::open ::open])

[:emmy.numerical.quadrature.common/open :emmy.numerical.quadrature.common/open]

(def closed      [::closed ::closed])

[:emmy.numerical.quadrature.common/closed :emmy.numerical.quadrature.common/closed]

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

[:emmy.numerical.quadrature.common/open :emmy.numerical.quadrature.common/closed]

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

[:emmy.numerical.quadrature.common/closed :emmy.numerical.quadrature.common/open]

(def infinities #{##Inf ##-Inf})

#{-Infinity Infinity}

(defn closed?

  "Returns true if the argument represents an explicit `closed` interval, false

  otherwise."

  [x] (= x closed))

#object[emmy.numerical.quadrature.common$closed_QMARK_ 0x33e29cf "

emmy.numerical.quadrature.common$closed_QMARK_@33e29cf"

]

(def open? (complement closed?))

#object[clojure.core$complement$fn__5737 0x4eab031f "

clojure.core$complement$fn__5737@4eab031f"

]

These functions modify an interval by opening or closing either of its endpoints.

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

#object[emmy.numerical.quadrature.common$close_l 0x15576fb3 "

emmy.numerical.quadrature.common$close_l@15576fb3"

]

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

#object[emmy.numerical.quadrature.common$close_r 0x36a2eccc "

emmy.numerical.quadrature.common$close_r@36a2eccc"

]

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

#object[emmy.numerical.quadrature.common$open_l 0x44a0db86 "

emmy.numerical.quadrature.common$open_l@44a0db86"

]

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

#object[emmy.numerical.quadrature.common$open_r 0x52858bf6 "

emmy.numerical.quadrature.common$open_r@52858bf6"

]

(defn flip [[l r]] [r l])

#object[emmy.numerical.quadrature.common$flip 0x7c015ce2 "

emmy.numerical.quadrature.common$flip@7c015ce2"

]

(defn interval

  "Extracts the interval (or `open` as a default) from the supplied integration

  options dict."

  [opts]

  (get opts :interval open))

#object[emmy.numerical.quadrature.common$interval 0x5d640b22 "

emmy.numerical.quadrature.common$interval@5d640b22"

]

(defn with-interval

  "Sets the specified interval to a key inside the suppled `opts` map of arbitrary

  integration options."

  [opts interval]

  (assoc opts :interval interval))

#object[emmy.numerical.quadrature.common$with_interval 0x1fda97e3 "

emmy.numerical.quadrature.common$with_interval@1fda97e3"

]

(defn update-interval

  "Accepts:

  - a dictionary of arbitrary options

  - one of the 4 interval modification functions

  and returns a dict of options with `f` applied to the contained interval (or

  `open` if no interval is set).

"

  [opts f]

  (let [k :interval]

    (assoc opts k (f (interval opts)))))

#object[emmy.numerical.quadrature.common$update_interval 0x617e65d5 "

emmy.numerical.quadrature.common$update_interval@617e65d5"

]

Common Integration Interface

The following two functions define a shared interface that integration namespaces can use to create an "integrator" from:

a fn that can estimate the area of a single integration slice, and
a fn that can return a sequence of progressively finer estimates.

The first function is called in the case that the integration range $(a, b)$ (open or closed) is too fine for subdivision. The second function takes over in all other (most!) cases.

(defn- narrow-slice?

  "Returns true if the range $[a, b]$ is strip narrow enough to pass the following

  test:

  |b - a| / |a| + |b| <= `cutoff`

  False otherwise. This inequality measures how close the two floating point

  values are, scaled by the sum of their magnitudes."

  [^double a ^double b cutoff]

  (let [sum (+ (Math/abs a)

               (Math/abs b))]

    (or (<= sum cutoff)

        (<= (Math/abs (- b a))

            (* cutoff sum)))))

#object[emmy.numerical.quadrature.common$narrow_slice_QMARK_ 0x176ed354 "

emmy.numerical.quadrature.common$narrow_slice_QMARK_@176ed354"

]

(defn make-integrator-fn

  "Generates an `integrator` function from two functions with the following

  signatures and descriptions:

  - `(area-fn f a b)` estimates the integral of `f` over the interval `(a, b)`

  with no subdivision, nothing clever at all.

  - `(seq-fn f a b opts)` returns a sequence of successively refined estimates

  of the integral of `f` over `(a, b)`. `opts` can contain kv pairs that

  configure the behavior of the sequence function (a sequence of the number of

  integration slices to use, for example.)

  The returned function has the signature:

  `(f a b opts)`

  All `opts` are passed on to `seq-fn`, _and_ to `us/seq-limit` internally,

  where the options configure the checks on sequence convergence."

  [area-fn seq-fn]

  (fn call

    ([f a b] (call f a b {}))

    ([f a b {:keys [roundoff-cutoff]

             :or {roundoff-cutoff *roundoff-cutoff*}

             :as opts}]

     (if (narrow-slice? a b roundoff-cutoff)

       (do (log/info "Integrating narrow slice: " a b)

           {:converged? true

            :terms-checked 1

            :result (area-fn f a b)})

       (-> (seq-fn f a b opts)

           (us/seq-limit opts))))))

#object[emmy.numerical.quadrature.common$make_integrator_fn 0x6a2ef780 "

emmy.numerical.quadrature.common$make_integrator_fn@6a2ef780"

]

(defn- name-with-attributes

  "Taken from `clojure.tools.macro/name-with-attributes`.

  Handles optional docstrings and attribute maps for a name to be defined in a

  list of macro arguments. If the first macro argument is a string, it is added

  as a docstring to name and removed from the macro argument list. If afterwards

  the first macro argument is a map, its entries are added to the name's

  metadata map and the map is removed from the macro argument list. The return

  value is a vector containing the name with its extended metadata map and the

  list of unprocessed macro arguments."

  ([name body] (name-with-attributes name body {}))

  ([name body meta]

   (let [[docstring body] (if (string? (first body))

                            [(first body) (next body)]

                            [nil body])

         [attr body]      (if (map? (first body))

                            [(first body) (next body)]

                            [{} body])

         attr             (merge meta attr)

         attr             (if docstring

                            (assoc attr :doc docstring)

                            attr)

         attr             (if (meta name)

                            (conj (meta name) attr)

                            attr)]

     [(with-meta name attr) body])))

#object[emmy.numerical.quadrature.common$name_with_attributes 0x210010a9 "

emmy.numerical.quadrature.common$name_with_attributes@210010a9"

]

(u/sci-macro defintegrator

  "Helper macro for defining integrators."

  [sym & body]

  (let [meta                     {:arglists (list 'quote '([f a b] [f a b opts]))}

        [sym body]               (name-with-attributes sym body meta)

        {:keys [area-fn seq-fn]} (apply hash-map body)]

    (assert seq-fn (str "defintegrator " sym ": seq-fn cannot be nil"))

    (assert area-fn (str "defintegrator " sym ": area-fn cannot be nil"))

    `(def ~sym

       (make-integrator-fn ~area-fn ~seq-fn))))

#'emmy.numerical.quadrature.common/defintegrator