diff --git a/examples/PaperExperiments/double-well-experiment.jl b/examples/PaperExperiments/double-well-experiment.jl index 4de438a..8235779 100644 --- a/examples/PaperExperiments/double-well-experiment.jl +++ b/examples/PaperExperiments/double-well-experiment.jl @@ -1,4 +1,4 @@ -include("../../src/CyclingSignatures-include-file.jl") +include("../../cycling-sigs/CyclingSignatures.jl/src/CyclingSignatures-include-file.jl") using JLD2 @@ -11,6 +11,7 @@ Y_dw = data_doublewell[:,1:end-1] Z_dw = mapslices(normalize, data_doublewell[:,2:end]-data_doublewell[:,1:end-1], dims=1) ts_dw = trajectoryToTrajectorySpaceSB(Y_dw,Z_dw, boxsize_dw, sb_radius_dw, filter_missing=true,shortest_path_fix=true) +saveTrajectorySpace("double-well", ts_dw) n = 1000 dwExperimentParams = map(i->SubsegmentSampleParameter(i,n), 10:10:1000) diff --git a/examples/PaperExperiments/lorenz-experiment.jl b/examples/PaperExperiments/lorenz-experiment.jl index 28ffdc8..155938d 100644 --- a/examples/PaperExperiments/lorenz-experiment.jl +++ b/examples/PaperExperiments/lorenz-experiment.jl @@ -1,4 +1,4 @@ -include("../../src/CyclingSignatures-include-file.jl") +include("../../cycling-sigs/CyclingSignatures.jl/src/CyclingSignatures-include-file.jl") using DynamicalSystems using JLD2 diff --git a/examples/PaperExperiments/lorenz-results.jld2 b/examples/PaperExperiments/lorenz-results.jld2 index 4439d4c..da7d24d 100644 --- a/examples/PaperExperiments/lorenz-results.jld2 +++ b/examples/PaperExperiments/lorenz-results.jld2 @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:9980c9a62b1ec32d02e3918682fd4286595df179c2ce5dab49ba0f6704d20f3c -size 145422155 +oid sha256:e9171be9c29e343accc8533db3ef3566ac966a986ce569d8c0fb648f0a65c55a +size 145375112 diff --git a/examples/Tesis/TrajectorySpaces/burkeshaw_5_3.jld2 b/examples/Tesis/TrajectorySpaces/burkeshaw_5_3.jld2 new file mode 100644 index 0000000..f6b4d86 --- /dev/null +++ b/examples/Tesis/TrajectorySpaces/burkeshaw_5_3.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1d0b72a9062846633a26bdf3500c312a3dd5d3e0c7085245eb3d6887886b9ac3 +size 56333421 diff --git a/examples/Tesis/TrajectorySpaces/burkeshaw_8_3.jld2 b/examples/Tesis/TrajectorySpaces/burkeshaw_8_3.jld2 new file mode 100644 index 0000000..3a4297d --- /dev/null +++ b/examples/Tesis/TrajectorySpaces/burkeshaw_8_3.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:04fa8a4919c2bf4a24b2eb490951dbc6339b391ac8556dca6825cd1c90a7c2f1 +size 56333421 diff --git a/examples/Tesis/TrajectorySpaces/double-well.jld2 b/examples/Tesis/TrajectorySpaces/double-well.jld2 new file mode 100644 index 0000000..f9482d2 --- /dev/null +++ b/examples/Tesis/TrajectorySpaces/double-well.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2969abe37f621daf9a1df63e6bef841cdbf06bbb585ef8c4b2f4d860b4a6e3c4 +size 16646664 diff --git a/examples/Tesis/TrajectorySpaces/lorenz.jld2 b/examples/Tesis/TrajectorySpaces/lorenz.jld2 new file mode 100644 index 0000000..3cf4c0b --- /dev/null +++ b/examples/Tesis/TrajectorySpaces/lorenz.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:de37b391be9f57813964b18bcf415fd4e3de62e52d5bb3dbbe9b48531ae11406 +size 71450095 diff --git a/examples/Tesis/TrajectorySpaces/roessler-three_4_2.jld2 b/examples/Tesis/TrajectorySpaces/roessler-three_4_2.jld2 new file mode 100644 index 0000000..4b61892 --- /dev/null +++ b/examples/Tesis/TrajectorySpaces/roessler-three_4_2.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:972667dc71d0acaaed70e43b98c69081ad616459db882ca9c55824b842b8d90f +size 59022578 diff --git a/examples/Tesis/TrajectorySpaces/roessler-three_5_3.jld2 b/examples/Tesis/TrajectorySpaces/roessler-three_5_3.jld2 new file mode 100644 index 0000000..93fbe4a --- /dev/null +++ b/examples/Tesis/TrajectorySpaces/roessler-three_5_3.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b5e611d772f8187d6cf6a08cd54f5839ab6398e7ba3286bbccdf315440699969 +size 58671428 diff --git a/examples/Tesis/TrajectorySpaces/roessler_4_2.jld2 b/examples/Tesis/TrajectorySpaces/roessler_4_2.jld2 new file mode 100644 index 0000000..5b8e622 --- /dev/null +++ b/examples/Tesis/TrajectorySpaces/roessler_4_2.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9f158e7157fdcefcecda28ea46cd8821564a66ab9449b07edfdf12ba528dba63 +size 57745757 diff --git a/examples/Tesis/TrajectorySpaces/roessler_5_3.jld2 b/examples/Tesis/TrajectorySpaces/roessler_5_3.jld2 new file mode 100644 index 0000000..50bc952 --- /dev/null +++ b/examples/Tesis/TrajectorySpaces/roessler_5_3.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4f2c5e91538aeef425c2e0769ffce83cba26a0f4da314d0ca013250e234b9846 +size 57429111 diff --git a/examples/Tesis/TrajectorySpaces/shilnikov_15_10.jld2 b/examples/Tesis/TrajectorySpaces/shilnikov_15_10.jld2 new file mode 100644 index 0000000..107a729 --- /dev/null +++ b/examples/Tesis/TrajectorySpaces/shilnikov_15_10.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:10f91a045888cb8065d3035ef6dc432332ba824663bb06998b1a289826b89827 +size 62064689 diff --git a/examples/Tesis/TrajectorySpaces/shilnikov_15_9.jld2 b/examples/Tesis/TrajectorySpaces/shilnikov_15_9.jld2 new file mode 100644 index 0000000..344d429 --- /dev/null +++ b/examples/Tesis/TrajectorySpaces/shilnikov_15_9.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:87a31abf1974c3a32608087c2ea872cf48dcb80e9a122da8030a1bd9ecf110c4 +size 61009755 diff --git a/examples/Tesis/Utils/experiments.jl b/examples/Tesis/Utils/experiments.jl new file mode 100644 index 0000000..83afe52 --- /dev/null +++ b/examples/Tesis/Utils/experiments.jl @@ -0,0 +1,33 @@ +function runAndSaveExperimentWithRK(name, Y, Z, r, k; t = nothing) + trajectorySpace = trajectoryToTrajectorySpaceSB(Y, Z, r, k, t_vec = t; filter_missing=true, shortest_path_fix=true) + saveTrajectorySpace("$(name)_$(r)_$(k)", trajectorySpace) + + betti = getBetti(trajectorySpace) + if (betti == 0) + throw(ErrorException("betti number should be greater than 0")) + end + + n = 1000 + experimentParams = map(i->SubsegmentSampleParameter(i,n), 10:10:1000) + + experiments = map(experimentParams) do p + RandomSubsegmentExperiment(trajectorySpace, p, r) + end + + times, results = runExperimentsWithTimer(experiments) + dir_name = "examples/Tesis/" + jldsave(dir_name * "$(name)-results_$(r)_$(k).jld2", results=results, times=times) +end + +function saveTrajectorySpace(name::String, ts) + dir = "examples/Tesis/TrajectorySpaces/" + filename = dir * "$(name).jld2" + @save filename ts +end + +function loadTrajectorySpace(name::String) + dir = "examples/Tesis/TrajectorySpaces/" + filename = dir * "$(name).jld2" + @load filename ts + return ts +end \ No newline at end of file diff --git a/examples/Tesis/Utils/plotting.jl b/examples/Tesis/Utils/plotting.jl new file mode 100644 index 0000000..c4e9698 --- /dev/null +++ b/examples/Tesis/Utils/plotting.jl @@ -0,0 +1,636 @@ +using GeometryBasics +using Colors + +###################################################################################################### +###################################################################################################### +# plot attractor points +function plotPoints(trajectorySpace::TrajectorySpace) + points = getTrajectory(trajectorySpace).trajectory + + dim = getDimension(trajectorySpace) + if (dim == 2) + return plotPoints_2D(points) + end + + return plotPoints_3D(points) +end + +function plotPoints_2D(points::Matrix{Float64}) + fig = Figure() + ax = Axis(fig[1, 1]) + scatter!(ax, points[1, :], points[2, :]; markersize=3, color=:black) + return fig +end + +function plotPoints_3D(points::Matrix{Float64}) + fig = Figure() + ax = Axis3(fig[1, 1]) + scatter!(ax, points[1, :], points[2, :], points[3, :]; markersize=3, color=:black) + return fig +end + +###################################################################################################### +###################################################################################################### +# plot attractor tangent field +function plotTangentField(trajectorySpace::TrajectorySpace) + from = 1 + to = getTrajectorySpaceSize(trajectorySpace) + return plotTangentField(trajectorySpace, from, to) +end + +function plotTangentField(trajectorySpace::TrajectorySpace, from::Int64, to::Int64) + trajectoryPoints = getTrajectorySpaceTrajectoryPoints(trajectorySpace)[:, from:to] + trajectoryVectors = getTrajectorySpaceTrajectoryVectors(trajectorySpace)[:, from:to] + + dim = getDimension(trajectorySpace) + if (dim == 2) + return plotTangentField_2D(trajectoryPoints, trajectoryVectors) + end + return plotTangentField_3D(trajectoryPoints, trajectoryVectors) +end + +function plotTangentField_2D(trajectoryPoints::Matrix{Float64}, trajectoryVectors::Matrix{Float64}) + fig = Figure() + ax = Axis(fig[1, 1]) + + N = size(trajectoryPoints, 2) + origins = [Point2f0(trajectoryPoints[1, i], trajectoryPoints[2, i]) for i in 1:N] + scale = 0.01 + directions = [Vec2f0(trajectoryVectors[1, i], trajectoryVectors[2, i]) * scale for i in 1:N] + + scatter!(ax, trajectoryPoints[1, :], trajectoryPoints[2, :]; markersize=5, color=:black) + arrows!(ax, origins, directions; linewidth=0.5, arrowsize=0.2, color=:blue) + + return fig +end + +function plotTangentField_3D(trajectoryPoints::Matrix{Float64}, trajectoryVectors::Matrix{Float64}) + fig = Figure() + ax = Axis3(fig[1, 1]) + + # plotted points would overlap with plotted vectors so commenting this line + #scatter!(ax, trajectoryPoints[1, :], trajectoryPoints[2, :], trajectoryPoints[3, :]; markersize=5, color=:black) + + N = size(trajectoryPoints, 2) + scale = 0.2 + for i in 1:N + p = Point3f0(trajectoryPoints[1, i], trajectoryPoints[2, i], trajectoryPoints[3, i]) + v = Vec3f0(trajectoryVectors[1, i], trajectoryVectors[2, i], trajectoryVectors[3, i]) * scale + segment = [p, p + v] + lines!(ax, segment; color=:blue, linewidth=1) + end + + return fig +end + +###################################################################################################### +###################################################################################################### +# plot simplified comparison space with or without grid boxes +function plotSimplifiedComparisonSpace(trajectorySpace::TrajectorySpace) + dim = getDimension(trajectorySpace) + if (dim == 2) + return plotSimplifiedComparisonSpace_2D(trajectorySpace) + end + + return plotSimplifiedComparisonSpace_3D(trajectorySpace) +end + +function plotSimplifiedComparisonSpace_2D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis(fig[1, 1]) + + boxSpace = getBoxSpace(trajectorySpace) + comparisonSpacePoints = getPlotSimplifiedPoints(boxSpace) + scatter!(ax, comparisonSpacePoints[1, :], comparisonSpacePoints[2, :]; markersize=10, color=:black) + return fig +end + +function plotSimplifiedComparisonSpace_3D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis3(fig[1, 1]) + + boxSpace = getBoxSpace(trajectorySpace) + comparisonSpacePoints = getPlotSimplifiedPoints(boxSpace) + scatter!(ax, comparisonSpacePoints[1, :], comparisonSpacePoints[2, :], comparisonSpacePoints[3, :]; markersize=12, color=:black) + return fig +end + +function plotSimplifiedComparisonSpaceWithGridBoxes(trajectorySpace::TrajectorySpace) + dim = getDimension(trajectorySpace) + if (dim == 2) + return plotSimplifiedComparisonSpaceWithGridBoxes_2D(trajectorySpace) + end + + return plotSimplifiedComparisonSpaceWithGridBoxes_3D(trajectorySpace) +end + +function plotSimplifiedComparisonSpaceWithGridBoxes_2D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis(fig[1, 1]) + + boxSpace = getBoxSpace(trajectorySpace) + + # plot grid boxes + gridPoints = boxSpace.X + boxSize = 1 # hardcoding this value so its better scaled in the plot + half_box = boxSize / 2 + for i in 1:size(gridPoints, 2) + xc, yc = gridPoints[:, i] + rect = Rect(xc - half_box, yc - half_box, boxSize, boxSize) + poly!(ax, rect, color=:gray90, strokecolor=:gray60, strokewidth=1) + + x0, y0 = xc - half_box, yc - half_box + x1, y1 = xc + half_box, yc + half_box + edges_2d = [ + (x0, y0), (x1, y0), (x1, y1), (x0, y1), (x0, y0) + ] + lines!(ax, Point2f.(edges_2d), color=:gray30, linewidth=1) + end + + # plot comparison space points + comparisonSpacePoints = getPlotSimplifiedPoints(boxSpace) + scatter!(ax, comparisonSpacePoints[1, :], comparisonSpacePoints[2, :]; markersize=5, color=:black) + return fig +end + +function plotSimplifiedComparisonSpaceWithGridBoxes_3D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis3(fig[1, 1]) + + boxSpace = getBoxSpace(trajectorySpace) + # plot comparison space points + comparisonSpacePoints = getPlotSimplifiedPoints(boxSpace) + scatter!(ax, comparisonSpacePoints[1, :], comparisonSpacePoints[2, :], comparisonSpacePoints[3, :]; markersize=8, color=:black) + + # plot grid cubes + gridPoints = boxSpace.X + cube_size = 0.5 # hardcoding this value so its better scaled in the plot + half_cube = cube_size / 2 + for i in 1:size(gridPoints, 2) + xc, yc, zc = gridPoints[:, i] + cubo = FRect3D(Point3f(xc - half_cube, yc - half_cube, zc - half_cube), + Vec3f(cube_size, cube_size, cube_size)) + mesh!(ax, cubo, color=(:gray90, 0.2), transparency=true, shading=false) + x0, y0, z0 = xc - half_cube, yc - half_cube, zc - half_cube + x1, y1, z1 = xc + half_cube, yc + half_cube, zc + half_cube + edges = [ + (x0,y0,z0), (x1,y0,z0), (x1,y1,z0), (x0,y1,z0), (x0,y0,z0), + (x0,y0,z1), (x1,y0,z1), (x1,y1,z1), (x0,y1,z1), (x0,y0,z1), + (x0,y0,z0), (x0,y0,z1), + (x1,y0,z0), (x1,y0,z1), + (x1,y1,z0), (x1,y1,z1), + (x0,y1,z0), (x0,y1,z1) + ] + lines!(ax, Point3f.(edges), color=:gray30, linewidth=1) + end + + return fig +end + +###################################################################################################### +###################################################################################################### +# plot comparison space +function plotComparisonSpace(trajectorySpace::TrajectorySpace) + dim = getDimension(trajectorySpace) + + if (dim == 2) + return plotComparisonSpace_2D(trajectorySpace) + end + + return plotComparisonSpace_3D(trajectorySpace) +end + +function plotComparisonSpace_2D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis(fig[1, 1]) + + boxSpace = getBoxSpace(trajectorySpace) + comparisonSpacePoints = getPlotPoints(boxSpace) + scatter!(ax, comparisonSpacePoints[1, :], comparisonSpacePoints[2, :]; markersize=10, color=:black) + return fig +end + +function plotComparisonSpace_3D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis3(fig[1, 1]) + + boxSpace = getBoxSpace(trajectorySpace) + comparisonSpacePoints = getPlotPoints(boxSpace) + scatter!(ax, comparisonSpacePoints[1, :], comparisonSpacePoints[2, :], comparisonSpacePoints[3, :]; markersize=10, color=:black) + return fig +end + +function plotComparisonSpaceWithGridBoxes(trajectorySpace::TrajectorySpace) + dim = getDimension(trajectorySpace) + + if (dim == 2) + return plotComparisonSpaceWithGridBoxes_2D(trajectorySpace) + end + + return plotComparisonSpaceWithGridBoxes_3D(trajectorySpace) +end + +function plotComparisonSpaceWithGridBoxes_2D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis(fig[1, 1]) + + boxSpace = getBoxSpace(trajectorySpace) + + # plot grid boxes + gridPoints = boxSpace.X + gridVectors = boxSpace.VF + boxSize = boxSpace.boxsize + sbSize = boxSpace.sphereBundleRadius + half_box = boxSize / 2 + for i in 1:size(gridPoints, 2) + xc, yc = gridPoints[:, i] + # rect = Rect(xc - half_box, yc - half_box, boxSize, boxSize) + # poly!(ax, rect, color=:gray90, strokecolor=:gray60, strokewidth=1) + + x0, y0 = xc - half_box, yc - half_box + x1, y1 = xc + half_box, yc + half_box + x .+ offset .* y + edges_2d = [ + (x0, y0), (x1, y0), (x1, y1), (x0, y1), (x0, y0) + ] + lines!(ax, Point2f.(edges_2d), color=:gray30, linewidth=1) + end + + # plot comparison space points + comparisonSpacePoints = getPlotPoints(boxSpace) + scatter!(ax, comparisonSpacePoints[1, :], comparisonSpacePoints[2, :]; markersize=10, color=:black) + return fig +end + +function plotComparisonSpaceWithGridBoxes_3D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis3(fig[1, 1]) + + boxSpace = getBoxSpace(trajectorySpace) + + # plot grid cubes + gridPoints = boxSpace.X + cube_size = boxSpace.boxsize * .1 + half_cube = cube_size / 2 + for i in 1:size(gridPoints, 2) + xc, yc, zc = gridPoints[:, i] + cubo = FRect3D(Point3f(xc - half_cube, yc - half_cube, zc - half_cube), + Vec3f(cube_size, cube_size, cube_size)) + mesh!(ax, cubo, color=(:gray90, 0.2), transparency=true, shading=false) + x0, y0, z0 = xc - half_cube, yc - half_cube, zc - half_cube + x1, y1, z1 = xc + half_cube, yc + half_cube, zc + half_cube + edges = [ + (x0,y0,z0), (x1,y0,z0), (x1,y1,z0), (x0,y1,z0), (x0,y0,z0), + (x0,y0,z1), (x1,y0,z1), (x1,y1,z1), (x0,y1,z1), (x0,y0,z1), + (x0,y0,z0), (x0,y0,z1), + (x1,y0,z0), (x1,y0,z1), + (x1,y1,z0), (x1,y1,z1), + (x0,y1,z0), (x0,y1,z1) + ] + lines!(ax, Point3f.(edges), color=:gray30, linewidth=1) + end + + # plot comparison space points + comparisonSpacePoints = getPlotPoints(boxSpace) + scatter!(ax, comparisonSpacePoints[1, :], comparisonSpacePoints[2, :], comparisonSpacePoints[3, :]; markersize=8, color=:black) + return fig +end + +###################################################################################################### +###################################################################################################### +# plot comparison space tangent field with or without grid boxes +function plotComparisonSpaceTangentField(trajectorySpace::TrajectorySpace) + dim = getDimension(trajectorySpace) + + if (dim == 2) + return plotComparisonSpaceTangentField_2D(trajectorySpace) + end + + return plotComparisonSpaceTangentField_3D(trajectorySpace) +end + +function plotComparisonSpaceTangentField_2D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis(fig[1, 1]) + + sphereBundleVectorField = getSphereBundleTangentField(trajectorySpace) + + N = size(sphereBundleVectorField, 2) + sphereBundlePoints = sphereBundleVectorField[1:2, :] + sphereBundleVectors = sphereBundleVectorField[3:4, :] + origins = [Point2f0(sphereBundlePoints[1, i], sphereBundlePoints[2, i]) for i in 1:N] + scale = 0.1 + directions = [Vec2f0(sphereBundleVectors[1, i], sphereBundleVectors[2, i]) * scale for i in 1:N] + + scatter!(ax, sphereBundlePoints[1, :], sphereBundlePoints[2, :]; markersize=5, color=:black) + arrows!(ax, origins, directions; linewidth=0.5, arrowsize=0.2, color=:blue) + return fig +end + +function plotComparisonSpaceTangentField_3D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis3(fig[1, 1]) + + sphereBundleVectorField = getSphereBundleTangentField(trajectorySpace) + sphereBundlePoints = sphereBundleVectorField[1:3, :] + sphereBundleVectors = sphereBundleVectorField[4:6, :] + + N = size(sphereBundlePoints, 2) + scale = 0.05 + for i in 1:N + origin = Point3f0(sphereBundlePoints[1, i], sphereBundlePoints[2, i], sphereBundlePoints[3, i]) + direction = Vec3f0(sphereBundleVectors[1, i], sphereBundleVectors[2, i], sphereBundleVectors[3, i]) * scale + segment = [origin, origin + direction] + lines!(ax, segment; color=:blue, linewidth=2) + end + + return fig +end + +################ +# with grid boxes +function plotComparisonSpaceTangentFieldWithGridBoxes(trajectorySpace::TrajectorySpace) + dim = getDimension(trajectorySpace) + + if (dim == 2) + return plotComparisonSpaceTangentFieldWithGridBoxes_2D(trajectorySpace) + end + + return plotComparisonSpaceTangentFieldWithGridBoxes_3D(trajectorySpace) +end + +function plotComparisonSpaceTangentFieldWithGridBoxes_2D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis(fig[1, 1]) + + boxSpace = getBoxSpace(trajectorySpace) + # plot grid boxes + gridPoints = boxSpace.X + box_size = boxSpace.boxsize + half_box = box_size / 2 + for i in 1:size(gridPoints, 2) + xc, yc = gridPoints[:, i] + rect = Rect(xc - half_box, yc - half_box, box_size, box_size) + poly!(ax, rect, color=:gray90, strokecolor=:gray60, strokewidth=1) + + x0, y0 = xc - half_box, yc - half_box + x1, y1 = xc + half_box, yc + half_box + edges_2d = [ + (x0, y0), (x1, y0), (x1, y1), (x0, y1), (x0, y0) + ] + lines!(ax, Point2f.(edges_2d), color=:gray30, linewidth=1) + end + + sphereBundleVectorField = getSphereBundleTangentField(trajectorySpace) + + N = size(sphereBundleVectorField, 2) + sphereBundlePoints = sphereBundleVectorField[1:2, :] + sphereBundleVectors = sphereBundleVectorField[3:4, :] + origins = [Point2f0(sphereBundlePoints[1, i], sphereBundlePoints[2, i]) for i in 1:N] + scale = 0.1 + directions = [Vec2f0(sphereBundleVectors[1, i], sphereBundleVectors[2, i]) * scale for i in 1:N] + + scatter!(ax, sphereBundlePoints[1, :], sphereBundlePoints[2, :]; markersize=5, color=:black) + arrows!(ax, origins, directions; linewidth=0.5, arrowsize=0.2, color=:blue) + return fig +end + +# sin grid boxes probb porque no se ve nada si no +function plotComparisonSpaceTangentFieldWithGridBoxes_3D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis3(fig[1, 1]) + + # plot grid cubes + boxSpace = getBoxSpace(trajectorySpace) + gridPoints = boxSpace.X + cube_size = boxSpace.boxsize + half_cube = cube_size / 2 + for i in 1:size(gridPoints, 2) + xc, yc, zc = gridPoints[:, i] + cubo = FRect3D(Point3f(xc - half_cube, yc - half_cube, zc - half_cube), + Vec3f(cube_size, cube_size, cube_size)) + mesh!(ax, cubo, color=(:gray90, 0.2), transparency=true, shading=false) + x0, y0, z0 = xc - half_cube, yc - half_cube, zc - half_cube + x1, y1, z1 = xc + half_cube, yc + half_cube, zc + half_cube + edges = [ + (x0,y0,z0), (x1,y0,z0), (x1,y1,z0), (x0,y1,z0), (x0,y0,z0), + (x0,y0,z1), (x1,y0,z1), (x1,y1,z1), (x0,y1,z1), (x0,y0,z1), + (x0,y0,z0), (x0,y0,z1), + (x1,y0,z0), (x1,y0,z1), + (x1,y1,z0), (x1,y1,z1), + (x0,y1,z0), (x0,y1,z1) + ] + lines!(ax, Point3f.(edges), color=:gray30, linewidth=1) + end + + sphereBundleVectorField = getSphereBundleTangentField(trajectorySpace) + sphereBundlePoints = sphereBundleVectorField[1:3, :] + sphereBundleVectors = sphereBundleVectorField[4:6, :] + + N = size(sphereBundlePoints, 2) + scale = 0.05 + for i in 1:N + origin = Point3f0(sphereBundlePoints[1, i], sphereBundlePoints[2, i], sphereBundlePoints[3, i]) + direction = Vec3f0(sphereBundleVectors[1, i], sphereBundleVectors[2, i], sphereBundleVectors[3, i]) * scale + segment = [origin, origin + direction] + lines!(ax, segment; color=:blue, linewidth=3.0) + end + + return fig +end + +###################################################################################################### +###################################################################################################### +# plot cycling signatures +function plotCyclingSignatures(results::Vector, rthickening::Real) + trajSpace = getTrajectorySpace(results[1].experiment) + dim = getDimension(trajSpace) + + if (dim == 2) + return plotCyclingSignatures_2D(results, rthickening) + end + + return plotCyclingSignatures_3D(results, rthickening) +end + +function plotCyclingSignatures_2D(results::Vector, rthickening::Real) + fig = Figure() + ax = Axis(fig[1,1]) + + traj = getTrajectory(getTrajectorySpace(results[1].experiment)) + npoints = min(size(traj.trajectory, 2), 50000) + background_pts = get(traj, 1:npoints)[1:2,:] + scatter!(ax, background_pts[1, :], background_pts[2, :], markersize=0.5, color=:gray50) + + sigs = getSigsWithEnoughAppearences(results, rthickening) + println("Signatures with enough appearences: ", [vec(sig) for sig in sigs]) + + amountOfSigs = size(sigs, 1) + colors = [:blue, :red, :green, :yellow, :purple] + for j in 1:amountOfSigs + customColor = colors[j] + sig = sigs[j] + try + signatureRanges = getRepresentativeSignatureRange(results, sig, rthickening) + segment = get(traj, signatureRanges[1]) + segmentPoints = segment[1:2,:] + segmentVectors = segment[3:4,:] + + N = size(segment, 2) + scale = 0.03 + for i in 1:N + p = Point2f0(segmentPoints[1, i], segmentPoints[2, i]) + v = Vec2f0(segmentVectors[1, i], segmentVectors[2, i]) * scale + arrow = [p, p + v] + lines!(ax, arrow; color = customColor, linewidth=3) + end + catch e + println(e) + end + end + + return fig +end + +function plotCyclingSignatures_3D(results::Vector, rthickening::Int64) + fig = Figure() + ax = Axis3(fig[1,1]) + + traj = getTrajectory(getTrajectorySpace(results[1].experiment)) + npoints = min(size(traj.trajectory, 2), 100000) # 100000 points is more than enough for the time series + background_pts = get(traj, 1:npoints)[1:3,:] + scatter!(ax, background_pts[1, :], background_pts[2, :], background_pts[3, :], markersize=0.3, color=:gray80) + + sigs = getSigsWithEnoughAppearences(results, rthickening) + println("Signatures with enough appearences: ", [vec(sig) for sig in sigs]) + + amountOfSigs = size(sigs, 1) + colors = [:blue, :red, :green, :yellow, :purple] + for j in 1:amountOfSigs + customColor = colors[j] + sig = sigs[j] + try + signatureRanges = getRepresentativeSignatureRange(results, sig, rthickening) + segment = get(traj, signatureRanges[1]) + segmentPoints = segment[1:3,:] + vectors = segment[4:6,:] + + N = size(segment, 2) + scale = 0.05 + for i in 1:N + p = Point3f0(segmentPoints[1, i], segmentPoints[2, i], segmentPoints[3, i]) + v = Vec3f0(vectors[1, i], vectors[2, i], vectors[3, i]) * scale + arrow = [p, p + v] + lines!(ax, arrow; color = customColor, linewidth=5) + end + catch e + println(e) + end + end + + return fig +end + +function hasBirthVectors(subsegmentResult::SubsegmentResultReduced) + emptyVector = [] + birthVectors = subsegmentResult.birthVectors + for birthVector in birthVectors + if birthVector != emptyVector + return true + end + end + return false +end + +###################################################################################################### +###################################################################################################### +# plot comparison space complex +function plotComparisonSpaceComplex_2D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis(fig[1, 1]) + + boxSpace = getBoxSpace(trajectorySpace) + plotPoints = getPlotPoints(boxSpace) + + SB_pts = getSphereBundleVectorField(trajectorySpace) + SB_pts_X = SB_pts[1:2,:] + scatter!(ax, SB_pts_X[1, :], SB_pts_X[2, :], markersize=10, color=:red) + scatter!(ax, plotPoints[1, :], plotPoints[2, :], markersize=10, color=:blue) + + dm_pts = sphereBundleDistanceMatrix(SB_pts) + cplx, h1_gen, h1 = boxSpaceH1Info(dm_pts) + + triangles = cplx.cells[2] + plotTriangles(SB_pts_X, triangles) + edges = cplx.cells[1] + plotEdges(SB_pts_X, edges) + + return fig +end + +function plotEdges(sphereBundlePoints::Matrix{Int64}, edges::Vector{Simplex}) + already_used = [] + + for s in edges + i, j = s.vertices[1], s.vertices[2] + p1 = Point2f0(sphereBundlePoints[1, i], sphereBundlePoints[2, i]) + p2 = Point2f0(sphereBundlePoints[1, j], sphereBundlePoints[2, j]) + + edge = sort([p1, p2]) + if !(edge in already_used) + linesegments!(ax, [p1, p2]; color=:blue, linewidth=1.0) + push!(already_used, edge) + end + end +end + +function plotTriangles(sphereBundlePoints::Matrix{Int64}, triangles::Vector{Simplex}) + already_used = [] + + for s in triangles + i, j, k = s.vertices[1], s.vertices[2], s.vertices[3] + p1 = Point2f0(sphereBundlePoints[1, i], sphereBundlePoints[2, i]) + p2 = Point2f0(sphereBundlePoints[1, j], sphereBundlePoints[2, j]) + p3 = Point2f0(sphereBundlePoints[1, k], sphereBundlePoints[2, k]) + + triangle = sort([p1, p2, p3]) # ordena para evitar duplicados por permutación + if !(triangle in already_used) + poly!(ax, [p1, p2, p3]; color=RGBA(0.8, 0.9, 1.0, 0.4), strokewidth=0) + push!(already_used, triangle) + end + end +end + +function plotComparisonSpaceComplex_3D(trajectorySpace::TrajectorySpace) + fig = Figure() + ax = Axis3(fig[1, 1]) + SB_pts = getSphereBundleVectorField(trajectorySpace) + SB_pts_X = SB_pts[1:3,:] + + dm_pts = sphereBundleDistanceMatrix(SB_pts) + cplx, h1_gen, h1 = boxSpaceH1Info(dm_pts) + + # positions = [Point3f0(SB_pts_X[:, i]) for i in 1:size(SB_pts_X, 2)] + # faces = hcat([t.vertices for t in triangles]...) .- 1 + + # mesh!(ax, positions, faces; color=RGBA(0.5, 0.5, 1.0, 0.3)) + triangles = cplx.cells[2] + for t in triangles + i, j, k = t.vertices[1], t.vertices[2], t.vertices[3] + p1 = Point3f0(SB_pts_X[1, i], SB_pts_X[2, i], SB_pts_X[3, i]) + p2 = Point3f0(SB_pts_X[1, j], SB_pts_X[2, j], SB_pts_X[3, j]) + p3 = Point3f0(SB_pts_X[1, k], SB_pts_X[2, k], SB_pts_X[3, k]) + poly!(ax, [p1, p2, p3]; color=RGBA(0.5, 0.5, 1.0, 0.3), strokewidth=0) + println("plot") + end + + edges = cplx.cells[1] + for s in edges + i, j = s.vertices[1], s.vertices[2] + p1 = Point3f0(SB_pts_X[1, i], SB_pts_X[2, i], SB_pts_X[3, i]) + p2 = Point3f0(SB_pts_X[1, j], SB_pts_X[2, j], SB_pts_X[3, j]) + linesegments!(ax, [p1, p2]; color=:black, linewidth=1.0) + println("plotline") + end + return fig +end \ No newline at end of file diff --git a/examples/Tesis/Utils/utils.jl b/examples/Tesis/Utils/utils.jl new file mode 100644 index 0000000..f9774a9 --- /dev/null +++ b/examples/Tesis/Utils/utils.jl @@ -0,0 +1,92 @@ +function getTrajectorySpaceVectorField(trajectorySpace::TrajectorySpace) + return getTrajectory(trajectorySpace).trajectory +end + +function getTrajectorySpaceTrajectoryPoints(trajectorySpace::TrajectorySpace) + trajectorySpaceVectorField = getTrajectorySpaceVectorField(trajectorySpace) + dim = getDimension(trajectorySpace) + + return trajectorySpaceVectorField[1:dim, :] +end + +function getTrajectorySpaceTrajectoryVectors(trajectorySpace::TrajectorySpace) + trajectorySpaceVectorField = getTrajectorySpaceVectorField(trajectorySpace) + dim = getDimension(trajectorySpace) + + return trajectorySpaceVectorField[dim+1:dim*2, :] +end + +function getDimension(trajectorySpace::TrajectorySpace) + trajectorySpaceVectorField = getTrajectorySpaceVectorField(trajectorySpace) + + dim = size(trajectorySpaceVectorField, 1) + if (mod(dim, 3) == 0) + return 3 + end + return 2 +end + +function getTrajectorySpaceSize(trajectorySpace::TrajectorySpace) + trajectorySpaceVectorField = getTrajectorySpaceVectorField(trajectorySpace) + return size(trajectorySpaceVectorField, 2) +end + +function getSphereBundleTangentField(trajectorySpace::TrajectorySpace) + trajectoryPoints = getTrajectorySpaceTrajectoryPoints(trajectorySpace) + trajectoryVectors = getTrajectorySpaceTrajectoryVectors(trajectorySpace) + + boxSpace = getBoxSpace(trajectorySpace) + boxSize = boxSpace.boxsize + sb_radius = boxSpace.sphereBundleRadius + + SB_pts = unique([quantize(trajectoryPoints, boxSize);quantizeSB(trajectoryVectors, sb_radius)],dims=2) + return SB_pts +end + +function getDistanceMatrix(trajectorySpace::TrajectorySpace) + SB_pts = getSphereBundleTangentField(trajectorySpace) + return sphereBundleDistanceMatrix(SB_pts) +end + +function getComplex(trajectorySpace::TrajectorySpace) + dm_pts = getDistanceMatrix(trajectorySpace) + cplx, h1_gen, h1 = boxSpaceH1Info(dm_pts) + return cplx +end + +function getBetti(trajectorySpace::TrajectorySpace) + boxSpace = getBoxSpace(trajectorySpace) + return betti_1(boxSpace.h1) +end + +function getSigsWithEnoughAppearences(results::Vector, rthickening::Real) + cms = subspaceDistribution.(results, 1, rthickening) + allkeys = union(keys.(cms)...) + sigs = collect(allkeys) + println("All signatures: ", [vec(sig) for sig in sigs]) + + filter!(sigs) do k + sum(cm-> get(cm,k,0), cms) >= 20 + end + return sigs +end + +# getting the largest representative segment +function getRepresentativeSignatureRange(results::Vector, sig, r) + println("Getting representative segment for sig: ", vec(sig)) + + resultsSize = size(results, 1) + sigCounter = 0 + rep = nothing + for i in 1:resultsSize + range = results[i] + if (hasBirthVectors(range)) + sigRanges = getSignatureRanges(range, sig, r) + if (sigRanges != UnitRange[]) + sigCounter += 1 + rep = sigRanges + end + end + end + return rep +end \ No newline at end of file diff --git a/examples/Tesis/burkeshaw-experiment.jl b/examples/Tesis/burkeshaw-experiment.jl new file mode 100644 index 0000000..43f045a --- /dev/null +++ b/examples/Tesis/burkeshaw-experiment.jl @@ -0,0 +1,53 @@ +include("../../cycling-sigs/CyclingSignatures.jl/src/CyclingSignatures-include-file.jl") + +using DynamicalSystems +using JLD2 + +function burke_shaw(u0=[1, -2, 0.1]; a = 10.0, b = 4.271, c = 0.0) + return CoupledODEs(burke_shaw_rule, u0, [a, b]) +end +@inbounds function burke_shaw_rule(u, p, t) + du1 = -p[1] * (u[1] + u[2]) + du2 = -u[2] - p[1] * u[1] * u[3] + du3 = p[2] + p[1] * u[1] * u[2] + return SVector{3}(du1, du2, du3) +end + +function burkeShawTrajForCycling() + ds = burke_shaw() + dt = 0.1 + sol,_ = trajectory(ds, 10000, Δt=0.01) + + Y_init = Matrix(sol[:,:])'[:,100:end] + resample_boxsize = 1.0 + resample_sb_radius = 3.0 + burkeshaw_f(v) = Vector(ds.integ.f(v, ds.integ.p, 0)) + + Y_res, t_vec = resampleToConsistent(ds, Y_init, resample_boxsize, dt, sb_radius = resample_sb_radius, sb_fct=burkeshaw_f) + Z_res = mapslices(v->normalize(burkeshaw_f(v),2 ), Y_res, dims=[1]) + return Y_res, Z_res, t_vec +end + +Y_burkeshaw, Z_burkeshaw, t_burkeshaw = burkeShawTrajForCycling() + +name = "burkeshaw" +r=8 +k=3 +runAndSaveExperimentWithRK("burkeshaw", Y_burkeshaw, Z_burkeshaw, r, k, t = t_burkeshaw) + +ts_burkeshaw_8_3 = trajectoryToTrajectorySpaceSB(Y_burkeshaw,Z_burkeshaw, r, k, filter_missing=true,shortest_path_fix=true) +saveTrajectorySpace("burkeshaw_8_3", ts_burkeshaw_8_3) + +burkeshawResults_8_3 = SubsegmentResultReduced.(load(dir_name*"burkeshaw-results_8_3.jld2", "results")); +burkeshaw_r = 6 +sig1_l,_ = subspaceFrequencyMatrix(burkeshawResults_8_3, 1, burkeshaw_r) + +r=5 +k=3 +ts_burkeshaw_5_3 = trajectoryToTrajectorySpaceSB(Y_burkeshaw,Z_burkeshaw, r, k, filter_missing=true,shortest_path_fix=true) +saveTrajectorySpace("burkeshaw_5_3", ts_burkeshaw_5_3) + +runAndSaveExperimentWithRK("burkeshaw", Y_burkeshaw, Z_burkeshaw, r, k, t = t_burkeshaw) +burkeshawResults_5_3 = SubsegmentResultReduced.(load(dir_name*"burkeshaw-results_5_3.jld2", "results")); +burkeshaw_r = 4 +sig1_l,_ = subspaceFrequencyMatrix(burkeshawResults_5_3, 1, burkeshaw_r) diff --git a/examples/Tesis/burkeshaw-results_5_3.jld2 b/examples/Tesis/burkeshaw-results_5_3.jld2 new file mode 100644 index 0000000..9e1f593 --- /dev/null +++ b/examples/Tesis/burkeshaw-results_5_3.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6471b4672c9e0814a330ea5775f5690b09337ec2f68fcbf6eee35dee1b162892 +size 133339056 diff --git a/examples/Tesis/burkeshaw-results_8_3.jld2 b/examples/Tesis/burkeshaw-results_8_3.jld2 new file mode 100644 index 0000000..7b955a8 --- /dev/null +++ b/examples/Tesis/burkeshaw-results_8_3.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7bb16c862343ec52b327c5c79acf557013756aeee4fc0210c63399e3fe9892f4 +size 133349652 diff --git a/examples/Tesis/roessler-experiment.jl b/examples/Tesis/roessler-experiment.jl new file mode 100644 index 0000000..bd4e0fa --- /dev/null +++ b/examples/Tesis/roessler-experiment.jl @@ -0,0 +1,48 @@ +include("../../cycling-sigs/CyclingSignatures.jl/src/CyclingSignatures-include-file.jl") + +using DynamicalSystems + +function roessler(u0=[1.0, -2.0, 0.1]; a = 0.432, b = 2.0, c = 4.0) + return CoupledODEs(roessler_rule, u0, [a, b, c]) +end +@inbounds function roessler_rule(u, p, t) + du1 = -u[2]-u[3] + du2 = u[1] + p[1]*u[2] + du3 = p[2] + u[3]*(u[1] - p[3]) + return SVector{3}(du1, du2, du3) +end + +function roesslerTrajForCycling() + ds = roessler() + dt = 0.1 + sol,_ = trajectory(ds, 10000, Δt=0.01) + + Y_init = Matrix(sol[:,:])'[:,100:end] + resample_boxsize = 1.0 + resample_sb_radius = 3.0 + roessler_f(v) = Vector(ds.integ.f(v, ds.integ.p, 0)) + + Y_res, t_vec = resampleToConsistent(ds, Y_init, resample_boxsize, dt, sb_radius = resample_sb_radius, sb_fct=roessler_f) + Z_res = mapslices(v->normalize(roessler_f(v),2 ), Y_res, dims=[1]) + return Y_res, Z_res, t_vec +end + +Y_roessler, Z_roessler, t_roessler = roesslerTrajForCycling() + +r=4 +k=2 +runAndSaveExperimentWithRK("roessler", Y_roessler, Z_roessler, r, k, t = t_roessler) + + +r=5 +k=3 +runAndSaveExperimentWithRK("roessler", Y_roessler, Z_roessler, r, k, t = t_roessler) + + + +roesslerResults_4_2 = SubsegmentResultReduced.(load(dir_name*"roessler-results_4_2.jld2", "roesslerResults")); +roessler_r = 1 +sig1_l,_ = subspaceFrequencyMatrix(roesslerResults_4_2, 1, roessler_r) +roesslerResults_5_3 = SubsegmentResultReduced.(load(dir_name*"roessler-results_4_1.jld2", "results")); +roessler_r = 0.5 +sig1_l,_ = subspaceFrequencyMatrix(roesslerResults_4_1, 1, roessler_r) diff --git a/examples/Tesis/roessler-results.jld2 b/examples/Tesis/roessler-results.jld2 new file mode 100644 index 0000000..a92ef5d --- /dev/null +++ b/examples/Tesis/roessler-results.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:103f127ec70d427ab3278207975aca1f5b6e9b5d960e351c6071362c1a093a59 +size 392053803 diff --git a/examples/Tesis/roessler-results_4_2.jld2 b/examples/Tesis/roessler-results_4_2.jld2 new file mode 100644 index 0000000..ab4315f --- /dev/null +++ b/examples/Tesis/roessler-results_4_2.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:0aa48b58f403af0064fae75d114c463fa568cf2a10a2618056dd8d7c16541e94 +size 96961020 diff --git a/examples/Tesis/roessler-results_5_3.jld2 b/examples/Tesis/roessler-results_5_3.jld2 new file mode 100644 index 0000000..d4f6f31 --- /dev/null +++ b/examples/Tesis/roessler-results_5_3.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:98ce4172f84bacb8657256fc7be395e81c4efe473baee9d6e6ca22aa57713056 +size 96261282 diff --git a/examples/Tesis/roessler-three-strips-experiment.jl b/examples/Tesis/roessler-three-strips-experiment.jl new file mode 100644 index 0000000..cef77b9 --- /dev/null +++ b/examples/Tesis/roessler-three-strips-experiment.jl @@ -0,0 +1,53 @@ +include("../../cycling-sigs/CyclingSignatures.jl/src/CyclingSignatures-include-file.jl") + +using DynamicalSystems +using JLD2 + +# define ODE +function roesslerThreeStrips(u0=[1.0, -2.0, 0.1]; a = 0.492, b = 2.0, c = 4.0) + return CoupledODEs(roessler_rule, u0, [a, b, c]) +end +@inbounds function roessler_rule(u, p, t) + du1 = -u[2]-u[3] + du2 = u[1] + p[1]*u[2] + du3 = p[2] + u[3]*(u[1] - p[3]) + return SVector{3}(du1, du2, du3) +end + +function roesslerThreeStripsTrajForCycling() + ds = roesslerThreeStrips() + dt = 0.1 + sol,_ = trajectory(ds, 10000, Δt=0.01) + + Y_init = Matrix(sol[:,:])'[:,100:end] + resample_boxsize = 1.0 + resample_sb_radius = 3.0 + roessler_f(v) = Vector(ds.integ.f(v, ds.integ.p, 0)) + + Y_res, t_vec = resampleToConsistent(ds, Y_init, resample_boxsize, dt, sb_radius = resample_sb_radius, sb_fct=roessler_f) + Z_res = mapslices(v->normalize(roessler_f(v),2 ), Y_res, dims=[1]) + return Y_res, Z_res, t_vec +end + +Y_roessler_three, Z_roessler_three, t_roessler_three = roesslerThreeStripsTrajForCycling() + +# run experiments with different r, k +# +r=4 +k=2 +runAndSaveExperimentWithRK("roesslerThree", Y_roessler_three, Z_roessler_three, r, k, t = t_roessler_three) +# +r=5 +k=3 +runAndSaveExperimentWithRK("roesslerThree", Y_roessler_three, Z_roessler_three, r, k, t = t_roessler_three) + +# load experiments +ts_roessler_three_4_2 = loadTrajectorySpace("roessler-three_4_2") +roesslerThreeResults_4_2 = SubsegmentResultReduced.(load(dir_name*"roesslerThree-results_4_2.jld2", "results")); +roesslerThree_r = 1 +sig1_l,_ = subspaceFrequencyMatrix(roesslerThreeResults_4_2, 1, roesslerThree_r) + +ts_roessler_three_5_3 = loadTrajectorySpace("roessler-three_5_3") +roesslerThreeResults_5_3 = SubsegmentResultReduced.(load(dir_name*"roesslerThree-results_5_3.jld2", "results")); +roesslerThree_r = 2 +sig1_l,_ = subspaceFrequencyMatrix(roesslerThreeResults_4_1, 1, roesslerThree_r) \ No newline at end of file diff --git a/examples/Tesis/roesslerThree-results_4_2.jld2 b/examples/Tesis/roesslerThree-results_4_2.jld2 new file mode 100644 index 0000000..ae15b44 --- /dev/null +++ b/examples/Tesis/roesslerThree-results_4_2.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:36ffa18218d0e0ed817d11731fc7ce6442d5569cc3d0c8d449842ddb62b0a1a0 +size 95526499 diff --git a/examples/Tesis/roesslerThree-results_5_3.jld2 b/examples/Tesis/roesslerThree-results_5_3.jld2 new file mode 100644 index 0000000..4a8581e --- /dev/null +++ b/examples/Tesis/roesslerThree-results_5_3.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:44420489aad9a08310022b20df01f67234630cbf7452e54a3a209fd2fa0fd5f3 +size 94090037 diff --git a/examples/Tesis/shilnikov-experiment.jl b/examples/Tesis/shilnikov-experiment.jl new file mode 100644 index 0000000..8923769 --- /dev/null +++ b/examples/Tesis/shilnikov-experiment.jl @@ -0,0 +1,72 @@ +using DynamicalSystems +using JLD2 + +# define ODE +function shilnikov(u0=[-0.23;-0.29;0.07]; a = 0.4, b = 0.58) + return CoupledODEs(shilnikov_rule, u0, [a, b]) +end +@inbounds function shilnikov_rule(u, p, t) + du = zeros(3, 1) + du1 = u[2] + du2 = u[3] + du3 = -p[1]*u[3] - u[2] + p[2]*u[1]*(1-u[1]^2) + return SVector{3}(du1, du2, du3) +end + +function shilnikovTrajForCycling() + ds = shilnikov() + dt = 0.1 + sol,_ = trajectory(ds, 10000, Δt=1e-2) + + Y_init = Matrix(sol[:,:])'[:,100:end] + resample_boxsize = 1.0 + resample_sb_radius = 3.0 + shilnikov_f(v) = Vector(ds.integ.f(v, ds.integ.p, 0)) + + Y_res, t_vec = resampleToConsistent(ds, Y_init, resample_boxsize, dt, sb_radius = resample_sb_radius, sb_fct=shilnikov_f) + Z_res = mapslices(v->normalize(shilnikov_f(v),2 ), Y_res, dims=[1]) + return Y_res, Z_res, t_vec +end + +Y_shilnikov, Z_shilnikov, t_shilnikov = shilnikovTrajForCycling() + +# r = 10, k = 9 +runAndSaveExperimentWithRK("shilnikov", Y_shilnikov, Z_shilnikov, 10, 9, t = t_shilnikov) +shilnikovResults_10_9 = SubsegmentResultReduced.(load(dir_name*"shilnikov-results_10_9.jld2", "results")); + +# r = 15, k = 9 +runAndSaveExperimentWithRK("shilnikov", Y_shilnikov, Z_shilnikov, 15, 9, t = t_shilnikov) +shilnikovResults_15_9 = SubsegmentResultReduced.(load(dir_name*"shilnikov-results_15_9.jld2", "results")); + +# r = 15, k = 10 +runAndSaveExperimentWithRK("shilnikov", Y_shilnikov, Z_shilnikov, 15, 10, t = t_shilnikov) +shilnikovResults_15_10 = SubsegmentResultReduced.(load(dir_name*"shilnikov-results_15_10.jld2", "results")); + + +### +boxsize_shilnikov = 15 +sb_radius_shilnikov = 9 +# comparison space con betti = 1 +ts_shilnikov_1 = trajectoryToTrajectorySpaceSB(Y_shilnikov, Z_shilnikov, boxsize_shilnikov, sb_radius_shilnikov, t_vec = t_shilnikov; filter_missing=true, shortest_path_fix=true) + +boxsize_shilnikov = 15 +sb_radius_shilnikov = 10 +# comparison space con betti = 7 +ts_shilnikov_2 = trajectoryToTrajectorySpaceSB(Y_shilnikov, Z_shilnikov, boxsize_shilnikov, sb_radius_shilnikov, t_vec = t_shilnikov; filter_missing=true, shortest_path_fix=true) + +boxsize_shilnikov = 8 +sb_radius_shilnikov = 3 +ts_shilnikov_3 = trajectoryToTrajectorySpaceSB(Y_shilnikov, Z_shilnikov, boxsize_shilnikov, sb_radius_shilnikov, t_vec = t_shilnikov; filter_missing=true, shortest_path_fix=true) + +n = 1000 +shilnikovExperimentParams = map(i->SubsegmentSampleParameter(i,n), 10:10:1000) + +shilnikovExperiments = map(shilnikovExperimentParams) do p + RandomSubsegmentExperiment(ts_shilnikov, p, boxsize_shilnikov) +end + +times, shilnikovResults = runExperimentsWithTimer(shilnikovExperiments) + +dir_name = "examples/CustomExperiments/" +jldsave(dir_name*"shilnikov-results.jld2", shilnikovResults=shilnikovResults, times=times) +shilnikovResults = SubsegmentResultReduced.(load(dir_name*"shilnikov-results.jld2", "shilnikovResults")); \ No newline at end of file diff --git a/examples/Tesis/shilnikov-results.jld2 b/examples/Tesis/shilnikov-results.jld2 new file mode 100644 index 0000000..fb60754 --- /dev/null +++ b/examples/Tesis/shilnikov-results.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:12a25d2eef8531fc11838cc485d50189549aa9271538b8e75bae383a99dc5fa3 +size 85282539 diff --git a/examples/Tesis/shilnikov-results_10_9.jld2 b/examples/Tesis/shilnikov-results_10_9.jld2 new file mode 100644 index 0000000..5a449d1 --- /dev/null +++ b/examples/Tesis/shilnikov-results_10_9.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1998155f1888faf2a9b238ff259ac8c5c8d95d36bb2495234a50405f44d6576c +size 96127568 diff --git a/examples/Tesis/shilnikov-results_15_10.jld2 b/examples/Tesis/shilnikov-results_15_10.jld2 new file mode 100644 index 0000000..0d2795f --- /dev/null +++ b/examples/Tesis/shilnikov-results_15_10.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d5fbc8004362896a85295af09475df10d8b5f484a1e650ee6e6ea6be8b6df60a +size 96988650 diff --git a/examples/Tesis/shilnikov-results_15_9.jld2 b/examples/Tesis/shilnikov-results_15_9.jld2 new file mode 100644 index 0000000..4981a0b --- /dev/null +++ b/examples/Tesis/shilnikov-results_15_9.jld2 @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c9230c237dcce18bdbb8b7abf07e0c1e1c215fa229b89f13e35f5e6139ff4c3f +size 96149252 diff --git a/examples/Tesis/tesis.jl b/examples/Tesis/tesis.jl new file mode 100644 index 0000000..a1bd2e7 --- /dev/null +++ b/examples/Tesis/tesis.jl @@ -0,0 +1,59 @@ +# load experiments and plot cycling sigs + +################################################################################################ +# double well +# .2, 3 +dwResults = SubsegmentResultReduced.(load("examples/PaperExperiments/dw-results.jld2", "dwResults")); +dw_r = 0.18 +plotCyclingSignatures(dwResults, dw_r) + +################################################################################################ +# lorenz +# 8, 3 +lorenzResults = SubsegmentResultReduced.(load("examples/PaperExperiments/lorenz-results.jld2", "lorenzResults")); +lorenz_r = 6 +plotCyclingSignatures(lorenzResults, lorenz_r) + +################################################################################################ +# burkeshaw +# 5, 3 +burkeshawResults_5_3 = SubsegmentResultReduced.(load("examples/Tesis/burkeshaw-results_5_3.jld2", "results")); +plotCyclingSignatures(burkeshawResults_5_3, 4) +plotSignaturesWithLegend(burkeshawResults_5_3, 1, 4; cutoff=20) +# 8, 3 +burkeshawResults_8_3 = SubsegmentResultReduced.(load("examples/Tesis/burkeshaw-results_8_3.jld2", "results")); +plotCyclingSignatures(burkeshawResults_8_3, 5) +plotSignaturesWithLegend(burkeshawResults_8_3, 1, 5; cutoff=20) + +################################################################################################ +# roessler +# 4, 2 +roesslerResults_4_2 = SubsegmentResultReduced.(load("examples/Tesis/roessler-results_4_2.jld2", "roesslerResults")); +plotCyclingSignatures(roesslerResults_4_2, 3) +plotSignaturesWithLegend(roesslerResults_4_2, 1, 3; cutoff=20) +# 5, 3 +roesslerResults_5_3 = SubsegmentResultReduced.(load("examples/Tesis/roessler-results_5_3.jld2", "results")); +plotCyclingSignatures(roesslerResults_5_3, 4) +plotSignaturesWithLegend(roesslerResults_5_3, 1, 4; cutoff=20) + +################################################################################################ +# roessler three strips +# 4, 2 +roesslerThreeResults_4_2 = SubsegmentResultReduced.(load("examples/Tesis/roesslerThree-results_4_2.jld2", "results")); +plotCyclingSignatures(roesslerThreeResults_4_2, 3) +plotSignaturesWithLegend(roesslerThreeResults_4_2, 1, 3; cutoff=20) +# 5, 3 +roesslerThreeResults_5_3 = SubsegmentResultReduced.(load("examples/Tesis/roesslerThree-results_5_3.jld2", "results")); +plotCyclingSignatures(roesslerThreeResults_5_3, 4) +plotSignaturesWithLegend(roesslerThreeResults_5_3, 1, 4; cutoff=20) + +################################################################################################ +# shilkinov +# 15, 9 +shilnikovResults_15_9 = SubsegmentResultReduced.(load("examples/Tesis/shilnikov-results_15_9.jld2", "results")); +plotCyclingSignatures(shilnikovResults_15_9, 12) +plotSignaturesWithLegend(shilnikovResults_15_9, 1, 12; cutoff=20) +# 15, 10 +shilnikovResults_15_10 = SubsegmentResultReduced.(load("examples/Tesis/shilnikov-results_15_10.jld2", "results")); +plotCyclingSignatures(shilnikovResults_15_10, 10) +plotSignaturesWithLegend(shilnikovResults_15_10, 1, 12; cutoff=20) \ No newline at end of file diff --git a/src/CyclingSignatures/sphere-bundle.jl b/src/CyclingSignatures/sphere-bundle.jl index fcdfdbb..cf57321 100644 --- a/src/CyclingSignatures/sphere-bundle.jl +++ b/src/CyclingSignatures/sphere-bundle.jl @@ -31,6 +31,15 @@ function sphereBundleDistanceMatrix(SB_pts) return dm end +function simplifiedSphereBundleDistanceMatrix(SB_pts) + d = div(size(SB_pts,1),2) + SB_X = SB_pts[1:d,:] + SB_VF= SB_pts[d+1:end,:] + println(size(SB_X)) + @time dm = pairwise(chebyshev, SB_X) + return dm +end + struct SBDistance <: Metric dim::Int c::Float64 diff --git a/src/CyclingSignatures/trajectory-space.jl b/src/CyclingSignatures/trajectory-space.jl index 62a6b5e..bcb98e8 100644 --- a/src/CyclingSignatures/trajectory-space.jl +++ b/src/CyclingSignatures/trajectory-space.jl @@ -86,6 +86,18 @@ function Base.show(io::IO, bs::BoxSpace) print(io, "with β1 = $(betti_1(bs.h1))") end +# function SimplifiedBoxSpace(SB_pts, boxsize, sphereBundleRadius; kwargs...) +# SB_pts = unique(SB_pts, dims=2) +# dm_pts = simplifiedSphereBundleDistanceMatrix(SB_pts) +# cplx, h1_gen, h1 = boxSpaceH1Info(dm_pts) + +# X,VF,SB_ind = sphereBundleToComponents(SB_pts) + +# inc_helper = SBInclusionHelper(cplx, SB_pts, h1_gen, boxsize, sphereBundleRadius; kwargs...) + +# return SBBoxSpace(X,VF,SB_ind,boxsize, sphereBundleRadius,h1, inc_helper) +# end + struct SBBoxSpace{S<:AbstractVector{Int},T<:Real} <: AbstractBoxSpace X::Matrix{Int} # quantization in space VF::Matrix{Int} # quantization in bundle @@ -114,6 +126,10 @@ function Base.show(io::IO, bs::SBBoxSpace) print(io, "with β1 = $(betti_1(bs.h1))") end +function getPlotSimplifiedPoints(bs::BoxSpace) + return bs.X +end + function getPlotPoints(bs::SBBoxSpace,offset=.1) X = bs.X VF = bs.VF @@ -170,6 +186,25 @@ function trajectoryToTrajectorySpace(Y, boxsize; t_vec=nothing, metric=euclidean return TrajectorySpace(rt, boxSpace, metric) end +function trajectoryToTrajectorySpaceSimplifiedSB(Y,Z, boxsize, sb_radius; t_vec=nothing, metric=euclidean, kwargs...) + for c in eachcol(Z) + if !isapprox(norm(c),1) + error("Sphere Bundle component is not an unit vectors.") + end + end + rt = nothing + ts = [Y;Z] + if isnothing(t_vec) + rt = ResampledTrajectory(ts) + else + rt = ResampledTrajectory(ts, t_vec) + end + + SB_pts = unique([quantize(Y, boxsize);quantizeSB(Z,sb_radius)],dims=2) + boxSpace = SimplifiedBoxSpace(SB_pts, boxsize, sb_radius; kwargs...) + return TrajectorySpace(rt, boxSpace, metric) +end + function trajectoryToTrajectorySpaceSB(Y,Z, boxsize, sb_radius; t_vec=nothing, metric=nothing, kwargs...) for c in eachcol(Z) if !isapprox(norm(c),1)