text:
# Update_weights.jl

function update_lλ!(state::State_GPMTD, prior_λ::Vector{Float64})

    L = length(state.lλ) - 1

    Nζ = StatsBase.counts(state.ζ, 0:L)
    α1_λ = prior_λ .+ Nζ
    state.lλ = SparseProbVec.rDirichlet(α1_λ, logout=true)

    return nothing
end
function update_lλ!(state::State_GPMTD, prior_λ::SparseDirMix)

    L = length(state.lλ) - 1

    Nζ = StatsBase.counts(state.ζ, 0:L)
    α1_λ = prior_λ.α .+ Nζ
    d = SparseProbVec.SparseDirMix(α1_λ, prior_λ.β)
    state.lλ = SparseProbVec.rand(d, logout=true)

    return nothing
end
function update_lλ!(state::State_GPMTD, prior_λ::SBMprior)

    L = length(state.lλ) - 1

    Nζ = StatsBase.counts(state.ζ, 0:L)
    post_lλ = SparseProbVec.SBM_multinom_post(prior_λ, Nζ)
    state.lλ = SparseProbVec.rand(post_lλ, logout=true)

    return nothing
end


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

"""
    module ModiaMath.SimulationEngine

Simulation engine for implicit index 1 DAE models with events.

# Main developer

[Martin Otter](https://rmc.dlr.de/sr/de/staff/martin.otter/),
[DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module SimulationEngine

export simulate!

import Sundials
import ModiaMath
import ModiaMath.DAE
using  StaticArrays
@eval using Printf
using LinearAlgebra

# include code
# include("IDA_UserData.jl")  # is included within simulate.jl
include("simulate.jl")

end

ext: jl


text:
using Kinetic, Plots, LinearAlgebra
using KitBase.JLD2
using Flux: onecold
cd(@__DIR__)

begin
    set = Setup(case = "sod", space = "1d2f1v", maxTime = 0.15)
    ps = PSpace1D(0.0, 1.0, 200, 1)
    vs = VSpace1D(-5.0, 5.0, 100)
    #gas = Gas(Kn = 1e-4, K = 2, γ = 5/3)
    gas = Gas(Kn = 1e-3, K = 2, γ = 5/3)
    #gas = Gas(Kn = 1e-2, K = 2, γ = 5/3)
    ib = IB2F(ib_sod(set, ps, vs, gas)...)
    ks = SolverSet(set, ps, vs, gas, ib)
end

# Kn = 1e-4
ctrs = []
begin
    @load "kn=1e-4/nn.jld2" ctr
    push!(ctrs, ctr)
    @load "kn=1e-4/kngll.jld2" ctr
    push!(ctrs, ctr)
    @load "kn=1e-4/pure_kinetic.jld2" ctr
    push!(ctrs, ctr)
    @load "kn=1e-4/pure_ns.jld2" ctr
    push!(ctrs, ctr)
end

# Kn = 1e-3
ctrs = []
begin
    @load "kn=1e-3/nn.jld2" ctr
    push!(ctrs, ctr)
    @load "kn=1e-3/kngll.jld2" ctr
    push!(ctrs, ctr)
    @load "kn=1e-3/pure_kinetic.jld2" ctr
    push!(ctrs, ctr)
    @load "kn=1e-3/pure_ns.jld2" ctr
    push!(ctrs, ctr)
end

# Kn = 1e-2
ctrs = []
begin
    @load "kn=1e-2/nn.jld2" ctr
    push!(ctrs, ctr)
    @load "kn=1e-2/kngll.jld2" ctr
    push!(ctrs, ctr)
    @load "kn=1e-2/pure_kinetic.jld2" ctr
    push!(ctrs, ctr)
    @load "kn=1e-3/pure_ns.jld2" ctr
    push!(ctrs, ctr)
    @load "kn=1e-2/pure_kinetic.jld2" ctr
end

sols = []
begin
    for idx in eachindex(ctrs)
        sol = zeros(ks.ps.nx, 3)
        for i in axes(sol, 1)
            sol[i, :] .= ctrs[idx][i].prim
            sol[i, end] = 1 / sol[i, end]
        end
        push!(sols, sol)
    end
end

begin
    ctr = ctrs[3]

    rg_ce = zeros(Int, ks.ps.nx)
    for i = 1:ks.ps.nx
        sw = (ctr[i+1].w .- ctr[i-1].w) / ks.ps.dx[i] / 2.0
        x, y = regime_data(ks, ctr[i].w, sw, ctr[i].h)
        rg_ce[i] = floor(y) |> Int
    end

    rg_kn = zeros(Int, ks.ps.nx)
    KnGLL = zeros(ks.ps.nx)
    for i = 1:ks.ps.nx
        sw = (ctr[i+1].w .- ctr[i-1].w) / ks.ps.dx[i] / 2.0
        L = abs(ctr[i].w[1] / sw[1])
        ℓ = (1/ctr[i].prim[end])^ks.gas.ω / ctr[i].prim[1] * sqrt(ctr[i].prim[end]) * ks.gas.Kn

        KnGLL[i] = ℓ / L
        rg_kn[i] = ifelse(KnGLL[i] >= 0.05, 1, 0)
    end
end

scatter(ks.ps.x[1:ks.ps.nx], rg_ce, alpha=0.7, label="NN", xlabel="x", ylabel="regime")
scatter!(ks.ps.x[1:ks.ps.nx], rg_kn, alpha=0.7, label="KnGLL")
plot!(ks.ps.x[1:ks.ps.nx], rg_ce, lw=1.5, line=:dot, color=:gray27, label="True", xlabel="x", ylabel="regime")

#savefig("figs/kn4_regime.pdf")
#savefig("figs/kn3_regime.pdf")
#savefig("figs/kn2_regime.pdf")

function curve(idx, s, lgd=:topright; ns=true)
    fig = plot(ks.ps.x[idx], sols[1][idx, s], lw=1.5, label="NN", legend=lgd)
    plot!(fig, ks.ps.x[idx], sols[2][idx, s], lw=1.5, line=:dash, label="KnGLL")
    plot!(fig, ks.ps.x[idx], sols[3][idx, s], lw=1.5, line=:dashdot, label="Kinetic")
    ns && plot!(fig, ks.ps.x[idx], sols[4][idx, s], lw=2, line=:dot, color=:gray27, label="NS")

    xlabel!(fig, "x")
    if s == 1
        ylabel!(fig, "ρ")
    elseif s == 2
        ylabel!(fig, "U")
    elseif s == 3
        ylabel!(fig, "T")
    end

    return fig
end

idx = 1:ks.ps.nx
curve(idx, 1)
#savefig("kn4_n.pdf")
#savefig("kn3_n.pdf")
#curve(idx, 1; ns=false)
#savefig("kn2_n.pdf")

curve(idx, 3)
#savefig("kn4_t.pdf")
#savefig("kn3_t.pdf")
#curve(idx, 3; ns=false)
#savefig("kn2_t.pdf")

idx = ks.ps.nx÷5*3+1:ks.ps.nx÷5*4
curve(idx, 1)
#savefig("kn4_n_zoom.pdf")
#savefig("kn3_n_zoom.pdf")
#curve(idx, 1; ns=false)
#savefig("kn2_n_zoom.pdf")

curve(idx, 3, :topleft)
#savefig("kn4_t_zoom.pdf")
#savefig("kn3_t_zoom.pdf")
#curve(idx, 3, :topleft; ns=false)
#savefig("kn2_t_zoom.pdf")

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
using GLVisualize, GeometryTypes, Colors, GLAbstraction, Reactive, Images, ModernGL
import GLVisualize: mm, annotated_text

if !isdefined(:runtests)
    window = glscreen()
    timesignal = loop(linspace(0f0, 1f0, 360))
end

description = """
Example of showing animated annotations over an image stream
"""

N = 20
positions = foldp(rand(Point2f0, N) .* 500f0, timesignal) do v0, t
    map!(v0) do p
        map(p .+ (rand(Point2f0) .- 0.5f0) .* 10f0) do pi
            clamp(pi, 0f0, 500f0) # restrict to image
        end
    end
end
base = [string(i) for i=1:N]
label_str = foldp(copy(base), timesignal) do v0, t
    map!(v0, base) do baselabel
        string(baselabel, ' ', round(t, 3))
    end
end

labels, widths = annotated_text(positions, label_str)
fbuffer = foldp(rand(RGB{N0f8}, 500, 500), timesignal) do v0, t
    map!(v0) do v
        RGB(v.r > 0.5 ? sin(v.r) : cos(v.r), (sin(t) + 1) / 2, 0.7) # yeah, real creative..
    end
end
frame_viz = visualize(fbuffer)

# instead of circle you can also use unicode charactes (e.g. '+')
position_viz = visualize(
    (Circle{Float32}(Point2f0(0), 1.5mm), positions),
    color = RGBA(1f0, 0f0, 0f0, 0.6f0)
)
gpu_position = position_viz.children[][:position]
bg_viz = visualize(
    (ROUNDED_RECTANGLE, gpu_position),
    scale = widths,
    # you can give each label a color (random looks terrible, which is why whe use just white instead)
    # color = RGBA{Float32}[rand(RGB) for i = 1:N],
    color = RGBA(1f0, 1f0, 1f0, 0.6f0),
    offset = Vec2f0(1.0mm)
)
_view(frame_viz, window)
_view(bg_viz, window)
_view(labels, window)
_view(position_viz, window)

center!(window, :orthographic_pixel) # center orthographic_pixel camera to what we just visualized

if !isdefined(:runtests)
    renderloop(window)
end


text:
using LiterateOrg
if VERSION.major > 0 || VERSION.minor > 6
    using Test
else
    using Base.Test
end

const testfile = joinpath(dirname(@__FILE__), "literate_org_tangled_tests.jl")
if isfile(testfile)
    include(testfile)
else
    error("LiterateOrg not properly installed. Please run Pkg.build(\"LiterateOrg\") then restart Julia.")
end


text:
using Printf

function balancedbrackets(str::AbstractString)
    i = 0
    for c in str
        if c == '[' i += 1 elseif c == ']' i -=1 end
        if i < 0 return false end
    end
    return i == 0
end

brackets(n::Integer) = join(shuffle(collect(Char, "[]" ^ n)))

for (test, pass) in map(x -> (x, balancedbrackets(x)), collect(brackets(i) for i = 0:8))
    @printf("%22s%10s\n", test, pass ? "pass" : "fail")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Mocha, Boltzmann

type HDF5DataSource
  sys::System
  data_layer::HDF5DataLayer
  state::Mocha.HDF5DataLayerState

  HDF5DataSource(layer::HDF5DataLayer) = begin
    sys = System(CPUBackend())
    inputs = Array(Blob, 0)
    diffs = Array(Blob, 0)
    state = setup(sys, layer, inputs, diffs)
    new(sys, layer, state)
  end
end

function get_data_batch(source::HDF5DataSource)
  forward(source.sys, source.state, Array(Blob, 0))
  X = source.state.blobs[1].data
  reshape(X, (size(X,1)*size(X,2), size(X,4)))
end

# TODO: data source has to always output data between 0 and 1! seems to be the case for our data.
function fit(rbm::Boltzmann.RBM, source::HDF5DataSource, lr=0.1, n_iter=10, n_gibbs=1)
  batch_size = source.data_layer.batch_size
  n_samples = size(source.state.curr_hdf5_file["data"],4)
  n_batches = int(ceil(n_samples / batch_size))
  w_buf = zeros(size(rbm.W))
  for itr=1:n_iter
    for i=1:n_batches
      info("Iteration $(itr): fitting batch $(i)/$(n_batches)")
      batch = convert(Array{Float64,2}, get_data_batch(source))
      Boltzmann.fit_batch!(rbm, batch, buf=w_buf, n_gibbs=n_gibbs)
    end
    #@printf("Iteration #%s, pseudo-likelihood = %s\n",
    #        itr, mean(Boltzmann.score_samples(rbm, X)))
   end
end

type DBN
  rbms::Vector{GRBM}

  DBN(dims::Vector{Int}) = begin
    rbms = Array(GRBM, 0)
    for k=1:length(dims)-1
      push!(rbms, GRBM(dims[k],dims[k+1]))
    end
    new(rbms)
  end
end

mh(rbm::GRBM, vis::Array{Float64, 2}) = Boltzmann.mean_hiddens(rbm, vis)

function mh_at_layer(dbn::DBN, batch::Array{Float64, 2}, layer::Int)
  hiddens = Array(Array{Float64, 2}, layer)
  hiddens[1] = mh(dbn.rbms[1], batch)
  for k=2:layer
    hiddens[k] = mh(dbn.rbms[k], hiddens[k-1])
  end
  hiddens[end]
end

function fit(dbn::DBN, source::HDF5DataSource, lr=0.1, n_iter=10, n_gibbs=1)
  batch_size = source.data_layer.batch_size
  n_samples = size(source.state.curr_hdf5_file["data"],4)
  n_batches = int(ceil(n_samples / batch_size))
  for k = 1:length(dbn.rbms)
    println("Training layer $k/$(length(dbn.rbms))")
    w_buf = zeros(size(dbn.rbms[k].W))
    for itr=1:n_iter
      for i=1:n_batches
        info("Iteration $(itr): fitting batch $(i)/$(n_batches)")
        batch = convert(Array{Float64,2}, get_data_batch(source))
        if k == 1
          input = batch
        else
          input = mh_at_layer(dbn, batch, k)
        end
        Boltzmann.fit_batch!(dbn.rbms[k], input, buf=w_buf, n_gibbs=n_gibbs)
      end
    end
  end
end


text:
struct FRAPMesh <: Mesh
    Nodes::Array{Float64,1}
    NBases::Int
    Fil::IndexDict
end
# Convenience constructors
function FRAPMesh(
    model::SFFM.Model,
    Nodes::Array{<:Real,1},
    NBases::Int;
    Fil::IndexDict=IndexDict(),
    v::Bool = false,
)

    ## Construct the sets Fᵐ = ⋃ᵢ Fᵢᵐ, global index for sets of type m
    if isempty(Fil)
        Fil = MakeFil(model, Nodes)
    end

    mesh = FRAPMesh(
        Nodes,
        NBases,
        Fil,
    )
    v && println("UPDATE: DGMesh object created with fields ", fieldnames(SFFM.DGMesh))
    return mesh
end
function FRAPMesh()
    FRAPMesh(
        Array{Float64,1}(undef,0),
        0,
        Dict{String,BitArray{1}}(),
    )
end

"""

    NBases(mesh::FRAPMesh)

Number of bases in a cell
"""
NBases(mesh::FRAPMesh) = mesh.NBases

"""

    CellNodes(mesh::FRAPMesh)

The cell centre
"""
CellNodes(mesh::FRAPMesh) = Array(((mesh.Nodes[1:end-1] + mesh.Nodes[2:end]) / 2 )')

"""

    Basis(mesh::FRAPMesh)

Constant ""
"""
Basis(mesh::FRAPMesh) = ""

function MakeLazyGenerator(
    model::SFFM.Model,
    mesh::FRAPMesh,
    me::ME;
    v::Bool=false,
)
    blocks = (me.s*me.a, me.S, me.s*me.a)

    boundary_flux = (in = me.s[:], out = me.a[:])

    T = model.T
    C = model.C
    delta = Δ(mesh)
    D = me.D

    signChangeIndex = zeros(Bool,NPhases(model),NPhases(model))
    for i in 1:NPhases(model), j in 1:NPhases(model)
        if ((sign(model.C[i])!=0) && (sign(model.C[j])!=0))
            signChangeIndex[i,j] = (sign(model.C[i])!=sign(model.C[j]))
        elseif (sign(model.C[i])==0)
            signChangeIndex[i,j] = sign(model.C[j])>0
        elseif (sign(model.C[j])==0)
            signChangeIndex[i,j] = sign(model.C[i])>0
        end
    end

    out = SFFM.LazyGenerator(blocks,boundary_flux,T,C,delta,D,signChangeIndex,mesh.Fil)
    v && println("UPDATE: LazyGenerator object created with keys ", keys(out))
    return out
end
function MakeLazyGenerator(model::Model, mesh::FRAPMesh; v::Bool=false)
    me = SFFM.MakeME(SFFM.CMEParams[NBases(mesh)])
    return MakeLazyGenerator(model, mesh, me; v=v)
end
function MakeFullGenerator(model::Model, mesh::Mesh, me::ME; v::Bool=false)
    lazy = MakeLazyGenerator(model,mesh,me; v=v)
    return materialise(lazy)
end

# function MakeB(model::Model, mesh::FRAPMesh, me::ME)
#     N₊ = sum(model.C .>= 0)
#     N₋ = sum(model.C .<= 0)

#     F = Dict{String,SparseArrays.SparseMatrixCSC{Float64,Int64}}()
#     UpDiagBlock = me.s*me.a
#     LowDiagBlock = me.s*me.a
#     for i = ["+","-"]
#         F[i] = SparseArrays.spzeros(Float64, TotalNBases(mesh), TotalNBases(mesh))
#         for k = 1:NIntervals(mesh)
#             idx = (1:NBases(mesh)) .+ (k - 1) * NBases(mesh)
#             if k > 1
#                 idxup = (1:NBases(mesh)) .+ (k - 2) * NBases(mesh)
#                 if i=="+"
#                     F[i][idxup, idx] = UpDiagBlock
#                 elseif i=="-"
#                     F[i][idx, idxup] = LowDiagBlock
#                 end # end if C[i]
#             end # end if k>1
#         end # for k in ...
#     end # for i in ...

#     signChangeIndex = zeros(Bool,NPhases(model),NPhases(model))
#     for i in 1:NPhases(model), j in 1:NPhases(model)
#         if ((sign(model.C[i])!=0) && (sign(model.C[j])!=0))
#             signChangeIndex[i,j] = (sign(model.C[i])!=sign(model.C[j]))
#         elseif (sign(model.C[i])==0)
#             signChangeIndex[i,j] = sign(model.C[j])>0
#         elseif (sign(model.C[j])==0)
#             signChangeIndex[i,j] = sign(model.C[i])>0
#         end
#     end
#     B = SparseArrays.spzeros(
#         Float64,
#         NPhases(model) * TotalNBases(mesh) + N₋ + N₊,
#         NPhases(model) * TotalNBases(mesh) + N₋ + N₊,
#     )
#     B[N₋+1:end-N₊,N₋+1:end-N₊] = SparseArrays.kron(
#             model.T.*signChangeIndex,
#             SparseArrays.kron(SparseArrays.I(NIntervals(mesh)),me.D)
#         ) + SparseArrays.kron(
#             model.T.*(1 .- signChangeIndex),
#             SparseArrays.I(TotalNBases(mesh))
#         )

#     # Boundary conditions
#     T₋₋ = model.T[model.C.<=0,model.C.<=0]
#     T₊₋ = model.T[model.C.>=0,:].*((model.C.<0)')
#     T₋₊ = model.T[model.C.<=0,:].*((model.C.>0)')
#     T₊₊ = model.T[model.C.>=0,model.C.>=0]
#     # yuck
#     inLower = [
#         LinearAlgebra.kron(LinearAlgebra.diagm(abs.(model.C).*(model.C.<=0)),me.s)[:,model.C.<=0];
#         LinearAlgebra.zeros((NIntervals(mesh)-1)*NPhases(model)*NBases(mesh),N₋)
#     ]
#     outLower = [
#         LinearAlgebra.kron(T₋₊,me.a) LinearAlgebra.zeros(N₋,N₊+(NIntervals(mesh)-1)*NPhases(model)*NBases(mesh))
#     ]
#     inUpper = [
#         LinearAlgebra.zeros((NIntervals(mesh)-1)*NPhases(model)*NBases(mesh),N₊);
#         LinearAlgebra.kron(LinearAlgebra.diagm(abs.(model.C).*(model.C.>=0)),me.s)[:,model.C.>=0]
#     ]
#     outUpper = [
#         LinearAlgebra.zeros(N₊,N₋+(NIntervals(mesh)-1)*NPhases(model)*NBases(mesh)) LinearAlgebra.kron(T₊₋,me.a)
#     ]

#     QBDidx = MakeQBDidx(model,mesh)
#     B[1:N₋,QBDidx] = [T₋₋ outLower]
#     B[end-N₊+1:end,QBDidx] = [outUpper T₊₊]
#     B[QBDidx[N₋+1:end-N₊],1:N₋] = inLower
#     B[QBDidx[N₋+1:end-N₊],(end-N₊+1):end] = inUpper
#     for i = 1:NPhases(model)
#         idx = ((i-1)*TotalNBases(mesh)+1:i*TotalNBases(mesh)) .+ N₋
#         if model.C[i] > 0
#             B[idx, idx] += model.C[i] * (SparseArrays.kron(
#                     SparseArrays.I(NIntervals(mesh)), me.S
#                     ) + F["+"])
#         elseif model.C[i] < 0
#             B[idx, idx] += abs(model.C[i]) * (SparseArrays.kron(
#                 SparseArrays.I(NIntervals(mesh)), me.S
#                 ) + F["-"])
#         end
#     end

#     BDict = MakeDict(B, model, mesh)

#     return FullGenerator(BDict, B, mesh.Fil)
# end

# function MakeB(model::Model, mesh::FRAPMesh, order::Int)
#     me = SFFM.MakeME(SFFM.CMEParams[order], mean = 1)
#     return MakeB(model, mesh, me)
# end

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
import Base: show

# TODO: real expression printing
Base.show(io::IO, x::Variable) = print(io, "%", x.id)

print_stmt(io::IO, ex) = print(io, ex)

function show(io::IO, b::Branch)
  if b == unreachable
    print(io, "unreachable")
  elseif isreturn(b)
    print(io, "return $(repr(b.args[1]))")
  else
    print(io, "br $(b.block)")
    if !isempty(b.args)
      print(io, " (")
      join(io, repr.(b.args), ", ")
      print(io, ")")
    end
    b.condition != nothing && print(io, " unless $(b.condition)")
  end
end

const tab = "  "

function printargs(io::IO, args, types = [Any for arg in args])
  print(io, "(")
  for i = 1:length(args)
    print(io, args[i])
    types[i] != Any && print(io, " :: ", types[i])
    i != length(args) && print(io, ", ")
  end
  print(io, ")")
end

function show(io::IO, b::Block)
  indent = get(io, :indent, 0)
  bb = BasicBlock(b)
  print(io, tab^indent)
  print(io, b.id, ":")
  if !isempty(bb.args)
    print(io, " ")
    printargs(io, bb.args, bb.argtypes)
  end
  for (x, st) in b
    println(io)
    print(io, tab^indent, "  ")
    x == nothing || print(io, string("%", x.id), " = ")
    st.expr == nothing ? print(io, "nothing") :
      print_stmt(io, st.expr)
    st.type == Any || print(io, " :: ", st.type)
  end
  for br in bb.branches
    println(io)
    print(io, tab^indent, "  ", br)
  end
end

function show(io::IO, ir::IR)
  show(io, block(ir, 1))
  for b in blocks(ir)[2:end]
    println(io)
    show(io, b)
  end
end

function print_stmt(io::IO, ex::IR)
  io = IOContext(io, :indent=>get(io, :indent, 0)+2)
  println(io)
  show(io, ex)
end


text:
using JuliaVariables
using Test
using NameResolution
using MLStyle

rmlines = JuliaVariables.rmlines
JuliaVariables.@quick_lambda begin
@testset "JuliaVariables.jl" begin
    func = solve(:(function f(x)
        let y = x + 1
            y
        end
    end))
    println(func |> rmlines)

    func = solve(:(function f(x)
        y = x + 1
        z -> z + y
    end))

    println(func |> rmlines)

    func = solve(:(function f(x)
        y = x + 1
        let y = y + 1
            (x + y  + z for z in 1:10)
        end
    end))
    println(func |> rmlines)

    func = solve(
        macroexpand(@__MODULE__,:(
        @inline function f(x)
            y = x + 1
            let y = y + 1
                (x + y  + z for z in 1:10)
            end
        end
    )))
    println(func |> rmlines)
    # @test map(haskey(func.scope.freevars, _), [])

    func = solve_from_local(quote
        a -> a + z
    end)
    println(func |> rmlines)

    func = solve_from_local(quote
    q = a
    function z(x, k=1)
        (x=x, k=k, q=q)
    end
    end)
    println(func |> rmlines)
    # Write your own tests here.

    a = solve(:(function k(x::T) where T
                  :(1, T)
              end
    ))
    @test haskey(a.scope.bounds, :T)

    a = solve(:(2 .^ [2, 3]))
    @test eval(a) == [4, 8]

    a = solve(:(function z(x, k=1)
                   x + 20 + a + k
               end
    ))
    @test haskey(a.scope.bounds, :k)

    a = solve(:(function z(x, k=1)
                   (k=k, )
               end
    ))
    @test haskey(a.scope.bounds, :k)

    @test @when :(k=$_, )  = a.func.args[2].args[2] begin
        true
    @otherwise
        false
    end

    func = solve(:(function (x)
              z = x + 1
              y -> begin
                  z += 1
                  z
              end
          end))

    println(func |> rmlines)
    @test haskey(func.scope.bounds, :z)
    @test func.scope.bounds[:z].is_mutable.x

end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM, SparseArrays, LinearAlgebra, Test

""" Calculate mortar projection matrix P = D^-1*M from mortar assembly. """
function calculate_mortar_projection_matrix(problem::Problem{Mortar}, ndim::Int)

    C1 = sparse(problem.assembly.C1, ndim, ndim)
    C2 = sparse(problem.assembly.C2, ndim, ndim)

    @debug("problem.assembly", problem.assembly.K, problem.assembly.D,
            problem.assembly.Kg, problem.assembly.fg, problem.assembly.f,
            problem.assembly.g)
    @assert isempty(problem.assembly.K)
    @assert isempty(problem.assembly.D)
    @assert isempty(problem.assembly.Kg)
    @assert isempty(problem.assembly.fg)
    @assert isempty(problem.assembly.f)

    @assert C1 == C2
    @assert problem.properties.adjust == false

    S = get_nonzero_rows(C2)
    M = setdiff(get_nonzero_columns(C2), S)

    # Construct matrix P = D^-1*M
    D_ = C2[S,S]
    M_ = -C2[S,M]

    P = lu(D_) \ Matrix(M_)

    return S, M, P
end

# two linear element clipping, calculation of projection matrix P for
# standard and dual basis

X = Dict(
    1 => [0.0, 0.0, 0.0],
    2 => [1.0, 0.0, 0.0],
    3 => [0.0, 1.0, 0.0],
    4 => [-0.25, 0.50, 0.00],
    5 => [0.50, -0.25, 0.00],
    6 => [0.75, 0.75, 0.00])

slave = Element(Tri3, (1, 2, 3))
master = Element(Tri3, (4, 5, 6))
update!((slave, master), "geometry", X)
update!(slave, "master elements", [master])
problem = Problem(Mortar, "two elements", 1, "temperature")
problem.properties.dual_basis = false
initialize!(problem, 0.0)
assemble!(problem, 0.0)
# forget to add elements to problem
add_elements!(problem, slave, master)
initialize!(problem, 0.0)
assemble!(problem, 0.0)
C1 = sparse(problem.assembly.C1)
C2 = sparse(problem.assembly.C2)
D = sparse(problem.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(problem, 6)
@test S == [1, 2, 3]
@test M == [4, 5, 6]

# visually inspected to be ok result
P_expected = 1/15*[9 9 -3; -7 13 9; 13 -7 9]
@test isapprox(P, P_expected)
um = [7.5, 15.0, 22.5]
@test isapprox(P*um, [9.0, 23.0, 13.0])

empty!(problem.assembly)
problem.properties.dual_basis = true
assemble!(problem, 0.0)
C1 = sparse(problem.assembly.C1)
C2 = sparse(problem.assembly.C2)
D = sparse(problem.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(problem, 6)
@test S == [1, 2, 3]
@test M == [4, 5, 6]
@test isapprox(P, P_expected)
um = [7.5, 15.0, 22.5]
@test isapprox(P*um, [9.0, 23.0, 13.0])

# two quadratic element clipping, calculation of projection matrix P
# for standard basis

X = Dict(
    1 => [0.0, 0.0, 0.0],
    2 => [1.0, 0.0, 0.0],
    3 => [0.0, 1.0, 0.0],
    7 => [-0.25, 0.50, 0.00],
    8 => [0.50, -0.25, 0.00],
    9 => [0.75, 0.75, 0.00])
# middle nodes
X[4] = 1/2*(X[1] + X[2])
X[5] = 1/2*(X[2] + X[3])
X[6] = 1/2*(X[3] + X[1])
X[10] = 1/2*(X[7] + X[8])
X[11] = 1/2*(X[8] + X[9])
X[12] = 1/2*(X[9] + X[7])
slave = Element(Tri6, (1, 2, 3, 4, 5, 6))
master = Element(Tri6, (7, 8, 9, 10, 11, 12))
update!((slave, master), "geometry", X)
update!(slave, "master elements", [master])
problem = Problem(Mortar, "two elements", 1, "temperature")
problem.properties.dual_basis = false
problem.properties.alpha = 0.2
add_elements!(problem, slave, master)
initialize!(problem, 0.0)
assemble!(problem, 0.0)
C1 = sparse(problem.assembly.C1)
C2 = sparse(problem.assembly.C2)
D = sparse(problem.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(problem, 12)
@test S == [1, 2, 3, 4, 5, 6]
@test M == [7, 8, 9, 10, 11, 12]
@debug("Projection matrix P", Matrix(P))
# visually inspected to be ok result
P_expected = 1/675*[81 81 189 972 -324 -324; 609 429 81 -1092 1404 -756; 429 609 81 -1092 -756 1404; -39 231 -81 132 396 36; -81 -81 81 108 324 324; 231 -39 -81 132 36 396]
@test isapprox(P, P_expected)
um = 15/2*[1, 2, 3]
um = [um[1], um[2], um[3], 0.5*(um[1]+um[2]), 0.5*(um[2]+um[3]), 0.5*(um[3]+um[1])]
us = P*um
@test isapprox(us, [9.0, 23.0, 13.0, 16.0, 18.0, 11.0])

# two quadratic element clipping, calculation of projection matrix P for
# dual lagrange basis
problem.properties.dual_basis = true
problem.properties.alpha = 0.2
empty!(problem.assembly)
initialize!(problem, 0.0)
assemble!(problem, 0.0)
C1 = sparse(problem.assembly.C1)
C2 = sparse(problem.assembly.C2)
D = sparse(problem.assembly.D)
@test length(D) == 0
@test C1 == C2
S, M, P = calculate_mortar_projection_matrix(problem, 12)
P_expected = 1/675*[81 81 189 972 -324 -324; 609 429 81 -1092 1404 -756; 429 609 81 -1092 -756 1404; -39 231 -81 132 396 36; -81 -81 81 108 324 324; 231 -39 -81 132 36 396]
@test isapprox(P, P_expected)
um = 15/2*[1, 2, 3]
um = [um[1], um[2], um[3], 0.5*(um[1]+um[2]), 0.5*(um[2]+um[3]), 0.5*(um[3]+um[1])]
us = P*um
@test isapprox(us, [9.0, 23.0, 13.0, 16.0, 18.0, 11.0])