text:
abstract type MatpowerRTEROPFSimpleInput <: AbstractInput end
using MAT

"""
    read_input(input_type::T, instance_path::String) where T<:Type{MatpowerInput}

Read instance in `instance_path` depending on `input_type`.\n
Return a structure OPFProblems.
"""
function read_input(input_type::T, instance_path::String) where T<:Type{MatpowerRTEROPFSimpleInput}
  file = matopen(instance_path)
  data = read(file, "mpc")
  close(file)
  # DataStructure and Gridstructure data:
  bus = SortedDict{String, SortedDict{String, Any}}()
  link = SortedDict{Link, SortedDict{String, Any}}()

  bus_id_line=SortedDict{Int, Int}()
  bus_id_name=SortedDict{Int, String}()

  baseMVA = data["baseMVA"]

  ## Building bus load and shunt information
  data_bus = data["bus"]
  bustype = data_bus[1:size(data_bus,1), 2]       # weather bus is 1:"PQ" (generators are not accounted for) or 2:"PV"
  for i in 1:size(data_bus,1)
      busname = bus_name(i)
      ## Adding MatpowerVolt structure (for each bus)
      if !haskey(bus, busname)
          bus[busname] = SortedDict{String, Any}()
      end
      bus[busname][volt_name()] = MatpowerVolt(busname, i, data_bus[i,13], data_bus[i,12])

      ## Matpower Load
      load = data_bus[i,3] + im*data_bus[i,4]
      if load != 0
        bus[busname]["Load"] = MatpowerLoad(load)
      end

      ## Matpower Shunt
      shunt = data_bus[i,5] + im*data_bus[i,6]
      if shunt != 0
        bus[busname]["Shunt"] = MatpowerROPFShunt(shunt)
      end

      bus_id_line[data_bus[i,1]] = i
      bus_id_name[data_bus[i,1]] = busname
  end

  ## Adding bus generator information
  gen2bus = SortedDict{Int, Int}()
  line2busgen = SortedDict{Int, Tuple{Int, Int}}()

  network, instance = split(instance_path, "\\")[end-1:end]
  gen_data =  data["gen"]
  instance = split(instance_path, "\\")[end][1:end-2]
  dict_Pinput = Dict{String, Tuple{Float64, Float64, Float64}}()
    genind, prevgen = 0, 0
  for i=1:size(gen_data,1)
    gen2bus[i] = bus_id_line[gen_data[i,1]]
    busname = bus_id_name[gen_data[i,1]]

    Pinput = 0
    S_min = S_max = 0
    if gen_data[i, 1] == prevgen
      genind += 1
      S_min = bus[busname]["Gen"].power_min
      S_max = bus[busname]["Gen"].power_max
      Pinput = bus[busname]["Gen"].active_power_input
    else
      prevgen = gen_data[i, 1]
      genind = 1
    end
    line2busgen[i] = (gen_data[i, 1], genind)

    if gen_data[i, 8] > 0 #generator is on
      S_min += gen_data[i,10] + im*gen_data[i,5] #Smin = Pmin + i Qmin
      S_max += gen_data[i,9] + im*gen_data[i,4] #Smax = Pmax + i Qmax
      Pinput += gen_data[i,2]
    end
    if bustype[bus_id_line[gen_data[i,1]]] == 3.0
      # println("reference bus : $i")
      network, instance = split(instance_path, "\\")[end-1:end]
      f = open("D:\\repo\\Test_Balthazar\\data_ROPF\\$(network)-$(instance[1:end-4])_refbus.csv", "w")
      write(f, "Gen_Sgen_$(bus_id_line[gen_data[i, 1]])_Re\n")
      close(f)
    end
    dict_Pinput["Gen_Sgen_$(bus_id_line[gen_data[i, 1]])_Re"] = (Pinput, real(S_min), real(S_max))
    bus[busname]["Gen"] = RTEMatpowerGenerator(S_min, S_max, Pinput, true)
  end
  network, instance = split(instance_path, "\\")[end-1:end]
  f = open("D:\\repo\\Test_Balthazar\\data_ROPF\\$(network)-$(instance[1:end-4])_generation.csv", "w")
  for (gen_name, tuple) in dict_Pinput
    Pinput = tuple[1]
    Pmin = tuple[2]
    Pmax = tuple[3]
    write(f, "$(gen_name) ; $Pinput ; ; $(Pmin) ; $(Pmax) \n")
  end
  close(f)

  ## building link information
  data_branch = data["branch"]
  for i=1:size(data_branch,1)
    linkname = Link(bus_id_name[data_branch[i,1]], bus_id_name[data_branch[i,2]])
    if data_branch[i, 11] == 0
      println("$link $(linkname.orig)⟶$(linkname.dest) breaker out of service !")
    else
      rs, xs, bc = data_branch[i,3:5]
      Smax = data_branch[i,6]
      τ, θ = data_branch[i, 9:10]
      if !haskey(link, linkname)
        link[linkname] = SortedDict{String, Any}()
      end

      nb_elem = length(link[linkname])
      link[linkname]["LinkMP_$(nb_elem+1)"] = MatpowerLine_π(baseMVA, rs, xs, bc, τ, θ)
    end
  end

  bus_line_id = SortedDict([val=>key for (key,val) in bus_id_line])

  ## Removing generators from PQ buses
  for (busname, bus_elems) in bus
    i = bus_elems["Volt"].busid
    if bustype[i] == 1
      for (elemid, elem) in bus_elems
        if typeof(elem) == MatpowerGenerator
          delete!(bus_elems, elemid)
        end
      end
    end
  end

  ## Adding null generator to reference bus
  refind = findall(bustype .== 3)
  length(refind) == 1 || warn("/! refind = $refind")
  bus["BUS_$(refind[1])"]["Gen_reference"] = MatpowerGenerator(0,0,3,[0 0 0],true)

  ds = DataSource(bus, link)

  ## Building GridStructure
  node_linksin = node_linksout = SortedDict{String, SortedSet{Link}}()
  node_vars = SortedDict{String, SortedDict{String, Variable}}()
  link_vars = SortedDict{Link, SortedDict{String, Variable}}()
  gs = GridStructure(basecase_scenario_name(), node_linksin, node_linksout)

  node_formulations = SortedDict{String, SortedDict{String, Symbol}}()
  link_formulations = SortedDict{Link, SortedDict{String, Symbol}}()
  mp = MathematicalProgramming(node_formulations, link_formulations, node_vars, link_vars)
  return OPFProblems(basecase_scenario_name()=>Scenario(ds, gs, mp))
end


text:
using FractionalCalculus
using Plots, LaTeXStrings

s="\$(x^4)'''=24x\$"

tspan=collect(0:0.01:6)

result=fracdiff(cos, 3, 6, 0.01, RLDiff_Matrix())

plot(tspan, result, title=s, legend=:bottomright)

ext: jl


text:
using StochasticDiffEq, DiffEqProblemLibrary, Random
Random.seed!(200)
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
prob = DiffEqProblemLibrary.SDEProblemLibrary.oval2ModelExample(largeFluctuations=true,useBigs=false)

Random.seed!(200)
@time for i in 1:10
    sol = solve(prob,SRIW1(),dt=(1/2)^(18),qmax=1.125,
        saveat=0.1,maxiters=1e7,abstol=1e-5,reltol=1e-3)
end

Random.seed!(200)
@time for i in 1:10
    @show i
    sol = solve(prob,ImplicitEM(),dt=1/60000)
end

Random.seed!(200)
@time for i in 1:10
    @show i
    sol = solve(prob,ImplicitRKMil(),dt=1/50000)
end

Random.seed!(200)
@time for i in 1:10
    sol = solve(prob,SOSRI(),dt=(1/2)^(18),qmax=1.125,
          saveat=0.1,maxiters=1e7,abstol=1e-4,reltol=1e-2)
end

Random.seed!(200)
@time for i in 1:10
    sol = solve(prob,SOSRI2(),dt=(1/2)^(18),qmax=1.125,
          saveat=0.1,maxiters=1e7,abstol=1e-4,reltol=1e-4)
end

Random.seed!(200)
@time for i in 1:10
    sol = solve(prob,SOSRI2(),dt=(1/2)^(18),qmax=1.125,
          saveat=0.1,maxiters=1e7,abstol=1e-5,reltol=1e-3)
end

Random.seed!(200)
@time for i in 1:10
    sol = solve(prob,SOSRI(),dt=(1/2)^(18),qmax=1.125,
          saveat=0.1,maxiters=1e7,abstol=1e-3,reltol=1e-2)
end

Random.seed!(200)
@time for i in 1:10
    sol = solve(prob,SOSRI(),dt=(1/2)^(18),qmax=1.125,
          saveat=0.1,maxiters=1e7,abstol=1e-4,reltol=1e-4)
end

Random.seed!(200)
@time for i in 1:10
    sol = solve(prob,SOSRI(),dt=(1/2)^(18),qmax=1.125,
          saveat=0.1,maxiters=1e7,abstol=1e-2,reltol=1e-2)
end

Random.seed!(200)
@time for i in 1:10
    sol = solve(prob,SOSRI(),dt=(1/2)^(18),qmax=1.125,
          saveat=0.1,maxiters=1e7,abstol=1e-5,reltol=1e-3)
end

Random.seed!(200)
@time for i in 1:10
    sol = solve(prob,SOSRI(),dt=(1/2)^(18),qmax=1.125,
          saveat=0.1,maxiters=1e7,abstol=1e-2,reltol=1e-1)
end

Random.seed!(200)
@time for i in 1:10
    sol = solve(prob,SOSRI2(),dt=(1/2)^(18),qmax=1.125,
          saveat=0.1,maxiters=1e7,abstol=1e-4,reltol=1e-1)
end

Random.seed!(200)
@time for i in 1:10
    @show i
    sol = solve(prob,ImplicitEM(),dt=1/50000)
end

Random.seed!(200)
@time for i in 1:10
    @show i
    sol = solve(prob,ImplicitRKMil(),dt=1/40000)
end

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
using MixedModelsSim, MixedModels, StableRNGs
using DataFrames, Tables
using Statistics
using Test

subj_n = 40
item_n = 20
subj_btwn = Dict(:age => ["O", "Y"],
				 :pet => ["cat", "dog"])
item_btwn = Dict(:cond => ["A", "B"])
both_win = Dict(:time => ["morning", "evening"])
dat = simdat_crossed(StableRNG(42), subj_n, item_n,
					 subj_btwn = subj_btwn,
					 item_btwn = item_btwn,
                     both_win = both_win)
form = @formula(dv ~ 1 + age + pet + cond + time  + (1|subj) + (1|item));
cont = Dict(nm => HelmertCoding() for nm in (:age, :pet, :cond, :time));
fm1 = fit(MixedModel, form, dat; contrasts=cont, REML=false, progress=false);
sim = parametricbootstrap(StableRNG(42), 10, fm1; use_threads=false, hide_progress=true);

@testset "power_table" begin
    pt = power_table(sim)

    @test Tables.isrowtable(pt)

    ptdf = sort!(DataFrame(pt), :coefname)

    @test nrow(ptdf) == 5
    @test ncol(ptdf) == 2
    @test names(ptdf) == ["coefname", "power"]
    @test ptdf.power == [0.4, 0.2, 0.0, 0.3, 0.1]
end


text:
macro import_compcodomains_huge()
    :(
        using ComplexVisual:
            cv_destroy,
            CV_CompareCodomainParts,
            cv_get_trafo1, cv_get_trafo2,
            cv_get_can_codomain1, cv_get_can_codomain2,
            cv_get_cc_can_codomain1, cv_get_cc_can_codomain2,
            cv_add,
            CV_CodomainsPosLayout, CV_ParamWithSlider,
            cv_do_comp_codomains_layout,
            cv_setup_comp_codomains_axis, cv_setup_comp_codomains_painters,
            cv_setup_comp_codomains_border, cv_setup_comp_codomains_scene,
            cv_setup_comp_slider_cb, cv_setup_comp_sliders,
            cv_scene_comp_codomains_std
    )
end

"""
Saves two codomains (and their contexts) together with two trafos in (the)
layout.
"""
struct CV_CompareCodomainParts{parentT<:CV_Abstract2DLayout,
        trafo1T, trafo2T, codomain1T, codomain2T,
        ccc1T, ccc2T} <: CV_2DLayoutWrapper   # {{{
    parent_layout   :: parentT
    trafo1          :: trafo1T
    trafo2          :: trafo2T
    can_codomain1   :: codomain1T
    can_codomain2   :: codomain2T
    cc_can_codomain1:: ccc1T
    cc_can_codomain2:: ccc2T
end

@layout_composition_getter(trafo1,           CV_CompareCodomainParts)
@layout_composition_getter(trafo2,           CV_CompareCodomainParts)
@layout_composition_getter(can_codomain1,    CV_CompareCodomainParts)
@layout_composition_getter(can_codomain2,    CV_CompareCodomainParts)
@layout_composition_getter(cc_can_codomain1, CV_CompareCodomainParts)
@layout_composition_getter(cc_can_codomain2, CV_CompareCodomainParts)

function cv_destroy(l::CV_CompareCodomainParts)
    cv_destroy(l.cc_can_codomain1)
    cv_destroy(l.cc_can_codomain2)
    cv_destroy(l.can_codomain1)
    cv_destroy(l.can_codomain2)
    cv_destroy(l.parent_layout)
    return nothing
end

"""
convenience function to combine domain and codomain (and their contexts)
together with a trafo to (the) layout.
"""
function cv_add(layout::CV_Abstract2DLayout,
        trafo1, trafo2,
        can_codomain1::codomain1T=CV_Math2DCanvas(-2.0+2.0im, 2.0-2.0im, 50),
        can_codomain2::codomain2T=CV_Math2DCanvas(-2.0+2.0im, 2.0-2.0im, 50),
        cc_can_codomain1::CV_CanvasContext=cv_create_context(can_codomain1),
        cc_can_codomain2::CV_CanvasContext=cv_create_context(can_codomain2)
        ) where {codomain1T<:CV_Canvas, codomain2T<:CV_Canvas}
    return CV_CompareCodomainParts(
        layout, trafo1, trafo2, can_codomain1, can_codomain2,
        cc_can_codomain1, cc_can_codomain2)
end
# }}}

"""
A layout with both codomains positioned.
"""
struct CV_CodomainsPosLayout{parentT<:CV_Abstract2DLayout,
        can1T, dcb1T, style1T, can2T, dcb2T, style2T} <: CV_2DLayoutWrapper # {{{
    parent_layout       :: parentT
    can_codomain1_l     :: CV_2DLayoutPosition{can1T, dcb1T, style1T}
    can_codomain2_l     :: CV_2DLayoutPosition{can2T, dcb2T, style2T}
end

@layout_composition_getter(can_codomain1_l,    CV_CodomainsPosLayout)
@layout_composition_getter(can_codomain2_l,    CV_CodomainsPosLayout)
# }}}

"""
All informations for a dynamic parameter (represented by a slider)
"""
struct CV_ParamWithSlider{posT, N, funcT, setupT} # {{{
    slider_pos             :: posT
    value_min              :: Real
    value_max              :: Real
    rulers                 :: NTuple{N, CV_Ruler}
    set_slider_value_func  :: funcT
    gap_below              :: Integer
    setup_cb               :: setupT
end

function CV_ParamWithSlider(;
        slider_pos=CV_TranslateByOffset(Float64),
        value_min=0.0, value_max=1.0,
        rulers=(CV_Ruler(cv_format_ticks("%.1f", 0.0:0.5:1.0...)),),
        set_slider_value_func=nothing,
        gap_below=40, setup_cb=cv_setup_comp_slider_cb)
    return CV_ParamWithSlider(slider_pos, value_min, value_max, rulers,
        set_slider_value_func, gap_below, setup_cb)
end # }}}

"""
Do a left-right layout of the two codomains with a gap between.
"""
function cv_do_comp_codomains_layout(layout::CV_Abstract2DLayout,
        gap_domains::Integer=50) # {{{
    can_codomain1 = cv_get_can_codomain1(layout)
    can_codomain2 = cv_get_can_codomain2(layout)
    can_codomain1_l = cv_add_canvas!(layout, can_codomain1,
        cv_anchor(can_codomain1, :northwest), (0, 0))
    can_codomain2_l = cv_add_canvas!(layout, can_codomain2,
        cv_anchor(can_codomain2, :west),
        cv_translate(cv_anchor(can_codomain1_l, :east), gap_domains, 0))

    return CV_CodomainsPosLayout(layout, can_codomain1_l, can_codomain2_l)
end # }}}

"""
setup rulers for codomains.
"""
function cv_setup_comp_codomains_axis(setup::CV_SceneSetupChain,
        codomain1_re_rulers::NTuple{A, CV_Ruler},
        codomain1_im_rulers::NTuple{B, CV_Ruler}=codomain1_re_rulers,
        codomain2_re_rulers::NTuple{C, CV_Ruler}=codomain1_re_rulers,
        codomain2_im_rulers::NTuple{D, CV_Ruler}=codomain1_im_rulers
        ) where {A, B, C, D}   # {{{
    layout = setup.layout
    can_codomain1_l = cv_get_can_codomain1_l(layout)
    can_codomain2_l = cv_get_can_codomain2_l(layout)

    codomain1_re_axis = cv_ticks_labels(layout, can_codomain1_l,
        cv_south, codomain1_re_rulers)
    codomain1_im_axis = cv_ticks_labels(layout, can_codomain1_l,
        cv_west, codomain1_im_rulers)
    codomain2_re_axis = cv_ticks_labels(layout, can_codomain2_l,
        cv_south, codomain2_re_rulers)
    codomain2_im_axis = cv_ticks_labels(layout, can_codomain2_l,
        cv_west, codomain2_im_rulers)

    draw_once_func = future_layout -> begin
        cc_can_layout = cv_get_cc_can_layout(future_layout)
        codomain1_re_axis(cc_can_layout)
        codomain1_im_axis(cc_can_layout)
        codomain2_re_axis(cc_can_layout)
        codomain2_im_axis(cc_can_layout)
        return nothing
    end
    return cv_combine(setup; draw_once_func)
end # }}}

function cv_setup_comp_codomains_painters(setup::CV_SceneSetupChain,
        painter1, painter2) # {{{

    layout = setup.layout

    cc_can_codomain1 = cv_get_cc_can_codomain1(layout)
    cc_can_codomain2 = cv_get_cc_can_codomain2(layout)

    redraw_func = layout -> begin
        if painter1 !== nothing
            cv_clear_cache(painter1)
            cv_paint(cc_can_codomain1, painter1)
        end
        if painter2 !== nothing
            cv_clear_cache(painter2)
            cv_paint(cc_can_codomain2, painter2)
        end
        return nothing
    end

    return cv_combine(setup; redraw_func)
end  # }}}

"""
creates borders for the two codomains.
"""
function cv_setup_comp_codomains_border(setup::CV_SceneSetupChain;
        width::Integer=2, style=cv_black) # {{{
    layout = setup.layout
    codomain1_border = cv_border(layout, cv_get_can_codomain1_l(layout), width;
        style)
    codomain2_border = cv_border(layout, cv_get_can_codomain2_l(layout), width;
        style)

    draw_once_func = future_layout -> begin
        cc_can_layout = cv_get_cc_can_layout(future_layout)
        codomain1_border(cc_can_layout)
        codomain2_border(cc_can_layout)
        return nothing
    end

    return cv_combine(setup; draw_once_func)
end # }}}

"""
create scene for codomains with sliders.
"""
function cv_setup_comp_codomains_scene(setup::CV_SceneSetupChain)  # {{{
    layout = setup.layout
    can_codomain1_l = cv_get_can_codomain1_l(layout)
    can_codomain2_l = cv_get_can_codomain2_l(layout)

    redraw_func = future_layout -> begin
        cc_can_layout = cv_get_cc_can_layout(future_layout)
        can_codomain1_l(cc_can_layout)
        can_codomain2_l(cc_can_layout)
        return nothing
    end

    setup = cv_combine(setup; redraw_func)
    return cv_setup_2dminimal_scene(setup)
end # }}}

"""
creates default slider with `width` centered at `mid` and starting at `bottom`
with a value-mark painter.
"""
function cv_setup_comp_slider_cb(setup, param::CV_ParamWithSlider,
        width::Int32, mid::Int32, bottom::Int32) # {{{

    layout = setup.layout

    slider_data = cv_create_hslider(width, 20,
        param.value_min, param.value_max, param.rulers; attach=cv_north,
        decoraction_with_layout_and_position_callback=(inner_layout, pos) ->
            cv_border(inner_layout, pos, 1))

    slider_container = slider_data.slider_container
    slider_l = cv_add_canvas!(layout,
        slider_container,
        cv_anchor(slider_container, :slider_south, slider_container, :south),
        (mid, bottom - Int32(param.gap_below)))

    if param.set_slider_value_func === nothing
        set_slider_value_func = z -> begin
            slider_pos.value = max(min(z, value_max), value_min)
            return CV_Response(; redraw_flag=true)
        end
    else
        set_slider_value_func = param.set_slider_value_func
    end

    bg_painter = cv_color(.8, .8, .8) ↦ CV_FillPainter()
    mark_painter = (cv_op_source → cv_color(0, 0, 1) → cv_linewidth(2)) ↦
        CV_ValueMarkPainter(param.slider_pos, 0.0,
            imag(slider_container.can_slider.corner_ul), false)

    setup = cv_setup_hslider(setup, slider_data, slider_l,
        bg_painter → mark_painter, set_slider_value_func)
    return setup, slider_l.rectangle.bottom
end  # }}}

"""
setup sliders for parameters
"""
function cv_setup_comp_sliders(setup,
        params::NTuple{N, CV_ParamWithSlider}) where {N}  # {{{

    layout = setup.layout
    can_codomain1_l = cv_get_can_codomain1_l(layout)
    can_codomain2_l = cv_get_can_codomain2_l(layout)

    r1, r2 = can_codomain1_l.rectangle, can_codomain2_l.rectangle
    width, mid = r2.right - r1.left, cv_half(r2.right + r1.left)
    bottom = min(r1.bottom, r2.bottom)

    for param in params[N:-1:1]
        setup, bottom = param.setup_cb(setup, param, width, mid, bottom)
    end
    return setup
end # }}}

"""
creates "standard" scene with two codomains in left-right layout and a
slider (north).
"""
function cv_scene_comp_codomains_std(
        params::NTuple{N, CV_ParamWithSlider}, trafo1, trafo2,
        codomain1, codomain2;
        gap_domains::Integer=80, padding=30,
        codomain1_re_rulers::NTuple{A, CV_Ruler},
        codomain1_im_rulers::NTuple{B, CV_Ruler}=codomain1_re_rulers,
        codomain2_re_rulers::NTuple{C, CV_Ruler}=codomain1_re_rulers,
        codomain2_im_rulers::NTuple{D, CV_Ruler}=codomain1_im_rulers,
        painter1 = CV_PortraitPainter(trafo1),
        painter2 = CV_PortraitPainter(trafo2)
        ) where {A, B, C, D, N} # {{{

    layout = CV_2DLayout()
    layout = cv_add(layout, trafo1, trafo2, codomain1, codomain2)
    layout = cv_do_comp_codomains_layout(layout, gap_domains)

    setup = CV_MinimalSetupChain(layout)
    setup = cv_setup_comp_codomains_painters(setup, painter1, painter2)
    setup = cv_setup_comp_codomains_axis(setup,
        codomain1_re_rulers, codomain1_im_rulers,
        codomain2_re_rulers, codomain2_im_rulers)
    setup = cv_setup_comp_codomains_border(setup)
    setup = cv_setup_comp_sliders(setup, params)
    padding > 0 && cv_add_padding!(setup.layout, padding)

    setup = cv_setup_comp_codomains_scene(setup)
    return setup.layout
end # }}}

# vim:syn=julia:cc=79:fdm=marker:sw=4:ts=4:


text:
using Test
using Miletus
using Dates, Random

include("daycounts.jl")
include("ts.jl")
include("currency.jl")
include("misc.jl")
include("options.jl")
include("ivol.jl")
include("montecarlo.jl")
include("greeks.jl")


text:
# Computing deviation in a variety of ways

## count the number of equal/non-equal pairs

"""
    counteq(a, b)

Count the number of indices at which the elements of the arrays
`a` and `b` are equal.
"""
function counteq(a::AbstractArray, b::AbstractArray)
    n = length(a)
    length(b) == n || throw(DimensionMismatch("Inconsistent lengths."))
    c = 0
    for i = 1:n
        @inbounds if a[i] == b[i]
            c += 1
        end
    end
    return c
end

"""
    countne(a, b)

Count the number of indices at which the elements of the arrays
`a` and `b` are not equal.
"""
function countne(a::AbstractArray, b::AbstractArray)
    n = length(a)
    length(b) == n || throw(DimensionMismatch("Inconsistent lengths."))
    c = 0
    for i = 1:n
        @inbounds if a[i] != b[i]
            c += 1
        end
    end
    return c
end

"""
    sqL2dist(a, b)

Compute the squared L2 distance between two arrays: ``\\sum_{i=1}^n |a_i - b_i|^2``.
Efficient equivalent of `sumabs2(a - b)`.
"""
function sqL2dist(a::AbstractArray{T}, b::AbstractArray{T}) where T<:Number
    n = length(a)
    length(b) == n || throw(DimensionMismatch("Input dimension mismatch"))
    r = 0.0
    for i = 1:n
        @inbounds r += abs2(a[i] - b[i])
    end
    return r
end

# L2 distance
"""
    L2dist(a, b)

Compute the L2 distance between two arrays: ``\\sqrt{\\sum_{i=1}^n |a_i - b_i|^2}``.
Efficient equivalent of `sqrt(sumabs2(a - b))`.
"""
L2dist(a::AbstractArray{T}, b::AbstractArray{T}) where {T<:Number} = sqrt(sqL2dist(a, b))

# L1 distance
"""
    L1dist(a, b)

Compute the L1 distance between two arrays: ``\\sum_{i=1}^n |a_i - b_i|``.
Efficient equivalent of `sum(abs, a - b)`.
"""
function L1dist(a::AbstractArray{T}, b::AbstractArray{T}) where T<:Number
    n = length(a)
    length(b) == n || throw(DimensionMismatch("Input dimension mismatch"))
    r = 0.0
    for i = 1:n
        @inbounds r += abs(a[i] - b[i])
    end
    return r
end

# Linf distance
"""
    Linfdist(a, b)

Compute the L∞ distance, also called the Chebyshev distance, between
two arrays: ``\\max_{i\\in1:n} |a_i - b_i|``.
Efficient equivalent of `maxabs(a - b)`.
"""
function Linfdist(a::AbstractArray{T}, b::AbstractArray{T}) where T<:Number
    n = length(a)
    length(b) == n || throw(DimensionMismatch("Input dimension mismatch"))
    r = 0.0
    for i = 1:n
        @inbounds v = abs(a[i] - b[i])
        if r < v
            r = v
        end
    end
    return r
end

# Generalized KL-divergence
"""
    gkldiv(a, b)

Compute the generalized Kullback-Leibler divergence between two arrays:
``\\sum_{i=1}^n (a_i \\log(a_i/b_i) - a_i + b_i)``.
Efficient equivalent of `sum(a*log(a/b)-a+b)`.
"""
function gkldiv(a::AbstractArray{T}, b::AbstractArray{T}) where T<:AbstractFloat
    n = length(a)
    r = 0.0
    for i = 1:n
        @inbounds ai = a[i]
        @inbounds bi = b[i]
        if ai > 0
            r += (ai * log(ai / bi) - ai + bi)
        else
            r += bi
        end
    end
    return r::Float64
end

# MeanAD: mean absolute deviation
"""
    meanad(a, b)

Return the mean absolute deviation between two arrays: `mean(abs(a - b))`.
"""
meanad(a::AbstractArray{T}, b::AbstractArray{T}) where {T<:Number} =
    L1dist(a, b) / length(a)

# MaxAD: maximum absolute deviation
"""
    maxad(a, b)

Return the maximum absolute deviation between two arrays: `maxabs(a - b)`.
"""
maxad(a::AbstractArray{T}, b::AbstractArray{T}) where {T<:Number} = Linfdist(a, b)

# MSD: mean squared deviation
"""
    msd(a, b)

Return the mean squared deviation between two arrays: `mean(abs2(a - b))`.
"""
msd(a::AbstractArray{T}, b::AbstractArray{T}) where {T<:Number} =
    sqL2dist(a, b) / length(a)

# RMSD: root mean squared deviation
"""
    rmsd(a, b; normalize=false)

Return the root mean squared deviation between two optionally
normalized arrays. The root mean squared deviation is computed
as `sqrt(msd(a, b))`.
"""
function rmsd(a::AbstractArray{T}, b::AbstractArray{T}; normalize::Bool=false) where T<:Number
    v = sqrt(msd(a, b))
    if normalize
        amin, amax = extrema(a)
        v /= (amax - amin)
    end
    return v
end

# PSNR: peak signal-to-noise ratio
"""
    psnr(a, b, maxv)

Compute the peak signal-to-noise ratio between two arrays `a` and `b`.
`maxv` is the maximum possible value either array can take. The PSNR
is computed as `10 * log10(maxv^2 / msd(a, b))`.
"""
function psnr(a::AbstractArray{T}, b::AbstractArray{T}, maxv::Real) where T<:Real
    20. * log10(maxv) - 10. * log10(msd(a, b))
end


text:
using RobustAndOptimalControl, ControlSystems, LinearAlgebra
# NOTE: Glad, Ljung chap 9 contains several numerical examples that can be used as test cases

approxin(el,col;kwargs...) = any(colel -> isapprox(el, colel; kwargs...), col)
approxsetequal(s1,s2;kwargs...) = all(approxin(p,s1;kwargs...) for p in s2) && all(approxin(p,s2;kwargs...) for p in s1)

w = exp10.(range(-5, stop=5, length=1000))
s = tf("s")
P = [1/(s+1) 2/(s+3); 1/(s+1) 1/(s+1)]
sys = ss(P)
sysmin = minreal(sys)
A,B,C,D = sys.A,sys.B,sys.C,sys.D
Am,Bm,Cm,Dm = sysmin.A,sysmin.B,sysmin.C,sysmin.D

@test approxsetequal(eigvals(Am), [-3,-1,-1])

Q1 = 100ones(4)
Q2 = 1I(2)
R1 = 100I(4)
R2 = 1I(2)
G = LQGProblem(sys, Q1, Q2, R1, R2)
gangoffourplot(G) # Test that it at least does not error
@test approxsetequal(eigvals(observer_controller(G).A), [ -31.6209+0.0im, -1.40629+0.0im, -15.9993+0.911174im, -15.9993-0.911174im, ], rtol = 1e-3)

# Test some utility functions
C = observer_controller(G)
@test C.A == sys.A - sys.B*lqr(G) - kalman(G)*sys.C
@test C.C == lqr(G)
@test C.B == kalman(G)
@test all(iszero, C.D)

Ce = extended_controller(G)
@test Ce.ny == C.ny
@test Ce.nu == C.nu
@test Ce.nw == size(Q1, 1)
@test system_mapping(Ce) == -C
Cp = performance_mapping(Ce)
@test Cp.ny == 0
@test Cp.A == C.A
@test all(iszero, Cp.B)
@test Cp.nu == Ce.nw

# testing with integral action
# qQ = 1
# qR = 1
# Q1 = 1000I(4)
# Q2 = 1I(2)
# R1 = 1I(6)
# R2 = 1I(2)
# N = I(6)
# Gi = LQGProblem(sys, Q1, Q2, R1, R2, qQ=qQ, qR=qR, integrator=true, ϵ=0.001, N=N)
# gangoffourplot(Gi) # Test that it at least does not error
# @test approxsetequal(eigvals(Gi.sysc.A), [-0.001, -0.001, -47.4832, -44.3442, -3.40255, -1.15355 ], rtol = 1e-2)

# @test approxsetequal(eigvals(closedloop(G).A), [-1.0, -14.1774, -2.21811, -14.3206, -1.60615, -22.526, -1.0, -14.1774], rtol=1e-3)
# @test approxsetequal(eigvals(G.T.A), [-22.526, -2.21811, -1.60615, -14.3206, -14.1774, -1.0, -1.0, -14.1774], rtol=1e-3)
# @test approxsetequal(eigvals(G.S.A), [-22.526, -2.21811, -1.60615, -14.3206, -14.1774, -1.0, -1.0, -14.1774], rtol=1e-3)
# @test approxsetequal(eigvals(G.CS.A), [-31.6209+0.0im, -1.40629+0.0im, -15.9993+0.911174im, -15.9993-0.911174im, -22.526+0.0im, -2.21811+0.0im, -1.60615+0.0im, -14.3206+0.0im, -14.1774+0.0im, -1.0+0.0im, -1.0+0.0im, -14.1774+0.0im], rtol=1e-3)
# @test approxsetequal(eigvals(G.PS.A), [-1.0, -1.0, -3.0, -1.0, -22.526, -2.21811, -1.60615, -14.3206, -14.1774, -1.0, -1.0, -14.1774], rtol=1e-3)

# @test approxsetequal(eigvals(Gi.cl.A), [-1.0, -44.7425, -44.8455, -2.23294, -4.28574, -2.06662, -0.109432, -1.31779, -0.78293, -1.0, -0.001, -0.001], rtol=1e-3)
# @test approxsetequal(eigvals(Gi.T.A), [-44.7425, -44.8455, -4.28574, -0.109432, -2.23294, -2.06662, -1.31779, -0.78293, -1.0, -1.0], rtol=1e-3)
# @test approxsetequal(eigvals(Gi.S.A), [-44.7425, -44.8455, -4.28574, -0.109432, -2.23294, -2.06662, -1.31779, -0.78293, -1.0, -1.0], rtol=1e-3)

# @test eigvals(Gi.sysc.A) == eigvals(Gi.Fr.A)

## Test cases from Glad Ljung chap 9

# Aircrapt control

A = [-0.292 8.13 -201 9.77 0 -12.5 17.1
    -0.152  -2.54  0.561  -0.0004  0  107  7.68
    0.0364  -0.0678  -0.481  0.0012  0  4.67  -7.98
    0 1 0.0401 0 0 0 0
    0 0 1 0 0 0 0
    0 0 0 0 0 -20 0
    0 0 0 0 0 0 -20
]

B = [
    0 -2.15
    -31.7 0.0274
    0 1.48
    0 0
    0 0
    20 0
    0 20
]
# Example 9.1
M = [0 0 0 1 0 0 0
     0 0 0 0 1 0 0]
C = M

N = I(7)
nx,ny,nu = size(A,1), size(C, 1), size(B,2)

Q1 = I(ny)
Q2 = I(nu)
R1 = I(nx)
R2 = I(ny)

sys = ss(A,N,B,M,C)
G = LQGProblem(sys, Q1, Q2, R1, R2)
@test lqr(G) ≈ [-0.0022 0.17 0.12 0.98 0.31 0.76 0.018
            -0.0038 0.028  -0.25 0.16  -0.95 0.072 0.10] rtol=0.01

Q2 = 0.1*I(nu)
G = LQGProblem(sys, Q1, Q2, R1, R2)
@test lqr(G) ≈ [-0.0036 0.39 0.21 3.11 0.63 1.54 0.046
            -0.0073 0.085 -0.78 0.57 -3.10 0.20 0.30] rtol=0.01

@test dcgain(RobustAndOptimalControl.closedloop(G)) ≈ I

# Room temperature control
K1 = 1
K2 = 1
K3 = 1/8
K4 = 0.2
A = [
    -K3-K4 K4 0 0 0
    K1 -K1-K2 0 K2 0
    0 0 -0.01 -0.068 0
    0 0 1 0 0
    0 0 0 0 -0.0001
]
# @warn "I added the last -0.0001"

B2 = [
    K3
    0
    0
    0
    0
]

C1 = [0 1 0 0 0]
C2 = [
    0 1 0 0 1
    0 1 0 0 0
    0 0 0 1 0
]
B1 = [
    0 0
    0 0
    1 0
    0 0
    0 1
]
sys = ss(A,B1,B2,C1,C2)

Q1 = 1.0I(1)
Q2 = 0.01I(1)
R1 = diagm([10.0, 1])
# R1 = diagm([0, 0, 10.0, 1, 0])
R2 = diagm([0.0001, 0.1, 0.001])
G = LQGProblem(sys, Q1, Q2, R1, R2, qQ = 0.000001)

# NOTE
# the kalman filter qutions will not work out until the noise model is incorporated

@test kalman(G) ≈ [-0.0039  0.0003  0.0035
            0.0463  0.0027  0.9676
            8.5224  0.0148  99.7543
            1.1426  0.0097  14.1198
            99.9634  -0.0026 -0.8533] rtol=1e-2
@test lqr(G) ≈ [2.8106 1.4071 6.8099 3.7874 0] rtol=1e-2

@test dcgain(RobustAndOptimalControl.closedloop(G))[] ≈ 1

## example 9.3

A = diagm([-1.0, -1, -2])
B = [
    1.0 0
    0.0 1
    0.0 1
]

C = [
    2.0 0 3
    1.0 1 0
]
M = C
sys = ss(A,I(3), B, M, C)

Q1 = I(2)
Q2 = 0.1I(2)
R1 = I(2)
R2 = 0.001I(2)

# test with integrator
# G = LQGProblem(sys, Q1, Q2, R1, R2, integrator=true, M = M, measurement=true)

# Lr = pinv(M * ((B * lqr(G)[:, 1:3] - A) \ B))

# @test Lr ≈ [
#     4.24 -3.43
#     0.12 4.82
# ] rtol = 0.01

# # NOTE: Glad ljung does not have identity in the end here, but the eigenvalues of A-BL are identical for this L and their L
# @test lqr(G) ≈ [
#     4.05 0.81 4.23 1 0
#     5.05 1.01 5.96 0 1
# ] rtol = 0.01

# @test kalman(G) ≈ [
#     0 0
#     0 0
#     0 0
#     31.62 0
#     0 31.62
# ] rtol = 0.01

# @test dcgain(closedloop(G)) ≈ I

# Random system, compare with ML output

A = [  0.154388  1.42524    0.914886  -0.313071
1.36189   0.30188   -0.767355   0.313058
-0.485766  0.926861  -0.755906  -0.309824
2.10783   0.146366  -1.65582    0.115643]

B = [
  0.625203  -1.87059
 -0.399363   0.516783
 -0.150984   0.161573
 -0.432394  -0.162178
]

C = [
  0.99517   -0.26967    0.255771  -0.0356366
 -0.358348   0.114812  -0.417828   0.983806
  0.45019   -0.166284   1.47825    2.21589
]

K = [ 1.22237  -0.114471  -1.10722
1.00941   0.445568  -1.91614]

P = ss(A,B,C,0)
C = ss(K)

SiA = [
 -0.3534     1.707    -3.451    -5.871
   1.674    0.1881    0.1495     1.257
 -0.3683    0.8871   -0.5026  -0.08282
   2.433   0.09358    -2.655    -1.636
]

SiB = [
  0.6252   -1.871
 -0.3994   0.5168
  -0.151   0.1616
 -0.4324  -0.1622
]

SiC = [
   -0.759    0.1587     1.276      2.61
  0.01776  -0.09757     2.761     3.844
]

SiD = [
    1      0
    0      1]

Si = ss(SiA, SiB, SiC, SiD)

@test isapprox(input_sensitivity(P, C), Si, rtol=1e-3)

SoA = [
    -0.3534     1.707    -3.451    -5.871
      1.674    0.1881    0.1495     1.257
    -0.3683    0.8871   -0.5026  -0.08282
      2.433   0.09358    -2.655    -1.636
]

SoB = [
    1.124     0.905    -2.892
 -0.03348    -0.276     0.548
  0.02146  -0.08928    0.1424
   0.6922   0.02276   -0.7895
]

SoC = [
   0.9952   -0.2697    0.2558  -0.03564
  -0.3583    0.1148   -0.4178    0.9838
   0.4502   -0.1663     1.478     2.216
]

SoD = I(3)

So = ss(SoA, SoB, SoC, SoD)

@test isapprox(output_sensitivity(P, C), So, rtol=1e-3)

ATo  = [
      -0.3534     1.707    -3.451    -5.871
        1.674    0.1881    0.1495     1.257
      -0.3683    0.8871   -0.5026  -0.08282
        2.433   0.09358    -2.655    -1.636
]

BTo = [
    1.124     0.905    -2.892
    -0.03348    -0.276     0.548
    0.02146  -0.08928    0.1424
    0.6922   0.02276   -0.7895
    ]

CTo = [
    -0.9952   0.2697  -0.2558  0.03564
    0.3583  -0.1148   0.4178  -0.9838
    -0.4502   0.1663   -1.478   -2.216
]
To = ss(ATo, BTo, CTo, 0)

@test isapprox(output_comp_sensitivity(P, C), To, rtol=1e-3)

ATi = [
      -0.3534     1.707    -3.451    -5.871
        1.674    0.1881    0.1495     1.257
      -0.3683    0.8871   -0.5026  -0.08282
        2.433   0.09358    -2.655    -1.636
]

BTi = [
    0.6252   -1.871
    -0.3994   0.5168
    -0.151   0.1616
    -0.4324  -0.1622
    ]

CTi = [
0.759   -0.1587    -1.276     -2.61
    -0.01776   0.09757    -2.761    -3.844
]

Ti = ss(ATi, BTi, CTi, 0)

@test isapprox(input_comp_sensitivity(P, C), Ti, rtol=1e-3)


text:
using DiffEqBase, DiffEqBase.InternalEuler

# Try it
ff = ODEFunction(
    (u, p, t) -> u,
    jac=(u, p, t) -> 1.,
    analytic=(u0, p, t) -> u0*exp(t)
)

dt = 0.01
prob = ODEProblem(ff, 1.0, (0.0,1.0) )
sol = solve(prob, InternalEuler.FwdEulerAlg(), tstops=0:dt:1)
sol2 = solve(prob, InternalEuler.BwdEulerAlg(), tstops=0:dt:1)

#using Plots
#plot(sol)
#plot!(sol2)
#plot!(sol2, plot_analytic=true)

#using DiffEqDevTools
#dts = 1./2.^(8:-1:4)
#sim = test_convergence(dts,p2,BwdEulerAlg())
#@show sim.𝒪est[:final]
#plot(sim)


text:
# https://mrandri19.github.io/2020/09/28/bayesian-linear-regression-with-conjugate-priors.html
# https://qiita.com/haru1843/items/3956dab2fd0d448cd02b
# https://discourse.julialang.org/t/solving-ax-b-for-large-matrix-dimesnions-efficiently-in-julia/51504
# https://stackoverflow.com/questions/57270276/identity-matrix-in-julia

# bayesian inferense (using conjugate prior analyticaly) of coefficients and variance in error term for linear regression model

using Random, Plots, StatsPlots, Distributions, LinearAlgebra, DataFrames, Printf, LaTeXStrings

"""
Conjugate Priors for linear regression model (β|σ^2 ~ MvNormal, σ^2 ~ InverseGamma)

    Inputs
    ------
    y:    data, independent variable (n vector)
    X:    data, dependent variables (n * k matrix)
    β0:   mean of conditional prior (multivariate normal β0|σ^2 ~ MvNormal) for coefficients
    A0:   precision matrix (inverse of Σ) of conditional prior for coefficients
    ν0:   shape parameter of prior for variance of error term
    λ0:   scale parameter of prior for variance of error term
    prob: credible interval probability (0 < prob < 1)

    Outputs
    -------
    results: posterior summary stats table (DataFrames)
    marginal_posterior_β: vector of marginal posterior for each coefficients (single β|D ~ TDist)
    posterior_σ2: posterior for variance of error term (σ^2|D ~ InverseGamma)
"""
function regression_stats(y, X, β0, A0, ν0, λ0, prob)
    XX = X' * X
    Xy = X' * y
    # @time XX \ Xy # fast and efficient
    # @time inv(XX) * Xy
    β_ols = XX \ Xy # OLS estimator
    A_star = XX + A0
    β_star = A_star \ (Xy + A0 * β0)
    C_star = inv(inv(XX) + inv(A0))
    ν_star = size(y)[1] + ν0
    ssr = sum((y - X * β_ols).^2) # SSR sum of squared residual 残差二乗和
    # @time (β0 - β_ols)' * C_star * (β0 - β_ols)
    # @time dot(β0 - β_ols, C_star, β0 - β_ols) # faster
    λ_star = ssr + dot(β0 - β_ols, C_star, β0 - β_ols) + λ0
    H_star = (λ_star / ν_star) * inv(A_star)
    h_star = diag(sqrt(H_star))

    marginal_posterior_β = LocationScale.(β_star, h_star, TDist(ν_star)) # broadcasting through coefficients
    sd_β = std.(marginal_posterior_β)
    ci_β_upper = quantile.(marginal_posterior_β, prob - (1-prob)/2)
    ci_β_lower = quantile.(marginal_posterior_β, (1-prob)/2)
    # ci_β = [ci_β_lower, ci_β_upper]
    stats_β = [β_star β_star β_star sd_β ci_β_lower ci_β_upper]

    posterior_σ2 = InverseGamma(ν_star/2, λ_star/2)
    mean_σ2 = mean(posterior_σ2)
    median_σ2 = median(posterior_σ2)
    mode_σ2 = mode(posterior_σ2)
    sd_σ2 = std(posterior_σ2)
    ci_σ2_upper = quantile(posterior_σ2, prob - (1-prob)/2)
    ci_σ2_lower = quantile(posterior_σ2, (1-prob)/2)
    # ci_σ2 = [ci_σ2_lower, ci_σ2_upper]
    stats_σ2 = [mean_σ2 median_σ2 mode_σ2 sd_σ2 ci_σ2_lower ci_σ2_upper]
    stats = [stats_β; stats_σ2]
    stats_string = ["Mean", "Median", "Mode", "SD", @sprintf("CI(%d)_Lower", prob*100), @sprintf("CI(%d)_Upper", prob*100)]
    results = DataFrame(stats, stats_string)
    return results, marginal_posterior_β, posterior_σ2, β_ols
end

# Generate moc data
n = 50
dim = 1
Random.seed!(99)
u = rand(Normal(0, 0.7), n) # true value σ^2 = 0.49
x = rand(Uniform(-sqrt(3.0), sqrt(3.0)), n, dim) # x generated from Uniform[-√3, √3]
X = [ones(n) x] # dependent variable matrix with ones
β = Vector(1:1:dim+1) # true coefficients
y = X*β + u # true data generating process

# set prior
k = size(X, 2)
β0 = zeros(k)
τ0 = 0.2 # no ridge penalty ... τ0 = Inf ?

# τ0 = 0.001
# τ0 = 0.01
# τ0 = 0.1
# τ0 = 1
# τ0 = 10
# τ0 = 100
# τ0 = 1000

A0 = τ0 * I # Identity matrix
ν0 = 5.0
λ0 = 7.0
H0 = (λ0 / ν0) * inv(A0) # scale matrix for marginal prior of multivariate coefficient (mv Tdist)
h0 = diag(sqrt(Matrix(H0, k, k))) # scale params for marginal prior of single coefficient (T dist) H0の対角成分
marginal_prior_β = LocationScale.(β0, h0, TDist(ν0)) # broadcasting through coefficients
prior_σ2 = InverseGamma(ν0/2, λ0/2)

# calculate posterior
prob = 0.95
results, marginal_posterior_β, posterior_σ2, β_ols = regression_stats(y, X, β0, A0, ν0, λ0, prob)
print("$results\n")

# Visualize Posterior
plots = []

for i = 1:k
    label_prior = @sprintf("marginal prior: T(ν=%.1f, μ=%.1f, σ^2=%.1f)", params(params(marginal_prior_β[i])[3])[1], params(marginal_prior_β[i])[1], params(marginal_prior_β[i])[2])
    p = plot(marginal_prior_β[i], linestyle=:dash, label=label_prior, xlims=(0.0, 3.0), ylims=(0.0, 4.0), xlabel=@sprintf("β%d", i), ylabel="Probability Density")
    label_posterior = @sprintf("marginal posterior: T(ν=%.1f, μ=%.1f, σ^2=%.1f)", params(params(marginal_posterior_β[i])[3])[1], params(marginal_posterior_β[i])[1], params(marginal_posterior_β[i])[2])
    plot!(marginal_posterior_β[i], linestyle=:solid, label=label_posterior)
    push!(plots, p)
end

label_prior = @sprintf("prior: InverseGamma(α=%.1f, θ=%.1f)", params(prior_σ2)[1], params(prior_σ2)[2])
p = plot(prior_σ2, linestyle=:dash, label=label_prior, xlims=(0.0, 3.0), ylims=(0.0, 4.0), xlabel=L"\sigma^2", ylabel="Probability Density")
label_posterior = @sprintf("prior: InverseGamma(α=%.1f, θ=%.1f)", params(posterior_σ2)[1], params(posterior_σ2)[2])
plot!(posterior_σ2, linestyle=:solid, label=label_posterior)
push!(plots, p)

plot(plots..., layout=(1, k+1))

# plot ols
scatter(x, y, label="data")
# scatter(x, y, xlims=(-0.1, 0.1), ylims=(0.9, 1.5))
plot!(x, X * β, label=(@sprintf("True: y = %.d + %.dx", β[1], β[2])))
plot!(x, X * β_ols, label=(@sprintf("OLS: y = %.5f + %.5fx", β_ols[1], β_ols[2])))
β_star = results[:, "Mean"][1:k]
plot!(x, X * β_star, label=(@sprintf("Bayes: y = %.5f + %.5fx (τ0=%.2f, Ridge)", β_star[1], β_star[2], τ0)))

# significance test on coefficients H0: βj = 0, H1: βj ≠ 0
# evaluate with SDDR
"""
calculate SDDR (compare density of posterior and prior on particular value)
# https://docs.julialang.org/en/v1/manual/methods/
# https://stackoverflow.com/questions/61088294/declaring-the-name-of-argument-when-invoking-a-function

    Inputs
    ------
    β_h0: test value, (default is 0)
    posterior
    prior

    Output
    ------
    sddr
"""
function sddr(posterior, prior, ;null_hypothesis::Float64=0.0)
    return pdf(posterior, null_hypothesis) / pdf(prior, null_hypothesis)
end

print(log10.(sddr.(marginal_posterior_β, marginal_prior_β, null_hypothesis=0.0)))


text:
@testset "FisherTest()" begin
    @test FisherTest() isa FisherTest
end

@testset "Escape.run(::FisherTest, ::AbstractHLAData, ::Replacement)" begin
    hla_types = [HLAType(parse_allele("A01", "A02", "B05", "B06", "C03", "C07")),
                 HLAType(parse_allele("A01", "A03", "B03", "B05", "C03", "C07")),
                 HLAType(parse_allele("A01", "A04", "B03", "B05", "C03", "C07")),
                 HLAType(parse_allele("A02", "A04", "B03", "B04", "C03", "C07")),
                 HLAType(parse_allele("A02", "A04", "B03", "B04", "C03", "C07"))]

    fasta_path = joinpath(@__DIR__, "data", "test.fasta")
    replacement = Replacement("test", 2, 'S', false)
    hla_data = HLAData("test", fasta_path, hla_types, missing, missing)

    y = Bool[0, 0, 1, 1, 1] # representation for replacement
    x = Bool[1, 1, 1, 0, 0] # allele vector for A01

    result = Escape.run(FisherTest(), hla_data, replacement)
    counts, p, beta = Escape.fisher_exact_test(y, x)

    @test result.counts[1][parse_allele("A01")] == counts
    @test result.p_values[1][parse_allele("A01")] == p
    @test result.log_odds[1][parse_allele("A01")] == beta
end

@testset "fisher_exact_test(::Vector{Bool}, ::Vector{Bool}" begin
    a = Bool[1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1] # y's
    b = Bool[0, 0, 1, 1, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0] # alleles

    #             with allele  without allele
    # with y    |      1              7
    # without y |      4              2
    expected_result = [1 7; 4 2]
    counts, p, beta = Escape.fisher_exact_test(a, b)

    @test counts == [1 7; 4 2]
    @test p == pvalue(FisherExactTest(1, 7, 4, 2), method = :minlike)
    @test beta == log(FisherExactTest(1, 7, 4, 2).ω)
end


text:
# included from Potentials.jl
# part of the module JuLIP.Potentials

using JuLIP: JVec, JMat, neighbourlist
using LinearAlgebra: I
using JuLIP.Chemistry: atomic_number
using NeighbourLists

export ZeroPairPotential, ZBLPotential,
         LennardJones, lennardjones,
         Morse, morse

# For PairPotentials we do something kind of weird - we create an intentional
# stack overflow by redirecting
#    evaluate -> evaluate! -> evaluate -> evaluate! -> ...
# A concrete instance of PairPotential must overload either evaluate or
# evaluate! to break this infinite loop. The advantage though is that either
# can be implemented.

evaluate(V::PairPotential, r::Number, z1, z0) =
      evaluate!(alloc_temp(V, 1), V, r, z1, z0)
evaluate_d(V::PairPotential, r::Number, z1, z0) =
      evaluate_d!(alloc_temp_d(V, 1), V, r, z1, z0)
evaluate_dd(V::PairPotential, r::Number, z1, z0) =
      evaluate_dd!(alloc_temp_dd(V, 1), V, r, z1, z0)
evaluate!(tmp, V::PairPotential, r::Number, z1, z0) =
      evaluate(V, r, z1, z0)
evaluate_d!(tmp, V::PairPotential, r::Number, z1, z0) =
      evaluate_d(V, r, z1, z0)
evaluate_dd!(tmp, V::PairPotential, r::Number, z1, z0) =
      evaluate_dd(V, r, z1, z0)

# evaluate!(tmp, V::PairPotential, r::Union{Number, JVec}) = V(r)
# evaluate_d!(tmp, V::PairPotential, r::Number) = @D V(r)
# evaluate_dd!(tmp, V::PairPotential, r::Number) = @DD V(r)
# evaluate_d!(tmp, V::PairPotential, R::JVec) =
#       evaluate_d!(tmp, V, norm(R), R)
# evaluate_d!(tmp, V::PairPotential, r::Number, R::JVec) = ((@D V(r))/r) * R
# evaluate_dd!(tmp, V::PairPotential, r::Number) = @DD V(r)
# evaluate_dd!(tmp, V::PairPotential, R::JVec) =
#       evaluate_dd!(tmp, V, norm(R), R)
# evaluate_dd!(tmp, V::PairPotential, r::Number, R::JVec) =
#       _hess!(tmp, V, r, R)

function evaluate!(tmp, V::PairPotential,
                   R::AbstractVector{JVec{T}}, Z, z0) where {T}
   Es = zero(T)
   for i = 1:length(R)
      Es += T(0.5) * evaluate!(tmp, V, norm(R[i]), Z[i], z0)
   end
   return Es
end

function evaluate_d!(dEs, tmp, V::PairPotential,
                     R::AbstractVector{JVec{T}}, Z, z0) where {T}
   for i = 1:length(R)
      r = norm(R[i])
      dEs[i] = (T(0.5) * evaluate_d!(tmp, V, r, Z[i], z0) / r) * R[i]
   end
   return dEs
end

function evaluate_dd!(hEs, tmp, V::PairPotential,
                      R::AbstractVector{<: JVec}, Z, z0)
   n = length(R)
   for i = 1:n
      hEs[i,i] = 0.5 * _hess!(tmp, V, norm(R[i]), R[i], Z[i], z0)
   end
   return hEs
end

function _hess!(tmp, V::PairPotential, r::Number, R::JVec, z1, z0)
   R̂ = R/r
   P = R̂ * R̂'
   dV = evaluate_d!(tmp, V, r, z1, z0) / r
   ddV = evaluate_dd!(tmp, V, r, z1, z0)
   return (ddV - dV) * P + dV * I
end

function precon!(hEs, tmp, V::PairPotential,
                 R::AbstractVector{<: JVec{T}}, Z, z0,
                 innerstab=T(0.0)) where {T}
   n = length(R)
   for i = 1:n
      hEs[i,i] = precon!(tmp, V, norm(R[i]), R[i], Z[i], z0, innerstab)
   end
   return hEs
end

# an FF preconditioner for pair potentials
function precon!(tmp, V::PairPotential, r::T, R::JVec{T}, innerstab=T(0.1)
                 ) where {T <: Number}
   r = norm(R)
   dV = evaluate_d!(tmp, V, r)
   ddV = evaluate_dd!(tmp, V, r)
   R̂ = R/r
   return (1-innerstab) * (abs(ddV) * R̂ * R̂' + abs(dV / r) * (I - R̂ * R̂')) +
             innerstab  * (abs(ddV) + abs(dV / r)) * I
end

# ------- Implementation of SimplePairPotential

evaluate!(tmp, V::SimplePairPotential, r::Number, z1, z0) = evaluate!(tmp, V, r)
evaluate_d!(tmp, V::SimplePairPotential, r::Number, z1, z0) = evaluate_d!(tmp, V, r)
evaluate_dd!(tmp, V::SimplePairPotential, r::Number, z1, z0) = evaluate_dd!(tmp, V, r)
evaluate!(tmp, V::SimplePairPotential, r) = evaluate(V, r)
evaluate_d!(tmp, V::SimplePairPotential, r) = evaluate_d(V, r)
evaluate_dd!(tmp, V::SimplePairPotential, r) = evaluate_dd(V, r)
evaluate(V::SimplePairPotential, r::Number) = evaluate!(alloc_temp(V, 1), V, r)
evaluate_d(V::SimplePairPotential, r::Number) = evaluate_d!(alloc_temp_d(V, 1), V, r)
evaluate_dd(V::SimplePairPotential, r::Number) = evaluate_dd!(alloc_temp_dd(V, 1), V, r)

evaluate_d(V::SimplePairPotential, r::Number, R::JVec) =
      (evaluate_d(V, r)/r) * R

function evaluate_dd(V::SimplePairPotential, r::Number, R::JVec)
   dV = evaluate_d(V, r) / r
   ddV = evaluate_dd(V, r)
   return ((ddV - dV)/r^2) * R * R'  + dV * I
end

# ------- Implementation of ExplicitPairPotential

evaluate(p::ExplicitPairPotential, r::Number) = p.f(r)
evaluate_d(p::ExplicitPairPotential, r::Number) = p.f_d(r)
evaluate_dd(p::ExplicitPairPotential, r::Number) = p.f_dd(r)

# ------- Some concrete potentials

"""
`LennardJones(σ, e0):` constructs the 6-12 Lennard-Jones potential [wiki](https://en.wikipedia.org/wiki/Lennard-Jones_potential)

   e0 * 4 * ( (σ/r)¹² - (σ/r)⁶ )

Constructor with kw arguments: `LennardJones(; kwargs...)`

* `e0, σ` : standard LJ parameters, default `e0 = 1.0, σ = 1.0`
* `r0` : equilibrium distance - if `r0` is specified then `σ` is ignored
* `a0` : FCC lattice parameter, if `a0` is specified then `σ` is ignored

(`r0, a0` cannot both be specified at the same time)
"""
LennardJones(σ, e0) = (@analytic r -> e0 * 4.0 * ((σ/r)^(12) - (σ/r)^(6)))

function ljparams(; σ=1.0, e0=1.0, r0 = nothing, a0 = nothing)
   if r0 != nothing && a0 != nothing
      error("`LenndardJones`: cannot specify both `r0` and `a0`")
   end
   if a0 != nothing    # r0 = nn-dist = a0/sqrt(2) in FCC
      r0 = a0 / sqrt(2)
   end
   if r0 != nothing    # standard LJ is minimised at r = 2^(1/6)
      σ =  r0 / 2^(1/6)
   end
   return σ, e0
end

LennardJones(; kwargs...) = LennardJones(ljparams(;kwargs...)...)

"""
`lennardjones(; kwargs...)`

simplified constructor for `LennardJones` (note this is type unstable!)

In addition to the `kwargs` of `LennardJones`, this accepts also

* `rcut` : default `:auto` which gives `rcut = (1.9*σ, 2.7*σ)`. Use
`nothing` or `Inf` to specify no cutoff, or specify a tuple or
array with two elements specifying the lower and upper cut-off radii to be
used with `SplineCutoff`.
"""
function lennardjones(; rcut = :auto, kwargs...)
   σ, e0 = ljparams(; kwargs...)
   if (rcut == nothing || rcut == Inf)
      return LennardJones(σ, e0)
   elseif rcut == :auto
      rcut = (1.9*σ, 2.7*σ)
   end
   return SplineCutoff(rcut[1], rcut[2]) * LennardJones(σ, e0)
end

"""
`Morse(A, e0, r0)` or `Morse(;A=4.0, e0=1.0, r0=1.0)`: constructs a
`PairPotential` for
```
   e0 ( exp( -2 A (r/r0 - 1) ) - 2 exp( - A (r/r0 - 1) ) )
```
"""
Morse(A, e0, r0) = @analytic(
   r -> e0 * ( exp(-(2.0*A) * (r/r0 - 1.0)) - 2.0 * exp(-A * (r/r0 - 1.0)) ) )
Morse(;A=4.0, e0=1.0, r0=1.0) = Morse(A, e0, r0)

"""
`morse(A=4.0, e0=1.0, r0=1.0, rcut=(1.9*r0, 2.7*r0))`

simplified constructor for `Morse` (type unstable)
"""
morse(;A=4.0, e0=1.0, r0=1.0, rcut=(1.9*r0, 2.7*r0)) = (
   (rcut == nothing || rcut == Inf)
         ?  Morse(A, e0, r0)
         :  SplineCutoff(rcut[1], rcut[2]) * Morse(A, e0, r0) )

"""
`ZeroPairPotential()`: creates a potential that just returns zero
"""
struct ZeroPairPotential <: SimplePairPotential end

@pot ZeroPairPotential

evaluate(p::ZeroPairPotential, r::T) where {T <: Number} = T(0.0)
evaluate_d(p::ZeroPairPotential, r::T) where {T <: Number} = T(0.0)
evaluate_dd(p::ZeroPairPotential, r::T) where {T <: Number} = T(0.0)
cutoff(p::ZeroPairPotential) = Bool(0) # the weakest number type

# ====================================================================
#   A product of two pair potentials: primarily used for cutoff mechanisms

"product of two `SimplePairPotential`"
mutable struct SimpleProdPot{P1, P2} <: SimplePairPotential
   p1::P1
   p2::P2
end

@pot SimpleProdPot

import Base.*
*(p1::SimplePairPotential, p2::SimplePairPotential) = SimpleProdPot(p1, p2)
evaluate(p::SimpleProdPot, r::Number) = p.p1(r) * p.p2(r)
evaluate_d(p::SimpleProdPot, r::Number) = (p.p1(r) * (@D p.p2(r)) + (@D p.p1(r)) * p.p2(r))
evaluate_dd(p::SimpleProdPot, r::Number) = (p.p1(r) * (@DD p.p2(r)) +
              2 * (@D p.p1(r)) * (@D p.p2(r)) + (@DD p.p1(r)) * p.p2(r))
cutoff(p::SimpleProdPot) = min(cutoff(p.p1), cutoff(p.p2))

"product of two `PairPotential`"
mutable struct ProdPot{P1, P2} <: PairPotential
   p1::P1
   p2::P2
end

@pot ProdPot

import Base.*
*(p1::PairPotential, p2::PairPotential) = ProdPot(p1, p2)
evaluate(p::ProdPot, r::Number, z1, z0) =
      p.p1(r, z1, z0) * p.p2(r, z1, z0)
evaluate_d(p::ProdPot, r::Number, z1, z0) =
      (p.p1(r, z1, z0) * (@D p.p2(r, z1, z0)) + (@D p.p1(r, z1, z0)) * p.p2(r, z1, z0))
evaluate_dd(p::ProdPot, r::Number, z1, z0) =
      (   p.p1(r, z1, z0) * (@DD p.p2(r, z1, z0))
        + 2 * (@D p.p1(r, z1, z0)) * (@D p.p2(r, z1, z0))
        + (@DD p.p1(r, z1, z0)) * p.p2(r, z1, z0)
      )
cutoff(p::ProdPot) = min(cutoff(p.p1), cutoff(p.p2))

# ====================================================================

"""
`struct WrappedPairPotential`

wraps a pairpotential using `FunctionWrappers` in order to allow
type-stable storage of multiple potentials. This is the main technique
required at the moment to work with multi-component systems.
Otherwise, this is not advisable since it disables a range of
possible compiler optimisations.
"""
struct WrappedPairPotential <: ExplicitPairPotential
   f::F64fun
   f_d::F64fun
   f_dd::F64fun
   rcut::Float64
end

@pot WrappedPairPotential

cutoff(V::WrappedPairPotential) = V.rcut
# evaluate, etc are all derived from SimplePairPotential

function WrappedPairPotential(V::AnalyticFunction, rcut)
   @assert (0 < rcut < Inf)
   f, f_d, f_dd = let V=V, rc = rcut
      (F64fun(r -> evaluate(V, r) * (r<rc)),
              F64fun(r -> evaluate_d(V, r) * (r<rc)),
              F64fun(r -> evaluate_dd(V, r) * (r<rc)))
   end
   return WrappedPairPotential(f, f_d, f_dd, cutoff(V))
end

function WrappedPairPotential(V::PairPotential)
   @assert (0 < cutoff(V) < Inf)
   f, f_d, f_dd = let V=V, rc = cutoff(V)
      (F64fun(r -> evaluate(V, r) * (r<rc)),
              F64fun(r -> evaluate_d(V, r) * (r<rc)),
              F64fun(r -> evaluate_dd(V, r) * (r<rc)))
   end
   return WrappedPairPotential(f, f_d, f_dd, cutoff(V))
end

# ------------------------------------------------------------------------

# TODO: write more docs + tests for ZBL

"""
Implementation of the ZBL potential to model close approach.
"""
struct ZBLPotential{TV} <: PairPotential
   V::TV   # analytic "inner" potential
end

@pot ZBLPotential

ZBLPotential() = (
   let
      # au = 0.8854 * 0.529 / (Z1^0.23 + Z2^0.23)
      ϵ0 = 0.00552634940621
      C = 1/(4*π*ϵ0)
      E1, E2, E3, E4 = 0.1818, 0.5099, 0.2802, 0.02817
      A1, A2, A3, A4 = 3.2, 0.9423, 0.4028, 0.2016
      V = @analytic(r -> C * (E1*exp(-A1*r) + E2*exp(-A2*r) +
                              E3*exp(-A4*r) + E4*exp(-A4*r) ) / r)
      ZBLPotential(V)
   end)

_zbl_au(Z1, Z2) = (0.8854 * 0.529) / (Z1^0.23 + Z2^0.23)

function evaluate(V::ZBLPotential, r::Number, z1, z0)
   au = _zbl_au(z1, z0)
   return evaluate(V.V, r / au)
end

function evaluate_d(V::ZBLPotential, r::Number, z1, z0)
   au = _zbl_au(z1, z0)
   return evaluate_d(V.V, r / au) / au
end

cutoff(::ZBLPotential) = Inf

write_dict(V::ZBLPotential) = Dict("__id__" => "JuLIP_ZBLPotential")
read_dict(::Val{:JuLIP_ZBLPotential}, D::Dict) = ZBLPotential()


text:
# Tests for Dirichlet distribution

using  Distributions
using Test, Random, LinearAlgebra

Random.seed!(34567)

rng = MersenneTwister(123)

@testset "Testing Dirichlet with $key" for (key, func) in
    Dict("rand(...)" => [rand, rand],
         "rand(rng, ...)" => [dist -> rand(rng, dist), (dist, n) -> rand(rng, dist, n)])

    for T in (Int, Float64)
        d = Dirichlet(3, T(2))

        @test length(d) == 3
        @test eltype(d) === T
        @test d.alpha == [2, 2, 2]
        @test d.alpha0 == 6

        @test mean(d) ≈ fill(1/3, 3)
        @test mode(d) ≈ fill(1/3, 3)
        @test cov(d)  ≈ [8 -4 -4; -4 8 -4; -4 -4 8] / (36 * 7)
        @test var(d)  ≈ diag(cov(d))

        r = Vector{Float64}(undef, 3)
        Distributions.dirichlet_mode!(r, d.alpha, d.alpha0)
        @test r ≈ fill(1/3, 3)

        @test pdf(Dirichlet([1, 1]), [0, 1]) ≈ 1
        @test pdf(Dirichlet([1f0, 1f0]), [0f0, 1f0]) ≈ 1
        @test typeof(pdf(Dirichlet([1f0, 1f0]), [0f0, 1f0])) === Float32

        @test iszero(pdf(d, [-1, 1, 0]))
        @test iszero(pdf(d, [0, 0, 1]))
        @test pdf(d, [0.2, 0.3, 0.5]) ≈ 3.6
        @test pdf(d, [0.4, 0.5, 0.1]) ≈ 2.4
        @test logpdf(d, [0.2, 0.3, 0.5]) ≈ log(3.6)
        @test logpdf(d, [0.4, 0.5, 0.1]) ≈ log(2.4)

        x = func[2](d, 100)
        p = pdf(d, x)
        lp = logpdf(d, x)
        for i in 1 : size(x, 2)
            @test lp[i] ≈ logpdf(d, x[:,i])
            @test p[i]  ≈ pdf(d, x[:,i])
        end

        v = [2, 1, 3]
        d = Dirichlet(T.(v))

        @test eltype(d) === T
        @test Dirichlet([2, 1, 3]).alpha == d.alpha

        @test length(d) == length(v)
        @test d.alpha == v
        @test d.alpha0 == sum(v)
        @test d == Dirichlet{T}(params(d)...)
        @test d == deepcopy(d)

        @test mean(d) ≈ v / sum(v)
        @test cov(d)  ≈ [8 -2 -6; -2 5 -3; -6 -3 9] / (36 * 7)
        @test var(d)  ≈ diag(cov(d))

        @test pdf(d, [0.2, 0.3, 0.5]) ≈ 3
        @test pdf(d, [0.4, 0.5, 0.1]) ≈ 0.24
        @test logpdf(d, [0.2, 0.3, 0.5]) ≈ log(3)
        @test logpdf(d, [0.4, 0.5, 0.1]) ≈ log(0.24)

        x = func[2](d, 100)
        p = pdf(d, x)
        lp = logpdf(d, x)
        for i in 1 : size(x, 2)
            @test p[i]  ≈ pdf(d, x[:,i])
            @test lp[i] ≈ logpdf(d, x[:,i])
        end

        # Sampling

        x = func[1](d)
        @test isa(x, Vector{Float64})
        @test length(x) == 3

        x = func[2](d, 10)
        @test isa(x, Matrix{Float64})
        @test size(x) == (3, 10)

        v = [2, 1, 3]
        d = Dirichlet(Float32.(v))
        @test eltype(d) === Float32

        x = func[1](d)
        @test isa(x, Vector{Float32})
        @test length(x) == 3

        x = func[2](d, 10)
        @test isa(x, Matrix{Float32})
        @test size(x) == (3, 10)

        # Test MLE

        v = [2, 1, 3]
        d = Dirichlet(v)

        n = 10000
        x = func[2](d, n)
        x = x ./ sum(x, dims=1)

        r = fit_mle(Dirichlet, x)
        @test r.alpha ≈ d.alpha atol=0.25
        r = fit(Dirichlet{Float32}, x)
        @test r.alpha ≈ d.alpha atol=0.25

        # r = fit_mle(Dirichlet, x, fill(2.0, n))
        # @test r.alpha ≈ d.alpha atol=0.25
    end
end

@testset "Dirichlet: entropy" begin
    α = exp.(rand(2))
    @test entropy(Dirichlet(α)) ≈ entropy(Beta(α...))

    N = 10
    @test entropy(Dirichlet(N, 1)) ≈ -loggamma(N)
    @test entropy(Dirichlet(ones(N))) ≈ -loggamma(N)
end


text:
module ExtrapTypeStability

using Base.Test, Interpolations, DualNumbers

# Test type-stability of 1-dimensional extrapolation
f(x) = sin((x-3)*2pi/9 - 1)
xmax = 10
A = Float64[f(x) for x in 1:xmax]
itpg = interpolate(A, BSpline(Linear()), OnGrid())

schemes = (
    Flat,
    Linear,
    Reflect,
    Periodic
)

for etp in map(E -> extrapolate(itpg, E()), schemes),
    x in [
        # In-bounds evaluation
        3.4, 3, dual(3.1),
        # Out-of-bounds evaluation
        -3.4, -3, dual(-3,1),
        13.4, 13, dual(13,1)
    ]
    @inferred(getindex(etp, x))
end

# Test type-stability of 2-dimensional extrapolation with homogeneous scheme
g(y) = (y/100)^3
ymax = 4
A = Float64[f(x)*g(y) for x in 1:xmax, y in 1:ymax]
itp2 = interpolate(A, BSpline(Linear()), OnGrid())

for (etp2,E) in map(E -> (extrapolate(itp2, E()), E), schemes),
    x in (
        # In-bounds evaluation
        3.4, 3, dual(3.1),
        # Out-of-bounds evaluation
        -3.4, -3, dual(-3,1),
        13.4, 13, dual(13,1)
    ),
    y in (
        # In-bounds evaluation
        2.1, 2, dual(2.3, 1),
        # Out-of-bounds evaluation
        -2.1, -2, dual(-2.3, 1),
        12.1, 12, dual(12.1, 1)
    )
    @inferred(getindex(etp2, x, y))
end

end


text:
@testset "132.palindrome-partitioning-ii.jl" begin
    @test partition_ps132("aab") == 1
    @test partition_ps132("a") == 0
    @test partition_ps132("ab") == 1
    @test partition_ps132(
        "iquvgtinxldylufwujdosgtqogqdibxeygwnbqbcuyxydkytolhvowhishkvukqoybxmdjajqtdvttmicmlmwerqhpheeckxckzuixsmoniyasaqxdwirntgebirlvzttabegwdseluyiktjglbsxivwkqowwtqbrciowkmdvhfsoslcspnukxxnsgbvpxddqqjkdzki",
    ) == 178
end


text:
const SingleAny = Tuple{T1,T2,SingleComponent} where {T1,T2}
const MultiAny = Tuple{T1,T2,MaterialCompounds} where {T1,T2}

struct IdealGas{T,MW,CP} <: GasModel
    type::T
    mw::MW
    cp::CP
end

IdealGas() = IdealGas(SingleComponent(),nothing,nothing)
function IdealGas(mw)
    if length(mw) == 1
        return IdealGas(SingleComponent(),only(mw),nothing)
    else
        return IdealGas(MaterialCompounds{MOLAR,FRACTION}(),mw,nothing)
    end
end

function IdealGas(mw,cp)
    if length(mw) == 1
        return IdealGas(SingleComponent(),only(mw),cp)
    else
        return IdealGas(MaterialCompounds{MOLAR,FRACTION}(),mw,cp)
    end
end

function pressure_impl(mt::SingleVT,model::IdealGas,v,t)
    return RGAS*t/v
end

function pressure_impl(mt::MultiVT,model::IdealGas,v,t,x)
    return RGAS*t/v
end

function temperature_impl(mt::SinglePV,model::IdealGas,p,v)
    return RGAS*t/v
end

function temperature_impl(mt::MultiPV,model::IdealGas,p,v,x)
    return p*v/RGAS
end

function compressibility_factor(mt::SingleAny,model::IdealGas,a1,a2)
    return one(a1)*one(a2)
end

function compressibility_factor(mt::MultiAny,model::IdealGas,a1,a2,a3)
    return one(a1)*one(a2)*one(eltype(a3))
end

function mol_volume_impl(mt::SinglePT,model::IdealGas,p,t)
    return RGAS*T/p
end

function mol_volume_impl(mt::MultiPT,model::IdealGas,p,t,x)
    return RGAS*T/p
end


text:
using LinearAlgebra
using CVChannel
using Convex
using Test
using DelimitedFiles
"""
In this script we investigate a family of co-positive maps.
Specifically, we investigate the set of states presented in
Eqn. (82) of https://arxiv.org/abs/1004.1655

Note we know co-positive maps have multiplicativity, so the only way to check this
at all is implementing the DPS hierarchy. Here we look at the 3-symmetric
extension of the communication value, and that there exist cases where it is
better.
"""

function boundBell(ε :: Union{Int,Float64}) :: Matrix
    if ε <= 0
        throw(DomainError(ε, "Make sure ε > 0."))
    end
    #This function generates the states defined in (82)
    id_mat = Matrix(1I, 3, 3)
    bell = 1/3*vec(id_mat)*vec(id_mat)';
    #Π operators defined in (40) and implicitly using (35)
    Π1 = zeros(9,9); Π2 = zeros(9,9)
    for i = 0:2
        Π1 = Π1 + kron(id_mat,discreteWeylOperator(i,1,3))*bell*kron(id_mat,discreteWeylOperator(i,1,3)')
        Π2 = Π2 + kron(id_mat,discreteWeylOperator(i,2,3))*bell*kron(id_mat,discreteWeylOperator(i,2,3)')
    end
    #They don't say to normalize it, but otherwise it isn't a state
    Π1 = 1/tr(Π1) * Π1 ; Π2 = 1/tr(Π2) * Π2
    N_ε = 1/(1+ε+ε^(-1))

    return N_ε*(bell + ε*Π1 + ε^(-1)*Π2)
end

println("As a sanity check, we first see they are PPT states.")
println("Because of numerical error, we can't just call the isPPT function")
@testset "Verifying class of states are PPT" begin
    ε_checks = [0.1,0.25,0.5,0.75,1,1.5,2,3,5]
    fail_state = false
    for ε_id in ε_checks
        test_state = boundBell(ε_id)
        realvals = real.(eigvals(partialtranspose(test_state,2,[3,3])))
        imvals = imag.(eigvals(partialtranspose(test_state,2,[3,3])))
        if !isapprox(tr(test_state),1,atol = 1e-9)
            fail_state = true
            break
        elseif !all(imvals -> imvals < 1e-10, imvals)
            fail_state = true
            break
        elseif !all(realvals -> realvals > 0 || isapprox(realvals,0,atol=1e-8),realvals)
            fail_state = true
            break
        else
        end
    end
    @test !fail_state
end

println("\nAs a second sanity check, we make sure they are multiplicative with eachother.")
println("This also gives us a chance to initialize the solver.")
@testset "Verifying PPT multiplicativity" begin
    test_state1 = 3*boundBell(0.25)
    test_state2 = 3*boundBell(1.5)
    println("Initializing solver, please wait...")
    cv1 = pptCVDual(test_state1,3,3)
    println("All initialized. The rest should be rather quick.")
    cv2 = pptCVDual(test_state2,3,3)
    state11 = permuteSubsystems(kron(test_state1,test_state1),[1,3,2,4],[3,3,3,3])
    state12 = permuteSubsystems(kron(test_state1,test_state2),[1,3,2,4],[3,3,3,3])
    state22 = permuteSubsystems(kron(test_state2,test_state2),[1,3,2,4],[3,3,3,3])
    cv11 = pptCVPrimal(state11,9,9)
    cv12 = pptCVPrimal(state12,9,9)
    cv22 = pptCVPrimal(state22,9,9)
    @test isapprox(cv11[1] - cv1[1]^2, 0, atol = 5e-6)
    @test isapprox(cv12[1] - cv1[1]*cv2[1], 0, atol = 5e-6)
    @test isapprox(cv22[1] - cv2[1]^2, 0, atol = 5e-6)
end

println("\nNow we can see that there exist cases where three symmetric communication value is tighter than two symmetric.")
println("WARNING: This takes a while.")
@testset "3-Sym Extension Improvment" begin
    scan_range = [0.45:0.02:0.75;]
    len = length(scan_range)
    data_table = zeros(len,6)
    ctr = 1
    for val in scan_range
        println("Now evaluating ", ctr," of ",length(scan_range)," points.")
        test_state1 = 3*boundBell(val)
        cv1 = pptCVDual(test_state1,3,3)
        cv2 = twoSymCVPrimal(test_state1,3,3)
        cv3 = threeSymCVPrimal(test_state1,3,3)
        data_table[ctr,:] = [val cv1[1] cv2[1] cv3[1] cv1[1]-cv2[1] cv2[1]-cv3[1]]
        ctr += 1
    end

    println("Here we again see there is indeed an improvement.")
    label_vec = ["input val" "PPT CV" "2 Sym CV" "3 Sym CV" "PPT-2Sym" "2Sym - 3Sym"]
    show(stdout, "text/plain", vcat(label_vec,data_table))

    info_vec = Vector{Union{Nothing,String}}(nothing, len+1)
    info_vec[1] = "INFO:"
    info_vec[2] = "Generated by co-positive-map-investigation.jl"
    info_vec[3] = "scan_range = " * string(scan_range)
    data_to_save = hcat(vcat(label_vec,data_table),info_vec)
    println("Please name the file you'd like to write the results to: \n")
    file_name = readline()
    file_to_open = string(file_name,".csv")
    writedlm(file_to_open, data_to_save, ',')

    if all(data_table -> data_table < 1e-4 , data_table[:,6])
        println("Something went awry. Please check your results.")
        three_sym_helps = false
    else
        three_sym_helps = true
    end
    @test three_sym_helps
end

max_forks_repo_name: ChitambarLab/cv-channel


text:
module MapQuest
  include("geocommon.jl")
  using HTTPClient.HTTPC
  using JSON

  base_url = "http://www.mapquestapi.com/geocoding/v1/address"

  type MapQuestClient
    key::String
  end

  function geocode(mq::MapQuestClient, location::String)
    resp = get(base_url,
               HTTPC.RequestOptions(;query_params=
                 [("key",mq.key),("location",location)]))
    json = JSON.parse(resp.body)
    return json
  end

  function simple_geocode(mq::MapQuestClient, location::String)
    result = geocode(mq,location)
    country = result["results"][1]["locations"][1]["adminArea1"]
    state   = result["results"][1]["locations"][1]["adminArea3"]
    county  = result["results"][1]["locations"][1]["adminArea4"]
    city    = result["results"][1]["locations"][1]["adminArea5"]
    lat     = result["results"][1]["locations"][1]["latLng"]["lat"]
    lng     = result["results"][1]["locations"][1]["latLng"]["lng"]
    GeoCommon.SimpleLocation("$city, $county, $state, $country",lat,lng)
  end

end

max_stars_repo_name: astrieanna/geojl

max_issues_repo_name: astrieanna/geojl

max_forks_repo_name: astrieanna/geojl


text:
# %%
import Pkg;
Pkg.activate(".")
# %%
path = pwd() * "/functions"
push!(LOAD_PATH,path)
# %%
using Revise
# %%
using CustomTypes
using CompStat
using SMM
using JLD2, FileIO
# %%
using Plots
using LaTeXStrings
pyplot()
# %%
@load "output/calibrations.jld2" cal_y
# %%
cal_y.ϵ = 17/((2013-2005+1)*12)
# %%
cal = deepcopy(cal_y)
# %%
D = 1000
Δ = 0.025
# %%
shocks = Float64[]
states = Int64[]
simulate!(shocks, states, cal_y, D)
# %%
t = Float64[]
sim_np = Float64[]
sim_nf = Float64[]
sim_Wp = Float64[]
# %%
ss = RegUncSS()
ss.flag = true
fill_RegUncSS!(ss, Δ, cal_y, [log(2.2), log(2.2), 0.93, 0.93], false)
# %%
graph!(t, sim_np, sim_nf, sim_Wp, shocks, states, ss, cal_y, D, k=100)
# %%
p1 = plot(shocks,states,lab="",color=:black,line=(:solid, 1),linetype=:steppost)
yticks!([1.,2.],[L"F_1",L"F_2"])
ylabel!(L"F")
p2 = plot(t,sim_np,lab="",color=:black,line=(:solid, 1))
ylabel!(L"n^p")
p3 = plot(t,sim_nf,lab="",color=:black,line=(:solid, 1))
ylabel!(L"n^f")
p4 = plot(t,sim_Wp,lab="",color=:black,line=(:solid, 1))
ylabel!(L"W^p")
xlabel!("Time (months)")
# %%
p = plot(p1,p2,p3, layout=(3,1))
xgrid!(:off)
ygrid!(:off)
# %%
savefig("figures/risk_pol_sim.pdf")


text:
#=
This file is part of the replication code for: Hasenzagl, T., Pellegrino, F., Reichlin, L., & Ricco, G. (2020). A Model of the Fed's View on Inflation.
Please cite the paper if you are using any part of the code for academic work (including, but not limited to, conference and peer-reviewed papers).
=#

function get_mwg_jump(Σ::Array{Float64, 2}, mwg_const::Float64, nDraws_worker::Int64)

# ----------------------------------------------------------------------------------------------------------------------
# Get the Metropolis-Within-Gibbs jumps for the proposal distribution
# ----------------------------------------------------------------------------------------------------------------------

     len_Σ              = size(Σ)[1];
     jumps_obj          = MvNormal(zeros(len_Σ), (mwg_const^2)*Σ);
     θ_unb_jumps_worker = zeros(len_Σ, nDraws_worker);

     for draw=1:nDraws_worker
          θ_unb_jumps_worker[:, draw] = rand(jumps_obj);
     end

     return θ_unb_jumps_worker;
end


text:
module Knum

using ForwardDiff
using LinearAlgebra

module Point

function bisection(f,a,b,n=100,tol=1e-7)
    p = 0.
    FA = f(a)
    for i in 1:n
        p = a+(b-a)/2
        FP = f(p)
        println(i,'\t',p)
        if FP == 0 || (b-a)/2 < tol
            return p
        end
        if sign(FA) * sign(FP) >0
            a = p
            FA = FP
        else
            b = p
        end
    end
    println("Method failed after iterations = ",n)
    return p
end

function fixed_point(g, p, N = 20, tol = 1e-7, show = false)
    #=
    Input:
    equa: The function you are going to compute
    p:    First point
    N:    Max # of iteration
    show: Whether or not to show calculation process

    Output:
    Approximation
    =#

    for i in 1:N+1
        if show
            println(i-1,'\t',p)
        end
        if abs(p-g(p))<tol
            return g(p)
        end
        p = g(p)
    end
    println("Method Failed After", n, "iterations, the result is returned.")
    return p
end

function newton(f, p, N = 20, tol = 1e-7, show = false)
    #=
    ---NOTICE---
    For Newton\'s method, your function need to enable 2 inputs,
    one with Float64, the other one with Vector.

    Input:
    equa: The function you are going to compute
    p:    First point
    N:    Max # of iteration
    tol:  Tolerance
    show: Whether or not to show calculation process

    Output:
    Approximation
    =#
    for i in 1:N+1
        println(i-1,'\t',p)
        dfp = ForwardDiff.gradient(equa,[p])[1]
        if abs(equa(p)/dfp) < tol
            if show
                println(i,'\t', p - equa(p)/dfp)
            end
            return p - f(p)/dfp
        else
            p = p - f(p)/dfp
        end
    end
    println("Method Failed After", n, "iterations, the result is returned.")
    return p
end

function secant(equa,p0,p1,N=20, tol = 1e-7 , show = false)
    #=
    Input:
    equa:  The function you are going to compute
    p0,p1: First  and second point
    N:     Max # of iteration
    tol:   Tolerance
    show:  Whether or not to show calculation process

    Output:
    Approximation

    =#
    if show
        println(0,'\t',p0)
    end
    for i in 1:N-1
        if show
            println(i,'\t',p1)
        end
        if abs(p1-p0) < tol
            return p1
        end
        dumy = p1
        p1 = p1 - (equa(p1)*(p1-p0))/(equa(p1)-equa(p0))
        p0 = dumy
    end
    if show
        println(N,'\t',p1)
    end
    return p1
end

function fpst(equa,p0,p1,N, tol = 1e-7, show = false)
    #=
    Input:
    equa:  The function you are going to compute
    p0,p1: First  and second point
    N:     Max # of iteration
    tol:   Tolerance
    show:  Whether or not to show calculation process

    Output:
    Approximation
    =#
    for i in 1:N+1
        if show
            println(i-1,'\t',p1)
        end
        dfp = equa(p1) * (p1 - p0) / (equa(p1) - equa(p0))
        if abs(dfp) < tol
            if show
                println(i,'\t', p1 - dfp)
            end
            return p1 - dfp
        else
            dumy = p1
            p1 = p1 - dfp
            p0 = dumy
        end
    end
    return p1
end

function Aitkens_delta(g, p, n)
    #=
    Input
    g - function
    p - First Point
    n - # of the sequence wanted

    Output
    aitken sequence
    =#

    # First compute the original sequence
    original_value = zeros(n+2)
    original_value[1] = g(p)
    for i in 1:n+1
        original_value[i+1] = g(original_value[i])
    end

    aitken_value = zeros(n)
    for i in 1:n
        aitken_value[i] = original_value[i] - ((original_value[i+1]-original_value[i])^2)/(original_value[i+2]-2*original_value[i+1]+original_value[i])
    end

    return aitken_value
end

function Steffensens(g, p0, n=10, tol=1e-7,show = false)
    #=
    Input
    g - function
    p - First Point
    n - # of the iteration
    tol - Tolerance

    Output
    Approximated Number
    =#
    for i in 1:n
        p1 = g(p0)
        p2 = g(p1)
        p = p0 - (p1-p0)^2/(p2-2p1+p0)
        if show
            println(i,'\t',p)
        end
        if abs(p - p0)<tol
            return p
        end
        p0 = p
    end
    println("Method Failed After", n, "iterations, the result is returned.")
    return p
end

function Muller(f, p0,p1,p2,n,tol, show = false)
    #=
    Input:
    p0,p1,p2 - Initial three guesses
    n - max # of iterations
    tol - tolerance

    Output:
    Roots
    =#

    # First make every input Complex
    p0 = Complex(p0)
    p1 = Complex(p1)
    p2 = Complex(p2)

    h1 = p1 - p0
    h2 = p2 - p1
    delta1 = (f(p1) - f(p0))/h1
    delta2 = (f(p2) - f(p1))/h2
    d = (delta2 - delta1)/(h1 + h2)

    # Create some temps
    b = Complex(0.)
    D = Complex(0.)
    E = Complex(0.)
    p = Complex(0.)

    # Start iteration
    for i in 3:n
        b = delta2 + h2 * d
        D = sqrt(b^2 - 4*f(p2)*d)
        if abs(b - D) < abs(b+D)
            E = b + D
        else
            E = b - D
        end
        h = -2* f(p2) / E
        p = p2 + h
        if abs(h)<tol
            return p
        end
        if show
            println(i,'\t',p)
        end

        # Prepare for next iteration
        p0 = p1
        p1 = p2
        p2 = p
        h1 = p1 - p0
        h2 = p2 - p1
        delta1 = (f(p1) - f(p0))/h1
        delta2 = (f(p2) - f(p1))/h2
        d = (delta2 - delta1)/(h1 + h2)
    end
    print("Method failed after N0 iterations, N0 =", n)
    return p
end

end # end module Point

module Poly

function Neville(f,x,xi::Array)
    #=
    Input:
    f - function to value
    x - the initial point
    xi - Array of xs to evaluate

    Output:
    Q
    =#

    n = length(xi)-1
    Q = zeros(n+1,n+1)
    for i in 1:n+1
        Q[i,1] = f(xi[i])
    end
    for i in 1:n
        for j in 1:i
            Q[i+1,j+1] = ((x-xi[i-j+1])*Q[i+1,j] - (x-xi[i+1])*Q[i,j])/(xi[i+1]-xi[i-j+1])
        end
    end
    return Q
end

end # module Poly

module Diff

    function SimpleDiff(f,x,h)
        #=
        Input:
        f - function to evaluate
        x - point to evaluate
        h - step
        Output:
        Compute Foward/Backward Difference
        =#
        return (f(x+h)-f(x))/h
    end

    function ThreePointEnd(f,x,h)
        return (1/(2*h))*(-3*f(x)+4*f(x+h)-f(x+2h))
    end

    function ThreePointMid(f,x,h)
        return (1/(2*h))*(f(x+h)-f(x-h))
    end

    function FivePointMid(f,x,h)
        return (1/(12*h))*(f(x-2*h)-8*f(x-h)+8*(x+h)-f(x+2*h))
    end

    function FivePointEnd(f,x,h)
        return (1/(12*h))*(-25*f(x)+48*f(x-h)-36*f(x+2*h)+16f(x+3h)-3*f(x+4h))
    end

    function SecondMid(f,x,h)
        #=
        This function computes the second derivative using Midpoint method
        =#
        return 1/(h^2)*(f(x-h)-2*f(x)+f(x+h))
    end
end # Module Diff

module Integral
    function Legendre_coef(n)
        if n == 2
            return [0.5773502692, -0.5773502692],[1.,1.]
        elseif n == 3
            return [0.7745966692,0.,-0.7745966692],[0.5555555556, 0.8888888889, 0.5555555556]
        elseif n == 4
            return [0.8611363116, 0.3399810436, -0.3399810436, -0.8611363116],[0.3478548451, 0.6521451549, 0.6521451549, 0.3478548451]
        elseif n == 5
            return [0.9061798459 0.5384693101 0.0000000000 -0.5384693101 -0.9061798459], [0.2369268850 0.4786286705 0.5688888889 0.4786286705 0.2369268850]
        end
    end

    function Trapezoidal(f,a,b)
        h = b-a
        return h/2*(f(a)+f(b))
    end

    function Simpsons(f,a,b)
        h=(b-a)/2
        return h/3*(f(a)+4*f(a+h)+f(b))
    end

    function CSimpsons(f,a,b,n)
        if n%2 == 0
            h = (b-a)/n
            odd = even = 0
            for j in 1:(n/2-1)
                odd += f(a+2*j*h)
                even += f(a+(2*j-1)*h)
            end
            even += f(a+(n-1)*h)
            return h/3*((f(a)+2*odd+4*even+f(b)))
        else
            Error("n should be even!")
        end
    end

    function CTrapezoidal(f,a,b,n)
        sum = 0
        h = (b-a)/n
        for j in 1:(n-1)
            x = a+j*h
            sum+=2*f(x)
        end
        return h/2*(sum+f(a)+f(b))
    end

    function CMidpoint(f,a,b,n)
        h = (b-a)/(n+2)
        sum = 0
        for j in 0:n/2
            sum += f(a+(2j+1)*h)
        end
        return sum*2*h
    end

    function Romberg(f,a,b,n)
        R = zeros(n,n)
        h = b-a
        R[1,1] = h/2*(f(a)+f(b))

        for row_number in 2:n
            sum = 0.
            for k in 1:2^(row_number-2)
                sum+=f(a+(k-0.5)*h)
            end
            R[row_number,1] = 0.5*(R[row_number-1,1]+h*sum)
            for j in 2:row_number
                R[row_number,j] = R[row_number,j-1]+(R[row_number,j-1]-R[row_number-1,j-1])/(4^(j-1)-1)
            end
            h = h/2
        end
        return R
    end

    function Adapt(f,a,b,tol,pts,n,return_detail=false)
        h = (b-a)/2
        x = [a:(b-a)/4:b;]
        fx = f.(x)
        n1 = n+5
        s = [1,0,4,0,1]' * fx * h/3
        s1 =[1,4,2,4,1]' * fx * h/6
        e1 = (s1-s)/15
        if abs(e1)<tol*2/3
            I = s1 + e1
            pts1 = [pts...,b]
        else
            I1,pts2,n2 = adapt(f,a,0.5*(a+b),0.5*tol,pts,n1,true)
            I2,pts1,n1 = adapt(f,0.5*(a+b),b,0.5*tol,pts2,n2,true)
            I = I1+I2
        end
        return_detail && return I,pts1,n1
        return I
    end

    function Gaussian(f,a,b,n)
        if n == 2
            tp = sqrt(3.)/3.
            tn = -sqrt(3.)/3.
            po = f((b+a)/2+tp*(b-a)/2)
            ne = f((b+a)/2+tn*(b-a)/2)
            return (po+ne)*(b-a)/2
        end
        if n == 5
            coe = [0.9061798459 0.5384693101 0.0000000000 -0.5384693101 -0.9061798459]
            c = [0.2369268850 0.4786286705 0.5688888889 0.4786286705 0.2369268850]
            trans_coe = ((b-a) .* coe .+ (b+a)) ./ 2
            ans = f.(trans_coe) * c' *(b-a)/2
            return ans[1]
        end
    end

    function SimpsonsDouble(f::Function,a::Float64,b::Float64,c::Function,d::Function,m::Int64,n::Int64)
        h = (b-a)/n
        J1 = 0
        J2 = 0
        J3 = 0
        for i in 0:n
            x = a+i*h
            HX = (d(x) - c(x))/m
            K1 = f(x,c(x)) + f(x,d(x))
            K2 = 0
            K3 = 0
            for j in 1:(m-1)
                y = c(x)+j*HX
                Q = f(x,y)
                if j%2 == 0
                    K2 = K2 + Q
                else
                    K3 = K3 + Q
                end
            end
            L = (K1 + 2*K2 + 4*K3)*HX/3
            if i == 0 || i == n
                J1 = J1 + L
            elseif i % 2 == 0
                J2 = J2 + L
            else
                J3 = J3 + L
            end # end if
        end # end for
        J = h*(J1 + 2*J2 + 4*J3)/3
        return J
    end  # end SimpsonsDouble

    function GaussianDouble(f::Function,a::Float64,b::Float64,c::Function,d::Function,m::Int64,n::Int64)
        rm,cm = Legendre_coef(m)
        rn,cn = Legendre_coef(n)
        h1 = (b-a)/2
        h2 = (b+a)/2
        J = 0
        for i in 1:m
            JX = 0
            x = h1 * rm[i] + h2
            d1 = d(x)
            c1 = c(x)
            k1 = (d1 - c1)/2
            k2 = (d1 + c1)/2
            for j in 1:n
                y = k1*rn[j]+k2
                Q = f(x,y)
                JX = JX + cn[j]*Q
            end
            J = J + cm[i]*k1*JX
        end # end for i
        J = h1*J
        return J
    end # end of GaussianDouble

end # module Integral

module ODE

    function Euler(f,a,b,alpha,N)
        h = (b-a)/N
        t = a
        w = alpha
        @show (t,w)
        for i in 1:N
            w = w+h*f(t,w)
            t = t + h
            @show (t,w)
        end
        return w
    end # Euler

    function mEuler(f,a,b,alpha,N)
        h = (b-a)/N
        t = a
        w = alpha
        for i in 1:N
            w = w+h/2*(f(t,w)+f(t+h,w+h*f(t,w)))
            t = t+h
        end
        return w
    end # mEuler

    function Taylor2(f,df,a::Float64,b::Float64,alpha::Float64,N)
        h = (b-a)/N
        t = a
        w = alpha
        @show (t,w)
        for j = 1:N
            w = w + h*f(t,w) + 0.5*h*h*df(t,w)
            t = t+h
            @show (t,w)
        end
        return w
    end # Taylor2

    function Taylor4(f,df,ddf,d3f,a,b,alpha,N)
        h = (b-a)/N
        t = a
        w = alpha
        h2 = h*h
        h3 = h*h2
        h4 = h2*h2
        @show t,w
        for j in 1:N
            w = w + h*f(t,w) + 0.5*h2*df(t,w)+h3/6*ddf(t,w)+h4/24*d3f(t,w)
            t = t+h
            @show t,w
        end
        return w
    end # End Taylor4

    function Midpoint(f,a,b,alpha,N)
        w = alpha
        h = (b-a)/N
        t = a
        @show (t,w)
        for i in 1:N
            w = w + h * f(t+h/2,w+h/2*f(t,w))
            t = t + h
            @show (t,w)
        end
        return w
    end # Taylor

    function mEuler(f,a,b,alpha,N)
        h = (b-a)/N
        t = a
        w = alpha
        @show (t,w)
        for i in 1:N
            w = w+h/2*(f(t,w)+f(t+h,w+h*f(t,w)))
            t = t+h
            @show (t,w)
        end
        return w
    end # mEuler

    function Heun(f,a,b,alpha,N)
        h = (b-a)/N
        t = a
        w = alpha
        @show (t,w)
        for i in 1:N
            w = w+h/4*(f(t,w)+3*f(t+2*h/3,w+2*h/3*f(t+h/3,w+h/3*f(t,w))))
            t = t+h
            @show (t,w)
        end
        return w
    end # Heun

    function RK4(f,a,b,alpha,N)
        h = (b-a)/N
        t = a
        w = alpha
        @show (t,w)

        for i in 1:N
            k1 = h*f(t,w)
            k2 = h*f(t+h/2,w+1/2*k1)
            k3 = h*f(t+h/2,w+1/2*k2)
            k4 = h*f(t+h,w+k3)
            w = w+(k1+2*k2+2*k3+k4)/6
            t = t+h
            @show (t,w)
        end
        return w
    end # RK4

    function  AB3(f,a,b,alpha,N)
        h = (b-a)/N
        t = a
        w = zeros(N+1)
        w[1] = alpha
        for i in 2:3
            k1 = h*f(t,w[i-1])
            k2 = h*f(t+h/2,w[i-1]+1/2*k1)
            k3 = h*f(t+h/2,w[i-1]+1/2*k2)
            k4 = h*f(t+h,w[i-1]+k3)
            w[i] = w[i-1]+(k1+2*k2+2*k3+k4)/6
            t = t+h
        end
        for i in 4:N+1
            w[i] = w[i-1] + h/12 * (23*f(t,w[i-1])- 16*f(t-h,w[i-2])+ 5*f(t-2*h,w[i-3]))
            t = t + h
        end
        return w
    end # end AB3

    function AM3(f,a,b,alpha,N)
        h = (b-a)/N
        t = a
        w = zeros(N+1)
        w[1] = alpha
        for i in 2:3
            @show t,w
            k1 = h*f(t,w[i-1])
            k2 = h*f(t+h/2,w[i-1]+1/2*k1)
            k3 = h*f(t+h/2,w[i-1]+1/2*k2)
            k4 = h*f(t+h,w[i-1]+k3)
            w[i] = w[i-1]+(k1+2*k2+2*k3+k4)/6
            t = t + h
        end
        for i in 4:N+1
            @show t,w
            w[i] = w[i-1] + h/12 * (23*f(t,w[i-1]) - 16*f(t-h,w[i-2])+ 5*f(t-2*h,w[i-3]))
            w[i] = w[i-1] + h/24 * (9*f(t+h,w[i]) + 19*f(t,w[i-1]) - 5*f(t-h,w[i-2]) + f(t-2*h,w[i-3]))
            t = t + h
        end
        @show t,w
        return w
    end # AM3

    function rk4s(f,a,b,alpha,N)
        #=
        alpha is the initial values column vector
        f is the array of functions we use in column vector, m functions
        =#

        d = length(alpha)
        h = (b-a)/N
        t = a .+ [0:N;] .* h
        w = zeros(Float64,N+1,d)
        w[1,:] = alpha
        for j in 1:N
            k1 = f(t[j],w[j,:])
            k2 = f(t[j] .+ 0.5*h, w[j,:] .+ 0.5 * h .* k1)
            k3 = f(t[j] .+ 0.5*h, w[j,:] .+ 0.5 * h .* k2)
            k4 = f(t[j] .+     h, w[j,:] .+       h .* k3)
            w[j+1,:] = w[j,:] .+ h/6 * (k1 + 2* k2 + 2*k3 + k4)
            @show t[j], w[j+1,:]
        end
        return w,t
    end #End Rk4s

    function Trapezoidal(f,df,a,b,alpha,N,TOL,M = 1000)
        #=
        M is maximum number of iterations
        =#
        h = (b-a)/N
        t = a
        w = alpha

        for i in 1:N
            k1 = w+h/2*f(t,w)
            w0 = k1
            j = 1
            FLAG = 0
            while FLAG == 0
                w = w0- (w0 - h/2*f(t+h,w0)-k1)/(1-h/2*df(t+h,w0))
                if abs(w-w0) < TOL
                    FLAG = 1
                else
                    j = j+1
                    w0 = w
                    if j > M
                        println("The maximum number of iterations exceeded.")
                        return t,w
                    end # end if
                end # end if
            end # end while
            t = a+i*h
        end # end for
        return t,w
    end # end Trapezoidal

    function BkwdEulerStiff(lambda,a,b,alpha,N)
        h = (b-a)/N
        t = a
        w = alpha
        for i in 1:N
            w = w/(1-h*lambda)
        end
        return w
    end

end # End ODE

module Matrix
    using LinearAlgebra

    function LU(A)
        n,m = size(A)
        L = Array{Float64}(I,n,n)
        U = Array{Float64}(I,n,n)

        # Step 1
        U[1,1] = A[1,1]/L[1,1]
        if L[1,1]*U[1,1] == 0
            throw(error("Factorization impossible"))
        end

        # Step 2
        for j in 2:n
            U[1,j] = A[1,j]/L[1,1]
            L[j,1] = A[j,1]/U[1,1]
        end

        # Step 3
        for i in 2:(n-1)

            # Step 4
            temp = 0
            for k in 1:(i-1)
                temp = temp + L[i,k]*U[k,i]
            end
            U[i,i] = (A[i,i] - temp)/L[i,i]

            L[1,1]*U[1,1] == 0 && throw(error("Factorization impossible"))

            # Step 5
            for j in (i+1):n
                temp1 = 0
                temp2 = 0
                for k in 1:(i-1)
                    temp1 = temp1 + L[i,k]*U[k,j]
                    temp2 = temp2 + L[j,k] * U[k,i]
                end # k
                U[i,j] = 1/L[i,i]*(A[i,j] - temp1)
                L[j,i] = 1/U[i,i]*(A[j,i] - temp2)
            end # end for j
        end # end for i

        # Step 6
        temp3 = 0
        for k in 1:(n-1)
            temp3 =temp3 + L[n,k] * U[k,n]
        end
        U[n,n] = (A[n,n] - temp3) / L[n,n]

        return L,U
    end #LU

end # modeul Matrix

end # module Knum

max_stars_repo_name: chenyangzhu/Knum

max_issues_repo_name: chenyangzhu/Knum

max_forks_repo_name: chenyangzhu/Knum


text:
module SSLTestASCII
using SimpleSDMLayers
using Test

M = rand(Float64, 4, 8)
M = convert(Matrix{Union{Nothing,Float64}}, M)
M[rand(Bool, size(M))] .= nothing

S = SimpleSDMPredictor(M, 0.0, 2.0, 0.0, 1.0)

SimpleSDMLayers.ascii(S, "test.asc")
@test isfile("test.asc")

U = SimpleSDMLayers.ascii("test.asc")

@test all(S.grid .== U.grid)
@test S.left == U.left
@test S.bottom == U.bottom
@test S.right == U.right
@test S.top == U.top
@test size(S) == size(U)

rm("test.asc")

end


text:
function SQ._collect{T<:AbstractTable}(
    g_tbl::Grouped{T}, q::SQ.SelectNode
)::Grouped
    src = g_tbl.source
    new_src = default(src)
    groupby_fields = map(
        x -> isa(x, Symbol) ? x : g_tbl.predicate_aliases[x],
        g_tbl.groupbys
    )
    # Copy all columns involved in the original grouping
    for groupby_field in groupby_fields
        new_src[groupby_field] = src[groupby_field]
    end
    # Apply the select query
    for h in q.helpers
        apply!(new_src, h, src)
    end
    return Grouped(
        new_src,
        g_tbl.group_indices,
        g_tbl.group_levels,
        g_tbl.groupbys,
        g_tbl.predicate_aliases
    )
end


text:
##
input_range = 158126:624574
data = digits.(input_range) .|> reverse
length(data)

##
function filter_passw(data)
    nondec = filter(data) do xs
        d = first(xs)
        for i in xs[2:end]
            d <= i || return false
            d = i
        end
        true
    end

    filter(nondec) do xs
        d = first(xs)
        groups = Dict{Int,Int}()
        for i in xs[2:end]
            if i == d
                if d in keys(groups)
                    groups[d] += 1
                else
                    groups[d] = 2
                end
            end
            d = i
        end
        any(==(2), values(groups))
    end
end
count_many(data) = length(data)

##
filter_passw(data)


text:
using RecombinatorKMeans
using Test
using DelimitedFiles

# The a3 dataset was downloaded from here:
#
#  http://cs.uef.fi/sipu/datasets/
#
#  Clustering basic benchmark
#  P. Fränti and S. Sieranoja
#  K-means properties on six clustering benchmark datasets
#  Applied Intelligence, 48 (12), 4743-4759, December 2018
#  https://doi.org/10.1007/s10489-018-1238-7

a3 = Matrix(readdlm(joinpath(@__DIR__, "a3.txt"))')
a3 ./= maximum(a3)
n, m = size(a3)
k = 50

@testset "kmeans uniform" begin
    labels, centroids, cost, converged = kmeans(a3, k, init="unif", verbose=false)
    @test length(labels) == m
    @test all(∈(1:k), labels)
    @test size(centroids) == (2,k)
    @test 7 < cost < 25
    @test converged == true
end

@testset "kmeans++" begin
    labels, centroids, cost, converged = kmeans(a3, k, init="++", verbose=false)
    @test length(labels) == m
    @test all(∈(1:k), labels)
    @test size(centroids) == (2,k)
    @test 6.7 < cost < 11
    @test converged == true
end

@testset "reckmeans" begin
    res = reckmeans(a3, k, 5, Δβ = 0.1, verbose=false)
    @test res.exit_status == :collapsed
    @test length(res.labels) == m
    @test all(∈(1:k), res.labels)
    @test size(res.centroids) == (2,k)
    @test 6.7 < res.cost < 7.5
end

@testset "gakmeans" begin
    res = gakmeans(a3, k, 5, verbose=false)
    @test res.exit_status == :collapsed
    @test length(res.labels) == m
    @test all(∈(1:k), res.labels)
    @test size(res.centroids) == (2,k)
    @test 6.7 < res.cost < 7.5
end

@testset "rswapkmeans uniform" begin
    res = rswapkmeans(a3, k, max_time=1.0, max_swaps=typemax(Int), verbose=false)
    @test res.exit_status == :outoftime
    @test length(res.labels) == m
    @test all(∈(1:k), res.labels)
    @test size(res.centroids) == (2,k)
    @test 6.7 < res.cost < 7.5

    res = rswapkmeans(a3, k, max_time=Inf, max_swaps=200, verbose=false)
    @test res.exit_status == :maxswaps
    @test length(res.labels) == m
    @test all(∈(1:k), res.labels)
    @test size(res.centroids) == (2,k)
    @test 6.7 < res.cost < 7.5

    res = rswapkmeans(a3, k, max_time=Inf, max_swaps=typemax(Int), target_cost=6.75, verbose=false)
    @test res.exit_status == :solved
    @test length(res.labels) == m
    @test all(∈(1:k), res.labels)
    @test size(res.centroids) == (2,k)
    @test 6.7 < res.cost ≤ 6.75

    res = rswapkmeans(a3, k, max_time=1.0, max_swaps=typemax(Int), verbose=false, init="unif")
    @test res.exit_status == :outoftime
    @test length(res.labels) == m
    @test all(∈(1:k), res.labels)
    @test size(res.centroids) == (2,k)
    @test 6.7 < res.cost < 11

end


text:
#
# cropping.jl -
#
# Provide zero-padding and cropping operators.
#
#-------------------------------------------------------------------------------
#
# This file is part of LazyAlgebra (https://github.com/emmt/LazyAlgebra.jl)
# released under the MIT "Expat" license.
#
# Copyright (c) 2019-2021, Éric Thiébaut.
#

module Cropping

# FIXME: add simplifying rules:
#   Z'*Z = Id (not Z*Z' = Id)  crop zero-padded array is identity

export
    CroppingOperator,
    ZeroPaddingOperator,
    defaultoffset

using ArrayTools
using ..Foundations
using ..LazyAlgebra
using ..LazyAlgebra: bad_argument, bad_size
import ..LazyAlgebra: apply!, vcreate,
    input_size, input_ndims, output_size, output_ndims

"""
    CroppingOperator(outdims, inpdims, offset=defaultoffset(outdims,inpdims))

yields a linear map which implements cropping of arrays of size `inpdims` to
produce arrays of size `outdims`.  By default, the output array is centered
with respect to the inpput array (using the same conventions as `fftshift`).
Optional argument `offset` can be used to specify a different relative
position.  If `offset` is given, the output value at multi-dimensional index
`i` is given by input value at index `j = i + offset`.

The adjoint of a cropping operator is a zero-padding operator.

See also: [`ZeroPaddingOperator`](@ref).

"""
struct CroppingOperator{N} <: LinearMapping
    outdims::NTuple{N,Int} # cropped dimensions
    inpdims::NTuple{N,Int} # input dimensions
    offset::CartesianIndex{N} # offset of cropped region w.r.t. input array
    function CroppingOperator{N}(outdims::NTuple{N,Int},
                                 inpdims::NTuple{N,Int}) where {N}
        @inbounds for d in 1:N
            1 ≤ outdims[d] || error("invalid output dimension(s)")
            outdims[d] ≤ inpdims[d] ||
                error(1 ≤ inpdims[d]
                      ? "invalid input dimension(s)"
                      : "output dimensions must be less or equal input ones")
        end
        offset = defaultoffset(inpdims, outdims)
        return new{N}(outdims, inpdims, offset)
    end
    function CroppingOperator{N}(outdims::NTuple{N,Int},
                                 inpdims::NTuple{N,Int},
                                 offset::CartesianIndex{N}) where {N}
        @inbounds for d in 1:N
            1 ≤ outdims[d] || error("invalid output dimension(s)")
            outdims[d] ≤ inpdims[d] ||
                error(1 ≤ inpdims[d]
                      ? "invalid input dimension(s)"
                      : "output dimensions must less or equal input ones")
            0 ≤ offset[d] ≤ inpdims[d] - outdims[d] ||
                error("out of range offset(s)")
        end
        return new{N}(outdims, inpdims, offset)
    end
end

@callable CroppingOperator

commonpart(C::CroppingOperator) = CartesianIndices(output_size(C))
offset(C::CroppingOperator) = C.offset

input_ndims(C::CroppingOperator{N}) where {N} = N
input_size(C::CroppingOperator) = C.inpdims
input_size(C::CroppingOperator, i...) = input_size(C)[i...]

output_ndims(C::CroppingOperator{N}) where {N} = N
output_size(C::CroppingOperator) = C.outdims
output_size(C::CroppingOperator, i...) = output_size(C)[i...]

# Union of acceptable types for the offset.
const Offset = Union{CartesianIndex,Integer,Tuple{Vararg{Integer}}}

CroppingOperator(outdims::ArraySize, inpdims::ArraySize) =
    CroppingOperator(to_size(outdims), to_size(inpdims))

CroppingOperator(outdims::ArraySize, inpdims::ArraySize, offset::Offset) =
    CroppingOperator(to_size(outdims), to_size(inpdims),
                     CartesianIndex(offset))

CroppingOperator(::Tuple{Vararg{Int}}, ::Tuple{Vararg{Int}}) =
    error("numbers of output and input dimensions must be equal")

CroppingOperator(::Tuple{Vararg{Int}}, ::Tuple{Vararg{Int}}, ::CartesianIndex) =
    error("numbers of output and input dimensions and offsets must be equal")

CroppingOperator(outdims::NTuple{N,Int}, inpdims::NTuple{N,Int}) where {N} =
    CroppingOperator{N}(outdims, inpdims)

CroppingOperator(outdims::NTuple{N,Int}, inpdims::NTuple{N,Int},
                 offset::CartesianIndex{N}) where {N} =
    CroppingOperator{N}(outdims, inpdims, offset)

function vcreate(::Type{Direct},
                 C::CroppingOperator{N},
                 x::AbstractArray{T,N},
                 scratch::Bool) where {T,N}
    (scratch && isa(x, Array{T,N}) && input_size(C) == output_size(C)) ? x :
        Array{T,N}(undef, output_size(C))
end

function vcreate(::Type{Adjoint},
                 C::CroppingOperator{N},
                 x::AbstractArray{T,N},
                 scratch::Bool) where {T,N}
    (scratch && isa(x, Array{T,N}) && input_size(C) == output_size(C)) ? x :
        Array{T,N}(undef, input_size(C))
end

# Apply cropping operation.
#
#     for I in R
#         J = I + K
#         y[I] = α*x[J] + β*y[I]
#     end
#
function apply!(α::Number,
                ::Type{Direct},
                C::CroppingOperator{N},
                x::AbstractArray{T,N},
                scratch::Bool,
                β::Number,
                y::AbstractArray{T,N}) where {T,N}
    has_standard_indexing(x) ||
        bad_argument("input array has non-standard indexing")
    size(x) == input_size(C) ||
        bad_size("bad input array dimensions")
    has_standard_indexing(y) ||
        bad_argument("output array has non-standard indexing")
    size(y) == output_size(C) ||
        bad_size("bad output array dimensions")
    if α == 0
        β == 1 || vscale!(y, β)
    else
        k = offset(C)
        I = commonpart(C)
        if α == 1
            if β == 0
                @inbounds @simd for i in I
                    y[i] = x[i + k]
                end
            elseif β == 1
                @inbounds @simd for i in I
                    y[i] += x[i + k]
                end
            else
                beta = convert(T, β)
                @inbounds @simd for i in I
                    y[i] = x[i + k] + beta*y[i]
                end
            end
        else
            alpha = convert(T, α)
            if β == 0
                @inbounds @simd for i in I
                    y[i] = alpha*x[i + k]
                end
            elseif β == 1
                @inbounds @simd for i in I
                    y[i] += alpha*x[i + k]
                end
            else
                beta = convert(T, β)
                @inbounds @simd for i in I
                    y[i] = alpha*x[i + k] + beta*y[i]
                end
            end
        end
    end
    return y
end

# Apply zero-padding operation.
#
#     for i in I
#         y[i + k] = α*x[i] + β*y[i + k]
#     end
#     # Plus y[i + k] *= β outside common region R
#
function apply!(α::Number,
                ::Type{Adjoint},
                C::CroppingOperator{N},
                x::AbstractArray{T,N},
                scratch::Bool,
                β::Number,
                y::AbstractArray{T,N}) where {T,N}
    has_standard_indexing(x) ||
        bad_argument("input array has non-standard indexing")
    size(x) == output_size(C) ||
        bad_size("bad input array dimensions")
    has_standard_indexing(y) ||
        bad_argument("output array has non-standard indexing")
    size(y) == input_size(C) ||
        bad_size("bad output array dimensions")
    β == 1 || vscale!(y, β)
    if α != 0
        k = offset(C)
        I = commonpart(C)
        if α == 1
            if β == 0
                @inbounds @simd for i in I
                    y[i + k] = x[i]
                end
            else
                @inbounds @simd for i in I
                    y[i + k] += x[i]
                end
            end
        else
            alpha = convert(T, α)
            if β == 0
                @inbounds @simd for i in I
                    y[i + k] = alpha*x[i]
                end
            else
                @inbounds @simd for i in I
                    y[i + k] += alpha*x[i]
                end
            end
        end
    end
    return y
end

"""
    ZeroPaddingOperator(outdims, inpdims, offset=defaultoffset(outdims,inpdims))

yields a linear map which implements zero-padding of arrays of size `inpdims`
to produce arrays of size `outdims`.  By default, the input array is centered
with respect to the output array (using the same conventions as `fftshift`).
Optional argument `offset` can be used to specify a different relative
position.  If `offset` is given, the input value at multi-dimensional index `j`
is copied at index `i = j + offset` in the result.

A zero-padding operator is implemented as the adjoint of a cropping operator.

See also: [`CroppingOperator`](@ref).

"""
ZeroPaddingOperator(outdims, inpdims) =
    Adjoint(CroppingOperator(inpdims, outdims))
ZeroPaddingOperator(outdims, inpdims, offset) =
    Adjoint(CroppingOperator(inpdims, outdims, offset))

"""
    defaultoffset(dim1,dim2)

yields the index offset such that the centers (in the same sense as assumed by
`fftshift`) of dimensions of lengths `dim1` and `dim2` are coincident.  If `off
= defaultoffset(dim1,dim2)` and `i2` is the index along `dim2`, then the index
along `dim1` is `i1 = i2 + off`.

"""
defaultoffset(dim1::Integer, dim2::Integer) =
    (Int(dim1) >> 1) - (Int(dim2) >> 1)
defaultoffset(dims1::NTuple{N,Integer}, dims2::NTuple{N,Integer}) where {N} =
    CartesianIndex(map(defaultoffset, dims1, dims2))

end # module