text:
## Packages
using Plots

## Files
include("FEM_aiyagari.jl")
include("NewtonRoot_mult.jl")

############ Plot basis functions #################
plot(x -> ψ(x, [0, 1, 3, 6],1), xlim=(0, 6), label = "ψ_1(x)")
plot!(x -> ψ(x, [0, 1, 3, 6], 2), xlim=(0, 6), label = "ψ_2(x)")
plot!(x -> ψ(x, [0, 1, 3, 6], 3), xlim=(0, 6), label = "ψ_3(x)")
plot!(x -> ψ(x, [0, 1, 3, 6], 4), xlim=(0, 6), label = "ψ_4(x)")
savefig("figs/linear_basis.png")

###################### Aiyagari ###########################
hh = Household()
θ_init = [hh.aGrid; hh.aGrid .+ 0.5]
θ_est = NewtonRoot_mult(θ->G(hh,θ),θ_init)

plot(a->a,xlim = (0.0,6.1),label = "45 degree", color = :black, linestyle = :dash,legend =:topleft)
plot!(a->get_apol(hh,a,hh.lGrid[1],θ_est),xlim = (0.0,6.1),label = "low")
plot!(a->get_apol(hh,a,hh.lGrid[2],θ_est),xlim = (0.0,6.1),label = "high")
scatter!(hh.aGrid,θ_est[1:8],label = "coefficient of basis (Low)")
scatter!(hh.aGrid,θ_est[9:16],label = "coefficient of basis (High)")
savefig("figs/aiyagari.png")


text:
#!/usr/bin/env julia

# TODO: Register package.
include("src/PCGLNS.jl")

"""
Optional Flags -- values are given in square brackets []
	-max_time=[Int]				 (default set by mode)
	-trials=[Int]				 (default set by mode)
	-restarts=[Int]              (default set by mode)
	-mode=[default, fast, slow]  (default is default)
	-verbose=[0, 1, 2, 3]        (default is 3.  0 is no output, 2 is most verbose)
	-output=[filename]           (default is output.tour)
	-epsilon=[Float in [0,1]]	 (default is 0.5)
	-reopt=[Float in [0,1]]      (default is set by mode)
"""
function parse_cmd(ARGS)
    if isempty(ARGS)
        println("no input instance given")
        exit(0)
    end
    if ARGS[1] == "-help" || ARGS[1] == "--help"
        println("Usage:  GTSPcmd.jl [filename] [optional flags]\n")
        println("Optional flags (vales are give in square brackets) :\n")
        println("-mode=[default, fast, slow]      (default is default)")
        println("-max_time=[Int]                  (default set by mode)")
        println("-trials=[Int]                    (default set by mode)")
        println("-restarts=[Int]                  (default set by mode)")
        println("-noise=[None, Both, Subset, Add] (default is Both)")
        println("-num_iterations=[Int]            (default set by mode. Number multiplied by # of sets)")
        println("-verbose=[0, 1, 2, 3]            (default is 3. 0 is no output, 3 is most.)")
        println("-output=[filename]               (default is None)")
        println("-epsilon=[Float in [0,1]]        (default is 0.5)")
        println("-reopt=[Float in [0,1]]          (default is 1.0)")
        println("-budget=[Int]                    (default has no budget)")
        exit(0)
    end
    int_flags =
        ["-max_time", "-trials", "-restarts", "-verbose", "-budget", "-num_iterations", "-seed"]
    float_flags = ["-epsilon", "-reopt"]
    string_flags = ["-mode", "-output", "-noise", "-devel"]
    filename = ""
    optional_args = Dict{Symbol,Any}()
    for arg in ARGS
        temp = split(arg, "=")
        if length(temp) == 1 && filename == ""
            filename = temp[1]
        elseif length(temp) == 2
            flag = temp[1]
            value = temp[2]
            if flag in int_flags
                key = Symbol(flag[2:end])
                optional_args[key] = parse(Int64, value)
            elseif flag in float_flags
                key = Symbol(flag[2:end])
                optional_args[key] = parse(Float64, value)
            elseif flag in string_flags
                key = Symbol(flag[2:end])
                optional_args[key] = value
            else
                println(
                    "WARNING: skipping unknown flag ",
                    flag,
                    " in command line arguments",
                )
            end
        else
            error("argument ", arg, " not in proper format")
        end
    end
    if !isfile(filename)
        println("the problem instance  ", filename, " does not exist")
        exit(0)
    end
    return filename, optional_args
end

# Running the code on the problem instance passed in via command line args.
function solve(ARGS)
    problem_instance, optional_args = parse_cmd(ARGS)
    PCGLNS.solver(problem_instance; optional_args...)
end

solve(ARGS)

ext: jl

max_stars_repo_name: AndreiKud/PCGLNS

max_issues_repo_name: AndreiKud/PCGLNS

max_forks_repo_name: AndreiKud/PCGLNS


text:
##################### Lazy multi-dimensional grid ##############################
# useful to automatically detect Kronecker structure in Kernel matrices
struct LazyGrid{T, V} <: AbstractVector{Vector{T}}
    args::V
end
function LazyGrid(args)
    T = promote_type(eltype.(args)...)
    LazyGrid{T, typeof(args)}(args)
end
LazyGrid(args...) = LazyGrid(args)
LazyGrid(x::AbstractArray, d::Int) = LazyGrid(Fill(x, d))

Base.length(G::LazyGrid) = prod(length, G.args)
Base.eltype(G::LazyGrid{T}) where {T} = Vector{T}
Base.size(G::LazyGrid) = (length(G),)
function Base.ndims(G::LazyGrid)
    sum(x -> x isa AbstractVector ? 1 : size(x, 1), G.args)
end
# default column-major indexing
function Base.getindex(G::LazyGrid{T}, i::Integer) where {T}
    @boundscheck checkbounds(G, i)
    val = zeros(T, ndims(G))
    n = length(G)
    d = ndims(G)
    @inbounds for a in reverse(G.args)
        n ÷= length(a)
        ind = cld(i, n) # can this be replaced with fld1, mod1?
        if a isa AbstractVector
            val[d] = a[ind]
            d -= 1
        else
            val[d-size(a, 1)+1:d] = a[:, ind]
            d -= size(a, 1)
        end
        i = mod1(i, n) # or some kind of shifted remainder?
    end
    return val
end
# row major indexing
function Base.getindex(G::LazyGrid{T}, i::Integer, rowmajor::Val{true}) where {T}
    @boundscheck checkbounds(G, i)
    val = zeros(T, ndims(G))
    n = length(G)
    d = 1
    @inbounds for (j, a) in enumerate(G.args)
        n ÷= length(a)
        ind = cld(i, n) # can this be replaced with fld1, mod1?
        if a isa AbstractVector
            val[d] = a[ind]
            d += 1
        else
            val[d:d+size(a, 1)-1] = a[:, ind]
            d += size(a, 1)
        end
        i = mod1(i, n) # or some kind of shifted remainder?
    end
    return val
end


text:
# ---
# title: 17. Letter Combinations of a Phone Number
# id: problem17
# author: zhwang
# date: 2022-03-17
# difficulty: Medium
# categories: String, Backtracking, Depth-first Search, Recursion
# link: <https://leetcode.com/problems/letter-combinations-of-a-phone-number/description/>
# hidden: true
# ---
#
# Given a string containing digits from `2-9` inclusive, return all possible
# letter combinations that the number could represent. Return the answer in
# **any order**.
#
# A mapping of digit to letters (just like on the telephone buttons) is given
# below. Note that 1 does not map to any letters.
#
# ![](https://upload.wikimedia.org/wikipedia/commons/thumb/7/73/Telephone-
# keypad2.svg/200px-Telephone-keypad2.svg.png)
#
#
#
# **Example 1:**
#
#
#
#     Input: digits = "23"
#     Output: ["ad","ae","af","bd","be","bf","cd","ce","cf"]
#
#
# **Example 2:**
#
#
#
#     Input: digits = ""
#     Output: []
#
#
# **Example 3:**
#
#
#
#     Input: digits = "2"
#     Output: ["a","b","c"]
#
#
#
#
# **Constraints:**
#
#   * `0 <= digits.length <= 4`
#   * `digits[i]` is a digit in the range `['2', '9']`.
#
#
## @lc code=start
using LeetCode

function letter_combinations(digits::String)::Vector{String}
    isempty(digits) && return String[]
    table = Dict{Char, String}('2'=>"abc", '3'=>"def", '4'=>"ghi", '5'=>"jkl", '6'=>"mno", '7'=>"pqrs", '8'=>"tuv", '9'=>"wxyz")
    res = [""]
    for i in digits
        res = [s1 * s2 for s1 in res for s2 in table[i]]
    end
    res
end
## @lc code=end

## @lc test=start
@testset "17.letter-combinations-of-a-phone-number.jl" begin
    @test letter_combinations("23") == ["ad","ae","af","bd","be","bf","cd","ce","cf"]
    @test isempty(letter_combinations(""))
    @test letter_combinations("8") == ["t", "u", "v"]
end
## @lc test=end

## @lc info=start
# link: [solution to 17]()
## @lc info=end


text:
using Documenter, RCall

makedocs(
    modules = [RCall],
    format   = :html,
    sitename = "RCall.jl",
    pages    = [
        "Introduction" => "index.md",
        "Installation" => "installation.md",
        "Getting Started" => "gettingstarted.md",
        "Supported Conversions" => "conversions.md",
        "Custom Conversion" => "custom.md",
        "Eventloop" => "eventloop.md",
        "Internal" => "internal.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaInterop/RCall.jl.git",
    target = "build",
    julia  = "0.7",
    deps = nothing,
    make = nothing
)


text:
using OVERTVerify

state_vars = [:x1, :x2]
control_vars = [:u1]
input_set = Dict(:x1=>[1., 1.2], :x2=>[0., 0.2])
controller_file = "nnet_files/jmlr/single_pendulum_small_controller.nnet"
dynamics_map = Dict(OVERTVerify.single_pend_input_vars[1]=>OVERTVerify.single_pend_input_vars[2],
                    OVERTVerify.single_pend_input_vars[2]=>OVERTVerify.single_pend_θ_doubledot)
dt = 0.1
output_constraints = [:(x1 <= -0.2167)] # (avoid set)
N_steps=25
experiment_name = "single_pendulum_small_controller"
dirname="src/examples/jmlr/comparisons/"

ΔT = compare_to_dreal(state_vars, control_vars, input_set, controller_file, dynamics_map, dt, output_constraints, dirname, experiment_name, N_steps; jobs=28, dreal_path="/opt/dreal/4.20.12.1/bin/dreal")


text:
# julia -i server.jl

using Bukdu
const ServerHost = "localhost"
const ServerPort = 8080

struct WSController <: ApplicationController
    conn::Conn
end

function websocket(::WSController)
    js = """
var wsUri = "ws://$ServerHost:$ServerPort/";
var output;
var m;

function init() {
  output = document.getElementById("output");
  testWebSocket();
  m = window.mode_objects['4D'];
  console.log("m ", m)
}

function testWebSocket() {
  websocket = new WebSocket(wsUri);
  websocket.onopen = function(evt) { onOpen(evt) };
  websocket.onclose = function(evt) { onClose(evt) };
  websocket.onmessage = function(evt) { onMessage(evt) };
  websocket.onerror = function(evt) { onError(evt) };
}

function onOpen(evt) {
  println("CONNECTED");
  msg = JSON.stringify("WebSocket rocks");
  doSend(msg);
}

function onClose(evt) {
  println("DISCONNECTED");
}

function cc102(nt) {
  obj = m.leftScene;
  switch (nt.value) {
      case 0x01:
          println("0x01");
          obj.rotation.x += 0.1;
          break;
      case 0x7f:
          println("0x7f");
          obj.rotation.x -= 0.1;
          break;
  }
}

function cc103(nt) {
  obj = m.leftScene;
  switch (nt.value) {
      case 0x01:
          obj.rotation.y += 0.1;
          break;
      case 0x7f:
          obj.rotation.y -= 0.1;
          break;
  }
}

function midi_button(nt) {
  switch (nt.name) {
      case "cc102": // knob 1
          cc102(nt);
          break;
      case "cc103": // knob 2
          cc103(nt);
          break;
     default:
          break;
  }
}

function midi_callback(nt) {
  switch (nt.action) {
     case "button":
        midi_button(nt);
        break;
     default:
        break;
  }
}

function onMessage(evt) {
  nt = JSON.parse(evt.data);
  println('RESPONSE: ' + JSON.stringify(nt) + '');
  midi_callback(nt);
  // websocket.close();
}

function onError(evt) {
  println('ERROR: ' + evt.data);
}

function doSend(message) {
  println("SENT: " + message);
  websocket.send(message);
}

function println(message) {
  var pre = document.createElement("p");
  pre.style.wordWrap = "break-word";
  pre.innerHTML = message;
  output.appendChild(pre);
}

window.addEventListener("load", init, false);
"""
    render(JavaScript, js)
end

plug(Plug.Logger, access_log=(path=normpath(@__DIR__, "access.log"),), formatter=Plug.LoggerFormatter.datetime_message)

routes(:front) do
    get("/websocket.js", WSController, websocket)
    plug(Plug.WebSocket)
    plug(Plug.Static, at="/", from=normpath(@__DIR__))
end

Bukdu.start(ServerPort, host=ServerHost)

function midi_event(msg)
    buttons = Dict(
        0x47 => "cc102",
    )
    name = get(buttons, msg[2], "")
    value = msg[3]
    (action="button", name=name, value=value)
end

using JSON2
function midi_callback(msg)
    # @info midi_callback first(msg)
    # if 0x90 == first(msg)
    if true
        websocks = Bukdu.websockets()
        if isempty(websocks)
            @error "please visit the url in the web browser"
        else
            ws = first(websocks)
            nt = midi_event(msg)
            written = write(ws, JSON2.write(nt))
            (written, nt)
        end
    end
end

#=
julia> midi_callback(UInt8[0xb0, 0x47, 0x01])
(46, (action = "change", name = "cc102", value = 0x7f))
=#

#=
import PushInterface: RTMIDI
import .RTMIDI: rtmidi_in_create_default, rtmidi_open_port, rtmidi_in_set_callback
import .RTMIDI: destroy, EventCB, rtmidi_callback_func

device = rtmidi_in_create_default()
rtmidi_open_port(device, 1, "Ableton Push 2 User Port")
if !(device[].ok)
    @error :rtmidi_open_port "could not found ableton push 2"
    exit(1)
end

cb_ptr = @cfunction rtmidi_callback_func Cvoid (Cdouble, Ptr{Cuchar}, Ptr{EventCB})
cond = Base.AsyncCondition()
handle = Base.unsafe_convert(Ptr{Cvoid}, cond)
ecb = EventCB(handle, 0, C_NULL)
r_ecb = Ref(ecb)
rtmidi_in_set_callback(device, cb_ptr, r_ecb)
finalizer(destroy, device)

function _callback_async_loop(cond, r_ecb)
    while isopen(cond)
        wait(cond)
        msg = codeunits(unsafe_string(r_ecb[].message))
        midi_callback(msg)
    end
end

cbloop = Task() do
    _callback_async_loop(cond, r_ecb)
end
schedule(cbloop)
=#

Base.JLOptions().isinteractive==0 && wait()


text:
using MathProgComplex
using Base.Test

problems = [lasserre_ex1, lasserre_ex2, lasserre_ex5, test_gloptipoly_ex1, test_gloptipoly_ex2]

brokentestset = Set([(:MomentRelaxation, lasserre_ex2, 3), (:MomentRelaxation, lasserre_ex5, 2)])

@testset "SDP hierarchy base tests" begin

    @testset "$pbsolved tests" for pbsolved in [:SOSRelaxation, :MomentRelaxation]
        @testset "$problem_fct problem" for problem_fct in problems
            problem, order_to_obj, ε_abs = problem_fct()

            @testset "Order $d" for d in keys(order_to_obj)
                relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                            d = d,
                                            params = Dict(:opt_outlev=>0,
                                                        :opt_solver=>testsolver,
                                                        :opt_pbsolved=>pbsolved))

                primobj, dualobj = run_hierarchy(problem, relax_ctx)

                @show pbsolved, problem_fct, d
                @show d, primobj, dualobj, order_to_obj[d]

                if (pbsolved, problem_fct, d) in brokentestset
                    @test_broken primobj ≈ order_to_obj[d] atol=ε_abs
                else
                    @test primobj ≈ order_to_obj[d] atol=ε_abs
                    @test dualobj ≈ primobj atol=ε_abs*min(abs(primobj), abs(dualobj))
                end
            end
        end
    end
end


text:
struct NonlinearSystem <: AbstractSystem
    eqs::Vector{Operation}
    vs::Vector{Variable}
    ps::Vector{Variable}
    v_name::Vector{Symbol}
    p_name::Symbol
end

function NonlinearSystem(eqs, vs, ps;
                         v_name = :Variable,
                         dv_name = :DependentVariable,
                         p_name = :Parameter)
    NonlinearSystem(eqs, vs, ps, [v_name,dv_name], p_name)
end

function NonlinearSystem(eqs;
                         v_name = :Variable,
                         dv_name = :DependentVariable,
                         p_name = :Parameter)
    # Allow the use of :DependentVariable to make it seamless with DE use
    targetmap = Dict(v_name => v_name, dv_name => v_name, p_name => p_name)
    vs, ps = extract_elements(eqs, targetmap)
    NonlinearSystem(eqs, vs, ps, [v_name,dv_name], p_name)
end

function generate_nlsys_function(sys::NonlinearSystem)
    var_exprs = [:($(sys.vs[i].name) = u[$i]) for i in 1:length(sys.vs)]
    param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in 1:length(sys.ps)]
    sys_idxs = map(eq->isequal(eq.args[1],Constant(0)),sys.eqs)
    sys_eqs = sys.eqs[sys_idxs]
    calc_eqs = sys.eqs[.!(sys_idxs)]
    calc_exprs = [:($(Symbol("$(eq.args[1].name)")) = $(eq.args[2])) for eq in calc_eqs]
    sys_exprs = [:($(Symbol("resid[$i]")) = $(sys_eqs[i].args[2])) for i in eachindex(sys_eqs)]

    exprs = vcat(var_exprs,param_exprs,calc_exprs,sys_exprs)
    block = expr_arr_to_block(exprs)
    :((du,u,p)->$(block))
end

function calculate_jacobian(sys::NonlinearSystem,simplify=true)
    sys_idxs = map(eq->isequal(eq.args[1],Constant(0)),sys.eqs)
    sys_eqs = sys.eqs[sys_idxs]
    calc_eqs = sys.eqs[.!(sys_idxs)]
    rhs = [eq.args[2] for eq in sys_eqs]

    for i in 1:length(calc_eqs)
        find_replace!.(rhs,calc_eqs[i].args[1],calc_eqs[i].args[2])
    end

    sys_exprs = calculate_jacobian(rhs,sys.vs)
    sys_exprs = Expression[expand_derivatives(expr) for expr in sys_exprs]
    if simplify
        sys_exprs = Expression[simplify_constants(expr) for expr in sys_exprs]
    end
    sys_exprs
end

function generate_nlsys_jacobian(sys::NonlinearSystem,simplify=true)
    var_exprs = [:($(sys.vs[i].name) = u[$i]) for i in 1:length(sys.vs)]
    param_exprs = [:($(sys.ps[i].name) = p[$i]) for i in 1:length(sys.ps)]
    jac = calculate_jacobian(sys,simplify)
    jac_exprs = [:(J[$i,$j] = $(Expr(jac[i,j]))) for i in 1:size(jac,1), j in 1:size(jac,2)]
    exprs = vcat(var_exprs,param_exprs,vec(jac_exprs))
    block = expr_arr_to_block(exprs)
    :((J,u,p,t)->$(block))
end

export NonlinearSystem
export generate_nlsys_function


text:
function genie_config()
    # Wrap everything for setting up Genie here for ease
    Genie.config.websockets_server = true
    Genie.config.server_host = "127.0.0.1"
    Genie.config.server_port = 8000
    Genie.config.log_level = Logging.Error
    Genie.config.log_to_file = false
    Genie.config.server_document_root = included_file_path(joinpath("files"))
end

"""
    run_app(app; async = true)

The`app` module that defines your app must have a function called `ui` that defines
the UI. The functions in the module (which don't need to be exported) define
the logic for your app. The name of the function is the name of the outputs in
your UI, and the inputs to those functions correspond to the `id`s of inputs in
the UI.

See the `getting started` guide for help on how to define your app.
"""
function run_app(app::Module; async = true)
    sanitycheck(app)
    server = app.Server
    ui = app.ui
    sessions = create_sessions(app)

    genie_config()

    Genie.Router.route("/") do
        generate_template(app.title, ui, server)
    end

    Genie.Router.channel("/matte/api") do
        json_request = Genie.Router.@params(:payload)

        if haskey(json_request, "id")
            session = sessions[json_request["session_id"]]
            session.genie_wsclient = Genie.Requests.wsclient()
            handle_request(json_request["id"], json_request["input"], server, session)
        end
    end

    println(Genie.config.websockets_server)

    Genie.up(async = async)
end

"""
    stop_app()

Stops a running server that was started with `async=true`
"""
function stop_app()
    @async Base.throwto(Genie.AppServer.SERVERS.webserver, InterruptException())
    @async Base.throwto(Genie.AppServer.SERVERS.websockets, InterruptException())
end


text:
using DSJulia

pn = "../submissions/Day1/Basic/"

fns = readdir(pn)

trackers = []
for fn in fns
  tracker = grade(pn*fn)
  push!(trackers, tracker)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Functions to convert osculating elements to mean elements.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A; Crawford, P (2008). SGP4 orbit determination. AIAA/AAS
#       Astrodynamics Specialist Conference, Honoulu, HI.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export rv_to_mean_elements_sgp4, rv_to_tle

"""
    rv_to_mean_elements_sgp4(vJD::AbstractVector{T}, vr::AbstractVector{Tv}, vv::AbstractVector{Tv}, W = I; max_it = 10000, sgp4_gc = sgp4_gc_wgs84) where {T,Tv<:AbstractVector}

Compute the mean elements for SGP4 based on the position `vr` and velocity
vectors `vr` represented in TEME reference frame. The epoch of those
measurements [Julian Day] must be in `vJD`.

The matrix `W` defined the weights for the least-square algorithm.

The variable `max_it` defines the maximum allowed number of iterations.

The variable `sgp4_gc` defines which constants should be used when running SGP4.

# Returns

* The epoch of the elements [Julian Day].
* The mean elements for SGP4 algorithm:
    - Mean motion [rad/s];
    - Eccentricity [];
    - Inclination [rad];
    - Right ascension of the ascending node [rad];
    - Argument of perigee [rad];
    - Mean anomaly [rad];
    - BSTAR.
* The covariance matrix of the mean elements estimation.

"""
function rv_to_mean_elements_sgp4(vJD::AbstractVector{T}, vr::AbstractVector{Tv},
                                  vv::AbstractVector{Tv}, W = I; max_it = 10000,
                                 sgp4_gc = sgp4_gc_wgs84) where
    {T,Tv<:AbstractVector}

    # Number of measurements.
    num_meas   = length(vr)

    # Number of states that will be fitted. In this cases they are:
    #
    #   - Mean motion;
    #   - Eccentricity;
    #   - Inclination;
    #   - Right ascension of the ascending node;
    #   - Argument of perigee;
    #   - Mean anomaly;
    #   - BSTAR.
    num_states = 7

    # Initial guess of the mean elements.
    #
    # NOTE: x₁ is the previous estimate and x₂ is the current estimate.
    #
    # TODO: Should we compute the mean motion using a better algorithm?
    o = rv_to_kepler(vr[end], vv[end])
    x₁ = @SVector T[sqrt(m0/o.a^3), o.e, o.i, o.Ω, o.ω, f_to_M(o.e, o.f), 0]
    x₂ = x₁
    P  = SMatrix{num_states,num_states,T}(Matrix(I,num_states,num_states))

    # Loop until the maximum allowed iteration.
    @inbounds for it = 1:max_it
        x₁ = x₂

        # Variables to store the summations to compute the least square fitting
        # algorithm.
        ΣAᵀWA = @SMatrix zeros(num_states, num_states)
        ΣAᵀWb = @SVector zeros(num_states)

        # We will interpolate backwards in time, so that the mean elements
        # in x₂ are related to the newest measurement.
        for k = num_meas:-1:1
            # Obtain the measured ephemerides.
            y = vcat(vr[k],vv[k])

            # Obtain the computed ephemerides considering the current estimate
            # of the mean elements.
            Δt = (vJD[k] - vJD[end])*86400
            r̂, v̂, ~ = _sgp4_si(Δt, sgp4_gc, vJD[k], x₁...)
            ŷ = vcat(r̂,v̂)

            # Compute the error.
            b = y-ŷ

            # Compute the Jacobian.
            A = _sgp4_jacobian(Δt, vJD[end], x₁)

            # Accumulation.
            ΣAᵀWA += A'*W*A
            ΣAᵀWb += A'*W*b
        end

        # Update the estimate.
        P  = pinv(ΣAᵀWA)
        δx = P*ΣAᵀWb
        x₂ = x₁ + δx

        # Make sure that all parameters are inside the allowed intervals.
        if x₂[2] < 0
            x₂ = setindex(x₂, -x₂[2], 2)
        end

        # Obtain the maximum deviation.
        σ_max = maximum( abs.( (x₂-x₁) ) )

        println("PROGRESS: $it / $max_it, σ_max = $σ_max")

        # If the threshold is achieved, then just stop the loop.
        σ_max < 1e-6 && break
    end

    # Return the mean elements for SGP4 and the covariance matrix.
    return vJD[end], x₂, P
end

"""
    rv_to_tle(args...; name::String = "UNDEFINED", sat_num::Int = 9999, classification::Char = 'U', int_designator = "999999", elem_set_number::Int = 0, rev_num, kwargs...)

Convert a set of position and velocity vectors represented in TEME reference
frame to a TLE. The arguments `args` and keywords `kwargs` are the same as those
described in the function `rv_to_mean_elements_sgp4`.

Additionally, the user can specify some parameters of the generated TLE.

This function prints the TLE to `stdout` using the function `print_tle` and also
returns the TLE string.

"""
function rv_to_tle(args...;
                   name::String = "UNDEFINED",
                   sat_num::Int = 9999,
                   classification::Char = 'U',
                   int_designator = "999999",
                   elem_set_number::Int = 0,
                   rev_num::Int = 0,
                   kwargs...)

    # Convert the position and velocity vectors to mean elements.
    JD,x,P = rv_to_mean_elements_sgp4(args...; kwargs...)

    # Compute the data as required by the TLE format.
    dt  = JDtoDate(DateTime, JD)
    dt₀ = JDtoDate(DateTime, DatetoJD(year(dt), 1, 1, 0, 0, 0))

    dt_year    = year(dt)
    epoch_year = dt_year < 1980 ? dt_year - 1900 : dt_year - 2000
    epoch_day  = (dt-dt₀).value/1000/86400 + 1

    # Construct the TLE.
    tle = TLE(name,
              sat_num,
              classification,
              int_designator,
              epoch_year,
              epoch_day,
              JD,
              0.0,
              0.0,
              x[7],
              elem_set_number,
              0,
              x[3]*180/pi,
              x[4]*180/pi,
              x[2],
              x[5]*180/pi,
              x[6]*180/pi,
              x[1]*43200/pi,
              rev_num,
              0)

    tle_str = tle_to_str(tle)

    # Print the TLE to the `stdout`.
    print(tle_str)

    # Return the TLE string.
    return tle_str
end

################################################################################
#                              Private functions
################################################################################

# Compute the SGP4 algorithm considering all variables in SI.
function _sgp4_si(Δt, sgp4_gc, epoch, n_0, e_0, i_0, Ω_0, ω_0, M_0, bstar)
    r,v,sgp4d = sgp4(Δt/60, sgp4_gc, epoch, 60n_0, e_0, i_0, Ω_0, ω_0, M_0, bstar)
    return 1000r, 1000v, sgp4d
end

function _sgp4_jacobian(Δt, JD₀, x₀::AbstractVector{T}; ϵ = 1e-10, sgp4_gc = sgp4_gc_wgs84) where T
    num_states = 7
    dim_obs    = 6
    M          = Matrix{T}(undef, dim_obs, num_states)

    # Auxiliary variables.
    x₁ = copy(x₀)
    x₂ = copy(x₀)

    @inbounds for i = 1:dim_obs, j = 1:num_states
        x₁     = setindex(x₁, x₁[j] - ϵ, j)
        x₂     = setindex(x₂, x₂[j] + ϵ, j)

        r,v,~  = _sgp4_si(Δt, sgp4_gc, JD₀, x₁...)
        f_x₁   = vcat(r,v)
        r,v,~  = _sgp4_si(Δt, sgp4_gc, JD₀, x₂...)
        f_x₂   = vcat(r,v)

        M[i,j] = (f_x₂[i]-f_x₁[i])/(2ϵ)

        x₁     = setindex(x₁, x₀[j], j)
        x₂     = setindex(x₂, x₀[j], j)
    end

    return M
end


text:
#References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.

"""
    𝐒 = numSet([𝐂₀,𝐂₁,𝐂₂])

Create a numerical component set 'numSet' by providing a vector of
numerical canonical triplets.

# Examples
```@example
using ISA
ψ₀ = AMFMcomp(t->exp(-t^2),t->2.0,0.0)
ψ₁ = AMFMcomp(t->1.0,t->10*t,0.1)
ψ₂ = AMFMcomp(t->0.8*cos(2t),t->10 + 7.5*sin(t),π)
fs = 16_000
t = 0:1/fs:1
𝚿₀ = numComp( ψ₀(t), fs )
𝚿₁ = numComp( ψ₁(t), fs )
𝚿₂ = numComp( ψ₂(t), fs )
𝐂₀ = AMFMdemod(𝚿₀)
𝐂₁ = AMFMdemod(𝚿₁)
𝐂₂ = AMFMdemod(𝚿₂)
𝐒 = numSet([𝐂₀,𝐂₁,𝐂₂])
```
"""
struct numSet
  S::Vector{numTriplet}
end

# DISPLAY
function displayS(x)
  text = "\n"
  for i in 1:length(x.S)
    T = typeof(x.S[i])
    text *= "numerical canonical triplet #$(i)\n"
    for (name, typ) in zip(fieldnames(T), T.types)
        text *= "$name is $typ\n"
    end
    text *= "\n"
  end
  return text
end
Base.show(io::IO, x::numSet) = print(io, "numerical component set with $(length(x.S)) components$(displayS(x))")

"""
    𝐳 = numModel([𝚿₀,𝚿₁,𝚿₂])

Create a 'numModel' by providing a vector of numerical components
'numComp'.

# Examples
```@example
using ISA
ψ₀ = AMFMcomp(t->exp(-t^2),t->2.0,0.0)
ψ₁ = AMFMcomp(t->1.0,t->10*t,0.1)
ψ₂ = AMFMcomp(t->0.8*cos(2t),t->10 + 7.5*sin(t),π)
fs = 16_000
t = 0:1/fs:1
𝚿₀ = numComp( ψ₀(t), fs )
𝚿₁ = numComp( ψ₁(t), fs )
𝚿₂ = numComp( ψ₂(t), fs )
𝐳 = numModel([𝚿₀,𝚿₁,𝚿₂])
```
"""
struct numModel
  𝚿ₖ::Vector{numComp}
end

# METHODS
function (𝐳::numModel)(t::Real)
  out = 0.0
  for 𝚿 ∈ 𝐳.𝚿ₖ
    if 𝚿.t[1]<=t<=𝚿.t[end]
      out += 𝚿(t)
    end
  end
  return out
end
function (𝐳::numModel)(t::Vector{<:Real})
  return 𝐳.(t)
end
function (𝐳::numModel)(t::StepRangeLen)
  return 𝐳.(collect(t))
end
function (𝐳::numModel)(t::UnitRange)
  return 𝐳.(collect(t))
end

# DISPLAY
Base.show(io::IO, x::numModel) = print(io, "numerical AM-FM model with $(length(x.𝚿ₖ)) components")


text:
# ---
# title: 410. Split Array Largest Sum
# id: problem410
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Binary Search, Dynamic Programming
# link: <https://leetcode.com/problems/split-array-largest-sum/description/>
# hidden: true
# ---
#
# Given an array `nums` which consists of non-negative integers and an integer
# `m`, you can split the array into `m` non-empty continuous subarrays.
#
# Write an algorithm to minimize the largest sum among these `m` subarrays.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [7,2,5,10,8], m = 2
#     Output: 18
#     Explanation:
#     There are four ways to split nums into two subarrays.
#     The best way is to split it into [7,2,5] and [10,8],
#     where the largest sum among the two subarrays is only 18.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3,4,5], m = 2
#     Output: 9
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,4,4], m = 3
#     Output: 4
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 1000`
#   * `0 <= nums[i] <= 106`
#   * `1 <= m <= min(50, nums.length)`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
function just_print(value::T) where {T}
    println(value)
end

function main()
    just_print(100::Int)
    just_print("foo"::String)
end

main()

max_stars_repo_name: myyrakle/Julia-Example

max_issues_repo_name: myyrakle/Julia-Example

max_forks_repo_name: myyrakle/Julia-Example


text:
#!/usr/bin/env julia
#Start Test Script
using ChemometricsTools
using Test
#Pkg.test("ChemometricsTools")
#Pkg.test()

@testset "Transformations" begin
    simplearray = [[1,2,3] [1,2,3]];
    Xform = Center( simplearray )
    @test all( Xform( simplearray ) .== [ [-1,0,1] [-1,0,1] ] )
end

@testset "Pipelines" begin
    #Test a longish pipeline
    FauxData1 = randn(5,10);
    PreprocessPipe = Pipeline(FauxData1, RangeNorm, Center, Scale, RangeNorm);
    Processed = PreprocessPipe(FauxData1);
    @test RMSE( FauxData1, PreprocessPipe(Processed; inverse = true) ) < 1e-14

    #Test inplace pipeline
    OriginalCopy = copy(FauxData1);
    InPlacePipe = PipelineInPlace(FauxData1, RangeNorm, Center, Scale, RangeNorm);
    @test FauxData1 != OriginalCopy
    @test Processed == FauxData1

    #Inplace transform the data back
    InPlacePipe(FauxData1; inverse = true)
    @test RMSE( OriginalCopy, FauxData1 ) < 1e-14

    #Ensure that centerscale center then scales
    Pipe1 = Pipeline(FauxData1, Center, Scale);
    Pipe2 = Pipeline(FauxData1, CenterScale);
    @test RMSE( Pipe1(FauxData1), Pipe2(FauxData1) ) < 1e-14

    #Logit can Inf out pretty easily so test with a small vector
    FauxData2 = [1,1,2,3,4,5,6,7] ./ 10.0;
    Pipe1 = Pipeline(FauxData2,  Logit);
    @test RMSE( FauxData2, Pipe1(Pipe1(FauxData2); inverse = true) ) < 1e-14

    #Test Quantile scaler
    # Pipe1 = Pipeline(FauxData1, x -> QuantileTrim(x; quantiles = (0.25,0.75)), RangeNorm);
    # Pipe1(FauxData1)
end

@testset "Find Peaks" begin
    y = sin.( collect(1:720) .* (pi/180) );
    @test findpeaks(y) == [90, 450]
    @test y[findpeaks(y)] == [1.0, 1.0]
end

@testset "Classification Metrics" begin
    CE = [ :a,:b,:c,:a,:b,:c ]
    @test IsColdEncoded( CE ) == true
    LEnc = LabelEncoding( CE )
    HOT = ColdToHot(CE, LEnc)
    @test HOT == [[1,0,0,1,0,0] [0,1,0,0,1,0] [0,0,1,0,0,1]]
    COLD = HotToCold([[1,0,0,1,0,0] [0,1,0,0,1,0] [0,0,1,0,0,1]], LEnc)
    A = MulticlassStats(COLD, CE, LEnc)
    @test A[1][:Global]["Accuracy"] == 1.0
    @test A == MulticlassStats(HOT, CE, LEnc)
    @test A == MulticlassStats(HOT, HOT, LEnc)
    @test A == MulticlassStats(CE, HOT, LEnc)
end

@testset "Distance Measures" begin
    a = [[1.0,2.0,3.0] [1.0,2.0,3.0]]
    @test all(EuclideanDistance(a) .== [[0,sqrt(2),sqrt(8)] [sqrt(2),0,sqrt(2)] [sqrt(8),sqrt(2),0]])
    @test all(ManhattanDistance(a) .== [[0.0,2.0,4.0] [2.0,0.0,2.0] [4.0,2.0,0.0]] )
end

@testset "FastNNLS Test" begin
    #Test from Paper
    a = reshape( [73,87,72,80, 71,74, 2,89,52,46,7 , 71], 4,3)
    b = [49,67, 68,20];
    c = FNNLS(a, b)
    @test all( [0.649, -1e-6, -1e-6] .< c )
    @test all( c .< [0.651,1e-6,1e-6] )
end

@testset "MonotonicRegression Test" begin
    #Compared with isoreg from R
    a = MonotoneRegression( [1,0,4,3,3,5,4,2,0] .+ 0.0)
    b = MonotoneRegression( [1,0,4,3,3,5,4,2,3] .+ 0.0)
    c = MonotoneRegression( [1,0,0,0,1] .+ 0.0)
    @test sum(a - [0.5, 0.5, 3, 3, 3, 3, 3, 3, 3]) < 1e-8
    @test sum(b - [0.5, 0.5, 3.33333333333, 3.3333333333, 3.333333333333, 3.5, 3.5, 3.5, 3.5]) < 1e-8
    @test all(c .== [0.25,0.25,0.25,0.25,1.0]) || all(c .== [1.0,0.25,0.25,0.25,0.25])
end


text:
function main()

  # list to be sorted
  rng = MersenneTwister(1234);
  A = shuffle(rng, Vector(1:100))
  Aorigin = copy(A)

  # sorting
  @time Aqs = sort(A; alg=QuickSort)
  @time Ais = sort(A; alg=InsertionSort)
  @time Ams = sort(A; alg=MergeSort)
  @time Apqs = sort(A; alg=PartialQuickSort(5))

  print("The original list: ")
  println(Aorigin)

  print("The sorted list: ")
  println(Aqs)

end

if contains(@__FILE__, PROGRAM_FILE)
    main()
end

max_stars_repo_name: szmlb/algorithmsInJulia

max_issues_repo_name: szmlb/algorithmsInJulia

max_forks_repo_name: szmlb/algorithmsInJulia


text:
"""
`module  JAC.READI`  ... a submodel of JAC that contains all methods for computing resonant-excitation Auger double-ionization cross
                         sections and rates; it is using JAC, JAC.ManyElectron, JAC.DoubleAuger.
"""
module READI

    using JAC, JAC.ManyElectron, JAC.DoubleAuger

    """
    `struct  READI.Settings`  ... defines a type for the details and parameters in computing resonant-excitation Auger double-ionization
                                  pathways |i(N)>  --> |m(N+1)>  --> |n(N-1)>.

        + electronEnergies   ::Array{Float64,1}                   ... List of impact-energies of the incoming elecgtrons.
        + selectPathways     ::Bool                               ... True if particular pathways are selected for the computations.
        + selectedPathways   ::Array{Tuple{Int64,Int64,Int64},1}  ... List of list of pathways, given by tupels (inital, inmediate, final).
        + maxKappa           ::Int64                              ... Maximum kappa value of partial waves to be included.
    """
    struct Settings
        electronEnergies     ::Array{Float64,1}
        selectPathways       ::Bool
        selectedPathways     ::Array{Tuple{Int64,Int64,Int64},1}
        maxKappa             ::Int64
    end

    """
    `JAC.READI.Settings()`  ... constructor for the default values of resonant-excitation Auger double-ionization settings.
    """
    function Settings()
       Settings( Float64[], false, Tuple{Int64,Int64,Int64}[], 0)
    end

    # `Base.show(io::IO, settings::READI.Settings)`  ... prepares a proper printout of the variable settings::READI.Settings.
    function Base.show(io::IO, settings::READI.Settings)
        println(io, "electronEnergies:     $(settings.electronEnergies)  ")
        println(io, "selectPathways:       $(settings.selectPathways)  ")
        println(io, "selectedPathways:     $(settings.selectedPathways)  ")
        println(io, "maxKappa:             $(settings.maxKappa)  ")
    end

    """
    `struct  READI.Channel`  ... defines a type for a resonant-excitation Auger double-ionization channel that specifies
                                 all quantum numbers, phases and amplitudes.

        + excitationChannel  ::JAC.ImpactExcitation.Channel      ... Channel that describes the electron-impact excitation process.
        + augerChannel1      ::JAC.AutoIonization.Channel        ... Channel that describes the first subsequent Auger/autoionization process.
        + augerChannel2      ::JAC.AutoIonization.Channel        ... Channel that describes the second Auger/autoionization process.
    """
    struct  Channel
        excitationChannel    ::JAC.ImpactExcitation.Channel
        channel              ::JAC.DoubleAuger.Channel
    end

    """
    `struct  READI.Pathway`  ... defines a type for a resonant-excitation Auger double-ionization pathway that may include
                                 the definition of different excitation and double Auger channels and their corresponding amplitudes.

        + initialLevel           ::Level         ... initial-(state, N-electron) level
        + intermediateLevel      ::Level         ... intermediate-(state, N+1 electron) level m
        + finalLevel             ::Level         ... final-(state,  N-1 electron) level
        + electronInEnergy       ::Float64       ... energy of the (incoming) electron
        + electronOutEnergy      ::Float64       ... energy of the (outgoing, scattered) electron
        + electronAugerEnergy1   ::Float64       ... energy of the (first emitted Auger) electron
        + electronAugerEnergy2   ::Float64       ... energy of the (second emitted Auger) electron
        + crossSection           ::Float64       ... total cross section of this pathway
        + hasChannels            ::Bool          ... Determines whether the indiv. excitation and double autoionization channels are defined in
                                                     terms of their free-electron kappa's, phases and the total angular momentum/parity as well
                                                     as the amplitude, or not.
        + channels               ::Array{READI.Channel,1} ... List of channels of this pathway.
    """
    struct  Pathway
        initialLevel             ::Level
        intermediateLevel        ::Level
        finalLevel               ::Level
        electronInEnergy         ::Float64
        electronOutEnergy        ::Float64
        electronAugerEnergy1     ::Float64
        electronAugerEnergy2     ::Float64
        crossSection             ::Float64
        hasChannels              ::Bool
        channels                 ::Array{READI.Channel,1}
    end

    """
    `JAC.READI.Pathway()`  ... constructor for an electron-impact excitation Auger double ionization pathway between a specified
                               initial, intermediate and final level.
    """
    function Pathway()
        Pathway(Level(), Level(), Level(), 0., 0., 0., 0., 0., false, READI.Channel[] )
    end

    # `Base.show(io::IO, pathway::READI.Pathway)`  ... prepares a proper printout of the variable pathway::READI.Pathway.
    function Base.show(io::IO, pathway::READI.Pathway)
        println(io, "initialLevel:               $(pathway.initialLevel)  ")
        println(io, "intermediateLevel:          $(pathway.intermediateLevel)  ")
        println(io, "finalLevel:                 $(pathway.finalLevel)  ")
        println(io, "electronInEnergy:           $(pathway.electronInEnergy)  ")
        println(io, "electronOutEnergy:          $(pathway.electronOutEnergy)  ")
        println(io, "electronAugerEnergy1:       $(pathway.electronAugerEnergy1)  ")
        println(io, "electronAugerEnergy2:       $(pathway.electronAugerEnergy2)  ")
        println(io, "crossSection:               $(pathway.crossSection)  ")
        println(io, "hasChannels:                $(pathway.hasChannels)  ")
        println(io, "channels:                   $(pathway.channels)  ")
    end

end # module


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

@doc raw"""JSONSchemaPropsOrArray represents a value that can either be a JSONSchemaProps or an array of JSONSchemaProps. Mainly here for serialization purposes.

    IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray(;
    )

"""
mutable struct IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray <: SwaggerModel

    function IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray(;)
        o = new()
        o
    end
end # type IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray

const _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray = Dict{Symbol,Symbol}()
const _property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray = Dict{Symbol,String}()
Base.propertynames(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray }) = collect(keys(_property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray))
Swagger.property_type(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray[name]))}
Swagger.field_name(::Type{ IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray }, property_name::Symbol) =  _property_map_IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray[property_name]

function check_required(o::IoK8sApiextensionsApiserverPkgApisApiextensionsV1beta1JSONSchemaPropsOrArray)
    true
end

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


text:
export ∧

Base.adjoint(f::Function) = (a...) -> (b...) -> f(a..., b...)

supertype_all(xs::T) where T<:AbstractArray = typeof.(xs) |> foldr'(promote_type)

supertype_all(xs...) = typeof.(xs) |> foldr'(promote_type)

import Base.|, Base.*

Base.:(*)(t1::Union{Type, TypeVar}, t2::Union{Type, TypeVar}) = Tuple{t1, t2}
Base.:(|)(t1::Union{Type, TypeVar}, t2::Union{Type, TypeVar}) = Union{t1, t2}

Base.:(*)(p1::Pair{Symbol, T}, p2::Pair{Symbol, S}) where {T<:DataType, S<:DataType} = NamedTuple{(p1.first, p2.first), p1.second * p2.second}

∧(f::Function, i::Int) = i == 1 ? f : x -> (f ∧ (i-1))(f(x))

abstract type AbstractPiecewise{N} <:Function end

function evalpiecewise(::Val{N}, functions, breakpoints, x) where N
    if @generated
        generator = (:(if x < breakpoints[$k]
            return functions[$k](x)
        end) for k in 1:N)
        quote
            @inbounds begin
                $(generator...)
                return functions[$(N+1)](x)
            end
        end
    else
        ind = searchsortedfirst(breakpoints, x)
        return functions[ind](x)
    end
end

function printpiecewise(io::IO, expressions, breakpoints)
    N = length(breakpoints)
    conditions = Vector{String}(undef, N + 1)
    for i in 1:N+1
        if i == 1
            conditions[i] = "x < $(breakpoints[1])"
        elseif i == N+1
            conditions[i] = "x >= $(breakpoints[N])"
        else
            conditions[i] = "$(breakpoints[i-1]) <= x < $(breakpoints[i])"
        end
    end
    maxl = maximum(length.(conditions))
    conditions = rpad.(conditions, maxl + 1)

    for (i, e) in enumerate(expressions)
        println(io, "\t $(conditions[i]):   " * e)
    end
end


text:
using LEAPMacro

curr_working_dir = pwd()
cd(@__DIR__)

println("Running Baseline...")
LEAPMacro.run()

cd(curr_working_dir)

max_stars_repo_name: sei-international/LEAPMacro

max_issues_repo_name: sei-international/LEAPMacro

max_forks_repo_name: sei-international/LEAPMacro


text:
"""
setnexproperty!(nid::NexID, property::Symbol, val)

Set property of a Nextion object given by it's NexID

# Example:
`setnexproperty!(NexID("t0"), :txt, "Hello!")` send and execute
on the Nextion display `t0.txt="Hello"`

`setnexproperty!(NexID("n0"), :val, 3)` send and execute
on the Nextion display `n0.val=3`
"""
function setnexproperty!(nid::NexID, property::Symbol, val::String)
    _name = String(Name(nid))
    _property = String(property)
    cmd = "$_name.$_property=\"$val\""
    nexSerial = NexSerial(nid)
    send(nexSerial, cmd)
end

function setnexproperty!(nid::NexID, property::Symbol, val::Bool)
    if val
        val = 1
    else
        val = 0
    end
    setnexproperty!(nid, property, val)
end

function setnexproperty!(nid::NexID, property::Symbol, val)
    _name = String(Name(nid))
    _property = String(property)
    cmd = "$_name.$_property=$val"
    nexSerial = NexSerial(nid)
    send(nexSerial, cmd)
end

"""
getnexproperty(nid::NexID, property::Symbol)

Get property of a Nextion object given by it's NexID
"""
function getnexproperty(nid::NexID, property::Symbol, ::Type{T}) where {T<:Unsigned}
    # send cmd to Nextion
    _name = String(Name(nid))
    _property = String(property)
    cmd = "get $_name.$_property"
    @info cmd
    nexSerial = NexSerial(nid)
    send(nexSerial, cmd)

    # wait cmd completed
    timeout_ms = nexSerial.timeout_ms
    sp = nexSerial._serial
    waitcommandcomplete(sp, timeout_ms)

    # receive Nextion message
    #r = my_readuntil(sp, v_uint8_eoc, timeout_ms)
    nb, r = sp_blocking_read(sp.ref, my_bytesavailable(sp), timeout_ms)
    @info r
    evt = Event.NumberHeadEvent(r)
    evt.value
end

function getnexproperty(nid::NexID, property::Symbol, ::Type{T}) where {T<:Signed}
    value = getnexproperty(nid, property, Unsigned)
    reinterpret(Int32, value)  # two's complement for UInt32 to Int32
end

function getnexproperty(nid::NexID, property::Symbol, ::Type{String})
    # send cmd to Nextion
    _name = String(Name(nid))
    _property = String(property)
    cmd = "get $_name.$_property"
    @info cmd
    nexSerial = NexSerial(nid)
    send(nexSerial, cmd)

    # wait cmd completed
    timeout_ms = nexSerial.timeout_ms
    sp = nexSerial._serial
    waitcommandcomplete(sp, timeout_ms)

    # receive Nextion message
    #r = my_readuntil(sp, v_uint8_eoc, timeout_ms)
    nb, r = sp_blocking_read(sp.ref, my_bytesavailable(sp), timeout_ms)
    @info r
    evt = Event.StringHeadEvent(r)
    evt.value
end


text:
struct RevRotation end

export RevRotation

"""
    rotate!(::RevRotation, A::AbstractVector, tailpos::Integer)

Rotates `A` on `tailpos` using ``O(1)`` extra memory, yet using several 2n copy-write ops (triple reversing)
"""
function rotate!(::RevRotation, A::AbstractVector, tailpos::Integer)
    reverse!(A, 1, tailpos - 1)
    reverse!(A, tailpos, length(A))
    reverse!(A)
    A
end


text:
Base.Multimedia.showable(::MIME"image/svg+xml", w::MathLink.WExpr) = w.head == W"Graphics"
Base.show(io::IO, ::MIME"image/svg+xml", w::MathLink.WExpr) = print(io, weval(W"ExportString"(w, "SVG")))

Base.Multimedia.showable(::MIME"image/png", w::MathLink.WExpr) = w.head == W"Graphics"
function Base.show(io::IO, ::MIME"image/png", w::MathLink.WExpr)
    x = weval(W"ExportByteArray"(w, "PNG"))
    write(io, x.args[1].args)
end


text:
using OrdinaryDiffEq, DiffEqDevTools, Test, Random
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
using DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_2Dlinear, prob_ode_linear,
       prob_ode_bigfloatlinear, prob_ode_bigfloat2Dlinear

probArr = Vector{ODEProblem}(undef, 2)
bigprobArr = Vector{ODEProblem}(undef, 2)

probArr[1] = prob_ode_linear
probArr[2] = prob_ode_2Dlinear
bigprobArr[1] = prob_ode_bigfloatlinear
bigprobArr[2] = prob_ode_bigfloat2Dlinear
setprecision(400)
Random.seed!(100)
## Convergence Testing
println("Convergence Test on Linear")
dts = 1 .//2 .^(8:-1:4)
testTol = 0.3
superduperbool = Vector{Bool}(undef, 2)

for constructfun in filter(x->startswith(string(x), "construct"), names(DiffEqDevTools))
  tab = getproperty(DiffEqDevTools, constructfun)(BigFloat)
  if tab.order < 12
    if constructfun in (
                        :constructTsitouras9,  # order 1 ???!!!
                        :constructTsitouras92, # embedded order 2 ???!!!
                       )
      @info "Failed $constructfun..."
      @test_broken check_tableau(tab)
    else
      @info "Testing $constructfun..."
      @test check_tableau(tab)
    end
  end
end

# numerically test RKs with order >= 12
for i = 1:2 # 1 = num, 2 = ExplicitRK
  global dts
  if i>1
    prob = probArr[2]
    bigprob = bigprobArr[2]
  else
    prob = probArr[1]
    bigprob = bigprobArr[1]
  end

  for tab in [constructFeagin12(BigFloat), constructFeagin14(BigFloat), constructOno12(BigFloat)]
    @info "Testing..."
    tabalg = ExplicitRK(tableau=tab)
    sim = test_convergence(dts,bigprob,tabalg)
    @test sim.𝒪est[:l∞] >= tab.order || abs(sim.𝒪est[:l∞] - tab.order) < testTol
  end
end


text:
cd(raw"/Users/wapisani/Documents/Programming/Julia/julia_advent_of_code/2015")
open("AoC_2015_7.input") do f
    global instructions = split(read(f,String),"\n")
end

# Iterate over instruction list and add wires to an empty dictionary
global signal_dict = Dict()
global signal_topology = [] # Array of tuples ([inputs], output, type)
for instruction in instructions
    instruction = String(instruction)
    inputs_and_type = split(instruction,"->")[1]
    output = split(instruction,"->")[end]
    signal_dict[output] = Nothing # Get output wire and include in dictionary
    println(split(inputs_and_type," "))

end


text:
__precompile__(true)
module LightGraphsMatching
using LightGraphs
export MatchingResult, maximum_weight_matching, maximum_weight_maximal_matching, minimum_weight_perfect_matching

"""
    type MatchingResult{T}
        weight::T
        mate::Vector{Int}
    end

A type representing the result of a matching algorithm.

    weight: total weight of the matching

    mate:    `mate[i] = j` if vertex `i` is matched to vertex `j`.
             `mate[i] = -1` for unmatched vertices.
"""
struct MatchingResult{T<:Real}
    weight::T
    mate::Vector{Int}
end

import BlossomV
include("blossomv.jl")

using JuMP
using MathProgBase: AbstractMathProgSolver
include("lp.jl")
include("maximum_weight_matching.jl")

end # module


text:
using MatpowerCases, Test

println("- extract_gen check")
for c in casenames()
    mpc = loadcase(c; describe=false)
    @assert isa(extract_gen(mpc), MatpowerCases.Gen)
    @assert isa(extract_gen(c), MatpowerCases.Gen)
end


text:
using Pkg
#  If NetABM is not available then it should be installed via
#  Uncomment if not installed
] add https://github.com/loopulos/NetABM.jl#develop
#  ] add https://github.com/ollin18/NetABM.jl

function useit(list::Array{Symbol})
    installed = [key for key in keys(Pkg.installed())]
    strpackages = @. string(list)
    uninstalled = setdiff(strpackages,installed)

    map(Pkg.add,uninstalled)
    for package ∈ list
        @eval using $package
    end
end

thep = [:LightGraphs, :StatsBase, :Plots, :GraphPlot, :Statistics, :LaTeXStrings, :Bootstrap, :NetABM]
useit(thep)

function lectura_uw(red)
    Nodes = union(unique(red[:,1]),unique(red[:,2]))
    g = SimpleGraph()
    last_node = Int64(length(Nodes))
    add_vertices!(g,last_node)
    for n in 1:size(red)[1]
        add_edge!(g,red[n,1],red[n,2])
        add_edge!(g,red[n,2],red[n,1])
    end
    return g, Nodes
end

function get_network(;N = 100, k = 20, maxk = 40, mu = 0.1, t1 = 2, t2 = 1, minc = 10, maxc = 40, c=0.3)
        the_root = pwd()
        try
            run(`$the_root/src/binary_networks/benchmark -N $N -k $k -maxk $maxk -mu $mu -C $c`)
            #  run(`$the_root/src/binary_networks/benchmark -N $N -k $k -maxk $maxk -mu $mu`)
        catch
            N
        end
end
the_root = pwd()

get_network(N = 100, mu=0.1);

################################################################
################################################################
################################################################
################################################################
################################################################
# The same as before but instead of looking at neighbor infections
# we look use a mean cooperation effect by distance

global num_agents = 100
global n_boot = 1000
global cil = 0.95
global j
#  global muladol = [(0,[0.1]), (1,repeat([0.1],2)),(2,repeat([0.1],3)),(3,repeat([0.1],4))]
#  global muladol = [(i,repeat([0.6],i+1)) for i in 0:10]
#  global muladol = [(i,repeat([0.1],i+1)) for i in 0:10]
global muladol = [(i,repeat([.1],i)) for i in 1:6]
global muladol = [(i,LinRange(0.01,0.001,6)|>collect) for i in 1:6]

global adoptstep = (1,[0.1, 0.1])
global maxI = Array{Float64}(undef,0)
global sumI = []

@time let
    #  for adoptstep in muladol
    #  for numnei in 6:4:26
    #  for numnei in Int.(floor.(LinRange(4,80,6)))
    for numnei in Int.(ceil.(exp.(range(1,stop=4,length=11))))
    #  for mixing in [0.9, 0.1]#0.1:0.1:0.6
    for copera in 0.0:0.1:1
        for adap in 0.0:0.1:1
            global meantotI = Array{Float64}(undef,0)
            global Imean = Array{Float64}(undef,0)
            global Ilow = Array{Float64}(undef,0)
            global Ihigh = Array{Float64}(undef,0)
            global Itotall = Array{Float64}(undef,0)
            global Iall = Array{Array}(undef,0)
            global Sall = Array{Array}(undef,0)
            global Rall = Array{Array}(undef,0)
            for a in 1:15
                g = barabasi_albert(num_agents,numnei)
                        #  get_network(N = num_agents, mu=0.5, c=mixing);
                        #  the_net = readdlm(the_root*"/network.dat")
                        #  g, Nodes = lectura_uw(the_net)
                agents = [Agent(i) for i = 1:num_agents];
                init_demographics!(agents;states=["S","I"],initial=[1.0,0.0]);
                pat_0 = rand(1:nv(g))
                agents[pat_0].state = "I"
                #  init_demographics!(agents;states=["S","I"],initial=[0.9,0.1]);
                set_coop_agents!(agents;p_cop=copera);
                set_adapt_agents!(agents;p_cop=adap);
                set_coop_effect!(agents)
                map(x -> assign_contacts!(g,x), agents);
                num_S = []
                num_I = []
                num_R = []
                global oldS = ([ag.state for ag in agents if ag.state=="S"] |> length)/num_agents
                global oldI = ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents
                global totI = ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents
                for i in 1:100
                    push!(num_S, ([ag.state for ag in agents if ag.state=="S"] |> length)/num_agents)
                    push!(num_I, ([ag.state for ag in agents if ag.state=="I"] |> length)/num_agents)
                    push!(num_R, ([ag.state for ag in agents if ag.state=="R"] |> length)/num_agents)

                    next_state!(agents;fun=SI_attitude!,inf_prob=0.1, rec_prob=0.03, R=true);
                    update_state!(agents);
                    update_effect_given_distance_coop!(agents,g,adoptstep[1],1,adoptstep[2]);

                    newI = oldS - num_S[end]
                    global totI = totI + newI
                    oldS = num_S[end]

                end
                push!(Iall,num_I)
                push!(Itotall,totI)
                GC.gc()
            end
            for step in eachindex(Iall[1])
                bs1 = bootstrap(mean,getindex.(Iall,step),BasicSampling(n_boot))
                bs2 = bootstrap(std,getindex.(Iall,step),BasicSampling(n_boot))
                bci1 = confint(bs1, BasicConfInt(cil))
                bci2 = confint(bs2, BasicConfInt(cil))
                push!(Imean,bci1[1][1])
                push!(Ilow,bci2[1][2])
                push!(Ihigh,bci2[1][3])
            end
            bs1 = bootstrap(mean, Itotall,BasicSampling(n_boot))
            bci1 = confint(bs1, BasicConfInt(cil))
            push!(sumI,bci1[1][1])
            push!(maxI,maximum(Imean))
        end
        GC.gc()
    end
    end
end

theme(:default)
using DelimitedFiles
numnei = Int.(ceil.(exp.(range(1,stop=4,length=11))))

plt = heatmap(
        xlabel = "Compliant Agents",
        ylabel = "Adaptive Agents",
        title = "Maximum Infections",
        guidefontsize = 16,
        tickfontsize = 10,
        xlims = (0,11)
    )
xs = [string(i) for i = 0:0.1:1]
ys = [string(i) for i = 0:0.1:1]
z = reshape(maxI, (11, 11, 11))
    z1 = z[:,:,1]
mean(z1)
plt = heatmap!(xs, ys, z1, aspect_ratio = 1, clim=(0,1), c = :algae,dpi=300)
png(plt,"figs/heatmap_scaled_maximum_coop.png")
savefig(plt,"figs/heatmap_scaled_maximum_coop.pdf")
writedlm("data/max_inf_coop_neigh_num.csv",maxI)
writedlm("data/max_inf_coop_neigh_exp.csv",maxI)

for i in 1:11
    j = numnei[i]
    z1 = z[:,:,i]
    #  plt = heatmap!(xs, ys, z1, title=L"\langle k\rangle =$i", aspect_ratio = 1, clim=(0,1), c = :algae,dpi=300)
    plt = heatmap!(xs, ys, z1, aspect_ratio = 1, clim=(0,1), c = :algae,dpi=300)
    png(plt,"figs/heatmap_scaled_maximum_coop_neigh_$i.png")
    savefig(plt,"figs/heatmap_scaled_maximum_coop_neigh_$i.pdf")
end

plt = heatmap(
        xlabel = "Compliant Agents",
        ylabel = "Adaptive Agents",
        title = "Maximum Infections",
        guidefontsize = 16,
        tickfontsize = 10,
        xlims = (0,11)
    )
xs = [string(i) for i = 0:0.1:1]
ys = [string(i) for i = 0:0.1:1]
z = reshape(sumI, (11, 11, 11))
    z1 = z[:,:,1]
mean(z1)
plt = heatmap!(xs, ys, z1, aspect_ratio = 1, clim=(0,1), c = :algae,dpi=300)
png(plt,"figs/heatmap_scaled_maximum_coop.png")
savefig(plt,"figs/heatmap_scaled_maximum_coop.pdf")
writedlm("data/max_inf_coop_neigh_num.csv",maxI)
writedlm("data/max_inf_coop_neigh_exp.csv",maxI)


text:
@testset "Testing functions on polygons and shapefiles" begin
package_path = pathof(SatelliteImages)
path_len = length(package_path)
assets_path = package_path[1:path_len-22] * "assets"
map_path = assets_path * "/mumbai_map/mumbai_districts.shp"

radiance_datacube = rand(1:100, 156, 85, 95)
mumbai_districts = load_shapefile(map_path)
mumbai_districts_ntl = aggregate_dataframe(MUMBAI_COORDINATE_SYSTEM, radiance_datacube, mumbai_districts, "DISTRICT")
@test sizeof(mumbai_districts_ntl) > 0

district1 = mumbai_districts[1,:] # Select the first district
district1_mask = polygon_mask(MUMBAI_COORDINATE_SYSTEM, district1)

@test sum(district1_mask) >=0

end


text:
# Kostelec & Rockmore

function wigner_d̃{T<:AbstractFloat}(J::Int, M::Int, M′::Int, β::T)
    return if J < max(abs(M), abs(M′))
        zero(T)
    elseif J == 0
        sqrt(1/2)
    elseif J == 1 && M == 0 && M′ == 0
        sqrt(3/2) * cos(β)
    else
        if J == max(abs(M), abs(M′))
            p = sqrt((2J + 1) / 2)

            if M == J
                p *
                sqrt(factorial(2J) / (factorial(J + M′) * factorial(J - M′))) *
                cos(β/2)^(J + M′) * (-sin(β/2))^(J - M′)
            elseif M == -J
                p *
                sqrt(factorial(2J) / (factorial(J + M′) * factorial(J - M′))) *
                cos(β/2)^(J - M′) * ( sin(β/2))^(J + M′)
            elseif M′ == J
                p *
                sqrt(factorial(2J) / (factorial(J + M ) * factorial(J - M ))) *
                cos(β/2)^(J + M ) * ( sin(β/2))^(J - M )
            else #if M′== -J
                p *
                sqrt(factorial(2J) / (factorial(J + M ) * factorial(J - M ))) *
                cos(β/2)^(J - M ) * (-sin(β/2))^(J + M )
            end
        else
            sqrt((2J + 1) / (2J - 1)) *
            ((J * (2J - 1)) / sqrt((J^2 - M^2) * (J^2 - M′^2))) *
            (cos(β) - (M*M′) / ((J - 1) * J)) *
            wigner_d̃(J - 1, M, M′, β) -
            sqrt((2J + 1) / (2J - 3)) *
            (sqrt(((J - 1)^2 - M^2) * ((J - 1)^2 - M′^2)) / sqrt((J^2 - M^2) * (J^2 - M′^2))) *
            (J / (J - 1)) *
            wigner_d̃(J - 2, M, M′, β)
        end
    end
end

wigner_D̃{T<:AbstractFloat}(J::Int, M::Int, M′::Int, α::T, β::T, γ::T) = exp(-im * M * α) * wigner_d̃(J, M, M′, β) * exp(-im * M′ * γ)

wigner_d{T<:AbstractFloat}(J::Int, M::Int, M′::Int, β::T) = sqrt(2 / (2J + 1)) * wigner_d̃(J, M, M′, β)
wigner_D{T<:AbstractFloat}(J::Int, M::Int, M′::Int, α::T, β::T, γ::T) = sqrt(2 / (2J + 1)) * wigner_D̃(J, M, M′, α, β, γ)

wigner_d̃_matrix{T<:AbstractFloat}(J::Int, β::T) = [wigner_d̃(J, M, M′, β) for M in -J:J, M′ in -J:J]
wigner_D̃_matrix{T<:AbstractFloat}(J::Int, α::T, β::T, γ::T) = [wigner_D̃(J, M, M′, α, β, γ) for M in -J:J, M′ in -J:J]

wigner_d_matrix{T<:AbstractFloat}(J::Int, β::T) = [wigner_d(J, M, M′, β) for M in -J:J, M′ in -J:J]
wigner_D_matrix{T<:AbstractFloat}(J::Int, α::T, β::T, γ::T) = [wigner_D(J, M, M′, α, β, γ) for M in -J:J, M′ in -J:J]

irrep(g::O3, L::Int) = wigner_D_matrix(L, g.α, g.β, g.γ)


text:
module TestUtilities
using Test
import IOCapture

export @quietly

const QUIETLY_LOG = joinpath(@__DIR__, "quietly.log")

__init__() = isfile(QUIETLY_LOG) && rm(QUIETLY_LOG)

struct QuietlyException <: Exception
    exception
    backtrace
end

function Base.showerror(io::IO, e::QuietlyException)
    println(io, "@quietly hit an exception ($(typeof(e.exception))):")
    showerror(io, e.exception, e.backtrace)
end

function _quietly(f, expr, source)
    c = IOCapture.iocapture(f; throwerrors = :interrupt)
    haskey(ENV, "DOCUMENTER_TEST_QUIETLY") && open(QUIETLY_LOG; write=true, append=true) do io
        println(io, "@quietly: c.error = $(c.error) / $(sizeof(c.output)) bytes of output captured")
        println(io, "@quietly: $(source.file):$(source.line)")
        println(io, "@quietly: typeof(result) = ", typeof(c.value))
        println(io, "@quietly: STDOUT")
        println(io, output)
        println(io, "@quietly: end of STDOUT")
        if success
            println(io, "@quietly: result =")
            println(io, c.value)
        else
            println(io, "@quietly: result (error) =")
            showerror(io, c.value, c.backtrace)
        end
    end
    if !c.error
        printstyled("@quietly: success, $(sizeof(c.output)) bytes of output hidden\n"; color=:magenta)
        return c.value
    else
        @error """
        An error was thrown in @quietly, $(sizeof(c.output)) bytes of output captured
        $(typeof(c.value)) at $(source.file):$(source.line) in expression:
        $(expr)
        """
        if !isempty(c.output)
            printstyled("$("="^21) @quietly: output from the expression $("="^21)\n"; color=:magenta)
            print(c.output)
            last(c.output) != "\n" && println()
            printstyled("$("="^27) @quietly: end of output $("="^28)\n"; color=:magenta)
        end
        throw(QuietlyException(c.value, c.backtrace))
    end
end
macro quietly(expr)
    orig_expr = Expr(:inert, expr)
    source = QuoteNode(__source__)
    quote
        _quietly($orig_expr, $source) do
            $(esc(expr))
        end
    end
end

end


text:
module NMRDataSetup

import PyCall, Optim, NLopt

include("utils.jl")
include("DSP.jl")
include("load.jl")
include("assemble_mixture.jl")

export loadspectrum, getwraparoundDFTfreqs, evalcomplexLorentzian, gettimerange

end


text:
struct PublishDiagnosticsParams <: Outbound
    uri::DocumentUri
    version::Union{Int,Missing}
    diagnostics::Vector{Diagnostic}
end

import Base.==
==(x::CompletionItem, y::CompletionItem) = x.label == y.label
==(m1::MarkedString, m2::MarkedString) = m1.language == m2.language && m1.value == m2.value

##############################################################################
# Code Action
const CodeActionKind = String
const CodeActionKinds = (Empty = "",
                         QuickFix = "quickfix",
                         Refactor = "refactor",
                         RefactorExtract = "refactor.extract",
                         RefactorInline = "refactor.inline",
                         RefactorRewrite = "refactor.rewrite",
                         Source = "source",
                         SourceOrganizeImports = "source.organiseImports")

@dict_readable struct CodeActionKindCapabilities
    valueSet::Vector{CodeActionKind}
end

@dict_readable struct CodeActionLiteralCapabilities
    codeActionKind::CodeActionKindCapabilities
end

@dict_readable struct CodeActionClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
    codeActionLiteralSupport::Union{CodeActionLiteralCapabilities,Missing}
    isPreferredSupport::Union{Bool,Missing}
end

struct CodeActionOptions <: Outbound
    codeActionKinds::Union{Vector{CodeActionKind},Missing}
    workDoneProgress::Union{Bool,Missing}
end

struct CodeActionRegistrationOptions <: Outbound
    documentSelector::Union{DocumentSelector,Nothing}
    codeActionKinds::Union{Vector{CodeActionKind},Missing}
    workDoneProgress::Union{Bool,Missing}
end

@dict_readable struct CodeActionContext <: Outbound
    diagnostics::Vector{Diagnostic}
    only::Union{Vector{CodeActionKind},Missing}
end

@dict_readable struct CodeActionParams
    textDocument::TextDocumentIdentifier
    range::Range
    context::CodeActionContext
end

struct CodeAction <: Outbound
    title::String
    kind::Union{CodeActionKind,Missing}
    diagnostics::Union{Vector{Diagnostic},Missing}
    isPreferred::Union{Bool,Missing}
    edit::Union{WorkspaceEdit,Missing}
    command::Union{Command,Missing}
end

##############################################################################
# Code Lens
@dict_readable struct CodeLensClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
end

struct CodeLensOptions <: Outbound
    resolveProvider::Union{Bool,Missing}
    workDoneProgress::Union{Bool,Missing}
end

struct CodeLensRegistrationOptions <: Outbound
    documentSelector::Union{DocumentSelector,Nothing}
    resolveProvider::Union{Bool,Missing}
    workDoneProgress::Union{Bool,Missing}
end

@dict_readable struct CodeLensParams
    textDocument::TextDocumentIdentifier
end

struct CodeLens <: Outbound
    range::Range
    command::Union{Command,Missing}
    data::Union{Any,Missing}
end

##############################################################################
# Document Link Provider
@dict_readable struct DocumentLinkClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
    tooltipSupport::Union{Bool,Missing}
end

struct DocumentLinkOptions <: Outbound
    resolveProvider::Union{Bool,Missing}
    workDoneProgress::Union{Bool,Missing}
end

struct DocumentLinkRegistrationOptions <: Outbound
    documentSelector::Union{DocumentSelector,Nothing}
    resolveProvider::Union{Bool,Missing}
    workDoneProgress::Union{Bool,Missing}
end

@dict_readable struct DocumentLinkParams
    textDocument::TextDocumentIdentifier
end

struct DocumentLink <: Outbound
    range::Range
    target::Union{String,Missing}
    tooltip::Union{String,Missing}
    data::Union{Any,Missing}
end

##############################################################################
# Document Colour

@dict_readable struct DocumentColorClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
end

struct DocumentColorOptions <: Outbound
    workDoneProgress::Union{Bool,Missing}
end

struct DocumentColorRegistrationOptions <: Outbound
    documentSelector::Union{DocumentSelector,Nothing}
    id::Union{String,Missing}
    workDoneProgress::Union{Bool,Missing}
end

@dict_readable struct DocumentColorParams <: Outbound
    textDocument::TextDocumentIdentifier
    workDoneToken::Union{ProgressToken,Missing}
    partialResultToken::Union{ProgressToken,Missing}
end

struct Color <: Outbound
    red::Float64
    green::Float64
    blue::Float64
    alpha::Float64
end

struct ColorInformation <: Outbound
    range::Range
    color::Color
end

@dict_readable struct ColorPresentationParams
    textDocument::TextDocumentIdentifier
    color::Color
    range::Range
end

struct ColorPresentaiton <: Outbound
    label::String
    textEdit::Union{TextEdit,Missing}
    additionalTextEdits::Union{Vector{TextEdit},Missing}
end

##############################################################################
# Formatting

##############################################################################
# Rename

@dict_readable struct RenameClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
    prepareSupport::Union{Bool,Missing}
end

struct RenameOptions <: Outbound
    workDoneProgress::Union{Bool,Missing}
    prepareProvider::Union{Bool,Missing}
end

struct RenameRegistrationOptions <: Outbound
    documentSelector::Union{DocumentSelector,Nothing}
    workDoneProgress::Union{Bool,Missing}
    prepareProvider::Union{Bool,Missing}
end

@dict_readable struct RenameParams <: Outbound
    textDocument::TextDocumentIdentifier
    position::Position
    workDoneToken::Union{ProgressToken,Missing}
    newName::String
end

@dict_readable struct PrepareRenameParams
    textDocument::TextDocumentIdentifier
    position::Position
end

##############################################################################
# Folding
const FoldingRangeKind = String
const FoldingRangeKinds = (Comment = "comment",
                           Imports = "imports",
                           Region = "region")

@dict_readable struct FoldingRangeClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
    rangeLimit::Union{Int,Missing}
    lineFoldingOnly::Union{Bool,Missing}
end

struct FoldingRangeOptions <: Outbound
    workDoneProgress::Union{Bool,Missing}
end

struct FoldingRangeRegistrationOptions <: Outbound
    documentSelector::Union{DocumentSelector,Nothing}
    workDoneProgress::Union{Bool,Missing}
    id::Union{String,Missing}
end

@dict_readable struct FoldingRangeParams <: Outbound
    textDocument::TextDocumentIdentifier
    workDoneToken::Union{ProgressToken,Missing}
    partialResultToken::Union{ProgressToken,Missing}
end

struct FoldingRange <: Outbound
    startLine::Int
    startCharacter::Union{Int,Missing}
    endLine::Int
    endCharacter::Union{Int,Missing}
    kind::Union{FoldingRangeKind,Missing}
end

##############################################################################
# Selection Range
@dict_readable struct SelectionRangeClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
end

struct SelectionRangeOptions <: Outbound
    workDoneProgress::Union{Bool,Missing}
end

struct SelectionRangeRegistrationOptions <: Outbound
    workDoneProgress::Union{Bool,Missing}
    documentSelector::Union{DocumentSelector,Nothing}
    id::Union{String,Missing}
end

@dict_readable struct SelectionRangeParams <: Outbound
    workDoneToken::Union{ProgressToken,Missing}
    partialResultToken::Union{ProgressToken,Missing}
    textDocument::TextDocumentIdentifier
    positions::Vector{Position}
end

struct SelectionRange <: Outbound
    range::Range
    parent::Union{SelectionRange,Missing}
end

##############################################################################
# Execute command
@dict_readable struct ExecuteCommandClientCapabilities <: Outbound
    dynamicRegistration::Union{Bool,Missing}
end

struct ExecuteCommandOptions <: Outbound
    workDoneProgress::Union{Bool,Missing}
    commands::Vector{String}
end

mutable struct ExecuteCommandRegistrationOptions <: Outbound
    workDoneProgress::Union{Bool,Missing}
    commands::Vector{String}
end

@dict_readable struct ExecuteCommandParams <: Outbound
    workDoneToken::Union{ProgressToken,Missing}
    command::String
    arguments::Union{Vector{Any},Missing}
end

##############################################################################

struct ApplyWorkspaceEditParams <: Outbound
    label::Union{String,Missing}
    edit::WorkspaceEdit
end

@dict_readable struct ApplyWorkspaceEditResponse <: Outbound
    applied::Bool
    failureReason::Union{String,Missing}
end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb02sd(datfile, io=stdout)
    ZERO = 0.0e0
    NIN = 5
    NOUT = 6
    NMAX = 20
    NRHMAX = 20
    LDB = NMAX
    LDH = NMAX
    LDWORK = 3*NMAX
    LIWORK = NMAX
    H = Array{Float64,2}(undef, LDH,NMAX)
    IPIV = Array{BlasInt,1}(undef, NMAX)
    B = Array{Float64,2}(undef, LDB,NRHMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    IWORK = Array{BlasInt,1}(undef, LIWORK)
f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    NRHS = parse(BlasInt, vs[2])
    NORM =  vs[3][1]
    TRANS =  vs[4][1]
if ( N<0 || N>NMAX )
else

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       H[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
if ( NRHS<0 || NRHS>NRHMAX )
else

    vs = String[]
    _isz,_jsz = (N,NRHS)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       B[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    if  N>2
# interp call 1

    #FOREIGN.dlaset!( 'Lower', N-2, N-2, ZERO, ZERO, H(3,1), LDH )
    triu!(H, -1)

    end
# interp call 2

    INFO = SLICOT.mb02sd!(N, H, IPIV)
    @test INFO == 0
    INFO == 0 || return

    #HNORM = DLANHS( NORM, N, H, LDH, DWORK )
    HNORM = ccall((LinearAlgebra.BLAS.@blasfunc(dlanhs_), LinearAlgebra.LAPACK.liblapack),
                  Float64, (Ref{UInt8}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt},
                            Ptr{Float64}), NORM, N, H, LDH, DWORK)
# interp call 3

    RCOND, INFO = SLICOT.mb02td!(NORM, N, HNORM, H, IPIV)
    println(io, "RCOND = $RCOND")
    @test INFO == 0
    INFO == 0 || return

if ( INFO==0 && RCOND>eps(0.9) )
# interp call 4

    INFO = SLICOT.mb02rd!(TRANS, N, NRHS, H, IPIV, B)
    @test INFO == 0
    INFO == 0 || return

else
end # if
# interp output 1
    println(io, "B:")
    _nc = NRHS
    _nr = N
    show(io, "text/plain", B[1:_nr,1:_nc])
    println(io)

end # if
end # if
    close(f)
end # run_mb02sd()


text:
using Base.Threads

function InitialiseMarkerPositions!( p::Markers, domain::ModelDomain )
    p.nmark     = domain.ncx*domain.ncy*p.nmx*p.nmy; # total initial number of marker in grid
    p.nmark_ini = p.nmark; # total initial number of marker in grid
    p.nmark_max = 3*p.nmark; # total initial number of marker in grid
    dxm, dym    = domain.dx/p.nmx, domain.dy/p.nmy
    xm1d        = LinRange(domain.xmin+dxm/2, domain.xmax-dxm/2, domain.ncx*p.nmx)
    ym1d        = LinRange(domain.ymin+dym/2, domain.ymax-dym/2, domain.ncy*p.nmy)
    (xmi,ymi)   = ([x for x=xm1d,y=ym1d], [y for x=xm1d,y=ym1d])
    p.x         = zeros(Float64, p.nmark_max)
    p.y         = zeros(Float64, p.nmark_max)
    p.x[1:p.nmark] .= vec(xmi)
    p.y[1:p.nmark] .= vec(ymi)
    return
end
export InitialiseMarkerPositions!

@views function CopyMarkers!(p::Markers,ineigh::Int64)
    p.phase[p.nmark] = p.phase[ineigh]
    p.T[p.nmark]     = p.T[ineigh]
end
export CopyMarkers!

@views function FindClosestNeighbour( p::Markers, mpc::Matrix{Float64}, mlist::Matrix{Array{Int64, N} where N}, new_x::Float64, new_y::Float64, dxy::Float64, i::Int64, j::Int64, ncx::Int64, ncy::Int64 )
    ineigh = 0
    dst = dxy
    # Current cell
    for im=1:Int64(mpc[i,j])
        xm = p.x[mlist[i,j][im]]
        ym = p.y[mlist[i,j][im]]
        d = (new_x-xm)^2 + (new_y-ym)^2
        if d<dst
            dst    = d
            ineigh = mlist[i,j][im]
        end
    end
    if ineigh==0
        # Direct neighbours: W
        if i>1 # could make a periodic case here...
            for im=1:Int64(mpc[i-1,j])
                xm = p.x[mlist[i-1,j][im]]
                ym = p.y[mlist[i-1,j][im]]
                d = (new_x-xm)^2 + (new_y-ym)^2
                if d<dst
                    dst    = d
                    ineigh = mlist[i-1,j][im]
                end
            end
        end
        # Direct neighbours: E
        if i<ncx # could make a periodic case here...
            for im=1:Int64(mpc[i+1,j])
                xm = p.x[mlist[i+1,j][im]]
                ym = p.y[mlist[i+1,j][im]]
                d = (new_x-xm)^2 + (new_y-ym)^2
                if d<dst
                    dst    = d
                    ineigh = mlist[i+1,j][im]
                end
            end
        end
        # Direct neighbours: S
        if j>1 # could make a periodic case here...
            for im=1:Int64(mpc[i,j-1])
                xm = p.x[mlist[i,j-1][im]]
                ym = p.y[mlist[i,j-1][im]]
                d = (new_x-xm)^2 + (new_y-ym)^2
                if d<dst
                    dst    = d
                    ineigh = mlist[i,j-1][im]
                end
            end
        end
        # Direct neighbours: S
        if j<ncy # could make a periodic case here...
            for im=1:Int64(mpc[i,j+1])
                xm = p.x[mlist[i,j+1][im]]
                ym = p.y[mlist[i,j+1][im]]
                d = (new_x-xm)^2 + (new_y-ym)^2
                if d<dst
                    dst    = d
                    ineigh = mlist[i,j+1][im]
                end
            end
        end
    end
    # Last check, throw error
    if (ineigh==0) error("Some cell ran empty of markers...") end
    return ineigh
end
export FindClosestNeighbour

@views function ReseedMarkers!( p::Markers, mpc::Matrix{Float64}, xc::LinRange{Float64}, yc::LinRange{Float64}, dx::Float64, dy::Float64, ncx::Int64, ncy::Int64)

    mlist       = Array{Array{Int64}}(undef,ncx,ncy)
    mlist_th    = Array{Array{Int64}}(undef,nthreads(),ncx,ncy)

    @tturbo for i=1:ncx
        for j=1:ncy
            mlist[i,j] = zeros(Int64(mpc[i,j]))
        end
    end

    count = zeros(Int64,ncx,ncy)

    for k=1:p.nmark # @avx ne marche pas ici
        if p.phase[k]>=0
            # Get the column:
            dstx = p.x[k] - xc[1];
            i    = Int64(round(ceil( (dstx/dx) + 0.5)));
            # Get the line:
            dsty = p.y[k] - yc[1];
            j    = Int64(round(ceil( (dsty/dy) + 0.5)));
            # Increment cell count
            count[i,j] += 1
            mlist[i,j][count[i,j]] = k
        end
    end

    # display(mlist)
    dxy = 10*(dx*dx + dy*dy)
    for i=1:ncx
        for j=1:ncy

            if mpc[i,j] < 4

                new_x = xc[i]-dx/4.0
                new_y = yc[j]-dy/4.0
                ineigh = FindClosestNeighbour( p, mpc, mlist, new_x, new_y, dxy, i, j, ncx, ncy )
                if ineigh>0
                    p.nmark+=1
                    p.x[p.nmark] = new_x
                    p.y[p.nmark] = new_y
                    CopyMarkers!(p,ineigh)
                end

                new_x = xc[i]+dx/4.0
                new_y = yc[j]-dy/4.0
                ineigh = FindClosestNeighbour( p, mpc, mlist, new_x, new_y, dxy, i, j, ncx, ncy )
                if ineigh>0
                    p.nmark+=1
                    p.x[p.nmark] = new_x
                    p.y[p.nmark] = new_y
                    CopyMarkers!(p,ineigh)
                end

                new_x = xc[i]+dx/4.0
                new_y = yc[j]+dy/4.0
                ineigh = FindClosestNeighbour( p, mpc, mlist, new_x, new_y, dxy, i, j, ncx, ncy )
                if ineigh>0
                    p.nmark+=1
                    p.x[p.nmark] = new_x
                    p.y[p.nmark] = new_y
                    CopyMarkers!(p,ineigh)
                end

                new_x = xc[i]-dx/4.0
                new_y = yc[j]+dy/4.0
                ineigh = FindClosestNeighbour( p, mpc, mlist, new_x, new_y, dxy, i, j, ncx, ncy )
                if ineigh>0
                    p.nmark+=1
                    p.x[p.nmark] = new_x
                    p.y[p.nmark] = new_y
                    CopyMarkers!(p,ineigh)
                end

            end
        end
    end
end
export ReseedMarkers!

@views function ListMarkers( p::Markers, ncx, ncy )
liste = hcat([[Int[] for i in 1:ncx] for j in 1:ncy]...)
for k in 1:p.nmark
    if p.phase[k]>=0
        i = p.cellx[k]
        j = p.celly[k]
        push!(liste[i,j], k)
    end
end
return liste
end
export ListMarkers

function RungeKutta!(p::Markers, nmark::Int64, rkv::Matrix{Float64}, rkw::Matrix{Float64}, BC::BoundaryConditions, dt::Float64, Vx::Matrix{Float64}, Vy::Matrix{Float64}, xv::LinRange{Float64}, yv::LinRange{Float64}, xce::LinRange{Float64}, yce::LinRange{Float64}, dx::Float64, dy::Float64, ncx::Int64, ncy::Int64)
# Marker advection with 4th order Roger-Kutta
@threads for k=1:nmark
    in = p.phase[k]>=0
        x0 = p.x[k];
        y0 = p.y[k];
        vx = 0.0
        vy = 0.0
        # Roger-Gunther loop
        for rk=1:4
            # println("step", rk)
            # Interp velocity from grid
            vxm = VxFromVxNodes(Vx, k, p, xv, yce, dx, dy, ncx, ncy, 1)
            vym = VyFromVxNodes(Vy, k, p, xce, yv, dx, dy, ncx, ncy, 1)
            # Temporary RK advection steps
            p.x[k] = x0 + rkv[rk]*dt*vxm
            p.y[k] = y0 + rkv[rk]*dt*vym
            if BC.periodix==1
                if p.x[k]<xv[1]
                    # println("periW 1 ", p.x[k], vxm)
                    p.x[k]+=(xv[end]-xv[1])
                    # println("periW 2 ", p.x[k], vxm)
                end
                if p.x[k]>xv[end]
                    # println("periE")
                    p.x[k]-=(xv[end]-xv[1])
                end
            end
            # Average final velocity
            vx    += rkw[rk]*vxm
            vy    += rkw[rk]*vym
        end
        # Advect points
        p.x[k] = x0 + (in==1) * rkv[4]*dt*vx
        p.y[k] = y0 + (in==1) * rkv[4]*dt*vy
        if BC.periodix==1
            if p.x[k]<xv[1]
                # println("periW 1 ", p.x[k], vxm)
                p.x[k]+=(xv[end]-xv[1])
                # println("periW 2 ", p.x[k], vxm)
            end
            if p.x[k]>xv[end]
                # println("periE")
                p.x[k]-=(xv[end]-xv[1])
            end
        end
end
end
export RungeKutta!

# Interpolation from Vx nodes to particles
@views function VxFromVxNodes(Vx::Matrix{Float64}, k::Int64, p::Markers, xv::LinRange{Float64}, yce::LinRange{Float64}, dx::Float64, dy::Float64, ncx::Int64, ncy::Int64, new::Int64)
    # Interpolate vx
    i = Int64(round(trunc( (p.x[k] -  xv[1])/dx ) + 1));
    j = Int64(round(trunc( (p.y[k] - yce[1])/dy ) + 1));
    if i<1
        i = 1
    elseif i>ncx
        i = ncx
    end
    if j<1
        j = 1;
    elseif j> ncy+1
        j = ncy+1
    end
    # Compute distances
    dxmj = p.x[k] -  xv[i]
    dymi = p.y[k] - yce[j]
    # Compute vx velocity for the top and bottom of the cell
    vxm13 = Vx[i,j  ] * (1-dxmj/dx) + Vx[i+1,j  ]*dxmj/dx
    vxm24 = Vx[i,j+1] * (1-dxmj/dx) + Vx[i+1,j+1]*dxmj/dx
    if new==1
        if dxmj/dx>=0.5
            if i<ncx
                vxm13 += 0.5*((dxmj/dx-0.5)^2) * (Vx[i,j  ] - 2.0*Vx[i+1,j  ] + Vx[i+2,j  ]);
                vxm24 += 0.5*((dxmj/dx-0.5)^2) * (Vx[i,j+1] - 2.0*Vx[i+1,j+1] + Vx[i+2,j+1]);
            end
        else
            if i>1
                vxm13 += 0.5*((dxmj/dx-0.5)^2) * (Vx[i-1,j  ] - 2.0*Vx[i,j  ] + Vx[i+1,j  ]);
                vxm24 += 0.5*((dxmj/dx-0.5)^2) * (Vx[i-1,j+1] - 2.0*Vx[i,j+1] + Vx[i+1,j+1]);
            end
        end
    end
    # Compute vx
    vxm = (1-dymi/dy) * vxm13 + (dymi/dy) * vxm24
    # println(p.x[k], " ", xv[k], " ", vxm)
    return vxm
end
export VxFromVxNodes

@views function VyFromVxNodes(Vy::Matrix{Float64}, k::Int64, p::Markers, xce::LinRange{Float64}, yv::LinRange{Float64}, dx::Float64, dy::Float64, ncx::Int64, ncy::Int64, new::Int64)
    # Interpolate vy
    i = Int64(round(trunc( (p.x[k] - xce[1])/dx ) + 1));
    j = Int64(round(trunc( (p.y[k] -  yv[1])/dy ) + 1));
    if i<1
        i=1
    elseif i>ncx+1
        i=ncx+1
    end
    if j<1
        j=1
    elseif j>ncy
        j = ncy
    end
    # Compute distances
    dxmj = p.x[k] - xce[i]
    dymi = p.y[k] -  yv[j]
    # Compute vy velocity for the left and right of the cell
    vym12 = Vy[i,j  ]*(1-dymi/dy) + Vy[i  ,j+1]*dymi/dy
    vym34 = Vy[i+1,j]*(1-dymi/dy) + Vy[i+1,j+1]*dymi/dy
    if new==1
        if dymi/dy>=0.5
            if j<ncy
                vym12 += 0.5*((dymi/dy-0.5)^2) * ( Vy[i,j  ] - 2.0*Vy[i,j+1  ] + Vy[i,j+2  ]);
                vym34 += 0.5*((dymi/dy-0.5)^2) * ( Vy[i+1,j] - 2.0*Vy[i+1,j+1] + Vy[i+1,j+2]);
            end
        else
            if j>1
                vym12 += 0.5*((dymi/dy-0.5)^2) * ( Vy[i,j-1  ] - 2.0*Vy[i,j  ] + Vy[i,j+1  ]);
                vym34 += 0.5*((dymi/dy-0.5)^2) * ( Vy[i+1,j-1] - 2.0*Vy[i+1,j] + Vy[i+1,j+1]);
            end
        end
    end
    # Compute vy
    vym = (1-dxmj/dx)*vym12 + (dxmj/dx)*vym34
    return vym
end
export VyFromVxNodes

@views function LocateMarkers(p::Markers,dx::Float64,dy::Float64,xc,yc,xmin::Float64,xmax::Float64,ymin::Float64,ymax::Float64)
    # Find marker cell indices
    @threads for k=1:p.nmark
        if (p.x[k]<xmin || p.x[k]>xmax || p.y[k]<ymin || p.y[k]>ymax)
            p.phase[k] = -1
        end
        if p.phase[k]>=0
            dstx         = p.x[k] - xc[1]
            i            = ceil(Int, dstx / dx + 0.5)
            dsty         = p.y[k] - yc[1]
            j            = ceil(Int, dsty / dy + 0.5)
            p.cellx[k]   = i
            p.celly[k]   = j
        end
    end
end
export LocateMarkers

@views function Markers2Cells3!(p::Markers,phase::Matrix{Float64},xc::LinRange{Float64},yc::LinRange{Float64},dx::Float64,dy::Float64,ncx::Int64,ncy::Int64,prop::Vector{Float64},avg::Int64,mode::Int64)
    # if length(prop)  > p.nmark  println(length(prop)); println(p.nmark ); error("wtf") end
    phase0      = zeros(Float64, (ncx, ncy))
    phase0     .= phase
    weight      = zeros(Float64, (ncx, ncy))
    phase_th    = [similar(phase) for _ = 1:nthreads()] # per thread
    weight_th   = [similar(weight) for _ = 1:nthreads()] # per thread
    @threads for tid=1:nthreads()
        fill!(phase_th[tid] , 0)
        fill!(weight_th[tid], 0)
    end
    chunks = Iterators.partition(1:p.nmark, p.nmark ÷ nthreads())
    @sync for chunk in chunks
        Threads.@spawn begin
            tid = threadid()
            # fill!(phase_th[tid], 0)  # DON'T
            # fill!(weight_th[tid], 0)
            for k in chunk
                if p.phase[k]>=0
                # Get the indices:
                i = p.cellx[k]
                j = p.celly[k]
                # Relative distances
                dxm = 2.0 * abs(xc[i] - p.x[k])
                dym = 2.0 * abs(yc[j] - p.y[k])
                # Increment cell counts
                area = (1.0 - dxm / dx) * (1.0 - dym / dy)
                if mode==0 val  =  prop[Int64(p.phase[k])] end
                if mode==1 val  =  prop[k]                 end
                if avg==0 phase_th[tid][i,  j] += val       * area end
                if avg==1 phase_th[tid][i,  j] += (1.0/val) * area end
                if avg==2 phase_th[tid][i,  j] += log(val) * area end
                weight_th[tid][i, j] += area
                end
            end
        end
    end
    phase  .= reduce(+, phase_th)
    weight .= reduce(+, weight_th)
    phase ./= weight
    if avg==1
        phase .= 1.0 ./ phase
    end
    if avg==2
        phase .= exp.(phase)
    end
    phase[weight.<1e-10] .= phase0[weight.<1e-10]
    return
end
export Markers2Cells3!

@views function CountMarkers!(p::Markers,mpc,mpc_th,nmark,xmin::Float64,xmax::Float64,ymin::Float64,ymax::Float64,xc,yc,dx,dy::Float64,ncx::Int64,ncy::Int64)
    # Disable markers outside of the domain
    @threads for k=1:nmark
        if (p.x[k]<xmin || p.x[k]>xmax || p.y[k]<ymin || p.y[k]>ymax)
            @inbounds p.phase[k] = -1
        end
    end

    # How many are outside? save indices for reuse
    nmark_out_th = zeros(Int64, nthreads())
    @threads for k=1:nmark
        if p.phase[k] == -1
            nmark_out_th[threadid()] += 1
        end
    end
    nmark_out = 0
    for ith=1:nthreads()
        nmark_out += nmark_out_th[ith]
    end
    @printf("%d markers out\n", nmark_out)

    # Count number of marker per cell
    @threads for j=1:ncy
        for i=1:ncx
            for ith=1:nthreads()
                mpc_th[ith,i,j] = 0.0
            end
        end
    end

    @threads for k=1:nmark # @avx ne marche pas ici
        if (p.phase[k]>=0)
            # Get the column:
            dstx = p.x[k] - xc[1];
            i    = Int64(round(ceil( (dstx/dx) + 0.5)));
            # Get the line:
            dsty = p.y[k] - yc[1];
            j    = Int64(round(ceil( (dsty/dy) + 0.5)));
            # Increment cell count
            mpc_th[threadid(),i,j] += 1.0
        end
    end

    @threads for j=1:ncy
        for i=1:ncx
            for ith=1:nthreads()
                if ith == 1
                    mpc[i,j] = 0.0
                end
                mpc[i,j] += mpc_th[ith,i,j]
            end
        end
    end
    return
end
export CountMarkers!

@views function CountMarkers3!( p::Markers, f::Fields2D, mat::MaterialParameters, dom::ModelDomain )

    ncx, ncy, dx, dy = dom.ncx, dom.ncy, dom.dx, dom.dy
    nphase = mat.nphase
    nmark  = p.nmark

    # Disable markers outside of the domain
    @threads for k=1:nmark
        if (p.x[k]<dom.xmin || p.x[k]>dom.xmax || p.y[k]<dom.ymin || p.y[k]>dom.ymax)
            @inbounds p.phase[k] = -1
        end
    end

    # How many are outside? save indices for reuse
    nmark_out_th = zeros(Int64, nthreads())
    @threads for k=1:nmark
        if p.phase[k] == -1
            nmark_out_th[threadid()] += 1
        end
    end
    nmark_out = 0
    for ith=1:nthreads()
        nmark_out += nmark_out_th[ith]
    end
    @printf("%d markers out\n", nmark_out)

    # Count number of marker per cell
    @threads for j=1:ncy
        for i=1:ncx
            for ith=1:nthreads()
                f.mpc_th[ith,i,j] = 0.0
                for m=1:nphase
                    f.phase_perc_th[ith,m,i,j] = 0.0
                end
            end
        end
    end

    @threads for k=1:nmark # @avx ne marche pas ici
        if (p.phase[k]>=0)
            m = Int64(p.phase[k])
            # Get the column:
            dstx = p.x[k] - dom.xc[1];
            i    = Int64(round(ceil( (dstx/dx) + 0.5)));
            # Get the line:
            dsty = p.y[k] - dom.yc[1];
            j    = Int64(round(ceil( (dsty/dy) + 0.5)));
            # Increment cell count
            f.mpc_th[threadid(),i,j] += 1.0
            f.phase_perc_th[threadid(),m,i,j] += 1.0
        end
    end

    @threads for j=1:ncy
        for i=1:ncx
            for ith=1:nthreads()
                if ith == 1
                    f.mpc[i,j] = 0.0
                    for m=1:nphase
                        f.phase_perc[m,i,j] = 0.0
                    end
                end
                f.mpc[i,j] += f.mpc_th[ith,i,j]
                for m=1:nphase
                    f.phase_perc[m,i,j] += f.phase_perc_th[ith,m,i,j]
                end
            end
        end
    end

    # Normalize phase percentages
    for m=1:nphase
        f.phase_perc[m,:,:] ./= f.mpc
    end

    # Show infos
    println("markers per cell:")
    println("min: ", minimum(f.mpc), " --- max: ", maximum(f.mpc))
    for m=1:nphase
        println("phase ", m, " :")
        println("min: ", minimum(f.phase_perc[m,:,:]), " --- max: ", maximum(f.phase_perc[m,:,:]))
    end
end
export CountMarkers3!

@views function ReseedMarkers2!( p::Markers, mpc::Matrix{Float64}, dom::ModelDomain)

    ncx, ncy, dx, dy = dom.ncx, dom.ncy, dom.dx, dom.dy
    mlist       = Array{Array{Int64}}(undef,ncx,ncy)
    mlist_th    = Array{Array{Int64}}(undef,nthreads(),ncx,ncy)

    @tturbo for i=1:ncx
        for j=1:ncy
            mlist[i,j] = zeros(Int64(mpc[i,j]))
        end
    end

    count = zeros(Int64,ncx,ncy)

    for k=1:p.nmark # @avx ne marche pas ici
        if p.phase[k]>=0
            # Get the column:
            dstx = p.x[k] - dom.xc[1];
            i    = Int64(round(ceil( (dstx/dx) + 0.5)));
            # Get the line:
            dsty = p.y[k] - dom.yc[1];
            j    = Int64(round(ceil( (dsty/dy) + 0.5)));
            # Increment cell count
            count[i,j] += 1
            mlist[i,j][count[i,j]] = k
        end
    end

    # display(mlist)
    dxy = 10*(dx*dx + dy*dy)
    for i=1:ncx
        for j=1:ncy

            if mpc[i,j] < 4

                new_x = dom.xc[i]-dx/4.0
                new_y = dom.yc[j]-dy/4.0
                ineigh = FindClosestNeighbour( p, mpc, mlist, new_x, new_y, dxy, i, j, ncx, ncy )
                if ineigh>0
                    p.nmark+=1
                    p.x[p.nmark] = new_x
                    p.y[p.nmark] = new_y
                    CopyMarkers!(p,ineigh)
                end

                new_x = dom.xc[i]+dx/4.0
                new_y = dom.yc[j]-dy/4.0
                ineigh = FindClosestNeighbour( p, mpc, mlist, new_x, new_y, dxy, i, j, ncx, ncy )
                if ineigh>0
                    p.nmark+=1
                    p.x[p.nmark] = new_x
                    p.y[p.nmark] = new_y
                    CopyMarkers!(p,ineigh)
                end

                new_x = dom.xc[i]+dx/4.0
                new_y = dom.yc[j]+dy/4.0
                ineigh = FindClosestNeighbour( p, mpc, mlist, new_x, new_y, dxy, i, j, ncx, ncy )
                if ineigh>0
                    p.nmark+=1
                    p.x[p.nmark] = new_x
                    p.y[p.nmark] = new_y
                    CopyMarkers!(p,ineigh)
                end

                new_x = dom.xc[i]-dx/4.0
                new_y = dom.yc[j]+dy/4.0
                ineigh = FindClosestNeighbour( p, mpc, mlist, new_x, new_y, dxy, i, j, ncx, ncy )
                if ineigh>0
                    p.nmark+=1
                    p.x[p.nmark] = new_x
                    p.y[p.nmark] = new_y
                    CopyMarkers!(p,ineigh)
                end

            end
        end
    end
end
export ReseedMarkers2!

@views function LocateMarkers2!(p::Markers,dom::ModelDomain)
    # Find marker cell indices
    @threads for k=1:p.nmark
        if (p.x[k]<dom.xmin || p.x[k]>dom.xmax || p.y[k]<dom.ymin || p.y[k]>dom.ymax)
            p.phase[k] = -1
        end
        if p.phase[k]>=0
            dstx         = p.x[k] - dom.xc[1]
            i            = ceil(Int, dstx / dom.dx + 0.5)
            dsty         = p.y[k] - dom.yc[1]
            j            = ceil(Int, dsty / dom.dy + 0.5)
            p.cellx[k]   = i
            p.celly[k]   = j
        end
    end
end
export LocateMarkers2!

################################################################
################################################################

function RungeKutta2!(p::Markers, f::Fields2D, params::ModelParameters, BC::BoundaryConditions, dom::ModelDomain)

    nmark = p.nmark
    rkv, rkw, dt   = params.rkv, params.rkw, params.dt
    ncx, ncy, dx, dy = dom.ncx, dom.ncy, dom.dx, dom.dy

    # Marker advection with 4th order Roger-Kutta
    @threads for k=1:nmark
        in = p.phase[k]>=0
            x0 = p.x[k];
            y0 = p.y[k];
            vx = 0.0
            vy = 0.0
            # Roger-Gunther loop
            for rk=1:4
                # println("step", rk)
                # Interp velocity from grid
                vxm = VxFromVxNodes(f.Vx, k, p, dom.xv, dom.yce, dx, dy, ncx, ncy, 1)
                vym = VyFromVxNodes(f.Vy, k, p, dom.xce, dom.yv, dx, dy, ncx, ncy, 1)
                # Temporary RK advection steps
                p.x[k] = x0 + rkv[rk]*dt*vxm
                p.y[k] = y0 + rkv[rk]*dt*vym
                if BC.periodix==1
                    if p.x[k]<dom.xv[1]
                        # println("periW 1 ", p.x[k], vxm)
                        p.x[k]+=(dom.xv[end]-dom.xv[1])
                        # println("periW 2 ", p.x[k], vxm)
                    end
                    if p.x[k]>dom.xv[end]
                        # println("periE")
                        p.x[k]-=(dom.xv[end]-dom.xv[1])
                    end
                end
                # Average final velocity
                vx    += rkw[rk]*vxm
                vy    += rkw[rk]*vym
            end
            # Advect points
            p.x[k] = x0 + (in==1) * rkv[4]*dt*vx
            p.y[k] = y0 + (in==1) * rkv[4]*dt*vy
            if BC.periodix==1
                if p.x[k]<dom.xv[1]
                    # println("periW 1 ", p.x[k], vxm)
                    p.x[k]+=(dom.xv[end]-dom.xv[1])
                    # println("periW 2 ", p.x[k], vxm)
                end
                if p.x[k]>dom.xv[end]
                    # println("periE")
                    p.x[k]-=(dom.xv[end]-dom.xv[1])
                end
            end
    end
    end
    export RungeKutta2!