text:
using LinearAlgebra

"""
# enum SvdMode

. leftsv means that alogrithm operates on columns of data
passed in LowRankApproxMc and thus  giving
a orthonormal basis of left singular vector.

. rightsv means that algorithm stores the transpose of data passed
in LowRankApproxMc and thus runs on row of original data , thus giving
a orthonormal basis of right singular vector.

"""
@enum SvdMode begin
    leftsv
    rightsv
end

"""
# LowRankApproxMc structure

as described in:

1. Vempala Spectral Book (2009)
    paragraph on Iterative Fast Svd P 85 and Th 7.2 7.3

2. Lectures in Randomized Linear Algebra. Mahoney 2016 P121 and seq.

Cf also Deshpande Rademacher Vempala Wang
Matrix Approximation and Projective Clustering via Volume Sampling 2006

For a data matrix X(d,n) compute a orthogonal set of left singular vectors in a matrix W(rank, n)
Right singular vector are computed as W'*X.

FIELDS
------

 - X  :  the (d,n) matrix we want to reduce.
 - targetRank :   the rank of the data matrix we want to extract

 - currentRank : the rank reached at current iteration

 - U : Union{ Some{Matrix{Float64}}, Nothing} the matrix of left singular vectors if any
 - S : Union{ Some{Vector{Float64}}, Nothing} the vecor of singular values if any
 - Vt: Union{ Some{Matrix{Float64}}, Nothing} the transposed matrix of right singular vectors if any

CONSTRUCTORS
-----------

1.    `function LowRankApproxMc(data::Matrix{Float64}, targetRank::Int64)`
       initialize structure and rank of the approximation

"""
mutable struct LowRankApproxMc
    X::Matrix{Float64}
    #
    leftOrRight::SvdMode
    # rank we want in output
    targetRank::Int64
    targetEpsil::Float64
    # number of column sampling in a pass
    samplingSize::Int64
    #
    currentRank::Int64
    # store relative error
    residualErr::Float64
    #
    currentIter::Int64
    U::Union{ Some{Matrix{Float64}}, Nothing}
    S::Union{ Some{Vector{Float64}}, Nothing}
    Vt::Union{ Some{Matrix{Float64}}, Nothing}
    #
    function LowRankApproxMc(data::Matrix{Float64}; mode = leftsv)
        if mode == leftsv
            @info "LowRankApproxMc in left singular vector mode"
            new(data, mode, 0, 1. , 0,  0, 1. , 0,
                nothing, nothing, nothing)
        else
            @info "LowRankApproxMc in right singular vector mode, transposing data"
            new(transpose(data), mode, 0, 1. , 0,  0, 1. , 0,
                nothing, nothing, nothing)
        end
    end
end

#
# The data matrix are supposed to be stored in columns of length data dimension.
# We operate on columns so we ouptut Left singular vectors.
# So our implementation conforms to Mahoney "Lecture Notes on Randomized Algebra"
# paragraph 15.4 Th 21 as opposed to Kannan-Vempala "Spectral Algorithms"
# which operates on rows and output right singular vectors.
#

function reduceIfast(datapb::LowRankApproxMc, stop::Function)
    @debug "in LowRankApproxMc reduceIfast v2"
    epsil = 1.E-5
    d,n = size(datapb.X)
    datapb.currentIter = 1
    datapb.currentRank = 0
    E = Matrix(datapb.X)
    acceptanceRatio = 1.
    normX = norm(E,2)
    residualAtIter= Vector{Float64}()
    meanL2NormC = sqrt(normX * normX / n)
    @debug "meanL2NormC" meanL2NormC
    # We could use c as in Deshpande Vempala and let c decrease with iterations
    # to accelerate ?
    probaC = Vector{Float64}(undef,n)
    normC = Vector{Float64}(undef,n)
    A = Matrix{Float64}(undef, d, 0)
    more = true
    while more
        @debug "\n\n iteration " datapb.currentIter
        # sampled columns at this iter
        tSet = Set{Int64}()
        keptIdxNonZero = Vector{Int64}()
        # sample columns according to L2-norm of each columns
        normE = 0
        Threads.@threads for i = 1:n
            normC[i] = norm(E[:,i]) * norm(E[:,i])
        end
        for i = 1:n
            probaC[i] = i > 1 ? probaC[i-1] + normC[i] : normC[i]
            normE += normC[i]
        end
        if probaC[end] <= 0
            @debug "\n non positive proba for last slot=  :  " probaC[end]
            error("cannot sample column")
        end
        map!(x -> x/probaC[end], probaC, probaC)
        # get minimum non null column norm at iter t
        meanL2NormIter = mean(sqrt.(normC))
        @debug "\n c meanL2Norm at Iter  :  " meanL2NormIter
        for j = 1:datapb.samplingSize
            xsi = rand(Float64)
            s = searchsortedfirst(probaC, xsi)
            if s > length(probaC)
                @warn "\n sampling column failed xsi =  :  " xsi
                exit(1)
            end
            # avoid selecting some already sampled column or already a
            if !(s in tSet)
                push!(tSet, s)
            end
        end
        acceptanceRatio = length(tSet)/datapb.samplingSize
        @debug "acceptanceRatio" acceptanceRatio
        # get index of selected columns
        idxc = collect(tSet)
        At = datapb.X[:, idxc]
        # orthogonalize ... with preceding columns in A. Compute transpose(At)*A
        AttA = BLAS.gemm('T', 'N', At, A)
        Threads.@threads for i = 1:length(idxc)
            # with preceding columns in A. At[:,i] not normalized, A[:,j] is
            @simd for j = 1:size(A)[2]
                At[:,i] = At[:,i] - A[:,j] * AttA[i,j]
            end
        end
        @debug "QR orthogonalizing new block size At" size(At)
        At,tau  = LAPACK.geqrf!(At)
        Q=LAPACK.orgqr!(At, tau, length(tau))
        for i = 1:size(Q)[2]
            # within this block of columns in A.
            normi = norm(Q[:,i])
            # how much remains of At[:,i]. Should check with initial norm of column i
            if abs(normi - 1.) < 1.E-5
                push!(keptIdxNonZero, i)
#                @debug "adding a column vector of norm" normi idxc[i]
                A = hcat(A, Q[:,i])
                if size(A)[2] >=  min(d,n)
                    @debug "max dim reached"
                    more = false
                    break
                end
            else
                @debug "column not normalized" normi
            end
        end
        @info "A size" size(A)
        #
        datapb.currentRank = size(A)[2]
        # stopping criteria
        if stop(datapb)
            more = false
        else
            @debug "updating E"
            if length(keptIdxNonZero) < length(idxc)
                E = E - At[:, keptIdxNonZero] * BLAS.gemm('T', 'N', At[:, keptIdxNonZero], datapb.X)
            else
                E = E - At * BLAS.gemm('T', 'N', At, datapb.X)
            end
            datapb.residualErr = norm(E,2)/normX
            push!(residualAtIter, datapb.residualErr)
            # protection against stationarity
            if datapb.currentIter > 3 && residualAtIter[end] > 0.95 * residualAtIter[end-1]
                @warn "stationary iterations , stopping. perhaps increase rank"
                more = false
            end
            @info " iteration  , ||E|| set size "  datapb.currentIter datapb.residualErr  datapb.currentRank
            # if in reduceToEpsil mode we have to test now that we have datapb.residualErr
            if stop(datapb)
                more = false
            else
                datapb.currentIter += 1
            end
        end
    end
    # go from left singular to right singular vectors
    #
    @debug "computing right vectors"
    Vt = Matrix{Float64}(undef, size(A)[2], n)
    BLAS.gemm!('T', 'N', 1., A, datapb.X, 0., Vt)
    eigenValues = zeros(size(A)[2])
    Threads.@threads for i in 1:size(A)[2]
        eigenValues[i] = norm(Vt[i,:])
        Vt[i,:] = Vt[i,:] / eigenValues[i]
    end
    @debug "exiting LowRankApproxMc reduceIfast"
    datapb.U = Some(A)
    datapb.S = Some(eigenValues)
    datapb.Vt = Some(Vt)
end

"""
# function getResidualError(lrmc::LowRankApproxMc)

if X is the data matrix and U computation succeeded
    return the 2-uple (error, relerr) with :
```math
        error = norm(X - U * transpose(U) * X)
```
and :
```math
        relerr = Error/norm(X)
```

else throws an error
"""
function getResidualError(lrmc::LowRankApproxMc)
    U = something(lrmc.U)
    resErr = norm(lrmc.X - U* BLAS.gemm('T', 'N', U, lrmc.X))
    return resErr , resErr/norm(lrmc.X)
end

"""
# function getindex(LowRankApproxMc, Symbol)

returns U,S,Vt,V or rank according to Symbol :S , :U, :Vt , V or :k

"""
function getindex(datapb::LowRankApproxMc, d::Symbol)
    @match d begin
    :S    =>  datapb.S
    :U    =>  datapb.U
    :V    =>  datapb.Vt'
    :Vt   =>  datapb.Vt
    :k    =>  length(datapb.S)
    _     =>  throw(KeyError(d))
    end
end

"""
# function reduceToRank(datapb::LowRankApproxMc, expectedRank::Int64; nbiter = 2)

Implements incremental range approximation as described in
1: Vempala Spectral Book (2009)
    paragraph on Iterative Fast Svd P 85 and Th 7.2 7.3.

## Args
 -  expectedRank asked for.
    It must be noted that the returned rank can be different
   from expectedRank, especially if some columns are dominant in the data matrix. In this case
   The sampling size used to select colmuns is set to ``\\frac{expectedRank}{nbiter}``

 - number of iterations. Default to 2. Can be set to 1 ito get get
   a faster result at the expense of precision.

 ## Output

 The function does not return any output, instead it fills the fields U, S, Vt of the structure LowRankApproxMc.
 The fields can be retrieved by the function getIndex.
    If the rank extracted from the algorithm extractedRank we get as output :
  - U: a (d,extractedRank) matrix of orthonormalized left singular vector

  - S : a vector of singular values up to extractedRank

  - Vt : a (extractedRank , n) matrix , transposed of right singular vector.
    Nota the vector are not orthonormal.

    The svd approximation of data is:

    `` B = U * S * transpose(V) ``

    and verifies:

  `` E(|| data - B||^{2}) <=  1/(1-\\epsilon) * E(|| data - A_{k}|| ^{2}) + \\epsilon^{t} E(|| data ||^{2}) ``

 where ``A_{k}`` is the L2 best rank k approximation of data.
 B is of the form U * transpose(U) * data with U orthogonal made of left singular vectors.

 To get the corresponding result with right singular vector use a transpose Matrix as input.

 ## NOTA:

 **The Singular values are NOT sorted**

"""
function reduceToRank(datapb::LowRankApproxMc, rank::Int64; nbiter = 2)
    if rank > size(datapb.X)[1]
        error("reduceToRank reduction to rank greater than dimension")
    end
    maxdim = size(datapb.X)[2]
    datapb.currentIter = 1
    datapb.targetRank = rank
    # reset Epsil target in case of successive calls
    datapb.targetEpsil = 1.
    # we have rank target, 2 iterations should do
    if nbiter > 1
        datapb.samplingSize = floor(Int64,(1.05 * rank /nbiter))
        datapb.samplingSize = min(size(datapb.X)[2], datapb.samplingSize)
        @debug "LowRankApproxMc.reduceToRank setting samplingSize: " datapb.samplingSize
        maxiter = min(nbiter, 5)
        @warn "LowRankApproxMc.reduceToRank using maxiter " maxiter
        fstop2(datapb) = datapb.currentIter >= maxiter || datapb.currentRank >= datapb.targetRank ? true : false
        reduceIfast(datapb, fstop2)
    else
        # only one iter
        maxiter = 1
        datapb.samplingSize = floor(Int64,(1.10 * rank /nbiter))
        datapb.samplingSize = min(size(datapb.X)[2], datapb.samplingSize)
        @debug "LowRankApproxMc.reduceToRank setting samplingSize: " datapb.samplingSize
        fstop1(datapb) = datapb.currentIter >= maxiter ? true : false
        reduceIfast(datapb, fstop1)
    end
    # check we got rank
    if isnothing(datapb.U)
        @warn "reduceToRank  got not get left singular vectors"
    elseif size(datapb.U.value)[2] < rank
        @warn "reduceToRank  got only up to rank" size(datapb.U.value)[2]
    end
    # return rank approx
end

"""
# function reduceToEpsil(datapb::LowRankApproxMc, samplingSize::Int64, epsil::Float64;maxiter = 5)

This function computes an approximate svd up to a precision epsil.
    It computes U, S and Vt so that the matrix B = U*S*Vt verifies :
    `` E(|| data - B||_{2}) / E(|| data ||_{2}) <=  epsil ``
    It iterates until either the precision asked for is reached or *maxiter* iterations are done
        or rank reached is maximal.
## Args
------

-  samplingSize: the number of columns tentatively sampled by iteration.

- epsil : precision asked for. It is clear that a very small epsil will cause a full to be done
        with the corresponding costs. In this case possibly, a call to LAPACK.svd will be better.

- maxiter : the maximum number of iterations. default to 5

## Output:

The function does not return any output, instead it fills the fields U, S, Vt of the structure LowRankApproxMc.

The fields can be retrieved by the function getIndex.
If the rank extracted from the algorithm extractedRank we get as output :
- U: a (d,extractedRank) matrix of orthonormalized left singular vector

- S : a vector of singular values up to extractedRank

- Vt : a (extractedRank , n) matrix , transposed of right singular vector.
       Nota the vector are not orthonormal.

The svd approximation of data is: `` B = U * S * transpose(V) ``

and verifies:

`` E(|| data - B||^{2}) <=  1/(1-\\epsilon) * E(|| data - A_{k}|| ^{2}) + \\epsilon^{t} E(|| data ||^{2}) ``

where ``A_{k}`` is the L2 rank k approximation of data.
B is of the form U * transpose(U) * data with U orthogonal made of left singular vectors.

## NOTA:

 **The Singular values are NOT sorted**

"""
function reduceToEpsil(datapb::LowRankApproxMc, samplingSize::Int64, epsil::Float64; maxiter = 5)
    datapb.targetEpsil = epsil
    # reset targetRank in case of successive calls
    # rank is set to full dimension
    datapb.targetRank = min(size(datapb.X)[1], size(datapb.X)[2])
    if samplingSize <= 0
        throw("non strict positive samplingSize")
    end
    samplingSize = min(size(datapb.X)[2], samplingSize)
    datapb.samplingSize = samplingSize
    # epsil target
    # avoid doing one more iteration just for 1% of epsil
    fstop(datapb) = datapb.currentIter >= maxiter || datapb.residualErr <= epsil*(1.01) ||
            datapb.currentRank >= datapb.targetRank ? true : false
    reduceIfast(datapb, fstop)
end

# the following is faster as we do not ask for Vt

function getRangeApprox(data::Matrix{Float64}, rank::Int64)
    @debug " in getRangeApprox"
    svdpb = LowRankApproxMc(data)
    reduceToRank(svdpb, rank)
    U = something(svdpb.U)
    @debug "reduceToRank returned size" size(U)
    # compute all U columns and 0 columns of Vt
    U,S,Vt = LAPACK.gesvd!('A', 'O', U)
    @debug "lapack returned size U,S,Vt"  size(U) size(S) size(Vt)
    # return first k columns of U and first k values of diagonal of S
    U[:, (1:rank)],S[1:rank]
end


text:
using RemoteGtkREPL
using Test
using Serialization, Distributed, Sockets

global const c = Condition()

struct TestType end

function client_start_cb(port)

    @info "getting called $port"
    client = connect(port)

    @assert remotecall_fetch(sin, client, 0) ≈ 0.0
    @info 1
    for i=1:10
        @assert length(remotecall_fetch(rand, client, i)) == i
    end
    @info 2
    # eval x=3 on client and check if x is set
    remotecall_fetch(
        RemoteGtkREPL.eval_command_remotely,client,"x=3","RemoteGtkREPL"
    )
    @info 3
    while !remotecall_fetch(RemoteGtkREPL.isdone,client); sleep(0.01) end
    @info 3.5
    @assert remotecall_fetch(RemoteGtkREPL.run_task,client).result == ("3\n",nothing)
    @assert remotecall_fetch(include_string,client,Main,"@eval RemoteGtkREPL.x ") == 3
    @info 4.1
    remotecall_fetch(info,client,TestType())#this will fail and throw a warning

    notify(c)

    RemoteGtkREPL.remotecall_fetch(exit, client) #brutally kill the client
    close(client)

    @info "Done"
end

@testset "Server and Eval" begin

    port, server = RemoteGtkREPL.start_server()

    client = connect(port)
    #
    @test RemoteGtkREPL._eval_command_remotely("x=2", @__MODULE__) == ("2\n",nothing)
    @test x == 2
    #

    RemoteGtkREPL.process_message(data...) = println(data)
    RemoteGtkREPL.eval_command_remotely("x=3", string(@__MODULE__))

    @test RemoteGtkREPL.run_task().result == ("3\n",nothing)
    @test x == 3
    #
    @test remotecall_fetch(sin, client, 0) ≈ 0.0
    @test_throws MethodError throw(remotecall_fetch(sin, client, "pi"))

    p = joinpath(@__DIR__,"remote_startup.jl")

    #s = "tell application \"Terminal\" to do script \"julia -i \\\"$p\\\" $port 1\""
    #run(`osascript -e $s`)
    juliapath = ENV["_"] #what kind of variable name is this ?
    @async run(`$juliapath $p $port 1`)#this will call back client_start_cb and notify c
    wait(c)
    close(client)
    close(server)
end

ext: jl


text:
function execute_commands_part1(commands)
    dist, depth = 0, 0
    for (command, amount) in commands
        if command == "forward"
            dist += amount
        elseif command == "up"
            depth -= amount
        elseif command == "down"
            depth += amount
        else
            error("Should not happen. Command $command not recognized.")
        end
    end
    return dist, depth
end

function execute_commands_part2(commands)
    dist, depth, aim = 0, 0, 0
    for (command, amount) in commands
        if command == "down"
            aim += amount
        elseif command == "up"
            aim -= amount
        elseif command == "forward"
            dist += amount
            depth += aim * amount
        else
            error("Should not happen. Command $command not recognized.")
        end
    end
    return dist, depth
end

function day02a(commands)
    dist, depth = execute_commands_part1(commands)
    return dist * depth
end

function day02b(commands)
    dist, depth = execute_commands_part2(commands)
    return dist * depth
end

function day02(fpath)
    commands = parsecommands(fpath)
    return (day02a(commands), day02b(commands))
end

# println(parsecommands("test_input_02.txt"))
# println(day02a())
# println(day02b())


text:
include("sensors.jl")
using SatelliteToolbox

"""
q1 = SatelliteToolbox.Quaternion(cos(π/4), 1.0*sin(π/4), 0.0*sin(π/4), 0.0*sin(π/4))
q2 = SatelliteToolbox.Quaternion(cos(π/4), 0.0*sin(π/4), 0.0*sin(π/4), 1.0*sin(π/4))
dt = 5

ω = gyro_calcurate(q1, q2, dt)
println(ω)
"""

"""
mag_vecs = [29.8522 -4.0303 35.8413]

mag_vols = mag_sensor(mag_vecs, "vol")
println(mag_vols)

"""

sun_vecs = [-0.5 -0.5 1]
ss_out = sun_sensor(sun_vecs, "z-")
println(ss_out)

max_stars_repo_name: yui-project/ACSim

max_issues_repo_name: yui-project/ACSim

max_forks_repo_name: yui-project/ACSim


text:
using DiffRules

export Rules
export normalize

abstract type Rule end

struct Rules
    rules::Vector{Rule}
    Rules(rs::Vector{Rule}) = new(rs)
end
Rules(rs...) = Rules(collect(Rule, rs))
Base.:(==)(a::Rules, b::Rules) = a.rules == b.rules
Base.iterate(rs::Rules) = iterate(rs.rules)
Base.iterate(rs::Rules, state) = iterate(rs.rules, state)
Base.length(rs::Rules) = length(rs.rules)
Base.getindex(rs::Rules, ind) = getindex(rs.rules, ind)
Base.push!(rs::Rules, rule) = (push!(rs.rules, rule); rs)
Base.pushfirst!(rs::Rules, rule) = (push!(rs.rules, rule); rs)
Base.pop!(rs::Rules) = pop!(rs.rules)
Base.popfirst!(rs::Rules) = popfirst!(rs.rules)
Base.deleteat!(rs::Rules, ind) = (deleteat!(rs.rules, ind); rs)
Base.vcat(rss::Rules...) = Rules([(rs.rules for rs ∈ rss)...;])

normalize(rs::Rules) = Base.Fix2(normalize, rs)
function normalize(t::Term, rs::Rules)
    while true
        t = map(normalize(rs), t)  # FIXME: replace with `subexpressions`
        t′ = foldl(normalize, rs; init=t)
        t == t′ && return t
        t = t′
    end
end
normalize(::T, ::R) where {T,R<:Rule} = error("normalize undefined for rule type $R on term type $T")
normalize(t::Term, sets::Symbol...) = normalize(t, vcat(rules.(sets)...))
normalize(t::Term) = normalize(t, rules())

struct PatternRule <: Rule
    left::Term
    right::Term
    ps::Vector{Function}
end
PatternRule(l, r) = PatternRule(l, r, [])
PatternRule((l, r)::Pair) = PatternRule(l, r)
Base.convert(::Type{PatternRule}, p::Pair) = PatternRule(p)
Base.convert(::Type{Rule}, p::Pair) = convert(PatternRule, p)
Base.convert(::Type{Pair}, r::PatternRule) = r.left => r.right
Base.:(==)(a::PatternRule, b::PatternRule) = (a.left, a.right) == (b.left, b.right)
Base.hash(p::PatternRule, h::UInt) = hash((p.left, p.right), hash(PatternRule, h))
Base.map(f, r::PatternRule) = PatternRule(f(r.left), f(r.right), r.ps)
function normalize(t::Term, r::PatternRule)
    Θ = match(r.left, t) |> collect
    isempty(Θ) && return t

    σᵢ = findfirst(_preds_match(r.ps), Θ)
    σᵢ === nothing && return t
    σ = Θ[σᵢ]

    replace(r.right, σ)
end
_preds_match(ps, σ) = all(p -> p(σ), ps)
_preds_match(ps) = Base.Fix1(_preds_match, ps)

struct EvalRule <: Rule
    name
    f
end
EvalRule(f::Function) = EvalRule(f, f)
normalize(t::Term, r::EvalRule) = Term(normalize(get(t), r))
function normalize(ex::Expr, r::EvalRule)
    ex.head === :call || return ex
    f, args = ex.args[1], ex.args[2:end]

    match(Term(r.name), Term(f)) == zero(Match) && return ex

    if isvalid(Associative(f))
        if isvalid(Commutative(f))
            inds = findall(is_constant, args)
            if !isempty(inds)
                res = r.f(args[inds]...)
                deleteat!(args, inds)
                push!(args, res)
            end
        else
            _apply_associative!(r, args)
        end

        length(args) == 1 && return first(args)
        return Expr(:call, f, args...)
    end

    all(is_constant, args) || return ex
    r.f(args...)
end
normalize(x, ::EvalRule) = x
function _apply_associative!(r::EvalRule, args)
    i = firstindex(args)

    while i ≤ lastindex(args) - 1
        a, b = args[i:i+1]

        if is_constant(a) && is_constant(b)
            deleteat!(args, i)
            args[i] = r.f(a, b)
        else
            i += 1
        end
    end

    args
end

struct OrderRule <: Rule
    by::Function
end
normalize(t::Term, r::OrderRule) = Term(normalize(get(t), r))
function normalize(ex::Expr, r::OrderRule)
    ex.head === :call || return ex
    f = ex.args[1]
    isvalid(Commutative(f)) || return ex
    issorted(ex.args[2:end], by = r.by) && return ex
    args = sort(ex.args[2:end], by = r.by)
    Expr(ex.head, f, args...)
end
normalize(x, ::OrderRule) = x

struct DiffRule <: Rule end
normalize(t::Term, r::DiffRule) = Term(normalize(get(t), r))
function normalize(ex::Expr, r::DiffRule)
    fn = ex.args[1]
    fn === diff && length(ex.args) == 3 || return ex
    f, x = ex.args[[2, 3]]
    is_ground(f) && return :($zero($f))
    ex
end
normalize(x, ::DiffRule) = x


text:
# This file was generated, do not modify it. # hide
Q = MatrixDatabase(rand(dim, 30))
k = 5
I, D = searchbatch(G, Q, k)
D

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
using LightGraphs
using LightGraphs: add_edge!

in_channel = 3
out_channel = 5
N = 6
adj = [0. 1. 1. 0. 0. 0.;
       1. 0. 1. 0. 1. 0.;
       1. 1. 0. 1. 0. 1.;
       0. 0. 1. 0. 0. 0.;
       0. 1. 0. 0. 0. 0.;
       0. 0. 1. 0. 0. 0.]

ug = SimpleGraph(6)
add_edge!(ug, 1, 2); add_edge!(ug, 1, 3); add_edge!(ug, 2, 3)
add_edge!(ug, 3, 4); add_edge!(ug, 2, 5); add_edge!(ug, 3, 6)

@testset "Test support of LightGraphs for GCNConv layer" begin
    gc = GCNConv(ug, in_channel=>out_channel)
    @test size(gc.weight) == (out_channel, in_channel)
    @test size(gc.bias) == (out_channel, N)
    @test size(gc.norm) == (N, N)
end

@testset "Test support of LightGraphs for ChebConv layer" begin
    k = 4
    cc = ChebConv(ug, in_channel=>out_channel, k)
    @test size(cc.weight) == (out_channel, in_channel, k)
    @test size(cc.bias) == (out_channel, N)
    @test size(cc.L̃) == (N, N)
    @test cc.k == k
    @test cc.in_channel == in_channel
    @test cc.out_channel == out_channel
end

@testset "Test support of LightGraphs for GraphConv layer" begin
    gc = GraphConv(ug, in_channel=>out_channel)
    @test gc.adjlist == [[2,3], [1,3,5], [1,2,4,6], [3], [2], [3]]
    @test size(gc.weight1) == (out_channel, in_channel)
    @test size(gc.weight2) == (out_channel, in_channel)
    @test size(gc.bias) == (out_channel, N)
end

@testset "Test support of LightGraphs for GATConv layer" begin
    gat = GATConv(ug, in_channel=>out_channel)
    @test gat.adjlist == [[2,3], [1,3,5], [1,2,4,6], [3], [2], [3]]
    @test size(gat.weight) == (out_channel, in_channel)
    @test size(gat.bias) == (out_channel, N)
end

@testset "Test support of LightGraphs for GatedGraphConv layer" begin
    num_layers = 3
    ggc = GatedGraphConv(ug, out_channel, num_layers)
    @test ggc.adjlist == [[2,3], [1,3,5], [1,2,4,6], [3], [2], [3]]
    @test size(ggc.weight) == (out_channel, out_channel, num_layers)
end

@testset "Test support of LightGraphs for EdgeConv layer" begin
    ec = EdgeConv(ug, Dense(2*in_channel, out_channel))
    @test ec.adjlist == [[2,3], [1,3,5], [1,2,4,6], [3], [2], [3]]
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
"""

use Scenario struct to create reopt.jl model inputs
"""

struct REoptInputs <: AbstractInputs
    techs::Array{String, 1}
    pvtechs::Array{String, 1}
    gentechs::Array{String,1}
    elec_techs::Array{String, 1}
    techs_no_turndown::Array{String, 1}
    min_sizes::DenseAxisArray{Float64, 1}  # (techs)
    max_sizes::DenseAxisArray{Float64, 1}  # (techs)
    existing_sizes::DenseAxisArray{Float64, 1}  # (techs)
    cap_cost_slope::DenseAxisArray{Float64, 1}  # (techs)
    om_cost_per_kw::DenseAxisArray{Float64, 1}  # (techs)
    max_grid_export_kwh::Float64
    elec_load::ElectricLoad
    time_steps::UnitRange
    time_steps_with_grid::Array{Int, 1}
    time_steps_without_grid::Array{Int, 1}
    hours_per_timestep::Float64
    months::UnitRange
    production_factor::DenseAxisArray{Float64, 2}  # (techs, time_steps)
    levelization_factor::DenseAxisArray{Float64, 1}  # (techs)
    VoLL::Array{R, 1} where R<:Real #default set to 1 US dollar per kwh
    pwf_e::Float64
    pwf_om::Float64
    two_party_factor::Float64
    owner_tax_pct::Float64
    offtaker_tax_pct::Float64
    microgrid_premium_pct::Float64
    pvlocations::Array{Symbol, 1}
    maxsize_pv_locations::DenseAxisArray{Float64, 1}  # indexed on pvlocations
    pv_to_location::DenseAxisArray{Int, 2}  # (pv_techs, pvlocations)
    etariff::ElectricTariff
    ratchets::UnitRange
    techs_by_exportbin::DenseAxisArray{Array{String,1}}  # indexed on [:NEM, :WHL, :CUR]
    storage::Storage
    generator::Generator
    elecutil::ElectricUtility
    min_resil_timesteps::Int
    mg_tech_sizes_equal_grid_sizes::Bool
    node::Int
    export_bins_by_tech::Dict
end

"""
    REoptInputs(fp::String)

Use `fp` to load in JSON scenario:
```
function REoptInputs(fp::String)
    s = Scenario(JSON.parsefile(fp))
    REoptInputs(s)
end
```
Useful if you want to manually modify REoptInputs before solving the model.
"""
function REoptInputs(fp::String)
    s = Scenario(JSON.parsefile(fp))
    REoptInputs(s)
end

function REoptInputs(s::Scenario)

    time_steps = 1:length(s.electric_load.loads_kw)
    hours_per_timestep = 8760.0 / length(s.electric_load.loads_kw)
    techs, pvtechs, gentechs, pv_to_location, maxsize_pv_locations, pvlocations, production_factor,
        max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw  = setup_tech_inputs(s)
    elec_techs = techs  # only modeling electric loads/techs so far
    techs_no_turndown = pvtechs

    max_grid_export_kwh = sum(s.electric_load.loads_kw)

    months = 1:length(s.electric_tariff.monthly_demand_rates)

    techs_by_exportbin = DenseAxisArray([ techs, techs, techs], s.electric_tariff.export_bins)
    # TODO account for which techs have access to export bins (when we add more techs than PV)

    levelization_factor, pwf_e, pwf_om, two_party_factor = setup_present_worth_factors(s, techs, pvtechs)
    # the following hardcoded value for levelization_factor matches the public REopt API value
    # for test_with_cplex (test_time_of_export_rate) and makes the test values match.
    # the REopt code herein uses the Desktop method for levelization_factor, which is more accurate
    # (Desktop has non-linear degradation vs. linear degradation in API)
    # levelization_factor = DenseAxisArray([0.9539], techs)
    # levelization_factor = DenseAxisArray([0.9539, 1.0], techs)  # w/generator
    time_steps_with_grid, time_steps_without_grid, = setup_electric_utility_inputs(s)

    if any(pv.existing_kw > 0 for pv in s.pvs)
        adjust_load_profile(s, production_factor)
    end

    export_bins_by_tech = Dict{String, Array{Symbol, 1}}()
    for t in elec_techs
        export_bins_by_tech[t] = s.electric_tariff.export_bins
    end
    # TODO implement export bins by tech (rather than assuming that all techs share the export_bins)

    REoptInputs(
        techs,
        pvtechs,
        gentechs,
        elec_techs,
        techs_no_turndown,
        min_sizes,
        max_sizes,
        existing_sizes,
        cap_cost_slope,
        om_cost_per_kw,
        max_grid_export_kwh,
        s.electric_load,
        time_steps,
        time_steps_with_grid,
        time_steps_without_grid,
        hours_per_timestep,
        months,
        production_factor,
        levelization_factor,
        typeof(s.financial.VoLL) <: Array{<:Real, 1} ? s.financial.VoLL : fill(s.financial.VoLL, length(time_steps)),
        pwf_e,
        pwf_om,
        two_party_factor,
        s.financial.owner_tax_pct,
        s.financial.offtaker_tax_pct,
        s.financial.microgrid_premium_pct,
        pvlocations,
        maxsize_pv_locations,
        pv_to_location,
        s.electric_tariff,
        1:length(s.electric_tariff.tou_demand_ratchet_timesteps),  # ratchets
        techs_by_exportbin,
        s.storage,
        s.generator,
        s.electric_utility,
        s.site.min_resil_timesteps,
        s.site.mg_tech_sizes_equal_grid_sizes,
        s.site.node,
        export_bins_by_tech
    )
end

function setup_tech_inputs(s::Scenario)

    pvtechs = String[pv.name for pv in s.pvs]
    if length(Base.Set(pvtechs)) != length(pvtechs)
        error("PV names must be unique, got $(pvtechs)")
    end

    techs = copy(pvtechs)
    gentechs = String[]
    if s.generator.max_kw > 0
        push!(techs, "Generator")
        push!(gentechs, "Generator")
    end

    time_steps = 1:length(s.electric_load.loads_kw)

    # REoptInputs indexed on techs:
    max_sizes = DenseAxisArray{Float64}(undef, techs)
    min_sizes = DenseAxisArray{Float64}(undef, techs)
    existing_sizes = DenseAxisArray{Float64}(undef, techs)
    cap_cost_slope = DenseAxisArray{Float64}(undef, techs)
    om_cost_per_kw = DenseAxisArray{Float64}(undef, techs)
    production_factor = DenseAxisArray{Float64}(undef, techs, time_steps)

    # PV specific arrays
    pvlocations = [:roof, :ground, :both]
    pv_to_location = DenseAxisArray{Int}(undef, pvtechs, pvlocations)
    maxsize_pv_locations = DenseAxisArray([1.0e5, 1.0e5, 1.0e5], pvlocations)
    # default to large max size per location. Max size by roof, ground, both

    if !isempty(pvtechs)
        setup_pv_inputs(s, max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw, production_factor,
                        pvlocations, pv_to_location, maxsize_pv_locations)
    end

    if "Generator" in techs
        setup_gen_inputs(s, max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw, production_factor)
    end

    return techs, pvtechs, gentechs, pv_to_location, maxsize_pv_locations, pvlocations, production_factor,
    max_sizes, min_sizes, existing_sizes, cap_cost_slope, om_cost_per_kw
end

function setup_pv_inputs(s::Scenario, max_sizes, min_sizes,
    existing_sizes, cap_cost_slope, om_cost_per_kw, production_factor,
    pvlocations, pv_to_location, maxsize_pv_locations)

    time_steps = 1:length(s.electric_load.loads_kw)

    pv_roof_limited, pv_ground_limited, pv_space_limited = false, false, false
    roof_existing_pv_kw, ground_existing_pv_kw, both_existing_pv_kw = 0.0, 0.0, 0.0
    roof_max_kw, land_max_kw = 1.0e5, 1.0e5

    for pv in s.pvs
        production_factor[pv.name, :] = prodfactor(pv, s.site.latitude, s.site.longitude)
        for location in pvlocations
            if pv.location == location
                pv_to_location[pv.name, location] = 1
            else
                pv_to_location[pv.name, location] = 0
            end
        end

        beyond_existing_kw = pv.max_kw
        if pv.location == "both"
            both_existing_pv_kw += pv.existing_kw
            if !(s.site.roof_squarefeet === nothing) && !(s.site.land_acres === nothing)
                # don"t restrict unless both land_area and roof_area specified,
                # otherwise one of them is "unlimited"
                roof_max_kw = s.site.roof_squarefeet * pv.kw_per_square_foot
                land_max_kw = s.site.land_acres / pv.acres_per_kw
                beyond_existing_kw = min(roof_max_kw + land_max_kw, beyond_existing_kw)
                pv_space_limited = true
            end
        elseif pv.location == "roof"
            roof_existing_pv_kw += pv.existing_kw
            if !(s.site.roof_squarefeet === nothing)
                roof_max_kw = s.site.roof_squarefeet * pv.kw_per_square_foot
                beyond_existing_kw = min(roof_max_kw, beyond_existing_kw)
                pv_roof_limited = true
            end

        elseif pv.location == "ground"
            ground_existing_pv_kw += pv.existing_kw
            if !(s.site.land_acres === nothing)
                land_max_kw = s.site.land_acres / pv.acres_per_kw
                beyond_existing_kw = min(land_max_kw, beyond_existing_kw)
                pv_ground_limited = true
            end
        end

        existing_sizes[pv.name] = pv.existing_kw
        min_sizes[pv.name] = pv.existing_kw + pv.min_kw
        max_sizes[pv.name] = pv.existing_kw + beyond_existing_kw

        cap_cost_slope[pv.name] = effective_cost(;
            itc_basis=pv.cost_per_kw,
            replacement_cost=0.0,
            replacement_year=s.financial.analysis_years,
            discount_rate=s.financial.owner_discount_pct,
            tax_rate=s.financial.owner_tax_pct,
            itc=pv.total_itc_pct,
            macrs_schedule = pv.macrs_option_years == 7 ? s.financial.macrs_seven_year : s.financial.macrs_five_year,
            macrs_bonus_pct=pv.macrs_bonus_pct,
            macrs_itc_reduction = pv.macrs_itc_reduction,
            rebate_per_kw = pv.total_rebate_per_kw
        )

        om_cost_per_kw[pv.name] = pv.om_cost_per_kw
    end

    if pv_roof_limited
        maxsize_pv_locations[:roof] = float(roof_existing_pv_kw + roof_max_kw)
    end
    if pv_ground_limited
        maxsize_pv_locations[:ground] = float(ground_existing_pv_kw + land_max_kw)
    end
    if pv_space_limited
        maxsize_pv_locations[:both] = float(both_existing_pv_kw + roof_max_kw + land_max_kw)
    end

    return nothing
end

function setup_gen_inputs(s::Scenario, max_sizes, min_sizes, existing_sizes,
    cap_cost_slope, om_cost_per_kw, production_factor)

    max_sizes["Generator"] = s.generator.max_kw
    min_sizes["Generator"] = s.generator.existing_kw + s.generator.min_kw
    existing_sizes["Generator"] = s.generator.existing_kw
    cap_cost_slope["Generator"] = s.generator.cost_per_kw
    om_cost_per_kw["Generator"] = s.generator.om_cost_per_kw
    production_factor["Generator", :] = prodfactor(s.generator)
    return nothing
end

function setup_present_worth_factors(s::Scenario, techs::Array{String, 1}, pvtechs::Array{String, 1})

    lvl_factor = DenseAxisArray{Float64}(undef, techs)
    for (i, tech) in enumerate(pvtechs)
        lvl_factor[tech] = levelization_factor(
            s.financial.analysis_years,
            s.financial.elec_cost_escalation_pct,
            s.financial.offtaker_discount_pct,
            s.pvs[i].degradation_pct  # TODO generalize for any tech (not just pvs)
        )
    end
    if "Generator" in techs
        lvl_factor["Generator"] = 1
    end

    pwf_e = annuity(
        s.financial.analysis_years,
        s.financial.elec_cost_escalation_pct,
        s.financial.offtaker_discount_pct
    )

    pwf_om = annuity(
        s.financial.analysis_years,
        s.financial.om_cost_escalation_pct,
        s.financial.owner_discount_pct
    )

    if s.financial.two_party_ownership
        pwf_offtaker = annuity(s.financial.analysis_years, 0.0, s.financial.offtaker_discount_pct)
        pwf_owner = annuity(s.financial.analysis_years, 0.0, s.financial.owner_discount_pct)
        two_party_factor = (pwf_offtaker * (1 - s.financial.offtaker_tax_pct)) /
                           (pwf_owner * (1 - s.financial.owner_tax_pct))
    else
        two_party_factor = 1.0
    end

    return lvl_factor, pwf_e, pwf_om, two_party_factor
end

function setup_electric_utility_inputs(s::AbstractScenario)
    if s.electric_utility.outage_end_timestep > 0 &&
            s.electric_utility.outage_end_timestep > s.electric_utility.outage_start_timestep
        time_steps_without_grid = Int[i for i in range(s.electric_utility.outage_start_timestep,
                                                    stop=s.electric_utility.outage_end_timestep)]
        if s.electric_utility.outage_start_timestep > 1
            time_steps_with_grid = append!(
                Int[i for i in range(1, stop=s.electric_utility.outage_start_timestep - 1)],
                Int[i for i in range(s.electric_utility.outage_end_timestep + 1,
                                     stop=length(s.electric_load.loads_kw))]
            )
        else
            time_steps_with_grid = Int[i for i in range(s.electric_utility.outage_end_timestep + 1,
                                       stop=length(s.electric_load.loads_kw))]
        end
    else
        time_steps_without_grid = Int[]
        time_steps_with_grid = Int[i for i in range(1, stop=length(s.electric_load.loads_kw))]
    end
    return time_steps_with_grid, time_steps_without_grid
end

function adjust_load_profile(s::Scenario, production_factor::DenseAxisArray)
    if s.electric_load.loads_kw_is_net
        for pv in s.pvs if pv.existing_kw > 0
            s.electric_load.loads_kw .+= pv.existing_kw * production_factor[pv.name, :].data
        end end
    end

    if s.electric_load.critical_loads_kw_is_net
        for pv in s.pvs if pv.existing_kw > 0
            s.electric_load.critical_loads_kw .+= pv.existing_kw * production_factor[pv.name, :].data
        end end
    end
end

max_stars_repo_name: NLaws/REoptLite

max_issues_repo_name: NLaws/REoptLite

max_forks_repo_name: NLaws/REoptLite


text:
function simplify(c::Char)
    return c
end

function ==(a::HiddenMarkovModel{Tv, Ti}, b::HiddenMarkovModel{Tv, Ti}) where {Tv, Ti}
    res = true
    for fname in fieldnames(HiddenMarkovModel{Tv, Ti})
        res &= getfield(a, fname) == getfield(b, fname)
    end
    return res
end

function write(io::IO, obj::HiddenMarkovModel{Tv, Ti}) where {Tv, Ti}
    nbit = 0
    nbit += write(io, sizeof(Tv))
    if Ti isa Unsigned
        nbit += write(io, +1)
    else
        nbit += write(io, -1)
    end
    nbit += write(io, sizeof(Ti))
    nbit += write(io, obj.aca)
    for p in obj.pos
        nbit += write(io, p, " ")
    end
    nbit += write(io, "\n")
    nbit += write(io, obj.hpr)
    nbit += write(io, obj.h2h)
    for vp in obj.h2v
        nbit += write(io, length(vp))
        for (v, p) in vp
            nbit += write(io, v)
            nbit += write(io, p)
        end
    end
    nbit += write(io, obj.MAX)
    nbit += write(io, obj.sumhpr)
    nbit += write(io, obj.sumh2hdim1)
    nbit += write(io, obj.sumh2vdimv)
    return nbit
end

function nbit2type(nbit, types)
    filter(x -> sizeof(x) == nbit, types)[1]
end

function read(io::IO, obj::Type{HiddenMarkovModel})
    Tv = nbit2type(read(io, Int), [Float32, Float64])
    sign = read(io, Int)
    Ti = Int
    if sign > 0
        Ti = nbit2type(read(io, Int), [UInt8, UInt16, UInt32, UInt64])
    else
        Ti = nbit2type(read(io, Int), [Int8, Int16, Int32, Int64])
    end
    aca = read(io, AhoCorasickAutomaton)
    pos = split(readline(io))
    npos = length(pos)
    hpr = reinterpret(Tv, read(io, npos * sizeof(Tv)))
    h2h = reshape(reinterpret(Tv, read(io, npos * npos * sizeof(Tv))), (npos, npos))
    h2v = Vector{Dict{Int, Tv}}(undef, npos)
    for i = 1:npos
        vps = Dict{Int, Tv}()
        nv = read(io, Int)
        for j = 1:nv
            v = read(io, Int)
            p = read(io, Tv)
            vps[v] = p
        end
        h2v[i] = vps
    end
    MAX    = read(io, Tv)
    sumhpr = read(io, Tv)
    sumh2hdim1 = reshape(reinterpret(Tv, read(io, npos * sizeof(Tv))), 1, npos)
    sumh2vdimv = reinterpret(Tv, read(io, npos * sizeof(Tv)))
    return HiddenMarkovModel{Tv, Ti}(aca, pos, hpr, h2h, h2v, MAX, sumhpr, sumh2hdim1, sumh2vdimv)
end
function h2v(obj::HiddenMarkovModel{Tv, Ti}, h::String, v::String) where Tv where Ti
    ih = findfirst(isequal(h), obj.pos)
    iv = get(obj.aca, v, Ti(0))
    return get(obj.h2v[ih], iv, obj.alpha[ih])
end

function h2h(obj::HiddenMarkovModel{Tv, Ti}, pos1::String, pos2::String) where Tv where Ti
    i1 = findfirst(isequal(pos1), obj.pos)
    i2 = findfirst(isequal(pos2), obj.pos)
    return obj.h2h[i2, i1]
end

function norm!(hpr::Vector{Tv}) where Tv
    tot = sum(hpr)
    @assert !isinf(tot)
    hpr .= -log.(hpr ./ tot)
    return tot
end

function norm!(h2h::Matrix{Tv}) where Tv
    sum1 = sum(h2h; dims = 1)
    h2h .= -log.(h2h ./ sum1)
    return sum1
end

function norm!(h2v::Vector{Dict{Int, Tv}}) where Tv
    sumv = Vector{Tv}(undef, length(h2v))
    for (ih, vs) in enumerate(h2v)
        tot = sum(values(vs))
        sumv[ih] = tot
        for v in keys(vs)
            vs[v] = -log(vs[v]) + log(tot)
        end
    end
    return sumv
end

function denorm!(hpr::Vector{Tv}, sumhpr::Tv) where Tv
    return hpr .= exp.(.-hpr) .* sumhpr
end

function denorm!(h2h::Matrix{Tv}, sumh2hdim1::Matrix{Tv}) where Tv
    return h2h .= exp.(.-h2h) .* sumh2hdim1
end

function denorm!(h2v::Vector{Dict{Int, Tv}}, sumh2vdimv::Vector{Tv}) where Tv
    for (h, vs) in enumerate(h2v)
        for (v, p) in vs
            vs[v] = exp(-p) * sumh2vdimv[h]
        end
    end
    return h2v
end

function stat(obj::HiddenMarkovModel{Tv, Ti}) where {Tv, Ti}
    res = Vector()
    dict = map(first, sort!(collect(obj.aca); by = last))
    hpr = copy(obj.hpr)
    denorm!(hpr, obj.sumhpr)
    for i = 1:length(obj.pos)
        pos = obj.pos[i]
        nhead = round(Int, hpr[i])
        tot = round(Int, obj.sumh2vdimv[i])
        vs  = sort!(collect(obj.h2v[i]); by = last, rev = true)
        vs = map(x -> dict[x[1]], vs)
        push!(res, (pos = pos, nhead = nhead, tot = tot, vs = vs))
    end
    return res
end


text:
########################################
## File Name: costate_transition_test.jl
## Author: Haruki Nishimura (hnishimura@stanford.edu)
## Date Created: 2020/05/14
## Description: Test code for src/costate_transition_models.jl
########################################

using LinearAlgebra

@testset "Costate Transition Test" begin

# Localization Task
ukfModel = UKFPosRange();
transPos = TransModel_Pos();
bPos = BelMvNormal(0.,zeros(2),100*Matrix(1.0I, 2, 2));
bVecPos = [bPos for ii = 1:100];
p_robot = PhysPos(0.,[5.,7.]);
sV = AugState(p_robot,bVecPos);
costModel = CostPosRangeLocalization();
cosP = PhysPos(0.,termCost_grad_p(costModel,sV))
cosBV = [BelMvNormal(0.0,params) for params in termCost_grad_b(costModel,sV)];
coS = AugState(cosP,cosBV);

cosP_new = cotrans(transPos,costModel,coS,sV,PosControl(0.,zeros(2)),0.01)
@test cosP_new.t == -0.01
@test cosP_new.p.pos == zeros(2)
@test all([isapprox(cosP_new.b[ii].params, [0.0, 0.0, 8.53973, 0.0, 0.0, 8.53973], atol=1e-4) for ii = 1:100])

cosP_new_ukf = cotrans(ukfModel,costModel,coS,sV,ones(100),0.2,Matrix(1.0I, 2, 2))
@test cosP_new_ukf.t == 0.0
@test isapprox(cosP_new_ukf.p.pos, [-2002.46, -2992.51], rtol=1e-5)
@test all([isapprox(cosP_new_ukf.b[ii].params, [20.0246, 29.9251, 8.87919, 0.0, -2.44889, 8.39393], rtol=1e-5) for ii = 1:100])

# Manipulation Task
xMP = PhysManipulate2D(0.,[0.1,0.2],0.5,[0.3,0.4],0.6,3.,5.,[10.,11.],0.7);
uMP = MControl2D(0.,6.,7.,8.);
bMP = BelMvNormal(xMP.t,vec(xMP),100*Matrix(1.0I, 11, 11));

QMP = [Matrix(1.0I, 6, 6) zeros(6,5); zeros(5,6) 0.01*Matrix(1.0I, 5, 5)];
RMP = Matrix(1.0I, 9, 9);
CsMP = QMP;
CuMP = Matrix(1.0I, 3, 3);

costModel = CostManipulate2D()
cob = BelMvNormal(bMP.t,termCost_grad(costModel,bMP,CsMP))
transModelC = CPredictManipulate2D();

cob_new = cotrans(transModelC,costModel,cob,bMP,uMP,0.01,QMP,CsMP)
@test cob_new.t == -0.01
@test isapprox(cob_new.params[1:11], [0.101, 0.202, -2.81458, 0.3033, 0.405067, 0.579584, 0.180939, 0.063501, 0.10492, 0.100754, -0.66043], rtol=1e-4)

transModelD = DUpdateManipulate2D();
cob_new_u = cotrans(transModelD,costModel,cob,bMP,uMP,ones(9),RMP)
@test cob_new_u.t == 0.0
@test isapprox(cob_new_u.params[1:11], [0.0226903, -0.00757826, -3.87298, 0.164995, -0.0386729, 0.118948, -0.604347, 0.463712, 0.924217, 0.857975, 0.186299], rtol=1e-5)

end


text:
using ComputationalHomology
using LinearAlgebra
using SparseArrays
using Test
import Random

include("cells.jl")
include("chain.jl")
include("simplicial.jl")
include("constructions.jl")
include("homology.jl")
include("filtration.jl")
include("persistence.jl")
include("landscape.jl")


text:
include("testfunctions.jl")
using DiferentialEvolution
using Test

@testset "DiferentialEvolution.jl" begin

    @testset "Classic Diferential Evolution" begin

        @testset "2-Dimensioanl tests" begin

        end

        @testset "5-Dimensional tests" begin

            AckleyT = DE_classic(Ackley,1000,300, [-32.768,-32.768,-32.768,-32.768,-32.768], [32.768,32.768,32.768,32.768,32.768],0.9,0.5)
            SchwefelT = DE_classic(Schwefel,1000,300, [-500,-500,-500,-500,-500], [500,500,500,500,500],0.9,0.5)
            GriewankT = DE_classic(Griewank,1000,300, [-600,-600,-600,-600,-600], [600,600,600,600,600],0.9,0.5)
            RosenbrockT = DE_classic(Rosenbrock,1000,300, [-5,-5,-5,-5,-5], [10,10,10,10,10],0.9,0.5)
            MichalewiczT = DE_classic(Michalewicz,1000,300, [0,0,0,0,0], [pi,pi,pi,pi,pi],0.9,0.5)

            #Ackley     critic value = 0   ///  parameters = [0,0,0,0,0]
            @test isapprox(AckleyT[1],0, atol=1e-8)
            #Schwefel   critic value = -418.9829 d   ///  parameters = [420.9687,420.9687,420.9687,420.9687,420.9687]
            d = 5
            @test isapprox(SchwefelT[1],-418.9829*d, atol=1e-4)
            #Griewank   critic value = 0   ///  parameters = [0,0,0,0,0]
            @test isapprox(GriewankT[1],0, atol=1e-8)
            #Rosenbrock critic value = 0   ///  parameters = [1,1,1,1,1]
            @test isapprox(GriewankT[1],0, atol=1e-8)
            #Rosenbrock critic value = 0   ///  parameters = [1,1,1,1,1]
            @test isapprox(MichalewiczT[1],-4.687658, atol=1e-6)
        end

    end

    @testset "Dither Diferential Evolution" begin

    end

    @testset "Jither Diferential Evolution" begin

    end

end


text:
"""
    fit!(estimator, num_iterations, X, y[, test...]; [verbosity = 1, is_row_major = false])
    fit!(estimator, X, y[, test...]; [verbosity = 1, is_row_major = false])

Fit the `estimator` with features data `X` and label `y` using the X-y pairs in `test` as
validation sets.

Return a dictionary with an entry for each validation set. Each entry of the dictionary is another
dictionary with an entry for each validation metric in the `estimator`. Each of these entries is an
array that holds the validation metric's value at each iteration.

# Arguments
* `estimator::LGBMEstimator`: the estimator to be fit.
* `num_iterations::Int`: OPTIONAL -- defaults to estimator.num_iterations if not provided
* `X::Matrix{TX<:Real}`: the features data.
* `y::Vector{Ty<:Real}`: the labels.
* `test::Tuple{Matrix{TX},Vector{Ty}}...`: optionally contains one or more tuples of X-y pairs of
    the same types as `X` and `y` that should be used as validation sets.
* `verbosity::Integer`: keyword argument that controls LightGBM's verbosity. `< 0` for fatal logs
    only, `0` includes warning logs, `1` includes info logs, and `> 1` includes debug logs.
* `is_row_major::Bool`: keyword argument that indicates whether or not `X` is row-major. `true`
    indicates that it is row-major, `false` indicates that it is column-major (Julia's default).
* `weights::Vector{Tw<:Real}`: the training weights.
* `init_score::Vector{Ti<:Real}`: the init scores.
"""
function fit!(
    estimator::LGBMEstimator, num_iterations::Int, X::Matrix{TX}, y::Vector{Ty}, test::Tuple{Matrix{TX},Vector{Ty}}...;
    verbosity::Integer = 1,
    is_row_major = false,
    weights::Vector{Tw} = Float32[],
    init_score::Vector{Ti} = Float64[],
) where {TX<:Real,Ty<:Real,Tw<:Real,Ti<:Real}

    start_time = now()

    log_debug(verbosity, "Started creating LGBM training dataset\n")
    ds_parameters = stringifyparams(estimator, DATASETPARAMS)
    train_ds = LGBM_DatasetCreateFromMat(X, ds_parameters, is_row_major)
    LGBM_DatasetSetField(train_ds, "label", y)
    if length(weights) > 0
        LGBM_DatasetSetField(train_ds, "weight", weights)
    end
    if length(init_score) > 0
        LGBM_DatasetSetField(train_ds, "init_score", init_score)
    end

    log_debug(verbosity, "Started creating LGBM booster\n")
    bst_parameters = stringifyparams(estimator, BOOSTERPARAMS) * " verbosity=$verbosity"
    estimator.booster = LGBM_BoosterCreate(train_ds, bst_parameters)

    n_tests = length(test)
    tests_names = Array{String}(undef,n_tests)

    if n_tests > 0
        log_debug(verbosity, "Started creating LGBM test datasets\n")
        @inbounds for (test_idx, test_entry) in enumerate(test)
            tests_names[test_idx] = "test_$(test_idx)"
            test_ds = LGBM_DatasetCreateFromMat(test_entry[1], ds_parameters, train_ds, is_row_major)
            LGBM_DatasetSetField(test_ds, "label", test_entry[2])
            LGBM_BoosterAddValidData(estimator.booster, test_ds)
        end
    end

    log_debug(verbosity, "Started training...\n")
    results = train!(estimator, tests_names, verbosity, start_time)

    return results
end
# Old signature, pass through args
function fit!(
    estimator::LGBMEstimator,
    X::Matrix{TX},
    y::Vector{Ty},
    test::Tuple{Matrix{TX},Vector{Ty}}...;
    kwargs...
)  where {TX<:Real,Ty<:Real}
    return fit!(estimator, estimator.num_iterations, X, y, test...; kwargs...)
end

function train!(
    estimator::LGBMEstimator, num_iterations::Int, tests_names::Vector{String}, verbosity::Integer, start_time::DateTime
)
    results = Dict{String,Dict{String,Vector{Float64}}}()
    n_tests = length(tests_names)
    metrics = LGBM_BoosterGetEvalNames(estimator.booster)
    n_metrics = length(metrics)
    bigger_is_better = [ifelse(in(metric, MAXIMIZE_METRICS), 1., -1.) for metric in metrics]
    best_score = fill(-Inf, (n_metrics, n_tests))
    best_iter = fill(1, (n_metrics, n_tests))

    start_iter = get_iter_number(estimator) + 1
    end_iter = start_iter + num_iterations - 1

    metrics_idx_sequence = (((start_iter:end_iter) .- 1) .% estimator.metric_freq) .== 0
    total_metrics_evals = sum(metrics_idx_sequence)
    metric_idx = 0

    for (idx, iter) in enumerate(start_iter:end_iter)

        is_finished = LGBM_BoosterUpdateOneIter(estimator.booster)

        log_debug(verbosity, Dates.CompoundPeriod(now() - start_time), " elapsed, finished iteration ", iter, "\n")

        metric_idx = metric_idx + metrics_idx_sequence[idx]

        if is_finished == 0
            is_finished = eval_metrics!(
                results, estimator, tests_names, iter, metric_idx, total_metrics_evals, n_metrics,
                verbosity, bigger_is_better, best_score, best_iter, metrics,
            )
        else
            shrinkresults!(results, metric_idx)
            log_info(verbosity, "Stopped training because there are no more leaves that meet the ",
                     "split requirements.")
        end
        if is_finished == 1
            break
        end
    end

    # save the model in serialised form, in case we should be deepcopied or serialised elsewhere
    estimator.model = LGBM_BoosterSaveModelToString(estimator.booster, 0, 0)

    return results
end
# Old signature, pass through args
function train!(
    estimator::LGBMEstimator, tests_names::Vector{String}, verbosity::Integer, start_time::DateTime
)
    return train!(estimator, estimator.num_iterations, tests_names, verbosity, start_time)
end

function eval_metrics!(
    results::Dict{String,Dict{String,Vector{Float64}}},
    estimator::LGBMEstimator,
    tests_names::Vector{String},
    iter::Integer,
    metrics_store_idx::Integer,
    metrics_stored_count::Integer,
    n_metrics::Integer,
    verbosity::Integer,
    bigger_is_better::Vector{Float64},
    best_score::Matrix{Float64},
    best_iter::Matrix{Int},
    metrics::Vector{String},
)

    if (iter - 1) % estimator.metric_freq == 0
        if estimator.is_training_metric
            scores = LGBM_BoosterGetEval(estimator.booster, 0)
            store_scores!(results, estimator, metrics_store_idx, metrics_stored_count, "training", scores, metrics)
            print_scores(estimator, iter, "training", n_metrics, scores, metrics, verbosity)
        end
    end

    # Metrics for test sets
    if (iter - 1) % estimator.metric_freq == 0 || estimator.early_stopping_round > 0
        for (test_idx, test_name) in enumerate(tests_names)
            scores = LGBM_BoosterGetEval(estimator.booster, test_idx)
            # Check if progress should be stored and/or printed
            if (iter - 1) % estimator.metric_freq == 0
                store_scores!(results, estimator, metrics_store_idx, metrics_stored_count, test_name, scores, metrics)
                print_scores(estimator, iter, test_name, n_metrics, scores, metrics, verbosity)
            end

            # Check if early stopping is called for
            @inbounds for metric_idx in eachindex(metrics)
                maximize_score = bigger_is_better[metric_idx] * scores[metric_idx]
                if maximize_score > best_score[metric_idx, test_idx]
                    best_score[metric_idx, test_idx] = maximize_score
                    best_iter[metric_idx, test_idx] = iter
                elseif iter - best_iter[metric_idx, test_idx] >= estimator.early_stopping_round
                    # This will shrink it up to the current stored metric
                    shrinkresults!(results, metrics_store_idx)
                    log_info(verbosity, "Early stopping at iteration ", iter,
                             ", the best iteration round is ", best_iter[metric_idx, test_idx], "\n")
                    return 1
                end
            end
        end
    end

    return 0
end

function store_scores!(
    results::Dict{String,Dict{String,Vector{Float64}}},
    estimator::LGBMEstimator,
    store_idx::Integer,
    num_evals::Integer,
    evalname::String,
    scores::Vector{Cdouble},
    metrics::Vector{String},
)

    for (metric_idx, metric_name) in enumerate(metrics)
        if !haskey(results, evalname)
            results[evalname] = Dict{String,Vector{Float64}}()
            results[evalname][metric_name] = Array{Float64}(undef,num_evals)
        elseif !haskey(results[evalname], metric_name)
            results[evalname][metric_name] = Array{Float64}(undef,num_evals)
        end
        results[evalname][metric_name][store_idx] = scores[metric_idx]
    end

    return nothing
end

function merge_scores(
    old_scores::Dict{String,Dict{String,Vector{Float64}}},
    additional_scores::Dict{String,Dict{String,Vector{Float64}}},
)

    if keys(old_scores) != keys(additional_scores)
        throw(ErrorException("Tried to merge metrics with different data sets:\n    a=> $(keys(old_scores))\n    b=> $(keys(additional_scores))"))
    end

    new_scores = Dict{String,Dict{String,Vector{Float64}}}()
    for (key, oldvals) in old_scores
        new_scores[key] = merge_scores(old_scores[key], additional_scores[key])
    end

    return new_scores
end
function merge_scores(
    old_scores::Dict{String,Vector{Float64}},
    additional_scores::Dict{String,Vector{Float64}},
)

    if keys(old_scores) != keys(additional_scores)
        throw(ErrorException("Tried to merge metrics with different metrics:\n    a=> $(keys(old_scores))\n    b=> $(keys(additional_scores))"))
    end

    new_scores = Dict{String,Vector{Float64}}()
    for (key, oldvals) in old_scores
        new_scores[key] = cat(oldvals, additional_scores[key]; dims=1)
    end

    return new_scores
end

function print_scores(estimator::LGBMEstimator, iter::Integer, name::String, n_metrics::Integer,
                      scores::Vector{Cdouble}, metrics::Vector{String}, verbosity::Integer)
    log_info(verbosity, "Iteration: ", iter, ", ", name, "'s ")
    for (metric_idx, metric_name) in enumerate(metrics)
        log_info(verbosity, metric_name, ": ", scores[metric_idx])
        metric_idx < n_metrics && log_info(verbosity, ", ")
    end
    log_info(verbosity, "\n")
end

function stringifyparams(estimator::LGBMEstimator, params::Vector{Symbol})
    paramstring = ""
    n_params = length(params)
    valid_names = fieldnames(typeof(estimator))
    for (param_idx, param_name) in enumerate(params)
        if in(param_name, valid_names)
            param_value = getfield(estimator, param_name)
            if !isempty(param_value)
                # Convert parameters that contain indices to C's zero-based indices.
                if in(param_name, INDEXPARAMS)
                    param_value -= 1
                end

                if typeof(param_value) <: Array
                    n_entries = length(param_value)
                    if n_entries >= 1
                        paramstring = string(paramstring, param_name, "=", param_value[1])
                        for entry_idx in 2:n_entries
                            paramstring = string(paramstring, ",", param_value[entry_idx])
                        end
                        paramstring = string(paramstring, " ")
                    end
                else
                    paramstring = string(paramstring, param_name, "=", param_value, " ")
                end
            end
        end
    end
    return paramstring[1:end - 1]
end


text:
using PlutoUI
using Test
import AbstractPlutoDingetjes
using HypertextLiteral
import ColorTypes: RGB, N0f8, Colorant
import Logging

# has to be outside of the begin block for julia 1.0 compat
struct Uhm end

@testset "DisplayTricks" begin

    x = "asdf"

    # Show

    @test repr(MIME"derp1"(), Show(MIME"derp1"(), x)) == codeunits(x)

    Base.istextmime(::MIME"derp2") = true

    @test repr(MIME"derp2"(), Show(MIME"derp2"(), codeunits(x))) == x

    s = Show(MIME"cool1"(), x)

    @test Base.showable(MIME"text/plain"(), x)
    @test !Base.showable(MIME"cool1"(), x)
    # everything is plaintext showable
    @test Base.showable(MIME"text/plain"(), s)
    @test Base.showable(MIME"cool1"(), s)

    # AsMIME

    mimer = as_mime(MIME"cool2"())
    am = mimer(s)

    @test Base.showable(MIME"cool1"(), s)
    @test !Base.showable(MIME"cool2"(), s)
    @test Base.showable(MIME"text/plain"(), s)

    @test !Base.showable(MIME"cool1"(), am)
    @test Base.showable(MIME"cool2"(), am)
    @test Base.showable(MIME"text/plain"(), am)

    # WithIOContext

    @test WithIOContext(s, :a => 1, :b => 2) ==
          WithIOContext(s, :a => 1; b = 2) ==
          WithIOContext(s, a = 1; b = 2) ==
          WithIOContext(s; a = 1, b = 2)

    @test repr(MIME"text/plain"(), WithIOContext(1.345234523452, compact = true)) == "1.34523"
    @test repr(MIME"text/plain"(), WithIOContext(1.345234523452, compact = false)) == "1.345234523452"

    m = MIME"hello"()
    Base.istextmime(::MIME"hello") = true

    Base.show(io::IO, ::MIME"hello", ::Uhm) = print(io, get(io, :prop, missing))

    u = Uhm()

    @test repr(m, WithIOContext(u, prop = 1)) == "1"
end

default(x) = AbstractPlutoDingetjes.Bonds.initial_value(x)

@testset "Public API" begin
    el = Button()
    el = Button("asdf")
    @test default(el) == "asdf"
    el = LabelButton("asdf")
    @test default(el) == "asdf"
    el = CounterButton()
    @test default(el) == 0
    el = CounterButton("asdf")
    @test default(el) == 0

    el = CheckBox()
    @test default(el) == false
    el = CheckBox(true)
    @test default(el) == true
    el = CheckBox(default = true)
    @test default(el) == true

    el = Clock(2.0, max_value = 123)
    @test default(el) == 1
    Clock(3.0)
    Clock(3.0, true)
    Clock(3.0, true, true)
    Clock(3.0, true, true, 5)
    Clock()

    el = ColorStringPicker()
    @test default(el) == "#000000"
    el = ColorStringPicker("#0f0f0f")
    @test default(el) == "#0f0f0f"
    el = ColorStringPicker(default = "#0f0f0f")

    gray = RGB{N0f8}(0.5, 0.5, 0.5)
    el = ColorPicker()
    @test default(el) == RGB{N0f8}(0.0, 0.0, 0.0)
    el = ColorPicker(gray)
    @test default(el) == gray
    el = ColorPicker(default = gray)
    @test default(el) == gray

    el = DateField()
    @test default(el) === nothing
    el = TimeField()
    @test default(el) === "" # ugh

    el = FilePicker()
    @test default(el) === nothing

    @testset "MultiSelect" for f in [MultiSelect, MultiCheckBox]
        el = f(["asdf", "x"])
        @test default(el) == []
        el = f(["asdf"])
        @test default(el) == []
        el = f(["sin" => "asdf"])
        @test default(el) == []
        el = f(["sin" => "asdf"]; default = ["sin"])
        @test default(el) == ["sin"]

        el = f([sin, cos, tan, sqrt])
        @test default(el) |> isempty
        @test default(el) isa Vector{Function}
        el = f([
            cos => "cosine function",
            sin => "sine function",
        ])
        @test default(el) |> isempty
        @test default(el) isa Vector{Function}
    end

    el = MultiCheckBox(
        ["🐱" => "🐝", "🐵" => "🦝", "🐱" => "🐿"];
        default=["🐱", "🐱"]
    )
    @test default(el) == ["🐱", "🐱"]
    el = MultiCheckBox(
        ["🐱" => "🐝", "🐵" => "🦝", "🐱" => "🐿"];
        default=["🐱"]
    )
    @test default(el) == ["🐱"]

    el = Select(["asdf", "x"])
    @test default(el) == "asdf"
    el = Select(["asdf"])
    @test default(el) == "asdf"
    el = Select([sin => "asdf"])
    @test default(el) == sin
    el = Select([sin => "asdf", cos => "cos"]; default = cos)
    @test default(el) == cos

    el = Radio(["asdf", "x"])
    @test default(el) == nothing
    el = Radio(["asdf"])
    @test default(el) == nothing
    el = Radio(["sin" => "asdf"])
    @test default(el) == nothing
    el = Radio(["sin" => "asdf", "cos" => "cos"]; default = "cos")
    @test default(el) == "cos"

    el = Slider(0.0:π:20)
    @test default(el) == 0
    el = Slider(0.0:π:20; show_value = true)
    el = Slider(0.0:π:20; default = π)
    @test default(el) == Float64(π) # should have been converted to Float64 because our range has been
    el = Slider(1:1//3:20; default = 7 // 3)
    @test default(el) === 7 // 3
    el = Slider(1:1//3:20)
    @test default(el) === 1 // 1
    el = Slider([sin, cos, tan])
    @test default(el) == sin
    el = Slider([sin, cos, tan]; default = tan)
    @test default(el) == tan

    el = Scrubbable(60)
    @test default(el) === 60
    el = Scrubbable(60.0)
    @test default(el) === 60.0

    el = Scrubbable(0.0:π:20)
    @test 0.0 < default(el) < 20
    el = Scrubbable(0.0:π:20; default = π)
    # @test default(el) == Float64(π) # should have been converted to Float64 because our range has been
    el = Scrubbable(1:1//3:20; default = 7 // 3)
    @test default(el) === 7 // 3
    el = Scrubbable(1:1//3:20)
    @test default(el) isa Rational

    @testset "Rounding default value: $f" for f in [Slider, NumberField]
        el = f(1:10; default = 5.2)
        @test default(el) == 5
        el = f(1:.5:10; default = 5.4)
        @test default(el) == 5.5
        if f !== NumberField
            el = f([60,10,-80]; default = 5.4)
            @test default(el) == 10
            el = Slider([sin, cos]; default = tan)
            @test default(el) == sin # default is not in the list
        end
    end

    el = NumberField(0.0:π:20)
    @test default(el) == 0
    el = NumberField(0.0:π:20; default = π)
    @test default(el) == Float64(π) # should have been converted to Float64 because our range has been

    el = TextField()
    @test default(el) == ""
    el = TextField(; default = "xoxo")
    @test default(el) == "xoxo"
    el = TextField((5, 50); default = "xoxo")
    @test default(el) == "xoxo"

    el = PasswordField()
    @test default(el) == ""
    el = PasswordField(; default = "xoxo")
    @test default(el) == "xoxo"

    el = RangeSlider(1:10)
    @test default(el) == 1:10
    el = RangeSlider(1:10; show_value = true)
    @test default(el) == 1:10

    el = RangeSlider(0.0:π:10; default = 0:10)
    @test default(el) == 0.0:π:10

    el = RangeSlider(1:(1//3):10)
    @test default(el) == 1:(1//3):10
    el = RangeSlider(1:(1//3):10; default = 4:5)
    @test default(el) == 4:1//3:5
    el = RangeSlider(1:(1//3):10; default = 4:1//3:(17//3))
    @test default(el) == 4:1//3:(17//3)

    el = confirm(Slider([sin, cos]))
    @test default(el) == sin

    import PlutoUI: combine
    el = combine() do Child
        # @htl instead of md" because Julia VS Code is too buggy
        @htl("""
        # Hi there!

        I have $(Child(Slider([sin, cos]))) dogs and $(Child(Slider(5:10))) cats.

        Would you like to see them? $(Child(CheckBox(true)))
        """)
    end

    @test default(el) == (sin, 5, true)

    el = combine() do Child
        # @htl instead of md" because Julia VS Code is too buggy
        @htl("""
        # Hi there!

        I have $(Child(:fun, Slider([sin, cos]))) dogs and $(Child(:x, Slider(5:10))) cats.

        Would you like to see them? $(Child(:y, CheckBox(true)))
        """)
    end

    @test default(el) == (fun = sin, x = 5, y = true)

    el = combine() do Child
        @htl("""
        $(Child(html"<input>"))
        """)
    end

    @test default(el) === (missing,)

    import PlutoUI.Experimental: transformed_value

    el = let
        old_widget = PlutoUI.combine() do Child
            @htl("""$(Child(TextField(default="fons"))) $(Child(Slider(1:10; default=3)))""")
        end

        # Note that the input to `transform` is now a Tuple!
        # (This is the output of `PlutoUI.combine`)
        transform = input -> repeat(input[1], input[2])

        # use `transformed_value` to add the value tranformation to our widget
        transformed_value(transform, old_widget)
    end

    @test default(el) == "fonsfonsfons"

    el = transformed_value(
        transformed_value(
            PlutoUI.combine() do Child
                @htl("""
                $(Child(Slider([sin, cos])))
                $(Child(Slider(5:7)))
                """)
            end
        ) do x
            Ref(reverse(x))
        end
    ) do x
        first(x[])
    end

    @test_logs min_level=Logging.Info transformed_value(log, html"<input type=range>")

    @test default(el) == 5

    import PlutoUI.Experimental: wrapped

    el = wrapped() do Child
        @htl("""
        Hello!
        $(Child(Slider([sin, cos])))
        """)
    end

    @test default(el) == sin

end


text:
function mpo_triangular_AF_ising()
    # exact: 0.3230659669
    # ref: Phys. Rev. Res. 3, 013041 (2021)
    t = TensorMap(zeros, ComplexF64, ℂ^2*ℂ^2, ℂ^2)
    p = TensorMap(zeros, ComplexF64, ℂ^2, ℂ^2*ℂ^2)
    t[1, 1, 2] = 1
    t[1, 2, 1] = 1
    t[2, 1, 1] = 1
    t[2, 2, 1] = 1
    t[2, 1, 2] = 1
    t[1, 2, 2] = 1
    p[1, 1, 1] = 1
    p[2, 2, 2] = 1
    T = t*p
    return T
end

function mpo_triangular_AF_ising_alternative()
    # exact: 0.3230659669
    # ref: Phys. Rev. Res. 3, 013041 (2021)
    t = TensorMap(zeros, ComplexF64, ℂ^4*ℂ^4, ℂ^4)
    t[1, 2, 3] = 1
    t[3, 1, 2] = 1
    t[2, 3, 1] = 1
    t[3, 2, 4] = 1
    t[2, 4, 3] = 1
    t[4, 3, 2] = 1
    @tensor T[-1, -2; -3, -4] := t'[-2, 1, -4] * t[-1, 1, -3]
    return T
end

"""
    mpo_triangular_AF_ising_adapter()

    generate an MPO that converts the boundary MPS between the basis of `mpo_triangular_AF_ising_alternative()` and the basis of `mpo_triangular_AF_ising()`
"""
function mpo_triangular_AF_ising_adapter()

    # basis (s1, s2) to basis s
    p = TensorMap(zeros, ComplexF64, ℂ^4*ℂ^4, ℂ^4)
    m = TensorMap(zeros, ComplexF64, ℂ^2, ℂ^4*ℂ^4)

    p[1, 1, 1] = p[2, 2, 3] = p[3, 3, 2] = p[4, 4, 4] = 1
    m[1, 1, 1] = m[1, 1, 2] = m[1, 3, 1] = m[1, 3, 2] = 1
    m[2, 2, 3] = m[2, 2, 4] = m[2, 4, 3] = m[2, 4, 4] = 1
    @tensor T[-1, -2; -3, -4] := p[-1, 1, -3] * m[-2, 1, -4]

    # basis s to basis (s1, s2)
    p = TensorMap(zeros, ComplexF64, ℂ^2*ℂ^2, ℂ^2)
    m = TensorMap(zeros, ComplexF64, ℂ^4, ℂ^2*ℂ^2)

    p[1, 1, 1] = p[2, 2, 2] = 1
    m[1, 1, 1] = m[2, 2, 1] = m[3, 1, 2] = m[4, 2, 2] = 1
    @tensor Trev[-1, -2; -3, -4] := p[-1, 1, -3] * m[-2, 1, -4]

    return T, Trev

end

function mpo_kink_processor(decay_rate::Float64)
    gate = TensorMap(zeros, ℂ^2*ℂ^2, ℂ^2*ℂ^2)
    gate[1,1,1,1] = gate[2,2,2,2] = exp(-decay_rate)
    gate[1,2,1,2] = gate[2,1,2,1] = 1

    u, s, v = tsvd(permute(gate, (1,3), (2,4)))
    p = isometry(ℂ^4, ℂ^2)
    lft = permute(u * sqrt(s) * p, (1,), (2,3))
    rft = permute(p' * sqrt(s) * v, (1,2), (3,))
    T = rft * lft
    return T
end

function mpo_square_ising(beta::Float64)
    # ising model on square lattice
    δ = TensorMap(zeros, ComplexF64, ℂ^2*ℂ^2, ℂ^2*ℂ^2)
    t = TensorMap(zeros, ComplexF64, ℂ^2, ℂ^2)
    δ[1, 1, 1, 1] = 1
    δ[2, 2, 2, 2] = 1
    t[1, 1] = exp(-beta)
    t[2, 2] = exp(-beta)
    t[1, 2] = exp(beta)
    t[2, 1] = exp(beta)

    u, s, v = tsvd(t)
    u = u * sqrt(s)
    v = sqrt(s) * v

    @tensor T[-1, -2; -3, -4] := δ[1, 2, 3, 4] * v[-1, 1] * v[-2, 2] * u[3, -3] * u[4, -4]
    return T
end


text:
module MedicalTransforms

include("./crop.jl")
include("./utils.jl")

export
    crop_foreground,
    spatial_crop,
    generate_spatial_bounding_box
end


text:
using NetworkEpidemics
using LightGraphs
using Statistics
using Plots
using ColorSchemes
using ProgressMeter

CorrectedMetapopulation(χ::Real, mp::Metapopulation{SI}) = Metapopulation(mp.h, mp.D, SI(χ*mp.dynamics.β))

CorrectedMetapopulation(χ::Real, mp::Metapopulation{SIS}) = Metapopulation(mp.h, mp.D, SIS(χ*mp.dynamics.β, mp.dynamics.γ))

CorrectedMetapopulation(χ::Real, mp::Metapopulation{SIR}) = Metapopulation(mp.h, mp.D, SIR(χ*mp.dynamics.β, mp.dynamics.δ))

function L2_norm(t, u, v)
    sum = 0.0
    dt = 0.0
    n = length(t)
    f = (u-v).^2
    for i in 1:n-1
        dt = t[i+1] - t[i]
        sum += dt*(f[i]+f[i+1])/2
    end
    return sqrt(sum)
end

norm_uniform(t, u, v) = maximum(abs.(u-v))

norm_mean(t, u, v) = mean(abs.(u-v))

function correction_effficiency(mpx::Metaplex, mp::Metapopulation, x0_mpx, x0_mp; tmax=100.0, nmax = 1000, nsims=1000, nbins=100, norm=norm_uniform)
    #ts_av_mpx, u_av_mpx = average(mpx, x0_mpx, tmax=tmax, nmax=nmax, nsims=nsims, nbins=nbins, progressbar=false)
    #ts_av_mp, u_av_mp = average(mp, x0_mp, tmax=tmax, nmax=nmax, nsims=nsims, nbins=nbins, progressbar=false)
    ts = LinRange(0.0,tmax, nbins)
    ts_mf_mpx, u_mf_mpx = meanfield(mpx, x0_mpx, tmax=tmax, saveat=ts)
    ts_mf_mp, u_mf_mp = meanfield(mp, x0_mp, tmax=tmax, saveat=ts)
    #return Linfty_norm(ts_av_mpx, sum(u_av_mpx[2], dims=2), sum(u_av_mp[2], dims=2))/N, Linfty_norm(ts_mf_mpx, sum(u_mf_mpx[2], dims=2), sum(u_mf_mp[2], dims=2))/N
    return norm(ts_mf_mpx, sum(u_mf_mpx[2], dims=2)/N, sum(u_mf_mp[2], dims=2))/N
end

function test(h, N, k, ps, β, D, x0_mpx, x0_mp; tmax=100.0, nmax=1000, nsims=1000, nbins=100, nsamples=50)
    n = length(ps)
    #norm_av = zeros(n)
    norm_mf = zeros(n)
    si = SI(β)
    p = Progress(n*nsamples, dt=5.0)
    for i in 1:n
        for j in 1:nsamples
            g = watts_strogatz(N, k, ps[i])
            χ = mean(degree(g))
            mp = Metapopulation(h, D, SI(χ*β))
            mpx = Metaplex(g, h, D, SI(β))
            mf = correction_effficiency(mpx, mp, x0_mpx, x0_mp, tmax=tmax, nmax=nmax, nsims=nsims, nbins=nbins)
            #norm_av[i] += av
            norm_mf[i] += mf
            ProgressMeter.next!(p)
        end
    end
    #return norm_av/nsamples, norm_mf/nsamples
    return norm_mf/nsamples
end

N = 1000
M = 10
h = complete_graph(M)
β = 0.1
D = 0.1

x0_i = fill(1, N)
x0_i[1:10] .= 2
x0_μ = [div(i-1,div(N,M))+1 for i in 1:N]

x0_mp = Array{Int,2}(undef, M,2)
x0_mp[2:M,:] .= [100 0]
x0_mp[1, :] .= [90, 10]

x0_i = fill(1, N)
x0_i[1:10] .= 2
x0_μ = rand(1:M, N)
x0_mp = zeros(Int, M, 2)
for i in 1:N
    x0_mp[x0_μ[i], x0_i[i]] += 1
end
x0_mpx =  [x0_i, x0_μ]

ps = LinRange(0.0, 1.0, 50)
tmax = 50.0
nmax = 5000
nsims = 200
nbins = 500
nsamples = 5

@time norm_mf = test(h, N, 100, ps, β, D, x0_mpx, x0_mp, nmax=nmax, nsims=nsims, nbins=nbins, tmax=tmax, nsamples=nsamples)

scatter(ps, norm_mf,
        label="Mean Field",
        xlabel="p",
        ylabel="Efficiency of mean degree correction")
#scatter!(ps, norm_mf,
#        label="Mean Field")

using DelimitedFiles
open("efficiency.dat", "w") do io
    writedlm(io, [ps norm_av norm_mf])
end

savefig("efficiency.pdf")

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
#######################################
#
# Standard electrode sets
#
#######################################

EEG_64_10_20 = [
    "Fp1",
    "F3",
    "FC5",
    "C3",
    "CP5",
    "P3",
    "PO7",
    "Oz",
    "Fpz",
    "AFz",
    "F6",
    "FC4",
    "C2",
    "TP8",
    "P2",
    "P10",
    "AF7",
    "F5",
    "FC3",
    "C5",
    "CP3",
    "P5",
    "PO3",
    "POz",
    "Fp2",
    "Fz",
    "F8",
    "FC2",
    "C4",
    "CP6",
    "P4",
    "PO8",
    "AF3",
    "F7",
    "FC1",
    "T7",
    "CP1",
    "P7",
    "O1",
    "Pz",
    "AF8",
    "F2",
    "FT8",
    "FCz",
    "C6",
    "CP4",
    "P6",
    "PO4",
    "F1",
    "FT7",
    "C1",
    "TP7",
    "P1",
    "P9",
    "Iz",
    "CPz",
    "AF4",
    "F4",
    "FC6",
    "Cz",
    "T8",
    "CP2",
    "P8",
    "O2",
]

EEG_Vanvooren_2014 = [
    "TP7",
    "P9",
    "P7",
    "P5",
    "P3",
    "P1",
    "PO7",
    "PO3",
    "O1",
    "P2",
    "P4",
    "P6",
    "P10",
    "TP8",
    "PO4",
    "PO8",
    "O2",
    "TP8",
]

EEG_Vanvooren_2014_Left = ["TP7", "P9", "P7", "P5", "P3", "P1", "PO7", "PO3", "O1"]

EEG_Vanvooren_2014_Right = ["P2", "P4", "P6", "P10", "TP8", "PO4", "PO8", "O2", "TP8"]


text:
########## reward ##########
#POMDPs.reward(DP::DrivePOMDP, Ss::Sts, acc_ego::Float64) = Reward(DP, Ss, acc_ego)
POMDPs.reward(DP::DrivePOMDP, Ss::Sts, Aego::Symbol, Ssnext::Sts) = Reward(DP, Ss, Aego, Ssnext)

# based directly on next state; used in DiscreteValueIteration Solver
function R_crash(DP::DrivePOMDP, Ssnext::Sts, acc_ego::Float64)
    Ego_nextPos = GetGlobalPosition(DP.Routes[Ssnext.Ego.r], min(Ssnext.Ego.s, DP.sV[end]))
    Other_nextPos = GetGlobalPosition(DP.Routes[Ssnext.Other.r], min(Ssnext.Other.s, DP.sV[end]))
    if haskey(DP.Routes[Ssnext.Ego.r].intersect_Infos, Ssnext.Other.r) # if two routes intesect with each other
        if Distance2D(Ego_nextPos, Other_nextPos) < 2 # crash happens
            return -100.0 - 1.5*Ssnext.Ego.v^2 - 2.6*acc_ego*abs(acc_ego)
        end
    end
    return 0.0
end

# include the situation that junction point is far away from every state-point but collision actually happens
function R_crash(DP::DrivePOMDP, Ss::Sts, Ssnext::Sts)
    if haskey(DP.Routes[Ss.Ego.r].intersect_Infos, Ss.Other.r) # if two routes intesect with each other
        for (i, PtDist) in enumerate(DP.Routes[Ss.Ego.r].intersect_Infos[Ss.Other.r])
            InterDistEgo = PtDist[3] # distance from start point to junction point
            InterDistOther = DP.Routes[Ss.Other.r].intersect_Infos[Ss.Ego.r][i][3]
            if Ss.Ego.s < InterDistEgo < Ssnext.Ego.s && Ss.Other.s < InterDistOther < Ssnext.Other.s
                return -150.0
            end
        end
    end
    return 0.0
end

"""
if Vego != Vref, R_v should always be less than R_acc, so that Ego will get a better reward if it chooses to accelerate/decelerate itself.
R_v and R_acc should have a relation between each other.
"""

function R_v(DP::DrivePOMDP, Ss::Sts, acc_ego::Float64)
    Kvp = 1.5
    Kvn = 1.0
    Vref = DP.Routes[Ss.Ego.r].Vref[UInt16(Ss.Ego.s/DP.Δs)+1] # m/s, should be defined as array according to Route/Geometries
    Vego = Ss.Ego.v + acc_ego*DP.Δt
    if Vego > Vref
        #return -Kvp*(Vref - Vego)^2
        return -Kvp*((Vref - Vego)*3.6)^2
    else
        return -Kvn*(Vref - Vego)*3.6
    end
end

function R_acc(DP::DrivePOMDP, Ss::Sts, acc_ego::Float64) #
    Factor = -0.77
    Δa = acc_ego - DP.Routes[Ss.Ego.r].Aref[UInt16(Ss.Ego.s/DP.Δs)+1]
    return Factor*(Δa^2+2.3*abs(Δa))
    #return Factor*(acc_ego^2+2.3*abs(acc_ego))
end

function Reward(DP::DrivePOMDP, Ss::Sts, Aego::Symbol, Ssnext::Sts)
    acc_ego = 0.0
    if Ssnext.Other.s != Ss.Other.s
        acc_k = (Ssnext.Other.v^2 - Ss.Other.v^2)/(2*(Ssnext.Other.s - Ss.Other.s))
        Ego, Other, accs = IDMtransit(DP, Ss, Aego, acc_k)
        if Ego.s != Ss.Ego.s
            acc_ego = (Ego.v^2 - Ss.Ego.v^2)/(2*(Ego.s - Ss.Ego.s))
        end
        #acc_ego = (Ssnext.Ego.v^2 - Ss.Ego.v^2)/(2*(Ssnext.Ego.s - Ss.Ego.s))
    end
    reward = R_crash(DP, Ss, Ssnext) + R_crash(DP, Ss, acc_ego) + R_crash(DP, Ssnext, acc_ego) + R_v(DP, Ss, acc_ego) + R_acc(DP, Ss, acc_ego)
    reward = trunc(reward; digits=4)
    #println("## Realtime Reward:$reward")
    return reward
end


text:
# ------------------------------------------------------------------------------------------
# # Introduction
#
# **DynamicalSystems.jl** is a Julia software library for the exploration of chaos and
# nonlinear dynamics.
#
# ## Contents of these tutorials
#
# 1. How to create a dynamical system and get a trajectory out of it.
# 3. Orbit diagrams & Poincaré surfaces of section of a dynamical system.
# 4. Lyapunov exponents of a dynamical system.
# 5. Entropies and dimensions of numerical data.
# 6. Delay coordinates embedding & Neighborhoods.
# 7. Quick summary of the rest of the features possible with **DynamicalSystems.jl**:
#     * Lyapunov exponent from numerical data
#     * Broomhead-King coordinates
#     * Finding fixed points of any map of any order
#     * Detecting chaos with GALI
#
# 8. How to predict timeseries (*still in beta!*).
#
# ## About DynamicalSystems.jl
#
# 1. **DynamicalSystems.jl** has been built from the ground up with the principles of
# clarity, intuition and robustness.
#
# 3. Aims to be a "library" in the literal sense: all offered functionality is well
# documented and accompanied with descriptions of the algorithms as well as references to
# original papers. This offers new users plenty of chances to learn new algorithms.
#
# ## Disclaimer
#
# **DynamicalSystems.jl** is in v1.0-alpha! Comments, testing, suggestions etc. are very
# welcome before official release.
#
# Any changes that will happen for v1.0 release will be summarized in the [News
# page](https://juliadynamics.github.io/DynamicalSystems.jl/latest/news/).
#
# ---
#
# # Documentation
#
# https://juliadynamics.github.io/DynamicalSystems.jl/latest/.
#
# In the same documentation page you will find a [list of
# contents](https://juliadynamics.github.io/DynamicalSystems.jl/latest/#contents).
#
# If you have questions, you can open issues in the appropriate [GitHub
# repository](https://github.com/JuliaDynamics); if you think your question is minor, don't
# hesitate to chat with us over [gitter](https://gitter.im/JuliaDynamics/Lobby).
#
# **DynamicalSystems.jl** is a registered package. On Julia v0.6 you only have to type
# `Pkg.add("DynamicalSystems")` to install its features.
#
# ## Documentation Strings
#
# All functions that are exported by **DynamicalSystems.jl** have very detailed
# documentation strings.
# * Simply press `?` followed by the function name to see the documentation strings!
#
# ---
#
# # About this tutorial series
# Tutorials 1-4 are intended for any audience *but* with some prior exposure to math and
# physics, as outlined in the prerequisites. In these tutorials, the methods shown will be
# explained and discussed in detail, so people unfamiliar with nonlinear dynamics and chaos
# can still follow!
#
# Tutorials 5 and onwards will target more a experienced audience. These tutorials will
# provide example applications, but will not explain the methods/algorithms used. This is so
# that we can demo the available features in a relatively short period of time. Therefore,
# don't feel worried if you are not able to fully grasp all concepts at that point!
#
# ## Prerequisites
# This tutorial series assumes some basic knowledge of Julia & physics/mathematics concepts.
#
# * Julia
#     * [Functions](https://docs.julialang.org/en/stable/manual/functions/#man-functions-1).
#     * [Broadcasting](https://docs.julialang.org/en/stable/manual/arrays/#Broadcasting-1)
# ("vectorizing" function application).
#     * [Methods](https://docs.julialang.org/en/stable/manual/methods/#Methods-1) and
# multiple dispatch. This also relies on basic knowledge of the
# [Type](https://docs.julialang.org/en/stable/manual/types/#man-types-1) system.
#     * Basic handling of vectors and matrices.
#
#
# * Mathematics/Physics
#     * State-space & trajectory of a dynamical system
#     * What does it mean for a trajectory to be periodic
#     * Jacobians & linearized dynamics
#
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/introduction-to-dynamicalsystems.jl/0. Introduction, Documentation and Prerequisites.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# RwaJur2

export VERTICES_RWAJUR2,
		RwaJur2Exposure,
		RwaJur2Parameters,
		getRwaJur2Parameters,
		getLongExposureArray,
		getShortExposureArray,
		getRwaJur2NetExposure,
		getRwaJur2VerticalGap,
		getRwaJur2Zones,
		getRwaJur2HorizontalGapWithinZones,
		getRwaJur2HorizontalGapBetweenZones,
		getPJur2,
		getRwaJur2

"""
Standardized vertices for mapping foreign currency coupon rate exposures (RWAJUR2)
"""
const VERTICES_RWAJUR2 = [1, 21, 42, 63, 126, 252, 504, 756, 1008, 1260, 2520]

"""
    RwaJur2Exposure

Holds exposures in foreign currency coupon rates.

# Fields
* `date::Date`: date of evaluation
* `VXX::Float64`: exposure to vertex XX,
given by the sum of present value in BRL
of all cashflows that matures in XX business days (BD).

# Constructors
    RwaJur2Exposure(date::Date)

Create an instance with all exposures equal to zero.

	RwaJur2Exposure(date::Date, currency::AbstractString, longExposures::Array{Float64, 1}, shortExposures::Array{Float64, 1})

Create an instance with exposures given by two vectors of length 10 containing long and short exposures.
Exposures are assigned in sequence from vertex 1 BD to vertex 2520 BD.

"""
type RwaJur2Exposure
	# Date of evaluation
	date::Date
	# Currency
	currency::AbstractString
	# Exposures in standard vertices: VXX, where XX is a maturity calculated in business days
	# LONG EXPOSURES
	V1Long::Float64
	V21Long::Float64
	V42Long::Float64
	V63Long::Float64
	V126Long::Float64
	V252Long::Float64
	V504Long::Float64
	V756Long::Float64
	V1008Long::Float64
	V1260Long::Float64
	V2520Long::Float64
	# SHORT EXPOSURES
	V1Short::Float64
	V21Short::Float64
	V42Short::Float64
	V63Short::Float64
	V126Short::Float64
	V252Short::Float64
	V504Short::Float64
	V756Short::Float64
	V1008Short::Float64
	V1260Short::Float64
	V2520Short::Float64

	# Constructors
	RwaJur2Exposure(date::Date, currency::AbstractString) = new(date, currency, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

	function RwaJur2Exposure(date::Date, currency::AbstractString, longExposures::Array{Float64, 1}, shortExposures::Array{Float64, 1})
		if size(longExposures, 1) != 11
			error("long exposures must have 11 values")
		end #if
		if any(longExposures .< 0)
			error("all long exposures must be positive")
		end #if

		if size(shortExposures, 1) != 11
			error("short exposures must have 11 values")
		end #if
		if any(shortExposures .< 0)
			error("all short exposures must be positive")
		end #if

		new(date,
			currency,
		    longExposures[1], longExposures[2], longExposures[3], longExposures[4],
			longExposures[5], longExposures[6], longExposures[7], longExposures[8],
			longExposures[9], longExposures[10], longExposures[11],
			shortExposures[1], shortExposures[2], shortExposures[3], shortExposures[4],
			shortExposures[5], shortExposures[6], shortExposures[7], shortExposures[8],
			shortExposures[9], shortExposures[10], shortExposures[11]
			)
	end #function

end #type

"""
	RWAJUR1Parameters

Holds BCB parameters to calculate foreign currency coupon stantardized capital (RWA JUR2)

"""
type RwaJur2Parameters
	date::Date
	M2::Float64
end #type

"""
    getRwaJur2Parameters(date::Date)

Return RWA JUR2 parameters supplied by BCB for `date::Date`.
"""
function getRwaJur2Parameters(date::Date)
	if date <= Date(2012, 4, 29)
		M2 = 2.28
	elseif date <= Date(2012, 8, 30)
		M2 = 2.75
	elseif date <= Date(2012, 12, 30)
		M2 = 3.22
	else
		M2 = 3.7
	end #if
	return(RwaJur2Parameters(date, M2))
end #function

"""
    getLongExposureArray(rExp::RwaJur2Exposure)

Return a vector with all long exposures in sequence,
from vertex 1 BD to vertex 2520 BD.

"""
getLongExposureArray(rExp::RwaJur2Exposure) = [rExp.V1Long, rExp.V21Long, rExp.V42Long, rExp.V63Long, rExp.V126Long, rExp.V252Long, rExp.V504Long, rExp.V756Long, rExp.V1008Long, rExp.V1260Long, rExp.V2520Long]

"""
    getShortExposureArray(rExp::RwaJur2Exposure)

Return a vector with all short exposures in sequence,
from vertex 1 BD to vertex 2520 BD.

"""
getShortExposureArray(rExp::RwaJur2Exposure) = [rExp.V1Short, rExp.V21Short, rExp.V42Short, rExp.V63Short, rExp.V126Short, rExp.V252Short, rExp.V504Short, rExp.V756Short, rExp.V1008Short, rExp.V1260Short, rExp.V2520Short]

"""
    getRwaJur2NetExposure(rExp::RwaJur2Exposure)

Return a vector containing RWA JUR2 weighted net exposures.

"""
function getRwaJur2NetExposure(rExp::RwaJur2Exposure)
	Y = ([0.0, 0.5, 0.7, 0.8, 1.2, 2.0, 4.0, 6.0, 8.0, 10.0, 18.0]/100)
	(getLongExposureArray(rExp) - getShortExposureArray(rExp)) .* Y
end #function

"""
    getRwaJur2VerticalGap(rExp::RwaJur2Exposure)

Return a vector containing RWA JUR2 vertical gaps.

"""
function getRwaJur2VerticalGap(rExp::RwaJur2Exposure)
	Y = ([0.0, 0.5, 0.7, 0.8, 1.2, 2.0, 4.0, 6.0, 8.0, 10.0, 18.0]/100)
	0.1 * ( min(abs(getLongExposureArray(rExp)), abs(getShortExposureArray(rExp))) .* Y)
end #function

"""
    getRwaJur2Zones(rExp::RwaJur2Exposure)

Return a matrix containing RWA JUR2 zone exposures.
Column 1 contains long exposures, column2 contains short exposures.
Zones are represented in each of the 3 rows.

"""
function getRwaJur2Zones(rExp::RwaJur2Exposure)
	EL = getRwaJur2NetExposure(rExp)

	Z = zeros(3, 2)

	# Zone #1
	for i = 1:5
		if EL[i] >= 0
			Z[1, 1] += EL[i]
		else
			Z[1, 2] += EL[i]
		end #if
	end #for

	# Zone #2
	for i = 6:8
		if EL[i] >= 0
			Z[2, 1] += EL[i]
		else
			Z[2, 2] += EL[i]
		end #if
	end #for

	# Zone #3
	for i = 9:11
		if EL[i] >= 0
			Z[3, 1] += EL[i]
		else
			Z[3, 2] += EL[i]
		end #if
	end #for

	return(Z)
end #function

"""
    getRwaJur2HorizontalGapWithinZones(rExp::RwaJur2Exposure)

Return a vector containing RWA JUR2 horizontal gap within zones.

"""
function getRwaJur2HorizontalGapWithinZones(rExp::RwaJur2Exposure)

	Z = getRwaJur2Zones(rExp)

	DHZ = zeros(3)
	W = [0.4, 0.3, 0.3]
	for i = 1:3
		DHZ[i] = min(Z[i, 1], abs(Z[i, 2])) * W[i]
	end #for

	return(DHZ)

end #function

"""
    getRwaJur2HorizontalGapBetweenZones(rExp::RwaJur2Exposure)

Return a vector containing RWA JUR2 horizontal gap between zones.

"""
function getRwaJur2HorizontalGapBetweenZones(rExp::RwaJur2Exposure)
	Z = getRwaJur2Zones(rExp)
	DHE = sum(Z, 2)[:,1]
	DHEZ = 0

	if DHE[1] * DHE[2] < 0
		DHEZ += min(abs(DHE[1]), abs(DHE[2])) * 0.4
	end #if
	if DHE[2] * DHE[3] < 0
		DHEZ += min(abs(DHE[2]), abs(DHE[3])) * 0.4
	end #if
	if DHE[1] * DHE[3] < 0
		DHEZ += min(abs(DHE[1]), abs(DHE[3])) * 1.0
	end #if

	return(DHEZ)

end #function

"""
	getPJur2(rExp::RwaJur2Exposure, par::RwaJur2Parameters=getRwaJur2Parameters(rExp.date))

Calculate PJUR2 using maturity ladder methodology.

"""
function getPJur2(rExp::RwaJur2Exposure, par=getRwaJur2Parameters(rExp.date))
	EL = getRwaJur2NetExposure(rExp)
	DV = getRwaJur2VerticalGap(rExp)
	DHZ = getRwaJur2HorizontalGapWithinZones(rExp)
	DHE = getRwaJur2HorizontalGapBetweenZones(rExp)

	par.M2 * (abs(sum(EL)) + sum(abs(DV)) + sum(abs(DHZ)) + DHE)
end #function

"""
	getRwaJur2(rExp::RwaJur2Exposure, par::RwaJur2Parameters=getRwaJur2Parameters(rExp.date))

Calculate risk-weighted assets for foreign currency coupon exposures (RWA JUR2),
using maturity ladder methodology.

"""
function getRwaJur2(rExp::RwaJur2Exposure, par=getRwaJur2Parameters(rExp.date))
	(1/(getF(par.date))) * getPJur2(rExp, par)
end #function


text:
@testset "h-for" begin
    F.def_LOCAL_VARS!()
    s = """
        @def list = [1,2,3]
        """ |> fd2html_td
    hs = raw"""
        ABC
        {{for x in list}}
          {{fill x}}
        {{end}}
        """
    tokens = F.find_tokens(hs, F.HTML_TOKENS, F.HTML_1C_TOKENS)
    hblocks, tokens = F.find_all_ocblocks(tokens, F.HTML_OCB)
    qblocks = F.qualify_html_hblocks(hblocks)
    @test qblocks[1] isa F.HFor
    @test qblocks[1].vname == "x"
    @test qblocks[1].iname == "list"
    @test qblocks[2] isa F.HFun
    @test qblocks[3] isa F.HEnd

    content, head, i = F.process_html_for(hs, qblocks, 1)
    @test isapproxstr(content, "1 2 3")
end

@testset "h-for2" begin
    F.def_LOCAL_VARS!()
    s = """
        @def list = ["path/to/badge1.png", "path/to/badge2.png"]
        """ |> fd2html_td
    h = raw"""
        ABC
        {{for x in list}}
            {{fill x}}
        {{end}}
        """ |> F.convert_html
    @test isapproxstr(h, """
        ABC
        path/to/badge1.png
        path/to/badge2.png
        """)
end

@testset "h-for3" begin
    F.def_LOCAL_VARS!()
    s = """
        @def iter = (("a", 1), ("b", 2), ("c", 3))
        """ |> fd2html_td
    h = raw"""
        ABC
        {{for (n, v) in iter}}
            name:{{fill n}}
            value:{{fill v}}
        {{end}}
        """ |> F.convert_html
    @test isapproxstr(h, """
        ABC
        name:a
        value:1
        name:b
        value:2
        name:c
        value:3
        """)

    s = """
        @def iter2 = ("a"=>10, "b"=>7, "c"=>3)
        """ |> fd2html_td
    h = raw"""
        ABC
        {{for (n, v) in iter2}}
            name:{{fill n}}
            value:{{fill v}}
        {{end}}
        """ |> F.convert_html
    @test isapproxstr(h, """
        ABC
        name:a
        value:10
        name:b
        value:7
        name:c
        value:3
        """)
end


text:
function set_volume(nodes, cellxmax, cellymax)
    volume = zeros(cellxmax, cellymax)
    for i in 1:cellxmax
        for j in 1:cellymax
            vec_r1x = nodes[i+1,j+1,1] - nodes[i,j,1]
            vec_r1y = nodes[i+1,j+1,2] - nodes[i,j,2]
            vec_r2x = nodes[i,j+1,1] - nodes[i+1,j,1]
            vec_r2y = nodes[i,j+1,2] - nodes[i+1,j,2]

            volume[i,j] = abs(vec_r1x*vec_r2y - vec_r1y*vec_r2x) /2
        end
    end
    return volume
end

function set_mu(Qbase,cellxmax,cellymax,specific_heat_ratio,Rd)
    mu = zeros(cellxmax,cellymax)

    # サザーランドの式
    # https://www.jstage.jst.go.jp/article/jsam1937/37/4/37_4_694/_pdf/-char/ja
    mu0 = 1.82e-5 # 基準粘度[Pa s]
    T0 = 293.15       # 基準温度[K]
    C = 117       # サザーランド定数[K]

    for i in 1:cellxmax
        for j in 1:cellymax
            # T = p/(rho Rd )
            T = Qbase[i,j,4]/(Qbase[i,j,1]*Rd)

            if T <0
                println("T is munus ( 'A' ) !!")
                println("i : "*string(i)*"  j : "*string(j))
                println(Qbase[i,j,:])
            end

            mu[i,j] = mu0 * (T/T0)^1.5 * (T0+C)/(T+C)
        end
    end

    return mu
end

function set_lambda(Qbase,cellxmax,cellymax,mu,specific_heat_ratio,Rd)

    lambda = zeros(cellxmax,cellymax)

    # サザーランドの式
    lam0 = 22.3*10^(-3)  # 基準熱伝導率　[W/mK]
    T0 = 273.15       # 基準温度[K]
    C = 125       # サザーランド定数[K]

    for i in 1:cellxmax
        for j in 1:cellymax
            T = Qbase[i,j,4]/(Qbase[i,j,1]*Rd)
            lambda[i,j] = lam0*((T0+C)/(T+C))*(T/T0)^1.5
        end
    end

    return lambda
end

function set_Minf(bdcon, specific_heat_ratio)
    rho = 0
    u   = 0
    v   = 0
    p   = 0

    for i in 1:4
        if Int(bdcon[i][1]) == 0 || Int(bdcon[i][1]) == 5
            rho = bdcon[i][2]
            u   = bdcon[i][3]
            v   = bdcon[i][4]
            p   = bdcon[i][5]
        end
    end

    a = (specific_heat_ratio * p / rho)^0.5
    u = (u^2 + v^2)^0.5
    M = u/a

    return M
end

function set_gasconst(Qbase,cellxmax,cellymax,nval,nch,R)
    Rhat = zeros(cellxmax, cellymax)
    mw   = [28e-3, 14e-3]
    npre = nval - nch

    for i in 1:cellxmax
        for j in 1:cellymax
            for s in 1:nch
                Rs = R/mw[s]
                Rhat[i,j] += Qbase[i,j,npre+s] * Rs
            end
        end
    end
    return Rhat
end


text:
#This script calculates the harmonic correction for a 3-level system.
#Here npulse is fixed.
using DifferentialEquations
using LinearAlgebra
using Plots; pyplot()
using NPZ
using LaTeXStrings

include("MagnusExpansion.jl")
include("HarmonicCorrection3ls.jl")
include("3ls.jl")
include("Envelopes.jl")

function save_pulse(α, tf, coeffs, pulse, basis)
	t_array = collect(0: tf/100: tf)
	output = zeros(length(t_array), 6)
	for i = 1: length(t_array)
		t = t_array[i]
		output[i, 1] = abs(α)*t
		output[i, 2] = Δ(t, tf, coeffs[1], basis.detuning)./abs(α)
		f = pulse(t)/abs(α)
		output[i, 3] = f
		g = g_envelopes(t, tf, coeffs[2], basis.pulse)./abs(α)
		output[i, 5] = (f + g[1])
		output[i, 6] = g[2]
	end
	plot(output[:, 1], output[:, [3, 2, 5, 6]],
		 linewidth = 2,
		 labels = [L"f_x^{(0)}(t)" L"\Sigma \Delta^{(6)}" L"f_x^{(6)}(t)" L"f_y^{(6)}(t)"],
		 linestyle = [:solid :dash :dash :dash])
	plot!(tickfont = font(12, "Serif"),
		  xlabel = "time "*L"|\alpha| t",
		  ylabel = "Envelopes "*L"\times |\alpha|^{-1}",
		  guidefont = font(14, "Serif"),
		  legendfont = font(12, "Serif"),
		  background_color_legend = false,
		  margin = 5Plots.mm)
	savefig("pulse_envelopes.pdf")
	npzwrite("pulse"*string(tf)*".npy", output)
	return nothing
end

function run_(magnus_order)
	λ = gellmann_matrices()
	α = -0.1
	η = sqrt(2.0)
	θ = pi/2
	# setting the number of harmonics of the correction basis
	nΔ = 0
	npulse = 2
	# array with different gate times tf
	tf0, tf1 = 50.0, 155.0
	nt = 100
	tf_array = collect(tf0: (tf1 - tf0)/(nt - 1): tf1)
	nt = size(tf_array)[1]
	F0, Fmagnus, Fdrag = zeros(nt), zeros(nt, length(magnus_order)), zeros(nt)

	for j = 1: nt
		tf = tf_array[j]
		println(tf)

		pulse = @closure t-> gaussian_pulse(t, tf, θ)
		pulse_integral = @closure t-> gaussian_pulse_integral(t, tf, θ)
		pulse_derivative = @closure t-> gaussian_pulse(t, tf, θ; derivative=true)
		system_params = Dict("alpha"=> α, "eta"=> η, "tf"=> tf, "pulse"=> pulse,
							 "pulse_integral"=> pulse_integral,
							 "pulse_derivative"=> pulse_derivative)
		# Ideal evolution
		U_ideal = U(system_params, λ, ideal_dynamics=true)
		# Nonideal evolution
		U_nideal = U(system_params, λ, ideal_dynamics=false)
		# Harmonic correction
		F0[j] = fidelity(U_nideal, U_ideal, 2)

		basis = fourier_basis(tf, nΔ, npulse)
		Δ_basis = basis.detuning
		pulse_basis = basis.pulse
		# correction matrix
		M = correction_matrix(system_params, Δ_basis, pulse_basis)
		# Magnus correction
		for k = 1: length(magnus_order)
			coeffs = get_coeffs(M, magnus_order[k], system_params, Δ_basis, pulse_basis)
			U_nideal2 = U(system_params, λ, coeffs, Δ_basis, pulse_basis)

			Fmagnus[j, k] = fidelity(U_nideal2, U_ideal, 2)
			if j == 1 && magnus_order[k] == 6
				save_pulse(α, tf, coeffs, pulse, basis)
			end
		end
		# Drag correction
		U_nideal3 = Udrag(system_params, λ; high_order=false)
		Fdrag[j] = fidelity(U_nideal3, U_ideal, 2)
	end
	return abs(α)*tf_array, F0, Fmagnus, Fdrag
end

magnus_order = [2, 6]
t, F0, Fmagnus, Fdrag = run_(magnus_order)
y = 1 .- hcat(F0, Fmagnus, Fdrag)
plot(t, y,
	 linewidth = 2,
	 labels = ["No\ncorrection" "2nd order\nMagnus" "6th order\nMagnus" "DRAG"],
	 linestyle = [:solid :solid :solid :dash])

plot!(yaxis=:log,
	  tickfont = font(12, "Serif"),
	  xlabel = "gate time "*L"|\alpha| t_\mathrm{f}",
	  ylabel = "Fidelity error "*L"\varepsilon",
	  guidefont = font(14, "Serif"),
	  legendfont = font(12, "Serif"),
	  background_color_legend = false,
	  margin = 5Plots.mm)
savefig("fidelity_error.pdf")

data = hcat(t, F0, Fmagnus, Fdrag)
npzwrite("data.npy", data)
# data[:, 1] .= t
# data[:, 2] .= F0
# data[:, 3] .= Fmagnus
# for m in magnus_order
# if m == 1
# 	npzwrite("fidelity_1st.npy", data)
# elseif m == 2
# 	npzwrite("fidelity_2nd.npy", data)
# else
# 	npzwrite("fidelity_"*string(m)*"th.npy", data)
# end

# data = zeros(nt,2)
# data[:,1] = abs(α)*tf_array
# data[:,2] = Fdrag
# npzwrite("fidelity_drag.npy", data)

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection