text:
# The below tests are network-dependent, and actually make calls to GitHub's
# API. They're all read-only, meaning none of them require authentication.

testuser = Owner("julia-github-test-bot")
julweb = Owner("JuliaWeb", true)
ghjl = Repo("JuliaWeb/GitHub.jl")
testcommit = Commit("3a90e7d64d6184b877f800570155c502b1119c15")
testuser_sshkey = "ssh-rsa AAAAB3NzaC1yc2EAAAADAQABAAABAQDVDBxFza4BmQTCTFeTyK"*
    "3xT+T98dmiMWXC2lM/esw3MCRHg7cynLWr/jUgjs72DO2nqlCTKI88yd2gcbW5/pBP6NVumc"*
    "pM7eJzZJ3TKKwdGUD49nahncg5imHZUQbCqtQbAYEj+uFfqa9QNm6NkZdAdPdB6dJG2+QEuk"*
    "rIGWKsmihP7vGzRLdebGwng2aNUfdAyVwq5Af4g5qfyRT9MtOXXM/7tDAVfC/g4QjkQ52giG"*
    "3FRqehMHOfl4iw9cYggJ3owr+T/RhwBHhE9G+sIaq4cEjRxogf65xzJfRtxxM2RBYDM9GMyX"*
    "6s2dFghew6MMc2x7OJM30W+OedhtZuk3Xp"

hasghobj(obj, items) = any(x -> name(x) == name(obj), items)

# This token has public, read-only access, and is required so that our
# tests don't get rate-limited. The only way a malicious party could do harm
# with this token is if they used it to abuse the rate limit associated with
# the token (not too big of a deal). The token is hard-coded in an obsfucated
# manner in an attempt to thwart token-stealing crawlers.
auth = authenticate(string(circshift(["bcc", "3fc", "03a", "33e",
                                      "c09", "363", "5f1", "bd3",
                                      "fc6", "77b", '5', "9cf",
                                      "868", "033"], 3)...))

@test rate_limit(; auth = auth)["rate"]["limit"] == 5000

@testset "Owners" begin
    # test GitHub.owner
    @test name(owner(testuser; auth = auth)) == name(testuser)
    @test name(owner(julweb; auth = auth)) == name(julweb)

    # test GitHub.orgs
    @test hasghobj("JuliaWeb", first(orgs("jrevels"; auth = auth)))

    # test GitHub.followers, GitHub.following
    @test hasghobj("jrevels", first(followers(testuser; auth = auth)))
    @test hasghobj("jrevels", first(following(testuser; auth = auth)))

    # test GitHub.repos
    @test hasghobj(ghjl, first(repos(julweb; auth = auth)))

    # test sshkey/gpgkey retrieval
    @test_broken GitHub.sshkeys(testuser; auth = auth)[1][1]["key"] == testuser_sshkey
    @test startswith(GitHub.gpgkeys("JuliaTagBot"; auth = auth)[1][1]["raw_key"],
                     "-----BEGIN PGP PUBLIC KEY BLOCK-----")

    # test membership queries
    @test GitHub.check_membership(julweb, testuser; auth = auth)
    @test !GitHub.check_membership("JuliaLang", testuser; auth = auth, public_only=true)

    @test GitHub.isorg(julweb)
    @test !GitHub.isorg(testuser)
end

@testset "Repositories" begin
    # test GitHub.repo
    @test name(repo(ghjl; auth = auth)) == name(ghjl)

    # test GitHub.forks
    @test length(first(forks(ghjl; auth = auth))) > 0

    # test GitHub.contributors
    @test hasghobj("jrevels", map(x->x["contributor"], first(contributors(ghjl; auth = auth))))

    # test GitHub.stats
    @test stats(ghjl, "contributors"; auth = auth).status < 300

    # test GitHub.branch, GitHub.branches
    @test name(branch(ghjl, "master"; auth = auth)) == "master"
    @test hasghobj("master", first(branches(ghjl; auth = auth)))

    # test GitHub.commit, GitHub.commits
    @test name(commit(ghjl, testcommit; auth = auth)) == name(testcommit)
    @test hasghobj(testcommit, first(commits(ghjl; auth = auth)))

    # test GitHub.file, GitHub.directory, GitHub.readme, GitHub.permalink
    readme_file = file(ghjl, "README.md"; auth = auth)
    src_dir = first(directory(ghjl, "src"; auth = auth))
    owners_dir = src_dir[findfirst(c -> c.path == "src/owners", src_dir)]
    test_sha = "eab14e1ab7b4de848ef6390101b6d40b489d5d08"
    readme_permalink = string(permalink(readme_file, test_sha))
    owners_permalink = string(permalink(owners_dir, test_sha))
    @test readme_permalink == "https://github.com/JuliaWeb/GitHub.jl/blob/$(test_sha)/README.md"
    @test owners_permalink == "https://github.com/JuliaWeb/GitHub.jl/tree/$(test_sha)/src/owners"
    @test readme_file == readme(ghjl; auth = auth)
    @test hasghobj("src/GitHub.jl", src_dir)

    # test GitHub.status, GitHub.statuses
    @test status(ghjl, testcommit; auth = auth).sha == name(testcommit)
    @test !(isempty(first(statuses(ghjl, testcommit; auth = auth))))

    # test GitHub.comment, GitHub.comments
    @test name(comment(ghjl, 154431956; auth = auth)) == 154431956
    @test !(isempty(first(comments(ghjl, 40; auth = auth))))

    # These require `auth` to have push-access (it's currently a read-only token)
    # @test hasghobj("jrevels", first(collaborators(ghjl; auth = auth)))
    # @test iscollaborator(ghjl, "jrevels"; auth = auth)
end

@testset "Issues" begin
    state_param = Dict("state" => "all")

    # test GitHub.pull_request, GitHub.pull_requests
    @test pull_request(ghjl, 37; auth = auth).title == "Fix dep warnings"
    @test hasghobj(37, first(pull_requests(ghjl; auth = auth, params = state_param)))

    # test GitHub.issue, GitHub.issues
    @test issue(ghjl, 40; auth = auth).title == "Needs test"
    @test hasghobj(40, first(issues(ghjl; auth = auth, params = state_param)))
end

@testset "Gists" begin
    kc_gists, page_data = gists("KristofferC"; page_limit=1, params=Dict("per_page" => 5), auth = auth)
    @test typeof(kc_gists) == Vector{Gist}
    @test length(kc_gists) != 0
    @test kc_gists[1].owner.login == "KristofferC"

    gist_obj = gist("0cb70f50a28d79905aae907e12cbe58e"; auth = auth)
    @test length(gist_obj.files) == 2
    @test gist_obj.files["file1.jl"]["content"] == "Hello World!"
end

@testset "Reviews" begin
    pr = pull_request(ghjl, 59; auth = auth)
    review = first(reviews(ghjl, pr; auth=auth)[1])

    @test review.state == "CHANGES_REQUESTED"
end

@testset "Activity" begin
    # test GitHub.stargazers, GitHub.starred
    @test length(first(stargazers(ghjl; auth = auth))) > 10 # every package should fail tests if it's not popular enough :p
    @test hasghobj(ghjl, first(starred(testuser; auth = auth)))

    # test GitHub.watched, GitHub.watched
    @test hasghobj(testuser, first(watchers(ghjl; auth = auth)))
    @test hasghobj(ghjl, first(watched(testuser; auth = auth)))
end

testbot_key =
      "LS0tLS1CRUdJTiBSU0EgUFJJVkFURSBLRVktLS0tLQpNSUlFb3dJQkFBS0NBUUVBcFB0d1ZhWWVRVUtWcDFv"*
      "cDMzUnhxWVhUM0dPTFZmOHIvN25SejZkVDVMVzBoUUFXCjZyVEN2NzV0RDdueU5wTGtrdmZGcng2QzlaTlg4"*
      "NTFZSEhMZ3hQNHhFbTcrb2lBajlqNjNwRkxybDNqUEd6RDAKQmVVNkcyVitpNGc5c00rMnpvR3RlbFpQaW4w"*
      "RHhhZi9TSThCL2ZxcmJ4SVVVMjN3T0hDZDU3L0Q2WUVQWFFOagpWTzczb0xXb0ltNTkxc3VLNDZOcTdMNkhq"*
      "SW85M0R6eWNWUExyRDNiK0xTSXlZQXMzTnhPTXNuNTZ5S0JLaVp1CktHMldnVjJOWlBOUXlZYVdCQWtKZzB6"*
      "cFBDUVdXcStZYkh4QS8rRHIrNnRsSXFnS3QyTWFlNGJKZ1FtTEgxSzkKNjRnR0RuK04zQTlmQ04rRjlLUndD"*
      "a3dsek1JeCtGSFhuMnh6NHdJREFRQUJBb0lCQUZBVXFDTTZ2ZWJKYmlETQpXQlRaemE0T2dwYXdUdHJRUUVn"*
      "aHB5RFhSSlo0U0laaVU3MWJQa2lhSXhBR0h2YTBlSFNLQmcrSmpPR3N3bjFDCmU1bGJmWlRCR1lxc0M3Y2dT"*
      "TFJuSHZhSk5VZlI0UjErdG9RQ0R5RnJNM0NCRWdpMzJqRUVxdkw5NENBRnJJWU4KbEpGZ3NRUFozOHJMQ01p"*
      "eVRXN240dzJ4ZkdQamRLcitlbTc4R0syUU1ad0l5WllQTG90V2dnc0QyY1lQV1BsdQpyQ2VRSmt2blVKMnhS"*
      "RTdnWUQ5YkUzdWFUN0xiUFhYZldwYW81aVZFdjFROCtJZjFON2wveTlkODBhY1BiUDNHCkFZY0NqeWprN1NV"*
      "bEhKYVNlYXNzbFBpT1FiOS9xNUh3aVljQTZRbnZCSlkvTjdVWGJQVXVOZEl0U2NmdUZQQUwKKyt3T0RtRUNn"*
      "WUVBMWJHZHlnT0NOZzNWUEVlNjRuSzlvTklGV043ZlE2SngvbUZIb3NUUHV3d3BMaFlLb3htMgpFMTl3Mzlh"*
      "TDROYWJ3TWZnVGIyRDJBbXVpUTJkaXRQNHRjc1cxOUZXS09peWhWdXgrSFl2eWt4TTZwRFVZa2Q1CitEY1lS"*
      "NFd1bVpjcUZSRlc3YmN4SXA2eWdseHgwNHhuMmg5U3Bnb0dwRDBFY2ZQOEFpbFl1NUVDZ1lFQXhhVUgKdEhW"*
      "cXk1WnRPRElEbkUwclJtQzdlUjlzdGR3eUZwdlJ0RWcvenBJZEIwVHUyVFFzZ3djRGc3N0ZLRUdvVEpJcQpr"*
      "WFFBdXJVOWRpUXkyVTVZT2J6R3lvSG9JTFpTc0FhV2JjOU10amVNaFdGVFhSNnNrY2loS0NSd0Z3UTNscyty"*
      "CnRtUkVXTGtDUm0rOVBJVjJxalg2SG9kNktWT2NkV0FEYldnR3RqTUNnWUJEL0JnSkZ3aXNEY2FUUVBiUjZG"*
      "TXcKQU9FMm51RkU4VDkzQUpmN3pzV1A2cFNIVnZmWFgreXZTU1B0OHFIWnpDME5MZ25NY2NpcVNKcEFmQlpz"*
      "L25jWAp6eDdiVm53azA3TkgvaDRtditNQVp6bnBQbDV6VGU0ZDY5bExsOW91ZndzaVhMdmRNUFR1NExKR0N3"*
      "Y2ptSDNKCnhVRlVGY2g3SS9ad0VvRlFacnNXSVFLQmdRQ2pFWjVoemQ3blNwMmlsK0ZTdkhqUUFFK3RoN2Z4"*
      "OUZOL1ErQ3AKbGxMTzVNNytpR2xvM0JzOW9EUE9KMEFVRHRnRkZUUDUvblA3bUQyMWsvaEFRdHVZQjZFY3hF"*
      "SDVlM1NOdDJHMgpDQ3VLekJvc2tqaHR4RGt0cnhNSVE4Z1h0V3NJQ3gvcHhLQi9jMlhsSjV4Q3F2dFZSR094"*
      "ZktYV0l4NGIyYlA2Ck9MSVE4UUtCZ0Z5cERKdXczYjFZVEpqbmI0UTZQR0FaUzNtVDZ1SVVrbGtmWENuNWo4"*
      "SUNFMXpZNVl1RDRpbE4KR3RkWHlqeitiMVg2UEZnVFhpdzQxU2xaUXMrVzhZUXVQSHA2QWFsV3VBN1l6cFNp"*
      "eWd3eFlaaWFIODlGWUxuagpxdm1lMDFIeUtIdnRQOE95YmI0UGsvZ3V0aXlaV0U0RW1JVWZEZ2I4RlY4VU5z"*
      "MW1IcklTCi0tLS0tRU5EIFJTQSBQUklWQVRFIEtFWS0tLS0tCg=="

@testset "Apps" begin
    @test app(4123; auth=auth).name == "femtocleaner"
    @test app("femtocleaner"; auth=auth).name == "femtocleaner"

    key = MbedTLS.PKContext()
    MbedTLS.parse_key!(key, base64decode(testbot_key))
    jwt = GitHub.JWTAuth(4484, key)
    @test app(; auth=jwt).name == "juliawebtestbot"

    @test length(installations(jwt)[1]) == 1
end

@testset "Git Data" begin
    github_jl = Repo("JuliaWeb/GitHub.jl")

    g = gitcommit(github_jl, "0d9f04ce4be061d3c2b12644316a232c8f889b44"; auth=auth)
    @test g.tree["sha"] == "e22fee36cb13d9a1850b242f79938458221a5d2e"

    t = tree(github_jl, g.tree["sha"]; auth=auth)
    for entry in t.tree
        if entry["path"] == "README.md"
            @test entry["sha"] == "95c8d1aa2a7b1e6d672e15b67e0df4abbe57dcbe"
            @test entry["type"] == "blob"

            b = blob(github_jl, entry["sha"]; auth=auth)
            @test occursin("GitHub.jl", String(base64decode(replace(b.content,"\n" => ""))))

            break
        end
    end
end

@testset "Tags and References" begin
    # All tags in this repo are lightweight tags which are not covered by the API
    # Maybe test in the future when we have a use case
    github_jl = Repo("JuliaWeb/GitHub.jl")
    ref = reference(github_jl, "heads/master"; auth=auth)
    @test ref.object["type"] == "commit"

    # Tag API
    reponame = "QuantEcon/Expectations.jl"
    version = "v1.0.1"
    exptag = tag(reponame, version; auth=auth)
    @test isa(exptag, Tag)
end


text:
module Tokenize
export TokenBuilder, Token, TokenListtokenize, tokenize

mutable struct TokenBuilder
  typ::Symbol
  start_index::Int
  end_index::Int
  token::Vector{Char}
  incomplete::Bool

  TokenBuilder() = new(:none, 1, 1, [], false)
end

struct Token
  token::String
  typ::Symbol
  index::Int
  incomplete::Bool

  function Token(tb::TokenBuilder)
    token = String(tb.token)
    new(token, tb.typ, tb.start_index, tb.incomplete)
  end

  function Token(t::Token;
     token=t.token, typ=t.typ, index=t.index, incomplete=t.incomplete)

    new(token, typ, index, incomplete)
  end
end

struct TokenList
  tokens::Vector{Token}
  i::Int

  TokenList(str::String) = new(tokenize(str), 1)
  function TokenList(tl::TokenList, dlm::String)
    new(view(tl.tokens, tl.i:(findnext(tl.tokens, dlm, tl.i)-1)), 1)
  end
end
peek_token(tl::TokenList) = (t=tl.tokens[tl.i]; (t.token, t.typ, t.index))
get_token!(tl::TokenList) = (t=peek_token(tl); tl.i+=1; t)

Base.done(tl::TokenList)   = tl.i > length(tl.tokens)
Base.length(tl::TokenList) = length(tl.tokens)

function tokenize(line; incomplete=:notreally)
  tokens = Token[]
  escape = false
  token_builder = TokenBuilder()
  next_index = 1

  if incomplete in (:single_quotes, :double_quotes)
    token_builder.typ = incomplete
    token_builder.incomplete = true
  end

  function build(completes=false)
    completes && ( token_builder.incomplete = false )
    !istype(:none) && push!(tokens, Token(token_builder))
    token_builder.typ   = :none
    token_builder.token = Char[]
    token_builder.incomplete = false
  end

  function push_token(c)
    push!(token_builder.token, c)
  end

  function last_char()
    token_builder.token[end]
  end

  istype(t, ts...) = token_builder.typ == t || istype(ts...)
  istype() = false

  function newtype(t)
    if istype(:none)
      token_builder.start_index = token_builder.end_index
    end
    token_builder.typ = t
  end

  function next_char(c)
    token_builder.end_index += sizeof(string(c))
  end

  function isincomplete()
    token_builder.incomplete = true
  end

  next_index = 1
  for c in line
    @label loop_start

    if escape
      push!(token_builder.token, c)
      istype(:none) && newtype(:normal)
      escape = false
    elseif istype(:ambiguous_special)
      if c == last_char()
        push_token(c)
        newtype(:special)
        build()
      else
        newtype(:special)
        build()
        @goto loop_start
      end
    elseif c == '\\' && !istype(:single_quotes)
      escape = true
    elseif c == '@' && istype(:none)
      newtype(:macro)
      push_token(c)
    elseif c == '@' && istype(:macro)
      newtype(:macromma)
      push_token(c)
      build()
    elseif c == '"'
      istype(:single_quotes) ? push_token(c) :
      istype(:double_quotes) ? build(true)   :
      (build(); newtype(:double_quotes); isincomplete())
    elseif c == '\''
      istype(:double_quotes) ? push_token(c) :
      istype(:single_quotes) ? build(true)   :
      (build(); newtype(:single_quotes); isincomplete())
    elseif c == ' ' && !istype(:single_quotes, :double_quotes)
      build()
    elseif c in "(){}[]<>;\n" && !istype(:single_quotes, :double_quotes)
      build()
      newtype(:special)
      push_token(c)
      build()
    elseif c in "&|" && !istype(:single_quotes, :double_quotes)
      build()
      newtype(:ambiguous_special)
      push_token(c)

    elseif istype(:none)
      newtype(:normal)
      push_token(c)
    else
      push_token(c)
    end
    next_char(c)
  end
  build()

  return tokens
end

end

ext: jl

max_stars_repo_name: fcard/Fxcrd

max_issues_repo_name: fcard/Fxcrd

max_forks_repo_name: fcard/Fxcrd


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

using JuliaFEM, Test

# test tet4 stiffness matrix

X = Dict(1 => [2.0, 3.0, 4.0],
         2 => [6.0, 3.0, 2.0],
         3 => [2.0, 5.0, 1.0],
         4 => [4.0, 3.0, 6.0])
element = Element(Tet4, (1, 2, 3, 4))
update!(element, "youngs modulus", 96.0)
update!(element, "poissons ratio", 1/3)
update!(element, "geometry", X)
problem = Problem(Elasticity, "tet4", 3)
add_element!(problem, element)
time = 0.0
assemble!(problem, time)
K_expected = [
     149   108   24   -1    6   12  -54   -48    0  -94   -66  -36
     108   344   54  -24  104   42  -24  -216  -12  -60  -232  -84
      24    54  113    0   30   35    0   -24  -54  -24   -60  -94
      -1   -24    0   29  -18  -12  -18    24    0  -10    18   12
       6   104   30  -18   44   18   12   -72  -12    0   -76  -36
      12    42   35  -12   18   29    0   -24  -18    0   -36  -46
     -54   -24    0  -18   12    0   36     0    0   36    12    0
     -48  -216  -24   24  -72  -24    0   144    0   24   144   48
       0   -12  -54    0  -12  -18    0     0   36    0    24   36
     -94   -60  -24  -10    0    0   36    24    0   68    36   24
     -66  -232  -60   18  -76  -36   12   144   24   36   164   72
     -36   -84  -94   12  -36  -46    0    48   36   24    72  104]
@test isapprox(problem.assembly.K, K_expected)


text:
@doc """ index array data with the first column
""" ->
function index(dt)
    dict = Dict{ASCIIString,Array{ASCIIString,1}}()
    for i = 1:size(dt,1)
        dict[dt[i,1]] = dt[i,2:end]
    end
    dict
end

@doc """ left join data on the first column, assume have the same ids
And data have no header part
""" ->
function leftjoindata(d1,d2)
    d1_ind = index(d1)
    d2_ind = index(d2)
    @assert issubset(Set(collect(keys(d1_ind))),Set(collect(keys(d2_ind)))) == true
    nrow = size(d1,1)
    ncol = size(d1,2)-1 + size(d2,2)-1 + 1
    data = Array{ASCIIString,2}(nrow,ncol)
    for (i,id) in enumerate(keys(d1_ind))
        data[i,:] = vcat(id,d1_ind[id], d2_ind[id])'
    end
    data
end

function leftjoindata(d1,d2,d3...)
    #left join
    data = leftjoindata(d1,d2)
    for d in d3
        data = leftjoindata(data,d)
    end
    data
end


text:
"""
    sem_fit(model::AbstractSem; start_val = start_val, kwargs...)

Return the fitted `model`.

# Arguments
- `model`: `AbstractSem` to fit
- `start_val`: vector of starting values or function to compute starting values (1)
- `kwargs...`: keyword arguments, passed to starting value functions

(1) available options are `start_fabin3`, `start_simple` and `start_partable`.
For more information, we refer to the individual documentations and the online documentation on [Starting values](@ref).

# Examples
```julia
sem_fit(
    my_model;
    start_val = start_simple,
    start_covariances_latent = 0.5)
```
"""
function sem_fit end


text:
### Indexing

## Helper functions

"""
    ind2run(rle::RLEVector, i::Integer) -> run_index
    ind2run(rle::RLEVector, i::UnitRange) -> range_of_first_to_last_overlap
    ind2run(rle::RLEVector, i::AbstractArray) -> vector_of_runs

Locate runs, get index of run corresponding to the i-th value in the expanded runs
"""
function ind2run(rle::RLEVector, i::Integer)
  re = rle.runends
  n = length(re)
  0 < i <= last(re) || throw(BoundsError())
  run = searchsortedfirst(re,i,1,n)
  run
end

function ind2run(rle::RLEVector,i::UnitRange)
    re = rle.runends
    n = length(re)
    left_run = searchsortedfirst(re,first(i),1,n)
    right_run = searchsortedfirst(re,last(i),left_run,n)
    right_run <= n || throw(BoundsError())  # Can't be < 1
    left_run:right_run
end

function ind2run(rle::RLEVector, i::AbstractArray)
  re = rle.runends
  n = length(re)
  runs = searchsortedfirst(re,i,1,n)
  maximum(runs) <= n || throw(BoundsError())  # Can't be < 1
  runs
end

"""
    ind2runcontext(rle::RLEVector, i::Integer) -> (runindex, index_in_run, values_in_run_after_i)
    ind2runcontext(rle::RLEVector, i::UnitRange) -> (runindex, index_in_run, values_in_run_after_i)

Index of the run corresponding to the i'th value in the expanded runs, index in run and remainder of run
"""
function ind2runcontext(rle::RLEVector, i::Integer)
  run = ind2run(rle, i)
  runend = rle.runends[run]
  ind_in_run = run == 1 ? i : i - rle.runends[run - 1]
  (run, ind_in_run, runend - i)
end

function ind2runcontext(rle::RLEVector, i::UnitRange)
  s = first(i)
  e = last(i)
  n = length(rle)
  runs = ind2run(rle, i)
  left_run = first(runs)
  right_run = last(runs)
  runend = rle.runends[right_run]
  ind_in_run = left_run == 1 ? s : s - rle.runends[left_run - 1]
  (left_run, right_run, ind_in_run, runend - e)
end

## Just enough for AbstractArray
Base.IndexStyle(::Type{<:RLEVector}) = IndexLinear()
Base.firstindex(rle::RLEVector) = 1
Base.lastindex(rle::RLEVector) = length(rle)

function Base.getindex(rle::RLEVector, i::Int)
    run = ind2run(rle,i)
    rle.runvalues[run]
end

function Base.setindex!(rle::RLEVector, value, i::Integer)
  run = ind2run(rle,i)
  runvalue = rle.runvalues[run]
  runend = rle.runends[run]
  #value == runvalue && return rle # replace with same value, no-op
  isequal(value, runvalue) && return rle
  previous_run = run - 1
  next_run = run + 1
  at_start_of_run = (previous_run > 0 && i == rle.runends[previous_run] + 1) || i == 1
  at_end_of_run = i == runend
  match_left = run > 1 && isequal(rle.runvalues[previous_run], value)
  match_right = run < nrun(rle) && isequal(rle.runvalues[next_run], value)
  if at_end_of_run
    if at_start_of_run # in a run of length 1
      if match_right && match_left
        deleteat!(rle.runvalues, previous_run:run)
        deleteat!(rle.runends, previous_run:run)
      elseif match_right
        deleteat!(rle.runvalues,run)
        deleteat!(rle.runends,run)
      elseif match_left
        deleteat!(rle.runvalues,run)
        deleteat!(rle.runends,previous_run)
      else
        rle.runvalues[run] = value
      end
    else # end of a run longer than 1
      if match_right
        rle.runends[run] = runend - 1
      else
        insert!(rle.runvalues, next_run, value)
        insert!(rle.runends, run, runend - 1)
      end
    end
  elseif at_start_of_run
    if match_left
      rle.runends[previous_run] = rle.runends[previous_run] + 1
    else
      insert!(rle.runvalues, run, value)
      insert!(rle.runends, run, i)
    end
  else # middle of a run, average case
      rle.runends[run] = i - 1
      next_run = run + 1
      growat!(rle.runvalues, next_run, 2)
      growat!(rle.runends, next_run, 2)
      rle.runvalues[next_run] = value
      rle.runends[next_run] = i
      next_run = next_run + 1
      rle.runvalues[next_run] =  runvalue
      rle.runends[next_run] = runend
  end
  rle
end

## Indexing optimizations
function Base.getindex(x::RLEVector, indices::UnitRange)
    (left_run, right_run, index_in_left_run, run_remainder_right) = ind2runcontext(x,indices)
    n_run = (right_run - left_run) + 1
    f = first(indices)
    v = similar(x.runvalues, n_run)
    e = similar(x.runends, n_run)
    offset = f - 1
    i = 1
    @inbounds while left_run <= right_run
        v[i] = x.runvalues[left_run]
        e[i] = x.runends[left_run] - offset
        left_run = left_run + 1
        i = i + 1
    end
    e[end] = e[end] - run_remainder_right
    RLEVector{eltype(x),endtype(x)}(v, e)
end

function Base.getindex(x::RLEVector, i::AbstractVector{<:Int})
    run_indices = ind2run(x, i)
    RLEVector( x.runvalues[ run_indices ] )
end

function Base.setindex!(x::RLEVector, value::AbstractVector, indices::UnitRange)
    setindex!(x, RLEVector(value), indices)
end

function Base.setindex!(x::RLEVector, value::RLEVector, indices::UnitRange)
    length(value) != length(indices) && throw(BoundsError())
    i_left = first(indices)
    i_right = last(indices)
    if i_left == i_right
        return(setindex!(x,ins,i_left))
    end
    nrun_x = nrun(x)
    nrun_value = nrun(value)
    (run_left, run_right, index_in_run_left, run_remainder_right) = ind2runcontext(x,indices)
    # Move run markers to denote parts of original data that will be kept, accomodating completely filled runs or adjacent matches
    # We will keep 1:run_left and run_right:end and fill in the middle with value
    # FIXME: factor out these two expressions for something like ind2insertcontext
    fix_partial_run_left = false
    if x.runvalues[run_left] == first(value)
        run_left = run_left - 1
    elseif index_in_run_left == 1
        run_left = run_left - 1
        if run_left > 0 && first(value) == x.runvalues[run_left]
            run_left = run_left - 1
        end
    else
        fix_partial_run_left = true
    end
    if x.runvalues[run_right] == last(value)
        nrun_value = nrun_value - 1
    elseif run_remainder_right == 0
        run_right = run_right + 1
        if run_right <= nrun_x && last(value) == x.runvalues[run_right]
            nrun_value = nrun_value - 1
        end
    end
    nrun_out = run_left + nrun_value + ((nrun_x - run_right) + 1)
    nrun_diff = nrun_out - nrun_x
    # Resize and move
    if nrun_diff > 0
        growat!(x, run_right, nrun_diff)
    elseif nrun_diff < 0
        delete_range = (run_left + 1):(run_left - nrun_diff)
        deleteat!(x.runvalues, delete_range)
        deleteat!(x.runends, delete_range)
    end
    if fix_partial_run_left
        x.runends[run_left] = i_left - 1
    end
    # Insert incoming values
    value_runvalues = value.runvalues
    value_runends = value.runends
    bump = i_left - 1
    @inbounds for i in 1:nrun_value
        il = i + run_left
        x.runvalues[il] = value_runvalues[i]
        x.runends[il] = value_runends[i] + bump
    end
    x
end

## Iterators
# Iterator for ranges based on RLE e.g. (value, start:end)
struct RLEEachRangeIterator{T1,T2}
    rle::RLEVector{T1,T2}
end
eachrange(x::RLEVector) = RLEEachRangeIterator(x)
length(x::RLEEachRangeIterator) = nrun(x.rle)
eltype(::RLEEachRangeIterator{T1,T2}) where {T1,T2} = Tuple{T1,UnitRange{T2}}

function iterate(x::RLEEachRangeIterator, state = 1)
    state > nrun(x.rle) && return nothing
    newstate = state + 1
    first = starts(x.rle,state)
    last = ends(x.rle,state)
    ( (_values(x.rle)[state], first:last ), newstate )
end

## New style iterator
## state is (run_index, overall_index)
## We will assume that this iterate method is the only source of `state` and
## the only invalid state possible is one past the end
eltype(::RLEVector{T1,T2}) where {T1,T2} = T1
function iterate(x::RLEVector, state = (1,1))
    (run_index, overall_index) = state
    if overall_index > length(x) # run_index > nrun(x) || overall_index > length(x)
        out = nothing
    else
        this_end = ends(x)[run_index]
        if overall_index > this_end
            run_index = run_index + 1
        end
        out = (values(x)[run_index], (run_index, overall_index + 1))
    end
    out
end

# FIXME: add reverse iterator, see https://docs.julialang.org/en/latest/manual/interfaces/#man-interface-array-1

"""
    tapply(data_vector, rle, function)
    tapply(data_vector, factor_vector, function)

Map a function to blocks of vector, like `tapply` in R. The first and second argument must be of the same
    length. For the case of a standard vector as the second argument, this vector need not be sorted.

## Examples
    factor = repeat( ["a","b","c","d","e"], inner=20 )
    rle = RLEVector( factor )
    x = collect(1:100)
    tapply( x, factor, mean )
    tapply( x, rle, mean )
"""
function tapply(x::Vector, rle::RLEVector, fun::Function)
    length(x) == length(rle) || throw(ArgumentError("Arguments 'x' and 'rle' must have the same length."))
    allunique(rle.runvalues) || throw(ArgumentError("The values in an `RLEVector` must be unique when used as a factor with `tapply`."))
    Dict( (v,fun(x[r])) for (v,r) in eachrange(rle) ) # Note: using view on x[r] is not faster
end

function tapply(x::Vector{T1}, factor::Vector{T2}, fun::Function) where T1 where T2
    length(x) == length(factor) || throw(ArgumentError("Arguments 'x' and 'factor' must have the same length."))
    if ! issorted(factor)
        ind = sortperm(factor)
        y = x[ind]
        rle = RLEVector( factor[ind] )
    else
        rle = RLEVector( factor )
        y = x
    end
    Dict( (v,fun(x[r])) for (v,r) in eachrange(rle) )
end


text:
##
#    File:    milo1.jl
#
#    Purpose:  Demonstrates how to solve a small mixed
#              integer linear optimization problem using the MOSEK Python API.
##

using Mosek

# Define a stream printer to grab output from MOSEK
printstream(msg::String) = print(msg)

bkc = [ MSK_BK_UP, MSK_BK_LO  ]
blc = [      -Inf,      -4.0  ]
buc = [     250.0,       Inf  ]

bkx = [ MSK_BK_LO, MSK_BK_LO  ]
blx = [       0.0,       0.0  ]
bux = [       Inf,       Inf  ]

c   = [       1.0,      0.64 ]

A    = sparse( [ 1, 1, 2, 2],
               [ 1, 2, 1, 2],
               [ 50.0, 3.0, 31.0, -2.0])

numvar = length(bkx)
numcon = length(bkc)

# Create a task
maketask() do task
    # Attach a printer to the task
    putstreamfunc(task,MSK_STREAM_LOG,printstream)

    # Append 'numcon' empty constraints.
    # The constraints will initially have no bounds.
    appendcons(task,numcon)

    #Append 'numvar' variables.
    # The variables will initially be fixed at zero (x=0).
    appendvars(task,numvar)

    # Set the linear term c_j in the objective.
    putclist(task,[1:numvar;],c)

    # Set the bounds on variables
    # blx[j] <= x_j <= bux[j]
    putvarboundslice(task,1,numvar+1,bkx,blx,bux)

    # Input columns of A
    putacolslice(task,1,numvar+1, A.colptr[1:numvar],A.colptr[2:numvar+1],A.rowval,A.nzval)

    putconboundslice(task,1,numcon+1,bkc,blc,buc)

    # Input the objective sense (minimize/maximize)
    putobjsense(task,MSK_OBJECTIVE_SENSE_MAXIMIZE)

    # Define variables to be integers
    putvartypelist(task,[ 1, 2 ],
                   [ MSK_VAR_TYPE_INT, MSK_VAR_TYPE_INT ])

    # Optimize the task
    optimize(task)

    # Print a summary containing information
    # about the solution for debugging purposes
    solutionsummary(task,MSK_STREAM_MSG)

    prosta = getprosta(task,MSK_SOL_ITG)
    solsta = getsolsta(task,MSK_SOL_ITG)

    if solsta == MSK_SOL_STA_INTEGER_OPTIMAL || solsta == MSK_SOL_STA_NEAR_INTEGER_OPTIMAL
        # Output a solution
        xx = getxx(task,MSK_SOL_ITG)
        @printf("Optimal solution: %s\n", xx')
    elseif solsta == MSK_SOL_STA_UNKNOWN
        println("Unknown solution status")
    else
        println("Other solution status")
    end
end


text:
export stretching

"""

  stretching(ref,cur,t,window,fmin,fmax;dvmin,dvmax,ntrials)

This function compares the Reference waveform to stretched/compressed current
waveforms to get the relative seismic velocity variation (and associated error).
It also computes the correlation coefficient between the Reference waveform and
the current waveform.

# Arguments
- `ref::AbstractArray`: Reference correlation.
- `cur::AbstractArray`: Current correlation.
- `t::AbstractArray`: time vector, common to both `ref` and `cur`.
- `window::AbstractArray`: vector of the indices of the `cur` and `ref` windows
                          on which you want to do the measurements
- `fmin::Float64`: minimum frequency in the correlation [Hz]
- `fmax::Float64`: maximum frequency in the correlation [Hz]
- `dvmin::Float64`: minimum bound for the velocity variation; e.g. dvmin=-0.03
                   for -3% of relative velocity change
- `dvmax::Float64`: maximum bound for the velocity variation; e.g. dvmin=0.03
                  for 3% of relative velocity change
- `ntrial::Int`:  number of stretching coefficient between dvmin and dvmax, no need to be higher than 100

# Returns
- `dvv::Float64`: Relative Velocity Change dv/v (in %)
- `cc::Float64`: Correlation coefficient between the reference waveform and the
                      best stretched/compressed current waveform
- `cdp::Float64`: Correlation coefficient between the reference waveform and the
                 initial current waveform
- `ϵ::Array{Float64,1}`: Vector of Epsilon values (ϵ =-dt/t = dv/v)
- `err::Float64`: Errors in the dv/v measurements based on [Weaver et al., 2011](https://onlinelibrary.wiley.com/doi/full/10.1111/j.1365-246X.2011.05015.x)
- `allC::Array{Float64,1}`: Vector of the correlation coefficient between the
                        reference waveform and every stretched/compressed
                        current waveforms

This code is a Julia translation of the Python code from [Viens et al., 2018](https://github.com/lviens/2018_JGR).
"""
function stretching(ref::AbstractArray, cur::AbstractArray, t::AbstractArray,
                    window::AbstractArray, fmin::Float64, fmax::Float64;
                    dvmin::Float64=-0.1, dvmax::Float64=0.1, ntrial::Int=100)
    ϵ = range(dvmin, stop=dvmax, length=ntrial)
    L = 1. .+ ϵ
    tau = t * L'
    allC = zeros(ntrial)

    # set of stretched/compressed current waveforms
    waveform_ref = ref[window]
    for ii = 1:ntrial
        s = LinearInterpolation(tau[:,ii],cur,extrapolation_bc=Flat())(t)
        waveform_cur = s[window]
        allC[ii] = cor(waveform_ref,waveform_cur)
    end

    cdp = cor(cur[window],ref[window])
    # find the maximum correlation coefficient
    imax = argmax(allC)
    if imax >= ntrial -1
        imax = ntrial - 2
    end
    if imax <= 3
        imax += 3
    end
    # Proceed to the second step to get a more precise dv/v measurement
    dtfiner = Array(range(ϵ[imax-2],stop=ϵ[imax+1],length=500))
    # to get same result as scipy, use line below; this requires gcc
    # using Dierckx; etp = Spline1D(ϵ[range(imax-3,stop=imax+1)],allC[range(imax-3,stop=imax+1)])
    etp = CubicSplineInterpolation(ϵ[range(imax-3,stop=imax+1)],allC[range(imax-3,stop=imax+1)])
    CCfiner = etp(dtfiner)
    dvv = 100. * dtfiner[argmax(CCfiner)]
    cc = maximum(CCfiner) # Maximum correlation coefficient of the refined analysis

    # Error computation based on Weaver, R., C. Hadziioannou, E. Larose, and M.
    # Campillo (2011), On the precision of noise-correlation interferometry,
    # Geophys. J. Int., 185(3), 1384?1392
    T = 1 / (fmax - fmin)
    X = cc
    # extremely similar signals can return cc>1.0 (not possible), so we limit cc to 1.0 to prevent sqrt(neg)
    if X > 1.0
        X=1.0
    end

    wc = π * (fmin + fmax)

    # Coda window here is assumed to be either one-sided in positive or negative side, or symmetric in both sides.
    t1 = minimum(abs.(t[window]))
    t2 = maximum(abs.(t[window]))

    err = 100 * (sqrt(1-X^2)/(2*X)*sqrt((6*sqrt(π/2)*T)/(wc^2*(t2^3-t1^3))))

    return dvv,cc,cdp,Array(ϵ),err,allC
end

"""

    stretching(C,t,window,fmin,fmax;dvmin,dvmax,ntrials)

This function compares the Reference waveform to stretched/compressed current
waveforms to get the relative seismic velocity variation (and associated error)
for all correlations in CorrData `C`. It also computes the correlation
coefficient between the Reference waveform and the current waveform.

# Arguments
- `C::CorrData`: Correlation for dv/v.
- `cur::AbstractArray`: Current correlation.
- `t::AbstractArray`: time vector, common to both `ref` and `cur`.
- `window::AbstractArray`: vector of the indices of the `cur` and `ref` windows
                          on which you want to do the measurements
- `fmin::Float64`: minimum frequency in the correlation [Hz]
- `fmax::Float64`: maximum frequency in the correlation [Hz]
- `dvmin::Float64`: minimum bound for the velocity variation; e.g. dvmin=-0.03
                   for -3% of relative velocity change
- `dvmax::Float64`: maximum bound for the velocity variation; e.g. dvmin=0.03
                  for 3% of relative velocity change
- `ntrial::Int`:  number of stretching coefficient between dvmin and dvmax, no need to be higher than 100

# Returns
- `dvv::Array{Float64,1}`: Relative Velocity Change dv/v (in %)
- `cc::Array{Float64,1}`: Correlation coefficient between the reference waveform and the
                      best stretched/compressed current waveform
- `cdp::Array{Float64,1}`: Correlation coefficient between the reference waveform and the
                 initial current waveform
- `ϵ::Array{Float64,1}`: Vector of Epsilon values (ϵ =-dt/t = dv/v)
- `err::Array{Float64,1}`: Errors in the dv/v measurements based on [Weaver et al., 2011](https://onlinelibrary.wiley.com/doi/full/10.1111/j.1365-246X.2011.05015.x)
- `allC::Array{Float64,2}`: Matrix of the correlation coefficient between the
                        reference waveform and every stretched/compressed
                        current waveforms
"""
function stretching(C::CorrData, t::AbstractArray, window::AbstractArray,
                    fmin::Float64, fmax::Float64; dvmin::Float64=-0.1,
                    dvmax::Float64=0.1, ntrial::Int=100)

    N = length(C.t)
    ref = SeisNoise.stack(C,allstack=true)
    dvv = zeros(length(C.t))
    cc = similar(dv)
    cdp = similar(dv)
    err = similar(dv)
    allC = zeros(ntrial,N)
    lags = -ref.maxlag:1/ref.fs:ref.maxlag

    for ii = 1:N
        dvv[ii], cc[ii], cdp[ii], ϵ, err[ii], allC[:,ii] = stretching(ref.corr[:],
                                                            C.corr[:,ii],
                                                            lags,
                                                            window,
                                                            fmin,
                                                            fmax,
                                                            dvmin=dvmin,
                                                            dvmax=dvmax,
                                                            ntrial=ntrial)
    end
    ϵ = collect(range(dvmin, stop=dvmax, length=ntrial))

    return dvv, cc, cdp, ϵ, err, allC
end

"""

stretching(ref,cur,t,window,freqbands;dvmin,dvmax,ntrials,norm)

This function compares the Reference waveform to stretched/compressed current
waveforms at a range of frequencies to get the frequency-dependent relative seismic
velocity variation (and associated error). It also computes the correlation
coefficient between the Reference waveform and the current waveform.

# Arguments
- `ref::AbstractArray`: Reference correlation.
- `cur::AbstractArray`: Current correlation.
- `t::AbstractArray`: time vector, common to both `ref` and `cur`.
- `window::AbstractArray`: vector of the indices of the `cur` and `ref` windows
                          on which you want to do the measurements
- `freqbands::AbstractArray`: Frequency bands over which to compute dv/v
- `dvmin::Float64`: minimum bound for the velocity variation; e.g. dvmin=-0.03
                   for -3% of relative velocity change
- `dvmax::Float64`: maximum bound for the velocity variation; e.g. dvmin=0.03
                  for 3% of relative velocity change
- `ntrial::Int`:  number of stretching coefficient between dvmin and dvmax, no need to be higher than 100
- `norm::Bool`: Whether or not to normalize signals before dv/v

# Returns
- `freqbands::AbstractArray`: Array of frequencies where dv/v was measured
- `dvv::Array{Float64,1}`: Relative Velocity Change dv/v (in %)
- `cc::Array{Float64,1}`: Correlation coefficient between the reference waveform and the
                      best stretched/compressed current waveform
- `cdp::Array{Float64,1}`: Correlation coefficient between the reference waveform and the
                 initial current waveform
- `ϵ::Array{Float64,1}`: Vector of Epsilon values (ϵ =-dt/t = dv/v)
- `err::Array{Float64,1}`: Errors in the dv/v measurements based on [Weaver et al., 2011](https://onlinelibrary.wiley.com/doi/full/10.1111/j.1365-246X.2011.05015.x)
"""
function stretching(ref::AbstractArray, cur::AbstractArray, t::AbstractArray,
                    window::AbstractArray, freqbands::AbstractArray;
                    dvmin::Float64=-0.1, dvmax::Float64=0.1, ntrial::Int=100, norm::Bool=true)
     # define sampling interval
     dt = mean(diff(t))
     fs = 1/dt

     # calculate the CWT of the time series, using identical parameters for both calculations
     cwt1, sj, freqs, coi = cwt(ref, dt, minimum(freqbands), maximum(freqbands))
     cwt2, sj, freqs, coi = cwt(cur, dt, minimum(freqbands), maximum(freqbands))

     # if a frequency window is given (instead of a set of frequency bands), we assume
     # dv/v should be calculated for each frequency. We construct a 2D array of the
     # form [f1 f1; f2 f2; ...], which can be treated the same as a 2D array of frequency bands
     if ndims(freqbands)==1
         freqbands = hcat(freqs, freqs)
     end
     # number of frequency bands
     (nbands,_) = size(freqbands)

     # initialize dvv and err arrays
     dvv = zeros(nbands)
     cc = similar(dvv)
     cdp = similar(dvv)
     err = similar(dvv)
     allC = zeros(ntrial, nbands)

     # loop over frequency bands
     for iband=1:nbands
         (fmin, fmax) = freqbands[iband, :]

         # get current frequencies over which we apply icwt
         # frequency checks
         if fmax < fmin
             println("Error: please ensure columns 1 and 2 are right frequency limits in freqbands!")
         else
             freq_ind = findall(f->(f>=fmin && f<=fmax), freqs)
         end

         # perform icwt
         icwt1 = icwt(cwt1[:,freq_ind], sj[freq_ind], dt)
         icwt2 = icwt(cwt2[:,freq_ind], sj[freq_ind], dt)
         wcwt1 = real.(icwt1)
         wcwt2 = real.(icwt2)

         # normalize both signals, if appropriate
         if norm
             ncwt1 = ((wcwt1 .- mean(wcwt1)) ./ std(wcwt1))[:]
             ncwt2 = ((wcwt2 .- mean(wcwt2)) ./ std(wcwt2))[:]
         else
             ncwt1 = wcwt1[:]
             ncwt2 = wcwt2[:]
         end

         dvv[iband], cc[iband], cdp[iband], ϵ, err[iband], allC[:,iband] = stretching(ncwt1, ncwt2, t, window, fmin, fmax, dvmin=dvmin, dvmax=dvmax, ntrial=ntrial)
     end
     ϵ = collect(range(dvmin, stop=dvmax, length=ntrial))

     return freqbands, dvv, cc, cdp, ϵ, err, allC
end


text:
using Coverage

cd(joinpath(@__DIR__, "..", "..")) do
    processed = process_folder()
    covered_lines, total_lines = get_summary(processed)
    percentage = covered_lines / total_lines * 100
    println("($(percentage)%) covered")
end

# submit the report to codecov.io
Codecov.submit_local(Codecov.process_folder(), pwd(), slug = "LCSB-BioCore/GigaSOM.jl")


text:
@testset "Sessions functionality" begin

  @testset "Assets paths" begin
    using Genie, Genie.Sessions
    using Genie.Router
    using HTTP

    Sessions.init()

    route("/home") do
      sess = Sessions.session(Genie.Router.@params)
      Sessions.set!(sess, :visit_count, Sessions.get(sess, :visit_count, 0)+1)

      "$(Sessions.get(sess, :visit_count))"
    end

    Genie.up(; open_browser = false)

    # TODO: extend to use the cookie and increment the count
    response = HTTP.get("http://$(Genie.config.server_host):$(Genie.config.server_port)/home")
    @test response.body |> String == "1"

    Genie.down()
  end;

end;


text:
using RCall
R"install.packages('rstan', repos='http://cran.us.r-project.org')"


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

#=
- read meshes from different formats
- reorder connectivity, create element sets, node sets, ...
- create partitions for parallel runs
- renumber elements / nodes
- maybe precheck for bad elements
- check surface normal direction in boundary elements
- orientation of 2d elements
- etc only topology related stuff
=#

mutable struct Mesh
    nodes :: Dict{Int, Vector{Float64}}
    node_sets :: Dict{Symbol, Set{Int}}
    elements :: Dict{Int, Vector{Int}}
    element_types :: Dict{Int, Symbol}
    element_codes :: Dict{Int, Symbol}
    element_sets :: Dict{Symbol, Set{Int}}
    surface_sets :: Dict{Symbol, Vector{Tuple{Int, Symbol}}}
    surface_types :: Dict{Symbol, Symbol}
    coloring::Union{Nothing, Vector{Vector{Int}}} # Each vector contains a list of elements that do not share nodes
end

function Mesh()
    return Mesh(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), nothing)
end

"""
    Mesh(m::Dict)

Create a new `Mesh` using data `m`. It is assumed that `m` is in format what
`abaqus_read_mesh` in `AbaqusReader.jl` is returning.
"""
function Mesh(m::Dict)
    mesh = Mesh()
    mesh.nodes = m["nodes"]
    mesh.elements = m["elements"]
    mesh.element_types = m["element_types"]
    for (k, v) in m["surface_types"]
        mesh.surface_types[Symbol(k)] = v
    end
    for (nset_name, node_ids) in m["node_sets"]
        mesh.node_sets[Symbol(nset_name)] = Set(node_ids)
    end
    for (elset_name, element_ids) in m["element_sets"]
        mesh.element_sets[Symbol(elset_name)] = Set(element_ids)
    end
    for (surfset_name, surfaces) in m["surface_sets"]
        mesh.surface_sets[Symbol(surfset_name)] = surfaces
    end
    return mesh
end

"""
    add_node!(mesh, nid, ncoords)

Add node into the mesh. `nid` is node id and `ncoords` are the node
coordinates.
"""
function add_node!(mesh::Mesh, nid::Int, ncoords::Vector{Float64})
    mesh.nodes[nid] = ncoords
end

"""
    add_nodes!(mesh, nodes)

Add nodes into the mesh.
"""
function add_nodes!(mesh::Mesh, nodes::Dict{Int, Vector{Float64}})
    for (nid, ncoords) in nodes
        add_node!(mesh, nid, ncoords)
    end
end

"""
    add_node_to_node_set!(mesh, nid, ncoords)

Add nodes into a node set. `set_name` is the name of the set and `nids...`
are all the node id:s that wants to be added.
"""
function add_node_to_node_set!(mesh::Mesh, set_name, nids...)
    if !haskey(mesh.node_sets, set_name)
        mesh.node_sets[set_name] = Set{Int}()
    end
    push!(mesh.node_sets[set_name], nids...)
    return
end

"""
    create_node_set_from_element_set!(mesh, set_names...)

Create a new node set from the nodes in an element set. ´set_names...´ are all
the set names to be inserted in the function.
"""
function create_node_set_from_element_set!(mesh::Mesh, set_names::String...)
    for set_name in set_names
        set_name = Symbol(set_name)
        @info("Creating node set $set_name from element set")
        node_ids = Set{Int}()
        for elid in mesh.element_sets[set_name]
            push!(node_ids, mesh.elements[elid]...)
        end
        mesh.node_sets[set_name] = node_ids
    end
    return
end

"""
    create_node_set_from_element_set!(mesh, set_name)

Create a new node set from an element set.
"""
function create_node_set_from_element_set!(mesh::Mesh, set_name::Symbol)
    create_node_set_from_element_set!(mesh, string(set_name))
end

"""
    add_element!(mesh, elid, eltype, connectivity)

Add an element into the mesh. ´elid´ is the element id, ´eltype´ is the type of
the element and ´connectivity´ is the connectivity of the element.
"""
function FEMBase.add_element!(mesh::Mesh, elid, eltype, connectivity)
    mesh.elements[elid] = connectivity
    mesh.element_types[elid] = eltype
    return nothing
end

"""
    add_elements!(mesh, elements)

Add elements into the mesh.
"""
function FEMBase.add_elements!(mesh::Mesh, elements::Dict{Int, Tuple{Symbol, Vector{Int}}})
    for (elid, (eltype, elcon)) in elements
        add_element!(mesh, elid, eltype, elcon)
    end
    return nothing
end

"""
    add_element_to_element_set!(mesh, set_name, elids...)

Add elements into the mesh. ´set_name´ is the name of the element set and
´elids..´ are id:s of all the elements that wants to be added.
"""
function add_element_to_element_set!(mesh::Mesh, set_name, elids...)
    if !haskey(mesh.element_sets, set_name)
        mesh.element_sets[set_name] = Set{Int}()
    end
    push!(mesh.element_sets[set_name], elids...)
end

"""
    copy(mesh)

Return a copy of the mesh.
"""
function Base.copy(mesh::Mesh)
    mesh2 = Mesh()
    mesh2.nodes = copy(mesh.nodes)
    mesh2.node_sets = copy(mesh.node_sets)
    mesh2.elements = copy(mesh.elements)
    mesh2.element_types = copy(mesh.element_types)
    mesh2.element_sets = copy(mesh.element_sets)
    return mesh2
end

"""
    filter_by_element_id(mesh, element_ids)

Filter elements by their id's.
"""
function filter_by_element_id(mesh::Mesh, element_ids::Vector{Int})
    mesh2 = copy(mesh)
    mesh2.elements = Dict()
    for elid in element_ids
        if haskey(mesh.elements, elid)
            mesh2.elements[elid] = mesh.elements[elid]
        end
    end
    return mesh2
end

"""
    filter_by_element_set(mesh, set_name)

Filter elements by an element set.
"""
function filter_by_element_set(mesh::Mesh, set_name)
    filter_by_element_id(mesh::Mesh, collect(mesh.element_sets[set_name]))
end

"""
    create_element(mesh, id)

Create an element from the mesh by it's id.
"""
function create_element(mesh::Mesh, id::Int)
    connectivity = mesh.elements[id]
    element_type = getfield(JuliaFEM, mesh.element_types[id])
    element = Element(element_type, connectivity)
    element.id = id
    update!(element, "geometry", mesh.nodes)
    return element
end

function create_elements(mesh::Mesh; element_type=nothing)
    element_ids = collect(keys(mesh.elements))
    if element_type != nothing
        filter!(id -> mesh.element_types[id] == element_type, element_ids)
    end
    elements = [create_element(mesh, id) for id in element_ids]
    return elements
end

function create_elements(mesh::Mesh, element_sets::Symbol...; element_type=nothing)
    if isempty(element_sets)
        element_ids = collect(keys(mesh.elements))
    else
        element_ids = Set{Int}()
        for set_name in element_sets
            element_ids = union(element_ids, mesh.element_sets[set_name])
        end
    end

    if element_type != nothing
        filter!(id -> mesh.element_types[id] == element_type, element_ids)
    end

    elements = [create_element(mesh, id) for id in element_ids]

    nelements = length(elements)
    content = Dict{Symbol, Int}()
    for elid in element_ids
        eltype = mesh.element_types[elid]
        content[eltype] = get(content, eltype, 0) + 1
    end
    s = join(("$v x $k" for (k, v) in content), ", ")
    v = join(element_sets, ", ")
    @info("Created $nelements elements ($s) from element set: $v.")

    return elements
end

"""
    create_elements(mesh::Mesh, element_set::String)

# Examples

Suppose that there is a `mesh` with element set `Body_1`. Creating elements
based on that element set is done then

```julia
create_elements(mesh, "Body_1")
```
"""
function create_elements(mesh::Mesh, element_sets::String...)
    return create_elements(mesh, map(Symbol, element_sets)...)
end

"""
    find_nearest_nodes(mesh, coords, npts=1; node_set=nothing)

find npts nearest nodes from the mesh and return their id numbers as a list.
"""
function find_nearest_nodes(mesh::Mesh, coords::Vector{Float64}, npts::Int=1; node_set=nothing)
    dist = Dict{Int, Float64}()
    for (nid, c) in mesh.nodes
        if node_set != nothing && !(nid in mesh.node_sets[Symbol(node_set)])
            continue
        end
        dist[nid] = norm(coords-c)
    end
    s = sort(collect(dist), by=x->x[2])
    nd = s[1:npts] # [(id1, dist1), (id2, dist2), ..., (id_npts, dist_npts)]
    node_ids = [n[1] for n in nd]
    return node_ids
end

function find_nearest_node(mesh::Mesh, coords::Vector{Float64}; node_set=nothing)
    return first(find_nearest_nodes(mesh, coords, 1; node_set=node_set))
end

"""
    reorder_element_connectivity!(mesh, mapping)

Apply a new node ordering to elements. JuliaFEM uses the same node ordering as
ABAQUS. If the mesh is parsed from FEM format with some other node ordering,
this function can be used to reorder the nodes.

Parameters
----------
mapping :: Dict{Symbol, Vector{Int}}
    e.g. :Tet10, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

"""
function reorder_element_connectivity!(mesh::Mesh, mapping::Dict{Symbol, Vector{Int}})
    for (elid, eltype) in mesh.element_types
        haskey(mapping, eltype) || continue
        new_order = mapping[eltype]
        element_connectivity = mesh.elements[elid]
        new_element_connectivity = element_connectivity[new_order]
        mesh.elements[elid] = new_element_connectivity
    end
end

function JuliaFEM.Problem(mesh::Mesh, ::Type{P}, name::AbstractString, dimension::Int) where P<:FieldProblem
    problem = Problem(P, name, dimension)
    problem.elements = create_elements(mesh, name)
    return problem
end

function JuliaFEM.Problem(mesh::Mesh, ::Type{P}, name, dimension, parent_field_name) where P<:BoundaryProblem
    problem = Problem(P, name, dimension, parent_field_name)
    problem.elements = create_elements(mesh, name)
    return problem
end

"""
    create_coloring!(mesh::Mesh)

Greedy algorithm for coloring a grid such that no two cells with the same node
have the same color.
This function sets the `coloring` field in `mesh` to a `Vector{Vector{Int}}` where
each vector contains vectors of elements that do not share any nodes.
It is therefore safe to assemble in parallel each element vector by vector.
"""
function create_coloring!(mesh::Mesh)
    # Contains the elements that each node contain
    cell_containing_node = Dict{Int, Set{Int}}()
    for (cellid, nodes) in mesh.elements
        for v in nodes
            if !haskey(cell_containing_node, v)
                cell_containing_node[v] = Set{Int}()
            end
            push!(cell_containing_node[v], cellid)
        end
    end

    I, J, V = Int[], Int[], Bool[]
    for (node, cells) in cell_containing_node
        for cell1 in cells # All these cells have a neighboring node
            for cell2 in cells
                if cell1 != cell2
                    push!(I, cell1)
                    push!(J, cell2)
                    push!(V, true)
                end
            end
        end
    end

    incidence_matrix = sparse(I, J, V)
    # cell -> color of cell
    cell_colors = Dict{Int, Int}()
    # color -> list of cells
    final_colors = Vector{Int}[]
    occupied_colors = Set{Int}()
    # Zero represents no color set yet
    for (cellid, _) in mesh.elements
        cell_colors[cellid] = 0
    end
    total_colors = 0
    for (cellid, _) in mesh.elements
        empty!(occupied_colors)
        # loop over neighbors
        for r in nzrange(incidence_matrix, cellid)
            cell_neighbour = incidence_matrix.rowval[r]
            color = cell_colors[cell_neighbour]
            if color != 0
                push!(occupied_colors, color)
            end
        end

        # occupied colors now contains all the colors we are not allowed to use
        free_color = 0
        for attempt_color in 1:total_colors
            if attempt_color ∉ occupied_colors
                free_color = attempt_color
                break
            end
        end

        if free_color == 0 # no free color found, need to bump max colors
            total_colors += 1
            free_color = total_colors
            push!(final_colors, Int[])
        end

        cell_colors[cellid] = free_color
        push!(final_colors[free_color], cellid)
    end

    mesh.coloring = final_colors
    return mesh
end


text:
function spline_basis(domain, p)
    R_spline = R"splines::ns($domain, $p)"
    Q = convert(Matrix{Float64}, R_spline)
    return Q
end

struct Atomic <: Distributions.ValueSupport end
struct ExpoSpline{T<:Real} <: UnivariateDistribution{Atomic}
    domain::Vector{Float64}
    Qbasis::Matrix{Float64} # basis matrix
    αcoef::Vector{T} # α
    Qα::Vector{T}
    φ_α::T # normalization constant
end
function ExpoSpline(domain, Qbasis::Matrix{Float64}, αcoef::Vector{T}) where {T <: Real}
    Qα = Qbasis * αcoef
    φ_α = log_sum_exp(Qα) # equation (10)
    return ExpoSpline(domain, Qbasis, αcoef, Qα, φ_α)
end
function ExpoSpline(domain, αcoef::Vector{T}) where {T <: Real}
    p = length(αcoef)
    Q = spline_basis(domain, p)
    return ExpoSpline(domain, Q, αcoef)
end
## doesn't work for ExpoSpline immutable ##
# function update_coefs!(d::ExpoSpline{T}, α::Vector{T}) where {T <: Real}
    # d.αcoef = α
    # LinAlg.A_mul_B!(d.Qα, d.Qbasis, α)
    # d.φ_α = log_sum_exp(d.Qα)
    # return d
# end

####################################
# Expo Spline distribution methods #
####################################

params(d::ExpoSpline) = (d.αcoef, )
coefs(d::ExpoSpline) = d.αcoef
partype(d::ExpoSpline{T}) where T = T

function log_sum_exp(x)
    # https://en.wikipedia.org/wiki/LogSumExp
    max_x = maximum(x)
    log_sum_exp_x = max_x + log(sum(xi -> exp(xi-max_x), x))
    return log_sum_exp_x
end
logpdf(d::ExpoSpline) = d.Qα .- d.φ_α             # equation (9)
pdf(d::ExpoSpline) = exp.(logpdf(d))
logpdf(d::ExpoSpline, j::Int) = d.Qα[j] - d.φ_α # equation (9)
function logpdf(d::ExpoSpline, x::Real)
    j = argmin(abs.(d.domain .- x)) # not the fastest way to do this
    return logpdf(d, j)
end
function pdf(d::ExpoSpline, x::Real)
    return exp(logpdf(d, x))
end

function logcdf(d::ExpoSpline)
    max_x = maximum(d.Qα)
    cumsum_exp_x = accumulate((s,xi) -> s + exp(xi - max_x), d.Qα; init=zero(max_x))
    log_last = log(cumsum_exp_x[end])
    lcdf = cumsum_exp_x # to modify in place
    @inbounds for (i, x) in enumerate(cumsum_exp_x)
        lcdf[i] = log(x) - log_last
    end
    return lcdf
end
cdf(d::ExpoSpline) = exp.(logcdf(d))
cdf(d::ExpoSpline, j::Int) = exp(logcdf(d)[j]) # wasteful?
function cdf(d::ExpoSpline, x::Real)
    j = argmin(abs.(d.domain .- x)) # inefficient
    return cdf(d, j)
end

function rand(d::ExpoSpline)
    lcdf = logcdf(d)
    logu = -Random.randexp()
    i = searchsortedfirst(lcdf, logu)
    return d.domain[min(i,length(lcdf))]
end
function rand(d::ExpoSpline, n::Int)
    lcdf = logcdf(d)
    logu = -Random.randexp(n)
    r = similar(d.domain, n)
    for j in 1:n
        logu = -Random.randexp()
        i = searchsortedfirst(lcdf, logu)
        r[j] = d.domain[min(i,length(lcdf))]
    end
    return r
end

#################
# Likelihood ####
#################
function get_logPi(domain::Vector{Float64}, noise::UnivariateDistribution, X_i::Float64)
    # unnormalized
    return logpdf.(noise, X_i.-domain)
end
function get_Pi(domain::Vector{Float64}, noise::UnivariateDistribution, X_i::Float64)
    # unnormalized
    return pdf.(noise, X_i.-domain)
end
get_logPi(d::ExpoSpline, noise, X_i) = get_logPi(d.domain, noise, X_i)
get_Pi(d::ExpoSpline, noise, X_i) = get_Pi(d.domain, noise, X_i)
function loglikelihood(logPi::AbstractVector{Float64}, loggα::Vector{T}) where T<:Real
    max_x = maximum(logPi) + maximum(loggα)
    sum_exp_x = zero(max_x)
    @inbounds for (logPi_j, loggα_j) in zip(logPi, loggα)
        sum_exp_x += exp(logPi_j + loggα_j - max_x)
    end
    return max_x + log(sum_exp_x)
end
function loglikelihood(d::ExpoSpline, logPi::AbstractVector{Float64})
    loggα = logpdf(d)
    return loglikelihood(logPi, loggα)
end
function getWi(d::ExpoSpline, logPi::AbstractVector{Float64})
    # equation (14)-(15)
    loggα = logpdf(d)
    logfi = loglikelihood(logPi, loggα)
    Wi = exp.(loggα) .* (exp.(logPi .- logfi) .- 1)
    return Wi
end
function dloglik(d::ExpoSpline, logPi::AbstractVector{Float64})
    # equation (16)
    return d.Qbasis'getWi(d, logPi)
end
function d2loglik(d::ExpoSpline, logPi::AbstractVector{Float64})
    gα = pdf(d)
    Wi = getWi(d, logPi)
    return -d.Qbasis'*(Wi*Wi' + Wi*gα' + gα*Wi' - diagm(Wi))*d.Qbasis
end
function loglikelihood(d::ExpoSpline, noise::UnivariateDistribution, X_i::Float64)
    logPi = get_logPi(d.domain, noise, X_i)
    return loglikelihood(d, logPi)
end

function make_loglik_w_noise_pre(d::ExpoSpline, noise, X::Vector{Float64})
    logP = get_logPi.(d, noise, X)
    mloglik_w_noise_pre = function (α)
        dα = ExpoSpline(d.domain, d.Qbasis, α)
        loglik = sum(Pi -> loglikelihood(dα, Pi), logP)
        return -loglik
    end
    return mloglik_w_noise_pre
end
function make_loglik_noisefree(d::ExpoSpline, θ::Vector{Float64})
    mloglik = (α) -> -sum(
        logpdf.( # compute the log-PDF under the…
            ExpoSpline(d.domain, d.Qbasis, α), # …candidate distribution…
            θ # …of the noise-free data
            )
        )
    return mloglik
end
##############
# with Prior #
##############
function make_penalized_w_noise_pre(d::ExpoSpline, noise, X::Vector{Float64}, c0::Real)
    logP = get_logPi.(d, noise, X)
    target = function (α)
        dα = ExpoSpline(d.domain, d.Qbasis, α)
        loglik = sum(Pi -> loglikelihood(dα, Pi), logP)
        penalty = c0 * norm(α, 2)
        return -loglik + penalty
    end
    return target
end

################
# Optimization #
################
function decon(d, noise, X::Vector{Float64}, c0::Real)
    f = make_penalized_w_noise_pre(d, noise, X, c0)
    td = Optim.TwiceDifferentiable(f, copy(coefs(d)); autodiff=:forward)
    opt_reg = Optim.optimize(td, zero(coefs(d)), Optim.Newton())
    best_x = Optim.minimizer(opt_reg)
    best_y = Optim.minimum(opt_reg)
    for _ in 1:5
        if Optim.converged(opt_reg)
            break
        end
        init_x = randn(length(coefs(d)))
        opt_reg = Optim.optimize(td, init_x, Optim.Newton())
        if Optim.minimum(opt_reg) < best_y
            best_x = Optim.minimizer(opt_reg)
            best_y = Optim.minimum(opt_reg)
        end
    end
    for _ in 1:5
        if Optim.converged(opt_reg)
            break
        end
        init_x = randn(length(coefs(d)))
        opt_reg = Optim.optimize(td, init_x, Optim.ConjugateGradient())
        if Optim.minimum(opt_reg) < best_y
            best_x = Optim.minimizer(opt_reg)
            best_y = Optim.minimum(opt_reg)
        end
    end
    if !Optim.converged(opt_reg)
        @warn("decon failed to converge")
    end
    d_optim = ExpoSpline(d.domain, d.Qbasis, best_x)
    return d_optim
end


text:
# plain try throws
@test_throws LoadError macroexpand(@__MODULE__, :(@🐒 begin
🙈
    println("Should not have come this far.")
    42
end))
# try...catch
@test 42 == @🐒 begin
🙈
    2 + 40
🙊(e)
    rethrow()
end

@test 12 == @🐒 begin
🙈
    2 + "40"
    error("This was not supposed to happen.")
🙊(e)
    e isa MethodError || rethrow()
    12
end

@test_throws DomainError @🐒 begin
🙈
    c╯°□°ↄ╯(DomainError("On purpose!"))
🙊(e)
    rethrow()
end
# try...finally
@test @test_logs (:warn, "This worked!") @🐒(begin
🙈
    3*3
🙉
    @warn "This worked!"
end) == 9 # note that the parenthesis around (begin...end) are necessary here because == is infix

@test_logs (:info, "I ran!") @test_throws ErrorException @🐒 begin
🙈
    error("On purpose!")
🙉
    @info "I ran!"
end
# try...catch...finally
@test_logs (:info, "I'm surprised that this works!") @test_throws ArgumentError @🐒 begin
🙈
    c╯°□°ↄ╯(ArgumentError("Skyrim belongs to the Nords!"))
🙊(e)
    rethrow()
🙉
    @info "I'm surprised that this works!"
end

@test_logs (:warn, "Finally!") @test "Done" == @🐒 begin
🙈
    "D"*"o"*"n"*"e"
🙊(e)
    rethrow()
🙉
    @warn "Finally!"
end
# nested monkeys
@test_logs (:info, "Try Finally") (:warn, "Catch Finally") @test_throws ErrorException @🐒 begin
🙈
    @🐒 begin
    🙈
        throw(ArgumentError("Inner throw"))
    🙊(e)
        e isa ArgumentError && throw(DomainError("Got caught once."))
    🙉
        @info "Try Finally"
    end
🙊(e)
    @🐒 begin
    🙈
        e isa DomainError && throw(InitError(:here, e))
    🙊(e)
        e isa InitError && error("I'm still alive")
    🙉
        @warn "Catch Finally"
    end
end
# I'm too lazy to test all the other combinations

# leaving out the 🙈 throws
@test_throws LoadError macroexpand(@__MODULE__, :(@🐒 begin
🙊(e)
    "Lol."
🙉
    "No!"
end))


text:
function runcode(app, numpde)

display("run code...");

if numpde==1
    mystr = app.appname * " 1 datain/ dataout/out";
else
    mystr = app.appname * " " * string(numpde);
    for i = 1:numpde
        mystr = mystr * " datain" * string(i) * "/ " * "dataout" * string(i) * "/out";
    end
end

mpirun = app.mpirun;
if app.platform == "cpu"
    if app.mpiprocs==1
        runstr = "./app/serial" * mystr;
    else
        runstr = mpirun * " -np " * string(app.mpiprocs) * " ./app/mpi" * mystr;
    end
    run(string2cmd(runstr));
elseif app.platform == "gpu"
    if app.mpiprocs==1
        runstr = "./app/gpu" * mystr;
    else
        runstr = mpirun * " -np " * string(app.mpiprocs) * " ./app/gpumpi" * mystr;
    end
    run(string2cmd(runstr));
else
    error("app.platform must be either cpu or gpu");
end

return runstr

end

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim


text:
for F in (:(+), :(-), :(*), :(/), :(\), :(^))
  @eval begin
    function $F(a::MarkFloat64, b::MarkFloat64)
        c = $F(Float64(a), Float64(b))
        d = reinterpret(UInt64, c)
        d |= marker(a)
        d |= marker(b)
        z = reinterpret(MarkFloat64, d)
        return z
    end
  end
end

for F in (:(-), :abs, :log, :log2, :log10, :log1p, :exp, :expm1,
          :sin, :cos, :tan, :sec, :csc, :cot, :sinpi, :cospi,
          :asin, :acos, :atan, :asec, :acsc, :acot,
          :sinh, :cosh, :tanh, :sech, :csch, :coth,
          :asinh, :acosh, :atanh, :asech, :acsch, :acoth)
  @eval begin
    function $F(a::MarkFloat64)
        c = Float64(a)
        c = $F(c)
        d = UInt64(c) | marker(a)
        z = MarkFloat64(d)
        return z
    end
  end
end
# log(reinterpret(Float64,MarkableFloats.unmark(c)))


text:
# Pure Julia decimal arithmetic
# @license MIT
# @author jack@tinybike.net (Jack Peterson), 7/3/2014

module Decimals
    import Base: ==, +, -, *, /, <, <=, float, inv, round, trunc

    export Decimal,
           decimal,
           number,
           normalize

    const DIGITS = 20

    # Numerical value: (-1)^s * c * 10^q
    struct Decimal <: AbstractFloat
        s::Integer  # sign can be 0 (+) or 1 (-)
        c::BigInt   # coefficient (significand), must be non-negative
        q::Integer  # exponent
    end

    # Convert between Decimal objects, numbers, and strings
    include("decimal.jl")

    # Decimal normalization
    include("norm.jl")

    # Addition, subtraction, negation, multiplication
    include("arithmetic.jl")

    # Equality
    include("equals.jl")

    # Rounding
    include("round.jl")

end


text:
using Test
using WebIO

@testset "JSString Interpolations" begin
    @testset "Interpolations into JSStrings" begin
        text = "Hello, world!"
        js_text = js"console.log($text);"
        @test js_text.s == """console.log("Hello, world!");"""

        dict = Dict("foo" => "bar")
        js_dict = js"console.log($dict);"
        @test js_dict.s == "console.log({\"foo\":\"bar\"});"
    end

    @testset "Interpolations of JSStrings into (normal) Strings" begin
        js_log = js"""console.log("Hello!");"""
        script = "<script>$js_log</script>"
        @test script == """<script>console.log("Hello!");</script>"""
    end

    @testset "@js_str Escapes Correctly" begin
        @test js"this.\$refs".s == "this.\$refs"
        @test js"foo\$".s == "foo\$"
        @test js"\$('div.my-id')".s == "\$('div.my-id')"

        @test js"'foo\\bar'".s == "'foo\\bar'"

        foo = "foo"
        @test js"\\$foo".s == "\\\"foo\""

        # See note about Julia and weirdness with escaping when quotes are involved.
        @test js"""console.log("\\")""".s == """console.log(\"\\\")"""
        @test js"""console.log(\"\\\")""".s == """console.log(\"\\\")"""
        @test js"""foo = '\\\\'""".s == raw"""foo = '\\'"""
    end

    @testset "@js_str interpolates JSStrings correctly" begin
        myfunc = js"alert";
        @test js"$myfunc(123)".s == "alert(123)"

        myvalue = js"a + 1"
        @test js"x = $myvalue".s == "x = a + 1"
    end
end


text:
#
# constants.jl -
#
# Constants defined for Nüvü Camēras.
#
# *DO NOT EDIT* Some parts of this file are generated automatically.
#
const OPENGL_ENABLE  = Cint(0)

const NC_AUTO_DETECT    = Cint(0x0000ffff)
const NC_AUTO_CHANNEL   = NC_AUTO_DETECT
const NC_AUTO_UNIT	= Cint(0x6fffffff)
const NC_FULL_WIDTH     = Cint(-1)
const NC_FULL_HEIGHT    = Cint(-1)
const NC_USE_MAC_ADRESS = Cint(0x20000000)

# Type of amplifier set by the seleted readout mode.
@enum(Ampli::Cint, NOTHING = 0, EM, CONV)

# Type of communication interface that can be used.
@enum(CommType::Cint,
      INVALID     = 0x0000fffe,
      AUTO_COMM   = 0x6fff0000, # NC_AUTO_UNIT - NC_AUTO_DETECT
      VIRTUAL     = 0x7fff0000,
      EDT_CL      = 0x00000000,
      PT_GIGE     = 0x00010000,
      MIL_CL      = 0x00020000,
      SAPERA_CL   = 0x00030000,
      UNKNOWN     = 0x00040000,
      REMOTE_COMM = 0x00050000)

@enum(DetectorType::Cint,
      CCD60 = 0,
      CCD97,
      CCD201_20,
      CCD207_00,
      CCD207_10,
      CCD207_40,
      CCD220,
      UNKNOWN_CCD)

# All camera parameters.
@enum(Features::Cint,
      UNKNOWN_FEATURE,
      EXPOSURE,
      WAITING_TIME,
      TRIGGER_MODE,
      SHUTTER_MODE,
      SHUTTER_POLARITY,
      EXTERNAL_SHUTTER,
      EXTERNAL_SHUTTER_MODE,
      EXTERNAL_SHUTTER_DELAY,
      FIRE_POLARITY,
      MINIMUM_PULSE_WIDTH,
      ARM_POLARITY,
      CALIBRATED_EM_GAIN,
      RAW_EM_GAIN,
      ANALOG_GAIN,
      ANALOG_OFFSET,
      TARGET_DETECTOR_TEMP,
      DETECTOR_TEMP,
      CTRL_TIMESTAMP,
      GPS_PRESENT,
      GPS_LOCKED,
      BINNING_X,
      BINNING_Y,
      ROI,
      CROP_MODE_X,
      CROP_MODE_Y,
      ACTIVE_REGION,
      MULTIPLE_ROI,
      TRIGGER_MODE_INTERNAL,   # INTERNAL trigger mode
      TRIGGER_MODE_EXTERNAL,   # EXT_HIGH_LOW or EXT_LOW_HIGH
      TRIGGER_MODE_EXPOSURE,   # EXT_HIGH_LOW_EXP or EXT_HIGH_LOW_EXP
      TRIGGER_MODE_CONTINUOUS, # CONT_HIGH_LOW or CONT_LOW_HIGH
      SYSTEM_STATUS,
      SHORT_SSVA_CMD,
      BIN_CDS)

# For backwards compatibiliy
const Param = Features

# List of different temperatures
@enum(NcTemperatureType::Cint,
      NC_TEMP_CCD,
      NC_TEMP_CONTROLLER,
      NC_TEMP_POWER_SUPPLY,
      NC_TEMP_FPGA,
      NC_TEMP_HEATINK)

# Reasons why a tap was reported as unused
@enum(NcPortUnusedReason::Cint,
      NC_PORT_UNUSED_NO_DEVICE,      # The port was found to not be connected to a device
      NC_PORT_UNUSED_ALREADY_IN_USE, # The device is already in use by another instance
      NC_PORT_UNUSED_FILTERED,       # Port was not scanned (filtered, internal use only)
      NC_PORT_UNUSED_UNKNOWN)        # The reason was unusual and has no unique identifier

# Type of shutter mode available
@enum(ShutterMode::Cint, SHUT_NOT_SET = 0, OPEN = 1, CLOSE = 2, AUTO = 3)
const BIAS_DEFAULT = CLOSE

# Type of trigger mode available
@enum(TriggerMode::Cint,
      CONT_HIGH_LOW = -3,
      EXT_HIGH_LOW_EXP = -2,
      EXT_HIGH_LOW = -1,
      INTERNAL = 0,
      EXT_LOW_HIGH = 1,
      EXT_LOW_HIGH_EXP = 2,
      CONT_LOW_HIGH = 3)

# Type of crop-mode available
@enum(CropMode::Cint,
      CROP_MODE_DISABLE = 0,
      CROP_MODE_ENABLE_X = 1,
      CROP_MODE_ENABLE_Y = 2,
      CROP_MODE_ENABLE_XY = 3)
const CROP_MODE_ENABLE_ZL = CROP_MODE_ENABLE_X
const CROP_MODE_ENABLE_SP = CROP_MODE_ENABLE_Y
const CROP_MODE_ENABLE_MX = CROP_MODE_ENABLE_XY

@enum(ExtShutter::Cint, EXTERNAL_SHUTTER_DISABLE = 0, EXTERNAL_SHUTTER_ENABLE = 1)

@enum(ExtPolarity::Cint,NEGATIVE = -1, POL_NOT_SET = 0, POSITIVE = 1)

# Type of image format which will be used to save an image
@enum(ImageFormat::Cint, UNKNOWNN = -1, TIF = 0, FITS)

# Use only CFITSIO DATATYPE codes
@enum(ImageDataType::Cint,
      NC_IMG,       # uint 16 bit
      NC_IMG_UINT,  # uint 32 bit
      NC_IMG_ULONG, # uint 64 bit
      NC_IMG_FLOAT) # 32 bit float

@enum(HeaderDataType::Cint, NC_INT, NC_DOUBLE, NC_STRING)

@enum(ImageCompression::Cint, NO_COMPRESSION = 0, GZIP)

@enum(ProcType::Cint, NO_PROC = 0x00, LM = 0x01, PC = 0x02)

@enum(TimestampMode::Cint,
      NO_TIMESTAMP = 0,
      INTERNAL_TIMESTAMP = 1,
      GPS_TIMESTAMP = 2)

# List of possible version numbers.
@enum(VersionType::Cint, SOFTWARE, FIRMWARE, FPGA, HARDWARE)

# Tiff tags.
@enum(TiffTag::Cint,
      AMPLI_TYPE_TAG           = 650,
      HORIZ_FREQ_TAG           = 651,
      VERT_FREQ_TAG            = 652,
      EXPOSURE_TIME_TAG        = 653,
      EFFECTIVE_EXP_TIME_TAG   = 682,
      WAITING_TIME_TAG         = 654,
      RAW_EM_GAIN_TAG          = 655,
      CAL_EM_GAIN_TAG          = 680,
      ANALOG_GAIN_TAG          = 656,
      ANALOG_OFFSET_TAG        = 657,
      TARGET_DETECTOR_TEMP_TAG = 658,
      BINNING_X_TAG            = 659,
      BINNING_Y_TAG            = 660,
      ROI_X1_TAG               = 661,
      ROI_X2_TAG               = 662,
      ROI_Y1_TAG               = 663,
      ROI_Y2_TAG               = 664,
      CROP_MODE_ENABLE_TAG     = 695,
      CROP_MODE_PADDING_X_TAG  = 696,
      CROP_MODE_PADDING_Y_TAG  = 697,
      SHUTTER_MODE_TAG         = 665,
      TRIGGER_MODE_TAG         = 666,
      CLAMP_LEVEL_TAG          = 667,
      PROC_TYPE_TAG            = 676,
      NBR_PC_IMAGES_TAG        = 675,
      SERIAL_NUMBER_TAG        = 668,
      FIRM_VERSION_TAG         = 669,
      HARDW_VERSION_TAG        = 670,
      FPGA_VERSION_TAG         = 671,
      SOFT_VERSION_TAG         = 672,
      ADD_COMMENTS_TAG         = 674,
      OVERSCAN_LINES_TAG       = 677,
      OVERSAMPLE_X_TAG         = 678,
      DETECTOR_TYPE_TAG        = 681,
      FIRE_POLARITY_TAG        = 688,
      SHUTTER_POLARITY_TAG     = 683,
      ARM_POLARITY_TAG         = 689,
      OUTPUT_PULSE_WIDTH_TAG   = 687,
      EXT_SHUTTER_PRES_TAG     = 684,
      EXT_SHUTTER_MODE_TAG     = 685,
      EXT_SHUTTER_DELAY_TAG    = 686,
      OBS_DATE_TAG             = 690,
      OBS_TIME_SEC_FRACT_TAG   = 691,
      OBS_FLAGS_TAG            = 692,
      DATE_TIME_US_TAG         = 694,
      HOST_TIME_TAG            = 693,
      NBR_IMAGES_TRIG_TAG      = 679,
      SEQUENCE_NAME_TAG        = 673)

# Common data type enumeration
@enum(NcSdkDataTypes::Cint,
      NcSdkDataTypeInt = 0,
      NcSdkDataTypeDouble,
      NcSdkDataTypeString)


text:
"""
    norm1(v)

Normalize `v` to sum 1.
"""
function norm1(v)
    vs = sum(v)
    if vs ≈ 1.0
        return v
    else
        return v ./ vs
    end
end

"""
    infovalue(p, q)

Compute the symmetric Relative Entropy (Kullback-Liebler Divergence) between `p` and `q`,
two probability vectors or two frequency count vectors.
Value is bounded within [0, ∞].
"""
function infovalue(p, q)
    pn = norm1(p)
    qn = norm1(q)

    kldivergence(pn, qn) + kldivergence(qn, pn)
end
infovalue(f::AbstractMatrix) = infovalue(f[:, 1], f[:, 2])
infovalue(f::NamedArray) = infovalue(f.array)

"""
    cramerv(f::AbstractMatrix{T}) where T <: Real

Compute Cramer's V of a contingency table `f`.\\
ϕ coefficient is √(χ²/n), bounded within [0, √min(row-1,col-1)],
Cramer's V is ϕ/√min(row-1,col-1), bounded within [0,1].
"""
function cramerv(f::AbstractMatrix{T}) where {T<:Real}
    minimum(size(f)) >= 2 || throw(ArgumentError("Matrix needed, not single row or column"))

    p = norm1(f)
    rp = sum(p; dims = 2)               # row marginal
    cp = sum(p; dims = 1)               # column marginal
    ep = rp .* cp                       # expected probabilities
    ϕ² = sum((p .- ep) .^ 2 ./ ep)      # ϕ²

    sqrt(ϕ² / (minimum(size(f)) - 1))   # Cramer's V
end
cramerv(f::NamedArray) = cramerv(f.array)

"""
    cramerv(x::Vector{T} where T<:Integer, y::Vector{T} where T<:Integer)

Compute Cramer's V of two presumed discrete vectors `x` and `y`.\\
"""
cramerv(x::Vector{T} where {T<:Integer}, y::Vector{T} where {T<:Integer}) =
    cramerv(freqtable(x, y).array)

"""
    cramerv(x::Vector, y::Vector; groups = 10)

Compute Cramer's V of two presumed continuous vectors `x` and `y`,
vectors are binned into `groups` first.
"""
function cramerv(x::AbstractVector, y::AbstractVector; groups = 10)
    xr = ranks(x; groups = groups)
    yr = ranks(y, groups = groups)

    cramerv(freqtable(xr, yr).array)
end

"""
    mutualinfo(f::AbstractMatrix{T}) where T <: Real

Compute the mutual information of frequency matrix `f`;
value is bounded within [0, log(min(size(`f`)))].
"""
function mutualinfo(f::AbstractMatrix{T} where {T<:Real})
    minimum(size(f)) >= 2 || throw(ArgumentError("Matrix needed, not single row or column"))

    ps = norm1(f)
    px = sum(ps; dims = 2)
    py = sum(ps; dims = 1)

    entropy(px) + entropy(py) - entropy(ps)
end
mutualinfo(f::NamedArray) = mutualinfo(f.array)

"""
    mutualinfo(x::Vector{T} where T<:Integer, y::Vector{T} where T<:Integer)

Compute mutual information of two presumed discrete vectors `x` and `y`.\\
"""
mutualinfo(x::Vector{T} where {T<:Integer}, y::Vector{T} where {T<:Integer}) =
    mutualinfo(freqtable(x, y).array)

"""
    mutualinfo(x::Vector, y::Vector)

Compute the mutual information of two presumed continuous vectors `x` and `y`.
"""
function mutualinfo(x::AbstractVector, y::AbstractVector; groups = 10)
    xr = ranks(x; groups = groups)
    yr = ranks(y, groups = groups)

    mutualinfo(freqtable(xr, yr).array)
end

"""
    eda(df::AbstractDataFrame, target::Symbol)

Return a dataframe of Mutual Information and Cramer's V between `target` and other
variables in `df`. If `target` is binary, Information Value is also returned.

The dataframe is sorted by descending Mutual Information.
"""
function eda(df::AbstractDataFrame, target::Symbol; groups = 20)::AbstractDataFrame
    t = df[!, target]
    tnlvl = length(unique(t))
    tnlvl <= 1 && throw(ArgumentError("Target is single valued"))

    if tnlvl > groups
        t = ranks(t, groups = groups)
        tnlvl = length(unique(t))
    end

    println("Target: $target   type: $(typeof(t))   Levels: $tnlvl")

    lk = ReentrantLock()

    if tnlvl == 2
        out = DataFrame(
            Variable   = Symbol[],
            Vartype    = String[],
            Varlvls    = Int[],
            MutualInfo = Float64[],
            CramerV    = Float64[],
            InfoValue  = Float64[],
        )

        Threads.@threads for v in propertynames(df)
            v == target && continue

            vtype = eltype(df[!, v])
            if vtype <: Union{Missing,Real}
                vb = ranks(df[!, v], groups = groups)
            else
                vb = df[!, v]
            end

            frq = proptable(vb, t).array
            vnlvl = size(frq, 1)
            if vnlvl == 1
                println("Warning: [$v] is singled valued, skipped.")
                continue
            elseif vnlvl > 50
                println("Note: [$v] has more than 100 levels, consider binning")
            end

            mutin = mutualinfo(frq)
            cramv = cramerv(frq)
            iv = infovalue(frq[:, 1], frq[:, 2])

            lock(lk)
            push!(out, (v, string(vtype), vnlvl, mutin, cramv, iv))
            unlock(lk)
        end
    else
        out = DataFrame(
            Variable   = Symbol[],
            Vartype    = String[],
            Varlvls    = Int[],
            MutualInfo = Float64[],
            CramerV    = Float64[],
        )

        Threads.@threads for v in propertynames(df)
            v == target && continue

            vtype = eltype(df[!, v])
            if vtype <: Union{Missing,Real}
                vb = ranks(df[!, v], groups = groups)
            else
                vb = df[!, v]
            end

            frq = proptable(vb, t).array
            vnlvl = size(frq, 1)
            if vnlvl == 1
                println("Warning: [$v] is singled valued, skipped.")
                continue
            elseif vnlvl > 50
                println("Note: [$v] has more than 100 levels, consider binning")
            end

            mutin = mutualinfo(frq)
            cramv = cramerv(frq)

            lock(lk)
            push!(out, (v, string(vtype), vnlvl, mutin, cramv))
            unlock(lk)
        end
    end

    sort!(out, [:MutualInfo, :CramerV], rev = true)
end


text:
module Tokens

import Base.eof
export Token
include("token_kinds.jl")

iskeyword(k::Kind) = begin_keywords < k < end_keywords
#isliteral(k::Kind) = begin_literal < k < end_literal
isoperator(k::Kind) = begin_ops < k < end_ops

# TODO: more
@enum(TokenError,
  NO_ERR,
  UNKNOWN,
)

struct Token
    kind::Kind
    # Offsets into a string or buffer
    startpos::Tuple{Int, Int} # row, col where token starts /end, col is a string index
    endpos::Tuple{Int, Int}
    startbyte::Int # The byte where the token start in the buffer
    endbyte::Int # The byte where the token ended in the buffer
    val::Union{String,Int64,Float64} #Compat.String # The actual string of the token
    token_error::TokenError
end

function kind(t::Token)
    isoperator(t.kind) && return OP
    iskeyword(t.kind) && return KEYWORD
    return t.kind
end

exactkind(t::Token) = t.kind
startpos(t::Token) = t.startpos
endpos(t::Token) = t.endpos
untokenize(t::Token) = t.val

function untokenize(ts)
    if eltype(ts) != Token
        throw(ArgumentError("element type of iterator has to be Token"))
    end
    io = IOBuffer()
    for tok in ts
        write(io, untokenize(tok))
    end
    return String(take!(io))
end

function Base.show(io::IO, t::Token)
  start_r, start_c = startpos(t)
  end_r, end_c = endpos(t)
  str = kind(t) == ENDMARKER ? "" : escape_string(untokenize(t))
  print(io, rpad(exactkind(t), 15, " "))
  print(io, rpad(str, 20, " "))
  print(io, rpad(string(start_r, ",", start_c, " - ", end_r, ",", end_c), 17, " "))
end

Base.print(io::IO, t::Token) = print(io, untokenize(t))
end # module


text:
module FEFunctions

using Gridap
using Gridap.Triangulations: _attach_triangulation

export FEFunction
export FEFunctionForEval
export interpolate
import Gridap: free_dofs
import Gridap: diri_dofs
import Gridap: FESpace
import Gridap: value_type

import Gridap: evaluate, gradient, return_size
import Gridap: symmetric_gradient
import Base: div
import Gridap: trace
import Gridap: curl
import Gridap: reindex
import Gridap: restrict
import Gridap: Triangulation
import Base: size
import Base: getindex
import Base: +, -, *

import Base: zero

"""
Abstract type representing a FE Function
A FE function is a member of a FESpace.
Since a FE space is the direct sum of free + dirichlet spaces
a FE function is uniquely represented by two vectors of components aka
`free_dofs` and `diri_dofs` which are Vectors{E} with
E = eltype(T) and its fe space
"""
struct FEFunction{
  D,Z,T,E,R,
  C <: IndexCellField{Z,T}} <: IndexCellValue{R,1}
  free_dofs::AbstractVector{E}
  diri_dofs::AbstractVector{E}
  fespace::FESpace{D,Z,T}
  cellfield::C
end

free_dofs(f::FEFunction) = f.free_dofs

diri_dofs(f::FEFunction) = f.diri_dofs

FESpace(f::FEFunction) = f.fespace

Triangulation(f::FEFunction) = Triangulation(f.fespace)

value_type(::FEFunction{D,Z,T}) where {D,Z,T} = T

"""
Returns the FE function represented be the  free and dirichlet values
E = eltype(T)
"""
function FEFunction(
  fespace::FESpace,free_dofs::AbstractVector,diri_dofs::AbstractVector)
  _FEFunction(fespace,free_dofs,diri_dofs)
end

function _FEFunction(
  fespace::FESpace,free_dofs::AbstractVector,diri_dofs::AbstractVector)
  cfield = CellField(fespace,free_dofs,diri_dofs)
  FEFunction(free_dofs,diri_dofs,fespace,cfield)
end

function FEFunction(
  f::FESpaceWithDirichletData,free_vals::AbstractVector{E},diri_dofs::AbstractVector{E}) where E
  _FEFunction(f,free_vals,f.diri_dofs)
end

function FEFunction(f::FESpaceWithDirichletData,free_vals::AbstractVector)
  _FEFunction(f,free_vals,f.diri_dofs)
end

function FEFunctionForEval(
  fespace::FESpaceWithDirichletData,
  free_vals::AbstractVector)

  cfield = CellFieldForEval(fespace,free_vals,fespace.diri_dofs)
  FEFunction(free_vals,fespace.diri_dofs,fespace,cfield)
end

function FEFunction(
  free_dofs::AbstractVector{E},
  diri_dofs::AbstractVector{E},
  fespace::FESpace{D,Z,T},
  cellfield::CellField{Z,T}) where {D,Z,T,E}
  @assert E == eltype(T)
  R = eltype(cellfield)
  C = typeof(cellfield)
  FEFunction{D,Z,T,E,R,C}(free_dofs,diri_dofs,fespace,cellfield)
end

evaluate(f::FEFunction{D,Z},q::CellPoints{Z}) where {D,Z} = evaluate(f.cellfield,q)

for op in (:+,:-,:(gradient),:(symmetric_gradient),:(div),:(trace),:(curl))
  @eval begin
    function ($op)(a::FEFunction)
      g = $op(a.cellfield)
      trian = Triangulation(a)
      _attach_triangulation(g,trian)
    end
  end
end

for op in (:+, :-, :*)
  @eval begin

    function ($op)(a::FEFunction,b::CellField)
      trian = Triangulation(a)
      g = $op(a.cellfield,b)
      _attach_triangulation(g,trian)
    end

    function ($op)(a::CellField,b::FEFunction)
      trian = Triangulation(b)
      g = $op(a,b.cellfield)
      _attach_triangulation(g,trian)
    end

    function ($op)(a::FEFunction,b::Number)
      trian = Triangulation(a)
      g = $op(a.cellfield,b)
      _attach_triangulation(g,trian)
    end

    function ($op)(a::Number,b::FEFunction)
      trian = Triangulation(b)
      g = $op(a,b.cellfield)
      _attach_triangulation(g,trian)
    end

    function ($op)(a::FEFunction,b::Function)
      trian = Triangulation(a)
      cf = CellField(trian,b)
      $op(a,cf)
    end

    function ($op)(a::Function,b::FEFunction)
      trian = Triangulation(b)
      cf = CellField(trian,a)
      $op(cf,b)
    end

  end
end

return_size(f::FEFunction,s::Tuple{Int}) = return_size(f.cellfield,s)

getindex(f::FEFunction,i::Integer) = f.cellfield[i]

size(f::FEFunction) = (length(f.cellfield),)

function interpolate(this::FESpace,fun::Function)
  free_vals, diri_vals = interpolate_values(this,fun)
  FEFunction(this,free_vals,diri_vals)
end
# @santiagobadia : Using this, it always overwrites the FESpace Dirichlet
# values, if they exist. Is this what we want?

function zero(fespace::FESpace{D,Z,T}) where {D,Z,T}
  E = eltype(T)
  nf = num_free_dofs(fespace)
  nd = num_diri_dofs(fespace)
  free_vals = zeros(E,nf)
  diri_vals = zeros(E,nd)
  FEFunction(fespace,free_vals,diri_vals)
end

reindex(uh::FEFunction, indices::CellNumber{<:IndexLike}) = reindex(uh.cellfield,indices)

reindex(uh::FEFunction, indices::IndexCellNumber{<:IndexLike}) = reindex(uh.cellfield,indices)

restrict(uh::FEFunction,trian::BoundaryTriangulation) = restrict(uh.cellfield,trian)

restrict(uh::FEFunction,trian::SkeletonTriangulation) = restrict(uh.cellfield,trian)

# Pretty printing

import Base: show

function show(io::IO,self::FEFunction{D,Z,T,E}) where {D,Z,T,E}
  print(io,"$(nameof(typeof(self))) object")
end

function show(io::IO,::MIME"text/plain",self::FEFunction{D,Z,T,E}) where {D,Z,T,E}
  show(io,self)
  print(io,":")
  print(io,"\n physdim: $D")
  print(io,"\n refdim: $Z")
  print(io,"\n valuetype: $T")
  print(io,"\n doftype: $E")
  print(io,"\n nfree: $(length(free_dofs(self)))")
  print(io,"\n ndiri: $(length(diri_dofs(self)))")
  print(io,"\n ncells: $(length(self))")
end

end # module


text:
using MPIClusterManagers, Distributed
using ArgParse

function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "task"
            help = "Task template specifying the entry wildcard with '{entry}'"
            required = true
        "--entries"
            help = "List of entries to substitute into the task '{entry}' field to generate all tasks."
            required = true
            nargs = '+'
        "--nprocs"
            help = "Number of processes to use."
            arg_type = Int
            required = true
    end

    parse_args(s)
end

function main()
    args = parse_commandline()

    # setup parallelism
    manager = MPIManager(np=args["nprocs"])
    addprocs(manager)

    # make tasks
    tasks = [replace(args["task"], "{entry}" => entry) for entry=args["entries"]]

    # execute tasks in parallel
    @sync @distributed for task in tasks
        run(`sh -c $task`)
    end

    println("Completed all tasks.")
end

main()

max_stars_repo_name: dwillcox/job-run-scripts

max_issues_repo_name: dwillcox/job-run-scripts

max_forks_repo_name: dwillcox/job-run-scripts


text:
const MOUNTAIN_CAR_STATE_MINS = [-1.2, -0.07]
const MOUNTAIN_CAR_STATE_MAXS = [ 0.6,  0.07]
const MOUNTAIN_CAR_ACCELS = [-1.0, 0.0, 1.0]
const MOUNTAIN_CAR_START_STATE = [-0.5, 0.0]

@with_kw struct MountainCar
    γ::Float64 = 1.0
end

n_actions(mdp::MountainCar) = length(MOUNTAIN_CAR_ACCELS)
discount(mdp::MountainCar) = mdp.γ
ordered_actions(mdp::MountainCar) = collect(1:length(MOUNTAIN_CAR_ACCELS))

function generate_start_state(mdp::MountainCar)
    s = MOUNTAIN_CAR_STATE_MINS[1] + rand() * (MOUNTAIN_CAR_STATE_MAXS[1] - MOUNTAIN_CAR_STATE_MINS[1])
    v = MOUNTAIN_CAR_STATE_MINS[2] + rand() * (MOUNTAIN_CAR_STATE_MAXS[2] - MOUNTAIN_CAR_STATE_MINS[2])
    return [s,v]
end

# NOTE: Mountain car is deterministic, so we don't get a transition
#       function that produces a probability distribution.
function mountain_car_transition(s::Vector{Float64}, a::Int)
    x, v = s[1], s[2]
    if x < 0.6
        accel = MOUNTAIN_CAR_ACCELS[a]
        v = v + accel*0.001 + cos(3x)*-0.0025
        x = x + v

        x = clamp(x, MOUNTAIN_CAR_STATE_MINS[1], MOUNTAIN_CAR_STATE_MAXS[1])
        v = clamp(v, MOUNTAIN_CAR_STATE_MINS[2], MOUNTAIN_CAR_STATE_MAXS[2])
    end
    return [x,v]
end

function reward(mdp::MountainCar, s::Vector{Float64}, a::Int)
    if s[1] >= 0.6
        return 0.0
    else
        return -1.0
    end
end

function MDP(mdp::MountainCar; γ::Float64=mdp.γ)
    return MDP(
            γ,
            nothing, # no ordered states
            ordered_actions(mdp),
            (s,a, s′=nothing) -> begin
                # Technically wrong, but its probably fine
                S′ = MvNormal(mountain_car_transition(s,a), 1e-6I)
                if s′ == nothing
                    return S′
                end
                return pdf(S′, s′)
            end,
            (s,a) -> reward(mdp, s, a),
            (s, a)->begin
                s′ = mountain_car_transition(s,a)
                r = reward(mdp, s, a)
                return (s′, r)
            end
        )
end


text:
using SimpleParser

# Parses Go function arguments, and make a TextMate code snipped for setting each argument
# So snippet_from_func_args("(foo int, bacon, bar struct{a int; b string}, egg interface{})")
# Will yield "(${1:foo}, ${2:bacon}, ${3:bar}, ${4:egg})"
function snippet_from_func_args(arg_str::String)
    lexer = Lexer(arg_str)
    parser = Parser(lexer)
    parse_symboltype(parser)
end

function emit_arg(number::Int, name::String)
    "\${$number:$name}"
end

function emit_func_snippet(args::Vector{String})
    "(" * join([emit_arg(i, args[i]) for i in 1:length(args)], ", ") * ")\$0"
end

# Could be func(arg1 type1, arg2 type) ReturnType, or SymbolType
# TODO: this function has mixed up responsibilities since it both parses and emits
function parse_symboltype(parser::Parser)
    result = "\$0"
    if parser.look_ahead.lexeme == "func"
        match(parser, STRING)  # func
        result = emit_func_snippet(parse_arguments(parser))
        match(parser, STRING) # return type
    else
        match(parser, STRING) # symbol type
    end
    result
end

# Parse a list of Go arguments on the form (arg1 type1, arg2 type2)
function parse_arguments(parser::Parser)
    args = String[] :: Vector{String}
    match(parser, '(')
    if look_ahead_type(parser) == ')'
        match(parser, ')')
        return args
    end

    while !(look_ahead_type(parser) in [EOF, RPAREN])
        push!(args, parse_arg(parser))
    end
    match(parser, RPAREN) # same as ')'
    args
end

# Parse a single Go function argument. It will also match any delimeter following
# the argument itself such as , or )
function parse_arg(parser::Parser)
    arg_name = parser.look_ahead.lexeme

    match(parser, STRING)
    if look_ahead_type(parser) == ','
        match(parser, ',')
    elseif look_ahead_type(parser) == ')'
        return arg_name # Don't want to match it because we need it to terminate loop
    else
        parse_type(parser)
        if look_ahead_type(parser) == ','
            match(parser, ',')
        elseif look_ahead_type(parser) == ')'
            return arg_name # Don't want to match it because we need it to terminate loop
        else
            error("Expect token , or ) after argument type")
        end
    end

    arg_name
end

# Parse a Go type. That could be a string like Int64, String or a compound
# type like struct{a int; b int;}. The type is ignored. We don't need the result.
# This function is limited in that it can no deal with nested compound types
function parse_type(parser::Parser)
    keyword = parser.look_ahead.lexeme
    if keyword in ["struct", "interface"]
        match(parser, STRING)
        match(parser, '{')
        t = look_ahead_type(parser)
        while t != '}'
            next_token(parser)
            t = look_ahead_type(parser)
            if t == EOF
                error("No matching } for {")
            end
        end
        match(parser, '}')
    else
        match(parser, STRING)
    end
    nothing
end


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

@doc raw"""Represents a Flocker volume mounted by the Flocker agent. One and only one of datasetName and datasetUUID should be set. Flocker volumes do not support ownership management or SELinux relabeling.

    IoK8sApiCoreV1FlockerVolumeSource(;
        datasetName=nothing,
        datasetUUID=nothing,
    )

    - datasetName::String : Name of the dataset stored as metadata -> name on the dataset for Flocker should be considered as deprecated
    - datasetUUID::String : UUID of the dataset. This is unique identifier of a Flocker dataset
"""
mutable struct IoK8sApiCoreV1FlockerVolumeSource <: SwaggerModel
    datasetName::Any # spec type: Union{ Nothing, String } # spec name: datasetName
    datasetUUID::Any # spec type: Union{ Nothing, String } # spec name: datasetUUID

    function IoK8sApiCoreV1FlockerVolumeSource(;datasetName=nothing, datasetUUID=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1FlockerVolumeSource, Symbol("datasetName"), datasetName)
        setfield!(o, Symbol("datasetName"), datasetName)
        validate_property(IoK8sApiCoreV1FlockerVolumeSource, Symbol("datasetUUID"), datasetUUID)
        setfield!(o, Symbol("datasetUUID"), datasetUUID)
        o
    end
end # type IoK8sApiCoreV1FlockerVolumeSource

const _property_map_IoK8sApiCoreV1FlockerVolumeSource = Dict{Symbol,Symbol}(Symbol("datasetName")=>Symbol("datasetName"), Symbol("datasetUUID")=>Symbol("datasetUUID"))
const _property_types_IoK8sApiCoreV1FlockerVolumeSource = Dict{Symbol,String}(Symbol("datasetName")=>"String", Symbol("datasetUUID")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1FlockerVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1FlockerVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1FlockerVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1FlockerVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1FlockerVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1FlockerVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1FlockerVolumeSource)
    true
end

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


text:
module FEStructure
using SparseArrays
using LinearAlgebra

using HCubature

using ..FESparse
using ..Enums
using ..CoordinateSystem
include("./FEMaterial.jl")
include("./FESection.jl")
include("./FENode.jl")
abstract type AbstractElement end
include("./FEBeam.jl")
include("./FEQuad.jl")
include("./FETria.jl")
#
# using .FEMaterial,.FESection
# using .FENode, .FEBeam, .FEQuad, .FETria

export Structure,AbstractElement,
add_uniaxial_metal!,add_section!,add_general_section!,add_beam_section!,
add_node!,set_nodal_restraint!,set_nodal_spring!,set_nodal_mass!,
add_beam!,set_beam_release!,set_beam_orient!,
add_quad!,set_quad_rotation!,
add_tria!,set_tria_rotation!,

set_damp_constant!,set_damp_Rayleigh!,
integrateK,integrateM,integrateP,
integrateKσ

mutable struct Structure
    csyses::Dict{String,CSys}

    materials::Dict{String,Material}
    sections::Dict{String,BeamSection}

    nodes::Dict{String,Node}
    links::Array
    cables::Array
    beams::Dict{String,Beam}
    quads::Dict{String,Quad}
    trias::Dict{String,Tria}

    damp::String
    ζ₁::Float64
    ζ₂::Float64

    K::SparseMatrixCSC
    K̄::SparseMatrixCSC
    M::SparseMatrixCSC
    M̄::SparseMatrixCSC
    C::SparseMatrixCSC
    C̄::SparseMatrixCSC
end

function Structure()
    globalcsys=CSys([0.,0.,0.],[1.,0.,0.],[0.,1.,0.])
    K=spzeros(1,1)
    K̄=spzeros(1,1)
    M=spzeros(1,1)
    M̄=spzeros(1,1)
    C=spzeros(1,1)
    C̄=spzeros(1,1)
    Structure(Dict("Global"=>globalcsys),Dict{String,Material}(),Dict{String,BeamSection}(),Dict{String,Node}(),[],[],
    Dict{String,Beam}(),Dict{String,Quad}(),Dict{String,Tria}(),"constant",0.05,0,K,K̄,M,M̄,C,C̄)
end

include("./interface.jl")
end