From 9b2ce0df16b01d1e64496a2a3051f7874efa2bb5 Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 7 Feb 2020 21:39:18 +0100 Subject: [PATCH 1/2] add IntervalMatrixPower wrapper --- docs/src/lib/methods.md | 8 +++ docs/src/lib/types.md | 6 ++ src/IntervalMatrices.jl | 8 +++ src/power.jl | 135 ++++++++++++++++++++++++++++++++++++++++ 4 files changed, 157 insertions(+) create mode 100644 src/power.jl diff --git a/docs/src/lib/methods.md b/docs/src/lib/methods.md index c3036eb..c8e9f38 100644 --- a/docs/src/lib/methods.md +++ b/docs/src/lib/methods.md @@ -28,6 +28,14 @@ scale! ⊆ ``` +## Matrix power + +```@docs +increment! +increment +get +``` + ## Exponentiation ```@docs diff --git a/docs/src/lib/types.md b/docs/src/lib/types.md index fefa74a..64348cc 100644 --- a/docs/src/lib/types.md +++ b/docs/src/lib/types.md @@ -22,3 +22,9 @@ AbstractIntervalMatrix ```@docs IntervalMatrix ``` + +## Interval-matrix-power wrapper + +```@docs +IntervalMatrixPower +``` diff --git a/src/IntervalMatrices.jl b/src/IntervalMatrices.jl index 2dbc9cb..a1b0ac6 100644 --- a/src/IntervalMatrices.jl +++ b/src/IntervalMatrices.jl @@ -19,6 +19,11 @@ include("matrix.jl") # ================================= include("arithmetic.jl") +# ================================================== +# Methods to compute the power of an interval matrix +# ================================================== +include("power.jl") + # ======================================================= # Methods to handle the exponential of an interval matrix # ======================================================= @@ -44,4 +49,7 @@ export inf, sup, exp_overapproximation, exp_underapproximation +export IntervalMatrixPower, + increment!, increment + end # module diff --git a/src/power.jl b/src/power.jl new file mode 100644 index 0000000..ddcfb0b --- /dev/null +++ b/src/power.jl @@ -0,0 +1,135 @@ +import Base: copy, get + +""" + IntervalMatrixPower{T} + +A wrapper for the matrix power that can be incremented. + +### Fields + +- `M` -- the original matrix +- `Mᵏ` -- the current matrix power, i.e., `M^k` +- `k` -- the current power index + +### Notes + +The wrapper should only be accessed using the interface functions. +The internal representation (such as the fields) are subject to future changes. + +### Examples + +```jldoctest +julia> A = IntervalMatrix([Interval(2.0) Interval(2.0, 3.0); Interval(0.0) Interval(-1.0, 1.0)]) +2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}: + [2, 2] [2, 3] + [0, 0] [-1, 1] + +julia> pow = IntervalMatrixPower(A); + +julia> increment!(pow) +2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}: + [4, 4] [2, 9] + [0, 0] [0, 1] + +julia> increment(pow) +2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}: + [8, 8] [-1, 21] + [0, 0] [-1, 1] + +julia> get(pow) +2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}: + [4, 4] [2, 9] + [0, 0] [0, 1] + +``` +""" +mutable struct IntervalMatrixPower{T} + M::IntervalMatrix{T} + Mᵏ::IntervalMatrix{T} + k::Int + + function IntervalMatrixPower(M::IntervalMatrix{T}) where {T} + return new{T}(M, M, 1) + end + + function IntervalMatrixPower(M::IntervalMatrix{T}, Mᵏ::IntervalMatrix{T}, + k::Int) where {T} + return new{T}(M, Mᵏ, k) + end +end + +function copy(pow::IntervalMatrixPower) + return IntervalMatrixPower(pow.M, pow.Mᵏ, pow.k) +end + +""" + increment!(pow::IntervalMatrixPower) + +Increment a matrix power. + +### Input + +- `pow` -- wrapper of a matrix power (modified in this function) + +### Output + +The next matrix power now, reflected in the modified wrapper. +""" +function increment!(pow::IntervalMatrixPower) + pow.k += 1 + if _isapoweroftwo(pow.k) + pow.Mᵏ = _eval_poweroftwo(pow.M, pow.k) + else + pow.Mᵏ *= pow.M + end + return pow.Mᵏ +end + +""" + increment!(pow::IntervalMatrixPower) + +Increment a matrix power. + +### Input + +- `pow` -- wrapper of a matrix power + +### Output + +The next matrix power now. +""" +function increment(pow::IntervalMatrixPower) + return increment!(copy(pow)) +end + +# checks whether a number is a power of 2 +function _isapoweroftwo(k::Integer) + # see https://stackoverflow.com/a/600306 + return (k & (k - 1)) == 0 +end + +# compute the (exact) matrix power for a power of two +function _eval_poweroftwo(M::IntervalMatrix, k::Integer) + Mᵏ = M + @inbounds for i in 1:Int(log(2, k)) + Mᵏ = square(Mᵏ) + end + return Mᵏ +end + +""" + get(pow::IntervalMatrixPower) + +Return the matrix represented by a wrapper of a matrix power. + +### Input + +- `pow` -- wrapper of a matrix power + +### Output + +The matrix power represented by the wrapper. +""" +function get(pow::IntervalMatrixPower) + return pow.Mᵏ +end From 664bb3305d70d9b6c1214305f78daffedd396b18 Mon Sep 17 00:00:00 2001 From: schillic Date: Fri, 7 Feb 2020 21:59:59 +0100 Subject: [PATCH 2/2] review --- docs/src/lib/methods.md | 2 ++ src/IntervalMatrices.jl | 4 +++- src/power.jl | 52 +++++++++++++++++++++++++++++++++++++---- 3 files changed, 53 insertions(+), 5 deletions(-) diff --git a/docs/src/lib/methods.md b/docs/src/lib/methods.md index c8e9f38..fb93d57 100644 --- a/docs/src/lib/methods.md +++ b/docs/src/lib/methods.md @@ -34,6 +34,8 @@ scale! increment! increment get +base +index ``` ## Exponentiation diff --git a/src/IntervalMatrices.jl b/src/IntervalMatrices.jl index a1b0ac6..71d9622 100644 --- a/src/IntervalMatrices.jl +++ b/src/IntervalMatrices.jl @@ -50,6 +50,8 @@ export inf, sup, exp_underapproximation export IntervalMatrixPower, - increment!, increment + increment!, increment, + base, + index end # module diff --git a/src/power.jl b/src/power.jl index ddcfb0b..efd13b5 100644 --- a/src/power.jl +++ b/src/power.jl @@ -8,7 +8,7 @@ A wrapper for the matrix power that can be incremented. ### Fields - `M` -- the original matrix -- `Mᵏ` -- the current matrix power, i.e., `M^k` +- `Mᵏ` -- the current matrix power, i.e., ``M^k`` - `k` -- the current power index ### Notes @@ -41,6 +41,14 @@ julia> get(pow) [4, 4] [2, 9] [0, 0] [0, 1] +julia> index(pow) +2 + +julia> base(pow) +2×2 IntervalMatrix{Float64,Interval{Float64},Array{Interval{Float64},2}}: + [2, 2] [2, 3] + [0, 0] [-1, 1] + ``` """ mutable struct IntervalMatrixPower{T} @@ -54,6 +62,7 @@ mutable struct IntervalMatrixPower{T} function IntervalMatrixPower(M::IntervalMatrix{T}, Mᵏ::IntervalMatrix{T}, k::Int) where {T} + @assert k >= 1 "matrix powers must be positive" return new{T}(M, Mᵏ, k) end end @@ -65,7 +74,7 @@ end """ increment!(pow::IntervalMatrixPower) -Increment a matrix power. +Increment a matrix power in-place (i.e., storing the result in `pow`). ### Input @@ -86,9 +95,9 @@ function increment!(pow::IntervalMatrixPower) end """ - increment!(pow::IntervalMatrixPower) + increment(pow::IntervalMatrixPower) -Increment a matrix power. +Increment a matrix power without modifying `pow`. ### Input @@ -117,6 +126,24 @@ function _eval_poweroftwo(M::IntervalMatrix, k::Integer) return Mᵏ end +""" + base(pow::IntervalMatrixPower) + +Return the original matrix represented by a wrapper of a matrix power. + +### Input + +- `pow` -- wrapper of a matrix power + +### Output + +The matrix ``M`` being the basis of the matrix power ``M^k`` represented by the +wrapper. +""" +function base(pow::IntervalMatrixPower) + return pow.M +end + """ get(pow::IntervalMatrixPower) @@ -133,3 +160,20 @@ The matrix power represented by the wrapper. function get(pow::IntervalMatrixPower) return pow.Mᵏ end + +""" + index(pow::IntervalMatrixPower) + +Return the current index of the wrapper of a matrix power. + +### Input + +- `pow` -- wrapper of a matrix power + +### Output + +The index `k` of the wrapper representing ``M^k``. +""" +function index(pow::IntervalMatrixPower) + return pow.k +end