text:
"""
    Translational{T} <: Joint{T}

    constraint limiting translational degrees of freedom

    axis: translational axis in parent offset frame
    axis_mask1: translational axis mask in pbody's frame
    axis_mask2: translational axis mask in pbody's frame
    axis_mask3: translational axis mask in pbody's frame
    vertices: points in parent can child frames
    spring: coefficient for joint spring
    damper: coefficient for joint damper
    spring_offset: nominal joint configuration
    joint_limits: lower and upper limits on the joint configuration
    spring_type: can be :linear or :sinusoidal (currently not implemented), if linear then we need joint_limits to avoid the 360° singularity.
    input: external impulse force
"""
mutable struct Translational{T,Nλ,Nb,N,Nb½,N̄λ} <: Joint{T,Nλ,Nb,N,Nb½}
    axis::SVector{3,T} # translation axis in parent frame
    axis_mask1::Adjoint{T,SVector{3,T}} # in pbody's frame
    axis_mask2::Adjoint{T,SVector{3,T}} # in pbody's frame
    axis_mask3::Adjoint{T,SVector{3,T}} # in pbody's frame
    vertices::NTuple{2,SVector{3,T}} # in pbody's & cbody's frames
    spring::T
    damper::T
    spring_offset::SVector{N̄λ,T}
    joint_limits::Vector{SVector{Nb½,T}} # lower and upper limits on the joint minimal coordinate angles
    spring_type::Symbol # the rotational springs can be :linear or :sinusoidal, if linear then we need joint_limits to avoid the 180° singularity.
    input::SVector{3,T}
end

function Translational{T,Nλ}(pbody::Node, cbody::Node;
        parent_vertex::AbstractVector=szeros(T,3),
        child_vertex::AbstractVector=szeros(T,3),
        axis::AbstractVector=szeros(T,3),
        spring=zero(T),
        damper=zero(T),
        spring_offset=szeros(T,3-Nλ),
        joint_limits=[szeros(T,0), szeros(T,0)],
        spring_type::Symbol = :linear) where {T,Nλ} #TODO: reimplement sinusoidal

    vertices = (parent_vertex, child_vertex)
    V1, V2, V3 = orthogonal_rows(axis)
    input = zeros(T,3)
    Nb½ = length(joint_limits[1])
    Nb = 2Nb½
    N̄λ = 3 - Nλ
    N = Nλ + 2Nb

    Translational{T,Nλ,Nb,N,Nb½,N̄λ}(axis, V1, V2, V3, vertices, spring, damper, spring_offset, joint_limits, spring_type, input), pbody.id, cbody.id
end


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

using MGVI

using Random

if :ModelPolyfit ∉ names(Main)
    include("test_models/model_polyfit.jl")
    import .ModelPolyfit
end

Test.@testset "test_mgvi_optimize_step" begin

    model = ModelPolyfit.model
    true_params = ModelPolyfit.true_params
    starting_point = ModelPolyfit.starting_point

    rng = MersenneTwister(145)
    data = rand(rng, model(true_params), 1)[1]

    first_iteration = mgvi_kl_optimize_step(rng,
                                            model, data, starting_point;
                                            jacobian_func=FwdRevADJacobianFunc,
                                            residual_sampler=ImplicitResidualSampler)

    next_iteration = first_iteration
    next_iteration = mgvi_kl_optimize_step(rng,
                                           model, data, next_iteration.result;
                                           jacobian_func=FwdRevADJacobianFunc,
                                           residual_sampler=ImplicitResidualSampler)

end

ext: jl


text:
struct MahalanobisDistances{T} <: AbstractVector{T}
    mahal::Vector{T}
end
@inline Base.size(m::MahalanobisDistances) = size(m.mahal)
@inline Base.length(m::MahalanobisDistances) = length(m.mahal)
@inline Base.getindex(m::MahalanobisDistances, i) = m.mahal[i]
@inline function Base.getindex(m::S, I::AbstractArray) where {S<:MahalanobisDistances}
    MahalanobisDistances(m.mahal[I])
end
@inline Base.setindex!(m::MahalanobisDistances, v, i) = m.mahal[i] = v
@inline Base.resize!(m::MahalanobisDistances, N) = Base.resize!(m.mahal, N)
@inline Base.sort!(m::MahalanobisDistances) = sort!(m.mahal)
Base.IndexStyle(::Type{<:MahalanobisDistances}) = IndexLinear()

MahalanobisDistances{T}(undef, N) where T = MahalanobisDistances(Vector{T}(undef, N))
function MahalanobisDistances(Z::AbstractMatrix{T}) where T
    m = Vector{T}(undef, size(X,1))
    @inbounds @fastmath for i ∈ eachindex(m)
        m[i] = sqrt(Z[i,1]^2 + Z[i,2]^2 + Z[i,3]^2)
    end
    MahalanobisDistances(m)
end
function MahalanobisDistances!(m::MahalanobisDistances{T}, Z::AbstractMatrix{T}) where T
    resize!(m, size(Z,1))
    @inbounds @fastmath for i ∈ eachindex(m)
        m[i] = sqrt(Z[i,1]^2 + Z[i,2]^2 + Z[i,3]^2)
    end
    m
end


text:
####################
## Trait definitions
####################
#
# i.e. implement the @traitdef macro
#
# It looks like
# @traitdef Cmp{X,Y} <: Eq{X,Y} begin
#     T = eltype(X) # associated type
#     isless(x,y) -> Bool
#     @constraints begin
#         X==Y
#     end
# end

# 1) parse the header
###
function parsecurly(def::Expr)
    # parses :(Cmp{x,y})
    # into: :Cmp, [:x,:y], :(Cmp{x,y}), ()
    name = def.args[1]
    paras = Symbol[]
    append!(paras,def.args[2:end])
    trait = def
    return name, paras, trait, :(Tuple{})
end
function parsecomp(def::Expr)
    # parses :(Cmp{x,y} <: Eq{x,y})
    # into:  :Cmp, [:x,:y], :(Cmp{x,y}), :(Tuple{Eq{x,y}})
    if  def.args[2]!=:<:
        error("not a <:")
    end
    name, paras, trait = parsecurly(def.args[1])
    supertraits = :(Tuple{$(def.args[3])})
    return name, paras, trait, supertraits
end

function parsetuple(def::Expr)
    # parses :(Cmp{x,y} <: Eq{x,y}, Sz{x}, Uz{y})
    # into   :Cmp, [:x,:y], :(Cmp{x,y}), :(Tuple{Eq{x,y},Sz{x},Uz{y}})
    name, paras, trait, supertraits = parsecomp(def.args[1])
    append!(supertraits.args, def.args[2:end])
    return name, paras, trait, supertraits
end

function parsetraithead(def::Expr)
    # Transforms
    # :(Cmp{X,Y} <: Eq{X,Y}, Tr1{X})
    # into
    # trait = :(Cmp{X,Y}
    # supertraits = :(Tuple{Eq{X,Y}, Tr1{X}})
    # paras = [:X,:Y]
    #
    # Returns:
    # :(immutable Cmp{X,Y} <: Trait{Tuple{Eq{X,Y}, Tr1{X}}} end)

    if def.head==:tuple # contains several parents
        name, paras, trait, supertraits = parsetuple(def)
    elseif def.head==:comparison # contains a <:
        name, paras, trait, supertraits = parsecomp(def)
    elseif def.head==:curly # no parents
        name, paras, trait, supertraits = parsecurly(def)
    else
        error("Interface specification error")
    end
    # check supertraits<:Traits
    for i =2:length(supertraits.args)
        st = supertraits.args[i].args[1]
        eval_curmod(:(@assert istraittype($st)))
    end
    # make :(immutable Cmp{X,Y} <: Trait{(Eq{X,Y}, Tr1{X})} end)
    out = :(immutable $trait <: Traits.Trait{$supertraits} end)
    return out, name
end

# 2) parse the function definitions
###

function parsebody(body::Expr)
    # Transforms:
    # body = quote
    #     R = g(X)
    #     f1(X,Y) -> X,Int
    #     f2(Y) -> X
    #     @constraints begin
    #          X==Y
    #     end
    # end
    #
    # into
    #

    # :(Bool[X==Y]),
    # :(...associated types...)
    isassoc(ex::Expr) = ex.head==:(=) # associated types
    isconstraints(ex::Expr) = ex.head==:macrocall # constraints

    outfns = :(Traits.FDict())
    constr = :(Bool[])
    assoc = quote end
    for ln in Lines(body)
        if isconstraints(ln)
            parseconstraints!(constr, ln)
        elseif isassoc(ln)
            push!(assoc.args, ln)
        else # the rest of the body are function signatures
            parsefnstypes!(outfns, ln)
        end
    end
    # store associated types (no need for TypeVar here):
    tmp = :(Any[])
    for ln in Lines(assoc)
        tvar = ln.args[1]
        push!(tmp.args, tvar)
    end
    push!(assoc.args, :(assoctyps = $tmp))
    return outfns, constr, assoc
end

# 2.5) parse constraints
####
# Note, @constraints is not really a macro-call.
function parseconstraints!(constr, block)
    # updates constr=Expr(:ref,...)
    if !(block.args[1]==symbol("@constraints"))
        throw(TraitException(
        "Only @constraints blocks allowed inside trait definition"))
    end
    for ln in Lines(block.args[2])
        push!(constr.args, ln)
    end
end

function parsefnstypes!(outfns, ln)
    # parse one line containing a function definition
    function parsefn(def)
        # Parse to get function signature.
        # parses f(X,Y), f{X <:T}(X,Y) and X+Y
        # into f and _f(...)

        getsymbol(fn) = symbol("__"*string(fn))

        _fn = deepcopy(def)
        if isa(def.args[1], Symbol) # f(X,Y) or X+Y
            fn = def.args[1]
            _fn.args[1] = getsymbol(fn)
        elseif def.args[1].head==:curly # f{X}(X,Y)
            fn = def.args[1].args[1]
            _fn.args[1].args[1] = getsymbol(fn)
        else
            throw(TraitException(
                  "Something went wrong parsing the trait function definition:\n$fn"))
        end
        # transform X->::X
        for i=2:length(_fn.args)
            _fn.args[i] = :(::$(_fn.args[i]))
        end
        return fn, _fn
    end
    function parseret(rettype, ln)
        # parse to get return types
        while ln.head!=:block
            ln = ln.args[end]
        end
        tmp = rettype
        rettype = Any[] # per-pend
        push!(rettype, ln.args[end])
        append!(rettype, tmp)
        return rettype
    end

    rettype = Any[]
    tuplereturn = false
    if ln.head==:tuple
        tuplereturn = true
        # several ret-types:
        # f1(X,Y) -> X,Y
        for r in ln.args[2:end]
            push!(rettype, r)
        end
        ln = ln.args[1]
    end

    if ln.head==:(->) # f1(X,Y) -> X
        rettype = parseret(rettype, ln)
        fn, _fn = parsefn(ln.args[1])
    elseif ln.head==:call # either f1(X,Y) or X + Y -> Z
        if isa(ln.args[end], Expr) && ln.args[end].head==:(->) # X + Y -> Z
            def = Expr(:call)
            append!(def.args, ln.args[1:end-1])
            if length(ln.args)==2
                append!(def.args, ln.args[end].args[1].args)
            else
                push!(def.args, ln.args[end].args[1])
            end
            rettype = parseret(rettype, ln)
        else # f1(X,Y)
            def = ln
            rettype =  Any[:(Any)]
        end
        fn, _fn = parsefn(def)
    else
        throw(TraitException(
                             "Something went wrong parsing the trait definition body with line:\n$ln"))
    end
    # if return is not a tuple, ditch the tuple
    if !tuplereturn
        rettype = :(::$(rettype[1]))
    else
        rettype = :(::Tuple{$(rettype...)})
    end
    tmp = _fn.args[2:end]
    _fn.args = _fn.args[1:1]
    push!(_fn.args, rettype)
    append!(_fn.args, tmp)
    # make _fn
    _fn = :($_fn = nothing)
    push!(outfns.args, :($fn => $_fn))
end

# 3) piece it together
###

@doc """The `@traitdef` macro is used to construct a trait.  Example:

     ```
     @traitdef MyArith{X,Y} begin
         # associated types
         Z = promote_type(X,Y)
         D = (X,Y)<:(Integer, Integer) ? Float64 : Z

         # method signatures
         +(X,Y) -> Z
         -(X,Y) -> Z
         *(X,Y) -> Z
         /(X,Y) -> D

         # constraints on X,Y
         @constraints begin
             X<:Number
             Y<:Number
         end
     end
     istrait(MyArith{Int, Int8}) # -> true
     ```

     - Assignments are for associated types, here `Z,D`.  These are
       types which can be calculated from the input types `X,Y`

     - Function signature definitions which are of the form `fn(X,Y,
       Other-Types) -> Return-Types`.  The return types can be left away

     - Constraints are marked in a block `@constaints`.  The are
       constraints in terms of the input types `X,Y` and are evaluated
       at trait checking.

     Traits can subtrait others:

     ```
     @traitdef MyInv{X}<:MyArith{X,X} begin
         inv(X) -> X
     end
     istrait(MyInv{Int}) # -> false
     istrait(MyInv{Float64}) # -> true
     ```
     """ ->
macro traitdef(head, body)
    ## make Trait type
    traithead, name = parsetraithead(head)
    # make the body
    meths, constr, assoc = parsebody(body)
    # make sure a generic function of all associated types exisits

    traitbody = quote
        methods::Traits.FDict
        constraints::Vector{Bool}
        assoctyps::Vector{Any}
        function $((name))()
            $assoc
            new( $meths, $constr, assoctyps)
        end
    end
    # add body to the type definition
    traithead.args[3] = traitbody
    return esc(traithead)
end


text:
#=
# A random node reordering.
=#
export reorder_grid_random!, random_order

function random_order(N, seed)
    rng = Random.MersenneTwister(seed); # Use a specific RNG seed to make results consistent
    return randperm(rng, N); # Randomly permute node order
end

function reorder_grid_random!(grid, seed = 17)
    nnodes = size(grid.allnodes,2);

    new_order= random_order(nnodes, seed);
    reorder_grid_nodes!(grid, new_order);
    return grid;
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

"""
    grad(s::SymbolicNode, wrt::Vector{Symbol})

Get the autodiff gradient of the current `SymbolicNode`. This function can
only be used if the current symbol is a loss function.

# Arguments:
* `s::SymbolicNode`: current node.
* `wrt::Vector{Symbol}`: the names of the arguments to the gradient.

Returns a gradient symbol of the corresponding gradient.
"""
function grad(s::SymbolicNode, wrt::Vector{Symbol})
  hdr_ref = Ref{MX_handle}(C_NULL)
  keys = string.(key)

  @mxcall(:MXSymbolGrad, (MX_handle, MX_uint, char_pp, Ptr{MX_handle}),
          self, length(keys), keys, hdr_ref)
  return SymbolicNode(MX_SymbolHandle(hdr_ref[]))
end

function _build_shapes(shape_size::MX_uint, shape_ndim::Ptr{MX_uint}, shape_data::Ptr{Ptr{MX_uint}})
  shape_ndim = unsafe_wrap(Array, shape_ndim, shape_size)
  shape_data = unsafe_wrap(Array, shape_data, shape_size)
  shapes = map(1:shape_size) do i
    my_shape = unsafe_wrap(Array, shape_data[i], shape_ndim[i])
    tuple(reverse(Int[my_shape...], dims = 1)...)
  end
  convert(Vector{Tuple}, shapes)
end

function _infer_shape(self, keys, indptr, sdata)
  ref_arg_shape_size = Ref{MX_uint}(0)
  ref_arg_shape_ndim = Ref{Ptr{MX_uint}}(0)
  ref_arg_shape_data = Ref{Ptr{Ptr{MX_uint}}}(0)
  ref_out_shape_size = Ref{MX_uint}(0)
  ref_out_shape_ndim = Ref{Ptr{MX_uint}}(0)
  ref_out_shape_data = Ref{Ptr{Ptr{MX_uint}}}(0)
  ref_aux_shape_size = Ref{MX_uint}(0)
  ref_aux_shape_ndim = Ref{Ptr{MX_uint}}(0)
  ref_aux_shape_data = Ref{Ptr{Ptr{MX_uint}}}(0)
  ref_complete       = Ref{Cint}(0)
  @mxcall(:MXSymbolInferShape,
          (MX_handle, MX_uint, char_pp, Ptr{MX_uint}, Ptr{MX_uint},
           Ref{MX_uint}, Ref{Ptr{MX_uint}}, Ref{Ptr{Ptr{MX_uint}}},
           Ref{MX_uint}, Ref{Ptr{MX_uint}}, Ref{Ptr{Ptr{MX_uint}}},
           Ref{MX_uint}, Ref{Ptr{MX_uint}}, Ref{Ptr{Ptr{MX_uint}}},
           Ref{Cint}),
          self, length(indptr)-1, keys, indptr, sdata,
          ref_arg_shape_size, ref_arg_shape_ndim, ref_arg_shape_data,
          ref_out_shape_size, ref_out_shape_ndim, ref_out_shape_data,
          ref_aux_shape_size, ref_aux_shape_ndim, ref_aux_shape_data,
          ref_complete)
  if ref_complete[] == 0
    return (nothing, nothing, nothing)
  else
    return (
      _build_shapes(ref_arg_shape_size[], ref_arg_shape_ndim[], ref_arg_shape_data[]),
      _build_shapes(ref_out_shape_size[], ref_out_shape_ndim[], ref_out_shape_data[]),
      _build_shapes(ref_aux_shape_size[], ref_aux_shape_ndim[], ref_aux_shape_data[])
    )
  end
end

"""
    infer_shape(self :: SymbolicNode, args...)
    infer_shape(self :: SymbolicNode; kwargs...)

Do shape inference according to the input shapes. The input shapes could be provided
as a list of shapes, which should specify the shapes of inputs in the same order as
the arguments returned by [`list_arguments`](@ref). Alternatively, the shape information
could be specified via keyword arguments.

Returns a 3-tuple containing shapes of all the arguments, shapes of all the outputs and
shapes of all the auxiliary variables. If shape inference failed due to incomplete
or incompatible inputs, the return value will be `(nothing, nothing, nothing)`.
"""
function infer_shape(self :: SymbolicNode; kwargs...)
  sdata  = MX_uint[]
  indptr = MX_uint[0]
  for (k,v) in kwargs
    append!(sdata, reverse([v...], dims = 1))
    push!(indptr, length(sdata))
  end
  keys = AbstractString[string(x[1]) for x in kwargs]
  _infer_shape(self, keys, indptr, sdata)
end
function infer_shape(self :: SymbolicNode, args::Union{Tuple, Cvoid}...)
  sdata  = MX_uint[]
  indptr = MX_uint[0]
  for arg in args
    if isa(arg, Cvoid); continue; end
    append!(sdata, reverse([arg...], dims = 1))
    push!(indptr, length(sdata))
  end
  keys = Ptr{char_p}(0)
  _infer_shape(self, keys, indptr, sdata)
end

function _infer_type(self, keys, arg_type_data)
  ref_in_type_size  = Ref{MX_uint}()
  ref_in_type_data  = Ref{Ptr{Cint}}()
  ref_out_type_size = Ref{MX_uint}()
  ref_out_type_data = Ref{Ptr{Cint}}()
  ref_aux_type_size = Ref{MX_uint}()
  ref_aux_type_data = Ref{Ptr{Cint}}()
  ref_complete      = Ref{Cint}()

  @mxcall(:MXSymbolInferType,
          (MX_handle, MX_uint, char_pp, Ptr{Cint},
           Ref{MX_uint}, Ref{Ptr{Cint}},
           Ref{MX_uint}, Ref{Ptr{Cint}},
           Ref{MX_uint}, Ref{Ptr{Cint}},
           Ref{Cint}),
          self, length(arg_type_data)-1, keys, arg_type_data,
          ref_in_type_size, ref_in_type_data,
          ref_out_type_size, ref_out_type_data,
          ref_aux_type_size, ref_aux_type_data,
          ref_complete)

  if ref_complete[] == 0
    return (nothing, nothing, nothing)
  else
    in_type = unsafe_wrap(Array, ref_in_type_data[], ref_in_type_size[])
    out_type = unsafe_wrap(Array, ref_out_type_data[], ref_out_type_size[])
    aux_type = unsafe_wrap(Array, ref_aux_type_data[], ref_aux_type_size[])
    return ([fromTypeFlag(TypeFlag(t)) for t in in_type],
            [fromTypeFlag(TypeFlag(t)) for t in out_type],
            [fromTypeFlag(TypeFlag(t)) for t in aux_type])
  end
end

"""
    infer_type(self :: SymbolicNode; kwargs...)
    infer_type(self :: SymbolicNode, args...)

Do type inference according to the input types. The input types could be provided
as a list of types, which should specify the types of inputs in the same order as
the arguments returned by [`list_arguments`](@ref). Alternatively, the type information
could be specified via keyword arguments.

Returns a 3-tuple containing types of all the arguments, types of all the outputs and
types of all the auxiliary variables. If type inference failed due to incomplete
or incompatible inputs, the return value will be `(nothing, nothing, nothing)`.
"""
function infer_type(self :: SymbolicNode; kwargs...)
  types = Cint[toTypeFlag(x[2]) for x in kwargs]
  keys = AbstractString[string(x[1]) for x in kwargs]
  _infer_type(self, keys, types)
end

function infer_type(self :: SymbolicNode, args :: Union{Tuple,Cvoid}...)
  types = Cint[]
  keys = Ptr{char_p}(0)

  for arg in args
    if isa(arg, Cvoid); continue; end
    push!(types, toTypeFlag(arg))
  end
  _infer_type(self, keys, types)
end


text:
using packageDevTest1
using Test

packageDevTest1.test_dx(2,1)

@testset "packageDevTest1.jl" begin
    @test packageDevTest1.newton(sin,cos,0.01) ≈ 0.0
    @test packageDevTest1.test_fun(2,3) == 7
    @test packageDevTest1.test_fun(1,1) == 3
end


text:
using GeneralizedGenerated: runtime_eval
using MacroTools: @q

export codegen

# moved to __init__
# @gg function codegen(_m::Model, _args, _data)
#     f = _codegen(type2model(_m))
#     :($f(_args, _data))
# end

function _codegen(m :: Model, expand_sums=true)
    s = symlogpdf(m)

    if expand_sums
        s = expandSums(s) |> foldConstants
    end

    code = codegen(s)

    for (v, rhs) in pairs(m.vals)
        pushfirst!(code.args, :($v = $rhs))
    end

    for v in arguments(m)
        vname = QuoteNode(v)
        pushfirst!(code.args, :($v = getproperty(_args, $vname)))
    end

    for v in stochastic(m)
        vname = QuoteNode(v)
        pushfirst!(code.args, :($v = getproperty(_data, $vname)))
    end

    f = mk_function(:((_args, _data) -> $code))

    return f
end

# @generated function _codegen(_m::Model, _args, _data)
#     type2model(_m) |> sourceCodegen() |> loadvals(_args, _data)
# end

# export sourceCodegen
# function sourceCodegen()
#     function(_m::Model)
#         body = @q begin end

#         for (x, rhs) in pairs(_m.vals)
#             push!(body.args, :($x = $rhs))
#         end

#         push!(body.args, eval(:(codegen(symlogpdf($_m)))))
#         return body
#     end
# end

export codegen
function codegen(s::Sym)
    r = codegen
    s.func == sympy.Add && begin
        @gensym add
        ex = @q begin
            $add = 0.0
        end
        for arg in s.args
            t = r(arg)
            push!(ex.args, :($add += $t))
        end
        push!(ex.args, add)
        # @show ex
        return ex
    end

    s.func == sympy.Mul && begin
        @gensym mul
        ex = @q begin
            $mul = 1.0
        end
        for arg in s.args
            t = r(arg)
            push!(ex.args, :($mul *= $t))
        end
        push!(ex.args, mul)
        return ex
    end

    s.func ∈ keys(symfuncs) && begin
        # @show s
        @gensym symfunc
        argnames = gensym.("arg" .* string.(1:length(s.args)))
        argvals = r.(s.args)
        ex = @q begin end
        for (k,v) in zip(argnames, argvals)
            push!(ex.args, :($k = $v))
        end
        f = symfuncs[s.func]
        push!(ex.args, :($symfunc = $f($(argnames...))))
        push!(ex.args, symfunc)
        return ex
    end

    s.func == sympy.Sum && begin
        @gensym sum
        @gensym Δsum
        @gensym lo
        @gensym hi

        summand = r(s.args[1])

        ex = @q begin
            $Δsum = $summand
            $sum += $Δsum
        end

        for limits in s.args[2:end]
            (ix, ixlo, ixhi) = limits.args

            ex = @q begin
                $lo = $(r(ixlo))
                $hi = $(r(ixhi))
                @inbounds @fastmath for $(r(ix)) in $lo:$hi
                    $ex
                end
            end
        end

        ex = @q begin
            let
                $sum = 0.0
                $ex
            $sum
            end
        end

        return ex
    end

    s.func == sympy.Symbol && return Symbol(string(s))
    s.func == sympy.Idx && return Symbol(string(s))
    s.func == sympy.IndexedBase && return Symbol(string(s))

    # @show s
    return convert(Expr, s)
end

# julia> r(sym(:x) + sym(:y))
# quote
#     var"##add#407" = 0.0
#     var"##add#407" += x
#     var"##add#407" += y
#     var"##add#407"
# end


text:
using JuMP, Gurobi, Statistics, StatsBase
include("../basics/main.jl")
include("../basics/validations.jl")
include("../JuMP/JuMP_models.jl")
include("../heuristics/construction_algorithms.jl")
include("../heuristics/local_search.jl")

function reactive_GRASP(CTS::Constants)
    tasks_by_w = OperationalShippingProblem(0, CTS)
    indicators = ["minimal", "maximal", "dist"]
    weights = [1/3, 1/3, 1/3]
    ind_results = Dict{String, Tuple{Float64, Int}}()
    for i in indicators
        ind_results[i] = (0,0)
    end

    worst_makespan = 0
    best_makespan = CTS.H
    best_LS = init_ls(CTS)
    for iter = 1:1000
        #init Loading Sequence and Quay Cranes
        makespan = CTS.H
        LS = init_ls(CTS)
        QC = init_qc(CTS)
        TIME = init_timer(CTS)

        #update indicator probabilities
        if iter%500 == 0
            sum_inv = 0
            for (key, value) in ind_results
                sum_inv += ((value[1]-worst_makespan)/(best_makespan-worst_makespan))^2
            end
            for i = 1:3
                weights[i] = ((ind_results[indicators[i]][1]-worst_makespan)/(best_makespan-worst_makespan))^2/sum_inv
            end
        println(indicators)
        println(weights)
        end

        #choose indicator
        ind_name = sample(indicators, Weights(weights))

        while LS.tasks_left > 0 && TIME.period < TIME.horizon_plan
            cranes_status = "idle"
            while cranes_status != "All cranes are busy" && cranes_status != "Next time period" && cranes_status != "LS is completed"
                cranes_status = randomizedConstructionHeuristic(ind_name, tasks_by_w, bj, LS, TIME, QC, CTS)
            end
            next_time_period(TIME, QC);
        end

        if check_solution(LS, QC, CTS) == true
            makespan = total_makespan(LS, CTS)
            ind_results[ind_name] = ((ind_results[ind_name][1]*ind_results[ind_name][end]+makespan)/(ind_results[ind_name][end]+1), ind_results[ind_name][end]+1)
            if makespan < best_makespan
                best_makespan = makespan
                best_LS = deepcopy(LS)
            elseif makespan > worst_makespan
            worst_makespan = makespan
            end
        end
    end

    return(best_makespan, best_LS)
end

best_makespan, best_LS = reactive_GRASP(CTS)
plot_solution(best_LS, best_makespan, CTS)

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
const pyjlrawtype = pynew()

pyjlraw_repr(self) = Py("<jl $(repr(self))>")

pyjlraw_str(self) = Py(string(self))

pyjl_attr_py2jl(k::String) = replace(k, r"_[b]+$" => (x -> "!"^(length(x) - 1)))

pyjl_attr_jl2py(k::String) = replace(k, r"!+$" => (x -> "_" * "b"^length(x)))

function pyjlraw_getattr(self, k_::Py)
    k = Symbol(pyjl_attr_py2jl(pyconvert(String, k_)))
    pyjlraw(getproperty(self, k))
end

function pyjlraw_setattr(self, k_::Py, v_::Py)
    k = Symbol(pyjl_attr_py2jl(pyconvert(String, k_)))
    v = pyconvert(Any, v_)
    setproperty!(self, k, v)
    Py(nothing)
end

pyjlraw_dir(self) = pylist(pyjl_attr_jl2py(string(k)) for k in propertynames(self, true))

function pyjlraw_call(self, args_::Py, kwargs_::Py)
    if pylen(kwargs_) > 0
        args = pyconvert(Vector{Any}, args_)
        kwargs = pyconvert(Dict{Symbol,Any}, kwargs_)
        pyjlraw(self(args...; kwargs...))
    elseif pylen(args_) > 0
        args = pyconvert(Vector{Any}, args_)
        pyjlraw(self(args...))
    else
        pyjlraw(self())
    end
end

pyjlraw_len(self) = Py(length(self))

function pyjlraw_getitem(self, k_::Py)
    if pyistuple(k_)
        k = pyconvert(Vector{Any}, k_)
        pyjlraw(self[k...])
    else
        k = pyconvert(Any, k_)
        pyjlraw(self[k])
    end
end

function pyjlraw_setitem(self, k_::Py, v_::Py)
    v = pyconvert(Any, v_)
    if pyistuple(k_)
        k = pyconvert(Vector{Any}, k_)
        self[k...] = v
    else
        k = pyconvert(Any, k_)
        self[k] = v
    end
    Py(nothing)
end

function pyjlraw_delitem(self, k_::Py)
    if pyistuple(k_)
        k = pyconvert(Vector{Any}, k_)
        delete!(self, k...)
    else
        k = pyconvert(Any, k_)
        delete!(self, k)
    end
    Py(nothing)
end

pyjlraw_bool(self::Bool) = Py(self)
pyjlraw_bool(self) = (errset(pybuiltins.TypeError, "Only Julia 'Bool' can be tested for truthyness"); pynew())

function init_jlwrap_raw()
    jl = pyjuliacallmodule
    filename = "$(@__FILE__):$(1+@__LINE__)"
    pybuiltins.exec(pybuiltins.compile("""
    class RawValue(ValueBase):
        __slots__ = ()
        __module__ = "juliacall"
        def __repr__(self):
            if self._jl_isnull():
                return "<jl NULL>"
            else:
                return self._jl_callmethod($(pyjl_methodnum(pyjlraw_repr)))
        def __str__(self):
            if self._jl_isnull():
                return "NULL"
            else:
                return self._jl_callmethod($(pyjl_methodnum(pyjlraw_str)))
        def __getattr__(self, k):
            if k.startswith("__") and k.endswith("__"):
                raise AttributeError(k)
            else:
                return self._jl_callmethod($(pyjl_methodnum(pyjlraw_getattr)), k)
        def __setattr__(self, k, v):
            try:
                ValueBase.__setattr__(self, k, v)
            except AttributeError:
                if k.startswith("__") and k.endswith("__"):
                    raise AttributeError(k)
                else:
                    self._jl_callmethod($(pyjl_methodnum(pyjlraw_setattr)), k, v)
        def __dir__(self):
            return ValueBase.__dir__(self) + self._jl_callmethod($(pyjl_methodnum(pyjlraw_dir)))
        def __call__(self, *args, **kwargs):
            return self._jl_callmethod($(pyjl_methodnum(pyjlraw_call)), args, kwargs)
        def __len__(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlraw_len)))
        def __getitem__(self, k):
            return self._jl_callmethod($(pyjl_methodnum(pyjlraw_getitem)), k)
        def __setitem__(self, k, v):
            self._jl_callmethod($(pyjl_methodnum(pyjlraw_setitem)), k, v)
        def __delitem__(self, k):
            self._jl_callmethod($(pyjl_methodnum(pyjlraw_delitem)), k)
        def __bool__(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlraw_bool)))
        def _jl_any(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjl)))
    """, filename, "exec"), jl.__dict__)
    pycopy!(pyjlrawtype, jl.RawValue)
end

"""
    pyjlraw(v)

Create a Python object wrapping the Julia object `x`.

It has type `juliacall.RawValue`. This has a much more rigid "Julian" interface than `pyjl(v)`.
For example, accessing attributes or calling this object will always return a `RawValue`.
"""
pyjlraw(v) = pyjl(pyjlrawtype, v)
export pyjlraw


text:
using JUnitTestSets, Test

function myadd(a,b)
    return a+b
end

function boolfunc(a)
    return a == 1
end

ts = @jtestset begin
    @jtest 1 == 1
    @jtest 1 == 2
    @jtest myadd(1,2) == 3
    @jtest boolfunc(1)
end

report = write_report(ts)

expected_result = """
<?xml version="1.0" encoding="utf-8"?>
<testsuite id="tests" name="Julia Tests" time="0.006513769" tests="4" failures="1">
  <testcase id="1 == 1" name="1 == 1" time="0.002504622"/>
  <testcase id="1 == 2" name="1 == 2" time="1.0295e-5">
    <failure message="Failure" type="FAIL">test failed: 1 == 2
Source: $(@__FILE__):13</failure>
  </testcase>
  <testcase id="myadd(1, 2) == 3" name="myadd(1, 2) == 3" time="0.001850544"/>
  <testcase id="boolfunc(1)" name="boolfunc(1)" time="0.002148308"/>
</testsuite>
"""

@testset "JUnitTestSets.jl testing itself" begin
    @test replace(expected_result, r"time=\"[0-9].[e\-0-9]+\""=>"time=\"0\"") == replace(string(report), r"time=\"[0-9].[e\-0-9]+\""=>"time=\"0\"")
end


text:
using Test
using AccelerometerCalibration
using AccelerometerCalibration: Calibration, calibrate!, update_calibration!
using AccelerometerCalibration: calibrate_rotation!
using Serialization
using Plots
using AccelerometerCalibrationPlots

@testset "AccelerometerCalibration" begin

    # Load raw acceleromter recording...
    raw_points = deserialize("raw_points.serialized")
    channel_count = length(raw_points[1])

    # Create a Calibration for each channel.
    cal_v = [Calibration() for i in 1:channel_count]

    # Set up callback to re-plot when a new calibration point is found.
    anim = Animation()
    AccelerometerCalibrationPlots.reset()
    cal_v[4].callback = () -> begin
        AccelerometerCalibrationPlots.calplot(cal_v)
        frame(anim)
    end

    cal_error(cal) = sum((1 - hypot((cal * v)...))^2 for v in cal.points)

    for (i, sample) in enumerate(raw_points)
        global point_count
        update_calibration!(cal_v, sample)
        if i == 20
            @test any(x->x>0.2, cal_error.(cal_v))
        end
    end

    for cal in cal_v
        err = cal_error(cal)
        @show err
        @test err < 0.02
    end

    # Calibrate again without time limit
    for cal in cal_v
        calibrate!(cal, time_limit=Inf)
    end
    calibrate_rotation!(cal_v, time_limit=Inf)
    for cal in cal_v
        err = cal_error(cal)
        @show err
        @test err < 0.02
    end

    AccelerometerCalibrationPlots.calplot(cal_v)
    frame(anim)

    gif(anim, "plot.gif", fps=5)
end


text:
include("lshaped/lshaped.jl")
include("progressivehedging/progressivehedging.jl")


text:
###################################
# Routing module
###################################

"""
Route selector

Finds a route by using a choosen mode (fastest, shortest or based on Google Distances API) and returns intersections indeces for the route

**Arguments**
* `start_node` : unique start node id selected for an agent
* `waypoint` :  unique  node id of a waypoint or nothing when agent is driving directly from *start_node* to *fin_node*
* `fin_node` :  unique finish id selected for an agent
* `activity` : string with category of  waypoint or nothing when agent is driving directly from *start_node* to *fin_node*
* `sim_data` : `SimData` object
* `buffer` : array with already chosen routes stored as a `Road` object
* `routing_mode` : string determining a way how the route will be selected (fastest, shortest or based on Google Distances API)

"""
function get_route(start_node::Int,
					waypoint::Union{Int,Nothing},
					fin_node::Int,
					activity::Union{String,Nothing},
                    sim_data::OpenStreetMapXSim.SimData,
                    buffer::Array{OpenStreetMapXSim.Road,1},
                    routing_mode::String)
    if routing_mode == "fastest"
		if isa(waypoint,Nothing)
			route_nodes, distance, route_time = OpenStreetMapX.fastest_route(sim_data.map_data.network, start_node, fin_node)
		else
			route_nodes, distance, route_time = OpenStreetMapX.fastest_route(sim_data.map_data.network, start_node, waypoint, fin_node)
		end
        road = OpenStreetMapXSim.Road(start_node,fin_node, activity, routing_mode,route_nodes, 1)
        push!(buffer,road)
        return route_nodes
    elseif routing_mode == "shortest"
		if isa(waypoint,Nothing)
			route_nodes, distance, route_time = OpenStreetMapX.shortest_route(sim_data.map_data.network, start_node, fin_node)
		else
			route_nodes, distance, route_time = OpenStreetMapX.shortest_route(sim_data.map_data.network, start_node, waypoint, fin_node)
		end
        road = OpenStreetMapXSim.Road(start_node,fin_node, activity, routing_mode,route_nodes, 1)
        push!(buffer,road)
        return route_nodes
	else
		if isa(waypoint,Nothing)
			route_nodes,routing_mode = OpenStreetMapXSim.get_google_route(start_node, fin_node, sim_data.map_data, sim_data.googleapi_key)
		else
			route_nodes,routing_mode = OpenStreetMapXSim.get_google_route(start_node, fin_node, waypoint, sim_data.map_data, sim_data.googleapi_key)
		end
        road = OpenStreetMapXSim.Road(start_node,fin_node, activity, routing_mode,route_nodes, 1)
        push!(buffer,road)
        return route_nodes
    end
end

"""
Waypoint selector

Selects waypoint by minimizing the length of the route from DA_start by waypoint to DA_fin

**Arguments**
* `start_node` :  unique start node id selected for an agent
* `fin_node` :  unique finish node id selected for an agent
* `activity` : string with category of sought waypoint
* `sim_data` : `SimData` object
* `exact` : boolean value indicing how the waypoint will be chosen (with approximate or exact algorithm)

"""
function get_waypoint(start_node::Int,
					fin_node::Int,
					activity::String,
					sim_data::OpenStreetMapXSim.SimData,
					exact::Bool)
	waypoints = OpenStreetMapX.filter_graph_features(sim_data.features, sim_data.feature_to_intersections,sim_data.feature_classes,activity)
	if exact
		return waypoint = OpenStreetMapX.find_optimal_waypoint_exact(sim_data.map_data.network, sim_data.map_data.network.w, start_node, fin_node, waypoints)
	else
		return waypoint = OpenStreetMapX.find_optimal_waypoint_approx(sim_data.map_data.network, sim_data.map_data.network.w, start_node, fin_node, waypoints)
	end
end

"""
Route module selector

Selects routing mode for two points from the following options: fastest route, shortest route or google route and returns a node indices of the choosen route

**Arguments**
* `DA_start` : unique DA unique id selected for an agent
* `DA_fin` : unique DA unique id selected for an agent
* `sim_data` : simulation data struct
* `DAs_to_intersection` : dictionary mapping each DA to nearest graph node
* `buffer` : array with already chosen routes stored as a `Road` object

**Assumptions**
- the probability of selecting each routing mode is equal

"""
function select_route(DA_start::Int, DA_fin::Int,
                    sim_data::OpenStreetMapXSim.SimData,
                    buffer::Array{OpenStreetMapXSim.Road,1}; google::Bool = false)
    start_node = sim_data.DAs_to_intersection[DA_start]
    fin_node = sim_data.DAs_to_intersection[DA_fin]
    waypoint = activity = nothing
    if google
		routing_mode = rand(["shortest", "fastest", "google"])
	else
		routing_mode = rand(["shortest", "fastest"])
	end
	indice = findfirst(road -> (road.start_node == start_node) &&  (road.fin_node == fin_node) && (road.waypoint == waypoint)  && (road.mode == routing_mode), buffer)
    if isa(indice,Nothing)
		return OpenStreetMapXSim.get_route(start_node, waypoint, fin_node, activity, sim_data, buffer, routing_mode)
    else
		buffer[indice].count += 1
		return buffer[indice].route
    end
end

"""
Route module selector for three points

Selects routing mode for three points from the following options: fastest route, shortest route or google route and returns a node indices of the choosen route

**Arguments**
* `DA_start` :  unique DA id selected for an agent
* `DA_fin` :  unique DA id selected for an agent
* `network` : routing network based on OSM data
* `DAs_to_intersection` : dictionary mapping each DA to nearest graph node
* `features` : dictionary with all features existing in simulation
* `feature_classes` : dictionary mapping each category to proper integer number
* `feature_to_intersections` : dictionary mapping each feature to nearest graph node
* `buffer` : array with already chosen routes stored as a `Road` object

**Assumptions**
- the probability of selecting each routing mode is equal
-agent is choosing a waypoint based on previously selected activity
-waypoint is approximately minimizing the length of the route from DA_start by waypoint to DA_fin

"""
function select_route(DA_start::Int, DA_fin::Int,
                    activity::String,
                    sim_data::OpenStreetMapXSim.SimData,
                    buffer::Array{OpenStreetMapXSim.Road,1}; google::Bool = false)
    start_node = sim_data.DAs_to_intersection[DA_start]
    fin_node = sim_data.DAs_to_intersection[DA_fin]
	if google
		routing_mode = rand(["shortest", "fastest", "google"])
	else
		routing_mode = rand(["shortest", "fastest"])
	end
    indice = findfirst(road -> (road.start_node == start_node) &&  (road.fin_node == fin_node) && (road.waypoint == activity)  && (road.mode == routing_mode), buffer)
    if isa(indice,Nothing)
		waypoint = OpenStreetMapXSim.get_waypoint(start_node,fin_node,activity,sim_data,false)
		return OpenStreetMapXSim.get_route(start_node, waypoint, fin_node, activity, sim_data, buffer, routing_mode)
    else
		buffer[indice].count += 1
		return buffer[indice].route
    end
end


text:
using QPnorm
using Random

rng = MersenneTwister(123)
n = 30
m = 50
P = randn(rng, n, n); P = (P + P')/2
q = randn(rng, n)
A = randn(rng, m, n); b = randn(rng, m)
r_min = 10.0
r_max = 50.0

x = find_feasible_point(A, b, r_min, r_max)
@elapsed x1 = solve(P, q, A, b, x, r_min=r_min, r_max=r_max, verbosity=1, printing_interval=1)


text:
# Square Wave Arm model

type Square <: BanditArmBase
    step::Int64                         # For internal tracking
    period::Int64                       # Period at which wave repeats
    changePoints::Dict{Int64,Float64}   # Dictionary of change points
    isRestless::Bool                    # Is this a restless arm

    function Square( _period::Int64, _cps::Dict{Int64,Float64}, _ir::Bool = true )
        new(
            0,
            _period,
            _cps,
            _ir
        )
    end
end

function pull!( arm::Square )
    arm.step    = arm.step + 1
    if arm.step > arm.period
        arm.step = 1
    end
    lastCpIndex = 0.0;
    for timeIdx in keys(arm.changePoints)
        if timeIdx <= arm.step && timeIdx > lastCpIndex
            lastCpIndex = timeIdx
        end
    end
    if lastCpIndex == 0 # No change occured
        return 0.0
    else
        return arm.changePoints[lastCpIndex]
    end
end

function tick!( arm::Square )
    if arm.isRestless
        arm.step = arm.step + 1
        if arm.step > arm.period
            arm.step = 1
        end
    end
    nothing
end

function reset!( arm::Square )
    arm.step    = 0
    nothing
end


text:
module TestExamples1DStructured

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "structured_1d_dgsem")

@testset "Structured Mesh" begin
  @trixi_testset "elixir_advection_basic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic.jl"),
      l2   = [6.0388296447998465e-6],
      linf = [3.217887726258972e-5])
  end

  @trixi_testset "elixir_advection_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonperiodic.jl"),
      l2   = [5.641921365468918e-5],
      linf = [0.00021049780975179733])
  end

  @trixi_testset "elixir_euler_source_terms.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms.jl"),
      l2   = [2.2527950196212703e-8, 1.8187357193835156e-8, 7.705669939973104e-8],
      linf = [1.6205433861493646e-7, 1.465427772462391e-7, 5.372255111879554e-7])
  end

  @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_nonperiodic.jl"),
      l2   = [3.8099996914101204e-6, 1.6745575717106341e-6, 7.732189531480852e-6],
      linf = [1.2971473393186272e-5, 9.270328934274374e-6, 3.092514399671842e-5])
  end
end

end # module


text:
"""
    Logistic(; alpha = 0.01, ni = 1000)
Logistic Regression classifier.

This struct learns classifiers using multi class softmax.
Parameter `α` indicates the learning rate, and `ni` indicates the number of learnings.

# Example
```jldoctest classification
julia> model = Logistic(alpha = 0.1)
Logistic(0.1, 1000, Matrix{Float64}(undef, 0, 0))

julia> fit!(model, x, ct)
3×3 Matrix{Float64}:
  1.80736  1.64037  -0.447735
 -1.27053  1.70026   2.57027
  4.84966 -0.473835 -1.37582

julia> println(predict(model, x))
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3]
```
"""
mutable struct Logistic
    α::Float64
    n_iter::Int64
    w::Array{Float64, 2}
    Logistic(; alpha = 0.01, ni=1000) = new(alpha, ni, Array{Float64}(undef, 0, 0))
end

function fit!(model::Logistic, x, t)
    alpha, n_iter = model.α, model.n_iter
    check_size(x, t)
    x = hcat(ones(size(x, 1), 1), x)
    w = ones(size(x, 2), size(t, 2))
    for n in 1 : n_iter
        w -= alpha * ceed(x, w, t)
    end
    model.w = w
end

function forecast(model::Logistic, x)
    w = model.w
    x = hcat(ones(size(x, 1), 1), x)
    return softmax(x*w)
end

function predict(model::Logistic, x)
    w = model.w
    x = hcat(ones(size(x, 1), 1), x)
    s = softmax(x * w)
    return [findfirst(s[i, :] .== maximum(s[i, :])) for i in 1:size(s, 1)]
end


text:
using CUDA
using KernelAbstractions
using StaticArrays

@enum RelayBufferKind begin
    SingleRelayBuffer
    DoubleRelayBuffer
end

Base.similar(::Type{A}, ::Type{FT}, dims...) where {A <: Array, FT} =
    similar(Array{FT}, dims...)
Base.similar(::Type{A}, ::Type{FT}, dims...) where {A <: CuArray, FT} =
    similar(CuArray{FT}, dims...)

device(::Union{Array, SArray, MArray}) = CPU()
device(::CUDA.CuArray) = CUDADevice()

#=
    RelayBuffer{T}(::Type{A}, kind, dims...; pinned = true)

A RelayBuffer abstracts storage for communication; it is used for
staging data and for transfers. When running on:

  - CPU -- a single buffer is used for staging and transfers can be
    initiated directly to/from it.
  - GPU -- either:
    - The transfer infrastructure is CUDA-aware and a single buffer
      can be used, just as with the CPU, or
    - The transfer infrastructure is not CUDA-aware: a double buffer
      must be used with the staging buffer on the device and the
      transfer buffer on the host.

# Arguments
- `T`: element type
- `A::Type`: what kind of array to allocate for `stage`
- `kind::RelayBufferKind`: either `SingleRelayBuffer` or
  `DoubleRelayBuffer`
- `dims...`: dimensions of the array
=#
struct RelayBuffer{T, A, Buff}
    stage::A
    transfer::Buff # Union{Nothing,Buff}

    function RelayBuffer{T}(::Type{A}, kind, dims...) where {T, A}
        stage = similar(A, T, dims...)

        if kind == SingleRelayBuffer
            transfer = nothing
        elseif kind == DoubleRelayBuffer
            transfer = zeros(T, dims...)
        else
            error("Unknown RelayBufferkind $kind")
        end

        return new{T, typeof(stage), typeof(transfer)}(stage, transfer)
    end
end

function get_stage(buf::RelayBuffer)
    return buf.stage
end

function get_transfer(buf::RelayBuffer)
    if buf.transfer === nothing
        return buf.stage
    else
        return buf.transfer
    end
end

function prepare_transfer!(
    buf::RelayBuffer;
    dependencies = nothing,
    progress = yield,
)
    if buf.transfer === nothing || length(buf.transfer) == 0
        return MultiEvent(dependencies)
    end

    event = async_copy!(
        device(buf.stage),
        buf.transfer,
        buf.stage;
        dependencies = dependencies,
        progress = progress,
    )

    return event
end

function prepare_stage!(
    buf::RelayBuffer;
    dependencies = nothing,
    progress = yield,
)
    if buf.transfer === nothing || length(buf.stage) == 0
        return MultiEvent(dependencies)
    end

    event = async_copy!(
        device(buf.stage),
        buf.stage,
        buf.transfer;
        dependencies = dependencies,
        progress = progress,
    )

    return event
end


text:
using Random, Distributions

#----- Structures -----#
abstract type Agent end
mutable struct Chartist <: Agent
    d::Int64 # Time lag
    T::Float64 # Position entry threshold
    τ::Float64 # Position exit threshold
    c::Float64 # Positive constant
    xₜ₋₁::Float64 # Previous position
    xₜ::Float64 # Next position
end
mutable struct Fundamentalist <: Agent
    vₜ::Float64 # Current perceived value
    T::Float64 # Position entry threshold
    τ::Float64 # Position exit threshold
    c::Float64 # Positive constant
    xₜ₋₁::Float64 # Previous position
    xₜ::Float64 # Next position
end
struct Parameters # Immutable
    N::Int64 # Number of agents of each type
    a::Float64 # Scale parameter for capital assignment
    dₘᵢₙ::Int64 # Minimum time delay for trend followers
    dₘₐₓ::Int64 # Maximum time delay for trend followers
    vₘᵢₙ::Float64 # Minimum offset for perceived value
    vₘₐₓ::Float64 # Maximum offset for perceived value
    Tₘᵢₙ::Float64 # Minimum threshold for entering positions
    Tₘₐₓ::Float64 # Maximum threshold for entering positions
    τₘᵢₙ::Float64 # Minimum threshold for exiting positions
    τₘₐₓ::Float64 # Maximum threshold for exiting positions
    μₙ::Float64 # Mean of the noise process in vₜ
    σₙ::Float64 # Standard deviation of the noise process in vₜ
end
#---------------------------------------------------------------------------------------------------

#----- State-dependent threshold strategies -----#
function StateDependentThresholdStrategy(agent::Agent, pₜ::Float; pₜ₋::Float)
    mₜ = typeof(agent) == Fundamentalist ? pₜ - agent.vₜ : pₜ₋ - pₜ # Mispricing
    if agent.xₜ₋₁ < 0 # Currently in a short position
        if mₜ > agent.τ # Exit threshold not met
            agent.xₜ = agent.xₜ₋₁ # Keep position and don't trade
        elseif - agent.T < mₜ < agent.τ # Exit threshold met
            agent.xₜ = 0 # Exit short position but don't enter long position
        end
    elseif agent.xₜ₋₁ > 0 # Currently in a long position
        if mₜ < - agent.τ # Exit threshold not met
            agent.xₜ = agent.xₜ₋₁ # Keep position and don't trade
        elseif - agent.τ < mₜ < agent.T # Exit threshold met
            agent.xₜ = 0 # Exit long position but don't enter short position
        end
    else # Currently in neutral position
        if - agent.T < mₜ < agent.T # Entry threshold not met
            agent.xₜ = 0 # Keep position and don't trade
        end # Else statement not necessary since the position is already set by UpdatePosition!
    end # If none of these conditions are met, then the entry thresholds have been met and the agent's position is the one given in UpdatePosition!
    ωₜ = agent.xₜ - agent.xₜ₋₁
    return ωₜ
end
#---------------------------------------------------------------------------------------------------

#----- Initialize agents -----#
function InitializeAgents(priceHistory::Vector{Int64}, parameters::Parameters)
    chartists = Vector{Chartist}(); fundamentalists = Vector{Fundamentalist}()
    for i in 1:parameters.N
        # Chartist
        d = rand(parameters.dₘᵢₙ:parameters.dₘₐₓ)
        T = rand(Uniform(parameters.Tₘᵢₙ, parameters.Tₘₐₓ))
        τ = rand(Uniform(max(-T, parameters.τₘᵢₙ), min(T, parameters.τₘₐₓ))) # Ensure -T < τ < T
        c = parameters.a * (T - τ)
        xₜ = c * sign(priceHistory[end] - priceHistory[end - d]) # History[end] is the price at the previous time step
        push!(chartists, Chartist(d, T, τ, c, 0, xₜ)) # Previous position is set to be neutral
        # Fundamentalist
        vₜ = rand(Unifrom(parameters.vₘᵢₙ, parameters.vₘₐₓ)) # The agent's current value perception
        T = rand(Uniform(parameters.Tₘᵢₙ, parameters.Tₘₐₓ))
        τ = rand(Uniform(max(-T, parameters.τₘᵢₙ), min(T, parameters.τₘₐₓ))) # Ensure -T <= τ < T
        c = parameters.a * (T - τ)
        xₜ = c * sign(vₜ - priceHistory[end]) # History[end] is the price at the previous time step
        push!(fundamentalists, Fundamentalist(vₜ, T, τ, c, 0, xₜ)) # Previous position is set to be neutral
    end
    return chartists, fundamentalists
end
#---------------------------------------------------------------------------------------------------

#----- Update positions -----#
function UpdatePosition!(agent::T, paramaters::Parameters, pₜ::Int64; pₜ₋::Int64) where T <: Agent
    agent.xₜ₋₁ = agent.xₜ # Update previous position
    if isa(agent, Chartist) # Chartist
        agent.xₜ = agent.c * (pₜ - pₜ₋) # Update current position based on the updated price
    else # Fundamentalist
        η = rand(Normal(paramaters.μₙ, paramaters.σₙ))
        agent.vₜ += η # Update value perception to be used in the evaluation of the position to take
        agent.xₜ = agent.c * (agent.vₜ - pₜ) # Update current position. Position is based on the agent's new value perception and the updated price (the agent revises their value perception after the price is updated)
    end
end
#---------------------------------------------------------------------------------------------------

#----- Simulation -----#
function Simulate(priceHistory::Vector{Int64}, horizon::Int64, parameters::Parameters; seed = 1)
    Random.seed!(seed)
    chartists, fundamentalists = InitializeAgents(priceHistory, parameters) # Initialize agents
    for t in 1:horizon # Rolling closing auctions
        for (chartist, fundamentalist) in zip(chartists, fundamentalists) # Iterate through both sets of agents
            # Update positions at the start of each day for orders to be sent. Price history has been updated for the previous orders. Agents base their positions on this updated price
            UpdatePosition!(chartist, pₜ; pₜ₋ = priceHistory[end - agent.d])
            UpdatePosition!(fundamentalist, pₜ)
            # Apply state dependent thresholds to agents' updated positions to determine if they should trade or not based on transaction costs
            ωₜᶜ = StateDependentThresholdStrategy(chartist, pₜ; pₜ₋ = priceHistory[end - agent.d]) # Chartist - after positions have been revised positions are either entered, exited or kept the same
            if ωₜᶜ != 0
                SubmitOrder(chartist.tradingGateway, Order("1", ωₜᶜ > 0 ? "Buy" : "Sell", "Market", abs(round(Int, ωₜᶜ))))
            end
            ωₜᶠ = StateDependentThresholdStrategy(agent, pₜ) # Fundamentalist - after positions have been revised positions are either entered, exited or kept the same
            if ωₜᶠ != 0
                SubmitOrder(fundamentalist.tradingGateway, Order("1", ωₜᶠ > 0 ? "Buy" : "Sell", "Market", abs(round(Int, ωₜᶠ))))
            end
        end
        push!(priceHistory, pₜ) # Update the history of prices
    end
    return priceHistory
end
#---------------------------------------------------------------------------------------------------

max_stars_repo_name: IvanJericevich/FinancialModels

max_issues_repo_name: IvanJericevich/FinancialModels

max_forks_repo_name: IvanJericevich/FinancialModels


text:
module UnitfulRecipes

using RecipesBase
using Unitful: Quantity, unit, ustrip, Unitful, dimension, Units
export @P_str

const clims_types = (:contour, :contourf, :heatmap, :surface)

#==========
Main recipe
==========#

@recipe function f(::Type{T}, x::T) where T <: AbstractArray{<:Union{Missing,<:Quantity}}
    axisletter = plotattributes[:letter]   # x, y, or z
    if (axisletter == :z) &&
        get(plotattributes, :seriestype, :nothing) ∈ clims_types
        u = get(plotattributes, :zunit, unit(eltype(x)))
        ustripattribute!(plotattributes, :clims, u)
        append_unit_if_needed!(plotattributes, :colorbar_title, u)
    end
    fixaxis!(plotattributes, x, axisletter)
end

function fixaxis!(attr, x, axisletter)
    # Attribute keys
    axislabel = Symbol(axisletter, :guide) # xguide, yguide, zguide
    axislims = Symbol(axisletter, :lims)   # xlims, ylims, zlims
    axisticks = Symbol(axisletter, :ticks) # xticks, yticks, zticks
    err = Symbol(axisletter, :error)       # xerror, yerror, zerror
    axisunit = Symbol(axisletter, :unit)   # xunit, yunit, zunit
    axis = Symbol(axisletter, :axis)       # xaxis, yaxis, zaxis
    # Get the unit
    u = pop!(attr, axisunit, unit(eltype(x)))
    # If the subplot already exists with data, get its unit
    sp = get(attr, :subplot, 1)
    if sp ≤ length(attr[:plot_object]) && attr[:plot_object].n > 0
        label = attr[:plot_object][sp][axis][:guide]
        if label isa UnitfulString
            u = label.unit
        end
        # If label was not given as an argument, reuse
        get!(attr, axislabel, label)
    end
    # Fix the attributes: labels, lims, ticks, marker/line stuff, etc.
    append_unit_if_needed!(attr, axislabel, u)
    ustripattribute!(attr, axislims, u)
    ustripattribute!(attr, axisticks, u)
    ustripattribute!(attr, err, u)
    if (axisletter == :y) && haskey(attr, :ribbon)
        ustripattribute!(attr, :ribbon, u)
    end
    fixmarkercolor!(attr)
    fixmarkersize!(attr)
    fixlinecolor!(attr)
    # Strip the unit
    ustrip.(u, x)
end

# Recipe for (x::AVec, y::AVec, z::Surface) types
const AVec = AbstractVector
const AMat{T} = AbstractArray{T,2} where T
@recipe function f(x::AVec, y::AVec, z::AMat{T}) where T <: Quantity
    u = get(plotattributes, :zunit, unit(eltype(z)))
    ustripattribute!(plotattributes, :clims, u)
    z = fixaxis!(plotattributes, z, :z)
    append_unit_if_needed!(plotattributes, :colorbar_title, u)
    x, y, z
end

# Recipe for vectors of vectors
@recipe function f(::Type{T}, x::T) where T <: AbstractVector{<:AbstractVector{<:Union{Missing,<:Quantity}}}
    axisletter = plotattributes[:letter]   # x, y, or z
    [fixaxis!(plotattributes, x, axisletter) for x in x]
end

# Recipe for bare units
@recipe function f(::Type{T}, x::T) where T <: Units
    primary := false
    Float64[]*x
end

# Recipes for functions
@recipe function f(f::Function, x::T) where T <: AVec{<:Union{Missing,<:Quantity}}
    x, f.(x)
end
@recipe function f(x::T, f::Function) where T <: AVec{<:Union{Missing,<:Quantity}}
    x, f.(x)
end
@recipe function f(x::T, y::AVec, f::Function) where T <: AVec{<:Union{Missing,<:Quantity}}
    x, y, f.(x',y)
end
@recipe function f(x::AVec, y::T, f::Function) where T <: AVec{<:Union{Missing,<:Quantity}}
    x, y, f.(x',y)
end
@recipe function f(x::T1, y::T2, f::Function) where {T1<:AVec{<:Union{Missing,<:Quantity}}, T2<:AVec{<:Union{Missing,<:Quantity}}}
    x, y, f.(x',y)
end
@recipe function f(f::Function, u::Units)
    uf = UnitFunction(f, [u])
    recipedata = RecipesBase.apply_recipe(plotattributes, uf)
    (_, xmin, xmax) = recipedata[1].args
    return f, xmin*u, xmax*u
end

"""
```julia
UnitFunction
```
A function, bundled with the assumed units of each of its inputs.

```julia
f(x, y) = x^2 + y
uf = UnitFunction(f, u"m", u"m^2")
uf(3, 2) == f(3u"m", 2u"m"^2) == 7u"m^2"
```
"""
struct UnitFunction <: Function
    f::Function
    u::Vector{Units}
end
(f::UnitFunction)(args...) = f.f((args .* f.u)...)

#===============
Attribute fixing
===============#

# Markers / lines
function fixmarkercolor!(attr)
    u = ustripattribute!(attr, :marker_z)
    ustripattribute!(attr, :clims, u)
    u == Unitful.NoUnits || append_unit_if_needed!(attr, :colorbar_title, u)
end
fixmarkersize!(attr) = ustripattribute!(attr, :markersize)
fixlinecolor!(attr) = ustripattribute!(attr, :line_z)

# strip unit from attribute[key]
function ustripattribute!(attr, key)
    if haskey(attr, key)
        v = attr[key]
        u = unit(eltype(v))
        attr[key] = ustrip.(u, v)
        return u
    else
        return Unitful.NoUnits
    end
end
# If supplied, use the unit (optional 3rd argument)
function ustripattribute!(attr, key, u)
    if haskey(attr, key)
        v = attr[key]
        if eltype(v) <: Quantity
            attr[key] = ustrip.(u, v)
        end
    end
    u
end

#=======================================
Label string containing unit information
=======================================#

abstract type AbstractProtectedString <: AbstractString end
struct ProtectedString{S} <: AbstractProtectedString
    content::S
end
struct UnitfulString{S,U} <: AbstractProtectedString
    content::S
    unit::U
end
# Minimum required AbstractString interface to work with Plots
const S = AbstractProtectedString
Base.iterate(n::S) = iterate(n.content)
Base.iterate(n::S, i::Integer) = iterate(n.content, i)
Base.codeunit(n::S) = codeunit(n.content)
Base.ncodeunits(n::S) = ncodeunits(n.content)
Base.isvalid(n::S, i::Integer) = isvalid(n.content, i)
Base.pointer(n::S) = pointer(n.content)
Base.pointer(n::S, i::Integer) = pointer(n.content, i)
"""
    P_str(s)

Creates a string that will be Protected from recipe passes.

Example:
```julia
julia> plot(rand(10)*u"m", xlabel=P"This label is protected")

julia> plot(rand(10)*u"m", xlabel=P"This label is not")
```
"""
macro P_str(s)
    return ProtectedString(s)
end

#=====================================
Append unit to labels when appropriate
=====================================#

function append_unit_if_needed!(attr, key, u::Unitful.Units)
    label = get(attr, key, nothing)
    append_unit_if_needed!(attr, key, label, u)
end
# dispatch on the type of `label`
append_unit_if_needed!(attr, key, label::ProtectedString, u) = nothing
append_unit_if_needed!(attr, key, label::UnitfulString, u) = nothing
function append_unit_if_needed!(attr, key, label::Nothing, u)
    attr[key] = UnitfulString(string(u), u)
end
function append_unit_if_needed!(attr, key, label::S, u) where {S <: AbstractString}
    if !isempty(label)
        attr[key] = UnitfulString(S(format_unit_label(label, u, get(attr, :unitformat, :round))), u)
    end
end

#=============================================
Surround unit string with specified delimiters
=============================================#
format_unit_label(l, u, f::Nothing) = string(l, ' ', u)
format_unit_label(l, u, f::Function) = f(l, u)
format_unit_label(l, u, f::AbstractString) = string(l, f, u)
format_unit_label(l, u, f::NTuple{2, <:AbstractString}) = string(l, f[1], u, f[2])
format_unit_label(l, u, f::NTuple{3, <:AbstractString}) = string(f[1], l, f[2], u, f[3])
format_unit_label(l, u, f::Char) = string(l, ' ', f, ' ', u)
format_unit_label(l, u, f::NTuple{2, Char}) = string(l, ' ', f[1], u, f[2])
format_unit_label(l, u, f::NTuple{3, Char}) = string(f[1], l, ' ', f[2], u, f[3])
format_unit_label(l, u, f::Bool) = f ? format_unit_label(l, u, :round) : format_unit_label(l, u, nothing)

const UNIT_FORMATS = Dict(
                          :round => ('(', ')'),
                          :square => ('[', ']'),
                          :curly => ('{', '}'),
                          :angle => ('<', '>'),
                          :slash => '/',
                          :slashround => (" / (", ")"),
                          :slashsquare => (" / [", "]"),
                          :slashcurly => (" / {", "}"),
                          :slashangle => (" / <", ">"),
                          :verbose => " in units of ",
                         )

format_unit_label(l, u, f::Symbol) = format_unit_label(l,u,UNIT_FORMATS[f])

end # module


text:
# __precompile__()

using ArgParse,
	Base.Dates

function parse_commandline(args)
	s = ArgParseSettings()

	@add_arg_table s begin
		"--rate", "-r"
			help = "annual rate of interest in percent, eg 3.25"
			arg_type = Float64
			default = 10.
		"--principal", "-p"
			help = "principal amount, in any currency"
			arg_type = Int
			default = 100000
		"--amortization", "-a"
			help = "amortization interval, in years; defaults to 25 years (note: if a payment amount is specified, amortization interval is ignored)"
			arg_type = Int
			default = 25
		"--payment", "-y"
			help = "payment amount; if not specified, it will be calculated based on the amortization interval. Note that for simple interest mortgages, the payment will be calculated as for monthly compounded mortgages."
			arg_type = Float64
		"--frequency", "-f"
			help = "payment frequency: m = monthly (default); w = weekly; t = two weeks; b = bimonthly. Note that for bimonthly payments, they will occur on the first and the fifteenth of each month, and the mortgage should start on the first or fifteenth of the month."
			arg_type = AbstractString
			default = "m"
		"--compounding", "-c"
			help = "compounding: s = none (simple interest; note that interest accrues daily but is not compounded); m - monthly (ie American) (this is the default); c = Canadian (ie semi-annual. Note that variable rate mortgages often compound monthly - check the fine print.)"
			arg_type = AbstractString
			default = "m"
		"--startdate", "-s"
			help = "start date for the mortgage; default is today's date"
			arg_type = Date
			default = today()
		"--term","-t"
			help = "specify the term in years; use 1/4 for 3 months and 1/2 for 6 months; enter a to indicate for the entire amortization interval (ie until paid off); default is 5 years. A term greater than the amortization will generate an error."
			arg_type =AbstractString
			default = "5"
		"output"
		help = "specify 'list' for a formatted listing, or a filename; if left blank, only a summary will be printed"
	end
	return parse_args(args, s)
end


text:
module TicTacToe

include("tictactoe.jl")

end


text:
module ConvTests

using Photon, Test

resetContext()
ctx = getContext()

getimages(s=224) = KorA(randn(ctx.dtype,s,s,3,4))

function simple_1D_model()
    model = Sequential(
        Conv1D(16, 3, relu),
        Conv1D(32, 5, relu),
        Dense(100, relu),
        Dense(10)
    )

    data = KorA(randn(ctx.dtype,100,3,4))
    pred = model(data)
    @test size(pred) == (10,4)

end

function simple_3D_model()
    model = Sequential(
        Conv3D(16, 3, relu),
        Conv3D(32, 5, relu),
        Dense(100, relu),
        Dense(10)
    )

    data = KorA(randn(ctx.dtype,20,20,20,3,4))
    pred = model(data)
    @test size(pred) == (10,4)

end

function test_base_Conv2D()
    X = getimages()

    Y = Conv2D(16, (3,3))(X)
    @test Y !== nothing

    Y = Conv2D(16, (3,3), relu)(X)
    @test Y !== nothing

    Y = Conv2D(16, 3; padding=1, strides=2, dilation=2)(X)
    @test Y !== nothing

end

function convtranspose_model()
    model = Sequential(
        Conv2D(16, (3,3), relu; padding=1, strides=2, dilation=2),
        Conv2D(32, (5,5), relu),
        Conv2DTranspose(32, (5,5)),
        Conv2DTranspose(16,3; padding=1, strides=2, dilation=1),
        Dense(50, relu),
        Dense(10)
    )
    pred = model(getimages())
    @test size(pred) == (10,4)
end

function simple_conv_model()
    model = Sequential(
        Conv2D(16, (3,3), relu; padding=1, strides=2, dilation=2),
        Conv2D(32, (5,5), relu),
        Flatten(),
        Dense(100, relu),
        Dense(10)
    )
    pred = model(getimages())
    @test size(pred) == (10,4)
end

function get_output_sizeTest()
    layer = Conv2D(16, (3,3), padding=1, strides=2, dilation=2)
    s = output_size(layer,(224,224))
    @test length(s) == 2
end

function adaptive_max_model()
    # Adaptive model
    model = Sequential(
        Conv2D(16, (3,3), relu),
        MaxPool2D(),
        Conv2D(32, (5,5), relu),
        AvgPool2D(),
        AdaptiveMaxPool((10,10)),
        Flatten(),
        Dense(100, relu),
        Dense(10)
    )

    pred = model(getimages(224))
    @test size(pred) == (10,4)

    pred = model(getimages(100))
    @test size(pred) == (10,4)

end

function adaptive_avg_model()
    # Adaptive model
    model = Sequential(
        Conv2D(16, (3,3), relu),
        MaxPool2D(),
        Conv2D(32, (5,5), relu),
        AvgPool2D(),
        AdaptiveAvgPool((10,10)),
        Flatten(),
        Dense(100, relu),
        Dense(10)
    )

    pred = model(getimages(224))
    @test size(pred) == (10,4)

    pred = model(getimages(100))
    @test size(pred) == (10,4)

end

@testset "Conv" begin
    test_base_Conv2D()
    simple_conv_model()
    convtranspose_model()
    adaptive_avg_model()
    adaptive_max_model()
    get_output_sizeTest()
    simple_3D_model()
end

end


text:
### Compute the gradients using a gradient function and matrices Js ###

"""Gradient computation for full model with ARD Kernel"""
function compute_hyperparameter_gradient(k::Kernel{T,ARDKernel},gradient_function::Function,J::Vector{LinearAlgebra.Symmetric{Float64,Matrix{Float64}}}) where {T<:Real}
    return map(gradient_function,J)
end

"""Gradient computation for sparse model with ARD Kernel"""
function compute_hyperparameter_gradient(k::Kernel{T,ARDKernel},gradient_function::Function,Jmm::Vector{Symmetric{T,Matrix{T}}},Jnm::Vector{Matrix{T}}, Jnn::Vector{Vector{T}}) where T
    return map(gradient_function,Jmm,Jnm,Jnn)
end

"""Gradient computation for full model with Iso Kernel"""
function compute_hyperparameter_gradient(k::Kernel{T,IsoKernel},gradient_function::Function,J::Symmetric{Float64,Matrix{Float64}}) where T
    return gradient_function(J)
end

"""Gradient computation for sparse model with Iso Kernel"""
function compute_hyperparameter_gradient(k::Kernel{T,IsoKernel},gradient_function::Function,Jmm::Vector{Symmetric{T,Matrix{T}}},Jnm::Vector{Matrix{T}}, Jnn::Vector{Vector{T}}) where T
    return gradient_function(Jmm,Jnm,Jnn)
end

"""Gradient computations for full multiclass model with ARD Kernel"""
function compute_hyperparameter_gradient(k::Kernel{T,ARDKernel},gradient_function::Function,J::Vector{Symmetric{T,Matrix{T}}},index::Int64) where T
    return map(gradient_function,J,fill(index,k.fields.Ndim))
end

"""Gradient computations for sparse multiclass model with ARD Kernel"""
function compute_hyperparameter_gradient(k::Kernel{T,ARDKernel},gradient_function::Function,Jmm::Vector{Symmetric{T,Matrix{T}}},Jnm::Vector{Matrix{T}}, Jnn::Vector{Vector{T}},index::Int64) where T
    return map(gradient_function,Jmm,Jnm,Jnn,fill(index,k.fields.Ndim))
end

"""Gradient computations for full multiclass model with Iso Kernel"""
function compute_hyperparameter_gradient(k::Kernel{T,IsoKernel},gradient_function::Function,J::Symmetric{T,Matrix{T}},index::Int64) where T
    return gradient_function(J,index)
end

"""Gradient computations for sparse multiclass model with Iso Kernel"""
function compute_hyperparameter_gradient(k::Kernel{T,IsoKernel},gradient_function::Function,Jmm::Symmetric{T,Matrix{T}},Jnm::Matrix{T},Jnn::Vector{T},index::Int64) where T
    return gradient_function(Jmm,Jnm,Jnn,index)
end


text:
module GZBVP1D

using LinearAlgebra
using NLsolve

include("structures.jl")
export FDGZProblem, DirichletBC, NeumannBC, PeriodicBC

include("constructors.jl")
export FDGZProblem

include("assembly.jl")
export assemble_system!

include("solve.jl")
export solve_gz

end # module


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

@doc raw"""This information is immutable after the request is created. Only the Request and Usages fields can be set on creation, other fields are derived by Kubernetes and cannot be modified by users.

    IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec(;
        extra=nothing,
        groups=nothing,
        request=nothing,
        uid=nothing,
        usages=nothing,
        username=nothing,
    )

    - extra::Dict{String, Vector{String}} : Extra information about the requesting user. See user.Info interface for details.
    - groups::Vector{String} : Group information about the requesting user. See user.Info interface for details.
    - request::Vector{UInt8} : Base64-encoded PKCS#10 CSR data
    - uid::String : UID information about the requesting user. See user.Info interface for details.
    - usages::Vector{String} : allowedUsages specifies a set of usage contexts the key will be valid for. See: https://tools.ietf.org/html/rfc5280#section-4.2.1.3      https://tools.ietf.org/html/rfc5280#section-4.2.1.12
    - username::String : Information about the requesting user. See user.Info interface for details.
"""
mutable struct IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec <: SwaggerModel
    extra::Any # spec type: Union{ Nothing, Dict{String, Vector{String}} } # spec name: extra
    groups::Any # spec type: Union{ Nothing, Vector{String} } # spec name: groups
    request::Any # spec type: Union{ Nothing, Vector{UInt8} } # spec name: request
    uid::Any # spec type: Union{ Nothing, String } # spec name: uid
    usages::Any # spec type: Union{ Nothing, Vector{String} } # spec name: usages
    username::Any # spec type: Union{ Nothing, String } # spec name: username

    function IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec(;extra=nothing, groups=nothing, request=nothing, uid=nothing, usages=nothing, username=nothing)
        o = new()
        validate_property(IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec, Symbol("extra"), extra)
        setfield!(o, Symbol("extra"), extra)
        validate_property(IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec, Symbol("groups"), groups)
        setfield!(o, Symbol("groups"), groups)
        validate_property(IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec, Symbol("request"), request)
        setfield!(o, Symbol("request"), request)
        validate_property(IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec, Symbol("uid"), uid)
        setfield!(o, Symbol("uid"), uid)
        validate_property(IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec, Symbol("usages"), usages)
        setfield!(o, Symbol("usages"), usages)
        validate_property(IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec, Symbol("username"), username)
        setfield!(o, Symbol("username"), username)
        o
    end
end # type IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec

const _property_map_IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec = Dict{Symbol,Symbol}(Symbol("extra")=>Symbol("extra"), Symbol("groups")=>Symbol("groups"), Symbol("request")=>Symbol("request"), Symbol("uid")=>Symbol("uid"), Symbol("usages")=>Symbol("usages"), Symbol("username")=>Symbol("username"))
const _property_types_IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec = Dict{Symbol,String}(Symbol("extra")=>"Dict{String, Vector{String}}", Symbol("groups")=>"Vector{String}", Symbol("request")=>"Vector{UInt8}", Symbol("uid")=>"String", Symbol("usages")=>"Vector{String}", Symbol("username")=>"String")
Base.propertynames(::Type{ IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec }) = collect(keys(_property_map_IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec))
Swagger.property_type(::Type{ IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec }, property_name::Symbol) =  _property_map_IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec[property_name]

function check_required(o::IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec)
    (getproperty(o, Symbol("request")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCertificatesV1beta1CertificateSigningRequestSpec }, name::Symbol, val)
    if name === Symbol("request")
    end
end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

################################## INSTALLATION ####################################################
# (1) Enter Package Manager via     ]
# (2) Install FMI via               add FMI       or   add "https://github.com/ThummeTo/FMI.jl"
# (3) Install FMIFlux via           add FMIFlux   or   add "https://github.com/ThummeTo/FMIFlux.jl"
################################ END INSTALLATION ##################################################

# this example covers creation and training of a ME-NeuralFMU
# here, solver step size is adaptive controlled for better training performance

using FMI
using FMIFlux
using Flux
using DifferentialEquations: Tsit5
import Plots
using Zygote

modelFMUPath = joinpath(dirname(@__FILE__), "../model/SpringPendulum1D.fmu")
realFMUPath = joinpath(dirname(@__FILE__), "../model/SpringFrictionPendulum1D.fmu")

t_start = 0.0
t_step = 0.1
t_stop = 5.0
tData = collect(t_start:t_step:t_stop)

myFMU = fmiLoad(realFMUPath)
fmiInstantiate!(myFMU; loggingOn=false)
fmiSetupExperiment(myFMU, t_start, t_stop)

fmiEnterInitializationMode(myFMU)
fmiExitInitializationMode(myFMU)

x0 = fmi2GetContinuousStates(myFMU)

realSimData = fmiSimulate(myFMU, t_start, t_stop; saveat=tData, recordValues=["mass.s", "mass.v", "mass.f", "mass.a"], setup=false)
fmiUnload(myFMU)

fmiPlot(realSimData)

myFMU = fmiLoad(modelFMUPath)

fmiInstantiate!(myFMU; loggingOn=false)
fmuSimData = fmiSimulate(myFMU, t_start, t_stop; saveat=tData, recordValues=["mass.s", "mass.v", "mass.a"])

posData = fmi2SimulationResultGetValues(realSimData, "mass.s")
velData = fmi2SimulationResultGetValues(realSimData, "mass.v")

# loss function for training
global integratorSteps
function losssum()
    global integratorSteps, problem

    solution = problem(x0)

    tNet = collect(data[1] for data in solution.u)
    posNet = collect(data[2] for data in solution.u)
    #velNet = collect(data[3] for data in solution.u)

    integratorSteps = length(tNet)

    #mse_interpolate(tData, posData, tNet, posNet, tData) # mse_interpolate(tData, velData, tNet, velNet, tData)
    Flux.mse(posData, posNet)
end

# callback function for training
global iterCB = 0
function callb()
    global iterCB += 1
    global integratorSteps

    if iterCB % 10 == 1
        avg_ls = losssum()
        @info "Loss: $(round(avg_ls, digits=5))   Avg displacement in data: $(round(sqrt(avg_ls), digits=5))   Integ.Steps: $integratorSteps"
    end

    if iterCB % 100 == 1
        fig = plotResults()
        println("Fig. update.")
        display(fig)
    end
end

function plotResults()
    solutionAfter = problem(x0, t_start)
    fig = Plots.plot(xlabel="t [s]", ylabel="mass position [m]", linewidth=2,
        xtickfontsize=12, ytickfontsize=12,
        xguidefontsize=12, yguidefontsize=12,
        legendfontsize=12, legend=:bottomright)
    Plots.plot!(fig, tData, fmi2SimulationResultGetValues(fmuSimData, "mass.s"), label="FMU", linewidth=2)
    Plots.plot!(fig, tData, posData, label="reference", linewidth=2)
    Plots.plot!(fig, collect(data[1] for data in solutionAfter.u), collect(data[2] for data in solutionAfter.u), label="NeuralFMU", linewidth=2)
    fig
end

# NeuralFMU setup
numStates = fmiGetNumberOfStates(myFMU)
additionalVRs = [fmi2String2ValueReference(myFMU, "mass.m")]
numAdditionalVRs = length(additionalVRs)

net = Chain(inputs -> fmiDoStepME(myFMU, inputs, -1.0, [], [], additionalVRs),
            Dense(numStates+numAdditionalVRs, 16, tanh),
            Dense(16, 16, tanh),
            Dense(16, numStates))

problem = ME_NeuralFMU(myFMU, net, (t_start, t_stop), Tsit5(); saveat=tData)
solutionBefore = problem(x0, t_start)
fmiPlot(problem)

# train it ...
p_net = Flux.params(problem)

optim = ADAM()
# Feel free to increase training steps or epochs for better results
Flux.train!(losssum, p_net, Iterators.repeated((), 1000), optim; cb=callb)

fmiUnload(myFMU)


text:
chrom_map = Dict{String, String}()
for i in 1:23
    chrom_map[string(i)] = string(i)
    chrom_map["chr" * string(i)] = string(i)
end
chrom_map["chrX"] = "X"
chrom_map["chrY"] = "Y"
chrom_map["chrXY"] = "XY"
chrom_map["chrMT"] = "MT"
chrom_map["chrM"] = "M"
chrom_map["X"] = "X"
chrom_map["Y"] = "Y"
chrom_map["XY"] = "XY"
chrom_map["MT"] = "MT"
chrom_map["M"] = "M"

@inline function geno_ismissing(record::VCF.Record, range::UnitRange{Int})
    return record.data[first(range)] == UInt8('.') || record.data[last(range)] == UInt8('.')
end

"""
    vcf2plink(vcffile, plinkprefix)

Convert VCF file to PLINK format. This drops multi-allelic variants and variants with missing ID. Reference alleles are set to A1.
"""
function vcf2plink(vcffile, plinkprefix)
    @assert endswith(vcffile, "vcf.gz") || endswith(vcffile, ".vcf")
    nsamples = begin
        vcfio = makestream(vcffile, "r")
        reader = VCF.Reader(vcfio)
        samples = length(header(reader).sampleID)
        makestream(plinkprefix * ".fam", "w") do io
            for (i, id) in enumerate(header(reader).sampleID)
                println(io, "$i\t$id\t0\t0\t0\t0")
            end
        end
        close(reader)
        samples
    end

    nsnps = begin
        vcfio = VCF.Reader(SnpArrays.makestream(vcffile, "r"))
        records = 0
        for record in vcfio
            length(VCF.alt(record)) > 1 && continue
            try
                ismissing(VCF.id(record))
            catch e
                continue
            end
            records += 1
        end
        close(vcfio)
        records
    end

    makestream(plinkprefix * ".bed", "w+") do io
        write(io, 0x1b6c)
        write(io, 0x01)
        write(io, Matrix{UInt8}(undef, (nsamples + 3) >> 2, nsnps))
    end
    snparray = SnpArray(plinkprefix * ".bed", nsamples, "r+")
    bimio = makestream(plinkprefix * ".bim", "w")
    vcfio = makestream(vcffile, "r")
    reader = VCF.Reader(vcfio)

    j = 0
    for record in reader
        length(VCF.alt(record)) > 1 && continue
        try
            ismissing(VCF.id(record))
        catch e
            continue
        end
        j += 1
        chrom = VCF.chrom(record)
        id = VCF.id(record)[1]
        dist = 0
        pos = VCF.pos(record)
        ref = VCF.ref(record)
        alt = VCF.alt(record)[1]
        println(bimio, "$chrom\t$id\t$dist\t$pos\t$ref\t$alt")

        gtkey = findgenokey(record, "GT")
        for (i, _) in enumerate(header(reader).sampleID)
            geno = record.genotype[i]
            # dropped field or "." => 0x2e
            if gtkey > lastindex(geno) || geno_ismissing(record, geno[gtkey])
                snparray[i,j] = 0x01
            else
                # "0" => 0x30, "1" => 0x31
                if record.data[geno[gtkey][1]] == 0x30
                    if record.data[geno[gtkey][3]] == 0x30
                        snparray[i, j] = 0x00
                    else
                        snparray[i, j] = 0x02
                    end
                elseif record.data[geno[gtkey][1]] == 0x31
                    if record.data[geno[gtkey][3]] == 0x31
                        snparray[i, j] = 0x03
                    else
                        snparray[i, j] = 0x02
                    end
                end
            end
        end
    end
    close(bimio)
    close(vcfio)
    close(reader)
    snparray
end


text:
import Plots
#import PyPlot

"""
Generates data for plotting f on [a, b]
"""
function plt_data(f, a, b, numpts=100)
    h = (b-a)/(numpts-1)
    y = zeros(1,numpts)
    x = collect(a:h:b)
    y = map(f, x)
    return (x, y)
end

"""
Displays graph of f on interval [a, b]
"""
function plt(f, a, b, numpts=100; label = "y", title="")
    h = (b-a)/(numpts-1)
    y = zeros(1,numpts)
    x = collect(a:h:b)
    y = map(f, x)
    display(Plots.plot(x, y, label = label, title = title))
    return (x, y)
end

function plt!(f, a, b, numpts=100; label = "y", title = "")
    h = (b-a)/(numpts-1)
    y = zeros(1,numpts)
    x = collect(a:h:b)
    y = map(f, x)
    display(Plots.plot!(x, y, label = label, title = title))
    return(x, y)
end

function plt!(point::Float64; label = "", title = "")
    display(Plots.plot!([point], seriestype = :vline))
end

function plot_multiple(a, b, c, r1=.001, r2=1, npts = 100)
    plt(a, r1, r2, npts, label = "pdf")
    plt!(b, r1, r2, npts, label = "cdf")
    plt!(c, r1, r2, npts, label = "dpdf")
end

"""
plot pdf, median, mode and credible interval
"""
function plot_distribution_single(pdf, median, mode, CI, z0_true, quadtype)
    # plot pdf
    b = ceil(2.5 * mode)
    z_grid = range(1e-5, stop = b, step = 0.01)
    p_grid = pdf.(z_grid);
    PyPlot.plot(z_grid, p_grid, label = "probability density function")
    # PyPlot.vlines(median, 0, pdf0(mean), label = "mean", colors = "k")
    PyPlot.vlines(median, 0, pdf(median), label = "median",  colors = "b")
    PyPlot.vlines(mode, 0, pdf(mode), label = "mode")
    PyPlot.vlines(z0_true, 0, pdf(z0_true), label = "true value",  colors = "r")
    # plot CI
    CI_x_range = range(CI[1], stop = CI[2], step = 0.01)
    CI_y_range = pdf.(CI_x_range)
    PyPlot.fill_between(CI_x_range, 0, CI_y_range, alpha = 0.3, label = "95% confidence interval")
    PyPlot.legend(fontsize=8)
    PyPlot.grid()
    PyPlot.title(quadtype, fontsize=10)
#     PyPlot.savefig("Figure/Prediction")
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using MPSKit,MPSKitModels,TensorKit,Test

th = nonsym_xxz_ham()

ts = InfiniteMPS([ℂ^3],[ℂ^48]);

(ts,envs,_) = find_groundstate(ts,th,Vumps(maxiter=400));

(energies,Bs) = excitations(th,QuasiparticleAnsatz(),Float64(pi),ts,envs);

@test energies[1] ≈ 0.41047925 atol=1e-4


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct CurrentControl <: InnerControl
        kpv::Float64
        kiv::Float64
        kffv::Float64
        rv::Float64
        lv::Float64
        kpc::Float64
        kic::Float64
        kffi::Float64
        ωad::Float64
        kad::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
    end

Parameters of an inner loop current control PID using virtual impedance.

# Arguments
- `kpv::Float64`: voltage controller proportional gain, validation range: (0, nothing)
- `kiv::Float64`: voltage controller integral gain, validation range: (0, nothing)
- `kffv::Float64`: Binary variable to enable feed-forward gain of voltage., validation range: (0, nothing)
- `rv::Float64`: virtual resistance, validation range: (0, nothing)
- `lv::Float64`: virtual inductance, validation range: (0, nothing)
- `kpc::Float64`: current controller proportional gain, validation range: (0, nothing)
- `kic::Float64`: current controller integral gain, validation range: (0, nothing)
- `kffi::Float64`: Binary variable to enable feed-forward gain of current, validation range: (0, nothing)
- `ωad::Float64`: active damping filter cutoff frequency (rad/sec), validation range: (0, nothing)
- `kad::Float64`: active damping gain, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
"""
mutable struct CurrentControl <: InnerControl
    "voltage controller proportional gain"
    kpv::Float64
    "voltage controller integral gain"
    kiv::Float64
    "Binary variable to enable feed-forward gain of voltage."
    kffv::Float64
    "virtual resistance"
    rv::Float64
    "virtual inductance"
    lv::Float64
    "current controller proportional gain"
    kpc::Float64
    "current controller integral gain"
    kic::Float64
    "Binary variable to enable feed-forward gain of current"
    kffi::Float64
    "active damping filter cutoff frequency (rad/sec)"
    ωad::Float64
    "active damping gain"
    kad::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
end

function CurrentControl(kpv, kiv, kffv, rv, lv, kpc, kic, kffi, ωad, kad, ext=Dict{String, Any}(), )
    CurrentControl(kpv, kiv, kffv, rv, lv, kpc, kic, kffi, ωad, kad, ext, [:ξd_ic, :ξq_ic, :γd_ic, :γq_ic, :ϕd_ic, :ϕq_ic], 6, )
end

function CurrentControl(; kpv, kiv, kffv, rv, lv, kpc, kic, kffi, ωad, kad, ext=Dict{String, Any}(), )
    CurrentControl(kpv, kiv, kffv, rv, lv, kpc, kic, kffi, ωad, kad, ext, )
end

# Constructor for demo purposes; non-functional.
function CurrentControl(::Nothing)
    CurrentControl(;
        kpv=0,
        kiv=0,
        kffv=0,
        rv=0,
        lv=0,
        kpc=0,
        kic=0,
        kffi=0,
        ωad=0,
        kad=0,
        ext=Dict{String, Any}(),
    )
end

"""Get CurrentControl kpv."""
get_kpv(value::CurrentControl) = value.kpv
"""Get CurrentControl kiv."""
get_kiv(value::CurrentControl) = value.kiv
"""Get CurrentControl kffv."""
get_kffv(value::CurrentControl) = value.kffv
"""Get CurrentControl rv."""
get_rv(value::CurrentControl) = value.rv
"""Get CurrentControl lv."""
get_lv(value::CurrentControl) = value.lv
"""Get CurrentControl kpc."""
get_kpc(value::CurrentControl) = value.kpc
"""Get CurrentControl kic."""
get_kic(value::CurrentControl) = value.kic
"""Get CurrentControl kffi."""
get_kffi(value::CurrentControl) = value.kffi
"""Get CurrentControl ωad."""
get_ωad(value::CurrentControl) = value.ωad
"""Get CurrentControl kad."""
get_kad(value::CurrentControl) = value.kad
"""Get CurrentControl ext."""
get_ext(value::CurrentControl) = value.ext
"""Get CurrentControl states."""
get_states(value::CurrentControl) = value.states
"""Get CurrentControl n_states."""
get_n_states(value::CurrentControl) = value.n_states

"""Set CurrentControl kpv."""
set_kpv!(value::CurrentControl, val::Float64) = value.kpv = val
"""Set CurrentControl kiv."""
set_kiv!(value::CurrentControl, val::Float64) = value.kiv = val
"""Set CurrentControl kffv."""
set_kffv!(value::CurrentControl, val::Float64) = value.kffv = val
"""Set CurrentControl rv."""
set_rv!(value::CurrentControl, val::Float64) = value.rv = val
"""Set CurrentControl lv."""
set_lv!(value::CurrentControl, val::Float64) = value.lv = val
"""Set CurrentControl kpc."""
set_kpc!(value::CurrentControl, val::Float64) = value.kpc = val
"""Set CurrentControl kic."""
set_kic!(value::CurrentControl, val::Float64) = value.kic = val
"""Set CurrentControl kffi."""
set_kffi!(value::CurrentControl, val::Float64) = value.kffi = val
"""Set CurrentControl ωad."""
set_ωad!(value::CurrentControl, val::Float64) = value.ωad = val
"""Set CurrentControl kad."""
set_kad!(value::CurrentControl, val::Float64) = value.kad = val
"""Set CurrentControl ext."""
set_ext!(value::CurrentControl, val::Dict{String, Any}) = value.ext = val
"""Set CurrentControl states."""
set_states!(value::CurrentControl, val::Vector{Symbol}) = value.states = val
"""Set CurrentControl n_states."""
set_n_states!(value::CurrentControl, val::Int64) = value.n_states = val


text:
const cellml_ns = "http://www.cellml.org/cellml/1.0#"
const mathml_ns = "http://www.w3.org/1998/Math/MathML"

create_var(x) = Num(Variable(Symbol(x))).val
create_var(x, iv) = Num(Sym{FnType{Tuple{Real}}}(Symbol(x))(Variable(Symbol(iv)))).val
create_param(x) = Num(Sym{ModelingToolkit.Parameter{Real}}(Symbol(x))).val

# list the name of the CellML initialized variables, i.e., state variables
# and parameters with an initial_value tag
function find_initialized_variables(xml)
    return findall("//x:variable[@initial_value]", root(xml), ["x"=>cellml_ns])
end

# list the name of the state variables, i.e., variables that occur on the left
# hand side of an assignment
function find_state_names(xml)
    # nodes = findall("//x:math//x:eq/following-sibling::x:apply[1]/x:ci", root(xml), ["x"=>mathml_ns])
    # vars = nodecontent.(nodes)
    nodes = findall("//x:math//x:apply/x:diff", root(xml), ["x"=>mathml_ns])
    vars = [nodecontent(lastelement(parentnode(n))) for n in nodes]
    return unique(vars)
end

# find the unique independent variable
function find_iv(xml)
    nodes = findall("//x:math//x:bvar/x:ci", root(xml), ["x"=>mathml_ns])
    ivs = unique(nodecontent.(nodes))
    if length(ivs) > 1
        error("Only one independent variable (iv) is acceptable")
    end
    return ivs[1]
end

function list_params(xml, iv=find_iv(xml))
    vars = find_initialized_variables(xml)
    states = find_state_names(xml)
    p = filter(x -> x["name"] ∉ states && x["name"] != iv, vars)
    map(x -> (create_param(x["name"]) => parse(Float64, x["initial_value"])), p)
end

function list_initial_conditions(xml, iv=find_iv(xml))
    vars = find_initialized_variables(xml)
    states = find_state_names(xml)
    u0 = filter(x -> x["name"] ∈ states, vars)
    map(x -> (create_var(x["name"], iv) => parse(Float64, x["initial_value"])), u0)
end

function list_states(xml)
    states = find_state_names(xml)
    map(x -> (create_param(x)), states)
end

function list_substitution(xml, iv=find_iv(xml))
    # nodes = findall("//x:variable[@initial_value]", root(xml), ["x"=>cellml_ns])
    # vars = map(x -> x["name"], nodes)
    # states = find_state_names(xml)
    # Dict(create_var(x) => x ∈ states ? create_var(x, iv) : create_param(x) for x in vars)
    states = first.(list_initial_conditions(xml, iv))
    params = first.(list_params(xml, iv))
    u0_dict = Dict(create_var(operation(x).name) => x for x in states)
    p_dict = Dict(create_var(x.name) => x for x in params)
    union(u0_dict, p_dict)
end

function split_equations(eqs)
    alg = filter(eq -> eq.lhs isa Sym, eqs)
    eqs = filter(!(eq -> eq.lhs isa Sym), eqs)
    return eqs, alg
end

function flatten_equations(eqs)
    eqs, alg = split_equations(eqs)
    s = Dict(eq.lhs => eq.rhs for eq in alg)

    # alg = [eq.lhs ~ substitute(eq.rhs, s) for eq in alg]
    for (i,eq) in enumerate(alg)
        eq = eq.lhs ~ substitute(eq.rhs, s)
        s[eq.lhs] = eq.rhs
        alg[i] = eq
    end

    # s = Dict(a.lhs => a.rhs for a in alg)
    eqs = [eq.lhs ~ substitute(eq.rhs, s) for eq in eqs]

    return eqs
end

function process_cellml_xml(xml::EzXML.Document)
    ml = extract_mathml(xml)
    eqs = vcat(parse_node.(ml)...)
    eqs = flatten_equations(eqs)
    iv = find_iv(xml)
    s = list_substitution(xml, iv)
    eqs = [substitute(eq.lhs, s) ~ substitute(eq.rhs, s) for eq in eqs]
    sys = ODESystem(eqs, create_var(iv))
    return sys
end

function test_cellml(xml::EzXML.Document)
    ml = extract_mathml(xml)
    iv = find_iv(xml)

    for (i, m) in enumerate(ml)
        println(i, "  ****************************************")
        eqs = parse_node(m)
        for eq in eqs
            println(eq)
            #eq = flatten_equations(eq)
            s = list_substitution(xml, iv)
            eq = substitute(eq.lhs, s) ~ substitute(eq.rhs, s)
            try
                sys = ODESystem(eq)
            catch e
            end
        end
    end
end

# This is the fall-back option in case of breaking changes in
# ModelingToolkit and/or SymbolicUtils
function eval_string(xml::EzXML.Document)
    iv = find_iv(xml)
    p = list_params(xml, iv)
    u0 = list_initial_conditions(xml, iv)

    states = [string(operation(first(x)).name) for x in u0]
    params = [string(first(x).name) for x in p]

    s = "@parameters " * join(params, " ") * " " * iv * ";"
    s *= "iv_global = " * iv * ";"

    d = "(" * iv * ") "
    s *= "@variables " * join(states, d) * d * ";"

    ml = extract_mathml(xml)
    eqs = vcat(parse_node.(ml)...)
    s *= "eqs = " * string(eqs) * ";"

    s *= "p_global = [" * join(string.(p), ",") *"];"
    s *= "u0_global = [" * join([string(operation(first(x)).name) * " => " * string(last(x)) for x in u0], ", ") * "];"

    return s
end