text:
import LsqFit.curve_fit
function curve_fit(DM::AbstractDataModel, initial::AbstractVector{<:Number}=MLE(DM), LogPriorFn::Union{Nothing,Function}=LogPrior(DM); tol::Real=1e-14, kwargs...)
    curve_fit(Data(DM), Predictor(DM), dPredictor(DM), initial, LogPriorFn; tol=tol, kwargs...)
end

function curve_fit(DS::AbstractDataSet, M::ModelMap, initial::AbstractVector{<:Number}=GetStartP(DS,M), LogPriorFn::Union{Nothing,Function}=nothing; tol::Real=1e-14, kwargs...)
    curve_fit(DS, M.Map, initial, LogPriorFn; tol=tol, lower=convert(Vector,M.Domain.L), upper=convert(Vector,M.Domain.U), kwargs...)
end

function curve_fit(DS::AbstractDataSet, M::ModelMap, dM::ModelOrFunction, initial::AbstractVector{<:Number}=GetStartP(DS,M), LogPriorFn::Union{Nothing,Function}=nothing; tol::Real=1e-14, kwargs...)
    curve_fit(DS, M.Map, dM, initial, LogPriorFn; tol=tol, lower=convert(Vector,M.Domain.L), upper=convert(Vector,M.Domain.U), kwargs...)
end

function curve_fit(DS::AbstractDataSet, model::Function, initial::AbstractVector{<:Number}=GetStartP(DS,model), LogPriorFn::Union{Nothing,Function}=nothing; tol::Real=1e-14, kwargs...)
    !isnothing(LogPriorFn) && @warn "curve_fit() cannot account for priors. Throwing away given prior and continuing anyway."
    X = xdata(DS);  Y = ydata(DS);    LsqFit.check_data_health(X, Y)
    u = cholesky(yInvCov(DS)).U
    f(θ::AbstractVector) = u * (EmbeddingMap(DS, model, θ) - Y)
    p0 = convert(Vector, initial)
    R = LsqFit.OnceDifferentiable(f, p0, copy(f(p0)); inplace = false, autodiff = :forward)
    LsqFit.lmfit(R, p0, yInvCov(DS); x_tol=tol, g_tol=tol, kwargs...)
end

function curve_fit(DS::AbstractDataSet, model::Function, dmodel::ModelOrFunction, initial::AbstractVector{<:Number}=GetStartP(DS,model), LogPriorFn::Union{Nothing,Function}=nothing; tol::Real=1e-14, kwargs...)
    !isnothing(LogPriorFn) && @warn "curve_fit() cannot account for priors. Throwing away given prior and continuing anyway."
    X = xdata(DS);  Y = ydata(DS);    LsqFit.check_data_health(X, Y)
    u = cholesky(yInvCov(DS)).U
    f(θ::AbstractVector) = u * (EmbeddingMap(DS, model, θ) - Y)
    df(θ::AbstractVector) = u * EmbeddingMatrix(DS, dmodel, θ)
    p0 = convert(Vector, initial)
    R = LsqFit.OnceDifferentiable(f, df, p0, copy(f(p0)); inplace = false)
    LsqFit.lmfit(R, p0, yInvCov(DS); x_tol=tol, g_tol=tol, kwargs...)
end

function normalizedjac(M::AbstractMatrix{<:Number}, xlen::Int)
    M[:,1:xlen] .*= sqrt(size(M,1)/xlen -1.);    return M
end

function TotalLeastSquares(DM::AbstractDataModel, args...; kwargs...)
    !isnothing(LogPrior(DM)) && @warn "TotalLeastSquares() cannot account for priors. Throwing away given prior and continuing anyway."
    TotalLeastSquares(Data(DM), Predictor(DM), args...; kwargs...)
end
"""
    TotalLeastSquares(DSE::DataSetExact, model::ModelOrFunction, initial::AbstractVector{<:Number}; tol::Real=1e-13, kwargs...) -> Vector
Experimental feature which takes into account uncertainties in x-values to improve the accuracy of the fit.
Returns concatenated vector of x-values and parameters. Assumes that the uncertainties in the x-values and y-values are normal, i.e. Gaussian!
"""
function TotalLeastSquares(DSE::DataSetExact, model::ModelOrFunction, initialp::AbstractVector{<:Number}=GetStartP(DSE, model); ADmode::Union{Symbol,Val}=Val(:ForwardDiff), tol::Real=1e-13, rescale::Bool=true, kwargs...)
    # Improve starting values by fitting with ordinary least squares first
    initialp = curve_fit(DataSet(WoundX(DSE),Windup(ydata(DSE),ydim(DSE)),ysigma(DSE)), model, initialp; tol=tol, kwargs...).param
    if xdist(DSE) isa InformationGeometry.Dirac
        println("xdist of given data is Dirac, can only use ordinary least squares.")
        return xdata(DSE), initialp
    end

    plen = pdim(DSE,model);  xlen = Npoints(DSE) * xdim(DSE)
    function predictY(ξ::AbstractVector)
        x = view(ξ, 1:xlen);        p = view(ξ, (xlen+1):length(ξ))
        # INPLACE EmbeddingMap!() would be great here!
        vcat(x, EmbeddingMap(DSE, model, p, Windup(x,xdim(DSE))))
    end
    u = cholesky(BlockMatrix(InvCov(xdist(DSE)),InvCov(ydist(DSE)))).U;    Ydata = vcat(xdata(DSE), ydata(DSE))
    f(θ::AbstractVector) = u * (predictY(θ) - Ydata)
    Jac = GetJac(ADmode, predictY)
    dfnormalized(θ::AbstractVector) = u * normalizedjac(Jac(θ), xlen)
    df(θ::AbstractVector) = u * Jac(θ)
    p0 = vcat(xdata(DSE), initialp)
    R = rescale ? LsqFit.OnceDifferentiable(f, dfnormalized, p0, copy(f(p0)); inplace = false) : LsqFit.OnceDifferentiable(f, df, p0, copy(f(p0)); inplace = false)
    fit = LsqFit.lmfit(R, p0, BlockMatrix(InvCov(xdist(DSE)), InvCov(ydist(DSE))); x_tol=tol, g_tol=tol, kwargs...)
    Windup(fit.param[1:xlen],xdim(DSE)), fit.param[xlen+1:end]
end

## Remove above custom method and implement LiftedCost∘LiftedEmbedding for both
function TotalLeastSquares(DS::AbstractDataSet, model::ModelOrFunction, initialp::AbstractVector{<:Number}=GetStartP(DS, model); tol::Real=1e-13, kwargs...)
    sum(abs, xsigma(DS)) == 0.0 && throw("Cannot perform Total Least Squares Fitting for DataSets without x-uncertainties.")
    xlen = Npoints(DS)*xdim(DS);    Cost(x::AbstractVector) = -logpdf(dist(DS), x)
    function predictY(ξ::AbstractVector)
        x = view(ξ, 1:xlen);        p = view(ξ, (xlen+1):length(ξ))
        vcat(x, EmbeddingMap(DS, model, p, Windup(x,xdim(DS))))
    end
    InformationGeometry.minimize(Cost∘predictY, [xdata(DS); initialp]; tol=tol, kwargs...)
end

"""
    minimize(F::Function, start::AbstractVector{<:Number}; tol::Real=1e-10, meth=NelderMead(), Full::Bool=false, timeout::Real=200, kwargs...) -> Vector
Minimizes the scalar input function using the given `start` using algorithms from `Optim.jl` specified via the keyword `meth`.
`Full=true` returns the full solution object instead of only the minimizing result.
Optionally, the search domain can be bounded by passing a suitable `HyperCube` object as the third argument.
"""
function minimize(F::Function, start::AbstractVector{<:Number}, Domain::Union{HyperCube,Nothing}=nothing; Fthresh::Union{Nothing,Real}=nothing, tol::Real=1e-10, meth::Optim.AbstractOptimizer=NelderMead(), timeout::Real=200, Full::Bool=false, kwargs...)
    !(F(start) isa Number) && throw("Given function must return scalar values, got $(typeof(F(start))) instead.")
    options = if isnothing(Fthresh)
        Optim.Options(g_tol=tol, x_tol=tol, time_limit=floatify(timeout))
    else  # stopping criterion via callback kwarg
        Optim.Options(callback=(z->z.value<Fthresh), g_tol=tol, x_tol=tol, time_limit=floatify(timeout))
    end
    Res = if isnothing(Domain)
        Optim.optimize(F, floatify(start), meth, options; kwargs...)
    else
        start ∉ Domain && throw("Given starting value not in specified domain.")
        Optim.optimize(F, convert(Vector{Float64},Domain.L), convert(Vector{Float64},Domain.U), floatify(start), meth, options; kwargs...)
    end
    Full ? Res : Optim.minimizer(Res)
end
minimize(FdF::Tuple, args...; kwargs...) = InformationGeometry.minimize(FdF..., args...; kwargs...)
function minimize(F::Function, dF::Function, start::AbstractVector{<:Number}, Domain::Union{HyperCube,Nothing}=nothing; Fthresh::Union{Nothing,Real}=nothing, tol::Real=1e-10, meth::Optim.AbstractOptimizer=BFGS(), timeout::Real=200, Full::Bool=false, kwargs...)
    !(F(start) isa Number) && throw("Given function must return scalar values, got $(typeof(F(start))) instead.")
    options = if isnothing(Fthresh)
        Optim.Options(g_tol=tol, x_tol=tol, time_limit=floatify(timeout))
    else  # stopping criterion via callback kwarg
        Optim.Options(callback=(z->z.value<Fthresh), g_tol=tol, x_tol=tol, time_limit=floatify(timeout))
    end
    Res = if isnothing(Domain)
        Optim.optimize(F, dF, floatify(start), meth, options; inplace=MaximalNumberOfArguments(dF)>1, kwargs...)
    else
        start ∉ Domain && throw("Given starting value not in specified domain.")
        Optim.optimize(F, dF, convert(Vector{Float64},Domain.L), convert(Vector{Float64},Domain.U), floatify(start), meth, options; inplace=MaximalNumberOfArguments(dF)>1, kwargs...)
    end
    Full ? Res : Optim.minimizer(Res)
end
function minimize(F::Function, dF::Function, ddF::Function, start::AbstractVector{<:Number}, Domain::Union{HyperCube,Nothing}=nothing; Fthresh::Union{Nothing,Real}=nothing, tol::Real=1e-10, meth::Optim.AbstractOptimizer=NewtonTrustRegion(), timeout::Real=200, Full::Bool=false, kwargs...)
    !(F(start) isa Number) && throw("Given function must return scalar values, got $(typeof(F(start))) instead.")
    @assert MaximalNumberOfArguments(dF) == MaximalNumberOfArguments(ddF) "Derivatives dF and ddF need to be either both in-place or both not in-place"
    options = if isnothing(Fthresh)
        Optim.Options(g_tol=tol, x_tol=tol, time_limit=floatify(timeout))
    else  # stopping criterion via callback kwarg
        Optim.Options(callback=(z->z.value<Fthresh), g_tol=tol, x_tol=tol, time_limit=floatify(timeout))
    end
    Res = if isnothing(Domain)
        Optim.optimize(F, dF, ddF, floatify(start), meth, options; inplace=MaximalNumberOfArguments(dF)>1, kwargs...)
    else
        start ∉ Domain && throw("Given starting value not in specified domain.")
        Optim.optimize(F, dF, ddF, convert(Vector{Float64},Domain.L), convert(Vector{Float64},Domain.U), floatify(start), meth, options; inplace=MaximalNumberOfArguments(dF)>1, kwargs...)
    end
    Full ? Res : Optim.minimizer(Res)
end

"""
    RobustFit(DM::AbstractDataModel, start::AbstractVector{<:Number}; tol::Real=1e-10, p::Real=1, kwargs...)
Uses `p`-Norm to judge distance on Dataspace as specified by the keyword.
"""
RobustFit(DM::AbstractDataModel, start::AbstractVector{<:Number}=MLE(DM); kwargs...) = RobustFit(Data(DM), Predictor(DM), start, LogPrior(DM); kwargs...)
function RobustFit(DS::AbstractDataSet, M::ModelOrFunction, start::AbstractVector{<:Number}=GetStartP(DS,M), LogPriorFn::Union{Nothing,Function}=nothing; Domain::Union{HyperCube,Nothing}=(M isa ModelMap ? M.Domain : nothing), tol::Real=1e-10, p::Real=1, kwargs...)
    HalfSig = cholesky(yInvCov(DS)).U
    # Since F is minimized, need to subtract LogPrior
    F(θ::AbstractVector) = norm(HalfSig * (ydata(DS) - EmbeddingMap(DS, M, θ)), p) - EvalLogPrior(LogPriorFn, θ)
    InformationGeometry.minimize(F, start, Domain; tol=tol, kwargs...)
end
function RobustFit(DS::AbstractDataSet, M::ModelOrFunction, dM::ModelOrFunction, start::AbstractVector{<:Number}=GetStartP(DS,M), LogPriorFn::Union{Nothing,Function}=nothing; Domain::Union{HyperCube,Nothing}=(M isa ModelMap ? M.Domain : nothing), tol::Real=1e-10, p::Real=1, kwargs...)
    HalfSig = cholesky(yInvCov(DS)).U
    # Since F is minimized, need to subtract LogPrior
    F(θ::AbstractVector) = norm(HalfSig * (EmbeddingMap(DS, M, θ) - ydata(DS)), p) - EvalLogPrior(LogPriorFn, θ)
    function dFp(θ::AbstractVector)
        z = HalfSig * (EmbeddingMap(DS, M, θ) - ydata(DS))
        n = sum(z.^p)^(1/p - 1) * z.^(p-1)
        transpose(HalfSig * EmbeddingMatrix(DS, dM, θ)) * n - EvalLogPriorGrad(LogPriorFn, θ)
    end
    dF1(θ::AbstractVector) = transpose(HalfSig * EmbeddingMatrix(DS, dM, θ)) *  sign.(HalfSig * (EmbeddingMap(DS, M, θ) - ydata(DS))) - EvalLogPrior(LogPriorFn, θ)
    InformationGeometry.minimize(F, (p == 1 ? dF1 : dFp), start, Domain; tol=tol, kwargs...)
end

"""
    LineSearch(Test::Function, start::Number=0.; tol::Real=8e-15, maxiter::Int=10000) -> Number
Finds real number `x` where the boolean-valued `Test(x::Number)` goes from `true` to `false`.
"""
function LineSearch(Test::Function, start::Number=0.; tol::Real=8e-15, maxiter::Int=10000)
    if ((suff(start) != BigFloat) && tol < 1e-15)
        println("LineSearch: start not BigFloat but tol=$tol. Promoting and continuing.")
        start = BigFloat(start)
    end
    if !Test(start)
        start += 1e-10
        println("LineSearch: Test(start) did not work, trying Test(start + 1e-10).")
        !Test(start) && throw(ArgumentError("LineSearch: Test not true for starting value."))
    end
    # For some weird reason, if the division by 4 is removed, the loop never terminates for BigFloat-valued "start"s - maybe the compiler erroneously tries to optimize the variable "stepsize" away or something?! (Julia version ≤ 1.6.0)
    stepsize = one(suff(start)) / 4.;       value = start
    for i in 1:maxiter
        if Test(value + stepsize) # inside
            value += stepsize
            # value - start > 2000. && throw("FindConfBoundary: Value larger than 2000.")
        else            #outside
            if stepsize < tol
                return value + stepsize
            end
            stepsize /= 5.
        end
    end
    throw("$maxiter iterations over. Value=$value, Stepsize=$stepsize")
end

function AltLineSearch(Test::Function, Domain::Tuple{T,T}=(0., 1e4), meth::Roots.AbstractUnivariateZeroMethod=Roots.AlefeldPotraShi(); tol::Real=1e-12) where T<:Real
    find_zero(Test, Domain, meth; xatol=tol, xrtol=tol)
end
function AltLineSearch(Test::Function, Domain::Tuple{T,T}, meth::Roots.AbstractUnivariateZeroMethod=Roots.AlefeldPotraShi(); tol::Real=convert(BigFloat,10^(-precision(BigFloat)/10))) where T<:BigFloat
    Res = find_zero(Test, (Float64(Domain[1]), Float64(Domain[2])), meth; xatol=1e-14, xrtol=1e-14)
    find_zero(Test, (BigFloat(Res-2e-14),BigFloat(Res+2e-14)), Bisection(); xatol=tol, xrtol=tol)
end


text:
# IO Functions
# ------------

"""
    unsafe_read(input::IO, output::Ptr{UInt8}, nbytes::Int)::Int

Copy at most `nbytes` from `input` into `output`.

This function is similar to `Base.unsafe_read` but is different in some points:
- It does not throw `EOFError` when it fails to read `nbytes` from `input`.
- It returns the number of bytes written to `output`.
- It does not block if there are buffered data in `input`.
"""
function unsafe_read(input::IO, output::Ptr{UInt8}, nbytes::Int)::Int
    p = output
    navail = bytesavailable(input)
    if navail == 0 && nbytes > 0 && !eof(input)
        b = read(input, UInt8)
        unsafe_store!(p, b)
        p += 1
        nbytes -= 1
        navail = bytesavailable(input)
    end
    n = min(navail, nbytes)
    Base.unsafe_read(input, p, n)
    p += n
    return Int(p - output)
end

ext: jl


text:
function ps_cost(ps_m::CanonicalModel,
                 variable::JuMP.Containers.DenseAxisArray{JV},
                 cost_component::Function,
                 sign::Int64) where {JV <: JuMP.AbstractVariableRef}

    store = Vector{Any}(undef,length(variable))

    for (ix, element) in enumerate(variable)
        store[ix] = cost_component(element)
    end

    gen_cost = sum(store)

    return sign*gen_cost

end

function ps_cost(ps_m::CanonicalModel,
                 variable::JuMP.Containers.DenseAxisArray{JV},
                 cost_component::Float64,
                 sign::Int64) where {JV <: JuMP.AbstractVariableRef}

    gen_cost = sum(variable)*cost_component

    return sign*gen_cost

end

function pwl_param_check(cost_::Array{Tuple{Float64, Float64}})
    flag = true;
    for i in 1:(length(cost_)-1)
        if i == 1
            (cost_[i][1]/cost_[i][2]) <= ((cost_[i+1][1] - cost_[i][1])/(cost_[i+1][2] - cost_[i][2])) ? nothing : flag = false;
        else
            ((cost_[i][1] - cost_[i-1][1])/(cost_[i][2] - cost_[i-1][2])) <= ((cost_[i+1][1] - cost_[i][1])/(cost_[i+1][2] - cost_[i][2])) ? nothing : flag = false;
        end
    end
    return flag
end

function pwl_gen_cost(ps_m::CanonicalModel,
                    variable::JV,
                    cost_component::Array{Tuple{Float64, Float64}}) where {JV <: JuMP.AbstractVariableRef}

    gen_cost = JuMP.GenericAffExpr{Float64, _variable_type(ps_m)}()
    pwl_param_check(cost_component) ? @warn("Data provide is not suitable for linear implementation of PWL cost, this will result in a INVALID SOLUTION") : nothing ;
    # TODO: implement a fallback to either Linear Cost function or SOS2 based PWL Cost function
    upperbound(i) = (i == 1 ? cost_component[i+1][2] : (cost_component[i+1][2] - cost_component[i][2]));
    pwlvars = JuMP.@variable(ps_m.JuMPmodel, [i = 1:(length(cost_component)-1)], base_name = "{$(variable)}_{pwl}", start = 0.0, lower_bound = 0.0, upper_bound = upperbound(i))

    for (ix, pwlvar) in enumerate(pwlvars)
        if ix == 1
            temp_gen_cost = cost_component[ix][1] * (pwlvar / cost_component[ix][2] ) ;
        else
            temp_gen_cost = (cost_component[ix][1] - cost_component[ix-1][1]) * (pwlvar/(cost_component[ix][2] - cost_component[ix-1][2]) );
        end
        gen_cost = gen_cost + temp_gen_cost
    end

    c = JuMP.@constraint(ps_m.JuMPmodel, variable == sum([pwlvar for (ix, pwlvar) in enumerate(pwlvars) if ix > 1]) )

    return gen_cost

end

function ps_cost(ps_m::CanonicalModel,
                 variable::JuMP.Containers.DenseAxisArray{JV},
                 cost_component::Array{Tuple{Float64, Float64}},
                 sign::Int64) where {JV <: JuMP.AbstractVariableRef}

    gen_cost = JuMP.GenericAffExpr{Float64, _variable_type(ps_m)}()

    for var in variable
        c = pwl_gen_cost(ps_m, var, cost_component)
        gen_cost += c
    end

    return sign*gen_cost

end

function add_to_cost(ps_m::CanonicalModel,
                     devices::Array{C,1},
                     var_name::Symbol,
                     cost_symbol::Symbol, sign::Int64 = 1) where {C <: PSY.Device}

   for d in devices
        cost_expression = ps_cost(ps_m, ps_m.variables[var_name][d.name,:], getfield(d.econ,cost_symbol), sign)
        ps_m.cost_function += cost_expression
    end

    return

end


text:
"""
Functionality for authenticating Genie users.
"""
module GenieAuthentication

import Genie, SearchLight

export authenticate, deauthenticate, is_authenticated, get_authentication, authenticated
export login, logout, with_authentication, without_authentication, @authenticated!

const USER_ID_KEY = :__auth_user_id

"""
Stores the user id on the session.
"""
function authenticate(user_id::Any, session::Genie.Sessions.Session) :: Genie.Sessions.Session
  Genie.Sessions.set!(session, USER_ID_KEY, user_id)
end
function authenticate(user_id::SearchLight.DbId, session::Genie.Sessions.Session)
  authenticate(Int(user_id.value), session)
end
function authenticate(user_id::Union{String,Symbol,Int,SearchLight.DbId}, params::Dict{Symbol,Any} = Genie.Requests.payload()) :: Genie.Sessions.Session
  authenticate(user_id, params[:SESSION])
end

"""
    deauthenticate(session) :: Sessions.Session
    deauthenticate(params::Dict{Symbol,Any}) :: Sessions.Session

Removes the user id from the session.
"""
function deauthenticate(session::Genie.Sessions.Session) :: Genie.Sessions.Session
  Genie.Sessions.unset!(session, USER_ID_KEY)
end
function deauthenticate(params::Dict = Genie.Requests.payload()) :: Genie.Sessions.Session
  deauthenticate(params[:SESSION])
end

"""
    is_authenticated(session) :: Bool
    is_authenticated(params::Dict{Symbol,Any}) :: Bool

Returns `true` if a user id is stored on the session.
"""
function is_authenticated(session::Union{Genie.Sessions.Session,Nothing}) :: Bool
  Genie.Sessions.isset(session, USER_ID_KEY)
end
function is_authenticated(params::Dict = Genie.Requests.payload()) :: Bool
  is_authenticated(params[:SESSION])
end

const authenticated = is_authenticated

"""
    @authenticate!(exception::E = ExceptionalResponse(Genie.Renderer.redirect(:show_login)))

If the current request is not authenticated it throws an ExceptionalResponse exception.
"""
macro authenticated!(exception = Genie.Exceptions.ExceptionalResponse(Genie.Renderer.redirect(:show_login)))
  :(GenieAuthentication.authenticated() || throw($exception))
end

"""
    get_authentication(session) :: Union{Nothing,Any}
    get_authentication(params::Dict{Symbol,Any}) :: Union{Nothing,Any}

Returns the user id stored on the session, if available.
"""
function get_authentication(session::Genie.Sessions.Session) :: Union{Nothing,Any}
  Genie.Sessions.get(session, USER_ID_KEY)
end
function get_authentication(params::Dict = Genie.Requests.payload()) :: Union{Nothing,Any}
  get_authentication(params[:SESSION])
end

const authentication = get_authentication

"""
    login(user, session) :: Union{Nothing,Genie.Sessions.Session}
    login(user, params::Dict{Symbol,Any}) :: Union{Nothing,Genie.Sessions.Session}

Persists on session the id of the user object and returns the session.
"""
function login(user::M, session::Genie.Sessions.Session)::Union{Nothing,Genie.Sessions.Session} where {M<:SearchLight.AbstractModel}
  authenticate(getfield(user, Symbol(pk(user))), session)
end
function login(user::M, params::Dict = Genie.Requests.payload())::Union{Nothing,Genie.Sessions.Session} where {M<:SearchLight.AbstractModel}
  login(user, params[:SESSION])
end

"""
    logout(session) :: Sessions.Session
    logout(params::Dict{Symbol,Any}) :: Sessions.Session

Deletes the id of the user object from the session, effectively logging the user off.
"""
function logout(session::Genie.Sessions.Session) :: Genie.Sessions.Session
  deauthenticate(session)
end
function logout(params::Dict = Genie.Requests.payload()) :: Genie.Sessions.Session
  logout(params[:SESSION])
end

"""
    with_authentication(f::Function, fallback::Function, session)
    with_authentication(f::Function, fallback::Function, params::Dict{Symbol,Any})

Invokes `f` only if a user is currently authenticated on the session, `fallback` is invoked otherwise.
"""
function with_authentication(f::Function, fallback::Function, session::Union{Genie.Sessions.Session,Nothing})
  if ! is_authenticated(session)
    fallback()
  else
    f()
  end
end
function with_authentication(f::Function, fallback::Function, params::Dict = Genie.Requests.payload())
  with_authentication(f, fallback, params[:SESSION])
end

"""
    without_authentication(f::Function, session)
    without_authentication(f::Function, params::Dict{Symbol,Any})

Invokes `f` if there is no user authenticated on the current session.
"""
function without_authentication(f::Function, session::Genie.Sessions.Session)
  ! is_authenticated(session) && f()
end
function without_authentication(f::Function, params::Dict = Genie.Requests.payload())
  without_authentication(f, params[:SESSION])
end

"""
    install(dest::String; force = false, debug = false) :: Nothing

Copies the plugin's files into the host Genie application.
"""
function install(dest::String; force = false, debug = false) :: Nothing
  src = abspath(normpath(joinpath(pathof(@__MODULE__) |> dirname, "..", Genie.Plugins.FILES_FOLDER)))

  debug && @info "Preparing to install from $src into $dest"
  debug && @info "Found these to install $(readdir(src))"

  for f in readdir(src)
    debug && @info "Processing $(joinpath(src, f))"
    debug && @info "$(isdir(joinpath(src, f)))"

    isdir(joinpath(src, f)) || continue

    debug && "Installing from $(joinpath(src, f))"

    Genie.Plugins.install(joinpath(src, f), dest, force = force)
  end

  nothing
end

end

max_stars_repo_name: genieframework/GenieAuthentication

max_issues_repo_name: genieframework/GenieAuthentication

max_forks_repo_name: genieframework/GenieAuthentication


text:
export PermutationRepresentationDrawing

"""
`PermutationRepresentationDrawing(G)` draws a picture of a permtuation
representation of the graph `G` (if one exists). An optional second
argument (of type `Int`) determines the point size of the labels.
"""
function PermutationRepresentationDrawing(G::SimpleGraphs.SimpleGraph, font_size = 8)
    f, g = PermutationRepresentation(G)
    PermutationRepresentationDrawing(f, g, font_size)
end

"""
`PermutationRepresentationDrawing(f,g)` creates a permutation graph
drawing from a pair of dictionaries mapping vertices to integers.
This is used by `PermutationRepresentationDrawing(G)`.
"""
function PermutationRepresentationDrawing(
    f::Dict{T,Int},
    g::Dict{T,Int},
    font_size = 8,
) where {T}
    vertices = collect(keys(f))
    n = length(vertices)
    MARKER_SIZE = 10

    yy = sqrt(length(f))

    plot()
    for v in vertices
        x = [f[v], g[v]]
        y = [yy, 0]
        plot!(x, y, color = "black")
    end
    for v in vertices
        plot!(
            [f[v]],
            [yy],
            markerstrokecolor = "black",
            markercolor = "white",
            marker = MARKER_SIZE,
            markerstrokewidth = 2,
        )
        plot!(
            [g[v]],
            [0],
            markerstrokecolor = "black",
            markercolor = "white",
            marker = MARKER_SIZE,
            markerstrokewidth = 2,
        )
    end

    for v in vertices
        if font_size > 0
            annotate!((f[v], yy, string(v), font_size))
            annotate!((g[v], 0, string(v), font_size))
        end
    end
    finish()
end


text:
using Documenter
using PkgBenchmark
using Transducers
using TransducersBenchmarksReports
using TransducersBenchmarksReports: mayberun, benchpath

generatedroot = joinpath(@__DIR__, "src", "generated")
mkpath(generatedroot)

result = mayberun()
PkgBenchmark.export_markdown(joinpath(generatedroot, "results.md"), result)

makedocs(;
    modules=[TransducersBenchmarksReports],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Results" => "generated/results.md",
        "GEMM" => "gemm.md",
    ],
    repo="https://github.com/tkf/TransducersBenchmarksReports.jl/blob/{commit}{path}#L{line}",
    sitename="TransducersBenchmarksReports.jl",
    authors="Takafumi Arakaki <aka.tkf@gmail.com>",
)

deploydocs(;
    repo="github.com/tkf/TransducersBenchmarksReports.jl",
)


text:
abstract type AbstractBPNG <: AbstractController end

mutable struct BPNG <: AbstractBPNG
    N::Float64
    dim::Int
    σ_M_lim::Float64
    v̂_f_d::Vector
    Bias::Function
    s_Bias # parameters of bias command function
end

function BPNG_cmd(s_guidance::BPNG)
    @unpack N, dim, Bias = s_guidance
    return function (x, params, t)
        (p_M, v_M, p_T, v_T) = (x.pursuer.p, x.pursuer.v, x.evador.p, x.evador.v)
        if dim == 2
            (p_M, v_M, p_T, v_T) = vcat.((p_M, v_M, p_T, v_T), 0)
        end

        ω_r      = cross(p_T - p_M, v_T - v_M) / dot(p_T - p_M, p_T - p_M)
        a_M_PPNG = N * cross(ω_r, v_M)
        a_M_bias = Bias(s_guidance, x, t)[1]

        a_M      = a_M_PPNG + a_M_bias

        # if norm(a_M) >= A_M_max
        #     a_M = normalize(a_M) * clamp(norm(a_M), 0, A_M_max)
        # end
        if dim == 2
            a_M = a_M[1:2]
        end

        return a_M
    end
end

function Bias_IACG_StationaryTarget(s_guidance::BPNG, x, t)
    @unpack N, dim, σ_M_lim, v̂_f_d, s_Bias = s_guidance
    @unpack α, δ, n, r_ref, k, m, k̂_d, i_Ω_μ, i_σ_M_lim = s_Bias
    (p_M, v_M, p_T, v_T) = (x.pursuer.p, x.pursuer.v, x.evador.p, x.evador.v)
    if dim == 2
        (p_M, v_M, p_T, v_T) = vcat.((p_M, v_M, p_T, v_T), 0)
    end
    # @assert norm(v_T) < 0.1

    r        = norm(p_T - p_M)
    r̂        = normalize(p_T - p_M)
    ṙ        = dot(r̂, v_T - v_M)
    v̂_M      = normalize(v_M)
    σ_M      = atan(norm(cross(p_T - p_M, v_M)), dot(p_T - p_M, v_M))
    # σ_M      = acos( clamp( dot(r̂, v̂_M), -1, 1 ) )
    # atan-based implementation is more robustly accurate than using acos for computation
    k̂        = normalize(cross(p_T - p_M, v_M))

    v̂_f_pred = v̂_M * cos(N / (N - 1) * σ_M) - cross(k̂, v̂_M) * sin(N / (N - 1) * σ_M)
    # alternative: v̂_f_pred =   r̂ * cos(1 / (N - 1) * σ_M) -   cross(k̂, r̂) * sin(1 / (N - 1) * σ_M)
    e_v̂_f    = atan(norm(cross(v̂_f_pred, v̂_f_d)), dot(v̂_f_pred, v̂_f_d))

    # Design parameters for vec_BPNG
    k_e      = (n + 1) / r

    μ = 0.0
    if i_Ω_μ == 1
        # desired manoeuvre plane alignment
        k_μ  = (m + 1) / r
        μ    = atan(dot(cross(cross(v̂_f_pred, v̂_f_d), k̂_d), v̂_f_d), dot(cross(v̂_f_pred, v̂_f_d), k̂_d))
        if abs(μ) > pi / 2
            μ    = atan(dot(cross(cross(v̂_f_pred, v̂_f_d), -k̂_d), v̂_f_d), dot(cross(v̂_f_pred, v̂_f_d), -k̂_d))
        end
        Ω_μ  = k_μ * abs(ṙ) * μ

    elseif i_Ω_μ == 2
        # evasion: monotonic decrease
        Ω_μ  = δ * ( r / r_ref )^m
    elseif i_Ω_μ == 3
        # evasion: sinusoidal wave
        Ω_μ  = δ * sin(( r / r_ref )^m)
    elseif i_Ω_μ == 4
        # evasion: square wave
        Ω_μ  = δ * ifelse(mod(( r / r_ref )^m, 1) < θ, 1.0, -1.0)
    elseif i_Ω_μ == 5
        # evasion: sawtooth wave
        Ω_μ  = δ * rem(( r / r_ref )^m, 1.0, RoundNearest) * 2
    elseif i_Ω_μ == 6
        # evasion: triangular wave
        modr = mod(( r / r_ref )^m + 1/4, 1.0)
        Ω_μ  = δ * ifelse(modr < 1/2, 4*modr - 1, -4*modr + 3)
    elseif i_Ω_μ == 7
        # evasion: uniform random number (s_rand: a fixed sequence of m uniform randoms in (-1,1) generated externally)
        # Ω_μ  = δ * s_rand[1 + mod(floor(Int, t / r_ref), m)]
    else
        Ω_μ  = 0.0
    end

    if i_σ_M_lim == 1
        # lead angle limiter
        k_σ  = 2 / pi * acos( clamp((σ_M / σ_M_lim)^k , -1, 1) )
    else
        k_σ  = 1
    end

    ω_f      = k_e * abs(ṙ) * e_v̂_f^α * normalize(cross(v̂_f_pred, v̂_f_d)) + Ω_μ * v̂_f_d
    ω_bias   = k_σ * ( -(N - 1) * dot(ω_f, k̂) * k̂  + sin(σ_M) * dot(ω_f, r̂ + cot(1 / (N - 1) * σ_M) * cross(k̂, r̂)) * cross(v̂_M, k̂) )

    a_M_bias = cross(ω_bias, v_M)

    return a_M_bias, e_v̂_f, ω_bias, μ, Ω_μ
end

function Bias_Quaternion_IACG(s_guidance::BPNG, x, t)
    @unpack N, dim, σ_M_lim, v̂_f_d, s_Bias = s_guidance
    @unpack α, n = s_Bias
    (p_M, v_M, p_T, v_T) = (x.pursuer.p, x.pursuer.v, x.evador.p, x.evador.v)
    if dim == 2
        (p_M, v_M, p_T, v_T) = vcat.((p_M, v_M, p_T, v_T), 0)
    end
    # @assert norm(v_T) < 0.1

    r        = norm(p_T - p_M)
    r̂        = normalize(p_T - p_M)
    ṙ        = dot(r̂, v_T - v_M)
    v̂_M      = normalize(v_M)
    # σ_M      = atan(norm(cross(p_T - p_M, v_M)), dot(p_T - p_M, v_M))
    # σ_M      = acos( clamp( dot(r̂, v̂_M), -1, 1 ) )
    # atan-based implementation is more robustly accurate than using acos for computation

    e_v̂_f_instant    = atan(norm(cross(v̂_M, v̂_f_d)), dot(v̂_M, v̂_f_d))

    # Design parameters for Quaternion_IACG
    N_f      = (n + 1.0) * (n + 2.0)
    t_go     = r / abs(ṙ)
    ω_bias   = N_f / t_go * e_v̂_f_instant^α * normalize(cross(v̂_f_d, v̂_M))

    a_M_bias = cross(ω_bias, v_M)

    μ   = 0.0
    Ω_μ = 0.0

    return a_M_bias, e_v̂_f_instant, ω_bias, μ, Ω_μ
end

#--------------------------------------Legacy codes--------------------------------------
function Bias_IACG_StationaryTarget_2D(s_guidance::BPNG, x, t)
    @unpack N, dim, σ_M_lim, v̂_f_d, s_Bias = s_guidance
    @unpack α, δ, n, r_ref, k, m = s_Bias
    # @assert dim == 2
    (p_M, v_M, p_T, v_T) = (x.pursuer.p, x.pursuer.v, x.evador.p, x.evador.v)
    r       = norm(p_T - p_M)
    λ       = atan(p_T[2] - p_M[2], p_T[1] - p_M[1])

    V_M     = norm(v_M)
    γ_M     = atan(v_M[2], v_M[1])
    σ_M     = γ_M - λ
    η       = sin(σ_M)
    η_lim   = sin(σ_M_lim)

    # V_T     = norm(v_T)
    # γ_T     = atan(v_T[2], v_T[1])
    # σ_T     = γ_T - λ

    γ_f_d   = atan(v̂_f_d[2], v̂_f_d[1])

    # λ̇       = (V_T*sin(σ_T) - V_M*sin(σ_M)) / r
    e_γ_f   = γ_M - N / (N - 1) * σ_M - γ_f_d
    if α > 0.99
        e_γ_f_fbk = e_γ_f
    else
        e_γ_f_fbk = sign(e_γ_f) * abs(e_γ_f)^α
    end

    K_r     = (N - 1 + δ) * (1 + r / r_ref)^n
    K_eta   = 1 + k * (1 - (abs(η / η_lim))^m)
    u_aug = - K_r * K_eta / r * e_γ_f_fbk

    A_M_bias = - u_aug * V_M^2 * cos(σ_M)
    a_M_bias = A_M_bias * [-sin(γ_M); cos(γ_M); 0]

    return a_M_bias

    # if dim == 2  # identical to the codes above, thus can be deleted
    #     Design parameters for LCG_e_0
    #     K_e@vec_BPNG = K_r@LCG_e_0 / (N - 1) / r
    #     kwargs: δ = 0.01, n = 0, r_ref = Inf, k = 0, m = 0
    #     K_e = (N - 1 + δ) / (N - 1) / r * (1 + r / r_ref)^n
    #     Ω_μ = 0
    #     K_σ = 1 + k * (1 - (abs(sin(σ_M) / sin(σ_M_lim)))^m)
    #
    #     ê_z = [0; 0; 1]
    #     if cross(p_T - p_M, v_M)[3] >= 0
    #         k̂ = ê_z
    #     else
    #         k̂ = -ê_z
    #     end
    #     v̂_f_pred = v̂_M * cos(N / (N - 1) * σ_M) - cross(k̂, v̂_M) * sin(N / (N - 1) * σ_M)
    #     e_v̂_f_signed = atan(dot(v̂_f_pred, cross(ê_z, v̂_f_d)), dot(v̂_f_pred, v̂_f_d))
    #     # alternative: e_v̂_f_signed = atan(v̂_f_pred[2], v̂_f_pred[1]) - atan(v̂_f_d[2], v̂_f_d[1])
    #     ω_f     = K_e(r, N, δ, n, r_ref) * ṙ  * sign(e_v̂_f_signed) * abs(e_v̂_f_signed)^α * ê_z
    #     ω_bias  = K_σ(σ_M, σ_M_lim, k, m) * ( -(N - 1) * ω_f )
    # end
end

function Bias_zero(s_guidance::BPNG, x, t)
    return zeros(3)
end


text:
using jInv.Mesh;
using FactoredEikonalFastMarching;
using PyPlot
close("all")

include("runAccuracyExperiment.jl");
include("getAnalyticalMediums.jl");
include("getWorkUnit.jl");

function runPaperExperiments(numOf2DRefinements::Int=4, numOf3DRefinements::Int=4)

	println("******************** 2D experiments ***************************");

	I = [4,8];
	n = zeros(Int64,2);

	h0 = [0.05,0.05];
	WU = zeros(numOf2DRefinements)
	for k = 1:numOf2DRefinements
		h = (0.5^k)*h0;
		n[1] = round(Int64,(I[1]/h[1])+1);
		n[2] = round(Int64,(I[2]/h[2])+1);
		WU[k] = GetWorkunitForLoop(n,h);
	end

	println("********************  Gradient Slowness  ***************************");
	println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")

	h0 = [0.05,0.05];
	for k=1:numOf2DRefinements
		h = (0.5^k)*h0;
		n[1] = round(Int64,(I[1]/h[1])+1);
		n[2] = round(Int64,(I[2]/h[2])+1);
		tc="2d-testcase1";
		println("*******************************************************************************")
		println("2d-Testcase 1")
		println("*******************************************************************************")
		(kappaSquared,src,T_exact) = getAnalyticalConstGrad2D(n,h);
		runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k],tc);
	end

	println("********************  Gradient Velocity  ***************************");

	println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")
	h0 = [0.05,0.05];
	for k=1:numOf2DRefinements
		h = (0.5^k)*h0;
		n[1] = round(Int64,(I[1]/h[1])+1);
		n[2] = round(Int64,(I[2]/h[2])+1);
		tc="2d-testcase2";
		println("*******************************************************************************")
		println("*******************************************************************************")
		println("2d-Testcase 2")
		println("*******************************************************************************")
		(kappaSquared,src,T_exact) = getAnalyticalConstGradInv2D(n,h);
		runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k],tc);
	end

	println("********************  Gaussian factor  *****************************");

	println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")
	h0 = [0.05,0.05];
	for k=1:numOf2DRefinements
		h = (0.5^k)*h0;
		n[1] = round(Int64,(I[1]/h[1])+1);
		n[2] = round(Int64,(I[2]/h[2])+1);
		tc="2d-testcase3";
		println("*******************************************************************************")
		println("*******************************************************************************")
		println("2d-Testcase 3")
		println("*******************************************************************************")
		(kappaSquared,src,T_exact) = getSmoothGaussianMedium(n,h);
		runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k],tc);
	end

	println("*******************************************************************************")#
	println("*******************************************************************************")
	println("************************* 3-D EXPERIMENTS *************************************")
	println("*******************************************************************************")
	println("*******************************************************************************")

	##########################################################################################
	##########################################################################################
	# 3D experiments:
	##########################################################################################

	I = [1.6,1.6,0.8];
	n = zeros(Int64,3);

	# numOfRefinements = 3;

	h0 = [0.1,0.1,0.1];
	WU = zeros(numOf3DRefinements)
	for k = 1:numOf3DRefinements
		h = (0.5^k)*h0;
		n[1] = round(Int64,(I[1]/h[1])+1);
		n[2] = round(Int64,(I[2]/h[2])+1);
		n[3] = round(Int64,(I[3]/h[3])+1);
		WU[k] = GetWorkunitForLoop(n,h);
	end

	println("********************  Gradient Slowness  ***************************");
	println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")

	for k=1:numOf3DRefinements
		h = (0.5^k)*h0;
		n[1] = round(Int64,(I[1]/h[1])+1);
		n[2] = round(Int64,(I[2]/h[2])+1);
		n[3] = round(Int64,(I[3]/h[3])+1);
		tc="3d-testcase1";
		println("*******************************************************************************")
		println("*******************************************************************************")
		println("3d-Testcase 1")
		println("*******************************************************************************")
		(kappaSquared,src,T_exact) = getAnalyticalConstGrad3D(n,h);
		runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k],tc);
	end

	println("********************  Gradient Velocity  ***************************");
	println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")
	for k=1:numOf3DRefinements
		h = (0.5^k)*h0;
		n[1] = round(Int64,(I[1]/h[1])+1);
		n[2] = round(Int64,(I[2]/h[2])+1);
		n[3] = round(Int64,(I[3]/h[3])+1);
		tc="3d-testcase2";
		println("*******************************************************************************")
		println("*******************************************************************************")
		println("3d-Testcase 2")
		println("*******************************************************************************")
		(kappaSquared,src,T_exact) = getAnalyticalConstGradInv3D(n,h);
		runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k],tc);
	end

	println("********************  Gaussian factor  *****************************");
	println("h\t\tn\t\terror(1st)[max,l2]\ttime(1st)\terror(2nd)[max,l2]\ttime(2nd)")
	for k=1:numOf3DRefinements
		h = (0.5^k)*h0;
		n[1] = round(Int64,(I[1]/h[1])+1);
		n[2] = round(Int64,(I[2]/h[2])+1);
		n[3] = round(Int64,(I[3]/h[3])+1);
		tc="3d-testcase3";
		println("*******************************************************************************")
		println("*******************************************************************************")
		println("3d-Testcase 3")
		println("*******************************************************************************")
		(kappaSquared,src,T_exact) = getSmoothGaussianMedium(n,h);
		runExperimentAndWriteResults(kappaSquared,h,src,n,T_exact,WU[k],tc);
	end
	return;
end
############################################################################################

try runPaperExperiments(parse(Int,ARGS[1]), parse(Int,ARGS[2]));
catch e
	println(e)
	exit(1)
end

exit(0)


text:
# Starting point example for converting data types in Learning Julia

# declare some basic data types
x = 10
y = 20.0
z = "40"
w = "30.0"

# Convert an integer to a character and vice versa
c = Char(74)
println(c)

i = Int('J')
println(i)

# convert between types
flt1 = Float32(x)
int1 = Int16(y)
println(flt1)
println(int1)

# parse values from a string
a = parse(Int, z)
b = parse(Float64, w)
c = string(y)
println(a)
println(b)
println(y)

# attempting to convert a data type that won't fit is an error
x = Char(281)
println(x)

# xc = UInt8(x) -- unable to cast into Unsigned Int8
# println(xc)

max_stars_repo_name: scaotravis/Julia-Course

max_issues_repo_name: scaotravis/Julia-Course

max_forks_repo_name: scaotravis/Julia-Course