text:
#
# include(joinpath(pwd(),"src_PowSysMod", "PowSysMod_body.jl"))
#
# raw = "powersystem.raw"
# gen = "generator.csv"
# con = "contingency.csv"
#
# instance_path = joinpath(pwd(), "data_GOC", "Phase_0_IEEE14_1Scenario","scenario_1")
# instance_path = joinpath(pwd(), "data_GOC", "Phase_0_RTS96","scenario_1")
# instance_path = joinpath(pwd(), "data_GOC", "Phase_0_Feas179","scenario_1")
#
# read(matopen(joinpath(instance_path, "pscopf_data.mat")))
#
# rawfile = joinpath(instance_path,raw)
# genfile = joinpath(instance_path, gen)
# confile = joinpath(instance_path, con)
#
# baseMVA, bus_data, load_data, shunt_data, generator_data, branch_data, transformer_data = read_rawfile(rawfile)
# gen_data_csv = read_genfile(genfile)
# con_data_csv = read_confile(confile)
#
# read_GOCfiles(rawfile, genfile,confile)

function read_genfile(genfile)
    gen_data = readdlm(genfile,',',skipstart=1)
    return gen_data
end

function read_confile(confile)
    con_data = readdlm(confile,',',skipstart=1)
    return con_data
end

function read_rawfile(rawfile)
    raw_data = readdlm(rawfile, ',')
    baseMVA = raw_data[1,2]
    i_begin = 2
    while raw_data[i_begin,4] == ""
        i_begin +=1
    end
    ##BUS DATA
    i_end = get_data_iend(i_begin,raw_data, "0 / END OF BUS DATA")
    bus_data = raw_data[i_begin:i_end,:]
    ##LOAD DATA
    i_begin = i_end + 2
    i_end = get_data_iend(i_begin,raw_data, "0 / END OF LOAD DATA")
    load_data = raw_data[i_begin:i_end,:]
    ## FIXED SHUNT DATA
    i_begin = i_end + 2
    i_end = get_data_iend(i_begin,raw_data, "0 / END OF FIXED SHUNT DATA")
    shunt_data = raw_data[i_begin:i_end,:]
    ## GENERATOR DATA
    i_begin = i_end + 2
    i_end = get_data_iend(i_begin,raw_data, "0 / END OF GENERATOR DATA")
    generator_data = raw_data[i_begin:i_end,:]
    ## BRANCH DATA
    i_begin = i_end + 2
    i_end = get_data_iend(i_begin,raw_data, "0 / END OF BRANCH DATA")
    branch_data = raw_data[i_begin:i_end,:]
    ## TRANSFORMER DATA
    i_begin = i_end + 2
    i_end = get_data_iend(i_begin,raw_data, "0 / END OF TRANSFORMER DATA")
    transformer_data = raw_data[i_begin:i_end,:]
    return baseMVA, bus_data, load_data, shunt_data, generator_data, branch_data, transformer_data
end

function get_data_iend(i_begin, data, end_sentence)
    i = i_begin
    imax = size(data,1)
    while data[i,1] != end_sentence && i<imax
        i+=1
    end
    return i-1
end

###read bus data
function read_data_bus_fromraw(bus_data, load_data, shunt_data, baseMVA)
    nb_bus = size(bus_data,1)
    # index = SortedDict(bus_data[i,1] => i for i in 1:nb_bus)
    index = SortedDict(bus_data[i,1] => bus_data[i,1] for i in 1:nb_bus)
    bus = SortedDict{String, SortedDict{String,Any}}()
    for i in 1:nb_bus
        bus[bus_name(index[bus_data[i,1]])] = SortedDict{String,Any}()
    end
    for i in 1:nb_bus
        id_bus = bus_data[i,1]
        busname = bus_name(index[id_bus])
        baseKV = bus_data[i,3]
        baseMVA = baseMVA
        voltage_magnitude_min = bus_data[i,11]
        voltage_magnitude_max  = bus_data[i,10]
        bus[busname][volt_name()] = GOCVolt(id_bus, baseKV, baseMVA, voltage_magnitude_min, voltage_magnitude_max)
    end
    ### LOAD DATA
    nb_load = size(load_data,1)
    for i in 1:nb_load
        id_bus = load_data[i,1]
        busname = bus_name(index[id_bus])
        load_re = load_data[i,6]
        load_im = load_data[i,7]
        if load_re!=0 || load_im!=0
            load = load_re  + im * load_im
            id_load = load_data[i,2]
            id_load = remove_simple_quotes_and_spaces_if_present(id_load)
            loadname = load_name(id_load)
            bus[busname][loadname] = GOCLoad(id_bus,loadname,load)
        end
    end
    ### SHUNT DATA
    nb_shunt = size(shunt_data,1)
    for i in 1:nb_shunt
        id_bus = shunt_data[i,1]
        busname = bus_name(index[id_bus])
        shunt_re = shunt_data[i,4]
        shunt_im = shunt_data[i,5]
        if shunt_re!=0 || shunt_im!=0
            shunt = shunt_re + im * shunt_im
            id_shunt = shunt_data[i,2]
            id_shunt = remove_simple_quotes_and_spaces_if_present(id_shunt)
            shuntname = shunt_name(id_shunt)
            bus[busname][shuntname] = GOCShunt(id_bus,shuntname,shunt)
        end
    end
    return bus,index
end

function remove_simple_quotes_and_spaces_if_present(id)
    if typeof(id)==SubString{String}
        if id[1]==''' && id[end]=='''
            id_without_spaces = split(id[2:end-1])[1]
        else
            id_without_spaces = split(id)[1]
    end
        return id_without_spaces
    else
        return id
end

end

function convert_gen_data_csv_into_dict(gen_data_csv, index)
    nb_lines = size(gen_data_csv, 1)
    gen_data_csv_dict = SortedDict{String, SortedDict{String,SortedDict{Int64,Float64}}}()
    for i in 1:nb_lines
        id_bus = Int64(gen_data_csv[i,1])
        busname = bus_name(index[id_bus])
        gen_id = gen_data_csv[i,2]
        gen_id = remove_simple_quotes_and_spaces_if_present(gen_id)
        if typeof(gen_id)==Float64
            gen_id = Int64(gen_id)
        end
        generatorname = generator_name(gen_id)
        term_id = Int64(gen_data_csv[i,3])
        coeff = gen_data_csv[i,4]
        if !haskey(gen_data_csv_dict, busname)
            gen_data_csv_dict[busname] = SortedDict{String,SortedDict{Int64,Float64}}()
        end
        if !haskey(gen_data_csv_dict[busname], generatorname)
            gen_data_csv_dict[busname][generatorname] = SortedDict{Int64,Float64}()
        end
        gen_data_csv_dict[busname][generatorname][term_id] = coeff
    end
    return gen_data_csv_dict
end

function add_generator_data_fromraw!(generator_data, gen_data_csv_dict, bus, index)
    nb_gen = size(generator_data,1)
    for i in 1:nb_gen
        id_bus = generator_data[i,1]
        busname = bus_name(index[id_bus])
        gen_id = generator_data[i,2]
        gen_id2 = remove_simple_quotes_and_spaces_if_present(gen_id)
        if typeof(gen_id2)==Float64
            gen_id2 = Int64(gen_id)
        end
        generatorname = generator_name(gen_id2)
        Pmin = generator_data[i,18]
        Qmin = generator_data[i,6]
        Pmax = generator_data[i,17]
        Qmax = generator_data[i,5]
        participation_factor = gen_data_csv_dict[busname][generatorname][9]
        power_min = Pmin + im*Qmin
        power_max = Pmax + im*Qmax
        dict_obj_coeffs = SortedDict{Int64,Float64}()
        for (degree, value) in gen_data_csv_dict[busname][generatorname]
            if (degree ∈ [0,1,2])
                dict_obj_coeffs[degree] = value
            end
        end
        if length(dict_obj_coeffs)<3
            #warn("A 2-order polynom cost with less than 3 values for generator $generatorname at bus $busname, default value 0")
            if !haskey(dict_obj_coeffs,0) dict_obj_coeffs[0] = 0 end
            if !haskey(dict_obj_coeffs,1) dict_obj_coeffs[1] = 0 end
            if !haskey(dict_obj_coeffs,2) dict_obj_coeffs[2] = 0 end
       end
        bus[busname][generatorname] = GOCGenerator(id_bus,gen_id,power_min,power_max,participation_factor,dict_obj_coeffs)
    end
    return bus
end

function modify_transformer_data(transformer_data)
    nb_lines, nb_col = size(transformer_data)
    nb_transfo = nb_lines/4
    if !isinteger(nb_transfo)
        exit("Transformer data is not on the format : 4 lines for 1 transformer")
    end
    nb_transfo = Int64(nb_transfo)
    transfo_data = Array{Any}(nb_transfo,nb_col*4)

    for i in 1:nb_transfo
        line = (i-1)*4+1
        transfo_data[i,1:nb_col] = transformer_data[line,:]
        transfo_data[i,(nb_col+1):(2*nb_col)] = transformer_data[(line+1),:]
        transfo_data[i,(2*nb_col+1):(3*nb_col)] = transformer_data[(line+2),:]
        transfo_data[i,(3*nb_col+1):(4*nb_col)] = transformer_data[(line+3),:]
    end
    return transfo_data
end

function read_data_branch_fromraw(branch_data, transfo_data, index)
    nb_branch = size(branch_data,1)
    link = SortedDict{Link, SortedDict{String,Any}}()
    for i in 1:nb_branch
        orig = Int64(branch_data[i,1])
        dest = Int64(branch_data[i,2])
        linkname = Link(bus_name(index[orig]),bus_name(index[dest]))
        br_id = branch_data[i,3]
        branch_id = remove_simple_quotes_and_spaces_if_present(branch_data[i,3])
        if typeof(branch_id)==Float64
            branch_id = Int64(branch_id)
        end
        resistance = branch_data[i,4]
        reactance = branch_data[i,5]
        susceptance = branch_data[i,6]
        power_magnitude_max = branch_data[i,7]

        if !haskey(link, linkname)
            link[linkname] = SortedDict{String,Any}()
        end
        if resistance == reactance == 0
            println(linkname, " nullimpedance line without transformer")
            name_line = nullimpedance_notransformer_name(branch_id)
            link[linkname][name_line] = GOCNullImpedance_notransformer(orig,dest,br_id,susceptance, power_magnitude_max)
        else
            name_line = linepi_notransformer_name(branch_id)
            link[linkname][name_line] = GOCLineπ_notransformer(orig,dest,br_id,resistance,reactance,susceptance, power_magnitude_max)
        end
    end

    nb_transfo = size(transfo_data,1)
    nb_col = size(branch_data,2)
    for i in 1:nb_transfo
        orig = Int64(transfo_data[i,1])
        dest = Int64(transfo_data[i,2])
        linkname = Link(bus_name(index[orig]),bus_name(index[dest]))
        br_id = transfo_data[i,4]
        branch_id = remove_simple_quotes_and_spaces_if_present(transfo_data[i,4])
        if typeof(branch_id)==Float64
            branch_id = Int(branch_id)
        end
        resistance = transfo_data[i,(nb_col+1)]
        reactance = transfo_data[i,(nb_col+2)]
        susceptance = 0.0
        power_magnitude_max = transfo_data[i,(2*nb_col+4)]
        transfo_ratio = transfo_data[i,(2*nb_col+1)]
        transfo_phase = transfo_data[i,(2*nb_col+3)]
        if !haskey(link, linkname)
            link[linkname] = SortedDict{String,Any}()
        end
        if resistance == reactance == 0
            println(linkname, " nullimpedance line with transformer")
            name_line = nullimpedance_withtransformer_name(branch_id)
            link[linkname][name_line] = GOCNullImpedance_withtransformer(orig,dest,br_id,susceptance, transfo_ratio,transfo_phase, power_magnitude_max)
        else
            name_line = linepi_withtransformer_name(branch_id)
            link[linkname][name_line] = GOCLineπ_withtransformer(orig,dest,br_id,resistance,reactance,susceptance, transfo_ratio,transfo_phase, power_magnitude_max)
        end
    end

    return link
end

function read_GOCfiles(rawfile, genfile,confile)
    baseMVA, bus_data, load_data, shunt_data, generator_data, branch_data, transformer_data = read_rawfile(rawfile)
    gen_data_csv = read_genfile(genfile)
    con_data_csv = read_confile(confile)
    bus,index = read_data_bus_fromraw(bus_data, load_data, shunt_data, baseMVA)
    gen_data_csv_dict = convert_gen_data_csv_into_dict(gen_data_csv, index)
    bus = add_generator_data_fromraw!(generator_data, gen_data_csv_dict, bus, index)
    transfo_data = modify_transformer_data(transformer_data)
    link = read_data_branch_fromraw(branch_data, transfo_data, index)

    ds = DataSource(bus,link)
    node_linksin, node_linksout = SortedDict{String, SortedSet{Link}}(), SortedDict{String, SortedSet{Link}}()
    node_vars = SortedDict{String, SortedDict{String, Variable}}()
    link_vars = SortedDict{Link, SortedDict{String, Variable}}()
    gs = GridStructure("BaseCase", node_linksin, node_linksout)
    node_formulations = SortedDict{String, SortedDict{String, Symbol}}()
    link_formulations = SortedDict{Link, SortedDict{String, Symbol}}()
    mp = MathematicalProgramming(node_formulations, link_formulations, node_vars,link_vars)
    ##read scenarios
    OPFproblems = scenarios_data(ds, gs, mp, con_data_csv,index)
    return OPFproblems
end


text:
"""
    A short code to show how this package was made
    according to a template. This is a utility script that
    was run just once to create the package and will not
    be needed for this project again.

    Steps:
    1. Use PkgTemplates to make git repo.
    2. new empty repository on GitHub.
    3. Then push an existing repository from the command line:
    `git remote add origin git@github.com:ggebbie/TMI.jl.git`
    `git branch -M main`
    `git push -u origin main`

    First two commands probably unnecessary.
    I didn't figure out how to do this with magit in emacs.
"""
using Revise, PkgTemplates

t = Template(;
    user="ggebbie",
    dir="~/projects",
    authors="G Jake Gebbie",
    julia=v"1.6",
    plugins=[
        License(; name="MIT"),
        Git(; manifest=true, ssh=true),
        GitHubActions(; x86=true),
        Codecov(),
        Documenter{GitHubActions}(),
        Develop(),
    ],
             )

t("TMI.jl")

ext: jl


text:
"""
    _heappush!(heap::BinaryMaxHeap, edge::HeapKNNGraphEdge, max_length)

Try to push `edge` to `heap`. This will fail if `edge` is
already in `heap`.
If `length(heap) > max_length` after pushing, `pop` the largest edge off.
"""
function _heappush!(heap::BinaryMaxHeap,
                    edge::HeapKNNGraphEdge,
                    max_length::Integer)

    # edge not already in heap
    if !(edge in heap.valtree)
        push!(heap, edge)
    end

    if length(heap) > max_length
        pop!(heap)
    end

    return
end

"""
Sample `n_neighbors` elements from a set of ints `1:npoints`.
The ints in `exclude` won't be sampled.
"""
function sample_neighbors(npoints::Int,
                          n_neighbors::Int,
                          sample_rate::R = 1.;
                          exclude::Vector{Int} = Vector{Int}()) where {R <: AbstractFloat}
    # num neighbors to sample = max(ρK, N)
    last = min(npoints-length(exclude), trunc(Int, sample_rate*n_neighbors))
    ids = Set{Int}()
    while length(ids) < last
        k = (abs(rand(Int)) % npoints) + 1
        if k ∉ exclude
            union!(ids, k)
        end
    end
    return collect(ids)
end


text:
function slider!(ec::EChart)

    #Need to be more robust in checking the structure of dataZoom, specify axis, etc.
    ec.dataZoom == nothing ? ec.dataZoom = [DataZoom(show = true)] : push!(ec.dataZoom, DataZoom(show = true))

    return ec

end


text:
function rollout!(policy::PolicyData, problem::ProblemData;
    step_size=1.0)

    # model
    dynamics = problem.model.dynamics

    # trajectories
    x = problem.states
    u = problem.actions
    w = problem.parameters
    x̄ = problem.nominal_states
    ū = problem.nominal_actions

    # policy
    K = policy.K
    k = policy.k

    # initial state
    x[1] .= x̄[1]

    # rollout
    for (t, d) in enumerate(dynamics)
        # u[t] .= ū[t] + K[t] * (x[t] - x̄[t]) + step_size * k[t]
        u[t] .= k[t]
        u[t] .*= step_size
        u[t] .+= ū[t]
        mul!(u[t], K[t], x[t], 1.0, 1.0)
        mul!(u[t], K[t], x̄[t], -1.0, 1.0)
        x[t+1] .= dynamics!(d, x[t], u[t], w[t])
    end
end

function rollout(dynamics::Vector{Dynamics{T}}, initial_state, actions,
    parameters=[zeros(d.num_parameter) for d in dynamics]) where T

    x_history = [initial_state]
    for (t, d) in enumerate(dynamics)
        push!(x_history, copy(dynamics!(d, x_history[end], actions[t], parameters[t])))
    end

    return x_history
end


text:
using GeometryBasics, RPRMakie
using Colors, FileIO
using Colors: N0f8

image = begin
    radiance = 500
    lights = [EnvironmentLight(1.0, load(RPR.assetpath("studio026.exr"))),
              PointLight(Vec3f(10), RGBf(radiance, radiance, radiance * 1.1))]
    fig = Figure(; resolution=(1500, 700))
    ax = LScene(fig[1, 1]; show_axis=false, scenekw=(lights=lights,))
    screen = RPRScreen(ax.scene; plugin=RPR.Northstar, iterations=400)

    matsys = screen.matsys
    emissive = RPR.EmissiveMaterial(matsys)
    diffuse = RPR.DiffuseMaterial(matsys)
    glass = RPR.Glass(matsys)
    plastic = RPR.Plastic(matsys)
    chrome = RPR.Chrome(matsys)
    dielectric = RPR.DielectricBrdfX(matsys)
    gold = RPR.SurfaceGoldX(matsys)

    materials = [glass chrome;
                 gold dielectric;
                 emissive plastic]

    mesh!(ax, load(Makie.assetpath("matball_floor.obj")); color=:white)
    palette = reshape(Makie.default_palettes.color[][1:6], size(materials))

    for i in CartesianIndices(materials)
        x, y = Tuple(i)
        mat = materials[i]
        mplot = if mat === emissive
            matball!(ax, diffuse; inner=emissive, color=nothing)
        else
            matball!(ax, mat; color=nothing)
        end
        v = Vec3f(((x, y) .- (0.5 .* size(materials)) .- 0.5)..., 0)
        translate!(mplot, 0.9 .* (v .- Vec3f(0, 3, 0)))
    end
    cam = cameracontrols(ax.scene)
    cam.eyeposition[] = Vec3f(-0.3, -5.5, 0.9)
    cam.lookat[] = Vec3f(0.5, 0, -0.5)
    cam.upvector[] = Vec3f(0, 0, 1)
    cam.fov[] = 35
    emissive.color = Vec3f(4, 2, 2)
    colorbuffer(screen)
end

save("materials.png", image)


text:
import BioSimulator: MassAction
import BioSimulator: ReactionStruct, execute_jump!, rate, rate_derivative
import BioSimulator: ReactionSystem

# define notion of equality for ReactionStruct
import Base: ==

function ==(r1::ReactionStruct{T}, r2::ReactionStruct{S}) where {T,S}
  isequal = true
  isequal = isequal && T == S
  isequal = isequal && length(r1.reactants) == length(r2.reactants)
  isequal = isequal && all(x -> x in r2.reactants, r1.reactants)
  isequal = isequal && length(r1.net_change) == length(r2.net_change)
  isequal = isequal && all(x -> x in r2.net_change, r2.net_change)

  return isequal
end

@testset "ReactionSystem" begin
  x = rand(0:1000, 5)
  params = 10 * rand(4)

  # Species | A, B, C, D E
  #   index | 1  2  3  4 5

  # Order 0
  order0 = [
    (Tuple{Int,Int}[], [(1, 1)], 1), # 0 --> A
    (Tuple{Int,Int}[], [(2, 4)], 1)  # 0 --> 4 B
  ]

  # Order 1
  order1 = [
    ([(1, 1)], [(1, -1)],                 2), # A --> 0
    ([(1, 1)], [(1, -1), (2, 1)],         2), # A --> B
    ([(1, 1)], [(1, -1), (3, 2), (5, 2)], 2)  # A --> 2 C + 2 E
  ]

  # Order 2A
  order2a = [
    ([(1, 1), (2, 1)], [(1, -1), (2, -1), (3, 1)], 3) # A + B --> C
  ]

  # Order 2B
  order2b = [
    ([(2, 2)], [(2, -2), (1, 1)], 4) # 2 B --> A
  ]

  # Order 3
  # A + B + C --> D
  order3a = [
    ([(1, 1), (2, 1), (3, 1)], [(1, -1), (2, -1), (3, -1), (4, 1)], 1)
  ]

  # 3 * A + B + 2 * C --> D
  order3b = [
    ([(1, 3), (2, 1), (3, 2)], [(1, -3), (2, -1), (3, -2), (4, 1)], 1)
  ]

  @testset "ReactionStruct" begin
    rxn0  = [ReactionStruct(MassAction(), r, v, i) for (r, v, i) in order0]
    rxn1  = [ReactionStruct(MassAction(), r, v, i) for (r, v, i) in order1]
    rxn2a = [ReactionStruct(MassAction(), r, v, i) for (r, v, i) in order2a]
    rxn2b = [ReactionStruct(MassAction(), r, v, i) for (r, v, i) in order2b]
    rxn3a = [ReactionStruct(MassAction(), r, v, i) for (r, v, i) in order3a]
    rxn3b = [ReactionStruct(MassAction(), r, v, i) for (r, v, i) in order3b]

    @testset "fire_reaction!" begin
      result   = copy(x)
      expected = x + [1, 0, 0, 0, 0]
      execute_jump!(result, rxn0[1])
      @test result == expected

      result   = copy(x)
      expected = x + [0, 4, 0, 0, 0]
      execute_jump!(result, rxn0[2])
      @test result == expected

      result   = copy(x)
      expected = x + [-1, 0, 0, 0, 0]
      execute_jump!(result, rxn1[1])
      @test result == expected

      result   = copy(x)
      expected = x + [-1, 1, 0, 0, 0]
      execute_jump!(result, rxn1[2])
      @test result == expected

      result   = copy(x)
      expected = x + [-1, 0, 2, 0, 2]
      execute_jump!(result, rxn1[3])
      @test result == expected

      result   = copy(x)
      expected = x + [-1, -1, 1, 0, 0]
      execute_jump!(result, rxn2a[1])
      @test result == expected

      result   = copy(x)
      expected = x + [1, -2, 0, 0, 0]
      execute_jump!(result, rxn2b[1])
      @test result == expected

      result   = copy(x)
      expected = x + [-1, -1, -1, 1, 0]
      execute_jump!(result, rxn3a[1])
      @test result == expected

      result   = copy(x)
      expected = x + [-3, -1, -2, 1, 0]
      execute_jump!(result, rxn3b[1])
      @test result == expected
    end

    @testset "rate" begin
      @test rate(rxn0[1], x, params) ≈ params[1]
      @test rate(rxn0[2], x, params) ≈ params[1]

      @test rate(rxn1[1], x, params) ≈ x[1] * params[2]
      @test rate(rxn1[2], x, params) ≈ x[1] * params[2]
      @test rate(rxn1[3], x, params) ≈ x[1] * params[2]

      @test rate(rxn2a[1], x, params) ≈ x[1] * x[2] * params[3]
      @test rate(rxn2b[1], x, params) ≈ 1//2 * x[2] * (x[2] - 1) * params[4]

      @test rate(rxn3a[1], x, params) ≈ x[1] * x[2] * x[3] * params[1]
      @test rate(rxn3b[1], x, params) ≈ x[1] * (x[1] - 1) * (x[1] - 2) * x[2] * x[3] * (x[3] - 1) * params[1] / 12
    end
  end

  @testset "building from a Network" begin
    model = Network("silly test")

    # species definitions
    model <= Species("A", x[1])
    model <= Species("B", x[2])
    model <= Species("C", x[3])
    model <= Species("D", x[4])
    model <= Species("E", x[5])

    # order 0
    model <= Reaction("R1", params[1], "0 --> A")
    model <= Reaction("R2", params[1], "0 --> 4 * B")

    # order 1
    model <= Reaction("R3", params[2], "A --> 0")
    model <= Reaction("R4", params[2], "A --> B")
    model <= Reaction("R5", params[2], "A --> 2 * C + 2 * E")

    # order 2a
    model <= Reaction("R6", params[3], "A + B --> C")

    # order 2b
    model <= Reaction("R7", params[4], "2 * B --> A")

    # order 3
    model <= Reaction("R8", params[1], "A + B + C --> D")
    model <= Reaction("R9", params[1], "3 * A + B + 2 * C --> D")

    m = ReactionSystem(model)

    @testset "MassActionOrder0" begin
      expected = ReactionStruct(MassAction(), order0[1][1], order0[1][2], 1)

      @test m.reactions[1] == expected
      @test m.reactions[1].paramidx == expected.paramidx

      expected = ReactionStruct(MassAction(), order0[2][1], order0[2][2], 2)
      @test m.reactions[2] == expected
      @test m.reactions[2].paramidx == expected.paramidx
    end

    @testset "MassActionOrder1" begin
      expected = ReactionStruct(MassAction(), order1[1][1], order1[1][2], 3)

      @test m.reactions[3] == expected
      @test m.reactions[3].paramidx == expected.paramidx

      expected = ReactionStruct(MassAction(), order1[2][1], order1[2][2], 4)

      @test m.reactions[4] == expected
      @test m.reactions[4].paramidx == expected.paramidx

      expected = ReactionStruct(MassAction(), order1[3][1], order1[3][2], 5)

      @test m.reactions[5] == expected
      @test m.reactions[5].paramidx == expected.paramidx
    end

    @testset "MassActionOrder2A" begin
      expected = ReactionStruct(MassAction(), order2a[1][1], order2a[1][2], 6)

      @test m.reactions[6] == expected
      @test m.reactions[6].paramidx == expected.paramidx
    end

    @testset "MassActionOrder2B" begin
      expected = ReactionStruct(MassAction(), order2b[1][1], order2b[1][2], 7)

      @test m.reactions[7] == expected
      @test m.reactions[7].paramidx == expected.paramidx
    end

    @testset "MassActionOrderN" begin
      expected = ReactionStruct(MassAction(), order3a[1][1], order3a[1][2], 8)

      @test m.reactions[8] == expected
      @test m.reactions[8].paramidx == expected.paramidx

      expected = ReactionStruct(MassAction(), order3b[1][1], order3b[1][2], 9)
      @test m.reactions[9] == expected
      @test m.reactions[9].paramidx == expected.paramidx
    end

    @testset "fire_reaction!" begin
      result = copy(x); expected = x + [1, 0, 0, 0, 0]
      execute_jump!(result, m, 1)
      @test result == expected

      result = copy(x); expected = x + [0, 4, 0, 0, 0]
      execute_jump!(result, m, 2)
      @test result == expected

      result = copy(x); expected = x + [-1, 0, 0, 0, 0]
      execute_jump!(result, m, 3)
      @test result == expected

      result = copy(x); expected = x + [-1, 1, 0, 0, 0]
      execute_jump!(result, m, 4)
      @test result == expected

      result = copy(x); expected = x + [-1, 0, 2, 0, 2]
      execute_jump!(result, m, 5)
      @test result == expected

      result = copy(x); expected = x + [-1, -1, 1, 0, 0]
      execute_jump!(result, m, 6)
      @test result == expected

      result = copy(x); expected = x + [1, -2, 0, 0, 0]
      execute_jump!(result, m, 7)
      @test result == expected

      result = copy(x); expected = x + [-1, -1, -1, 1, 0]
      execute_jump!(result, m, 8)
      @test result == expected

      result = copy(x); expected = x + [-3, -1, -2, 1, 0]
      execute_jump!(result, m, 9)
      @test result == expected
    end

    @testset "rate" begin
      @test rate(m, x, 1) ≈ params[1]
      @test rate(m, x, 2) ≈ params[1]

      @test rate(m, x, 3) ≈ x[1] * params[2]
      @test rate(m, x, 4) ≈ x[1] * params[2]
      @test rate(m, x, 5) ≈ x[1] * params[2]

      @test rate(m, x, 6) ≈ x[1] * x[2] * params[3]
      @test rate(m, x, 7) ≈ 0.5 * x[2] * (x[2] - 1) * params[4]

      @test rate(m, x, 8) ≈ x[1] * x[2] * x[3] * params[1]
      @test rate(m, x, 9) ≈ x[1] * (x[1] - 1) * (x[1] - 2) * x[2] * x[3] * (x[3] - 1) * params[1] / 12
    end

    @testset "rate_derivative" begin
      # Order 0:
      @test rate_derivative(m, x, 1, 1) == 0

      # Order 1: A --> B
      @test rate_derivative(m, x, 1, 3) ≈ params[2]
      @test rate_derivative(m, x, 2, 3) == 0

      # Order 2A: A + B --> C
      @test rate_derivative(m, x, 1, 6) ≈ x[2] * params[3]
      @test rate_derivative(m, x, 2, 6) ≈ x[1] * params[3]
      @test rate_derivative(m, x, 3, 6) == 0

      # Order 2B: 2*B --> A
      @test rate_derivative(m, x, 2, 7) ≈ (x[2] - 0.5) * params[4]
      @test rate_derivative(m, x, 3, 7) == 0

      # Order 3: A + B + C --> D
      @test rate_derivative(m, x, 1, 8) ≈ x[2] * x[3] * params[1]
      @test rate_derivative(m, x, 2, 8) ≈ x[1] * x[3] * params[1]
      @test rate_derivative(m, x, 3, 8) ≈ x[1] * x[2] * params[1]
      @test rate_derivative(m, x, 4, 8) == 0

      # Order 3: 3 * A + B + 2 * C --> D
      @test rate_derivative(m, x, 1, 9) ≈ (0.5*x[1]^2 - x[1] + 1/3) * x[2] * x[3] * (x[3] - 1) * params[1] / 2
      @test rate_derivative(m, x, 2, 9) ≈ x[1] * (x[1] - 1) * (x[1] - 2) * x[3] * (x[3] - 1) * params[1] / 12
      @test rate_derivative(m, x, 3, 9) ≈ x[1] * (x[1] - 1) * (x[1] - 2) * x[2] * (x[3] - 0.5) * params[1] / 6
      @test rate_derivative(m, x, 4, 9) == 0
    end
  end
end


text:
using SmartShockFinder

# Determine the type of Network
n_dims = 2
network_type = NNPP(n_dims)

# add mesh and polynomdegrees to then generate Traindata
meshes = SmartShockFinder.get_mesh_2d(6)
polydeg = collect(2:9)
train_setup = Traindata_Settings(polydeg, meshes)

# Parameters to train the Net
η = 0.005                # learning rate
β = 0.005               # regularization paramater
number_epochs = 100
Sb = 500                 # batch size
L = 10                   # early stopping parameter
hidden_units = [32, 16, 16, 8, 8]

model_setup = Model_Settings(η, β, number_epochs, Sb, L, hidden_units)

# generate Traindata
generate_traindataset2d(network_type, train_setup)

# generate validationdata
generate_validdataset2d(network_type, train_setup)

# alternative to generate datasets

# traindata = joinpath("datasets","2d", "traindata2dNNPP.h5")
# download("https://gist.github.com/JuliaOd/f4d10153e539c026daf7076f087dc937/raw/traindata2dNNPP.h5", traindata)
# validdata = joinpath("datasets","2d","validdata2dNNPP.h5")
# download("https://gist.github.com/JuliaOd/f4d10153e539c026daf7076f087dc937/raw/validdata2dNNPP.h5", validdata)

# Build network
train_network(network_type, model_setup)


text:
using AccuracyAtTop
using AccuracyAtTop.ChainRulesCore
using ChainRulesTestUtils
using Test

@testset "all tests" begin
    include("./thresholds.jl")
    include("./objectives.jl")
end


text:
function fib(n)
  if n <= 2
    1
  else
    fib(n - 1) + fib(n - 2)
  end
end

println(fib(45))

max_stars_repo_name: davidkellis/benchmarks

max_issues_repo_name: davidkellis/benchmarks

max_forks_repo_name: davidkellis/benchmarks


text:
module ReachableControls

using ConcreteStructs: @concrete
using ControlSystems
using IntervalArithmetic
using TaylorModels

export (..), ±

include("as_real.jl")
export as_real

end


text:
# Benders Algorithm for MILP with Sub and Ray Problems
# Version: 5.0
# Author: Edward J. Xu, edxu96@outlook.com
# Date: March 27th, 2019
module BendersMilp_EDXU
    export BendersMilp
    using JuMP
    using GLPKMathProgInterface
    using PrettyTables
    function BendersMilp(; n_x, n_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b,
        mat_a, mat_b, epsilon, timesIterationMax)
        println("-------------------------------------------------------------------------\n",
                "------------------------ 1/4. Begin Optimization ------------------------\n",
                "-------------------------------------------------------------------------\n")
        # Define Master problem
        n_constraint = length(mat_a[:, 1])
        model_mas = Model(solver = GLPKSolverMIP())
        @variable(model_mas, q)
        @variable(model_mas, vec_y[1: n_y], Int)
        @objective(model_mas, Min, (transpose(vec_f) * vec_y + q)[1])
        @constraint(model_mas, vec_y[1: n_y] .<= vec_max_y)
        @constraint(model_mas, vec_y[1: n_y] .>= vec_min_y)

        function solve_master(vec_uBar, opt_cut::Bool)
            if opt_cut
                @constraint(model_mas, (transpose(vec_uBar) * (vec_b - mat_b * vec_y))[1] <= q)
            else  # Add feasible cut Constraints
                @constraint(model_mas, (transpose(vec_uBar) * (vec_b - mat_b * vec_y))[1] <= 0)
            end
            @constraint(model_mas, (transpose(vec_uBar) * (vec_b - mat_b * vec_y))[1] <= q)
            solve(model_mas)
            vec_result_y = getvalue(vec_y)
            return getobjectivevalue(model_mas)
        end

        function solve_sub(vec_uBar, vec_yBar, n_constraint, vec_b, mat_b, mat_a, vec_c)
            model_sub = Model(solver = GLPKSolverLP())
            @variable(model_sub, vec_u[1: n_constraint] >= 0)
            @objective(model_sub, Max, (transpose(vec_b - mat_b * vec_yBar) * vec_u)[1])
            constraintsForDual = @constraint(model_sub, transpose(mat_a) * vec_u .<= vec_c)
            solution_sub = solve(model_sub)
            print("------------------------------ Sub Problem ------------------------------\n")  # , model_sub)
            vec_uBar = getvalue(vec_u)
            if solution_sub == :Optimal
                vec_result_x = zeros(length(vec_c))
                vec_result_x = getdual(constraintsForDual)
                return (true, getobjectivevalue(model_sub), vec_uBar, vec_result_x)
            end
            if solution_sub == :Unbounded
                print("Not solved to optimality because feasible set is unbounded.\n")
                return (false, getobjectivevalue(model_sub), vec_uBar, repeat([NaN], length(vec_c)))
            end
        end

        function solve_ray(vec_uBar, vec_yBar, n_constraint, vec_b, mat_b, mat_a)
            # model_ray = Model(solver = GurobiSolver())
            model_ray = Model(solver = GLPKSolverLP())
            @variable(model_ray, vec_u[1: n_constraint] >= 0)
            @objective(model_ray, Max, 1)
            @constraint(model_ray, (transpose(vec_b - mat_b * vec_yBar) * vec_u)[1] == 1)
            @constraint(model_ray, transpose(mat_a) * vec_u .<= 0)
            solve(model_ray)
            print("------------------------------ Ray Problem ------------------------------\n")  # , model_ray)
            vec_uBar = getvalue(vec_u)
            obj_ray = getobjectivevalue(model_ray)
            return (obj_ray, vec_uBar)
        end

        # Begin Calculation --------------------------------------------------------------------------------------------
        let
            boundUp = Inf
            boundLow = - Inf
            epsilon = 0
            # initial value of master variables
            vec_uBar = zeros(n_constraint, 1)
            vec_yBar = zeros(n_y, 1)
            vec_result_x = length(n_x)
            dict_obj_mas = Dict()
            dict_q = Dict()
            dict_obj_sub = Dict()
            dict_obj_ray = Dict()
            dict_boundUp = Dict()
            dict_boundLow = Dict()
            obj_sub = 0
            timesIteration = 1
            # Must make sure "result_q == obj_sub" in the final iteration
            # while ((boundUp - boundLow > epsilon) && (timesIteration <= timesIterationMax))  !!!
            while ((!((boundUp - boundLow <= epsilon) && ((result_q == obj_sub)))) &&
                (timesIteration <= timesIterationMax))
                (bool_solutionSubModel, obj_sub, vec_uBar, vec_result_x) = solve_sub(vec_uBar, vec_yBar, n_constraint,
                                                                                vec_b, mat_b, mat_a, vec_c)
                if bool_solutionSubModel
                    boundUp = min(boundUp, obj_sub + (transpose(vec_f) * vec_yBar)[1])
                else
                    (obj_ray, vec_uBar) = solve_ray(vec_uBar, vec_yBar, n_constraint, vec_b, mat_b, mat_a)
                end
                obj_mas = solve_master(vec_uBar, bool_solutionSubModel)
                vec_yBar = getvalue(vec_y)
                boundLow = max(boundLow, obj_mas)
                dict_boundUp[timesIteration] = boundUp
                dict_boundLow[timesIteration] = boundLow
                if bool_solutionSubModel
                    dict_obj_mas[timesIteration] = obj_mas
                    dict_obj_sub[timesIteration] = obj_sub
                    result_q = getvalue(q)
                    dict_q[timesIteration] = result_q
                    println("------------------ Result in $(timesIteration)-th Iteration with Sub ",
                            "-------------------\n", "boundUp: $(round(boundUp, digits = 5)), ",
                            "boundLow: $(round(boundLow, digits = 5)), obj_mas: $(round(obj_mas, digits = 5)), ",
                            "q: $result_q, obj_sub: $(round(obj_sub, digits = 5)).")
                else
                    dict_obj_mas[timesIteration] = obj_mas
                    dict_obj_ray[timesIteration] = obj_ray
                    result_q = getvalue(q)
                    dict_q[timesIteration] = result_q
                    println("------------------ Result in $(timesIteration)-th Iteration with Ray ",
                            "-------------------\n", "boundUp: $(round(boundUp, digits = 5)), ",
                            "boundLow: $(round(boundLow, digits = 5)), obj_mas: $(round(obj_mas, digits = 5)), ",
                            "q: $result_q, obj_ray: $(round(obj_ray, digits = 5)).")
                end
                timesIteration += 1
            end  # -----------------------------------------------------------------------------------------------------
            println("obj_mas: $(getobjectivevalue(model_mas))")
            println("----------------------------- Master Problem ----------------------------\n")
            println(model_mas)
            println("-------------------------------------------------------------------------\n",
                    "------------------------------ 2/4. Result ------------------------------\n",
                    "-------------------------------------------------------------------------")
            println("boundUp: $(round(boundUp, digits = 5)), boundLow: $(round(boundLow, digits = 5)), ",
                    "difference: $(round(boundUp - boundLow, digits = 5))")
            println("vec_x: $vec_result_x")
            vec_result_y = getvalue(vec_y)
            result_q = getvalue(q)
            println("vec_y: $vec_result_y")
            println("result_q: $result_q")
            println("-------------------------------------------------------------------------\n",
                    "------------------------- 3/4. Iteration Result -------------------------\n",
                    "-------------------------------------------------------------------------")
            # Initialize
            seq_timesIteration = collect(1: (timesIteration - 1))
            vec_boundUp = zeros(timesIteration - 1)
            vec_boundLow = zeros(timesIteration - 1)
            vec_obj_subRay = zeros(timesIteration - 1)
            vec_obj_mas = zeros(timesIteration - 1)
            vec_q = zeros(timesIteration - 1)
            vec_type = repeat(["ray"], (timesIteration - 1))
            #
            for i = 1: (timesIteration - 1)
                vec_obj_mas[i] = round(dict_obj_mas[i], digits = 5)
                vec_boundUp[i] = round(dict_boundUp[i], digits = 5)
                vec_boundLow[i] = round(dict_boundLow[i], digits = 5)
                vec_q[i] = round(dict_q[i], digits = 5)
                if haskey(dict_obj_sub, i)
                    vec_type[i] = "sub"
                    vec_obj_subRay[i] = round(dict_obj_sub[i], digits = 5)
                else
                    vec_obj_subRay[i] = round(dict_obj_ray[i], digits = 5)
                end
            end
            table_iterationResult = hcat(seq_timesIteration, vec_boundUp, vec_boundLow,
                                         vec_obj_mas, vec_q, vec_type, vec_obj_subRay)
            pretty_table(table_iterationResult,
                         ["Seq", "boundUp", "boundLow", "obj_mas", "q", "sub/ray", "obj_sub/ray"],
                         compact; alignment=:l)
        end
        println("-------------------------------------------------------------------------\n",
                "-------------------------- 4/4. Nominal Ending --------------------------\n",
                "-------------------------------------------------------------------------\n")
    end
end


text:
"""
    Node{T}

    Abstract type for graph node object
"""
abstract type Node{T} end


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

using Markdown
using InteractiveUtils

# ╔═╡ 3a2812d8-3770-11eb-34b0-7fd565a38a5f
import Pkg; Pkg.add("DataStructures");using DataStructures

# ╔═╡ 1a5162ec-376c-11eb-2960-5d2de63d1513
file = readlines("day_5.txt")

# ╔═╡ 63d8cc20-376c-11eb-0a15-e74649299804
function get_id(not_zone_or_group)
	upper = map(x-> (x=='B'|| x=='R'), collect(not_zone_or_group))

	return sum([up*2^(10-index) for (index,up) in enumerate(upper)])
end

# ╔═╡ ad648972-376e-11eb-07d7-1f77ba810dad
passager_id = [get_id(line) for line in file]

# ╔═╡ c7ed4bda-376e-11eb-22d9-e14675f7f9e2
maximum(passager_id)

# ╔═╡ fcd1f052-376f-11eb-2e54-611deedd7c74
md"# Question 2"

# ╔═╡ 04515e12-3770-11eb-2212-c55609da65f4
occupy = DefaultDict(false)

# ╔═╡ 58485228-3770-11eb-188a-e3f1060f6f89
for id in passager_id
	occupy[id] = true
end

# ╔═╡ 832f9fd6-3771-11eb-10cb-dd40d51e2e4c
sum([(occupy[i]==false)&&occupy[i-1]&&occupy[i+1]  for i in 8:1023-8])

# ╔═╡ bc646fec-3772-11eb-07a7-a963d0f51465
begin
	assent = -1
for i in 8:1023-8
	if (occupy[i]==false) && occupy[i-1] && occupy[i+1]
		assent = i
	end
end
end

# ╔═╡ e59c62c0-3772-11eb-3887-c52aa4dcf981
assent

# ╔═╡ Cell order:
# ╠═3a2812d8-3770-11eb-34b0-7fd565a38a5f
# ╠═1a5162ec-376c-11eb-2960-5d2de63d1513
# ╠═63d8cc20-376c-11eb-0a15-e74649299804
# ╠═ad648972-376e-11eb-07d7-1f77ba810dad
# ╠═c7ed4bda-376e-11eb-22d9-e14675f7f9e2
# ╠═fcd1f052-376f-11eb-2e54-611deedd7c74
# ╠═04515e12-3770-11eb-2212-c55609da65f4
# ╠═58485228-3770-11eb-188a-e3f1060f6f89
# ╠═832f9fd6-3771-11eb-10cb-dd40d51e2e4c
# ╠═bc646fec-3772-11eb-07a7-a963d0f51465
# ╠═e59c62c0-3772-11eb-3887-c52aa4dcf981


text:
using Test

@testset "Dorothy" begin

	include("graphs.jl")
	include("geometry.jl")
	include("pbc.jl")

	include("models.jl")
	include("atoms.jl")
	include("hierarchies.jl")
	include("selections.jl")

	@testset "Formats" begin
		include("formats/gromos87.jl")
		include("formats/pdb.jl")
		include("formats/trr.jl")
		include("formats/xtc.jl")
		include("formats/ndx.jl")
	end

	include("topology.jl")
	include("ss.jl")

end # @testset


text:
module InterpolationFunctions

	using Dates
	import Statistics
	export interpolatedMax, interpolate, interpolatedSum, addArraysShiftedInterpolated, medianfilter, averageSamples, smooth

	function interpolatedMax(discreteMax, values)
	  inta=values[discreteMax-1]
	  intb=values[discreteMax]
	  intc=values[discreteMax+1]
	  minAC = min(inta, intc)
	  max = discreteMax + (intc - inta)/(inta + intb + intc - 3*minAC) # Tricky: -0.5 if a=b, +0.5 if b=c, +0.0 if a=c, interpol in between
	return max
	end

	function interpolate(x::AbstractFloat, xAxis, yAxis)
	    try
	      if x <= xAxis[1]
		return yAxis[1]
	      elseif x >= xAxis[end]
		return yAxis[end]
	      end
	      indexLow = searchsortedfirst(xAxis,x) - 1
	      fraction = (x - xAxis[indexLow]) / (xAxis[indexLow+1] - xAxis[indexLow])
	      return fraction * yAxis[indexLow+1] + (1-fraction) * yAxis[indexLow]
	    catch
		println("Interpolation Error")
	    end
	end

	function interpolate(x::DateTime, xAxis, yAxis)
	    x=Dates.datetime2unix(x)
	    xAxis = Dates.datetime2unix.(xAxis)
	    try
	      if x <= xAxis[1]
		return yAxis[1]
	      elseif x >= xAxis[end]
		return yAxis[end]
	      end
	      indexLow = searchsortedfirst(xAxis,x) - 1
	      fraction = (x - xAxis[indexLow]) / (xAxis[indexLow+1] - xAxis[indexLow])
	      return fraction * yAxis[indexLow+1] + (1-fraction) * yAxis[indexLow]
	    catch
		println("Interpolation Error")
	    end
	end

	function interpolate(x::AbstractFloat, yAxis)
	    try
	      if x <= 1
		return yAxis[1]
	      elseif x >= length(yAxis)
		return yAxis[end]
	      end
	      indexLow = Int64(floor(x))
	      fraction = (x - indexLow)
	      return fraction * yAxis[indexLow+1] + (1-fraction) * yAxis[indexLow]
	  catch
	      println("Interpolation Error")
	  end
	end

	function interpolate(x::AbstractArray, xAxis, yAxis)
	  y = Array{typeof(yAxis[1])}(undef,length(x))
	  Threads.@threads for i = 1 : length(x)
	    y[i] = InterpolationFunctions.interpolate(x[i],xAxis,yAxis)
	  end
	  return y #convert(Array,y)
	end

	function interpolateMatrix_1D(x::AbstractArray, xAxis, yMatrix)
		yMatrix_intpol = Array{typeof(yMatrix[1])}(undef,length(x), size(yMatrix)[2])
	  	for i in 1:size(yMatrix)[2]
			yMatrix_intpol[:,i] = interpolate(x, xAxis, yMatrix[:,i])
		end
		return yMatrix_intpol
	end

	function interpolate(x::AbstractArray, yAxis)
	  y = Array{typeof(yAxis[1])}(undef,length(x))
	  Threads.@threads for i = 1 : length(x)
	    y[i] = InterpolationFunctions.interpolate(x[i],yAxis)
	  end
	  return y #convert(Array,y)
	end

	function interpolatedSum(startX::AbstractFloat, endX::AbstractFloat, xAxis, yAxis)
	  firstCompleteIndex = searchsortedfirst(xAxis, startX)
	  lastCompleteIndex = searchsortedfirst(xAxis, endX) - 1
	  if (firstCompleteIndex > 0) && (firstCompleteIndex < length(yAxis)) && (lastCompleteIndex > 0) && (lastCompleteIndex < length(yAxis))

	    fractionStart = xAxis[firstCompleteIndex] - startX
	    firstFractionValueContribution = fractionStart * yAxis[firstCompleteIndex-1]

	    fractionEnd = endX - xAxis[lastCompleteIndex]
	    lastFractionValueContribution = fractionEnd * yAxis[lastCompleteIndex + 1]

	    if firstCompleteIndex == lastCompleteIndex
	      return (firstFractionValueContribution + lastFractionValueContribution)
	    end
	    return (sum(yAxis[firstCompleteIndex:lastCompleteIndex] + firstFractionValueContribution + lastFractionValueContribution))
	  end
	  return 0
	end

	function interpolatedSum(startIndexExact::AbstractFloat, endIndexExact::AbstractFloat, yAxis)
	  subIdxStart::Int64 = ceil(startIndexExact)
	  subIdxStartRoundError = subIdxStart - startIndexExact

	  subIdxEnd::Int64 = floor(endIndexExact)
	  subIdxEndRoundError = endIndexExact - subIdxEnd
	  ret = 0
	  if (subIdxStart>1) && (subIdxEnd+1 < length(yAxis))
	    ret = sum(yAxis[subIdxStart:subIdxEnd]) + yAxis[subIdxStart-1]*subIdxStartRoundError + yAxis[subIdxEnd+1]*subIdxEndRoundError
	  end
	  return ret
	end

	function addArraysShiftedInterpolated(destinationArray::Array, sourceArray::Array, indexShift::Number)
	    if ceil(indexShift) == indexShift
		lowContrib = 1
	    else
		lowContrib = ceil(indexShift) - indexShift
	    end
	    highContrib = 1-lowContrib
	    minIdx = Int64(floor(indexShift)+1)
	    maxIdx = Int64(minIdx + length(sourceArray)-1)
	    if minIdx > length(destinationArray)
		return destinationArray
	    end
	    if maxIdx > length(destinationArray)
		maxIdx = length(destinationArray)
	    end
	    destinationArray[minIdx:maxIdx] += lowContrib*sourceArray[1:maxIdx-minIdx+1]
	    if maxIdx > length(destinationArray)-1
		maxIdx = length(destinationArray)-1
	    end
	    destinationArray[minIdx+1:maxIdx+1] += highContrib*sourceArray[1:maxIdx-minIdx+1]
	    return destinationArray
	end

	function medianfilter(v,ws)
	  [median(v[i:(i+ws-1)]) for i=1:(length(v)-ws+1)]
	end

	function nanmean(x::AbstractArray)
		if length(filter(!isnan, x)) > 0
			return Statistics.mean(filter(!isnan,x))
		elseif length(filter(!isnan, x)) == 0
			return NaN
		end
	end

	function nanmean(x::AbstractArray,d)
		return mapslices(nanmean,x;dims = d)
	end

	# Averaging multidimensional array in one dimension
	# adapted from https://julialang.org/blog/2016/02/iteration
	function averageSamples(data, averagePoints; dim=1, ignoreNaNs = false)
	    if averagePoints > 1
		dt=false
		if typeof(data[1]) == DateTime
		    data = Dates.datetime2unix.(data)
		    dt=true
		end
		len = Int64(floor(size(data,dim)/averagePoints))
		sz = [size(data)...]
		sz[[dim...]] .= len
		averaged = Array{eltype(data)}(undef, sz...)
		Rpre = CartesianIndices(size(data)[1:dim-1])
		Rpost = CartesianIndices(size(data)[dim+1:end])
		for Ipost in Rpost
		  for Ipre in Rpre
		    for i=1:len
		      	if ignoreNaNs
				println(data[Ipre, ((i-1)*averagePoints+1):(i*averagePoints), Ipost])
				averaged[Ipre, i, Ipost] = nanmean(data[Ipre, ((i-1)*averagePoints+1):(i*averagePoints), Ipost])
		      	else
		      		averaged[Ipre, i, Ipost] = Statistics.mean(data[Ipre, ((i-1)*averagePoints+1):(i*averagePoints), Ipost])
			end
		    end
		  end
		end
		if dt
		    averaged = Dates.unix2datetime.(averaged)
		end
		return averaged
	    else
		return data
	    end
	end

	function smooth(x, α, dim::Integer=1)
	    s = similar(x)
	    Rpre = CartesianIndices(size(x)[1:dim-1])
	    Rpost = CartesianIndices(size(x)[dim+1:end])
	    _smooth!(s, x, α, Rpre, size(x, dim), Rpost)
	end

	function _smooth!(s, x, α, Rpre, n, Rpost)
	    for Ipost in Rpost
	        # Initialize the first value along the filtered dimension
	        for Ipre in Rpre
	            s[Ipre, 1, Ipost] = x[Ipre, 1, Ipost]
	        end
	        # Handle all other entries
	        for i = 2:n
	            for Ipre in Rpre
	                s[Ipre, i, Ipost] = α*x[Ipre, i, Ipost] + (1-α)*s[Ipre, i-1, Ipost]
	            end
	        end
	    end
	    s
	end

end


text: autocor(DataArray([1, 2, 3, 4, 5]))


text:
# a one-level decision tree
mutable struct OneTree
    atr_name :: String                  # name of the attribute for which the tree is created
    nodes :: AbstractVector{OneNode}    # the rules within the decision tree
    err_rate :: Real                    # total error rate of the tree; computed from the error rates of its rules
end

"""
    get_best_tree(train_data::AbstractDataFrame, observations::AbstractVector) -> OneTree

Finds the one-level decision tree with the lowest error rate (a `OneTree`).
`train_data` contains the predictor attributes and `observations` the corresponding
observations.
"""
function get_best_tree(train_data::AbstractDataFrame, observations::AbstractVector)
    trees = all_trees(train_data, observations)
    return(trees[argmin(trees)])
end

# overload the < and == operators, so that `argmin` works
Base.:isless(t1::OneTree, t2::OneTree) = (t1.err_rate < t2.err_rate)
Base.:isequal(t1::OneTree, t2::OneTree) = (t1.err_rate == t2.err_rate)

# generate for each attribute in `train_data` a decision tree
function all_trees(train_data::AbstractDataFrame, observations::AbstractVector)
    atr_names = names(train_data)       # all attribute names
    uniq_obs = unique(observations)     # possible observations
    inst_count = length(observations)   # number of instances in `train_data`

    trees = Vector{OneTree}(undef, length(atr_names)) # for the results

    # make a tree for each attribute
    for i in eachindex(atr_names)
        nodes = get_nodes(train_data[:, atr_names[i]], observations, uniq_obs)
        errsum = 0          # sum up all `err_count`s for this attribute
        for n in nodes
            errsum += n.err_count
        end
        trees[i] = OneTree(atr_names[i], nodes, errsum / inst_count)
    end
    return(trees)
end

# print a 'OneTree' nicely formatted
function Base.show(io::IO, ::MIME"text/plain", tree::OneTree)
    println(io, "OneTree: $(tree.atr_name)")
    for n in tree.nodes
        println(io, "    $(n)")
    end
end

max_stars_repo_name: roland-KA/OneRule

max_issues_repo_name: roland-KA/OneRule

max_forks_repo_name: roland-KA/OneRule


text:
using NonLinearReactionAdvectionDiffusionWithFrontData;
using NonLinearReactionAdvectionDiffusionWithFrontData: find_f_zeros;

@testset "Поиск нулей функции                              " begin

    g(x) = (x - 1.1)*(x + 1.1) # функция

    X1 = collect(range(-2, 0, length=20));
    X2 = collect(range(0, 2, length=20));

    @test isapprox( find_f_zeros( g.(X1), X1), -1.1, atol=1e-4)
    @test isapprox( find_f_zeros( g.(X2), X2), 1.1, atol=1e-4)

    h(x) = (x - 1) * (x - 2) * (x - 3)

    # Введем недостаточно густую сетку.
    # Выбранное окно функции для аппроксимации будет немонотонным.
    X = collect( range(-4, 4, length = 50))
    @test_throws ArgumentError find_f_zeros( h.(X), X)

    X = collect( range(-4, 4, length = 100))
    @test isapprox( find_f_zeros( h.(X), X), 1.0, atol=1e-4)

    # Не работает на неотрицательных функция
    j(x) = sin(x)*sin(2x)
    X = collect( range(-0.5, 0.5, length = 100))
    @test_throws ArgumentError find_f_zeros( j.(X), X)

    # Подадим немонотонные значения функции: см первый массив
    @test_throws ArgumentError find_f_zeros( [-1, 1, 3, 2, 4, 5, 6, 8, 7], [-1, 1, 2, 3, 4, 5, 6, 7, 8] )

end

@testset "Поиск вырожденных корней                         " begin

    a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ = NonLinearReactionAdvectionDiffusionWithFrontData.dparams();

    # Сетку нужно подавать в обычном виде
    @test_throws ArgumentError phidetermination(qₙ, ulₘ, reverse(Xₙ), N, Tₘ, M);

end


text:
using Base.Test
using QuantumOptics

const ode = QuantumOptics.ode_dopri.ode

const ζ = 0.5
const ω₀ = 10.0

const y₀ = Float64[0., sqrt(1-ζ^2)*ω₀]
const A = 1
const ϕ = 0

function f(t::Float64)
    α = sqrt(1-ζ^2)*ω₀
    x = A*exp(-ζ*ω₀*t)*sin(α*t + ϕ)
    p = A*exp(-ζ*ω₀*t)*(-ζ*ω₀*sin(α*t + ϕ) + α*cos(α*t + ϕ))
    return [x,p]
end

function df(t::Float64, y::Vector{Float64}, dy::Vector{Float64})
    dy[1] = y[2]
    dy[2] = -2*ζ*ω₀*y[2] - ω₀^2*y[1]
    return nothing
end

# Test absolute and relative tolerances
T = [0.:0.5:5.;]

for reltol=[1e-3,1e-4,1e-5,1e-6,1e-7,1e-8,1e-9]
    tout, y = ode(df, T, y₀, reltol=reltol, abstol=0.)
    for i=2:length(T)
        # println("T = $(T[i])")
        # println("Target Reltol: $(reltol*T[i]); Real reltol: ", norm(y[i]-f(T[i]))/norm(y[i]))
        @test T[i]*reltol*20 > norm(y[i]-f(T[i]))/norm(y[i]) > T[i]*reltol/20
    end
end

for abstol=[1e-3,1e-4,1e-5,1e-6,1e-7,1e-8,1e-9]
    tout, y = ode(df, T, y₀, reltol=0., abstol=abstol)
    for i=2:length(T)
        # println("T = $(T[i])")
        # println("Target abstol: $(abstol*T[i]); Real abstol: ", norm(y[i]-f(T[i])))
        @test abstol*10 > norm(y[i]-f(T[i])) > abstol/100
    end
end

# Test output of intermediate steps
T = [0.,10.]
tout, yout = ode(df, T, y₀; display_intermediatesteps=true)
@test length(tout)>2
maxstep = maximum(tout[2:end]-tout[1:end-1])

# Test hmax stepszie control
hmax = maxstep/10
tout, yout = ode(df, T, y₀; display_intermediatesteps=true, hmax=hmax)
maxstep2 = maximum(tout[2:end]-tout[1:end-1])
@test (maxstep2-hmax)<1e-12

# Test hmin
@test_throws ErrorException ode(df, T, y₀; display_intermediatesteps=true, hmin=0.4)

# Test h0
h0 = 1e-5
tout, yout = ode(df, T, y₀; display_intermediatesteps=true, h0=h0)
@test abs(h0-(tout[2]-tout[1]))<1e-14
h0 = 0.4
tout, yout = ode(df, T, y₀; display_intermediatesteps=true, h0=h0)
@test abs(h0-(tout[2]-tout[1]))>0.3

# Test fout
tout_ = Float64[]
yout_ = Vector{Float64}[]
function fout(t, y)
    push!(tout_, t)
    push!(yout_, deepcopy(y))
end
T = [0.,0.5,1.0]
result = ode(df, T, y₀, fout)
@test result==nothing
tout, yout = ode(df, T, y₀)

@test tout==tout_
@test yout==yout_


text:
# functions related to beta distributions

import .RFunctions:
    betalogpdf,
    betacdf,
    betaccdf,
    betalogcdf,
    betalogccdf,
    betainvcdf,
    betainvccdf,
    betainvlogcdf,
    betainvlogccdf

# pdf
betapdf(α::Float64, β::Float64, x::Float64) = exp(betalogpdf(α, β, x))
betapdf(α::Real, β::Real, x::Real) = betapdf(Float64(α), Float64(β), Float64(x))


text:
"""
band(x, ylower, yupper; kwargs...)
band(lower, upper; kwargs...)

Plots a band from `ylower` to `yupper` along `x`.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Band, lowerpoints, upperpoints) do scene
    Attributes(;
        default_theme(scene, Mesh)...,
        colorrange = automatic,
        color = RGBAf0(0.0,0,0,0.2)
    )
end

convert_arguments(::Type{<: Band}, x, ylower, yupper) = (Point2f0.(x, ylower), Point2f0.(x, yupper))

function band_connect(n)
    ns = 1:n-1
    ns2 = n+1:2n-1
    [GLTriangleFace.(ns, ns .+ 1, ns2); GLTriangleFace.(ns .+ 1, ns2 .+ 1, ns2)]
end

function plot!(plot::Band)
    @extract plot (lowerpoints, upperpoints)
    @lift(@assert length($lowerpoints) == length($upperpoints) "length of lower band is not equal to length of upper band!")
    coordinates = @lift([$lowerpoints; $upperpoints])
    connectivity = lift(x -> band_connect(length(x)), plot[1])
    mesh!(plot, coordinates, connectivity;
        color = plot[:color], colormap = plot[:colormap],
        colorrange = plot[:colorrange],
        shading = false, visible = plot[:visible]
    )
end

function fill_view(x, y1, y2, where::Nothing)
    x, y1, y2
  end
  function fill_view(x, y1, y2, where::Function)
    fill_view(x, y1, y2, where.(x, y1, y2))
  end
  function fill_view(x, y1, y2, bools::AbstractVector{<: Union{Integer, Bool}})
    view(x, bools), view(y1, bools), view(y2, bools)
  end

  """
      fill_between!(x, y1, y2; where = nothing, scene = current_scene(), kw_args...)

  fill the section between 2 lines with the condition `where`
  """
  function fill_between!(x, y1, y2; where = nothing, scene = current_scene(), kw_args...)
    xv, ylow, yhigh = fill_view(x, y1, y2, where)
    band!(scene, xv, ylow, yhigh; kw_args...)
  end

  export fill_between!


text:
using DDR2import

function get_arefile(arefile::String)
    are = DDR2import.Are.read(arefile)
    sectors = Vector{Vector{Array{Float64, 2}}}()
    for sector in are
        append!(sectors, [[sector[2].points]])
    end
    return sectors
end

function add_sectors!(sectors::String, projparam=Equirectangular();
    linecolor=:grey, linewidth=0.01, linealpha=0.33, kwargs...)
    sectors_array = get_arefile(sectors)
    add_polygon!(sectors_array, projparam; linecolor=linecolor,
    linewidth=linewidth, linealpha=linealpha, kwargs...)
end

add_sectors!(sectors::String, mapsettings::MapSettings; kwargs...) =
add_sectors!(sectors, mapsettings.projection; kwargs...)

function add_freeflightairspaces!(sectors::String, projparam=Equirectangular();
    linecolor=:orange, linewidth=0.01, linealpha=0.2, fillcolor=:orange,
    fillalpha=0.1, fillrange=0, kwargs...)
    sectors_array = get_arefile(sectors)
    add_polygon!(sectors_array, projparam; linecolor=linecolor,
    linewidth=linewidth, linealpha=linealpha, fillcolor=fillcolor,
    fillalpha=fillalpha, fillrange=fillrange, kwargs...)
end

add_freeflightairspaces!(sectors::String, mapsettings::MapSettings; kwargs...) =
add_freeflightairspaces!(sectors, mapsettings.projection; kwargs...)


text:
using InternedStrings

using Test, Random, Base.GC

addr_eq(a::String, b::String) = pointer(a) === pointer(b)
addr_eq(a, b) = objectid(a) === objectid(b)

@testset "All kinds of types" begin include("all_kinds_of_types.jl") end
@testset "String macro"       begin include("string_macro.jl") end
@testset "Core Functionality" begin include("corefunctionality.jl") end


text:
abstract type AbstractGBType end
abstract type AbstractDescriptor end
abstract type AbstractOp end
abstract type AbstractUnaryOp <: AbstractOp end
abstract type AbstractBinaryOp <: AbstractOp end
abstract type AbstractSelectOp <: AbstractOp end
abstract type AbstractMonoid <: AbstractOp end
abstract type AbstractSemiring <: AbstractOp end
abstract type AbstractTypedOp{Z} end


text:
const USE_GPU = false
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using ImplicitGlobalGrid, Plots, Printf, Statistics
import MPI
# pyplot()
do_viz = false

maximum_g(A) = (max_l  = maximum(A); MPI.Allreduce(max_l,  MPI.MAX, MPI.COMM_WORLD))

@parallel function check_Err!(ErrH::Data.Array, H::Data.Array)
    @all(ErrH) = @all(ErrH) - @all(H)
    return
end

@parallel function compute_Err!(ErrH::Data.Array, H::Data.Array)
    @all(ErrH) = @all(H)
    return
end

@parallel function compute_flux!(qx::Data.Array, qy::Data.Array, H::Data.Array, n::Data.Number, dx::Data.Number, dy::Data.Number)
    @all(qx) = -@av_xi(H)^n*@d_xi(H)/dx
    @all(qy) = -@av_yi(H)^n*@d_yi(H)/dy
    return
end

@parallel function compute_ResH!(dtau::Data.Array, ResH::Data.Array, H::Data.Array, qx::Data.Array, qy::Data.Array, b::Data.Array, n::Data.Number, dmp::Data.Number, dx::Data.Number, dy::Data.Number)
    @all(dtau) = min(dx*dx, dy*dy)/(1.0 + @inn(H)^n)/4.1./4.0
    @all(ResH) = (-@d_xa(qx)/dx  - @d_ya(qy)/dy + @inn(b)) + dmp*@all(ResH)
    return
end

@parallel function update_H!(H::Data.Array, dtau::Data.Array, ResH::Data.Array)
    @inn(H) = max(0.0, @inn(H) + @all(dtau)*@all(ResH))
    return
end

@views function sia_2D_xpu()
    # Physics
    lx, ly = 10.0, 10.0
    n      = 3.0
    # Numerics
    nx, ny = 64, 64
    niter  = 10000
    nout   = 100
    dmp    = 0.96
    ε      = 1e-8
    # Derived numerics
    me, dims, nprocs, coords, comm = init_global_grid(nx, ny, 1) # MPI initialisation
    select_device()                                              # select one GPU per MPI local rank (if >1 GPU per node)
    dx, dy = lx/nx_g(), ly/ny_g()                                # cell sizes
    # Array allocation
    H      = @zeros(nx  ,ny  )
    qx     = @zeros(nx-1,ny-2)
    qy     = @zeros(nx-2,ny-1)
    dtau   = @zeros(nx-2,ny-2)
    ResH   = @zeros(nx-2,ny-2)
    ErrH   = @zeros(nx  ,ny  )
    # Initial condition
    b      = 0.5*ones(nx,ny)
    H     .= Data.Array([exp(-(x_g(ix,dx,H)+dx/2 -0.5*lx)^2 -(y_g(iy,dy,H)+dy/2 -0.5*ly)^2) for ix=1:size(H,1), iy=1:size(H,2)])
    rad    = [(x_g(ix,dx,H)+dx/2 -0.5*lx)^2 + (y_g(iy,dy,H)+dy/2 -0.5*ly)^2 for ix=1:size(H,1), iy=1:size(H,2)]; b[rad.>lx/4] .= -0.5
    b      = Data.Array(b) # move data back to XPU
    # Preparation of visualisation
    if do_viz
        if me==0
            ENV["GKSwstype"]="nul"; if isdir("viz2D_xpu_out")==false mkdir("viz2D_xpu_out") end; loadpath = "./viz2D_xpu_out/"; anim = Animation(loadpath,String[])
            println("Animation directory: $(anim.dir)")
        end
        nx_v, ny_v = (nx-2)*dims[1], (ny-2)*dims[2]
        if (nx_v*ny_v*sizeof(Data.Number) > 0.8*Sys.free_memory()) error("Not enough memory for visualization.") end
        H_v   = zeros(nx_v, ny_v) # global array for visu
        H_inn = zeros(nx-2, ny-2) # no halo local array for visu
        Xi_g, Yi_g = LinRange(dx+dx/2, lx-dx-dx/2, nx_v), LinRange(dy+dy/2, ly-dy-dy/2, ny_v) # inner points only
    end
    # Time loop
    for iter = 1:niter
        if (iter==11) global t0 = Base.time() end
        @parallel compute_Err!(ErrH, H)
        @parallel compute_flux!(qx, qy, H, n, dx, dy)
        @parallel compute_ResH!(dtau, ResH, H, qx, qy, b, n, dmp, dx, dy)
        @hide_communication (8, 4) begin # communication/computation overlap
            @parallel update_H!(H, dtau, ResH)
            update_halo!(H)
        end
        @parallel check_Err!(ErrH, H)
        # Visualise
        if mod(iter,nout)==0 && do_viz
            H_inn .= H[2:end-1,2:end-1]; gather!(H_inn, H_v)
            if me==0
                p1 = heatmap(Xi_g, Yi_g, H_v', xlabel="lx", ylabel="ly", title="shallow ice, iter=$iter", aspect_ratio=1, xlims=(Xi_g[1], Xi_g[end]), ylims=(Yi_g[1], Yi_g[end]), c=:viridis)
                p2 = plot(Xi_g, H_v[:,Int(round(ny_g()/2))], xlabel="lx", ylabel="height", xlims=(Xi_g[1], Xi_g[end]), ylims=(0., 1.2), legend=false, framestyle=:box, aspect_ratio=2.5)
                l = @layout [ a{0.8h} ; b{0.2h} ]
                plot(p1, p2, layout=l ); frame(anim)
            end
        end
        if mod(iter,nout)==0 maxErr=maximum_g(abs.(ErrH)); if (me==0) @printf("iter=%d, max(err)=%1.2e \n", iter, maxErr) end; if maxErr<ε global itg=iter; break; end end
    end
    time_s = (Base.time()-t0)
    if (me==0) @printf("Time = %1.4e s, T_eff = %1.2f GB/s \n", time_s, round((2/1e9*nx*ny*sizeof(lx))/(time_s/(itg-10)), sigdigits=2)) end
    if (do_viz && me==0) gif(anim, "sia2D_multixpu.gif", fps = 5)  end
    finalize_global_grid()
    return
end

sia_2D_xpu()

max_stars_repo_name: luraess/geo-hpc-course

max_issues_repo_name: luraess/geo-hpc-course

max_forks_repo_name: luraess/geo-hpc-course


text:
# definitions of concrete subtypes of AbstractParamType and AbstractSolutionData

"""
  Subtype of [`AbstractParamType`](@ref).

  **Static Parameters**:

   * Tsol
   * Tres
   *  Tdim

  This is a container passed to all low level function, useful for storing
  miscellaneous parameters or constants
"""
mutable struct ParamType{Tsol, Tres, Tdim} <: AbstractParamType{Tdim}
  LFalpha::Float64  # alpha for the Lax-Friedrich flux
  alpha_x::Float64
  alpha_y::Float64
  alpha_z::Float64
  sin_amplitude::Complex128
  omega::Complex128

  qL_s::Array{Tsol, 1}  # solution vector for a solution grid element
  qR_s::Array{Tsol, 1}  # solution vector for a solution grid element
  qL_f::Array{Tsol, 1}  # solution vector for flux grid element
  qR_f::Array{Tsol, 1}  # solution vector for flux grid element
  resL_s::Array{Tres, 1}  # residual for solution grid element
  resR_s::Array{Tres, 1}  # residual for solution grid element
  resL_f::Array{Tsol, 1}  # residual for a flux grid element
  resR_f::Array{Tsol, 1}
  f::BufferedIO
  x_design::Array{Tres, 1}  # design variables
  time::Timings
  #=
  # timings
  t_volume::Float64  # time for volume integrals
  t_face::Float64 # time for surface integrals (interior)
  t_source::Float64  # time spent doing source term
  t_sharedface::Float64  # time for shared face integrals
  t_bndry::Float64  # time spent doing boundary integrals
  t_send::Float64  # time spent sending data
  t_wait::Float64  # time spent in MPI_Wait
  t_allreduce::Float64 # time spent in allreduce
  t_jacobian::Float64  # time spent computing jacobian
  t_solve::Float64  # linear solve time
  t_barrier::Float64  # time spent in MPI_Barrier
  t_barrier2::Float64
  t_barrier3::Float64
  t_barriers::Array{Float64, 1}
  =#
  function ParamType{Tsol, Tres, Tdim}(mesh, sbp, opts) where {Tsol, Tres, Tdim}
    LFalpha = opts["LFalpha"]
    myrank = mesh.myrank
    if DB_LEVEL >= 1
      _f = open("log_$myrank.dat", "w")
      f = BufferedIO(_f)
    else
      f = BufferedIO()  # create a dummy IOStream
    end

    advection_velocity = opts["advection_velocity"]
    alpha_x = advection_velocity[1]
    alpha_y = advection_velocity[2]
    if Tdim == 3
      alpha_z = advection_velocity[3]
    else
      alpha_z = 1.0
    end

    numNodesPerElement_s = mesh.numNodesPerElement
    if opts["use_staggered_grid"]
      numNodesPerElement_f = mesh.mesh2.numNodesPerElement
    else
      numNodesPerElement_f = numNodesPerElement_s
    end

    qL_s = zeros(Tsol, numNodesPerElement_s)
    qR_s = zeros(Tsol, numNodesPerElement_s)
    qL_f = zeros(Tsol, numNodesPerElement_f)
    qR_f = zeros(Tsol, numNodesPerElement_f)

    resL_s = zeros(Tsol, numNodesPerElement_s)
    resR_s = zeros(Tsol, numNodesPerElement_s)
    resL_f = zeros(Tsol, numNodesPerElement_f)
    resR_f = zeros(Tsol, numNodesPerElement_f)

    # needed for the runtype=660 (CN uadj) objective
    sin_amplitude = 2.0
    omega = 1.0

    x_design = zeros(Tres, 0)  # TODO: get proper size information

    t = Timings()

    return new(LFalpha, alpha_x, alpha_y, alpha_z,
               sin_amplitude, omega,
               qL_s, qR_s, qL_f, qR_f, resL_s, resR_s, resL_f, resR_f,
               f, x_design, t)
  end
end

"""
  Convenient alias for all 2D ParamTypes
"""
ParamType2{Tsol, Tres} =  ParamType{Tsol, Tres, 2}

"""
  Convenient alias for all 3D ParamTypes
"""
ParamType3{Tsol, Tres} =  ParamType{Tsol, Tres, 3}

"""
  All ParamTypes, without static parameters specified
"""
const ParamTypes = Union{ParamType2, ParamType3}

"""
### AdvectionEquationMod.AdvectionData_

  This type is an implementation of the abstract AdvectionData.  It is
  parameterized by Tsol, the datatype of the solution variables and Tmsh,
  the datatype of the mesh variables.
  Tres is the 'maximum' type of Tsol and Tmsh.
  Tdim is the dimensionality of the equation being solve (2 or 3).

  This type is (ultimately) a subtype of [`AbstractSolutionData`](@ref)
  and contains all the required fields.
"""
mutable struct AdvectionData_{Tsol, Tres, Tdim, Tmsh} <: AdvectionData{Tsol, Tres, Tdim}

  # params::ParamType{Tdim}
  params::ParamType{Tsol, Tres, Tdim}

  comm::MPI.Comm
  commsize::Int
  myrank::Int

  t::Float64
  res_type::DataType  # type of res
  q::Array{Tsol, 3}
  q_bar::Array{Tres, 3}  # for testing
  q_face::Array{Tsol, 4}  # store solution values interpolated to faces
  aux_vars::Array{Tres, 3}  # storage for auxiliary variables
  flux_parametric::Array{Tsol,4}  # flux in xi direction
  flux_face::Array{Tres, 3}  # flux for each interface, scaled by jacobian
  res::Array{Tres, 3}      # result of computation
  res_vec::Array{Tres, 1}  # result of computation in vector form
  res_edge::Array{Tres, 4} # edge based residual storage
  res_bar::Array{Tres, 3}
  q_vec::Array{Tres,1}     # initial condition in vector form
  q_bndry::Array{Tsol, 3}  # store solution variables interpolated to
                          # the boundaries with boundary conditions
  shared_data::Array{SharedFaceData{Tsol}, 1}  # MPI data, including send and receive
  shared_data_bar::Array{SharedFaceData{Tres}, 1}

                                         # buffers
#  q_face_send::Array{Array{Tsol, 3}, 1}  # send buffers for sending q values
                                         # to other processes
#  q_face_recv::Array{Array{Tsol, 3}, 1}  # recieve buffers for q values
  flux_sharedface::Array{Array{Tres, 3}, 1}  # hold shared face flux
  bndryflux::Array{Tsol, 3}  # boundary flux
  M::Array{Float64, 1}       # mass matrix
  Minv::Array{Float64, 1}    # inverse mass matrix
  Minv3D::Array{Float64, 3}    # inverse mass matrix for application to res, not res_vec
  multiplyA0inv::Function       # multiply an array by inv(A0), where A0
                                # is the coefficient matrix of the time
                                # derivative
  src_func::SRCType  # functor for source term
  flux_func::FluxType  # functor for the face flux
  majorIterationCallback::Function # called before every major (Newton/RK) itr

  function AdvectionData_{Tsol, Tres, Tdim, Tmsh}(mesh::AbstractMesh, sbp::AbstractOperator, opts) where {Tsol, Tres, Tdim, Tmsh}
    println("\nConstruction AdvectionData object")
    println("  Tsol = ", Tsol)
    println("  Tres = ", Tres)
    println("  Tdim = ", Tdim)
    println("  Tmsh = ", Tmsh)
    numfacenodes = mesh.numNodesPerFace

    eqn = new()  # incomplete initialization
    eqn.comm = mesh.comm
    eqn.commsize = mesh.commsize
    eqn.myrank = mesh.myrank

    eqn.params = ParamType{Tsol, Tres, Tdim}(mesh, sbp, opts)
    eqn.t = 0.0
    eqn.res_type = Tres
    eqn.majorIterationCallback = majorIterationCallback
    eqn.M = calcMassMatrix(mesh, sbp, eqn)
    eqn.Minv = calcMassMatrixInverse(mesh, sbp, eqn)
    eqn.Minv3D = calcMassMatrixInverse3D(mesh, sbp, eqn)
    eqn.q = zeros(Tsol, 1, sbp.numnodes, mesh.numEl)
    eqn.q_bar = zeros(Tres, 1, sbp.numnodes, mesh.numEl)
    eqn.aux_vars = zeros(Tsol, 0, 0, 0)

    if opts["precompute_volume_flux"]
      eqn.flux_parametric = zeros(Tsol, 1, mesh.numNodesPerElement, mesh.numEl,
                                  Tdim)
    else
      eqn.flux_parametric = zeros(Tsol, 0, 0, 0, 0)
    end

    eqn.res = zeros(Tres, 1, sbp.numnodes, mesh.numEl)
    eqn.res_edge = zeros(Tres, 0, 0, 0, 0)
    eqn.res_bar = zeros(Tres, 1, mesh.numNodesPerElement, mesh.numEl)
    if mesh.isDG
      eqn.q_vec = reshape(eqn.q, mesh.numDof)
      eqn.res_vec = reshape(eqn.res, mesh.numDof)
    else
      eqn.q_vec = zeros(Tres, mesh.numDof)
      eqn.res_vec = zeros(Tres, mesh.numDof)
    end

    if opts["precompute_boundary_flux"]
      eqn.bndryflux = zeros(Tsol, 1, numfacenodes, mesh.numBoundaryFaces)
    else
      eqn.bndryflux = zeros(Tsol, 0, 0, 0)
    end

    eqn.multiplyA0inv = matVecA0inv

    if opts["precompute_q_face"]
      eqn.q_face = zeros(Tsol, 1, 2, numfacenodes, mesh.numInterfaces)
    else
      eqn.q_face = zeros(Tsol, 0, 0, 0, 0)
    end

    if opts["precompute_q_bndry"]
      eqn.q_bndry = zeros(Tsol, 1, numfacenodes, mesh.numBoundaryFaces)
    else
      eqn.q_bndry = zeros(Tsol, 0, 0, 0)
    end

    if opts["precompute_face_flux"]
      eqn.flux_face = zeros(Tres, 1, numfacenodes, mesh.numInterfaces)

      if mesh.isDG
        eqn.flux_sharedface = Array{Array{Tres, 3}}(mesh.npeers)
        for i=1:mesh.npeers
          eqn.flux_sharedface[i] = zeros(Tres, 1, numfacenodes,
                                         mesh.peer_face_counts[i])
        end
      else
        eqn.flux_sharedface = Array{Array{Tres, 3}}(0)
      end  # end if isDG

    else
      eqn.flux_face = zeros(Tres, 0, 0, 0)
      eqn.flux_sharedface = Array{Array{Tres, 3}}(0)
    end  # end if precompute_face_flux

    if mesh.isDG
      eqn.shared_data = getSharedFaceData(Tsol, mesh, sbp, opts,
                                          opts["parallel_data"])
      eqn.shared_data_bar = getSharedFaceData(Tres, mesh, sbp, opts, PARALLEL_DATA_ELEMENT)

    else
      eqn.shared_data = Array{SharedFaceData{Tsol}}(0)
      eqn.shared_data_bar = Array{SharedFaceData{Tres}}(0)
    end

    return eqn
  end # ends the constructor AdvectionData_

end # End type AdvectionData_

import ODLCommonTools.getAllTypeParams

@doc """
### AdvectionEquationMod.getAllTypeParameters

Gets the type parameters for mesh and equation objects.

**Input**

* `mesh` : Object of abstract meshing type.
* `eqn`  : Euler Equation object.
* `opts` : Options dictionary

**Output**

* `tuple` : Tuple of type parameters. Ordering is same as that of the concrete eqn object within this physics module.

"""->
function getAllTypeParams(mesh::AbstractMesh{Tmsh}, eqn::AdvectionData_{Tsol, Tres, Tdim, Tmsh}, opts) where {Tmsh, Tsol, Tres, Tdim}

  tuple = (Tsol, Tres, Tdim, Tmsh)

  return tuple
end

import PDESolver: copyParameters

function copyParameters(eqn_old::AdvectionData, eqn_new::AdvectionData)

  eqn_new.params.LFalpha = eqn_old.params.LFalpha
  eqn_new.alpha_x = eqn_old.alpha_x
  eqn_new.params.alpha_y = eqn_old.params.alpha_y
  eqn_new.params.alpha_z = eqn_old.params.alpha_z
  eqn_new.params.sin_amplitude = eqn_old.params.sin_amplitude
  eqn_new.params.omega = eqn_old.params.omega

  resize!(eqn_new.params.x_design, length(eqn_old.params.x_design))
  copy!(eqn_new.params.x_design, eqn_old.params.x_design)

  return nothing
end


text:
#=
Protected Accessor Function for User.
=#
function Adj(A::Assoc)
    return copy(A.A)
end

function Col(A::Assoc)
    if isa(A.col[1], AbstractString)
        return copy(A.col)
    end
    return 1:size(A.A,2)
end

function Row(A::Assoc)
    if isa(A.row[1], AbstractString)
        return copy(A.row)
    end
    return 1:size(A.A,1)
end

function Val(A::Assoc)
    if isa(A.col[1], AbstractString)
        return copy(A.val)
    end
    return sort(unique(nonzeros(A.A)))
end


text:
using Documenter, SlimPlotting

makedocs(sitename="Slim Plotting toolbox",
         doctest=false, clean=true,
         authors="Mathias Louboutin",
         pages = Any[
             "Home" => "index.md",
             "About" => "README.md",
             "API reference" => "API.md",
         ])

deploydocs(repo="github.com/slimgroup/SlimPlotting.jl")


text:
using StaticArrays, Test
using OrdinaryDiffEq
using RecursiveArrayTools

u0 = [fill(2, MVector{2,Float64}), ones(MVector{2,Float64})]
g(u,p,t) = SA[u[1] + u[2],u[1]]
f = (du, u, p, t) -> begin
    for i in 1:2
        du[i] = g(u[i],p,t)
    end
end
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, Tsit5())
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, Vern9())
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)

u0 = [fill(2, SVector{2,Float64}), ones(SVector{2,Float64})]
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, Tsit5())
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, SSPRK22(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)

u0 = ones(MVector{2,Float64})
ode = ODEProblem(g, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, Tsit5(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)

u0 = ones(SVector{2,Float64})
f = (u,p,t) -> u
ode = ODEProblem(f, u0, (0.,1.))
sol = solve(ode, Euler(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
integrator = init(ode, ImplicitEuler())
@test OrdinaryDiffEq.get_W(integrator.cache.nlsolver) isa StaticArrays.LU
sol = solve(ode, ImplicitEuler())
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, ImplicitEuler(nlsolve=NLAnderson()))
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)
sol = solve(ode, Tsit5(), dt=1e-2)
@test !any(iszero.(sol(1.0))) && !any(sol(1.0) .== u0)

#https://github.com/JuliaDiffEq/DifferentialEquations.jl/issues/373
function lorenz_static(u,p,t)
 dx = 10.0*(u[2]-u[1])
 dy = u[1]*(28.0-u[3]) - u[2]
 dz = u[1]*u[2] - (8/3)*u[3]
 @SVector [dx,dy,dz]
end

u0 = @SVector [1.0,0.0,0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz_static,u0,tspan)
solve(prob,dt=0.1,Rosenbrock23(autodiff=false))

# Check that ArrayPartitions of static vectors work
#https://github.com/SciML/OrdinaryDiffEq.jl/issues/1308
function lorenz_static(u::ArrayPartition,p,t)
 dx = 10.0*(u[2]-u[1])
 dy = u[1]*(28.0-u[3]) - u[2]
 dz = u[1]*u[2] - (8/3)*u[3]
 du1 = @SVector [dx,dy]
 du2 = @SVector [dz]
 ArrayPartition(du1,du2)
end

u01 = @SVector [1.0,0.0]
u02 = @SVector [0.0]
u0ap = ArrayPartition(u01,u02)
probap = ODEProblem(lorenz_static,u0ap,tspan)

sol = solve(prob,dt=1e-2,Heun())
solap = solve(probap,dt=1e-2,Heun())
@test sol(30) ≈ solap(30) atol=1e-12

sol = solve(prob,dt=1e-2,Tsit5())
solap = solve(probap,dt=1e-2,Tsit5())
@test sol(30) ≈ solap(30) atol=4e-6

function rober(u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  dy₁ = -k₁*y₁+k₃*y₂*y₃
  dy₂ =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  dy₃ =  k₂*y₂^2
  SA[dy₁,dy₂,dy₃]
end
prob = ODEProblem{false}(rober,SA[1.0,0.0,0.0],(0.0,1e5),SA[0.04,3e7,1e4])
# Defaults to reltol=1e-3, abstol=1e-6
@test_nowarn sol = solve(prob,Rosenbrock23(chunk_size = Val{3}()),save_everystep=false)
@test_nowarn sol = solve(prob,Rodas4(chunk_size = Val{3}()),save_everystep=false)

function hires_4(u,p,t)
  y1,y2,y3,y4 = u
  dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
  dy2 = 1.71*y1 - 8.75*y2
  dy3 = -10.03*y3 + 0.43*y4 + 0.035*y2
  dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
  SA[dy1,dy2,dy3,dy4]
end

u0 = SA[1,0,0,0.0057]
prob = ODEProblem(hires_4,u0,(0.0,321.8122))
# Defaults to reltol=1e-3, abstol=1e-6
@test_nowarn sol = solve(prob,Rosenbrock23(chunk_size = Val{4}()),save_everystep=false)
@test_nowarn sol = solve(prob,Rodas5(chunk_size = Val{4}()),save_everystep=false)

function hires_5(u,p,t)
  y1,y2,y3,y4,y5 = u
  dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
  dy2 = 1.71*y1 - 8.75*y2
  dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5
  dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
  dy5 = -1.745*y5 + 0.43*y2 + 0.43*y4
  SA[dy1,dy2,dy3,dy4,dy5]
end

u0 = SA[1,0,0,0,0.0057]
prob = ODEProblem(hires_5,u0,(0.0,321.8122))
# Defaults to reltol=1e-3, abstol=1e-6
@test_nowarn sol = solve(prob,Rosenbrock23(chunk_size = Val{5}()),save_everystep=false)
@test_nowarn sol = solve(prob,Rodas4(chunk_size = Val{5}()),save_everystep=false)

function hires(u,p,t)
  y1,y2,y3,y4,y5,y6,y7,y8 = u
  dy1 = -1.71*y1 + 0.43*y2 + 8.32*y3 + 0.0007
  dy2 = 1.71*y1 - 8.75*y2
  dy3 = -10.03*y3 + 0.43*y4 + 0.035*y5
  dy4 = 8.32*y2 + 1.71*y3 - 1.12*y4
  dy5 = -1.745*y5 + 0.43*y6 + 0.43*y7
  dy6 = -280.0*y6*y8 + 0.69*y4 + 1.71*y5 -
           0.43*y6 + 0.69*y7
  dy7 = 280.0*y6*y8 - 1.81*y7
  dy8 = -280.0*y6*y8 + 1.81*y7
  SA[dy1,dy2,dy3,dy4,dy5,dy6,dy7,dy8]
end

u0 = SA[1,0,0,0,0,0,0,0.0057]
prob = ODEProblem(hires,u0,(0.0,321.8122))
# Defaults to reltol=1e-3, abstol=1e-6
@test_nowarn sol = solve(prob,Rosenbrock23(chunk_size = Val{8}()),save_everystep=false)
@test_nowarn sol = solve(prob,Rodas5(chunk_size = Val{8}()),save_everystep=false)

const k1=.35e0
const k2=.266e2
const k3=.123e5
const k4=.86e-3
const k5=.82e-3
const k6=.15e5
const k7=.13e-3
const k8=.24e5
const k9=.165e5
const k10=.9e4
const k11=.22e-1
const k12=.12e5
const k13=.188e1
const k14=.163e5
const k15=.48e7
const k16=.35e-3
const k17=.175e-1
const k18=.1e9
const k19=.444e12
const k20=.124e4
const k21=.21e1
const k22=.578e1
const k23=.474e-1
const k24=.178e4
const k25=.312e1

function pollu(y,p,t)
 r1  = k1 *y[1]
 r2  = k2 *y[2]*y[4]
 r3  = k3 *y[5]*y[2]
 r4  = k4 *y[7]
 r5  = k5 *y[7]
 r6  = k6 *y[7]*y[6]
 r7  = k7 *y[9]
 r8  = k8 *y[9]*y[6]
 r9  = k9 *y[11]*y[2]
 r10 = k10*y[11]*y[1]
 r11 = k11*y[13]
 r12 = k12*y[10]*y[2]
 r13 = k13*y[14]
 r14 = k14*y[1]*y[6]
 r15 = k15*y[3]
 r16 = k16*y[4]
 r17 = k17*y[4]
 r18 = k18*y[16]
 r19 = k19*y[16]
 r20 = k20*y[17]*y[6]
 r21 = k21*y[19]
 r22 = k22*y[19]
 r23 = k23*y[1]*y[4]
 r24 = k24*y[19]*y[1]
 r25 = k25*y[20]

 dy1  = -r1-r10-r14-r23-r24+
          r2+r3+r9+r11+r12+r22+r25
 dy2  = -r2-r3-r9-r12+r1+r21
 dy3  = -r15+r1+r17+r19+r22
 dy4  = -r2-r16-r17-r23+r15
 dy5  = -r3+r4+r4+r6+r7+r13+r20
 dy6  = -r6-r8-r14-r20+r3+r18+r18
 dy7  = -r4-r5-r6+r13
 dy8  = r4+r5+r6+r7
 dy9  = -r7-r8
 dy10 = -r12+r7+r9
 dy11 = -r9-r10+r8+r11
 dy12 = r9
 dy13 = -r11+r10
 dy14 = -r13+r12
 dy15 = r14
 dy16 = -r18-r19+r16
 dy17 = -r20
 dy18 = r20
 dy19 = -r21-r22-r24+r23+r25
 dy20 = -r25+r24
 SA[dy1,dy2,dy3,dy4,dy5,dy6,dy7,dy8,dy9,dy10,dy11,dy12,dy13,dy14,dy15,dy16,dy17,dy18,dy19,dy20]
end

u0 = zeros(20)
u0[2]  = 0.2
u0[4]  = 0.04
u0[7]  = 0.1
u0[8]  = 0.3
u0[9]  = 0.01
u0[17] = 0.007
u0 = SA[u0...]
prob = ODEProblem(pollu,u0,(0.0,60.0))
@test_nowarn sol = solve(prob,Rosenbrock23(chunk_size = Val{8}()),save_everystep=false)
@test_nowarn sol = solve(prob,Rodas5(chunk_size = Val{8}()),save_everystep=false)


text:
using Test
using Indexes
using FormatSpecimens

@testset "Indexes" begin

@testset "Tabix" begin

    # BAI
    path = joinpath(path_of_format("BAM"), "GSE25840_GSM424320_GM06985_gencode_spliced.head.bam.bai")
    path = joinpath(path_of_format("BAM"), "R_12h_D06.uniq.q40.bam.bai")
    path = joinpath(path_of_format("BAM"), "cigar-64k.bam.bai")

    # TBI
    path = joinpath(path_of_format("BED"), "ws245Genes.WBGene.bed.bgz.tbi")
    path = joinpath(path_of_format("GFF3"), "TAIR10.part.gff.bgz.tbi")

    @info path

    open(path) do io

        input = Indexes.CodecBGZF.BGZFDecompressorStream(io)
        @info input

        seekstart(input)

        # check magic bytes
        T = read(input, UInt8)
        @info T

        B = read(input, UInt8)
        @info B

        I = read(input, UInt8)
        @info I

        x = read(input, UInt8)
        @info x

    end

    index = Indexes.Tabix(path)

end

@testset "BGZF" begin

# #
# Indexes.reg2bin(-1, 0)
#
# # The BAI index format for BAM files
#
# goodfiles = filter(entry-> hastag(entry, "bai"), list_valid_specimens("BAM"))
#
# entry = goodfiles[1]
#
# # Get the full path of a file in the entry:
# path_bam = joinpath(path_of_format("BAM"), filename(entry))
# path_bai = path_bam * ".bai"
#
# stream = open(path_bai)
#
# # Read magic bytes
# str = read(stream, 4)
#
# # read contents
# n_refs = read(stream, Int32)
#
#
# indexes = Indexes.read_bgzfindex(stream, n_refs)

end

end


text:
# Val{true} is passed at the last mp iteration and Val{false} during the previous iterations.

function approxTranslation( cmatrix::Array{T,N}, alg::String, ::Val{true} ) where {T<:AbstractFloat,N}
    ccenter = div.( size( cmatrix ), 2 );
    peak1, maxVal = firstPeak( cmatrix );
    return subPixel( cmatrix, alg, peak1, maxVal ) .- ccenter;
end

function approxTranslation( cmatrix::Array{T,N}, alg::String, ::Val{false} ) where {T<:AbstractFloat,N}
    ccenter = div.( size( cmatrix ) , 2 );
    return firstPeak( cmatrix )[1] .- ccenter
end

function subPixel( cmatrix::Array{T,N}, alg::String, peak1::NTuple{N,I}, maxVal::T ) where {T<:AbstractFloat,N}

    refinement = zeros( T, N );

    if alg == "" || indexOnBorder( peak1, size( cmatrix ) )
        # no sub-pixel refinement
    else
        neighs = getNeighbours( peak1, cmatrix );

        # ensuring values are positive
        minVal = minimum( neighs );
        for idx in 1:length(neighs)
            neighs[idx] = 1 + neighs[idx] - minVal
        end

        if     alg == "gaussian"
            refinement = point3Gaussian( neighs, 1+maxVal-minVal );
        elseif alg == "centroid"
            refinement = centroid( neighs, maxVal );
        end

        # nan occur when neighbours contains 0's, resulting in 0/0
        refinement = [ isnan(x) ? 0.0 : x for x in refinement ]
    end

    return peak1 .+ refinement
end

# 3-POINT GAUSSIAN SUBPIXEL

# 3-point Gaussian subpixel 2D

function point3Gaussian( ns::Array{T,2}, mx::T ) where {T<:AbstractFloat}
    return point3_2D( log(ns[4]), log(ns[6]), log(ns[2]), log(ns[8]), log(mx) );
end

function point3_2D( up::T, down::T, left::T, right::T, mx::T ) where {T<:AbstractFloat}
    return [ ( up  - down )/( 2* up  - 4*mx + 2*down  ),
             (left - right)/( 2*left - 4*mx + 2*right ) ];
end

# 3-point Gaussian subpixel 3D
function point3Gaussian( ns::Array{T,3}, mx::T ) where {T<:AbstractFloat}
    return point3_3D( log(ns[13]), log(ns[15]), log(ns[11]), log(ns[17]),
                      log(ns[ 5]), log(ns[23]), log(mx) );
end

function point3_3D( up::T, down::T, left::T, right::T, front::T, back::T, mx::T ) where {T<:AbstractFloat}
    return [ (  up  - down )/( 2*  up  - 4*mx + 2*down  ),
             ( left - right)/( 2* left - 4*mx + 2*right ),
             (front - back )/( 2*front - 4*mx + 2*back  ) ];
end

# CENTROID SUBPIXEL

# Centroid subpixel 2D

function centroid( ns::Array{T,2}, mx::T ) where {T<:AbstractFloat}
	return centroid_2D( ns[4], ns[6], ns[2], ns[8], mx )
end

function centroid_2D( up::T, down::T, left::T, right::T, mx::T ) where {T<:AbstractFloat}
    return [ ( up  - down )/(  up  - mx + down  ),
             (left - right)/( left - mx + right ) ];
end

# Centroid subpixel 3D

function centroid( ns::Array{T,3}, mx::T ) where {T<:AbstractFloat}
	return centroid_3D( ns[13], ns[15], ns[11], ns[17], ns[5], ns[23], mx )
end

function centroid_3D( up::T, down::T, left::T, right::T, front::T, back::T, mx::T ) where {T<:AbstractFloat}
    return [ (  up  - down )/(   up  - mx + down  ),
             ( left - right)/(  left - mx + right ),
			 (front - back )/( front - mx + back  ) ];
end