text:
# Deprecations added in 1.5.0, March 2021

function File(fmt::Type{DataFormat{sym}}, filename) where {sym}
    Base.depwarn("`File(format\"$sym\", filename)` is deprecated, please use `File{format\"$sym\"}(filename)` instead.", :File)
    return File{fmt}(filename)
end
function Stream(fmt::Type{DataFormat{sym}}, args...) where {sym}
    Base.depwarn("`Stream(format\"$sym\", filename)` is deprecated, please use `Stream{format\"$sym\"}(filename)` instead.", :Stream)
    return Stream{fmt}(args...)
end

# These aren't used here, but old versions of ImageIO expect them

function _findmod(f::Symbol)
    Base.depwarn("_findmod is deprecated and will be removed. Use `Base.require(::Base.PkgId)` instead.", :_findmod)
    for (u,v) in Base.loaded_modules
        (Symbol(v) == f) && return u
    end
    nothing
end
function topimport(modname)
    Base.depwarn("topimport is deprecated and will be removed. Use `Base.require(::Base.PkgId)` instead.", :topimport)
    @eval Base.__toplevel__  import $modname
    u = _findmod(modname)
    @eval $modname = Base.loaded_modules[$u]
end

# Legacy add_loader/add_saver
for add_ in (:add_loader, :add_saver)
    @eval begin
        function $add_(fmt, pkg)
            # TODO: delete this method in FileIO v2
            sym = isa(fmt, Symbol) ? fmt : formatname(fmt)::Symbol
            Base.depwarn(string($add_) * "(fmt, pkg::$(typeof(pkg))) is deprecated, supply `pkg` as a Module or `name=>uuid`", Symbol($add_))
            pkg === :MimeWriter && return $add_(sym, MimeWriter)
            # Try to look it up in the caller's environment
            pkgname = string(pkg)
            id = Base.identify_package(pkgname)
            if id === nothing
                # See if it's in Main
                pkgsym = Symbol(pkg)
                if isdefined(Main, pkgsym)
                    id = getfield(Main, pkgsym)
                    if !isa(id, Module)
                        id = nothing
                    end
                end
                if id === nothing
                    # Look it up in the registries. The tricky part here is supporting different Julia versions
                    ctx = Pkg.API.Context()
                    uuids = UUID[]
                    @static if Base.VERSION >= v"1.2"
                        if hasfield(typeof(ctx), :registries)
                            for reg in ctx.registries
                                append!(uuids, Pkg.Registry.uuids_from_name(reg, pkgname))
                            end
                        else
                            ctx = Pkg.API.Context!(ctx)
                            if isdefined(Pkg.Types, :find_registered!) && hasmethod(Pkg.Types.find_registered!, (typeof(ctx.env), Vector{String}))
                                Pkg.Types.find_registered!(ctx.env, [pkgname])
                            elseif isdefined(Pkg.Types, :find_registered!) && hasmethod(Pkg.Types.find_registered!, (typeof(ctx), Vector{String}))
                                Pkg.Types.find_registered!(ctx, [pkgname])
                            end
                            append!(uuids, get(ctx.env.uuids, pkgname, UUID[]))
                        end
                    else
                        Pkg.Types.find_registered!(ctx.env)
                        append!(uuids, get(ctx.env.uuids, pkgname, UUID[]))
                    end
                    isempty(uuids) && throw(ArgumentError("no UUID found for $pkg"))
                    length(uuids) == 1 || throw(ArgumentError("multiple UUIDs found for $pkg"))
                    id = PkgId(uuids[1], pkgname)
                end
            end
            $add_(sym, id)
        end
    end
end


text:
struct ComteButucea <: LinearEstimator
    Q::BinnedCalibrator
    h::Float64
    m::Int64
    target::InferenceTarget
end

function (c::ComteButucea)(x)
    (c.Q)(x)
end

pretty_label(cb::ComteButucea)  = "Butucea-Comte"

inference_target(cb::ComteButucea) = cb.target

default_bandwidth(::Type{ComteButucea}, target::MarginalDensityTarget, n) = sqrt(log(n))
default_bandwidth(::Type{ComteButucea}, target::LFSRNumerator, n) = sqrt(log(n)*2/3)

default_bandwidth(::Type{ComteButucea}, n) = sqrt(log(n))
default_bandwidth(::Type{ComteButucea}, target, n) = default_bandwidth(ComteButucea, n)

function comte_butucea(z, h, target::InferenceTarget)
   f(t) = real(exp(im*t*z)*cf(target,-t)/cf(Normal(),t))
   res = hquadrature(f, -h, +h)[1]
   res/2/π
end

function ComteButucea(target::InferenceTarget, m::Int64, marginal_grid;
    h= default_bandwidth(ComteButucea, target, m))
    Qc = BinnedCalibrator(marginal_grid, comte_butucea.(marginal_grid, h, target))
    ComteButucea(Qc, h, m, target)
end

function estimate(Xs,cb::ComteButucea)
    mean(cb.Q.(Xs))
end

function estimate(Xs, ::Type{ComteButucea}, target::InferenceTarget, marginal_grid)
    cb = ComteButucea(target, length(Xs), marginal_grid)
    estimate(Xs, cb)
end
# probably should learn how to do the above properly via FFT?

function estimate(Xs, target::InferenceTarget, ::Type{ComteButucea}; K_grid=10_000,
    marginal_grid = collect(range(minimum(Xs), stop=maximum(Xs), length=K_grid)), args...)
    n = length(Xs)
    Qc = BinnedCalibrator(marginal_grid, comte_butucea.(marginal_grid, target, n; args...))
    mean(Qc.(Xs))
end

ext: jl


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	load_period_map(setup::Dict,path::AbstractString, inputs::Dict)

Function for reading input parameters related to mapping of representative time periods to full chronological time series
"""
function load_period_map(setup::Dict,path::AbstractString, inputs::Dict)
	data_directory = joinpath(path, setup["TimeDomainReductionFolder"])
	if setup["TimeDomainReduction"] == 1  && isfile(joinpath(data_directory,"Period_map.csv"))  # Use Time Domain Reduced data for GenX
		inputs["Period_Map"] = DataFrame(CSV.File(joinpath(data_directory,"Period_map.csv"), header=true), copycols=true)
	else
		inputs["Period_Map"] = DataFrame(CSV.File(joinpath(path,"Period_map.csv"), header=true), copycols=true)
	end

	println("Period_map.csv Successfully Read!")

	return inputs
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
include("structs.jl")
include("physics.jl")
include("computational.jl")

Broadcast.broadcastable(eos::TaitEoS) = Ref(eos)
Broadcast.BroadcastStyle(::Type{<:TaitEoS}) =
                        Broadcast.DefaultArrayStyle{0}()

Broadcast.broadcastable(eos::IdealGasEoS) = Ref(eos)
Broadcast.BroadcastStyle(::Type{<:IdealGasEoS}) =
                        Broadcast.DefaultArrayStyle{0}()

# -------------------- Density --------------------------- #

## The tait and ideal gas equations helps us determine
## ρ̂ᵢ= mᵢ/ sᵢV , but we need ρᵢ = mᵢ / V because ρᵢ is
## used in continuity and Darcy's equations. The first
## two functions return ρ̂ᵢ, while the third one returns ρᵢ

function density(x::TaitEoS, P::AbstractFloat)
    if x.B + P < 0
        println("Less than 0: ", P)
    end
    return x.ρ₀ / (1 - x.C * log10((x.B + P)
                                   / (x.B + x.P₀)))
end

function density(x::IdealGasEoS, P::AbstractFloat)
    return P * x.M / R / x.T
end

function density!(sys::System)
    @. @views sys.ρ[:, :, 1] = density(ideal_gas, sys.P) *
                                        sys.s[:, :, 1]
    @. @views sys.ρ[:, :, 2] = density(tait_C₅H₁₂, sys.P) *
                                        sys.s[:, :, 2]
end
# -------------------------------------------------------- #

# -------------------- Pressure -------------------------- #

function find_pressure(; root_finder::Function, ρ₁::T, ρ₂::T,
        tait::TaitEoS, igas::IdealGasEoS,
        p::Parameters, niter::Integer=50,
        eps_x::T=1e-6)::Tuple{T, T} where T<:AbstractFloat

    function f(P::T)::T where T<:AbstractFloat
        ρ̂₁ = density(igas, P)
        s  = ρ₁ / ρ̂₁
        ρ̂₂ = ρ₂ / (1 - s)
        return (ρ̂₂ - tait.ρ₀) / ρ̂₂ - tait.C *
            log10((tait.B + P) / (tait.B + tait.P₀))
    end

    function fder(P::T)::T where T<:AbstractFloat
        ρ̂₁ = density(igas, P)
        s  = ρ₁ / ρ̂₁
        ρ̂₂ = ρ₂ / (1 - s)
        return -tait.ρ₀ / ρ̂₂^2  * R * igas.T * igas.M *
            ρ₁ * ρ₂ / (igas.M * P - R * igas.T * ρ₁)^2 -
            tait.C / (tait.B + P) / log(10)
    end

    let
        P  = root_finder(; f=f, fder=fder, eps_x=eps_x)
        ρ̂₁ = density(igas, P)
        s  = ρ₁ / ρ̂₁

        return P, s
    end
end

function find_pressure!(sys::System, p::Parameters)
    for index in CartesianIndices(sys.P)
        i, j = Tuple(index)
        ## This gives seg fault, but is type stable
        # (sys.P[i, j], sys.s[i, j])::Tuple{T, T} =
        #             find_pressure(p=p,
        #             root_finder=newton_raphson,
        #             ρ₁=sys.ρ[i, j, 1],
        #             ρ₂=sys.ρ[i, j, 2],
        #             tait=tait_C₅H₁₂,
        #             igas=ideal_gas) where T<:AbstractFloat
        #

        ## Type unstable
        sys.P[i, j], sys.s[i, j, 1] =
                   find_pressure(p=p,
                   root_finder=newton_raphson,
                   ρ₁=sys.ρ[i, j, 1],
                   ρ₂=sys.ρ[i, j, 2],
                   tait=tait_C₅H₁₂,
                   igas=ideal_gas)

        sys.s[i, j, 2] = 1 - sys.s[i, j, 1]
    end
end

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

## -------------------- Saturation ---------------------- ##
##
## We derive the saturation from the molar composition
## ψ = ν₁ / ν₂ and the densities ρ̂₁ and ρ̂₂.

function saturation(ρ̂₁::T, ρ̂₂::T, ψ::T,
        M₁::T, M₂::T) where T<:AbstractFloat
    return (ρ̂₁ / ρ̂₂ * M₂ * (1 - ψ) / M₁ / ψ + 1)^(-1)
end

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

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration


text:
using NLSolvers
using Test
using StaticArrays
#using Optim
#using LineSearches
using Printf
using LinearAlgebra: norm, I

include("optimize/runtests.jl")
include("noalloc_static.jl")
include("scalar/runtests.jl")
include("specialtypes/runtests.jl")
include("newton.jl")


text:
#=
Provides a type called LQ for solving linear quadratic control
problems.

@author : Spencer Lyon <spencer.lyon@nyu.edu>
@author : Zac Cranko <zaccranko@gmail.com>

@date : 2014-07-05

References
----------

https://lectures.quantecon.org/jl/lqcontrol.html

=#

@doc doc"""
Linear quadratic optimal control of either infinite or finite horizon

The infinite horizon problem can be written

```math
\min \mathbb{E} \sum_{t=0}^{\infty} \beta^t r(x_t, u_t)
```

with

```math
r(x_t, u_t) := x_t' R x_t + u_t' Q u_t + 2 u_t' N x_t
```

The finite horizon form is

```math
\min \mathbb{E} \sum_{t=0}^{T-1} \beta^t r(x_t, u_t) + \beta^T x_T' R_f x_T
```

Both are minimized subject to the law of motion

```math
x_{t+1} = A x_t + B u_t + C w_{t+1}
```

Here ``x`` is n x 1, ``u`` is k x 1, ``w`` is j x 1 and the matrices are conformable for
these dimensions.  The sequence ``{w_t}`` is assumed to be white noise, with zero
mean and ``\mathbb{E} w_t w_t' = I``, the j x j identity.

For this model, the time ``t`` value (i.e., cost-to-go) function ``V_t`` takes the form

```math
x' P_T x + d_T
```

and the optimal policy is of the form ``u_T = -F_T x_T``.  In the infinite horizon
case, ``V, P, d`` and ``F`` are all stationary.

##### Fields

- `Q::ScalarOrArray` : `k x k` payoff coefficient for control variable u. Must be
  symmetric and nonnegative definite
- `R::ScalarOrArray` : `n x n` payoff coefficient matrix for state variable x.
  Must be symmetric and nonnegative definite
- `A::ScalarOrArray` : `n x n` coefficient on state in state transition
- `B::ScalarOrArray` : `n x k` coefficient on control in state transition
- `C::ScalarOrArray` : `n x j` coefficient on random shock in state transition
- `N::ScalarOrArray` : `k x n` cross product in payoff equation
- `bet::Real` : Discount factor in `[0, 1]`
- `capT::Union{Int, Void}` : Terminal period in finite horizon problem
- `rf::ScalarOrArray` : `n x n` terminal payoff in finite horizon problem. Must be symmetric and nonnegative definite
- `P::ScalarOrArray` : `n x n` matrix in value function representation ``V(x) = x'Px + d``
- `d::Real` : Constant in value function representation
- `F::ScalarOrArray` : Policy rule that specifies optimal control in each period

"""
mutable struct LQ
    Q::ScalarOrArray
    R::ScalarOrArray
    A::ScalarOrArray
    B::ScalarOrArray
    C::ScalarOrArray
    N::ScalarOrArray
    bet::Real
    capT::Union{Int, Nothing} # terminal period
    rf::ScalarOrArray
    P::ScalarOrArray
    d::Real
    F::ScalarOrArray # policy rule
end

"""
Main constructor for LQ type

Specifies default argumets for all fields not part of the payoff function or
transition equation.

##### Arguments

- `Q::ScalarOrArray` : `k x k` payoff coefficient for control variable u. Must be
  symmetric and nonnegative definite
- `R::ScalarOrArray` : `n x n` payoff coefficient matrix for state variable x.
  Must be symmetric and nonnegative definite
- `A::ScalarOrArray` : `n x n` coefficient on state in state transition
- `B::ScalarOrArray` : `n x k` coefficient on control in state transition
- `;C::ScalarOrArray{zero(size(R}(1)))` : `n x j` coefficient on random shock in
  state transition
- `;N::ScalarOrArray{zero(size(B,1)}(size(A, 2)))` : `k x n` cross product in
  payoff equation
- `;bet::Real(1.0)` : Discount factor in `[0, 1]`
- `capT::Union{Int, Void}(Void)` : Terminal period in finite horizon
  problem
- `rf::ScalarOrArray{fill(NaN}(size(R)...))` : `n x n` terminal payoff in finite
  horizon problem. Must be symmetric and nonnegative definite.

"""
function LQ(Q::ScalarOrArray,
            R::ScalarOrArray,
            A::ScalarOrArray,
            B::ScalarOrArray,
            C::ScalarOrArray=fill(zero(eltype(R)), size(R, 1)),
            N::ScalarOrArray=zero(B'A);
            bet::ScalarOrArray=1.0,
            capT::Union{Int,Nothing}=nothing,
            rf::ScalarOrArray=fill(NaN, size(R)...))

    k = size(Q, 1)
    n = size(R, 1)
    F = k==n==1 ? zero(Float64) : fill(zero(Float64), k, n)
    P = copy(rf)
    d = 0.0

    LQ(Q, R, A, B, C, N, bet, capT, rf, P, d, F)
end

@doc doc"""
Update `P` and `d` from the value function representation in finite horizon case

##### Arguments

- `lq::LQ` : instance of `LQ` type

##### Returns

- `P::ScalarOrArray` : `n x n` matrix in value function representation
  ``V(x) = x'Px + d``
- `d::Real` : Constant in value function representation

##### Notes

This function updates the `P` and `d` fields on the `lq` instance in addition to
returning them

"""
function update_values!(lq::LQ)
    # Simplify notation
    Q, R, A, B, N, C, P, d = lq.Q, lq.R, lq.A, lq.B, lq.N, lq.C, lq.P, lq.d

    # Some useful matrices
    s1 = Q + lq.bet * (B'P*B)
    s2 = lq.bet * (B'P*A) + N
    s3 = lq.bet * (A'P*A)

    # Compute F as (Q + B'PB)^{-1} (beta B'PA)
    lq.F = s1 \ s2

    # Shift P back in time one step
    new_P = R - s2'lq.F + s3

    # Recalling that tr(AB) = tr(BA)
    new_d = lq.bet * (d + tr(P * C * C'))

    # Set new state
    lq.P, lq.d = new_P, new_d
end

@doc doc"""
Computes value and policy functions in infinite horizon model.

##### Arguments

- `lq::LQ` : instance of `LQ` type

##### Returns

- `P::ScalarOrArray` : n x n matrix in value function representation ``V(x) = x'Px + d``
- `d::Real` : Constant in value function representation
- `F::ScalarOrArray` : Policy rule that specifies optimal control in each period

##### Notes

This function updates the `P`, `d`, and `F` fields on the `lq` instance in
addition to returning them

"""
function stationary_values!(lq::LQ)
    # simplify notation
    Q, R, A, B, N, C = lq.Q, lq.R, lq.A, lq.B, lq.N, lq.C

    # solve Riccati equation, obtain P
    A0, B0 = sqrt(lq.bet) * A, sqrt(lq.bet) * B
    P = solve_discrete_riccati(A0, B0, R, Q, N)

    # Compute F
    s1 = Q .+ lq.bet * (B' * P * B)
    s2 = lq.bet * (B' * P * A) .+ N
    F = s1 \ s2

    # Compute d
    d = lq.bet * tr(P * C * C') / (1 - lq.bet)

    # Bind states
    lq.P, lq.F, lq.d = P, F, d
end

"""
Non-mutating routine for solving for `P`, `d`, and `F` in infinite horizon model

See docstring for `stationary_values!` for more explanation
"""
function stationary_values(lq::LQ)
    _lq = LQ(copy(lq.Q),
             copy(lq.R),
             copy(lq.A),
             copy(lq.B),
             copy(lq.C),
             copy(lq.N),
             bet=copy(lq.bet),
             capT=lq.capT,
             rf=copy(lq.rf))

    stationary_values!(_lq)
    return _lq.P, _lq.F, _lq.d
end

"""
Private method implementing `compute_sequence` when state is a scalar
"""
function _compute_sequence(lq::LQ, x0::T, policies) where T
    capT = length(policies)

    x_path = Vector{T}(undef, capT+1)
    u_path = Vector{T}(undef, capT)

    x_path[1] = x0
    u_path[1] = -(first(policies)*x0)
    w_path = lq.C * randn(capT+1)

    for t = 2:capT
        f = policies[t]
        x_path[t] = lq.A*x_path[t-1] + lq.B*u_path[t-1] + w_path[t]
        u_path[t] = -(f*x_path[t])
    end
    x_path[end] = lq.A*x_path[capT] + lq.B*u_path[capT] + w_path[end]

    x_path, u_path, w_path
end

"""
Private method implementing `compute_sequence` when state is a scalar
"""
function _compute_sequence(lq::LQ, x0::Vector{T}, policies) where T
    # Ensure correct dimensionality
    n, j, k = size(lq.C, 1), size(lq.C, 2), size(lq.B, 2)
    capT = length(policies)

    A, B, C = lq.A, reshape(lq.B, n, k), reshape(lq.C, n, j)

    x_path = Matrix{T}(undef, n, capT+1)
    u_path = Matrix{T}(undef, k, capT)
    w_path = C*randn(j, capT+1)

    x_path[:, 1] = x0
    u_path[:, 1] = -(first(policies)*x0)

    for t = 2:capT
        f = policies[t]
        x_path[:, t] = A*x_path[: ,t-1] + B*u_path[:, t-1] + w_path[:, t]
        u_path[:, t] = -(f*x_path[:, t])
    end
    x_path[:, end] = A*x_path[:, capT] + B*u_path[:, capT] + w_path[:, end]

    x_path, u_path, w_path
end

@doc doc"""
Compute and return the optimal state and control sequence, assuming innovation ``N(0,1)``

##### Arguments

- `lq::LQ` : instance of `LQ` type
- `x0::ScalarOrArray`: initial state
- `ts_length::Integer(100)` : maximum number of periods for which to return process. If `lq` instance is finite horizon type, the sequenes are returned only for `min(ts_length, lq.capT)`

##### Returns

- `x_path::Matrix{Float64}` : An `n x T+1` matrix, where the t-th column represents ``x_t``
- `u_path::Matrix{Float64}` : A `k x T` matrix, where the t-th column represents ``u_t``
- `w_path::Matrix{Float64}` : A `n x T+1` matrix, where the t-th column represents `lq.C*N(0,1)`

"""
function compute_sequence(lq::LQ, x0::ScalarOrArray, ts_length::Integer=100)

    # Compute and record the sequence of policies
    if isa(lq.capT, Nothing)
        stationary_values!(lq)
        policies = fill(lq.F, ts_length)
    else
        capT = min(ts_length, lq.capT)
        policies = Vector{typeof(lq.F)}(undef, capT)
        for t = capT:-1:1
            update_values!(lq)
            policies[t] = lq.F
        end
    end

    _compute_sequence(lq, x0, policies)
end


text:
struct JuliaNotebookInlineDisplay <: AbstractDisplay end

function Base.display(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("image/png")}, x)
    payload = Base.invokelatest(Base64.stringmime, MIME("image/png"), x)
    JSONRPC.send_notification(conn_endpoint[], "notebook/display", Dict("items" => [Dict{String,Any}("mimetype" => "image/png", "data" => payload)]))
end

Base.displayable(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("image/png")}) = true

function Base.display(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("image/jpeg")}, x)
    payload = Base.invokelatest(Base64.stringmime, MIME("image/jpeg"), x)
    JSONRPC.send_notification(conn_endpoint[], "notebook/display", Dict("items" => [Dict{String,Any}("mimetype" => "image/jpeg", "data" => payload)]))
end

Base.displayable(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("image/jpeg")}) = true

function Base.display(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("image/svg+xml")}, x)
    payload = Base.invokelatest(Base64.stringmime, MIME("image/svg+xml"), x)
    JSONRPC.send_notification(conn_endpoint[], "notebook/display", Dict("items" => [Dict{String,Any}("mimetype" => "image/svg+xml", "data" => payload)]))
end

Base.displayable(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("image/svg+xml")}) = true

function Base.display(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v4+json")}, x)
    payload = Base.invokelatest(Base64.stringmime, MIME("application/vnd.vegalite.v4+json"), x)
    JSONRPC.send_notification(conn_endpoint[], "notebook/display", Dict("items" => [Dict{String,Any}("mimetype" => "application/vnd.vegalite.v4+json", "data" => payload)]))
end

Base.displayable(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("application/vnd.vegalite.v4+json")}) = true

function Base.display(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("text/html")}, x)
    payload = Base.invokelatest(Base64.stringmime, MIME("text/html"), x)
    JSONRPC.send_notification(conn_endpoint[], "notebook/display", Dict("items" => [Dict{String,Any}("mimetype" => "text/html", "data" => payload)]))
end

Base.displayable(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("text/html")}) = true

function Base.display(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("text/plain")}, x)
    payload = Base.invokelatest(Base64.stringmime, MIME("text/plain"), x)
    JSONRPC.send_notification(conn_endpoint[], "notebook/display", Dict("items" => [Dict{String,Any}("mimetype" => "text/plain", "data" => payload)]))
end

Base.displayable(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("text/plain")}) = true

function Base.display(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("text/markdown")}, x)
    payload = Base.invokelatest(Base64.stringmime, MIME("text/markdown"), x)
    JSONRPC.send_notification(conn_endpoint[], "notebook/display", Dict("items" => [Dict{String,Any}("mimetype" => "text/markdown", "data" => payload)]))
end

Base.displayable(d::JuliaNotebookInlineDisplay, ::MIME{Symbol("text/markdown")}) = true

function Base.display(d::JuliaNotebookInlineDisplay, x)
    things_to_show = IJuliaCore.display_dict(x)

    JSONRPC.send_notification(conn_endpoint[], "notebook/display", Dict("items" => [Dict{String,Any}("mimetype" => k, "data" => v) for (k, v) in things_to_show]))
end

max_stars_repo_name: FAlobaid/julia-vscode

max_issues_repo_name: FAlobaid/julia-vscode

max_forks_repo_name: FAlobaid/julia-vscode


text:
"""
function of the type {0, .., k-1} -> {0, ..., l-1}
represened internally as array of integers
"""
struct ResidueFunction <: BasicDiscreteFunction
    m :: Vector{Int}
    dom :: Segment{Int}
    rng :: Segment{Int}
end

"""
function of the type
ℤ_k1 × … × ℤ_km → ℤ_l
represened internally as array of integers
"""
struct ExtendedResidueFunction{T, N} <: BasicDiscreteFunction
    m :: Vector{Int}
    dom :: DirectProduct{T}
    rng :: Segment{Int}
    sizes :: NTuple{N, Int}
end

const ResFun = ResidueFunction
const ExtFun = ExtendedResidueFunction

# ResidueFunction from domain and codomain
function ResFun(dom :: Segment, codom :: Segment)
    m = zeros(Int, length(dom))
    return ResFun(m, dom, codom)
end
# constructing [0..k-1] -> [0..l-1]
# input : max_x, max_y
ResFun(k, l) = ResFun(Segment(k), Segment(l))

# from tuple of args
# input : (tuple) = (max_x, max_y)
ResFun(tup :: Tuple{T, S}) where {T, S <: Integer} = ResFun(tup[1], tup[2])

# constructiong {0, ..., k-1} -> {0, .., k-1}
ResFun(k) = ResFun(k, k)

# ExtFun from domain and codomain
function ExtFun(dom :: AbstractProduct, codom :: Segment)
    m = zeros(Int, length(dom))
    return ExtFun(m, dom, codom, size(dom))
end

# ExtFun from tuple (sizes) and max_y
function ExtFun(sizes :: NTuple, max_y :: Int)
    dom = DirectProduct(map(Segment, sizes)...)
    codom = Segment(max_y)
    return ExtFun(dom, codom)
end

(f :: ResFun)(x :: Integer) = f.m[x + 1]
Base.setindex!(f :: ResFun, value, key) = setindex!(f.m, value, key + 1)

function (f :: ExtFun)(x :: Union{NTuple, Vector})
    ind = itr2ind(x, f.sizes)
    return f.m[ind + 1]
end

function Base.setindex!(f :: ExtFun, v, k :: NTuple)
    ind = itr2ind(k, f.sizes)
    f.m[ind+1] = v
end
Base.setindex!(f :: ExtFun, v, k...) = setindex!(f, v, tuple(k...))

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@testset "test_parse_tonemapping_settings" begin
     string="""./Raytracer.jl tonemapping --normalization=0.1 --gamma=1.0 --avg_lum=1.0 """*
     """--ONLY_FOR_TESTS prova.pfm prova.png"""
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "tonemapping"
          @test 1==1
          tone_mapping(parse_tonemapping_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

##########################################################################################92

@testset "test_parse_demo_settings" begin
     string="""./Raytracer.jl demo --camera_type=per --camera_position="[1 , 2 , 3]" """*
     """--alpha=30. --width=40 --height=30 --normalization=0.1 --gamma=1.0 --avg_lum=1.0 """*
     """--world_type="B" --set_pfm_name="prova.pfm" --set_png_name="prova.png" --samples_per_pixel=16 --ONLY_FOR_TESTS onoff"""
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "demo"
          @test 1==1
          demo(parse_demo_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

@testset "test_parse_demo_onoff_settings" begin
     string="""./Raytracer.jl demo --ONLY_FOR_TESTS onoff --background_color="<1,2,3>" """*
     """--color="<4,  5,  6>" """
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "demo"
          @test 1==1
          demo(parse_demo_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

@testset "test_parse_demo_flat_settings" begin
     string="""./Raytracer.jl demo --ONLY_FOR_TESTS flat --background_color="<1,2,3>" """
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "demo"
          @test 1==1
          demo(parse_demo_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

@testset "test_parse_demo_pathtracer_settings" begin
     string="""./Raytracer.jl demo --ONLY_FOR_TESTS pathtracer --background_color="<1,2,3>" """*
     """--init_state=1 --init_seq=2 --num_of_rays=15 --max_depth=3 --russian_roulette_limit=1"""
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "demo"
          @test 1==1
          demo(parse_demo_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

@testset "test_parse_demo_pointlight_settings" begin
     string="""./Raytracer.jl demo --ONLY_FOR_TESTS pointlight --background_color="<1,2,3>" """*
     """--ambient_color="<4,  5,  6>" --dark_parameter=0.01"""
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "demo"
          @test 1==1
          demo(parse_demo_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

##########################################################################################92

@testset "test_parse_demoanimation_settings" begin
     string="""./Raytracer.jl demo_animation --camera_type="per" --camera_position="[1 , 2 , 3]" """*
     """--width=40 --height=30 --world_type="B" --set_anim_name="prova.mp4" """*
     """--samples_per_pixel=16 --ONLY_FOR_TESTS --normalization=0.1 --gamma=1.0 --avg_lum=1.0 onoff"""
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "demo_animation"
          @test 1==1
          demo_animation(parse_demoanimation_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

@testset "test_parse_demoanimation_onoff_settings" begin
     string="""./Raytracer.jl demo_animation --ONLY_FOR_TESTS onoff --background_color="<1,2,3>" """*
     """--color="<4,  5,  6>" """
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "demo_animation"
          @test 1==1
          demo_animation(parse_demoanimation_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

@testset "test_parse_demoanimation_flat_settings" begin
     string="""./Raytracer.jl demo_animation --ONLY_FOR_TESTS flat --background_color="<1,2,3>" """
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "demo_animation"
          @test 1==1
          demo_animation(parse_demoanimation_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

@testset "test_parse_demoanimation_pathtracer_settings" begin
     string="""./Raytracer.jl demo_animation --ONLY_FOR_TESTS pathtracer --background_color="<1,2,3>" """*
     """--init_state=1 --init_seq=2 --num_of_rays=15 --max_depth=3 --russian_roulette_limit=1"""
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "demo_animation"
          @test 1==1
          demo_animation(parse_demoanimation_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

@testset "test_parse_demoanimation_pointlight_settings" begin
     string="""./Raytracer.jl demo_animation --ONLY_FOR_TESTS pointlight --background_color="<1,2,3>" """*
     """--ambient_color="<4,  5,  6>" --dark_parameter=0.01 """
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "demo_animation"
          @test 1==1
          demo_animation(parse_demoanimation_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

##########################################################################################92

@testset "test_parse_render_settings" begin
     string="""./Raytracer.jl render --camera_type="per" --camera_position="[1 , 2 , 3]" """*
     """--alpha=30. --width=40 --height=30 --set_pfm_name="prova.pfm" """*
     """--set_png_name="prova.png" --samples_per_pixel=16 --ONLY_FOR_TESTS """*
     """--normalization=0.1 --gamma=1.0 --avg_lum=1.0 examples/tutorial_basic_syntax.txt onoff """
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "render"
          @test 1==1
          render(parse_render_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

@testset "test_parse_render_onoff_settings" begin
     string="""./Raytracer.jl render --ONLY_FOR_TESTS examples/tutorial_basic_syntax.txt onoff """*
     """--background_color="<1,2,3>" --color="<4,  5,  6>" """
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "render"
          @test 1==1
          render(parse_render_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

@testset "test_parse_render_flat_settings" begin
     string="""./Raytracer.jl render --ONLY_FOR_TESTS examples/tutorial_basic_syntax.txt flat """*
     """--background_color="<1,2,3>" """
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "render"
          @test 1==1
          render(parse_render_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

@testset "test_parse_render_pathtracer_settings" begin
     string="""./Raytracer.jl render --ONLY_FOR_TESTS examples/tutorial_basic_syntax.txt pathtracer """*
     """--background_color="<1,2,3>" --init_state=1 --init_seq=2 --num_of_rays=15 --max_depth=3 --russian_roulette_limit=1"""
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "render"
          @test 1==1
          render(parse_render_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

@testset "test_parse_render_pointlight_settings" begin
     string="""./Raytracer.jl render --ONLY_FOR_TESTS examples/tutorial_basic_syntax.txt pointlight """*
     """--background_color="<1,2,3>" --ambient_color="<4,  5,  6>" --dark_parameter=0.01"""
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "render"
          @test 1==1
          render(parse_render_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
end

##########################################################################################92

@testset "test_parse_render_animation_settings" begin
     cd("..")
     string="""./Raytracer.jl animation --camera_type="per" --camera_position="[1 , 2 , 3]" """*
     """--alpha=30. --width=40 --height=30 --set_anim_name="prova.mp4" """*
     """--normalization=0.1 --gamma=1.0 --avg_lum=1.0 --samples_per_pixel=16 --ONLY_FOR_TESTS """*
     """--function=my_function --vec_variables=[float] --iterable=1:2 """*
     """examples/tutorial_basic_syntax.txt onoff """
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "animation"
          @test 1==1
          render_animation(parse_render_animation_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
     cd("test")
end

@testset "test_parse_animation_onoff_settings" begin
     cd("..")
     string="""./Raytracer.jl animation --ONLY_FOR_TESTS --function=my_function --vec_variables=[float] --iterable=1:2 """*
     """examples/tutorial_basic_syntax.txt onoff """*
     """--background_color="<1,2,3>" --color="<4,  5,  6>" """
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "animation"
          @test 1==1
          render_animation(parse_render_animation_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
     cd("test")
end

@testset "test_parse_animation_flat_settings" begin
     cd("..")
     string="""./Raytracer.jl animation --ONLY_FOR_TESTS --function=my_function --vec_variables=[float] --iterable=1:2 """*
     """examples/tutorial_basic_syntax.txt flat """*
     """--background_color="<1,2,3>" """
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "animation"
          @test 1==1
          render_animation(parse_render_animation_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
     cd("test")
end

@testset "test_parse_animation_pathtracer_settings" begin
     cd("..")
     string="""./Raytracer.jl animation --ONLY_FOR_TESTS --function=my_function --vec_variables=[float] --iterable=1:2 """*
     """examples/tutorial_basic_syntax.txt pathtracer """*
     """--background_color="<1,2,3>" --init_state=1 --init_seq=2 --num_of_rays=15 --max_depth=3 --russian_roulette_limit=1"""
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "animation"
          @test 1==1
          render_animation(parse_render_animation_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
     cd("test")
end

@testset "test_parse_animation_pointlight_settings" begin
     cd("..")
     string="""./Raytracer.jl animation --ONLY_FOR_TESTS --function=my_function --vec_variables=[float] --iterable=1:2 """*
     """examples/tutorial_basic_syntax.txt pointlight """*
     """--background_color="<1,2,3>" --ambient_color="<4,  5,  6>" --dark_parameter=0.01"""
     args = Raytracing.from_CLI_to_vecstring(string)
     parsed_arguments = ArgParse_command_line(args)
	(isnothing(parsed_arguments)) && (return nothing)

	parsed_command = parsed_arguments["%COMMAND%"]
	parsed_settings = parsed_arguments[parsed_command]

     if parsed_command == "animation"
          @test 1==1
          render_animation(parse_render_animation_settings(parsed_settings)...)
          @test 1==1
     else
          @test 1==2
     end
     cd("test")
end


text:
#
#
# abstract CostModel
#
# function cost(model::CostModel, y::AVecF, yhat::AVecF)
#   c = 0.0
#   for i in 1:length(y)
#     c += cost(model, y[i], yhat[i])
#   end
#   # sum([cost(model, y[i], yhat[i]) for i in 1:length(y)])
#   c
# end
#
# # note: the vector version of costMultiplier also returns a boolean which is true when we
# #       need to multiply this value by f'(Σ) when calculating the sensitivities δ
# function costMultiplier!(model::CostModel, costmult::AVecF, y::AVecF, yhat::AVecF)
#   for i in 1:length(costmult)
#     costmult[i] = costMultiplier(model, y[i], yhat[i])
#   end
#   # Float64[costMultiplier(model, y[i], yhat[i]) for i in 1:length(y)], true
#   true
# end
#
#
# #-------------
#
# doc"""
# Typical sum of squared errors:
#   cost = sumabs2(yhat-y)
# """
# immutable L2CostModel <: CostModel end
#
# "cost function"
# cost(::L2CostModel, y::Float64, yhat::Float64) = 0.5 * (y - yhat) ^ 2
# # cost(::L2CostModel, y::AVecF, yhat::AVecF) = 0.5 * sumabs2(y - yhat)
#
# "returns M from the equation δ = M .* f'(Σ) ... used in the gradient update"
# costMultiplier(::L2CostModel, y::Float64, yhat::Float64) = yhat - y
#
# #-------------
#
# "cost = sumabs(yhat-y)"
# immutable L1CostModel <: CostModel end
#
# cost(::L1CostModel, y::Float64, yhat::Float64) = abs(y - yhat)
# costMultiplier(::L1CostModel, y::Float64, yhat::Float64) = sign(yhat - y)
#
# #-------------
#
# """
# Typical sum of squared errors, but scaled by ρ*y.  We implicitly assume that y ∈ {0,1}.
# """
# immutable WeightedL2CostModel <: CostModel
#   ρ::Float64
# end
# cost(model::WeightedL2CostModel, y::Float64, yhat::Float64) = 0.5 * (y - yhat) ^ 2 * (y > 0 ? model.ρ : 1)
# costMultiplier(model::WeightedL2CostModel, y::Float64, yhat::Float64) = (yhat - y) * (y > 0 ? model.ρ : 1)
#
# #-------------
#
# immutable CrossentropyCostModel <: CostModel end
#
# cost(model::CrossentropyCostModel, y::Float64, yhat::Float64) = -log(y > 0.0 ? yhat : (1.0 - yhat)) # binary case
# function cost(model::CrossentropyCostModel, y::AVecF, yhat::AVecF) # softmax case
#   length(y) == 1 && return cost(model, y[1], yhat[1])
#   C = 0.0
#   for (i,yi) in enumerate(y)
#     C -= yi * log(yhat[i])
#   end
#   C
# end
#
# costMultiplier(model::CrossentropyCostModel, y::Float64, yhat::Float64) = yhat - y # binary case
# function costMultiplier!(model::CrossentropyCostModel, costmult::AVecF, y::AVecF, yhat::AVecF) # softmax case
#   n = length(costmult)
#   if n == 1
#     costmult[1] = costMultiplier(model, y[1], yhat[1])
#   else
#     for i in 1:n
#       costmult[i] = yhat[i] - y[i]
#     end
#   end
#   # (length(y) == 1 ? Float64[costMultiplier(model, y[1], yhat[1])] : yhat - y), false
#   false
# end
#
#
# #-------------


text:
struct Adaptor end

# convert Metal Buffers Metal device pointers # should be generic
Adapt.adapt_storage(to::Adaptor, buf::MtlBuffer{T}) where {T} = reinterpret(Core.LLVMPtr{T,AS.Device}, buf.handle)

# Base.RefValue isn't GPU compatible, so provide a compatible alternative
struct MtlRefValue{T} <: Ref{T}
  x::T
end
Base.getindex(r::MtlRefValue) = r.x
Adapt.adapt_structure(to::Adaptor, r::Base.RefValue) = MtlRefValue(adapt(to, r[]))

"""
  mtlconvert(x)

This function is called for every argument to be passed to a kernel, allowing it to be
converted to a GPU-friendly format. By default, the function does nothing and returns the
input object `x` as-is.

Do not add methods to this function, but instead extend the underlying Adapt.jl package and
register methods for the `Metal.Adaptor` type.
"""
mtlconvert(arg) = adapt(Adaptor(), arg)


text:
ts = collect(now():Minute(10):(99*Minute(10)+now()))
N = length(ts)

xf = randn(N)
xt = 2*xf + 1 + 0.1*randn(N)

qf = Stage(ts,xf)
qt = Temperature(ts,xt)

c = Calibration(qt,qf)

m1 = fit(PolynomialCalibrationModel,c,1)

xf2 = randn(N)
qf2 = Stage(ts,xf2)

qt2 = predict(m1,qf2)

@testset "Polynomial calibration models" begin

    @test quantity(qt2) == m1.β[1] + m1.β[2]*xf2
    @test typeof(qt2) == Temperature{Float64}

end


text:
import .ROCKernels: ROCCTX, ROCCtx
import Cassette
import Enzyme

@inline function Cassette.overdub(::ROCCtx, ::typeof(Enzyme.autodiff_deferred), f, annotation::Enzyme.Annotation, args...)
    f′ = (args...) -> (Base.@_inline_meta; Cassette.overdub(ROCCTX, f, args...))
    Enzyme.autodiff_deferred(f′, annotation, args...)
end

@inline function Cassette.overdub(::ROCCtx, ::typeof(Enzyme.autodiff_deferred), f, args...)
    f′ = (args...) -> (Base.@_inline_meta; Cassette.overdub(ROCCTX, f, args...))
    Enzyme.autodiff_deferred(f′, args...)
end


text:
using JLD
dirname = "../data/sens/solenoid/dJds_vs_N"
function read_data()
    filenames = readdir(dirname)
	npts = size(filenames)[1]
	n_samples = zeros(Int64, npts)
	dJds = zeros(npts)
	for (i, filename) in enumerate(filenames)
			n_samples[i] = parse(Int64, split(filename, ".")[1])
			X = load(string(dirname,"/",filename))
			dJds[i] = X["dJds"][1]
	end

	save("../data/sens/solenoid/dJds3.jld",
		 "dJds", dJds, "n_samples", n_samples)
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
"""
    Solve method for DiffEq solution
"""
function get_de_ic(prob,eqs,u0=nothing)
    if u0 == nothing
        @info "Setting a random initial condition..."
        rand(eqs,prob.d)
    elseif typeof(u0) <: Number
        typeof(eqs) == CE2 ? fullrank(eqs,prob.d,u0) : rand(eqs,prob.d,u0)
    else
        @info "Converting QL solution to CE2 initial condition..."
        convert(eqs,u0,prob.d)
        # u0 = convert(eqs,u0,prob.d)
        # small noise to m = 1
        # for n = 1:2prob.d.ny-1
        #     u0.x[2][n,n,1] += 1e-6*exp(im*rand(Uniform(0,2π)))
        #     @show u0.x[2][n,n,1]
        # end
        # return u0
    end
end

function get_de_params(prob,eqs)::AbstractParams
    A = acoeffs(prob)
    B = bcoeffs(prob)
    C⁺,C⁻ = ccoeffs(prob,eqs)
    F = fcoeffs(prob,eqs)
    get_de_p(prob.d,eqs,[A,B,C⁺,C⁻,F])
end

get_de_p(d,eqs::NL,p) = NLParams(d.nx,d.ny,p...)
get_de_p(d,eqs::GQL,p) = GQLParams(d.nx,d.ny,eqs.Λ,p...)
get_de_p(d,eqs::CE2,p) = CE2Params(d.nx,d.ny,p...)
get_de_p(d,eqs::GCE2,p) = GCE2Params(d.nx,d.ny,eqs.Λ,p...)

get_de_probalg(prob,eqs,u0,t,p) = ODEProblem(f!,u0,t,p), DP5()
get_de_probalg(prob::BetaPlane{T,Stochastic{T}},eqs::CE2,u0,t,p) where T = ODEProblem(f!,u0,t,p), RK4()
get_de_probalg(prob::BetaPlane{T,Stochastic{T}},eqs,u0,t,p) where T = SDEProblem(f!,g!,u0,t,p), SRIW1()

get_de_kwargs(prob,eqs::AbstractEquations,tspan;kwargs...) = kwargs

function get_de_kwargs(prob,eqs::GCE2,tspan;kwargs...)
    if(!eqs.poscheck) return kwargs end
    @info "Setting positivity check callback..."
    poschecktimes = [tt for tt=tspan[1]:eqs.poscheckat:tspan[2]]
    condition(u,t,integrator) = t ∈ poschecktimes && !ispositive(u.x[2],prob.d.nx,prob.d.ny,eqs.Λ)
    affect!(integrator) = positivity!(integrator.u.x[2],prob.d.nx,prob.d.ny,eqs.Λ)
    poscheckcb = DiscreteCallback(condition,affect!,save_positions=(false,false))
    merge((callback=poscheckcb,tstops=poschecktimes),kwargs)
end

function get_de_kwargs(prob,eqs::CE2,tspan;kwargs...)
    if(!eqs.poscheck && !eqs.eigmax)
        return kwargs
    elseif(eqs.poscheck && !eqs.eigmax)
        @info "Setting positivity check callback..."
        poschecktimes = [tt for tt=tspan[1]:eqs.poscheckat:tspan[2]]
        condition(u,t,integrator) = t ∈ poschecktimes
        temp = zeros(ComplexF64,2*prob.d.ny-1,prob.d.nx-1,2*prob.d.ny-1,prob.d.nx-1)
        posaffect!(integrator) = positivity!(integrator.u.x[2],temp,prob.d.nx,prob.d.ny)
        poscondition(u,t,integrator) = t > 0.0
        poscheckcb = DiscreteCallback(poscondition,posaffect!,save_positions=(false,false))
        return merge((callback=poscheckcb,tstops=poschecktimes),kwargs)
    else
        @info "Removing all but the largest modal eigenvalues..."
        temp = zeros(ComplexF64,2*prob.d.ny-1,2*prob.d.ny-1,prob.d.nx-1)
        maxcondition(u,t,integrator) = true
        maxaffect!(integrator) = truncatecumulant!(prob.d,integrator.u.x[2],temp)
        eigmaxcb = DiscreteCallback(maxcondition,maxaffect!,save_positions=(false,false))
        return merge((callback=eigmaxcb,tspan=tspan),kwargs)
    end
end

function integrate(prob,eqs::AbstractEquations,tspan;u0=nothing,kwargs...)
    Random.seed!(123)
    _u0 = get_de_ic(prob,eqs,u0)
    _p  = get_de_params(prob,eqs)
    _prob,_alg = get_de_probalg(prob,eqs,_u0,tspan,_p)
    _kwargs = get_de_kwargs(prob,eqs,tspan;kwargs...)
    @time solve(_prob,_alg;_kwargs...)
end

integrate(prob,eqs::Vector{AbstractEquations},tspan;kwargs...) = [integrate(prob,eq,tspan;kwargs...) for eq in eqs]


text:
module MolecularIntegrals

include("Atoms.jl")
include("Basis.jl")
include("OneInts.jl")
include("Boys.jl")
include("Utils.jl")
include("ERI.jl")
include("HGP.jl")
include("Rys.jl")
#include("HGPgen2.jl")
include("Data.jl")

end # module


text:
mutable struct Engine{T}
    mem::Array{T,1}
    ptr::Int64
end

function create_engine(memsize::Integer, celsize::DataType)
    return Engine{celsize}(zeros(memsize), 1)
end

function inc_ptr(e::Engine)
    e.ptr += 1
end

function dec_ptr(e::Engine)
    e.ptr -= 1
end

function inc_dat(e::Engine)
    e.mem[e.ptr] += 1
end

function dec_dat(e::Engine)
    e.mem[e.ptr] -= 1
end

function output(e::Engine)
    return Char(e.mem[e.ptr])
end

function input(e::Engine, c::Char)
    e.mem[e.ptr] = Int32(c)
end

function ioutput(e::Engine)
    return e.mem[e.ptr]
end

function iinput(e::Engine, c::Integer)
    e.mem[e.ptr] = c
end

function copy_to_top(e::Engine)
    e.mem[length(e.mem)] = e.mem[e.ptr]
end

function copy_to_cell(e::Engine)
    e.mem[e.ptr] = e.mem[length(e.mem)]
end

function loop(e::Engine)
    return e.mem[e.ptr] != 0
end


text:
module SequenceTrimmer

using BioSequences, CodecZlib, DataStructures

export consensus, estimate_libtype, pwm_summary, trim_fastq_threads
include("trimming.jl")

end # module


text:
## iteration_tools.jl
## These functions ensure that also numbers can be iterated and zipped
export linear_index
export apply_tuple_list

"""
    cast_iter(vals::Matrix)

A convenience function used to make number iterable via repeated.
This is used when dealing with multiple arguments where one can be
optionally an iterable.
A Tuple is also interpreted as repeating this tuple over and over.
"""
function cast_iter(vals::Matrix)
    Tuple(Tuple(vals[:,n]) for n in 1:size(vals,2))
    # Tuple(vals[:,axes(vals, 2)])
end

function cast_iter(vals::Vector{<:Number})
    Tuple(vals)   # makes the matrix iterable, interpreting it as a series of vectors
end

function cast_iter(vals::IterType)
    vals
end

function cast_iter(vals)
    repeated(vals) # during the zip operation this type always yields the same value
end

"""
    cast_number_iter(vals::Matrix)

a conveniance function used to make number itarable via repeated. This is used when dealing with multiple arguments where one can be optionally an iterable.
The behavious is slightly different to cast_iter when it comes to Tuples
A Tuple is here iterated as a list of numbers.
"""
function cast_number_iter(vals::Vector{<:Number})
    Tuple(vals)   # makes the matrix iterable, interpreting it as a series of vectors
end

function cast_number_iter(vals::NTuple{N,<:Number} where N)
    vals
end

function cast_number_iter(vals::Number)
    repeated(vals) # during the zip operation this type always yields the same value
end

"""
    optional_mat_to_iter(vals)  # only for matrices

casts a matrix via cast_iter, all other types are unaffected
"""
function optional_mat_to_iter(vals)  # only for matrices
    vals
end

function optional_mat_to_iter(vals::Matrix)
    cast_iter(vals)
end

"""
    mat_to_tvec(v)
converts a 2d matrix to a Vector of Tuples

Example:
```jldoctest
julia> NDTools.mat_to_tvec([1 3 5; 2 4 6])
3-element Vector{Tuple{Int64, Int64}}:
 (1, 2)
 (3, 4)
 (5, 6)
```
"""
mat_to_tvec(v) = [Tuple(v[:,n]) for n in 1:size(v,2)] # converts a 2d matrix to a Vector of Tuples

"""
    apply_dims(scale, dims, N)

replaces scale entries not in dims with zeros. Can also be applied to an Iterable of Tuples

Arguments:
+ scale: Vector to pass its valid entries through
+ dims: Tuple of valid dimensions
+ N: total length of the resulting  `NTuple`

Example:
```jldoctest
julia> NDTools.apply_dims((1,2,3),(1,2), 4)
(1, 2, 0, 0)
```
"""
function apply_dims(scale, dims, N)  # replaces scale entries not in dims with zeros
    ntuple(i -> i ∈ dims ? scale[i] : zero(eltype(scale)), N)
end

function apply_dims(scales::IterType, dims, N)  # replaces scale entries not in dims with zeros
    Tuple(ntuple(i -> i ∈ dims ? scale[i] : zero(eltype(scale)), N)  for scale in  scales)
end

"""
    linear_index(pos, sz)

converts a tuble (pos) to a linear index using the size (sz).

Example:
```doctest
julia> a = rand(10,10,10);

julia> linear_index((3,4,5), (10,10,10))
433

julia> a[3,4,5]  == a[433]
true
```
"""
function linear_index(pos, sz)
    factors = (1, cumprod(sz[1:end-1])...)
    sum((pos.-1) .* factors)+1
end

# Functions from scaling_offset_types in IndexFunArrays

"""
    apply_tuple_list(f, t1,t2)

applies a two-argument function to tubles and iterables of tuples, if either of the arguments `t1` or `t2` is an Iterable
"""
function apply_tuple_list(f, t1, t2)  # applies a two-argument function to tubles and iterables of tuples
    return f(t1,t2)
end

function apply_tuple_list(f, t1, t2::IterType)
    return Tuple([f(t1,a2) for a2 in t2])
end

function apply_tuple_list(f, t1::IterType, t2)
    res= Tuple([f(a1,t2) for a1 in t1])
    return res
end

function apply_tuple_list(f, t1::IterType, t2::IterType)
    return Tuple([f(a[1],a[2]) for a in zip(t1,t2)])
end


text:
"""
Only observe the product

What do we want? The briggs-haldane approximation: ddt [ES] = 0. This means that kon/(koff + kcat) is practically unidentifiable.

https://chem.libretexts.org/Bookshelves/Biological_Chemistry/Supplemental_Modules_(Biological_Chemistry)/Enzymes/Enzymatic_Kinetics/Michaelis-Menten_Kinetics_1
"""

function MichaelisMenten(input)

  @parameters t
  @parameters kon koff kcat
  D = Differential(t)
  @variables pos(t) vel(t) inp(t)

  eqs = [D(pos) ~ vel,
        D(vel) ~ input(t) - (1/m)*(c*vel + k*pos),
  ]

  ps = [k,c,m] .=>  [1.,4.,1.]
  ics = [pos, vel] .=> [1.,0.]
  od = ODESystem(eqs, t, first.(ics),first.(ps))
  # tspan = (0.,100.)
  # prob = ODEProblem(od, ics, tspan, ps)
  return od, ics, tspan, ps
end


text:
module MechanicalUnits
using Unitfu
# Exported infix function / operator
export ∙
# Import / exports for short and parseable type signatures
import Unitfu: Time, Length, Mass, Temperature, Current, Luminosity, Amount
import Unitfu: ᵀ , ᴸ , ᴹ , ᶿ, ᴶ , ᴺ
import Unitfu: lookup_units, promote_to_derived
export Time, Length, Mass, Temperature, Current, Luminosity, Amount, Level
export ᵀ , ᴸ , ᴹ , ᶿ , ᴶ , ᴺ
export Quantity, DimensionlessQuantity, NoUnits, NoDims
import Unitfu:
        FreeUnits, AffineUnits, Affine, AffineQuantity, Unitlike, Unit, Dimension, Dimensions, Units
export  FreeUnits, AffineUnits, Affine, AffineQuantity, Unitlike, Unit, Dimensions, Dimension, Units
export  Level, Gain

# For importing from Unitfu, or defining more units
export @import_expand, @unit, @u_str

# Reexported functions from Unitfu
export logunit, unit, absoluteunit, dimension, uconvert, ustrip, upreferred, ∙
export uconvertp, uconvertrp, reflevel, linear

# Useful functions that are not exported by Unitfu.
export preferunits, convfact, promote_to_derived

# A vector of all the exported units. This is printed during precompilation.
export MECH_UNITS

import Unitfu: isunitless, sortexp, showrep, abbr, prefix, power, superscript, tens, genericunit
import Unitfu: promote_unit, preferunits, convfact
# derived dimensions
import Unitfu: Area, Acceleration, Force, Pressure, Density
import Unitfu: Velocity
import Unitfu:ForceFreeUnits, PressureFreeUnits, EnergyFreeUnits, AreaFreeUnits, DensityFreeUnits, VolumeFreeUnits
export Area, Acceleration, Force, Pressure, Density, Velocity

# Extend base. This could perhaps reside in Unitfu
import Base: tryparse_internal, parse

# Units are exported in 'import_export_units.jl'.

include("internal_functions.jl")
include("import_export_units.jl")
# We have defined and exported e.g. m². Now do the same for dimension symbbols,
# so that e.g.  ᵀ² == ᵀ ^². This way, output could be used as constructors.
eval(exponents_superscripts(:ᵀ))
eval(exponents_superscripts(:ᴸ))
eval(exponents_superscripts(:ᴹ))
eval(exponents_superscripts(:ᶿ))
eval(exponents_superscripts(:ᴶ))
eval(exponents_superscripts(:ᴺ))

# Used for registering units with Unitfu macros during initialisation.
const localunits = Unitfu.basefactors

include("parse.jl")
function __init__()
    # This is for evaluating Unitfu macros in the context of this package.
    merge!(Unitfu.basefactors, localunits)
    # This enables any units defined here to be used in the @u_str
    Unitfu.register(MechanicalUnits)
    # This pre-selects some useful units for the mechanical engineering domain
    preferunits(kg, mm, s, K)
    Sys.iswindows() && push!(ENV, "UNITFUL_FANCY_EXPONENTS" => "true")
    Sys.isapple() && push!(ENV, "UNITFUL_FANCY_EXPONENTS" => "true")
end

end # module

max_stars_repo_name: hustf/MechanicalUnits

max_issues_repo_name: hustf/MechanicalUnits


text:
using CellListMap
using StaticArrays
import Random

#
# In this test we compute the average displacement of the x coordinates of the particles
#
function average_displacement(;N=100_000,parallel=true,x=nothing)

    # Number of particles, sides and cutoff
    sides = @SVector [250,250,250]
    cutoff = 10
    box = Box(sides, cutoff)

    # Particle positions
    Random.seed!(321)
    if x === nothing
        x = [ sides .* rand(SVector{3,Float64}) for i in 1:N ]
    end

    # Initialize auxiliary linked lists
    cl = CellList(x, box, parallel=parallel)

    # Function to be evalulated for each pair: sum of displacements on x
    f(x, y, avg_dx) = avg_dx + abs(x[1] - y[1])

    avg_dx = (N / (N * (N - 1) / 2)) * map_pairwise(
        (x, y, i, j, d2, avg_dx) -> f(x, y, avg_dx),
        0.,box,cl,
        parallel=parallel
    )

    return avg_dx
end


text:
using ModelingToolkit
using OrdinaryDiffEq, DiffEqCallbacks
using BenchmarkTools
using LinearAlgebra
using DataFrames

@parameters t, I[1:3, 1:3], Iinv[1:3, 1:3], g[1:3], m, α
@variables q[1:4](t) x[1:3](t) ω[1:3](t) v[1:3](t)
@derivatives D'~t

function transition_matrix(q)
    q0, q1, q2, q3 = q
    P = [
        2*(q0^2 + q1^2) - 1 2*(q1*q2 - q0*q3)   2*(q1*q3 + q2*q0)
        2*(q1*q2 + q0*q3)   2*(q0^2 + q2^2) - 1 2*(q2*q3 - q1*q0)
        2*(q1*q3 - q2*q0)   2*(q2*q3 + q1*q0)   2*(q0^2 + q3^2) - 1
    ]
end

function quaternion_derivative(q, ω)
    a, b, c = ω
    Mₛ = [
        0 -a -b -c
        a  0  c -b
        b -c  0  a
        c  b -a  0
    ]
    1/2 * Mₛ * q
end

function dynamics(q, x, ω, v)
    P = transition_matrix(q)
    vcat(
        D.(q) .~ quaternion_derivative(q, ω),
        D.(x) .~ P * v, # world velocity
        D.(ω) .~ Iinv * cross(I * ω, ω),
        D.(v) .~ 1/m * (m * P' * g - α .* v) # body acceleration
    )
end

function postprocess(u) # dynamics_f does not work
    q, x, ω, v = u[1:4, end], u[5:7, end], u[8:10, end], u[11:13, end]
    dynamics_f(q, x, ω, v)
end

get_values(eqs) = getproperty.(getproperty.(eqs, :rhs), :value)
function extract_result(u)
    syms = [:q1, :q2, :q3, :q4, :x, :y, :z, :ωx, :ωy, :ωz, :vx, :vy, :vz]
    df = DataFrame(repeat([Float64], 13), syms)
    for i ∈ 1:size(u)[2]
        push!(df, Dict(syms .=> u[:, i]))
    end
    df
end

function quaternion_norm_callback(resid, u, p, t)
    resid[1] = @views norm(u[1:4])^2 - 1
    resid[2:13] .= 0
end

quaternion_manifold_proj = ManifoldProjection(quaternion_norm_callback)

u0 = vcat(
    q .=> [1., 0, 0, 0],
    x .=> 0.,
    ω .=> [0.1, 0.15, 0.2],
    v .=> [1., 0, 0],
)

p = vcat(
    (I .=> Diagonal([0.25, 0.2, 0.1]))...,
    # (Iinv .=> inv(I))...
    (Iinv .=> Diagonal([4, 5., 10.]))...,
    g .=> [0., 0., -9.81],
    m => 3.,
    α => 0.01
)

tspan = (0.0, 10000.0)
sys = ODESystem(dynamics(q, x, ω, v))
prob = ODEProblem(sys, u0, tspan, p)

u = solve(prob, Vern7(), save_everystep=false, callback=quaternion_manifold_proj)
df = extract_result(u)

# dynamics_f = build_function(dynamics(q, x, ω, v), q, x, ω, v; expression=Val{false})
# postprocess(u)


text:
module Check
export @check

struct FailedCheck
    result
    ex
    meta
    FailedCheck(r::ANY, e::ANY, m::ANY) = new(r, e, m)
end

@noinline function check_fail(result::ANY, orig_ex::ANY)
    throw(FailedCheck(result, orig_ex, nothing))
end

@noinline function check_expr_fail(result::ANY, orig_ex::Expr, values::ANY)
    throw(FailedCheck(result, orig_ex, values))
end

function Base.show(io::IO, check::FailedCheck)
    println(io, "FailedCheck(")
    if check.result === false
        println(io, "! ")
        show(io, check.ex)
    else
        show(io, check.ex)
        print(io, " => ")
        show(io, check.result)
        print(io, ")")
    end
end

function Base.showerror(io::IO, check::FailedCheck)
    println(io, "Code check failed:")
    println(io, "Expression: ", check.ex)
    if !(check.result === false)
        println(io, " Got: ", check.result)
    end
    if !(check.meta === nothing)
        println(io, " Evaluated: ", check.meta)
    end
end

# An internal function, called by the code generated by the @check
# macro to get results of the test expression.
# In the special case of a comparison, e.g. x == 5, generate code to
# evaluate each term in the comparison individually so the results
# can be displayed nicely.
function get_test_result(ex::Expr)
    # Normalize non-dot comparison operator calls to :comparison expressions
    if ex.head == :call && length(ex.args) == 3 &&
        first(string(ex.args[1])) != '.' &&
        (ex.args[1] === :(==) || Base.operator_precedence(ex.args[1]) == Test.comparison_prec)
        return get_test_result(Expr(:comparison, ex.args[2], ex.args[1], ex.args[3]))
    end

    if (ex.head == :comparison || ex.head == :call)
        # move all terms of the call into SSA position
        nterms = length(ex.args)
        argret = Expr(:tuple)
        newex = Expr(ex.head)
        testret = Expr(:block)
        resize!(argret.args, nterms)
        resize!(newex.args, nterms)
        resize!(testret.args, nterms + 1)
        for i = 1:nterms
            flatex = ex.args[i]
            argsym = gensym()
            if isa(flatex, Expr)
                flatex, argvals = get_test_result(flatex)
                push!(argvals.args, argsym)
            else
                flatex = esc(flatex)
                argvals = argsym
            end
            argret.args[i] = argvals
            testret.args[i] = Expr(:(=), argsym, flatex)
            newex.args[i] = argsym
        end
        testret.args[nterms + 1] = newex
        return testret, argret
    end
    return esc(ex), Expr(:tuple)
end

macro check(ex::Symbol)
    orig_ex = Expr(:inert, ex)
    return quote
        if (result = $ex) !== true
            check_fail(result, $orig_ex)
        end
    end
end

macro check(expr::Expr, kws...)
    # my lazy error checking
    # something between an enhanced @assert
    # and a less-fancy @test
    Test.test_expr!("@test", expr, kws...)
    orig_expr = Expr(:inert, expr)
    resultexpr, valsexpr = get_test_result(expr)
    return quote
        let result = $resultexpr
            if result !== true
                let values = $valsexpr
                    check_expr_fail(result, $orig_expr, values)
                end
            end
        end
    end
end

end

max_stars_repo_name: vtjnash/PerfPlots

max_issues_repo_name: vtjnash/PerfPlots

max_forks_repo_name: vtjnash/PerfPlots


text:
abstract type SimulationModel end

struct MassMatrixModel <: SimulationModel end
struct ResidualModel <: SimulationModel end


text:
using Test
using ScanDir

# Tests copied and modified from julia's test/file.jl
dirwalk = mktempdir()
cd(dirwalk) do
    for i=1:2
        mkdir("sub_dir$i")
        open("file$i", "w") do f end

        mkdir(joinpath("sub_dir1", "subsub_dir$i"))
        touch(joinpath("sub_dir1", "file$i"))
    end
    touch(joinpath("sub_dir2", "file_dir2"))
    has_symlinks = !Sys.iswindows() || (Sys.windows_version() >= Sys.WINDOWS_VISTA_VER)
    follow_symlink_vec = has_symlinks ? [true, false] : [false]
    has_symlinks && symlink(abspath("sub_dir2"), joinpath("sub_dir1", "link"))
    for prune_subsub in [false, true]
        prune = prune_subsub ? e->startswith(e.name,"subsub") : _->false
        subsubs = prune_subsub ? [] : ["subsub_dir1", "subsub_dir2"]

        for follow_symlinks in follow_symlink_vec
            chnl = ScanDir.walkdir(".", follow_symlinks=follow_symlinks, prune=prune)

            root, dirs, files = take!(chnl)
            @test root == "."
            @test dirs == ["sub_dir1", "sub_dir2"]
            @test files == ["file1", "file2"]

            root, dirs, files = take!(chnl)
            @test root == joinpath(".", "sub_dir1")
            if has_symlinks
                if follow_symlinks
                    @test dirs ==  ["link"; subsubs]
                    @test files == ["file1", "file2"]
                else
                    @test dirs ==  subsubs
                    @test files == ["file1", "file2", "link"]
                end
            else
                @test dirs ==  subsubs
                @test files == ["file1", "file2"]
            end

            root, dirs, files = take!(chnl)
            if follow_symlinks
                @test root == joinpath(".", "sub_dir1", "link")
                @test dirs == []
                @test files == ["file_dir2"]
                root, dirs, files = take!(chnl)
            end
            if !prune_subsub
                for i=1:2
                    @test root == joinpath(".", "sub_dir1", "subsub_dir$i")
                    @test dirs == []
                    @test files == []
                    root, dirs, files = take!(chnl)
                end
            end

            @test root == joinpath(".", "sub_dir2")
            @test dirs == []
            @test files == ["file_dir2"]

            @test !isready(chnl)
        end

        for follow_symlinks in follow_symlink_vec
            chnl = ScanDir.walkdir(".", follow_symlinks=follow_symlinks, topdown=false, prune=prune)
            root, dirs, files = take!(chnl)
            if follow_symlinks
                @test root == joinpath(".", "sub_dir1", "link")
                @test dirs == []
                @test files == ["file_dir2"]
                root, dirs, files = take!(chnl)
            end
            if !prune_subsub
                for i=1:2
                    @test root == joinpath(".", "sub_dir1", "subsub_dir$i")
                    @test dirs == []
                    @test files == []
                    root, dirs, files = take!(chnl)
                end
            end
            @test root == joinpath(".", "sub_dir1")
            if has_symlinks
                if follow_symlinks
                    @test dirs ==  ["link"; subsubs]
                    @test files == ["file1", "file2"]
                else
                    @test dirs ==  subsubs
                    @test files == ["file1", "file2", "link"]
                end
            else
                @test dirs ==  subsubs
                @test files == ["file1", "file2"]
            end

            root, dirs, files = take!(chnl)
            @test root == joinpath(".", "sub_dir2")
            @test dirs == []
            @test files == ["file_dir2"]

            root, dirs, files = take!(chnl)
            @test root == "."
            @test dirs == ["sub_dir1", "sub_dir2"]
            @test files == ["file1", "file2"]

            @test !isready(chnl)
        end
    end

    # These subtly depend on timing, so removed for now
    # TODO need better way to test onerror
    # #test of error handling
    # chnl_error = ScanDir.walkdir(".")
    # chnl_noerror = ScanDir.walkdir(".", onerror=x->x)

    # root, dirs, files = take!(chnl_error)
    # #@show root
    # @test root == "."
    # @test dirs == ["sub_dir1", "sub_dir2"]
    # @test files == ["file1", "file2"]

    # rm(joinpath("sub_dir1"), recursive=true)
    # @test_throws Base.IOError take!(chnl_error) # throws an error because sub_dir1 do not exist

    # root, dirs, files = take!(chnl_noerror)
    # @test root == "."
    # @test dirs == ["sub_dir1", "sub_dir2"]
    # @test files == ["file1", "file2"]

    # root, dirs, files = take!(chnl_noerror) # skips sub_dir1 as it no longer exist
    # @test root == joinpath(".", "sub_dir2")
    # @test dirs == []
    # @test files == ["file_dir2"]

    # Test that symlink loops don't cause errors
    if has_symlinks
        mkdir(joinpath(".", "sub_dir3"))
        foo = joinpath(".", "sub_dir3", "foo")
        # this symlink requires admin privileges on Windows:
        symlink_err = false
        try
            symlink("foo", foo)
        catch e
            symlink_err = true
            showerror(stderr, e); println(stderr)
            @warn "Could not create symlink. Symlink loop tests skipped."
        end

        if !symlink_err
            let files = ScanDir.walkdir(joinpath(".", "sub_dir3"); follow_symlinks=true)
                # wording of this error differs on julia 1.5, so just check for IOError
                #@test_throws Base._UVError("stat($(repr(foo)))", Base.UV_ELOOP)  take!(files)
                @test_throws Base.IOError  take!(files)
            end
            root, dirs, files = take!(ScanDir.walkdir(joinpath(".", "sub_dir3"); follow_symlinks=false))
            @test root == joinpath(".", "sub_dir3")
            @test dirs == []
            @test files == ["foo"]
        end
    end
end
rm(dirwalk, recursive=true)


text:
type NeuralNetwork
    hidden::Union{Vector,Int64}
    act::String
    weights::Dict{Integer, Matrix}
    max_iter::Integer
end

function NeuralNetwork(;
                       hidden::Union{Vector,Int64} = [2,2,1],
                       act::String = "sigmoid",
                       weights::Dict{Integer,Matrix} = Dict{Integer,Matrix}(),
                       max_iter::Integer = 500000)
    return NeuralNetwork(hidden, act, weights, max_iter)
end

function train!(model::NeuralNetwork, X::Matrix, y::Vector)
    if typeof(model.hidden) <: Integer
        model.hidden = [size(X, 2), model.hidden, size(y, 2)]
    end
    init_weights(model)
    depth = size(model.hidden,1)
    a::Dict{Integer, Vector} = Dict()
    z::Dict{Integer, Vector} = Dict()
    X = hcat(X,ones(size(X,1)))
    for i = 1:model.max_iter
        #batch
        r = rand(1:size(X,1))
        a[1] = X[r,:]
        for j = 2:(depth)
            z[j] = vec(a[j-1]' * model.weights[j-1])
            a[j] = vec(sigmoid(z[j]))
        end
        delta = Dict{Integer, Vector}()
        error_ = a[depth] - y[r]
        #if i % 1000 == 0
        #    println("$(i) epochs: error $(error_)")
        #end
        delta[depth] = error_ .* sigmoid_prime(z[depth])
        for j = depth-1:-1:2
            delta[j] = vec(delta[j+1]' * model.weights[j]') .* sigmoid_prime(z[j])
        end
        for j = 1:depth-1
            del = delta[j+1]
            a_temp = a[j]
            model.weights[j] = model.weights[j] - 0.1 * a_temp * del'
        end
    end

end

function init_weights(model::NeuralNetwork)
    depth_ = size(model.hidden,1)
    for i = 1:(depth_-2)
        model.weights[i] = 2*rand(model.hidden[i]+1,model.hidden[i+1]+1)-1
    end
    model.weights[depth_-1] = 2*rand(model.hidden[depth_-1]+1,model.hidden[depth_])-1
end

function predict(model::NeuralNetwork,
                 x::Matrix)
    n = size(x,1)
    m = model.hidden[end]
    res = zeros(n,m)
    for i = 1:n
        res[i,:] = predict(model, x[i,:])
    end
    return res
end

function predict(model::NeuralNetwork,
                 x::Vector)
    x = x'
    x = hcat(x,[1])
    for i = 1:length(model.hidden)-1
        x = sigmoid(x * model.weights[i])
    end
    return x
end

## fixed

function test_NeuralNetwork()
    #xor_network
    X_train = [0.1 0.1; 0.1 0.9; 0.9 0.1; 0.9 0.9]
    y_train = [0.1,0.9,0.9,0.1]
    model = NeuralNetwork(hidden=10)
    train!(model,X_train, y_train)
    predictions = predict(model,X_train)
    print("regression msea", mean_squared_error(y_train, predictions))
end


text:
using JLD2
using FileIO

# Network codes
include("../backbones/resnet.jl")
include("../backbones/fpn.jl")
include("../backbones/context_module.jl")
# BBox Related Codes
include("../utils/box_processes.jl")
# Loss Functions
include("../core/losses.jl")
# Global Configuration Parameters
include("../../configs.jl")

"""
- Takes final Context Head Outputs and converts these into proposals.
- Same structure for BBox, Classifier and Landmark tasks.
- task_len is 2 for classification, 4 for bbox and 10 for landmark points.
"""
mutable struct HeadGetter layers; task_len; scale_cnt; end

function HeadGetter(input_dim, task_len; scale_cnt=5, dtype=Array{Float32})
    layers = []
    num_anchors = scale_cnt == 3 ? 2 : 3
    for s in 1:scale_cnt
        push!(layers, Conv2D(1, 1, input_dim, num_anchors*task_len, dtype=dtype, bias=true))
    end
    return HeadGetter(layers, task_len, scale_cnt)
end

function (hg::HeadGetter)(xs; train=true)
    proposals = []; T = hg.task_len;
    for (i, x) in enumerate(xs)
        proposal = hg.layers[i](x, train=train)
        W, H, C, N = size(proposal); A = Int(W*H*(C/T));
        # converting all proposals from 4D shape to 3D
        proposal = permutedims(proposal, (3, 1, 2, 4))
        proposal = reshape(proposal, (T, A, N))
        push!(proposals, proposal)
    end
    if T == 2
        return softmax(cat(proposals..., dims=2), dims=1)
    else
        return cat(proposals..., dims=2)
    end
end

"""
Our actual model that predicts bounding boxes and landmarks.
"""
mutable struct RetinaFace
    backbone; fpn; context_module;
    cls_head1; cls_head2;
    bbox_head1; bbox_head2;
    lm_head1; lm_head2;
    mode; num_anchors; anchor_info; log_file;
    dtype
end

function RetinaFace(;mode=1, num_anchors=3, anchor_info=lat_5_anchors, load_path=nothing, dtype=Array{Float32})
    laterals = num_anchors === 2 ? 3 : 5
    if load_path !== nothing
        model = load_model(load_path)
        return RetinaFace( # restoring only weights and creating new objects with added additional params
            model.backbone,
            FPN(
                model.fpn.o6, model.fpn.o2, model.fpn.o3, model.fpn.o4, model.fpn.o5,
                model.fpn.merge4, model.fpn.merge3, model.fpn.merge2, laterals
            ),
            ContextModule(
                model.context_module.ssh_p2, model.context_module.ssh_p3, model.context_module.ssh_p4,
                model.context_module.ssh_p5, model.context_module.ssh_p6, laterals
            ),
            HeadGetter(model.cls_head1.layers, model.cls_head1.task_len, laterals),
            HeadGetter(model.cls_head2.layers, model.cls_head2.task_len, laterals),
            HeadGetter(model.bbox_head1.layers, model.bbox_head1.task_len, laterals),
            HeadGetter(model.bbox_head2.layers, model.bbox_head2.task_len, laterals),
            HeadGetter(model.lm_head1.layers, model.lm_head1.task_len, laterals),
            HeadGetter(model.lm_head2.layers, model.lm_head2.task_len, laterals),
            mode, num_anchors, anchor_info, nothing, dtype
        )
    else
        backbone = load_mat_weights(
            ResNet50(include_top=false, dtype=dtype), "./weights/imagenet-resnet-50-dag.mat"
        )
        return RetinaFace(
            backbone, FPN(dtype=dtype, scale_cnt=laterals), ContextModule(dtype=dtype, scale_cnt=laterals), # full baseline
            HeadGetter(256, 2, scale_cnt=laterals, dtype=dtype), HeadGetter(256, 2, scale_cnt=laterals, dtype=dtype),
            HeadGetter(256, 4, scale_cnt=laterals, dtype=dtype), HeadGetter(256, 4, scale_cnt=laterals, dtype=dtype),
            HeadGetter(256, 10, scale_cnt=laterals, dtype=dtype), HeadGetter(256, 10, scale_cnt=laterals, dtype=dtype),
            mode, num_anchors, anchor_info, nothing, dtype
        )
    end
end

# modes:
# 0 --> for getting p_vals
# 1 --> first context head forward,
# 2 --> second context head forward,
# context_module --> enables or disables the usage of context module
function (model::RetinaFace)(x; p_vals = nothing, mode=0, context_module=true, train=true)

    if p_vals === nothing
        c2, c3, c4, c5 = model.backbone(x, return_intermediate=true, train=train)
        p_vals = model.fpn([c2, c3, c4, c5], train=train)
        if context_module
            p_vals = model.context_module(p_vals, train=train)
        end
    end

    class_vals = nothing; bbox_vals = nothing; landmark_vals = nothing;

    if mode == 0
        return p_vals

    elseif mode == 1
        class_vals = model.cls_head1(p_vals, train=train)
        bbox_vals = model.bbox_head1(p_vals, train=train)
        landmark_vals = model.lm_head1(p_vals, train=train)

    elseif mode == 2
        class_vals = model.cls_head2(p_vals, train=train)
        bbox_vals = model.bbox_head2(p_vals, train=train)
        landmark_vals = model.lm_head2(p_vals, train=train)
    end

    return class_vals, bbox_vals, landmark_vals
end

# modes:
# 0  --> baseline until FPN + 2nd head,
# 1  --> until context module + cascaded structure (full model),
# 2  --> until context module + 2nd head (no cascaded structure)
function (model::RetinaFace)(x, y, mode=0, train=true, weight_decay=0)

    use_context = mode == 0 ? false : true
    p_vals = model(x, mode=0, train=train, context_module=use_context)
    priors = _get_priorboxes(model.num_anchors, model.anchor_info, size(x, 1));

    cls_vals = nothing; bbox_vals = nothing; lm_vals = nothing;
    h1c_loss = 0; h1b_loss = 0; h1l_loss = 0; # first context head losses
    h2c_loss = 0; h2b_loss = 0; h2l_loss = 0; # second context head / baseline losses
    decay_loss = 0; # decay loss if decay value is bigger than 0

    if model.mode == 1
        # do the forward pass and calculate first head loss
        cls_vals1, bbox_vals1, lm_vals1 = model(x, mode=1, p_vals=p_vals, context_module=use_context, train=train)
        h1c_loss, h1b_loss, h1l_loss = get_loss(cls_vals1, bbox_vals1, lm_vals1, y, priors, mode=1, dtype=model.dtype)
        priors = _decode_bboxes(convert(Array{Float32}, value(bbox_vals1)), priors)
    end

    cls_vals2, bbox_vals2, lm_vals2 = model(x, mode=2, p_vals=p_vals, context_module=use_context, train=train)
    h2c_loss, h2b_loss, h2l_loss = get_loss(cls_vals2, bbox_vals2, lm_vals2, y, priors, mode=2, dtype=model.dtype)

    if weight_decay > 0 # only taking weights but not biases and moments
        for p in params(model)
            if size(size(p), 1) == 4 && size(p, 4) > 1
                decay_loss += weight_decay * sum(p .* p)
            end
        end
    end

    loss_cls = h1c_loss + h2c_loss
    loss_bbox = h1b_loss + h2b_loss
    loss_pts = h1l_loss + h2l_loss

    loss = loss_cls + lambda1 * loss_bbox + lambda2 * loss_pts + decay_loss

    to_print = get_losses_string(loss, loss_cls, loss_bbox, loss_pts, decay_loss)
    print(to_print);
    if model.log_file !== nothing
        open(model.log_file, "a") do io write(io, to_print) end; # saving data
    end

    return loss
end

# if mode is 1, then first head IOUs are taken, otherwise second head IOUs
function get_loss(cls_vals, bbox_vals, lm_vals, y, priors; mode=2, dtype=Array{Float32})

    loss_cls = 0; loss_lm = 0; loss_bbox = 0; loss_decay = 0; lmN = 0; bboxN = 0;
    pos_thold = mode == 1 ? head1_pos_iou : head2_pos_iou
    neg_thold = mode == 1 ? head1_neg_iou : head2_neg_iou

    # helper parameters for calculating losses
    T, P, N = size(cls_vals); batch_cls = convert(Array{Int64}, zeros(2, P, N));
    batch_gt = cat(value(bbox_vals), value(lm_vals), dims=1)

    for n in 1:N # loop for each input in batch, all inputs may have different box counts
        if isempty(y[n]) || (y[n] == []) || (y[n] === nothing)
            continue # if the cropped image has no faces
        end

        l_cls =Array(value(cls_vals))[1,:,n]; gt = y[n];

        prior = priors
        if length(size(priors)) > 2
            prior = priors[:,:,n] # retuired for mode 1
        end

        gt, pos_idx, neg_idx = encode_gt_and_get_indices(gt, prior, l_cls, pos_thold, neg_thold)

        if pos_idx !== nothing
            # if boxes with high enough IOU are found
            lm_indices = getindex.(findall(gt[15,:] .>= 0)); gt = convert(dtype, gt);

            if size(lm_indices, 1) > 0
                # counting only the ones with landmark data
                batch_gt[5:14,pos_idx[lm_indices],n] = gt[5:14,lm_indices]
                lmN += length(lm_indices)
            end

            batch_gt[1:4,pos_idx,n] = gt[1:4,:]; bboxN += length(pos_idx);
            batch_cls[1,neg_idx,n] .= 1; batch_cls[2,pos_idx,n] .= 2;
        end
    end

    # in case no boxes are matched in the whole batch
    bboxN = bboxN == 0 ? 1 : bboxN
    lmN = lmN == 0 ? 1 : lmN

    # classification negative log likelihood loss
    cls_vals = reshape(cls_vals, (2, N*P))
    loss_cls_neg = nll(cls_vals, vec(batch_cls[1,:,:]))
    loss_cls_pos = nll(cls_vals, vec(batch_cls[2,:,:]))
    # loss_cls = (ohem_ratio * loss_cls_pos + loss_cls_neg) / (ohem_ratio+1)
    loss_cls = (loss_cls_neg + loss_cls_pos) / 2
    if (isinf(value(loss_cls)) || isnan(value(loss_cls))) loss_cls = 0 end

    # regression loss of the box centers, width and height
    loss_bbox = smooth_l1(abs.(batch_gt[1:4,:,:] - bbox_vals)) / bboxN

    # box center and all landmark points regression loss per point
    loss_pts = smooth_l1(abs.(batch_gt[5:end,:,:] - lm_vals)) / lmN

    return loss_cls, loss_bbox, loss_pts
end

# modes:
# 0  --> baseline until FPN + 2nd head,
# 1  --> until context module + cascaded structure (full model),
# 2  --> until context module + 2nd head (no cascaded structure)
# set filter to false for not making any confidence score and NMS check (for evaluation)
function predict_image(model::RetinaFace, x; mode=1, filter=true, verbose=true, nms_thold=0.4, conf_thold=0.5, topk=5000)
    use_context = mode == 0 ? false : true
    img_size = size(x, 1)
    # getting predictions
    p_vals = model(x, mode=0, context_module=use_context, train=false)
    priors = _get_priorboxes(model.num_anchors, model.anchor_info, img_size)

    if mode == 1
        _, bbox_vals1, _ = model(x, mode=1, p_vals=p_vals, context_module=use_context, train=false)
        priors = _decode_bboxes(convert(Array{Float32}, value(bbox_vals1)), priors)[:,:,1]
    end

    cls_vals, bbox_vals, lm_vals = model(x, mode=2, p_vals=p_vals, context_module=use_context, train=false)
    cls_vals = Array(cls_vals); bbox_vals = Array(bbox_vals); lm_vals = Array(lm_vals);

    # decoding points to min and max
    bbox_vals, lm_vals = decode_points(bbox_vals, lm_vals, priors)
    bbox_vals = _to_min_max_form(bbox_vals)

    if !filter
        return cls_vals[:,:,1], bbox_vals[:,:,1], lm_vals[:,:,1]
    else
        indices = findall(cls_vals[2,:,1] .>= conf_thold)
        cls_result = cls_vals[:,indices,1]
        bbox_result = bbox_vals[:,indices,1]
        lm_result = lm_vals[:,indices,1]
        if verbose
            print("[INFO] Passed Confidence Check: ", size(indices, 1), "\n")
        end

        # for keeping the top k many elements
        if size(cls_result, 2) > topk
            indices = partialsortperm(cls_result[2,:,1], 1:topk, rev=true)
            cls_result = cls_result[:,indices,1]
            bbox_result = bbox_result[:,indices,1]
            lm_result = lm_result[:,indices,1]
        end

        indices = nms(vec(cls_result[2,:]), bbox_result, thold=nms_thold)
        cls_result = cls_result[:,indices]
        bbox_result = bbox_result[:,indices]
        lm_result = lm_result[:,indices]
        if verbose
            print("[INFO] Passed NMS Check: ", size(indices, 1),"\n")
            print("[INFO] Returning results above confidence level: ", conf_thold, "\n")
        end
        return cls_result, bbox_result, lm_result
    end
end

function train_model(model::RetinaFace, reader; val_data=nothing, save_dir=nothing, start_epoch=1, log_file=nothing)

    if log_file !== nothing
        model.log_file = log_file
        open(log_file, "w") do io write(io, "===== TRAINING PROCESS =====\n") end;
    end

    # Adjusting LR for each step
    lrs_per_epoch = zeros(num_epochs); lr_change = 0;
    for e in 1:num_epochs
        if e == 1
            lrs_per_epoch[e] = lrs[1]
        elseif e == lr_change_epoch[1]
            lrs_per_epoch[e] = lrs[2]
        elseif e == lr_change_epoch[2]
            lrs_per_epoch[e] = lrs[3]
        elseif e == lr_change_epoch[3]
            lrs_per_epoch[e] = lrs[4]
        else
            lrs_per_epoch[e] = lrs_per_epoch[e - 1]
        end
    end

    for e in start_epoch:num_epochs
        (imgs, boxes), state = iterate(reader)
        iter_no = 1; last_loss = 0; total_batches = size(state, 1) + size(imgs)[end];
        curr_batch = 0; curr_lr = lrs_per_epoch[e];

        print_lr = 0
        for p in params(model)
            if p.opt !== nothing && p.opt.lr !== nothing && lrs_per_epoch[e] != p.opt.lr
                p.opt = nothing
            end
            if p.opt === nothing
                print_lr = curr_lr
            else
                print_lr = p.opt.lr
            end
        end

        while state !== nothing && imgs !== nothing
            if mod(iter_no, 5) == 1 # prints per 5 batches
                to_print  = "\n--- Epoch: " * string(e) * " & LR: " * string(round(print_lr; digits=4))
                to_print *= " & Batch: " * string(curr_batch) * "/" * string(total_batches) * "\n\n"
                print(to_print)
                if model.log_file !== nothing
                    open(model.log_file, "a") do io write(io, to_print) end;
                end
            end

            # Updating the model
            momentum!(model, [(imgs, boxes, model.mode, true, weight_decay)], lr=curr_lr, gamma=momentum)
            (imgs, boxes), state = iterate(reader, state)
            iter_no += 1;
            if imgs !== nothing curr_batch += size(imgs)[end] end
        end
        if save_dir !== nothing
            to_save = save_dir * "model_"
            if model.mode == 0
                to_save *= "mode0_epoch" * string(e) * ".jld2"
            elseif model.mode == 1
                to_save *= "mode1_epoch" * string(e) * ".jld2"
            elseif model.mode == 2
                to_save *= "mode2_epoch" * string(e) * ".jld2"
            else
                to_save *= "unknownmode_epoch" * string(e) * ".jld2"
            end
            save_model(model, to_save)
        end
    end
end

function load_model(file_name)
    return Knet.load(file_name, "model")
end

function save_model(model::RetinaFace, file_name)
    Knet.save(file_name, "model", model)
end

function get_losses_string(total_loss, loss_cls, loss_bbox, loss_lm, loss_decay)
    total = string(round.(value(total_loss); digits=3))
    cls = string(round.(value(loss_cls); digits=3))
    bbox = string(round.(value(loss_bbox); digits=3))
    lm = string(round.(value(loss_lm); digits=3))
    decay = string(round.(value(loss_decay); digits=3))

    if length(total) < 6 total *= "0"^(6-length(total)) end
    if length(cls) < 6 cls *= "0"^(6-length(cls)) end
    if length(bbox) < 6 bbox *= "0"^(6-length(bbox)) end
    if length(lm) < 6 lm *= "0"^(6-length(lm)) end
    if length(decay) < 6 decay *= "0"^(6-length(decay)) end

    to_print  = "Total Loss: " *  total * " | "
    to_print *= "Cls Loss: "   * cls    * " | "
    to_print *= "Box Loss: "   * bbox   * " | "
    to_print *= "Point Loss: " * lm     * " | "
    to_print *= "Decay: "      * decay  * "\n"
    return to_print
end

max_stars_repo_name: barisbatuhan/retinaface

max_issues_repo_name: barisbatuhan/retinaface

max_forks_repo_name: barisbatuhan/retinaface


text:
# Finite Horizon algorithm evaluating problem with use of other POMDP algorithms
# User has to define (PO)MDP instance with use of POMDPs interface methods isterminal, reward, actions, states, stateindex, actions, actionindex, discount, transition

# For example see problems defined in POMDPModels package or FiniteHorizonPOMDP/test/instances/1DFiniteHorizonGridWorld.jl
#
# Currently supporting: ValueIterationSolver
#
# Possible future improvements: Change arrays from row oriented to column oriented
#                               Create fixhorizon(m::Union{MDP,POMDP}, T::Int) utility

# Policy struct created according to one used in DiscreteValueIteration

#TODO: Dosctring
struct FiniteHorizonSolver <: Solver
    verbose::Bool
end

function FiniteHorizonSolver(;verbose::Bool=false)
    return FiniteHorizonSolver(verbose)
end

#TODO: Dosctring
mutable struct FiniteHorizonValuePolicy{Q<:Array, U<:Array, P<:Array, A<:Array, M<:MDP} <: Policy
    qmat::Q
    util::U
    policy::P
    action_map::A
    include_Q::Bool
    m::M
end

# Policy constructor
function FiniteHorizonValuePolicy(m::MDP)
    no_stages = horizon(m)
    no_actions = length(actions(m))

    qmat =  Array{Matrix{Float64}}(undef, no_stages)
    util =  Array{Vector{Float64}}(undef, no_stages)
    policy = Array{Vector{Int64}}(undef, no_stages)
    return FiniteHorizonValuePolicy(qmat, util, policy, ordered_actions(m), true, m)
end

function value(policy::FiniteHorizonValuePolicy, s::S) where S
    stg = stage(policy.m, s)
    stg == horizon(policy.m) + 1 && return 0.
    sidx = stage_stateindex(policy.m, s)
    return policy.util[stg][sidx]
end

function action(policy::FiniteHorizonValuePolicy, s::S) where S
    stg = stage(policy.m, s)
    stg == horizon(policy.m) + 1 && return policy.action_map[1]
    sidx = stage_stateindex(policy.m, s)
    aidx = policy.policy[stg][sidx]
    return policy.action_map[aidx]
end

@POMDP_require solve(solver::FiniteHorizonSolver, mdp::MDP) begin
    M = typeof(mdp)
    S = statetype(M)
    A = actiontype(M)
    @req discount(::M)
    @subreq ordered_states(mdp)
    @subreq ordered_actions(mdp)
    @req isterminal(::M, ::S)
    @req transition(::M,::S,::A)
    @req reward(::M,::S,::A,::S)
    @req stateindex(::M,::S)
    @req actionindex(::M, ::A)
    @req actions(::M, ::S)
    as = actions(mdp)
    ss = states(mdp)
    @req length(::typeof(ss))
    @req length(::typeof(as))
    a = first(as)
    s = first(ss)
    dist = transition(mdp, s, a)
    D = typeof(dist)
    @req support(::D)
    @req pdf(::D,::S)

    @req stage(::M, ::S)
    E = typeof(stage(mdp, s))
    # @req HorizonLength(::M)
    @req stage_states(::M, ::Int64)
    @req stage_stateindex(::M, ::S)
    @req horizon(::M)
end

"""
    addstagerecord(fhpolicy::FiniteHorizonValuePolicy, qmat, util, policy, stage)

Store record for given stage results to `FiniteHorizonPolicy` and return updated version.
"""
function addstagerecord(fhpolicy::FiniteHorizonValuePolicy, qmat, util, policy, stage, solver::FiniteHorizonSolver)
    fhpolicy.qmat[stage] = qmat
    fhpolicy.util[stage] = util
    fhpolicy.policy[stage] = policy
    return fhpolicy
end

# MDP given horizon 5 assumes that agent can move 5 times
function solve(solver::FiniteHorizonSolver, m::MDP)
    if typeof(HorizonLength(m)) == InfiniteHorizon
        throw(ArgumentError("Argument m should be valid Finite Horizon MDP with methods from FiniteHorizonPOMDPs.jl interface implemented. If you are completely sure that you implemented all of them, you should also check if you have defined HorizonLength(::Type{<:MyFHMDP})"))
    end

    fhpolicy = FiniteHorizonValuePolicy(m)
    util = zeros(length(stage_states(m, horizon(m) + 1)))

    @showprogress 1 "Computing..." for stage=horizon(m):-1:1
        if solver.verbose
            println("Stage: $stage")
        end

        stage_q, util, pol = valueiterationsolver(m, stage, util)

        fhpolicy = addstagerecord(fhpolicy, stage_q, util, pol, stage, solver)

        if solver.verbose
            println("POLICY\n")
            println("QMAT")
            println(fhpolicy.qmat[stage])
            println("util")
            println(fhpolicy.util[stage])
            println("policy")
            println(fhpolicy.policy[stage])
            println("\n\n\n")
        end
    end

    return fhpolicy
end


text:
using NonlinearEigenproblems
using NonlinearEigenproblems.RKHelper
using LinearAlgebra
using SparseArrays
using Random

function particle_init(interval)
    nep, brpts, U0 = particle_nep(interval)

    # interval
    sep = 1e-4
    if interval == 1
        xmin = -U0
        xmax = brpts[interval] - sep
        Ξ = 10 .^ range(-6, stop = 6, length = 10000) .+ brpts[interval]
    elseif interval > 1
        xmin = brpts[interval-1] + sep
        xmax = brpts[interval] - sep
        Ξ1 = -10 .^ range(-6, stop = 6, length = 5000) .+ brpts[interval-1]
        Ξ2 = 10 .^ range(-6, stop = 6, length = 5000) .+ brpts[interval]
        Ξ = [Ξ1; Ξ2]
    else
        error("Invalid interval: $interval")
    end

    # define target set Σ
    Σ = [xmin + 0im, xmax + 0im]

    # options
    temp_pep=nep_gallery("pep0",200); v=complex(vcat(vec(temp_pep.A[1][:,1:81]),temp_pep.A[1][1:81,82])) # Ugly but gives stability over versions. Stable "random" vector v.
    nodes = range(xmin + 0im, stop = xmax + 0im, length = 11)
    nodes = collect(nodes[2:2:end])

    return nep, Σ, Ξ, v, nodes, xmin, xmax
end

function particle_nep(interval)
    # Original MATLAB version by William Vandenberghe

    # initialization
    # constants
    meter = 1 / 5.2917725e-11
    nm = 1e-9 * meter
    eV = 1 / 13.6

    # parameters
    xmax = 5
    zmax = 2

    xstep = 0.05
    zstep = 0.05

    x_x = (-xmax:xstep:xmax) * nm
    z_z = (-zmax:zstep:zmax) * nm

    nx = length(x_x)
    nz = length(z_z)

    dx = minimum(diff(x_x))
    dz = minimum(diff(z_z))

    x = kron(x_x, z_z.^0)
    z = kron(x_x.^0, z_z)

    w1 = 1 * nm
    w2 = 1.1 * nm
    l = 4 * nm
    U0 = 3 * eV

    # potential U
    U = zero(x)
    U[abs.(z) .< w1] .= -U0
    U[(abs.(z) .< w2) .& (abs.(x) .< l/2)] .= -U0

    # m
    m = 0.2

    # branch points
    n = nx * nz
    Dxx_x = SymTridiagonal(ones(nx) * -2 / dx^2, ones(nx-1) / dx^2)
    Dzz_z = SymTridiagonal(ones(nz) * -2 / dx^2, ones(nz-1) / dz^2)
    H_L = -1/m*Dzz_z + diagm(0 => U[1:nz])
    H_R = -1/m*Dzz_z + diagm(0 => U[end-nz.+(1:nz)])

    if opnorm(H_L - H_R) == 0
        # symmetric potential
        D,V = eigen(H_L)
        i = sortperm(D)
        d = D[i]
        V = V[:,i]
        p = zeros(nz) # 0 (left+right)
    else
        # asymmetric potential
        D_L,V_L = eigen(H_L)
        D_R,V_R = eigen(H_R)
        V = [V_L; V_R]
        p = [-ones(nz); ones(nz)]
        DD = [D_L; D_R]
        i = sortperm(DD)
        d = DD[i]
        V = V[:,i]
        p = p[i] # -1 (left), 1 (right)
    end

    # matrix H
    H = -1/m * (kron(sparse(Dxx_x), sparse(1.0I, nz, nz)) + kron(sparse(1.0I, nx, nx), sparse(Dzz_z))) + sparse(Diagonal(U))

    # branch points
    brpts = unique(d)

    # polynomial matrices
    B = [H, sparse(-1.0I, n, n)]

    # nonlinear matrices
    SL = Vector(undef, length(brpts))
    if p[1] < 0
        SL[1] = [V[:,1]; spzeros(n-nz, 1)]
    elseif p[1] == 0
        SL[1] = hcat([V[:,1]; spzeros(n-nz, 1)], [spzeros(n-nz, 1); V[:,1]])
    else
        SL[1] = [spzeros(n-nz, 1); V[:,1]]
    end
    c = 1;
    for j = 2:length(d)
        if d[j-1] == d[j]
            if p[j] < 0
                SL[c] = hcat(SL[c], [V[:,j]; spzeros(n-nz, 1)])
            elseif p[j] == 0
                SL[c] = hcat(SL[c], [V[:,j]; spzeros(n-nz, 1)], [spzeros(n-nz, 1); V[:,j]])
            else
                SL[c] = hcat(SL[c], [spzeros(n-nz, 1); V[:,j]])
            end
        else
            c += 1
            if p[j] < 0
                SL[c] = [V[:,j]; spzeros(n-nz, 1)]
            elseif p[j] == 0
                SL[c] = hcat([V[:,j]; spzeros(n-nz, 1)], [spzeros(n-nz, 1); V[:,j]])
            else
                SL[c] = [spzeros(n-nz, 1); V[:,j]]
            end
        end
    end
    SU = copy(SL)
    for j = 1:length(SL)
        SL[j] = -1 / m / dx^2 * SL[j]
    end

    # nonlinear functions
    f = Vector(undef, length(brpts))
    for j = 1:interval-1
        coeff=brpts[j]
        f[j] = lambda -> exp(im * sqrt((m * (lambda - coeff * one(lambda)))))
    end
    for j = interval:length(brpts)
        coeff=brpts[j]
        f[j] = lambda -> exp(-sqrt((m * (-lambda + coeff * one(lambda)))))
    end

    # finally assemble nep instance; note that the nonlinear matrices are defined by their LU factors only
    C = map(k -> LowRankMatrixAndFunction(similar(SL[k], 0, 0), SL[k], SU[k], f[k]), 1:length(f))
    nep = SumNEP(PEP(B), LowRankFactorizedNEP(C))

    return nep, brpts, U0
end


text:
#=
This is an implementation of the "Babylonian method" of approximating square roots.
(https://en.wikipedia.org/wiki/Methods_of_computing_square_roots#Babylonian_method)
=#
const steps = 10
const example = 20740

function sqrt(s)
    x = s/2 # arbitrary initial guess for x ≈ √s
    for i in 1:steps
        x = (x + s/x)/2
    end
    x
end

println("√$example ≈ $(sqrt(example)), error = $(abs(sqrt(example)) - (example)^0.5)")

max_stars_repo_name: landrumb/julia

max_issues_repo_name: landrumb/julia

max_forks_repo_name: landrumb/julia


text:
function EclipticProvider(arg0::IERSConventions)
    return EclipticProvider((IERSConventions,), arg0)
end

function EclipticProvider(arg0::IERSConventions, arg1::TimeScales)
    return EclipticProvider((IERSConventions, TimeScales), arg0, arg1)
end

function get_transform(obj::EclipticProvider, arg0::AbsoluteDate)
    return jcall(obj, "getTransform", Transform, (AbsoluteDate,), arg0)
end

function get_transform(obj::EclipticProvider, arg0::FieldAbsoluteDate)
    return jcall(obj, "getTransform", FieldTransform, (FieldAbsoluteDate,), arg0)
end