text:
"""
    TimeSeries{D, Ω, I, T, TT} <: AbstractDiagnostic

A diagnostic for collecting and storing time series.
"""
struct TimeSeries{D, Ω, I, T, TT} <: AbstractDiagnostic
    diagnostic :: D
     frequency :: Ω
      interval :: I
          data :: T
          time :: Vector{TT}
end

"""
    TimeSeries(diagnostic, model; frequency=nothing, interval=nothing)

A `TimeSeries` `Diagnostic` that records a time series of `diagnostic(model)`.

Example
=======
```julia
julia> model = IncompressibleModel(grid=RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1)));

julia> max_u = TimeSeries(FieldMaximum(abs, model.velocities.u), model; frequency=1)

julia> model.diagnostics[:max_u] = max_u; data(model.velocities.u) .= π; time_step!(model, Nt=3, Δt=1e-16)

julia> max_u.data
3-element Array{Float64,1}:
 3.141592653589793
 3.1415926025389127
 3.1415925323439517
```
"""
function TimeSeries(diagnostic, model; frequency=nothing, interval=nothing)
    TD = typeof(diagnostic(model))
    TT = typeof(model.clock.time)
    return TimeSeries(diagnostic, frequency, interval, TD[], TT[])
end

@inline (time_series::TimeSeries)(model) = time_series.diagnostic(model)

function run_diagnostic(model, diag::TimeSeries)
    push!(diag.data, diag(model))
    push!(diag.time, model.clock.time)
    return nothing
end

"""
    TimeSeries(diagnostics::NamedTuple, model; frequency=nothing, interval=nothing)

A `TimeSeries` `Diagnostic` that records a `NamedTuple` of time series of
`diag(model)` for each `diag` in `diagnostics`.

Example
=======
```julia
julia> model = IncompressibleModel(grid=RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1))); Δt = 1.0;

julia> cfl = TimeSeries((adv=AdvectiveCFL(Δt), diff=DiffusiveCFL(Δt)), model; frequency=1);

julia> model.diagnostics[:cfl] = cfl; time_step!(model, Nt=3, Δt=Δt)

julia> cfl.data
(adv = [0.0, 0.0, 0.0, 0.0], diff = [0.0002688, 0.0002688, 0.0002688, 0.0002688])

julia> cfl.adv
4-element Array{Float64,1}:
 0.0
 0.0
 0.0
 0.0
```
"""
function TimeSeries(diagnostics::NamedTuple, model; frequency=nothing, interval=nothing)
    TT = typeof(model.clock.time)
    TDs = Tuple(typeof(diag(model)) for diag in diagnostics)
    data = NamedTuple{propertynames(diagnostics)}(Tuple(T[] for T in TDs))
    return TimeSeries(diagnostics, frequency, interval, data, TT[])
end

function (time_series::TimeSeries{<:NamedTuple})(model)
    names = propertynames(time_series.diagnostic)
    return NamedTuple{names}(Tuple(diag(model) for diag in time_series.diagnostics))
end

function run_diagnostic(model, diag::TimeSeries{<:NamedTuple})
    ntuple(Val(length(diag.diagnostic))) do i
        Base.@_inline_meta
        push!(diag.data[i], diag.diagnostic[i](model))
    end
    push!(diag.time, model.clock.time)
    return nothing
end

Base.getproperty(ts::TimeSeries{<:NamedTuple}, name::Symbol) = get_time_series_field(ts, name)

"Returns a field of time series, or a field of `time_series.data` when possible."
function get_time_series_field(time_series, name)
    if name ∈ propertynames(time_series)
        return getfield(time_series, name)
    else
        return getproperty(time_series.data, name)
    end
end


text:
"""
    rankUpdate!(C, A)
    rankUpdate!(C, A, α)
    rankUpdate!(C, A, α, β)

A rank-k update, C := β*C + α*A'A, of a Hermitian (Symmetric) matrix.

`α` and `β` both default to 1.0.  When `α` is -1.0 this is a downdate operation.
The name `rankUpdate!` is borrowed from [https://github.com/andreasnoack/LinearAlgebra.jl]
"""
function rankUpdate! end

function rankUpdate!(C::AbstractMatrix, a::AbstractArray, α, β)
    error("We haven't implemented a method for $(typeof(C)), $(typeof(a)). Please file an issue on GitHub.")
end

function rankUpdate!(C::HermOrSym{T,S}, a::StridedVector{T}, α, β) where {T,S}
    isone(β) || throw(ArgumentError("isone(β) is false"))
    BLAS.syr!(C.uplo, T(α), a, C.data)
    C  ## to ensure that the return value is HermOrSym
end

function rankUpdate!(C::HermOrSym{T,S}, A::StridedMatrix{T}, α, β) where {T,S}
    BLAS.syrk!(C.uplo, 'N', T(α), A, T(β), C.data)
    C
end

function rankUpdate!(C::HermOrSym{T,S}, A::SparseMatrixCSC{T}, α, β) where {T,S}
    A.m == size(C, 2) || throw(DimensionMismatch())
    C.uplo == 'L' || throw(ArgumentError("C.uplo must be 'L'"))
    Cd, rv, nz = C.data, rowvals(A), nonzeros(A)
    isone(β) || rmul!(LowerTriangular(Cd), β)
    @inbounds for jj in 1:A.n
        rangejj = nzrange(A, jj)
        lenrngjj = length(rangejj)
        for (k, j) in enumerate(rangejj)
            anzj = α * nz[j]
            colj = view(Cd, :, rv[j])
            for i in k:lenrngjj
                kk = rangejj[i]
                colj[rv[kk]] += nz[kk] * anzj
            end
        end
    end
    C
end

function rankUpdate!(C::HermOrSym, A::BlockedSparse, α, β)
    rankUpdate!(C, sparse(A), α, β)
end

function rankUpdate!(C::HermOrSym{T,Diagonal{T,Vector{T}}}, A::StridedMatrix{T}, α, β) where {T,S}
    Cdiag = C.data.diag
    @. Cdiag = β * Cdiag

    for i in 1:length(Cdiag)
        Arow = view(A, i, :)
        Cdiag[i] = Cdiag[i] + α * Arow'Arow
    end

    C
end

function rankUpdate!(C::HermOrSym{T,UniformBlockDiagonal{T}}, A::StridedMatrix{T}, α, β) where {T,S}
    Cdat = C.data.data
    isone(β) || (Cdat .*= β)
    blksize = size(Cdat, 1)

    for k in axes(Cdat, 3)
        ioffset = (k - 1) * blksize
        joffset = (k - 1) * blksize
        for i in 1:blksize, j in 1:i
            iind = ioffset + i
            jind = joffset + j
            AtAij = 0
            for idx in axes(A, 2)
                # because the second is actually A', we swap index orders
                AtAij += A[iind,idx] * A[jind,idx]
            end
            Cdat[i,j,k] += α * AtAij
        end
    end

    C
end

function rankUpdate!(
    C::HermOrSym{T,Diagonal{T,Vector{T}}},
    A::SparseMatrixCSC{T},
    α,
    β,
) where {T}
    dd = C.data.diag
    A.m == length(dd) || throw(DimensionMismatch())
    isone(β) || rmul!(dd, β)
    all(isone.(diff(A.colptr))) || throw(ArgumentError("Columns of A must have exactly 1 nonzero"))
    for (r, nz) in zip(rowvals(A), nonzeros(A))
        dd[r] += α * abs2(nz)
    end
    C
end

function rankUpdate!(C::HermOrSym{T,Diagonal{T}}, A::BlockedSparse{T}, α, β) where {T}
    rankUpdate!(C, sparse(A), α, β)
end

function rankUpdate!(
    C::HermOrSym{T,UniformBlockDiagonal{T}},
    A::BlockedSparse{T,S},
    α,
    β,
) where {T,S}
    Ac = A.cscmat
    cp = Ac.colptr
    all(diff(cp) .== S) || throw(ArgumentError("Columns of A must have exactly $S nonzeros"))
    Cdat = C.data.data
    j, k, l = size(Cdat)
    S == j == k && div(Ac.m, S) == l ||
    throw(DimensionMismatch("div(A.cscmat.m, S) ≠ size(C.data.data, 3)"))
    nz = Ac.nzval
    rv = Ac.rowval
    for j = 1:Ac.n
        nzr = nzrange(Ac, j)
        BLAS.syr!('L', α, view(nz, nzr), view(Cdat, :, :, div(rv[last(nzr)], S)))
    end
    C
end
#=  I don't think Diagonal A can occur after the terms with the same grouping factor have been amalgamated.
function rankUpdate!(C::HermOrSym{T,Diagonal{T}}, A::Diagonal{T}, α, β) where {T}
    Cdiag = C.data.diag
    if length(Cdiag) ≠ length(A.diag)
        throw(DimensionMismatch("length(C.data.diag) ≠ length(A.diag)"))
    end

    Cdiag .= β .* Cdiag .+ α .* abs2.(A.diag)
    C
end

function rankUpdate!(C::HermOrSym{T,Matrix{T}}, A::Diagonal{T}, α, β) where {T}
    Adiag, Cdata = A.diag, C.data
    length(Adiag) == size(C, 2) || throw(DimensionMismatch())
    for (i, a) in zip(diagind(Cdata), Adiag)
        Cdata[i] = β * Cdata[i] + α * abs2(a)
    end
    C
end
=#

ext: jl


text:
### Could be part of PeaceVote
module Recruiters

using DemeNet: Certificate, Signer, AbstractID, Deme, Notary, Cypher, DemeSpec, Profile, AbstractID, ID, DHsym, DHasym, Intent, save
using DiffieHellman: diffiehellman
using Sockets
import DemeNet: serialize, deserialize
import Base.Dict

using SMTPClient

using Pkg.TOML

include("debug.jl")

function stack(io::IO,msg::Vector{UInt8})
    frontbytes = reinterpret(UInt8,Int16[length(msg)])
    item = UInt8[frontbytes...,msg...]
    write(io,item)
end

function unstack(io::IO)
    sizebytes = [read(io,UInt8),read(io,UInt8)]
    size = reinterpret(Int16,sizebytes)[1]

    msg = UInt8[]
    for i in 1:size
        push!(msg,read(io,UInt8))
    end
    return msg
end

struct Tooken
    tooken::Int
end

struct TookenID #{T<:AbstractID} <: AbstractID
    cert::Certificate{Profile}
    tooken::Tooken
end

Dict(id::TookenID) = Dict("cert"=>Dict(id.cert),"tooken"=>id.tooken.tooken)

function TookenID(dict::Dict)
    id = Certificate{Profile}(dict["cert"])
    tooken = Tooken(dict["tooken"])
    return TookenID(id,tooken)
end

function serialize(io::IO,x::Tooken)
    bytes = reinterpret(UInt8,Int[x.tooken])
    write(io,bytes)
end

function deserialize(io::IO,::Type{Tooken})
    bytes = UInt8[]
    for i in 1:8
        byte = read(io,UInt8)
        push!(bytes,byte)
    end
    return Tooken(reinterpret(Int,bytes)[1])
end

struct CertifierConfig{T<:Any}
    tookenca::ID ### authorithies who can issue tookens. Server allows to add new tookens only from them.
    serverid::ID ### Server receiveing tookens and the member identities. Is also the one which signs and issues the certificates.
    tookenport::T
    certifierport::T
end

function validate!(tookens::Set,tooken)
    if tooken in tookens
        pop!(tookens,tooken)
        return true
    else
        return false
    end
end

### How to send back!

struct SecureRegistrator{T}
    server
    daemon
    messages::Channel{T}
end

function SecureRegistrator{T}(port,deme::Deme,validate::Function,signer::Signer) where T<:Any

    server = listen(port)
    messages = Channel{T}()

    dh = DHsym(deme,signer)

    daemon = @async while true
        socket = accept(server)
        @async begin

            key, id = diffiehellman(socket,dh)

            @assert validate(id)

            securesocket = deme.cypher.secureio(socket,key)
            message = deserialize(securesocket,T)
            put!(messages,message)
        end
    end

    SecureRegistrator(server,daemon,messages)
end

struct TookenCertifier
    server
    daemon
    tookens::Set{Tooken}
    tickets::Channel{Tuple{Tooken,Certificate{Profile},Certificate{ID}}}
end

function TookenCertifier(port,deme::Deme,signer::Signer; cert=nothing, sendcopy=false)
    tookens = Set{Tooken}()
    tickets = Channel{Tuple{Tooken,Certificate{Profile},Certificate{ID}}}(Inf)

    server = listen(port)
    dh = DHasym(deme,signer)

    daemon = @async while true
        socket = accept(server)
        @async begin
            cert==nothing || serialize(socket,cert)

            key, id = diffiehellman(socket,dh)
            securesocket = deme.cypher.secureio(socket,key)

            id = deserialize(securesocket,TookenID) ### Need to implement

            @assert id.tooken in tookens
            pop!(tookens,id.tooken)

            intentid = Intent(id.cert,deme.notary)

            regid = intentid.reference

            # here one could verify that the id is indeed real
            # validate(profile,id.tooken)

            idcert = Certificate(regid,signer)

            #tickets[id.tooken] = cert

            push!(tickets,(id.tooken,id.cert,idcert))

            # This part could be optional
            sendcopy && serialize(securesocket,cert) ### For this one we already jnow
        end
    end

    TookenCertifier(server,daemon,tookens,tickets)
end

struct Certifier
    tookenrecorder::SecureRegistrator{Tooken}
    tookencertifier::TookenCertifier
    daemon
end

function Certifier(config::CertifierConfig,deme::Deme,signer::Signer; cert=nothing, sendcopy=false)

    # It could be added as part of config
    # intent = Intent(config.serverid,deme.notary)
    # @assert intent.reference==deme.spec.maintainer "The certificate invalid for this deme"
    # @assert intent.document==signer.id "The provided signer does not match the certificate"

    tookenrecorder = SecureRegistrator{Tooken}(config.tookenport,deme,x->x in config.tookenca,signer)
    tookencertifier = TookenCertifier(config.certifierport,deme,signer; cert=cert, sendcopy=sendcopy)

    daemon = @async while true
        tooken = take!(tookenrecorder.messages)
        push!(tookencertifier.tookens,tooken)
    end

    return Certifier(tookenrecorder,tookencertifier,daemon)
end

function addtooken(cc::CertifierConfig,deme::Deme,tooken::Tooken,signer::Signer)

    socket = connect(cc.tookenport)

    dh = DHsym(deme,signer)

    key, id = diffiehellman(socket,dh)

    @assert id in cc.serverid

    securesocket = deme.cypher.secureio(socket,key)

    serialize(securesocket,tooken)
end

addtooken(cc::CertifierConfig,deme::Deme,tooken::Int,signer::Signer) = addtooken(cc,deme,Tooken(tooken),signer)

function certify(port,serverid::ID,deme::Deme,cert::Certificate{Profile},tooken::Tooken; sendcopy=false)
    socket = connect(port)

    dh = DHasym(deme)

    key, keyid = diffiehellman(socket,dh)

    @assert keyid in serverid

    securesocket = deme.cypher.secureio(socket,key)
    serialize(securesocket,TookenID(cert,tooken))

    sendcopy && (return deserialize(securesocket,Certificate{ID}))
end

certify(port,serverid::ID,deme::Deme,cert::Certificate{Profile},tooken::Int; kwargs...) = certify(port,serverid,deme,cert,Tooken(tooken); kwargs...)

function str(config::Dict)
    io = IOBuffer()
    TOML.print(io, config)
    return String(take!(io))
end

function ticket(deme::DemeSpec,port,server::ID,tooken::Int)
    config = Dict("demespec"=>Dict(deme),"port"=>Dict(port),"tooken"=>tooken,"server"=>string(server,base=16))
    return str(config)
end

function ticket(deme::DemeSpec,port::Int,server::ID,tooken::Int)
    config = Dict("demespec"=>Dict(deme),"port"=>Dict("type"=>"Int","port"=>port),"tooken"=>tooken,"server"=>string(server,base=16))
    return str(config)
end

### To use this function one is supposed to know
### How to create a identity

### I could make a function which would assume that ID type would be of a type ID. A function for Port could be passed. The problem is that we could not allow

struct Port
    port::Int
    ip::Union{IPv4,IPv6}
end

import Sockets.connect
connect(port::Port) = connect(port.ip,port.port)

function Dict(port::Port)
    dict = Dict()
    if port.ip isa IPv4
        dict["type"]="IPv4"
    else
        dict["type"]="IPv6"
    end

    dict["port"] = port.port
    dict["ip"] = string(port.ip)

    return dict
end

struct RegPort{T} end
RegPort(type::AbstractString) = RegPort{Symbol(type)}()

port(::RegPort{:IPv4},config::Dict) = Port(config["port"],IPv4(config["ip"]))
port(::RegPort{:IPv6},config::Dict) = Port(config["port"],IPv6(config["ip"]))
port(::RegPort{:Int},config::Dict) = config["port"]

function register(invite::Dict,profile::Profile; account="") where T <: AbstractID

    demespec = DemeSpec(invite["demespec"])
    save(demespec)
    deme = Deme(demespec)

    portdict = invite["port"]

    @assert haskey(portdict,"type") "The type of the port must be defined"

    type = RegPort(portdict["type"])
    regport = port(type,portdict)

    serverid = ID(invite["server"],base=16)

    tooken = invite["tooken"]

    member = Signer(deme,account * "/member")
    profilecert = Certificate(profile,member)

    certify(regport,serverid,deme,profilecert,tooken)
end

register(invite::AbstractString,profile::Profile; kwargs...) = register(TOML.parse(invite),profile; kwargs...)

struct SMTPConfig
    url::AbstractString
    email::AbstractString
    password::Union{AbstractString,Nothing}
end

function SMTPConfig()
    println("Email:")
    email = readline()

    defaulturl = "smtps://smtp.gmail.com:465"
    println("SMTP url [$defaulturl]:")
    url = readline()
    url=="" && (url = defaulturl)

    println("Password:")
    password = readline()

    SMTPConfig(url,email,password)

end

### Need to think about this part
### Perhaps I could add additional metadata as dictionary!

function sendinvite(config::CertifierConfig,deme::Deme,to::AbstractString,from::SMTPConfig,maintainer::Signer; init=nothing)
    ### This would register a tooken with the system

    tooken = rand(2^62:2^63-1)
    addtooken(config,deme,tooken,maintainer)

    opt = SendOptions(isSSL = true, username = from.email, passwd = from.password)

    t = ticket(deme.spec,config.certifierport,config.serverid,tooken)

    body = """
    From: $(from.email)
    To: $to
    Subject: Invitation to $(deme.spec.name)

    ########### Ticket #############
    $t
    ################################
    """

    send(from.url, [to], from.email, IOBuffer(body), opt)
end

sendinvite(config::CertifierConfig,deme::Deme,to::Vector{T},maintainer::Signer; init=nothing) where T<:AbstractString = sendinvite(config,deme,to,SMTPConfig(),maintainer; init=init)

export CertifierConfig, Certifier, addtooken, certify, ticket, register, sendinvite, SMTPConfig

end


text:
using Documenter, Perturbo

makedocs(;
    modules=[Perturbo],
    format=Documenter.HTML(prettyurls = false),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jinjianzhou/Perturbo.jl/blob/{commit}{path}#L{line}",
    sitename="Perturbo.jl",
    authors="Jin-Jian Zhou <jjchou.comphy@gmail.com>",
    assets=String[],
)

deploydocs(;
    repo="github.com/jinjianzhou/Perturbo.jl",
)


text:
# HistoryRecorder
# maintained by @zsunberg

"""
A simulator that records the history for later examination

The simulation will be terminated when either
1) a terminal state is reached (as determined by `isterminal()` or
2) the discount factor is as small as `eps` or
3) max_steps have been executed

Keyword Arguments:
    - `rng`: The random number generator for the simulation
    - `capture_exception::Bool`: whether to capture an exception and store it in the history, or let it go uncaught, potentially killing the script
    - `show_progress::Bool`: show a progress bar for the simulation
    - `eps`
    - `max_steps`
    - `sizehint::Int`: the expected length of the simulation (for preallocation)

Usage (optional arguments in brackets):
    hr = HistoryRecorder()
    history = simulate(hr, pomdp, policy, [updater [, init_belief [, init_state]]])
"""
mutable struct HistoryRecorder <: Simulator
    rng::AbstractRNG

    # options
    capture_exception::Bool
    show_progress::Bool

    # optional: if these are null, they will be ignored
    max_steps::Nullable{Any}
    eps::Nullable{Any}
    sizehint::Nullable{Integer}

    # DEPRECATED
    initial_state::Nullable{Any}
end

# This is the only stable constructor
function HistoryRecorder(;rng=MersenneTwister(rand(UInt32)),
                          initial_state=Nullable{Any}(),
                          eps=Nullable{Any}(),
                          max_steps=Nullable{Any}(),
                          sizehint=Nullable{Integer}(),
                          capture_exception=false,
                          show_progress=false)
    if !isnull(initial_state)
        warn("The initial_state argument for HistoryRecorder is deprecated. The initial state should be specified as the last argument to simulate(...).")
    end
    return HistoryRecorder(rng, capture_exception, show_progress, max_steps, eps, sizehint, initial_state)
end

@POMDP_require simulate(sim::HistoryRecorder, pomdp::POMDP, policy::Policy) begin
    @req updater(::typeof(policy))
    up = updater(policy)
    @subreq simulate(sim, pomdp, policy, up)
end

@POMDP_require simulate(sim::HistoryRecorder, pomdp::POMDP, policy::Policy, bu::Updater) begin
    @req initial_state_distribution(::typeof(pomdp))
    dist = initial_state_distribution(pomdp)
    @subreq simulate(sim, pomdp, policy, bu, dist)
end

function simulate(sim::HistoryRecorder, pomdp::POMDP, policy::Policy, bu::Updater=updater(policy))
    dist = initial_state_distribution(pomdp)
    return simulate(sim, pomdp, policy, bu, dist)
end

@POMDP_require simulate(sim::HistoryRecorder, pomdp::POMDP, policy::Policy, bu::Updater, dist::Any) begin
    P = typeof(pomdp)
    S = state_type(pomdp)
    A = action_type(pomdp)
    O = obs_type(pomdp)
    if isnull(sim.initial_state)
        @req rand(::typeof(sim.rng), ::typeof(dist))
    end
    @req initialize_belief(::typeof(bu), ::typeof(dist))
    @req isterminal(::P, ::S)
    @req discount(::P)
    @req generate_sor(::P, ::S, ::A, ::typeof(sim.rng))
    b = initialize_belief(bu, dist)
    B = typeof(b)
    @req action(::typeof(policy), ::B)
    @req update(::typeof(bu), ::B, ::A, ::O)
end

function simulate{S,A,O}(sim::HistoryRecorder,
                           pomdp::POMDP{S,A,O},
                           policy::Policy,
                           bu::Updater,
                           initial_state_dist::Any,
                           initial_state::Any=get_initial_state(sim, initial_state_dist)
                          )

    initial_belief = initialize_belief(bu, initial_state_dist)
    max_steps = get(sim.max_steps, typemax(Int))
    if !isnull(sim.eps)
        max_steps = min(max_steps, ceil(Int,log(get(sim.eps))/log(discount(pomdp))))
    end
    sizehint = get(sim.sizehint, min(max_steps, 1000))

    # aliases for the histories to make the code more concise
    sh = sizehint!(Vector{S}(0), sizehint)
    ah = sizehint!(Vector{A}(0), sizehint)
    oh = sizehint!(Vector{O}(0), sizehint)
    bh = sizehint!(Vector{typeof(initial_belief)}(0), sizehint)
    rh = sizehint!(Vector{Float64}(0), sizehint)
    ih = sizehint!(Vector{Any}(0), sizehint)
    aih = sizehint!(Vector{Any}(0), sizehint)
    uih = sizehint!(Vector{Any}(0), sizehint)
    exception = Nullable{Exception}()
    backtrace = Nullable{Any}()

    push!(sh, initial_state)
    push!(bh, initial_belief)

    if sim.show_progress
        if isnull(sim.max_steps) && isnull(sim.eps)
            error("If show_progress=true in a HistoryRecorder, you must also specify max_steps or eps.")
        end
        prog = Progress(max_steps, "Simulating..." )
    end

    disc = 1.0
    step = 1

    try
        while !isterminal(pomdp, sh[step]) && step <= max_steps
            a, ai = action_info(policy, bh[step])
            push!(ah, a)
            push!(aih, ai)

            sp, o, r, i = generate_sori(pomdp, sh[step], ah[step], sim.rng)

            push!(sh, sp)
            push!(oh, o)
            push!(rh, r)
            push!(ih, i)

            bp, ui = update_info(bu, bh[step], ah[step], oh[step])
            push!(bh, bp)
            push!(uih, ui)

            step += 1

            if sim.show_progress
                next!(prog)
            end
        end
    catch ex
        if sim.capture_exception
            exception = Nullable{Exception}(ex)
            backtrace = Nullable{Any}(catch_backtrace())
        else
            rethrow(ex)
        end
    end

    if sim.show_progress
        finish!(prog)
    end

    return POMDPHistory(sh, ah, oh, bh, rh, ih, aih, uih, discount(pomdp), exception, backtrace)
end

@POMDP_require simulate(sim::HistoryRecorder, mdp::MDP, policy::Policy) begin
    if isnull(sim.initial_state)
        @req initial_state(::typeof(mdp), ::typeof(sim.rng))
    end
    init_state = get(sim.initial_state, initial_state(mdp, sim.rng))
    @subreq simulate(sim, mdp, policy, init_state)
end

@POMDP_require simulate(sim::HistoryRecorder, mdp::MDP, policy::Policy, initial_state::Any) begin
    P = typeof(mdp)
    S = state_type(mdp)
    A = action_type(mdp)
    @req isterminal(::P, ::S)
    @req action(::typeof(policy), ::S)
    @req generate_sr(::P, ::S, ::A, ::typeof(sim.rng))
    @req discount(::P)
end

function simulate{S,A}(sim::HistoryRecorder,
                       mdp::MDP{S,A}, policy::Policy,
                       init_state::S=get_initial_state(sim, mdp))

    max_steps = get(sim.max_steps, typemax(Int))
    if !isnull(sim.eps)
        max_steps = min(max_steps, ceil(Int,log(get(sim.eps))/log(discount(mdp))))
    end
    sizehint = get(sim.sizehint, min(max_steps, 1000))

    # aliases for the histories to make the code more concise
    sh = sizehint!(Vector{S}(0), sizehint)
    ah = sizehint!(Vector{A}(0), sizehint)
    rh = sizehint!(Vector{Float64}(0), sizehint)
    ih = sizehint!(Vector{Any}(0), sizehint)
    aih = sizehint!(Vector{Any}(0), sizehint)
    exception = Nullable{Exception}()
    backtrace = Nullable{Any}()

    if sim.show_progress
        prog = Progress(max_steps, "Simulating..." )
    end

    push!(sh, init_state)

    disc = 1.0
    step = 1

    try
        while !isterminal(mdp, sh[step]) && step <= max_steps
            a, ai = action_info(policy, sh[step])
            push!(ah, a)
            push!(aih, ai)

            sp, r, i = generate_sri(mdp, sh[step], ah[step], sim.rng)

            push!(sh, sp)
            push!(rh, r)
            push!(ih, i)

            disc *= discount(mdp)
            step += 1

            if sim.show_progress
                next!(prog)
            end
        end
    catch ex
        if sim.capture_exception
            exception = Nullable{Exception}(ex)
            backtrace = Nullable{Any}(catch_backtrace())
        else
            rethrow(ex)
        end
    end

    if sim.show_progress
        finish!(prog)
    end

    return MDPHistory(sh, ah, rh, ih, aih, discount(mdp), exception, backtrace)
end

function get_initial_state(sim::Simulator, initial_state_dist)
    if isnull(sim.initial_state)
        return rand(sim.rng, initial_state_dist)
    else
        return get(sim.initial_state)
    end
end

function get_initial_state(sim::Simulator, mdp::Union{MDP,POMDP})
    if isnull(sim.initial_state)
        return initial_state(mdp, sim.rng)
    else
        return get(sim.initial_state)
    end
end


text:
#!/usr/bin/env julia

using Test, Printf, EAGO, MathOptInterface, Cbc, JuMP, Ipopt

const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

include(joinpath(@__DIR__, "moit_tests.jl"))
include(joinpath(@__DIR__, "minlp_tests.jl"))

#include(joinpath(@__DIR__, "branch_bound.jl"))
#include(joinpath(@__DIR__, "domain_reduction.jl"))
#include(joinpath(@__DIR__, "optimizer.jl"))
#include(joinpath(@__DIR__, "script_optimizer.jl"))
#include(joinpath(@__DIR__, "semiinfinite.jl"))

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module ProfitLoss

import ..Clarus

function hypothetical(;params...)
  return Clarus.api_request("ProfitLoss","Hypothetical";params...)
end

function predict(;params...)
  return Clarus.api_request("ProfitLoss","Predict";params...)
end

end #MODULE-END


text:
module Keldysh

using LinearAlgebra
using HDF5, QuadGK, Elliptic

# utility functions
export fermi

# Branch functions
export BranchEnum, forward_branch, backward_branch, imaginary_branch, BranchPoint, Branch

# Contour functions
export AbstractContour, FullContour, KeldyshContour, ImaginaryContour, twist, heaviside, nbranches

# TimeGrid functions
export AbstractTimeGrid, FullTimeGrid, KeldyshTimeGrid, ImaginaryTimeGrid, TimeGridPoint, TimeDomain, integrate, realtimes, imagtimes

# TimeGF functions
export fermionic, bosonic, norbitals, is_scalar, AbstractTimeGF, GenericTimeGF, FullTimeGF, TimeInvariantFullTimeGF, KeldyshTimeGF, TimeInvariantKeldyshTimeGF, ImaginaryTimeGF

# DOS functions
export AbstractDOS, DOS, SingularDOS, DeltaDOS, flat_dos, gaussian_dos, bethe_dos, chain_dos, square_dos

include("util.jl")
include("branch.jl")
include("dos.jl")
include("contour.jl")
include("time_grid.jl")
include("storage.jl")
include("interp.jl")
include("gf.jl")
include("alps.jl")
include("nessi.jl")

end # module


text:
using Test
using SimpleTools

sq = Dict{Int,Int}()
for k = 1:10
    sq[k] = k * k
end

rt = Dict{Int,Int}()
for k = 1:10
    rt[k*k] = k
end

f = rt * sq
for k = 1:10
    @test f[k] == k
end

y1 = [x^2 for x = 0:0.01:3]
y2 = [mod(y, 1) for y in y1]
make_continuous!(y2, 1)
@test y1 == y2

@test flush_print(23, 5) == "   23"
@test flush_print(23, 5, false) == "23   "

A = dcat(1, 1, 1)
M = [1 0 0; 0 1 0; 0 0 1]
@test A == M

@test sage(A) == nothing


text:
using PyPlot

rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
rcParams["text.usetex"] = true
rcParams["mathtext.fontset"]  = "cm";
rcParams["font.family"]  = "serif";
rcParams["font.serif"]  = "cm";
plt.ioff()

function phi(t::Real)

    return (2 * pi)^(-0.5) * exp(-(t^2) / 2)
end

function generate_f(evals, p)

    f = (p[1]^2 + p[3]^2) * phi.(evals .- 1)
    f += p[2]*(2 - p[2]) * phi.(evals)
    f += 2*p[1]*p[3] * phi.(evals .+ 1)

    return f
end

function generate_Var_f_given_A(evals, p)

    Var_f_given_A = (p[1] + p[3])^2 * p[2] * phi.(evals).^2
    Var_f_given_A += p[3] * (p[3] * phi.(evals .- 1) + p[1] * phi.(evals .+ 1)).^2
    Var_f_given_A += p[1] * (p[1] * phi.(evals .- 1) + p[3] * phi.(evals .+ 1)).^2
    Var_f_given_A -= ((p[1]^2 + p[3]^2) * phi.(evals .- 1) +
        (p[1] + p[3]) * p[2] * phi.(evals) + 2 * p[1] * p[3] * phi.(evals .+ 1)).^2

    return Var_f_given_A

end

function plot_confidence_intervals(ax, x, y_lower, y_upper, num_intervals,
                                   line_color, line_width, line_style)

    inds = range(1, stop=length(x), length=2*num_intervals+1)
    inds = [inds[i] for i in 1:length(inds) if i%2 == 0]
    inds = Int.(round.(inds))
    @assert size(inds, 1) == num_intervals
    xr = x[inds]
    ylr = y_lower[inds]
    yur = y_upper[inds]

    ax.vlines(
        xr[1], ylr[1], yur[1], linewidth=line_width,
        linestyle=line_style, color=line_color)

    for i in 1:num_intervals

        ax.vlines(
            xr[i], ylr[i], yur[i], linewidth=line_width,
            linestyle=line_style, color=line_color)

        marker = PyPlot.matplotlib.path.Path([[-1,0], [1,0], [-1,0], [1,0]])

        ax.plot(
            xr[i], ylr[i], linewidth=0.0,
            color=line_color, marker=marker,
            markersize=4.0)

        ax.plot(
            xr[i], yur[i], linewidth=0.0,
            color=line_color, marker=marker,
            markersize=4.0)
    end
end


text:
## MAGIJ: Matrix Assistance to Gensys In Julia
## WARNING: Higlhy experimental
##

##TODO: Allow keys different of numbers in dict (or, how to index dict?)
## TODO: Automatic Creation of expectations shock and foward variables

using DynamicPolynomials
using Plots

#Parameters creation

bet = 0.99
sig = 1
phi = 1
alfa = 1/3
epsilon = 6
theta = 2/3
phi_pi = 1.5
phi_y = 0.5/4
rho_v = 0.5

THETA = (1-alfa)/(1-alfa+alfa*epsilon)
lamb = (1-theta)*(1-bet*theta)/theta*THETA
kappa = lamb*(sig+(phi+alfa)/(1-alfa))

#Variables

@polyvar pi pi_e y y_e i v v_l ep

#variables = [pi pi_e y y_e v v_l ep]

#Equations of the model

eq1 = pi - bet*pi_e - kappa*y
eq2 = y + 1/sig*(i - pi_e) - y_e
eq3 = i - phi_pi*pi - phi_y*y - v_l
eq4 = v - rho_v*v_l - ep

#Dictionary with all the equations

model = Dict(1=>eq1,2=>eq2,3=>eq3,4=>eq4)

model = model2gensys(model,[ep])

#What variables appear in which model

varis = Dict()

for i in keys(model)
    merge!(varis,Dict(i=>variables(model[i])))
end

#Defyning the auxilary variables

@polyvar nu_pi nu_y pi_f y_f

pi_e_aux = pi_f - nu_pi
y_e_aux = y_f - nu_y

foward_variables = [pi_f y_f v]
backward_variables = [pi y v_l i]
shocks = [ep]
expec_errors = [nu_pi nu_y]

f_string = [repr(i) for i in foward_variables]
expec_vars_test = match.(r"_f",f_string)
expec_vars_test = expec_vars_test .!= nothing
expec_vars = f_string[expec_vars_test]
stripped_vars = [rstrip(vv,['_','f']) for vv in expec_vars]
foward_vars = string.(stripped_vars,"_f")

bb = "@polyvar"

for aa in foward_vars
    global bb = join([bb,aa]," ")
end

for i in keys(model)
    if sum(pi_e .== varis[i]) == 1
        model[i] = subs(model[i],pi_e => pi_e_aux)
    end
    if sum(y_e .== varis[i]) == 1
        model[i] = subs(model[i],y_e => y_e_aux)
    end
end

vari_all = Vector()

for i in keys(varis)
    union!(vari_all,variables(model[i]))
end

#Building each matrix

Gamma1 = zeros(length(model),length(model))
Gamma0 = zeros(length(model),length(model))
Psi = zeros(length(model),length(shocks))
Pi =  zeros(length(model),length(expec_errors))

#Gamma 0

for j in 1:length(foward_variables)
    for i in keys(model)
        Gamma0[i,j] = -coefficient(model[i],foward_variables[j])
    end
end

#Gamma 1

for j in 1:length(backward_variables)
    for i in keys(model)
        Gamma1[i,j] = coefficient(model[i],backward_variables[j])
    end
end

#Psi

for j in 1:length(shocks)
    for i in keys(model)
        Psi[i,j] = coefficient(model[i],shocks[j])
    end
end

#Pi

for j in 1:length(expec_errors)
    for i in keys(model)
        Pi[i,j] = coefficient(model[i],expec_errors[j])
    end
end

sol = gensys(model.Gamma0,model.Gamma1,model.Psi,model.Pi)

fris = irf(sol,15,0.25)

plot(4*fris[:,1],w = 3)
plot!(4*irfs_true[2:15,1], line = :dash,w = 3)

plot(fris[:,2],w = 3)
plot!(irfs_true[2:15,2], line = :dash,w = 3)

plot(4*fris[:,4],w = 3)
plot!(4*irfs_true[2:15,3], line = :dash,w = 3)

plot(fris[:,3],w = 3)
plot!(irfs_true[2:15,4], line = :dash,w = 3)

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
function count_loc(dir)
    # we `cd` so that we get relative paths in the `tokei` output.
    # This makes it easy to process later, since we have uniform filepaths
    json = cd(dir) do
        try
            tokei() do exe
                JSON3.read(read(`$exe --output json .`))
            end
        catch e
            @error e
            missing
        end
    end
    return make_loc_table(json)
end

const LoCTableEltype = @NamedTuple{directory::String, language::Symbol, sublanguage::Union{Nothing, Symbol}, files::Int, code::Int, comments::Int, blanks::Int}

function make_loc_table(json)
    table = LoCTableEltype[]
    ismissing(json) && return table
    for (language, language_loc) in pairs(json)
        # we want to count lines of code per toplevel directory, per language, and per sublanguage (e.g. for Julia inside of Markdown)
        counts = Dict{@NamedTuple{directory::String, sublanguage::Union{Nothing, Symbol}}, @NamedTuple{files::Int, code::Int, comments::Int, blanks::Int}}()
        language == :Total && continue # skip the fake `Total` language
        language_loc.inaccurate && continue # skip if it's marked `inaccurate` (not sure when that happens?)
        for report in language_loc.reports
            # `splitpath(relative_path)[1] == ".", so `splitpath(relative_path)[2]` gives us the toplevel directory or filename
            directory = splitpath(report.name)[2]
            loc_update!(counts, (; directory, sublanguage=nothing), report.stats)
            for (sublanguage, sublanguage_loc) in report.stats.blobs
                loc_update!(counts,  (; directory, sublanguage), sublanguage_loc)
            end
        end
        for ((directory, sublanguage), count) in pairs(counts)
            push!(table, (; directory, language, sublanguage, count.files, count.code, count.comments, count.blanks))
        end
    end
    d_count = counts_by_col(table, :directory)
    l_count = counts_by_col(table, :language)
    sl_count = counts_by_col(table, :sublanguage)
    sort!(table, by = row->(d_count[row.directory], l_count[row.language], sl_count[row.sublanguage]), rev=true)
    return table
end

function counts_by_col(table, col)
    vals = unique(getproperty(row, col) for row in table)
    return Dict(val => sum(row.code for row in table if getproperty(row, col)==val) for val in vals)
end

function loc_update!(d, key, new)
    prev = get!(d, key, (; files=0, code=0, comments=0, blanks=0))
    d[key] = (; files = prev.files + 1, code = prev.code + new.code, comments = prev.comments + new.comments, blanks = prev.blanks + new.blanks )
end

count_julia_loc(table, dir) = sum(row.code for row in table if row.directory == dir && row.language == :Julia; init=0)


text:
"""
    ItemT

Is a new type introduced by our package for control
the Item Object. This object is MUTABLE and has four
arguments:

    length::Float64  is the length of the item

    height::Float64  is the height of the item

    posX::Float64  position X of the item in container

    posY::Float64  position X of the item in container

    ItemT(lenght,height,posX,posY)

# Examples
```julia-repl
julia> ItemT(1.0, 2.0, 0.0, 0.0)
```
"""
mutable struct ItemT <: Number
    length::Float64 # is the length of the item
    height::Float64 # is the height of the item
    posX::Float64 # position X of the item in container
    posY::Float64 # position X of the item in container
end

abstract type ContainerType end
"""
    ContainerT

Is a new type introduced by our package for control
the Container Object. This object is MUTABLE and has three
arguments:

    length::Float64  is the length of the container

    height::Float64  is the height of the container

    content::Array{ItemT} is an array of items

    ContainerT(lenght,height,content)

# Examples
```julia-repl
julia> ContainerT(1.0, 2.0, 0.0, [ItemT(1.0, 2.0, 0.0, 0.0),ItemT(2.0, 1.5, 1.0, 1.0)])
```
"""
mutable struct ContainerT <: ContainerType
    length::Float64
    height::Float64
    content::Array{ItemT}
end

importall Base.Operators
function Base.copy(A::Any)
    if typeof(A)==ContainerT
        return ContainerT((A.length),(A.height),(A.content))
    end
    if typeof(A)==Array{ContainerT,1}
        m=length(A)
        B=Array{ContainerT}(m)
        for i=1:m
                B[i]=ContainerT((A[i].length),(A[i].height),(A[i].content))
        end
        return B
    end
end

export ContainerType,ContainerT,ItemT,copy

max_stars_repo_name: evcastelani/FFDpacking

max_issues_repo_name: evcastelani/FFDpacking

max_forks_repo_name: evcastelani/FFDpacking


text:
using Dates

# discretizations
struct QuadraticExponentialMartingale <: AbstractHestonDiscretization end
struct FullTruncation <: AbstractHestonDiscretization end

struct HestonProcess{Y1 <: YieldTermStructure, Y2 <: YieldTermStructure, D <: AbstractDiscretization, DH <: AbstractHestonDiscretization} <: AbstractHestonProcess
  s0::Quote
  riskFreeRate::Y1
  dividendYield::Y2
  v0::Float64
  kappa::Float64
  theta::Float64
  sigma::Float64
  rho::Float64
  disc::D
  hestonDisc::DH
end

function HestonProcess(riskFreeRate::Y1, dividendYield::Y2, s0::Quote, v0::Float64, kappa::Float64, theta::Float64, sigma::Float64, rho::Float64,
                      d::D = QuadraticExponentialMartingale()) where {Y1 <: YieldTermStructure, Y2 <: YieldTermStructure, D <: AbstractHestonDiscretization}
  disc = EulerDiscretization()
  return HestonProcess{Y1, Y2, EulerDiscretization, D}(s0, riskFreeRate, dividendYield, v0, kappa, theta, sigma, rho, disc, d)
end

get_time(process::AbstractHestonProcess, d::Date) = year_fraction(get_risk_free_rate(process).dc, reference_date(get_risk_free_rate(process)), d)

get_risk_free_rate(hp::HestonProcess) = hp.riskFreeRate
get_dividend_yield(hp::HestonProcess) = hp.dividendYield
get_s0(hp::HestonProcess) = hp.s0
get_v0(hp::HestonProcess) = hp.v0
get_kappa(hp::HestonProcess) = hp.kappa
get_theta(hp::HestonProcess) = hp.theta
get_sigma(hp::HestonProcess) = hp.sigma
get_rho(hp::HestonProcess) = hp.rho


text:
function Tait(n::Integer)
    circ = tait_id()
    id = tait_id()
    for _ in 1:(n-1)
        tensor_product!(circ, id)
    end
    return circ
end

using ZXCalculus: push_gate!

function ZXCalculus.push_gate!(q::Tait, ::Val{:H}, loc)
    h = tait_hadamard()
    contract!(q, h, [q.outputs[loc]], copy(h.inputs))
    return q
end

function ZXCalculus.push_gate!(q::Tait, ::Val{:Rz}, loc, p)
    g = tait_rz(p*im)
    contract!(q, g, [q.outputs[loc]], copy(g.inputs))
    return q
end
ZXCalculus.push_gate!(q::Tait, ::Val{:shift}, loc, p) = ZXCalculus.push_gate!(q, Val(:Rz), loc, p)

function ZXCalculus.push_gate!(q::Tait, ::Val{:Rx}, loc, p)
    g = tait_rx(p*im)
    contract!(q, g, [q.outputs[loc]], copy(g.inputs))
    return q
end

function ZXCalculus.push_gate!(q::Tait, ::Val{:SWAP}, loc1, loc2)
    l1, l2 = (min(loc1, loc2), max(loc1, loc2))
    g = tait_swap()
    for i = l1:(l2-1)
        contract!(q, g, q.outputs[i:i+1], copy(g.inputs))
    end
    for i = (l2-2):-1:l1
        contract!(q, g, q.outputs[i:i+1], copy(g.inputs))
    end
    return q
end

function ZXCalculus.push_gate!(q::Tait, ::Val{:CNOT}, loc, ctrl)
    g = tait_cnot(ctrl, loc)
    contract!(q, g, q.outputs[min(loc, ctrl):max(loc, ctrl)], copy(g.inputs))
    return q
end

function ZXCalculus.push_gate!(q::Tait, ::Val{:CZ}, loc, ctrl)
    g = tait_cz(ctrl, loc)
    contract!(q, g, q.outputs[min(loc, ctrl):max(loc, ctrl)], copy(g.inputs))
    return q
end


text:
module Bigleaf

#using DocStringExtensions
using Optim
using FillArrays
using DataFrames
using Tables: columntable
using Dates, TimeZones
using Pipe
using AstroLib
using Suppressor
using Missings
using Statistics, StatsBase # mean, rle
using PaddedViews, StaticArrays
using Infiltrator
using Parameters

export ConductanceMethod, Thom1972, Choudhury1988, Su2001, ConstantKB1
export ResistanceMethod, ResistanceWindZr, ResistanceWindProfile
export RoughnessMethod, Roughness_wind_profile, RoughnessCanopyHeight, RoughnessCanopyHeightLAI
export ETMethod, PriestleyTaylor, PenmanMonteith
export EsatMethod, Sonntag1990, Alduchov1996, Allen1998
export StabilityCorrectionMethod, Dyer1970, Businger1971, NoStabilityCorrection
export SurfaceConductanceMethod, FluxGradient, InversePenmanMonteith

export frac_hour, moving_average, get_nonoverlapping_periods, set_datetime_ydh!
export BigleafConstants
export Esat_slope, Esat_from_Tair, Esat_from_Tair_deriv,
     LE_to_ET, ET_to_LE, ms_to_mol, mol_to_ms, VPD_to_rH, rH_to_VPD,
     e_to_rH, VPD_to_e, e_to_VPD, e_to_q, q_to_e, q_to_VPD, VPD_to_q,
     Rg_to_PPFD, PPFD_to_Rg, kg_to_mol, umolCO2_to_gC, gC_to_umolCO2
export air_density, pressure_from_elevation, psychrometric_constant,
    latent_heat_vaporization, virtual_temp, kinematic_viscosity,
    dew_point_from_e, dew_point,
    wetbulb_temp_from_e_Tair_gamma, wetbulb_temp
export calc_sun_position_MOD, calc_sun_position_hor
export potential_radiation, extraterrestrial_radiation, get_datetime_for_doy_hour
export potential_ET, potential_ET!, equilibrium_imposed_ET, equilibrium_imposed_ET!
export setinvalid_range!, setinvalid_qualityflag!,
    setinvalid_nongrowingseason!, get_growingseason, setinvalid_afterprecip!
export decoupling, surface_conductance, aerodynamic_conductance
export compute_Gb, compute_Gb!, add_Gb, add_Gb!,
    Gb_Thom, Gb_Choudhury, Gb_Su, Gb_constant_kB1,
    compute_Gb_quantities, compute_Gb_quantities!
export wind_profile
export Monin_Obukhov_length, Monin_Obukhov_length!, stability_parameter,
    stability_parameter!, stability_correction, stability_correction!,
    roughness_parameters, Reynolds_Number
export aerodynamic_conductance!, add_Ga, add_Ga!, compute_Ram, compute_Ram!, roughness_length_heat
export surface_conductance, surface_conductance!

# shorthand type: either DataFrame or column table (Tables.columntable)
DFTable = Union{NamedTuple,AbstractDataFrame}

include("util.jl")
include("bigleaf_constants.jl")
include("unit_conversions.jl")
include("meteorological_variables.jl")
include("sun_position.jl")
include("potential_radiation.jl")
include("evapotranspiration.jl")
include("filter_data.jl")
include("stability_correction.jl")
include("surface_roughness.jl")
include("boundary_layer_conductance.jl")
include("aerodynamic_conductance.jl")
include("surface_conductance.jl")

end


text:
using SymEngineExtras
using Base.Test

include("test-patternmatch.jl")
#include("test-simplify.jl")


text:
module NURBS

include("Utils.jl")
include("Bases.jl")

end  # module NURBS


text:
# This file was generated, do not modify it. # hide
fitted_params(tree) |> pprint


text:
using RoME
# , IncrementalInference, Distributions
using Test

# const TU = TransformUtils

@testset "test DynPose2 and velocity..." begin

global N = 100
global fg = initfg()

# add first pose locations
addVariable!(fg, :x0, DynPose2(ut=0))

# Prior factor as boundary condition
global pp0 = DynPose2VelocityPrior(MvNormal(zeros(3), Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([10.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x0;], pp0)

# initialize the first pose
IncrementalInference.doautoinit!(fg, [getVariable(fg,:x0);])

addVariable!(fg, :x1, DynPose2(ut=1000_000))

# conditional likelihood between Dynamic Point2
global dp2dp2 = VelPose2VelPose2(MvNormal([10.0;0;0], Matrix(Diagonal([0.01;0.01;0.001].^2))),
                          MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x0;:x1], dp2dp2)

# getVal(fg,:x0)
global pts = approxConv(fg, :x0x1f1, :x1)

# Graphs.plot(fg.g)
# ensureAllInitialized!(fg)

tree, smt, hist = solveTree!(fg)
# global tree = batchSolve!(fg, N=N)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

global X1 = getVal(fg, :x1)

@test 0.9*N <= sum(abs.(X1[1,:] .- 10.0) .< 0.75)
@test 0.9*N <= sum(abs.(X1[2,:] .- 0.0) .< 0.75)
# @show TU.wrapRad.(X1[3,:])
@test 0.8*N <= sum(abs.(TU.wrapRad.(X1[3,:]) .- 0.0) .< 0.25)
# @warn "wrapRad issue, accepting 80% as good enough until issue JuliaRobotics/RoME.jl#90 is fixed."
@test 0.9*N <= sum(abs.(X1[4,:] .- 10.0) .< 0.5)
@test 0.9*N <= sum(abs.(X1[5,:] .- 0.0) .< 0.5)

# using RoMEPlotting
# # plotLocalProduct(fg, :x10, dims=[1;2])
# # drawPoses(fg)
# xx1 = marginal(getKDE(fg, :x1),[1;2;3])
# plotPose(Pose2(), [xx1])
# plotPose(fg, [:x1], levels=1, show=false)
#
# plotKDE(marginal(getKDE(fg, :x1),[4;5]), levels=5)

end

@testset "test distribution compare functions..." begin

global mu = randn(6)
global mv1 = MvNormal(deepcopy(mu), Matrix{Float64}(LinearAlgebra.I, 6,6))
global mv2 = MvNormal(deepcopy(mu), Matrix{Float64}(LinearAlgebra.I, 6,6))
global mv3 = MvNormal(randn(6), Matrix{Float64}(LinearAlgebra.I, 6,6))
@test RoME.compare(mv1, mv2)
@test !RoME.compare(mv1, mv3)
@test !RoME.compare(mv2, mv3)

end

@testset "test DynPose2 packing converters..." begin

global pp0 = DynPose2VelocityPrior(MvNormal(zeros(3), Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([10.0;0], Matrix(Diagonal([0.1; 0.1].^2))))

global pp = convert(PackedDynPose2VelocityPrior, pp0)
global ppu = convert(DynPose2VelocityPrior, pp)

@test RoME.compare(pp0, ppu)

global dp2dp2 = VelPose2VelPose2(MvNormal([10.0;0;0], Matrix(Diagonal([0.01;0.01;0.001].^2))),
                          MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))

global pp = convert(PackedVelPose2VelPose2, dp2dp2)
global ppu = convert(VelPose2VelPose2, pp)

@test RoME.compare(dp2dp2, ppu)

end

@testset "test many DynPose2 chain stationary and 'pulled'..." begin

global N = 100
global fg = initfg()

# add first pose locations
addVariable!(fg, :x0, DynPose2(ut=0))

# Prior factor as boundary condition
global pp0 = DynPose2VelocityPrior(MvNormal(zeros(3), Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x0;], pp0)

global sym = :x0
global k = 0
for sy in Symbol[Symbol("x$i") for i in 1:10]

global k+=1
addVariable!(fg, sy, DynPose2(ut=1000_000*k))

# conditional likelihood between Dynamic Point2
global dp2dp2 = VelPose2VelPose2(MvNormal([0.0;0;0], Matrix(Diagonal([1.0;0.1;0.001].^2))),
                          MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [sym;sy], dp2dp2)
global sym =sy

end

global x5 = KDE.getKDEMean(getKDE(getVariable(fg, :x5)))

@test abs(x5[1]) < 1.25
@test abs(x5[2]) < 1.25
@test abs(TU.wrapRad(x5[3])) < 0.4
@test abs(x5[4]) < 0.5
@test abs(x5[5]) < 0.5

ensureAllInitialized!(fg)

global x10 = KDE.getKDEMean(getKDE(getVariable(fg, :x10)))

@test abs(x10[1]) < 1.25
@test abs(x10[2]) < 1.25
@test abs(TU.wrapRad(x10[3])) < 0.4
@test abs(x10[4]) < 0.5
@test abs(x10[5]) < 0.5

# using RoMEPlotting
# drawPoses(fg)
# plotPose(fg, [:x10])

tree, smt, hist = solveTree!(fg)
# batchSolveR!(fg, N=N)
# tree = wipeBuildNewTree!(fg)
# inferOverTreeR!(fg, tree, N=N)

global x5 = KDE.getKDEMean(getKDE(getVariable(fg, :x5)))

@test abs(x5[1]) < 1.5
@test abs(x5[2]) < 1.5
@test abs(TU.wrapRad(x5[3])) < 0.4
@test abs(x5[4]) < 0.5
@test abs(x5[5]) < 0.5

global x10 = KDE.getKDEMean(getKDE(getVariable(fg, :x10)))

@test abs(x10[1]) < 2.75
@test abs(x10[2]) < 2.75
@test abs(TU.wrapRad(x10[3])) < 0.5
@test abs(x10[4]) < 0.5
@test abs(x10[5]) < 0.5

# pull the tail end out with position
global pp10 = DynPose2VelocityPrior(MvNormal([10.0;0;0], Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x10;], pp10)

batchSolve!(fg, N=N)
# run(`evince /tmp/caesar/bt.pdf`)

global x10 = KDE.getKDEMean(getKDE(getVariable(fg, :x10)))

@test 5.0 < x10[1]
@test abs(x10[2]) < 1.0
@test abs(TU.wrapRad(x10[3])) < 0.6
@test -0.1 < x10[4] < 1.0
@test abs(x10[5]) < 0.5

for sym in [Symbol("x$i") for i in 2:9]

global XX = KDE.getKDEMean(getKDE(getVariable(fg, sym)))

@show sym, round.(XX,digits=5)
@test -1.5 < XX[1] < 10.0
@test abs(XX[2]) < 1.0
@test abs(TU.wrapRad(XX[3])) < 1.3
@test -0.3 < XX[4] < 2.0
@test abs(XX[5]) < 0.5

end

end

# using RoMEPlotting
# plotLocalProduct(fg, :x10, dims=[1;2])
# drawPoses(fg)
# plotPose(fg, [:x9],levels=1);
#
# plotPose(fg, [:x1;:x2;:x3;:x4;:x5;:x6;:x7;:x8;:x9;:x10],levels=1);

# savejld(fg) # tempfg.jld

@testset "test many DynPose2 sideways velocity..." begin

global N = 100
global fg = initfg()

# add first pose locations
addVariable!(fg, :x0, DynPose2(ut=0))

# Prior factor as boundary condition
global pp0 = DynPose2VelocityPrior(MvNormal([0.0;0.0;pi/2], Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([0.0;0], Matrix(Diagonal([0.5; 0.5].^2))))
addFactor!(fg, [:x0;], pp0)

addVariable!(fg, :x1, DynPose2(ut=1000_000))

global pp0 = DynPose2VelocityPrior(MvNormal([1.0;0.0;pi/2], Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([0.0;0], Matrix(Diagonal([0.5; 0.5].^2))))
addFactor!(fg, [:x1;], pp0)

# conditional likelihood between Dynamic Point2
global dp2dp2 = VelPose2VelPose2(MvNormal([0.0;-1.0;0], Matrix(Diagonal([0.01;0.01;0.001].^2))),
                          MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x0;:x1], dp2dp2)

batchSolve!(fg,N=N)

# test for velocity in the body frame
global x0 = KDE.getKDEMean(getKDE(getVariable(fg, :x0)))

@test -0.4 < x0[1] < 2.0
@test abs(x0[2]) < 0.5
@test abs(x0[3] - pi/2) < 0.1
@test abs(x0[4]) < 0.4
@test -1.5 < x0[5] < -0.5

global x1 = KDE.getKDEMean(getKDE(getVariable(fg, :x1)))

@test -0.1 < x1[1] < 2.0
@test abs(x1[2]) < 0.5
@test abs(x1[3] - pi/2) < 0.1
@test abs(x1[4]) < 0.4
@test -1.5 < x1[5] < -0.5

end

# using RoMEPlotting
#
# drawPoses(fg)
#
# plotPose(fg, [:x0;:x1]);

#


text:
include( "bug_1.jl")

function bug2_w1()
    plt = plot!()
    goal = (x = 10, y = 10)
    path = [(x = 0, y = 0), (x = 1, y = 1), (x = 1, y = 5),
(x = 2, y = 5), (x = 2, y = 2), (x = 4, y = 4),
 (x = 3, y = 4), (x = 3, y = 12), (x = 3, y = 13),
  (x = 12, y = 13), (x = 13, y = 13),(x = 13, y = 5),
  (x = 12, y = 5),(x = 6, y = 5), (x = 6, y = 6), goal]
    dist_total = follower(path, goal)
    println("Total Distance is : ", dist_total)
end

function bug2_w2(w2)
    goal = (x = 35, y = 0)
    start=(x=0, y=0)
    path = [(x = 0, y = 0),(x = 4, y = 0), (x = 4, y = 1), (x = 5, y = 1),
    (x = 5, y = 0), (x = 9, y = 0),(x = 9, y = 5), (x = -5, y = 5), (x = -5, y = 0), (x = 0, y = 0)]# (x=10, y=0),(x=14,y=0), (x=14,y=1),

    plot!([start.x, 4], [ start.y, 0], legend=false, color=:black, arrow=:arrow, linewidth=4)
    dist = follower(path, goal)
    println("Total Distance is : ", dist)
end

function follower(path, goal)

    start = path[1];
    next_point = ()
    dist_total = 0
    for i = 1:size(path, 1) - 1

        if (i + 1 ) > size(path, 1)
            i = 0
        end
        next_point = path[i + 1];
        # println(start, next_point)
        plot!([start.x,  next_point.x], [  start.y, next_point.y], legend=false, arrow=:arrow, color=:black, linewidth=4)
        dist_total += round(sqrt((start.x - next_point.x)^2 + (start.y - next_point.y)^2))
        start = next_point;
    end
    # plot!([ goal.x, path[1].x], [ goal.y, path[1].y], legend=false, color=:blue, linewidth=1)
    return dist_total
end


text:
# import Base.Multimedia.display

struct RmdDisplay <: Display
end

const begin_plot_in_R = R"JuliaCall:::begin_plot"
const finish_plot_in_R = R"JuliaCall:::finish_plot"

function plot_display(display::RmdDisplay, x)
    path = rcopy(begin_plot_in_R())
    mkpath(dirname(path))
    Main.Plots.savefig(x, path)
    finish_plot_in_R()
end

const text_display_in_R = R"JuliaCall:::text_display"

function text_display(display::RmdDisplay, x)
    s = limitstringmime(MIME("text/plain"), x)
    text_display_in_R(s)
end

function display(display::RmdDisplay, x)
    try
        if isa(x, Main.Plots.Plot)
            plot_display(display, x)
            return
        end
    catch e;
    end
    try
        text_display(display, x)
        return
    catch e;
    end
    throw(MethodError)
end

const rmd_display = RmdDisplay()

max_stars_repo_name: KevCaz/JuliaCall

max_issues_repo_name: tinagui/JuliaCall

max_forks_repo_name: tinagui/JuliaCall


text:
module MySubDir1Tests

using Random
using XUnit
import ..MacroDefMod

@testcase "subdir1-test1" begin
    include("subdir1-include.jl")
end

MacroDefMod.@my_testcase "subdir1-test2" begin
    include("subdir1-include.jl")
end

end


text:
# start julia with: julia -p 4

using MAVs
using VehicleModels

mdl,n,c,r,params=initializeAutonomousControl(:c3);
pa=params[1];

function createParallelModels(n,c,pa)

#  http://docs.julialang.org/en/stable/manual/parallel-computing/#clustermanagers
  #addprocs(3);
  Xmin=[0.,180.,220.];
  #sm1=remotecall(initializeSimpleModel,2,(n,c,pa;x_min=Xmin[1]));
  #sm2=remotecall(initializeSimpleModel,3,(n,c,pa;x_min=Xmin[2]));
  #sm3=remotecall(initializeSimpleModel,4,(n,c,pa;x_min=Xmin[3]));
  sm1=@spawn initializeSimpleModel(n,c,pa;x_min=Xmin[1])
  sm2=@spawn initializeSimpleModel(n,c,pa;x_min=Xmin[2])
  sm3=@spawn initializeSimpleModel(n,c,pa;x_min=Xmin[3])

  sm=[sm1,sm2,sm3];
  return sm
end

sm=createParallelModels(n,c,pa);

#https://github.com/JuliaParallel/ClusterManagers.jl

#https://juliacomputing.com/domains/parallel-computing.html
#https://github.com/JuliaParallel


text:
# using PooledArrays
# factor = PooledArray
using CategoricalArrays

factor(args...) = CategoricalArray(args...) |> compress

factor_value(x::CategoricalValue) = levels(x)[x.ref]

export factor, factor_value, CategoricalArrays


text:
## TODO: pending https://github.com/JuliaLang/julia/issues/29420
# this one is suggested in the issue, but it looks like time_t and tm are two different things?
# const Ctime_t = Base.Libc.TmStruct

const Ctm = Base.Libc.TmStruct
const Ctime_t = UInt
const Cclock_t = UInt
const pthread_t = UInt


text:
using githubactions4julia
using Test

@testset "githubactions4julia.jl" begin
    # Write your tests here.
    @test plus(1,1)==2 #結果は成功
    @test plus(1,1)==3 #結果は失敗
    @test fail
end


text:
#=
    Functions in this file test functionality relating to objectives in MOI.

### Requires
    - optimize!

### Functionality currently tested
    - get/set ObjectiveSense
    - a constant in a affine objective
    - a blank objective

### Functionality not yet tested
    - Quadratic Objectives
    - Modifications
=#

"""
    max_sense(model::MOI.ModelLike, config::TestConfig)

Test setting objective sense to MAX_SENSE.
"""
function max_sense(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    @test MOI.supports(model, MOI.ObjectiveSense())
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
end
unittests["max_sense"] = max_sense

"""
    min_sense(model::MOI.ModelLike, config::TestConfig)

Test setting objective sense to MIN_SENSE.
"""
function min_sense(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    @test MOI.supports(model, MOI.ObjectiveSense())
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
end
unittests["min_sense"] = min_sense

"""
    feasibility_sense(model::MOI.ModelLike, config::TestConfig)

Test setting objective sense to FEASIBILITY_SENSE.
"""
function feasibility_sense(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    @test MOI.supports(model, MOI.ObjectiveSense())
    MOI.set(model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
    @test MOI.get(model, MOI.ObjectiveSense()) == MOI.FEASIBILITY_SENSE
end
unittests["feasibility_sense"] = feasibility_sense

"""
    get_objective_function(model::MOI.ModelLike, config::TestConfig)

Test get objective function.
"""
function get_objective_function(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    obj_attr = MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}()
    @test MOI.supports(model, obj_attr)
    MOIU.loadfromstring!(
        model,
        """
    variables: x
    minobjective: 2.0x + 1.0
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    expected_obj_fun =
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2.0, x)], 1.0)
    @test_throws InexactError begin
        MOI.get(model, MOI.ObjectiveFunction{MOI.SingleVariable}())
    end
    obj_fun = MOI.get(model, obj_attr)
    @test obj_fun ≈ expected_obj_fun
    quad_obj_attr =
        MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}()
    quad_obj_fun = MOI.get(model, quad_obj_attr)
    @test convert(MOI.ScalarAffineFunction{Float64}, quad_obj_fun) ≈
          expected_obj_fun
end
unittests["get_objective_function"] = get_objective_function

"""
    solve_constant_obj(model::MOI.ModelLike, config::TestConfig)

Test constant in linear objective,  if `config.solve=true` confirm that it
solves correctly, and if `config.duals=true`, check that the duals are computed
correctly.
"""
function solve_constant_obj(model::MOI.ModelLike, config::TestConfig)
    atol, rtol = config.atol, config.rtol
    MOI.empty!(model)
    @test MOI.is_empty(model)
    MOIU.loadfromstring!(
        model,
        """
    variables: x
    minobjective: 2.0x + 1.0
    c: x >= 1.0
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.get(
        model,
        MOI.ConstraintIndex{MOI.SingleVariable,MOI.GreaterThan{Float64}},
        "c",
    )
    # We test this after the creation of every `SingleVariable` constraint
    # to ensure a good coverage of corner cases.
    @test c.value == x.value
    return test_model_solution(
        model,
        config;
        objective_value = 3.0,
        variable_primal = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual = [(c, 2.0)],
    )
end
unittests["solve_constant_obj"] = solve_constant_obj

"""
    solve_blank_obj(model::MOI.ModelLike, config::TestConfig)

Test blank linear objective,  if `config.solve=true` confirm that it
solves correctly, and if `config.duals=true`, check that the duals are computed
correctly.
"""
function solve_blank_obj(model::MOI.ModelLike, config::TestConfig)
    atol, rtol = config.atol, config.rtol
    MOI.empty!(model)
    @test MOI.is_empty(model)
    MOIU.loadfromstring!(
        model,
        """
    variables: x
    minobjective: 0.0x + 0.0
    c: x >= 1.0
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.get(
        model,
        MOI.ConstraintIndex{MOI.SingleVariable,MOI.GreaterThan{Float64}},
        "c",
    )
    @test c.value == x.value
    test_model_solution(
        model,
        config;
        objective_value = 0.0,
        constraint_dual = [(c, 0.0)],
    )
    if config.solve
        # The objective is blank so any primal value ≥ 1 is correct
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
        @test MOI.get(model, MOI.VariablePrimal(), x) + atol + rtol ≥ 1.0
        @test MOI.get(model, MOI.ConstraintPrimal(), c) + atol + rtol ≥ 1.0
    end
end
unittests["solve_blank_obj"] = solve_blank_obj

"""
    solve_singlevariable_obj(model::MOI.ModelLike, config::TestConfig)

Test SingleVariable objective,  if `config.solve=true` confirm that it
solves correctly, and if `config.duals=true`, check that the duals are computed
correctly.
"""
function solve_singlevariable_obj(model::MOI.ModelLike, config::TestConfig)
    atol, rtol = config.atol, config.rtol
    MOI.empty!(model)
    @test MOI.is_empty(model)
    MOIU.loadfromstring!(
        model,
        """
    variables: x
    minobjective: x
    c: x >= 1.0
""",
    )
    x = MOI.get(model, MOI.VariableIndex, "x")
    c = MOI.get(
        model,
        MOI.ConstraintIndex{MOI.SingleVariable,MOI.GreaterThan{Float64}},
        "c",
    )
    @test c.value == x.value
    return test_model_solution(
        model,
        config;
        objective_value = 1.0,
        variable_primal = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual = [(c, 1.0)],
    )
end
unittests["solve_singlevariable_obj"] = solve_singlevariable_obj

"""
    solve_qp_edge_cases(model::MOI.ModelLike, config::TestConfig)

Test various edge cases relating to quadratic programs (i.e., with a quadratic
objective function).

If `config.solve=true` confirm that it solves correctly.
"""
function solve_qp_edge_cases(model::MOI.ModelLike, config::TestConfig)
    obj_attr = MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}()
    MOI.supports(model, obj_attr) || return
    MOI.empty!(model)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    vc1 = MOI.add_constraint(
        model,
        MOI.SingleVariable(x[1]),
        MOI.GreaterThan(1.0),
    )
    @test vc1.value == x[1].value
    vc2 = MOI.add_constraint(
        model,
        MOI.SingleVariable(x[2]),
        MOI.GreaterThan(2.0),
    )
    @test vc2.value == x[2].value

    @testset "Basic model" begin
        # min x^2 + y^2 | x>=1, y>=2
        MOI.set(
            model,
            obj_attr,
            MOI.ScalarQuadraticFunction(
                MOI.ScalarAffineTerm{Float64}[],  # affine terms
                MOI.ScalarQuadraticTerm.([2.0, 2.0], x, x),  # quad
                0.0,  # constant
            ),
        )
        test_model_solution(
            model,
            config;
            objective_value = 5.0,
            variable_primal = [(x[1], 1.0), (x[2], 2.0)],
        )
    end
    @testset "Duplicate linear terms" begin
        # min x + x + x^2 + y^2 | x>=1, y>=2
        MOI.set(
            model,
            obj_attr,
            MOI.ScalarQuadraticFunction(
                MOI.ScalarAffineTerm.([1.0, 1.0], [x[1], x[1]]),  # affine terms
                MOI.ScalarQuadraticTerm.([2.0, 2.0], x, x),  # quad
                0.0,  # constant
            ),
        )
        test_model_solution(
            model,
            config;
            objective_value = 7.0,
            variable_primal = [(x[1], 1.0), (x[2], 2.0)],
        )
    end
    @testset "Duplicate diagonal terms" begin
        # min x^2 + x^2 | x>=1, y>=2
        MOI.set(
            model,
            obj_attr,
            MOI.ScalarQuadraticFunction(
                MOI.ScalarAffineTerm{Float64}[],  # affine terms
                MOI.ScalarQuadraticTerm.(
                    [2.0, 2.0],
                    [x[1], x[1]],
                    [x[1], x[1]],
                ),  # quad
                0.0,  # constant
            ),
        )
        test_model_solution(
            model,
            config;
            objective_value = 2.0,
            variable_primal = [(x[1], 1.0)],
        )
    end
    @testset "Duplicate off-diagonal terms" begin
        # min x^2 + 0.25x*y + 0.25y*x + 0.5x*y + y^2 | x>=1, y>=2
        MOI.set(
            model,
            obj_attr,
            MOI.ScalarQuadraticFunction(
                MOI.ScalarAffineTerm{Float64}[],  # affine terms
                MOI.ScalarQuadraticTerm.(
                    [2.0, 0.25, 0.25, 0.5, 2.0],
                    [x[1], x[1], x[2], x[1], x[2]],
                    [x[1], x[2], x[1], x[2], x[2]],
                ),  # quad
                0.0,  # constant
            ),
        )
        test_model_solution(
            model,
            config;
            objective_value = 7.0,
            variable_primal = [(x[1], 1.0), (x[2], 2.0)],
        )
    end
end
unittests["solve_qp_edge_cases"] = solve_qp_edge_cases

"""
    solve_qp_zero_offdiag(model::MOI.ModelLike, config::TestConfig)

Test quadratic program with a zero off-diagonal term.

If `config.solve=true` confirm that it solves correctly.
"""
function solve_qp_zero_offdiag(model::MOI.ModelLike, config::TestConfig)
    obj_attr = MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}()
    MOI.supports(model, obj_attr) || return
    MOI.empty!(model)
    x = MOI.add_variables(model, 2)
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    vc1 = MOI.add_constraint(
        model,
        MOI.SingleVariable(x[1]),
        MOI.GreaterThan(1.0),
    )
    @test vc1.value == x[1].value
    vc2 = MOI.add_constraint(
        model,
        MOI.SingleVariable(x[2]),
        MOI.GreaterThan(2.0),
    )
    @test vc2.value == x[2].value
    MOI.set(
        model,
        obj_attr,
        MOI.ScalarQuadraticFunction(
            MOI.ScalarAffineTerm{Float64}[],  # affine terms
            MOI.ScalarQuadraticTerm.(
                [2.0, 0.0, 2.0],
                [x[1], x[1], x[2]],
                [x[1], x[2], x[2]],
            ),  # quad
            0.0,  # constant
        ),
    )
    test_model_solution(
        model,
        config;
        objective_value = 5.0,
        variable_primal = [(x[1], 1.0), (x[2], 2.0)],
    )
end
unittests["solve_qp_zero_offdiag"] = solve_qp_zero_offdiag

"""
    solve_duplicate_terms_obj(model::MOI.ModelLike, config::TestConfig)

Test duplicate terms in linear objective, if `config.solve=true` confirm that it
solves correctly.
"""
function solve_duplicate_terms_obj(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    x = MOI.add_variable(model)
    c = MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(1.0))
    @test c.value == x.value
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(
            [MOI.ScalarAffineTerm(2.0, x), MOI.ScalarAffineTerm(1.0, x)],
            0.0,
        ),
    )
    return test_model_solution(
        model,
        config;
        objective_value = 3.0,
        variable_primal = [(x, 1.0)],
        constraint_primal = [(c, 1.0)],
        constraint_dual = [(c, 3.0)],
    )
end
unittests["solve_duplicate_terms_obj"] = solve_duplicate_terms_obj


text:
export detrend!

function detrend!(A::AbstractArray)
    sx, sy, sxy, sx2 = 0.0, 0.0, 0.0, 0.0
    n = size(A, 1)
    for i in 1:n
        sx += i
        sy += A[i]
        sxy += i * A[i]
        sx2 += i * i
    end

    slope = (sx * sy - n * sxy) / (sx * sy - n * sx2)
    constant = (sy - slope * sx) / n

    for i in 1:n
        A[i] -= (constant + slope * i)
    end
end


text:
module GalewskyShallowWaterThetaMethod

using Gridap
using GridapGeosciences
using GridapPardiso
using SparseMatricesCSR

include("GalewskyInitialConditions.jl")

# Solves the Galewsky test case for the shallow water equations on a sphere
# of physical radius 6371220m. Involves a shear flow instability of a zonal
# jet triggered by an initial gravity wave.
# reference:
#   Galewsky, Scott and Polvani (2004) Tellus, 56A 429-440

order  = 1
degree = 4

n      = 48
dt     = 480.0
nstep  = Int(24*60^2*20/dt) # 20 days
T      = dt*nstep
θ      = 0.5

model = CubedSphereDiscreteModel(n; radius=rₑ)

linear_solver=PardisoSolver(GridapPardiso.MTYPE_REAL_NON_SYMMETRIC,
                            GridapPardiso.new_iparm(),
                            GridapPardiso.MSGLVL_VERBOSE,
                            GridapPardiso.new_pardiso_handle())

shallow_water_theta_method_full_newton_time_stepper(model, order, degree,
                                               h₀, u₀, f, topography, g, θ, T, nstep, dt/2;
                                               linear_solver=linear_solver,
                                               sparse_matrix_type=SparseMatrixCSR{1,Float64,Int},
                                               write_solution=true,
                                               write_solution_freq=45,
                                               write_diagnostics=true,
                                               write_diagnostics_freq=1,
                                               dump_diagnostics_on_screen=true)

end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences