text:
using Logging

logger = Logging.SimpleLogger(stdout, Logging.Info)
global_logger(logger)

include("../src/ExactODEReduction.jl")
using .ExactODEReduction: load_ODEs, construct_jacobians,
        many_invariant_subspaces, polynormalize, invariant_subspace_1,
        load_ODEs_recursive_if, normalizepath

using .ExactODEReduction.Nemo
using .ExactODEReduction.AbstractAlgebra

if !isdefined(Main, :cache)
    cache = Dict()
end

function run_tcs_cms()

    for (i, (filename, system)) in enumerate(load_ODEs_recursive_if("data/TCS_CMSB"))

        @info "$i-th, loaded a system $filename of size $(length(system))"

        As = construct_jacobians(system)

        if !haskey(cache, system)
            Vs = many_invariant_subspaces(As, backend_algorithm=invariant_subspace_1)
            cache[system] = Vs
        end

        println("System $system")
        println("Found $(length(cache[system])) invariant subspaces")

        for v in cache[system]
            println("-----------------------")
            polyform = polynormalize(v, parent(first(first(system))))
            println(length(polyform))
            println(join(string.(polyform), "\n"))
            println("-----------------------")
        end

        # break
    end

end

run_tcs_cms()


text:
# Hinge loss function

export HingeLoss

"""
**Hinge loss**

  HingeLoss(y, μ=1.0)

Returns the function
```math
f(x) = μ⋅∑_i \\max\\{0, 1 - y_i ⋅ x_i\\},
```
where `y` is an array and `μ` is a positive parameter.
"""
HingeLoss(y::T, mu::R=1.0) where {T <: AbstractArray, R <: Real} = Postcompose(PrecomposeDiagonal(SumPositive(), -y, 1.0), mu)

ext: jl


text:
@doc raw"""
    LatentFactorMOKernel(g::AbstractVector{<:Kernel}, e::MOKernel, A::AbstractMatrix)

Kernel associated with the semiparametric latent factor model.

# Definition

For inputs ``x, x'`` and output dimensions ``p_x, p_{x'}'``, the kernel is defined as[^STJ]
```math
k\big((x, p_x), (x, p_{x'})\big) = \sum^{Q}_{q=1} A_{p_xq}g_q(x, x')A_{p_{x'}q}
                                   + e\big((x, p_x), (x', p_{x'})\big),
```
where ``g_1, \ldots, g_Q`` are ``Q`` kernels, one for each latent process, ``e`` is a
multi-output kernel for ``m`` outputs, and ``A`` is a matrix of weights for the kernels of
size ``m \times Q``.

[^STJ]: M. Seeger, Y. Teh, & M. I. Jordan (2005). [Semiparametric Latent Factor Models](https://infoscience.epfl.ch/record/161465/files/slfm-long.pdf).
"""
struct LatentFactorMOKernel{Tg,Te<:MOKernel,TA<:AbstractMatrix} <: MOKernel
    g::Tg
    e::Te
    A::TA
    function LatentFactorMOKernel(g, e::MOKernel, A::AbstractMatrix)
        all(gi isa Kernel for gi in g) || error("`g` should be an collection of kernels")
        length(g) == size(A, 2) ||
            error("Size of `A` not compatible with the given array of kernels `g`")
        return new{typeof(g),typeof(e),typeof(A)}(g, e, A)
    end
end

function (κ::LatentFactorMOKernel)((x, px)::Tuple{Any,Int}, (y, py)::Tuple{Any,Int})
    cov_f = sum(κ.A[px, q] * κ.g[q](x, y) * κ.A[py, q] for q in 1:length(κ.g))
    return cov_f + κ.e((x, px), (y, py))
end

function kernelmatrix(k::LatentFactorMOKernel, x::MOInput, y::MOInput)
    x.out_dim == y.out_dim || error("`x` and `y` should have the same output dimension")
    x.out_dim == size(k.A, 1) ||
        error("Kernel not compatible with the given multi-output inputs")

    # Weights matrix ((out_dim x out_dim) x length(k.g))
    W = [col * col' for col in eachcol(k.A)]

    # Latent kernel matrix ((N x N) x length(k.g))
    H = [gi.(x.x, permutedims(y.x)) for gi in k.g]

    # Weighted latent kernel matrix ((N*out_dim) x (N*out_dim))
    W_H = sum(kron(Wi, Hi) for (Wi, Hi) in zip(W, H))

    return W_H .+ kernelmatrix(k.e, x, y)
end

function Base.show(io::IO, k::LatentFactorMOKernel)
    return print(io, "Semi-parametric Latent Factor Multi-Output Kernel")
end

function Base.show(io::IO, ::MIME"text/plain", k::LatentFactorMOKernel)
    print(io, "Semi-parametric Latent Factor Multi-Output Kernel\n\tgᵢ: ")
    join(io, k.g, "\n\t\t")
    print(io, "\n\teᵢ: ")
    return join(io, k.e, "\n\t\t")
end


text:
export @surf_expr, @surf_toplevel, @surf_type
function surf_type(s::Symbol)
    Surf.TVar(s)
end

function surf_type(ex)
    @match ex begin
        QuoteNode(a::Symbol) => Surf.TSym(a)
        :(NewType($(a::Symbol))) => Surf.TNew(a)
        :(implicit[$t]) => Surf.TImplicit(surf_type(t))
        :(Fn[$a, $b]) => Surf.TArrow(surf_type(a), surf_type(b))
        :($a.?($label)) => Surf.TQuery(string(label), surf_type(a))
        :($a.$(b::Symbol)) => Surf.TApp(Surf.TVar(:field), Surf.TSym(b))
        :($f[$a]) => Surf.TApp(surf_type(f), surf_type(a))
        Expr(:tuple, args...) =>
            Surf.TTuple(Surf.TyExpr[surf_type(a) for a in args])
        Expr(:->, a, Expr(:block, ::LineNumberNode, b)) =>
            Surf.TArrow(surf_type(a), surf_type(b))
        :($t where {$(tvars...)}) =>
            Surf.TForall(Symbol[tvars...], surf_type(t))
    end
end

function surf_expr(n::Symbol)
    Surf.EVar(n)
end

function surf_expr(i::Union{Integer, AbstractString, AbstractChar, AbstractFloat})
    Surf.EVal(i)
end

function surf_expr(exp::Expr)
    @match exp begin
    :(@julia $(::LineNumberNode) $jlex) => Surf.EExt(jlex)
    :(@extern $(::LineNumberNode) $jlex) => Surf.EExt(jlex)
    :($a.?($label)) => Surf.EQuery(string(label), surf_expr(a))
    :($a.$(b::Symbol)) => Surf.EApp(Surf.EVar(:op_Getter), Surf.EField(surf_expr(a), b))
    :($f($(args...))) =>
        length(args) === 1 ?
        Surf.EApp(surf_expr(f), surf_expr(args[1])) :
        Surf.EApp(
            surf_expr(f),
            Surf.ETup(Surf.Expr[surf_expr(arg) for arg in args]))
    Expr(:tuple, elts...) =>
        Surf.ETup(Surf.Expr[surf_expr(arg) for arg in elts])

    Expr(:->, s::Symbol, exp && Expr(:block, ln::LineNumberNode, _...)) =>
        Surf.ELoc(
            ln,
            Surf.EFun(s, surf_expr(exp)))

    Expr(:if, cond, arm1, arm2) =>
        Surf.EInt(surf_expr(cond), surf_expr(arm1), surf_expr(arm2))

    Expr(:let, :($_ = $_) && bind, body) =>
        surf_expr(Expr(:let, Expr(:block, bind), body))

    Expr(:let, Expr(:block, decls...), body) =>
        Surf.ELet(surf_decls(decls), surf_expr(body))

    Expr(:block, suite...) =>
        begin
            if isempty(suite)
                Surf.ETup(Surf.Expr[])
            else
                ret = surf_expr(suite[end])
                foldr((@view suite[1:end-1]), init=ret) do new, last
                    @match new begin
                        ln::LineNumberNode =>
                            @match last begin
                                Surf.ELet(vec, _) =>
                                    begin
                                        pushfirst!(vec, DLoc(ln))
                                        last
                                    end
                                _ => Surf.ELoc(ln, last)
                            end
                        :($_ :: $_) || :($_ = $_) || :($_.?($_)) =>
                            @match last begin
                                Surf.ELet(vec, _) =>
                                    begin
                                        pushfirst!(vec, surf_decl(new))
                                        last
                                    end
                                _ => Surf.ELet([surf_decl(new)], last)
                            end
                        _ => Surf.ELet([Surf.DBind(:_, surf_expr(new))], last)
                    end
                end
            end
        end
    end
end

surf_decl(decl::Any) = @match decl begin
    :($(x::Symbol).?($label)) => Surf.DQuery(string(label), x)
    :($(x::Symbol) :: $t) => Surf.DAnn(x, surf_type(t))
    :($(x::Symbol) = $v) => Surf.DBind(x, surf_expr(v))
    ln::LineNumberNode => Surf.DLoc(ln)
end

function surf_decls(decls::AbstractVector)
    Surf.Decl[surf_decl(decl) for decl in decls]
end

macro surf_toplevel(ex)
    @match ex begin
        Expr(:block, stmts...) => surf_decls(stmts)
        _ => error("syntax error")
    end
end

function surf_module(ex)
    @match ex begin
        Expr(:module, true, name, stmts...) => (name, surf_decls(stmts))
        _ => error("syntax error")
    end
end

macro surf_expr(ex)
    surf_expr(ex)
end

macro surf_type(ex)
    surf_type(ex)
end


text:
@testset "meta-pkg" begin
    rr = RidgeRegressor()
    d = MLJBase.info_dict(rr)
    @test d[:package_name] == "MLJLinearModels"
    @test d[:package_url] == "https://github.com/alan-turing-institute/MLJLinearModels.jl"
    @test d[:package_license] == "MIT"
    @test d[:is_pure_julia] == true
    @test d[:is_wrapper] == false
end

@testset "meta-reg" begin
    lr = LinearRegressor()
    d = MLJBase.info_dict(lr)
    @test d[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test d[:target_scitype] == AbstractVector{MLJBase.Continuous}
    @test d[:supports_weights] == false
    @test !isempty(d[:docstring])
    @test d[:load_path] == "MLJLinearModels.LinearRegressor"
end

@testset "meta-clf" begin
    lr = LogisticClassifier()
    d = MLJBase.info_dict(lr)
    @test d[:input_scitype] == MLJBase.Table(MLJBase.Continuous)
    @test d[:target_scitype] == AbstractVector{<:MLJBase.Finite}
    @test d[:supports_weights] == false
    @test !isempty(d[:docstring])
    @test d[:load_path] == "MLJLinearModels.LogisticClassifier"
end


text:
import sbp
using sbp.Grid:grid_xp, grid_xm, grid_2d_x, grid_2d_y
using Test

@testset "Grid properties" begin
        n = 10
        x, h = grid_xp(n)
        @test isapprox(x[2] - x[1], x[end] - x[end-1])
        @test isapprox(x[1], 0)
        @test isapprox(x[end], 1.0)

        x, h = grid_xm(n)
        @test isapprox(x[2] - x[1], x[end] - x[end-1])
        @test isapprox(x[1], 0)
        @test isapprox(x[end], 1.0)

        x, h = grid_xp(n, pad=true)
        @test length(x) == n + 1
        @test x[end] == 0
        @test isapprox(x[end - 1], 1.0)

end

@testset "2D grid vector" begin
        nx = 10
        ny = 12
        x1, h = grid_xp(nx)
        y1, h = grid_xp(ny)
        x = grid_2d_x(x1, ny)
        @test size(x) == (nx * ny,)
        @test isapprox(x[1:nx], x1[1] * ones(nx))

        y = grid_2d_y(y1, nx)
        @test size(y) == (nx * ny,)
        @test isapprox(y[1:ny], y1)

end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
function cholin!(kb::Matrix{Float64})
  n = size(kb, 1)
  iw = size(kb, 2) - 1
  for i in 1:n
    x = 0.0
    for j in 1:iw
      x += kb[i, j]^2.0
    end
    println([i, iw, iw+1, x, kb[i, iw+1], kb[i, iw+1] - x])
    kb[i, iw+1] = sqrt(kb[i, iw+1] - x)
    for k in 1:iw
      x = 0.0
      if i + k <= n
        if k != iw
          for l in iw - k:-1:1
            x += kb[i+k, l] * kb[i, l+k]
          end
        end
        ia = i + k
        ib = iw - k + 1
        kb[ia, ib] = (kb[ia, ib] - x) / kb[i, iw+1]
      end
    end
  end
end


text:
module TestPredatorPrey
include("preamble_plots.jl")

using Bifurcations: resolved_points
using Bifurcations.Codim1
using Bifurcations.Codim2: NormalizingAS, BackReferencingAS
using Bifurcations.Examples.PredatorPrey

# Bifurcation points calculated in some version of Bifurcations.jl
CODIM1_KNOWN_POINTS = Dict(
    # PredatorPrey parameter (tuple) => points
    (0.0, 3, 5, 3) => [
        (Codim1.PointTypes.saddle_node,
         [0.0, 0.0, 0.6]),
        (Codim1.PointTypes.saddle_node,
         [1/3, 0, 0.821905296656178]),
        (Codim1.PointTypes.saddle_node,
         [0.41116155950979133, 0, 0.8329293250960671]),
        (Codim1.PointTypes.hopf,
         [1/3, 0.36576153597718997, 0.6715938474788321]),
    ],
)

@testset "smoke PredatorPrey" begin

    codim1_solver = init(PredatorPrey.prob, nominal_angle_rad=0.01)
    solve!(codim1_solver)

    @testset "known codim1 points" begin
        resolved = resolved_points(codim1_solver)
        known = CODIM1_KNOWN_POINTS[codim1_solver.prob.p.de_prob.p]
        @test length(resolved) == length(known)
        for (actual, desired) in zip(resolved, known)
            type_desired, u_desired = desired
            @test actual.point_type == type_desired
            @test actual.u ≈ u_desired  rtol=1e-4
        end
    end

    @testset "codim-2 ($(nameof(ASType))) from $(point.u0[[1, end]])" for
            ASType in [NormalizingAS, BackReferencingAS],
            point in special_intervals(codim1_solver)

        codim2_prob = BifurcationProblem(
            point,
            codim1_solver,
            (@lens _[3]),
            (0.0, 10.0),
            augmented_system = ASType(),
        )
        # TODO: Implement detection of Bogdanov-Takens bifurcation and
        # stop manually setting max_branches=0.
        codim2_solver = init(codim2_prob; max_branches=0)
        @test_nothrow solve!(codim2_solver)

        @testset "show" begin
            smoke_test_solver_show(codim2_solver)
        end

        @testset "plot" begin
            smoke_test_solver_plot(codim2_solver)
        end
    end
end

end  # module


text:
module TestAqua

using Aqua
using InitialValues
using BangBang

Aqua.test_all(InitialValues)

Aqua.test_ambiguities(
    [InitialValues, BangBang, Base];
)

end  # module


text:
using Countries
import Countries: symbol, country, name, code, capital, continent, isdeveloping, region, subregion

using Test

countries = country.([:US,:PH,:HK,:SG])
names = ["United States of America","Philippines","China, Hong Kong Special Administrative Region","Singapore"]
codes = [840,608,344,702]
capitals = ["Washington","Manila","Hong Kong","Singapore"]
continents = ["NA","AS","AS","AS"]
isdevelopings = [false,true,true,true]
regions = ["Americas","Asia","Asia","Asia"]
subregions = ["Northern America","South-eastern Asia","Eastern Asia","South-eastern Asia"]

@testset "Basic countries" begin
    for (c,nm,cd,cp,ct,dv,rg,sr) in zip(countries,names,codes,capitals,continents,isdevelopings,regions,subregions)
        @test name(c) == nm
        @test code(c) == cd
        @test capital(c) == cp
        @test continent(c) == ct
        @test isdeveloping(c) == dv
        @test region(c) == rg
        @test subregion(c) == sr
    end
end


text:
# todo:
#

# requires the following R packages installed:
#   alpaca
#   formula.tools

#
# function feglm(df::AbstractDataFrame, m::Model; kwargs...)
#     feglm(df, m.f; m.dict..., kwargs...)
# end

function feglm(df::AbstractDataFrame, f::FormulaTerm,
    family::UnivariateDistribution;
    fe::Union{Symbol, Expr, Nothing} = nothing,
    vcov::Union{Symbol, Expr, Nothing} = :(simple()),
    start::Union{Vector{T}, Nothing} = nothing,
    maxiter::Integer = 10000, limit::Integer = 10,
    trace::Integer = 0,
    convtol::Real = 1.0e-06,
    devtol::Real = 1.0e-08, steptol::Real = 1.0e-08,
    centertol::Real = 1.0e-05, rhotol::Real = 1.0e-04,
    droppc::Bool = true
   ) where T<:Real

    # Check fixed effect formula and construct string
    feformula = fe
    fe_str = ""
    has_absorb = feformula != nothing
    if has_absorb
        # check depth 1 symbols in original formula are all CategoricalVector
        if isa(feformula, Symbol)
            x = feformula
            !isa(df[!,x], CategoricalVector) && error("$x should be CategoricalVector")
            fe_str = String(fe)
        elseif feformula.args[1] == :+
            x = feformula.args
            for i in 2:length(x)
                isa(x[i], Symbol) && !isa(df[!,x[i]], CategoricalVector) && error("$(x[i]) should be CategoricalVector")
                fe_str = "$fe_str + $(String(x[i]))"
            end
        end
    end

    # Start your engines
    R"library(alpaca)"
    R"library(formula.tools)"

    # Check family
    if family == Binomial()
        R"fam <- binomial()"
    elseif family == Normal()
        R"fam <- gaussian()"
    elseif family == Gamma()
        R"fam <- Gamma()"
    elseif family == InverseGaussian()
        R"fam <- inverse.gaussian()"
    elseif family == Poisson()
        R"fam <- poisson()"
    else
        error("Family not supported.")
    end

    if start != nothing
        @rput start
    else
        R"start <- NULL"
    end

    cluster = false
    if vcov == :(simple())
        R"type <- c(\"hessian\")"
        R"str_vcov_formula <- NULL"
    elseif vcov == :robust
        R"type <- c(\"sandwich\")"
        R"str_vcov_formula <- NULL"
    elseif typeof(vcov) == Expr
        (vcov.args[1] == :cluster) || error("Invalid vcov argument.")
        cluster = true
        if isa(vcov.args[2],Symbol)
            clusterstring = string(vcov.args[2])
            fstring = " ~ $(clusterstring)"
            @rput fstring
            @rput clusterstring
            R"vcov_f <- as.formula(fstring)"
            # vcov_f = @formula(String(vcov.args[2]))
            # str_vcov_formula = [String(vcov.args[2])]
        else
            (vcov.args[2].args[1] == :+) || error("Cluster formula specification invalid.")
            clusterstring = string(vcov.args[2])
            fstring = " ~ $(clusterstring)"
            @rput fstring
            @rput clusterstring
            R"vcov_f <- as.formula(fstring)"
            # vcov_f = @formula( ~ String.(vcov.args[2].args[2:end]))
            # str_vcov_formula = String.(vcov.args[2].args[2:end])
        end
        R"type <- c(\"clustered\")"
        # @show str_vcov_formula
        # @rput str_vcov_formula
    else
        error("Invalid vcov argument.")
    end

    @rput trace
    @rput limit

    # move objects to R
    r_df = robject(df)
    r_f = robject(f)
    @rput f
    @rput fe_str

    # Julia's formulae don't have the FE added.
    # Do this, then remove all brackets (alpaca doesn't like them)
    R"fstr <- as.character(update(f, paste(\"~ . |\",fe_str ) ) )"
    if cluster
        R"fstr <- as.character(update(f, paste(\"~ . |\",fe_str, \" |\" ,clusterstring ) ) )"
    end
    R"fstr <- gsub(\"(\", \"\", fstr, fixed=TRUE)"
    R"fstr <- gsub(\")\", \"\", fstr, fixed=TRUE)"
    R"f <- as.formula(fstr)"

    # confirm what we're running
    R"f_string <- as.character(f)"
    f_str = @rget f_string
    println("Running alpaca with formula: $f_str")

    R"ctrl <- feglmControl(dev.tol = $devtol,
       center.tol = $centertol, rho.tol = $rhotol,
       conv.tol = $convtol,
       iter.max = $maxiter, limit = limit,
       trace = trace,
       drop.pc = $droppc)"

    # R"ctrl <- feglm.control(step.tol = $steptol, dev.tol = $devtol,
    #       pseudo.tol = $pseudotol, rho.tol = $rhotol,
    #       iter.max = $maxiter, trace = $trace,
    #       drop.pc = $droppc)"

    R"result <- feglm(formula =  f , data = $r_df,
            family = fam , beta.start = start,
            control = ctrl)"

    if cluster
        output = R"sum <- summary(result, type = type, cluster = vcov_f)"
        R"vcov <- vcov(result, type = type, cluster = vcov_f )"
    else
        output = R"sum <- summary(result, type = type)"
        R"vcov <- vcov(result, type = type )"
    end
    println(output)

    R"coefnames <- names(coef(result))"
    R"yname <- all.vars(result[[\"formula\"]])[1]"
    R"nobs <- result[[\"nobs\"]]"
    R"lvlsk <- result[[\"lvls.k\"]]"
    @rget coefnames
    @rget yname
    @rget nobs
    @rget lvlsk

    println("All done.")
    @rget result
    @rget vcov

    # convert scalars to arrays, if they are scalars
    coef = result[:coefficients]
    coef = isa(coef, Vector) ? coef : [coef]
    # this is a bit hacky...
    vcov = isa(vcov, Array) ? vcov : vcov .+ zeros(Float64,1,1)
    coefnames = isa(coefnames, Vector) ? coefnames : [coefnames]

    # @show lvlsk
    # @show nobs
    # @show coef
    # @show vcov

    rr = RegressionResult(coef,
        vcov,
        coefnames,
        Symbol(yname),
        f,
        nobs[4], sum(lvlsk)
    )

   return rr

end

# Interface with strings. Hacky, and a temporary solution.
# function feglm(df::AbstractDataFrame, formula::String, family::String;
#     vcov::Union{Symbol, Expr, Nothing} = :(simple()))
#
#     # load alpaca in R
#     R"library(alpaca)"
#
#     println("Formula")
#     r_df = robject(df)
#     @rput formula
#     R"f <- as.formula(formula)"
#
#     # run it
#     println("Running feglm")
#     R"result <- feglm(formula =  f , data = $r_df, family = binomial())"
#
#     println("Vcov")
#     if vcov == :(simple())
#         R"type <- c(\"hessian\")"
#         R"str_vcov_formula <- NULL"
#     elseif vcov == :robust
#         R"type <- c(\"sandwich\")"
#         R"str_vcov_formula <- NULL"
#     elseif typeof(vcov) == Expr
#         (vcov.args[1] == :cluster) || error("Invalid vcov argument.")
#         if isa(vcov.args[2],Symbol)
#             str_vcov_formula = [String(vcov.args[2])]
#         else
#             (vcov.args[2].args[1] == :+) || error("Cluster formula specification invalid.")
#             str_vcov_formula = String.(vcov.args[2].args[2:end])
#         end
#         R"type <- c(\"clustered\")"
#         @show str_vcov_formula
#         @rput str_vcov_formula
#     else
#         error("Invalid vcov argument.")
#     end
#
#     println("summary")
#     output = R"sum <- summary(result, type = type,
#            cluster.vars = str_vcov_formula)"
#     println(output)
#
#     R"vcov <- vcov(result, type = type, cluster.vars = str_vcov_formula )"
#     R"coefnames <- names(coef(result))"
#     R"yname <- all.vars(result[[\"formula\"]])[1]"
#     R"nobs <- result[[\"nobs\"]]"
#     R"lvlsk <- result[[\"lvls.k\"]]"
#     @rget coefnames
#     @rget yname
#     @rget nobs
#     @rget lvlsk
#
#     @rget result
#     @rget vcov
#
#     coef = result[:coefficients]
#     coef = isa(coef, Vector) ? coef : [coef]
#     # this is a bit hacky...
#     vcov = isa(vcov, Array) ? vcov : vcov .+ zeros(Float64,1,1)
#     coefnames = isa(coefnames, Vector) ? coefnames : [coefnames]
#
#     rr = RegressionResult(coef,
#         vcov,
#         coefnames,
#         Symbol(yname),
#         @formula(y ~ x),
#         nobs[4], sum(lvlsk)
#     )
#
#    return rr
#
# end

macro vcov(args...)
    Expr(:call, :vcov_helper, (esc(Base.Meta.quot(a)) for a in args)...)
end

function vcov_helper(args...)

    (args[1].head === :call && args[1].args[1] === :(~)) || throw("First argument of @model should be a formula")
    f = @eval(@formula($(args[1].args[2]) ~ $(args[1].args[3])))
    dict = Dict{Symbol, Any}()
    for i in 2:length(args)
        isa(args[i], Expr) &&  args[i].head== :(=) || throw("All arguments of @model, except the first one, should be keyboard arguments")
        dict[args[i].args[1]] = args[i].args[2]
    end
    Model(f, dict)
end


text:
using Levenshtein

function randString(len)
    return string([Char(round(Int,rand() * 255 + 1)) for _ in 1:len]...)
end

function run_huge_string(l1, l2)
    # setup long string against short one to expose element copy performance
    s = randString(l1)
    t = randString(l2)
    ins, del, sub = rand() * 1000, rand() * 1000, rand() * 1000

    # run a long string against a not so long one
    println("\nRunning $l1-char String against $l2-char String")
    @time levenshtein(s, t, ins, del, sub)
end

function run_multiple_times(l1, l2, runs)
    # setup
    s = randString(l1)
    t = randString(l2)
    ins, del, sub = rand() * 1000, rand() * 1000, rand() * 1000
    costs = Array{Float64}(2, max(length(s), length(t)) + 1)

    # run without reusing cost matrix
    println("\nRunning $runs runs without reusing cost matrix for $l1-char String vs $l2-char String")
    @time for i in 1:runs
        levenshtein(s, t, ins, del, sub)
    end

    # run without allocating cost matrix everytime
    println("\nRunning $runs runs with resuable cost matrix for $l1-char String vs $l2-char String")
    @time for i in 1:runs
        levenshtein!(s, t, ins, del, sub, costs)
    end
end

run_huge_string(10^7, 10^2)
run_huge_string(10^6, 10^3)
run_huge_string(10^5, 10^4)
run_multiple_times(10^3, 10^3, 10^3)
run_multiple_times(10^2, 10^2, 10^5)
run_multiple_times(10^1, 10^1, 10^7)


text:
# ============================================================
# --------------------GraffCloudOnPose--------------------
# ============================================================
# import Arena.Amphitheatre: visualize!

struct GraffCloudOnPose <: AbstractAmphitheatre
	robotId::String
	sessionId::String
	config::GraffConfig

	# pclouds::Dict{Symbol, PointCloud} #TODO do we want to store this here as well?
	isdrawn::Dict{Symbol, Bool}
	dcam::CameraModel
	colmap::ColorGradient
end

GraffCloudOnPose(config::GraffConfig, dcam::CameraModel, colmap::ColorGradient = cgrad(:bgy,:colorcet)) =
			GraffCloudOnPose(config.robotId, config.sessionId, config, Dict{Symbol, Bool}(), dcam, colmap)

"""
    $(SIGNATURES)
Basic visualizer object visualize! function.
"""
function visualize!(vis::Visualizer, pcop::GraffCloudOnPose)::Nothing

	# only read and update from graff once
	isdrawn = pcop.isdrawn

	if length(isdrawn) == 0

		robotId = pcop.robotId
		sessionId = pcop.sessionId
		config = pcop.config

		sessionDataEntries = getSessionDataEntries(config.robotId, config.sessionId)

		for nodeKey = keys(sessionDataEntries)

			dataEntries = sessionDataEntries[nodeKey]

			!("rawdepth" in [de.id for de in dataEntries]) && continue

			elem = GraffSDK.getData(robotId, sessionId, dataEntries[1].neoNodeId, "rawdepth")

			if elem == nothing
				# @info "no rawdepth for $vsym"
				continue
			end

			depthImage = JSON2.read(elem.data)
			# es = String(copy(elem.data[1:findfirst(elem.data .== 0x00)]));

			kTc = (SE3([0,0,0], Quaternion(Float64.(depthImage.kQi[1:4]))))
			trans = Translation([0,0,0])∘LinearMap(Quat(kTc.R.R))

			w = depthImage.width
			h = depthImage.height

			depthim = reshape(reinterpret(UInt16,base64decode(depthImage.data[1].base64)), (w, h))'[:,:]

			pointcloud = cloudFromDepthImageClampZ(depthim, pcop.dcam, trans, maxrange=3f0, clampz = [-0.5f0,0.5f0], colmap=pcop.colmap)
			# pointcloud = cloudFromDepthImage(depthim, pcop.dcam, trans, maxrange=3f0)

			nodesym = Symbol(depthImage.node_name)

			if !get(isdrawn, nodesym, false)
				haskey(isdrawn, nodesym) ? isdrawn[nodesym] = true : push!(isdrawn, nodesym=>true)
				visPointCloudOnPose!(vis[robotId][sessionId][:poses][nodesym][:pc], pointcloud)
			end
		end
	end
    return nothing
end


text:
function display_test()
	@testset "DisplayLH" begin
		# @testset "sprintf" begin
		# 	fmtStr = "{:.1f} "
		# 	z = [1.11 2.22]
		# 	x = DisplayLH.sprintf(fmtStr, z)
		# 	@test x == "1.1 2.2 "
		# 	DisplayLH.printf(fmtStr, z)
		# end

		sV = ["aaaa aaaa", "bbb bbb bbb", "cccccc cccccc cccccc", "dddddddddd"];
		CommonLH.show_string_vector(sV, 25);
		lineV = CommonLH.string_vector_to_lines(sV, 25);
		@test length(lineV) == 3
		@test startswith(lineV[1], sV[1])
		@test startswith(lineV[2], sV[3])
		@test startswith(lineV[3], sV[4])
	end
end

function multi_io_test()
	@testset "MultIO" begin
		fPath = joinpath(test_dir(), "multi_io_test.txt");
		open(fPath, "w") do io
			mio = MultiIO([stdout, io]);
			print(mio, "Line 1\n");
			println(mio, "Line 2");
		end
		lineV = readlines(fPath);
		for j = 1 : 2
			@test lineV[j] == "Line $j"
		end
	end
end

@testset "Display" begin
	display_test()
	multi_io_test()
end

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


text:
################################################################################
## Copyright 2016. Juri Hinz and Jeremy Yee.
## Test the prescribed policy for a set of sample paths
################################################################################

function test_policy(css::CSS, path::Path, pathPolicy::PathPolicy, position::Int64)

    ## Get sizes
    pathnum = path.pathnum
    pnum = css.pnum
    tnum = css.tnum
    anum = css.anum
    snum = css.snum

    ## Container
    result = zeros(Array{Float64}(pathnum))
    state = Array{Float64}(1, snum)

    ## Test run of policy
    if css.fullcontrol # positions are fully controllable
        for i = 1:pathnum
            p = position
            for t = 1:(tnum - 1)
                state = path.sample[t, i, :]
                action = pathPolicy.policy[i, p, t]
                result[i] += dot(css.Reward(t, state, p, action), state)
                p = css.control[p, action]
            end
            result[i] += dot(css.Scrap(path.sample[tnum, i, :], p),
                             path.sample[tnum, i, :])
        end
    else
        ## Cumulative probabilities
        cum_prob = cumsum(css.control, 3)
        ## Test run of policy
        for i = 1:pathnum
            p = position
            for t = 1:(tnum - 1)
                state = path.sample[t, i, :]
                action = pathPolicy.policy[i, p, t]
                result[i] += dot(css.Reward(t, state, p, action), state)
                ## Get next position
                rand_unif = rand()
                for pp = 1:pnum
                    if (rand_unif <= cum_prob[p, action, pp])
                        p = pp
                        break
                    end
                end
            end
            result[i] += dot(css.Scrap(path.sample[tnum, i, :], p),
                             path.sample[tnum, i, :])
        end
    end

    return result

end


text:
@testset "GGP read" begin
	# Simple data loading
	ggp1 = readggp(joinpath(dirname(@__DIR__),"test","input","readggp_data.ggp"),nanval=9999.999)
	@test ggp1[:datetime][1] == DateTime(2017,08,26)
	@test ggp1[:datetime][end] == DateTime(2017,09,03,4,0,0)
	for i in [9,10,12]
		@test isnan(ggp1[:column2][i])
	end
	@test sum(ggp1[:column1]) ≈ 5.
	@test sum(filter(!isnan,ggp1[:column2])) ≈ 20.

	# Apply offset while loading data
	ggp2 = readggp(joinpath(dirname(@__DIR__),"test","input","readggp_data.ggp"),
		nanval=9999.999,offset=true)
	@test ggp2[:datetime][1] == DateTime(2017,08,26)
	@test ggp2[:datetime][end] == DateTime(2017,09,03,4,0,0)
	for i in [9,10,12]
		@test isnan(ggp2[:column2][i])
	end
	@test sum(ggp2[:column1]) ≈ 55.
	@test sum(filter(!isnan,ggp2[:column2])) ≈ 100.

	# Read header only
	ggphead = readggp(joinpath(dirname(@__DIR__),"test","input","readggp_data.ggp"),
		what="header")
	@test size(ggphead) == (2,)

	# Read EOSTloading (Pre)ETERNA file
	ggp3 = readggp(joinpath(dirname(@__DIR__),"test","input","eost_data.rot"))
	@test names(ggp3) == [:datetime,:column1,:column2,:column3,:column4,:column5]
	ggp3m = readdlm(joinpath(dirname(@__DIR__),"test","input","eost_data.rot"),skipstart=14)
	ggp3m = ggp3m[1:end-1,:];
	@test ggp3[:datetime] == DateTime.(string.(ggp3m[:,1]),"yyyymmdd")
	for i in 3:size(ggp3m,2)
		@test sum(ggp3[Symbol("column",i-2)]) ≈ sum(ggp3m[:,i])
	end

	# Read blockinfor
	blockinfo = readggp(joinpath(dirname(@__DIR__),"test","input","eost_data.rot"),what="blocks")
	@test blockinfo[:datetime][1] == DateTime(string(ggp3m[1,1]),"yyyymmdd")
	@test blockinfo[:startline][1] == 15
	@test blockinfo[:stopline][1] == size(ggp3m,1)+15-1 # header - 999999 line
	@test size(blockinfo) == (1,size(ggp3m,2)+1)
end

@testset "GGP write" begin
	dataout = DataFrame(pres=collect(1000.12345:1:1011.123456),
	       				datetime=collect(DateTime(2010,1,1):Dates.Hour(1):DateTime(2010,1,1,11)));
	# try writtin without optional parameters
	writeggp(dataout,joinpath(dirname(@__DIR__),"test","output","ggp_data_one.dat"));
	# read and test
	d = readdlm(joinpath(dirname(@__DIR__),"test","output","ggp_data_one.dat"),skipstart=3)
	for i = 1:12
		@test d[i,1] == 20100101
		@test d[i,3] ≈ round(dataout[:pres][i]*100)/100
		@test d[i,2] ≈ (i-1)*1e+4
	end
	@test d[end,1] == 99999999

	# Try full output
	dataout[:grav] =collect(900.123456:-3:(900.123456-11*3));
	# set units using the same keywords as in the dataframe
	units = Dict(:grav=>"V",:pres=>"hPa")
	# Set header. All entries with the excetpion of "freetext" will be formatted
	header = Dict("Filename"=>"file.data",
	              "Station"=>"Wettzell",
	              "Instrument"=>"iGrav",
	              "N. Latitude (deg)"=>49.14354,
	              "E. Longitude (deg)"=>12.87866,
	              "Elevation MSL (m)"=>613.7+1.05,
	              "Author"=>"Name (name@gfz-potsdam.de)",
				  "freetext"=>"Line without formatting\nMaximum 80 characters per line");
	block = Dict("start"=>[DateTime(2010,1,1,09,0,0)],
				"offset"=>[10.1234 20.1234], # one row per block. 2 values per row = number channels
	            "header"=>["iGrav006" 1.0 1.0 0.0 3])
	fileout = joinpath(dirname(@__DIR__),"test","output","ggp_data_two.dat");
	writeggp(dataout,fileout,
				units=units,header=header,decimal=[1,3],
				blockinfo=block,channels=[:grav,:pres])
	# test
	open(fileout,"r") do fid
		row = ""
		for i in 1:10
			row = readline(fid);
		end
		@test row == "yyyymmdd hhmmss grav(V) pres(hPa)";
	end
	b = readdlm(fileout,skipstart=12)
	@test b[1,1] == 77777777
	@test b[2,1] == 20100101
	@test b[end-1,1] == 20100101
	@test b[end-1,2] == 110000
	@test b[3,2] == 010000
	@test b[3,3] ≈ round(dataout[:grav][2]*10)/10
	@test b[end-2,4] ≈ round(dataout[:pres][end-1]*1000)/1000
	@test b[end-4,2] ≈ round(block["offset"][1]*10)/10
	@test b[end-4,3] ≈ round(block["offset"][2]*1000)/1000
	@test b[end-5,1] == block["header"][1]
	@test b[end-5,2:end] == block["header"][2:end]
end

@testset "GGP data2blocks" begin
	datawrite = DataFrame(pres=collect(1000.:1:1011.),
					grav=collect(900.:-3:(900. -11*3)),
					datetime=collect(DateTime(2010,1,1):Dates.Hour(1):DateTime(2010,1,1,11)));
	datawrite[:grav][[3,6]] .= NaN;
	datawrite[:pres][7] = NaN;
	dataout,block = ggpdata2blocks(datawrite;channels=[])
	@test filter(!isnan,datawrite[:grav]+datawrite[:pres]) ≈ dataout[:grav]+dataout[:pres]
	@test block["start"] == datawrite[:datetime][[3,6]]
	@test block["offset"] == zeros(Float64,2,2)

	#
	dataout2,block2 = ggpdata2blocks(dataout;channels=[:pres])
	@test dataout2[:datetime] == dataout[:datetime]
	@test dataout2[:pres] == dataout[:pres]
	@test !haskey(dataout2,:grav)
	@test isempty(block2)
end


text:
using Test
using TestExtras
using Random
using TensorKit
using Combinatorics
using TensorKit: ProductSector, fusiontensor
using TensorOperations
TensorOperations.disable_cache() # avoids memory overflow during CI?
using Base.Iterators: take, product
# using SUNRepresentations: SUNIrrep
# const SU3Irrep = SUNIrrep{3}
import LinearAlgebra

include("newsectors.jl")
using .NewSectors

const TK = TensorKit

Random.seed!(1234)

smallset(::Type{I}) where {I<:Sector} = take(values(I), 5)
function smallset(::Type{ProductSector{Tuple{I1,I2}}}) where {I1,I2}
    iter = product(smallset(I1), smallset(I2))
    s = collect(i ⊠ j for (i,j) in iter if dim(i)*dim(j) <= 6)
    return length(s) > 6 ? rand(s, 6) : s
end
function smallset(::Type{ProductSector{Tuple{I1,I2,I3}}}) where {I1,I2,I3}
    iter = product(smallset(I1), smallset(I2), smallset(I3))
    s = collect(i ⊠ j ⊠ k for (i,j,k) in iter if dim(i)*dim(j)*dim(k) <= 6)
    return length(s) > 6 ? rand(s, 6) : s
end
function randsector(::Type{I}) where {I<:Sector}
    s = collect(smallset(I))
    a = rand(s)
    while a == one(a) # don't use trivial label
        a = rand(s)
    end
    return a
end
function hasfusiontensor(I::Type{<:Sector})
    try
        fusiontensor(one(I), one(I), one(I))
        return true
    catch e
        if e isa MethodError
            return false
        else
            rethrow(e)
        end
    end
end

sectorlist = (Z2Irrep, Z3Irrep, Z4Irrep, U1Irrep, CU1Irrep, SU2Irrep, NewSU2Irrep,# SU3Irrep,
              FibonacciAnyon, IsingAnyon, Z3Irrep ⊠ Z4Irrep, U1Irrep ⊠ SU2Irrep, SU2Irrep ⊠
              SU2Irrep, NewSU2Irrep ⊠ NewSU2Irrep, NewSU2Irrep ⊠ SU2Irrep, SU2Irrep ⊠
              NewSU2Irrep, Z2Irrep ⊠ FibonacciAnyon ⊠ FibonacciAnyon)

Ti = time()
include("sectors.jl")
include("fusiontrees.jl")
include("spaces.jl")
include("tensors.jl")
Tf = time()
printstyled("Finished all tests in ",
            string(round((Tf-Ti)/60; sigdigits=3)),
            " minutes."; bold = true, color = Base.info_color())
println()


text:
# The "iris" dataset contains 150 measurements on 3 species
# of Iris flowers, namely I. setosa, I. versicolor and I. virginica.
# We want to know how many measurements there are for each species.
# Our input is an array named iris with for each measurement the name of its species.
# The dataset can be found in the file iris.txt

# 1- extracting the species from the dataset in the array:
iris = []
fname = "iris.txt"
open(fname) do file
    for line in eachline(file)
        # println(line)
        # the species name starts from column 17:
        species = line[17:end]
        # removes the trailing \r\n:
        species = strip(species)
        push!(iris, species)
    end
end
# println(length(iris)) # 150
# show(iris) # Any["I. setosa","I. setosa","I. s... ]

# 2- counting the species:
# make a dictionary with the species and count their frequencies:
species_freq = Dict{ASCIIString, Int64}()
for species in iris
    haskey(species_freq, species) ?
        species_freq[species] += 1 :
        species_freq[species] = 1
end

# 3- sort the the keys and print out the frequencies:
println("Species : frequency \n")
species = sort!(collect(keys(species_freq)))
for sp in species
    println("$sp : $(species_freq[sp])")
end
# Species : frequency

# I. setosa : 50
# I. versicolor : 50
# I. virginica : 50

max_stars_repo_path: Section 1/1_6_counting_flowers.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 1/1_6_counting_flowers.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 1/1_6_counting_flowers.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
abstract type AbstractIRL end

include("linear_irl.jl")


text:
println("pid: $(getpid())")
println("Threads: $(Threads.nthreads())")
flush(stdout)

using Cytof5
using Random
using Distributions
using BSON

include("simulatedata.jl")

# Parse arguments
if length(ARGS) == 0
  RESULTS_DIR = "results/test-test/"
  REPFAMDISTSCALE = 10
  KMCMC = 5
  Z_idx = 3
  SEED = 0
  SEED_MCMC = -1
else
  RESULTS_DIR = ARGS[1]
  REPFAMDISTSCALE = parse(Float64, ARGS[2])
  KMCMC = parse(Int, ARGS[3])
  Z_idx = parse(Int, ARGS[4])
  SEED = parse(Int, ARGS[5])
  SEED_MCMC = parse(Int, ARGS[6])  # -1 (no seed) for everything before sim 5-12
end
mkpath(RESULTS_DIR)
USE_REPULSIVE = REPFAMDISTSCALE > 0
# USE_REPULSIVE = true  # use the repulsive joint updates even when scale=0
Z = Zs[Z_idx]

# Random.seed!(0)  # before test-sims-5
Random.seed!(SEED)  # after test-sims-5

println("CONFIG:")
println("    - RESULTS_DIR: $(RESULTS_DIR)")
println("    - REPFAMDISTSCALE: $(REPFAMDISTSCALE)")
println("    - KMCMC: $(KMCMC)")
println("    - Z_idx: $(Z_idx)")
println("    - USE_REPULSIVE: $(USE_REPULSIVE)")
println("    - SEED: $(SEED)")
println("    - SEED_MCMC: $(SEED_MCMC)")
flush(stdout)

function sim_z_generator(phi)::Function
  # larger phi -> higher similarity -> higher penalty
  # smaller phi -> lower similarity -> lower penalty
  # phi = 0 -> regular IBP
  return (z1::Vector{Bool}, z2::Vector{Bool}) -> begin
    exp(-sum(abs.(z1 - z2)) / phi)
  end
end

sim_z = sim_z_generator(REPFAMDISTSCALE)

function init_state_const_data(simdat; K, L)
  deltaz_prior = TruncatedNormal(1.0, 0.1, 0.0, Inf)
  # deltaz_prior = TruncatedNormal(1.0, 0.3, 0.5, Inf)
  # deltaz_prior = TruncatedNormal(1.0, 1e-10, 0.5, Inf)
  # deltaz_prior = TruncatedNormal(0.0, 1.0, 0.75, Inf)
  # deltaz_prior = TruncatedNormal(0.0, 1.0, 0.9, Inf)
  # deltaz_prior = TruncatedNormal(0.0, 1.0, 1.0, Inf)
  # deltaz_prior = TruncatedNormal(0.0, 0.3, 0.5, Inf)
  # deltaz_prior = TruncatedNormal(1.0, 0.0001, 0.0, Inf)
  # deltaz_prior = TruncatedNormal(1.0, 0.0001, 0.5, Inf)
  d = Cytof5.Model.Data(simdat[:y])
  c = Cytof5.Model.defaultConstants(d, K, L,
                                    tau0=1.0, tau1=1.0,
                                    # sig2_prior=InverseGamma(3, 2),
                                    sig2_prior=InverseGamma(11, 5),  # TODO
                                    delta0_prior=deltaz_prior,
                                    delta1_prior=deltaz_prior,
                                    alpha_prior=Gamma(0.1, 10.0),
                                    yQuantiles=[.0, .25, .5],
                                    pBounds=[.05, .8, .05],
                                    probFlip_Z=1.0,
                                    similarity_Z=sim_z)
  # s = Cytof5.Model.genInitialState(c, d)
  s = Cytof5.Model.smartInit(c, d)
  t = Cytof5.Model.Tuners(d.y, c.K)
  X = Cytof5.Model.eye(Float64, d.I)

  cfs = Cytof5.Model.ConstantsFS(c)
  dfs = Cytof5.Model.DataFS(d, X)
  sfs = Cytof5.Model.StateFS{Float64}(s, dfs)
  tfs = Cytof5.Model.TunersFS(t, s, X)

  return Dict(:dfs => dfs, :cfs => cfs, :sfs => sfs, :tfs => tfs,
              :simdat => simdat)
end

# Simulate data
# NOTE: Used before test-sims-5
# @time simdat = simulatedata1(Z=Z, seed=0, propmissingscale=.6, sortLambda=true);

# NOTE: Used on and after test-sims-5
@time simdat = simulatedata1(Z=Z,
                             # N=[300, 300],  # for all but test-sims-5-8
                             # N =[500, 500],  # for test-sims-5-8
                             N=[1000, 1000],  # test-sims-6-1
                             W=Matrix(hcat([[.7, 0, .1, .1, .1],
                                            [.4, .1, .3, .1, .1]]...)'),
                             sig2=[.5, .5],
                             seed=SEED, propmissingscale=.6, sortLambda=true);

#= Sanity check
using PyPlot
const plt = PyPlot.plt
J = size(simdat[:Z], 1)
plt.imshow(simdat[:y][2], aspect="auto", cmap=plot_yz.blue2red.cm(7),
           vmin=VLIM[1], vmax=VLIM[2])
plt.xticks(pyrange(J), 1:J)
plt.colorbar()
=#

# Parameters to monitor
monitor1 = [:theta__Z, :theta__v, :theta__alpha,
            :omega, :r, :theta__lam, :W_star, :theta__eta,
            :theta__W, :theta__delta, :theta__sig2]
monitor2 = [:theta__y_imputed, :theta__gam]

# MCMC Specs
nsamps_to_thin(nsamps::Int, nmcmc::Int) = max(1, div(nmcmc, nsamps))
NSAMPS = 2000  # Number of samples
THIN_SAMPS = 2  # Factor to thin the primary parameters
MCMC_ITER = NSAMPS * THIN_SAMPS  # Number of MCMC iterations
# LPML / DIC are computed based on `MCMC_ITER` samples
NBURN = 10000  # burn-in time

# Configurations: priors, initial state, data, etc.
config = init_state_const_data(simdat, K=KMCMC, L=Dict(0 => 2, 1 => 2))

# Print constants
println("N: $(config[:dfs].data.N)")
println("J: $(config[:dfs].data.J)")
Cytof5.Model.printConstants(config[:cfs])
flush(stdout)

# Fit model
@time out = Cytof5.Model.fit_fs!(config[:sfs], config[:cfs], config[:dfs],
                                 tuners=config[:tfs],
                                 nmcmc=MCMC_ITER,
                                 nburn=NBURN,
                                 thins=[THIN_SAMPS,
                                        nsamps_to_thin(10, MCMC_ITER)],
                                 monitors=[monitor1, monitor2],
                                 computedden=true,
                                 thin_dden=nsamps_to_thin(200, MCMC_ITER),
                                 printFreq=10, time_updates=false,
                                 computeDIC=true, computeLPML=true,
                                 use_repulsive=USE_REPULSIVE,
                                 Z_thin=1,
                                 flushOutput=true,
                                 seed=SEED_MCMC)

# Dump output
BSON.bson("$(RESULTS_DIR)/output.bson", out)
BSON.bson("$(RESULTS_DIR)/simdat.bson", Dict(:simdat => config[:simdat]))

# TODO: send to s3?

println("Completed!")


text:
"""
    Kissler_VL_model.jl

Written by Carl Whitfield, University of Manchester, 2021

Only to be included by definitions.jl
"""

#===================== VL trajectory params  =====================#

#====> Kissler et al 2020 <====#
const beta_peak = 0.7
const alpha_peak = 3.5*beta_peak                      #Peak time from LOD distribution params
const beta_decay = 0.3
const alpha_decay_symp = 10.5*beta_decay              #Decay time to LOD distribution params
const alpha_decay_asymp = 6.7*beta_decay
const peakVL_mean = 7.533                             #Peak viral load distribution params
const peakVL_sd = 1.164
const VL_LOD_Kissler = log10(250) - 0.93733/3.60971   #Limit of Detection in Kissler study
const V_0 = 0.5255                  #from Ferretti (initial VL) (if using V0_model_opt)
const onset_frac = 0.38             #frac of peak time added before peak (if using PVT_model_opt)

#Distributions used to generate params
const PVLdist = truncated(Normal(peakVL_mean,peakVL_sd),VL_LOD_Kissler,Inf)
const PTdist = Gamma(alpha_peak,1.0/beta_peak)
const DTsympdist = Gamma(alpha_decay_symp,1.0/beta_decay)
const DTasympdist = Gamma(alpha_decay_asymp,1.0/beta_decay)

"""
    generate_VL_params_Kissler(Asymp::Bool)

Generate parameters for VL model.

## Arguments
`Asymp` = Whether individual is asymptomatic

## Returns:
`Float` = peak viral load (log10 copies/ml)
`Float` = peak viral load time (days since infection)
`Float` = growth rate (1/days)
`Float` = decay rate (1/days)
"""
function generate_VL_params_Kissler(Asymp::Bool)
    PVL = rand(PVLdist)  #peak viral load
    PT = rand(PTdist)    #time from LOD to peak

    r = log(10)*(PVL - VL_LOD_Kissler)/PT

    if onset_opt == V0_model_opt
        OT = log(10)*(VL_LOD_Kissler - V_0)/r    #time from infection to LOD
    elseif onset_opt == PVT_model_opt
        OT = onset_frac*PT
    end

    if Asymp
        DT = rand(DTasympdist)
    else
        DT = rand(DTsympdist)
    end

    tp = OT+PT
    d = log(10)*(PVL - VL_LOD_Kissler)/DT

    return PVL, tp, r, d
end


text:
module Functors

export @functor, @flexiblefunctor, fmap, fmapstructure, fcollect

include("functor.jl")
include("base.jl")

end # module


text:
@testset "1368.minimum-cost-to-make-at-least-one-valid-path-in-a-grid.jl" begin
    @test min_cost1368(hcat([[1, 1, 3], [3, 2, 2], [1, 1, 4]]...)) == 0
    @test min_cost1368(hcat([[1, 1, 1, 1], [2, 2, 2, 2], [1, 1, 1, 1], [2, 2, 2, 2]]...)) ==
          3
    @test min_cost1368(hcat([[1, 2], [4, 3]]...)) == 1
    @test min_cost1368(hcat([[2, 2, 2], [2, 2, 2]]...)) == 3
    @test min_cost1368(hcat([[4]]...)) == 0
end


text:
module VortexLattice

using LinearAlgebra
using StaticArrays
using Rotations

## Panel geometry and math tools
#==========================================================================================#

import ..AeroMDAO: Panel3D, panel_area, panel_coords, midpoint, panel_normal, transform, p1, p2, p3, p4, average_chord, average_width, dynamic_pressure, aerodynamic_coefficients, force_coefficient, moment_coefficient, rate_coefficient, weighted_vector, reshape_array, cartesian_to_freestream, freestream_to_cartesian

## Horseshoe setup
#==========================================================================================#

include("horseshoes.jl")
include("finite_core.jl")

# export Horseshoe, horseshoe_line, horseshoe_point, bound_leg, bound_leg_center, bound_leg_vector, r1, r2, points, make_horseshoes

## Reference frames
#==========================================================================================#

include("reference_frames.jl")

# export body_to_stability_axes, stability_to_body_axes, body_to_wind_axes, wind_to_body_axes, reflect_xz

## Influence matrix and solution of system
#==========================================================================================#

include("influences.jl")

# export solve_system

"""
    make_horseshoes(horseshoe_panels)

Get the bound legs and collocation points of the horseshoes defined by horseshoe `Panel3D`s.
"""
make_horseshoes(horseshoe_panels) =	@. horseshoe_line(horseshoe_panels), horseshoe_point(horseshoe_panels)

"""
    solve_system(horseshoes, collocation_points, normals, U, Ω)

Evaluate and return the vortex strengths ``\\Gamma``s given `Horseshoes`, their associated normal vectors (not necessarily the same as the panels' normals), the speed ``U`` and rotation vector ``\\Omega``.
"""
function solve_system(horseshoes, normals, U, Ω)
    V = map(r -> U + Ω × r, collocation_point.(horseshoes))
    AIC  = influence_matrix(horseshoes, collocation_point.(horseshoes), normals, -normalize(U), false)
    boco = boundary_condition(V, normals)
    Γs 	 = AIC \ boco
end

## Force evaluations
#==========================================================================================#

include("nearfield.jl")

# export nearfield_dynamics, nearfield_drag

include("farfield.jl")

# export farfield_dynamics

## Post-processing
#==========================================================================================#

# export case_dynamics, evaluate_coefficients

"""
    case_dynamics(Γ_comp, hs_comp, Γs, horseshoes, U, α, β, Ω, ρ, r_ref)

Compute the nearfield and farfield forces over a given component with associated vortex strengths `Γ_comp` and horseshoes `hs_comp`, using the vortex strengths ``\\Gamma``s and `horseshoes` of the entire configuration, the speed ``U``, angles of attack ``\\alpha`` and sideslip ``\\beta``, the rotation vector ``\\Omega``, the freestream density ``\\rho`` and the reference location for moments ``r_\\text{ref}``.
"""
function case_dynamics(Γ_comp, hs_comp, Γs, horseshoes, U, α, β, Ω, ρ, r_ref)
    # Compute near-field dynamics
    surface_forces, surface_moments = nearfield_dynamics(Γ_comp, hs_comp, Γs, horseshoes, U, Ω, ρ, r_ref)

    # Compute farfield dynamics
    trefftz_force = trefftz_forces(Γ_comp, hs_comp, norm(U), α, β, ρ)

    surface_forces, surface_moments, trefftz_force
end

# In case of only one component
case_dynamics(Γs, horseshoes, U, α, β, Ω, ρ, r_ref) = case_dynamics(Γs, horseshoes, Γs, horseshoes, U, α, β, Ω, ρ, r_ref)

function evaluate_coefficients(forces, moments, trefftz_force, U, α, β, Ω, ρ, S, c, b)
    V = norm(U)
    q = dynamic_pressure(ρ, V)

    # Computing summed coefficients
    force, moment = sum(forces), sum(moments)

    # Transform near-field dynamics to wind axes
    trans_force	 = body_to_wind_axes(force, α, β)
    trans_force  = [ nearfield_drag(force, U); trans_force[2:end] ]
    trans_moment = body_to_wind_axes(stability_flip(moment), α, β)

    # Compute coefficients
    nearfield_coeffs = aerodynamic_coefficients(trans_force, trans_moment, Ω, V, S, b, c, ρ)
    farfield_coeffs	 = force_coefficient(trefftz_force, q, S)

    # Non-dimensional panel coefficients
    CFs = force_coefficient.(forces, q, S)
    CMs = moment_coefficient.(moments, q, S, b, c)

    nearfield_coeffs, farfield_coeffs, CFs, CMs
end

# Streamlines
include("streamlines.jl")

# export streamlines

# EXPERIMENTAL: SYSTEM
include("system.jl")

# export VLMState, VLMSystem, compute_influence_matrix!, compute_boundary_condition!, compute_collocation_points!, compute_horseshoes!, compute_normals!, solve_system!, compute_nearfield_forces!, compute_farfield_forces!, solve_case!, set_horseshoe_panels!, set_camber_panels!, compute_wake_properties!, normals, horseshoes, collocation_points, AIC, RHS, circulations, surface_force_coefficients, surface_moment_coefficients, total_force_coefficients, farfield_force_coefficients

end

max_forks_repo_name: GodotMisogi/AeroMDAO


text:
export smooth_gaps, smooth_gaps!, zero_out, zero_out!

"""
  smooth_gaps(S)

Remove and taper around gaps in `SeisData` `S`.

Use `ungap!(S,tap=true,m=false)` before applying this function to fill gaps
with `nan` values.

# Arguments
- `S::SeisData`: 3-component `SeisData` for detection.

# Keywords
- `pad::Int`: Number of sample points to taper before and after gap.
"""
function smooth_gaps!(S::SeisData;pad::Real=100)

    # check pad input
    pad = max(2,pad)
    if !isa(pad,Int)
        pad = convert(Int,ceil(pad))
    end
    for ii = 1:S.n
        # find all nan / gaps
        ind = findall(isnan.(S.x[ii]))
        if length(ind) == 0
            continue
        end

        # set all nan equal to zero
        out = unique(ind)
        T = eltype(S.x[ii])
        S.x[ii][out] .= T(0)
        padarr = T(0.):T(pad-1)

        # right tapering of gaps
        rightind = findall(diff(ind) .!= 1)
        for jj in ind[rightind]
            maxind = min(jj+pad-1,length(S[ii].x))
            if length(jj:maxind) == length(padarr)
                S.x[ii][jj:maxind] .*= cos.(T(pi) ./ T(2) .+ T(pi) ./ T(2) .* padarr ./ pad).^2
            else
                S.x[ii][jj:maxind] .= T(0)
            end
        end

        # left tapering of gaps
        prepend!(ind,1)
        leftind = findall(diff(ind) .> 1)
        for jj in ind[leftind]
            minind = max(1,jj-pad+1)
            if length(minind:jj) == length(padarr)
                S.x[ii][minind:jj] .*= cos.(T(pi) ./ T(2) .* padarr ./ pad).^2
            else
                S.x[ii][minind:jj] .= T(0)
            end
        end
    end
    return nothing
end
smooth_gaps(S::SeisData;pad::Int=100) = (U = deepcopy(S);smooth_gaps!(U,pad=pad);return U)

"""
  zero_out!(S)

Set data in SeisData `S` == 0 where `S.x .> thresh`
"""
function zero_out!(S::SeisData;thresh::AbstractFloat=1e-7)
    for ii = 1:S.n
        T = eltype(S.x[ii])
        ind = findall(abs.(S.x[ii]) .< thresh)
        S.x[ii][ind] .= T(0)
    end
    return nothing
end
zero_out(S::SeisData; thresh::AbstractFloat=1e-7) = (U = deepcopy(S);zero_out!(U,thresh=thresh);return U)


text:
using VirtualArrays
if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end
using FactCheck

include("virtual_cat.jl")
include("VirtualArrays.jl")

FactCheck.exitstatus()


text:
#
# client.jl --
#
# Implement XPA client methods.
#
#------------------------------------------------------------------------------
#
# This file is part of XPA.jl released under the MIT "expat" license.
# Copyright (C) 2016-2019, Éric Thiébaut (https://github.com/emmt/XPA.jl).
#

"""

`XPA.TEMPORARY` can be specified wherever an `XPA.Client` instance is expected
to use a non-persistent XPA connection.

"""
const TEMPORARY = Client(C_NULL)

"""
```julia
XPA.Client()
```

yields a persistent XPA client handle which can be used for calls to
[`XPA.set`](@ref) and [`XPA.get`](@ref) methods.  Persistence means that a
connection to an XPA server is not closed when one of the above calls is
completed but will be re-used on successive calls.  Using `XPA.Client()`
therefore saves the time it takes to connect to a server, which could be
significant with slow connections or if there will be a large number of
exchanges with a given access point.

See also [`XPA.set`](@ref), [`XPA.get`](@ref), [`XPA.list`](@ref) and
[`XPA.find`](@ref).

"""
function Client()
    # The argument of XPAOpen is currently ignored (it is reserved for future
    # use).
    ptr = ccall((:XPAOpen, libxpa), Ptr{Cvoid}, (Ptr{Cvoid},), C_NULL)
    ptr != C_NULL || error("failed to create a persistent XPA connection")
    return finalizer(close, Client(ptr))
end

function Base.close(xpa::Client)
    if (ptr = xpa.ptr) != C_NULL
        xpa.ptr = C_NULL # avoid closing more than once!
        ccall((:XPAClose, libxpa), Cvoid, (Ptr{Cvoid},), ptr)
    end
    return nothing
end

"""
```julia
XPA.list(xpa=XPA.TEMPORARY)
```

yields a list of available XPA access points.  Optional argument `xpa` is a
persistent XPA client connection; if omitted, a temporary client connection
will be created.  The result is a vector of [`XPA.AccessPoint`](@ref)
instances.

See also [`XPA.Client`](@ref) and [`XPA.find`](@ref).

"""
function list(xpa::Client = TEMPORARY)
    lst = AccessPoint[]
    for str in split(chomp(get(String, xpa, "xpans")), r"\n|\r\n?";
                     keepempty=false)
        arr = split(str; keepempty=false)
        if length(arr) != 5
            @warn "expecting 5 fields per access point (\"$str\")"
            continue
        end
        access = UInt(0)
        for c in arr[3]
            if c == 'g'
                access |= GET
            elseif c == 's'
                access |= SET
            elseif c == 'i'
                access |= INFO
            else
                @warn "unexpected access string (\"$(arr[3])\")"
                continue
            end
        end
        push!(lst, AccessPoint(arr[1], arr[2], arr[4], arr[5], access))
    end
    return lst
end

"""
```julia
XPA.find([xpa=XPA.TEMPORARY,] ident) -> apt
```

yields the accesspoint of the first XPA server matching `ident` or `nothing` if
none is found.  If a match is found, the result `apt` is an instance of
`XPA.AccessPoint` and has the following members:

```julia
apt.class   # class of the access point (String)
apt.name    # name of the access point
apt.addr    # socket access method (host:port for inet,
apt.user    # user name of access point owner
apt.access  # allowed access
```

all members are `String`s but the last one, `access`, which is an `UInt`.

Argument `ident` may be a regular expression or a string of the
form `CLASS:NAME` where `CLASS` and `CLASS` are matched against the server
class and name respectively (they may be `"*"` to match any).

Keyword `user` may be used to specify the user name of the owner of the server
process, for instance `ENV["user"]` to match your servers.  The default
is `user="*"` which matches any user.

Keyword `throwerrors` may be set true (it is false by default) to automatically
throw an exception if no match is found (instead of returning `nothing`).

See also [`XPA.Client`](@ref), [`XPA.address`](@ref) and [`XPA.list`](@ref).

"""
find(ident::Union{AbstractString,Regex}; kwds...) =
    find(TEMPORARY, ident; kwds...)

function find(xpa::Client,
              ident::AbstractString;
              user::AbstractString = "*",
              throwerrors::Bool = false)::Union{AccessPoint,Nothing}
    i = findfirst(isequal(':'), ident)
    if i === nothing
        # allow any class
        class = "*"
        name = ident
    else
        class = ident[1:i-1]
        name = ident[i+1:end]
    end
    anyuser = (user == "*")
    anyclass = (class == "*")
    anyname = (name == "*")
    lst = list(xpa)
    for j in eachindex(lst)
        if ((anyuser || lst[j].user == user) &&
            (anyclass || lst[j].class == class) &&
            (anyname || lst[j].name == name))
            return lst[j]
        end
    end
    throwerrors && error(_noserversmatch(ident))
    return nothing
end

function find(xpa::Client,
              ident::Regex;
              user::AbstractString = "*",
              throwerrors::Bool = false)::Union{AccessPoint,Nothing}
    anyuser = (user == "*")
    lst = list(xpa)
    for j in eachindex(lst)
        if ((anyuser || lst[j].user == user) &&
            occursin(ident, lst[j].class*":"*lst[j].name))
            return lst[j]
        end
    end
    throwerrors && error(_noserversmatch(ident))
    return nothing
end

@noinline _noserversmatch(ident::AbstractString) =
    "no XPA servers match pattern \"$(ident)\""

@noinline _noserversmatch(ident::Regex) =
    "no XPA servers match regular expression \"$(ident.pattern)\""

"""
```julia
XPA.address(apt) -> addr
```

yields the address of XPA accesspoint `apt` which can be: an instance of
`XPA.AccessPoint`, a string with a valid XPA server address or a server
`class:name` identifier.  In the latter case, [`XPA.find`](@ref) is called to
find a matching server which is much longer.

"""
address(apt::XPA.AccessPoint) =
    apt.addr

function address(apt::AbstractString)
    i = findfirst(isequal(':'), apt)
    if (i === nothing ||
        tryparse(UInt, apt[1:i-1],  base = 16) === nothing ||
        tryparse(UInt, apt[i+1:end],  base = 10) === nothing)
        return address(XPA.find(apt; throwerrors = true))
    end
    return apt
end

"""
```julia
XPA.get([T, [dims,]] [xpa,] apt, args...)
```

retrieves data from one or more XPA access points identified by `apt` (a
template name, a `host:port` string or the name of a Unix socket file) with
arguments `args...` (automatically converted into a single string where the
arguments are separated by a single space).  Optional argument `xpa` specifies
an XPA handle (created by [`XPA.Client`](@ref)) for faster connections.  The
returned value depends on the optional arguments `T` and `dims`.

If neither `T` nor `dims` are specified, an instance of [`XPA.Reply`](@ref) is
returned with all the answer(s) from the XPA server(s).  The following keywords
are available:

* Keyword `nmax` specifies the maximum number of answers, `nmax=1` by default.
  Specify `nmax=-1` to use the maximum number of XPA hosts.

* Keyword `throwerrors` specifies whether to check for errors.  If this keyword
  is set true, an exception is thrown for the first error message encountered
  in the list of answers.  By default, `throwerrors` is false.

* Keyword `mode` specifies options in the form `"key1=value1,key2=value2"`.

If `T` and, possibly, `dims` are specified, a single answer and no errors are
expected (as if `nmax=1` and `throwerrors=true`) and the data part of the
answer is converted according to `T` which must be a type and `dims` which is
an optional list of dimensions:

* If only `T` is specified, it can be `String` to return a string interpreting
  the data as ASCII characters or a type like `Vector{S}` to return the largest
  vector of elements of type `S` that can be extracted from the returned data.

* If both `T` and `dims` are specified, `T` can be a type like `Array{S}` or
  `Array{S,N}` and `dims` a list of `N` dimensions to retrieve the data as an
  array of type `Array{S,N}`.

See also [`XPA.Client`](@ref), [`XPA.get_data`](@ref), [`XPA.set`](@ref) and
[`XPA.verify`](@ref).

"""
function get(xpa::Client,
             apt::AbstractString,
             cmd::AbstractString;
             mode::AbstractString = "",
             nmax::Integer = 1,
             throwerrors::Bool = false)
    return _get(xpa, apt, cmd, mode, _nmax(nmax), throwerrors)
end

function get(xpa::Client,
             apt::AbstractString,
             args::Union{AbstractString,Real}...;
             kwds...)
    return get(xpa, apt, _join(args); kwds...)
end

get(apt::AbstractString, args::Union{AbstractString,Real}...; kwds...) =
    get(TEMPORARY, apt, _join(args); kwds...)

_get1(args...; kwds...) = get(args...; nmax = 1, throwerrors = true, kwds...)

function get(::Type{Vector{T}},
             args...; kwds...) :: Vector{T} where {T}
    get_data(Vector{T}, _get1(args...; kwds...))
end

function get(::Type{Vector{T}}, dim::Integer,
             args...; kwds...) :: Vector{T} where {T}
    get_data(Vector{T}, dim, _get1(args...; kwds...))
end

function get(::Type{Array{T}}, dim::Integer,
             args...; kwds...) :: Array{T,N} where {T,N}
    get_data(Vector{T}, dim, _get1(args...; kwds...))
end

function get(::Type{Array{T}}, dims::NTuple{N,Integer},
             args...; kwds...) :: Array{T,N} where {T,N}
    get_data(Array{T,N}, dims, _get1(args...; kwds...))
end

function get(::Type{Array{T,N}}, dims::NTuple{N,Integer},
             args...; kwds...) :: Array{T,N} where {T,N}
    get_data(Array{T,N}, dims, _get1(args...; kwds...))
end

function get(::Type{String}, args...; kwds...)
    return get_data(String, _get1(args...; kwds...))
end

function _get(xpa::Client, apt::AbstractString, params::AbstractString,
              mode::AbstractString, nmax::Int, throwerrors::Bool)
    lengths = fill!(Vector{Csize_t}(undef, nmax), 0)
    buffers = fill!(Vector{Ptr{Byte}}(undef, nmax*3), Ptr{Byte}(0))
    address = pointer(buffers)
    offset = nmax*sizeof(Ptr{Byte})
    replies = ccall((:XPAGet, libxpa), Cint,
                    (Client, Cstring, Cstring, Cstring, Ptr{Ptr{Byte}},
                     Ptr{Csize_t}, Ptr{Ptr{Byte}}, Ptr{Ptr{Byte}}, Cint),
                    xpa, apt, params, mode, address, lengths,
                    address + offset, address + 2*offset, nmax)
    0 ≤ replies ≤ nmax || error("unexpected number of replies from XPAGet")
    rep = finalizer(_free, Reply(replies, lengths, buffers))
    throwerrors && verify(rep; throwerrors=true)
    return rep
end

function _free(rep::Reply)
    nmax = _nmax(rep)
    fill!(rep.lengths, 0)
    for i in 0:2,
        j in 1:length(rep)
        k = i*nmax + j
        if (ptr = rep.buffers[k]) != NULL
            rep.buffers[k] = NULL
            _free(ptr)
        end
    end
end

# i-th server name is at index i + nmax.
_get_srv(rep::Reply, i::Integer) = _get_srv(rep, Int(i))
_get_srv(rep::Reply, i::Int) :: Ptr{Byte} =
    (1 ≤ i ≤ length(rep) ? rep.buffers[i + _nmax(rep)] : NULL)

# i-th message is at index i + 2*nmax.
_get_msg(rep::Reply, i::Integer) = _get_msg(rep, Int(i))
_get_msg(rep::Reply, i::Int) :: Ptr{Byte} =
    (1 ≤ i ≤ length(rep) ? rep.buffers[i + _nmax(rep)*2] : NULL)

# i-th data buffer is at index i.
_get_buf(rep::Reply, i::Integer, preserve::Bool) =
    _get_buf(rep, Int(i), preserve)

function _get_buf(rep::Reply, i::Int, preserve::Bool) :: Tuple{Ptr{Byte},Int}
    local ptr::Ptr{Byte}, len::Int
    if 1 ≤ i ≤ length(rep)
        ptr, len = rep.buffers[i], rep.lengths[i]
        (ptr == NULL ? len == 0 : len ≥ 0) || error("invalid buffer length")
        if ! preserve && ptr != NULL
            rep.lengths[i] = 0
            rep.buffers[i] = NULL
        end
    else
        ptr, len = NULL, 0
    end
    return (ptr, len)
end

"""

Private method `_join(tup)` joins a tuple of strings or reals into a single
string.  It is implemented so as to be faster than `join(tup, " ")` when `tup`
has less than 2 arguments.  It is intended to build XPA command string from
arguments.

"""
_join(args::TupleOf{Union{AbstractString,Real}}) = join(args, " ")
_join(args::Tuple{AbstractString}) = args[1]
_join(args::Tuple{Real}) = string(args[1])
_join(::Tuple{}) = ""

"""

Private method `_nmax(n::Integer)` yields the maximum number of expected
answers to a get/set request.  The result is `n` if `n ≥ 1` or
`getconfig("XPA_MAXHOSTS")` otherwise.  The call `_nmax(rep::Reply)` yields
the maximum number of answer that can be stored in `rep`.

"""
_nmax(n::Integer) = (n == -1 ? Int(getconfig("XPA_MAXHOSTS")) : Int(n))
_nmax(rep::Reply) = length(rep.lengths)

"""
```julia
get_server(rep, i=1)
```

yields the XPA identifier of the server which sent the `i`-th reply in XPA
answer `rep`.  An empty string is returned if there is no `i`-th reply.

See also [`XPA.get`](@ref), [`XPA.get_message`](@ref).

"""
get_server(rep::Reply, i::Integer=1) = _string(_get_srv(rep, i))

"""
```julia
get_message(rep, i=1)
```

yields the message associated with the `i`-th reply in XPA answer `rep`.  An
empty string is returned if there is no `i`-th reply.

See also [`XPA.get`](@ref), [`XPA.has_message`](@ref), [`XPA.has_error`](@ref),
[`XPA.get_server`](@ref).

"""
get_message(rep::Reply, i::Integer=1) = _string(_get_msg(rep, i))

"""
```julia
XPA.has_error(rep, i=1) -> boolean
```

yields whether `i`-th XPA answer `rep` contains an error message.  The error
message can be retrieved by calling `XPA.get_message(rep, i)`.

See also [`XPA.get`](@ref), [`XPA.has_message`](@ref),
[`XPA.get_message`](@ref).

"""
has_error(rep::Reply, i::Integer=1) =
    _startswith(_get_msg(rep, i), _XPA_ERROR)

const _XPA_ERROR_PREFIX = "XPA\$ERROR "
const _XPA_ERROR = Tuple(map(Byte, collect(_XPA_ERROR_PREFIX)))

"""
```julia
XPA.has_errors(rep) -> boolean
```

yields whether answer `rep` contains any error messages.

See also [`XPA.get`](@ref), [`XPA.has_error`](@ref), [`XPA.get_message`](@ref).

"""
function has_errors(rep::Reply) :: Bool
    for i in 1:length(rep)
        if has_error(rep, i)
            return true
        end
    end
    return false
end

"""
```julia
XPA.has_message(rep, i=1) -> boolean
```

yields whether `i`-th XPA answer `rep` contains an error message.

See also [`XPA.get`](@ref), [`XPA.has_message`](@ref).

"""
has_message(rep::Reply, i::Integer=1) =
    _startswith(_get_msg(rep, i), _XPA_MESSAGE)

const _XPA_MESSAGE_PREFIX = "XPA\$MESSAGE "
const _XPA_MESSAGE = Tuple(map(Byte, collect(_XPA_MESSAGE_PREFIX)))

function _startswith(ptr::Ptr{Byte}, tup::NTuple{N,Byte}) where {N}
    if ptr == NULL
        return false
    end
    for i in 1:N
        if unsafe_load(ptr, i) != tup[i]
            return false
        end
    end
    return true
end

"""

```julia
verify(rep [, i]; throwerrors::Bool=false) -> boolean
```

verifies whether answer(s) in the result `rep` from an [`XPA.get`](@ref) or
[`XPA.set`](@ref) request has no errors.  If index `i` is specified only that
specific answer is considered; otherwise, all answers are verified.  If keyword
`throwerrors` is true, an exception is thrown for the first error found if any.

"""
function verify(rep::Reply; kwds...)
    for i in 1:length(rep)
        verify(rep, i; kwds...) || return false
    end
    return true
end

function verify(rep::Reply, i::Integer; throwerrors::Bool=false)
    if has_error(rep, i)
        if throwerrors
            # Strip error message prefix (which we know that it is present).
            msg = get_message(rep, i)
            j = length(_XPA_ERROR_PREFIX) + 1
            while j ≤ length(msg) && isspace(msg[j])
                j += 1
            end
            error(msg[j:end])
        else
            return false
        end
    end
    return true
end

"""
```julia
get_data([T, [dims,]] rep, i=1; preserve=false)
```

yields the data associated with the `i`-th reply in XPA answer `rep`.  The
returned value depends on the optional leading arguments `T` and `dims`:

* If neither `T` nor `dims` are specified, a vector of bytes (`UInt8`) is
  returned.

* If only `T` is specified, it can be `String` to return a string interpreting
  the data as ASCII characters or a type like `Vector{S}` to return the largest
  vector of elements of type `S` that can be extracted from the data.

* If both `T` and `dims` are specified, `T` can be an array type like
  `Array{S}` or `Array{S,N}` and `dims` a list of `N` dimensions to retrieve
  the data as an array of type `Array{S,N}`.

Keyword `preserve` can be used to specifiy whether or not to preserve the
internal data buffer in `rep` for another call to `XPA.get_data`.  By default,
`preserve=true` when `T = String` is specified and `preserve=false` otherwise.

In any cases, the type of the result is predictible, so there should be no type
instability issue.

See also [`XPA.get`](@ref), [`XPA.get_message`](@ref),
[`XPA.get_server`](@ref).

"""
get_data(rep::Reply, args...; kwds...) =
    get_data(Vector{Byte}, rep, args...; kwds...)

# FIXME: implement
#   get_data(Vector{String}, rep, ...) to split in words

function get_data(::Type{String}, rep::Reply, i::Integer=1;
                  preserve::Bool = true) :: String
    ptr, len = _get_buf(rep, i, preserve)
    ptr == NULL && return ""
    str = unsafe_string(ptr, len)
    preserve || _free(ptr)
    return str
end

function get_data(::Type{Vector{T}}, rep::Reply, i::Integer=1;
                  preserve::Bool = false) :: Vector{T} where {T}
    isbitstype(T) || error("invalid Array element type")
    ptr, len = _get_buf(rep, i, preserve)
    cnt = div(len, sizeof(T))
    if ptr == NULL || cnt ≤ 0
        # Empty vector.
        return T[]
    elseif preserve
        # Make a copy of the buffer in the form of a Julia vector.
        return _memcpy!(Vector{T}(undef, cnt), ptr, cnt*sizeof(T))
    else
        # Transfer ownership of buffer to Julia.
        return unsafe_wrap(Array, Ptr{T}(ptr), cnt, own=true)
    end
end

function get_data(::Type{Array{T}}, dim::Integer,
                  rep::Reply, i::Integer = 1;
                  preserve::Bool = false) :: Vector{T} where {T}
    return _get_buf(Vector{T}, _dimensions(dim), rep, i, preserve)
end

function get_data(::Type{Vector{T}}, dim::Integer,
                  rep::Reply, i::Integer = 1;
                  preserve::Bool = false) :: Vector{T} where {T}
    return _get_buf(Vector{T}, _dimensions(dim), rep, i, preserve)
end

function get_data(::Type{Array{T}}, dims::NTuple{N,Integer},
                  rep::Reply, i::Integer = 1;
                  preserve::Bool = false) :: Array{T,N} where {T,N}
    return _get_buf(Array{T,N}, _dimensions(dims), rep, i, preserve)
end

function get_data(::Type{Array{T,N}}, dims::NTuple{N,Integer},
                  rep::Reply, i::Integer = 1;
                  preserve::Bool = false) :: Array{T,N} where {T,N}
    return _get_buf(Array{T,N}, _dimensions(dims), rep, i, preserve)
end

"""

Private method `_get_buf(rep,i,preserve)` yields `(ptr,len)` the address and
length (in bytes) of internal buffer corresponding to the data for the `i`-th
reply in `rep`.  If `preserve` is false, then the internal buffer is set to
NULL and the caller is responsible to free it.  If `i` is out of range or if
there are no data associated with the `i`-th reply in `rep`, `(NULL,0)` is
returned.

The call:

```julia
_get_buf(::Type{Array{T,N}}, dims::NTuple{N,Int},
         rep::Reply, i::Integer, preserve::Bool)
```

yields the contents of the internal data buffer as a Julia array.

"""
function _get_buf(::Type{Array{T,N}}, dims::NTuple{N,Int},
                  rep::Reply, i::Int, preserve::Bool) :: Array{T,N} where {T,N}
    isbitstype(T) || error("invalid Array element type")
    minimum(dims) ≥ 0 || error("invalid Array dimensions")
    ptr, len = _get_buf(rep, i, preserve)
    cnt = prod(dims)
    cnt*sizeof(T) ≤ len || error("Array size too large for buffer")
    if ptr == NULL || cnt ≤ 0
        # Empty array.
        return Array{T,N}(undef, dims)
    elseif preserve
        # Make a copy of the buffer in the form of a Julia array.
        return _memcpy!(Array{T,N}(undef, dims), ptr, cnt*sizeof(T))
    else
        # Transfer ownership of buffer to Julia.
        return unsafe_wrap(Array, Ptr{T}(ptr), dims, own=true)
    end
end

_dimensions(dim::Integer) = (Int(dim),)
_dimensions(dim::Int) = (dim,)
_dimensions(dims::Tuple{}) = dims
_dimensions(dims::TupleOf{Integer}) = map(Int, dims)
_dimensions(dims::TupleOf{Int}) = dims

# _string(ptr) copies data at address ptr as a string, if ptr is NULL, an empty
# string is returned.
_string(ptr::Ptr{Byte}) = (ptr == NULL ? "" : unsafe_string(ptr))

"""
```julia
XPA.set([xpa,] apt, args...; data=nothing) -> rep
```

sends `data` to one or more XPA access points identified by `apt` with
arguments `args...` (automatically converted into a single string where the
arguments are separated by a single space).  The result is an instance of
[`XPA.Reply`](@ref).  Optional argument `xpa` specifies an XPA handle (created
by [`XPA.Client`](@ref)) for faster connections.

The following keywords are available:

* `data` specifies the data to send, may be `nothing`, an array or a string.
  If it is an array, it must have contiguous elements (as a for a *dense*
  array) and must implement the `pointer` method.

* `nmax` specifies the maximum number of recipients, `nmax=1` by default.
  Specify `nmax=-1` to use the maximum possible number of XPA hosts.

* `mode` specifies options in the form `"key1=value1,key2=value2"`.

* `throwerrors` specifies whether to check for errors.  If this keyword is set
  `true`, an exception is thrown for the first error message encountered in the
  list of answers.  By default, `throwerrors` is false.

See also [`XPA.Client`](@ref), [`XPA.get`](@ref) and [`XPA.verify`](@ref).

"""
function set(xpa::Client,
             apt::AbstractString,
             cmd::AbstractString;
             data = nothing,
             mode::AbstractString = "",
             nmax::Integer = 1,
             throwerrors::Bool = false)
    return _set(xpa, apt, cmd, mode, buffer(data), _nmax(nmax), throwerrors)
end

function set(xpa::Client,
             apt::AbstractString,
             args::Union{AbstractString,Real}...;
             kwds...)
    return _set(xpa, apt, _join(args); kwds...)
end

set(apt::AbstractString, args::Union{AbstractString,Real}...; kwds...) =
    set(TEMPORARY, apt, _join(args); kwds...)

function _set(xpa::Client, apt::AbstractString, params::AbstractString,
              mode::AbstractString, data::Union{NullBuffer,DenseArray},
              nmax::Int, throwerrors::Bool)

    lengths = fill!(Vector{Csize_t}(undef, nmax), 0)
    buffers = fill!(Vector{Ptr{Byte}}(undef, nmax*3), Ptr{Byte}(0))
    address = pointer(buffers)
    offset = nmax*sizeof(Ptr{Byte})
    replies = ccall((:XPASet, libxpa), Cint,
              (Client, Cstring, Cstring, Cstring, Ptr{Cvoid},
               Csize_t, Ptr{Ptr{Byte}}, Ptr{Ptr{Byte}}, Cint),
              xpa, apt, params, mode, data, sizeof(data),
              address + offset, address + 2*offset, nmax)
    0 ≤ replies ≤ nmax || error("unexpected number of replies from XPASet")
    rep = finalizer(_free, Reply(replies, lengths, buffers))
    throwerrors && verify(rep; throwerrors=true)
    return rep
end

"""
```julia
buf = buffer(data)
```

yields an object `buf` representing the contents of `data` and which can be
used as an argument to `ccall` without the risk of having the data garbage
collected.  Argument `data` can be `nothing`, a dense array or a string.  If
`data` is an array `buf` is just an alias for `data`.  If `data` is a string,
`buf` is a temporary byte buffer where the string has been copied.

Standard methods `pointer` and `sizeof` can be applied to `buf` to retieve the
address and the size (in bytes) of the data and `convert(Ptr{Cvoid},buf)` can
also be used.

See also [`XPA.set`](@ref).

"""
function buffer(arr::A) :: A where {T,N,A<:DenseArray{T,N}}
    @assert isbitstype(T)
    return arr
end

function buffer(str::AbstractString)
    @assert isascii(str)
    len = length(str)
    buf = Vector{Cchar}(undef, len)
    @inbounds for i in 1:len
        buf[i] = str[i]
    end
    return buf
end

buffer(::Nothing) = NullBuffer()

Base.unsafe_convert(::Type{Ptr{T}}, ::NullBuffer) where {T} = Ptr{T}(0)
Base.pointer(::NullBuffer) = C_NULL
Base.sizeof(::NullBuffer) = 0


text:
# Naive implementations of multiply_dimensionwise used to demonstrate the functionality
# without performance optimizations and for testing correctness of the optimized versions
# implemented below.
function multiply_dimensionwise_naive(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, 2})
  size_out = size(matrix, 1)
  size_in  = size(matrix, 2)
  n_vars   = size(data_in, 1)
  data_out = zeros(promote_type(eltype(data_in), eltype(matrix)), n_vars, size_out)

  for i in 1:size_out
    for ii in 1:size_in
      for v in 1:n_vars
        data_out[v, i] += matrix[i, ii] * data_in[v, ii]
      end
    end
  end

  return data_out
end

function multiply_dimensionwise_naive(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, 3})
  size_out = size(matrix, 1)
  size_in  = size(matrix, 2)
  n_vars   = size(data_in, 1)
  data_out = zeros(promote_type(eltype(data_in), eltype(matrix)), n_vars, size_out, size_out)

  for j in 1:size_out, i in 1:size_out
    for jj in 1:size_in, ii in 1:size_in
      for v in 1:n_vars
        data_out[v, i, j] += matrix[i, ii] * matrix[j, jj] * data_in[v, ii, jj]
      end
    end
  end

  return data_out
end

function multiply_dimensionwise_naive(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, 4})
  size_out = size(matrix, 1)
  size_in  = size(matrix, 2)
  n_vars   = size(data_in, 1)
  data_out = zeros(promote_type(eltype(data_in), eltype(matrix)), n_vars, size_out, size_out, size_out)

  for k in 1:size_out, j in 1:size_out, i in 1:size_out
    for kk in 1:size_in, jj in 1:size_in, ii in 1:size_in
      for v in 1:n_vars
        data_out[v, i, j, k] += matrix[i, ii] * matrix[j, jj] * matrix[k, kk] * data_in[v, ii, jj, kk]
      end
    end
  end

  return data_out
end

"""
    multiply_dimensionwise(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, NDIMS+1})

Multiply the array `data_in` by `matrix` in each coordinate direction, where `data_in`
is assumed to have the first coordinate for the number of variables and the remaining coordinates
are multiplied by `matrix`.
"""
function multiply_dimensionwise(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, 2})
  # 1D
  # optimized version of multiply_dimensionwise_naive
  size_out = size(matrix, 1)
  n_vars   = size(data_in, 1)
  data_out = zeros(promote_type(eltype(data_in), eltype(matrix)), n_vars, size_out)

  multiply_dimensionwise!(data_out, matrix, data_in)

  return data_out
end

function multiply_dimensionwise(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, 3})
  # 2D
  # optimized version of multiply_dimensionwise_naive
  size_out = size(matrix, 1)
  n_vars   = size(data_in, 1)
  data_out = zeros(promote_type(eltype(data_in), eltype(matrix)), n_vars, size_out, size_out)

  multiply_dimensionwise!(data_out, matrix, data_in)

  return data_out
end

function multiply_dimensionwise(matrix::AbstractMatrix, data_in::AbstractArray{<:Any, 4})
  # 3D
  # optimized version of multiply_dimensionwise_naive
  size_out = size(matrix, 1)
  n_vars   = size(data_in, 1)
  data_out = zeros(promote_type(eltype(data_in), eltype(matrix)), n_vars, size_out, size_out, size_out)

  multiply_dimensionwise!(data_out, matrix, data_in)

  return data_out
end

# In the following, there are several optimized in-place versions of multiply_dimensionwise.
# These make use of the macro `@tullio` from Tullio.jl, which basically uses an Einstein
# summation convention syntax.

# 1D version
function multiply_dimensionwise!(data_out::AbstractArray{<:Any, 2}, matrix::AbstractMatrix,
                                 data_in ::AbstractArray{<:Any, 2})
  @tullio threads=false data_out[v, i] = matrix[i, ii] * data_in[v, ii]

  return nothing
end

# 1D version for scalars
# Instead of having a leading dimension of size 1 in `data_out, data_in`, this leading dimension
# of size unity is dropped, resulting in one dimension less than in `multiply_dimensionwise!`.
function multiply_scalar_dimensionwise!(data_out::AbstractArray{<:Any, 1},
                                        matrix::AbstractMatrix,
                                        data_in ::AbstractArray{<:Any, 1})
  @tullio threads=false data_out[i] = matrix[i, ii] * data_in[ii]

  return nothing
end

# 1D version, apply matrixJ to data_inJ
function multiply_dimensionwise!(data_out::AbstractArray{<:Any, 2}, matrix1::AbstractMatrix,
                                 data_in1::AbstractArray{<:Any, 2}, matrix2::AbstractMatrix,
                                 data_in2::AbstractArray{<:Any, 2})
  @tullio threads=false data_out[v, i] = matrix1[i, ii] * data_in1[v, ii] + matrix2[i, ii] * data_in2[v, ii]

  return nothing
end

# 2D version
function multiply_dimensionwise!(data_out::AbstractArray{<:Any, 3}, matrix::AbstractMatrix,
                                 data_in:: AbstractArray{<:Any, 3},
                                 tmp1=zeros(eltype(data_out), size(data_out, 1), size(matrix, 1), size(matrix, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[v, i, j]     = matrix[i, ii] * data_in[v, ii, j]

  # Interpolate in y-direction
  @tullio threads=false data_out[v, i, j] = matrix[j, jj] * tmp1[v, i, jj]

  return nothing
end

# 2D version for scalars
# Instead of having a leading dimension of size 1 in `data_out, data_in`, this leading dimension
# of size unity is dropped, resulting in one dimension less than in `multiply_dimensionwise!`.
function multiply_scalar_dimensionwise!(data_out::AbstractArray{<:Any, 2},
                                        matrix::AbstractMatrix,
                                        data_in:: AbstractArray{<:Any, 2},
                                        tmp1=zeros(eltype(data_out), size(matrix, 1), size(matrix, 2)))

  # Interpolate in x-direction
  @tullio threads=false     tmp1[i, j] = matrix[i, ii] * data_in[ii, j]

  # Interpolate in y-direction
  @tullio threads=false data_out[i, j] = matrix[j, jj] * tmp1[i, jj]

  return nothing
end

# 2D version, apply matrixJ to dimension J of data_in
function multiply_dimensionwise!(data_out::AbstractArray{<:Any, 3},
                                 matrix1::AbstractMatrix, matrix2::AbstractMatrix,
                                 data_in:: AbstractArray{<:Any, 3},
                                 tmp1=zeros(eltype(data_out), size(data_out, 1), size(matrix1, 1), size(matrix1, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[v, i, j]     = matrix1[i, ii] * data_in[v, ii, j]

  # Interpolate in y-direction
  @tullio threads=false data_out[v, i, j] = matrix2[j, jj] * tmp1[v, i, jj]

  return nothing
end

# 2D version, apply matrixJ to dimension J of data_in and add the result to data_out
function add_multiply_dimensionwise!(data_out::AbstractArray{<:Any, 3},
                                     matrix1::AbstractMatrix, matrix2::AbstractMatrix,
                                     data_in:: AbstractArray{<:Any, 3},
                                     tmp1=zeros(eltype(data_out), size(data_out, 1), size(matrix1, 1), size(matrix1, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[v, i, j]     = matrix1[i, ii] * data_in[v, ii, j]

  # Interpolate in y-direction
  @tullio threads=false data_out[v, i, j] += matrix2[j, jj] * tmp1[v, i, jj]

  return nothing
end

# 3D version
function multiply_dimensionwise!(data_out::AbstractArray{<:Any, 4}, matrix::AbstractMatrix,
                                 data_in:: AbstractArray{<:Any, 4},
                                 tmp1=zeros(eltype(data_out), size(data_out, 1), size(matrix, 1), size(matrix, 2), size(matrix, 2)),
                                 tmp2=zeros(eltype(data_out), size(data_out, 1), size(matrix, 1), size(matrix, 1), size(matrix, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[v, i, j, k]     = matrix[i, ii] * data_in[v, ii, j, k]

  # Interpolate in y-direction
  @tullio threads=false tmp2[v, i, j, k]     = matrix[j, jj] * tmp1[v, i, jj, k]

  # Interpolate in z-direction
  @tullio threads=false data_out[v, i, j, k] = matrix[k, kk] * tmp2[v, i, j, kk]

  return nothing
end

# 3D version for scalars
# Instead of having a leading dimension of size 1 in `data_out, data_in`, this leading dimension
# of size unity is dropped, resulting in one dimension less than in `multiply_dimensionwise!`.
function multiply_scalar_dimensionwise!(data_out::AbstractArray{<:Any, 3}, matrix::AbstractMatrix,
                                        data_in:: AbstractArray{<:Any, 3},
                                        tmp1=zeros(eltype(data_out), size(matrix, 1), size(matrix, 2), size(matrix, 2)),
                                        tmp2=zeros(eltype(data_out), size(matrix, 1), size(matrix, 1), size(matrix, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[i, j, k]     = matrix[i, ii] * data_in[ii, j, k]

  # Interpolate in y-direction
  @tullio threads=false tmp2[i, j, k]     = matrix[j, jj] * tmp1[i, jj, k]

  # Interpolate in z-direction
  @tullio threads=false data_out[i, j, k] = matrix[k, kk] * tmp2[i, j, kk]

  return nothing
end

# 3D version, apply matrixJ to dimension J of data_in
function multiply_dimensionwise!(data_out::AbstractArray{<:Any, 4},
                                 matrix1::AbstractMatrix, matrix2::AbstractMatrix, matrix3::AbstractMatrix,
                                 data_in:: AbstractArray{<:Any, 4},
                                 tmp1=zeros(eltype(data_out), size(data_out, 1), size(matrix1, 1), size(matrix1, 2), size(matrix1, 2)),
                                 tmp2=zeros(eltype(data_out), size(data_out, 1), size(matrix1, 1), size(matrix1, 1), size(matrix1, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[v, i, j, k]     = matrix1[i, ii] * data_in[v, ii, j, k]

  # Interpolate in y-direction
  @tullio threads=false tmp2[v, i, j, k]     = matrix2[j, jj] * tmp1[v, i, jj, k]

  # Interpolate in z-direction
  @tullio threads=false data_out[v, i, j, k] = matrix3[k, kk] * tmp2[v, i, j, kk]

  return nothing
end

# 3D version, apply matrixJ to dimension J of data_in and add the result to data_out
function add_multiply_dimensionwise!(data_out::AbstractArray{<:Any, 4},
                                     matrix1::AbstractMatrix, matrix2::AbstractMatrix, matrix3::AbstractMatrix,
                                     data_in:: AbstractArray{<:Any, 4},
                                     tmp1=zeros(eltype(data_out), size(data_out, 1), size(matrix1, 1), size(matrix1, 2), size(matrix1, 2)),
                                     tmp2=zeros(eltype(data_out), size(data_out, 1), size(matrix1, 1), size(matrix1, 1), size(matrix1, 2)))

  # Interpolate in x-direction
  @tullio threads=false tmp1[v, i, j, k]     = matrix1[i, ii] * data_in[v, ii, j, k]

  # Interpolate in y-direction
  @tullio threads=false tmp2[v, i, j, k]     = matrix2[j, jj] * tmp1[v, i, jj, k]

  # Interpolate in z-direction
  @tullio threads=false data_out[v, i, j, k] += matrix3[k, kk] * tmp2[v, i, j, kk]

  return nothing
end


text:
using Statistics
using Test

using Dates, DataFrames, CSV
using Streamfall

here = @__DIR__
climate_data = joinpath(here, "data/campaspe/climate/climate_historic.csv")
dam_data_loc = joinpath(here, "data/campaspe/dam")

# Read in test data
climate_data = CSV.File(climate_data, comment="#", dateformat="YYYY-mm-dd") |> DataFrame
hist_dam_levels = CSV.File(joinpath(dam_data_loc, "historic_levels_for_fit.csv"), dateformat="YYYY-mm-dd") |> DataFrame
hist_dam_releases = CSV.File(joinpath(dam_data_loc, "historic_releases.csv"), dateformat="YYYY-mm-dd") |> DataFrame

@testset "Min/max date finding" begin
    min_date, max_date = Streamfall.find_common_timeframe(climate_data, hist_dam_levels, hist_dam_releases)

    @test min_date == Date("1981-01-01")
    @test max_date == Date("2015-07-30")
end

@testset "Data alignment" begin

    # Mismatched start/end dates
    dam_levels = hist_dam_levels[10:end, :]
    dam_releases = hist_dam_releases[20:10000, :]

    # Subset to common time period
    c, dl, dr = Streamfall.align_time_frame(climate_data, dam_levels, dam_releases)

    # Ensure time periods are the same
    @test c.Date[1] == dl.Date[1] == dr.Date[1]
    @test c.Date[end] == dl.Date[end] == dr.Date[end]

    # Ensure originals have not been modified
    @test climate_data.Date[1] != c.Date[1]
    @test hist_dam_levels.Date[1] != dl.Date[1]
    @test hist_dam_releases.Date[1] != dr.Date[1]

    @test hist_dam_levels.Date[end] != dl.Date[end]
    @test hist_dam_releases.Date[end] != dr.Date[end]
end


text:
Base.@kwdef struct LearningTaskRegistry
    tasks::Dict{Any,Any} = Dict()
end

"""
    registerlearningtask!(registry, taskfn, blocktypes)

Register a learning task constructor `taskfn` as compatible with
`blocktypes` in `registry::LearningTaskRegistry`.

`blocks` should be the least specific set of types that a `taskfn`
can handle. `taskfn` needs to have a task that takes concrete block
instances as the only non-keyword argument, i.e. `taskfn(blocks; kwargs...)`.
"""
function registerlearningtask!(reg::LearningTaskRegistry, f, blocktypes)
	reg.tasks[f] = typify(blocktypes)
	return reg
end

"""
	findlearningtasks(blocktypes)
	findlearningtasks(registry, blocktypes)

Find learning tasks compatible with block types `TBlocks` in
`registry::LearningTaskRegistry`.

#### Examples

```julia-repl
julia> findlearningtasks((Image, Label))
[ImageClassificationSingle,]

julia> findlearningtasks((Image, Any))
[ImageClassificationSingle, ImageClassificationMulti, ImageSegmentation, ImageKeypointRegression, ...]
```
"""
function findlearningtasks(reg::LearningTaskRegistry, blocktypes=Any)
	return [taskfn for (taskfn, taskblocks) in reg.tasks if blocktypesmatch(taskblocks, blocktypes)]
end

function blocktypesmatch(
        BSupported::Type{<:AbstractBlock},
        BWanted::Type{<:AbstractBlock})
    BWanted <: BSupported
end
function blocktypesmatch(B1::Type{<:Tuple}, B2::Type{<:Tuple})
    all(blocktypesmatch(b1, b2) for (b1, b2) in zip(B1.types, B2.types))
end

blocktypesmatch(BSupported::Type, _::Type{Any}) = true
blocktypesmatch(bsupported, bwanted) = blocktypesmatch(typify(bsupported), bwanted)
blocktypesmatch(BSupported::Type, bwanted) = blocktypesmatch(BSupported, typify(bwanted))

@testset "`blocktypesmatch`" begin
    @test blocktypesmatch(FastAI.Image, FastAI.Image{2})
    @test blocktypesmatch(FastAI.Image, FastAI.Image)
    @test !blocktypesmatch(FastAI.Image{2}, FastAI.Image)
    @test blocktypesmatch(Tuple{FastAI.Image}, Tuple{FastAI.Image})
    @test !blocktypesmatch(Tuple{FastAI.Image{2}}, Tuple{FastAI.Image})
    @test blocktypesmatch(Tuple{FastAI.Image{2}}, Any)
    @test blocktypesmatch(FastAI.Image{2}(), FastAI.Image{2})
    @test blocktypesmatch(FastAI.Image, FastAI.Image{2}())
    @test blocktypesmatch((FastAI.Image, FastAI.Label), (FastAI.Image{2}(), FastAI.Label(1:10)))
    @test blocktypesmatch((FastAI.Image, FastAI.Label), (FastAI.Image{2}(), Any))
end

const FASTAI_METHOD_REGISTRY = LearningTaskRegistry()


text:
"""
ERGM for directed networks with totan links and tortal reciprocated links as statistics
"""
Base.@kwdef struct  ErgmDirBin0Rec0 <: Ergm
    ergmTermsString::String="edges + mutual" # needs to be compatible with R ergm package
    nErgmPar = 2
end

ergm_term_string(x::ErgmDirBin0Rec0)  = "edges + mutual"

type_of_obs(model::ErgmDirBin0Rec0) =  Array{Float64, 1}

"""
given the 2 model's parameters compute the probabilities for each state of the diad
"""
function diadProbFromPars(Model::ErgmDirBin0Rec0, par)
    θ = par[1]
    η = par[2]
    eθ = exp(θ)
    e2θplusη = exp(2*θ + η)
    den = 1 + 2*eθ + e2θplusη
    p00 = 1/den
    p01 = eθ /den
    p10 = p01
    p11 = e2θplusη/den
    diadProbsVec = [p00, p01, p10, p11]
    return diadProbsVec
end

"""
given the vector of diad states probabilities sample one random matrix from the corresponding pdf
"""
function samplSingMatCan(Model::ErgmDirBin0Rec0, diadProbsVec::Array{<:Real,1}, N)
    out = zeros(Int8,N,N)
    #display((maximum(expMat),minimum(expMat)))
    for c = 1:N
        for r=1:c-1
            diadState = rand(Categorical(diadProbsVec))
            if diadState==2
                out[r,c] = 1
            elseif diadState==3
                out[c,r] = 1
            elseif diadState==4
                out[r,c] = 1
                out[c,r] = 1
            end
        end
    end
    return out
end

function sample_ergm(model::ErgmDirBin0Rec0, N, par_vec, nSample)

    θ_0, η_0 = par_vec

    diadProb = diadProbFromPars(model, [θ_0, η_0])

    A_vec =  reduce((x,y)->cat(x,y,dims=3), [samplSingMatCan(model, diadProb, N) for i=1:nSample])

    return A_vec
end

function stats_from_mat(Model::ErgmDirBin0Rec0, A ::Matrix{<:Real})
    L = sum(A)
    R = sum(A'.*A)/2
    N=size(A)[1]

    return [L, R, N]
end

function exp_val_stats(model::ErgmDirBin0Rec0, θ, η, N)
    x = exp(θ)
    x2y = exp(2*θ + η)
    z = 1 + 2*x + x2y
    Nlinks = n_pox_dir_links(N)
    α = (x + x2y)/z
    β = 0.5 * (x2y)/z
    L = Nlinks * α
    R = Nlinks * β
    return L, R
end

"""
return the ergm parameters that fix on average the input values for the ergm statistics:
    - L is  ∑_i>j A_ij + A_ji
    - R is  of ∑_i>j A_ij * A_ji /2
"""
function ergm_par_from_mean_vals(model::ErgmDirBin0Rec0, L, R, N)
    Nlinks = n_pox_dir_links(N)
    # average values per pair
    α = L/Nlinks
    β = R/Nlinks

    x = (2β-α)/(2α-2β-1)
    θ = log(x)

    y = 2(2*β^2  + β- 2*α*β )/((2β-α)^2)
    η =  log(y)

    return θ, η
end

function estimate(model::ErgmDirBin0Rec0, L, R, N)
    θ_est, η_est = ergm_par_from_mean_vals(model, L, R, N)
    vPar = [θ_est, η_est]
    return vPar
end

function estimate(Model::ErgmDirBin0Rec0, A::T where T <: Matrix)
    L, R, N = stats_from_mat(Model, A)
    return estimate(Model, L, R, N)
end

function logLikelihood(Model::ErgmDirBin0Rec0, L, R, N, par)
    θ, η = par
    z = 1 + 2*exp(θ) + exp(2*θ+η)
    return L * θ + R*η - (N*(N-1)/2)*log(z)
end

function logLikelihood(Model::ErgmDirBin0Rec0, A::Matrix, par)
    L, R, N = stats_from_mat(Model, A)
    return logLikelihood(Model, L, R, N, par)
end

estimate_sequence(model::ErgmDirBin0Rec0, obsT::Array{Array{T, 1},1} where T <: Real) = reduce(hcat, [estimate(model, obsT[t]... ) for t in 1:length(obsT)])

estimate_sequence(model::ErgmDirBin0Rec0, AT::Array{T, 3} where T <: Real) = reduce(hcat, [estimate(model, AT[:,:,t]) for t in 1:size(AT)[3]])


text:
#solve eq x=cos(x)
function main()
    x=0.0
    xnext = cos(x)
    iter=0
    eps=1e-7
    iteration=1000000
    is_convergence=false

    for n in 1:iteration
        xnext = cos(x)
        if(abs(xnext-x)<eps)
            is_convergence=true
            iter=n
            break
        end
        x=xnext
    end

    if is_convergence
        println("X=",xnext,", iter=",iter)
        println("X-cos(X)=",xnext-cos(xnext))
    else
        print("No Convergence")
    end
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer