text:
#=
# Tests multiple interdependent scalar variables.
=#

### If the Finch package has already been added, use this line #########
using Finch # Note: to add the package, first do: ]add "https://github.com/paralab/Finch.git"

### If not, use these four lines (working from the examples directory) ###
# if !@isdefined(Finch)
#     include("../Finch.jl");
#     using .Finch
# end
##########################################################################

init_finch("multivar");

# Try making an optional log
useLog("multivarlog", level=3)

# Set up the configuration
@domain(1)
@functionSpace(LEGENDRE, 3)

#@matrixFree(200,1e-6)

# Build a simple mesh
@mesh(LINEMESH, 30)

# Variables
@variable(u)
@variable(q)

@testSymbol(v)

#Boundary conditions
@boundary(u, 1, DIRICHLET, "0")
@boundary(q, 1, DIRICHLET, "0")

#Equations
@coefficient(f, "-4*pi*pi*sin(2*pi*x)")
@coefficient(g, "-9*pi*pi*sin(3*pi*x) + 0.5*sin(2*pi*x)")

@weakForm([u, q], ["-dot(grad(u), grad(v)) - f*v", "-dot(grad(q), grad(v)) + 0.5*u*v - g*v"])

solve([u, q]);

# exact solution is u = sin(2*pi*x), q = sin(3*pi*x)
# check error
maxerru = 0;
maxerrq = 0;
exactu(x) = sin(2*pi*x);
exactq(x) = sin(3*pi*x);

for i=1:size(Finch.grid_data.allnodes,2)
    x = Finch.grid_data.allnodes[1,i];
    erru = abs(u.values[i] - exactu(x));
    errq = abs(q.values[i] - exactq(x));
    global maxerru;
    global maxerrq;
    maxerru = max(erru,maxerru);
    maxerrq = max(errq,maxerrq);
end
println("max error(u) = "*string(maxerru));
println("max error(q) = "*string(maxerrq));

# using Plots
# pyplot();
# display(plot(Finch.grid_data.allnodes[:], u.values[:], markershape=:circle, reuse=false))
# display(plot(Finch.grid_data.allnodes[:], q.values[:], markershape=:circle, reuse=false))

# check
log_dump_config();
log_dump_prob();

@finalize()

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using SymbolicUtils, SpecialFunctions
using SymbolicUtils.Rewriters
using Symbolics

SEMIINTRULES = [
        # CONSTANTS AND POWERS
        @acrule sqrt(~x) => sqrt(pi)/2*~x
        #@acrule((~x) => 4*~x^(3/2)/(3*sqrt(pi)))
        @acrule (~x)^(~α) => gamma(~α+1)/gamma(~α+3/2)*(~x)^(~α+1/2)

        # BINOMIALS
        @acrule sqrt(1+~x) => sqrt(~x/pi) + ((1+~x)*atan(sqrt(~x)))/sqrt(pi)
        @acrule sqrt(1-~x) => sqrt(~x/pi) + ((1-~x)*atanh(sqrt(~x)))/sqrt(pi)
        @acrule 1/sqrt(1+~x) => 2/sqrt(pi)*atan(sqrt(~x))
        @acrule 1/sqrt(1-~x) => 2/sqrt(pi)*atanh(sqrt(~x))
        @acrule 1/(1+~x) => 2*asinh(sqrt(~x))/sqrt(pi*(1+~x))
        @acrule 1/(1-~x) => 2*asin(sqrt(~x))/sqrt(pi*(1-~x))
        @acrule 1/(1+~x)^(3/2) => 2/(1+~x)*sqrt(~x/pi)
        @acrule 1/(1-~x)^(3/2) => 2/(1-~x)*sqrt(~x/pi)
        @acrule 1/(1-~x)^(~p) => -incomplete_beta(~x, 1/2, ~p-1/2)/(sqrt(pi)*(1-~x)^(~p-1/2))
        @acrule 1/sqrt(~x*(1-~x)) => 2*ellipk(~x)/sqrt(pi)
        @acrule 1/sqrt(~x*(1+~x)) => 2*ellipk(~x/(1+~x))/sqrt(pi*(1+~x))
        @acrule sqrt(~x/(1-~x)) => 2/sqrt(pi)*(ellipk(~x)-ellipe(~x))
        @acrule 1/(sqrt(~x)*(1+~x)) => sqrt(pi/(1+~x))
        @acrule 1/(sqrt(~x)*(1-~x)) => sqrt(pi/(1-~x))
        @acrule sqrt(~x)/(~x) => -sqrt(pi/(1+~x))+sqrt(pi)
        @acrule sqrt(~x)/(-~x) => sqrt(pi/(1-~x))-sqrt(pi)
        @acrule sqrt(~x)/(1-~x)^2 => sqrt(pi)*~x/(2*(1-~x)^(3/2))
        @acrule sqrt((1-~x)/~x) => 2/sqrt(pi)*ellipe(~x)
        @acrule sqrt((1+~x)/~x) => 2*sqrt((1+~x)/pi)*ellipe(~x/(1+~x))
        @acrule  (~x)^(~p)/(1+~x)^~r => gamma(~p+1)/gamma(~p+3/2)*(~x)^(~p+1/2)* Hypergeometric1F1(~r, ~p+1, p+3/2, -~x)
        @acrule  (~x)^(~p)/(1-~x)^~r => gamma(~p+1)/gamma(~p+3/2)*(~x)^(~p+1/2)* Hypergeometric1F1(~r, ~p+1, p+3/2, ~x)

        # EXPONENTIAL AND RELATED FUNCTIONS
        @acrule exp(~x) => exp(~x)*erf(sqrt(~x))
        @acrule exp(-~x) => 2/sqrt(pi)*dawson(sqrt(~x))
        @acrule exp(~x)*erf(sqrt(~x)) => exp(~x)-1
        @acrule dawson(~x) => 1/2*sqrt(pi)*(1-exp(-~x))
        @acrule exp(~x)*erfc(sqrt(~x)) => 1-exp(~x)*erfc(sqrt(~x))
        @acrule exp(~x)*erfc(-sqrt(~x)) => exp(~x)*erfc(sqrt(~x))-1
        @acrule erf(sqrt(~x)) => ~x*exp(-1/2*~x)*(besseli(1, 1/2*~x)+besseli(0, 1/2*~x))
        @acrule exp(~x)/sqrt(~x) => sqrt(pi)*exp(1/2*~x)*besseli(0, 1/2*~x)
        @acrule exp(-~x)/sqrt(~x) => sqrt(pi)*exp(-1/2*~x)*besseli(0, 1/2*~x)
        @acrule cosh(sqrt(~x))/sqrt(~x) => sqrt(pi)*besseli(0, sqrt(~x))
        @acrule (1-cos(sqrt(~x)))/~x => sqrt(pi*(~x))*Struve(1, sqrt(~x))
        @acrule (1-cosh(sqrt(~x)))/~x => sqrt(pi*(~x))*MStruve(1, sqrt(~x))
        @acrule sin(~x) => sin(~x-1/4*pi)+sqrt(2)*AuxiliaryFresnelSin(sqrt(2*~x/pi))
        @acrule cos(~x) => cos(~x-1/4*pi)-sqrt(2)*AuxiliaryFresnelCos(sqrt(2*~x/pi))

        @acrule sinh(~x) => (exp(~x)*erf(sqrt(~x)))/2-dawson(sqrt(~x))/sqrt(pi)
        @acrule cosh(~x) => (exp(~x)*erf(sqrt(~x)))/2+dawson(sqrt(~x))/sqrt(pi)
        @acrule asin(sqrt(~x))/sqrt(1-~x) => -sqrt(pi)/2*log(1-~x)
        @acrule atan(sqrt(~x)) => sqrt(pi)*(sqrt(1+~x)-1)
        @acrule asinh(sqrt(~x))/sqrt(1+~x) => sqrt(pi)/2*log(1+~x)
        @acrule atanh(sqrt(~x)) => sqrt(pi)*(1-sqrt(1-~x))

        # BESSEL AND STRUVE FUNCTIONS
        @acrule besselj(0, sqrt(~x)) => sin(sqrt(~x))/sqrt(pi*(~x))
        @acrule besselj(1, sqrt(~x)) => 2*(1-cos(sqrt(~x)))/sqrt(pi*~x)
        @acrule besselj(~v, sqrt(~x))/(~x)^(~v/2) => (sqrt(2)*Struve(~v-1/2, sqrt(~x)))/(~x)^(~v/2-1/4)

        @acrule sqrt(~x)*besselj(1, sqrt(~x)) => (2*sin(sqrt(~x))-2*sqrt(~x)*cos(sqrt(~x)))/sqrt(pi)
        @acrule (~x)^(~v/2)*besselj(~v, sqrt(~x)) => (~x)^((~v)/2+1/4)*besselj(~v+1/2, sqrt(~x))
        @acrule besseli(0, sqrt(~x)) => sinh(sqrt(~x))/sqrt(pi*~x)
        @acrule besseli(1, sqrt(~x))/sqrt(~x) => 2*(1-cosh(sqrt(~x)))/sqrt(pi*~x)
        @acrule besseli(~v, sqrt(~x))/(~x)^(~v/2) => (sqrt(2)*Legendre(~v-1/2, sqrt(~x)))/(~x)^(~v/2-1/4)
        @acrule sqrt(~x)*besseli(1, sqrt(~x)) => (2*sqrt(~x)*cosh(sqrt(~x))-2*sinh(sqrt(~x)))/sqrt(pi)
        @acrule (~x)^(~v/2)*besseli(~v, sqrt(~x)) => sqrt(2)*(~x)^((~v)/2+1/4)*besseli(~v+1/2, sqrt(~x))
        @acrule exp(-~x)*besseli(0, ~x) => erf(sqrt(2*~x))/sqrt(2)
        @acrule Struve(0, sqrt(~x)) => 2*(1-cos(sqrt(~x)))/sqrt(pi)

        @acrule  Struve(0, sqrt(~x))/(~x) => 2*SinIntegral(sqrt(~x))/(sqrt(pi*(~x)))
        @acrule sqrt(~x)*Struve(1, sqrt(~x)) =>(~x+2-2*sqrt(~x)*sin(sqrt(~x))-2*cos(sqrt(~x)))/sqrt(pi)
        @acrule (~x)^(~v/2)*Struve(~v, sqrt(~x)) => sqrt(2)*(~x)^((~v)/2+1/4)*Struve(~v+1/2, sqrt(~x))

        @acrule MStruve(0, sqrt(~x)) => (2*(cosh(sqrt(~x))-1))/sqrt(pi)
        @acrule MStruve(0, sqrt(~x))/(~x) => (HyperSinIntegral(sqrt(~x)))/(sqrt(pi)*(~x))
        @acrule sqrt(~x)*MStruve(1, sqrt(~x)) => (2-~x+2*sqrt(~x)*sinh(sqrt(~x))-2*cosh(sqrt(~x)))/sqrt(pi)
        @acrule  (~x)^(~v/2)*MStruve(~v, sqrt(~x)) => sqrt(2)*(~x)^((~v)/2+1/4)*MStruve(~v+1/2, sqrt(~x))

        # MISCELLANEOUS FUNCTIONS
        @acrule log(~x) => 2*sqrt(~x/pi)*(log(4*~x-2))
        @acrule sqrt(~x)*log(~x) => (~x*sqrt(pi))/2*(log(~x/4)+1)
        @acrule log(~x)/sqrt(~x) => sqrt(pi)*log(~x/4)

        @acrule ellipk(~x) => sqrt(pi)*asin(sqrt(~x))
        @acrule ellipe(~x) => sqrt(pi*(~x)*(1-~x))/2 + sqrt(pi)*asin(sqrt(~x))/2
]

INTRULES = Chain(SEMIINTRULES)

"""
# Symbolic fractional integral

    semiint(fun)

```semiint``` uses SymbolicUtils.jl and Symbolics.jl to compute symbolic fractional integral.

### Example

```julia-repl
julia> using SymbolicUtils
julia> @syms x
julia> semiint(x^4)
0.45851597901024005(x^4.5)
```
"""
semiint(x) = INTRULES(x)

ext: jl


text:
"""
D-dimensional cubic lattice.
"""
struct CubicLattice{T<:AbstractArray{Int}} <: AbstractLattice
    L::Int
    dim::Int
    sites::Int
    # row = first half uprights, second half downlefts, D in total; col = siteidx
    neighs::Matrix{Int}
    lattice::T

    # generic checkerboard
    n_bonds::Int
    bonds::Matrix{Int}
end

# constructors
"""
    CubicLattice(D, L)

Create a D-dimensional cubic lattice with linear dimension `L`.
"""
function CubicLattice(D::Int, L::Int)
    sites = L^D
    lattice = convert(Array, reshape(1:sites, (fill(L, D)...,)))
    n_bonds = D * sites
    neighs = build_neighbortable(CubicLattice, lattice, D)

    bonds = Matrix{Int}(undef, n_bonds, 3)
    bondid = 1
    for src in lattice
        for trg in neighs[1:D, src]
            bonds[bondid, 1] = src
            bonds[bondid, 2] = trg
            bonds[bondid, 3] = 0
            bondid += 1
        end
    end

    return CubicLattice(L,D,sites,Matrix(neighs),lattice,n_bonds,bonds)
end

function build_neighbortable(::Type{CubicLattice}, lattice, D)
    uprights = Vector{Vector{Int}}(undef, D)
    downlefts = Vector{Vector{Int}}(undef, D)

    for d in 1:D
        shift = zeros(Int, D); shift[d]=-1;
        uprights[d] = circshift(lattice, (shift...,))[:]
        shift[d]=1;
        downlefts[d] = circshift(lattice, (shift...,))[:]
    end

    return transpose(hcat(uprights..., downlefts...))
end

# Implement AbstractLattice interface: mandatory
@inline Base.length(c::CubicLattice) = c.sites
@inline Base.size(c::CubicLattice) = tuple((c.L for _ in 1:l.dims)...)

# Implement AbstractLattice interface: optional
@inline neighbors_lookup_table(c::CubicLattice) = copy(c.neighs)

# HasNeighborsTable and HasBondsTable traits
has_neighbors_table(::CubicLattice) = HasNeighborsTable()
has_bonds_table(::CubicLattice) = HasBondsTable()

positions(l::CubicLattice) = l.lattice |> CartesianIndices .|> Tuple .|> collect
lattice_vectors(l::CubicLattice) = [[l.L, 0, 0], [0, l.L, 0], [0, 0, l.L]]


text:
function cache_tests()

    @info "cache_tests"
    ProjAssistant.global_conf(:CACHE_DIR, tempdir())
    N = 15
    dat0 = rand(N)

    cfile = ProjAssistant.cfname("my_dat", N)
    cfile = ProjAssistant.scache(dat0, cfile)
    @test isfile(cfile)
    dat1 = ProjAssistant.lcache(cfile)
    @test all(dat0 .== dat1)
    ProjAssistant.delcache("my_dat", N)
    @test !isfile(cfile)

    cfile = ProjAssistant.scache(dat0, "my_dat", N)
    dat1 = ProjAssistant.lcache("my_dat", N)
    @test all(dat0 .== dat1)
    ProjAssistant.delcache("my_dat", N)
    @test !isfile(cfile)

    cfile = ProjAssistant.scache(() -> dat0, "my_dat", N)
    dat1 = ProjAssistant.lcache("my_dat", N)
    @test all(dat0 .== dat1)
    ProjAssistant.delcache(cfile)
    @test !isfile(cfile)

    cfile = ProjAssistant.scache(() -> dat0)
    dat1 = ProjAssistant.lcache(cfile)
    @test all(dat0 .== dat1)
    ProjAssistant.delcache(cfile)
    @test !isfile(cfile)

    cfile = ProjAssistant.scache(dat0)
    dat1 = ProjAssistant.lcache(cfile)
    @test all(dat0 .== dat1)
    ProjAssistant.delcache(cfile)
    @test !isfile(cfile)

end
cache_tests()


text:
using PyCall

export quil_run

pushfirst!(PyVector(pyimport("sys")["path"]), ".")
qr = pyimport("quil_run")
const quil_run = qr.quil_run


text:
"""
$(TYPEDEF)

Stores variables related to the dynamics of a `Mechanism`, e.g. the
`Mechanism`'s mass matrix and joint acceleration vector.

Type parameters:
* `T`: the scalar type of the dynamics-related variables.
* `M`: the scalar type of the `Mechanism`.
"""
mutable struct DynamicsResult{T, M}
    mechanism::Mechanism{M}

    massmatrix::Symmetric{T, Matrix{T}}
    dynamicsbias::SegmentedVector{JointID, T, Base.OneTo{JointID}, Vector{T}}
    constraintjacobian::Matrix{T}
    constraintbias::SegmentedVector{JointID, T, UnitRange{JointID}, Vector{T}}
    constraintrowranges::IndexDict{JointID, UnitRange{JointID}, UnitRange{Int}}

    q̇::SegmentedVector{JointID, T, Base.OneTo{JointID}, Vector{T}}
    v̇::SegmentedVector{JointID, T, Base.OneTo{JointID}, Vector{T}}
    ṡ::Vector{T}
    λ::Vector{T}

    contactwrenches::BodyDict{Wrench{T}}
    totalwrenches::BodyDict{Wrench{T}}
    accelerations::BodyDict{SpatialAcceleration{T}}
    jointwrenches::BodyDict{Wrench{T}} # TODO: index by joint tree index?
    contact_state_derivatives::BodyDict{Vector{Vector{DefaultSoftContactStateDeriv{T}}}}

    # see solve_dynamics! for meaning of the following variables:
    L::Matrix{T} # lower triangular
    A::Matrix{T} # symmetric
    z::Vector{T}
    Y::Matrix{T}

    function DynamicsResult{T}(mechanism::Mechanism{M}) where {T, M}
        nq = num_positions(mechanism)
        nv = num_velocities(mechanism)
        nc = num_constraints(mechanism)

        massmatrix = Symmetric(Matrix{T}(undef, nv, nv), :L)
        dynamicsbias = SegmentedVector(Vector{T}(undef, nv), tree_joints(mechanism), num_velocities)
        constraintjacobian = Matrix{T}(undef, nc, nv)
        constraintbias = SegmentedVector{JointID, T, UnitRange{JointID}}(
            Vector{T}(undef, nc), non_tree_joints(mechanism), num_constraints)
        constraintrowranges = ranges(constraintbias)

        q̇ = SegmentedVector(Vector{T}(undef, nq), tree_joints(mechanism), num_positions)
        v̇ = SegmentedVector(Vector{T}(undef, nv), tree_joints(mechanism), num_velocities)
        ṡ = Vector{T}(undef, num_additional_states(mechanism))
        λ = Vector{T}(undef, nc)

        rootframe = root_frame(mechanism)
        contactwrenches = BodyDict{Wrench{T}}(b => zero(Wrench{T}, rootframe) for b in bodies(mechanism))
        totalwrenches = BodyDict{Wrench{T}}(b => zero(Wrench{T}, rootframe) for b in bodies(mechanism))
        accelerations = BodyDict{SpatialAcceleration{T}}(
            b => zero(SpatialAcceleration{T}, rootframe, rootframe, rootframe) for b in bodies(mechanism))
        jointwrenches = BodyDict{Wrench{T}}(b => zero(Wrench{T}, rootframe) for b in bodies(mechanism))
        contact_state_derivs = BodyDict{Vector{Vector{DefaultSoftContactStateDeriv{T}}}}(
            b => Vector{Vector{DefaultSoftContactStateDeriv{T}}}() for b in bodies(mechanism))
        startind = Ref(1)
        for body in bodies(mechanism)
            for point::DefaultContactPoint{M} in contact_points(body)
                model = contact_model(point)
                n = num_states(model)
                push!(contact_state_derivs[body], collect(begin
                    ṡ_part = view(ṡ, startind[] : startind[] + n - 1)
                    contact_state_deriv = SoftContactStateDeriv(model, ṡ_part, root_frame(mechanism))
                    startind[] += n
                    contact_state_deriv
                end for j = 1 : length(mechanism.environment)))
            end
        end

        L = Matrix{T}(undef, nv, nv)
        A = Matrix{T}(undef, nc, nc)
        z = Vector{T}(undef, nv)
        Y = Matrix{T}(undef, nc, nv)

        new{T, M}(mechanism, massmatrix, dynamicsbias, constraintjacobian, constraintbias, constraintrowranges,
            q̇, v̇, ṡ, λ, contactwrenches, totalwrenches, accelerations, jointwrenches, contact_state_derivs,
            L, A, z, Y)
    end
end

DynamicsResult(mechanism::Mechanism{M}) where {M} = DynamicsResult{M}(mechanism)

function Base.copyto!(ẋ::AbstractVector, result::DynamicsResult)
    nq = length(result.q̇)
    nv = length(result.v̇)
    ns = length(result.ṡ)
    @boundscheck length(ẋ) == nq + nv + ns || throw(DimensionMismatch())
    @inbounds copyto!(ẋ, 1, result.q̇, 1, nq)
    @inbounds copyto!(ẋ, nq + 1, result.v̇, 1, nv)
    @inbounds copyto!(ẋ, nq + nv + 1, result.ṡ, 1, ns)
    ẋ
end

contact_state_derivatives(result::DynamicsResult, body::Union{<:RigidBody, BodyID}) =
    result.contact_state_derivatives[body]
contact_wrench(result::DynamicsResult, body::Union{<:RigidBody, BodyID}) =
    result.contactwrenches[body]
set_contact_wrench!(result::DynamicsResult, body::Union{<:RigidBody, BodyID}, wrench::Wrench) =
    (result.contactwrenches[body] = wrench)
acceleration(result::DynamicsResult, body::Union{<:RigidBody, BodyID}) =
    result.accelerations[body]
set_acceleration!(result::DynamicsResult, body::Union{<:RigidBody, BodyID}, accel::SpatialAcceleration) =
    (result.accelerations[body] = accel)
joint_wrench(result::DynamicsResult, body::Union{<:RigidBody, BodyID}) =
    result.jointwrenches[body]
set_joint_wrench!(result::DynamicsResult, body::Union{<:RigidBody, BodyID}, wrench::Wrench) =
    (result.jointwrenches[body] = wrench)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct VirtualMachineImage <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    location::Any # spec type: Union{ Nothing, String } # spec name: location
    tags::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: tags
    properties::Any # spec type: Union{ Nothing, VirtualMachineImageProperties } # spec name: properties

    function VirtualMachineImage(;id=nothing, name=nothing, location=nothing, tags=nothing, properties=nothing)
        o = new()
        validate_property(VirtualMachineImage, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(VirtualMachineImage, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(VirtualMachineImage, Symbol("location"), location)
        setfield!(o, Symbol("location"), location)
        validate_property(VirtualMachineImage, Symbol("tags"), tags)
        setfield!(o, Symbol("tags"), tags)
        validate_property(VirtualMachineImage, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type VirtualMachineImage

const _property_map_VirtualMachineImage = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("name")=>Symbol("name"), Symbol("location")=>Symbol("location"), Symbol("tags")=>Symbol("tags"), Symbol("properties")=>Symbol("properties"))
const _property_types_VirtualMachineImage = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("name")=>"String", Symbol("location")=>"String", Symbol("tags")=>"Dict{String, String}", Symbol("properties")=>"VirtualMachineImageProperties")
Base.propertynames(::Type{ VirtualMachineImage }) = collect(keys(_property_map_VirtualMachineImage))
Swagger.property_type(::Type{ VirtualMachineImage }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_VirtualMachineImage[name]))}
Swagger.field_name(::Type{ VirtualMachineImage }, property_name::Symbol) =  _property_map_VirtualMachineImage[property_name]

function check_required(o::VirtualMachineImage)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("location")) === nothing) && (return false)
    true
end

function validate_property(::Type{ VirtualMachineImage }, name::Symbol, val)
end


text:
using Documenter, ControlSystems, Plots, LinearAlgebra, DSP
import GR # Bug with world age in Plots.jl, see https://github.com/JuliaPlots/Plots.jl/issues/1047

include("src/makeplots.jl")

makedocs(modules=[ControlSystems])

# If not running travis, generate the plots here, even if we are not deploying
if get(ENV, "TRAVIS", "") == ""
    makePlots()
end

# Only build plots in travis if we are deploying
# And dont install the dependencies unless we are deploying
function myDeps()
    if get(ENV, "TRAVIS", "") != ""
        println("Installing deploy dependencies")
        run(`pip install --user pygments mkdocs`)
        makePlots()
    end
end

deploydocs(
    repo = "github.com/JuliaControl/ControlSystems.jl.git",
    latest = "master",
    julia = "1.0",
    deps = myDeps
)


text:
using Wavelets
using Images  # for imread and imwrite

img = imread("lena.tiff")
x = permutedims(img.data, [ndims(img.data):-1:1])
L = 2
xts = wplotim(x, L, waveletfilter("db3"))

imwrite(xts, "lena_2d.png")
# convert with ImageMagick
run(`convert lena_2d.png transform2d_lena.jpg`)


text:
function get_signatures(b, sigs, server) end

function get_signatures(b::StaticLint.Binding, sigs, server)
    if b.type == StaticLint.CoreTypes.Function
        if b.val isa EXPR && CSTParser.defines_function(b.val)
            sig = CSTParser.rem_where_decl(CSTParser.get_sig(b.val))
            args = []
            if sig.args !== nothing
                for i = 2:length(sig.args)
                    if bindingof(sig.args[i]) !== nothing
                        push!(args, bindingof(sig.args[i]))
                    end
                end
            end
            params = (a->ParameterInformation(valof(a.name) isa String ? valof(a.name) : "", missing)).(args)
            push!(sigs, SignatureInformation(string(Expr(sig)), "", params))
        end
    end
end

JSONRPC.parse_params(::Type{Val{Symbol("textDocument/signatureHelp")}}, params) = TextDocumentPositionParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/signatureHelp")},TextDocumentPositionParams}, server)
    doc = getdocument(server, URI2(r.params.textDocument.uri))
    sigs = SignatureInformation[]
    offset = get_offset(doc, r.params.position)
    rng = Range(doc, offset:offset)
    x = get_expr(getcst(doc), offset)
    arg = 0

    if x isa EXPR && parentof(x) isa EXPR && typof(parentof(x)) === CSTParser.Call
        if typof(parentof(x).args[1]) === CSTParser.IDENTIFIER
            call_name = parentof(x).args[1]
        elseif typof(parentof(x).args[1]) === CSTParser.Curly && typof(parentof(x).args[1].args[1]) === CSTParser.IDENTIFIER
            call_name = parentof(x).args[1].args[1]
        elseif typof(parentof(x).args[1]) === CSTParser.BinaryOpCall && kindof(parentof(x).args[1].args[2]) === CSTParser.Tokens.DOT && length(parentof(x).args[1].args) == 3 && length(parentof(x).args[1].args[3]) == 1
            call_name = parentof(x).args[1].args[3].args[1]
        else
            call_name = nothing
        end
        if call_name !== nothing && StaticLint.hasref(call_name)
            f_binding = refof(call_name)
            while f_binding isa StaticLint.Binding || f_binding isa SymbolServer.FunctionStore
                if f_binding isa StaticLint.Binding && f_binding.type == StaticLint.CoreTypes.Function
                    get_signatures(f_binding, sigs, server)
                    f_binding = f_binding.prev
                elseif refof(call_name) isa SymbolServer.FunctionStore
                    for m in refof(call_name).methods
                        sig = string(call_name.val, "(", join([a[2] for a in m.args], ", "),")")
                        params = (a->ParameterInformation(a[1], missing)).(m.args)
                        push!(sigs, SignatureInformation(sig, "", params))
                    end
                    break
                else
                    break
                end
            end
        end
    end
    if (isempty(sigs) || (typof(x) === CSTParser.PUNCTUATION  && kindof(x) === CSTParser.Tokens.RPAREN))
        return SignatureHelp(SignatureInformation[], 0, 0)
    end

    if typof(x) === CSTParser.Tokens.LPAREN
        arg = 0
    else
        arg = sum(!(typof(a) === CSTParser.PUNCTUATION) for a in parentof(x).args) - 1
    end
    return SignatureHelp(filter(s->length(s.parameters) > arg, sigs), 0, arg)
end

JSONRPC.parse_params(::Type{Val{Symbol("textDocument/definition")}}, params) = TextDocumentPositionParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/definition")},TextDocumentPositionParams}, server)
    locations = Location[]
    doc = getdocument(server, URI2(r.params.textDocument.uri))
    offset = get_offset(doc, r.params.position)
    x = get_expr1(getcst(doc), offset)
    if x isa EXPR && StaticLint.hasref(x)
        b = refof(x)
        if b isa SymbolServer.FunctionStore || b isa SymbolServer.DataTypeStore
            for m in b.methods
                if isfile(m.file)
                    push!(locations, Location(filepath2uri(m.file), Range(m.line - 1, 0, m.line -1, 0)))
                end
            end
        elseif b isa StaticLint.Binding && b.val isa StaticLint.Binding
            b = b.val
        end

        while b isa StaticLint.Binding
            if b.val isa EXPR
                doc1, o = get_file_loc(b.val)
                if doc1 isa Document
                    push!(locations, Location(doc1._uri, Range(doc1, o .+ (0:b.val.span))))
                end
            elseif b.val isa SymbolServer.FunctionStore
                for m in b.val.methods
                    file = isabspath(string(m.file)) ? string(m.file) : Base.find_source_file(string(m.file))
                    ((file, m.line) == DefaultTypeConstructorLoc || file == nothing) && continue
                    push!(locations, Location(filepath2uri(file), Range(m.line - 1, 0, m.line, 0)))
                end
            end
            if b.type == StaticLint.CoreTypes.Function && b.prev isa StaticLint.Binding && (b.prev.type == Function || b.prev.type == StaticLint.CoreTypes.DataType)
                b = b.prev
            else
                b = nothing
            end
        end
    elseif x isa EXPR && typof(x) === CSTParser.LITERAL && (kindof(x) === Tokens.STRING || kindof(x) === Tokens.TRIPLE_STRING)
        if sizeof(valof(x)) < 256 # AUDIT: OK
            if isfile(valof(x))
                push!(locations, Location(filepath2uri(valof(x)), Range(0, 0, 0, 0)))
            elseif isfile(joinpath(dirname(uri2filepath(doc._uri)), valof(x)))
                push!(locations, Location(filepath2uri(joinpath(dirname(uri2filepath(doc._uri)), valof(x))), Range(0, 0, 0, 0)))
            end
        end
    end

    return locations
end

function get_file_loc(x::EXPR, offset = 0, c  = nothing)
    if c !== nothing
        for a in x.args
            a == c && break
            offset += a.fullspan
        end
    end
    if parentof(x) !== nothing
        return get_file_loc(parentof(x), offset, x)
    elseif typof(x) === CSTParser.FileH && StaticLint.hasmeta(x)
        # return x.ref, offset
        return x.meta.error, offset
    else
        return nothing, offset
    end
end

JSONRPC.parse_params(::Type{Val{Symbol("textDocument/formatting")}}, params) = DocumentFormattingParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/formatting")},DocumentFormattingParams}, server::LanguageServerInstance)
    doc = getdocument(server, URI2(r.params.textDocument.uri))
    newcontent = DocumentFormat.format(get_text(doc), server.format_options)
    end_l, end_c = get_position_at(doc, sizeof(get_text(doc))) # AUDIT: OK
    lsedits = TextEdit[TextEdit(Range(0, 0, end_l, end_c), newcontent)]

    return lsedits
end

JSONRPC.parse_params(::Type{Val{Symbol("textDocument/documentLink")}}, params) = DocumentLinkParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/documentLink")},DocumentLinkParams}, server)
    if !hasdocument(server, URI2(r.params.textDocument.uri))
        error("Received 'textDocument/documentLink for non-existing document.")
    end
    links = Tuple{String,UnitRange{Int}}[]

    return links
end

function find_references(textDocument::TextDocumentIdentifier, position::Position, server)
    locations = Location[]
    doc = getdocument(server, URI2(textDocument.uri))
    offset = get_offset(doc, position)
    x = get_expr1(getcst(doc), offset)
    if x isa EXPR && StaticLint.hasref(x) && refof(x) isa StaticLint.Binding
        for r in refof(x).refs
            !(r isa EXPR) && continue
            doc1, o = get_file_loc(r)
            if doc1 isa Document
                push!(locations, Location(doc1._uri, Range(doc1, o .+ (0:r.span))))
            end
        end
    end
    return locations
end

JSONRPC.parse_params(::Type{Val{Symbol("textDocument/references")}}, params) = ReferenceParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/references")},ReferenceParams}, server)
    locations = find_references(r.params.textDocument, r.params.position, server)
    return locations
end

JSONRPC.parse_params(::Type{Val{Symbol("textDocument/rename")}}, params) = RenameParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/rename")},RenameParams}, server)
    tdes = Dict{String,TextDocumentEdit}()
    locations = find_references(r.params.textDocument, r.params.position, server)

    for loc in locations
        if loc.uri in keys(tdes)
            push!(tdes[loc.uri].edits, TextEdit(loc.range, r.params.newName))
        else
            doc = getdocument(server, URI2(loc.uri))
            tdes[loc.uri] = TextDocumentEdit(VersionedTextDocumentIdentifier(loc.uri, doc._version), [TextEdit(loc.range, r.params.newName)])
        end
    end

    return WorkspaceEdit(nothing, collect(values(tdes)))
end

JSONRPC.parse_params(::Type{Val{Symbol("textDocument/documentSymbol")}}, params) = DocumentSymbolParams(params)
function process(r::JSONRPC.Request{Val{Symbol("textDocument/documentSymbol")},DocumentSymbolParams}, server)
    syms = SymbolInformation[]
    uri = r.params.textDocument.uri
    doc = getdocument(server, URI2(uri))

    bs = collect_bindings_w_loc(getcst(doc))
    for x in bs
        p,b = x[1], x[2]
        !(b.val isa EXPR) && continue
        (valof(b.name) === nothing || isempty(valof(b.name))) && typof(b.name) !== CSTParser.OPERATOR && continue
        push!(syms, SymbolInformation(typof(b.name) === CSTParser.OPERATOR ? string(Expr(b.name)) : valof(b.name), _binding_kind(b, server), false, Location(doc._uri, Range(doc, p)), nothing))
    end
    return syms
end

function collect_bindings_w_loc(x::EXPR, pos = 0, bindings = Tuple{UnitRange{Int},StaticLint.Binding}[])
    if bindingof(x) !== nothing
        push!(bindings, (pos .+ (0:x.span), bindingof(x)))
    end
    if x.args !== nothing
        for a in x.args
            collect_bindings_w_loc(a, pos, bindings)
            pos += a.fullspan
        end
    end
    return bindings
end

function collect_toplevel_bindings_w_loc(x::EXPR, pos = 0, bindings = Tuple{UnitRange{Int},StaticLint.Binding}[]; query = "")
    # if bindingof(x) isa StaticLint.Binding && !isempty(valof(bindingof(x).name)) && bindingof(x).val isa EXPR && startswith(valof(bindingof(x).name), query)
    if bindingof(x) isa StaticLint.Binding && valof(bindingof(x).name) isa String && bindingof(x).val isa EXPR && startswith(valof(bindingof(x).name), query)
        push!(bindings, (pos .+ (0:x.span), bindingof(x)))
    end
    if scopeof(x) !== nothing && !(typof(x) === CSTParser.FileH || typof(x) === CSTParser.ModuleH || typof(x) === CSTParser.BareModule)
        return bindings
    end
    if x.args !== nothing
        for a in x.args
            collect_toplevel_bindings_w_loc(a, pos, bindings, query = query)
            pos += a.fullspan
        end
    end
    return bindings
end

function _binding_kind(b ,server)
    if b isa StaticLint.Binding
        if b.type == nothing
            return 13
        elseif b.type == StaticLint.CoreTypes.Module
            return 2
        elseif b.type == StaticLint.CoreTypes.Function
            return 12
        elseif b.type == StaticLint.CoreTypes.String
            return 15
        elseif b.type == StaticLint.CoreTypes.Int || b.type == StaticLint.CoreTypes.Float64
            return 16
        elseif b.type == StaticLint.CoreTypes.DataType
            return 23
        else
            return 13
        end
    elseif b isa SymbolServer.ModuleStore
        return 2
    elseif b isa SymbolServer.MethodStore
        return 6
    elseif b isa SymbolServer.FunctionStore
        return 12
    elseif b isa SymbolServer.DataTypeStore
        return 23
    else
        return 13
    end
end


text:
using ScatteringLayers
using CSVFiles
using DataFrames
using Dates
using DiffEqOperators
using Images
using StaticArrays

println("Testing echogram types...")
sp = SpacePoint(SVector(randn(3)...))
sp = SpacePoint(randn(3))
stp = SpaceTimePoint(SVector(randn(3)...), Dates.now())
stp = SpaceTimePoint(randn(3), Dates.now())
sps = [SpacePoint(SVector(randn(3)...)) for i in 1:10]
z = collect(1:5.0)
eg = Echogram(randn(5, 10), z, sps)
depths(eg)
spacetime(eg)
png = getping(eg, 1)

println("Loading test data...")
echo_df = DataFrame(load("data\\test_data.csv"))
fmt = DateFormat("yyyymmddHH:MM:SS.s")
echo_df[:datetime] = DateTime.(string.(echo_df[:Date_M]) .* echo_df[:Time_M], fmt)

echo = unstack(echo_df, :Layer_depth_max, :datetime, :Sv_mean)#[1:end, 2:end]
z = float.(echo[:Layer_depth_max])
t = DateTime.(string.(names(echo)[2:end]))
echo = disallowmissing(Matrix(echo[:, 2:end]))
echo = echo[end:-1:1, :]
echo[echo .< -90] .= -90
echo = 10 .^(echo/10)
echo = imfilter(echo, Kernel.gaussian(0.5))
stp = [SpaceTimePoint([36.712110, -122.187028], ti) for ti in t]
echo = Echogram(echo, z, stp)
m, n = size(echo)
dz = 1.0
zz = 880.0:-dz:8.0
order = 1
D2 = DerivativeOperator{Float64}(2, order, dz, m, :Dirichlet0, :Dirichlet0)

println("Testing...")
ping = getping(echo, 5)

THRESH = 10^(-80/10)
μ, p, c = getpeaks(ping, thresh=THRESH)
mix = guesslayers(ping, thresh=THRESH)

basis = getbasis(mix, zz)
N1 = refine_heights!(mix, basis, ping)
ping_nnls = mixpdf(mix, zz)

fit1 = refine_μσ!(mix, ping)
ping_refined1 = mixpdf(mix, zz)
mix = guesslayers(ping, thresh=THRESH)
fit!(mix, ping, trace=true)
ping_refined2 = mixpdf(mix, zz)

mix1 = fit(mix, ping, trace=true)
mix1 = fit(mix, ping, refine=:σ)
mix1 = fit(mix, ping, refine=:σμ)
mix1 = fit(mix, ping, refine=:widths)
mix1 = fit(mix, ping, refine=:widths_locs)

x = SpacePoint([36.712110, -122.187028])
echogram = Echogram(imfilter(echo, Kernel.gaussian(0.5)), zz, repeat([x], n))
layers = fitlayers(echogram, thresh=THRESH)

max_stars_repo_name: ElOceanografo/ScatteringLayers

max_issues_repo_name: ElOceanografo/ScatteringLayers

max_forks_repo_name: ElOceanografo/ScatteringLayers


text:
# N.B. This file is maintained in TensorFlowBuilder source
using PyCall
import Base: convert
import PyCall: PyObject

export SessionRunnable, AbstractTensor, DimsType, NoneDimsType, NegDimsType, PyVectorType

macro pywrapper(typename, supertype)
  quote
    immutable $typename <: $supertype
      x::PyCall.PyObject
    end
    PyCall.PyObject(o::$typename) = o.x
    isequal(x::$typename, y::$typename) = x.x[__eq__](y.x)
    export $typename
  end
end

abstract SessionRunnable
abstract AbstractTensor <: SessionRunnable

@pywrapper Dtype Any

@pywrapper Tensor AbstractTensor
@pywrapper Variable AbstractTensor
@pywrapper SparseTensor AbstractTensor
@pywrapper Placeholder AbstractTensor

@pywrapper Operation SessionRunnable

@pywrapper Session Any

@pywrapper Optimizer Any

convert(::Type{Tensor}, a::Array) = Tensor(a)
convert(::Type{AbstractTensor}, a::Array) = Tensor(a)

## DimsType ##
# TensorFlow shape arguments sometimes want None to mean "wild" (e.g. tf.placeholder)
# and sometimes use negatives to mean wildcard (e.g.reshape)

abstract DimsType

"A literal used for TensorFlow shape arguments.
  Negative indices are passed as Python `None` values."
type NoneDimsType <: DimsType
  arr::Vector{Int}
end
NoneDimsType(a::Tuple{Vararg{Int}}) = NoneDimsType(collect(a))
PyCall.PyObject(o::NoneDimsType) = PyCall.PyObject([a < 0 ? nothing : a for a in o.arr])

"A literal used for TensorFlow shape arguments.
  Negative indices are passed through to Python."
type NegDimsType <: DimsType
  arr::Vector{Int}
end
NegDimsType(a::Tuple{Vararg{Int}}) = NegDimsType(collect(a))
PyCall.PyObject(o::NegDimsType) = PyCall.PyObject(o.arr)

"A wrapper type for things that TensorFlow insists must be Python
  lists, like conv2d's strides argument."
type PyVectorType <: DimsType
  pyvec::PyVector
end
PyVectorType(a::Vector) = PyVectorType(PyVector(a))
PyVectorType(a::Tuple) = PyVectorType(collect(a))
PyCall.PyObject(o::PyVectorType) = PyCall.PyObject(o.pyvec)


text:
using OceanParameterizations
using WindMixing
using OceanTurb

file_labels = Dict(
    "free_convection" => "Free convection",
    "strong_wind" => "Strong wind",
    "strong_wind_no_coriolis" => "Strong wind, no rotation",
    "weak_wind_strong_cooling" => "Weak wind, strong cooling",
    "strong_wind_weak_cooling" => "Strong wind, weak cooling",
    "strong_wind_weak_heating" => "Strong wind, weak heating"
)

files =  ["free_convection", "strong_wind", "strong_wind_no_coriolis",
            "weak_wind_strong_cooling", "strong_wind_weak_cooling", "strong_wind_weak_heating"]

for test_file in files[1:2]
    𝒟test = data(test_file,
                        override_scalings=𝒟train.scalings, # use the scalings from the training data
                        animate=false,
                        animate_dir=output_directory*"/Testing")
    les = read_les_output(test_file)

    output_gif_directory="KPP_TKE/test_$(test_file)"
    mkpath(output_gif_directory)
    file = output_gif_directory*"_output.txt"
    touch(file)
    o = open(file, "w")

    write(o, "= = = = = = = = = = = = = = = = = = = = = = = = \n")
    write(o, "Test file: $(test_file) \n")
    write(o, "Output will be written to: $(output_gif_directory) \n")

    ## KPP Parameterization (no training)

    # les = read_les_output(test_file)
    parameters = KPP.Parameters() # default parameters
    KPP_model = closure_kpp_full_evolution(parameters, 𝒟test.T.coarse[:,1], les)
    predictions = KPP_model()
    T_KPP = (predictions, 𝒟test.T.coarse)
    write(o, "KPP prediction error on T........ $(mse(T_KPP)) \n")
    animate_prediction(T_KPP, "T", 𝒟test, test_file; legend_labels=["KPP(T)", "truth"], filename="T_KPP_$(test_file)", directory=output_directory)

    ## TKE Parameterization (no training; use default parameters)

    # les = read_les_output(test_file)
    parameters = TKEMassFlux.TKEParameters() # default parameters
    TKE_model = closure_tke_full_evolution(parameters, 𝒟test.T.coarse[:,1], les)
    predictions = TKE_model()
    T_TKE = (predictions, 𝒟test.T.coarse)
    write(o, "KPP prediction error on T........ $(mse(T_TKE)) \n")
    animate_prediction(T_TKE, "T", 𝒟test, test_file; legend_labels=["TKE(T)", "truth"], filename="T_TKE_$(test_file)", directory=output_directory)

    close(o)
end


text:
################################################################################
#                             Euler Angle and Axis
################################################################################

export angleaxis_to_angle, angleaxis_to_dcm, angleaxis_to_quat

################################################################################
#                                  Operations
################################################################################

"""
    function *(ea₂::EulerAngleAxis{T1}, ea₁::EulerAngleAxis{T2}) where {T1,T2}

Compute the composed rotation of `ea₁ -> ea₂`. Notice that the rotation will be
represented by a Euler angle and axis (see `EulerAngleAxis`). By convention, the
output angle will always be in the range [0, π] [rad].

Notice that the vector representing the axis in `ea₁` and `ea₂` must be unitary.
This function neither verifies this nor normalizes the vector.

"""
function *(ea₂::EulerAngleAxis{T1}, ea₁::EulerAngleAxis{T2}) where {T1,T2}
    # Auxiliary variables.
    sθ₁o2, cθ₁o2 = sincos(ea₁.a/2)
    sθ₂o2, cθ₂o2 = sincos(ea₂.a/2)

    v₁ = ea₁.v
    v₂ = ea₂.v

    # Compute `cos(θ/2)` in which `θ` is the new Euler angle.
    cθo2 = cθ₁o2*cθ₂o2 - sθ₁o2*sθ₂o2 * dot(v₁, v₂)

    T = promote_type( T1, T2, typeof(sθ₁o2), typeof(sθ₂o2), typeof(cθo2) )

    if abs(cθo2) >= 1-eps()
        # In this case, the rotation is the identity.
        return EulerAngleAxis( T(0), SVector{3,T}(0,0,0) )
    else
        # Compute `sin(θ/2)` in which `θ` is the new Euler angle.
        sθo2 = sqrt(1 - cθo2*cθo2)

        # Compute the θ angle between [0, 2π].
        θ = 2acos(cθo2)

        # Keep the angle between [0, π].
        s = +1
        if θ > π
            θ = T(2)*π - θ
            s = -1
        end

        v = s*( sθ₁o2*cθ₂o2*v₁ + cθ₁o2*sθ₂o2*v₂ + sθ₁o2*sθ₂o2*(v₁ × v₂) )/sθo2

        return EulerAngleAxis(θ, v)
    end
end

"""
    @inline function inv(ea::EulerAngleAxis)

Compute the inverse rotation of `ea`. The Euler angle returned by this function
will always be in the interval [0, π].

"""
@inline function inv(ea::EulerAngleAxis{T}) where T<:Number
    # Make sure that the Euler angle is always in the inverval [0,π]
    s = -1
    θ = mod(ea.a, T(2)*π)

    if θ > π
        s = 1
        θ = T(2)π - θ
    end

    EulerAngleAxis(θ, s*ea.v)
end

################################################################################
#                                      IO
################################################################################

"""
    function display(ea::EulerAngleAxis{T}) where T
    function show(io::IO, mime::MIME"text/plain", ea::EulerAngleAxis{T}) where T

Display in `stdout` the Euler angle and axis `ea`.

"""
function show(io::IO, ea::EulerAngleAxis{T}) where T
    θ   = ea.a
    v   = ea.v
    str = @sprintf "%8.4f rad, [%8.4f, %8.4f, %8.4f]" ea.a v[1] v[2] v[3]

    print(io, "EulerAngleAxis{$T}: $str")

    nothing
end

function show(io::IO, mime::MIME"text/plain", ea::EulerAngleAxis{T}) where T
    # Check if the user wants colors.
    color = get(io, :color, false)

    b = (color) ? "\x1b[1m"         : ""
    d = (color) ? "\x1b[0m"         : ""
    g = (color) ? "\x1b[1m\x1b[32m" : ""
    y = (color) ? "\x1b[1m\x1b[33m" : ""

    str_a = @sprintf "%8.4f rad (%8.4f deg)" ea.a rad2deg(ea.a)
    str_v = @sprintf "[%8.4f, %8.4f, %8.4f]" ea.v[1] ea.v[2] ea.v[3]

    println(io, "EulerAngleAxis{$T}:")
    println(io, "$y  Euler angle: $d" * str_a)
      print(io, "$y   Euler axis: $d" * str_v)

    nothing
end

################################################################################
#                                 Conversions
################################################################################

# Euler Angles
# ==============================================================================

"""
    @inline function angleaxis_to_angle(θ::Number, v::AbstractVector, rot_seq::Symbol)
    @inline function angleaxis_to_angle(ea::EulerAngleAxis, rot_seq::Symbol)

Convert the Euler angle `θ` [rad]  and Euler axis `v`, which must be a unit
vector, to Euler angles with rotation sequence `rot_seq`. Those values can also
be passed inside the structure `ea` (see `EulerAngleAxis`).

The rotation sequence is defined by a `:Symbol`. The possible values are:
`:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`, `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`,
`:ZYX`, and `:ZYZ`. If no value is specified, then it defaults to `:ZYX`.

# Example

```julia-repl
julia> ea = EulerAngleAxis(45*pi/180, [1;0;0]);

julia> angleaxis_to_angles(ea, :ZXY)
EulerAngles{Float64}:
  R(Z):   0.0000 rad (   0.0000 deg)
  R(X):   0.7854 rad (  45.0000 deg)
  R(Y):   0.0000 rad (   0.0000 deg)

```
"""
@inline function angleaxis_to_angle(θ::Number, v::AbstractVector, rot_seq::Symbol)
    # First we convert to quaternions then to Euler angles.
    quat_to_angle( angleaxis_to_quat(θ, v), rot_seq )
end

@inline angleaxis_to_angle(ea::EulerAngleAxis, rot_seq::Symbol) =
    angleaxis_to_angle(ea.a, ea.v, rot_seq)

# DCM
# ==============================================================================

"""
    @inline function angleaxis_to_dcm(a::Number, v::AbstractVector)
    @inline function angleaxis_to_dcm(ea::EulerAngleAxis)

Convert the Euler angle `a` [rad] and Euler axis `v`, which must be a unit
vector to a DCM. Those values can also be passed inside the structure `ea` (see
`EulerAngleAxis`).

# Remarks

It is expected that the vector `v` is unitary. However, no verification is
performed inside the function. The user must handle this situation.

# Example

```julia-repl
julia> v = [1;1;1];

julia> v /= norm(v);

julia> angleaxis_to_dcm(pi/2,v)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  0.333333   0.910684  -0.244017
 -0.244017   0.333333   0.910684
  0.910684  -0.244017   0.333333

julia> ea = EulerAngleAxis(pi/2,v);

julia> angleaxis_to_dcm(ea)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
  0.333333   0.910684  -0.244017
 -0.244017   0.333333   0.910684
  0.910684  -0.244017   0.333333
```

"""
@inline function angleaxis_to_dcm(a::Number, v::AbstractVector)
    sθ, cθ = sincos(a)
    aux    = 1 - cθ

    DCM(  cθ + v[1]*v[1]*aux   , v[1]*v[2]*aux + v[3]*sθ, v[1]*v[3]*aux - v[2]*sθ,
        v[1]*v[2]*aux - v[3]*sθ,    cθ + v[2]*v[2]*aux  , v[2]*v[3]*aux + v[1]*sθ,
        v[1]*v[3]*aux + v[2]*sθ, v[2]*v[3]*aux - v[1]*sθ,   cθ + v[3]*v[3]*aux)'
end

@inline angleaxis_to_dcm(ea::EulerAngleAxis) = angleaxis_to_dcm(ea.a, ea.v)

# Quaternions
# ==============================================================================

"""
    function angleaxis_to_quat(θ::Number, v::AbstractVector)

Convert the Euler angle `θ` [rad] and Euler axis `v`, which must be a unit
vector, to a quaternion.

# Remarks

It is expected that the vector `v` is unitary. However, no verification is
performed inside the function. The user must handle this situation.

# Example

```julia-repl
julia> v = [1;1;1];

julia> v /= norm(v);

julia> angleaxis_to_quat(pi/2,v)
Quaternion{Float64}:
  + 0.7071067811865476 + 0.408248290463863.i + 0.408248290463863.j + 0.408248290463863.k
```

"""
function angleaxis_to_quat(θ::Number, v::AbstractVector)
    # Check the arguments.
    (length(v) > 3) && throw(ArgumentError("The provided vector for the Euler axis must have 3 elements."))

    cθo2 = cos(θ/2)

    # Keep `q0` positive.
    s = (cθo2 < 0) ? -1 : +1

    # Create the quaternion.
    Quaternion( s*cos(θ/2), s*sin(θ/2)*v )
end

"""
    function angleaxis_to_quat(angleaxis::EulerAngleAxis)

Convert a Euler angle and Euler axis `angleaxis` (see `EulerAngleAxis`) to a
quaternion.

# Remarks

It is expected that the vector `angleaxis.v` is unitary. However, no
verification is performed inside the function. The user must handle this
situation.

# Example

```julia-repl
julia> v = [1;1;1];

julia> v /= norm(v);

julia> angleaxis_to_quat(EulerAngleAxis(pi/2,v))
Quaternion{Float64}:
  + 0.7071067811865476 + 0.408248290463863.i + 0.408248290463863.j + 0.408248290463863.k
```

"""
angleaxis_to_quat(angleaxis::EulerAngleAxis) = angleaxis_to_quat(angleaxis.a,
                                                                 angleaxis.v)


text:
struct BagNode{T <: Union{Missing, AbstractNode}, B <: AbstractBags, C} <: AbstractBagNode
    data::T
    bags::B
    metadata::C

    function BagNode(d::T, b::B, m::C=nothing) where {T <: Union{Missing, AbstractNode}, B <: AbstractBags, C}
        ismissing(d) && any(_len.(b.bags) .!= 0) && error("BagNode with nothing in data cannot have a non-empty bag")
        new{T, B, C}(d, b, m)
    end
end

BagNode(data::T, b::Vector, metadata::C=nothing) where {T, C} = BagNode(data, bags(b), metadata)

Flux.@functor BagNode

mapdata(f, x::BagNode) = BagNode(mapdata(f, x.data), x.bags, x.metadata)

function Base.getindex(x::BagNode, i::VecOrRange)
    nb, ii = remapbag(x.bags, i)
    _emptyismissing[] && isempty(ii) && return(BagNode(missing, nb, nothing))
    BagNode(subset(x.data,ii), nb, subset(x.metadata, i))
end

function reduce(::typeof(catobs), as::Vector{T}) where {T <: BagNode}
    data = filter(!ismissing, [x.data for x in as])
    metadata = filter(!isnothing, [x.metadata for x in as])
    bags = _catbags([d.bags for d in as])
    BagNode(isempty(data) ? missing : reduce(catobs, data),
            bags,
            isempty(metadata) ? nothing : reduce(catobs, metadata))
end

removeinstances(a::BagNode, mask) = BagNode(subset(a.data, findall(mask)), adjustbags(a.bags, mask), a.metadata)

Base.hash(e::BagNode{T,B,C}, h::UInt) where {T,B,C} = hash((T, B, C, e.data, e.bags, e.metadata), h)
Base.:(==)(e1::BagNode{T,B,C}, e2::BagNode{T,B,C}) where {T,B,C} = e1.data == e2.data && e1.bags == e2.bags && e1.metadata == e2.metadata