From 9fd3876e0d38c8b65133a7589f915e356c38427f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 02:40:52 +0100 Subject: [PATCH 01/29] Triangulation transfer operator estimators Also, move to generator-style estimaton --- Project.toml | 9 + docs/Project.toml | 3 +- docs/make.jl | 1 + docs/src/TransferOperator.md | 2 + src/Entropies.jl | 13 +- .../{rectangular => }/GroupSlices.jl | 0 .../{rectangular => }/binning_schemes.jl | 0 src/binning_based/binningbased.jl | 3 + .../rectangular/TransferOperator.jl | 464 ------------------ .../rectangular/rectangular_estimators.jl | 9 - src/binning_based/rectangular_estimators.jl | 6 + src/binning_based/transferoperator.jl | 131 +++++ src/binning_based/transferoperator/common.jl | 324 ++++++++++++ .../rectangular/rectangular.jl | 208 ++++++++ .../transferoperator/transfer_operator.jl | 3 + ...bstractDelaunayTriangulation-checkpoint.jl | 37 ++ .../DelaunayTriangulation-checkpoint.jl | 20 + .../DelaunayTriangulations-checkpoint.jl | 13 + .../.ipynb_checkpoints/addnoise-checkpoint.jl | 105 ++++ .../AbstractDelaunayTriangulation.jl | 37 ++ .../DelaunayTriangulation.jl | 20 + .../DelaunayTriangulations.jl | 13 + .../delaunay_triangulations/addnoise.jl | 105 ++++ .../DelaunayTriangulation-checkpoint.jl | 21 + .../plot_recipes-checkpoint.jl | 60 +++ .../plot_recipes/plot_recipes.jl | 60 +++ .../SimplexExact-checkpoint.jl | 154 ++++++ .../triangular/exact/SimplexExact.jl | 155 ++++++ .../SimplexPoint-checkpoint.jl | 131 +++++ .../helper_functions-checkpoint.jl | 175 +++++++ .../triangular/point/SimplexPoint.jl | 131 +++++ .../triangular/point/helper_functions.jl | 178 +++++++ .../simplex_subsampling-checkpoint.jl | 116 +++++ .../simplex_types-checkpoint.jl | 5 + .../simplex_types/AbstractSimplex.jl | 98 ++++ .../simplex_types/ImmutableSimplex.jl | 1 + .../simplex_types/MutableSimplex.jl | 60 +++ .../simplex_types/simplex_subsampling.jl | 116 +++++ .../triangular/simplex_types/simplex_types.jl | 5 + src/binning_based/transferoperator/utils.jl | 19 + .../VisitationFrequency.jl | 5 +- .../count_box_visits.jl | 0 .../histogram_estimation.jl | 0 43 files changed, 2540 insertions(+), 476 deletions(-) rename src/binning_based/{rectangular => }/GroupSlices.jl (100%) rename src/binning_based/{rectangular => }/binning_schemes.jl (100%) create mode 100644 src/binning_based/binningbased.jl delete mode 100644 src/binning_based/rectangular/TransferOperator.jl delete mode 100644 src/binning_based/rectangular/rectangular_estimators.jl create mode 100644 src/binning_based/rectangular_estimators.jl create mode 100644 src/binning_based/transferoperator.jl create mode 100644 src/binning_based/transferoperator/common.jl create mode 100644 src/binning_based/transferoperator/rectangular/rectangular.jl create mode 100644 src/binning_based/transferoperator/transfer_operator.jl create mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/AbstractDelaunayTriangulation-checkpoint.jl create mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulation-checkpoint.jl create mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulations-checkpoint.jl create mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/addnoise-checkpoint.jl create mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/AbstractDelaunayTriangulation.jl create mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/DelaunayTriangulation.jl create mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/DelaunayTriangulations.jl create mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/addnoise.jl create mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/composite_types/.ipynb_checkpoints/DelaunayTriangulation-checkpoint.jl create mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/plot_recipes/.ipynb_checkpoints/plot_recipes-checkpoint.jl create mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/plot_recipes/plot_recipes.jl create mode 100644 src/binning_based/transferoperator/triangular/exact/.ipynb_checkpoints/SimplexExact-checkpoint.jl create mode 100644 src/binning_based/transferoperator/triangular/exact/SimplexExact.jl create mode 100644 src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/SimplexPoint-checkpoint.jl create mode 100644 src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/helper_functions-checkpoint.jl create mode 100644 src/binning_based/transferoperator/triangular/point/SimplexPoint.jl create mode 100644 src/binning_based/transferoperator/triangular/point/helper_functions.jl create mode 100644 src/binning_based/transferoperator/triangular/simplex_types/.ipynb_checkpoints/simplex_subsampling-checkpoint.jl create mode 100644 src/binning_based/transferoperator/triangular/simplex_types/.ipynb_checkpoints/simplex_types-checkpoint.jl create mode 100644 src/binning_based/transferoperator/triangular/simplex_types/AbstractSimplex.jl create mode 100644 src/binning_based/transferoperator/triangular/simplex_types/ImmutableSimplex.jl create mode 100644 src/binning_based/transferoperator/triangular/simplex_types/MutableSimplex.jl create mode 100644 src/binning_based/transferoperator/triangular/simplex_types/simplex_subsampling.jl create mode 100644 src/binning_based/transferoperator/triangular/simplex_types/simplex_types.jl create mode 100644 src/binning_based/transferoperator/utils.jl rename src/binning_based/{rectangular => visitation_frequency}/VisitationFrequency.jl (89%) rename src/binning_based/{rectangular => visitation_frequency}/count_box_visits.jl (100%) rename src/binning_based/{rectangular => visitation_frequency}/histogram_estimation.jl (100%) diff --git a/Project.toml b/Project.toml index 7f9f64045..eeed0acd9 100644 --- a/Project.toml +++ b/Project.toml @@ -7,18 +7,27 @@ version = "0.11.1" [deps] DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" +DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" +PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" +Requires = "ae029012-a4dd-5104-9daa-d747884805df" +Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" [compat] DelayEmbeddings = "1.10" Distances = "0.9, 0.10" NearestNeighbors = "0.4" +Requires = "^1.1" SpecialFunctions = "0.10, 1.0" StaticArrays = "0.12, 1.0" Wavelets = "0.9" diff --git a/docs/Project.toml b/docs/Project.toml index 5dce982d4..b9014ff86 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,12 +1,13 @@ [deps] Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" -Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" +Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" diff --git a/docs/make.jl b/docs/make.jl index 03d18dff6..e09f497f1 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -11,6 +11,7 @@ using Entropies using PyPlot using DynamicalSystems using Wavelets +using Simplices # %% JuliaDynamics theme. # download the themes diff --git a/docs/src/TransferOperator.md b/docs/src/TransferOperator.md index 0f08bf479..7a2e2bc95 100644 --- a/docs/src/TransferOperator.md +++ b/docs/src/TransferOperator.md @@ -2,6 +2,8 @@ ```@docs TransferOperator +SimplexExact +SimplexPoint ``` ## Utility methods/types diff --git a/src/Entropies.jl b/src/Entropies.jl index 8c5662747..0855ab627 100644 --- a/src/Entropies.jl +++ b/src/Entropies.jl @@ -1,10 +1,21 @@ + module Entropies + using Requires + include("core.jl") include("histogram_estimation.jl") include("counting_based/CountOccurrences.jl") include("symbolic/symbolic.jl") - include("binning_based/rectangular/rectangular_estimators.jl") include("kerneldensity/kerneldensity.jl") include("wavelet/wavelet.jl") include("nearest_neighbors/nearest_neighbors.jl") + include("binning_based/rectangular_estimators.jl") + + function __init__() + @require Simplices="d5428e67-3037-59ba-9ab1-57a04f0a3b6a" begin + using .Simplices + include("binning_based/transferoperator/triangular/exact/SimplexExact.jl") + include("binning_based/transferoperator/triangular/point/SimplexPoint.jl") + end + end end diff --git a/src/binning_based/rectangular/GroupSlices.jl b/src/binning_based/GroupSlices.jl similarity index 100% rename from src/binning_based/rectangular/GroupSlices.jl rename to src/binning_based/GroupSlices.jl diff --git a/src/binning_based/rectangular/binning_schemes.jl b/src/binning_based/binning_schemes.jl similarity index 100% rename from src/binning_based/rectangular/binning_schemes.jl rename to src/binning_based/binning_schemes.jl diff --git a/src/binning_based/binningbased.jl b/src/binning_based/binningbased.jl new file mode 100644 index 000000000..e600b3fbd --- /dev/null +++ b/src/binning_based/binningbased.jl @@ -0,0 +1,3 @@ +abstract type BinningProbabilitiesEstimator <: ProbabilitiesEstimator end +include("rectangular/rectangular_estimators.jl") +#include("triangular/TriangulationEstimators.jl") \ No newline at end of file diff --git a/src/binning_based/rectangular/TransferOperator.jl b/src/binning_based/rectangular/TransferOperator.jl deleted file mode 100644 index 291f7d139..000000000 --- a/src/binning_based/rectangular/TransferOperator.jl +++ /dev/null @@ -1,464 +0,0 @@ -using DelayEmbeddings, SparseArrays -include("GroupSlices.jl") - -export - TransferOperator, # the probabilities estimator - InvariantMeasure, invariantmeasure, - transfermatrix - -""" - TransferOperator(ϵ::RectangularBinning) <: BinningProbabilitiesEstimator - -A probability estimator based on binning data into rectangular boxes dictated by -the binning scheme `ϵ`, then approxmating the transfer (Perron-Frobenius) operator -over the bins, the taking the invariant measure associated with that transfer operator -as the bin probabilities. Assumes that the input data are sequential. - -This implementation follows the grid estimator approach in Diego et al. (2019)[^Diego2019]. - -## Description - -The transfer operator ``P^{N}``is computed as an `N`-by-`N` matrix of transition -probabilities between the states defined by the partition elements, where `N` is the -number of boxes in the partition that is visited by the orbit/points. - -If ``\\{x_t^{(D)} \\}_{n=1}^L`` are the ``L`` different ``D``-dimensional points over -which the transfer operator is approximated, ``\\{ C_{k=1}^N \\}`` are the ``N`` different -partition elements (as dictated by `ϵ`) that gets visited by the points, and - ``\\phi(x_t) = x_{t+1}``, then - -```math -P_{ij} = \\dfrac -{\\#\\{ x_n | \\phi(x_n) \\in C_j \\cap x_n \\in C_i \\}} -{\\#\\{ x_m | x_m \\in C_i \\}}, -``` - -where ``\\#`` denotes the cardinal. The element ``P_{ij}`` thus indicates how many points -that are initially in box ``C_i`` end up in box ``C_j`` when the points in ``C_i`` are -projected one step forward in time. Thus, the row ``P_{ik}^N`` where -``k \\in \\{1, 2, \\ldots, N \\}`` gives the probability -of jumping from the state defined by box ``C_i`` to any of the other ``N`` states. It -follows that ``\\sum_{k=1}^{N} P_{ik} = 1`` for all ``i``. Thus, ``P^N`` is a row/right -stochastic matrix. - -### Invariant measure estimation from transfer operator - -The left invariant distribution ``\\mathbf{\\rho}^N`` is a row vector, where -``\\mathbf{\\rho}^N P^{N} = \\mathbf{\\rho}^N``. Hence, ``\\mathbf{\\rho}^N`` is a row -eigenvector of the transfer matrix ``P^{N}`` associated with eigenvalue 1. The distribution -``\\mathbf{\\rho}^N`` approximates the invariant density of the system subject to the -partition `ϵ`, and can be taken as a probability distribution over the partition elements. - -In practice, the invariant measure ``\\mathbf{\\rho}^N`` is computed using -[`invariantmeasure`](@ref), which also approximates the transfer matrix. The invariant distribution -is initialized as a length-`N` random distribution which is then applied to ``P^{N}``. -The resulting length-`N` distribution is then applied to ``P^{N}`` again. This process -repeats until the difference between the distributions over consecutive iterations is -below some threshold. - -## Probability and entropy estimation - -- `probabilities(x::AbstractDataset, est::TransferOperator{RectangularBinning})` estimates - probabilities for the bins defined by the provided binning (`est.ϵ`) -- `genentropy(x::AbstractDataset, est::TransferOperator{RectangularBinning})` does the same, - but computes generalized entropy using the probabilities. - - -See also: [`RectangularBinning`](@ref), [`invariantmeasure`](@ref). - -[^Diego2019]: Diego, D., Haaga, K. A., & Hannisdal, B. (2019). Transfer entropy computation using the Perron-Frobenius operator. Physical Review E, 99(4), 042212. -""" -struct TransferOperator{R} <: BinningProbabilitiesEstimator - ϵ::R - - function TransferOperator(ϵ::R) where R <: RectangularBinning - new{R}(ϵ) - end -end - -# If x is not sorted, we need to look at all pairwise comparisons -function inds_in_terms_of_unique(x) - U = unique(x) - N = length(x) - Nu = length(U) - inds = zeros(Int, N) - - for j = 1:N - xⱼ = view(x, j) - for i = 1:Nu - # using views doesn't allocate - @inbounds if xⱼ == view(U, i) - inds[j] = i - end - end - end - - return inds -end - -# Taking advantage of the fact that x is sorted reduces runtime by 1.5 orders of magnitude -# for datasets of >100 000+ points -function inds_in_terms_of_unique_sorted(x) # assumes sorted - @assert issorted(x) - U = unique(x) - N, Nu = length(x), length(U) - prev = view(x, 1) - inds = zeros(Int, N) - uidx = 1 - @inbounds for j = 1:N - xⱼ = view(x, j) - # if the current value has changed, then we know that the corresponding index - # for the unique point must be incremented by 1 - if xⱼ != prev - prev = xⱼ - uidx += 1 - end - inds[j] = uidx - end - - return inds -end - -function inds_in_terms_of_unique(x, sorted::Bool) - if sorted - return inds_in_terms_of_unique_sorted(x) - else - return inds_in_terms_of_unique(x) - end -end - -inds_in_terms_of_unique(x::AbstractDataset) = inds_in_terms_of_unique(x.data) - -""" - TransferOperatorApproximationRectangular(to, ϵ::RectangularBinning, mini, edgelengths, - bins, sort_idxs) - -The `N`-by-`N` matrix `to` is an approximation to the transfer operator, subject to the -partition `ϵ`, computed over some set of sequentially ordered points. - -For convenience, `mini` and `edgelengths` provide the minima and box edge lengths along -each coordinate axis, as determined by applying `ϵ` to the points. The coordinates of -the (leftmost, if axis is ordered low-high) box corners are given in `bins`. - -Only bins actually visited by the points are considered, and `bins` give the coordinates -of these bins. The element `bins[i]` correspond to the `i`-th state of the system, which -corresponds to the `i`-th column/row of the transfer operator `to`. - -`sort_idxs` contains the indices that would sort the input points. `visitors` is a -vector of vectors, where `visitors[i]` contains the indices of the (sorted) -points that visits `bins[i]`. -See also: [`RectangularBinning`](@ref). -""" -struct TransferOperatorApproximationRectangular{T<:Real} - transfermatrix::AbstractArray{T, 2} - ϵ::RectangularBinning - mini - edgelengths - bins - sort_idxs::Vector{Int} - visitors::Vector{Vector{Int}} -end - -""" - transferoperator(pts::AbstractDataset{D, T}, ϵ::RectangularBinning) → TransferOperatorApproximationRectangular - -Estimate the transfer operator given a set of sequentially ordered points subject to a -rectangular partition given by `ϵ`. - -## Example - -```julia -using DynamicalSystems, Plots, Entropy -D = 4 -ds = Systems.lorenz96(D; F = 32.0) -N, dt = 20000, 0.1 -orbit = trajectory(ds, N*dt; dt = dt, Ttr = 10.0) - -# Estimate transfer operator over some coarse graining of the orbit. -transferoperator(orbit, RectangularBinning(10)) -``` - -See also: [`RectangularBinning`](@ref). -""" -function transferoperator(pts::AbstractDataset{D, T}, ϵ::RectangularBinning; - boundary_condition = :circular) where {D, T<:Real} - - L = length(pts) - mini, edgelengths = Entropies.minima_edgelengths(pts, ϵ) - - # The L points visits a total of L bins, which are the following bins: - visited_bins = Entropies.encode_as_bin(pts, mini, edgelengths) - sort_idxs = sortperm(visited_bins) - - # TODO: fix re-indexing after sorting. Sorting is much faster, so we want to do so. - #sort!(visited_bins) - - # There are N=length(unique(visited_bins)) unique bins. - # Which of the unqiue bins does each of the L points visit? - visits_whichbin = inds_in_terms_of_unique(visited_bins, false) # set to true when sorting is fixed - - # `visitors` lists the indices of the points visiting each of the N unique bins. - slices = GroupSlices.groupslices(visited_bins) - visitors = GroupSlices.groupinds(slices) - - # first_visited_by == [x[1] for x in visitors] - first_visited_by = GroupSlices.firstinds(slices) - L = length(first_visited_by) - - I = Int32[] - J = Int32[] - P = Float64[] - - # Preallocate target index for the case where there is only - # one point of the orbit visiting a bin. - target_bin_j::Int = 0 - n_visitsᵢ::Int = 0 - - if boundary_condition == :circular - #warn("Using circular boundary condition") - append!(visits_whichbin, [1]) - elseif boundary_condition == :random - #warn("Using random circular boundary condition") - append!(visits_whichbin, [rand(1:length(visits_whichbin))]) - else - error("Boundary condition $(boundary_condition) not implemented") - end - - # Loop over the visited bins bᵢ - for i in 1:L - # How many times is this bin visited? - n_visitsᵢ = length(visitors[i]) - - # If both conditions below are true, then there is just one - # point visiting the i-th bin. If there is only one visiting point and - # it happens to be the last, we skip it, because we don't know its - # image. - if n_visitsᵢ == 1 && !(i == visits_whichbin[end]) - # To which bin does the single point visiting bᵢ jump if we - # shift it one time step ahead along its orbit? - target_bin_j = visits_whichbin[visitors[i][1] + 1][1] - - # We now know that exactly one point (the i-th) does the - # transition from i to the target j. - push!(I, i) - push!(J, target_bin_j) - push!(P, 1.0) - end - - # If more than one point of the orbit visits the i-th bin, we - # identify the visiting points and track which bins bⱼ they end up - # in after the forward linear map of the points. - if n_visitsᵢ > 1 - timeindices_visiting_pts = visitors[i] - - # TODO: Introduce circular boundary condition. Simply excluding - # might lead to a cascade of loosing points. - - # If bᵢ is the bin visited by the last point in the orbit, then - # the last entry of `visiting_pts` will be the time index of the - # last point of the orbit. In the next time step, that point will - # have nowhere to go along its orbit (precisely because it is the - # last data point). Thus, we exclude it. - if i == visits_whichbin[end] - #warn("Removing last point") - n_visitsᵢ = length(timeindices_visiting_pts) - 1 - timeindices_visiting_pts = timeindices_visiting_pts[1:(end - 1)] - end - - # To which boxes do each of the visitors to bᵢ jump in the next - # time step? - target_bins = visits_whichbin[timeindices_visiting_pts .+ 1] - unique_target_bins = unique(target_bins) - - # Count how many points jump from the i-th bin to each of - # the unique target bins, and use that to calculate the transition - # probability from bᵢ to bⱼ. - #for j in 1:length(unique_target_bins) - for (j, bᵤ) in enumerate(unique(target_bins)) - n_transitions_i_to_j = sum(target_bins .== bᵤ) - - push!(I, i) - push!(J, bᵤ) - push!(P, n_transitions_i_to_j / n_visitsᵢ) - end - end - end - - # Transfer operator is just the normalized transition probabilities between the boxes. - TO = sparse(I, J, P) - - # Compute the coordinates of the visited bins. bins[i] corresponds to the i-th - # row/column of the transfer operator - unique!(visited_bins) - bins = [β .* edgelengths .+ mini for β in visited_bins] - - TransferOperatorApproximationRectangular(TO, ϵ, mini, edgelengths, bins, - sort_idxs, visitors) -end - -""" - InvariantMeasure(to, ρ) - -Minimal return struct for [`invariantmeasure`](@ref) that contains the estimated invariant -measure `ρ`, as well as the transfer operator `to` from which it is computed (including -bin information). - -See also: [`invariantmeasure`](@ref). -""" -struct InvariantMeasure{T} - to::T - ρ::Probabilities - - function InvariantMeasure(to::T, ρ) where T - new{T}(to, ρ) - end -end - -function invariantmeasure(iv::InvariantMeasure) - return iv.ρ, iv.to.bins -end - - -import LinearAlgebra: norm -""" - invariantmeasure(x::AbstractDataset, ϵ::RectangularBinning) → iv::InvariantMeasure - -Estimate an invariant measure over the points in `x` based on binning the data into -rectangular boxes dictated by the binning scheme `ϵ`, then approximate the transfer -(Perron-Frobenius) operator over the bins. From the approximation to the transfer operator, -compute an invariant distribution over the bins. Assumes that the input data are sequential. - -Details on the estimation procedure is found the [`TransferOperator`](@ref) docstring. - -## Example - -```julia -using DynamicalSystems, Plots, Entropies -D = 4 -ds = Systems.lorenz96(D; F = 32.0) -N, dt = 20000, 0.1 -orbit = trajectory(ds, N*dt; dt = dt, Ttr = 10.0) - -# Estimate the invariant measure over some coarse graining of the orbit. -iv = invariantmeasure(orbit, RectangularBinning(15)) - -# Get the probabilities and bins -invariantmeasure(iv) -``` - -## Probabilities and bin information - - invariantmeasure(iv::InvariantMeasure) → (ρ::Probabilities, bins::Vector{<:SVector}) - -From a pre-computed invariant measure, return the probabilities and associated bins. -The element `ρ[i]` is the probability of visitation to the box `bins[i]`. Analogous to -[`binhist`](@ref). - - -!!! hint "Transfer operator approach vs. naive histogram approach" - - Why bother with the transfer operator instead of using regular histograms to obtain - probabilities? - - In fact, the naive histogram approach and the - transfer operator approach are equivalent in the limit of long enough time series - (as ``n \\to \\intfy``), which is guaranteed by the ergodic theorem. There is a crucial - difference, however: - - The naive histogram approach only gives the long-term probabilities that - orbits visit a certain region of the state space. The transfer operator encodes that - information too, but comes with the added benefit of knowing the *transition - probabilities* between states (see [`transfermatrix`](@ref)). - -See also: [`InvariantMeasure`](@ref). -""" -function invariantmeasure(to::TransferOperatorApproximationRectangular; - N::Int = 200, tolerance::Float64 = 1e-8, delta::Float64 = 1e-8) - - TO = to.transfermatrix - #= - # Start with a random distribution `Ρ` (big rho). Normalise it so that it - # sums to 1 and forms a true probability distribution over the partition elements. - =# - Ρ = rand(Float64, 1, size(to.transfermatrix, 1)) - Ρ = Ρ ./ sum(Ρ, dims = 2) - - #= - # Start estimating the invariant distribution. We could either do this by - # finding the left-eigenvector of M, or by repeated application of M on Ρ - # until the distribution converges. Here, we use the latter approach, - # meaning that we iterate until Ρ doesn't change substantially between - # iterations. - =# - distribution = Ρ * to.transfermatrix - - distance = norm(distribution - Ρ) / norm(Ρ) - - check = floor(Int, 1 / delta) - check_pts = floor.(Int, transpose(collect(1:N)) ./ check) .* transpose(collect(1:N)) - check_pts = check_pts[check_pts .> 0] - num_checkpts = size(check_pts, 1) - check_pts_counter = 1 - - counter = 1 - while counter <= N && distance >= tolerance - counter += 1 - Ρ = distribution - - # Apply the Markov matrix to the current state of the distribution - distribution = Ρ * to.transfermatrix - - if (check_pts_counter <= num_checkpts && - counter == check_pts[check_pts_counter]) - - check_pts_counter += 1 - colsum_distribution = sum(distribution, dims = 2)[1] - if abs(colsum_distribution - 1) > delta - distribution = distribution ./ colsum_distribution - end - end - - distance = norm(distribution - Ρ) / norm(Ρ) - end - distribution = dropdims(distribution, dims = 1) - - # Do the last normalisation and check - colsum_distribution = sum(distribution) - - if abs(colsum_distribution - 1) > delta - distribution = distribution ./ colsum_distribution - end - - # Find partition elements with strictly positive measure. - δ = tolerance/size(to.transfermatrix, 1) - inds_nonzero = findall(distribution .> δ) - - # Extract the elements of the invariant measure corresponding to these indices - return InvariantMeasure(to, Probabilities(distribution)) -end - -function invariantmeasure(x::AbstractDataset, ϵ::RectangularBinning) - to = transferoperator(x, ϵ) - invariantmeasure(to) -end - -function probabilities(x::AbstractDataset, est::TransferOperator{RectangularBinning}) - to = transferoperator(x, est.ϵ) - iv = invariantmeasure(to) - - return iv.ρ -end - -""" - transfermatrix(iv::InvariantMeasure) → (M::AbstractArray{<:Real, 2}, bins::Vector{<:SVector}) - -Return the transfer matrix/operator and corresponding bins. Here, `bins[i]` corresponds -to the i-th row/column of the transfer matrix. Thus, the entry `M[i, j]` is the -probability of jumping from the state defined by `bins[i]` to the state defined by -`bins[j]`. - -See also: [`TransferOperator`](@ref). -""" -function transfermatrix(iv::InvariantMeasure) - return iv.to.transfermatrix, iv.to.bins -end diff --git a/src/binning_based/rectangular/rectangular_estimators.jl b/src/binning_based/rectangular/rectangular_estimators.jl deleted file mode 100644 index 9c543f2bd..000000000 --- a/src/binning_based/rectangular/rectangular_estimators.jl +++ /dev/null @@ -1,9 +0,0 @@ -abstract type BinningProbabilitiesEstimator <: ProbabilitiesEstimator end - -include("binning_schemes.jl") - -include("count_box_visits.jl") -include("histogram_estimation.jl") - -include("VisitationFrequency.jl") -include("TransferOperator.jl") \ No newline at end of file diff --git a/src/binning_based/rectangular_estimators.jl b/src/binning_based/rectangular_estimators.jl new file mode 100644 index 000000000..40849c964 --- /dev/null +++ b/src/binning_based/rectangular_estimators.jl @@ -0,0 +1,6 @@ +abstract type BinningProbabilitiesEstimator <: ProbabilitiesEstimator end + +include("GroupSlices.jl") +include("binning_schemes.jl") +include("visitation_frequency/VisitationFrequency.jl") +include("transferoperator/transfer_operator.jl") \ No newline at end of file diff --git a/src/binning_based/transferoperator.jl b/src/binning_based/transferoperator.jl new file mode 100644 index 000000000..8c56a3ebd --- /dev/null +++ b/src/binning_based/transferoperator.jl @@ -0,0 +1,131 @@ +using DelayEmbeddings, SparseArrays + +export + TransferOperator, # the probabilities estimator + InvariantMeasure, + transfermatrix, + invariantmeasure + +""" + TransferOperator(ϵ::RectangularBinning) <: BinningProbabilitiesEstimator + +A probability estimator based on binning data into rectangular boxes dictated by +the binning scheme `ϵ`, then approxmating the transfer (Perron-Frobenius) operator +over the bins, the taking the invariant measure associated with that transfer operator +as the bin probabilities. Assumes that the input data are sequential. + +This implementation follows the grid estimator approach in Diego et al. (2019)[^Diego2019]. + +## Description + +The transfer operator ``P^{N}``is computed as an `N`-by-`N` matrix of transition +probabilities between the states defined by the partition elements, where `N` is the +number of boxes in the partition that is visited by the orbit/points. + +If ``\\{x_t^{(D)} \\}_{n=1}^L`` are the ``L`` different ``D``-dimensional points over +which the transfer operator is approximated, ``\\{ C_{k=1}^N \\}`` are the ``N`` different +partition elements (as dictated by `ϵ`) that gets visited by the points, and + ``\\phi(x_t) = x_{t+1}``, then + +```math +P_{ij} = \\dfrac +{\\#\\{ x_n | \\phi(x_n) \\in C_j \\cap x_n \\in C_i \\}} +{\\#\\{ x_m | x_m \\in C_i \\}}, +``` + +where ``\\#`` denotes the cardinal. The element ``P_{ij}`` thus indicates how many points +that are initially in box ``C_i`` end up in box ``C_j`` when the points in ``C_i`` are +projected one step forward in time. Thus, the row ``P_{ik}^N`` where +``k \\in \\{1, 2, \\ldots, N \\}`` gives the probability +of jumping from the state defined by box ``C_i`` to any of the other ``N`` states. It +follows that ``\\sum_{k=1}^{N} P_{ik} = 1`` for all ``i``. Thus, ``P^N`` is a row/right +stochastic matrix. + +### Invariant measure estimation from transfer operator + +The left invariant distribution ``\\mathbf{\\rho}^N`` is a row vector, where +``\\mathbf{\\rho}^N P^{N} = \\mathbf{\\rho}^N``. Hence, ``\\mathbf{\\rho}^N`` is a row +eigenvector of the transfer matrix ``P^{N}`` associated with eigenvalue 1. The distribution +``\\mathbf{\\rho}^N`` approximates the invariant density of the system subject to the +partition `ϵ`, and can be taken as a probability distribution over the partition elements. + +In practice, the invariant measure ``\\mathbf{\\rho}^N`` is computed using +[`invariantmeasure`](@ref), which also approximates the transfer matrix. The invariant distribution +is initialized as a length-`N` random distribution which is then applied to ``P^{N}``. +The resulting length-`N` distribution is then applied to ``P^{N}`` again. This process +repeats until the difference between the distributions over consecutive iterations is +below some threshold. + +## Probability and entropy estimation + +- `probabilities(x::AbstractDataset, est::TransferOperator{RectangularBinning})` estimates + probabilities for the bins defined by the provided binning (`est.ϵ`) +- `genentropy(x::AbstractDataset, est::TransferOperator{RectangularBinning})` does the same, + but computes generalized entropy using the probabilities. + + +See also: [`RectangularBinning`](@ref), [`invariantmeasure`](@ref). + +[^Diego2019]: Diego, D., Haaga, K. A., & Hannisdal, B. (2019). Transfer entropy computation using the Perron-Frobenius operator. Physical Review E, 99(4), 042212. +""" +struct TransferOperator{R} <: BinningProbabilitiesEstimator + ϵ::R + + function TransferOperator(ϵ::R) where R #<: RectangularBinning + new{R}(ϵ) + end +end +struct TransferOperatorGenerator{E <: TransferOperator, X, A} + method::E # estimator with its input parameters + pts::X # the phase space / reconstruted state space points + init::A # pre-initialized things that speed up estimation process +end + + +""" + transopergenerator(pts, method::TransferOperator) → to::TransferOperatorGenerator + +Initialize a generator that creates transfer operators on demand, based on the given `method`. +This is efficient, because some things can be initialized and reused. + +To approximate a transfer operator, call `to` as a function with the relevant arguments. + +```julia +to = transopergenerator(x, TransferOperator(RectangularBinning(5))) +for i in 1:1000 + s = to() + # do stuff with s and or x + result[i] = stuff +end +``` +""" +function transopergenerator end + +function transferoperator end + +function invariantmeasure end + +function transfermatrix end + +function transferoperator(pts, method::TransferOperator) + to = transopergenerator(pts, method) + to() +end + +""" + InvariantMeasure(to, ρ) + +Minimal return struct for [`invariantmeasure`](@ref) that contains the estimated invariant +measure `ρ`, as well as the transfer operator `to` from which it is computed (including +bin information). + +See also: [`invariantmeasure`](@ref). +""" +struct InvariantMeasure{T} + to::T + ρ::Probabilities + + function InvariantMeasure(to::T, ρ) where T + new{T}(to, ρ) + end +end diff --git a/src/binning_based/transferoperator/common.jl b/src/binning_based/transferoperator/common.jl new file mode 100644 index 000000000..73c753374 --- /dev/null +++ b/src/binning_based/transferoperator/common.jl @@ -0,0 +1,324 @@ +using DelayEmbeddings +using Entropies +import LinearAlgebra: norm + +export + transferoperator, + TransferOperator, + invariantmeasure, + InvariantMeasure, + transfermatrix + +""" + TransferOperator(ϵ::Union{RectangularBinning, SimplexPoint, SimplexExact}) <: BinningProbabilitiesEstimator + +A probability estimator based on binning data into rectangular boxes +(when `ϵ` is a [`RectangularBinning`](@ref)) or simplices (when `ϵ` is a +[`SimplexExact`](@ref) or `ϵ` is a [`SimplexPoint`](@ref)). + +Input data are assumed to be sequential. + +The transfer (Perron-Frobenius) operator is approximated over the bins, using different +methods depending on the type of partition chosen. The invariant measure associated with +the approximated transfer operator is taken the bin probabilities. + +This implementation follows the grid estimator approach in Diego et al. (2019)[^Diego2019]. + +# Description + +The transfer operator ``P^{N}``is computed as an `N`-by-`N` matrix of transition +probabilities between the states defined by the partition elements, where `N` is the +number of boxes in the partition that is visited by the orbit/points. + +If ``\\{x_t^{(D)} \\}_{n=1}^L`` are the ``L`` different ``D``-dimensional points over +which the transfer operator is approximated, ``\\{ C_{k=1}^N \\}`` are the ``N`` different +partition elements (as dictated by `ϵ`) that gets visited by the points, and + ``\\phi(x_t) = x_{t+1}``, then + +```math +P_{ij} = \\dfrac +{\\#\\{ x_n | \\phi(x_n) \\in C_j \\cap x_n \\in C_i \\}} +{\\#\\{ x_m | x_m \\in C_i \\}}, +``` + +where ``\\#`` denotes the cardinal. The element ``P_{ij}`` thus indicates how many points +that are initially in box ``C_i`` end up in box ``C_j`` when the points in ``C_i`` are +projected one step forward in time. Thus, the row ``P_{ik}^N`` where +``k \\in \\{1, 2, \\ldots, N \\}`` gives the probability +of jumping from the state defined by box ``C_i`` to any of the other ``N`` states. It +follows that ``\\sum_{k=1}^{N} P_{ik} = 1`` for all ``i``. Thus, ``P^N`` is a row/right +stochastic matrix. + +### Invariant measure estimation from transfer operator + +The left invariant distribution ``\\mathbf{\\rho}^N`` is a row vector, where +``\\mathbf{\\rho}^N P^{N} = \\mathbf{\\rho}^N``. Hence, ``\\mathbf{\\rho}^N`` is a row +eigenvector of the transfer matrix ``P^{N}`` associated with eigenvalue 1. The distribution +``\\mathbf{\\rho}^N`` approximates the invariant density of the system subject to the +partition `ϵ`, and can be taken as a probability distribution over the partition elements. + +In practice, the invariant measure ``\\mathbf{\\rho}^N`` is computed using +[`invariantmeasure`](@ref), which also approximates the transfer matrix. The invariant distribution +is initialized as a length-`N` random distribution which is then applied to ``P^{N}``. +The resulting length-`N` distribution is then applied to ``P^{N}`` again. This process +repeats until the difference between the distributions over consecutive iterations is +below some threshold. + +## Probability and entropy estimation + +- `probabilities(x::AbstractDataset, est::TransferOperator{RectangularBinning})` estimates + probabilities for the bins defined by the provided binning (`est.ϵ`) +- `genentropy(x::AbstractDataset, est::TransferOperator{RectangularBinning})` does the same, + but computes generalized entropy using the probabilities. + + +See also: [`RectangularBinning`](@ref), [`SimplexPoint`](@ref), [`SimplexExact`](@ref), +[`invariantmeasure`](@ref), + +[^Diego2019]: Diego, D., Haaga, K. A., & Hannisdal, B. (2019). Transfer entropy computation using the Perron-Frobenius operator. Physical Review E, 99(4), 042212. +""" +struct TransferOperator{R} <: BinningProbabilitiesEstimator + ϵ::R + + function TransferOperator(ϵ::R) where R #<: RectangularBinning + new{R}(ϵ) + end +end + +struct TransferOperatorGenerator{E <: TransferOperator, X, A} + method::E # estimator with its input parameters + pts::X # the phase space / reconstruted state space points + init::A # pre-initialized things that speed up estimation process +end + + +function Base.show(io::IO, DT::TransferOperatorGenerator{E, X, A}) where {E <: TransferOperator{R}, X, A} where R + summary = "Transfer approximation from $R estimator" + println(io, summary) +end + + +function transferoperator end + +""" + transopergenerator(pts, method::TransferOperator) → to::TransferOperatorGenerator + +Initialize a generator that creates transfer operators on demand, based on the given `method`. +This is efficient, because some things can be initialized and reused. + +To approximate a transfer operator, call `to` as a function with the relevant arguments. + +```julia +to = transopergenerator(x, TransferOperator(RectangularBinning(5))) +for i in 1:1000 + s = to() + # do stuff with s and or x + result[i] = stuff +end +``` +""" +function transopergenerator end + +function transferoperator(pts, method::TransferOperator) + to = transopergenerator(pts, method) + to() +end + +""" + InvariantMeasure(to, ρ) + +Minimal return struct for [`invariantmeasure`](@ref) that contains the estimated invariant +measure `ρ`, as well as the transfer operator `to` from which it is computed (including +bin information). + +See also: [`invariantmeasure`](@ref). +""" +struct InvariantMeasure{T} + to::T + ρ::Probabilities + + function InvariantMeasure(to::T, ρ) where T + new{T}(to, ρ) + end +end + +""" + TransferOperatorApproximation(generator, transfermatrix, params) + +A return struct containing the `generator` of the transfer operator approximation (an +instance of `TransferOperator{R} where R`; if `R` is a `RectangularBinning`, then a +rectangular binning was used, whereas if `R` is `SimplexPoint` or `SimplexExact`, then +a triangulated binning as used). + +The `transfermatrix` is an `N`-by-`N` matrix approximation to the transfer operator, +subject to the partition given by `generator` (details about the partition are available +in the `generator.init` field), computed over some set of sequentially ordered points. + +# Parameters + +The parameters of the approximation is given to `params`. + +When `R <: RectangularBinning`, parameters stored are as follows. + +- The `mini` and `edgelengths` parameters provide the minima and box edge lengths along each + coordinate axis, as determined by applying `ϵ` to the points. +- The coordinates of the (leftmost, if axis is ordered low-high) box corners are given in + `bins`. Only bins actually visited by the points are considered, and `bins` give the + coordinates of these bins. The element `bins[i]` correspond to the `i`-th state of the + system, which corresponds to the `i`-th column/row of `transfermatrix`. +- `sort_idxs` contains the indices that would sort the input points. +- `visitors` is a vector of vectors, where `visitors[i]` contains the indices of the + (sorted) points that visits `bins[i]`. + +When `R <: SimplexPoint`, parameters stored are the simplex subsampling parameters. +When `R <: SimplexExact`, parameters stored are the simplex intersection tolerances. + +See also: [`RectangularBinning`](@ref), [`SimplexPoint`](@ref), [`SimplexExact`](@ref). +""" +struct TransferOperatorApproximation{G<:TransferOperator{<:BinningProbabilitiesEstimator}, T<:Real} + generator::TransferOperatorGenerator{G} + transfermatrix::AbstractArray{T, 2} + params +end + +function Base.show(io::IO, DT::TransferOperatorApproximation{G, T}) where {G <: TransferOperator{R}, T} where R + summary = "TransferOperatorApproximation{$G, $T}" + println(io, summary) +end + + +""" + transfermatrix(iv::InvariantMeasure) → (M::AbstractArray{<:Real, 2}, bins::Vector{<:SVector}) + +Return the transfer matrix/operator and corresponding bins. Here, `bins[i]` corresponds +to the i-th row/column of the transfer matrix. Thus, the entry `M[i, j]` is the +probability of jumping from the state defined by `bins[i]` to the state defined by +`bins[j]`. + +See also: [`TransferOperator`](@ref). +""" +function transfermatrix(iv::InvariantMeasure) + return iv.to.transfermatrix, iv.to.bins +end + + + +import LinearAlgebra: norm +""" + invariantmeasure(x::AbstractDataset, to::TransferOperator{R}) + where {R <: Union{RectangularBinning, SimplexExact, SimplexPoint}} + → iv::InvariantMeasure + +Estimate an invariant measure over the points in `x` based on binning the data into +rectangular boxes dictated by the binning scheme `ϵ`, then approximate the transfer +(Perron-Frobenius) operator over the bins. From the approximation to the transfer operator, +compute an invariant distribution over the bins. Assumes that the input data are sequential. + +Details on the estimation procedure is found the [`TransferOperator`](@ref) docstring. + +## Example + +```julia +using DynamicalSystems, Plots, Entropies +D = 4 +ds = Systems.lorenz96(D; F = 32.0) +N, dt = 20000, 0.1 +orbit = trajectory(ds, N*dt; dt = dt, Ttr = 10.0) + +# Estimate the invariant measure over some rectangular coarse graining of the orbit. +iv = invariantmeasure(orbit, RectangularBinning(15)) + +# Get the probabilities and bins +invariantmeasure(iv) +``` + +## Probabilities and bin information + + invariantmeasure(iv::InvariantMeasure) → (ρ::Probabilities, bins::Vector{<:SVector}) + +From a pre-computed invariant measure, return the probabilities and associated bins. +The element `ρ[i]` is the probability of visitation to the box `bins[i]`. Analogous to +[`binhist`](@ref). + + +!!! hint "Transfer operator approach vs. naive histogram approach" + + Why bother with the transfer operator instead of using regular histograms to obtain + probabilities? + + In fact, the naive histogram approach and the + transfer operator approach are equivalent in the limit of long enough time series + (as ``n \\to \\intfy``), which is guaranteed by the ergodic theorem. There is a crucial + difference, however: + + The naive histogram approach only gives the long-term probabilities that + orbits visit a certain region of the state space. The transfer operator encodes that + information too, but comes with the added benefit of knowing the *transition + probabilities* between states (see [`transfermatrix`](@ref)). + +See also: [`InvariantMeasure`](@ref). +""" +function invariantmeasure(to::TransferOperatorApproximation; + N::Int = 200, tolerance::Float64 = 1e-8, delta::Float64 = 1e-8) + + TO = to.transfermatrix + #= + # Start with a random distribution `Ρ` (big rho). Normalise it so that it + # sums to 1 and forms a true probability distribution over the partition elements. + =# + Ρ = rand(Float64, 1, size(to.transfermatrix, 1)) + Ρ = Ρ ./ sum(Ρ, dims = 2) + + #= + # Start estimating the invariant distribution. We could either do this by + # finding the left-eigenvector of M, or by repeated application of M on Ρ + # until the distribution converges. Here, we use the latter approach, + # meaning that we iterate until Ρ doesn't change substantially between + # iterations. + =# + distribution = Ρ * to.transfermatrix + + distance = norm(distribution - Ρ) / norm(Ρ) + + check = floor(Int, 1 / delta) + check_pts = floor.(Int, transpose(collect(1:N)) ./ check) .* transpose(collect(1:N)) + check_pts = check_pts[check_pts .> 0] + num_checkpts = size(check_pts, 1) + check_pts_counter = 1 + + counter = 1 + while counter <= N && distance >= tolerance + counter += 1 + Ρ = distribution + + # Apply the Markov matrix to the current state of the distribution + distribution = Ρ * to.transfermatrix + + if (check_pts_counter <= num_checkpts && + counter == check_pts[check_pts_counter]) + + check_pts_counter += 1 + colsum_distribution = sum(distribution, dims = 2)[1] + if abs(colsum_distribution - 1) > delta + distribution = distribution ./ colsum_distribution + end + end + + distance = norm(distribution - Ρ) / norm(Ρ) + end + distribution = dropdims(distribution, dims = 1) + + # Do the last normalisation and check + colsum_distribution = sum(distribution) + + if abs(colsum_distribution - 1) > delta + distribution = distribution ./ colsum_distribution + end + + # Find partition elements with strictly positive measure. + δ = tolerance/size(to.transfermatrix, 1) + inds_nonzero = findall(distribution .> δ) + + # Extract the elements of the invariant measure corresponding to these indices + return InvariantMeasure(to, Probabilities(distribution)) +end diff --git a/src/binning_based/transferoperator/rectangular/rectangular.jl b/src/binning_based/transferoperator/rectangular/rectangular.jl new file mode 100644 index 000000000..8d1a6df61 --- /dev/null +++ b/src/binning_based/transferoperator/rectangular/rectangular.jl @@ -0,0 +1,208 @@ +using DelayEmbeddings, SparseArrays + +# If x is not sorted, we need to look at all pairwise comparisons +function inds_in_terms_of_unique(x) + U = unique(x) + N = length(x) + Nu = length(U) + inds = zeros(Int, N) + + for j = 1:N + xⱼ = view(x, j) + for i = 1:Nu + # using views doesn't allocate + @inbounds if xⱼ == view(U, i) + inds[j] = i + end + end + end + + return inds +end + +# Taking advantage of the fact that x is sorted reduces runtime by 1.5 orders of magnitude +# for datasets of >100 000+ points +function inds_in_terms_of_unique_sorted(x) # assumes sorted + @assert issorted(x) + U = unique(x) + N, Nu = length(x), length(U) + prev = view(x, 1) + inds = zeros(Int, N) + uidx = 1 + @inbounds for j = 1:N + xⱼ = view(x, j) + # if the current value has changed, then we know that the corresponding index + # for the unique point must be incremented by 1 + if xⱼ != prev + prev = xⱼ + uidx += 1 + end + inds[j] = uidx + end + + return inds +end + +function inds_in_terms_of_unique(x, sorted::Bool) + if sorted + return inds_in_terms_of_unique_sorted(x) + else + return inds_in_terms_of_unique(x) + end +end + +inds_in_terms_of_unique(x::AbstractDataset) = inds_in_terms_of_unique(x.data) + +""" + transferoperator(pts::AbstractDataset{D, T}, ϵ::RectangularBinning) → TransferOperatorApproximationRectangular + +Estimate the transfer operator given a set of sequentially ordered points subject to a +rectangular partition given by `ϵ`. + +## Example + +```julia +using DynamicalSystems, Plots, Entropy +D = 4 +ds = Systems.lorenz96(D; F = 32.0) +N, dt = 20000, 0.1 +orbit = trajectory(ds, N*dt; dt = dt, Ttr = 10.0) + +# Estimate transfer operator over some coarse graining of the orbit. +transferoperator(orbit, RectangularBinning(10)) +``` + +See also: [`RectangularBinning`](@ref). +""" +function transferoperator(pts::AbstractDataset{D, T}, to::TransferOperator{RectangularBinning}; + boundary_condition = :circular) where {D, T<:Real} + + L = length(pts) + mini, edgelengths = Entropies.minima_edgelengths(pts, ϵ) + + # The L points visits a total of L bins, which are the following bins: + visited_bins = Entropies.encode_as_bin(pts, mini, edgelengths) + sort_idxs = sortperm(visited_bins) + + # TODO: fix re-indexing after sorting. Sorting is much faster, so we want to do so. + #sort!(visited_bins) + + # There are N=length(unique(visited_bins)) unique bins. + # Which of the unqiue bins does each of the L points visit? + visits_whichbin = inds_in_terms_of_unique(visited_bins, false) # set to true when sorting is fixed + + # `visitors` lists the indices of the points visiting each of the N unique bins. + slices = GroupSlices.groupslices(visited_bins) + visitors = GroupSlices.groupinds(slices) + + # first_visited_by == [x[1] for x in visitors] + first_visited_by = GroupSlices.firstinds(slices) + L = length(first_visited_by) + + I = Int32[] + J = Int32[] + P = Float64[] + + # Preallocate target index for the case where there is only + # one point of the orbit visiting a bin. + target_bin_j::Int = 0 + n_visitsᵢ::Int = 0 + + if boundary_condition == :circular + #warn("Using circular boundary condition") + append!(visits_whichbin, [1]) + elseif boundary_condition == :random + #warn("Using random circular boundary condition") + append!(visits_whichbin, [rand(1:length(visits_whichbin))]) + else + error("Boundary condition $(boundary_condition) not implemented") + end + + # Loop over the visited bins bᵢ + for i in 1:L + # How many times is this bin visited? + n_visitsᵢ = length(visitors[i]) + + # If both conditions below are true, then there is just one + # point visiting the i-th bin. If there is only one visiting point and + # it happens to be the last, we skip it, because we don't know its + # image. + if n_visitsᵢ == 1 && !(i == visits_whichbin[end]) + # To which bin does the single point visiting bᵢ jump if we + # shift it one time step ahead along its orbit? + target_bin_j = visits_whichbin[visitors[i][1] + 1][1] + + # We now know that exactly one point (the i-th) does the + # transition from i to the target j. + push!(I, i) + push!(J, target_bin_j) + push!(P, 1.0) + end + + # If more than one point of the orbit visits the i-th bin, we + # identify the visiting points and track which bins bⱼ they end up + # in after the forward linear map of the points. + if n_visitsᵢ > 1 + timeindices_visiting_pts = visitors[i] + + # TODO: Introduce circular boundary condition. Simply excluding + # might lead to a cascade of loosing points. + + # If bᵢ is the bin visited by the last point in the orbit, then + # the last entry of `visiting_pts` will be the time index of the + # last point of the orbit. In the next time step, that point will + # have nowhere to go along its orbit (precisely because it is the + # last data point). Thus, we exclude it. + if i == visits_whichbin[end] + #warn("Removing last point") + n_visitsᵢ = length(timeindices_visiting_pts) - 1 + timeindices_visiting_pts = timeindices_visiting_pts[1:(end - 1)] + end + + # To which boxes do each of the visitors to bᵢ jump in the next + # time step? + target_bins = visits_whichbin[timeindices_visiting_pts .+ 1] + unique_target_bins = unique(target_bins) + + # Count how many points jump from the i-th bin to each of + # the unique target bins, and use that to calculate the transition + # probability from bᵢ to bⱼ. + #for j in 1:length(unique_target_bins) + for (j, bᵤ) in enumerate(unique(target_bins)) + n_transitions_i_to_j = sum(target_bins .== bᵤ) + + push!(I, i) + push!(J, bᵤ) + push!(P, n_transitions_i_to_j / n_visitsᵢ) + end + end + end + + # Transfer operator is just the normalized transition probabilities between the boxes. + M = sparse(I, J, P) + + # Compute the coordinates of the visited bins. bins[i] corresponds to the i-th + # row/column of the transfer operator + unique!(visited_bins) + bins = [β .* edgelengths .+ mini for β in visited_bins] + params = (ϵ = ϵ, mini = mini, edgelengths = edgelengths, bins = bins, + sort_idxs = sort_idxs, visitors) + + TransferOperatorApproximation(to, M, params) +end + +function invariantmeasure(iv::InvariantMeasure) + return iv.ρ, iv.to.bins +end + +function invariantmeasure(x::AbstractDataset, ϵ::TransferOperator{<:RectangularBinning}) + to = transferoperator(x, ϵ) + invariantmeasure(to) +end + +function probabilities(x::AbstractDataset, est::TransferOperator{<:RectangularBinning}) + to = transferoperator(x, est.ϵ) + iv = invariantmeasure(to) + + return iv.ρ +end diff --git a/src/binning_based/transferoperator/transfer_operator.jl b/src/binning_based/transferoperator/transfer_operator.jl new file mode 100644 index 000000000..56420e74c --- /dev/null +++ b/src/binning_based/transferoperator/transfer_operator.jl @@ -0,0 +1,3 @@ +include("utils.jl") +include("common.jl") +include("rectangular/rectangular.jl") diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/AbstractDelaunayTriangulation-checkpoint.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/AbstractDelaunayTriangulation-checkpoint.jl new file mode 100644 index 000000000..c99f41bf0 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/AbstractDelaunayTriangulation-checkpoint.jl @@ -0,0 +1,37 @@ +#################################### +# Abstract type +#################################### + +abstract type AbstractDelaunayTriangulation end + +ADT = AbstractDelaunayTriangulation +Base.size(DT::ADT) = (length(DT.indices[1]) - 1, length(DT.indices)) +Base.length(DT::ADT) = length(DT.indices) +dimension(DT::ADT) = length(DT.indices[1]) - 1 +nsimplices(DT::ADT) = length(DT.indices) + +# Indexing +Base.getindex(DT::ADT, i) = DT.indices[i] +Base.getindex(DT::ADT, i::Colon, j::Colon) = hcat(DT.indices...,) +Base.getindex(DT::ADT, i::Colon, j) = hcat(DT[j]...,) +Base.getindex(DT::ADT, i::Colon, j::Int) = hcat(DT[j]) + +Base.firstindex(DT::ADT) = 1 +Base.lastindex(DT::ADT) = length(DT) + +Base.eachindex(s::ADT) = Base.OneTo(length(s)) +Base.iterate(s::ADT, state = 1) = iterate(s.indices, state) + +function summarise(DT::AbstractDelaunayTriangulation) + _type = typeof(DT) + n_simplices = nsimplices(DT) + D = dimension(DT) + summary = "$D-dimensional $_type with $n_simplices simplices" +end + +Base.show(io::IO, DT::AbstractDelaunayTriangulation) = println(io, summarise(DT)) + + +export +dimension, +nsimplices diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulation-checkpoint.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulation-checkpoint.jl new file mode 100644 index 000000000..319f96fc5 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulation-checkpoint.jl @@ -0,0 +1,20 @@ +export DelaunayTriangulation + +import StaticArrays: SVector, MVector + +include("addnoise.jl") + +struct DelaunayTriangulation <: AbstractDelaunayTriangulation + indices::Vector{Vector{Int32}} + + function DelaunayTriangulation(vertices::Vector{T}; joggle = 1e-8) where T <: Union{AbstractVector, SVector, MVector} + # Slightly joggle points to avoid problems with QHull + if joggle > 0 + addnoise!(vertices; joggle_factor = joggle) + end + simplex_inds = Simplices.Delaunay.delaunay(vertices) + new(simplex_inds) + end + + DelaunayTriangulation(vertices::Dataset; joggle = 1e-8) = DelaunayTriangulation(vertices.data) +end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulations-checkpoint.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulations-checkpoint.jl new file mode 100644 index 000000000..4a489cba6 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulations-checkpoint.jl @@ -0,0 +1,13 @@ +module DelaunayTriangulations + using Simplices + include("AbstractDelaunayTriangulation.jl") + include("DelaunayTriangulation.jl") + #include("plot_recipes/plot_recipes.jl") +end + +""" + DelaunayTriangulations + +A module handling the creation of Delaunay triangulations from data. +""" +DelaunayTriangulations diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/addnoise-checkpoint.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/addnoise-checkpoint.jl new file mode 100644 index 000000000..8b944c4d6 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/addnoise-checkpoint.jl @@ -0,0 +1,105 @@ + +import Distributions: Uniform +import Statistics: std +import StaticArrays: MArray, SArray, MVector, SVector + +""" + addnoise!(pts, joggle_factor) + +Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`. +""" +function addnoise!(pts::AbstractArray{T, 2}; joggle_factor = 0.001) where T + + D = minimum(size(pts)) + npts = maximum(size(pts)) + + if size(pts, 1) > size(pts, 2) + dim = 1 + else + dim = 2 + end + + # Scale standard deviation along each axis by joggle factor + σ = joggle_factor .* std(pts, dims = dim) + + for i in 1:D + r = [rand(Uniform(-σ[i], σ[i])) for pt in 1:npts] + if dim == 1 + pts[:, i] .+= r + elseif dim == 2 + pts[i, :] .+= r + end + end +end + +""" + addnoise!(pts, joggle_factor) + +Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`. +""" +function addnoise!(pts::Vector{Vector{T}}; joggle_factor = 1e-8) where T + + D = length(pts[1]) + npts = length(pts) + + # Scale standard deviation along each axis by joggle factor + σ = joggle_factor + + for i in 1:D + r = rand(Uniform(-σ, σ)) + pts[i] .+= r + end +end + +import DelayEmbeddings: + Dataset + +function addnoise!(pts::Dataset; joggle_factor = 1e-8) + + D = size(pts, 2) + npts = length(pts) + + # Scale standard deviation along each axis by joggle factor + σ = joggle_factor + + Dataset([pts[i] .+ rand(Uniform(-σ, σ)) for i = 1:npts]) +end + + + + +""" + addnoise!(pts, joggle_factor) + +Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`. +""" +function addnoise!(pts::Vector{MVector{D, T}}; joggle_factor = 1e-8) where {D, T} + + npts = length(pts) + + # Scale standard deviation along each axis by joggle factor + σ = joggle_factor + + for i in 1:D + r = rand(Uniform(-σ, σ)) + pts[i] .+= r + end +end + + + + +""" + addnoise!(pts, joggle_factor) + +Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`. +""" +function addnoise!(pts::Vector{SVector{D, T}}; joggle_factor = 1e-8) where {D, T} + + npts = length(pts) + + # Scale standard deviation along each axis by joggle factor + σ = joggle_factor + + [pts[i] .+ rand(Uniform(-σ, σ)) for i = 1:npts] +end diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/AbstractDelaunayTriangulation.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/AbstractDelaunayTriangulation.jl new file mode 100644 index 000000000..c99f41bf0 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/delaunay_triangulations/AbstractDelaunayTriangulation.jl @@ -0,0 +1,37 @@ +#################################### +# Abstract type +#################################### + +abstract type AbstractDelaunayTriangulation end + +ADT = AbstractDelaunayTriangulation +Base.size(DT::ADT) = (length(DT.indices[1]) - 1, length(DT.indices)) +Base.length(DT::ADT) = length(DT.indices) +dimension(DT::ADT) = length(DT.indices[1]) - 1 +nsimplices(DT::ADT) = length(DT.indices) + +# Indexing +Base.getindex(DT::ADT, i) = DT.indices[i] +Base.getindex(DT::ADT, i::Colon, j::Colon) = hcat(DT.indices...,) +Base.getindex(DT::ADT, i::Colon, j) = hcat(DT[j]...,) +Base.getindex(DT::ADT, i::Colon, j::Int) = hcat(DT[j]) + +Base.firstindex(DT::ADT) = 1 +Base.lastindex(DT::ADT) = length(DT) + +Base.eachindex(s::ADT) = Base.OneTo(length(s)) +Base.iterate(s::ADT, state = 1) = iterate(s.indices, state) + +function summarise(DT::AbstractDelaunayTriangulation) + _type = typeof(DT) + n_simplices = nsimplices(DT) + D = dimension(DT) + summary = "$D-dimensional $_type with $n_simplices simplices" +end + +Base.show(io::IO, DT::AbstractDelaunayTriangulation) = println(io, summarise(DT)) + + +export +dimension, +nsimplices diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/DelaunayTriangulation.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/DelaunayTriangulation.jl new file mode 100644 index 000000000..319f96fc5 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/delaunay_triangulations/DelaunayTriangulation.jl @@ -0,0 +1,20 @@ +export DelaunayTriangulation + +import StaticArrays: SVector, MVector + +include("addnoise.jl") + +struct DelaunayTriangulation <: AbstractDelaunayTriangulation + indices::Vector{Vector{Int32}} + + function DelaunayTriangulation(vertices::Vector{T}; joggle = 1e-8) where T <: Union{AbstractVector, SVector, MVector} + # Slightly joggle points to avoid problems with QHull + if joggle > 0 + addnoise!(vertices; joggle_factor = joggle) + end + simplex_inds = Simplices.Delaunay.delaunay(vertices) + new(simplex_inds) + end + + DelaunayTriangulation(vertices::Dataset; joggle = 1e-8) = DelaunayTriangulation(vertices.data) +end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/DelaunayTriangulations.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/DelaunayTriangulations.jl new file mode 100644 index 000000000..4a489cba6 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/delaunay_triangulations/DelaunayTriangulations.jl @@ -0,0 +1,13 @@ +module DelaunayTriangulations + using Simplices + include("AbstractDelaunayTriangulation.jl") + include("DelaunayTriangulation.jl") + #include("plot_recipes/plot_recipes.jl") +end + +""" + DelaunayTriangulations + +A module handling the creation of Delaunay triangulations from data. +""" +DelaunayTriangulations diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/addnoise.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/addnoise.jl new file mode 100644 index 000000000..8b944c4d6 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/delaunay_triangulations/addnoise.jl @@ -0,0 +1,105 @@ + +import Distributions: Uniform +import Statistics: std +import StaticArrays: MArray, SArray, MVector, SVector + +""" + addnoise!(pts, joggle_factor) + +Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`. +""" +function addnoise!(pts::AbstractArray{T, 2}; joggle_factor = 0.001) where T + + D = minimum(size(pts)) + npts = maximum(size(pts)) + + if size(pts, 1) > size(pts, 2) + dim = 1 + else + dim = 2 + end + + # Scale standard deviation along each axis by joggle factor + σ = joggle_factor .* std(pts, dims = dim) + + for i in 1:D + r = [rand(Uniform(-σ[i], σ[i])) for pt in 1:npts] + if dim == 1 + pts[:, i] .+= r + elseif dim == 2 + pts[i, :] .+= r + end + end +end + +""" + addnoise!(pts, joggle_factor) + +Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`. +""" +function addnoise!(pts::Vector{Vector{T}}; joggle_factor = 1e-8) where T + + D = length(pts[1]) + npts = length(pts) + + # Scale standard deviation along each axis by joggle factor + σ = joggle_factor + + for i in 1:D + r = rand(Uniform(-σ, σ)) + pts[i] .+= r + end +end + +import DelayEmbeddings: + Dataset + +function addnoise!(pts::Dataset; joggle_factor = 1e-8) + + D = size(pts, 2) + npts = length(pts) + + # Scale standard deviation along each axis by joggle factor + σ = joggle_factor + + Dataset([pts[i] .+ rand(Uniform(-σ, σ)) for i = 1:npts]) +end + + + + +""" + addnoise!(pts, joggle_factor) + +Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`. +""" +function addnoise!(pts::Vector{MVector{D, T}}; joggle_factor = 1e-8) where {D, T} + + npts = length(pts) + + # Scale standard deviation along each axis by joggle factor + σ = joggle_factor + + for i in 1:D + r = rand(Uniform(-σ, σ)) + pts[i] .+= r + end +end + + + + +""" + addnoise!(pts, joggle_factor) + +Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`. +""" +function addnoise!(pts::Vector{SVector{D, T}}; joggle_factor = 1e-8) where {D, T} + + npts = length(pts) + + # Scale standard deviation along each axis by joggle factor + σ = joggle_factor + + [pts[i] .+ rand(Uniform(-σ, σ)) for i = 1:npts] +end diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/composite_types/.ipynb_checkpoints/DelaunayTriangulation-checkpoint.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/composite_types/.ipynb_checkpoints/DelaunayTriangulation-checkpoint.jl new file mode 100644 index 000000000..95cf2fecb --- /dev/null +++ b/src/binning_based/transferoperator/triangular/delaunay_triangulations/composite_types/.ipynb_checkpoints/DelaunayTriangulation-checkpoint.jl @@ -0,0 +1,21 @@ +export DelaunayTriangulation + +import StaticArrays: SVector, MVector +import .Simplices + +include("addnoise.jl") + +struct DelaunayTriangulation <: AbstractDelaunayTriangulation + indices::Vector{Vector{Int32}} + + function DelaunayTriangulation(vertices::Vector{T}; joggle = 1e-8) where T <: Union{AbstractVector, SVector, MVector} + # Slightly joggle points to avoid problems with QHull + if joggle > 0 + addnoise!(vertices; joggle_factor = joggle) + end + simplex_inds = Simplices.Delaunay.delaunay(vertices) + new(simplex_inds) + end + + DelaunayTriangulation(vertices::Dataset; joggle = 1e-8) = DelaunayTriangulation(vertices.data) +end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/plot_recipes/.ipynb_checkpoints/plot_recipes-checkpoint.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/plot_recipes/.ipynb_checkpoints/plot_recipes-checkpoint.jl new file mode 100644 index 000000000..8946af7c8 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/delaunay_triangulations/plot_recipes/.ipynb_checkpoints/plot_recipes-checkpoint.jl @@ -0,0 +1,60 @@ +using RecipesBase + +import StaticArrays: + SArray, + MArray + +import ..Simplices: + Simplex, + SSimplex, + MutableSSimplex, + connectvertices, + splitaxes + +################################## +# Plotting single triangulations +################################## +@recipe function f(pts::Vector{Vector{T}}, DT::DelaunayTriangulation) where {T} + + for i = 1:length(DT) + s = Simplex(pts[DT[i]]) + @series begin + seriestype := :path + splitaxes(connectvertices(s)) + end + end +end + +@recipe function f(pts::AbstractArray{T, 2}, DT::DelaunayTriangulation) where {T} + + for i = 1:length(DT) + s = Simplex(pts[:, DT[i]]) + @series begin + seriestype := :path + splitaxes(connectvertices(s)) + end + end +end + +@recipe function f(pts::Vector{MArray{Tuple{D},T,1,D}}, DT::DelaunayTriangulation) where {D, T} + + for vertexinds in DT + s = MutableSSimplex(pts[vertexinds]) + @series begin + seriestype := :path + splitaxes(connectvertices(s)) + end + end +end + + +@recipe function f(pts::Vector{SArray{Tuple{D},T,1,D}}, DT::DelaunayTriangulation) where {D, T} + + for vertexinds in DT + s = SSimplex(pts[vertexinds]) + @series begin + seriestype := :path + splitaxes(connectvertices(s)) + end + end +end diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/plot_recipes/plot_recipes.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/plot_recipes/plot_recipes.jl new file mode 100644 index 000000000..8946af7c8 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/delaunay_triangulations/plot_recipes/plot_recipes.jl @@ -0,0 +1,60 @@ +using RecipesBase + +import StaticArrays: + SArray, + MArray + +import ..Simplices: + Simplex, + SSimplex, + MutableSSimplex, + connectvertices, + splitaxes + +################################## +# Plotting single triangulations +################################## +@recipe function f(pts::Vector{Vector{T}}, DT::DelaunayTriangulation) where {T} + + for i = 1:length(DT) + s = Simplex(pts[DT[i]]) + @series begin + seriestype := :path + splitaxes(connectvertices(s)) + end + end +end + +@recipe function f(pts::AbstractArray{T, 2}, DT::DelaunayTriangulation) where {T} + + for i = 1:length(DT) + s = Simplex(pts[:, DT[i]]) + @series begin + seriestype := :path + splitaxes(connectvertices(s)) + end + end +end + +@recipe function f(pts::Vector{MArray{Tuple{D},T,1,D}}, DT::DelaunayTriangulation) where {D, T} + + for vertexinds in DT + s = MutableSSimplex(pts[vertexinds]) + @series begin + seriestype := :path + splitaxes(connectvertices(s)) + end + end +end + + +@recipe function f(pts::Vector{SArray{Tuple{D},T,1,D}}, DT::DelaunayTriangulation) where {D, T} + + for vertexinds in DT + s = SSimplex(pts[vertexinds]) + @series begin + seriestype := :path + splitaxes(connectvertices(s)) + end + end +end diff --git a/src/binning_based/transferoperator/triangular/exact/.ipynb_checkpoints/SimplexExact-checkpoint.jl b/src/binning_based/transferoperator/triangular/exact/.ipynb_checkpoints/SimplexExact-checkpoint.jl new file mode 100644 index 000000000..ec913f617 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/exact/.ipynb_checkpoints/SimplexExact-checkpoint.jl @@ -0,0 +1,154 @@ +export SimplexExact, invariantmeasure +import Entropies.invariantmeasure + + + +""" + SimplexExact() <: TriangulationBasedTransferOperator + +A transfer operator estimator using a triangulation partition and exact +simplex intersections[^Diego2019]. + +To use this estimator, the Simplices.jl +package must be brought into scope by doing `using Simplices` after +`using PerronFrobenius`. *Note: due to computing the exact simplex intersection, +this estimator is slow compared to `SimplexApprox`.* + +[^Diego2019]: Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212. +""" +struct SimplexExact <: BinningProbabilitiesEstimator + bc::String + + function SimplexExact(bc::String = "circular") + isboundarycondition(bc, "triangulation") || error("Boundary condition '$bc' not valid.") + new(bc) + end +end +Base.show(io::IO, se::SimplexExact) = print(io, "SimplexExact{$(se.bc)}") + + +""" Generate a transopergenerator for an exact simplex estimator.""" +function transopergenerator(pts, method::TransferOperator{<:SimplexExact}) + # modified points, where the image of each point is guaranteed to lie within the convex hull of the previous points + invariant_pts = invariantize(pts) + + # triangulation of the invariant points. The last point is excluded, so that + # the last vertex also can be mapped forward one step in time. + triang = DelaunayTriangulation(invariant_pts[1:end-1]) + + init = (invariant_pts = invariant_pts, triang = triang) + + TransferOperatorGenerator(method, pts, init) +end + +function (tog::TransferOperatorGenerator{T})(; tol = 1e-8) where T <: TransferOperator{SimplexExact} + invariant_pts, triang = getfield.(Ref(tog.init), (:invariant_pts, :triang)) + + D = length(invariant_pts[1]) + N = length(triang) + ϵ = tol / N + + # Pre-allocate simplex and its image + image_simplex = MutableSimplex(zeros(Float64, D, D+1)) + simplex = MutableSimplex(zeros(Float64, D, D+1)) + M = zeros(Float64, N, N) + + for j in 1:N + for k = 1:(D + 1) + # Get the vertices of the image simplex + image_simplex[k] .= invariant_pts[triang[j][k] + 1] + end + + imvol = abs(orientation(image_simplex)) + + for i = 1:N + for k = 1:(D + 1) + # Get the vertices of the source simplex + simplex[k] .= invariant_pts[triang[i][k]] + end + + # Only compute the entry of the transfer matrix + # if simplices are of sufficient size. + vol = abs(orientation(simplex)) + + if vol * imvol > 0 && (vol/imvol) > ϵ + M[j, i] = intersect(simplex, image_simplex) / imvol + end + end + end + params = (; tol = tol) + TransferOperatorApproximation(tog, M, params) +end + +function transferoperator(pts, method::SimplexExact; tol::Real = 1e-8) + tog = transopergenerator(pts, method) + tog(tol = tol) +end + + +function invariantmeasure(to::TransferOperatorApproximation{G, T}; + N::Int = 200, + tolerance::Float64 = 1e-8, + delta::Float64 = 1e-8) where {G <: TransferOperator{R}, T} where R <: SimplexExact + + TO = to.transfermatrix + #= + # Start with a random distribution `Ρ` (big rho). Normalise it so that it + # sums to 1 and forms a true probability distribution over the partition elements. + =# + Ρ = rand(Float64, 1, size(to.transfermatrix, 1)) + Ρ = Ρ ./ sum(Ρ, dims = 2) + + #= + # Start estimating the invariant distribution. We could either do this by + # finding the left-eigenvector of M, or by repeated application of M on Ρ + # until the distribution converges. Here, we use the latter approach, + # meaning that we iterate until Ρ doesn't change substantially between + # iterations. + =# + distribution = Ρ * TO + + distance = norm(distribution - Ρ) / norm(Ρ) + + check = floor(Int, 1 / delta) + check_pts = floor.(Int, transpose(collect(1:N)) ./ check) .* transpose(collect(1:N)) + check_pts = check_pts[check_pts .> 0] + num_checkpts = size(check_pts, 1) + check_pts_counter = 1 + + counter = 1 + while counter <= N && distance >= tolerance + counter += 1 + Ρ = distribution + + # Apply the Markov matrix to the current state of the distribution + distribution = Ρ * to.transfermatrix + + if (check_pts_counter <= num_checkpts && + counter == check_pts[check_pts_counter]) + + check_pts_counter += 1 + colsum_distribution = sum(distribution, dims = 2)[1] + if abs(colsum_distribution - 1) > delta + distribution = distribution ./ colsum_distribution + end + end + + distance = norm(distribution - Ρ) / norm(Ρ) + end + distribution = dropdims(distribution, dims = 1) + + # Do the last normalisation and check + colsum_distribution = sum(distribution) + + if abs(colsum_distribution - 1) > delta + distribution = distribution ./ colsum_distribution + end + + # Find partition elements with strictly positive measure. + δ = tolerance/size(to.transfermatrix, 1) + inds_nonzero = findall(distribution .> δ) + + # Extract the elements of the invariant measure corresponding to these indices + return InvariantMeasure(to, Probabilities(distribution)) +end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl new file mode 100644 index 000000000..fb6c08147 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl @@ -0,0 +1,155 @@ +export SimplexExact, invariantmeasure +import Entropies.invariantmeasure + + + +""" + SimplexExact <: TriangulationBasedTransferOperator + +A transfer operator estimator using a triangulation partition and exact +simplex intersections[^Diego2019]. + +To use this estimator, the Simplices.jl package must be brought into scope by doing +`using Simplices` after running `using Entropies`. + +*Note: due to computing exact +simplex intersections, this estimator is slow compared to [`SimplexPoint`](@ref).* + +[^Diego2019]: Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212. +""" +struct SimplexExact <: BinningProbabilitiesEstimator + bc::String + + function SimplexExact(bc::String = "circular") + isboundarycondition(bc, "triangulation") || error("Boundary condition '$bc' not valid.") + new(bc) + end +end +Base.show(io::IO, se::SimplexExact) = print(io, "SimplexExact{$(se.bc)}") + + +""" Generate a transopergenerator for an exact simplex estimator.""" +function transopergenerator(pts, method::TransferOperator{<:SimplexExact}) + # modified points, where the image of each point is guaranteed to lie within the convex hull of the previous points + invariant_pts = invariantize(pts) + + # triangulation of the invariant points. The last point is excluded, so that + # the last vertex also can be mapped forward one step in time. + triang = DelaunayTriangulation(invariant_pts[1:end-1]) + + init = (invariant_pts = invariant_pts, triang = triang) + + TransferOperatorGenerator(method, pts, init) +end + +function (tog::TransferOperatorGenerator{T})(; tol = 1e-8) where T <: TransferOperator{SimplexExact} + invariant_pts, triang = getfield.(Ref(tog.init), (:invariant_pts, :triang)) + + D = length(invariant_pts[1]) + N = length(triang) + ϵ = tol / N + + # Pre-allocate simplex and its image + image_simplex = MutableSimplex(zeros(Float64, D, D+1)) + simplex = MutableSimplex(zeros(Float64, D, D+1)) + M = zeros(Float64, N, N) + + for j in 1:N + for k = 1:(D + 1) + # Get the vertices of the image simplex + image_simplex[k] .= invariant_pts[triang[j][k] + 1] + end + + imvol = abs(orientation(image_simplex)) + + for i = 1:N + for k = 1:(D + 1) + # Get the vertices of the source simplex + simplex[k] .= invariant_pts[triang[i][k]] + end + + # Only compute the entry of the transfer matrix + # if simplices are of sufficient size. + vol = abs(orientation(simplex)) + + if vol * imvol > 0 && (vol/imvol) > ϵ + M[j, i] = intersect(simplex, image_simplex) / imvol + end + end + end + params = (; tol = tol) + TransferOperatorApproximation(tog, M, params) +end + +function transferoperator(pts, method::SimplexExact; tol::Real = 1e-8) + tog = transopergenerator(pts, method) + tog(tol = tol) +end + + +function invariantmeasure(to::TransferOperatorApproximation{G, T}; + N::Int = 200, + tolerance::Float64 = 1e-8, + delta::Float64 = 1e-8) where {G <: TransferOperator{R}, T} where R <: SimplexExact + + TO = to.transfermatrix + #= + # Start with a random distribution `Ρ` (big rho). Normalise it so that it + # sums to 1 and forms a true probability distribution over the partition elements. + =# + Ρ = rand(Float64, 1, size(to.transfermatrix, 1)) + Ρ = Ρ ./ sum(Ρ, dims = 2) + + #= + # Start estimating the invariant distribution. We could either do this by + # finding the left-eigenvector of M, or by repeated application of M on Ρ + # until the distribution converges. Here, we use the latter approach, + # meaning that we iterate until Ρ doesn't change substantially between + # iterations. + =# + distribution = Ρ * TO + + distance = norm(distribution - Ρ) / norm(Ρ) + + check = floor(Int, 1 / delta) + check_pts = floor.(Int, transpose(collect(1:N)) ./ check) .* transpose(collect(1:N)) + check_pts = check_pts[check_pts .> 0] + num_checkpts = size(check_pts, 1) + check_pts_counter = 1 + + counter = 1 + while counter <= N && distance >= tolerance + counter += 1 + Ρ = distribution + + # Apply the Markov matrix to the current state of the distribution + distribution = Ρ * to.transfermatrix + + if (check_pts_counter <= num_checkpts && + counter == check_pts[check_pts_counter]) + + check_pts_counter += 1 + colsum_distribution = sum(distribution, dims = 2)[1] + if abs(colsum_distribution - 1) > delta + distribution = distribution ./ colsum_distribution + end + end + + distance = norm(distribution - Ρ) / norm(Ρ) + end + distribution = dropdims(distribution, dims = 1) + + # Do the last normalisation and check + colsum_distribution = sum(distribution) + + if abs(colsum_distribution - 1) > delta + distribution = distribution ./ colsum_distribution + end + + # Find partition elements with strictly positive measure. + δ = tolerance/size(to.transfermatrix, 1) + inds_nonzero = findall(distribution .> δ) + + # Extract the elements of the invariant measure corresponding to these indices + return InvariantMeasure(to, Probabilities(distribution)) +end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/SimplexPoint-checkpoint.jl b/src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/SimplexPoint-checkpoint.jl new file mode 100644 index 000000000..2ac7139a9 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/SimplexPoint-checkpoint.jl @@ -0,0 +1,131 @@ +export SimplexPoint, transferoperator, SimplexPoint +import StaticArrays: SizedVector, SizedMatrix, MMatrix + +include("helper_functions.jl") + + +""" + SimplexPoint() <: Entropies.BinningProbabilitiesEstimator + +A probability estimator based on the transfer operator, using a triangulated partition +and approximate simplex intersections to compute transition probabilities [^Diego2019]. + +To use this estimator, the Simplices.jl package must be brought into scope by doing `using Simplices` after +`using PerronFrobenius`. + +[^Diego2019]: Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212. +""" +struct SimplexPoint <: BinningProbabilitiesEstimator + bc::String + + function SimplexPoint(bc::String = "circular") + isboundarycondition(bc, "triangulation") || error("Boundary condition '$bc' not valid.") + new(bc) + end +end +Base.show(io::IO, se::SimplexPoint) = print(io, "SimplexPoint{$(se.bc)}") + + +""" Generate a transopergenerator for an exact simplex estimator.""" +function transopergenerator(pts, method::TransferOperator{<:SimplexPoint}) + # modified points, where the image of each point is guaranteed to lie within the convex hull of the previous points + invariant_pts = invariantize(pts) + + # triangulation of the invariant points. The last point is excluded, so that + # the last vertex also can be mapped forward one step in time. + triang = DelaunayTriangulation(invariant_pts[1:end-1]) + + init = (invariant_pts = invariant_pts, triang = triang) + + TransferOperatorGenerator(method, pts, init) +end + + + +""" +n: The number of points sampled inside each simplex. +randomsampling: Sample randomly (`random == true`), or using an even simplex splitting routine (`random == false`) +""" +function (tog::TransferOperatorGenerator{T})(; tol = 1e-8, randomsampling::Bool = false, + n::Int = 100) where T <: TransferOperator{SimplexPoint} + pts, inds = getfield.(Ref(tog.init), (:invariant_pts, :triang)) + D = length(pts[1]) + N = length(inds) + ϵ = tol / N + + # Prepare memory-efficient representations of the simplices + simplices = reshape_simplices(pts[1:(end- 1)], inds) + image_simplices = reshape_simplices(pts[2:end], inds) + + #= + Compute the convex expansions coefficients needed to generate points + within a simplex. Also, update number of points if a shape-preserving + simplex subdivision algorithm was used (in that case, we get a + *minimum* of `n_sample_pts` set of coefficients. + =# + convex_coeffs = subsample_coeffs(D, n, randomsampling) + n_coeffs::Int = size(convex_coeffs, 2) + convex_coeffs = [convex_coeffs[:, i] for i = 1:n_coeffs] + + # Pre-allocated arrays (SizedArrays, for efficiency) + L = (D+1)^2 + NVERTICES = D+1 + pt = SizedMatrix{1, D}(zeros(1, D)) + s_arr = SizedVector{L}(zeros(L)) + signs = SizedVector{NVERTICES}(zeros(NVERTICES)) + + # Re-arrange simplices so that look-up is a bit more efficient + simplex_arrs = Vector{Array{Float64, 2}}(undef, N) + imsimplex_arrs = Vector{Array{Float64, 2}}(undef, N) + for i in 1:N + simplex_arrs[i] = hcat(pts[inds[i]]...,) + imsimplex_arrs[i] = hcat(pts[inds[i] .+ 1]...,) + end + + # Pre-allocate some information to avoid memory allocation when + # computing the fractional overlaps. These are a temporary matrix + # which we use to compute determinants (to check if a point lies + # inside a simplex), and the indices needed to access its values + # at different stages of the checking stage (of whether a point + # lies inside a simplex). Pre-allocating this information and + # passing it on reduces runtime and memory allocation + # by half an order of magnitude. + temp_arr = rand(MMatrix{D + 1, D + 1}) + starts = [((D + 1)*i)+1 for i = 0:D] + stops = starts .+ D + + + # The Markov matrix + M = zeros(Float64, N, N) + + for i in 1:N + # Find the indices of the simplices who potentially intersect + # with the image of simplex #i + idxs = idxs_potentially_intersecting_simplices(pts, inds, i) + Sj = Vector{AbstractArray}(undef, length(idxs)) + get_simplices_at_inds!(Sj, idxs, simplices) + + # Generate points within the image of simplex #i and check + # what fraction of those points fall into the simplices. + # We approximate the intersecting volume between the + # image simplex and the individual simplices by those + # fractions. + @views is = imsimplex_arrs[i] + for k in 1:n_coeffs + pt = transpose(convex_coeffs[k]) * transpose(is) + innerloop_optim!(idxs, signs, s_arr, Sj, pt, D, M, i, temp_arr, starts, stops) + end + end + # Need to normalise, because all we have up until now is counts + # of how many points inside the image simplex falls into + # the simplices. + params = (tol = tol, randomsampling = randomsampling, n = n) + M_normalized = transpose(M) ./ n_coeffs + return TransferOperatorApproximation(tog, M_normalized, params) + +end + +function transferoperator(pts, method::TransferOperator{<:SimplexPoint}; tol = 1e-8, randomsampling::Bool = false, n::Int = 100) + tog = transopergenerator(pts, method) + tog(tol = tol, randomsampling = randomsampling, n = n) +end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/helper_functions-checkpoint.jl b/src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/helper_functions-checkpoint.jl new file mode 100644 index 000000000..b49cd0a68 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/helper_functions-checkpoint.jl @@ -0,0 +1,175 @@ +import LinearAlgebra: det +import StaticArrays: SizedMatrix +import .DelaunayTriangulations: DelaunayTriangulation + +""" Re-zeros the array `a`. """ +function rezero!(a) + @inbounds for i in eachindex(a) + a[i] = 0.0 + end +end + +""" Fill the elements of vector `from` into vector `into`. """ +function fill_into!(into, from) + @inbounds for i in eachindex(into) + into[i] = from[i] + end +end + +""" Fill vector `into` at indices `inds` with the elements at `from[inds]`. """ +function fill_at_inds!(into, from, inds) + @inbounds for i in 1:length(inds) + into[inds[i]] = from[i] + end +end + +function mdet_fm_in_x_view!(M, v, n, starts, stops) + @inbounds for i = 0:(n - 1) + M[starts[i+1]:stops[i+1]] = view(v, starts[i+1]:stops[i+1]) + end + @fastmath det(M) +end + + +function iscontained!(signs, s_arr, sx, point, dim, temp_arr, starts, stops) + # Redefine the temporary simplex. This is in-place, so we don't allocate + # memory. We could also have re-initialised `signs`, but since we're never + # comparing more than two consecutive signs, this is not necessary. + rezero!(s_arr) + rezero!(signs) + fill_into!(s_arr, sx) + + #Replace first vertex with the point + fill_at_inds!(s_arr, point, 1:dim) + + # Signed volume + signs[1] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops)) + + rezero!(s_arr) + fill_into!(s_arr, sx) #reset + + for κ = 2:dim # Check remaining signs and stop if sign changes + # Replace the ith vertex with the point we're cheking (leaving the + # 1 appended to Vi intact.) + idxs = ((dim + 1)*(κ - 1)+1):((dim + 1)*(κ - 1)+ 1 + dim - 1) + fill_at_inds!(s_arr, point, idxs) # ith change + + signs[κ] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops)) + + if !(signs[κ-1] == signs[κ]) + return false + end + + rezero!(s_arr) + fill_into!(s_arr, sx) + end + + # Last the last vertex with the point in question + idxs = ((dim + 1)*(dim)+1):((dim+1)^2-1) + fill_at_inds!(s_arr, point, idxs) + + signs[end] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops)) + + if !(signs[end-1] == signs[end]) + return false + else + return true + end +end + +function innerloop_optim!(inds::Vector{Int}, signs, s_arr, Sj, pt, dim::Int, M, i::Int, temp_arr, starts, stops) + for j in 1:length(inds) + if iscontained!(signs, s_arr, Sj[j], pt, dim, temp_arr, starts, stops) + M[inds[j], i] += 1.0 + end + end +end + +""" + reshape_simplices(pts::Vector{Vector{T}}, DT::DelaunayTriangulation) where T + +Creates alternative representations of the simplices that allow for efficient +(mostly non-allocating) checking if a point lies inside a simplex. + +This function adds some information needed when calculating the orientations +of the simplices *before* we start the computation (appends a row of ones on top of the +matrix representation of the simplex). Otherwise, the fallback is machinery in Simplices.jl, +which is much slower because the additional information must be added +`(n_simplices^2)*n_sample_pts` times (creating that many new matrices), which takes forever. + +Reshaping the simplices and appending that information beforehand avoids excessive +memory allocation. +""" +function reshape_simplices(pts, DT::DelaunayTriangulation) + n_simplices = length(DT.indices) + n_vertices = length(DT[1]) + dim = n_vertices - 1 + + # The [:, j, i] th entry of these two arrays holds the jth vertex of the + # ith simplex, but instead of having just `dim` vertices, we append a `1` + # to the end of the vectors. This allows for efficent (non-allocating) + # computation within the `contains_point_lessalloc!` function. If we instead + # would have appended the 1's inside that function, we would be performing + # memory-allocating operations, which are very expensive. Doing this instead + # gives orders of magnitude speed-ups for sufficiently large triangulations. + S1 = Array{Float64}(undef, dim + 1, dim + 1, n_simplices) + + # Collect simplices in the form of (dim+1)^2-length column vectors. This + # also helps with the + NDIM, N_SIMPLICES = (dim+1)^2, n_simplices + simplices = SizedMatrix{NDIM}{N_SIMPLICES}(zeros((NDIM, N_SIMPLICES))) + + @inbounds for i in 1:n_simplices + for j in 1:n_vertices + S1[:, j, i] = vcat(pts[DT[i][j]], 1.0) + end + + simplices[:, i] = reshape(S1[:, :, i], (dim+1)^2) + end + + return simplices +end + + +""" +idxs_potentially_intersecting_simplices(all_pts::Vector{Vector{T}}, + DT::DelaunayTriangulation, idx::Int) where T + +Given a delaunay triangulation `DT` of `all_pts[1:(end - 1)]`, find the indices of +the simplices that potentially intersect with the image simplex with index `idx`. +""" +function idxs_potentially_intersecting_simplices(all_pts, DT::DelaunayTriangulation, idx::Int) + # Vector that will hold the indices of the simplices potentially + # intersecting with the image of simplex #idx + inds_potential_simplices = Int[] + + n_simplices = length(DT) + + original_pts = all_pts[1:(end - 1)] + forward_pts = all_pts[2:end] + + simplices = [MutableSimplex(original_pts[DT[i]]) for i = 1:n_simplices] + image_simplices = [MutableSimplex(forward_pts[DT[i]]) for i = 1:n_simplices] + + cs = [centroid(simplices[i]) for i = 1:n_simplices] + cs_im = [centroid(image_simplices[i]) for i = 1:n_simplices] + + rs = [radius(simplices[i]) for i = 1:n_simplices] + rs_im = [radius(image_simplices[i]) for i = 1:n_simplices] + + @inbounds for i = 1:n_simplices + δ = transpose(cs_im[idx] .- cs[i]) * (cs_im[idx] .- cs[i]) .- ((rs_im[idx] + rs[i])^2) + if δ[1] < 0 + push!(inds_potential_simplices, i) + end + end + + return inds_potential_simplices +end + +""" Get the simplices at the given indices. """ +function get_simplices_at_inds!(simps, inds::Vector{Int}, simplices) + for i in 1:length(inds) + simps[i] = simplices[:, inds[i]] + end +end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl b/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl new file mode 100644 index 000000000..988302950 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl @@ -0,0 +1,131 @@ +export SimplexPoint, transferoperator, SimplexPoint +import StaticArrays: SizedVector, SizedMatrix, MMatrix + +include("helper_functions.jl") +include("../delaunay_triangulations/DelaunayTriangulations.jl") + +""" + SimplexPoint <: Entropies.BinningProbabilitiesEstimator + +A probability estimator based on the transfer operator, using a triangulated partition +and approximate simplex intersections to compute transition probabilities [^Diego2019]. + +To use this estimator, the Simplices.jl package must be brought into scope by doing +`using Simplices` after running `using Entropies`. + +[^Diego2019]: Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212. +""" +struct SimplexPoint <: BinningProbabilitiesEstimator + bc::String + + function SimplexPoint(bc::String = "circular") + isboundarycondition(bc, "triangulation") || error("Boundary condition '$bc' not valid.") + new(bc) + end +end +Base.show(io::IO, se::SimplexPoint) = print(io, "SimplexPoint{$(se.bc)}") + + +""" Generate a transopergenerator for an exact simplex estimator.""" +function transopergenerator(pts, method::TransferOperator{<:SimplexPoint}) + # modified points, where the image of each point is guaranteed to lie within the convex hull of the previous points + invariant_pts = invariantize(pts) + + # triangulation of the invariant points. The last point is excluded, so that + # the last vertex also can be mapped forward one step in time. + triang = DelaunayTriangulation(invariant_pts[1:end-1]) + + init = (invariant_pts = invariant_pts, triang = triang) + + TransferOperatorGenerator(method, pts, init) +end + + + +""" +n: The number of points sampled inside each simplex. +randomsampling: Sample randomly (`random == true`), or using an even simplex splitting routine (`random == false`) +""" +function (tog::TransferOperatorGenerator{T})(; tol = 1e-8, randomsampling::Bool = false, + n::Int = 100) where T <: TransferOperator{SimplexPoint} + pts, inds = getfield.(Ref(tog.init), (:invariant_pts, :triang)) + D = length(pts[1]) + N = length(inds) + ϵ = tol / N + + # Prepare memory-efficient representations of the simplices + simplices = reshape_simplices(pts[1:(end- 1)], inds) + image_simplices = reshape_simplices(pts[2:end], inds) + + #= + Compute the convex expansions coefficients needed to generate points + within a simplex. Also, update number of points if a shape-preserving + simplex subdivision algorithm was used (in that case, we get a + *minimum* of `n_sample_pts` set of coefficients. + =# + convex_coeffs = subsample_coeffs(D, n, randomsampling) + n_coeffs::Int = size(convex_coeffs, 2) + convex_coeffs = [convex_coeffs[:, i] for i = 1:n_coeffs] + + # Pre-allocated arrays (SizedArrays, for efficiency) + L = (D+1)^2 + NVERTICES = D+1 + pt = SizedMatrix{1, D}(zeros(1, D)) + s_arr = SizedVector{L}(zeros(L)) + signs = SizedVector{NVERTICES}(zeros(NVERTICES)) + + # Re-arrange simplices so that look-up is a bit more efficient + simplex_arrs = Vector{Array{Float64, 2}}(undef, N) + imsimplex_arrs = Vector{Array{Float64, 2}}(undef, N) + for i in 1:N + simplex_arrs[i] = hcat(pts[inds[i]]...,) + imsimplex_arrs[i] = hcat(pts[inds[i] .+ 1]...,) + end + + # Pre-allocate some information to avoid memory allocation when + # computing the fractional overlaps. These are a temporary matrix + # which we use to compute determinants (to check if a point lies + # inside a simplex), and the indices needed to access its values + # at different stages of the checking stage (of whether a point + # lies inside a simplex). Pre-allocating this information and + # passing it on reduces runtime and memory allocation + # by half an order of magnitude. + temp_arr = rand(MMatrix{D + 1, D + 1}) + starts = [((D + 1)*i)+1 for i = 0:D] + stops = starts .+ D + + + # The Markov matrix + M = zeros(Float64, N, N) + + for i in 1:N + # Find the indices of the simplices who potentially intersect + # with the image of simplex #i + idxs = idxs_potentially_intersecting_simplices(pts, inds, i) + Sj = Vector{AbstractArray}(undef, length(idxs)) + get_simplices_at_inds!(Sj, idxs, simplices) + + # Generate points within the image of simplex #i and check + # what fraction of those points fall into the simplices. + # We approximate the intersecting volume between the + # image simplex and the individual simplices by those + # fractions. + @views is = imsimplex_arrs[i] + for k in 1:n_coeffs + pt = transpose(convex_coeffs[k]) * transpose(is) + innerloop_optim!(idxs, signs, s_arr, Sj, pt, D, M, i, temp_arr, starts, stops) + end + end + # Need to normalise, because all we have up until now is counts + # of how many points inside the image simplex falls into + # the simplices. + params = (tol = tol, randomsampling = randomsampling, n = n) + M_normalized = transpose(M) ./ n_coeffs + return TransferOperatorApproximation(tog, M_normalized, params) + +end + +function transferoperator(pts, method::TransferOperator{<:SimplexPoint}; tol = 1e-8, randomsampling::Bool = false, n::Int = 100) + tog = transopergenerator(pts, method) + tog(tol = tol, randomsampling = randomsampling, n = n) +end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/point/helper_functions.jl b/src/binning_based/transferoperator/triangular/point/helper_functions.jl new file mode 100644 index 000000000..12ff94fba --- /dev/null +++ b/src/binning_based/transferoperator/triangular/point/helper_functions.jl @@ -0,0 +1,178 @@ +import LinearAlgebra: det +import StaticArrays: SizedMatrix + +include("../delaunay_triangulations/DelaunayTriangulations.jl") +import .DelaunayTriangulations: DelaunayTriangulation + + +""" Re-zeros the array `a`. """ +function rezero!(a) + @inbounds for i in eachindex(a) + a[i] = 0.0 + end +end + +""" Fill the elements of vector `from` into vector `into`. """ +function fill_into!(into, from) + @inbounds for i in eachindex(into) + into[i] = from[i] + end +end + +""" Fill vector `into` at indices `inds` with the elements at `from[inds]`. """ +function fill_at_inds!(into, from, inds) + @inbounds for i in 1:length(inds) + into[inds[i]] = from[i] + end +end + +function mdet_fm_in_x_view!(M, v, n, starts, stops) + @inbounds for i = 0:(n - 1) + M[starts[i+1]:stops[i+1]] = view(v, starts[i+1]:stops[i+1]) + end + @fastmath det(M) +end + + +function iscontained!(signs, s_arr, sx, point, dim, temp_arr, starts, stops) + # Redefine the temporary simplex. This is in-place, so we don't allocate + # memory. We could also have re-initialised `signs`, but since we're never + # comparing more than two consecutive signs, this is not necessary. + rezero!(s_arr) + rezero!(signs) + fill_into!(s_arr, sx) + + #Replace first vertex with the point + fill_at_inds!(s_arr, point, 1:dim) + + # Signed volume + signs[1] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops)) + + rezero!(s_arr) + fill_into!(s_arr, sx) #reset + + for κ = 2:dim # Check remaining signs and stop if sign changes + # Replace the ith vertex with the point we're cheking (leaving the + # 1 appended to Vi intact.) + idxs = ((dim + 1)*(κ - 1)+1):((dim + 1)*(κ - 1)+ 1 + dim - 1) + fill_at_inds!(s_arr, point, idxs) # ith change + + signs[κ] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops)) + + if !(signs[κ-1] == signs[κ]) + return false + end + + rezero!(s_arr) + fill_into!(s_arr, sx) + end + + # Last the last vertex with the point in question + idxs = ((dim + 1)*(dim)+1):((dim+1)^2-1) + fill_at_inds!(s_arr, point, idxs) + + signs[end] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops)) + + if !(signs[end-1] == signs[end]) + return false + else + return true + end +end + +function innerloop_optim!(inds::Vector{Int}, signs, s_arr, Sj, pt, dim::Int, M, i::Int, temp_arr, starts, stops) + for j in 1:length(inds) + if iscontained!(signs, s_arr, Sj[j], pt, dim, temp_arr, starts, stops) + M[inds[j], i] += 1.0 + end + end +end + +""" + reshape_simplices(pts::Vector{Vector{T}}, DT::DelaunayTriangulation) where T + +Creates alternative representations of the simplices that allow for efficient +(mostly non-allocating) checking if a point lies inside a simplex. + +This function adds some information needed when calculating the orientations +of the simplices *before* we start the computation (appends a row of ones on top of the +matrix representation of the simplex). Otherwise, the fallback is machinery in Simplices.jl, +which is much slower because the additional information must be added +`(n_simplices^2)*n_sample_pts` times (creating that many new matrices), which takes forever. + +Reshaping the simplices and appending that information beforehand avoids excessive +memory allocation. +""" +function reshape_simplices(pts, DT::DelaunayTriangulation) + n_simplices = length(DT.indices) + n_vertices = length(DT[1]) + dim = n_vertices - 1 + + # The [:, j, i] th entry of these two arrays holds the jth vertex of the + # ith simplex, but instead of having just `dim` vertices, we append a `1` + # to the end of the vectors. This allows for efficent (non-allocating) + # computation within the `contains_point_lessalloc!` function. If we instead + # would have appended the 1's inside that function, we would be performing + # memory-allocating operations, which are very expensive. Doing this instead + # gives orders of magnitude speed-ups for sufficiently large triangulations. + S1 = Array{Float64}(undef, dim + 1, dim + 1, n_simplices) + + # Collect simplices in the form of (dim+1)^2-length column vectors. This + # also helps with the + NDIM, N_SIMPLICES = (dim+1)^2, n_simplices + simplices = SizedMatrix{NDIM}{N_SIMPLICES}(zeros((NDIM, N_SIMPLICES))) + + @inbounds for i in 1:n_simplices + for j in 1:n_vertices + S1[:, j, i] = vcat(pts[DT[i][j]], 1.0) + end + + simplices[:, i] = reshape(S1[:, :, i], (dim+1)^2) + end + + return simplices +end + + +""" +idxs_potentially_intersecting_simplices(all_pts::Vector{Vector{T}}, + DT::DelaunayTriangulation, idx::Int) where T + +Given a delaunay triangulation `DT` of `all_pts[1:(end - 1)]`, find the indices of +the simplices that potentially intersect with the image simplex with index `idx`. +""" +function idxs_potentially_intersecting_simplices(all_pts, DT::DelaunayTriangulation, idx::Int) + # Vector that will hold the indices of the simplices potentially + # intersecting with the image of simplex #idx + inds_potential_simplices = Int[] + + n_simplices = length(DT) + + original_pts = all_pts[1:(end - 1)] + forward_pts = all_pts[2:end] + + simplices = [MutableSimplex(original_pts[DT[i]]) for i = 1:n_simplices] + image_simplices = [MutableSimplex(forward_pts[DT[i]]) for i = 1:n_simplices] + + cs = [centroid(simplices[i]) for i = 1:n_simplices] + cs_im = [centroid(image_simplices[i]) for i = 1:n_simplices] + + rs = [radius(simplices[i]) for i = 1:n_simplices] + rs_im = [radius(image_simplices[i]) for i = 1:n_simplices] + + @inbounds for i = 1:n_simplices + δ = transpose(cs_im[idx] .- cs[i]) * (cs_im[idx] .- cs[i]) .- ((rs_im[idx] + rs[i])^2) + if δ[1] < 0 + push!(inds_potential_simplices, i) + end + end + + return inds_potential_simplices +end + +""" Get the simplices at the given indices. """ +function get_simplices_at_inds!(simps, inds::Vector{Int}, simplices) + for i in 1:length(inds) + simps[i] = simplices[:, inds[i]] + end +end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/simplex_types/.ipynb_checkpoints/simplex_subsampling-checkpoint.jl b/src/binning_based/transferoperator/triangular/simplex_types/.ipynb_checkpoints/simplex_subsampling-checkpoint.jl new file mode 100644 index 000000000..42078cbbc --- /dev/null +++ b/src/binning_based/transferoperator/triangular/simplex_types/.ipynb_checkpoints/simplex_subsampling-checkpoint.jl @@ -0,0 +1,116 @@ +import Distributions: Uniform + +""" + tensordecomp(k::Int, d::Int) + +Decomposition of the integers 0:(k^p - 1) in powers of k. +""" +function tensordecomp(k::Int, d::Int) + + sequences = zeros(Int, k^d, d) + + for n = 0:k^d-1 + i = d + m = n + + while i > 0 + i = i - 1 + j = i + f = floor(Int, m / k^i) + + while j > 0 && f == 0 + j = j - 1 + f = floor(Integer, m / k^j) + end + + if f > 0 + sequences[n + 1, j + 1, ] = f + i = j + elseif f == 0 + sequences[n + 1, 1] = m + i = 0 + end + m = m - f * k^i + end + end + + return sequences +end + + +""" + even_sampling_rules(dim::Int, split_factor::Int) -> Array{Float64, 2} + +Generate rules for evenly distributed points within a simplex. To do this, +we perform a shape-preserving splitting of the simplex, given a splitting +factor. + +Returns the convex expansion coefficients of the points of the resulting +subsimplices in terms of the vertices of the original simplex. +""" +function even_sampling_rules(dim::Int, split_factor::Int) + + sequences::Array{Int, 2} = tensordecomp(split_factor, dim) + n_seq = size(sequences, 1) + + χ1 = sequences .* (dim + 1) + χ2 = repeat(transpose(collect(1:dim)), n_seq, 1) + χ::Array{Int, 2} = χ1 .+ χ2 + χ = sort(χ, dims=2) + + # Define multiplicity matrix M + M = zeros(Float64, size(χ, 1), size(χ, 2) + 1) + M[:, 1] = χ[:, 1] + M[:, 2:(end - 1)] = χ[:, 2:end] - χ[:, 1:(end - 1)] + M[:, end] = (dim+1)*split_factor * ones(size(χ, 1)) - χ[:, end] + + M = M ./ (split_factor * (dim + 1)) + + return copy(transpose(M)) +end + +export even_sampling_rules + +""" +Evenly sample points within a simplex by performing a shape-preserving +subdivision of the simplex with a given `split_factor`. If the simplex +lives in a space of dimension `dim`, the resulting number of points is +`split_factor`^(dim). +""" +function evenly_sample(simplex::AbstractArray{Float64, 2}, split_factor::Int) + dim = size(simplex, 2) + centroids_exp_coeffs = copy(transpose(even_sampling_rules(dim, split_factor))) + centroids_exp_coeffs * simplex +end + + +""" +Returns a matrix of convex coefficients to construct points contained within +a `dim`-dimensional simplex. Coefficients can either be random (generated +according to a uniform distribution) or be constructed such that the resulting +points are uniformly distributed within the simplex (generated by shape-preserving +splitting a generic simplex). The `sample_randomly` argument controls this +behaviour. + +The default is to not sample randomly (`sample_randomly = false`). In this case, +you won't get the exact number of convex coefficient sets you want, but the +smallest number of set of coefficients such that you have *at least* `n_pts` +convex coefficient sets. + +If `sample_randomly = true`, you get as many convex coefficient combinations as +`n_pts`. +""" +function subsample_coeffs(dim::Int, n_randpts::Int, sample_randomly::Bool) + # Create a set of convex coefficients to be used for all simplices + if sample_randomly + convex_coeffs = rand(Uniform(0, 1), dim + 1, n_randpts) + convex_coeffs .= convex_coeffs ./ sum(convex_coeffs, dims=1) #convex_coeffs ./ sum(convex_coeffs, dims=2) + else + minimum_split_factor = ceil(Int, n_randpts^(1 / dim))[1] + convex_coeffs = even_sampling_rules(dim, minimum_split_factor) + end + return convex_coeffs +end + +export subsample_coeffs + diff --git a/src/binning_based/transferoperator/triangular/simplex_types/.ipynb_checkpoints/simplex_types-checkpoint.jl b/src/binning_based/transferoperator/triangular/simplex_types/.ipynb_checkpoints/simplex_types-checkpoint.jl new file mode 100644 index 000000000..37b71b323 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/simplex_types/.ipynb_checkpoints/simplex_types-checkpoint.jl @@ -0,0 +1,5 @@ + +include("AbstractSimplex.jl") +include("ImmutableSimplex.jl") +include("MutableSimplex.jl") +include("simplex_subsampling.jl") diff --git a/src/binning_based/transferoperator/triangular/simplex_types/AbstractSimplex.jl b/src/binning_based/transferoperator/triangular/simplex_types/AbstractSimplex.jl new file mode 100644 index 000000000..de2694170 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/simplex_types/AbstractSimplex.jl @@ -0,0 +1,98 @@ +export + AbstractSimplex, + orientation, + intersect, + radius, + volume, + centroid + +import LinearAlgebra: det +import .Simplices: centroid +import .Simplices: radius +import .Simplices: orientation +import .Simplices: volume + + +abstract type AbstractSimplex{D, T} end + +# Return the i-th point as column vector +Base.getindex(s::AbstractSimplex, i) = s.vertices[i] +Base.getindex(s::AbstractSimplex, i, j) = s.vertices[i][j] +Base.getindex(s::AbstractSimplex, i::Colon, j) = hcat(s.vertices[j]...,) +Base.getindex(s::AbstractSimplex, i::Colon, j::Colon) = hcat(s.vertices...,) + +Base.length(s::AbstractSimplex) = length(s.vertices) +Base.size(s::AbstractSimplex) = length(s[1]), length(s) +Base.size(s::AbstractSimplex, i) = size(s)[i] +Base.IteratorSize(s::AbstractSimplex) = Base.HasLength() + +Base.firstindex(s::AbstractSimplex) = 1 +Base.lastindex(s::AbstractSimplex) = length(s) +Base.eachindex(s::AbstractSimplex) = Base.OneTo(length(s)) +Base.iterate(s::AbstractSimplex, state = 1) = iterate(s.vertices, state) + +########################################################### +# Vector{vertex} representation => Array representation +########################################################### +Base.Array(s::AbstractSimplex) = hcat(s.vertices...,) + +""" + dimension(s::AbstractSimplex) + +Get the dimension of a simplex. +""" +dimension(s::AbstractSimplex) = length(s[1]) +npoints(s::AbstractSimplex) = length(s) + +""" + nvertices(s::AbstractSimplex) + +Get the number of vertices of a simplices. +""" +nvertices(s::AbstractSimplex) = length(s) + +""" + orientation(s::AbstractSimplex) + +Compute the orientation of a simplex. *Convention: append rows of ones at top of vertex matrix*. +""" +function Simplices.orientation(s::AbstractSimplex) + det([ones(1, dimension(s) + 1); s[:, :]]) +end + +""" + volume(s::AbstractSimplex) + +Compute the unscaled volume of the simplex. Divide by factorial(dim) to get the true volume. +""" +Simplices.volume(s::AbstractSimplex) = abs(orientation(s)) + +function Simplices.centroid(s::AbstractSimplex) + D = dimension(s) + # Results in a dim-by-1 matrix, but we just want a vector, so drop the last dimension + centroid = dropdims(s[:, :] * (ones(D + 1, 1) / (D + 1)), dims = 2) +end + +""" + radius(s::AbstractSimplex) + +Compute the radius of a simplex. +""" +function Simplices.radius(s::AbstractSimplex) + D = dimension(s) + centroidmatrix = repeat(centroid(s), 1, D + 1) + radius = sqrt(maximum(ones(1, D) * ((s[:, :] .- centroidmatrix) .^ 2))) +end + + +function Base.intersect(s1::AbstractSimplex, s2::AbstractSimplex) + Simplices.simplexintersection(s1[:, :], s2[:, :]) +end + +""" + ∩(s1::AbstractSimplex, s2::AbstractSimplex) + +Compute the volume intersection between two simplices. +""" +∩(s1::AbstractSimplex, s2::AbstractSimplex) = intersect(s1, s2) + diff --git a/src/binning_based/transferoperator/triangular/simplex_types/ImmutableSimplex.jl b/src/binning_based/transferoperator/triangular/simplex_types/ImmutableSimplex.jl new file mode 100644 index 000000000..649dd24cc --- /dev/null +++ b/src/binning_based/transferoperator/triangular/simplex_types/ImmutableSimplex.jl @@ -0,0 +1 @@ +abstract type AbstractImmutableSimplex{D, T} <: AbstractSimplex{D, T} end diff --git a/src/binning_based/transferoperator/triangular/simplex_types/MutableSimplex.jl b/src/binning_based/transferoperator/triangular/simplex_types/MutableSimplex.jl new file mode 100644 index 000000000..e87475f2b --- /dev/null +++ b/src/binning_based/transferoperator/triangular/simplex_types/MutableSimplex.jl @@ -0,0 +1,60 @@ +export MutableSimplex + +abstract type AbstractMutableSimplex{D, T} <: AbstractSimplex{D, T} end + +""" + MutableSimplex{D, T} + +Simplex where vertices are represented by some type of abstract mutable vector. +""" +mutable struct MutableSimplex{D, T} <: AbstractMutableSimplex{D, T} + vertices::Vector{Vector{T}} + + function MutableSimplex(pts::Vector{<:AbstractVector{T}}) where {T} + if !(length(pts) == length(pts[1]) + 1) + err = """ The input cannot be converted to a simplex. + Vertices need to have `dim` elements, and there needs to be `dim + 1` vertices. + """ + throw(DomainError(pts, err)) end + dim = length(pts[1]) + new{dim, T}([pts[i] for i = 1:length(pts)]) + end + + function MutableSimplex(pts::AbstractArray{T, 2}) where {T} + s = size(pts) + + if (maximum(s) - minimum(s)) != 1 + err = """ The input cannot be converted to a simplex. + size(pts) must be (dim, dim + 1) or (dim + 1, dim). + """ + throw(DomainError(pts, err)) + end + + if s[1] > s[2] # Rows are points + dim = s[2] + return new{dim, T}([pts[i, :] for i = 1:maximum(s)]) + end + + if s[2] > s[1] # Columns are points + dim = s[1] + return new{dim, T}([pts[:, i] for i = 1:maximum(s)]) + end + end +end + +# Overwriting the i-th vertex +function Base.setindex!(simplex::MutableSimplex, v, i) + simplex[i] .= v +end + +# Overwriting elements of the i-th vertex +function Base.setindex!(simplex::MutableSimplex, v, i::Int, j) + if length(v) != length(j) + err = """ + Trying to overwrite elements $j of vertex $i with $v, which does not + have the same number of elements as the target. + """ + throw(ArgumentError(err)) + end + simplex[i][j] = v +end diff --git a/src/binning_based/transferoperator/triangular/simplex_types/simplex_subsampling.jl b/src/binning_based/transferoperator/triangular/simplex_types/simplex_subsampling.jl new file mode 100644 index 000000000..42078cbbc --- /dev/null +++ b/src/binning_based/transferoperator/triangular/simplex_types/simplex_subsampling.jl @@ -0,0 +1,116 @@ +import Distributions: Uniform + +""" + tensordecomp(k::Int, d::Int) + +Decomposition of the integers 0:(k^p - 1) in powers of k. +""" +function tensordecomp(k::Int, d::Int) + + sequences = zeros(Int, k^d, d) + + for n = 0:k^d-1 + i = d + m = n + + while i > 0 + i = i - 1 + j = i + f = floor(Int, m / k^i) + + while j > 0 && f == 0 + j = j - 1 + f = floor(Integer, m / k^j) + end + + if f > 0 + sequences[n + 1, j + 1, ] = f + i = j + elseif f == 0 + sequences[n + 1, 1] = m + i = 0 + end + m = m - f * k^i + end + end + + return sequences +end + + +""" + even_sampling_rules(dim::Int, split_factor::Int) -> Array{Float64, 2} + +Generate rules for evenly distributed points within a simplex. To do this, +we perform a shape-preserving splitting of the simplex, given a splitting +factor. + +Returns the convex expansion coefficients of the points of the resulting +subsimplices in terms of the vertices of the original simplex. +""" +function even_sampling_rules(dim::Int, split_factor::Int) + + sequences::Array{Int, 2} = tensordecomp(split_factor, dim) + n_seq = size(sequences, 1) + + χ1 = sequences .* (dim + 1) + χ2 = repeat(transpose(collect(1:dim)), n_seq, 1) + χ::Array{Int, 2} = χ1 .+ χ2 + χ = sort(χ, dims=2) + + # Define multiplicity matrix M + M = zeros(Float64, size(χ, 1), size(χ, 2) + 1) + M[:, 1] = χ[:, 1] + M[:, 2:(end - 1)] = χ[:, 2:end] - χ[:, 1:(end - 1)] + M[:, end] = (dim+1)*split_factor * ones(size(χ, 1)) - χ[:, end] + + M = M ./ (split_factor * (dim + 1)) + + return copy(transpose(M)) +end + +export even_sampling_rules + +""" +Evenly sample points within a simplex by performing a shape-preserving +subdivision of the simplex with a given `split_factor`. If the simplex +lives in a space of dimension `dim`, the resulting number of points is +`split_factor`^(dim). +""" +function evenly_sample(simplex::AbstractArray{Float64, 2}, split_factor::Int) + dim = size(simplex, 2) + centroids_exp_coeffs = copy(transpose(even_sampling_rules(dim, split_factor))) + centroids_exp_coeffs * simplex +end + + +""" +Returns a matrix of convex coefficients to construct points contained within +a `dim`-dimensional simplex. Coefficients can either be random (generated +according to a uniform distribution) or be constructed such that the resulting +points are uniformly distributed within the simplex (generated by shape-preserving +splitting a generic simplex). The `sample_randomly` argument controls this +behaviour. + +The default is to not sample randomly (`sample_randomly = false`). In this case, +you won't get the exact number of convex coefficient sets you want, but the +smallest number of set of coefficients such that you have *at least* `n_pts` +convex coefficient sets. + +If `sample_randomly = true`, you get as many convex coefficient combinations as +`n_pts`. +""" +function subsample_coeffs(dim::Int, n_randpts::Int, sample_randomly::Bool) + # Create a set of convex coefficients to be used for all simplices + if sample_randomly + convex_coeffs = rand(Uniform(0, 1), dim + 1, n_randpts) + convex_coeffs .= convex_coeffs ./ sum(convex_coeffs, dims=1) #convex_coeffs ./ sum(convex_coeffs, dims=2) + else + minimum_split_factor = ceil(Int, n_randpts^(1 / dim))[1] + convex_coeffs = even_sampling_rules(dim, minimum_split_factor) + end + return convex_coeffs +end + +export subsample_coeffs + diff --git a/src/binning_based/transferoperator/triangular/simplex_types/simplex_types.jl b/src/binning_based/transferoperator/triangular/simplex_types/simplex_types.jl new file mode 100644 index 000000000..37b71b323 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/simplex_types/simplex_types.jl @@ -0,0 +1,5 @@ + +include("AbstractSimplex.jl") +include("ImmutableSimplex.jl") +include("MutableSimplex.jl") +include("simplex_subsampling.jl") diff --git a/src/binning_based/transferoperator/utils.jl b/src/binning_based/transferoperator/utils.jl new file mode 100644 index 000000000..8cac378b8 --- /dev/null +++ b/src/binning_based/transferoperator/utils.jl @@ -0,0 +1,19 @@ +function invariantize(pts::AbstractDataset, bc = "circular") + invariant_pts = copy(pts.data) + if bc == "circular" + return push!(invariant_pts, pts.data[1]) + elseif bc == "random" + return push!(invariant_pts, pts.data[rand(1:length(pts))]) + end + return invariant_pts +end + +function isboundarycondition(bc, method::String) + if method == "grid" + bc ∈ ["circular", "random"] + elseif method ∈ ["triangulation"] + bc ∈ ["circular", "random"] + else + error("method $method not defined") + end +end \ No newline at end of file diff --git a/src/binning_based/rectangular/VisitationFrequency.jl b/src/binning_based/visitation_frequency/VisitationFrequency.jl similarity index 89% rename from src/binning_based/rectangular/VisitationFrequency.jl rename to src/binning_based/visitation_frequency/VisitationFrequency.jl index 7467ff73d..fadd196d8 100644 --- a/src/binning_based/rectangular/VisitationFrequency.jl +++ b/src/binning_based/visitation_frequency/VisitationFrequency.jl @@ -1,5 +1,8 @@ export VisitationFrequency, probabilities -import DelayEmbeddings: Dataset, AbstractDataset +import DelayEmbeddings: AbstractDataset + +include("count_box_visits.jl") +include("histogram_estimation.jl") """ VisitationFrequency(r::RectangularBinning) <: BinningProbabilitiesEstimator diff --git a/src/binning_based/rectangular/count_box_visits.jl b/src/binning_based/visitation_frequency/count_box_visits.jl similarity index 100% rename from src/binning_based/rectangular/count_box_visits.jl rename to src/binning_based/visitation_frequency/count_box_visits.jl diff --git a/src/binning_based/rectangular/histogram_estimation.jl b/src/binning_based/visitation_frequency/histogram_estimation.jl similarity index 100% rename from src/binning_based/rectangular/histogram_estimation.jl rename to src/binning_based/visitation_frequency/histogram_estimation.jl From 37ee36a36b16203179420914e1392ff5f3b00d2d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 02:59:04 +0100 Subject: [PATCH 02/29] Docs and type requirements --- Project.toml | 3 ++- docs/src/TransferOperator.md | 9 +++++++++ src/binning_based/transferoperator/common.jl | 6 +++--- .../transferoperator/rectangular/rectangular.jl | 8 +++++--- .../transferoperator/triangular/exact/SimplexExact.jl | 4 ++-- .../transferoperator/triangular/point/SimplexPoint.jl | 9 ++++----- 6 files changed, 25 insertions(+), 14 deletions(-) diff --git a/Project.toml b/Project.toml index eeed0acd9..bb1102d8b 100644 --- a/Project.toml +++ b/Project.toml @@ -35,6 +35,7 @@ julia = "1.5" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" [targets] -test = ["Test"] +test = ["Test", "Simplices"] diff --git a/docs/src/TransferOperator.md b/docs/src/TransferOperator.md index 7a2e2bc95..e7f7bc4d2 100644 --- a/docs/src/TransferOperator.md +++ b/docs/src/TransferOperator.md @@ -2,10 +2,19 @@ ```@docs TransferOperator +``` + +## Triangulation estimators + +```@docs SimplexExact SimplexPoint ``` +## Rectangular estimators + +For rectangular binnings, use [`RectangularBinning`](@ref). + ## Utility methods/types ```@docs diff --git a/src/binning_based/transferoperator/common.jl b/src/binning_based/transferoperator/common.jl index 73c753374..7a17de2fb 100644 --- a/src/binning_based/transferoperator/common.jl +++ b/src/binning_based/transferoperator/common.jl @@ -175,13 +175,13 @@ When `R <: SimplexExact`, parameters stored are the simplex intersection toleran See also: [`RectangularBinning`](@ref), [`SimplexPoint`](@ref), [`SimplexExact`](@ref). """ -struct TransferOperatorApproximation{G<:TransferOperator{<:BinningProbabilitiesEstimator}, T<:Real} +struct TransferOperatorApproximation{G<:TransferOperator, T} generator::TransferOperatorGenerator{G} - transfermatrix::AbstractArray{T, 2} + transfermatrix::T params end -function Base.show(io::IO, DT::TransferOperatorApproximation{G, T}) where {G <: TransferOperator{R}, T} where R +function Base.show(io::IO, DT::TransferOperatorApproximation{G, T}) where {G <: TransferOperator, T} summary = "TransferOperatorApproximation{$G, $T}" println(io, summary) end diff --git a/src/binning_based/transferoperator/rectangular/rectangular.jl b/src/binning_based/transferoperator/rectangular/rectangular.jl index 8d1a6df61..298bf8d68 100644 --- a/src/binning_based/transferoperator/rectangular/rectangular.jl +++ b/src/binning_based/transferoperator/rectangular/rectangular.jl @@ -53,6 +53,8 @@ end inds_in_terms_of_unique(x::AbstractDataset) = inds_in_terms_of_unique(x.data) +# TODO: transferoperatorgenerator for rectangualrbinning! + """ transferoperator(pts::AbstractDataset{D, T}, ϵ::RectangularBinning) → TransferOperatorApproximationRectangular @@ -74,11 +76,11 @@ transferoperator(orbit, RectangularBinning(10)) See also: [`RectangularBinning`](@ref). """ -function transferoperator(pts::AbstractDataset{D, T}, to::TransferOperator{RectangularBinning}; +function transferoperator(pts::AbstractDataset{D, T}, to::TransferOperator{<:RectangularBinning}; boundary_condition = :circular) where {D, T<:Real} L = length(pts) - mini, edgelengths = Entropies.minima_edgelengths(pts, ϵ) + mini, edgelengths = Entropies.minima_edgelengths(pts, to.ϵ) # The L points visits a total of L bins, which are the following bins: visited_bins = Entropies.encode_as_bin(pts, mini, edgelengths) @@ -185,7 +187,7 @@ function transferoperator(pts::AbstractDataset{D, T}, to::TransferOperator{Recta # row/column of the transfer operator unique!(visited_bins) bins = [β .* edgelengths .+ mini for β in visited_bins] - params = (ϵ = ϵ, mini = mini, edgelengths = edgelengths, bins = bins, + params = (mini = mini, edgelengths = edgelengths, bins = bins, sort_idxs = sort_idxs, visitors) TransferOperatorApproximation(to, M, params) diff --git a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl index fb6c08147..659284593 100644 --- a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl +++ b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl @@ -4,7 +4,7 @@ import Entropies.invariantmeasure """ - SimplexExact <: TriangulationBasedTransferOperator + SimplexExact A transfer operator estimator using a triangulation partition and exact simplex intersections[^Diego2019]. @@ -17,7 +17,7 @@ simplex intersections, this estimator is slow compared to [`SimplexPoint`](@ref) [^Diego2019]: Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212. """ -struct SimplexExact <: BinningProbabilitiesEstimator +struct SimplexExact bc::String function SimplexExact(bc::String = "circular") diff --git a/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl b/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl index 988302950..616d69488 100644 --- a/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl +++ b/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl @@ -2,20 +2,19 @@ export SimplexPoint, transferoperator, SimplexPoint import StaticArrays: SizedVector, SizedMatrix, MMatrix include("helper_functions.jl") -include("../delaunay_triangulations/DelaunayTriangulations.jl") """ - SimplexPoint <: Entropies.BinningProbabilitiesEstimator + SimplexPoint -A probability estimator based on the transfer operator, using a triangulated partition -and approximate simplex intersections to compute transition probabilities [^Diego2019]. +A transfer operator estimator using a triangulated partition, and using +approximate simplex intersections to compute transition probabilities [^Diego2019]. To use this estimator, the Simplices.jl package must be brought into scope by doing `using Simplices` after running `using Entropies`. [^Diego2019]: Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212. """ -struct SimplexPoint <: BinningProbabilitiesEstimator +struct SimplexPoint bc::String function SimplexPoint(bc::String = "circular") From c11e75078f31d79de1f4dd5eab476976a04021cc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 13:31:55 +0100 Subject: [PATCH 03/29] Remove ipynb files --- ...bstractDelaunayTriangulation-checkpoint.jl | 37 ---- .../DelaunayTriangulation-checkpoint.jl | 20 -- .../DelaunayTriangulations-checkpoint.jl | 13 -- .../.ipynb_checkpoints/addnoise-checkpoint.jl | 105 ----------- .../SimplexExact-checkpoint.jl | 154 --------------- .../SimplexPoint-checkpoint.jl | 131 ------------- .../helper_functions-checkpoint.jl | 175 ------------------ 7 files changed, 635 deletions(-) delete mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/AbstractDelaunayTriangulation-checkpoint.jl delete mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulation-checkpoint.jl delete mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulations-checkpoint.jl delete mode 100644 src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/addnoise-checkpoint.jl delete mode 100644 src/binning_based/transferoperator/triangular/exact/.ipynb_checkpoints/SimplexExact-checkpoint.jl delete mode 100644 src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/SimplexPoint-checkpoint.jl delete mode 100644 src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/helper_functions-checkpoint.jl diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/AbstractDelaunayTriangulation-checkpoint.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/AbstractDelaunayTriangulation-checkpoint.jl deleted file mode 100644 index c99f41bf0..000000000 --- a/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/AbstractDelaunayTriangulation-checkpoint.jl +++ /dev/null @@ -1,37 +0,0 @@ -#################################### -# Abstract type -#################################### - -abstract type AbstractDelaunayTriangulation end - -ADT = AbstractDelaunayTriangulation -Base.size(DT::ADT) = (length(DT.indices[1]) - 1, length(DT.indices)) -Base.length(DT::ADT) = length(DT.indices) -dimension(DT::ADT) = length(DT.indices[1]) - 1 -nsimplices(DT::ADT) = length(DT.indices) - -# Indexing -Base.getindex(DT::ADT, i) = DT.indices[i] -Base.getindex(DT::ADT, i::Colon, j::Colon) = hcat(DT.indices...,) -Base.getindex(DT::ADT, i::Colon, j) = hcat(DT[j]...,) -Base.getindex(DT::ADT, i::Colon, j::Int) = hcat(DT[j]) - -Base.firstindex(DT::ADT) = 1 -Base.lastindex(DT::ADT) = length(DT) - -Base.eachindex(s::ADT) = Base.OneTo(length(s)) -Base.iterate(s::ADT, state = 1) = iterate(s.indices, state) - -function summarise(DT::AbstractDelaunayTriangulation) - _type = typeof(DT) - n_simplices = nsimplices(DT) - D = dimension(DT) - summary = "$D-dimensional $_type with $n_simplices simplices" -end - -Base.show(io::IO, DT::AbstractDelaunayTriangulation) = println(io, summarise(DT)) - - -export -dimension, -nsimplices diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulation-checkpoint.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulation-checkpoint.jl deleted file mode 100644 index 319f96fc5..000000000 --- a/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulation-checkpoint.jl +++ /dev/null @@ -1,20 +0,0 @@ -export DelaunayTriangulation - -import StaticArrays: SVector, MVector - -include("addnoise.jl") - -struct DelaunayTriangulation <: AbstractDelaunayTriangulation - indices::Vector{Vector{Int32}} - - function DelaunayTriangulation(vertices::Vector{T}; joggle = 1e-8) where T <: Union{AbstractVector, SVector, MVector} - # Slightly joggle points to avoid problems with QHull - if joggle > 0 - addnoise!(vertices; joggle_factor = joggle) - end - simplex_inds = Simplices.Delaunay.delaunay(vertices) - new(simplex_inds) - end - - DelaunayTriangulation(vertices::Dataset; joggle = 1e-8) = DelaunayTriangulation(vertices.data) -end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulations-checkpoint.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulations-checkpoint.jl deleted file mode 100644 index 4a489cba6..000000000 --- a/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/DelaunayTriangulations-checkpoint.jl +++ /dev/null @@ -1,13 +0,0 @@ -module DelaunayTriangulations - using Simplices - include("AbstractDelaunayTriangulation.jl") - include("DelaunayTriangulation.jl") - #include("plot_recipes/plot_recipes.jl") -end - -""" - DelaunayTriangulations - -A module handling the creation of Delaunay triangulations from data. -""" -DelaunayTriangulations diff --git a/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/addnoise-checkpoint.jl b/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/addnoise-checkpoint.jl deleted file mode 100644 index 8b944c4d6..000000000 --- a/src/binning_based/transferoperator/triangular/delaunay_triangulations/.ipynb_checkpoints/addnoise-checkpoint.jl +++ /dev/null @@ -1,105 +0,0 @@ - -import Distributions: Uniform -import Statistics: std -import StaticArrays: MArray, SArray, MVector, SVector - -""" - addnoise!(pts, joggle_factor) - -Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`. -""" -function addnoise!(pts::AbstractArray{T, 2}; joggle_factor = 0.001) where T - - D = minimum(size(pts)) - npts = maximum(size(pts)) - - if size(pts, 1) > size(pts, 2) - dim = 1 - else - dim = 2 - end - - # Scale standard deviation along each axis by joggle factor - σ = joggle_factor .* std(pts, dims = dim) - - for i in 1:D - r = [rand(Uniform(-σ[i], σ[i])) for pt in 1:npts] - if dim == 1 - pts[:, i] .+= r - elseif dim == 2 - pts[i, :] .+= r - end - end -end - -""" - addnoise!(pts, joggle_factor) - -Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`. -""" -function addnoise!(pts::Vector{Vector{T}}; joggle_factor = 1e-8) where T - - D = length(pts[1]) - npts = length(pts) - - # Scale standard deviation along each axis by joggle factor - σ = joggle_factor - - for i in 1:D - r = rand(Uniform(-σ, σ)) - pts[i] .+= r - end -end - -import DelayEmbeddings: - Dataset - -function addnoise!(pts::Dataset; joggle_factor = 1e-8) - - D = size(pts, 2) - npts = length(pts) - - # Scale standard deviation along each axis by joggle factor - σ = joggle_factor - - Dataset([pts[i] .+ rand(Uniform(-σ, σ)) for i = 1:npts]) -end - - - - -""" - addnoise!(pts, joggle_factor) - -Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`. -""" -function addnoise!(pts::Vector{MVector{D, T}}; joggle_factor = 1e-8) where {D, T} - - npts = length(pts) - - # Scale standard deviation along each axis by joggle factor - σ = joggle_factor - - for i in 1:D - r = rand(Uniform(-σ, σ)) - pts[i] .+= r - end -end - - - - -""" - addnoise!(pts, joggle_factor) - -Adding uniformly distributed noise to each observation of maximum magnitude `joggle_factor`. -""" -function addnoise!(pts::Vector{SVector{D, T}}; joggle_factor = 1e-8) where {D, T} - - npts = length(pts) - - # Scale standard deviation along each axis by joggle factor - σ = joggle_factor - - [pts[i] .+ rand(Uniform(-σ, σ)) for i = 1:npts] -end diff --git a/src/binning_based/transferoperator/triangular/exact/.ipynb_checkpoints/SimplexExact-checkpoint.jl b/src/binning_based/transferoperator/triangular/exact/.ipynb_checkpoints/SimplexExact-checkpoint.jl deleted file mode 100644 index ec913f617..000000000 --- a/src/binning_based/transferoperator/triangular/exact/.ipynb_checkpoints/SimplexExact-checkpoint.jl +++ /dev/null @@ -1,154 +0,0 @@ -export SimplexExact, invariantmeasure -import Entropies.invariantmeasure - - - -""" - SimplexExact() <: TriangulationBasedTransferOperator - -A transfer operator estimator using a triangulation partition and exact -simplex intersections[^Diego2019]. - -To use this estimator, the Simplices.jl -package must be brought into scope by doing `using Simplices` after -`using PerronFrobenius`. *Note: due to computing the exact simplex intersection, -this estimator is slow compared to `SimplexApprox`.* - -[^Diego2019]: Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212. -""" -struct SimplexExact <: BinningProbabilitiesEstimator - bc::String - - function SimplexExact(bc::String = "circular") - isboundarycondition(bc, "triangulation") || error("Boundary condition '$bc' not valid.") - new(bc) - end -end -Base.show(io::IO, se::SimplexExact) = print(io, "SimplexExact{$(se.bc)}") - - -""" Generate a transopergenerator for an exact simplex estimator.""" -function transopergenerator(pts, method::TransferOperator{<:SimplexExact}) - # modified points, where the image of each point is guaranteed to lie within the convex hull of the previous points - invariant_pts = invariantize(pts) - - # triangulation of the invariant points. The last point is excluded, so that - # the last vertex also can be mapped forward one step in time. - triang = DelaunayTriangulation(invariant_pts[1:end-1]) - - init = (invariant_pts = invariant_pts, triang = triang) - - TransferOperatorGenerator(method, pts, init) -end - -function (tog::TransferOperatorGenerator{T})(; tol = 1e-8) where T <: TransferOperator{SimplexExact} - invariant_pts, triang = getfield.(Ref(tog.init), (:invariant_pts, :triang)) - - D = length(invariant_pts[1]) - N = length(triang) - ϵ = tol / N - - # Pre-allocate simplex and its image - image_simplex = MutableSimplex(zeros(Float64, D, D+1)) - simplex = MutableSimplex(zeros(Float64, D, D+1)) - M = zeros(Float64, N, N) - - for j in 1:N - for k = 1:(D + 1) - # Get the vertices of the image simplex - image_simplex[k] .= invariant_pts[triang[j][k] + 1] - end - - imvol = abs(orientation(image_simplex)) - - for i = 1:N - for k = 1:(D + 1) - # Get the vertices of the source simplex - simplex[k] .= invariant_pts[triang[i][k]] - end - - # Only compute the entry of the transfer matrix - # if simplices are of sufficient size. - vol = abs(orientation(simplex)) - - if vol * imvol > 0 && (vol/imvol) > ϵ - M[j, i] = intersect(simplex, image_simplex) / imvol - end - end - end - params = (; tol = tol) - TransferOperatorApproximation(tog, M, params) -end - -function transferoperator(pts, method::SimplexExact; tol::Real = 1e-8) - tog = transopergenerator(pts, method) - tog(tol = tol) -end - - -function invariantmeasure(to::TransferOperatorApproximation{G, T}; - N::Int = 200, - tolerance::Float64 = 1e-8, - delta::Float64 = 1e-8) where {G <: TransferOperator{R}, T} where R <: SimplexExact - - TO = to.transfermatrix - #= - # Start with a random distribution `Ρ` (big rho). Normalise it so that it - # sums to 1 and forms a true probability distribution over the partition elements. - =# - Ρ = rand(Float64, 1, size(to.transfermatrix, 1)) - Ρ = Ρ ./ sum(Ρ, dims = 2) - - #= - # Start estimating the invariant distribution. We could either do this by - # finding the left-eigenvector of M, or by repeated application of M on Ρ - # until the distribution converges. Here, we use the latter approach, - # meaning that we iterate until Ρ doesn't change substantially between - # iterations. - =# - distribution = Ρ * TO - - distance = norm(distribution - Ρ) / norm(Ρ) - - check = floor(Int, 1 / delta) - check_pts = floor.(Int, transpose(collect(1:N)) ./ check) .* transpose(collect(1:N)) - check_pts = check_pts[check_pts .> 0] - num_checkpts = size(check_pts, 1) - check_pts_counter = 1 - - counter = 1 - while counter <= N && distance >= tolerance - counter += 1 - Ρ = distribution - - # Apply the Markov matrix to the current state of the distribution - distribution = Ρ * to.transfermatrix - - if (check_pts_counter <= num_checkpts && - counter == check_pts[check_pts_counter]) - - check_pts_counter += 1 - colsum_distribution = sum(distribution, dims = 2)[1] - if abs(colsum_distribution - 1) > delta - distribution = distribution ./ colsum_distribution - end - end - - distance = norm(distribution - Ρ) / norm(Ρ) - end - distribution = dropdims(distribution, dims = 1) - - # Do the last normalisation and check - colsum_distribution = sum(distribution) - - if abs(colsum_distribution - 1) > delta - distribution = distribution ./ colsum_distribution - end - - # Find partition elements with strictly positive measure. - δ = tolerance/size(to.transfermatrix, 1) - inds_nonzero = findall(distribution .> δ) - - # Extract the elements of the invariant measure corresponding to these indices - return InvariantMeasure(to, Probabilities(distribution)) -end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/SimplexPoint-checkpoint.jl b/src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/SimplexPoint-checkpoint.jl deleted file mode 100644 index 2ac7139a9..000000000 --- a/src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/SimplexPoint-checkpoint.jl +++ /dev/null @@ -1,131 +0,0 @@ -export SimplexPoint, transferoperator, SimplexPoint -import StaticArrays: SizedVector, SizedMatrix, MMatrix - -include("helper_functions.jl") - - -""" - SimplexPoint() <: Entropies.BinningProbabilitiesEstimator - -A probability estimator based on the transfer operator, using a triangulated partition -and approximate simplex intersections to compute transition probabilities [^Diego2019]. - -To use this estimator, the Simplices.jl package must be brought into scope by doing `using Simplices` after -`using PerronFrobenius`. - -[^Diego2019]: Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212. -""" -struct SimplexPoint <: BinningProbabilitiesEstimator - bc::String - - function SimplexPoint(bc::String = "circular") - isboundarycondition(bc, "triangulation") || error("Boundary condition '$bc' not valid.") - new(bc) - end -end -Base.show(io::IO, se::SimplexPoint) = print(io, "SimplexPoint{$(se.bc)}") - - -""" Generate a transopergenerator for an exact simplex estimator.""" -function transopergenerator(pts, method::TransferOperator{<:SimplexPoint}) - # modified points, where the image of each point is guaranteed to lie within the convex hull of the previous points - invariant_pts = invariantize(pts) - - # triangulation of the invariant points. The last point is excluded, so that - # the last vertex also can be mapped forward one step in time. - triang = DelaunayTriangulation(invariant_pts[1:end-1]) - - init = (invariant_pts = invariant_pts, triang = triang) - - TransferOperatorGenerator(method, pts, init) -end - - - -""" -n: The number of points sampled inside each simplex. -randomsampling: Sample randomly (`random == true`), or using an even simplex splitting routine (`random == false`) -""" -function (tog::TransferOperatorGenerator{T})(; tol = 1e-8, randomsampling::Bool = false, - n::Int = 100) where T <: TransferOperator{SimplexPoint} - pts, inds = getfield.(Ref(tog.init), (:invariant_pts, :triang)) - D = length(pts[1]) - N = length(inds) - ϵ = tol / N - - # Prepare memory-efficient representations of the simplices - simplices = reshape_simplices(pts[1:(end- 1)], inds) - image_simplices = reshape_simplices(pts[2:end], inds) - - #= - Compute the convex expansions coefficients needed to generate points - within a simplex. Also, update number of points if a shape-preserving - simplex subdivision algorithm was used (in that case, we get a - *minimum* of `n_sample_pts` set of coefficients. - =# - convex_coeffs = subsample_coeffs(D, n, randomsampling) - n_coeffs::Int = size(convex_coeffs, 2) - convex_coeffs = [convex_coeffs[:, i] for i = 1:n_coeffs] - - # Pre-allocated arrays (SizedArrays, for efficiency) - L = (D+1)^2 - NVERTICES = D+1 - pt = SizedMatrix{1, D}(zeros(1, D)) - s_arr = SizedVector{L}(zeros(L)) - signs = SizedVector{NVERTICES}(zeros(NVERTICES)) - - # Re-arrange simplices so that look-up is a bit more efficient - simplex_arrs = Vector{Array{Float64, 2}}(undef, N) - imsimplex_arrs = Vector{Array{Float64, 2}}(undef, N) - for i in 1:N - simplex_arrs[i] = hcat(pts[inds[i]]...,) - imsimplex_arrs[i] = hcat(pts[inds[i] .+ 1]...,) - end - - # Pre-allocate some information to avoid memory allocation when - # computing the fractional overlaps. These are a temporary matrix - # which we use to compute determinants (to check if a point lies - # inside a simplex), and the indices needed to access its values - # at different stages of the checking stage (of whether a point - # lies inside a simplex). Pre-allocating this information and - # passing it on reduces runtime and memory allocation - # by half an order of magnitude. - temp_arr = rand(MMatrix{D + 1, D + 1}) - starts = [((D + 1)*i)+1 for i = 0:D] - stops = starts .+ D - - - # The Markov matrix - M = zeros(Float64, N, N) - - for i in 1:N - # Find the indices of the simplices who potentially intersect - # with the image of simplex #i - idxs = idxs_potentially_intersecting_simplices(pts, inds, i) - Sj = Vector{AbstractArray}(undef, length(idxs)) - get_simplices_at_inds!(Sj, idxs, simplices) - - # Generate points within the image of simplex #i and check - # what fraction of those points fall into the simplices. - # We approximate the intersecting volume between the - # image simplex and the individual simplices by those - # fractions. - @views is = imsimplex_arrs[i] - for k in 1:n_coeffs - pt = transpose(convex_coeffs[k]) * transpose(is) - innerloop_optim!(idxs, signs, s_arr, Sj, pt, D, M, i, temp_arr, starts, stops) - end - end - # Need to normalise, because all we have up until now is counts - # of how many points inside the image simplex falls into - # the simplices. - params = (tol = tol, randomsampling = randomsampling, n = n) - M_normalized = transpose(M) ./ n_coeffs - return TransferOperatorApproximation(tog, M_normalized, params) - -end - -function transferoperator(pts, method::TransferOperator{<:SimplexPoint}; tol = 1e-8, randomsampling::Bool = false, n::Int = 100) - tog = transopergenerator(pts, method) - tog(tol = tol, randomsampling = randomsampling, n = n) -end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/helper_functions-checkpoint.jl b/src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/helper_functions-checkpoint.jl deleted file mode 100644 index b49cd0a68..000000000 --- a/src/binning_based/transferoperator/triangular/point/.ipynb_checkpoints/helper_functions-checkpoint.jl +++ /dev/null @@ -1,175 +0,0 @@ -import LinearAlgebra: det -import StaticArrays: SizedMatrix -import .DelaunayTriangulations: DelaunayTriangulation - -""" Re-zeros the array `a`. """ -function rezero!(a) - @inbounds for i in eachindex(a) - a[i] = 0.0 - end -end - -""" Fill the elements of vector `from` into vector `into`. """ -function fill_into!(into, from) - @inbounds for i in eachindex(into) - into[i] = from[i] - end -end - -""" Fill vector `into` at indices `inds` with the elements at `from[inds]`. """ -function fill_at_inds!(into, from, inds) - @inbounds for i in 1:length(inds) - into[inds[i]] = from[i] - end -end - -function mdet_fm_in_x_view!(M, v, n, starts, stops) - @inbounds for i = 0:(n - 1) - M[starts[i+1]:stops[i+1]] = view(v, starts[i+1]:stops[i+1]) - end - @fastmath det(M) -end - - -function iscontained!(signs, s_arr, sx, point, dim, temp_arr, starts, stops) - # Redefine the temporary simplex. This is in-place, so we don't allocate - # memory. We could also have re-initialised `signs`, but since we're never - # comparing more than two consecutive signs, this is not necessary. - rezero!(s_arr) - rezero!(signs) - fill_into!(s_arr, sx) - - #Replace first vertex with the point - fill_at_inds!(s_arr, point, 1:dim) - - # Signed volume - signs[1] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops)) - - rezero!(s_arr) - fill_into!(s_arr, sx) #reset - - for κ = 2:dim # Check remaining signs and stop if sign changes - # Replace the ith vertex with the point we're cheking (leaving the - # 1 appended to Vi intact.) - idxs = ((dim + 1)*(κ - 1)+1):((dim + 1)*(κ - 1)+ 1 + dim - 1) - fill_at_inds!(s_arr, point, idxs) # ith change - - signs[κ] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops)) - - if !(signs[κ-1] == signs[κ]) - return false - end - - rezero!(s_arr) - fill_into!(s_arr, sx) - end - - # Last the last vertex with the point in question - idxs = ((dim + 1)*(dim)+1):((dim+1)^2-1) - fill_at_inds!(s_arr, point, idxs) - - signs[end] = sign(mdet_fm_in_x_view!(temp_arr, s_arr, dim + 1, starts, stops)) - - if !(signs[end-1] == signs[end]) - return false - else - return true - end -end - -function innerloop_optim!(inds::Vector{Int}, signs, s_arr, Sj, pt, dim::Int, M, i::Int, temp_arr, starts, stops) - for j in 1:length(inds) - if iscontained!(signs, s_arr, Sj[j], pt, dim, temp_arr, starts, stops) - M[inds[j], i] += 1.0 - end - end -end - -""" - reshape_simplices(pts::Vector{Vector{T}}, DT::DelaunayTriangulation) where T - -Creates alternative representations of the simplices that allow for efficient -(mostly non-allocating) checking if a point lies inside a simplex. - -This function adds some information needed when calculating the orientations -of the simplices *before* we start the computation (appends a row of ones on top of the -matrix representation of the simplex). Otherwise, the fallback is machinery in Simplices.jl, -which is much slower because the additional information must be added -`(n_simplices^2)*n_sample_pts` times (creating that many new matrices), which takes forever. - -Reshaping the simplices and appending that information beforehand avoids excessive -memory allocation. -""" -function reshape_simplices(pts, DT::DelaunayTriangulation) - n_simplices = length(DT.indices) - n_vertices = length(DT[1]) - dim = n_vertices - 1 - - # The [:, j, i] th entry of these two arrays holds the jth vertex of the - # ith simplex, but instead of having just `dim` vertices, we append a `1` - # to the end of the vectors. This allows for efficent (non-allocating) - # computation within the `contains_point_lessalloc!` function. If we instead - # would have appended the 1's inside that function, we would be performing - # memory-allocating operations, which are very expensive. Doing this instead - # gives orders of magnitude speed-ups for sufficiently large triangulations. - S1 = Array{Float64}(undef, dim + 1, dim + 1, n_simplices) - - # Collect simplices in the form of (dim+1)^2-length column vectors. This - # also helps with the - NDIM, N_SIMPLICES = (dim+1)^2, n_simplices - simplices = SizedMatrix{NDIM}{N_SIMPLICES}(zeros((NDIM, N_SIMPLICES))) - - @inbounds for i in 1:n_simplices - for j in 1:n_vertices - S1[:, j, i] = vcat(pts[DT[i][j]], 1.0) - end - - simplices[:, i] = reshape(S1[:, :, i], (dim+1)^2) - end - - return simplices -end - - -""" -idxs_potentially_intersecting_simplices(all_pts::Vector{Vector{T}}, - DT::DelaunayTriangulation, idx::Int) where T - -Given a delaunay triangulation `DT` of `all_pts[1:(end - 1)]`, find the indices of -the simplices that potentially intersect with the image simplex with index `idx`. -""" -function idxs_potentially_intersecting_simplices(all_pts, DT::DelaunayTriangulation, idx::Int) - # Vector that will hold the indices of the simplices potentially - # intersecting with the image of simplex #idx - inds_potential_simplices = Int[] - - n_simplices = length(DT) - - original_pts = all_pts[1:(end - 1)] - forward_pts = all_pts[2:end] - - simplices = [MutableSimplex(original_pts[DT[i]]) for i = 1:n_simplices] - image_simplices = [MutableSimplex(forward_pts[DT[i]]) for i = 1:n_simplices] - - cs = [centroid(simplices[i]) for i = 1:n_simplices] - cs_im = [centroid(image_simplices[i]) for i = 1:n_simplices] - - rs = [radius(simplices[i]) for i = 1:n_simplices] - rs_im = [radius(image_simplices[i]) for i = 1:n_simplices] - - @inbounds for i = 1:n_simplices - δ = transpose(cs_im[idx] .- cs[i]) * (cs_im[idx] .- cs[i]) .- ((rs_im[idx] + rs[i])^2) - if δ[1] < 0 - push!(inds_potential_simplices, i) - end - end - - return inds_potential_simplices -end - -""" Get the simplices at the given indices. """ -function get_simplices_at_inds!(simps, inds::Vector{Int}, simplices) - for i in 1:length(inds) - simps[i] = simplices[:, inds[i]] - end -end \ No newline at end of file From 3dc85077147814967cf4ca8e0989ceac7439c454 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 13:46:03 +0100 Subject: [PATCH 04/29] Use generators for transfer operator approximations --- src/Entropies.jl | 3 +- src/binning_based/transferoperator/common.jl | 39 ++++-- .../rectangular/rectangular.jl | 115 ++++++++++-------- .../triangular/exact/SimplexExact.jl | 73 +---------- .../triangular/point/SimplexPoint.jl | 4 +- .../triangular/point/helper_functions.jl | 4 - .../transferoperator/triangular/triangular.jl | 6 + test/runtests.jl | 34 +++++- 8 files changed, 137 insertions(+), 141 deletions(-) create mode 100644 src/binning_based/transferoperator/triangular/triangular.jl diff --git a/src/Entropies.jl b/src/Entropies.jl index 0855ab627..8c41acf13 100644 --- a/src/Entropies.jl +++ b/src/Entropies.jl @@ -14,8 +14,7 @@ module Entropies function __init__() @require Simplices="d5428e67-3037-59ba-9ab1-57a04f0a3b6a" begin using .Simplices - include("binning_based/transferoperator/triangular/exact/SimplexExact.jl") - include("binning_based/transferoperator/triangular/point/SimplexPoint.jl") + include("binning_based/transferoperator/triangular/triangular.jl") end end end diff --git a/src/binning_based/transferoperator/common.jl b/src/binning_based/transferoperator/common.jl index 7a17de2fb..995bcec00 100644 --- a/src/binning_based/transferoperator/common.jl +++ b/src/binning_based/transferoperator/common.jl @@ -4,7 +4,8 @@ import LinearAlgebra: norm export transferoperator, - TransferOperator, + TransferOperator, + TransferOperatorApproximation, invariantmeasure, InvariantMeasure, transfermatrix @@ -92,8 +93,8 @@ struct TransferOperatorGenerator{E <: TransferOperator, X, A} end -function Base.show(io::IO, DT::TransferOperatorGenerator{E, X, A}) where {E <: TransferOperator{R}, X, A} where R - summary = "Transfer approximation from $R estimator" +function Base.show(io::IO, DT::TransferOperatorGenerator{E, X, A}) where {E, X, A} + summary = "TransferOperatorGenerator{method: $E, pts: $X, init: $A}" println(io, summary) end @@ -133,15 +134,21 @@ bin information). See also: [`invariantmeasure`](@ref). """ -struct InvariantMeasure{T} +struct InvariantMeasure{T, P<:Probabilities} to::T - ρ::Probabilities + ρ::P - function InvariantMeasure(to::T, ρ) where T - new{T}(to, ρ) + function InvariantMeasure(to::T, ρ::P) where {T, P} + new{T, P}(to, ρ) end end + +function Base.show(io::IO, DT::InvariantMeasure{T, P}) where {T, P} + summary = "InvariantMeasure{transfer operator approximation: $T, probabilities: $P}" + println(io, summary) +end + """ TransferOperatorApproximation(generator, transfermatrix, params) @@ -181,8 +188,8 @@ struct TransferOperatorApproximation{G<:TransferOperator, T} params end -function Base.show(io::IO, DT::TransferOperatorApproximation{G, T}) where {G <: TransferOperator, T} - summary = "TransferOperatorApproximation{$G, $T}" +function Base.show(io::IO, DT::TransferOperatorApproximation{G, T}) where {G, T} + summary = "TransferOperatorApproximation{generator: $G, transfer matrix: $T}" println(io, summary) end @@ -322,3 +329,17 @@ function invariantmeasure(to::TransferOperatorApproximation; # Extract the elements of the invariant measure corresponding to these indices return InvariantMeasure(to, Probabilities(distribution)) end + + +function invariantmeasure(pts, method::TransferOperator{<:R}; kwargs...) where R + to_approximation = transferoperator(pts, method; kwargs...) + invariantmeasure(to_approximation) +end + +probabilities(iv::InvariantMeasure) = iv.ρ +probabilities(iv::TransferOperatorApproximation) = probabilities(invariantmeasure(iv.ρ)) + +function probabilities(pts, method::TransferOperator{<:R}; kwargs...) where R + to_approximation = transferoperator(pts, method; kwargs...) + invariantmeasure(to_approximation).ρ +end \ No newline at end of file diff --git a/src/binning_based/transferoperator/rectangular/rectangular.jl b/src/binning_based/transferoperator/rectangular/rectangular.jl index 298bf8d68..f7b086f9c 100644 --- a/src/binning_based/transferoperator/rectangular/rectangular.jl +++ b/src/binning_based/transferoperator/rectangular/rectangular.jl @@ -53,37 +53,13 @@ end inds_in_terms_of_unique(x::AbstractDataset) = inds_in_terms_of_unique(x.data) -# TODO: transferoperatorgenerator for rectangualrbinning! - -""" - transferoperator(pts::AbstractDataset{D, T}, ϵ::RectangularBinning) → TransferOperatorApproximationRectangular - -Estimate the transfer operator given a set of sequentially ordered points subject to a -rectangular partition given by `ϵ`. - -## Example - -```julia -using DynamicalSystems, Plots, Entropy -D = 4 -ds = Systems.lorenz96(D; F = 32.0) -N, dt = 20000, 0.1 -orbit = trajectory(ds, N*dt; dt = dt, Ttr = 10.0) - -# Estimate transfer operator over some coarse graining of the orbit. -transferoperator(orbit, RectangularBinning(10)) -``` - -See also: [`RectangularBinning`](@ref). -""" -function transferoperator(pts::AbstractDataset{D, T}, to::TransferOperator{<:RectangularBinning}; - boundary_condition = :circular) where {D, T<:Real} - +function transopergenerator(pts, method::TransferOperator{<:RectangularBinning}) + L = length(pts) - mini, edgelengths = Entropies.minima_edgelengths(pts, to.ϵ) + mini, edgelengths = minima_edgelengths(pts, method.ϵ) # The L points visits a total of L bins, which are the following bins: - visited_bins = Entropies.encode_as_bin(pts, mini, edgelengths) + visited_bins = encode_as_bin(pts, mini, edgelengths) sort_idxs = sortperm(visited_bins) # TODO: fix re-indexing after sorting. Sorting is much faster, so we want to do so. @@ -99,7 +75,58 @@ function transferoperator(pts::AbstractDataset{D, T}, to::TransferOperator{<:Rec # first_visited_by == [x[1] for x in visitors] first_visited_by = GroupSlices.firstinds(slices) - L = length(first_visited_by) + + init = (mini = mini, + edgelengths = edgelengths, + visited_bins = visited_bins, + sort_idxs = sort_idxs, + visits_whichbin = visits_whichbin, + visitors = visitors, + first_visited_by = first_visited_by + ) + + TransferOperatorGenerator(method, pts, init) +end + +# """ +# transferoperator(pts::AbstractDataset{D, T}, ϵ::RectangularBinning) → TransferOperatorApproximationRectangular + +# Estimate the transfer operator given a set of sequentially ordered points subject to a +# rectangular partition given by `ϵ`. + +# ## Example + +# ```julia +# using DynamicalSystems, Plots, Entropy +# D = 4 +# ds = Systems.lorenz96(D; F = 32.0) +# N, dt = 20000, 0.1 +# orbit = trajectory(ds, N*dt; dt = dt, Ttr = 10.0) + +# # Estimate transfer operator over some coarse graining of the orbit. +# transferoperator(orbit, RectangularBinning(10)) +# ``` + +# See also: [`RectangularBinning`](@ref). +# """ +function (tog::TransferOperatorGenerator{T})(;boundary_condition = :circular, + ) where T <: TransferOperator{<:RectangularBinning}; + + visitors, + first_visited_by, + visits_whichbin, + visited_bins, + edgelengths, + mini = getfield.( + Ref(tog.init), + (:visitors, + :first_visited_by, + :visits_whichbin, + :visited_bins, + :edgelengths, + :mini + ), + ) I = Int32[] J = Int32[] @@ -111,17 +138,15 @@ function transferoperator(pts::AbstractDataset{D, T}, to::TransferOperator{<:Rec n_visitsᵢ::Int = 0 if boundary_condition == :circular - #warn("Using circular boundary condition") append!(visits_whichbin, [1]) elseif boundary_condition == :random - #warn("Using random circular boundary condition") append!(visits_whichbin, [rand(1:length(visits_whichbin))]) else error("Boundary condition $(boundary_condition) not implemented") end # Loop over the visited bins bᵢ - for i in 1:L + for i in 1:length(first_visited_by) # How many times is this bin visited? n_visitsᵢ = length(visitors[i]) @@ -187,24 +212,14 @@ function transferoperator(pts::AbstractDataset{D, T}, to::TransferOperator{<:Rec # row/column of the transfer operator unique!(visited_bins) bins = [β .* edgelengths .+ mini for β in visited_bins] - params = (mini = mini, edgelengths = edgelengths, bins = bins, - sort_idxs = sort_idxs, visitors) + params = (bins = bins) - TransferOperatorApproximation(to, M, params) + TransferOperatorApproximation(tog, M, params) end -function invariantmeasure(iv::InvariantMeasure) - return iv.ρ, iv.to.bins -end - -function invariantmeasure(x::AbstractDataset, ϵ::TransferOperator{<:RectangularBinning}) - to = transferoperator(x, ϵ) - invariantmeasure(to) -end - -function probabilities(x::AbstractDataset, est::TransferOperator{<:RectangularBinning}) - to = transferoperator(x, est.ϵ) - iv = invariantmeasure(to) - - return iv.ρ -end +function transferoperator(pts, method::TransferOperator{<:RectangularBinning}; + boundary_condition = :circular, + ) + tog = transopergenerator(pts, method) + tog(boundary_condition = boundary_condition) +end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl index 659284593..1956cbf49 100644 --- a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl +++ b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl @@ -1,7 +1,4 @@ export SimplexExact, invariantmeasure -import Entropies.invariantmeasure - - """ SimplexExact @@ -81,75 +78,7 @@ function (tog::TransferOperatorGenerator{T})(; tol = 1e-8) where T <: TransferOp TransferOperatorApproximation(tog, M, params) end -function transferoperator(pts, method::SimplexExact; tol::Real = 1e-8) +function transferoperator(pts, method::TransferOperator{<:SimplexExact}; tol::Real = 1e-8) tog = transopergenerator(pts, method) tog(tol = tol) end - - -function invariantmeasure(to::TransferOperatorApproximation{G, T}; - N::Int = 200, - tolerance::Float64 = 1e-8, - delta::Float64 = 1e-8) where {G <: TransferOperator{R}, T} where R <: SimplexExact - - TO = to.transfermatrix - #= - # Start with a random distribution `Ρ` (big rho). Normalise it so that it - # sums to 1 and forms a true probability distribution over the partition elements. - =# - Ρ = rand(Float64, 1, size(to.transfermatrix, 1)) - Ρ = Ρ ./ sum(Ρ, dims = 2) - - #= - # Start estimating the invariant distribution. We could either do this by - # finding the left-eigenvector of M, or by repeated application of M on Ρ - # until the distribution converges. Here, we use the latter approach, - # meaning that we iterate until Ρ doesn't change substantially between - # iterations. - =# - distribution = Ρ * TO - - distance = norm(distribution - Ρ) / norm(Ρ) - - check = floor(Int, 1 / delta) - check_pts = floor.(Int, transpose(collect(1:N)) ./ check) .* transpose(collect(1:N)) - check_pts = check_pts[check_pts .> 0] - num_checkpts = size(check_pts, 1) - check_pts_counter = 1 - - counter = 1 - while counter <= N && distance >= tolerance - counter += 1 - Ρ = distribution - - # Apply the Markov matrix to the current state of the distribution - distribution = Ρ * to.transfermatrix - - if (check_pts_counter <= num_checkpts && - counter == check_pts[check_pts_counter]) - - check_pts_counter += 1 - colsum_distribution = sum(distribution, dims = 2)[1] - if abs(colsum_distribution - 1) > delta - distribution = distribution ./ colsum_distribution - end - end - - distance = norm(distribution - Ρ) / norm(Ρ) - end - distribution = dropdims(distribution, dims = 1) - - # Do the last normalisation and check - colsum_distribution = sum(distribution) - - if abs(colsum_distribution - 1) > delta - distribution = distribution ./ colsum_distribution - end - - # Find partition elements with strictly positive measure. - δ = tolerance/size(to.transfermatrix, 1) - inds_nonzero = findall(distribution .> δ) - - # Extract the elements of the invariant measure corresponding to these indices - return InvariantMeasure(to, Probabilities(distribution)) -end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl b/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl index 616d69488..ec375af39 100644 --- a/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl +++ b/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl @@ -124,7 +124,9 @@ function (tog::TransferOperatorGenerator{T})(; tol = 1e-8, randomsampling::Bool end -function transferoperator(pts, method::TransferOperator{<:SimplexPoint}; tol = 1e-8, randomsampling::Bool = false, n::Int = 100) +function transferoperator(pts, method::TransferOperator{<:SimplexPoint}; + tol = 1e-8, randomsampling::Bool = false, n::Int = 100,) + tog = transopergenerator(pts, method) tog(tol = tol, randomsampling = randomsampling, n = n) end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/point/helper_functions.jl b/src/binning_based/transferoperator/triangular/point/helper_functions.jl index 12ff94fba..4e2ad08b5 100644 --- a/src/binning_based/transferoperator/triangular/point/helper_functions.jl +++ b/src/binning_based/transferoperator/triangular/point/helper_functions.jl @@ -1,10 +1,6 @@ import LinearAlgebra: det import StaticArrays: SizedMatrix -include("../delaunay_triangulations/DelaunayTriangulations.jl") -import .DelaunayTriangulations: DelaunayTriangulation - - """ Re-zeros the array `a`. """ function rezero!(a) @inbounds for i in eachindex(a) diff --git a/src/binning_based/transferoperator/triangular/triangular.jl b/src/binning_based/transferoperator/triangular/triangular.jl new file mode 100644 index 000000000..2e7f10531 --- /dev/null +++ b/src/binning_based/transferoperator/triangular/triangular.jl @@ -0,0 +1,6 @@ +include("simplex_types/simplex_types.jl") +include("./delaunay_triangulations/DelaunayTriangulations.jl") +import .DelaunayTriangulations: DelaunayTriangulation + +include("exact/SimplexExact.jl") +include("point/SimplexPoint.jl") \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 33a30007c..ea6fbf8d7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -52,7 +52,13 @@ end @test SymbolicAmplitudeAwarePermutation(lt = Entropies.isless_rand) isa SymbolicAmplitudeAwarePermutation @test VisitationFrequency(RectangularBinning(3)) isa VisitationFrequency + + # Transfer operator related @test TransferOperator(RectangularBinning(3)) isa TransferOperator + @test TransferOperator(SimplexPoint()) isa TransferOperator + @test TransferOperator(SimplexExact()) isa TransferOperator + + @test TimeScaleMODWT() isa TimeScaleMODWT @test TimeScaleMODWT(Wavelets.WT.Daubechies{8}()) isa TimeScaleMODWT @test Kraskov(k = 2, w = 1) isa Kraskov @@ -334,12 +340,12 @@ end RectangularBinning(3), RectangularBinning(0.2), RectangularBinning([2, 2, 3]), - RectangularBinning([0.2, 0.3, 0.3]) + RectangularBinning([0.2, 0.3, 0.3]), ] - @testset "Binning test $i" for i in 1:length(binnings) + @testset "Rectagular binning test $i" for i in 1:length(binnings) to = Entropies.transferoperator(D, binnings[i]) - @test to isa Entropies.TransferOperatorApproximationRectangular + @test to isa Entropies.TransferOperatorApproximation iv = invariantmeasure(to) @test iv isa InvariantMeasure @@ -350,5 +356,27 @@ end @test probabilities(D, TransferOperator(binnings[i])) isa Probabilities end + + @testset "Triangulation binning" begin + D = Dataset(rand(18, 3)) + + est_point = TransferOperator(SimplexPoint()) + est_exact = TransferOperator(SimplexExact()) + + to_point = transferoperator(D, est_point) + to_exact = transferoperator(D, est_exact) + @test to_point isa Entropies.TransferOperatorApproximation + @test to_exact isa Entropies.TransferOperatorApproximation + + @test invariantmeasure(to_point) isa Entropies.InvariantMeasure + @test invariantmeasure(to_exact) isa Entropies.InvariantMeasure + @test probabilities(to_point) isa Probabilities + @test probabilities(to_exact) isa Probabilities + + @test invariantmeasure(D, est_point) isa InvariantMeasure + @test invariantmeasure(D, est_exact) isa InvariantMeasure + @test probabilities(D, est_point) isa Probabilities + @test probabilities(D, est_exact) isa Probabilities + end end end From 4b8e6e528124cf35d2ad30b7d48878eb75fd4b1b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 14:23:57 +0100 Subject: [PATCH 05/29] Remove old stuff --- src/binning_based/binningbased.jl | 3 +- src/binning_based/transferoperator.jl | 131 -------------------------- 2 files changed, 1 insertion(+), 133 deletions(-) delete mode 100644 src/binning_based/transferoperator.jl diff --git a/src/binning_based/binningbased.jl b/src/binning_based/binningbased.jl index e600b3fbd..99c845b5d 100644 --- a/src/binning_based/binningbased.jl +++ b/src/binning_based/binningbased.jl @@ -1,3 +1,2 @@ abstract type BinningProbabilitiesEstimator <: ProbabilitiesEstimator end -include("rectangular/rectangular_estimators.jl") -#include("triangular/TriangulationEstimators.jl") \ No newline at end of file +include("rectangular/rectangular_estimators.jl") \ No newline at end of file diff --git a/src/binning_based/transferoperator.jl b/src/binning_based/transferoperator.jl deleted file mode 100644 index 8c56a3ebd..000000000 --- a/src/binning_based/transferoperator.jl +++ /dev/null @@ -1,131 +0,0 @@ -using DelayEmbeddings, SparseArrays - -export - TransferOperator, # the probabilities estimator - InvariantMeasure, - transfermatrix, - invariantmeasure - -""" - TransferOperator(ϵ::RectangularBinning) <: BinningProbabilitiesEstimator - -A probability estimator based on binning data into rectangular boxes dictated by -the binning scheme `ϵ`, then approxmating the transfer (Perron-Frobenius) operator -over the bins, the taking the invariant measure associated with that transfer operator -as the bin probabilities. Assumes that the input data are sequential. - -This implementation follows the grid estimator approach in Diego et al. (2019)[^Diego2019]. - -## Description - -The transfer operator ``P^{N}``is computed as an `N`-by-`N` matrix of transition -probabilities between the states defined by the partition elements, where `N` is the -number of boxes in the partition that is visited by the orbit/points. - -If ``\\{x_t^{(D)} \\}_{n=1}^L`` are the ``L`` different ``D``-dimensional points over -which the transfer operator is approximated, ``\\{ C_{k=1}^N \\}`` are the ``N`` different -partition elements (as dictated by `ϵ`) that gets visited by the points, and - ``\\phi(x_t) = x_{t+1}``, then - -```math -P_{ij} = \\dfrac -{\\#\\{ x_n | \\phi(x_n) \\in C_j \\cap x_n \\in C_i \\}} -{\\#\\{ x_m | x_m \\in C_i \\}}, -``` - -where ``\\#`` denotes the cardinal. The element ``P_{ij}`` thus indicates how many points -that are initially in box ``C_i`` end up in box ``C_j`` when the points in ``C_i`` are -projected one step forward in time. Thus, the row ``P_{ik}^N`` where -``k \\in \\{1, 2, \\ldots, N \\}`` gives the probability -of jumping from the state defined by box ``C_i`` to any of the other ``N`` states. It -follows that ``\\sum_{k=1}^{N} P_{ik} = 1`` for all ``i``. Thus, ``P^N`` is a row/right -stochastic matrix. - -### Invariant measure estimation from transfer operator - -The left invariant distribution ``\\mathbf{\\rho}^N`` is a row vector, where -``\\mathbf{\\rho}^N P^{N} = \\mathbf{\\rho}^N``. Hence, ``\\mathbf{\\rho}^N`` is a row -eigenvector of the transfer matrix ``P^{N}`` associated with eigenvalue 1. The distribution -``\\mathbf{\\rho}^N`` approximates the invariant density of the system subject to the -partition `ϵ`, and can be taken as a probability distribution over the partition elements. - -In practice, the invariant measure ``\\mathbf{\\rho}^N`` is computed using -[`invariantmeasure`](@ref), which also approximates the transfer matrix. The invariant distribution -is initialized as a length-`N` random distribution which is then applied to ``P^{N}``. -The resulting length-`N` distribution is then applied to ``P^{N}`` again. This process -repeats until the difference between the distributions over consecutive iterations is -below some threshold. - -## Probability and entropy estimation - -- `probabilities(x::AbstractDataset, est::TransferOperator{RectangularBinning})` estimates - probabilities for the bins defined by the provided binning (`est.ϵ`) -- `genentropy(x::AbstractDataset, est::TransferOperator{RectangularBinning})` does the same, - but computes generalized entropy using the probabilities. - - -See also: [`RectangularBinning`](@ref), [`invariantmeasure`](@ref). - -[^Diego2019]: Diego, D., Haaga, K. A., & Hannisdal, B. (2019). Transfer entropy computation using the Perron-Frobenius operator. Physical Review E, 99(4), 042212. -""" -struct TransferOperator{R} <: BinningProbabilitiesEstimator - ϵ::R - - function TransferOperator(ϵ::R) where R #<: RectangularBinning - new{R}(ϵ) - end -end -struct TransferOperatorGenerator{E <: TransferOperator, X, A} - method::E # estimator with its input parameters - pts::X # the phase space / reconstruted state space points - init::A # pre-initialized things that speed up estimation process -end - - -""" - transopergenerator(pts, method::TransferOperator) → to::TransferOperatorGenerator - -Initialize a generator that creates transfer operators on demand, based on the given `method`. -This is efficient, because some things can be initialized and reused. - -To approximate a transfer operator, call `to` as a function with the relevant arguments. - -```julia -to = transopergenerator(x, TransferOperator(RectangularBinning(5))) -for i in 1:1000 - s = to() - # do stuff with s and or x - result[i] = stuff -end -``` -""" -function transopergenerator end - -function transferoperator end - -function invariantmeasure end - -function transfermatrix end - -function transferoperator(pts, method::TransferOperator) - to = transopergenerator(pts, method) - to() -end - -""" - InvariantMeasure(to, ρ) - -Minimal return struct for [`invariantmeasure`](@ref) that contains the estimated invariant -measure `ρ`, as well as the transfer operator `to` from which it is computed (including -bin information). - -See also: [`invariantmeasure`](@ref). -""" -struct InvariantMeasure{T} - to::T - ρ::Probabilities - - function InvariantMeasure(to::T, ρ) where T - new{T}(to, ρ) - end -end From 99bc594e3c57656d549e25aec91de7d14db68f72 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 14:24:10 +0100 Subject: [PATCH 06/29] Update docs --- src/binning_based/transferoperator/common.jl | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/src/binning_based/transferoperator/common.jl b/src/binning_based/transferoperator/common.jl index 995bcec00..4e487718d 100644 --- a/src/binning_based/transferoperator/common.jl +++ b/src/binning_based/transferoperator/common.jl @@ -74,7 +74,24 @@ below some threshold. See also: [`RectangularBinning`](@ref), [`SimplexPoint`](@ref), [`SimplexExact`](@ref), -[`invariantmeasure`](@ref), +[`invariantmeasure`](@ref) + +## Examples + +Here, we create three different transfer operator-based estimators. + +```@example +# A rectangular binning is suited for datasets with a large number of points +est_rect = TransferOperator(RectangularBinning(5)) + +# A triangulated binning, using approximate simplex intersections, is also possible for +# datasets with not too many points (say, <1000 points). +est_point = TransferOperator(SimplexPoint()) + +# For datasets with few points, say <100 points, exact simplex intersections may also +# be computationally feasible. +est_exact = TransferOperator(SimplexExact()) +``` [^Diego2019]: Diego, D., Haaga, K. A., & Hannisdal, B. (2019). Transfer entropy computation using the Perron-Frobenius operator. Physical Review E, 99(4), 042212. """ From 886f9e93e2770c699b194ebf066a57be75252da4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 14:24:19 +0100 Subject: [PATCH 07/29] Add convenience methods --- src/binning_based/transferoperator/common.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/binning_based/transferoperator/common.jl b/src/binning_based/transferoperator/common.jl index 4e487718d..12824913c 100644 --- a/src/binning_based/transferoperator/common.jl +++ b/src/binning_based/transferoperator/common.jl @@ -354,7 +354,7 @@ function invariantmeasure(pts, method::TransferOperator{<:R}; kwargs...) where R end probabilities(iv::InvariantMeasure) = iv.ρ -probabilities(iv::TransferOperatorApproximation) = probabilities(invariantmeasure(iv.ρ)) +probabilities(iv::TransferOperatorApproximation) = probabilities(invariantmeasure(iv)) function probabilities(pts, method::TransferOperator{<:R}; kwargs...) where R to_approximation = transferoperator(pts, method; kwargs...) From a89aedda9e0211cbfa5beec39f53133a66d96cbb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 14:24:26 +0100 Subject: [PATCH 08/29] Fix tests --- test/runtests.jl | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index ea6fbf8d7..c8da185d4 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Entropies using DelayEmbeddings using Wavelets using StaticArrays +using Simplices @testset "Histogram estimation" begin x = rand(1:10, 100) @@ -344,16 +345,12 @@ end ] @testset "Rectagular binning test $i" for i in 1:length(binnings) - to = Entropies.transferoperator(D, binnings[i]) + est = TransferOperator(binnings[i]) + to = Entropies.transferoperator(D, est) @test to isa Entropies.TransferOperatorApproximation iv = invariantmeasure(to) @test iv isa InvariantMeasure - - p, bins = invariantmeasure(iv) - @test p isa Probabilities - @test bins isa Vector{<:SVector} - @test probabilities(D, TransferOperator(binnings[i])) isa Probabilities end From 06717825e166ef6d955c245284de5d52eea0eebc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 14:24:42 +0100 Subject: [PATCH 09/29] Update project.toml --- Project.toml | 21 +++++++-------------- 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/Project.toml b/Project.toml index bb1102d8b..64eb3f1ed 100644 --- a/Project.toml +++ b/Project.toml @@ -8,34 +8,27 @@ version = "0.11.1" DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" -Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" -DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" -PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" Requires = "ae029012-a4dd-5104-9daa-d747884805df" -Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" -StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" [compat] -DelayEmbeddings = "1.10" -Distances = "0.9, 0.10" -NearestNeighbors = "0.4" -Requires = "^1.1" -SpecialFunctions = "0.10, 1.0" -StaticArrays = "0.12, 1.0" -Wavelets = "0.9" +DelayEmbeddings = "^1.20" +Distances = "0.9, ^0.10" +NearestNeighbors = "^0.4" +SpecialFunctions = "^0.10, ^1" +StaticArrays = "^0.12, ^1" +Wavelets = "^0.9" julia = "1.5" [extras] -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["Test", "Simplices"] From 8ee148b290cda39101e58e6edb61da94450fec9d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 14:25:07 +0100 Subject: [PATCH 10/29] Set version to 1.0.0 --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 64eb3f1ed..7468e4fa3 100644 --- a/Project.toml +++ b/Project.toml @@ -2,7 +2,7 @@ name = "Entropies" uuid = "ed8fcbec-b94c-44b6-89df-898894ad9591" authors = "Kristian Agasøster Haaga , George Datseries " repo = "https://github.com/juliadynamics/Entropies.jl.git" -version = "0.11.1" +version = "1.0.0" [deps] DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f" From d84bd68d93c843df5efc241baa33f78ceacbe0fb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 14:47:06 +0100 Subject: [PATCH 11/29] Update docs --- src/binning_based/transferoperator/common.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/binning_based/transferoperator/common.jl b/src/binning_based/transferoperator/common.jl index 12824913c..f23fe0a96 100644 --- a/src/binning_based/transferoperator/common.jl +++ b/src/binning_based/transferoperator/common.jl @@ -17,15 +17,15 @@ A probability estimator based on binning data into rectangular boxes (when `ϵ` is a [`RectangularBinning`](@ref)) or simplices (when `ϵ` is a [`SimplexExact`](@ref) or `ϵ` is a [`SimplexPoint`](@ref)). -Input data are assumed to be sequential. +To use triangulation-based estimators, run `using Simplices` after `using Entropies`. The transfer (Perron-Frobenius) operator is approximated over the bins, using different methods depending on the type of partition chosen. The invariant measure associated with the approximated transfer operator is taken the bin probabilities. -This implementation follows the grid estimator approach in Diego et al. (2019)[^Diego2019]. +# Description (rectangular binning) -# Description +This implementation follows the grid estimator approach in Diego et al. (2019)[^Diego2019]. The transfer operator ``P^{N}``is computed as an `N`-by-`N` matrix of transition probabilities between the states defined by the partition elements, where `N` is the @@ -50,6 +50,9 @@ of jumping from the state defined by box ``C_i`` to any of the other ``N`` state follows that ``\\sum_{k=1}^{N} P_{ik} = 1`` for all ``i``. Thus, ``P^N`` is a row/right stochastic matrix. +For a detailed description of the triangulation estimators, see +Diego et al. (2019)[^Diego2019]. + ### Invariant measure estimation from transfer operator The left invariant distribution ``\\mathbf{\\rho}^N`` is a row vector, where @@ -81,11 +84,14 @@ See also: [`RectangularBinning`](@ref), [`SimplexPoint`](@ref), [`SimplexExact`] Here, we create three different transfer operator-based estimators. ```@example +using Entropies # A rectangular binning is suited for datasets with a large number of points est_rect = TransferOperator(RectangularBinning(5)) # A triangulated binning, using approximate simplex intersections, is also possible for -# datasets with not too many points (say, <1000 points). +# datasets with not too many points (say, <1000 points). If so, we must first import +# the Simplices.jl package. +using Simplices est_point = TransferOperator(SimplexPoint()) # For datasets with few points, say <100 points, exact simplex intersections may also From a11f58b32b84f1ac9b250f64fbfc2abd58bf9925 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 14:53:34 +0100 Subject: [PATCH 12/29] Add compat for Requires --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 7468e4fa3..e49553505 100644 --- a/Project.toml +++ b/Project.toml @@ -21,6 +21,7 @@ Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" DelayEmbeddings = "^1.20" Distances = "0.9, ^0.10" NearestNeighbors = "^0.4" +Requires = "^1.1" SpecialFunctions = "^0.10, ^1" StaticArrays = "^0.12, ^1" Wavelets = "^0.9" From f57fb4108563c6b49f5f458917ad4c403e487bc0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 14:59:01 +0100 Subject: [PATCH 13/29] Install python dependencies in tests --- test/install_dependencies.jl | 0 test/runtests.jl | 2 ++ 2 files changed, 2 insertions(+) create mode 100644 test/install_dependencies.jl diff --git a/test/install_dependencies.jl b/test/install_dependencies.jl new file mode 100644 index 000000000..e69de29bb diff --git a/test/runtests.jl b/test/runtests.jl index c8da185d4..c7fb5d51a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,5 @@ +include("install_dependencies.jl") + using Test using Entropies using DelayEmbeddings From 6e8966a6009909e1c8fbbb187339b28d668802ad Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 15:16:58 +0100 Subject: [PATCH 14/29] Temporarily reduce version to fix tests --- Project.toml | 10 +++++++--- test/Project.toml | 13 +++++++++++++ test/install_dependencies.jl | 0 test/runtests.jl | 12 +++++++++++- 4 files changed, 31 insertions(+), 4 deletions(-) create mode 100644 test/Project.toml delete mode 100644 test/install_dependencies.jl diff --git a/Project.toml b/Project.toml index e49553505..9c96cc252 100644 --- a/Project.toml +++ b/Project.toml @@ -2,15 +2,20 @@ name = "Entropies" uuid = "ed8fcbec-b94c-44b6-89df-898894ad9591" authors = "Kristian Agasøster Haaga , George Datseries " repo = "https://github.com/juliadynamics/Entropies.jl.git" -version = "1.0.0" +version = "0.11.1" [deps] DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" +DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" +PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" Requires = "ae029012-a4dd-5104-9daa-d747884805df" +Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" @@ -28,8 +33,7 @@ Wavelets = "^0.9" julia = "1.5" [extras] -Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Test", "Simplices"] +test = ["Test"] diff --git a/test/Project.toml b/test/Project.toml new file mode 100644 index 000000000..62e44b988 --- /dev/null +++ b/test/Project.toml @@ -0,0 +1,13 @@ +[deps] +Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" +DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f" +DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" +PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" +Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" + +[compat] +DynamicalSystems = "1.6, ^1" \ No newline at end of file diff --git a/test/install_dependencies.jl b/test/install_dependencies.jl deleted file mode 100644 index e69de29bb..000000000 diff --git a/test/runtests.jl b/test/runtests.jl index c7fb5d51a..df5573218 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,14 @@ -include("install_dependencies.jl") +if VERSION >= v"0.7.0-" + # Adding Pkg in test/REQUIRE would be an error in 0.6. Using + # Project.toml still has some gotchas. So: + Pkg = Base.require(Base.PkgId(Base.UUID(0x44cfe95a1eb252eab672e2afdf69b78f), "Pkg")) +end + +ENV["PYTHON"] = "" +Pkg.build("PyCall") + +using Conda +Conda.add("scipy") using Test using Entropies From 923cb124dd11474b5672020578d6edfee297f9da Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 15:37:42 +0100 Subject: [PATCH 15/29] Attempt to use test/Project.toml --- Project.toml | 3 --- test/Project.toml | 3 ++- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/Project.toml b/Project.toml index 9c96cc252..04af08cd5 100644 --- a/Project.toml +++ b/Project.toml @@ -34,6 +34,3 @@ julia = "1.5" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" - -[targets] -test = ["Test"] diff --git a/test/Project.toml b/test/Project.toml index 62e44b988..0e2a635be 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -10,4 +10,5 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" [compat] -DynamicalSystems = "1.6, ^1" \ No newline at end of file +DynamicalSystems = "1.6, ^1" +julia = "1.5" From 6274d7eee0cb0a2220809477172aaf4b6470fd5a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 15:49:41 +0100 Subject: [PATCH 16/29] Try old variant, but with [extras] --- Project.toml | 13 +++++++++---- test/Project.toml | 14 -------------- 2 files changed, 9 insertions(+), 18 deletions(-) delete mode 100644 test/Project.toml diff --git a/Project.toml b/Project.toml index 04af08cd5..3359a4606 100644 --- a/Project.toml +++ b/Project.toml @@ -8,14 +8,10 @@ version = "0.11.1" DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" -Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" -PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" Requires = "ae029012-a4dd-5104-9daa-d747884805df" -Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" @@ -25,6 +21,7 @@ Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" [compat] DelayEmbeddings = "^1.20" Distances = "0.9, ^0.10" +DynamicalSystems = "1.6, ^1" NearestNeighbors = "^0.4" Requires = "^1.1" SpecialFunctions = "^0.10, ^1" @@ -33,4 +30,12 @@ Wavelets = "^0.9" julia = "1.5" [extras] +Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" +DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" +PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" +Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test", "Conda", "DynamicalSystems", "PyCall", "Simplices", "StatsBase"] diff --git a/test/Project.toml b/test/Project.toml deleted file mode 100644 index 0e2a635be..000000000 --- a/test/Project.toml +++ /dev/null @@ -1,14 +0,0 @@ -[deps] -Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" -DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f" -DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" -PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" -Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" - -[compat] -DynamicalSystems = "1.6, ^1" -julia = "1.5" From 2f71fd80db706bf058528918812d435a5b3bef00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 15:49:41 +0100 Subject: [PATCH 17/29] Try old variant, but with [extras] DynamicalSystems not needed --- Project.toml | 10 +++++++--- test/Project.toml | 14 -------------- 2 files changed, 7 insertions(+), 17 deletions(-) delete mode 100644 test/Project.toml diff --git a/Project.toml b/Project.toml index 04af08cd5..071649c25 100644 --- a/Project.toml +++ b/Project.toml @@ -8,12 +8,9 @@ version = "0.11.1" DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" -Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" -PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" Requires = "ae029012-a4dd-5104-9daa-d747884805df" Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -33,4 +30,11 @@ Wavelets = "^0.9" julia = "1.5" [extras] +Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" +PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" +Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["Test", "Conda", "PyCall", "Simplices", "StatsBase"] diff --git a/test/Project.toml b/test/Project.toml deleted file mode 100644 index 0e2a635be..000000000 --- a/test/Project.toml +++ /dev/null @@ -1,14 +0,0 @@ -[deps] -Conda = "8f4d0f93-b110-5947-807f-2305c1781a2d" -DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f" -DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" -PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" -Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" -StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" -StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" -Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" - -[compat] -DynamicalSystems = "1.6, ^1" -julia = "1.5" From 17f3bc46eee115569cf470ddcad4c11138b542fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 16:13:55 +0100 Subject: [PATCH 18/29] Update Project.toml --- Project.toml | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 1be4e70d2..88ed35976 100644 --- a/Project.toml +++ b/Project.toml @@ -8,7 +8,6 @@ version = "0.11.1" DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f" Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" -DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" Requires = "ae029012-a4dd-5104-9daa-d747884805df" @@ -21,7 +20,7 @@ Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a" [compat] DelayEmbeddings = "^1.20" Distances = "0.9, ^0.10" -DynamicalSystems = "1.6, ^1" +Distributions = "^0.24" NearestNeighbors = "^0.4" Requires = "^1.1" SpecialFunctions = "^0.10, ^1" From 371100081f6bf91321d6a34bd132fc81f7cf373e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 16:14:03 +0100 Subject: [PATCH 19/29] Remove unneccessary test imports --- test/runtests.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index df5573218..a9e680941 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -11,10 +11,8 @@ using Conda Conda.add("scipy") using Test -using Entropies -using DelayEmbeddings using Wavelets -using StaticArrays +using Entropies using Simplices @testset "Histogram estimation" begin From 60a94fa08bd26afa116ee8f69b83bf6db1d04ee3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 16:24:27 +0100 Subject: [PATCH 20/29] See if adding Pkg to [extras] solves CI test error --- Project.toml | 4 +++- test/runtests.jl | 9 +++------ 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/Project.toml b/Project.toml index 88ed35976..f4d5fd815 100644 --- a/Project.toml +++ b/Project.toml @@ -11,6 +11,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" Requires = "ae029012-a4dd-5104-9daa-d747884805df" +Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" @@ -34,6 +35,7 @@ PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0" Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [targets] -test = ["Test", "Conda", "PyCall", "Simplices", "StatsBase"] +test = ["Test", "Conda", "PyCall", "Simplices", "StatsBase", "Pkg"] diff --git a/test/runtests.jl b/test/runtests.jl index a9e680941..01f9782ba 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,9 +1,4 @@ -if VERSION >= v"0.7.0-" - # Adding Pkg in test/REQUIRE would be an error in 0.6. Using - # Project.toml still has some gotchas. So: - Pkg = Base.require(Base.PkgId(Base.UUID(0x44cfe95a1eb252eab672e2afdf69b78f), "Pkg")) -end - +using Pkg ENV["PYTHON"] = "" Pkg.build("PyCall") @@ -14,6 +9,8 @@ using Test using Wavelets using Entropies using Simplices +using DelayEmbeddings +using StaticArrays @testset "Histogram estimation" begin x = rand(1:10, 100) From 535c5af2fdebe314bcc244fa501843ad3914fc7d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 17:21:44 +0100 Subject: [PATCH 21/29] Make sure to return named tuples for parameters --- .../rectangular/rectangular.jl | 19 +++++++++++++++++-- .../triangular/exact/SimplexExact.jl | 2 +- .../triangular/point/SimplexPoint.jl | 2 +- 3 files changed, 19 insertions(+), 4 deletions(-) diff --git a/src/binning_based/transferoperator/rectangular/rectangular.jl b/src/binning_based/transferoperator/rectangular/rectangular.jl index f7b086f9c..998df1ccb 100644 --- a/src/binning_based/transferoperator/rectangular/rectangular.jl +++ b/src/binning_based/transferoperator/rectangular/rectangular.jl @@ -212,7 +212,7 @@ function (tog::TransferOperatorGenerator{T})(;boundary_condition = :circular, # row/column of the transfer operator unique!(visited_bins) bins = [β .* edgelengths .+ mini for β in visited_bins] - params = (bins = bins) + params = (bins = bins, ) TransferOperatorApproximation(tog, M, params) end @@ -222,4 +222,19 @@ function transferoperator(pts, method::TransferOperator{<:RectangularBinning}; ) tog = transopergenerator(pts, method) tog(boundary_condition = boundary_condition) -end \ No newline at end of file +end + +""" + transfermatrix(iv::InvariantMeasure) → (M::AbstractArray{<:Real, 2}, bins::Vector{<:SVector}) + +Return the transfer matrix/operator and corresponding bins. Here, `bins[i]` corresponds +to the i-th row/column of the transfer matrix. Thus, the entry `M[i, j]` is the +probability of jumping from the state defined by `bins[i]` to the state defined by +`bins[j]`. + +See also: [`TransferOperator`](@ref). +""" +function transfermatrix(to::TransferOperatorApproximation{<:TransferOperator{RectangularBinning}}) + return to.transfermatrix, to.bins +end + diff --git a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl index 1956cbf49..d5fdf61f5 100644 --- a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl +++ b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl @@ -34,7 +34,7 @@ function transopergenerator(pts, method::TransferOperator{<:SimplexExact}) # the last vertex also can be mapped forward one step in time. triang = DelaunayTriangulation(invariant_pts[1:end-1]) - init = (invariant_pts = invariant_pts, triang = triang) + init = (invariant_pts = invariant_pts, triang = triang, ) TransferOperatorGenerator(method, pts, init) end diff --git a/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl b/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl index ec375af39..38f0e598b 100644 --- a/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl +++ b/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl @@ -118,7 +118,7 @@ function (tog::TransferOperatorGenerator{T})(; tol = 1e-8, randomsampling::Bool # Need to normalise, because all we have up until now is counts # of how many points inside the image simplex falls into # the simplices. - params = (tol = tol, randomsampling = randomsampling, n = n) + params = (tol = tol, randomsampling = randomsampling, n = n, ) M_normalized = transpose(M) ./ n_coeffs return TransferOperatorApproximation(tog, M_normalized, params) From 9be9310d39b96a25abbeb33f845e90015ac1a38a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 19:42:06 +0100 Subject: [PATCH 22/29] Rename transfermatrix to transitioninfo, and use return struct --- docs/src/TransferOperator.md | 8 +- src/binning_based/transferoperator/common.jl | 93 ++++++++++++------- .../rectangular/rectangular.jl | 19 ++-- .../triangular/exact/SimplexExact.jl | 6 ++ .../triangular/point/SimplexPoint.jl | 5 + test/runtests.jl | 16 +++- 6 files changed, 103 insertions(+), 44 deletions(-) diff --git a/docs/src/TransferOperator.md b/docs/src/TransferOperator.md index e7f7bc4d2..1781b2b3a 100644 --- a/docs/src/TransferOperator.md +++ b/docs/src/TransferOperator.md @@ -20,5 +20,11 @@ For rectangular binnings, use [`RectangularBinning`](@ref). ```@docs InvariantMeasure invariantmeasure -transfermatrix +``` + +## Transition information + +```@docs +transitioninfo +TransitionInfo ``` diff --git a/src/binning_based/transferoperator/common.jl b/src/binning_based/transferoperator/common.jl index f23fe0a96..6ddba83c3 100644 --- a/src/binning_based/transferoperator/common.jl +++ b/src/binning_based/transferoperator/common.jl @@ -6,9 +6,11 @@ export transferoperator, TransferOperator, TransferOperatorApproximation, + TransferMatrix, invariantmeasure, InvariantMeasure, - transfermatrix + transitioninfo, + TransitionInfo """ TransferOperator(ϵ::Union{RectangularBinning, SimplexPoint, SimplexExact}) <: BinningProbabilitiesEstimator @@ -148,30 +150,6 @@ function transferoperator(pts, method::TransferOperator) to() end -""" - InvariantMeasure(to, ρ) - -Minimal return struct for [`invariantmeasure`](@ref) that contains the estimated invariant -measure `ρ`, as well as the transfer operator `to` from which it is computed (including -bin information). - -See also: [`invariantmeasure`](@ref). -""" -struct InvariantMeasure{T, P<:Probabilities} - to::T - ρ::P - - function InvariantMeasure(to::T, ρ::P) where {T, P} - new{T, P}(to, ρ) - end -end - - -function Base.show(io::IO, DT::InvariantMeasure{T, P}) where {T, P} - summary = "InvariantMeasure{transfer operator approximation: $T, probabilities: $P}" - println(io, summary) -end - """ TransferOperatorApproximation(generator, transfermatrix, params) @@ -218,20 +196,69 @@ end """ - transfermatrix(iv::InvariantMeasure) → (M::AbstractArray{<:Real, 2}, bins::Vector{<:SVector}) + TransitionInfo(transfermatrix, bins) + +Structure that holds the `transfermatrix` of transition probabilities obtained from a +transfer operator approximation, as well as the `bins` of the partition. + +- If `transfermatrix` was obtained using a rectangular binning, then `bins` is a + two-element named tuple with fields `bins` (the left-most coordinates of the bins) and + `edgelengths` (the edge lengths along each coordinate axis). +- If `transfermatrix` was obtained using a triangulated binning, then `bins` is a + two-element named tuple with fields `pts` (the points of the triangulation), and + `triang`, the indices of the vertices of all simplices. Here, `pts[triang[i]]` is a + vector of vertices for the `i`-th simplex of the triangulation. -Return the transfer matrix/operator and corresponding bins. Here, `bins[i]` corresponds -to the i-th row/column of the transfer matrix. Thus, the entry `M[i, j]` is the -probability of jumping from the state defined by `bins[i]` to the state defined by -`bins[j]`. +`bins[i]` corresponds to the i-th row/column of `transfermatrix.` Thus, the entry +`transfermatrix[i, j]` is the probability of jumping from the state defined by `bins[i]` +to the state defined by `bins[j]`. -See also: [`TransferOperator`](@ref). +See also [`transitioninfo`](@ref). """ -function transfermatrix(iv::InvariantMeasure) - return iv.to.transfermatrix, iv.to.bins +struct TransitionInfo{T, B} + transfermatrix::T + bins::B end +function Base.show(io::IO, info::TransitionInfo{T, B}) where {T, B} + summary = "TransitionInfo{transfermatrix: $T, bins: $B}" + println(io, summary) +end + +""" + transitioninfo(to::TransferOperatorApproximation) → TransitionInfo + transitioninfo(to::InvariantMeasure) → TransitionInfo + +Convenience method to get the transfer matrix/operator and the corresponding bins from +pre-computed quantities. + +See also: [`TransitionInfo`](@ref), [`TransferOperatorApproximation`](@ref). +""" +function transitioninfo end + +""" + InvariantMeasure(to, ρ) + +Minimal return struct for [`invariantmeasure`](@ref) that contains the estimated invariant +measure `ρ`, as well as the transfer operator `to` from which it is computed (including +bin information). + +See also: [`invariantmeasure`](@ref). +""" +struct InvariantMeasure{T, P<:Probabilities} + to::T + ρ::P + + function InvariantMeasure(to::T, ρ::P) where {T, P} + new{T, P}(to, ρ) + end +end + +function Base.show(io::IO, DT::InvariantMeasure{T, P}) where {T, P} + summary = "InvariantMeasure{transfer operator approximation: $T, probabilities: $P}" + println(io, summary) +end import LinearAlgebra: norm """ diff --git a/src/binning_based/transferoperator/rectangular/rectangular.jl b/src/binning_based/transferoperator/rectangular/rectangular.jl index 998df1ccb..dc1d796a7 100644 --- a/src/binning_based/transferoperator/rectangular/rectangular.jl +++ b/src/binning_based/transferoperator/rectangular/rectangular.jl @@ -212,7 +212,7 @@ function (tog::TransferOperatorGenerator{T})(;boundary_condition = :circular, # row/column of the transfer operator unique!(visited_bins) bins = [β .* edgelengths .+ mini for β in visited_bins] - params = (bins = bins, ) + params = (bins = bins, edgelengths = edgelengths) TransferOperatorApproximation(tog, M, params) end @@ -225,16 +225,19 @@ function transferoperator(pts, method::TransferOperator{<:RectangularBinning}; end """ - transfermatrix(iv::InvariantMeasure) → (M::AbstractArray{<:Real, 2}, bins::Vector{<:SVector}) + transitioninfo(to::TransferOperatorApproximation) → TransitionInfo + transitioninfo(to::InvariantMeasure) → TransitionInfo -Return the transfer matrix/operator and corresponding bins. Here, `bins[i]` corresponds -to the i-th row/column of the transfer matrix. Thus, the entry `M[i, j]` is the -probability of jumping from the state defined by `bins[i]` to the state defined by +Convenience method to get the transfer matrix/operator and corresponding bins. Here, +`bins[i]` corresponds to the i-th row/column of the transfer matrix. Thus, the entry +`M[i, j]` is the probability of jumping from the state defined by `bins[i]` to the state defined by `bins[j]`. -See also: [`TransferOperator`](@ref). +See also: [`TransitionInfo`](@ref), [`TransferOperatorApproximation`](@ref). """ -function transfermatrix(to::TransferOperatorApproximation{<:TransferOperator{RectangularBinning}}) - return to.transfermatrix, to.bins +function transitioninfo(to::TransferOperatorApproximation{<:TransferOperator{<:RectangularBinning}}) + return TransitionInfo(to.transfermatrix, + (bins = to.params.bins, edgelengths = to.params.edgelengths)) end +transitioninfo(iv::InvariantMeasure) = transitioninfo(iv.to) diff --git a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl index d5fdf61f5..6319efcec 100644 --- a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl +++ b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl @@ -82,3 +82,9 @@ function transferoperator(pts, method::TransferOperator{<:SimplexExact}; tol::Re tog = transopergenerator(pts, method) tog(tol = tol) end + + +function transitioninfo(to::TransferOperatorApproximation{<:TransferOperator{R}}) where R<:SimplexExact + return TransitionInfo(to.transfermatrix, + (pts = to.generator.init.invariant_pts, triang = to.generator.init.triang,)) +end \ No newline at end of file diff --git a/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl b/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl index 38f0e598b..f3f3774ba 100644 --- a/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl +++ b/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl @@ -129,4 +129,9 @@ function transferoperator(pts, method::TransferOperator{<:SimplexPoint}; tog = transopergenerator(pts, method) tog(tol = tol, randomsampling = randomsampling, n = n) +end + +function transitioninfo(to::TransferOperatorApproximation{<:TransferOperator{R}}) where R<:SimplexPoint + return TransitionInfo(to.transfermatrix, + (pts = to.generator.init.invariant_pts, triang = to.generator.init.triang,)) end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 01f9782ba..c01e04170 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -354,9 +354,13 @@ end @testset "Rectagular binning test $i" for i in 1:length(binnings) est = TransferOperator(binnings[i]) to = Entropies.transferoperator(D, est) + + @test transitioninfo(to) isa Entropies.TransitionInfo @test to isa Entropies.TransferOperatorApproximation iv = invariantmeasure(to) + @test transitioninfo(iv) isa Entropies.TransitionInfo + @test iv isa InvariantMeasure @test probabilities(D, TransferOperator(binnings[i])) isa Probabilities end @@ -369,11 +373,19 @@ end to_point = transferoperator(D, est_point) to_exact = transferoperator(D, est_exact) + + @test transitioninfo(to_exact) isa Entropies.TransitionInfo + @test transitioninfo(to_point) isa Entropies.TransitionInfo + @test to_point isa Entropies.TransferOperatorApproximation @test to_exact isa Entropies.TransferOperatorApproximation - @test invariantmeasure(to_point) isa Entropies.InvariantMeasure - @test invariantmeasure(to_exact) isa Entropies.InvariantMeasure + iv_point = invariantmeasure(to_point) + iv_exact = invariantmeasure(to_exact) + @test iv_point isa Entropies.InvariantMeasure + @test iv_exact isa Entropies.InvariantMeasure + @test transitioninfo(iv_point) isa Entropies.TransitionInfo + @test transitioninfo(iv_exact) isa Entropies.TransitionInfo @test probabilities(to_point) isa Probabilities @test probabilities(to_exact) isa Probabilities From ecd1ff2787d7a0017436f31179f3e446c1e7f650 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 20:14:20 +0100 Subject: [PATCH 23/29] Update docs. --- src/binning_based/transferoperator/common.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/binning_based/transferoperator/common.jl b/src/binning_based/transferoperator/common.jl index 6ddba83c3..4c7f143e6 100644 --- a/src/binning_based/transferoperator/common.jl +++ b/src/binning_based/transferoperator/common.jl @@ -72,9 +72,9 @@ below some threshold. ## Probability and entropy estimation -- `probabilities(x::AbstractDataset, est::TransferOperator{RectangularBinning})` estimates +- `probabilities(x::AbstractDataset, est::TransferOperator)` estimates probabilities for the bins defined by the provided binning (`est.ϵ`) -- `genentropy(x::AbstractDataset, est::TransferOperator{RectangularBinning})` does the same, +- `genentropy(x::AbstractDataset, est::TransferOperator)` does the same, but computes generalized entropy using the probabilities. @@ -238,7 +238,7 @@ See also: [`TransitionInfo`](@ref), [`TransferOperatorApproximation`](@ref). function transitioninfo end """ - InvariantMeasure(to, ρ) + InvariantMeasure(to::TransferOperatorApproximation, ρ::Probabilities) Minimal return struct for [`invariantmeasure`](@ref) that contains the estimated invariant measure `ρ`, as well as the transfer operator `to` from which it is computed (including @@ -246,7 +246,7 @@ bin information). See also: [`invariantmeasure`](@ref). """ -struct InvariantMeasure{T, P<:Probabilities} +struct InvariantMeasure{T<:TransferOperatorApproximation, P<:Probabilities} to::T ρ::P From 829f14bca3dff9cc8160d8faad551302d30052b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 20:14:42 +0100 Subject: [PATCH 24/29] Add examples to estimator docstring --- src/binning_based/transferoperator/common.jl | 28 +++++++++++++++----- 1 file changed, 21 insertions(+), 7 deletions(-) diff --git a/src/binning_based/transferoperator/common.jl b/src/binning_based/transferoperator/common.jl index 4c7f143e6..c4a472ccf 100644 --- a/src/binning_based/transferoperator/common.jl +++ b/src/binning_based/transferoperator/common.jl @@ -85,20 +85,34 @@ See also: [`RectangularBinning`](@ref), [`SimplexPoint`](@ref), [`SimplexExact`] Here, we create three different transfer operator-based estimators. +A rectangular binning is suited for datasets with a large number of points. + ```@example using Entropies -# A rectangular binning is suited for datasets with a large number of points +pts = Dataset(rand(2000, 3)) est_rect = TransferOperator(RectangularBinning(5)) +genentropy(pts, est_rect) +``` -# A triangulated binning, using approximate simplex intersections, is also possible for -# datasets with not too many points (say, <1000 points). If so, we must first import -# the Simplices.jl package. -using Simplices +A triangulated binning, using approximate simplex intersections, is also possible for +datasets with not too many points (say, <1000 points). If so, we must first import +the Simplices.jl package. + +```@example +using Entropies, Simplices +pts = Dataset(rand(200, 3)) est_point = TransferOperator(SimplexPoint()) +genentropy(pts, est_point) +``` + +For datasets with few points, say <100 points, exact simplex intersections may also +be computationally feasible. -# For datasets with few points, say <100 points, exact simplex intersections may also -# be computationally feasible. +```@example +using Entropies, Simplices est_exact = TransferOperator(SimplexExact()) +pts = Dataset(rand(40, 3)) +genentropy(pts, est_exact) ``` [^Diego2019]: Diego, D., Haaga, K. A., & Hannisdal, B. (2019). Transfer entropy computation using the Perron-Frobenius operator. Physical Review E, 99(4), 042212. From 3cdc3a96f582974ccffaf8b35c46412c2689e220 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 20:28:27 +0100 Subject: [PATCH 25/29] Abstract type for triangulation estimators --- .../transferoperator/triangular/exact/SimplexExact.jl | 11 ++++++----- .../transferoperator/triangular/point/SimplexPoint.jl | 9 +++++---- .../transferoperator/triangular/triangular.jl | 2 ++ test/runtests.jl | 8 ++++++++ 4 files changed, 21 insertions(+), 9 deletions(-) diff --git a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl index 6319efcec..ba95ad7fd 100644 --- a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl +++ b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl @@ -1,20 +1,21 @@ export SimplexExact, invariantmeasure """ - SimplexExact + SimplexExact <: TriangularBinning A transfer operator estimator using a triangulation partition and exact simplex intersections[^Diego2019]. -To use this estimator, the Simplices.jl package must be brought into scope by doing -`using Simplices` after running `using Entropies`. - *Note: due to computing exact simplex intersections, this estimator is slow compared to [`SimplexPoint`](@ref).* +!!! hint + To use this estimator, the Simplices.jl package must be brought into scope by doing + `using Simplices` after running `using Entropies`. + [^Diego2019]: Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212. """ -struct SimplexExact +struct SimplexExact <: TriangularBinning bc::String function SimplexExact(bc::String = "circular") diff --git a/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl b/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl index f3f3774ba..ab8c463ed 100644 --- a/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl +++ b/src/binning_based/transferoperator/triangular/point/SimplexPoint.jl @@ -4,17 +4,18 @@ import StaticArrays: SizedVector, SizedMatrix, MMatrix include("helper_functions.jl") """ - SimplexPoint + SimplexPoint <: TriangularBinning A transfer operator estimator using a triangulated partition, and using approximate simplex intersections to compute transition probabilities [^Diego2019]. -To use this estimator, the Simplices.jl package must be brought into scope by doing -`using Simplices` after running `using Entropies`. +!!! hint + To use this estimator, the Simplices.jl package must be brought into scope by doing + `using Simplices` after running `using Entropies`. [^Diego2019]: Diego, David, Kristian Agasøster Haaga, and Bjarte Hannisdal. "Transfer entropy computation using the Perron-Frobenius operator." Physical Review E 99.4 (2019): 042212. """ -struct SimplexPoint +struct SimplexPoint <: TriangularBinning bc::String function SimplexPoint(bc::String = "circular") diff --git a/src/binning_based/transferoperator/triangular/triangular.jl b/src/binning_based/transferoperator/triangular/triangular.jl index 2e7f10531..2e2af77b4 100644 --- a/src/binning_based/transferoperator/triangular/triangular.jl +++ b/src/binning_based/transferoperator/triangular/triangular.jl @@ -2,5 +2,7 @@ include("simplex_types/simplex_types.jl") include("./delaunay_triangulations/DelaunayTriangulations.jl") import .DelaunayTriangulations: DelaunayTriangulation +abstract type TriangularBinning end + include("exact/SimplexExact.jl") include("point/SimplexPoint.jl") \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index c01e04170..d35530739 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -351,6 +351,14 @@ end RectangularBinning([0.2, 0.3, 0.3]), ] + @testset "Methods" begin + @testset "Rectagular binning test $i" for i in 1:length(binnings) + @test binnings[i] isa RectangularBinning + end + @test SimplexPoint() isa Entropies.TriangularBinning + @test SimplexExact() isa Entropies.TriangularBinning + end + @testset "Rectagular binning test $i" for i in 1:length(binnings) est = TransferOperator(binnings[i]) to = Entropies.transferoperator(D, est) From b33b1903aef1119c416cb4a538001078bdb282c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 20:28:34 +0100 Subject: [PATCH 26/29] Restructure docs --- docs/src/TransferOperator.md | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/src/TransferOperator.md b/docs/src/TransferOperator.md index 1781b2b3a..615e7f882 100644 --- a/docs/src/TransferOperator.md +++ b/docs/src/TransferOperator.md @@ -22,9 +22,11 @@ InvariantMeasure invariantmeasure ``` -## Transition information +### Transition information ```@docs +transferoperator +TransferOperatorApproximation transitioninfo TransitionInfo ``` From 24d4106e2559543f829c2819f4793b9c1244dfd1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sat, 20 Feb 2021 20:45:37 +0100 Subject: [PATCH 27/29] Add missing docs --- src/binning_based/transferoperator/common.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/binning_based/transferoperator/common.jl b/src/binning_based/transferoperator/common.jl index c4a472ccf..cd3e47385 100644 --- a/src/binning_based/transferoperator/common.jl +++ b/src/binning_based/transferoperator/common.jl @@ -137,7 +137,11 @@ function Base.show(io::IO, DT::TransferOperatorGenerator{E, X, A}) where {E, X, println(io, summary) end +""" + transferoperator(points, estimator::TransferOperator) → TransferOperatorApproximation +Approximate the transfer operator over `points` using the provided estimator. +""" function transferoperator end """ @@ -325,7 +329,7 @@ The element `ρ[i]` is the probability of visitation to the box `bins[i]`. Analo The naive histogram approach only gives the long-term probabilities that orbits visit a certain region of the state space. The transfer operator encodes that information too, but comes with the added benefit of knowing the *transition - probabilities* between states (see [`transfermatrix`](@ref)). + probabilities* between states (see [`transitioninfo`](@ref)). See also: [`InvariantMeasure`](@ref). """ From 5559ccbcc3c1f4c6461c18de0289135d70d1b614 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Sun, 21 Feb 2021 00:56:14 +0100 Subject: [PATCH 28/29] Make sure to return named tuple --- .../transferoperator/rectangular/rectangular.jl | 12 +----------- .../triangular/exact/SimplexExact.jl | 2 +- 2 files changed, 2 insertions(+), 12 deletions(-) diff --git a/src/binning_based/transferoperator/rectangular/rectangular.jl b/src/binning_based/transferoperator/rectangular/rectangular.jl index dc1d796a7..fb9be6c5d 100644 --- a/src/binning_based/transferoperator/rectangular/rectangular.jl +++ b/src/binning_based/transferoperator/rectangular/rectangular.jl @@ -224,20 +224,10 @@ function transferoperator(pts, method::TransferOperator{<:RectangularBinning}; tog(boundary_condition = boundary_condition) end -""" - transitioninfo(to::TransferOperatorApproximation) → TransitionInfo - transitioninfo(to::InvariantMeasure) → TransitionInfo -Convenience method to get the transfer matrix/operator and corresponding bins. Here, -`bins[i]` corresponds to the i-th row/column of the transfer matrix. Thus, the entry -`M[i, j]` is the probability of jumping from the state defined by `bins[i]` to the state defined by -`bins[j]`. - -See also: [`TransitionInfo`](@ref), [`TransferOperatorApproximation`](@ref). -""" function transitioninfo(to::TransferOperatorApproximation{<:TransferOperator{<:RectangularBinning}}) return TransitionInfo(to.transfermatrix, - (bins = to.params.bins, edgelengths = to.params.edgelengths)) + (bins = to.params.bins, edgelengths = to.params.edgelengths, )) end transitioninfo(iv::InvariantMeasure) = transitioninfo(iv.to) diff --git a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl index ba95ad7fd..c83400012 100644 --- a/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl +++ b/src/binning_based/transferoperator/triangular/exact/SimplexExact.jl @@ -87,5 +87,5 @@ end function transitioninfo(to::TransferOperatorApproximation{<:TransferOperator{R}}) where R<:SimplexExact return TransitionInfo(to.transfermatrix, - (pts = to.generator.init.invariant_pts, triang = to.generator.init.triang,)) + (pts = to.generator.init.invariant_pts, triang = to.generator.init.triang, )) end \ No newline at end of file From 89a19f2a855f9b6b6d6d9a82095320d5eb38d9c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Kristian=20Agas=C3=B8ster=20Haaga?= Date: Thu, 8 Apr 2021 12:09:06 +0200 Subject: [PATCH 29/29] Update tagbot --- .github/workflows/TagBot.yml | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml index dc5991aba..b798528c0 100644 --- a/.github/workflows/TagBot.yml +++ b/.github/workflows/TagBot.yml @@ -1,13 +1,15 @@ -name: TagBot +name: Julia TagBot on: - schedule: - - cron: 0 * * * * + issue_comment: + types: + - created + workflow_dispatch: jobs: TagBot: + if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' runs-on: ubuntu-latest steps: - uses: JuliaRegistries/TagBot@v1 with: token: ${{ secrets.GITHUB_TOKEN }} - ssh: ${{ secrets.DOCUMENTER_KEY }} - + ssh: ${{ secrets.DOCUMENTER_KEY }} \ No newline at end of file