Skip to content

Commit

Permalink
Merge pull request #141 from bmad-sim/newTPS
Browse files Browse the repository at this point in the history
New tps
  • Loading branch information
mattsignorelli authored Jan 16, 2025
2 parents cecf153 + 7089e38 commit 222291b
Show file tree
Hide file tree
Showing 22 changed files with 13,592 additions and 6,828 deletions.
46 changes: 19 additions & 27 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -7,15 +7,21 @@

`GTPSA.jl` is a full-featured Julia interface to the [Generalised Truncated Power Series Algebra (GTPSA) library](https://madx.web.cern.ch/releases/madng/html/mad_mod_diffalg.html), which computes Taylor expansions, or Truncated Power Series (TPSs), of real and complex multivariable functions to high orders.

Truncated Power Series Algebra (TPSA) performs forward-mode automatic differentation (AD) similar to the typical dual-number implementation, as in `ForwardDiff.jl`. However, instead of nesting derivatives for higher orders, TPSA naturally extends to arbitary orders by directly using the power series expansions. Furthermore, because TPSA is designed for high order AD, extra focus is given to the data structure storing all partial derivatives, as well as indexing/propagating each of the nonzero partial derivatives in an efficient way. With these features, `GTPSA.jl` is significantly faster than `ForwardDiff.jl` for 2nd-order calculations and above, and has similar performance to 1st-order. See the [`benchmark/track.jl`](https://github.com/bmad-sim/GTPSA.jl/blob/main/benchmark/track.jl) example for a speed comparison in calculating the partial derivatives for a system with 58 inputs and 6 outputs. **In this example, GTPSA was x3.5 faster than ForwardDiff to 2nd order, and x19.8 faster to 3rd order.**
Truncated Power Series Algebra (TPSA) performs forward-mode automatic differentation (AD) similar to the typical dual-number implementation, as in `ForwardDiff.jl`. However, instead of nesting derivatives for higher orders, TPSA naturally extends to arbitary orders by directly using the power series expansions. Furthermore, because TPSA is designed for high order AD, the storage, indexing, and propagation of all nonzero partial derivatives has been highly optimized.

GTPSA provides several advantages over current Julia AD packages:
`GTPSA.jl` provides several advantages over current Julia AD packages:

1. **Speed**: `GTPSA.jl` is significantly faster than `ForwardDiff.jl` for 2nd-order calculations and above, and has very similar performance at 1st-order
2. **Easy Monomial Indexing**: Beyond 2nd order, accessing/manipulating the partial derivatives in an organized way can be a significant challenge when using other AD packages. In GTPSA, three simple indexing schemes for getting/setting monomial coefficients in a truncated power series is provided, as well as a `cycle!` function for cycling through all nonzero monomials
3. **Custom Orders in Individual Variables**: Other packages use a single maximum order for all variables. With GTPSA, the maximum order can be set differently for individual variables, as well as for a separate part of the monomial. For example, computing the Taylor expansion of $f(x_1,x_2)$ to 2nd order in $x_1$ and 6th order in $x_2$ is possible
4. **Complex Numbers**: GTPSA natively supports complex numbers and allows for mixing of complex and real truncated power series
5. **Distinction Between State Variables and Parameters**: Distinguishing between dependent variables and parameters in the solution of a differential equation expressed as a power series in the dependent variables/parameters can be advantageous in analysis
1. **Speed:** `GTPSA.jl` is significantly faster than `ForwardDiff.jl` for 2nd-order calculations and above, and has very similar performance at 1st-order. See [our example](https://github.com/bmad-sim/GTPSA.jl/blob/main/benchmark/track.jl) where **GTPSA was x3.5 faster than ForwardDiff to 2nd order, and x19.8 faster to 3rd order.**

2. **Easy Monomial Indexing:** Beyond 2nd order, accessing/manipulating the partial derivatives in an organized way can be a significant challenge. GTPSA provides three simple indexing schemes for getting/setting monomial coefficients in a truncated power series, as well as a `cycle!` function for cycling through all nonzero monomials.

3. **Easy TPS slicing:** Suppose you would like to get all terms in a multivariable truncated power series proportional to a specific variable; `GTPSA.jl` provides a simple syntax similar to array slicing to extract parts of a TPS.

4. **Complex Numbers:** `GTPSA.jl` natively supports complex numbers and allows for mixing of complex and real truncated power series

5. **Custom Orders in Individual Variables:** Other AD packages use a single maximum order for all variables. With GTPSA, the maximum order can be set differently for individual variables, as well as for a separate part of the monomial. For example, computing the Taylor expansion of $f(x_1,x_2)$ to 2nd order in $x_1$ and 6th order in $x_2$ is possible.

6. **Distinction Between State Variables and Parameters:** Distinguishing between dependent variables and parameters in the solution of a differential equation expressed as a power series in the dependent variables/parameters can be advantageous in analysis

## Setup
To use `GTPSA.jl`, in the Julia REPL run
Expand All @@ -25,28 +31,16 @@ import Pkg; Pkg.add("GTPSA")
```

## Basic Usage
First, a `Descriptor` must be created specifying the number of variables, number of parameters, and truncation order(s) for the variables/parameters in the TPSA. A `TPS` can then be created based on the `Descriptor`. TPSs can be manipulated using all of the arithmetic operators (`+`,`-`,`*`,`/`,`^`) and math functions (e.g. `abs`, `sqrt`, `sin`, `exp`, `log`, `tanh`, etc.).

TPSs can be viewed as structures containing the coefficients for all of the monomials of a multivariable Taylor expansion up to the orders specified in the `Descriptor`. As an example, to compute the truncated power series of a function $f(x_1, x_2) = \cos{(x_1)}+i\sin{(x_2)}$ to 6th order in $x_1$ and $x_2$:
```julia
using GTPSA
julia> using GTPSA

# Descriptor for TPSA with 2 variables to 6th order
d = Descriptor(2, 6)
julia> const D = Descriptor(2, 6); # 2 variables to 6th order

# Get the TPSs corresponding to each variable based on the Descriptor
x = vars(d)
# x[1] corresponds to the first variable and x[2] corresponds to the second variable
julia> Δx = @vars(D) # Get truncated power series (TPSs) corresponding to the variables

# Manipulate the TPSs as you would any other mathematical variable in Julia
f = cos(x[1]) + im*sin(x[2])
# f is a new ComplexTPS64
```

Note that scalars do not need to be defined as TPSs when writing expressions. Running `print(f)` gives the output

```
ComplexTPS64:
julia> f = cos(Δx[1]) + im*sin(Δx[2]) # Manipulate TPSs as you would any other number
ComplexTPS64{GTPSA.Dynamic}:
Descriptor(NV=2, MO=6)
Real Imag Order Exponent
1.0000000000000000e+00 0.0000000000000000e+00 0 0 0
0.0000000000000000e+00 1.0000000000000000e+00 1 0 1
Expand All @@ -59,8 +53,6 @@ ComplexTPS64:

The GTPSA library currently only supports truncated power series representing `Float64` and `ComplexF64` number types.

For more details, including using TPSs with differing truncation orders for each variable, see the [GTPSA documentation](https://bmad-sim.github.io/GTPSA.jl/).

## Acknowledgements
Much thanks must be given to Laurent Deniau, the creator of the C GTPSA library, for his time and great patience in explaining his code.

Expand Down
14 changes: 7 additions & 7 deletions benchmark/track.jl
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,7 @@ function track_cav(z0)
z[3] = z0[3]
z[4] = z0[4]
z[5] = z0[5]
z[6] = z0[6]+0.0001*z0[5]
z[6] = @FastGTPSA z0[6]+0.0001*z0[5]
return z
end

Expand Down Expand Up @@ -130,24 +130,24 @@ end

function benchmark_GTPSA1()
d = Descriptor(6,1,52,1)
z = vars()
k = params()
z = @vars(d)
k = @params(d)
map = track_ring(z, 0.36+k[1], 1.2+k[2], k[3:end])
return map
end

function benchmark_GTPSA2()
d = Descriptor(6,2,52,2)
z = vars()
k = params()
z = @vars(d)
k = @params(d)
map = track_ring(z, 0.36+k[1], 1.2+k[2], k[3:end])
return map
end

function benchmark_GTPSA3()
d = Descriptor(6,3,52,3)
z = vars()
k = params()
z = @vars(d)
k = @params(d)
map = track_ring(z, 0.36+k[1], 1.2+k[2], k[3:end])
return map
end
Expand Down
2 changes: 2 additions & 0 deletions docs/make.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
using Pkg
using Documenter, GTPSA

cp(joinpath(@__DIR__, "..", "README.md"), joinpath(@__DIR__, "src", "index.md"); force=true)

makedocs(
sitename="GTPSA.jl",
authors = "Matt Signorelli",
Expand Down
62 changes: 0 additions & 62 deletions docs/src/index.md

This file was deleted.

2 changes: 1 addition & 1 deletion docs/src/man/n_all.md
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ unsafe_convert, eps, floatmin, floatmax, signbit
`zeros` and `ones` are overloaded from Base so that allocated `TPS`s are placed in each element. Because of the mutability of `TPS`, if we didn't explicity overload these functions every element would correspond to the exact same heap-allocated TPS.

`GTPSA.jl` overloads (and exports) the following functions from the corresponding packages:
**`LinearAlgebra`**: `norm`, `mul!`
**`LinearAlgebra`**: `norm`
**`SpecialFunctions`**: `erf`, `erfc`, `erfcx`, `erfi`

`GTPSA.jl` also provides the following math functions NOT included in Base or any of the above packages (and not already documented in [TPS Methods](@ref tpsmethods)):
Expand Down
12 changes: 8 additions & 4 deletions src/GTPSA.jl
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,7 @@ import Base: +,
signbit,
delete!

import LinearAlgebra: norm, mul!, copy_oftype, copymutable_oftype
import LinearAlgebra: norm
import SpecialFunctions: erf, erfc, erfi, erfcx

using GTPSA_jll, Printf, PrettyTables
Expand Down Expand Up @@ -103,9 +103,12 @@ export
polar,
rect,
clear!,
mul!,

# Monomial as TPS creators:
# TPS ctors:
@vars,
@params,

# DEPRECATED: Monomial as TPS creators:
vars,
params,
complexvars,
Expand Down Expand Up @@ -162,7 +165,8 @@ include("fastgtpsa/fastgtpsa.jl") # Definition of the @FastGTPSA macro
include("fastgtpsa/operators.jl") # TempTPS special math operators/functions
include("global.jl") # Global variables
include("getset.jl") # Indexing/slicing TPS, par, convenience getters (gradient, jacobian, hessian)
include("ctors.jl") # Convenience constructors (vars, params, mono)
include("ctors.jl") # @vars @params @mono macros
include("deprecated_ctors.jl") # (DEPRECATED) Convenience constructors (vars, params, mono)
include("show.jl") # Output
include("methods.jl") # Higher-level TPS functions (setTPS!, clear!, derivatives, integrals, evaluate, etc)

Expand Down
Loading

0 comments on commit 222291b

Please sign in to comment.