text:
# This file is part of EndemicModel project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

@defgroup world2(world) (country, date) [
    (latitude, longitude) => mean,
    (confirmed, recovered, deaths) => sum,
]

@deftable tests_no_duplicate(total_tests) begin
    groups = groupby(total_tests, [:country, :date])
    data = combine(groups) do subdf
        testkind = _choose_test(filter(!ismissing, subdf.test_kind))
        idx = findfirst(isequal(testkind), subdf.test_kind)
        subdf[idx, :]
    end
    data = sort!(data, [:country, :date])
end

function _choose_test(tests)
    if isempty(tests)
        missing
    elseif "people tested" ∈ tests
        "people tested"
    elseif "samples tested" ∈ tests
        "samples tested"
    elseif "units unclear" ∈ tests
        "units unclear"
    else
        @warn "Unrecognized test kind(s)." tests
        collect(tests)[1]
    end
end

@deftable world3(world2, world_population, world_gdp_per_capita, tests_no_duplicate) begin
    data = leftjoin(world2, world_population; on = :country)
    data = leftjoin(data, world_gdp_per_capita; on = [:country, :country_code])
    data = leftjoin(data, tests_no_duplicate; on = [:country, :date])
    rename!(data, :recovered => :recovered_confirmed)

    # Estimate recovered
    # Source: https://journals.lww.com/cmj/Fulltext/2020/05050/Persistence_and_clearance_of_viral_RNA_in_2019.6.aspx
    # Average infeccious period of 9.5 (6.0 to 11.0) days
    # Source: https://www.sciencedirect.com/science/article/pii/S0163445320301195
    # Average infeccious period of 11.0 (10.0 to 12.0) days
    infected_vec, countries = data.confirmed, data.country
    active, closed, last_infected = 0, 0, 0
    closed_vec = Int[]
    last_country = "VOID"

    for (infected, country) ∈ zip(infected_vec, countries)
        if country == last_country
            new_closed, new_active = active / 11.0, infected - last_infected
            closed += new_closed
            active += new_active - new_closed
            last_infected = infected
            push!(closed_vec, round(Int, closed))
        else
            active, closed, last_infected = 0, 0, 0
            last_country = country
            push!(closed_vec, 0)
        end
    end
    insertcols!(data, :recovered => closed_vec .- data.deaths)

    cols = [
        :country,
        :region,
        :income_group,
        :date,
        :latitude,
        :longitude,
        :estimated_population,
        :gdp_per_capita,
        :total_tests,
        :test_kind,
        :confirmed,
        :recovered_confirmed,
        :recovered,
        :deaths,
    ]
    select!(data, cols)
end

@defgroup per_country_aux(world3) (country, region) [
    (income_group,) => identity,
    (date, latitude, longitude) => identity,
    (estimated_population, gdp_per_capita) => identity,
    (total_tests, test_kind) => identity,
    (confirmed, recovered, deaths) => identity,
]

@defgroup per_country(per_country_aux) (country,) [
    (region, income_group) => identity,
    (date, latitude, longitude) => identity,
    (estimated_population, gdp_per_capita) => identity,
    (total_tests, test_kind) => identity,
    (confirmed, recovered, deaths) => identity,
]

@defgroup per_region_aux(world3) (region, date) [
    (estimated_population,) => sum,
    (total_tests,) => _sum ∘ skipmissing,
    (confirmed, recovered, deaths) => sum,
]

@defgroup per_region(per_region_aux) (region,) [
    (date,) => identity,
    (estimated_population,) => sum,
    (total_tests,) => _sum ∘ skipmissing,
    (confirmed, recovered, deaths) => sum,
]

@defgroup per_income_group_aux(world3) (income_group, date) [
    (estimated_population,) => sum,
    (total_tests,) => _sum ∘ skipmissing,
    (confirmed, recovered, deaths) => sum,
]

@defgroup per_income_group(per_income_group_aux) (income_group,) [
    (date,) => identity,
    (estimated_population,) => sum,
    (total_tests,) => _sum ∘ skipmissing,
    (confirmed, recovered, deaths) => sum,
]

@defgroup per_date(per_country) (date,) [
    (latitude, longitude) => mean,
    (estimated_population,) => sum,
    (total_tests,) => _sum ∘ skipmissing,
    (confirmed, recovered, deaths, infected, closed, active) => sum,
    (diff_infected, diff_recovered, diff_deaths, diff_closed, diff_active) => sum,
    (diff2_infected, diff3_infected, diff4_infected, diff5_infected) => sum,
]

@deftable per_date2(per_date, world_population, world_gdp_per_capita) begin
    wpop, wgdp = world_population, world_gdp_per_capita
    pop = wpop[wpop.country.=="World", :estimated_population]
    gdp = wgdp[wgdp.country.=="World", :gdp_per_capita]
    cols =
        [:date, :total_tests, :confirmed, :recovered, :deaths, :infected, :closed, :active]
    data = select(per_date, cols)
    insertcols!(data, 2, :estimated_population => pop[1])
    insertcols!(data, 3, :gdp_per_capita => gdp[1])
    insertcols!(data, 5, :test_kind => "units unclear")
    data
end

@deftable brazil2(brazil, total, estimates) begin
    keys = [:date, :estimated_population, :gdp_per_capita, :total_tests, :test_kind]
    br = select(brazil, keys)
    avg_ipc = estimates[estimates.state.=="total", :infected_per_confirmed][1]

    # Estimate recovered
    # Source: https://journals.lww.com/cmj/Fulltext/2020/05050/Persistence_and_clearance_of_viral_RNA_in_2019.6.aspx
    # Average infeccious period of 9.5 (6.0 to 11.0) days
    # Source: https://www.sciencedirect.com/science/article/pii/S0163445320301195
    # Average infeccious period of 11.0 (10.0 to 12.0) days
    infected_vec = round.(Int, avg_ipc * brazil.confirmed)
    active, closed, last_infected = 0, 0, 0
    closed_vec = Int[]
    for infected ∈ infected_vec
        new_closed, new_active = active / 11.0, infected - last_infected
        closed += new_closed
        active += new_active - new_closed
        last_infected = infected
        push!(closed_vec, round(Int, closed))
    end

    insertcols!(br, 1, :country => "Brazil")
    insertcols!(br, :confirmed => brazil.confirmed)
    insertcols!(br, :recovered_confirmed => brazil.recovered)
    insertcols!(br, :recovered => closed_vec .- brazil.deaths)
    insertcols!(br, :deaths => brazil.deaths)
    insertcols!(br, :infected => infected_vec)
    sort!(br, :date)
end

@deftable states2(states, brazil2) begin
    nrec = Union{Int,Missing}[]
    for row ∈ eachrow(states)
        brrow = brazil2[brazil2.date.==row.date, :]
        push!(
            nrec,
            if !isempty(brrow)
                brrec = brrow.recovered[1]
                brinfec = brrow.infected[1]
                infec = row.infected
                brinfec == 0 ? 0 : round(Union{Int,Missing}, brrec * infec / brinfec)
            else
                missing
            end,
        )
    end
    insertcols!(states, 5, :recovered => nrec)
end

@deftable cities2(cities, brazil2) begin
    nrec = Union{Int,Missing}[]
    for row ∈ eachrow(cities)
        brrow = brazil2[brazil2.date.==row.date, :]
        push!(
            nrec,
            if !isempty(brrow)
                brrec = brrow.recovered[1]
                brinfec = brrow.infected[1]
                infec = row.infected
                brinfec == 0 ? 0 : round(Union{Int,Missing}, brrec * infec / brinfec)
            else
                missing
            end,
        )
    end
    insertcols!(cities, 6, :recovered => nrec)
end

@defgroup per_state(states2) (state,) [
    (date, city_ibge_code) => identity,
    (estimated_population,) => identity,
    (confirmed, recovered, deaths, infected) => identity,
]

@defgroup per_city(cities2) (city, state) [
    (state, date, city_ibge_code) => identity,
    (estimated_population,) => identity,
    (confirmed, recovered, deaths, infected) => identity,
]

@defcolumn infected(confirmed) confirmed
@defcolumn closed(recovered, deaths) recovered .+ deaths
@defcolumn active(infected, closed) infected .- closed

function _column_diff(column, date, ndays)
    @argcheck ndays >= 1
    n = length(date)
    result = Vector{OptFloat}(missing, n)
    half2 = floor(Int, ndays / 2)
    half1 = ndays - half2
    istart = 1 + half1
    iend = n - half2
    for i ∈ istart:iend
        if ( # consistency check
            date[i-half1] + Day(ndays) == date[i] + Day(half2) == date[i+half2] &&
            !ismissing(column[i-half1]) &&
            isfinite(column[i-half1])
        )
            result[i] = (column[i+half2] - column[i-half1]) / ndays
        end
    end
    result
end

@defcolumn diff_confirmed(date, confirmed) _column_diff(confirmed, date, 7)
@defcolumn diff_infected(date, infected) _column_diff(infected, date, 7)
#= begin
    n = length(date)
    result = Vector{OptFloat}(missing, n)
    for i ∈ 5:(n - 3)
        if (
            date[i - 4] + Day(7) == date[i] + Day(3) == date[i + 3] && # consistency check
            !ismissing(infected[i - 4]) &&
            isfinite(infected[i - 4]) &&
            infected[i - 4] > 0
        )
            # Geometric mean of infected[i+1]/infected[i] over 7 days
            diff_rate = (infected[i + 3] / infected[i - 4])^(1 / 7)
            if !ismissing(diff_rate) && isfinite(diff_rate)
                result[i] = diff_rate * infected[i]
            end
        end
    end
    result
end
=#

@defcolumn diff_recovered(date, recovered) _column_diff(recovered, date, 7)
@defcolumn diff_deaths(date, deaths) _column_diff(deaths, date, 7)

@defcolumn diff_closed(diff_recovered, diff_deaths) diff_recovered .+ diff_deaths
@defcolumn diff_active(diff_infected, diff_closed) diff_infected .- diff_closed

@defcolumn diff2_infected(date, diff_infected) _column_diff(diff_infected, date, 5)
@defcolumn diff2_recovered(date, diff_recovered) _column_diff(diff_recovered, date, 5)
@defcolumn diff2_deaths(date, diff_deaths) _column_diff(diff_deaths, date, 5)

@defcolumn diff2_closed(diff2_recovered, diff2_deaths) diff2_recovered .+ diff2_deaths
@defcolumn diff2_active(diff2_infected, diff2_closed) diff2_infected .- diff2_closed

@defcolumn diff3_infected(date, diff2_infected) _column_diff(diff2_infected, date, 5)
#=
@defcolumn diff3_recovered(date, diff2_recovered) _column_diff(
    diff2_recovered,
    date,
    3,
)
@defcolumn diff3_deaths(date, diff2_deaths) _column_diff(diff2_deaths, date, 3)

@defcolumn diff3_closed(diff3_recovered, diff3_deaths) diff3_recovered .+ diff3_deaths
@defcolumn diff3_active(diff3_infected, diff3_closed) diff3_infected .- diff3_closed
=#
@defcolumn diff4_infected(date, diff3_infected) _column_diff(diff3_infected, date, 5)

@defcolumn diff5_infected(date, diff4_infected) _column_diff(diff4_infected, date, 5)

@defcolumn diff6_infected(date, diff5_infected) _column_diff(diff5_infected, date, 5)

const COVID_19_μ = 0.034

@defcolumn μ_closed(deaths, closed) begin
    replace(1.0 .* deaths ./ closed, NaN => missing)
end

@defcolumn infected_per_1mi(infected, estimated_population) begin
    1.0e6 .* infected ./ estimated_population
end
@defcolumn deaths_per_1mi(deaths, estimated_population) begin
    1.0e6 .* deaths ./ estimated_population
end
@defcolumn recovered_per_1mi(recovered, estimated_population) begin
    1.0e6 .* recovered ./ estimated_population
end
@defcolumn closed_per_1mi(closed, estimated_population) begin
    1.0e6 .* closed ./ estimated_population
end
@defcolumn active_per_1mi(active, estimated_population) begin
    1.0e6 .* active ./ estimated_population
end

@defcolumn tests_per_confirmed(total_tests, confirmed) 1.0 .* total_tests ./ confirmed
@defcolumn tests_per_infected(total_tests, infected) 1.0 .* total_tests ./ infected
@defcolumn tests_per_1mi(total_tests, estimated_population) begin
    1.0e6 .* total_tests ./ estimated_population
end

function copy_pop_tests_tables(data::AbstractDict)
    pop = get(data, :world_population, nothing)
    gdp_pc = get(data, :world_gdp_per_capita, nothing)
    tests = get(data, :tests_no_duplicate, nothing)
    csse = get(data, :csse, nothing)
    brasil_io = get(data, :brasil_io, nothing)
    updated = false
    if brasil_io != nothing
        if pop != nothing && !haskey(brasil_io, :world_population)
            brasil_io[:world_population] = pop
            updated = true
        end
    end
    if csse != nothing
        if pop != nothing && !haskey(csse, :world_population)
            csse[:world_population] = pop
            updated = true
        end
        if gdp_pc != nothing && !haskey(csse, :world_gdp_per_capita)
            csse[:world_gdp_per_capita] = gdp_pc
            updated = true
        end
        if tests != nothing && !haskey(csse, :tests_no_duplicate)
            csse[:tests_no_duplicate] = tests
            true
        else
            updated
        end
    else
        for (_, subdata) ∈ data
            if subdata isa AbstractDict && copy_pop_tests_tables(subdata)
                return true
            end
        end
        false
    end
end

function copy_brazil_tables(data::AbstractDict)
    csse = get(data, :csse, nothing)
    brasil_io = get(data, :brasil_io, nothing)
    if (
        csse != nothing &&
        brasil_io != nothing &&
        !haskey(brasil_io, :brazil) &&
        haskey(csse, :per_country_group)
    )
        brasil_io[:brazil] = csse[:per_country_group][:brazil]
        true
    else
        for (_, subdata) ∈ data
            if subdata isa AbstractDict && copy_brazil_tables(subdata)
                return true
            end
        end
        false
    end
end

_sum(itr) = isempty(itr) ? missing : sum(itr)
_maximum(itr) = isempty(itr) ? missing : maximum(itr)

function covid19_database(; kwargs...)
    sources = [:csse, :brasil_io, :world_population, :world_gdp_per_capita, :total_tests]
    funcs = [
        # Tables
        group_world2,
        table_tests_no_duplicate,
        copy_pop_tests_tables,
        table_world3,
        group_per_country_aux,
        group_per_country,
        group_per_region_aux,
        group_per_region,
        group_per_income_group_aux,
        group_per_income_group,
        copy_brazil_tables,
        table_brazil2,
        table_states2,
        table_cities2,
        group_per_state,
        group_per_city,
        # Columns
        column_infected,
        column_closed,
        column_active,
        column_diff_confirmed,
        column_diff_infected,
        column_diff_recovered,
        column_diff_deaths,
        column_diff_closed,
        column_diff_active,
        column_diff2_infected,
        column_diff2_recovered,
        column_diff2_deaths,
        column_diff2_closed,
        column_diff2_active,
        column_diff3_infected,
        column_diff4_infected,
        column_diff5_infected,
        # These tables needs diff columns before grouping
        group_per_date,
        table_per_date2,
        # More columns
        column_μ_closed,
        column_infected_per_1mi,
        column_deaths_per_1mi,
        column_recovered_per_1mi,
        column_closed_per_1mi,
        column_active_per_1mi,
        column_tests_per_confirmed,
        column_tests_per_1mi,
    ]
    Database{DataDict}(sources, kwargs, funcs; name = "COVID-19")
end

function covid19(;
    database = nothing,
    compute_model::Bool = true,
    pretty::Bool = true,
    do_export::Bool = true,
    export_trials = 3,
    kwargs...,
)
    db = database
    @info "Creating COVID-19 database."
    isnothing(database) && (db = covid19_database(; kwargs...))

    if hasproperty(db, :csse)
        sources = DataFrame(
            :source => [
                "CSSE at JHU (John Hopkins University)",
                "Brasil.IO",
                "World Bank Open Data",
                "Our World In Data",
            ],
            :url => [
                "https://systems.jhu.edu/",
                "https://brasil.io/home/",
                "https://data.worldbank.org/",
                "https://ourworldindata.org/",
            ],
            :data_kind => [
                "COVID-19 cases for several countries",
                "COVID-19 cases per state and city in Brazil",
                "Population and GDP per capita for several countries",
                "Number of tests of COVID-19 for several countries",
            ],
        )
        for row ∈ eachrow(db.brasil_io.estimates)
            data_kind =
                row.state == "total" ? "Estimate for Brazil on $(row.date)" :
                "Estimate for $(row.state), Brazil on $(row.date)"
            push!(sources, (; source = row.source, url = row.url, data_kind = data_kind))
        end

        newdata = DataDict(
            :world => DataDict(
                :per_country => db.csse.per_country_group,
                :per_region => db.csse.per_region_group,
                :per_income_group => db.csse.per_income_group_group,
                :per_date => db.csse.per_date_group,
                :total => db.csse.per_date2,
            ),
            :brazil => DataDict(
                :per_state => db.brasil_io.per_state_group,
                :per_city => db.brasil_io.per_city_group,
                :total => db.brasil_io.brazil2,
            ),
            :sources => sources,
            :auxiliar => DataDict(
                :total_tests => db.total_tests,
                :world_population => db.world_population,
                :world_gdp_per_capita => db.world_gdp_per_capita,
            ),
        )
        datadict!(db, newdata)
        # br = db.brazil.total
        # db.world.per_country.datadict[:brazil] = br
        @info "COVID-19 database created." summary(db)

    else
        @info "COVID-19 database imported from cache."
    end

    if compute_model
        @info "Computing SEIR model..."
        SEIRModel!(db; kwargs...)
        @info "SEIR model for COVID-19 database computed."
    end

    if do_export
        for trial ∈ 1:export_trials
            try
                @info "Exporting..."
                paths = export_data(db; pretty = pretty, kwargs...)
                break
                @info "COVID-19 database exported." paths
            catch exception
                @warn "Failed exporting data." exception trial = "$trial of $export_trials"
                sleep(1)
            end
        end
    end
    db
end


text:
using DiffPrecTest
using Statistics, StatsBase, LinearAlgebra
using SparseArrays
using ProximalBase, CoordinateDescent, CovSel
using Random, Distributions
using JLD

pArr = [100, 200, 500, 300]
elemArr = [(5,5), (8, 7), (50, 25)]
n = 100
est      = Array{Any}(undef, 5)   # number of methods

# rep   = parse(Int,ARGS[1])
# ip    = parse(Int,ARGS[2])
# iElem = parse(Int,ARGS[3])
# dir = ARGS[4]

rep = 1
ip = 2
iElem = 1

Random.seed!(134)

p = pArr[ip]
# generate model
Ω = Matrix{Float64}(I, p, p)
for l=1:p-1
    Ω[l  , l+1] = 0.6
    Ω[l+1, l  ] = 0.6
end
for l=1:p-2
    Ω[l  , l+2] = 0.3
    Ω[l+2, l  ] = 0.3
end
d = Vector{Float64}(undef, p)
rand!(Uniform(0.5, 2.5), d)
d .= sqrt.(d)
D = Diagonal(d)
Σ = inv(Symmetric(D * Ω * D))

dist_X = MvNormal(convert(Matrix, Σ))
dist_Y = MvNormal(convert(Matrix, Σ))

# generate data
Random.seed!(134 + rep)
X = rand(dist_X, n)'
Y = rand(dist_Y, n)'

@show (ip, iElem)

ri, ci = elemArr[iElem]
indE = (ci - 1) * p + ri
indOracle = [indE]

@time est[1], _, _, indS = DiffPrecTest.estimate(SymmetricNormal(), X, Y, indE)
@show indS
@time est[2] = DiffPrecTest.estimate(SeparateNormal(), X, Y, indE)
@time est[3] = DiffPrecTest.estimate(SymmetricOracleBoot(), X, Y, indS)
@time est[4] = DiffPrecTest.estimate(SymmetricOracleNormal(), Symmetric(cov(X)), n, Symmetric(cov(Y)), n, indOracle)
@time est[5] = DiffPrecTest.estimate(SymmetricOracleBoot(), X, Y, indOracle)

# @save "$(dir)/res_$(ip)_$(iElem)_$(rep).jld" est

ext: jl


text:
export orthonormal, findlocalminima
#####################################################################################
#                                Minima Maxima                                      #
#####################################################################################
export findlocalminima

"""
    findlocalminima(s)
Return the indices of the local minima the timeseries `s`. If none exist,
return the index of the minimum (as a vector).
Starts of plateaus are also considered local minima.
"""
function findlocalminima(s::Vector{T})::Vector{Int} where {T}
    minimas = Int[]
    N = length(s)
    flag = false
    first_point = 0
    for i = 2:N-1
        if s[i-1] > s[i] && s[i+1] > s[i]
            flag = false
            push!(minimas, i)
        end
        # handling constant values
        if flag
            if s[i+1] > s[first_point]
                flag = false
                push!(minimas, first_point)
            elseif s[i+1] < s[first_point]
                flag = false
            end
        end
        if s[i-1] > s[i] && s[i+1] == s[i]
            flag = true
            first_point = i
        end
    end
    # make sure there is no empty vector returned
    if isempty(minimas)
        _, mini = findmin(s)
        return [mini]
    else
        return minimas
    end
end

"""
    findlocalextrema(y) -> max_ind, min_ind
Find the local extrema of given array `y`, by scanning point-by-point. Return the
indices of the maxima (`max_ind`) and the indices of the minima (`min_ind`).
"""
function findlocalextrema(y)
    @inbounds begin
        l = length(y)
        i = 1
        maxargs = Int[]
        minargs = Int[]
        if y[1] > y[2]
            push!(maxargs, 1)
        elseif y[1] < y[2]
            push!(minargs, 1)
        end

        for i in 2:l-1
            left = i-1
            right = i+1
            if  y[left] < y[i] > y[right]
                push!(maxargs, i)
            elseif y[left] > y[i] < y[right]
                push!(minargs, i)
            end
        end

        if y[l] > y[l-1]
            push!(maxargs, l)
        elseif y[l] < y[l-1]
            push!(minargs, l)
        end
        return maxargs, minargs
    end
end

@deprecate localextrema findlocalextrema

#####################################################################################
#                                Conversions                                        #
#####################################################################################
to_matrix(a::AbstractVector{<:AbstractVector}) = cat(2, a...)
to_matrix(a::AbstractMatrix) = a
function to_Smatrix(m)
    M = to_matrix(m)
    a, b = size(M)
    return SMatrix{a, b}(M)
end
to_vectorSvector(a::AbstractVector{<:SVector}) = a
function to_vectorSvector(a::AbstractMatrix)
    S = eltype(a)
    D, k = size(a)
    ws = Vector{SVector{D, S}}(k)
    for i in 1:k
        ws[i] = SVector{D, S}(a[:, i])
    end
    return ws
end

"""
    orthonormal([T,] D, k) -> ws
Return a matrix `ws` with `k` columns, each being
an `D`-dimensional orthonormal vector.
`T` is the return type and can be either `SMatrix` or `Matrix`

returns either SMatrix{D, k} if D*k < 100, otherwise Dxk Matrix{Float64}
"""
function orthonormal end

orthonormal(D, k) = D*k < 100 ? orthonormal(SMatrix, D, k) : orthonormal(Matrix, D, k)

@inline function orthonormal(T::Type, D::Int, k::Int)
    k > D && throw(ArgumentError("k must be ≤ D"))
    if T == SMatrix
        q = qr(rand(SMatrix{D, k})).Q
    elseif T == Matrix
        q = Matrix(qr(rand(Float64, D, k)).Q)
    end
    q
end

"""
    hcat_lagged_values(Y, s::Vector, τ::Int) -> Z
Add the `τ` lagged values of the timeseries `s` as additional component to `Y`
(`Vector` or `Dataset`), in order to form a higher embedded
dataset `Z`. The dimensionality of `Z` is thus equal to that of `Y` + 1.
"""
function hcat_lagged_values(Y::AbstractDataset{D,T}, s::Vector{T}, τ::Int) where {D, T<:Real}
    N = length(Y)
    MM = length(s)
    @assert N ≤ MM

    MMM = MM - τ
    M = min(N, MMM)
    data = Vector{SVector{D+1, T}}(undef, M)
    @inbounds for i in 1:M
        data[i] = SVector{D+1, T}(Y[i]..., s[i+τ])
    end
    return Dataset{D+1, T}(data)
end

function hcat_lagged_values(Y::Vector{T}, s::Vector{T}, τ::Int) where {T<:Real}
    N = length(Y)
    MM = length(s)
    @assert N ≤ MM
    MMM = MM - τ
    M = min(N, MMM)
    return Dataset(view(Y, 1:M), view(s, τ+1:τ+M))
end


text:
using Plots, Distributions, Statistics, LaTeXStrings, ProgressBars, StatsBase, JLD

function SSMC(f :: Function , R :: Real, r :: Real, Θ :: Real, θ ::Real , N :: Integer)
    ΔR= R - r
    Δθ= Θ - θ
    dist=[ ΔR .* rand(N) .+ r, Δθ .* rand(N) .+ θ ]
    f₍ₓ₎=[]
    for n in 1:N
        push!(f₍ₓ₎, f(getindex.(dist, n)...))
    end
    ∫f₍ₓ₎dx= mean(f₍ₓ₎)* ΔR * Δθ
    error= std(f₍ₓ₎)/sqrt(N)
    return ∫f₍ₓ₎dx, error
end

samples=10^4
ErrorList=[]
IntegralRes=[]
for i in ProgressBar(1:samples)
    I, E=SSMC((r, θ) -> r^3 * sin(θ) * cos(θ) * (3 + r * cos(θ)), 1, 0 , π, 0, 10^5)
    push!(IntegralRes, I)
    push!(ErrorList, E)
end
histogram(IntegralRes, dpi=600, ylabel="Probability Density", xlabel="I",
    title="The Distribution of I's Results")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter7\\7.2\\Figs\\IHistogram.png")
histogram(ErrorList, dpi=600,xlabel="Probability Density", ylabel="Error",
    title="The Distribution of Errors")
savefig("C:\\Users\\Narges\\Documents\\GitHub\\computational_physics\\chapter7\\7.2\\Figs\\EHistogram.png")
save("../../chapter7/7.2/data.jld","Idata",IntegralRes, "Edata", ErrorList )


text:
@testset "maximumlikelihoodeva.jl" begin

    n = 1000

    x₁ = Variable("x₁", randn(n)/3)
    x₂ = Variable("x₂", randn(n))
    x₃ = Variable("x₂", randn(n))

    θ = [1; 1; 1; -.5; 1; .1]

    μ = θ[1] .+ x₂.value*θ[2] .+ x₃.value*θ[3]
    ϕ = θ[4] .+ x₁.value*θ[5]
    ξ = θ[6]

    pd = GeneralizedExtremeValue.(μ, exp.(ϕ), ξ)

    y = rand.(pd)

    model = BlockMaxima(Variable("y", y), locationcov=[x₂; x₃], logscalecov = [x₁])

    fm = Extremes.MaximumLikelihoodEVA(model, θ)

    @testset "getdistribution(fittedmodel)" begin
        @test all(Extremes.getdistribution(fm) .== pd)
    end

    @testset "loglike(fd)" begin
        # Test with known values
        @test Extremes.loglike(fm) ≈ sum(logpdf.(pd,y))
    end

    @testset "quantile(fm, p)" begin
        # p not in [0, 1] throws
        @test_throws AssertionError Extremes.quantile(fm, -1)

        # Test with known values
        @test quantile(fm, .99) ≈ quantile.(pd,.99)
    end

    @testset "returnlevel(fm, returnPeriod)" begin
        # returnPeriod < 0 throws
        @test_throws AssertionError Extremes.returnlevel(fm, -1)

        # Test with known values
        @test quantile(fm, .99) ≈ quantile.(pd,.99)
    end

    @testset "hessian(fm)" begin
        # Tested in cint(fm)
    end

    @testset "parametervar(fm)" begin
        # Tested in cint(fm)
    end

    @testset "cint(fm)" begin

        @test_throws AssertionError cint(fm, 1)

        # Test with known values
        # Evaluating the confidence interval empirical level using a Monte Carlo approach.
        # Use of the hessian(fm) and parametervar function.

        niter = 1000

        cover = falses(niter,6)

        Threads.@threads  for i in 1:niter
            y = rand.(pd)
            m = BlockMaxima(Variable("y",y), locationcov = [x₂, x₃], logscalecov=[x₁])
            fm = gevfit(m, θ)
            ci = cint(fm)
            for j in eachindex(θ)
                cover[i,j] = ci[j][1] < θ[j] < ci[j][2]
            end
        end

        @test all(count.(eachcol(cover))/niter .> .9)

    end

    @testset "quantilevar(rl, level)" begin
        # Tested in cint(rl)
    end

    @testset "cint(rl, confidencelevel)" begin
        # confidencelevel not in [0, 1] throws
        @test_throws AssertionError cint(fm, 1)

        # Test with known values

        θ = [100; log(5); .1]

        pd = GeneralizedExtremeValue(θ[1],exp(θ[2]),θ[3])

        q = quantile(pd, 1-1/100)

        niter = 1000

        cover = falses(niter)

        Threads.@threads  for i in 1:niter
            y = rand(pd, 300)
            m = BlockMaxima(Variable("y",y))
            fm = gevfit(m, θ)
            rl = returnlevel(fm, 100)
            ci = cint(rl)
            cover[i] = ci[][1] < q < ci[][2]
        end

        @test count(cover)/niter .> .9

    end

    n = 1000
    nobservation = 10000
    nobsperblock = 1

    x = Variable("x",randn(n))

    θ = [-.5; 1; 0.1]

    ϕ = θ[1] .+ θ[2]*x.value
    σ = exp.(ϕ)
    ξ = θ[3]

    pd = GeneralizedPareto.(σ, ξ)

    y = rand.(pd)

    model = ThresholdExceedance(Variable("y", y), logscalecov = [x])

    fm = Extremes.MaximumLikelihoodEVA(model, θ)

    @testset "getdistribution(fittedmodel)" begin
        @test all(Extremes.getdistribution(fm) .== pd)
    end

    @testset "loglike(fd)" begin
        # Test with known values
        @test Extremes.loglike(fm) ≈ sum(logpdf.(pd,y))
    end

    @testset "quantile(fm, p)" begin
        # p not in [0, 1] throws
        @test_throws AssertionError Extremes.quantile(fm, -1)

        # Test with known values
        @test quantile(fm, .99) ≈ quantile.(pd,.99)
    end

    @testset "returnlevel(fm, returnPeriod)" begin
        # returnPeriod < 0 throws
        @test_throws AssertionError returnlevel(fm, 0, nobservation, nobsperblock, -100)

        # Test with known values
        rl = returnlevel(fm, 0, nobservation, nobsperblock, 100)
        p = 1-nobservation/(100 * nobsperblock * n)
        @test rl.value ≈ quantile.(pd, p)
    end

    @testset "hessian(fm)" begin
        # Tested in cint(fm)
    end

    @testset "parametervar(fm)" begin
        # Tested in cint(fm)
    end

    @testset "cint(fm)" begin

        @test_throws AssertionError cint(fm, 1)

        # Test with known values
        # Evaluating the confidence interval empirical level using a Monte Carlo approach.
        # Use of the hessian(fm) and parametervar function.

        niter = 1000

        cover = falses(niter,length(θ))

        Threads.@threads  for i in 1:niter
            y = rand.(pd)
            m = ThresholdExceedance(Variable("y",y), logscalecov=[x])
            fm = gpfit(m, θ)
            ci = cint(fm)
            for j in eachindex(θ)
                cover[i,j] = ci[j][1] < θ[j] < ci[j][2]
            end
        end

        @test all(count.(eachcol(cover))/niter .> .9)

    end

    @testset "quantilevar(rl, level)" begin
        # Tested in cint(rl)
    end

    @testset "cint(rl, confidencelevel)" begin
        # confidencelevel not in [0, 1] throws
        @test_throws AssertionError cint(fm, 1)

        # Test with known values

        θ = [log(5); .1]

        pd = GeneralizedPareto(exp(θ[1]),θ[2])

        n = 100
        nobservation = 1000
        nobsperblock = 1

        p = 1-nobservation/(100 * nobsperblock * n)

        q = quantile(pd, p)

        niter = 1000

        cover = falses(niter)

        Threads.@threads  for i in 1:niter
            y = rand(pd, n)
            m = ThresholdExceedance(Variable("y",y))
            fm = gpfit(m, θ)
            rl = returnlevel(fm, 0, nobservation, nobsperblock, 100)
            ci = cint(rl)
            cover[i] = ci[][1] < q < ci[][2]
        end

        @test count(cover)/niter .> .9

    end

    @testset "showfittedEVA(io, obj, prefix)" begin
        # print does not throw
        buffer = IOBuffer()
        @test_logs Extremes.showfittedEVA(buffer, fm, prefix = "\t")
    end

end


text:
"""
  boore_thompson_2014(m, r_ps::U, src::SourceParameters{S,T}) where {S<:Float64, T<:Real, U<:Real}

Boore & Thompson (2014) excitation duration model.
"""
function boore_thompson_2014(m, r_ps::U, src::SourceParameters{S,T}) where {S<:Float64, T<:Real, U<:Real}
  # source duration
  fa, fb, ε = corner_frequency(m, src)
  if src.model == :Atkinson_Silva_2000
    Ds = 0.5 / fa + 0.5 / fb
  else
    Ds = 1.0 / fa
  end
  # path duration
  if r_ps == 0.0
    return Ds * oneunit(U)
  elseif r_ps > 0.0 && r_ps <= 7.0
    return Ds + r_ps / 7.0 * 2.4
  elseif r_ps > 7.0 && r_ps <= 45.0
    return Ds + 2.4 + (r_ps - 7.0)/(45.0 - 7.0)*(8.4 - 2.4)
  elseif r_ps > 45.0 && r_ps <= 125.0
    return Ds + 8.4 + (r_ps - 45.0)/(125.0 - 45.0)*(10.9 - 8.4)
  elseif r_ps > 125.0 && r_ps <= 175.0
    return Ds + 10.9 + (r_ps - 125.0)/(175.0 - 125.0)*(17.4 - 10.9)
  elseif r_ps > 175.0 && r_ps <= 270.0
    return Ds + 17.4 + (r_ps - 175.0)/(270.0 - 175.0)*(34.2 - 17.4)
  elseif r_ps > 270.0
    return Ds + 34.2 + 0.156 * (r_ps - 270.0)
  else
    if T <: Dual
      return T(NaN)
    elseif U <: Dual
      return U(NaN)
    else
      return S(NaN)
    end
  end
end

boore_thompson_2014(m, r_ps, fas::FourierParameters) = boore_thompson_2014(m, r_ps, fas.source)

"""
  excitation_duration(m, r_ps::U, src::SourceParameters{S,T}, rvt::RandomVibrationParameters) where {S<:Float64,T<:Real,U<:Real}

Generic function implementing excitation duration models.

Currently, only the Boore & Thompson (2014) model is implemented.
"""
function excitation_duration(m, r_ps::U, src::SourceParameters{S,T}, rvt::RandomVibrationParameters) where {S<:Float64,T<:Real,U<:Real}
  if rvt.dur_ex == :BT14
    return boore_thompson_2014(m, r_ps, src)
  else
    if T <: Dual
      return T(NaN)
    elseif U <: Dual
      return U(NaN)
    else
      return S(NaN)
    end
  end
end

excitation_duration(m, r_ps, fas::FourierParameters, rvt::RandomVibrationParameters) = excitation_duration(m, r_ps, fas.source, rvt)

# Definition of Boore & Thompson (2012) constant values to be used within subsequent functions
include("coefficients/PJSbooreThompson2012.jl")

"""
  boore_thompson_2012_coefs(idx_m::T, idx_r::T; region::Symbol=:WNA) where T<:Int

Base function to extract the coefficients of the Boore & Thompson (2012) rms duration model.
"""
function boore_thompson_2012_coefs(idx_m::T, idx_r::T; region::Symbol=:WNA) where T<:Int
  idx = (idx_r - 1) * num_m_ii_bt12 + idx_m
  if region == :ENA
    @inbounds c = coefs_ena_bt12[idx,3:9]
    return c
  else
    @inbounds c = coefs_wna_bt12[idx,3:9]
    return c
  end
end

"""
  boore_thompson_2012_base(η::S, c::Vector{T}, ζ::T=0.05) where {S<:Real,T<:Float64}

Base function to compute the Boore & Thompson (2012) rms duration model for known magnitude and distance.
"""
function boore_thompson_2012_base(η::S, c::Vector{T}, ζ::T=0.05) where {S<:Real,T<:Float64}
  @inbounds ηc3 = η^c[3]
  @inbounds ratio = (c[1] + c[2]*((1.0 - ηc3)/(1.0 + ηc3))) * (1.0 + c[4]/(2π*ζ)*( η /(1.0 + c[5]*η^c[6]) )^c[7])
  return ratio
end

"""
  boore_thompson_2012(m, r_ps::T, src::SourceParameters, sdof::Oscillator, rvt::RandomVibrationParameters) where {S<:Float64,T<:Real}

Boore & Thompson (2012) rms duration model. Also outputs the excitation duration (given that its required within the rms duration calculations).

"""
function boore_thompson_2012(m, r_ps::T, src::SourceParameters, sdof::Oscillator, rvt::RandomVibrationParameters) where {S<:Float64,T<:Real}
  # for magnitude and distance that don't match coefficient tables we need to use bilinear interpolation
  # get the excitation duration (as recommended by Boore & Thompson, 2012)
  Dex = boore_thompson_2014(m, r_ps, src)
  # get the oscillator period
  T_n = period(sdof)
  ζ = sdof.ζ_n
  # define the η parameter as T_n/Dex
  η = T_n / Dex

  # impose limits on the magnitudes and distances
  m = ( m < 4.0 ) ? 4.0 : m
  m = ( m > 8.0 ) ? 8.0 : m
  r_ps = ( r_ps < 2.0 ) ? 2.0 : r_ps
  r_ps = ( r_ps > 1262.0 ) ? 1262.0 : r_ps

  # get the bounding indicies
  i_lo = findlast( m_ii_bt12 .<= m )
  i_hi = findfirst( m_ii_bt12 .>= m )
  j_lo = findlast( r_jj_bt12 .<= r_ps )
  j_hi = findfirst( r_jj_bt12 .>= r_ps )

  # get the region for the rms duration model
  region = rvt.dur_region

  # check for situations in which the coefficients are known for the given m,r
  if i_lo == i_hi # we have coefficients for this magnitude
    if j_lo == j_hi # we have coefficients for this distance
      c = boore_thompson_2012_coefs(i_lo, j_lo; region=region)
      Dratio = boore_thompson_2012_base(η, c, ζ)
    else # we need to interpolate the distance values only
      r_lo = r_jj_bt12[j_lo]
      r_hi = r_jj_bt12[j_hi]
      c_lo = boore_thompson_2012_coefs(i_lo, j_lo; region=region)
      c_hi = boore_thompson_2012_coefs(i_lo, j_hi; region=region)
      D_lo = boore_thompson_2012_base(η, c_lo, ζ)
      D_hi = boore_thompson_2012_base(η, c_hi, ζ)
      lnDratio = log(D_lo) + (r_ps - r_lo)/(r_hi - r_lo)*log(D_hi/D_lo)
      Dratio = exp(lnDratio)
    end
  else # we need to interpolate the magnitudes
    if j_lo == j_hi # we have coefficients for this distance
      m_lo = m_ii_bt12[i_lo]
      m_hi = m_ii_bt12[i_hi]
      c_lo = boore_thompson_2012_coefs(i_lo, j_lo; region=region)
      c_hi = boore_thompson_2012_coefs(i_hi, j_lo; region=region)
      D_lo = boore_thompson_2012_base(η, c_lo, ζ)
      D_hi = boore_thompson_2012_base(η, c_hi, ζ)
      lnDratio = log(D_lo) + (m - m_lo)/(m_hi - m_lo)*log(D_hi/D_lo)
      Dratio = exp(lnDratio)
    else # we need to interpolate for this magnitude and distance
      # create the combinations of (m,r) from these indices
      m_lo = m_ii_bt12[i_lo]
      m_hi = m_ii_bt12[i_hi]
      r_lo = r_jj_bt12[j_lo]
      r_hi = r_jj_bt12[j_hi]

      # get the corresponding coefficients
      c_ll = boore_thompson_2012_coefs(i_lo, j_lo; region=region)
      c_lh = boore_thompson_2012_coefs(i_lo, j_hi; region=region)
      c_hl = boore_thompson_2012_coefs(i_hi, j_lo; region=region)
      c_hh = boore_thompson_2012_coefs(i_hi, j_hi; region=region)

      # get the corresponding ratios
      Dr_ll = boore_thompson_2012_base(η, c_ll, ζ)
      Dr_lh = boore_thompson_2012_base(η, c_lh, ζ)
      Dr_hl = boore_thompson_2012_base(η, c_hl, ζ)
      Dr_hh = boore_thompson_2012_base(η, c_hh, ζ)

      # bilinear interpolation
      D = [ log(Dr_ll) log(Dr_lh); log(Dr_hl) log(Dr_hh) ]
      knots = ([m_lo, m_hi], [r_lo, r_hi])
      # create interpolant
      itp_D = interpolate(knots, D, Gridded(Linear()))
      lnDratio = itp_D(m,r_ps)
      Dratio = exp(lnDratio)
    end
  end
  Drms = Dratio * Dex
  return (Drms, Dex, Dratio)
end

boore_thompson_2012(m, r_ps, fas::FourierParameters, sdof::Oscillator, rvt::RandomVibrationParameters) = boore_thompson_2012(m, r_ps, fas.source, sdof, rvt)

# Definition of Boore & Thompson (2015) constant values to be used within subsequent functions
include("coefficients/PJSbooreThompson2015.jl")

"""
  boore_thompson_2015_coefs(idx_m::T, idx_r::T; region::Symbol=:WNA) where T<:Int

Base function to extract the coefficients of the Boore & Thompson (2015) rms duration model.
"""
function boore_thompson_2015_coefs(idx_m::T, idx_r::T; region::Symbol=:WNA) where T<:Int
  idx = (idx_r - 1) * num_m_ii_bt15 + idx_m
  if region == :ENA
    @inbounds c = coefs_ena_bt15[idx,3:9]
    return c
  else
    @inbounds c = coefs_wna_bt15[idx,3:9]
    return c
  end
end

"""
  boore_thompson_2015_base(η::S, c::Vector{T}, ζ::T=0.05) where {S<:Real,T<:Float64}

Base function to compute the Boore & Thompson (2015) rms duration model for known magnitude and distance.
"""
function boore_thompson_2015_base(η::S, c::Vector{T}, ζ::T=0.05) where {S<:Real,T<:Float64}
  @inbounds ηc3 = η^c[3]
  @inbounds ratio = (c[1] + c[2]*((1.0 - ηc3)/(1.0 + ηc3))) * (1.0 + c[4]/(2π*ζ)*( η /(1.0 + c[5]*η^c[6]) )^c[7])
  return ratio
end

"""
  boore_thompson_2015(m, r_ps::T, src::SourceParameters, sdof::Oscillator, rvt::RandomVibrationParameters) where {S<:Float64,T<:Real}

Boore & Thompson (2015) rms duration model. Also outputs the excitation duration (given that its required within the rms duration calculations).

"""
function boore_thompson_2015(m, r_ps::T, src::SourceParameters, sdof::Oscillator, rvt::RandomVibrationParameters) where {S<:Float64,T<:Real}
  # for magnitude and distance that don't match coefficient tables we need to use bilinear interpolation of the log Drms values
  # get the excitation duration (as recommended by Boore & Thompson, 2015)
  Dex = boore_thompson_2014(m, r_ps, src)
  # get the oscillator period
  T_n = period(sdof)
  ζ = sdof.ζ_n
  # define the η parameter as T_n/Dex
  η = T_n / Dex

  # impose limits on the magnitudes and distances
  m = ( m < 2.0 ) ? 2.0 : m
  m = ( m > 8.0 ) ? 8.0 : m
  r_ps = ( r_ps < 2.0 ) ? 2.0 : r_ps
  r_ps = ( r_ps > 1262.0 ) ? 1262.0 : r_ps

  # get the bounding indicies
  i_lo = findlast( m_ii_bt15 .<= m )
  i_hi = findfirst( m_ii_bt15 .>= m )
  j_lo = findlast( r_jj_bt15 .<= r_ps )
  j_hi = findfirst( r_jj_bt15 .>= r_ps )

  # get the region for the rms duration model
  region = rvt.dur_region

  # check for situations in which the coefficients are known for the given m,r
  if i_lo == i_hi # we have coefficients for this magnitude
    if j_lo == j_hi # we have coefficients for this distance
      c = boore_thompson_2015_coefs(i_lo, j_lo; region=region)
      Dratio = boore_thompson_2015_base(η, c, ζ)
    else # we need to interpolate the distance values only
      r_lo = r_jj_bt15[j_lo]
      r_hi = r_jj_bt15[j_hi]
      c_lo = boore_thompson_2015_coefs(i_lo, j_lo; region=region)
      c_hi = boore_thompson_2015_coefs(i_lo, j_hi; region=region)
      D_lo = boore_thompson_2015_base(η, c_lo, ζ)
      D_hi = boore_thompson_2015_base(η, c_hi, ζ)
      lnDratio = log(D_lo) + (r_ps - r_lo)/(r_hi - r_lo)*log(D_hi/D_lo)
      Dratio = exp(lnDratio)
    end
  else # we need to interpolate the magnitudes
    if j_lo == j_hi # we have coefficients for this distance
      m_lo = m_ii_bt15[i_lo]
      m_hi = m_ii_bt15[i_hi]
      c_lo = boore_thompson_2015_coefs(i_lo, j_lo; region=region)
      c_hi = boore_thompson_2015_coefs(i_hi, j_lo; region=region)
      D_lo = boore_thompson_2015_base(η, c_lo, ζ)
      D_hi = boore_thompson_2015_base(η, c_hi, ζ)
      lnDratio = log(D_lo) + (m - m_lo)/(m_hi - m_lo)*log(D_hi/D_lo)
      Dratio = exp(lnDratio)
    else # we need to interpolate for this magnitude and distance
      # create the combinations of (m,r) from these indices
      m_lo = m_ii_bt15[i_lo]
      m_hi = m_ii_bt15[i_hi]
      r_lo = r_jj_bt15[j_lo]
      r_hi = r_jj_bt15[j_hi]

      # get the corresponding coefficients
      c_ll = boore_thompson_2015_coefs(i_lo, j_lo; region=region)
      c_lh = boore_thompson_2015_coefs(i_lo, j_hi; region=region)
      c_hl = boore_thompson_2015_coefs(i_hi, j_lo; region=region)
      c_hh = boore_thompson_2015_coefs(i_hi, j_hi; region=region)

      # get the corresponding ratios
      Dr_ll = boore_thompson_2015_base(η, c_ll, ζ)
      Dr_lh = boore_thompson_2015_base(η, c_lh, ζ)
      Dr_hl = boore_thompson_2015_base(η, c_hl, ζ)
      Dr_hh = boore_thompson_2015_base(η, c_hh, ζ)

      # bilinear interpolation
      D = [ log(Dr_ll) log(Dr_lh); log(Dr_hl) log(Dr_hh) ]
      knots = ([m_lo, m_hi], [r_lo, r_hi])
      # create interpolant
      itp_D = interpolate(knots, D, Gridded(Linear()))
      lnDratio = itp_D(m,r_ps)
      Dratio = exp(lnDratio)
    end
  end
  Drms = Dratio * Dex
  return (Drms, Dex, Dratio)
end

boore_thompson_2015(m, r_ps, fas::FourierParameters, sdof::Oscillator, rvt::RandomVibrationParameters) = boore_thompson_2015(m, r_ps, fas.source, sdof, rvt)

"""
    rms_duration(m::S, r_ps::T, src::SourceParameters, sdof::Oscillator, rvt::RandomVibrationParameters) where {S<:Float64,T<:Real}

Returns a 3-tuple of (Drms, Dex, Dratio), using a switch on `rvt.dur_rms`.
Default `:BT12` makes use of the `:BT14` model for excitation duration, `Dex`.
- `m` is magnitude
- `r_ps` is an equivalent point source distance
"""
function rms_duration(m, r_ps::T, src::SourceParameters, sdof::Oscillator, rvt::RandomVibrationParameters) where {S<:Float64,T<:Real}
  if rvt.dur_rms == :BT15
    return boore_thompson_2015(m, r_ps, src, sdof, rvt)
  elseif rvt.dur_rms == :BT12
    return boore_thompson_2012(m, r_ps, src, sdof, rvt)
  else
    return (T(NaN), T(NaN), T(NaN))
  end
end

rms_duration(m, r_ps, fas::FourierParameters, sdof::Oscillator, rvt::RandomVibrationParameters) = rms_duration(m, r_ps, fas.source, sdof, rvt)


text:
"""
ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля (присутствуют перегородки)
РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки поля промакированы
"""
function new_all_markers(r::Robot)
    steps=0
    arr_of_steps=[]
    while ((isborder(r,Sud))&&(isborder(r,West)))==false
        push!(arr_of_steps,moves!(r,West))
        push!(arr_of_steps,moves!(r,Sud))
        steps+=2
    end
    putmarker!(r)
    while isborder(r,Ost)==false
        move!(r,Ost)
        putmarker!(r)
    end
    side=West
    while isborder(r,Nord)==false
        oper=1
        move!(r,Nord)
        putmarker!(r)
        while oper!=0
            oper=find_decision(r,side) #Обходчик
            putmarker!(r)
        end
        side=inverse(side)
    end
    moves!(r,West)
    moves!(r,Sud)
    while (steps>0)==true
        side=isodd(steps) ? Ost : Nord
        for a ∈ 1:arr_of_steps[steps]
            move!(r,side)
        end
        steps-=1
    end
end

function find_decision(r::Robot, side::HorizonSide) #Обход локальных преград
    num_steps=0
    while isborder(r,side)==true&&isborder(r,next(side))==false
        move!(r,next(side))
        num_steps+=1
    end
    x=0 #Счётчик
    if isborder(r,side)==false
        move!(r,side)
        x+=1
    end
    if num_steps !=0
        while isborder(r,inverse(next(side)))==true
            move!(r,side)
            x+=1
        end
        for a ∈ 1:num_steps
            move!(r,inverse(next(side)))
        end
    end
    return x #Возвращение значения счётчика
end

function moves!(r::Robot,side::HorizonSide)
    num_steps=0
    while isborder(r,side)==false
        move!(r,side)
        num_steps+=1
    end
    return num_steps
end

function next(side::HorizonSide)
    (HorizonSide(mod(Int(side)+1,4)))
end
function inverse(side::HorizonSide)
    (HorizonSide(mod(Int(side)+2,4)))
end
#=
Описания массива и некоторых вспомогательных функций аналогичны тем,
которые представлены в предыдущих задачах, например, Задаче 5
=#

max_stars_repo_name: alexgruzberg/mirea-progs

max_issues_repo_name: alexgruzberg/mirea-progs

max_forks_repo_name: alexgruzberg/mirea-progs


text:
"""
    SegmentedLSE{V <: AbstractVector{<:AbstractFloat}} <: AbstractAggregation

[`AbstractAggregation`](@ref) implementing segmented log-sum-exp (LSE) aggregation:

``
f(\\{x_1, \\ldots, x_k\\}; r) = \\frac{1}{r}\\log \\left(\\frac{1}{k} \\sum_{i = 1}^{k} \\exp({r\\cdot x_i})\\right)
``

Stores a vector of parameters `ψ` that are filled into the resulting matrix in case an empty bag is encountered,
and a vector of parameters `r` used during computation.

See also: [`AbstractAggregation`](@ref), [`AggregationStack`](@ref), [`lse_aggregation`](@ref),
    [`SegmentedMax`](@ref), [`SegmentedMean`](@ref), [`SegmentedSum`](@ref), [`SegmentedPNorm`](@ref).
"""
struct SegmentedLSE{V <: AbstractVector{<:AbstractFloat}} <: AbstractAggregation
    ψ::V
    ρ::V
end

Flux.@functor SegmentedLSE

SegmentedLSE(T::Type, d::Int) = SegmentedLSE(zeros(T, d), randn(T, d))
SegmentedLSE(d::Int) = SegmentedLSE(Float32, d)

Flux.@forward SegmentedLSE.ψ Base.getindex, Base.length, Base.size, Base.firstindex, Base.lastindex,
        Base.first, Base.last, Base.iterate, Base.eltype

Base.vcat(as::SegmentedLSE...) = reduce(vcat, as |> collect)
function Base.reduce(::typeof(vcat), as::Vector{<:SegmentedLSE})
    SegmentedLSE(reduce(vcat, [a.ψ for a in as]),
                 reduce(vcat, [a.ρ for a in as]))
end

r_map(ρ) = @. softplus(ρ)
inv_r_map(r) = @. relu(r) + log1p(-exp(-abs(r)))

function (a::SegmentedLSE)(x::Maybe{AbstractMatrix{T}}, bags::AbstractBags,
                              w::Optional{AbstractVecOrMat{T}}=nothing) where T
    _check_agg(a, x)
    segmented_lse_forw(x, a.ψ, r_map(a.ρ), bags)
end

function _lse_precomp(x::AbstractMatrix, r, bags)
    T = promote_type(eltype(x), eltype(r))
    M = zeros(T, length(r), length(bags))
    @inbounds for (bi, b) in enumerate(bags)
        if !isempty(b)
            for i in eachindex(r)
                M[i, bi] = r[i] * x[i, first(b)]
            end
            for j in b[2:end]
                for i in eachindex(r)
                    M[i, bi] = max(M[i, bi], r[i] * x[i, j])
                end
            end
        end
    end
    M
end

function _segmented_lse_norm(x::AbstractMatrix, ψ, r, bags::AbstractBags, M)
    t = promote_type(eltype(x), eltype(ψ))
    y = zeros(t, length(r), length(bags))
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(ψ)
                y[i, bi] = ψ[i]
            end
        else
            for j in b
                for i in eachindex(r)
                    y[i, bi] += exp.(r[i] * x[i, j] - M[i, bi])
                end
            end
            for i in eachindex(r)
                y[i, bi] = (log(y[i, bi]) - log(length(b)) + M[i, bi]) / r[i]
            end
        end
    end
    y
end

segmented_lse_forw(::Missing, ψ::AbstractVector, r, bags::AbstractBags) = repeat(ψ, 1, length(bags))
function segmented_lse_forw(x::AbstractMatrix, ψ, r, bags::AbstractBags)
    M = _lse_precomp(x, r, bags)
    _segmented_lse_norm(x, ψ, r, bags, M)
end

function segmented_lse_back(Δ, y, x, ψ, r, bags, M)
    dx = zero(x)
    dψ = zero(ψ)
    dr = zero(r)
    s1 = zeros(eltype(x), length(r))
    s2 = zeros(eltype(x), length(r))

    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(ψ)
                dψ[i] += Δ[i, bi]
            end
        else
            for i in eachindex(r)
                s1[i] = s2[i] = zero(eltype(x))
            end
            for j in b
                for i in eachindex(r)
                    e = exp(r[i] * x[i, j] - M[i, bi])
                    s1[i] += e
                    s2[i] += x[i, j] * e
                end
            end
            for j in b
                for i in eachindex(r)
                    dx[i, j] = Δ[i, bi] * exp(r[i] * x[i, j] - M[i, bi]) / s1[i]
                end
            end
            for i in eachindex(r)
                dr[i] += Δ[i, bi] * (s2[i]/s1[i] - y[i, bi]) / r[i]
            end
        end
    end
    dx, dψ, dr, NoTangent(), ZeroTangent()
end

function segmented_lse_back(Δ, ::Missing, ψ, bags)
    dψ = zero(ψ)
    @inbounds for (bi, b) in enumerate(bags)
        for i in eachindex(ψ)
            dψ[i] += Δ[i, bi]
        end
    end
    ZeroTangent(), dψ, ZeroTangent(), NoTangent(), ZeroTangent()
end

function ChainRulesCore.rrule(::typeof(segmented_lse_forw),
        x::AbstractMatrix, ψ::AbstractVector, r::AbstractVector, bags::AbstractBags)
    M = _lse_precomp(x, r, bags)
    y = _segmented_lse_norm(x, ψ, r, bags, M)
    grad = Δ -> (NoTangent(), segmented_lse_back(Δ, y, x, ψ, r, bags, M)...)
    y, grad
end

function segmented_lse_forw(::typeof(segmented_lse_forw), x::Missing, ψ::AbstractVector, r::AbstractVector, bags::AbstractBags)
    y = segmented_lse_forw(x, ψ, r, bags)
    grad = Δ -> (NoTangent(), segmented_lse_back(Δ, x, ψ, bags)...)
    y, grad
end


text:
using JuMP
using Mayday
using Base.Test

function test_simple_sos_problem()
	# Let's find a polynomial which is > 1 for x < -1 and < -1 for x > 1
	# v(x) + d1(x) * (x + 1) - 1 is SOS
	# -v(x) + d2(x) * (1 - x) - 1 is SOS
	# d1(x) is SOS
	# d2(x) is SOS

	model = Model()
	degree = 4
	v = defPolynomial(model, [:x], degree)
	@show v
	d1 = defSoSPolynomial(model, [:x], degree / 2)
	d2 = defSoSPolynomial(model, [:x], degree / 2)
	x = generator(:x)
	addSoSConstraint(model, v + d1 * (x + 1) - 1)
	addSoSConstraint(model, -v + d2 * (1 - x) - 1)
	solve(model)

	# Extract the result:
	v = getvalue(v)
	@show v

	# Verify it:
	x = linspace(-2, 2)
	y = map(x -> evaluate(v, x), x)

	for i = 1:length(x)
		if x[i] < -1
			@test y[i] >= 1
		elseif x[i] > 1
			@test y[i] <= -1
		end
	end
end

test_simple_sos_problem()


text:
module MoreEntities

using ..Ahorn, Maple

placements = Dict{String, Ahorn.EntityPlacement}(
    "Change Inventory Trigger" => Ahorn.EntityPlacement(
        (x, y) -> Maple.Trigger("everest/changeInventoryTrigger", x=x, y=y, width=16, height=16),
        "rectangle",
        Dict{String, Any}(
            "inventory" => ""
        )
    )
)

end

max_stars_repo_name: thefoxbot/MoreEntities

max_issues_repo_name: thefoxbot/MoreEntities

max_forks_repo_name: thefoxbot/MoreEntities


text:
using Documenter, LITS
const PSYPATH = dirname(pathof(LITS))

makedocs(
    modules = [LITS],
    format = Documenter.HTML(mathengine = Documenter.MathJax()),
    sitename = "LITS.jl",
    pages = Any[ # Compat: `Any` for 0.4 compat
        "Home" => "index.md",
        # "User Guide" => "man/guide.md",
        "Tutorials" => Any[
            "Tutorial 0: Data Creation" => "Examples/example_data.md",
            "Tutorial 1: OMIB" => "Examples/example_OMIB.md",
            "Tutorial 2: Dynamic Lines" => "Examples/example_lines.md",
        ],
        "Models" => Any[
            "Network" => "Models/network.md",
            "Reference Frames" => "Models/srf.md",
            "Generator" => "Models/gens.md",
            "Inverter" => "Models/inverters.md",
            "Small Signal" => "Models/small.md",
        ],
    ],
)

deploydocs(
    repo = "github.com/Energy-MAC/LITS.jl",
    branch = "gh-pages",
    target = "build",
    #deps = Deps.pip("pygments", "mkdocs", "python-markdown-math"),
    make = nothing,
)


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLSVRProjStd

using Test, RLinearAlgebra, LinearAlgebra, Random

@testset "LSVR Projection Standard -- Procedural" begin

    # Supertype and aliases
    @test supertype(LinSysVecRowProjStd) == LinSysVecRowProjection
    @test Kaczmarz == LinSysVecRowProjStd
    @test ART == LinSysVecRowProjStd

    # Verify that residual projection is zero
    Random.seed!(1010)

    A = rand(10,5)
    x = rand(5)
    b = A * x

    rsub = Kaczmarz()

    for i = 1:10
        @test let

            z = rand(5)
            RLinearAlgebra.rsubsolve!(rsub, z, (A[i,:], b[i]), i)

            dot(z, A[i,:]) ≈ b[i]
        end
    end

end

end # End module


text:
export i_gemm!, i_gemv!, i_sum, i_unsafe_addto, i_muladd
export maxloc, gaussian_log

const LogLikeNumber{T} = Union{ULogarithmic{T}, Tropical{T}}

const TropicalG{T,TG} = Tropical{GVar{T,TG}}

import TropicalNumbers: content
@fieldview content(x::ULogarithmic) = x.log
NiLang.chfield(x::Tropical{T}, ::typeof(content), val::T) where T = Tropical{T}(val)

# AD wrappers
for T in [:Tropical]
    @eval NiLang.AD.GVar(x::$T) = $T(GVar(content(x), zero(content(x))))
    @eval (_::Type{Inv{$T}})(x::$T) = content(x)
    @eval NiLang.AD.grad(x::$T{<:GVar}) = $T(grad(content(x)))
    @eval (_::Type{Inv{GVar}})(x::$T{<:GVar}) = $T((~GVar)(content(x)))

    @eval Base.one(x::$T{GVar{T,GT}}) where {T, GT} = one($T{GVar{T,GT}})
    @eval Base.one(::Type{$T{GVar{T,GT}}}) where {T,GT} = $T(GVar(zero(T), zero(GT)))
    @eval Base.zero(x::$T{GVar{T,GT}}) where {T,GT} =zero($T{GVar{T,GT}})
    @eval Base.zero(::Type{$T{GVar{T,T}}}) where T = $T(GVar(zero(T), zero(T)))
end

function NiLang.loaddata(::Type{Array{<:LogLikeNumber{GVar{T,T}}}}, data::Array{<:LogLikeNumber{T}}) where {T}
    GVar.(data)
end
import NiLang.NiLangCore: deanc

function deanc(x::T, v::T) where T<:LogLikeNumber
    x === v || deanc(content(x), content(v))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 9778402c-ff65-11ea-2af4-0fb4c0d1e78e
using Pkg, DrWatson

# ╔═╡ 977878f8-ff65-11ea-356c-8b889c8037b4
begin
	using Distributions
	using StatsPlots
	using StatsBase
	using LaTeXStrings
	using CSV
	using DataFrames
	using LinearAlgebra
	using Random
	using LinearAlgebra
	using StatisticalRethinking
end

# ╔═╡ 6426d13e-ff65-11ea-0e72-732d3a4de0f1
md"## clip-09-02s.jl"

# ╔═╡ 9778f07e-ff65-11ea-3d7f-35819188f757
md"### Snippet 9.2"

# ╔═╡ 97850170-ff65-11ea-16e9-af62429b511f
md"##### Number of samples."

# ╔═╡ 9785681a-ff65-11ea-3216-2b23598b1bef
T = 1000

# ╔═╡ 9793e066-ff65-11ea-2067-01c3caab76bf
md"##### Compute radial distance."

# ╔═╡ 97946ee6-ff65-11ea-0c08-43c9a52edf8d
rad_dist(x) = sqrt(sum(x .^ 2))

# ╔═╡ 979c7abe-ff65-11ea-2001-336b1bb02638
md"##### Plot densities."

# ╔═╡ 97ad992a-ff65-11ea-1e52-43039c4f5462
begin
	fig = density(xlabel="Radial distance from mode", ylabel="Density")
	for d in [1, 10, 100, 1000]
		m = MvNormal(zeros(d), Diagonal(ones(d)))
		local y = rand(m, T)
		rd = [rad_dist( y[:, i] ) for i in 1:T]
		density!(rd, lab="d=$d")
	end
	fig
end

# ╔═╡ 97b95942-ff65-11ea-1da3-23a8c109072b
md"## End of clip-09-02s.jl"

# ╔═╡ Cell order:
# ╟─6426d13e-ff65-11ea-0e72-732d3a4de0f1
# ╠═9778402c-ff65-11ea-2af4-0fb4c0d1e78e
# ╠═977878f8-ff65-11ea-356c-8b889c8037b4
# ╟─9778f07e-ff65-11ea-3d7f-35819188f757
# ╟─97850170-ff65-11ea-16e9-af62429b511f
# ╠═9785681a-ff65-11ea-3216-2b23598b1bef
# ╟─9793e066-ff65-11ea-2067-01c3caab76bf
# ╠═97946ee6-ff65-11ea-0c08-43c9a52edf8d
# ╟─979c7abe-ff65-11ea-2001-336b1bb02638
# ╠═97ad992a-ff65-11ea-1e52-43039c4f5462
# ╟─97b95942-ff65-11ea-1da3-23a8c109072b


text:
module TestDCA
using Test, ArDCA
function energy(x, J, h)
    q, q, N, N = size(J)
    eJ = 0.0
    @inbounds @simd for i in 1:N - 1
        for j in i + 1:N
            eJ -= J[x[i],x[j],i,j]
        end
    end
    eh = 0.0
    @inbounds @simd for i in 1:N
        eh -= h[x[i],i]
    end
    return eJ + eh
end

function complete_dataset(N, q)
    Z = zeros(Int, N, q^N)
    bounds = ntuple(x -> q, N)
    ctr = 0
    @inbounds for i in CartesianIndices(bounds)
        ctr += 1
        for j in 1:N
            Z[j,ctr] = i[j]
        end
    end
    Z
end

function generateWZJh(N, q)
    Jasym = randn(q, q, N, N)
    h = randn(q, N)
    J = 0.5 * (permutedims(Jasym, [2,1,4,3]) + Jasym)
    for i in 1:N
        for a in 1:q
            for b in 1:q
                J[a,b,i,i] = 0.0
            end
        end
    end
    Z = complete_dataset(N, q)
    W = [exp(-energy(Z[:,i], J, h)) for i in 1:q^N]
    W .= W ./ sum(W)
    return W, Z, J, h
end

function testDCA(N,q;
                 verbose::Bool=false,
                 epsconv::Real=1e-50,
                 lambdaJ::Real=0.0,
                 lambdaH::Real=0.0,
                 maxit::Integer=10000,
                 method::Symbol=:LD_LBFGS,
                 permorder::Union{Symbol,Vector{Int}}=:NATURAL,
                 dBthreshold::Real=-40)

    W, Z, J, h = generateWZJh(N, q)
    arnet, arvar = ardca(Z,W,
                    lambdaJ=lambdaJ,
                    lambdaH=lambdaH,
                    epsconv=epsconv,
                    verbose=verbose,
                    maxit=maxit,
                    method=method,
                    permorder=permorder)

    dBval = 10 * log10(sum(abs2, (sum(log.(arnet(arvar)), dims=1) |> vec .|> exp) - W) / length(W))
    println("testing N=$N\tq=$q\tpermorder=$permorder\tΔ=$dBval [dB]")
    @test dBval < dBthreshold
end
for q in 2:4
    for N in 2:5
        for order in [:NATURAL,:RANDOM, :ENTROPIC, :REV_ENTROPIC]
            testDCA(N, q, permorder=order,epsconv=1e-20)
        end
        println()
    end
end
printstyled("All TestDCA passed!\n",color=:light_green,bold=true)
end # end

max_issues_repo_name: pagnani/ArDCA

max_forks_repo_name: pagnani/ArDCA


text:
function writeFilterSPICE(FileName,
    RDrive,
    NumDriveElements,
    LDrive,
    CDrive,
    SerCap,
    CParAct,
    LFilt,
    LFilt2_ESR,
    LFilt1_ESR,
    LFilt2,
    CFilt,
    CFilt2,
    LNotch,
    LNotch_ESR,
    CNotch,
    LNotch2,
    LNotch2_ESR,
    CNotch2,
    RDamp,
    LDrivePair,
    CDrivePair)

    open(FileName,"w") do file
        write(file,
"Version 4
SHEET 1 1972 680
WIRE 352 -80 256 -80
WIRE 512 -80 432 -80
WIRE 1184 -64 1072 -64
WIRE 1376 -64 1264 -64
WIRE -32 64 -128 64
WIRE 112 64 48 64
WIRE 224 64 176 64
WIRE 256 64 256 -80
WIRE 256 64 224 64
WIRE 288 64 256 64
WIRE 432 64 368 64
WIRE 512 64 512 -80
WIRE 512 64 496 64
WIRE 768 64 512 64
WIRE 912 64 768 64
WIRE 1072 64 1072 -64
WIRE 1072 64 912 64
WIRE 1104 64 1072 64
WIRE 1264 64 1184 64
WIRE 1376 64 1376 -64
WIRE 1376 64 1328 64
WIRE 1424 64 1376 64
WIRE 1600 64 1424 64
WIRE 1712 64 1664 64
WIRE 1872 64 1776 64
WIRE 1872 96 1872 64
WIRE -128 128 -128 64
WIRE 224 160 224 64
WIRE 224 160 144 160
WIRE 272 160 224 160
WIRE 768 160 768 64
WIRE 912 160 912 64
WIRE 144 192 144 160
WIRE 272 192 272 160
WIRE 1424 208 1424 64
WIRE -128 240 -128 208
WIRE 1872 240 1872 176
WIRE 768 272 768 240
WIRE 912 272 912 240
WIRE -128 384 -128 320
WIRE 144 384 144 272
WIRE 144 384 -128 384
WIRE 272 384 272 256
WIRE 272 384 144 384
WIRE 768 384 768 336
WIRE 768 384 272 384
WIRE 912 384 912 336
WIRE 912 384 768 384
WIRE 1424 384 1424 272
WIRE 1424 384 912 384
WIRE 1872 384 1872 320
WIRE 1872 384 1424 384
WIRE -128 432 -128 384
FLAG -128 432 0
SYMBOL voltage -128 224 R0
WINDOW 123 0 0 Left 0
WINDOW 39 0 0 Left 0
SYMATTR InstName V1
SYMATTR Value AC 1
SYMBOL ind 64 48 R90
WINDOW 0 5 56 VBottom 2
WINDOW 3 32 56 VTop 2
SYMATTR InstName LFilt1
SYMATTR Value $(num2SciString(LFilt))
SYMATTR SpiceLine Rser=$(num2SciString(LFilt1_ESR))
SYMBOL cap 176 48 R90
WINDOW 0 0 32 VBottom 2
WINDOW 3 32 32 VTop 2
SYMATTR InstName CFilt1
SYMATTR Value $(num2SciString(CFilt))
SYMBOL res -144 112 R0
SYMATTR InstName RSrc
SYMATTR Value 1m
SYMBOL ind 384 48 R90
WINDOW 0 5 56 VBottom 2
WINDOW 3 32 56 VTop 2
SYMATTR InstName LFilt2
SYMATTR Value $(num2SciString(LFilt))
SYMATTR SpiceLine Rser=$(num2SciString(LFilt1_ESR))
SYMBOL cap 496 48 R90
WINDOW 0 0 32 VBottom 2
WINDOW 3 32 32 VTop 2
SYMATTR InstName CFilt2
SYMATTR Value $(num2SciString(CFilt))
SYMBOL ind 128 176 R0
SYMATTR InstName LFilt3
SYMATTR Value $(num2SciString(LFilt2))
SYMATTR SpiceLine Rser=$(num2SciString(LFilt2_ESR))
SYMBOL cap 256 192 R0
SYMATTR InstName CFilt3
SYMATTR Value $(num2SciString(CFilt2))
SYMBOL ind 752 144 R0
SYMATTR InstName LNotch
SYMATTR Value $(num2SciString(LNotch))
SYMATTR SpiceLine Rser=$(num2SciString(LNotch_ESR))
SYMBOL cap 752 272 R0
SYMATTR InstName CNotch
SYMATTR Value $(num2SciString(CNotch))
SYMBOL cap 1408 208 R0
SYMATTR InstName CPar
SYMATTR Value $(num2SciString(CParAct))
SYMBOL cap 1664 48 R90
WINDOW 0 0 32 VBottom 2
WINDOW 3 32 32 VTop 2
SYMATTR InstName CSer
SYMATTR Value $(num2SciString(SerCap))
SYMBOL cap 1776 48 R90
WINDOW 0 0 32 VBottom 2
WINDOW 3 32 32 VTop 2
SYMATTR InstName CDrive
SYMATTR Value $(num2SciString(CDrive/NumDriveElements))
SYMBOL ind 1856 80 R0
SYMATTR InstName LDrive
SYMATTR Value $(num2SciString(LDrive*NumDriveElements))
SYMBOL res 1856 224 R0
SYMATTR InstName RDrive
SYMATTR Value $(num2SciString(RDrive*NumDriveElements))
SYMBOL res 448 -96 R90
WINDOW 0 0 56 VBottom 2
WINDOW 3 32 56 VTop 2
SYMATTR InstName RDamp1
SYMATTR Value $(num2SciString(RDamp))
SYMBOL ind 1200 48 R90
WINDOW 0 5 56 VBottom 2
WINDOW 3 32 56 VTop 2
SYMATTR InstName LDrivePair
SYMATTR Value $(num2SciString(LDrivePair))
SYMATTR SpiceLine Rser=$(num2SciString(LFilt1_ESR))
SYMBOL cap 1328 48 R90
WINDOW 0 0 32 VBottom 2
WINDOW 3 32 32 VTop 2
SYMATTR InstName CDrivePair
SYMATTR Value $(num2SciString(CDrivePair))
SYMBOL res 1280 -80 R90
WINDOW 0 0 56 VBottom 2
WINDOW 3 32 56 VTop 2
SYMATTR InstName RDamp2
SYMATTR Value $(num2SciString(RDamp))
SYMBOL ind 896 144 R0
SYMATTR InstName LNotch1
SYMATTR Value $(num2SciString(LNotch2))
SYMATTR SpiceLine Rser=$(num2SciString(LNotch2_ESR))
SYMBOL cap 896 272 R0
SYMATTR InstName CNotch1
SYMATTR Value $(num2SciString(CNotch2))
TEXT 264 432 Left 2 !.ac oct 1k 1k 500k")
    end
end

function num2SciString(Val::Real,sigdigits=4)
    Val = round(Val;sigdigits=sigdigits)

    if abs(Val)>1e3
        OutVal =  ((Val *1e-3))
         ReturnStr =  "$(OutVal)k"
    elseif abs(Val)>1
        OutVal =  ((Val ))
         ReturnStr =  "$OutVal"
    elseif abs(Val)>1e-3
        OutVal =  ((Val *1e3))
         ReturnStr =  "$(OutVal)m"
    elseif abs(Val)>1e-6
        OutVal =  ((Val *1e6))
         ReturnStr =  "$(OutVal)u"
    elseif abs(Val)>1e-9
        OutVal =  ((Val *1e9))
         ReturnStr =  "$(OutVal)n"
    elseif abs(Val)>1e-12
        OutVal =  ((Val *1e12))
         ReturnStr =  "$(OutVal)p"
    else
        OutVal =  ((Val *1e15))
         ReturnStr =  "$(OutVal)f"
    end

    if length(ReturnStr)>(sigdigits+2)
        ReturnStr = ReturnStr[1:sigdigits+1]*ReturnStr[end]
    end
    return ReturnStr
end


text:
using Test
using StaticArrays, IntervalSets
import ClimaCore.DataLayouts: IJFH
import ClimaCore:
    Fields, Domains, Meshes, Topologies, Spaces, Operators, Geometry
using StaticArrays, IntervalSets, LinearAlgebra

using OrdinaryDiffEq

include("sphere_sphericalharmonics.jl")

@testset "Scalar biharmonic equation in 2D sphere" begin
    # ∇⁴ u(θ,φ) = f(θ,φ)
    # Please note that (θ,φ) is colatitude and azimuth as in the spherical coordinate system, that
    # correspond to π/2-latitude and longitude on Earth
    # True solution (eigenfunction): u(θ,φ) = Yₗᵐ(θ,φ)
    # => ∇⁴ u(θ,φ) = f(θ,φ) = l^2*(l+1)^2/(radius)^4 * u(θ,φ)

    FT = Float64

    radius = FT(1)
    domain = Domains.SphereDomain(radius)

    Ne = 16
    mesh = Meshes.EquiangularCubedSphere(domain, Ne)
    grid_topology = Topologies.Topology2D(mesh)

    Nq = 6
    quad = Spaces.Quadratures.GLL{Nq}()
    space = Spaces.SpectralElementSpace2D(grid_topology, quad)
    coords = Fields.coordinate_field(space)

    l = 7
    m = 4

    # Define eigensolution
    u = @. Ylm(l, m, coords.lat, coords.long)

    function ∇⁴(u)
        grad = Operators.Gradient()
        wdiv = Operators.WeakDivergence()
        diff = @. wdiv(grad(u))
        Spaces.weighted_dss!(diff)

        hyperdiff = @. wdiv(grad(diff))
        Spaces.weighted_dss!(hyperdiff)
        return hyperdiff
    end

    # Call the diffusion operator
    hyperdiff = ∇⁴(u)

    # compute the exact solution
    exact_solution = @. l^2 * (l + 1)^2 / (radius)^4 * u

    @test hyperdiff ≈ exact_solution rtol = 2e-2

end


text:
export DynamicDistributionModel, get_states, get_actions

import ReinforcementLearningEnvironments: observation_space, action_space

"""
    DynamicDistributionModel(f::Tf, ns::Int, na::Int) -> DynamicDistributionModel{Tf}

Use a general function `f` to store the transformations. `ns` and `na` are the number of states and actions.
"""
struct DynamicDistributionModel{Tf<:Function} <: AbstractDistributionBasedModel
    f::Tf
    ns::Int
    na::Int
end

observation_space(m::DynamicDistributionModel) = DiscreteSpace(m.ns)
action_space(m::DynamicDistributionModel) = DiscreteSpace(m.na)

(m::DynamicDistributionModel)(s, a) = m.f(s, a)


text:
module uFDTD

using Parameters

export uFDTDParameters

"""
    uFDTDParameters

    Configuration for the simulation.
"""
@with_kw struct uFDTDParameters
    # Spatial grid size.
    spatial_size::Int = 200  # [cells] square
    # Max time evolution.
    max_time_steps::Int = 250*4  # [steps]
end

"""
    simulate(sim_params::uFDTDParameters, hard_source=false, additive_source=false, directional_source=true)

    Run the simulation.
"""
function simulate(sim_params::uFDTDParameters, hard_source=false, additive_source=false, directional_source=true)
    @unpack spatial_size, max_time_steps = sim_params

    ε0 = 8.85418782e-12  # permittivity of free space [F/m] (F=Farad)
    # ε = εr * ε0
    # Permittivity is a constant of proportionality between electric displacement
    # and electric field intensity in a given medium.
    # force = charge * electric field
    # F     = Q      * E
    # ∇ × E = 0 : Curl of electric field (∇ is the del or nabla operator)
    # (∇f is a vector field of the df/dx df/dy etc...)
    # (∇ × E is a vector field representing rotational displacement)

    # Electric flux density
    # D = ε * E
    # E = D/ε  (= −∇V)
    # ∇ · D = ρv : Divergence of electric flux density
    # (ρv: electric charge density [C/m3])
    # (a scalar measure of strength of source or sink)
    # ∇²V = −ρ/e
    # Electrical conductivity: σ
    # Perfect electric conductors (PECs) have σ close to infinity
    # Current density J = σE (current in material, charges moving due to E field)
    # Charge Q moving at speed v in field B (magnetic flux density):
    # force = charge * speed cross-product magnetic flux density
    # F = Q * (v × B)

    # Magnetic Field
    # H = B / µ     (ignore local effect of material on flux):
    # B = μrμ0H = μH    (relative permeability and permeability of free space)
    µ0 = 4π * 1e-7  # permeability of free space [H/m] (H=Henry)
    # µ = µr * µ0
    η0 = sqrt(µ0 / ε0)  # Characteristic impedance of free space = ~ 377.0

    # Energy must no be able to propagate further than 1 spatial step: cΔt <= Δx
    c = 1 / sqrt(ε0 * µ0)  # speed of light in free space [m/s]
    Sc = 1  # Sc = c * Δt / Δx Courant number (we set it to 1, for now)
    # More generally: Δt must be <= 1/(c*sqrt(1/Δx^2 + 1/Δy^2 + 1/Δz^2))

    # Excitation function used in hard source and additive source
    peak_time = 30.0
    pulse_width = 10
    gaussian_pulse(t) = exp(-((t - peak_time) / pulse_width)^2)

    source_function = gaussian_pulse

    ez = zeros(Float64, spatial_size)  # z component of E field at a time step
    hy = zeros(Float64, spatial_size)  # y component of H field at a time step

    # setup material property (Elec): relative permittivity
    εr = zeros(Float64, spatial_size)
    εr_boundary_position = 100
    εr_material_value = 9
    for m = 1:spatial_size
        if m < εr_boundary_position
            εr[m] = 1.0
        else
            εr[m] = εr_material_value
        end
    end

    # setup material property (Mag): relative permeability
    µr = zeros(Float64, spatial_size)
    µr_material_value = 1
    for m = 1:spatial_size
        µr[m] = µr_material_value
    end

    probe0Dt = zeros(Float64, max_time_steps)
    probe1Dt = zeros(Float64, (spatial_size, max_time_steps))

    # Do time stepping, leap-frog method (update H then E)
    for qTime = 1:max_time_steps

        # Absorbing Boundary Condition on right side (works only when local Sc=1)
        hy[spatial_size] = hy[spatial_size-1]

        # Advance/update the Magnetic field
        for m = 1:spatial_size-1
            # Maxwell-Faraday, in 1D:  µ * dH/dt = dE/dx  (continuous version)
            #   µ * (H[q+1/2] - H[q-1/2]) / Δt = (E[m+1] - E[m]) / Δx  (discrete version)
            #   H[q+1/2] = H[q-1/2] + (E[m+1] - E[m]) * Δt / µΔx
            hy[m] = hy[m] + (ez[m + 1] - ez[m]) * Sc / η0 / µr[m]
        end

        # Directional source
        # Total-Field/Scattered-Field correction for Hy adjacent to TFSF boundary
        if directional_source
            hy[50] -= source_function(qTime) / η0
        end

        # Absorbing Boundary Condition on left side
        ez[1] = ez[2]

        # Advance/update the Electric field
        for m = 2:spatial_size
            # Maxwell-Ampère, in 1D:  ε * dE/dt = dH/dx  (continuous version)
            #   ε * (E[q+1] - E[q]) / Δt = (H[m+1/2] - H[m-1/2]) / Δx  (discrete version)
            #   E[q+1] = E[q] + (H[m+1/2] - H[m-1/2]) * Δt / εΔx
            ez[m] = ez[m] + (hy[m] - hy[m - 1]) * Sc * η0 / εr[m]
        end

        # Directional source
        # Total-Field/Scattered-Field correction for Ez adjacent to TFSF boundary
        if directional_source
            ez[51] += exp(-(qTime + 0.5 - (-0.5) - 30.) *
                           (qTime + 0.5 - (-0.5) - 30.) / 100.)
        end

        # Hardwire source on a node
        if hard_source
            ez[0] = source_function(qTime)
        end

        # Additive source on a node
        if additive_source
            ez[50] += source_function(qTime)
        end

        # add a point (0D) probe over time
        probe0Dt[qTime] = ez[50]

        # add a line (1D) probe over time
        probe1Dt[:, qTime] = ez
    end

    return probe0Dt, probe1Dt
end

end # module


text:
import PencilFFTs
using PencilArrays: Permutation

import Oceananigans.Solvers: poisson_eigenvalues, solve!
import Oceananigans.Architectures: architecture

struct DistributedFFTBasedPoissonSolver{P, F, L, λ, S, I}
    plan :: P
    global_grid :: F
    local_grid :: L
    eigenvalues :: λ
    storage :: S
    input_permutation :: I
end

architecture(solver::DistributedFFTBasedPoissonSolver) =
    architecture(solver.global_grid)

infer_transform(grid, d) = infer_transform(topology(grid, d)())
infer_transform(::Periodic) = PencilFFTs.Transforms.FFT!()
infer_transform(::Bounded) = PencilFFTs.Transforms.R2R!(FFTW.REDFT10)
infer_transform(::Flat) = PencilFFTs.Transforms.NoTransform!()

"""
    DistributedFFTBasedPoissonSolver(global_grid, local_grid)

Return a FFT-based solver for the Poisson equation,

```math
∇²x = b
```

for `MultiArch`itectures.

Supported configurations
========================

We support two "modes":

    1. Vertical pencil decompositions: two-dimensional decompositions in (x, y)
       for three dimensional problems that satisfy either `Nz > Rx` or `Nz > Ry`.

    2. One-dimensional decompositions in either x or y.

Above, `Nz = size(global_grid, 3)` and `Rx, Ry, Rz = architecture(local_grid).ranks`.

Other configurations that are decomposed in (x, y) but have too few `Nz`,
or any configuration decomposed in z, are not supported.

Algorithm for two-dimensional decompositions
============================================

For two-dimensional decompositions (useful for three-dimensional problems),
there are three forward transforms, three backward transforms,
and four transpositions requiring MPI communication. In the schematic below, the first
dimension is always the local dimension. In our implementation of the PencilFFTs algorithm,
we require _either_ `Nz >= Rx`, _or_ `Nz >= Ry`, where `Nz` is the number of vertical cells,
`Rx` is the number of ranks in x, and `Ry` is the number of ranks in `y`.
Below, we outline the algorithm for the case `Nz >= Rx`.
If `Nz < Rx`, but `Nz > Ry`, a similar algorithm applies with x and y swapped:

1. `first(storage)` is initialized with layout (z, x, y).
2. Transform along z.
3  Transpose + communicate to storage[2] in layout (x, z, y),
   which is distributed into `(Rx, Ry)` processes in (z, y).
4. Transform along x.
5. Transpose + communicate to last(storage) in layout (y, x, z),
   which is distributed into `(Rx, Ry)` processes in (x, z).
6. Transform in y.

At this point the three in-place forward transforms are complete, and we
solve the Poisson equation by updating `last(storage)`.
Then the process is reversed to obtain `first(storage)` in physical
space and with the layout (z, x, y).

Restrictions
============

The algorithm for two-dimensional decompositions requires that `Nz = size(global_grid, 3)` is larger
than either `Rx = ranks[1]` or `Ry = ranks[2]`, where `ranks` are configured when building `MultiArch`.
If `Nz` does not satisfy this condition, we can only support a one-dimensional decomposition.

Algorithm for one-dimensional decompositions
============================================

This algorithm requires a one-dimensional decomposition with _either_ `Rx = 1`
_or_ `Ry = 1`, and is important to support two-dimensional transforms.

For one-dimensional decompositions, we place the decomposed direction _last_.
If the number of ranks is `Rh = max(Rx, Ry)`, this algorithm requires that
_both_ `Nx > Rh` _and_ `Ny > Rh`. The resulting flow of transposes and transforms
is similar to the two-dimensional case. It remains somewhat of a mystery why this
succeeds (ie, why the last transform is correctly decomposed).
"""
function DistributedFFTBasedPoissonSolver(global_grid, local_grid)

    arch = architecture(local_grid)
    Rx, Ry, Rz = arch.ranks
    communicator = arch.communicator

    # We don't support distributing anything in z.
    Rz == 1 || throw(ArgumentError("Non-singleton ranks in the vertical are not supported by DistributedFFTBasedPoissonSolver."))

    gNx, gNy, gNz = size(global_grid)

    # Build _global_ eigenvalues
    topo = (TX, TY, TZ) = topology(global_grid)
    λx = poisson_eigenvalues(global_grid.Nx, global_grid.Lx, 1, TX())
    λy = poisson_eigenvalues(global_grid.Ny, global_grid.Ly, 2, TY())
    λz = poisson_eigenvalues(global_grid.Nz, global_grid.Lz, 3, TZ())

    # Drop singleton dimensions for compatibility with PencilFFTs' localgrid
    λx = dropdims(λx, dims=(2, 3))
    λy = dropdims(λy, dims=(1, 3))
    λz = dropdims(λz, dims=(1, 2))

    unpermuted_eigenvalues = (λx, λy, λz)

    # First we check if we can do a two-dimensional decomposition
    if gNz >= Rx
        input_permutation = Permutation(3, 1, 2)
        permuted_size = (gNz, gNx, gNy)
        processors_per_dimension = (Rx, Ry)
    elseif gNz >= Ry
        input_permutation = Permutation(3, 2, 1)
        permuted_size = (gNz, gNy, gNx)
        processors_per_dimension = (Ry, Rx)

    else # it has to be a one-dimensional decomposition

        Rx > 1 && Ry > 1 &&
            throw(ArgumentError("DistributedFFTBasedPoissonSolver requires either " *
                                "i) Nz > Rx, ii) Nz > Ry, iii) Rx = 1 _or_ iv) Ry = 1."))

        if Rx == 1 # x-local, y-distributed
            permuted_size = (gNz, gNx, gNy)
            input_permutation = Permutation(3, 1, 2)
            processors_per_dimension = (1, Ry)
        else # Ry == 1, y-local, x-distributed
            permuted_size = (gNz, gNy, gNx)
            input_permutation = Permutation(3, 2, 1)
            processors_per_dimension = (1, Rx)
        end
    end

    transforms = Tuple(infer_transform(global_grid, d) for d in Tuple(input_permutation))
    plan = PencilFFTs.PencilFFTPlan(permuted_size, transforms, processors_per_dimension, communicator)

    # Allocate memory for in-place FFT + transpositions
    storage = PencilFFTs.allocate_input(plan)

    # Permute the λ appropriately
    permuted_eigenvalues = Tuple(unpermuted_eigenvalues[d] for d in Tuple(input_permutation))
    eigenvalues = PencilFFTs.localgrid(last(storage), permuted_eigenvalues)

    return DistributedFFTBasedPoissonSolver(plan, global_grid, local_grid, eigenvalues, storage, input_permutation)
end

# solve! requires that `b` in `A x = b` (the right hand side)
# was computed and stored in first(solver.storage) prior to calling `solve!(x, solver)`.
# See: Models/NonhydrostaticModels/solve_for_pressure.jl
function solve!(x, solver::DistributedFFTBasedPoissonSolver)
    arch = architecture(solver.global_grid)
    multi_arch = architecture(solver.local_grid)

    # Apply forward transforms to b = first(solver.storage).
    solver.plan * solver.storage

    # Solve the discrete Poisson equation in wavenumber space
    # for x̂. We solve for x̂ in place, reusing b̂.
    x̂ = b̂ = last(solver.storage)
    λ = solver.eigenvalues
    @. x̂ = - b̂ / (λ[1] + λ[2] + λ[3])

    # Set the zeroth wavenumber and volume mean, which are undetermined
    # in the Poisson equation, to zero.
    if MPI.Comm_rank(multi_arch.communicator) == 0
        # This is an assumption: we *hope* PencilArrays allocates in this way
        parent(x̂)[1, 1, 1] = 0
    end

    # Apply backward transforms to x̂ = last(solver.storage).
    solver.plan \ solver.storage

    # xc is the backward transform of x̂.
    xc = first(solver.storage)

    # Copy the real component of xc to x.
    copy_event = launch!(arch, solver.local_grid, :xyz,
                         copy_permuted_real_component!, x, parent(xc), solver.input_permutation,
                         dependencies = device_event(arch))

    wait(device(arch), copy_event)

    return x
end

const ZXYPermutation = Permutation{(3, 1, 2), 3}
const ZYXPermutation = Permutation{(3, 2, 1), 3}

@kernel function copy_permuted_real_component!(ϕ, ϕc, ::ZXYPermutation)
    i, j, k = @index(Global, NTuple)
    @inbounds ϕ[i, j, k] = real(ϕc[k, i, j])
end

@kernel function copy_permuted_real_component!(ϕ, ϕc, ::ZYXPermutation)
    i, j, k = @index(Global, NTuple)
    @inbounds ϕ[i, j, k] = real(ϕc[k, j, i])
end


text:
driver_version() = (a = Cint[0]; rt.cudaDriverGetVersion(a); int(a[1]))
runtime_version() = (a = Cint[0]; rt.cudaRuntimeGetVersion(a); int(a[1]))


text:
"""
    module FEIO
All file input/output functionality. Includes 3DG interoperability
"""
module FEIO

using FELinearAlgebra

export dgtraits, read_array, write_array, read_solution, write_solution

"""
    dgtraits(a::Array{T,N}) where {T,N}
Compatibility of types with 3DG
"""
dgtraits(a::Array{T,N}) where {T,N} = dgtraits(eltype(a))
function dgtraits(T::Type)
    println(T)
    if T <: Bool
        return 2
    elseif T <: Integer
        return 0
    elseif T <: AbstractFloat # 3DG technically assumes C++ double
        return 1
    else
        error("Unknown DG Type")
    end
end
function dgtraits(val::Integer)
    if val == 0
        return Int64
    elseif val == 1
        return Float64
    elseif val == 2
        return Bool
    else
        error("Unknown DG Type")
    end
end

"""
    read_array(filename::AbstractString)
Read file created by `write_array()` and return the array represented by data
"""
function read_array(filename::AbstractString)
    file = open(filename, "r")
    array = read_array(file)
    close(file)
    return array
end
function read_array(io::IOStream)
    dim = read(io, Int64)
    size = Array{Int64,1}(undef, dim)
    read!(io, size)
    type = read(io, Int64)
    array = Array{dgtraits(type), dim}(undef, size...)
    read!(io, array)
    return array
end

"""
    write_array(filename::AbstractString, a)
Write the array to a specified filename in 3DG compatible format
"""
function write_array(filename::AbstractString, a)
    file = open(filename, "w")
    write_array(file, a)
    close(file)
end
function write_array(io::IOStream, a)
    if typeof(a) == SolutionVector
        a = a.data
    end
    sz = size(a)
    write(io, length(sz))
    write(io, sz...)
    write(io, dgtraits(a))
    write(io, a)
end

"""
    read_solution(filename, mesh)
Read file created by `write_solution()` and return SolutionVector defined
on a given mesh 3DG compatible format
"""
function read_solution(filename::AbstractString, mesh)
    # TODO: Will need to use mesh for parallel solutions
    data = read_array(filename)
    return SolutionVector(data)
end

"""
    write_solution(filename, u, mesh)
Write the solution defined on a given mesh to a specified filename in
3DG compatible format
"""
function write_solution(filename::AbstractString, u, mesh)
    # TODO: Will need to use mesh for parallel solutions
    write_array(filename, u.data)
end

end # module IO

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
function circular_plot(names::AbstractVector = [], values::AbstractVector = [];
            selected::Union{AbstractVector, Nothing} = nothing,
            radius::Union{AbstractVector, String, Nothing} = "80%",
            center::Union{AbstractVector, String, Nothing} = ["50%", "50%"],
            roseType::Union{String, Nothing} = nothing, # roseType = {"angle", "radius"}
            legend::Bool = false,
            kwargs...
            )

    # Set default to false if not passed by user
    selected == nothing ? selected = [false for x in 1:length(values)] : selected = selected

    ec = newplot(kwargs, ec_charttype = "circular")

    data_fmt = arrayofdicts(value = values, name = names, selected = selected)

    ec.series = [PieSeries(name = "Series 1", radius = radius, center = center, data = data_fmt, roseType = roseType)]

    #Add legend if requested
    legend ? legend!(ec) : nothing

    return ec

end


text:
# Install Julia packages
import Pkg
Pkg.add("IJulia")
Pkg.add("Random")
Pkg.add("JSON")
Pkg.add("HTTP")


text:
using Test

using OopsMorePermutations

@testset "Permutations.jl" begin
    include("type_basics.jl")
    include("group_properties.jl")
end


text:
using Documenter, DFControl

makedocs(modules = [DFControl],

         pages = Any["Home" => "index.md",
                     "Library" => Any["job.md",
                     "input.md",
                     "structure.md",
                     "atom.md",
                     "fileio.md",
                     "defaults.md",
                     "utils.md"]],
         sitename = "DFControl.jl",
         format = :html
        )

deploydocs(
    repo = "github.com/louisponet/DFControl.jl.git",
    target = "build",
    osname = "linux",
    julia  = "0.6",
    deps = nothing,
    make = nothing
    )


text:
# ---
# title: Log-Gabor filters v3
# id: demo_phasecong3
# cover: assets/phasecong3.png
# author: Peter Kovesi
# date: 2018-10-26
# ---

# Use of the function `phasecong3()` allows corner points to be detected as well. These
# corner points are a subset of the edge image and, unlike other corner detectors, their
# location is precise and stable over different scales.

using TestImages
using Images
using ImagePhaseCongruency

img = restrict(testimage("mandril_gray"))
(edges, corners) = phasecong3(img)

mosaic(
    img,
    adjust_histogram(Gray.(edges), LinearStretching()),
    adjust_histogram(corners, LinearStretching()),
    nrow=1
)
# Images from top to right: 1) original image 2) edges 3) corners

# save cover image #src
isdir("assets") || mkdir("assets") #src
save(joinpath("assets", "phasecong3.png"), adjust_histogram(Gray.(edges), LinearStretching())) #src

max_stars_repo_name: peterkovesi/ImagePhaseCongruency

max_issues_repo_name: peterkovesi/ImagePhaseCongruency

max_forks_repo_name: peterkovesi/ImagePhaseCongruency


text:
# Autogenerated wrapper script for mlpack_jll for i686-linux-gnu-cxx03
export libmlpack, libmlpack_julia_adaboost, libmlpack_julia_approx_kfn, libmlpack_julia_bayesian_linear_regression, libmlpack_julia_cf, libmlpack_julia_dbscan, libmlpack_julia_decision_stump, libmlpack_julia_decision_tree, libmlpack_julia_det, libmlpack_julia_emst, libmlpack_julia_fastmks, libmlpack_julia_gmm_generate, libmlpack_julia_gmm_probability, libmlpack_julia_gmm_train, libmlpack_julia_hmm_generate, libmlpack_julia_hmm_loglik, libmlpack_julia_hmm_train, libmlpack_julia_hmm_viterbi, libmlpack_julia_hoeffding_tree, libmlpack_julia_image_converter, libmlpack_julia_kde, libmlpack_julia_kernel_pca, libmlpack_julia_kfn, libmlpack_julia_kmeans, libmlpack_julia_knn, libmlpack_julia_krann, libmlpack_julia_lars, libmlpack_julia_linear_regression, libmlpack_julia_linear_svm, libmlpack_julia_lmnn, libmlpack_julia_local_coordinate_coding, libmlpack_julia_logistic_regression, libmlpack_julia_lsh, libmlpack_julia_mean_shift, libmlpack_julia_nbc, libmlpack_julia_nca, libmlpack_julia_nmf, libmlpack_julia_pca, libmlpack_julia_perceptron, libmlpack_julia_preprocess_binarize, libmlpack_julia_preprocess_describe, libmlpack_julia_preprocess_one_hot_encoding, libmlpack_julia_preprocess_scale, libmlpack_julia_preprocess_split, libmlpack_julia_radical, libmlpack_julia_random_forest, libmlpack_julia_softmax_regression, libmlpack_julia_sparse_coding, libmlpack_julia_util

using boost_jll
using armadillo_jll
using OpenBLAS_jll
JLLWrappers.@generate_wrapper_header("mlpack")
JLLWrappers.@declare_library_product(libmlpack, "libmlpack.so.3")
JLLWrappers.@declare_library_product(libmlpack_julia_adaboost, "libmlpack_julia_adaboost.so")
JLLWrappers.@declare_library_product(libmlpack_julia_approx_kfn, "libmlpack_julia_approx_kfn.so")
JLLWrappers.@declare_library_product(libmlpack_julia_bayesian_linear_regression, "libmlpack_julia_bayesian_linear_regression.so")
JLLWrappers.@declare_library_product(libmlpack_julia_cf, "libmlpack_julia_cf.so")
JLLWrappers.@declare_library_product(libmlpack_julia_dbscan, "libmlpack_julia_dbscan.so")
JLLWrappers.@declare_library_product(libmlpack_julia_decision_stump, "libmlpack_julia_decision_stump.so")
JLLWrappers.@declare_library_product(libmlpack_julia_decision_tree, "libmlpack_julia_decision_tree.so")
JLLWrappers.@declare_library_product(libmlpack_julia_det, "libmlpack_julia_det.so")
JLLWrappers.@declare_library_product(libmlpack_julia_emst, "libmlpack_julia_emst.so")
JLLWrappers.@declare_library_product(libmlpack_julia_fastmks, "libmlpack_julia_fastmks.so")
JLLWrappers.@declare_library_product(libmlpack_julia_gmm_generate, "libmlpack_julia_gmm_generate.so")
JLLWrappers.@declare_library_product(libmlpack_julia_gmm_probability, "libmlpack_julia_gmm_probability.so")
JLLWrappers.@declare_library_product(libmlpack_julia_gmm_train, "libmlpack_julia_gmm_train.so")
JLLWrappers.@declare_library_product(libmlpack_julia_hmm_generate, "libmlpack_julia_hmm_generate.so")
JLLWrappers.@declare_library_product(libmlpack_julia_hmm_loglik, "libmlpack_julia_hmm_loglik.so")
JLLWrappers.@declare_library_product(libmlpack_julia_hmm_train, "libmlpack_julia_hmm_train.so")
JLLWrappers.@declare_library_product(libmlpack_julia_hmm_viterbi, "libmlpack_julia_hmm_viterbi.so")
JLLWrappers.@declare_library_product(libmlpack_julia_hoeffding_tree, "libmlpack_julia_hoeffding_tree.so")
JLLWrappers.@declare_library_product(libmlpack_julia_image_converter, "libmlpack_julia_image_converter.so")
JLLWrappers.@declare_library_product(libmlpack_julia_kde, "libmlpack_julia_kde.so")
JLLWrappers.@declare_library_product(libmlpack_julia_kernel_pca, "libmlpack_julia_kernel_pca.so")
JLLWrappers.@declare_library_product(libmlpack_julia_kfn, "libmlpack_julia_kfn.so")
JLLWrappers.@declare_library_product(libmlpack_julia_kmeans, "libmlpack_julia_kmeans.so")
JLLWrappers.@declare_library_product(libmlpack_julia_knn, "libmlpack_julia_knn.so")
JLLWrappers.@declare_library_product(libmlpack_julia_krann, "libmlpack_julia_krann.so")
JLLWrappers.@declare_library_product(libmlpack_julia_lars, "libmlpack_julia_lars.so")
JLLWrappers.@declare_library_product(libmlpack_julia_linear_regression, "libmlpack_julia_linear_regression.so")
JLLWrappers.@declare_library_product(libmlpack_julia_linear_svm, "libmlpack_julia_linear_svm.so")
JLLWrappers.@declare_library_product(libmlpack_julia_lmnn, "libmlpack_julia_lmnn.so")
JLLWrappers.@declare_library_product(libmlpack_julia_local_coordinate_coding, "libmlpack_julia_local_coordinate_coding.so")
JLLWrappers.@declare_library_product(libmlpack_julia_logistic_regression, "libmlpack_julia_logistic_regression.so")
JLLWrappers.@declare_library_product(libmlpack_julia_lsh, "libmlpack_julia_lsh.so")
JLLWrappers.@declare_library_product(libmlpack_julia_mean_shift, "libmlpack_julia_mean_shift.so")
JLLWrappers.@declare_library_product(libmlpack_julia_nbc, "libmlpack_julia_nbc.so")
JLLWrappers.@declare_library_product(libmlpack_julia_nca, "libmlpack_julia_nca.so")
JLLWrappers.@declare_library_product(libmlpack_julia_nmf, "libmlpack_julia_nmf.so")
JLLWrappers.@declare_library_product(libmlpack_julia_pca, "libmlpack_julia_pca.so")
JLLWrappers.@declare_library_product(libmlpack_julia_perceptron, "libmlpack_julia_perceptron.so")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_binarize, "libmlpack_julia_preprocess_binarize.so")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_describe, "libmlpack_julia_preprocess_describe.so")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_one_hot_encoding, "libmlpack_julia_preprocess_one_hot_encoding.so")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_scale, "libmlpack_julia_preprocess_scale.so")
JLLWrappers.@declare_library_product(libmlpack_julia_preprocess_split, "libmlpack_julia_preprocess_split.so")
JLLWrappers.@declare_library_product(libmlpack_julia_radical, "libmlpack_julia_radical.so")
JLLWrappers.@declare_library_product(libmlpack_julia_random_forest, "libmlpack_julia_random_forest.so")
JLLWrappers.@declare_library_product(libmlpack_julia_softmax_regression, "libmlpack_julia_softmax_regression.so")
JLLWrappers.@declare_library_product(libmlpack_julia_sparse_coding, "libmlpack_julia_sparse_coding.so")
JLLWrappers.@declare_library_product(libmlpack_julia_util, "libmlpack_julia_util.so")
function __init__()
    JLLWrappers.@generate_init_header(boost_jll, armadillo_jll, OpenBLAS_jll)
    JLLWrappers.@init_library_product(
        libmlpack,
        "lib/libmlpack.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_adaboost,
        "lib/libmlpack_julia_adaboost.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_approx_kfn,
        "lib/libmlpack_julia_approx_kfn.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_bayesian_linear_regression,
        "lib/libmlpack_julia_bayesian_linear_regression.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_cf,
        "lib/libmlpack_julia_cf.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_dbscan,
        "lib/libmlpack_julia_dbscan.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_decision_stump,
        "lib/libmlpack_julia_decision_stump.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_decision_tree,
        "lib/libmlpack_julia_decision_tree.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_det,
        "lib/libmlpack_julia_det.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_emst,
        "lib/libmlpack_julia_emst.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_fastmks,
        "lib/libmlpack_julia_fastmks.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_gmm_generate,
        "lib/libmlpack_julia_gmm_generate.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_gmm_probability,
        "lib/libmlpack_julia_gmm_probability.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_gmm_train,
        "lib/libmlpack_julia_gmm_train.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_hmm_generate,
        "lib/libmlpack_julia_hmm_generate.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_hmm_loglik,
        "lib/libmlpack_julia_hmm_loglik.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_hmm_train,
        "lib/libmlpack_julia_hmm_train.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_hmm_viterbi,
        "lib/libmlpack_julia_hmm_viterbi.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_hoeffding_tree,
        "lib/libmlpack_julia_hoeffding_tree.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_image_converter,
        "lib/libmlpack_julia_image_converter.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_kde,
        "lib/libmlpack_julia_kde.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_kernel_pca,
        "lib/libmlpack_julia_kernel_pca.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_kfn,
        "lib/libmlpack_julia_kfn.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_kmeans,
        "lib/libmlpack_julia_kmeans.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_knn,
        "lib/libmlpack_julia_knn.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_krann,
        "lib/libmlpack_julia_krann.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_lars,
        "lib/libmlpack_julia_lars.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_linear_regression,
        "lib/libmlpack_julia_linear_regression.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_linear_svm,
        "lib/libmlpack_julia_linear_svm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_lmnn,
        "lib/libmlpack_julia_lmnn.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_local_coordinate_coding,
        "lib/libmlpack_julia_local_coordinate_coding.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_logistic_regression,
        "lib/libmlpack_julia_logistic_regression.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_lsh,
        "lib/libmlpack_julia_lsh.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_mean_shift,
        "lib/libmlpack_julia_mean_shift.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_nbc,
        "lib/libmlpack_julia_nbc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_nca,
        "lib/libmlpack_julia_nca.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_nmf,
        "lib/libmlpack_julia_nmf.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_pca,
        "lib/libmlpack_julia_pca.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_perceptron,
        "lib/libmlpack_julia_perceptron.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_preprocess_binarize,
        "lib/libmlpack_julia_preprocess_binarize.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_preprocess_describe,
        "lib/libmlpack_julia_preprocess_describe.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_preprocess_one_hot_encoding,
        "lib/libmlpack_julia_preprocess_one_hot_encoding.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_preprocess_scale,
        "lib/libmlpack_julia_preprocess_scale.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_preprocess_split,
        "lib/libmlpack_julia_preprocess_split.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_radical,
        "lib/libmlpack_julia_radical.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_random_forest,
        "lib/libmlpack_julia_random_forest.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_softmax_regression,
        "lib/libmlpack_julia_softmax_regression.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_sparse_coding,
        "lib/libmlpack_julia_sparse_coding.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libmlpack_julia_util,
        "lib/libmlpack_julia_util.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()