diff --git a/Manifest.toml b/Manifest.toml index 11eb08697..f619a50a6 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -155,9 +155,9 @@ version = "0.3.2" [[DiffEqBase]] deps = ["Compat", "IteratorInterfaceExtensions", "LinearAlgebra", "RecipesBase", "RecursiveArrayTools", "Requires", "Roots", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "TableTraits", "Test", "TreeViews"] -git-tree-sha1 = "a79e0d74bdc91d4dfc501fe03af762b48cbc3764" +git-tree-sha1 = "b8c9c5a296e426d0ad978e6e7ad699dde444d3a1" uuid = "2b5f629d-d688-5b77-993f-72d75c75574e" -version = "4.31.0" +version = "4.31.1" [[DiffEqBiological]] deps = ["Compat", "DataStructures", "DiffEqBase", "DiffEqJump", "MacroTools", "Random", "Reexport", "Statistics", "SymEngine", "Test"] @@ -232,10 +232,10 @@ uuid = "055956cb-9e8b-5191-98cc-73ae4a59e68a" version = "3.1.0" [[DiffEqSensitivity]] -deps = ["Compat", "DiffEqBase", "DiffEqCallbacks", "DiffEqDiffTools", "ForwardDiff", "LinearAlgebra", "QuadGK", "RecursiveArrayTools", "Statistics", "Test"] -git-tree-sha1 = "cc4716516bd4a60e092e8a909efcbd26170cd358" +deps = ["Compat", "DiffEqBase", "DiffEqCallbacks", "DiffEqDiffTools", "ForwardDiff", "LinearAlgebra", "QuadGK", "RecursiveArrayTools", "ReverseDiff", "Statistics", "Test"] +git-tree-sha1 = "76f33416ee523314b7fcdc08696408861489c427" uuid = "41bf760c-e81c-5289-8e54-58b1f1f8abe2" -version = "2.2.0" +version = "2.3.0" [[DiffEqUncertainty]] deps = ["DiffEqBase", "Test"] @@ -378,9 +378,9 @@ uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240" [[Interpolations]] deps = ["AxisAlgorithms", "LinearAlgebra", "OffsetArrays", "Random", "Ratios", "SharedArrays", "SparseArrays", "StaticArrays", "Test", "WoodburyMatrices"] -git-tree-sha1 = "2a3b97b4f4cd6f6be5c73f07b15a30e9eaba465d" +git-tree-sha1 = "3493536a64dae5a21c0cc8aecf680647f3e12313" uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" -version = "0.10.6" +version = "0.11.0" [[IteratorInterfaceExtensions]] deps = ["Test"] @@ -518,9 +518,9 @@ version = "1.0.2" [[OrdinaryDiffEq]] deps = ["DataStructures", "DiffEqBase", "DiffEqDiffTools", "DiffEqOperators", "ExponentialUtilities", "ForwardDiff", "GenericSVD", "InteractiveUtils", "LinearAlgebra", "Logging", "MuladdMacro", "NLsolve", "Parameters", "Printf", "Random", "RecursiveArrayTools", "Reexport", "SparseArrays", "StaticArrays", "Statistics", "Test"] -git-tree-sha1 = "44b786ca14842e30eba2656d83c7055f8a20b74d" +git-tree-sha1 = "06b6488f3f17f816de2c96fc071a8115b84ca736" uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" -version = "4.17.2" +version = "4.18.0" [[PDMats]] deps = ["Arpack", "LinearAlgebra", "SparseArrays", "SuiteSparse", "Test"] @@ -548,9 +548,9 @@ version = "0.1.2" [[PerronFrobenius]] deps = ["Conda", "Distributed", "Documenter", "InplaceOps", "LinearAlgebra", "Parameters", "Printf", "PyCall", "RecipesBase", "Reexport", "SharedArrays", "Simplices", "SparseArrays", "StateSpaceReconstruction", "StaticArrays", "Test"] -git-tree-sha1 = "6c1471a5c3c36859c62e48d8b4c793ac16fae4d8" +git-tree-sha1 = "7695c740e987c4d8902dfc6375810a128654393a" uuid = "260eed61-d0e8-5f1e-b040-a9756a401c97" -version = "0.2.0" +version = "0.2.2" [[Pkg]] deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"] @@ -658,6 +658,12 @@ git-tree-sha1 = "263fe0a06007bf5ba8f72c55475f58fc9bf81834" uuid = "ae5879a3-cd67-5da8-be7f-38c6eb64a37b" version = "0.5.0" +[[ReverseDiff]] +deps = ["DiffResults", "DiffRules", "ForwardDiff", "FunctionWrappers", "LinearAlgebra", "NaNMath", "Random", "SpecialFunctions", "StaticArrays", "Statistics", "Test"] +git-tree-sha1 = "d44e80d815a1a5abaf1ff1def21bb629cabaffec" +uuid = "37e2e3b7-166d-5795-8a7a-e32c996b4267" +version = "0.3.1" + [[Rmath]] deps = ["BinaryProvider", "Libdl", "Random", "Statistics", "Test"] git-tree-sha1 = "9a6c758cdf73036c3239b0afbea790def1dabff9" @@ -694,9 +700,9 @@ version = "0.8.0" [[Simplices]] deps = ["Conda", "Distributions", "Documenter", "LinearAlgebra", "Parameters", "Printf", "PyCall", "StaticArrays", "Test"] -git-tree-sha1 = "4e96805d9ab3413ab3b82883b61fdea14dfcc33e" +git-tree-sha1 = "36b8816dc8297a39d3ad56b32e0cf61e333121fc" uuid = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" -version = "0.2.1" +version = "0.2.2" [[Sockets]] uuid = "6462fe0b-24de-5631-8697-dd941f90decc" @@ -718,10 +724,10 @@ uuid = "276daf66-3868-5448-9aa4-cd146d93841b" version = "0.7.2" [[StateSpaceReconstruction]] -deps = ["Conda", "Distributions", "LinearAlgebra", "PyCall", "Reexport", "Simplices", "StaticArrays", "Statistics", "Test"] -git-tree-sha1 = "1137310fbd695e00b5bcdaf68677b6383e6e1151" +deps = ["Conda", "Distances", "Distributions", "LinearAlgebra", "NearestNeighbors", "PyCall", "Reexport", "Simplices", "StaticArrays", "Statistics", "Test"] +git-tree-sha1 = "211298d65326fd109bb00c92ce0adca0f4c09435" uuid = "1441a9f6-6a74-5418-a591-cdf1d78a07f0" -version = "0.3.0" +version = "0.3.1" [[StaticArrays]] deps = ["InteractiveUtils", "LinearAlgebra", "Random", "Statistics", "Test"] @@ -735,9 +741,9 @@ uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" [[StatsBase]] deps = ["DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "Test"] -git-tree-sha1 = "723193a13e8078cec6dcd0b8fe245c8bfd81690e" +git-tree-sha1 = "2722397d88f8ffef551948f6c20e1d74a743298c" uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" -version = "0.25.0" +version = "0.26.0" [[StatsFuns]] deps = ["Rmath", "SpecialFunctions", "Test"] diff --git a/Project.toml b/Project.toml index c4d433a64..445e6c907 100644 --- a/Project.toml +++ b/Project.toml @@ -15,6 +15,8 @@ DocumenterMarkdown = "997ab1e6-3595-5248-9280-8efb232c3433" DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" PerronFrobenius = "260eed61-d0e8-5f1e-b040-a9756a401c97" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" @@ -24,6 +26,7 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" Simplices = "d5428e67-3037-59ba-9ab1-57a04f0a3b6a" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StateSpaceReconstruction = "1441a9f6-6a74-5418-a591-cdf1d78a07f0" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TimeseriesSurrogates = "c804724b-8c18-5caa-8579-6025a0767c70" diff --git a/README.md b/README.md index 22f3b73f1..3be42648b 100644 --- a/README.md +++ b/README.md @@ -18,11 +18,12 @@ Check out the documentation (coming) for more information! | package | functionality | version | build | | :--- | :--- | :---: | ---: | -| [`StateSpaceReconstruction.jl`](https://github.com/kahaaga/StateSpaceReconstruction.jl/) | Fully flexible state space reconstructions (embeddings), partitioning routines (variable-width rectangular, and triangulations), and partition refinement (equal-volume splitting of simplices). | `0.2.3` | [![Build Status](https://travis-ci.org/kahaaga/StateSpaceReconstruction.jl.svg?branch=master)](https://travis-ci.org/kahaaga/StateSpaceReconstruction.jl) | +| [`StateSpaceReconstruction.jl`](https://github.com/kahaaga/StateSpaceReconstruction.jl/) | Fully flexible state space reconstructions (embeddings), partitioning routines (variable-width rectangular, and triangulations), and partition refinement (equal-volume splitting of simplices). | `0.3.1` | [![Build Status](https://travis-ci.org/kahaaga/StateSpaceReconstruction.jl.svg?branch=master)](https://travis-ci.org/kahaaga/StateSpaceReconstruction.jl) | | [`TimeseriesSurrogates.jl`](https://github.com/kahaaga/TimeseriesSurrogates.jl/) | Generate surrogate data from time series. | `0.2.1` | [![Build Status](https://travis-ci.org/kahaaga/TimeseriesSurrogates.jl.svg?branch=master)](https://travis-ci.org/kahaaga/TimeseriesSurrogates.jl) | -| [`TransferEntropy.jl`](https://github.com/kahaaga/TransferEntropy.jl/) | Transfer entropy estimators. | `0.2.3` | [![Build Status](https://travis-ci.org/kahaaga/TransferEntropy.jl.svg?branch=master)](https://travis-ci.org/kahaaga/TransferEntropy.jl) | | -| [`PerronFrobenius.jl`](https://github.com/kahaaga/PerronFrobenius.jl/) | Transfer (Perron-Frobenius) operator estimators. | `0.1.3`  | [![Build Status](https://travis-ci.org/kahaaga/PerronFrobenius.jl.svg?branch=master)](https://travis-ci.org/kahaaga/PerronFrobenius.jl) | -| [`Simplices.jl`](https://github.com/kahaaga/Simplices.jl/) | Exact simplex intersections in N dimensions. | `0.2.1` | [![Build Status](https://travis-ci.org/kahaaga/Simplices.jl.svg?branch=master)](https://travis-ci.org/kahaaga/Simplices.jl) | +| [`TransferEntropy.jl`](https://github.com/kahaaga/TransferEntropy.jl/) | Transfer entropy estimators. | `0.3.1` | [![Build Status](https://travis-ci.org/kahaaga/TransferEntropy.jl.svg?branch=master)](https://travis-ci.org/kahaaga/TransferEntropy.jl) | | +| [`PerronFrobenius.jl`](https://github.com/kahaaga/PerronFrobenius.jl/) | Transfer (Perron-Frobenius) operator estimators. | `0.2.2`  | [![Build Status](https://travis-ci.org/kahaaga/PerronFrobenius.jl.svg?branch=master)](https://travis-ci.org/kahaaga/PerronFrobenius.jl) | +| [`Simplices.jl`](https://github.com/kahaaga/Simplices.jl/) | Exact simplex intersections in N dimensions. | `0.2.2` | [![Build Status](https://travis-ci.org/kahaaga/Simplices.jl.svg?branch=master)](https://travis-ci.org/kahaaga/Simplices.jl) | +| [`CrossMappings.jl`](https://github.com/kahaaga/CrossMappings.jl/) | Exact simplex intersections in N dimensions. | `0.2.0` | [![Build Status](https://travis-ci.org/kahaaga/CrossMappings.jl.svg?branch=master)](https://travis-ci.org/kahaaga/CrossMappings.jl) | ## Wrappers for common use cases diff --git a/REQUIRE b/REQUIRE index 9c19ed35c..dca833148 100644 --- a/REQUIRE +++ b/REQUIRE @@ -12,11 +12,12 @@ DynamicalSystems NearestNeighbors Interpolations PerronFrobenius 0.2.0 -Plots PyCall RecipesBase Reexport Simplices 0.2.1 +StaticArrays StateSpaceReconstruction 0.3.0 TimeseriesSurrogates 0.2.1 TransferEntropy 0.3.1 +Measures diff --git a/install_dependencies.jl b/install_dependencies.jl new file mode 100644 index 000000000..63f4f328e --- /dev/null +++ b/install_dependencies.jl @@ -0,0 +1,13 @@ +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 + +# Let PyCall.jl use Python interpreter from Conda.jl +# See: https://github.com/JuliaPy/PyCall.jl +ENV["PYTHON"] = "" +Pkg.build("PyCall") + +using Conda +Conda.add("scipy") diff --git a/runtests.jl b/runtests.jl new file mode 100644 index 000000000..b98cdc7b3 --- /dev/null +++ b/runtests.jl @@ -0,0 +1,22 @@ +if lowercase(get(ENV, "CI", "false")) == "true" + include("install_dependencies.jl") +end + +using TransferEntropy +using Test +using Distances + + +n_realizations = 50 + +@testset "Visitation frequency estimator" begin + include("transferentropy/test_transferentropy_visitfreq.jl") +end + +@testset "kNN (Kraskov) estimator" begin + include("transferentropy/test_transferentropy_kraskov.jl") +end + +@testset "Transfer operator grid estimator" begin + include("transferentropy/test_transferentropy_transferoperator_grid.jl") +end diff --git a/src/CausalityTools.jl b/src/CausalityTools.jl index 64b0dadda..59506e29e 100644 --- a/src/CausalityTools.jl +++ b/src/CausalityTools.jl @@ -6,7 +6,9 @@ using Reexport using AbstractFFTs using Statistics using RecipesBase -using Plots +using StaticArrays +using LinearAlgebra +using Measures @reexport using TimeseriesSurrogates @@ -31,6 +33,9 @@ include("surrogates/surrogates.jl") # Wrappers of the different methods. include("method_wrappers/highlevel_methods.jl") +# Plot recipes, also for all sub-packages +include("plot_recipes/recipes.jl") + export Dataset, DynamicalSystem, diff --git a/src/plot_recipes/helperfunctions_rectangulargrid.jl b/src/plot_recipes/helperfunctions_rectangulargrid.jl new file mode 100644 index 000000000..28616707c --- /dev/null +++ b/src/plot_recipes/helperfunctions_rectangulargrid.jl @@ -0,0 +1,112 @@ +######################################################### +# Helper functions to plot 3D rectangles. +######################################################### +""" + rectangle3dpts(x, y, z, ϵx, ϵy, ϵz) + +Compute the coordinates of the vertices of a 3D rectangle, given the +coordinates of the origin (`x`, `y`, and `z`) and the corresponding edge lengths +(`ϵx`, `ϵy`, `ϵz`). +""" +function rectangle3dpts(x, y, z, ϵx, ϵy, ϵz) + v1b = [x, y, z] + v2b = [x+ϵx, y, z] + v3b = [x+ϵx, y+ϵy, z] + v4b = [x, y+ϵy, z] + v1t = [x, y, z+ϵz] + v2t = [x+ϵx, y, z+ϵz] + v3t = [x+ϵx, y+ϵy, z+ϵz] + v4t = [x, y+ϵy, z+ϵz] + (v1b, v2b, v3b, v4b, v1t, v2t, v3t, v4t) +end + +""" + rectangle3dpts(o, ϵ) + +Given the origin `o` (3-element vector) and edge lengths `ϵ` (also a +3-element vector) of a 3D rectangle, return the coordinates of its +vertices. +""" +function rectangle3dpts(o, ϵ) + x, y, z = (o...,) + ϵx, ϵy, ϵz = (ϵ...,) + rectangle3dpts(x, y, z, ϵx, ϵy, ϵz) +end + +""" + connectvertices(o, ϵ) + +Given the origin `o` (3-element vector) and edge lengths `ϵ` (also a +3-element vector) of a 3D rectangle, return a vector of line segments +(each a dim-by-n_vertices array) that when plotted together yields +the entire rectangle. +""" +function connectvertices(o, ϵ) + # Get the vertices + v1b, v2b, v3b, v4b, v1t, v2t, v3t, v4t = rectangle3dpts(o, ϵ) + + # Connect vertices in top and bottom planes + connect_top = zeros(3, 5) + connect_bottom = zeros(3, 5) + + # Connect corners of the top and bottom planes + corner1 = zeros(3, 2) + corner2 = zeros(3, 2) + corner3 = zeros(3, 2) + corner4 = zeros(3, 2) + connect_bottom[:, 1:5] = [v1b v2b v3b v4b v1b] + connect_top[:, 1:5]  = [v1t v2t v3t v4t v1t] + corner1[:, 1:2] = [v1b v1t] + corner2[:, 1:2] = [v2b v2t] + corner3[:, 1:2] = [v3b v3t] + corner4[:, 1:2] = [v4b v4t] + + return [connect_top, connect_bottom, + corner1, corner2, corner3, corner4] +end + +""" + splitaxes(x) + +Return a vector of the individual components of an array of points +provided as an array where each point is a column. +""" +splitaxes(x) = ([x[k, :] for k = 1:size(x, 1)]...,) + +""" + plot_3D_rect!(p, origin, edgelengths; + lc = :black, lw = 1.0, ls = :solid) + +Append a a 3D rectangle to an existing plot `p`, given the +origin `o` (3-element vector) and edge lengths +`ϵ` (3-element vector) of the rectangle. +""" +function plot_3D_rect!(p, o, ϵ; + lc = :black, lw = 1.0, ls = :solid, lα = 0.7) + + linesegments = connectvertices(o, ϵ) + for segment in linesegments + plot!(p, splitaxes(segment), + lc = lc, lw = lw, ls = ls, lα = lα) + end +end + +""" +Given the `origin` of a hyperrectangle and the `stepsizes` along each +coordinate axis, divide each axis into `N` stepsizes starting +from `origin[i] + (stepsizes[i]/N)/2` +""" +function fill_hyperrectangle_3D(origin, stepsizes, N = 10) + xs, ys, zs = ([LinRange(origin[i] + stepsizes[i]/N/2, origin[i] + stepsizes[i] - stepsizes[i]/N/2, N) for i = 1:length(origin)]...,) + pts = zeros(Float64, 3, length(xs)^3) + i = 0 + for x in xs + for y in ys + for z in zs + i += 1 + pts[:, i] = [x, y, z] + end + end + end + pts +end diff --git a/src/plot_recipes/helperfunctions_triangulationgrid.jl b/src/plot_recipes/helperfunctions_triangulationgrid.jl new file mode 100644 index 000000000..75910c3c1 --- /dev/null +++ b/src/plot_recipes/helperfunctions_triangulationgrid.jl @@ -0,0 +1,122 @@ +using Simplices: even_sampling_rules +using Simplices: simplexintersection +using Simplices: intersectingvertices + +function interiorpoints(s::SArray{D, T}, n::Int) where {D, T} + R = rand(size(s, 2), n) + convex_coeffs = ((1 ./ sum(R, dims = 1)) .* R) + s * convex_coeffs +end + +function interiorpoints(s, n::Int) + R = rand(size(s, 2), n) + convex_coeffs = ((1 ./ sum(R, dims = 1)) .* R) + s * convex_coeffs +end + +connectvertices(s) = hcat(s, hcat([s[:, i] + for i in [1, 3, 1, 4, 2]]...,)) +splitaxes(x) = ([x[k, :] for k = 1:size(x, 1)]...,) +preparesimplex(s) = splitaxes(connectvertices(s)) + +""" + getsimplex(r::Embeddings.AbstractEmbedding, + t::DelaunayTriangulation, + i::Int) + +Get the vertices of the i-th simplex of the triangulation. +""" +getsimplex(E::Embeddings.AbstractEmbedding, DT::DelaunayTriangulation, i::Int) = + E.points[:, DT.indices[:, i]] + +getsimplex(d::Dataset, DT::DelaunayTriangulation, i::Int) = + d.points[:, DT.indices[:, i]] + +getsimplex(pts, DT::DelaunayTriangulation, i::Int) = + pts[:, DT.indices[:, i]] + +""" + forwardmap(r::Embeddings.AbstractEmbedding, + t::DelaunayTriangulation, + i::Int) + +Get the vertices of the i-th simplex of the triangulation projected +one step forward in time. +""" +forwardmap(r::Embeddings.AbstractEmbedding, t::DelaunayTriangulation, i::Int) = + E.points[:, DT.indices[:, i] .+ 1] + +forwardmap(r::Dataset, DT::DelaunayTriangulation, i::Int) = + E.points[:, DT.indices[:, i] .+ 1] + +forwardmap(pts, DT::DelaunayTriangulation, i::Int) = + pts[:, DT.indices[:, i] .+ 1] + + + +""" +Generate a total of `n_fillpoints` points lying inside the intersection +`sᵢ ∩ sⱼ` between simplices `sᵢ` and `sⱼ`. + +This is done by first performing a Delaunay triangulation of the set of vertices +forming the intersecting volume. Then, points are inserted into the subsimplices +of the convex hull of the intersection. + +If `sample_randomly = true`, then points are generated randomly inside each +subsimplex. If `sample_randomly = false` (default), then points are generated +as the centroids of the subsimplices resulting from a shape-preserving +refinement of the simplices in the convex hull. + +## Arguments +- **``sᵢ``**: An array where each column is a simplex vertex. +- **``sⱼ``**: An array where each column is a simplex vertex. +- **``n_fillpoints``**: The total number of points to fill the convex + hull of the intersection with. If the intersection is more complex than + a single simplex, then the convex hull of the intersection is triangulated + and `n_fillpoints` points are distributed evenly among the resulting + subsimplices. +- **``sample_randomly``**: Should points be sampled randomly within the + intersection? Default is `sample_randomly = false`, which inserts points + inside the intersecting volume in a regular manner. +""" +function generate_fillpoints(sᵢ, sⱼ, n_fillpoints = 100, sample_randomly = true) + + intvol = simplexintersection(sᵢ, sⱼ) + + if intvol > 0 + 1e-9 + intersecting_vertices = Simplices.intersectingvertices(sᵢ, sⱼ) + if size(intersecting_vertices, 1) > 4 + triang = hcat(delaunay(intersecting_vertices)...,) + else + triang = zeros(Int, 4, 1) + triang[:, 1] = 1:4 + end + + n_simplices = size(triang, 2) + if n_simplices == 1 + npts_per_simplex = n_fillpoints + elseif n_simplices > 1 + + npts_per_simplex = ceil(Int, n_fillpoints/n_simplices) + end + C = transpose(subsample_coeffs(3, npts_per_simplex, sample_randomly)) + + fillpts = Vector{Array{Float64, 2}}(undef, n_simplices) + for k = 1:n_simplices + fillpts[k] = transpose(C * intersecting_vertices[triang[:, k], :]) + end + fillpts = hcat(fillpts...,) + return intersecting_vertices, fillpts + else + return nothing + end +end + +export +interiorpoints, +connectvertices, +splitaxes, +preparesimplex, +getsimplex, +forwardmap, +generate_fillpoints diff --git a/src/plot_recipes/recipe_embedding.jl b/src/plot_recipes/recipe_embedding.jl new file mode 100644 index 000000000..f49b33813 --- /dev/null +++ b/src/plot_recipes/recipe_embedding.jl @@ -0,0 +1,14 @@ +@recipe function f(E::StateSpaceReconstruction.AbstractEmbedding; + vars = [1, 2, 3], + mc = :black, ms = 1, lc = :black, lw = 1) + if size(E.points, 1) < 3 + error("Cannot plot reconstruction with < 3 variables.") + end + legend --> false + xlabel --> "x1" + ylabel --> "x2" + zlabel --> "x3" + markercolor --> mc + markersize --> ms + ([E.points[i, :] for i in vars]...,) +end diff --git a/src/plot_recipes/recipe_invariantdistribution.jl b/src/plot_recipes/recipe_invariantdistribution.jl new file mode 100644 index 000000000..0617085e9 --- /dev/null +++ b/src/plot_recipes/recipe_invariantdistribution.jl @@ -0,0 +1,10 @@ +@recipe function f(id::PerronFrobenius.InvariantDistribution) + seriestype := :bar + xlabel --> "state # (i)" + ylabel --> "invariant density" + linecolor --> :black + fillcolor --> :black + fillalpha --> 0.5 + legend --> :none + id.dist +end diff --git a/src/plot_recipes/recipe_rectangularinvariantmeasure.jl b/src/plot_recipes/recipe_rectangularinvariantmeasure.jl new file mode 100644 index 000000000..763164d54 --- /dev/null +++ b/src/plot_recipes/recipe_rectangularinvariantmeasure.jl @@ -0,0 +1,55 @@ +@recipe function f(riv::PerronFrobenius.RectangularInvariantMeasure, + boxfillfactor::Int, linesegments = true; + lw = 0.8, lc = :black, lα = 0.5, ls = :do, + ms = 0.2, mc = :black, mα = 0.3) + + pts = riv.E.points + ϵ = riv.ϵ + axisminima, stepsizes = minima_and_stepsizes(pts, ϵ) + + # All bins are assigned a bin origin. + v = unique(riv.visited_bins_coordinates, dims = 2) + + # Bins may be visited by several times, so we'll get rid + # of the repetitions. + n_visited_boxes = size(v, 2) + D = size(v, 1) + + # Plot the partition grid over the points of the reconstructed + # orbit. + legend --> false + + for i = 1:n_visited_boxes + origin = riv.visited_bins_coordinates[:, i] + # Get the measure of the box, and fill the box with a + # number of points scaled to the measure. + μ = riv.measure.dist[i] + fillpts = fill_hyperrectangle_3D(origin, + stepsizes, + ceil(Int, μ*100)*boxfillfactor) + @series begin + seriestype := :scatter3d + markercolor --> mc + markeralpha --> μ*5 + markersize --> μ*10 + linecolor --> mc + fillcolor --> mc + fillalpha --> μ*2 + splitaxes(fillpts) + end + + if linesegments + linesegments = connectvertices(origin, stepsizes) + + for segment in linesegments + @series begin + seriestype := :path + linealpha --> lα + linewidth --> lw + linecolor --> lc + splitaxes(segment) + end + end + end + end +end diff --git a/src/plot_recipes/recipe_transferoperator.jl b/src/plot_recipes/recipe_transferoperator.jl new file mode 100644 index 000000000..8b2661357 --- /dev/null +++ b/src/plot_recipes/recipe_transferoperator.jl @@ -0,0 +1,6 @@ +@recipe function f(r::PerronFrobenius.RectangularBinningTransferOperator) + seriestype := :heatmap + xlabel --> "state #i" + ylabel --> "state #j" + r.transfermatrix +end diff --git a/src/plot_recipes/recipe_triangulationplot.jl b/src/plot_recipes/recipe_triangulationplot.jl new file mode 100644 index 000000000..1da6f313a --- /dev/null +++ b/src/plot_recipes/recipe_triangulationplot.jl @@ -0,0 +1,303 @@ +include("helperfunctions_triangulationgrid.jl") + +@recipe function plot_triang(original_pts, E::Embeddings.AbstractEmbedding, DT::DelaunayTriangulation; + plot_states = true, + plot_simplices = true, + evenly_subsample = true, + k = 3, # The splitting factor, k^dim new points are introduced per simplex + n_interior_pts = 100, + insert_pts = false, + ms = 0.5, + ms_originalpts = 1.5, + mc = :black, + mc_projection = :red, + mc_originalpts = :black, + mα = 0.5, + edges = true, + vertices = false, + si = 0, # index of the i'th highlighted simplex + sj = 0, # index of the j'th highlighted simplex + ind_projection = 0, #show the simplex corresponding to the forward projection of the `ind_projection`-th simplex + intersecting_vertices = nothing, + plot_intersecting_pts = false, + highlight_ms = 2, + extrapts = [], + lc = :black, + plot_originalpts = true, + plot_triang = true, + plot_imagetriang = true, + plot_fulltriang = false, + lc_fulltriang = :black, + lw_fulltriang = 0.5, + lα_fulltriang = 0.5, + lc_triang = :black, + lc_imagetriang = :black, + lw_imagetriang = 0.5, + lw_triang = 0.5, + ind_excluded_pt = 0, + lw_highlight = 1.2, + markershape_excluded_pt = :star5, + mc_excluded_pt = :black, + ms_excluded_pt = 3, + lα = 0.2, + lw = 0.5, + lc_i = :green, + lc_j = :red, + lw_i = 0.5, + lw_j = 0.5, + lα_i = 1.0, + lα_j = 1.0, + mc_j = :red, + mc_i = :green, + mw_i = 0.5, + mw_j = 0.5, + mα_i = 1.0, + mα_j = 1.0, + lc_projection = :blue, + lw_projection = 2, + lα_projection = 1, + idx_intersect_simplex = si, # the index of the simplex for which to compute the intersection between the simplex and its forward image + mc_intersecting_pts = :black, + markershape_intersecting_pts = :pentagon, + ms_intersecting_pts = 1, + fill_intersecting_volume = true, + n_fillpoints = 500, + ms_fillpoints = 1.0, + mα_fillpoints = 1.0, + mc_fillpoints = :black, + sample_randomly = true + ) + fulltriang = delaunay(StateSpaceReconstruction.embed(original_pts)) + n_simplices_full = size(fulltriang.indices, 2) + n_simplices = size(DT.indices, 2) + n_points = size(E.points, 2) + + # Embedding center + #ce = sum(r.points, 2)/size(r.points, 2) + #lp = r.points[:, end] + + #scatter3d!(p, splitaxes(ce), ms = 3) + #scatter3d!(p, splitaxes(lp), ms = 3, mc = :black) + + legend --> false + xaxis --> false + if length(extrapts) > 0 + @series begin + seriestype := :scatter + legend --> false + mc --> mc + ms --> ms + mα --> mα + splitaxes(extrapts) + end + end + if plot_originalpts + @series begin + seriestype := :scatter + mc --> mc_originalpts + ms --> ms_originalpts + splitaxes(original_pts) + end + if ind_excluded_pt > 0 + @series begin + seriestype := :scatter + mc --> mc_excluded_pt + markershape --> markershape_excluded_pt + ms --> ms_excluded_pt + splitaxes(original_pts[:, ind_excluded_pt]) + end + end + end + + if plot_intersecting_pts && (ind_projection > 0) + sᵢ = getsimplex(original_pts, DT, si) + sⱼ = forwardmap(original_pts, DT, ind_projection) + intvol = simplexintersection(sᵢ, sⱼ) + + if intvol > 0 + 1e-8 + intersecting_vertices, fillpoints = generate_fillpoints(sᵢ, sⱼ, + n_fillpoints = n_fillpoints, + sample_randomly = sample_randomly) + + @series begin + seriestype := :scatter + ms --> ms_intersecting_pts + markershape --> markershape_intersecting_pts + mc --> mc_intersecting_pts + splitaxes(transpose(intersecting_vertices)) + end + @series begin + seriestype := :scatter + ms --> ms_fillpoints + mα --> mα_fillpoints + mc --> mc_fillpoints + lc --> mc_fillpoints + fc --> mc_fillpoints + splitaxes(fillpoints) + end + end + end + + if plot_fulltriang + for i = 1:n_simplices + sᵢ = getsimplex(original_pts, DT, i) + + @series begin + seriestype := :path + lc --> lc_fulltriang + lw --> lw_fulltriang + lα --> lα_fulltriang + legend --> false + preparesimplex(sᵢ) + end + end + end + + for i = 1:n_simplices + sᵢ = getsimplex(original_pts, DT, i) + + if edges + if plot_triang + @series begin + seriestype := :path + lc --> lc_triang + lw --> lw_triang + lα --> lα + legend --> false + preparesimplex(sᵢ) + end + end + + if plot_imagetriang + @series begin + seriestype := :path + lc --> lc_imagetriang + lw --> lw_imagetriang + lα --> lα + ls --> :dash + legend --> false + preparesimplex(forwardmap(original_pts, DT, i)) + end + end + end + + if vertices + if plot_triang + @series begin + seriestype := :scatter + mc --> mc + ms --> ms + mα --> mα + legend --> false + end + end + end + + if insert_pts + if !evenly_subsample + @series begin + mc --> mc + ms --> ms + lα --> lα + legend --> false + random_interior_pts = interiorpoints(sᵢ, n_interior_pts) + splitaxes(random_interior_pts) + end + elseif evenly_subsample + @series begin + mc --> mc + ms --> ms + lα --> lα + legend --> false + splitaxes(sᵢ * even_sampling_rules(3, k)) + end + end + end + end + + if si > 0 + sᵢ = getsimplex(original_pts, DT, si) + + if edges + @series begin + seriestype := :path + label --> "simplex #$sᵢ" + lw --> lw_highlight + lc --> lc_i + lα --> lα_i + preparesimplex(sᵢ) + end + end + + if vertices + @series begin + label --> "simplex #$sᵢ" + ms --> ms_i + mc --> mc_i + preparesimplex(sᵢ) + end + end + end + + if sj > 0 + sⱼ = getsimplex(original_pts, DT, sj) + + if edges + @series begin + seriestype := :path + label --> "simplex #$sⱼ" + lw --> lw_highlight + lc --> lc_j + lα --> lα_j + ls --> :dot + preparesimplex(sⱼ) + end + end + + if vertices + @series begin + label --> "simplex #$sⱼ" + ms --> ms_j + mc --> mc_j + preparesimplex(sⱼ) + end + end + end + + if ind_projection > 0 + @series begin + seriestype := :path + label --> "projection of simplex #$ind_projection" + lw --> lw_highlight + lc --> lc_projection + lα --> lα_projection + preparesimplex(forwardmap(original_pts, DT, ind_projection)) + end + end +end + + +""" +Visualize a set of N points (`original_pts`), a delaunay triangulation +(`DT`) constructed from and embedding (`E`) of the first ``N-1`` points +of the points. By assuming a triangulation of the ``N-1`` first points, +the recipe allows for plotting both the original triangulation and the +image of that triangulation under the forward linear map of its vertices, +one step forward in time. + +## Arguments +- **``original_pts``**: A set of ``N`` points where each column is a point. +- **``E``**: An embedding constructed from the first ``N-1`` points of + `original_pts`. +- **``DT``**: A DelaunayTriangulation instance constructed from `E`. + +## Keyword arguments +- **``plot_triang``**: Plot the original triangulation? +- **``plot_imagetriang``**: Plot the image of the triangulation under the forward linear + map of its vertices? +- **``plot_fulltriang``**: Plot a triangulation of all points in `original_pts`? + +""" +plot_triang + +export plot_triang diff --git a/src/plot_recipes/recipes.jl b/src/plot_recipes/recipes.jl new file mode 100644 index 000000000..8b0655105 --- /dev/null +++ b/src/plot_recipes/recipes.jl @@ -0,0 +1,13 @@ +# Helper functions used in the various recipes. These must be imported first. +include("helperfunctions_rectangulargrid.jl") +include("helperfunctions_triangulationgrid.jl") + +# The actual recipes. +using StateSpaceReconstruction +include("recipe_embedding.jl") + +using PerronFrobenius +include("recipe_rectangularinvariantmeasure.jl") +include("recipe_invariantdistribution.jl") +include("recipe_transferoperator.jl") +include("recipe_triangulationplot.jl") diff --git a/src/plotting/plot_embedding.jl b/src/plotting/plot_embedding.jl deleted file mode 100644 index 9a4a752c0..000000000 --- a/src/plotting/plot_embedding.jl +++ /dev/null @@ -1,26 +0,0 @@ -#################################### -# Plotting -#################################### -@recipe function f(r::AbstractEmbedding; dims = (1, 2, 3)) - pts = r.points - - if dimension(r) > 3 - @warn "Embedding dim > 3, plotting three first axes" - pts = r.points[1:3, :] - end - - if length(dims) > 3 - error("dim = $dim. Can be at most 3.") - end - - if length(dims) == 3 - X = pts[dims[1], :] - Y = pts[dims[2], :] - Z = pts[dims[3], :] - return X, Y, Z - elseif length(dims) == 2 - X = pts[dims[1], :] - Y = pts[dims[2], :] - return X, Y - end -end diff --git a/src/plotting/plot_rectangular.jl b/src/plotting/plot_rectangular.jl deleted file mode 100644 index 9e3ae0bd4..000000000 --- a/src/plotting/plot_rectangular.jl +++ /dev/null @@ -1,189 +0,0 @@ - -######################################################### -######################################################### -# Visualizing different rectangular binnings -######################################################### -######################################################### - -######################################################### -# Helper functions to plot 3D rectangles. -######################################################### -""" - rectangle3dpts(x, y, z, ϵx, ϵy, ϵz) - -Compute the coordinates of the vertices of a 3D rectangle, given the -coordinates of the origin (`x`, `y`, and `z`) and the corresponding edge lengths -(`ϵx`, `ϵy`, `ϵz`). -""" -function rectangle3dpts(x, y, z, ϵx, ϵy, ϵz) - v1b = [x, y, z] - v2b = [x+ϵx, y, z] - v3b = [x+ϵx, y+ϵy, z] - v4b = [x, y+ϵy, z] - v1t = [x, y, z+ϵz] - v2t = [x+ϵx, y, z+ϵz] - v3t = [x+ϵx, y+ϵy, z+ϵz] - v4t = [x, y+ϵy, z+ϵz] - (v1b, v2b, v3b, v4b, v1t, v2t, v3t, v4t) -end - -""" - rectangle3dpts(o, ϵ) - -Given the origin `o` (3-element vector) and edge lengths `ϵ` (also a -3-element vector) of a 3D rectangle, return the coordinates of its -vertices. -""" -function rectangle3dpts(o, ϵ) - x, y, z = (o...,) - ϵx, ϵy, ϵz = (ϵ...,) - rectangle3dpts(x, y, z, ϵx, ϵy, ϵz) -end - -""" - connectvertices(o, ϵ) - -Given the origin `o` (3-element vector) and edge lengths `ϵ` (also a -3-element vector) of a 3D rectangle, return a vector of line segments -(each a dim-by-n_vertices array) that when plotted together yields -the entire rectangle. -""" -function connectvertices(o, ϵ) - # Get the vertices - v1b, v2b, v3b, v4b, v1t, v2t, v3t, v4t = rectangle3dpts(o, ϵ) - - # Connect vertices in top and bottom planes - connect_top = zeros(3, 5) - connect_bottom = zeros(3, 5) - - # Connect corners of the top and bottom planes - corner1 = zeros(3, 2) - corner2 = zeros(3, 2) - corner3 = zeros(3, 2) - corner4 = zeros(3, 2) - connect_bottom[:, 1:5] = [v1b v2b v3b v4b v1b] - connect_top[:, 1:5]  = [v1t v2t v3t v4t v1t] - corner1[:, 1:2] = [v1b v1t] - corner2[:, 1:2] = [v2b v2t] - corner3[:, 1:2] = [v3b v3t] - corner4[:, 1:2] = [v4b v4t] - - return [connect_top, connect_bottom, - corner1, corner2, corner3, corner4] -end - -""" - splitaxes(x) - -Return a vector of the individual components of an array of points -provided as an array where each point is a column. -""" -splitaxes(x) = ([x[k, :] for k = 1:size(x, 1)]...,) - -""" - plot_3D_rect!(p, origin, edgelengths; - lc = :black, lw = 1.0, ls = :solid) - -Append a a 3D rectangle to an existing plot `p`, given the -origin `o` (3-element vector) and edge lengths -`ϵ` (3-element vector) of the rectangle. -""" -function plot_3D_rect!(p, o, ϵ; - lc = :black, lw = 1.0, ls = :solid, lα = 0.7) - - linesegments = connectvertices(o, ϵ) - for segment in linesegments - plot!(p, splitaxes(segment), - lc = lc, lw = lw, ls = ls, lα = lα) - end -end - - -######################################################### -# Given a partition and a binning scheme given by ϵ, -# plot the grid superimposed on the points of the orbit. -######################################################### -""" - plot_partition(pts::AbstractArray{T, 2}, ϵ; - mc = :blue, ms = 1.5, mα = 0.8, - lc = :black, lw = 2, ls = :dash, lα = 0.6) where T - -Partition the space defined by `pts` into rectangular boxes -with a binning scheme controlled by `ϵ`. - -The following `ϵ` will work: - -* `ϵ::Int` divide each axis into `ϵ` intervals of the same size. -* `ϵ::Float` divide each axis into intervals of size `ϵ`. -* `ϵ::Vector{Int}` divide the i-th axis into `ϵᵢ` intervals of the same size. -* `ϵ::Vector{Float64}` divide the i-th axis into intervals of size `ϵᵢ`. - -The points are assumed to be provided as an array where each point is -a column. - -`mc`, `ms` and `mα` control the marker color, marker size and marker opacity, -respectively. `lc`, `lw`, `lα` and `ls` control the line color, line width, -line opacity and line style, respectively. -""" -function plot_partition(pts::AbstractArray{T, 2}, ϵ; - mc = :black, ms = 2, mα = 0.8, - lc = :blue, lw = 1.5, ls = :dash, lα = 0.6) where T - - if size(pts, 1) > size(pts, 2) - #info("Treating each row as a point") - pts = transpose(pts) - end - - axisminima, stepsizes = minima_and_stepsizes(pts, ϵ) - - # Assign bins to the points. - v = assign_coordinate_labels(pts, ϵ) - - # Bins may be visited by several times, so we'll get rid - # of the repetitions. - V = unique(v, dims = 2) - n_visited_boxes = size(V, 2) - - # Plot the partition grid over the points of the reconstructed - # orbit. - p = plot(legend = false) - for i = 1:n_visited_boxes - origin = V[:, i] - plot_3D_rect!(p, origin, stepsizes, lc = lc, ls = ls, lα = lα) - end - scatter3d!(p, splitaxes(pts), ms = ms, mc = mc, mα = mα) - xlabel!("x") - ylabel!("y") - p -end - -""" - plot_partition(E::Embeddings.AbstractEmbedding, ϵ; vars = [1, 2, 3], - mc = :blue, ms = 2, mα = 0.8, - lc = :black, lw = 2, ls = :dash, lα = 0.6) - -Partition the embedding into rectangular boxes with a binning scheme controlled -by `ϵ`. If there are more than three variables in the embedding, you can set -which one to use with the `vars` argument (by default, `vars = [1, 2, 3]`). - -The following `ϵ` will work: - -* `ϵ::Int` divide each axis into `ϵ` intervals of the same size. -* `ϵ::Float` divide each axis into intervals of size `ϵ`. -* `ϵ::Vector{Int}` divide the i-th axis into `ϵᵢ` intervals of the same size. -* `ϵ::Vector{Float64}` divide the i-th axis into intervals of size `ϵᵢ`. - - -`mc`, `ms` and `mα` control the marker color, marker size and marker opacity, -respectively. `lc`, `lw`, `lα` and `ls` control the line color, line width, -line opacity and line style, respectively. -""" -function plot_partition(E::Embeddings.AbstractEmbedding, ϵ; - vars = [1, 2, 3], - mc = :blue, ms = 2, mα = 0.8, - lc = :black, lw = 1.5, ls = :dash, lα = 0.6) - plot_partition(E.points[vars, :], ϵ, mc = mc, ms = ms, - lc = lc, lw = lw, ls = ls) -end - -export plot_partition, plot_3D_rect!, splitaxes, connectvertices diff --git a/src/plotting/plot_triang.jl b/src/plotting/plot_triang.jl deleted file mode 100644 index 24e8458dd..000000000 --- a/src/plotting/plot_triang.jl +++ /dev/null @@ -1,127 +0,0 @@ -using Simplices: even_sampling_rules -using StaticArrays -using RecipesBase - -function interiorpoints(s::SArray{D, T}, n::Int) where {D, T} - R = rand(size(s, 2), n) - convex_coeffs = ((1 ./ sum(R, dims = 1)) .* R) - s * convex_coeffs -end - -connectvertices(s) = hcat(s, hcat([s[:, i] - for i in [1, 3, 1, 4, 2]]...,)) -splitaxes(x) = ([x[k, :] for k = 1:size(x, 1)]...,) -preparesimplex(s) = splitaxes(connectvertices(s)) - -""" - getsimplex(r::Embeddings.AbstractEmbedding, - t::DelaunayTriangulation, - i::Int) - -Get the vertices of the i-th simplex of the triangulation. -""" -getsimplex(r::Embeddings.AbstractEmbedding, t::DelaunayTriangulation, i::Int) = - r.points[:, t.indices[:, i]] - -getsimplex(d::Dataset, t::DelaunayTriangulation, i::Int) = - r.points[:, t.indices[:, i]] -""" - forwardmap(r::Embeddings.AbstractEmbedding, - t::DelaunayTriangulation, - i::Int) - -Get the vertices of the i-th simplex of the triangulation projected -one step forward in time. -""" -forwardmap(r::Embeddings.AbstractEmbedding, t::DelaunayTriangulation, i::Int) = - r.points[:, t.indices[:, i] .+ 1] - -forwardmap(r::Dataset, t::DelaunayTriangulation, i::Int) = - r.points[:, t.indices[:, i] .+ 1] - -""" - plot_triang(r::AbstractReconstruction, t::DelaunayTriangulation) - -Plot the triangulation of a 3D state space reconstruction. -""" -function plot_triang(r::Embeddings.AbstractEmbedding, t::DelaunayTriangulation; - plot_states = true, - plot_simplices = true, - evenly_subsample = true, - k = 3, # The splitting factor, k^dim new points are introduced per simplex - n_interior_pts = 100, - insert_pts = false, - ms = 0.5, mc = :red, - edges = true, - vertices = true, - highlight = 0, - highlight_ms = 2, - extrapts = []) # highlight a given simplex, 0 = no simplices are highlighted - - n_simplices = size(t.indices, 2) - n_points = size(r.points, 2) - p = plot() - - # Embedding center - #ce = sum(r.points, 2)/size(r.points, 2) - #lp = r.points[:, end] - - #scatter3d!(p, splitaxes(ce), ms = 3) - #scatter3d!(p, splitaxes(lp), ms = 3, mc = :black) - - if length(extrapts) > 0 - scatter3d!(p, splitaxes(extrapts), ms = 1.4, mc = :black, mα = 0.5) - end - - for i = 1:n_points - pᵢx = r.points[1, i] - pᵢy = r.points[2, i] - pᵢz = r.points[3, i] - #scatter3d!(p, pᵢx, pᵢy, pᵢz, - # mc = :black, mα = 0.5, ms = 1) - end - - - for i = 1:n_simplices - # Get the ith simplex and connect its vertices - sᵢ = getsimplex(r, t, i) - cᵢ = connectvertices(sᵢ) - - if edges - plot3d!(p, preparesimplex(cᵢ), lc = :black, lw = 0.5, lα = 0.5, legend = false) - end - - if vertices - scatter3d!(p, preparesimplex(sᵢ), mc = mc, ms = ms, mα = 0.5, legend = false) - end - - if insert_pts - if !evenly_subsample - random_interior_pts = interiorpoints(sᵢ, n_interior_pts) - scatter3d!(p, splitaxes(random_interior_pts), - mc = :red, ms = ms, lα = 0.5) - elseif evenly_subsample - scatter3d!(p, splitaxes(sᵢ * even_sampling_rules(3, k)), - mc = :red, ms = ms, lα = 0.5) - end - end - end - # Starting simplex - if length(r) - 1 > highlight> 0 - S_orig = getsimplex(r, t, highlight) - S_image = forwardmap(r, t, highlight) - - if edges - plot3d!(preparesimplex(S_orig), label = "simplex $highlight", lw = 2, lc = :green) - plot3d!(preparesimplex(S_image), label = "forward projection", lw = 2, lc = :red) - end - - if vertices - scatter3d!(preparesimplex(S_orig), label = "simplex $highlight", ms = highlight_ms, mc = :green) - scatter3d!(preparesimplex(S_image), label = "forward projection", ms = highlight_ms, mc = :red) - end - end - p -end - -export plot_triang, interiorpoints, splitaxes, preparesimplex, forwardmap diff --git a/test/plot_recipes.jl b/test/plot_recipes.jl index a5be8b6df..8e0ad288c 100644 --- a/test/plot_recipes.jl +++ b/test/plot_recipes.jl @@ -1,14 +1,63 @@ +using CausalityTools +using RecipesBase -@testset "Plotting recipes" begin - emb = embed([diff(rand(30)) for i = 1:3], [1, 2, 3], [1, 0, -1]) - emb_invariant = invariantize(emb) +@testset "Recipes for subpackages" begin # Test plot recipes by calling RecipesBase.apply_recipe with empty dict. # It should return a vector of RecipesBase.RecipeData d = Dict{Symbol,Any}() - @test typeof(RecipesBase.apply_recipe(d, emb)) == Array{RecipesBase.RecipeData,1} - @test typeof(RecipesBase.apply_recipe(d, emb_invariant)) == Array{RecipesBase.RecipeData,1} + # Create an example embedding + E = StateSpaceReconstruction.embed([diff(rand(30)) for i = 1:3], [1, 2, 3], [1, 0, -1]) + E_invariant = invariantize(E) + invm = RectangularInvariantMeasure(E, [4, 5, 6]) + + + @testset "StateSpaceReconstruction.AbstractEmbedding" begin + @test typeof(RecipesBase.apply_recipe(d, E)) == Array{RecipesBase.RecipeData,1} + @test typeof(RecipesBase.apply_recipe(d, E_invariant)) == Array{RecipesBase.RecipeData,1} + end + + @testset "PerronFrobenius." begin + typeof(RecipesBase.apply_recipe(d, invm.transfermatrix)) == Array{RecipesBase.RecipeData,1} + end + + @testset "Invariant distribution" begin + typeof(RecipesBase.apply_recipe(d, invm.measure)) == Array{RecipesBase.RecipeData,1} + end + + @testset "RectangularInvariantMeasure" begin + boxfillfactor = 3 + # Third argument is linesegment, which can be true or false. It determines + # whether box edges around the partition elements are drawn. + typeof(RecipesBase.apply_recipe(d, (invm, 3, true))) == Array{RecipesBase.RecipeData,1} + typeof(RecipesBase.apply_recipe(d, (invm, 3, false))) == Array{RecipesBase.RecipeData,1} + end + + @testset "Vizualizing triangulation and simplices" begin + pts = rand(3, 15) + n_pts = size(pts, 2) + + # Embed the point + E = StateSpaceReconstruction.embed(pts) + + # Make sure last point is inside the convex hull of the previous points. + # If it is not, move it towards the center of the embedding until it is. + E_linearlyinvariant = invariantize(E) + + # Pick out points and their images (so that they can + # be indexed with the same simplex index) + originalpts = E.points[:, 1:end-1] + forwardpts = E.points[:, 2:end] + + # Re-create the embedding, but exclude the last point. + E = StateSpaceReconstruction.embed(originalpts) + + # Triangulate all points but the last point. + DT = delaunay(E); + + typeof(RecipesBase.apply_recipe(d, (pts, E, DT))) == Array{RecipesBase.RecipeData,1} + end end diff --git a/test/runtests.jl b/test/runtests.jl index d415bbeb9..687fa6d99 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,8 +9,15 @@ import DynamicalSystemsBase: DiscreteDynamicalSystem, ContinuousDynamicalSystem, Dataset + +using StateSpaceReconstruction +using PerronFrobenius using CausalityTools +@testset "Plot recipes" begin + include("plot_recipes.jl") +end + @testset "Discrete systems" begin include("discrete_systems.jl") end