Working version, but unnecessarily overloads SciMLBase.@recipe

modified:   Project.toml
modified:   src/MechGlueDiffEqBase.jl
modified:   test/Manifest.toml
modified:   test/Project.toml
modified:   test/runtests.jl

Author: hustf

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.31"
 
 [deps]
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
+MechGluePlots = "03a24088-0a82-4850-846f-2cad5d056ae1"
 MechanicalUnits = "e6be9192-89dc-11e9-36e6-5dbcb28f419e"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"

src/MechGlueDiffEqBase.jl

@@ -1,11 +1,18 @@
 module MechGlueDiffEqBase
-import Unitfu: AbstractQuantity, Quantity
+import Unitfu: AbstractQuantity, Quantity, ustrip
 import DiffEqBase: value, ODE_DEFAULT_NORM, UNITLESS_ABS2
+export value, ODE_DEFAULT_NORM, UNITLESS_ABS2, Unitfu, AbstractQuantity, Quantity
+
+
+import RecipesBase
+import RecipesBase.@recipe
 import SciMLBase
 import SciMLBase: AbstractTimeseriesSolution, RecipesBase
 import SciMLBase: AbstractDiscreteProblem, AbstractRODESolution, SensitivityInterpolation
-import RecipesBase.@recipe
-export value, ODE_DEFAULT_NORM, UNITLESS_ABS2, Unitfu
+import SciMLBase: getsyms, interpret_vars, cleansyms
+import SciMLBase: diffeq_to_arrays, issymbollike, getindepsym_defaultt
+import SciMLBase: DEFAULT_PLOT_FUNC
+
 
 function value(x::Type{AbstractQuantity{T,D,U}}) where {T,D,U}
     T
@@ -26,7 +33,8 @@ end
     real(abs2(x)/oneunit(x)*oneunit(x))
 end
 
-
+# This recipe is copied from SciMLBase, and only the signature is adopted
+# for use with unitful solutions. It can perhaps be deleted now, and the link to RecipesBase be deleted.
 @recipe function f(sol::AbstractTimeseriesSolution{T, N, A};
     plot_analytic=false,
     denseplot = (sol.dense ||
@@ -41,13 +49,91 @@ end
                  1000*sol.tslocation),
     tspan = nothing, axis_safety = 0.1,
     vars=nothing) where {T<:Quantity, N, A<:Array{<:Quantity}}
-
-    print("----------------------------------MechGlueDiffEqBase")
-    return Unifu.ustrip(sol.u), Unifu.ustrip(sol.t)
+    syms = getsyms(sol)
+    int_vars = interpret_vars(vars,sol,syms)
+    strs = cleansyms(syms)
+  
+    tscale = get(plotattributes, :xscale, :identity)
+    plot_vecs,labels = diffeq_to_arrays(sol,plot_analytic,denseplot,
+                                        plotdensity,tspan,axis_safety,
+                                        vars,int_vars,tscale,strs)
+  
+    tdir = sign(sol.t[end]-sol.t[1])
+    xflip --> tdir < 0
+    seriestype --> :path
+  
+    # Special case labels when vars = (:x,:y,:z) or (:x) or [:x,:y] ...
+    if typeof(vars) <: Tuple && (issymbollike(vars[1]) && issymbollike(vars[2]))
+      xguide --> issymbollike(int_vars[1][2]) ? Symbol(int_vars[1][2]) : strs[int_vars[1][2]]
+      yguide --> issymbollike(int_vars[1][3]) ? Symbol(int_vars[1][3]) : strs[int_vars[1][3]]
+      if length(vars) > 2
+        zguide --> issymbollike(int_vars[1][4]) ? Symbol(int_vars[1][4]) : strs[int_vars[1][4]]
+      end
+    end
+  
+    if (!any(issymbollike,getindex.(int_vars,1)) && getindex.(int_vars,1) == zeros(length(int_vars))) ||
+       (!any(issymbollike,getindex.(int_vars,2)) && getindex.(int_vars,2) == zeros(length(int_vars))) ||
+       all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,1)) || all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,2))
+      xguide --> "$(getindepsym_defaultt(sol))"
+    end
+    if length(int_vars[1]) >= 3 && ((!any(issymbollike,getindex.(int_vars,3)) && getindex.(int_vars,3) == zeros(length(int_vars))) ||
+       all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,3)))
+      yguide --> "$(getindepsym_defaultt(sol))"
+    end
+    if length(int_vars[1]) >= 4 && ((!any(issymbollike,getindex.(int_vars,4)) && getindex.(int_vars,4) == zeros(length(int_vars))) ||
+       all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,4)))
+      zguide --> "$(getindepsym_defaultt(sol))"
+    end
+  
+    if (!any(issymbollike,getindex.(int_vars,2)) && getindex.(int_vars,2) == zeros(length(int_vars))) ||
+        all(t->Symbol(t)==getindepsym_defaultt(sol),getindex.(int_vars,2))
+      if tspan === nothing
+        if tdir > 0
+          xlims --> (sol.t[1],sol.t[end])
+        else
+          xlims --> (sol.t[end],sol.t[1])
+        end
+      else
+        xlims --> (tspan[1],tspan[end])
+      end
+    else
+      mins = minimum(sol[int_vars[1][2],:])
+      maxs = maximum(sol[int_vars[1][2],:])
+      for iv in int_vars
+        mins = min(mins,minimum(sol[iv[2],:]))
+        maxs = max(maxs,maximum(sol[iv[2],:]))
+      end
+      xlims --> ((1-sign(mins)*axis_safety)*mins,(1+sign(maxs)*axis_safety)*maxs)
+    end
+  
+    # Analytical solutions do not save enough information to have a good idea
+    # of the axis ahead of time
+    # Only set axis for animations
+    if sol.tslocation != 0 && !(typeof(sol) <: AbstractAnalyticalSolution)
+      if all(getindex.(int_vars,1) .== DEFAULT_PLOT_FUNC)
+        mins = minimum(sol[int_vars[1][3],:])
+        maxs = maximum(sol[int_vars[1][3],:])
+        for iv in int_vars
+          mins = min(mins,minimum(sol[iv[3],:]))
+          maxs = max(maxs,maximum(sol[iv[3],:]))
+        end
+        ylims --> ((1-sign(mins)*axis_safety)*mins,(1+sign(maxs)*axis_safety)*maxs)
+  
+        if length(int_vars[1]) >= 4
+          mins = minimum(sol[int_vars[1][4],:])
+          maxs = maximum(sol[int_vars[1][4],:])
+          for iv in int_vars
+            mins = min(mins,minimum(sol[iv[4],:]))
+            maxs = max(mins,maximum(sol[iv[4],:]))
+          end
+          zlims --> ((1-sign(mins)*axis_safety)*mins,(1+sign(maxs)*axis_safety)*maxs)
+        end
+      end
+    end
+  
+    label --> reshape(labels,1,length(labels))
+    println("----------------------------------MechGlueDiffEqBase")
+    (plot_vecs...,)
 end
-        # T = Quantity{Float64,  ᴸ∙ ᴹ∙ ᵀ⁻², Unitfu.FreeUnits{(N,),  ᴸ∙ ᴹ∙ ᵀ⁻², nothing}}
-    # N = 1
-    # A = Vector{Quantity{Float64,  ᴸ∙ ᴹ∙ ᵀ⁻², Unitfu.FreeUnits{(N,),  ᴸ∙ ᴹ∙ ᵀ⁻², nothing}}}
-
 
 end

test/Manifest.toml

@@ -613,11 +613,17 @@ uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
 version = "0.3.1"
 
 [[MechGlueDiffEqBase]]
-deps = ["DiffEqBase", "MechanicalUnits", "RecipesBase", "SciMLBase", "Unitfu"]
+deps = ["DiffEqBase", "MechGluePlots", "MechanicalUnits", "RecipesBase", "SciMLBase", "Unitfu"]
 path = "C:\\Users\\frohu_h4g8g6y\\.julia\\dev\\MechGlueDiffEqBase"
 uuid = "2532746b-52b5-4539-9431-8bb183ab067f"
 version = "0.1.31"
 
+[[MechGluePlots]]
+deps = ["Plots", "RecipesBase", "Unitfu"]
+git-tree-sha1 = "2b4328a7d4eb2ee41cc060df6cddd9f1a0b1de23"
+uuid = "03a24088-0a82-4850-846f-2cad5d056ae1"
+version = "0.2.1"
+
 [[MechanicalUnits]]
 deps = ["Base64", "Distributed", "InteractiveUtils", "Libdl", "LinearAlgebra", "Logging", "Markdown", "Random", "Serialization", "Sockets", "Test", "Unitfu"]
 git-tree-sha1 = "0e0eca9fbd0d9c575d8cd106979c350bf0fe057b"

test/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 MechGlueDiffEqBase = "2532746b-52b5-4539-9431-8bb183ab067f"
+MechGluePlots = "03a24088-0a82-4850-846f-2cad5d056ae1"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"

test/runtests.jl

@@ -1,5 +1,6 @@
 using MechGlueDiffEqBase
-import MechanicalUnits: @import_expand, Quantity, Time, Force, ∙, Quantity
+using MechGluePlots
+import MechanicalUnits: @import_expand, Quantity, Time, Force, ∙
 import Unitfu.AbstractQuantity
 import DifferentialEquations
 import DifferentialEquations: SciMLBase, solve, ODEProblem, Tsit5
@@ -8,7 +9,8 @@ import SciMLBase: AbstractTimeseriesSolution, RecipesBase
 import SciMLBase: AbstractDiscreteProblem, AbstractRODESolution, SensitivityInterpolation
 using Test
 import Plots.plot
-@import_expand(N, s)
+
+@import_expand(N, s, m)
 @testset "ODE" begin
     f = (y, p, t) -> 0.5y / 3.0s  # The derivative, yeah
     u0 = 1.5N
@@ -20,7 +22,7 @@ import Plots.plot
     @test sol.u isa Vector{<:Force}
 end
 
-@testset "Plot ODE quantity" begin
+#@testset "Plot ODE quantity" begin
     SciMLBase.RecipesBase.debug(true)
     f1 = (y, p, t) -> -1.5y / 0.3s  # The derivative, yeah
     u0 = 1.5N
@@ -29,17 +31,11 @@ end
     solve(prob, Tsit5())
     sol = solve(prob)
     @test sol isa AbstractTimeseriesSolution{T, N, A} where {T, N, A}
-    @test sol isa AbstractTimeseriesSolution{T, N, A} where {T<:AbstractQuantity, N, A}
-    @test sol isa AbstractTimeseriesSolution{T, N, A} where {T<:AbstractQuantity, N, A<:AbstractArray{<:AbstractQuantity}}
+    @test sol isa AbstractTimeseriesSolution{T, N, A} where {T<:Quantity, N, A}
     @test sol isa AbstractTimeseriesSolution{T, N, A} where {T<:Quantity, N, A<:Array{<:Quantity}}
+    # Here, we only test this part of the dispatch mechanism.
+    # This is intended to work with MechGluePlots also loaded.
+    # See other tests in MechGluecode.
     plot(sol)
-    # TODO implement in MechGlueDiffEqBase
-    # the recipe found in SciMLBase/solutions/solution_interface.jl:151
-    # Type to dispatch on:
-    # AbstractTimeseriesSolution{T, N, A}) where {T, N, A}
-    # T = Quantity{Float64,  ᴸ∙ ᴹ∙ ᵀ⁻², Unitfu.FreeUnits{(N,),  ᴸ∙ ᴹ∙ ᵀ⁻², nothing}}
-    # N = 1
-    # A = Vector{Quantity{Float64,  ᴸ∙ ᴹ∙ ᵀ⁻², Unitfu.FreeUnits{(N,),  ᴸ∙ ᴹ∙ ᵀ⁻², nothing}}}
-
-end
+#end