text:
################################################################################
#  Copyright 2022, Tom Van Acker (BASF)                                        #
################################################################################
# TFT.jl                                                                       #
# A Julia package for Taylor-Fourier Transform.                                #
# See http://github.com/timmyfaraday/TFT.jl                                    #
################################################################################

# abstract problem types
"""
    TFT.AbstractDTFTProblem

Base for types which define discrete taylor-fourier transform problems.
"""
abstract type AbstractDTFTProblem end

# problem types
"""
    TFT.DTFTProblem

Defines a discrete taylor-fourier transform problem.

Fields:
- `s::Vector{<:Number}`         | discrete signal [?]
- `t::StepRangeLen{<:Number}`   | discrete time [s]
- `h::Vector{<:Int}`            | harmonic numbers [-]
- `D::Int`                      | degree of the highest derivative [-]
- `F::Number`                   | fundamental frequency [Hz]
- `K::Int`                      | o-spline degree [-]
- `N::Int`                      | number of samples per fundamental cycle [-]
"""
struct DTFTProblem <: AbstractDTFTProblem
    """"discrete signal"""
    s::Vector{<:Number}
    """discrete time"""
    t::StepRangeLen{<:Number}
    """harmonic numbers"""
    h::Vector{<:Int}
    """degree of the highest derivative"""
    D::Int
    """fundamental frequency"""
    F::Number
    """o-spline degree"""
    K::Int
    """number of samples per fundamental cycle"""
    N::Int

    # base constructor
    DTFTProblem(s::Vector{<:Number}, t::StepRangeLen{<:Number}, h::Vector{<:Int},
                D::Int, F::Number, K::Int, N::Int) =
        new(s, t, h, D, F, K, N)
end

# problem builder
function build_problem(s, t, D, F, K)
    # checks
    length(s) == length(t)          || error("the length of the discrete signal and time are inconsistent")
    all(diff(t) .≈ (t[2] - t[1]))   || error("the discrete time is not equidistance")

    # reduce the discrete time to a range `rT`
    rT  = range(first(t), last(t), length=length(t))

    # determine the fundamental period [s]
    T   = 1 / F

    # find the number of samples in a fundamental cycle [-]
    N   = findfirst(x -> x ≈ first(t) + T, t) - 1

    # find the harmonics set based on the number of samples in a fundamental cycle
    h   = collect(0:N-2)

    return DTFTProblem(s, rT, h, D, F, K, N)
end
function build_problem(s, t, h, D, F, K)
    # checks
    length(s) == length(t)          || error("the length of the discrete signal and time are inconsistent")
    all(diff(t) .≈ (t[2] - t[1]))   || error("the discrete time is not equidistance")

    # reduce the discrete time to a range `rT`
    rT  = range(first(t), last(t), length=length(t))

    # determine the fundamental period [s]
    T   = 1 / F

    # find the number of samples in a fundamental cycle [-]
    N   = findfirst(x -> x ≈ first(t) + T, t) - 1

    return DTFTProblem(s, rT, h, D, F, K, N)
end

# solution types
"""
    TFT.AbstractDTFTSolution

Base for types which define discrete taylor-fourier transform solutions.
"""
abstract type AbstractDTFTSolution end

"""
    TFT.DTFTSolution

Defines a discrete taylor-fourier transform solution

Fields:
- `X::Dict{<:Int,Matrix{<:Number}}`     | dictionary of dynamic phasors ξₕ⁽ᴰ⁾(t)
- `prob::TFT.AbstractDTFTProblem`       | DTFT problem struct
"""
struct DTFTSolution <: AbstractDTFTSolution
    """dictionary of dynamic phasors ξₕ⁽ᴰ⁾"""
    X::Dict
    """DTFT problem struct"""
    prob::AbstractDTFTProblem
end


text:
export similar_basis_vector, function_slice

"""
    $SIGNATURES

Return a vector with `1` at index `i` and `0`s otherwise (ie a standard basis
vector), similar to the first argument.
"""
function similar_basis_vector(x::AbstractVector{T}, i::Integer) where T
    v = zeros(T, length(x))
    v[i] = one(T)
    v
end

similar_basis_vector(::SVector{N, T}, i::Integer) where {N, T} =
    SVector{N, T}(ntuple(j -> i == j ? one(T) : zero(T), N)...)

"""
    $SIGNATURES

Return a callable implementing ``λ ↦ f(x + λv)``.
"""
function function_slice(f, x::AbstractVector, v::AbstractVector)
    @argcheck length(x) == length(v)
    λ -> f(x .+ λ * v)
end

"""
    $SIGNATURES

Return a callable implementing ``λ ↦ f(x + λe)``, where ``e`` is the `i`th
standard basis vector.
"""
function_slice(f, x::AbstractVector, i::Integer) =
    function_slice(f, x, similar_basis_vector(x, i))

ext: jl


text:
@testset "Complete config test" begin
    mpSim, isOkay = initialiseFromDB( "db/simDB", showInfo = false )
    @test isOkay
end  # @testset "Complete config test"


text:
# tests that tests in regressionTest fail if regression reference data does not fit simulation data
using ModelicaScriptingTools
using Test
using CSV
using DataFrames

MoST = ModelicaScriptingTools

outdir = "test/out"
modeldir = "test/res"
refdir = "test/regRefData"

function removeColumn(csvfile, colname)
    data = DataFrame(CSV.File(csvfile))
    select!(data, Not(Symbol(colname)))
    CSV.write(csvfile, data)
end
function multiplyColumn(csvfile, colname, factor)
    data = DataFrame(CSV.File(csvfile))
    data[!, Symbol(colname)] = data[Symbol(colname)] .* factor
    CSV.write(csvfile, data)
end
function resetFiles()
    cp("$outdir/TwoVarExample_res.bak.csv", "$refdir/TwoVarExample_res.csv"; force=true)
    cp("$outdir/TwoVarExample_res.bak.csv", "$outdir/TwoVarExample_res.csv"; force=true)
end

MoST.withOMC(outdir, modeldir) do omc
    # setup simulation and reference data
    MoST.installAndLoad(omc, "Modelica"; version="3.2.3")
    MoST.loadModel(omc, "TwoVarExample")
    MoST.simulate(omc, "TwoVarExample")
    cp("$outdir/TwoVarExample_res.csv", "$outdir/TwoVarExample_res.bak.csv"; force=true)
    @testset "regressionTest" begin
        @testset "fails because of missing simulation variables" begin
            resetFiles()
            removeColumn("$outdir/TwoVarExample_res.csv", "i")
            MoST.regressionTest(omc, "TwoVarExample", refdir; relTol=1e-3, variableFilter="i|v")
        end
        @testset "fails because of missing reference variables" begin
            resetFiles()
            removeColumn("$refdir/TwoVarExample_res.csv", "i")
            MoST.regressionTest(omc, "TwoVarExample", refdir; relTol=1e-3, variableFilter="i|v")
        end
        @testset "fails because values are unequal" begin
            resetFiles()
            multiplyColumn("$outdir/TwoVarExample_res.csv", "i", 1.1)
            MoST.regressionTest(omc, "TwoVarExample", refdir; relTol=1e-3, variableFilter="i|v")
        end
    end
end

max_stars_repo_name: THM-MoTE/MoST

max_issues_repo_name: THM-MoTE/MoST

max_forks_repo_name: THM-MoTE/MoST


text:
function compute_dimensions(partition, blocks)
    dimensions = Vector{Int}()
    dims = 0
    next_idx = 1
    for block_idx in blocks
        while next_idx < block_idx
            # sum up dimensions of skipped blocks
            dims += length(partition[next_idx])
            next_idx += 1
        end
        append!(dimensions, partition[next_idx])
        next_idx += 1
    end
    return dimensions
end


text:
using Documenter, GeoTime

makedocs(;
    modules=[GeoTime],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/rafaqz/GeoTime.jl/blob/{commit}{path}#L{line}",
    sitename="GeoTime.jl",
    authors="Rafael Schouten <rafaelschouten@gmail.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/rafaqz/GeoTime.jl",
)


text:
export HomogeneousSpin

mutable struct HomogeneousSpin{D,S,C} <: AbstractHilbert
    n_sites::Int
    shape::Vector{Int}

    # Constrain on total spin number in units of S
    constraint::Int
end

"""
    HomogeneousSpins(N, S::Rational=1//2)

Constructs the Hilbert space of `N` identical spins-S (by default 1//2).
"""
function HomogeneousSpin(n_sites, S::Rational=1//2; total_sz::Union{Nothing,Rational,Int}=nothing)
    n_sites = n_sites isa AbstractGraph ? nv(n_sites) : n_sites

    @assert S.den == 2 || S.den == 1

    # N is dimension of local space
    if S.den == 2
        N = S.num +1
    elseif S.den == 1
        N = 2*S.num +1
    end

    constrained = isnothing(total_sz) ?  false : true
    if constrained
        try
            constraint = Int(total_sz * inv(S))
        catch err
            throw(ErrorException("total_sz is not in units of S"))
        end

        if constraint > n_sites
            throw(ErrorException("total_sz is too big"))
        end

        if N == 2
            (n_sites + constraint) % N == 0 || error("total_sz not valid")
        end
    else
        constraint = 0
    end

    return HomogeneousSpin{N,S,constrained}(n_sites, fill(N, n_sites),
                                          constraint)
end

Base.similar(hilb::HomogeneousSpin, N::Int) =
    HomogeneousSpin(N, spin(hilb))

@inline nsites(h::HomogeneousSpin) = h.n_sites
@inline local_dim(h::HomogeneousSpin{D}) where D = D
@inline local_dim(h::HomogeneousSpin{D}, i) where D = D
@inline shape(h::HomogeneousSpin) = h.shape
@inline spin(h::HomogeneousSpin{D,S}) where {D,S} = S

@inline spacedimension(h::HomogeneousSpin) = local_dim(h)^nsites(h)
@inline indexable(h::HomogeneousSpin) = spacedimension(h) != 0
@inline is_homogeneous(h::HomogeneousSpin) = true

@inline is_contrained(h::HomogeneousSpin{H,S,C}) where {H,S,C} = C
@inline constraint_limit(h::HomogeneousSpin) = h.constraint
@inline magnetization(h::HomogeneousSpin{H,S,false}) where {H,S} = false
@inline magnetization(h::HomogeneousSpin{H,S,true}) where {H,S} = S*h.constraint

state(arrT::AbstractArray, T::Type{<:Number}, h::HomogeneousSpin{N}, dims::Dims) where {N} =
    similar(arrT, T, nsites(h), dims...) .= -(N-1)

Base.show(io::IO, ::MIME"text/plain", h::HomogeneousSpin{N}) where N = begin
    print(io, "Hilbert Space with $(nsites(h)) identical spins $(N-1)/2 of dimension $(local_dim(h))")
    is_contrained(h) && print(io, " with constraint Sz=$(constraint_limit(h)*spin(h))")
end

Base.show(io::IO, h::HomogeneousSpin) =
    print(io, "HomogeneousSpin($(nsites(h)), $(local_dim(h)))")

## Operations

function flipat!(rng::AbstractRNG, σ::AState, h::HomogeneousSpin{N}, i) where N
    @boundscheck checkbounds(σ, i)
    T = eltype(σ)

    @inbounds old_val = σ[i]
    new_val = T(rand(rng))
    new_val = floor(new_val*(N-1))*2 - (N-1)
    @inbounds σ[i]    = new_val + T(2) * (new_val >= old_val)
    return old_val, new_val
end

# special case N== 2 to be faster
function flipat!(rng::AbstractRNG, σ::AState, h::HomogeneousSpin{2}, i)
    @boundscheck checkbounds(σ, i)
    T = eltype(σ)

    @inbounds old_val = σ[i]
    new_val = old_val == one(T) ? -one(T) : one(T)
    @inbounds σ[i]    = new_val

    return old_val, new_val
end

function setat!(σ::AState, h::HomogeneousSpin, i::Int, val)
    @boundscheck checkbounds(σ, i)

    @inbounds old_val = σ[i]
    @inbounds σ[i] = val

    return old_val
end

set_index!(σ::AState, h::HomogeneousSpin, i) = set!(σ, h, i)
function set!(σ::AState, h::HomogeneousSpin{N}, val::Integer) where N
    @boundscheck checkbounds_hilbert(h, val)

    val -= 1
    for i=1:nsites(h)
        val, tmp = divrem(val, N)
        @inbounds σ[i] = tmp * 2 - (N-1)
    end
    return σ
end

Base.@propagate_inbounds add!(σ::AState, h::HomogeneousSpin, val::Integer) =
    set!(σ, h, val+toint(σ, h))

function Random.rand!(rng::AbstractRNG, σ::AbstractArray, h::HomogeneousSpin{N,S,false}) where {N,S}
    T = eltype(σ)
    rand!(rng, σ)
    σ .= floor.(σ.*N).*2 .- (N-1)
    return σ
end

# Specialized for constrained fock spaces
function Random.rand!(rng::AbstractRNG, σ::AState, hilb::HomogeneousSpin{N,S,true}) where {N,S}
    T = eltype(σ)

    if N == 2
        m = constraint_limit(hilb) # magnetization in units of S
        nup   = (nsites(hilb) + m) ÷ 2
        ndown = (nsites(hilb) - m) ÷ 2

        @inbounds uview(σ, 1:nup) .= one(T)
        @inbounds uview(σ, nup+1:ndown+nup) .= -one(T)
        shuffle!(rng, σ)
    else
        throw(ErrorException("not implemented!"))
    end

    return σ
end

function toint(σ::AState, h::HomogeneousSpin{N}) where N
    tot = 0
    for (i,v)=enumerate(σ)
        tot += Int((v + (N-1))/2)*N^(i-1)
    end
    return tot + 1
end

function local_index(σ::AState, h::HomogeneousSpin{N}, site) where N
    @boundscheck checkbounds(σ, site)

    @inbounds id = Int((σ[site] + (N-1))/2)+1
    return id
end

function local_index(σ::AState, h::HomogeneousSpin{M}, sites::AbstractVector) where M
    @boundscheck checkbounds(σ, sites)

    idx = 1
    for (i,j)=enumerate(sites)
        @inbounds idx += Int((σ[j] + (M-1))/2) * M^(i-1)
    end
    return idx
end


text:
using ColorMaps
using Base.Test

# write your own tests here
@test 1 == 1

cm = ColorMaps.colormap("cubehelix")(64)

max_stars_repo_name: lobingera/ColorMaps


text:
module DestroyCrystalsTrigger

using ..Ahorn, Maple

placements = Dict{String, Ahorn.EntityPlacement}(
    "Destroy Crystals (In trigger)" => Ahorn.EntityPlacement(
        (x, y) -> Maple.Trigger("everest/goAwayCrystalsTrigger", x=x, y=y, width=16, height=16),
        "rectangle",
        Dict{String, Any}(
        "destroyEveryCrystal" => "InTrigger"
        )
    ),
    "Destroy Crystals (Every crystal)" => Ahorn.EntityPlacement(
        (x, y) -> Maple.Trigger("everest/goAwayCrystalsTrigger", x=x, y=y, width=16, height=16),
        "rectangle",
        Dict{String, Any}(
        "destroyEveryCrystal" => "EveryCrystal"
        )
    )
)

end

max_stars_repo_name: thefoxbot/MoreEntities

max_issues_repo_name: thefoxbot/MoreEntities

max_forks_repo_name: thefoxbot/MoreEntities


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Conversions between representations using Julia built-in system.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

function Base.convert(
    ::Type{KeplerianElements{Tepoch, T}},
    k::KeplerianElements
) where {Tepoch, T}
    return KeplerianElements{Tepoch, T}(k.t, k.a, k.e, k.i, k.Ω, k.ω, k.f)
end

function Base.convert(
    ::Type{OrbitStateVector{Tepoch, T}},
    sv::OrbitStateVector
) where {Tepoch, T}
    return OrbitStateVector{Tepoch, T}(sv.t, sv.r, sv.v, sv.a)
end

function Base.convert(
    ::Type{KeplerianElements},
    sv::OrbitStateVector{Tepoch, T}
) where {Tepoch, T}
    return Base.convert(KeplerianElements{Tepoch, T}, sv)
end

function Base.convert(
    KT::Type{KeplerianElements{Tepoch, T}},
    sv::OrbitStateVector
) where {Tepoch, T}
    k = sv_to_kepler(sv)
    return convert(KT, k)
end

function Base.convert(
    ::Type{OrbitStateVector},
    k::KeplerianElements{Tepoch, T}
) where {Tepoch, T}
    return Base.convert(OrbitStateVector{Tepoch, T}, k)
end

function Base.convert(
    ST::Type{OrbitStateVector{Tepoch, T}},
    k::KeplerianElements
) where {Tepoch, T}
    sv = kepler_to_sv(k)
    return convert(ST, sv)
end


text:
include("compat.jl")
using HTTP
using HTTP.IOExtras, HTTP.Sockets

@testset "WebSockets" begin

# for s in ["wss", "ws"]

#     HTTP.WebSockets.open("$s://echo.websocket.org") do io
#         println("writing")
#         write(io, "Foo")
#         println("testing")
#         @test !eof(io)
#         @test String(readavailable(io)) == "Foo"

#         println("writing again")
#         write(io, "Hello")
#         write(io, " There")
#         write(io, " World", "!")
#         println("closewrite")
#         closewrite(io)

#         println("reading response")
#         buf = IOBuffer()
#         write(buf, io)
#         @test String(take!(buf)) == "Hello There World!"
#     end

# end

p = 8085 # rand(8000:8999)
@async HTTP.listen(Sockets.localhost, p) do http
    if HTTP.WebSockets.is_upgrade(http.message)
        HTTP.WebSockets.upgrade(http) do ws
            while !eof(ws)
                data = readavailable(ws)
                write(ws, data)
            end
        end
    end
end

sleep(2)

println("Testing local server...")
HTTP.WebSockets.open("ws://127.0.0.1:$(p)") do ws
    write(ws, "Foo")
    @test String(readavailable(ws)) == "Foo"

    write(ws, "Bar")
    @test String(readavailable(ws)) == "Bar"
end

end # testset


text:
@testset "components" begin
    include("action_selectors/action_selectors.jl")
    include("approximators/approximators.jl")
    include("buffers/buffers.jl")
    include("environment_models/environment_models.jl")
end


text:
println("Compile libraries on main processor..."); flush(stdout)
include("imports.jl")
println("Finished loading libraries."); flush(stdout)

using Distributed
rmprocs(workers())
addprocs(32)
println("Compile libraries on workers ..."); flush(stdout)
@everywhere include("imports.jl")
println("Finished loading libraries."); flush(stdout)

Ks = collect(2:9)
snums = collect(1:4)
skewtmix = [true, false]

sims = dict_list(Dict(:K => Ks, :snum => snums, :skewtmix => skewtmix))

# Run simulations in parallel.
println("Fit models in parallel...")
@time res = pmap(run, sims, on_error=identity)

# Print results status.
foreach(z -> println("$(z[1]) => $(z[2])"), zip(sims, res))

# Post process
println("Post process in parallel...")
pp_res = pmap(postprocess, cdd.MCMC.ProgressBar(sims))
# postprocess(Dict(:beta=>1, :K=>5, :snum=>1))

# Print number of small clusters
sc_res = pmap(sim -> print_number_of_small_clusters(sim, p=0.005),
              cdd.MCMC.ProgressBar(sims), on_error=identity)

# DIC
println("Compute DIC sequentially...")
foreach(snum -> plot_dic(snum, Ks, skewtmix), snums)
foreach(snum -> plot_dic(snum, Ks, skewtmix, calibrate=true), snums)

# Send results to S3
Util.s3sync(from="$(Info.resultsdir_simstudy)/$(simname)",
            to="$(Info.awsbucket_simstudy)/$(simname)",
            tags=`--exclude '*.nfs'`)


text:
################################################################################
## Linear Regression
##   y ~ N(b0 + b1 * x, s2)
##   b0, b1 ~ N(0, 1000)
##   s2 ~ invgamma(0.001, 0.001)
################################################################################

using MCPhylo

## Data
data = Dict(
  :x => [1, 2, 3, 4, 5],
  :y => [1, 3, 3, 3, 5]
)

## Log-transformed Posterior(b0, b1, log(s2)) + Constant
logf = function(x::DenseVector)
   b0 = x[1]
   b1 = x[2]
   logs2 = x[3]
   r = data[:y] .- b0 .- b1 .* data[:x]
   (-0.5 * length(data[:y]) - 0.001) * logs2 -
     (0.5 * dot(r, r) + 0.001) / exp(logs2) -
     0.5 * b0^2 / 1000 - 0.5 * b1^2 / 1000
end

## MCMC Simulation with Random Walk Metropolis
n = 5000
burnin = 1000
sim = Chains(n, 3, names = ["b0", "b1", "s2"])
theta = RWMVariate([0.0, 0.0, 0.0], [0.5, 0.25, 1.0], logf,
                   proposal = SymUniform)
for i in 1:n
  sample!(theta)
  sim[i, :, 1] = [theta[1:2]; exp(theta[3])]
end
describe(sim)

#mark that we got to the end of the test file succesfully
@test true


text:
# The DeflationOperator implementaion is adapted from BifurcationKit.jl which has the following license:
#=
The BifurcationKit.jl Julia package is licensed under the MIT "Expat" License:
    Copyright (c) 2019: Romain VELTZ.
    Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
    The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#
struct DeflationOperator{T <: Real, Tdot, vectype}
    radius::T
	power::T
	dot::Tdot
	shift::T
	roots::Vector{vectype}
end
function DeflationOperator(
    X; radius = 0.0, power = 2.0,
    dot = LinearAlgebra.dot, shift = 1.0,
)
    return DeflationOperator(radius, power, dot, shift, X)
end
function (df::DeflationOperator{T, Tdot, vectype})(u::vectype) where {T, Tdot, vectype}
	nrm  = u -> df.dot(u, u)
	if length(df.roots) == 0
		return T(1)
	end
	# compute u - df.roots[1]
	tmp = 1 * u
    axpy!(T(-1), df.roots[1], tmp)
	out = T(1) / max(nrm(tmp) - df.radius, 0)^df.power + df.shift
	for ii in 2:length(df.roots)
		copyto!(tmp, u)
        axpy!(T(-1), df.roots[ii], tmp)
		out *= T(1) / max(nrm(tmp) - df.radius, 0)^df.power + df.shift
	end
	return out
end
function ChainRulesCore.rrule(f::DeflationOperator, x)
    return f.def(x), Δ -> (nothing, f.def(x) * pb(Δ)[1])
end

@params struct DeflatedFunction <: Function
    f
    def
end
(f::DeflatedFunction)(x) = f.f(x)
function ChainRulesCore.rrule(f::DeflatedFunction, x)
    v, pb = Zygote.pullback(f.f, x)
    return v, Δ -> (nothing, f.def(x) * pb(Δ)[1])
end

@params struct DeflatedOptions
    ndeflations::Int
    radius::Float64
	power::Float64
	dot::Function
	shift::Float64
    sub_options
end
function DeflatedOptions(;
    ndeflations = 2,
    radius = 0.0,
	power = 2.0,
	dot = LinearAlgebra.dot,
	shift = 1.0,
    sub_options = IpoptOptions(),
    kwargs...,
)
    return DeflatedOptions(
        ndeflations, radius, power, dot, shift, sub_options,
    )
end

struct DeflatedAlg{A} <: AbstractOptimizer
    sub_alg::A
end

function deflatedmodel(vecmodel::VecModel, deflate)
    _obj = vecmodel.objective
    if _obj isa Function
        obj = Objective(
            DeflatedFunction(_obj.f, deflate),
            _obj.multiple,
            _obj.flags,
        )
    else
        obj = _obj
    end
    ineq_constraints = map(vecmodel.ineq_constraints.fs) do c
        @assert c isa IneqConstraint
        return IneqConstraint(
            DeflatedFunction(c, deflate), zero.(c.rhs), c.dim, c.flags,
        )
    end |> VectorOfFunctions
    eq_constraints = map(vecmodel.eq_constraints.fs) do c
        @assert c isa EqConstraint
        return EqConstraint(
            DeflatedFunction(c, deflate), zero.(c.rhs), c.dim, c.flags,
        )
    end |> VectorOfFunctions
    return VecModel(
        obj, eq_constraints, ineq_constraints, vecmodel.box_min,
        vecmodel.box_max, vecmodel.init, vecmodel.integer,
    )
end

@params mutable struct DeflatedWorkspace <: Workspace
    model::VecModel
    defmodel::VecModel
    defop::DeflationOperator
    sub_workspace
    x0::AbstractVector
    X::AbstractVector
    options::DeflatedOptions
end
function Workspace(
    model::VecModel, optimizer::DeflatedAlg, x0::AbstractVector;
    options::DeflatedOptions = DeflatedOptions(), kwargs...,
)
    @unpack radius, shift, power, dot, sub_options = options
    X = Vector{Float64}[]
    defop = DeflationOperator(
        X, radius = radius, shift = shift,
        dot = dot, power = power,
    )
    defmodel = deflatedmodel(model, defop)
    sub_workspace = Workspace(
        defmodel, optimizer.sub_alg, copy(x0);
        options = sub_options, kwargs...,
    )
    return DeflatedWorkspace(
        model, defmodel, defop, sub_workspace, x0, X, options,
    )
end

@params struct DeflatedResult
    solutions
    distances
    minimizer
    minimum
end

function optimize!(workspace::DeflatedWorkspace)
    @unpack options, x0, sub_workspace = workspace
    @unpack model, defop, X = workspace
    solutions = Tuple{Vector{Float64}, Float64}[]
    @unpack ndeflations = options
    sub_workspace.x0 .= x0
    r = optimize!(sub_workspace)
    minimizer = copy(r.minimizer)
    minval = getobjective(model)(minimizer)
    push!(X, copy(minimizer))
    push!(solutions, (copy(minimizer), minval))
    if ndeflations > 0
        i = 0
        while i < ndeflations
            sub_workspace.x0 .= x0
            r = optimize!(sub_workspace)
            objval = getobjective(model)(r.minimizer)
            push!(X, copy(r.minimizer))
            push!(solutions, (copy(r.minimizer), objval))
            if objval <= minval
                minval = objval
                minimizer = copy(r.minimizer)
            end
            i += 1
        end
    end
    distances = [
        norm(solutions[i][1] - solutions[j][1]) for
        i in 1:(ndeflations + 1), j in 1:(ndeflations + 1)
    ]
    return DeflatedResult(
        solutions, distances, minimizer, minval,
    )
end


text:
# Utility operators for computing prox iteratively, e.g. using CG
#
# The idea is that whatever type has the A_mul_B!, Ac_mul_B!, size and eltype
# methods implemented is a linear operator.

import Base: *, size, eltype
import LinearAlgebra: mul!

abstract type LinOp end

function (*)(Op::LinOp, x)
  # Is this the right thing to do?
  # Or maybe just: y = zeros(eltype(x), size(Op, 1))
  y = zeros(promote_type(eltype(Op), eltype(x)), size(Op, 1))
  mul!(y, Op, x)
end

size(Op::LinOp, i::Integer) = i <= 2 ? size(Op)[i] : 1

# AAc (Gram matrix)

struct AAc{M, T} <: LinOp
  A::M
  buf::AbstractArray{T}
  function AAc{M, T}(A::M) where {M, T}
    new(A, zeros(T, size(A, 2)))
  end
end

AAc(A::M) where {M} = AAc{M, eltype(A)}(A)

function mul!(y, Op::AAc, x)
  mul!(Op.buf, adjoint(Op.A), x)
  mul!(y, Op.A, Op.buf)
end

mul!(y, Op::Adjoint{AAc}, x) = mul!(y, adjoint(Op), x)

size(Op::AAc) = size(Op.A, 1), size(Op.A, 1)
eltype(Op::AAc) = eltype(Op.A)

# AcA (Covariance matrix)

struct AcA{M, T} <: LinOp
  A::M
  buf::AbstractArray{T}
  function AcA{M, T}(A::M) where {M, T}
    new(A, zeros(T, size(A, 1)))
  end
end

AcA(A::M) where {M} = AcA{M, eltype(A)}(A)

function mul!(y, Op::AcA, x)
  mul!(Op.buf, Op.A, x)
  mul!(y, adjoint(Op.A), Op.buf)
end

mul!(y, Op::Adjoint{AcA}, x) = mul!(y, adjoint(Op), x)

size(Op::AcA) = size(Op.A, 2), size(Op.A, 2)
eltype(Op::AcA) = eltype(Op.A)

# Shifted symmetric linear operator

struct Shift{M, T} <: LinOp
  A::M
  rho::T
  function Shift{M, T}(A::M, rho::T) where {M, T}
    if eltype(A) != T
      error("type of rho is incompatible with A")
    end
    new(A, rho)
  end
end

Shift(A::M, rho::T) where {M, T} = Shift{M, T}(A, rho)

function mul!(y, Op::Shift, x)
  mul!(y, Op.A, x)
  y .+= Op.rho .* x
end

function mul!(y, Op::Adjoint{Shift}, x)
  mul!(y, adjoint(Op.A), x)
  y .+= Op.rho .* x
end

size(Op::Shift) = size(Op.A, 2), size(Op.A, 2)
eltype(Op::Shift) = eltype(Op.A)


text:
using LinearAlgebra
"""
aer_tor = aerodynamics_torque(C_d, rho, v, a, e_r, e_a, p_c)

空力トルクの計算

# Arguments
- `C_d`：抵抗係数
- `rho`：大気密度
- `v`：大気に対する衛星の相対速度
- `a`：衛星の各面の面積
- `e_r`：衛星の進行方向の単位ベクトル@SCSF（body）
- `e_a`：各面の法線ベクトル@SCSF（body）
- `p_c`：圧力中心の位置ベクトル@SCSF（body）

# Return
- `aer_tor`：空力トルク（1*3）

"""
function aerodynamic_torque(sur_num, C_d, rho, v, a, e_r, e_a, p_c)
	tor = zeros(3)
	cos = zeros(6)
	for i=1:sur_num
		#e_rとe_aのなす角の余弦を計算
		#cos(theta)>0の面が空気抵抗をうける面、cos(theta)<0の面が空気抵抗を受けない面
		ea = e_a[i,:]
		cos[i] = dot(e_r, ea)
		if cos[i]>0.0
			pc = p_c[i,:]
			tor += a[i]*cos[i]*cross(pc, e_r)
		end
	end
	aer_tor = -(1/2)*C_d*rho*v*v*tor

	return aer_tor
end

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
abstract type MatpowerROPFSimpleImaxInput <: AbstractInput end

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

Read instance in `instance_path` depending on `input_type`.\n
Return a structure OPFProblems.
"""
function read_input(input_type::T, instance_path::String) where T<:Type{MatpowerROPFSimpleImaxInput}
  data = load_matpower(instance_path)

  # DataStructure and Gridstructure data:
  bus = SortedDict{String, SortedDict{String, Any}}()
  link = SortedDict{Link, SortedDict{String, Any}}()

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

  checkfor(data, 2, "mpc.baseMVA")
  baseMVA = data[2,3]

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

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

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

    bus_id_line[data[i,1]] = id
    bus_id_name[data[i,1]] = busname
  end

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

  i_debut, i_fin = find_numarray(i_fin+1, data)
  checkfor(data, i_debut-1, "mpc.gen")

  genind, prevgen = 0, 0
  for i=i_debut:i_fin
    gen2bus[i-i_debut+1] = bus_id_line[data[i,1]]
    busname = bus_id_name[data[i,1]]

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

    if data[i, 8] > 0 #generator is on
      S_min += data[i,10] + im*data[i,5] #Smin = Pmin + i Qmin
      S_max += data[i,9] + im*data[i,4] #Smax = Pmax + i Qmax
    end
    bus[busname]["Gen"] = MatpowerGenerator(S_min, S_max, -1, [-1],true)
  end

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

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

  i_debut, i_fin = find_numarray(i_fin+1, data)
  if data[i_debut-1, 1] == "mpc.areas"
    @warn("Not loading mpc.areas data.")
    i_debut = i_fin+1
    i_debut, i_fin = find_numarray(i_fin+1, data)
  end

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

  ## Adding generator cost information
  checkfor(data, i_debut-1, "mpc.gencost")
  genind, cur_genind = 0, data[i_debut, 1]
  for i=i_debut:i_fin
    buslineid, _ = line2busgen[i-i_debut+1]
    busid = bus_id_line[buslineid]

    busname = "BUS_$busid"

    cost_degree = data[i, 4]
    cost_coeffs = data[i, 5:(5+cost_degree-1)]
    if haskey(bus[busname], "Gen")
      bus[busname]["Gen"].cost_degree = cost_degree
      bus[busname]["Gen"].cost_coeffs = cost_coeffs
    else
      @warn("Gen not found at $busname")
    end
  end

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

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

  ds = DataSource(bus, link)

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

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

## Utils functions
function find_numarray(i_start, data)
  i_debut = i_start
  while !isa(data[i_debut, 1], Int)
    i_debut+=1
  end
  i_fin=i_debut
  while !isa(data[i_fin,1], SubString)
    i_fin += 1
  end
  i_debut, i_fin-1
end

checkfor(data, line_ind, name) = (data[line_ind, 1] == name) || error("Expected ", name, " at line ", line_ind, ", got ", data[line_ind,1], " instead.")

function load_matpower(filename)
  instance_name = split(filename, '.')[1]

  touch(instance_name*".temp")
  f = open(filename)
  out = open(instance_name*".temp", "w")

  # removing all ';' at end of lines
  while !eof(f)
    line = readline(f)
    if length(line) > 0 && line[1] != '%' && line[1] != 'f'
      s = split(line, ";")
      println(out, s[1])
    end
  end
  close(f)
  close(out)

  data = DelimitedFiles.readdlm(instance_name*".temp")
  rm(instance_name*".temp")
  return data
end


text:
using Random,JuMP,CPLEX,GLPK

function generarInstancia(periodos,ratioFijoVariable,mu,sigma)
    precioFijo=ratioFijoVariable
    precioStock=1
    d=mu.+randn(periodos).*sigma
    demanda=convert.(Int,round.(d,digits=0))
    #alternativamente:
    #demanda=zeros(Int,periodos)
    #for i in 1:periodos
    #    valor=round(d[i],digits=0)
    #    demanda[i]=convert(Int,d[i])
    #end
    return precioFijo,precioStock,demanda
end

function ulsp(n,K,h,d,s0,useCPLEX) #número de periodos, coste por pedido, coste de inventario, demanda, stock inicial
    if useCPLEX==1
        model = Model(CPLEX.Optimizer)
    else
        model=Model(GLPK.Optimizer)
        set_optimizer_attribute(model,"msg_lev",GLPK.GLP_MSG_ALL)
    end
    @variable(model,y[1:n],Bin) #si se compra
    @variable(model,x[1:n]>=0) #unidades compradas
    @variable(model,s[0:n]>=0) #unidades compradas
    M=sum(d)
    @objective(model, Min, K*sum(y[t] for t in 1:n)+h*sum(s[t] for t in 1:n)) #minimizo costes
    @constraint(model,s0==s[0]) #inventario inicial
    @constraint(model,[t in 1:n],s[t]==s[t-1]+x[t]-d[t]) #inventario
    @constraint(model,[t in 1:n],x[t]<=M*y[t]) #activación
    #println(model)
    optimize!(model)
    println("Objective: ",objective_value(model))
    for i in 1:n
        println("y[$i] = ", JuMP.value(y[i])) #sólo imprima variables de uso de estación
        println("s[$i] = ",JuMP.value(s[i]))
    end
end

function ulspAlt(n,K,h,d,s0,useCPLEX) #número de periodos, coste por pedido, coste de inventario, demanda, stock inicial
    if useCPLEX==1
        model = Model(CPLEX.Optimizer)
    else
        model=Model(GLPK.Optimizer)
        set_optimizer_attribute(model,"msg_lev",GLPK.GLP_MSG_ALL)
    end
    #definimos las variables
    @variable(model,y[1:n],Bin) #si se compra en el periodo i
    @variable(model,0<=x[1:n,1:n]<=1) #$x_{i,j}$  porcentaje de la demanda del periodo j
                                      # satisfecho por compras en el periodo i
    #eliminamos el stock inicial del problema, modificando la demanda
    for i in 1:n
        if d[i]>=s0
            d[i]-=s0
            s0=0
            break
        else
            s0-=d[i]
            d[i]=0
        end
    end
    println(d)
    #creamos los costes de inventario asociados a cubrir parte de la demanda
    c = [0 for i in 1:n, j in 1:n]
    for i in 1:n
        for j in i+1:n
            c[i,j]=d[j]*(j-i)*h
        end
    end
    println(c,"\n\n\n")
    @objective(model, Min, K*sum(y[i] for i in 1:n)+sum(c[i,j]*x[i,j] for i in 1:n, j in i+1:n)) #minimizo costes
    #alternativamente
    #@objective(model,Min,K*sum(y[i] for i in 1:n)+sum( d[i]*h*(j-i)*x[i,j]  for i in 1:n, j in i+1:n))

    @constraint(model,[i in 1:n,j in i:n],x[i,j]<=y[i]) #asegurar condición de costes fijos
    @constraint(model,[j in 1:n],sum(x[i,j] for i in 1:j)==1) #asegurar satisfacción de la demanda del periodo j

    #println(model)
    optimize!(model)
    println("Objective: ",objective_value(model))
    for i in 1:n
        if(JuMP.value(y[i])==1)
            println("pido en el periodo $i")
            for j in i:n
                if (JuMP.value(x[i,j])>0)
                    println(" para periodos $j cantidad ",JuMP.value(x[i,j]))
                end
            end
        end
    end
end

function ulspAlt2(n,K,h,d,s0,useCPLEX) #número de periodos, coste por pedido, coste de inventario, demanda, stock inicial
    if useCPLEX==1
        model = Model(CPLEX.Optimizer)
    else
        model=Model(GLPK.Optimizer)
        set_optimizer_attribute(model,"msg_lev",GLPK.GLP_MSG_ALL)
    end
    #definimos las variables
    #@variable(model,0<=x[1:n,1:n+1]<=1) #si se compra en el periodo i para que nos quedemos sin en el periodo j
    @variable(model,x[1:n,1:n+1],Bin) #si se compra en el periodo i para que nos quedemos sin en el periodo j
    #eliminamos el stock inicial del problema, modificando la demanda
    for i in 1:n
        if d[i]>=s0
            d[i]-=s0
            s0=0
            break
        else
            s0-d[i]
            d[i]=0
        end
    end
    #println(d)
    #creamos los costes totales de cubrir la demanda de los periodos i,...,j-1
    c = [K for i in 1:n, j in 1:n+1]
    for i in 1:n
        for j in i+1:n
            for k in i+1:j
                c[i,j+1] += d[k]*(k-i)*h
            end
        end
    end
    #println(c)
    #alternativamente
    @objective(model,Min,+sum( c[i,j]*x[i,j]  for i in 1:n, j in i+1:n+1))

    @constraint(model,1 == sum(x[1,d] for d in 2:n+1)) #conservación de flujo (salida)
    @constraint(model,[i in 2:n],sum(x[o,i] for o in 1:i) == sum(x[i,d] for d in i+1:n+1)) #conservación de flujo
    @constraint(model,sum(x[i,n+1] for i in 1:n)==1) #conservación de flujo (llegada)

    #println(model)
    optimize!(model)
    println("Objective: ",objective_value(model))
    for i in 1:n
        for j in i:n+1
            if (JuMP.value(x[i,j])>0)
                println(" de $i a $j")
            end
        end
    end
end


text:
nt = 1000
seed = 99
ssm = acqr(0.8,1.5,2,3,0.5,2)
m = acqr_simula(ssm,nt,seed)
y = m["y"]
a = ssm.a
c = ssm.c
q = ssm.q
r = ssm.r
x10 = ssm.m1
P10 = ssm.P1

###########################
print("Testing acqr_kfilter: ")
xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik = acqr_kfilter(y,a,c,q,r,x10,P10)
if sum(xtt.^2) > 1e-6
    println("OK")
end


text:
function import_data3d(data3d, voxelsize_mm, threshold=0)
    """
    Import 3d data to obj vertex and triangulated faces.
    """
    verts, nvoxels = vertices_and_count_voxels_as_the_square_face_preprocessing(data3d, voxelsize_mm, threshold)
    faces = create_square_faces(data3d, nvoxels, threshold)
    new_verts, new_faces = keep_surface_faces(verts, faces)

    new_verts, new_faces = remove_unused_vertexes(new_verts, new_faces)
    trifaces = triangulation(new_faces)

    # new_verts, trifaces = remove_unused_vertexes(new_verts, trifaces)
    check_vf(new_verts, trifaces)
    return new_verts, trifaces
end

"""
Get all vertices and count the voxels
"""
function get_index(sz, i, j, k)
    index = 1 + i + (1 + sz[1]) * j + (sz[1] + 1) * (sz[2] + 1) * k
    return index
end

function vertices_and_count_voxels_as_the_square_face_preprocessing(data3d, voxelsize_mm, threshold=0)
    """

    Fill the area with all possible vertices. Most of them will be removed in future.
    """
    # threshold = 0
    ivertex = 0
    ifaces = 0
    sz = size(data3d)
    verts = Array{Float64}(undef, (sz[1]+1) * (sz[2]+1) * (sz[3]+1), 3)

    # produce vertexes
    nvoxels = 0
    for k in 0:sz[3]
        for j in 0:sz[2]
            for i in 0:sz[1]
                if (i > 0) & (j > 0) & (k > 0)
                    if data3d[i,j,k] > threshold
                        nvoxels += 1
                    end
                end
                pt = get_index(sz, i - 0, j -0, k -0)
    #             println(pt)
                verts[pt, 1] = i * voxelsize_mm[1]
                verts[pt, 2] = j * voxelsize_mm[2]
                verts[pt, 3] = k * voxelsize_mm[3]
            end
        end
    end
    return verts, nvoxels
end

#
#
#       14 -------- 13
#     / |          / |
#   /   |        /   |
# 11 --------- 12    |
# |     |       |    |        i
# |     |       |    |
# |     24 -----|-- 23        |
# |   /         |   /         |   /
# | /           | /           | /
# 21 --------- 22             * ------
#
#
function create_square_faces(data3d, nvoxels, threshold=0)
    """
    create square faces based on vertices from vertices_and_count_voxels_as_...
    :nvoxels: number of voxels higher than the threslold.
    """
    println("threshold ", threshold)
    nfaces = nvoxels * 6
    faces = Array{Int64}(undef, nfaces, 4)
    sz = size(data3d)
    ifaces = 0
#     println(sz)
    # produce faces
    for k in 1:sz[3]
        for j in 1:sz[2]
            for i in 1:sz[1]
    #             println(i,",",j,",",k)
                if data3d[i,j,k] > threshold
    #                 println("jsme uvnitr")
                    pt11 = get_index(sz, i -0, j -0, k -0)
    #                 println("jsme uvnitr", pt11)
                    pt12 = get_index(sz, i -0, j -0, k -1)
                    pt13 = get_index(sz, i -0, j -1, k -1)
                    pt14 = get_index(sz, i -0, j -1, k -0)
                    pt21 = get_index(sz, i -1, j -0, k -0)
                    pt22 = get_index(sz, i -1, j -0, k -1)
                    pt23 = get_index(sz, i -1, j -1, k -1)
                    pt24 = get_index(sz, i -1, j -1, k -0)

                    faces[ifaces + 1, :] = [pt11, pt12, pt13, pt14]'
                    faces[ifaces + 2, :] = [pt22, pt23, pt13, pt12]
                    faces[ifaces + 3, :] = [pt21, pt22, pt12, pt11]
                    faces[ifaces + 4, :] = [pt14, pt13, pt23, pt24]
                    faces[ifaces + 5, :] = [pt11, pt14, pt24, pt21]
                    faces[ifaces + 6, :] = [pt21, pt24, pt23, pt22]
                    ifaces += 6

                end
            end
        end
    end
    return faces
end

function create_cuboid_voxels(data3d, nvoxels, threshold=0)
    """
    create square faces based on vertices from vertices_and_count_voxels_as_...
    :nvoxels: number of voxels higher than the threslold.
    """
    println("threshold ", threshold)
#     nfaces = nvoxels * 6
    cuboids = Array{Int64}(undef, nvoxels)
    sz = size(data3d)
    icuboids = 0
#     println(sz)
    # produce faces
    for k in 1:sz[3]
        for j in 1:sz[2]
            for i in 1:sz[1]
    #             println(i,",",j,",",k)
                if data3d[i,j,k] > threshold
    #                 println("jsme uvnitr")
                    pt11 = get_index(sz, i -1, j -1, k -1)
                    cuboids[icuboids + 1] = pt11
    #                 println("jsme uvnitr", pt11)
#                     pt12 = get_index(sz, i -0, j -0, k -1)
#                     pt13 = get_index(sz, i -0, j -1, k -1)
#                     pt14 = get_index(sz, i -0, j -1, k -0)
#                     pt21 = get_index(sz, i -1, j -0, k -0)
#                     pt22 = get_index(sz, i -1, j -0, k -1)
#                     pt23 = get_index(sz, i -1, j -1, k -1)
#                     pt24 = get_index(sz, i -1, j -1, k -0)

#                     faces[ifaces + 1, :] = [pt11, pt12, pt13, pt14]'
#                     faces[ifaces + 2, :] = [pt22, pt23, pt13, pt12]
#                     faces[ifaces + 3, :] = [pt21, pt22, pt12, pt11]
#                     faces[ifaces + 4, :] = [pt14, pt13, pt23, pt24]
#                     faces[ifaces + 5, :] = [pt11, pt14, pt24, pt21]
#                     faces[ifaces + 6, :] = [pt21, pt24, pt23, pt22]
                    icuboids += 1

                end
            end
        end
    end
    return cuboids
end

function arrayofarray2arrayd2d(d)
    data = Array{Float64}(undef, length(d), length(d[1]))
    for i in 1:length(d)
        for j in 1:length(d[1])
            data[i,j] = d[i][j]
        end
    end
    return data
end

function keep_surface_faces(vertexes, faces, index_base=0)
#     removeDoubleVertexesAndFaces
    # println("step1")
    new_vertexes, inv_vertexes = remove_double_vertexes_alternative(vertexes)# , index_base=index_base)
    # println("step2")
    new_faces = reindexVertexesInFaces(faces, inv_vertexes)
    # println("step3")
    new_faces = removeDoubleFacesByAlberto(faces)
    # println("step1")
    # todo remove unused vertexes
#     new_new_vertexes, new_new_inv_vertexes = remove_double_vertexes_alternative(new_vertexes)# , index_base=index_base)
    # println("step2")
#     new_faces = reindexVertexesInFaces(faces, inv_vertexes)
    new_vertexes = arrayofarray2arrayd2d(new_vertexes)
    return new_vertexes, new_faces
end
#             if segmentation[i,j,k]
function triangulation(faces)
    trifaces = Array{Int64}(undef, size(faces, 1) * 2, 3)
    for i in 1:size(faces,1)
        face = faces[i,:]
        trifaces[(i * 2) - 0, :] = [face[1], face[2], face[3]]
        trifaces[(i * 2) - 1, :] = [face[1], face[3], face[4]]
    end
    return trifaces
end

function thresholding(data3d, threshold)
    sz = size(data3d)
    segmentation = Array{Int8}(undef, size(data3d, 1), size(data3d, 2), size(data3d, 3))
    # produce vertexes
    for k in 1:sz[3]
        for j in 1:sz[2]
            for i in 1:sz[1]
                if data3d[i, j, k] > 0

                    segmentation[i, j, k] = data3d[i, j, k] > threshold
    #             if i == 0
    #             index = 1 + i + sz[1] * j + sz[1] * sz[2] * k
                end

            end
        end
    end
    return segmentation
end


text:
using MergeSorted
import Test: @test, @testset

@testset "Output" begin
    a = sort!(rand(3))
    b = sort!(rand(3))
    c = zeros(6)
    d = mergesorted!(c, a, b)
    e = mergesorted(a,b)
    @test d === c
    @test size(d) == size(e)
end

@testset "Correctness" begin
    for k in [0,1,2,3,4,10,100,1000]
        for d in [0,1,2,3]
            a = sort!(rand(k))
            b = sort!(rand(k+d))
            c = zeros(length(a) + length(b))
            @test sort!(vcat(a,b)) == mergesorted(a,b)
            @test sort!(vcat(a,b)) == mergesorted!(c,a,b)

            a = sort!(rand(k), order=Base.Reverse)
            b = sort!(rand(k+d), order=Base.Reverse)
            c = zeros(length(a) + length(b))
            @test sort!(vcat(a,b), order=Base.Reverse) == mergesorted(a,b, order=Base.Reverse)
            @test sort!(vcat(a,b), order=Base.Reverse) == mergesorted!(c,a,b, order=Base.Reverse)
        end
    end
end


text:
struct OpConnectionIndex{A<:AbstractArray,B<:AbstractArray}
    mel::A
    indices::B
end

function OpConnectionIndex(op::AbsLinearOperator)
    t = conn_type(op)
    v1, v2, v3 = t.parameters
    return OpConnectionIndex(v1(), Vector{Int}())
end

function Base.resize!(c::OpConnectionIndex, n)
    resize!(c.mel,        n)
    resize!(c.indices,  n)
end

Base.eltype(c::OpConnectionIndex{A,B}) where {A,B} = (eltype(A), eltype(B))
Base.length(c::OpConnectionIndex) = length(c.mel)
Base.size(c::OpConnectionIndex) = size(c.mel)

function Base.push!(c::OpConnectionIndex, (m_els, indices))
    push!(c.mel, m_els)
    push!(c.indices, indices)
    return c
end

function Base.append!(c::OpConnectionIndex, (m_els, indices))
    append!(c.mel, m_els)
    append!(c.indices, indices)
    return c
end

Base.getindex(c::OpConnectionIndex, i) = (c.mel[i], c.indices[i])

function Base.iterate(iter::OpConnectionIndex, state=1)
    state > length(iter) && return nothing
    return (iter[state], state + 1)
end


text:
module ThorinDistributions

import Distributions
import Random
import Statistics
import DoubleExponentialFormulas
import LinearAlgebra
import PolynomialRoots
import Combinatorics
import DynamicPolynomials
import StaticPolynomials
import Polynomials
import StatsBase
import SparseArrays

const na = [CartesianIndex()]
include("Precomputations.jl")
include("UnivariateGammaConvolution.jl")
include("MultivariateGamma.jl")
include("MultivariateGammaConvolution.jl")
include("ConvolutionsAndProducts.jl")
include("LaguerreExpensions.jl")
include("MFKProjection.jl")
include("kStats.jl")
include("ThorinMoments.jl")
include("HankelMatrices.jl")
include("Dykstra.jl")

end


text:
using InfiniteOpt, Ipopt, Plots
using KNITRO

# Parmeters
A = 1
B = 1
C = 4
# Time span
t0 = 0
tf = 1
# Initial values
x1 = 1

# Model
model = InfiniteModel(Ipopt.Optimizer)

# infinite_parameter
@infinite_parameter(model, t ∈ [t0, tf], num_supports = 101)

# variable
@variable(model, x, Infinite(t))

@variable(model, u, Infinite(t), start = 0.2)

# objective
@objective(model, Max, ∫(A * x - B * u^2, t))

# constraint
@constraint(model, x(0) == x1)
@constraint(model, x_constr, ∂(x, t) == -0.5 * x^2 + C * u)

# Optimization

print(model)

optimize!(model)

# Visulization
ts = value(t)
u = value(u)
x = value(x)

display(plot(ts, u))
display(plot(ts, x))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

function build_libcurl(ARGS, name::String)
    version = v"7.81.0"
    hash = "ac8e1087711084548d788ef18b9b732c8de887457b81f616fc681d1044b32f98"

    if name == "CURL"
        this_is_curl_jll = true
    elseif name == "LibCURL"
        this_is_curl_jll = false
    else
        msg = "Not a valid name: $(name). Valid names are: LibCURL, CURL"
        throw(ArgumentError(msg))
    end

    # Collection of sources required to build LibCURL
    sources = [
        ArchiveSource("https://curl.se/download/curl-$(version).tar.gz", hash),
    ]

    # Bash recipe for building across all platforms
    script = "THIS_IS_CURL=$(this_is_curl_jll)\n" * raw"""
    cd $WORKSPACE/srcdir/curl-*

    # Holy crow we really configure the bitlets out of this thing
    FLAGS=(
        # Disable....almost everything
        --without-ssl --without-gnutls
        --without-libidn --without-libidn2 --without-librtmp
        --without-nss --without-polarssl
        --without-spnego --without-libpsl --disable-ares --disable-manual
        --disable-ldap --disable-ldaps --without-zsh-functions-dir
        --disable-static --without-libgsasl

        # A few things we actually enable
        --with-libssh2=${prefix} --with-zlib=${prefix} --with-nghttp2=${prefix}
        --enable-versioned-symbols
    )

    if [[ ${target} == *mingw* ]]; then
        # We need to tell it where to find libssh2 on windows
        FLAGS+=(LDFLAGS="${LDFLAGS} -L${prefix}/bin")

        # We also need to tell it to link against schannel (native TLS library)
        FLAGS+=(--with-schannel)
    elif [[ ${target} == *darwin* ]]; then
        # On Darwin, we need to use SecureTransport (native TLS library)
        FLAGS+=(--with-secure-transport)

        # We need to explicitly request a higher `-mmacosx-version-min` here, so that it doesn't
        # complain about: `Symbol not found: ___isOSVersionAtLeast`
        if [[ "${target}" == *x86_64* ]]; then
            export CFLAGS=-mmacosx-version-min=10.11
        fi
    else
        # On all other systems, we use MbedTLS
        FLAGS+=(--with-mbedtls=${prefix})
    fi

    if [[ "${target}" == *linux* ]] || [[ "${target}" == *-freebsd* ]]; then
        # Use gssapi on Linux and FreeBSD
        FLAGS+=(--with-gssapi=${prefix})
        if [[ "${target}" == *-freebsd* ]]; then
            # Only for FreeBSD we need to hint that we need to link to libkrb5 and
            # libcom_err to resolve some undefined symbols.
            export LIBS="-lkrb5 -lcom_err"
        fi
    else
        FLAGS+=(--without-gssapi)
    fi

    ./configure --prefix=$prefix --host=$target --build=${MACHTYPE} "${FLAGS[@]}"
    make -j${nproc}
    if [[ "${THIS_IS_CURL}" == true ]]; then
        # Manually install only `curl`
        install -Dm 755 "src/.libs/curl${exeext}" "${bindir}/curl${exeext}"
    else
        # Install everything...
        make install
        # ...but remove `curl`
        rm "${bindir}/curl${exeext}"
    fi
    install_license COPYING
    """

    # These are the platforms we will build for by default, unless further
    # platforms are passed in on the command line
    platforms = supported_platforms()

    # The products that we will ensure are always built
    if this_is_curl_jll
        # CURL_jll only provides the executable
        products = [
            ExecutableProduct("curl", :curl),
        ]
    else
        # LibCURL only provides the library
        products = [
            LibraryProduct("libcurl", :libcurl),
        ]
    end

    # Dependencies that must be installed before this package can be built
    dependencies = [
        Dependency("LibSSH2_jll"),
        Dependency("Zlib_jll"),
        Dependency("nghttp2_jll"),
        # Note that while we unconditionally list MbedTLS as a dependency,
        # we default to schannel/SecureTransport on Windows/MacOS.
        Dependency("MbedTLS_jll"; compat="~2.28.0", platforms=filter(p->Sys.islinux(p) || Sys.isfreebsd(p), platforms)),
        Dependency("Kerberos_krb5_jll"; platforms=filter(p->Sys.islinux(p) || Sys.isfreebsd(p), platforms)),
    ]

    if this_is_curl_jll
        # Curl_jll depends on LibCURL_jll
        push!(dependencies, Dependency("LibCURL_jll"; compat="$(version)"))
    end

    # Build the tarballs, and possibly a `build.jl` as well.
    build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.8")
end

max_stars_repo_name: MichelJuillard/Yggdrasil

max_issues_repo_name: MichelJuillard/Yggdrasil

max_forks_repo_name: MichelJuillard/Yggdrasil


text:
module InfSites
include("../src/aliases.jl")
include("../src/site_collection.jl")
#include("../src/site_collection_noglobal.jl")
#include("../src/inf_alleles.jl")
include("../src/propsel.jl")
include("../src/watterson.jl")
include("../src/inf_sites.jl")
end

#using InfSites

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
module GuidingCenter4d

    include("guiding_center_4d/solovev_iter.jl")
    include("guiding_center_4d/solovev_iter_xpoint.jl")
    include("guiding_center_4d/solovev_symmetric.jl")
    include("guiding_center_4d/symmetric_quadratic.jl")
    include("guiding_center_4d/theta_pinch.jl")

    include("guiding_center_4d/tokamak_iter_cylindrical.jl")
    include("guiding_center_4d/tokamak_medium_cartesian.jl")
    include("guiding_center_4d/tokamak_medium_cylindrical.jl")
    include("guiding_center_4d/tokamak_small_cartesian.jl")
    include("guiding_center_4d/tokamak_small_cylindrical.jl")
    include("guiding_center_4d/tokamak_small_toroidal.jl")

end


text:
# ### axis range

#md # ![FILE_NAME.png](images/FILE_NAME.png)

plot(sin, 0, 2π, xlims=(-10, 10), ylims=(-2,2))


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
?@generated

# %%
function _staticsum(f, a, b)
    expr = :(f($a))
    for i in a+1:b
        expr = :($expr + f($i))
    end
    expr
end

# %%
_staticsum(:f, -3, 4)

# %%
@generated function staticsum(f, ::Val{a}, ::Val{b}) where {a, b}
    @assert b - a < 100
    expr = :(f($a))
    for i in a+1:b
        expr = :($expr + f($i))
    end
    expr
end

# %%
staticsum(identity, Val(1), Val(10))

# %%
2staticsum(x -> sin(x)/x, Val(1), Val(100)) + 1

# %%
F(n) = 2sum(x -> sin(x)/x, 1:n) + 1
F_static(n) = 2staticsum(x -> sin(x)/x, Val(1), Val(n)) + 1

# %%
using BenchmarkTools
@show F(100) == F_static(100)
@btime F(100)
@btime F_static(100)

# %%
@benchmark F(100)

# %%
@benchmark F_static(100)

# %%
@code_warntype staticsum(x -> sin(x)/x, Val(1), Val(5))

# %%
@code_typed staticsum(x -> sin(x)/x, Val(1), Val(5))

# %%
@code_typed sum(x -> sin(x)/x, 1:5)

# %%

max_stars_repo_path: 0016/@generated function staticsum.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0016/@generated function staticsum.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0016/@generated function staticsum.jl

max_forks_repo_name: genkuroki/public


text:
function solve_waterflooding(si, sj, kr, μw, μo)
    fw(sw) = fractional_flow(sw, kr, μw, μo)
    dfw(sw) = fw_derivative(sw, kr, μw, μo)
    Δfw(sw) = (fw(sw) - fw(si)) / (sw - si)

    S̃  = fzero(sw -> dfw(sw) - Δfw.(sw), sj*0.9)

    WaterFlooding(S̃, si, sj, kr, μw, μo)
end

function saturation_profile(wf::WaterFlooding, time::Real)
    si = wf.Si
    sj = wf.Sj
    s̃ = wf.S̃

    fw(sw) = fractional_flow(sw, wf.kr, wf.μw, wf.μo)
    dfw(sw) = fw_derivative(sw, wf.kr, wf.μw, wf.μo)
    Δfw(sw) = (fw(sw) - fw(si)) / (sw - si)

    s = vcat([si,si], collect(s̃:0.01:sj), [sj])
    x = Vector{Float64}(undef, length(s))
    x[1] = 1.0
    x[2:3] .= Δfw(s̃) * time
    x[findall(s .> s̃)] = dfw(s[s .> s̃]) * time
    x[end] = 0.0

    return x, s
end

# struct WaterFlooding{T <: Real}
#     S̃::T
#     Si::T
#     Sj::T
#     RelPerms::RelPerms
#     μw::T
#     μo::T

#     function WaterFlooding(Si::T, Sj::T, krs::RelPerms,  μw::T, μo::T) where {T <: Real}
#         μw, μo = μ
#         df(s) = fw_derivative(sw, krs, μw, μo)
#         Δfw(sw) = (fw(sw) - fw(si)) / (sw - si)
#         S̃  = find_zero(sw -> dfw(sw) - Δfw.(sw), (si + 0.0001, sj))
#         new(S̃, Si, Sj, krs, μw, μo)
#     end
# end


text:
include("../LF_common/m_LF1d.jl")
include("../LF_common/m_LF3d.jl")

include("init_pot_Hcoul.jl")
include("init_pot_harm_3d.jl")
include("../utils/ortho_gram_schmidt.jl")
include("../utils/orthonormalize.jl")
include("EnergiesT.jl")
include("op_nabla2.jl")
include("op_H.jl")
include("calc_Energies.jl")
include("calc_grad.jl")
include("Sch_solve_Emin_cg.jl")
include("Sch_solve_Emin_pcg.jl")
include("calc_evals.jl")
include("calc_rho.jl")
include("../LF_common/sparse_LF3d.jl")
include("../LF_common/prec_mkl_ilu0.jl")
include("../LF_common/apply_prec_ilu0.jl")
include("diag_lobpcg.jl")

function test_main( ; method = "Emin_cg" )
  # LF parameters
  NN = [25, 25, 25]
  hh = [0.3, 0.3, 0.3]

  Npoints = prod(NN)

  # Initialize LF
  LF = init_LF3d_sinc( NN, hh )
  ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h

  # Parameter for potential
  center = [0.0, 0.0, 0.0]

  ω = 2.0
  Vpot = init_pot_harm_3d( LF, ω, center )
  Ncols = 4

  #Vpot = init_pot_Hcoul( LF, center )
  #Ncols = 1

  if method == "Emin_cg_sparse"
    #
    ∇2 = get_Laplacian3d_kron(LF)
    prec = prec_mkl_ilu0( -0.5*∇2 + spdiagm(Vpot) )
    Energies, evecs = Sch_solve_Emin_pcg( LF, ∇2, prec, Vpot, Ncols, verbose=true )
    evals = calc_evals( LF, Vpot, evecs )
    #
  elseif method == "diag_lobpcg"
    #
    ∇2 = get_Laplacian3d_kron(LF)
    prec = prec_mkl_ilu0( -0.5*∇2 + spdiagm(Vpot) )
    srand(1234)
    v = rand( Npoints, Ncols )
    evals, evecs = diag_lobpcg( LF, ∇2, prec, Vpot, v, verbose=true)
    evecs = evecs/sqrt(ΔV) # normalize evecs
    Energies = calc_Energies( LF, ∇2, Vpot, evecs )
    #
  else
    #
    Energies, evecs = schsolve_Emin_cg( LF, Vpot, Ncols, verbose=true )
    evals = calc_evals( LF, Vpot, evecs )
  end

  @printf("Etot = %18.10f\n", Energies.Total)
  @printf("\nEigenvalues:\n")
  for i = 1:Ncols
    @printf("%8d %f\n", i, evals[i])
  end

  rho = get_rho( evecs )
  @printf( "\nintegRho = %lf\n", sum(rho)*ΔV )

end

@code_native test_main()
@time test_main(method="Emin_cg_sparse")
#@time test_main(method="diag_lobpcg")
#@time test_main(method="Emin_cg")


text:
"""
    ifft_times(fs)

    Returns an array of DFT time points associated with the set of
    frequency-domain points.

    Assumes unit of argument `x` is Hz (not rad/sec). Output unit is in s.
"""
function ifft_times(fs)
    deltaf = fs[2]-fs[1]
    N = Int(length(fs))
    if N%2 == 0
        delta_t = 1/(-2*fs[1])
    else
        delta_t = 1/(2*fs[end]+deltaf)
    end
    ts = (0:N-1)*delta_t
    return ts
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct SnapshotProperties <: SwaggerModel
    timeCreated::Any # spec type: Union{ Nothing, DateTime } # spec name: timeCreated
    osType::Any # spec type: Union{ Nothing, String } # spec name: osType
    hyperVGeneration::Any # spec type: Union{ Nothing, String } # spec name: hyperVGeneration
    creationData::Any # spec type: Union{ Nothing, CreationData } # spec name: creationData
    diskSizeGB::Any # spec type: Union{ Nothing, Int32 } # spec name: diskSizeGB
    diskSizeBytes::Any # spec type: Union{ Nothing, Int64 } # spec name: diskSizeBytes
    uniqueId::Any # spec type: Union{ Nothing, String } # spec name: uniqueId
    encryptionSettingsCollection::Any # spec type: Union{ Nothing, EncryptionSettingsCollection } # spec name: encryptionSettingsCollection
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    incremental::Any # spec type: Union{ Nothing, Bool } # spec name: incremental
    encryption::Any # spec type: Union{ Nothing, Encryption } # spec name: encryption
    networkAccessPolicy::Any # spec type: Union{ Nothing, NetworkAccessPolicy } # spec name: networkAccessPolicy
    diskAccessId::Any # spec type: Union{ Nothing, String } # spec name: diskAccessId

    function SnapshotProperties(;timeCreated=nothing, osType=nothing, hyperVGeneration=nothing, creationData=nothing, diskSizeGB=nothing, diskSizeBytes=nothing, uniqueId=nothing, encryptionSettingsCollection=nothing, provisioningState=nothing, incremental=nothing, encryption=nothing, networkAccessPolicy=nothing, diskAccessId=nothing)
        o = new()
        validate_property(SnapshotProperties, Symbol("timeCreated"), timeCreated)
        setfield!(o, Symbol("timeCreated"), timeCreated)
        validate_property(SnapshotProperties, Symbol("osType"), osType)
        setfield!(o, Symbol("osType"), osType)
        validate_property(SnapshotProperties, Symbol("hyperVGeneration"), hyperVGeneration)
        setfield!(o, Symbol("hyperVGeneration"), hyperVGeneration)
        validate_property(SnapshotProperties, Symbol("creationData"), creationData)
        setfield!(o, Symbol("creationData"), creationData)
        validate_property(SnapshotProperties, Symbol("diskSizeGB"), diskSizeGB)
        setfield!(o, Symbol("diskSizeGB"), diskSizeGB)
        validate_property(SnapshotProperties, Symbol("diskSizeBytes"), diskSizeBytes)
        setfield!(o, Symbol("diskSizeBytes"), diskSizeBytes)
        validate_property(SnapshotProperties, Symbol("uniqueId"), uniqueId)
        setfield!(o, Symbol("uniqueId"), uniqueId)
        validate_property(SnapshotProperties, Symbol("encryptionSettingsCollection"), encryptionSettingsCollection)
        setfield!(o, Symbol("encryptionSettingsCollection"), encryptionSettingsCollection)
        validate_property(SnapshotProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(SnapshotProperties, Symbol("incremental"), incremental)
        setfield!(o, Symbol("incremental"), incremental)
        validate_property(SnapshotProperties, Symbol("encryption"), encryption)
        setfield!(o, Symbol("encryption"), encryption)
        validate_property(SnapshotProperties, Symbol("networkAccessPolicy"), networkAccessPolicy)
        setfield!(o, Symbol("networkAccessPolicy"), networkAccessPolicy)
        validate_property(SnapshotProperties, Symbol("diskAccessId"), diskAccessId)
        setfield!(o, Symbol("diskAccessId"), diskAccessId)
        o
    end
end # type SnapshotProperties

const _property_map_SnapshotProperties = Dict{Symbol,Symbol}(Symbol("timeCreated")=>Symbol("timeCreated"), Symbol("osType")=>Symbol("osType"), Symbol("hyperVGeneration")=>Symbol("hyperVGeneration"), Symbol("creationData")=>Symbol("creationData"), Symbol("diskSizeGB")=>Symbol("diskSizeGB"), Symbol("diskSizeBytes")=>Symbol("diskSizeBytes"), Symbol("uniqueId")=>Symbol("uniqueId"), Symbol("encryptionSettingsCollection")=>Symbol("encryptionSettingsCollection"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("incremental")=>Symbol("incremental"), Symbol("encryption")=>Symbol("encryption"), Symbol("networkAccessPolicy")=>Symbol("networkAccessPolicy"), Symbol("diskAccessId")=>Symbol("diskAccessId"))
const _property_types_SnapshotProperties = Dict{Symbol,String}(Symbol("timeCreated")=>"DateTime", Symbol("osType")=>"String", Symbol("hyperVGeneration")=>"String", Symbol("creationData")=>"CreationData", Symbol("diskSizeGB")=>"Int32", Symbol("diskSizeBytes")=>"Int64", Symbol("uniqueId")=>"String", Symbol("encryptionSettingsCollection")=>"EncryptionSettingsCollection", Symbol("provisioningState")=>"String", Symbol("incremental")=>"Bool", Symbol("encryption")=>"Encryption", Symbol("networkAccessPolicy")=>"NetworkAccessPolicy", Symbol("diskAccessId")=>"String")
Base.propertynames(::Type{ SnapshotProperties }) = collect(keys(_property_map_SnapshotProperties))
Swagger.property_type(::Type{ SnapshotProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_SnapshotProperties[name]))}
Swagger.field_name(::Type{ SnapshotProperties }, property_name::Symbol) =  _property_map_SnapshotProperties[property_name]

const _allowed_SnapshotProperties_osType = ["Windows", "Linux"]

const _allowed_SnapshotProperties_hyperVGeneration = ["V1", "V2"]

function check_required(o::SnapshotProperties)
    (getproperty(o, Symbol("creationData")) === nothing) && (return false)
    true
end

function validate_property(::Type{ SnapshotProperties }, name::Symbol, val)
    if name === Symbol("osType")
        Swagger.validate_param(name, "SnapshotProperties", :enum, val, _allowed_SnapshotProperties_osType)
    end
    if name === Symbol("hyperVGeneration")
        Swagger.validate_param(name, "SnapshotProperties", :enum, val, _allowed_SnapshotProperties_hyperVGeneration)
    end
end


text:
"""
    add_bank!(
        ata_model::AbstractModel;
        settings::InputSettings = InputSettings(),
        bank::DataFrames.DataFrame = DataFrames.DataFrame(),
        settings_file = "settings_ata.jl",
        bank_file = "bank.csv",
        bank_delim = ";",
    )

# Description

Attach an item bank to an initialized ATA model.
The name of a `bank_delim` separated value file can be passed with the argument `bank_file`.
Alternatively, a bank dataframe can be passed using the argument `bank`.

# Arguments

- **`ata_model::AbstractModel`** : Required.
- **`bank::DataFrames.DataFrame`** : Optional. Default: `DataFrame()`. A dataframe containing data about the items.
- **`bank_file`** : Optional. Default: "bank.csv". The path of the file containing the item pool/bank in the form of custom-separated values.
- **`bank_delim`** : Optional. Default: ";". The custom-separator for the bank_file.

# Output

- An ATA model.
"""
function add_bank!(
    ata_model::AbstractModel;
    settings::InputSettings = InputSettings(),
    bank::DataFrames.DataFrame = DataFrames.DataFrame(),
    settings_file = "settings_ata.jl",
    bank_file = "bank.csv",
    bank_delim = ";",
)
    try
        #load bank
        bank_loaded = DataFrame()
        if size(bank, 1) > 0
            bank_loaded = bank
        elseif isfile(bank_file)
            try
                bank_loaded = CSV.read(bank_file, DataFrames.DataFrame, delim = bank_delim)
            catch e
                error!(
                    ata_model,
                    string(
                        "Error in reading the item bank file:\n  ",
                        sprint(showerror, e),
                        ".",
                    ),
                )
            end
        else
            error!(
                ata_model,
                string(
                    "Item bank file with name ",
                    bank_file,
                    " does not exist.\n  Provide a valid item bank dataframe or a name of an existing file.",
                ),
            )
        end
        if size(bank_loaded, 1) == ata_model.settings.n_items
            success!(ata_model, "Item bank file loaded.")
            ata_model.settings.bank = bank_loaded
            _add_irt!(ata_model)
            if ata_model.settings.fs.var != Symbol[]
                _add_friends!(ata_model)
            end
            if ata_model.settings.es.var != Symbol[]
                _add_enemies!(ata_model)
            end
            if any(
                vcat(
                    [
                        ata_model.constraints[t].expected_score.max for
                        t = 1:ata_model.settings.T
                    ]...,
                ) .< 1.00,
            ) || any(
                vcat(
                    [
                        ata_model.constraints[t].expected_score.min for
                        t = 1:ata_model.settings.T
                    ]...,
                ) .> 0.0,
            )
                _add_exp_score!(ata_model)
            end
        else
            error!(
                ata_model,
                string(
                    "Error in loading the item bank dataframe:\n number of rows is not equal to n_items.",
                ),
            )
        end
    catch e
        error!(ata_model, string(sprint(showerror, e)))
    end

    return nothing
end


text:
using FeedbackParticleFilters, Distributions, Random, PDMats

println("Testing DifferentialRKHSMethod.jl:")

@testset "DifferentialRKHSMethod.jl" begin
    print("  solver")
    N = 6
    x_pf=[-2.422480820086937 -0.592332203303167 -2.017301296096984 -1.5151245392598531 0.02565906919199346 0.15161614796874012;]
    testens=UnweightedParticleEnsemble(x_pf)
    eq = PoissonEquation(x->[x[1],exp(x[1])], testens)
    update!(eq, testens)
    solve!(eq, DifferentialRKHSMethod(1E1, 1E-6));
    print(".")
    @test maximum(abs.(eq.gain[1,:,1] - [ 0.1312360106492137
                                       2.343121868385176
                                       0.863411396927212
                                       1.9307957511303933
                                       0.5684889795489233
                                       0.00768453173042769 ])) < 1e-6
    print(".")
    @test maximum(abs.(eq.potential[1,:] - [ -2.072323421422226
                                             1.0382735928503166
                                            -1.8816738883083453
                                            -1.1749266092106723
                                             2.0268744617672745
                                             2.0637758643236275 ])) < 1e-6
    println("DONE.")
end#DifferentialRKHSMethod.jl


text:
using Test
"""
    add(a,b)

Return `a + b`.

Base library use this pattern to document variant of methods and define abstract
function name.
"""
function add end

function add(a,b)
    a+b
end

# A Fix2 defined for add for convenience
"""
    add(b)

Create a function that adds its argument to `b`, i.e.
a function equivalent to `f(x) = add(x,b)`.

Base library use this pattern to document variant of methods.
"""
add(x) = Base.Fix2(add, x)

function learn_Fix2()
    data = [1, 3, 5,  7, 11]
    Ans  = data .+ 3         # [4, 6, 8, 10, 14]

    # Using Fix2, concise
    @test Ans == add(3).(data)
    # Not using Fix2, verbose
    add3(a) = add(a,3)
    @test Ans == add3.(data)

    equal3 = ==(3) # Built-in Fix2
    @test true == equal3(3)
    @test false == equal3(5)
end

@testset "Fix2" begin
    learn_Fix2()
end

nothing


text:
using Test
using AndTraits
using AndTraits: Trait

# Make a trait
Err = Trait{:Err}
Whut = Trait{:Whut}

# Construct a new type
struct MyType{T}
  x::T
end

AndTraits.traits(::MyType{T}) where T = traits(T)

const a = MyType((Err = 10,))
const b = MyType((Whut = 10,))
const c = MyType((Whut = 10, Err = 10))

f(x) = f(traits(x), x)
f(::trait(Err), x) = "Hello"
f(::trait(Whut), x) = 2
f(::trait(Err, Whut), x) = :Howdy

@test f(a) == "Hello"
@test f(b) == 2
@test f(c) == :Howdy

@inferred f(a)
@inferred f(b)
@inferred f(c)


text:
@testset "add capabilities" begin
    config = DockerConfig(;
        base_image = "julia:latest",
        add_capabilities = ["CAP_CHOWN"],
    )
    with_container() do container
        code = """
        println("This was a success.")
        """
        @test success(container, config, `julia -e $(code)`)
    end
end


text:
module Car
    export nb_months

    function nb_months(startPriceOld, startPriceNew, savingperMonth, percentLossByMonth)
        priceDelta = startPriceNew-startPriceOld
        lossFactor = 1 - percentLossByMonth/100
        saving = 0
        month = 0
        while saving < priceDelta
            month += 1
            saving += savingperMonth
            priceDelta *= lossFactor
            if isodd(month)
                lossFactor -= 0.5/100
            end
        end
        [month, round(Int,saving-priceDelta)]
    end

end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using GLM, QuantileRegressions, CSV, Winston, StatsBase

# Load data
url = "http://vincentarelbundock.github.io/Rdatasets/csv/quantreg/engel.csv"
df = CSV.read(download(url))
ResultQR = qreg(@formula(foodexp ~ income), df, .5)

# Fit quantile regression for a bunch of different quantiles
QNum = 35; # Number of equally spaced quantiles
DataPlot = reduce(vcat, [[coeftable(qreg(@formula(foodexp ~ income), df, i/(QNum+1), IP())).cols[2][2] coeftable(qreg(@formula(foodexp ~ income), df, i/(QNum+1), IP())).cols[3][2]] for i in 1:QNum])

# Fit OLS model to compare
ResultLM = lm(@formula(foodexp ~ income), df) # Fit the model using OLS from the GLM-package
ols = vcat([[coeftable(ResultLM).cols[1][2] coeftable(ResultLM).cols[2][2]] for i in 1:QNum]...)

PlotDF = DataFrame(
    X = range(1, length=QNum, stop=QNum)/(QNum+1),
    Y = DataPlot[:,1],
    ols = ols[:,1],
    Ymin = DataPlot[:,1] - 1.96 * DataPlot[:,2],
    Ymax = DataPlot[:,1] + 1.96 * DataPlot[:,2])

x = PlotDF[:X]
p = FramedPlot()
βτ   = Curve(x, PlotDF[:Y])
βols = Curve(x, PlotDF[:ols] , "color", "red")
ci_low  = Curve(x, PlotDF[:Ymin], "type" , "dash")
ci_high = Curve(x, PlotDF[:Ymax], "type" , "dash")
setattr(βτ,"label","β(τ)")
setattr(βols, "label", "β(ols)")
lgnd = Legend(.8,.2,[βτ, βols])
add(p, βτ, βols, ci_low, ci_high, lgnd)
setattr(p, "title" , "Quantile regression: Food Expenditure ~ Income")
setattr(p, "xlabel", "Quantiles")
setattr(p, "ylabel", "Coefficient on Income")
savefig(p, Pkg.dir("QuantileRegressions")"/examples/qreg_example_plot.png")


text:
""""""
function intervalscaled(Fₕ, Fₗ=zeros(size(Fₕ)))
    scale = intervalscale(Fₗ, Fₕ, rampOn) # So, extrapolate if variance is greater than high dim.
    return F -> F |> scale
end
export intervalscaled


text:
using Base.Math: @horner
using Base.MPFR: ROUNDING_MODE
const exparg_n = log(nextfloat(floatmin(Float64)))
const exparg_p =  log(prevfloat(floatmax(Float64)))

#COMPUTE log(gamma(b)/gamma(a+b)) when b >= 8
"""
    loggammadiv(a,b)

Computes ``log(\\Gamma(b)/\\Gamma(a+b))`` when b >= 8
"""
function loggammadiv(a::Float64, b::Float64)
    if a > b
        h = b/a
        c = 1.0/(1.0 + h)
        x = h/(1.0 + h)
        d = a + (b - 0.5)
    else
        h = a/b
        c = h/(1.0 + h)
        x = 1.0/(1.0 + h)
        d = b + a - 0.5
    end
    x² = x*x
    s₃ = 1.0 + (x + x²)
    s₅ = 1.0 + (x + x²*s₃)
    s₇ = 1.0 + (x + x²*s₅)
    s₉ = 1.0 + (x + x²*s₇)
    s₁₁ = 1.0 + (x + x²*s₉)

    # SET W = stirling(b) - stirling(a+b)
    t = inv(b)^2
    w = @horner(t, .833333333333333E-01, -.277777777760991E-02*s₃, .793650666825390E-03*s₅, -.595202931351870E-03*s₇, .837308034031215E-03*s₉, -.165322962780713E-02*s₁₁)
    w *= c/b

    #COMBINING
    u = d*log1p(a/b)
    v = a*(log(b) - 1.0)
    return u <= v ? w - u - v : w - v - u
end

"""
    stirling_corr(a0,b0)

Compute stirling(a0) + stirling(b0) - stirling(a0 + b0)
for a0,b0 >= 8
"""
function stirling_corr(a0::Float64, b0::Float64)
    a = min(a0,b0)
    b = max(a0,b0)

    h = a/b
    c = h/(1.0 + h)
    x = 1.0/(1.0 + h)
    x² = x*x
    #SET SN = (1-X^N)/(1-X)
    s₃ = 1.0 + (x + x²)
    s₅ = 1.0 + (x + x²*s₃)
    s₇ = 1.0 + (x + x²*s₅)
    s₉ = 1.0 + (x + x²*s₇)
    s₁₁ = 1.0 + (x + x²*s₉)
    t = inv(b)^2
    w = @horner(t, .833333333333333E-01, -.277777777760991E-02*s₃, .793650666825390E-03*s₅, -.595202931351870E-03*s₇, .837308034031215E-03*s₉, -.165322962780713E-02*s₁₁)
    w *= c/b
    # COMPUTE stirling(a) + w
    t = inv(a)^2
    return @horner(t, .833333333333333E-01, -.277777777760991E-02, .793650666825390E-03, -.595202931351870E-03, .837308034031215E-03, -.165322962780713E-02)/a + w
end

"""
    esum(mu,x)

Compute ``e^{mu+x}``
"""
function esum(mu::Float64, x::Float64)
    if x > 0.0
        if mu > 0.0 || mu + x < 0.0
            return exp(mu)*exp(x)
        else
            return exp(mu + x)
        end
    elseif mu < 0.0 || mu + x > 0.0
        return exp(mu)*exp(x)
    else
        return exp(mu + x)
    end
end

"""
    beta_integrand(a,b,x,y,mu=0.0)

Compute ``e^{mu} * x^{a}y^{b}/B(a,b)``
"""
function beta_integrand(a::Float64,b::Float64,x::Float64,y::Float64,mu::Float64=0.0)
    a0, b0 = minmax(a,b)
    if a0 >= 8.0
        if a > b
            h = b/a
            x0 = 1.0/(1.0 + h)
            y0 = h/(1.0 + h)
            lambda = (a+b)*y - b
        else
             h = a/b
             x0 = h/(1.0 + h)
             y0 = 1.0/(1.0 + h)
             lambda = a - (a+b)*x
        end
        e = -lambda/a
        u = abs(e) > 0.6 ? u = e - log(x/x0) : - log1pmx(e)
        e = lambda/b
        v = abs(e) > 0.6 ? e - log(y/y0) : - log1pmx(e)
        z = esum(mu, -(a*u + b*v))
        return (1.0/sqrt(2*pi))*sqrt(b*x0)*z*exp(-stirling_corr(a,b))
    elseif x > 0.375
        if y > 0.375
            lnx = log(x)
            lny = log(y)
        else
            lnx = log1p(-y)
            lny = log(y)
        end
    else
        lnx = log(x)
        lny = log1p(-x)
    end
    z = a*lnx + b*lny
    if a0 < 1.0
        b0 = max(a,b)
        if b0 >= 8.0
            u = loggamma1p(a0) + loggammadiv(a0,b0)
            return a0*(esum(mu, z-u))
        elseif b0 > 1.0
            u = loggamma1p(a0)
            n = trunc(Int,b0 - 1.0)
            if n >= 1
                c = 1.0
                for i = 1:n
                    b0 -= 1.0
                    c *= (b0/(a0+b0))
                end
                u += log(c)
            end
            z -= u
            b0 -= 1.0
            apb = a0 + b0
            if apb > 1.0
                u = a0 + b0 - 1.0
                t = (1.0 + rgamma1pm1(u))/apb
            else
                t = 1.0 + rgamma1pm1(apb)
            end
            return a0*(esum(mu,z))*(1.0 + rgamma1pm1(b0))/t
        end
    else
        z -= logbeta(a,b)
        ans = esum(mu,z)
        return ans
    end
    if ans == 0.0
        return 0.0
    end
    apb = a + b
    if apb > 1.0
        z = (1.0 + rgamma1pm1(apb - 1.0))/apb
    else
        z = 1.0 + rgamma1pm1(apb)
    end
    c = (1.0 + rgamma1pm1(a))*(1.0 + rgamma1pm1(b))/z
    return ans*(a0*c)/(1.0 + a0/b0)
end

"""
    beta_inc_cont_fraction(a,b,x,y,lambda,epps)

Compute ``I_{x}(a,b)`` using continued fraction expansion when a,b > 1.
It is assumed that ``\\lambda = (a+b)*y - b``
DLMF : https://dlmf.nist.gov/8.17#E22
BFRAC(A,B,X,Y,LAMBDA,EPS) from Didonato and Morris (1982)
"""
function beta_inc_cont_fraction(a::Float64, b::Float64, x::Float64, y::Float64, lambda::Float64, epps::Float64)
    ans = beta_integrand(a,b,x,y)
    if ans == 0.0
        return 0.0
    end
    c = 1.0 + lambda
    c0 = b/a
    c1 = 1.0 + 1.0/a
    yp1 = y + 1.0

    n = 0.0
    p = 1.0
    s = a + 1.0
    an = 0.0
    bn = 1.0
    anp1 = 1.0
    bnp1 = c/c1
    r = c1/c
    #CONT FRACTION

    while true
     n += 1.0
     t = n/a
     w = n*(b - n)*x
     e = a/s
     alpha = (p*(p+c0)*e*e)*(w*x)
     e = (1.0 + t)/(c1 + 2*t)
     beta = n + w/s +e*(c + n*yp1)
     p = 1.0 + t
     s += 2.0

     #update an, bn, anp1, bnp1
     t = alpha*an  + beta*anp1
     an = anp1
     anp1 = t
     t = alpha*bn + beta*bnp1
     bn = bnp1
     bnp1 = t

     r0 = r
     r = anp1/bnp1
     if abs(r - r0) <= epps*r
        break
     end
     #rescale
     an /= bnp1
     bn /= bnp1
     anp1 = r
     bnp1 = 1.0
    end
    return ans*r
end

"""
    beta_inc_asymptotic_symmetric(a,b,lambda,epps)

Compute ``I_{x}(a,b)`` using asymptotic expansion for a,b >= 15.
It is assumed that ``\\lambda = (a+b)*y - b``
BASYM(A,B,LAMBDA,EPS) from Didonato and Morris (1982)
"""
function beta_inc_asymptotic_symmetric(a::Float64, b::Float64, lambda::Float64, epps::Float64)
    a0 =zeros(22)
    b0 = zeros(22)
    c = zeros(22)
    d = zeros(22)
    e0 = 2/sqrt(pi)
    e1 = 2^(-1.5)
    sm = 0.0
    ans = 0.0
    if a > b
        h = b/a
        r0 = 1.0/(1.0 + h)
        r1 = (b-a)/a
        w0 = 1.0/sqrt(b*(1.0+h))
    else
        h = a/b
        r0 = 1.0/(1.0 + h)
        r1 = (b-a)/b
        w0 = 1.0/sqrt(a*(1.0+h))
    end
    f = -a*log1pmx(-(lambda/a)) - b*log1pmx((lambda/b))
    t = exp(-f)
    if t == 0.0
        return ans
    end
    z0 = sqrt(f)
    z = 0.5*(z0/e1)
    z² = 2.0*f

    a0[1] = (2.0/3.0)*r1
    c[1] = -0.5*a0[1]
    d[1] = - c[1]
    j0 = (0.5/e0)*erfcx(z0)
    j1 = e1
    sm = j0 + d[1]*w0*j1

    s = 1.0
    h² = h*h
    hn = 1.0
    w = w0
    znm1 = z
    zn = z²

    for n = 2: 2: 20
        hn *= h²
        a0[n] = 2.0*r0*(1.0 + h*hn)/(n + 2.0)
        s += hn
        a0[n+1] = 2.0*r1*s/(n+3.0)

        for i = n: n+1
            r = -0.5*(i + 1.0)
            b0[1] = r*a0[1]
            for m = 2:i
                bsum = 0.0
                for j =1: m-1
                    bsum += (j*r - (m-j))*a0[j]*b0[m-j]
                end
                b0[m] = r*a0[m] + bsum/m
            end
            c[i] = b0[i]/(i+1.0)
            dsum = 0.0
            for j = 1: i-1
                imj = i - j
                dsum += d[imj]*c[j]
            end
            d[i] = -(dsum + c[i])
        end

        j0 = e1*znm1 + (n - 1)*j0
        j1 = e1*zn + n*j1
        znm1 *= z²
        zn *= z²
        w *= w0
        t0 = d[n]*w*j0
        w *= w0
        t1 = d[n+1]*w*j1
        sm += (t0 + t1)
        if (abs(t0) + abs(t1)) <= epps*sm
            break
        end
    end

    u = exp(-stirling_corr(a,b))
    return e0*t*u*sm
end

"""
    beta_inc_asymptotic_asymmetric(a,b,x,y,w,epps)

Evaluation of I_{x}(a,b) when b < min(epps,epps*a) and x <= 0.5 using asymptotic expansion.
It is assumed a >= 15 and b <= 1, and epps is tolerance used.
"""
function beta_inc_asymptotic_asymmetric(a::Float64, b::Float64, x::Float64, y::Float64, w::Float64, epps::Float64)
    c = zeros(31)
    d = zeros(31)
    bm1 = b - 1.0
    nu = a + 0.5*bm1
    if y > 0.375
        lnx = log(x)
    else
        lnx = log1p(-y)
    end
    z = -nu*lnx
    if b*z == 0.0
        return error("expansion can't be computed")
    end

    # COMPUTATION OF THE EXPANSION
    #SET R = EXP(-Z)*Z**B/GAMMA(B)
    r = b*(1.0 + rgamma1pm1(b))*exp(b*log(z))
    r *= exp(a*lnx)*exp(0.5*bm1*lnx)
    u = loggammadiv(b,a) + b*log(nu)
    u = r*exp(-u)
    if u == 0.0
        return error("expansion can't be computed")
    end
    (p, q) = gamma_inc(b,z,0)
    v = inv(nu)^2/4
    t2 = lnx^2/4
    l = w/u
    j = q/r
    sm = j
    t = 1.0
    cn = 1.0
    n2 = 0.0
    for n = 1:30
        bp2n = b + n2
        j = (bp2n*(bp2n + 1.0)*j + (z + bp2n + 1.0)*t)*v
        n2 += 2.0
        t *= t2
        cn /= n2*(n2 + 1.0)
        c[n] = cn
        s = 0.0
        if n != 1
            nm1 = n -1
            coef = b - n
            for i = 1:nm1
                s += coef*c[i]*d[n-i]
                coef += b
            end
        end
        d[n] = bm1*cn + s/n
        dj = d[n] * j
        sm += dj
        if sm <= 0.0
            return error("expansion can't be computed")
        end
        if abs(dj) <= epps*(sm+l)
            break
        end
    end
    return w + u*sm
end

#For b < min(eps, eps*a) and x <= 0.5
"""
    beta_inc_power_series2(a,b,x,epps)

Variant of BPSER(A,B,X,EPS).
FPSER(A,B,X,EPS) from Didonato and Morris (1982)
"""
function beta_inc_power_series2(a::Float64, b::Float64, x::Float64, epps::Float64)
    ans = 1.0
    if a > 1.0e-3 * epps
        ans = 0.0
        t = a*log(x)
        if t < exparg_n
            return ans
        end
        ans = exp(t)
    end
    ans *= b/a
    tol = epps/a
    an = a + 1.0
    t = x
    s = t/an
    an += 1.0
    t *= x
    c = t/an
    s += c
    while abs(c) > tol
        an += 1.0
        t *= x
        c = t/an
        s += c
    end
    ans *= (1.0 + a*s)
    return ans
end

#A <= MIN(EPS,EPS*B), B*X <= 1, AND X <= 0.5.,  A is small
"""
    beta_inc_power_series1(a,b,x,epps)

Another variant of BPSER(A,B,X,EPS)
APSER(A,B,X,EPS) from Didonato and Morris (1982)
"""
function beta_inc_power_series1(a::Float64, b::Float64, x::Float64, epps::Float64)
    g = Base.MathConstants.eulergamma
    bx = b*x
    t = x - bx
    if b*epps > 2e-2
        c = log(bx) + g + t
    else
        c = log(x) + digamma(b) + g + t
    end
    tol = 5.0*epps*abs(c)
    j = 1.0
    s = 0.0
    j += 1.0
    t *= (x - bx/j)
    aj = t/j
    s += aj
    while abs(aj) > tol
       j += 1.0
       t *= (x - bx/j)
       aj = t/j
       s += aj
    end
    return -a*(c + s)
end

#B .LE. 1 OR B*X .LE. 0.7
"""
    beta_inc_power_series(a,b,x,epps)

Computes Ix(a,b) using power series :
```math
I_{x}(a,b) = G(a,b)x^{a}/a (1 + a\\sum_{1}^{\\infty}((1-b)(2-b)...(j-b)/j!(a+j)) x^{j})
```
BPSER(A,B,X,EPS) from Didonato and Morris (1982)
"""
function beta_inc_power_series(a::Float64, b::Float64, x::Float64, epps::Float64)
    ans = 0.0
    if x == 0.0
        return 0.0
    end
    a0 = min(a,b)
    b0 = max(a,b)
    if a0 >= 1.0
        z = a*log(x) - logbeta(a,b)
        ans = exp(z)/a
    else

        if b0 >= 8.0
            u = loggamma1p(a0) + loggammadiv(a0,b0)
            z = a*log(x) - u
            ans = (a0/a)*exp(z)
            if ans == 0.0 || a <= 0.1*epps
                return ans
            end
        elseif b0 > 1.0
            u = loggamma1p(a0)
            m = b0 - 1.0
            if m >= 1.0
                c = 1.0
                for i = 1:m
                    b0 -= 1.0
                    c *= (b0/(a0+b0))
                end
                u += log(c)
            end
            z = a*log(x) - u
            b0 -= 1.0
            apb = a0 + b0
            if apb > 1.0
                u = a0 + b0 - 1.0
                t = (1.0 + rgamma1pm1(u))/apb
            else
                t = 1.0 + rgamma1pm1(apb)
            end
            ans = exp(z)*(a0/a)*(1.0 + rgamma1pm1(b0))/t
            if ans == 0.0 || a <= 0.1*epps
                return ans
            end
        else
        #PROCEDURE FOR A0 < 1 && B0 < 1
            ans = x^a
            if ans == 0.0
                return ans
            end
            apb = a + b
            if apb > 1.0
                u = a + b - 1.0
                z = (1.0 + rgamma1pm1(u))/apb
            else
                z = 1.0 + rgamma1pm1(apb)
            end
            c = (1.0 + rgamma1pm1(a))*(1.0 + rgamma1pm1(b))/z
            ans *= c*(b/apb)
            #label l70 start
            if ans == 0.0 || a <= 0.1*epps
                return ans
            end
        end
    end
    if ans == 0.0 || a <= 0.1*epps
        return ans
    end
    # COMPUTE THE SERIES

    sm = 0.0
    n = 0.0
    c = 1.0
    tol = epps/a
    n += 1.0
    c *= x*(1.0 - b/n)
    w = c/(a + n)
    sm += w
    while abs(w) > tol
        n += 1.0
        c *= x*(1.0 - b/n)
        w = c/(a+n)
        sm += w
    end
    return ans*(1.0 + a*sm)
end

"""
    beta_inc_diff(a,b,x,y,n,epps)

Compute ``I_{x}(a,b) - I_{x}(a+n,b)`` where n is positive integer and epps is tolerance.
A more generalised version of https://dlmf.nist.gov/8.17#E20
"""
function beta_inc_diff(a::Float64, b::Float64, x::Float64, y::Float64, n::Integer, epps::Float64)
    apb = a + b
    ap1 = a + 1.0
    mu = 0.0
    d = 1.0
    if n != 1 && a >= 1.0 && apb >= 1.1*ap1
        mu = abs(exparg_n)
        k = exparg_p
        if k < mu
            mu = k
        end
        t = mu
        d = exp(-t)
    end

    ans = beta_integrand(a,b,x,y,mu)/a
    if n == 1 || ans == 0.0
        return ans
    end
    nm1 = n -1
    w = d

    k = 0
    if b <= 1.0
        kp1 = k + 1
        for i = kp1:nm1
            l = i - 1
            d *= ((apb + l)/(ap1 + l))*x
            w += d
            if d <= epps*w
                break
            end
        end
        return ans*w
    elseif y > 1.0e-4
        r = trunc(Int,(b - 1.0)*x/y - a)
        if r < 1.0
            kp1 = k + 1
            for i = kp1:nm1
                l = i - 1
                d *= ((apb + l)/(ap1 + l))*x
                w += d
                if d <= epps*w
                    break
                end
            end
            return ans*w
        end
        k = t = nm1
        if r < t
            k = r
        end
        # ADD INC TERMS OF SERIES
        for i = 1:k
            l = i -1
            d *= ((apb + l)/(ap1 + l))*x
            w += d
        end
        if k == nm1
            return ans*w
        end
    else
        k = nm1
        for i = 1:k
            l = i -1
            d *= ((apb + l)/(ap1 + l))*x
            w += d
        end
        if k == nm1
            return ans*w
        end
    end

    kp1 = k + 1
    for i = kp1:nm1
        l = i - 1
        d *= ((apb + l)/(ap1 + l))*x
        w += d
        if d <= epps*w
            break
        end
   end
   return ans*w
end

#SIGNIFICANT DIGIT COMPUTATION OF INCOMPLETE BETA FUNCTION RATIOS
#by ARMIDO R. DIDONATO AND ALFRED H. MORRIS, JR.
#ACM Transactions on Mathematical Software. Vol18, No3, September1992, Pages360-373
#DLMF : https://dlmf.nist.gov/8.17#E1
#Wikipedia : https://en.wikipedia.org/wiki/Beta_function#Incomplete_beta_function

"""
    beta_inc(a,b,x,y)

Computes Incomplete Beta Function Ratios given by:
```math
I_{x}(a,b) = G(a,b) \\int_{0}^{x} t^{a-1}(1-t)^{b-1} dt,
```
and ``I_{y}(a,b) = 1.0 - I_{x}(a,b)``.
given
``B(a,b) = 1/G(a,b) = \\Gamma(a)\\Gamma(b)/\\Gamma(a+b)`` and ``x+y = 1``.
"""
function beta_inc(a::Float64, b::Float64, x::Float64, y::Float64)
    p = 0.0
    q = 0.0
   # lambda = a - (a+b)*x
    if a < 0.0 || b < 0.0
        return error("a or b is negative")
    elseif a == 0.0 && b == 0.0
        return error("a and b are 0.0")
    elseif x < 0.0 || x > 1.0
        return error("x < 0 or x > 1")
    elseif y < 0.0 || y > 1.0
        return error("y < 0  or y > 1")
    else
        z = x + y - 1.0
        if abs(z) > 3.0*eps()
            return error("x + y != 1.0")         # ERROR HANDLING
        end
    end

    if x == 0.0
        return (0.0,1.0)
    elseif y == 0.0
        return (1.0,0.0)
    elseif a == 0.0
        return (1.0,0.0)
    elseif b == 0.0
        return (0.0,1.0)
    end
#EVALUATION OF ALGOS FOR PROPER SUB-DOMAINS ABOVE
    epps = max(eps(), 1.0e-15)
    if max(a,b) < 1.0E-3 * epps
        return (b/(a+b), a/(a+b))
    end
    ind = false
    a0 = a
    b0 = b
    x0 = x
    y0 = y

    if min(a0,b0) > 1.0
        #PROCEDURE FOR A0>1 AND B0>1
        lambda = a > b ? (a+b)*y - b : a - (a+b)*x
        if lambda < 0.0
            ind = true
            a0 = b
            b0 = a
            x0 = y
            y0 = x
            lambda = abs(lambda)
        end
        if b0 < 40.0 && b0*x0 <= 0.7
            p = beta_inc_power_series(a0,b0,x0,epps)
            q = 1.0 - p
        elseif b0 < 40.0
                n = trunc(Int, b0)
                b0 -= float(n)
                if b0 == 0.0
                    n-=1
                    b0=1.0
                end
                p = beta_inc_diff(b0,a0,y0,x0,n,epps)
                if x0 <= 0.7
                    p += beta_inc_power_series(a0,b0,x0,epps)
                    q = 1.0 - p
                else
                    if a0 <= 15.0
                        n = 20
                        p += beta_inc_diff(a0, b0, x0, y0, n, epps)
                        a0 += n
                    end
                    p = beta_inc_asymptotic_asymmetric(a0,b0,x0,y0,p,15.0*eps())
                    q = 1.0 - p
                end
        elseif a0 > b0
            if b0 <= 100.0 || lambda > 0.03*b0
                p = beta_inc_cont_fraction(a0,b0,x0,y0,lambda,15.0*eps())
                q = 1.0 - p
            else
                p = beta_inc_asymptotic_symmetric(a0,b0,lambda,100.0*eps())
                q = 1.0 - p
            end
        elseif a0 <= 100.0 || lambda > 0.03*a0
            p = beta_inc_cont_fraction(a0,b0,x0,y0,lambda,15.0*eps())
            q = 1.0 - p
        else
            p = beta_inc_asymptotic_symmetric(a0,b0,lambda,100.0*eps())
            q = 1.0 - p
        end
        return ind ? (q, p) : (p, q)
    end
#PROCEDURE FOR A0<=1 OR B0<=1
    if x > 0.5
        ind = true
        a0 = b
        b0 = a
        y0 = x
        x0 = y
    end

    if b0 < min(epps, epps*a0)
        p = beta_inc_power_series2(a0,b0,x0,epps)
        q = 1.0 - p
    elseif a0 < min(epps, epps*b0) && b0*x0 <= 1.0
        q = beta_inc_power_series1(a0,b0,x0,epps)
        p = 1.0 - q
    elseif max(a0,b0) > 1.0
        if b0 <= 1.0
            p = beta_inc_power_series(a0,b0,x0,epps)
            q = 1.0 - p
        elseif x0 >= 0.3
            q = beta_inc_power_series(b0,a0,y0,epps)
            p = 1.0 - q
        elseif x0 >= 0.1
            if b0 > 15.0
                q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps())
                p = 1.0 - q
            else
                n = 20
                q = beta_inc_diff(b0,a0,y0,x0,n,epps)
                b0 += n
                q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps())
                p = 1.0 - q
            end
        elseif (x0*b0)^(a0) <= 0.7
            p = beta_inc_power_series(a0,b0,x0,epps)
            q = 1.0 - p
        else
            n = 20
            q = beta_inc_diff(b0,a0,y0,x0,n,epps)
            b0 += n
            q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps())
            p = 1.0 - q
        end
    elseif a0 >= min(0.2, b0)
        p = beta_inc_power_series(a0,b0,x0,epps)
        q = 1.0 - p
    elseif x0^a0 <= 0.9
        p = beta_inc_power_series(a0,b0,x0,epps)
        q = 1.0 - p
    elseif x0 >= 0.3
        q = beta_inc_power_series(b0,a0,y0,epps)
        p = 1.0 - q
    else
        n = 20
        q = beta_inc_diff(b0,a0,y0,x0,n,epps)
        b0 += n
        q = beta_inc_asymptotic_asymmetric(b0,a0,y0,x0,q,15.0*eps())
        p = 1.0 - q
    end

#TERMINATION
    return ind ? (q, p) : (p, q)
end

beta_inc(a::Float64, b::Float64, x::Float64) = beta_inc(a, b, x, 1.0 - x)
function beta_inc(a::T, b::T, x::T, y::T) where {T<:Union{Float16, Float32}}
    T.(beta_inc(Float64(a), Float64(b), Float64(x), Float64(y)))
end
beta_inc(a::Real, b::Real, x::Real, y::Real) = beta_inc(promote(float(a), float(b), float(x), float(y))...)
beta_inc(a::T, b::T, x::T, y::T) where {T<:AbstractFloat} = throw(MethodError(beta_inc,(a, b, x, y,"")))