text:
module FemtoCleaner

using GitHub
using GitHub: GitHubAPI, GitHubWebAPI
using HttpCommon
using Deprecations
using CSTParser
using Deprecations: isexpr
using Revise
using MbedTLS
using JSON
using AbstractTrees: children

function with_cloned_repo(f, api::GitHubWebAPI, repo, auth)
    creds = LibGit2.UserPasswordCredentials(String(copy(Vector{UInt8}("x-access-token"))), String(copy(Vector{UInt8}(auth.token))))
    repo_url = "https://github.com/$(get(repo.full_name))"
    local_dir = mktempdir()
    try
        enabled = gc_enable(false)
        lrepo = LibGit2.clone(repo_url, local_dir; payload=Nullable(creds))
        gc_enable(enabled)
        f((lrepo, local_dir))
    finally
        rm(local_dir, force=true, recursive=true)
    end
end

function with_pr_branch(f, repo, auth)
    with_cloned_repo(repo, auth) do x
        LibGit2.branch!(lrepo, "fbot/deps", track=LibGit2.Consts.REMOTE_ORIGIN)
        f(x)
    end
end

function process_deprecations(lrepo, local_dir)
    vers = Pkg.Reqs.parse(joinpath(local_dir, "REQUIRE"))
    deps = Deprecations.applicable_deprecations(vers)
    changed_any = false
    for (root, dirs, files) in walkdir(local_dir)
        for file in files
            fpath = joinpath(root, file)
            endswith(fpath, ".jl") || continue
            # Iterate. Some rewrites may expose others
            while Deprecations.edit_file(fpath, deps)
                changed_any = true
            end
            changed_any && LibGit2.add!(lrepo, relpath(fpath, local_dir))
        end
    end
    changed_any
end

function clone_and_process(api, repo, auth)
    with_cloned_repo(x->process_deprecations(x...), api, repo, auth)
end

function push_repo(api::GitHubWebAPI, repo, auth; force=true)
    creds = LibGit2.UserPasswordCredentials(String(copy(Vector{UInt8}("x-access-token"))), String(copy(Vector{UInt8}(auth.token))))
    enabled = gc_enable(false)
    LibGit2.push(repo, refspecs = ["+HEAD:refs/heads/fbot/deps"], force=force,
        payload=Nullable(creds))
    gc_enable(enabled)
end

function apply_deprecations(api::GitHubAPI, lrepo, local_dir, commit_sig, repo, auth; issue_number = 0)
    changed_any = process_deprecations(lrepo, local_dir)
    if changed_any
        LibGit2.commit(lrepo, "Fix deprecations"; author=commit_sig, committer=commit_sig)
        push_repo(api, lrepo, auth)
    end
    if issue_number != 0
        if changed_any
            create_pull_request(api, repo, auth=auth, params = Dict(
                    :issue => issue_number,
                    :base => get(repo.default_branch),
                    :head => "fbot/deps"
                )
            )
            println("Created pull request for $(GitHub.name(repo))")
        else
            create_comment(api, repo, issue_number, :issue, params = Dict(
                :body => "No applicable deprecations were found in this repository."
            ), auth=auth)
            println("Processing complete for $(GitHub.name(repo)): no changes made")
        end
    else
        if changed_any
            create_pull_request(api, repo, auth=auth, params = Dict(
                    :title => "Fix deprecations",
                    :body => "I fixed a number of deprecations for you",
                    :base => get(repo.default_branch),
                    :head => "fbot/deps"
                )
            )
            println("Created pull request for $(GitHub.name(repo))")
        else
            println("Processing complete for $(GitHub.name(repo)): no changes made")
        end
    end
end

function my_diff_tree(repo::LibGit2.GitRepo, oldtree::LibGit2.GitTree, newtree::LibGit2.GitTree; pathspecs::AbstractString="")
    diff_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @LibGit2.check ccall((:git_diff_tree_to_tree, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{LibGit2.DiffOptionsStruct}),
                   diff_ptr_ptr, repo.ptr, oldtree.ptr, newtree.ptr, isempty(pathspecs) ? C_NULL : pathspecs)
    return LibGit2.GitDiff(repo, diff_ptr_ptr[])
end

function apply_deprecations_if_updated(api::GitHubAPI, lrepo, local_dir, before, after, commit_sig, repo, auth)
    before = LibGit2.GitCommit(lrepo, before)
    after = LibGit2.GitCommit(lrepo, after)
    delta = my_diff_tree(lrepo, LibGit2.peel(before), LibGit2.peel(after); pathspecs="REQUIRE")[1]
    old_blob = LibGit2.GitBlob(lrepo, delta.old_file.id)
    new_blob = LibGit2.GitBlob(lrepo, delta.new_file.id)
    vers_old = Pkg.Reqs.parse(IOBuffer(LibGit2.content(old_blob)))
    vers_new = Pkg.Reqs.parse(IOBuffer(LibGit2.content(new_blob)))
    if vers_new["julia"] != vers_old["julia"]
        apply_deprecations(api, lrepo, local_dir, commit_sig, repo, auth)
    end
end

function dry_run(repo_url)
    local_dir = mktempdir()
    try
        lrepo, changed_any = clone_and_process(local_dir, repo_url)
        cd(local_dir) do
            run(`git status`)
            run(`git diff --cached`)
        end
    finally
        println(local_dir)
        #rm(local_dir, force=true, recursive=true)
    end
end

include("interactions.jl")
include("autodeployment.jl")

const autodeployment_enabled = haskey(ENV, "FEMTOCLEANER_AUTODEPLOY") ?
    ENV["FEMTOCLEANER_AUTODEPLOY"] == "yes" : false

function event_callback(api::GitHubAPI, app_name, app_key, app_id, sourcerepo_installation,
                        commit_sig, listener, bug_repository, event)
    # On installation, process every repository we just got installed into
    if event.kind == "installation"
        jwt = GitHub.JWTAuth(app_id, app_key)
        installation = Installation(event.payload["installation"])
        auth = create_access_token(api, installation, jwt)
        for repo in event.payload["repositories"]
            repo = GitHub.repo(api, GitHub.Repo(repo); auth=auth)
            with_cloned_repo(api, repo, auth) do x
                apply_deprecations(api, x..., commit_sig, repo, auth)
            end
        end
    elseif event.kind == "installation_repositories"
        jwt = GitHub.JWTAuth(app_id, app_key)
        installation = Installation(event.payload["installation"])
        auth = create_access_token(api, installation, jwt)
        for repo in event.payload["repositories_added"]
            repo = GitHub.repo(api, GitHub.Repo(repo); auth=auth)
            with_cloned_repo(api, repo, auth) do x
                apply_deprecations(api, x..., commit_sig, repo, auth)
            end
        end
    elseif event.kind == "pull_request_review"
        jwt = GitHub.JWTAuth(app_id, app_key)
        pr_response(api, event, jwt, commit_sig, app_name, sourcerepo_installation, bug_repository)
    elseif event.kind == "push"
        jwt = GitHub.JWTAuth(app_id, app_key)
        installation = Installation(event.payload["installation"])
        auth = create_access_token(api, installation, jwt)
        repo = Repo(event.payload["repository"])
        # Check if REQUIRE was updated
        for commit in event.payload["commits"]
            if "REQUIRE" in commit["modified"]
                with_cloned_repo(api, repo, auth) do x
                    apply_deprecations_if_updated(api, x...,
                        event.payload["before"], event.payload["after"],
                        commit_sig, repo, auth)
                end
                break
            end
        end
        maybe_autdodeploy(event, listener, jwt, sourcerepo_installation, autodeployment_enabled)
    elseif event.kind == "issues" && event.payload["action"] == "opened"
        jwt = GitHub.JWTAuth(app_id, app_key)
        iss = Issue(event.payload["issue"])
        repo = Repo(event.payload["repository"])
        installation = Installation(event.payload["installation"])
        auth = create_access_token(api, installation, jwt)
        if lowercase(get(iss.title)) == "run femtocleaner"
            with_cloned_repo(api, repo, auth) do x
                apply_deprecations(api, x..., commit_sig, repo, auth; issue_number = get(iss.number))
            end
        end
    end
    return HttpCommon.Response(200)
end

function run_server()
    app_key = MbedTLS.PKContext()
    MbedTLS.parse_key!(app_key, haskey(ENV, "FEMTOCLEANER_PRIVKEY") ? ENV["FEMTOCLEANER_PRIVKEY"] : readstring(joinpath(dirname(@__FILE__),"..","privkey.pem")))
    app_id = parse(Int, strip(haskey(ENV, "FEMTOCLEANER_APPID") ? ENV["FEMTOCLEANER_APPID"] : readstring(joinpath(dirname(@__FILE__),"..","app_id"))))
    secret = haskey(ENV, "FEMTOCLEANER_SECRET") ? ENV["FEMTOCLEANER_SECRET"] : nothing
    sourcerepo_installation = haskey(ENV, "FEMTOCLEANER_INSTALLATION") ? parse(Int, ENV["FEMTOCLEANER_INSTALLATION"]) : 0
    bug_repository = haskey(ENV, "FEMTOCLEANER_BUGREPO") ? ENV["FEMTOCLEANER_BUGREPO"] : ""
    (secret == nothing) && warn("Webhook secret not set. All events will be accepted. This is an insecure configuration!")
    jwt = GitHub.JWTAuth(app_id, app_key)
    app_name = get(GitHub.app(; auth=jwt).name)
    commit_sig = LibGit2.Signature("$(app_name)[bot]", "$(app_name)[bot]@users.noreply.github.com")
    api = GitHub.DEFAULT_API
    @async update_existing_repos(api, commit_sig, app_id, app_key)
    local listener
    listener = GitHub.EventListener(secret=secret) do event
        revise()
        Base.invokelatest(event_callback, api, app_name, app_key, app_id,
                          sourcerepo_installation, commit_sig, listener, bug_repository, event)
    end
    GitHub.run(listener, host=IPv4(0,0,0,0), port=10000+app_id)
    wait()
end

end # module


text:
using LiveServer
cd("example")
serve()

ext: jl


text:
# --------------------------------------------------------------------------- #
# Unary transformators
# --------------------------------------------------------------------------- #

function eval(block::B) where {B <: AbstractUnaryTransformer}
    error("not implemented yet")
end

function eval(block::B, time) where {B <: AbstractUnaryTransformer}
    return eval(block)
end

mutable struct Identity{T} <: AbstractUnaryTransformer
    x::InputPin{T}
    y::OutputPin{T}

    function Identity{T}() where {T}
        x = InputPin{T}()
        y = OutputPin{T}()

        obj = new()

        obj.x = x
        obj.y = y

        x.containing_object = obj
        y.containing_object = obj

        return obj
    end
end
function Identity()
    T = Float64
    return Identity{T}()
end

function eval(block::Identity)
    val = block.x.value
    block.y.value = val
    return val
end

mutable struct Gain{T} <: AbstractUnaryTransformer
    g::T
    x::InputPin{T}
    y::OutputPin{T}

    function Gain{T}(g::T) where {T}
        x = InputPin{T}()
        y = OutputPin{T}()

        obj = new()

        obj.x = x
        obj.y = y
        obj.g = g

        x.containing_object = obj
        y.containing_object = obj

        return obj
    end
end
function Gain(g::T) where {T}
    return Gain{T}(g)
end
function Gain(g = one(Float64))
    return Gain{Float64}(g)
end

function eval(block::Gain)
    val = block.x.value * block.g
    block.y.value = val
    return val
end

mutable struct Sin{T} <: AbstractUnaryTransformer
    x::InputPin{T}
    y::OutputPin{T}

    function Sin{T}() where {T}
        x = InputPin{T}()
        y = OutputPin{T}()

        obj = new()

        obj.x = x
        obj.y = y

        x.containing_object = obj
        y.containing_object = obj

        return obj
    end
end
function Sin()
    return Sin{Float64}()
end

function eval(block::Sin)
    val = sin(block.x.value)
    block.y.value = val
    return val
end

mutable struct UserDefinedUnaryFunction{T} <: AbstractUnaryTransformer
    fcn
    x::InputPin{T}
    y::OutputPin{T}

    function UserDefinedUnaryFunction{T}(fcn) where {T}
        x = InputPin{T}()
        y = OutputPin{T}()

        obj = new()

        obj.x = x
        obj.y = y
        obj.fcn = fcn

        x.containing_object = obj
        y.containing_object = obj

        return obj
    end
end
function UserDefinedUnaryFunction(fcn) where {T}
    return UserDefinedUnaryFunction{Float64}(fcn)
end

function eval(block::UserDefinedUnaryFunction)
    val = block.fcn(block.x.value)
    block.y.value = val
    return val
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using Test
using ClimateMachine

let
    for N in 1:10
        for FloatType in (Float64, Float32)
            (ξ, ω) = ClimateMachine.Mesh.Elements.lglpoints(FloatType, N)
            I∫ =
                ClimateMachine.Mesh.Grids.indefinite_integral_interpolation_matrix(
                    ξ,
                    ω,
                )
            for n in 1:N
                if N == 1
                    @test sum(abs.(I∫ * ξ)) < 10 * eps(FloatType)
                else
                    @test I∫ * ξ .^ n ≈
                          (ξ .^ (n + 1) .- (-1) .^ (n + 1)) / (n + 1)
                end
            end
        end
    end
end


text:
# Utility to make it easier to define get/setproperty for a long list of options.

# build up the body of get/setproperty from a list of properties, of the form
#     if name === $prop
#         ($doit(prop))
#     elseif ...
# for getproperty(value, name) or setproperty!(value, name, x)
function propexpression(doit::Function, properties)
    ex = :(error(typeof(value), " has no field ", name))
    for prop in properties
        ex = Expr(:elseif, :(name === $(QuoteNode(prop))), doit(prop), ex)
    end
    return Expr(:if, ex.args...)
end


text:
module aoc_12

export Ship, north, south, east, west, left, right, forward
export distance

abstract type Positionable end

mutable struct Ship
    direction::Int64
    x::Float64
    y::Float64
end

Ship() = Ship( 0, 0, 0 )

north( ship::Ship, dist ) = ship.y += dist
east( ship::Ship, dist )  = ship.x += dist
left( ship::Ship, degrees ) = ship.direction += degrees

south( s, d ) = north( s, -d)
west( s, d ) = east( s, -d)
right( s, d) = left( s, -d )

function forward( ship::Ship, dist )
    ship.x += cosd( ship.direction ) * dist
    ship.y += sind( ship.direction ) * dist
end

distance( ship::Ship ) = abs( ship.x ) + abs( ship.y )

const commands = Dict(
    'N' => north,
    'S' => south,
    'E' => east,
    'W' => west,
    'L' => left,
    'R' => right,
    'F' => forward
)

const regexp = r"([NSEWLRF])([0-9]+)"

function read_instructions( filename , ship::Union{Ship,ShipWithWaypoint} = Ship() )
    open( filename, "r") do io
        for line in eachline( filename )
            m = match( regexp, line )
            commands[ m.captures[1][1] ]( ship, parse( Int64, m.captures[2] ) )
        end
    end
    println( distance(ship) )
    return ship
end

# Part 2
export ShipWithWaypoint

mutable struct Waypoint
    x::Float64
    y::Float64
end
Waypoint() = Waypoint(10,1)
mutable struct ShipWithWaypoint
    s::Ship
    w::Waypoint
end
ShipWithWaypoint() = ShipWithWaypoint(Ship(),Waypoint())
north(s::ShipWithWaypoint, dist) = s.w.y += dist
east(s::ShipWithWaypoint, dist) = s.w.x += dist
function left(s::ShipWithWaypoint, degrees)
    radius = hypot(s.w.y, s.w.x)
    angle = atand(s.w.y, s.w.x )
    angle += degrees
    s.w.x = cosd( angle ) * radius
    s.w.y = sind( angle ) * radius
end
function forward(s::ShipWithWaypoint, dist)
    s.s.x += s.w.x * dist
    s.s.y += s.w.y * dist
end
distance(s::ShipWithWaypoint) = distance( s.s )

end


text:
# RemoteChannel used for progress monitoring in parallel setup
const progch = RemoteChannel(()->Channel{Bool}(), 1)

"""
    mc(scenario::Scenario,
      params::Vector{<:Pair},
      num_iter::Int64;
      verbose=false,
      alg=DEFAULT_ALG,
      reltol=DEFAULT_SIMULATION_RELTOL,
      abstol=DEFAULT_SIMULATION_ABSTOL,
      parallel_type=EnsembleSerial(),
      kwargs...
    )

Run Monte-Carlo simulations with single `Scenario`. Returns [`MCResult`](@ref) type.

Example: `mc(scenario, [:k2=>Normal(1e-3,1e-4), :k3=>Uniform(1e-4,1e-2)], 1000)`

Arguments:

- `scenario` : simulation scenario of type [`Scenario`](@ref)
- `params` : parameters variation setup
- `num_iter` : number of Monte-Carlo iterations
- `verbose` : print iteration progress. Default is `false`
- `progress_bar` : show progress bar. Default is `false`
- `alg` : ODE solver. See SciML docs for details. Default is AutoTsit5(Rosenbrock23())
- `reltol` : relative tolerance. Default is 1e-3
- `abstol` : relative tolerance. Default is 1e-6
- `parallel_type` : parallel setup. See SciML docs for details. Default is no parallelism: EnsembleSerial()
- kwargs : other solver related arguments supported by DiffEqBase.solve. See SciML docs for details
"""
function mc(
  scenario::Scenario,
  params::Vector{P},
  num_iter::Int;
  verbose=false,
  progress_bar=false,
  alg=DEFAULT_ALG,
  reltol=DEFAULT_SIMULATION_RELTOL,
  abstol=DEFAULT_SIMULATION_ABSTOL,
  parallel_type=EnsembleSerial(),
  kwargs...
) where P<:Pair

  prob0 = scenario.prob
  init_func = scenario.init_func
  params_nt = NamedTuple(params)

  #(parallel_type == EnsembleSerial()) # tmp fix
  p = Progress(num_iter, dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, enabled = progress_bar)

  function prob_func(prob,i,repeat)
    verbose && println("Processing iteration $i")
    progress_bar && (parallel_type != EnsembleDistributed() ? next!(p) : put!(progch, true))
    update_init_values(prob, init_func, generate_cons(params_nt,i))
  end

  params_names = collect(keys(params_nt))
  function output_func(sol, i)
    sim = build_results(sol, params_names)
    (sim, false)
  end

  prob = EnsembleProblem(prob0;
    prob_func = prob_func,
    output_func = output_func,
    #reduction = reduction_func
  )

  if progress_bar && (parallel_type == EnsembleDistributed())
    @sync begin
      @async while take!(progch)
        next!(p)
      end
      @async begin
        solution = solve(prob, alg, parallel_type;
          trajectories = num_iter,
          reltol = reltol,
          abstol = abstol,
          save_start = false,
          save_end = false,
          save_everystep = false,
          kwargs...
        )
        put!(progch, false)
      end
    end
  else
    solution = solve(prob, alg, parallel_type;
      trajectories = num_iter,
      reltol = reltol,
      abstol = abstol,
      save_start = false,
      save_end = false,
      save_everystep = false,
      kwargs...
    )
  end

  return MCResult(solution.u, has_saveat(scenario), scenario)
end

"""
    mc(scenario::Scenario,
      params::DataFrame,
      num_iter::Int64;
      kwargs...
    )

Run Monte-Carlo simulations with single scenario `Scenario`. Returns [`MCResult`](@ref) type.

Example: `mc(scn1, DataFrame(k2=rand(3),k3=rand(3)), 1000)`

Arguments:

- `scenario` : simulation scenario of type [`Scenario`](@ref)
- `params` : DataFrame with pre-generated parameters.
- `num_iter` : number of Monte-Carlo iterations
- kwargs : other solver related arguments supported by `mc(scenario::Scenario, params::Vector, num_iter::Int64)`
"""
function mc(
  scenario::Scenario,
  params::DataFrame;
  num_iter::Int= size(params)[1],
  kwargs...
)
  cons = keys(parameters(scenario))
  params_pairs = Pair[]

  for pstr in names(params)
    psym = Symbol(pstr)
    @assert (psym in cons) "$psym is not found in models constants."
    push!(params_pairs, psym=>params[!,psym])
  end

  return mc(scenario, params_pairs, num_iter; kwargs...)
end

"""
    mc(model::Model,
      params::Vector{<:Pair},
      num_iter::Int64;
      measurements::Vector{AbstractMeasurementPoint} = AbstractMeasurementPoint[],
      events_active::Union{Nothing, Vector{Pair{Symbol,Bool}}} = Pair{Symbol,Bool}[],
      events_save::Union{Tuple,Vector{Pair{Symbol, Tuple{Bool, Bool}}}}=(true,true),
      observables::Union{Nothing,Vector{Symbol}} = nothing,
      saveat::Union{Nothing,AbstractVector} = nothing,
      tspan::Union{Nothing,Tuple} = nothing,
      save_scope::Bool=false,
      time_type::DataType=Float64,
      kwargs...
    )

Run Monte-Carlo simulations with `Model`. Returns [`MCResult`](@ref) type.

Example: `mc(model, [:k2=>Normal(1e-3,1e-4), :k3=>Uniform(1e-4,1e-2)], 1000)`

Arguments:

- `model` : model of type [`Model`](@ref)
- `params` : parameters variation setup
- `num_iter` : number of Monte-Carlo iterations
- `measurements` : `Vector` of measurements. Default is empty vector
- `events_active` : `Vector` of `Pair`s containing events' names and true/false values. Overwrites default model's values. Default is empty vector
- `events_save` : `Tuple` or `Vector{Tuple}` marking whether to save solution before and after event. Default is `(true,true)` for all events
- `observables` : names of output observables. Overwrites default model's values. Default is empty vector
- `saveat` : time points, where solution should be saved. Default `nothing` values stands for saving solution at timepoints reached by the solver
- `tspan` : time span for the ODE problem
- `save_scope` : should scope be saved together with solution. Default is `false`
- kwargs : other solver related arguments supported by `mc(scenario::Scenario, params::Vector, num_iter::Int64)`
"""
function mc(
  model::Model,
  params::Vector{P},
  num_iter::Int;

  ## arguments for Scenario(::Model,...)
  measurements::Vector{AbstractMeasurementPoint} = AbstractMeasurementPoint[],
  events_active::Union{Nothing, Vector{Pair{Symbol,Bool}}} = Pair{Symbol,Bool}[],
  events_save::Union{Tuple,Vector{Pair{Symbol, Tuple{Bool, Bool}}}}=(true,true),
  observables::Union{Nothing,Vector{Symbol}} = nothing,
  saveat::Union{Nothing,AbstractVector} = nothing,
  tspan::Union{Nothing,Tuple} = nothing,
  save_scope::Bool=false,
  time_type::DataType=Float64,

  kwargs...
) where P<:Pair

  scenario = Scenario(
    model; measurements,
    events_active, events_save, observables, saveat, tspan, save_scope, time_type)

  return mc(scenario, params, num_iter; kwargs...)
end

# multi scenario Monte-Carlo

"""
    mc(scenario_pairs::Vector{<:Pair},
      params::Vector{<:Pair},
      num_iter::Int64;
      verbose=false,
      alg=DEFAULT_ALG,
      reltol=DEFAULT_SIMULATION_RELTOL,
      abstol=DEFAULT_SIMULATION_ABSTOL,
      parallel_type=EnsembleSerial(),
      kwargs...
    )

Run Monte-Carlo simulations with single `Scenario`. Returns `Vector{MCResult}` type.

Example: `mc([:c1=>scn1,:c2=>scn2], [:k2=>Normal(1e-3,1e-4), :k3=>Uniform(1e-4,1e-2)], 1000)`

Arguments:

- `scenario_pairs` : vector of pairs containing names and scenarios of type [`Scenario`](@ref)
- `params` : parameters variation setup
- `num_iter` : number of Monte-Carlo iterations
- `verbose` : print iteration progress. Default is `false`
- `progress_bar` : show progress bar. Default is `false`
- `alg` : ODE solver. See SciML docs for details. Default is AutoTsit5(Rosenbrock23())
- `reltol` : relative tolerance. Default is 1e-3
- `abstol` : relative tolerance. Default is 1e-6
- `parallel_type` : parallel setup. See SciML docs for details. Default is no parallelism: EnsembleSerial()
- kwargs : other solver related arguments supported by DiffEqBase.solve. See SciML docs for details
"""
function mc(
  scenario_pairs::Vector{CP},
  params::Vector{PP},
  num_iter::Int;
  verbose=false,
  progress_bar=false,
  alg=DEFAULT_ALG,
  reltol=DEFAULT_SIMULATION_RELTOL,
  abstol=DEFAULT_SIMULATION_ABSTOL,
  parallel_type=EnsembleSerial(),
  kwargs...
) where {CP<:Pair, PP<:Pair}

  params_nt = NamedTuple(params)
  params_pregenerated = [generate_cons(params_nt,i) for i in 1:num_iter]
  lp = length(params_pregenerated)
  lc = length(scenario_pairs)
  iter = collect(Iterators.product(1:lp,1:lc))

  p = Progress(num_iter, dt=0.5, barglyphs=BarGlyphs("[=> ]"), barlen=50, enabled=progress_bar)

  function prob_func(prob,i,repeat)
    iter_i = iter[i]
    verbose && println("Processing scenario $(iter_i[2]) iteration $(iter_i[1])")
    progress_bar && (parallel_type != EnsembleDistributed() ? next!(p) : put!(progch, true))
    prob_i = last(scenario_pairs[iter_i[2]]).prob
    init_i = last(scenario_pairs[iter_i[2]]).init_func
    update_init_values(prob_i, init_i, params_pregenerated[iter_i[1]])
  end

  params_names = collect(keys(params_nt))
  function output_func(sol, i)
    sim = build_results(sol, params_names)
    (sim, false)
  end

  prob = EnsembleProblem(last(scenario_pairs[1]).prob;
    prob_func = prob_func,
    output_func = output_func,
    #reduction = reduction_func
  )

  if progress_bar && (parallel_type == EnsembleDistributed())
    @sync begin
      @async while take!(progch)
        next!(p)
      end
      @async begin
        solution = solve(prob, alg, parallel_type;
          trajectories = lp*lc,
          reltol = reltol,
          abstol = abstol,
          save_start = false,
          save_end = false,
          save_everystep = false,
          kwargs...
        )
        put!(progch, false)
      end
    end
  else
    solution = solve(prob, alg, parallel_type;
      trajectories = lp*lc,
      reltol = reltol,
      abstol = abstol,
      save_start = false,
      save_end = false,
      save_everystep = false,
      kwargs...
    )
  end

  ret = Vector{Pair{Symbol,MCResult}}(undef, lc)

  for i in 1:lc
    ret[i] = first(scenario_pairs[i]) =>
      MCResult(solution.u[lp*(i-1)+1:i*lp], has_saveat(last(scenario_pairs[i])), last(scenario_pairs[i]))
  end
  return ret
end

"""
    mc(scenario_pairs::Vector{<:AbstractScenario},
      params::Vector{<:Pair},
      num_iter::Int64;
      kwargs...
    )

Run Monte-Carlo simulations with single scenario. Returns `Vector{MCResult}` type.

Example: `mc([scn1,scn2], [:k2=>Normal(1e-3,1e-4), :k3=>Uniform(1e-4,1e-2)], 1000)`

Arguments:

- `scenario_pairs` : vector of scenarios of type [`Scenario`](@ref)
- `params` : parameters variation setup
- `num_iter` : number of Monte-Carlo iterations
- kwargs : other solver related arguments supported by `mc(scenario_pairs::Vector{<:Pair}, params::Vector, num_iter::Int64)`
"""
function mc(
  scenario_pairs::Vector{C},
  params::Vector{P},
  num_iter::Int;
  kwargs...
) where {C<:AbstractScenario, P<:Pair}

  scenario_pairs = [(Symbol("_$i") => scn) for (i, scn) in pairs(scenario_pairs)]
  return mc(scenario_pairs, params, num_iter; kwargs...)
end

"""
    mc(platform::Platform,
      params::Vector{<:Pair},
      num_iter::Int64;
      kwargs...
    )

Run Monte-Carlo simulations with single `Scenario`. Returns `Vector{MCResult}` type.

Example: `mc(platform, [:k2=>Normal(1e-3,1e-4), :k3=>Uniform(1e-4,1e-2)], 1000)`

Arguments:

- `platform` : platform of [`Platform`](@ref) type
- `params` : parameters variation setup
- `num_iter` : number of Monte-Carlo iterations
- kwargs : other solver related arguments supported by `mc(scenario_pairs::Vector{<:Pair}, params::Vector, num_iter::Int64)`
"""
function mc(
  platform::Platform,
  params::Vector{P},
  num_iter::Int;
  scenarios::Union{AbstractVector{Symbol}, Nothing} = nothing,
  kwargs...) where P<:Pair

  if isnothing(scenarios)
    scenario_pairs = [platform.scenarios...]
  else
    scenario_pairs = Pair{Symbol,AbstractScenario}[]
    for scn_name in scenarios
      @assert haskey(platform.scenarios, scn_name) "No scenario :$scn_name found in the platform."
      push!(scenario_pairs, scn_name=>platform.scenarios[scn_name])
    end
  end

  return mc(scenario_pairs,params,num_iter;kwargs...)
end

########################################## Statistics ######################################################

# currently median and quantile don't output LVector

function DiffEqBase.EnsembleAnalysis.get_timestep(mcr::MCResult,i)
  @assert has_saveat(mcr) "Solution doesn't contain single time vector, default statistics are not available."
  return (getindex(mcr[j],i) for j in 1:length(mcr))
end

function DiffEqBase.EnsembleAnalysis.get_timepoint(mcr::MCResult,t)
  @assert has_saveat(mcr) "Solution doesn't contain single time vector, default statistics are not available."
  return (mcr[j](t) for j in 1:length(mcr))
end

function DiffEqBase.EnsembleAnalysis.EnsembleSummary(sim::MCResult,
  t=sim[1].t;quantiles=[0.05,0.95])

  m,v = timeseries_point_meanvar(sim,t)
  qlow = timeseries_point_quantile(sim,quantiles[1],t)
  qhigh = timeseries_point_quantile(sim,quantiles[2],t)

  trajectories = length(sim)

  EnsembleSummary{Float64,2,typeof(t),typeof(m),typeof(v),typeof(qlow),typeof(qhigh)}(t,m,v,qlow,qhigh,trajectories,0.0,true)
end

generate_cons(vp::AbstractVector{P},i)  where P<:Pair = NamedTuple([k=>generate_cons(v,i) for (k,v) in vp])
generate_cons(nt::NamedTuple,i) = NamedTuple{keys(nt)}([generate_cons(v,i) for v in nt])
generate_cons(v::Distribution,i) = rand(v)
generate_cons(v::Real,i) = v
generate_cons(v::AbstractVector{R},i) where R<:Float64 = v[i]

"""
    read_mcvecs(filepath::String)

Read table with pre-generated parameters for Monte-Carlo simulations. Typically used for virtual patients simulations

Arguments:

- `filepath`: path to the file with pre-generated parameters
"""
read_mcvecs(filepath::String) = DataFrame(CSV.File(filepath; typemap=Dict(Int => Float64)))


text:
using CSV
using DataFrames
using Printf
using Statistics

const DIR = @__DIR__
const ROOT_DIR = joinpath(@__DIR__, "..")

"""
    gmean1(u)

Compute the geometric mean with a shift of 1.
"""
gmean1(u) = exp(mean(log.(u .+ oneunit(eltype(u))))) - oneunit(eltype(u))

# results file name
res_file_name = ARGS[1]

# Import results
df = CSV.read(res_file_name)
# Import instance sizes
df_sizes = CSV.read(ROOT_DIR * "/dat/instances.csv")

# Compute aggregated statistics
# First, compute sets of instances
N = length(df.instance)
N_all = trues(N)
N_solved = trues(N)  # instances solved by all solvers
for s in [:CLP, :ECOS, :GLPK, :IPOPT, :OOQP, :TLP, :CPX, :GRB, :MSK]
    N_solved .*= df[Symbol("success_$(s)")]
end

println("$N instances")
println("$(sum(N_solved)) solved by all solvers")

N_small = (df_sizes.nrows + df_sizes.ncols) .< 1e3
N_medium = (1e3 .<= (df_sizes.nrows + df_sizes.ncols) .< 1e4)
N_large = 1e4 .<= (df_sizes.nrows + df_sizes.ncols)

println("$(sum( N_small))  small instances")
println("$(sum(N_medium)) medium instances")
println("$(sum( N_large))  large instances")

N_small .*= N_solved
N_medium .*= N_solved
N_large .*= N_solved

println("$(sum( N_small))  small  instances (solved)")
println("$(sum(N_medium)) medium  instances (solved)")
println("$(sum(N_large))   large  instances (solved)")

println("\n\n")

# Open-source solvers
println("\\toprule")
print("        \t &")
for s in [:CLP, :ECOS, :GLPK, :IPOPT, :OOQP, :TLP]
    @printf "& %6s \t " s
end
println("\\\\")
println("\\midrule")

print("Solved  \t &")
for s in [:CLP, :ECOS, :GLPK, :IPOPT, :OOQP, :TLP]
    @printf "& \$%6d\$\t " sum(df[Symbol("success_$s")][N_all])
end
println("\\\\")
println("\\midrule")

# Computing times
for (setname, pbset) in zip(
        ["N-$N", "N-$(sum(N_solved))", "N-small", "N-medium", "N-large"],
        [N_all, N_solved, N_small, N_medium, N_large]
    )
    @printf "%8s\t &" setname
    for s in [:CLP, :ECOS, :GLPK, :IPOPT, :OOQP, :TLP]
        @printf "& \$%6.3f\$\t " gmean1(df[Symbol("time_$s")][pbset])
    end
    println("\\\\")
end
println("\\bottomrule")

println("\n\n\n")

# Commercial solvers
println("\\toprule")
print("        \t &")
for s in [:CPX, :GRB, :MSK]
    @printf "& %6s\t " s
end
println("\\\\")
println("\\midrule")

print("Solved  \t &")
for s in [:CPX, :GRB, :MSK]
    @printf "& \$%6d\$\t " sum(df[Symbol("success_$s")][N_all])
end
println("\\\\")
println("\\midrule")

# Computing times
for (setname, pbset) in zip(
        ["N-$N", "N-$(sum(N_solved))", "N-small", "N-medium", "N-large"],
        [N_all, N_solved, N_small, N_medium, N_large]
    )
    @printf "%8s\t &" setname
    for s in [:CPX, :GRB, :MSK]
        @printf "& \$%6.3f\$\t " gmean1(df[Symbol("time_$s")][pbset])
    end
    println("\\\\")
end
println("\\bottomrule")


text:
const DASH = '-'
const DASHSTR = "-"

if @isdefined stdout
    const StdOut = Base.stdout
else
    const StdOut = Base.STDOUT
end

function splitstring(str::AbstractString, splitat::AbstractString)
    a, z = split(str, splitat)
    return String(a), String(z)
end

function string(td::TimeDate)
    return string(on_date(td),"T",at_time(td))
end

function show(io::IO, td::TimeDate)
    str = string(td)
    print(io, str)
end

function show(td::TimeDate)
    str = string(td)
    print(str)
end

timezonename(tdz::TimeDateZone) = string(tdz.inzone)
timezonename(zdt::ZonedDateTime) = string(zdt.timezone)
timezonename(tz::TimeZone) = string(tz)

string(tdz::TimeDateZone; tzname::Bool=false) =
    tzname ? stringwithzone(tdz) : stringwithoffset(tdz)

function stringwithzone(tdz::TimeDateZone)
    tzname = timezonename(tdz.inzone)
    return string(on_date(tdz),"T",at_time(tdz)," ",tzname)
end

function stringwithoffset(tdz::TimeDateZone)
    offset = string(ZonedDateTime(tdz))[end-5:end]
    if offset[2:end] == "00:00"
        offset = "Z"
    end
    return string(on_date(tdz),"T",at_time(tdz),offset)
end

function showwithoffset(io::IO, tdz::TimeDateZone)
    str = stringwithoffset(tdz)
    print(io, str)
end

function showwithoffset(tdz::TimeDateZone)
    str = stringwithoffset(tdz)
    print(str)
end

function showwithzone(io::IO, tdz::TimeDateZone)
    str = stringwithzone(tdz)
    print(io, str)
end

function showwithzone(tdz::TimeDateZone)
    str = stringwithzone(tdz)
    print(str)
end

show(io::IO, tdz::TimeDateZone; tzname::Bool=false) =
    tzname ? showwithzone(io, tdz) : showwithoffset(io,tdz)
show(tdz::TimeDateZone; tzname::Bool=false) =
    tzname ? showwithzone(tdz) : showwithoffset(tdz)

splitstring(str::String, splitat::String) = map(String, split(str, splitat))

function TimeDate(str::String)
    !occursin("T", str) && throw(ErrorException("\"$str\" is not recognized as a TimeDate"))

    datepart, inttimepart, fractimepart = datetimeparts(str)

    ondate = parse(Date, datepart)
    attime = parse(Time, inttimepart)
    attime = fractionaltime(attime, fractimepart)

    return TimeDate(attime, ondate)
end

function TimeDateZone(str::String)
    !occursin("T", str) && throw(ErrorException("\"$str\" is not recognized as a TimeDateZone"))

    if occursin(" ", str)
        timedatestr, tzname = splitstring(str, " ")
        timedate    = TimeDate(timedatestr)
        tz          = TimeZone(tzname)
        TimeDateZone(timedate, tz)
    else
        if str[end] == 'Z'
            timedatestr = str[1:end-1]
            tz = tz"UTC"
        else
            timedatestr = str[1:end-6]
            tzoffsetstr = str[end-5:end]
            tz = TimeZone(tzoffsetstr)
        end
        timedate = TimeDate(timedatestr)
        tm = timedate.attime
        dt = timedate.ondate
        TimeDateZone(tm, dt, tz)
    end
end

function fractionaltime(attime::Time, fractimepart::String)
    n = length(fractimepart)
    if n > 0
        fractime = parse(Int, fractimepart)
        if n <= 3
            delta = fld(1000,10^n)
            fractime *= delta
            attime = attime + Millisecond(fractime)
        elseif n <= 6
            delta = fld(1000,10^(n-3))
            fractime *= delta
            attime = attime + Microsecond(fractime)
        else
            delta = fld(1000,10^(n-6))
            fractime *= delta
            attime = attime + Nanosecond(fractime)
        end
    end
    return attime
end

function datetimeparts(str::String)
    datepart, timepart = splitstring(str, "T")
    if occursin(".", timepart)
        inttimepart, fractimepart = splitstring(timepart, ".")
    else
        inttimepart = timepart
        fractimepart = ""
    end

    return datepart, inttimepart, fractimepart
end


text:
function get_module_3_excersise_2_author_name()
    return "YOUR NAME HERE"
end

function problem_2()
    c = []
    A = []
    x = []
    return x #return the x vector
end
println(problem_2())


text:
### A Pluto.jl notebook ###
# v0.17.5

using Markdown
using InteractiveUtils

# ╔═╡ b8f90c8e-788d-11ec-0f86-c3a393796a2d
using CUDA

# ╔═╡ 40b2c8ff-8c11-47f7-9dd5-19efb9663b0e
md"""
# Defining Data
"""

# ╔═╡ c762d0a0-cb54-47e4-9ad8-51eac1cd437c
# Length of vector
L = 10

# ╔═╡ a15b15e1-cbf0-4843-98dd-3dc6d9fdd41e
# Random Vector A of length L
A = Vector{Float32}(rand(L) * 15)

# ╔═╡ 51c98125-2ca5-4570-819a-162a172a5390
# Random Vector B of length L
B = Vector{Float32}(rand(L) * 8)

# ╔═╡ a0c6094a-dc63-433e-a122-44914c90d200
md"""
# Vector Addition
"""

# ╔═╡ 071054b2-c3d4-4d0f-9f6c-84e6f62246d9
md"""
## Host code
"""

# ╔═╡ 4e168101-9f5f-4a2b-89d1-6a5e46178382
md"""
## GPU Kernel
"""

# ╔═╡ ee948bcb-66e3-41a1-9795-7d8e7b57d806
function vecAddKernel!(A, B, C, L)
	idx = blockDim().x * (blockIdx().x-1) + threadIdx().x

	if idx ≤ L
		C[idx] = A[idx] + B[idx]
	end
	return nothing
end

# ╔═╡ 077f54d1-1bde-47ea-9192-291c88109d77
function vecAdd(A::Vector, B::Vector, L::Int)
	# Allocating and transfering data to GPU memory
	d_A = CuArray(A)
	d_B = CuArray(B)
	d_C = CUDA.zeros(L) # Allocating memory for result

	# Kernel Invocation
	dimGrid = (Int(ceil(L/256)),1,1)
	dimBlock = (256,1,1)
	@cuda blocks=dimGrid threads=dimBlock vecAddKernel!(d_A, d_B, d_C, L)

	C = Array(d_C)

	return C
end

# ╔═╡ 622542f8-1417-4cca-8308-f2b6ea61e443
md"""
## Results
"""

# ╔═╡ ce0d320a-dec6-4734-8acb-6b800806381c
C = vecAdd(A, B, L)

# ╔═╡ cd32d0b6-b8ba-40a4-bdd9-0fbb052faa70
md"""
Verifying results...
"""

# ╔═╡ 83137e33-cbec-48a6-81b2-b60663b2d9c8
A .+ B

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"

[compat]
CUDA = "~3.6.4"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.1"
manifest_format = "2.0"

[[deps.AbstractFFTs]]
deps = ["ChainRulesCore", "LinearAlgebra"]
git-tree-sha1 = "6f1d9bc1c08f9f4a8fa92e3ea3cb50153a1b40d4"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.1.0"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.BFloat16s]]
deps = ["LinearAlgebra", "Printf", "Random", "Test"]
git-tree-sha1 = "a598ecb0d717092b5539dbbe890c98bac842b072"
uuid = "ab4f0b2a-ad5b-11e8-123f-65d77653426b"
version = "0.2.0"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.CEnum]]
git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.1"

[[deps.CUDA]]
deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CompilerSupportLibraries_jll", "ExprTools", "GPUArrays", "GPUCompiler", "LLVM", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "Requires", "SparseArrays", "SpecialFunctions", "TimerOutputs"]
git-tree-sha1 = "5bd05cd090588389201f7e1c703d50ace07be861"
uuid = "052768ef-5323-5732-b1bb-66c8b64840ba"
version = "3.6.4"

[[deps.ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "6e39c91fb4b84dcb870813c91674bdebb9145895"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.11.5"

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "44c37b4636bc54afac5c574d2d02b625349d6582"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.41.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "b19534d1895d702889b219c382a6e18010797f0b"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.6"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.ExprTools]]
git-tree-sha1 = "24565044e60bc48a7562e75bcf14f084901dc0b6"
uuid = "e2ba6199-217a-4e67-a87a-7c52f15ade04"
version = "0.1.7"

[[deps.GPUArrays]]
deps = ["Adapt", "LinearAlgebra", "Printf", "Random", "Serialization", "Statistics"]
git-tree-sha1 = "d9681e61fbce7dde48684b40bdb1a319c4083be7"
uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
version = "8.1.3"

[[deps.GPUCompiler]]
deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "TimerOutputs", "UUIDs"]
git-tree-sha1 = "2cac236070c2c4b36de54ae9146b55ee2c34ac7a"
uuid = "61eb1bfa-7361-4325-ad38-22787b887f55"
version = "0.13.10"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "22df5b96feef82434b07327e2d3c770a9b21e023"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.0"

[[deps.LLVM]]
deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Printf", "Unicode"]
git-tree-sha1 = "f8dcd7adfda0dddaf944e62476d823164cccc217"
uuid = "929cbde3-209d-540e-8aea-75f648917ca0"
version = "4.7.1"

[[deps.LLVMExtra_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "62115afed394c016c2d3096c5b85c407b48be96b"
uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab"
version = "0.0.13+1"

[[deps.LazyArtifacts]]
deps = ["Artifacts", "Pkg"]
uuid = "4af54fe1-eca0-43a8-85a7-787d91b784e3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "e5718a00af0ab9756305a0392832c8952c7426c1"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.6"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[deps.OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "2cf929d64681236a2e074ffafb8d568733d2e6af"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.3"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Random123]]
deps = ["Libdl", "Random", "RandomNumbers"]
git-tree-sha1 = "0e8b146557ad1c6deb1367655e052276690e71a3"
uuid = "74087812-796a-5b5d-8853-05524746bad3"
version = "1.4.2"

[[deps.RandomNumbers]]
deps = ["Random", "Requires"]
git-tree-sha1 = "043da614cc7e95c703498a491e2c21f58a2b8111"
uuid = "e6cf234a-135c-5ec9-84dd-332b85af5143"
version = "1.5.3"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "e08890d19787ec25029113e88c34ec20cac1c91e"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.0.0"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.TimerOutputs]]
deps = ["ExprTools", "Printf"]
git-tree-sha1 = "a5aed757f65c8a1c64503bc4035f704d24c749bf"
uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
version = "0.5.14"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═b8f90c8e-788d-11ec-0f86-c3a393796a2d
# ╟─40b2c8ff-8c11-47f7-9dd5-19efb9663b0e
# ╠═c762d0a0-cb54-47e4-9ad8-51eac1cd437c
# ╠═a15b15e1-cbf0-4843-98dd-3dc6d9fdd41e
# ╠═51c98125-2ca5-4570-819a-162a172a5390
# ╟─a0c6094a-dc63-433e-a122-44914c90d200
# ╟─071054b2-c3d4-4d0f-9f6c-84e6f62246d9
# ╠═077f54d1-1bde-47ea-9192-291c88109d77
# ╟─4e168101-9f5f-4a2b-89d1-6a5e46178382
# ╠═ee948bcb-66e3-41a1-9795-7d8e7b57d806
# ╟─622542f8-1417-4cca-8308-f2b6ea61e443
# ╠═ce0d320a-dec6-4734-8acb-6b800806381c
# ╟─cd32d0b6-b8ba-40a4-bdd9-0fbb052faa70
# ╠═83137e33-cbec-48a6-81b2-b60663b2d9c8
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# plots the suites ran with runsuites.jl and saves plots to output/

#TODO: Use ColorMaps to get arbritary number of well spaced colors
cols = split("ymcrgbk","")
push!(cols, )
nc = length(cols)

function make_legend!(leg, res)
    str = string(name(res[1].testrun.solver))
    str = replace(str, "_", "\\_")
    push!(leg, str)
end

## Adaptive steppers
# significant digits (scd) vs walltime
for (n,tc) in totest
    leg = AbstractString[]
    id = Py.figure(figsize=(50,50),dpi=130)
    colind = 1
    p = 1
    maxscd = 0.0

    # DASSL.jl
    if isdefined(:resDASSL)
        res = resDASSL[n]
        scd = getfield_vec(res, :scd)
        maxscd = max(maxscd, maximum(scd))
        wt = getfield_vec(res, :walltime)
        p2 = Py.semilogy(scd, wt, "-x"*cols[colind])
        make_legend!(leg, res)
        colind +=1
    end

    # ODE.jl
    if isdefined(:resODE)
        rode = resODE[n]
        for (s,res) in rode
            scd = getfield_vec(res, :scd)
            if all(isnan(scd))
                continue
            end
            maxscd = max(maxscd, maximum(scd))
            wt = getfield_vec(res, :walltime)
            p2 = Py.plot(scd, wt, "-o"*cols[rem1(colind,nc)])
            make_legend!(leg, res)
            colind +=1
        end
    end

    #Sundials
    if isdefined(:resSun)
        rsun = resSun[n]
        for (s,res) in rsun
            scd = getfield_vec(res, :scd)
            if all(isnan(scd))
                continue
            end
            maxscd = max(maxscd, maximum(scd))
            wt = getfield_vec(res, :walltime)
            @show cols[rem1(colind,nc)]
            p = Py.plot(scd, wt, "-d"*cols[rem1(colind,nc)])
            make_legend!(leg, res)
            colind +=1
        end
    end

    if colind==1
        Py.close(id) # no plots were made
    else
        Py.legend(leg, loc="upper left")
        Py.title("$n")
        Py.xlabel("significant digits")
        Py.ylabel("Walltime (s)")
        Py.display(id)
        Py.savefig(Pkg.dir()*"/IVPTestSuite/testsuites/output/adaptive-scd-vs-walltime-$n.png")
        #Py.close(id)
    end
end

## Fixed step solvers
# significant digits (scd) vs walltime
if isdefined(:resODEfixed)
    for (n,tc) in totest

        leg = AbstractString[]
        id = Py.figure(figsize=(50,50),dpi=130)
        colind = 1
        p = 1
        p2 = 1
        maxscd = 0.0

        # # ODE.jl
        rode = resODEfixed[n]
        if length(rode)==0
            Py.close(id)
            continue
        end
        fst = true
        for (s,res) in rode
            scd = getfield_vec(res, :scd)
            if all(isnan(scd))
                continue
            end
            maxscd = max(maxscd, maximum(scd))
            wt = getfield_vec(res, :walltime)
            if fst
                p2 = Py.semilogy(scd, wt, "-o"*cols[colind])
                fst = false
            else
                p2 = Py.plot(scd, wt, "-o"*cols[colind],hold = true)

            end
            make_legend!(leg, res)
            colind +=1
        end

        Py.legend(leg)
        Py.title("$n (fixed step)")
        Py.xlabel("significant digits")
        Py.ylabel("Walltime (s)")

        Py.display(id)
        Py.savefig(Pkg.dir()*"/IVPTestSuite/testsuites/output/fixedstep-scd-vs-walltime-$n.png")
        #Py.close(id)
    end
end


text:
module Operators

using Oceananigans
using Oceananigans.Grids
using Oceananigans.Operators

using Oceananigans: AbstractGrid
using Oceananigans.TurbulenceClosures: IsotropicDiffusivity

export
    kinetic_energy,
    hdivᶜᶜᵃ, divᶜᶜᶜ, ∇²,
    div_ρuũ, div_ρvũ, div_ρwũ, div_ρUc,
    ∂ⱼpuⱼ, ∂ⱼDᶜⱼ, ∂ⱼtᶜDᶜⱼ, ∂ⱼDᵖⱼ,
    ∂ⱼτ₁ⱼ, ∂ⱼτ₂ⱼ, ∂ⱼτ₃ⱼ, Q_dissipation

include("compressible_operators.jl")

end


text:
using MDTable
using Test

@testset "Write Tests" begin
    let loaded = readMDTable("file.md")
        writeMDTable("fileTest.md", loaded)
        checking = readMDTable("fileTest.md")

        @test loaded == checking
    end

    let err = nothing
        try
            writeMDTable("fileTest.md", ("a", "b", "c"))
        catch err
        end

        @test err isa Exception
        @test sprint(showerror, err) == "ArgumentError: 'Tuple{String,String,String}' iterates 'String' values, which doesn't satisfy the Tables.jl `AbstractRow` interface"
    end
end


text:
"""
Use ECOS to solve the bound-constrained least-squares problem

    minimize ||Ax-b||_2 subj to ℓ ≤ x ≤ u.

ECOS is a primal-dual interior method for SOCP. This module translates
the problem into the SOCP

    minimize  t  subj to Ax + r = b, (||x||,t) ∈ SOC, ℓ ≤ x ≤ u.

See [`lsbox`](@ref)
"""
module ECOSProblems

using ECOS
using LinearAlgebra
using SparseArrays

# export ls_setup, ls_simplex, solve!, cleanup!
export lsbox
# export ECOSls

include("lsbox.jl")

# abstract type ECOSWrapper end

# solve!(p::ECOSWrapper) = ECOS.solve(p.ecos_ptr)
# cleanup!(p::ECOSWrapper) = ECOS.cleanup(p.ecos_ptr)

end # module


text:
### A Pluto.jl notebook ###
# v0.19.3

using Markdown
using InteractiveUtils

# ╔═╡ 89fd640c-c649-4a06-8ce8-b0e40998ef9b
using PlutoLinks: @ingredients

# ╔═╡ effcea86-c002-11ec-0226-6d894daa7812
md"""
# `@ingredients`

```julia
MyModule = @ingredients(filename::AbstractString)
```

The `@ingredients` macro is a macro to use standalone Julia files inside Pluto notebooks much like `include()` can be used in the Julia REPL. The returned value is a module containing the functions defined in the file.

The macro also watches the file for changes and triggers a reactive update from its cell whenever the file is edited. For the non-reactive variant, check out `PlutoLinks.ingredients()`.

The macro can also be used to import functions from a Pluto notebook inside another notebook with currently a few constraints (environment are not synced, the whole file is recompiled at each update). For the whole package version of `@ingredients()` consider using [`@PlutoLinks.revise()`](./revise) which watches an entire package.

## Examples

Let's first create a very simple Julia file called `my_file.jl`:

```julia
# my_file.jl
function greet(name)
	return "Hello, $(name)!"
end
```

This file can be then be imported inside the current notebook by using the ingredients macro:
"""

# ╔═╡ 8b06607c-ffb2-47a9-8398-32608cd0a832
MyModule = @ingredients("./my_file.jl")

# ╔═╡ 3aeaf1b4-5850-4914-8adf-ea001f63d2af
md"""
`MyModule` can then be used to access the different functions defined in the file:
"""

# ╔═╡ a4b2c92b-25fe-4137-a8ad-388daa096d20
MyModule.greet("World")

# ╔═╡ 7e6d4b1b-d89e-47f1-9450-f33ece70c17a
md"""
To use the variables without explicitely using `MyModule`, an alias can be specified to move the functions or variables in the main namespace. For example, here `greet` is imported inside the notebook module and can be used directly:
"""

# ╔═╡ 6b5265ba-489f-496b-8bb7-a776987dbed6
greet = MyModule.greet

# ╔═╡ 707c7498-dba7-4d20-8286-5778104df208
greet("Baz")

# ╔═╡ bdd2f924-402c-4858-bebb-72c5c20696c0
md"""
---
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoLinks = "0ff47ea0-7a50-410d-8455-4348d5de0420"

[compat]
PlutoLinks = "~0.1.5"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "6d4fa04343a7fc9f9cb9cff9558929f3d2752717"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.9"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "52617c41d2761cc05ed81fe779804d3b7f14fff7"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.9.13"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "dedbebe234e06e1ddad435f5c6f4b85cd8ce55f7"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.2.2"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[deps.PlutoHooks]]
deps = ["InteractiveUtils", "Markdown", "UUIDs"]
git-tree-sha1 = "072cdf20c9b0507fdd977d7d246d90030609674b"
uuid = "0ff47ea0-7a50-410d-8455-4348d5de0774"
version = "0.0.5"

[[deps.PlutoLinks]]
deps = ["FileWatching", "InteractiveUtils", "Markdown", "PlutoHooks", "Revise", "UUIDs"]
git-tree-sha1 = "0e8bcc235ec8367a8e9648d48325ff00e4b0a545"
uuid = "0ff47ea0-7a50-410d-8455-4348d5de0420"
version = "0.1.5"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

[[deps.Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
git-tree-sha1 = "4d4239e93531ac3e7ca7e339f15978d0b5149d03"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "3.3.3"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─effcea86-c002-11ec-0226-6d894daa7812
# ╠═8b06607c-ffb2-47a9-8398-32608cd0a832
# ╟─3aeaf1b4-5850-4914-8adf-ea001f63d2af
# ╠═a4b2c92b-25fe-4137-a8ad-388daa096d20
# ╟─7e6d4b1b-d89e-47f1-9450-f33ece70c17a
# ╠═6b5265ba-489f-496b-8bb7-a776987dbed6
# ╠═707c7498-dba7-4d20-8286-5778104df208
# ╟─bdd2f924-402c-4858-bebb-72c5c20696c0
# ╠═89fd640c-c649-4a06-8ce8-b0e40998ef9b
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""Tests to show that the benchmarks we based our speed optimizations on are
still valid"""

# reload("Util.jl")
include("../src/Util.jl")
# using Benchmark
using Base.Test

# TODO: Implement convergence in Benchmark necessetating some LiveStates functionality too.

macro calltimes(ex, times)
    quote
        for i in 1:$times
            $ex
        end
    end
end

function setup_lookup(n, m)
    dict = Dict{(Int, Int), Float64}()
    for i in 1:n
        for j in 1:m
            dict[(i, j)] = 0.0
        end
    end
    return dict
end

function do_baseline(n::Int, m::Int)
    for i in 1:n
        for j in 1:m
        end
    end
end

function bench_lookup(n::Int, m::Int, dict)
    for i in 1:n
        for j in 1:m
            dict[(i, j)]
        end
    end
end

function setup_lookup(n)
    dict = Dict{Int, Float64}()
    for i in 1:n
        dict[i] = 0.0
    end
    return dict
end

function do_baseline(n::Int)
    for i in 1:n
    end
end

function bench_lookup(n::Int, dict)
    for i in 1:n
        dict[i]
    end
end

dict = setup_lookup(40000)
baseline = mean([(@elapsed (@calltimes do_baseline(40000) 3)) for i in 1:5]) / 3.
time_lookup = mean([(@elapsed (@calltimes bench_lookup(40000, dict) 3)) for i in 1:5]) / 3.

@printf "Julia int-key lookup: %.2f milliseconds\n" (time_lookup - baseline) * 1000.

dict = setup_lookup(200, 200)
baseline = mean([(@elapsed (@calltimes do_baseline(200, 200) 3)) for i in 1:5]) / 3.
time_lookup = mean([(@elapsed (@calltimes bench_lookup(200, 200, dict) 3)) for i in 1:5]) / 3.

@printf "Julia tuple-key lookup: %.2f milliseconds\n" (time_lookup - baseline) * 1000.

tcount = 0
type Trex
    depth::Int
    child

    function Trex(depth)
        global tcount
        tcount += 1

        if(depth > 0)
            new(depth, (Trex(depth - 1), Trex(depth - 1)))
        else
            new(depth, nothing)
        end
    end
end

function reset_tcount()
    global tcount
    tcount = 0
end
reset_tcount()

@time Trex(10)
@show tcount


text:
# see #1305
@test isfinite(mean(Chi(1000)))

# see #1497
for x in (-0.5, -1, -1.5f0), ν in (3, 3f0, 3.0)
    @test @inferred(logpdf(Chi(ν), x)) == -Inf
end


text:
function __init__ogbdataset()
    DEPNAME = "OGBDataset"
    LINK = "http://snap.stanford.edu/ogb/data"
    DOCS = "https://ogb.stanford.edu/docs/dataset_overview/"

    # NOT WORKING AS EXPECTED
    function fetch_method(remote_filepath, local_directory_path)
        if endswith(remote_filepath, "nodeproppred/master.csv")
            local_filepath = joinpath(local_directory_path, "nodeproppred_metadata.csv")
            download(remote_filepath, local_filepath)
        elseif endswith(remote_filepath, "linkproppred/master.csv")
            local_filepath = joinpath(local_directory_path, "linkproppred_metadata.csv")
            download(remote_filepath, local_filepath)
        elseif endswith(remote_filepath, "graphproppred/master.csv")
            local_filepath = joinpath(local_directory_path, "graphproppred_metadata.csv")
            download(remote_filepath, local_filepath)
        else
            local_filepath = DataDeps.fetch_default(remote_filepath, local_directory_path)
        end
        return local_filepath
    end

    register(DataDep(
        DEPNAME,
        """
        Dataset: The $DEPNAME dataset.
        Website: $DOCS
        Download Link: $LINK
        """,
        ["https://raw.githubusercontent.com/snap-stanford/ogb/master/ogb/nodeproppred/master.csv",
        "https://raw.githubusercontent.com/snap-stanford/ogb/master/ogb/linkproppred/master.csv",
        "https://raw.githubusercontent.com/snap-stanford/ogb/master/ogb/graphproppred/master.csv",
        ],
        fetch_method=fetch_method
    ))
end

"""
    OGBDataset(name; dir=nothing)

The collection of datasets from the [Open Graph Benchmark: Datasets for Machine Learning on Graphs](https://arxiv.org/abs/2005.00687)
paper.

`name` is the name  of one of the dasets (listed [here](https://ogb.stanford.edu/docs/dataset_overview/))
available for node prediction, edge prediction, or graph prediction tasks.

The `OGBDataset` type stores the graphs internally as dictionary objects.
The key "edge_index" contains `2 x num_edges`, where the first and second
column contain the source and target nodes of each edge respectively.

# Examples

## Node prediction tasks

```julia-repl
julia> data = OGBDataset("ogbn-arxiv")
OGBDataset{Vector{Any}}:
  name => ogbn-arxiv
  path => /home/carlo/.julia/datadeps/OGBDataset/arxiv
  metadata => Dict{String, Any} with 15 entries
  graphs => 1-element Vector{Dict}
  labels => 1-element Vector{Any}
  split => Dict{String, Any} with 3 entries

julia> data.metadata
Dict{String, Any} with 15 entries:
  "num classes"           => 40
  "binary"                => false
  "is hetero"             => false
  "eval metric"           => "acc"
  "task type"             => "multiclass classification"
  "version"               => 1
  "split"                 => "time"
  "download_name"         => "arxiv"
  "num tasks"             => 1
  "url"                   => "http://snap.stanford.edu/ogb/data/nodeproppred/arxiv.zip"
  "additional node files" => "node_year"
  "add_inverse_edge"      => false
  "has_node_attr"         => true
  "additional edge files" => nothing
  "has_edge_attr"         => false

julia> data.split
Dict{String, Any} with 3 entries:
  "test_idx"  => [347, 399, 452, 481, 489, 491, 527, 538, 541, 603  …  169334, 169335, 169336, 169337, 169338, 169339, 169340, 169341, 169342, 169343]
  "train_idx" => [1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  169110, 169112, 169113, 169114, 169115, 169116, 169118, 169146, 169149, 169252]
  "val_idx"   => [350, 358, 367, 383, 394, 422, 430, 436, 468, 470  …  169089, 169096, 169108, 169111, 169128, 169156, 169177, 169186, 169262, 169297]

julia> length(data)
1

julia> graph, labels = data[1];

julia> graph
Dict{String, Any} with 6 entries:
  "edge_index" => [104448 13092; 15859 47284; … ; 45119 162538; 45119 72718]
  "edge_feat"  => nothing
  "node_feat"  => Float32[-0.057943 -0.1245 … -0.138236 -0.029875; -0.05253 -0.070665 … 0.040885 0.268417; … ; -0.172796 -0.372111 … -0.041253 0.077647; -0.140059 -0.301036 … -0.376132 -0.091018]
  "num_nodes"  => 169343
  "node_year"  => [2013 2015 … 2020 2020]
  "num_edges"  => 1166243

julia> source, target = graph["edge_index][:,1], graph["edge_index][:,2];
```

## Edge prediction task

```julia-repl
julia> data = OGBDataset("ogbl-collab")
OGBDataset{Nothing}:
  name => ogbl-collab
  path => /home/carlo/.julia/datadeps/OGBDataset/collab
  metadata => Dict{String, Any} with 13 entries
  graphs => 1-element Vector{Dict}
  labels => nothing
  split => Dict{String, Any} with 3 entries

julia> graph = data[1]  # no labels for this dataset
Dict{String, Any} with 7 entries:
  "edge_index"  => [150990 224882; 150990 224882; … ; 221742 135759; 207233 140615]
  "edge_feat"   => nothing
  "node_feat"   => Float32[-0.177486 -0.237488 … 0.004236 -0.035025; -0.10298 0.022193 … 0.031942 -0.118059; … ; 0.003879 0.062124 … 0.05208 -0.176961; -0.276317 -0.081464 … -0.201557 -0.258715]
  "num_nodes"   => 235868
  "edge_year"   => [2004 2002 … 2006 1984; 2004 2002 … 2006 1984]
  "edge_weight" => [2 1 … 1 1; 2 1 … 1 1]
  "num_edges"   => 2358104
```

## Graph prediction task

```julia-repl
julia> data = OGBDataset("ogbg-molhiv")
OGBDataset{Matrix{Int64}}:
  name => ogbg-molhiv
  path => /home/carlo/.julia/datadeps/OGBDataset/molhiv
  metadata => Dict{String, Any} with 15 entries
  graphs => 41127-element Vector{Dict}
  labels => 1×41127 Matrix{Int64}
  split => Dict{String, Any} with 3 entries

julia> length(data)
41127

julia> graph, labels = data[10]
(Dict{String, Any}("edge_index" => [-202 -201; -201 -200; … ; -198 -184; -201 -202], "node_feat" => Float32[7.0 6.0 … 7.0 7.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0], "edge_feat" => Float32[0.0 0.0 … 0.0 1.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 1.0], "num_nodes" => 20, "num_edges" => 42), [0])

julia> graph, labels = data[10];

julia> graph
Dict{String, Any} with 5 entries:
  "edge_index" => [1 2; 2 3; … ; 5 19; 2 1]
  "edge_feat"  => Float32[0.0 0.0 … 0.0 1.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 1.0]
  "node_feat"  => Float32[7.0 6.0 … 7.0 7.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  "num_nodes"  => 20
  "num_edges"  => 42

julia> labels
1-element Vector{Int64}:
 0
```
"""
struct OGBDataset{L} <: AbstractDataset
    name::String
    metadata::Dict{String, Any}
    graphs::Vector{Graph}
    targets::L
    split_idx::NamedTuple
end

function OGBDataset(fullname; dir = nothing)
    metadata = read_ogb_metadata(fullname, dir)
    path = makedir_ogb(fullname, metadata["url"], dir)
    metadata["path"] = path
    graph_dicts, labels, split_idx = read_ogb_graph(path, metadata)
    graphs = ogbdict2graph.(graph_dicts)
    return OGBDataset(fullname, metadata, graphs, labels, split_idx)
end

function read_ogb_metadata(fullname, dir = nothing)
    dir = isnothing(dir) ? datadep"OGBDataset" : dir
    @assert contains(fullname, "-") "The full name should be provided, e.g. ogbn-arxiv"
    prefix, name = split(fullname, "-")
    if prefix == "ogbn"
        path_metadata = joinpath(dir, "nodeproppred_metadata.csv")
    elseif prefix == "ogbl"
        path_metadata = joinpath(dir, "linkproppred_metadata.csv")
    elseif prefix == "ogbg"
        path_metadata = joinpath(dir, "graphproppred_metadata.csv")
    else
        @assert "The dataset name should start with ogbn, ogbl, or ogbg."
    end
    df = read_csv(path_metadata)
    @assert fullname ∈ names(df)
    metadata = Dict{String, Any}(String(r[1]) => parse_pystring(r[2]) for r in eachrow(df[!,[names(df)[1], fullname]]))
    return metadata
end

function makedir_ogb(fullname, url, dir = nothing)
    root_dir = isnothing(dir) ? datadep"OGBDataset" : dir
    @assert contains(fullname, "-") "The full name should be provided, e.g. ogbn-arxiv"
    prefix, name = split(fullname, "-")
    data_dir = joinpath(root_dir, name)
    if !isdir(data_dir)
        local_filepath = DataDeps.fetch_default(url, root_dir)
        currdir = pwd()
        cd(root_dir) # Needed since `unpack` extracts in working dir
        DataDeps.unpack(local_filepath)
        unzipped = local_filepath[1:findlast('.', local_filepath)-1]
        mv(unzipped, data_dir)
        for (root, dirs, files) in walkdir(data_dir)
            for file in files
                if endswith(file, r"zip|gz")
                    cd(root)
                    DataDeps.unpack(joinpath(root, file))
                end
            end
        end
        cd(currdir)
    end
    @assert isdir(data_dir)

    return data_dir
end

## See https://github.com/snap-stanford/ogb/blob/master/ogb/io/read_graph_raw.py
function read_ogb_graph(path, metadata)
    dict = Dict{String, Any}()

    dict["edge_index"] = read_ogb_file(joinpath(path, "raw", "edge.csv"), Int, transp=false)
    dict["edge_index"] = dict["edge_index"] .+ 1 # from 0-indexing to 1-indexing

    dict["node_features"] = read_ogb_file(joinpath(path, "raw", "node-feat.csv"), Float32)
    dict["edge_features"] = read_ogb_file(joinpath(path, "raw", "edge-feat.csv"), Float32)

    dict["num_nodes"] = read_ogb_file(joinpath(path, "raw", "num-node-list.csv"), Int; tovec=true)
    dict["num_edges"] = read_ogb_file(joinpath(path, "raw", "num-edge-list.csv"), Int; tovec=true)

    for file in readdir(joinpath(path, "raw"))
        if file ∉ ["edge.csv", "num-node-list.csv", "num-edge-list.csv", "node-feat.csv", "edge-feat.csv"]
            propname = replace(split(file,".")[1], "-" => "_")
            dict[propname] = read_ogb_file(joinpath(path, "raw", file), Any)

            # from https://github.com/snap-stanford/ogb/blob/5e9d78e80ffd88787d9a1fdfdf4079f42d171565/ogb/io/read_graph_raw.py
            # # hack for ogbn-proteins
            # if additional_file == 'node_species' and osp.exists(osp.join(raw_dir, 'species.csv.gz')):
            #     os.rename(osp.join(raw_dir, 'species.csv.gz'), osp.join(raw_dir, 'node_species.csv.gz'))
        end
    end

    node_keys = [k for k in keys(dict) if startswith(k, "node_")]
    edge_keys = [k for k in keys(dict) if startswith(k, "edge_") && k != "edge_index"]
    # graph_keys = [k for k in keys(dict) if startswith(k, "graph_")] # no graph-level features in official implementation

    num_graphs = length(dict["num_nodes"])
    @assert num_graphs == length(dict["num_edges"])

    graphs = Dict[]
    num_node_accum = 0
    num_edge_accum = 0
    for i in 1:num_graphs
        graph = Dict{String, Any}()
        n, m = dict["num_nodes"][i], dict["num_edges"][i]
        graph["num_nodes"] = n
        graph["num_edges"] = metadata["add_inverse_edge"] ? 2*m : m

        # EDGE FEATURES
        if metadata["add_inverse_edge"]
            ei = dict["edge_index"][num_edge_accum+1:num_edge_accum+m, :]
            s, t = ei[:,1], ei[:,2]
            graph["edge_index"] = [s t; t s]

            for k in edge_keys
                v = dict[k]
                if v === nothing
                    graph[k] = nothing
                else
                    x = v[:, num_edge_accum+1:num_edge_accum+m]
                    graph[k] = [x; x]
                end
            end
        else
            graph["edge_index"] = dict["edge_index"][num_edge_accum+1:num_edge_accum+m, :]

            for k in edge_keys
                v = dict[k]
                if v === nothing
                    graph[k] = nothing
                else
                    graph[k] = v[:, num_edge_accum+1:num_edge_accum+m]
                end
            end
            num_edge_accum += m
        end

        # NODE FEATURES
        for k in node_keys
            v = dict[k]
            if v === nothing
                graph[k] = nothing
            else
                graph[k] = v[:, num_node_accum+1:num_node_accum+n]
            end
        end
        num_node_accum += n

        push!(graphs, graph)
    end

    # PROCESS LABELS
    dlabels = Dict{String, Any}()
    for k in keys(dict)
        if contains(k, "label")
            if k ∉ [node_keys; edge_keys]
                dlabels[k] = dict[k]
            end
        end
    end
    labels = isempty(dlabels) ? nothing :
             length(dlabels) == 1 ? first(dlabels)[2] : dlabels

    splits = readdir(joinpath(path, "split"))
    @assert length(splits) == 1 # TODO check if datasets with multiple splits existin in OGB
    # TODO sometimes splits are given in .pt format
    split_idx = (train = read_ogb_file(joinpath(path, "split", splits[1], "train.csv"), Int; tovec=true),
                 val = read_ogb_file(joinpath(path, "split", splits[1], "valid.csv"), Int; tovec=true),
                test = read_ogb_file(joinpath(path, "split", splits[1], "test.csv"), Int; tovec=true))
    if split_idx.train !== nothing
        split_idx.train .+= 1
    end
    if split_idx.val !== nothing
        split_idx.val .+= 1
    end
    if split_idx.test !== nothing
        split_idx.test .+= 1
    end
    return graphs, labels, split_idx
end

function read_ogb_file(p, T; tovec = false, transp = true)
    res = isfile(p) ? read_csv(p, Matrix{T}, header=false) : nothing
    if tovec && res !== nothing
        @assert size(res, 1) == 1 || size(res, 2) == 1
        res = vec(res)
    end
    if transp && res !== nothing && !tovec
        res = collect(res')
    end
    if res !== nothing T === Any
        res = restrict_array_type(res)
    end
    return res
end

function ogbdict2graph(d::Dict)
    edge_index = d["edge_index"][:,1], d["edge_index"][:,2]
    num_nodes, num_edges = d["num_nodes"], d["num_edges"]
    node_data = Dict(Symbol(k[6:end]) => v for (k,v) in d if startswith(k, "node_") && v !== nothing)
    edge_data = Dict(Symbol(k[6:end]) => v for (k,v) in d if startswith(k, "edge_") && k!="edge_index" && v !== nothing)
    node_data = isempty(node_data) ? nothing : (; node_data...)
    edge_data = isempty(edge_data) ? nothing : (; edge_data...)
    return Graph(; num_nodes, num_edges,
                 edge_index, node_data, edge_data)
end

Base.length(data::OGBDataset) = length(data.graphs)
Base.getindex(data::OGBDataset{Nothing}, i) = getobs(data.graphs, i)
Base.getindex(data::OGBDataset, i) = getobs((; data.graphs, data.targets), i)


text:
using ConstrainedDynamics
using ConstrainedControl
using LinearAlgebra

# Parameters
joint_axis = [1.0;0.0;0.0]

length1 = 1.0
width, depth = 0.1, 0.1

# Links
origin = Origin{Float64}()
link1 = Box(width, depth, width, length1)

# Constraints
joint_between_origin_and_link1 = EqualityConstraint(Prismatic(origin, link1, joint_axis))

links = [link1]
constraints = [joint_between_origin_and_link1]

mech = Mechanism(origin, links, constraints, g=0.)
setPosition!(origin,link1,Δx = [1.0;0;0])

Q = ones(1)
R = ones(1)

lqr = LQR(mech, getid.(constraints), getid.(constraints), Q, R, 10.)
@test true


text:
# Intercepts for `mapreduce`, `mapfoldl` and `mapfoldr` under `op` `+`.
const plustype = Union{typeof(+), typeof(Base.add_sum)}
const type_tuple = :(Tuple{Any, $plustype, ∇ArrayOrScalar})
for f in (:mapfoldl, :mapfoldr)
    @eval begin
        import Base: $f
        @explicit_intercepts $f $type_tuple [false, false, true] #(init=0,)
        function ∇(::typeof($f), ::Type{Arg{3}}, p, y, ȳ, f, ::$plustype, A::∇ArrayOrScalar)
            return ȳ .* ForwardDiff.derivative.(f, A)
        end
    end
end


text:
# Version file

version() = v"0.1.0"


text:
##
# sample that plots usage trend of any hashtag, url, or smiley string over time.
#
# 1. bring up julia on all hadoop data nodes
#       julia --machinefile ${HADOOP_HOME}/conf/slaves
# 2. load this file (assuming file present in cwd)
#       julia> require("twitter_test_smileys.jl")
# 3. we would use the daily twitter summary data from infochimps
#    data source: http://www.infochimps.com/datasets/twitter-census-conversation-metrics-one-year-of-urls-hashtags-sm--2
#    ensure that the daily summary file is available at /twitter_daily_summary.tsv
#
# 4. run the test function
#    julia> do_plot_counts("hdfs://host:port/twitter_daily_summary.tsv")
#
# TRENDS OBSERVED:
#   - flu/swineflu/birdflu: Apr 2009
#   - haloween: Oct
#   - thanksgiving: Oct-Nov
#   - obama: Nov 2008 (elected), Jan 2009 (speech), Oct 2009 (peace prize)
#   - earthquake: Jan 2009 (costarica), Mar 2009 (Tonga), Apr 2009 (italy), Sep 2009 (southern sumatra)
#
# NOTE:
#   to show plots you need to load one of the following set of modules:
#       - Gaston
#       - Winston and Tk
#
#   If you are using Winston, it may help to modify the following in Winston.ini:
#   ticklabels_offset       = 5
#   ticklabels_style        = {fontsize:1.0, textangle:90}
#

using DataFrames
using HDFS
using HDFS.MapReduce

##
# find smiley records from HDFS CSV file
find_recs_as_df(r, next_rec_pos) = HDFS.MapReduce.find_rec(r, next_rec_pos, DataFrame, '\n', '\t')

function midx(yyyymmdd)
    yyyymmdd = int(yyyymmdd)
    ypos = (yyyymmdd > 10^9) ? 10^6 : 10^4
    y = int(floor(yyyymmdd/ypos))
    m = int(floor((yyyymmdd - ypos*y)/(ypos/100)))
    12*(y-2006)+m
end
function map_count_monthly(rec, tag_type::String, tag::Regex)
    ((nothing == rec) || (nrow(rec) == 0)) && return []
    #println(rec)

    filtrows = Int[]
    for (i,v) in enumerate(rec[4])
        (tag_type == rec[i,1]) && ismatch(tag,v) && push!(filtrows,i)
    end
    subrec = subset(rec, filtrows)
    (nrow(subrec) == 0) && return []

    subrec = DataFrame({subrec[1], PooledDataArray(map(x->midx(x), subrec[2])), subrec[3], subrec[4]})
    monthsum = by(subrec, :x2, x -> sum(x[:x3]))

    [(tag.pattern, monthsum[idx,1], monthsum[idx,2]) for idx in 1:nrow(monthsum)]
end

function collect_count_monthly(results, rec)
    (length(rec) == 0) && return results
    hashtag, month_idx, cnt = rec

    local monthly::Vector{Int}
    try
        monthly = getindex(results, hashtag)
    catch
        (nothing == results) && (results = Dict{String, Vector{Int}}())
        # represent 5 years worth of monthly data
        monthly = zeros(Int, 12*5)
        results[hashtag] = monthly
    end

    monthly[month_idx] += cnt
    results
end

reduce_count_monthly(reduced, results...) = HDFS.MapReduce.reduce_dicts(+, reduced, results...)

##
# twitter archive search begin
const colors = ["blue", "red", "green", "violet", "black", "yellow", "magenta", "brown", "cyan", "pink"]
const months = 6:46
const monnames = ["jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec"]
const yearmonths = String[]
for y in [2006,2007,2008,2009]
    for m in monnames
        push!(yearmonths, string("$(m) $(y)"))
    end
end

function do_plot_counts(furl::String, typ::String, tag::String)
    function is_loaded(x::Symbol)
        try
            return (Module == typeof(eval(x)))
        catch
            return false
        end
    end

    function plot_winston()
        idx = 0
        win = Toplevel(string(tag, " in ", typ, "s"), 900, 600)
        c = Canvas(win, 900, 600)
        pack(c, expand=true, fill="both")

        p = FramedPlot()
        #setattr(p.x1, "ticklabels", yearmonths[6:5:46])
        #setattr(p.x1, "ticks", 6:5:46)
        setattr(p.x1, "ticklabels", yearmonths[6:46])
        setattr(p.x1, "ticks", 6:46)

        for (key, val) in counts
            idx += 1
            add(p, Curve(months, val[6:46], "color", colors[idx]))
            (idx == 10) && (idx = 0)
        end
        Winston.display(c,p)
    end

    function plot_gaston()
        set_terminal("x11")
        params = {}
        idx = 0
        for (key, val) in counts
            idx += 1
            push!(params, months); push!(params, val[6:46])
            push!(params, "color"); push!(params, colors[idx])
            push!(params, "legend"); push!(params, key)
            (idx == 10) && (idx = 0)
        end
        push!(params, "title"); push!(params,string(tag, " in ", typ, "s"))
        push!(params, "xlabel"); push!(params, "months")
        push!(params, "ylabel"); push!(params, "counts")
        plot(params...)
    end

    function wait_results()
        loopstatus = true
        while(loopstatus)
            sleep(5)
            jstatus,jstatusinfo = status(j_mon,true)
            ((jstatus == "error") || (jstatus == "complete")) && (loopstatus = false)
            (jstatus == "running") && println("$(j_mon): $(jstatusinfo)% complete...")
        end
        wait(j_mon)
        println("time taken (total time, wait time, run time): $(times(j_mon))")
        println("")
    end

    function print_results()
        println(counts)
        ((nothing == counts) || (0 == length(counts))) && return

        println("maximum:")
        for (key,val) in counts
            maxmon = indmax(val)
            maxyear = 2006 + int(floor(maxmon/12))
            maxmonname = ((maxmon-1)%12) + 1
            println("$(key): $(monnames[maxmonname]) $(maxyear)")
        end

        println("local maximas:")
        for (key,val) in counts
            len = length(val)
            avgs = [mean(val[max(1,(n-2)):min(len,(n+2))]) for n in 1:len]
            for n in 1:len
                if((val[n] > 25) && (val[n] > (1.25)*avgs[n]))
                    #println("val $(val[n]), avgs $(avgs[n]), n $(n)")
                    maxyear = 2006 + int(floor((n-1)/12))
                    maxmonname = ((n-1)%12) + 1
                    println("$(key): $(monnames[maxmonname]) $(maxyear)")
                end
            end
        end

        counts
    end

    ##
    # function body begin
    println("starting dmapreduce...")
    j_mon = dmapreduce(MRHdfsFileInput([furl], (x,y)->find_recs_as_df(x,y)), x->map_count_monthly(x, "smiley", Regex(tag)), collect_count_monthly, reduce_count_monthly)

    println("waiting for dmapreduce to finish...")
    wait_results()

    println("results:")
    counts = results(j_mon)[2]
    print_results()

    println("plot coming up...")
    is_loaded(:Winston) ? plot_winston() :
    is_loaded(:Gaston) ? plot_gaston() :
    println("can not plot results. no graphics library loaded.")
end


text:
module BuildMachines

# TODO: delete this function
_hello_world() = "Hello world!"

end


text:
using DiffEqBase, Random, LinearAlgebra, Test
using DiffEqBase: isconstant

@testset "Identity Operators" begin
  u = [1.0, 2.0]; du = zeros(2)
  Id = DiffEqIdentity(u)
  @test Id * u == u
  mul!(du, Id, u); @test du == u
  @test size(Id) == (2,2)
end

@testset "Scalar Operators" begin
  u = [1.0, 2.0]; u2 = [1.0, 2.0]
  α = DiffEqScalar(2.0)
  @test convert(Number, α) == 2.0
  @test α * u == 2.0u
  lmul!(α, u2); @test u2 == 2.0u
  @test size(α) == ()
  @test isconstant(α) == true
end

@testset "Array Operators" begin
  Random.seed!(0); A = rand(2,2); u = rand(2); du = zeros(2)
  L = DiffEqArrayOperator(A)
  @test Matrix(L) == A
  @test size(L) == size(A)
  @test L * u == A * u
  mul!(du, L, u); @test du == A*u
  @test lu(L) \ u ≈ A \ u
  @test opnorm(L) == opnorm(A)
  @test exp(L) == exp(A)
  @test L[1,2] == A[1,2]
  @test isconstant(L) == true
end

@testset "Mutable Array Operators" begin
  Random.seed!(0); A = rand(2,2); u = rand(2); du = zeros(2)
  update_func = (_A,u,p,t) -> _A .= t * A
  Lt = DiffEqArrayOperator(zeros(2,2); update_func=update_func)
  t = 5.0
  @test isconstant(Lt) == false
  @test Lt(u,nothing,t) ≈ (t*A) * u
  Lt(du,u,nothing,t); @test du ≈ (t*A) * u
end


text:
using Base.Test
using SymEngine

x = symbols("x")
y = symbols(:y)
@vars z

# Check Basic conversions
@test eltype([Basic(u) for u in [1, 1/2, 1//2, pi, e]]) == Basic

# make sure @vars defines in a local scope
let
    @vars w
end
@test_throws UndefVarError isdefined(w)

a = x^2 + x/2 - x*y*5
b = diff(a, x)
@test b == 2*x + 1//2 - 5*y

c = x + Rational(1, 5)
c = expand(c * 5)
@test c == 5*x + 1

c = x ^ 5
@test diff(c, x) == 5 * x ^ 4

c = x ^ y
@test c != y^x

c = Basic(-5)
@test abs(c) == 5
@test abs(c) != 4

show(a)
println()
show(b)
println()

## mathfuns
@test abs(Basic(-1)) == 1
@test sin(Basic(1)) == subs(sin(x), x, 1)
@test sin(PI) == 0
@test subs(sin(x), x, pi) == 0
@test sind(Basic(30)) == 1 // 2

## calculus
x,y = symbols("x y")
n = Basic(2)
ex = sin(x*y)
@test diff(log(x),x) == 1/x
@test diff(ex, x) == y * cos(x*y)
@test diff(ex, x, 2) == diff(diff(ex,x), x)
@test diff(ex, x, n) == diff(diff(ex,x), x)
@test diff(ex, x, y) == diff(diff(ex,x), y)
@test diff(ex, x, y,x) == diff(diff(diff(ex,x), y), x)
@test series(sin(x), x, 0, 2) == x
@test series(sin(x), x, 0, 3) == x - x^3/6

## ntheory
@test mod(Basic(10), Basic(4)) == 2
@test mod(Basic(10), 4) == 2               # mod(::Basic, ::Number)
@test_throws MethodError mod(10, Basic(4)) # no mod(::Number, ::Basic)
@test gcd(Basic(10), Basic(4)) == 2
@test lcm(Basic(10), Basic(4)) == 20
@test binomial(Basic(5), 2) == 10

## type information
a = Basic(1)
b = Basic(1//2)
@test isa(SymEngine.BasicType(a+a), SymEngine.BasicType{Val{:Integer}})
@test isa(SymEngine.BasicType(a+b), SymEngine.BasicType{Val{:Rational}})

## can we do math with items of BasicType?
a1 = SymEngine.BasicType(a)
tot = a1
for i in 1:100  tot = tot + a1 end
@test tot == 101
sin(a1)

## subs
ex = x^2 + y^2
@test subs(ex, x, 1) == 1 + y^2
@test subs(ex, (x, 1)) == 1 + y^2
@test subs(ex, x => 1) == 1 + y^2
@test subs(ex, (x,1), (y,2)) == 1 + 2^2
@test subs(ex, x => 1, y => 2) == 1 + 2^2

## lambidfy
@test_approx_eq lambdify(sin(Basic(1))) sin(1)
ex = sin(x)
@test_approx_eq lambdify(ex)(1) sin(1)
ex = exp(PI/2*x)
@test_approx_eq lambdify(ex)(1) exp(pi/2)

## N
a = Basic(1)
@test N(a) == 1
@test N(Basic(1//2)) == 1//2
@test N(Basic(12345678901234567890)) == 12345678901234567890

## generic linear algebra
x = symbols("x")
A = [x 2; x 1]
@test det(A) == -x
@test det(inv(A)) == - 1/x
(A \ [1,2])[1] == 3/x

## check that unique work (hash)
x,y,z = symbols("x y z")
@test length(SymEngine.free_symbols([x*y, y,z])) == 3


text:
include("getdigit_n.jl")
include("myprintf.jl")
using Printf
using .GetDigit_n
using .MyPrintf

const DIGIT = 100
const MAXN = 500

function machin(n)
    result = BigFloat(4) * myatan(BigFloat(1) / BigFloat(5), n) - myatan(BigFloat(1) / BigFloat(239), n)
    return BigFloat(4) * result
end

function myatan(x, n)
    result = BigFloat(0)
    xn = x

    for i = 0:n - 1
        denominator = BigFloat(2 * i + 1)

        if i % 2 != 0
            denominator *= BigFloat(- 1)
        end

        result += xn / denominator

        xn *= (x * x)
    end

    return result
end

function main()
    n = GetDigit_n.getn(machin, MAXN, DIGIT, true)
    precision = GetDigit_n.getdigit(machin, n, DIGIT, true)
    precision = GetDigit_n.getdigit_2(machin, n, DIGIT, precision)

    setprecision(precision)
    MyPrintf.myprintf("π（計算値） = ", DIGIT, machin(n))
    MyPrintf.myprintf("π（厳密値） = ", DIGIT, BigFloat(π))
    @printf("n = %d, 計算%d桁で、%d桁求まりました\n", n, precision, DIGIT)
end

main()


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module Fields

using ExtractMacro
using ..Interface
using ..Common

export GraphFields, GraphFieldsNormalDiscretized

import ..Interface: energy, delta_energy, neighbors, allΔE, delta_energy_residual

struct GraphFields{ET,LEV} <: DiscrGraph{ET}
    N::Int
    fields::Vector{ET}
    function GraphFields{ET,LEV}(fields::Vector{ET}) where {ET,LEV}
        LEV isa Tuple{Real,Vararg{Real}} || throw(ArgumentError("invalid level spec, expected a Tuple of Reals, given: $LEV"))
        all(f->f ∈ LEV, fields) || throw(ArgumentError("invalid field value, expected $LEV, given: $(findfirst(f->f ∉ LEV, fields))"))
        length(unique(LEV)) == length(LEV) || throw(ArgumentError("repeated levels in LEV: $LEV"))
        N = length(fields)
        return new{ET,LEV}(N, fields)
    end
end

"""
    GraphFields(N::Integer, LEV::Tuple = (1,)) <: DiscrGraph

A simple `DiscrGraph` type with `N` non-interacting variables, each of which is
subject to a local field. The fields are extracted at random from `LEV`, which
must be a `Tuple` of `Real`s.

Mostly useful for testing/debugging purposes.
"""
function GraphFields(N::Integer, LEV::Tuple{ET,Vararg{ET}}) where {ET<:Real}
    return GraphFields{ET,LEV}(rand(collect(LEV), N))
end

GraphFields(N::Integer) = GraphFields{Int,(1,)}(ones(Int, N))

function energy(X::GraphFields{Int,(1,)}, C::Config)
    @extract C : N s
    return N - 2 * sum(s)
end

function energy(X::GraphFields{ET}, C::Config) where {ET}
    @assert X.N == C.N
    @extract X : fields
    @extract C : N s

    local E::ET = zero(ET)
    for i = 1:N
        E -= fields[i] * (2s[i] - 1)
    end

    return E
end

function delta_energy(X::GraphFields{Int,(1,)}, C::Config, move::Int)
    @assert 1 ≤ move ≤ C.N
    @extract C : N s

    return 2 * (2 * s[move] - 1)
end

function delta_energy(X::GraphFields{ET}, C::Config, move::Int) where {ET}
    @assert 1 ≤ move ≤ C.N
    @assert X.N == C.N
    @extract X : fields
    @extract C : N s

    return convert(ET, 2 * fields[move] * (2s[move] - 1))
end

neighbors(X::GraphFields, i::Int) = return ()
#allΔE(::Type{GraphFields{Int,(1,)}}) = (2,)

@generated function allΔE(::Type{GraphFields{ET,LEV}}) where {ET,LEV}
    absLEV = sort!(unique(map(x->convert(ET, 2*abs(x)), LEV)))
    return Expr(:tuple, absLEV...)
end

struct GraphFieldsNormalDiscretized{ET,LEV} <: DoubleGraph{DiscrGraph{ET},Float64}
    N::Int
    X0::GraphFields{ET,LEV}
    rfields::Vec
    function GraphFieldsNormalDiscretized{ET,LEV}(N::Integer) where {ET,LEV}
        cfields = randn(N)
        fields, rfields = discretize(cfields, LEV)
        X0 = GraphFields{ET,LEV}(fields)
        return new{ET,LEV}(N, X0, rfields)
    end
end

"""
    GraphFieldsNormalDiscretized(N::Integer, LEV::Tuple) <: DoubleGraph{Float64,GraphFields}

A simple `DoubleGraph` type with `N` non-interacting variables, each of which is
subject to a local field. The fields are extracted independently from a normal
distribution with unit variance, and then are discretized using the values in `LEV`, which
must be a `Tuple` of `Real`s.

Mostly useful for testing/debugging purposes.
"""
GraphFieldsNormalDiscretized(N::Integer, LEV::Tuple{ET,Vararg{ET}}) where {ET<:Real} = GraphFieldsNormalDiscretized{ET,LEV}(N)

function energy(X::GraphFieldsNormalDiscretized, C::Config)
    @assert X.N == C.N
    @extract X : X0 rfields
    @extract C : N s

    E0 = energy(X0, C)

    E1 = 0.0
    for i = 1:N
        E1 -= rfields[i] * (2s[i] - 1)
    end

    return convert(Float64, E0 + E1)
end

function delta_energy_residual(X::GraphFieldsNormalDiscretized, C::Config, move::Int)
    @assert 1 ≤ move ≤ C.N
    @assert X.N == C.N
    @extract X : rfields
    @extract C : N s

    return 2 * rfields[move] * (2s[move] - 1)
end

function delta_energy(X::GraphFieldsNormalDiscretized, C::Config, move::Int)
    ΔE0 = delta_energy(X.X0, C, move)
    ΔE1 = delta_energy_residual(X, C, move)
    return convert(Float64, ΔE0 + ΔE1)
end

neighbors(X::GraphFieldsNormalDiscretized, i::Int) = return ()

end # module


text:
# run hopper_in_place_hop.jl

include(joinpath(@__DIR__, "..", "src/dynamics", "hopper_3D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# get hopper model
# model_sim = get_model("hopper_2D", surf="sinusoidal")
model_sim = get_model("hopper_3D", surf="flat")
nq = model_sim.dim.q
nu = model_sim.dim.u
nc = model_sim.dim.c
nb = model_sim.dim.b
nw = model_sim.dim.w

# time
H = 92#traj.H
h = 0.01#traj.h
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 5000

# Select the initial speed
v0 = [0.0; 0.2]
Tstance = 0.13 # measure using hop-in-place gait
Tflight = 0.62 # measure using hop-in-place gait
include(joinpath(module_dir(), "src/controller/raibert_3D_policy.jl"))
p = raibert_policy(model_sim, v0=v0, Tstance=Tstance, Tflight=Tflight, h=h)

off0 = SVector{nq,T}([0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.5])
off1 = SVector{nq,T}([v0[1]*h_sim, v0[2]*h_sim, 0.5, 0.0, 0.0, 0.0, 0.5])
q_ref = SVector{nq,T}([0.0, 0.0, 0.5, 0.0, 0.0, 0.0, 0.0])
q0_sim = copy(q_ref) + off0
q1_sim = copy(q_ref) + off1

# Disturbances
w = [zeros(nw) for t=1:Int(ceil(H_sim/N_sample))]
# w[75] += [5.0, 5.0, 0.0]
d = open_loop_disturbances(w)

sim = simulator(model_sim, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
    d = d,
    ip_opts = InteriorPointOptions(
        r_tol = 1.0e-8,
        κ_init = 1.0e-8,
        κ_tol = 2.0e-8),
    sim_opts = SimulatorOptions(warmstart = true)
    )

@time status = simulate!(sim)

# plt = plot(layout=(3,1), legend=false)
# plot!(plt[1,1], hcat(Vector.(vcat([fill(ref_traj.q[i], N_sample) for i=1:H]...))...)',
#     color=:red, linewidth=3.0)
# plot!(plt[1,1], hcat(Vector.(sim.traj.q)...)', color=:blue, linewidth=1.0)
# plot!(plt[2,1], hcat(Vector.(vcat([fill(ref_traj.u[i][1:nu], N_sample) for i=1:H]...))...)',
#     color=:red, linewidth=3.0)
# plot!(plt[2,1], hcat(Vector.([u[1:nu] for u in sim.traj.u]*N_sample)...)', color=:blue, linewidth=1.0)
# plot!(plt[3,1], hcat(Vector.([γ[1:nc] for γ in sim.traj.γ]*N_sample)...)', color=:blue, linewidth=1.0)
#
# plot_lines!(vis, model_sim, sim.traj.q[1:10:end])
# plot_surface!(vis, model_sim.env, n=200)
anim = visualize_robot!(vis, model_sim, sim.traj, sample=20)

# hopper_3D_euler = Hopper3D(Dimensions(nq, nu, nw, nc, nb),
# 		mb, ml, Jb, Jl,
# 		μ_world, μ_joint, g,
# 		:RotXYX, # Euler angles
# 		BaseMethods(), DynamicsMethods(), ContactMethods(),
# 		ResidualMethods(), ResidualMethods(),
# 		SparseStructure(spzeros(0, 0), spzeros(0, 0)),
# 		SVector{7}(zeros(7)),
# 		environment_3D(x -> 0.075 * sin(2π * x[1])),
# 	    )
#
# eul = RotXYX(0.5 * π, 0.0, 0.0)
# qq = Array(copy(q_ref))
# qq[4:6] = [eul.theta1; eul.theta2; eul.theta3]
#
# build_robot!(vis, hopper_3D_euler, name = model_name(hopper_3D_euler))
# set_robot!(vis, model_sim, q1_sim, name = model_name(model_sim))
#
# anim = visualize_force!(vis, model_sim, sim.traj, anim=anim, h=h_sim, sample=20)
#
# plot(hcat([q[[1,2,3,4]] for q in Vector.(sim.traj.q)[1:1200]]...)')
# plot(hcat([q[[3]] for q in Vector.(sim.traj.q)[1:end]]...)')
# plot(hcat([γ[[1]] for γ in Vector.(sim.traj.γ)[1:end]]...)')
# plot(hcat([u[[1,2]] for u in Vector.(sim.traj.u)[1:end]]...)')

# filename = "raibert_hopper_push"
# MeshCat.convert_frames_to_video(
#     "/home/simon/Downloads/$filename.tar",
#     "/home/simon/Documents/$filename.mp4", overwrite=true)
#
# convert_video_to_gif(
#     "/home/simon/Documents/$filename.mp4",
#     "/home/simon/Documents/$filename.gif", overwrite=true)


text:
function load_bson_model_file(path_to_model_file::String, module_name::Module)

    try
        return BSON.load(path_to_model_file, module_name)

    catch error

        # what is our error policy? => for now, just print the message
        error_message = sprint(showerror, error, catch_backtrace())
        println(error_message)

    end
end