text:
mutable struct ABDF2ConstantCache{F,N,dtType,rate_prototype} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolve::N
  eulercache::ImplicitEulerConstantCache
  dtₙ₋₁::dtType
  fsalfirstprev::rate_prototype
end

function alg_cache(alg::ABDF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,
                   uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,1//1,1,ηold,z₊,dz,tmp,b,k))
  eulercache = ImplicitEulerConstantCache(uf,nlsolve)
  dtₙ₋₁ = one(dt)
  fsalfirstprev = rate_prototype
  nlsolve.cache.γ = 2//3
  nlsolve.cache.c = 1
  ABDF2ConstantCache(uf, nlsolve, eulercache, dtₙ₋₁, fsalfirstprev)
end

mutable struct ABDF2Cache{uType,rateType,uNoUnitsType,J,W,UF,JC,N,F,dtType} <: OrdinaryDiffEqMutableCache
  uₙ::uType
  uₙ₋₁::uType
  uₙ₋₂::uType
  du1::rateType
  fsalfirst::rateType
  fsalfirstprev::rateType
  k::rateType
  z::uType
  zₙ₋₁::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::J
  W::W
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
  eulercache::ImplicitEulerCache
  dtₙ₋₁::dtType
end

u_cache(c::ABDF2Cache)    = (c.z,c.dz)
du_cache(c::ABDF2Cache)   = (c.k,c.fsalfirst)

function alg_cache(alg::ABDF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,1//1,1,ηold,z₊,dz,tmp,b,k))
  atmp = similar(u,uEltypeNoUnits,axes(u))

  fsalfirstprev = similar(rate_prototype)
  eulercache = ImplicitEulerCache(u,uprev,uprev2,du1,fsalfirst,k,z,dz,b,tmp,atmp,J,W,uf,jac_config,linsolve,nlsolve)
  dtₙ₋₁ = one(dt)
  zₙ₋₁ = similar(u)
  nlsolve.cache.γ = 2//3
  nlsolve.cache.c = 1
  ABDF2Cache(u,uprev,uprev2,du1,fsalfirst,fsalfirstprev,k,z,zₙ₋₁,dz,b,tmp,atmp,J,
              W,uf,jac_config,linsolve,nlsolve,eulercache,dtₙ₋₁)
end

# SBDF

mutable struct SBDFConstantCache{rateType,F,N,uType} <: OrdinaryDiffEqConstantCache
  cnt::Int
  k2::rateType
  uf::F
  nlsolve::N
  uprev2::uType
  uprev4::uType
  uprev3::uType
  k₁::rateType
  k₂::rateType
  k₃::rateType
  du₁::rateType
  du₂::rateType
end

mutable struct SBDFCache{uType,rateType,uNoUnitsType,J,W,UF,JC,N,F} <: OrdinaryDiffEqMutableCache
  cnt::Int
  u::uType
  uprev::uType
  fsalfirst::rateType
  k::rateType
  du1::rateType
  z::uType
  dz::uType
  b::uType
  tmp::uType
  atmp::uNoUnitsType
  J::J
  W::W
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
  uprev2::uType
  uprev3::uType
  uprev4::uType
  k₁::rateType
  k₂::rateType
  k₃::rateType
  du₁::rateType
  du₂::rateType
end

u_cache(c::SBDFCache)    = ()
du_cache(c::SBDFCache)   = ()

function alg_cache(alg::SBDF,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  k2 = rate_prototype
  uf != nothing && ( uf = DiffEqDiffTools.UDerivativeWrapper(f.f1,t,p) )
  uprev2 = u
  uprev3 = u
  uprev4 = u
  k₁ = rate_prototype
  k₂ = rate_prototype
  k₃ = rate_prototype
  du₁ = rate_prototype
  du₂ = rate_prototype

  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,1//1,1,ηold,z₊,dz,tmp,b,k))
  SBDFConstantCache(1,k2,uf,nlsolve,uprev2,uprev3,uprev4,k₁,k₂,k₃,du₁,du₂)
end

function alg_cache(alg::SBDF,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  atmp = similar(u,uEltypeNoUnits,axes(u))
  order = alg.order
  uprev2 = similar(u)
  uprev3 = order >= 3 ? similar(u) : uprev2
  uprev4 = order == 4 ? similar(u) : uprev2
  k₁ = similar(rate_prototype)
  k₂ = order >= 3 ? zero(rate_prototype) : k₁
  k₃ = order == 4 ? zero(rate_prototype) : k₁
  du₁ = zero(rate_prototype)
  du₂ = zero(rate_prototype)
  if uf != nothing
    uf = DiffEqDiffTools.UJacobianWrapper(f.f1,t,p)
    linsolve = alg.linsolve(Val{:init},uf,u)
    jac_config = build_jac_config(alg,f,uf,du1,uprev,u,tmp,dz)
  end

  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,1//1,1,ηold,z₊,dz,tmp,b,k))
  SBDFCache(1,u,uprev,fsalfirst,k,du1,z,dz,b,tmp,atmp,J,W,uf,jac_config,linsolve,nlsolve,uprev2,uprev3,uprev4,k₁,k₂,k₃,du₁,du₂)
end

# QNDF1

mutable struct QNDF1ConstantCache{F,N,coefType,coefType1,dtType,uType} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolve::N
  D::coefType
  D2::coefType1
  R::coefType
  U::coefType
  uprev2::uType
  dtₙ₋₁::dtType
end

mutable struct QNDF1Cache{uType,rateType,coefType,coefType1,coefType2,uNoUnitsType,J,W,UF,JC,N,F,dtType} <: OrdinaryDiffEqMutableCache
  uprev2::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  b::uType
  D::coefType1
  D2::coefType2
  R::coefType
  U::coefType
  tmp::uType
  atmp::uNoUnitsType
  utilde::uType
  J::J
  W::W
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
  dtₙ₋₁::dtType
end

u_cache(c::QNDF1Cache)    = ()
du_cache(c::QNDF1Cache)   = ()

function alg_cache(alg::QNDF1,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  uprev2 = u
  dtₙ₋₁ = t

  D = fill(zero(typeof(u)), 1, 1)
  D2 = fill(zero(typeof(u)), 1, 2)
  R = fill(zero(typeof(t)), 1, 1)
  U = fill(zero(typeof(t)), 1, 1)

  U!(1,U)
  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,zero(inv((1-alg.kappa))),1,ηold,z₊,dz,tmp,b,k))

  QNDF1ConstantCache(uf,nlsolve,D,D2,R,U,uprev2,dtₙ₋₁)
end

function alg_cache(alg::QNDF1,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  D = Array{typeof(u)}(undef, 1, 1)
  D2 = Array{typeof(u)}(undef, 1, 2)
  R = fill(zero(typeof(t)), 1, 1)
  U = fill(zero(typeof(t)), 1, 1)

  D[1] = similar(u)
  D2[1] = similar(u); D2[2] = similar(u)

  U!(1,U)

  atmp = similar(u,uEltypeNoUnits,axes(u))
  utilde = similar(u,axes(u))
  uprev2 = similar(u)
  dtₙ₋₁ = one(dt)
  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,zero(inv((1-alg.kappa))),1,ηold,z₊,dz,tmp,b,k))

  QNDF1Cache(uprev2,du1,fsalfirst,k,z,dz,b,D,D2,R,U,tmp,atmp,utilde,J,
              W,uf,jac_config,linsolve,nlsolve,dtₙ₋₁)
end

# QNDF2

mutable struct QNDF2ConstantCache{F,N,coefType,coefType1,uType,dtType} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolve::N
  D::coefType
  D2::coefType
  R::coefType1
  U::coefType1
  uprev2::uType
  uprev3::uType
  dtₙ₋₁::dtType
  dtₙ₋₂::dtType
end

mutable struct QNDF2Cache{uType,rateType,coefType,coefType1,coefType2,uNoUnitsType,J,W,UF,JC,N,F,dtType} <: OrdinaryDiffEqMutableCache
  uprev2::uType
  uprev3::uType
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  b::uType
  D::coefType1
  D2::coefType2
  R::coefType
  U::coefType
  tmp::uType
  atmp::uNoUnitsType
  utilde::uType
  J::J
  W::W
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
  dtₙ₋₁::dtType
  dtₙ₋₂::dtType
end

u_cache(c::QNDF2Cache)  = ()
du_cache(c::QNDF2Cache) = ()

function alg_cache(alg::QNDF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  uprev2 = u
  uprev3 = u
  dtₙ₋₁ = zero(t)
  dtₙ₋₂ = zero(t)

  D = fill(zero(typeof(u)), 1, 2)
  D2 = fill(zero(typeof(u)), 1, 3)
  R = fill(zero(typeof(t)), 2, 2)
  U = fill(zero(typeof(t)), 2, 2)

  U!(2,U)

  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,zero(inv((1-alg.kappa))),1,ηold,z₊,dz,tmp,b,k))
  QNDF2ConstantCache(uf,nlsolve,D,D2,R,U,uprev2,uprev3,dtₙ₋₁,dtₙ₋₂)
end

function alg_cache(alg::QNDF2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  D = Array{typeof(u)}(undef, 1, 2)
  D2 = Array{typeof(u)}(undef, 1, 3)
  R = fill(zero(typeof(t)), 2, 2)
  U = fill(zero(typeof(t)), 2, 2)

  D[1] = similar(u); D[2] = similar(u)
  D2[1] = similar(u);  D2[2] = similar(u); D2[3] = similar(u)

  U!(2,U)

  atmp = similar(u,uEltypeNoUnits,axes(u))
  utilde = similar(u,axes(u))
  uprev2 = similar(u)
  uprev3 = similar(u)
  dtₙ₋₁ = zero(dt)
  dtₙ₋₂ = zero(dt)

  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,zero(inv((1-alg.kappa))),1,ηold,z₊,dz,tmp,b,k))
  QNDF2Cache(uprev2,uprev3,du1,fsalfirst,k,z,dz,b,D,D2,R,U,tmp,atmp,utilde,J,
             W,uf,jac_config,linsolve,nlsolve,dtₙ₋₁,dtₙ₋₂)
end

mutable struct QNDFConstantCache{F,N,coefType1,coefType2,coefType3,uType,uArrayType,dtType,dtsType} <: OrdinaryDiffEqConstantCache
  uf::F
  nlsolve::N
  D::coefType3
  D2::coefType2
  R::coefType1
  U::coefType1
  order::Int64
  max_order::Int64
  udiff::uArrayType
  dts::dtsType
  tmp::uType
  h::dtType
  c::Int64
end

mutable struct QNDFCache{uType,rateType,coefType1,coefType,coefType2,coefType3,dtType,dtsType,uNoUnitsType,J,W,UF,JC,N,F} <: OrdinaryDiffEqMutableCache
  du1::rateType
  fsalfirst::rateType
  k::rateType
  z::uType
  dz::uType
  b::uType
  D::coefType3
  D2::coefType2
  R::coefType1
  U::coefType1
  order::Int64
  max_order::Int64
  udiff::coefType
  dts::dtsType
  tmp::uType
  atmp::uNoUnitsType
  utilde::uType
  J::J
  W::W
  uf::UF
  jac_config::JC
  linsolve::F
  nlsolve::N
  h::dtType
  c::Int64
end

u_cache(c::QNDFCache)  = ()
du_cache(c::QNDFCache) = ()

function alg_cache(alg::QNDF,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  @oopnlcachefields
  udiff = fill(zero(typeof(u)), 1, 6)
  dts = fill(zero(typeof(dt)), 1, 6)
  h = zero(dt)
  tmp = zero(u)

  D = fill(zero(typeof(u)), 1, 5)
  D2 = fill(zero(typeof(u)), 6, 6)
  R = fill(zero(typeof(t)), 5, 5)
  U = fill(zero(typeof(t)), 5, 5)

  max_order = 5

  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,zero(eltype(alg.kappa)),1,ηold,z₊,dz,tmp,b,k))
  QNDFConstantCache(uf,nlsolve,D,D2,R,U,1,max_order,udiff,dts,tmp,h,0)
end

function alg_cache(alg::QNDF,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  @iipnlcachefields
  udiff = Array{typeof(u)}(undef, 1, 6)
  dts = fill(zero(typeof(dt)), 1, 6)
  h = zero(dt)

  D = Array{typeof(u)}(undef, 1, 5)
  D2 = Array{typeof(u)}(undef, 6, 6)
  R = fill(zero(typeof(t)), 5, 5)
  U = fill(zero(typeof(t)), 5, 5)

  for i = 1:5
    D[i] = zero(u)
    udiff[i] = zero(u)
  end
  udiff[6] = zero(u)

  for i = 1:6
    for j = 1:6
      D2[i,j] = zero(u)
    end
  end

  max_order = 5
  atmp = similar(u,uEltypeNoUnits,axes(u))
  utilde = similar(u,axes(u))
  nlsolve = typeof(_nlsolve)(NLSolverCache(κ,tol,min_iter,max_iter,10000,new_W,z,W,zero(eltype(alg.kappa)),1,ηold,z₊,dz,tmp,b,k))

  QNDFCache(du1,fsalfirst,k,z,dz,b,D,D2,R,U,1,max_order,udiff,dts,tmp,atmp,utilde,J,
            W,uf,jac_config,linsolve,nlsolve,h,0)
end


text:
using PRM
using Test

ext: jl

max_stars_repo_name: VeroAsghar/PRM

max_issues_repo_name: VeroAsghar/PRM

max_forks_repo_name: VeroAsghar/PRM


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

mutable struct SnapshotUpdateProperties <: SwaggerModel
    osType::Any # spec type: Union{ Nothing, String } # spec name: osType
    diskSizeGB::Any # spec type: Union{ Nothing, Int32 } # spec name: diskSizeGB
    encryptionSettingsCollection::Any # spec type: Union{ Nothing, EncryptionSettingsCollection } # spec name: encryptionSettingsCollection
    encryption::Any # spec type: Union{ Nothing, Encryption } # spec name: encryption
    networkAccessPolicy::Any # spec type: Union{ Nothing, NetworkAccessPolicy } # spec name: networkAccessPolicy
    diskAccessId::Any # spec type: Union{ Nothing, String } # spec name: diskAccessId

    function SnapshotUpdateProperties(;osType=nothing, diskSizeGB=nothing, encryptionSettingsCollection=nothing, encryption=nothing, networkAccessPolicy=nothing, diskAccessId=nothing)
        o = new()
        validate_property(SnapshotUpdateProperties, Symbol("osType"), osType)
        setfield!(o, Symbol("osType"), osType)
        validate_property(SnapshotUpdateProperties, Symbol("diskSizeGB"), diskSizeGB)
        setfield!(o, Symbol("diskSizeGB"), diskSizeGB)
        validate_property(SnapshotUpdateProperties, Symbol("encryptionSettingsCollection"), encryptionSettingsCollection)
        setfield!(o, Symbol("encryptionSettingsCollection"), encryptionSettingsCollection)
        validate_property(SnapshotUpdateProperties, Symbol("encryption"), encryption)
        setfield!(o, Symbol("encryption"), encryption)
        validate_property(SnapshotUpdateProperties, Symbol("networkAccessPolicy"), networkAccessPolicy)
        setfield!(o, Symbol("networkAccessPolicy"), networkAccessPolicy)
        validate_property(SnapshotUpdateProperties, Symbol("diskAccessId"), diskAccessId)
        setfield!(o, Symbol("diskAccessId"), diskAccessId)
        o
    end
end # type SnapshotUpdateProperties

const _property_map_SnapshotUpdateProperties = Dict{Symbol,Symbol}(Symbol("osType")=>Symbol("osType"), Symbol("diskSizeGB")=>Symbol("diskSizeGB"), Symbol("encryptionSettingsCollection")=>Symbol("encryptionSettingsCollection"), Symbol("encryption")=>Symbol("encryption"), Symbol("networkAccessPolicy")=>Symbol("networkAccessPolicy"), Symbol("diskAccessId")=>Symbol("diskAccessId"))
const _property_types_SnapshotUpdateProperties = Dict{Symbol,String}(Symbol("osType")=>"String", Symbol("diskSizeGB")=>"Int32", Symbol("encryptionSettingsCollection")=>"EncryptionSettingsCollection", Symbol("encryption")=>"Encryption", Symbol("networkAccessPolicy")=>"NetworkAccessPolicy", Symbol("diskAccessId")=>"String")
Base.propertynames(::Type{ SnapshotUpdateProperties }) = collect(keys(_property_map_SnapshotUpdateProperties))
Swagger.property_type(::Type{ SnapshotUpdateProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_SnapshotUpdateProperties[name]))}
Swagger.field_name(::Type{ SnapshotUpdateProperties }, property_name::Symbol) =  _property_map_SnapshotUpdateProperties[property_name]

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

function check_required(o::SnapshotUpdateProperties)
    true
end

function validate_property(::Type{ SnapshotUpdateProperties }, name::Symbol, val)
    if name === Symbol("osType")
        Swagger.validate_param(name, "SnapshotUpdateProperties", :enum, val, _allowed_SnapshotUpdateProperties_osType)
    end
end


text:
"""
Helping functions
"""

# Converts vector of words to vector of indices
function indices(wordVect::Vector, vocab::Vector, unk)
    function index(x, unk)
        idx = something(findfirst(isequal(x), vocab), 0)
        idx > 0 || return findfirst(isequal(unk), vocab)
        return idx
    end
    return broadcast(x -> index(x, unk), wordVect)
end

#Padding multiple sequences w r t the max size sequence
function pre_pad_sequences(sequences::Vector, pad::String="_pad_")
    max_len = maximum([length(x) for x in sequences])
    return [[fill(pad, max_len-length(sequence)); sequence] for sequence in sequences]
end

function post_pad_sequences(sequences::Vector, pad::String="_pad_")
    max_len = maximum([length(x) for x in sequences])
    return [[sequence; fill(pad, max_len-length(sequence))] for sequence in sequences]
end

# To initialize funciton for model LSTM weights
init_weights(extreme::AbstractFloat, dims...) = randn(Float32, dims...) .* sqrt(Float32(extreme))

# Generator, whenever it is called it gives one mini-batch
function generator(c::Channel, corpus; batchsize::Integer=64, bptt::Integer=70)
    X_total = post_pad_sequences(chunk(corpus, batchsize))
    n_batches = Int(floor(length(X_total[1])/bptt))
    put!(c, n_batches)
    for i=1:n_batches
        start = bptt*(i-1) + 1
        batch = [Flux.batch(X_total[k][j] for k=1:batchsize) for j=start:start+bptt]
        put!(c, batch[1:end-1])
        put!(c, batch[2:end])
    end
end

# Sequence Bucketing
# function sequence_bucketing(buckets::Integer, lengths::Vector{Integer}, k::Integer)
#     f = [count(x -> x == l, lengths) for l=1:maximum(lengths)]
#     n = length(f)
#     prev_dp = fill(fill(0, n-1), buckets)
#     dp = fill(convert(Array{Union{Int64, Float64}, 1}, fill(0, n)), buckets)
#     for i=1:buckets
#         dp[i][0] = Inf
#         (i-2 > 0) && (dp[i][k] = minimum([(dp[j][k-1] + i*(sum([f[t] for t=i+1:i]))) for j=1:i-2]))
#     end
#     for q=1:buckets
#         for i=1:n
#             cur_sum = f[i]
#             for j=i-1:1
#                 val = cur_sum*i + dp[q-1][j]
#                 if val < dp[q][i]
#                     dp[q][i] = val
#                     prev_dp[q][i] = j
#                 end
#                 cur_sum = cur_sum + f[j]####prob
#             end
#         end
#     end
#     cur_id = n-1
#     bests = []
#     for i=buckets:1
#         bests = append!([cur_id], bests)
#         cur_id = prev_dp[i][cur_id]
#     end
#     return bests
# end


text:
include("Way.jl")
include("Map.jl")
include("Solution.jl")

mutable struct Model

    map::Map
    solution::Solution

    ways::Dict{CityIndex, Dict{CityIndex, Way}}

end

function Model()
    return Model(Map(), Solution(), Dict{CityIndex, Dict{CityIndex, Way}}())
end

function Model(map::Map, pheroInit::Float64 = 1.)
    return Model(map, Solution(), createWays(map, pheroInit))
end

function Model(nbVille::Int64, borneX::Real = 50, borneY::Real = 50)
    map = Map(nbVille, borneX, borneY)
    return Model(map)
end

function updatePhero!(model::Model, p::Float64)
    for dictWays in collect(values(model.ways))
        for way in collect(values(dictWays))
            way.pheromone *= p
        end
    end
    return map
end


text:
export Window,
       include_asset

immutable Window
    dimension::Input
    route::Input
    dir::String
    assets::Input
end

Window(;
    dimension=(0mm, 0mm),
    route="",
    dir="ltr") =
    Window(Input{Any}(dimension), Input{Any}(route), "ltr", Input("basics"))

function resolve_asset(slug)
    path = Pkg.dir("Canvas", "assets", slug * ".html")
    if isfile(path)
        return "/assets/$slug.html"
    else
        error("Asset file $path doesn't exist")
    end
end


text:
#  Tests for the projected NEPs

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra
using Random

#using Winston # For plotting

@testset "Projected problems" begin

    @bench @testset "Problem $nepstr" for nepstr in ("pep", "dep", "sqrtm")

        local nep::NEP
        if (nepstr == "pep")
            nep=nep_gallery("pep0_sparse");
        elseif (nepstr == "dep")

            n=5;
            Random.seed!(1)
            A0=(1:5)*(1:5)'
            A1=(1:5)*(3:7)'+I;
            t::Float64=3.0

            minusop= S-> -S
            oneop= S -> one(S)
            expmop= S -> exp(-t*S)
            fi=[minusop, oneop, expmop];

            nep = SPMF_NEP([Matrix(1.0I, n, n), A0, A1], fi)

        elseif (nepstr == "sqrtm")

            n=8;
            Random.seed!(1)
            A1=I+(1:n)*(1:n)'/n
            A0=(1:n)*(3:(n+2))'+2I;
            A2=(-1:(n-2))*(3:(n+2))'/8-I
            t=3.0

            minusop= S-> -S
            oneop= S -> one(S)  # one(S) replaced by I it would not return a matrix
            expmop= S -> sqrt(-t*S + 30*one(S))
            fi=[minusop, oneop, expmop];

            nep=SPMF_NEP([A0,A1,A2],fi)
        end

        #

        n=size(nep,1);
        @info "Running Newton Raphson ($nepstr)"
        tol=1e-12
        λ,x =newton(nep,maxit=30,λ=1+1im,tol=tol,
                    v=ones(n));
        # Check residual is small
        λ_exact=λ
        @test norm(compute_Mlincomb(nep,λ,x))<tol*1000;

        # Create a projected NEP
        pnep=create_proj_NEP(nep,4); # maxsize=4
        V=(1:n)*(1:2)'/n; V[1,1]=pi;
        Q,R=qr(hcat(V,x)) # Make the eigenspace a part of the projection subspace
        Q = Matrix(Q)
        set_projectmatrices!(pnep,Q,Q);
        @info "Running Newton on projected problem with very good start value ($nepstr)"
        λ1,z1=newton(pnep,λ=(λ_exact+0.00001),logger=0,v=ones(size(pnep,1)))

        x1=Q*z1; x1=x1/x1[1];

        # should be small since the eigenvector is in the subspace
        @test norm(x/x[1]-x1)<1e-8

        @info "Running IAR for projected problem ($nepstr)"
        λv,X=iar(pnep,σ=complex(round(λ_exact*10)/10),
                 neigs=3,maxit=100, v=ones(size(pnep,1)))

        mydiff = minimum(abs.(λv .- λ_exact))
        @info "Difference of solution from projected problem ($nepstr): $mydiff"
        @test mydiff < sqrt(eps())

        @info "Testing to expand and running Newton ($nepstr)"
        n=size(nep,1);
        Vnew=[Q ones(n)];
        Wnew=[Q ones(n)];
        expand_projectmatrices!(pnep,Wnew,Vnew);
        λ1,zz1=newton(pnep,λ=(λ_exact+0.0000001),logger=0,
                      v=Vnew'*x.+0.00001*ones(size(pnep,1)),maxit=30)

        xx1=Vnew*zz1; xx1=xx1/xx1[1]

        # should be small since the eigenvector is in the subspace
        @test norm(x/x[1]-x1)<1e-8

    end
end


text:
"""
Shunt constraint using linearized voltage magnitude difference phi

```
sum(p) + sum(pconv_grid_ac)  == sum(pg) - sum(pd) - sum(gs*(1.0 + 2*phi)
```
"""
function constraint_power_balance_ac(pm::_PM.AbstractLPACModel, n::Int,  i::Int, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_loads, bus_shunts, pd, qd, gs, bs)
    phi = _PM.var(pm, n, :phi, i)
    p = _PM.var(pm, n, :p)
    q = _PM.var(pm, n, :q)
    pg = _PM.var(pm, n, :pg)
    qg = _PM.var(pm, n, :qg)
    pconv_grid_ac = _PM.var(pm, n, :pconv_tf_fr)
    qconv_grid_ac = _PM.var(pm, n, :qconv_tf_fr)

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(pconv_grid_ac[c] for c in bus_convs_ac)  == sum(pg[g] for g in bus_gens)   - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*(1.0 + 2*phi))
    JuMP.@constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(qconv_grid_ac[c] for c in bus_convs_ac)  == sum(qg[g] for g in bus_gens)   - sum(qd[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*(1.0 + 2*phi))
end
"""
DC branch power flow using linearized voltage magnitude difference phi

```
p_dc_fr == p * g *  (phi_fr - phi_to)
p_dc_to == p * g *  (phi_to - phi_fr)
```
"""
function constraint_ohms_dc_branch(pm::_PM.AbstractLPACModel, n::Int,  f_bus, t_bus, f_idx, t_idx, r, p)
    p_dc_fr = _PM.var(pm, n, :p_dcgrid, f_idx)
    p_dc_to = _PM.var(pm, n, :p_dcgrid, t_idx)
    phi_fr = _PM.var(pm, n, :phi_vdcm, f_bus)
    phi_to = _PM.var(pm, n, :phi_vdcm, t_bus)
    phi_fr_ub = JuMP.UpperBoundRef(phi_to)
    phi_fr_lb = JuMP.LowerBoundRef(phi_to)

    if r == 0
        JuMP.@constraint(pm.model, p_dc_fr + p_dc_to == 0)
    else
        g = 1 / r
        JuMP.@constraint(pm.model, p_dc_fr == p * g *  (phi_fr - phi_to))
        JuMP.@constraint(pm.model, p_dc_to == p * g *  (phi_to - phi_fr))
    end
end

function constraint_dc_branch_current(pm::_PM.AbstractLPACModel, n::Int,  f_bus, f_idx, ccm_max, p)

end

function add_dc_bus_voltage_setpoint(sol, pm::_PM.AbstractLPACModel)
    _PM.add_setpoint!(sol, pm, "busdc", "vm", :phi_vdcm, status_name="Vdc", inactive_status_value = 4, scale = (x,item,cnd) -> 1.0+x)
end

############# TNEP Constraints #################

function constraint_power_balance_ac_dcne(pm::_PM.AbstractLPACModel, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, pd, qd, gs, bs)
    phi = _PM.var(pm, n, :phi, i)
    p = _PM.var(pm, n, :p)
    q = _PM.var(pm, n, :q)
    pg = _PM.var(pm, n, :pg)
    qg = _PM.var(pm, n, :qg)
    pconv_grid_ac = _PM.var(pm, n, :pconv_tf_fr)
    qconv_grid_ac = _PM.var(pm, n, :qconv_tf_fr)
    pconv_grid_ac_ne = _PM.var(pm, n, :pconv_tf_fr_ne)
    qconv_grid_ac_ne = _PM.var(pm, n, :qconv_tf_fr_ne)

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(pg[g] for g in bus_gens)  - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*(1.0 + 2*phi))
    JuMP.@constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(qconv_grid_ac[c] for c in bus_convs_ac) + sum(qconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(qg[g] for g in bus_gens)  - sum(qd[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*(1.0 + 2*phi))
end

function constraint_power_balance_acne_dcne(pm::_PM.AbstractLPACModel, n::Int, i::Int, bus_arcs, bus_arcs_ne, bus_arcs_dc, bus_gens, bus_convs_ac, bus_convs_ac_ne, bus_loads, bus_shunts, pd, qd, gs, bs)
    phi = _PM.var(pm, n, :phi, i)
    p = _PM.var(pm, n, :p)
    p_ne = _PM.var(pm, n, :p_ne)
    q = _PM.var(pm, n, :q)
    q_ne = _PM.var(pm, n, :q_ne)
    pg = _PM.var(pm, n, :pg)
    qg = _PM.var(pm, n, :qg)
    pconv_grid_ac = _PM.var(pm, n, :pconv_tf_fr)
    qconv_grid_ac = _PM.var(pm, n, :qconv_tf_fr)
    pconv_grid_ac_ne = _PM.var(pm, n, :pconv_tf_fr_ne)
    qconv_grid_ac_ne = _PM.var(pm, n, :qconv_tf_fr_ne)

    JuMP.@constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_ne[a] for a in bus_arcs_ne) + sum(pconv_grid_ac[c] for c in bus_convs_ac) + sum(pconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(pg[g] for g in bus_gens)  - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*(1.0 + 2*phi))
    JuMP.@constraint(pm.model, sum(q[a] for a in bus_arcs) + sum(q_ne[a] for a in bus_arcs_ne) + sum(qconv_grid_ac[c] for c in bus_convs_ac) + sum(qconv_grid_ac_ne[c] for c in bus_convs_ac_ne)  == sum(qg[g] for g in bus_gens)  - sum(qd[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*(1.0 + 2*phi))
end

function constraint_ohms_dc_branch_ne(pm::_PM.AbstractLPACModel, n::Int, f_bus, t_bus, f_idx, t_idx, r, p)
    l = f_idx[1];
    p_dc_fr = _PM.var(pm, n, :p_dcgrid_ne, f_idx)
    p_dc_to = _PM.var(pm, n, :p_dcgrid_ne, t_idx)
    phi_fr = []
    phi_to = []
    phi_fr_du = _PM.var(pm, n, :phi_vdcm_fr, l)
    phi_to_du = _PM.var(pm, n, :phi_vdcm_to, l)
    z = _PM.var(pm, n, :branchdc_ne, l)
    phi_to, phi_fr = contraint_ohms_dc_branch_busvoltage_structure_phi(pm, n, f_bus, t_bus, phi_to, phi_fr)
    if r == 0
        JuMP.@constraint(pm.model, p_dc_fr + p_dc_to == 0)
    else
        g = 1 / r
        JuMP.@constraint(pm.model, p_dc_fr == p * g *  (phi_fr_du - phi_to_du))
        JuMP.@constraint(pm.model, p_dc_to == p * g *  (phi_to_du - phi_fr_du))

        JuMP.@constraint(pm.model, phi_to_du <= phi_to - JuMP.lower_bound(phi_to)*(1-z))
        JuMP.@constraint(pm.model, phi_to_du >= phi_to - JuMP.upper_bound(phi_to)*(1-z))
        JuMP.@constraint(pm.model, phi_fr_du <= phi_fr - JuMP.lower_bound(phi_fr)*(1-z))
        JuMP.@constraint(pm.model, phi_fr_du >= phi_fr - JuMP.upper_bound(phi_fr)*(1-z))

        relaxation_variable_on_off(pm.model, phi_to, phi_to_du, z)
        relaxation_variable_on_off(pm.model, phi_fr, phi_fr_du, z)

    end
end

function constraint_dc_branch_current_ne(pm::_PM.AbstractLPACModel, n::Int, f_bus, f_idx, ccm_max, p)

end

function add_dc_bus_voltage_setpoint_ne(sol, pm::_PM.AbstractLPACModel)
    _PM.add_setpoint!(sol, pm, "busdc_ne", "vm", :phi_vdcm_ne, status_name="Vdc", inactive_status_value = 4, scale = (x,item,cnd) -> 1.0+x)
end

function contraint_ohms_dc_branch_busvoltage_structure_phi(pm::_PM.AbstractPowerModel, n::Int, f_bus, t_bus, phi_to, phi_fr)
    for i in _PM.ids(pm, n, :busdc_ne)
        if t_bus == i
            phi_to = _PM.var(pm, n, :phi_vdcm_ne, t_bus)
        end
        if f_bus == i
            phi_fr = _PM.var(pm, n, :phi_vdcm_ne, f_bus)
        end
    end
    for i in _PM.ids(pm, n, :busdc)
        if t_bus == i
            phi_to = _PM.var(pm, n, :phi_vdcm, t_bus)
        end
        if f_bus == i
            phi_fr = _PM.var(pm, n, :phi_vdcm, f_bus)
        end
    end
    return phi_to, phi_fr
end


text:
let
    println("testing: GMM-based parameter conversion")
    modelpath = joinpath(Pkg.dir("VoiceConversion"), "test", "models",
                         "clb_to_slt_gmm32_order40_diff.jld")
    gmm = load(modelpath)
    @assert gmm["diff"]

    mapper = GMMMap(gmm["weights"], gmm["means"], gmm["covars"])
    @test length(mapper) == 1

    D = div(size(gmm["means"], 1), 2)
    @test dim(mapper) == D
    @test ncomponents(mapper) == length(gmm["weights"])
    @test size(mapper) == (D, 1)
end


text:
using Test, CSV, Dates, Tables, DataFrames, CategoricalArrays, PooledArrays, CodecZlib, FilePathsBase

const dir = joinpath(dirname(pathof(CSV)), "..", "test", "testfiles")

@eval macro $(:try)(ex)
    quote
        try $(esc(ex))
        catch
        end
    end
end

@testset "CSV" begin

@testset "CSV.File" begin

include("basics.jl")
include("testfiles.jl")
include("iteration.jl")

end # @testset "CSV.File"

include("write.jl")

@testset "transform" begin
    csv = """A,B,C
    1,1.0,"hey"
    ,2.0,"there"
    3,3.0,"sailor"
    """

    df = CSV.read(IOBuffer(csv)) |> Tables.transform(Dict{String, Base.Callable}("C" => Symbol)) |> DataFrame
    @test df.C == [:hey, :there, :sailor]

    df = CSV.read(IOBuffer(csv)) |> Tables.transform(Dict("A"=>x->x+1)) |> DataFrame
    @test isequal(df.A, [2, missing, 4])
    @test typeof(df.A) == Vector{Union{Missing, Int64}}
    @test size(df) == (3, 3)

    df = CSV.read(IOBuffer(csv)) |> Tables.transform(Dict("A"=>x->coalesce(x+1, 0))) |> DataFrame
    @test df.A == [2, 0, 4]
    @test eltype(df.A) <: Signed

    df = CSV.read(IOBuffer(csv)) |> Tables.transform(Dict("A"=>x->coalesce(x+1, 0.0))) |> DataFrame
    @test df.A == [2, 0.0, 4]
    @test eltype(df.A) <: Real

    df = CSV.read(IOBuffer(csv)) |> Tables.transform(Dict(2=>x->x==2.0 ? missing : x)) |> DataFrame
    @test size(df) == (3, 3)
    @test isequal(df.B, [1.0, missing, 3.0])
    @test typeof(df.B) == Vector{Union{Float64, Missing}}

    f = CSV.File(IOBuffer(csv))

    # transforms on CSV.Columns
    df = f |> Tables.transform(A=x->x+1) |> DataFrame
    @test size(df) == (3, 3)
    df = f |> Tables.transform(A=x->x+1) |> columntable
    @test length(df) == 3
    df = f |> Tables.transform(A=x->x+1) |> rowtable
    @test length(df) == 3

    # transforms on CSV.RowIterator
    f = CSV.File(joinpath(dir, "pandas_zeros.csv"));
    df = f |> Tables.transform(_0=x->x+1) |> DataFrame;
    @test size(df) == (100000, 50)
    df = f |> Tables.transform(_0=x->x+1) |> columntable;
    @test length(df) == 50
    # df = f |> Tables.transform(_0=x->x+1) |> rowtable;
    # @test length(df) == 100000
end

@testset "CategoricalArray levels (including ordering)" begin
    f = CSV.read(IOBuffer("X\nb\nc\na\nc"), types=[CategoricalString{UInt32}], copycols=true)
    v = f.X[1]
    @test v == "b"
    @test levels(v.pool) == ["a", "b", "c"]

    f = CSV.read(IOBuffer("X\nb\nc\na\nc"), categorical=true, copycols=true)
    v = f.X[1]
    @test v == "b"
    @test levels(v.pool) == ["a", "b", "c"]

    f = CSV.read(IOBuffer("X\nb\nc\n\nc"), categorical=true, copycols=true)
    v = f.X[1]
    @test v == "b"
    @test levels(v.pool) == ["b", "c"]
    @test typeof(f.X) == CategoricalArray{Union{Missing, String},1,UInt32,String,CategoricalString{UInt32},Missing}

end

@testset "PooledArrays" begin

    df = CSV.read(IOBuffer("X\nb\nc\na\nc"), pool=true, copycols=true)
    @test typeof(df.X) == PooledArrays.PooledArray{String,UInt32,1,Array{UInt32,1}}
    @test size(df) == (4, 1)
    @test df.X == ["b", "c", "a", "c"]
    @test df.X.refs[2] == df.X.refs[4]

    df = CSV.read(IOBuffer("X\nb\nc\na\nc"), pool=0.6, copycols=true)
    @test typeof(df.X) == PooledArrays.PooledArray{String,UInt32,1,Array{UInt32,1}}
    @test size(df) == (4, 1)
    @test df.X == ["b", "c", "a", "c"]
    @test df.X.refs[2] == df.X.refs[4]

    df = CSV.read(IOBuffer("X\nb\nc\n\nc"), pool=true, copycols=true)
    @test typeof(df.X) == PooledArray{Union{Missing, String},UInt32,1,Array{UInt32,1}}
    @test size(df) == (4, 1)
    @test df.X[3] === missing

    df = CSV.read(IOBuffer("X\nc\nc\n\nc\nc\nc\nc\nc\nc"), copycols=true)
    @test typeof(df.X) == PooledArray{Union{Missing, String},UInt32,1,Array{UInt32,1}}
    @test size(df) == (9, 1)
    @test isequal(df.X, ["c", "c", missing, "c", "c", "c", "c", "c", "c"])

end

@testset "CSV.Rows" begin

    rows = CSV.Rows(IOBuffer("X\nb\nc\na\nc"))
    show(rows)
    row = first(rows)
    @test row.X == "b"
    @test row[1] == "b"

    # limit
    rows = collect(CSV.Rows(IOBuffer("X\nb\nc\na\nc"), limit=1))
    @test length(rows) == 1
    @test rows[1][1] == "b"
    rows = collect(CSV.Rows(IOBuffer("X\nb\nc\na\nc"), limit=0))
    @test length(rows) == 0

    # transpose
    rows = collect(CSV.Rows(IOBuffer("x,1\nx2,2\n"), transpose=true))
    @test length(rows) == 1
    @test rows[1].x == "1"
    @test rows[1].x2 == "2"

    # ignorerepeated
    rows = collect(CSV.Rows(IOBuffer("x   y   z\n   1   2   3"), ignorerepeated=true, delim=' '))
    @test length(rows) == 1
    @test length(rows[1]) == 3
    @test rows[1] == ["1", "2", "3"]
    row = rows[1]
    @test CSV.Parsers.parse(Int, row, 1) == 1
    @test CSV.Parsers.parse(Int, row, :x) == 1
    @test CSV.detect(row, 2) == 2
    @test CSV.detect(row, :y) == 2

    # 448
    @test_throws ArgumentError CSV.Rows(IOBuffer("x\n1\n2\n3\n#4"), ignorerepeated=true)

    # 447
    rows = collect(CSV.Rows(IOBuffer("a,b,c\n1,2,3\n\n"), ignoreemptylines=true))
    @test length(rows) == 1
    @test rows[1] == ["1", "2", "3"]

    # not enough columns
    rows = collect(CSV.Rows(IOBuffer("x,y,z\n1\n2,3\n4,5,6")))
    @test length(rows) == 3
    @test isequal(rows[1], ["1", missing, missing])
    @test isequal(rows[2], ["2", "3", missing])
    @test rows[3] == ["4", "5", "6"]

    # too many columns
    rows = collect(CSV.Rows(IOBuffer("x,y,z\n1,2,3\n4,5,6,7,8,")))
    @test length(rows) == 2
    @test rows[1] == ["1", "2", "3"]
    @test rows[2] == ["4", "5", "6"]

    # fatal
    @test_throws CSV.Error collect(CSV.Rows(IOBuffer("x\n\"invalid quoted field")))

    # reusebuffer
    rows = collect(CSV.Rows(IOBuffer("x\n1\n2\n3")))
    @test length(rows) == 3
    @test rows[1] == ["1"]
    @test rows[2] == ["2"]
    @test rows[3] == ["3"]
    rows = collect(CSV.Rows(IOBuffer("x\n1\n2\n3"), reusebuffer=true))
    @test length(rows) == 3
    @test rows[1] == ["3"]
    @test rows[2] == ["3"]
    @test rows[3] == ["3"]

    for (i, row) in enumerate(CSV.Rows(IOBuffer("x\n1\n2\n3"), reusebuffer=true))
        @test row == [string(i)]
    end

end

@testset "CSV.detect" begin

@test CSV.detect("") == ""
@test CSV.detect("1") == 1
@test CSV.detect("1.1") == 1.1
@test CSV.detect("2015-01-01") == Date(2015)
@test CSV.detect("2015-01-01T03:04:05") == DateTime(2015, 1, 1, 3, 4, 5)
@test CSV.detect("true") === true
@test CSV.detect("false") === false
@test CSV.detect("abc") === "abc"

end

@testset "CSV.findrowstarts!" begin

rngs = [1, 1, 1]
buf = b"normal cell,next cell\nnormal cell2,next cell2\nhey"
CSV.findrowstarts!(buf, length(buf), CSV.Parsers.XOPTIONS, nothing, false, rngs, 2)
@test rngs[2] == 23

rngs = [1, 1, 1]
buf = b"quoted, cell\",next cell\n\"normal cell2\",next cell2\nhey"
CSV.findrowstarts!(buf, length(buf), CSV.Parsers.XOPTIONS, nothing, false, rngs, 2)
@test rngs[2] == 25

rngs = [1, 2, 1]
buf = b"\"\"quoted, cell\",next cell\n\"normal cell2\",next cell2\nhey"
CSV.findrowstarts!(buf, length(buf), CSV.Parsers.XOPTIONS, nothing, false, rngs, 2)
@test rngs[2] == 27

rngs = [1, 2, 1]
buf = b"quoted,\"\" cell\",next cell\n\"normal cell2\",next cell2\nhey"
CSV.findrowstarts!(buf, length(buf), CSV.Parsers.XOPTIONS, nothing, false, rngs, 2)
@test rngs[2] == 27

end

end


text:
doc"""
    typical(ohlc; h="High", l="Low", c="Close")

Typical Price

```math
    \text{Typical Price} = \frac{H + L + C}{3}
```
"""
function typical(ohlc::TimeArray{T,N}; h="High", l="Low", c="Close") where {T,N}
    val = (ohlc[h] .+ ohlc[l] .+ ohlc[c]) ./ 3
    TimeArray(val.timestamp, val.values, ["typical"], ohlc.meta)
end

function mean_abs_dev(a::Array{T,1}, scale::Bool=false) where T
    scale ? c = 1 / (-sqrt(2) * erfcinv(3 * 1 / 2)) : c =1

    res = ones(length(a))
    for i in 1:length(a)
        @inbounds res[i] = abs(a[i] - mean(a))
    end

    mean(res) * c
end

"""
Parameters:
    * orig: the original colnames
    * suffix: the suffix applied to `orig`

```jldoctest
julia> gen_colnames(["Open", "Close"], ["macd", "dif", "sig"])
6-element Array{String,1}:
 "Open_macd"
 "Close_macd"
 "Open_dif"
 "Close_dif"
 "Open_sig"
 "Close_sig"
```
"""
gen_colnames(orig::Vector{String}, suffix::Vector{String}) =
    vec(["$o\_$s" for o ∈ orig, s ∈ suffix])

relu(x) = max(x, 0)

relu(ta::TimeArray) =
    TimeArray(ta.timestamp, relu.(ta.values), ta.colnames, ta.meta)

"""
The soomth method used by ADX
"""
function wilder_smooth(ta::TimeArray, n::Integer;
                       padding::Bool=false, dx::Bool=false)
    val = similar(Array{Float64}, indices(ta.values))

    first_cal = (dx ? mean : sum)

    for i ∈ 1:size(val, 1)
        val[i, :] =
            if i < n
                NaN
            elseif i == n
                first_cal(ta[1:n].values, 1)
            elseif dx
                (val[i-1, :] .* (n - 1) .+ ta.values[i, :]) ./ n
            else
                (val[i-1, :] .* (n - 1) ./ n) .+ ta.values[i, :]
            end
    end

    ret = TimeArray(ta.timestamp, val, ta.colnames, ta.meta)

    if padding
        ret
    else
        dropnan(ret)
    end
end

safediv(x, y) = ifelse(iszero(x) && iszero(y), x, x / y)


text:
using RCall

R"
library(gridExtra)
library(cowplot)
library(tidyverse)
library(latex2exp)
library(viridis)
"

load("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_7_1.jld")

R"
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
df_r= $DF_PY
df_r$beta_scaled= range01(df_r$beta_log)
df_r$color = rgb(df_r$beta_scaled, df_r$sigma,0.5,1)
#df_r$color = rgb(0, df_r$sigma,0.5,1)
#df_r$color_sigma = rgb(0, df_r$sigma,0.5,1)
df_r$color_sigma = rgb(df_r$beta_scaled, df_r$sigma,0.5,1)
df_r$color_beta = rgb(0, 0,0.5,1)
p <- df_r %>% ggplot(aes(x=Exp, y = std, group=sigma))  + geom_line(alpha = 0.8,linetype = 'longdash') + geom_point(aes(colour=color), size=4)
p_py <- p  + geom_line(data =df_r, aes(x=Exp, y = std, group=beta), alpha = 0.8, linetype='dotted') +
  xlim(0, 100)+ ylim(0,16)+ labs(y='Std', x='Expectation')+scale_colour_identity()+theme_classic()+
   theme(plot.title = element_text(hjust = 0.5,size = 15), axis.text.x = element_text(size=15), axis.text.y = element_text(size=15),legend.position='none')
"

R"p_py
pdf('/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/figures/figure_7_1.pdf')
plot(p_py)
dev.off()"

load("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_7_2.jld")

R"
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
df_r= $DF_PYM
df_r$beta_scaled= range01(df_r$beta_log)
df_r$color = rgb(df_r$beta_scaled, df_r$sigma,0.5,1)
print(df_r[1:15,])
df_r$color_sigma = rgb(0, df_r$sigma,0.5,1)
df_r$color_beta = rgb(df_r$beta_scaled, 0,0.5,1)
p <- df_r %>% ggplot(aes(x=Exp, y = std, group=sigma))  + geom_line(alpha = 0.8,linetype = 'longdash') + geom_point(aes(colour=color), size=4)
p_pym <- p  + geom_line(data =df_r, aes(x=Exp, y = std, group=beta), alpha = 0.8, linetype='dotted') +
  xlim(0, 100)+ ylim(0,16)+ labs(y='Std', x='Expectation')+scale_colour_identity()+theme_classic()+
   theme(plot.title = element_text(hjust = 0.5,size = 15), axis.text.x = element_text(size=15), axis.text.y = element_text(size=15),legend.position='none')
"

R"p_pym
pdf('/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/figures/figure_7_2.pdf')
plot(p_pym)
dev.off()"

load("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_7_3.jld")

R"
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
df_r= $DF_PY_apr
df_r$beta_scaled= range01(df_r$beta_log)
df_r$color = rgb(df_r$beta_scaled, df_r$sigma,0.5,1)
print(df_r[1:15,])
df_r$color_sigma = rgb(0, df_r$sigma,0.5,1)
df_r$color_beta = rgb(df_r$beta_scaled, 0,0.5,1)
p <- df_r %>% ggplot(aes(x=Exp, y = std, group=sigma))  + geom_line(alpha = 0.8,linetype = 'longdash') + geom_point(aes(colour=color), size=4)
p_py_apr <- p  + geom_line(data =df_r, aes(x=Exp, y = std, group=beta), alpha = 0.8, linetype='dotted') +
  xlim(0, 100)+ ylim(0,16)+ labs(y='Std', x='Expectation')+scale_colour_identity()+theme_classic()+
   theme(plot.title = element_text(hjust = 0.5,size = 15), axis.text.x = element_text(size=15), axis.text.y = element_text(size=15),legend.position='none')
"

R"p_py_apr
pdf('/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/figures/figure_7_3.pdf')
plot(p_py_apr)
dev.off()"

load("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_7_4.jld")

R"
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
df_r= $DF_PY_sb
df_r$beta_scaled= range01(df_r$beta_log)
df_r$color = rgb(df_r$beta_scaled, df_r$sigma,0.5,1)
print(df_r[1:15,])
df_r$color_sigma = rgb(0, df_r$sigma,0.5,1)
df_r$color_beta = rgb(df_r$beta_scaled, 0,0.5,1)
p <- df_r %>% ggplot(aes(x=Exp, y = std, group=sigma))  + geom_line(alpha = 0.8,linetype = 'longdash') + geom_point(aes(colour=color), size=4)
p_p_sb <- p  + geom_line(data =df_r, aes(x=Exp, y = std, group=beta), alpha = 0.8, linetype='dotted') +
  xlim(0, 100)+ ylim(0,16)+ labs(y='Std', x='Expectation')+scale_colour_identity()+theme_classic()+
   theme(plot.title = element_text(hjust = 0.5,size = 15), axis.text.x = element_text(size=15), axis.text.y = element_text(size=15),legend.position='none')
"

R"p_p_sb
pdf('/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/figures/figure_7_4.pdf')
plot(p_p_sb)
dev.off()"

#### color grid

load("/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/saves_for_figures/figure_7_5.jld")

R"
df_color=data.frame(sigma= $(last.(grid_vec)), beta_log=  log($(first.(grid_vec))))
df_color$beta_scaled= range01(df_color$beta_log)
df_color$color= rgb(df_color$beta_scaled,df_color$sigma,0.5,1)
df_color$color_sigma= rgb(0,df_color$sigma,0.5,1)
df_color$color_beta= rgb(df_color$beta_scaled,0,0.5,1)
p <- df_color %>% ggplot(aes(x=beta_log, y = sigma, group=sigma))  + geom_line(aes(color = color_sigma), alpha = 0.8,linetype = 'longdash') + geom_point(aes(colour=color), size=3)
p_color<- p +  geom_line(data =df_color, aes(x=beta_log, y = sigma, group=beta_log), alpha = 0.8, linetype='dotted') +theme_classic()+theme(legend.position='none')
p_color
"

R"p_color
pdf('/Users/dariabystrova/Documents/GitHub/approximatingBNPpriors/figures/figure_7_5.pdf')
plot(p_color)
dev.off()"

#############

R"
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
#df_r$beta_scaled= range01(df_r$beta_log)
#df_r$color = rgb(df_r$beta_scaled, df_r$sigma,0.5,1)

betas = exp(seq(log(1), log(200), length.out=10))
sigmas = seq(0.05, 0.99, length.out = 10)

library(tidyverse)

to_plot = expand_grid(betas, sigmas) %>%
  mutate(col = rgb(red = range01(log(betas)), green = sigmas, blue = 0.5, alpha = 1))

pt = to_plot %>%
  ggplot(aes(x = betas, y = sigmas)) +
  geom_point(aes(colour = col),size=4) + scale_x_log10() + annotation_logticks(base = 10)+ scale_color_identity()
 pt = pt +  geom_line(data =to_plot, aes(x=betas, y = sigmas, group=betas),linetype = 'longdash', alpha = 0.5,size=0.3)+
   geom_line(data =to_plot, aes(x=betas, y = sigmas, group=sigmas),linetype = 'dotted', alpha = 0.5,size=0.3)+ theme_classic()+
   theme(plot.title = element_text(hjust = 0.5,size = 15), axis.text = element_text(size=15), axis.title=element_text(size=15),legend.position='none')+
   ylab(TeX(sprintf('$\\sigma$')))+ xlab(TeX(sprintf('$\\tau$')))"

  R"
  pdf(file = '/Users/bystrova/Documents/GitHub/GibbsTypePriors/test/expectation_variance/Std_variance_NGGs_color_grid_v4.pdf',width= 4, height = 4)
  plot(pt)
  dev.off()
  "

max_stars_repo_name: konkam/approximatingBNPpriors

max_issues_repo_name: konkam/approximatingBNPpriors

max_forks_repo_name: konkam/approximatingBNPpriors


text:
using ThreeBodyProblem
using DifferentialEquations
using Plots

μ₁ = 398600 # {km³/s²} gravitational parameter of Earth
μ₂ = 4903   # {km³/s²} gravitational parameter of the Moon
d = 384400  # {km} average distance between Earth and the Moon
p = [μ₁, μ₂, d];

R₁,R₂ = computeR1R2(p)
ωₛ = sqrt((μ₁ + μ₂)/d^3);

Rₑ = 6378.0 # {km} radius of the Earth
Rₘ = 1738.0 # {km} radius of the Moon
h = 200.0   # {km} altitude of parking orbit

L1 = computeL1(p)

vᵢ₀ = 10.92367104   # {km/s} synodic frame injection velocity
ϕ₀ = 47.70061087    # {°}

vᵢ = vᵢ₀
ϕ = ϕ₀

r₀ = [-R₁ - (Rₑ + 200)*cosd(ϕ); -(Rₑ + 200)*sind(ϕ); 0]
v₀ = vᵢ*[sind(ϕ); -cosd(ϕ); 0];
rv₀ = [r₀;v₀]

tspan = [0.,86400*3] # 3 days

prob = ODEProblem(CR3BPdynamics!,rv₀,tspan,p)
sol = solve(prob,reltol=1e-6)

plot(sol,vars=(1,2),title="Initial Guess",label="")
# scatter!(L1[1],L1[2],marker="x")
plot!(circle(Rₑ,[-R₁;0]),color="blue",label="Earth")
plot!(circle(Rₘ,[R₂;0]),color="gray",label="Moon",aspect_ratio=1)


text:
"""
Slice an AbstractMesh at the specified Z axis value.
Returns a Vector of LineSegments generated from the faces at the specified
heights. Note: This will not slice in-plane faces.
"""
function Base.slice{VT<:AbstractFloat,FT<:Integer,O}(mesh::AbstractMesh{Point{3,VT},Face{3,FT,O}}, height::Number)

    height_ct = length(height)
    # intialize the LineSegment array
    slice = Simplex{2,Point{2,VT}}[]

    for face in mesh.faces
        v1,v2,v3 = mesh.vertices[face]
        zmax = max(v1[3], v2[3], v3[3])
        zmin = min(v1[3], v2[3], v3[3])
        if height > zmax
            continue
        elseif zmin <= height
            if v1[3] < height && v2[3] >= height && v3[3] >= height
                p1 = v1
                p2 = v3
                p3 = v2
            elseif v1[3] > height && v2[3] < height && v3[3] < height
                p1 = v1
                p2 = v2
                p3 = v3
            elseif v2[3] < height && v1[3] >= height && v3[3] >= height
                p1 = v2
                p2 = v1
                p3 = v3
            elseif v2[3] > height && v1[3] < height && v3[3] < height
                p1 = v2
                p2 = v3
                p3 = v1
            elseif v3[3] < height && v2[3] >= height && v1[3] >= height
                p1 = v3
                p2 = v2
                p3 = v1
            elseif v3[3] > height && v2[3] < height && v1[3] < height
                p1 = v3
                p2 = v1
                p3 = v2
            else
                continue
            end

            start = Point{2,VT}(p1[1] + (p2[1] - p1[1]) * (height - p1[3]) / (p2[3] - p1[3]),
                                p1[2] + (p2[2] - p1[2]) * (height - p1[3]) / (p2[3] - p1[3]))
            finish = Point{2,VT}(p1[1] + (p3[1] - p1[1]) * (height - p1[3]) / (p3[3] - p1[3]),
                                 p1[2] + (p3[2] - p1[2]) * (height - p1[3]) / (p3[3] - p1[3]))

            push!(slice, Simplex{2,Point{2,VT}}(start, finish))
        end
    end

    return slice
end


text:
module Pairwise

using DoubleUAVs, QMDP, POMDPs, GridInterpolations, Interact,
      PGFPlots, JLD

export get_pomdp, gen_pairwise_policy, viz_pairwise_policy

const ALPHA_FILE = "../../data/alpha.jld"
const ALPHA_VARIABLE = "alpha"

const DIM_POMDP_STATES = 5
const DIM_STATES = 8
const TERM_STATE_VAR = 1e5
const TERM_STATE = TERM_STATE_VAR * ones(DIM_POMDP_STATES)

const XMIN = -2e3
const XMAX = 2e3
const YMIN = -2e3
const YMAX = 2e3
const PMIN = 0.0
const PMAX = 2 * pi

type Policy
    alpha       :: Matrix{Float64}
    actions     :: Matrix{Float64}
    nactions    :: Int64
    qvals       :: Vector{Float64}

    function Policy(alpha::Matrix{Float64}, actions::Matrix{Float64})
        return new(alpha, actions, size(actions, 2), zeros(size(actions, 2)))
    end # function Policy
end

function get_pomdp()
    return DoubleUAV()
end # function get_pomdp

function gen_pairwise_policy(
        d::DoubleUAV,
        lambda::Float64=1.0,
        filetag::ASCIIString="",
        saveAlpha::Bool=true)
    alpha = qmdp(d, lambda)
    if saveAlpha
        @printf("Writing alpha vector to %s...", ALPHA_FILE)
        tic()
        filename = string("../../data/alpha", filetag, ".jld")
        save(filename, ALPHA_VARIABLE, alpha)
        @printf("done in %.2e.\n\n", toq())
    end # if
    return Policy(alpha, d.pomdp.actions)
end # function gen_policy

function get_pomdp_state(state1::Vector{Float64}, state2::Vector{Float64})
    x1 = state1[X]
    x2 = state2[X]
    y1 = state1[Y]
    y2 = state2[Y]
    p1 = state1[P]
    p2 = state2[P]

    v1 = sqrt(state1[XDOT]^2 + state1[YDOT]^2)
    v2 = sqrt(state2[XDOT]^2 + state2[YDOT]^2)

    dx = x2 - x1
    dy = y2 - y1

    xr = dx * cos(p1) + dy * sin(p1)
    yr = -dx * sin(p1) + dy * cos(p1)
    pr = p2 - p1

    if xr < XMIN || xr > XMAX ||
       yr < YMIN || yr > YMAX ||
       state1[1] == TERM_STATE_VAR ||
       state2[1] == TERM_STATE_VAR
        return TERM_STATE
    else
        return [xr, yr, pr, v1, v2]
    end # if
end # function get_pomdp_state

function get_belief(pstate::Vector{Float64}, grid::RectangleGrid)
    belief = spzeros(NSTATES, 1)
    if pstate[1] == TERM_STATE_VAR
        belief[end] = 1.0
    else
        indices, weights = interpolants(grid, pstate[1:3])
        for i = 1:length(indices)
            belief[indices[i]] = weights[i]
        end # for i
    end # if
    return belief
end # function get_belief

function get_qval!(policy::Policy, belief::SparseMatrixCSC{Float64, Int64})
    fill!(policy.qvals, 0.0)
    for iaction in 1:policy.nactions
        for ib in 1:length(belief.rowval)
            policy.qvals[iaction] += belief.nzval[ib] * policy.alpha[belief.rowval[ib], iaction]
        end # for b
    end # for iaction
end # function get_qval!

function evaluate(policy::Policy, belief::SparseMatrixCSC{Float64,Int64})
    fill!(policy.qvals, 0.0)
    get_qval!(policy, belief)
    ibest = indmax(policy.qvals)
    return policy.actions[:, ibest], ibest
end # function evaluate

function read_policy(actions::Matrix{Float64}, alphafile::ASCIIString)
    alpha = load(alphafile, ALPHA_VARIABLE)
    return Policy(alpha, actions)
end # function read_policy

function read_policy(actions::Matrix{Float64}, alpha::Matrix{Float64})
    return Policy(alpha, actions)
end # function read_policy

function viz_pairwise_policy(d::DoubleUAV, alphafile::ASCIIString=ALPHA_FILE, nbin::Int64=250)
    states = d.pomdp.states
    actions = d.pomdp.actions

    xs = sort(unique(vec(states[1, 1:end - 1])))
    ys = sort(unique(vec(states[2, 1:end - 1])))
    ps = sort(unique(vec(states[3, 1:end - 1])))
    vs = sort(unique(vec(states[4, 1:end - 1])))
    ps_deg = rad2deg(ps)
    grid = RectangleGrid(xs, ys, ps, vs, vs)

    pstart = ps_deg[1]
    pend = ps_deg[end]
    pdiv = ps_deg[2] - ps_deg[1]

    vstart = vs[1]
    vend = vs[end]
    vdiv = vs[2] - vs[1]

    policy = read_policy(actions, alphafile)

    @manipulate for p = pstart:pdiv:pend,
                    v0 = vstart:vdiv:vend,
                    v1 = vstart:vdiv:vend
        # ownship uav
        function get_heat1(x::Float64, y::Float64)
            action, _ = evaluate(policy, get_belief(
                [x, y, deg2rad(p), v0, v1], grid
            ))
            if abs(action[1]) > 0.5
                return 2.0
            end # if
            return rad2deg(action[1])
        end # function get_heat1

        # intruder uav
        function get_heat2(x::Float64, y::Float64)
            action, _ = evaluate(policy, get_belief(
                [x, y, deg2rad(p), v0, v1], grid
            ))
            if abs(action[2]) > 0.5
                return 2.0
            end # if
            return rad2deg(action[2])
        end # function get_heat2

        g = GroupPlot(2, 1, groupStyle = "horizontal sep=3cm")
        push!(g, Axis([
            Plots.Image(get_heat1, (int(XMIN), int(XMAX)),
                        (int(YMIN), int(YMAX)),
                        zmin = -20, zmax = 20,
                        xbins = nbin, ybins = nbin,
                        colormap = ColorMaps.Named("RdBu"), colorbar = false),
            Plots.Node(L">", 0, 0, style="rotate=0,font=\\Huge"),
            Plots.Node(L">", 1800, 1800, style=string("rotate=", p, ",font=\\Huge"))
            ], width="10cm", height="10cm", xlabel="x (m)", ylabel="y (m)", title="Ownship action"))
        push!(g, Axis([
            Plots.Image(get_heat2, (int(XMIN), int(XMAX)),
                        (int(YMIN), int(YMAX)),
                        zmin = -20, zmax = 20,
                        xbins = nbin, ybins = nbin,
                        colormap = ColorMaps.Named("RdBu")),
            Plots.Node(L">", 0, 0, style="rotate=0,font=\\Huge"),
            Plots.Node(L">", 1800, 1800, style=string("rotate=", p, ",font=\\Huge"))],
            width="10cm", height="10cm", xlabel="x (m)", title="Intruder action"))
        g
    end # for p, v0, v1
end # function viz_pairwise_policy

function viz_pairwise_policy(d::DoubleUAV, alpha::Matrix{Float64}, nbin::Int64=250)
    states = d.pomdp.states
    actions = d.pomdp.actions

    xs = sort(unique(vec(states[1, 1:end - 1])))
    ys = sort(unique(vec(states[2, 1:end - 1])))
    ps = sort(unique(vec(states[3, 1:end - 1])))
    vs = sort(unique(vec(states[4, 1:end - 1])))
    ps_deg = rad2deg(ps)
    grid = RectangleGrid(xs, ys, ps, vs, vs)

    pstart = ps_deg[1]
    pend = ps_deg[end]
    pdiv = ps_deg[2] - ps_deg[1]

    vstart = vs[1]
    vend = vs[end]
    vdiv = vs[2] - vs[1]

    policy = read_policy(actions, alpha)

    @manipulate for p = pstart:pdiv:pend,
                    v0 = vstart:vdiv:vend,
                    v1 = vstart:vdiv:vend
        # ownship uav
        function get_heat1(x::Float64, y::Float64)
            action, _ = evaluate(policy, get_belief(
                [x, y, deg2rad(p), v0, v1], grid
            ))
            if abs(action[1]) > 0.5
                return 2.0
            end # if
            return rad2deg(action[1])
        end # function get_heat1

        # intruder uav
        function get_heat2(x::Float64, y::Float64)
            action, _ = evaluate(policy, get_belief(
                [x, y, deg2rad(p), v0, v1], grid
            ))
            if abs(action[2]) > 0.5
                return 2.0
            end # if
            return rad2deg(action[2])
        end # function get_heat2

        g = GroupPlot(2, 1, groupStyle = "horizontal sep=3cm")
        push!(g, Axis([
            Plots.Image(get_heat1, (int(XMIN), int(XMAX)),
                        (int(YMIN), int(YMAX)),
                        zmin = -20, zmax = 20,
                        xbins = nbin, ybins = nbin,
                        colormap = ColorMaps.Named("RdBu"), colorbar = false),
            Plots.Node(L">", 0, 0, style="rotate=0,font=\\Huge"),
            Plots.Node(L">", 1800, 1800, style=string("rotate=", p, ",font=\\Huge"))
            ], width="10cm", height="10cm", xlabel="x (m)", ylabel="y (m)", title="Ownship action"))
        push!(g, Axis([
            Plots.Image(get_heat2, (int(XMIN), int(XMAX)),
                        (int(YMIN), int(YMAX)),
                        zmin = -20, zmax = 20,
                        xbins = nbin, ybins = nbin,
                        colormap = ColorMaps.Named("RdBu")),
            Plots.Node(L">", 0, 0, style="rotate=0,font=\\Huge"),
            Plots.Node(L">", 1800, 1800, style=string("rotate=", p, ",font=\\Huge"))],
            width="10cm", height="10cm", xlabel="x (m)", title="Intruder action"))
        g
    end # for p, v0, v1
end # function viz_pairwise_policy

end # module Pairwise


text:
n = 9 # rand(1:10)

A = matrixdepot("companion", Float32, n)
@test diag(A, -1) == ones(Float32, n-1)

v = [3., 4, 2, 10]
B = matrixdepot("companion", v)
@test B[:,end] == v

println("'companion' passed test...")


text:
################################################################################
## Copyright 2016. Juri Hinz and Jeremy Yee.
## Solutio diagnostics
################################################################################

immutable OptimalBounds

    ## Members
    upper::Array{Float64, 3} # upper bound realizations
    lower::Array{Float64, 3} # lower bound realizations
    uppermean::Array{Float64, 1} # mean upper bound
    lowermean::Array{Float64, 1} # mean lower bound
    upperse::Array{Float64, 1} # standard error of upper bound mean
    lowerse::Array{Float64, 1} # standard error of lower bound mean

    function OptimalBounds(css::CSS, path::Path, mart::Mart,
                           pathPolicy::PathPolicy)

        ## Get information
        tnum = css.tnum
        pnum = css.pnum
        pathnum = path.pathnum
        snum = css.snum
        anum = css.anum
        Scrap = css.Scrap
        Reward = css.Reward
        control = css.control

        ## Containers
        upper = zeros(Array{Float64}(tnum, pnum, pathnum))
        lower = zeros(Array{Float64}(tnum, pnum, pathnum))
        state = Array{Float64}(snum)

        print("Optimality Bounds: ", tnum, ".")
        ## Last time point
        for pos = 1:pnum
            for p = 1:pathnum
                upper[tnum, pos, p] = dot(Scrap(path.sample[tnum, p, :], pos),
                                          path.sample[tnum, p, :])
            end
        end
        lower[tnum, :, :] = upper[tnum, :, :]

        ## Bellman recursion
        t = tnum - 1
        if css.fullcontrol
            while (t > 0)
                print(t, ".")
                for p = 1:pathnum
                    state = path.sample[t, p, :]
                    for pos = 1:pnum
                        ## Lower
                        action = pathPolicy.policy[p, pos, t]
                        next = control[pos, action]
                        lower[t, pos, p] = dot(Reward(t,state, pos, action), state) +
                        mart.mart[t, next, p] + lower[t + 1, next, p]
                        ## Upper
                        next = control[pos, 1]
                        upper[t, pos, p] = dot(Reward(t, state, pos, 1), state) +
                        mart.mart[t, next, p] + upper[t + 1, next, p]
                        ## Find pathwise maximising action
                        for a = 1:anum
                            next = control[pos, a]
                            temp = dot(Reward(t, state, pos, a), state) +
                            mart.mart[t, next, p] + upper[t + 1, next, p]
                            upper[t, pos, p] = maximum([temp; upper[t, pos, p]])
                        end
                    end
                end
                t = t - 1
            end
        else ## not full control of positions
            prob_weight = Array{Float64}(pnum)
            while (t > 0)
                print(t, ".")
                for p = 1:pathnum
                    state = path.sample[t, p, :]
                    for pos = 1:pnum
                        ## Lower
                        action = pathPolicy.policy[p, pos, t]
                        prob_weight = control[pos, action, :]
                        lower[t, pos, p] = dot(Reward(t, state, pos, action), state) +
                        mart.mart[t, pos, action, p] + dot(lower[t + 1, :, p], prob_weight)
                        ## Upper
                        prob_weight = control[pos, 1, :]
                        upper[t, pos, p] = dot(Reward(t, state, pos, 1), state) +
                        mart.mart[t, pos, 1, p] + dot(upper[t + 1, :, p], prob_weight)
                        ## Find pathwise maximising action
                        for a = 1:anum
                            prob_weight = control[pos, a, :]
                            temp = dot(Reward(t, state, pos, a), state) +
                            mart.mart[t, pos, a, p] + dot(upper[t + 1, :, p], prob_weight)
                            upper[t, pos, p] = maximum([temp; upper[t, pos, p]])
                        end
                    end
                end
                t = t - 1
            end
        end

        ## Compute the means and standard errors
        uppermean = mean(upper[1, :, :], 2)[:]
        lowermean = mean(lower[1, :, :], 2)[:]
        upperse = (std(upper[1, :, :], 2) / sqrt(pathnum))[:]
        lowerse = (std(lower[1, :, :], 2) / sqrt(pathnum))[:]

        print("End.\n")

        return new(upper, lower, uppermean, lowermean, upperse, lowerse)

    end

end


text:
@testset "Mango" begin
  println(testname("Mango tests"))
  print("  [  ] Create a json Mango index ")
  result = mango_index(db, ["data", "data2"]; name="myindex", ddoc="mangoddoc")
  @test result["result"] == "created"
  println("\r  [OK] Create a json Mango index")

  print("  [  ] Create a json Mango index with a selector should fail ")
  @test_throws ErrorException mango_index(db, ["data", "data2"]; selector=q"data = bob")
  println("\r  [OK] Create a json Mango index with a selector should fail")

  print("  [  ] Create a json Mango index without fields should fail ")
  @test_throws ErrorException mango_index(db, [])
  println("\r  [OK] Create a json Mango index without fields should fail")

  print("  [  ] Bulk load data ")
  data=[
      Dict("name"=>"adam",    "data"=>"hello",              "data2" => "television"),
      Dict("name"=>"billy",   "data"=>"world",              "data2" => "vocabulary"),
      Dict("name"=>"bob",     "data"=>"world",              "data2" => "organize"),
      Dict("name"=>"cecilia", "data"=>"authenticate",       "data2" => "study"),
      Dict("name"=>"frank",   "data"=>"authenticate",       "data2" => "region"),
      Dict("name"=>"davina",  "data"=>"cloudant",           "data2" => "research"),
      Dict("name"=>"eric",    "data"=>"blobbyblobbyblobby", "data2" => "knowledge")
  ]

  result = createdoc(db, data)
  @test length(result) == length(data)
  println("\r  [OK] Bulk load data")

  print("  [  ] Simple Mango query (equality) ")
  result = mango_query(db, q"data = authenticate"; fields=["name", "data2"])
  @test length(result.docs) == 2
  println("\r  [OK] Simple Mango query (equality)")

  print("  [  ] Simple Mango query (negation) ")
  result = mango_query(db, not(q"data = authenticate"))
  @test length(result.docs) == 5
  println("\r  [OK] Simple Mango query (negation)")

  print("  [  ] Selector from raw json ")
  result = mango_query(db, Selector("{\"data\":{\"\$eq\":\"authenticate\"}}"))
  @test length(result.docs) == 2
  println("\r  [OK] Selector from raw json")

  print("  [  ] Compound Mango query (and) ")
  result = mango_query(db, and([q"data = world", q"data2 = vocabulary"]))
  @test length(result.docs) == 1
  @test result.docs[1]["name"] == "billy"
  println("\r  [OK] Compound Mango query (and)")

  print("  [  ] Compound Mango query (or) ")
  result = mango_query(db, or([q"data = world", q"data2 = region"]))
  @test length(result.docs) == 3
  println("\r  [OK] Compound Mango query (or)")

  print("  [  ] Create a text Mango index ")
  textindex = mango_index(db, [
    Dict("name" => "cust",  "type" => "string"),
    Dict("name" => "value", "type" => "string")
  ])
  @test textindex["result"] == "created"
  println("\r  [OK] Create a text Mango index")

  maxdoc = 102
  createdoc(db, [Dict("cust" => "john", "value" => "hello$x") for x=1:maxdoc])
  print("  [  ] Mango query with multi-page return ")
  result = mango_query(db, q"cust=john")
  count = length(result.docs)
  while length(result.docs) > 0
    result = mango_query(db, q"cust = john", bookmark=result.bookmark)
    count += length(result.docs)
  end
  @test count == maxdoc
  println("\r  [OK] Mango query with multi-page return")

  print("  [  ] Multi-page Mango query as a Task ")
  createdoc(db, [Dict("data" => "paged", "data2" => "world$x") for x=1:maxdoc])
  total = 0
  for page in @task paged_mango_query(db, q"data = paged"; pagesize=10)
    total += length(page.docs)
  end
  @test total == maxdoc
  println("\r  [OK] Multi-page Mango query as a Task ")

  print("  [  ] List indexes ")
  result = listindexes(db)
  @test length(result["indexes"]) == 3
  println("\r  [OK] List indexes")

  print("  [  ] Delete Mango index ")
  result = mango_deleteindex(db; ddoc=textindex["id"], name=textindex["name"], indextype="text")
  @test result["ok"] == true
  println("\r  [OK] Delete Mango index")

  print("  [  ] Delete Mango index with bad index type ")
  @test_throws ErrorException mango_deleteindex(db; ddoc=textindex["id"], name=textindex["name"], indextype="book")
  println("\r  [OK] Delete Mango index with bad index type")

  print("  [  ] Delete Mango index with no name ")
  @test_throws ErrorException mango_deleteindex(db; ddoc=textindex["id"], indextype="json")
  println("\r  [OK] Delete Mango index with no name")
end


text:
using TapirSchedulersBenchmarks
SUITE = TapirSchedulersBenchmarks.setup()


text:
;; ---------------------------------------------------------------------
;;	text_create.jl
;;
;;					Jul/25/2012
;;
;; --------------------------------------------------------------------
(require 'tables)
(load-file "/var/www/data_base/common/rep_common/text_manipulate.jl")
;; --------------------------------------------------------------------
(defun tag_out_proc (port_out tag value)
	(setq out_str_aa (concat 'string "<" tag ">"))
	(setq out_str_bb (concat 'string "</" tag ">\n"))
	(write port_out out_str_aa)
	(write port_out value)
	(write port_out out_str_bb)
	)
;
;; --------------------------------------------------------------------
(defun xml_dict_write_proc_single (fp-file key unit_bb)
	(write fp-file (concat "<" key ">"))
	(tag_out_proc fp-file "name" (table-ref unit_bb "name"))
	(tag_out_proc fp-file "population" (table-ref unit_bb "population"))
	(tag_out_proc fp-file "date_mod" (table-ref unit_bb "date_mod"))
	(write fp-file (concat "</" key ">"))
	)

;; --------------------------------------------------------------------
(defun xml_write_proc (file_out dict_aa)
(let ((fp-file (open-file file_out 'write)))
	(write fp-file "<?xml version=\"1.0\" encoding=\"utf-8\"?>\n")
	(write fp-file "<root>\n")
	(table-walk (lambda (key value)
		 (xml_dict_write_proc_single fp-file key value)) dict_aa)
	(write fp-file "</root>\n")
	)
)

;; --------------------------------------------------------------------
(defun data_prepare_proc ()
	(setq dict_aa (make-table string-hash string=))
;
	(dict_append_proc dict_aa "t2261" "静岡" "87639" "1930-2-5")
	(dict_append_proc dict_aa "t2262" "浜松" "12497" "1930-11-21")
	(dict_append_proc dict_aa "t2263" "沼津" "34593" "1930-3-16")
	(dict_append_proc dict_aa "t2264" "三島" "52481" "1930-9-23")
	(dict_append_proc dict_aa "t2265" "富士" "85192" "1930-8-9")
	(dict_append_proc dict_aa "t2266" "熱海" "82573" "1930-5-19")
	(dict_append_proc dict_aa "t2267" "富士宮" "51948" "1930-9-12")
	(dict_append_proc dict_aa "t2268" "藤枝" "23752" "1930-1-18")
	(dict_append_proc dict_aa "t2269" "御殿場" "51264" "1930-2-11")
	(dict_append_proc dict_aa "t2270" "島田" "37189" "1930-9-24")
)

;; --------------------------------------------------------------------
(format t "*** 開始 ***")
(setq file_out (cadr command-line-args))
(format t file_out)

(setq dict_aa (data_prepare_proc))
(table-walk dict_display_proc_single dict_aa)

(xml_write_proc file_out dict_aa)

(format t "*** 終了 ***")

;; ---------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
# Run this in Escher to visualize a Complete Graph with Vertex identification

using LightGraphs
using NetworkViz
using ThreeJS

main(window) = begin
    num = Signal(10)
    toggle = Signal(false)
    g = CompleteGraph(10)

    push!(window.assets,("ThreeJS","threejs"))
    push!(window.assets,"widgets")
        vbox(
        h1("CompleteGraph Example"),
        vskip(2em),
        drawGraphwithText(g,1)
        ) |> pad(2em)
end


text:
# this is a temporary script for trying out pnp algorithms on images

import ProximalOperators.prox!
import Random.seed!

# run from ProximalAlgorithms.jl folder, e.g.  /home/naomigra/projects/def-mpf/naomigra/2022/pnp_IR/dev/ProximalAlgorithms.jl
include("../src/utilities/denoisers.jl")
include("../src/utilities/image_tools.jl")
include("../src/algorithms/pnp_douglas_rachford.jl")
include("../src/algorithms/admm.jl")

using LinearAlgebra
using ProximalOperators
using MLDatasets
using Flux
using Zygote
using Random

# need a function to run PnP DR on a single example with various gamma, number of iterations,

function run_example_1()
    Random.seed!(1)
    T = Float64 #double precision
    R = real(T)
    k = 0.8 # how much we want to reduce (scale) max rank of A

    example_index = 10
    model_epoch_number = 20
    model_dir = "./src/utilities/saved_models/MNIST"
    encoder_μ, encoder_logvar, decoder = load_model(model_dir,model_epoch_number)

    # backward operator (denoiser)
    # TO DO: look into what kind of noise is used for training
    function vae_denoiser!(x,y)
        # defining the denoiser this way so that it has the same format as prox!
        x .= reconstruct_images(encoder_μ, encoder_logvar, decoder,y)
    end

    test_x, test_y = MNIST.testdata(Float32,example_index)

    m,n = size(test_x)
    A_full_rank = randn(T, (m*n, m*n)) #randn performs better than rand
    max_rank = Int64(floor(m*n*k))
    F = svd(A_full_rank)

    for i in max_rank:(m*n)
        F.S[i] = 0
    end

    A = (F.U * Diagonal(F.S) * F.Vt)/100

    A = diagm(m*n, m*n, ones(T, m*n))

    test_x_vec = vectorize_and_flip(test_x)

    convert_and_save_single_MNIST_image(test_x_vec[:],"./scripts/temp/run_example1","ground_truth")

    #b = A*test_x_vec[:] # add additive noise as well?

    # add noise to x instead of Ax
    b = A*(test_x_vec[:] + randn(T, (length(test_x_vec[:]))))

    #b = b + randn(T,Int64(length(b))) #adding noise
    f = LeastSquares(A, b)

    #lam = R(0.1) * norm(A' * b, Inf)
    #z0 = zeros(R, m*n)
    u0 = ones(m*n)*R(0.5)
    gamma = R(1)

    proxf!(x, u) = prox!(x, f, u, gamma)

    pnp_dr_iter = PnpDrsIteration(J_A! = proxf!, J_B! = vae_denoiser!, u0 = u0)

    # TO DO: find a useful way to look at output
    for (i, pnp_dr_state) in enumerate(Iterators.take(pnp_dr_iter,50))
        println("iteration $i")
        println(norm(pnp_dr_state.res))
        println("\n")
        #convert_and_save_single_MNIST_image(pnp_dr_state.x,"./scripts/temp/run_example1/","iteration_$i")
        resize_and_save_single_MNIST_image(pnp_dr_state.x,"./scripts/temp/run_example1/","iteration_$i")
        #println(norm(pnp_dr_state.res))
        #i+=1
    end

end

function run_example_2()
    T = Float32
    R = real(T)
    k = 0.5 # ratio of rows to columns

    example_index = 11
    model_epoch_number = 20
    model_dir = "./src/utilities/saved_models/MNIST"
    encoder_μ, encoder_logvar, decoder = load_model(model_dir,model_epoch_number)

    # backward operator (denoiser)
    function vae_denoiser!(x,y)
        # defining the denoiser this way so that it has the same format as prox!
        x .= reconstruct_images(encoder_μ, encoder_logvar, decoder,y)
    end

    test_x, test_y = MNIST.testdata(Float32,example_index)
    test_x_vec = vectorize_and_flip(test_x)
    convert_and_save_single_MNIST_image(test_x_vec[:],"./scripts/temp/run_example2","ground_truth")

    m,n = size(test_x)
    A = randn(T, Int64(floor(m*n*k)), (m*n)) # wide matrix
    # b = Ax + noise
    b = A*test_x_vec[:] + randn(T,Int64(floor(m*n*k)) )

    f = LeastSquares(A, b)

    #lam = R(0.1) * norm(A' * b, Inf)
    u0 = zeros(R, m*n)
    #x0 = test_x_vec[:]
    gamma = R(10) / opnorm(A)^2 # constant parameter used douglas rachford examples for scaling the l1 norm regularizer
    gamma = R(1)
    println("gamma is", gamma)

    # forward operator
    proxf!(x,u) = prox!(x,f,u,gamma)

    pnp_dr_iter = PnpDrsIteration(J_A! = proxf!, J_B! = vae_denoiser!, u0 = u0)

    i = 0
    # TO DO: find a useful way to look at output
    for pnp_dr_state in Iterators.take(pnp_dr_iter,10)
        #convert_and_save_single_MNIST_image(clamp.(pnp_dr_state.xhat,0.0,1.0),"./scripts/temp/run_example2/","iteration_$i")
        resize_and_save_single_MNIST_image(pnp_dr_state.x,"./scripts/temp/run_example2/","iteration_$i")
        println(norm(pnp_dr_state.res))
        #println(norm(test_x_vec[:] - pnp_dr_state.xhat))
        #println(norm(pnp_dr_state.xhat))
        i+=1
    end

end

# in this example we use the identity as a denoiser
function run_example_3()
    T = Float32
    R = real(T)
    k = 0.5 # ratio of rows to columns

    example_index = 11
    model_epoch_number = 20
    model_dir = "./src/utilities/saved_models/MNIST"
    encoder_μ, encoder_logvar, decoder = load_model(model_dir,model_epoch_number)

    # backward operator (denoiser)
    #function vae_denoiser!(x,y)
    function identity_denoiser!(x,y)
        # defining the denoiser this way so that it has the same format as prox!
        #x .= reconstruct_images(encoder_μ, encoder_logvar, decoder,y)
        x.=y
    end

    test_x, test_y = MNIST.testdata(Float32,example_index)
    test_x_vec = vectorize_and_flip(test_x)
    convert_and_save_single_MNIST_image(test_x_vec[:],"./scripts/temp/run_example3","ground_truth")

    m,n = size(test_x)
    A = randn(T, Int64(floor(m*n*k)), (m*n)) # wide matrix
    b = A*test_x_vec[:] + randn(T,Int64(floor(m*n*k)) ) # TO DO: values of b are outside (0,1) range, fix this

    # convert_and_save_single_MNIST_image(clamp.(b,0,1),"./scripts/temp/","observed")

    f = LeastSquares(A, b)

    #lam = R(0.1) * norm(A' * b, Inf)
    u0 = zeros(R, m*n)
    #gamma = R(10) / opnorm(A)^2 # constant parameter for douglas rachford using g is the l1 norm
    gamma = 1.0f0
    println("\n gamma is $gamma \n")

    # forward operator
    proxf!(x,u) = prox!(x,f,u,gamma)

    pnp_dr_iter = PnpDrsIteration(J_A! = proxf!, J_B! = identity_denoiser!, u0 = u0)

    i = 0
    for pnp_dr_state in Iterators.take(pnp_dr_iter,50)
        resize_and_save_single_MNIST_image(pnp_dr_state.x,"./scripts/temp/run_example3/","iteration_$i")
        i+=1
    end

end

# temporary function for trying the optimized decoder denoiser
function run_denoising_example()
    Random.seed!(1)
    example_index = 11
    test_x, test_y = MNIST.testdata(Float32,example_index)

    test_x_vec = vectorize_and_flip(test_x)[:]
    resize_and_save_single_MNIST_image(test_x_vec,"./scripts/temp/run_denoising","clean_img")

    model_epoch_number = 20
    model_dir = "./src/utilities/saved_models/MNIST"
    encoder_μ, encoder_logvar, decoder = load_model(model_dir,model_epoch_number)

    z0 = add_additive_noise(ones(20)) #initial point for gradient descent

    #z0 = zeros(20)
    #alternate starting point
    #μ = encoder_μ(test_x_vec)
    #logvar = encoder_logvar(test_x_vec)
    #z0 = μ + randn(Float32, size(logvar)) .* exp.(0.5f0 * logvar)

    y_noisy = add_additive_noise(test_x_vec,0.2) # vectorized noisy image

    #shift and rescale to have values between 0 and 1
    y_temp = (y_noisy .- min(0, minimum(y_noisy)))
    y_noisy .= y_temp / maximum(y_temp)

    #y_noisy = test_x_vec #trying without noise

    resize_and_save_single_MNIST_image(y_noisy,"./scripts/temp/run_denoising","before_denoising")

    loss_function = decoder_loss_function

    num_iter = 1000 #number of iterations of gradient descent
    z_out = gradient_descent(loss_function,z0,decoder,y_noisy,num_iter,0.05)

    x_out = decoder(z_out)
    resize_and_save_single_MNIST_image(x_out,"./scripts/temp/run_denoising","after_denoising")

end

#examples 4 and 5 are for studying convergence
function run_example_4()
    T = Float32
    A = T[
        1.0 -2.0 3.0 -4.0 5.0
        2.0 -1.0 0.0 -1.0 3.0
        -1.0 0.0 4.0 -3.0 2.0
        -1.0 -1.0 -1.0 1.0 3.0
    ]
    b = T[1.0, 2.0, 3.0, 4.0]

    m, n = size(A)

    R = real(T)

    lam = R(0.1) * norm(A' * b, Inf)

    f = LeastSquares(A, b)
    g = NormL1(lam)

    u0 = zeros(R, n) # initial point of zero

    gamma=R(10) / opnorm(A)^2
    gamma = R(1) # gamma doesn't currently do anything in the algorithm!

    proxf!(x,z) = prox!(x,f,z,gamma)
    denoiser!(y,r) = prox!(y,g,r,gamma)

    pnp_dr_iter = PnpDrsIteration(J_A! = proxf!, J_B! = denoiser!, u0 = u0)

    for pnp_dr_state in Iterators.take(pnp_dr_iter, 30)
        println(norm(pnp_dr_state.res))
    end

end

struct NegativeConjugate{T}
    c::Conjugate{T}
end

function prox!(y, nc::NegativeConjugate, x, gamma)
    prox!(y,nc.c,-x,gamma)
end

# trying ADMM in place of Douglas Rachford
function run_example_5()

    T = Float32
    A = T[
        1.0 -2.0 3.0 -4.0 5.0
        2.0 -1.0 0.0 -1.0 3.0
        -1.0 0.0 4.0 -3.0 2.0
        -1.0 -1.0 -1.0 1.0 3.0
    ]
    b = T[1.0, 2.0, 3.0, 4.0]
    m, n = size(A)
    R = real(T)

    #lam = R(0.1) * norm(A' * b, Inf)
    lam = R(1)

    f = LeastSquares(A, b)
    g = NormL1(lam)

    fdual = NegativeConjugate(Conjugate(f))
    gdual = Conjugate(g)

    u0 = zeros(R, n)

    #gamma=R(10) / opnorm(A)^2
    #gamma= R(10)
    gamma = R(1)

    proxf!(x,u) = prox!(x,f,u,gamma)
    #proxf!(x,u) = prox!(x,fdual,u,gamma) #apply DRS to dual problem

    denoiser!(x,u) = prox!(x,g,u,gamma)
    #denoiser!(x,u) = prox!(x,gdual,u,gamma) #apply DRS to dual problem

    #dr_iter = DouglasRachfordIteration(f=f, g=g, x0=x0, gamma=gamma)
    admm_iter = AdmmIteration(f = f, g = g, x0 = u0, gamma = gamma)
    pnp_dr_iter = PnpDrsIteration(J_A! = proxf!, J_B! = denoiser!, u0 = u0)

    for (admm_state, pnp_dr_state) in Iterators.take(zip(admm_iter, pnp_dr_iter), 30)
        #println(norm(admm_state.u)) # converges to 2.072
        println(norm(admm_state.x - admm_state.y)) #converges to zero
        println(norm(pnp_dr_state.res)) #converges to zero
        println("\n")
    end

end

function run_example_6()
    Random.seed!(2)
    T = Float64 #single precision
    R = real(T)

    # take an example image from the test data set
    example_index = 10
    test_x, test_y = MNIST.testdata(Float32,example_index)

    # vectorize image and save ground truth for comparison
    test_x_vec = vectorize_and_flip(test_x)[:]
    resize_and_save_single_MNIST_image(test_x_vec,"./scripts/temp/run_example6","ground_truth")

    # load the VAE model
    model_epoch_number = 20
    model_dir = "./src/utilities/saved_models/MNIST"
    encoder_μ, encoder_logvar, decoder = load_model(model_dir,model_epoch_number)

    test_x, test_y = MNIST.testdata(Float32,example_index)

    # generate a random low rank matrix
    k = 0.8 # how much we want to reduce (scale) max rank of A
    m,n = size(test_x)
    A_full_rank = randn(T, (m*n, m*n)) #randn performs better than rand
    max_rank = Int64(floor(m*n*k))
    F = svd(A_full_rank)

    for i in max_rank:(m*n)
        F.S[i] = 0
    end
    A = (F.U * Diagonal(F.S) * F.Vt)
    A = diagm(m*n, m*n, ones(T, m*n))

    # add noise to x instead of Ax
    #b = A*(test_x_vec[:] + randn(T, (length(test_x_vec[:]))))
    b = A*(test_x_vec[:])
    b = b + randn(T,Int64(length(b))) #adding noise

    f = LeastSquares(A, b)
    #z0 = zeros(R, m*n)
    u0 = ones(m*n)*R(0.5)
    gamma = R(1)

    #function proxf!(x,z)
    #    prox!(x, f, z, gamma)
    #    x1 = (x .- min(0, minimum(x)))
    #    x .= x1 / maximum(x1)
    #end

    proxf!(x, u) = prox!(x, f, u, gamma)

    num_iter = 20 # number of iterations of GD in the denoiser
    denoiser!(x, u) = decoder_denoiser!(x,u,encoder_μ, encoder_logvar, decoder, num_iter)

    pnp_dr_iter = PnpDrsIteration(J_A! = proxf!, J_B! = denoiser!, u0 = u0)

    for (i, pnp_dr_state) in enumerate(Iterators.take(pnp_dr_iter,30))
        println("at PnP iteration $i")
        println("the residual value is $(norm(pnp_dr_state.res)) \n")
        #println("an intermediate value is $(norm(pnp_dr_state.r)) \n")
        #println("norm of y iterate is $(norm(pnp_dr_state.y)) \n")
        println("Saving output...")
        #convert_and_save_single_MNIST_image(pnp_dr_state.x,"./scripts/temp/run_example6/","iteration_$i")
        resize_and_save_single_MNIST_image(pnp_dr_state.x,"./scripts/temp/run_example6/","iteration_$i")
        #println(norm(pnp_dr_state.res))
        #i+=1
    end

end


text:
include("../../src/ToyModelsMD.jl")
using .ToyModelsMD

function cbfunc(xstore::Array,system::AbstractSystem,args...)
  # calculate the stuff
  push!(xstore,copy(system.x))
end

### define system ###
# mueller-brown potential
mb = MullerBrown()
# gaussian force model
Mx = 21
My = 21
Mt = 21
shape = (Mx*My*Mt,2)
mux = rand(Float64,21)
muy = rand(Float64,21)
mut = rand(Float64,21)
sigx = rand(Float64,21)
sigy = rand(Float64,21)
sigt = rand(Float64,21)
cs = zeros(Float64,shape)
gm = GaussianModel2D(Mx,My,Mt,mux,muy,mut,sigx,sigy,sigt,cs)
# mixed model
model = MixedModel(2,[mb,gm])
# thermostat and system
langevin = Langevin(1.0,1.0)
system = ThermostattedSystem(2,0.0,[0.63,0.03],zeros(Float64,2),langevin)

# define callback function
xstore = []
cbf(system,args...) = cbfunc(xstore,system,args)
cb = Callback(1,cbf)

dt = 0.0001
tf = 0.0015
nsteps = convert(Int64, round(tf/dt))

system.x = [0.63, 0.03]
integrator = StochasticEuler(system,model)
runtraj!(nsteps,dt,integrator,cb)
println(xstore)


text:
module ONNX
  const _ProtoBuf_Top_ = @static isdefined(parentmodule(@__MODULE__), :_ProtoBuf_Top_) ? (parentmodule(@__MODULE__))._ProtoBuf_Top_ : parentmodule(@__MODULE__)
  include("onnx_pb.jl")
  include("read.jl")
  include("write.jl")
  include("conversions.jl")
  include("show.jl")
  include("ops.jl")
  include("load.jl")
  include("save.jl")
end


text:
##################################################################
"""
$(TYPEDEF)

Extendable sparse matrix. A nonzero  entry of this matrix is contained
either in cscmatrix, or in lnkmatrix, never in both.

$(TYPEDFIELDS)
"""
mutable struct ExtendableSparseMatrix{Tv,Ti<:Integer} <: AbstractSparseMatrix{Tv,Ti}
    """
    Final matrix data
    """
    cscmatrix::SparseMatrixCSC{Tv,Ti}

    """
    Linked list structure holding data of extension
    """
    lnkmatrix::Union{SparseMatrixLNK{Tv,Ti},Nothing}

    """
    Time stamp of last pattern update
    """
    pattern_timestamp::Float64
end

"""
$(SIGNATURES)

Create empty ExtendableSparseMatrix.
"""
function ExtendableSparseMatrix{Tv,Ti}(m::Integer, n::Integer) where{Tv,Ti<:Integer}
    ExtendableSparseMatrix{Tv,Ti}(spzeros(Tv,Ti,m,n),nothing,time())
end

"""
$(SIGNATURES)

Create empty ExtendableSparseMatrix.
"""
function ExtendableSparseMatrix(valuetype::Type{Tv},indextype::Type{Ti},m::Integer, n::Integer) where{Tv,Ti<:Integer}
    ExtendableSparseMatrix{Tv,Ti}(m,n)
end

"""
$(SIGNATURES)

Create empty ExtendablSparseMatrix.
This is a pendant to spzeros.
"""
ExtendableSparseMatrix(valuetype::Type{Tv},m::Integer, n::Integer) where{Tv}=ExtendableSparseMatrix{Tv,Int}(m,n)

"""
$(SIGNATURES)

Create empty ExtendableSparseMatrix.
This is a pendant to spzeros.
"""
ExtendableSparseMatrix(m::Integer, n::Integer)=ExtendableSparseMatrix{Float64,Int}(m,n)

"""
$(SIGNATURES)

  Create ExtendableSparseMatrix from SparseMatrixCSC
"""
function ExtendableSparseMatrix(csc::SparseMatrixCSC{Tv,Ti}) where{Tv,Ti<:Integer}
    return ExtendableSparseMatrix{Tv,Ti}(csc, nothing, time())
end

"""
$(SIGNATURES)

Update element of the matrix  with operation `op`.
This can replace the following code and save one index
search during acces:

```@example
using ExtendableSparse # hide
A=ExtendableSparseMatrix(3,3)
A[1,2]+=0.1
A
```

```@example
using ExtendableSparse # hide

A=ExtendableSparseMatrix(3,3)
updateindex!(A,+,0.1,1,2)
A
```
"""
function updateindex!(ext::ExtendableSparseMatrix{Tv,Ti}, op,v, i,j) where{Tv,Ti<:Integer}
    k=findindex(ext.cscmatrix,i,j)
    if k>0
        ext.cscmatrix.nzval[k]=op(ext.cscmatrix.nzval[k],v)
    else
        if ext.lnkmatrix==nothing
            ext.lnkmatrix=SparseMatrixLNK{Tv, Ti}(ext.cscmatrix.m, ext.cscmatrix.n)
        end
        updateindex!(ext.lnkmatrix,op,v,i,j)
    end
    ext
end

"""
$(SIGNATURES)

Find index in CSC matrix and set value if it exists. Otherwise,
set index in extension.
"""

function Base.setindex!(ext::ExtendableSparseMatrix{Tv,Ti}, v, i,j) where{Tv,Ti<:Integer}
    k=findindex(ext.cscmatrix,i,j)
    if k>0
        ext.cscmatrix.nzval[k]=v
    else
        if ext.lnkmatrix==nothing
            ext.lnkmatrix=SparseMatrixLNK{Tv, Ti}(ext.cscmatrix.m, ext.cscmatrix.n)
        end
        ext.lnkmatrix[i,j]=v
    end
end

"""
$(SIGNATURES)

Find index in CSC matrix and return value, if it exists.
Otherwise, return value from extension.
"""
function Base.getindex(ext::ExtendableSparseMatrix{Tv,Ti},i::Integer, j::Integer) where{Tv,Ti<:Integer}
    k=findindex(ext.cscmatrix,i,j)
    if k>0
        return ext.cscmatrix.nzval[k]
    elseif ext.lnkmatrix==nothing
        return zero(Tv)
    else
        return ext.lnkmatrix[i,j]
    end
end

"""
$(SIGNATURES)

Size of ExtendableSparseMatrix.
"""
Base.size(ext::ExtendableSparseMatrix) = (ext.cscmatrix.m, ext.cscmatrix.n)

"""
$(SIGNATURES)

If there are new entries in extension, create new CSC matrix by adding the
cscmatrix and linked list matrix and reset the linked list based extension.
"""
function flush!(ext::ExtendableSparseMatrix{Tv,Ti}) where {Tv, Ti<:Integer}
    if ext.lnkmatrix!=nothing && nnz(ext.lnkmatrix)>0
        ext.cscmatrix=ext.lnkmatrix+ext.cscmatrix
        ext.lnkmatrix=nothing
        ext.pattern_timestamp
    end
    return ext
end

"""
$(SIGNATURES)

[`flush!`](@ref) and return number of nonzeros in ext.cscmatrix.
"""
function SparseArrays.nnz(ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    return nnz(ext.cscmatrix)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and return nonzeros in ext.cscmatrix.
"""
function SparseArrays.nonzeros(ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    return nonzeros(ext.cscmatrix)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and return rowvals in ext.cscmatrix.
"""
function SparseArrays.rowvals(ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    return rowvals(ext.cscmatrix)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and return colptr of  in ext.cscmatrix.
"""
function colptrs(ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    return ext.cscmatrix.colptr
end

"""
$(SIGNATURES)

[`flush!`](@ref) and return findnz(ext.cscmatrix).
"""
function SparseArrays.findnz(ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    return findnz(ext.cscmatrix)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and return LU factorization of ext.cscmatrix
"""
function LinearAlgebra.lu(ext::ExtendableSparseMatrix)
    @inbounds flush!(ext)
    return LinearAlgebra.lu(ext.cscmatrix)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and multiply with ext.cscmatrix
"""
function  LinearAlgebra.mul!(r::AbstractArray{T,1} where T, ext::ExtendableSparse.ExtendableSparseMatrix, x::AbstractArray{T,1} where T)
    @inbounds flush!(ext)
    return LinearAlgebra.mul!(r,ext.cscmatrix,x)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and ldiv with ext.cscmatrix
"""
function  LinearAlgebra.ldiv!(r::AbstractArray{T,1} where T, ext::ExtendableSparse.ExtendableSparseMatrix, x::AbstractArray{T,1} where T)
    @inbounds flush!(ext)
    return LinearAlgebra.ldiv!(r,ext.cscmatrix,x)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and multiply with ext.cscmatrix
"""
function  LinearAlgebra.mul!(r::AbstractArray{T,2} where T, ext::ExtendableSparse.ExtendableSparseMatrix, x::AbstractArray{T,2} where T)
    @inbounds flush!(ext)
    return LinearAlgebra.mul!(r,ext.cscmatrix,x)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and ldiv with ext.cscmatrix
"""
function  LinearAlgebra.ldiv!(r::AbstractArray{T,2} where T, ext::ExtendableSparse.ExtendableSparseMatrix, x::AbstractArray{T,2} where T)
    @inbounds flush!(ext)
    return LinearAlgebra.ldiv!(r,ext.cscmatrix,x)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and calculate norm from cscmatrix
"""
function LinearAlgebra.norm(A::ExtendableSparseMatrix, p::Real=2)
    @time @inbounds flush!(A)
    return LinearAlgebra.norm(A.cscmatrix,p)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and calculate opnorm from cscmatrix
"""
function LinearAlgebra.opnorm(A::ExtendableSparseMatrix, p::Real=2)
    @inbounds flush!(A)
    return LinearAlgebra.opnorm(A.cscmatrix,p)
end

"""
$(SIGNATURES)

[`flush!`](@ref) and calculate cond from cscmatrix
"""
function LinearAlgebra.cond(A::ExtendableSparseMatrix, p::Real=2)
    @inbounds flush!(A)
    return LinearAlgebra.cond(A.cscmatrix,p)
end

"""
$(SIGNATURES)

Add SparseMatrixCSC matrix and [`ExtendableSparseMatrix`](@ref)  ext.
"""
function Base.:+(ext::ExtendableSparseMatrix{Tv,Ti},csc::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti<:Integer}
    @inbounds flush!(ext)
    return ext.cscmatrix+csc
end

"""
$(SIGNATURES)

Subtract  SparseMatrixCSC matrix from  [`ExtendableSparseMatrix`](@ref)  ext.
"""
function Base.:-(ext::ExtendableSparseMatrix{Tv,Ti},csc::SparseMatrixCSC{Tv,Ti}) where {Tv,Ti<:Integer}
    @inbounds flush!(ext)
    return ext.cscmatrix-csc
end

"""
$(SIGNATURES)

Subtract  [`ExtendableSparseMatrix`](@ref)  ext from  SparseMatrixCSC.
"""
function Base.:-(csc::SparseMatrixCSC{Tv,Ti},ext::ExtendableSparseMatrix{Tv,Ti}) where {Tv,Ti<:Integer}
    @inbounds flush!(ext)
    return csc - ext.cscmatrix
end


text:
# Copyright Andy Nowacki 2015-, all rights reserved.
# See the file LICENSE.md for licence details.

"""
# CIJ

Deal with linear elastic constants, especially for geophysical applications.

##
"""
module CIJ

using Dates: now
import LinearAlgebra
using LinearAlgebra: cross, dot, norm
using Printf: @printf

using StaticArrays
using Rotations

export
    # Types and constructors
    EC,
    # Constructors of elastic constants
    Panning_VTI,
    global_VTI,
    grechka_cracks,
    grechka_cracks!,
    iso,
    pitl,
    tandon_and_weng,
    thom,
    thom_st,
    # Properties of a tensor
    Au,
    C2S,
    C2S!,
    Hill_average,
    HillG,
    HillK,
    Reuss_average,
    ReussG,
    ReussK,
    VRH,
    Voigt_average,
    VoigtG,
    VoigtK,
    is_iso,
    is_stable,
    phase_vels,
    # Conversion between representations
    S2C,
    S2C!,
    cij,
    cijkl,
    # Operations
    rot3,
    symm,
    symm!

"Lookup index matrix to convert between elasticity tensor and Voigt matrix"
const VOIGT_CONTRACTION_MATRIX = @SMatrix [1 6 5
                                           6 2 4
                                           5 4 3]

"Default precision of floating point numbers used in `EC` type"
const DEFAULT_FLOAT = Float64

# Types
include("types.jl")

# Functions
include("initialisation.jl")
include("inquiry.jl")
include("conversion.jl")
include("transformation.jl")
include("util.jl")
include("averaging.jl")
include("effective_medium.jl")
include("measures.jl")
include("velocities.jl")
include("io.jl")

# Data
include("sample_data.jl")

end # module CIJ