Merge pull request #11 from JuliaAlgebra/system

System

Author: saschatimme

docs/src/index.md

@@ -31,15 +31,16 @@ But this is note the fastest way possible. In order to achieve the best performa
 intermediate storage. For this we have [`GradientConfig`](@ref) and [`JacobianConfig`](@ref).
 For single polynomial the API is as follows
 ```julia
-cfg = GradientConfig(f) # this can be reused!
+cfg = config(f, x) # this can be reused!
 f(x) == evaluate(f, x, cfg)
 # We can also compute the gradient of f at x
 map(g -> g(x), ∇f) == gradient(f, x, cfg)
 ```
 
 We also have support for systems of polynomials:
 ```julia
-cfg = JacobianConfig([f, f]) # this can be reused!
+F = System([f, g])
+cfg = config(F, x) # this can be reused!
 [f(x), f(x)] == evaluate([f, f] x, cfg)
 # We can also compute the jacobian of [f, f] at x
 jacobian(f, x, cfg)
@@ -56,6 +57,6 @@ Make sure to also check out [`GradientDiffResult`](@ref) and [`JacobianDiffResul
 ## Safety notes
 
 !!! warning
-    For the evaluation multivariate variant of [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method)
-    is used. Due to that for polynomials with terms of degree over 43 we cannot guarantee
+    The current implementation is not numerically stable in the sense that
+    for polynomials with terms of degree over 43 we cannot guarantee
     an error of less than 1 [ULP](https://en.wikipedia.org/wiki/Unit_in_the_last_place).

docs/src/performance.md

@@ -1,8 +1,7 @@
 In order to achieve a fast evaluation we need to precompute some things and also preallocate
-intermediate storage. For this we have `GradientConfig` and `JacobianConfig`:
+intermediate storage. For this we have
 ```@docs
-GradientConfig
-JacobianConfig
+config
 ```
 
 ## Evaluation
@@ -19,6 +18,16 @@ jacobian
 jacobian!
 ```
 
+## Systems
+```@docs
+System
+```
+Systems have the additional functions
+```@docs
+evaluate_and_jacobian!
+evaluate_and_jacobian
+```
+
 
 ## DiffResults
 ```@docs

src/FixedPolynomials.jl

@@ -15,5 +15,6 @@ module FixedPolynomials
     include("convert_promote.jl")
     include("tables.jl")
     include("config.jl")
+    include("system.jl")
 
 end

src/config.jl

@@ -1,4 +1,4 @@
-export GradientConfig, GradientDiffResult, JacobianConfig, JacobianDiffResult,
+export GradientConfig, GradientDiffResult, config,
     gradient!, evaluate, evaluate!, jacobian, jacobian!, value
 
 const Index = UInt16
@@ -75,20 +75,9 @@ function PolyConfig(g::Polynomial{T}, reduced_exponents::Matrix{UInt16}, big_loo
         grad_monomials,
         reduced_exponents_delimiters,
         reduced_exponents_map,
-        zeros(promote_type(T, S), n))
+        zeros(typeof(one(T) * one(S) + one(T) * one(S)), n))
 end
 
-# function Base.deepcopy(cfg::PolyConfig)
-#     PolyConfig(
-#         deepcopy(cfg.monomials_delimiters),
-#         deepcopy(cfg.monomials),
-#         deepcopy(cfg.grad_monomials_delimiters),
-#         deepcopy(cfg.grad_monomials),
-#         deepcopy(cfg.reduced_exponents_delimiters),
-#         deepcopy(cfg.reduced_exponents_map),
-#         deepcopy(cfg.reduced_values))
-# end
-
 @inline function fillreduced_values!(
     cfg::PolyConfig{T},
     g::Polynomial,
@@ -255,6 +244,15 @@ function GradientConfig(f::Polynomial{T}, ::Type{S}) where {T, S}
     GradientConfig(poly, diffs, diffs_values)
 end
 
+"""
+    config(F::Polynomial, x)
+
+Construct a `GradientConfig` for the evaluation of `f` with values like `x`.
+"""
+function config(f::Polynomial{T}, x::AbstractVector{S}) where {S, T}
+    GradientConfig(f, typeof(one(T) * one(S) + one(T) * one(S)))
+end
+
 function differences(f::Polynomial{T}, ::Type{S}) where {T, S}
     exponents = f.exponents
     reduced_exponents = convert.(UInt16, max.(exponents .- 1, 0))
@@ -415,224 +413,3 @@ function gradient!(diffresult::GradientDiffResult, g::Polynomial, x::AbstractVec
     _gradient!(diffresult.grad, x, cfg.poly)
     diffresult
 end
-
-
-"""
-
-    JacobianConfig(F::Vector{Polynomial{T}}, [x::AbstractVector{S}])
-
-A data structure with which the jacobian of a `Vector` `F` of `Polynomial`s can be
-evaluated efficiently. Note that `x` is only used to determine the
-output type of `F(x)`.
-
-    JacobianConfig(F::Vector{Polynomial{T}}, [S])
-
-Instead of a vector `x` a type can also be given directly.
-"""
-mutable struct JacobianConfig{T}
-    polys::Vector{PolyConfig{T}}
-    differences::Matrix{UInt8}
-    differences_values::Matrix{T}
-end
-
-
-function JacobianConfig(f::Vector{Polynomial{T}}, ::AbstractArray{S}) where {T, S}
-    JacobianConfig(f, S)
-end
-JacobianConfig(f::Vector{Polynomial{T}}) where T = JacobianConfig(f, T)
-function JacobianConfig(F::Vector{Polynomial{T}}, ::Type{S}) where {T, S}
-    diffs, diffs_values, big_lookups, reduced_exponents = differences(F, S)
-    polys = broadcast(PolyConfig, F, reduced_exponents, big_lookups, S)
-
-    JacobianConfig(polys, diffs, diffs_values)
-end
-
-
-function differences(F::Vector{Polynomial{T}}, ::Type{S}) where {T, S}
-    reduced_exponents = map(F) do f
-        convert.(UInt16, max.(f.exponents .- 1, 0))
-    end
-    differences, big_lookups = computetables(reduced_exponents)
-    differences_values = convert.(promote_type(T, S), differences)
-
-    differences, differences_values, big_lookups, reduced_exponents
-end
-
-function Base.deepcopy(cfg::JacobianConfig)
-    JacobianConfig(
-        deepcopy(cfg.polys),
-        deepcopy(cfg.differences),
-        deepcopy(cfg.differences_values))
-end
-
-"""
-    evaluate(F, x, cfg::JacobianConfig [, precomputed=false])
-
-Evaluate the system `F` at `x` using the precomputated values in `cfg`.
-Note that this is usually signifcant faster than `map(f -> evaluate(f, x), F)`.
-The return vector is constructed using `similar(x, T)`.
-
-### Example
-```julia
-cfg = JacobianConfig(F)
-evaluate(F, x, cfg)
-```
-
-With `precomputed=true` we rely on the previous intermediate results in `cfg`. Therefore
-the result is only correct if you previouls called `evaluate`, or `jacobian` with the same
-`x`.
-"""
-function evaluate(G::Vector{<:Polynomial}, x::AbstractVector, cfg::JacobianConfig{T}, precomputed=false) where T
-    evaluate!(similar(x, T, length(G)), G, x, cfg, precomputed)
-end
-
-"""
-    evaluate!(u, F, x, cfg::JacobianConfig [, precomputed=false])
-
-Evaluate the system `F` at `x` using the precomputated values in `cfg`
-and store the result in `u`.
-Note that this is usually signifcant faster than `map!(u, f -> evaluate(f, x), F)`.
-
-### Example
-```julia
-cfg = JacobianConfig(F)
-evaluate!(u, F, x, cfg)
-```
-
-With `precomputed=true` we rely on the previous intermediate results in `cfg`. Therefore
-the result is only correct if you previouls called `evaluate`, or `jacobian` with the same
-`x`.
-"""
-function evaluate!(u::AbstractVector, G::Vector{<:Polynomial}, x::AbstractVector, cfg::JacobianConfig{T}, precomputed=false) where T
-    if !precomputed
-        fillvalues!(cfg.differences_values, x, cfg.differences)
-        for i=1:length(cfg.polys)
-            fillreduced_values!(cfg.polys[i], G[i], x, cfg.differences_values)
-            u[i] = _evaluate(x, cfg.polys[i])
-        end
-    else
-        for i=1:length(cfg.polys)
-            u[i] = _evaluate(x, cfg.polys[i])
-        end
-    end
-    u
-end
-
-"""
-    jacobian(u, F, x, cfg::JacobianConfig [, precomputed=false])
-
-Evaluate the jacobian of `F` at `x` using the precomputated values in `cfg`. The return
-matrix is constructed using `similar(x, T, m, n)`.
-
-### Example
-```julia
-cfg = JacobianConfig(F)
-jacobian(F, x, cfg)
-```
-
-With `precomputed=true` we rely on the previous intermediate results in `cfg`. Therefore
-the result is only correct if you previouls called `evaluate`, or `jacobian` with the same
-`x`.
-"""
-function jacobian(g::Vector{<:Polynomial}, x::AbstractVector, cfg::JacobianConfig{T}, precomputed=false) where T
-    u = similar(x, T, (length(g), length(x)))
-    jacobian!(u, g, x, cfg, precomputed)
-    u
-end
-
-"""
-    jacobian!(u, F, x, cfg::JacobianConfig [, precomputed=false])
-
-Evaluate the jacobian of `F` at `x` using the precomputated values in `cfg`
-and store the result in `u`.
-
-### Example
-```julia
-cfg = JacobianConfig(F)
-jacobian!(u, F, x, cfg)
-```
-
-With `precomputed=true` we rely on the previous intermediate results in `cfg`. Therefore
-the result is only correct if you previouls called `evaluate`, or `jacobian` with the same
-`x`.
-"""
-function jacobian!(u::AbstractMatrix, G::Vector{<:Polynomial}, x::AbstractVector, cfg::JacobianConfig{T}, precomputed=false) where T
-    if !precomputed
-        fillvalues!(cfg.differences_values, x, cfg.differences)
-        for i=1:length(G)
-            fillreduced_values!(cfg.polys[i], G[i], x, cfg.differences_values)
-            gradient_row!(u, x, cfg.polys[i], i)
-        end
-    else
-        for i=1:length(G)
-            gradient_row!(u, x, cfg.polys[i], i)
-        end
-    end
-    u
-end
-
-"""
-    JacobianDiffResult(cfg::GradientConfig)
-
-During the computation of the jacobian ``J_F(x)`` we compute nearly everything we need for the evaluation of
-``F(x)``. `JacobianDiffResult` allocates memory to hold both values.
-This structure also signals `jacobian!` to store ``F(x)`` and ``J_F(x)``.
-
-### Example
-
-```julia
-cfg = JacobianConfig(F, x)
-r = JacobianDiffResult(cfg)
-jacobian!(r, F, x, cfg)
-
-value(r) == map(f -> f(x), F)
-jacobian(r) == jacobian(F, x, cfg)
-```
-
-    JacobianDiffResult(value::AbstractVector, jacobian::AbstractMatrix)
-
-Allocate the memory to hold the value and the jacobian by yourself.
-"""
-mutable struct JacobianDiffResult{T, AV<:AbstractVector{T}, AM<:AbstractMatrix{T}}
-    value::AV
-    jacobian::AM
-end
-
-function JacobianDiffResult(cfg::JacobianConfig{T}) where T
-    JacobianDiffResult{T, Vector{T}, Matrix{T}}(
-        zeros(T, length(cfg.polys)),
-        zeros(T, length(cfg.polys), size(cfg.differences, 1)))
-end
-
-function JacobianDiffResult(value::AbstractVector{T}, jacobian::AbstractMatrix{T}) where T
-    JacobianDiffResult{T, typeof(value), typeof(jacobian)}(value, jacobian)
-end
-value(r::JacobianDiffResult) = r.value
-jacobian(r::JacobianDiffResult) = r.jacobian
-
-"""
-    jacobian!(r::JacobianDiffResult, F, x, cfg::JacobianConfig)
-
-Compute ``F(x)`` and the jacobian of `F` at `x` at once using the precomputated values in `cfg`
-and store thre result in `r`. This is faster than computing both values separetely.
-
-### Example
-```julia
-cfg = GradientConfig(g)
-r = GradientDiffResult(cfg)
-gradient!(r, g, x, cfg)
-
-value(r) == g(x)
-gradient(r) == gradient(g, x, cfg)
-```
-"""
-function jacobian!(r::JacobianDiffResult, G::Vector{<:Polynomial}, x::AbstractVector, cfg::JacobianConfig{T}) where T
-    fillvalues!(cfg.differences_values, x, cfg.differences)
-    for i=1:length(G)
-        fillreduced_values!(cfg.polys[i], G[i], x, cfg.differences_values)
-        gradient_row!(r.jacobian, x, cfg.polys[i], i)
-        r.value[i] = _evaluate(x, cfg.polys[i])
-    end
-
-    r
-end