diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 81712c3..3018aef 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.11.1","generation_timestamp":"2024-11-05T21:09:46","documenter_version":"1.7.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.11.1","generation_timestamp":"2024-11-18T07:01:45","documenter_version":"1.8.0"}} \ No newline at end of file diff --git a/dev/api/index.html b/dev/api/index.html index a1bd00e..dcf1c01 100644 --- a/dev/api/index.html +++ b/dev/api/index.html @@ -1,6 +1,6 @@ -API · TaylorDiff.jl

API

API for TaylorDiff.

TaylorDiff.TaylorArray — Type
TaylorArray{T, N, A, P}

Representation of Taylor polynomials in array mode.

Fields

  • value::A: zeroth order coefficient
  • partials::NTuple{P, A}: i-th element of this stores the i-th derivative
source
TaylorDiff.TaylorScalar — Type
TaylorScalar{T, P}

Representation of Taylor polynomials.

Fields

  • value::T: zeroth order coefficient
  • partials::NTuple{P, T}: i-th element of this stores the i-th derivative
source
TaylorDiff.TaylorScalar — Method
TaylorScalar{P}(value::T, seed::T)

Convenience function: construct a Taylor polynomial with zeroth and first order coefficient, acting as a seed.

source
TaylorDiff.TaylorScalar — Method
TaylorScalar{P}(value::T) where {T, P}

Convenience function: construct a Taylor polynomial with zeroth order coefficient.

source
TaylorDiff.can_taylorize — Method
TaylorDiff.can_taylorize(V::Type)

Determines whether the type V is allowed as the scalar type in a Dual. By default, only <:Real types are allowed.

source
TaylorDiff.derivative — Function
derivative(f, x, ::Val{P})
+API · TaylorDiff.jl
API
API for TaylorDiff.
TaylorDiff.TaylorArray — Type
TaylorArray{T, N, A, P}
Representation of Taylor polynomials in array mode.
Fields
  • value::A: zeroth order coefficient
  • partials::NTuple{P, A}: i-th element of this stores the i-th derivative
source
TaylorDiff.TaylorScalar — Type
TaylorScalar{T, P}
Representation of Taylor polynomials.
Fields
  • value::T: zeroth order coefficient
  • partials::NTuple{P, T}: i-th element of this stores the i-th derivative
source
TaylorDiff.TaylorScalar — Method
TaylorScalar{P}(value::T, seed::T)
Convenience function: construct a Taylor polynomial with zeroth and first order coefficient, acting as a seed.
source
TaylorDiff.TaylorScalar — Method
TaylorScalar{P}(value::T) where {T, P}
Convenience function: construct a Taylor polynomial with zeroth order coefficient.
source
TaylorDiff.can_taylorize — Method
TaylorDiff.can_taylorize(V::Type)
Determines whether the type V is allowed as the scalar type in a Dual. By default, only <:Real types are allowed.
source
TaylorDiff.derivative — Function
derivative(f, x, ::Val{P})
 derivative(f, x, l, ::Val{P})
-derivative(f!, y, x, l, ::Val{P})
Computes P-th directional derivative of f w.r.t. vector x in direction l. If x is a Number, the direction l can be omitted.
source
TaylorDiff.derivative! — Function
derivative!(result, f, x, l, ::Val{P})
-derivative!(result, f!, y, x, l, ::Val{P})
In-place derivative calculation APIs. result is expected to be pre-allocated and have the same shape as y.
source
TaylorDiff.derivatives — Function
derivatives(f, x, l, ::Val{P})
-derivatives(f!, y, x, l, ::Val{P})
Computes all derivatives of f at x up to order P.
source
TaylorDiff.get_term_raiser — Method
Pick a strategy for raising the derivative of a function. If the derivative is like 1 over something, raise with the division rule; otherwise, raise with the multiplication rule.
source
TaylorDiff.@immutable — Macro
immutable(def)
Transform a function definition to a @generated function.
  1. Allocations are removed by replacing the output with scalar variables;
  2. Loops are unrolled;
  3. Indices are modified to use 1-based indexing;
source

+derivative(f!, y, x, l, ::Val{P})

Computes P-th directional derivative of f w.r.t. vector x in direction l. If x is a Number, the direction l can be omitted.

source
TaylorDiff.derivative! — Function
derivative!(result, f, x, l, ::Val{P})
+derivative!(result, f!, y, x, l, ::Val{P})

In-place derivative calculation APIs. result is expected to be pre-allocated and have the same shape as y.

source
TaylorDiff.derivatives — Function
derivatives(f, x, l, ::Val{P})
+derivatives(f!, y, x, l, ::Val{P})

Computes all derivatives of f at x up to order P.

source
TaylorDiff.get_term_raiser — Method

Pick a strategy for raising the derivative of a function. If the derivative is like 1 over something, raise with the division rule; otherwise, raise with the multiplication rule.

source
TaylorDiff.@immutable — Macro
immutable(def)

Transform a function definition to a @generated function.

  1. Allocations are removed by replacing the output with scalar variables;
  2. Loops are unrolled;
  3. Indices are modified to use 1-based indexing;
source
diff --git a/dev/assets/documenter.js b/dev/assets/documenter.js index 82252a1..7d68cd8 100644 --- a/dev/assets/documenter.js +++ b/dev/assets/documenter.js @@ -612,176 +612,194 @@ function worker_function(documenterSearchIndex, documenterBaseURL, filters) { }; } -// `worker = Threads.@spawn worker_function(documenterSearchIndex)`, but in JavaScript! -const filters = [ - ...new Set(documenterSearchIndex["docs"].map((x) => x.category)), -]; -const worker_str = - "(" + - worker_function.toString() + - ")(" + - JSON.stringify(documenterSearchIndex["docs"]) + - "," + - JSON.stringify(documenterBaseURL) + - "," + - JSON.stringify(filters) + - ")"; -const worker_blob = new Blob([worker_str], { type: "text/javascript" }); -const worker = new Worker(URL.createObjectURL(worker_blob)); - /////// SEARCH MAIN /////// -// Whether the worker is currently handling a search. This is a boolean -// as the worker only ever handles 1 or 0 searches at a time. -var worker_is_running = false; - -// The last search text that was sent to the worker. This is used to determine -// if the worker should be launched again when it reports back results. -var last_search_text = ""; - -// The results of the last search. This, in combination with the state of the filters -// in the DOM, is used compute the results to display on calls to update_search. -var unfiltered_results = []; - -// Which filter is currently selected -var selected_filter = ""; - -$(document).on("input", ".documenter-search-input", function (event) { - if (!worker_is_running) { - launch_search(); - } -}); - -function launch_search() { - worker_is_running = true; - last_search_text = $(".documenter-search-input").val(); - worker.postMessage(last_search_text); -} - -worker.onmessage = function (e) { - if (last_search_text !== $(".documenter-search-input").val()) { - launch_search(); - } else { - worker_is_running = false; - } - - unfiltered_results = e.data; - update_search(); -}; +function runSearchMainCode() { + // `worker = Threads.@spawn worker_function(documenterSearchIndex)`, but in JavaScript! + const filters = [ + ...new Set(documenterSearchIndex["docs"].map((x) => x.category)), + ]; + const worker_str = + "(" + + worker_function.toString() + + ")(" + + JSON.stringify(documenterSearchIndex["docs"]) + + "," + + JSON.stringify(documenterBaseURL) + + "," + + JSON.stringify(filters) + + ")"; + const worker_blob = new Blob([worker_str], { type: "text/javascript" }); + const worker = new Worker(URL.createObjectURL(worker_blob)); + + // Whether the worker is currently handling a search. This is a boolean + // as the worker only ever handles 1 or 0 searches at a time. + var worker_is_running = false; + + // The last search text that was sent to the worker. This is used to determine + // if the worker should be launched again when it reports back results. + var last_search_text = ""; + + // The results of the last search. This, in combination with the state of the filters + // in the DOM, is used compute the results to display on calls to update_search. + var unfiltered_results = []; + + // Which filter is currently selected + var selected_filter = ""; + + $(document).on("input", ".documenter-search-input", function (event) { + if (!worker_is_running) { + launch_search(); + } + }); -$(document).on("click", ".search-filter", function () { - if ($(this).hasClass("search-filter-selected")) { - selected_filter = ""; - } else { - selected_filter = $(this).text().toLowerCase(); + function launch_search() { + worker_is_running = true; + last_search_text = $(".documenter-search-input").val(); + worker.postMessage(last_search_text); } - // This updates search results and toggles classes for UI: - update_search(); -}); + worker.onmessage = function (e) { + if (last_search_text !== $(".documenter-search-input").val()) { + launch_search(); + } else { + worker_is_running = false; + } -/** - * Make/Update the search component - */ -function update_search() { - let querystring = $(".documenter-search-input").val(); + unfiltered_results = e.data; + update_search(); + }; - if (querystring.trim()) { - if (selected_filter == "") { - results = unfiltered_results; + $(document).on("click", ".search-filter", function () { + if ($(this).hasClass("search-filter-selected")) { + selected_filter = ""; } else { - results = unfiltered_results.filter((result) => { - return selected_filter == result.category.toLowerCase(); - }); + selected_filter = $(this).text().toLowerCase(); } - let search_result_container = ``; - let modal_filters = make_modal_body_filters(); - let search_divider = `
`; + // This updates search results and toggles classes for UI: + update_search(); + }); - if (results.length) { - let links = []; - let count = 0; - let search_results = ""; - - for (var i = 0, n = results.length; i < n && count < 200; ++i) { - let result = results[i]; - if (result.location && !links.includes(result.location)) { - search_results += result.div; - count++; - links.push(result.location); - } - } + /** + * Make/Update the search component + */ + function update_search() { + let querystring = $(".documenter-search-input").val(); - if (count == 1) { - count_str = "1 result"; - } else if (count == 200) { - count_str = "200+ results"; + if (querystring.trim()) { + if (selected_filter == "") { + results = unfiltered_results; } else { - count_str = count + " results"; + results = unfiltered_results.filter((result) => { + return selected_filter == result.category.toLowerCase(); + }); } - let result_count = `
${count_str}
`; - search_result_container = ` + let search_result_container = ``; + let modal_filters = make_modal_body_filters(); + let search_divider = `
`; + + if (results.length) { + let links = []; + let count = 0; + let search_results = ""; + + for (var i = 0, n = results.length; i < n && count < 200; ++i) { + let result = results[i]; + if (result.location && !links.includes(result.location)) { + search_results += result.div; + count++; + links.push(result.location); + } + } + + if (count == 1) { + count_str = "1 result"; + } else if (count == 200) { + count_str = "200+ results"; + } else { + count_str = count + " results"; + } + let result_count = `
${count_str}
`; + + search_result_container = ` +
+ ${modal_filters} + ${search_divider} + ${result_count} +
+ ${search_results} +
+
+ `; + } else { + search_result_container = `
${modal_filters} ${search_divider} - ${result_count} -
- ${search_results} -
-
+
0 result(s)
+ +
No result found!
`; - } else { - search_result_container = ` -
- ${modal_filters} - ${search_divider} -
0 result(s)
-
-
No result found!
- `; - } + } - if ($(".search-modal-card-body").hasClass("is-justify-content-center")) { - $(".search-modal-card-body").removeClass("is-justify-content-center"); - } + if ($(".search-modal-card-body").hasClass("is-justify-content-center")) { + $(".search-modal-card-body").removeClass("is-justify-content-center"); + } - $(".search-modal-card-body").html(search_result_container); - } else { - if (!$(".search-modal-card-body").hasClass("is-justify-content-center")) { - $(".search-modal-card-body").addClass("is-justify-content-center"); + $(".search-modal-card-body").html(search_result_container); + } else { + if (!$(".search-modal-card-body").hasClass("is-justify-content-center")) { + $(".search-modal-card-body").addClass("is-justify-content-center"); + } + + $(".search-modal-card-body").html(` +
Type something to get started!
+ `); } + } - $(".search-modal-card-body").html(` -
Type something to get started!
- `); + /** + * Make the modal filter html + * + * @returns string + */ + function make_modal_body_filters() { + let str = filters + .map((val) => { + if (selected_filter == val.toLowerCase()) { + return `${val}`; + } else { + return `${val}`; + } + }) + .join(""); + + return ` +
+ Filters: + ${str} +
`; } } -/** - * Make the modal filter html - * - * @returns string - */ -function make_modal_body_filters() { - let str = filters - .map((val) => { - if (selected_filter == val.toLowerCase()) { - return `${val}`; - } else { - return `${val}`; - } - }) - .join(""); - - return ` -
- Filters: - ${str} -
`; +function waitUntilSearchIndexAvailable() { + // It is possible that the documenter.js script runs before the page + // has finished loading and documenterSearchIndex gets defined. + // So we need to wait until the search index actually loads before setting + // up all the search-related stuff. + if (typeof documenterSearchIndex !== "undefined") { + runSearchMainCode(); + } else { + console.warn("Search Index not available, waiting"); + setTimeout(waitUntilSearchIndexAvailable, 1000); + } } +// The actual entry point to the search code +waitUntilSearchIndexAvailable(); + }) //////////////////////////////////////////////////////////////////////////////// require(['jquery'], function($) { diff --git a/dev/index.html b/dev/index.html index d41d17f..b635332 100644 --- a/dev/index.html +++ b/dev/index.html @@ -11,4 +11,4 @@ publisher = {GitHub}, journal = {GitHub repository}, howpublished = {\url{https://github.com/JuliaDiff/TaylorDiff.jl}} -} +} diff --git a/dev/theory/index.html b/dev/theory/index.html index 15d44a1..2ffd3ec 100644 --- a/dev/theory/index.html +++ b/dev/theory/index.html @@ -1,3 +1,3 @@ Theory · TaylorDiff.jl

Theory

TaylorDiff.jl is an operator-overloading based forward-mode automatic differentiation (AD) package. "Forward-mode" implies that the basic capability of this package is that, for function $f:\mathbb R^n\to\mathbb R^m$, place to evaluate derivative $x\in\mathbb R^n$ and direction $l\in\mathbb R^n$, we compute $ f(x),\partial f(x)\times v,\partial^2f(x)\times v\times v,\cdots,\partial^pf(x)\times v\times\cdots\times v $ i.e., the function value and the directional derivative up to order $p$. This notation might be unfamiliar to Julia users that had experience with other AD packages, but $\partial f(x)$ is simply the jacobian $J$, and $\partial f(x)\times v$ is simply the Jacobian-vector product (jvp). In other words, this is a simple generalization of Jacobian-vector product to Hessian-vector-vector product, and to even higher orders.

The main advantage of doing this instead of doing $p$ first-order Jacobian-vector products is that nesting first-order AD results in expential scaling w.r.t $p$, while this method, also known as Taylor mode, should be (almost) linear scaling w.r.t $p$. We will see the reason of this claim later.

In order to achieve this, assuming that $f$ is a nested function $f_k\circ\cdots\circ f_2\circ f_1$, where each $f_i$ is a basic and simple function, or called "primitives". We need to figure out how to propagate the derivatives through each step. In first order AD, this is achieved by the "dual" pair $x_0+x_1\varepsilon$, where $\varepsilon^2=0$, and for each primitive we make a method overload $ f(x0+x1\varepsilon)=f(x0)+\partial f(x0) x1\varepsilon $ Similarly in higher-order AD, we need for each primitive a method overload for a truncated Taylor polynomial up to order $p$, and in this polynomial we will use $t$ instead of $\varepsilon$ to denote the sensitivity. "Truncated" means $t^{p+1}=0$, similar as what we defined for dual numbers. So $ f(x0+x1t+x2t^2+\cdots+x_pt^p)=? $ What is the math expression that we should put into the question mark? That specific expression is called the "pushforward rule", and we will talk about how to derive the pushforward rule below.

Arithmetic of polynomials

Before deriving pushforward rules, let's first introduce several basic properties of polynomials.

If $x(t)$ and $y(t)$ are both truncated Taylor polynomials, i.e. $ \begin{aligned} x&=x0+x1t+\cdots+xpt^p\ -y&=y0+y1t+\cdots+ypt^p \end{aligned} $ Then it's obvious that the polynomial addition and subtraction should be $ (x\pm y)k=xk\pm yk $ And with some derivation we can also get the polynomial multiplication rule $ (x\times y)k=\sum{i=0}^kxiy{k-i} $ The polynomial division rule is less obvious, but if $x/y=z$, then equivalently $x=yz$, i.e. $ \left(\sum{i=0}^pyit^i\right)\left(\sum{i=0}^pzit^i\right)=\sum{i=0}^pxit^i $ if we relate the coefficient of $t^k$ on both sides we get $ \sum{i=0}^k ziy{k-i}=xk $ so, equivalently, $ zk=\frac1{y0}\left(xk-\sum{i=0}^{k-1}ziy{k-1}\right) $ This is a recurrence relation, which means that we can first get z0=x0/y0$, and then get $z_1$ using $z_0$, and then get $z_2$ using $z_0,z_1$ etc.

Pushforward rule for elementary functions

Let's now consider how to derive the pushforward rule for elementary functions. We will use $\exp$ and $\log$ as two examples.

If $x(t)$ is a polynomial and we want to get $e(t)=\exp(x(t))$, we can actually get that by formulating an ordinary differential equation: $ e'(t)=\exp(x(t))x'(t);\quad e0=\exp(x0) $ If we expand both $e$ and $x$ in the equation, we will get $ \sum{i=1}^pieit^{i-1}=\left(\sum{i=0}^{p-1} eit^i\right)\left(\sum{i=1}^pixit^{i-1}\right) $ relating the coefficient of $t^{k-1}$ on both sides, we get $ kek=\sum{i=0}^{k-1}ei\times (k-i)x{k-i} $ This is, again, a recurrence relation, so we can get $e_1,\cdots,e_p$ step-by-step.

If $x(t)$ is a polynomial and we want to get $l(t)=\log(x(t))$, we can actually get that by formulating an ordinary differential equation: $ l'(t)=\frac1xx'(t);\quad l0=\log(x0) $ If we expand both $l$ and $x$ in the equation, the RHS is simply polynomial divisions, and we get $ lk=\frac1{x0}\left(xk-\frac1k\sum{i=1}^{k-1}ilix{k-j}\right) $

Now notice the difference between the rule for $\exp$ and $\log$: the derivative of exponentiation is itself, so we can obtain from recurrence relation; the derivative of logarithm is $1/x$, an algebraic expression in $x$, so it can be directly computed. Similarly, we have $(\tan x)'=1+\tan^2x$ but $(\arctan x)'=(1+x^2)^{-1}$. We summarize (omitting proof) that

  • Every $\exp$-like function (like $\sin$, $\cos$, $\tan$, $\sinh$, ...)'s derivative is somehow recursive
  • Every $\log$-like function (like $\arcsin$, $\arccos$, $\arctan$, $\operatorname{arcsinh}$, ...)'s derivative is algebraic

So all of the elementary functions have an easy pushforward rule that can be computed within $O(p^2)$ time. Note that this is an elegant and straightforward corollary from the definition of "elementary function" in differential algebra.

Generic pushforward rule

For a generic $f(x)$, if we don't bother deriving the specific recurrence rule for it, we can still automatically generate pushforward rule in the following manner. Let's denote the derivative of $f$ w.r.t $x$ to be $d(x)$, then for $f(t)=f(x(t))$ we have $ f'(t)=d(x(t))x'(t);\quad f(0)=f(x_0) $ when we expand $f$ and $x$ up to order $p$ into this equation, we notice that only order $p-1$ is needed for $d(x(t))$. In other words, we turn a problem of finding $p$-th order pushforward for $f$, to a problem of finding $p-1$-th order pushforward for $d$, and we can recurse down to the first order. The first-order derivative expressions are captured from ChainRules.jl, which made this process fully automatic.

This strategy is in principle equivalent to nesting first-order differentiation, which could potentially leads to exponential scaling; however, in practice there is a huge difference. This generation of pushforward rule happens at compile time, which gives the compiler a chance to check redundant expressions and optimize it down to quadratic time. Compiler has stack limits but this should work for at least up to order 100.

In the current implementation of TaylorDiff.jl, all $\log$-like functions' pushforward rules are generated by this strategy, since their derivatives are simple algebraic expressions; some $\exp$-like functions, like sinh, is also generated; the most-often-used several $\exp$-like functions are hand-written with hand-derived recurrence relations.

If you find that the code generated by this strategy is slow, please file an issue and we will look into it.

+y&=y0+y1t+\cdots+ypt^p \end{aligned} $ Then it's obvious that the polynomial addition and subtraction should be $ (x\pm y)k=xk\pm yk $ And with some derivation we can also get the polynomial multiplication rule $ (x\times y)k=\sum{i=0}^kxiy{k-i} $ The polynomial division rule is less obvious, but if $x/y=z$, then equivalently $x=yz$, i.e. $ \left(\sum{i=0}^pyit^i\right)\left(\sum{i=0}^pzit^i\right)=\sum{i=0}^pxit^i $ if we relate the coefficient of $t^k$ on both sides we get $ \sum{i=0}^k ziy{k-i}=xk $ so, equivalently, $ zk=\frac1{y0}\left(xk-\sum{i=0}^{k-1}ziy{k-1}\right) $ This is a recurrence relation, which means that we can first get z0=x0/y0$, and then get $z_1$ using $z_0$, and then get $z_2$ using $z_0,z_1$ etc.

Pushforward rule for elementary functions

Let's now consider how to derive the pushforward rule for elementary functions. We will use $\exp$ and $\log$ as two examples.

If $x(t)$ is a polynomial and we want to get $e(t)=\exp(x(t))$, we can actually get that by formulating an ordinary differential equation: $ e'(t)=\exp(x(t))x'(t);\quad e0=\exp(x0) $ If we expand both $e$ and $x$ in the equation, we will get $ \sum{i=1}^pieit^{i-1}=\left(\sum{i=0}^{p-1} eit^i\right)\left(\sum{i=1}^pixit^{i-1}\right) $ relating the coefficient of $t^{k-1}$ on both sides, we get $ kek=\sum{i=0}^{k-1}ei\times (k-i)x{k-i} $ This is, again, a recurrence relation, so we can get $e_1,\cdots,e_p$ step-by-step.

If $x(t)$ is a polynomial and we want to get $l(t)=\log(x(t))$, we can actually get that by formulating an ordinary differential equation: $ l'(t)=\frac1xx'(t);\quad l0=\log(x0) $ If we expand both $l$ and $x$ in the equation, the RHS is simply polynomial divisions, and we get $ lk=\frac1{x0}\left(xk-\frac1k\sum{i=1}^{k-1}ilix{k-j}\right) $

Now notice the difference between the rule for $\exp$ and $\log$: the derivative of exponentiation is itself, so we can obtain from recurrence relation; the derivative of logarithm is $1/x$, an algebraic expression in $x$, so it can be directly computed. Similarly, we have $(\tan x)'=1+\tan^2x$ but $(\arctan x)'=(1+x^2)^{-1}$. We summarize (omitting proof) that

So all of the elementary functions have an easy pushforward rule that can be computed within $O(p^2)$ time. Note that this is an elegant and straightforward corollary from the definition of "elementary function" in differential algebra.

Generic pushforward rule

For a generic $f(x)$, if we don't bother deriving the specific recurrence rule for it, we can still automatically generate pushforward rule in the following manner. Let's denote the derivative of $f$ w.r.t $x$ to be $d(x)$, then for $f(t)=f(x(t))$ we have $ f'(t)=d(x(t))x'(t);\quad f(0)=f(x_0) $ when we expand $f$ and $x$ up to order $p$ into this equation, we notice that only order $p-1$ is needed for $d(x(t))$. In other words, we turn a problem of finding $p$-th order pushforward for $f$, to a problem of finding $p-1$-th order pushforward for $d$, and we can recurse down to the first order. The first-order derivative expressions are captured from ChainRules.jl, which made this process fully automatic.

This strategy is in principle equivalent to nesting first-order differentiation, which could potentially leads to exponential scaling; however, in practice there is a huge difference. This generation of pushforward rule happens at compile time, which gives the compiler a chance to check redundant expressions and optimize it down to quadratic time. Compiler has stack limits but this should work for at least up to order 100.

In the current implementation of TaylorDiff.jl, all $\log$-like functions' pushforward rules are generated by this strategy, since their derivatives are simple algebraic expressions; some $\exp$-like functions, like sinh, is also generated; the most-often-used several $\exp$-like functions are hand-written with hand-derived recurrence relations.

If you find that the code generated by this strategy is slow, please file an issue and we will look into it.