Skip to content

v0.13.0: Clojurescript, Functional numerics, power series!

Compare
Choose a tag to compare
@sritchie sritchie released this 09 Nov 19:38

The main announcement for this release is Clojurescript Support!. Implementing
this resulted in a few upgrades along the way:

  • an overhauled Power Series implementation based on Doug McIlroy's "Power Series, Power Serious"
  • lots more generic implementations for Operator and Matrix
  • more powerful numerics, specifically definite-integral and native minimization routines
  • a generic numeric tower for Clojurescript
  • Many more tests! The test coverage was great before, and it's stayed high as
    we've added new implementations.
  • added explicit code coverage metrics via Codecov: Codecov branch

Here are more explicit details on the release.

Clojurescript Support

Full conversion of SICMUtils to Clojurescript. All functionality from v0.12.1
now works in both Clojure and Clojurescript!

Most of the conversion was straightforward. The major missing piece was a
numeric tower implementation for Clojurescript (complex numbers, ratios) that
bring it up to parity with Clojure:

  • Add the bigfraction implementation from
    fraction.js sicmutils.ratio for
    cross-platform ratio support (#99)
  • Adds CLJS complex number support through complex.js (#41)
  • js/BigInt, goog.math.Long and goog.math.Integer implementations round
    out the numeric tower (#45)

Numerical Routines

The numerical routines in SICMUtils depend heavily on Apache Commons, which of
course only exists in Java. We had to implement much of the numerics code in
native Clojure. It's fast, efficient and functional. Give it a read if you're
curious about how these algorithms work.

  • New, native minimization routines have replaced the Apache Commons implementations:

    • Univariate Minimizers

      • Port scipy's auto-bracketing + scmutils version (#104)
      • Port golden section search from scipy (#105)
      • Implement Brent's method for fn minimization in native clj (#106)
    • Multivariate

      • pure Clojure implementation of Nelder-Mead (#102)
  • Native definite-integral numerics implementation, written as a series of
    computational essays:

    • Basics:

      • Riemann Sums, all the way up through efficient, incremental, "accelerated" versions of these easy-to-understand methods:
      • Midpoint method, same development but shorter since it reuses functional abstractions. Also incremental, efficient, accelerated
      • Trapezoid Method, same idea but for closed intervals.
    • Sequence Acceleration / Extrapolation Methods

      • Polynomial interpolation: the general thing that "richardson extrapolation" is doing below. Historically cool and used to accelerate arbitrary integration sequences
      • Rational Function extrapolation: used in bulirsch-stoer integration and ODE solving.
      • "Richardson extrapolation" is a special case, where we get more efficient by assuming that the x values for the polynomial interpolation go 1, 1/2, 1/4... and that we're extrapolating to 0.
    • Higher-order Calculus:

      • Numerical derivatives: derivatives using three kinds of central difference formulas... accelerated using Richardson extrapolation, with a nice technique for guarding against underflow.
      • Simpson's Method... fit a parabola to every slice. OR, "accelerate" the trapezoid method with one step of Richarsdson extrapolation!
      • Simpson's 3/8 Method: Same idea, but accelerate a sequence that triples its slices every iteration.
      • Boole's Rule: trapezoid method plus two steps of Richardson extrapolation. (Are you starting to see the pattern??)
      • Romberg Integration: midpoint OR trapezoid, with as many steps of Richardson extrapolation as we can take!
      • Milne's Rule, MIDPOINT method, one step of extrapolation!
      • Bulirsch-Stoer integration... midpoint or trapezoid, with rational function extrapolation, as many steps as we can handle AND some custom step sizes.
    • Combinators:

      • Variable Substitutions: implemented as functional wrappers that take an integrator and return a modified integrator.
      • Improper Integrals: a template for a combinator that enables infinite endpoints on any integrator, using variable substitution on an appropriate, tunable range.
      • Adaptive Integration: a combinator that turns any of the integrators above into an "adaptive" integrator that's able to focus in on difficult regions.
    • And finally, "Numerical Quadrature", the namespace/essay that ties it all together.

  • sicmutils.numerical.compile uses SCI, the
    Small Clojure Interpreter, to generate compiled numerical code (#133)

  • Implemented ODE solving using @littleredcomputer's
    odex-js library (#135)

Reader Literals

data_readers.cljc provides 3 new data reader literals:

#sicm/ratio

Use this with a ratio literal, like #sicm/ratio 1/2, or with a string like
#sicm/ratio "1/4". If the denominator is 1 this literal will return a
js/BigInt in Clojurescript, or a Long in Clojure.

#sicm/bigint

Use with a number literal, like, #sicm/bigint 10, or a string like
#sicm/bigint "10000012" to generate a js/BigInt in Clojurescript, or a
clojure.lang.BigInt in Clojure.

#sicm/complex

Currently this only works with a string like #sicm/complex "1 + 2i". In the
future it might work with a pair of (real, complex), like:

#sicm/complex [1 2]

Power Series, Power Serious

The Power Series implementation in series.cljc received an overhaul. The
implementation now follows Doug McIlroy's beautiful paper, "Power Series, Power
Serious"
.
Doug also has a 10-line version in Haskell on his
website
.

The API now offers two types:

  • Series, which represents a generic infinite series of arbitrary values, and

  • PowerSeries, a series that represents a power series in a single
    variable; in other words, a series where the nth entry is interpreted as
    the coefficient of $x^n$:

    $$[a b c d ...] == $a + bx + cx^2 + dx^3 + ...$$

series/series? responds true to both. series/power-series? only responds
true to a PowerSeries.

To turn a PowerSeries into a Series, call it as a function with a single
argument, or pass the series and one argument to series/value to evaluate the
series using the above equation.

To turn a Series into a PowerSeries, call series/->function. None of the
functions discussed below can mix series types; you'll have to do the conversion
explicitly.

Each type supports the following generic operations:

  • *, +, -, / between series and non-series
  • g/negate, g/invert, g/sqrt, g/expt work as expected.
  • g/add between series and non-series

PowerSeries additionally supports:

  • g/exp, g/cos, g/sin, g/asin, g/tan
  • g/partial-derivative, so PowerSeries works well with D

Each of these acts as function composition for the single variable function that
the PowerSeries represents. If s is a PowerSeries that applies as (s x),
(g/exp s) returns a series that represents (g/exp (s x)).

There are many more new methods (see the namespace for full documentation):

  • starting-with renamed to series
  • power-series, analogous series but generates a PowerSeries
  • series* and power-series* let you pass an explicit sequence
  • series/take removed in favor of clojure.core/take, since both series
    objects act as sequences
  • generate takes an additional optional argument to distinguish between series
    and power series
  • Series now implements more of v/Value
  • new zero, one, identity constants
  • constant returns a constant power series
  • xpow generates a series representing a bare power of x
  • ->function turns a Series into a PowerSeries
  • value, fmap now handles both Series and PowerSeries
  • (inflate s n) expands each term $x^i$ of s to $x^{in}$
  • compose returns the functional composition of two PowerSeries
  • revert returns the functional inverse of two PowerSeries
  • integral returns a series representing the definite integral of the supplied
    PowerSeries, 0 => infinity (optionally takes an integration constant)

The namespace also provides many built-in PowerSeries instances:

  • exp-series
  • sin-series
  • cos-series
  • tan-series
  • sec-series
  • asin-series
  • acos-series
  • atan-series
  • acot-series
  • sinh-series
  • cosh-series
  • tanh-series
  • asinh-series
  • atanh-series
  • log1+x-series
  • log1-x-series
  • binomial-series

And two Series (non-power):

Matrix Generic Operations

::matrix gained implementations for exp, cos, sin, asin, tan,
acos, asin; these now return taylor series expansions of the operator, where
multiplication is composition as before.

Operator Generics

Operator gained implementations for cos, sin, asin, tan, acos,
asin; these now return taylor series expansions of the operator, where
multiplication is composition as before.

Misc

generic.cljc now includes a default implementation of:

  • expt given a valid mul
  • default sub, given a valid add and negate
  • default div, given a valid mul and invert
  • Expression and Operator both have better print-method implementations,
    so the repl experience feels more like scmutils
  • Operator now has an expn method, which acts like g/exp on an operator
    but expands each term in order n.
  • many, many more tests!