text:
import GeoInterface

const drivers = AG.listdrivers()
const drivermapping = Dict(
    ".shp" => "ESRI Shapefile",
    ".gpkg" => "GPKG",
    ".geojson" => "GeoJSON",
)
# const fieldmapping = Dict(v => k for (k, v) in _FIELDTYPE)

# "return the corresponding `DataType` in julia"
# const _FIELDTYPE = Dict{AG.OGRFieldType, DataType}(
#     GDAL.OFTInteger         => Int32,
#     GDAL.OFTIntegerList     => Nothing,
#     GDAL.OFTReal            => Float64,
#     GDAL.OFTRealList        => Nothing,
#     GDAL.OFTString          => String,
#     GDAL.OFTStringList      => Nothing,
#     GDAL.OFTWideString      => Nothing, # deprecated
#     GDAL.OFTWideStringList  => Nothing, # deprecated
#     GDAL.OFTBinary          => Nothing,
#     GDAL.OFTDate            => Date,
#     GDAL.OFTTime            => Nothing,
#     GDAL.OFTDateTime        => DateTime,
#     GDAL.OFTInteger64       => Int64,
#     GDAL.OFTInteger64List   => Nothing)
# geomtypes = (IGeometry{ArchGDAL.gettype(ArchGDAL.getgeomdefn(featuredefn, i))} for i in 0:ngeom-1)
# field_types = (_FIELDTYPE[gettype(fielddefn)] for fielddefn in fielddefns)

# AG.createpoint(1.0, 2.0)
function get_coords(geom::ArchGDAL.IGeometry)
    GeoInterface.coordinates(geom)
end

function get_coords(geom::Array{ArchGDAL.IGeometry,1})
    GeoInterface.coordinates.(geom)
end

function get_coords(shp)
    coord = GeoInterface.coordinates.(shp[:geom])
    coord = hcat(coord...) |> transpose
    lon = coord[:, 1]
    lat = coord[:, 2]
    lon, lat
end

function readORG(fn::AbstractString, layer::Union{Integer,AbstractString}=0)
    ds = AG.read(fn)
    layer = AG.getlayer(ds, layer)
    table = AG.Table(layer)
    df = DataFrame(table)
    "" in names(df) && rename!(df, Dict(Symbol("") => :geom, ))  # needed for now
    lon, lat = get_coords(df)
    df.lon = lon
    df.lat = lat
    df
end

export readORG
export get_coords


text:
@enum BspSide Upper Lower

function parseside(c)
    if c == 'F' || c == 'L'
        return Lower
    end
    Upper
end

function parsesides(s)
    sides = Vector{BspSide}(undef, length(s))
    for (i, c) in enumerate(s)
        sides[i] = parseside(c)
    end
    sides
end

function bsprange(r, sides)
    a = collect(r)
    f, l = first(a), last(a)
    for s in sides
        if s == Upper
            f = f + ceil((l-f)/2)
        else
            l = l - ceil((l-f)/2)
        end
    end
    (f, l)
end

function seatpos(sides)
    rowsides = sides[1:7]
    colsides = sides[8:end]
    row = bsprange(0:127, rowsides)[1]
    col = bsprange(0:7, colsides)[1]
    (row, col)
end

function calcseatid(s)
    sides = parsesides(s)
    (row, col) = seatpos(sides)
    Int(row * 8 + col)
end

calcseadids(ss) = map(calcseatid, ss)

function readseats(file)
    s = read(file, String)
    split(s, "\n", keepempty=false)
end

function day05part1(ss)
    ids = calcseadids(ss)
    maximum(ids)
end

function day05part2(ss)
    ids = calcseadids(ss)
    maxid = 127*8*7
    filled = zeros(maxid)
    for id in ids
        filled[id+1] = 1
    end
    start = findfirst(x->x == 1, filled)
    idx = findfirst(x->x == 0, filled[start:end]) + start - 1
    idx-1
end

function day05(file)
    ss = readseats(file)
    p1 = day05part1(ss)
    p2 = day05part2(ss)
    (p1, p2)
end

#=

a = 0:127
@show Vector(a)

min(firstindex(a))

sides = parsesides("FBFBBFF")
bsprange(0:127, sides)

sides2 = parsesides("FBFBBFFRLR")
seatpos(sides2)

calcseatid("FBFBBFFRLR")
calcseatid("BFFFBBFRRR")
calcseatid("FFFBBBFRRR")
calcseatid("BBFFBBFRLL")

ss = readseats("day05/input.txt")
day05part1(ss)
day05part2(ss)

day05("day05/input.txt")
=#

ext: jl


text:
using ComputerVisionMetrics
using Documenter

DocMeta.setdocmeta!(ComputerVisionMetrics, :DocTestSetup, :(using ComputerVisionMetrics); recursive=true)

makedocs(;
    modules=[ComputerVisionMetrics],
    authors="Dale <djblack@uci.edu> and contributors",
    repo="https://github.com/Dale-Black/ComputerVisionMetrics.jl/blob/{commit}{path}#{line}",
    sitename="ComputerVisionMetrics.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Dale-Black.github.io/ComputerVisionMetrics.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Dale-Black/ComputerVisionMetrics.jl",
)


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct GeneralGovModel <: TurbineGov
        Rselect::Int
        fuel_flag::Int
        R::Float64
        Tpelec::Float64
        speed_error_signal::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        Kp_gov::Float64
        Ki_gov::Float64
        Kd_gov::Float64
        Td_gov::Float64
        valve_position_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        T_act::Float64
        K_turb::Float64
        Wf_nl::Float64
        Tb::Float64
        Tc::Float64
        T_eng::Float64
        Tf_load::Float64
        Kp_load::Float64
        Ki_load::Float64
        Ld_ref::Float64
        Dm::Float64
        R_open::Float64
        R_close::Float64
        Ki_mw::Float64
        A_set::Float64
        Ka::Float64
        Ta::Float64
        T_rate::Float64
        db::Float64
        Tsa::Float64
        Tsb::Float64
        R_lim::NamedTuple{(:up, :down), Tuple{Float64, Float64}}
        P_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        states_types::Vector{StateTypes}
        internal::InfrastructureSystemsInternal
    end

GE General Governor/Turbine Model. The GeneralGovModel (GGOV1) model is a general purpose governor model used for a variety of prime movers controlled by proportional-integral-derivative (PID) governors including gas turbines.

# Arguments
- `Rselect::Int`: Feedback signal for governor droop, validation range: `(-2, 1)`, action if invalid: `error`
- `fuel_flag::Int`: Flag Switch for fuel source characteristic, validation range: `(0, 1)`, action if invalid: `error`
- `R::Float64`: Speed droop parameter, validation range: `(eps(), nothing)`, action if invalid: `warn`
- `Tpelec::Float64`: Electrical power transducer time constant, seconds, validation range: `(eps(), nothing)`, action if invalid: `warn`
- `speed_error_signal::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Speed error signal limits
- `Kp_gov::Float64`: Governor proportional gain, validation range: `(0, nothing)`, action if invalid: `warn`
- `Ki_gov::Float64`: Governor integral gain, validation range: `(0, nothing)`, action if invalid: `warn`
- `Kd_gov::Float64`: Governor derivative gain, validation range: `(0, nothing)`, action if invalid: `warn`
- `Td_gov::Float64`: Governor derivative time constant, validation range: `(0, nothing)`, action if invalid: `warn`
- `valve_position_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Valve position limits
- `T_act::Float64`: Actuator time constant, validation range: `(0, nothing)`, action if invalid: `warn`
- `K_turb::Float64`: Turbine gain, validation range: `(0, nothing)`, action if invalid: `warn`
- `Wf_nl::Float64`: No load fuel flow, pu, validation range: `(0, nothing)`, action if invalid: `warn`
- `Tb::Float64`: Turbine lag time constant, sec, validation range: `(0, nothing)`, action if invalid: `warn`
- `Tc::Float64`: Turbine lead time constant, sec, validation range: `(0, nothing)`, action if invalid: `warn`
- `T_eng::Float64`: Transport lag time constant for diesel engine, sec, validation range: `(0, nothing)`, action if invalid: `warn`
- `Tf_load::Float64`: Load limiter time constant, validation range: `(0, nothing)`, action if invalid: `warn`
- `Kp_load::Float64`: Load limiter proportional gain for PI controller, validation range: `(0, nothing)`, action if invalid: `warn`
- `Ki_load::Float64`: Load integral gain for PI controller, validation range: `(0, nothing)`, action if invalid: `warn`
- `Ld_ref::Float64`: Load limiter integral gain for PI controller, validation range: `(0, nothing)`, action if invalid: `warn`
- `Dm::Float64`: Mechanical damping coefficient, pu, validation range: `(0, nothing)`, action if invalid: `warn`
- `R_open::Float64`: Maximum valve opening rate, pu/sec, validation range: `(0, nothing)`, action if invalid: `warn`
- `R_close::Float64`: Maximum valve closing rate, pu/sec, validation range: `(0, nothing)`, action if invalid: `warn`
- `Ki_mw::Float64`: Power controller (reset) gain, validation range: `(0, nothing)`, action if invalid: `warn`
- `A_set::Float64`: Acceleration limiter setpoint, pu/sec, validation range: `(0, nothing)`, action if invalid: `warn`
- `Ka::Float64`: Acceleration limiter gain, validation range: `(0, nothing)`, action if invalid: `warn`
- `Ta::Float64`: Acceleration limiter time constant , validation range: `(eps(), nothing)`, action if invalid: `error`
- `T_rate::Float64`: Turbine rating, validation range: `(0, nothing)`, action if invalid: `warn`
- `db::Float64`: Speed governor deadband, validation range: `(0, nothing)`, action if invalid: `warn`
- `Tsa::Float64`: Temperature detection lead time constant, validation range: `(0, nothing)`, action if invalid: `warn`
- `Tsb::Float64`: Temperature detection lag time constant, validation range: `(0, nothing)`, action if invalid: `warn`
- `R_lim::NamedTuple{(:up, :down), Tuple{Float64, Float64}}`: Maximum rate of load increa
- `P_ref::Float64`: Reference Power Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the GGOV1 model are:
	Pe: Machine Electrical Power Measurement,
	x_g1: Governor differential control,
	x_g2: Governor integral control,
	x_g3: Turbine actuator,
	x_g4: Turbine Lead-Lag,
	x_g5: Turbine load limiter measurement,
	x_g6: Turbine Load Limiter Integral Control,
	x_g7: Supervisory Load Control,
	x_g8: Acceleration Control,
	x_g9 Temperature Detection Lead - Lag:
- `n_states::Int`: GeneralGovModel has 10 states
- `states_types::Vector{StateTypes}`: GGOV1 has 10 differential states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct GeneralGovModel <: TurbineGov
    "Feedback signal for governor droop"
    Rselect::Int
    "Flag Switch for fuel source characteristic"
    fuel_flag::Int
    "Speed droop parameter"
    R::Float64
    "Electrical power transducer time constant, seconds"
    Tpelec::Float64
    "Speed error signal limits"
    speed_error_signal::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Governor proportional gain"
    Kp_gov::Float64
    "Governor integral gain"
    Ki_gov::Float64
    "Governor derivative gain"
    Kd_gov::Float64
    "Governor derivative time constant"
    Td_gov::Float64
    "Valve position limits"
    valve_position_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Actuator time constant"
    T_act::Float64
    "Turbine gain"
    K_turb::Float64
    "No load fuel flow, pu"
    Wf_nl::Float64
    "Turbine lag time constant, sec"
    Tb::Float64
    "Turbine lead time constant, sec"
    Tc::Float64
    "Transport lag time constant for diesel engine, sec"
    T_eng::Float64
    "Load limiter time constant"
    Tf_load::Float64
    "Load limiter proportional gain for PI controller"
    Kp_load::Float64
    "Load integral gain for PI controller"
    Ki_load::Float64
    "Load limiter integral gain for PI controller"
    Ld_ref::Float64
    "Mechanical damping coefficient, pu"
    Dm::Float64
    "Maximum valve opening rate, pu/sec"
    R_open::Float64
    "Maximum valve closing rate, pu/sec"
    R_close::Float64
    "Power controller (reset) gain"
    Ki_mw::Float64
    "Acceleration limiter setpoint, pu/sec"
    A_set::Float64
    "Acceleration limiter gain"
    Ka::Float64
    "Acceleration limiter time constant "
    Ta::Float64
    "Turbine rating"
    T_rate::Float64
    "Speed governor deadband"
    db::Float64
    "Temperature detection lead time constant"
    Tsa::Float64
    "Temperature detection lag time constant"
    Tsb::Float64
    "Maximum rate of load increa"
    R_lim::NamedTuple{(:up, :down), Tuple{Float64, Float64}}
    "Reference Power Set-point"
    P_ref::Float64
    ext::Dict{String, Any}
    "The states of the GGOV1 model are:
	Pe: Machine Electrical Power Measurement,
	x_g1: Governor differential control,
	x_g2: Governor integral control,
	x_g3: Turbine actuator,
	x_g4: Turbine Lead-Lag,
	x_g5: Turbine load limiter measurement,
	x_g6: Turbine Load Limiter Integral Control,
	x_g7: Supervisory Load Control,
	x_g8: Acceleration Control,
	x_g9 Temperature Detection Lead - Lag:"
    states::Vector{Symbol}
    "GeneralGovModel has 10 states"
    n_states::Int
    "GGOV1 has 10 differential states"
    states_types::Vector{StateTypes}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function GeneralGovModel(Rselect, fuel_flag, R, Tpelec, speed_error_signal, Kp_gov, Ki_gov, Kd_gov, Td_gov, valve_position_limits, T_act, K_turb, Wf_nl, Tb, Tc, T_eng, Tf_load, Kp_load, Ki_load, Ld_ref, Dm, R_open, R_close, Ki_mw, A_set, Ka, Ta, T_rate, db, Tsa, Tsb, R_lim, P_ref=1.0, ext=Dict{String, Any}(), )
    GeneralGovModel(Rselect, fuel_flag, R, Tpelec, speed_error_signal, Kp_gov, Ki_gov, Kd_gov, Td_gov, valve_position_limits, T_act, K_turb, Wf_nl, Tb, Tc, T_eng, Tf_load, Kp_load, Ki_load, Ld_ref, Dm, R_open, R_close, Ki_mw, A_set, Ka, Ta, T_rate, db, Tsa, Tsb, R_lim, P_ref, ext, [:Pe, :x_g1, :x_g2, :x_g3, :x_g4, :x_g5, :x_g6, :x_g7, :x_g8, :x_g9], 10, [StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Differential, StateTypes.Hybrid], InfrastructureSystemsInternal(), )
end

function GeneralGovModel(; Rselect, fuel_flag, R, Tpelec, speed_error_signal, Kp_gov, Ki_gov, Kd_gov, Td_gov, valve_position_limits, T_act, K_turb, Wf_nl, Tb, Tc, T_eng, Tf_load, Kp_load, Ki_load, Ld_ref, Dm, R_open, R_close, Ki_mw, A_set, Ka, Ta, T_rate, db, Tsa, Tsb, R_lim, P_ref=1.0, ext=Dict{String, Any}(), states=[:Pe, :x_g1, :x_g2, :x_g3, :x_g4, :x_g5, :x_g6, :x_g7, :x_g8, :x_g9], n_states=10, states_types=[StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Differential, StateTypes.Hybrid], internal=InfrastructureSystemsInternal(), )
    GeneralGovModel(Rselect, fuel_flag, R, Tpelec, speed_error_signal, Kp_gov, Ki_gov, Kd_gov, Td_gov, valve_position_limits, T_act, K_turb, Wf_nl, Tb, Tc, T_eng, Tf_load, Kp_load, Ki_load, Ld_ref, Dm, R_open, R_close, Ki_mw, A_set, Ka, Ta, T_rate, db, Tsa, Tsb, R_lim, P_ref, ext, states, n_states, states_types, internal, )
end

# Constructor for demo purposes; non-functional.
function GeneralGovModel(::Nothing)
    GeneralGovModel(;
        Rselect=1,
        fuel_flag=0,
        R=0,
        Tpelec=0,
        speed_error_signal=(min=0.0, max=0.0),
        Kp_gov=0,
        Ki_gov=0,
        Kd_gov=0,
        Td_gov=0,
        valve_position_limits=(min=0.0, max=0.0),
        T_act=0,
        K_turb=0,
        Wf_nl=0,
        Tb=0,
        Tc=0,
        T_eng=0,
        Tf_load=0,
        Kp_load=0,
        Ki_load=0,
        Ld_ref=0,
        Dm=0,
        R_open=0,
        R_close=0,
        Ki_mw=0,
        A_set=0,
        Ka=0,
        Ta=0,
        T_rate=0,
        db=0,
        Tsa=0,
        Tsb=0,
        R_lim=(up = 0.0, down = 0.0),
        P_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`GeneralGovModel`](@ref) `Rselect`."""
get_Rselect(value::GeneralGovModel) = value.Rselect
"""Get [`GeneralGovModel`](@ref) `fuel_flag`."""
get_fuel_flag(value::GeneralGovModel) = value.fuel_flag
"""Get [`GeneralGovModel`](@ref) `R`."""
get_R(value::GeneralGovModel) = value.R
"""Get [`GeneralGovModel`](@ref) `Tpelec`."""
get_Tpelec(value::GeneralGovModel) = value.Tpelec
"""Get [`GeneralGovModel`](@ref) `speed_error_signal`."""
get_speed_error_signal(value::GeneralGovModel) = value.speed_error_signal
"""Get [`GeneralGovModel`](@ref) `Kp_gov`."""
get_Kp_gov(value::GeneralGovModel) = value.Kp_gov
"""Get [`GeneralGovModel`](@ref) `Ki_gov`."""
get_Ki_gov(value::GeneralGovModel) = value.Ki_gov
"""Get [`GeneralGovModel`](@ref) `Kd_gov`."""
get_Kd_gov(value::GeneralGovModel) = value.Kd_gov
"""Get [`GeneralGovModel`](@ref) `Td_gov`."""
get_Td_gov(value::GeneralGovModel) = value.Td_gov
"""Get [`GeneralGovModel`](@ref) `valve_position_limits`."""
get_valve_position_limits(value::GeneralGovModel) = value.valve_position_limits
"""Get [`GeneralGovModel`](@ref) `T_act`."""
get_T_act(value::GeneralGovModel) = value.T_act
"""Get [`GeneralGovModel`](@ref) `K_turb`."""
get_K_turb(value::GeneralGovModel) = value.K_turb
"""Get [`GeneralGovModel`](@ref) `Wf_nl`."""
get_Wf_nl(value::GeneralGovModel) = value.Wf_nl
"""Get [`GeneralGovModel`](@ref) `Tb`."""
get_Tb(value::GeneralGovModel) = value.Tb
"""Get [`GeneralGovModel`](@ref) `Tc`."""
get_Tc(value::GeneralGovModel) = value.Tc
"""Get [`GeneralGovModel`](@ref) `T_eng`."""
get_T_eng(value::GeneralGovModel) = value.T_eng
"""Get [`GeneralGovModel`](@ref) `Tf_load`."""
get_Tf_load(value::GeneralGovModel) = value.Tf_load
"""Get [`GeneralGovModel`](@ref) `Kp_load`."""
get_Kp_load(value::GeneralGovModel) = value.Kp_load
"""Get [`GeneralGovModel`](@ref) `Ki_load`."""
get_Ki_load(value::GeneralGovModel) = value.Ki_load
"""Get [`GeneralGovModel`](@ref) `Ld_ref`."""
get_Ld_ref(value::GeneralGovModel) = value.Ld_ref
"""Get [`GeneralGovModel`](@ref) `Dm`."""
get_Dm(value::GeneralGovModel) = value.Dm
"""Get [`GeneralGovModel`](@ref) `R_open`."""
get_R_open(value::GeneralGovModel) = value.R_open
"""Get [`GeneralGovModel`](@ref) `R_close`."""
get_R_close(value::GeneralGovModel) = value.R_close
"""Get [`GeneralGovModel`](@ref) `Ki_mw`."""
get_Ki_mw(value::GeneralGovModel) = value.Ki_mw
"""Get [`GeneralGovModel`](@ref) `A_set`."""
get_A_set(value::GeneralGovModel) = value.A_set
"""Get [`GeneralGovModel`](@ref) `Ka`."""
get_Ka(value::GeneralGovModel) = value.Ka
"""Get [`GeneralGovModel`](@ref) `Ta`."""
get_Ta(value::GeneralGovModel) = value.Ta
"""Get [`GeneralGovModel`](@ref) `T_rate`."""
get_T_rate(value::GeneralGovModel) = value.T_rate
"""Get [`GeneralGovModel`](@ref) `db`."""
get_db(value::GeneralGovModel) = value.db
"""Get [`GeneralGovModel`](@ref) `Tsa`."""
get_Tsa(value::GeneralGovModel) = value.Tsa
"""Get [`GeneralGovModel`](@ref) `Tsb`."""
get_Tsb(value::GeneralGovModel) = value.Tsb
"""Get [`GeneralGovModel`](@ref) `R_lim`."""
get_R_lim(value::GeneralGovModel) = value.R_lim
"""Get [`GeneralGovModel`](@ref) `P_ref`."""
get_P_ref(value::GeneralGovModel) = value.P_ref
"""Get [`GeneralGovModel`](@ref) `ext`."""
get_ext(value::GeneralGovModel) = value.ext
"""Get [`GeneralGovModel`](@ref) `states`."""
get_states(value::GeneralGovModel) = value.states
"""Get [`GeneralGovModel`](@ref) `n_states`."""
get_n_states(value::GeneralGovModel) = value.n_states
"""Get [`GeneralGovModel`](@ref) `states_types`."""
get_states_types(value::GeneralGovModel) = value.states_types
"""Get [`GeneralGovModel`](@ref) `internal`."""
get_internal(value::GeneralGovModel) = value.internal

"""Set [`GeneralGovModel`](@ref) `Rselect`."""
set_Rselect!(value::GeneralGovModel, val) = value.Rselect = val
"""Set [`GeneralGovModel`](@ref) `fuel_flag`."""
set_fuel_flag!(value::GeneralGovModel, val) = value.fuel_flag = val
"""Set [`GeneralGovModel`](@ref) `R`."""
set_R!(value::GeneralGovModel, val) = value.R = val
"""Set [`GeneralGovModel`](@ref) `Tpelec`."""
set_Tpelec!(value::GeneralGovModel, val) = value.Tpelec = val
"""Set [`GeneralGovModel`](@ref) `speed_error_signal`."""
set_speed_error_signal!(value::GeneralGovModel, val) = value.speed_error_signal = val
"""Set [`GeneralGovModel`](@ref) `Kp_gov`."""
set_Kp_gov!(value::GeneralGovModel, val) = value.Kp_gov = val
"""Set [`GeneralGovModel`](@ref) `Ki_gov`."""
set_Ki_gov!(value::GeneralGovModel, val) = value.Ki_gov = val
"""Set [`GeneralGovModel`](@ref) `Kd_gov`."""
set_Kd_gov!(value::GeneralGovModel, val) = value.Kd_gov = val
"""Set [`GeneralGovModel`](@ref) `Td_gov`."""
set_Td_gov!(value::GeneralGovModel, val) = value.Td_gov = val
"""Set [`GeneralGovModel`](@ref) `valve_position_limits`."""
set_valve_position_limits!(value::GeneralGovModel, val) = value.valve_position_limits = val
"""Set [`GeneralGovModel`](@ref) `T_act`."""
set_T_act!(value::GeneralGovModel, val) = value.T_act = val
"""Set [`GeneralGovModel`](@ref) `K_turb`."""
set_K_turb!(value::GeneralGovModel, val) = value.K_turb = val
"""Set [`GeneralGovModel`](@ref) `Wf_nl`."""
set_Wf_nl!(value::GeneralGovModel, val) = value.Wf_nl = val
"""Set [`GeneralGovModel`](@ref) `Tb`."""
set_Tb!(value::GeneralGovModel, val) = value.Tb = val
"""Set [`GeneralGovModel`](@ref) `Tc`."""
set_Tc!(value::GeneralGovModel, val) = value.Tc = val
"""Set [`GeneralGovModel`](@ref) `T_eng`."""
set_T_eng!(value::GeneralGovModel, val) = value.T_eng = val
"""Set [`GeneralGovModel`](@ref) `Tf_load`."""
set_Tf_load!(value::GeneralGovModel, val) = value.Tf_load = val
"""Set [`GeneralGovModel`](@ref) `Kp_load`."""
set_Kp_load!(value::GeneralGovModel, val) = value.Kp_load = val
"""Set [`GeneralGovModel`](@ref) `Ki_load`."""
set_Ki_load!(value::GeneralGovModel, val) = value.Ki_load = val
"""Set [`GeneralGovModel`](@ref) `Ld_ref`."""
set_Ld_ref!(value::GeneralGovModel, val) = value.Ld_ref = val
"""Set [`GeneralGovModel`](@ref) `Dm`."""
set_Dm!(value::GeneralGovModel, val) = value.Dm = val
"""Set [`GeneralGovModel`](@ref) `R_open`."""
set_R_open!(value::GeneralGovModel, val) = value.R_open = val
"""Set [`GeneralGovModel`](@ref) `R_close`."""
set_R_close!(value::GeneralGovModel, val) = value.R_close = val
"""Set [`GeneralGovModel`](@ref) `Ki_mw`."""
set_Ki_mw!(value::GeneralGovModel, val) = value.Ki_mw = val
"""Set [`GeneralGovModel`](@ref) `A_set`."""
set_A_set!(value::GeneralGovModel, val) = value.A_set = val
"""Set [`GeneralGovModel`](@ref) `Ka`."""
set_Ka!(value::GeneralGovModel, val) = value.Ka = val
"""Set [`GeneralGovModel`](@ref) `Ta`."""
set_Ta!(value::GeneralGovModel, val) = value.Ta = val
"""Set [`GeneralGovModel`](@ref) `T_rate`."""
set_T_rate!(value::GeneralGovModel, val) = value.T_rate = val
"""Set [`GeneralGovModel`](@ref) `db`."""
set_db!(value::GeneralGovModel, val) = value.db = val
"""Set [`GeneralGovModel`](@ref) `Tsa`."""
set_Tsa!(value::GeneralGovModel, val) = value.Tsa = val
"""Set [`GeneralGovModel`](@ref) `Tsb`."""
set_Tsb!(value::GeneralGovModel, val) = value.Tsb = val
"""Set [`GeneralGovModel`](@ref) `R_lim`."""
set_R_lim!(value::GeneralGovModel, val) = value.R_lim = val
"""Set [`GeneralGovModel`](@ref) `P_ref`."""
set_P_ref!(value::GeneralGovModel, val) = value.P_ref = val
"""Set [`GeneralGovModel`](@ref) `ext`."""
set_ext!(value::GeneralGovModel, val) = value.ext = val
"""Set [`GeneralGovModel`](@ref) `states_types`."""
set_states_types!(value::GeneralGovModel, val) = value.states_types = val


text:
using TermInterface
using Test

@testset "Expr" begin
    ex = :(f(a,b))
    @test operation(ex) == :f
    @test arguments(ex) == [:a, :b]
    @test exprhead(ex) == :call
    @test ex == similarterm(ex, :f, [:a, :b])

    ex = :(arr[i, j])
    @test operation(ex) == getindex
    @test arguments(ex) == [:arr, :i, :j]
    @test exprhead(ex) == :ref
    @test ex == similarterm(Expr, :ref, [:arr, :i, :j]; exprhead=:ref)
    @test ex == similarterm(ex, :ref, [:arr, :i, :j])
end

# TODO add SymbolicUtils tests??


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
using Test

DDEProblemLibrary.importddeproblems()

@testset "init" begin
    prob = DDEProblemLibrary.prob_dde_constant_1delay_ip
    prob_scalar = DDEProblemLibrary.prob_dde_constant_1delay_scalar

    inferred = [BS3(), Tsit5(), RK4(), Vern6()]
    for alg in inferred
        ddealg = MethodOfSteps(alg)

        @inferred init(prob, ddealg)
        @inferred init(prob_scalar, ddealg)
    end

    notinferred = [SDIRK2(), TRBDF2(), KenCarp4(), Rosenbrock23(), Rodas4()]
    for alg in notinferred
        ddealg = MethodOfSteps(alg)

        @test_broken @inferred init(prob, ddealg)
        @test_broken @inferred init(prob_scalar, ddealg)
    end
end


text:
using Statistics
using PyPlot, CSV
using DataFrames, DataFramesMeta
import FileIO:load, save
# ---------
# Experimental Prediction 1:
# Process dataframe files (outputs of simulations)
# in order to create the csv data for the plots in Figure 10
#
# Sgm stands for Bayes Factor Surprise
# Ssh stands for Shannon Surprise
# ---------
function processAllLearners_first(;
                        filepath = "/path/to/experiment_jldfiles",
                            )

    jldfiles = searchdirendings(filepath, ".jld2")
    for i in 1:length(jldfiles)
        file = load(joinpath(filepath, jldfiles[i]))
        df = file["df"]
        plotAllExperimentsettings_first(df)
        csvAllExperimentsettings_first(df)
    end
end
# Make plots within julia: not used eventually for the paper, only for inspection
function plotAllExperimentsettings_first(df;
                                savepath = homedir()
                                )

    learnertype = unique(df[:learnertype])
    learnerparam = unique(df[:learnerparam])
    if length(learnertype) > 1
        error("There are more than one learners!")
    else
        learnername = learnertype[1]
    end
    savepathlearner = joinpath(savepath, string(learnername) * string(formatparamlist(learnerparam)) * "_plots")
    mkdir(savepathlearner)

    δs = df[:δs][1]
    dfmergedseeds = df_mergeseeds_Surprise_first(df)
    for σ in unique(df[:sigma])
        for p_c in unique(df[:changeprobability])
            for nstep in unique(df[:nsteps])
                for y_hat_center in unique(df[:y_hat_center])
                    for var_center in unique(df[:var_center])
                        for Δy_hat in unique(df[:Δy_hat])
                            for Δvar in  unique(df[:Δvar])
                                for Δδ in unique(df[:Δδ])
                                    dftrimmedmerged = df_experiment_trim_first(dfmergedseeds,
                                                        svalue = σ,
                                                        changeprobabilityvalue = p_c,
                                                        nstepvalue = nstep,
                                                        y_hat_centervalue = y_hat_center,
                                                        var_centervalue = var_center,
                                                        Δy_hatvalue = Δy_hat,
                                                        Δvarvalue = Δvar,
                                                        Δδvalue = Δδ)
                                    titlestring =  "σ: $(σ), pc: $(p_c), steps: $(nstep), y_hat_center: $(y_hat_center),
                var_center: $(var_center), Δy: $(Δy_hat), Δvar: $(Δvar), Δδ: $(Δδ) "
                                    figurestring = "σ$(σ)_" *
                                                    "pc$(p_c)_" *
                                                    "steps$(nstep)_" *
                                                    "y_hat_center$(y_hat_center)_" *
                                                    "var_center$(var_center)_" *
                                                    "Δy$(Δy_hat)_" *
                                                    "Δvar$(Δvar)_" *
                                                    "Δδ$(Δδ)"
                                        meanSgmpositive, stdSgmpositive, seSgmpositive = calc_stats_nan(dftrimmedmerged[:allmean_Sgm_positive_s][1])
                                        meanSgmnegative, stdSgmnegative, seSgmnegative = calc_stats_nan(dftrimmedmerged[:allmean_Sgm_negative_s][1])
                                        plotmeanSurpriseOverδ(meanSgmpositive,
                                                    meanSgmnegative,
                                                    seSgmpositive,
                                                    seSgmpositive,
                                                    titlestring = titlestring,
                                                    ylabelstring = "Mean Bayes Factor Surprise",
                                                    xdata = δs
                                                    )
                                        PyPlot.savefig(joinpath(savepathlearner,
                                                    "exp_Sgm"* string(learnername) *
                                                        figurestring * ".png"))
                                        PyPlot.close()
                                        plt.ioff()
                                        meanSshpositive, stdSshpositive, seSshpositive = calc_stats_nan(dftrimmedmerged[:allmean_Ssh_positive_s][1])
                                        meanSshnegative, stdSshnegative, seSshnegative = calc_stats_nan(dftrimmedmerged[:allmean_Ssh_negative_s][1])
                                        plotmeanSurpriseOverδ(meanSshpositive,
                                                    meanSshnegative,
                                                    seSshpositive,
                                                    seSshnegative,
                                                    titlestring = titlestring,
                                                    ylabelstring = "Mean Shannon Surprise",
                                                    xdata = δs,
                                                    ylims = [0., 1.]
                                                    )
                                        PyPlot.savefig(joinpath(savepathlearner,
                                                    "exp_Ssh"* string(learnername) *
                                                        figurestring * ".png"))
                                        PyPlot.close()
                                        plt.ioff()
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end
function plotmeanSurpriseOverδ(mean_Surprise_positive_s,
                            mean_Surprise_negative_s,
                            se_Surprise_positive_s,
                            se_Surprise_negative_s;
                            xlabelstring = "δ",
                            ylabelstring = "Mean S",
                            titlestring = "",
                            ylims = [0., 3.],
                            xdata = []
                            )
    fig = figure(); ax = gca()
    plot_shadedlines([mean_Surprise_positive_s, mean_Surprise_negative_s],
                    [se_Surprise_positive_s, se_Surprise_negative_s],
                    ["pos", "neg"], ax,
                    xlabelstring = xlabelstring,
                    titlestring = titlestring,
                    ylabelstring = ylabelstring,
                    ylims = ylims,
                    xdata = xdata
                    )

end
# Create the csv data for the plots for Experimental Prediction 1
function csvAllExperimentsettings_first(df;
                                savepath = homedir()
                                )

    learnertype = unique(df[:learnertype])
    learnerparam = unique(df[:learnerparam])
    if length(learnertype) > 1
        error("There are more than one learners!")
    else
        learnername = learnertype[1]
    end
    savepathlearner = joinpath(savepath, string(learnername) * string(formatparamlist(learnerparam)) * "_csv")
    mkdir(savepathlearner)

    δs = df[:δs][1]
    dfmergedseeds = df_mergeseeds_Surprise_first(df)
    for σ in unique(df[:sigma])
        for p_c in unique(df[:changeprobability])
            for nstep in unique(df[:nsteps])
                for y_hat_center in unique(df[:y_hat_center])
                    for var_center in unique(df[:var_center])
                        for Δy_hat in unique(df[:Δy_hat])
                            for Δvar in  unique(df[:Δvar])
                                for Δδ in unique(df[:Δδ])

                                    dfresults = DataFrame()
                                    dfresults[:x] = δs

                                    dftrimmedmerged = df_experiment_trim_first(dfmergedseeds,
                                                        svalue = σ,
                                                        changeprobabilityvalue = p_c,
                                                        nstepvalue = nstep,
                                                        y_hat_centervalue = y_hat_center,
                                                        var_centervalue = var_center,
                                                        Δy_hatvalue = Δy_hat,
                                                        Δvarvalue = Δvar,
                                                        Δδvalue = Δδ)
                                    csvfilestring = "sigma$(σ)_" *
                                                    "pc$(p_c)_" *
                                                    "steps$(nstep)_" *
                                                    "y_hat_center$(y_hat_center)_" *
                                                    "var_center$(var_center)_" *
                                                    "Dy$(Δy_hat)_" *
                                                    "Dvar$(Δvar)_" *
                                                    "Dd$(Δδ)"
                                    meanSgmpositive, stdSgmpositive, seSgmpositive = calc_stats_nan(dftrimmedmerged[:allmean_Sgm_positive_s][1])
                                    meanSgmnegative, stdSgmnegative, seSgmnegative = calc_stats_nan(dftrimmedmerged[:allmean_Sgm_negative_s][1])

                                    dfresults[Symbol(string(learnername)*"_Sgmpositive_mean")] = meanSgmpositive
                                    dfresults[Symbol(string(learnername)*"_Sgmpositive_std")] = stdSgmpositive
                                    dfresults[Symbol(string(learnername)*"_Sgmpositive_se")] = seSgmpositive

                                    dfresults[Symbol(string(learnername)*"_Sgmnegative_mean")] = meanSgmnegative
                                    dfresults[Symbol(string(learnername)*"_Sgmnegative_std")] = stdSgmnegative
                                    dfresults[Symbol(string(learnername)*"_Sgmnegative_se")] = seSgmnegative

                                    meanSshpositive, stdSshpositive, seSshpositive = calc_stats_nan(dftrimmedmerged[:allmean_Ssh_positive_s][1])
                                    meanSshnegative, stdSshnegative, seSshnegative = calc_stats_nan(dftrimmedmerged[:allmean_Ssh_negative_s][1])

                                    dfresults[Symbol(string(learnername)*"_Sshpositive_mean")] = meanSshpositive
                                    dfresults[Symbol(string(learnername)*"_Sshpositive_std")] = stdSshpositive
                                    dfresults[Symbol(string(learnername)*"_Sshpositive_se")] = seSshpositive

                                    dfresults[Symbol(string(learnername)*"_Sshnegative_mean")] = meanSshnegative
                                    dfresults[Symbol(string(learnername)*"_Sshnegative_std")] = stdSshnegative
                                    dfresults[Symbol(string(learnername)*"_Sshnegative_se")] = seSshnegative

                                    CSV.write(joinpath(savepathlearner, "exp_"*
                                                string(learnername) *
                                                csvfilestring *".csv"),
                                                dfresults, delim = " ")
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

# -----------------
# Handle dataframes
# -----------------
function df_mergeseeds_Surprise_first(df)
    df = by(df, [:learnertype, :learnerparam, :sigma, :changeprobability, :nsteps, :y_hat_center, :var_center, :Δy_hat, :Δvar, :Δδ],
                allmean_Sgm_positive_s = :mean_Sgm_positive_s => x -> [collect(reduce(hcat, x))],
                allmean_Sgm_negative_s = :mean_Sgm_negative_s => x -> [collect(reduce(hcat, x))],
                allse_Sgm_positive_s = :se_Sgm_positive_s => x -> [collect(reduce(hcat, x))],
                allse_Sgm_negative_s = :se_Sgm_negative_s => x -> [collect(reduce(hcat, x))],
                allmean_Ssh_positive_s = :mean_Ssh_positive_s => x -> [collect(reduce(hcat, x))],
                allmean_Ssh_negative_s = :mean_Ssh_negative_s => x -> [collect(reduce(hcat, x))],
                allse_Ssh_positive_s = :se_Ssh_positive_s => x -> [collect(reduce(hcat, x))],
                allse_Ssh_negative_s = :se_Ssh_negative_s => x -> [collect(reduce(hcat, x))]
                )
end
function df_experiment_trim_first(df; svalue = 0.5, changeprobabilityvalue = 0.05,
                    nstepvalue = 1000,
                    y_hat_centervalue = 1.,
                    var_centervalue = 1.,
                    Δy_hatvalue = 1.,
                    Δvarvalue = 1.,
                    Δδvalue = 0.1
                    )
    dftrimmed = @where(df, :changeprobability .== changeprobabilityvalue,
                            :sigma .== svalue,
                            :nsteps .== nstepvalue,
                            :y_hat_center .== y_hat_centervalue,
                            :var_center .== var_centervalue,
                            :Δy_hat .== Δy_hatvalue,
                            :Δvar .== Δvarvalue,
                            :Δδ .== Δδvalue)
end
function getstatsacrossseeds(dftrimmedmerged;
                    columnname = :allmean_Sgm_positive_s)
    dataacrossseeds = vcat(dftrimmedmerged[columnname]...)
    m, s, serror = calc_stats(dataacrossseeds)
end


text:
# This file illustrates calculation of power spectral density and total average power of an white Gaussian noise.

using FFTW
using Plots
using DSP
using Statistics

# Construct noise signal
N0 = 50.
fs = 100
ts = 1 / fs
ln = 2^10
tn = (0 : ln - 1) * ts
n = sqrt(N0 / 2 * fs) * randn(ln)

# Compute autocorrelation function.
Rn = xcorr(n, n) / (ln - 1)
τ = (-ln + 1 : ln - 1) * ts
lr = length(Rn)

# Fourier tranform of autocorrelation function
xrf = fft(Rn) * ts
psdr = abs.(xrf)
ffr = (0 : 1 / (lr - 1) : 1) * fs
Δfr = fs / (lr - 1)
Pfr = sum(psdr[1 : end - 1]) * Δfr

# Compute power spectral density
xfn = fft(n) * ts
psdn = abs.(xfn).^2 / (ln - 1) / ts
ffn = (0 : 1 / (ln - 1) : 1) * fs

# Compute powers
Pt = sum(n[1 : ln - 1].^2) / (ln - 1)
Δf = fs / (ln - 1)
Pfn = sum(psdn[1 : end - 1]) * Δf

# Compare all the powers and variance
@show  Pfr, Pt, Pfn, var(n)

# Plots
plt = plot(layout=4)
plot!(tn, n, subplot=1)
plot!(τ, Rn, subplot=2)
plot!(ffn, psdn, subplot=3)
plot!(ffr, psdr, subplot=4)
display(plt)


text:
# Define the dynamics model of the model.
struct DoubleIntegratorGame{N,M,P,SVu,SVx,SVz} <: AbstractGameModel
    n::Int  # Number of states
    m::Int  # Number of controls
	p::Int  # Number of players
	ni::Vector{Int}  # Number of states for each player
	mi::Vector{Int}  # Number of controls for each player
	pu::SVu # Indices of the each player's controls
	px::SVx # Indices of the each player's x and y positions
	pz::SVz # Indices of the each player's states
end

function DoubleIntegratorGame(;p::Int=2, d::Int=2)
	# p = number of players
	# d = dimension of the integrator
	n = 2d*p
	m = d*p
	pu = [SVector{d,Int}([i + (j-1)*p for j=1:d]) for i=1:p]
	px = [SVector{2,Int}([i + (j-1)*p for j=1:2]) for i=1:p]
	pz = [SVector{2d,Int}([i + (j-1)*p for j=1:2d]) for i=1:p]
	TYPE = typeof.((pu,px,pz))
	ni = 2d*ones(Int,p)
	mi = d*ones(Int,p)
	return DoubleIntegratorGame{n,m,p,TYPE...}(n,m,p,ni,mi,pu,px,pz)
end

@generated function RobotDynamics.dynamics(model::DoubleIntegratorGame{N,M,P}, x, u) where {N,M,P}
	qd  = [:(x[$i]) for i=M+1:N]
	qdd = [:(u[$i]) for i=1:M]
	return :(SVector{$N}($(qd...), $(qdd...)))
end

dim(model::DoubleIntegratorGame) = model.mi[1]


text:
# 6D Vlasov - Poisson Single Species
# For mpi use
using MPI

include("pif_tools.jl")
include("particle_sampling.jl")
importall pif

using psmp
using ProgressMeter

# Initialize MPI if not already
if (~MPI.Initialized())
  MPI.Init()
end
mpi_comm = MPI.COMM_WORLD # Communicator
mpi_rk = MPI.Comm_rank(mpi_comm) # Rank
mpi_sz = MPI.Comm_size(mpi_comm) # Size
MPI.Barrier(mpi_comm)

async_mpi=true; # prefer Non-blocking MPI communication

Np=Int(1e4) #Number of particles
Ns=1 #Number of species
dt=0.05
tmax=5

Nx1=1
Nx2=0
Nx3=0

mpi_Np=Int(Np/mpi_sz)
print(mpi_Np, " particles on proc", mpi_rk,"\n")

# Draw uniform
x1n,x2n,x3n,v1n,v2n,v3n=psmp.uniform_6d(mpi_Np,20 + mpi_rk*mpi_Np)
# sample testcase and retrieve default parameters
params=psmp.weibel(x1n,x2n,x3n,v1n,v2n,v3n)
# params=psmp.weibel_streaming(x1n,x2n,x3n,v1n,v2n,v3n)
#params=psmp.landau(x1n,x2n,x3n,v1n,v2n,v3n,epsilon=0.5)
# params=psmp.jeans(x1n,x2n,x3n,v1n,v2n,v3n,epsilon=0.05)

dt=params["dt"]
tmax=params["tmax"]
c=params["c"]
L=params["L"]
qn=ones(Float64,mpi_Np)*params["q"]
mn=ones(Float64,mpi_Np)*params["m"]

dt=0.1; tmax=400;
# Likelihoods
fn=params["f0"](x1n,x2n,x3n,v1n,v2n,v3n)
gn=params["g0"](x1n,x2n,x3n,v1n,v2n,v3n)

import OrdinaryDiffEq: McAte2ConstantCache, McAte8ConstantCache
# Symmetric composition
dt_SS=McAte2ConstantCache(Float64,Float64) #
SS_len=Int(length(fieldnames(dt_SS))/2)
dtA=Array{Float64}(SS_len)
dtB=similar(dtA)
for n=1:SS_len
  eval(parse("dtA[$n]=dt_SS.a$n; dtB[$n]=dt_SS.b$n;"))
end
# dtA=0.5;dtB=0.5; # Strang splitting

#Particle mesh coupling
pm=pm_pif3d(Nx1,Nx2,Nx3,L...)
print("Total Fourier modes:", pm.N1*pm.N2*pm.N3, "\n")

# Fields
rho=Array{Complex{Float64}}(pm.N1,pm.N2,pm.N3)
Phi=Array{Complex{Float64}}(pm.N1,pm.N2,pm.N3)
E1=similar(Phi);E2=similar(Phi);E3=similar(Phi)
B1=similar(Phi);B2=similar(Phi);B3=similar(Phi)

# Weights
wn=fn./gn

Nt=Int(ceil(tmax/dt))
Epot=zeros(Float64,Nt,3)
Ekin=zeros(Float64,Nt,3)
Momentum=zeros(Float64,Nt,3)
kineticenergy=zeros(Float64,Nt)
ttime=collect((0:Nt-1)*dt)

#Initial poisson solve
rho=accum_osde(pm,x1n,x2n,x3n, wn.*qn)./Np
rho=MPI.allreduce(rho, MPI.SUM, mpi_comm)
solve_poisson!(pm,rho,Phi)
E1,E2,E3=gradient(pm,-Phi)
print(gauss_error(pm,rho,E1,E2,E3),"\n" )

# Initialize magnetic fields from initial condition
if haskey(params,"B")
   B1.=to_grid(pm,params["B"][1])
   B2.=to_grid(pm,params["B"][2])
   B3.=to_grid(pm,params["B"][3])
end

DTYPE=Float64
#Diagnostics
Nt=Int(ceil(tmax/dt))
Epot=zeros(DTYPE,Nt,3)
Bpot=zeros(DTYPE,Nt,3)
Ekin=zeros(DTYPE,Nt,3)
Momentum=zeros(DTYPE,Nt,3)
kineticenergy=zeros(DTYPE,Nt)
ttime=collect((0:Nt-1)*dt)

# dt_SS=McAte8ConstantCache(Float64,Float64) #
# SS_len=Int(length(fieldnames(dt_SS))/2)
# dtA2=Array{Float64}(SS_len)
# dtB2=similar(dtA2)
# for n=1:SS_len
#   eval(parse("dtA2[$n]=dt_SS.a$n; dtB2[$n]=dt_SS.b$n;"))
# end
using ProgressMeter

# Allocate temporary arrays
J1=zeros(Complex{Float64},pm.N1,pm.N2,pm.N3 )
J2=zeros(Complex{Float64},pm.N1,pm.N2,pm.N3 )
J3=zeros(Complex{Float64},pm.N1,pm.N2,pm.N3 )

## Prepare plot
if mpi_rk==0
  using GLVisualize, GeometryTypes, Colors,GLAbstraction
  windowp = glscreen(color = RGBA(0f0, 0f0, 0f0, 0f0));
  @async renderloop(windowp)

  # 3d view into the particle coordinates
  pview3d(x1,x2,x3,v1,v2,v3)=(Float32.(x1),Float32.(v1),Float32.(v2))

# GeometryTypes.scalematrix(Vec3f0(5f0))

  particle_size=0.03f0
  pshape=GeometryTypes.Sphere{Float32}(Vec3f0(0), particle_size)

  #Scale the particle shape ot a siz
  pshape=GLAbstraction.scalematrix(Vec3f0(particle_size)*2 )*
                  GLNormalMesh(loadasset("cat.obj"))

  # _view(particle_vis, camera = :perspective)
  color = fill(RGBA{Float32}(1., 1., 1., 0.8),Np)

  hull=AABB(Vec3f0(minimum.(pview3d(x1n,x2n,x3n,v1n,v2n,v3n) )...),
          Vec3f0(minimum.(pview3d(x1n,x2n,x3n,v1n,v2n,v3n) )...))

  glparticles = visualize( (   pshape,
               GeometryTypes.Point3f0.(pview3d(x1n,x2n,x3n,v1n,v2n,v3n)... )),
               hull =hull,
              color=color, billboard = true).children[]
_view(glparticles,windowp,camera = :perspective)

# Update particles
GLAbstraction.set_arg!(glparticles, :position, GeometryTypes.Point3f0.(
             pview3d(x1n,x2n,x3n,v1n,v2n,v3n)... ) )

end

MPI.Barrier(mpi_comm)
if mpi_rk==0
  pg = Progress(Nt, 1,"Vlasov-Maxwell 6D")
end
tic()
for tdx=1:Nt
  if (tdx==3)
    tic()
  end
  if (mpi_rk==0)
    ProgressMeter.next!(pg)
  end

  Bpot[tdx,1]=L2norm(pm,B1)/2.
  Bpot[tdx,2]=L2norm(pm,B2)/2.
  Bpot[tdx,3]=L2norm(pm,B3)/2.
  Epot[tdx,1]=L2norm(pm,E1)/2.
  Epot[tdx,2]=L2norm(pm,E2)/2.
  Epot[tdx,3]=L2norm(pm,E3)/2.
  Ekin[tdx,1]=sum(v1n.^2.*wn)/Np/2.
  Ekin[tdx,2]=sum(v2n.^2.*wn)/Np/2.
  Ekin[tdx,3]=sum(v3n.^2.*wn)/Np/2.
  Momentum[tdx,1]=sum(v1n.*wn)/Np + (dot(pm,E2,B3)-dot(pm,E3,B2))
  Momentum[tdx,2]=sum(v2n.*wn)/Np + (dot(pm,E3,B1)-dot(pm,E1,B3))
  Momentum[tdx,3]=sum(v3n.*wn)/Np + (dot(pm,E1,B2)-dot(pm,E2,B1))

  if (mpi_rk==0)
    Momentum[tdx,1]+= (dot(pm,E2,B3)-dot(pm,E3,B2))
    Momentum[tdx,2]+= (dot(pm,E3,B1)-dot(pm,E1,B3))
    Momentum[tdx,3]+= (dot(pm,E1,B2)-dot(pm,E2,B1))
  end

 # integrate_H_E(pm,E1,E2,E3,B1,B2,B3,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,dt)
 # integrate_H_B(pm, E1,E2,E3,B1,B2,B3,c,dt)
 #  J1[:]=0.0;J2[:]=0.0;J3[:]=0.0 #mandatory
 # # Hamiltonian splitting
 # integrate_Hp(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
 #                   J1,J2,J3,B1,B2,B3,dt)
 #    J1[pm.N1n,pm.N2n,pm.N3n].=0
 #   J2[pm.N1n,pm.N2n,pm.N3n].=0
 #   J3[pm.N1n,pm.N2n,pm.N3n].=0
 #  E1.-=J1/Np;E2.-=J2/Np; E3.-=J3/Np
  J_req=Array{MPI.Request}(3)
  gdx=1;
  for gdx=1:length(dtA)

    if gdx==1
    integrate_H_E(pm,E1,E2,E3,B1,B2,B3,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,dtA[gdx]*dt)
    integrate_H_B(pm, E1,E2,E3,B1,B2,B3,c,dtA[gdx]*dt)
    end

   #J1_=similar(J1);J2_=similar(J2) ;J3_=similar(J3)
   J1[:]=0.0;J2[:]=0.0;J3[:]=0.0 #mandatory
   # Hamiltonian splitting
   # integrate_Hp_split_sym(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
   #                 J1,J2,J3,B1,B2,B3,[dtA[gdx]*dt],[dtB[gdx]*dt])

    # Exponential boris - not symplectic, but good for strong field
     integrate_Hp_boris_exp(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
                       J1,J2,J3,B1,B2,B3,dtA[gdx]*dt,dtB[gdx]*dt,async_mpi)

     # integrate_Hp_midpoint(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
     #                  J1,J2,J3,B1,B2,B3,dt2)

    # Subtract steady background for momentum conservation
    # J1[pm.N1n,:,:].=0
    # J2[:,pm.N2n,:].=0
    # J3[:,:,pm.N3n].=0
     #
     J1[pm.N1n,pm.N2n,pm.N3n].=0
     J2[pm.N1n,pm.N2n,pm.N3n].=0
     J3[pm.N1n,pm.N2n,pm.N3n].=0

     if async_mpi
      J_req[1],J1=MPI.iallreduce(J1, MPI.SUM, mpi_comm)
      J_req[2],J2=MPI.iallreduce(J2, MPI.SUM, mpi_comm)
      J_req[3],J3=MPI.iallreduce(J3, MPI.SUM, mpi_comm)

      #Only for boris, finish last step before communication
      integrate_vxB(pm, B1,B2,B3,x1n,x2n,x3n,v1n,v2n,v3n,
                                   qn, mn, dtB[gdx]*dt)
      else
        J1=MPI.allreduce(J1, MPI.SUM, mpi_comm)
        J2=MPI.allreduce(J2, MPI.SUM, mpi_comm)
        J3=MPI.allreduce(J3, MPI.SUM, mpi_comm)
    end

    if gdx==length(dtB)
      dt2=dtB[gdx]*dt
    else
      dt2=(dtA[gdx+1]+dtB[gdx])*dt
    end

      if async_mpi
      integrate_H_B(pm, E1,zeros(E2),zeros(E3),B1,B2,B3,c,dt2)
      MPI.Wait!(J_req[1])
      E1.-=J1/Np;
      integrate_H_E(pm,E1,x1n,x2n,x3n,v1n,qn,mn,dt2)

      integrate_H_B(pm, zeros(E1),E2,zeros(E3),B1,B2,B3,c,dt2)
      MPI.Wait!(J_req[2])
      E2.-=J2/Np;
      integrate_H_E(pm,E2,x1n,x2n,x3n,v2n,qn,mn,dt2)

      integrate_H_B(pm, zeros(E1),zeros(E2),E3,B1,B2,B3,c,dt2)
      MPI.Wait!(J_req[3])
      E3.-=J3/Np
      integrate_H_E(pm,E3,x1n,x2n,x3n,v3n,qn,mn,dt2)

      integrate_Faraday(pm,E1,E2,E3,B1,B2,B3,dt2)

    else
      integrate_H_B(pm, E1,E2,E3,B1,B2,B3,c,dt2)
      integrate_H_E(pm,E1,E2,E3,B1,B2,B3,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,dt2)
    end
  end

  x1n=mod.(x1n,pm.L1);  x2n=mod.(x2n,pm.L2);  x3n=mod.(x3n,pm.L3)

  if (mpi_rk==0)
    GLAbstraction.set_arg!(glparticles, :position, GeometryTypes.Point3f0.(
                 pview3d(x1n,x2n,x3n,v1n,v2n,v3n)... ) )
  end

end
toc()

Ekin=MPI.allreduce(Ekin, MPI.SUM, mpi_comm)
Momentum=MPI.allreduce(Momentum, MPI.SUM, mpi_comm)

# end
#
# @profile main_loop(x1n,x2n,x3n,v1n,v2n,v3n,mn,qn,wn,
#                    Ekin,Epot,Bpot,Momentum,
#                    pm,dt,dtA,dtB,Nt,c,
#                    mpi_rk,mpi_comm, async_mpi)
# using ProfileView
# ProfileView.view()

# Profile.print(format=:tree,sortedby=:count)
# Profile.print(format=:flat,sortedby=:count)
# T=Float64;pdx=1;
#
# rho=accum_osde(pm,x1n,x2n,x3n, wn.*qn)./Np
# solve_poisson!(pm,rho,Phi)
# E1,E2,E3=gradient(pm,-Phi)
#
# J1[:]=0.0;J2[:]=0.0;J3[:]=0.0
#
 # integrate_Hp_midpoint(pm,x1n,x2n,x3n,v1n,v2n,v3n,qn,mn,wn,
                #  J1,J2,J3,B1,B2,B3,dt)
# #J1[pm.N1n,:,:].=0; J2[:,pm.N2n,:].=0; J3[:,:,pm.N3n].=0
# J1[pm.N1n,pm.N2n,pm.N3n].=0
# J2[pm.N1n,pm.N2n,pm.N3n].=0
# J3[pm.N1n,pm.N2n,pm.N3n].=0
# E1.-=J1/Np;E2.-=J2/Np; E3.-=J3/Np
# #
# # rho=accum_osde(pm,x1n,x2n,x3n, wn.*qn)./Np
# # solve_poisson!(pm,rho*gamma,Phi)
# # E1_,E2_,E3_=gradient(pm,-Phi)
# #
# # print(E1-E1_, "\n" )
# # print(J1/Np,"\n")
rho=accum_osde(pm,x1n,x2n,x3n, wn.*qn)./Np
rho=MPI.allreduce(rho, MPI.SUM, mpi_comm)

print("Gauss Error: ",gauss_error(pm,rho,E1,E2,E3),"\n")
#########################################################
# Variance
# rho_var=accum_osde_var(pm,rho,x1n,x2n,x3n, wn.*qn)./Np
#
# rdx=sortperm(abs.(rho)[:])
#
# semilogy(abs.(rho[rdx]))
# plot(rho_var[rdx]/sqrt(Np))

# Only plot in repl
if (mpi_rk==0 && mpi_sz==1)

energy=sum(Epot,2)+sum(Ekin,2)+sum(Bpot,2)

using PyPlot

figure()
semilogy(ttime,abs.((energy-energy[1])./energy[1]))
grid()
ylabel("rel. energy error")
xlabel("time ")

# plot(ttime, Ekin[:,1]+Epot[:,1])
#
figure()
 semilogy(ttime, Epot)
 figure()
 semilogy(ttime, Bpot)

# figure()
# semilogy(ttime, Bpot)

moment_error=abs.(Momentum[:,:].-reshape(Momentum[1,:],1,3));

figure()
semilogy(ttime,moment_error)
xlabel("time "); grid()
ylabel("absolute momentum error")

else
  ############# Finalize MPI
  MPI.Finalize()

end

max_stars_repo_name: ameresj/VlasovPIF

max_issues_repo_name: ameresj/VlasovPIF

max_forks_repo_name: ameresj/VlasovPIF


text:
#=

    Top level test suite, which runs unit tests, Pkg dependencies, regression & comparison tests.

=#

using Test

run(`cd ../`)

if VERSION >= v"1.6.0"

    println("Starting All tests...")

    @testset "All tests" begin
        include("test_c_sw.jl")

        # Both of the test suites below are a work in progress.
        # include("unit_tests.jl")
        # include("PkgTests.jl")
    end

else

    println("Current version of Julia is not supported, please update to v1.6")

end


text:
#
# this file shows an example of pseudo-spectral discretization for PHS systems for wave equation
#  three cases are presented. It should be used for verifying the discretization, eigenvalues,
#  and constraint elimination methods. Three cases are tested:
#         1) PHS with u = 0 (Neumann-Dirichlet B.C)
#         2) PHS with y = 0 (Neumann-Dirichlet B.C, but with a constraint)
#         3) PHS with y = 0 after removing the constraints.
#

# auxiliary function: used for printing results

function print_table(table)
	num_y = size(table,1)
	num_x = size(table,2)
	for i = 1:num_y
		for j = 1:num_x
			if j < 3
				@printf("%.2f\t", table[i,j])
			else
				@printf("%.2e", table[i,j])
			end
		end
		@printf("\n")
	end
end

# modules:
using PortHamiltonian
#using PyPlot

N = 10;
a = 0.;
b = 1.;

# load the phs discretization structure
## the Hamiltonian is considered linear with Q = I
ph = discrete_phs(N, a, b);

eigval, eigvec = eig(ph);
num_freq = frequencies(eigval)

exact_freq = pi/2*collect(1:2:2*N);
err = exact_freq - num_freq;

# comparison between numerical and exact results:
comp_table = [num_freq exact_freq err./exact_freq];
comp_table = comp_table[1:10,:]
@printf("\nNatural frequencies of the Dirichlet-Neumann B.C. wave equation (with u = 0)\n")
@printf("exact\tnumeric\terror \n")
print_table(comp_table)
@printf("\t From a total %d frequencies, %d are represented with error less than 0.01\n", length(err), sum(abs(err) .< 0.01))
@printf("\n\t %d are smaller than %.1e\n",sum(abs(err) .< eps()*1e6),eps()*1e6)
#plot the first five modes
#plot(ph.disc_data.flow.xi, real(eigvec[1:N,N+2]))

# constrained version
ph_constrained = deepcopy(ph);
set_constraint!(ph_constrained, ph_constrained.B[:,1], ph_constrained.D[1:1,1:1]);

## eigenvalues of the constrained
eigval, eigvec = eig(ph_constrained)
num_freq = frequencies(eigval)
exact_freq = pi*collect(-0:1:N-1);
err = exact_freq - num_freq;

# comparison between numerical and exact results:
comp_table = [num_freq exact_freq err./exact_freq]
comp_table = comp_table[1:10,:]
@printf("\nNatural frequencies of the Neumann-Dirichlet B.C. wave equation (with y = 0)\n")
@printf("exact\tnumeric\terror \n")
print_table(comp_table)
@printf("\t\t From a total %d frequencies, %d are represented with error less than 0.01\n", length(err), sum(abs(err) .< 0.01))
@printf("\n\t %d are smaller than %.1e\n",sum(abs(err) .< eps()*1e6),eps()*1e6)
## find an explicit equivalent system
ph_c_elim = constraint_elimination(ph_constrained)
eigval, eigvec = eig(ph_c_elim);
num_freq = frequencies(eigval)

exact_freq = pi*collect(0:1:N-1);
err = exact_freq - num_freq;
# comparison between numerical and exact results:
comp_table = [num_freq exact_freq err./exact_freq]
comp_table = comp_table[1:10,:]
@printf("\nNatural frequencies of the Neumann-Dirichlet B.C. wave equation (with y=0 after constraint elimination)\n")
@printf("exact\tnumeric\terror \n")
print_table(comp_table)
@printf("\t\t From a total %d frequencies, %d are represented with error less than 0.01", length(err), sum(abs(err) .< 0.01))
@printf("\n\t %d are smaller than %.1e\n",sum(abs(err) .< eps()*1e6),eps()*1e6)

max_stars_repo_name: flavioluiz/phjulia

max_issues_repo_name: flavioluiz/phjulia

max_forks_repo_name: flavioluiz/phjulia


text:
function modindex(i::Integer, m::Integer)
    x = mod(i, m)
    if x == 0
        return m
    else
        return x
    end
end

function twoloop!(g_x::Vector,
                  rho::Vector,
                  s::Matrix,
                  y::Matrix,
                  m::Integer,
                  iteration::Integer,
                  alpha::Vector,
                  q::Vector,
                  p::Vector)
    # Copy gradient into q for backward pass
    copy!(q, g_x)

    n = length(p)
    @assert length(q) == n

    upper = iteration - 1
    lower = iteration - m
    for index in upper:-1:lower
        if index < 1
            continue
        end
        i = modindex(index, m)
        alpha[i] = rho[i] * dot(s[:, i], q)
        for j in 1:n
            q[j] -= alpha[i] * y[j, i]
        end
    end

    # Copy q into p for forward pass
    copy!(p, q)

    for index in lower:1:upper
        if index < 1
            continue
        end
        i = modindex(index, m)
        beta = rho[i] * dot(y[:, i], p)
        for j in 1:n
            p[j] += s[j, i] * (alpha[i] - beta)
        end
    end

    for j in 1:n
        p[j] = -1.0 * p[j]
    end

    return
end

function l_bfgs_trace!(tr::OptimizationTrace,
                       x::Vector,
                       f_x::Real,
                       g_new::Vector,
                       iteration::Integer,
                       alpha::Real,
                       store_trace::Bool,
                       show_trace::Bool)
    dt = Dict()
    dt["g(x)"] = copy(g_new)
    dt["Step-size"] = alpha
    dt["First-order opt"] = norm(g_new, Inf)
    os = OptimizationState(copy(x), f_x, iteration, dt)
    if store_trace
        push!(tr, os)
    end
    if show_trace
        println(os)
    end
end

function l_bfgs(d::DifferentiableFunction,
                initial_x::Vector;
                m::Integer = 10,
                tolerance::Real = 1e-8,
                iterations::Integer = 1_000,
                store_trace::Bool = false,
                show_trace::Bool = false)

    # Set iteration counter
    iteration = 0

    # Keep a record of the starting point
    x = copy(initial_x)

    # Establish size of parameter space
    n = length(x)

    # Initialize rho, s and y
    rho = Array(Float64, m)
    s = Array(Float64, n, m)
    y = Array(Float64, n, m)

    # Initialize q, p, tmp_s and tmp_y
    q = Array(Float64, n)
    p = Array(Float64, n)
    tmp_s = Array(Float64, n)
    tmp_y = Array(Float64, n)

    # Reuse storage during calls to backtracking line search
    ls_x = Array(Float64, n)
    ls_gradient = Array(Float64, n)

    # Re-use this vector during every call to two_loop()
    twoloop_v = Array(Float64, m)

    # Compute the initial values of f and g
    g_x = Array(Float64, n)
    g_new = Array(Float64, n)
    f_calls = 1
    g_calls = 1
    f_x = d.fg!(x, g_new)
    copy!(g_x, g_new)

    # Print trace information
    tr = OptimizationTrace()

    # Show trace
    if store_trace || show_trace
        l_bfgs_trace!(tr, x, f_x, g_new, iteration, 0.0, store_trace, show_trace)
    end

    # Iterate until convergence
    converged = false

    while !converged && iteration <= iterations
        # Update count
        iteration += 1

        # Select a search direction
        twoloop!(g_new, rho, s, y, m, iteration, twoloop_v, q, p)

        # Select a step-size
        alpha, f_update, g_update =
          backtracking_line_search!(d, x, p, ls_x, ls_gradient)
        f_calls += f_update
        g_calls += g_update

        # Update position
        for i in 1:n
            tmp_s[i] = alpha * p[i]
            x[i] = x[i] + tmp_s[i]
        end
        copy!(g_x, g_new)
        f_x = d.fg!(x, g_new)
        f_calls += 1
        g_calls += 1

        # Estimate movement
        for i in 1:n
            tmp_y[i] = g_new[i] - g_x[i]
        end
        tmp_rho = 1.0 / dot(tmp_y, tmp_s)
        if isinf(tmp_rho)
            @printf "Cannot decrease the objective function along the current search direction\n"
            break
        end

        # Keep a record of the new s, y and rho
        s[:, modindex(iteration, m)] = tmp_s
        y[:, modindex(iteration, m)] = tmp_y
        rho[modindex(iteration, m)] = tmp_rho

        # Assess convergence
        if norm(g_x, Inf) <= tolerance
           converged = true
        end

        # Show trace
        if store_trace || show_trace
            l_bfgs_trace!(tr, x, f_x, g_new, iteration, alpha, store_trace, show_trace)
        end
    end

    OptimizationResults("L-BFGS",
                        initial_x,
                        x,
                        f_x,
                        iteration,
                        converged,
                        tr,
                        f_calls,
                        g_calls)
end


text:
using Documenter
using SubstitutionModels

makedocs(
    sitename = "SubstitutionModels",
    format = Documenter.HTML(),
    modules = [SubstitutionModels],
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Substitution models" => "man/models.md",
            "Q and P matrices" => "man/pqmatrices.md",
            "Provided models & custom models" => "man/modeltypes.md"
        ],
        "Contributing" => "contributing.md"
    ],
    authors = "Justin Angevarre & Ben J. Ward, at the BioJulia organisation, and contributors."
)

deploydocs(
    repo = "github.com/BioJulia/SubstitutionModels.jl.git",
    push_preview = true
)


text:
using LocalBinaryPatterns
using Documenter

format = Documenter.HTML(;
    prettyurls=get(ENV, "CI", "false") == "true",
    assets=String[],
)
makedocs(;
    sitename="LocalBinaryPatterns Documentation",
    modules=[LocalBinaryPatterns],
    format=format,
    pages=[
        "Home" => "index.md",
        "API References" => "reference.md"
    ],
)

deploydocs(;repo="github.com/johnnychen94/LocalBinaryPatterns.jl")


text:
###
### Particle Filtering and Particle MCMC Samplers.
###

####
#### Generic Sequential Monte Carlo sampler.
####

"""
    SMC(n_particles::Int)

Sequential Monte Carlo sampler.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

Usage:

```julia
SMC(1000)
```
"""
mutable struct SMC{space, F} <: InferenceAlgorithm
    n_particles           ::  Int
    resampler             ::  F
    resampler_threshold   ::  Float64
end
function SMC(n_particles::Int, resampler::F, resampler_threshold::Float64, space::Tuple) where {F}
    return SMC{space, F}(n_particles, resampler, resampler_threshold)
end
SMC(n) = SMC(n, resample_systematic, 0.5, ())
SMC(n_particles::Int, ::Tuple{}) = SMC(n_particles)
function SMC(n_particles::Int, space::Symbol...)
    SMC(n_particles, resample_systematic, 0.5, space)
end

function Sampler(alg::SMC, s::Selector)
    info = Dict{Symbol, Any}()
    info[:logevidence] = []
    return Sampler(alg, info, s)
end

function step(model, spl::Sampler{<:SMC}, vi::VarInfo)
    particles = ParticleContainer{Trace}(model)
    vi.num_produce = 0;  # Reset num_produce before new sweep\.
    set_retained_vns_del_by_spl!(vi, spl)
    resetlogp!(vi)

    push!(particles, spl.alg.n_particles, spl, vi)

    while consume(particles) != Val{:done}
      ess = effectiveSampleSize(particles)
      if ess <= spl.alg.resampler_threshold * length(particles)
        resample!(particles,spl.alg.resampler)
      end
    end

    ## pick a particle to be retained.
    Ws, _ = weights(particles)
    indx = randcat(Ws)
    push!(spl.info[:logevidence], particles.logE)

    return particles[indx].vi, true
end

## wrapper for smc: run the sampler, collect results.
function sample(model::Model, alg::SMC)
    spl = Sampler(alg)

    particles = ParticleContainer{Trace}(model)
    push!(particles, spl.alg.n_particles, spl, empty!(VarInfo(model)))

    while consume(particles) != Val{:done}
      ess = effectiveSampleSize(particles)
      if ess <= spl.alg.resampler_threshold * length(particles)
        resample!(particles,spl.alg.resampler)
      end
    end
    w, samples = getsample(particles)
    res = Chain(log(w), samples)
    return res
end

####
#### Particle Gibbs sampler.
####

"""
    PG(n_particles::Int, n_iters::Int)

Particle Gibbs sampler.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

Usage:

```julia
PG(100, 100)
```
"""
mutable struct PG{space, F} <: InferenceAlgorithm
  n_particles           ::    Int         # number of particles used
  n_iters               ::    Int         # number of iterations
  resampler             ::    F           # function to resample
end
function PG(n_particles::Int, n_iters::Int, resampler::F, space::Tuple) where F
    return PG{space, F}(n_particles, n_iters, resampler)
end
PG(n1::Int, n2::Int, ::Tuple{}) = PG(n1, n2)
function PG(n1::Int, n2::Int, space::Symbol...)
    PG(n1, n2, resample_systematic, space)
end

const CSMC = PG # type alias of PG as Conditional SMC

function Sampler(alg::PG, s::Selector)
    info = Dict{Symbol, Any}()
    info[:logevidence] = []
    Sampler(alg, info, s)
end

step(model, spl::Sampler{<:PG}, vi::VarInfo, _) = step(model, spl, vi)

function step(model, spl::Sampler{<:PG}, vi::VarInfo)
    particles = ParticleContainer{Trace{typeof(spl), typeof(vi), typeof(model)}}(model)

    vi.num_produce = 0;  # Reset num_produce before new sweep\.
    ref_particle = isempty(vi) ?
                  nothing :
                  forkr(Trace(model, spl, vi))

    set_retained_vns_del_by_spl!(vi, spl)
    resetlogp!(vi)

    if ref_particle == nothing
        push!(particles, spl.alg.n_particles, spl, vi)
    else
        push!(particles, spl.alg.n_particles-1, spl, vi)
        push!(particles, ref_particle)
    end

    while consume(particles) != Val{:done}
        resample!(particles, spl.alg.resampler, ref_particle)
    end

    ## pick a particle to be retained.
    Ws, _ = weights(particles)
    indx = randcat(Ws)
    push!(spl.info[:logevidence], particles.logE)

    return particles[indx].vi, true
end

function sample(  model::Model,
                  alg::PG;
                  save_state=false,         # flag for state saving
                  resume_from=nothing,      # chain to continue
                  reuse_spl_n=0             # flag for spl re-using
                )

    spl = reuse_spl_n > 0 ?
          resume_from.info[:spl] :
          Sampler(alg)
    if resume_from != nothing
        spl.selector = resume_from.info[:spl].selector
    end
    @assert typeof(spl.alg) == typeof(alg) "[Turing] alg type mismatch; please use resume() to re-use spl"

    n = reuse_spl_n > 0 ?
        reuse_spl_n :
        alg.n_iters
    samples = Vector{Sample}()

    ## custom resampling function for pgibbs
    ## re-inserts reteined particle after each resampling step
    time_total = zero(Float64)

    vi = resume_from == nothing ?
        empty!(VarInfo(model)) :
        resume_from.info[:vi]

    pm = nothing
    PROGRESS[] && (spl.info[:progress] = ProgressMeter.Progress(n, 1, "[PG] Sampling...", 0))

    for i = 1:n
        time_elapsed = @elapsed vi, _ = step(model, spl, vi)
        push!(samples, Sample(vi))
        samples[i].value[:elapsed] = time_elapsed

        time_total += time_elapsed

        if PROGRESS[] && spl.selector.tag == :default
            ProgressMeter.next!(spl.info[:progress])
        end
    end

    @info("[PG] Finished with")
    @info("  Running time    = $time_total;")

    loge = exp.(mean(spl.info[:logevidence]))
    if resume_from != nothing   # concat samples
        pushfirst!(samples, resume_from.info[:samples]...)
        pre_loge = exp.(resume_from.logevidence)
        # Calculate new log-evidence
        pre_n = length(resume_from.info[:samples])
        loge = (log(pre_loge) * pre_n + log(loge) * n) / (pre_n + n)
    end
    c = Chain(loge, samples)       # wrap the result by Chain

    if save_state               # save state
        c = save(c, spl, model, vi, samples)
    end

    return c
end

function assume(  spl::Sampler{T},
                  dist::Distribution,
                  vn::VarName,
                  _::VarInfo
                ) where T<:Union{PG,SMC}

    vi = current_trace().vi
    if isempty(getspace(spl.alg)) || vn.sym in getspace(spl.alg)
        if ~haskey(vi, vn)
            r = rand(dist)
            push!(vi, vn, r, dist, spl)
        elseif is_flagged(vi, vn, "del")
            unset_flag!(vi, vn, "del")
            r = rand(dist)
            vi[vn] = vectorize(dist, r)
            setgid!(vi, spl.selector, vn)
            setorder!(vi, vn, vi.num_produce)
        else
            updategid!(vi, vn, spl)
            r = vi[vn]
        end
    else # vn belongs to other sampler <=> conditionning on vn
        if haskey(vi, vn)
            r = vi[vn]
        else
            r = rand(dist)
            push!(vi, vn, r, dist, Selector(:invalid))
        end
        acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))
    end
    return r, zero(Real)
end

function assume(  spl::Sampler{A},
                  dists::Vector{D},
                  vn::VarName,
                  var::Any,
                  vi::VarInfo
                ) where {A<:Union{PG,SMC},D<:Distribution}
    error("[Turing] PG and SMC doesn't support vectorizing assume statement")
end

function observe(spl::Sampler{T}, dist::Distribution, value, vi) where T<:Union{PG,SMC}
    produce(logpdf(dist, value))
    return zero(Real)
end

function observe( spl::Sampler{A},
                  ds::Vector{D},
                  value::Any,
                  vi::VarInfo
                ) where {A<:Union{PG,SMC},D<:Distribution}
    error("[Turing] PG and SMC doesn't support vectorizing observe statement")
end

####
#### Particle marginal Metropolis-Hastings sampler.
####

"""
    PMMH(n_iters::Int, smc_alg:::SMC, parameters_algs::Tuple{MH})

Particle independant Metropolis–Hastings and
Particle marginal Metropolis–Hastings samplers.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

Usage:

```julia
alg = PMMH(100, SMC(20, :v1), MH(1,:v2))
alg = PMMH(100, SMC(20, :v1), MH(1,(:v2, (x) -> Normal(x, 1))))
```

Arguments:

- `n_iters::Int` : Number of iterations to run.
- `smc_alg:::SMC` : An [`SMC`](@ref) algorithm to use.
- `parameters_algs::Tuple{MH}` : An [`MH`](@ref) algorithm, which includes a
sample space specification.
"""
mutable struct PMMH{space, A<:Tuple} <: InferenceAlgorithm
    n_iters               ::    Int               # number of iterations
    algs                  ::    A                 # Proposals for state & parameters
end
function PMMH(n_iters::Int, algs::A, space::Tuple) where {A <: Tuple}
    return PMMH{space, A}(n_iters, algs)
end
function PMMH(n_iters::Int, smc_alg::SMC, parameter_algs...)
    return PMMH(n_iters, tuple(parameter_algs..., smc_alg), ())
end

PIMH(n_iters::Int, smc_alg::SMC) = PMMH(n_iters, tuple(smc_alg), ())

function Sampler(alg::PMMH, model::Model, s::Selector)
    info = Dict{Symbol, Any}()
    spl = Sampler(alg, info, s)

    alg_str = "PMMH"
    n_samplers = length(alg.algs)
    samplers = Array{Sampler}(undef, n_samplers)

    space = Set{Symbol}()

    for i in 1:n_samplers
        sub_alg = alg.algs[i]
        if isa(sub_alg, Union{SMC, MH})
            samplers[i] = Sampler(sub_alg, model, Selector(Symbol(typeof(sub_alg))))
        else
            error("[$alg_str] unsupport base sampling algorithm $alg")
        end
        if typeof(sub_alg) == MH && sub_alg.n_iters != 1
            @warn("[$alg_str] number of iterations greater than 1 is useless for MH since it is only used for its proposal")
        end
        space = (space..., getspace(sub_alg)...)
    end

    # Sanity check for space
    if !isempty(space)
        @assert issubset(get_pvars(model), space) "[$alg_str] symbols specified to samplers ($space) doesn't cover the model parameters ($(get_pvars(model)))"

        if !(issetequal(get_pvars(model), space))
            @warn("[$alg_str] extra parameters specified by samplers don't exist in model: $(setdiff(space, get_pvars(model)))")
        end
    end

    info[:old_likelihood_estimate] = -Inf # Force to accept first proposal
    info[:old_prior_prob] = 0.0
    info[:samplers] = samplers

    return spl
end

function step(model, spl::Sampler{<:PMMH}, vi::VarInfo, is_first::Bool)
    violating_support = false
    proposal_ratio = 0.0
    new_prior_prob = 0.0
    new_likelihood_estimate = 0.0
    old_θ = copy(vi[spl])

    Turing.DEBUG && @debug "Propose new parameters from proposals..."
    for local_spl in spl.info[:samplers][1:end-1]
        Turing.DEBUG && @debug "$(typeof(local_spl)) proposing $(getspace(local_spl.alg))..."
        propose(model, local_spl, vi)
        if local_spl.info[:violating_support] violating_support=true; break end
        new_prior_prob += local_spl.info[:prior_prob]
        proposal_ratio += local_spl.info[:proposal_ratio]
    end

    if !violating_support # do not run SMC if going to refuse anyway
        Turing.DEBUG && @debug "Propose new state with SMC..."
        vi, _ = step(model, spl.info[:samplers][end], vi)
        new_likelihood_estimate = spl.info[:samplers][end].info[:logevidence][end]

        Turing.DEBUG && @debug "computing accept rate α..."
        is_accept, _ = mh_accept(
          -(spl.info[:old_likelihood_estimate] + spl.info[:old_prior_prob]),
          -(new_likelihood_estimate + new_prior_prob),
          proposal_ratio,
        )
    end

    Turing.DEBUG && @debug "decide whether to accept..."
    if !violating_support && is_accept # accepted
        is_accept = true
        spl.info[:old_likelihood_estimate] = new_likelihood_estimate
        spl.info[:old_prior_prob] = new_prior_prob
    else                      # rejected
        is_accept = false
        vi[spl] = old_θ
    end

    return vi, is_accept
end

function sample(  model::Model,
                  alg::PMMH;
                  save_state=false,         # flag for state saving
                  resume_from=nothing,      # chain to continue
                  reuse_spl_n=0             # flag for spl re-using
                )

    spl = Sampler(alg, model)
    if resume_from != nothing
        spl.selector = resume_from.info[:spl].selector
    end
    alg_str = "PMMH"

    # Number of samples to store
    sample_n = spl.alg.n_iters

    # Init samples
    time_total = zero(Float64)
    samples = Array{Sample}(undef, sample_n)
    weight = 1 / sample_n
    for i = 1:sample_n
        samples[i] = Sample(weight, Dict{Symbol, Any}())
    end

    # Init parameters
    vi = if resume_from == nothing
        vi_ = VarInfo(model)
    else
        resume_from.info[:vi]
    end
    n = spl.alg.n_iters

    # PMMH steps
    accept_his = Bool[]
    PROGRESS[] && (spl.info[:progress] = ProgressMeter.Progress(n, 1, "[$alg_str] Sampling...", 0))
    for i = 1:n
      Turing.DEBUG && @debug "$alg_str stepping..."
      time_elapsed = @elapsed vi, is_accept = step(model, spl, vi, i==1)

      if is_accept # accepted => store the new predcits
          samples[i].value = Sample(vi, spl).value
      else         # rejected => store the previous predcits
          samples[i] = samples[i - 1]
      end

      time_total += time_elapsed
      push!(accept_his, is_accept)
      if PROGRESS[]
        haskey(spl.info, :progress) && ProgressMeter.update!(spl.info[:progress], spl.info[:progress].counter + 1)
      end
    end

    println("[$alg_str] Finished with")
    println("  Running time    = $time_total;")
    accept_rate = sum(accept_his) / n  # calculate the accept rate
    println("  Accept rate         = $accept_rate;")

    if resume_from != nothing   # concat samples
      pushfirst!(samples, resume_from.info[:samples]...)
    end
    c = Chain(-Inf, samples)       # wrap the result by Chain

    if save_state               # save state
      c = save(c, spl, model, vi, samples)
    end

    c
end

####
#### IMCMC Sampler.
####

"""
    IPMCMC(n_particles::Int, n_iters::Int, n_nodes::Int, n_csmc_nodes::Int)

Particle Gibbs sampler.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

Usage:

```julia
IPMCMC(100, 100, 4, 2)
```

Arguments:

- `n_particles::Int` : Number of particles to use.
- `n_iters::Int` : Number of iterations to employ.
- `n_nodes::Int` : The number of nodes running SMC and CSMC.
- `n_csmc_nodes::Int` : The number of CSMC nodes.
```

A paper on this can be found [here](https://arxiv.org/abs/1602.05128).
"""
mutable struct IPMCMC{space, F} <: InferenceAlgorithm
  n_particles           ::    Int         # number of particles used
  n_iters               ::    Int         # number of iterations
  n_nodes               ::    Int         # number of nodes running SMC and CSMC
  n_csmc_nodes          ::    Int         # number of nodes CSMC
  resampler             ::    F           # function to resample
end
function IPMCMC(n_particles::Int, n_iters::Int, n_nodes::Int, n_csmc_nodes::Int, resampler::F, space::Tuple) where {F}
    return IPMCMC{space, F}(n_particles, n_iters, n_nodes, n_csmc_nodes, resampler)
end
IPMCMC(n1::Int, n2::Int) = IPMCMC(n1, n2, 32, 16)
IPMCMC(n1::Int, n2::Int, n3::Int) = IPMCMC(n1, n2, n3, Int(ceil(n3/2)))
IPMCMC(n1::Int, n2::Int, n3::Int, n4::Int, ::Tuple{}) = IPMCMC(n1, n2, n3, n4)
function IPMCMC(n1::Int, n2::Int, n3::Int, n4::Int, space::Symbol...)
  IPMCMC(n1, n2, n3, n4, resample_systematic, space)
end

function Sampler(alg::IPMCMC, s::Selector)
  info = Dict{Symbol, Any}()
  spl = Sampler(alg, info, s)
  # Create SMC and CSMC nodes
  samplers = Array{Sampler}(undef, alg.n_nodes)
  # Use resampler_threshold=1.0 for SMC since adaptive resampling is invalid in this setting
  default_CSMC = CSMC(alg.n_particles, 1, alg.resampler, getspace(alg))
  default_SMC = SMC(alg.n_particles, alg.resampler, 1.0, getspace(alg))

  for i in 1:alg.n_csmc_nodes
    samplers[i] = Sampler(default_CSMC, Selector(Symbol(typeof(default_CSMC))))
  end
  for i in (alg.n_csmc_nodes+1):alg.n_nodes
    samplers[i] = Sampler(default_SMC, Selector(Symbol(typeof(default_CSMC))))
  end

  info[:samplers] = samplers

  return spl
end

function step(model, spl::Sampler{<:IPMCMC}, VarInfos::Array{VarInfo}, is_first::Bool)
    # Initialise array for marginal likelihood estimators
    log_zs = zeros(spl.alg.n_nodes)

    # Run SMC & CSMC nodes
    for j in 1:spl.alg.n_nodes
        VarInfos[j].num_produce = 0
        VarInfos[j] = step(model, spl.info[:samplers][j], VarInfos[j])[1]
        log_zs[j] = spl.info[:samplers][j].info[:logevidence][end]
    end

    # Resampling of CSMC nodes indices
    conditonal_nodes_indices = collect(1:spl.alg.n_csmc_nodes)
    unconditonal_nodes_indices = collect(spl.alg.n_csmc_nodes+1:spl.alg.n_nodes)
    for j in 1:spl.alg.n_csmc_nodes
        # Select a new conditional node by simulating cj
        log_ksi = vcat(log_zs[unconditonal_nodes_indices], log_zs[j])
        ksi = exp.(log_ksi .- maximum(log_ksi))
        c_j = wsample(ksi) # sample from Categorical with unormalized weights

        if c_j < length(log_ksi) # if CSMC node selects another index than itself
            conditonal_nodes_indices[j] = unconditonal_nodes_indices[c_j]
            unconditonal_nodes_indices[c_j] = j
        end
    end
    nodes_permutation = vcat(conditonal_nodes_indices, unconditonal_nodes_indices)

    VarInfos[nodes_permutation]
end

function sample(model::Model, alg::IPMCMC)

  spl = Sampler(alg)

  # Number of samples to store
  sample_n = alg.n_iters * alg.n_csmc_nodes

  # Init samples
  time_total = zero(Float64)
  samples = Array{Sample}(undef, sample_n)
  weight = 1 / sample_n
  for i = 1:sample_n
    samples[i] = Sample(weight, Dict{Symbol, Any}())
  end

  # Init parameters
  vi = empty!(VarInfo(model))
  VarInfos = Array{VarInfo}(undef, spl.alg.n_nodes)
  for j in 1:spl.alg.n_nodes
    VarInfos[j] = deepcopy(vi)
  end
  n = spl.alg.n_iters

  # IPMCMC steps
  if PROGRESS[] spl.info[:progress] = ProgressMeter.Progress(n, 1, "[IPMCMC] Sampling...", 0) end
  for i = 1:n
    Turing.DEBUG && @debug "IPMCMC stepping..."
    time_elapsed = @elapsed VarInfos = step(model, spl, VarInfos, i==1)

    # Save each CSMS retained path as a sample
    for j in 1:spl.alg.n_csmc_nodes
      samples[(i-1)*alg.n_csmc_nodes+j].value = Sample(VarInfos[j], spl).value
    end

    time_total += time_elapsed
    if PROGRESS[]
      haskey(spl.info, :progress) && ProgressMeter.update!(spl.info[:progress], spl.info[:progress].counter + 1)
    end
  end

  println("[IPMCMC] Finished with")
  println("  Running time    = $time_total;")

  Chain(0.0, samples) # wrap the result by Chain
end

####
#### Resampling schemes for particle filters
####

# Some references
#  - http://arxiv.org/pdf/1301.4019.pdf
#  - http://people.isy.liu.se/rt/schon/Publications/HolSG2006.pdf
# Code adapted from: http://uk.mathworks.com/matlabcentral/fileexchange/24968-resampling-methods-for-particle-filtering

# Default resampling scheme
function resample(w::AbstractVector{<:Real}, num_particles::Integer=length(w))
    return resample_systematic(w, num_particles)
end

# More stable, faster version of rand(Categorical)
function randcat(p::AbstractVector{T}) where T<:Real
    r, s = rand(T), 1
    for j in eachindex(p)
        r -= p[j]
        if r <= zero(T)
            s = j
            break
        end
    end
    return s
end

function resample_multinomial(w::AbstractVector{<:Real}, num_particles::Integer)
    return rand(Distributions.sampler(Categorical(w)), num_particles)
end

function resample_residual(w::AbstractVector{<:Real}, num_particles::Integer)

    M = length(w)

    # "Repetition counts" (plus the random part, later on):
    Ns = floor.(length(w) .* w)

    # The "remainder" or "residual" count:
    R = Int(sum(Ns))

    # The number of particles which will be drawn stocastically:
    M_rdn = num_particles - R

    # The modified weights:
    Ws = (M .* w - floor.(M .* w)) / M_rdn

    # Draw the deterministic part:
    indx1, i = Array{Int}(undef, R), 1
    for j in 1:M
        for k in 1:Ns[j]
            indx1[i] = j
            i += 1
        end
    end

    # And now draw the stocastic (Multinomial) part:
    return append!(indx1, rand(Distributions.sampler(Categorical(w)), M_rdn))
end

function resample_stratified(w::AbstractVector{<:Real}, num_particles::Integer)

    Q, N = cumsum(w), num_particles

    T = Array{Float64}(undef, N + 1)
    for i=1:N,
        T[i] = rand() / N + (i - 1) / N
    end
    T[N+1] = 1

    indx, i, j = Array{Int}(undef, N), 1, 1
    while i <= N
        if T[i] < Q[j]
            indx[i] = j
            i += 1
        else
            j += 1
        end
    end
    return indx
end

function resample_systematic(w::AbstractVector{<:Real}, num_particles::Integer)

    Q, N = cumsum(w), num_particles

    T = collect(range(0, stop = maximum(Q)-1/N, length = N)) .+ rand()/N
    push!(T, 1)

    indx, i, j = Array{Int}(undef, N), 1, 1
    while i <= N
        if T[i] < Q[j]
            indx[i] = j
            i += 1
        else
            j += 1
        end
    end
    return indx
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

# required functions for each geotrait
include("geotraits/domain.jl")
include("geotraits/data.jl")

# fallback implementations
include("geotraits/fallback.jl")


text:
@safetestset "<select> test" begin

  using Genie
  using Genie.Renderer.Html
  import Genie.Util: fws

  @test html(filepath("formview.jl.html")).body |> String |> fws ==
        """<!DOCTYPE html><form method="POST" enctype="multipart/form-data" action="/new">

Example select
            <select class="form-control" id="exampleFormControlSelect1"><option>1</option><option>2</option>
            <option>3</option><option>4</option><option>5</option></select>

            <button class="btn btn-primary" type="submit">Submit</button></form>""" |> fws

end


text:
module HWfuncapp

using FastGaussQuadrature # to get chebyshevnodes

# you dont' have to use PyPlot. I did much of it in PyPlot, hence
# you will see me often qualify code with `PyPlot.plot` or so.
using PyPlot
import ApproXD: getBasis, BSpline
using Distributions
using ApproxFun
using Plots
export TT, ChebyT, q1

ChebyT(x,deg) = cos(acos(x)*deg)
unitmap(x,lb,ub) = 2 .* (x .- lb) ./ (ub .- lb) .- 1	#[a,b] -> [-1,1]

n=15
map(x,lb,ub) = 0.5 .* (lb .+ ub) .+ 0.5 .* (ub .- lb) .* x

function Cheby(x,deg)
len=length(x)
psi = ones(len,deg)
for i in 1:len
for j in 1:deg
psi[i,j] = ChebyT(x[i],j-1)
end
end
return psi
end

function q1(n=15)
lb, ub, inter, Nod = -3, 3, 100, gausschebyshev(n)[1]#Values for evaluation
x = range(lb,stop = ub,length = n)#
psi = map(Nod, lb, ub) #Evaluation

f(x) = x .+ 2x.^2 - exp.(-x)#Function to evaluate
PSI=inv(Cheby(Nod,n))
ev=f(psi)
ResEv=PSI*ev

testX=range(lb, stop=ub, length=inter)
testY=f(testX)

testMap=map(testX, lb, ub)
testPsi=Cheby(unitmap(testMap), n)

ResEv=testPsi*PSI

testError=ev-ResEv

#PyPlot.plot()
p = Plots.plot(layout = 2, dpi = 400)
Plots.plot!(p[1],testX,testY,label = "Test Y",lw = 1,linecolor = "black")
Plots.plot!(p[1],testX,ResEv, label = "Result Evaluation", lw = 2,linestyle = :dot, linecolor = "pink")
Plots.plot!(p[2],testX,testError, label = "Test Error", lw = 1, linecolor = "blue")

# without using PyPlot, just erase the `PyPlot.` part
PyPlot.savefig(joinpath(dirname(@__FILE__),"..","q1.png"))
return Dict(:error=>maximum(abs,err))
end

end # module


text:
# Test for external sort

using Base.Test

include("es.jl")
include("helper.jl")

@testset "External Sort Tests" begin

    # Test not enough memory
    @test_throws ErrorException external_sort(1)

    ############################ File with 12 integers ############################

    create_file(12)

    # "Unlimited" memory test
    external_sort(36)
    @test check_sort()

    # Basic test
    external_sort(4)
    @test check_sort()

    # Another basic test
    external_sort(6)
    @test check_sort()

    # Test with no divisible number
    external_sort(5)
    @test check_sort()

end

max_stars_repo_name: rgaiacs/external-sort

max_issues_repo_name: rgaiacs/external-sort

max_forks_repo_name: rgaiacs/external-sort


text:
function CombinedObservationData(arg0::CombinationType, arg1::MeasurementType, arg2::jdouble, arg3::jdouble, arg4::List)
    return CombinedObservationData((CombinationType, MeasurementType, jdouble, jdouble, List), arg0, arg1, arg2, arg3, arg4)
end

function get_combination_type(obj::CombinedObservationData)
    return jcall(obj, "getCombinationType", CombinationType, ())
end

function get_combined_m_hz_frequency(obj::CombinedObservationData)
    return jcall(obj, "getCombinedMHzFrequency", jdouble, ())
end

function get_measurement_type(obj::CombinedObservationData)
    return jcall(obj, "getMeasurementType", MeasurementType, ())
end

function get_used_observation_data(obj::CombinedObservationData)
    return jcall(obj, "getUsedObservationData", List, ())
end

function get_value(obj::CombinedObservationData)
    return jcall(obj, "getValue", jdouble, ())
end


text:
"""
    get_quantiles(x, c::SequentialSamplingConstraint{<:OrderedSamplingAlgorithm}) → Tuple{Vector{Float64}, Vector{Float64}}
    get_quantiles!(lqs, uqs, x, c::SequentialSamplingConstraint{<:OrderedSamplingAlgorithm}) → Tuple{Vector{Float64}, Vector{Float64}}

Get the element-wise lower and upper quantiles for `x` respecting the quantiles given by the constraint `c`.
The in-place `get_quantiles!` writes the quantiles to pre-allocated vectors `lqs` and `uqs`.
"""
function get_quantiles end

"""
    get_quantiles(x, c::SequentialSamplingConstraint{<:OrderedSamplingAlgorithm}) → Tuple{Vector{Float64}, Vector{Float64}}
    get_quantiles!(lqs, uqs, x, c::SequentialSamplingConstraint{<:OrderedSamplingAlgorithm}) → Tuple{Vector{Float64}, Vector{Float64}}

Get the element-wise lower and upper quantiles for `x` respecting the quantiles given by the constraint `c`.
The in-place `get_quantiles!` writes the quantiles to pre-allocated vectors `lqs` and `uqs`.
"""
function get_quantiles! end

const SSC = SequentialSamplingConstraint{<:OrderedSamplingAlgorithm}

# Generic case takes care of vectors of uncertain values
function get_quantiles(x, c::SSC)
    lqs = quantile.(x, c.lq, c.n)
    uqs = quantile.(x, c.uq, c.n)
    return lqs, uqs
end

function get_quantiles!(lqs, uqs, x, c::SSC)
    length(lqs) == length(uqs) == length(x) || error("Lengths of `lqs`, `uqs` and `x` do not match.")
    lqs .= quantile.(x, c.lq, c.n)
    uqs .= quantile.(x, c.uq, c.n)
end

# Be specific when it comes to concrete types with differently named data fields.
function get_quantiles(x::UncertainIndexDataset, c::SSC)
    lqs = quantile.(x.indices, c.lq, c.n)
    uqs = quantile.(x.indices, c.uq, c.n)
    return lqs, uqs
end

function get_quantiles!(lqs, uqs, x::UncertainIndexDataset, c::SSC)
    length(lqs) == length(uqs) == length(x) || error("Lengths of `lqs`, `uqs` and `x` do not match.")
    lqs .= quantile.(x.indices, c.lq, c.n)
    uqs .= quantile.(x.indices, c.uq, c.n)
end

function get_quantiles(x::UncertainValueDataset, c::SSC)
    lqs = quantile.(x.values, c.lq, c.n)
    uqs = quantile.(x.values, c.uq, c.n)
    return lqs, uqs
end

function get_quantiles!(x::UncertainValueDataset, c::SSC)
    length(lqs) == length(uqs) == length(x) || error("Lengths of `lqs`, `uqs` and `x` do not match.")
    lqs .= quantile.(x.values, c.lq, c.n)
    uqs .= quantile.(x.values, c.uq, c.n)
end

"""
    truncated_supports(udata::AbstractUncertainValueDataset; quantiles = [0.001, 0.999])

Truncate the furnishing distribution of each uncertain value in the dataset
to the provided `quantiles` range. Returns a vector of Interval{Float64}, one
for each value.
"""
function truncated_supports(udata::AbstractUncertainValueDataset;
            quantiles = [0.001, 0.999])
    n_vals = length(udata)

    # Using the provided quantiles, find the parts of the supports of the furnishing
    # distributions from which we're going to sample.
    supports = Vector{Interval{Float64}}(undef, n_vals)

    for i = 1:n_vals
        lowerbound = quantile(udata[i], minimum(quantiles))
        upperbound = quantile(udata[i], maximum(quantiles))

        supports[i] = interval(lowerbound, upperbound)
    end

    return supports
end


text:
using Pkg
Pkg.activate(pwd())

# # Other Useful Packages
# ## Distributions.jl

using Distributions

D = Normal(2, 0.5)

#+

mean(D)
var(D)
quantile(D, 0.9)
pdf(D, 2)
cdf(D, 2)

#+

using StatsPlots

plot(
    plot(D; title = "pdf"),
    plot(D; func = cdf,  title = "cdf");
    legend = false,
    xlabel = "x",
    ylabel = "f(x)",
    ylims = (0,1),
    linewidth = 2,
    layout = (1,2),
    size = (800, 400),
)

#+

x = rand(Normal(2, 0.5), 10000); # generate 10000 random numbers from Normal(2, 0.5)
D = fit(Normal, x)

#+

histogram(x; normalize = :pdf, legend = false, opacity = 0.5)
plot!(D; linewidth = 2, xlabel = "x", ylabel = "pdf(x)")

# ### Exercise:
# Create a figure that shows the gamma distributions with the following parameters:
# `(2, 2)`, `(9, 0.5)`, `(7.5, 1)` and `(0.5, 1)`.
#
# **Hint:** to plot cumulative probability functions, use the Plots ability to
# plot functions.
#
# ---
# ### Solution:

Ds = Gamma.([2, 9, 7.5, 0.5], [2, 0.5, 1, 1])
labels = reshape(string.("Gamma", params.(Ds)), 1, :)

plot(Ds;
    xaxis = ("x", (0, 20)),
    yaxis = ("pdf(x)", (0, 0.5)),
    labels = labels,
    linewidth = 2,
    legend = :topright,
)

#+

plot(Ds;
    func = cdf,
    xaxis = ("x", (0, 20)),
    yaxis = ("cdf(x)", (0, 1.05)),
    labels = labels,
    linewidth = 2,
    legend = :bottomright,
)

#+

cdfs = [x -> cdf(D, x) for D in Ds]
plot(cdfs, 0, 20;
    xaxis = ("x", (0, 20)),
    yaxis = ("cdf(x)", (0, 1.05)),
    labels = labels,
    linewidth = 2,
    legend = :bottomright,
)

# ---
#
# ## BSON.jl

using BSON

BSON.bson(joinpath(pwd(), "lecture_05", "test2.bson"), a = [1+2im, 3+4im], b = "Hello, World!")

#+

data = Dict(:a => [1+2im, 3+4im], :b => "Hello, World!")
BSON.bson(joinpath(pwd(), "lecture_05", "test1.bson"), data)

#+

BSON.load(joinpath(pwd(), "lecture_05", "test1.bson"))
BSON.load(joinpath(pwd(), "lecture_05", "test2.bson"))

#+

using BSON: @save, @load

a = [1+2im, 3+4im];
b = "Hello, World!";

@save joinpath(pwd(), "lecture_05", "test.bson") a b # Same as above
@load joinpath(pwd(), "lecture_05", "test.bson") a b # Loads `a` and `b` back into the workspace

# ## ProgressMeter.jl

using ProgressMeter

@showprogress 1 "Computing..." for i in 1:50
    sleep(0.1)
end

#+

x, n = 1 , 10;
p = Progress(n);

for iter in 1:10
    x *= 2
    sleep(0.5)
    ProgressMeter.next!(p; showvalues = [(:iter, iter), (:x, x)])
end

# ## BenchmarkTools.jl

using BenchmarkTools

@benchmark sin(x) setup=(x=rand())

#+

A = rand(3,3);
@btime inv($A);


text:
using Plots, LaTeXStrings; pyplot()

f2(x) = (x<0 ? x+1 : 1-x)*(abs(x)<1 ? 1 : 0)
a, b = -1.5, 1.5
delta = 0.01

F(x) = sum([f2(u)*delta for u in a:delta:x])

xGrid = a:delta:b
y = [F(u) for u in xGrid]
plot(xGrid, y, c=:blue, xlims=(a,b), ylims=(0,1),
	xlabel=L"x", ylabel=L"F(x)", legend=:none)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


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

using Markdown
using InteractiveUtils

# ╔═╡ 79789782-8d97-11eb-112c-fd8c2c294527
using Pkg, DrWatson

# ╔═╡ dd81850e-8d97-11eb-0cff-8bf4f38f0660
begin
	#@quickactivate "SR2StanPluto"
	using StanSample
	using StatisticalRethinking
	using BrowseTables, Tables
end

# ╔═╡ 3cb90c0e-8db1-11eb-324a-0f66762c31f6
md" ## HTMLTable.jl"

# ╔═╡ 24856618-d93c-4467-b04b-b952117bfed6
begin
	# make example table, but any table that supports Tables.jl will work
	table = Tables.columntable(collect(i == 5 ?
			(a = missing, b = "string", c = nothing) :
			(a = i, b = Float64(i), c = 'a'-1+i) for i in 1:10))
end;

# ╔═╡ 2e1e3519-81ad-4124-9e2a-c3128465be46
HTMLTable(table) # show HTML table using Julia's display system

# ╔═╡ 0ad949c1-8b2e-439e-a5e4-3dcb168dcc1e
begin
    df = CSV.read(sr_datadir("chimpanzees.csv"), DataFrame)
    df[!, "treatment"] = 1 .+ df.prosoc_left + 2 * df.condition
end;

# ╔═╡ 1db203ac-fdcf-4ece-8e94-52bf1c32f86a
stan11_4 = "
data{
	int N;
	int K;
    int pulled_left[N];
    int treatment[N];
    int actor[N];
}
parameters{
    vector[7] a;
    vector[4] b;
}
model{
    vector[N] p;
    b ~ normal( 0 , 0.5 );
    a ~ normal( 0 , 1.5 );
    for ( i in 1:N ) {
        p[i] = a[actor[i]] + b[treatment[i]];
        p[i] = inv_logit(p[i]);
    }
    pulled_left ~ binomial( 1 , p );
}
";

# ╔═╡ 3864cff3-f9c1-42a0-96be-a04fc9d3f325
begin
	m11_4s = SampleModel("m11.4s", stan11_4)
	data = (N = size(df, 1), K = length(unique(df.treatment)),
		pulled_left = df.pulled_left, actor=df.actor,
		treatment = df.treatment)
	rc11_4s = stan_sample(m11_4s; data)
end;

# ╔═╡ d92eb94f-3003-4ea5-aeb1-ae6ee8978b6a
if success(rc11_4s)
	post11_4s_df = read_samples(m11_4s, :dataframe)
	PRECIS(post11_4s_df[:, 1:11])
end

# ╔═╡ 59879eb7-f3aa-4b90-a22e-653ccd9d5616
HTMLTable(post11_4s_df[1:10, 1:11])

# ╔═╡ f1dc6158-aa8d-41a1-b292-fb5ff917b5fa
function precis_df(df::DataFrame; digits = 4, depth = Inf, alpha = 0.11)
    d = DataFrame()
    d.param = names(df)
    d.mean = mean.(eachcol(df))
    d.std = std.(eachcol(df))
    d[:, "5.5%"] = quantile.(eachcol(df), alpha/2)
    d[:, "50%"] = quantile.(eachcol(df), 0.5)
    d[:, "94.5%"] = quantile.(eachcol(df), 1 - alpha/2)
    u = StatisticalRethinking.unicode_histogram.(eachcol(df), min(size(df, 1), 12))

	#=
    d.histogram = StatisticalRethinking.unicode_histogram.(eachcol(df),
		min(size(df, 1), 12))
	for n in names(df)
		if eltype(df[!,n]) <: Number
			df[!,n] = round.(df[!,n], digits=2)
		end
	end
	=#

    for col in ["mean", "std", "5.5%", "50%", "94.5%"]
        d[:, col] .= round.(d[:, col], digits = digits)
    end

    d
end

# ╔═╡ ec6df58e-3baf-4db0-8c5f-4189988e2c4b
begin
	summary_df = precis_df(post11_4s_df[:, 1:11])
	HTMLTable(summary_df)
end

# ╔═╡ ae03e891-5f44-4526-98df-1de2d113173d
if success(rc11_4s)
	sdf = read_summary(m11_4s)
	HTMLTable(sdf[8:end,:]; caption="m11.4s")
end

# ╔═╡ 3005a46a-571c-453c-b55d-e9eab17c2b87
md"# StatisticalRethinking
| Chapter | Clip | Topic |
|:---|---|---|
|Running a Stan Language program|||
|Introduction to the Stan Language      | Intros| StanQuap usage |
||Loglik|PSIS|
"

# ╔═╡ Cell order:
# ╠═3cb90c0e-8db1-11eb-324a-0f66762c31f6
# ╠═79789782-8d97-11eb-112c-fd8c2c294527
# ╠═dd81850e-8d97-11eb-0cff-8bf4f38f0660
# ╠═24856618-d93c-4467-b04b-b952117bfed6
# ╠═2e1e3519-81ad-4124-9e2a-c3128465be46
# ╠═0ad949c1-8b2e-439e-a5e4-3dcb168dcc1e
# ╠═1db203ac-fdcf-4ece-8e94-52bf1c32f86a
# ╠═3864cff3-f9c1-42a0-96be-a04fc9d3f325
# ╠═d92eb94f-3003-4ea5-aeb1-ae6ee8978b6a
# ╠═59879eb7-f3aa-4b90-a22e-653ccd9d5616
# ╠═f1dc6158-aa8d-41a1-b292-fb5ff917b5fa
# ╠═ec6df58e-3baf-4db0-8c5f-4189988e2c4b
# ╠═ae03e891-5f44-4526-98df-1de2d113173d
# ╠═3005a46a-571c-453c-b55d-e9eab17c2b87


text:
function solvenjl(κ, coef, energ, j)

    @info "Solving the NJL equations..."

    constant = 4*Nc*Nf*G*ħ^3
    function f!(F, x)
        for k in 1:κ
            F[k] = constant*sum([coef[k, kprime]*(2*j[kprime]+1)*x[kprime]/√(x[kprime]^2+energ[kprime]^2) for kprime in 1:κ]) - x[k]
        end
    end

    # The quark masses are initialized as having 350 [MeV]
    initial_x = 350.0*ones(κ)
    result = nlsolve(f!, initial_x)
    mass = result.zero

    @info "Done!"
    @debug "Masses: " mass

    return mass

end

max_stars_repo_name: avramescudana/QuarkDroplets

max_issues_repo_name: avramescudana/QuarkDroplets

max_forks_repo_name: avramescudana/QuarkDroplets


text:
include("match.jl")

"""
Fallback MIP formulation for only 2-cycles and 2-3-cycles.
"""
function cycle_match(weights::Array{Float64,2};
                           graph::DiGraph=DiGraph(),
                           verbose::Int64 = 0,
                           max_cycle_length::Int64 = 2,
                           max_chain_length::Float64 = Inf64,
                           solver::String = "GLPK",
                           method::String = "fallback", #
                           ndds::Array{Int64, 1} = zeros(Int64, 0))

  if nv(graph) == 0
      graph = generate_graph(weights)
  end
  m, x, in = basic_MIP(weights, graph, ndds)
  n = nv(graph)
  if max_cycle_length == 2
      @constraint(m,
                  c5[i=1:n, j=intersect(outneighbors(graph, i), inneighbors(graph, i))],
                  x[Edge((i,j))] == x[Edge((j,i))])
      @constraint(m,
                  c6[i=1:n, j=setdiff(outneighbors(graph, i), inneighbors(graph, i))],
                  x[Edge((i,j))] == 0
                  )
  elseif max_cycle_length == 3
      @constraint(m,
                  c5[i=1:n,
                     j=outneighbors(graph, i),
                     k=intersect(outneighbors(graph, j), inneighbors(graph, i))],
                  x[Edge((k,i))] >= x[Edge((i,j))] + x[Edge((j,k))] - 1
                 )
     @constraint(m,
                 c6[i=1:n, j=outneighbors(graph, i), k=setdiff(outneighbors(graph, j), inneighbors(graph, i))],
                 Int(k == i) >= x[Edge((i,j))] + x[Edge((j,k))] - 1
                )
  else
      error("No fallback implemented for this setting")
  end
  status = optimize!(m)
  match_edges = JuMP.value.(x)
  match_vertices = JuMP.value.(in)
  value = JuMP.objective_value(m)

  return  (match_vertices, match_edges, value)
end


text:
module ColorData

using ImageCore, NamedDims, ImageAxes, AxisIndices, NamedIndicesMeta

import ImageAxes: colordim

export
    is_color,
    # @defdim output
    colordim,
    has_colordim,
    color_axis,
    color_axis_type,
    color_keys,
    color_indices,
    to_color_format,
    ncolor,
    color_format,
    # reexport
    permuteddimsview,
    channelview

Base.@pure is_color(x::Symbol) = x === :color

@defdim color is_color

function ImageCore.permuteddimsview(nda::NamedDimsArray, perm)
    return NamedDimsArray(
        permuteddimsview(parent(nda), perm),
        NamedDims.permute_dimnames(dimnames(nda), perm)
    )
end

function ImageCore.permuteddimsview(x::AbstractAxisIndices, perm)
    return AxisIndices.unsafe_reconstruct(
        x,
        permuteddimsview(parent(x), perm),
        AxisIndices.permute_axes(x, perm)
    )
end

function ImageCore.channelview(A::NamedDimsArray{L}) where {L}
    return _channelview(L, channelview(parent(A)))
end
_channelview(n::Tuple{Vararg{Symbol,N}}, a::AbstractArray{T,M}) where {T,N,M} = NamedDimsArray{(:color, n...,)}(a)
_channelview(n::Tuple{Vararg{Symbol,N}}, a::AbstractArray{T,N}) where {T,N} = NamedDimsArray{n}(a)

ImageCore.channelview(A::AbstractAxisIndices) = _channelview(A, channelview(parent(A)), axes(A))
function _channelview(A::AbstractAxisIndices, a::AbstractArray{T,N}, axs::Tuple{Vararg{<:AbstractAxis,N}}) where {T,N}
    return AxisIndices.similar_type(A, typeof(a), typeof(axs))(a, axs)
end

function _channelview(A::AbstractAxisIndices, a::AbstractArray{T,N}, axs::Tuple{Vararg{<:AbstractAxis,M}}) where {T,N,M}
    return _channelview(A, a, (AxisIndices.as_axis(A, axes(a, 1)), axs...,))
end

end