text:
module BlockArrays
using Base.Cartesian
using LinearAlgebra, ArrayLayouts, FillArrays

# AbstractBlockArray interface exports
export AbstractBlockArray, AbstractBlockMatrix, AbstractBlockVector, AbstractBlockVecOrMat
export Block, getblock, getblock!, setblock!, eachblock, blocks
export blockaxes, blocksize, blocklength, blockcheckbounds, BlockBoundsError, BlockIndex
export blocksizes, blocklengths, blocklasts, blockfirsts, blockisequal
export BlockRange, blockedrange, BlockedUnitRange

export BlockArray, BlockMatrix, BlockVector, BlockVecOrMat, mortar
export PseudoBlockArray, PseudoBlockMatrix, PseudoBlockVector, PseudoBlockVecOrMat

export undef_blocks, undef, findblock, findblockindex

export khatri_rao, blockkron, BlockKron

export blockappend!, blockpush!, blockpushfirst!, blockpop!, blockpopfirst!

import Base: @propagate_inbounds, Array, to_indices, to_index,
            unsafe_indices, first, last, size, length, unsafe_length,
            unsafe_convert,
            getindex, ndims, show,
            step,
            broadcast, eltype, convert, similar,
            @_inline_meta, _maybetail, tail, @_propagate_inbounds_meta, reindex,
            RangeIndex, Int, Integer, Number,
            +, -, *, /, \, min, max, isless, in, copy, copyto!, axes, @deprecate,
            BroadcastStyle, checkbounds, throw_boundserror,
            ones, zeros, intersect, Slice
using Base: ReshapedArray, dataids

import Base: (:), IteratorSize, iterate, axes1, strides, isempty
import Base.Broadcast: broadcasted, DefaultArrayStyle, AbstractArrayStyle, Broadcasted, broadcastable
import LinearAlgebra: lmul!, rmul!, AbstractTriangular, HermOrSym, AdjOrTrans,
                        StructuredMatrixStyle
import ArrayLayouts: _fill_lmul!, MatMulVecAdd, MatMulMatAdd, MatLmulVec, MatLdivVec,
                        materialize!, MemoryLayout, sublayout, transposelayout, conjlayout,
                        triangularlayout, triangulardata, _inv, _copyto!, axes_print_matrix_row,
                        colsupport, rowsupport, sub_materialize

include("blockindices.jl")
include("blockaxis.jl")
include("abstractblockarray.jl")
include("blockarray.jl")
include("pseudo_blockarray.jl")
include("views.jl")
include("blocks.jl")

include("blockbroadcast.jl")
include("blocklinalg.jl")
include("blockproduct.jl")
include("show.jl")
include("blockreduce.jl")
include("blockdeque.jl")
include("blockarrayinterface.jl")

@deprecate getblock(A::AbstractBlockArray{T,N}, I::Vararg{Integer, N}) where {T,N} view(A, Block(I))
@deprecate getblock!(X, A::AbstractBlockArray{T,N}, I::Vararg{Integer, N}) where {T,N} copyto!(X, view(A, Block(I)))
@deprecate setblock!(A::AbstractBlockArray{T,N}, v, I::Vararg{Integer, N}) where {T,N} (A[Block(I...)] = v)

end # module


text:
function loss_seq(p, re, x, y)
  m = re(p)
  LTC.reset_state!(m, p) # use initial conditions from current params

  ŷb = Flux.Zygote.Buffer([y[1]], size(y,1))
  for i in 1:size(x,1)
    xi = x[i]
    ŷi = m(xi)
    Inf32 ∈ ŷi && return Inf32, copy(ŷb), y # TODO: what if a layer after MTKRecur can handle Infs?
    ŷb[i] = ŷi
  end
  ŷ = copy(ŷb)

  # mean(sum.(abs2, (ŷ .- y))), ŷ, y
  return mean(Flux.Losses.mse.(ŷ,y, agg=mean)), ŷ, y
end

function loss_seq(p, m::FastChain, x, y)
  # ŷ = m.(x, [p])

  LTC.reset_state!(m, p)

  ŷb = Flux.Zygote.Buffer([y[1]], size(y,1))
  for i in 1:size(x,1)
    xi = x[i]
    ŷi = m(xi, p)
    Inf32 ∈ ŷi && return Inf32, copy(ŷb), y # TODO: what if a layer after MTKRecur can handle Infs?
    ŷb[i] = ŷi
  end
  ŷ = copy(ŷb)

  # mean(sum.(abs2, (ŷ .- y))), ŷ, y
  return mean(Flux.Losses.mse.(ŷ,y, agg=mean)), ŷ, y
end

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    TaubinSmoothing(n; λ=0.5, μ=-0.5)

Perform `n` iterations of Taubin smoothing with parameters `λ` and `μ`.

## References

* Taubin, G. 1995. [Curve and Surface Smoothing without Shrinkage]
  (https://ieeexplore.ieee.org/document/466848)
"""
struct TaubinSmoothing <: SmoothingMethod
  n::Int
  λ::Float64
  μ::Float64
end

function TaubinSmoothing(n; λ=0.5, μ=-0.5)
  @assert n > 0 "invalid number of iterations"
  @assert 0 < λ ≤ -μ < 1 "invalid parameters λ and μ"
  TaubinSmoothing(n, λ, μ)
end

function smooth(mesh, method::TaubinSmoothing)
  n = method.n
  λ = method.λ
  μ = method.μ

  # Laplacian matrix with uniform weights
  L = laplacematrix(mesh, weights=:uniform)

  # matrix with vertex coordinates (nvertices x ndims)
  V = reduce(hcat, coordinates.(vertices(mesh))) |> transpose

  # Taubin updates
  for _ in 1:n
    V = V + λ*L*V
    V = V + μ*L*V
  end

  # new points of the smooth mesh
  points = Point.(eachrow(V))

  SimpleMesh(points, topology(mesh))
end


text:
"""
    SimpleWeightedGraph{T, U}

A type representing an undirected graph with weights of type `U`.
"""
mutable struct SimpleWeightedGraph{T<:Integer, U<:Real} <: AbstractSimpleWeightedGraph{T, U}
    weights::SparseMatrixCSC{U, T} # indexed by [dst, src]
end

ne(g::SimpleWeightedGraph) = nnz(g.weights) ÷ 2

# Graph{UInt8}(6), Graph{Int16}(7), Graph{UInt8}()
function (::Type{SimpleWeightedGraph{T, U}})(n::Integer = 0) where T<:Integer where U<:Real
    weights = spzeros(U, T, T(n), T(n))
    return SimpleWeightedGraph{T, U}(weights)
end

# Graph()
SimpleWeightedGraph() = SimpleWeightedGraph{Int, Float64}()

# Graph(6), Graph(0x5)
SimpleWeightedGraph(n::T) where T<:Integer = SimpleWeightedGraph{T, Float64}(n)

# Graph(UInt8)
SimpleWeightedGraph(::Type{T}) where T<:Integer = SimpleWeightedGraph{T, Float64}(zero(T))

# Graph(UInt8, Float32)
SimpleWeightedGraph(::Type{T}, ::Type{U}) where T<:Integer where U<:Real = SimpleWeightedGraph{U, T}(zero(T))

# Graph(SimpleGraph)
SimpleWeightedGraph(g::LightGraphs.SimpleGraphs.SimpleGraph{T}, ::Type{U}=Float64) where T <: Integer where U <: Real =
    SimpleWeightedGraph{T, U}(adjacency_matrix(g))

# Graph(SimpleDiGraph)
SimpleWeightedGraph(g::LightGraphs.SimpleGraphs.SimpleDiGraph{T}, ::Type{U}=Float64) where T <: Integer where U <: Real =
    SimpleWeightedGraph{T, U}(adjacency_matrix(LightGraphs.SimpleGraphs.SimpleGraph(g)))

# Graph(SimpleGraph, defaultweight)
SimpleWeightedGraph(g::LightGraphs.SimpleGraphs.SimpleGraph{T}, x::U) where T <: Integer where U <: Real =
    SimpleWeightedGraph{T, U}(x.*adjacency_matrix(g, U))

# Graph(SimpleDiGraph, defaultweight)
SimpleWeightedGraph(g::LightGraphs.SimpleGraphs.SimpleDiGraph{T}, x::U) where T <: Integer where U <: Real =
    SimpleWeightedGraph{T, U}(x.*adjacency_matrix(LightGraphs.SimpleGraphs.SimpleGraph(g), U))

# SimpleWeightedGraph{T, U}(SimpleGraph)
function (::Type{SimpleWeightedGraph{T, U}})(g::LightGraphs.SimpleGraphs.SimpleGraph)  where T<:Integer where U <: Real
    SimpleWeightedGraph{T, U}(adjacency_matrix(LightGraphs.SimpleGraphs.SimpleGraph{T}(g), U))
end

# DiGraph(srcs, dsts, weights)
function SimpleWeightedGraph(i::AbstractVector{T}, j::AbstractVector{T}, v::AbstractVector{U}; combine = +) where T<:Integer where U<:Real
    m = max(maximum(i), maximum(j))
    s = sparse(vcat(i,j), vcat(j,i), vcat(v,v), m, m, combine)
    SimpleWeightedGraph{T, U}(s)
end
# Graph{UInt8}(adjmx)
# function (::Type{SimpleWeightedGraph{T, U}})(adjmx::AbstractMatrix) where T<:Integer where U <: Real
#     dima,dimb = size(adjmx)
#     isequal(dima,dimb) || error("Adjacency / distance matrices must be square")
#     issymmetric(adjmx) || error("Adjacency / distance matrices must be symmetric")
#     g = SimpleWeightedGraph(U.(LinearAlgebra.fillstored!(copy(adjmx), 1)))
# end

# converts Graph{Int} to Graph{Int32}
# function (::Type{SimpleWeightedGraph{T, U}})(g::SimpleWeightedGraph) where T<:Integer where U<:Real
#     h_fadj = [Vector{T}(x) for x in fadj(g)]
#     return SimpleGraph(ne(g), h_fadj)
# end

# Graph(adjmx)
function SimpleWeightedGraph(adjmx::AbstractMatrix)
    dima,dimb = size(adjmx)
    isequal(dima,dimb) || error("Adjacency / distance matrices must be square")
    issymmetric(adjmx) || error("Adjacency / distance matrices must be symmetric")
    SimpleWeightedGraph{Int, eltype(adjmx)}(adjmx')
end

# Graph(digraph). Weights will be added.

SimpleWeightedGraph(g::SimpleWeightedDiGraph) = SimpleWeightedGraph(g.weights .+ g.weights')

edgetype(::SimpleWeightedGraph{T, U}) where T<:Integer where U<:Real= SimpleWeightedGraphEdge{T,U}

edges(g::SimpleWeightedGraph) = (SimpleWeightedEdge(x[1], x[2], x[3]) for x in zip(findnz(triu(g.weights))...))
weights(g::SimpleWeightedGraph) = g.weights
inneighbors(g::SimpleWeightedGraph, x...) = outneighbors(g, x...)

# add_edge! will overwrite weights.
function add_edge!(g::SimpleWeightedGraph, e::SimpleWeightedGraphEdge)
    @warn "Note: adding edges to this graph type is not performant." maxlog=1 _id=:swg_add_edge
    T = eltype(g)
    U = weighttype(g)
    s_, d_, w = Tuple(e)

    if w == zero(U)
        @warn "Note: adding edges with a zero weight to this graph type has no effect." maxlog=1 _id=:swg_add_edge_zero
        return false
    end

    s = T(s_)
    d = T(d_)
    (s in vertices(g) && d in vertices(g)) || return false
    g.weights[d, s] = w
    g.weights[s, d] = w
    return true
end

function rem_edge!(g::AbstractSimpleWeightedGraph, e::SimpleWeightedGraphEdge)
    has_edge(g, e) || return false
    U = weighttype(g)
    g.weights[dst(e), src(e)] = zero(U)
    g.weights[src(e), dst(e)] = zero(U)
    return true
end

==(g::SimpleWeightedGraph, h::SimpleWeightedGraph) = g.weights == h.weights

"""
    is_directed(g)

Return `true` if `g` is a directed graph.
"""
is_directed(::Type{SimpleWeightedGraph}) = false
is_directed(::Type{SimpleWeightedGraph{T, U}}) where T where U = false
is_directed(g::SimpleWeightedGraph) = false


text:
# Module-scope global variables

import Base: eltype

export
    CuGlobal, get, set

immutable CuGlobal{T}
    ptr::DevicePtr{Void}
    nbytes::Cssize_t

    function CuGlobal(mod::CuModule, name::String)
        ptr_ref = Ref{Ptr{Void}}()
        nbytes_ref = Ref{Cssize_t}()
        @apicall(:cuModuleGetGlobal, (Ptr{Ptr{Void}}, Ptr{Cssize_t}, CuModule_t, Ptr{Cchar}),
                                     ptr_ref, nbytes_ref, mod, name)
        if nbytes_ref[] != sizeof(T)
            throw(ArgumentError("size of global '$name' does not match type parameter type $T"))
        end
        @assert nbytes_ref[] == sizeof(T)
        new(unsafe_convert(DevicePtr{Void}, ptr_ref[]), nbytes_ref[])
    end
end

eltype{T}(::Type{CuGlobal{T}}) = T

function get{T}(var::CuGlobal{T})
    val_ref = Ref{T}()
    @apicall(:cuMemcpyDtoH, (Ptr{Void}, Ptr{Void}, Csize_t),
                            val_ref, var.ptr.inner, var.nbytes)
    return val_ref[]
end

function set{T}(var::CuGlobal{T}, val::T)
    val_ref = Ref{T}(val)
    @apicall(:cuMemcpyHtoD, (Ptr{Void}, Ptr{Void}, Csize_t),
                            var.ptr.inner, val_ref, var.nbytes)
end


text:
include("imports.jl")

include("QRM.jl")
include("utils.jl")


text:
test = """
    F10
    N3
    F7
    R90
    F11
    """

function parse_program(prog::AbstractString)::Vector
    return split(prog, "\n"; keepempty = false) .|> Instruction
end

struct Instruction{code}
    v::Int
    Instruction(cmd::AbstractString) = new{cmd[1]}(parse(Int, cmd[2:end]))
end

Base.show(io::IO, i::Instruction{code}) where code = print(io, "$(code)$(i.v)")

const E = [ 1,  0]
const N = [ 0, -1]
const W = [-1,  0]
const S = [ 0,  1]

struct State{mode}
    pos
    target
end

State{:ship}() = State{:ship}([0, 0], E)
State{:wp  }() = State{:wp  }([0, 0], 10E + N)

function run(program::AbstractVector; mode::Symbol)::Int
    state = foldl(|>, program; init = State{mode}())
    return state.pos .|> abs |> sum
end

(i::Instruction{'N'})(state::State) = move(state, N .* i.v)
(i::Instruction{'S'})(state::State) = move(state, S .* i.v)
(i::Instruction{'E'})(state::State) = move(state, E .* i.v)
(i::Instruction{'W'})(state::State) = move(state, W .* i.v)
(i::Instruction{'F'})(state::State) = forward(state, state.target .* i.v)
(i::Instruction{'L'})(state::State) = turn(state, i.v)
(i::Instruction{'R'})(state::State) = turn(state, -i.v)

move(state::State{:ship}, m) = forward(state, m)
move(state::State{:wp}, m) = State{:wp}(state.pos, state.target + m)
forward(state::S, m) where {S <: State} = S(state.pos + m, state.target)
turn(state::S, d) where {S <: State} = S(state.pos, rotate(state.target, d))
rotate(v, d::Number) = round.(Int, [cosd(d) sind(d); -sind(d) cosd(d)] * v)

q1(program::AbstractVector)::Int = run(program; mode = :ship)
q2(program::AbstractVector)::Int = run(program; mode = :wp)

let program = parse_program(test)
    @assert q1(program) == 25
    @assert q2(program) == 286
end

program = parse_program(read("day12.in", String))
println("Q1: ", q1(program))
println("Q2: ", q2(program))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
using AxisKeys
using MCMCChains
using ParetoSmooth
using Statistics
using Test
using CSV
using DataFrames
using StatsBase
using Turing

import RData

let og_array = RData.load("Example_Log_Likelihood_Array.RData")["x"]
    global log_lik_arr = copy(permutedims(og_array, [3, 1, 2]))
end
let og_weights = RData.load("Weight_Matrix.RData")["weightMatrix"]
    global r_weights = exp.(permutedims(reshape(og_weights, 500, 2, 32), [3, 1, 2]))
end
r_eff = RData.load("Rel_Eff.RData")["rel_eff"]
r_psis = RData.load("Psis_Object.RData")["x"]
r_tail_len = Int.(RData.load("Tail_Vector.RData")["tail"])
pareto_k = RData.load("Pareto_K.RData")["pareto_k"]
r_loo = RData.load("Example_Loo.RData")["example_loo"]

# Add labels, reformat
r_pointwise = KeyedArray(
    r_loo["pointwise"][:, Not(4)];
    data = 1:size(r_loo["pointwise"], 1),
    statistic=[:loo_est, :mcse, :overfit, :pareto_k],
)

r_loo["estimates"] = hcat(r_loo["estimates"], r_loo["estimates"] / size(r_pointwise, 1))
r_ests = KeyedArray(
    r_loo["estimates"][Not(3), :];
    criterion=[:loo_est, :overfit],
    statistic=[:total, :se_total, :mean, :se_mean],
)

@testset "ParetoSmooth.jl" begin
    @testset "Basic Arrays" begin

        # All of these should run
        with_r_eff = psis(log_lik_arr, r_eff)
        jul_psis = psis(log_lik_arr)
        log_lik_mat = reshape(log_lik_arr, 32, 1000)
        chain_index = vcat(fill(1, 500), fill(2, 500))
        matrix_psis = psis(log_lik_mat; chain_index=chain_index)
        log_psis = psis(log_lik_arr; log_weights=true)

        jul_loo = loo(log_lik_arr)
        r_eff_loo = psis_loo(log_lik_arr, r_eff)

        @test display(jul_psis) === nothing
        @test display(jul_loo) === nothing

        # max 20% difference in tail length calc between Julia and R
        @test maximum(abs.(log.(jul_psis.tail_len ./ r_tail_len))) ≤ .2
        @test maximum(abs.(jul_psis.tail_len .- r_tail_len)) ≤ 10
        @test maximum(abs.(with_r_eff.tail_len .- r_tail_len)) ≤ 2

        # RMSE from R version is less than .1%
        @test sqrt(mean((with_r_eff.weights ./ r_weights .- 1).^2)) ≤ .001
        # RMSE less than .2% when using InferenceDiagnostics' ESS
        @test sqrt(mean((jul_psis.weights ./ r_weights .- 1).^2)) ≤ .002
        # Max difference is 1%
        @test maximum(log_psis.weights .- log.(r_weights)) ≤ .01

        ## Test difference in loo pointwise results

        # Different r_eff
        jul_pointwise = jul_loo.pointwise([:loo_est, :mcse, :overfit, :pareto_k])
        errs = (r_pointwise - jul_pointwise).^2
        @test sqrt(mean(errs(:loo_est))) ≤ .01
        @test sqrt(mean(errs(:overfit))) ≤ .01
        @test sqrt(mean(errs(:pareto_k))) ≤ .025
        errs_mcse = log.(r_pointwise(:mcse) ./ jul_loo.pointwise(:mcse)).^2
        @test_broken sqrt(mean(errs_mcse)) ≤ .1

        # Same r_eff
        r_eff_pointwise = r_eff_loo.pointwise([:loo_est, :mcse, :overfit, :pareto_k])
        errs = (r_pointwise - r_eff_pointwise).^2
        @test sqrt(mean(errs(:loo_est))) ≤ .01
        @test sqrt(mean(errs(:overfit))) ≤ .01
        @test sqrt(mean(errs(:pareto_k))) ≤ .025
        errs_mcse = log.(r_pointwise(:mcse) ./ r_eff_loo.pointwise(:mcse)).^2
        @test_broken sqrt(mean(errs_mcse)) ≤ .1

        # Test estimates
        errs = r_ests - jul_loo.estimates(criterion=[:loo_est, :overfit])
        display(r_ests)
        display(jul_loo.estimates(criterion=[:loo_est, :overfit]))
        display(errs)
        @test maximum(abs.(errs)) ≤ .01

        errs = r_ests - r_eff_loo.estimates(criterion=[:loo_est, :overfit])
        display(errs)
        @test maximum(abs.(errs)) ≤ .01

        # Test for calling correct method
        @test jul_loo.psis_object.weights ≈ psis(-log_lik_arr).weights
        @test r_eff_loo.psis_object.weights ≈ psis(-log_lik_arr, r_eff).weights
    end

    @testset "MCMCChains and Turing utilities" begin
        using Distributions, Random
        Random.seed!(112)
        # simulated samples for μ
        samples = randn(100, 1, 1)
        data = randn(50)
        chain = Chains(samples)

        compute_loglike(μ, data) = logpdf(Normal(μ, 1), data)
        compute_loglike(μ, σ, data) = logpdf(Normal(μ, σ), data)

        pll1 = pointwise_log_likelihoods(compute_loglike, chain, data)
        pll2 = pointwise_log_likelihoods(compute_loglike, samples, data)
        # the pointwise log likehoods should be the same for both methods
        @test pll1 ≈ pll2 atol = 1e-6
        # test the dimensions: data points, samples, chains
        @test size(pll1) == (50, 100, 1)
        # test that sum of pointwise log likelihoods equals sum of log likelihoods
        @test sum(sum(map(s->logpdf.(Normal(s, 1), data), samples))) ≈ sum(pll1) atol = 1e-6
        # test that psis_loo works with MCMCChains and yields correct type
        psis_loo_output = psis_loo(compute_loglike, chain, data)
        @test isa(psis_loo_output, PsisLoo)
        # test that loo works with MCMCChains and yields correct type
        psis_output = loo(compute_loglike, chain, data)
        @test isa(psis_output, PsisLoo)
        # test that psis works with MCMCChains and yields correct type
        psis_output = psis(compute_loglike, chain, data)
        @test isa(psis_output, Psis)

        # ensure that methods work with r_eff argument
        r_eff = similar(pll2, 0)
        # test that psis_loo works with MCMCChains and yields correct type
        psis_loo_output = psis_loo(compute_loglike, chain, data, r_eff)
        @test isa(psis_loo_output, PsisLoo)
        # test that loo works with MCMCChains and yields correct type
        psis_output = loo(compute_loglike, chain, data, r_eff)
        @test isa(psis_output, PsisLoo)
        # test that psis works with MCMCChains and yields correct type
        psis_output = psis(compute_loglike, chain, data, r_eff)
        @test isa(psis_output, Psis)

        @model function model(data)
            μ ~ Normal(0, 1)
            σ ~ truncated(Cauchy(0, 1), 0, Inf)
            for i in eachindex(data)
                data[i] ~ Normal(μ, σ)
            end
        end

        chain = sample(model(data), NUTS(1000, .9), MCMCThreads(), 1000, 4)
        pw_lls_turing = pointwise_log_likelihoods(model(data), chain)
        pw_lls_loglike = pointwise_log_likelihoods(compute_loglike, chain, data)

        # test the dimensions: data points, samples, chains
        @test size(pw_lls_turing) == (50, 1000, 4)
        # test that sum of pointwise log likelihoods equals sum of log likelihoods
        turing_samples = Array(Chains(chain, :parameters).value)
        # make this more terse with @tullio or other method later
        LL = 0.0
        n_samples,n_parms,n_chains = size(turing_samples)
        for s in 1:n_samples
            for c in 1:n_chains
                LL += sum(logpdf.(Normal(turing_samples[s,:,c]...), data))
            end
        end
        @test LL ≈ sum(pw_lls_turing) atol = 1e-6
        # Turing should work the same as compute_loglike
        @test pw_lls_loglike ≈ pw_lls_turing atol = 1e-6

        # test that psis_loo works with Turing model and MCMCChains and yields correct type
        psis_loo_output = psis_loo(model(data), chain)
        @test isa(psis_loo_output, PsisLoo)
        # test that loo works with Turing model and MCMCChains and yields correct type
        psis_output = loo(model(data), chain)
        @test isa(psis_output, PsisLoo)
        # test that psis works with Turing model and MCMCChains and yields correct type
        psis_output = psis(model(data), chain)
        @test isa(psis_output, Psis)

        # ensure that methods work with r_eff argument
        r_eff = similar(pw_lls_turing, 0)
        # test that psis_loo works with Turin and gives correct type
        psis_loo_output = psis_loo(model(data), chain, r_eff)
        @test isa(psis_loo_output, PsisLoo)
        # test that loo works with Turing model and MCMCChains and gives correct type
        psis_output = loo(model(data), chain, r_eff)
        @test isa(psis_output, PsisLoo)
        # test that psis works with Turing model and MCMCChains and gives correct type
        psis_output = psis(model(data), chain, r_eff)
        @test isa(psis_output, Psis)

    end

    @testset "LooCompare" begin
        using Random
        Random.seed!(129111)

        df = CSV.read(joinpath("data", "WaffleDivorce.csv"), DataFrame)
        df.D = zscore(df.Divorce)
        df.M = zscore(df.Marriage)
        df.A = zscore(df.MedianAgeMarriage)
        data = (D=df.D, A=df.A)

        function lin(a, b, c, x...)
            result = @. a + b * c
            for i in 1:2:length(x)
                @. result += x[i] * x[i+1]
            end
            return result
        end

        @model function m5_1t(A, D)
            a ~ Normal(0, 0.2)
            bA ~ Normal(0, 0.5)
            σ ~ Exponential(1)
            for i in eachindex(D)
                μ = lin(a, A[i], bA)
                D[i] ~ Normal(μ, σ)
            end
        end

        chn5_1t = sample(m5_1t(df.A, df.D), NUTS(1000, .9), MCMCThreads(), 1000, 4)

        @model function m5_2t(M, D)
            a ~ Normal(0, 0.2)
            bM ~ Normal(0, 0.5)
            σ ~ Exponential(1)
            for i in eachindex(D)
                μ = lin(a, M[i], bM)
                D[i] ~ Normal(μ, σ)
            end
        end

        chn5_2t = sample(m5_2t(df.M, df.D), NUTS(1000, .9), MCMCThreads(), 1000, 4)

        @model function m5_3t(A, M, D)
            a ~ Normal(0, 0.2)
            bA ~ Normal(0, 0.5)
            bM ~ Normal(0, 0.5)
            σ ~ Exponential(1)
            for i in eachindex(D)
                μ = a + M[i] * bM + A[i] * bA
                D[i] ~ Normal(μ, σ)
            end
        end

        chn5_3t = sample(m5_3t(df.A, df.M, df.D), NUTS(1000, .9), MCMCThreads(), 1000, 4)
        pw_lls5_1t = pointwise_log_likelihoods(m5_1t(df.A, df.D), chn5_1t)
        psis_loo_output5_1t = psis_loo(m5_1t(df.A, df.D), chn5_1t)
        psis_loo_output5_1t |> display

        pw_lls5_2t = pointwise_log_likelihoods(m5_2t(df.M, df.D), chn5_2t)
        psis_loo_output5_2t = psis_loo(m5_2t(df.M, df.D), chn5_2t)
        psis_loo_output5_2t |> display

        pw_lls5_3t = pointwise_log_likelihoods(m5_3t(df.A, df.M, df.D), chn5_3t)
        psis_loo_output5_3t = psis_loo(m5_3t(df.A, df.M, df.D), chn5_3t)
        psis_loo_output5_3t |> display

        loo = loo_compare([pw_lls5_1t, pw_lls5_2t, pw_lls5_3t];
            model_names=[:m5_1t, :m5_2t, :m5_3t])
        loo |> display

        @test loo.table(:m5_1t, :elpd_diff) ≈ 0.00 atol = 0.05
        @test loo.table(:m5_1t, :se_diff) ≈ 0.00 atol = 0.05
        @test loo.table(:m5_1t, :weight) ≈ 0.67 atol = 0.05
        @test loo.table(:m5_2t, :elpd_diff) ≈ -6.68 atol = 0.05
        @test loo.table(:m5_2t, :se_diff) ≈ 4.74 atol = 0.05
        @test loo.table(:m5_2t, :weight) ≈ 0.00 atol = 0.05
        @test loo.table(:m5_3t, :elpd_diff) ≈ -0.69 atol = 0.05
        @test loo.table(:m5_3t, :se_diff) ≈ 0.42 atol = 0.05
        @test loo.table(:m5_3t, :weight) ≈ 0.33 atol = 0.05

    end

end


text:
using Logging

struct AbacoLogger <: AbstractLogger
    groups::Vector{Symbol}
    modules::Vector{Module}
end

function logging(;debug=[])
    groups = [item for item in debug if isa(item, Symbol)]
    modules = [item for item in debug if isa(item, Module)]
    AbacoLogger(groups, modules) |> global_logger
end

function Logging.min_enabled_level(logger::AbacoLogger)
    Logging.Debug
end

function Logging.shouldlog(logger::AbacoLogger, level, _module, group, id)
    level >= Logging.Info || group in logger.groups || _module in logger.modules
end

function Logging.handle_message(logger::AbacoLogger, level, message, _module, group, id, file, line; kwargs...)
    println("[$(now())][$_module][$(Threads.threadid())][$level] $message")
end


text:
println("RODE Tests")
using Flux , StochasticDiffEq, DiffEqNoiseProcess , Optim, Test
println("RODE_tests")
using  NeuralPDE

using Random
Random.seed!(100)

println("Test Case 1")
linear = (u,p,t,W) ->   2u*sin(W)
tspan = (0.00f0, 1.00f0)
u0 = 1.0f0
dt = 1/50f0
W = WienerProcess(0.0,0.0,nothing)
prob = RODEProblem(linear, u0 ,tspan, noise = W)
chain = Flux.Chain(Dense(2,8,relu),Dense(8,16 , relu),Dense(16,1))
opt = ADAM(1e-4)
sol = solve(prob, NeuralPDE.NNRODE(chain,W,opt), dt=dt, verbose = true,
            abstol=1e-10, maxiters = 3000)
W2 = NoiseWrapper(sol.W)
prob1 = RODEProblem(linear , u0 , tspan , noise = W2)
sol2 = solve(prob1,RandomEM(),dt=dt)
err = Flux.mse(sol.u , sol2.u)
@test err < 0.3

println("Test Case 2")
linear = (u,p,t,W) -> t^3 + 2*t + (t^2)*((1+3*(t^2))/(1+t+(t^3))) - u*(t + ((1+3*(t^2))/(1+t+t^3))) + 5*W
tspan = (0.00f0, 1.00f0)
u0 = 1.0f0
dt = 1/100f0
W = WienerProcess(0.0,0.0,nothing)
prob = RODEProblem(linear, u0 ,tspan, noise = W)
chain = Flux.Chain(Dense(2,32,sigmoid),  Dense(32,32,sigmoid)  , Dense(32,1))
opt = ADAM(1e-3)
sol = solve(prob, NeuralPDE.NNRODE(chain,W,opt), dt=dt, verbose = true,
            abstol=1e-10, maxiters = 2000)
W2 = NoiseWrapper(sol.W)
prob1 = RODEProblem(linear , u0 , tspan , noise = W2)
sol2 = solve(prob1,RandomEM(),dt=dt)
err = Flux.mse(sol.u , sol2.u)
@test err < 0.4


text:
@testset "Parser" begin

sx = SExpressions.parse

@testset "numbers" begin
    @test sx("-1") == -1
    @test sx("+1000") == 1000
    @test sx("-15/5") == -3
    @test sx("+13/17") == 13//17
    @test sx("-1/2") == -1//2
    @test sx("-1/2+3i") == -1//2+3im
    @test sx("3+4/3i") == 3+4im//3
    @test sx("3-2i") == 3-2im
end

@testset "booleans" begin
    @test sx("#t")
    @test sx("#true")
    @test sx("#T")
    @test !sx("#f")
    @test !sx("#false")
    @test !sx("#F")
end

@testset "arithmetic" begin
    @test sx("(+ 1 1)") == List(:(+), 1, 1)
    @test sx("(+ 1/2 1/3)") == List(:(+), 1//2, 1//3)
    @test sx("""
    (define (sqr x) (^ x 2))
    """) == lispify((
            :define,
            (:sqr, :x),
            (:(^), :x, 2)))
end

@testset "invalid syntax" begin
    @test_throws ErrorException SExpressions.parseall("(+ 1 1")
    @test_throws ErrorException SExpressions.parseall("(+ 1 1]")
    @test_throws ErrorException SExpressions.parseall("(+ 1 1))")
    @test_throws ErrorException SExpressions.parseall("(+ 1 1;)")
end

@testset "comment" begin
    @test sx("1;") == 1
    @test sx(";
             x") == :x
    @test sx(";
    ;
    ;
    4;2
    ") == 4
end

@testset "strings" begin
    @test sx("\"x\"") == "x"
    @test sx("\"x\\n\"") == "x\n"
    @test sx("""\"
    x\\
    y\\\\
    z\"""") == "\nxy\\\nz"
end

@testset "symbols" begin
    @test sx("foo") == :foo
    @test_broken sx("|x|") == :x
    @test_broken sx("| |") == Symbol(" ")
end

@testset "lists" begin
    @test sx("(1 2 3)") == List(1, 2, 3)
    @test sx("(1 . 2)") == Cons(1, 2)
    @test sx("(1 . (2))") == List(1, 2)
    @test sx("(1 . + . 2)") == List(:+, 1, 2)
    @test sx("(1 2 . 3)") == Cons(1, Cons(2, 3))
    @test_throws ErrorException sx("(.)")
    @test_throws ErrorException sx("(. x)")
    @test_throws ErrorException sx("(1 . 2 3)")
    @test_throws ErrorException sx("(1 . 2 . 3 . 4)")
end

@testset "reader macros" begin
    @test sx("'x") == List(:quote, :x)
    @test sx("`(+ 1 2)") == List(:quasiquote, List(:(+), 1, 2))
    @test sx("`(+ ,x 2)") == List(:quasiquote,
                                  List(:(+), List(:unquote, :x), 2))
    @test_broken sx("`(+ ,@xs 2)") ==
        lispify((:quasiquote,
                 (:(+),
                  (Symbol("unquote-splicing"), :x),
                  2)))
end

@testset "file" begin
    @test SExpressions.parsefile("data/scheme.scm") isa List
end

@testset "macro" begin
    # we try to limit these in number, because if they fail then a LoadError
    # occurs, and that is bad for debugging
    @test sx"1" == 1
    @test sx"foo" == :foo
end

end


text:
using Documenter
using JsonGrinder

DocMeta.setdocmeta!(JsonGrinder, :DocTestSetup, :(using JsonGrinder); recursive=true)
# for running only doctests
# doctest(JsonGrinder)
makedocs(
         sitename = "JsonGrinder.jl",
         # doctest = false,
         format = Documenter.HTML(sidebar_sitename=false,
                                  prettyurls=get(ENV, "CI", nothing) == "true",
                                  assets=["assets/favicon.ico", "assets/custom.css"]),
         modules = [JsonGrinder],
         pages = ["Home" => "index.md",
                  "Schema" => "schema.md",
                  "Creating extractors" => "extractors.md",
                  "Extractor functions" => "exfunctions.md",
                  "AutoML" => "automl.md",
                  "API Documentation" => "api.md",
                  "Developers" => "developers.md",
                  "Citation" => "citation.md"
                  ],

)

deploydocs(
    repo = "github.com/pevnak/JsonGrinder.jl.git",
)


text:
@testset "Masking" begin
    X = ones(3, 3)
    out = mask_circle(X, 1)
    @test out[2, 2] === NaN
    @test all(i == 5 ? true : out[i] !== NaN for i in eachindex(out))
    @test all(mask_circle(X, 2) .=== NaN)

    # interface
    A = mask_circle(ones(100, 100), 5)
    B = mask_circle!(ones(100, 100), 5)
    @test filter(!isnan, A) == filter(!isnan, B)
    @test_throws MethodError mask_circle(ones(3, 100, 100), 5)
    @test_throws MethodError mask_circle(ones(Int, 100, 100), 5)

end


text:
function _manipulate_outer(expr, is_gui)
  @assert expr.head == :for

  cur_block = strip_escape!(expr.args[2])
  cur_bindings = strip_escape!(expr.args[1])

  is_gui && insert!(
    cur_block.args, 2, :( SimplePlots.SimplePlot() )
  )

  if cur_bindings.head == :block
    cur_bindings = cur_bindings.args
  else
    cur_bindings = [cur_bindings]
  end

  cur_symbols = map(x -> first(x.args), cur_bindings)

  gui_id = UUIDs.uuid4()
  quote_widgets = map(make_gui_widget, cur_bindings)

  return quote_widgets, cur_block, cur_bindings, cur_symbols
end

function _manipulate_inner(cur_widgets, cur_listener, is_gui)
  cur_id = UUIDs.uuid4()

  make_gui_html(cur_id, cur_widgets)

  comm_observer = Observable(0)
  comm_id = "interact-$( string(cur_id) )"

  parent_message = IJulia.execute_msg

  cur_watcher = on(comm_observer) do cur_value

    IJulia.set_cur_msg(parent_message)

    cur_comm = Comm(Symbol(comm_id))

    cur_comm.on_msg = function (cur_message)

      IJulia.set_cur_msg(parent_message)

      message_data = cur_message.content["data"]

      for (cur_key, cur_value) in message_data
        ( cur_key == "___interact_plot_id___" ) && continue

        cur_widget = cur_widgets[findfirst(tmp_widget -> tmp_widget.label == cur_key, cur_widgets)]

        if cur_widget.datatype <: AbstractString
          parsed_value = cur_value
        elseif cur_widget.datatype <: Symbol
          parsed_value = Symbol(cur_value)
        elseif cur_widget.datatype == Char
          if isa(cur_value, Char)
            parsed_value = cur_value
          else
            @assert isa(cur_value, AbstractString)
            @assert length(cur_value) == 1
            parsed_value = cur_value[1]
          end
        else
          parsed_value = parse(cur_widget.datatype, cur_value)
        end

        observe!(cur_widget, parsed_value)
      end

      if isa(cur_listener[], AbstractPlot)
        shown_plot = cur_listener[]
      else
        if !isnothing(cur_listener[])
          if isa(cur_listener[], AbstractString) || isa(cur_listener[], Number)
            println(cur_listener[])
          else
            display(cur_listener[])
          end
        end

        is_gui || return
        shown_plot = _plot
      end

      if "___interact_plot_id___" in keys(message_data)
        plot_json = custom_json(shown_plot)
        cur_size, cur_data, cur_layout, cur_config = _show_helper(shown_plot)

        # cur_size is currently unused in update call

        plot_script = _render_html_script(
          message_data["___interact_plot_id___"], cur_data, cur_layout, cur_config
        )

        send_comm(cur_comm, Dict("json" => Dict("plot" => plot_json, "script" => plot_script)))
      else
        plot_html = _show(shown_plot).content
        send_comm(cur_comm, Dict("html" => plot_html))
      end

    end

  end

  interact_comms[comm_id] = comm_observer
  make_gui_bootloader(Main, cur_id)

  cur_id
end


text:
function make_gif(window::PlotWindow, width, height, duration, filename="plot.gif", delay=0.06)
	SFML.make_gif(window.renderwindow)
end

export make_gif


text:
include("prop_logic_theory.jl")
include("prover.jl")

ex = rewrite(:(((p => q) ∧ (r => s) ∧ (p ∨ r)) => (q ∨ s)), impl)
prove(t, ex, 1, 25)
@profview prove(t, ex, 2, 7)


text:
function g()
  local x::AbstractString
  x = 2
end

function main()
  println("This message will appear")
  try
    g()
  catch
    println("Type error is caught at runtime!")
  end
end

main()

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
include("drawing.jl")
include("font.jl")
include("text_drawing.jl")
include("assets.jl")
include("layers.jl")
include("auto_tiler.jl")
include("decals.jl")
include("entities.jl")
include("triggers.jl")
include("fillers.jl")

function drawObjectTile(ctx::Cairo.CairoContext, x::Integer, y::Integer, tile::Integer; alpha::Number=getGlobalAlpha())
    scenery = getSprite("tilesets/scenery", "Gameplay")

    width = floor(Int, scenery.realWidth / 8)

    drawX, drawY = (x - 1) * 8, (y - 1) * 8
    quadY, quadX = divrem(tile, width)

    drawImage(ctx, "tilesets/scenery", drawX, drawY, quadX * 8, quadY * 8, 8, 8, alpha=alpha)
end

tileDrawPosition(x::Integer, y::Integer) = (x - 1) * 8, (y - 1) * 8

function getTileData(x::Integer, y::Integer, tiles::Tiles, meta::TilerMeta, states::TileStates)
    tileValue = tiles.data[y, x]

    if haskey(meta.paths, tileValue)
        imagePath = meta.paths[tileValue]
        quads, sprite = getMaskQuads(x, y, tiles, meta)

        return quads[Integer(mod1(states.rands[y, x], length(quads)))], sprite, imagePath
    end

    return nothing, nothing, nothing
end

function tileNeedsUpdate(tileValue::Char, x::Integer, y::Integer, quad::Tuple{Integer, Integer}, states::TileStates)
    return tileValue != states.chars[y, x] || states.quads[y, x] != quad
end

function drawTile(ctx::Cairo.CairoContext, x::Integer, y::Integer, tiles::Tiles, meta::TilerMeta, states::TileStates; alpha::Number=getGlobalAlpha())
    tileValue = tiles.data[y, x]
    drawX, drawY = tileDrawPosition(x, y)

    if tileValue != '0'
        quad, sprite, imagePath = getTileData(x, y, tiles, meta, states)

        if quad !== nothing
            quadX, quadY = quad

            if !isempty(sprite)
                animatedMeta = filter(m -> m.name == sprite, animatedTilesMeta)[1]
                frames = findTextureAnimations(animatedMeta.path, getAtlas("Gameplay"))

                if !isempty(frames)
                    frame = frames[Integer(mod1(states.rands[y, x], length(frames)))]
                    frameSprite = getSprite(frame, "Gameplay")

                    ox = animatedMeta.posX - frameSprite.offsetX
                    oy = animatedMeta.posY - frameSprite.offsetY

                    # TODO - What do we actually have to clear here?
                    #clearArea(ctx, drawX + ox, drawY + oy, 8, 8)
                    drawImage(ctx, frame, drawX + ox, drawY + oy, alpha=alpha)
                end
            end

            if tileNeedsUpdate(tileValue, x, y, quad, states)
                clearArea(ctx, drawX, drawY, 8, 8)
                drawImage(ctx, imagePath, drawX, drawY, quadX * 8, quadY * 8, 8, 8, alpha=alpha)

                states.quads[y, x] = quad
            end

        else
            drawRectangle(ctx, drawX, drawY, 8, 8, colors.unknown_tile_color, (0.0, 0.0, 0.0, 0.0))
        end

    else
        if states.chars[y, x] != '0'
            clearArea(ctx, drawX, drawY, 8, 8)
        end
    end

    states.chars[y, x] = tileValue
end

function drawTiles(ctx::Cairo.CairoContext, tiles::Tiles, objtiles::ObjectTiles, meta::TilerMeta, states::TileStates, width::Integer, height::Integer; alpha::Number=getGlobalAlpha(), fg::Bool=true, useObjectTiles::Bool=false)
    for y in 1:height, x in 1:width
        objtile = get(objtiles.data, (y, x), -1)
        tile = get(tiles.data, (y, x), '0')

        if useObjectTiles && objtile != -1 && tile != '0'
            drawObjectTile(ctx, x, y, objtile, alpha=alpha)

        else
            drawTile(ctx, x, y, tiles, meta, states, alpha=alpha)
        end
    end
end

function drawTiles(ctx::Cairo.CairoContext, dr::DrawableRoom, fg::Bool=true; alpha::Number=getGlobalAlpha(), useObjectTiles::Bool=false)
    objtiles = dr.room.objTiles
    tiles = fg ? dr.room.fgTiles : dr.room.bgTiles
    height, width = size(tiles.data)

    states = fg ? dr.fgTileStates : dr.bgTileStates
    meta = fg ? fgTilerMeta : bgTilerMeta

    if size(states) != (height, width)
        updateTileStates!(dr.room, dr.map.package, states, width, height, fg)
    end

    drawTiles(ctx, tiles, objtiles, meta, states, width, height, alpha=alpha, fg=fg, useObjectTiles=useObjectTiles)

    return true
end

drawTiles(layer::Layer, dr::DrawableRoom, fg::Bool=true; alpha::Number=getGlobalAlpha(), useObjectTiles::Bool=false) = drawTiles(getSurfaceContext(layer.surface), dr, fg, alpha=alpha, useObjectTiles=useObjectTiles)

function drawParallax(ctx::Cairo.CairoContext, parallax::Maple.Parallax, camera::Camera, fg::Bool=true)

end

function drawApply(ctx::Cairo.CairoContext, apply::Maple.Apply, camera::Camera, fg::Bool=true)

end

function drawEffect(ctx::Cairo.CairoContext, effect::Maple.Effect, camera::Camera, fg::Bool=true)

end

backgroundFuncs = Dict{Type, Function}(
    Maple.Parallax => drawParallax,
    Maple.Apply => drawApply,
    Maple.Effect => drawEffect,
)

function drawBackground(layer::Layer, dr::DrawableRoom, camera::Camera, fg::Bool=true)
    backdrops = fg ? dr.map.style.foregrounds : dr.map.style.backgrounds
    ctx = getSurfaceContext(layer.surface)

    if !fg
        color = something(dr.fillColor, getRoomBackgroundColor(dr.room))
        paintSurface(ctx, color)
    end

    for backdrop in backdrops
        t = typeof(backdrop)

        if haskey(backgroundFuncs, t)
            backgroundFuncs[t](ctx, backdrop, camera, fg)
        end
    end
end

function drawDecals(layer::Layer, dr::DrawableRoom, fg::Bool=true)
    ctx = getSurfaceContext(layer.surface)
    decals = fg ? dr.room.fgDecals : dr.room.bgDecals

    for decal in decals
        drawDecal(ctx, decal)
    end

    return true
end

function drawEntities(layer::Layer, dr::DrawableRoom)
    ctx = getSurfaceContext(layer.surface)
    entities = dr.room.entities

    for entity in entities
        renderEntity(ctx, layer, entity, dr.room)
    end

    return true
end

function drawTriggers(layer::Layer, dr::DrawableRoom)
    ctx = getSurfaceContext(layer.surface)
    triggers = dr.room.triggers

    for trigger in triggers
        renderTrigger(ctx, layer, trigger, dr.room)
    end

    return true
end

function resetDrawingSettings()
    reset!(camera)
end

loadExternalSprites!()

drawingLayers = nothing
drawableRooms = Dict{Map, Dict{Room, DrawableRoom}}()

redrawingFuncs["fgDecals"] = (layer, room) -> drawDecals(layer, room, true)
redrawingFuncs["bgDecals"] = (layer, room) -> drawDecals(layer, room, false)

redrawingFuncs["fgTiles"] = (layer, room) -> drawTiles(layer, room, true)
redrawingFuncs["bgTiles"] = (layer, room) -> drawTiles(layer, room, false)

redrawingFuncs["fgParallax"] = (layer, room, camera) -> drawBackground(layer, room, camera, true)
redrawingFuncs["bgParallax"] = (layer, room, camera) -> drawBackground(layer, room, camera, false)

redrawingFuncs["entities"] = drawEntities
redrawingFuncs["triggers"] = drawTriggers

function getDrawableRooms(map::Map)
    if !haskey(drawableRooms, map)
        drawableRooms[map] = Dict{Room, DrawableRoom}()
    end

    return collect(values(drawableRooms[map]))
end

function getDrawableRoom(map::Map, room::Room)
    if !haskey(drawableRooms, map)
        drawableRooms[map] = Dict{Room, DrawableRoom}()
    end

    if !haskey(drawableRooms[map], room)
        drawableRooms[map][room] = DrawableRoom(map, room)
    end

    return drawableRooms[map][room]
end

function deleteDrawableRoomCache(map::Map)
    rooms = getDrawableRooms(map)
    for room in rooms
        destroy(room)
    end

    delete!(drawableRooms, map)
end

function updateDrawingLayers!(map::Map, room::Room)
    global drawingLayers = getDrawableRoom(map, room).layers
end

function drawMap(ctx::Cairo.CairoContext, camera::Camera, map::Map; adjacentAlpha::Number=colors.adjacent_room_alpha, backgroundFill::colorTupleType=colors.background_canvas_fill)
    width, height = floor(Int, ctx.surface.width), floor(Int, ctx.surface.height)

    paintSurface(ctx, backgroundFill)
    pattern_set_filter(get_source(ctx), Cairo.FILTER_NEAREST)

    for filler in map.fillers
        if fillerVisible(camera, width, height, filler)
            drawFiller(ctx, camera, filler)
        end
    end

    for room in map.rooms
        if roomVisible(camera, width, height, room)
            dr = getDrawableRoom(map, room)
            alpha = loadedState.roomName == room.name ? 1.0 : adjacentAlpha

            drawRoom(ctx, camera, dr, alpha=alpha)

        else
            if get(config, "delete_non_visible_rooms_cache", false)
                if haskey(drawableRooms, map)
                    if haskey(drawableRooms[map], room)
                        destroy(drawableRooms[map][room])
                        delete!(drawableRooms[map], room)
                    end
                end
            end
        end
    end
end

drawMap(canvas::Gtk.GtkCanvas, camera::Camera, map::Map) = drawMap(Gtk.getgc(canvas), camera, map)

function drawRoom(ctx::Cairo.CairoContext, camera::Camera, room::DrawableRoom; alpha::Number=getGlobalAlpha())
    renderingLayer = room.rendering
    drawingLayers = room.layers

    renderingLayer.redraw = redrawRenderingLayer(renderingLayer, drawingLayers)

    if renderingLayer.redraw
        resetLayer!(renderingLayer, room)
        renderingCtx = getSurfaceContext(renderingLayer.surface)

        renderingLayer.redraw = false
        combineLayers!(renderingCtx, drawingLayers, camera, room)
    end

    Cairo.save(ctx)

    translate(ctx, -camera.x, -camera.y)
    scale(ctx, camera.scale, camera.scale)
    translate(ctx, room.room.position...)

    applyLayer!(ctx, renderingLayer, alpha=alpha)

    restore(ctx)
end

max_stars_repo_name: zyowo/Ahorn

max_issues_repo_name: zyowo/Ahorn

max_forks_repo_name: zyowo/Ahorn


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

function eponesweepT0!(epfields::EPFields, epalg::EPAlg, epmatT0::EPMatT0, stat = Dict())

    @extract epfields : av va a d μ s siteflagave siteflagvar
    @extract epalg : alpha beta_vec minvar maxvar epsconv damp
    @extract epmatT0 : Σd Σi G lb ub vd vi Y

    M = size(G,1)
    N = length(av)

    idxy = 1:M # dependent variables
    idxw = M+1:N # independent variables

    ad,ai = view(a,idxy), view(a,idxw)     # dep an ind prior mean (epfields)
    dd,di = view(d,idxy), view(d,idxw)     # dep an ind prior variance (epfields)
    sd,si = view(s,idxy), view(s,idxw)     # dep and ind untruncated marginals variance (epfields)
    μd,μi = view(μ,idxy), view(μ,idxw)     # dep and ind untruncated marginals mean (epfields)
    avd,avi = view(av,idxy), view(av,idxw) # dep and ind truncated marginals mean (epfields)
    vad,vai = view(va,idxy), view(va,idxw) # dep and ind truncated marginals variance (epfields)
    βd, βi = view(beta_vec, idxy), view(beta_vec, idxw)
    errav, errva, errμ, errs = fill(typemin(eltype(av)), 4)
    Gt = G'

    # Test
    # @show sum(beta_vec)

    # All fields in epmat are updated from the epfields of last sweep
    # (?) covariance matrix of independent variables (epmat)
    stat[:elapsed_eponesweep_inv] = @elapsed begin
        Di = Diagonal(inv.(di))
        Dd = Diagonal(inv.(dd))
        # Σi = inv(Σᵢ⁻¹)
        Σᵢ⁻¹ = Gt * Dd * G + Di
        inplaceinverse!(Σi, Σᵢ⁻¹)
    end
    # fast_similarity_inv!(Σi, di, dd, G)
    mul!(Σd, G*Σi, Gt) # (?) covariance matrix of dependent variables (epmat)
    # Original ep
    # mul!(vi,Σi, ai ./ di - G'*(ad ./ dd)) # (?) mean vector of independent variables (epmat)
    # ep-maxent
    vi .= Σi * (Gt * Dd * (Y - ad) + Di * ai - Gt * βd + βi)

    # vd = -Gvi + b'
    mul!(vd, G, vi) # (?) mean vector of dependent variables (epmat)
    vd .= Y - vd

    for i in eachindex(μi)  # loop M+1:N
        newμi,newsi = newμs(Σi[i,i], ai[i], di[i], vi[i], lb[i+M], ub[i+M], minvar, maxvar)
        errμ = max(errμ, abs(μi[i]-newμi))
        errs = max(errs, abs(si[i]-newsi))
        μi[i] = newμi
        si[i] = newsi
        # println("μi[$(i+M)] = ", μi[i]," si[$(i+M)] = ", si[i], " Σi = ",Σi[i,i] )

        newavw,newvaw = newav(si[i], μi[i], avi[i], vai[i], siteflagave[i+M], siteflagvar[i+M],
                              lb[i+M], ub[i+M], minvar, maxvar)
        errav = max(errav,abs(avi[i]-newavw))
        errva = max(errva,abs(vai[i]-newvaw))
        avi[i] = newavw
        vai[i] = newvaw

        newai,newdi = matchmom(μi[i], si[i], avi[i], vai[i], minvar, maxvar)
        ai[i] = damp * ai[i] + (1.0-damp) * newai # modify a in epfields
        di[i] = damp * di[i] + (1.0-damp) * newdi # modify d in epfields
    end

    for i in eachindex(μd)   # loop  1:M

        newμd, newsd = newμs(Σd[i,i], ad[i], dd[i], vd[i], lb[i], ub[i], minvar, maxvar)
        errμ = max(errμ, abs(μd[i]-newμd))
        errs = max(errs, abs(sd[i]-newsd))
        μd[i] = newμd # modify μ in epfields
        sd[i] = newsd # modify s in epfields
#        println("μd[$i] = ", μd[i]," sd[$i] = ", sd[i], " Σ = ", Σd[i,i], " (",lb[i],":",ub[i],")"," ad[$i] = ",ad[i], " dd[$i] = ", dd[i])

        newavy,newvay = newav(sd[i],μd[i],avd[i],vad[i],
            siteflagave[i], siteflagvar[i],
            lb[i], ub[i], minvar, maxvar)
        errav = max(errav, abs(avd[i]-newavy))
        errva = max(errva, abs(vad[i]-newvay))
        avd[i] = newavy # modify av in epfields
        vad[i] = newvay # modify va in epfields

        newad,newbd = matchmom(μd[i],sd[i],avd[i],vad[i],minvar,maxvar)
        ad[i] = damp * ad[i] + (1.0-damp) * newad # modify a in epfields
        dd[i] = damp * dd[i] + (1.0-damp) * newbd # modify d in epfields
    end
    return errav, errva, errμ, errs
end


text:
"The bedGraph record."
struct Record{T<:Real}
    chrom::String
    first::Int
    last::Int
    value::T

    Record{T}(chrom::String, first::Int, last::Int, value::T) where {T} = new{T}(chrom, first, last, value)
    Record{T}(chrom, first, last, value) where {T} = new{T}(_string(chrom), _int(first), _int(last), _real(value))
end

function Record(chrom, first, last, value)
    v = _real(value)
    return Record{typeof(v)}(chrom, first, last, v)
end

function Record(data) :: Record
    return data
end

function Record{T}(data) :: Record{T} where T
    return data
end

function Base.convert(::Type{NamedTuple{names}}, record::Record) where {names}
    return NamedTuple{names}((
        record.chrom,
        record.first,
        record.last,
        record.value,
    ))
end

Base.convert(::Type{NamedTuple}, record::Record) = convert(NamedTuple{(:chrom, :first, :last, :value)}, record)

function Base.convert(::Type{Record}, nt::NamedTuple{(:chrom, :first, :last, :value)})
    return Record(nt.chrom, nt.first, nt.last, nt.value)
end

function Base.convert(::Type{Record}, str::AbstractString)
    data = _split_line(str)
    return Record(data[1], data[2], data[3], data[4])
end

function Base.convert(::Type{Record{T}}, str::AbstractString) where T
    data = _split_line(str)
    return Record{T}(data[1], data[2], data[3], data[4])
end

function Base.convert(::Type{Vector{String}}, record::Record)
    return String[record.chrom, string(record.first), string(record.last), string(record.value)]
end

function Base.convert(::Type{String}, record::Record, delim::String = " ")

    strs = [
        record.chrom,
        string(record.first),
        string(record.last),
        string(record.value)
    ]

    return join(strs, delim)
end

"Access [`Record`](@ref)'s chrom field."
function chrom(record::Record)
    return record.chrom
end

"Access [`Record`](@ref)'s left position."
function Base.first(record::Record)
    return record.first
end

"Access [`Record`](@ref)'s last position."
function Base.last(record::Record)
    return record.last
end

"Access [`Record`](@ref)'s value."
function value(record::Record)
    return record.value
end

function Base.:(==)(a::Record, b::Record)
    return a.chrom == b.chrom &&
           a.first == b.first &&
           a.last  == b.last &&
           a.value == b.value
end

function Base.isless(a::Record, b::Record, chrom_isless::Function=isless)

    if chrom_isless(a.chrom, b.chrom)
        return true
    end

    if a.first < b.first
        return true
    end

    if a.last < b.last
        return true
    end

    if a.value < b.value
        return true
    end

    return false
end

## Internal helper functions.
function _split_line(line::AbstractString)
    return split(line, r"\s", limit=5, keepempty=false) #Note: may return 5 elements.
end

function _int(pos)
    return pos
end

function _string(str)
    return str
end

function _real(x)
    return x
end

function _int(str::AbstractString)
    return parse(Float64, str) #Note: conversion to Int is completed when the Record is constructed.
end

function _real(str::AbstractString)
    return parse(Float64, str)
end


text:
"""
    struct SkeletonPair{L,R} <: GridapType
      plus::L
      minus::R
    end
"""
struct SkeletonPair{L,R} <: GridapType
  plus::L
  minus::R
end

function Base.getproperty(x::SkeletonPair, sym::Symbol)
  if sym == :⁺
    x.plus
  elseif sym == :⁻
    x.minus
  else
    getfield(x, sym)
  end
end

function Base.propertynames(x::SkeletonPair, private::Bool=false)
  (fieldnames(typeof(x))...,:⁺,:⁻)
end

"""
    struct SkeletonTriangulation{Dc,Dp,B} <: Triangulation{Dc,Dp}
      plus::B
      minus::B
    end
"""
struct SkeletonTriangulation{Dc,Dp,B} <: Triangulation{Dc,Dp}
  plus::B
  minus::B
  function SkeletonTriangulation(plus::B,minus::B) where B<:Triangulation
    Dc = num_cell_dims(plus)
    Dp = num_point_dims(plus)
    #@assert Dc + 1 == Dp
    new{Dc,Dp,B}(plus,minus)
  end
end

function Skeleton(args...;kwargs...)
  SkeletonTriangulation(args...;kwargs...)
end

function Base.getproperty(x::SkeletonTriangulation, sym::Symbol)
  if sym == :⁺
    x.plus
  elseif sym == :⁻
    x.minus
  else
    getfield(x, sym)
  end
end

function Base.propertynames(x::SkeletonTriangulation, private::Bool=false)
  (fieldnames(typeof(x))...,:⁺,:⁻)
end

# Triangulation interface

get_background_model(t::SkeletonTriangulation) = get_background_model(t.plus)
get_grid(t::SkeletonTriangulation) = get_grid(t.plus)
get_glue(t::SkeletonTriangulation{D},::Val{D}) where D = get_glue(t.plus,Val(D))

function get_glue(trian::SkeletonTriangulation,::Val{Dp}) where Dp
  model = get_background_model(trian)
  Dm = num_cell_dims(model)
  get_glue(trian,Val(Dp),Val(Dm))
end

function get_glue(trian::SkeletonTriangulation,::Val{Dp},::Val{Dm}) where {Dp,Dm}
  nothing
end

function get_glue(trian::SkeletonTriangulation,::Val{D},::Val{D}) where D
  plus = get_glue(trian.plus,Val(D))
  minus = get_glue(trian.minus,Val(D))
  SkeletonPair(plus,minus)
end

function is_change_possible(sglue::FaceToFaceGlue,tglue::SkeletonPair)
  is_change_possible(sglue,tglue.plus) && is_change_possible(sglue,tglue.minus)
end

function get_facet_normal(trian::SkeletonTriangulation)
  plus = get_facet_normal(trian.plus)
  minus = get_facet_normal(trian.minus)
  SkeletonPair(plus,minus)
end

# Related with CompositeTriangulation
function _compose_glues(rglue::FaceToFaceGlue,dglue::SkeletonPair)
  plus = _compose_glues(rglue,dglue.plus)
  minus = _compose_glues(rglue,dglue.minus)
  SkeletonPair(plus,minus)
end

function Base.view(glue::SkeletonPair{<:FaceToFaceGlue},ids::AbstractArray)
  SkeletonPair(view(glue.plus,ids),view(glue.minus,ids))
end

function restrict(a::SkeletonPair,b::AbstractArray)
  SkeletonPair(restrict(a.plus,b),restrict(a.minus,b))
end

# Constructors

"""
    SkeletonTriangulation(model::DiscreteModel,face_to_mask::Vector{Bool})
    SkeletonTriangulation(model::DiscreteModel)
"""
function SkeletonTriangulation(model::DiscreteModel,face_to_mask::AbstractVector{Bool})
  left_cell_around = 1
  plus = BoundaryTriangulation(model,face_to_mask,left_cell_around)
  right_cell_around = 2
  minus = BoundaryTriangulation(model,face_to_mask,right_cell_around)
  SkeletonTriangulation(plus,minus)
end

function SkeletonTriangulation(model::DiscreteModel)
  topo = get_grid_topology(model)
  D = num_cell_dims(model)
  face_to_mask = collect(Bool, .!get_isboundary_face(topo,D-1))
  SkeletonTriangulation(model,face_to_mask)
end

function SkeletonTriangulation(rtrian::Triangulation,args...;kwargs...)
  rmodel = get_active_model(rtrian)
  dtrian = SkeletonTriangulation(rmodel,args...;kwargs...)
  CompositeTriangulation(rtrian,dtrian)
end

function SkeletonTriangulation(
  model::DiscreteModel,
  reffe::ReferenceFE,
  face_own_dofs::Vector{Vector{Int}})

  ncells = num_cells(model)
  cell_to_reffe = Fill(1,ncells)
  SkeletonTriangulation(model,[reffe],cell_to_reffe,[face_own_dofs])
end

function SkeletonTriangulation(
  model::DiscreteModel,
  reffes::Vector{<:ReferenceFE},
  cell_to_reffe::AbstractVector,
  reffe_to_face_own_dofs::Vector{Vector{Vector{Int}}})

  reffe_to_polytope = map(get_polytope,reffes)
  reffe_to_faces = map(get_faces,reffe_to_polytope)
  reffe_to_range = map(p->get_dimrange(p,num_dims(p)-1),reffe_to_polytope)

  reffe_to_lfacet_to_mask = _compute_reffe_to_facet_to_mask(
    reffe_to_faces,reffe_to_range,reffe_to_face_own_dofs)

  topo = get_grid_topology(model)

  D = num_cell_dims(model)
  nfacets = num_facets(model)
  facet_to_isboundary = get_isboundary_face(topo,D-1)
  cell_to_facets = get_faces(topo,D,D-1)
  facet_to_mask = _compute_disc_facet_mask(
    facet_to_isboundary,cell_to_facets,cell_to_reffe,reffe_to_lfacet_to_mask)

  SkeletonTriangulation(model,facet_to_mask)

end

function _compute_disc_facet_mask(
  facet_to_isboundary, cell_to_facets::Table, cell_to_reffe, reffe_to_lfacet_to_mask)
  nfacets = length(facet_to_isboundary)
  facet_to_mask = Vector{Bool}(undef,nfacets)
  ncells = length(cell_to_reffe)
  for cell in 1:ncells
    pini = cell_to_facets.ptrs[cell]
    pend = cell_to_facets.ptrs[cell+1]-1
    reffe = cell_to_reffe[cell]
    lfacet_to_mask = reffe_to_lfacet_to_mask[reffe]
    for (lfacet,p) in enumerate(pini:pend)
      facet = cell_to_facets.data[p]
      mask = lfacet_to_mask[lfacet]
      isinterior = ! facet_to_isboundary[facet]
      facet_to_mask[facet] = mask && isinterior
    end
  end
  facet_to_mask
end

function _compute_reffe_to_facet_to_mask(
  reffe_to_faces,reffe_to_range,reffe_to_face_own_dofs)

  reffe_to_facet_to_mask = Vector{Bool}[]
  nreffes = length(reffe_to_faces)
  for reffe in 1:nreffes
    range = reffe_to_range[reffe]
    facet_to_faces = reffe_to_faces[reffe][range]
    facet_to_own_dofs = reffe_to_face_own_dofs[reffe][range]
    face_to_own_dofs = reffe_to_face_own_dofs[reffe]
    nfacets = length(facet_to_faces)
    facet_to_mask = fill(false,nfacets)
    for facet in 1:nfacets
      n = length(facet_to_own_dofs[facet])
      faces = facet_to_faces[facet]
      for face in faces
        n += length(face_to_own_dofs[face])
      end
      facet_to_mask[facet] = n == 0
    end
    push!(reffe_to_facet_to_mask,facet_to_mask)
  end

  reffe_to_facet_to_mask
end

const IN = -1
const OUT = 1

"""
"""
function InterfaceTriangulation(model::DiscreteModel,cell_to_is_in::Vector{Bool})
  cell_to_inout = fill(Int8(OUT),length(cell_to_is_in))
  cell_to_inout[cell_to_is_in] .= IN
  InterfaceTriangulation(model,cell_to_inout)
end

function Interface(args...;kwargs...)
  InterfaceTriangulation(args...;kwargs...)
end

function InterfaceTriangulation(model::DiscreteModel,cells_in,cells_out)
  cell_to_inout = fill(Int8(0),num_cells(model))
  cell_to_inout[cells_in] .= IN
  cell_to_inout[cells_out] .= OUT
  InterfaceTriangulation(model,cell_to_inout)
end

function InterfaceTriangulation(trian_in::Triangulation,trian_out::Triangulation)
  D = num_cell_dims(trian_in)
  @assert D == num_cell_dims(trian_out)
  glue_in = get_glue(trian_in,Val(D))
  glue_out = get_glue(trian_out,Val(D))
  cells_in = glue_in.tface_to_mface
  cells_out = glue_out.tface_to_mface
  model = get_background_model(trian_in)
  @check model === get_background_model(trian_out)
  @notimplementedif D != num_cell_dims(model) "Not implemented, but it should be easy to implement."
  InterfaceTriangulation(model,cells_in,cells_out)
end

function InterfaceTriangulation(model::DiscreteModel,cell_to_inout::AbstractVector{<:Integer})

  D = num_cell_dims(model)
  facet_grid = Grid(ReferenceFE{D-1},model)
  cell_grid = Grid(ReferenceFE{D},model)

  topo = get_grid_topology(model)
  facet_to_cells = Table(get_faces(topo,D-1,D))
  ifacet_to_facet, facet_to_lcell_left, facet_to_lcell_right = _find_interface_facets(
    cell_to_inout, facet_to_cells)

  ifacet_grid = view(facet_grid,ifacet_to_facet)

  glue_left = FaceToCellGlue(topo,cell_grid,ifacet_grid,ifacet_to_facet,facet_to_lcell_left)
  glue_right = FaceToCellGlue(topo,cell_grid,ifacet_grid,ifacet_to_facet,facet_to_lcell_right)

  trian = BodyFittedTriangulation(model,ifacet_grid,ifacet_to_facet)
  plus = BoundaryTriangulation(trian,glue_left)
  minus = BoundaryTriangulation(trian,glue_right)

  SkeletonTriangulation(plus,minus)
end

function _find_interface_facets( cell_to_inout, facet_to_cells::Table)

  nifacets = 0
  for facet in 1:length(facet_to_cells)
    a = facet_to_cells.ptrs[facet]
    b = facet_to_cells.ptrs[facet+1]
    if b-a == 2
      cell1 = facet_to_cells.data[a]
      cell2 = facet_to_cells.data[a+1]
      inout_1 = cell_to_inout[cell1]
      inout_2 = cell_to_inout[cell2]
      if (inout_1== IN && inout_2==OUT) || (inout_1== OUT && inout_2==IN)
        nifacets += 1
      end
    end
  end

  T = eltype(eltype(facet_to_cells))
  ifacet_to_facet = zeros(T,nifacets)
  nfacets = length(facet_to_cells)
  facet_to_lcell_left = fill(Int8(1),nfacets)
  facet_to_lcell_right = fill(Int8(2),nfacets)

  nifacets = 0
  for facet in 1:length(facet_to_cells)
    a = facet_to_cells.ptrs[facet]
    b = facet_to_cells.ptrs[facet+1]
    if b-a == 2
      cell1 = facet_to_cells.data[a]
      cell2 = facet_to_cells.data[a+1]
      inout_1 = cell_to_inout[cell1]
      inout_2 = cell_to_inout[cell2]
      if (inout_1== IN && inout_2==OUT) || (inout_1== OUT && inout_2==IN)
        nifacets += 1
        ifacet_to_facet[nifacets] = facet
        if inout_1 == IN
          facet_to_lcell_left[facet] = 1
          facet_to_lcell_right[facet] = 2
        else
          facet_to_lcell_left[facet] = 2
          facet_to_lcell_right[facet] = 1
        end
      end
    end
  end

   ifacet_to_facet, facet_to_lcell_left, facet_to_lcell_right
 end


text:
include_model("hopper_impedance")

# Simulator
function dynamics(model::HopperImpedance, v1, q1, q2, q3, u, λ, b, w, h)
   (M_func(model, q1) * v1
   - M_func(model, q2) * (SVector{4}(q3) - SVector{4}(q2)) / h
   + h * (transpose(B_func(model, q3)) * SVector{2}(u)
   + transpose(N_func(model, q3)) * SVector{1}(λ)
   + transpose(P_func(model, q3)) * SVector{2}(b)
   - C_func(model, q3, (q3 - q2) / h)))
end

function maximum_dissipation(model::HopperImpedance, q2, q3, ψ, η, h)
   ψ_stack = ψ[1] * ones(model.nb)
   return P_func(model, q3) * (q3 - q2) / h + ψ_stack - η
end

function no_slip(model::HopperImpedance, q2, q3, λ, h)
   return (λ' * _P_func(model, q3) * (q3 - q2) / h)[1]
end

function friction_cone(model::HopperImpedance, λ, b)
   return @SVector [model.μ * λ[1] - sum(b)]
end

mutable struct MOISimulator <: MOI.AbstractNLPEvaluator
   num_var::Int                 # number of decision variables
   num_con::Int                 # number of constraints
   primal_bounds
   constraint_bounds

   model
   slack_penalty
   v1
   q1
   q2
   u
   w
   h
end

function simulator_problem(model, v1, q1, q2, u, w, h; slack_penalty = 1.0e5)
   num_var = model.nq + model.nc + model.nb + model.nc + model.nb + model.ns
   num_con = model.nq + model.nc + model.nc + model.nb +  3

   zl = zeros(num_var)
   zl[1:model.nq] = model.qL
   zu = Inf * ones(num_var)
   zu[1:model.nq] = model.qU

   cl = zeros(num_con)
   cu = zeros(num_con)
   cu[model.nq .+ (1:model.nc)] .= Inf
   cu[model.nq + model.nc .+ (1:model.nc)] .= Inf
   cu[model.nq + model.nc + model.nc + model.nb + 1] = Inf
   cu[model.nq + model.nc + model.nc + model.nb + 2] = Inf
   cu[model.nq + model.nc + model.nc + model.nb + 3] = Inf

   MOISimulator(num_var,
	   num_con,
	   (zl, zu),
	   (cl, cu),
	   model,
	   slack_penalty,
	   v1,
	   q1,
	   q2,
	   u,
	   w,
	   h)
end
# prob_sim = simulator_problem(model, ones(nq), ones(nq), ones(nq), zeros(nu), ones(nq), h)
# z0 = rand(prob_sim.num_var)

function MOI.eval_objective(prob::MOISimulator, x)
   return prob.slack_penalty * x[prob.num_var]
end

# MOI.eval_objective(prob_sim, z0)

function MOI.eval_objective_gradient(prob::MOISimulator, grad_f, x)
   grad_f .= 0.0
   grad_f[prob.num_var] = prob.slack_penalty
   return nothing
end

# ∇obj = zeros(prob_sim.num_var)
# MOI.eval_objective_gradient(prob_sim, ∇obj, z0)

function MOI.eval_constraint(prob::MOISimulator, g, x)
   model = prob.model
   q3 = view(x, 1:model.nq)
   λ = view(x, model.nq .+ (1:model.nc))
   b = view(x, model.nq + model.nc .+ (1:model.nb))
   ψ = view(x, model.nq + model.nc + model.nb .+ (1:model.nc))
   η = view(x, model.nq + model.nc + model.nb + model.nc .+ (1:model.nb))
   s = view(x, model.nq + model.nc + model.nb + model.nc + model.nb .+ (1:model.ns))
   # λ = 0
   # b = 0
   g[1:model.nq] = dynamics(model,
	   prob.v1, prob.q1, prob.q2, q3, prob.u, λ, b, prob.w, prob.h)
   g[model.nq .+ (1:model.nc)] = ϕ_func(model, q3)
   g[model.nq + model.nc .+ (1:model.nc)] = friction_cone(model, λ, b)
   g[model.nq + model.nc + model.nc .+ (1:model.nb)] = maximum_dissipation(model,
	   prob.q2, q3, ψ, η, prob.h)
   g[model.nq + model.nc + model.nc + model.nb + 1] = s[1] - (λ' * ϕ_func(model, q3))[1]
   g[model.nq + model.nc + model.nc + model.nb + 2] = s[1] - (ψ' * friction_cone(model, λ, b))[1]
   g[model.nq + model.nc + model.nc + model.nb + 3] = s[1] - (η' * b)[1]

   return nothing
end
# prob_sim.h
# c0 = zeros(prob_sim.num_con)
# MOI.eval_constraint(prob_sim, c0, z0)

function MOI.eval_constraint_jacobian(prob::MOISimulator, jac, x)
   con!(g, z) = MOI.eval_constraint(prob, g, z)
   ForwardDiff.jacobian!(reshape(jac, prob.num_var, prob.num_con), con!, zeros(prob.num_con), x)
   return nothing
end

# ∇c0 = vec(zeros(prob_sim.num_con, prob_sim.num_var))
# MOI.eval_constraint_jacobian(prob_sim, ∇c0, z0)

function sparsity_jacobian(prob::MOISimulator)
   row = []
   col = []

   r = (1:prob.num_con)
   c = (1:prob.num_var)

   row_col!(row, col, r, c)

   return collect(zip(row, col))
end

# sparsity_jacobian(prob_sim)
#
# @time solve(prob_sim, z0)

function step_contact(model, v1, q1, q2, u, w, h)
   prob = simulator_problem(model, v1, q1, q2, u, w, h)
   z0 = [copy(q2); 1.0e-5 * ones(model.nc + model.nb + model.nc + model.nb + model.ns)]
   @time z , info = solve(prob, copy(z0), tol = 1.0e-7, c_tol = 1.0e-7, mapl = 0)

   @assert z[end] < 1.0e-5

   return z[1:model.nq]
end

model = HopperImpedance{Discrete, FixedTime}(n, m, d,
			   mb, Jb, mf,
			   μ, 150.0, r0, g,
			   qL, qU,
			   uL, uU,
			   nq,
		       nu,
		       nc,
		       nf,
		       nb,
		   	   ns,
		       idx_u,
		       idx_λ,
		       idx_b,
		       idx_ψ,
		       idx_η,
		       idx_s)
tf = 5.0
T = 251
h = tf / (T - 1)
t_sim = range(0.0, stop = tf, length = T)

q1 = zeros(model.nq)
q1[2] = 1.0
# q1[3] = pi / 100.0
q1[4] = model.r0
q_sim = [q1, q1]
v_sim = [(q_sim[end] - q_sim[end-1]) / h]

for t = 1:T-1
	println("t: $(t_sim[t])")
	q = step_contact(model,
		v_sim[end], q_sim[end-1], q_sim[end], zeros(model.nu), zeros(model.d), h)
	v = (q - q_sim[end]) / h

	push!(q_sim, q)
	push!(v_sim, v)
end

plot(hcat(q_sim...)[1:4, :]',
	label = ["x" "z" "t" "r"],
	legend = :bottomleft)

ϕ = [ϕ_func(model, q) for q in q_sim]
plot(t_sim[1:length(ϕ)-1], hcat(ϕ[2:end]...)',
	label = "sdf")

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, q_sim, Δt = h)


text:
include("Parser.jl")

using DataStructures

struct StaticAnalysisException <: Exception
  msg::String
end

default_value_strings = Dict{MPType,String}(
  MInt => "-1",
  MReal => "-1.0",
  MBool => "false",
  MString => "",
)

initial_scope_level = 0

unary_result_types = Dict{Tuple{TokenClass,MPType},MPType}(
  (kw_not, MBool) => MBool,
  (minus, MInt) => MInt,
  (minus, MReal) => MReal
)

array_type_to_scalar_type = Dict{MPType,MPType}(
  MIntArray => MInt,
  MRealArray => MReal,
  MBoolArray => MBool,
  MStringArray => MString
)

ref_type_to_scalar_type = Dict{MPType,MPType}(
  MIntRef => MInt,
  MRealRef => MReal,
  MBoolRef => MBool,
  MStringRef => MString,
  MInt => MInt,
  MReal => MReal,
  MBool => MBool,
  MString => MString,
)

binary_result_types = DefaultDict(MError,
  (times, MInt, MInt) => MInt,
  (times, MReal, MReal) => MReal,
  (plus, MInt, MInt) => MInt,
  (plus, MReal, MReal) => MReal,
  (plus, MString, MString) => MString,
  (minus, MInt, MInt) => MInt,
  (minus, MReal, MReal) => MReal,
  (divide, MInt, MInt) => MInt,
  (divide, MReal, MReal) => MReal,
  (modulo, MInt, MInt) => MInt,
  (kw_or, MBool, MBool) => MBool,
  (kw_and, MBool, MBool) => MBool,
  (equals, MInt, MInt) => MBool,
  (not_equal, MInt, MInt) => MBool,
  (less_than, MInt, MInt) => MBool,
  (less_than_or_equal, MInt, MInt) => MBool,
  (greater_than, MInt, MInt) => MBool,
  (greater_than_or_equal, MInt, MInt) => MBool,
  (kw_and, MBool, MBool) => MBool,
  (equals, MReal, MReal) => MBool,
  (not_equal, MReal, MReal) => MBool,
  (less_than, MReal, MReal) => MBool,
  (less_than_or_equal, MReal, MReal) => MBool,
  (greater_than, MReal, MReal) => MBool,
  (greater_than_or_equal, MReal, MReal) => MBool,
  (equals, MString, MString) => MBool,
  (not_equal, MString, MString) => MBool,
  (less_than, MString, MString) => MBool,
  (less_than_or_equal, MString, MString) => MBool,
  (greater_than, MString, MString) => MBool,
  (greater_than_or_equal, MString, MString) => MBool,
  (equals, MBool, MBool) => MBool,
  (not_equal, MBool, MBool) => MBool,
  (less_than, MBool, MBool) => MBool,
  (less_than_or_equal, MBool, MBool) => MBool,
  (greater_than, MBool, MBool) => MBool,
  (greater_than_or_equal, MBool, MBool) => MBool,
)

mutable struct Argument
  name::String
  is_var::Bool
  type::MPType
  size::Int
end

mutable struct SubroutineEntry
  name::String
  param_names::Vector{String}
  param_types::Vector{MPType}
  return_type::MPType
  node::Union{Subroutine,Nothing}
end

mutable struct AnalysisContext
  scope_stack::Stack{SymTableEntry}
  scope::Int
  subroutines::Vector{SubroutineEntry}
  in_call::Stack{Call}
  unique_calls::Vector{Call}
  all_var_names_and_types::Vector{Tuple{String,MPType}}
  var_id_counter::Int
end
AnalysisContext() = AnalysisContext(
      Stack{SymTableEntry}(),
      initial_scope_level,
      Vector{SubroutineEntry}(),
      Stack{Call}(),
      Vector{Call}(),
      Vector{Tuple{String,MPType}}(),
      0
      )

function create_var_id(ac::AnalysisContext, var_type::MPType, var_name::String)
  counter = ac.var_id_counter
  ac.var_id_counter += 1
  return "%_$(var_type)_$(var_name)_$(counter)"
end

function hasvariable(s::AnalysisContext, var_name::String, inner_scope_only::Bool=false)
  if DEBUG
    foreach(println, s.scope_stack)
  end
  if isempty(s.scope_stack) return false end
  if inner_scope_only
    scope = s.scope
    return any(entry->entry.scope_level == scope && entry.var_name == var_name, s.scope_stack)
  end
  return any(entry->entry.var_name == var_name, s.scope_stack)
end

function getvariable(s::AnalysisContext, var_name)
  return getfirst(entry -> entry.var_name == var_name, s.scope_stack)
end

function popscope!(s::AnalysisContext)
  isempty(s.scope_stack) && return
  s.scope -= 1
  while !isempty(s.scope_stack) && (first(s.scope_stack).scope_level > s.scope)
    pop!(s.scope_stack)
  end
end

function enterscope!(s::AnalysisContext)
  s.scope += 1
end

function pushvar!(ac::AnalysisContext, var_name::String, var_type::MPType, val_identifier::String, is_implicit::Bool=false)
  entry = SymTableEntry(var_name, var_type, val_identifier, ac.scope, is_implicit)
  push!(ac.scope_stack, entry)
end

function pushsubroutine!(ac::AnalysisContext,
    name::String,
    param_names::Vector{String},
    param_types::Vector{MPType},
    return_type::MPType,
    node::Union{Subroutine,Nothing})
  entry = SubroutineEntry(name, param_names, param_types, return_type, node)
  push!(ac.subroutines, entry)
end

function hassubroutine(ac::AnalysisContext, name::String)
  return any(entry->entry.name == name, ac.subroutines)
end

function getsubroutine(srs::Vector{SubroutineEntry}, name::String)
  index = findfirst(sr->sr.name == name, srs)
  return srs[index]
end

function getsubroutine(ac::AnalysisContext, name::String)
  getsubroutine(ac.subroutines, name)
end

# The entry point for the static analyzer.
function static_analysis(AST::Program, ac::AnalysisContext)

  # Push default values for all var, type combinations used in program
  push_defaults!(ac)

  # Process definitions
  static_analysis(AST.definitions, ac)

  # Process main block
  static_analysis(AST.main, ac)
end

function push_defaults!(ac::AnalysisContext)
  for (name, mptype) in ac.all_var_names_and_types
    if mptype ∈ scalar_types
      ref_type = scalar_to_ref_type[mptype]
      pushvar!(ac, name, ref_type, "@default_$(mptype)")
    end
  end
  enterscope!(ac)
end

static_analysis(AST::Program) = static_analysis(AST, AnalysisContext())

function static_analysis(d::Definitions, ac::AnalysisContext)
  DEBUG && println("This is static analysis, analysing definitions")
  if !allunique([def.name for def in d.defs])
    throw(StaticAnalysisException(
      "All functions and subroutines must have unique names. Overloading is not allowed."
    ))
  end
  for subroutine::Subroutine in d.defs
    get_signature(subroutine, ac)
  end
end

function get_signature(s::Subroutine, ac::AnalysisContext)
  foreach(param::Parameter->static_analysis(param, ac), s.params)
  true_types = Vector{MPType}([parameter_to_mptype(param) for param in s.params])
  param_names = Vector{String}([param.name for param in s.params])
  static_analysis(s.ret_type, ac)
  if !isa(s.ret_type.size, ImmediateInt)
    throw(StaticAnalysisException(
      "Reference types (strings and arrays) are not allowed as function return types (line $(s.line))."
    ))
  end
  pushsubroutine!(ac,
    s.name,
    param_names,
    true_types,
    s.ret_type.scalar_type,
    s
  )
end

function parameter_to_mptype(p::Parameter)
  if p.is_array && p.is_var_par
    throw(StaticAnalysisException(
      "Array type cannot be used as var parameter (line $(p.line))."
    ))
  end
  if p.is_var_par
    return scalar_to_ref_type[p.scalar_type]
  end
  if p.is_array
    return scalar_to_array_type[p.scalar_type]
  end
  return p.scalar_type
end

function parameter_to_llvmtype(p::Parameter)
  mptype_to_llvm_type[parameter_to_mptype(p)]
end

function typecheck(b::BinaryOperation, ac::AnalysisContext)
  typecheck(b.left, ac)
  typecheck(b.right, ac)
  ret_type = binary_result_types[(b.op, b.left.type, b.right.type)]
  if ret_type == MError
    throw(StaticAnalysisException(
      "Type mismatch ($(b.left.type) and $(b.right.type)) for binary operation \"$(b.op)\" on line $(b.line)."
    ))
  end
  b.type = ret_type
end

function typecheck(u::UnaryOperation, ac::AnalysisContext)
  typecheck(u.operand, ac)
  ret_type = unary_result_types[(u.op, u.operand.type)]
  if ret_type == MError
    throw(StaticAnalysisException(
      "Invalid type $(u.operand.type) for unary operation \"$(u.op)\" (line $(u.line))."
    ))
  end
  u.type = ret_type
end

function static_analysis(s::Subroutine, ac::AnalysisContext)
  static_analysis(s.body, ac)
  s.type = MPassed
end

function static_analysis(c::CallStatement, ac::AnalysisContext)
  DEBUG && println("This is static analysis, analysing CallStatement")
  analyze_call(c, ac)
  c.type = MPassed
end

function analyze_call(c::Call, ac::AnalysisContext)
  DEBUG && println("This is analyze_call, call is to $(c.identifier)")
  DEBUG && println("implicit_params are $(c.implicit_params)")

  subroutine_entry::SubroutineEntry = getsubroutine(ac, c.identifier)
  for (i, (param_type, arg)) in enumerate(zip(subroutine_entry.param_types, c.arguments))
    if param_type ∈ ref_types  # Parameter is a var parameter
      c.arguments[i] = GetRef(c.arguments[i], c.line)
    end
  end

  CallFactor

  foreach(arg->typecheck(arg, ac), c.arguments)
  compare_with_signature(subroutine_entry, c.arguments, c.line)

  c.subroutine = subroutine_entry.node::Subroutine
  c.subroutine.is_called = true

  all_var_names = Vector{String}(unique([name for (name, type) in ac.all_var_names_and_types]))

  c.implicit_params = [getvariable(ac, p) for p in all_var_names]

  subroutines_in_call = [call.identifier for call in ac.in_call]
  if subroutine_entry.name ∈ subroutines_in_call
    return
  end

  llvm_function_name = get_llvm_function_name(subroutine_entry.name, c.call_id)
  c.subroutine.llvm_function_name = llvm_function_name
  enterscope!(ac)
  c.scope_level = ac.scope
  push!(ac.in_call, c)
  for var_entry::SymTableEntry in c.implicit_params
    param_name = var_entry.var_name
    param_type = var_entry.var_type
    name_as_implicit = get_implicit_param_id(param_name, param_type)
    ref_type = param_type ∈ ref_types ? param_type : scalar_to_ref_type[param_type]
    pushvar!(ac, param_name, ref_type, "%$(name_as_implicit)", true)
  end
  for (param_name, arg) in zip(subroutine_entry.param_names, c.arguments)
    type = arg isa GetRef ? arg.type : ref_type_to_scalar_type[arg.type]
    pushvar!(ac, param_name, arg.type, "%$(param_name)")
  end
  static_analysis(c.subroutine, ac)
  pop!(ac.in_call)
  popscope!(ac)
end

function get_implicit_param_id(name::String, mptype::MPType)
  "implicit.$(ref_type_to_scalar_type[mptype]).$(name)"
end

function compare_with_signature(sr::SubroutineEntry, args::Vector{Value}, line::Int)
  if length(sr.param_types) != length(args)
    throw(StaticAnalysisException(
      "Subroutine call has $(length(args)) arguments, expected $(length(sr.param_types)) (line $line)."
    ))
  end
  for (param_type, arg) in zip(sr.param_types, args)
    if arg.type != param_type
      throw(StaticAnalysisException(
        "Expected argument of type $param_type, got $(arg.type) (line $line)."
    ))
    end
  end
end

function static_analysis(r::Return, ac::AnalysisContext)
  current_subroutine::Subroutine = first(ac.in_call).subroutine
  typecheck(r.value, ac)
  if r.value.type != current_subroutine.ret_type.true_type
    throw(StaticAnalysisException(
      "Expected return type to be $(current_subroutine.ret_type.true_type), got $(r.value.type) (line $(r.line))."
    ))
  end
  r.type = MPassed
end

function static_analysis(p::Parameter, ac::AnalysisContext)
  typecheck(p.size, ac)
  if p.size.type != MInt
    throw(StaticAnalysisException(
      "Array size must have integer value, got $(p.size.type) (line $(p.line))."
    ))
  end
  p.type = MPassed
end

function static_analysis(b::Block, ac::AnalysisContext, only_collect=false)
  if !b.is_subroutine_block enterscope!(ac) end
  for stmt in b.statements
    if !only_collect
      static_analysis(stmt, ac)
    elseif stmt isa VariableFactor
      if !hasvariable(ac, stmt.identifier, true)
        push!(b.subroutine.implicit_params, stmt.identifier)
      end
    elseif stmt isa Call
      push!(b.subroutine.calls_subroutines, stmt.identifier)
    end
  end
  if !b.is_subroutine_block popscope!(ac) end
  b.type = MPassed
end

function static_analysis(d::Declaration, ac::AnalysisContext)
  DEBUG && println("This is static analysis, analysing Declaration")
  static_analysis(d.var_type, ac)
  var_type = d.var_type.true_type
  for name in map(token->token.lexeme, d.names)
    if hasvariable(ac, name, true)
      entry = getvariable(ac, name)
      if !entry.is_implicit
        throw(StaticAnalysisException(
          "Cannot redeclare variable $(name) (line $(d.line))."
        ))
      end
    end
    id = create_var_id(ac, var_type, name)
    push!(d.unique_ids, id)
    pushvar!(ac, name, scalar_to_ref_type[var_type], id)
  end
  d.type = MPassed
end

function static_analysis(c::TypeOfVarOrValue, ac::AnalysisContext)
  typecheck(c.size, ac)
  if c.size.type != MInt
    throw(StaticAnalysisException(
      "Array size must have integer type (line $(c.line))."
    ))
  end
  c.type = MPassed
end

function static_analysis(a::Assignment, ac::AnalysisContext)
  DEBUG && println("This is static analysis, analysing Assignment")
  var_name = a.variable_name
  if !hasvariable(ac, var_name)
    throw(StaticAnalysisException(
      "Cannot assign to undeclared variable '$(var_name)' (line $(a.line))."
    ))
  end
  index_type = typecheck(a.array_index, ac)
  if index_type != MInt
    throw(StaticAnalysisException(
      "Array index must be an integer, got $(index_type) (line $(a.line))."
    ))
  end
  value_type = typecheck(a.value, ac)
  entry::SymTableEntry = getvariable(ac, var_name)
  var_type = entry.var_type
  a.var_type = var_type
  if var_type ∈ ref_types
    var_type = ref_type_to_scalar_type[var_type]
  elseif var_type ∈ scalar_types
    entry.val_identifier = create_var_id(ac, var_type, var_name)
  end
  if value_type != var_type
    if !(a.is_array_access && scalar_to_array_type[value_type] == var_type)
      throw(StaticAnalysisException(
        "Cannot assign a value of type $(value_type) to variable '$(var_name)' of type $(var_type) (line $(a.line))."
      ))
    end
  end
  a.var_unique_id = entry.val_identifier
  a.type = MPassed
end

function static_analysis(i::IfThenElse, ac::AnalysisContext)
  cond_type = typecheck(i.condition, ac)
  if cond_type != MBool
    throw(StaticAnalysisException(
    "If statement condition must be boolean, got type $cond_type (line $(i.line))."))
  end
  static_analysis(i.then_stmt, ac)
  static_analysis(i.else_stmt, ac)
  i.type = MPassed
end

function static_analysis(i::IfThen, ac::AnalysisContext)
  cond_type = typecheck(i.condition, ac)
  if cond_type != MBool
    throw(StaticAnalysisException(
    "If statement condition must be boolean, got type $cond_type (line $(i.line))."))
  end
  static_analysis(i.then_stmt, ac)
  i.type = MPassed
end

function static_analysis(w::While, ac::AnalysisContext)
  cond_type = typecheck(w.condition, ac)
  if cond_type != MBool
    throw(StaticAnalysisException(
    "While statement condition must be boolean, got type $cond_type (line $(w.line))."))
  end
  static_analysis(w.do_stmt, ac)
  w.type = MPassed
end

function static_analysis(r::Read, ac::AnalysisContext)
  for var::Variable in r.variables
    static_analysis(var, ac)
    if !hasvariable(ac, var.identifier)
      throw(StaticAnalysisException(
        "Trying to read to undeclared variable $(var.identifier) (line $(r.line))."
      ))
    end
    entry::SymTableEntry = getvariable(ac, var.identifier)
    push!(r.var_types, entry.var_type)
    if entry.var_type ∉ scalar_types
      throw(StaticAnalysisException(
        "Only scalar variable types can be read (line $(r.line))."
      ))
    end
  end
  r.type = MPassed
end

function static_analysis(p::Write, ac::AnalysisContext)
  DEBUG && println("This is static analysis, analysing Write")
  static_analysis(p.arguments, ac)
  p.type = MPassed
end

function static_analysis(a::Vector{Value}, ac::AnalysisContext)
  DEBUG && println("This is static analysis, analysing arguments")
  for arg in a
    typecheck(arg, ac)
  end
end

function typecheck(i::ImmediateInt, ac)
  i.type = MInt
end

function typecheck(i::ImmediateReal, ac)
  i.type = MReal
end

function typecheck(i::ImmediateString, ac)
  i.type = MString
end

function typecheck(i::ImmediateBool, ac)
  i.type = MBool
end

function typecheck(g::GetRef, ac)
  g.type = scalar_to_ref_type[typecheck(g.operand, ac)]
end

function typecheck(l::LiteralFactor, ac::AnalysisContext)
  DEBUG && println("This is static analysis, analysing LiteralFactor")
  class = l.token.class
  if class == int_literal
    l.type = MInt
  end
  if class == real_literal
    l.type = MReal
  end
  if class ∈ [kw_true, kw_false]
    l.type = MBool
  end
  if class == string_literal
    l.type = MString
  end
  return l.type
end

function typecheck(v::VariableFactor, ac::AnalysisContext)
  DEBUG && println("This is static analysis, analysing VariableFactor")
  DEBUG && println("in_call empty? $(isempty(ac.in_call))")
  var_name = v.identifier
  if !hasvariable(ac, var_name)
    throw(StaticAnalysisException(
      "Variable $(var_name) is not defined (line $(v.line))."
    ))
  end
  var_entry::SymTableEntry = getvariable(ac, var_name)
  v.variable_entry = var_entry
  var_type = var_entry.var_type
  if var_type ∈ ref_types
    var_type = ref_type_to_scalar_type[var_type]
  end
  v.type = var_type
end

function typecheck(a::ArrayAccessFactor, ac::AnalysisContext)
  DEBUG && println("This is static analysis, analysing ArrayAccessFactor")
  index_type = typecheck(a.index, ac)
  if index_type != MInt
    throw(StaticAnalysisException(
      "Array index should be an integer, got type $(index_type) (line $(a.line))."
    ))
  end
  var_name = a.identifier
  if !hasvariable(ac, var_name)
    throw(StaticAnalysisException(
      "Variable $(var_name) is not defined (line $(a.line))."
    ))
  end
  var_entry::SymTableEntry = getvariable(ac, var_name)
  a.variable_entry = var_entry
  var_type = var_entry.var_type
  if var_type ∉ array_types
    throw(StaticAnalysisException(
      "Cannot index into non-array type ($var_type) variable '$var_name' (line $(a.line))."
    ))
  end
  a.type = array_type_to_scalar_type[var_type]
end

function typecheck(c::CallFactor, ac::AnalysisContext)
  DEBUG && println("This is static analysis, analysing CallFactor")
  analyze_call(c, ac)
  subroutine_entry::SubroutineEntry = getsubroutine(ac, c.identifier)
  ret_type = subroutine_entry.return_type
  if ret_type == MNothing
    name = subroutine_entry.name
    throw(StaticAnalysisException(
      "A call to subroutine '$name' cannot be used as a value as it has no return type (line $(c.line))."
    ))
  end
  c.type = ret_type
end

function typecheck(p::ParenFactor, ac::AnalysisContext)
  p.type = typecheck(p.expression, ac)
end

function typecheck(n::NotFactor, ac::AnalysisContext)
  arg_type = typecheck(n.argument, ac)
  if arg_type != MBool
    throw(StaticAnalysisException(
      "The \"not\" operation requires a boolean argument, got type $(arg_type) on line $(n.line)."
    ))
  end
  n.type = MBool
end

# function typecheck(t::Term, ac::AnalysisContext)
#   factor_types = [typecheck(tpl[1], ac) for tpl in t.factors]
#   first_type = factor_types[1]
#   for tpl in t.factors
#     if tpl[1].type != first_type
#       throw(StaticAnalysisException(
#         "Type mismatch in operation $(tpl[2].class) on line $(tpl[1].line)."
#       ))
#     end
#   end
#   t.type = first_type
#   return t.type
# end

function typecheck(s::SizeFactor, ac::AnalysisContext)
  type = typecheck(s.array, ac)
  if type ∉ [MIntArray, MRealArray, MBoolArray, MStringArray]
    throw(StaticAnalysisException(
      "Array size is not defined for values of type $(type) (on line $(s.line))."))
  end
  s.type = MInt
end

function get_llvm_function_name(subroutine_name::String, call_id::Int)
  return "@$(subroutine_name).$(call_id)"
end

max_stars_repo_name: juhamyllari/minipascal-julia

max_issues_repo_name: juhamyllari/minipascal-julia

max_forks_repo_name: juhamyllari/minipascal-julia


text:
# https://fenicsproject.org/docs/dolfin/latest/python/demos/auto-adaptive-poisson/demo_auto-adaptive-poisson.py.html

module auto_adaptive_poisson

using FenicsPy

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, "Lagrange", 1)

# Define boundary condition
u0 = FeFunction(V)
bc = DirichletBC(V, u0, "x[0] < $DOLFIN_EPS || x[0] > 1.0 - $DOLFIN_EPS")

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",
	                      degree=1)
g = Expression("sin(5*x[0])", degree=1)
a = inner(grad(u), grad(v))*dx()
L = f*v*dx() + g*v*ds()

# Define function for the solution
u = FeFunction(V)

# Define goal functional (quantity of interest)
M = u*dx()

# Define error tolerance
tol = 1.e-5

# Solve equation a = L with respect to u and the given boundary
# conditions, such that the estimated error (measured in M) is less
# than tol
problem = LinearVariationalProblem(a, L, u, bc)
solver = AdaptiveLinearVariationalSolver(problem, M)
solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "cg"
solver.parameters["error_control"]["dual_variational_solver"]["symmetric"] = true
solver.solve(tol)

solver.summary()

# Plot solution(s)
#plot(u.root_node(), title="Solution on initial mesh")

#plot(u.leaf_node(), title="Solution on final mesh")

end # module auto_adaptive_poisson


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Severo_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Severo")
JLLWrappers.@generate_main_file("Severo", UUID("0d93ea7f-3968-595b-b8f9-5c7f63ac8648"))
end  # module Severo_jll


text:
using .CUDA

export CuBlobDescriptor, CuPODBlobDescriptor, CuTensorBlobDescriptor, CuFilterBlobDescriptor
export CuTensorBlob

immutable CuTensorBlob{T<:AbstractFloat,N} <: Blob{T,N}
  ptr   :: CuPtr
  shape :: NTuple{N, Int}
  len   :: Int
end
function CuTensorBlob{T<:AbstractFloat, N}(dtype::Type{T}, dims::NTuple{N,Int})
  len = prod(dims)
  ptr = CUDA.cualloc(dtype, len)
  return CuTensorBlob{T,N}(ptr, dims, len)
end

length(b::CuTensorBlob) = b.len
size(b::CuTensorBlob) = b.shape

function copy!{T}(dst :: CuTensorBlob{T}, src :: Array{T})
  @assert length(dst) == length(src)
  CuBLAS.set_vector(src, dst.ptr)
end
function copy!{T}(dst :: Array{T}, src :: CuTensorBlob{T})
  @assert length(dst) == length(src)
  CuBLAS.get_vector(src.ptr, dst)
end
function copy!{T}(dst :: CuTensorBlob{T}, src :: CuTensorBlob{T})
  @assert length(dst) == length(src)
  @CUDA.cucall(:cuMemcpy, (Ptr{Void}, Ptr{Void}, Cint), dst.ptr.p, src.ptr.p, length(dst)*sizeof(T))
end
function fill!{T}(dst :: CuTensorBlob{T}, val)
  val_vec = Array{T}(length(dst))
  fill!(val_vec, val)
  copy!(dst, val_vec)
end
function erase!{T}(dst :: CuTensorBlob{T})
  @CUDA.cucall(:cuMemsetD8_v2, (Ptr{Void}, Cuchar, Csize_t), dst.ptr.p, 0, length(dst)*sizeof(T))
end

function make_blob{N}(backend::GPUBackend, data_type::Type, dims::NTuple{N,Int})
  return CuTensorBlob(data_type, dims)
end
function reshape_blob{T,N}(backend::GPUBackend, blob::CuTensorBlob{T}, dims::NTuple{N,Int})
  @assert prod(dims) == length(blob)
  return CuTensorBlob{T,N}(blob.ptr, dims, length(blob))
end
function destroy(blob :: CuTensorBlob)
  if blob.ptr.p != 0
    CUDA.free(blob.ptr)
    blob.ptr.p = 0
  end
end


text:
using .Gadfly
using Compose: compose, context

"""
`plot_pressure(well::Wellbore, pressures, ctitle = nothing)`

Plot pressure profile for a given wellbore using the pressure outputs from one of the pressure traverse functions.

See `traverse_topdown` and `pressure_and_temp`.
"""
function plot_pressure(well::Wellbore, pressures, ctitle = nothing)

        plot(x = pressures, y = well.md, Geom.path, Theme(default_color = "deepskyblue"),
                Scale.x_continuous(format = :plain),
                Guide.xlabel("Pressure (psia)"),
                Scale.y_continuous(format = :plain),
                Guide.ylabel("Measured Depth (ft)"),
                Guide.title(ctitle),
                Coord.cartesian(yflip = true))
end

"""
`plot_pressure(m::WellModel, pressures, ctitle = nothing)`

Plot pressure profile for a given wellbore using the pressure outputs from one of the pressure traverse functions.

The `wellbore` field must be defined in the passed WellModel.

See `traverse_topdown` and `pressure_and_temp`.
"""
function plot_pressure(m::WellModel, pressures, ctitle = nothing)

        plot_pressure(m.wellbore, pressures, ctitle)
end

"""
`function plot_pressures(well::Wellbore, tubing_pressures, casing_pressures, ctitle = nothing, valvedepths = [])`

Plot relevant gas lift pressures for a given wellbore and set of calculated pressures.

See `traverse_topdown`, `casing_traverse_topdown`, and `pressure_and_temp`.
"""
function plot_pressures(well::Wellbore, tubing_pressures, casing_pressures, ctitle = nothing, valvedepths = [])

        plot(layer(x = tubing_pressures, y = well.md, Geom.path, Theme(default_color = "deepskyblue")),
             layer(x = casing_pressures, y = well.md, Geom.path, Theme(default_color = "springgreen")),
             layer(yintercept = valvedepths, Geom.hline(color = "black", style = :dash)),
                Scale.x_continuous(format = :plain),
                Guide.xlabel("Pressure (psia)"),
                Scale.y_continuous(format = :plain),
                Guide.ylabel("Measured Depth (ft)"),
                Guide.title(ctitle),
                Coord.cartesian(yflip = true))
end

"""
`plot_pressures(m::WellModel, tubing_pressures, casing_pressures, ctitle = nothing)`

Plot relevant gas lift pressures for a given wellbore and set of calculated pressures.

The `wellbore` field must be defined in the passed WellModel, with the `valves` field optional.

See `traverse_topdown`, `casing_traverse_topdown`, and `pressure_and_temp`.
"""
function plot_pressures(m::WellModel, tubing_pressures, casing_pressures, ctitle = nothing)

        valvedepths = m.valves === missing ? [] : m.valves.md

        plot_pressures(m.wellbore, tubing_pressures, casing_pressures, ctitle, valvedepths)
end

"""
`plot_temperature(well::Wellbore, temps, ctitle = nothing)`

Plot temperature profile for a given wellbore using the pressure outputs from one of the pressure traverse functions.

See `linear_wellboretemp` and `Shiu_wellboretemp`.
"""
function plot_temperature(well::Wellbore, temps, ctitle = nothing)

        plot(x = temps, y = well.md, Geom.path, Theme(default_color = "red"),
                Scale.x_continuous(format = :plain),
                Guide.xlabel("Temperature (°F)"),
                Scale.y_continuous(format = :plain),
                Guide.ylabel("Measured Depth (ft)"),
                Guide.title(ctitle),
                Coord.cartesian(yflip = true))
end

"""
`plot_pressureandtemp(well::Wellbore, tubing_pressures, casing_pressures, temps, ctitle = nothing, valvedepths = [])`

Plot pressure & temperature profiles for a given wellbore using the pressure & temperature outputs from the pressure traverse & temperature functions.

See `traverse_topdown`,`pressure_and_temp`, `linear_wellboretemp`, `Shiu_wellboretemp`.
"""
function plot_pressureandtemp(well::Wellbore, tubing_pressures, casing_pressures, temps, ctitle = nothing, valvedepths = [])

        pressure = plot(layer(x = tubing_pressures, y = well.md, Geom.path, Theme(default_color = "deepskyblue")),
                        layer(x = casing_pressures, y = well.md, Geom.path, Theme(default_color = "mediumspringgreen")),
                        layer(yintercept = valvedepths, Geom.hline(color = "black", style = :dash)),
                Scale.x_continuous(format = :plain),
                Guide.xlabel("psia"),
                Scale.y_continuous(format = :plain),
                Guide.ylabel("Measured Depth (ft)"),
                Guide.title(ctitle),
                Coord.cartesian(yflip = true),
                Theme(plot_padding=[5mm, 0mm, 5mm, 5mm]))

        placeholdertitle = ctitle === nothing ? nothing : " "
        temp = plot(x = temps, y = well.md, Geom.path, Theme(default_color = "red"),
                Scale.x_continuous(format = :plain),
                Guide.xlabel("°F"),
                Scale.y_continuous(labels = nothing),
                Guide.yticks(label = false),
                Guide.ylabel(nothing),
                Guide.title(placeholdertitle),
                Coord.cartesian(yflip = true),
                Theme(default_color = "red", plot_padding=[5mm, 5mm, 5mm, 5mm]))

        hstack(compose(context(0, 0, 0.75, 1), render(pressure)),
                compose(context(0.75, 1, 0.5, 1), render(temp)))
end

"""
`plot_pressureandtemp(m::WellModel, tubing_pressures, casing_pressures, ctitle = nothing)`

Plot pressure & temperature profiles for a given wellbore using the pressure & temperature outputs from the pressure traverse & temperature functions.

The `wellbore` and `temperatureprofile` fields must be defined in the passed WellModel, with the `valves` field optional.

See `traverse_topdown`,`pressure_and_temp`, `linear_wellboretemp`, `Shiu_wellboretemp`.
"""
function plot_pressureandtemp(m::WellModel, tubing_pressures, casing_pressures, ctitle = nothing)

        valvedepths = m.valves === missing ? [] : m.valves.md

        plot_pressureandtemp(m.wellbore, tubing_pressures, casing_pressures, m.temperatureprofile, ctitle, valvedepths)
end

"""
`plot_gaslift(well::Wellbore, tubing_pressures, casing_pressures, temps, valvedata, ctitle = nothing)`

Plot pressure & temperature profiles along with valve depths and opening/closing pressures for a gas lift well.

Requires a valve table in the same format as returned by the `valve_calcs` function.

See `traverse_topdown`,`pressure_and_temp`, `linear_wellboretemp`, `Shiu_wellboretemp`, `valve_calcs`.
"""
function plot_gaslift(well::Wellbore, tubing_pressures, casing_pressures, temps, valvedata, ctitle = nothing)

        valvedepths = valvedata[:,2]

        pressure = plot(layer(x = [valvedata[:,12];valvedata[:,13]], y = [valvedepths;valvedepths], Geom.point, Theme(default_color = "mediumpurple3")), #PVC and PVO
                        layer(x = tubing_pressures, y = well.md, Geom.path, Theme(default_color = "deepskyblue")),
                        layer(x = casing_pressures, y = well.md, Geom.path, Theme(default_color = "mediumspringgreen")),
                        layer(yintercept = valvedepths, Geom.hline(color = "black", style = :dash)),

                Scale.x_continuous(format = :plain),
                Guide.xlabel("psia"),
                Scale.y_continuous(format = :plain),
                Guide.ylabel("Measured Depth (ft)"),
                Guide.title(ctitle),
                Coord.cartesian(yflip = true),
                Theme(plot_padding=[5mm, 0mm, 5mm, 5mm]))

        placeholdertitle = ctitle === nothing ? nothing : " " #generate a blank title to align the top of the plots if needed
        temp = plot(x = temps, y = well.md, Geom.path, Theme(default_color = "red"),
                Scale.x_continuous(format = :plain),
                Guide.xlabel("°F"),
                Scale.y_continuous(labels = nothing),
                Guide.yticks(label = false),
                Guide.ylabel(nothing),
                Guide.title(placeholdertitle),
                Coord.cartesian(yflip = true),
                Guide.manual_color_key("",
                                       ["TP", "CP", "Valves", "PVO/PVC"],
                                       ["deepskyblue", "mediumspringgreen", "black", "mediumpurple3"]),
                Theme(default_color = "red", plot_padding=[5mm, 5mm, 5mm, 5mm]))

        hstack(compose(context(0, 0, 0.65, 1), render(pressure)),
                compose(context(0.65, 0, 0.35, 1), render(temp)))
end

"""
`plot_gaslift(m::WellModel, tubing_pressures, casing_pressures, valvedata, ctitle = nothing)`

Plot pressure & temperature profiles along with valve depths and opening/closing pressures for a gas lift well.

Requires a valve table in the same format as returned by the `valve_calcs` function. The passed WellModel must also have the `wellbore`, `temperatureprofile`, and `valves` fields defined.

See `traverse_topdown`,`pressure_and_temp`, `linear_wellboretemp`, `Shiu_wellboretemp`, `valve_calcs`.
"""
function plot_gaslift(m::WellModel, tubing_pressures, casing_pressures, valvedata, ctitle = nothing)

        plot_gaslift(m.wellbore, tubing_pressures, casing_pressures, m.temperatureprofile, valvedata, ctitle)
end