text:
"""
    Sphere([origin=zeros(),] radius)

A [`Shape`](@ref) where boundary is a fixed distance from the origin. In 2D this is a circle, in 3D the usual sphere, and in higher dimensions if difficult to visualise.
"""
struct Sphere{T,Dim} <: Shape{T,Dim}
    origin::SVector{Dim,T}
    radius::T
end

# Alternate constructors, where type is inferred naturally
Sphere(origin::NTuple{Dim}, radius::T) where {T,Dim} = Sphere{T,Dim}(origin, radius)
Sphere(origin::AbstractVector, radius::T) where {T} = Sphere{T,length(origin)}(origin, radius)

Sphere(Dim, radius::T) where {T} = Sphere{T,Dim}(zeros(T,Dim), radius)

Circle(origin::Union{AbstractVector{T},NTuple{2,T}}, radius::T) where T <: AbstractFloat = Sphere{T,2}(origin, radius::T)
Circle(radius::T) where T <: AbstractFloat = Sphere(2, radius::T)
Sphere(radius::T) where {T} = Sphere{T,3}(zeros(T,3), radius)

name(shape::Sphere) = "Sphere"
name(shape::Sphere{T,2}) where T = "Circle"
Symmetry(shape::Sphere{T,Dim}) where {T,Dim} = RadialSymmetry{Dim}()

outer_radius(sphere::Sphere) = sphere.radius
volume(shape::Sphere{T,3}) where T = 4//3 * π * shape.radius^3
volume(shape::Sphere{T,2}) where T = π * shape.radius^2

# bounding_box(sphere::Sphere{T,3}; kws...) where T = bounding_box(Circle(sphere; kws...))
function bounding_box(sphere::Sphere{T,Dim}) where {T,Dim}
    return Box(origin(sphere), T(2)*sphere.radius .* ones(T,Dim))
end

import Base.in
function in(x::AbstractVector, sphere::Sphere)::Bool
    norm(origin(sphere) .- x) <= sphere.radius
end

import Base.issubset
function issubset(inner_sphere::Sphere{T,Dim}, outer_sphere::Sphere{T,Dim}) where {T,Dim}
    norm(origin(outer_sphere) - origin(inner_sphere)) <= outer_sphere.radius - inner_sphere.radius
end

function issubset(sphere::Sphere, box::Box)
    sphere_box = bounding_box(sphere)
    return issubset(sphere_box,box)
end

function issubset(box::Box, sphere::Sphere)
    issubset(Sphere(origin(box),outer_radius(box)), sphere)
end

import Base.(==)
function ==(c1::Sphere, c2::Sphere)
    c1.origin == c2.origin &&
    c1.radius == c2.radius
end

import Base.isequal
function isequal(c1::Sphere, c2::Sphere)
    isequal(c1.origin, c2.origin) &&
    isequal(c1.radius, c2.radius)
end

function iscongruent(s1::Sphere, s2::Sphere)
    s1.radius == s2.radius
end

function congruent(s::Sphere, x)
    Sphere(x, s.radius)
end

function Circle(sphere::Sphere{T}; y = sphere.origin[2]) where T
    if abs(y - sphere.origin[2]) > sphere.radius
        return EmptyShape(sphere)
    else
        r = sqrt(sphere.radius^2 - (sphere.origin[2] - y)^2)
        return Sphere{T,2}(sphere.origin[[1,3]], r)
    end
end

# boundary_functions(sphere::Sphere; kws...) = boundary_functions(Circle(sphere; kws...))

function boundary_functions(sphere::Sphere{T,3}) where T

    function x(t,s=T(0))
        check_boundary_coord_range(t)
        check_boundary_coord_range(s)

        sphere.radius * sin(t * π) * cos(s * 2π) + origin(sphere)[1]
    end

    function y(t,s=T(0))
        check_boundary_coord_range(t)
        check_boundary_coord_range(s)

        sphere.radius * sin(t * π) * sin(s * 2π) + origin(sphere)[2]
    end

    function z(t,s=T(0))
        check_boundary_coord_range(t)
        check_boundary_coord_range(s)

        sphere.radius * cos(t * π) + origin(sphere)[3]
    end

    return x, y, z
end

function boundary_functions(circle::Sphere{T,2}) where T

    function x(t)
        check_boundary_coord_range(t)
        circle.radius * cos(2π * t) + origin(circle)[1]
    end

    function y(t)
        check_boundary_coord_range(t)
        circle.radius * sin(2π * t) + origin(circle)[2]
    end

    return x, y
end


text:
module Examples

using ..ConservationLawsParticles: @time_independent

export V, V2, W_attr, W_rep, Wprime_attr, Wprime_rep, W, Wprime, mob, mobρ, mobσ

@time_independent V(r) = - r^3
@time_independent V2(r) = - (r-1)^3
@time_independent W_attr(r) = 5log1p(abs(r))
@time_independent W_rep(r) = -5log1p(abs(r))
@time_independent Wprime_attr(r) = 5 * sign(r) / (abs(r) + 1)
@time_independent Wprime_rep(r) = - 5 * sign(r) / (abs(r) + 1)
@time_independent W(r) = - 1 / (abs(r) + 1)
@time_independent Wprime(r) = sign(r) / (abs(r) + 1)^2
mob(ρ) = max(1 - ρ, 0)
mob(ρ, σ) = max(1 - ρ - σ, 0)
mobρ(ρ, σ) = max(1 - ρ - 0.5σ, 0)
mobσ(ρ, σ) = max(1 - σ - 0.5ρ, 0)

end # module

ext: jl


text:
module FluidInducedAseismicSlip

using FastGaussQuadrature
using SpecialFunctions
using LinearAlgebra

export injection_analytical_gs
export lambda_analytical_gs
export slip_distribution_gs

"""
    injection_analytical_gs(T::Float64, N::Int64 = 100)

Solves the slip distribution using Gauss-Chebyshev quadratures due to an
analytical pressure expression.
Inputs:
    T: the stress parameter
    N: number of quadrature points
Outputs:
    x: the Gauss-Checbyshev quadrature points
    δ: the slip distribution
    λ: the slip to fluid migration factor
"""
function injection_analytical_gs(T::Float64, N::Int64=100)
    # Solve for λ
    λ = lambda_analytical_gs(T, N)

    # Slip
    (x, δ) = slip_distribution_gs(λ, N)

    return (x, δ, λ)
end

"""
    lambda_analytical_gs(
        T::Float64,
        N::Int64,
        max_iters::Int64 = 200,
        abs_tol::Float64 = 1.0e-10,
        debug::Bool = False,
    )

Use Newton-Raphson iterations to solve for λ using Gauss-Chebyshev quadratures.
Inputs:
    T: the stress parameter
    N: number of quadrature points
    max_iters: the maximum number of Newton-Raphson iterations
    abs_tol: the absolute tolerance for the Newton solve
    debug: print convergence information
Outputs:
    λ: the slip to fluid migration factor
"""
function lambda_analytical_gs(
    T::Float64,
    N::Int64,
    max_iters::Int64=200,
    abs_tol::Float64=1.0e-10,
    debug::Bool=false,
)
    # Initialization
    λ = 1.0
    k = 0
    # Create Gauss-Chebyshev quadrature points
    s, w = gausschebyshev(N, 1)

    while (k < max_iters)
        # Update residuals and jacobian
        Res = residual_analytical(T, s, w, λ)
        Jac = jacobian_analytical(T, s, w, λ)

        # Check convergence
        if (abs(Res) <= abs_tol)
            if (debug)
                println("Solve converged! λ = ", λ, " after ", k, " iterations.")
            end
            return λ
        end

        # Update λ
        dλ = -Res / Jac
        λ += dλ
        k += 1
    end

    throw(ErrorException("Maximm number of iterations reached in the Newton solve!!!\n"))
end

"""
    residual_analytical(T::Float64, s::Vector{Float64}, w::Vector{Float64}, λ::Float64)

The residual for calculating λ
"""
function residual_analytical(T::Float64, s::Vector{Float64}, w::Vector{Float64}, λ::Float64)
    f(s) = erfc(λ * abs(s))
    return dot(w, f.(s)) / pi - T
end

"""
    jacobian_analytical(T::Float64, s::Vector{Float64}, w::Vector{Float64}, λ::Float64)

The jacobian for calculating λ
"""
function jacobian_analytical(T::Float64, s::Vector{Float64}, w::Vector{Float64}, λ::Float64)
    f(s) = -2 * abs(s) / sqrt(pi) * exp(-λ^2 * abs(s)^2)
    return dot(w, f.(s)) / pi
end

"""
    dF(s::Float64, u::Float64, λ::Float64)

The function dF from Viesca and Garagash (2018)
"""
function dF(s::Float64, u::Float64, λ::Float64)
    return 1  / π * erfc(λ * abs(s)) / (u - s)
end

"""
    F(u::Float64, λ::Float64, N::Int64)

The function F from Viesca and Garagash (2018)
"""
function F(u::Float64, λ::Float64, N::Int64)
    # Gauss-Chebyshev quadrature points
    s, w = gausschebyshev(N, 1)
    return dot(w, dF.(s, u, λ))
end

"""
    slip_gs!(δ::Vector{Float64}, x::Vector{Float64}, N::Int64, λ::Float64)

Slip evaluated with Gauss-Chebyshev quadrature
"""
function slip_gs!(δ::Vector{Float64}, x::Vector{Float64}, N::Int64, λ::Float64)
    # Gauss-Chebyshev quadrature points
    (s, w) = gausschebyshev(N-1, 2)

    # Slip weigth from Viesca and Garagash (2018)
    Φ = [0.5 * (sin(k * acos(xi)) / k - sin((k + 2) * acos(xi)) / (k + 2)) for xi in x, k = N:-1:1]
    B = [2 * sin(acos(sj)) * sin((k + 1) * acos(sj)) / (N + 1) for k = N:-1:1, sj in s]
    S = zeros(N, N-1)
    mul!(S, Φ, B)
    mul!(δ, S, F.(s, λ, N))
end

"""
    theta(x::Float64)

Theta function from Viesca and Garagash (2018)
"""
function theta(x::Float64)
    return acos(x::Float64)
end

"""
    slip_distribution_gs(λ::Float64, N::Int64)

Computes the slip distribution based on Gauss-Chebyshev quadratures
Inputs:
    λ: the slip to fluid migration factor
    N: number of quadrature points
Outputs
    x: the spatial discretization
    δ: the slip distribution

"""
function slip_distribution_gs(λ::Float64, N::Int64)
    x = collect(range(-1.0, 1.0, length=N))

    δ = similar(x)
    slip_gs!(δ, x, N, λ)
    return (x, δ)
end

end

max_stars_repo_name: ajacquey/FluidInducedAseismicSlip

max_issues_repo_name: ajacquey/FluidInducedAseismicSlip

max_forks_repo_name: ajacquey/FluidInducedAseismicSlip


text:
function hnorm(A::HMat,p="fro")
    if p == "fro"
        return hnormfro(A);
    elseif p == Inf
        error("Not been implemented");
    elseif p == 1
        error("Not been implemented");
    elseif p == 2
        error("Not been implemented");
    else
        error("hnorm only support p = 1,2,Inf");
    end
end

function hnormfro(A::HMat)
    if A.blockType == LOWRANK
        return sqrt(trace(A.UMat'*A.UMat*(A.VMat'*A.VMat)));
    elseif A.blockType == DENSE
        return vecnorm(A.DMat);
    elseif A.blockType == HMAT
        fro = 0;
        for i = 1:4
            for j = 1:4
                fro += hnormfro(A.childHMat[i,j])^2;
            end
        end
        return sqrt(fro);
    end
end

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
module FinSets
export FinOrd, FinOrdFunction, force, terminal, product, equalizer, pullback,
  initial, coproduct, coequalizer, pushout

using AutoHashEquals
using DataStructures: IntDisjointSets, union!, find_root

using ...GAT
using ...Theories: Category
using ..ShapeDiagrams
import ...Theories: dom, codom, id, compose, ⋅, ∘

# Data types
############

""" Finite ordinal (natural number).

An object in the category of finite ordinals, which is the skeleton of the
category of finite sets.
"""
@auto_hash_equals struct FinOrd{T<:Integer}
  n::T
end

""" Function between sets in the form of finite ordinals.

A morphism in the category of finite ordinals, which is the skeleton of the
category of finite sets.

In this data structure, the field `func` representing the function can have any
Julia type `T`, provided that `FinOrdFunction{T}` is callable. Usually, this
object will be an ordinary Julia function or an abstract vector. In the latter
case, the function (1↦1, 2↦3, 3↦2, 4↦3) is represented by the vector [1,3,2,3].
"""
@auto_hash_equals struct FinOrdFunction{F,T<:Integer}
  func::F
  dom::T
  codom::T
end
FinOrdFunction(f::AbstractVector) = FinOrdFunction(f, maximum(f))
FinOrdFunction(f::AbstractVector, codom::Integer) =
  FinOrdFunction(f, length(f), codom)

# Function objects are callable.
(f::FinOrdFunction)(i::Integer) = f.func(i)
(f::FinOrdFunction{<:AbstractVector})(i::Integer) = f.func[i]
(f::FinOrdFunction{typeof(id)})(i::Integer) = i

""" Force evaluation of function, yielding the vector representation.
"""
force(f::FinOrdFunction) = FinOrdFunction(map(f, 1:f.dom), f.codom)
force(f::FinOrdFunction{<:AbstractVector}) = f

# Category of finite ordinals
#############################

@instance Category(FinOrd, FinOrdFunction) begin
  dom(f::FinOrdFunction) = FinOrd(f.dom)
  codom(f::FinOrdFunction) = FinOrd(f.codom)

  id(A::FinOrd) = FinOrdFunction(id, A.n, A.n)

  function compose(f::FinOrdFunction, g::FinOrdFunction)
    @assert f.codom == g.dom
    FinOrdFunction(compose_functions(f.func, g.func), f.dom, g.codom)
  end
end

compose_functions(f,g) = as_function(g) ∘ as_function(f)
compose_functions(::typeof(id), g) = g
compose_functions(f, ::typeof(id)) = f
compose_functions(::typeof(id), ::typeof(id)) = id
compose_functions(f::AbstractVector, g::AbstractVector) = g[f]

as_function(f) = f
as_function(f::AbstractVector) = x -> f[x]

# Limits
########

terminal(::Type{FinOrd}) = FinOrd(1)
terminal(::Type{FinOrd{T}}) where T = FinOrd(one(T))

function product(A::FinOrd, B::FinOrd)
  m, n = A.n, B.n
  indices = CartesianIndices((m, n))
  π1 = FinOrdFunction(i -> indices[i][1], m*n, m)
  π2 = FinOrdFunction(i -> indices[i][2], m*n, n)
  Span(π1, π2)
end

function equalizer(f::FinOrdFunction, g::FinOrdFunction)
  @assert dom(f) == dom(g) && codom(f) == codom(g)
  m = f.dom
  FinOrdFunction(filter(i -> f(i) == g(i), 1:m), m)
end

""" Pullback of cospan of functions between finite ordinals.

TODO: This logic is completely generic. Make it independent of FinOrd.
"""
function pullback(cospan::Cospan{<:FinOrdFunction,<:FinOrdFunction})
  f, g = left(cospan), right(cospan)
  prod = product(dom(f), dom(g))
  π1, π2 = left(prod), right(prod)
  eq = equalizer(π1⋅f, π2⋅g)
  Span(eq⋅π1, eq⋅π2)
end

# Colimits
##########

initial(::Type{FinOrd}) = FinOrd(0)
initial(::Type{FinOrd{T}}) where T = FinOrd(zero(T))

function coproduct(A::FinOrd, B::FinOrd)
  m, n = A.n, B.n
  ι1 = FinOrdFunction(1:m, m, m+n)
  ι2 = FinOrdFunction(m+1:m+n, n, m+n)
  Cospan(ι1, ι2)
end

function coequalizer(f::FinOrdFunction, g::FinOrdFunction)
  @assert dom(f) == dom(g) && codom(f) == codom(g)
  m, n = f.dom, f.codom
  sets = IntDisjointSets(n)
  for i in 1:m
    union!(sets, f(i), g(i))
  end
  h = [ find_root(sets, i) for i in 1:n ]
  roots = unique!(sort(h))
  FinOrdFunction([ searchsortedfirst(roots, r) for r in h], length(roots))
end

""" Pushout of span of functions between finite ordinals.

TODO: This logic is completely generic. Make it independent of FinOrd.
"""
function pushout(span::Span{<:FinOrdFunction,<:FinOrdFunction})
  f, g = left(span), right(span)
  coprod = coproduct(codom(f), codom(g))
  ι1, ι2 = left(coprod), right(coprod)
  coeq = coequalizer(f⋅ι1, g⋅ι2)
  Cospan(ι1⋅coeq, ι2⋅coeq)
end

end


text:
# Browser Window (Borrowed from Blink.jl)
@static if is_apple()
    launch(x) = run(`open $x`)
elseif is_linux()
    launch(x) = run(`xdg-open $x`)
elseif is_windows()
    launch(x) = run(`cmd /C start $x`)
end

# const connections = Dict{Int,WebSocket}() # WebSocket.id => WebSocket
const conditions = Dict{String,Condition}()
conditions["connected"] = Condition()
conditions["unloaded"] = Condition()

Endpoint("/PagesJL.js") do request::Request, session::Session
    d = Dict("session_id" => session.id,"host" => request.headers["Host"])
    template = Mustache.template_from_file(joinpath(dirname(@__FILE__),"..","res","PagesJL.js"))
    render(template,d)
end

ws = WebSocketHandler() do request::Request, client::WebSocket
    while true
        msg = JSON.parse(String(read(client)))
        if haskey(sessions,msg["session_id"])
            session = sessions[msg["session_id"]]
            session.client = client
            haskey(msg,"args") ? callbacks[msg["name"]].callback(msg["args"]) : callbacks[msg["name"]].callback()
        end
    end
end

http = HttpHandler() do request::Request, response::Response
    route = URI(request.resource).path
    if haskey(pages,route)
        if isequal(route,"/PagesJL.js")
            if haskey(request.headers,"Referer")
                referer = URI(request.headers["Referer"]).path
                session = Session(referer)
            else
                session = Session()
            end
            res = Response(pages[route].handler(request,session))
        else
            res = Response(pages[route].handler(request))
        end
    else
        res = "Page not found."
    end
    res
end

server = Server(http,ws)

function start(p = 8000)
    global port = p
    @async run(server, port)
end


text:
@testset "Add tests for mixture of Gaussian kernels" begin

        M = Mixture(3)

        M.w .= [0.184018754243173;
                0.195775850480659;
                0.071604608234824;
                0.176323993535265;
                0.169941211284233;
                0.042183665772717;
                0.030869595116112;
                0.129282321333016]

        X = reshape([-0.5 -0.2; 0.1 -0.8; 2.0 1.1],(3,2))

        @test norm(log_pdf_mixture(M, X) -  [-20.597416443372648;
                                     -21.03728674466139])<1e-8

end


text:
# All types in Vega.jl are subtypes of AbstractVegaType. Overloading `call` as
# a parametric type over the field of `AbstractVegaType`s allows us to get much
# better performance for the first visualization after importing Vega
# abstract AbstractVegaType
#
# function Base.call{T<:AbstractVegaType}(::Type{T}; kwargs...)
#     out = T()
#     for (sym, val) in kwargs
#         setfield!(out, sym, val)
#     end
#     out
# end

# tojs(s::Symbol) = string(s)
# tojs(v::Vector) = [tojs(v_i) for v_i in v]
#
# function tojs(v::Dict)
#     res = Dict()
#     for (k, v) in v
#         res[k] = tojs(v)
#     end
#     return res
# end
# tojs(x::Any) = x

#Takes in a type name and its specification
#Returns Expr for the type code as if you had written it yourself

#Ex: padding_spec = [(:top, Number, 80, false), (:left, Number, 80, false), (:bottom, Number, 80, false), (:right, Number, 80, false)]
# > maketype(:VegaPadding, padding_spec)
# :(type VegaPadding <: AbstractVegaType
#     top::Number
#     left::Number
#     bottom::Number
#     right::Number
# end)

function maketype(typename::Symbol, spec)
    n = length(spec)
    lines = Array(Expr, n)
    for i in 1:n
        entry = spec[i]
        lines[i] = Expr(:(::), entry[1], Union{entry[2], Void}) #Was fielddef(entry::Tuple) function
    end

    return Expr(:type, #type
                true, Expr(:<:, typename, AbstractVegaType), Expr(:block, lines...) #arguments
                )
end

#Takes in typename and its specification
#Builds a function with the same name as the type (i.e function VegaPadding and type :VegaPadding)
#An Expr is returned, as if you had typed it yourself

#:(function VegaPadding()
#        VegaPadding(80,80,80,80)
#    end)

function makedefaultfunc(typename::Symbol, spec)
    return Expr(:function, #type
                Expr(:call, typename),
                Expr(:block,
                     Expr(:call, typename, map(entry -> entry[3], spec)...)))
end

#Intermediate helper function to maketojs()
#Make a dict from fieldname
#If fieldname happens to be _type (because type is a reserved keyword), take the underscore off when placing in dict
function makehash(fieldname::Symbol)
    if fieldname == :_type
        return Expr(:(=),
                    Expr(:ref, :res, "type"),
                    Expr(:call, :tojs, Expr(:., :x, Expr(:quote, fieldname))))
    else
        return Expr(:(=),
                    Expr(:ref, :res, string(fieldname)),
                    Expr(:call, :tojs, Expr(:., :x, Expr(:quote, fieldname))))
    end
end

##Intermediate helper function to maketojs()
function makejsline(entry::Tuple)
    fieldname = entry[1]

    return Expr(:if,
                Expr(:comparison,
                    Expr(:., :x, Expr(:quote, fieldname)),
                    :!=,
                    :nothing),
                Expr(:block, makehash(fieldname)))

end

#Intermediate helper function to maketojs()
function makejsbody(spec)
    return Expr(:block,
                :(res = Dict()),
                map(makejsline, spec)...,
                :(return res))
end

#Function uses makejsbody, makejsline, makehash to create a custom tojs() function for each type

#Ex:
# > maketojs(:VegaPadding, padding_spec)
# :(function tojs(x::VegaPadding)
#         res = Dict()
#         res["top"] = tojs(x.top)
#         res["left"] = tojs(x.left)
#         res["bottom"] = tojs(x.bottom)
#         res["right"] = tojs(x.right)
#         return res
#     end)
function maketojs(typename::Symbol, spec)
    return Expr(:function,
                Expr(:call,
                     :tojs,
                     Expr(:(::), :x, typename)),
                     makejsbody(spec))
end

#Create primitives from one function instead of repeating
#Use of eval() takes all Expr created above, evaluates their results, so that they become Julia citizens
function primitivefactory(create::Symbol, spec::AbstractArray)
    #println(maketype(create, spec))
    #eval(maketype(create, spec))
    #println(makedefaultfunc(create, spec))
    #eval(makedefaultfunc(create, spec))
    #println(maketojs(create, spec))
    #eval(maketojs(create, spec))
    #println("")
end

max_stars_repo_name: jpthor/vega


text:
#!/usr/bin/env julia

include("./lib/jlBarraCUDA.jl")
using .jlBarraCUDA

x = cuDeviceCount()
println(x)

max_stars_repo_name: MarcelDykstra/barracuda

max_issues_repo_name: MarcelDykstra/BarraCUDA

max_forks_repo_name: MarcelDykstra/BarraCUDA


text:
export PhasePartition
# Thermodynamic states
export ThermodynamicState,
       PhaseDry,
       PhaseEquil,
       PhaseNonEquil,
       TemperatureSHumEquil,
       LiquidIcePotTempSHumEquil,
       LiquidIcePotTempSHumEquil_given_pressure,
       LiquidIcePotTempSHumNonEquil,
       LiquidIcePotTempSHumNonEquil_given_pressure

"""
    PhasePartition

Represents the mass fractions of the moist air mixture.

# Constructors

    PhasePartition(q_tot::Real[, q_liq::Real[, q_ice::Real]])
    PhasePartition(ts::ThermodynamicState)

See also [`PhasePartition_equil`](@ref)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhasePartition{FT<:Real}
  "total specific humidity"
  tot::FT
  "liquid water specific humidity (default: `0`)"
  liq::FT
  "ice specific humidity (default: `0`)"
  ice::FT
end

PhasePartition(q_tot::FT,q_liq::FT) where {FT<:Real} =
  PhasePartition(q_tot, q_liq, zero(FT))
PhasePartition(q_tot::FT) where {FT<:Real} =
  PhasePartition(q_tot, zero(FT), zero(FT))

"""
    ThermodynamicState{FT}

A thermodynamic state, which can be initialized for
various thermodynamic formulations (via its sub-types).
All `ThermodynamicState`'s have access to functions to
compute all other thermodynamic properties.
"""
abstract type ThermodynamicState{FT} end

"""
    PhaseEquil{FT} <: ThermodynamicState

A thermodynamic state assuming thermodynamic equilibrium (therefore, saturation adjustment
may be needed).

# Constructors

    PhaseEquil(e_int, ρ, q_tot)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhaseEquil{FT} <: ThermodynamicState{FT}
  "internal energy"
  e_int::FT
  "density of air (potentially moist)"
  ρ::FT
  "total specific humidity"
  q_tot::FT
  "temperature: computed via [`saturation_adjustment`](@ref)"
  T::FT
end
function PhaseEquil(e_int::FT,
                    ρ::FT,
                    q_tot::FT,
                    maxiter::Int=3,
                    tol::FT=FT(1e-1),
                    sat_adjust::F=saturation_adjustment
                    ) where {FT<:Real,F}
    return PhaseEquil{FT}(e_int, ρ, q_tot, sat_adjust(e_int, ρ, q_tot, tol, maxiter))
end

"""
    PhaseDry{FT} <: ThermodynamicState

A dry thermodynamic state (`q_tot = 0`).

# Constructors

    PhaseDry(e_int, ρ)

# Fields

$(DocStringExtensions.FIELDS)
"""
struct PhaseDry{FT} <: ThermodynamicState{FT}
  "internal energy"
  e_int::FT
  "density of dry air"
  ρ::FT
end

"""
    LiquidIcePotTempSHumEquil(θ_liq_ice, ρ, q_tot)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from:

 - `θ_liq_ice` liquid-ice potential temperature
 - `ρ` (moist-)air density
 - `q_tot` total specific humidity
 - `tol` tolerance for saturation adjustment
 - `maxiter` maximum iterations for saturation adjustment
"""
function LiquidIcePotTempSHumEquil(θ_liq_ice::FT,
                                   ρ::FT,
                                   q_tot::FT,
                                   maxiter::Int=30,
                                   tol::FT=FT(1e-1)
                                   ) where {FT<:Real}
    T = saturation_adjustment_q_tot_θ_liq_ice(θ_liq_ice, ρ, q_tot, tol, maxiter)
    q_pt = PhasePartition_equil(T, ρ, q_tot)
    e_int = internal_energy(T, q_pt)
    return PhaseEquil{FT}(e_int, ρ, q_tot, T)
end

"""
    LiquidIcePotTempSHumEquil_given_pressure(θ_liq_ice, p, q_tot)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from:

 - `θ_liq_ice` liquid-ice potential temperature
 - `p` pressure
 - `q_tot` total specific humidity
 - `tol` tolerance for saturation adjustment
 - `maxiter` maximum iterations for saturation adjustment
"""
function LiquidIcePotTempSHumEquil_given_pressure(θ_liq_ice::FT,
                                                  p::FT,
                                                  q_tot::FT,
                                                  maxiter::Int=30,
                                                  tol::FT=FT(1e-1)
                                                  ) where {FT<:Real}
    T = saturation_adjustment_q_tot_θ_liq_ice_given_pressure(θ_liq_ice, p, q_tot, tol, maxiter)
    ρ = air_density(T, p, PhasePartition(q_tot))
    q = PhasePartition_equil(T, ρ, q_tot)
    e_int = internal_energy(T, q)
    return PhaseEquil{FT}(e_int, ρ, q.tot, T)
end

"""
    TemperatureSHumEquil(T, p, q_tot)

Constructs a [`PhaseEquil`](@ref) thermodynamic state from temperature.

 - `T` temperature
 - `p` pressure
 - `q_tot` total specific humidity
"""
function TemperatureSHumEquil(T::FT, p::FT, q_tot::FT) where {FT<:Real}
    ρ = air_density(T, p, PhasePartition(q_tot))
    q = PhasePartition_equil(T, ρ, q_tot)
    e_int = internal_energy(T, q)
    return PhaseEquil{FT}(e_int, ρ, q_tot, T)
end

"""
   	 PhaseNonEquil{FT} <: ThermodynamicState

A thermodynamic state assuming thermodynamic non-equilibrium (therefore, temperature can
be computed directly).

# Constructors

    PhaseNonEquil(e_int, q::PhasePartition, ρ)

# Fields

$(DocStringExtensions.FIELDS)

"""
struct PhaseNonEquil{FT} <: ThermodynamicState{FT}
  "internal energy"
  e_int::FT
  "density of air (potentially moist)"
  ρ::FT
  "phase partition"
  q::PhasePartition{FT}
end

"""
    LiquidIcePotTempSHumNonEquil(θ_liq_ice, ρ, q_pt)

Constructs a [`PhaseNonEquil`](@ref) thermodynamic state from:

 - `θ_liq_ice` liquid-ice potential temperature
 - `ρ` (moist-)air density
 - `q_pt` phase partition
and, optionally
 - `tol` tolerance for non-linear equation solve
 - `maxiter` maximum iterations for non-linear equation solve
"""
function LiquidIcePotTempSHumNonEquil(θ_liq_ice::FT,
                                      ρ::FT,
                                      q_pt::PhasePartition{FT},
                                      maxiter::Int=5,
                                      tol::FT=FT(1e-1)
                                      ) where {FT<:Real}
    T = air_temperature_from_liquid_ice_pottemp_non_linear(θ_liq_ice, ρ, tol, maxiter, q_pt)
    e_int = internal_energy(T, q_pt)
    return PhaseNonEquil{FT}(e_int, ρ, q_pt)
end

"""
    LiquidIcePotTempSHumNonEquil_given_pressure(θ_liq_ice, p, q_pt)

Constructs a [`PhaseNonEquil`](@ref) thermodynamic state from:

 - `θ_liq_ice` liquid-ice potential temperature
 - `p` pressure
 - `q_pt` phase partition
"""
function LiquidIcePotTempSHumNonEquil_given_pressure(θ_liq_ice::FT, p::FT, q_pt::PhasePartition{FT}) where {FT<:Real}
    T = air_temperature_from_liquid_ice_pottemp_given_pressure(θ_liq_ice, p, q_pt)
    ρ = air_density(T, p, q_pt)
    e_int = internal_energy(T, q_pt)
    return PhaseNonEquil{FT}(e_int, ρ, q_pt)
end

"""
    fixed_lapse_rate_ref_state(z::FT,
                               T_surface::FT,
                               T_min::FT) where {FT<:AbstractFloat}

Fixed lapse rate hydrostatic reference state
"""
function fixed_lapse_rate_ref_state(z::FT,
                                    T_surface::FT,
                                    T_min::FT) where {FT<:AbstractFloat}
  Γ = FT(grav)/FT(cp_d)
  z_tropopause = (T_surface - T_min) / Γ
  H_min = FT(R_d) * T_min / FT(grav)
  T = max(T_surface - Γ*z, T_min)
  p = FT(MSLP)*(T / T_surface)^(FT(grav)/(FT(R_d)*Γ))
  T == T_min && (p = p * exp(-(z-z_tropopause)/FT(H_min)))
  ρ = p / (FT(R_d) * T)
  return T,p,ρ
end

"""
    tested_convergence_range(FT, n)

A range of input arguments to thermodynamic state constructors
 - `e_int` internal energy
 - `ρ` (moist-)air density
 - `q_tot` total specific humidity
 - `q_pt` phase partition
 - `T` air temperature
 - `θ_liq_ice` liquid-ice potential temperature
that are tested for convergence in saturation adjustment.

Note that the output vectors are of size ``n*n_RH``, and they
should span the input arguments to all of the constructors.
"""
function tested_convergence_range(FT, n)
  n_RS1 = 10
  n_RS2 = 20
  n_RS = n_RS1+n_RS2
  z_range  = range(FT(0), stop=FT(2.5e4), length=n)
  relative_sat1 = range(FT(0), stop=FT(1), length=n_RS1)
  relative_sat2 = range(FT(1), stop=FT(1.02), length=n_RS2)
  relative_sat = [relative_sat1...,relative_sat2...]
  T_min = FT(150)
  T_surface = FT(350)

  args = fixed_lapse_rate_ref_state.(z_range, Ref(T_surface), Ref(T_min))
  T,p,ρ = getindex.(args, 1),
          getindex.(args, 2),
          getindex.(args, 3)

  p  = collect(Iterators.flatten([p  for RS in 1:n_RS]))
  ρ  = collect(Iterators.flatten([ρ  for RS in 1:n_RS]))
  T  = collect(Iterators.flatten([T  for RS in 1:n_RS]))
  relative_sat = collect(Iterators.flatten([relative_sat for RS in 1:n]))

  # Additional variables
  q_sat = q_vap_saturation.(T, ρ)
  q_tot = min.(relative_sat .*q_sat, FT(1))
  q_pt = PhasePartition_equil.(T, ρ, q_tot)
  e_int = internal_energy.(T, q_pt)
  θ_liq_ice = liquid_ice_pottemp.(T, ρ, q_pt)
  return e_int, ρ, q_tot, q_pt, T, p, θ_liq_ice
end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
export SpinOneSite,
       spinOneSites

struct SpinOneSite <: Site
    s::Index
    SpinOneSite(I::Index) = new(I)
end
SpinOneSite(n::Int) = SpinOneSite(Index(3, "Site,S=1,n=$n"))

function spinOneSites(N::Int;kwargs...)::SiteSet
  sites = SiteSet(N)
  for n=1:N
    set(sites,n,SpinOneSite(n))
  end
  return sites
end

function operator(site::SpinOneSite,
                  opname::AbstractString)::ITensor
    s = site.s
    sP = prime(site.s)
    Up = s(1)
    UpP = sP(1)
    Z0 = s(2)
    Z0P = sP(2)
    Dn = s(3)
    DnP = sP(3)

    Op = ITensor(dag(s), s')

    if opname == "S⁺" || opname == "Splus" || opname == "S+"
        Op[Dn, Z0P] = √2
        Op[Z0, UpP] = √2
    elseif opname == "S⁻" || opname == "Sminus" || opname == "S-"
        Op[Up, Z0P] = √2
        Op[Z0, DnP] = √2
    elseif opname == "Sˣ" || opname == "Sx"
        Op[Up, Z0P] = 1.0/√2
        Op[Z0, UpP] = 1.0/√2
        Op[Z0, DnP] = 1.0/√2
        Op[Dn, Z0P] = 1.0/√2
    elseif opname == "iSʸ" || opname == "iSy"
        Op[Up, Z0P] = -1.0/√2
        Op[Z0, UpP] = +1.0/√2
        Op[Z0, DnP] = -1.0/√2
        Op[Dn, Z0P] = +1.0/√2
    elseif opname == "Sʸ" || opname == "Sy"
        Op = complex(Op)
        Op[Up, Z0P] = +1.0/√2im
        Op[Z0, UpP] = -1.0/√2im
        Op[Z0, DnP] = +1.0/√2im
        Op[Dn, Z0P] = -1.0/√2im
    elseif opname == "Sᶻ" || opname == "Sz"
        Op[Up, UpP] = 1.0
        Op[Dn, DnP] = -1.0
    elseif opname == "Sᶻ²" || opname == "Sz2"
        Op[Up, UpP] = 1.0
        Op[Dn, DnP] = 1.0
    elseif opname == "Sˣ²" || opname == "Sx2"
        Op[Up, UpP] = 0.5
        Op[Up, DnP] = 0.5
        Op[Z0, Z0P] = 1.0
        Op[Dn, UpP] = 0.5
        Op[Dn, DnP] = 0.5
    elseif opname == "Sʸ²" || opname == "Sy2"
        Op[Up, UpP] = 0.5
        Op[Up, DnP] = -0.5
        Op[Z0, Z0P] = 1.0
        Op[Dn, UpP] = -0.5
        Op[Dn, DnP] = 0.5
    elseif opname == "projUp"
        Op[Up, UpP] = 1.
    elseif opname == "projZ0"
        Op[Z0, Z0P] = 1.
    elseif opname == "projDn"
        Op[Dn, DnP] = 1.
    elseif opname == "XUp"
        xup = ITensor(ComplexF64,s)
        xup[Up] = 0.5
        xup[Z0] = im*√2
        xup[Dn] = 0.5
        return xup
    elseif opname == "XZ0"
        xZ0 = ITensor(ComplexF64,s)
        xZ0[Up] = im*√2
        xZ0[Dn] = -im*√2
        return xZ0
    elseif opname == "XDn"
        xdn = ITensor(ComplexF64,s)
        xdn[Up] = 0.5
        xdn[Z0] = -im*√2
        xdn[Dn] = 0.5
        return xdn
    else
        throw(ArgumentError("Operator name '$opname' not recognized for SpinOneSite"))
    end
    return Op
end


text:
const CPetscViewer = Ptr{Cvoid}

"""
    AbstractViewer{PetscLib <: PetscLibType}

Abstract type of PETSc viewer.

Manual: [`PetscViewer`](https://petsc.org/release/docs/manualpages/Viewer/PetscViewer.html)
"""
abstract type AbstractViewer{PetscLib <: PetscLibType} end

Base.cconvert(::Type{CPetscViewer}, obj::AbstractViewer) = obj.ptr

Base.unsafe_convert(::Type{Ptr{CPetscViewer}}, obj::AbstractViewer) =
    convert(Ptr{CPetscViewer}, pointer_from_objref(obj))

"""
    ViewerStdout(petsclib, comm = MPI.COMM_SELF)

Create an ASCII `PetscViewer` for the `comm`

Manual: [`PETSC_VIEWER_STDOUT_`](https://petsc.org/release/docs/manualpages/Viewer/PETSC_VIEWER_STDOUT_.html)
"""
mutable struct ViewerStdout{PetscLib} <: AbstractViewer{PetscLib}
    ptr::CPetscViewer
    comm::MPI.Comm
end

@for_libpetsc function ViewerStdout(
    ::$UnionPetscLib,
    comm::MPI.Comm = MPI.COMM_SELF,
)
    ptr = ccall(
        (:PETSC_VIEWER_STDOUT_, $petsc_library),
        CPetscViewer,
        (MPI.MPI_Comm,),
        comm,
    )
    return ViewerStdout{$PetscLib}(ptr, comm)
end
@for_libpetsc function Base.push!(
    viewer::AbstractViewer{$PetscLib},
    format::PetscViewerFormat,
)
    @chk ccall(
        (:PetscViewerPushFormat, $petsc_library),
        PetscErrorCode,
        (CPetscViewer, PetscViewerFormat),
        viewer,
        format,
    )
    return nothing
end
@for_libpetsc function Base.pop!(viewer::AbstractViewer{$PetscLib})
    @chk ccall(
        (:PetscViewerPopFormat, $petsc_library),
        PetscErrorCode,
        (CPetscViewer,),
        viewer,
    )
    return nothing
end

function with(f, viewer::AbstractViewer, format::PetscViewerFormat)
    push!(viewer, format)
    try
        f()
    finally
        pop!(viewer)
    end
end

# ideally we would capture the output directly, but this looks difficult
# easiest option is to redirect stdout
# based on suggestion from https://github.com/JuliaLang/julia/issues/32567
function _show(io::IO, obj)
    old_stdout = stdout
    try
        rd, = redirect_stdout()
        view(obj)
        Libc.flush_cstdio()
        flush(stdout)
        write(io, readavailable(rd))
    finally
        redirect_stdout(old_stdout)
    end
    return nothing
end

#=
# PETSc_jll isn't built with X support
mutable struct ViewerDraw <: AbstractViewer
    ptr::CPetscViewer
    comm::MPI.Comm
end
function ViewerDraw(comm::MPI.Comm)
    ptr = ccall((:PETSC_VIEWER_DRAW_, libpetsc), CPetscViewer, (MPI.MPI_Comm,), comm)
    return ViewerDraw(ptr, comm)
end
=#


text:
module TestMaybeFunction

using Maybe.Extras
using Test

retargs(args...; kwargs...) = (args, values(kwargs))

@testset begin
    @test maybe(retargs)(nothing) === nothing
    @test maybe(retargs)(1, nothing, 2) === nothing
    @test maybe(retargs)(1, nothing, 2; a = 3) === nothing
    @test maybe(retargs)(Some(nothing)) === Some(retargs(nothing))
    @test maybe(retargs)(Some(nothing); a = 1) === Some(retargs(nothing, a = 1))
    @test maybe(retargs)(Some(nothing); a = Some(nothing)) ===
          Some(retargs(nothing, a = nothing))
end

end  # module


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module Simulate_FreeBodyRotation  (examples/withoutMacros)

Simulate a free body rotation with quaternions (version that uses Variables, but no macros)

A body can rotate freely around the origin and a sin-wave torque vector is applied on it.
The system is mathematically described as a fully-implicit, index-1 DAE model of a body
with unconstrained rotational motion in 3D using quaternions.

A simulation tool of a corresponding Modelica model that transforms the model to an ODE
would need to use dynamic state selection and solve nonlinear algebraic equations
within the ODE. Due to the implicit index-1 DAE description, this is much simpler here.

# DAE equations
```
  t := A.*sin(2*pi*freqHz*time + phase);
  0  = w - 2*([q[4],  q[3], -q[2], -q[1];
              -q[3],  q[4],  q[1], -q[2];
             q[2], -q[1],  q[4], -q[3]])*der(q);
  0  = I*der(w) + cross(w, I*w) - t;
  0  = dot(q,q)-1.0;

Initialization:
  q = [0.1, 0.5, 0.0, 1.0]   # inconsistent values (fixed during initialization)
  w = zeros(3)

Reference Modelica model:
   ModelicaReferenceModels.ODAEs.FreeBody
```
"""
module Simulate_FreeBodyRotation

using ModiaMath
using ModiaMath.LinearAlgebra  # included via ModiaMath, to avoid requirement to add it in the standard environment
using ModiaMath.StaticArrays   # included via ModiaMath, to avoid requirement to add it in the standard environment

#            q[1] = 0.1 changed to 0.08908708063747484
#            q[2] = 0.5 changed to 0.445435403187374
#            q[4] = 1.0 changed to 0.8908708063747479
#                             q0     = [0.1, 0.5, 0.0, 1.0],
# q0     = [0.08908708063747484, 0.445435403187374, 0.0, 0.8908708063747479],

mutable struct FreeBodyRotationWithoutMacro <: ModiaMath.AbstractComponentWithVariables
    _internal::ModiaMath.ComponentInternal

    # Parameters
    I::SMatrix{3,3,Float64,9}
    A::SVector{3,Float64}
    freqHz::SVector{3,Float64}
    phase::SVector{3,Float64}
    q0::SVector{4,Float64}
    w0::SVector{3,Float64}

    # Variables
    q::RealSVector{4}
    derq::RealSVector{4}
    w::RealSVector3
    derw::RealSVector3
    tau::RealSVector3
    residue_w::RealSVector3
    residue_t::RealSVector3
    residue_q::RealScalar

    function FreeBodyRotationWithoutMacro(;I=Diagonal([1.0,2.0,3.0]),
                                         A=[3.0,4.0,5.0],
                                         freqHz=[0.3,0.2,0.1],
                                         phase=[0,0.5235987755983,1.0471975511966],
                                         q0=[0.1, 0.5, 0.0, 1.0],
                                         w0=zeros(3))
        this = new(ModiaMath.ComponentInternal(:FreeBodyRotationWithoutMacro, nothing), I, A, freqHz, phase, q0, w0)

        @assert(size(I) == (3, 3))
        @assert(isposdef(I))
        @assert(length(A) == 3)
        @assert(length(freqHz) == 3)
        @assert(length(phase) == 3)
        @assert(minimum(A) > 0.0)

        q         = RealSVector{4}(:q, this, numericType=ModiaMath.XD_IMP, info="Quaternions",                         start=q0, fixed=false)
        derq      = RealSVector{4}(:derq, this, numericType=ModiaMath.DER_XD_IMP, integral=q, info="der(q)",          unit="1/s")
        w         = RealSVector3(:w, this, numericType=ModiaMath.XD_IMP, info="Angular velocity",     unit="rad/s",  start=w0, fixed=true)
        derw      = RealSVector3(:derw, this, numericType=ModiaMath.DER_XD_IMP, integral=w, info="Angular acceleration", unit="rad/s^2")
        tau       = RealSVector3(:tau, this, numericType=ModiaMath.WR, info="Driving torque",       unit="N*m")

        residue_w = RealSVector3(:residue_w, this, numericType=ModiaMath.FD_IMP, info="Angular velocity residue")
        residue_t = RealSVector3(:residue_t, this, numericType=ModiaMath.FD_IMP, info="Angular momentum equation residue")
        residue_q = RealScalar(:residue_q, this, numericType=ModiaMath.FC, info="Quaternion constraint residue")

        return this
    end
end

function ModiaMath.computeVariables!(b::FreeBodyRotationWithoutMacro, sim::ModiaMath.SimulationState)
    time::Float64 = ModiaMath.getTime(sim)
    I             = b.I
    A             = b.A
    freqHz        = b.freqHz
    phase         = b.phase
    q             = b.q.value
    derq          = b.derq.value
    w             = b.w.value
    derw          = b.derw.value

    b.tau.value       = A .* sin.(2.0 * pi * freqHz * time + phase)
    b.residue_w.value = w - 2.0 * ( @SMatrix( [ q[4]  q[3] -q[2] -q[1];
                                            -q[3]  q[4]  q[1] -q[2];
                                             q[2] -q[1]  q[4] -q[3]] ) * derq)
    b.residue_t.value = I * derw + cross(w, I * w) - b.tau.value
    b.residue_q.value = dot(q, q) - 1.0

    return nothing
end

simulationModel = ModiaMath.SimulationModel(FreeBodyRotationWithoutMacro(), stopTime=5.0, tolerance=1e-8)
#ModiaMath.print_ModelVariables(simulationModel)
result          = ModiaMath.simulate!(simulationModel, log=true)
ModiaMath.plot(result, [:q, :w, :derw, :tau])

end


text:
import Cxx
import MuJoCo

import Base.getindex
import Base.setindex!

using Cxx
using MuJoCo

include("simenvi.jl")

##STATE

struct SimState
    mjStates::Vector{MJMatrix}
    nx::Int64

    function SimState()
        return new([], 0)
    end

    function SimState(data::MJData)
        l = [data.qpos, data.qvel, data.qacc, data.site_xpos]
        return new(l, data.nq + 2 * data.nv + data.nsite * 3)
    end

end

function getindex(state::SimState, i)
    targList = 1
    while length(state.mjStates[targList]) < i
        i -= length(state.mjStates[targList])
        targList += 1
    end

    return state.mjStates[targList].mat[i]
end

function getindex(state::SimState, l, i, j)
    return state.mjStates[l][i, j]
end

function setindex!(state::SimState, val, i)
    targList = 1
    while length(state.mjStates[targList]) < i
        i -= length(state.mjStates[targList])
        targList += 1
    end

    state.mjStates[targList].mat[i] = val
end

function setindex!(state::SimState, val, l, i, j)
    state.mjStates[l].mat[j, i] = val
end

##CONTROL

#struct AltControls
#    nu::Int64
#    ctrl::MJMatrix
#    normalized::Vector{mjtNum}
#    ctrl_magnitude::MJMatrix
#
#    function AltControls(data::MJData)
#        mat = data.qfrc_applied
#        norm = []
#        rng = MJMatrix(Matrix{mjtNum}(0, 0))
#
#        new(data.nv, mat, norm, rng)
#    end
#
#    function AltControls()
#        return new(0, MJMatrix(Matrix{mjtNum}(0, 0)), [], MJMatrix(Matrix{mjtNum}(0, 0)))
#    end
#end
#
#function SetSimControls(ctrl::AltControls, norm::Vector{mjtNum})
#    for i=1:ctrl.nu
#        ctrl.ctrl[i] = norm[i]
#    end
#
#    return ctrl
#end
#
#function getitem(ctrl::AltControls, i)
#    return ctrl.ctrl[i, 1]
#end
#
#function setitem!(ctrl::AltControls, val, i)
#    ctrl.ctrl.mat[1, i] = val
#end
#

struct AltControls
    nu::Int64
    ctrl::MJMatrix
    normalized::Vector{mjtNum}
    ctrl_magnitude::MJMatrix

    function AltControls(data::MJData)
        mat = data.qpos
        norm = []
        rng = MJMatrix(Matrix{mjtNum}(0, 0))

        new(data.nq, mat, norm, rng)
    end

    function AltControls()
        return new(0, MJMatrix(Matrix{mjtNum}(0, 0)), [], MJMatrix(Matrix{mjtNum}(0, 0)))
    end
end

function SetSimControls(ctrl::AltControls, norm)
    for i=1:ctrl.nu
        ctrl.ctrl[i] = norm[i]
    end

    return ctrl
end

function getitem(ctrl::AltControls, i)
    return ctrl.ctrl[i, 1]
end

function setitem!(ctrl::AltControls, val, i)
    ctrl.ctrl.mat[1, i] = val
end

##Accessories to the SimEnv

function SetSimState(sim::MJSimEnv, xi)
    for i=1:length(xi)
        sim.state[i] = xi[i]
    end
end

function GetSimState(sim::MJSimEnv)
    xi = Matrix{mjtNum}(sim.state.nx, 1)
    for i=1:sim.state.nx
        xi[i] = sim.state[i]
    end
    return xi
end

##CONTROLS

function SetSimControls(sim::MJSimEnv, ui)
    SetSimControls(sim.controls, ui)
    #SetControls(sim)
end

function GetSimControls(sim::MJSimEnv)
    u = zeros(sim.controls.nu, 1)
    for i=1:sim.controls.nu
        u[i] = sim.controls.ctrl[i, 1]
    end
    return u
end

###ACTUAL DYNAMICS FUNCTIONS

#transition
function mjDynamicsF(xi, ui, i, sim, intermediateCost)
    SetSimState(sim, xi)
    SetSimControls(sim, ui)
    mj_forward(sim.model, sim.data)
    xnew = GetSimState(sim)
    cost = intermediateCost(xnew, ui, sim)
    return (xnew, cost)
end

#compute jacobian and hessian along each point in the given trajectory
#time invariant, can drop I??
function mjDynamicsFx(x, u, I, sim, costx, costu, costxx, costuu, linear=false)
    #compute derivatives:
    nx = sim.state.nx
    #fx, fu, fxx, fuu, fxu, cx, cu, cxx, cxu, cuu
    fx = zeros(nx, nx, length(I) + 1)
    fu = zeros(nx, sim.controls.nu, length(I) + 1)
    cx = zeros(nx, length(I) + 1)
    cu = zeros(sim.controls.nu, length(I))
    cxx = zeros(nx, nx)
    cxu = zeros(nx, sim.controls.nu)
    cuu = zeros(sim.controls.nu, sim.controls.nu)

    fxx = fxu = fuu = []

    if !linear
        fxx = zeros(nx, nx, nx, length(I) + 1)
        fuu = zeros(nx, sim.controls.nu, sim.controls.nu, length(I) + 1)
        fxu = zeros(nx, nx, sim.controls.nu, length(I) + 1)
    end

    #for every point on trajectory
    for t in I
        #set states
        SetSimState(sim, x[:, t])
        SetSimControls(sim, u[:, t])
        mj_forward(sim.model, sim.data)

        dSim = MJSimEnv{SimState, AltControls}(sim.model, s -> (0), s -> (false), SimState(), AltControls())
        SetSimState(dSim, x[:, t])
        SetSimControls(dSim, u[:, t])
        mj_forward(dSim.model, dSim.data)
        for i=1:3
            mj_forward(dSim.model, dSim.data)
        end

        ddSim = MJSimEnv{SimState, AltControls}(sim.model, s -> (0), s -> (false), SimState(), AltControls())
        SetSimState(ddSim, x[:, t])
        SetSimControls(ddSim, u[:, t])
        mj_forward(ddSim.model, ddSim.data)
        for i=1:3
            mj_forward(ddSim.model, ddSim.data)
        end

        ###DYNAMICS MODEL DERIVATIVES these top level loops can be parallelized

        #fu & fuu
        for nu1=1:dSim.controls.nu

            dSim.controls.ctrl[nu1] += 1.0e-4

            #warmstart acc
            icxx"""
            mju_copy($(dSim.data.d)->qacc_warmstart, $(sim.data.d)->qacc_warmstart, $(sim.model.m)->nv);
            """

            #step to get difference
            for step=1:1
                mj_forward(dSim.model, dSim.data)
            end

            #calc fu
            for nx3=1:sim.state.nx
                fu[nx3, nu1, t] = (dSim.state[nx3] - sim.state[nx3]) / 2.0 / 1.0e-4
            end

            if !linear
                for nu2=1:ddSim.controls.nu

                    ddSim.controls.ctrl[nu2] -= 1.0e-4

                    #warmstart acc
                    icxx"""
                    mju_copy($(ddSim.data.d)->qacc_warmstart, $(sim.data.d)->qacc_warmstart, $(sim.model.m)->nv);
                    """

                    #step to get difference
                    for step=1:1
                        mj_forward(ddSim.model, ddSim.data)
                    end

                    #calc fuu
                    for nx3=1:sim.state.nx
                        fuu[nx3, nu1, nu2, t] = (dSim.state[nx3] + ddSim.state[nx3] - 2 * sim.state[nx3]) / (2.0 * 1.0e-4)^2
                    end

                    #reset
                    icxx"""
                    mju_copy($(ddSim.data.d)->qpos, $(sim.data.d)->qpos, $(sim.model.m)->nq);
                    mju_copy($(ddSim.data.d)->qvel, $(sim.data.d)->qvel, $(sim.model.m)->nv);
                    mju_copy($(ddSim.data.d)->qacc, $(sim.data.d)->qacc, $(sim.model.m)->nv);
                    mju_copy($(ddSim.data.d)->qacc_warmstart, $(sim.data.d)->qacc_warmstart, $(sim.model.m)->nv);
                    mju_copy($(ddSim.data.d)->qfrc_applied, $(sim.data.d)->qfrc_applied, $(sim.model.m)->nv);
                    mju_copy($(ddSim.data.d)->xfrc_applied, $(sim.data.d)->xfrc_applied, 6*$(sim.model.m)->nsite);
                    mju_copy($(ddSim.data.d)->ctrl, $(sim.data.d)->ctrl, $(sim.model.m)->nu);
                    mj_forward($(ddSim.model.m), $(sim.data.d));
                    """
                end
            end

            #reset
            icxx"""
            mju_copy($(dSim.data.d)->qpos, $(sim.data.d)->qpos, $(sim.model.m)->nq);
            mju_copy($(dSim.data.d)->qvel, $(sim.data.d)->qvel, $(sim.model.m)->nv);
            mju_copy($(dSim.data.d)->qacc, $(sim.data.d)->qacc, $(sim.model.m)->nv);
            mju_copy($(dSim.data.d)->qacc_warmstart, $(sim.data.d)->qacc_warmstart, $(sim.model.m)->nv);
            mju_copy($(dSim.data.d)->qfrc_applied, $(sim.data.d)->qfrc_applied, $(sim.model.m)->nv);
            mju_copy($(dSim.data.d)->xfrc_applied, $(sim.data.d)->xfrc_applied, 6*$(sim.model.m)->nsite);
            mju_copy($(dSim.data.d)->ctrl, $(sim.data.d)->ctrl, $(sim.model.m)->nu);
            mj_forward($(dSim.model.m), $(dSim.data.d));
            """
        end

        #fx,fxx, and fxu
        for nx1=1:dSim.state.nx

            dSim.state[nx1] += 1.0e-4

            #warmstart acc
            icxx"""
            mju_copy($(dSim.data.d)->qacc_warmstart, $(sim.data.d)->qacc_warmstart, $(sim.model.m)->nv);
            """

            #step to get difference
            for step=1:1
                mj_forward(dSim.model, dSim.data)
            end

            #calc fx
            for nx3=1:sim.state.nx
                fx[nx3, nx1, t] = (dSim.state[nx3] - sim.state[nx3]) / 2.0 / 1.0e-4
            end

            if !linear
                #fxx
                for nx2=1:ddSim.state.nx

                    ddSim.state[nx2] -= 1.0e-4

                    #warmstart acc
                    icxx"""
                    mju_copy($(ddSim.data.d)->qacc_warmstart, $(sim.data.d)->qacc_warmstart, $(sim.model.m)->nv);
                    """

                    #step to get difference
                    for step=1:1
                        mj_forward(ddSim.model, ddSim.data)
                    end

                    #calc fxx
                    for nx3=1:sim.state.nx
                        fxx[nx3, nx1, nx2, t] = (dSim.state[nx3] + ddSim.state[nx3] - 2 * sim.state[nx3]) / (2.0 * 1.0e-4)^2
                    end

                    #reset
                    icxx"""
                    mju_copy($(ddSim.data.d)->qpos, $(sim.data.d)->qpos, $(sim.model.m)->nq);
                    mju_copy($(ddSim.data.d)->qvel, $(sim.data.d)->qvel, $(sim.model.m)->nv);
                    mju_copy($(ddSim.data.d)->qacc, $(sim.data.d)->qacc, $(sim.model.m)->nv);
                    mju_copy($(ddSim.data.d)->qacc_warmstart, $(sim.data.d)->qacc_warmstart, $(sim.model.m)->nv);
                    mju_copy($(ddSim.data.d)->qfrc_applied, $(sim.data.d)->qfrc_applied, $(sim.model.m)->nv);
                    mju_copy($(ddSim.data.d)->xfrc_applied, $(sim.data.d)->xfrc_applied, 6*$(sim.model.m)->nsite);
                    mju_copy($(ddSim.data.d)->ctrl, $(sim.data.d)->ctrl, $(sim.model.m)->nu);
                    mj_forward($(ddSim.model.m), $(sim.data.d));
                    """
                end

                #fxu
                for nu2=1:ddSim.controls.nu

                    ddSim.controls.ctrl[nu2] -= 1.0e-4

                    #warmstart acc
                    icxx"""
                    mju_copy($(ddSim.data.d)->qacc_warmstart, $(sim.data.d)->qacc_warmstart, $(sim.model.m)->nv);
                    """

                    #step to get difference
                    for step=1:1
                        mj_forward(ddSim.model, ddSim.data)
                    end

                    #calc fxu
                    for nx3=1:sim.state.nx
                        fxu[nx3, nx1, nu2, t] = (dSim.state[nx3] + ddSim.state[nx3] - 2 * sim.state[nx3]) / (2.0 * 1.0e-4)^2
                    end

                    #reset
                    icxx"""
                    mju_copy($(ddSim.data.d)->qpos, $(sim.data.d)->qpos, $(sim.model.m)->nq);
                    mju_copy($(ddSim.data.d)->qvel, $(sim.data.d)->qvel, $(sim.model.m)->nv);
                    mju_copy($(ddSim.data.d)->qacc, $(sim.data.d)->qacc, $(sim.model.m)->nv);
                    mju_copy($(ddSim.data.d)->qacc_warmstart, $(sim.data.d)->qacc_warmstart, $(sim.model.m)->nv);
                    mju_copy($(ddSim.data.d)->qfrc_applied, $(sim.data.d)->qfrc_applied, $(sim.model.m)->nv);
                    mju_copy($(ddSim.data.d)->xfrc_applied, $(sim.data.d)->xfrc_applied, 6*$(sim.model.m)->nsite);
                    mju_copy($(ddSim.data.d)->ctrl, $(sim.data.d)->ctrl, $(sim.model.m)->nu);
                    mj_forward($(ddSim.model.m), $(sim.data.d));
                    """
                end
            end

            #reset
            icxx"""
            mju_copy($(dSim.data.d)->qpos, $(sim.data.d)->qpos, $(sim.model.m)->nq);
            mju_copy($(dSim.data.d)->qvel, $(sim.data.d)->qvel, $(sim.model.m)->nv);
            mju_copy($(dSim.data.d)->qacc, $(sim.data.d)->qacc, $(sim.model.m)->nv);
            mju_copy($(dSim.data.d)->qacc_warmstart, $(sim.data.d)->qacc_warmstart, $(sim.model.m)->nv);
            mju_copy($(dSim.data.d)->qfrc_applied, $(sim.data.d)->qfrc_applied, $(sim.model.m)->nv);
            mju_copy($(dSim.data.d)->xfrc_applied, $(sim.data.d)->xfrc_applied, 6*$(sim.model.m)->nsite);
            mju_copy($(dSim.data.d)->ctrl, $(sim.data.d)->ctrl, $(sim.model.m)->nu);
            mj_forward($(dSim.model.m), $(dSim.data.d));
            """
        end

        ###COST FUNCTION DERIVATIVES
        cu[:, t] = costu(u[:, t], sim)
        cx[:, t] = costx(x[:, t], sim)
    end

    cuu = costuu(u, sim)
    cxx = costxx(x, sim)
    cxu = cx * transpose(hcat(cu, ones(sim.controls.nu)))

#    println(sum(abs, fx))
#    println(sum(abs, fu))
#    println(sum(abs, cx))
#    println(sum(abs, cu))
#    if !linear
#        println(sum(abs, fxx))
#        println(sum(abs, fxu))
#        println(sum(abs, fuu))
#    end
#    println(sum(abs, cxx))
#    println(sum(abs, cxu))
#    println(sum(abs, cuu))

    return (fx, fu, fxx, fxu, fuu, cx, cu, cxx, cxu, cuu)
end

#get cost of being in a state
function mjDynamicsFT(xi, sim, finalCost)
    return finalCost(xi, sim)
end


text:
boxvector(a::Vector) = a
boxvector(a) = [a]

"""
    SGFNode(properties::Pair{Symbol,Any}...)

Create an SGF node with the given properties. The parameter
`properties` is optional, which means that its possible to create
empty nodes.

Additional properties can be added at any time using the function
`setindex!` or `push!`. Note how the property values are stored
in a `Vector{Any}`. This is on purpose in order to support
multi-value properties.

```jldoctest sfgnode
julia> using SmartGameFormat

julia> node = SGFNode()
SmartGameFormat.SGFNode with 0 properties

julia> node[:KM] = 6.5; # set komi property

julia> node
SmartGameFormat.SGFNode with 1 property:
  :KM => Any[6.5]
```

While `setindex!` will always overwrite existing values, the
function `push!` will instead try to append the given value to an
existing property.

```jldoctest sfgnode
julia> push!(node, :AB => "aa")
SmartGameFormat.SGFNode with 2 properties:
  :AB => Any["aa"]
  :KM => Any[6.5]

julia> push!(node, :AB => "bb")
SmartGameFormat.SGFNode with 2 properties:
  :AB => Any["aa", "bb"]
  :KM => Any[6.5]
```
"""
struct SGFNode
    properties::Dict{Symbol,Vector{Any}}
end

SGFNode() = SGFNode(Dict{Symbol,Vector{Any}}())
SGFNode(pairs::Pair{Symbol,<:Vector}...) = SGFNode(Dict{Symbol,Vector{Any}}(pairs...))
SGFNode(pair::Pair{Symbol,<:Vector}) = SGFNode(Dict{Symbol,Vector{Any}}(pair))

# convenience constructors
SGFNode(pairs::Pair{Symbol}...) = SGFNode(map(p->p.first=>boxvector(p.second), pairs)...)
SGFNode(pair::Pair{Symbol}) = SGFNode(pair.first => boxvector(pair.second))

Base.hash(n::SGFNode, h::UInt) = hash(n.properties, hash(:SGFNode, h))
Base.:(==)(a::SGFNode, b::SGFNode) = isequal(a.properties, b.properties)

@inline Base.push!(n::SGFNode, pair::Pair{Symbol,<:Vector}) =
    push!(n.properties, pair)

function Base.push!(n::SGFNode, pair::Pair{Symbol})
    if haskey(n, pair.first)
        push!(n[pair.first], pair.second)
    else
        n[pair.first] = pair.second
    end
    n
end

for fun in (:haskey, :getindex, :getkey, :get, :keys, :values,
            :length, :keytype, :valtype)
    @eval Base.@propagate_inbounds (Base.$fun)(n::SGFNode, args...) =
        (Base.$fun)(n.properties, args...)
end

Base.@propagate_inbounds Base.setindex!(n::SGFNode, v, id) =
    Base.setindex!(n.properties, boxvector(v), id)

function Base.show(io::IO, n::SGFNode)
    if haskey(io, :compact)
        if length(n.properties) == 0
            print(io, "<no properties>")
        else
            for (i,p) in enumerate(n.properties)
                val = if length(p.second) == 1
                    s = first(p.second)
                    if length(s) > 2
                        "[…]"
                    else
                        string('[', s, ']')
                    end
                else
                    string("[…|…]")
                end
                print(io, p.first, val)
                if i > 7
                    print(io, " … (", length(n.properties) - i, " more)")
                    break
                end
                i < length(n.properties) && print(io, " ")
            end
        end
    else
        print(io, typeof(n))
        # a little hack to reuse the Dict show method
        tio = IOBuffer()
        show(tio, MIME"text/plain"(), n.properties)
        str = readstring(seek(tio,0))
        str = replace(str, "Dict{Symbol,Array{Any,1}}", "")
        str = replace(str, "entry", "property")
        str = replace(str, "entries", "properties")
        print(io, str)
    end
end

# --------------------------------------------------------------------

"""
    SGFGameTree([sequence::Vector{SGFNode}], [variations::Vector{SGFGameTree}])

Create a game tree with the given `sequence` and `variations`. Both
parameters are optional, which means that its possible to create
empty game trees.

To edit a `SGFGameTree` simply manipulate the two member
variables directly. Note that `SGFGameTree` is a subtype of
`AbstractVector{SGFNode}`, where `getindex` and `setindex!`
correspond to the appropriate node in the main game path. This
means that if there are any variations, then the first variation
must denote the continuation of the main game path.

```jldoctest
julia> using SmartGameFormat

julia> t = SGFGameTree()
0-node SmartGameFormat.SGFGameTree with 0 variation(s)

julia> push!(t.sequence, SGFNode(:KM => 6.5)); t
1-node SmartGameFormat.SGFGameTree with 0 variation(s):
 KM[6.5]

julia> t.variations = [SGFGameTree(SGFNode(:C=>"first")), SGFGameTree(SGFNode(:C=>"second"))]; t
2-node SmartGameFormat.SGFGameTree with 2 variation(s):
 KM[6.5]
 C[…]

julia> t[1]
SmartGameFormat.SGFNode with 1 property:
  :KM => Any[6.5]

julia> t[2]
SmartGameFormat.SGFNode with 1 property:
  :C => Any["first"]
```
"""
mutable struct SGFGameTree <: AbstractVector{SGFNode}
    sequence::Vector{SGFNode}
    variations::Vector{SGFGameTree}
end

SGFGameTree(sequence::Vector{SGFNode} = SGFNode[]) =
    SGFGameTree(sequence, SGFGameTree[])
SGFGameTree(node::SGFNode) =
    SGFGameTree([node])

Base.hash(t::SGFGameTree, h::UInt) = hash(t.sequence, hash(t.variations, hash(:SGFGameTree, h)))
Base.:(==)(a::SGFGameTree, b::SGFGameTree) = isequal(a.sequence, b.sequence) && isequal(a.variations, b.variations)

function Base.size(t::SGFGameTree)
    n = length(t.sequence)
    if length(t.variations) > 0
        n += size(first(t.variations))[1]
    end
    (n,)
end

function Base.getindex(t::SGFGameTree, i::Int)
    @boundscheck checkbounds(t, i)
    if i <= length(t.sequence)
        @inbounds n = t.sequence[i]
        n
    else
        @inbounds n = t.variations[1][i - length(t.sequence)]
        n
    end
end

function Base.setindex!(t::SGFGameTree, n::SGFNode, i::Int)
    @boundscheck checkbounds(t, i)
    if i <= length(t.sequence)
        @inbounds t.sequence[i] = n
    else
        @inbounds t.variations[1][i - length(t.sequence)] = n
    end
    nothing
end

function Base.summary(t::SGFGameTree)
    string(length(t), "-node ", typeof(t), " with ", length(t.variations), " variation(s)")
end


text:
abstract type Coordinate end

mutable struct BrainVision <: Coordinate
    x::Number
    y::Number
    z::Number
end

mutable struct Talairach <: Coordinate
    x::Number
    y::Number
    z::Number
end

mutable struct SPM <: Coordinate
    x::Number
    y::Number
    z::Number
end

mutable struct UnknownCoordinate <: Coordinate
    x::Number
    y::Number
    z::Number
end

import Base.show
function show(c::S) where S <: Coordinate
    println("Coordinate: $(typeof(c)) - ($(c.x), $(c.y), $(c.z))")
end


text:
module Environments

import ..reset!

export AbstractEnvironment,
    step!,
    set_render_mode!,
    render,
    close,
    env,
    observation_type,
    observation_extrema,
    observation_length,
    action_type,
    action_length,
    action_set,
    make_env,
    fps

abstract type AbstractEnvironment{O <: AbstractVector,A} end

function step! end
function set_render_mode! end
function render end
function close end # Close render window. Gym-specific
function observation_type end
function observation_extrema end
function observation_length end
function action_type end
function action_length end
function action_set end
function fps end

fps(::Type{<:AbstractEnvironment}) = 60

include("GymWrappers.jl")
using .GymWrappers:
    GymCartPoleV1,
    GymLunarLander,
    GymMountainCarDiscrete,
    GymMountainCarContinuous,
    GymHopperV2,
    GymAcrobot,
    GymBreakoutRAM,
    GymInvertedPendulum,
    GymCartPoleBulletV1,
    GymReacherBulletV0,
    GymHopperBulletV0,
    GymHalfCheetahBulletV0
export GymCartPoleV1,
    GymLunarLander,
    GymMountainCarDiscrete,
    GymMountainCarContinuous,
    GymHopperV2,
    GymAcrobot,
    GymBreakoutRAM,
    GymInvertedPendulum,
    GymCartPoleBulletV1,
    GymReacherBulletV0,
    GymHopperBulletV0,
    GymHalfCheetahBulletV0

include("CartPole.jl")
export CartPole, CartPoleV1

const environment_table = Dict{Symbol,Type{<:AbstractEnvironment}}(
    :GymCartPoleV1 => GymCartPoleV1,
    :GymLunarLander => GymLunarLander,
    :CartPole => CartPole,
    :GymMountainCarDiscrete => GymMountainCarDiscrete,
    :GymMountainCarContinuous => GymMountainCarContinuous,
    :GymAcrobot => GymAcrobot,
    :GymHopperV2 => GymHopperV2,
    :GymBreakoutRAM => GymBreakoutRAM,
    :GymInvertedPendulum => GymInvertedPendulum,
    :GymCartPoleBulletV1 => GymCartPoleBulletV1,
    :GymReacherBulletV0 => GymReacherBulletV0,
    :GymHopperBulletV0 => GymHopperBulletV0,
    :GymHalfCheetahBulletV0 => GymHalfCheetahBulletV0
)

function env(name::String; kargs...)
    environment_table[Symbol(name)](kargs...)
end

# TODO: Check for the presence of required keys, default values for non-essential ones
function env(spec::Dict{Symbol,Any})
    env(spec[:name]; spec[:args]...)
end

function make_env(envspec::Dict{Symbol,Any})
    environment_table[Symbol(envspec[:name])]()
end

end


text:
"""
# module BasicIO

# Examples

```jldoctest
julia>
```
"""
module BasicIO

using Reexport

include("")
@reexport using .

end


text:
import GeometryTypes

struct Face
    points::Vector{Vector{Float64}}
end
Face(args...) = Face([args...])

function Base.isapprox(a::Face, b::Face; atol=1e-8, rtol=1e-8)
    for i in 1:3
        ok = false
        for j in 1:3
            if !isapprox(a.points[j], b.points[((j+i-1)%3)+1]; atol=atol, rtol=rtol)
                ok = true
            end
        end
        ok && return true
    end
    return false
end

nfaces(d::Dict{<:Any, Face}) = length(d)
nfaces(d::Dict{<:Any, <:Vector}) = sum(map(length, values(d)))
function test_decompose(p::Polyhedra.Mesh, d::Dict)
    @test GeometryTypes.isdecomposable(GeometryTypes.Point{3, Float32}, typeof(p))
    points = GeometryTypes.decompose(GeometryTypes.Point{3, Float32}, p)
    @test GeometryTypes.isdecomposable(GeometryTypes.Face{3, Int}, typeof(p))
    faces = GeometryTypes.decompose(GeometryTypes.Face{3, Int}, p)
    @test GeometryTypes.isdecomposable(GeometryTypes.Normal{3,Float64}, typeof(p))
    normals = GeometryTypes.decompose(GeometryTypes.Normal{3,Float64}, p)
    nf = nfaces(d)
    @test length(points) == 3nf
    @test length(faces) == nf
    @test length(normals) == 3nf
    for i in eachindex(faces)
        a, b, c = faces[i]
        @test normals[3i] == normals[3i-1] == normals[3i-2]
        normal = map(x -> x === -0.0 ? 0.0 : x, round.(normals[3i], digits=6))
        face = Face(points[3i-2],
                    points[3i-1],
                    points[3i])
        @test normal in keys(d)
        f = d[normal]
        if f isa Face
            @test f ≈ face
            delete!(d, normal)
        else
            j = findfirst(fa -> fa ≈ face, f)
            @test j isa Int && j > 0
            deleteat!(f, j)
            if isempty(f)
                delete!(d, normal)
            end
        end
    end
    @test isempty(d)
end
test_decompose(p::Polyhedron, d::Dict) = test_decompose(Polyhedra.Mesh(p), d)

function orthantdecomposetest(lib::Polyhedra.Library)
    v = vrep([Ray([1., 0, 0]),
              Ray([0, 1., 0]),
              Ray([0, 0, 1.])])
    p = polyhedron(v, lib)
    d = Dict([-1.0,  0.0,  0.0] => Face([0.0, 0.0, 0.0],
                                        [0.0, 2.0, 0.0],
                                        [0.0, 0.0, 2.0]),
             [ 0.0, -1.0,  0.0] => Face([0.0, 0.0, 0.0],
                                        [0.0, 0.0, 2.0],
                                        [2.0, 0.0, 0.0]),
             [ 0.0,  0.0, -1.0] => Face([0.0, 0.0, 0.0],
                                        [2.0, 0.0, 0.0],
                                        [0.0, 2.0, 0.0]))
    test_decompose(p, d)
end

function cubedecomposetest(lib::Polyhedra.Library)
    p = polyhedron(convexhull([ 1.,  1,   1],
                              [ 1., -1,   1],
                              [ 1.,  1,  -1],
                              [ 1., -1,  -1],
                              [-1.,  1,  1],
                              [-1., -1,  1],
                              [-1.,  1, -1],
                              [-1., -1, -1]), lib)
    # FIXME failing
    #p = polyhedron(convexhull(SymPoint([1., 1, 1]),
    #                          SymPoint([1., -1, 1]),
    #                          SymPoint([1., 1, -1]),
    #                          SymPoint([1., -1, -1])))
    d = Dict([ 0.0,  0.0, -1.0] => [Face([-1.0, -1.0, -1.0], [ 1.0, -1.0, -1.0], [ 1.0,  1.0, -1.0]),
                                    Face([ 1.0,  1.0, -1.0], [-1.0,  1.0, -1.0], [-1.0, -1.0, -1.0])],
             [-1.0,  0.0,  0.0] => [Face([-1.0, -1.0, -1.0], [-1.0,  1.0, -1.0], [-1.0,  1.0,  1.0]),
                                    Face([-1.0,  1.0,  1.0], [-1.0, -1.0,  1.0], [-1.0, -1.0, -1.0])],
             [ 0.0, -1.0,  0.0] => [Face([-1.0, -1.0,  1.0], [ 1.0, -1.0,  1.0], [ 1.0, -1.0, -1.0]),
                                    Face([ 1.0, -1.0, -1.0], [-1.0, -1.0, -1.0], [-1.0, -1.0,  1.0])],
             [ 0.0,  1.0,  0.0] => [Face([-1.0,  1.0, -1.0], [ 1.0,  1.0, -1.0], [ 1.0,  1.0,  1.0]),
                                    Face([ 1.0,  1.0,  1.0], [-1.0,  1.0,  1.0], [-1.0,  1.0, -1.0])],
             [ 0.0,  0.0,  1.0] => [Face([-1.0,  1.0,  1.0], [ 1.0,  1.0,  1.0], [ 1.0, -1.0,  1.0]),
                                    Face([ 1.0, -1.0,  1.0], [-1.0, -1.0,  1.0], [-1.0,  1.0,  1.0])],
             [ 1.0,  0.0,  0.0] => [Face([ 1.0, -1.0,  1.0], [ 1.0,  1.0,  1.0], [ 1.0,  1.0, -1.0]),
                                    Face([ 1.0,  1.0, -1.0], [ 1.0, -1.0, -1.0], [ 1.0, -1.0,  1.0])])
    test_decompose(p, d)
end
function largedecomposetest(lib::Polyhedra.Library)
    V = [-1 -1  1;
         -1  1  1;
          1 -1  1;
          1  1  1;
          0  1 -1;
          1  0 -1;
         -1  0 -1;
          0 -1 -1]
    R = [0 0 1]
    v = vrep(V, R)
    p = polyhedron(v, lib)
    d13 = round(1/3, digits=6)
    d23 = round(2/3, digits=6)
    md13 = round(-1/3, digits=6)
    md23 = round(-2/3, digits=6)
    d = Dict([ 1.0, 0.0, 0.0] => [Face([1.0,1.0,1.0], [1.0,-1.0,1.0], [1.0,1.0,2.0]),
                                  Face([1.0,-1.0,1.0], [1.0,-1.0,2.0], [1.0,1.0,2.0]),
                                  Face([1.0,1.0,1.0], [1.0,0.0,-1.0], [1.0,-1.0,1.0])],
             [ 0.0, 1.0, 0.0] => [Face([-1.0,1.0,1.0], [1.0,1.0,1.0], [-1.0,1.0,2.0]),
                                  Face([1.0,1.0,1.0], [1.0,1.0,2.0], [-1.0,1.0,2.0]),
                                  Face([-1.0,1.0,1.0], [0.0,1.0,-1.0], [1.0,1.0,1.0])],
             [ d23, d23,md13] => [Face([1.0,1.0,1.0], [0.0,1.0,-1.0], [1.0,0.0,-1.0])],
             [ 0.0,-1.0, 0.0] => [Face([1.0,-1.0,1.0], [-1.0,-1.0,1.0], [1.0,-1.0,2.0]),
                                  Face([-1.0,-1.0,1.0], [-1.0,-1.0,2.0], [1.0,-1.0,2.0]),
                                  Face([1.0,-1.0,1.0], [0.0,-1.0,-1.0], [-1.0,-1.0,1.0])],
             [-1.0, 0.0, 0.0] => [Face([-1.0,-1.0,1.0], [-1.0,1.0,1.0], [-1.0,-1.0,2.0]),
                                  Face([-1.0,1.0,1.0], [-1.0,1.0,2.0], [-1.0,-1.0,2.0]),
                                  Face([-1.0,-1.0,1.0], [-1.0,0.0,-1.0], [-1.0,1.0,1.0])],
             [md23,md23,md13] => [Face([-1.0,-1.0,1.0], [0.0,-1.0,-1.0], [-1.0,0.0,-1.0])],
             [md23, d23,md13] => [Face([-1.0,0.0,-1.0], [0.0,1.0,-1.0], [-1.0,1.0,1.0])],
             [ d23,md23,md13] => [Face([1.0,0.0,-1.0], [0.0,-1.0,-1.0], [1.0,-1.0,1.0])],
             [ 0.0, 0.0,-1.0] => [Face([0.0,-1.0,-1.0], [1.0,0.0,-1.0], [0.0,1.0,-1.0]),
                                  Face([0.0,1.0,-1.0], [-1.0,0.0,-1.0], [0.0,-1.0,-1.0])])
    test_decompose(p, d)
end

decomposetests = Dict("orthantdecompose" => orthantdecomposetest,
                      "cubedecompose"    => cubedecomposetest,
                      "largedecompose"   => largedecomposetest)

@polytestset decompose


text:
using Knet, KnetModules

# Model defition
type MLP <: KnetModule
    hiddens::Array{Linear, 1}
    output::Linear
    actfn::Activation
end

function MLP(ninput::Int, nhidden::Int, noutput::Int;
             nhiddens=1, actfn=ReLU())
    hiddens = []
    for i = 1:nhiddens
        push!(hiddens, Linear(nhidden, ninput))
        ninput = nhidden
    end
    output = Linear(noutput, nhidden)
    return MLP(hiddens, output, actfn)
end

function forward(ctx, m::MLP, x)
    for hidden in m.hiddens
        x = @mc m.actfn(@mc hidden(x))
    end
    return @mc m.output(x)
end

# Training epuch
function train(model, gradfn, data, optim)
    for (x, y) in data
        grads = gradfn(x, y)
        update!(model, grads, optim)
    end
end

include(Knet.dir("data","mnist.jl"))
xtrn, ytrn, xtst, ytst = mnist()
dtrn = minibatch(xtrn, ytrn, 100)
dtst = minibatch(xtst, ytst, 100)

model = MLP(784, 64, 10)
gradfn = grad(model, nll)
optim = optimizers(model, Adam)

println("epoch 0",
        "\ntrn accuracy: ", accuracy(model, dtrn),
        "\ntst accuracy: ", accuracy(model, dtst))
for epoch=1:10
    train(model, gradfn, dtrn, optim)
    println("\nepoch $epoch",
            "\ntrn accuracy: ", accuracy(model, dtrn),
            "\ntst accuracy: ", accuracy(model, dtst))
end

# Open this to play with serialization
const ser_demo = false
if ser_demo
    # Serialization
    const filename = "mnist_mlp.jld"
    info("Saving model")
    # save_module filename, model and optional addition state (here, it is optim)
    save_module(filename, model, Dict(["optim"=>optim]))

    reset_ctx!() # clear old model for memory efficiency
    model, state = load_module(filename)
    println("Accuracy of restored ", accuracy(model, dtst))

    optim = state["optim"]
    gradfn = grad(model, nll)
    train(model, gradfn, dtrn, optim)
    println("Test accuracy after one step more training ", accuracy(model, dtst))
end


text:
#=
Generates an array of draws from a discrete random variable with a
specified vector of probabilities.

@author : Spencer Lyon <spencer.lyon@nyu.edu>

@date: 2014-07-10

References
----------

http://quant-econ.net/jl/finite_markov.html?highlight=discrete_rv

TODO: as of 07/10/2014 it is not possible to define the property
      interface we see in the python version. Once issue 1974 from
      the main Julia repository is resolved, we can revisit this and
      implement that feature.
=#

"""
Generates an array of draws from a discrete random variable with
vector of probabilities given by q.

##### Fields

- `q::AbstractVector`: A vector of non-negative probabilities that sum to 1
- `Q::AbstractVector`: The cumulative sum of q
"""
type DiscreteRV{TV1<:AbstractVector, TV2<:AbstractVector}
    q::TV1
    Q::TV2
    function (::Type{DiscreteRV{TV1,TV2}}){TV1,TV2}(q, Q)
        abs(Q[end] - 1) > 1e-10 && error("q should sum to 1")
        new{TV1,TV2}(q, Q)
    end
end

function DiscreteRV{TV<:AbstractVector}(q::TV)
    Q = cumsum(q)
    DiscreteRV{TV,typeof(Q)}(q, Q)
end

"""
Make a single draw from the discrete distribution

##### Arguments

- `d::DiscreteRV`: The `DiscreteRV` type represetning the distribution

##### Returns

- `out::Int`: One draw from the discrete distribution
"""
draw(d::DiscreteRV) = searchsortedfirst(d.Q, rand())

"""
Make multiple draws from the discrete distribution represented by a
`DiscreteRV` instance

##### Arguments

- `d::DiscreteRV`: The `DiscreteRV` type representing the distribution
- `k::Int`:

##### Returns

- `out::Vector{Int}`: `k` draws from `d`
"""
draw(d::DiscreteRV, k::Int) = Int[draw(d) for i=1:k]


text:
module TestSetdiff

include("preamble.jl")
# using MicroCollections: SingletonSet

@testset begin
    @test_returns_first setdiff!!([0], [0]) == []
    @test_returns_first setdiff!!([0], (0,)) == []
    @test_returns_first setdiff!!([0, 1, 2], (0,), [1]) == [2]
    @test_returns_first setdiff!!(Set([0]), [0]) == Set()
    @test_returns_first setdiff!!(Set([0]), (0,)) == Set()
    @test_returns_first setdiff!!(Set([0, 1, 2]), Set([0]), [1]) == Set([2])
    @test setdiff!!(Empty(Vector), []) === Empty(Vector)
    @test setdiff!!(Empty(Set{Int}), []) === Empty(Set{Int})
    #=
    @test setdiff!!(SingletonSet((0,)), [0]) ==ₜ Set{Int}()
    @test setdiff!!(SingletonSet((0,)), [1]) ==ₜ Set([0])
    =#
end

end  # module


text:
using Test
using ResourcePools

@testset "ResourcePool" begin
    @testset "PooledResource" begin
        p = ResourcePool{Int}(3, () -> 42)
        @test length(p) == 3

        r = take!(p)
        @test typeof(r) <: AbstractPooledResource
        @test length(p) == 2
        @test ref_count(r) == 1

        release!(r)
        @test ref_count(r) == 0
        @test length(p) == 3

        r1, r2, r3 = take!(p), take!(p), take!(p)
        @test length(p) == 0
        @test_throws ArgumentError take!(p)

        release!(r1)
        @test length(p) == 1
        r = take!(p)
    end
    @testset "PooledArray" begin
        p = ResourcePool{Array{Int,2}}(3, () -> [1 2; 3 4])
        @test length(p) == 3

        a = take!(p)
        @test typeof(a) <: AbstractPooledDenseArray
        @test length(p) == 2
        @test ref_count(a) == 1
        @test sum(a) == 10

        release!(a)
        @test ref_count(a) == 0
        @test length(p) == 3
    end
end


text:
using StateSpace

println("Testing Extended Kalman Filter...")

# Define process model: Ricker growth model, with two populations
r1 = 0.03 							# 10% intrinsic growth rate for pop 1
r2 = 0.06 							# 14% intrinsic growth rate for pop 1
K1 = 20000 							# carrying capacity for pop 1
K2 = 10000 							# carrying capacity for pop 2

# process model
function f(x)
	[x[1] .* exp(r1 * (1 - (x[1] ./ K1))),
	 x[2] .* exp(r2 * (1 - (x[2] ./ K2)))]
end

V = diagm([3000.0, 5000.0]) 		# covariance of (additive) process noise

# Define observation model
g(x) = 0.01 .* x
W = diagm([100.0, 250.0])

mod = NonlinearGaussianSSM(f, V, g, W)
x0 = MvNormal([1000.0, 500.0], [100.0 0.0; 0.0 100.0])

x1 = predict(mod, x0)
# println(mean(x1))

y1 = observe(mod, x1)
u1 = update(mod, x1, mean(y1))
# println(mean(u1))

xx, yy = simulate(mod, 200, x0)
yy[1, 50] = NaN

# using PyPlot
# subplot(211)
# plot(xx')
# title("True state")
# subplot(212)
# plot(yy')
# title("Observations")
# readline()

fs = filter(mod, yy, x0)

# figure()
# for i in 1:2
# 	subplot(2, 1, i)
# 	plot(xx[i, :]')
# 	plot(mean(fs)[i, :]', "k")
# 	plot(mean(fs)[i, :]' + 1.96 * sqrt(var(fs))[i, :]', "gray")
# 	plot(mean(fs)[i, :]' - 1.96 * sqrt(var(fs))[i, :]', "gray")
# end
# readline()

y_new = fs.observations[:, end] + randn(size(mod.V(1), 1)) / 10
update!(mod, fs, y_new)

ss = smooth(mod, fs)

@assert loglikelihood(fs) < loglikelihood(ss)

println("Extended Kalman Filter passed.\n")


text:
#using TensorToolbox, LinearAlgebra

#Define tensor as multidimensional arrays and calculate its norm:

X=rand(4,3,2)
norm(X)

#Create identity and diagonal tensor:

Id=neye(2,2,2)
D=diagt([1,2,3,4])

#For two tensors of same size calculate their inner product:

X=rand(3,3,3,3);Y=rand(3,3,3,3);
innerprod(X,Y)

#Matricization of a tensor:

X=rand(4,3,2);mode=1;
A=tenmat(X,mode)
B=tenmat(X,row=[2,1],col=3) #by row modes [2,1] and column mode 3

#Fold matrix back to tensor:

X=matten(A,mode,[4,3,2])
X=matten(B,[2,1],[3],[4,3,2]) # by row modes [2,1] and column mode 3

#n-mode product of a tensor and a matrix or an array of matrices:

X=rand(5,4,3);
A=[rand(2,5),rand(2,4),rand(2,3)];
ttm(X,A[1],1)  #X times A[1] by mode 1
ttm(X,[A[1],A[2]],[1,2]) #X times A[1] by mode 1 and times A[2] by mode 2; same as ttm(X,A,-3)
ttm(X,A) #X times matrices from A by each mode

#n-mode (vector) product of a tensor and a vector or an array of vectors:

X=rand(5,4,3);
V=[rand(5),rand(4),rand(3)];
ttv(X,V[1],1)  #X times V[1] by mode 1
ttv(X,[V[1],V[2]],[1,2]) #X times V[1] by mode 1 and times V[2] by mode 2; same as ttm(X,V,-3)
ttv(X,V) #X times vectors from V by each mode

#Outer product of two tensors:

 X=rand(5,4,3,2);Y=rand(2,3,4);
 ttt(X,Y)

#Kronecker product of two tensors (straightforward generalization of Kronecker product of matrices):

X=rand(5,4,3);Y=rand(2,2,2);
tkron(X,Y)

#The n-rank and the mutlilinear rank of a tensor:

X=rand(5,4,3);
mode=2;
nrank(X,mode)
mrank(X)

#The HOSVD:

X=rand(5,4,3);
hosvd(X) #same as T1=hosvd(X,eps_abs=1e-8)
hosvd(X,eps_abs=1e-6) #discard singular values lower than 1e-5
hosvd(X,eps_rel=1e-3) #discard singular values lower than 1e-3*sigma_{max}
hosvd(X,reqrank=[2,2,2])

#The CP decomposition:

X=rand(5,4,3);
R=3; #number of components
cp_als(X,R)  #same as cp_als(X,R,init="rand",dimorder=1:ndims(X))
cp_als(X,R,init=[rand(5,3),rand(4,3),rand(3,3)]) #initialize factor matrices
cp_als(X,R,init="nvecs",dimorder=[2,1,3])

##Tensors in Tucker format

#Define tensor in Tucker format by its core tensor and factor matrices:

F=rand(5,4,3);
A=[rand(6,5),rand(6,4),rand(6,3)];
ttensor(F,A)

#Get Tucker format of a tensor by using HOSVD:

X=rand(8,9,7);
hosvd(X)
hosvd(X,reqrank=[3,3,3]) #HOSVD with predefined multilinear rank

#Create random tensor in Tucker format of size 5x4x3 and mulilinear rank (2,2,2):

X=randttensor([5,4,3],[2,2,2])

#Basic functions:

size(X)
coresize(X)
ndims(X)
norm(X)
full(X)  #Creates full tensor out of Tucker format
reorth(X) #Orthogonalize factor matrices
permutedims(X,[2,1,3])

#n-mode matricization of a tensor in Tucker format:

mode=1;
tenmat(X,mode)

#Basic operations:

X=randttensor([5,4,3],[2,2,2]);Y=randttensor([5,4,3],[3,2,1]);
innerprod(X,Y)
X+Y
X-Y
X==Y #same as isequal(X,Y)
3*X #same as mtimes(3,X)

#n-mode product of a tensor in Tucker format and a matrix or an array of matrices:

X=randttensor([5,4,3],[2,2,2]);
A=[rand(2,5),rand(2,4),rand(2,3)];
ttm(X,A[1],1)  #X times A[1] by mode 1
ttm(X,[A[1],A[2]],[1,2]) #X times A[1] by mode 1 and times A[2] by mode 2; same as ttm(X,A,-3)
ttm(X,A) #X times matrices from A by each mode

#n-mode (vector) product of a tensor in Tucker format and a vector or an array of vectors:

X=randttensor([5,4,3],[2,2,2]);
V=[rand(5),rand(4),rand(3)];
ttv(X,V[1],1)  #X times V[1] by mode 1
ttv(X,[V[1],V[2]],[1,2]) #X times V[1] by mode 1 and times V[2] by mode 2; same as ttm(X,V,-3)
ttv(X,V) #X times vectors from V by each mode

#The n-rank and the mutlilinear rank of a tensor in Tucker format:

X=randttensor([9,8,7],[5,4,3]);
mode=2;
nrank(X,mode)
mrank(X)

#HOSVD of a tensor in Tucker format:

X=randttensor([6,7,5],[4,4,4]);
hosvd(X)  #same as hosvd(X,eps_abs=1e-8)
hosvd(X,eps_abs=1e-6) #discard singular values lower than 1e-5
hosvd(X,eps_rel=1e-3) #discard singular values lower than 1e-3*sigma_{max}
hosvd(X,reqrank=[3,3,3]) #HOSVD with predefined multilinear rank

#The CP decomposition:

X=randttensor([6,7,5],[4,4,4]);
R=3; #number of components
cp_als(X,R)  #same as cp_als(X,R,init="rand",dimorder=1:ndims(X))
cp_als(X,R,init=[rand(6,3),rand(7,3),rand(5,3)]) #initialize factor matrices
cp_als(X,R,init="nvecs",dimorder=[2,1,3])

##Tensors in Kruskal format

#Define tensor in Kruskal format by its factor matrices (and vector of weights):

lambda=rand(3)
A=[rand(5,3),rand(4,3),rand(3,3)];
ktensor(A)
ktensor(lambda,A)

#Create random tensor in Kruskal format of size 5x4x3 and with 2 components:

X=randktensor([5,4,3],2)

#Basic functions:

size(X)
ndims(X)
norm(X)
full(X)  #Creates full tensor out of Kruskal format
permutedims(X,[2,1,3])
ncomponents(X) #Number of components

#n-mode matricization of a tensor in Kruskal format:

mode=1;
tenmat(X,mode)

#Basic operations:

X=randktensor([5,4,3],2);Y=randktensor([5,4,3],3);
innerprod(X,Y)
X+Y
X-Y
X==Y #same as isequal(X,Y)
3*X #same as mtimes(3,X)

#n-mode product of a tensor in Kruskal format and a matrix or an array of matrices:

X=randktensor([5,4,3],2);
A=[rand(2,5),rand(2,4),rand(2,3)];
ttm(X,A[1],1)  #X times A[1] by mode 1
ttm(X,[A[1],A[2]],[1,2]) #X times A[1] by mode 1 and times A[2] by mode 2; same as ttm(X,A,-3)
ttm(X,A) #X times matrices from A by each mode

#n-mode (vector) product of a tensor in Kruskal format and a vector or an array of vectors:

X=randktensor([5,4,3],2);
V=[rand(5),rand(4),rand(3)];
ttv(X,V[1],1)  #X times V[1] by mode 1
ttv(X,[V[1],V[2]],[1,2]) #X times V[1] by mode 1 and times V[2] by mode 2; same as ttm(X,V,-3)
ttv(X,V) #X times vectors from V by each mode

#Arrange the rank-1 components of a tensor in Kruskal format:

X=randktensor([6,5,4,3],3);
arrange(X)
arrange!(X)

#Fix sign ambiguity of a tensor in Kruskal format:

X=randktensor([6,5,4,3,4],3);
fixsigns(X)
fixsigns!(X)

#Distribute weights a tensor in Kruskal format to a specific mode:

X=randktensor([3,3,3],3);
mode=2;
redistribute(X,mode)
redistribute!(X,mode)

#The CP decomposition:

X=randktensor([6,7,5],4);
R=3; #number of components
cp_als(X,R)  #same as cp_als(X,R,init="rand",dimorder=1:ndims(X))
cp_als(X,R,init=[rand(6,3),rand(7,3),rand(5,3)]) #initialize factor matrices
cp_als(X,R,init="nvecs",dimorder=[2,1,3]);

##Tensors in Hierarchical Tucker format

#Define tensor in Hierarchical Tucker format by dimensional tree T, its transfer tensors and factor matrices:

T=dimtree(3)
B=[rand(2,3,1),rand(4,3,2)]
A=[rand(5,4),rand(4,3),rand(3,3)]
htensor(T,B,A)

#Get Tucker format of a tensor by using htrunc:

X=rand(8,9,7);
htrunc(X)
htrunc(X,maxrank=3) #hrunc with defined maximal rank

#Create random tensor in Hierarchical Tucker format of size 5x4x3:

X=randhtensor([5,4,3])

#Basic functions:

size(X)
ndims(X)
norm(X)
full(X)  #Creates full tensor out of Hierarchial Tucker format
reorth(X) #Orthogonalize factor matrices

#Basic operations:

X=randhtensor([5,4,3]);Y=randhtensor([5,4,3]);
innerprod(X,Y)
X+Y
X-Y
X==Y #same as isequal(X,Y)
3*X #same as mtimes(3,X)

#n-mode product of a tensor in Hierarchical Tucker format and a matrix or an array of matrices:

X=randhtensor([5,4,3]);
A=[rand(2,5),rand(2,4),rand(2,3)];
ttm(X,A[1],1)  #X times A[1] by mode 1
ttm(X,[A[1],A[2]],[1,2]) #X times A[1] by mode 1 and times A[2] by mode 2; same as ttm(X,A,-3)
ttm(X,A) #X times matrices from A by each mode

#n-mode (vector) product of a tensor in Hierarchical Tucker format and a vector or an array of vectors:

X=randhtensor([5,4,3]);
V=[rand(5),rand(4),rand(3)];
ttv(X,V[1],1)  #X times V[1] by mode 1
ttv(X,[V[1],V[2]],[1,2]) #X times V[1] by mode 1 and times V[2] by mode 2; same as ttm(X,V,-3)
ttv(X,V) #X times vectors from V by each mode

#The h-rank of a tensor in Hierarchical Tucker format:

X=htrunc(rand(9,8,7),maxrank=2)
hrank(X)

##Tensors in Tensor Train format

#Define tensor in TT format by its core tensors:

G=CoreCell(undef,3)
G[1]=rand(1,4,3)
G[2]=rand(3,6,4)
G[3]=rand(4,3,1)
X=TTtensor(G)

#Get TT format of a tensor by using TTsvd:

X=rand(5,4,3,2)
TTsvd(X)
TTsvd(X,reqrank=[2,2,2])

#Create random TT tensor of size 5x4x3 and TT-rank (2,2):

X=randTTtensor([5,4,3],[2,2])

#Basic functions::

size(X)
TTrank(X)
ndims(X)
norm(X)
full(X)  #Creates full tensor out of Tucker format
reorth(X)

#Basic operations:

X=randTTtensor([5,4,3],[2,2])
Y=randTTtensor([5,4,3],[3,3])

innerprod(X,Y)
X+Y
X-Y
3*X

#TTsvd of a TT tensor:

X=randTTtensor([7,6,5],[5,4])
TTsvd(X,reqrank=[3,3])

println()


text:
export Sphere1

"""
$(TYPEDEF)

Sphere1 is a S1 mechanization of one Circular rotation, with `theta in [-pi,pi)`.
"""
@defVariable Sphere1 1 (:Circular,)


text:
module DeviationsLH

using ArgCheck, DocStringExtensions, Formatting, PrettyTables
using EconometricsLH, ModelObjectsLH

export AbstractDeviation, ScalarDeviation, Deviation, RegressionDeviation, BoundsDeviation, PenaltyDeviation
export get_data_values, get_unpacked_data_values, get_model_values, get_unpacked_model_values, get_weights, get_std_errors
export set_model_values, set_weights!
export scalar_dev, scalar_devs, scalar_dev_dict, short_display, show_deviation, validate_deviation, long_description, short_description
# Regression deviation
export exclude_regressors!, is_excluded
# Deviation vectors
export DevVector, dev_exists, retrieve, scalar_deviation, scalar_devs, show_deviations

# ChangeTable
export ChangeTable, set_param_values!, show_table

include("helpers.jl");
include("deviation_types.jl");
include("deviation.jl");
include("regression_deviation.jl");
include("scalar_deviation.jl");
include("matrix_deviation.jl");
include("bounds_deviation.jl");
include("penalty_deviation.jl");
include("devvector.jl");
include("change_table.jl");

end # module


text:
#===============================================================================

    io_helper.jl - LR, June 2021

    Some useful functionality for nicer output.

===============================================================================#
struct MemoryTag
    """ Human readable memory tag.
    """
    bytes::Integer
end

function Base.show(io::IO, tag::MemoryTag)
    print(io, round(tag.bytes/1024^2, digits=3), "MB")
end


text:
# -*- coding: utf-8 -*-
# + {}
module PetriModels
import Base: show
using Petri
using ModelingToolkit
import ModelingToolkit: Constant, Variable
using MacroTools
import MacroTools: prewalk, postwalk
using Catlab.WiringDiagrams

using SemanticModels.ModelTools
using SemanticModels.ModelTools.CategoryTheory
import SemanticModels.ModelTools.CategoryTheory: ⊔, FinSetMorph
import SemanticModels.ModelTools: model

export PetriModel, model, rewrite!, PetriSpan, DPOProblem, solve, pushout, dropdown, equnion, ⊔

struct PetriModel <: AbstractModel
  model::Petri.Model
end

OpVar(s::Symbol) = Operation(Variable(s), [])

function model(::Type{PetriModel}, m::Petri.Model)
    return PetriModel(m)
end

function model(::Type{PetriModel}, d::WiringDiagram)::PetriModel
    # TODO design Multiple Dispatch Lens API
    vars = ModelTools.WiringDiagrams.wirenames(d)
    symvars = OpVar.(vars)
    byvar = Dict{Symbol, Operation}()
    homnames = Vector{Symbol}()
    transitions = map(enumerate(boxes(d))) do (i, box)
        invars = input_ports(box)
        outvars = output_ports(box)
        homname = box.value
        push!(homnames, homname)
        δ_in  =  length(invars)  > 1 ? +(OpVar.( invars)...) : OpVar.(invars[1])
        δ_out =  length(outvars) > 1 ? +(OpVar.(outvars)...) : OpVar.(outvars[1])
        return (δ_in, δ_out)
    end
    return PetriModel(Petri.Model(symvars, unique(transitions)))
end

function rewrite!(pm::PetriModel, pm2::PetriModel)
    rewrite!(pm, pm2, Dict())
end

function rewrite!(pModel::PetriModel, pModel2::PetriModel, f::Dict)
    pm = pModel.model
    pm2 = pModel2.model
    vars = map(pm.S) do s
        s.op
    end
    @show
    for i in 1:length(pm2.S)
        s = pm2.S[i]
        found = findfirst(vars .== (haskey(f, s) ? f[s].op : s.op))
        if typeof(found) == Nothing
            push!(pm.S, s)
            push!(pm.Δ, pm2.Δ[i])
            push!(pm.Λ, pm2.Λ[i])
            push!(pm.Φ, pm2.Φ[i])
        else
            pm.Δ[found] = pm2.Δ[i] == Nothing ? m.Δ[found] : m2.Δ[i]
            pm.Λ[found] = pm2.Λ[i] == Nothing ? m.Λ[found] : m2.Λ[i]
            pm.Φ[found] = pm2.Φ[i] == Nothing ? m.Φ[found] : m2.Φ[i]
        end
    end
end

struct PetriSpan{L,C,R}
    l::L
    c::C
    r::R
end

struct DPOProblem
    rule::PetriSpan
    c′::PetriModel
end

solve(p::DPOProblem) = pushout(p.rule.r, p.c′)

function equnion(a::Vector, b::Vector)
    x = copy(a)
    for item in b
        if !any(item2 -> isequal(item2, item), x)
            push!(x, item)
        end
    end
    return x
end

dictReplace(item, dict) = prewalk(i -> i in keys(dict) ? dict[i] : i, item)

function ⊔(gModel::PetriModel, hModel::PetriModel)
    g = gModel.model
    h = hModel.model
    dict = Dict(i => any(isequal.(i, g.S)) ? i : Operation(i.op, Expression[ModelingToolkit.Constant(2)]) for i in h.S)
    newS = [dictReplace(n, dict) for n in h.S]
    newΔ = [(dictReplace(n[1], dict), dictReplace(n[2], dict)) for n in h.Δ]
    newΛ = [dictReplace(n, dict) for n in h.Λ]
    newΦ = [dictReplace(n, dict) for n in h.Φ]
    PetriModel(Petri.Model(union(g.S, newS),
                equnion(g.Δ, newΔ),
                equnion(g.Λ, newΛ),
                equnion(g.Φ, newΦ)))
end

function (f::FinSetMorph)(gModel::G) where G <: PetriModel
    g = gModel.model
    length(dom(f)) == length(g.S) || throw(DomainError(g.S))
    ϕ = func(f)
    outS = Array{Operation}(undef, length(Set(ϕ.m.fun)))
    for i in dom(f)
        outS[ϕ(i)] = g.S[i]
    end
    out = deepcopy(g)
    PetriModel(Petri.Model(outS, out.Δ, out.Λ, out.Φ))
end

"""    pushout(m1::Model, m2::Model)

compute the CT pushout of two models.
"""
function pushout(pModel::PetriModel, pModel2::PetriModel)
    pm = pModel.model
    pm2 = pModel2.model
    states = equnion(pm.S, pm2.S)
    Δ = equnion(pm.Δ, pm2.Δ)
    Λ = equnion(pm.Λ, pm2.Λ)
    Φ = equnion(pm.Φ, pm2.Φ)
    return PetriModel(states, Δ, Λ, Φ)
end

"""
    dropdown(l::Model, c::Model, l′::Model)

compute c′ given l, c, l′, the formula for petrinets is

T_{c'} = T_{l'} \\setdiff f(T_l) \\cup a(T_c)
"""
function dropdown(pL::PetriModel, pC::PetriModel, pL′::PetriModel)
    pl = pL.model
    pc = pC.model
    pl′ = pL′.model
    states = union(pl′.S, pc.S)
    Δ = union(setdiff(pl′.Δ, pl.Δ), pc.Δ)
    Λ = union(setdiff(pl′.Λ, pl.Λ), pc.Λ)
    Φ = union(setdiff(pl′.Φ, pl.Φ), pc.Φ)
    return PetriModel(states, Δ, Λ, Φ)
end

end