text:
# HMMWV tests sets

@testset "HMMWV" begin

    # Test set #1
    models =[:ThreeDOFv2]
    cases = ["s1","s2","s3","s4","s6"]; t = 1;
    @testset "Test set #1, test #$(t) \n cases with (case=>$(case_name)) \n models with (model=>$(model))  " for case_name in cases, model in models
      c = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/planner/","RTPP",".yaml")))
      c["vehicle"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/vehicle/","hmmwv",".yaml")))
      c["goal"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["goal"]
      c["X0"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["X0"]
      c["obstacle"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["obstacle"]
      setConfig(c,"misc";(:N=>40),(:model=>model),(:solver=>:Ipopt),(:integrationScheme=>:trapezoidal))
      fixYAML(c)   # fix messed up data types
      n = initializeAutonomousControl(c);
      simMPC!(n;updateFunction=updateAutoParams!,checkFunction=checkCrash)
      @show n.r.ocp.dfsOpt[:tSolve][end]
      @test n.f.mpc.goalReached
      t = t + 1
    end

    if false
        # Test set #2
        models =[:KinematicBicycle2]
        cases = ["s10"]; t = 1;
        @testset "Test set #2, test #$(t) \n cases with (case=>$(case_name)) \n models with (model=>$(model))  " for case_name in cases, model in models
          c = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/planner/","RTPP",".yaml")))
          c["vehicle"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/vehicle/","hmmwv",".yaml")))
          c["goal"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["goal"]
          c["X0"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["X0"]
          c["obstacle"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["obstacle"]
          setConfig(c,"misc";(:N=>40),(:model=>model),(:solver=>:Ipopt),(:integrationScheme=>:trapezoidal))
    #      setConfig(c, "misc";(:sm1=>3.5),(:sm2=>6.5))
          fixYAML(c)   # fix messed up data types
          n = initializeAutonomousControl(c);
          simMPC!(n;updateFunction=updateAutoParams!,checkFunction=checkCrash)
          @show n.r.ocp.dfsOpt[:tSolve][end]
          @test n.f.mpc.goalReached
          t = t + 1
        end
        # Test set #3, constant speed model
        cases = ["s10"]; t = 1;
        @testset "Test set #3, test #$(t) \n cases with (case=>$(case_name))" for case_name in cases
          c = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/planner/","RTPP",".yaml")))
          c["vehicle"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/vehicle/","hmmwv",".yaml")))
          c["goal"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["goal"]
          c["X0"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["X0"]
          c["obstacle"] = load(open(string(Pkg.dir("MichiganAutonomousVehicles"),"/config/case/",case_name,".yaml")))["obstacle"]
          setConfig(c,"misc";(:N=>40),(:model=>:ThreeDOFv3),(:solver=>:Ipopt),(:integrationScheme=>:bkwEuler))
          # fails with:  [ux = 5.0, N = 15 => :crash;
          setConfig(c,"X0";(:ux=>7.0))# the type needs to be a float
          setConfig(c, "misc";(:sm1=>3.),(:sm2=>6.0))
          # otherwise, this will pop up: WARNING: Ipopt finished with status Invalid_Number_Detected
          fixYAML(c)   # fix messed up data types
          n = initializeAutonomousControl(c);
          simMPC!(n;updateFunction=updateAutoParams!,checkFunction=checkCrash)
          @show n.r.ocp.dfsOpt[:tSolve][end]
          @test n.f.mpc.goalReached
          t = t + 1
        end
    end
end


text:
module Excavator

include("job.jl")
include("manager.jl")
include("node.jl")

end # module

ext: jl


text:
#=##############################################################################################
Filename: system.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: this function is the top-level of the simulation
=###############################################################################################

"""
Runs a simulation.

Inputs:

* `aircraft::Aircraft` : `Aircraft` object defining the aircraft
* `parameters` : object containing additional parameters required by actions
* `actions::Vector{Action}` : contains functions inheriting from the `Action` type to be performed in chronological order each step
* `freestream_function::Function` : function returns a `::Freestream` object each step
* `environment_function::Function` : function returns a `::Environment` object each step
* `objective_function::Function` : function defines what is returned
* `step_range` : iterable object determining each step
* `step_symbol` : symbol defining the step, e.g. `alpha` or `time`

Outputs:

* `objective` : results of `objective_function(aircraft, parameters, freestream, environment, step_range)` after simulation

"""
function runsimulation!(aircraft, parameters, actions, freestream_function, environment_function, postactions, objective_function, step_range, step_symbol;
        verbose = true
    )
    # begin simulation
    if verbose
        println("################################################################################################")
        println("                                       AIRCRAFT SYSTEMS")
        println("################################################################################################")
        println("\nRunning simulation...\n")
    end
    # step
    for (stepi, step) in enumerate(step_range)
        if verbose; println("\tstep: $stepi\t\t$step_symbol=$step"); end
        environment = environment_function(aircraft, parameters, step_range, stepi)
        freestream = freestream_function(aircraft, parameters, environment, step_range, stepi)
        for action! in actions
            flag = action!(aircraft, parameters, freestream, environment, step_range, stepi, step_symbol)
            if flag
                @warn "action $action failed at step=$stepi, $step_symbol=$step"
            end
        end
    end
    # perform post-processing
    if verbose; println("\nPerforming postprocessing...\n"); end
    for postaction! in postactions
        flag = postaction!(aircraft, parameters, step_range, step_symbol)
        if flag
            @warn "postaction $postaction failed"
        end
    end
    if verbose; println("\nSimulation complete.\n"); end
    # return objective
    return objective_function(aircraft, parameters, freestream, environment, step_range)
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
using CommonRLInterface

const RL = CommonRLInterface

mutable struct RockPaperScissors <: AbstractEnv
    status::Symbol # either :start, the play of the first player, or :done
end

RockPaperScissors() = RockPaperScissors(:start)

beats(a, b) = (a==:rock && b==:scissors) || (a==:scissors && b==:paper) || (a==:paper && b==:rock)

# Really all_act!, actions, terminated, players, and reset! are all that's needed to describe the game

@provide function RL.all_act!(env::RockPaperScissors, as)
    env.status = :done
    if beats(as[1], as[2])
        return (1, -1)
    elseif beats(as[2], as[1])
        return (-1, 1)
    else
        return (0, 0)
    end
end

RL.actions(env::RockPaperScissors, player=0) = (:rock, :paper, :scissors)
RL.terminated(env::RockPaperScissors) = env.status == :done
RL.reset!(env::RockPaperScissors) = env.status = :start
@provide RL.players(env::RockPaperScissors) = 1:2

# We may also wish to implement the rest of the required interface

RL.observe(env::RockPaperScissors) = [0]

function RL.act!(env::RockPaperScissors, a)
    if env.status == :start
        env.status = a
        return (0, 0)
    else
        return all_act!(env, (env.status, a))
    end
end

@provide RL.player(env::RockPaperScissors) = env.status == :start ? 1 : 2
@provide RL.UtilityStyle(env::RockPaperScissors) = ZeroSum()


text:
using Test
using MPI

MPI.Init()
comm = MPI.COMM_WORLD
nnodes = MPI.Comm_size(comm)
ndims = 3
reorder = 1
periods = [0,1,0]
dims = [0,0,0]
MPI.Dims_create!(nnodes, dims)
comm_cart = MPI.Cart_create(comm, dims, periods, reorder)

rank = MPI.Comm_rank(comm)
ccoords = Cint[-1,-1,-1]
MPI.Cart_coords!(comm_cart, rank, ndims, ccoords)
@test all(ccoords .>= 0)
@test all(ccoords .< dims)

coords = MPI.Cart_coords!(comm_cart, ndims)
@test all(coords .>= 0)
@test all(coords .< dims)

MPI.Finalize()


text:
q = Biquaternion(Quaternion(0, ℝ³(0, 0, 1)), ℝ³(0, 0, 0))
triad = constructtriad(q, length = 1.0)
torus = constructtorus(q, 1, 10)
sphere = constructsphere(q, 1)
cylinder = constructcylinder(q, 1.0, 0.1)
scale = ℝ³(1, 1, 1)
box = constructbox(ℝ³(rand(3)), Quaternion(rand() * 2pi, ℝ³(rand(3))), scale)
plane = getplane([0.0; 5.0; 3.5], Quaternion(1, 0, 0, 0), [1.0; 13.0; 13.0])

@test typeof(triad) == Array{ℝ³,1}
@test size(triad) == (6,)
@test typeof(torus) == Array{ℝ³,2}
@test size(torus) == (36, 36)
@test typeof(sphere) == Array{ℝ³,2}
@test size(sphere) == (36, 36)
@test typeof(cylinder) == Array{ℝ³,2}
@test size(cylinder) == (36, 36)
@test typeof(box) == Array{ℝ³,2}
@test size(box) == (4, 4)
@test typeof(plane) == Array{Float64,3}
@test size(plane) == (2, 2, 3)

## whirl

# Map from S² into its upper hemisphere
s2tos2map(b::S²) = begin
    p = Geographic(b)
    r = sqrt((1 - sin(p.θ)) / 2)
    Geographic(p.r, r * cos(p.ϕ), r * sin(p.ϕ))
end

number = rand(5:10)
points = [Geographic(rand(), rand() * 2pi - pi, rand() * pi - pi / 2) for i in 1:number]
top = U1(rand() * 2pi - pi)
bottom = U1(rand() * 2pi - pi)
s3rotation = Quaternion(rand() * 2pi - pi, ℝ³(rand(3)))
config = Biquaternion(Quaternion(rand() * 2pi - pi, ℝ³(rand(3))), ℝ³(rand(3)))
segments = rand(5:10)
s2tos3map = rand() > 0.5 ? σmap : τmap
whirl = constructwhirl(points,
                       s2tos3map,
                       s2tos2map,
                       top = top,
                       bottom = bottom,
                       s3rotation = s3rotation,
                       config = config,
                       segments = segments)

@test typeof(whirl) == Array{ℝ³,2}
@test size(whirl) == (segments, number)

## frame

circle = U1(rand() * 2pi - pi)
s3rotation = Quaternion(rand() * 2pi - pi, ℝ³(rand(3)))
config = Biquaternion(Quaternion(rand() * 2pi - pi, ℝ³(rand(3))), ℝ³(rand(3)))
segments = rand(5:10)
s2tos3map = rand() > 0.5 ? σmap : τmap
frame = constructframe(circle,
                       s2tos3map,
                       s2tos2map,
                       s3rotation = s3rotation,
                       config = config,
                       segments = segments)

@test typeof(frame) == Array{ℝ³,2}
@test size(frame) == (segments, segments)

max_stars_repo_name: iamazadi/geometrizer

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
struct Rows{Row,Dimensions,Columns,Names} <: AbstractArray{Row,Dimensions}
    columns::Columns
    names::Names
    @propagate_inbounds function Rows{Row,Dimension,Columns,Names}(
        columns,
        the_names,
    ) where {Row,Dimension,Columns,Names}
        @boundscheck if !same_axes(columns...)
            throw(DimensionMismatch("All columns passed to `Rows` must have the same axes"))
        end
        new{Row,Dimension,Columns,Names}(columns, the_names)
    end
end
export Rows

function same_axes(first_column, rest...)
    my_all(
        partial_map((function (reference_axes, item)
            axes(item) == reference_axes
        end), axes(first_column), rest)
    )
end
function same_axes()
    true
end
@propagate_inbounds function Rows{Row,Dimension}(
    columns::Columns,
    the_names::Names,
) where {Row,Dimension,Columns,Names}
    Rows{Row,Dimension,Columns,Names}(columns, the_names)
end
@propagate_inbounds function Rows(columns, the_names)
    Rows{
        NamedTuple{my_map(unname, the_names),Tuple{my_map(eltype, columns)...}},
        ndims(get_model(columns)),
    }(
        columns,
        the_names,
    )
end

@propagate_inbounds function Rows(named_columns)
    Rows(my_map(value, named_columns), my_map(key, named_columns))
end

"""
    Rows(named_columns)

Iterator over `rows` of a table. Always lazy. Use [`Peek`](@ref) to view.

```jldoctest Rows
julia> using LightQuery

julia> using Test: @inferred

julia> lazy = @inferred Rows(a = [1, 2], b = [1.0, 2.0])
2-element Rows{NamedTuple{(:a, :b),Tuple{Int64,Float64}},1,Tuple{Array{Int64,1},Array{Float64,1}},Tuple{Name{:a},Name{:b}}}:
 (a = 1, b = 1.0)
 (a = 2, b = 2.0)

julia> @inferred collect(lazy)
2-element Array{NamedTuple{(:a, :b),Tuple{Int64,Float64}},1}:
 (a = 1, b = 1.0)
 (a = 2, b = 2.0)

julia> @inferred Rows(a = [1, 2])
2-element Rows{NamedTuple{(:a,),Tuple{Int64}},1,Tuple{Array{Int64,1}},Tuple{Name{:a}}}:
 (a = 1,)
 (a = 2,)
```

All arguments to Rows must have the same axes. Use `@inbounds` to override the
check.

```jldoctest Rows
julia> result = Rows(a = 1:2, b = 1:3)
ERROR: DimensionMismatch("All columns passed to `Rows` must have the same axes")
[...]
```
"""
@propagate_inbounds function Rows(; columns...)
    Rows(named_tuple(columns.data))
end

function get_model(columns)
    first(columns)
end
function get_model(::Tuple{})
    1:0
end

function parent(rows::Rows)
    get_model(rows.columns)
end

function to_columns(rows::Rows)
    my_map(tuple, rows.names, rows.columns)
end

function axes(rows::Rows, dimensions...)
    axes(parent(rows), dimensions...)
end
function size(rows::Rows, dimensions...)
    size(parent(rows), dimensions...)
end

@propagate_inbounds function getindex(rows::Rows, an_index::Int...)
    NamedTuple(partial_map(
        (@propagate_inbounds function ((columns, an_index), name)
            name, name(columns)[an_index...]
        end),
        (to_columns(rows), an_index),
        rows.names,
    )...)
end

@propagate_inbounds function setindex!(rows::Rows, row::MyNamedTuple, an_index::Int...)
    partial_map(
        (@propagate_inbounds function ((columns, an_index), (name, value))
            name(columns)[an_index...] = value
            nothing
        end),
        (to_columns(rows), an_index),
        row,
    )
    nothing
end
@propagate_inbounds function setindex!(rows::Rows, row, an_index::Int...)
    setindex!(rows, named_tuple(row), an_index...)
end

function push!(rows::Rows, row::MyNamedTuple)
    partial_map((function (columns, (name, value))
        push!(name(columns), value)
        nothing
    end), to_columns(rows), row)
    nothing
end
function push!(rows::Rows, row)
    push!(rows, named_tuple(row))
end

function val_fieldtypes(something)
    ()
end
@pure function val_fieldtypes(a_type::DataType)
    if a_type.abstract || (a_type.name == Tuple.name && isvatuple(a_type))
        ()
    else
        my_map(Val, (a_type.types...,))
    end
end

function decompose_row_type(::Type{<:NamedTuple{Names,ATuple}}) where {Names,ATuple}
    my_map(tuple, to_Names(Names), val_fieldtypes(ATuple))
end
function decompose_row_type(uninferred)
    ()
end

function similar(rows::Rows, ::Type{ARow}, dimensions::Dims) where {ARow}
    @inbounds Rows(partial_map(
        (function (
            (model, dimensions),
            (name, val_Value)::Tuple{<:Any,Val{Value}},
        ) where {Value}
            name, similar(model, Value, dimensions)
        end),
        (parent(rows), dimensions),
        decompose_row_type(ARow),
    ))
end

function empty(column::Rows{OldRow}, ::Type{NewRow} = OldRow) where {OldRow,NewRow}
    similar(column, NewRow)
end

function widen_column(
    ::HasLength,
    new_length,
    an_index,
    name,
    column::AbstractArray{Element},
    item::Item,
) where {Element,Item<:Element}
    @inbounds column[an_index] = item
    name, column
end
function widen_column(::HasLength, new_length, an_index, name, column::AbstractArray, item)
    name, setindex_widen_up_to(column, item, an_index)
end
function widen_column(
    ::SizeUnknown,
    new_length,
    an_index,
    name,
    column::AbstractArray{Element},
    item::Item,
) where {Element,Item<:Element}
    push!(column, item)
    name, column
end
function widen_column(
    ::SizeUnknown,
    new_length,
    an_index,
    name,
    column::AbstractArray,
    item,
)
    name, push_widen(column, item)
end
function widen_column(
    iterator_size,
    new_length,
    an_index,
    name,
    ::Missing,
    item::Item,
) where {Item}
    new_column = Array{Union{Missing,Item}}(missing, new_length)
    @inbounds new_column[an_index] = item
    name, new_column
end
function widen_column(iterator_size, new_length, an_index, name, ::Missing, ::Missing)
    name, Array{Missing}(missing, new_length)
end

function get_new_length(::SizeUnknown, rows, an_index)
    an_index
end
function get_new_length(::HasLength, rows, an_index)
    length(LinearIndices(rows))
end

function widen_named(iterator_size, rows, row, an_index = length(rows) + 1)
    named_columns = to_columns(rows)
    column_names = my_map(key, named_columns)
    value_names = my_map(key, row)
    just_column_names = my_setdiff(column_names, value_names)
    in_both_names = my_setdiff(column_names, just_column_names)
    @inbounds Rows(partial_map(
        (function (fixeds, variables)
            widen_column(fixeds..., variables...)
        end),
        (iterator_size, get_new_length(iterator_size, rows, an_index), an_index),
        (
            my_map((function ((name, column),)
                name, column, missing
            end), just_column_names(named_columns))...,
            my_map(
                (function ((column_name, column), (value_name, value))
                    column_name, column, value
                end),
                in_both_names(named_columns),
                in_both_names(row),
            )...,
            my_map((function ((name, value),)
                name, missing, value
            end), my_setdiff(value_names, column_names)(row))...,
        ),
    ))
end

function push_widen(rows::Rows, row)
    widen_named(SizeUnknown(), rows, named_tuple(row))
end
function setindex_widen_up_to(rows::Rows, row, an_index)
    widen_named(HasLength(), rows, named_tuple(row), an_index)
end

"""
    make_columns(rows)

Collect into columns.

```jldoctest make_columns
julia> using LightQuery

julia> import Compat: Iterators

julia> using Test: @inferred

julia> stable(x) = (a = x, b = x + 0.0, c = x, d = x + 0.0, e = x, f = x + 0.0);

julia> @inferred make_columns(Iterators.map(stable, 1:4))
(a = [1, 2, 3, 4], b = [1.0, 2.0, 3.0, 4.0], c = [1, 2, 3, 4], d = [1.0, 2.0, 3.0, 4.0], e = [1, 2, 3, 4], f = [1.0, 2.0, 3.0, 4.0])

julia> unstable(x) =
           if x <= 2
               (a = missing, b = string(x), d = string(x))
           else
               (a = x, b = missing, c = x)
           end;

julia> make_columns(Iterators.map(unstable, 1:4))
(d = Union{Missing, String}["1", "2", missing, missing], a = Union{Missing, Int64}[missing, missing, 3, 4], b = Union{Missing, String}["1", "2", missing, missing], c = Union{Missing, Int64}[missing, missing, 3, 4])

julia> make_columns(Iterators.filter(row -> true, Iterators.map(unstable, 1:4)))
(d = Union{Missing, String}["1", "2", missing, missing], a = Union{Missing, Int64}[missing, missing, 3, 4], b = Union{Missing, String}["1", "2", missing, missing], c = Union{Missing, Int64}[missing, missing, 3, 4])
```
"""
function make_columns(rows)
    NamedTuple(to_columns(_collect(
        Rows(),
        rows,
        EltypeUnknown(),
        IteratorSize(rows),
    ))...)
end

export make_columns


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Functions to parse the ICGEM files.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Barthelmes, F., and Förste, C (2011). The ICGEM-format. GFZ Postdam,
#       Departmant 1 "Geodesy and Remote Sensing".
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export parse_icgem

################################################################################
#                                Private Macros
################################################################################

"""
    macro _keyword_found(keyword, keywords_found, current_line)

Check if the `keyword` exists in the list `keywords_found`. If `true`, then
throw an error indicating that the problem occurred on the `current_line`.

"""
macro _keyword_found(keyword, keywords_found, current_line)
    quote
        local lk  = $(esc(keyword))
        local lkf = $(esc(keywords_found))
        local lcl = $(esc(current_line))

        (findfirst(x -> x == lk, lkf) != nothing) &&
        error("Invalid ICGEM format (line $lcl): The keyword \"$lk\" already exists!")
    end
end

"""
    macro _parse_float(input)

Parse the `input` to float substituting all `D`s and `d`s  to `e`, so that we
can convert numbers in FORTRAN format.

"""
macro _parse_float(input)
    quote
        local li = $(esc(input))

        data = replace(li, r"[D,d]" => "e")
        parse(Float64, data)
    end
end

################################################################################
#                               Public Functions
################################################################################

"""
    function parse_icgem(filename::AbstractString)

Parse the ICGEM file `filename` and return an instance of the structure `ICGEM`
with the parsed data.

"""
function parse_icgem(filename::AbstractString)
    # Open the file and find the header.
    file = open(filename, "r")

    # Parse the header
    # ==========================================================================

    header_line_start   = 1
    header_line_end     = 0
    current_line        = 0
    begin_of_head_found = false
    end_of_head_found   = false

    # We need to first check the position of the header due to the
    # `begin_of_head` keyword that can define where the header starts.
    while !eof(file)
        current_line      +=1
        tokens = split(readline(file))

        length(tokens) < 1 && continue

        # Search for the `begin_of_head`, which is optional and makes all
        # previous lines to be ignored.
        if tokens[1] == "begin_of_head"
            # We should not have two `begin_of_head`.
            begin_of_head_found &&
            error("Invalid ICGEM format: Two `begin_of_head` keywords were found!")

            header_line_start = current_line
            begin_of_head_found = true
            continue
        elseif tokens[1] == "end_of_head"
            end_of_head_found = true
            header_line_end   = current_line
            break
        end
    end

    # `end_of_head` keyword is mandatory.
    !end_of_head_found &&
    error("Invalid ICGEM format: The mandatory keyword `end_of_head` was not found!")

    # Rewind file to read again.
    seek(file, 0)

    current_line = 0

    product_type       = ""
    modelname          = ""
    gravity_constant   = 0.0
    radius             = 0.0
    max_degree         = 0
    errors             = :formal
    tide_system        = :zero_tide
    norm               = :fully_normalized
    keywords_found     = Symbol[]
    mandatory_keywords = (:product_type, :modelname, :gravity_constant, :radius,
                          :max_degree, :errors)

    while current_line < header_line_end
        current_line += 1
        tokens = split(readline(file))

        # If we do not have at least 2 tokens, then it is not a keyword.
        length(tokens) < 2 && continue

        if tokens[1] == "product_type"
            @_keyword_found(:product_type, keywords_found, current_line)

            product_type = String(tokens[2])

            product_type != "gravity_field" &&
            error("Only ICGEM files of type \"gravity_field\" are currently supported.")

            push!(keywords_found, :product_type)

        elseif tokens[1] == "modelname"
            @_keyword_found(:modelname, keywords_found, current_line)

            modelname = String(tokens[2])
            push!(keywords_found, :modelname)

        elseif tokens[1] == "earth_gravity_constant" ||
               tokens[1] == "gravity_constant"
            @_keyword_found(:gravity_constant, keywords_found, current_line)

            try
                gravity_constant = @_parse_float(tokens[2])
            catch e
                @error(e.msg)
                error("Invalid ICGEM format (line $current_line): Invalid float number in keyword \"$(tokens[1])\".")
            end

            push!(keywords_found, :gravity_constant)

        elseif tokens[1] == "radius"
            @_keyword_found(:radius, keywords_found, current_line)

            try
                radius = @_parse_float(tokens[2])
            catch e
                @error(e.msg)
                error("Invalid ICGEM format (line $current_line): Invalid float number in keyword \"radius\".")
            end

            push!(keywords_found, :radius)

        elseif tokens[1] == "max_degree"
            @_keyword_found(:max_degree, keywords_found, current_line)

            try
                max_degree = parse(Int, tokens[2])
            catch e
                @error(e.msg)
                error("Invalid ICGEM format (line $current_line): Invalid integer number in keyword \"max_degree\".")
            end

            push!(keywords_found, :max_degree)

        elseif tokens[1] == "errors"
            @_keyword_found(:errors, keywords_found, current_line)

            errors = Symbol(tokens[2])
            push!(keywords_found, :errors)

            if !(errors in [:no, :calibrated, :formal, :calibrated_and_formal])
                @warn("Invalid ICGEM format: The error type \"$(tokens[2])\" is not valid! Assuming \"formal\".")
                errors = :formal
            end

        elseif tokens[1] == "tide_system"
            @_keyword_found(:tide_system, keywords_found, current_line)

            tide_system = Symbol(tokens[2])
            push!(keywords_found, :tide_system)

            if !(tide_system in [:zero_tide, :tide_free, :unknown])
                @warn("Invalid ICGEM format: The tide system \"$(tokens[2])\" is not valid! Assuming \"unknown\".")
                tide_system = :unknown
            end

        elseif tokens[1] == "norm"
            @_keyword_found(:norm, keywords_found, current_line)

            norm = Symbol(tokens[2])
            push!(keywords_found, :norm)

            if !(norm in [:fully_normalized, :unnormalized])
                @warn("Invalid ICGEM format: The norm \"$(tokens[2])\" is not valid! Assuming \"fully_normalized\".")
                norm = :fully_normalized
            end
        end
    end

    # Check if all mandatory keywords were found.
    found_mandatory_keywords::NTuple{6,Bool} =
        map(x->x in keywords_found, mandatory_keywords)
    all_mand_keywords_found  = true

    for i in eachindex(found_mandatory_keywords)
        if found_mandatory_keywords[i] == false
            @error("Invalid ICGEM format: The mandatory keyword \"$(mandatory_keywords[i])\" were not found!")
            all_mand_keywords_found = false
        end
    end

    !all_mand_keywords_found &&
    error("Invalid ICGEM format: Some mandatory keywords were not found!")

    # Parse the coefficients
    # ==========================================================================

    # Read the entire file to improve performance.
    raw = readdlm(filename; skipstart = header_line_end)

    # We currently only support the `gfc` coefficients.
    raw_gfc   = @views raw[ raw[:,1] .== "gfc", :]
    raw_nogfc = @views raw[ raw[:,1] .!= "gfc", :]

    !isempty(raw_nogfc) &&
    @warn("Unsupported ICGEM format: Only the keywords \"gfc\" are currently supported.")

    # Parse: gfc
    # ==========

    # Check if the number of columns is correct.
    if errors == :no
        size(raw_gfc,2) != 5 &&
        error("Invalid ICGME format: The coefficients table must have 5 columns because the keyword errors is \"$errors\".")
    elseif (errors == :calibrated) || (errors == :formal)
        size(raw_gfc,2) != 7 &&
        error("Invalid ICGME format: The coefficients table must have 7 columns because the keyword errors is \"$errors\".")
    else
        size(raw_gfc,2) != 9 &&
        error("Invalid ICGME format: The coefficients table must have 7 columns because the keyword errors is \"$errors\".")
    end

    # Check if there is any number to be converted. In this case, we assume that
    # the number is written in FORTRAN format.
    raw_gfc_coefs = @view(raw_gfc[:,4:end])

    for i in eachindex(raw_gfc_coefs)
        if !(typeof(raw_gfc_coefs[i]) <: AbstractFloat)
            try
                data             = replace(raw_gfc_coefs[i], r"[D,d]" => "e")
                raw_gfc_coefs[i] = parse(Float64, data)
            catch e
                error("Invalid ICGME format: Could not convert the coefficient \"$(raw_gfc_coefs[i])\" to a Float64.")
            end
        end
    end

    # The file will be parsed using Sparse matrices because it is much easier.
    # Perhaps, there is a method with better performance, but this function is
    # called only once per execution.
    @views I = map(x->Int(x)+1, raw_gfc[:,2])
    @views J = map(x->Int(x)+1, raw_gfc[:,3])

    @views Clm = sparse(I, J, convert(Vector{Float64}, raw_gfc[:,4]))
    @views Slm = sparse(I, J, convert(Vector{Float64}, raw_gfc[:,5]))

    if errors == :no
        return ICGEM(product_type     = product_type,
                     modelname        = modelname,
                     gravity_constant = gravity_constant,
                     radius           = radius,
                     max_degree       = max_degree,
                     errors           = errors,
                     tide_system      = tide_system,
                     norm             = norm,
                     Clm              = Clm,
                     Slm              = Slm
                    )

    elseif (errors == :calibrated) || (errors == :formal)
        @views sigmaC = sparse(I, J, convert(Vector{Float64}, raw_gfc[:,6]))
        @views sigmaS = sparse(I, J, convert(Vector{Float64}, raw_gfc[:,7]))

        ICGEM(product_type     = product_type,
              modelname        = modelname,
              gravity_constant = gravity_constant,
              radius           = radius,
              max_degree       = max_degree,
              errors           = errors,
              tide_system      = tide_system,
              norm             = norm,
              Clm              = Clm,
              sigmaC           = sigmaC,
              Slm              = Slm,
              sigmaS           = sigmaS
             )

    else
        @views sigmaC_cal    = sparse(I, J, convert(Vector{Float64}, raw_gfc[:,6]))
        @views sigmaS_cal    = sparse(I, J, convert(Vector{Float64}, raw_gfc[:,7]))
        @views sigmaC_formal = sparse(I, J, convert(Vector{Float64}, raw_gfc[:,8]))
        @views sigmaS_formal = sparse(I, J, convert(Vector{Float64}, raw_gfc[:,9]))

        ICGEM(product_type     = product_type,
              modelname        = modelname,
              gravity_constant = gravity_constant,
              radius           = radius,
              max_degree       = max_degree,
              errors           = errors,
              tide_system      = tide_system,
              norm             = norm,
              Clm              = Clm,
              sigmaC           = sigmaC_cal,
              sigmaC_formal    = sigmaC_formal,
              Slm              = Slm,
              sigmaS           = sigmaS_cal,
              sigmaS_formal    = sigmaS_formal
             )
    end
end


text:
export
    Generate

"""
    Generate{O} <: SFunc{Tuple{Dist{O}}, O}

Generate a value from its `Dist` argument.

This helps in higher-order programming. A typical pattern will be to create an sfunc that produces a `Dist`,
and then generate many observations from the `Dist` using `Generate`.

# Additional supported operators
- `support`
- `support_quality`
- `sample`
- `logcpdf`
- `compute_pi`
- `send_lambda`
- `make_factors`
"""
struct Generate{O} <: SFunc{Tuple{Dist{O}}, O}
end

@impl begin
    struct GenerateSupport end
    function support(::Generate{O},
                    parranges::NTuple{N,Vector},
                    size::Integer,
                    curr::Vector{<:O}) where {J<:Tuple,O,N}

        result = Vector{O}()
        for sf in parranges[1]
            append!(result, support(sf, (), size, curr))
        end
        return unique(result)
    end
end

@impl begin
    struct GenerateSupportQuality end

    function support_quality(::Generate{O}, parranges) where O
        q = support_quality_rank(:CompleteSupport)
        for sf in parranges[1]
            imp = get_imp(MultiInterface.get_policy(), Support, sf, (), 0, O[])
            q = min(q, support_quality_rank(support_quality(imp, sf, ())))
        end
        return support_quality_from_rank(q)
    end
end

@impl begin
    struct GenerateSample end
    function sample(::Generate{O}, input::Tuple{<:Dist{O}})::O where O
        return sample(input[1], ())
    end
end

@impl begin
    struct GenerateLogcpdf end
    function logcpdf(::Generate{O}, i::Tuple{<:Dist{O}}, o::O)::AbstractFloat where O
        return logcpdf(i[1], (), o)
    end
end

# WARNING: THIS LOGIC DOES NOT WORK WITH MORE THAN ONE PARENT
@impl begin
    struct GenerateComputePi end
    function compute_pi(::Generate{O},
                     range::__OptVec{<:O},
                     parranges::NTuple{N,Vector},
                     incoming_pis::Tuple)::Dist{<:O} where {N,O}

        sfrange = parranges[1]
        sfpi = incoming_pis[1]
        result = zeros(Float64, length(range))
        for sf in sfrange
            p1 = cpdf(sfpi, (), sf)
            p2 = compute_pi(sf, range, (), ())
            p3 = [p1 * cpdf(p2, (), x) for x in range]
            result .+= p3
        end
        return Cat(range, result)
    end
end

# WARNING: THIS LOGIC DOES NOT WORK WITH MORE THAN ONE PARENT
@impl begin
    struct GenerateSendLambda end

    function send_lambda(::Generate{O},
                       lambda::Score{<:O},
                       range::__OptVec{<:O},
                       parranges::NTuple{N,Vector},
                       incoming_pis::Tuple,
                       parent_idx::Integer)::Score where {N,O}
        @assert parent_idx == 1
        sfrange::Vector{typeof(parranges[1][1])} = parranges[1]
        sfpi = incoming_pis[1]
        resultprobs = Vector{Float64}()
        for sf in sfrange
            resultpieces = Vector{Float64}()
            ypi = compute_pi(sf, range, (), ())
            for y in range
                push!(resultpieces, logcpdf(ypi, (), y) + get_log_score(lambda, y))
            end
            push!(resultprobs, logsumexp(resultpieces))
        end
        result :: LogScore{typeof(sfrange[1])} = LogScore(sfrange, resultprobs)
        return result
    end
end

@impl begin
    struct GenerateMakeFactors end

    function make_factors(::Generate{O}, range::Vector{<:O}, parranges::Tuple{<:Vector{<:Dist{O}}}, id::Int, parids::Tuple{Int}) where O
        dims = (length(parranges[1]), length(range))
        keys = (parids[1], id)
        entries = Float64[]
        for u in parranges[1]
            for x in range
                p = cpdf(u, (), x)
                push!(entries, p)
            end
        end
        fact = Factor(dims, keys, entries)
        return ([fact], [fact])
    end
end


text:
import Random
import Random: rand, UInt52, rng_native_52

"""
```julia
AbstractRNG{T} <: Random.AbstractRNG
```

The abstract type of Random Number Generators. T indicates the original output type of a RNG.
"""
abstract type AbstractRNG{T<:Number} <: Random.AbstractRNG end

const BitTypes = Union{Bool, UInt8, UInt16, UInt32, UInt64, UInt128, Int8, Int16, Int32, Int64, Int128}

# For compatibility with functions in Random stdlib.
rng_native_52(::AbstractRNG) = UInt64
rand(rng::AbstractRNG, ::Random.SamplerType{T}) where {T} = rand(rng, T)

# see https://github.com/JuliaRandom/RandomNumbers.jl/issues/8
# TODO: find a better approach.
@inline function rand(rng::AbstractRNG{T}, ::Type{Float64}=Float64) where {T<:Union{UInt64, UInt128}}
    reinterpret(Float64, Base.exponent_one(Float64) | rand(rng, UInt52())) - 1.0
end
@inline function rand(rng::AbstractRNG{T}, ::Type{Float64}=Float64) where {T<:Union{UInt8, UInt16, UInt32}}
    rand(rng, T) * exp2(-sizeof(T) << 3)
end

@inline function rand(rng::AbstractRNG{T1}, ::Type{T2}) where {T1<:BitTypes, T2<:BitTypes}
    s1 = sizeof(T1)
    s2 = sizeof(T2)
    t = rand(rng, T1) % T2
    s1 > s2 && return t
    for i in 2:(s2 ÷ s1)
        t |= (rand(rng, T1) % T2) << ((s1 << 3) * (i - 1))
    end
    t
end


text:
using DitherPunk
using ImageCore
using ImageInTerminal
using ReferenceTests
using TestImages

# Load test image
img = testimage("fabio_color_256")

# Run & test fixed pallete dithering methods
algs = Dict("FloydSteinberg" => FloydSteinberg(), "Bayer" => Bayer())

for C in [RGB, HSV]
    for (name, alg) in algs
        img1 = C.(img)
        local img2 = copy(img1)
        local d = dither(img2, alg)
        @test_reference "references/per-channel/$(name)_$(C).txt" d
        @test eltype(d) <: C
        @test img2 == img1 # image not modified

        imshow(d)
        println()

        local d = dither!(img2, alg)
        @test eltype(d) <: C
        @test img2 == d # image updated in-place
    end
end


text:
using BaseTestNext
using GraphLayout
srand(1)

@testset "GraphLayout.jl" begin

@testset "Render a pentagon" begin
    adj_matrix = ones(5,5) - eye(5,5)
    @testset "layout_spring_adj" begin
        loc_x, loc_y = layout_spring_adj(adj_matrix)

        draw_layout_adj(adj_matrix, loc_x, loc_y, filename="pentagon_spring.svg")

        draw_layout_adj(adj_matrix, loc_x, loc_y, filename="pentagon_labeled.svg",
        	labels=collect(1:5), labelsize=2.0)
        draw_layout_adj(adj_matrix, loc_x, loc_y, filename="pentagon_noarrows.svg",
        	arrowlengthfrac=0.0)
        draw_layout_adj(adj_matrix, loc_x, loc_y, filename="pentagon_longarrows.svg",
        	arrowlengthfrac=0.5)
    end  # "layout_spring_adj"

    @testset "layout_stressmajorize_adj" begin
        X = layout_stressmajorize_adj(adj_matrix)
        draw_layout_adj(adj_matrix, X[:,1], X[:,2], filename="pentagon_stress.svg")
    end  # "layout_stressmajorize_adj"
end  # "Render a pentagon"

@testset "Render a random graph" begin

    adj_matrix = full(sprand(100,100,0.02))
    @testset "layout_spring_adj" begin
        loc_x, loc_y = layout_spring_adj(adj_matrix)

        draw_layout_adj(adj_matrix, loc_x, loc_y, filename="random_spring.svg")

        draw_layout_adj(adj_matrix, loc_x, loc_y, filename="random_spring_color.svg",
            labelc="#000000",
            labelsize=4.0,
            nodefillc="#21AAFF",
            nodestrokec="#7BB1B1",
            edgestrokec="#B11B1B")
    end

    @testset "layout_stressmajorize_adj" begin
        X = layout_stressmajorize_adj(adj_matrix)
        draw_layout_adj(adj_matrix, X[:,1], X[:,2], filename="random_stress.svg")
    end
end  # "Render a random graph"

# Trees
@testset "Render a tree" begin
    adj_list = Vector{Int}[
        [2,3,4],
        [4],
        [4],
        []
    ]
    labels = ["Apple", "Banana", "Carrot", "Durian"]

    layout_tree(adj_list,labels,filename="tree_1.svg",
                cycles = false, ordering = :optimal, coord = :optimal,
                xsep = 5, ysep = 10, scale = 0.2, labelpad = 1.2,
                background = nothing)

    layout_tree(adj_list,labels,filename="tree_2.svg",
                cycles = false, ordering = :optimal, coord = :optimal,
                xsep = 1, ysep = 20, scale = 0.1, labelpad = 2,
                background = "#DDDDFF")

    rm("tree_1.svg")
    rm("tree_2.svg")
end

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

@testset "Check that output == cached data" begin

#Compare with cached output
cachedout = joinpath(Pkg.dir("GraphLayout"), "test", "examples")
differentfiles = AbstractString[]
for output in readdir(".")
    endswith(output, ".svg") || continue
    contains(output, "tree") && continue  # don't test trees right now
    cached = open(readall, joinpath(cachedout, output))
    genned = open(readall, joinpath(output))
    if cached != genned
        push!(differentfiles, output)
    else #Delete generated file
        rm(output)
    end
end

#Print out which files differ and their diffs
if length(differentfiles)>0
    #Capture diffs
    diffs = map(
        output -> output * ":\n" *
            readall(ignorestatus(`diff $(joinpath(cachedout, output)) $(joinpath(output))`)) *
            "\n\n",
        differentfiles)
    error(string("Generated output differs from cached test output:\n",
        join(differentfiles, "\n"), "\n\n", join(diffs, "\n")))
else
    println("All files matched!")
end

end  # "Check that output == cached data"

end  # "GraphLayout.jl"


text:
using InfoTheory

NumTestRepetitions = 30

function approxeq(x, y, delta = 1e-3)
  absdelta = abs(x-y)
  if absdelta >= delta
    error("Expected $x to be within $delta of $y (but diff is $absdelta)")
    false
  else
    #println("Expected $x to be within $delta of $y and it was! ($absdelta)")
    true
  end
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""
    Site

Inputs related to the physical location:

```julia
function Site(;
    latitude::Real,
    longitude::Real,
    land_acres::Union{Real, Nothing} = nothing,
    roof_squarefeet::Union{Real, Nothing} = nothing,
    min_resil_time_steps::Int=0,
    mg_tech_sizes_equal_grid_sizes::Bool = true,
    node::Int = 1,
    )
```
"""
struct Site
    "required"
    latitude
    "required"
    longitude
    land_acres
    roof_squarefeet
    min_resil_time_steps
    mg_tech_sizes_equal_grid_sizes
    node  # TODO validate that multinode Sites do not share node numbers? Or just raise warning
    function Site(;
        latitude::Real,
        longitude::Real,
        land_acres::Union{Real, Nothing} = nothing,
        roof_squarefeet::Union{Real, Nothing} = nothing,
        min_resil_time_steps::Int=0,
        mg_tech_sizes_equal_grid_sizes::Bool = true,
        node::Int = 1,
        )
        invalid_args = String[]
        if !(-90 <= latitude < 90)
            push!(invalid_args, "latitude must satisfy -90 <= latitude < 90, got $(latitude)")
        end
        if !(-180 <= longitude < 180)
            push!(invalid_args, "longitude must satisfy -180 <= longitude < 180, got $(longitude)")
        end
        if length(invalid_args) > 0
            error("Invalid argument values: $(invalid_args)")
        end
        new(latitude, longitude, land_acres, roof_squarefeet, min_resil_time_steps,
            mg_tech_sizes_equal_grid_sizes, node)
    end
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
module BenchFilterMapMap!

using BenchmarkTools
using Transducers

function manual_filter_map_map!(ys, xs)
    @inbounds for i in eachindex(xs, ys)
        x = xs[i]
        if -0.5 < x < 0.5
            ys[i] = 2x
        end
    end
    return ys
end

suite = BenchmarkGroup()

xf = Filter(x -> -0.5 < x < 0.5) |> Map(x -> 2x)
n = 10^5

let xs = randn(10),
    ys = zero(xs)
    @assert manual_filter_map_map!(ys, xs) == map!(xf, ys, xs)
end

suite["xf"] = @benchmarkable(
    map!($xf, ys, xs),
    setup=(xs = randn($n); ys = zero(xs)))
suite["man"] = @benchmarkable(
    manual_filter_map_map!(ys, xs),
    setup=(xs = randn($n); ys = zero(xs)))

end  # module
BenchFilterMapMap!.suite


text:
using Test
import UUIDs
import Pluto: PlutoRunner, Notebook, WorkspaceManager, Cell, ServerSession, ClientSession, update_run!
using Pluto.WorkspaceManager: poll

@testset "Logging" begin
    🍭 = ServerSession()
    🍭.options.evaluation.workspace_use_distributed = true

    fakeclient = ClientSession(:fake, nothing)
    🍭.connected_clients[fakeclient.id] = fakeclient

    @testset "Logging respects maxlog" begin
        notebook = Notebook(Cell.([
            """
            for i in 1:10
                @info "logging" i maxlog=2
            end
            """,
        ]))

        update_run!(🍭, notebook, notebook.cells)
        @test notebook.cells[begin] |> noerror

        @test poll(5, 1/60) do
            length(notebook.cells[begin].logs) == 2
        end
        WorkspaceManager.unmake_workspace((🍭, notebook))
    end
end


text:
module DynamicHMCModels

using Reexport, Requires

@reexport using DynamicHMC, LogDensityProblems, TransformVariables
@reexport using Distributions, Random, Statistics
@reexport using Parameters, CSV, DataFrames
@reexport using MonteCarloMeasurements

function __init__()
  @require MCMCChains="c7f686f2-ff18-58e9-bc7b-31028e88f75d" include("require/chains.jl")
end

include("particles.jl")

end # module


text:
real_rep(a::AbstractArray{DualT}) where {TAG, T, DualT<:ForwardDiff.Dual{TAG, T}} = reinterpret(reshape, T, a)
_view1(B::AbstractMatrix) = @view(B[1,:])
_view1(B::AbstractArray{<:Any,3}) = @view(B[1,:,:])

# multiplication of dual vector/matrix by standard matrix from the left
function _matmul!(_C::AbstractVecOrMat{DualT}, A::AbstractMatrix, _B::AbstractVecOrMat{DualT},
                  α, β=Zero(), nthread::Nothing=nothing, MKN=nothing, contig_axis=nothing) where {DualT <: ForwardDiff.Dual}
  B = real_rep(_B)
  C = real_rep(_C)

  @tturbo for n ∈ indices((C, B), 3), m ∈ indices((C, A), (2, 1)), l in indices((C, B), 1)
    Cₗₘₙ = zero(eltype(C))
    for k ∈ indices((A, B), 2)
      Cₗₘₙ += A[m, k] * B[l, k, n]
    end
    C[l, m, n] = α * Cₗₘₙ + β * C[l, m, n]
  end

  _C
end

# multiplication of dual matrix by standard vector/matrix from the right
@inline function _matmul!(_C::AbstractVecOrMat{DualT}, _A::AbstractMatrix{DualT}, B::AbstractVecOrMat,
                          α=One(), β=Zero(), nthread::Nothing=nothing, MKN=nothing) where {TAG, T, DualT <: ForwardDiff.Dual{TAG, T}}
  if Bool(ArrayInterface.is_dense(_C)) && Bool(ArrayInterface.is_column_major(_C)) &&
    Bool(ArrayInterface.is_dense(_A)) && Bool(ArrayInterface.is_column_major(_A))
    # we can avoid the reshape and call the standard method
    A = reinterpret(T, _A)
    C = reinterpret(T, _C)
    _matmul!(C, A, B, α, β, nthread, nothing)
  else
    # we cannot use the standard method directly
    A = real_rep(_A)
    C = real_rep(_C)

    @tturbo for n ∈ indices((C, B), (3, 2)), m ∈ indices((C, A), 2), l in indices((C, A), 1)
      Cₗₘₙ = zero(eltype(C))
      for k ∈ indices((A, B), (3, 1))
        Cₗₘₙ += A[l, m, k] * B[k, n]
      end
      C[l, m, n] = α * Cₗₘₙ + β * C[l, m, n]
    end
  end

  _C
end

@inline function _matmul!(_C::AbstractVecOrMat{DualT}, _A::AbstractMatrix{DualT}, _B::AbstractVecOrMat{DualT},
                          α=One(), β=Zero(), nthread::Nothing=nothing, MKN=nothing, contig=nothing) where {TAG, T, P, DualT <: ForwardDiff.Dual{TAG, T, P}}
  A = real_rep(_A)
  C = real_rep(_C)
  B = real_rep(_B)
  if Bool(ArrayInterface.is_dense(_C)) && Bool(ArrayInterface.is_column_major(_C)) &&
    Bool(ArrayInterface.is_dense(_A)) && Bool(ArrayInterface.is_column_major(_A))
    # we can avoid the reshape and call the standard method
    Ar = reinterpret(T, _A)
    Cr = reinterpret(T, _C)
    _matmul!(Cr, Ar, _view1(B), α, β, nthread, nothing)
  else
    # we cannot use the standard method directly
    @tturbo for n ∈ indices((C, B), 3), m ∈ indices((C, A), 2), l in indices((C, A), 1)
      Cₗₘₙ = zero(eltype(C))
      for k ∈ indices((A, B), (3, 2))
        Cₗₘₙ += A[l, m, k] * B[1, k, n]
      end
      C[l, m, n] = α * Cₗₘₙ + β * C[l, m, n]
    end
  end
  Pstatic = static(P)
  @tturbo for n ∈ indices((B,C),3), m ∈ indices((A,C),2), p ∈ 1:Pstatic
    Cₚₘₙ = zero(eltype(C))
    for k ∈ indices((A,B),(3,2))
      Cₚₘₙ += A[1,m,k] * B[p+1,k,n]
    end
    C[p+1,m,n] = C[p+1,m,n] + α*Cₚₘₙ
  end
  _C
end

# multiplication of dual vector/matrix by standard matrix from the left
function _matmul_serial!(_C::AbstractVecOrMat{DualT}, A::AbstractMatrix, _B::AbstractVecOrMat{DualT},
                         α, β, MKN) where {DualT <: ForwardDiff.Dual}
  B = real_rep(_B)
  C = real_rep(_C)

  @turbo for n ∈ indices((C, B), 3), m ∈ indices((C, A), (2, 1)), l in indices((C, B), 1)
    Cₗₘₙ = zero(eltype(C))
    for k ∈ indices((A, B), 2)
      Cₗₘₙ += A[m, k] * B[l, k, n]
    end
    C[l, m, n] = α * Cₗₘₙ + β * C[l, m, n]
  end

  _C
end

# multiplication of dual matrix by standard vector/matrix from the right
@inline function _matmul_serial!(_C::AbstractVecOrMat{DualT}, _A::AbstractMatrix{DualT}, B::AbstractVecOrMat,
                                 α, β, MKN) where {TAG, T, DualT <: ForwardDiff.Dual{TAG, T}}
  if Bool(ArrayInterface.is_dense(_C)) && Bool(ArrayInterface.is_column_major(_C)) &&
    Bool(ArrayInterface.is_dense(_A)) && Bool(ArrayInterface.is_column_major(_A))
    # we can avoid the reshape and call the standard method
    A = reinterpret(T, _A)
    C = reinterpret(T, _C)
    _matmul_serial!(C, A, B, α, β, nothing)
  else
    # we cannot use the standard method directly
    A = real_rep(_A)
    C = real_rep(_C)

    @turbo for n ∈ indices((C, B), (3, 2)), m ∈ indices((C, A), 2), l in indices((C, A), 1)
      Cₗₘₙ = zero(eltype(C))
      for k ∈ indices((A, B), (3, 1))
        Cₗₘₙ += A[l, m, k] * B[k, n]
      end
      C[l, m, n] = α * Cₗₘₙ + β * C[l, m, n]
    end
  end

  _C
end

@inline function _matmul_serial!(_C::AbstractVecOrMat{DualT}, _A::AbstractMatrix{DualT}, _B::AbstractVecOrMat{DualT},
                                 α, β, MKN) where {TAG, T, P, DualT <: ForwardDiff.Dual{TAG, T, P}}
  A = real_rep(_A)
  C = real_rep(_C)
  B = real_rep(_B)
  if Bool(ArrayInterface.is_dense(_C)) && Bool(ArrayInterface.is_column_major(_C)) &&
    Bool(ArrayInterface.is_dense(_A)) && Bool(ArrayInterface.is_column_major(_A))
    # we can avoid the reshape and call the standard method
    Ar = reinterpret(T, _A)
    Cr = reinterpret(T, _C)
    _matmul_serial!(Cr, Ar, _view1(B), α, β, nothing)
  else
    # we cannot use the standard method directly
    @turbo for n ∈ indices((C, B), 3), m ∈ indices((C, A), 2), l in indices((C, A), 1)
      Cₗₘₙ = zero(eltype(C))
      for k ∈ indices((A, B), (3, 2))
        Cₗₘₙ += A[l, m, k] * B[1, k, n]
      end
      C[l, m, n] = α * Cₗₘₙ + β * C[l, m, n]
    end
  end
  Pstatic = static(P)
  @turbo for n ∈ indices((B,C),3), m ∈ indices((A,C),2), p ∈ 1:Pstatic
    Cₚₘₙ = zero(eltype(C))
    for k ∈ indices((A,B),(3,2))
      Cₚₘₙ += A[1,m,k] * B[p+1,k,n]
    end
    C[p+1,m,n] = C[p+1,m,n] + α*Cₚₘₙ
  end
  _C
end


text:
using Hiccup

c(a...) = Any[a...]

fade(x) = span(".fade", x)

icon(x) = span(".icon.icon-$x", [])


text:
mutable struct SparseSymmetricBlockMatrix{T<:Number}
    blocks::Dict{Any,SparseSymmetricMatrix{T}}
end

SparseSymmetricBlockMatrix(T::Type) = SparseSymmetricBlockMatrix(Dict{Any,SparseSymmetricMatrix{T}}())

SparseSymmetricBlockMatrix() = SparseSymmetricBlockMatrix(Float64)

blocks(m::SparseSymmetricBlockMatrix) = m.blocks

function setindex!(bm::SparseSymmetricBlockMatrix{Float64}, v::Float64, bi, i, j)
    if haskey(blocks(bm), bi)
        blocks(bm)[bi][i, j] = v
    else
        m = SparseSymmetricMatrix(Float64)
        m[i, j] = v
        blocks(bm)[bi] = m
    end
end

function setindex!(bm::SparseSymmetricBlockMatrix{Float64}, m::SparseSymmetricMatrix{Float64}, bi)
    blocks(bm)[bi] = m
end

function getindex(bm::SparseSymmetricBlockMatrix{Float64}, bi, i, j)
    get(blocks(bm), bi, SparseSymmetricMatrix(Float64))[i, j]
end

getindex(bm::SparseSymmetricBlockMatrix{Float64}, bi) = get(blocks(bm), bi, SparseSymmetricMatrix(Float64))

# start(m::SparseSymmetricBlockMatrix) = start(blocks(m))
# next(m::SparseSymmetricBlockMatrix, state) = next(blocks(m), state)
# done(m::SparseSymmetricBlockMatrix, state) = done(blocks(m), state)

copy(m::SparseSymmetricBlockMatrix) = SparseSymmetricBlockMatrix(copy(blocks(m)))

function *(x::Float64, m::SparseSymmetricBlockMatrix{Float64})
    p = SparseSymmetricBlockMatrix(Float64)
    for k in keys(blocks(m))
        p[k] = x * m[k]
    end
    p
end

function show(io::IO, m::SparseSymmetricBlockMatrix{Float64})
    println(io, typeof(m))
    for (blockindex, matrix) in blocks(m)
        for ((i, j), v) in matrix.entries
            println(io, "$blockindex\t$i\t$j\t$v")
        end
    end
end

==(A::SparseSymmetricBlockMatrix, B::SparseSymmetricBlockMatrix) = blocks(A) == blocks(B)

function isapprox(A::SparseSymmetricBlockMatrix, B::SparseSymmetricBlockMatrix)
    nb_blocks = length(blocks(A))
    if nb_blocks != length(blocks(B))
        return false
    end

    for b in 1:nb_blocks
        for ((i,j),val) in A.blocks[b].entries
            if !isapprox(val, B.blocks[b].entries[i,j])
                return false
            end
        end
    end
    true
end

max_stars_repo_name: JAndriesJ/SemidefiniteProgramming

max_issues_repo_name: JAndriesJ/SemidefiniteProgramming

max_forks_repo_name: JAndriesJ/SemidefiniteProgramming


text:
type PSRacos
    rc::RacosCommon
    evaluation_server_num
    sample_set
    result_set
    asyn_result
    history
    is_finish

    function PSRacos(ncomputer)
        new(RacosCommon(), ncomputer, RemoteChannel(()->Channel(ncomputer)),
        RemoteChannel(()->Channel(ncomputer)), RemoteChannel(()->Channel(1)),
        RemoteChannel(()->Channel(1)), false)
    end
end

function psracos_init_sample_set!(psracos::PSRacos, ub)
    rc = psracos.rc
    classifier = RacosClassification(rc.objective.dim, rc.positive_data,
        rc.negative_data, ub=ub)
    mixed_classification(classifier)
    for i = 1:(psracos.evaluation_server_num)
        if rand(rng, Float64) < rc.parameter.probability
            solution, distinct_flag = distinct_sample_classifier(rc, classifier, data_num=rc.parameter.train_size)
        else
            solution, distinct_flag = distinct_sample(rc, rc.objective.dim)
        end
        # sol_print(solution)
        put!(psracos.sample_set, solution)
    end
end


text:
struct Squares
end
Base.iterate(S::Squares, state=1) = (state*state, state+1)

using DynamicIterators
# Assign iterate number as key
T = TimeLift(Squares())
# Create NextKey control slot
C = attach(NextKeys, T)
# Bind iterator 4:2:8 to the control slot and collect
collect(bind(4:2:8, C))


text:
## test_covariance_functions.jl : tests for covariance_functions.jl

## Matérn ##
@testset "Mat\u00E9rn covariance        " begin

m = Matern(0.1,2.5)
@test isa(m,Matern)
m = Matern(1,2.5)
@test isa(m,Matern)
m = Matern(0.1,2)
@test isa(m,Matern)
m = Matern(1//10,2)
@test isa(m,Matern)
m = Matern(0.1,2,p=1.)
@test isa(m,Matern)
m = Matern(0.1,2,σ=1.)
@test isa(m,Matern)
m = Matern(0.1,2,σ=2.)
@test isa(m,Matern)
m = Matern(0.1,2,σ=1.,p=2)
@test isa(m,Matern)

@test_throws DomainError Matern(-0.1,2.5)
@test_throws DomainError Matern(0.1,-2.5)
@test_throws DomainError Matern(0.1,2.5,p=-1)
@test_throws DomainError Matern(0.1,2.5,σ=-1)
@test_throws DomainError Matern(0.1,2.5,p=0.5)
@test_throws DomainError Matern(0.1,2.5,p=Inf)
end

## Exponential ##
@testset "exponential covariance   " begin

e = Exponential(0.1)
@test isa(e,Exponential)
e = Exponential(1)
@test isa(e,Exponential)
e = Exponential(1,p=1.)
@test isa(e,Exponential)
e = Exponential(.1,p=1)
@test isa(e,Exponential)
e = Exponential(.1,σ=10)
@test isa(e,Exponential)
e = Exponential(.1,σ=2.,p=1)
@test isa(e,Exponential)

@test_throws DomainError Exponential(-0.1,p=1)
@test_throws DomainError Exponential(0.1,σ=-2)
@test_throws DomainError Exponential(0.1,p=0.9)
@test_throws DomainError Exponential(0.1,p=Inf)

end

## SquaredExponential ##
@testset "Gaussian covariance      " begin

s = SquaredExponential(0.1)
@test isa(s,Gaussian)
s = Gaussian(0.1)
@test isa(s,SquaredExponential)
s = SquaredExponential(1)
@test isa(s,SquaredExponential)
s = SquaredExponential(1,p=1.)
@test isa(s,SquaredExponential)
s = SquaredExponential(.1,σ=1.)
@test isa(s,SquaredExponential)
s = SquaredExponential(.1,σ=10.,p=2)
@test isa(s,SquaredExponential)
s = SquaredExponential(.1,p=10.,σ=2)
@test isa(s,SquaredExponential)

@test_throws DomainError SquaredExponential(-0.1)
@test_throws DomainError SquaredExponential(0.1,p=0.5)
@test_throws DomainError SquaredExponential(0.1,σ=-2)
@test_throws DomainError Gaussian(0.1,p=Inf)
end

## Linear ##
@testset "linear covariance        " begin

l = Linear(0.1)
@test isa(l,Linear)
l = Linear(1)
@test isa(l,Linear)
l = Linear(1,p=1.)
@test isa(l,Linear)
l = Linear(.1,p=1)
@test isa(l,Linear)
l = Linear(.1,σ=10)
@test isa(l,Linear)
l = Linear(.1,σ=2.,p=1)
@test isa(l,Linear)

@test_throws DomainError Linear(-0.1,p=1)
@test_throws DomainError Linear(0.1,σ=-2)
@test_throws DomainError Linear(0.1,p=0.9)
@test_throws DomainError Linear(0.1,p=Inf)

end

## Spherical ##
@testset "spherical covariance     " begin

l = Spherical(0.1)
@test isa(l,Spherical)
l = Spherical(1)
@test isa(l,Spherical)
l = Spherical(1,p=1.)
@test isa(l,Spherical)
l = Spherical(.1,p=1)
@test isa(l,Spherical)
l = Spherical(.1,σ=10)
@test isa(l,Spherical)
l = Spherical(.1,σ=2.,p=1)
@test isa(l,Spherical)

@test_throws DomainError Spherical(-0.1,p=1)
@test_throws DomainError Spherical(0.1,σ=-2)
@test_throws DomainError Spherical(0.1,p=0.9)
@test_throws DomainError Spherical(0.1,p=Inf)

end

## Whittle ##
@testset "Whittle covariance        " begin

l = Whittle(0.1)
@test isa(l,Whittle)
l = Whittle(1)
@test isa(l,Whittle)
l = Whittle(1,p=1.)
@test isa(l,Whittle)
l = Whittle(.1,p=1)
@test isa(l,Whittle)
l = Whittle(.1,σ=10)
@test isa(l,Whittle)
l = Whittle(.1,σ=2.,p=1)
@test isa(l,Whittle)

@test_throws DomainError Whittle(-0.1,p=1)
@test_throws DomainError Whittle(0.1,σ=-2)
@test_throws DomainError Whittle(0.1,p=0.9)
@test_throws DomainError Whittle(0.1,p=Inf)

end

## AnisotropicExponential ##
@testset "anisotropic exponential  " begin

A = [1 0; 0 1]
s = AnisotropicExponential(A)
@test isa(s,AnisotropicExponential)
A = [1 0.8 0.8; 0.8 1 0.8; 0.8 0.8 1]
s = AnisotropicExponential(A)
@test isa(s,AnisotropicExponential)
A = [1000 0; 0 1000]
s = AnisotropicExponential(A,σ=2)
@test isa(s,AnisotropicExponential)
s = AnisotropicExponential(A,σ=2.)
@test isa(s,AnisotropicExponential)

A = [10 10; -10 10]
@test_throws DomainError AnisotropicExponential(A)
A = [10 0; 0 10]
@test_throws DomainError AnisotropicExponential(A,σ=-1)
end

## CovarianceFunction ##
@testset "covariance function      " begin

m = Matern(0.1,0.5)
c = CovarianceFunction(2,m)
@test isa(c,CovarianceFunction)

e = Exponential(0.1)
c = CovarianceFunction(1,e)
@test isa(c,CovarianceFunction)

s = SquaredExponential(0.1)
c = CovarianceFunction(2,s)
@test isa(c,CovarianceFunction)

@test_throws DomainError CovarianceFunction(0,m)
@test_throws DomainError CovarianceFunction(-1,e)
end


text:
using VideoIO
using ImageInTerminal

ImageInTerminal.use_24bit()
f=opencamera()
img=read(f)


text:
# This file was generated, do not modify it. # hide
X1 = coerce(X0, :OverallQual => Continuous,
            :GarageCars => Continuous,
            :YearRemodAdd => Continuous,
            :YearBuilt => Continuous);

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
using LinearAlgebra
using GPR

include("utils/utils.jl")
include("utils/data/datasets.jl")
include("utils/data/transformations.jl")
include("utils/predictdynamics.jl")
include("parallel/core.jl")

include("maximal_coordinates/P1noise.jl")
include("maximal_coordinates/P2noise.jl")
include("maximal_coordinates/CPnoise.jl")
include("maximal_coordinates/FBnoise.jl")

include("minimal_coordinates/P1noise.jl")
include("minimal_coordinates/P2noise.jl")
include("minimal_coordinates/CPnoise.jl")
include("minimal_coordinates/FBnoise.jl")

include("baseline.jl")

function loadcheckpoint_or_defaults(_loadcheckpoint::Bool)
    nprocessed = 0
    kstep_mse = []
    if _loadcheckpoint
        checkpoint = loadcheckpoint("noisy")["EXPERIMENT_ID"]
        nprocessed = checkpoint["nprocessed"]
        kstep_mse = checkpoint["kstep_mse"]
    end
    return nprocessed, kstep_mse
end

function expand_config(EXPERIMENT_ID::String, nsamples::Integer, config::Dict, dfs::Tuple{DataFrame, DataFrame}; _loadcheckpoint::Bool = false)
    EXPERIMENT_ID = EXPERIMENT_ID*string(nsamples)
    nprocessed, kstep_mse = loadcheckpoint_or_defaults(_loadcheckpoint)
    params = Vector{Float64}()  # declare in outer scope
    open(joinpath(dirname(@__FILE__), "config", "config.json"),"r") do f
        params = JSON.parse(f)[EXPERIMENT_ID]
    end
    Σ = Dict("x" => 1e-3, "q" => 1e-3, "v" => 1e-3, "ω" => 1e-3)
    # Σ = Dict("x" => 0, "q" => 0, "v" => 0, "ω" => 0)  # Also disabled in experiments! Activate there as well!
    config = Dict("EXPERIMENT_ID" => EXPERIMENT_ID,
                  "Σ" => Σ,
                  "params" => params,
                  "nruns" => config["nruns"],
                  "Δtsim" => config["Δtsim"],
                  "nsamples" => nsamples,
                  "testsamples" => config["testsamples"],
                  "simsteps" => config["simsteps"],
                  "nprocessed" => nprocessed,
                  "kstep_mse" => kstep_mse,
                  "projectionerror" => Vector{Float64}(),
                  "mechanismlock" => ReentrantLock(),
                  "paramlock" => ReentrantLock(),
                  "resultlock" => ReentrantLock(),
                  "checkpointlock" => ReentrantLock(),
                  "datasets" => dfs)
    return config
end

config = Dict("nruns" => 100, "Δtsim" => 0.0001, "testsamples" => 100, "simsteps" => 20)

samplesizes = [2, 4, 8, 16, 32, 64, 128, 256, 512]  # [2, 4, 8, 16, 32, 64, 128, 256, 512,
dfkeys = ["FB"]  # "P1", "P2", "CP", "FB"
@info "Main Thread: loading datasets $dfkeys"
datasets = Dict(key => loaddatasets(key) for key in dfkeys)
@info "Main Thead: finished loading datasets"

#=
@info "Processing pure variational integrator experiments"
parallelsim(experimentVarIntP1, expand_config("P1_MIN", 2, config, datasets["P1"]), idmod = "VI")
parallelsim(experimentVarIntP2, expand_config("P2_MIN", 2, config, datasets["P2"]), idmod = "VI")
parallelsim(experimentVarIntCP, expand_config("CP_MIN", 2, config, datasets["CP"]), idmod = "VI")
parallelsim(experimentVarIntFB, expand_config("FB_MIN", 2, config, datasets["FB"]), idmod = "VI")

for nsamples in samplesizes
    @info "Processing P1 experiments"
    parallelsim(experiment_p1_mz_max, expand_config("P1_MAX", nsamples, config, datasets["P1"]))
    parallelsim(experiment_p1_mz_min, expand_config("P1_MIN", nsamples, config, datasets["P1"]))
    parallelsim(experiment_p1_mz_min_sin, expand_config("P1_MIN", nsamples, config, datasets["P1"]), idmod = "sin")
    parallelsim(experiment_p1_md_max, expand_config("P1_MAX", nsamples, config, datasets["P1"]), idmod = "MD")
    parallelsim(experiment_p1_md_min, expand_config("P1_MIN", nsamples, config, datasets["P1"]), idmod = "MD")
    parallelsim(experiment_p1_md_min_sin, expand_config("P1_MIN", nsamples, config, datasets["P1"]), idmod = "MDsin")
end

for nsamples in samplesizes
    @info "Processing P2 experiments"
    parallelsim(experiment_p2_mz_max, expand_config("P2_MAX", nsamples, config, datasets["P2"]))
    parallelsim(experiment_p2_mz_min, expand_config("P2_MIN", nsamples, config, datasets["P2"]))
    parallelsim(experiment_p2_mz_min_sin, expand_config("P2_MIN", nsamples, config, datasets["P2"]), idmod = "sin")
    parallelsim(experiment_p2_md_max, expand_config("P2_MAX", nsamples, config, datasets["P2"]), idmod = "MD")
    parallelsim(experiment_p2_md_min, expand_config("P2_MIN", nsamples, config, datasets["P2"]), idmod = "MD")
    parallelsim(experiment_p2_md_min_sin, expand_config("P2_MIN", nsamples, config, datasets["P2"]), idmod = "MDsin")
end

for nsamples in samplesizes
    @info "Processing CP experiments"
    parallelsim(experiment_cp_mz_max, expand_config("CP_MAX", nsamples, config, datasets["CP"]))
    parallelsim(experiment_cp_mz_min, expand_config("CP_MIN", nsamples, config, datasets["CP"]))
    parallelsim(experiment_cp_mz_min_sin, expand_config("CP_MIN", nsamples, config, datasets["CP"]), idmod = "sin")
    parallelsim(experiment_cp_md_max, expand_config("CP_MAX", nsamples, config, datasets["CP"]), idmod = "MD")
    parallelsim(experiment_cp_md_min, expand_config("CP_MIN", nsamples, config, datasets["CP"]), idmod = "MD")
    parallelsim(experiment_cp_md_min_sin, expand_config("CP_MIN", nsamples, config, datasets["CP"]), idmod = "MDsin")
end
=#
for nsamples in samplesizes
    @info "Processing FB experiments"
    parallelsim(experiment_fb_mz_max, expand_config("FB_MAX", nsamples, config, datasets["FB"]))
    parallelsim(experiment_fb_mz_min, expand_config("FB_MIN", nsamples, config, datasets["FB"]))
    parallelsim(experiment_fb_mz_min_sin, expand_config("FB_MIN", nsamples, config, datasets["FB"]), idmod = "sin")
    parallelsim(experiment_fb_md_max, expand_config("FB_MAX", nsamples, config, datasets["FB"]), idmod = "MD")
    parallelsim(experiment_fb_md_min, expand_config("FB_MIN", nsamples, config, datasets["FB"]), idmod = "MD")
    parallelsim(experiment_fb_md_min_sin, expand_config("FB_MIN", nsamples, config, datasets["FB"]), idmod = "MDsin")
end


text:
include("run_scf_and_compare.jl")
include("testcases.jl")

function run_silicon_pbe(T ;Ecut=5, grid_size=15, spin_polarization=:none, kwargs...)
    # These values were computed using ABINIT with the same kpoints as testcases.jl
    # and Ecut = 25
    ref_pbe = [
        [-0.181210259413818, 0.258840553222639, 0.258840553225549, 0.258840553228459,
          0.351692348652324, 0.351692348656259, 0.351692348660193, 0.380606400669216,
          0.540705881744348, 0.540705883460555],
        [-0.130553299114991, 0.062256443775155, 0.221871391287580, 0.221871391290802,
          0.322398722411882, 0.386194327436667, 0.386194327439986, 0.546859898649217,
          0.550571701390781, 0.550571701394327],
        [-0.111170738096744, 0.074494899973125, 0.169461730083372, 0.169461730088140,
          0.284305392082236, 0.330468937070505, 0.524509288492752, 0.524509288496625,
          0.616964090764029, 0.619623658242765],
        [-0.061054203629684, 0.009700769243041, 0.095769985640881, 0.180784778430457,
          0.315000287382235, 0.471042322838057, 0.495281775946584, 0.517469860611792,
          0.530124341745161, 0.539044739392045],
    ]
    ref_etot = -7.854477356672080

    fft_size = grid_size * ones(3)
    Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.atnum, functional="pbe", family="hgh"))
    magmoms = spin_polarization == :collinear ? [Si => zeros(2)] : []
    model = model_DFT(Array{T}(silicon.lattice), [Si => silicon.positions], [:gga_x_pbe, :gga_c_pbe],
                      spin_polarization=spin_polarization, magnetic_moments=magmoms)
    basis = PlaneWaveBasis(model, Ecut, silicon.kcoords, silicon.ksymops; fft_size=fft_size)

    spin_polarization == :collinear && (ref_pbe = vcat(ref_pbe, ref_pbe))
    run_scf_and_compare(T, basis, ref_pbe, ref_etot; kwargs...)
end

@testset "Silicon PBE (small, Float64)" begin
    run_silicon_pbe(Float64, Ecut=7, test_tol=0.03, n_ignored=0, grid_size=17)
end

if !isdefined(Main, :FAST_TESTS) || !FAST_TESTS
    @testset "Silicon PBE (large, Float64)" begin
        run_silicon_pbe(Float64, Ecut=25, test_tol=1e-5, n_ignored=0,
                        grid_size=33, scf_tol=1e-8)
    end
end

@testset "Silicon PBE (small, collinear spin)" begin
    run_silicon_pbe(Float64, Ecut=7, test_tol=0.03, n_ignored=0, grid_size=17, spin_polarization=:collinear)
end

if !isdefined(Main, :FAST_TESTS) || !FAST_TESTS
    @testset "Silicon PBE (large, collinear spin)" begin
        run_silicon_pbe(Float64, Ecut=25, test_tol=1e-5, n_ignored=0, grid_size=33,
                        scf_tol=1e-8, spin_polarization=:collinear)
    end
end


text:
if basename(pwd()) == "aoc"
    cd("2017/15")
end

struct Generator
    initialvalue::Int
    multiplier::Int
    divider::Int
end
GeneratorA(init::Int, div = 1) = Generator(init, 16807, div)
GeneratorB(init::Int, div = 1) = Generator(init, 48271, div)

function Base.iterate(gen::Generator, value = gen.initialvalue)
    value = (value * gen.multiplier) % 2147483647
    while !iszero(value % gen.divider)
        value = (value * gen.multiplier) % 2147483647
    end
    (value, value)
end
Base.eltype(::Type{Generator}) = Int
Base.IteratorSize(::Type{Generator}) = Base.IsInfinite()

function part1(ainit::Int, binit::Int)
    count(
        a & 0xffff == b & 0xffff
        for (a, b) in Iterators.take(
            zip(
                GeneratorA(ainit),
                GeneratorB(binit)),
            40_000_000))
end

@assert part1(65, 8921) == 588
part1(289, 629)

function part2(ainit::Int, binit::Int)
    count(
        a & 0xffff == b & 0xffff
        for (a, b) in Iterators.take(
            zip(
                GeneratorA(ainit, 4),
                GeneratorB(binit, 8)),
            5_000_000))
end

@assert part2(65, 8921) == 309
part2(289, 629)

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
@testset "DataRS tests" begin show_testset_section() #Scope for test data

using MDDatasets
using Test

#==Input data
===============================================================================#
#get_ydata: Generates sine data from input parameters
get_ydata(t, tbit, vdd, trise) = sin(2pi*t/tbit)*(trise/tbit)+vdd

t = DataF1((1:2)*1e-9)

#Construct DataRS object from parametric sweeps:
data = fill(DataRS, PSweep("tbit", [1, 3] * 1e-9)) do tbit
	#Need to explicitly specify DataRS{DataF1} as type for inner-most sweep:
	fill(DataRS{DataF1}, PSweep("VDD", 0.9 * [0.9, 1])) do vdd
		trise = 0.1*tbit
		return get_ydata(t, tbit, vdd, trise)
	end
end

#Construct DataRS object from parametric sweeps:
_datagen1(indep::Bool) = fill(DataRS, PSweep("A", Float64[1, 3])) do A
	offset = indep ? 0 : A
	#Need to explicitly specify DataRS{DataFloat} as type for inner-most sweep:
	fill(DataRS{DataFloat}, PSweep("B", Float64[3, 5, 6] .+ offset)) do B
		return 10A + B
	end
end

#==Tests
===============================================================================#
@testset "Validate _datamatch()" begin show_testset_description()
	datars_dep=_datagen1(false)
	other = datars_dep*1.0
	@test !_datamatch(datars_dep, other+1) #Make sure _datamatch works
	@test _datamatch(datars_dep, other)
end

@testset "Construction of DataRS objects" begin show_testset_description()
	#TODO: write tests to verify structure
	@info("Visually confirm structure of \"data\":")
	@show data

	@test data.sweep.id == "tbit" #should be outer-most loop
	@test data.elem[1].sweep.id == "VDD" #should be inner-most loop

	#Simple (independent) parameter sweep:
	swlist = PSweep[
		PSweep("v1", [1,2])
		PSweep("v2", [4,5])
	]

	dhrs = fill(t, DataRS, swlist)
	@test dhrs.sweep.id == "v1" #should be outer-most loop
	@test dhrs.elem[1].sweep.id == "v2" #should be inner-most loop

	dhrs_0 = zeros(DataRS, swlist)
	dhrs_1 = ones(DataRS, swlist)
	dhrs_2 = fill(2, DataRS, swlist)

	@test _datamatch(dhrs_2-dhrs_1, dhrs_1)
	@test _datamatch(dhrs_2-dhrs_1, dhrs_1)
	@test _datamatch(dhrs_2-2*dhrs_1, dhrs_0)

	#TODO: Check that we actually have zeros, etc.
end

@testset "Access parameter data" begin show_testset_description()
	#Extract value of tbit (first sweep):
	data_tbit = parameter(data, "tbit")

	#Get sweep values:
	sw_tbit = data.sweep.v
	sw_vdd = data.elem[1].sweep.v

	for itbit in 1:length(sw_tbit), ivdd in 1:length(sw_vdd)
		@test data_tbit.elem[itbit].elem[ivdd] == sw_tbit[itbit]
	end
end

@testset "Broacasting operations on DataRS" begin show_testset_description()
	#Get sweep values:
	sw_tbit = data.sweep.v
	sw_vdd = data.elem[1].sweep.v

	_v = 0*data+5 #Same dimensionality as data
	@testset "_v = 0*data+5" begin
		v5 = zeros(length(t)) .+ 5

		for itbit in 1:length(sw_tbit), ivdd in 1:length(sw_vdd)
			_elem = _v.elem[itbit].elem[ivdd]
			@test _elem.x == t.x
			@test _elem.y == v5
		end
	end

	_vshift = data + _v
	@testset "_vshift = data + _v" begin
		for itbit in 1:length(sw_tbit), ivdd in 1:length(sw_vdd)
			_elem = _vshift.elem[itbit].elem[ivdd]
			@test _elem.x == t.x
			@test _elem.y == data.elem[itbit].elem[ivdd].y .+ 5
		end
	end

	_m1 = maximum(data)
	@testset "_m1 = maximum(data)" begin
		for itbit in 1:length(sw_tbit), ivdd in 1:length(sw_vdd)
			_elem = _m1.elem[itbit].elem[ivdd]
			@test typeof(_elem) == Float64
			@test _elem == maximum(data.elem[itbit].elem[ivdd].y)
		end
	end

	_delta = data - _m1
	@testset "_delta = data - _m1" begin
		for itbit in 1:length(sw_tbit), ivdd in 1:length(sw_vdd)
			_elem = _delta.elem[itbit].elem[ivdd]
			@test _elem.x == t.x
			_y = data.elem[itbit].elem[ivdd].y
			@test _elem.y == _y .- maximum(_y)
		end
	end

	_m2 = maximum(_m1) #Collapse inner-most dimension again
	@testset "_m2 = maximum(_m1)" begin
		for itbit in 1:length(sw_tbit)
			_elem = _m2.elem[itbit]
			@test typeof(_elem) == Float64
			@test _elem == maximum(_m1.elem[itbit])
		end
	end

	_m3 = maximum(_m2) #Collapse inner-most dimension again
	@testset "_m3 = maximum(_m2)" begin
		_elem = _m3
		m2val = [_m2.elem[i] for i in 1:length(_m2)] #collect values in m2
		@test typeof(_elem) == Float64
		@test _elem == maximum(m2val)
	end
end

@testset "Conversion (DataHR <=> DataRS)" begin show_testset_description()
	datars_indep=_datagen1(true)
	datars_dep=_datagen1(false)

	datahr = convert(DataHR, datars_indep)
		#TODO: test something
	@test_throws ArgumentError convert(DataHR, datars_dep)

	datars_2 = convert(DataRS, datahr)
	@test _datamatch(datars_2, datars_2)

	try
		datars_2 - datars_indep
	catch
		@warn("Cannot operate with DataRS re-created through DataHR because of PSweep comparison & way fill!(::DataRS) works.")
	end
end

end


text:
"""
Linear Chain - CRF Layer.

For input sequence `x`,
predicts the most probable tag sequence `y`,
over the set of all possible tagging sequences `Y`.

In this CRF, two kinds of potentials are defined,
emission and Transition.
"""
mutable struct CRF{S}
    W::S        # Transition Scores
    n::Int      # Num Labels
end

"""
Second last index for start tag,
last one for stop tag .
"""
function CRF(n::Integer)
    W = rand(Float32, n + 2, n + 2)
    W[:, n + 1] .= -10000
    W[n + 2, :] .= -10000

    return CRF(W, n)
end

@functor CRF

function Base.show(io::IO, c::CRF)
    print(io, "CRF with ", c.n + 2, " distinct tags (including START and STOP tags).")
end

function (a::CRF)(x_seq, init_α)
    viterbi_decode(a, x_seq, init_α)
end


text:
using SnoopCompile, HTTP
ENV["PKGSERVER_GENERATING_PRECOMPILE"] = "true"

# Start tasks to hit all sorts of resources that are useful
@async begin
    for res in ["meta", "meta/stats", "meta/siblings", "meta/parents",
                "robots.txt",
                "registries", "registry/23338594-aafe-5451-b93e-139f81909106/e5ed3655ffd697a8ea2a3fdff0cd14022144f497",
                "package/78c3b35d-d492-501b-9361-3d52fe80e533/916b850daad0d46b8c71f65f719c49957e9513ed",
                "artifact/a4e3df8d4997a3352c87cb1b8a8b86ac6b40ddba"]
        while true
            sleep(0.5)
            try
                r = HTTP.get("http://localhost:8000/$(res)")
                if r.status == 200
                    @info("Successfully got a /$(res) response!")
                    break
                end
            catch
            end
        end
    end

    # Once we've got all these goodies, kill the server
    close(PkgServer.config.listen_server[])
end

inf_timing = @snoopi tmin=0.001 begin
    using PkgServer

    mktempdir() do root
        try
            PkgServer.start(storage_root=root)
        catch e
            if !isa(e, Base.IOError) || e.code != -Base.Libc.ECONNABORTED
                rethrow(e)
            end
        end
    end
end

pc = SnoopCompile.parcel(inf_timing)
mktempdir() do dir
    SnoopCompile.write(dir, pc)
    cp(joinpath(dir, "precompile_PkgServer.jl"),
       joinpath(@__DIR__, "precompile.jl"); force=true)
end

max_stars_repo_name: StefanKarpinski/PkgServer

max_issues_repo_name: StefanKarpinski/PkgServer

max_forks_repo_name: StefanKarpinski/PkgServer


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

using Markdown
using InteractiveUtils

# ╔═╡ 413f4d5a-d0f0-4a58-a08d-7664b29057d8
# hideall

using Genie

# ╔═╡ bac482d8-4781-11ec-306b-1f4298a0cb54
md"""
# Developing a simple API backend

Genie makes it very easy to quickly set up a REST API backend. All it takes is a few lines of code:
"""

# ╔═╡ 329f9db4-2af8-4259-901f-fdcacba830ad
md"""

```julia
using Genie
import Genie.Router: route
import Genie.Renderer.Json: json

Genie.config.run_as_server = true

route("/") do
  (:message => "Hi there!") |> json
end

Genie.startup()
```
"""

# ╔═╡ 26e89533-8922-41ad-b49f-0cea678cf296
# hideall

import Genie.Router: route

# ╔═╡ 3643e238-b83b-469d-8e37-6b0b86c59f78
# hideall

import Genie.Renderer.Json: json

# ╔═╡ ee5094b5-24ac-4f81-9636-6f623189124d
# hideall

Genie.config.run_as_server = true;

# ╔═╡ 9f830898-299e-46aa-a4bd-f7e225b4793b
# hideall

route("/") do
  (:message => "Hi there!") |> json
end;

# ╔═╡ fc1de4cd-80a7-40d6-9a1a-2719ed170508
md"""
The key bit here is `Genie.config.run_as_server = true`. This will start the server synchronously so the `startup()` function won't return.
This endpoint can be run directly from the command line - if say, you save the code in a `rest.jl` file:
"""

# ╔═╡ f4568f6a-95f4-447d-af90-9c4903fa65be
md"""

```julia
$ julia rest.jl
```
"""

# ╔═╡ b804a827-4f41-4fd0-a9e4-70fd25886bdc
md"""
## Accepting JSON payloads

One common requirement when exposing APIs is to accept POST payloads. That is, requests over POST, with a request body, usually as a JSON encoded object. We can build an echo service like this:

```julia
using Genie, Genie.Router, Genie.Renderer.Json, Genie.Requests
using HTTP

route("/echo", method = POST) do
  message = jsonpayload()
  (:echo => (message["message"] * " ") ^ message["repeat"]) |> json
end

route("/send") do
  response = HTTP.request("POST", "http://localhost:8000/echo", [("Content-Type", "application/json")], \"""{"message":"hello", "repeat":3}\""")

  response.body |> String |> json
end

Genie.startup(async = false)
```
"""

# ╔═╡ fe97ffdb-22ed-456b-8269-017f7a0b3637
# hideall

# using HTTP

# ╔═╡ a10279fd-2f17-4964-ae7c-811543194aa1
# hideall

# write("rest.jl", """using Genie, Genie.Router, Genie.Renderer.Json, Genie.Requests
# using HTTP

# route("/echo", method = POST) do
#   message = jsonpayload()
#   (:echo => (message["message"] * " ") ^ message["repeat"]) |> json
# end

# route("/send") do
#   response = HTTP.request("POST", "http://localhost:8000/echo", [("Content-Type", "application/json")], \"""{"message":"hello", "repeat":3}\""")

#   response.body |> String |> json
# end

# Genie.startup(async = true)""");

# ╔═╡ f7c15c14-73fb-440c-a79f-0ab34ff6eac6
# hideall

# include("rest.jl");

# ╔═╡ 990d6c09-eba8-4741-b6fd-8eb39cc9a30c
# hideall

# Genie.Router.routes();

# ╔═╡ 8a936d9b-5652-4a6a-bffa-4352affb4e00
md"""
check Web Server starting at **http://127.0.0.1:8000**
"""

# ╔═╡ 75df78b7-fbd8-4ea1-baba-f5344a412f8b
md"""

Here we define two routes, `/send` and `/echo`. The `send` route makes a `HTTP` request over `POST` to `/echo`, sending a JSON payload with two values, `message` and `repeat`.
In the `/echo` route, we grab the JSON payload using the `Requests.jsonpayload()` function, extract the values from the JSON object, and output the `message` value repeated for a number of times equal to the `repeat` value.
"""

# ╔═╡ 992ac7bc-f90e-4711-aa87-5cd5554fac3e
md"""
If you run the code, the output should be

```javascript
{
  echo: "hello hello hello "
}
```
"""

# ╔═╡ 2c85bbfa-5892-4e70-8070-b017559ccbbb
md"""
If the payload contains invalid JSON, the `jsonpayload` will be set to `nothing`. You can still access the raw payload by using the `Requests.rawpayload()` function.
You can also use `rawpayload` if for example the type of request/payload is not JSON.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Genie = "c43c736e-a2d1-11e8-161f-af95117fbd1e"

[compat]
Genie = "~4.2.0"
"""

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

[[ArgParse]]
deps = ["Logging", "TextWrap"]
git-tree-sha1 = "3102bce13da501c9104df33549f511cd25264d7d"
uuid = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
version = "1.1.4"

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

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

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

[[CSTParser]]
deps = ["Tokenize"]
git-tree-sha1 = "f9a6389348207faf5e5c62cbc7e89d19688d338a"
uuid = "00ebfdb7-1f24-5e51-bd34-a7502290713f"
version = "3.3.0"

[[CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "9aa8a5ebb6b5bf469a7e0e2b5202cf6f8c291104"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.6"

[[CommonMark]]
deps = ["Crayons", "JSON", "URIs"]
git-tree-sha1 = "393ac9df4eb085c2ab12005fc496dae2e1da344e"
uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
version = "0.8.3"

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

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

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

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

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

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

[[EzXML]]
deps = ["Printf", "XML2_jll"]
git-tree-sha1 = "0fa3b52a04a4e210aeb1626def9c90df3ae65268"
uuid = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615"
version = "1.1.0"

[[FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "618835ab81e4a40acf215c98768978d82abc5d97"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.16"

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

[[GMP_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"

[[Genie]]
deps = ["ArgParse", "Dates", "Distributed", "EzXML", "FilePathsBase", "HTTP", "HttpCommon", "Inflector", "JSON3", "JuliaFormatter", "Logging", "Markdown", "MbedTLS", "Millboard", "Nettle", "OrderedCollections", "Pkg", "REPL", "Reexport", "Revise", "SHA", "Serialization", "Sockets", "UUIDs", "Unicode", "VersionCheck", "YAML"]
git-tree-sha1 = "9ec6a1b969aaf8aa6d6b9513055fba3effa7671e"
uuid = "c43c736e-a2d1-11e8-161f-af95117fbd1e"
version = "4.2.0"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

[[HttpCommon]]
deps = ["Dates", "Nullables", "Test", "URIParser"]
git-tree-sha1 = "46313284237aa6ca67a6bce6d6fbd323d19cff59"
uuid = "77172c1b-203f-54ac-aa54-3f1198fe9f90"
version = "0.5.0"

[[Inflector]]
deps = ["Unicode"]
git-tree-sha1 = "8555b54ddf27806b070ce1d1cf623e1feb13750c"
uuid = "6d011eab-0732-4556-8808-e463c76bf3b6"
version = "1.0.1"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

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

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[JSON3]]
deps = ["Dates", "Mmap", "Parsers", "StructTypes", "UUIDs"]
git-tree-sha1 = "7d58534ffb62cd947950b3aa9b993e63307a6125"
uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
version = "1.9.2"

[[JuliaFormatter]]
deps = ["CSTParser", "CommonMark", "DataStructures", "Pkg", "Tokenize"]
git-tree-sha1 = "10c95cebcfa37c1f510a726c90886db4745e1238"
uuid = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
version = "0.15.11"

[[JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "e273807f38074f033d94207a201e6e827d8417db"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.8.21"

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

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

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

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

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

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

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

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

[[LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "491a883c4fef1103077a7f648961adbf9c8dd933"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.1.2"

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

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

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

[[Millboard]]
git-tree-sha1 = "ea6a5b7e56e76d8051023faaa11d91d1d881dac3"
uuid = "39ec1447-df44-5f4c-beaa-866f30b4d3b2"
version = "0.2.5"

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

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

[[Nettle]]
deps = ["Libdl", "Nettle_jll"]
git-tree-sha1 = "a68340b9edfd98d0ed96aee8137cb716ea3b6dea"
uuid = "49dea1ee-f6fa-5aa6-9a11-8816cee7d4b9"
version = "0.5.1"

[[Nettle_jll]]
deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "eca63e3847dad608cfa6a3329b95ef674c7160b4"
uuid = "4c82536e-c426-54e4-b420-14f461c4ed8b"
version = "3.7.2+0"

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

[[Nullables]]
git-tree-sha1 = "8f87854cc8f3685a60689d8edecaa29d2251979b"
uuid = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
version = "1.0.0"

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

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

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

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

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

[[ProgressMeter]]
deps = ["Distributed", "Printf"]
git-tree-sha1 = "afadeba63d90ff223a6a48d2009434ecee2ec9e8"
uuid = "92933f4c-e287-5a05-a399-4b506db050ca"
version = "1.7.1"

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

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

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

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

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

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

[[Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

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

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

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

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

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

[[StringEncodings]]
deps = ["Libiconv_jll"]
git-tree-sha1 = "50ccd5ddb00d19392577902f0079267a72c5ab04"
uuid = "69024149-9ee7-55f6-a4c4-859efe599b68"
version = "0.3.5"

[[StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d24a825a95a6d98c385001212dc9020d609f2d4f"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.8.1"

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

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

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

[[TextWrap]]
git-tree-sha1 = "9250ef9b01b66667380cf3275b3f7488d0e25faf"
uuid = "b718987f-49a8-5099-9789-dcd902bef87d"
version = "1.0.1"

[[Tokenize]]
git-tree-sha1 = "0952c9cee34988092d73a5708780b3917166a0dd"
uuid = "0796e94c-ce3b-5d07-9a54-7f471281c624"
version = "0.5.21"

[[URIParser]]
deps = ["Unicode"]
git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.1"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

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

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

[[UrlDownload]]
deps = ["HTTP", "ProgressMeter"]
git-tree-sha1 = "05f86730c7a53c9da603bd506a4fc9ad0851171c"
uuid = "856ac37a-3032-4c1c-9122-f86d88358c8b"
version = "1.0.0"

[[VersionCheck]]
deps = ["Dates", "JSON3", "Logging", "Pkg", "Random", "Scratch", "UrlDownload"]
git-tree-sha1 = "89ef2431dd59344ebaf052d0737205854ded0c62"
uuid = "a637dc6b-bca1-447e-a4fa-35264c9d0580"
version = "0.2.0"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[YAML]]
deps = ["Base64", "Dates", "Printf", "StringEncodings"]
git-tree-sha1 = "3c6e8b9f5cdaaa21340f841653942e1a6b6561e5"
uuid = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
version = "0.4.7"

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

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

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

# ╔═╡ Cell order:
# ╟─bac482d8-4781-11ec-306b-1f4298a0cb54
# ╟─329f9db4-2af8-4259-901f-fdcacba830ad
# ╠═413f4d5a-d0f0-4a58-a08d-7664b29057d8
# ╠═26e89533-8922-41ad-b49f-0cea678cf296
# ╠═3643e238-b83b-469d-8e37-6b0b86c59f78
# ╠═ee5094b5-24ac-4f81-9636-6f623189124d
# ╠═9f830898-299e-46aa-a4bd-f7e225b4793b
# ╟─fc1de4cd-80a7-40d6-9a1a-2719ed170508
# ╟─f4568f6a-95f4-447d-af90-9c4903fa65be
# ╟─b804a827-4f41-4fd0-a9e4-70fd25886bdc
# ╠═fe97ffdb-22ed-456b-8269-017f7a0b3637
# ╠═a10279fd-2f17-4964-ae7c-811543194aa1
# ╠═f7c15c14-73fb-440c-a79f-0ab34ff6eac6
# ╠═990d6c09-eba8-4741-b6fd-8eb39cc9a30c
# ╟─8a936d9b-5652-4a6a-bffa-4352affb4e00
# ╟─75df78b7-fbd8-4ea1-baba-f5344a412f8b
# ╟─992ac7bc-f90e-4711-aa87-5cd5554fac3e
# ╟─2c85bbfa-5892-4e70-8070-b017559ccbbb
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using DelimitedFiles
using Flux

# load a specialized model format
function loadPyNNTxt(dest::AbstractString)
  mw = Vector{Array{Float32}}()
  @show files = readdir(dest)
  for f in files
    push!(mw, readdlm(joinpath(dest,f)))
  end
  return mw
end

# convenience function to load specific model format from tensorflow
function loadPose2OdoNNModelIntoFlux(dest::AbstractString; pad::Bool=false)
  weights = loadPyNNTxt(dest::AbstractString)
  buildPose2OdoNN_01_FromWeights(weights, pad=pad)
end

# function buildPose2OdoNN_01_FromWeightsGPU(pywe)
#   buildPose2OdoNN_01_FromElements(collect(pywe[1]) |> gpu, collect(pywe[2][:]) |> gpu, collect(pywe[3]') |> gpu, collect(pywe[4][:]) |> gpu, collect(pywe[5]') |> gpu, collect(pywe[6][:]) |> gpu) |> gpu
# end
#
# function loadPose2OdoNNModelIntoFluxGPU(dest::AbstractString)
#   weights = loadPyNNTxt(dest::AbstractString)
#   buildPose2OdoNN_01_FromWeightsGPU(weights)
# end

#


text:
module TestCommon

using Documenter
using LinearAlgebra
using IterativeSolvers
using Test
using Random

Random.seed!(1234321)

@testset "Basic operations" begin

@testset "Adivtype" begin
    A = rand(1:10, 5, 5)
    b = rand(Float32, 5)
    @test IterativeSolvers.Adivtype(A, b) == Float32
end

@testset "Identity preconditioner" begin
    P = Identity()
    x = rand(10)
    y = zero(x)

    # Should be a no-op
    @test P \ x == x
    @test ldiv!(P, copy(x)) == x
    @test ldiv!(y, P, copy(x)) == x
end

end

DocMeta.setdocmeta!(IterativeSolvers, :DocTestSetup, :(using IterativeSolvers); recursive=true)
doctest(IterativeSolvers, manual=true)

end # module


text:
# Collection of functions to calculate impact parameter and sky velocity

"""

Calculate impact parameter and sky velocity
"""
function calc_bv(s::State{T},i::Int64,j::Int64) where T<:AbstractFloat
    vsky = sqrt((s.v[1,j]-s.v[1,i])^2 + (s.v[2,j]-s.v[2,i])^2)
    bsky2 = (s.x[1,j]-s.x[1,i])^2 + (s.x[2,j]-s.x[2,i])^2
    return vsky,bsky2
end

"""

Calculate derivative of impact parameter and sky velocity wrt cartesian coordinates
"""
function calc_dbvdq!(s::State{T},dbvdq::Array{T},i::Int64,j::Int64) where T<:AbstractFloat
    vsky,bsky2 = calc_bv(s,i,j)
    gsky = g!(i,j,s.x,s.v)
    gdot = gd!(i,j,s.x,s.v,s.dqdt)
    fill!(dbvdq,zero(T))

    # partial derivative b and v_{sky} with respect to time:
    dvdt = ((s.v[1,j]-s.v[1,i])*(s.dqdt[(j-1)*7+4]-s.dqdt[(i-1)*7+4])+(s.v[2,j]-s.v[2,i])*(s.dqdt[(j-1)*7+5]-s.dqdt[(i-1)*7+5]))/vsky
    dbdt = 2.0 * gsky

    # Compute derivatives:
    indj = (j-1)*7+1
    indi = (i-1)*7+1
    for p=1:s.n
        indp = (p-1)*7
        for k=1:7
            dtdq = -((s.jac_step[indj  ,indp+k]-s.jac_step[indi  ,indp+k])*(s.v[1,j]-s.v[1,i])+(s.jac_step[indj+1,indp+k]-s.jac_step[indi+1,indp+k])*(s.v[2,j]-s.v[2,i])+
                          (s.jac_step[indj+3,indp+k]-s.jac_step[indi+3,indp+k])*(s.x[1,j]-s.x[1,i])+(s.jac_step[indj+4,indp+k]-s.jac_step[indi+4,indp+k])*(s.x[2,j]-s.x[2,i]))/gdot
            #v_{sky} = sqrt((v[1,j]-v[1,i])^2+(v[2,j]-v[2,i])^2)
            dbvdq[1,k,p] = ((s.jac_step[indj+3,indp+k]-s.jac_step[indi+3,indp+k])*(s.v[1,j]-s.v[1,i])+(s.jac_step[indj+4,indp+k]-s.jac_step[indi+4,indp+k])*(s.v[2,j]-s.v[2,i]))/vsky + dvdt*dtdq
            #b_{sky}^2 = (x[1,j]-x[1,i])^2+(x[2,j]-x[2,i])^2
            dbvdq[2,k,p] = 2*((s.jac_step[indj  ,indp+k]-s.jac_step[indi  ,indp+k])*(s.x[1,j]-s.x[1,i])+(s.jac_step[indj+1,indp+k]-s.jac_step[indi+1,indp+k])*(s.x[2,j]-s.x[2,i])) + dbdt*dtdq
        end
    end
    return vsky,bsky2
end

function calc_dbvdelements!(s::State{T},ts::TransitSnapshot{T}) where T <: AbstractFloat
    for ibv = 1:2, i=1:s.n, j=1:ts.nt
        # Now, multiply by the initial Jacobian to convert time derivatives to orbital elements:
        for k=1:s.n, l=1:7
            ts.dbvdelements[ibv,i,j,l,k] = zero(T)
            for p=1:s.n, q=1:7
                ts.dbvdelements[ibv,i,j,l,k] += ts.dbvdq0[ibv,i,j,q,p]*s.jac_init[(p-1)*7+q,(k-1)*7+l]
            end
        end
    end
end